test_freesurfer.py

 #!/usr/bin/env python
#
# test_freesurfer.py -
#
# Author: Paul McCarthy <pauldmccarthy@gmail.com>
#

import os
import os.path as op
import shutil

import numpy as np

import nibabel             as nib
import nibabel.freesurfer  as nibfs

import pytest

import fsl.data.freesurfer as fslfs

from .test_mesh import (CUBE_VERTICES, CUBE_TRIANGLES_CCW)

from . import tempdir, touch


def gen_freesurfer_geometry(fname, verts, tris):
    nibfs.write_geometry(fname, verts, tris)

def gen_freesurfer_morph(fname, vdata):
    nibfs.write_morph_data(fname, vdata)

def gen_freesurfer_label(fname, verts, vdata):
    with open(fname, 'wt') as f:
        f.write('# Comment\n')
        f.write('{}\n'.format(len(verts)))
        for v, vd in zip(verts, vdata):
            f.write('{} 0.0 0.0 0.0 {}\n'.format(v, vd))


def gen_freesurfer_annot(fname, labels, rgbal, names):
    nibfs.write_annot(fname, labels, rgbal, names)

def test_FreesurferMesh_create():

    verts = np.array(CUBE_VERTICES)
    tris  = np.array(CUBE_TRIANGLES_CCW)

    with tempdir():
        gen_freesurfer_geometry('lh.pial', verts, tris)

        mesh = fslfs.FreesurferMesh('lh.pial')

        assert mesh.name == 'lh.pial'
        assert np.all(np.isclose(mesh.vertices, verts))
        assert np.all(np.isclose(mesh.indices,  tris))


def test_FreesurferMesh_create_loadall():

    verts = np.array(CUBE_VERTICES)
    tris  = np.array(CUBE_TRIANGLES_CCW)

    with tempdir():
        gen_freesurfer_geometry('lh.pial', verts, tris)

        vertSets = ['lh.orig', 'lh.white', 'lh.inflated']
        for vs in vertSets:
            shutil.copy('lh.pial', vs)

        mesh = fslfs.FreesurferMesh('lh.pial', loadAll=True)

        assert list(sorted(mesh.vertexSets())) == \
            list(sorted([op.abspath(p) for p in ['lh.pial'] + vertSets]))


def test_loadVertices():

    verts = np.array(CUBE_VERTICES)
    tris  = np.array(CUBE_TRIANGLES_CCW)

    with tempdir():
        gen_freesurfer_geometry('lh.pial', verts,     tris)
        gen_freesurfer_geometry('rh.pial', verts * 2, tris)

        # bad
        gen_freesurfer_geometry('lh.orig', verts[:-1, :], tris)

        np.savetxt('verts.txt', verts * 3)

        mesh = fslfs.FreesurferMesh('lh.pial')

        assert np.all(np.isclose(mesh.loadVertices('rh.pial'),   verts * 2))
        assert np.all(np.isclose(mesh.loadVertices('verts.txt'), verts * 3))

        with pytest.raises(ValueError):
            mesh.loadVertices('lh.orig')


def test_loadVertexData_morph():

    verts = np.array(CUBE_VERTICES)
    tris  = np.array(CUBE_TRIANGLES_CCW)

    with tempdir():
        vdata1 = np.random.randint(1, 100, verts.shape[0])
        vdata2 = np.random.randint(1, 100, verts.shape[0])
        vdata3 = np.random.randint(1, 100, verts.shape[0])

        gen_freesurfer_geometry('lh.pial', verts, tris)

        gen_freesurfer_morph('lh.curv',      vdata1)
        gen_freesurfer_morph('lh.thickness', vdata2)
        np.savetxt(          'vdata.txt',    vdata3)

        # bad
        gen_freesurfer_morph('rh.sulc', vdata2[:-1])

        mesh = fslfs.FreesurferMesh('lh.pial')

        assert np.all(np.isclose(mesh.loadVertexData('lh.curv'),      vdata1.reshape(-1, 1)))
        assert np.all(np.isclose(mesh.loadVertexData('lh.thickness'), vdata2.reshape(-1, 1)))
        assert np.all(np.isclose(mesh.loadVertexData('vdata.txt'),    vdata3.reshape(-1, 1)))

        with pytest.raises(ValueError):
            mesh.loadVertexData('rh.sulc')


def test_loadVertexData_label():

    verts  = np.array(CUBE_VERTICES)
    tris   = np.array(CUBE_TRIANGLES_CCW)
    nverts = verts.shape[0]

    with tempdir():

        # Currently, vertex scalar data is
        # ignored by the FreesurferMesh class
        lverts = np.random.choice(np.arange(nverts), 4)
        vdata  = np.random.randint(1, 100, 4)

        gen_freesurfer_geometry('lh.pial', verts, tris)
        gen_freesurfer_label('lh.vdata.label', lverts, vdata)

        mesh = fslfs.FreesurferMesh('lh.pial')

        exp = np.zeros((nverts, 1))
        exp[lverts, :] = 1

        assert np.all(np.isclose(mesh.loadVertexData('lh.vdata.label'), exp))


