Fix subidx to enable computing the NoRMCorre template frame from subset of movie frames

[ludovicbellier]

Description

Hi team CaImAn! I'm Ludovic Bellier, a Senior Computational Research Scientist at Inscopix, a Bruker company. First, thank you so much for making available to the community such an amazing and powerful toolbox!

In several projects we process concatenated movies, and we've noticed that in some instances motion correction wasn't working well. Upon further investigation, it happened when there was a big shift in Field Of View (FOV) between two concatenated movies, and we tracked down the issue to the reference frame that is computed from a subset of frames uniformly sampled throughout the concatenated movies: with a big shift between movies, it ends up averaging frames at two FOVs relatively shifted from one another, and yields a "dual state", "echo" template frame. Here's an ASCII, 1D example, hoping it makes it clearer 😁

Movie 1 FOV:   ". . | . . ." (let's say | is a blood vessel, represented on pixel 3)
Movie 2 FOV:   ". . . | . ." (leftward shift, now | is on pixel 4)
Template frame: ". . | | . ." (averaging an equal number of frames from Movie 1 and 2 results in a blurry template frame, with the blood vessel both on pixel 3 and 4)

Then, during the motion correction iterations, frames from movie 1 will match the template frame well enough, as will frames from movie 2. In other words, there is no incentive for frames of each movie to be shifted to match the template, when the template has "echoed", duplicated landmarks.

After a few unsuccessful attempts at providing a single frame as template frame, we dug into the code and noticed there was an argument, subidx, meant to index a slice of movie frames. Only, it was not functioning properly.

In this present PR, we fixed this bug and restore the ability to compute the motion-correction template frame by uniformly sampling frames from a slice of the Tiff stack only. Using this code successfully fixed motion correction in all instances where we had substantial shifts between concatenated movies.

Type of change

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)

Has your PR been tested?

All motion correction tests in python caimanmanager.py test ran successfully.

Task list

• Decide on how to enable both temporal cropping (#1106) and template computation
• Implement for nonrigid MC
• Make the parameter(s) part of params.motion
• Write documentation in the args docs in MotionCorrect

[EricThomson]

Great thanks a lot for the contribution. I am about to go out of town for Scipy next week, so not sure when I will get a chance to review, but I look forward to checking this out!

[EricThomson]

Hi I'm back from Scipy, and will set aside some time this week to look this over. Thanks again for the contribution!

[ludovicbellier]

Hi Eric, sounds good, thank you for the update! Please don't hesitate if you have any questions.

[EricThomson]

This is a really helpful contribution, and relates to #1106, which is something I've been wanting to do since before I opened that issue.

I'm not sure how closely it implements the functionality in that issue, so I will want to test it out, and see if it works as intended in multiple use cases. It will probably take me a little while to run through this and stress test it, and I'm going on vacation next week (sorry! 😄), but I will start with it today, and probably tweak it in a few places, and may have some questions.

[EricThomson]

I have just had a little time to look it over, but a few things that hit me on first view.

On the usage of this param: I've always considered it to be the temporal analog of indices, which lets you do spatial cropping of your movie using a convenient slicing operation as input. I saw subidx as a similar temporal cropping internal that was never exposed on the front end (as discussed in #1106). As I play with your PR first iteration, it is not yet working that way (though it is good that it helps you solve the problem you have mentioned above: that seems like almost a secondary benefit compared to the more obvious temporal cropping/subsampling intention: e.g., see how it shapes the T variable (which is the time index)). Maybe I'm misunderstanding, or doing it wrong, though. If so, maybe you could discuss more or add some more use-case examples? My guess is it should just solve both simultaneously which would be great. Fingers crossed, but please let me know if I'm wrong.

Also, it looks like currently this PR has integrated only into rigid motion correction, but not nonrigid? Nonrigid also has the same param.

We should add documentation to the args docs in MotionCorrect (but this will depend on bigger picture question of usage I think).

This should become part of params.motion (just like indices is part of the parameter object). This is important for data lineage/provenance tracking, as it will ultimately likely change the data that is pulled after motion correction. So it has wider consequences beyond motion_correction.py (i.e., if I give a 3000 frame movie slice(0,1000) input for this, then after that my movie will only be 1000 frames long.

I will be out for a week+, but properly describing and exposing subidx on the front end is important and we need to get it right. I'll be back in a week and a half to help work on this. I appreciate you helping to take it on, as it is something I've sort of avoided as there are quite a few moving parts.

[ludovicbellier]

Hi Eric, you are absolutely right, subidx seemed to be originally intended for temporal cropping, which is a slightly different use than our current use case (i.e., computing the template frame from a subset of movie frames).

I can see two ways for enabling both use cases:

1. We create another parameter, e.g., subidx_template, for the template computation use case, and keep subidx for temporal cropping. A caveat is that then we need to decide whether subidx_template needs to be contained within subidx or not (e.g., on a 3k frame movie, we compute the template on frames 0 to 1k, and output a temporally cropped motion-corrected movie with frames 2k to 2999), but that's solvable.
2. We create another parameter, e.g., subidx_type, to specify how subidx should be interpreted ('crop' or 'template'). The issue with this is that we can't do both at the same time (as in the example in 1.).

What do you think? Happy to help either way!
Wishing you a wonderful vacation!

[EricThomson]

Thanks for these thoughts they help clarify the path forward: it seems Option 1 is the the right way to go.

I think subidx_template should just be a separate parameter from subidx: people will typically want their full movie motion corrected, and this is a way to have it done right. However, you are right to bring up the concern with how it might collide with subidx (handling when the intersection is null, or whatever).

I'm thinking maybe the right path forward is to implement subidx functionality as a baseline so we know that is working right and really test it, and then superimpose the subidx_template (or whatever) as an additional parameter once we better appreciate the full ramifications of exposing subidx on the front end (which as I said is something on the pile of things to do).

This would have the unfortunate consequence of slowing things down a bit, but I do think in the long run would be a very cool double feature. This sounds like two separate PRs to me: first just subidx and second, subidx_template, otherwise it will get too big and messy.

Does this sound like a reasonable path forward to you?

[EricThomson]

Note @ludovicbellier we've done quite a bit of work in the past week or so on getting subindices to work in a couple of out-of-the-way methods (load() and load_iter()) that are used for basic movie loading for motion correction. This was a lot more work than I thought it would be (it's not code that we play with a lot and had some issues to be worked out). It's in the dev_load_pims_fallback branch.

[ludovicbellier]

Hi @EricThomson, sorry for the late reply, I took a vacation in France for the first time in 5 years! =)
Your path forward sounds great, and thanks a lot for making so much progress on it!! Is there anything I can help with at this point? How far from merging do you estimate this subidx functionality is?

[EricThomson]

Hi @ludovicbellier, we are currently patching up a few things in dev and a couple of us are traveling. When we are all back in a couple of weeks we should hopefully merge. I'll report back once this is done.

[ludovicbellier]

Perfect, thanks @EricThomson for this update and looking forward to hearing back from you soon!

[EricThomson]

@ludovicbellier We've pushed the new version with subindexing fixed so I now have a better sense for how it is supposed to work 😄 Will now be able to loop around to this (probably after SFN unfortunately I have tons of things to prepare for the SFN workshop sorry this is taking so long!!).