MantisBT - ITK
View Issue Details
IDProjectCategoryView StatusDate SubmittedLast Update
0006558ITKpublic2008-03-07 05:202010-10-21 11:15
ReporterAndreas Keil 
Assigned ToMichelKitware 
PrioritynormalSeveritymajorReproducibilityalways
StatusclosedResolutionfixed 
PlatformOSOS Version
Product Version 
Target VersionFixed in Version 
Resolution Date
Sprint
Sprint Statusbacklog
Summary0006558: Physical coordinates of a pixel - Severe inconsistency and bug in ImageBase
DescriptionProblems:
1) The implementation of the conversion between physical coordinates and indexes differs from the documentation. The former uses corner-based coordinates whereas the latter describes center-based coordinates.
2) The implementation contains a severe bug since it uses truncation (instead of either flooring or rounding). This yields wrong indexes for negative coordinates.

Severity: This bug affects every application/filter using physical coordinates of ITK images (esp. registration and reconstruction). It has to be fixed immediately.

Solution: Decide about whether to use center-based or corner-based coordinates and correct the implementation (and perhaps documentation).

Summary:
- The code and documentation are inconsistent and one or the other has
  to be fixed.
- Backwards compatibility cannot be achieved - either the docs or the
  code have to be changed.
- ImageBase definitly contains a bug which has to be fixed since it
  uses truncation instead of floor().
- The two possibilities received the following supporters:
  + center-based: Andreas, Steve, Rupert, Edson, Simon
  + corner-based: Maarten, Eliana

Arguments for center-based:
+ More intuitive. People from image processing interprete voxels as smapled data which is best represented at the center of the PSF.
+ If a single coordinate is required for representing a voxel in an algorithm, the center is usually the better choice.
+ Spacing is not necessarily equal to the extent (not stored in ITK images) of a voxel.
+ Origin would not have to change when downsampling an image.
+ Natural expression of formulas for interpolation.
+ Lower-left(-back) coordinates do not remain lower-left(-back) coordinates after transformations like a 180-degrees rotation. Transformations mandatorily have to return center-based coordinates in order to work.
+ Voxels do not have to be interpreted as rectangular objects. They could be circles, gaussians, ...
+ Consistency with DICOM.

Arguments for corner-based:
+ More intuituve. People from graphics/rendering interprete voxels as cubes of color which are easier to handle when represented using corners.
+ There may be filters assuming corner-based physical coordinates.

Non-countig argument for corner-based:
- Backwards compatibility of the code.
  Backwards compatibility cannot be kept anyways since truncation is
  wrong for negative physical coordinates and has to be replaced by
  rounding in both proposed solutinos (round or floor). Negative
  coordinates are also not exotic, they appear as soon as one centers
  the coordinate frame in a volume, e.g.

Please also see the implications of changing the implementation which Rupert mentions in his last post (see additional information).
Steps To Reproduce
Additional InformationIn the following sections, I gathered all the messages on this issue posted on the ITK users list:
2007-12-14 to 2007-12-19 (4 messages)
2008-01-08 to 2008-01-09 (3 messages)
2008-01-29 to 2008-01-30 (6 messages)

---------------------------------------------------------------

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.

---------------------------------------------------------------

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

---------------------------------------------------------------

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.

---------------------------------------------------------------

Hi Andreas,

with respect to the "co-ordinate representation of a voxel", I'm not completely sure what you mean, but I think that for the conversion between index and physical point, the following functions are needed:

- given a physical point, there should be a function that returns the index of the pixel that contains this physical point.

- given a pixel index, there should be a function that returns the physical point at the center of this pixel.

In my opinion, for an 2x2 image with spacing: 1mm and origin: 0mm/0mm, the function PhysicalPointToIndex( 0.6, 0.6 ) should return the index (0,0). The function IndexToPhysicalPointAtPixelCenter( 0, 0 ) should return the point (0.5, 0.5). The extent of the image should be (0,0) - (2,2). (Of course, the function arguments should be itkPoint, I simplified the notation). Again, this is just my preference.

Best regards
Maarten

---------------------------------------------------------------

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.

---------------------------------------------------------------

Hello Andreas,

I think you're right: the docs are correct but the code is buggy.


On Tue, Jan 08, 2008 at 10:37:32AM +0100, Andreas Keil wrote:

> 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.

Right. I find the centre-based coordinate convention to be much more
natural; e.g. when writing out formulas for interpolation. If you're
taking votes, put me down for the centre.

