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

[Stefan Klein]

Hi Dan,

I see your point that the results are counterintuitive if N=1.

But I still think it is a bad idea to just assume some boundary 
condition and increase the valid region. I can imagine that for many 
problems, people would prefer [0,N-1] as a valid range (no need for 
boundary conditions, so simpler/faster code, and no need for any 
arbitrary boundary conditions). By inheriting, more advanced 
interpolators could be implemented, which allow for specifying a 
boundary condition.

Actually, the IsInsideBuffer function definition is always specific for 
an ImageFunction. Only the specific type of ImageFunction (linear 
interpolator, nearest neighbour interpolator, etc) itself knows the 
region on which it can generate a meaningful function value. So, maybe 
it's impossible to define one single implementation that always is 
appropriate.

Regards,
Stefan.

Blezek, Daniel J., Ph.D. wrote:
> 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