From fdea6f0b34e6a93c4bebc716ff9967e29f127c04 Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Wed, 24 Jan 2018 18:19:52 +0000
Subject: [PATCH] Unit tests for freesurfer module
---
tests/__init__.py | 5 +
tests/test_freesurfer.py | 285 +++++++++++++++++++++++++++++++++++++++
2 files changed, 290 insertions(+)
create mode 100644 tests/test_freesurfer.py
diff --git a/tests/__init__.py b/tests/__init__.py
index a837c9c28..f44dfea33 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -46,6 +46,11 @@ def tempdir():
shutil.rmtree(testdir)
+def touch(fname):
+ with open(fname, 'wt') as f:
+ pass
+
+
class CaptureStdout(object):
"""Context manager which captures stdout and stderr. """
diff --git a/tests/test_freesurfer.py b/tests/test_freesurfer.py
new file mode 100644
index 000000000..451406205
--- /dev/null
+++ b/tests/test_freesurfer.py
@@ -0,0 +1,285 @@
+ #!/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')
+ nib.save(img, '/Users/paulmc/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)))
+
+
+
+
+
+
+
+def test_loadVertexData_annot():
+
+ # nibabel 2.2.1 is broken w.r.t. .annot files.
+ # 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)
+
+ 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
--
GitLab