[Insight-users] image coordinates

Simon Warfield simon.warfield at childrens.harvard.edu
Wed Jan 30 10:16:44 EST 2008 

I think it is best to use the center of the voxel to be consistent with 
DICOM.

In DICOM files, the Image Position (Patient) tag, which is frequently 
used to define the origin when reading DICOM files into ITK, gives the 
location of the center of the first voxel.

C.7.6.2.1.1 Image Position And Image Orientation.
The Image Position (0020,0032) specifies the x, y, and z coordinates of the
upper left hand corner of the image; it is the center of the first voxel
transmitted.





> Date: Tue, 29 Jan 2008 18:23:24 -0500 From: "Rupert Brooks" 
> <rupe.brooks at gmail.com> Subject: RE: [Insight-users] ITK Image 
> Coordinates / Problem withPhysicalPoint to Index Conversion To: 
> insight-users at itk.org, andreas.keil at cs.tum.edu Message-ID: 
> <a6cc33f0801291523w291d6805h17a996ebc4f03e68 at mail.gmail.com> 
> Content-Type: text/plain; charset=ISO-8859-1 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

-- 
Simon K. Warfield, Ph.D.        simon.warfield at childrens.harvard.edu
Associate Professor of Radiology,             Harvard Medical School
Director, Computational Radiology Laboratory,
Wolbach 215, Department of Radiology,
Children's Hospital Boston
300 Longwood Avenue Boston MA 02115
Phone: 617-355-4566                                FAX: 617-730-0635

More information about the Insight-users mailing list