transforms

`navis.transforms.AffineTransform` #

Affine transformation of 3D spatial data.

PARAMETER	DESCRIPTION
`matrix`	`Affine matrix.` TYPE: `(4, 4) np.ndarray`

Examples:

A simple scaling transform

>>> from navis import transforms
>>> import numpy as np
>>> M = np.diag([1e3, 1e3, 1e3, 1])
>>> tr = transforms.affine.AffineTransform(M)
>>> points = np.array([[0, 0, 0], [1, 1, 1]])
>>> tr.xform(points)
array([[   0.,    0.,    0.],
       [1000., 1000., 1000.]])

Source code in navis/transforms/affine.py

class AffineTransform(BaseTransform):
    """Affine transformation of 3D spatial data.

    Parameters
    ----------
    matrix :        (4, 4) np.ndarray
                    Affine matrix.

    Examples
    --------
    A simple scaling transform

    >>> from navis import transforms
    >>> import numpy as np
    >>> M = np.diag([1e3, 1e3, 1e3, 1])
    >>> tr = transforms.affine.AffineTransform(M)
    >>> points = np.array([[0, 0, 0], [1, 1, 1]])
    >>> tr.xform(points)
    array([[   0.,    0.,    0.],
           [1000., 1000., 1000.]])

    """

    def __init__(self, matrix: np.ndarray, direction: str = 'forward'):
        """Initialize transform."""
        assert direction in ('forward', 'inverse')

        self.matrix = matrix

        if direction == 'inverse':
            self.matrix = np.linalg.inv(self.matrix)

    def __eq__(self, other: 'AffineTransform') -> bool:
        """Implements equality comparison."""
        if isinstance(other, AffineTransform):
            if np.all(self.matrix == other.matrix):
                return True
        return False

    def __neg__(self):
        """Invert direction."""
        x = self.copy()

        # Invert affine matrix
        x.matrix = np.linalg.inv(x.matrix)

        return x

    def __add__(self, other: 'AffineTransform') -> 'AffineTransform':
        """Add two affine transforms."""
        if isinstance(other, AffineTransform):
            x = self.copy()
            x.matrix = np.dot(self.matrix, other.matrix)
            return x
        raise ValueError('Can only add `AffineTransform` objects.')

    def copy(self) -> 'AffineTransform':
        """Return copy of transform."""
        # Attributes not to copy
        no_copy = []
        # Generate new empty transform
        x = self.__class__(None)
        # Override with this neuron's data
        x.__dict__.update({k: copy.copy(v) for k, v in self.__dict__.items() if k not in no_copy})

        return x

    def xform(self, points: np.ndarray, invert: bool = False) -> np.ndarray:
        """Apply transform to points.

        Parameters
        ----------
        points :        np.ndarray
                        (N, 3) array of x/y/z locations.
        invert :        bool
                        If True, will invert the transform.

        Returns
        -------
        pointsxf :      np.ndarray
                        The transformed points.

        """
        points = np.asarray(points)

        if points.ndim != 2 or points.shape[1] != 3:
            raise ValueError('`points` must be of shape (N, 3)')

        # Add a fourth column to points
        points_mat = np.ones((points.shape[0], 4))
        points_mat[:, :3] = points

        # Apply transform
        if not invert:
            mat = self.matrix
        else:
            mat = np.linalg.inv(self.matrix)

        return np.dot(mat, points_mat.T).T[:, :3]

`init` #

Initialize transform.

Source code in navis/transforms/affine.py

def __init__(self, matrix: np.ndarray, direction: str = 'forward'):
    """Initialize transform."""
    assert direction in ('forward', 'inverse')

    self.matrix = matrix

    if direction == 'inverse':
        self.matrix = np.linalg.inv(self.matrix)

`copy` #

Return copy of transform.

Source code in navis/transforms/affine.py

def copy(self) -> 'AffineTransform':
    """Return copy of transform."""
    # Attributes not to copy
    no_copy = []
    # Generate new empty transform
    x = self.__class__(None)
    # Override with this neuron's data
    x.__dict__.update({k: copy.copy(v) for k, v in self.__dict__.items() if k not in no_copy})

    return x

`xform` #

Apply transform to points.

PARAMETER	DESCRIPTION
`points`	`(N, 3) array of x/y/z locations.` TYPE: `np.ndarray`
`invert`	`If True, will invert the transform.` TYPE: `bool` DEFAULT: `False`

RETURNS	DESCRIPTION
`pointsxf`	The transformed points. TYPE: `np.ndarray`

Source code in navis/transforms/affine.py

def xform(self, points: np.ndarray, invert: bool = False) -> np.ndarray:
    """Apply transform to points.

    Parameters
    ----------
    points :        np.ndarray
                    (N, 3) array of x/y/z locations.
    invert :        bool
                    If True, will invert the transform.

    Returns
    -------
    pointsxf :      np.ndarray
                    The transformed points.

    """
    points = np.asarray(points)

    if points.ndim != 2 or points.shape[1] != 3:
        raise ValueError('`points` must be of shape (N, 3)')

    # Add a fourth column to points
    points_mat = np.ones((points.shape[0], 4))
    points_mat[:, :3] = points

    # Apply transform
    if not invert:
        mat = self.matrix
    else:
        mat = np.linalg.inv(self.matrix)

    return np.dot(mat, points_mat.T).T[:, :3]

`navis.transforms.AliasTransform` #

Helper transform that simply passes points through.

Useful for defining aliases.

Source code in navis/transforms/base.py

class AliasTransform(BaseTransform):
    """Helper transform that simply passes points through.

    Useful for defining aliases.
    """

    def __init__(self):
        """Initialize."""
        pass

    def __neg__(self) -> "AliasTransform":
        """Invert transform."""
        return self.copy()

    def __eq__(self, other):
        """Check if the same."""
        if isinstance(other, AliasTransform):
            return True
        return False

    def copy(self):
        """Return copy."""
        x = AliasTransform()
        x.__dict__.update(self.__dict__)
        return x

    def xform(self, points: np.ndarray) -> np.ndarray:
        """Pass through.

        Note that the returned points are NOT a copy but the originals.
        """
        return points

`init` #

Initialize.

Source code in navis/transforms/base.py

def __init__(self):
    """Initialize."""
    pass

`copy` #

Return copy.

Source code in navis/transforms/base.py

def copy(self):
    """Return copy."""
    x = AliasTransform()
    x.__dict__.update(self.__dict__)
    return x

`xform` #

Pass through.

Note that the returned points are NOT a copy but the originals.

Source code in navis/transforms/base.py

def xform(self, points: np.ndarray) -> np.ndarray:
    """Pass through.

    Note that the returned points are NOT a copy but the originals.
    """
    return points

`navis.transforms.CMTKtransform` #

CMTK transforms of 3D spatial data.

Requires either CMTK to be installed, or navis-fastcore - the latter implements CMTK transforms in Rust and needs no external binaries. See the backend parameter.

PARAMETER	DESCRIPTION
`regs`	`Path(s) to CMTK transformations(s).` TYPE: `str \| list of str`
`directions`	Direction of transformation. Must provide one direction per `reg`. TYPE: `"forward" \| "inverse" \| list thereof` DEFAULT: `'forward'`
`threads`	`Number of threads to use.` TYPE: `int` DEFAULT: `None`
`backend`	Which implementation to use. `None` (default) defers to `navis.config.default_transform_backend`. "binary" is deprecated and will be removed in 3.0. TYPE: `"auto" \| "fastcore" \| "binary"` DEFAULT: `None`

Examples:

>>> from navis import transforms
>>> tr = transforms.cmtk.CMTKtransform('/path/to/CMTK_directory.list')
>>> tr.xform(points)

Notes

TODO (when navis-fastcore becomes a hard dependency and the "binary" backend goes away): switch affine_fallback to fastcore's fallback_to_affine="hop".

We currently ask fastcore for fallback_to_affine=True, i.e. "re-run a failed point affine-only through the whole chain, from where it started". That is what streamxform --affine-only does, and while we support both backends they have to agree - so it is the right choice for now.

But it is the cruder of the two. Where several registrations are merged into one transform (which happens on about half of all bridging paths), a point that clears the first few warps and then leaves the domain loses those warps too. fastcore's "hop" keeps them and swaps the affine in only for the hop that actually failed. On a 2-registration chain the two differ on ~19% of points, by a median of ~2 world units.

Once there is no binary to stay consistent with, "hop" is the better answer.

Source code in navis/transforms/cmtk.py

class CMTKtransform(BackendMixin, BaseTransform):
    """CMTK transforms of 3D spatial data.

    Requires either [CMTK](https://www.nitrc.org/projects/cmtk/) to be installed,
    or navis-fastcore - the latter implements CMTK transforms in Rust and needs
    no external binaries. See the `backend` parameter.

    Parameters
    ----------
    regs :          str | list of str
                    Path(s) to CMTK transformations(s).
    directions :    "forward" | "inverse" | list thereof
                    Direction of transformation. Must provide one direction per
                    `reg`.
    threads :       int, optional
                    Number of threads to use.
    backend :       "auto" | "fastcore" | "binary", optional
                    Which implementation to use. `None` (default) defers to
                    `navis.config.default_transform_backend`. "binary" is
                    deprecated and will be removed in 3.0.

    Examples
    --------
    >>> from navis import transforms
    >>> tr = transforms.cmtk.CMTKtransform('/path/to/CMTK_directory.list')
    >>> tr.xform(points) # doctest: +SKIP

    Notes
    -----
    TODO (when navis-fastcore becomes a hard dependency and the "binary" backend
    goes away): switch `affine_fallback` to fastcore's `fallback_to_affine="hop"`.

    We currently ask fastcore for `fallback_to_affine=True`, i.e. "re-run a failed
    point affine-only through the *whole* chain, from where it started". That is
    what `streamxform --affine-only` does, and while we support both backends they
    have to agree - so it is the right choice for now.

    But it is the cruder of the two. Where several registrations are merged into
    one transform (which happens on about half of all bridging paths), a point that
    clears the first few warps and then leaves the domain loses those warps too.
    fastcore's `"hop"` keeps them and swaps the affine in only for the hop that
    actually failed. On a 2-registration chain the two differ on ~19% of points, by
    a median of ~2 world units.

    Once there is no binary to stay consistent with, `"hop"` is the better answer.

    """

    # A CMTK warp has no closed-form inverse: going backwards means solving for
    # each point iteratively. That is ~2-4x slower than going forwards and fails on
    # a different (usually larger) set of points, which come back as NaN.
    #
    # 2 means navis will take one extra forward hop rather than invert a CMTK
    # registration once. Kept modest deliberately: unlike elastix, CMTK
    # registrations rarely ship a purpose-built reverse, so the graph genuinely
    # depends on inverting them - roughly three quarters of all bridging paths do.
    # Weighting them heavily would just push routes into long detours.
    inverse_weight_factor = 2

    def __init__(
        self,
        regs: list,
        directions: str = "forward",
        threads: int = None,
        backend: str = None,
    ):
        self._backend = backend
        self.directions = list(utils.make_iterable(directions))
        for d in self.directions:
            assert d in ("forward", "inverse"), (
                '`direction` must be "foward"' f'or "inverse", not "{d}"'
            )

        self.regs = list(utils.make_iterable(regs))
        self.command = "streamxform"
        self.threads = threads

        # Broadcast a single direction to all registrations
        if len(self.directions) == 1 and len(self.regs) >= 1:
            self.directions = self.directions * len(self.regs)

        if len(self.regs) != len(self.directions):
            raise ValueError("Must provide one direction per regs")

    def __add__(self, other: "CMTKtransform") -> "CMTKtransform":
        """Implement addition operator to concatenate transforms."""
        if not isinstance(other, CMTKtransform):
            raise TypeError(f"Cannot add {type(other)} to CMTKtransform")

        if self.command != other.command:
            raise ValueError("Unable to merge CMTKtransforms using different commands.")

        if self.backend != other.backend:
            raise ValueError(
                "Unable to merge CMTKtransforms using different backends."
            )

        x = self.copy()
        x.regs += other.regs
        x.directions += other.directions

        return x

    def __eq__(self, other: "CMTKtransform") -> bool:
        """Implement equality comparison."""
        if isinstance(other, CMTKtransform):
            if len(self) == len(other):
                if all([self.regs[i] == other.regs[i] for i in range(len(self))]):
                    if all(
                        [
                            self.directions[i] == other.directions[i]
                            for i in range(len(self))
                        ]
                    ):
                        return True
        return False

    def __len__(self) -> int:
        return len(self.regs)

    def __neg__(self) -> "CMTKtransform":
        """Invert direction."""
        x = self.copy()

        # Swap directions
        x.directions = [
            {"forward": "inverse", "inverse": "forward"}[d] for d in x.directions
        ]

        # Reverse order
        x.regs = x.regs[::-1]
        x.directions = x.directions[::-1]

        return x

    def __str__(self):
        return self.__repr__()

    def __repr__(self):
        return f"CMTKtransform with {len(self)} transform(s)"

    @staticmethod
    def from_file(filepath: str, **kwargs) -> "CMTKtransform":
        """Generate CMTKtransform from file.

        Parameters
        ----------
        filepath :  str
                    Path to CMTK transform.
        **kwargs
                    Keyword arguments passed to CMTKtransform.__init__

        Returns
        -------
        CMTKtransform

        """
        defaults = {"directions": "forward"}
        defaults.update(kwargs)
        return CMTKtransform(str(filepath), **defaults)

    def make_args(self, affine_only: bool = False) -> list:
        """Generate arguments passed to subprocess."""
        # Generate the arguments
        # The actual command (i.e. streamxform)
        args = [str(_cmtkbin / self.command)]

        if affine_only:
            args.append("--affine-only")

        if self.threads:
            args.append(f"--threads {int(self.threads)}")

        # Add the regargs
        args += self.regargs

        return args

    @property
    def regargs(self) -> list:
        """Generate regargs."""
        regargs = []
        for i, (reg, dir) in enumerate(zip(self.regs, self.directions)):
            if dir == "inverse":
                # For the first transform we need to prefix "--inverse" with
                # a solitary "--"
                if i == 0:
                    regargs.append("--")
                regargs.append("--inverse")
            # Note no double quotes!
            regargs.append(f"{reg}")
        return regargs

    def append(self, transform: "CMTKtransform", direction: str = None):
        """Add another transform.

        Parameters
        ----------
        transform :     str | CMTKtransform
                        Either another CMTKtransform or filepath to registration.
        direction :     "forward" | "inverse"
                        Only relevant if transform is filepath.

        """
        if isinstance(transform, CMTKtransform):
            if self.command != transform.command:
                raise ValueError(
                    "Unable to merge CMTKtransforms using " "different commands."
                )

            # Note we must raise a NotImplementedError (rather than e.g. a
            # ValueError) if we can't merge: that's the signal TransformSequence
            # uses to chain the transforms instead.
            if self.backend != transform.backend:
                raise NotImplementedError(
                    "Unable to merge CMTKtransforms using different backends."
                )

            self.regs += transform.regs
            self.directions += transform.directions
        elif isinstance(transform, str):
            if not direction:
                raise ValueError("Must provide direction along with new transform")
            self.regs.append(transform)
            self.directions.append(direction)
        else:
            raise NotImplementedError(
                f"Unable to append {type(transform)} to {type(self)}"
            )

    def check_if_possible(self, on_error: str = "raise"):
        """Check if this transform is possible."""
        msg = None
        # Note the binaries are only required for the "binary" backend - the
        # whole point of the fastcore backend is that it doesn't need them.
        if self.backend == "binary" and not _cmtkbin:
            msg = (
                "Folder with CMTK binaries not found. Make sure the "
                "directory is in your PATH environment variable - or install "
                "navis-fastcore to transform points without CMTK."
            )
        else:
            for r in self.regs:
                if not os.path.isdir(r) and not os.path.isfile(r):
                    msg = f"Registration {r} not found."
                    break

        if msg and on_error == "raise":
            raise BaseException(msg)

        return msg

    def copy(self) -> "CMTKtransform":
        """Return copy."""
        return self.__class__(
            regs=copy.copy(self.regs),
            directions=copy.copy(self.directions),
            threads=self.threads,
            backend=self._backend,
        )

    def parse_cmtk_output(self, output: str, fail_value=np.nan) -> np.ndarray:
        r"""Parse CMTK output.

        Briefly, CMTK output will be a byte literal like this:

            b'311 63 23 \n275 54 25 \n'

        In case of failed transforms we will get something like this where the
        original coordinates are returned with a "FAILED" flag

            b'343 72 23 \n-10 -10 -10 FAILED \n'

        Parameter
        ---------
        output :        tuple of (b'', None)
                        Stdout of CMTK call.
        fail_value
                        Value to use for points that failed to transform. By
                        default we use `np.nan`.

        Returns
        -------
        pointsxf :      (N, 3) numpy array
                        The parse transformed points.

        """
        # The original stout is tuple where we care only about the second one
        if isinstance(output, tuple):
            output = output[0]

        pointsx = []
        # Split string into rows - lazily using a generator
        for row in (x.group(1) for x in re.finditer(r"(.*?) \n", output.decode())):
            # Split into values
            values = row.split(" ")

            # If this point failed
            if len(values) != 3:
                values = [fail_value] * 3
            else:
                values = [float(v) for v in values]

            pointsx.append(values)

        return np.asarray(pointsx)

    def xform(
        self,
        points: np.ndarray,
        affine_only: bool = False,
        affine_fallback: bool = False,
    ) -> np.ndarray:
        """Xform data.

        Parameters
        ----------
        points :            (N, 3) numpy array | pandas.DataFrame
                            Points to xform. DataFrame must have x/y/z columns.
        affine_only :       bool
                            Whether to apply only the non-rigid affine
                            transform. This is useful if points are outside
                            the deformation field and would therefore not
                            transform properly.
        affine_fallback :   bool
                            If True and some points did not transform during the
                            non-rigid part of the transformation, we will apply
                            only the affine transformation to those points.

        Returns
        -------
        pointsxf :      (N, 3) numpy array
                        Transformed points. Points that failed to transform will
                        be `np.nan`.

        """
        self.check_if_possible(on_error="raise")

        if self.backend == "fastcore":
            # One parsed object serves every direction - fastcore takes the
            # per-hop direction at xform time, and our `directions` map straight
            # onto its `invert` flags.
            #
            # `fallback_to_affine=True` means "re-run the whole chain affine-only
            # from the original point" - what streamxform --affine-only does, and
            # so what we need while the binary backend still exists.
            # See the class Notes: switch this to "hop" once it doesn't.
            reg = get_cmtk_reg(tuple(str(r) for r in self.regs))
            return reg.xform(
                parse_points(points),
                transform="affine" if affine_only else "warp",
                fallback_to_affine=bool(affine_fallback) and not affine_only,
                invert=[d == "inverse" for d in self.directions],
                n_cores=self.threads,
            )

        if isinstance(points, pd.DataFrame):
            # Make sure x/y/z columns are present
            if np.any([c not in points for c in ["x", "y", "z"]]):
                raise ValueError("points DataFrame must have x/y/z columns.")
        elif (
            isinstance(points, np.ndarray) and points.ndim == 2 and points.shape[1] == 3
        ):
            points = pd.DataFrame(points, columns=["x", "y", "z"])
        else:
            raise TypeError(
                "`points` must be numpy array of shape (N, 3) or "
                "pandas DataFrame with x/y/z columns"
            )

        # Generate the result
        args = self.make_args(affine_only=affine_only)
        proc = subprocess.Popen(args, stdin=subprocess.PIPE, stdout=subprocess.PIPE)

        # Pipe in the points
        points_str = points[["x", "y", "z"]].to_string(index=False, header=False)

        # Do not use proc.stdin.write to avoid output buffer becoming full
        # before we finish piping in stdin.
        # proc.stdin.write(points_str.encode())
        # output = proc.communicate()

        # Read out results
        # This is equivalent to e.g.:
        # $ streamxform -args <<< "10, 10, 10"
        output = proc.communicate(input=points_str.encode())

        # If no output, something went wrong
        if not output[0]:
            raise utils.CMTKError("CMTK produced no output. Check points?")

        # Xformed points
        xf = self.parse_cmtk_output(output, fail_value=np.nan)

        # Check if any points not xformed
        if affine_fallback and not affine_only:
            not_xf = np.any(np.isnan(xf), axis=1)
            if np.any(not_xf):
                xf[not_xf] = self.xform(points.loc[not_xf], affine_only=True)

        return xf

    def to_grid_transform(self, template, absolute: bool = True, verbose: bool = False):
        """Convert to GridTransform via dense deformation field.

        Parameters
        ----------
        template :  str | TemplateBrain | (Nx, Ny, Nz, dx, dy, dz) tuple
                    This defines the bounds and voxel size of the deformation
                    field. Typically, this would correspond to the source space
                    (i.e. the moving image). We can work with:
                      - str: a filepath to a NRRD file
                      - TemplateBrain: a navis TemplateBrain object
                      - a tuple/list/array with (Nx, Ny, Nz, dx, dy, dz)
                        where N is the number of voxels in each dimension
                        and d is the voxel size.
        absolute :  bool
                    Whether to return absolute coordinates or offsets.
        verbose :   bool
                    Whether to print CMTK output.

        Returns
        -------
        transform:  GridTransform
                    A few notes on using the resulting GridTransform:
                      1. The transform will expect input coordinates in voxel space.
                         See the returned `voxel_size` array!
                      2. The transformed coordinates, on the other hand, will already
                         be in physical space (e.g. microns).
        voxel_size :   (3,) numpy array
                    The voxel size of the input space. This is important because
                    the GridTransform will expect input coordinates in voxel
                    space.
        """
        from .grid import GridTransform

        dfield, header = self.to_dfield(
            template, out=None, absolute=absolute, verbose=verbose
        )

        # Generate the transform
        transform = GridTransform(
            np.transpose(dfield, (1, 2, 3, 0)),
            type="coordinates" if absolute else "offsets",
        )
        # Get voxel size from header - note that the first row is for the 4th dimension
        # (i.e. the vector dimension of the offsets/coordinates) which we don't care about
        space_dir = np.diagonal(header["space directions"][1:])

        return transform, space_dir

    @requires_cmtk
    def to_dfield(
        self, template, out=None, absolute: bool = True, verbose: bool = False
    ) -> np.ndarray:
        """Convert transform to dense deformation field.

        Note this always requires CMTK: navis-fastcore implements point but not
        image transforms.

        Parameters
        ----------
        template :  str | TemplateBrain | (Nx, Ny, Nz, dx, dy, dz) tuple
                    This defines the bounds and voxel size of the deformation
                    field. Typically, this would correspond to the source space
                    (i.e. the moving image). We can work with:
                      - str: a filepath to a NRRD file
                      - TemplateBrain: a navis TemplateBrain object
                      - a tuple/list/array with (Nx, Ny, Nz, dx, dy, dz)
                        where N is the number of voxels in each dimension
                        and d is the voxel size.
        out :       str, optional
                    NRRD filepath to save the deformation field. If None (default),
                    the deformation field will be returned as numpy array.
        absolute :  bool
                    Whether to return absolute coordinates or offsets.
        verbose :   bool
                    Whether to print CMTK output.

        Returns
        -------
        dfield :    np.ndarray
                    The dense deformation field as (3, Nx, Ny, Nz) numpy array
                    with either absolute coordinates or offsets where the
                    first dimension contains the x/y/z coordinates. Note that
                    the coordinates stored in the field are in physical space
                    (e.g. microns), not voxel space.
        header :    dict
                    The NRRD header associated with the deformation field.


        See Also
        --------
        CMTKtransform.to_grid_transform
                    Method to directly convert to GridTransform. Please see
                    that method's notes for details on how to handle the
                    deformation field.

        """
        from .templates import TemplateBrain  # avoid circular import

        # Translate template info into shape and voxel size
        if isinstance(template, TemplateBrain):
            if not hasattr(template, "dims") or not hasattr(template, "voxdims"):
                raise ValueError(
                    "TemplateBrain must have `dims` and `voxdims` attributes"
                )
            if len(template.dims) != 3 or len(template.voxdims) != 3:
                raise ValueError(
                    "TemplateBrain `dims` and `voxdims` must be of length 3"
                )
            template = list(template.dims) + list(template.voxdims)

        to_remove = []  # track temporary files to clean up later

        try:
            if isinstance(template, str):
                # Assume NRRD filepath
                if not os.path.isfile(template):
                    raise ValueError(f"Template file not found: {template}")
                ref_image = template
            elif isinstance(template, (tuple, list, np.ndarray)):
                if len(template) != 6:
                    raise ValueError(
                        "`template` tuple/list/array must be of length 6: "
                        "(Nx, Ny, Nz, dx, dy, dz)"
                    )

                # Create a temporary NRRD file to use as reference
                # (we're using zeros because that's highly compressible)
                im = np.zeros(np.array(template[:3]).astype(int), dtype=np.uint8)
                header = {
                    "space directions": np.diag(template[3:6]),
                    "kinds": ["domain", "domain", "domain"],
                    "space units": ["microns", "microns", "microns"],
                    "space origin": [0.0, 0.0, 0.0],
                    "labels": ["x", "y", "z"],
                    "space dimension": 3,
                }
                with tempfile.NamedTemporaryFile(suffix=".nrrd", delete=False) as tf:
                    nrrd.write(tf.name, im, header=header)
                    to_remove.append(tf.name)
                    ref_image = tf.name
            else:
                raise TypeError(
                    f"`template` must be str, TemplateBrain or tuple, not {type(template)}"
                )

            if out is None:
                outfile = tempfile.NamedTemporaryFile(suffix=".nrrd", delete=False).name
                to_remove.append(outfile)
            elif isinstance(out, (str, pathlib.Path)):
                outfile = pathlib.Path(out).resolve()
            else:
                raise ValueError(f"Invalid output type: {type(out)}")

            # Compile the command
            args = [str(_cmtkbin / "xform2dfield")]

            if absolute:
                args += ["--output-absolute"]

            args += [f"{outfile}"]
            args += [f"{ref_image}"]
            args += self.regargs

            # run the binary
            # avoid resourcewarnings with null
            with open(os.devnull, "w") as devnull:
                startupinfo = None
                if platform.system() == "Windows":
                    startupinfo = subprocess.STARTUPINFO()
                    startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW
                if verbose:
                    # in debug mode print the output
                    stdout = None
                else:
                    stdout = devnull

                if verbose:
                    config.logger.info("executing: {}".format(" ".join(args)))
                check_call(
                    args,
                    stdout=stdout,
                    stderr=subprocess.STDOUT,
                    startupinfo=startupinfo,
                )

            if out is None:
                # Return transformed image
                return nrrd.read(outfile)
            elif verbose:
                config.logger.info(f"Transformed image saved to {outfile}")

        except BaseException:
            raise
        finally:
            # Clean up temporary files
            for f in to_remove:
                os.remove(f)

    @requires_cmtk
    def xform_image(
        self,
        im,
        target,
        out=None,
        interpolation="linear",
        verbose=False,
    ):
        """Transform an image using CMTK's reformatx.

        Note this always requires CMTK: navis-fastcore implements point but not
        image transforms.

        Parameters
        ----------
        im :        3D numpy array | filepath
                    The floating image to transform.
        target :    str | TemplateBrain | (Nx, Ny, Nz, dx, dy, dz) | (Nx, Ny, Nz, dx, dy, dz, Ox, Oy, Oz)
                    Defines the target image: dimensions in voxels (N), the voxel size (d) and optionally
                    an origin (0) for the target image. Can be provided as a string (name of a template),
                    a TemplateBrain object, a tuple/list/array with the target specs.
        out :       str, optional
                    The filepath to save the transformed image. If None (default), will return the
                    transformed image as np.ndarray.
        interpolation : "linear" | "nn" | "cubic" | "pv" | "sinc-cosine" | "sinc-hamming"
                    The interpolation method to use.
        verbose :   bool
                    Whether to print CMTK output.

        Returns
        -------
        np.ndarray | None
                    If out is None, returns the transformed image as np.ndarray. Otherwise, None.

        """
        assert interpolation in (
            "linear",
            "nn",
            "cubic",
            "pv",
            "sinc-cosine",
            "sinc-hamming",
        )

        # `reformatx` expects this format:
        # ./reformatx --floating {INPUT_FILE} -o {OUTPUT_FILE} {REFERENCE_SPECS} {TRANSFORMS}
        # where:
        # - {INPUT_FILE} is the image to transform
        # - {OUTPUT_FILE} is where the output will be saved
        # - {REFERENCE_SPECS} defines the target space; this needs to be either a NRRD
        #   file from which CMTK can extract the target grid or the actual specs:
        #   "--target-grid Nx,Ny,Nz:dX,dY,dZ:[Ox,Oy,Oz]" where N is the number of
        #   voxels in each dimension and d is the voxel size. The optional O is the
        #   origin of the image. If not provided, it is assumed to be (0, 0, 0).
        # - {TRANSFORMS} are the CMTK transform(s) to apply; prefix with "--inverse" to invert
        # Below command works to convert JFRC2 to FCWB:
        # /opt/local/lib/cmtk/bin/reformatx --verbose --floating JFRC2.nrrd -o JFRC2_xf.nrrd FCWB.nrrd ~/flybrain-data/BridgingRegistrations/FCWB_JFRC2.list
        # This took 2min 28s - should check if that is actually faster than the look-up approach we use in `images.py`
        # Note that inversion of transforms always comes with an overhead! When using the inverse transform instead, the above command took 2h 40min!

        if verbose and any((d == "inverse " for d in self.directions)):
            config.logger.warning(
                "Using inverse CMTK transforms with reformatx can be very slow! If possible, consider using only forward transforms."
            )

        target_specs = parse_target_specs(target)

        to_remove = []
        if isinstance(im, (str, pathlib.Path)):
            floating = pathlib.Path(im)
            if not im.is_file():
                raise ValueError(f"Image file not found: {im}")
        elif isinstance(im, np.ndarray):
            assert im.ndim == 3
            # Save to temporary file
            with tempfile.NamedTemporaryFile(suffix=".nrrd", delete=False) as tf:
                nrrd.write(tf.name, im)
                floating = tf.name
                to_remove.append(tf.name)
        else:
            raise ValueError(f"Invalid image type: {type(im)}")

        if out is None:
            outfile = tempfile.NamedTemporaryFile(suffix=".nrrd", delete=False).name
            to_remove.append(outfile)
        elif isinstance(out, (str, pathlib.Path)):
            outfile = pathlib.Path(out).resolve()
        else:
            raise ValueError(f"Invalid output type: {type(out)}")

        # Compile the command
        args = [str(_cmtkbin / "reformatx")]
        args += [f"-o {outfile}"]
        args += [f"--floating {floating}"]
        args += [f"--{interpolation}"]
        args += [target_specs]

        # Add the regargs
        args += self.regargs

        try:
            # run the binary
            # avoid resourcewarnings with null
            with open(os.devnull, "w") as devnull:
                startupinfo = None
                if platform.system() == "Windows":
                    startupinfo = subprocess.STARTUPINFO()
                    startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW
                if verbose:
                    # in debug mode print the output
                    stdout = None
                else:
                    stdout = devnull

                if verbose:
                    config.logger.info("executing: {}".format(" ".join(args)))
                check_call(
                    args,
                    stdout=stdout,
                    stderr=subprocess.STDOUT,
                    startupinfo=startupinfo,
                )

            if out is None:
                # Return transformed image
                return nrrd.read(outfile)
            elif verbose:
                config.logger.info(f"Transformed image saved to {outfile}")
        except BaseException:
            raise
        finally:
            # Clean up temporary files
            for f in to_remove:
                os.remove(f)

