prealigning CT and MRI: incorrect transformation read from freesurfer .lta file

[mmagnuski]

I've outlined the problem in more detail on the discourse, but because I think it might be either a bug or somthing not taken into account in the instructions in the tutorial, I post it here too.

I have to follow the "Note" section in the "Locating intracranial electrode contacts" tutorial, because automatic CT-MRI alignment is not satisfactory in a few patients. I successfully align the CT and MRI by hand in freeview and save the .lta transformation file. However when I read the .lta and transform it from vox->vox to ras->ras, as outlined in the Note section of the tutorial, I see that the CT and MRI are not correctly aligned (this is before running the registration again with starting_affine):

as a result running the registration with starting_affine fails, with the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[122], line 1
----> 1 reg_affine, _ = mne.transforms.compute_volume_registration(
      2     CT_orig, T1, pipeline=['rigid'],
      3     starting_affine=manual_reg_affine)

File <decorator-gen-27>:12, in compute_volume_registration(moving, static, pipeline, zooms, niter, starting_affine, verbose)

File ~\miniconda3\envs\mne\lib\site-packages\mne\transforms.py:1666, in compute_volume_registration(moving, static, pipeline, zooms, niter, starting_affine, verbose)
   1628 @verbose
   1629 def compute_volume_registration(moving, static, pipeline='all', zooms=None,
   1630                                 niter=None, *, starting_affine=None,
   1631                                 verbose=None):
   1632     """Align two volumes using an affine and, optionally, SDR.
   1633 
   1634     Parameters
   (...)
   1664     .. versionadded:: 0.24
   1665     """
-> 1666     return _compute_volume_registration(
   1667         moving, static, pipeline, zooms, niter,
   1668         starting_affine=starting_affine)[:2]

File ~\miniconda3\envs\mne\lib\site-packages\mne\transforms.py:1743, in _compute_volume_registration(moving, static, pipeline, zooms, niter, starting_affine)
   1740 if step in ('translation', 'rigid'):
   1741     dist = np.linalg.norm(reg_affine[:3, 3])
   1742     angle = np.rad2deg(_angle_between_quats(
-> 1743         np.zeros(3), rot_to_quat(reg_affine[:3, :3])))
   1744     logger.info(f'    Translation: {dist:6.1f} mm')
   1745     if step == 'rigid':

File ~\miniconda3\envs\mne\lib\site-packages\mne\transforms.py:1274, in rot_to_quat(rot)
   1256 """Convert a set of rotations to quaternions.
   1257 
   1258 Parameters
   (...)
   1271 quat_to_rot
   1272 """
   1273 rot = rot.reshape(rot.shape[:-2] + (9,))
-> 1274 return np.apply_along_axis(_one_rot_to_quat, -1, rot)

File <__array_function__ internals>:180, in apply_along_axis(*args, **kwargs)

File ~\miniconda3\envs\mne\lib\site-packages\numpy\lib\shape_base.py:379, in apply_along_axis(func1d, axis, arr, *args, **kwargs)
    375 except StopIteration as e:
    376     raise ValueError(
    377         'Cannot apply_along_axis when any iteration dimensions are 0'
    378     ) from None
--> 379 res = asanyarray(func1d(inarr_view[ind0], *args, **kwargs))
    381 # build a buffer for storing evaluations of func1d.
    382 # remove the requested axis, and add the new ones on the end.
    383 # laid out so that each write is contiguous.
    384 # for a tuple index inds, buff[inds] = func1d(inarr_view[inds])
    385 buff = zeros(inarr_view.shape[:-1] + res.shape, res.dtype)

File ~\miniconda3\envs\mne\lib\site-packages\mne\transforms.py:1226, in _one_rot_to_quat()
   1224 det = np.linalg.det(np.reshape(rot, (3, 3)))
   1225 if np.abs(det - 1.) > 1e-3:
-> 1226     raise ValueError('Matrix is not a pure rotation, got determinant != 1')
   1227 t = 1. + rot[0] + rot[4] + rot[8]
   1228 if t > np.finfo(rot.dtype).eps:

ValueError: Matrix is not a pure rotation, got determinant != 1

I'll post a link to the files and a reproducible code example below.

[mmagnuski]

The data are here.

The code to reproduce the problem:

import os
import os.path as op

import numpy as np
import matplotlib.pyplot as plt

import mne
import nibabel as nib

