Functions

Tools for aligning spatial transcriptomics data

The LDDMM fits (LDDMM, LDDMM_3D_to_slice, and align_cells for two sets of cells), rasterization and the transforms of fitted mappings are implemented in private modules next to this file (_lddmm*.py, _cells.py, _raster.py, _transforms.py) and exposed here. The LDDMM fits use a compiled C++ engine when it is installed and pure torch otherwise; backend_info() reports which.

Functions:

LDDMM(xI, I, xJ, J[, pointsI, pointsJ, L, ...])

Run LDDMM between a pair of images.

LDDMM_3D_to_slice(xI, I, xJ, J[, pointsI, ...])

LDDMM for 3D to 2D slice mapping.

L_T_from_points(pointsI, pointsJ)

Compute an affine transformation from points.

align_cells(xI, yI, xJ, yJ[, dx, blur_um, ...])

Align two sets of cells with LDDMM on density images, coarse to fine.

analyze3Dalign(labelfile, xv, v, mat, xJ, ...)

Create dataframe with region annotations and 3D coordinates of target brain slice

backend_info()

Which LDDMM backends are available.

build_transform(xv, v, A[, direction, XJ, ...])

Create sample points to transform source to target from affine and velocity.

build_transform3D(xv, v, A[, direction, XJ, ...])

Create sample points to transform source to target from affine and velocity.

calculate_tre(pointsI, pointsJ)

Target registration error (TRE) between two sets of corresponding points

clip(I)

clip an arrays values between 0 and 1.

download_aba_image_labels(imageurl, ...)

Create 3D altas image and region annotations

download_aba_ontology(url, file_name)

Create 3D altas ontology

extent_from_x(xJ)

Given a set of pixel locations, returns an extent 4-tuple for use with np.imshow.

interp(x, I, phii, **kwargs)

Interpolate the 2D image I, with regular grid positions stored in x (1d arrays), at the positions stored in phii (3D arrays with first channel storing component)

interp3D(x, I, phii, **kwargs)

Interpolate the 3D image I, with regular grid positions stored in x (1d arrays), at the positions stored in phii (4D arrays with first channel storing component)

make_scree(W, name[, p])

Create a scree plot for a given set of eigenvalues.

normalize(arr[, t_min, t_max])

Linearly normalizes an array between two specifed values.

plot_brain_regions(df)

Plot all brain regions in target brain slice with different color.

plot_subset_brain_regions(df, brain_regions)

Plot subset of brain regions in target brain slice with different color.

rasterize(x, y[, g, dx, blur, expand, draw, ...])

Rasterize a spatial transcriptomics dataset into a density image

rasterizePCA(x, y, G)

Rasterize a spatial transcriptomics dataset into a density image and perform PCA on it.

rasterize_with_signal(x, y[, s, dx, blur, ...])

Rasterize a spatial transcriptomics dataset with a signal per cell

saveRasters(X, Y, W, V, Z, scree_fig, nrows, ...)

Create an RGB image principal components 1-3 and 4-6.

to_A(L, T)

Convert a linear transform matrix and a translation vector into an affine matrix.

to_A_3D(L, T)

Convert a linear transform matrix and a translation vector into an affine matrix.

transform_image_source_to_target(xv, v, A, xI, I)

Transform an image from source to target space, by sampling it at the backward map.

transform_image_source_with_A(A, XI, I, XJ)

Transform an image with an affine matrix

transform_image_target_to_source(xv, v, A, xJ, J)

Transform an image from target to source space, by sampling it at the forward map.

transform_points_source_to_target(xv, v, A, ...)

Transform points from source to target space.

transform_points_target_to_source(xv, v, A, ...)

Transform points from target to source space.

v_to_phii(xv, v)

Integrate a velocity field over time to return a position field (diffeomorphism).

v_to_phii_3D(xv, v)

Integrate a 3D velocity field over time to return a 3D position field (diffeomorphism).
STalign.STalign.LDDMM(xI, I, xJ, J, pointsI=None, pointsJ=None, L=None, T=None, A=None, v=None, xv=None, a=500.0, p=2.0, expand=2.0, nt=3, niter=None, diffeo_start=0, epL=None, epT=None, epV=None, sigmaM=1.0, sigmaB=2.0, sigmaA=5.0, sigmaR=None, sigmaP=20.0, device='cpu', dtype=torch.float64, muB=None, muA=None, *, draw='end', num_threads=None, backend='auto', affine_step_control=False, mixture=None, prefit='auto', sigma_disp=94.03, legacy_defaults=False)

Run LDDMM between a pair of images.

This jointly estimates an affine transform A, and a diffeomorphism phi. The map is of the form x -> A phi x

By default the fit starts with a coarse affine pre-alignment (prefit), and the regularization weight sigmaR and the step sizes epL, epT, epV are chosen from the images unless given (see their descriptions). Coordinates are assumed to be in micrometers. legacy_defaults=True turns these defaults off.

