[Insight-users] Image index limits
Oleksandr Dzyubak
adzyubak at gmail.com
Thu Jan 7 12:38:00 EST 2010
Hi Bill,
Many thanks for your quick response.
Of course, you are right and my implementation of scanning long OX is
incorrect.
It was meant to be "pixelIndex[0] += 1".
But it was nice that it became "pixelIndex[0] += pixelIndex[0]" since
it enforced me to check the "out-of-range" values. So I updated my code.
1) I changed "pixelIndex[0] += pixelIndex[0]" for "pixelIndex[0] += 1".
2) I scanned along the same line twice just to make sure that Set() and
Get()
methods both work correctly.
a) The entire image volume was filled up homogeneously with I=100;
b) Using Set() method some pixels along the line were assigned I=200.
3) Using Get() method, I got some out-of-range voxel values.
The results are below.
These results surprised me.
In perfect life after clean reboot the system memory should hold zeros.
Well, in real life the system memory holds "nobody knows" what numbers.
However, my out-of-range voxel values are still those I used to initialize
my 1340x1300x720 image. As you pointed out, the first version of my program
used indices which were even much farther away vs I used now but
my out-of-range voxel values were still those I used to initialize my
1340x1300x720 image.
Does the ITK allocate() method reserves much more space as it was requested?
Am I terribly missing something?
Thanks,
Alex
******** Start *******
dzyubak at helium: /Float_1340x1300x720$ ./ImageFill float.hdr
This is initial pixelIndex
pixelIndex = [670, 650, 360]
pixelValue = 100
Now add 100 to element with pixelIndex[0] + 5
pixelIndex = [675, 650, 360]
pixelValue = 200
Roll back to initial index and scan along OX
pixelIndex = [671, 650, 360]
pixelValue = 100
pixelIndex = [672, 650, 360]
pixelValue = 100
pixelIndex = [673, 650, 360]
pixelValue = 100
pixelIndex = [674, 650, 360]
pixelValue = 100
pixelIndex = [675, 650, 360]
pixelValue = 200
pixelIndex = [676, 650, 360]
pixelValue = 100
pixelIndex = [677, 650, 360]
pixelValue = 100
pixelIndex = [678, 650, 360]
pixelValue = 100
pixelIndex = [679, 650, 360]
pixelValue = 100
pixelIndex = [680, 650, 360]
pixelValue = 100
Now pick some out-of-range values
This is out-of-range values
pixelIndex = [2680, 650, 360]
pixelValue = 100
This is out-of-range values
pixelIndex = [5360, 650, 360]
pixelValue = 100
WARNING: In
/mnt/helium/Root_BUILDs/ITK_3_16/InsightToolkit-3.16.0/Code/IO/itkAnalyzeImageIO.cxx,
line 1296
AnalyzeImageIO (0x1da7040): ERROR: Analyze 7.5 File Format Only Allows
RPI, PIR, and RIP Orientation
********* End*****
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.
>
>
> On Tue, Jan 5, 2010 at 4:16 PM, Oleksandr Dzyubak <adzyubak at gmail.com> wrote:
>
>> 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
>>>>
>>>>
>> _____________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html
>>
>> Kitware offers ITK Training Courses, for more information visit:
>> http://www.kitware.com/products/protraining.html
>>
>> Please keep messages on-topic and check the ITK FAQ at:
>> http://www.itk.org/Wiki/ITK_FAQ
>>
>> Follow this link to subscribe/unsubscribe:
>> http://www.itk.org/mailman/listinfo/insight-users
>>
>>
More information about the Insight-users
mailing list