def test_loadVertexData_mgh():

    verts  = np.array(CUBE_VERTICES)
    tris   = np.array(CUBE_TRIANGLES_CCW)
    nverts = verts.shape[0]

    with tempdir():

        data = np.random.randint(1, 100, (nverts, 1, 1), dtype=np.int32)
        img  = nib.freesurfer.mghformat.MGHImage(data, np.eye(4))

        nib.save(img, 'lh.vdata.mgh')

        gen_freesurfer_geometry('lh.pial', verts, tris)

        mesh = fslfs.FreesurferMesh('lh.pial')
        assert np.all(np.isclose(mesh.loadVertexData('lh.vdata.mgh'), data.reshape(-1, 1)))

        del img
        del mesh
        img = None


def test_loadVertexData_annot():

    import nibabel.info as nibinfo

    # assume nibabel 2.*
    # nibabel 2.2.1 is broken w.r.t. .annot files.
    if nibinfo._version_minor == 2 and nibinfo._version_micro <= 1:
        return

    verts = np.array(CUBE_VERTICES)
    tris  = np.array(CUBE_TRIANGLES_CCW)

    with tempdir():
        nlabels      = 3
        names        = ['label {}'.format(l) for l in range(1, nlabels + 1)]
        rgba         = np.random.randint(0, 255, (nlabels, 4), dtype=np.int32)
        labels       = list(range(nlabels)) + list(np.random.randint(0, nlabels, verts.shape[0] - nlabels))
        labels       = np.array(labels, dtype=np.int32)

        np.random.shuffle(labels)

        gen_freesurfer_geometry('lh.pial', verts, tris)
        gen_freesurfer_annot('lh.aparc.annot', labels, rgba, names)

        mesh = fslfs.FreesurferMesh('lh.pial')

        vdfile = op.abspath('lh.aparc.annot')

        loaded = mesh.loadVertexData(vdfile)
        ergbal, enames = mesh.getVertexDataColourTable(vdfile)

        enames = [n.decode() for n in enames]

        assert np.all(np.isclose(loaded, labels.reshape(-1, 1)))
        assert list(enames) == list(names)
        assert np.all(np.isclose(ergbal[:, :4], rgba))


def test_relatedGeometryFiles():

    with tempdir():
        touch('not.geometry')
        touch('lh.pial')
        touch('lh.sphere')
        touch('lhinflated')
        touch('rh.blob')
        touch('rhsphere')
        touch('rh.pial')
        touch('rh.orig')
        touch('rh.white')
        touch('rh.blob')

        def a(paths):
            return [op.abspath(p) for p in paths]

        assert        fslfs.relatedGeometryFiles('not.geometry') == []
        assert        fslfs.relatedGeometryFiles('lh.blob')      == []
        assert        fslfs.relatedGeometryFiles('lh.pial')      == a(['lh.sphere'])
        assert        fslfs.relatedGeometryFiles('lh.sphere')    == a(['lh.pial'])
        assert        fslfs.relatedGeometryFiles('rh.blob')      == []
        assert sorted(fslfs.relatedGeometryFiles('rh.pial'))     == a(['rh.orig', 'rh.white'])
        assert sorted(fslfs.relatedGeometryFiles('rh.orig'))     == a(['rh.pial', 'rh.white'])
        assert sorted(fslfs.relatedGeometryFiles('rh.white'))    == a(['rh.orig', 'rh.pial'])


def test_relatedVertexDataFiles():
    with tempdir():
        touch('lh.pial')
        touch('rh.pial')
        touch('rh.area')
        touch('foo.mgz')
        touch('lh.blob.mgz')
        touch('lh.blob.mgh')
        touch('lhthickness')
        touch('lh.thickness')
        touch('lh.curv')
        touch('lh.area')
        touch('lh.sulc')
        touch('lh.what.label')
        touch('lh.aparc.annot')

        def a(paths):
            return [op.abspath(p) for p in paths]

        assert fslfs.relatedVertexDataFiles('rh.pial') == a(['rh.area'])
        assert fslfs.relatedVertexDataFiles('rh.area') == []

        exp = sorted(['lh.blob.mgz',
                      'lh.blob.mgh',
                      'lh.thickness',
                      'lh.curv',
                      'lh.area',
                      'lh.sulc',
                      'lh.what.label',
                      'lh.aparc.annot'])

        assert sorted(fslfs.relatedVertexDataFiles('lh.pial')) == a(exp)


def test_findReferenceImage():
    with tempdir():

        surf = op.join('surf', 'lh.pial')
        t1   = op.join('mri',  'T1.mgz')

        os.mkdir('surf')
        os.mkdir('mri')
        touch(surf)
        touch(t1)

        assert fslfs.findReferenceImage(surf) == op.abspath(t1)

        os.remove(t1)

        assert fslfs.findReferenceImage(surf) is None