3dNwarpApply or 3dAllineate with affine transforms for different grid sizes

AFNI Message Board
1

Hi AFNI experts,
I am running into an issue for applying some transforms. I have a T1 brain which I have rigidly aligned to my EPI scan (tx_1.1D) with the cmass flag in 3dAllineate. I have that same T1 brain also rigidly aligned to the MNI template for AC-PC alignment (tx_2.1D) using -cmass. I use cat_matvec to combine the transforms as follows:

cat_matvec tx_1.1D -I > tx_1_inv.1D
cat_matvec tx_1_inv.1D tx_2.1D > tx_EPItoACPC.1D

If I apply this transform using 3dAllineate to the T1 image (aligned to EPI after tx_1) it looks correct for APCP alignment with MNI template:
3dAllineate -base tpl.nii -source T1toEPI.nii -1Dmatrix_apply tx_EPItoACPC.1D

When I apply the same transform to the EPI scan it looks fairly shifted in the saggital plane by 15mm or so and also below the template by several millimeters.

  1. Does the matrix coming from 3dAllineate take into account the initial cmass operation?
  2. Could any of this be related to the fact that my T1 is 0.8mm iso while the EPI resolution is 2mm iso?

I also tried a few things with 3dNwarpApply

3dNwarpApply -nwarp ‘IDENT(EPI.nii) tx_EPItoACPC.1D’
3dNwarpApply -nwarp ‘IDENT(EPI.nii) tx_1_inv.1D IDENT(T1.brain.nii) tx_2.1D’

and a few other combinations but no success.

If I run it in two stages it works fine:

3dAllineate -1Dmatrix_apply tx_1_inv.1D -source EPI.nii -base T1.brain.nii -prefix temp1.nii
3dAllineate -1Dmatrix_apply tx_2.1D -source temp1.nii -base tpl.nii -prefix EPItotpl.nii

Any help for how to concatenate the warps properly would be appreciated.

Thanks,
Ajay

2

I think align_epi_anat.py includes all the transformations you want here. Use -epi2anat, -rigid_body and -tlrc_apar. The -tlrc_apar takes a dataset as input that has been affinely aligned to a standard space with @auto_tlrc. Similar operations can be done as aea_opts within afni_proc.py’s framework too.

3

Thanks for the suggestion Daniel I will try it out.

I had one general question regarding some of the outputs with 1Dmatrix_save and -1Dparam_save.

Specifically when looking at the rotation matrix (param 4-6) from 1Dmatrix_save and the x,y,z rotational components in 1Dparam_save I could not get them to match exactly. When I converted the param rotational components from degrees to radians it would get close to the 1Dmatrix_save rotational components which I assume is precision error. However, I thought when looking at some pages that both param and 1Dmatrix_save were in degrees which made me want to verify what is the correct conversion? Or is it that one is referenced to the center of the image and one to the origin in RAI format?

Thanks,
Ajay

4

The 1Dmatrix_save option saves an xyz to xyz (in RAI coordinates) matrix, not one that includes any angles at all directly. The last column is the translation vector, so that should match the translation parameters saved with the 1Dparam_save file.

5

Hi Daniel,
I have a follow up question of how the affine matrix is organized based on the documentation below:

DEFINITION OF AFFINE TRANSFORMATION PARAMETERS

The 3x3 spatial transformation matrix is calculated as [S][D][U],
where [S] is the shear matrix,
[D] is the scaling matrix, and
[U] is the rotation (proper orthogonal) matrix.
Thes matrices are specified in DICOM-ordered (x=-R+L,y=-A+P,z=-I+S)
coordinates as:

