Registration

Functions:

`apply_translations`(imagestack, corrections)	Apply (potentially non-integer) translations to each frame in an imagestack.
`distributed_translation_estimator`(D)	Compute registration among many frames, based on pairwise estimated registrations.
`dotproducts_at_translations`(X, Y[, cutin, ...])	Computes dot product between X and various translations of Y, after a little bit of preprocessing.
`find_translations_using_model`(imagestack, ...)	A method that uses the codebook and the model to find a translation of the imagestack which is more consistent with the observation model.
`pairwise_correlation_registration`(X, Y[, ...])	Estimate a translation that makes X and Y line up with each other.

bardensr.registration.apply_translations(imagestack, corrections, mode='valid', interpolation_method='linear')

Apply (potentially non-integer) translations to each frame in an imagestack.

Input:

imagestack (N x M0 x M1 x M2 numpy array)
corrections (N x 3 floating point numpy)
mode (‘valid’ or ‘full’; this indicates what to do with voxels for which not all frames have been measured. valid trims them out, full sets them to zero.)
interpolation_method (‘hermite’ or ‘linear’ or ‘nearest’; how to deal with cases where corrections are not integers)

Output:

imagestack2 (N x M0’ x M1’ x M2’)
trimmed_corrections (N x 3 array, indicating the coordinates in imagestack which are used to supply imagestack2[:,0,0,0]). This will be the same as corrections, up to an overall shift which may be applied to every row independently. In particular, trimmed_corrections[i,j]-trimmed_corrections[i+1,j] = corrections[i,j]-corrections[i+1,j].

imagestack2 is a translated version of imagestack which satisfies:

imagestack[f,trimmed_corrections[f,0],trimmed_corrections[f,1],trimmed_corrections[f,1]
   approx
imagestack2[f,0,0,0]

The shape of imagestack2 and the value of trimmed_corrections depend upon the value of ‘mode’ supplied. “valid” insists that every value of imagestack2 arises from a value in imagestack. When “full” is used, every value in imagestack is placed somewhere in imagestack2.

bardensr.registration.distributed_translation_estimator(D)

Compute registration among many frames, based on pairwise estimated registrations.

Input: D (… x F x F numpy array)
Output: T (… x F numpy array)

Solves the minimization problem:

argmin_t sum_ij (D_ij^2 - (t_i-t_j))^2

Ihe input should indicate the result of trying to find a translation that registers each pair of frames. For example, D_ij could represent the result of running pairwise_correlation_registration(data[i],data[j]). The output tries to find a translation for all of the frames that is as consistent as possible with all the pairwise estimates.

For more info, cf. E. Varol et al., “Decentralized Motion Inference and Registration of Neuropixel Data,” ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2021, pp. 1085-1089, doi: 10.1109/ICASSP39728.2021.9414145.

bardensr.registration.dotproducts_at_translations(X, Y, cutin=None, demean=True)

Computes dot product between X and various translations of Y, after a little bit of preprocessing.

Input:

X – M0 x M1 X … M(n-1) numpy array
Y – M0 x M1 x … M(n-1) numpy array
[optional] cutin – a scalar or n-vector indicating that we should only look at an inner portion of Y; this can help especially if X,Y aren’t demeaned
[optional] demean – whether or not to subtract mean from X and Y before estimating translations

Output:

V – an n-dimensional numpy array indicating dot product between X and various translations of Y.
offset – an n-vector

The offset indicates how the entries in V should be interpreted: V[i0,i1,i2…] indicates the dot product between X and Y after translating Y by (i0+offset0,i1+offset1,i2+offset2…). For example, If demean and cutin are False:

V[i0,i1,i2...] = sum_j X[j0,j1,j2...] * Y[j0-(i0+offset0),j1-(i1+offset1)...].

Here the sum is taken over all values such that the indices aren’t out-of-bounds.

NOTE: this is a thin wrapper around sp.signal.correlate, and it does not use GPUs.

bardensr.registration.find_translations_using_model(imagestack, codebook, maximum_wiggle=10, niter=50, use_tqdm_notebook=False, initial_guess=None)

A method that uses the codebook and the model to find a translation of the imagestack which is more consistent with the observation model. Before running this code, we generally advocate preprocessing by running bardensr.preprocess_minmax, running bardensr.preprocess_bgsubtraction and then running bardensr.preprocess_minmax again.

Input

imagestack (N x M0 x M1 x M2 numpy array)
codebook (N x J numpy array)
[optional] maximum_wiggle (tuple of 3 integers; default (10,10,10); maximum possible wiggle permitted along each spatial dimension)
[optional] niter (integer; default 50; number of rounds of gradient descent to run in estimating the registration)

Output: a FindTranslationsUsingModelResults, including

corrections (N x 3 numpy array, indicating how each imagestack should be shifted)
losses (loss computed at each iteration, indicating how well the gradient descent is proceeding)

bardensr.registration.pairwise_correlation_registration(X, Y, cutin=None, demean=True)

Estimate a translation that makes X and Y line up with each other.

Input:

X – M0 x M1 X … M(n-1) numpy array
Y – M0 x M1 x … M(n-1) numpy array
[optional] cutin – a scalar or n-vector indicating that we should only look at an inner portion of Y; this can help especially if X,Y aren’t demeaned
[optional] demean – whether or not to subtract mean from X and Y before estimating translations

Output: t, an n-vector indicating how Y should be shifted to look like X. Roughly speaking, t tries to make it so that:

X[i0,i1,i2...] approx Y[i0-t0,i1-t1,...]

by looking at dot products between X and various translations of Y. Put another way, if we run:

(Xp,Yp) = apply_translations([X,Y],[zeros(n),-t],'full')

then Xp and Yp should line up.

See dotproducts_at_translations for more information.

NOTE: this is a thin wrapper around sp.signal.correlate, and it does not use GPUs.