I think the centre-based coordinate convention is more common with
image processing folk, whereas the corner-based convention is more
common in graphics circles. For what it's worth, VTK's image data is
said to use the centre-based coordinate system [1]. On the other
hand, VTK uses corner convention for screen pixel coordinates, and I
suspect this confusion underlies the problems in VTK's coordinate
conversion routines [2].


[1] http://public.kitware.com/pipermail/vtk-developers/2005-November/003825.html [^]
[2] http://www.vtk.org/Bug/view.php?id=5111 [^]

[Message by Steve M. Robbins]

---------------------------------------------------------------

Hi,

Independent of using center-based or corner-based data location, the actual implementation of PhysicalPoint To Index uses TRUNCATION as a rounding mode, and this is plain wrong (it is implementing neither of the solutions). If you allow a negative Origin (e.g., (-0.4, -0.4)), then you could have one pixel with coordinate at -0.4 mapped to 0.0, and then the next pixel with coordinate at 0.6 mapped to 0.0 again. I suggest using Floor (corner-based approach) or proper Rounding (center-based approach).

Whatever solution is chosen, I think it should also be reflected in itk::OrientedImage. Whatever solution is chosen, the data location can be changed from center to corner and vice versa through a shift of the Origin (by Spacing/2 for itk::Image or by Direction*Spacing*(0.5,0.5,0.5) for itk::OrientedImage); because of this, I would suggest to take the approach that breaks less code.
 
[Message by Edson Tadeu]

---------------------------------------------------------------

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

---------------------------------------------------------------

Hello Maarten and Andreas, and it-users
I would vote to better change the documentation than the code. I don't
know for sure, but I suppose there are some fitlers that already assume
the origin to be set in the lower-left-back (llb) corner in the image.
Eliana

---------------------------------------------------------------

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.

---------------------------------------------------------------

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

---------------------------------------------------------------

I completely agree.. the "pixel" (or voxel) do not have even to be rectangular... or the same size of the Spacing. They could be tiny circles... or gaussians... or even tiny dots (dirac delta). E.g., in visualization, when you consider the common linear interpolation, they have a width of twice the Spacing (overlapping each other); considering bi-cubic or other kinds of interpolation (Lanczos, etc.), they are even larger.

I vote for the center-based location approach (rounding mode), but an easy way to switch to vertex-based should be available (maybe through a flag bit and a shift in the Origin), because there may be some applications where vertex-based information is more natural (finite elements, PDE, etc.).

As a note, there should always be an unambiguous (non-overlapping) pixel domain [e.g., equivalent to the Voronoi diagram]. Even in the center-based location, there are still an issue to be resolved: the domain of the boundary pixels.

For example, for a 3x3 image, origin (0.0, 0.0), spacing (1.0, 1.0), the domain of each pixel it would be something like this:

Pix [0,0]: [-Inf, -Inf] <= x < [0.5, 0.5]
Pix [1,0]: [-Inf, 0.5] <= x < [0.5, 1.5]
Pix [2,0]: [-Inf, 1.5] <= x < [0.5, Inf]
Pix [0,1]: [0.5, -Inf] <= x < [1.5, 0.5]
Pix [1,1]: [0.5, 0.5] <= x < [1.5, 1.5]
Pix [2,1]: [0.5, 1.5] <= x < [1.5, Inf]
Pix [0,2]: [1.5, -Inf] <= x < [Inf, 0.5]
Pix [1,2]: [1.5, 0.5] <= x < [Inf, 1.5]
Pix [2,2]: [1.5, 1.5] <= x < [Inf, Inf]

Note that there are three options to the domain of the boundary pixels:
1) Extend them to infinity (as I did)
2) Extend them so that they have the same "size" of the inner pixels. In this case, the "image domain" would be (-0.5,-0.5) <= x < (2.5, 2.5)
3) Extend them to half (or 1/4) the size of the inner pixels, towards them. In this case, the "image domain" would be (0.0, 0.0) <= x < (2.0, 2.0), but the boundary pixels would be at the *vertices* of their domains.

Here is a fixed-width font diagram exemplifying it:


                | |
                | |
   (-0.5,2.5) | | (2.5,2.5)
        +.......+.......+.......+
        . | | .
        . [0,2].|.[1,2].|.[2,2] .
        . . | | . .
   -----+-------+-------+-------+-----
        . . | | . .
        . [0,1] | [1,1] | [2,1] .
        . . | | . .
   -----+-------+-------+-------+-----
        . . | | . .
        . [0,0].|.[1,0].|.[2,0] .
        . | | .
        +.......+.......+.......+
   (-0.5,-0.5) | | (2.5,-0.5)
                | |
                | |


