Day 4: Transforming coordinates

Yesterday we isolated the right amygdala using nibabel and freesurfer.
The coordinates were in voxels. mne.setup_volume_source_space wants
the coordinates in meters. Luckily, nibabel provides information about
coordinate transformations. Expanding our code from yesterday ...

In [1]:

import nibabel as nib
import numpy as np
import mne
from mne.datasets import spm_face

# get the segment information
subjects_dir = spm_face.data_path()
aseg_fname = subjects_dir + '/subjects/spm/mri/aseg.mgz'

# load segment info using nibabel
aseg = nib.load(aseg_fname)
aseg_data = aseg.get_data()
ix = aseg_data==54 # index for the right amygdala

# get the indices in x, y, z space
iix = []
for i in range(ix.shape[0]):
    for j in range(ix.shape[1]):
        for k in range(ix.shape[2]):
            if ix[i, j, k]:
                iix.append([i, j, k])
iix = np.array(iix) # convert to array

# get the header information
aseg_hdr = aseg.get_header()

# get the transformation matrix
trans = aseg_hdr.get_vox2ras_tkr()

# convert using the transformation matrix
xyz = np.dot(iix, m[:3,:3].T)+m[:3,3]

# convert to meters
xyz /= 1000.

I should qualify this by saying that lots of trial and error
went into finding the correct way to transpose and order
these matrices but the end result checks out in freeview.

So now we have one half of the pos dictionary. Tomorrow we'll
go over orientations.