`regargs: list` `property` #

Generate regargs.

`append` #

Add another transform.

PARAMETER	DESCRIPTION
`transform`	`Either another CMTKtransform or filepath to registration.` TYPE: `str \| CMTKtransform`
`direction`	`Only relevant if transform is filepath.` TYPE: `"forward" \| "inverse"` DEFAULT: `None`

Source code in navis/transforms/cmtk.py

def append(self, transform: "CMTKtransform", direction: str = None):
    """Add another transform.

    Parameters
    ----------
    transform :     str | CMTKtransform
                    Either another CMTKtransform or filepath to registration.
    direction :     "forward" | "inverse"
                    Only relevant if transform is filepath.

    """
    if isinstance(transform, CMTKtransform):
        if self.command != transform.command:
            raise ValueError(
                "Unable to merge CMTKtransforms using " "different commands."
            )

        # Note we must raise a NotImplementedError (rather than e.g. a
        # ValueError) if we can't merge: that's the signal TransformSequence
        # uses to chain the transforms instead.
        if self.backend != transform.backend:
            raise NotImplementedError(
                "Unable to merge CMTKtransforms using different backends."
            )

        self.regs += transform.regs
        self.directions += transform.directions
    elif isinstance(transform, str):
        if not direction:
            raise ValueError("Must provide direction along with new transform")
        self.regs.append(transform)
        self.directions.append(direction)
    else:
        raise NotImplementedError(
            f"Unable to append {type(transform)} to {type(self)}"
        )

`check_if_possible` #

Check if this transform is possible.

Source code in navis/transforms/cmtk.py

def check_if_possible(self, on_error: str = "raise"):
    """Check if this transform is possible."""
    msg = None
    # Note the binaries are only required for the "binary" backend - the
    # whole point of the fastcore backend is that it doesn't need them.
    if self.backend == "binary" and not _cmtkbin:
        msg = (
            "Folder with CMTK binaries not found. Make sure the "
            "directory is in your PATH environment variable - or install "
            "navis-fastcore to transform points without CMTK."
        )
    else:
        for r in self.regs:
            if not os.path.isdir(r) and not os.path.isfile(r):
                msg = f"Registration {r} not found."
                break

    if msg and on_error == "raise":
        raise BaseException(msg)

    return msg

`copy` #

Return copy.

Source code in navis/transforms/cmtk.py

def copy(self) -> "CMTKtransform":
    """Return copy."""
    return self.__class__(
        regs=copy.copy(self.regs),
        directions=copy.copy(self.directions),
        threads=self.threads,
        backend=self._backend,
    )

`from_file` `staticmethod` #

Generate CMTKtransform from file.

PARAMETER	DESCRIPTION
`filepath`	`Path to CMTK transform.` TYPE: `str`
`**kwargs`	`Keyword arguments passed to CMTKtransform.__init__` DEFAULT: `{}`

RETURNS	DESCRIPTION
`CMTKtransform`

Source code in navis/transforms/cmtk.py

@staticmethod
def from_file(filepath: str, **kwargs) -> "CMTKtransform":
    """Generate CMTKtransform from file.

    Parameters
    ----------
    filepath :  str
                Path to CMTK transform.
    **kwargs
                Keyword arguments passed to CMTKtransform.__init__

    Returns
    -------
    CMTKtransform

    """
    defaults = {"directions": "forward"}
    defaults.update(kwargs)
    return CMTKtransform(str(filepath), **defaults)

`make_args` #

Generate arguments passed to subprocess.

Source code in navis/transforms/cmtk.py

def make_args(self, affine_only: bool = False) -> list:
    """Generate arguments passed to subprocess."""
    # Generate the arguments
    # The actual command (i.e. streamxform)
    args = [str(_cmtkbin / self.command)]

    if affine_only:
        args.append("--affine-only")

    if self.threads:
        args.append(f"--threads {int(self.threads)}")

    # Add the regargs
    args += self.regargs

    return args

`parse_cmtk_output` #

Parse CMTK output.

Briefly, CMTK output will be a byte literal like this:

b'311 63 23 \n275 54 25 \n'

In case of failed transforms we will get something like this where the original coordinates are returned with a "FAILED" flag

b'343 72 23 \n-10 -10 -10 FAILED \n'

Parameter

output : tuple of (b'', None) Stdout of CMTK call. fail_value Value to use for points that failed to transform. By default we use np.nan.

RETURNS	DESCRIPTION
`pointsxf`	The parse transformed points. TYPE: `(N, 3) numpy array`

Source code in navis/transforms/cmtk.py

def parse_cmtk_output(self, output: str, fail_value=np.nan) -> np.ndarray:
    r"""Parse CMTK output.

    Briefly, CMTK output will be a byte literal like this:

        b'311 63 23 \n275 54 25 \n'

    In case of failed transforms we will get something like this where the
    original coordinates are returned with a "FAILED" flag

        b'343 72 23 \n-10 -10 -10 FAILED \n'

    Parameter
    ---------
    output :        tuple of (b'', None)
                    Stdout of CMTK call.
    fail_value
                    Value to use for points that failed to transform. By
                    default we use `np.nan`.

    Returns
    -------
    pointsxf :      (N, 3) numpy array
                    The parse transformed points.

    """
    # The original stout is tuple where we care only about the second one
    if isinstance(output, tuple):
        output = output[0]

    pointsx = []
    # Split string into rows - lazily using a generator
    for row in (x.group(1) for x in re.finditer(r"(.*?) \n", output.decode())):
        # Split into values
        values = row.split(" ")

        # If this point failed
        if len(values) != 3:
            values = [fail_value] * 3
        else:
            values = [float(v) for v in values]

        pointsx.append(values)

    return np.asarray(pointsx)

`to_dfield` #

Convert transform to dense deformation field.

Note this always requires CMTK: navis-fastcore implements point but not image transforms.

PARAMETER	DESCRIPTION
`template`	`This defines the bounds and voxel size of the deformation field. Typically, this would correspond to the source space (i.e. the moving image). We can work with: - str: a filepath to a NRRD file - TemplateBrain: a navis TemplateBrain object - a tuple/list/array with (Nx, Ny, Nz, dx, dy, dz) where N is the number of voxels in each dimension and d is the voxel size.` TYPE: `str \| TemplateBrain \| (Nx, Ny, Nz, dx, dy, dz) tuple`
`out`	`NRRD filepath to save the deformation field. If None (default), the deformation field will be returned as numpy array.` TYPE: `str` DEFAULT: `None`
`absolute`	`Whether to return absolute coordinates or offsets.` TYPE: `bool` DEFAULT: `True`
`verbose`	`Whether to print CMTK output.` TYPE: `bool` DEFAULT: `False`

RETURNS	DESCRIPTION
`dfield`	The dense deformation field as (3, Nx, Ny, Nz) numpy array with either absolute coordinates or offsets where the first dimension contains the x/y/z coordinates. Note that the coordinates stored in the field are in physical space (e.g. microns), not voxel space. TYPE: `np.ndarray`
`header`	The NRRD header associated with the deformation field. TYPE: `dict`

`to_grid_transform` #

Convert to GridTransform via dense deformation field.

PARAMETER DESCRIPTION

template

    This defines the bounds and voxel size of the deformation
    field. Typically, this would correspond to the source space
    (i.e. the moving image). We can work with:
      - str: a filepath to a NRRD file
      - TemplateBrain: a navis TemplateBrain object
      - a tuple/list/array with (Nx, Ny, Nz, dx, dy, dz)
        where N is the number of voxels in each dimension
        and d is the voxel size.

TYPE: str | TemplateBrain | (Nx, Ny, Nz, dx, dy, dz) tuple

absolute

    Whether to return absolute coordinates or offsets.

TYPE: bool DEFAULT: True

verbose

    Whether to print CMTK output.

TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
`transform`	A few notes on using the resulting GridTransform: 1. The transform will expect input coordinates in voxel space. See the returned `voxel_size` array! 2. The transformed coordinates, on the other hand, will already be in physical space (e.g. microns). TYPE: `GridTransform`
`voxel_size`	The voxel size of the input space. This is important because the GridTransform will expect input coordinates in voxel space. TYPE: `(3,) numpy array`

Source code in navis/transforms/cmtk.py

def to_grid_transform(self, template, absolute: bool = True, verbose: bool = False):
    """Convert to GridTransform via dense deformation field.

    Parameters
    ----------
    template :  str | TemplateBrain | (Nx, Ny, Nz, dx, dy, dz) tuple
                This defines the bounds and voxel size of the deformation
                field. Typically, this would correspond to the source space
                (i.e. the moving image). We can work with:
                  - str: a filepath to a NRRD file
                  - TemplateBrain: a navis TemplateBrain object
                  - a tuple/list/array with (Nx, Ny, Nz, dx, dy, dz)
                    where N is the number of voxels in each dimension
                    and d is the voxel size.
    absolute :  bool
                Whether to return absolute coordinates or offsets.
    verbose :   bool
                Whether to print CMTK output.

    Returns
    -------
    transform:  GridTransform
                A few notes on using the resulting GridTransform:
                  1. The transform will expect input coordinates in voxel space.
                     See the returned `voxel_size` array!
                  2. The transformed coordinates, on the other hand, will already
                     be in physical space (e.g. microns).
    voxel_size :   (3,) numpy array
                The voxel size of the input space. This is important because
                the GridTransform will expect input coordinates in voxel
                space.
    """
    from .grid import GridTransform

    dfield, header = self.to_dfield(
        template, out=None, absolute=absolute, verbose=verbose
    )

    # Generate the transform
    transform = GridTransform(
        np.transpose(dfield, (1, 2, 3, 0)),
        type="coordinates" if absolute else "offsets",
    )
    # Get voxel size from header - note that the first row is for the 4th dimension
    # (i.e. the vector dimension of the offsets/coordinates) which we don't care about
    space_dir = np.diagonal(header["space directions"][1:])

    return transform, space_dir

`xform` #

Xform data.

PARAMETER DESCRIPTION

points

            Points to xform. DataFrame must have x/y/z columns.

TYPE: (N, 3) numpy array | pandas.DataFrame

affine_only

            Whether to apply only the non-rigid affine
            transform. This is useful if points are outside
            the deformation field and would therefore not
            transform properly.

TYPE: bool DEFAULT: False

affine_fallback

            If True and some points did not transform during the
            non-rigid part of the transformation, we will apply
            only the affine transformation to those points.

TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
`pointsxf`	Transformed points. Points that failed to transform will be `np.nan`. TYPE: `(N, 3) numpy array`

Source code in navis/transforms/cmtk.py

def xform(
    self,
    points: np.ndarray,
    affine_only: bool = False,
    affine_fallback: bool = False,
) -> np.ndarray:
    """Xform data.

    Parameters
    ----------
    points :            (N, 3) numpy array | pandas.DataFrame
                        Points to xform. DataFrame must have x/y/z columns.
    affine_only :       bool
                        Whether to apply only the non-rigid affine
                        transform. This is useful if points are outside
                        the deformation field and would therefore not
                        transform properly.
    affine_fallback :   bool
                        If True and some points did not transform during the
                        non-rigid part of the transformation, we will apply
                        only the affine transformation to those points.

    Returns
    -------
    pointsxf :      (N, 3) numpy array
                    Transformed points. Points that failed to transform will
                    be `np.nan`.

    """
    self.check_if_possible(on_error="raise")

    if self.backend == "fastcore":
        # One parsed object serves every direction - fastcore takes the
        # per-hop direction at xform time, and our `directions` map straight
        # onto its `invert` flags.
        #
        # `fallback_to_affine=True` means "re-run the whole chain affine-only
        # from the original point" - what streamxform --affine-only does, and
        # so what we need while the binary backend still exists.
        # See the class Notes: switch this to "hop" once it doesn't.
        reg = get_cmtk_reg(tuple(str(r) for r in self.regs))
        return reg.xform(
            parse_points(points),
            transform="affine" if affine_only else "warp",
            fallback_to_affine=bool(affine_fallback) and not affine_only,
            invert=[d == "inverse" for d in self.directions],
            n_cores=self.threads,
        )

    if isinstance(points, pd.DataFrame):
        # Make sure x/y/z columns are present
        if np.any([c not in points for c in ["x", "y", "z"]]):
            raise ValueError("points DataFrame must have x/y/z columns.")
    elif (
        isinstance(points, np.ndarray) and points.ndim == 2 and points.shape[1] == 3
    ):
        points = pd.DataFrame(points, columns=["x", "y", "z"])
    else:
        raise TypeError(
            "`points` must be numpy array of shape (N, 3) or "
            "pandas DataFrame with x/y/z columns"
        )

    # Generate the result
    args = self.make_args(affine_only=affine_only)
    proc = subprocess.Popen(args, stdin=subprocess.PIPE, stdout=subprocess.PIPE)

    # Pipe in the points
    points_str = points[["x", "y", "z"]].to_string(index=False, header=False)

    # Do not use proc.stdin.write to avoid output buffer becoming full
    # before we finish piping in stdin.
    # proc.stdin.write(points_str.encode())
    # output = proc.communicate()

    # Read out results
    # This is equivalent to e.g.:
    # $ streamxform -args <<< "10, 10, 10"
    output = proc.communicate(input=points_str.encode())

    # If no output, something went wrong
    if not output[0]:
        raise utils.CMTKError("CMTK produced no output. Check points?")

    # Xformed points
    xf = self.parse_cmtk_output(output, fail_value=np.nan)

    # Check if any points not xformed
    if affine_fallback and not affine_only:
        not_xf = np.any(np.isnan(xf), axis=1)
        if np.any(not_xf):
            xf[not_xf] = self.xform(points.loc[not_xf], affine_only=True)

    return xf

`xform_image` #

Transform an image using CMTK's reformatx.

Note this always requires CMTK: navis-fastcore implements point but not image transforms.

