[Insight-developers] Bug in ImageFunction causing incorrect results from ResampleImageFilter

[Blezek, Daniel J., Ph.D.]

Hi Stefan,

  You are indeed correct, the range is closed: [0, N-1].  But I disagree
that this is correct behavior.  Consider a single slice.  It is
difficult to have an index inside it because the valid range is: N = 1,
[0, N-1] == [0, 0].  Even if my slices are 5mm thick, I am treating them
as infinitely thin.  In my opinion, this is incorrect.  I should
reasonably expect to resample a 5mm slice into 5 1mm slices.  I may not
like the results depending on the boundary conditions, but I should get
"something" out of the resampler.  Currently, I would get the first
slice, and 4 blank slices.

  Having our slices go from -0.5 to N-0.5 would match DICOM's definition
of a slice.  This would be very sensible.

OK, forgive the ASCII art...  For a two slice volume, interpolated to 10
slices

Z           0-----1-----2-----3
ITK         [+....00000]
DICOM    [..+....+..]
Dan         [+....+....]

"." are interpolated slices, "+" are slice centers, "0" is where ITK
decided the slice was out side the buffer.  ITK essentially only gives
me 5 slices (of data) when I should have 10.

  99% of the time we don't care about this issue.  When it shows up
(like in registration and resampling), it's ugly and we really care.

  Another issue is backwards compatibility.  Do we dare not change this
fundamental behavior for fear of breaking existing code?  If indeed the
community decides to change ITK's definitions, how do we do it?  Perhaps
Bill Lorensen would care to comment?

Regards,
-dan

-----Original Message-----
From: Stefan Klein [mailto:s.klein at erasmusmc.nl] 
Sent: Monday, June 02, 2008 11:39 AM
To: Blezek, Daniel J., Ph.D.
Cc: insight-developers at itk.org
Subject: Re: [Insight-developers] Bug in ImageFunction causing incorrect
results from ResampleImageFilter

Hi Dan,

I'm not sure, but: it seems that with the current implementation the
valid range is:

[0, N-1]

instead of [0,N-1) as you said.

(if index==m_EndContinuousIndex the method returns true).

In my opinion, the current implementation is right. For example, I see
no meaningful default implementation for a linear interpolator at
indices > N-1. (unless some boundary condition has been set).

Moreover, related to the discussion on the voxel centre (IndexToPoint
conversion): wouldn't it be more logical to define the valid region as
going from -0.5 to N-0.5?
(but that would cause troubles again for some interpolators, so i would
vote to keep it like it is right now)

Kind regards!
Stefan


Blezek, Daniel J., Ph.D. wrote:
> I have some fairly thick CT slices (5mm) that I would like to resample

> to be isotropic.  A straightforward use of ResampleImageFilter leaves 
> 4 blank slices in the output volume.  I've tracked this down to a 
> bug(?) in ImageFunction.  This code is incorrect IMHO:
> 
>   virtual bool IsInsideBuffer( const ContinuousIndexType & index )
const
>     { 
>       for ( unsigned int j = 0; j < ImageDimension; j++ )
>         {
>         if ( index[j] < m_StartContinuousIndex[j] ) { return false; };
>         if ( index[j] > m_EndContinuousIndex[j] ) { return false; };
>         }
>       return true;
>     }
> 
> (Dealing with Z in 3D) This code assumes the last slice has zero 
> thickness.  For N slices, the range where IsInsideBuffer returns true 
> is [0, N-1).  This is incorrect, the correct result should be [0, N).

> [0,
> N) gives the last slice it's proper thickness.  The windowed Sinc 
> interpolator properly handles the last slice, i.e. [N-1, N), but I 
> have not checked others.  A proper implementation would look like this

> (note the vcl_floor calls):
> 
>           if ( vcl_floor(index[j]) < m_StartContinuousIndex[j] ) { 
> return false; };
>           if ( vcl_floor(index[j]) > m_EndContinuousIndex[j] ) { 
> return false; };
> 
> 
> I've run all the tests on a fresh check out with the above change.  As

> may be expected, almost all of the registration tests fail, while most

> of the resample tests pass.  Surprisingly, 
> ResampleVolumesToBeIsotropicTest passed...
> 
> Any comments on this bug?  Should the tests be made to accommodate the

> change, or should I override the method in a custom interpolator?
> 
> Regards,
> -dan
> 
> 
> 
> Daniel Blezek, PhD
> Medical Imaging Informatics Innovation Center Enterprise Imaging 
> Systems Unit
> 
> P 127 or (77) 8 8886
> T 507 538 8886
> E blezek.daniel at mayo.edu
> 
> Mayo Clinic
> 200 First St. S.W.
> Harwick SL-44
> Rochester, MN 55905
> mayoclinic.org
> 
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk.org
> http://www.itk.org/mailman/listinfo/insight-developers