Demo 01 — V1, standard multi-plane 2P (TC030)

A 7-plane two-photon recording of mouse visual cortex, run stage by stage.

This is the demo to read first: it needs no parameter overrides at all beyond the acquisition geometry, so what you see below is Suite3D on its shipped defaults.

Data: <data-root>/v1/raw/*.tif from the figshare archive (10 tifs, 21.1 GB). Expect ~845 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('v1')

# `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, 'v1')
print(f'{len(tifs)} tifs')

job = Job(OUT_DIR, 'demo-v1', 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.

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