Coordinate space API

[darenlee]

Hi -

I'd like to start this thread to address coordinate space transform support in the 3D Volume/Atlas API. I've come up with an approach very similar to the one used in the Nifti file format so hopefully that means it is on the right track.

First some definitions:
1) Voxel Space (i,j,k): indices into data array
i = 0 .. dim[1]-1
j = 0 .. dim[2]-1
k = 0 .. dim[3]-1
2) Coordinate Space (x,y,z): coordinates relative to some origin defined by a standard (ie: Talairach, Waxholm, etc)

The basic idea: given a point in voxel space (i,j,k), the transform gives the point in the coordinate space (x,y,z).

Implementation: In the JLayerImage3D class, it holds the transform represented by a 4x4 affine transform matrix. Internally, the JLayerImage3D class will convert from 2D image coordinates to voxel coordinates. By default, if the affine transform exists, it will then convert from voxel coords to coordinate space coords by matrix multiplication. Theoritcally, I'm using the standard computer graphics transformation matrix (4x1 row vector x 4x4 matrix) but for our case, the 4th column is ignored

Caveats: If any physical dimensions need to be included in the transform, it must be included in the affine transformation matrix (I adopted this from the Nifti format; see below for more details).

API:

// set the affine transform matrix
public void setCoordinateSystemTransform ( float[16] );

/**
 * Converts from 3D voxel coordinates to 3D coordinate system space by 
 * multiplying voxel coordinates by the coordinate transformation matrix.
 * If no coordinate system transform matrix has been specified, it is 
 * assumed to be the identity matrix. Method allows for inplace computation.
 * 
 * @param vCoordinateCoords output 3D coordinate space coordinates
 * @param vVoxelCoords input 3D voxel coordinates
 * @return true if conversion was performed; false if coordinate transform matrix (m_pMatrixCoordinateSystem) is null
 */
public boolean convertVoxelCoordsToCoordinateCoords ( dluVector3i vCoordinateCoords, dluVector3i vVoxelCoords );

NOTE: convertVoxelCoordsToCoordinateCoords() can always be overriden by any plugins.

Let me know if you foresee any problems with this approach.

Thanks,
Daren



---------

Nifti header documentation:

http://nifti.nimh.nih.gov/pub/dist/s...tilib/nifti1.h

