diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 2630bdf90fb8990aa2cb351c377482cbc11923ae..09d15c5fca9ef6fac5b8d8dc67705e816e8c8ed1 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,6 +1,25 @@
This document contains the ``fslpy`` release history in reverse chronological
order.
+3.3.1 (Thursday 8th October 2020)
+---------------------------------
+
+
+Changed
+^^^^^^^
+
+
+* The :func:`.affine.decompose` and :func:`.affine.compose` functions now
+ have the ability to return/accept shear components.
+
+
+Fixed
+^^^^^
+
+
+* Fixed a bug in the :func:`.affine.decompose` function which was corrupting
+ the scale estimates when given an affine containing shears.
+
3.3.0 (Tuesday 22nd September 2020)
-----------------------------------
diff --git a/fsl/transform/affine.py b/fsl/transform/affine.py
index c7fb71f2a0571f8a1318a7e7b84749cfa498a9f2..47bdd4b371d0d1a9a26b19177795c0b4290b74c3 100644
--- a/fsl/transform/affine.py
+++ b/fsl/transform/affine.py
@@ -120,7 +120,7 @@ def scaleOffsetXform(scales, offsets):
return xform
-def compose(scales, offsets, rotations, origin=None):
+def compose(scales, offsets, rotations, origin=None, shears=None):
"""Compose a transformation matrix out of the given scales, offsets
and axis rotations.
@@ -133,6 +133,8 @@ def compose(scales, offsets, rotations, origin=None):
:arg origin: Origin of rotation - must be scaled by the ``scales``.
If not provided, the rotation origin is ``(0, 0, 0)``.
+
+ :arg shears: Sequence of three shear values
"""
preRotate = np.eye(4)
@@ -154,6 +156,7 @@ def compose(scales, offsets, rotations, origin=None):
scale = np.eye(4, dtype=np.float64)
offset = np.eye(4, dtype=np.float64)
rotate = np.eye(4, dtype=np.float64)
+ shear = np.eye(4, dtype=np.float64)
scale[ 0, 0] = scales[ 0]
scale[ 1, 1] = scales[ 1]
@@ -164,10 +167,15 @@ def compose(scales, offsets, rotations, origin=None):
rotate[:3, :3] = rotations
- return concat(offset, postRotate, rotate, preRotate, scale)
+ if shears is not None:
+ shear[0, 1] = shears[0]
+ shear[0, 2] = shears[1]
+ shear[1, 2] = shears[2]
+
+ return concat(offset, postRotate, rotate, preRotate, scale, shear)
-def decompose(xform, angles=True):
+def decompose(xform, angles=True, shears=False):
"""Decomposes the given transformation matrix into separate offsets,
scales, and rotations, according to the algorithm described in:
@@ -175,8 +183,7 @@ def decompose(xform, angles=True):
320-323 in *Graphics Gems II*, James Arvo (editor), Academic Press, 1991,
ISBN: 0120644819.
- It is assumed that the given transform has no perspective components. Any
- shears in the affine are discarded.
+ It is assumed that the given transform has no perspective components.
:arg xform: A ``(3, 3)`` or ``(4, 4)`` affine transformation matrix.
@@ -184,6 +191,8 @@ def decompose(xform, angles=True):
as axis-angles, in radians. Otherwise, the rotation matrix
is returned.
+ :arg shears: Defaults to ``False``. If ``True``, shears are returned.
+
:returns: The following:
- A sequence of three scales
@@ -191,6 +200,7 @@ def decompose(xform, angles=True):
was a ``(3, 3)`` matrix)
- A sequence of three rotations, in radians. Or, if
``angles is False``, a rotation matrix.
+ - If ``shears is True``, a sequence of three shears.
"""
# The inline comments in the code below are taken verbatim from
@@ -199,7 +209,7 @@ def decompose(xform, angles=True):
# The next step is to extract the translations. This is trivial;
# we find t_x = M_{4,1}, t_y = M_{4,2}, and t_z = M_{4,3}. At this
# point we are left with a 3*3 matrix M' = M_{1..3,1..3}.
- xform = xform.T
+ xform = np.array(xform).T
if xform.shape == (4, 4):
translations = xform[ 3, :3]
@@ -214,6 +224,7 @@ def decompose(xform, angles=True):
# The process of finding the scaling factors and shear parameters
# is interleaved. First, find s_x = |M'_1|.
sx = np.sqrt(np.dot(M1, M1))
+ M1 = M1 / sx
# Then, compute an initial value for the xy shear factor,
# s_xy = M'_1 * M'_2. (this is too large by the y scaling factor).
@@ -230,7 +241,7 @@ def decompose(xform, angles=True):
# The second row is normalized, and s_xy is divided by s_y to
# get its final value.
M2 = M2 / sy
- sxy = sxy / sy
+ sxy = sxy / sx
# The xz and yz shear factors are computed as in the preceding,
sxz = np.dot(M1, M3)
@@ -245,8 +256,8 @@ def decompose(xform, angles=True):
# the third row is normalized, and the xz and yz shear factors are
# rescaled.
M3 = M3 / sz
- sxz = sxz / sz
- syz = sxz / sz
+ sxz = sxz / sx
+ syz = syz / sy
# The resulting matrix now is a pure rotation matrix, except that it
# might still include a scale factor of -1. If the determinant of the
@@ -266,7 +277,13 @@ def decompose(xform, angles=True):
if angles: rotations = rotMatToAxisAngles(R.T)
else: rotations = R.T
- return np.array([sx, sy, sz]), translations, rotations
+ retval = [np.array([sx, sy, sz]), translations, rotations]
+
+ if shears:
+ retval.append(np.array((sxy, sxz, syz)))
+
+ return tuple(retval)
+
def rotMatToAffine(rotmat, origin=None):
diff --git a/tests/test_transform/test_affine.py b/tests/test_transform/test_affine.py
index 79a17fba6f9a3895b124b06f279be1c9b9f1b0f9..469c6c13d98950e379a9a042467b0a4f4d517dcb 100644
--- a/tests/test_transform/test_affine.py
+++ b/tests/test_transform/test_affine.py
@@ -226,8 +226,10 @@ def test_compose_and_decompose():
xform = lines[i * 4: i * 4 + 4]
xform = np.genfromtxt(xform)
- scales, offsets, rotations = affine.decompose(xform)
- result = affine.compose(scales, offsets, rotations)
+ scales, offsets, rotations, shears = affine.decompose(
+ xform, shears=True)
+
+ result = affine.compose(scales, offsets, rotations, shears=shears)
assert np.all(np.isclose(xform, result, atol=1e-5))
@@ -235,8 +237,10 @@ def test_compose_and_decompose():
# different rotation origin, but we test
# explicitly passing the origin for
# completeness
- scales, offsets, rotations = affine.decompose(xform)
- result = affine.compose(scales, offsets, rotations, [0, 0, 0])
+ scales, offsets, rotations, shears = affine.decompose(
+ xform, shears=True)
+ result = affine.compose(
+ scales, offsets, rotations, origin=[0, 0, 0], shears=shears)
assert np.all(np.isclose(xform, result, atol=1e-5))
diff --git a/tests/test_transform/testdata/test_transform_test_compose.txt b/tests/test_transform/testdata/test_transform_test_compose.txt
index e61f31bdd1578e526d9c9c78142b0dd31b5a2381..497fba8606f8d4617a9c8998b741eff2ce2504ef 100644
--- a/tests/test_transform/testdata/test_transform_test_compose.txt
+++ b/tests/test_transform/testdata/test_transform_test_compose.txt
@@ -23,7 +23,7 @@
1.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0
0.0 0.0 -1.0 0.0
-0.0 0.0 0.0 1.0
+0.0 0.0 0.0 1.0
0.5 0.0 0.0 0.0
0.0 2.4 0.0 0.0
@@ -35,22 +35,22 @@
0.0 0.0 1.6 53.0
0.0 0.0 0.0 1.0
-0.5 0.0 0.0 10.0
+0.5 0.0 0.0 10.0
0.0 0.0 -1.6 -6.7877502441
0.0 2.4 0.0 29.7931518555
-0.0 0.0 0.0 1.0
+0.0 0.0 0.0 1.0
0.0 -2.4 -0.0 30.6160888672
0.0 0.0 -1.6 -6.7877502441
0.5 0.0 0.0 50.5961608887
-0.0 0.0 0.0 1.0
+0.0 0.0 0.0 1.0
-0.0 1.3765834472 1.3106432709 0.0736983719
0.0 -1.9659649063 0.9177222982 14.1062102814
0.5 0.0 0.0 50.5961608887
0.0 0.0 0.0 1.0
- 0.453766452 -0.7004337463 0.4832135801 16.789591468
+ 0.453766452 -0.7004337463 0.4832135801 16.789591468
-0.0159784155 1.6056689296 1.1880787388 -16.2943532298
-0.2093817023 -1.640493784 0.9565424959 66.1123321137
0.0 0.0 0.0 1.0
@@ -109,3 +109,8 @@
0.041017 -0.081580 -6.319922 72.090378
-0.000000 0.381914 -1.365036 -41.159451
0.000000 0.000000 0.000000 1.000000
+
+1.1441229 0.833406 3.765340 3
+0.8312539 2.459607 -0.804545 4
+1.4142136 1.554275 1.724796 5
+0.0 0.0 0.0 1