PARAMETER	DESCRIPTION
`im`	`The floating image to transform.` TYPE: `3D numpy array \| filepath`
`target`	`Defines the target image: dimensions in voxels (N), the voxel size (d) and optionally an origin (0) for the target image. Can be provided as a string (name of a template), a TemplateBrain object, a tuple/list/array with the target specs.` TYPE: `str \| TemplateBrain \| (Nx, Ny, Nz, dx, dy, dz) \| (Nx, Ny, Nz, dx, dy, dz, Ox, Oy, Oz)`
`out`	`The filepath to save the transformed image. If None (default), will return the transformed image as np.ndarray.` TYPE: `str` DEFAULT: `None`
`interpolation`	`The interpolation method to use.` TYPE: `linear \| nn \| cubic \| pv \| sinc - cosine \| sinc - hamming` DEFAULT: `'linear'`
`verbose`	`Whether to print CMTK output.` TYPE: `bool` DEFAULT: `False`

RETURNS	DESCRIPTION
`np.ndarray \| None`	If out is None, returns the transformed image as np.ndarray. Otherwise, None.

Source code in navis/transforms/cmtk.py

@requires_cmtk
def xform_image(
    self,
    im,
    target,
    out=None,
    interpolation="linear",
    verbose=False,
):
    """Transform an image using CMTK's reformatx.

    Note this always requires CMTK: navis-fastcore implements point but not
    image transforms.

    Parameters
    ----------
    im :        3D numpy array | filepath
                The floating image to transform.
    target :    str | TemplateBrain | (Nx, Ny, Nz, dx, dy, dz) | (Nx, Ny, Nz, dx, dy, dz, Ox, Oy, Oz)
                Defines the target image: dimensions in voxels (N), the voxel size (d) and optionally
                an origin (0) for the target image. Can be provided as a string (name of a template),
                a TemplateBrain object, a tuple/list/array with the target specs.
    out :       str, optional
                The filepath to save the transformed image. If None (default), will return the
                transformed image as np.ndarray.
    interpolation : "linear" | "nn" | "cubic" | "pv" | "sinc-cosine" | "sinc-hamming"
                The interpolation method to use.
    verbose :   bool
                Whether to print CMTK output.

    Returns
    -------
    np.ndarray | None
                If out is None, returns the transformed image as np.ndarray. Otherwise, None.

    """
    assert interpolation in (
        "linear",
        "nn",
        "cubic",
        "pv",
        "sinc-cosine",
        "sinc-hamming",
    )

    # `reformatx` expects this format:
    # ./reformatx --floating {INPUT_FILE} -o {OUTPUT_FILE} {REFERENCE_SPECS} {TRANSFORMS}
    # where:
    # - {INPUT_FILE} is the image to transform
    # - {OUTPUT_FILE} is where the output will be saved
    # - {REFERENCE_SPECS} defines the target space; this needs to be either a NRRD
    #   file from which CMTK can extract the target grid or the actual specs:
    #   "--target-grid Nx,Ny,Nz:dX,dY,dZ:[Ox,Oy,Oz]" where N is the number of
    #   voxels in each dimension and d is the voxel size. The optional O is the
    #   origin of the image. If not provided, it is assumed to be (0, 0, 0).
    # - {TRANSFORMS} are the CMTK transform(s) to apply; prefix with "--inverse" to invert
    # Below command works to convert JFRC2 to FCWB:
    # /opt/local/lib/cmtk/bin/reformatx --verbose --floating JFRC2.nrrd -o JFRC2_xf.nrrd FCWB.nrrd ~/flybrain-data/BridgingRegistrations/FCWB_JFRC2.list
    # This took 2min 28s - should check if that is actually faster than the look-up approach we use in `images.py`
    # Note that inversion of transforms always comes with an overhead! When using the inverse transform instead, the above command took 2h 40min!

    if verbose and any((d == "inverse " for d in self.directions)):
        config.logger.warning(
            "Using inverse CMTK transforms with reformatx can be very slow! If possible, consider using only forward transforms."
        )

    target_specs = parse_target_specs(target)

    to_remove = []
    if isinstance(im, (str, pathlib.Path)):
        floating = pathlib.Path(im)
        if not im.is_file():
            raise ValueError(f"Image file not found: {im}")
    elif isinstance(im, np.ndarray):
        assert im.ndim == 3
        # Save to temporary file
        with tempfile.NamedTemporaryFile(suffix=".nrrd", delete=False) as tf:
            nrrd.write(tf.name, im)
            floating = tf.name
            to_remove.append(tf.name)
    else:
        raise ValueError(f"Invalid image type: {type(im)}")

    if out is None:
        outfile = tempfile.NamedTemporaryFile(suffix=".nrrd", delete=False).name
        to_remove.append(outfile)
    elif isinstance(out, (str, pathlib.Path)):
        outfile = pathlib.Path(out).resolve()
    else:
        raise ValueError(f"Invalid output type: {type(out)}")

    # Compile the command
    args = [str(_cmtkbin / "reformatx")]
    args += [f"-o {outfile}"]
    args += [f"--floating {floating}"]
    args += [f"--{interpolation}"]
    args += [target_specs]

    # Add the regargs
    args += self.regargs

    try:
        # run the binary
        # avoid resourcewarnings with null
        with open(os.devnull, "w") as devnull:
            startupinfo = None
            if platform.system() == "Windows":
                startupinfo = subprocess.STARTUPINFO()
                startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW
            if verbose:
                # in debug mode print the output
                stdout = None
            else:
                stdout = devnull

            if verbose:
                config.logger.info("executing: {}".format(" ".join(args)))
            check_call(
                args,
                stdout=stdout,
                stderr=subprocess.STDOUT,
                startupinfo=startupinfo,
            )

        if out is None:
            # Return transformed image
            return nrrd.read(outfile)
        elif verbose:
            config.logger.info(f"Transformed image saved to {outfile}")
    except BaseException:
        raise
    finally:
        # Clean up temporary files
        for f in to_remove:
            os.remove(f)

`navis.transforms.ElastixTransform` #

Elastix transforms of 3D spatial data.

Requires either Elastix to be installed, or navis-fastcore - the latter implements elastix transforms in Rust and needs no external binaries. See the backend parameter.

