[Insight-users] ITK Image Coordinates / Problem withPhysicalPoint to Index Conversion

[Andreas Keil]

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.