emlddmm package

Submodules

emlddmm.curve_annotator module

emlddmm.emlddmm module

Run the EMLDDMM algorithm for deformable registration between two different imaging modalities with possible missing data in one of them

Details of this algorithm can be found in

  • [1] Tward, Daniel, et al. “Diffeomorphic registration with intensity transformation and missing data: Application to 3D digital pathology of Alzheimer’s disease.” Frontiers in neuroscience 14 (2020): 52.

  • [2] Tward, Daniel, et al. “3d mapping of serial histology sections with anomalies using a novel robust deformable registration algorithm.” Multimodal Brain Image Analysis and Mathematical Foundations of Computational Anatomy. Springer, Cham, 2019. 162-173.

  • [3] Tward, Daniel, et al. “Solving the where problem in neuroanatomy: a generative framework with learned mappings to register multimodal, incomplete data into a reference brain.” bioRxiv (2020).

  • [4] Tward DJ. An optical flow based left-invariant metric for natural gradient descent in affine image registration. Frontiers in Applied Mathematics and Statistics. 2021 Aug 24;7:718607.

Note all parameters are keyword arguments, but the first four are required.

param xI

xI[i] stores the location of voxels on the i-th axis of the atlas image I (REQUIRED)

type xI

list of arrays

param I

4D array storing atlas imaging data. Channels (e.g. RGB are stored on the first axis, and the last three are spatial dimensions. (REQUIRED)

type I

4D array (numpy or torch)

param xJ

xJ[i] stores the location of voxels on the i-th axis of the target image J (REQUIRED)

type xJ

list of arrays

param J

4D array storing target imaging data. Channels (e.g. RGB are stored on the first axis, and the last three are spatial dimensions. (REQUIRED)

type J

4D array (numpy or torch)

param nt

Number of timesteps for integrating a velocity field to yeild a position field (default 5).

type nt

int

param eA

Gradient descent step size for affine component (default 1e-5). It is strongly suggested that you test this value and not rely on defaults. Note linear and translation components are combined following [4] so only one stepsize is required.

type eA

float

param ev

Gradient descent step size for affine component (default 1e-5). It is strongly suggested that you test this value and not rely on defaults.

type ev

float

param order

Order of the polynomial used for contrast mapping. If using local contranst, only order 1 is supported.

type order

int

param n_draw

Draw a picture every n_draw iterations. 0 for do not draw.

type n_draw

int

param sigmaR

Amount of regularization of the velocity field used for diffeomorphic transformation, of the form 1/sigmaR^2 * (integral over time of norm velocity squared ).

type sigmaR

float

param n_iter

How many iterations of optimization to run.

type n_iter

int

param n_e_step

How many iterations of M step to run before another E step is ran in expectation maximization algorithm for detecting outliers.

type n_e_step

int

param v_start

What iteration to start optimizing velocity field. One may want to compute an affine transformation before beginning to compute a deformation (for example).

type v_start

int

param n_reduce_step

Simple stepsize reducer for gradient descent optimization. Every this number of steps, we check if objective function is oscillating. If so we reduce the step size.

type n_reduce_step

int

param v_expand_factor

How much bigger than the atlas image should the domain of the velocity field be? This is helpful to avoid wraparound effects caused by discrete Fourier domain calculations. 0.2 means 20% larger.

type v_expand_factor

float

param v_res_factor

How much lower resolution should the velocity field be sampled at than the atlas image. This is overrided if you specify dv.

type v_res_factor

float

param dv

Explicitly state the resolution of the sampling grid for the velocity field.

type dv

None or float or list of 3 floats

param a

Constant with units of length. In velocity regularization, its square is multiplied against the Laplacian. Regularization is of the form 1/2/sigmaR^2 int |(id - a^2 Delta)^p v_t|^2_{L2} dt.

type a

float

param p

Power of the Laplacian operator in regularization of the velocity field. Regularization is of the form 1/2/sigmaR^2 int |(id - a^2 Delta)^p v_t|^2_{L2} dt.

type p

float

returns

out (dict) – Returns a dictionary of outputs storing computing transforms. if full_outputs==True, then more data is output including figures.

raises Exception

If the initial velocity does not have three components.

raises Exception

Local contrast transform requires either order = 1, or order > 1 and 1D atlas.

raises Exception

If order > 1. Local contrast transform not implemented yet except for linear.

raises Exception

Amode must be 0 (normal), 1 (rigid), or 2 (rigid+scale), or 3 (rigid using XJ for projection).

emlddmm.emlddmm.__dir__()

Default dir() implementation.

emlddmm.emlddmm.__format__(format_spec, /)

Default object formatter.

emlddmm.emlddmm.__init_subclass__()

This method is called when a class is subclassed.

The default implementation does nothing. It may be overridden to extend subclasses.

emlddmm.emlddmm.__new__(*args, **kwargs)

Create and return a new object. See help(type) for accurate signature.

emlddmm.emlddmm.__reduce__()

Helper for pickle.

emlddmm.emlddmm.__reduce_ex__(protocol, /)

Helper for pickle.

emlddmm.emlddmm.__sizeof__()

Size of object in memory, in bytes.

emlddmm.emlddmm.__subclasshook__()

Abstract classes can override this to customize issubclass().

This is invoked early on by abc.ABCMeta.__subclasscheck__(). It should return True, False or NotImplemented. If it returns NotImplemented, the normal algorithm is used. Otherwise, it overrides the normal algorithm (and the outcome is cached).

emlddmm.histsetup module

emlddmm.load_image_folder module

emlddmm.manual_point_align module

emlddmm.point_mapper module

emlddmm.setup module

emlddmm.transformation_graph_v01 module

emlddmm.twoDreg module

Module contents

class emlddmm.Image(space, name, fpath, mask=None, x=None)[source]

Bases: object

space

name of the image space

Type

string

name

image name. This is provided when instantiating an Image object.

Type

string

x

image voxel coordinates

Type

list of numpy arrays

data

image data

Type

numpy array

title

image title passed by the read_data function

Type

string

names

information about image data dimensions

Type

list of strings

mask

image mask

Type

array

path

path to image file or directory containing the 2D image series

Type

string

normalize(norm='mean', q=0.99)

normalize image

downsample(down)[source]

downsample image, image coordinated, and mask

fnames()[source]

get filenames of 2D images in a series, or return the single filename of an image volume

downsample(down)[source]

Downsample image

Parameters

down (list of ints) – Factor by which to downsample along each dimension

Returns

  • x (list of numpy arrays) – Pixel locations where each element of the list identifies pixel locations in corresponding axis.

  • data (numpy array) – image data

  • mask (numpy array) – binary mask array

fnames()[source]

Get a list of image file names for 2D series, or a single file name for volume image.

Returns

fnames (list of strings) – List of image file names

class emlddmm.Transform(data, direction='f', domain=None, dtype=torch.float32, device='cpu', verbose=False, **kwargs)[source]

Bases: object

A simple class for storing and applying transforms

Note that the types of transforms we can support are

  1. Deformations stored as a displacement field loaded from a vtk file. These should be a 1x3xrowxcolxslice array.

  2. Deformations stored as a position field in a python variable. These are a 3xrowxcolxslice array.

  3. Velocity fields stored as a python variable. These are ntx3xrowxcolxslice arrays.

  4. a 4x4 affine transform loaded from a text file.

  5. a 4x4 affine transform stored in a python variable

  6. a nslices x 3 x 3 sequence of affine transforms stored in an array. * new and special case *

Note the data stores position fields or matrices. If it is a vtk file it will store a position field.

Raises

  • Exception – If transform is not a txt or vtk file or valid python variable.

  • Exception – if direction is not ‘f’ or ‘b’.

  • Exception – When inputting a velocity field, if the domain is not included.

  • Exception – When inputting a matrix, if its shape is not 3x3 or 4x4.

  • Exception – When specifying a mapping, if the direction is ‘b’ (backward).

apply(X)[source]
emlddmm.affine_from_figure(x, shift=1000.0, angle=5.0)[source]

Build small affine transforms by looking at figures generated in draw.

Parameters

  • x (str) –

    Two letters. “t”,”m”,”b” for top row middle row bottom row.

    Then ‘w’,’e’,’n’,’s’ for left right up down

    Or ‘r’,’l’ for turn right turn left

    or ‘id’ for identity.

  • shift (float) – How far to shift.

  • angle (float) – How far to rotate (in degrees)

Returns

A0 (numpy array) – 4x4 numpy array affine transform matrix

emlddmm.apply_transform_float(x, I, Xout, **kwargs)[source]

Apply transform to image Image points stored in x, data stored in I transform stored in Xout

There is an issue with numpy integer arrays, I’ll have two functions

emlddmm.apply_transform_from_file_to_points(q, tform_file)[source]

To transform points from spacei to spacej (example from fluoro to nissl_registerd) We look for the output folder called outputs/spacei/spacej_to_spacei/transforms (example outputs/fluoro/nissl_registered_to_fluoro/transforms) Note “spacej to spacei” is not a typo (even though it looks backwards, point data uses inverse of transform as compared to images). In the transforms folder, there are transforms of the form “spacei to spacej”. If applying transforms to slice data, you will have to find the appropriate slice.

Parameters

  • q (numpy array) – A Nx3 numpy array of coordinates in slice,row,col order

  • tform_file (str) – A string pointing to a transform file

Returns

Tq (numpy array) – The transformed set of points.

emlddmm.apply_transform_int(x, I, Xout, double=True, **kwargs)[source]

Apply transform to image Image points stored in x, data stored in I transform stored in Xout

There is an issue with numpy integer arrays, I’ll have two functions

Note that we often require double precision when converting to floats and back

Raises

Exception – If mode is not ‘nearest’ for ints.

emlddmm.atlas_free_reconstruction(**kwargs)[source]

Atlas free slice alignment

Uses an MM algorithm to align slices. Minimizes a Sobolev norm over rigid transformations of each slice.

All arguments are keword arguments

Keword Arguments

xJlist of arrays

Lists of pixel locations in slice row col axes

Jarray

A C x slice x row x col set of 2D image.

Warray

A slice x row x col set of weights of 2D images. 1 for good quality, to 0 for bad quality or out of bounds.

afloat

length scale (multiplied by laplacian) for smootheness averaging between slices. Defaults to twice the slice separation.

pfloat

Power of laplacian in somothing. Defaults to 2.

drawint

Draw figure (default false)

n_stepsint

Number of iterations of MM algorithm.

**kwargsdict

All other keword arguments are passed along to slice matching in the emlddmm algorithm.

returns

  • out (dictionary) – All outputs from emlddmm algorithm, plus the new image I, and reconstructed image Jr

  • out[‘I’] (numpy array) – A C x slice x row x col set of 2D images (same size as J), averaged between slices.

  • out[‘Jr’] (numpy array) – A C x slice x row x col set of 2D images (same size as J)

Todo

Think about how to do this with some slices fixed. Think about normalizing slice to slice contrast.

emlddmm.compose_sequence(transforms, Xin, direction='f', verbose=False)[source]

Compose a set of transformations, to produce a single position field, suitable for use with emlddmm.apply_transform_float() for example.

Parameters

  • transforms –

    Several types of inputs are supported.

    1. A list of transforms class.

    2. A list of filenames (single direction in argument)

    3. A list of a list of 2 tuples that specify direction (f,b)

    4. An output directory

  • Xin (3 x slice x row x col array) – The points we want to transform (e.g. sample points in atlas). Also supports input as a list of voxel locations, along each axis which will be reshaped as above using meshgrid.

  • direction (char) – Can be ‘f’ for foward or ‘b’ for bakward. f is default which maps points from atlas to target, or images from target to atlas.

Returns

Xout (3 x slicex rowxcol array) – Points from Xin that have had a sequence of transformations applied to them.

Note

Note, if the input is a string, we assume it is an output directory and get A and V. In this case we use the direction argument. If the input is a tuple of length 2, we assume it is an output directory and a direction

Otherwise, the input must be a list. It can be a list of strings, or transforms, or string-direction tuples.

We check that it is an instace of a list, so it should not be a tuple.

Raises

Exception – Transforms must be either output directory, or list of objects, or list of filenames, or list of tuples storing filename/direction.

Todo

  1. use os path join

  2. support direction as a list, right now direction only is used for a single direction

emlddmm.compute_atlas_from_slices(J, W, ooop, niter=10, I=None, draw=False)[source]

Construct an atlas image by averaging between slices.

This uses an MM (also could be interpreted as EM) algorithm for converting a weighted least squares problem to an ordinary least squares problem. The latter can be solved as a stationary problem, and updated iteratively.

Parameters

  • J (array) – an C x slice x row x col image array.

  • W (array) – An slice x row x col array of weights

  • ooop (array) – a 1D array with “slice” elements. This is a frequency domain one over operator for doing smoothing.

  • niter (int) – Number of iterations to update in MM algorithm. (default to 10)

  • I (array) – An C x slice x row x col image array representing an initial guess.

  • draw (int) – Draw the image every draw iterations. Do not draw if 0 (default).

Note

This code uses numpy, not torch, since it is not gradient based.

Todo

Better boundary conditions

emlddmm.convert_points_from_json(points, d_high, n_high=None, sidecar=None, z=None, verbose=False)[source]

We load points from a json produced by Samik at cold spring harbor.

These are indexed to pixels in a high res image, rather than any physical units.

To convert to proper points in 3D for transforming, we need information about pixel size and origin.

Pixel size of the high res image is a required input.

If we have a json sidecar file that was prepared for the registration dataset, we can get all the info from this.

If not, we get it elsewhere.

Origin information can be determined from knowing the number of pixels in the high res image.

Z coordinate information is not required if we are only applying 2D transforms, but for 3D it will have to be input manually if we do not have a sidecar file.

Parameters

  • points (str or numpy array) – either a geojson filename, or a Nx2 numpy array with coordinates loaded from such a file.

  • d_high (float) – pixel size of high resolution image where cells were detected

  • n_high (str or numpy array) – WIDTH x HEIGHT of high res image. Or the filename of the high res image.

  • sidecar (str) – Filename of sidecar file to get z and origin info

  • z (float) – z coordinate

Returns

q (numpy array) – A Nx3 array of points in physical units using our coordinate system convention (origin in center).

Notes

If we are applying 3D transforms, we need a z coordinate. This can be determined either by specifying it or by using a sidecar file. If we have neither, it ill be set to 0.

Todo

Consider setting z to nan instead of zero.

emlddmm.downmedian(xI, S_, down)[source]

Downsamples a 3D image by taking the median among rectangular neighborhoods. This is often appropriate when image pixels has a small number of outliers, or when pixel values are assumed to be ordered, but don’t otherwise belong to a vector space.

Note

2D images can be hanled by adding a singleton dimension no leading batch dimensions

Parameters

  • xI (list of 3 numpy arrays) – Locations of image pixels along each axis

  • S (numpy array) – Numpy array storing imaging data. Note there should not be a leading dimension for channels.

  • down (list of 3 ints) – downsample by this factor along each axis.

Returns

  • xd (list of 3 numpy arrays) – Locations of image pixels along each axis after downsampling.

  • Sd (numpy array) – The downsampled image.

emlddmm.downmode(xI, S_, down)[source]

Downsamples a 3D image by taking the mode among rectangular neighborhoods. This is appropriate for label images, where averaging pixel values is not meaningful.

Note

2D images can be hanled by adding a singleton dimension no leading batch dimensions

Parameters

  • xI (list of 3 numpy arrays) – Locations of image pixels along each axis

  • S (numpy array) – Numpy array storing imaging data. Note there should not be a leading dimension for channels.

  • down (list of 3 ints) – downsample by this factor along each axis.

Returns

  • xd (list of 3 numpy arrays) – Locations of image pixels along each axis after downsampling.

  • Sd (numpy array) – The downsampled image.

emlddmm.downsample(I, down, W=None)[source]

Downsample an image by an integer factor along each axis. Note extra data at the end will be truncated if necessary.

If the first axis is for image channels, downsampling factor should be 1 on this.

Parameters

  • I (array (numpy or torch)) – Imaging data to downsample

  • down (list of int) – List of downsampling factors for each axis.

  • W (array (numpy or torch)) – A weight of the same size as I but without the “channel” dimension

Returns

Id (array (numpy or torch as input)) – Downsampled imaging data.

emlddmm.downsample_ax(I, down, ax, W=None)[source]

Downsample imaging data along one of the first 5 axes.

Imaging data is downsampled by averaging nearest pixels. Note that data will be lost from the end of images instead of padding. This function is generally called repeatedly on each axis.

Parameters

  • I (array like (numpy or torch)) – Image to be downsampled on one axis.

  • down (int) – Downsampling factor. 2 means average pairs of nearest pixels into one new downsampled pixel

  • ax (int) – Which axis to downsample along.

  • W (np array) – A mask the same size as I, but without a “channel” dimension

Returns

Id (array like) – The downsampled image.

Raises

Exception – If a mask (W) is included and ax == 0.

emlddmm.downsample_image_domain(xI, I, down, W=None)[source]

Downsample an image as well as pixel locations

Parameters

  • xI (list of numpy arrays) – xI[i] is a numpy array storing the locations of each voxel along the i-th axis.

  • I (array like) – Image to be downsampled

  • down (list of ints) – Factor by which to downsample along each dimension

  • W (array like) – Weights the same size as I, but without a “channel” dimension

Returns

  • xId (list of numpy arrays) – New voxel locations in the same format as xI

  • Id (numpy array) – Downsampled image.

Raises

Exception – If the length of down and xI are not equal.

emlddmm.draw(J, xJ=None, fig=None, n_slices=5, vmin=None, vmax=None, disp=True, cbar=False, slices_start_end=[None, None, None], **kwargs)[source]

Draw 3D imaging data.

Images are shown by sampling slices along 3 orthogonal axes. Color or grayscale data can be shown.

Parameters

  • J (array like (torch tensor or numpy array)) – A 3D image with C channels should be size (C x nslice x nrow x ncol) Note grayscale images should have C=1, but still be a 4D array.

  • xJ (list) – A list of 3 numpy arrays. xJ[i] contains the positions of voxels along axis i. Note these are assumed to be uniformly spaced. The default is voxels of size 1.0.

  • fig (matplotlib figure) – A figure in which to draw pictures. Contents of the figure will be cleared. Default is None, which creates a new figure.

  • n_slices (int) – An integer denoting how many slices to draw along each axis. Default 5.

  • vmin – A minimum value for windowing imaging data. Can also be a list of size C for windowing each channel separately. Defaults to None, which corresponds to tha 0.001 quantile on each channel.

  • vmax – A maximum value for windowing imaging data. Can also be a list of size C for windowing each channel separately. Defaults to None, which corresponds to tha 0.999 quantile on each channel.

  • disp (bool) – Figure display toggle

  • kwargs (dict) – Other keywords will be passed on to the matplotlib imshow function. For example include cmap=’gray’ for a gray colormap

Returns

  • fig (matplotlib figure) – The matplotlib figure variable with data.

  • axs (array of matplotlib axes) – An array of matplotlib subplot axes containing each image.

Example

Here is an example:

>>> example test

Todo

Put interpolation=’none’ in keywords

emlddmm.emlddmm(**kwargs)[source]

Run the EMLDDMM algorithm for deformable registration between two different imaging modalities with possible missing data in one of them

Details of this algorithm can be found in

  • [1] Tward, Daniel, et al. “Diffeomorphic registration with intensity transformation and missing data: Application to 3D digital pathology of Alzheimer’s disease.” Frontiers in neuroscience 14 (2020): 52.

  • [2] Tward, Daniel, et al. “3d mapping of serial histology sections with anomalies using a novel robust deformable registration algorithm.” Multimodal Brain Image Analysis and Mathematical Foundations of Computational Anatomy. Springer, Cham, 2019. 162-173.

  • [3] Tward, Daniel, et al. “Solving the where problem in neuroanatomy: a generative framework with learned mappings to register multimodal, incomplete data into a reference brain.” bioRxiv (2020).

  • [4] Tward DJ. An optical flow based left-invariant metric for natural gradient descent in affine image registration. Frontiers in Applied Mathematics and Statistics. 2021 Aug 24;7:718607.

Note all parameters are keyword arguments, but the first four are required.

Parameters

  • xI (list of arrays) – xI[i] stores the location of voxels on the i-th axis of the atlas image I (REQUIRED)

  • I (4D array (numpy or torch)) – 4D array storing atlas imaging data. Channels (e.g. RGB are stored on the first axis, and the last three are spatial dimensions. (REQUIRED)

  • xJ (list of arrays) – xJ[i] stores the location of voxels on the i-th axis of the target image J (REQUIRED)

  • J (4D array (numpy or torch)) – 4D array storing target imaging data. Channels (e.g. RGB are stored on the first axis, and the last three are spatial dimensions. (REQUIRED)

  • nt (int) – Number of timesteps for integrating a velocity field to yeild a position field (default 5).

  • eA (float) – Gradient descent step size for affine component (default 1e-5). It is strongly suggested that you test this value and not rely on defaults. Note linear and translation components are combined following [4] so only one stepsize is required.

  • ev (float) – Gradient descent step size for affine component (default 1e-5). It is strongly suggested that you test this value and not rely on defaults.

  • order (int) – Order of the polynomial used for contrast mapping. If using local contranst, only order 1 is supported.

  • n_draw (int) – Draw a picture every n_draw iterations. 0 for do not draw.

  • sigmaR (float) – Amount of regularization of the velocity field used for diffeomorphic transformation, of the form 1/sigmaR^2 * (integral over time of norm velocity squared ).

  • n_iter (int) – How many iterations of optimization to run.

  • n_e_step (int) – How many iterations of M step to run before another E step is ran in expectation maximization algorithm for detecting outliers.

  • v_start (int) – What iteration to start optimizing velocity field. One may want to compute an affine transformation before beginning to compute a deformation (for example).

  • n_reduce_step (int) – Simple stepsize reducer for gradient descent optimization. Every this number of steps, we check if objective function is oscillating. If so we reduce the step size.

  • v_expand_factor (float) – How much bigger than the atlas image should the domain of the velocity field be? This is helpful to avoid wraparound effects caused by discrete Fourier domain calculations. 0.2 means 20% larger.

  • v_res_factor (float) – How much lower resolution should the velocity field be sampled at than the atlas image. This is overrided if you specify dv.

  • dv (None or float or list of 3 floats) – Explicitly state the resolution of the sampling grid for the velocity field.

  • a (float) – Constant with units of length. In velocity regularization, its square is multiplied against the Laplacian. Regularization is of the form 1/2/sigmaR^2 int |(id - a^2 Delta)^p v_t|^2_{L2} dt.

  • p (float) – Power of the Laplacian operator in regularization of the velocity field. Regularization is of the form 1/2/sigmaR^2 int |(id - a^2 Delta)^p v_t|^2_{L2} dt.

Returns

out (dict) – Returns a dictionary of outputs storing computing transforms. if full_outputs==True, then more data is output including figures.

Raises

  • Exception – If the initial velocity does not have three components.

  • Exception – Local contrast transform requires either order = 1, or order > 1 and 1D atlas.

  • Exception – If order > 1. Local contrast transform not implemented yet except for linear.

  • Exception – Amode must be 0 (normal), 1 (rigid), or 2 (rigid+scale), or 3 (rigid using XJ for projection).

emlddmm.emlddmm_multiscale(**kwargs)[source]

Run the emlddmm algorithm multiple times, restarting with the results of the previous iteration. This is intended to be used to register data from coarse to fine resolution.

Parameters

  • value (emlddmm parameters either as a list of length 1 (to use the same) –

  • at (at each iteration) or a list of length N (to use different values) –

  • iterations). (each of the N) –

Returns

A list of emlddmm outputs (see documentation for emlddmm)

emlddmm.extent_from_x(xJ)[source]

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

Note

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 (tuple) – (xmin, xmax, ymin, ymax) tuple

Example

Draw a 2D image stored in J, with pixel locations of rows stored in xJ[0] and pixel locations of columns stored in xJ[1].

>>> import matplotlib.pyplot as plt
>>> extent_from_x(xJ)
>>> fig,ax = plt.subplots()
>>> ax.imshow(J,extent=extentJ)
emlddmm.find_bounding_box(Ji, startoff=30, search=20, n_iter=10, bbox=None, fixsize=False)[source]

Computes a bounding box using an Otsu intraclass variance objective function.

Parameters

  • Ji (numpy array) – Row x col x n_channels numpy array to find bounding box

  • startoff (int) – How far away from boundary to put initial guess. Default 30.

  • search (int) – How far to move bounding box edges at each iteration. Default +/- 10.

  • n_iter (int) – How many loops through top, left, bottom, right.

Returns

bbox (list) – A tuple of 4 ints. [row0, col0, row1, col1].

emlddmm.fnames(path)[source]

Get a list of image file names for 2D series, or a single file name for volume image.

Returns

fnames (list of strings) – List of image file names

emlddmm.initialize_A2d_with_bbox(J, xJ)[source]

Use bounding boxes to find an initial guess of A2d.

On each slice we will compute a bounding box.

Then we will compute a translation vector which will move the slice to the center of the field of view.

Then we will return the inverse.

TODO

emlddmm.interp(x, I, phii, interp2d=False, **kwargs)[source]

Interpolate an image with specified regular voxel locations at specified sample points.

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

Parameters

  • x (list of numpy arrays) – x[i] is a numpy array storing the pixel locations of imaging data along the i-th axis. Note that this MUST be regularly spaced, only the first and last values are queried.

  • I (array) – Numpy array or torch tensor storing 2D or 3D imaging data. In the 3D case, I is a 4D array with channels along the first axis and spatial dimensions along the last 3. For 2D, I is a 3D array with spatial dimensions along the last 2.

  • phii (array) – Numpy array or torch tensor storing positions of the sample points. phii is a 3D or 4D array with components along the first axis (e.g. x0,x1,x1) and spatial dimensions along the last axes.

  • interp2d (bool, optional) – If True, interpolates a 2D image, otherwise 3D. Default is False (expects a 3D image).

  • kwargs (dict) – keword arguments to be passed to the grid sample function. For example to specify interpolation type like nearest. See pytorch grid_sample documentation.

Returns

out (torch tensor) – Array storing an image with channels stored along the first axis. This is the input image resampled at the points stored in phii.

emlddmm.labels_to_rgb(S, seed=0, black_label=0, white_label=255)[source]

Convert an integer valued label image into a randomly colored image for visualization with the draw function

Parameters

  • S (numpy array) – An array storing integer labels. Expected to be 4D (1 x slices x rows x columns), but can be 3D (slices x rows x columns).

  • seed (int) – Random seed for reproducibility

  • black_label (int) – Color to assign black. Usually for background.

emlddmm.load_slices(target_name, xJ=None)[source]

Load a slice dataset.

Load a slice dataset for histology registration. Slice datasets include pairs of images and json sidecar files, as well as one tsv file explaining the dataset. Note this code creates a 3D array by padding.

Parameters

  • target_name (string) – Name of a directory containing slice dataset.

  • xJ (list, optional) – list of numpy arrays containing voxel positions along each axis. Images will be resampled by interpolation on this 3D grid.

Returns

  • xJ (list of numpy arrays) – Location of v

  • J (numpy array) – Numpy array of size C x nslices x nrows x ncols where C is the number of channels e.g. C=3 for RGB.

  • W0 (numpy array) – A nslices x nrows x ncols numpy array containing weights. Weights are 0 where there was padding

Raises

Exception – If the first image is not present in the image series.

emlddmm.map_image(emlddmm_path, root_dir, from_space_name, to_space_name, input_image_fname, output_image_directory=None, from_slice_name=None, to_slice_name=None, use_detjac=False, verbose=False, **kwargs)[source]

This function will map imaging data from one space to another.

There are four cases:

  1. 3D to 3D mapping: A single displacement field is used to map data

  2. 3D to 2D mapping: A single displacement field is used to map data, a slice filename is needed in addition to a space

  3. 2D to 2D mapping: A single matrix is used to map data.

  4. 2D to 3D mapping: Currently not supported. Ideally this will output data, and weights for a single slice, so it can be averaged with other slices.

Warning

This function was built for a particular use case at Cold Spring Harbor, and is generally not used. It may be removed in the future.

Parameters

  • emlddmm_path (str) – Path to the emlddmm python library, used for io

  • root_dir (str) – The root directory of the output structure

  • from_space_name (str) – The name of the space we are mapping data from

  • to_space_name (str) – The name of the space we are mapping data to

  • input_image_fname (str) – Filename of the input image to be transformed

  • output_image_fname (str) – Filename of the output image after transformation. If None (default), it will be returned as a python variable but not written to disk.

  • from_slice_name (str) – When transforming slice based image data only, we also need to know the filename of the slice the data came from.

  • to_slice_name (str) – When transforming slice based image data only, we also need to know the filename of the slice the data came from.

  • use_detjac (bool) – If the image represents a density, it should be transformed and multiplied by the Jacobian of the transformation

Keyword Arguments

**kwargs (dict) – Arguments passed to torch interpolation (grid_resample), e.g. padding_mode,

Returns

phiI (array) – Transformed image

Raises

  • Exception – If use_detjac set to True for 3D to 2D mapping. Detjac not currently supported for 3D to 2D.

  • Exception – 2D to 3D not implemented yet, may not get implemented.

  • Exception – Jacobian not implemented yet for 2D slices.

Warns

DetJac is ignored if mapping is 2D to 2D.

emlddmm.map_points(emlddmm_path, root_dir, from_space_name, to_space_name, input_points_fname, output_points_directory=None, from_slice_name=None, to_slice_name=None, use_detjac=False, verbose=False, **kwargs)[source]

For points we need to get the transforms in the opposite folder to images.

This function will map imaging data from one space to another. There are four cases:

  1. 3D to 3D mapping: A single displacement field is used to map data

  2. 3D to 2D mapping: Currently not supported.

  3. 2D to 2D mapping: A single matrix is used to map data.

  4. 2D to 3D mapping: A single displacement field is used to map data, a slice filename is needed in addition to a space

Warning

This function was built for a particular use case at Cold Spring Harbor, and is generally not used. It may be removed in the future.

Parameters

  • emlddmm_path (str) – Path to the emlddmm python library, used for io

  • root_dir (str) – The root directory of the output structure

  • from_space_name (str) – The name of the space we are mapping data from

  • to_space_name (str) – The name of the space we are mapping data to

  • input_points_fname (str) – Filename of the input image to be transformed

  • output_directory_fname (str) – Filename of the output image after transformation. If None (default), it will be returned as a python variable but not written to disk.

  • from_slice_name (str) – When transforming slice based image data only, we also need to know the filename of the slice the data came from.

  • to_slice_name (str) – When transforming slice based image data only, we also need to know the filename of the slice the data came from.

  • use_detjac (bool) – If the image represents a density, it should be transformed and multiplied by the Jacobian of the transformation

  • **kwargs (dict) – Arguments passed to torch interpolation (grid_resample), e.g. padding_mode,

Returns

  • phiP (array) – Transformed points

  • connectivity (list of lists) – Same connectivity entries as loaded data

  • connectivity_type (str) – Same connectivity type as loaded data

Raises

  • Exception – 3D to 2D mapping is not implemented for points.

  • Exception – If use_detjac set to True for 2D to 3D mapping. Detjac not currently supported for 2D to 3D

  • Exception – If use_detjac set to True. Jacobian is not implemented yet.

  • Exception – If attempting to map points from 3D to 2D.

Warns

DetJac is ignored if mapping is 2D to 2D.

emlddmm.orientation_to_RAS(orientation, verbose=False)[source]

Compute a linear transform from a given orientation to RAS.

Orientations are specified using 3 letters, by selecting one of each pair for each image axis: R/L, A/P, S/I

Parameters

orientation (3-tuple) – orientation can be any iterable with 3 components. Each component should be one of R/L, A/P, S/I. There should be no duplicates

Returns

Ao (3x3 numpy array) – A linear transformation to transform your image to RAS

emlddmm.orientation_to_orientation(orientation0, orientation1, verbose=False)[source]

Compute a linear transform from one given orientation to another.

Orientations are specified using 3 letters, by selecting one of each pair for each image axis: R/L, A/P, S/I

This is done by computing transforms to and from RAS.

Parameters

orientation (3-tuple) – orientation can be any iterable with 3 components. Each component should be one of R/L, A/P, S/I. There should be no duplicates

Returns

Ao (3x3 numpy array) – A linear transformation to transform your image from orientation0 to orientation1

emlddmm.pad(xI, I, n, **kwargs)[source]

Pad an image and its domain.

Perhaps include here

Parameters

  • xI (list of arrays) – Location of pixels in I

  • I (array) – Image

  • n (list of ints or list of pairs of ints) – Pad on front and back

emlddmm.project_affine_to_rigid(A, XJ)[source]

This function finds the closest rigid transform to the given affine transform.

Close is defined in terms of the action of A^{-1} on XJ

That is, we find R to minimize || A^{-1}XJ - R^{-1}XJ||^2_F = || R (A^{-1}XJ) - XJ ||^2_F.

We use a standard procurstes method.

Parameters

  • A (torch tensor) – A 4x4 affine transformation matrix

  • XJ (torch tensor) – A 3 x slice x row x col array of coordinates in the space of the target image

Returns

R (torch tensor) – A 4x43 rigid affine transformation matrix.

emlddmm.read_data(fname, x=None, **kwargs)[source]

Read array data from several file types.

This function will read array based data of several types and output x,images,title,names. Note we prefer vtk legacy format, but accept some other formats as read by nibabel.

Parameters

  • fname (str) – Filename (full path or relative) of array data to load. Can be .vtk or nibabel supported formats (e.g. .nii)

  • x (list of arrays, optional) – Coordinates for 2D series space

  • **kwargs (dict) – Keyword parameters that are passed on to the loader function

Returns

  • x (list of numpy arrays) – Pixel locations where each element of the list identifies pixel locations in corresponding axis.

  • images (numpy array) – Imaging data of size channels x slices x rows x cols, or of size time x 3 x slices x rows x cols for velocity fields

  • title (str) – Title of the dataset (read from vtk files)

  • names (list of str) – Names of each dataset (channel or time point)

Raises

  • Exception – If file type is nrrd.

  • Exception – If data is a single slice, json reader does not support it.

  • Exception – If opening with Nibabel and the affine matrix is not diagonal.

emlddmm.read_matrix_data(fname)[source]

Read linear transforms as matrix text file. Note in python we work in zyx order, but text files are in xyz order

Parameters

fname (str) –

Returns

A (array) – matrix in zyx order

emlddmm.read_vtk_data(fname, normalize=False, endian='b')[source]

Read vtk structured points legacy format data.

Note endian should always be big, but we support little as well.

Parameters

  • fname (str) – Name of .vtk file to read.

  • normalize (bool) – Whether or not to divide an image by its mean absolute value. Defaults to True.

  • endian (str) – Endian of data, with ‘b’ for big (default and only officially supported format) or ‘l’ for little (for compatibility if necessary).

Returns

  • x (list of numpy arrays) – Location of voxels along each spatial axis (last 3 axes)

  • images (numpy array) – Image with last three axes corresponding to spatial dimensions. If 4D, first axis is channel. If 5D, first axis is time, and second is xyz component of vector field.

Raises

  • Exception – The first line should include vtk DataFile Version X.X

  • Exception – If the file contains data type other than BINARY.

  • Exception – If the dataset type is not STRUCTURED_POINTS.

  • Exception – If the dataset does not have either 3 or 4 axes.

  • Exception – If dataset does not contain POINT_DATA

  • Exception – If the file does not contain scalars or vectors.

Warns

If data not written in big endian – Note (b) symbol not in data name {name}, you should check that it was written big endian. Specify endian=”l” if you want little

Todo

Torch does not support negative strides. This has lead to an error where x has a negative stride. I should flip instead of negative stride, or copy afterward.

emlddmm.read_vtk_polydata(fname)[source]

Read ascii vtk polydata from simple legacy files. Assume file contains xyz order, they are converted to zyx for python

Parameters

fname (str) – The name of the file to read

Returns

  • points (numpy float array) – nx3 array storing locations of points

  • connectivity (list of lists) – list of indices containing connectivity elements

  • connectivity_type (str) – VERTICES or LINES or POLYGONS

  • name (str) – name of the dataset

emlddmm.registered_domain(x, A2d)[source]

Construct a new domain that fits all rigidly aligned slices.

Parameters

  • x (list of arrays) – list of numpy arrays containing voxel positions along each axis.

  • A2d (numpy array) – Nx3x3 array of affine transformations

Returns

xr (list of arrays) – new list of numpy arrays containing voxel positions along each axis

emlddmm.reshape_for_local(J, local_contrast)[source]

Reshapes an image into blocks for simple local contrast estimation.

Parameters

  • J (tensor) – 3D image data where first index stores the channel information (i.e. 4D array)

  • local_contrast (tensor) – 1D tensor storing the block size on each dimension

Returns

Jv (tensor) – Reshaped imaging data to be used for contrast estimation

emlddmm.reshape_from_local(Jv, local_contrast=None)[source]

After changing contrast, transform back TODO: this did not get used

emlddmm.rigid2D(xI, I, xJ, J, **kwargs)[source]

Rigid transformation between 2D slices.

emlddmm.sinc_resample_numpy(I, n)[source]

Perform sinc resampling of an image in numpy.

This function does sinc resampling using numpy rfft torch does not let us control behavior of fft well enough This is intended to be used to resample velocity fields if necessary Only intending it to be used for upsampling.

Parameters

  • I (numpy array) – An image to be resampled. Can be an arbitrary number of dimensions.

  • n (list of ints) – Desired dimension of output data.

Returns

Id (numpy array) – A resampled image of size n.

emlddmm.v_to_phii(xv, v, **kwargs)[source]

Use Euler’s method to construct a position field from a velocity field by integrating over time.

This method uses interpolation and subtracts and adds identity for better behavior outside boundaries. This method is sometimes refered to as the method of characteristics.

Parameters

  • xv (list of 1D tensors) – xv[i] is a tensor storing the location of the sample points along the i-th dimension of v

  • v (5D tensor) – 5D tensor where first axis corresponds to time, second corresponds to component, and 3rd to 5th correspond to spatial dimensions.

Returns

phii (4D tensor) – Inverse transformation is output with component on the first dimension and space on the last 3. Note that the whole timeseries is not output.

emlddmm.weighted_intraclass_variance(bbox, Ji, ellipse=True)[source]

Returns weighted intraclass variance (as in Otsu’s method) for an image with a given bounding box.

Parameters

  • bbox (list) – A tuple of 4 ints. [row0, col0, row1, col1].

  • Ji (numpy array) – Row x col x n_channels numpy array to find bounding box

Returns

E (float) – The weighted intraclass variance between inside and outside the bounding box.

emlddmm.write_data(fname, x, out, title, names=None)[source]

Write data

Raises

  • Exception – If data is in .nii format, it must be grayscale.

  • Exception – If output is not .vtk or .nii/.nii.gz format.

Warns

If data to be written uses extension .nii or .nii.gz – Writing image in nii fomat, no title or names saved

emlddmm.write_matrix_data(fname, A)[source]

Write linear transforms as matrix text file. Note that in python we use zyx order, but we write outputs in xyz order

Parameter

fnamestr

Filename to write

A2D array

Matrix data to write. Assumed to be in zyx order.

returns

None

emlddmm.write_outputs_for_pair(output_dir, outputs, xI, I, xJ, J, WJ=None, atlas_space_name=None, target_space_name=None, atlas_image_name=None, target_image_name=None)[source]

Write outputs in standard format for a pair of images

Parameters

  • output_dir (str) – Location to store output data.

  • outputs (dict) – Dictionary of outputs from the emlddmm python code

  • xI (list of numpy array) – Location of voxels in atlas

  • I (numpy array) – Atlas image

  • xJ (list of numpy array) – Location of voxels in target

  • J (numpy array) – Target image

  • atlas_space_name (str) – Name of atlas space (default ‘atlas’)

  • target_space_name (str) – Name of target space (default ‘target’)

  • atlas_image_name (str) – Name of atlas image (default ‘atlasimage’)

  • target_image_name (str) – Name of target image (default ‘targetimage’)

Todo

Implement QC figures.

Check device more carefully, probably better to put everything on cpu.

Check dtype more carefully.

Get determinant of jacobian. (done)

Write velocity out and affine

Get a list of files to name outputs.

emlddmm.write_qc_outputs(output_dir, output, I, J, xS=None, S=None)[source]

write out registration qc images

Parameters

  • output_dir (string) – Path to output parent directory

  • output (dict) – Output dictionary from emlddmm algorithm

  • I (emlddmm image) – source image

  • J (emlddmm image) – target image

  • xS (list of arrays, optional) – Label coordinates

  • S (array, optional) – Labels

Returns

None

emlddmm.write_transform_outputs(output_dir, output, I, J)[source]

Note, daniel is redoing the above slightly. on March 2023

Write transforms output from emlddmm. Velocity field, 3D affine transform, and 2D affine transforms for each slice if applicable.

Parameters

  • output_dir (str) – Directory to place output data (will be created of it does not exist)

  • output (dict) – Output dictionary from emlddmm algorithm

  • A2d_names (list) – List of file names for A2d transforms

Returns

None

Todo

Update parameter list.

emlddmm.write_transform_outputs_(output_dir, output, I, J)[source]

Note this version (whose function name ends in “_” ) is obsolete. Write transforms output from emlddmm. Velocity field, 3D affine transform, and 2D affine transforms for each slice if applicable.

Parameters

  • output_dir (str) – Directory to place output data (will be created of it does not exist)

  • output (dict) – Output dictionary from emlddmm algorithm

  • A2d_names (list) – List of file names for A2d transforms

Returns

None

Todo

Update parameter list.

emlddmm.write_vtk_data(fname, x, out, title, names=None)[source]

Write data as vtk file legacy format file. Note data is written in big endian.

inputs should be numpy, but will check for tensor only structured points supported, scalars or vectors data type each channel is saved as a dataset (time for velocity field, or image channel for images) each channel is saved as a structured points with a vector or a scalar at each point

Parameters

  • fname (str) – filename to write to

  • x (list of arrays) – Voxel locations along last three axes

  • out (numpy array) – Imaging data to write out. If out is size nt x 3 x slices x height x width we assume vector if out is size n x slices x height x width we assume scalar

  • title (str) – Name of the dataset

  • names (list of str or None) – List of names for each dataset or None to use a default.

Raises

Exception – If out is not the right size.

emlddmm.write_vtk_polydata(fname, name, points, connectivity=None, connectivity_type=None)[source]

points should by Nx3 in zyx order It will be written out in xyz order connectivity should be lists of indices or nothing to write only cell data

Parameters

  • fname (str) – Filename to write

  • name (str) – Dataset name

  • points (array) –

  • connectivity (str) – Array of arrays storing each connectivity element as integers that refer to the points, size number of points by number of dimensions (expected to be 3)

  • connectivity_type (str) – Can by VERTICES, or POLYGONS, or LINES

Returns

nothing