[Insight-developers] Bug in ImageFunction causing incorrect results from ResampleImageFilter
Blezek, Daniel J., Ph.D.
Blezek.Daniel at mayo.edu
Wed Jun 18 13:45:35 EDT 2008
Stefan,
At the last ITK tcon (just Jim Miller and I), we decided to have a go
and making points and indices uniform across ITK. Simon Warfield and I
will be making a branch in a month or so when travel schedules permit.
The goal is to produce green dashboards with a consistent view of points
and indices. This will require crawling through much of the
registration framework, and adapting each ImageFunction to understand
and appreciate the conventions.
I would welcome any developers interested in this issue.
Regards,
-dan
-----Original Message-----
From: Stefan Klein [mailto:s.klein at erasmusmc.nl]
Sent: Tuesday, June 03, 2008 12:30 PM
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 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
More information about the Insight-developers
mailing list