%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 extract_batch, get_params, find_tifs
from suite3d.job import Job

params = get_params('lbm')

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

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

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()

# All three are estimated during the init pass.
print(f"cavity crosstalk coefficient : {summary['crosstalk_coeff']:.3f}")
print(f"strip fusion shift (px)      : {summary['fuse_shift']}")
print(f"cavity_size                  : {params['cavity_size']} planes per cavity")
print(f"planes (acquisition order)   : {list(params['planes'])}")

fig, axs = plt.subplots(1, 2, figsize=(11, 3.2))
axs[0].plot(np.asarray(summary['crosstalk_planes']), 'o-')
axs[0].set_xlabel('plane pair')
axs[0].set_ylabel('estimated crosstalk')
axs[0].set_title(f"per-plane crosstalk (coeff = {summary['crosstalk_coeff']:.3f})")

ps = np.asarray(summary['plane_shifts'])
axs[1].plot(ps[:, 0], label='y')
axs[1].plot(ps[:, 1], label='x')
axs[1].set_xlabel('plane')
axs[1].set_ylabel('shift (px)')
axs[1].legend()
axs[1].set_title('plane-to-plane shifts')
plt.tight_layout()
plt.show()

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()

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()

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()

job.compute_npil_masks()

# Trace extraction sets this pipeline's peak memory. Suite3D's default of
# 500 volumes per batch costs ~56 GiB on a 22-plane LBM volume; snapping the
# batch to whole on-disk chunks is both cheaper and faster. See
# common/pipeline.py:extract_batch.
job.extract_and_deconvolve(batchsize_frames=extract_batch(job))

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()

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()

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

Demo 02 — Light-beads microscopy, ~40k neurons (SS004)¶

1. Initialization¶

What makes this recording LBM¶

2. Registration¶

3. Correlation map¶

4. Segmentation¶

5. Neuropil masks, traces, deconvolution¶

6. Browse the result¶