[Insight-users] Image index limits

[Oleksandr Dzyubak]

Hi Luis,

The problem of empty slices when processing large images
still bothers me and this time I decided to make direct tests.
I wrote a tiny program which allocates and fills a buffer 
(1340x1300x720) with
constant values I=100. A code snippet is attached below.

1) This program does not read an external file
so we do not use my "corrupted" images anymore plus
we exclude influence of the ImageRead class.

2) Next, I directly read intensities from a buffer before and after
a "magic boundary" (slice #104) and print the values to a console thus
we can see whether or not this failure is caused by the allocation 
procedure.

3) In addition to 2), I read the buffer content and write it to a file 
and measure
the intensities again to see whether the ImageWrite class causes a problem.

4) Finally, I did massive tests for all possible pixel types.

Believe it or not, but it works perfectly well and produces
correct intensities over the entire volume 1340x1300x720.

So I conclude that at least this part works correctly.

Next step. Now I use the example from the ITK library

InsightToolkit-3.16.0/Examples/IO/ImageReadRegionOfInterestWrite.cxx

and I try to extract couple slices before and after the "magic boundary".

Wow! After the "magic boundary", I get all the slices with zeros.

What bothers me the most is that I actually need the classes which use
the Gaussians and Gaussian derivatives. Surprisingly those algorithms
have the same problem: after the "magic boundary", all the slices
they produce are with zeros.

As to me, it looks like they all use the same technique to fly over indices
which results in a sort of "partial volume effect"
discarding intensities with large index values.

Unless I solve this problem, I cannot go any further with
using Gaussians and Gaussian derivatives for processing large images.

Luis, I would appreciate any suggestion or any crazy hint.

Regards,

Alex


******** BeginCodeSnippet *********

#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;
    }
 
//  typedef unsigned char   PixelType;
//  typedef unsigned int    PixelType;
//  typedef short          PixelType;
  typedef float          PixelType;

  const unsigned int      Dimension = 3;
  typedef itk::Image< PixelType, Dimension > ImageType;

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

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

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

  size[0]  = 1340;  // size along X
  size[1]  = 1300;  // size along Y
  size[2]  =  720;  // size along 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;
 
  pixelIndex[0] = 670;   // x position
  pixelIndex[1] = 650;   // y position
  pixelIndex[2] =  10;   // z position

  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] = 670;   // x position
  pixelIndex[1] = 650;   // y position
  pixelIndex[2] = 360;   // z position

  for (int i=0; i<10; i++)
    {
      pixelIndex[0] +=   pixelIndex[0];
      ImageType::PixelType   pixelValue = image->GetPixel( pixelIndex );
      std::cout << "pixelValue = " << pixelValue << std::endl;
    }

//  image->SetPixel(   pixelIndex,   pixelValue+1  );

  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 0;
}

******** EndCodeSnippet ********


Oleksandr Dzyubak wrote:
> Hi Luis,
>
> Yes, indeed.
>
> This time I tried both pixel types, "float" and "signed short"
> and compared. The results are the same. It goes well up to
> the slice # 104. From the slice # 105 and on, it's all zeros.
>
> It is some sort of "saturation" or something specific to my machine?
> Did you try on your machine? What's your result?
>
> Regards,
>
> Alex
>
> Luis Ibanez wrote:
>> Hi Alex,
>>
>>  Did you changed to "float" the pixel type of
>>
>>       ImageReadRegionOfInterestWrite.cxx    ?
>>
>>     by default the pixel type is "signed short".
>>
>>    and your image is of pixel type "float".
>>
>>    ITK will perform silent C-Style casting
>>    at the moment of reading the image.
>>
>>
>> Please let us know,
>>
>>      Thanks,
>>
>>           Luis
>>