Add Motion Correction to 3D demo

[marcusflisher]

For better support, please use the template below to submit your issue. When your issue gets resolved please remember to close it.

Sometimes errors while running CNMF occur during parallel processing which prevents the log to provide a meaningful error message. Please reproduce your error with setting dview=None.

If you need to upgrade CaImAn follow the instructions given in the documentation.

• Tell us a bit about your setup:

  1. Operating system (Linux/macOS/Windows): macOS Catalina (10.15.4)

  2. Python version (3.x): 3.7.6

  3. Working environment (Python IDE/Jupyter Notebook/other): Jupyter Notebook

  4. Which of the demo scripts you're using for your analysis (if applicable): demo_pipeline.ipynb

  5. CaImAn version*: Most recent version

  6. CaImAn installation process (pip install ./pip install -e ./conda): pip install -e .

*You can get the CaImAn version by creating a params object and then typing params.data['caiman_version']. If the field doesn't exist, type N/A and consider upgrading)

• Describe the issue that you are experiencing

Hello,

Thank you for providing us with the CaImAn toolbox. I am new to python and coding so please excuse any novice mistakes.

I would like to use rigid motion correction and the CMNF algorithm to analyze my 4D dataset (TXYZ).

The dataset (TXYZ) is 181 x 350 x 350 x 8 and stored as a tif file.

The issue that I am having occurs while attempting to perform motion correction. Here is the notebook I am using: https://github.com/marcusflisher/CaImAn/blob/master/demo_pipeline_3D_MF_v1.ipynb

How can I successfully perform motion correction for this 3D dataset?

Let me know if it is useful to share the file that I am working with.

• Copy error log below

IndexError                                Traceback (most recent call last)
<ipython-input-10-4d1e407438e4> in <module>
      1 #%% Run piecewise-rigid motion correction using NoRMCorre
----> 2 mc.motion_correct(save_movie=True)
      3 m_els = cm.load(mc.fname_tot_els)
      4 border_to_0 = 0 if mc.border_nan is 'copy' else mc.border_to_0
      5     # maximum shift to be used for trimming against NaNs

~/CaImAn/caiman/motion_correction.py in motion_correct(self, template, save_movie)
    252                                     np.max(np.abs(self.y_shifts_els))))
    253         else:
--> 254             self.motion_correct_rigid(template=template, save_movie=save_movie)
    255             b0 = np.ceil(np.max(np.abs(self.shifts_rig)))
    256         self.border_to_0 = b0.astype(np.int)

~/CaImAn/caiman/motion_correction.py in motion_correct_rigid(self, template, save_movie)
    303                 var_name_hdf5=self.var_name_hdf5,
    304                 is3D=self.is3D,
--> 305                 indices=self.indices)
    306             if template is None:
    307                 self.total_template_rig = _total_template_rig

~/CaImAn/caiman/motion_correction.py in motion_correct_batch_rigid(fname, max_shifts, dview, splits, num_splits_to_process, num_iter, template, shifts_opencv, save_movie_rigid, add_to_movie, nonneg_movie, gSig_filt, subidx, use_cuda, border_nan, var_name_hdf5, is3D, indices)
   2670 
   2671     if is3D:
-> 2672         m = m[:, indices[0], indices[1], indices[2]]
   2673     else:
   2674         m = m[:, indices[0], indices[1]]

IndexError: too many indices for array

Any help would be greatly appreciated.

• If you're not reporting an error, type your message below

[j-friedrich]

Hello @marcusflisher Please make sure you have indeed the most recent caiman version. I used your notebook with some simulated data of the same size as yours and the IndexError does not occur with a fresh pull from the master branch. When creating the CNMFParams objects opts = params.CNMFParams(params_dict=opts_dict) you should get a message setting key indices automatically to (slice(None, None, None), slice(None, None, None), slice(None, None, None)) Please verify that opts.motion['indices'] is indeed a 3-tuple.

Note that given that you have only 8 z-layers, these might be spaced so far apart that it could potentially make more sense to process each layer individually.

[marcusflisher]

Hello @j-friedrich

Thank you for your assistance.

I uninstalled and pulled the most recent instance of CaImAn and ran my data back through my modified notebook that I previously posted but received the same error log as in my original post. I also tried to feed into my modified notebook the ‘demoMovie3D.tif’ file that is generated in the ‘demo_caiman_cnmf_3D’ notebook, but received the same error log as in my original post.

Here is a google drive link to the ‘demoMovie3D.tif’ file that I am using: https://drive.google.com/file/d/1oQSIOzvrqvgE-bBavheWSRMW77KCGn5_/view?usp=sharing