/*---------------------------------------------------------------------------*/
/* 3D IMAGE (VOLUME) ORIENTATION AND LOCATION IN SPACE:
   ---------------------------------------------------
   There are 3 different methods by which continuous coordinates can
   attached to voxels.  The discussion below emphasizes 3D volumes, and
   the continuous coordinates are referred to as (x,y,z).  The voxel
   index coordinates (i.e., the array indexes) are referred to as (i,j,k),
   with valid ranges:
     i = 0 .. dim[1]-1
     j = 0 .. dim[2]-1  (if dim[0] >= 2)
     k = 0 .. dim[3]-1  (if dim[0] >= 3)
   The (x,y,z) coordinates refer to the CENTER of a voxel.  In methods
   2 and 3, the (x,y,z) axes refer to a subject-based coordinate system,
   with
     +x = Right  +y = Anterior  +z = Superior.
   This is a right-handed coordinate system.  However, the exact direction
   these axes point with respect to the subject depends on qform_code
   (Method 2) and sform_code (Method 3).

   N.B.: The i index varies most rapidly, j index next, k index slowest.
    Thus, voxel (i,j,k) is stored starting at location
      (i + j*dim[1] + k*dim[1]*dim[2]) * (bitpix/8)
    into the dataset array.

   N.B.: The ANALYZE 7.5 coordinate system is
      +x = Left  +y = Anterior  +z = Superior
    which is a left-handed coordinate system.  This backwardness is
    too difficult to tolerate, so this NIFTI-1 standard specifies the
    coordinate order which is most common in functional neuroimaging.

   N.B.: The 3 methods below all give the locations of the voxel centers
    in the (x,y,z) coordinate system.  In many cases, programs will wish
    to display image data on some other grid.  In such a case, the program
    will need to convert its desired (x,y,z) values into (i,j,k) values
    in order to extract (or interpolate) the image data.  This operation
    would be done with the inverse transformation to those described below.

   N.B.: Method 2 uses a factor 'qfac' which is either -1 or 1; qfac is
    stored in the otherwise unused pixdim[0].  If pixdim[0]=0.0 (which
    should not occur), we take qfac=1.  Of course, pixdim[0] is only used
    when reading a NIFTI-1 header, not when reading an ANALYZE 7.5 header.

   N.B.: The units of (x,y,z) can be specified using the xyzt_units field.

   METHOD 1 (the "old" way, used only when qform_code = 0):
   -------------------------------------------------------
   The coordinate mapping from (i,j,k) to (x,y,z) is the ANALYZE
   7.5 way.  This is a simple scaling relationship:

     x = pixdim[1] * i
     y = pixdim[2] * j
     z = pixdim[3] * k

   No particular spatial orientation is attached to these (x,y,z)
   coordinates.  (NIFTI-1 does not have the ANALYZE 7.5 orient field,
   which is not general and is often not set properly.)  This method
   is not recommended, and is present mainly for compatibility with
   ANALYZE 7.5 files.

   METHOD 2 (used when qform_code > 0, which should be the "normal" case):
   ---------------------------------------------------------------------
   The (x,y,z) coordinates are given by the pixdim[] scales, a rotation
   matrix, and a shift.  This method is intended to represent
   "scanner-anatomical" coordinates, which are often embedded in the
   image header (e.g., DICOM fields (0020,0032), (0020,0037), (0028,0030),
   and (0018,0050)), and represent the nominal orientation and location of
   the data.  This method can also be used to represent "aligned"
   coordinates, which would typically result from some post-acquisition
   alignment of the volume to a standard orientation (e.g., the same
   subject on another day, or a rigid rotation to true anatomical
   orientation from the tilted position of the subject in the scanner).
   The formula for (x,y,z) in terms of header parameters and (i,j,k) is:

     [ x ]   [ R11 R12 R13 ] [        pixdim[1] * i ]   [ qoffset_x ]
     [ y ] = [ R21 R22 R23 ] [        pixdim[2] * j ] + [ qoffset_y ]
     [ z ]   [ R31 R32 R33 ] [ qfac * pixdim[3] * k ]   [ qoffset_z ]

   The qoffset_* shifts are in the NIFTI-1 header.  Note that the center
   of the (i,j,k)=(0,0,0) voxel (first value in the dataset array) is
   just (x,y,z)=(qoffset_x,qoffset_y,qoffset_z).

   The rotation matrix R is calculated from the quatern_* parameters.
   This calculation is described below.

   The scaling factor qfac is either 1 or -1.  The rotation matrix R
   defined by the quaternion parameters is "proper" (has determinant 1).
   This may not fit the needs of the data; for example, if the image
   grid is
     i increases from Left-to-Right
     j increases from Anterior-to-Posterior
     k increases from Inferior-to-Superior
   Then (i,j,k) is a left-handed triple.  In this example, if qfac=1,
   the R matrix would have to be

     [  1   0   0 ]
     [  0  -1   0 ]  which is "improper" (determinant = -1).
     [  0   0   1 ]

   If we set qfac=-1, then the R matrix would be

     [  1   0   0 ]
     [  0  -1   0 ]  which is proper.
     [  0   0  -1 ]

   This R matrix is represented by quaternion [a,b,c,d] = [0,1,0,0]
   (which encodes a 180 degree rotation about the x-axis).

   METHOD 3 (used when sform_code > 0):
   -----------------------------------
   The (x,y,z) coordinates are given by a general affine transformation
   of the (i,j,k) indexes:

     x = srow_x[0] * i + srow_x[1] * j + srow_x[2] * k + srow_x[3]
     y = srow_y[0] * i + srow_y[1] * j + srow_y[2] * k + srow_y[3]
     z = srow_z[0] * i + srow_z[1] * j + srow_z[2] * k + srow_z[3]

   The srow_* vectors are in the NIFTI_1 header.  Note that no use is
   made of pixdim[] in this method.

   WHY 3 METHODS?
   --------------
   Method 1 is provided only for backwards compatibility.  The intention
   is that Method 2 (qform_code > 0) represents the nominal voxel locations
   as reported by the scanner, or as rotated to some fiducial orientation and
   location.  Method 3, if present (sform_code > 0), is to be used to give
   the location of the voxels in some standard space.  The sform_code
   indicates which standard space is present.  Both methods 2 and 3 can be
   present, and be useful in different contexts (method 2 for displaying the
   data on its original grid; method 3 for displaying it on a standard grid).

   In this scheme, a dataset would originally be set up so that the
   Method 2 coordinates represent what the scanner reported.  Later,
   a registration to some standard space can be computed and inserted
   in the header.  Image display software can use either transform,
   depending on its purposes and needs.

   In Method 2, the origin of coordinates would generally be whatever
   the scanner origin is; for example, in MRI, (0,0,0) is the center
   of the gradient coil.

   In Method 3, the origin of coordinates would depend on the value
   of sform_code; for example, for the Talairach coordinate system,
   (0,0,0) corresponds to the Anterior Commissure.

[cgusafs]

Couple of things here:

The explicit result of the convert voxel to coordinate method is a boolean. What's the meaning of this result? True if the conversion succeeded? If so, what constitutes failure? Having an identity matrix as the transform? Having no (or null) matrix as the transform?

I also think being able to convert from coordinate space to voxel space is important. My understanding of what's going to happen is that the user will blissfully work in his desired coordinate space (wachy-home, or whatever), and the viewers will have to do the grunt work of converting back and forth (which is as it should be). The Neuroterrain server only understands, at least at this instant in time, it's own atlas space, which is a scaled and translated version of what would be voxel space. So, to fulfil the user's desire for a view using a cutting plane specified in wachy-home space, the viewer will have to do the translation using something like convertCoordinateCoordsToVoxelCoords() which would use the inverse of the matrix used in convertVoxelCoordsToCoordinateCoords(), which is what the viewers would have to use to convert mouse clicks on it's graphics turf to wachy-home space.

(And at no time, during that last paragraph, did my hands leave the end of my arms).

[darenlee]

Originally Posted by cgusafs

Couple of things here:

The explicit result of the convert voxel to coordinate method is a boolean. What's the meaning of this result? True if the conversion succeeded? If so, what constitutes failure? Having an identity matrix as the transform? Having no (or null) matrix as the transform?

Sorry for the incomplete documentation. I've edited the above post with the return values. This preliminary code is committed to SVN so you also take a peek there in JLayerImage3D.

I also think being able to convert from coordinate space to voxel space is important.

Ah very good point. I'll have to think about this a little bit more and will get back to you.

[cgusafs]

Well, it's not really a big deal. Either require the client code to provide the inverse transform, which is reasonable, since the client code needs to come up with the forward transform, or sportingly compute the inverse matrix directly (there's gotta be a bunch of liberaries to do this kind of stuff), then just add a second conversion method. Both forward and reverse methods could just call a private method to do the actual work, providing the correct matrix...

[darenlee]

Hi Carl -

Thanks for your input. I was thinking along the same lines as you and have opted for your first option - client code must specify the voxel->coord and coord->voxel transforms.

I've changed the API a bit:

1. added the transform matrix as an argument
2. allow coordinate space values to be floating point
3. voxel coordinates are always integers

This latest version is committed to SVN.

/**
 * Converts from 3D voxel coordinates to 3D coordinate system space by 
 * multiplying voxel coordinates by the coordinate transformation matrix.
 * If the coordinate system transformation matrix is not specified or not 
 * the right dimensions, it is assumed to be the identity matrix.<br>
 * <br>
 * @see dlGraphics.util.dlUtils#dluTransformVector3f(dluVector3f, dluVector3f, float[])  
 * @param vCoordinateCoords output 3D coordinate space coordinates (floating point)
 * @param vVoxelCoords input 3D voxel coordinates (integer)
 * @param pMatrix 4x4 coordinate transformation matrix
 * @return true if conversion was performed; false if coordinate transformation matrix is null or not 4x4
 */
