[Insight-users] Re: SoftwareGuide.pdf is ready.

[Luis Ibanez]

Hi Simon,

Please seee comments below,

    Luis

----------------------------------------
Simon Warfield wrote:
> 
>   I read over it rather quickly trying to understand the relationship
> between a point in space and a voxel in an image. Perhaps you could clarify
> section 4.1.4 regarding origin and spacing  (page 34)
> ``These coordinates correspond to the position of the first pixel of 
 >
the image with respect to an arbitrary reference system in physical 
space...user's
> responsibility [to keep it consistent]'' 
>      (there is a typo in this section by the way)
>
Thanks for pointing this out. (I'll fix the typo)


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

Probably the best way to analyze the relationship between
voxels and space points in ITK is with the diagrams of the
ResampleImageFilter. This is section 5.7.1: pages 105-118.

The diagram in Figure 5.40 is a good start (page 116 in text)
(page 123 of the PDF file).

First of all, The space coordinates associated with one pixel are those
of the pixel center. The figure shows two image grids and for both of
them the centers of pixels are represented with circles. Two filled
circles illustrate the origin of the image. That is, the coordinates
of the first pixel in the image. The rectangles represent the area
of coverage of the pixel. (Note that by 'pixel' we mean N-D pixels,
so, 'voxels' are just 3D pixels).

The image to the right, for example has "origin" = (50,130). This means
that the center of the first pixel is at 50mm along x and 130mm along y
from the origin of space coordinates.

The actual origin of space coordinates is completly arbitrary...
There is noting we can do about it since this comes from the very
acquisition process when every scanner is calibrated and some arbitraty
origin of coordinates is selected on the table.

Image spacing is defined as the separation between the centers of
neighbor pixels. For example in the left image grid, spacing is
20mm in X and 30mm in Y.

The Index<> class in ITK provide the integer coordinates associated
with the image grid. So in the right image, the blue arrow is comming
out of the pixel with index = [1,2] (since the origin pixel is [0,0]).

Registration in ITK is performed in space coordinates, not pixels.
This allows to take into account the origin and spacing of each image.

A detailed discussion on the mapping of the pixel (1,2) is presented
in this section.

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

> 
>   Just where in the voxel does the registration and interpolation framework
> understand the coordinate to point to ?  the center of the voxel or one of
> the corners of the voxel ?
>   Or is this up to the user to choose ?
>   

The center of the voxel is the point.

This is consistent with the physics of the acquisition process.
The value of a pixel is the representation of a continuous function
that has been digitized. A nearest neighbor interpolator will then
assign the same value to all the points surrounding the space point
associated with the coordinates of the pixel (the pixel center).

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

> Chapter 5 doesn't seem to address it.
> 


Thanks for pointing this out.
It is definetly worth to expand on this on the section describing
image interpolators in the registration framework...
(added to the TODO list...)


--------------------------------------------------------------------
> 7.7.1 isn't completely clear.
> ``NearestNeigbhorInterpolateImageFunction simply use the intensity of
>         ^^^^^^^^ typo by the way
> the nearest grid position. That is, it assumes that the image intensity
> is piecewise constant with jumps mid-way between grid positions. This
> interpolation scheme is cheap as it does not require any floating point
> computations.''
> 
>   Does this imply that voxels extend + and - 0.5 times the spacing from 
> each grid position ? That is, each grid position is intended to be at the
> center of a voxel, and similarly the origin is intended to point to the
> center of a voxel.
> 

That's right. The space position of a pixel in ITK is the
point associated with the center of the pixel.

We may have to add a figure here to make it clearer.
The final meaning should be that given a pixel at space
coordinates (x,y), The nearest neighbor interpolator should
assign this same value to all the points having coordinates
x',y' such that

    x' is in [ x - Xspacing/2 , x + Xspacing/2 ]
    y' is in [ y - Yspacing/2 , y + Yspacing/2 ]


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

>   How does this compare to the VTK scheme ?  How are origins and coordinates
> mapped between ITK and VTK in the various ITK*VTK filters ?
> 
> 

Well,
As far as the scheme for interpolation goes, they are both
the same.

The difference is what is called a 'Pixel' in ITK and
what is called a 'Pixel' in VTK.

Let's talk about "data points" first.  If you have a file with a 2D
image made of 5 x 7 data points, VTK will load it in a 5x7 grid and
will say that your image is made of 4x6 pixels. Since VTK consider
the pixels to be the squares whose corners are defined by data points.

So, if you consider the grid of data points, the VTK 'pixels' are the
Delaunay regions of the grid.  ITK on the other hand will load the
file in a 5x7 grid and tell you that you have a 5x7 pixels image.
The ITK 'pixels' are the Voronoi regions of the grid defined by
data points.

Despite the superficial difference, the good news is that when it
comes to interpolation both points of view are conciliated since
for example, linear interpolation is done in the Delaunay region
using the 4 data points in the corners, both in ITK an VTK.

Probably the most remarkable difference in the interpolation approach
is the treatment of borders.  VTK will not let you go outside the
borders of the data point grid. ITK will allow you to go 1/2 spacing
outside the border.

A VTK image has physical extent = (Nx-1)*Xspacing, (Ny-1)*Yspacing

An ITK image has physical extent = Nx * Xspacing , Ny * Yspacing

Where

    Nx = number of data points along x
    Ny = number of data points along y


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

A diagram will explain this much  better....


Please let us know if you have any other
questions or comments


Thanks for your feedback


     Luis