[Insight-developers] A possible bug in neighborhood interator
or filter?
Joshua Cates
cates at sci . utah . edu
Wed, 28 May 2003 13:24:35 -0600 (MDT)
Hi Paul,
Haven't yet had time to check this out today, but I'll have a look later
toady and see if I can understand what the problem is.
Perhaps there is a problem with GenerateInputRequestedRegion.
Josh.
______________________________
Josh Cates
School of Computer Science
University of Utah
Email: cates at sci . utah . edu
Phone: (801) 587-7697
URL: http://www . sci . utah . edu/~cates
On Wed, 28 May 2003, Paul Yushkevich wrote:
> 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
> --------------------------------
>
>
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at public . kitware . com
> http://public . kitware . com/mailman/listinfo/insight-developers
>