Based on code by Jasper Phelps (https://github.com/jasper-tms/pytransformix).

Note that elastix transforms can only be inverted on the "fastcore" backend: the transformix binary has no way to compute an inverse.

PARAMETER DESCRIPTION

file

            Filepath to elastix transformation file.

TYPE: str

copy_files

            Any files that need to be copied into the temporary
            directory where we perform the transform. These are
            typically files supplemental to the main transform
            file (e.g. defining an additional affine transform).

            Only relevant for the "binary" backend, and only because
            `transformix` looks for chained transform files in the
            current directory. fastcore resolves them itself - it
            follows the recorded path and, failing that, looks for
            that file's basename next to the transform that named it,
            which is the same thing copying them into one directory
            achieves. On that backend this argument is ignored.

TYPE: filepath | list DEFAULT: None

backend

            Which implementation to use. `None` (default) defers to
            `navis.config.default_transform_backend`. "binary" is
            deprecated and will be removed in 3.0.

TYPE: "auto" | "fastcore" | "binary" DEFAULT: None

Examples:

>>> from navis import transforms
>>> tr = transforms.ElastixTransform('/path/to/transform/transform')
>>> tr.xform(points)

Notes

Inverting an elastix transform is a Gauss-Newton solve per point: ~80x the cost of going forwards (and ~4x a dedicated reverse registration), it drops the ~0.04% of points that have no preimage, and where the warp folds ~4% of points come back somewhere other than they started - inherent to inverting a non-injective map, not a solver defect.

So a purpose-built reverse registration should always beat an on-the-fly inversion. It does: _pick_edge prefers a forward registration over the inverse of its counterpart unconditionally, and inverse_weight_factor (above) keeps an inverse hop from looking like a shortcut when routing.

navis.config.elastix_invertible nonetheless defaults to off, and should stay that way while the deprecated "binary" backend is still selectable: it cannot invert at all, so turning this on would let the two backends find different routes. Flip it in 3.0, when that backend goes away. It is safe to enable - on flybrains it changes nothing at all, because every elastix registration there already ships with a reverse - so the only thing it buys today is a route where somebody registered an elastix transform and no reverse for it.

One thing worth knowing before anyone tunes this: a dedicated reverse registration is not the inverse of its forward twin - they are independent fits. Ours agree to a median of ~0.2 world units, but the dedicated reverse only round-trips through the forward map to a median of 0.17 (p95 128), where the numerical inverse round-trips to ~1e-13. Preferring the dedicated registration is a choice about which map you want - the one its authors fitted, and the one the field uses - not a claim that it is more accurate.

Source code in navis/transforms/elastix.py

class ElastixTransform(BackendMixin, BaseTransform):
    """Elastix transforms of 3D spatial data.

    Requires either [Elastix](https://github.com/SuperElastix/elastix/) to be
    installed, or navis-fastcore - the latter implements elastix transforms in
    Rust and needs no external binaries. See the `backend` parameter.

    Based on code by Jasper Phelps (<https://github.com/jasper-tms/pytransformix>).

    Note that elastix transforms can only be inverted on the "fastcore" backend:
    the `transformix` binary has no way to compute an inverse.

    Parameters
    ----------
    file :              str
                        Filepath to elastix transformation file.
    copy_files :        filepath | list, optional
                        Any files that need to be copied into the temporary
                        directory where we perform the transform. These are
                        typically files supplemental to the main transform
                        file (e.g. defining an additional affine transform).

                        Only relevant for the "binary" backend, and only because
                        `transformix` looks for chained transform files in the
                        current directory. fastcore resolves them itself - it
                        follows the recorded path and, failing that, looks for
                        that file's basename next to the transform that named it,
                        which is the same thing copying them into one directory
                        achieves. On that backend this argument is ignored.
    backend :           "auto" | "fastcore" | "binary", optional
                        Which implementation to use. `None` (default) defers to
                        `navis.config.default_transform_backend`. "binary" is
                        deprecated and will be removed in 3.0.

    Examples
    --------
    >>> from navis import transforms
    >>> tr = transforms.ElastixTransform('/path/to/transform/transform')
    >>> tr.xform(points) # doctest: +SKIP

    Notes
    -----
    Inverting an elastix transform is a Gauss-Newton solve per point: ~80x the cost
    of going forwards (and ~4x a dedicated reverse registration), it drops the ~0.04%
    of points that have no preimage, and where the warp folds ~4% of points come back
    somewhere other than they started - inherent to inverting a non-injective map,
    not a solver defect.

    So a purpose-built reverse registration should always beat an on-the-fly
    inversion. It does: `_pick_edge` prefers a forward registration over the inverse
    of its counterpart unconditionally, and `inverse_weight_factor` (above) keeps an
    inverse hop from looking like a shortcut when routing.

    `navis.config.elastix_invertible` nonetheless defaults to *off*, and should stay
    that way while the deprecated "binary" backend is still selectable: it cannot
    invert at all, so turning this on would let the two backends find different
    routes. Flip it in 3.0, when that backend goes away. It is safe to enable - on
    `flybrains` it changes nothing at all, because every elastix registration there
    already ships with a reverse - so the only thing it buys today is a route where
    somebody registered an elastix transform and no reverse for it.

    One thing worth knowing before anyone tunes this: a dedicated reverse registration
    is *not* the inverse of its forward twin - they are independent fits. Ours agree
    to a median of ~0.2 world units, but the dedicated reverse only round-trips through
    the forward map to a median of 0.17 (p95 128), where the numerical inverse
    round-trips to ~1e-13. Preferring the dedicated registration is a choice about
    *which map you want* - the one its authors fitted, and the one the field uses - not
    a claim that it is more accurate.

    """

    # Traversing an elastix transform backwards is a Gauss-Newton solve per point:
    # ~80x the cost of going forwards, and it is *lossy* - points whose preimage
    # does not exist come back as NaN (~0.04%), and where the warp folds another
    # ~4% land somewhere other than they started.
    #
    # 5 means navis will detour through up to four extra forward transforms rather
    # than invert an elastix registration once. Deliberately not the 80x runtime
    # ratio: routing trades accuracy, not seconds, and every extra hop carries its
    # own interpolation error and loses its own out-of-domain points. Five is a
    # judgement call - it says "this is a last resort", without licensing an absurd
    # detour to avoid it.
    inverse_weight_factor = 5

    _invert = False

    def __init__(self, file: str, copy_files=None, backend: str = None):
        self.file = pathlib.Path(file)
        self.copy_files = copy_files if copy_files is not None else []
        self._backend = backend

    def __eq__(self, other: "ElastixTransform") -> bool:
        """Implement equality comparison."""
        if isinstance(other, ElastixTransform):
            if self.file == other.file and self._invert == other._invert:
                return True
        return False

    def __neg__(self) -> "ElastixTransform":
        """Invert transform."""
        # Note this is defined unconditionally (so that `hasattr(tr, '__neg__')`
        # is True) but the honest answer lives in `.can_invert`/`.invertible` -
        # which is why the registry must ask *those* instead.
        if self.backend != "fastcore":
            raise NotImplementedError(
                "Inverting elastix transforms requires navis-fastcore "
                "(`pip install navis-fastcore`): the `transformix` binary "
                "cannot compute an inverse."
            )

        if not self.can_invert:
            raise NotImplementedError(
                f"Elastix transform {self.file} cannot be inverted: its chain "
                'combines transforms via "Add".'
            )

        x = self.copy()
        x._invert = not x._invert
        return x

    @property
    def can_invert(self) -> bool:
        """Whether `-transform` works.

        Needs the fastcore backend (`transformix` cannot invert) *and* a file
        that is actually invertible - a chain combined via "Add" is not. The
        latter is settled with fastcore's header-only probe, so this stays cheap
        enough to ask of every transform in the registry.

        """
        if self.backend != "fastcore" or not self.file.is_file():
            return False
        return elastix_is_invertible(str(self.file))

    @property
    def invertible(self) -> bool:
        """Whether the bridging graph may traverse this transform backwards.

        Deliberately *narrower* than `can_invert`: inverting works whenever we're
        on the fastcore backend with an invertible file, but we only let the
        registry route through an inverse if `navis.config.elastix_invertible`
        says so. Every elastix registration we know of ships with a purpose-built
        reverse registration, so the inverse buys no new connectivity - it just
        adds a cheaper-looking parallel edge that drags unrelated routes through
        it.

        Note the config check comes first, so that in the default case this costs
        a single lookup and never touches the disk - it is evaluated for every
        registered transform each time the bridging graph is built.

        """
        if not getattr(config, "elastix_invertible", False):
            return False
        return self.can_invert

    def check_if_possible(self, on_error: str = "raise"):
        """Check if this transform is possible."""
        msg = None
        if not self.file.is_file():
            msg = f"Transformation file {self.file} not found."
        elif self.backend == "binary":
            if not _elastixbin:
                msg = (
                    "Folder with elastix binaries not found. Make sure the "
                    "directory is in your PATH environment variable - or install "
                    "navis-fastcore to transform points without elastix."
                )
            elif self._invert:
                msg = (
                    "Inverse elastix transforms require navis-fastcore: the "
                    "`transformix` binary cannot compute an inverse."
                )
        elif self._invert and not elastix_is_invertible(str(self.file)):
            # Header-only probe - no need to parse the coefficients just to find
            # out we can't use them.
            msg = (
                f"Elastix transform {self.file} cannot be inverted: its "
                'chain combines transforms via "Add".'
            )

        if msg and on_error == "raise":
            raise BaseException(msg)

        return msg

    def copy(self) -> "ElastixTransform":
        """Return copy."""
        # Attributes not to copy
        no_copy = []
        # Generate new empty transform
        x = self.__class__(self.file)
        # Override with this neuron's data
        x.__dict__.update(
            {k: copy.copy(v) for k, v in self.__dict__.items() if k not in no_copy}
        )

        return x

    def write_input_file(self, points, filepath):
        """Write a numpy array in format required by transformix."""
        with open(filepath, "w") as f:
            f.write("point\n{}\n".format(len(points)))
            for x, y, z in points:
                f.write(f"{x:f} {y:f} {z:f}\n")

    def read_output_file(self, filepath) -> np.ndarray:
        """Load output file.

        Parameter
        ---------
        filepath :      str
                        Filepath to output file.

        Returns
        -------
        pointsxf :      (N, 3) numpy array
                        The parse transformed points.

        """
        points = []
        with open(filepath, "r") as f:
            for line in f.readlines():
                output = line.split("OutputPoint = [ ")[1].split(" ]")[0]
                points.append([float(i) for i in output.split(" ")])
        return np.array(points)

    def xform(self, points: np.ndarray, return_logs=False) -> np.ndarray:
        """Xform data.

        Parameters
        ----------
        points :        (N, 3) numpy array | pandas.DataFrame
                        Points to xform. DataFrame must have x/y/z columns.
        return_logs :   bool
                        If True, will return logs instead of transformed points.
                        Really only useful for debugging.

        Returns
        -------
        pointsxf :      (N, 3) numpy array
                        Transformed points.

        """
        self.check_if_possible(on_error="raise")

        if self.backend == "fastcore":
            if return_logs:
                raise ValueError(
                    "`return_logs=True` requires the `transformix` binary "
                    "(use `backend='binary'`)."
                )

            if self.copy_files:
                logger.debug(
                    "`copy_files` is ignored on the fastcore backend: it resolves "
                    "chained transform files itself."
                )

            # One parsed object serves both directions - fastcore takes `invert`
            # at xform time - so a transform and its inverse share a cache entry.
            tr = get_elastix_transform(str(self.file))
            return tr.xform(parse_points(points), invert=self._invert)

        if isinstance(points, pd.DataFrame):
            # Make sure x/y/z columns are present
            if np.any([c not in points for c in ["x", "y", "z"]]):
                raise ValueError("points DataFrame must have x/y/z columns.")
            points = points[["x", "y", "z"]].values
        elif not (
            isinstance(points, np.ndarray) and points.ndim == 2 and points.shape[1] == 3
        ):
            raise TypeError(
                "`points` must be numpy array of shape (N, 3) or "
                "pandas DataFrame with x/y/z columns"
            )

        # Everything happens in a temporary directory
        with tempfile.TemporaryDirectory() as tempdir:
            p = pathlib.Path(tempdir)

            # If required, copy additional files into the temporary directory
            if self.copy_files:
                for f in make_iterable(self.copy_files):
                    _ = pathlib.Path(shutil.copy(f, p))

            # Write points to file
            in_file = p / "inputpoints.txt"
            self.write_input_file(points, in_file)

            out_file = p / "outputpoints.txt"

            # Prepare the command
            command = [
                _elastixbin / "transformix",
                "-out",
                str(p),
                "-tp",
                str(self.file),
                "-def",
                str(in_file),
            ]

            # Keep track of current working directory
            cwd = os.getcwd()
            try:
                # Change working directory to the temporary directory
                # This is apparently required because elastix stupidly expects
                # any secondary transform files to be in the current directory
                # (as opposed to where the main transform is)
                os.chdir(p)
                # Run the actual transform
                proc = subprocess.run(command, stdout=subprocess.PIPE)
            except BaseException:
                raise
            finally:
                # This makes sure we land on our feet even in case of an error
                os.chdir(cwd)

            if return_logs:
                logfile = p / "transformix.log"
                if not logfile.is_file():
                    raise FileNotFoundError("No log file found.")
                with open(logfile) as f:
                    logs = f.read()
                return logs

            if not out_file.is_file():
                raise FileNotFoundError(
                    "Elastix transform did not produce any "
                    f"output:\n {proc.stdout.decode()}"
                )

            points_xf = self.read_output_file(out_file)

        return points_xf

`can_invert: bool` `property` #

Whether -transform works.

Needs the fastcore backend (transformix cannot invert) and a file that is actually invertible - a chain combined via "Add" is not. The latter is settled with fastcore's header-only probe, so this stays cheap enough to ask of every transform in the registry.

`invertible: bool` `property` #

Whether the bridging graph may traverse this transform backwards.

Deliberately narrower than can_invert: inverting works whenever we're on the fastcore backend with an invertible file, but we only let the registry route through an inverse if navis.config.elastix_invertible says so. Every elastix registration we know of ships with a purpose-built reverse registration, so the inverse buys no new connectivity - it just adds a cheaper-looking parallel edge that drags unrelated routes through it.

Note the config check comes first, so that in the default case this costs a single lookup and never touches the disk - it is evaluated for every registered transform each time the bridging graph is built.

`check_if_possible` #

Check if this transform is possible.

Source code in navis/transforms/elastix.py

def check_if_possible(self, on_error: str = "raise"):
    """Check if this transform is possible."""
    msg = None
    if not self.file.is_file():
        msg = f"Transformation file {self.file} not found."
    elif self.backend == "binary":
        if not _elastixbin:
            msg = (
                "Folder with elastix binaries not found. Make sure the "
                "directory is in your PATH environment variable - or install "
                "navis-fastcore to transform points without elastix."
            )
        elif self._invert:
            msg = (
                "Inverse elastix transforms require navis-fastcore: the "
                "`transformix` binary cannot compute an inverse."
            )
    elif self._invert and not elastix_is_invertible(str(self.file)):
        # Header-only probe - no need to parse the coefficients just to find
        # out we can't use them.
        msg = (
            f"Elastix transform {self.file} cannot be inverted: its "
            'chain combines transforms via "Add".'
        )

    if msg and on_error == "raise":
        raise BaseException(msg)

    return msg

`copy` #

Return copy.

Source code in navis/transforms/elastix.py

def copy(self) -> "ElastixTransform":
    """Return copy."""
    # Attributes not to copy
    no_copy = []
    # Generate new empty transform
    x = self.__class__(self.file)
    # Override with this neuron's data
    x.__dict__.update(
        {k: copy.copy(v) for k, v in self.__dict__.items() if k not in no_copy}
    )

    return x

`read_output_file` #

Load output file.

Parameter

filepath : str Filepath to output file.

RETURNS	DESCRIPTION
`pointsxf`	The parse transformed points. TYPE: `(N, 3) numpy array`

Source code in navis/transforms/elastix.py

def read_output_file(self, filepath) -> np.ndarray:
    """Load output file.

    Parameter
    ---------
    filepath :      str
                    Filepath to output file.

    Returns
    -------
    pointsxf :      (N, 3) numpy array
                    The parse transformed points.

    """
    points = []
    with open(filepath, "r") as f:
        for line in f.readlines():
            output = line.split("OutputPoint = [ ")[1].split(" ]")[0]
            points.append([float(i) for i in output.split(" ")])
    return np.array(points)

`write_input_file` #

Write a numpy array in format required by transformix.

Source code in navis/transforms/elastix.py

def write_input_file(self, points, filepath):
    """Write a numpy array in format required by transformix."""
    with open(filepath, "w") as f:
        f.write("point\n{}\n".format(len(points)))
        for x, y, z in points:
            f.write(f"{x:f} {y:f} {z:f}\n")

`xform` #

Xform data.

PARAMETER	DESCRIPTION
`points`	`Points to xform. DataFrame must have x/y/z columns.` TYPE: `(N, 3) numpy array \| pandas.DataFrame`
`return_logs`	`If True, will return logs instead of transformed points. Really only useful for debugging.` TYPE: `bool` DEFAULT: `False`

RETURNS	DESCRIPTION
`pointsxf`	Transformed points. TYPE: `(N, 3) numpy array`

Source code in navis/transforms/elastix.py

def xform(self, points: np.ndarray, return_logs=False) -> np.ndarray:
    """Xform data.

    Parameters
    ----------
    points :        (N, 3) numpy array | pandas.DataFrame
                    Points to xform. DataFrame must have x/y/z columns.
    return_logs :   bool
                    If True, will return logs instead of transformed points.
                    Really only useful for debugging.

    Returns
    -------
    pointsxf :      (N, 3) numpy array
                    Transformed points.

    """
    self.check_if_possible(on_error="raise")

    if self.backend == "fastcore":
        if return_logs:
            raise ValueError(
                "`return_logs=True` requires the `transformix` binary "
                "(use `backend='binary'`)."
            )

        if self.copy_files:
            logger.debug(
                "`copy_files` is ignored on the fastcore backend: it resolves "
                "chained transform files itself."
            )

        # One parsed object serves both directions - fastcore takes `invert`
        # at xform time - so a transform and its inverse share a cache entry.
        tr = get_elastix_transform(str(self.file))
        return tr.xform(parse_points(points), invert=self._invert)

    if isinstance(points, pd.DataFrame):
        # Make sure x/y/z columns are present
        if np.any([c not in points for c in ["x", "y", "z"]]):
            raise ValueError("points DataFrame must have x/y/z columns.")
        points = points[["x", "y", "z"]].values
    elif not (
        isinstance(points, np.ndarray) and points.ndim == 2 and points.shape[1] == 3
    ):
        raise TypeError(
            "`points` must be numpy array of shape (N, 3) or "
            "pandas DataFrame with x/y/z columns"
        )

    # Everything happens in a temporary directory
    with tempfile.TemporaryDirectory() as tempdir:
        p = pathlib.Path(tempdir)

        # If required, copy additional files into the temporary directory
        if self.copy_files:
            for f in make_iterable(self.copy_files):
                _ = pathlib.Path(shutil.copy(f, p))

        # Write points to file
        in_file = p / "inputpoints.txt"
        self.write_input_file(points, in_file)

        out_file = p / "outputpoints.txt"

        # Prepare the command
        command = [
            _elastixbin / "transformix",
            "-out",
            str(p),
            "-tp",
            str(self.file),
            "-def",
            str(in_file),
        ]

        # Keep track of current working directory
        cwd = os.getcwd()
        try:
            # Change working directory to the temporary directory
            # This is apparently required because elastix stupidly expects
            # any secondary transform files to be in the current directory
            # (as opposed to where the main transform is)
            os.chdir(p)
            # Run the actual transform
            proc = subprocess.run(command, stdout=subprocess.PIPE)
        except BaseException:
            raise
        finally:
            # This makes sure we land on our feet even in case of an error
            os.chdir(cwd)

        if return_logs:
            logfile = p / "transformix.log"
            if not logfile.is_file():
                raise FileNotFoundError("No log file found.")
            with open(logfile) as f:
                logs = f.read()
            return logs

        if not out_file.is_file():
            raise FileNotFoundError(
                "Elastix transform did not produce any "
                f"output:\n {proc.stdout.decode()}"
            )

        points_xf = self.read_output_file(out_file)

    return points_xf

`navis.transforms.FunctionTransform` #

Apply custom function as transform.

PARAMETER	DESCRIPTION
`func`	`Function that accepts and returns an (N, 3) array.` TYPE: `callable`

Source code in navis/transforms/base.py

class FunctionTransform(BaseTransform):
    """Apply custom function as transform.

    Parameters
    ----------
    func :      callable
                Function that accepts and returns an (N, 3) array.

    """

    def __init__(self, func):
        """Initialize."""
        if not callable(func):
            raise TypeError("`func` must be callable")
        self.func = func

    def __eq__(self, other):
        """Check if the same."""
        if not isinstance(other, FunctionTransform):
            return False
        if self.func != other.func:
            return False
        return True

    def copy(self):
        """Return copy."""
        x = self.__class__(self.func)
        x.__dict__.update(self.__dict__)
        return x

    def xform(self, points: np.ndarray) -> np.ndarray:
        """Xform data.

        Parameters
        ----------
        points :        (N, 3) numpy array | pandas.DataFrame
                        Points to xform. DataFrame must have x/y/z columns.

        Returns
        -------
        pointsxf :      (N, 3) numpy array
                        Transformed points.

        """
        if isinstance(points, pd.DataFrame):
            # Make sure x/y/z columns are present
            if np.any([c not in points for c in ["x", "y", "z"]]):
                raise ValueError("points DataFrame must have x/y/z columns.")
            points = points[["x", "y", "z"]].values
        elif not (
            isinstance(points, np.ndarray) and points.ndim == 2 and points.shape[1] == 3
        ):
            raise TypeError(
                "`points` must be numpy array of shape (N, 3) or "
                "pandas DataFrame with x/y/z columns"
            )
        return self.func(points.copy())

`init` #

Initialize.

Source code in navis/transforms/base.py

def __init__(self, func):
    """Initialize."""
    if not callable(func):
        raise TypeError("`func` must be callable")
    self.func = func

`copy` #

Return copy.

Source code in navis/transforms/base.py

def copy(self):
    """Return copy."""
    x = self.__class__(self.func)
    x.__dict__.update(self.__dict__)
    return x

`xform` #

Xform data.

PARAMETER	DESCRIPTION
`points`	`Points to xform. DataFrame must have x/y/z columns.` TYPE: `(N, 3) numpy array \| pandas.DataFrame`

RETURNS	DESCRIPTION
`pointsxf`	Transformed points. TYPE: `(N, 3) numpy array`

Source code in navis/transforms/base.py

def xform(self, points: np.ndarray) -> np.ndarray:
    """Xform data.

    Parameters
    ----------
    points :        (N, 3) numpy array | pandas.DataFrame
                    Points to xform. DataFrame must have x/y/z columns.

    Returns
    -------
    pointsxf :      (N, 3) numpy array
                    Transformed points.

    """
    if isinstance(points, pd.DataFrame):
        # Make sure x/y/z columns are present
        if np.any([c not in points for c in ["x", "y", "z"]]):
            raise ValueError("points DataFrame must have x/y/z columns.")
        points = points[["x", "y", "z"]].values
    elif not (
        isinstance(points, np.ndarray) and points.ndim == 2 and points.shape[1] == 3
    ):
        raise TypeError(
            "`points` must be numpy array of shape (N, 3) or "
            "pandas DataFrame with x/y/z columns"
        )
    return self.func(points.copy())

`navis.transforms.GridTransform` #

Deformation or coordinate field transform of 3D spatial data.

This is effectively a simpler version of the H5transform class and only supports deformation fields stored as numpy arrays in memory.

PARAMETER	DESCRIPTION
`field`	`Deformation/coordinate field. The last dimension must contain the x/y/z coordinates.` TYPE: `(Nx, Ny, Nz, 3) numpy array`
`type`	`Whether the field contains absolute coordinates or offsets/displacements. Offsets are added to the input coordinates, while coordinates returned as is.` TYPE: `"offsets" \| "coordinates"` DEFAULT: `'offsets'`
`spacing`	`Spacing of the deformation field in x/y/z. If not provided, spacing of 1 is assumed.` TYPE: `tuple \| list \| numpy array` DEFAULT: `None`
`offset`	`Offset of the deformation field in x/y/z. If not provided, offset of 0 is assumed. N.B. offsets are applied _after_ spacing, i.e. need to be provided in voxel space of the deformation field.` TYPE: `tuple \| list \| numpy array` DEFAULT: `None`

Source code in navis/transforms/grid.py

class GridTransform(BaseTransform):
    """Deformation or coordinate field transform of 3D spatial data.

    This is effectively a simpler version of the H5transform class and
    only supports deformation fields stored as numpy arrays in memory.

    Parameters
    ----------
    field :         (Nx, Ny, Nz, 3) numpy array
                    Deformation/coordinate field. The last dimension must
                    contain the x/y/z coordinates.
    type :          "offsets" | "coordinates"
                    Whether the field contains absolute coordinates or offsets/displacements.
                    Offsets are added to the input coordinates, while
                    coordinates returned as is.
    spacing :       tuple | list | numpy array, optional
                    Spacing of the deformation field in x/y/z. If not provided,
                    spacing of 1 is assumed.
    offset :        tuple | list | numpy array, optional
                    Offset of the deformation field in x/y/z. If not provided,
                    offset of 0 is assumed. N.B. offsets are applied _after_ spacing,
                    i.e. need to be provided in voxel space of the deformation field.

    """

    def __init__(
        self, field: np.ndarray, type: str = "offsets", spacing=None, offset=None
    ):
        """Init class."""
        assert (
            field.ndim == 4 and field.shape[-1] == 3
        ), "Field must be a 4D numpy array with the last dimension of size 3."
        assert type in (
            "coordinates",
            "offsets",
        ), "type must be 'coordinates' or 'offsets'."
        self.field = field
        self.type = type
        self.dtype = field.dtype
        self.spacing = spacing
        self.offset = offset

    def __eq__(self, other) -> bool:
        """Compare with other Transform."""
        if isinstance(other, GridTransform):
            if np.array_equal(self.field, other.field):
                return True
        return False

    def __neg__(self) -> "GridTransform":
        """Invert direction."""
        # Note to future self: we could implement this by computing the inverse
        # deformation field using fixed-point iteration, but that's
        # non-trivial and not needed right now.
        raise NotImplementedError("Inversion of GridTransform is not implemented.")

    @property
    def spacing(self):
        return self._spacing

    @spacing.setter
    def spacing(self, value):
        if value is None:
            self._spacing = None
        else:
            value = np.asarray(value)
            assert (
                value.ndim == 1 and value.size == 3
            ), "spacing must be a tuple/list/array of size 3."
            self._spacing = np.array(value, dtype=self.dtype)

    @property
    def offset(self):
        return self._offset

    @offset.setter
    def offset(self, value):
        if value is None:
            self._offset = None
        else:
            value = np.asarray(value)
            assert (
                value.ndim == 1 and value.size == 3
            ), "offset must be a tuple/list/array of size 3."
            self._offset = np.array(value, dtype=self.dtype)

    @property
    def affine(self) -> AffineTransform:
        """Return affine part of the transform."""
        # We're delaying the calculation of the affine part until it's needed
        if not hasattr(self, "_affine"):
            self.calculate_affine()
        return self._affine

    @property
    def shape(self) -> tuple:
        """Return shape of the deformation field."""
        return self.field.shape

    def calculate_affine(self) -> None:
        """Calculate affine part of the transform."""
        # The strategy here is this:
        # 1. Take the 8 corners of the deformation field
        # 2. Transform them using the deformation field
        # 3. Treat them as landmarks and compute the affine transform using morphops

        mx = np.array(self.shape) - 1  # max indices in each dimension
        points = np.array(
            [
                [0, 0, 0],
                [0, 0, mx[2]],
                [0, mx[1], 0],
                [0, mx[1], mx[2]],
                [mx[0], 0, 0],
                [mx[0], 0, mx[2]],
                [mx[0], mx[1], 0],
                [mx[0], mx[1], mx[2]],
            ]
        )

        if self.offset is not None:
            points = points - self.offset[np.newaxis, :]

        if self.spacing is not None:
            points = points * self.spacing[np.newaxis, :]

        points_xf = self.xform(points, affine_fallback=False)

        m = TPStransform(points, points_xf).matrix_affine

        # Calculate the affine part as the mean displacement across the field
        self._affine = AffineTransform(m)

    def copy(self, copy_data: bool = False) -> "GridTransform":
        """Return copy."""
        return GridTransform(
            self.field.copy() if copy_data else self.field,
            type=self.type,
            spacing=self.spacing,
            offset=self.offset
        )

    @classmethod
    def from_file(cls, filepath: str) -> "GridTransform":
        """Create GridTransform from a file.

        Parameters
        ----------
        filepath :      str
                        Path to file. Currently supported formats:
                          - NRRD files with deformation fields
                          - Numpy .npy files with deformation fields
                          - Nifti files (.nii, .nii.gz) with deformation fields

        Returns
        -------
        GridTransform instance.

        """
        filepath = Path(filepath)
        if not filepath.exists():
            raise FileNotFoundError(f"File {filepath} does not exist.")

        if filepath.suffix in [".npy"]:
            field = np.load(filepath)
        elif filepath.suffix in [".nii", ".nii.gz"]:
            try:
                import nibabel as nib
            except ModuleNotFoundError:
                raise ImportError(
                    "`nibabel` package is required to read Nifti files (.nii, .nii.gz)"
                )

            img = nib.load(str(filepath))
            field = img.get_fdata()
        elif filepath.suffix in [".nrrd"]:
            import nrrd

            field, _ = nrrd.read(str(filepath))
        else:
            raise ValueError(
                f"Unsupported file format: {filepath.suffix}. "
                "Supported formats are: .npy, .nii, .nii.gz, .nrrd"
            )

        return cls(field)

    @classmethod
    def from_warpfield(cls, warpfield):
        """Create GridTransform from a Warpfield deformation field.

        Parameters
        ----------
        warpfield :     warpfield.WarpMap | str
                        Warpfield WarpMap instance or path to a WarpMap h5 file.

        """
        if isinstance(warpfield, str):
            import h5py

            with h5py.File(warpfield, "r") as h5:
                wm = h5["warp_map"]
                field = wm["warp_field"][:]
                block_size = wm["block_size"][:]
                block_stride = wm["block_stride"][:]
                # mov_shape = wm["moving_shape"][:]
                # ref_shape = wm["ref_shape"][:]
        else:
            field = np.asarray(warpfield.warp_field)
            block_size = np.asarray(warpfield.block_size)
            block_stride = np.asarray(warpfield.block_stride)

        # Reshape the field from (3, X, Y, Z) to (X, Y, Z, 3)
        field = np.moveaxis(field, 0, -1)

        spacing = block_stride
        offset = -block_size / block_stride / 2
        # Note to self regarding the offset: this code is taken straight from warpfield.

        return cls(field, type="offsets", spacing=spacing, offset=offset)

    def xform(
        self,
        points: np.ndarray,
        affine_fallback: bool = False,
    ) -> np.ndarray:
        """Xform data.

        Parameters
        ----------
        points :            (N, 3) numpy array | pandas.DataFrame
                            Points to xform. DataFrame must have x/y/z columns.
        affine_fallback :   bool
                            If False, sample outside the deformation field is treated
                            as constant `np.nan`. If True, sample outside the deformation
                            field uses nearest-neighbor extrapolation, i.e. we're not
                            exactly applying a global affine but we're keeping the
                            parameter name for compatibility with other transforms.

        Returns
        -------
        pointsxf :          (N, 3) numpy.ndarray
                            Transformed points. Will contain `np.nan` for points
                            that did not transform.

        """
        if isinstance(points, pd.DataFrame):
            # Make sure x/y/z columns are present
            if np.any([c not in points for c in ["x", "y", "z"]]):
                raise ValueError("points DataFrame must have x/y/z columns.")
            points = points[["x", "y", "z"]].values

        if (
            not isinstance(points, np.ndarray)
            or points.ndim != 2
            or points.shape[1] != 3
        ):
            raise TypeError(
                "`points` must be numpy array of shape (N, 3) or "
                "pandas DataFrame with x/y/z columns"
            )

        points_vxl = points

        if self.spacing is not None:
            points_vxl = points_vxl / self.spacing[np.newaxis, :]

        if self.offset is not None:
            points_vxl = points_vxl + self.offset[np.newaxis, :]

        # map_coordinates expects coordinates as (ndim, n_points)
        coords = points_vxl.T

        mode = "nearest" if affine_fallback else "constant"
        cval = 0 if affine_fallback else np.nan

        # Prepare output array
        if self.type == "coordinates":
            points_xf = np.zeros(points.shape, dtype=self.dtype)
        else:  # offsets
            points_xf = points.astype(self.dtype, copy=True)

        # Interpolate coordinates, re-combine to an x/y/z array and
        # add to the input points (if coordinates, the points are zeroed out above)
        shifts = np.vstack(
            (
                map_coordinates(self.field[:, :, :, 0], coords, order=1, mode=mode, cval=cval),
                map_coordinates(self.field[:, :, :, 1], coords, order=1, mode=mode, cval=cval),
                map_coordinates(self.field[:, :, :, 2], coords, order=1, mode=mode, cval=cval),
            )
        ).T.astype(self.dtype)
        points_xf += shifts

        return points_xf

`affine: AffineTransform` `property` #

Return affine part of the transform.

`shape: tuple` `property` #

Return shape of the deformation field.

`init` #

Init class.

Source code in navis/transforms/grid.py

def __init__(
    self, field: np.ndarray, type: str = "offsets", spacing=None, offset=None
):
    """Init class."""
    assert (
        field.ndim == 4 and field.shape[-1] == 3
    ), "Field must be a 4D numpy array with the last dimension of size 3."
    assert type in (
        "coordinates",
        "offsets",
    ), "type must be 'coordinates' or 'offsets'."
    self.field = field
    self.type = type
    self.dtype = field.dtype
    self.spacing = spacing
    self.offset = offset

`calculate_affine` #

Calculate affine part of the transform.

Source code in navis/transforms/grid.py

def calculate_affine(self) -> None:
    """Calculate affine part of the transform."""
    # The strategy here is this:
    # 1. Take the 8 corners of the deformation field
    # 2. Transform them using the deformation field
    # 3. Treat them as landmarks and compute the affine transform using morphops

    mx = np.array(self.shape) - 1  # max indices in each dimension
    points = np.array(
        [
            [0, 0, 0],
            [0, 0, mx[2]],
            [0, mx[1], 0],
            [0, mx[1], mx[2]],
            [mx[0], 0, 0],
            [mx[0], 0, mx[2]],
            [mx[0], mx[1], 0],
            [mx[0], mx[1], mx[2]],
        ]
    )

    if self.offset is not None:
        points = points - self.offset[np.newaxis, :]

    if self.spacing is not None:
        points = points * self.spacing[np.newaxis, :]

    points_xf = self.xform(points, affine_fallback=False)

    m = TPStransform(points, points_xf).matrix_affine

    # Calculate the affine part as the mean displacement across the field
    self._affine = AffineTransform(m)

`copy` #

Return copy.

Source code in navis/transforms/grid.py

def copy(self, copy_data: bool = False) -> "GridTransform":
    """Return copy."""
    return GridTransform(
        self.field.copy() if copy_data else self.field,
        type=self.type,
        spacing=self.spacing,
        offset=self.offset
    )

`from_file` `classmethod` #

Create GridTransform from a file.

PARAMETER	DESCRIPTION
`filepath`	`Path to file. Currently supported formats: - NRRD files with deformation fields - Numpy .npy files with deformation fields - Nifti files (.nii, .nii.gz) with deformation fields` TYPE: `str`

RETURNS	DESCRIPTION
`GridTransform instance.`

Source code in navis/transforms/grid.py

@classmethod
def from_file(cls, filepath: str) -> "GridTransform":
    """Create GridTransform from a file.

    Parameters
    ----------
    filepath :      str
                    Path to file. Currently supported formats:
                      - NRRD files with deformation fields
                      - Numpy .npy files with deformation fields
                      - Nifti files (.nii, .nii.gz) with deformation fields

    Returns
    -------
    GridTransform instance.

    """
    filepath = Path(filepath)
    if not filepath.exists():
        raise FileNotFoundError(f"File {filepath} does not exist.")

    if filepath.suffix in [".npy"]:
        field = np.load(filepath)
    elif filepath.suffix in [".nii", ".nii.gz"]:
        try:
            import nibabel as nib
        except ModuleNotFoundError:
            raise ImportError(
                "`nibabel` package is required to read Nifti files (.nii, .nii.gz)"
            )

        img = nib.load(str(filepath))
        field = img.get_fdata()
    elif filepath.suffix in [".nrrd"]:
        import nrrd

        field, _ = nrrd.read(str(filepath))
    else:
        raise ValueError(
            f"Unsupported file format: {filepath.suffix}. "
            "Supported formats are: .npy, .nii, .nii.gz, .nrrd"
        )

    return cls(field)

`from_warpfield` `classmethod` #

Create GridTransform from a Warpfield deformation field.

PARAMETER	DESCRIPTION
`warpfield`	`Warpfield WarpMap instance or path to a WarpMap h5 file.` TYPE: `warpfield.WarpMap \| str`

Source code in navis/transforms/grid.py

@classmethod
def from_warpfield(cls, warpfield):
    """Create GridTransform from a Warpfield deformation field.

    Parameters
    ----------
    warpfield :     warpfield.WarpMap | str
                    Warpfield WarpMap instance or path to a WarpMap h5 file.

    """
    if isinstance(warpfield, str):
        import h5py

        with h5py.File(warpfield, "r") as h5:
            wm = h5["warp_map"]
            field = wm["warp_field"][:]
            block_size = wm["block_size"][:]
            block_stride = wm["block_stride"][:]
            # mov_shape = wm["moving_shape"][:]
            # ref_shape = wm["ref_shape"][:]
    else:
        field = np.asarray(warpfield.warp_field)
        block_size = np.asarray(warpfield.block_size)
        block_stride = np.asarray(warpfield.block_stride)

    # Reshape the field from (3, X, Y, Z) to (X, Y, Z, 3)
    field = np.moveaxis(field, 0, -1)

    spacing = block_stride
    offset = -block_size / block_stride / 2
    # Note to self regarding the offset: this code is taken straight from warpfield.

    return cls(field, type="offsets", spacing=spacing, offset=offset)

`xform` #

Xform data.

PARAMETER DESCRIPTION

points

            Points to xform. DataFrame must have x/y/z columns.

TYPE: (N, 3) numpy array | pandas.DataFrame

affine_fallback

            If False, sample outside the deformation field is treated
            as constant `np.nan`. If True, sample outside the deformation
            field uses nearest-neighbor extrapolation, i.e. we're not
            exactly applying a global affine but we're keeping the
            parameter name for compatibility with other transforms.

TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
`pointsxf`	Transformed points. Will contain `np.nan` for points that did not transform. TYPE: `(N, 3) numpy.ndarray`

Source code in navis/transforms/grid.py

def xform(
    self,
    points: np.ndarray,
    affine_fallback: bool = False,
) -> np.ndarray:
    """Xform data.

    Parameters
    ----------
    points :            (N, 3) numpy array | pandas.DataFrame
                        Points to xform. DataFrame must have x/y/z columns.
    affine_fallback :   bool
                        If False, sample outside the deformation field is treated
                        as constant `np.nan`. If True, sample outside the deformation
                        field uses nearest-neighbor extrapolation, i.e. we're not
                        exactly applying a global affine but we're keeping the
                        parameter name for compatibility with other transforms.

    Returns
    -------
    pointsxf :          (N, 3) numpy.ndarray
                        Transformed points. Will contain `np.nan` for points
                        that did not transform.

    """
    if isinstance(points, pd.DataFrame):
        # Make sure x/y/z columns are present
        if np.any([c not in points for c in ["x", "y", "z"]]):
            raise ValueError("points DataFrame must have x/y/z columns.")
        points = points[["x", "y", "z"]].values

    if (
        not isinstance(points, np.ndarray)
        or points.ndim != 2
        or points.shape[1] != 3
    ):
        raise TypeError(
            "`points` must be numpy array of shape (N, 3) or "
            "pandas DataFrame with x/y/z columns"
        )

    points_vxl = points

    if self.spacing is not None:
        points_vxl = points_vxl / self.spacing[np.newaxis, :]

    if self.offset is not None:
        points_vxl = points_vxl + self.offset[np.newaxis, :]

    # map_coordinates expects coordinates as (ndim, n_points)
    coords = points_vxl.T

    mode = "nearest" if affine_fallback else "constant"
    cval = 0 if affine_fallback else np.nan

    # Prepare output array
    if self.type == "coordinates":
        points_xf = np.zeros(points.shape, dtype=self.dtype)
    else:  # offsets
        points_xf = points.astype(self.dtype, copy=True)

    # Interpolate coordinates, re-combine to an x/y/z array and
    # add to the input points (if coordinates, the points are zeroed out above)
    shifts = np.vstack(
        (
            map_coordinates(self.field[:, :, :, 0], coords, order=1, mode=mode, cval=cval),
            map_coordinates(self.field[:, :, :, 1], coords, order=1, mode=mode, cval=cval),
            map_coordinates(self.field[:, :, :, 2], coords, order=1, mode=mode, cval=cval),
        )
    ).T.astype(self.dtype)
    points_xf += shifts

    return points_xf

`navis.transforms.H5transform` #

Hdf5 transform of 3D spatial data.

See here for specifications of the format.

PARAMETER	DESCRIPTION
`f`	`Path to Hdf5 transformation.` TYPE: `str`
`direction`	`Direction of transformation.` TYPE: `"forward" \| "inverse"` DEFAULT: `'forward'`
`level`	`For Hdf5 files with deformation fields at multiple resolutions: what level of detail to use. Negative values go backwards from the highest available resolution (-1 = highest, -2 = second highest, etc). Ignored if only a single deformation field present.` TYPE: `int` DEFAULT: `-1`
`cache`	`If True, we will cache the deformation field for subsequent future transforms. This will speed up future calculations in the future but comes at a memory cost.` TYPE: `bool` DEFAULT: `False`
`full_ingest`	`If True, will read and cache the full deformation field at initialization. This additional upfront cost can pay off if you are about to make many transforms across the volume.` TYPE: `bool` DEFAULT: `False`

Source code in navis/transforms/h5reg.py

class H5transform(BaseTransform):
    """Hdf5 transform of 3D spatial data.

    See [here](https://github.com/saalfeldlab/template-building/wiki/Hdf5-Deformation-fields)
    for specifications of the format.

    Parameters
    ----------
    f :             str
                    Path to Hdf5 transformation.
    direction :     "forward" | "inverse"
                    Direction of transformation.
    level :         int, optional
                    For Hdf5 files with deformation fields at multiple
                    resolutions: what level of detail to use. Negative values
                    go backwards from the highest available resolution
                    (-1 = highest, -2 = second highest, etc). Ignored if only a
                    single deformation field present.
    cache :         bool
                    If True, we will cache the deformation field for subsequent
                    future transforms. This will speed up future calculations
                    in the future but comes at a memory cost.
    full_ingest :   bool
                    If True, will read and cache the full deformation field at
                    initialization. This additional upfront cost can pay off if
                    you are about to make many transforms across the volume.

    """

    def __init__(
        self,
        f: str,
        direction: str = "forward",
        level: Optional[int] = -1,
        cache: bool = False,
        full_ingest: bool = False,
    ):
        """Init class."""
        assert direction in ("forward", "inverse"), (
            '`direction` must be "forward"' f'or "inverse", not "{direction}"'
        )

        self.file = f
        self.direction = direction
        self.field = {"forward": "dfield", "inverse": "invdfield"}[direction]

        # Trying to avoid the file repeatedly so we are making these initial
        # adjustments all in one go even though it would be more Pythonic to
        # delegate to property getter/setter methods
        with h5py.File(self.file, "r") as h5:
            # Get the available levels
            available_levels = []
            for k in h5.keys():
                try:
                    available_levels.append(int(k))
                except ValueError:
                    continue
            available_levels = sorted(available_levels)

            # Check if there are indeed deformation fields at various resolutions
            if available_levels:
                if isinstance(level, type(None)):
                    level = available_levels[0]
                elif level < 0:
                    ix = level * -1 - 1
                    level = available_levels[ix]

                # Set level
                self._level = str(level)

                # Shape of deformation field
                self._shape = h5[self.level][self.field].shape

                # Data type of deformation field
                self.dtype = h5[self.level][self.field].dtype
            elif self.field in h5.keys():
                # Set level
                self._level = None

                # Shape of deformation field
                self._shape = h5[self.field].shape

                # Data type of deformation field
                self.dtype = h5[self.field].dtype
            else:
                raise ValueError(
                    "Unable to parse deformation fields from " f" {self.file}."
                )

        # Prepare cache if applicable
        if full_ingest:
            # Ingest the whole deformation field
            self.full_ingest()
        elif cache:
            self.use_cache = True

    def __eq__(self, other) -> bool:
        """Compare with other Transform."""
        if isinstance(other, H5transform):
            if self.file == other.file:
                if self.direction == other.direction:
                    if self.level == other.level:
                        return True
        return False

    def __neg__(self) -> "H5transform":
        """Invert direction."""
        # Swap direction
        new_direction = {"forward": "inverse", "inverse": "forward"}[self.direction]
        # We will re-iniatialize
        x = H5transform(
            self.file,
            direction=new_direction,
            level=int(self.level) if self.level else None,
            cache=self.use_cache,
            full_ingest=False,
        )

        return x

    @property
    def level(self):
        return self._level

    @level.setter
    def level(self, value):
        raise ValueError("`level` cannot be changed after initialization.")

    @property
    def shape(self):
        """Shape of the deformation field.

        Note that the deformation field is likely to be (z, y, x, 3) where
        the last dimension contains the x/y/z offsets.
        """
        return self._shape

    @shape.setter
    def shape(self, value):
        raise ValueError("`shape` cannot be changed after initialization.")

    @property
    def spacing(self):
        """Voxel spacing of the deformation field (z, y, x)."""
        if not hasattr(self, "_spacing"):
            with h5py.File(self.file, "r") as h5:
                if self.level:
                    self._spacing = h5[self.level][self.field].attrs["spacing"]
                else:
                    self._spacing = h5[self.field].attrs["spacing"]
        return self._spacing

    @property
    def quantization_multiplier(self):
        """Quantization multiplier of the deformation field."""
        if not hasattr(self, "_quantization_multiplier"):
            with h5py.File(self.file, "r") as h5:
                if self.level:
                    self._quantization_multiplier = h5[self.level][
                        self.field
                    ].attrs.get("quantization_multiplier", 1)
                else:
                    self._quantization_multiplier = h5[self.field].attrs.get(
                        "quantization_multiplier", 1
                    )
        return self._quantization_multiplier

    @property
    def use_cache(self):
        """Whether to cache the deformation field."""
        if not hasattr(self, "_use_cache"):
            self._use_cache = False
        return self._use_cache

    @use_cache.setter
    def use_cache(self, value):
        """Set whether to cache the deformation field."""
        assert isinstance(value, bool)

        # If was False and now set to True, build the cache
        if not getattr(self, "_use_cache", False) and value:
            # This is the actual cache
            self.cache = np.zeros(self.shape, dtype=self.dtype)
            # This is a mask that tells us which values have already been cached
            self.cached = np.zeros(self.shape[:-1], dtype=bool)
            self._use_cache = True
        # If was True and now is set to False, deconstruct cache
        elif getattr(self, "_use_cache", False) and not value:
            del self.cache
            del self.cached
            self._use_cache = False

            if hasattr(self, "_fully_ingested"):
                del self._fully_ingested

        # Note: should explore whether we can use sparse N-dimensional
        # arrays for caching to save memory
        # See https://github.com/pydata/sparse/

    def copy(self, drop_cache=False):
        """Return copy.

        Please note that by default we're carrying over the cache *WITHOUT* copying it!
        Currently, `xform_brain` turns transforms into a `TransformSequence` which
        makes copies of each transform. If we didn't carry over the cache, we would never
        actually use it.
        """
        transform = H5transform(
            self.file,
            direction=self.direction,
            level=int(self.level) if self.level else None,
            cache=self.use_cache,
            full_ingest=False,
        )
        if not drop_cache and self.use_cache:
            transform.cache = self.cache
            transform.cached = self.cached
            if hasattr(self, "_fully_ingested"):
                transform._fully_ingested = self._fully_ingested

        return transform

    def full_ingest(self):
        """Fully ingest the deformation field."""
        # Skip if already ingested
        if getattr(self, "_fully_ingested", False):
            return

        with h5py.File(self.file, "r") as h5:
            # Read in the entire field
            if self.level:
                self.cache = h5[self.level][self.field][:, :, :]
            else:
                self.cache = h5[self.field][:, :, :]
            # Keep a flag of this
            self._fully_ingested = True
            # We set `cached` to True instead of using a mask
            self.cached = True
            # Keep track of the caching
            self._use_cache = True

    def precache(self, bbox: Union[list, np.ndarray], padding=True):
        """Cache deformation field for given bounding box.

        Parameters
        ----------
        bbox :      list | array
                    Must be `[[x1, x2], [y1, y2], [z1, z2]]`.
        padding :   bool
                    If True, will add the (required!) padding to the bounding
                    box.

        """
        bbox = np.asarray(bbox)

        if bbox.ndim != 2 or bbox.shape != (3, 2):
            raise ValueError(f"Expected (3, 2) bounding box, got {bbox.shape}")

        # Set use_cache=True -> this also prepares the cache array(s)
        self.use_cache = True

        with h5py.File(self.file, "r") as h5:
            if self.level:
                spacing = h5[self.level][self.field].attrs["spacing"]
            else:
                spacing = h5[self.field].attrs["spacing"]

            # Note that we invert because spacing is given in (z, y, x)
            bbox_vxl = (bbox.T / spacing[::-1]).T
            # Digitize into voxels
            bbox_vxl = bbox_vxl.round().astype(int)

            if padding:
                bbox_vxl[:, 0] -= 2
                bbox_vxl[:, 1] += 2

            # Make sure we are within bounds
            bbox_vxl = np.clip(bbox_vxl.T, 0, self.shape[:-1][::-1]).T

            # Extract values
            x1, x2, y1, y2, z1, z2 = bbox_vxl.flatten()

            # Cache values in this bounding box
            if self.level:
                read_from = h5[self.level][self.field]
            else:
                read_from = h5[self.field]

            self.cache[z1:z2, y1:y2, x1:x2] = read_from[z1:z2, y1:y2, x1:x2]
            self.cached[z1:z2, y1:y2, x1:x2] = True

    @staticmethod
    def from_file(filepath: str, **kwargs) -> "H5transform":
        """Generate H5transform from file.

        Parameters
        ----------
        filepath :  str
                    Path to H5 transform.
        **kwargs
                    Keyword arguments passed to H5transform.__init__

        Returns
        -------
        H5transform

        """
        return H5transform(str(filepath), **kwargs)

    def xform(
        self,
        points: np.ndarray,
        affine_fallback: bool = True,
        force_deform: bool = True,
    ) -> np.ndarray:
        """Xform data.

        Parameters
        ----------
        points :            (N, 3) numpy array | pandas.DataFrame
                            Points to xform. DataFrame must have x/y/z columns.
        affine_fallback :   bool
                            If False, points outside the deformation field will
                            be returned as `np.nan`. If True, points outside the
                            deformation field fall back to the baked-in affine
                            transform unless `force_deform=True`.
        force_deform :      bool
                            If True and `affine_fallback=True`, points outside
                            the deformation field are deformed using the closest
                            point inside the field, equivalent to nearest-neighbor
                            extrapolation of the deformation component. If False,
                            out-of-bounds points receive only the affine part.
                            Ignored if `affine_fallback=False`.

        Returns
        -------
        pointsxf :      (N, 3) numpy array
                        Transformed points. Out-of-bounds handling depends on
                        `affine_fallback` and `force_deform`.

        """
        if isinstance(points, pd.DataFrame):
            # Make sure x/y/z columns are present
            if np.any([c not in points for c in ["x", "y", "z"]]):
                raise ValueError("points DataFrame must have x/y/z columns.")
            points = points[["x", "y", "z"]].values

        if (
            not isinstance(points, np.ndarray)
            or points.ndim != 2
            or points.shape[1] != 3
        ):
            raise TypeError(
                "`points` must be numpy array of shape (N, 3) or "
                "pandas DataFrame with x/y/z columns"
            )

        # Read the file
        with h5py.File(self.file, "r") as h5:
            if self.level:
                field = h5[self.level][self.field]
            else:
                field = h5[self.field]

            # We need the field to be (z, y, x, offsets) with `offsets` being
            # three values - for example: (293, 470, 1010, 3)
            # If that's not the case, something is fishy!
            if field.shape[-1] != 3:
                logger.warning(
                    "Expected the deformation field to be of shape "
                    f"(z, y, x, 3), got {field.shape}."
                )

            if "affine" in field.attrs:
                # The affine part of the transform is a 4 x 4 matrix where the upper
                # 3 x 4 part (row x columns) is an attribute of the h5 dataset
                M = np.eye(4)
                M[:3, :4] = field.attrs["affine"].reshape(3, 4)
                affine = AffineTransform(M)
            else:
                affine = False

            # Get quantization multiplier for later use
            quantization_multiplier = field.attrs.get("quantization_multiplier", 1)

            # For inverse direction, the affine part is applied first
            if self.direction == "inverse" and affine:
                xf = affine.xform(points)
            else:
                xf = points

            affine_xf = xf.copy()

            # Translate points into voxel space
            spacing = field.attrs["spacing"]
            # Note that we invert because spacing is given in (z, y, x)
            xf_voxel = xf / spacing[::-1]
            # Digitize points into voxels
            xf_indices = xf_voxel.round().astype(int)
            # Determine the bounding box of the deformation vectors we need
            # Note that we are grabbing a bit more than required - this is
            # necessary for interpolation later down the line
            mn = xf_indices.min(axis=0) - 2
            mx = xf_indices.max(axis=0) + 2

            # Make sure we are within bounds
            # Note that we clip `mn` at 0 and `mx` at 2 at the lower end?
            # This is to make sure we have enough of the deformation field
            # to interpolate later on `offsets`
            mn = np.clip(mn, 2, np.array(self.shape[:-1][::-1])) - 2
            mx = np.clip(mx, 0, np.array(self.shape[:-1][::-1]) - 2) + 2

            # Check if we can use cached values
            if self.use_cache and (
                hasattr(self, "_fully_ingested")
                or np.all(self.cached[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]])
            ):
                offsets = self.cache[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]]
            else:
                # Load the deformation values for this bounding box
                # This is faster than grabbing individual voxels and
                offsets = field[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]]

                if self.use_cache:
                    # Write these offsets to cache
                    self.cache[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]] = offsets
                    self.cached[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]] = True

        # Before we interpolate check how many points are outside the
        # deformation field -> these will only receive the affine part of the
        # transform
        is_out = (xf_voxel.min(axis=1) < 0) | np.any(
            xf_voxel >= self.shape[:-1][::-1], axis=1
        )

        # If more than 20% (arbitrary number) of voxels are out, there is
        # something suspicious going on
        frac_out = is_out.sum() / xf_voxel.shape[0]
        if frac_out > 0.2 and not getattr(self, "_silenced_large_out_warning", False):
            logger.warning(
                f"A suspiciously large fraction ({frac_out:.1%}) "
                f"of {xf_voxel.shape[0]} points appear to be outside "
                "the H5 deformation field. Please make doubly sure "
                "that the input coordinates are in the correct "
                "space/units"
            )

        # map_coordinates expects coordinates in (z, y, x) order and relative
        # to the extracted offsets block.
        coords = xf_voxel[:, ::-1].T - mn[::-1].reshape(3, 1)
        mode = "nearest" if affine_fallback and force_deform else "constant"
        cval = 0 if mode == "nearest" else np.nan

        offset_vxl = np.vstack(
            (
                map_coordinates(offsets[:, :, :, 0], coords, order=1, mode=mode, cval=cval),
                map_coordinates(offsets[:, :, :, 1], coords, order=1, mode=mode, cval=cval),
                map_coordinates(offsets[:, :, :, 2], coords, order=1, mode=mode, cval=cval),
            )
        ).T

        # Turn offsets into real-world coordinates
        offset_real = offset_vxl * quantization_multiplier

        # Please note that we must not use += here to allow upcasting as needed.
        xf = xf + offset_real

        # For forward direction, the affine part is applied last
        if self.direction == "forward" and affine:
            xf = affine.xform(xf)
            affine_xf = affine.xform(affine_xf)

        if affine_fallback and not force_deform and is_out.any():
            xf[is_out, :] = affine_xf[is_out, :]

        # If no affine_fallback, set outside points to np.nan
        if not affine_fallback:
            xf[is_out, :] = np.nan

        return xf

    def xform_image(
        self,
        image: np.ndarray,
        image_res: Union[tuple, list],
        out_res: Optional[Union[tuple, list]] = None,
        out_shape: Optional[Union[tuple, list]] = None,
        order: int = 1,
        mode: str = "constant",
        cval: float = 0,
        cache: bool = True,
        progress: bool = True,
    ) -> np.ndarray:
        """Transform a 3D image using backward mapping.

        If available, this method will use `numba` to accelerate the
        transformation. This is highly recommended as it can speed up
        the transformation by several orders of magnitude:

          pip install numba

        Note though that the numba-accelerated path only supports linear
        interpolation and constant-mode boundary handling (default).

        Parameters
        ----------
        image :         (N, M, K) numpy array
                        Image to be transformed.
        image_res :     (3, ) tuple | list
                        Voxel resolution of the input image.
        out_res :       (3, ) tuple | list, optional
                        Voxel resolution of the output image. If None, assumed
                        to be the same as `image_res`.
        out_shape :     (3, ) tuple | list, optional
                        Shape of the output image in voxels. If None, uses the
                        shape of the deformation field. This works as long as
                        the deformation field fully samples the target space.
        order :         int
                        The order of the spline interpolation, default is 1
                        (linear). The order has to be in the range 0-5.
        mode :          str
                        How to handle values outside the image bounds.
                        Default is 'constant' (pad with cval).
        cval :          float
                        Value used for points outside the boundaries when
                        mode='constant'.
        cache :         bool
                        If True, we will cache the deformation field for
                        subsequent future transforms. This is generally
                        recommended unless memory is very tight or the
                        deformation field is huge.
        progress :      bool
                        Whether to show a progress bar during processing.
                        Not available (and probably not necessary) if
                        numba-accelerated path is used.

        Returns
        -------
        transformed :   (N, M, K) numpy array
                        Transformed image in target space. The shape is
                        determined by `out_shape` and `out_res`.

        Notes
        -----
        This method uses backward mapping: for each output chunk, we compute
        where the voxels came from in the source image using the transformation,
        then interpolate the source image at those locations. So if your
        transform goes from A -> B, this method actually uses the inverse transform
        under the hood to figure out where in A the voxels in B come from.

        """
        if not isinstance(image, np.ndarray) or image.ndim != 3:
            raise TypeError("`image` must be a 3D numpy array")

        # Because we are using backward mapping, we need to invert the transform,
        # i.e. we're taking points in the target space, mapping them back to the
        # source space and sampling the source image at those locations)
        if hasattr(self, "_inverted_copy"):
            # We're tracking the copy so that we can use the cache for subsequent
            # transformations instead of having to re-read the file and re-ingest
            # the deformation field
            reg_inv = self._inverted_copy
        else:
            reg_inv = -self
            # If cache, track the inverted copy (and the deformation field)
            if cache:
                reg_inv.use_cache = True
                self._inverted_copy = reg_inv

        if out_res is None:
            out_res = image_res

        if out_shape is None:
            # Bounds in physical space
            bounds = np.array(list(reg_inv.shape)[:-1][::-1]) * reg_inv.spacing
            # Back to voxel space in output resolution
            out_shape = np.ceil(bounds / out_res).astype(int)

        # Generate the empty output array
        out = np.zeros(out_shape, dtype=image.dtype)

        try:
            from .h5reg_numba import h5reg_warp_image_linear_constant

            nb_available = True
        except ModuleNotFoundError:
            nb_available = False

        # Numba path: linear interpolation + constant mode only but much faster
        if order == 1 and mode == "constant":
            if not nb_available:
                global NUMBA_WARNING
                if not NUMBA_WARNING:
                    logger.warning(
                        "For faster transforming of images, please install numba:\n"
                        "  pip install numba"
                    )
                    NUMBA_WARNING = True
            else:
                reg_inv.full_ingest()
                field = reg_inv.cache

                with h5py.File(reg_inv.file, "r") as h5:
                    if reg_inv.level:
                        ds = h5[reg_inv.level][reg_inv.field]
                    else:
                        ds = h5[reg_inv.field]

                    spacing = np.asarray(ds.attrs["spacing"], dtype=np.float64)
                    qmult = float(ds.attrs.get("quantization_multiplier", 1))

                    if "affine" in ds.attrs:
                        affine = np.eye(4, dtype=np.float64)
                        affine[:3, :4] = ds.attrs["affine"].reshape(3, 4)
                        apply_affine = 1 if reg_inv.direction == "inverse" else 2
                    else:
                        affine = np.eye(4, dtype=np.float64)
                        apply_affine = 0

                return h5reg_warp_image_linear_constant(
                    image,
                    field,
                    spacing,
                    qmult,
                    affine,
                    apply_affine,
                    out,
                    out_res,
                    image_res,
                    float(cval),
                )

        # Pre-compute prefilter if needed for higher order interpolation
        if order > 1:
            from scipy.ndimage import spline_filter

            image = spline_filter(image, order=order)

        for z in config.trange(
            out_shape[2],
            disable=not progress or config.pbar_hide,
            desc="Warping image",
        ):
            # Generate a grid of points in target voxel space for this z-slice
            chunk_coords = np.mgrid[
                0 : out_shape[0] : 1, 0 : out_shape[1] : 1, z : z + 1
            ]
            points = np.stack(
                [
                    chunk_coords[0].flatten(),
                    chunk_coords[1].flatten(),
                    chunk_coords[2].flatten(),
                ],
                axis=-1,
            )

            # Convert to physical space
            points = points.astype(np.float32) * out_res

            # Transform points to source space using the inverse transform
            try:
                reg_inv._silenced_large_out_warning = True
                points_xf = reg_inv.xform(points)
            except Exception as e:
                reg_inv._silenced_large_out_warning = False
                raise e

            # Convert back to voxel space in source image
            points_xf_vox = points_xf / image_res

            # Sample the source image at these transformed coordinates
            sampled_values = map_coordinates(
                image,
                [points_xf_vox[:, 0], points_xf_vox[:, 1], points_xf_vox[:, 2]],
                order=1,
                mode="constant",
                cval=0,
            )

            # Reshape back to 2D and assign to this slice of the output image
            out[:, :, z] = sampled_values.reshape(out_shape[0], out_shape[1])

        return out

