diff --git a/fsl/data/image.py b/fsl/data/image.py index 8a2465d2971710b396b2b947e82661fa7979c02d..678db9f1932ac1661f0952626e5c0030df20ef75 100644 --- a/fsl/data/image.py +++ b/fsl/data/image.py @@ -460,7 +460,7 @@ class Nifti(notifier.Notifier): @deprecation.deprecated(deprecated_in='1.1.0', - removed_in='1.2.0', + removed_in='2.0.0', details='Use ndims instead') def is4DImage(self): """Returns ``True`` if this image is 4D, ``False`` otherwise. """ @@ -554,7 +554,16 @@ class Nifti(notifier.Notifier): @memoize.Instanceify(memoize.memoize) + @deprecation.deprecated(deprecated_in='1.2.0', + removed_in='2.0.0', + details='Use voxToScaledVoxMat instead') def voxelsToScaledVoxels(self): + """See :meth:`voxToScaledVoxMat`.""" + return self.voxToScaledVoxMat + + + @property + def voxToScaledVoxMat(self): """Returns a transformation matrix which transforms from voxel coordinates into scaled voxel coordinates, with a left-right flip if the image appears to be stored in neurological order. @@ -563,6 +572,12 @@ class Nifti(notifier.Notifier): _format_of_the_matrix_used_by_FLIRT.2C_and_how_does_it_relate_to\ _the_transformation_parameters.3F """ + return self.__voxToScaledVoxMat() + + + @memoize.Instanceify(memoize.memoize) + def __voxToScaledVoxMat(self): + """See :meth:`voxToScaledVoxMat`. """ shape = list(self.shape[ :3]) pixdim = list(self.pixdim[:3]) @@ -576,6 +591,20 @@ class Nifti(notifier.Notifier): return voxToPixdimMat + @property + def scaledVoxToVoxMat(self): + """Returns a transformation matrix which transforms from scaled voxels + into voxels, the inverse of the :meth:`voxToScaledVoxMat` transform. + """ + return self.__scaledVoxToVoxMat() + + + @memoize.Instanceify(memoize.memoize) + def __scaledVoxToVoxMat(self): + """See :meth:`scaledVoxToVoxMat`. """ + return transform.invert(self.voxToScaledVoxMat) + + def sameSpace(self, other): """Returns ``True`` if the ``other`` image (assumed to be a :class:`Nifti` instance) has the same dimensions and is in the diff --git a/fsl/utils/transform.py b/fsl/utils/transform.py index 4fe7fe8e4cfd33e97c4e4cd72db90081c5ce8582..6ff409b272a54962f8f09742c4ef8046a4c4aea3 100644 --- a/fsl/utils/transform.py +++ b/fsl/utils/transform.py @@ -555,10 +555,10 @@ def flirtMatrixToSform(flirtMat, srcImage, refImage): :arg refImage: Reference :class:`.Image` """ - srcScaledVoxelMat = srcImage.voxelsToScaledVoxels() - refScaledVoxelMat = refImage.voxelsToScaledVoxels() + srcScaledVoxelMat = srcImage.voxToScaledVoxMat + refInvScaledVoxelMat = refImage.voxToScaledVoxMat + refInvScaledVoxelMat = refImage.scaledVoxToVoxMat refVoxToWorldMat = refImage.voxToWorldMat - refInvScaledVoxelMat = invert(refScaledVoxelMat) return concat(refVoxToWorldMat, refInvScaledVoxelMat, @@ -581,15 +581,15 @@ def sformToFlirtMatrix(srcImage, refImage, srcXform=None): :attr:`.Nifti.voxToWorldMat` """ - srcScaledVoxelsToVoxelsMat = invert(srcImage.voxelsToScaledVoxels()) - srcVoxToWorldMat = srcImage.voxToWorldMat - refWorldToVoxMat = invert(refImage.voxToWorldMat) - refVoxelsToScaledVoxelsMat = refImage.voxelsToScaledVoxels() + srcScaledVoxToVoxMat = srcImage.scaledVoxToVoxMat + srcVoxToWorldMat = srcImage.voxToWorldMat + refWorldToVoxMat = refImage.worldToVoxMat + refVoxToScaledVoxMat = refImage.voxToScaledVoxMat if srcXform is not None: srcVoxToWorldMat = srcXform - return concat(refVoxelsToScaledVoxelsMat, + return concat(refVoxToScaledVoxMat, refWorldToVoxMat, srcVoxToWorldMat, - srcScaledVoxelsToVoxelsMat) + srcScaledVoxToVoxMat) diff --git a/fsl/version.py b/fsl/version.py index 52bbf39c4a29bfe012b6abd64005339ae5a028c5..b8c0a1da236d9f5f7d9b818c405f79024691149e 100644 --- a/fsl/version.py +++ b/fsl/version.py @@ -41,7 +41,7 @@ import re import string -__version__ = '1.1.1.dev' +__version__ = '1.2.0.dev' """Current version number, as a string. """ diff --git a/tests/test_image.py b/tests/test_image.py index 2bab26993425730b49681cd1728784d2f7e41ec6..5dc7aee98e6c1fb7cb6a30acd2b6b0125c81640f 100644 --- a/tests/test_image.py +++ b/tests/test_image.py @@ -12,7 +12,6 @@ import os.path as op import itertools as it import tempfile import shutil -import glob import pytest @@ -856,9 +855,12 @@ def _test_Image_5D(imgtype): assert img.ndims == 5 -def test_Image_voxelsToScaledVoxels(): +def test_Image_voxToScaledVox_analyze(): _test_Image_voxToScaledVox(0) +def test_Image_voxToScaledVox_nifti1(): _test_Image_voxToScaledVox(1) +def test_Image_voxToScaledVox_nifti2(): _test_Image_voxToScaledVox(2) + +def _test_Image_voxToScaledVox(imgtype): - imgTypes = [0, 1, 2] dims = [(10, 10, 10)] pixdims = [(-1, 1, 1), ( 1, 1, 1), @@ -879,14 +881,15 @@ def test_Image_voxelsToScaledVoxels(): return xf - for imgType, dim, pixdim in it.product(imgTypes, dims, pixdims): - nimg = make_image(imgtype=imgType, dims=dim, pixdims=pixdim) + for dim, pixdim in it.product(dims, pixdims): + nimg = make_image(imgtype=imgtype, dims=dim, pixdims=pixdim) img = fslimage.Image(nimg) - expected = expect(imgType, dim, pixdim) - result = img.voxelsToScaledVoxels() + expected = expect(imgtype, dim, pixdim) + invexpected = npla.inv(expected) - assert np.all(np.isclose(result, expected)) + assert np.all(np.isclose(expected, img.voxToScaledVoxMat)) + assert np.all(np.isclose(invexpected, img.scaledVoxToVoxMat)) def test_Image_sameSpace(): diff --git a/tests/test_transform.py b/tests/test_transform.py index 609df61df5aa842718b6f1a5690c3e1f321b1129..0ed50c64e36d3f181201d5f5003a4ec1282436be 100644 --- a/tests/test_transform.py +++ b/tests/test_transform.py @@ -92,6 +92,67 @@ def test_concat(): assert np.all(np.isclose(result, expected)) +def test_veclength(seed): + + def l(v): + v = np.array(v, copy=False).reshape((-1, 3)) + x = v[:, 0] + y = v[:, 1] + z = v[:, 2] + l = x * x + y * y + z * z + return np.sqrt(l) + + vectors = -100 + 200 * np.random.random((200, 3)) + + for v in vectors: + + vtype = random.choice((list, tuple, np.array)) + v = vtype(v) + + assert np.isclose(transform.veclength(v), l(v)) + + # Multiple vectors in parallel + result = transform.veclength(vectors) + expected = l(vectors) + assert np.all(np.isclose(result, expected)) + + +def test_normalise(seed): + + vectors = -100 + 200 * np.random.random((200, 3)) + + def parallel(v1, v2): + v1 = v1 / transform.veclength(v1) + v2 = v2 / transform.veclength(v2) + + return np.isclose(np.dot(v1, v2), 1) + + for v in vectors: + + vtype = random.choice((list, tuple, np.array)) + v = vtype(v) + vn = transform.normalise(v) + vl = transform.veclength(vn) + + assert np.isclose(vl, 1.0) + assert parallel(v, vn) + + # normalise should also be able + # to do multiple vectors at once + results = transform.normalise(vectors) + lengths = transform.veclength(results) + pars = np.zeros(200) + for i in range(200): + + v = vectors[i] + r = results[i] + + pars[i] = parallel(v, r) + + assert np.all(np.isclose(lengths, 1)) + assert np.all(pars) + + def test_scaleOffsetXform(): # Test numerically @@ -155,9 +216,6 @@ def test_scaleOffsetXform(): assert np.all(np.isclose(result, expected)) - - - def test_compose_and_decompose(): testfile = op.join(datadir, 'test_transform_test_compose.txt') @@ -187,14 +245,70 @@ def test_compose_and_decompose(): rots = [np.pi / 5, np.pi / 4, np.pi / 3] rmat = transform.axisAnglesToRotMat(*rots) xform = transform.compose([1, 1, 1], [0, 0, 0], rmat) - sc, of, rot = transform.decompose(xform) - sc = np.array(sc) - of = np.array(of) - rot = np.array(rot) - assert np.all(np.isclose(sc, [1, 1, 1])) - assert np.all(np.isclose(of, [0, 0, 0])) - assert np.all(np.isclose(rot, rots)) + # And the angles flag should cause decompose + # to return the rotation matrix, instead of + # the axis angls + sc, of, rot = transform.decompose(xform) + scat, ofat, rotat = transform.decompose(xform, angles=True) + scaf, ofaf, rotaf = transform.decompose(xform, angles=False) + + sc, of, rot = np.array(sc), np.array(of), np.array(rot) + scat, ofat, rotat = np.array(scat), np.array(ofat), np.array(rotat) + scaf, ofaf, rotaf = np.array(scaf), np.array(ofaf), np.array(rotaf) + + assert np.all(np.isclose(sc, [1, 1, 1])) + assert np.all(np.isclose(of, [0, 0, 0])) + assert np.all(np.isclose(scat, [1, 1, 1])) + assert np.all(np.isclose(ofat, [0, 0, 0])) + assert np.all(np.isclose(scaf, [1, 1, 1])) + assert np.all(np.isclose(ofaf, [0, 0, 0])) + + assert np.all(np.isclose(rot, rots)) + assert np.all(np.isclose(rotat, rots)) + assert np.all(np.isclose(rotaf, rmat)) + + +def test_rotMatToAxisAngles(seed): + + pi = np.pi + pi2 = pi / 2 + + for i in range(100): + + rots = [-pi + 2 * pi * np.random.random(), + -pi2 + 2 * pi2 * np.random.random(), + -pi + 2 * pi * np.random.random()] + + rmat = transform.axisAnglesToRotMat(*rots) + gotrots = transform.rotMatToAxisAngles(rmat) + + assert np.all(np.isclose(rots, gotrots)) + + +def test_rotMatToAffine(seed): + + pi = np.pi + pi2 = pi / 2 + + for i in range(100): + + rots = [-pi + 2 * pi * np.random.random(), + -pi2 + 2 * pi2 * np.random.random(), + -pi + 2 * pi * np.random.random()] + + if np.random.random() < 0.5: origin = None + else: origin = np.random.random(3) + + rmat = transform.axisAnglesToRotMat(*rots) + mataff = transform.rotMatToAffine(rmat, origin) + rotaff = transform.rotMatToAffine(rots, origin) + + exp = np.eye(4) + exp[:3, :3] = rmat + exp[:3, 3] = origin + + assert np.all(np.isclose(mataff, rotaff)) def test_axisBounds(): @@ -348,6 +462,33 @@ def test_transform_vector(seed): assert np.all(np.isclose(ptExpected, ptResult)) +def test_transformNormal(seed): + + normals = -100 + 200 * np.random.random((50, 3)) + + def tn(n, xform): + + xform = npla.inv(xform[:3, :3]).T + return np.dot(xform, n) + + for n in normals: + + scales = -10 + np.random.random(3) * 10 + offsets = -100 + np.random.random(3) * 200 + rotations = -np.pi + np.random.random(3) * 2 * np.pi + origin = -100 + np.random.random(3) * 200 + + xform = transform.compose(scales, + offsets, + rotations, + origin) + + expected = tn(n, xform) + result = transform.transformNormal(n, xform) + + assert np.all(np.isclose(expected, result)) + + def test_flirtMatrixToSform(): testfile = op.join(datadir, 'test_transform_test_flirtMatrixToSform.txt') @@ -396,92 +537,3 @@ def test_sformToFlirtMatrix(): assert np.all(np.isclose(result1, expected)) assert np.all(np.isclose(result2, expected)) - - - -def test_normalise(seed): - - vectors = -100 + 200 * np.random.random((200, 3)) - - def parallel(v1, v2): - v1 = v1 / transform.veclength(v1) - v2 = v2 / transform.veclength(v2) - - return np.isclose(np.dot(v1, v2), 1) - - for v in vectors: - - vtype = random.choice((list, tuple, np.array)) - v = vtype(v) - vn = transform.normalise(v) - vl = transform.veclength(vn) - - assert np.isclose(vl, 1.0) - assert parallel(v, vn) - - # normalise should also be able - # to do multiple vectors at once - results = transform.normalise(vectors) - lengths = transform.veclength(results) - pars = np.zeros(200) - for i in range(200): - - v = vectors[i] - r = results[i] - - pars[i] = parallel(v, r) - - assert np.all(np.isclose(lengths, 1)) - assert np.all(pars) - - -def test_veclength(seed): - - def l(v): - v = np.array(v, copy=False).reshape((-1, 3)) - x = v[:, 0] - y = v[:, 1] - z = v[:, 2] - l = x * x + y * y + z * z - return np.sqrt(l) - - vectors = -100 + 200 * np.random.random((200, 3)) - - for v in vectors: - - vtype = random.choice((list, tuple, np.array)) - v = vtype(v) - - assert np.isclose(transform.veclength(v), l(v)) - - # Multiple vectors in parallel - result = transform.veclength(vectors) - expected = l(vectors) - assert np.all(np.isclose(result, expected)) - - -def test_transformNormal(seed): - - normals = -100 + 200 * np.random.random((50, 3)) - - def tn(n, xform): - - xform = npla.inv(xform[:3, :3]).T - return np.dot(xform, n) - - for n in normals: - - scales = -10 + np.random.random(3) * 10 - offsets = -100 + np.random.random(3) * 200 - rotations = -np.pi + np.random.random(3) * 2 * np.pi - origin = -100 + np.random.random(3) * 200 - - xform = transform.compose(scales, - offsets, - rotations, - origin) - - expected = tn(n, xform) - result = transform.transformNormal(n, xform) - - assert np.all(np.isclose(expected, result))