Image registration fails on low signal to noise data

[danielsf]

@morriscb and I have been experimenting with using Suite2P registration on low signal to noise 2-photon movies generated at the Allen Institute for Brain Science. We have discovered a peculiar failure mode. For some cases, register.compute_reference shifts the reference image so far that a significant fraction of the field of view gets wrapped around the edge of the image (see examples below). We believe this problem arises because, when the frames are iteratively shifted to create a new reference image in compute_reference, the shifts are added to the frames in place, allowing the shifts to accumulate over the 8 iterations

https://github.com/MouseLand/suite2p/blob/main/suite2p/registration/register.py#L148-L149

this renders the limitation implied by ops['maxregshift'] moot, since, for instance, a series of 0.1*Lx shifts could, in principle, add up to an 0.8*Lx shift over the course of the iterations, if they all pointed in the same direction. This is probably not a concern for high signal to noise data in which there are clear structures for the registration to latch on to but, in low signal to noise data, it is possible that the iteration will have a hard time converging and could easily "wander away" as described (again: see below).

The proposed fix would be something like this

[scott.daniel@synapse suite2p]$ git show 85b2ba
commit 85b2ba8dc6b95eaf65646e2a7f3a13e5a4dd1eec
Author: danielsf <[email protected]>
Date:   Thu Jan 6 22:46:46 2022 -0800

    option to turn off cumulatie shifts

diff --git a/suite2p/registration/register.py b/suite2p/registration/register.py
index 9b0ec50..97c2b27 100644
--- a/suite2p/registration/register.py
+++ b/suite2p/registration/register.py
@@ -97,7 +97,7 @@ def pick_initial_reference(frames: np.ndarray):
     return refImg
 
 
-def compute_reference(ops, frames):
+def compute_reference(ops, frames, cumulative_shifts=True):
     """ computes the reference image
 
     picks initial reference then iteratively aligns frames to create reference
@@ -126,8 +126,15 @@ def compute_reference(ops, frames):
         refImg = utils.spatial_high_pass(refImg, int(ops['spatial_hp_reg']))
         frames = utils.spatial_high_pass(frames, int(ops['spatial_hp_reg']))
 
+    if not cumulative_shifts:
+        input_frames = np.copy(frames)
+
     niter = 8
     for iter in range(0, niter):
+
+        if not cumulative_shifts:
+            frames = np.copy(input_frames)
+
         # rigid registration
         ymax, xmax, cmax = rigid.phasecorr(
             data=rigid.apply_masks(
@@ -529,4 +536,4 @@ def enhanced_mean_image(ops):
         ops['xrange'][0]:ops['xrange'][1]] = mimg0
     ops['meanImgE'] = mimg
     print('added enhanced mean image')
-    return ops
\ No newline at end of file
+    return ops

copy the frames passed to compute_reference and, at each iteration, return to using the copied frames rather than stacking the rigid frame transformations on top of each other. This may not be terribly memory efficient, since it involves carrying two copies of frames around at a time. This is not a huge concern in our use case, but maybe is for some of your other users, which is why we are not just issuing a direct pull request.

Examples In each of these cases, the leftmost image is the average projection of the raw (pre-registration) frames in the movie. The center image is the result of calling suite2p.register.compute_reference on the specified nimg_init number of frames. The rightmost image the result of calling suite2p.register.compute_reference with the above change applied (cumulative_shifts=False). Note that the reference images generated by Suite2P all wrap the image around one of its edges, whereas the reference images generated by the modified versions of Suite2P do not (or, at least, not as signficantly).

We are willing to talk with you more about this. We are probably pushing the edges of what Suite2P was intended for (we run registration with ops['maxregshift'] = 0.2, which probably doesn't help; we've found need for that based on past results).

PS While we were examining register.compute_reference, we discovered some choices that confused us. Why is isort indexed [1:nmax] here

https://github.com/MouseLand/suite2p/blob/main/suite2p/registration/register.py#L152

It effectively throws away the most correlated image when constructing the new reference image. We haven't experimented with this parameter, but were wondering what the rationale for this choice was. A similar choice appears to be made in register.pick_initial_reference here

https://github.com/MouseLand/suite2p/blob/main/suite2p/registration/register.py#L93

[carsen-stringer]

thanks for pointing that out, we need to fix that bug with [1:nmax], it's from when we converted the code from matlab :) . right now it breaks our regression tests so we haven't done that yet.

the keeping of the original frames seems like a reasonable fix. would you like to submit a pull request?