`quantization_multiplier` `property` #

Quantization multiplier of the deformation field.

`shape` `property` `writable` #

Shape of the deformation field.

Note that the deformation field is likely to be (z, y, x, 3) where the last dimension contains the x/y/z offsets.

`spacing` `property` #

Voxel spacing of the deformation field (z, y, x).

`use_cache` `property` `writable` #

Whether to cache the deformation field.

`init` #

Init class.

Source code in navis/transforms/h5reg.py

def __init__(
    self,
    f: str,
    direction: str = "forward",
    level: Optional[int] = -1,
    cache: bool = False,
    full_ingest: bool = False,
):
    """Init class."""
    assert direction in ("forward", "inverse"), (
        '`direction` must be "forward"' f'or "inverse", not "{direction}"'
    )

    self.file = f
    self.direction = direction
    self.field = {"forward": "dfield", "inverse": "invdfield"}[direction]

    # Trying to avoid the file repeatedly so we are making these initial
    # adjustments all in one go even though it would be more Pythonic to
    # delegate to property getter/setter methods
    with h5py.File(self.file, "r") as h5:
        # Get the available levels
        available_levels = []
        for k in h5.keys():
            try:
                available_levels.append(int(k))
            except ValueError:
                continue
        available_levels = sorted(available_levels)

        # Check if there are indeed deformation fields at various resolutions
        if available_levels:
            if isinstance(level, type(None)):
                level = available_levels[0]
            elif level < 0:
                ix = level * -1 - 1
                level = available_levels[ix]

            # Set level
            self._level = str(level)

            # Shape of deformation field
            self._shape = h5[self.level][self.field].shape

            # Data type of deformation field
            self.dtype = h5[self.level][self.field].dtype
        elif self.field in h5.keys():
            # Set level
            self._level = None

            # Shape of deformation field
            self._shape = h5[self.field].shape

            # Data type of deformation field
            self.dtype = h5[self.field].dtype
        else:
            raise ValueError(
                "Unable to parse deformation fields from " f" {self.file}."
            )

    # Prepare cache if applicable
    if full_ingest:
        # Ingest the whole deformation field
        self.full_ingest()
    elif cache:
        self.use_cache = True