[U] = [Rotate_y(param#6)] [Rotate_x(param#5)] [Rotate_z(param #4)]
(angles are in degrees)

[D] = diag( param#7 , param#8 , param#9 )

    [    1        0     0 ]        [ 1 param#10 param#11 ]

[S] = [ param#10 1 0 ] OR [ 0 1 param#12 ]
[ param#11 param#12 1 ] [ 0 0 1 ]

The shift vector comprises parameters #1, #2, and #3.

  1. In looking at this documentation, if there is only a rigid registration is it correct to say that for the 3x3 matrix (excluding translational components), that the lower diagonal is shearing terms (not present in rigid registration), the diagonal is the scaling term (should be 1), and then upper triangle are rotational terms (going from RAI to RAI).

  2. If there is shearing/scaling terms and rotation, then there is an interaction such that the shearing/scaling and rotation cannot be separated accurately. However, it is still organized such that shearing terms are in the lower triangle and scaling in the diagonal, and rotation iin the upper triangle. This would suggest that most of the rotational components are in the upper triangle and shearing in lower triangle, while one still must account for interactions in order to be precise.

When looking at cat_matvec there is a polar decomposition option:

-P = Do a polar decomposition on the 3x3 matrix part
of the mfile. This would result in an orthogonal
matrix (rotation only, no scaling) Q that is closest,
in the Frobenius distance sense, to the input matrix A.
Note: if A = R * S * E, where R, S and E are the Rotation,
Scale, and shEar matrices, respctively, Q does not
necessarily equal R because of interaction; Each of R,
S and E affects most of the columns in matrix A.

If looking at the affine matrix from the description above is it accurate to say the lower triangle of the resulting 3x3 matrix after -P are shearing terms, the diagonal is scaling terms (set to 1) and the upper triangle is rotational components after doing polar decomposition, even if this does not exactly equal the input matrix? The part that confused me is that the last sentence says R, S and E affects most of the columns in A so I am not 100% sure if the SDU are organized as columns or upper/lower triangle and diagonal (I am assuming the latter).

Thanks,
Ajay

6

I’m not sure about all that. I would say for a translation only, yes, you would have just 1’s on the diagonal and the right column would be the translation amounts. Use the -P option with cat_matvec to extract the rigid equivalent (rotation and translation) from the affine matrix. There is a new program to compute affine transformations from parameters, 1dApar2mat, that might be useful for exploring the relationship between parameters and affine transformation matrices.

7

Hi Daniel,
Thanks for the suggestion!
a) It seems that when there is only rotation the upper and lower triangle (3x3 matrix) are symmetric aside from the sign flip.
1dApar2mat 0 0 0 -0.081944 -0.103117 -0.169312 1.000000 1 1.000 0 0 0

mat44 1dApar2mat 0 0 0 -0.081944 -0.103117 -0.169312 1.000000 1 1.000 0 0 0 :

  0.999995     -0.001425      0.002955       0.000000
  0.001430      0.999997     -0.001800       0.000000
 -0.002952      0.001804      0.999994       0.000000

Polar decomposition is the same:
0.999995 -0.00142484 0.00295479 0
0.00143016 0.999997 -0.00179989 0
-0.00295221 0.00180411 0.999994 0

b) If there is shearing only then this is only in the lower triangle:

1dApar2mat 0 0 0 0 0 0 1.000000 1 1.000000 -0.000783 0.02345 0.000108

mat44 1dApar2mat 0 0 0 0 0 0 1.000000 1 1.000000 -0.000783 0.02345 0.000108 :

  1.000000      0.000000      0.000000       0.000000
 -0.000783      1.000000      0.000000       0.000000
  0.023450      0.000108      1.000000       0.000000

Polar decomposition would be the following:
0.999931 0.000391157 -0.0117242 0
-0.00039179 1 -5.17014e-05 0
0.0117242 5.62912e-05 0.999931 0

c) If there is shearing and rotation you get the upper and lower triangles not matching:
1dApar2mat 0 0 0 -0.081944 -0.103117 -0.169312 1.000000 1 1.000000 -0.000783 0.02345 0.000108 > test1.txt

mat44 1dApar2mat 0 0 0 -0.081944 -0.103117 -0.169312 1.000000 1 1.000000 -0.000783 0.02345 0.000108 :

  0.999995     -0.001425      0.002955       0.000000
  0.000647      0.999999     -0.001802       0.000000
  0.020498      0.001879      1.000063       0.000000

Polar Decomposition would be the following:
cat_matvec test1.txt -P
0.999961 -0.00105473 -0.00877026 0
0.00103852 0.999998 -0.0018528 0
0.00877219 0.00184362 0.99996 0

So in general it seems that doing the polar decomposition and applying the 3x3 matrix would give the rotational components, whereas if you use the original full matrix then the upper triangle is the rotational component while the lower triangle is a mixture of shearing and rotation.

This definitely helps a lot, thanks for pointing that new tool out.
Thanks,
Ajay