Exporting surface sources as NIfTI volume

💬 Support & Discussions
1

Hi,

I am computing two separate set of sources, one in volume (beamformer) and another on surfaces (eLORETA). I’d like to export these as NIfTI volumes. Of course, for this, I understand that surface sources will need to be somehow extrapolated from vertices to voxels. I wrote a snippet of code using a kdtree spatial indexing with a nearest neighbor lookup that kinda works for that purpose, but I’d rather rely on MNE built-in functionalities as much as possible. I found that it is pretty easy to export volume sources with something like

vol_stc.save_as_volume("src_vol_A52.nii", vol_fwd["src"])

which I think is great. I could not find the equivalence for surface sources. It does appear that some code was written as part as WIP: Plot glass brain #4496 that would accomplish most of this work (except the saving)

    elif isinstance(stc, SourceEstimate):
        n_vertices = sum(len(v) for v in stc.vertices)
        offset = 0
        surf_to_mri = 0.
        for hi, hemi in enumerate(['lh', 'rh']):
            ribbon = nib.load(op.join(subjects_dir, subject, 'mri',
                                      '%s.ribbon.mgz' % hemi))
            xfm = ribbon.header.get_vox2ras_tkr()
            mri_data = ribbon.get_data()
            ijk = np.array(np.where(mri_data)).T
            xyz = apply_trans(xfm, ijk) / 1000.
            row_ind = np.where(mri_data.ravel())[0]
            data = np.ones(row_ind.size)
            rr = read_surface(op.join(subjects_dir, subject, 'surf',
                                      '%s.white' % hemi))[0]
            rr /= 1000.
            rr = rr[stc.vertices[hi]]
            col_ind = _compute_nearest(rr, xyz) + offset
            surf_to_mri = surf_to_mri + sparse.csr_matrix(
                (data, (row_ind, col_ind)), shape=(mri_data.size, n_vertices))
            offset += len(stc.vertices[hi])

        data = surf_to_mri.dot(stc.data[:, time_idx])
        data.shape = mri_data.shape
        xfm = _read_talairach(subject, subjects_dir)

        affine = np.dot(xfm, ribbon.header.get_vox2ras())
        img = nib.Nifti1Image(data, affine)

but that PR was abandoned. Is there such a functionality somewhere in MNE base code (saving surface sources as volumes)? I could not find it but MNE has many crannies and nooks so I might just have missed it!

2

AFAIK we don’t have a built-in way to convert surface STCs to volume STCs (nor to export surface STCs directly to NIfTI, which is your ultimate goal). @larsoner is that right?

3

I think the ribbon.mgz is the right way to do this. We should probably make a PR to add a BaseSourceEstimate.as_volume with this code in it. We could compare to the output of mri_surf2vol in tests. @christian-oreilly do you want to give it a shot? If not, feel free to open an issue link back here

4

Yes, after writing this post, I ended up modifying existing code (from the abandoned PR if I remember well) and doing this

from mne.transforms import apply_trans
from mne.surface import read_surface, _compute_nearest
from scipy import sparse
import xarray as xr 
from nibabel.funcs import concat_images

def save_surface_src_volume(stc, subject, subjects_dir, time_downsample_factor=20, 
                            save_file_name="data.nii.gz"):
    n_vertices = sum(len(v) for v in stc.vertices)
    offset = 0
    surf_to_mri = 0.
    for hi, hemi in enumerate(['lh', 'rh']):
        ribbon = nib.load(Path(subjects_dir / subject / 'mri' / f'{hemi}.ribbon.mgz'))
        xfm = ribbon.header.get_vox2ras_tkr()
        mri_data = np.asanyarray(ribbon.dataobj)
        ijk = np.array(np.where(mri_data)).T
        xyz = apply_trans(xfm, ijk) / 1000.
        row_ind = np.where(mri_data.ravel())[0]
        data = np.ones(row_ind.size)
        rr = read_surface(Path(subjects_dir / subject / 'surf'/ f'{hemi}.white'))[0]
        rr /= 1000.
        rr = rr[stc.vertices[hi]]
        col_ind = _compute_nearest(rr, xyz) + offset
        surf_to_mri = surf_to_mri + sparse.csr_matrix(
            (data, (row_ind, col_ind)), shape=(mri_data.size, n_vertices))
        offset += len(stc.vertices[hi])

    source_data = xr.DataArray(stc.data, dims=["sources", "time"]) \
                    .rolling(time=time_downsample_factor, center=True) \
                    .mean()[:, time_downsample_factor//2:-time_downsample_factor//2:time_downsample_factor]

    data_all_times = []
    for time_data in source_data.transpose("time", "sources").values:
        data = surf_to_mri.dot(time_data)
        data = data.reshape(mri_data.shape).astype("float32")
        data_all_times.append(nib.Nifti1Image(data, ribbon.affine))

    nib.save(concat_images(data_all_times), save_file_name)


save_surface_src_volume(pruned_stc, subject, subjects_dir, 
                        save_file_name=f"eLORETA_surface_{subject}.nii.gz")

That worked fine. When viewed in FreeSurfer:

image
image2194×1476 319 KB

The challenge here is that with 256 X 256 X 256 standard MRI space, the amount of data fast becomes unmanageable. We had an ERP going from -0.2 to 1.0 s, sampled at 1000 Hz. Full 1200 time frame was requiring around 100 GB of RAM (would be significantly smaller once compressed on the HDD, but it still need to be generated first in an uncompressed form). The volume sources conversion to MRI solves this with a coarser spatial resolution, but this does not work well with the ribbon because it is too thin to give good results at lower resolution. So, I ended up putting a temporal resolution downsampling factor. With default value I set (20), I obtained an ~55 MB nii.gz file showing the sources at 50 Hz.

I could used this code to implement a BaseSourceEstimate.as_volume and submit a PR for that. I cannot commit to preparing a test that compares it with the output of mri_surf2vol though. I can check it, but in my experience these things that seem simple have a way to end up taking much more time than initially expected and I have limited time I can put on this.

5

This is still helpful! If one of us needs to take over the PR to get it thoroughly tested, that is still fine. Many hands make for light work, as they say.

1 Like