Skip to content

Commit b2e2c13

Browse files
matthew-brettmatthew-brett
committed
RF: refactor to restore previous inverse behavior
Restore the previous behavior of returning a floating point inverse (mostly) from an int input AffineTransform. Add flag to inverse method to specify new behavior of returning None when no integer inverse can be found.
1 parent 5dc0d33 commit b2e2c13

File tree

2 files changed

+124
-23
lines changed

2 files changed

+124
-23
lines changed

‎nipy/core/reference/coordinate_map.py

Lines changed: 67 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -573,16 +573,69 @@ def __init__(self, function_domain, function_range, affine):
573573
#
574574
###################################################################
575575

576-
def inverse(self):
576+
def inverse(self, preserve_dtype=False):
577577
""" Return coordinate map with inverse affine transform or None
578578
579-
Try to invert our ``affine``, and see if it can be cast to the needed
580-
data type, which is ``self.function_domain.coord_dtype``. We need this
581-
dtype in order for the inverse to preserve the coordinate system dtypes.
579+
Parameters
580+
----------
581+
preserve_dtype : bool
582+
If False, return affine mapping from inverting the ``affine``. The
583+
domain / range dtypes for the inverse may then change as a function
584+
of the dtype of the inverted ``affine``. If True, try to invert our
585+
``affine``, and see if it can be cast to the needed data type, which
586+
is ``self.function_domain.coord_dtype``. We need this dtype in
587+
order for the inverse to preserve the coordinate system dtypes.
588+
589+
Returns
590+
-------
591+
aff_cm_inv : ``AffineTransform`` instance or None
592+
``AffineTransform`` mapping from the *range* of input `self` to the
593+
*domain* of input `self` - the inverse of `self`. If
594+
``self.affine`` was not invertible return None. If `preserve_dtype`
595+
is True, and the inverse of ``self.affine`` cannot be cast to
596+
``self.function_domain.coord_dtype``, then return None. Otherwise
597+
return ``AffineTransform`` inverse mapping. If `preserve_dtype` is
598+
False, the domain / range dtypes of the return inverse may well be
599+
different from those of the input `self`.
600+
601+
Examples
602+
--------
603+
>>> input_cs = CoordinateSystem('ijk', coord_dtype=np.int)
604+
>>> output_cs = CoordinateSystem('xyz', coord_dtype=np.int)
605+
>>> affine = np.array([[1,0,0,1],
606+
... [0,1,0,1],
607+
... [0,0,1,1],
608+
... [0,0,0,1]])
609+
>>> affine_transform = AffineTransform(input_cs, output_cs, affine)
610+
>>> affine_transform([2,3,4]) #doctest: +IGNORE_DTYPE
611+
array([3, 4, 5])
612+
613+
The inverse transform, by default, generates a floating point inverse
614+
matrix and therefore floating point output:
615+
616+
>>> affine_transform_inv = affine_transform.inverse()
617+
>>> affine_transform_inv([2, 6, 12])
618+
array([ 1., 5., 11.])
619+
620+
You can force it to preserve the coordinate system dtype with the
621+
`preserve_dtype` flag:
622+
623+
>>> at_inv_preserved = affine_transform.inverse(preserve_dtype=True)
624+
>>> at_inv_preserved([2, 6, 12]) #doctest: +IGNORE_DTYPE
625+
array([ 1, 5, 11])
626+
627+
If you `preserve_dtype`, and there is no inverse affine preserving the
628+
dtype, the inverse is None:
629+
630+
>>> affine2 = affine.copy()
631+
>>> affine2[0, 0] = 2 # now inverse can't be integer
632+
>>> aff_t = AffineTransform(input_cs, output_cs, affine2)
633+
>>> aff_t.inverse(preserve_dtype=True) is None
634+
True
582635
"""
583636
aff_dt = self.function_range.coord_dtype
584637
try:
585-
m_inv_in = npl.inv(self.affine)
638+
m_inv = npl.inv(self.affine)
586639
except npl.LinAlgError:
587640
return None
588641
except TypeError:
@@ -591,10 +644,12 @@ def inverse(self):
591644
from sympy import Matrix, matrix2numpy
592645
sym_inv = Matrix(self.affine).inv()
593646
m_inv = matrix2numpy(sym_inv).astype(aff_dt)
594-
else: # linalg inverse succeeded - can we cast back?
595-
m_inv = m_inv_in.astype(aff_dt)
596-
if (aff_dt != np.object) and not np.allclose(m_inv_in, m_inv):
597-
return None
647+
else: # linalg inverse succeeded
648+
if preserve_dtype and aff_dt != np.object: # can we cast back?
649+
m_inv_orig = m_inv
650+
m_inv = m_inv.astype(aff_dt)
651+
if not np.allclose(m_inv_orig, m_inv):
652+
return None
598653
return AffineTransform(self.function_range,
599654
self.function_domain,
600655
m_inv)
@@ -900,9 +955,6 @@ def __call__(self, x):
900955
>>> affine_transform = AffineTransform(input_cs, output_cs, affine)
901956
>>> affine_transform([2,3,4]) #doctest: +IGNORE_DTYPE
902957
array([3, 4, 5])
903-
>>> affine_transform_inv = affine_transform.inverse()
904-
>>> affine_transform_inv([2, 6, 12])
905-
array([ 1, 5, 11])
906958
"""
907959
x = np.asanyarray(x)
908960
out_shape = (self.function_range.ndim,)
@@ -1567,7 +1619,9 @@ def _function(x):
15671619
value += b
15681620
return value
15691621

1570-
affine_transform_inv = affine_transform.inverse()
1622+
# Preserve dtype check because the CoordinateMap expects to generate the
1623+
# expected dtype and checks this on object creation
1624+
affine_transform_inv = affine_transform.inverse(preserve_dtype=True)
15711625
if affine_transform_inv:
15721626
Ainv, binv = to_matvec(affine_transform_inv.affine)
15731627
def _inverse_function(x):

‎nipy/core/reference/tests/test_coordinate_map.py

Lines changed: 57 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -938,22 +938,69 @@ def test_make_cmap():
938938
def test_dtype_cmap_inverses():
939939
# Check that we can make functional inverses of AffineTransforms, and
940940
# CoordinateMap versions of AffineTransforms
941-
dtypes = (np.sctypes['int'] + np.sctypes['float'] + np.sctypes['complex'] +
942-
[np.object])
941+
dtypes = (np.sctypes['int'] + np.sctypes['uint'] + np.sctypes['float']
942+
+ np.sctypes['complex'] + [np.object])
943943
arr_p1 = np.eye(4)[:, [0, 2, 1, 3]]
944-
in_arr = [0, 1, 2]
945-
out_arr = [0, 2, 1]
944+
in_list = [0, 1, 2]
945+
out_list = [0, 2, 1]
946946
for dt in dtypes:
947947
in_cs = CoordinateSystem('ijk', coord_dtype=dt)
948948
out_cs = CoordinateSystem('xyz', coord_dtype=dt)
949949
cmap = AffineTransform(in_cs, out_cs, arr_p1.astype(dt))
950+
coord = np.array(in_list, dtype=dt)
951+
out_coord = np.array(out_list, dtype=dt)
952+
# Expected output type of inverse, not preserving
953+
if dt in np.sctypes['int'] + np.sctypes['uint']:
954+
exp_i_dt = np.float64
955+
else:
956+
exp_i_dt = dt
957+
# Default inverse cmap may alter coordinate types
950958
r_cmap = cmap.inverse()
951-
assert_array_equal(r_cmap.affine, arr_p1)
959+
res = r_cmap(out_coord)
960+
assert_array_equal(res, coord)
961+
assert_equal(res.dtype, exp_i_dt)
962+
# Default behavior is preserve_type=False
963+
r_cmap = cmap.inverse(preserve_dtype=False)
964+
res = r_cmap(out_coord)
965+
assert_array_equal(res, coord)
966+
assert_equal(res.dtype, exp_i_dt)
967+
# Preserve_dtype=True - preserves dtype
968+
r_cmap = cmap.inverse(preserve_dtype=True)
969+
res = r_cmap(out_coord)
970+
assert_array_equal(res, coord)
971+
assert_equal(res.dtype, dt)
972+
# Preserve_dtype=True is default for conversion to CoordinateMap
952973
cm_cmap = _as_coordinate_map(cmap)
953-
coord = np.array(in_arr, dtype=dt)
954-
assert_array_equal(cm_cmap(coord), out_arr)
974+
assert_array_equal(cm_cmap(coord), out_list)
955975
rcm_cmap = cm_cmap.inverse()
956-
assert_array_equal(rcm_cmap(coord), out_arr)
976+
assert_array_equal(rcm_cmap(coord), out_list)
977+
res = rcm_cmap(out_coord)
978+
assert_array_equal(res, coord)
979+
assert_equal(res.dtype, dt)
980+
# For integer types, where there is no integer inverse, return floatey
981+
# inverse by default, and None for inverse when preserve_dtype=True
982+
arr_p2 = arr_p1 * 2
983+
arr_p2[-1, -1] = 1
984+
out_list = [0, 4, 2]
985+
for dt in np.sctypes['int'] + np.sctypes['uint']:
986+
in_cs = CoordinateSystem('ijk', coord_dtype=dt)
987+
out_cs = CoordinateSystem('xyz', coord_dtype=dt)
988+
cmap = AffineTransform(in_cs, out_cs, arr_p2.astype(dt))
989+
coord = np.array(in_list, dtype=dt)
990+
out_coord = np.array(out_list, dtype=dt)
991+
# Default
992+
r_cmap = cmap.inverse()
993+
res = r_cmap(out_coord)
994+
assert_array_equal(res, coord)
995+
assert_equal(res.dtype, np.float64)
996+
# Default is preserve_type=False
997+
r_cmap = cmap.inverse(preserve_dtype=False)
998+
res = r_cmap(out_coord)
999+
assert_array_equal(res, coord)
1000+
assert_equal(res.dtype, np.float64)
1001+
# preserve_dtype=True means there is no valid inverse for non integer
1002+
# affine inverses, as here
1003+
assert_equal(cmap.inverse(preserve_dtype=True), None)
9571004

9581005

9591006
def test_subtype_equalities():
@@ -977,8 +1024,8 @@ def test_cmap_coord_types():
9771024
# Check that we can use full range of coordinate system types. The inverse
9781025
# of an AffineTransform should generate coordinates in the input coordinate
9791026
# system dtype
980-
dtypes = (np.sctypes['int'] + np.sctypes['float'] + np.sctypes['complex'] +
981-
[np.object])
1027+
dtypes = (np.sctypes['int'] + np.sctypes['uint'] + np.sctypes['float']
1028+
+ np.sctypes['complex'] + [np.object])
9821029
arr_p1 = np.eye(4)
9831030
arr_p1[:3, 3] = 1
9841031
for dt in dtypes:

0 commit comments

Comments
 (0)