mghimage.py

#!/usr/bin/env python
#
# mghimage.py - The MGHImage class
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#
"""This module provides the :class:`MGHImage` class, which can be used to load
Freesurfer ``mgh``/``mgz`` image files.
"""


import os.path as op

import            six
import numpy   as np
import nibabel as nib

import fsl.utils.path       as fslpath
import fsl.transform.affine as affine
import fsl.data.image       as fslimage


ALLOWED_EXTENSIONS = ['.mgz', '.mgh']
"""List of file extensions interpreted as MGH image files.
"""


EXTENSION_DESCRIPTIONS = ['Compressed MGH image', 'MGH image']
"""A description for each of the :attr:`ALLOWED_EXTENSIONS`."""


class MGHImage(fslimage.Image):
    """The ``MGHImage`` class is a NIFTI :class:`Image` which has been converted
    from a Freesurfer ``.mgh`` file.

    See:
     - https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/MghFormat
     - http://nipy.org/nibabel/reference/nibabel.freesurfer.html
    """

    def __init__(self, image, **kwargs):
        """Create a ``MGHImage``.

        :arg image: Name of MGH file, or a
                    ``nibabel.freesurfer.mghformat.MGHImage`` instance.

        All other arguments are passed through to :meth:`Image.__init__`
        """

        if isinstance(image, six.string_types):
            filename = op.abspath(image)
            name     = op.basename(filename)
            image    = nib.load(image)
        else:
            name     = 'MGH image'
            filename = None

        data     = np.asanyarray(image.dataobj)
        xform    = image.affine
        pixdim   = image.header.get_zooms()
        vox2surf = image.header.get_vox2ras_tkr()

        # the image may have an affine which
        # transforms the data into some space
        # with a scaling that is different to
        # the pixdims. So we create a header
        # object with both the affine and the
        # pixdims, so they are both preserved.
        #
        # Note that we have to set the zooms
        # after the s/qform, otherwise nibabel
        # will clobber them with zooms gleaned
        # fron the affine.
        header = nib.nifti1.Nifti1Header()
        header.set_data_shape(data.shape)
        header.set_sform(xform)
        header.set_qform(xform)
        header.set_zooms(pixdim)

        fslimage.Image.__init__(self,
                                data,
                                header=header,
                                name=name,
                                dataSource=filename,
                                **kwargs)

        if filename is not None:
            self.setMeta('mghImageFile', filename)

        self.__voxToSurfMat   = vox2surf
        self.__surfToVoxMat   = affine.invert(vox2surf)
        self.__surfToWorldMat = affine.concat(xform, self.__surfToVoxMat)
        self.__worldToSurfMat = affine.invert(self.__surfToWorldMat)


    def save(self, filename=None):
        """Overrides :meth:`.Image.save`.  If a ``filename`` is not provided,
        converts the original (MGH) file name into a NIFTI filename, before
        passing it to the :meth:`.Image.save` method.
        """
        if filename is None:
            filename = self.dataSource

        filename = fslpath.removeExt(filename, ALLOWED_EXTENSIONS)

        return fslimage.Image.save(self, filename)


    @property
    def mghImageFile(self):
        """If this ``MGHImage`` was loaded from a file, returns the file
        name. Otherwise returns ``None``.
        """
        return self.getMeta('mghImageFile', None)


    @property
    def voxToSurfMat(self):
        """Returns an affine which can be used to transform voxel
        coordinates into the surface coordinate system for this image.

        See: http://www.grahamwideman.com/gw/brain/fs/coords/fscoords.htm
        See: https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems
        """
        return self.__voxToSurfMat


    @property
    def surfToVoxMat(self):
        """Returns an affine which can be used to transform surface
        coordinates into the voxel coordinate system for this image.
        """
        return self.__surfToVoxMat


    @property
    def surfToWorldMat(self):
        """Returns an affine which can be used to transform surface
        coordinates into the world coordinate system for this image.
        """
        return self.__surfToWorldMat


    @property
    def worldToSurfMat(self):
        """Returns an affine which can be used to transform world
        coordinates into the surface coordinate system for this image.
        """
        return self.__worldToSurfMat


def voxToSurfMat(img):
    """Generate an affine which can transform the voxel coordinates of
    the given image into a corresponding Freesurfer surface coordinate
    system (known as "Torig", or "vox2ras-tkr").

    See https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems

    :arg img: An :class:`.Image` object.

    :return:  A ``(4, 4)`` matrix encoding an affine transformation from the
              image voxel coordinate system to the corresponding Freesurfer
              surface coordinate system.
    """

    zooms = np.array(img.pixdim[:3])
    dims  = img.shape[ :3] * zooms / 2

    xform        = np.zeros((4, 4), dtype=np.float32)
    xform[ 0, 0] = -zooms[0]
    xform[ 1, 2] =  zooms[2]
    xform[ 2, 1] = -zooms[1]
    xform[ 3, 3] = 1
    xform[:3, 3] = [dims[0], -dims[2], dims[1]]

    return xform