def plot_overlay(image, compare, title, thresh=None):
    """Define a helper function for comparing plots."""
    image = nib.orientations.apply_orientation(
        np.asarray(image.dataobj), nib.orientations.axcodes2ornt(
            nib.orientations.aff2axcodes(image.affine))).astype(np.float32)
    compare = nib.orientations.apply_orientation(
        np.asarray(compare.dataobj), nib.orientations.axcodes2ornt(
            nib.orientations.aff2axcodes(compare.affine))).astype(np.float32)
    if thresh is not None:
        compare[compare < np.quantile(compare, thresh)] = np.nan
    fig, axes = plt.subplots(1, 3, figsize=(12, 4))
    fig.suptitle(title)
    for i, ax in enumerate(axes):
        ax.imshow(np.take(image, [image.shape[i] // 2], axis=i).squeeze().T,
                  cmap='gray')
        ax.imshow(np.take(compare, [compare.shape[i] // 2],
                          axis=i).squeeze().T, cmap='gist_heat', alpha=0.5)
        ax.invert_yaxis()
        ax.axis('off')
    fig.tight_layout()

data_dir = r'CHANGE TO THE CORRECT PATH ON YOUR COMPUTER'
ct_fname = 'sub-U06_ct_Tilt_1.nii'
mri_fname = 'T1.mgz'

T1 = nib.load(op.join(data_dir, mri_fname))
CT_orig = nib.load(op.join(data_dir, ct_fname))

# load transform
manual_reg_affine_vox = mne.read_lta(op.join(
    data_dir, 'sub-U06_ct_aligned_manual.nii.lta'))

# convert from vox->vox to ras->ras
manual_reg_affine = \
    CT_orig.affine @ np.linalg.inv(manual_reg_affine_vox) \
    @ np.linalg.inv(CT_orig.affine)

# show alignment
CT_aligned = mne.transforms.apply_volume_registration(
    CT_orig, T1, manual_reg_affine, cval='1%')
plot_overlay(T1, CT_aligned, 'Aligned CT Overlaid on T1', thresh=0.95)

But the alignment looks correct in freeview:

freeview T1.mgz sub-U06_ct_Tilt_1.nii:colormap=heat:opacity=0.6:reg=sub-U06_ct_aligned_manual.nii.lta

[larsoner]

@alexrockhill can you look?

[alexrockhill]

Hmm the lta looks normal but from the error I would infer that there might be some scaling in the lta that's not rigid. Is that possible? Otherwise, it's just for logging, so you could try commenting those lines out and check if the result works. I looked in my own code and realized I was using:

manual_reg_affine_vox = mne.read_lta(
    CT_fname.replace('.mgz', '_manual.mgz.lta'))
# convert from vox->vox to ras->ras
manual_reg_affine = \
    CT_orig.affine @ np.linalg.inv(manual_reg_affine_vox) \
    @ np.linalg.inv(CT_orig.affine)
CT_aligned_fix_img, reg_affine_fix = affine_registration(
    moving=np.array(CT_orig.dataobj), static=np.array(T1.dataobj),
    moving_affine=CT_orig.affine, static_affine=T1.affine,
    pipeline=['rigid'], starting_affine=manual_reg_affine,
    level_iters=[100], sigmas=[0], factors=[1])

for this step. It would be ideal if it were all encapsulated within MNE but this is a workaround until we figure it out. I can make a PR if we can determine if the lta is purely rigid.

[alexrockhill]

I also looks like you skipped this step from the tutorial:

reg_affine, _ = mne.transforms.compute_volume_registration(
    CT_orig, T1, pipeline=['rigid'],
    starting_affine=manual_reg_affine)

which re-aligns based on the manual affine for machine-precision alignment. Was that on purpose? Doing this will give basically the same lta and I have tested that it runs for the sample data. That would fix any non-rigid alignment if it exists or if it's a rounding error or something.

[mmagnuski]

@alexrockhill

It also looks like you skipped this step from the tutorial: (...)

I did not skip the starting_affine step, this is exactly the step that raises the error I pasted in the first post. And it most likely raises the error because the manual transformation is not correct (ie it is different in mne to what I see in freeview). For more screenshots take a look at the discourse post.

Hmm the lta looks normal but from the error I would infer that there might be some scaling in the lta that's not rigid. Is that possible?

I don't think so, the only steps I performed in freeview when pre-aligning were translations and rotations. The only potentially problematic part is that the nii file of the ct has some tilt correction applied. This was done when creating the nii from dicoms (I think I did that with dcm2niix, but I'll have to check) - could it be present in the CT_orig.affine?

[mmagnuski]

I have two other subjects that require manual pre-alignment of the CT, but that did not have and corrections applied to the CT (at least looking at the postfixes that dcm2niix normally creates when it applies these corrections). I can check if the lta for these files is read and displayed correctly.

[alexrockhill]

Can you try the dev version of mne? I fixed a read_lta bug recently

[mmagnuski]

oh, sure, I will!

[mmagnuski]

Ok, great, the dev mne version of read_lta does the trick! Sorry for the noise!