When setting the CNMFParams object to opts = params.CNMFParams(params_dict=opts_dict) I do indeed get the message:

4359372 [params.py: check_consistency():936] [15587] is3D=True, hence setting key indices automatically to (slice(None, None, None), slice(None, None, None), slice(None, None, None))

"Please verify that opts.motion['indices'] is indeed a 3-tuple." Calling opts.motion[‘indices’] returns (slice(None, None, None), slice(None, None, None), slice(None, None, None))

Although my original dataset is only 8 Z-slices thick, the sections are in close proximity to one another and a single cell will typically occupy multiple Z slices, which is why I want to use 3D motion correction and 3D CNMF analysis.

Do you have any other suggestions?

[j-friedrich]

Oh, indeed, loading subindices of a file goes awry for 3D tif files. No need to wait for a fix; please simply save your file in memory-mapped format first prior to motion correction by replacing the line mc = MotionCorrect(fnames, dview=dview, **opts.get_group('motion')) (where fnames is your tif file) with

fname_new = cm.save_memmap([fnames], is_3D=True)
mc = MotionCorrect(fname_new, dview=dview, **opts.get_group('motion'))            

[marcusflisher]

Thank you for your timely response.

After following your suggestion, I now get the following error when I attempt to run motion correction (for both my own 3D dataset and the ‘demoMovie3D.tif’ file that I posted above):

AttributeError                            Traceback (most recent call last)
<ipython-input-11-4d1e407438e4> in <module>
      1 #%% Run piecewise-rigid motion correction using NoRMCorre
      2 mc.motion_correct(save_movie=True)
----> 3 m_els = cm.load(mc.fname_tot_els)
      4 border_to_0 = 0 if mc.border_nan is 'copy' else mc.border_to_0
      5     # maximum shift to be used for trimming against NaNs

AttributeError: 'MotionCorrect' object has no attribute 'fname_tot_els'

How should I best go about assigning a value to 'fname_tot_els', including specific location of assignment and format?

Or do you have any other suggestions?

[j-friedrich]

For rigid motion correction the filename of the motion corrected file is
not mc.fname_tot_els but mc.fname_tot_rig, thus use the latter.

[marcusflisher]

Thank you.

I now receive the following error log after attempting to motion correct both the ‘demoMovie3D.tif’ file and my 3D dataset:

ValueError                                Traceback (most recent call last)
<ipython-input-14-bf1b52bfe074> in <module>
      1 #%% Run piecewise-rigid motion correction using NoRMCorre
      2 mc.motion_correct(save_movie=True)
----> 3 m_els = cm.load(mc.fname_tot_rig)
      4 border_to_0 = 0 if mc.border_nan is 'copy' else mc.border_to_0
      5     # maximum shift to be used for trimming against NaNs

~/CaImAn/caiman/base/movies.py in load(file_name, fr, start_time, meta_data, subindices, shape, var_name_hdf5, in_memory, is_behavior, bottom, top, left, right, channel, outtype, is3D)
   1490                                 outtype=outtype,
   1491                                 var_name_hdf5=var_name_hdf5,
-> 1492                                 is3D=is3D)
   1493 
   1494     elif isinstance(file_name,tuple):

~/CaImAn/caiman/base/movies.py in load_movie_chain(file_list, fr, start_time, meta_data, subindices, var_name_hdf5, bottom, top, left, right, z_top, z_bottom, is3D, channel, outtype)
   1801                 m = m[np.newaxis, :, :]
   1802 
-> 1803             _, h, w = np.shape(m)
   1804             m = m[:, top:h - bottom, left:w - right]
   1805         else:

ValueError: too many values to unpack (expected 3)

Do you have any suggestions for troubleshooting this?

[j-friedrich]

You have to use m_els = cm.load(mc.fname_tot_rig, is3D=True). Note that gSig has to be a 3-tuple, but in your shared notebook it isn't

[marcusflisher]

Thank you so much @j-friedrich for all of your help.

For the most part, I have successfully ran my 3D data through the 3D motion correction and 3D CNMF.

Here is a GitHub link to my updated notebook: https://github.com/marcusflisher/CaImAn/blob/master/demo_pipeline_3D_MF_v2.ipynb

The only problem now is that the X dimension and Z dimension were swapped (see screenshot of 3D component view below).

I am not sure where exactly this occurred. Do you know how I can fix this?

Once this is solved, I will clean up and update this version of my notebook, in case anyone else finds it useful for working with their 3D data.

[j-friedrich]

We should probably change the default axis along with to slice the volume into layers, but you can specify it using e.g. cnm.estimates.nb_view_components_3d(image_type='mean', dims=dims, axis=2)

