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

No description provided.

Could you give an example?

I created a 2D nipy Image that map 2D 'i,j' to 3D world coordinate 'x, y, z'. The corresponding AffineTransform has affine matrix size 3x4. However, it turns out the inverse() method in AffineTransform does not work for 3x4 matrix. As a result, i cannot use the amazing functions such as resample().

Following codes can fix the problem:

in AffineTransform.inverse():

def 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-squared 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:    
                ...
        ...

Please help to fix the issue such that AffineTransorm can be applied to 2D image embedded in 3D space.

Please let me know if you need more details. I am new for github, but very willing to support your great works!

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}")

Navigation Menu

Repository navigation

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

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

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

Description

Activity

matthew-brett commented on Apr 4, 2019

shawleez commented on Apr 4, 2019

shawleez commented on Apr 7, 2019

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions