[HazenBabcock]

Summary

Two defects in how the fitters update overlapping and width-clamped peaks.

2D width clamp (dao_fit.c). The 2D model has one width, but the peak is drawn from XWIDTH in x and YWIDTH in y:

 daoUpdate2D        XWIDTH -= delta;  YWIDTH = XWIDTH
 daoCheck2D         clamp XWIDTH
+                   YWIDTH = XWIDTH        # after the clamp; also mpDaoUpdate()'s only source of YWIDTH
 daoCheck3D
-  daoCheck2D()                            # would now tie YWIDTH to XWIDTH
+  mFitCheck(); clamp XWIDTH               # 3D keeps independent widths
   clamp YWIDTH

Before this change, a clamped 2D peak was drawn elliptical. Multi-plane DAO never wrote YWIDTH at all, so every peak kept its starting y width.

Shared-background Jacobian (19 functions in dao_fit.c, cubic_fit.c, pupil_fit.c, fft_fit.c):

 model:  f = f_data + bg_data / bg_counts     # background averaged over the ROIs covering a pixel
 Jacobian, background column:
-  df/db = rqe
+  df/db = rqe / bg_counts

Where ROIs overlap, the update was aimed at a point that is not the maximum likelihood fit. Its background/height steps also made so little progress that overlapping peaks were marked converged early. Isolated peaks have bg_counts = 1, so their fits are unchanged.

Evidence

Tests. 282 pass before and after. test_mfit_8 / test_mfit_10 round-tripped a pre-specified 2D peak with ysigma ≠ xsigma. They now expect ysigma to come back as the clamped xsigma.

Python replica of mFitIterateLM(), compared with the C library sweep by sweep on 240 simulated fits of 1–4 overlapping peaks:

• Before: 2 runs diverge, exactly where a peak first hits the width clamp.
  After: 240 / 240 match (worst relative difference 1.4e-12) once the replica uses df/db = rqe / bg_counts. The replica with the old derivative now disagrees in exactly the 180 overlapping runs and agrees on the 60 single-peak runs.
• Isolated peaks: bit-identical before and after, 400 / 400.

Simulated clusters (σ 1.5 px, 2000 photons, background 50). Below is the median distance from the joint maximum likelihood fit when doFit() stops, in CRLB units:

model	cluster	1.5σ before → after	2σ before → after	3σ before → after
2dfixed	pair	0.11 → 0.03	0.14 → 0.01	0.11 → 0.003
2dfixed	square of 4	0.42 → 0.08	0.18 → 0.03	0.20 → 0.005
2d	pair	0.12 → 0.08	0.09 → 0.04	0.08 → 0.008
2d	square of 4	1.65 → 0.92	0.64 → 0.21	0.24 → 0.03

End to end (3D-DAOSTORM find + fit on the same clusters), where the right number of peaks was found:

• Error against the truth now matches the joint optimum. xy RMS ratio is 0.96–1.05 (was up to 1.49) and height RMS ratio 0.97–1.03 (was up to 2.12). The exception is lines of three at ≤1.5σ, where the optimum is itself near-degenerate.
• Spurious extra localizations fall, e.g. pairs at 2σ: 28% → 18% (2dfixed), 56% → 42% (2d).

Diagnostics, run on identical simulated data before and after:

diagnostic	before	after
daostorm_3d, sCMOS	—	unchanged (ROIs don't overlap)
spliner	—	summary statistics unchanged
pupilfn, psf_fft	—	XYZ RMSE within ±1.5% (one realisation resolves ~2%); 8–10% fewer iterations
multiplane_dao	recall 0.702, XY RMSE 19.6 / 19.3 nm	recall 0.723, XY RMSE 19.2 / 18.5 nm

daostorm_3d with random placement. This is the "Uniform background, STORM" block in make_data.py: random_list.hdf5 with SimpleSTORM blinking, about 117 emitters on per frame. It is switched off in the repo, so I reproduced it with 500 frames, giving about 58,600 emitter-frames per photon level. Scored by distance to the nearest emitter that is on in the same frame (1000 photons; 500 is similar):

nearest neighbour	emitters	recall before → after	xy RMSE (nm) before → after
< 1σ	854	0.477 → 0.478	36.9 → 39.8
1–2σ	2470	0.483 → 0.485	68.6 → 68.8
2–3σ	3652	0.616 → 0.622	65.6 → 64.4
> 3σ	51572	0.757 → 0.759	28.0 → 28.2

Unmatched localizations fall 412 → 319 (500 photons: 285 → 236). The gains are much smaller here than in the toy clusters. At these separations most of the error is detection and merging (recall below 0.5), and blinking emitters that overlap are usually very unequal in brightness.

Merge Danger

Door: two-way. This is a code-only change with no file format or parameter changes, and reverting restores the old behaviour exactly.

Blast Radius: numerical

• Localizations from every fitter change wherever ROIs overlap. Single-plane isolated peaks do not change.
• Tight clusters take more iterations because peaks no longer stop early: 15–90% more per cluster in the simulations above.
• Multi-plane DAO results change for every peak, since the y width is now fitted. It needs about twice the iterations (mean 6.2 → 12.4 per localization).
• Pre-specified (text / hdf5) 2D peaks now report ysigma = xsigma from the moment they are added, including with no_fitting.
• In the 2D model, fits the data pulls outside the σ range now reach the clamp instead of stalling just short of it. In the random-placement run, 0.3–0.4% of localizations end at 0.75 or 7.5 px (was ~0.07%), while those between 4 and 7.5 px fall by about a third. The mean fitted width is unchanged, but these outliers raise its SD (e.g. 0.31 → 0.36 px for isolated emitters at 1000 photons).

🤖 Generated with Claude Code