Apply offset values from Suite2P to raw movie

[wbwakeman]

Issue 268 determined that the smooth_sigma_time parameter has the potential to help a lot of experiments that may be very noisy, or have certain artifacts introduced by hardware problems. Values between 0 and 4 were shown to not introduce any additional problems that can be seen, but we would like to verify that there aren't any hidden problems introduced. The first step to do that will be to simply take the offset values that come out of suite2p motion correction and applying them manually to the raw movie. This will avoid the possibility of suite2p altering our data in any way other than exactly what we want.

Use Suite2P's shift_frame function: https://github.com/AllenInstitute/ophys_etl_pipelines/blob/44be06d5832fa09f4ea2d461b5257930b03608ce/src/ophys_etl/modules/suite2p_registration/__main__.py#L85-L93

Tasks:

• [x] Remove code that reads Suite2P output TIFFs
• [x] Implement shift_frame function on raw data read from the HDF5
• [x] Preserve or incorporate shift_frame implementation for clipped data

Acceptance criteria:

• [x] Generate outputs that clearly demonstrate that offset values generated by suite2p have been applied directly to the raw movie, and that Suite2P is not doing any other "processing" on the data

Update: The above task list has been completed, however the projection images come out blurry and still with the vertical artifact lines. The following task list is intended for debugging this problem.

• [x] Verify that frames are being shifted in the correct directions along each axis
• [x] View raw videos and verify if the artifacts appear there
• [x] Also make projection images of raw movie
• [x] Check that the csv outputs are identical to the non-mnual application csv outputs to verify that nothing else was changed.
• [x] Read the Suite2P documentation to see what else they did to these images other than just applying offsets

[danielsf]

As a proof of concept, I ran our wrapped Suite2P motion correction module (with smooth_sigma_time=4) on a cartoon movie with three bright ROIs of known fiducial position and motion. I then took the resulting suite2p_rigid_transformation.csv file and used np.roll to apply the prescribed translations in x and y directly to the raw movie. This video (which focuses on the period in the movie in which the cartoon ROIs actually move) shows the results.

The left panel is the raw movie; the central panel is the movie produced by Suite2P; the right panel is the result of applying the x,y translations by hand. The red dots mark the true, motion-corrected centers of the ROIs. If Suite2P were just applying a rigid motion correction, you would expect that, in all frames, the relationship between the gray squares and the red dots would be the same for all three ROIs in the central panel. If you scroll through the video by hand, you can find examples where this is not true, see, for example, the behavior between timestamp 0:02 and 0:03

https://images.zenhubusercontent.com/274240166/2939b272-3dab-4b6b-b2a6-f60fb050b1ff/video_side_by_side.mp4

As further indication that Suite2P is doing more than rigid motion correction, this plot zooms in on one ROI and three different timestamps. The left column shows the position of the ROI in the raw movie, the central column shows the position of the ROI in the Suite2P corrected movie, and the right column shows the position of the ROI in the by-hand corrected movie. The reported rigid motion displacements are reported in the title of the left hand column. Note that, naively interpreting the values of dx and dy in the second and third rows, you would expect the ROI to be displaced in only one dimension relative to the raw movie. The Suite2P ROIs are clearly displaced in both dimensions.

Finally, this figure shows the dx, dy values found by Suite2P relative to the true motion applied to the cartoon video

it does appear the Suite2P does a fairly good job of finding the correct displacements, so maybe the fact that Suite2P is doing something extra should not be too alarming.

The next step will be to try this by-hand correction on real data and see what the results look like.

[danielsf]

I have interpreted the scope of this ticket to be to explore the effects of Suite2P motion correction with smooth_sigma_time=4 to determine whether or not Suite2P does something beyond rigidly transform the frames of the raw movie (it does) and whether or not we can take the displacements reported by Suite2P and apply them to a movie by hand (we can). The visualizations below compare raw, un-motion-corrected movies to motion-corrected movies output by Suite2P with smooth_sigma_time=0 and smooth_sigma_time=4 to movies that are motion-corrected by applying Suite2P's reported rigid transformations by hand. The function to apply a rigid transformation to a movie by hand is:

def shift_by_hand(raw_data, dxdy_lookup):
    """
   dxdy_lookup is a dict mapping frame index to a dict like
    {'dx': dx, 'dy': dy}
    """
    new_data = np.zeros(raw_data.shape, dtype=raw_data.dtype)
    for iframe in range(new_data.shape[0]):
        if iframe not in dxdy_lookup:
            new_data[iframe, :, : ] = raw_data[iframe, :, :]
        else:
            dxdy = dxdy_lookup[iframe]
            dx = dxdy['dx']
            dy = dxdy['dy']
            new_frame = np.roll(raw_data[iframe, :, :],
                                (-dy, -dx),
                                axis=(0,1))
            new_data[iframe, :, : ] = new_frame

    return new_data

The experiments visualized below (ophys_experiment_id is reported in the title of upper left hand plot) all failed motion correction with the default smooth_sigma_time=0. The plots are as follows

• The first row is the raw, un-motion-corrected movie
• The second row is the motion-corrected movie produced by Suite2P with smooth_sigma_time=0
• The third row is the motion-corrected movie produced by Suite2P with smooth_sigma_time=4
• The fourth row is the motion-corrected movie produced by applying Suite2P's reported rigid transformations directly to the raw movie by hand

The columns of the plots are:

• The first column is the max projection image
• The second column is the average projection image
• The third column shows the dx, dy displacement curves as a function of frame number
• To create the fourth column, I randomly selected 100 frames. I took all of the pixels in those frames and sorted them. For each row, I plotted the histogram of the ratio of those sorted pixel values with the sorted pixel values taken from the un-motion-corrected movie (same set of 100 frames). If Suite2P only translated the pixels without altering their flux values, we would expect his histogram to be a spike at 1.0. If Suite2P renormalized the pixels by a constant factor, we would expect the histogram to be a spike at some other value. The fact that the histogram has non-zero extent indicates that Suite2P is redistributing flux between pixels in the course of doing registration. This is not an enormous effect (note the logarithmic scale on the histogram's vertical axis), but it is an effect. Note: before computing this histogram, I subtracted off the global minimum of each movie's flux value, just in case Suite2P was zeroing and then renormalizing the data. The histogram still has a non-zero extent.
• The fifth column reports some summary statistics. flux_ratio min/med/max reports the minimum, median, and maximum of the histograms in the fourth column. video min/max report the minimum and maximum values of flux in the full movie (I am surprised that negative values are allowed).

You will note that, by construction, the fourth row preserves all pixel flux values. It only applies rigid translations to the pixels' positions in the field of view.

[danielsf]

Note that np.roll operates by shifting pixels that fall off one edge onto the opposite edge. This is the same function that Suite2P uses in its rigid transform implementation here

https://github.com/MouseLand/suite2p/blob/main/suite2p/registration/rigid.py#L111-L129

In this experiment, the extremely large displacements mean that some bright regions are probably getting moved from the upper right corner of the raw max projection image to the lower central region of the motion-corrected max projection image.

[danielsf]

[danielsf]

[danielsf]

Note the abnormally large translation of features in the smooth_sigma_time=0 iteration of this experiment

[danielsf]

I believe the next step is to determine whether or not we want to accept whatever redistribution of flux Suite2P is doing on its own and, if not, alter our production module to apply the rigid transformation to the movie by hand using some variant of the method above.