[Insight-users] Image index limits

[Oleksandr Dzyubak]

Hi Bill,

I am still working on a large image limitation problem.
This time let me get rid of all nuances
that keep drawing us away from the point.

First, dear ITK users, I am very sorry but to still stay on this forum list,
I have to cut attached answers. Otherwise the nabble system
will bounce me off.

Now let the party begin.

1) At this point I only use "itkImage.h" and "itkImageFileWriter.h".
I wrote a code (attached below) which allocates a 1340x1300x720 buffer
to hold a "float" pixel type, fill in some voxel values,
and write the final result down.

    a) This code does not use any external image files.
    b) This code first fills up the entire volume with Intensity I=100;
    c) Since I already know that the "problematic boundary" is
        slice# 104, I create two 100x100x100 cubes with Intensity=1000:
        first cube is before the critical boundary and the second one is 
after.

    d) Write the result down using the "itkImageFileWriter.h"

What is the result?
Both cubes are where they supposed to be with correct intensities.
So no more discussions about Debian Linux, gcc 4.3.2 limitations,
C++ container and memory limits.

All the above supports the volumes I need
and the ITK library can perfectly handle this.

Now lets try to find out where the real problem resides.

2) At this point I am going to use "itkImage.h",  "itkImageFileWriter.h".
and "itkImageFileReader.h".
I wrote a code (attached below) to make plain reading the image from a file
and writing it down to another file. Nothing fancy.

A test image is the image produced by the program above:
two cubes at the beginning and the end of the image.

What is the result?
Surprise! Second cube disappeared and all slices after
the critical boundary are empty.

Thus I conclude that these are not Gaussians and Gaussian derivatives
which generate a problem.

There is something wrong with the "itkImageFileReader.h" class.
Well, it is getting warmer.
I would predict that we can solve the problem very soon.

Bill, what would you say?

Regards,

Alex


****** Begin ImageFill_cube.cxx ***

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

int main(int argc, char * argv[])
{

  if( argc < 2 )
    {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0] << "  outputImageFile" << std::endl;
    return 1;
    }

  const unsigned int      Dimension = 3;

//  typedef unsigned char   PixelType;
//  typedef unsigned int    PixelType;
//  typedef short          PixelType;
  typedef float          PixelType;

  typedef itk::Image< PixelType, Dimension > ImageType;

  ImageType::Pointer image = ImageType::New();

  // The image region should be initialized
  ImageType::IndexType start;
  ImageType::SizeType  size;

  size[0]  = 1340;  // size along X
  size[1]  = 1300;  // size along Y
  size[2]  =  720;  // size along Z

  start[0] =   0;  // first index on X
  start[1] =   0;  // first index on Y
  start[2] =   0;  // first index on Z

  ImageType::RegionType region;
  region.SetSize( size );
  region.SetIndex( start );
 
  // Pixel data is allocated
  image->SetRegions( region );
  image->Allocate();

  // The image buffer is initialized to a particular value
  ImageType::PixelType  initialValue = 100;
  image->FillBuffer( initialValue );


  ImageType::IndexType pixelIndex;
  ImageType::IndexType initialPixelIndex;
  ImageType::PixelType   pixelValue = 1000;
 
  initialPixelIndex[0] = 670;   // x position
  initialPixelIndex[1] = 650;   // y position
  initialPixelIndex[2] =   0;   // z position

  std::cout << std::endl;
  std::cout << "pixelIndex = " << pixelIndex << std::endl;
  std::cout << std::endl;


// Lets draw two cubes 100x100x100 with intensity I=1000;
// Before the critical boundary pixelIndex[2]=104 and after.

  for (int i=0; i<100; i++)
    {
      pixelIndex[0] = initialPixelIndex[0] + i;
    for(int j=0; j<100; j++)
        {
          pixelIndex[1] = initialPixelIndex[1] + j;
        for(int k=0; k<100; k++)
            {
              pixelIndex[2] = initialPixelIndex[2] + k;
              image->SetPixel( pixelIndex, pixelValue );
//              std::cout << "pixelIndex = " << pixelIndex << std::endl;
//              std::cout << "pixelValue = " << image->GetPixel( 
pixelIndex ) << std::endl;
            }
        }
    }


  initialPixelIndex[0] = 670;   // x position
  initialPixelIndex[1] = 650;   // y position
  initialPixelIndex[2] = 360;   // z position

  for (int i=0; i<100; i++)
    {
      pixelIndex[0] = initialPixelIndex[0] + i;
    for(int j=0; j<100; j++)
        {
          pixelIndex[1] = initialPixelIndex[1] + j;
        for(int k=0; k<100; k++)
            {
              pixelIndex[2] = initialPixelIndex[2] + k;
              image->SetPixel( pixelIndex, pixelValue );
//              std::cout << "pixelIndex = " << pixelIndex << std::endl;
//              std::cout << "pixelValue = " << image->GetPixel( 
pixelIndex ) << std::endl;
            }
        }
    }


  typedef itk::ImageFileWriter< ImageType > WriterType;
  WriterType::Pointer writer = WriterType::New();

  writer->SetFileName( argv[1] );
  writer->SetInput( image );

  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & exp )
    {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << exp << std::endl;
    }

return EXIT_SUCCESS;
}
****** End ImageFill_cube.cxx *****


****** Begin JustReadWrite.cxx ***

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

int main( int argc, char * argv[] )
{

  if( argc < 3 )
    {
    std::cerr << "Usage: " << argv[0]  << std::endl;
    std::cerr << " inputImageFile  outputImageFile " << std::endl; 
    return EXIT_FAILURE;
    }

  const unsigned int      Dimension = 3;

//  typedef unsigned char   PixelType;
//  typedef unsigned int    PixelType;
//  typedef short          PixelType;
  typedef float          PixelType;


  typedef itk::Image< PixelType, Dimension > ImageType;
  typedef itk::ImageFileReader< ImageType >  ReaderType;
  typedef itk::ImageFileWriter< ImageType >  WriterType;

  ReaderType::Pointer reader = ReaderType::New();
  WriterType::Pointer writer = WriterType::New();


  reader->SetFileName( argv[1] );
  writer->SetFileName( argv[2] );

  writer->SetInput( reader->GetOutput() );
  writer->Update();

  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & exp )
    {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << exp << std::endl;
    }

  return EXIT_SUCCESS;
}

****** End JustReadWrite.cxx ***

Bill Lorensen wrote:
> This code is incorrect:
>  for (int i=0; i<10; i++)
>   {
>     pixelIndex[0] +=   pixelIndex[0];
>     ImageType::PixelType   pixelValue = image->GetPixel( pixelIndex );
>     std::cout << "pixelValue = " << pixelValue << std::endl;
>   }
>
> pixelIndex[0] rapidly gets out of range.
>
>