z-stack registration

[jkim0731]

We are planning to take cortical z-stacks and local z-stacks using 'looping' protocol from now on.
Previously we used 'stepping' protocol, but recently we found out that 'looping' is more resistant to the disproportionate projection of neurons in neighboring planes due to spontaneous activities.

'step': [1,1,1,...,1], [2,2,2,...,2], [3,3,3,...,3], ..., [80,80,80,,...,80]
'loop': [1,2,3,...,80], [1,2,3,...,80], [1,2,3,...,80], ..., [1,2,3,...,80]

But due to larger time interval between frames in 'looping' protocol, we need to apply motion correction, both within planes and between neighboring planes.

Here's a code that worked for the test z-stack taken by looping.

def get_registered_z_stack_loop(loop_images, num_steps=160, num_loops=50):
"""Within-and-aross registration of a z-stack taken by looping protocol,
and return plane-averaged z-stack. Using phase cross correlation.
'step': [1,1,1,...,1], [2,2,2,...,2], [3,3,3,...,3], ..., [80,80,80,,...,80]
'loop': [1,2,3,...,80], [1,2,3,...,80], [1,2,3,...,80], ..., [1,2,3,...,80]

Parameters
----------
loop_images : np.ndarray (3d array)
    z-stack image taken by looping protocol. 1st dimension is the number of images
num_steps : int, optional
    Number of planes (or steps), by default 160
num_loops : int, optional
    Number of loops, by default 50

Returns
-------
np.ndimage (3d array)
    plane-averaged registered z-stack
"""
import skimage, scipy
assert loop_images.shape[0] == num_steps * num_loops
step_registered_mean = np.zeros((num_steps, *loop_images.shape[1:]))
for i in range(num_steps):
    step_images = loop_images[range(i,8000, num_steps),:,:]
    mean_step = np.mean(step_images, axis=0)
    step_registered = np.zeros_like(step_images)
    for j in range(num_loops):
        shift,_,_ = skimage.registration.phase_cross_correlation(mean_step, step_images[j,:,:], normalization=None)
        step_registered[j,:,:] = scipy.ndimage.shift(step_images[j,:,:], shift)
    step_mean = np.mean(step_registered, axis=0)
    if i > 0:
        shift,_,_ = skimage.registration.phase_cross_correlation(step_registered_mean[i-1,:,:], step_mean, normalization=None)
        step_mean = scipy.ndimage.shift(step_mean, shift)
    step_registered_mean[i,:,:] = step_mean
return step_registered_mean

Please add this code or something works like this to the pipeline for z-stack data.