[marcusflisher]

Thanks so much @j-friedrich! I am just learning python, and realized that I can use ? to access the various arguments for each function.

I was able to get the CNMF to work really well with my 3D dataset! However, I was only able to get good results when I directly loaded, memory mapped, and ran CNMF on my dataset (without motion correcting).

I noticed that the correlation image of non-motion correcting vs. motion correcting suggests I either have inappropriate parameter settings or am missing a key argument setting. I have attached the two respective correlation images below so you can see what I am talking about. I will read over the NoRMCorre paper and relevant documentation for motion correction to see if I can self troubleshoot. However, if you know of any obvious cause for this result, I would love to hear your suggestions.

Regardless, I will post back if I figure out what is causing the poor resulting correlation image after performing motion correction.

Correlation image without motion correction:

Correlation image with motion correction:

[j-friedrich]

Such figures typically occur when the order of dimensions is mixed up when converting between the original movie (a 4D tensor) and the (2D) matrices goes awry. Please check whether the motion corrected movie looks fine, e.g. by displaying the first z-layer cm.load(fname, is3D=True)[...,0].play() where fname is the filename of the original or motion corrected movie respectively. There’s a flag swap_dim when calculating the correlation image, try changing that.

[marcusflisher]

When I try to specify a z-layer by calling cm.load(fname, is3D=True)[...,0].play() the resulting playback only seems to show a single Y-row of pixels as seen in the screen shot below. The time series of this single Y-row vs the entire X dimension seemingly plays through properly, but it is hard to tell with only one Y-row:

Changing the the flag swap_dimfrom True to False didn't fix the problem. I think it is now plotting the Y-axis against the T-dimension.

Also, do you know how I can copy/add a false z-layer below and above the bottom and top z-layer or include the bottom and top z-layer in the CNMF initialization/extraction of components? When I run CNMF using the following commands, I don't think the bottom and top z-layer are being incorporated.

# INIT
cnm = cnmf.CNMF(n_processes, k=K, gSig=gSig, merge_thresh=merge_thresh, p=p, dview=dview)

# FIT
images = np.reshape(Yr.T, [T] + list(dims), order='F')    # reshape data in Python format (T x X x Y x Z)
cnm = cnm.fit(images)

# View the results
cnm.estimates.nb_view_components_3d(image_type='mean', dims=dims)

[j-friedrich]

Your (aptly named) tif file is in TZXY order, hence you have to do cm.load(frames, is3D=True)[:,0].play() to select the first z-layer. You selected from the last dimension of the array, i.e. a y-layer.

[marcusflisher]

Thank you. Running cm.load(fnames, is3D=True)[:,0].play() worked for proper playback of the original movie.

However, I am still struggling with the motion corrected video playback and I still get the funky correlation image that I originally posted.

I don't get any errors after running:

%%capture
#%% Run rigid motion correction using NoRMCorre
mc.motion_correct(save_movie=True)

But when I try to playback the motion corrected video using the following function I get the error below:

m_els = cm.load(mc.fname_tot_rig, is3D=True).play()
border_to_0 = 0 if mc.border_nan is 'copy' else mc.border_to_0 
    # maximum shift to be used for trimming against NaNs

Error Log:

error                                     Traceback (most recent call last)
<ipython-input-9-605c1b0cdc28> in <module>
----> 1 m_els = cm.load(mc.fname_tot_rig, is3D=True).play()
      2 border_to_0 = 0 if mc.border_nan is 'copy' else mc.border_to_0
      3     # maximum shift to be used for trimming against NaNs

~/CaImAn/caiman/base/movies.py in play(self, gain, fr, magnification, offset, interpolation, backend, do_loop, bord_px, q_max, q_min, plot_text, save_movie, opencv_codec, movie_name)
   1356                                     thickness=1)
   1357 
-> 1358                     cv2.imshow('frame', frame)
   1359                     if save_movie:
   1360                         if frame.ndim < 3:

error: OpenCV(4.2.0) ../modules/imgproc/src/color.simd_helpers.hpp:92: error: (-2:Unspecified error) in function 'cv::impl::(anonymous namespace)::CvtHelper<cv::impl::(anonymous namespace)::Set<3, 4, -1>, cv::impl::(anonymous namespace)::Set<3, 4, -1>, cv::impl::(anonymous namespace)::Set<0, 2, 5>, cv::impl::(anonymous namespace)::NONE>::CvtHelper(cv::InputArray, cv::OutputArray, int) [VScn = cv::impl::(anonymous namespace)::Set<3, 4, -1>, VDcn = cv::impl::(anonymous namespace)::Set<3, 4, -1>, VDepth = cv::impl::(anonymous namespace)::Set<0, 2, 5>, sizePolicy = cv::impl::(anonymous namespace)::NONE]'
> Invalid number of channels in input image:
>     'VScn::contains(scn)'
> where
>     'scn' is 350