Parameters:

  • xI (list of torch tensor) – Location of voxels in source image I

  • I (torch tensor) – source image I, with channels along first axis

  • xJ (list of torch tensor) – Location of voxels in target image J

  • J (torch tensor) – Target image J, with channels along first axis

  • L (torch tensor) – Initial guess for linear transform (2x2 torch tensor). Defaults to None (identity).

  • T (torch tensor) – Initial guess for translation (2 element torch tensor). Defaults to None (identity)

  • A (torch tensor) – Initial guess for affine matrix. Either L and T can be specified, or A, but not both. Defaults to None (identity).

  • v (torch tensor) – Initial guess for velocity field, shape (nt, len(xv[0]), len(xv[1]), 2)

  • xv (torch tensor) – pixel locations for velocity field

  • a (float) – Smoothness scale of velocity field (default 500.0)

  • p (float) – Power of Laplacian in velocity regularization (default 2.0)

  • expand (float) – Factor to expand size of velocity field around image boundaries (default 2.0)

  • nt (int) – Number of timesteps for integrating velocity field (default 3). Ignored if you input v.

  • pointsI (torch tensor) – N x 2 set of corresponding points for matching in source image. Default None (no points).

  • pointsJ (torch tensor) – N x 2 set of corresponding points for matching in target image. Default None (no points).

  • niter (int) – Number of iterations of gradient descent optimization. Default None: 2000 (5000 with legacy_defaults=True).

  • diffeo_start (int) – Number of iterations of gradient descent optimization for affine only, before nonlinear deformation.

  • epL (float) – Gradient descent step size for linear part of affine.

  • epT (float) – Gradient descent step size of translation part of affine.

  • epV (float) – Gradient descent step size for velocity field. Default None for each of epL, epT, epV: the steps of earlier versions (2e-8, 2e-1, 2e3) times sigmaM^2 * 1.1 / G, where G estimates the curvature of the image matching term from the target’s variance and the source’s gradients under the initial affine (after the pre-fit). This adapts the steps to the image intensities and pixel size, where the old defaults could be far too large or too small. Each is then lowered where needed to keep gradient descent stable, from the curvature of the energy at the initial affine (the image term with the same assumptions as G, the landmarks and the regularizer): epV to at most 1 / (lambda_v + DV / sigmaR^2), where lambda_v is the largest eigenvalue of the kernel-smoothed curvature of the image and landmark terms in v (by subspace iteration) and DV the area of a velocity grid cell; epL and epT by a common factor, so that the largest eigenvalue of the affine curvature scaled by the steps (divided by 10, as during the deformation, unless diffeo_start > 0) is at most 0.5. This matters for small images (G grows with the tissue area, the stable epV does not), landmarks, a small sigma_disp and coordinates far from the origin. A given value is used as is.

  • sigmaM (float) – Standard deviation of image matching term for Gaussian mixture modeling in cost function. This term generally controls matching accuracy with smaller corresponding to more accurate. As an common example (rule of thumb), you could chose this parameter to be the variance of the pixels in your target image.

  • sigmaB (float) – Standard deviation of backtround term for Gaussian mixture modeling in cost function. If there is missing tissue in target, we may label some pixels in target as background, and not enforce matching here.

  • sigmaA (float) – Standard deviation of artifact term for Gaussian mixture modeling in cost function. If there are artifacts in target or other lack of corresponding between template and target, we may label some pixels in target as artifact, and not enforce matching here.

  • sigmaR (float) – Standard deviation for regularization. Smaller sigmaR means a smoother resulting transformation. Regularization is of the form 0.5/sigmaR^2 sum_{t=0}^{nt-1} int_X |L v_t|^2 dx with L = (1 - a^2 Laplacian)^p, i.e. a sum over the nt time steps without a dt = 1/nt factor (nt times 0.5/sigmaR^2 int_0^1 int_X |Lv|^2 dx dt). Changing nt at a fixed sigmaR therefore changes the strength of the regularization. Default None: computed from sigma_disp as sigma_disp * sqrt(12 pi) * a * sqrt(nt) * (dx / 30) * sigmaM / s_J, where dx is the target pixel size in source units (under the affine after the pre-fit) and s_J the typical tissue intensity of the target. The prior then allows displacements of about sigma_disp micrometers whatever the pixel size, intensity scale, a and nt (at the defaults a = 500, nt = 3, a 30 um target with s_J = sigmaM gets 5e5, the previous default). A given value is used as is, as in earlier versions (5e5 with legacy_defaults=True).

  • sigmaP (float) – Standard deviation for matching of points. Cost is of the form 0.5/sigmaP^2 sum_i (source_point_i - target_point_i)^2, where the source points are mapped to the target by the exact inverse of the fitted target-to-source map, i.e. as transform_points_source_to_target maps them after the fit. (Earlier versions used forward Euler steps, which invert that map only approximately; legacy_defaults=True uses them.)

  • device (str) – torch device. defaults to ‘cpu’. Can also be ‘cuda:0’ for example.

  • dtype (torch dtype) – torch data type. defaults to torch.float64

  • muA (torch tensor whose dimension is the same as the target image) – Defaults to None, which means we estimate this. If you provide a value, we will not estimate it. Give one value per channel of the target (size 3 for a RGB image) or a single value for all channels. Ints, floats, lists, numpy arrays and torch tensors are accepted.

  • muB (torch tensor whose dimension is the same as the target image) – Defaults to None, which means we estimate this. If you provide a value, we will not estimate it. Same format as muA.

  • draw ('end', bool or int, keyword only) – Progress figures (images, weights, velocity and energies). ‘end’ (default) draws them once after the last iteration, True every 10 iterations, an int N every N iterations, and False, None or 0 never. Drawing is slow; a drawing error only issues a warning.

  • num_threads (int, keyword only) – Maximum number of CPU threads. Defaults to the environment variable STALIGN_NUM_THREADS, else torch.get_num_threads(). Fewer threads are used for small images.

  • backend ({'auto', 'cpp', 'torch'}, keyword only) – ‘auto’ (default) uses the compiled C++ engine when it is installed and the problem runs on the CPU in float64, and the pure torch engine otherwise. ‘cpp’ raises an error if the compiled engine cannot run the problem; ‘torch’ always uses torch. Both give the same results up to rounding.

  • affine_step_control (bool, keyword only) – If True, the step of each affine parameter (each entry of L and T) is scaled by its own multiplier between 1e-3 and 1. The multiplier is halved whenever that parameter’s gradient changes sign from one iteration to the next (the previous step overshot) and grows by 5% per iteration (up to 1) while the sign stays the same, so steps are never larger than epL and epT. This stops the affine part from oscillating when epL or epT is too large. Default False, which gives the plain gradient steps of earlier versions. (LDDMM_3D_to_slice uses it by default.)

  • mixture (None, 'em' or dict, keyword only) – How the weights WM, WA, WB of the matching, artifact and background classes are found. None (default): as in earlier versions (WM starts at 0.5; from iteration 50, every 5 iterations, an E-step with the fixed variances sigmaM^2, sigmaA^2, sigmaB^2). ‘em’: a mixture EM whose matching variance vM is estimated from the residual. WM starts at 1. From iteration 50, every 5 iterations, each target pixel is weighted between the matching class (variance vM) and the background class (mean muB, the median of the target’s border unless muB is given; variance (sigmaB/sigmaM)^2 vM), with the matching prior kept >= 0.5. Pixels that look like background while the deformed source does not match them, as where the target section is incomplete, then get a low matching weight. A dict selects ‘em’ and overrides some of its options: start (50), every (5), classes (‘MB’; ‘MAB’ adds the artifact class with mean muA, the mean of the target unless muA is given, and variance (sigmaA/sigmaM)^2 vM; ‘M’ and ‘MA’ are also accepted), pi_min (0.5), pi_init ((0.9, 0.05, 0.05)), fix_muA (False), fix_muB (True: muB is not re-estimated). Only the weights change; the matching term keeps sigmaM. Worth trying where the target lacks tissue that the source has; it can make targets of uneven density worse. align_cells uses it by default.

  • prefit ('auto', True or False, keyword only) – Coarse affine pre-alignment before LDDMM. The initial affine (L, T or A; identity by default) is rotated about the target tissue’s centroid in 12 steps of 30 degrees, each start is fitted as an affine map on images pooled to about 80 pixels across, and the best one is refined at twice that resolution; starts whose energies are within 1% of the best count as ties, and the smallest rotation among them wins. LDDMM then starts from the result. With landmarks the pre-fit includes the landmark term and keeps the initial orientation. ‘auto’ (default) runs it unless landmarks (pointsI, pointsJ) or a nonzero initial velocity field v are given; True always; False never. It is deterministic, and what it did is returned as ‘prefit’. If the target has no tissue (a constant image) or the pre-fit fails, a warning is issued and the fit starts from the given affine.

  • sigma_disp (float, keyword only) – Prior displacement scale in micrometers (default 94.03, which reproduces the previous default sigmaR=5e5 on 30 um images). Used only when sigmaR is None; larger values allow larger deformations. It assumes that coordinates are in micrometers.

  • legacy_defaults (bool, keyword only) – If True, use the defaults of version 1.0 instead of those chosen from the data (see CHANGES.md): no pre-fit (unless prefit=True), niter=5000, epL=2e-8, epT=2e-1, epV=2e3, sigmaR=5e5, and the forward Euler map in the landmark term. The corrections listed in CHANGES.md still apply, so results can differ from those of version 1.0. Default False.

  • dictionary (Returns a)

  • -------

  • {

  • 'A' (torch tensor) – Affine transform

  • 'v' (torch tensor) – Velocity field

  • 'xv' (list of torch tensor) – Pixel locations in v

  • 'WM' (torch tensor) – Resulting weight 2D array (matching)

  • 'WB' (torch tensor) – Resulting weight 2D array (background)

  • 'WA' (torch tensor) – Resulting weight 2D array (artifact)

  • 'prefit' (dict) – What the pre-fit did: ‘ran’ (bool) and ‘A’ (the affine LDDMM started from: the pre-fit’s result, or the initial affine if it did not run); if it ran, ‘A_start’ (the initial affine), ‘energy’ (the final energy of every start), ‘starts’ (their rotations in degrees), ‘best’ (the index of the winner), ‘pool’ (the source and target pooling factors) and ‘seconds’; if not, ‘reason’.

  • 'settings' (dict) – The niter, sigmaR, epL, epT and epV the fit used. Passing them with A=prefit[‘A’] (in place of L, T or A) and prefit=False (and legacy_defaults=True if it was given), the other arguments unchanged, repeats the fit exactly, without the data-driven choices.

  • }

  • pass (A and v are the pair used by the last iteration's forward)

  • build_transform(xv (so)

  • v

  • A)

  • returned. (reproduces the final alignment. With niter=0 the initial state (after the pre-fit) is)

STalign.STalign.LDDMM_3D_to_slice(xI, I, xJ, J, pointsI=None, pointsJ=None, L=None, T=None, A=None, v=None, xv=None, a=500.0, p=2.0, expand=1.25, nt=3, niter=5000, diffeo_start=0, epL=1e-06, epT=10.0, epV=None, sigmaM=1.0, sigmaB=2.0, sigmaA=5.0, sigmaR=None, sigmaP=20.0, device='cpu', dtype=torch.float64, muA=None, muB=None, *, draw='end', num_threads=None, backend='auto', affine_step_control=True, mixture=None, prefit='auto', sigma_disp=94.03, legacy_defaults=False)

LDDMM for 3D to 2D slice mapping.

Jointly estimates a 3D affine transform A and a 3D diffeomorphism phi that map the 3D source I (for example an atlas) onto the 2D target slice J. The map is of the form x -> A phi x. The parameters are those of LDDMM, with the differences below. By default the fit starts with a coarse pre-alignment of the slice plane (prefit) and sigmaR is chosen from the images; coordinates are assumed to be in micrometers.

Parameters:

  • xI (list of 3 arrays) – Location of voxels in the source volume I, with channels along the first axis.

  • xJ (list of 2 or 3 arrays) – Location of pixels in the target slice J (row, column). With two arrays the slice lies in the plane z = 0 of the target space.

  • J (array) – Target slice with channels along the first axis, shape (C, H, W) or (C, 1, H, W).

  • L (arrays) – Initial 3x3 linear part and 3 element translation, or 4x4 affine A (not both).

  • T (arrays) – Initial 3x3 linear part and 3 element translation, or 4x4 affine A (not both).

  • A (arrays) – Initial 3x3 linear part and 3 element translation, or 4x4 affine A (not both).

  • v (arrays) – Initial velocity field (nt, len(xv[0]), len(xv[1]), len(xv[2]), 3) and its grid.

  • xv (arrays) – Initial velocity field (nt, len(xv[0]), len(xv[1]), len(xv[2]), 3) and its grid.

  • expand (float) – Factor to expand size of velocity field around image boundaries (default 1.25)

  • sigmaR (float) – Standard deviation for regularization. Smaller sigmaR means a smoother transformation. Regularization is of the form 0.5/sigmaR^2 sum_{t=0}^{nt-1} int_X |L v_t|^2 dx with L = (1 - a^2 Laplacian)^p over all three axes, a sum over time steps without dt as in LDDMM. Earlier versions took its Fourier transform over only two of the three axes, which made it much stronger (and uneven), with the default 1e8. Default None: computed from sigma_disp as sigma_disp * sqrt(64 pi) * a^1.5 * sqrt(nt) * (dx / 30) * sigmaM / s_J (dx: the target pixel size in atlas units under the affine after the pre-fit; s_J: the typical tissue intensity of the target), so that the prior allows displacements of about sigma_disp micrometers (1.1e7 and 1.85e6 in the MERFISH and STARmap tutorials). A given value is used as is (2e6 with legacy_defaults=True, which gives deformations of about the size of the old regularizer’s). See CHANGES.md.

  • niter – Defaults 5000, 1e-6, 1e1 and None: 1e3 (no data-driven steps in 3D), lowered to sigmaR^2 / DV where needed to keep the velocity step stable for the regularizer (epV * DV / sigmaR^2 <= 1, DV the volume of a velocity grid cell). That happens only for a very small sigma_disp (a few um) with sigmaM well below the target’s intensity scale. A given value is used as is.

  • epL – Defaults 5000, 1e-6, 1e1 and None: 1e3 (no data-driven steps in 3D), lowered to sigmaR^2 / DV where needed to keep the velocity step stable for the regularizer (epV * DV / sigmaR^2 <= 1, DV the volume of a velocity grid cell). That happens only for a very small sigma_disp (a few um) with sigmaM well below the target’s intensity scale. A given value is used as is.

  • epT – Defaults 5000, 1e-6, 1e1 and None: 1e3 (no data-driven steps in 3D), lowered to sigmaR^2 / DV where needed to keep the velocity step stable for the regularizer (epV * DV / sigmaR^2 <= 1, DV the volume of a velocity grid cell). That happens only for a very small sigma_disp (a few um) with sigmaM well below the target’s intensity scale. A given value is used as is.

  • epV – Defaults 5000, 1e-6, 1e1 and None: 1e3 (no data-driven steps in 3D), lowered to sigmaR^2 / DV where needed to keep the velocity step stable for the regularizer (epV * DV / sigmaR^2 <= 1, DV the volume of a velocity grid cell). That happens only for a very small sigma_disp (a few um) with sigmaM well below the target’s intensity scale. A given value is used as is.

  • pointsI – Not used: there is no landmark term in 3D, and passing points issues a warning. Use L_T_from_points to initialize L and T from landmarks instead.

  • pointsJ – Not used: there is no landmark term in 3D, and passing points issues a warning. Use L_T_from_points to initialize L and T from landmarks instead.

  • sigmaP – Not used: there is no landmark term in 3D, and passing points issues a warning. Use L_T_from_points to initialize L and T from landmarks instead.

  • muA (torch tensor whose dimension is the same as the target image) – Defaults to None, which means we estimate this. If you provide a value, we will not estimate it. Give one value per channel of the target or a single value. For a grayscale target this should be a single value.

  • muB (torch tensor whose dimension is the same as the target image) – Defaults to None, which means we estimate this. If you provide a value, we will not estimate it. Same format as muA.

  • draw (keyword only) – As in LDDMM.

  • num_threads (keyword only) – As in LDDMM.

  • backend (keyword only) – As in LDDMM.

  • affine_step_control (bool, keyword only) – As in LDDMM, but True by default: each affine parameter’s step is scaled by a multiplier between 1e-3 and 1 that is halved when its gradient changes sign and grows by 5% per iteration (up to 1) otherwise. With the default epT, plain gradient steps make the translation jump back and forth by tens of micrometers per iteration until the end, so the result depends on rounding: in the tutorials, changing the initial translation in its last digit moved the cells by 24-84 um on average without the control, and by about 1e-12 um with it. False gives the plain gradient steps of earlier versions.

  • mixture (None, keyword only) – The mixture EM of LDDMM (mixture=’em’) is not available in 3D; anything but None raises NotImplementedError.

  • prefit ('auto', True or False, keyword only) – Coarse pre-alignment of the slice plane before LDDMM. Starting from the initial affine, the plane is tilted about its two in-plane axes by -6, 0 and 6 degrees each and shifted along its normal by -300, 0 and 300 um (27 starts); each start fits a rotation, an in-plane affine map and a shift on images pooled to about 80 pixels across, and the best one is refined at twice that resolution (ties within 1% go to the start closest to the initial plane). ‘auto’ (default) runs it unless a nonzero initial velocity field v is given; True always; False never. It is deterministic, and what it did is returned as ‘prefit’. If the target has no tissue (a constant image) or the pre-fit fails, a warning is issued and the fit starts from the given affine.

  • sigma_disp (float, keyword only) – Prior displacement scale in micrometers (default 94.03), used when sigmaR is None.

  • legacy_defaults (bool, keyword only) – If True, use fixed defaults instead of those chosen from the data (see CHANGES.md): no pre-fit (unless prefit=True), sigmaR=2e6 and epV=1e3. The corrections still apply (the regularizer over all three axes, the affine step control and the others listed in CHANGES.md), so results differ from those of version 1.0, whose default sigmaR=1e8 regularizes much less with the corrected regularizer. Default False.

  • dictionary (Returns a)

  • -------

  • {

  • 'A' (torch tensor) – 4x4 affine transform

  • 'v' (torch tensor) – Velocity field

  • 'xv' (list of torch tensor) – Pixel locations in v

  • 'WM' (torch tensor) – Resulting weights (matching, background, artifact), shape (1, H, W)

  • 'WB' (torch tensor) – Resulting weights (matching, background, artifact), shape (1, H, W)

  • 'WA' (torch tensor) – Resulting weights (matching, background, artifact), shape (1, H, W)

  • 'Xs' (torch tensor) – Locations in the source space sampled for each target pixel, shape (1, H, W, 3)

  • 'prefit' (dict) – As in LDDMM; the pre-fit’s ‘starts’ are (tilt 1, tilt 2, offset) triples.

  • 'settings' (dict) – As in LDDMM; the pre-fit’s ‘starts’ are (tilt 1, tilt 2, offset) triples.

  • }

  • A

  • state (v and Xs come from the last iteration's forward pass. With niter=0 the initial)

  • returned. ((after the pre-fit) is)

STalign.STalign.L_T_from_points(pointsI, pointsJ)

Compute an affine transformation from points.

Note for an affine transformation (6dof) we need 3 points.

Outputs, L,T should be rconstructed blockwize like [L,T;0,0,1]

Parameters:

  • pointsI (array) – An Nx2 array of floating point numbers describing source points in ROW COL order (not xy)

  • pointsJ (array) – An Nx2 array of floating point numbers describing target points in ROW COL order (not xy)

Returns:

  • L (array) – A 2x2 linear transform array.

  • T (array) – A 2 element translation vector

STalign.STalign.align_cells(xI, yI, xJ, yJ, dx=15.0, blur_um=30.0, levels=((4, 500), (2, 500)), *, mixture='em', **lddmm_kwargs)

Align two sets of cells with LDDMM on density images, coarse to fine.

Both sets of cells are rasterized with rasterize at pixel size dx and a Gaussian kernel of standard deviation blur_um (each image covers its own cells). LDDMM is then run on these images mean-pooled by the factor of each level in turn, from coarse to fine. Each level starts from the previous level’s affine transform and velocity field, so the result is a single map from source to target that the transform functions use as LDDMM’s.

Parameters:

  • xI (numpy arrays of length N) – x and y locations of the source cells.

  • yI (numpy arrays of length N) – x and y locations of the source cells.

  • xJ (numpy arrays of length M) – x and y locations of the target cells, in the same units as the source cells (for example micrometers).

  • yJ (numpy arrays of length M) – x and y locations of the target cells, in the same units as the source cells (for example micrometers).

  • dx (float) – Pixel size of the density images, in the units of the cell locations (default 15.0).

  • blur_um (float) – Standard deviation of the rasterization kernel, in the units of the cell locations (default 30.0, i.e. 2 pixels of 15). See rasterize.

  • levels (sequence of (factor, niter) pairs) – The pyramid, coarsest level first: at each level LDDMM runs niter iterations on the images mean-pooled by the integer factor (the source by the factor that gives it about the same pixel size as the target under the current affine). The default ((4, 500), (2, 500)) runs 500 iterations at 4*dx = 60 and then 500 at 2*dx = 30.

  • mixture ('em', None or dict, keyword only) – LDDMM’s mixture argument, passed to every level. The default is ‘em’, unlike LDDMM’s (None): LDDMM’s mixture EM, which estimates the noise of the matching term from the residual and gives a low matching weight to target pixels that look like background where the deformed source has tissue, as where the target section lacks cells that the source has. mixture=None turns it off (LDDMM’s weights of earlier versions, with the fixed variances sigmaM^2, sigmaA^2 and sigmaB^2); a dict selects ‘em’ with some of its options changed (see LDDMM). The value you pass is used as given. On 64 synthetic cell-to-cell test problems (see CHANGES.md), the default cut the geometric mean error of the 8 partial sections from 15.2 to 4.3 um (all 8 better), but raised that of the 6 sections of uneven cell density from 7.9 to 12.8 um (by up to 99%) and made 5 of 8 serial heart sections 10-57% worse; the other 42 changed by less than 2%. Pass mixture=None where the target section is complete.

  • **lddmm_kwargs –

    Keyword arguments of LDDMM, for example pointsI and pointsJ, sigmaM, a, num_threads or backend. They apply at every level, except as noted below. Everything else takes LDDMM’s defaults at each level, so the parameters that LDDMM derives from the images when they are not given are derived at each level from its pooled images.

    • niter: set by levels (passing niter is an error).

    • L, T or A, v and xv (the initial state), prefit and diffeo_start: apply to the first level, which is an ordinary LDDMM call. Each later level starts from the previous level’s A and v on the same velocity grid xv, without a pre-fit, and deforms from its first iteration.

    • sigmaP (only with landmarks): multiplied by the level’s pooling factor, so that the balance between the landmark and the image terms does not change with the number of pixels.

    • draw: applies to the last level; the other levels are not drawn.

Returns:

out – LDDMM’s result for the last level (‘A’, ‘v’, ‘xv’, ‘WM’, ‘WB’, ‘WA’, …), with two more entries:

’rasters’dict

The density images and their pixel locations: ‘XI’, ‘YI’, ‘I’ for the source and ‘XJ’, ‘YJ’, ‘J’ for the target, as returned by rasterize. The fitted map takes [YI, XI] to [YJ, XJ] (row, column order, as LDDMM’s xI and xJ).

’levels’list of dict

One per level: ‘factor’ and ‘source_factor’ (the pooling factors of the target and the source images), ‘niter’, and the pooled pixel locations ‘xI’ and ‘xJ’. The weights ‘WM’, ‘WB’ and ‘WA’ are on the last level’s ‘xJ’.

If LDDMM reports its affine pre-fit (‘prefit’), the record is the first level’s.

Return type:

dict

Raises:

  • ValueError – If the cell coordinates do not match in length, levels is empty or invalid, or a level pools an image to fewer than 2 pixels along an axis.

  • TypeError – If niter is passed.

Notes

Re-rasterizing at a fine pixel size with a fixed physical kernel and fitting a two level pyramid was the most accurate configuration on our synthetic test problems (see CHANGES.md). Images from other sources (H&E, Visium, an atlas) go to LDDMM directly.

The matching weights are not carried over from one level to the next: each level starts from LDDMM’s initial weights.

Examples

>>> out = STalign.align_cells(xI, yI, xJ, yJ, L=L, T=T)
>>> tpointsI = STalign.transform_points_source_to_target(out['xv'], out['v'], out['A'],
...                                                      np.stack([yI, xI], 1))
STalign.STalign.analyze3Dalign(labelfile, xv, v, mat, xJ, dx, scale_x, scale_y, x, y, X_, Y_, namesdict, device='cpu')

Create dataframe with region annotations and 3D coordinates of target brain slice

Each cell is mapped into the atlas with the backward map fitted by LDDMM_3D_to_slice (build_transform3D, direction ‘b’), evaluated at the cell’s own position, and gets the label of the nearest atlas voxel at that position.

Parameters:

  • labelfile (str) – File name of 3D atlas with region annotations (nrrd).

  • xv (list of array) – Sample points for velocity

  • v (array) – time dependent velocity field

  • mat (array) – Affine transformation matrix

  • xJ (list of arrays) – Pixel locations of the rasterized target image. Not used (kept for compatibility).

  • dx (float) – Step of rasterized image. Not used (kept for compatibility).

  • scale_x (double) – Value that x positions of atlas were scaled by. Not used (kept for compatibility).

  • scale_y (double) – Value that y positions of atlas were scaled by. Not used (kept for compatibility).

  • x (array) – X positions of target brain slice

  • y (array) – Y positions of target brain slice

  • X (array) – Not used (kept for compatibility).

  • Y (array) – Not used (kept for compatibility).

  • namesdict (Dictionary) – Dictionary of brain regions corresponding to region numbers

  • device (str) – Not used (kept for compatibility); the computation runs on the device of v.

Returns:

df – One row per cell of the target brain slice: its 3D coordinates in the atlas (coord0, coord1, coord2), x, y, the region annotation (struct_id, 0 outside the atlas) and its acronym.

Return type:

DataFrame

STalign.STalign.backend_info()

Which LDDMM backends are available.

Returns:

info – ‘auto’: the backend that backend=’auto’ uses for CPU float64 problems (‘cpp’ or ‘torch’); ‘cpp’: whether the compiled extension can be imported; ‘cpp_error’: why not (or None); ‘num_threads’: the default maximum number of threads (see resolve_num_threads).

Return type:

dict

STalign.STalign.build_transform(xv, v, A, direction='b', XJ=None, exact_inverse=True)

Create sample points to transform source to target from affine and velocity.

Parameters:

  • xv (list of array) – Sample points for velocity

  • v (array) – time dependent velocity field, (nt,v0,v1,2). None for an affine only transform.

  • A (array) – Affine transformation matrix

  • direction (char) – ‘f’ for forward and ‘b’ for backward. ‘b’ is default and is used for transforming images. ‘f’ is used for transforming points.

  • XJ (array) – Sample points for target (meshgrid with ij index style). Defaults to None to keep sampling on the xv. Can also be a list of 1D arrays (one per axis), or an array of points with components along the last axis.

  • exact_inverse (bool) – Only used for direction ‘f’. If True (default), the forward map is the exact inverse of the backward map, to about 1e-6 of a velocity grid spacing. If False, it is the forward Euler integration used by earlier versions of STalign, which inverts the backward map only up to a time discretization error (up to tens of microns on the tutorial alignments).

Returns:

Xs – Sample points in mehsgrid format, (…,2).

Return type:

torch tensor

Notes

The backward map (‘b’, target to source) is the map LDDMM fits: the inverse affine transform, followed by explicit Euler steps of -v[t]/nt for t = nt-1 down to 0. The exact forward map inverts each step by fixed-point iteration (an implicit Euler step of the forward flow), starting from the forward Euler step. A RuntimeWarning is raised for points that do not converge, which can only happen if v is too large for nt time steps.

Inputs can be numpy arrays or torch tensors, with integer or floating point values. The computation uses the floating point dtype of A and v (float64 if both are integer) and the device of v.

STalign.STalign.build_transform3D(xv, v, A, direction='b', XJ=None, exact_inverse=True)

Create sample points to transform source to target from affine and velocity.

Parameters:

  • xv (list of array) – Sample points for velocity

  • v (array) – time dependent velocity field, (nt,v0,v1,v2,3). None for an affine only transform.

  • A (array) – Affine transformation matrix

  • direction (char) – ‘f’ for forward and ‘b’ for backward. ‘b’ is default and is used for transforming images. ‘f’ is used for transforming points.

  • XJ (array) – Sample points for target (meshgrid with ij index style). Defaults to None to keep sampling on the xv. Can also be a list of 1D arrays (one per axis), or an array of points with components along the last axis, for example (N,3).

  • exact_inverse (bool) – Only used for direction ‘f’. If True (default), the forward map is the exact inverse of the backward map, to about 1e-6 of a velocity grid spacing. If False, it is the forward Euler integration, which inverts the backward map only up to a time discretization error.

Returns:

Xs – Sample points in mehsgrid format, (…,3).

Return type:

torch tensor

Notes

After LDDMM_3D_to_slice the source is the atlas and the target is the slice: ‘b’ maps slice positions to the atlas and ‘f’ maps atlas positions to the slice. See build_transform for how the two maps are computed.

STalign.STalign.calculate_tre(pointsI, pointsJ)

Target registration error (TRE) between two sets of corresponding points

Parameters:

  • pointsI (array or torch tensor (N x 2)) – Points, e.g. the output of transform_points_source_to_target

  • pointsJ (array or torch tensor (N x 2)) – Corresponding points

Returns:

  • meanTRE (float) – Mean distance between corresponding points

  • stdTRE (float) – Standard deviation of the distances

STalign.STalign.clip(I)

clip an arrays values between 0 and 1. Useful for visualizing

Parameters:

I (torch tensor) – A torch tensor to clip.

Returns:

Ic – Clipped torch tensor.

Return type:

torch tensor

STalign.STalign.download_aba_image_labels(imageurl, labelurl, imagefile, labelfile)

Create 3D altas image and region annotations

Parameters:

  • imageurl (array) – Link containing atlas cell stained image.

  • labelurl (array) – Link containing atlas altas images for each voxel in imageurl.

  • imagefile (array) – File location to save imageurl information.

  • labelfile (array) – File location to save labelurl information.

Returns:

  • imagefile (array) – File location to save imageurl information.

  • labelfile (array) – File location to save labelurl information.

Raises:

requests.HTTPError – If a download fails (the file is not written then).

STalign.STalign.download_aba_ontology(url, file_name)

Create 3D altas ontology

Parameters:

  • url (Link to url contain atlas ontology)

  • file_name (File name to save atlas ontology)

Returns:

  • ontology_name (File name to store ontologies)

  • namesdict (Dictionary of all brain rgeion names)

Raises:

requests.HTTPError – If the download fails (nothing is written then).

STalign.STalign.extent_from_x(xJ)

Given a set of pixel locations, returns an extent 4-tuple for use with np.imshow.

Note inputs are locations of pixels along each axis, i.e. row column not xy.

Parameters:

xJ (list of torch tensors) – Location of pixels along each axis

Returns:

extent – (xmin, xmax, ymin, ymax) tuple

Return type:

tuple

Examples

>>> extent_from_x(xJ)
>>> fig,ax = plt.subplots()
>>> ax.imshow(J,extent=extentJ)
STalign.STalign.interp(x, I, phii, **kwargs)

Interpolate the 2D image I, with regular grid positions stored in x (1d arrays), at the positions stored in phii (3D arrays with first channel storing component)

Parameters:

  • x (list of arrays) – List of arrays storing the pixel locations along each image axis. convention is row column order not xy.

  • I (array) – Image array. First axis should contain different channels.

  • phii (array) – Sampling array. First axis should contain sample locations corresponding to each axis.

  • **kwargs (dict) – Other arguments fed into the torch interpolation function torch.nn.grid_sample

Returns:

out – The image I resampled on the points defined in phii.

Return type:

torch tensor

Notes

Convention is to use align_corners=True.

This uses the torch library.

STalign.STalign.interp3D(x, I, phii, **kwargs)

Interpolate the 3D image I, with regular grid positions stored in x (1d arrays), at the positions stored in phii (4D arrays with first channel storing component)

Parameters:

  • x (list of arrays) – List of arrays storing the pixel locations along each image axis. convention is row column order not xy.

  • I (array) – Image array. First axis should contain different channels.

  • phii (array) – Sampling array. First axis should contain sample locations corresponding to each axis.

  • **kwargs (dict) – Other arguments fed into the torch interpolation function torch.nn.grid_sample

Returns:

out – The image I resampled on the points defined in phii.

Return type:

torch tensor

Notes

Convention is to use align_corners=True.

This uses the torch library.

STalign.STalign.make_scree(W, name, p=6)

Create a scree plot for a given set of eigenvalues.

Parameters:

  • W (numpy array) – Eigenvalues in descending order

  • name (str) – Name of the dataset the eigenvalue originate from Will be included in plot title

  • p (int) – Number of eigenvalue from W to plot. Defaults to 6

Returns:

fig – Figure handle for scree plot

Return type:

matplotlib Figure

Raises:

ValueError – If p is larger than the number to eigenvalues in W.

STalign.STalign.normalize(arr, t_min=0, t_max=1)

Linearly normalizes an array between two specifed values.

Parameters:

  • arr (numpy array) – array to be normalized

  • t_min (int or float) – Lower bound of normalization range

  • t_max (int or float) – Upper bound of normalization range

Returns:

norm_arr – Array of the same shape as arr with normalized values. A constant arr (zero range) maps to t_min.

Return type:

numpy array

STalign.STalign.plot_brain_regions(df)

Plot all brain regions in target brain slice with different color.

Parameters:

df (DataFrame) – Dataframe containing each cell in original target brain slice, region annotations and 3D coordinates in terms of the altas.

Returns:

  • fig (matplotlib figure)

  • ax (matplotlib axes)

STalign.STalign.plot_subset_brain_regions(df, brain_regions)

Plot subset of brain regions in target brain slice with different color.

Parameters:

  • df (DataFrame) – Dataframe containing each cell in original target brain slice, region annotations and 3D coordinates in terms of the altas.

  • brain_regions (array) – Subset of brain regions of interest (i.e. [‘CA1’, ‘CP’, ‘DG-sg’])

Returns:

  • fig (matplotlib figure)

  • ax (matplotlib axes)

STalign.STalign.rasterize(x, y, g=array([1.]), dx=30.0, blur=None, expand=1.1, draw=10000, wavelet_magnitude=False, use_windowing=True, *, blur_um=None)

Rasterize a spatial transcriptomics dataset into a density image

Each cell adds a Gaussian kernel centered on the cell, normalized to unit sum over the pixels where it is evaluated and multiplied by the cell’s weight g.

Parameters:

  • x (numpy array of length N) – x location of cells

  • y (numpy array of length N) – y location of cells

  • g (numpy array of length N) – RNA count of cells If not given, density image is created. Unless g is all ones it is min-max normalized to [0, 1] (see normalize), so cells at the minimum of g contribute nothing, and a constant g gives an image of zeros. Use rasterize_with_signal for an image of unnormalized values.

  • dx (float) – Pixel size to rasterize data (default 30.0, in same units as x and y)

  • blur (float or list of floats) – Width of the Gaussian interpolation kernel, in pixels. The kernel standard deviation is 2*blur pixels (2*blur*dx in the units of x and y). Use a list to do multi scale (one channel per blur). Defaults to 1.0, unless blur_um is given.

  • expand (float) – Factor to expand sampled area beyond cells. Defaults to 1.1.

  • draw (int) – If nonzero, draw the image in a new figure and return its handle. Defaults to 10000. The figure is drawn once, after all cells are added (earlier versions redrew it every draw cells).

  • wavelet_magnitude (bool) – If True, take the absolute value of difference between scales for raster images. When using this option blur should be sorted from greatest to least.

  • use_windowing (bool) – If True (default), each cell’s kernel is evaluated on a (2r+1) x (2r+1) pixel window around the cell’s nearest pixel (clipped at the image border), r = ceil(4*max(blur)) pixels (about 2 standard deviations), and normalized to unit sum on that window. The window holds about 95% of the Gaussian, so the kernel acts like one about 4% narrower. If False, the untruncated kernel is evaluated on the whole image.

  • blur_um (float or list of floats, keyword only) – The kernel standard deviation in the units of x and y (for example micrometers), instead of blur: blur = blur_um/(2*dx). The kernel then has the same physical width at any pixel size dx. Give either blur or blur_um, not both.

Returns:

  • X (numpy array) – Locations of pixels along the x axis

  • Y (numpy array) – Locations of pixels along the y axis

  • M (numpy array) – A rasterized image with len(blur) channels along the first axis

  • fig (matplotlib figure handle) – Only returned if draw is nonzero: handle to the figure of the image.

Raises:

  • Exception – If wavelet_magnitude is set to true but blur is not sorted from greatest to least.

  • ValueError – If both blur and blur_um are given.

  • IndexError – If the image is empty (e.g. a single cell, or cells on a line).

Notes

Cells with NaN coordinates (possible with pandas inputs, whose min and max skip NaN) are ignored, with a warning.

Examples

Rasterize a dataset at 30 micron pixel size, with three kernels.

>>> X,Y,M,fig = STalign.rasterize(x,y,dx=30.0,blur=[2.0,1.0,0.5],draw=10000)

Rasterize a dataset at 30 micron pixel size, with three kernels, using difference between scales.

>>> X,Y,M,fig = STalign.rasterize(x,y,dx=30.0,blur=[2.0,1.0,0.5],draw=10000, wavelet_magnitude=True)

Rasterize a dataset at 15 micron pixel size with a kernel of standard deviation 30 microns (the same as blur=1.0).

>>> X,Y,M = STalign.rasterize(x,y,dx=15.0,blur_um=30.0,draw=0)
STalign.STalign.rasterizePCA(x, y, G)

Rasterize a spatial transcriptomics dataset into a density image and perform PCA on it.

Parameters:

  • x (numpy array of length N) – x location of cells

  • y (numpy array of length N) – y location of cells

  • G (pandas Dataframe of shape (N,M)) – gene expression level of cells for M genes

Returns:

  • X (numpy array) – A rasterized image for each gene with the channel along the first axis

  • Y (numpy array) – A rasterized image unraveled to 1-D for each gene; data is centered

  • W (numpy array) – The eigenvalues for the principal components in descending order

  • V (numpy array) – The normalized eigenvectors for the principal components V[:, i] corresponds to eigenvalue at W[i]

  • Z (numpy array) – Y rotated to align with the principal component axes (Z = V.T @ Y)

  • nrows (int) – Row dimension of rasterized image

  • ncols (int) – Column dimension of rasterized image

Notes

Each value/row at the same index in x, y, and G should all correspond to the same cell. x[i] <-> y[i] <-> G[i,:]

Each gene is rasterized with rasterize(x, y, g, dx=30.0, blur=1.0), which min-max normalizes its counts; a constant gene gives an image of zeros.

STalign.STalign.rasterize_with_signal(x, y, s=None, dx=30.0, blur=None, expand=1.1, draw=0, wavelet_magnitude=False, use_windowing=True, *, blur_um=None)

Rasterize a spatial transcriptomics dataset with a signal per cell

Each cell adds a Gaussian kernel centered on the cell, normalized to unit sum over the pixels where it is evaluated, divided by dx**2 and multiplied by the cell’s signal value(s), so the image is a signal density per unit area. See rasterize for the kernel and the window.

Parameters:

  • x (numpy array of length N) – x location of cells

  • y (numpy array of length N) – y location of cells

  • s (numpy array of length N, or N by M for M signals) – Signal value of each cell, used as is (not normalized). Signals may be signed, and the image keeps their sign. Each of the M columns gives one channel. If None, a density image is created.

  • dx (float) – Pixel size to rasterize data (default 30.0, in same units as x and y)

  • blur (float or list of floats) – Width of the Gaussian interpolation kernel, in pixels. The kernel standard deviation is 2*blur pixels (2*blur*dx in the units of x and y). A list does multi scale, which is only possible if s is None. Defaults to 1.0, unless blur_um is given.

  • expand (float) – Factor to expand sampled area beyond cells. Defaults to 1.1.

  • draw (int) – If nonzero, draw the image in a new figure and return its handle. Defaults to 0. The figure is drawn once, after all cells are added.

  • wavelet_magnitude (bool) – If True, take the absolute value of difference between scales for raster images. When using this option blur should be sorted from greatest to least. With a signal (a single blur) this only takes the absolute value.

  • use_windowing (bool) – If True (default), evaluate each kernel on a window of about 2 standard deviations; if False, on the whole image (see rasterize).

  • blur_um (float or list of floats, keyword only) – The kernel standard deviation in the units of x and y, instead of blur: blur = blur_um/(2*dx). Give either blur or blur_um, not both.

Returns:

  • X (numpy array) – Locations of pixels along the x axis

  • Y (numpy array) – Locations of pixels along the y axis

  • M (numpy array) – A rasterized image with channels along the first axis: one per column of s, or len(blur) if s is None

  • fig (matplotlib figure handle) – Only returned if draw is nonzero: handle to the figure of the image.

Raises:

  • Exception – If s is given with more than one blur, or if wavelet_magnitude is set to true but blur is not sorted from greatest to least.

  • ValueError – If both blur and blur_um are given.

  • IndexError – If s has fewer rows than there are cells, or if the image is empty.

Examples

Rasterize the total counts and the counts of two genes at 30 micron pixel size.

>>> X,Y,M = STalign.rasterize_with_signal(x,y,s=np.stack((total,gene1,gene2),-1),dx=30.0,blur=1.0)

Rasterize the density at three scales, using difference between scales.

>>> X,Y,M = STalign.rasterize_with_signal(x,y,dx=30.0,blur=[2.0,1.0,0.5],wavelet_magnitude=True)
STalign.STalign.saveRasters(X, Y, W, V, Z, scree_fig, nrows, ncols, name, path)

Create an RGB image principal components 1-3 and 4-6. Save numpy arrays necessary for rasterizing and performing PCA. Save scree plot and RGB image as jpeg, and their figure handles.

Parameters:

  • X (numpy array) – A rasterized image for each gene with the channel along the first axis

  • Y (numpy array) – A rasterized image unraveled to 1-D for each gene; data is centered

  • W (numpy array) – The eigenvalues for the principal components in descending order

  • V (numpy array) – The normalized eigenvectors for the principal components V[:, i] corresponds to eigenvalue at W[i]

  • Z (numpy array) – X rotated to align with the principal component axes

  • scree_fig (matplotlib Figure) – Figure handle for scree plot

  • nrows (int) – Row dimension of rasterized image

  • ncols (int) – Column dimension of rasterized image

  • name (str) – Name of the dataset the arrays originate from

  • path (str) – Path to folder where files should be saved (created if needed). Files are named <name>_screefig.jpg, <name>_rgbfig.jpg and <name>_arrays.npz.

Return type:

None

STalign.STalign.to_A(L, T)

Convert a linear transform matrix and a translation vector into an affine matrix.

Parameters:

  • L (torch tensor) – 2x2 linear transform matrix

  • T (torch tensor) – 2 element translation vector (note NOT 2x1)

Returns:

A – Affine transform matrix

Return type:

torch tensor

STalign.STalign.to_A_3D(L, T)

Convert a linear transform matrix and a translation vector into an affine matrix.

Parameters:

  • L (torch tensor) – 3x3 linear transform matrix

  • T (torch tensor) – 3 element translation vector (note NOT 2x1)

Returns:

A – Affine transform matrix

Return type:

torch tensor

STalign.STalign.transform_image_source_to_target(xv, v, A, xI, I, XJ=None)

Transform an image from source to target space, by sampling it at the backward map.

Parameters:

  • xv (list of array) – Sample points for velocity

  • v (array) – time dependent velocity field, (nt,v0,v1,2). None for an affine only transform.

  • A (array) – Affine transformation matrix

  • xI (list of arrays) – Locations of the pixels of I along each axis (row column order, not xy).

  • I (array) – Source image, channels along the first axis.

  • XJ (list of arrays) – Locations of the output pixels (see build_transform). Defaults to the velocity grid.

Returns:

phiI – I in target space, channels along the first axis.

Return type:

torch tensor

Notes

I is sampled with bilinear interpolation, clamped at its border. Inputs can be numpy arrays or torch tensors of any numeric dtype.

STalign.STalign.transform_image_source_with_A(A, XI, I, XJ)

Transform an image with an affine matrix

Parameters:

  • A (torch tensor) – Affine transform matrix

  • XI (list of numpy arrays) – List of arrays storing the pixel location in image I along each image axis. convention is row column order not xy. i.e, locations of pixels along the y axis (rows) followed by locations of pixels along the x axis (columns)

  • I (numpy array) – A rasterized image with len(blur) channels along the first axis

  • XJ (list of numpy arrays) – List of arrays storing the pixel location in image I along each image axis. convention is row column order not xy. i.e, locations of pixels along the y axis (rows) followed by locations of pixels along the x axis (columns)

Returns:

AI – image I after affine transformation A, with channels along first axis

Return type:

torch tensor

STalign.STalign.transform_image_target_to_source(xv, v, A, xJ, J, XI=None, exact_inverse=True)

Transform an image from target to source space, by sampling it at the forward map.

Parameters:

  • xv (list of array) – Sample points for velocity

  • v (array) – time dependent velocity field, (nt,v0,v1,2). None for an affine only transform.

  • A (array) – Affine transformation matrix

  • xJ (list of arrays) – Locations of the pixels of J along each axis (row column order, not xy).

  • J (array) – Target image, channels along the first axis.

  • XI (list of arrays) – Locations of the output pixels (see build_transform). Defaults to the velocity grid.

  • exact_inverse (bool) – If True (default), use the exact inverse of the backward map; if False, the forward Euler map of earlier versions (see build_transform).

Returns:

phiiJ – J in source space, channels along the first axis.

Return type:

torch tensor

Notes

J is sampled with bilinear interpolation, clamped at its border. Inputs can be numpy arrays or torch tensors of any numeric dtype.

STalign.STalign.transform_points_source_to_target(xv, v, A, pointsI, exact_inverse=True)

Transform points from source to target space. Note points are in row column order, not xy.

Parameters:

  • xv (list of array) – Sample points for velocity

  • v (array) – time dependent velocity field, (nt,*nv,d). None for an affine only transform.

  • A (array) – Affine transformation matrix, (d+1,d+1). With d=3 (after LDDMM_3D_to_slice) this maps atlas points to the slice.

  • pointsI (array) – (N,d) points in source space. Numpy array or torch tensor, integer or floating point; it is not modified.

  • exact_inverse (bool) – If True (default), use the exact inverse of the backward map LDDMM fits; if False, the forward Euler map of earlier versions (see build_transform).

Returns:

pointsIt – (N,d) points in target space.

Return type:

torch tensor

STalign.STalign.transform_points_target_to_source(xv, v, A, pointsI)

Transform points from target to source space. Note points are in row column order, not xy.

This is the backward map LDDMM fits, identical to build_transform(xv,v,A,direction=’b’) evaluated at the points.

Parameters:

  • xv (list of array) – Sample points for velocity

  • v (array) – time dependent velocity field, (nt,*nv,d). None for an affine only transform.

  • A (array) – Affine transformation matrix, (d+1,d+1).

  • pointsI (array) – (N,d) points in target space. Numpy array or torch tensor, integer or floating point; it is not modified.

Returns:

pointsIt – (N,d) points in source space.

Return type:

torch tensor

STalign.STalign.v_to_phii(xv, v)

Integrate a velocity field over time to return a position field (diffeomorphism).

Parameters:

  • xv (list of torch tensor) – List of 1D tensors describing locations of sample points in v

  • v (torch tensor) – 4D (nt,v0,v1,2) velocity field, the layout LDDMM returns. The component-first layout (nt,2,v0,v1) is also accepted (on a 2x2 grid, where the two are ambiguous, the LDDMM layout is assumed).

Returns:

phii – Inverse map (position field) computed by method of characteristics, (2,v0,v1)

Return type:

torch tensor

Notes

This semi-Lagrangian scheme approximates the map build_transform(xv,v,A,’b’) computes for an identity A, up to interpolation error. Use build_transform for exactly the map LDDMM fitted.

STalign.STalign.v_to_phii_3D(xv, v)

Integrate a 3D velocity field over time to return a 3D position field (diffeomorphism).

Parameters:

  • xv (list of torch tensor) – List of 1D tensors describing locations of sample points in v

  • v (torch tensor) – 5D (nt,v0,v1,v2,3) velocity field, the layout LDDMM_3D_to_slice returns. The component-first layout (nt,3,v0,v1,v2) is also accepted (on a 3x3x3 grid, where the two are ambiguous, the LDDMM layout is assumed).

Returns:

phii – Inverse map (position field) computed by method of characteristics, (3,v0,v1,v2)

Return type:

torch tensor

Notes

This semi-Lagrangian scheme approximates the map build_transform3D(xv,v,A,’b’) computes for an identity A, up to interpolation error. Use build_transform3D for exactly the map LDDMM_3D_to_slice fitted.