[Insight-users] ITK Image Coordinates / Problem withPhysicalPoint
to Index Conversion
Rupert Brooks
rupe.brooks at gmail.com
Tue Jan 29 18:23:24 EST 2008
Hi Andreas,
This is an important issue, particularly for anyone attempting to do
subpixel accuracy registration. I also agree that the coordinate of a
pixel should be the center of that pixel. This is because I tend to
think of images as sampled signals - ie functions convolved with PSF
and then a comb function - so they dont really have width as such. If
the coordinate is of the corner, for many reasons that you have
pointed out, it could lead to strange contradictions or shifts in the
effective image position.
I got very curious (afraid?) of what the ramifications of this might
be... did a little looking.
The sticking point is only the TransformPhysicalPointToIndex function
as far as i can see. Just to see what kind of horrible problems it
would cause if that was changed, i did a quick grep of the ITK 3.4
source code.... This method is really not used much at all - see the
bottom of this email. It would probably be possible to change it to
conform to the definition.
I took a look to see what the NearestNeighborInterpolator does - it
uses the "ConvertContinuousIndexToNearestIndex" function, which, you
guessed it, rounds.
Here it is (from itkImageFunction.h):
/** Convert continuous index to nearest index. */
void ConvertContinuousIndexToNearestIndex( const ContinuousIndexType & cindex,
IndexType & index ) const
{
typedef typename IndexType::IndexValueType ValueType;
for ( unsigned int j = 0; j < ImageDimension; j++ )
{ index[j] = static_cast<ValueType>( vnl_math_rnd( cindex[j] ) ); }
}
Its also interesting to note that TransformPhysicalPointToIndex is the
slowest of the World to Image coordinate conversions, at least on
Intel processors, because truncation is remarkably slower than
rounding on those machines. (See this article for discussion
http://www.mega-nerd.com/FPcast/). But, this is probably starting to
digress from the original topic.
Cheers,
Rupert B.
Heres the list of files that contain references to TransformPhysicalPointToIndex
sutton:/home/scratch/rbrook/InsightToolkit-CVS/src/Insight/Code> grep
-Rl 'TransformPhysicalPointToIndex' .
./Algorithms/itkBioCellularAggregate.txx
./Algorithms/itkEuclideanDistancePointMetric.txx
./BasicFilters/itkBloxBoundaryPointImageToBloxBoundaryProfileImageFilter.txx
./BasicFilters/itkBloxBoundaryPointToCoreAtomImageFilter.txx
./BasicFilters/itkBloxBoundaryProfileImageToBloxCoreAtomImageFilter.txx
./BasicFilters/itkGradientImageToBloxBoundaryPointImageFilter.txx
./BasicFilters/itkPointSetToImageFilter.txx
./BasicFilters/itkPolylineMask2DImageFilter.txx
./BasicFilters/itkPolylineMaskImageFilter.txx
./BasicFilters/itkWarpMeshFilter.txx
./Common/itkImage.h
./Common/itkImageAdaptor.h
./Common/itkImageTransformHelper.h
./Common/itkOrientedImage.h
./Common/itkPhasedArray3DSpecialCoordinatesImage.h
./Common/itkSpecialCoordinatesImage.h
./Common/itkVectorImage.h
./Numerics/FEM/itkFEMSolver.cxx
./Numerics/Statistics/itkJointDomainImageToListAdaptor.h
> Message: 1
> Date: Tue, 29 Jan 2008 18:16:31 +0100
> From: "Andreas Keil" <andreas.keil at cs.tum.edu>
> Subject: RE: [Insight-users] ITK Image Coordinates / Problem
> withPhysicalPoint to Index Conversion
> To: <insight-users at itk.org>
> Message-ID: <003b01c8629a$b13a9030$690a9f83 at 08keil>
> Content-Type: text/plain; charset="us-ascii"
>
> Dear all,
>
> let me elaborate why I am voting for a center-based coordinate system:
>
> Let's consider the simple example of resampling an image. The resample
> filter does the following for every output pixel:
> 1) Calculate the physical coordinate for the given output index.
> 2) Transform the physical coordinate by the given transform.
> 3) Calculate the continuous index for the transformed physical point.
> 4) Linearly interpolate the input value at the continuous index and
> write it to the output.
>
> Now, think about this process with corner-based coordinated:
> 1) Would work fine.
> 2) Transforms are there for converting physical coordinates. Giving them
> corner coordinates would, of course, work, too. However, the
> interpolation in step 4 is done based on a single point representing
> the whole pixel. In general, we don't have an idea on how the box of
> a pixel is transformed. We can just pass a single coordinate through
> through to the interpolator. Keep this in mind for now.
> 3) Works fine.
> 4) Now we have a problem with a corner-based physical point having been
> transformed: The interpolator has to decide, which closest pixels to
> take into account for interpolation (and how to weight them). This
> would be a problem even for a the linear interpolator. How should
> it know, which 4 pixels (in the 2D case) to take? The problem is,
> that the interpolator has to estimate the location of the backward
> warped input pixel. This is never really possible due to the unknown
> transformation. But the guess of a center-based rectangle is better
> than guessing a corner-based rectangle (where you don't know, which
> of the corners is given).
> A numerical example (sorry, I was too lazy to paint): Let's resample
> a 2D image with an input and output spacing of 1mm x 1mm. We want to
> interpolate the output pixel at (0mm..1mm,0mm..1mm) which has the
> corner coordinates (0mm,0mm).
> A) For a translation of (0.1mm,0.2mm), we get an input coordinate
> at (0.1mm,0.2mm) and now have to decide, which image range of
> 2x2 pixels to use for interpolation. Since we translated a bit
> in positive x and y direction, we should use the pixels at
> (0mm..2mm,0mm..2mm). Alright. So the rule should be to round the
> physical coordinate and this yields the lower left corner of our
> 2x2 pixels interpolation region.
> B) Let's try a different transformation: A 180 deg. rotation around
> (0mm,0mm) and a subsequent translation of (0.1mm,0.2mm). We again
> compute (0.1mm,0.2mm) as input coordinate. So let's follow the
> rule and use again (0mm..2mm,0mm..2mm) as interpolation region.
> Ups! We should have used (-1mm..1mm,-1mm..1mm) instead! (You can
> see this by also transforming the "upper-right" corner output
> pixel which is (1mm,1mm) giving (-0.9mm,-0.8mm) as input corner.
> Therfore, the output pixel lies at (-0.9mm..0.1mm,-0.8mm..0.2mm).
> This input pixel overlaps the 4 input pixels in the region
> (-1mm..1mm,-1mm..1mm) as stated above.)
>
> Conclusion: The transformed corner coordinate in reality is no
> longer a "lower-left" corner but how should the interpolater have
> guessed that? The only remedy would be that every class using a
> transformation would not only transform one single coordinate per
> pixel but all 4 (in 2D) or 8 (in 3D) corners. Nobody wants this.
> Guessing an input region after having backward warped a pixel's
> center is much easier and faster. (Although I would opt for taking
> circles instead of rectangles here, but that's another issue.)
>
> Hopefully, I made my point clear that a pixel's center coordinate is the
> more useful representation for it than its corner coordinate. I am pretty
> sure, the original developers or writer of the docs had this point in mind
> when the wrote the Software Guide. And since the implementation is flawed
> anyways, and since a lot of classes could not be adopted to corner-based
> coordinates, I still strongly vote for center-based ones. This "physical
> point of view" is also more in the style of ITK than the "computer science
> point of view". ITK developers should not think in pixels but rather in
> millimeters.
>
> Thanks for reading. Waiting for your replies,
> Andreas.
>
>
>
> -----Original Message-----
> From: Maarten Nieber
> Sent: Monday, January 28, 2008 14:45
> To: insight-users at itk.org
> Subject: Re: [Insight-users] ITK Image Coordinates / Problem
> withPhysicalPoint to Index Conversion
>
> Hi Andreas,
>
> I agree with you that this issue is of very high importance!
>
> -- quote from your previous post --
> I think that the personal preference
> for one of the two options we have basically depends on whether one is
> more interested in the overall dimensions of a dataset or whether one has
> to do voxel-wise operations depending on the physical coordinate.
> -- end quote --
>
> Can you explain why you think - for doing voxel-wise operations - it would
> be better if the image origin co-incides with the center of a pixel?
>
> Best regards,
> Maarten
>
>
> ----- Original Message ----
> From: Andreas Keil
> To: insight-users at itk.org
> Sent: Tuesday, January 8, 2008 10:37:32 AM
> Subject: RE: [Insight-users] ITK Image Coordinates / Problem with
> PhysicalPoint to Index Conversion
>
> Dear all,
>
> I would like to try a restart of the discussion related to the definition
> of physical coordinates if ITK images:
>
> As mentioned in my earlier post (see below), the documentation and
> implementation differ in this point. The two proposed solutions are either
> fixing the documentation (ITK Software Guide, p. 40) and defining the
> origin to lie in the "lower-left(-back)" corner of the pixel with index
> 0/0(/0) or fixing the implementation of the respective methods
> (IndexToPhysicalPoint, PhysicalPointToIndex, etc.) and defining the origin
> to lie in the center of this pixel.
>
> So far, only Maarten and I were discussing about this, with him favoring
> to fix the docs and me favoring to fix the implementation. Another
> argument I found for my preference is that the origin would not have to
> change when downsampling an image. I think that the personal preference
> for one of the two options we have basically depends on whether one is
> more interested in the overall dimensions of a dataset or whether one has
> to do voxel-wise operations depending on the physical coordinate.
>
> --> Therefore, I strongly ask other ITK users and developers to take part
> in this discussion so that we can make a decision. This problem is really
> of high importance since it is related to the core of the lib (namely the
> ImageBase class) and it has a big influence on reconstruction problems. As
> soon as we have come to a conclusion, I would file a bug report with the
> proposed solution.
>
> Thank you,
> Andreas.
>
>
>
> -----Original Message-----
> From: Andreas Keil
> Sent: Tuesday, December 18, 2007 13:30
> To: insight-users at itk.org
> Subject: RE: [Insight-users] ITK Image Coordinates / Problem with
> PhysicalPoint to Index Conversion
>
> Hi Maarten,
>
> thank you for your reply (which other list subscribers may find below). My
> initial posting was biased towards changing the ITK implementation rather
> than the documentation for the following reasons:
>
> Working with the physical coordinates of a voxel usually means that one
> needs a single coordinate representation of a voxel. Algorithms relying
> one this single physical coordinate would usually prefer the center of a
> voxel for their space-related computations.
>
> Moreover, ITK images only reflect the spacing between voxels which is not
> necessarily equal to the width of a voxel. I am pretty sure that there are
> cases where the voxel width is smaller than the spacing (in CT for
> example). In this case, the term "spacing" would also be ambigious if it
> would not refer to the distance between adjacent voxels' centers.
>
> What do others think about this? I would appreciate a clearification as
> the definitions of image origin and spacing really matter for my
> application in 3D reconstruction.
>
> Thank you,
> Andreas.
>
>
>
> -----Original Message-----
> From: Maarten Nieber
> Sent: Tuesday, December 18, 2007 10:11
> To: Andreas Keil
> Subject: Re: [Insight-users] ITK Image Coordinates / Problem with Physical
> Point to Index Conversion
>
> Hi Andreas,
>
> I agree with you that according to the software guide, mapping (0.6, 0.6)
> should yield an index of (0,0).
> On the other hand, I think that the definition of origin in the software
> guide is not intuitive.
> It says "the image origin is associated with the co-ordinates of the first
> pixel in the image". Probably, it would be more accurate to say that the
> image origin is associated with the >center< of the first pixel in the
> image. However, by using this definition, the extent of an image with
> origin (0,0) contains negative co-ordinates. I would find it more
> intuitive if the extent of such an image is something like (0,0) - (size1,
> size2). This would be the case if the origin of the image is associated
> with the bottom-left corner of the first pixel.
>
> If in the itk code, the origin of an itk image co-incides with the bottom
> left corner of the first pixel, then I would prefer to change the software
> guide and not the code.
>
> Best regards
> Maarten Nieber
>
>
>
> ----- Original Message ----
> From: Andreas Keil
> To: insight-users at itk.org
> Sent: Friday, December 14, 2007 6:10:25 PM
> Subject: [Insight-users] ITK Image Coordinates / Problem with Physical
> Point to Index Conversion
>
> Hi,
>
> I think I have discovered an inconsistency between the ITK Software Guide
> (p.40) and the implementation of itk::Image (all the methods taking
> physical points / continuous indices as arguments):
>
> The trunctation of continuous index coordinates to integers does not yield
> the correct pixel coordinates as expected by looking at the definition in
> the software guide.
>
> A simple example is:
> Image spacing: 1mm
> Image origin: 0mm/0mm
>
> The pixel with index (1/1) would (according to the software guide) have
> the following extents:
> 0.5mm/0.5mm to 1.5mm/1.5mm
>
> However, the physical point 0.6mm/0.6mm gets mapped to the index 0/0 by
> the method TransformPhysicalPointToIndex.
>
> The solution would be to check those conversion methods as well as others
> (like IsInside) and replace integer truncations with rounding.
>
> If my reasoning is correct, I'll file a bug report. However, I'd like to
> have some confirmation first.
>
> Thanks,
> Andreas.
>
--
--------------------------------------------------------------
Rupert Brooks
McGill Centre for Intelligent Machines (www.cim.mcgill.ca)
Ph.D Student, Electrical and Computer Engineering
http://www.cyberus.ca/~rbrooks
More information about the Insight-users
mailing list