Demo 03 — densely packed CA1, standard 2P (ATL020)

A 4-plane hippocampal recording. Densely packed somata make this the hardest of the three to segment, and the shallow stack makes it the most interesting to register.

The one override worth understanding is apply_z_shift=False. Suite3D measures rigid motion in 3D, including an axial component. On a 4-plane stack the axial search window is only +/-2 planes, and 6.6% of this recording's volumes peg at that cap -- the estimate saturates instead of converging. Applying it would translate those frames by 62% of the volume depth. So z is still measured (3d_reg=True), just not applied.

That is not the same as turning off 3D registration. 3d_reg stays on.

Data: <data-root>/hippocampus/raw/*.tif (10 tifs, 21.1 GB). Expect ~1,347 ROIs.

%matplotlib inline
# inline PNGs, never ipywidgets: nbconvert renders a stateless widget
# as a silently empty div, and these pages go on the public web.

import sys, os
sys.path.insert(0, os.path.abspath('..'))

import numpy as np
import matplotlib.pyplot as plt

from common import get_params, find_tifs
from suite3d.job import Job

params = get_params('hippocampus')

# `get_params` asserts that fs is a plausible *volume* rate.
# Never use io.get_vol_rate() -- it returns the per-plane rate.
nz = len(params['planes'])
print(f"{nz} planes, fs = {params['fs']:.4f} Hz (volume rate), "
      f"voxel {params['voxel_size_um']} um (dz, dy, dx)")

# Everything NOT listed here comes from suite3d.default_params.
for k, v in params.items():
    print(f'  {k:28s} {v}')

# Point SUITE3D_DEMO_DATA at your figshare download, e.g.
#   export SUITE3D_DEMO_DATA=~/figshare
# Absolute paths are kept out of the notebook source on purpose: the
# executed outputs are published, and must not leak a filesystem layout.
DATA_ROOT = os.environ.get('SUITE3D_DEMO_DATA', './data')
OUT_DIR   = os.environ.get('SUITE3D_DEMO_OUT', './results')

tifs = find_tifs(DATA_ROOT, 'hippocampus')
print(f'{len(tifs)} tifs')

job = Job(OUT_DIR, 'demo-hippocampus', tifs=tifs, params=params,
          create=True, overwrite=True, verbosity=1)

1. Initialization

Reads the first few tifs, builds a 3D mean image, and estimates the plane-to-plane shifts. This is the reference everything else registers to.

job.run_init_pass()
summary = job.load_summary()
ref = summary['ref_img_3d']            # (nz, ny, nx)
print('reference volume:', ref.shape)

nz = ref.shape[0]
fig, axs = plt.subplots(1, nz, figsize=(3 * nz, 3.4))
for z, ax in enumerate(np.atleast_1d(axs)):
    ax.imshow(ref[z], cmap='gray',
              vmin=np.percentile(ref[z], 1), vmax=np.percentile(ref[z], 99.5))
    ax.set_title(f'z = {z}'); ax.axis('off')
fig.suptitle('Reference volume (mean image per plane)')
plt.tight_layout(); plt.show()

2. Registration

3D rigid + non-rigid alignment of every volume to the reference, by FFT phase correlation. The rigid offsets below should look like slow drift plus small jitter; a trace that sits pinned at a constant value has hit the edge of its search window.

job.register()

reg = job.load_registration_results()
print('registration keys:', sorted(reg.keys()))

# Rigid registration offsets over time.
#
# `int_shift` is the integer voxel shift that was actually applied to each
# volume; `sub_pixel_shifts` is the raw phase-correlation estimate before
# rounding. Both are (n_volumes, 3) with columns ordered (z, y, x)
# -- see reg_3d.py.
#
# This demo sets apply_z_shift=False, so the z column below is *measured*
# but never applied to the movie.

shifts = np.asarray(reg['int_shift'])

fig, ax = plt.subplots(figsize=(9, 3))
for i, axis_name in enumerate(('z', 'y', 'x')):
    ax.plot(shifts[:, i], lw=0.8, label='%s shift' % axis_name)
ax.set_xlabel('volume')
ax.set_ylabel('shift (voxels)')
ax.legend()
ax.set_title('Rigid registration offsets')
plt.tight_layout()
plt.show()

3. Correlation map

The heart of Suite3D. Each voxel is scored by how correlated it is with its neighbours over time, after filtering at a cell scale and a neuropil scale. Cells light up as blobs; vessels and neuropil do not.

This is computed on the whole volume at once, which is what lets a cell spanning several planes be found as one object.

job.calculate_corr_map()
cm = job.load_corr_map_results()
vmap, mean_img = cm['vmap'], cm['mean_img']
print('corr map:', vmap.shape)

nz = vmap.shape[0]
fig, axs = plt.subplots(2, nz, figsize=(3 * nz, 6.5), squeeze=False)
for z in range(nz):
    axs[0, z].imshow(mean_img[z], cmap='gray',
                     vmax=np.percentile(mean_img[z], 99.5))
    axs[0, z].set_title(f'mean, z={z}'); axs[0, z].axis('off')
    axs[1, z].imshow(vmap[z], cmap='magma',
                     vmax=np.percentile(vmap[z], 99.5))
    axs[1, z].set_title(f'corr map, z={z}'); axs[1, z].axis('off')
plt.tight_layout(); plt.show()

4. Segmentation

Peaks in the correlation map seed ROIs; each seed is grown into a 3D footprint, then subtracted so the next iteration finds the next cell. Voxels join an ROI when they pass vox_snr_thresh.

job.segment_rois()
stats = job.load_segmentation_results(to_load=['stats'])
print(f'{len(stats)} ROIs detected')

# How far do ROIs extend in z? Multi-plane cells are the point of 3D detection.
zspan = np.array([s['coords'][0].max() - s['coords'][0].min() + 1 for s in stats])
npix  = np.array([len(s['lam']) for s in stats])

fig, axs = plt.subplots(1, 2, figsize=(10, 3.2))
axs[0].hist(zspan, bins=np.arange(0.5, zspan.max() + 1.5))
axs[0].set_xlabel('z-planes spanned'); axs[0].set_ylabel('# ROIs')
axs[1].hist(npix, bins=40)
axs[1].set_xlabel('voxels per ROI'); axs[1].set_ylabel('# ROIs')
plt.tight_layout(); plt.show()
print(f'median z-span {np.median(zspan):.0f} planes, '
      f'{(zspan > 1).mean() * 100:.0f}% of ROIs span more than one plane')

# ROI footprints painted onto the correlation map, plane by plane.
lbl = np.zeros(vmap.shape, dtype=int)
for i, s in enumerate(stats):
    zc, yc, xc = s['coords']
    lbl[zc, yc, xc] = i + 1

nz = vmap.shape[0]
fig, axs = plt.subplots(1, nz, figsize=(3 * nz, 3.4), squeeze=False)
for z in range(nz):
    axs[0, z].imshow(vmap[z], cmap='gray', vmax=np.percentile(vmap[z], 99.5))
    m = np.ma.masked_where(lbl[z] == 0, lbl[z] % 20 + 1)
    axs[0, z].imshow(m, cmap='tab20', alpha=0.65, interpolation='nearest')
    axs[0, z].set_title(f'z={z}  ({(lbl[z] > 0).sum()} vox)'); axs[0, z].axis('off')
plt.tight_layout(); plt.show()

5. Neuropil masks, traces, deconvolution

An annular neuropil mask is built around each ROI; F is the raw trace, Fneu the surrounding neuropil, and the corrected trace is F - npil_coeff * Fneu. spks is the deconvolved estimate.

job.compute_npil_masks()
job.extract_and_deconvolve()

res = job.load_segmentation_results(to_load=['F', 'Fneu', 'spks'])
F, Fneu, spks = res['F'], res['Fneu'], res['spks']
print('traces:', F.shape, '(n_rois, n_volumes)')

fs = params['fs']                      # volume rate -> seconds
t = np.arange(F.shape[1]) / fs

# The most active ROIs, by deconvolved activity.
order = np.argsort(spks.sum(1))[::-1][:5]
fig, axs = plt.subplots(len(order), 1, figsize=(11, 2 * len(order)), sharex=True)
for ax, i in zip(np.atleast_1d(axs), order):
    ax.plot(t, F[i] - 0.7 * Fneu[i], lw=0.7, label='F - 0.7*Fneu')
    ax.plot(t, spks[i], lw=0.7, alpha=0.7, label='spks')
    ax.set_ylabel(f'ROI {i}'); ax.legend(loc='upper right', fontsize=7)
np.atleast_1d(axs)[-1].set_xlabel('time (s)')
plt.tight_layout(); plt.show()

6. Browse the result

Export, and open the portable HTML viewer (or napari).

job.export_results(OUT_DIR, make_viewer=True) \
    if hasattr(job, 'make_html_viewer') else job.export_results(OUT_DIR)

# Or, in napari:
#   from suite3d.ui import create_napari_ui   # NB: ui.py, not curation.py
#   import napari
#   outputs = job.load_segmentation_results()
#   outputs.update(job.load_corr_map_results())
#   create_napari_ui(outputs, scale=params['voxel_size_um'])
#   napari.run()

Equivalent one-liner:

python run_pipeline.py --data-root /path/to/figshare --out-dir ./results