Options 1 seems more natural to me, and it is explicit that this "image domain" (Infinity) does not correspond to the real "image metrics" (0.0, 0.0) to (3.0, 3.0).
Option 2 seems natural too, but note that *many* people will assume that the "image domain" corresponds to the "image metrics", i.e., the image domain to start at the Origin: (0.0, 0.0) to (3.0, 3.0) in this case. A solution would be to shift the Origin by (-0.5, -0.5), but then, again, we would be moving the information location from the center to the vertices ;).

Edson

---------------------------------------------------------------

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.

[Message by Simon Warfield]

---------------------------------------------------------------
TagsNo tags attached.
Relationships
has duplicate 0000738closed Cory W Quammen Physical Point to Index : floor/ceil/round 
Attached Filespatch itk-bug6558.patch (8,647) 2008-11-27 10:54
https://public.kitware.com/Bug/file/1876/itk-bug6558.patch
patch CenteredPixelCoordinatesMay-5-2009.patch (24,119) 2009-05-05 16:26
https://public.kitware.com/Bug/file/2208/CenteredPixelCoordinatesMay-5-2009.patch
patch CenteredPixelCoordinatesMay-5-2009-B.patch (27,387) 2009-05-05 17:57
https://public.kitware.com/Bug/file/2209/CenteredPixelCoordinatesMay-5-2009-B.patch
patch PortableRound-May-11-2009.patch (40,668) 2009-05-11 18:41
https://public.kitware.com/Bug/file/2231/PortableRound-May-11-2009.patch
patch itk-portable-round-2009-05-12.patch (44,535) 2009-05-12 08:46
https://public.kitware.com/Bug/file/2236/itk-portable-round-2009-05-12.patch
patch itk-portable-round-2009-05-15.patch (45,362) 2009-05-15 12:50
https://public.kitware.com/Bug/file/2254/itk-portable-round-2009-05-15.patch
patch itk-portable-round-2009-05-16.patch (56,066) 2009-05-16 12:07
https://public.kitware.com/Bug/file/2255/itk-portable-round-2009-05-16.patch
patch itk-portable-round-2009-05-26-bis.patch (34,310) 2009-05-26 08:48
https://public.kitware.com/Bug/file/2269/itk-portable-round-2009-05-26-bis.patch
patch itk-templatedmathround-2009-07-27.patch (42,410) 2009-07-27 12:48
https://public.kitware.com/Bug/file/2383/itk-templatedmathround-2009-07-27.patch
patch itk-templatedmathround-2009-07-30.patch (45,131) 2009-07-30 08:41
https://public.kitware.com/Bug/file/2390/itk-templatedmathround-2009-07-30.patch
Issue History
Date ModifiedUsernameFieldChange
2008-03-07 05:20Andreas KeilNew Issue
2008-05-14 03:25mnieberNote Added: 0011885
2008-05-14 04:57Volker DaumNote Added: 0011890
2008-05-14 07:40Andreas KeilNote Added: 0011891
2008-05-16 09:43Andreas KeilNote Added: 0011961
2008-11-27 10:54Tom VercauterenFile Added: itk-bug6558.patch
2008-11-27 10:58Tom VercauterenNote Added: 0014214
2009-04-20 16:41Luis IbanezStatusnew => assigned
2009-04-20 16:41Luis IbanezAssigned To => MichelKitware
2009-05-05 16:26Luis IbanezFile Added: CenteredPixelCoordinatesMay-5-2009.patch
2009-05-05 16:27Luis IbanezNote Added: 0016302
2009-05-05 17:57Luis IbanezFile Added: CenteredPixelCoordinatesMay-5-2009-B.patch
2009-05-05 17:59Luis IbanezNote Added: 0016304
2009-05-11 18:41Luis IbanezFile Added: PortableRound-May-11-2009.patch
2009-05-12 08:46Tom VercauterenFile Added: itk-portable-round-2009-05-12.patch
2009-05-15 12:50Tom VercauterenFile Added: itk-portable-round-2009-05-15.patch
2009-05-16 12:06Luis IbanezNote Added: 0016502
2009-05-16 12:07Luis IbanezFile Added: itk-portable-round-2009-05-16.patch
2009-05-26 08:48Tom VercauterenFile Added: itk-portable-round-2009-05-26-bis.patch
2009-07-27 12:48Tom VercauterenFile Added: itk-templatedmathround-2009-07-27.patch
2009-07-30 08:41Tom VercauterenFile Added: itk-templatedmathround-2009-07-30.patch
2009-10-21 16:57Bradley LowekampNote Added: 0018149
2010-10-21 11:15Hans JohnsonSprint Status => backlog
2010-10-21 11:15Hans JohnsonNote Added: 0022566
2010-10-21 11:15Hans JohnsonStatusassigned => closed
2010-10-21 11:15Hans JohnsonResolutionopen => fixed
2010-11-04 23:06Cory W QuammenRelationship addedhas duplicate 0000738