`copy` #

Return copy.

Please note that by default we're carrying over the cache WITHOUT copying it! Currently, xform_brain turns transforms into a TransformSequence which makes copies of each transform. If we didn't carry over the cache, we would never actually use it.

Source code in navis/transforms/h5reg.py

def copy(self, drop_cache=False):
    """Return copy.

    Please note that by default we're carrying over the cache *WITHOUT* copying it!
    Currently, `xform_brain` turns transforms into a `TransformSequence` which
    makes copies of each transform. If we didn't carry over the cache, we would never
    actually use it.
    """
    transform = H5transform(
        self.file,
        direction=self.direction,
        level=int(self.level) if self.level else None,
        cache=self.use_cache,
        full_ingest=False,
    )
    if not drop_cache and self.use_cache:
        transform.cache = self.cache
        transform.cached = self.cached
        if hasattr(self, "_fully_ingested"):
            transform._fully_ingested = self._fully_ingested

    return transform

`from_file` `staticmethod` #

Generate H5transform from file.

PARAMETER	DESCRIPTION
`filepath`	`Path to H5 transform.` TYPE: `str`
`**kwargs`	`Keyword arguments passed to H5transform.__init__` DEFAULT: `{}`

RETURNS	DESCRIPTION
`H5transform`

Source code in navis/transforms/h5reg.py

@staticmethod
def from_file(filepath: str, **kwargs) -> "H5transform":
    """Generate H5transform from file.

    Parameters
    ----------
    filepath :  str
                Path to H5 transform.
    **kwargs
                Keyword arguments passed to H5transform.__init__

    Returns
    -------
    H5transform

    """
    return H5transform(str(filepath), **kwargs)

`full_ingest` #

Fully ingest the deformation field.

Source code in navis/transforms/h5reg.py

def full_ingest(self):
    """Fully ingest the deformation field."""
    # Skip if already ingested
    if getattr(self, "_fully_ingested", False):
        return

    with h5py.File(self.file, "r") as h5:
        # Read in the entire field
        if self.level:
            self.cache = h5[self.level][self.field][:, :, :]
        else:
            self.cache = h5[self.field][:, :, :]
        # Keep a flag of this
        self._fully_ingested = True
        # We set `cached` to True instead of using a mask
        self.cached = True
        # Keep track of the caching
        self._use_cache = True

`precache` #

Cache deformation field for given bounding box.

PARAMETER	DESCRIPTION
`bbox`	Must be `[[x1, x2], [y1, y2], [z1, z2]]`. TYPE: `list \| array`
`padding`	`If True, will add the (required!) padding to the bounding box.` TYPE: `bool` DEFAULT: `True`

Source code in navis/transforms/h5reg.py

def precache(self, bbox: Union[list, np.ndarray], padding=True):
    """Cache deformation field for given bounding box.

    Parameters
    ----------
    bbox :      list | array
                Must be `[[x1, x2], [y1, y2], [z1, z2]]`.
    padding :   bool
                If True, will add the (required!) padding to the bounding
                box.

    """
    bbox = np.asarray(bbox)

    if bbox.ndim != 2 or bbox.shape != (3, 2):
        raise ValueError(f"Expected (3, 2) bounding box, got {bbox.shape}")

    # Set use_cache=True -> this also prepares the cache array(s)
    self.use_cache = True

    with h5py.File(self.file, "r") as h5:
        if self.level:
            spacing = h5[self.level][self.field].attrs["spacing"]
        else:
            spacing = h5[self.field].attrs["spacing"]

        # Note that we invert because spacing is given in (z, y, x)
        bbox_vxl = (bbox.T / spacing[::-1]).T
        # Digitize into voxels
        bbox_vxl = bbox_vxl.round().astype(int)

        if padding:
            bbox_vxl[:, 0] -= 2
            bbox_vxl[:, 1] += 2

        # Make sure we are within bounds
        bbox_vxl = np.clip(bbox_vxl.T, 0, self.shape[:-1][::-1]).T

        # Extract values
        x1, x2, y1, y2, z1, z2 = bbox_vxl.flatten()

        # Cache values in this bounding box
        if self.level:
            read_from = h5[self.level][self.field]
        else:
            read_from = h5[self.field]

        self.cache[z1:z2, y1:y2, x1:x2] = read_from[z1:z2, y1:y2, x1:x2]
        self.cached[z1:z2, y1:y2, x1:x2] = True

`xform` #

Xform data.

PARAMETER DESCRIPTION

points

            Points to xform. DataFrame must have x/y/z columns.

TYPE: (N, 3) numpy array | pandas.DataFrame

affine_fallback

            If False, points outside the deformation field will
            be returned as `np.nan`. If True, points outside the
            deformation field fall back to the baked-in affine
            transform unless `force_deform=True`.

TYPE: bool DEFAULT: True

force_deform

            If True and `affine_fallback=True`, points outside
            the deformation field are deformed using the closest
            point inside the field, equivalent to nearest-neighbor
            extrapolation of the deformation component. If False,
            out-of-bounds points receive only the affine part.
            Ignored if `affine_fallback=False`.

TYPE: bool DEFAULT: True

RETURNS	DESCRIPTION
`pointsxf`	Transformed points. Out-of-bounds handling depends on `affine_fallback` and `force_deform`. TYPE: `(N, 3) numpy array`

Source code in navis/transforms/h5reg.py

def xform(
    self,
    points: np.ndarray,
    affine_fallback: bool = True,
    force_deform: bool = True,
) -> np.ndarray:
    """Xform data.

    Parameters
    ----------
    points :            (N, 3) numpy array | pandas.DataFrame
                        Points to xform. DataFrame must have x/y/z columns.
    affine_fallback :   bool
                        If False, points outside the deformation field will
                        be returned as `np.nan`. If True, points outside the
                        deformation field fall back to the baked-in affine
                        transform unless `force_deform=True`.
    force_deform :      bool
                        If True and `affine_fallback=True`, points outside
                        the deformation field are deformed using the closest
                        point inside the field, equivalent to nearest-neighbor
                        extrapolation of the deformation component. If False,
                        out-of-bounds points receive only the affine part.
                        Ignored if `affine_fallback=False`.

    Returns
    -------
    pointsxf :      (N, 3) numpy array
                    Transformed points. Out-of-bounds handling depends on
                    `affine_fallback` and `force_deform`.

    """
    if isinstance(points, pd.DataFrame):
        # Make sure x/y/z columns are present
        if np.any([c not in points for c in ["x", "y", "z"]]):
            raise ValueError("points DataFrame must have x/y/z columns.")
        points = points[["x", "y", "z"]].values

    if (
        not isinstance(points, np.ndarray)
        or points.ndim != 2
        or points.shape[1] != 3
    ):
        raise TypeError(
            "`points` must be numpy array of shape (N, 3) or "
            "pandas DataFrame with x/y/z columns"
        )

    # Read the file
    with h5py.File(self.file, "r") as h5:
        if self.level:
            field = h5[self.level][self.field]
        else:
            field = h5[self.field]

        # We need the field to be (z, y, x, offsets) with `offsets` being
        # three values - for example: (293, 470, 1010, 3)
        # If that's not the case, something is fishy!
        if field.shape[-1] != 3:
            logger.warning(
                "Expected the deformation field to be of shape "
                f"(z, y, x, 3), got {field.shape}."
            )

        if "affine" in field.attrs:
            # The affine part of the transform is a 4 x 4 matrix where the upper
            # 3 x 4 part (row x columns) is an attribute of the h5 dataset
            M = np.eye(4)
            M[:3, :4] = field.attrs["affine"].reshape(3, 4)
            affine = AffineTransform(M)
        else:
            affine = False

        # Get quantization multiplier for later use
        quantization_multiplier = field.attrs.get("quantization_multiplier", 1)

        # For inverse direction, the affine part is applied first
        if self.direction == "inverse" and affine:
            xf = affine.xform(points)
        else:
            xf = points

        affine_xf = xf.copy()

        # Translate points into voxel space
        spacing = field.attrs["spacing"]
        # Note that we invert because spacing is given in (z, y, x)
        xf_voxel = xf / spacing[::-1]
        # Digitize points into voxels
        xf_indices = xf_voxel.round().astype(int)
        # Determine the bounding box of the deformation vectors we need
        # Note that we are grabbing a bit more than required - this is
        # necessary for interpolation later down the line
        mn = xf_indices.min(axis=0) - 2
        mx = xf_indices.max(axis=0) + 2

        # Make sure we are within bounds
        # Note that we clip `mn` at 0 and `mx` at 2 at the lower end?
        # This is to make sure we have enough of the deformation field
        # to interpolate later on `offsets`
        mn = np.clip(mn, 2, np.array(self.shape[:-1][::-1])) - 2
        mx = np.clip(mx, 0, np.array(self.shape[:-1][::-1]) - 2) + 2

        # Check if we can use cached values
        if self.use_cache and (
            hasattr(self, "_fully_ingested")
            or np.all(self.cached[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]])
        ):
            offsets = self.cache[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]]
        else:
            # Load the deformation values for this bounding box
            # This is faster than grabbing individual voxels and
            offsets = field[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]]

            if self.use_cache:
                # Write these offsets to cache
                self.cache[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]] = offsets
                self.cached[mn[2] : mx[2], mn[1] : mx[1], mn[0] : mx[0]] = True

    # Before we interpolate check how many points are outside the
    # deformation field -> these will only receive the affine part of the
    # transform
    is_out = (xf_voxel.min(axis=1) < 0) | np.any(
        xf_voxel >= self.shape[:-1][::-1], axis=1
    )

    # If more than 20% (arbitrary number) of voxels are out, there is
    # something suspicious going on
    frac_out = is_out.sum() / xf_voxel.shape[0]
    if frac_out > 0.2 and not getattr(self, "_silenced_large_out_warning", False):
        logger.warning(
            f"A suspiciously large fraction ({frac_out:.1%}) "
            f"of {xf_voxel.shape[0]} points appear to be outside "
            "the H5 deformation field. Please make doubly sure "
            "that the input coordinates are in the correct "
            "space/units"
        )

    # map_coordinates expects coordinates in (z, y, x) order and relative
    # to the extracted offsets block.
    coords = xf_voxel[:, ::-1].T - mn[::-1].reshape(3, 1)
    mode = "nearest" if affine_fallback and force_deform else "constant"
    cval = 0 if mode == "nearest" else np.nan

    offset_vxl = np.vstack(
        (
            map_coordinates(offsets[:, :, :, 0], coords, order=1, mode=mode, cval=cval),
            map_coordinates(offsets[:, :, :, 1], coords, order=1, mode=mode, cval=cval),
            map_coordinates(offsets[:, :, :, 2], coords, order=1, mode=mode, cval=cval),
        )
    ).T

    # Turn offsets into real-world coordinates
    offset_real = offset_vxl * quantization_multiplier

    # Please note that we must not use += here to allow upcasting as needed.
    xf = xf + offset_real

    # For forward direction, the affine part is applied last
    if self.direction == "forward" and affine:
        xf = affine.xform(xf)
        affine_xf = affine.xform(affine_xf)

    if affine_fallback and not force_deform and is_out.any():
        xf[is_out, :] = affine_xf[is_out, :]

    # If no affine_fallback, set outside points to np.nan
    if not affine_fallback:
        xf[is_out, :] = np.nan

    return xf

`xform_image` #

Transform a 3D image using backward mapping.

If available, this method will use numba to accelerate the transformation. This is highly recommended as it can speed up the transformation by several orders of magnitude:

pip install numba

Note though that the numba-accelerated path only supports linear interpolation and constant-mode boundary handling (default).

PARAMETER	DESCRIPTION
`image`	`Image to be transformed.` TYPE: `(N, M, K) numpy array`
`image_res`	`Voxel resolution of the input image.` TYPE: `(3, ) tuple \| list`
`out_res`	Voxel resolution of the output image. If None, assumed to be the same as `image_res`. TYPE: `(3, ) tuple \| list` DEFAULT: `None`
`out_shape`	`Shape of the output image in voxels. If None, uses the shape of the deformation field. This works as long as the deformation field fully samples the target space.` TYPE: `(3, ) tuple \| list` DEFAULT: `None`
`order`	`The order of the spline interpolation, default is 1 (linear). The order has to be in the range 0-5.` TYPE: `int` DEFAULT: `1`
`mode`	`How to handle values outside the image bounds. Default is 'constant' (pad with cval).` TYPE: `str` DEFAULT: `'constant'`
`cval`	`Value used for points outside the boundaries when mode='constant'.` TYPE: `float` DEFAULT: `0`
`cache`	`If True, we will cache the deformation field for subsequent future transforms. This is generally recommended unless memory is very tight or the deformation field is huge.` TYPE: `bool` DEFAULT: `True`
`progress`	`Whether to show a progress bar during processing. Not available (and probably not necessary) if numba-accelerated path is used.` TYPE: `bool` DEFAULT: `True`

RETURNS	DESCRIPTION
`transformed`	Transformed image in target space. The shape is determined by `out_shape` and `out_res`. TYPE: `(N, M, K) numpy array`

Notes

This method uses backward mapping: for each output chunk, we compute where the voxels came from in the source image using the transformation, then interpolate the source image at those locations. So if your transform goes from A -> B, this method actually uses the inverse transform under the hood to figure out where in A the voxels in B come from.

Source code in navis/transforms/h5reg.py

def xform_image(
    self,
    image: np.ndarray,
    image_res: Union[tuple, list],
    out_res: Optional[Union[tuple, list]] = None,
    out_shape: Optional[Union[tuple, list]] = None,
    order: int = 1,
    mode: str = "constant",
    cval: float = 0,
    cache: bool = True,
    progress: bool = True,
) -> np.ndarray:
    """Transform a 3D image using backward mapping.

    If available, this method will use `numba` to accelerate the
    transformation. This is highly recommended as it can speed up
    the transformation by several orders of magnitude:

      pip install numba

    Note though that the numba-accelerated path only supports linear
    interpolation and constant-mode boundary handling (default).

    Parameters
    ----------
    image :         (N, M, K) numpy array
                    Image to be transformed.
    image_res :     (3, ) tuple | list
                    Voxel resolution of the input image.
    out_res :       (3, ) tuple | list, optional
                    Voxel resolution of the output image. If None, assumed
                    to be the same as `image_res`.
    out_shape :     (3, ) tuple | list, optional
                    Shape of the output image in voxels. If None, uses the
                    shape of the deformation field. This works as long as
                    the deformation field fully samples the target space.
    order :         int
                    The order of the spline interpolation, default is 1
                    (linear). The order has to be in the range 0-5.
    mode :          str
                    How to handle values outside the image bounds.
                    Default is 'constant' (pad with cval).
    cval :          float
                    Value used for points outside the boundaries when
                    mode='constant'.
    cache :         bool
                    If True, we will cache the deformation field for
                    subsequent future transforms. This is generally
                    recommended unless memory is very tight or the
                    deformation field is huge.
    progress :      bool
                    Whether to show a progress bar during processing.
                    Not available (and probably not necessary) if
                    numba-accelerated path is used.

    Returns
    -------
    transformed :   (N, M, K) numpy array
                    Transformed image in target space. The shape is
                    determined by `out_shape` and `out_res`.

    Notes
    -----
    This method uses backward mapping: for each output chunk, we compute
    where the voxels came from in the source image using the transformation,
    then interpolate the source image at those locations. So if your
    transform goes from A -> B, this method actually uses the inverse transform
    under the hood to figure out where in A the voxels in B come from.

    """
    if not isinstance(image, np.ndarray) or image.ndim != 3:
        raise TypeError("`image` must be a 3D numpy array")

    # Because we are using backward mapping, we need to invert the transform,
    # i.e. we're taking points in the target space, mapping them back to the
    # source space and sampling the source image at those locations)
    if hasattr(self, "_inverted_copy"):
        # We're tracking the copy so that we can use the cache for subsequent
        # transformations instead of having to re-read the file and re-ingest
        # the deformation field
        reg_inv = self._inverted_copy
    else:
        reg_inv = -self
        # If cache, track the inverted copy (and the deformation field)
        if cache:
            reg_inv.use_cache = True
            self._inverted_copy = reg_inv

    if out_res is None:
        out_res = image_res

    if out_shape is None:
        # Bounds in physical space
        bounds = np.array(list(reg_inv.shape)[:-1][::-1]) * reg_inv.spacing
        # Back to voxel space in output resolution
        out_shape = np.ceil(bounds / out_res).astype(int)

    # Generate the empty output array
    out = np.zeros(out_shape, dtype=image.dtype)

    try:
        from .h5reg_numba import h5reg_warp_image_linear_constant

        nb_available = True
    except ModuleNotFoundError:
        nb_available = False

    # Numba path: linear interpolation + constant mode only but much faster
    if order == 1 and mode == "constant":
        if not nb_available:
            global NUMBA_WARNING
            if not NUMBA_WARNING:
                logger.warning(
                    "For faster transforming of images, please install numba:\n"
                    "  pip install numba"
                )
                NUMBA_WARNING = True
        else:
            reg_inv.full_ingest()
            field = reg_inv.cache

            with h5py.File(reg_inv.file, "r") as h5:
                if reg_inv.level:
                    ds = h5[reg_inv.level][reg_inv.field]
                else:
                    ds = h5[reg_inv.field]

                spacing = np.asarray(ds.attrs["spacing"], dtype=np.float64)
                qmult = float(ds.attrs.get("quantization_multiplier", 1))

                if "affine" in ds.attrs:
                    affine = np.eye(4, dtype=np.float64)
                    affine[:3, :4] = ds.attrs["affine"].reshape(3, 4)
                    apply_affine = 1 if reg_inv.direction == "inverse" else 2
                else:
                    affine = np.eye(4, dtype=np.float64)
                    apply_affine = 0

            return h5reg_warp_image_linear_constant(
                image,
                field,
                spacing,
                qmult,
                affine,
                apply_affine,
                out,
                out_res,
                image_res,
                float(cval),
            )

    # Pre-compute prefilter if needed for higher order interpolation
    if order > 1:
        from scipy.ndimage import spline_filter

        image = spline_filter(image, order=order)

    for z in config.trange(
        out_shape[2],
        disable=not progress or config.pbar_hide,
        desc="Warping image",
    ):
        # Generate a grid of points in target voxel space for this z-slice
        chunk_coords = np.mgrid[
            0 : out_shape[0] : 1, 0 : out_shape[1] : 1, z : z + 1
        ]
        points = np.stack(
            [
                chunk_coords[0].flatten(),
                chunk_coords[1].flatten(),
                chunk_coords[2].flatten(),
            ],
            axis=-1,
        )

        # Convert to physical space
        points = points.astype(np.float32) * out_res

        # Transform points to source space using the inverse transform
        try:
            reg_inv._silenced_large_out_warning = True
            points_xf = reg_inv.xform(points)
        except Exception as e:
            reg_inv._silenced_large_out_warning = False
            raise e

        # Convert back to voxel space in source image
        points_xf_vox = points_xf / image_res

        # Sample the source image at these transformed coordinates
        sampled_values = map_coordinates(
            image,
            [points_xf_vox[:, 0], points_xf_vox[:, 1], points_xf_vox[:, 2]],
            order=1,
            mode="constant",
            cval=0,
        )

        # Reshape back to 2D and assign to this slice of the output image
        out[:, :, z] = sampled_values.reshape(out_shape[0], out_shape[1])

    return out

