[Insight-developers] fast image pixel access method

[Gao, Yi]

Hi Matt, Hi Kent,

Thank you for the suggestions!

I have tried the NeighborhoodIterator and then I did a little test
comparing that with directly calling GetPixel() for each of the 27
neighbors.

Basically, I randomly pick 10,000,000 points inside the "inner-region" of
the an 3D volume. For each point, take average of all its 27 neighbors, and
set its value to the average value.

Since the random points are all in [1 ,size[i]-2], range, accessing their
+/-1 neighbors are safe. So I did not check boundary when accessing.

I just used linux command "time" to time the code, as a rough test.

Both takes about 5.5seconds (Release build, again Release version of ITK
4.3.0). There is not dramatic difference...

When using NeighborhoodIterator, i'm not sure if the boundary checking is
automatically turned on. I guess if that is by default on, then there is
not much speed gain with it.

I paste the code below. Could I know if there's something I didn't do
right? Thank you!


The code are:

1. using NeighborhoodIterator

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkNeighborhoodIterator.h"

#include <vnl/vnl_random.h>

int main( int argc, char ** argv )
{
  if ( argc < 3 )
    {
      std::cerr << "Usage: " << std::endl;
      std::cerr << argv[0]
                << " inputImage numberOfPoints\n";
      return -1;
    }

  typedef float                                  PixelType;
  typedef itk::Image< PixelType, 3 >             ImageType;

  typedef itk::ImageFileReader< ImageType > ImageReaderType;
  typename ImageReaderType::Pointer reader = ImageReaderType::New();
  reader->SetFileName(argv[1]);
  reader->Update();
  typename ImageType::Pointer image = reader->GetOutput();

  ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();

  typedef itk::NeighborhoodIterator< ImageType > NeighborhoodIteratorType;

  NeighborhoodIteratorType::RadiusType radius;
  radius.Fill(1);
  NeighborhoodIteratorType it(radius, image,
image->GetLargestPossibleRegion());

  ImageType::IndexType index;
  long n = ::atol(argv[2]);
  std::cout<<n<<std::endl;

  vnl_random rg;

  PixelType s = static_cast<PixelType>(it.Size());
  std::cout<<s<<std::endl;

  for (long ip = 0; ip < n; ++ip)
    {
      index[0] = rg.lrand32(1, size[0] - 1);
      index[1] = rg.lrand32(1, size[1] - 1);
      index[2] = rg.lrand32(1, size[2] - 1);

      it.SetLocation(index);
      //std::cout<<index<<std::endl;

      PixelType sum = it.GetCenterPixel();
      for (unsigned i = 0; i < it.Size(); i++)
        {
          sum += it.GetPixel(i);
        }

      it.SetCenterPixel( sum/s );

      //std::cout<<sum/s<<std::endl;
    }

  typedef itk::ImageFileWriter< ImageType > WriterType;

  typename WriterType::Pointer writer = WriterType::New();
  writer->SetFileName( "o.nrrd" );
  writer->SetInput(image);
  writer->UseCompressionOff();
  writer->Update();

  return 0;
}




2. using direct access
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include <vnl/vnl_random.h>

int main( int argc, char ** argv )
{
  if ( argc < 3 )
    {
      std::cerr << "Usage: " << std::endl;
      std::cerr << argv[0]
                << " inputImage numberOfPoints\n";
      return -1;
    }

  typedef float                                  PixelType;
  typedef itk::Image< PixelType, 3 >             ImageType;

  typedef itk::ImageFileReader< ImageType > ImageReaderType;
  typename ImageReaderType::Pointer reader = ImageReaderType::New();
  reader->SetFileName(argv[1]);
  reader->Update();
  typename ImageType::Pointer image = reader->GetOutput();

  ImageType::SizeType size = image->GetLargestPossibleRegion().GetSize();

  ImageType::IndexType index;
  ImageType::IndexType nbhdindex;

  long n = ::atol(argv[2]);
  std::cout<<n<<std::endl;

  vnl_random rg;

  PixelType s = 27.0;

  for (long ip = 0; ip < n; ++ip)
    {
      PixelType sum = 0.0;

      int cx = rg.lrand32(1, size[0] - 2);
      int cy = rg.lrand32(1, size[1] - 2);
      int cz = rg.lrand32(1, size[2] - 2);

      index[0] = cx;
      index[1] = cy;
      index[2] = cz;

      for (int iz = cz - 1; iz <= cz + 1; ++iz)
        {
          nbhdindex[2] = iz;
          for (int iy = cy - 1; iy <= cy + 1; ++iy)
            {
              nbhdindex[1] = iy;
              for (int ix = cx - 1; ix <= cx + 1; ++ix)
                {
                  nbhdindex[0] = ix;

                  sum += image->GetPixel(nbhdindex);
                }
            }
        }

      image->SetPixel(index, sum/s);
    }

  typedef itk::ImageFileWriter< ImageType > WriterType;

  typename WriterType::Pointer writer = WriterType::New();
  writer->SetFileName( "p.nrrd" );
  writer->SetInput(image);
  writer->UseCompressionOff();
  writer->Update();

  return 0;
}











On Mon, May 13, 2013 at 10:02 AM, Williams, Norman K <
norman-k-williams at uiowa.edu> wrote:

> I would think an itk::NeighborhoodIterator with a radius of [1,1,1] would
> do.  Are there any requirements other than direct adjacency to the
> neighborhood you want to traverse?
>
> I haven't looked closely at the code but I think it's about as well
> optimized as it could be.
> --
> Kent Williams norman-k-williams at uiowa.edu
>
>
>
>
>
>
> On 5/12/13 8:43 PM, "Gao, Yi" <gaoyi.cn at gmail.com> wrote:
>
> >Dear all,
> >
> >
> >Could I have some suggestions on ways to fast access the nerighbors of a
> >set of random voxels? The scenario is as follows:
> >
> >
> >I have to access some voxels in the volume, these voxels are quite
> >randomly distributed. For each of them I know their (x, y, z) indexes.
> >
> >At each of these voxel, I need to access all its 26 neighbors. I want
> >this step to be faster.
> >
> >
> >So i'm hoping to have something like:
> >
> >
> >1. for each (x, y, z), I get the pointer to that voxel: ptr.
> >
> >2. For each of its neighbors, I use things like ptr[1], ptr[nx], ptr[1+
> >nx*ny] to access (x+1, y, z), (x, y+1, z) and (x+1, y, z+1)....
> >
> >
> >Could I know if there is some good way to do that?
> >
> >Thank you in advance for any hint!
> >
> >Best,
> >yi
> >
>
>
>
> ________________________________
> Notice: This UI Health Care e-mail (including attachments) is covered by
> the Electronic Communications Privacy Act, 18 U.S.C. 2510-2521, is
> confidential and may be legally privileged.  If you are not the intended
> recipient, you are hereby notified that any retention, dissemination,
> distribution, or copying of this communication is strictly prohibited.
>  Please reply to the sender that you have received the message in error,
> then delete it.  Thank you.
> ________________________________
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-developers/attachments/20130513/d46dcd5e/attachment.htm>