Also, did the last question in my previous post make sense? Without MC, I can run CNMF and get nice components that incorporate the z-layers but the activity/pixel information from the bottom and top z-layers seem to be ignored.

This is the question I am referring to:

"Do you know how I can copy/add a false z-layer below and above the bottom and top z-layer or include the bottom and top z-layer in the CNMF initialization/extraction of components? When I run CNMF using the following commands, I don't think the bottom and top z-layer are being incorporated."

# INIT
cnm = cnmf.CNMF(n_processes, k=K, gSig=gSig, merge_thresh=merge_thresh, p=p, dview=dview)

# FIT
images = np.reshape(Yr.T, [T] + list(dims), order='F')    # reshape data in Python format (T x X x Y x Z)
cnm = cnm.fit(images)

# View the results
cnm.estimates.nb_view_components_3d(image_type='mean', dims=dims)

[marcusflisher]

Okay, I just realized my mistake on the playback of the motion corrected video. I forgot to take your advice and add [:,0] to the following line that I had in the above post:

m_els = cm.load(mc.fname_tot_rig, is3D=True)[:,0].play()

I no longer get an error, but the video that plays back has mixed up dimensions compared to the original file. It seems like it is playing back either the X or Y axis and the first z-layer over time, but I am not sure.

Is there an argument that specifies the order of dimensions when creating the motion correct object or does it require a specific order for 3D datasets? I can't find it in the documentation for ? MotionCorrect

[marcusflisher]

@j-friedrich Sorry if that was confusing. These are my questions:

What order of dimensions does motion correction take for 3D datasets?

Can I copy/add a false z-layer below and above the bottom and top z-layer because I think all of the pixels in these bottom/top frames are being treated as undefined?

[j-friedrich]

The order of dimensions is TXYZ.

I started updating the demo to include motion correction. Check it out here I see your issue with empty first and last plane, not sure why that happens. Until fixed, you can pad a z-layer below and above, by just concatenating to the movie array using numpy.pad, as I did in the linked file.

[marcusflisher]

Thanks so much for the ongoing support. Thank you for updating the demo notebook!

For the following code, how do I go about defining 'shifts', because I am getting the proceeding error?

plt.figure(figsize=(12,3))
for i, s in enumerate((mc.shifts_rig, shifts)):
    plt.subplot(1,2,i+1)
    for k in (0,1,2):
        plt.plot(np.array(s)[:,k], label=('x','y','z')[k])
    plt.legend()
    plt.title(('inferred shifts', 'true shifts')[i])
    plt.xlabel('frames')
    plt.ylabel('pixels')

Error Log:

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-35-67648d4d1748> in <module>
      1 plt.figure(figsize=(12,3))
----> 2 for i, s in enumerate((mc.shifts_rig, shifts)):
      3     plt.subplot(1,2,i+1)
      4     for k in (0,1,2):
      5         plt.plot(np.array(s)[:,k], label=('x','y','z')[k])

NameError: name 'shifts' is not defined

<Figure size 864x216 with 0 Axes>

[j-friedrich]

I added artificial motion to the simulated data. For real data you won't have the ground truth shifts to compare to, thus just skip this cell (or only plot mc.shifts_rig)

[marcusflisher]

If I plot mc.shifts_rig does the returned list represent the inferred pixel shift for Z, X, Y for each frame?

For example, when I run:

mc.shifts_rig

And it returns this list:

[(-0.0, 0.6, -0.2),
 (0.1, 0.7, -0.1),
 (-0.0, 0.7, -0.0),
 (-0.1, 0.7, 0.2),
 (-0.0, 0.8, -0.0),...

Also, I noticed that motion correction only works for my dataset if the dimensions are ordered TZXY as opposed to your recommendation of TXYZ.

[marcusflisher]

To answer my own question, yes I think this is what the coordinates represent. Here is a cool graph of the motion corrected results. Thank you.

[marcusflisher]

Should the output file size after concatenating two files be additive? (ie, if both my input files are 2.84 GB, should the output file be 2.84 + 2.84 = 5.68 GB?)

My output ends up being double the expected size.

[pgunn]

You shouldn't worry about it; the output may be upscaled or compressed differently than the original file.

[j-friedrich]

Addressed in the recent 1.8.5 release.