ijk to xyz (voxel to mm) coordinate transformations

AFNI Message Board
1

Dear All,

I just started using Afni. My routine inputs a CT .nii volume and a Nx3 set of voxel (ijk) coordinates to Afni (extracted from the CT from an upstream external software). Afni then transforms them to RAI xyz (mm) coordinates (1st set), then a subroutine generates a 2nd RAI coordinates after some coregistrations inside Afni. I need to generate back voxel coordinates from this 2nd RAI set, since I need to move the set of point coordinates into another software package using a different coordinate system.

So I need to know how Afni transforms the ijk to xyz, and back (xyz to ijk). I tried the black-box version (3dUndump to construct a 3D vol from these points, then 3dmaskdump with -ijk option to get the voxel indices. However, there is a rounding of the voxel coordinates inside afni (fractional voxel coordinates get rounded with loss of precision after the forward and backward transforms. So, I hope to do these transforms the old way, by using the matrix from the CT header and apply some matrix operations. I would need to know:

  1. How one can retrieve such a 4x4 matrix from the header (3dinfo CT generates a .geom line but not sure of the order if trying to organize the line in a 4x4;
  2. How is afni internally switching from ijk to xyz, and back (please provide formulas). My CT vol is in RPI orientation, but I can switch to any other with 3dresample -orient if required by these formulas.

Thank you very much, Octavian

2

Hi, Octavian-

In the first paragraph, what parts of AFNI are you referring to? The GUI, or some other set of programs/functions?

(In the 2nd paragraph, I am not sure I would call 3dUndump+3dmaskdump a black box-- each specifies what it does.)

*notation note below: for simplicity, I sometimes abbreviate “(i, j, k)” as “ijk”, and similarly for xyz.

The matrix that maps a set of (i, j, k) values to an (x, y, z) location is internally 4x4, but the bottom row is constrained to be (0, 0, 0, 1), so it is effectly 3x4, where the first three columns for a 3x3 matrix that multiples a column vector of (i, j, k) values, and the last column is then added to the result.

If you want to see the 3x4 matrix that maps a set of (i, j, k) values to a physical location in scanner coordinates for a given DSET, you can run:


cat_matvec DSET::IJK_TO_DICOM_REAL

… (noting that “DSET” above is generic placeholder for a filename, but “IJK_TO_DICOM_REAL” is actually a specific label that happens to be written in all-capital letters) and if you want to see the 3x4 matrix that can map (i, j, k) to an associated “cardinal” set of coordinates, that the GUI uses, you can run:


cat_matvec DSET::IJK_TO_DICOM

You can add “-4x4” as an option to either to see the “full” 4x4 matrix. You could redirect the output of either to a file by putting “> FILE” or “>> FILE” at the end (overwriting or appending FILE, respectively).

You could invert one of these, if you would like:


cat_matvec DSET::IJK_TO_DICOM_REAL -I

… should map (x, y, z) coordinates to (i, j, k)

You can perform matrix multiplication of these, such as matrix multiplying the forward and inverse matrix of each, by naming more than one matrix in order, and including a “-I” after one, if you wish. Thus,


cat_matvec -4x4 DSET::IJK_TO_DICOM_REAL  DSET::IJK_TO_DICOM_REAL -I

… produces


            1             0             0             0
            0             1             0             0
            0             0             1             0
            0             0             0             1

As to your 2nd question, I don’t quite understand what you are asking. Are you saying you have a CT dset (DSET_CT) and another dset (DSET_OTHER) that are in the same space (so xyz correspondence is meaningful), and you want to specify an ijk value in DSET_CT and figure out the ijk value in DSET_OTHER that corresponds to the same xyz?

–pt

3

Thank you very much, this is very helpful, Octavian