[Insight-developers] A possible bug in neighborhood interator or filter?

[Luis Ibanez]

Yeap, my mistake,

I missed the line with the initialization of the region

 > ImageRegion<3> region = image->GetLargestPossibleRegion();

Sorry about that,


It is likely that the DiscreteGaussian filter will attempt
to blurr also along the dimension that is 1 pixel wide.
Since it probably need a minimum number of neighbors, that
may be triggering the error.

The RecursiveGaussian, for example, will need at least 5
pixels along each dimension.

In your case you could use two RecursiveGaussian filters
and explicitly say that you only want to filter along
dimension 0 and 1. In that way, the fact that the
third dimension size is 1 will not affect the filtering.
It will be just like filtering in 2D



   Luis


-----------------
Paul Yushkevich wrote:
> Hi Luis,
> 
> I believe that the line
> ImageRegion<3> region = image->GetLargestPossibleRegion();
> 
> initializes all the components of the region.  The subsequent calls to 
> region.SetSize() and region.SetIndex() shrink the region in one of the 
> dimensions to be one pixel thin.  When this dimension is in the X axis, 
> everything works, but if it's in a different axis, then the code crashes.
> 
> What I'm trying to do is to perform only as much blurring as is 
> necessary to produce one slice in the output image.
> 
> I think that there is some sort of a bug here.  The DiscreteGaussian 
> filter chains three Neighborhood filters together, and I think that the 
> regions produced by these filters might not be being calculated 
> correctly.  That's my intuition after some debugging.
> 
> Thank you!
> 
> Paul
> 
> Luis Ibanez wrote:
> 
>> Hi Paul,
>>
>> The methods
>>
>>   SetSize(  component, value  )
>>   SetIndex( component, value  )
>>
>> set only one component of the Size/Index.
>>
>> In the ImageRegion() default constructor
>> both m_Size and m_Index are filled up with
>> zeros.
>>
>> In your example, since you are initializing
>> only the second component of the size, the
>> first one is still zero. So, you seem to be
>> asking for a region of size:
>>
>>      size =  (0,1)
>>
>>      0 along X
>>      1 along Y
>>
>> Most iterators will be confused with this size.
>> They expect to have a value of at least 1 along
>> each dimension.
>>
>> It is actually surprising that it works with
>>
>>    region.SetSize( 0, 1 );
>>
>>
>> You may want to call something like
>>
>>    region.SetSize(0, sizeX );
>>    region.SetSize(1, sizeY );
>>    region.SetSize(2, sizeZ );
>>
>>    region.SetIndex(0, indexX );
>>    region.SetIndex(1, indexY );
>>    region.SetIndex(2, indexZ );
>>
>>
>> Since your image is 3D
>>
>>
>>
>>    Luis
>>
>>
>>
>> -----------------------
>> Paul Yushkevich wrote:
>>
>>> The following code throws an exception in line 36 of 
>>> itkZeroFluxNeumannBoundaryCondition.txx when filter->Update() is called.
>>>
>>> #include "itkImageFileReader.h"
>>> #include "itkImage.h"
>>> #include "itkDiscreteGaussianImageFilter.h"
>>>
>>> using namespace itk;
>>>
>>> void main(int argc, char *argv[])
>>> {
>>>
>>>  // Typedefs
>>>  typedef Image<float,3> ImageType;
>>>  typedef ImageFileReader<ImageType> ReaderType;
>>>  typedef DiscreteGaussianImageFilter<ImageType,ImageType> FilterType;
>>>
>>>  // Load image
>>>  ReaderType::Pointer reader = ReaderType::New();
>>>  reader->SetFileName("MRIcrop-orig.gipl");
>>>  reader->Update();
>>>  ImageType::Pointer image = reader->GetOutput();
>>>
>>>  // Create a requested region
>>>  ImageRegion<3> region = image->GetLargestPossibleRegion();
>>>  region.SetIndex(1,55);
>>>  region.SetSize(1,1);
>>>
>>>  // Create filter
>>>  FilterType::Pointer filter = FilterType::New();
>>>  filter->SetVariance(1.0f);
>>>  filter->SetInput(image);
>>>  filter->GetOutput()->SetRequestedRegion(region);
>>>  filter->Update();
>>> }
>>>
>>> I didn't attach the image because it's a bit large for posting to the 
>>> list.  However, it is of sufficient size (78x110x64). It seems that 
>>> the neighborhood iterator jumps outside of bounds, causing the 
>>> exception.  The funny thing is that it works with
>>>
>>>  region.SetIndex(0,55);
>>>  region.SetSize(0,1);
>>>
>>> but not with
>>>
>>>  region.SetIndex(2,55);
>>>  region.SetSize(2,1);
>>>
>>> Paul.
>>>
>>
>>
>>
>>
> 
>