`navis.transforms.MovingLeastSquaresTransform` #

Source code in navis/transforms/moving_least_squares.py

class MovingLeastSquaresTransform(BackendMixin, BaseTransform):
    _fallback_backend = "python"

    def __init__(
        self,
        landmarks_source: np.ndarray,
        landmarks_target: np.ndarray,
        direction: str = "forward",
        batch_size: int = 100_000,
        backend: str = None,
    ) -> None:
        """Moving Least Squares transforms of 3D spatial data.

        Runs on navis-fastcore's Rust implementation of the affine algorithm
        published in [Schaefer et al. 2006](https://dl.acm.org/doi/pdf/10.1145/1179352.1141920).
        The deprecated "python" backend runs on
        [molesq](https://github.com/clbarnes/molesq) instead, which packages the
        [implementation](https://github.com/ceesem/catalysis/blob/master/catalysis/transform.py)
        by Casey Schneider-Mizell of the same algorithm - see the `backend`
        parameter. `molesq` is imported on demand and is only needed for that
        backend.

        Notes
        -----
        At least in my hands, `TPStransforms` are significantly faster than
        `MovingLeastSquaresTransforms`. The results are similar but not identical,
        so make sure to use the one that works best for your use case. The two
        backends agree with each other to ~1e-13.

        Parameters
        ----------
        landmarks_source : np.ndarray
            Source landmarks as x/y/z coordinates.
        landmarks_target : np.ndarray
            Target landmarks as x/y/z coordinates.
        direction : str
            'forward' (default) or 'inverse' (treat the target as the source and vice versa)
        batch_size : int, optional
            Batch size for transforming points. At one point during the transformation,
            molesq generates a (N, M) distance matrix, where N is the number of landmarks
            and M is the number of locations, which can get prohibitively expensive.
            We avoid the issue by batching the transformation by default. Note that the
            overhead from batching seems negligible.
            Ignored on the "fastcore" backend, which reduces over landmarks
            instead of materialising that matrix.
        backend : "auto" | "fastcore" | "python", optional
            Which implementation to use. `None` (default) defers to
            `navis.config.default_transform_backend`. "python" is deprecated and
            will be removed in 3.0.

        Examples
        --------
        >>> from navis import transforms
        >>> import numpy as np
        >>> # Generate some mock landmarks
        >>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
        >>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
        >>> tr = transforms.MovingLeastSquaresTransform(src, trg)
        >>> points = np.array([[0, 0, 0], [50, 50, 50]])
        >>> tr.xform(points)                                        # doctest: +SKIP
        array([[  1.        ,  15.        ,   5.        ],
               [ 81.56361725, 155.32071504, 187.3147564 ]])
        >>> # The global affine part of the transform
        >>> tr.matrix_affine.shape
        (4, 4)

        """
        assert direction in ("forward", "inverse")

        self.source = np.asarray(landmarks_source)
        self.target = np.asarray(landmarks_target)

        # Same checks molesq's `Transformer` would do, kept here so that they
        # still happen up front now that it is only built on demand.
        if self.source.shape != self.target.shape:
            raise ValueError("Landmarks must have the same shape")
        if self.source.ndim != 2:
            raise ValueError("Landmarks must be a 2D array")

        self.reverse = direction == "inverse"
        self.batch_size = int(batch_size)
        self._backend = backend
        self._fc = None
        self._transformer = None

    @property
    def transformer(self):
        """The (lazily built) molesq transformer used by the "python" backend.

        molesq is imported here rather than at module level so that `import
        navis` does not pay for it and, more to the point, so that an install
        without it can still run moving least squares on the default (fastcore)
        backend - which needs nothing but the landmarks.

        """
        if self._transformer is None:
            try:
                from molesq import Transformer
            except ModuleNotFoundError:
                raise ModuleNotFoundError(
                    'The deprecated "python" backend for moving least squares '
                    "requires the `molesq` library:\n"
                    "  pip3 install molesq\n"
                    'Drop the explicit backend (or use "fastcore") to do '
                    "without it."
                )
            self._transformer = Transformer(self.source, self.target)
        return self._transformer

    @property
    def _fastcore_mls(self):
        """The (lazily built) fastcore transform.

        Moving least squares has no fit - construction just stores the
        landmarks - so building this is cheap. It is lazy only so that an
        install without fastcore never touches the attribute.

        Always built in the forward direction; `xform` passes `reverse` per
        call. That matters because navis keeps the direction in a plain
        `self.reverse` which `__neg__` (and users) flip in place, so a cached
        object with the direction baked in would go stale.

        """
        if self._fc is None:
            from .. import utils

            self._fc = utils.fastcore.MlsTransform(self.source, self.target)
        return self._fc

    def xform(self, points: np.ndarray) -> np.ndarray:
        """Transform points.

        Parameters
        ----------
        points :    (N, 3) array
                    Points to transform.

        Returns
        -------
        pointsxf :  (N, 3) array
                    Transformed points.

        """
        if isinstance(points, pd.DataFrame):
            if any([c not in points for c in ["x", "y", "z"]]):
                raise ValueError("DataFrame must have x/y/z columns.")
            points = points[["x", "y", "z"]].values

        if self.backend == "fastcore":
            # No batching here: everything but the result is a reduction over
            # landmarks, so peak memory is just the output array.
            return self._fastcore_mls.xform(
                np.asarray(points), reverse=self.reverse
            )

        batch_size = self.batch_size if self.batch_size else len(points)
        points_xf = []
        for i in range(0, len(points), batch_size):
            batch = points[i : i + batch_size]
            points_xf.append(self.transformer.transform(batch, reverse=self.reverse))

        return np.concatenate(points_xf, axis=0)

    @property
    def matrix_affine(self):
        """Return the affine transformation matrix.

        Note that moving least squares is a *locally* weighted affine transform:
        every point effectively gets its own affine matrix. This property returns
        the *global* affine, i.e. the least-squares fit of the source onto the
        target landmarks. That is the transform the moving least squares converges
        to far away from the landmarks, where the distance weights even out.

        """
        source, target = self._control_points()
        source = np.asarray(source, dtype=float)
        target = np.asarray(target, dtype=float)
        ndim = source.shape[1]

        # Least-squares fit of `source_homogeneous @ coefs = target`
        source_hom = np.ones((source.shape[0], ndim + 1))
        source_hom[:, :ndim] = source
        coefs = np.linalg.lstsq(source_hom, target, rcond=None)[0]

        # Combine into a typical (4, 4) transformation matrix
        # where the last row is [0, 0, 0, 1]
        m = np.eye(ndim + 1)
        m[:ndim, :ndim] = coefs[:ndim].T
        m[:ndim, ndim] = coefs[ndim]
        return m

    def __neg__(self) -> "MovingLeastSquaresTransform":
        """Invert direction"""
        out = self.copy()
        out.reverse = not self.reverse
        return out

    def __eq__(self, o: object) -> bool:
        if not isinstance(o, MovingLeastSquaresTransform):
            return False
        for cp_this, cp_that in zip(self._control_points(), o._control_points()):
            if not np.array_equal(cp_this, cp_that):
                return False
        return True

    def _control_points(self):
        cp1, cp2 = self.source, self.target
        if self.reverse:
            cp2, cp1 = cp1, cp2
        return cp1, cp2

    def copy(self):
        """Make a copy"""
        return deepcopy(self)

`matrix_affine` `property` #

Return the affine transformation matrix.

Note that moving least squares is a locally weighted affine transform: every point effectively gets its own affine matrix. This property returns the global affine, i.e. the least-squares fit of the source onto the target landmarks. That is the transform the moving least squares converges to far away from the landmarks, where the distance weights even out.

`transformer` `property` #

The (lazily built) molesq transformer used by the "python" backend.

molesq is imported here rather than at module level so that import navis does not pay for it and, more to the point, so that an install without it can still run moving least squares on the default (fastcore) backend - which needs nothing but the landmarks.

`init` #

Moving Least Squares transforms of 3D spatial data.

Runs on navis-fastcore's Rust implementation of the affine algorithm published in Schaefer et al. 2006. The deprecated "python" backend runs on molesq instead, which packages the implementation by Casey Schneider-Mizell of the same algorithm - see the backend parameter. molesq is imported on demand and is only needed for that backend.

Notes

At least in my hands, TPStransforms are significantly faster than MovingLeastSquaresTransforms. The results are similar but not identical, so make sure to use the one that works best for your use case. The two backends agree with each other to ~1e-13.

PARAMETER	DESCRIPTION
`landmarks_source`	Source landmarks as x/y/z coordinates. TYPE: `np.ndarray`
`landmarks_target`	Target landmarks as x/y/z coordinates. TYPE: `np.ndarray`
`direction`	'forward' (default) or 'inverse' (treat the target as the source and vice versa) TYPE: `str` DEFAULT: `'forward'`
`batch_size`	Batch size for transforming points. At one point during the transformation, molesq generates a (N, M) distance matrix, where N is the number of landmarks and M is the number of locations, which can get prohibitively expensive. We avoid the issue by batching the transformation by default. Note that the overhead from batching seems negligible. Ignored on the "fastcore" backend, which reduces over landmarks instead of materialising that matrix. TYPE: `int` DEFAULT: `100000`
`backend`	Which implementation to use. `None` (default) defers to `navis.config.default_transform_backend`. "python" is deprecated and will be removed in 3.0. TYPE: `auto \| fastcore \| python` DEFAULT: `None`

Examples:

>>> from navis import transforms
>>> import numpy as np
>>> # Generate some mock landmarks
>>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
>>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
>>> tr = transforms.MovingLeastSquaresTransform(src, trg)
>>> points = np.array([[0, 0, 0], [50, 50, 50]])
>>> tr.xform(points)
array([[  1.        ,  15.        ,   5.        ],
       [ 81.56361725, 155.32071504, 187.3147564 ]])
>>> # The global affine part of the transform
>>> tr.matrix_affine.shape
(4, 4)

Source code in navis/transforms/moving_least_squares.py

def __init__(
    self,
    landmarks_source: np.ndarray,
    landmarks_target: np.ndarray,
    direction: str = "forward",
    batch_size: int = 100_000,
    backend: str = None,
) -> None:
    """Moving Least Squares transforms of 3D spatial data.

    Runs on navis-fastcore's Rust implementation of the affine algorithm
    published in [Schaefer et al. 2006](https://dl.acm.org/doi/pdf/10.1145/1179352.1141920).
    The deprecated "python" backend runs on
    [molesq](https://github.com/clbarnes/molesq) instead, which packages the
    [implementation](https://github.com/ceesem/catalysis/blob/master/catalysis/transform.py)
    by Casey Schneider-Mizell of the same algorithm - see the `backend`
    parameter. `molesq` is imported on demand and is only needed for that
    backend.

    Notes
    -----
    At least in my hands, `TPStransforms` are significantly faster than
    `MovingLeastSquaresTransforms`. The results are similar but not identical,
    so make sure to use the one that works best for your use case. The two
    backends agree with each other to ~1e-13.

    Parameters
    ----------
    landmarks_source : np.ndarray
        Source landmarks as x/y/z coordinates.
    landmarks_target : np.ndarray
        Target landmarks as x/y/z coordinates.
    direction : str
        'forward' (default) or 'inverse' (treat the target as the source and vice versa)
    batch_size : int, optional
        Batch size for transforming points. At one point during the transformation,
        molesq generates a (N, M) distance matrix, where N is the number of landmarks
        and M is the number of locations, which can get prohibitively expensive.
        We avoid the issue by batching the transformation by default. Note that the
        overhead from batching seems negligible.
        Ignored on the "fastcore" backend, which reduces over landmarks
        instead of materialising that matrix.
    backend : "auto" | "fastcore" | "python", optional
        Which implementation to use. `None` (default) defers to
        `navis.config.default_transform_backend`. "python" is deprecated and
        will be removed in 3.0.

    Examples
    --------
    >>> from navis import transforms
    >>> import numpy as np
    >>> # Generate some mock landmarks
    >>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
    >>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
    >>> tr = transforms.MovingLeastSquaresTransform(src, trg)
    >>> points = np.array([[0, 0, 0], [50, 50, 50]])
    >>> tr.xform(points)                                        # doctest: +SKIP
    array([[  1.        ,  15.        ,   5.        ],
           [ 81.56361725, 155.32071504, 187.3147564 ]])
    >>> # The global affine part of the transform
    >>> tr.matrix_affine.shape
    (4, 4)

    """
    assert direction in ("forward", "inverse")

    self.source = np.asarray(landmarks_source)
    self.target = np.asarray(landmarks_target)

    # Same checks molesq's `Transformer` would do, kept here so that they
    # still happen up front now that it is only built on demand.
    if self.source.shape != self.target.shape:
        raise ValueError("Landmarks must have the same shape")
    if self.source.ndim != 2:
        raise ValueError("Landmarks must be a 2D array")

    self.reverse = direction == "inverse"
    self.batch_size = int(batch_size)
    self._backend = backend
    self._fc = None
    self._transformer = None

`copy` #

Make a copy

Source code in navis/transforms/moving_least_squares.py

def copy(self):
    """Make a copy"""
    return deepcopy(self)

`xform` #

Transform points.

PARAMETER	DESCRIPTION
`points`	`Points to transform.` TYPE: `(N, 3) array`

RETURNS	DESCRIPTION
`pointsxf`	Transformed points. TYPE: `(N, 3) array`

Source code in navis/transforms/moving_least_squares.py

def xform(self, points: np.ndarray) -> np.ndarray:
    """Transform points.

    Parameters
    ----------
    points :    (N, 3) array
                Points to transform.

    Returns
    -------
    pointsxf :  (N, 3) array
                Transformed points.

    """
    if isinstance(points, pd.DataFrame):
        if any([c not in points for c in ["x", "y", "z"]]):
            raise ValueError("DataFrame must have x/y/z columns.")
        points = points[["x", "y", "z"]].values

    if self.backend == "fastcore":
        # No batching here: everything but the result is a reduction over
        # landmarks, so peak memory is just the output array.
        return self._fastcore_mls.xform(
            np.asarray(points), reverse=self.reverse
        )

    batch_size = self.batch_size if self.batch_size else len(points)
    points_xf = []
    for i in range(0, len(points), batch_size):
        batch = points[i : i + batch_size]
        points_xf.append(self.transformer.transform(batch, reverse=self.reverse))

    return np.concatenate(points_xf, axis=0)

`navis.transforms.SimilarityTransform` #

Calculate and apply a similarity transform.

Similarity transforms are a special case of affine transforms that consist of a rotation, a translation and a single uniform scaling factor but no shear and no reflection.

The transform is the least-squares optimal fit between the two sets of landmarks and is found in closed form via a singular value decomposition (see Umeyama, 1991). Because a similarity transform has only 7 degrees of freedom (3 rotation, 3 translation, 1 scale), it will generally not map the source landmarks onto the target landmarks exactly - unlike, say, a navis.transforms.thinplate.TPStransform.

PARAMETER	DESCRIPTION
`landmarks_source`	`Source landmarks as x/y/z coordinates.` TYPE: `(M, 3) array`
`landmarks_target`	Target landmarks as x/y/z coordinates. Must have the same number of landmarks as `landmarks_source`. TYPE: `(M, 3) array`
`weights`	Relative weight for each landmark pair. If `None`, all landmarks are weighted equally. TYPE: `(M, ) array` DEFAULT: `None`
`scale`	`Whether to also estimate a uniform scaling factor. If False, the transform is a rigid transform (rotation + translation only).` TYPE: `bool` DEFAULT: `True`

Examples:

>>> from navis import transforms
>>> import numpy as np
>>> # Generate some mock landmarks
>>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
>>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
>>> tr = transforms.similarity.SimilarityTransform(src, trg)
>>> tr.xform(np.array([[0, 0, 0], [50, 50, 50]]))
array([[-0.27589955, 10.98503005,  4.26630007],
       [38.44936026, 53.04761397, 62.7346341 ]])

Note that - unlike a thin-plate spline transform - the landmarks are not mapped exactly because a similarity transform has only 7 degrees of freedom:

>>> round(tr.scale_factor, 4)
0.9443

Source code in navis/transforms/similarity.py

class SimilarityTransform(BaseTransform):
    """Calculate and apply a similarity transform.

    Similarity transforms are a special case of affine transforms that consist
    of a rotation, a translation and a single uniform scaling factor but no
    shear and no reflection.

    The transform is the least-squares optimal fit between the two sets of
    landmarks and is found in closed form via a singular value decomposition
    (see Umeyama, 1991). Because a similarity transform has only 7 degrees of
    freedom (3 rotation, 3 translation, 1 scale), it will generally *not* map
    the source landmarks onto the target landmarks exactly - unlike, say, a
    [`navis.transforms.thinplate.TPStransform`][].

    Parameters
    ----------
    landmarks_source :  (M, 3) array
                        Source landmarks as x/y/z coordinates.
    landmarks_target :  (M, 3) array
                        Target landmarks as x/y/z coordinates. Must have the
                        same number of landmarks as `landmarks_source`.
    weights :           (M, ) array, optional
                        Relative weight for each landmark pair. If `None`, all
                        landmarks are weighted equally.
    scale :             bool
                        Whether to also estimate a uniform scaling factor. If
                        False, the transform is a rigid transform (rotation +
                        translation only).

    Examples
    --------
    >>> from navis import transforms
    >>> import numpy as np
    >>> # Generate some mock landmarks
    >>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
    >>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
    >>> tr = transforms.similarity.SimilarityTransform(src, trg)
    >>> tr.xform(np.array([[0, 0, 0], [50, 50, 50]]))
    array([[-0.27589955, 10.98503005,  4.26630007],
           [38.44936026, 53.04761397, 62.7346341 ]])

    Note that - unlike a thin-plate spline transform - the landmarks are not
    mapped exactly because a similarity transform has only 7 degrees of freedom:

    >>> round(tr.scale_factor, 4)
    0.9443

    """

    def __init__(
        self,
        landmarks_source: np.ndarray,
        landmarks_target: np.ndarray,
        weights: np.ndarray = None,
        scale: bool = True,
    ):
        self.source = landmarks_source
        self.target = landmarks_target

        if self.source.ndim != 2 or self.source.shape[1] != 3:
            raise ValueError(f"Expected (M, 3) array, got {self.source.shape}")
        if self.target.ndim != 2 or self.target.shape[1] != 3:
            raise ValueError(f"Expected (M, 3) array, got {self.target.shape}")

        if self.source.shape[0] != self.target.shape[0]:
            raise ValueError(
                "Number of source landmarks must match number of target "
                f"landmarks: {self.source.shape[0]} != {self.target.shape[0]}"
            )

        if self.n <= 1:
            raise ValueError("Need at least two landmarks.")

        if weights is None:
            weights = np.ones(self.n)
        else:
            weights = np.asarray(weights, dtype=float)
            if weights.shape != (self.n,):
                raise ValueError(
                    f"Expected {self.n} weights, got {weights.shape[0]}"
                )
            if np.any(weights < 0):
                raise ValueError("Weights must not be negative.")
            if not weights.sum() > 0:
                raise ValueError("Weights must not sum to zero.")
        self.weights = weights

        self.scale = scale

        self.update()

    def __eq__(self, other) -> bool:
        """Implement equality comparison."""
        if not isinstance(other, SimilarityTransform):
            return False
        return (
            self.source.shape == other.source.shape
            and np.all(self.source == other.source)
            and np.all(self.target == other.target)
            and np.all(self.weights == other.weights)
            and self.scale == other.scale
        )

    def __neg__(self) -> "SimilarityTransform":
        """Invert direction."""
        # Simply switching source and target would re-fit the transform which
        # - for an imperfect fit - is not the same as inverting it. Instead we
        # invert the matrix itself.
        x = self.copy()
        x.source, x.target = self.target, self.source
        x.matrix = np.linalg.inv(self.matrix)
        return x

    @property
    def n(self) -> int:
        """Number of landmarks."""
        return self.source.shape[0]

    @property
    def source(self) -> np.ndarray:
        return self._source

    @source.setter
    def source(self, value):
        self._source = np.asarray(value, dtype=float)

    @property
    def target(self) -> np.ndarray:
        return self._target

    @target.setter
    def target(self, value):
        self._target = np.asarray(value, dtype=float)

    def update(self):
        """(Re-)calculate the transform from the current landmarks."""
        # Normalize weights so they sum to 1
        w = self.weights / self.weights.sum()

        # Weighted centroids
        mu_src = w @ self.source
        mu_trg = w @ self.target

        # Center the landmarks
        src_c = self.source - mu_src
        trg_c = self.target - mu_trg

        # Weighted cross-covariance (target x source)
        cov = (trg_c * w[:, None]).T @ src_c

        U, D, Vt = np.linalg.svd(cov)

        # Guard against the SVD producing a reflection instead of a rotation:
        # flipping the sign of the least significant singular value gives the
        # best proper rotation (det = +1).
        S = np.ones(3)
        if np.linalg.det(U) * np.linalg.det(Vt) < 0:
            S[-1] = -1

        R = U @ np.diag(S) @ Vt

        if self.scale:
            var_src = (w * (src_c**2).sum(axis=1)).sum()
            # If the source landmarks are all coincident there is no scale to
            # recover - fall back to 1 rather than dividing by zero.
            c = (D * S).sum() / var_src if var_src > 0 else 1.0
        else:
            c = 1.0

        matrix = np.eye(4)
        matrix[:3, :3] = c * R
        matrix[:3, 3] = mu_trg - c * R @ mu_src

        self.matrix = matrix

    @property
    def rotation(self) -> np.ndarray:
        """The (3, 3) rotation matrix."""
        return self.matrix[:3, :3] / self.scale_factor

    @property
    def scale_factor(self) -> float:
        """The uniform scaling factor."""
        # The matrix' linear part is `scale * R` and R is orthonormal, so each
        # of its columns has length `scale`.
        return float(np.linalg.norm(self.matrix[:3, 0]))

    @property
    def translation(self) -> np.ndarray:
        """The (3, ) translation vector."""
        return self.matrix[:3, 3]

    def copy(self) -> "SimilarityTransform":
        """Return copy of transform."""
        x = self.__class__.__new__(self.__class__)
        x.__dict__.update({k: copy.copy(v) for k, v in self.__dict__.items()})
        return x

    def xform(self, points: np.ndarray) -> np.ndarray:
        """Transform points.

        Parameters
        ----------
        points :    (N, 3) array | pandas.DataFrame
                    Points to transform. DataFrames must have x/y/z columns.

        Returns
        -------
        pointsxf :  (N, 3) array
                    Transformed points.

        """
        if isinstance(points, pd.DataFrame):
            if any(c not in points for c in ["x", "y", "z"]):
                raise ValueError("DataFrame must have x/y/z columns.")
            points = points[["x", "y", "z"]].values

        points = np.asarray(points, dtype=float)

        if points.ndim != 2 or points.shape[1] != 3:
            raise ValueError(f"Expected (N, 3) array, got {points.shape}")

        return points @ self.matrix[:3, :3].T + self.matrix[:3, 3]

