[Neuroimaging] Reading FreeSurfer m3z files?

[Tim Schäfer]

Dear list,

I have attached a Python/numpy function that works like the Matlab function 'mris_read_m3z' from the official FreeSurfer repo (https://github.com/freesurfer/freesurfer/blob/dev/matlab/mris_read_m3z.m) at the end of this email.

My latest version is always 'read_m3z_file' in https://github.com/dfsp-spirit/brainload/blob/master/src/brainload/freesurferdata.py

There is a unit test called 'test_read_m3z_file' in https://github.com/dfsp-spirit/brainload/blob/master/tests/brainload/test_freesurferdata.py. The values which are checked in the unit test have been retrieved by running the Matlab implementation linked above on my example m3z file (the file is also in the brainload repo, see the test).

The current implementation does not convert the result to double (like the Matlab implementation does), and the order of the return values is different. Most likely it's not the most beautiful implementation as I'm rather new to numpy, but it works. Feel free to use it/improve it/whatever (license is MIT).

Best,

Tim

> On February 25, 2019 at 11:41 AM Tim Schäfer <ts+ml at rcmd.org> wrote:
> 
> 
> True, I'll give it a try and report back.
> 
> > On February 25, 2019 at 11:26 AM Matthew Brett <matthew.brett at gmail.com> wrote:
> > 
> > 
> > Hi,
> > 
> > That bit of code is verrry short - I guess it's a very simple format.
> >  It should be pretty easy to implement.
> > 
> > Cheers,
> > 
> > Matthew
> > 
> > On Sun, Feb 24, 2019 at 8:09 PM Noam Peled <peled.noam at gmail.com> wrote:
> > >
> > > I don't know any python library that can read those files. I attached a Matlab file that reads m3z files. You can also find it here:
> > > https://www.mail-archive.com/freesurfer@nmr.mgh.harvard.edu/msg51239.html
> > >
> > >  If you convert it to Python, let me know :)
> > >
> > > Best,
> > > Noam
> > >
> > > On Fri, Feb 22, 2019 at 8:07 AM Tim Schäfer <ts+ml at rcmd.org> wrote:
> > >>
> > >> Dear list,
> > >>
> > >> does anybody know of a Python library which can read FreeSurfer m3z files, like <subject>/mri/transforms/talairach.m3z?
> > >>
> > >> These are binary files encoding a dense vector field for a 3D registration. The format seems to be a custom FreeSurfer format, but I'm not 100% sure on that. (The respective function in FreeSurfer is called GCAMwrite [1], and there is a Matlab implementation that is called mris_write_m3z [2], it seems.)
> > >>
> > >> Best,
> > >>
> > >> --
> > >> Tim
> > >>
> > >> [1] https://github.com/freesurfer/freesurfer/blob/dev/utils/gcamorph.c
> > >> [2] https://github.com/freesurfer/freesurfer/blob/dev/matlab/mris_write_m3z.m
> > >> _______________________________________________
> > >> Neuroimaging mailing list
> > >> Neuroimaging at python.org
> > >> https://mail.python.org/mailman/listinfo/neuroimaging
> > >
> > > _______________________________________________
> > > Neuroimaging mailing list
> > > Neuroimaging at python.org
> > > https://mail.python.org/mailman/listinfo/neuroimaging
> > _______________________________________________
> > Neuroimaging mailing list
> > Neuroimaging at python.org
> > https://mail.python.org/mailman/listinfo/neuroimaging
> 
> -- 
> Tim
> 0176 84273709
> _______________________________________________
> Neuroimaging mailing list
> Neuroimaging at python.org
> https://mail.python.org/mailman/listinfo/neuroimaging

-- 
Tim