public boolean convertVoxelToCoordinateCoords ( dluVector3f vCoordinateCoords, dluVector3i vVoxelCoords, float[] pMatrix )

/**
 * Converts from 3D coordinate system space to 3D voxel coordinates by 
 * multiplying voxel coordinates by the coordinate transformation matrix.
 * If the coordinate system transform matrix is not specified or not 
 * the right dimensions, it is assumed to be the identity matrix. The output
 * voxel coordinates are rounded to the nearest integer values.<br>
 * <br>
 * @see dlGraphics.util.dlUtils#dluTransformVector3f(dluVector3f, dluVector3f, float[])
 * 
 * @param vVoxelCoords output 3D voxel coordinates (nearest integer)
 * @param vCoordinateCoords input 3D coordinate space coordinates (floating point)
 * @param pMatrix 4x4 coordinate transformation matrix
 * @return true if conversion was performed; false if coordinate transformation matrix is null or not 4x4
 */
public boolean convertCoordinateToVoxelCoords ( dluVector3i vVoxelCoords, dluVector3f vCoordinateCoords, float[] pMatrix )

[cgusafs]

OK, so it looks like you now have two conversion methods - one takes a vector representation of a coordinate point in integer format as input and outputs a corresponding point using floats, and one takes a vector of floats as input, and outputs integers.