`n: int` `property` #

Number of landmarks.

`rotation: np.ndarray` `property` #

The (3, 3) rotation matrix.

`scale_factor: float` `property` #

The uniform scaling factor.

`translation: np.ndarray` `property` #

The (3, ) translation vector.

`copy` #

Return copy of transform.

Source code in navis/transforms/similarity.py

def copy(self) -> "SimilarityTransform":
    """Return copy of transform."""
    x = self.__class__.__new__(self.__class__)
    x.__dict__.update({k: copy.copy(v) for k, v in self.__dict__.items()})
    return x

`update` #

(Re-)calculate the transform from the current landmarks.

Source code in navis/transforms/similarity.py

def update(self):
    """(Re-)calculate the transform from the current landmarks."""
    # Normalize weights so they sum to 1
    w = self.weights / self.weights.sum()

    # Weighted centroids
    mu_src = w @ self.source
    mu_trg = w @ self.target

    # Center the landmarks
    src_c = self.source - mu_src
    trg_c = self.target - mu_trg

    # Weighted cross-covariance (target x source)
    cov = (trg_c * w[:, None]).T @ src_c

    U, D, Vt = np.linalg.svd(cov)

    # Guard against the SVD producing a reflection instead of a rotation:
    # flipping the sign of the least significant singular value gives the
    # best proper rotation (det = +1).
    S = np.ones(3)
    if np.linalg.det(U) * np.linalg.det(Vt) < 0:
        S[-1] = -1

    R = U @ np.diag(S) @ Vt

    if self.scale:
        var_src = (w * (src_c**2).sum(axis=1)).sum()
        # If the source landmarks are all coincident there is no scale to
        # recover - fall back to 1 rather than dividing by zero.
        c = (D * S).sum() / var_src if var_src > 0 else 1.0
    else:
        c = 1.0

    matrix = np.eye(4)
    matrix[:3, :3] = c * R
    matrix[:3, 3] = mu_trg - c * R @ mu_src

    self.matrix = matrix

`xform` #

Transform points.

PARAMETER	DESCRIPTION
`points`	`Points to transform. DataFrames must have x/y/z columns.` TYPE: `(N, 3) array \| pandas.DataFrame`

RETURNS	DESCRIPTION
`pointsxf`	Transformed points. TYPE: `(N, 3) array`

Source code in navis/transforms/similarity.py

def xform(self, points: np.ndarray) -> np.ndarray:
    """Transform points.

    Parameters
    ----------
    points :    (N, 3) array | pandas.DataFrame
                Points to transform. DataFrames must have x/y/z columns.

    Returns
    -------
    pointsxf :  (N, 3) array
                Transformed points.

    """
    if isinstance(points, pd.DataFrame):
        if any(c not in points for c in ["x", "y", "z"]):
            raise ValueError("DataFrame must have x/y/z columns.")
        points = points[["x", "y", "z"]].values

    points = np.asarray(points, dtype=float)

    if points.ndim != 2 or points.shape[1] != 3:
        raise ValueError(f"Expected (N, 3) array, got {points.shape}")

    return points @ self.matrix[:3, :3].T + self.matrix[:3, 3]

`navis.transforms.TPStransform` #

Thin Plate Spline transforms of 3D spatial data.

Runs on navis-fastcore's Rust implementation; the deprecated "python" backend runs on morphops instead (see the backend parameter). The spline is always fitted with morphops - numpy's LAPACK-backed solve is faster than fastcore's here - so the backend only changes how points are transformed, which is where fastcore is ~10-15x faster. The two agree to ~1e-13.

Notes

At least in my hands, TPStransforms are significantly faster than MovingLeastSquaresTransforms. The results are similar but not identical, so make sure to use the one that works best for your use case.

PARAMETER	DESCRIPTION
`landmarks_source`	`Source landmarks as x/y/z coordinates.` TYPE: `(M, 3) numpy array`
`landmarks_target`	`Target landmarks as x/y/z coordinates.` TYPE: `(M, 3) numpy array`
`batch_size`	`Batch size for transforming points. The thin-plate spline generating a (N, M) distance matrix, where N is the number of points and M is the number of source landmarks. Because this can get prohibitively expensive, we're batching the transformation by default. Please note that the the overhead from batching seems negligible. Ignored on the "fastcore" backend, which never materialises that matrix in the first place.` TYPE: `int` DEFAULT: `100000`
`backend`	Which implementation to use. `None` (default) defers to `navis.config.default_transform_backend`. "python" is deprecated and will be removed in 3.0. TYPE: `"auto" \| "fastcore" \| "python"` DEFAULT: `None`

Examples:

>>> from navis import transforms
>>> import numpy as np
>>> # Generate some mock landmarks
>>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
>>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
>>> tr = transforms.thinplate.TPStransform(src, trg)
>>> points = np.array([[0, 0, 0], [50, 50, 50]])
>>> tr.xform(points)
array([[ 1.        , 15.        ,  5.        ],
       [40.55555556, 54.        , 65.        ]])

Source code in navis/transforms/thinplate.py

class TPStransform(BackendMixin, BaseTransform):
    """Thin Plate Spline transforms of 3D spatial data.

    Runs on navis-fastcore's Rust implementation; the deprecated "python"
    backend runs on `morphops` instead (see the `backend` parameter). The spline
    is always *fitted* with `morphops` - numpy's LAPACK-backed solve is faster
    than fastcore's here - so the backend only changes how points are
    *transformed*, which is where fastcore is ~10-15x faster. The two agree to
    ~1e-13.

    Notes
    -----
    At least in my hands, `TPStransforms` are significantly faster than
    `MovingLeastSquaresTransforms`. The results are similar but not identical,
    so make sure to use the one that works best for your use case.

    Parameters
    ----------
    landmarks_source :  (M, 3) numpy array
                        Source landmarks as x/y/z coordinates.
    landmarks_target :  (M, 3) numpy array
                        Target landmarks as x/y/z coordinates.
    batch_size :        int, optional
                        Batch size for transforming points. The
                        thin-plate spline generating a (N, M) distance
                        matrix, where N is the number of points and M
                        is the number of source landmarks. Because
                        this can get prohibitively expensive, we're
                        batching the transformation by default.
                        Please note that the the overhead from batching
                        seems negligible.
                        Ignored on the "fastcore" backend, which never
                        materialises that matrix in the first place.
    backend :           "auto" | "fastcore" | "python", optional
                        Which implementation to use. `None` (default) defers to
                        `navis.config.default_transform_backend`. "python" is
                        deprecated and will be removed in 3.0.

    Examples
    --------
    >>> from navis import transforms
    >>> import numpy as np
    >>> # Generate some mock landmarks
    >>> src = np.array([[0, 0, 0], [10, 10, 10], [100, 100, 100], [80, 10, 30]])
    >>> trg = np.array([[1, 15, 5], [9, 18, 21], [80, 99, 120], [5, 10, 80]])
    >>> tr = transforms.thinplate.TPStransform(src, trg)
    >>> points = np.array([[0, 0, 0], [50, 50, 50]])
    >>> tr.xform(points)
    array([[ 1.        , 15.        ,  5.        ],
           [40.55555556, 54.        , 65.        ]])

    """

    _fallback_backend = "python"

    def __init__(
        self,
        landmarks_source: np.ndarray,
        landmarks_target: np.ndarray,
        batch_size: int = 100_000,
        backend: str = None,
    ):
        """Initialize class."""
        self.batch_size = batch_size
        self.source = np.asarray(landmarks_source)
        self.target = np.asarray(landmarks_target)
        self._backend = backend

        # Some checks
        if self.source.shape[1] != 3:
            raise ValueError(f"Expected (N, 3) array, got {self.source.shape}")
        if self.target.shape[1] != 3:
            raise ValueError(f"Expected (N, 3) array, got {self.target.shape}")

        if self.source.shape[0] != self.target.shape[0]:
            raise ValueError(
                "Number of source landmarks must match number of target landmarks."
            )

        self._W, self._A = None, None
        self._fc = None

    def __eq__(self, other) -> bool:
        """Implement equality comparison."""
        if isinstance(other, TPStransform):
            if self.source.shape[0] == other.source.shape[0]:
                if np.all(self.source == other.source):
                    if np.all(self.target == other.target):
                        return True
        return False

    def __neg__(self) -> "TPStransform":
        """Invert direction."""
        # Switch source and target
        return TPStransform(
            self.target,
            self.source,
            batch_size=self.batch_size,
            backend=self._backend,
        )

    @property
    def _fastcore_tps(self):
        """The (lazily built) fastcore transform used to apply the spline.

        Built from morphops-computed coefficients via `from_coefs` rather than
        by letting fastcore fit the spline itself. The two halves of a TPS scale
        oppositely: fastcore's fused-distance `xform` is ~10-15x faster than the
        morphops one, but its fit (a blocked LU) is several times slower than
        numpy's LAPACK-backed `solve`. Fitting with morphops and applying with
        fastcore takes the faster of each - so on this backend the fit is never
        a regression and the transform is much quicker.

        Deferred for the same reason the coefficients are: template packages
        such as `flybrains` construct every transform they ship at import time,
        and most are never used.

        """
        if self._fc is None:
            from .. import utils

            # `self.W`/`self.A` trigger the (morphops) fit if it hasn't happened
            self._fc = utils.fastcore.TpsTransform.from_coefs(
                self.source, self.W, self.A
            )
        return self._fc

    def _calc_tps_coefs(self):
        # Always fit with morphops, on either backend: numpy's LAPACK-backed
        # solve beats fastcore's blocked LU and the two agree to ~1e-13. On the
        # fastcore backend the coefficients are then handed to
        # `TpsTransform.from_coefs` (see `_fastcore_tps`), so the fast xform is
        # applied to a numpy-fitted spline.
        self._W, self._A = _morphops().tps_coefs(self.source, self.target)

    @property
    def W(self):
        if isinstance(self._W, type(None)):
            # Calculate coefficients
            self._calc_tps_coefs()
        return self._W

    @property
    def A(self):
        if isinstance(self._A, type(None)):
            # Calculate coefficients
            self._calc_tps_coefs()
        return self._A

    @property
    def matrix_affine(self):
        """Return the affine transformation matrix."""
        # The first row in self.A is the translation vector
        # The next 3x3 block is the rotation matrix
        # Let's combine these into a typical 4x4 transformation matrix
        # where the last row is [0, 0, 0, 1]
        m = np.zeros((4, 4))
        m[0:3, 0:3] = self.A[1:4, :].T
        m[0:3, 3] = self.A[0, :]
        m[3] = [0, 0, 0, 1]
        return m

    def copy(self):
        """Make copy."""
        x = TPStransform(self.source, self.target)

        # N.B. this also carries over `_backend` and any already-computed
        # coefficients (`_W`/`_A`/`_fc`). Sharing the latter is safe - they are
        # never mutated - and saves re-fitting, which is cubic in the number of
        # landmarks.
        x.__dict__.update(self.__dict__)

        return x

    def xform(self, points: np.ndarray) -> np.ndarray:
        """Transform points.

        Parameters
        ----------
        points :    (N, 3) array
                    Points to transform.

        Returns
        -------
        pointsxf :  (N, 3) array
                    Transformed points.

        """
        if isinstance(points, pd.DataFrame):
            if any(c not in points for c in ["x", "y", "z"]):
                raise ValueError("DataFrame must have x/y/z columns.")
            points = points[["x", "y", "z"]].values

        if self.backend == "fastcore":
            # No batching here: fastcore fuses the (N, M) distance matrix into
            # the accumulation instead of building it, so peak memory is just
            # the output array.
            return self._fastcore_tps.xform(np.asarray(points))

        mops = _morphops()
        batch_size = self.batch_size if self.batch_size else points.shape[0]
        points_xf = []
        for i in range(0, points.shape[0], batch_size):
            # Get the current batch of points
            batch = points[i : i + batch_size]

            # N.B. U is of shape (N, M) where N is the number of points and M is the
            # number of source landmarks. This can get fairly expensive
            # (which is precisely why we batch the transformation)!
            U = mops.K_matrix(batch, self.source)
            P = mops.P_matrix(batch)
            # The warped pts are the affine part + the non-uniform part
            points_xf.append(np.matmul(P, self.A) + np.matmul(U, self.W))

        # Concatenate all batches
        return np.concatenate(points_xf, axis=0)

`matrix_affine` `property` #

Return the affine transformation matrix.

`init` #

Initialize class.

Source code in navis/transforms/thinplate.py

def __init__(
    self,
    landmarks_source: np.ndarray,
    landmarks_target: np.ndarray,
    batch_size: int = 100_000,
    backend: str = None,
):
    """Initialize class."""
    self.batch_size = batch_size
    self.source = np.asarray(landmarks_source)
    self.target = np.asarray(landmarks_target)
    self._backend = backend

    # Some checks
    if self.source.shape[1] != 3:
        raise ValueError(f"Expected (N, 3) array, got {self.source.shape}")
    if self.target.shape[1] != 3:
        raise ValueError(f"Expected (N, 3) array, got {self.target.shape}")

    if self.source.shape[0] != self.target.shape[0]:
        raise ValueError(
            "Number of source landmarks must match number of target landmarks."
        )

    self._W, self._A = None, None
    self._fc = None

`copy` #

Make copy.

Source code in navis/transforms/thinplate.py

def copy(self):
    """Make copy."""
    x = TPStransform(self.source, self.target)

    # N.B. this also carries over `_backend` and any already-computed
    # coefficients (`_W`/`_A`/`_fc`). Sharing the latter is safe - they are
    # never mutated - and saves re-fitting, which is cubic in the number of
    # landmarks.
    x.__dict__.update(self.__dict__)

    return x

`xform` #

Transform points.

PARAMETER	DESCRIPTION
`points`	`Points to transform.` TYPE: `(N, 3) array`

RETURNS	DESCRIPTION
`pointsxf`	Transformed points. TYPE: `(N, 3) array`

Source code in navis/transforms/thinplate.py

def xform(self, points: np.ndarray) -> np.ndarray:
    """Transform points.

    Parameters
    ----------
    points :    (N, 3) array
                Points to transform.

    Returns
    -------
    pointsxf :  (N, 3) array
                Transformed points.

    """
    if isinstance(points, pd.DataFrame):
        if any(c not in points for c in ["x", "y", "z"]):
            raise ValueError("DataFrame must have x/y/z columns.")
        points = points[["x", "y", "z"]].values

    if self.backend == "fastcore":
        # No batching here: fastcore fuses the (N, M) distance matrix into
        # the accumulation instead of building it, so peak memory is just
        # the output array.
        return self._fastcore_tps.xform(np.asarray(points))

    mops = _morphops()
    batch_size = self.batch_size if self.batch_size else points.shape[0]
    points_xf = []
    for i in range(0, points.shape[0], batch_size):
        # Get the current batch of points
        batch = points[i : i + batch_size]

        # N.B. U is of shape (N, M) where N is the number of points and M is the
        # number of source landmarks. This can get fairly expensive
        # (which is precisely why we batch the transformation)!
        U = mops.K_matrix(batch, self.source)
        P = mops.P_matrix(batch)
        # The warped pts are the affine part + the non-uniform part
        points_xf.append(np.matmul(P, self.A) + np.matmul(U, self.W))

    # Concatenate all batches
    return np.concatenate(points_xf, axis=0)

`navis.transforms.clear_transform_cache` #

Drop all cached (parsed) fastcore transforms.

Mostly useful if a registration changed on disk - we do not stat the files on every call.

Source code in navis/transforms/backends.py

def clear_transform_cache():
    """Drop all cached (parsed) fastcore transforms.

    Mostly useful if a registration changed on disk - we do not stat the files
    on every call.

    """
    get_cmtk_reg.cache_clear()
    get_elastix_transform.cache_clear()
    elastix_is_invertible.cache_clear()

`navis.transforms.set_transform_backend` #

Set the backend used for point transforms.

Applies to CMTK and elastix transforms, and to the landmark transforms (navis.transforms.TPStransform and navis.transforms.MovingLeastSquaresTransform).

PARAMETER DESCRIPTION

backend

    - "auto" (default) and "fastcore" both use navis-fastcore's
      in-process Rust implementation
    - "binary"/"python" force the original implementation: the
      external binaries (`streamxform`, `transformix`) for
      CMTK/elastix, `morphops`/`molesq` for the landmark
      transforms. The two names are interchangeable - each
      transform reports whichever fits it. **Both are deprecated
      and will be removed in 3.0**; selecting one warns.

TYPE: "auto" | "fastcore" | "binary" | "python" DEFAULT: None

elastix_invertible

    Whether the bridging graph may traverse an elastix registration
    backwards (only possible on the fastcore backend). Off by
    default - see `navis.config.elastix_invertible`. Note this
    re-routes a substantial fraction of the bridging paths.

TYPE: bool DEFAULT: None

Examples:

>>> import navis
>>> navis.transforms.set_transform_backend('fastcore')
>>> navis.transforms.set_transform_backend('auto')

Source code in navis/transforms/backends.py

def set_transform_backend(backend: str = None, elastix_invertible: bool = None):
    """Set the backend used for point transforms.

    Applies to CMTK and elastix transforms, and to the landmark transforms
    ([`navis.transforms.TPStransform`][] and
    [`navis.transforms.MovingLeastSquaresTransform`][]).

    Parameters
    ----------
    backend :   "auto" | "fastcore" | "binary" | "python", optional
                - "auto" (default) and "fastcore" both use navis-fastcore's
                  in-process Rust implementation
                - "binary"/"python" force the original implementation: the
                  external binaries (`streamxform`, `transformix`) for
                  CMTK/elastix, `morphops`/`molesq` for the landmark
                  transforms. The two names are interchangeable - each
                  transform reports whichever fits it. **Both are deprecated
                  and will be removed in 3.0**; selecting one warns.
    elastix_invertible : bool, optional
                Whether the bridging graph may traverse an elastix registration
                backwards (only possible on the fastcore backend). Off by
                default - see `navis.config.elastix_invertible`. Note this
                re-routes a substantial fraction of the bridging paths.

    Examples
    --------
    >>> import navis
    >>> navis.transforms.set_transform_backend('fastcore')
    >>> navis.transforms.set_transform_backend('auto')

    """
    if backend is not None:
        if backend not in BACKENDS:
            raise ValueError(
                f'Unknown transform backend "{backend}". Must be one of {BACKENDS}.'
            )
        # Warn *here*, where the choice is actually made: `stacklevel=2` puts the
        # blame on the caller's frame, which is what Python's default filters need
        # to show it at all.
        warn_if_deprecated_backend(backend, stacklevel=3)
        config.default_transform_backend = backend

    if elastix_invertible is not None:
        config.elastix_invertible = bool(elastix_invertible)

    # Both of the above feed the transforms' `invertible` property, which in turn
    # decides the shape of the (cached) bridging graph - so both caches go.
    clear_transform_cache()

    from .templates import registry

    registry.clear_caches()

transforms

navis.transforms.AffineTransform #

__init__ #

copy #

xform #

navis.transforms.AliasTransform #

__init__ #

copy #

xform #

navis.transforms.CMTKtransform #

regargs: list property #

append #

check_if_possible #

copy #

from_file staticmethod #

make_args #

parse_cmtk_output #

to_dfield #

to_grid_transform #

xform #

xform_image #

navis.transforms.ElastixTransform #

can_invert: bool property #

invertible: bool property #

check_if_possible #

copy #

read_output_file #

write_input_file #

xform #

navis.transforms.FunctionTransform #

__init__ #

copy #

xform #

navis.transforms.GridTransform #

affine: AffineTransform property #

shape: tuple property #

__init__ #

calculate_affine #

copy #

from_file classmethod #

from_warpfield classmethod #

xform #

navis.transforms.H5transform #

quantization_multiplier property #

shape property writable #

spacing property #

use_cache property writable #

__init__ #

copy #

from_file staticmethod #

full_ingest #

precache #

xform #

xform_image #

navis.transforms.MovingLeastSquaresTransform #

matrix_affine property #

transformer property #

__init__ #

copy #

xform #

navis.transforms.SimilarityTransform #

n: int property #

rotation: np.ndarray property #

scale_factor: float property #

translation: np.ndarray property #

copy #

update #

xform #

navis.transforms.TPStransform #

matrix_affine property #

__init__ #

copy #

xform #

navis.transforms.clear_transform_cache #

navis.transforms.set_transform_backend #

`navis.transforms.AffineTransform` #

`init` #

`copy` #

`xform` #

`navis.transforms.AliasTransform` #

`init` #

`copy` #

`xform` #

`navis.transforms.CMTKtransform` #

`regargs: list` `property` #

`append` #

`check_if_possible` #

`copy` #

`from_file` `staticmethod` #

`make_args` #

`parse_cmtk_output` #

`to_dfield` #

`to_grid_transform` #

`xform` #

`xform_image` #

`navis.transforms.ElastixTransform` #

`can_invert: bool` `property` #

`invertible: bool` `property` #

`check_if_possible` #

`copy` #

`read_output_file` #

`write_input_file` #

`xform` #

`navis.transforms.FunctionTransform` #

`init` #

`copy` #

`xform` #

`navis.transforms.GridTransform` #

`affine: AffineTransform` `property` #

`shape: tuple` `property` #

`init` #

`calculate_affine` #

`copy` #

`from_file` `classmethod` #

`from_warpfield` `classmethod` #

`xform` #

`navis.transforms.H5transform` #

`quantization_multiplier` `property` #

`shape` `property` `writable` #

`spacing` `property` #

`use_cache` `property` `writable` #

`init` #

`copy` #

`from_file` `staticmethod` #

`full_ingest` #

`precache` #

`xform` #

`xform_image` #

`navis.transforms.MovingLeastSquaresTransform` #

`matrix_affine` `property` #

`transformer` `property` #

`init` #

`copy` #

`xform` #

`navis.transforms.SimilarityTransform` #

`n: int` `property` #

`rotation: np.ndarray` `property` #

`scale_factor: float` `property` #

`translation: np.ndarray` `property` #

`copy` #

`update` #

`xform` #

`navis.transforms.TPStransform` #

`matrix_affine` `property` #

`init` #

`copy` #

`xform` #

`navis.transforms.clear_transform_cache` #

`navis.transforms.set_transform_backend` #