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

[Paul Yushkevich]

Hi Luis,

I am sorry, I have not communicated very well what I am trying to do.

I have an NxMxK image, which I want to blur with an isotropic Gaussian.  
However, I want to see only one slice from the blurred image, i.e., I am 
interested in the region R of the output image, located at index 0,X,0 
and having size Nx1xK.  The computation of R in the output image 
requires convolution over the slab-shaped region R* of the input image, 
with R* having index 0,X-radius,0 and size Nx(1+2*radius)xK, where 
'radius' is the radius of the Gaussian operator.  This is quite 
different from blurring a NxK two-dimensional slice, which ignores the 
adjacent slices.

During PropagateRequestedRegion, DiscreteGaussianImageFilter (DGIF), 
correctly computes R* given R.  However, in GenerateData, it does not 
seem to propagate the region correctly between different stages of 
filtering, hence the exception.

Thank you!

Paul.




Luis Ibanez wrote:

> 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.
>>>>
>>>
>>>
>>>
>>>
>>
>>
>
>
>
>


-- 
--------------------------------
Paul A. Yushkevich, Ph.D.
President, Cognitica Corporation

17 Flemington Rd
Chapel Hill, NC 27517
Tel: 1-919-929-7652
--------------------------------