# python function follows
def read_m3z_file(m3z_file):
    """
    Read a file in FreeSurfer m3z format.

    Read a file in FreeSurfer m3z format, usually mri/transforms/talairach.m3z of a subject. An m3z file is a gzipped binary file containing a dense vector field that describes a 3D registration between two volumes/images. This implementation follows the Matlab implementation from the FreeSurfer source repository at github, see https://github.com/freesurfer/freesurfer/blob/dev/matlab/mris_read_m3z.m.

    Parameters
    ----------
    m3z_file: str
        Path to a file in m3z format

    Returns
    -------
    vol_orig: numpy 4D array
        Array of np.single, containing 3 numbers per voxel. The shape is (width, height, depth, 3), where width, height, and depth are the dimensions of the source volume.

    vol_dest: numpy 4D array
        Array of np.single, containing 3 numbers per voxel. The shape is (width, height, depth, 3), where width, height, and depth are the dimensions of the destination volume.

    vol_ind0: numpy 4D array
        Array of np.int32, containing 3 numbers per voxel. The shape is (width, height, depth, 3).

    meta_data: dictionary
        Contains meta data on the registration. The following keys are included:
            - version: float, File format version
            - width: int, Volume width
            - height: int, Volume height
            - depth: int, Volume depth
            - spacing: int, voxel spacing of the volume
            - exp_k: float, exp_k of the volume

    Examples
    --------
    Read a talairach.m3z file that has been generated by FreeSurfer's ```recon-all``` command.

    >>> import brainload.freesurferdata as fsd
    >>> import os
    >>> m3z_file = os.path.join(os.getenv('HOME'), 'my_study_data', 'subject1', 'mri', 'transforms', 'talairach.m3z')
    >>> vol_orig, vol_dest, vol_ind0, meta_data = fsd.read_m3z_file(m3z_file)
    """
    m3z = gzip.open(m3z_file, 'rb')
    fdata = m3z.read()      # read the whole file contents
    # Read the file header
    data_start_pos = 24
    (version, width, height, depth, spacing, exp_k) = struct.unpack(">fiiiif", fdata[:data_start_pos])
    meta_data = {}
    meta_data['version'] = version
    meta_data['width'] = width
    meta_data['height'] = height
    meta_data['depth'] = depth
    meta_data['spacing'] = spacing
    meta_data['exp_k'] = exp_k
    meta_data['data_start_pos'] = data_start_pos

    numbers_per_voxel = 9
    bytes_uint8 = 4
    numbers_to_read = width * height * depth * numbers_per_voxel
    bytes_to_read = numbers_to_read * bytes_uint8
    data_end_pos = data_start_pos + bytes_to_read
    vol_data = struct.unpack(">%dB" % (numbers_to_read * 4), fdata[data_start_pos:data_end_pos])    # vol_data is a tuple
    meta_data['data_end_pos'] = data_end_pos
    meta_data['bytes_to_read'] = bytes_to_read
    meta_data['numbers_to_read'] = numbers_to_read
    meta_data['remaining_data_tag'] = struct.unpack(">i", fdata[data_end_pos:data_end_pos+4])[0]    # the first 4 bytes after the data are an integer tag, that holds information on what kind of data follows. We do not read that data but check the tag here.
    meta_data['file_end_pos'] = len(fdata)

    # create indices into the vol_data where data resides
    indices = np.kron(np.ones((12,1), dtype=np.int), range(width * height * depth))
    vox_offset = 9 * 4 * indices
    to_repeat = np.array([3, 2, 1, 0, 7, 6, 5, 4, 11, 10, 9, 8])
    reshaped_offsets = np.reshape(vox_offset, -1, order='F').copy()
    repeated = np.kron(np.ones((width * height * depth, ), dtype=int), to_repeat)
    inds = repeated + reshaped_offsets;
    inds = inds.astype(int)

    # Extract the three interleaved volumes and permute the result to match the original looped read. Start with viol_orig.
    buf = np.array(vol_data, dtype=np.uint8)
    buf_at_inds = buf[inds]
    single_casted = buf_at_inds.view(dtype=np.single)
    vol_orig_step1 = np.reshape(single_casted, (3, depth, width, height), order='F').copy()
    vol_orig = np.transpose(vol_orig_step1, (3, 2, 1, 0))

    # OK, vol_orig is reconstructed. Now for vol_dest
    inds = inds + 12
    buf_at_inds = buf[inds]
    single_casted = buf_at_inds.view(dtype=np.single)
    vol_dest_step1 = np.reshape(single_casted, (3, depth, width, height), order='F').copy()
    vol_dest = np.transpose(vol_dest_step1, (3, 2, 1, 0))

    # now for vol_ind0
    inds = inds + 12
    buf_at_inds = buf[inds]
    casted = buf_at_inds.view(dtype=np.int32)
    vol_ind0_step1 = np.reshape(casted, (3, depth, width, height), order='F').copy()
    vol_ind0 = np.transpose(vol_ind0_step1, (3, 2, 1, 0))
    return vol_orig, vol_dest, vol_ind0, meta_data

--
Dr. Tim Schäfer
Postdoc Computational Neuroimaging
Department of Child and Adolescent Psychiatry, Psychosomatics and Psychotherapy
University Hospital Frankfurt, Goethe University Frankfurt am Main, Germany