From 4fe31c1d42cbeeb0f69e700a6675617274e8523c Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Sun, 18 Oct 2026 14:39:14 -0400 Subject: [PATCH 1/2] Tie YWIDTH to the clamped XWIDTH in the 2D Gaussian model daoUpdate2D() copied XWIDTH into YWIDTH before daoCheck2D() applied the width clamp, and the clamp only touched XWIDTH. Whenever it fired, the peak was drawn with the clamped width in x and the unclamped one in y, while daoCalcJH2D() differentiates a round peak of width XWIDTH, so the error that accepts or rejects the step belonged to a different model. Multi-plane DAO fitting was worse off. mpDaoUpdate() updates XWIDTH and calls daoCheck2D(), and nothing on that path ever wrote YWIDTH, so every peak kept its starting y width for the whole fit. daoCheck2D() now sets YWIDTH from XWIDTH after the clamp, which fixes both. daoCheck3D() used to start by calling daoCheck2D(); it now does its own margin check and clamps the two widths independently, so the 3D model is unchanged. Pre-specified 2D peaks ("text" / "hdf5") pass through the same check when they are added, so a requested ysigma that differs from xsigma is now replaced by xsigma straight away rather than after the first update. test_mfit_8 and test_mfit_10 checked the old round trip; they now check the new one. On identical simulated data the multiplane_dao diagnostic goes from recall 0.702 to 0.725 and XY RMSE 19.6 / 19.3 nm to 19.2 / 18.5 nm, at the cost of twice the fitting iterations (mean 6.2 to 12.4 per localization) now that the y width actually moves. The daostorm_3d, sCMOS, spliner, pupilfn and psf_fft diagnostics are unchanged to within run-to-run noise (~1e-13). Co-Authored-By: Claude Opus 5.5 --- storm_analysis/sa_library/dao_fit.c | 26 +++++++++++++++++++++++--- storm_analysis/test/test_multifit.py | 10 ++++++---- 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/storm_analysis/sa_library/dao_fit.c b/storm_analysis/sa_library/dao_fit.c index 220ca93a..dfab8563 100644 --- a/storm_analysis/sa_library/dao_fit.c +++ b/storm_analysis/sa_library/dao_fit.c @@ -1090,7 +1090,16 @@ int daoCheck2D(fitData *fit_data){ else if(peak->params[XWIDTH] > dao_fit->width_max){ peak->params[XWIDTH] = dao_fit->width_max; } - + + /* + * The 2D model has a single width, but the peak shape is drawn from + * XWIDTH in x and YWIDTH in y. Copy it after the clamp, otherwise a + * clamped peak is drawn with an unclamped YWIDTH. This is also the + * only place that multi-plane DAO fitting, mpDaoUpdate(), gets its + * YWIDTH from. + */ + peak->params[YWIDTH] = peak->params[XWIDTH]; + return 0; } @@ -1110,11 +1119,22 @@ int daoCheck3D(fitData *fit_data){ peak = fit_data->working_peak; dao_fit = (daoFit *)fit_data->fit_model; - - if(daoCheck2D(fit_data)){ + + if(mFitCheck(fit_data)){ return 1; } + /* + * Width range clamp, XWIDTH and YWIDTH independently. Not done by + * calling daoCheck2D() as that ties YWIDTH to XWIDTH. + */ + if(peak->params[XWIDTH] < dao_fit->width_min){ + peak->params[XWIDTH] = dao_fit->width_min; + } + else if(peak->params[XWIDTH] > dao_fit->width_max){ + peak->params[XWIDTH] = dao_fit->width_max; + } + if(peak->params[YWIDTH] < dao_fit->width_min){ peak->params[YWIDTH] = dao_fit->width_min; } diff --git a/storm_analysis/test/test_multifit.py b/storm_analysis/test/test_multifit.py index 47b5f2af..8e8f6856 100644 --- a/storm_analysis/test/test_multifit.py +++ b/storm_analysis/test/test_multifit.py @@ -337,9 +337,10 @@ def test_mfit_8(): "z" : numpy.array([16.0, 26.0])} mfit.newPeaks(peaks, "text") - # Round trip verification. + # Round trip verification. The 2D model has a single width, so ysigma + # comes back as the (clamped) xsigma, not as the ysigma asked for. for pname in peaks: - pvals = peaks[pname] + pvals = peaks["xsigma"] if (pname == "ysigma") else peaks[pname] mvals = mfit.getPeakProperty(pname) for i in range(pvals.size): assert(abs(pvals[i] - mvals[i]) < 1.0e-6) @@ -402,9 +403,10 @@ def test_mfit_10(): "z" : numpy.array([16.0, 26.0])} mfit.newPeaks(peaks, "hdf5") - # Round trip verification. + # Round trip verification. The 2D model has a single width, so ysigma + # comes back as the (clamped) xsigma, not as the ysigma asked for. for pname in peaks: - pvals = peaks[pname] + pvals = peaks["xsigma"] if (pname == "ysigma") else peaks[pname] mvals = mfit.getPeakProperty(pname) for i in range(pvals.size): assert(abs(pvals[i] - mvals[i]) < 1.0e-6) From 9b3e948395fe1bf6a9ff57f3b3393a5b810b0802 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Sun, 18 Oct 2026 14:49:53 -0400 Subject: [PATCH 2/2] Divide the background Jacobian column by bg_counts Overlapping peaks share their background. Each peak adds its background to bg_data and increments bg_counts over its ROI, and the model value at a pixel uses bg_data / bg_counts, the average background of the peaks whose ROIs cover it (see mFitCalcErr()). The Jacobian functions used df/db = rqe at every ROI pixel, which is only right where a single ROI covers the pixel. Where ROIs overlap the derivative is rqe / bg_counts. With the wrong column the per-peak update heads for a point that is not the maximum likelihood estimate of the model. More importantly, the background / height part of each step makes so little progress that the per-peak convergence test marks overlapping peaks converged early. All 19 Jacobian functions that use this model are changed: 3D-DAOSTORM / sCMOS (12, including the two Anscombe variants, whose column is t1 rather than rqei), Spliner (5), pupil function (1) and PSF FFT (1). Multi-plane fitting calls these per plane and picks the change up. Isolated peaks have bg_counts = 1 throughout, so their fits do not change at all (bit-identical on 400 single-peak fits). Measured on simulated clusters of 2 to 4 overlapping emitters (sigma 1.5 px, 2000 photons on a background of 50), against the joint maximum likelihood fit of the same model: - doFit() now stops a median 0.003 to 0.21 CRLB from the joint optimum at 2 to 3 sigma separation, was 0.08 to 0.64; at 1.5 sigma 0.03 to 0.92, was 0.11 to 1.65. - In the full find + fit pipeline, where the right number of peaks was found at 1.5 to 3 sigma, position and height errors against the truth now match the joint optimum: xy RMS ratio 0.96 to 1.05, was up to 1.49, and height RMS ratio 0.97 to 1.03, was up to 2.12. The exception is lines of three at 1.5 sigma and closer, where the joint optimum is itself near-degenerate and the pipeline still does better than it. Spurious extra localizations fall, e.g. pairs 2 sigma apart from 28% to 18% with 2dfixed and 56% to 42% with 2d. - The cost is more iterations in tight clusters, because peaks no longer stop early: 15 to 90% more LM iterations per cluster in those runs. The diagnostics, on identical data: daostorm_3d and sCMOS unchanged (their ROIs do not overlap); spliner summary statistics unchanged; pupilfn and psf_fft move by up to 1.5% in XYZ RMSE in both directions, within what one realisation resolves, and need 8 to 10% fewer iterations. Co-Authored-By: Claude Opus 5.5 --- storm_analysis/psf_fft/fft_fit.c | 3 ++- storm_analysis/pupilfn/pupil_fit.c | 3 ++- storm_analysis/sa_library/dao_fit.c | 36 +++++++++++++++++++---------- storm_analysis/spliner/cubic_fit.c | 15 ++++++++---- 4 files changed, 38 insertions(+), 19 deletions(-) diff --git a/storm_analysis/psf_fft/fft_fit.c b/storm_analysis/psf_fft/fft_fit.c index 779cfe62..e20a13c3 100644 --- a/storm_analysis/psf_fft/fft_fit.c +++ b/storm_analysis/psf_fft/fft_fit.c @@ -98,7 +98,8 @@ void ftFitCalcJH3D(fitData *fit_data, double *jacobian, double *hessian) jt[1] = rqei*height*(dx[l]); jt[2] = rqei*height*(dy[l]); jt[3] = rqei*height*(dz[l]); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate jacobian. */ t1 = 2.0*(1.0 - xi/fi); diff --git a/storm_analysis/pupilfn/pupil_fit.c b/storm_analysis/pupilfn/pupil_fit.c index 1a163e40..b9bf5836 100644 --- a/storm_analysis/pupilfn/pupil_fit.c +++ b/storm_analysis/pupilfn/pupil_fit.c @@ -104,7 +104,8 @@ void pfitCalcJH3D(fitData *fit_data, double *jacobian, double *hessian) jt[1] = rqei*2.0*height*(psf_r[l]*dx_r[l]+psf_c[l]*dx_c[l]); jt[2] = rqei*2.0*height*(psf_r[l]*dy_r[l]+psf_c[l]*dy_c[l]); jt[3] = rqei*2.0*height*(psf_r[l]*dz_r[l]+psf_c[l]*dz_c[l]); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate jacobian. */ t1 = 2.0*(1.0 - xi/fi); diff --git a/storm_analysis/sa_library/dao_fit.c b/storm_analysis/sa_library/dao_fit.c index dfab8563..7eaa2f6b 100644 --- a/storm_analysis/sa_library/dao_fit.c +++ b/storm_analysis/sa_library/dao_fit.c @@ -86,7 +86,8 @@ void daoCalcJH2DFixed(fitData *fit_data, double *jacobian, double *hessian) jt[0] = rqei*e_t; jt[1] = rqei*2.0*a1*width*xt*e_t; jt[2] = rqei*2.0*a1*width*yt*e_t; - jt[3] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[3] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t1 = 2.0*(1.0 - xi/fi); @@ -163,7 +164,8 @@ void daoCalcJH2DFixedALS(fitData *fit_data, double *jacobian, double *hessian) jt[0] = t1*e_t; jt[1] = t1*2.0*a1*width*xt*e_t; jt[2] = t1*2.0*a1*width*yt*e_t; - jt[3] = t1; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[3] = t1/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t2 = (fi - fit_data->as_xi[k]); @@ -238,7 +240,8 @@ void daoCalcJH2DFixedLS(fitData *fit_data, double *jacobian, double *hessian) jt[0] = rqei*e_t; jt[1] = rqei*2.0*a1*width*xt*e_t; jt[2] = rqei*2.0*a1*width*yt*e_t; - jt[3] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[3] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t1 = (fi - fit_data->x_data[k]); @@ -310,7 +313,8 @@ void daoCalcJH2DFixedDWLS(fitData *fit_data, double *jacobian, double *hessian) jt[0] = rqei*e_t; jt[1] = rqei*2.0*a1*width*xt*e_t; jt[2] = rqei*2.0*a1*width*yt*e_t; - jt[3] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[3] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t1 = (fi - xi)/xi; @@ -378,7 +382,8 @@ void daoCalcJH2DFixedFWLS(fitData *fit_data, double *jacobian, double *hessian) jt[0] = rqei*e_t; jt[1] = rqei*2.0*a1*width*xt*e_t; jt[2] = rqei*2.0*a1*width*yt*e_t; - jt[3] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[3] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t1 = (fi - xi)/fi; @@ -448,7 +453,8 @@ void daoCalcJH2D(fitData *fit_data, double *jacobian, double *hessian) jt[1] = rqei*2.0*a1*width*xt*e_t; jt[2] = rqei*2.0*a1*width*yt*e_t; jt[3] = rqei*(-a1*xt*xt*e_t-a1*yt*yt*e_t); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector */ t1 = 2.0*(1.0 - xi/fi); @@ -527,7 +533,8 @@ void daoCalcJH2DALS(fitData *fit_data, double *jacobian, double *hessian) jt[1] = t1*2.0*a1*width*xt*e_t; jt[2] = t1*2.0*a1*width*yt*e_t; jt[3] = t1*(-a1*xt*xt*e_t-a1*yt*yt*e_t); - jt[4] = t1; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = t1/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t2 = (fi - fit_data->as_xi[k]); @@ -632,7 +639,8 @@ void daoCalcJH2DLS(fitData *fit_data, double *jacobian, double *hessian) jt[1] = rqei*2.0*a1*width*xt*e_t; jt[2] = rqei*2.0*a1*width*yt*e_t; jt[3] = rqei*(-a1*xt*xt*e_t-a1*yt*yt*e_t); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t1 = (fi - fit_data->x_data[k]); @@ -706,7 +714,8 @@ void daoCalcJH2DDWLS(fitData *fit_data, double *jacobian, double *hessian) jt[1] = rqei*2.0*a1*width*xt*e_t; jt[2] = rqei*2.0*a1*width*yt*e_t; jt[3] = rqei*(-a1*xt*xt*e_t-a1*yt*yt*e_t); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t1 = (fi - xi)/xi; @@ -780,7 +789,8 @@ void daoCalcJH2DFWLS(fitData *fit_data, double *jacobian, double *hessian) jt[1] = rqei*2.0*a1*width*xt*e_t; jt[2] = rqei*2.0*a1*width*yt*e_t; jt[3] = rqei*(-a1*xt*xt*e_t-a1*yt*yt*e_t); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t1 = (fi - fit_data->x_data[k])/fi; @@ -852,7 +862,8 @@ void daoCalcJH3D(fitData *fit_data, double *jacobian, double *hessian) jt[2] = rqei*(-a1*xt*xt*e_t); jt[3] = rqei*2.0*a1*a5*yt*e_t; jt[4] = rqei*(-a1*yt*yt*e_t); - jt[5] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[5] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t1 = 2.0*(1.0 - xi/fi); @@ -946,7 +957,8 @@ void daoCalcJHZ(fitData *fit_data, double *jacobian, double *hessian) jt[1] = rqei*2.0*a1*a3*xt*e_t; jt[2] = rqei*2.0*a1*a5*yt*e_t; jt[3] = rqei*(-a1*xt*xt*gx*e_t-a1*yt*yt*gy*e_t); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate b vector. */ t1 = 2.0*(1.0 - xi/fi); diff --git a/storm_analysis/spliner/cubic_fit.c b/storm_analysis/spliner/cubic_fit.c index 096a76e0..eee137de 100644 --- a/storm_analysis/spliner/cubic_fit.c +++ b/storm_analysis/spliner/cubic_fit.c @@ -89,7 +89,8 @@ void cfCalcJH2D(fitData *fit_data, double *jacobian, double *hessian) jt[0] = rqei*peak->psf[j]; jt[1] = -rqei*height*dxfAt2D(spline_fit->spline_data,l+y_start,m+x_start); jt[2] = -rqei*height*dyfAt2D(spline_fit->spline_data,l+y_start,m+x_start); - jt[3] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[3] = rqei/((double)fit_data->bg_counts[k]); /* Calculate jacobian. */ t1 = 2.0*(1.0 - xi/fi); @@ -163,7 +164,8 @@ void cfCalcJH3D(fitData *fit_data, double *jacobian, double *hessian) jt[1] = -rqei*height*dxfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); jt[2] = -rqei*height*dyfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); jt[3] = rqei*height*dzfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate jacobian. */ t1 = 2.0*(1.0 - xi/fi); @@ -245,7 +247,8 @@ void cfCalcJH3DALS(fitData *fit_data, double *jacobian, double *hessian) jt[1] = -t1*height*dxfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); jt[2] = -t1*height*dyfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); jt[3] = t1*height*dzfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); - jt[4] = t1; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = t1/((double)fit_data->bg_counts[k]); /* Calculate jacobian. */ t2 = (fi - fit_data->as_xi[k]); @@ -325,7 +328,8 @@ void cfCalcJH3DLS(fitData *fit_data, double *jacobian, double *hessian) jt[1] = -rqei*height*dxfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); jt[2] = -rqei*height*dyfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); jt[3] = rqei*height*dzfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate jacobian. */ t1 = (fi - fit_data->x_data[k]); @@ -406,7 +410,8 @@ void cfCalcJH3DFWLS(fitData *fit_data, double *jacobian, double *hessian) jt[1] = -rqei*height*dxfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); jt[2] = -rqei*height*dyfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); jt[3] = rqei*height*dzfAt3D(spline_fit->spline_data,zi,l+y_start,m+x_start); - jt[4] = rqei; + /* The background is shared with any overlapping peaks, see mFitCalcErr(). */ + jt[4] = rqei/((double)fit_data->bg_counts[k]); /* Calculate jacobian. */ t1 = (fi - xi)/fi;