Could you give an example?

Hi,

Following is a simplified example. I defined new_inverse method that have the fix, then patch the old inverse method with the new method. then a simple example is provided to demonstrate it works.

import numpy as np
from nibabel.affines import to_matvec, from_matvec
from nipy.core.reference.coordinate_map import AffineTransform
from nipy.core.reference.coordinate_system import CoordinateSystem

# define the inverse() method with the proposed solution
def new_inverse(self, preserve_dtype=False):
    aff_dt = self.function_range.coord_dtype
    try:
        # shawleez: consider non-squared affine
        if self.affine.shape[0] == self.affine.shape[1]:
            m_inv = npl.inv(self.affine)
        else:
            # shawleez: inversion for non-square affine matrix
            mat, vec = to_matvec(self.affine)
            invmat = np.linalg.pinv(mat)
            invvec = invmat.dot(-vec)
            m_inv = from_matvec(invmat, invvec)
    except npl.LinAlgError:
        return None
    except TypeError:
        # Try using sympy for the inverse.  This might be needed for sympy
        # symbols in the affine, or Float128, or complex numbers.
        from sympy import Matrix, matrix2numpy
        sym_inv = Matrix(self.affine).inv()
        m_inv = matrix2numpy(sym_inv).astype(aff_dt)
    else:  # linalg inverse succeeded
        if preserve_dtype and aff_dt != np.object:  # can we cast back?
            m_inv_orig = m_inv
            m_inv = m_inv.astype(aff_dt)
            if not np.allclose(m_inv_orig, m_inv):
                return None
    return AffineTransform(self.function_range,
                           self.function_domain,
                           m_inv)

# patch the inverse() method by the new method
AffineTransform.inverse = new_inverse

# non-squared 4x3 affine matrix that maps 'ij' to 'xyz'
affine = np.array([[1,0,1],[0,1,1], [0,0,1],[0,0,1]])
inp_cs = CoordinateSystem('ij')
out_cs = CoordinateSystem('xyz')
cm = AffineTransform(inp_cs, out_cs, affine)
# the inverse AffineTransform that map 'xyz' back to 'ij'
invcm = cm.inverse()

# example original coordinates in 'ij'
orig_ij2d = [0.,0.]
# map to 'xyz'
xyz3d = cm(orig_ij2d)
# map back to 'ij', the coordinates shall be equal to the original 'ij' coordinates
ij2d = invcm(xyz3d)
print(f"orig_ij2d = {orig_ij2d}")
print(f"xyz3d = {xyz3d}")
print(f"ij2d = {ij2d}")

AffineTransform does not work for 2D image embedded in 3D space #450

Description

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions