Plan.md section 7 requires V1-V3 to pass before any phase/abnormality work begins. This document records what each check does, the closed-form derivations behind it, and one real modelling bug the process caught.
Plan.md references two prior scripts, oblique_section_model.py (sphere
pair, closed form, said to be validated to 0.06% against voxel integration)
and shape_sensitivity.py (oblate spheroid variant), with specific target
figures (median area ratio 0.87, 22% below 0.6, 26% "scoreable" at
r=2.58um, d=10um, t=3um for the sphere-pair case). Neither script is present
in this environment, and the exact sampling convention behind those figures
(what "scoreable" means, the offset-sampling range, whether ratio is
measured on the same axis-projection convention used here) cannot be
recovered from the plan text alone.
Rather than guess at a convention that reproduces three specific numbers, M1 re-derives the sphere-pair and oblate-spheroid slab-intersection geometry from first principles (below) and validates it the way the plan says the original was validated: closed form checked against an independent numeric measurement. If the original scripts resurface, V1 should be extended to diff against them directly.
Both primitives in mitosim/geometry/shapes.py expose slab_projected_area(s, t):
the area of the xy-projection of shape ∩ slab, where the slab is
z in [s - t/2, s + t/2] (docs/conventions.md).
For a sphere of radius r centred at world-z zc, the world-z cross-section
at height z is a disk of radius sqrt(r^2 - (z - zc)^2), always centred at
the same (x, y). The union of these disks over an interval is therefore
just the single widest disk in that interval, at the closest-to-equator
point z0 = clip(0, zlo, zhi) (local coordinates, zlo/zhi = slab bounds
clipped to [-r, r]):
area = pi * max(r^2 - z0^2, 0)
An oblate spheroid (a, a, c) with a > c, symmetry axis axis (unit
vector), is {x : x^T Q x <= 1} in world coordinates, Q = R M^-1 R^T with
M = diag(a^2, a^2, c^2) and R any orthonormal frame whose third column
is axis (the result doesn't depend on the other two columns, since the
solid is rotationally symmetric about axis).
Writing Q in xy/z block form (Qxy, q = Q[:2,2], qzz = Q[2,2]), the
cross-section at world-z = z is an ellipse in (x, y), obtained by
completing the square:
u = (x, y) + z * v, v = Qxy^-1 q
u^T Qxy u <= 1 - kappa * z^2, kappa = qzz - q . Qxy^-1 . q
First-pass model (wrong). The natural first guess is that, like the
sphere, the footprint over an interval is just the single widest slice: peak
half-extent rho = 1/sqrt(kappa), peak coefficient A0 = a^2 c / rho, area
= pi * A0 * (1 - z0^2/rho^2). This is what M1 initially shipped.
Why it's wrong. The ellipse is centred at (x, y) = -z * v, which
shifts with z whenever v != 0 (i.e. whenever the tilt couples the xy
and z directions -- true for essentially any tilt not aligned with a
principal axis). The single-widest-slice formula silently assumes the
cross-sections are concentric, which is only true for the sphere (or for
axis exactly +/-z_hat). A brute-force voxel cross-check
(tests/test_v1_analytic_regression.py) caught this directly: for tilted
configurations the true (voxel-measured) footprint came out 5-25% larger
than the single-slice formula, because off-centre slices poke out sideways
beyond where the widest centred slice would reach.
Corrected model. The footprint is the union, over z in the slab, of a
family of similar ellipses whose centre drifts linearly in z. Rotate Qxy
to its principal axes and rescale (w = Qxy^{1/2} u) so the family becomes
circles of radius sqrt(1 - z^2/rho^2) centred at mu * z along a fixed
direction in w-space (mu = |Qxy^{1/2} v|). For fixed w-space Y, the
covered X-interval as z ranges over [zlo, zhi] is
[min_z(mu*z - sqrt(g)), max_z(mu*z + sqrt(g))] with g(z, Y) = 1 - z^2/rho^2 - Y^2 (the interval is connected because g is concave in z).
The maximiser/minimiser has a closed form (stationarity of mu*z +/- sqrt(g)):
z_crit = mu * rho^2 * sqrt(1 - Y^2) / sqrt(1 + mu^2 * rho^2)
clipped to [zlo, zhi] (compare against the endpoints too, since the
critical point may fall outside the truncation window). Integrating the
resulting width over Y in [-1, 1] (via scipy.integrate.quad) and
dividing by sqrt(det(Qxy)) gives the true footprint area
(OblateSpheroid.slab_projected_area). This reduces exactly to the sphere
formula when a == c (v = 0 identically) --
test_oblate_spheroid_formula_reduces_to_sphere_when_a_equals_c checks
this directly.
Why this matters beyond M1. The same trap applies to any tilted,
non-axisymmetric-about-z geometry -- in particular the swept-tube
chromosome/chromatid model in Plan.md section 3.2. Anyone extending
slab_projected_area-style logic to chromatid geometry should re-derive
from the actual cross-section family rather than assuming the widest slice
suffices; validate against a brute-force numeric measurement before trusting
a closed form for a tilted shape.
tests/test_v1_analytic_regression.py draws random configurations (centre
offset, tilt, slab offset) for both primitives and checks the closed-form
slab_projected_area against a dedicated high-precision numeric reference:
an xy grid at 0.02um spacing, with z sampled on 401 points anchored exactly
at the true slab boundaries s +/- t/2 (not snapped to a pre-built voxel
grid's z-planes -- see the V2/production-grid caveat below). Agreement is
required to within 0.5%.
tests/test_v2_conservation_convergence.py checks that render/voxelize.py's
production voxel grid (a single isotropic spacing, built independently of
any particular slab) reproduces each shape's analytic volume to <1% at
0.05um spacing, and that the error shrinks monotonically as spacing is
refined (0.4 -> 0.2 -> 0.1um).
Caveat recorded here on purpose: this general-purpose grid is not used
for V1's area cross-check, because a grid built independently of (s, t)
snaps the slab boundary to the nearest voxel z-plane. Right at a boundary
where the cross-section is thin and changing fast with z (e.g. near a
sphere's pole), even a sub-voxel snap can produce a large relative area
error, though the absolute error stays of order spacing. This is an
expected property of any fixed-resolution production grid, not a bug --
M4's optical convergence study should keep it in mind when reasoning about
z-resolution near mass boundaries.
tests/test_v3_orientation_sampling.py checks: sampled orientation vectors
are unit length; sin(tilt) matches Uniform[0, 1] (Kolmogorov-Smirnov,
p > 0.01); and the fraction of spindles tilted more than 30 degrees is
0.5 +/- 0.01 (Plan.md section 6.1, section 11.1).
mitosim/render/optics.py implements Plan.md section 6.2's pipeline:
voxelise chromatin (and the cell-body sphere) on a fine isotropic grid,
integrate a defocus-dependent Gaussian blur per z-slice into a path-length
map, apply it through Beer-Lambert to get an optical-density-like alpha,
composite hematoxylin/eosin colour over a background field of placeholder
resting nuclei, add scanner noise and colour jitter, then resample to the
target pixel size and crop/pad to a fixed patch. render_patch(..., include_color=False) stops after the Beer-Lambert step and returns the
raw chromatin optical-density map at native (fine) resolution -- steps 5-7
sit behind that flag per the plan, since geometric experiments don't need
colour and are meaningfully faster without it.
Every constant in OpticsParams (extinction coefficients, background
nucleus density, noise/jitter scale) is a placeholder chosen to look
plausible, not a measured value -- exactly the situation V4 exists to fix.
The one constant derived from physics rather than guessed is the default
PSF sigma, from the Abbe diffraction limit for a 0.75 NA objective
(abbe_resolution_um -> FWHM -> Gaussian sigma).
Tests in tests/test_m4_optics.py are engineering checks, not a
calibration: Beer-Lambert monotonicity and bounds, energy conservation and
correct widening of the defocus blur away from the focal plane, patch
shape/range/reproducibility, and that an empty CellState renders as
literally no chromatin. None of this substitutes for V4.
V1-V3 pass (pytest tests/). M4's optical model is implemented and unit
tested (test_m4_optics.py), but V4 (morphometric calibration against real
patches) and V5 (pathologist face-validity) remain out of scope --
both need data and pathologist time this environment doesn't have access
to (Plan.md section 9). Until V4 lands, no figure produced by this
pipeline should be read as quantitatively realistic, only as structurally
plausible -- see the M4 demo's own disclaimer.