Notes
(0011885)
mnieber   
2008-05-14 03:25   
I noticed in the overview that my last post has been missed (no offense taken :-). It is in reply to the post starting with "let me elaborate why I am voting for a center-based coordinate system".

Best regards,
Maarten

Hi Andreas,

it seems to me that the problems you indicate are related to the question: what is the co-ordinate of a pixel?
The way I see it, a pixel has no co-ordinate (but it has an area).
I could rephrase your step1 one as follows:

1. For the given output index, calculate the co-ordinate >at the center< of the corresponding output pixel.

Ideally, this should be performed by a helper function: IndexToPhysicalCoordinateAtPixelCenter(). Note that the output of this function depends on whether you have the origin of the image at the bottom left or center of a pixel. For a 'bottom-left' representation, the function IndexToPhysicalCoordinateAtPixelCenter(0) would return 0.5 * imageSpacing, whereas for a 'center' representation, it would be zero.
However, the writer of the resample algorithm will not be concerned with this, because in both cases IndexToPhysicalCoordinateAtPixelCenter() will return the co-ordinate at the center of the pixel, and the algorithm will perform correctly.

Best regards, and looking forward to your replies too,
Maarten
(0011890)
Volker Daum   
2008-05-14 04:57   
To add to Ruperts observations:
itk::NearestNeighborInterpolateImageFunction uses rounding, therefore conforms to the documentation.
itk::LinearInterpolateImageFunction uses floor (or something like it) to generate the base index and then generates the remaining needed indices from this (by adding +1 along the different dimensions). This also conforms to the documentation.
itk::BSplineInterpolateImageFunction seems to (as far as I understood by looking at it just shortly) conform to the documentation as well.

The interpolate image functions (as a quick grep shows) seem to be more widely used in the ITK source than the Transform method, which makes sense, since anyone concerned about sub-pixel accuracy will probably want to use interpolation of some kind. They are also used in the ResampleImageFilter, which is probably used in most user code that does resampling. Therefore I would urge to fix the bug in the TransformPhysicalPointToIndex method to conform with the documentation and the rest of the interpolators, as this is the solution that will affect the least other code. (I suspect this is also true for most user code. It definitly is for mine).
(0011891)
Andreas Keil   
2008-05-14 07:40   
Another argument for center-based physical coordinates:

DICOM explicitly specifies that coordinates are center-based, see e.g. the Image Position Patient tag (0020,0032) (explained in section C.7.6.2.1.1) and the section on deformable registration C.20.3.1.

ITK is not DICOM but DICOM is the standard for medical imaging and I bet that DICOM tags are currently not converted to corner-based coordinates when reading DICOM with ITK (gdcm).
(0011961)
Andreas Keil   
2008-05-16 09:43   
The discussion on this topic has moved to the following ITK Wiki page:
http://www.itk.org/Wiki/Proposals:Refactoring_Index_Point_Coordinate_System [^]

I'll try to compile the information from this bug report into the Wiki page. Please feel free, to add your own comments there, too.
(0014214)
Tom Vercauteren   
2008-11-27 10:58   
I have upload a preliminary patch. It allows the use of centered pixels in a consistent manner but only applies to the very core of ITK. It can be switched on or off from CMake.

More work is required to make all unit tests pass.
(0016302)
Luis Ibanez   
2009-05-05 16:27   
Michel and Luis have continued fixing secondary effects of the original patch.

The current cumulative patch has been uploaded in the file:
CenteredPixelCoordinatesMay-5-2009.patch
(0016304)
Luis Ibanez   
2009-05-05 17:59   
An updated patch has been uploaded with the name

CenteredPixelCoordinatesMay-5-2009-B.patch

(Note the "-B")

It contains additional fixes for the Origin of the output image in file

itkGradientImageToBloxBoundaryPointImageFilter.txx
(0016502)
Luis Ibanez   
2009-05-16 12:06   
The patch itk-portable-round-2009-05-16.patch has just been committed.
(0018149)
Bradley Lowekamp   
2009-10-21 16:57   
The patch itk-templatedmathround-2009-07-30.patch has been committed with some modifications.
(0022566)
Hans Johnson   
2010-10-21 11:15   
If this is still an issue, please reopen the bug