Is there a reason for requiring voxel coordinates to be integers? There are several reasons why floating point representation of voxel coordinates are of interest to me:

The default coordinate system used by the Neuroterrain server is a normalized system with x values running between -1.0 and 1.0. So, when the client deals with the server, it has to speak floating point. I'd rather eliminate floating point-integer conversions until the last possible moment, if any conversion needs to be done at all.

The Neuroterrain server is capable of slicing at arbitrary angles, so to generate voxel coordinates here, one needs to do a fair bit of math to convert the two dimension slice coordinates to three dimensions; the result of this computation is also going to be in floating point (or doubles).

We have been working on the alignment/registration/whatever-is-the-right-term of the MBL data sets to our Neuroterrain atlas. Given those parameters, we can specify points on an MBL section in Neuroterrain internal coordinates (and thus converted by the plugin to voxel space coordinates), but by their nature, these are going to be floating point values, since it turns out that the sections are not at uniform orientation. I could round or truncate, but prefer to do it at the last possible moment.

Finally, having input and output coordinates both in floating point means one less method to deal with, especially if the transformation matrix is to be supplied on each invocation.

What sayest thou?

[darenlee]

The conversion methods above are actually wrappers to the more generic transform function in dlUtils:

public static boolean dluTransformVector3f ( dluVector3f vCoordsOut, dluVector3f vCoordsIn, float[] pMatrix )

Hopefully this is what you are looking for so you can call this method directly.

[cgusafs]

Yes, I can do that, that's no problem; on reflection, I probaly will.

Just for reference, I've got a TransformMatrix class which I originally dragged in from the Samurai Microtome. For this discussion, it's cardinal advantage is that it allows you to apply the transform to a point, and returns the transformed point as an explicit result - doesn't do any checking to make sure your matrix is reasonable or not, a legacy of my C heritage, I guess* - the user is supposed to be paying attention. It also allows you to build additional matrices by multiplying one with another, so for example, you can do a scale change based on zooming or conversion from server space to client space and multiply it into the voxel->waxholm conversion.

I don't know if it's worthwhile to adapt this class to handle arbitrary and/or affine transforms, thereby remaining self-contained, vs. the benefit of using a generally available scheme.

--
*Java - a battery operated circular saw with a blade guard, heavy work gloves, and goggles. Java users go through life without disasters, but it takes a while to make a cut.
C++ - a circular saw with a 10 amp motor, no blade guard, no gloves, no goggles. Experienced C++ users can't provide full sets of fingerprints, and tend to squint. But, they can do more than straight cuts, and boy, do they cut quick!

This useless analogy brought to you by the Band-Aids company.