SpikeInterface / spikeinterface

A Python-based module for creating flexible and robust spike sorting pipelines.
https://spikeinterface.readthedocs.io
MIT License
524 stars 187 forks source link

InterpolateMotionRecording returns zeros for rigid motion #1893

Open TomBugnon opened 1 year ago

TomBugnon commented 1 year ago

That one is nasty

Loading a raw sglx file:

from spikeinterface.full import read_spikeglx
from spikeinterface.sortingcomponents.motion_interpolation import InterpolateMotionRecording

sglx_path = "/Volumes/ceph-tononi/npx_archive/CNPIX2-Segundo/1-21-2020/SpikeGLX/1-21-2020_g0/1-21-2020_g0_imec0"
rec = read_spikeglx(sglx_path, stream_id="imec0.ap").select_segments(0).frame_slice(start_frame=0, end_frame=30000 * 10)
# FrameSliceRecording: 384 channels - 30.0kHz - 1 segments - 300,000 samples - 10.00s - int16 dtype 
#                     219.73 MiB

print(rec.get_traces(0, 0, 10))

memmap([[ 10, -19, -13, ..., 21, -3, 24], [ -5, -10, -21, ..., 26, 4, 21], [ 5, -12, -24, ..., 32, 6, 14], ..., [-15, 0, -18, ..., 39, 16, 18], [-18, -4, -25, ..., 36, 8, 16], [-13, -13, -21, ..., 28, 2, 18]], dtype=int16)

Apply any rigid motion interpolation and you end up with zeros

# Dummy rigid motion
motion = np.zeros((10, 1))
spatial_bins=np.array([3760])
temporal_bins=np.arange(10)

rec_rigid = InterpolateMotionRecording(
    rec,
    motion,
    temporal_bins,
    spatial_bins,
) # Same thing with method="idw"
# InterpolateMotionRecording: 384 channels - 30.0kHz - 1 segments - 300,000 samples - 10.00s 
                            int16 dtype - 219.73 MiB

print(rec_rigid.get_traces(0, 0, 10))

/home/tbugnon/miniconda3/envs/ecephys_dev/lib/python3.10/site-packages/scipy/interpolate/_interpolate.py:698: RuntimeWarning: invalid value encountered in divide slope = (y_hi - y_lo) / (x_hi - x_lo)[:, None]

array([[0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0], ..., [0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0]], dtype=int16)

While it looks fine for non-rigid

# Dummy non-rigid motion
motion = np.zeros((10, 2))
spatial_bins=np.array([1000, 4000])
temporal_bins=np.arange(10)

rec_nonrigid = InterpolateMotionRecording(
    rec,
    motion,
    temporal_bins,
    spatial_bins,
)

rec_nonrigid.get_traces(0, 0, 10)

array([[ 9, -18, -12, ..., 21, -2, 23], [ -5, -9, -20, ..., 25, 4, 21], [ 4, -11, -23, ..., 31, 6, 14], ..., [-14, 0, -17, ..., 38, 15, 18], [-18, -3, -24, ..., 35, 8, 16], [-13, -12, -20, ..., 28, 2, 18]], dtype=int16)

samuelgarcia commented 1 year ago

Salut Tom. Let this open I will try to have a look. Note that the bin have to be the bin center here your temporal_bins is the apparently the left.

TomBugnon commented 1 year ago

Yes I just made a quick illustration, I get the same behavior with motion/spatial_bins/temporal_bins actually returned by SI functions

samuelgarcia commented 1 year ago

Is 3760 the middle of your probe ?

TomBugnon commented 1 year ago

Actually the real value of spatial_bins returned for rigid motion estimation is array([3830.]), but that gives the same results in the dummy examples I provide