diff --git a/.agents/skills b/.agents/skills new file mode 120000 index 00000000..4ca0ec66 --- /dev/null +++ b/.agents/skills @@ -0,0 +1 @@ +../.skills \ No newline at end of file diff --git a/.claude/skills b/.claude/skills new file mode 120000 index 00000000..4ca0ec66 --- /dev/null +++ b/.claude/skills @@ -0,0 +1 @@ +../.skills \ No newline at end of file diff --git a/.gitignore b/.gitignore index c01684c4..f70484a9 100644 --- a/.gitignore +++ b/.gitignore @@ -156,9 +156,12 @@ waveorder/_version.py /examples/data_temp/* /logs/* scripts/ +!.skills/**/scripts/ +!.skills/**/scripts/** # benchmark results runs/ -# Claude Code session state -.claude/ +# Claude Code session state (but keep shared project skills symlink) +.claude/* +!.claude/skills diff --git a/.skills/visualize-otf-from-config/SKILL.md b/.skills/visualize-otf-from-config/SKILL.md new file mode 100644 index 00000000..f4251c73 --- /dev/null +++ b/.skills/visualize-otf-from-config/SKILL.md @@ -0,0 +1,49 @@ +--- +name: visualize-otf-from-config +description: Compute a waveorder transfer function (OTF) from a reconstruction config and open it in napari. Use when someone wants to inspect/visualize the transfer function or OTF for a config, check how NA / wavelength / pixel-size choices shape the transfer function, or debug a config before reconstructing. Runs `wo compute-tf` then `wo view`. +--- + +# Visualize the OTF from a config + +Compute the transfer function for a `waveorder` reconstruction config and open it +in napari. Exactly two CLI calls. + +## Inputs + +- **config** — a reconstruction config YAML (the same kind used by `wo rec`; e.g. + a `phase:` or `fluorescence:` config). This defines the optics (NA, wavelength, + pixel sizes, RI, tilt, …) that shape the transfer function. +- **input position** — a path to one OME-Zarr position (`.../row/col/fov`). + `wo compute-tf` needs it only to read the **array shape**; the values are not + used, so any dataset with the right ZYX shape works. + +If the user gave a `/visualize-otf-from-config ` argument, use it. +Otherwise ask for the config path and the input position path. + +## The two calls + +```bash +# 1. Compute the transfer function +wo compute-tf -i /0/0/0 -c -o ./otf.zarr + +# 2. Open it in napari (backgrounded so it doesn't block) +wo view ./otf.zarr & +``` + +Notes: +- The subcommand is `compute-tf` (alias for `compute-transfer-function`) — there + is no `calc-tf`. +- `wo view` auto-detects a transfer function (it has `settings` in its zarr + attrs) and shows real/imag parts with a diverging `bwr` colormap, `ifftshift`-ed + so DC is centered. **3D** reconstructions store the transfer function directly. + **2D** reconstructions store a singular system (`U`, `S`, `Vh`); `wo view` + reconstructs the transfer function from the SVD (`H = U @ diag(S) @ Vh`) and + shows that (e.g. the absorption and phase transfer functions for 2D phase) + rather than the raw singular vectors. +- `wo view` starts napari's event loop and **blocks the terminal until closed**. + When you run it, launch it in the background (`&`, or the Bash tool's + `run_in_background: true`), tell the user the window is open, and let them + inspect it — do not render a static PNG. +- If `wo` isn't on PATH, prefix both calls with `uv run` from the repo. + +That's it — report the `./otf.zarr` path and that napari is open. diff --git a/.skills/waveorder-tutorial/SKILL.md b/.skills/waveorder-tutorial/SKILL.md new file mode 100644 index 00000000..492ea30c --- /dev/null +++ b/.skills/waveorder-tutorial/SKILL.md @@ -0,0 +1,236 @@ +--- +name: waveorder-tutorial +description: Interactively walk a user through their first waveorder reconstruction from a raw imaging file via the CLI. Use when someone has a microscopy image (TIFF, zarr, or array) and wants to reconstruct phase (brightfield) or deconvolved fluorescence, or asks to "try waveorder", "reconstruct my data", or "set up a waveorder config". Guides environment setup, OME-Zarr conversion, config authoring, viewing in napari, and parameter sweeps. +--- + +# waveorder reconstruction tutorial + +Guide the user, step by step, from a raw imaging file to a first `waveorder` +reconstruction they can inspect in napari. This is an **interactive** skill: +ask the questions, wait for answers, run the commands, and interpret the +results together. Do not dump every question at once — go one stage at a time +and adapt to what the user reports. + +Supporting material in this skill: +- `references/conversion.md` — how to turn TIFF / zarr / array into OME-Zarr +- `references/config_reference.md` — every config knob explained +- `references/example_qpi_2d.yml` — ready config for the QPI-defocus sample +- `scripts/get_example_data.py` — fetch the QPI-defocus sample OME-Zarr +- `scripts/to_omezarr.py` — converter (TIFF/npy/array → viewable OME-Zarr) +- `scripts/sweep.py` — regularization / parameter sweep + side-by-side view + +## The `wo` command + +The CLI entrypoint is `wo` (alias for `waveorder`). Key subcommands: +- `wo rec -i /*/*/* -c -o ` — reconstruct +- `wo view [ ...]` — open datasets / transfer functions in napari +- `wo sim -c -o ` — simulate phantom + measurement + +Input paths use the glob `input.zarr/*/*/*` to expand HCS plate positions. +If plain `wo` isn't on PATH, prefix commands with `uv run` from the repo. + +## Viewing convention — always napari, never PNGs + +**Always show data and reconstructions in napari via `wo view`. Never render +static PNGs / matplotlib figures / screenshots to inspect or present results.** +napari is interactive — the *user* scrolls Z, adjusts contrast, and toggles +layers, which is exactly what's needed to judge focus, density sign, noise, and +ringing. A flat PNG can't support that and would defeat the point. + +`wo view` starts napari's event loop and **blocks the terminal until the window +is closed**. So when *you* (the assistant) launch it, run it in the background +so the session isn't stuck: + +```bash +wo view ./data.zarr ./recon.zarr & # background; napari opens for the user +``` + +(With the Bash tool, use `run_in_background: true`.) Then tell the user the +napari window is open, describe what to look for, and **ask them what they +observe** — do not try to see the result yourself. Reconstructions come back +squeezed so 2D results (single Z) display correctly, and grid view is enabled so +input and output sit side by side. + +--- + +## Step 0 — Check the environment + +Verify the CLI and napari before anything else: + +```bash +wo --help # is the CLI installed? +wo view 2>&1 | head -1 # will error, but confirms import +python -c "import napari; print(napari.__version__)" +``` + +- If `wo` is missing, install into the user's env with uv: + `uv pip install "waveorder[all]"` (the `all`/`visual` extra pulls in + `napari[pyqt6]` and `napari-ome-zarr`). If they're developing the repo: + `uv pip install -e ".[all]"` from the waveorder checkout. +- If napari imports but no Qt backend is present, `uv pip install "napari[pyqt6]"`. +- Confirm `wo --help` lists `reconstruct`, `view`, `simulate` before moving on. +- If the user has no data of their own to practice on, note that + `scripts/get_example_data.py` can fetch a zenodo sample OME-Zarr (Step 1). + +--- + +## Step 1 — Get the data into a viewable OME-Zarr + +Ask: **"What form is your data in?"** +1. **an OME-Zarr** (already in `waveorder`'s native format) — **preferred** +2. a TIFF / OME-TIFF (single file or a folder of tiles) +3. an in-memory / `.npy` array or a plain (non-OME) zarr + +The key thing `waveorder` needs is a **TCZYX** OME-Zarr with named channels and +a correct pixel scale. + +**No data of their own?** Offer the QPI-from-defocus sample (the dataset from +`docs/examples/demos/QPI_defocus`). It's already an OME-Zarr — a single `BF` +brightfield defocus stack — so it skips conversion and is set up for a **2D +phase-from-defocus** reconstruction: +```bash +python scripts/get_example_data.py # downloads + prints paths & commands +``` +Use the printed `raw_data.zarr` path in place of `./data.zarr` below. This +example's answers are: contrast = brightfield→phase, dimension = **2D**, and the +optical parameters are pre-filled in `references/example_qpi_2d.yml` +(BF channel, yx 0.325 µm, z 2.0 µm, n 1.0, λ 0.532 µm, NA_det 0.55, NA_ill 0.4). + +**Case 1 — OME-Zarr:** no conversion. Ask the user to paste the path, then just +confirm the axis order and channel names: +`wo view ` — the stack should focus/defocus as you scroll Z. + +**Case 2 — TIFF:** ask the user to **paste the path to their TIFF** (a single +file, or a folder of tiles). Then convert with `scripts/to_omezarr.py`: +```bash +python scripts/to_omezarr.py ./data.zarr \ + --channel-name Brightfield --axes ZYX \ + --yx-pixel-size 0.1 --z-pixel-size 0.25 +``` + +**Case 3 — array / plain zarr:** save to `.npy` if needed, then use +`scripts/to_omezarr.py` with `--axes` describing the dimension order. + +See `references/conversion.md` for per-case flags. Then always: +- Confirm the result visually: `wo view ./data.zarr` — check Z focus and that + channel names look right. +- Note the array shape and pixel sizes; you'll reuse them in the config. + +--- + +## Step 2 — Choose contrast type and dimension + +Ask two questions: + +**A. Contrast type?** +- **Brightfield → phase** (label-free density). Use the `phase:` config block. +- **Fluorescence → deconvolution** (denoise + sharpen). Use `fluorescence:`. +- (Polarization/birefringence is out of scope for a first tutorial; mention it + exists but steer to phase or fluorescence.) + +**B. Reconstruction dimension — 2D or 3D?** +- **2D** works best when the object is **thin compared to the depth of focus** — + i.e., the whole object appears to go in and out of focus at the same time as + you scroll Z. Output is a single in-focus plane. +- **3D** for thick samples where different depths focus at different Z. Output + is a volume. Start here if unsure and the sample is clearly volumetric. + +Record the answers → this picks the config template and +`reconstruction_dimension: 2` or `3`. + +--- + +## Step 3 — Fill in the configuration + +Copy the matching template from the waveorder repo +(`docs/examples/cli/configs/{phase,fluorescence}_{2d,3d}.yml`) or write one from +`references/config_reference.md`. Ask the user for each optical parameter and +fill it in. Always set `input_channel_names` to the exact channel name in +their OME-Zarr. + +Ask for: +- **`yx_pixel_size`** (µm) and **`z_pixel_size`** (µm) — from the acquisition; + reuse what Step 1 recorded. +- **`numerical_aperture_detection`** — the objective NA. +- **`index_of_refraction_media`** — 1.0 air, ~1.33 water, ~1.47 oil, etc. +- **wavelength** — `wavelength_illumination` (phase) or `wavelength_emission` + (fluorescence), in µm. +- Phase only: **`numerical_aperture_illumination`** — the condenser NA (must be + ≤ `index_of_refraction_media`). +- Fluorescence only: **`confocal_pinhole_diameter`** — `null` for widefield. + +**Flag `invert_phase_contrast` explicitly (phase only).** This is the knob most +likely to be wrong on the first try, and it's hard to set from theory — it +depends on the microscope's contrast convention. Plan to try **both `true` and +`false`** and pick the one with the correct density sign: +- **less-dense regions** (nuclei, vacuoles, lumen) should appear **darker than + background** +- **more-dense regions** (cytoplasm, membranes, dense organelles) should appear + **brighter than background** + +Start `tilt_angle_zenith` and `tilt_angle_azimuth` at **0.0** (see Step 4). + +Save the config next to the data, e.g. `./config.yml`. + +--- + +## Step 4 — Run it and inspect in napari + +```bash +wo rec -i ./data.zarr/*/*/* -c ./config.yml -o ./recon.zarr +wo view ./data.zarr ./recon.zarr & # napari, backgrounded (see convention above) +``` + +Open the input and reconstruction together in napari (never a PNG). Tell the +user the window is open, then review the result **with them** — ask what they +see. Look for: +- **Density sign (phase):** apply the nuclei-dark / cytoplasm-bright test above. + If inverted, flip `invert_phase_contrast` and re-run. This comes first — + everything else is easier to judge once the sign is right. +- **Noise / graininess:** speckly, high-frequency texture → regularization is + too low. Go to Step 5. +- **Ringing / halos:** dark/bright overshoot around edges → regularization too + low, or try the `TV` algorithm. +- **Over-smoothing / loss of detail:** regularization too high → Step 5. + +**Tilt:** keep `tilt_angle_zenith = 0.0` unless the user can *explicitly see* +directional "gradient- or shadow-casting" contrast in the in-focus data — an +image that looks lit from one side, like DIC. Only then add a small zenith tilt +and set the azimuth to point along the shadow direction. Don't add tilt to +chase artifacts that aren't shadow-cast. + +--- + +## Step 5 — Parameter sweeps + +When the user reports noise, ringing, or blur, sweep regularization instead of +guessing. `scripts/sweep.py` writes one config per value, reconstructs each, and +opens all results side by side in napari (grid view) for comparison — again, no +PNGs: + +```bash +python scripts/sweep.py -i ./data.zarr -c ./config.yml \ + --param regularization_strength --values 1e-4 1e-3 1e-2 1e-1 +``` + +The script launches napari at the end, which blocks. When *you* run it, use +`run_in_background: true`; or pass `--no-view` and open the results yourself in a +backgrounded `wo view ...`. Then ask the user which value looks best. + +Guidance: +- **regularization_strength** is the main dial. Sweep over decades + (e.g. 1e-4 → 1e-1). Higher = smoother/less noise but more blur; lower = + sharper but noisier with more ringing. +- If ringing persists at a good noise level, sweep + `reconstruction_algorithm` between `Tikhonov` and `TV` (TV suppresses ringing + while keeping edges, at higher compute cost; tune `TV_rho_strength` / + `TV_iterations`). +- Pick the value at the "elbow" — the smallest regularization that removes the + objectionable noise/ringing without visibly softening real structures. +- Other sweepable knobs for fine-tuning: `z_focus_offset`, + `numerical_aperture_illumination`. (waveorder can also *optimize* some of + these automatically — mention `wo rec` with an `lr` key / the `optimization:` + block if the user wants that later.) + +Wrap up by saving the final config and reconstruction, and summarizing the +chosen parameters for the user. diff --git a/.skills/waveorder-tutorial/references/config_reference.md b/.skills/waveorder-tutorial/references/config_reference.md new file mode 100644 index 00000000..7d20241f --- /dev/null +++ b/.skills/waveorder-tutorial/references/config_reference.md @@ -0,0 +1,90 @@ +# waveorder config reference + +A reconstruction config is a YAML file. Top-level keys plus one of the +`phase:` / `fluorescence:` blocks. Copy a starting template from the repo: +`docs/examples/cli/configs/{phase,fluorescence}_{2d,3d}.yml`. + +## Top level + +```yaml +input_channel_names: # list; must match channel names in the OME-Zarr +- Brightfield +time_indices: all # "all", an int, or a list of ints +reconstruction_dimension: 3 # 2 for thin samples, 3 for thick +device: null # null=cpu, "auto", "cuda:0", "mps" +``` + +Exactly one of `phase:` or `fluorescence:` for a first reconstruction +(they cannot coexist in one file; birefringence is separate). + +## `phase:` block + +```yaml +phase: + transfer_function: + wavelength_illumination: 0.532 # µm + yx_pixel_size: 0.1 # µm + z_pixel_size: 0.25 # µm + z_padding: 0 # z slices padded for axial boundaries + z_focus_offset: 0 # (sweepable) focus slice offset + index_of_refraction_media: 1.3 # 1.0 air / 1.33 water / 1.47 oil + numerical_aperture_detection: 1.2 # objective NA + numerical_aperture_illumination: 0.9 # condenser NA (≤ index_of_refraction_media) + tilt_angle_zenith: 0.0 # radians; keep 0 unless shadow-cast contrast + tilt_angle_azimuth: 0.0 # radians; shadow direction + invert_phase_contrast: false # TRY BOTH — see density-sign test below + apply_inverse: + reconstruction_algorithm: Tikhonov # Tikhonov (fast) or TV (edge-preserving) + regularization_strength: 0.001 # main noise/sharpness dial + TV_rho_strength: 0.001 # TV only: ADMM rho + TV_iterations: 1 # TV only: ADMM iterations +``` + +### `invert_phase_contrast` +The hardest knob. It sets the sign convention of density contrast and depends on +the microscope. Reconstruct with both `true` and `false`, then pick the one +where: +- **less-dense** regions (nuclei, vacuoles, lumen) are **darker** than background +- **more-dense** regions (cytoplasm, membranes) are **brighter** than background + +### tilt angles +Leave at `0.0`. Only add a nonzero `tilt_angle_zenith` if the in-focus raw data +visibly shows one-sided, DIC-like "gradient / shadow-casting" contrast; set +`tilt_angle_azimuth` to the shadow direction. + +## `fluorescence:` block + +Same optical keys, but with emission instead of illumination: + +```yaml +fluorescence: + transfer_function: + yx_pixel_size: 0.1 + z_pixel_size: 0.25 + z_padding: 0 + z_focus_offset: 0 + index_of_refraction_media: 1.3 + numerical_aperture_detection: 1.2 + tilt_angle_zenith: 0.0 + tilt_angle_azimuth: 0.0 + wavelength_emission: 0.532 # µm + confocal_pinhole_diameter: null # null = widefield; else diameter (AU/µm) + apply_inverse: + reconstruction_algorithm: Tikhonov + regularization_strength: 0.001 + TV_rho_strength: 0.001 + TV_iterations: 1 +``` + +## Which knobs to sweep + +| symptom | knob | +|-----------------------------|------------------------------------------| +| grainy / noisy | ↑ `regularization_strength` | +| over-smoothed / lost detail | ↓ `regularization_strength` | +| ringing / halos at edges | try `TV`; tune `TV_rho_strength` | +| wrong density sign (phase) | flip `invert_phase_contrast` | +| slightly defocused output | sweep `z_focus_offset` | + +`regularization_strength` is best swept over decades (1e-4 … 1e-1). Pick the +"elbow": smallest value that removes noise/ringing without softening structure. diff --git a/.skills/waveorder-tutorial/references/conversion.md b/.skills/waveorder-tutorial/references/conversion.md new file mode 100644 index 00000000..2574eca0 --- /dev/null +++ b/.skills/waveorder-tutorial/references/conversion.md @@ -0,0 +1,88 @@ +# Converting data to OME-Zarr + +`waveorder` reads **OME-Zarr** (HCS layout) with a **TCZYX** array, named +channels, and a pixel-size scale transform. If the data is already OME-Zarr, +no conversion is needed. `scripts/to_omezarr.py` handles the other raw inputs. +Below are the cases and the flags to use. + +## Example data (no data of your own) + +`scripts/get_example_data.py` downloads the **QPI-from-defocus** sample +(`recOrder_session.zip` from zenodo, the dataset behind +`docs/examples/demos/QPI_defocus`). It is already an **OME-Zarr**, so use its +path directly and skip conversion: + +```bash +python scripts/get_example_data.py # ~10 MB +``` + +The stack is a single `BF` brightfield channel with 11 defocus slices +(yx = 0.325 µm, z = 2.0 µm), imaged in air. Reconstruct **2D phase** from it +with `references/example_qpi_2d.yml`: + +```bash +wo rec -i /0/0/0 -c references/example_qpi_2d.yml -o ./qpi_2d_recon.zarr +wo view ./qpi_2d_recon.zarr +``` + +## The converter + +```bash +python scripts/to_omezarr.py \ + --channel-name \ + --axes \ + --yx-pixel-size \ + --z-pixel-size +``` + +- `` — a `.tif`/`.tiff`/`.ome.tif`, a `.npy`, or a plain zarr array dir. +- `--axes` — the axis order of the *input* array, any subset/order of `TCZYX` + (e.g. `ZYX`, `CZYX`, `YX`, `TZYX`). The script transposes/expands to TCZYX. +- `--channel-name` — one name per channel (repeat the flag for multiple). + Must match `input_channel_names` in the config later. +- Pixel sizes in micrometers; used for the scale transform (and thus napari's + physical axes and the reconstruction's optical model). + +Verify with `wo view ` afterwards. + +## Case 1 — OME-Zarr (preferred) + +Already in `waveorder`'s native format — **don't convert**. Ask the user to +paste the path. Confirm axis order and channel names with `wo view `; if +it opens and scrolls in Z correctly, use it directly. + +## Case 2 — TIFF / OME-TIFF + +Ask the user to paste the path to their TIFF (a single file or a folder of +tiles). Single 3D stack saved ZYX: +```bash +python scripts/to_omezarr.py ./data.zarr \ + --channel-name Brightfield --axes ZYX \ + --yx-pixel-size 0.1 --z-pixel-size 0.25 +``` +An OME-TIFF usually already encodes axes and scale — the script reads them when +present, and the flags override. For a **folder of tiles**, stitch or pick one +tile first; `wo tile-stitch` (`wo ts`) handles multi-position mosaics once each +tile is an OME-Zarr position. + +## Case 3 — in-memory / `.npy` array or plain (non-OME) zarr + +For a plain zarr array, point the converter at the array directory and pass +`--axes` describing its dimension order. For an in-memory array, save to `.npy` +first. + +Save to `.npy` first (`np.save("arr.npy", arr)`), then: +```bash +python scripts/to_omezarr.py arr.npy ./data.zarr \ + --channel-name GFP --axes CZYX \ + --yx-pixel-size 0.1 --z-pixel-size 0.25 +``` + +## Sanity checks + +- **Axis order** is the most common mistake — if the reconstruction looks + scrambled, re-check `--axes`. +- **Z scrolls through focus:** in `wo view`, scrolling the Z slider should move + the plane of focus through the sample. +- **Pixel sizes are in µm** and roughly match the objective (e.g. a 63×/1.4 + objective on a typical camera → ~0.1 µm yx). diff --git a/.skills/waveorder-tutorial/references/example_qpi_2d.yml b/.skills/waveorder-tutorial/references/example_qpi_2d.yml new file mode 100644 index 00000000..c1226b07 --- /dev/null +++ b/.skills/waveorder-tutorial/references/example_qpi_2d.yml @@ -0,0 +1,22 @@ +input_channel_names: # matches the 'BF' channel in the sample +- BF +time_indices: all +reconstruction_dimension: 2 # 2D phase from a defocus stack +phase: + transfer_function: + wavelength_illumination: 0.532 # µm + yx_pixel_size: 0.325 # µm (6.5 µm camera / 20x) + z_pixel_size: 2.0 # µm (defocus step) + z_padding: 0 + z_focus_offset: 0 + index_of_refraction_media: 1.0 # imaging in air + numerical_aperture_detection: 0.55 + numerical_aperture_illumination: 0.4 + tilt_angle_zenith: 0.0 # no shadow-cast contrast in this data + tilt_angle_azimuth: 0.0 + invert_phase_contrast: false # try both; pick correct density sign + apply_inverse: + reconstruction_algorithm: Tikhonov + regularization_strength: 0.01 + TV_rho_strength: 0.001 + TV_iterations: 1 diff --git a/.skills/waveorder-tutorial/scripts/get_example_data.py b/.skills/waveorder-tutorial/scripts/get_example_data.py new file mode 100644 index 00000000..704af1a0 --- /dev/null +++ b/.skills/waveorder-tutorial/scripts/get_example_data.py @@ -0,0 +1,56 @@ +"""Download the QPI-from-defocus sample dataset for the tutorial. + +Fetches ``recOrder_session.zip`` from zenodo (the dataset used by the +``docs/examples/demos/QPI_defocus`` demo) and prints the path to the brightfield +defocus stack. It is already an **OME-Zarr**, so it plugs straight into +``wo view`` / ``wo rec`` with no conversion. + +The stack is a single ``BF`` channel, 11 defocus slices (yx = 0.325 µm, +z = 2.0 µm) — ideal for a **2D phase-from-defocus** reconstruction. A ready-made +config lives at ``references/example_qpi_2d.yml``. + +Examples +-------- + python get_example_data.py +""" + +import argparse +import io +import zipfile +from pathlib import Path + +URL = "https://zenodo.org/record/8386856/files/recOrder_session.zip" + + +def main() -> None: + argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter).parse_args() + + import requests + + data_dir = Path.home() / ".waveorder_tutorial_data" + data_dir.mkdir(exist_ok=True) + session_dir = data_dir / "recOrder_session" + + if not session_dir.exists(): + print("Downloading QPI-defocus sample from zenodo (~10 MB)...") + resp = requests.get(URL) + resp.raise_for_status() + with zipfile.ZipFile(io.BytesIO(resp.content)) as z: + z.extractall(session_dir) + + raw_path = session_dir / "recOrder_session" / "phase_snap_0" / "raw_data.zarr" + if not raw_path.exists(): + raise SystemExit(f"expected {raw_path} after unzip, not found") + + config = Path(__file__).resolve().parent.parent / "references" / "example_qpi_2d.yml" + + print(f"\nRaw OME-Zarr (BF defocus stack): {raw_path}") + print(f"\nInspect with: wo view {raw_path}") + print("Already an OME-Zarr — skip conversion (Step 1, Case 1); channel 'BF', dim 2D.") + print("\n2D phase-from-defocus reconstruction:") + print(f" wo rec -i {raw_path}/0/0/0 -c {config} -o ./qpi_2d_recon.zarr") + print(f" wo view {raw_path} ./qpi_2d_recon.zarr") + + +if __name__ == "__main__": + main() diff --git a/.skills/waveorder-tutorial/scripts/sweep.py b/.skills/waveorder-tutorial/scripts/sweep.py new file mode 100644 index 00000000..372669e7 --- /dev/null +++ b/.skills/waveorder-tutorial/scripts/sweep.py @@ -0,0 +1,101 @@ +"""Sweep one config parameter, reconstruct each value, and view side by side. + +Writes one config per value into a ``sweep/`` folder, runs ``wo rec`` for each, +then opens all reconstructions in napari with ``wo view`` for comparison. + +Examples +-------- + python sweep.py -i ./data.zarr -c ./config.yml \ + --param regularization_strength --values 1e-4 1e-3 1e-2 1e-1 + + python sweep.py -i ./data.zarr -c ./config.yml \ + --param reconstruction_algorithm --values Tikhonov TV --no-view +""" + +import argparse +import shutil +import subprocess +from pathlib import Path + +import yaml + + +def _wo_cmd() -> list[str]: + """Return the invocation prefix for the waveorder CLI.""" + if shutil.which("wo"): + return ["wo"] + return ["uv", "run", "wo"] + + +def _set_nested(cfg: dict, param: str, value) -> bool: + """Set ``param`` anywhere it appears under an ``apply_inverse``/``transfer_function`` block.""" + found = False + for block in ("phase", "fluorescence", "birefringence"): + sub = cfg.get(block) + if not isinstance(sub, dict): + continue + for section in ("apply_inverse", "transfer_function"): + sec = sub.get(section) + if isinstance(sec, dict) and param in sec: + sec[param] = value + found = True + return found + + +def _coerce(value: str): + """Best-effort str → number, leaving non-numeric strings as-is.""" + for cast in (int, float): + try: + return cast(value) + except ValueError: + continue + return value + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + p.add_argument("-i", "--input", required=True, type=Path, help="input OME-Zarr") + p.add_argument("-c", "--config", required=True, type=Path, help="base config .yml") + p.add_argument("--param", default="regularization_strength", help="parameter name to sweep") + p.add_argument("--values", nargs="+", required=True, help="values to sweep over") + p.add_argument("--outdir", type=Path, default=Path("sweep"), help="output folder") + p.add_argument("--no-view", action="store_true", help="skip launching napari") + args = p.parse_args() + + base = yaml.safe_load(args.config.read_text()) + args.outdir.mkdir(parents=True, exist_ok=True) + wo = _wo_cmd() + # Expand HCS positions here — subprocess does not run a shell to glob. + positions = [str(p) for p in sorted(args.input.glob("*/*/*")) if p.is_dir()] + if not positions: + raise SystemExit(f"no positions found under {args.input}/*/*/*") + + outputs = [] + for raw in args.values: + value = _coerce(raw) + cfg = yaml.safe_load(args.config.read_text()) # fresh copy + if not _set_nested(cfg, args.param, value): + raise SystemExit(f"param {args.param!r} not found in any config block") + + tag = str(raw).replace(".", "p").replace("-", "m") + cfg_path = args.outdir / f"{args.param}_{tag}.yml" + out_path = args.outdir / f"recon_{args.param}_{tag}.zarr" + cfg_path.write_text(yaml.safe_dump(cfg, sort_keys=False)) + + print(f"\n=== {args.param} = {value} → {out_path} ===") + subprocess.run(wo + ["rec", "-i", *positions, "-c", str(cfg_path), "-o", str(out_path)], check=True) + outputs.append(str(out_path)) + + print("\nReconstructions:") + for o in outputs: + print(" ", o) + + if not args.no_view and outputs: + print("\nOpening in napari (grid view)...") + subprocess.run(wo + ["view", *outputs]) + else: + print("\nView with:\n " + " ".join(wo + ["view", *outputs])) + + +if __name__ == "__main__": + main() diff --git a/.skills/waveorder-tutorial/scripts/to_omezarr.py b/.skills/waveorder-tutorial/scripts/to_omezarr.py new file mode 100644 index 00000000..d03d1fb4 --- /dev/null +++ b/.skills/waveorder-tutorial/scripts/to_omezarr.py @@ -0,0 +1,93 @@ +"""Convert a TIFF / .npy / plain-zarr array into a viewable OME-Zarr. + +Writes an HCS OME-Zarr (``.../0/0/0``) with a TCZYX array, named channels, and a +pixel-size scale transform, so it can be read by ``wo rec`` and shown by +``wo view``. + +Examples +-------- + python to_omezarr.py stack.tif ./data.zarr \ + --channel-name Brightfield --axes ZYX \ + --yx-pixel-size 0.1 --z-pixel-size 0.25 + + python to_omezarr.py arr.npy ./data.zarr \ + --channel-name GFP --axes CZYX --yx-pixel-size 0.1 --z-pixel-size 0.25 +""" + +import argparse +from pathlib import Path + +import numpy as np + + +def _load(path: Path) -> np.ndarray: + """Load a TIFF, .npy, or plain zarr array as a numpy array.""" + suffix = path.suffix.lower() + if suffix in {".tif", ".tiff"}: + import tifffile + + return tifffile.imread(str(path)) + if suffix == ".npy": + return np.load(path) + # Assume a plain zarr array directory. + import zarr + + return np.asarray(zarr.open(str(path), mode="r")) + + +def _to_tczyx(arr: np.ndarray, axes: str) -> np.ndarray: + """Reorder/expand an array with the given axis labels into TCZYX.""" + axes = axes.upper() + if len(axes) != arr.ndim: + raise ValueError(f"--axes {axes!r} has {len(axes)} labels but array has {arr.ndim} dims") + if set(axes) - set("TCZYX"): + raise ValueError(f"--axes may only contain T, C, Z, Y, X; got {axes!r}") + + # Move each present axis into TCZYX order. + order = [axes.index(a) for a in "TCZYX" if a in axes] + arr = np.transpose(arr, order) + # Insert singleton dims for any missing axis. + for i, a in enumerate("TCZYX"): + if a not in axes: + arr = np.expand_dims(arr, i) + return arr + + +def main() -> None: + p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + p.add_argument("input", type=Path, help="input .tif/.tiff/.npy or plain zarr dir") + p.add_argument("output", type=Path, help="output .zarr path") + p.add_argument("--axes", default="ZYX", help="axis order of the input (subset of TCZYX)") + p.add_argument("--channel-name", action="append", default=None, help="channel name (repeat per channel)") + p.add_argument("--yx-pixel-size", type=float, default=0.1, help="lateral pixel size (µm)") + p.add_argument("--z-pixel-size", type=float, default=1.0, help="axial pixel size (µm)") + args = p.parse_args() + + from iohub.ngff import open_ome_zarr + from iohub.ngff.models import TransformationMeta + + arr = _to_tczyx(_load(args.input), args.axes) + n_channels = arr.shape[1] + names = args.channel_name or [f"Channel{i}" for i in range(n_channels)] + if len(names) != n_channels: + raise ValueError(f"{len(names)} channel names for {n_channels} channels") + + ds = open_ome_zarr(args.output, layout="hcs", mode="w", channel_names=names) + pos = ds.create_position("0", "0", "0") + pos.create_image( + "0", + arr.astype(np.float32), + transform=[ + TransformationMeta( + type="scale", + scale=[1, 1, args.z_pixel_size, args.yx_pixel_size, args.yx_pixel_size], + ) + ], + ) + ds.close() + print(f"Wrote {args.output} shape TCZYX={arr.shape} channels={names}") + print(f"Inspect with: wo view {args.output}") + + +if __name__ == "__main__": + main() diff --git a/waveorder/cli/view.py b/waveorder/cli/view.py index 825f7fc1..454be6d3 100644 --- a/waveorder/cli/view.py +++ b/waveorder/cli/view.py @@ -1,37 +1,72 @@ import click +def _add_transfer_function_image(viewer, name, arr, shift_axes=None): + """Add one transfer function array to the viewer as real/imag bwr layers.""" + import numpy as np + + # Remove leading singleton dims (T, C) + while arr.ndim > 3 and arr.shape[0] == 1: + arr = arr[0] + + # ifftshift to center the DC component (all axes by default) + arr = np.fft.ifftshift(arr, axes=shift_axes) + + lim = np.max(np.abs(arr)) + if lim == 0: + lim = 1.0 + + viewer.add_image(arr.real, name=f"Re({name})", colormap="bwr", contrast_limits=(-lim, lim)) + if np.iscomplexobj(arr): + viewer.add_image(arr.imag, name=f"Im({name})", colormap="bwr", contrast_limits=(-lim, lim)) + + def _open_transfer_function(viewer, path): - """Open a transfer function zarr, displaying real/imag parts with bwr colormap.""" + """Open a transfer function zarr in napari. + + 2D reconstructions store a singular system (``U``, ``S``, ``Vh``) instead of + a direct transfer function. In that case the transfer function is + reconstructed from the SVD, ``H = U @ diag(S) @ Vh`` at every lateral + frequency, and the resulting transfer function(s) are shown rather than the + raw singular vectors. All arrays are displayed as real/imag parts with a + diverging ``bwr`` colormap. + """ import numpy as np import zarr root = zarr.open(path, mode="r") - for name in root.keys(): - arr = np.array(root[name]) - # Remove leading singleton dims (T, C) - while arr.ndim > 3 and arr.shape[0] == 1: - arr = arr[0] - - # ifftshift to center the DC component - arr = np.fft.ifftshift(arr) - - lim = np.max(np.abs(arr)) - if lim == 0: - lim = 1.0 - viewer.add_image( - arr.real, - name=f"Re({name})", - colormap="bwr", - contrast_limits=(-lim, lim), - ) - if np.iscomplexobj(arr): - viewer.add_image( - arr.imag, - name=f"Im({name})", - colormap="bwr", - contrast_limits=(-lim, lim), - ) + names = list(root.keys()) + + svd_components = {"singular_system_U", "singular_system_S", "singular_system_Vh"} + if svd_components.issubset(names): + # Reconstruct the transfer function from the singular system. + # Stored shapes: U (1, s, k, Vy, Vx), S (1, 1, k, Vy, Vx), + # Vh (1, k, Z, Vy, Vx). H[s, z] = sum_k U[s, k] * S[k] * Vh[k, z]. + U = np.asarray(root["singular_system_U"])[0] # (s, k, Vy, Vx) + S = np.asarray(root["singular_system_S"])[0, 0] # (k, Vy, Vx) + Vh = np.asarray(root["singular_system_Vh"])[0] # (k, Z, Vy, Vx) + H = np.einsum("skyx,kyx,kzyx->szyx", U, S.astype(U.dtype), Vh) # (s, Z, Vy, Vx) + + # Name the output transfer functions by object type. + labels = { + 1: ["transfer_function"], + 2: ["absorption_transfer_function", "phase_transfer_function"], + }.get(H.shape[0], [f"transfer_function_{i}" for i in range(H.shape[0])]) + + for component, label in zip(H, labels): + # Z is a real-space defocus axis, so only ifftshift the lateral + # frequency axes. + _add_transfer_function_image(viewer, label, component, shift_axes=(-2, -1)) + + # Show any remaining (non-SVD) transfer function arrays, if present. + for name in names: + if name.startswith("singular_system_"): + continue + _add_transfer_function_image(viewer, name, np.asarray(root[name])) + return + + for name in names: + _add_transfer_function_image(viewer, name, np.asarray(root[name])) def _is_transfer_function(path):