From a518b134a3b39421cc272dd7780d04a5378b16bd Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 13:29:15 -0400 Subject: [PATCH 01/13] Make the Spliner README match the directory Cross-checking every entry against the files on disk turned up four problems: find_peaks_fista.py listed, renamed to find_peaks_decon.py by 2a3a8faf, "we are no longer limited to just FISTA". The description went stale with the name: the module now has both a FISTAPeakFinder and an ADMMPeakFinder. offset_to_Z.py listed, the file is offset_to_z.py. spline_info.py on disk since 01/18, never listed. It prints a spline's details, which is not the same thing as the listed print_psf.py, a debugging tool that dumps PSF values. measure_psf_beads.py described as taking "a text for input rather than a molecule list file (.bin file)". The contrast is still right, it does read bead positions from a text file, but measure_psf.py has read HDF5 through saH5Py for years. With these fixed the README and the directory agree exactly, 21 entries each way, nothing listed that is missing and nothing present that is unlisted. No test run: this is a text file that nothing in the test suite or the packaging reads. Co-Authored-By: Claude Opus 5 --- storm_analysis/spliner/README.txt | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/storm_analysis/spliner/README.txt b/storm_analysis/spliner/README.txt index d9e6ee2f..c240df5f 100644 --- a/storm_analysis/spliner/README.txt +++ b/storm_analysis/spliner/README.txt @@ -12,8 +12,8 @@ cubic_spline_c.py - The Python wrapper of the C spline library. Given the splines coefficients this can be used to calculate the spline and it's derivatives. -find_peaks_fista.py - Peak finding using a compressed sensing approach - and the FISTA solver. +find_peaks_decon.py - Peak finding using a compressed sensing + deconvolution approach, with either the FISTA or the ADMM solver. find_peaks_std.py - Peak finding using a similar approach to 3D-DAOSTORM, start with the brightest local maxima and the proceed to dimmer @@ -25,12 +25,13 @@ hdf5_to_beads.py - Convert an HDF5 file to a beads file for PSF measure_psf.py - Measure the PSF of the microscope given a z stack movie and the location of the emitters in each frame. -measure_psf_beads.py - This is similar to measure_psf except that it uses - a text for input rather than a molecule list file (.bin file). +measure_psf_beads.py - This is similar to measure_psf except that it + takes the bead positions from a text file rather than from an HDF5 + localization file. measure_psf_utils.py - Utility functions used for PSF measurement. -offset_to_Z.py - Converts a .off file (from storm-control) to a z_offset +offset_to_z.py - Converts a .off file (from storm-control) to a z_offset file that can be used for measuring the PSF. print_psf.py - Print out some information about a PSF file. @@ -45,6 +46,8 @@ spline3D.py - 3D cubic spline in Python. spline_analysis.py - Run spline analysis on a single STORM movie. +spline_info.py - Print out the details of a spline. + spline_to_psf.py - Generate a PSF at a particular z value given a spline. From 57cda57689823ddb9edeed8383166d820cc019b7 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 13:30:33 -0400 Subject: [PATCH 02/13] Stop measure_psf.py calling its input a .bin file Two Insight3 era references left in a script that has read HDF5 for years: parser = argparse.ArgumentParser(description = 'Measure PSF given a movie, a list.bin file and (optionally) a z_offset file') # Load dax file, z offset file and molecule list file. measurePSF() opens its localizations with saH5Py.SAH5Py(), and the --bin argument's own help text already says "The name of the localizations HDF5 file". So the usage line a user sees first contradicts the help text three lines below it, and names a file format the project has left. Renaming the --bin argument itself would be a user visible change to the command line, so that is left alone. This is the wording only. Co-Authored-By: Claude Opus 5 --- storm_analysis/spliner/measure_psf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/storm_analysis/spliner/measure_psf.py b/storm_analysis/spliner/measure_psf.py index d5b7a2d6..cd6fa29e 100644 --- a/storm_analysis/spliner/measure_psf.py +++ b/storm_analysis/spliner/measure_psf.py @@ -54,7 +54,7 @@ def measurePSF(movie_name, zfile_name, movie_h5_name, psf_name, want2d = False, # Create z scaling object. z_sclr = measurePSFUtils.ZScaler(z_range, z_step) - # Load dax file, z offset file and molecule list file. + # Load the movie, the z offset file and the localizations. dax_data = datareader.inferReader(movie_name) z_off = None if os.path.exists(zfile_name): @@ -207,7 +207,7 @@ def measurePSF(movie_name, zfile_name, movie_h5_name, psf_name, want2d = False, import argparse - parser = argparse.ArgumentParser(description = 'Measure PSF given a movie, a list.bin file and (optionally) a z_offset file') + parser = argparse.ArgumentParser(description = 'Measure PSF given a movie, a localizations HDF5 file and (optionally) a z_offset file') parser.add_argument('--movie', dest='movie', type=str, required=True, help = "The name of the movie to analyze, can be .dax, .tiff or .spe format.") From 67fdf93411b3b3bd1034ac51c4ffb2ba8572aa14 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 13:31:39 -0400 Subject: [PATCH 03/13] Stop psf_localizations.py calling its files molecule lists The same Insight3 era wording as 57cda576, in the multiplane script that does the equivalent job. The module doc string described its input as "a molecule list" and its output as "molecule lists", and the --frame help text referred to "the frame in .bin file". psfLocalizations() reads with saH5Py.SAH5Py() and writes with saH5Py.saveLocalizations(), producing files named _c1_psf.hdf5 and so on, one per channel. Nothing about it has involved a .bin file for years. Also fixes "Giving a mapping file" to "Given", and says what the outputs actually are, which is one localization file per channel rather than an unspecified number of lists. The --bin argument name is left alone, as renaming it would change the command line. Co-Authored-By: Claude Opus 5 --- storm_analysis/multi_plane/psf_localizations.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/storm_analysis/multi_plane/psf_localizations.py b/storm_analysis/multi_plane/psf_localizations.py index bbc9fe0b..bb3aeb4c 100644 --- a/storm_analysis/multi_plane/psf_localizations.py +++ b/storm_analysis/multi_plane/psf_localizations.py @@ -1,8 +1,8 @@ #!/usr/bin/env python """ -Giving a mapping file (from multi_plane.mapper), and a -molecule list, generate molecule lists to use for the -PSF extraction step. +Given a mapping file (from multi_plane.mapper) and a localizations +HDF5 file, generate the per channel localization files to use for +the PSF extraction step. Hazen 05/17 """ @@ -127,7 +127,7 @@ def psfLocalizations(h5_filename, mapping_filename, frame = 0, aoi_size = 8, min parser.add_argument('--map', dest='mapping', type=str, required=True, help = "The name of the mapping file. This is the output of multi_plane.mapper.") parser.add_argument('--frame', dest='frame', type=int, required=False, default=0, - help = "The frame in .bin file to get the localizations from. The default is 0.") + help = "The frame of the localizations file to get the localizations from. The default is 0.") parser.add_argument('--aoi_size', dest='aoi_size', type=int, required=False, default=8, help = "The size of the area of interest around the bead in pixels. The default is 8.") parser.add_argument('--min_height', dest='min_height', type=float, required=False, default = 0.0, From b38b1f4c7fbdc9b1c4a79f095635d8ddb72c4533 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 13:35:03 -0400 Subject: [PATCH 04/13] Stop mFitAnscombe documenting a parameter it does not take The comment block listed two parameters: * x - Signal / function value. * var - (Gaussian) variance. */ double mFitAnscombe(double x) There is no var. The function has taken a single argument since it was introduced in d0ab0dd8 (2018-05-04), so this was never a leftover from a signature change, it was wrong when it was written. A generalized Anscombe transform does take a variance, which is presumably where it came from, but this one is the plain 2*sqrt(x + 3/8) and the code matches that. Comment only, and provably so: multi_fit.o is byte identical before and after, md5 4b8058df50f0d1b692237950820df4aa. Checked by scanning every comment block in every .c file in the tree against the parameter list of the function below it. Excluding the "Return" lines, which are return value documentation rather than parameters, mFitAnscombe was one of four mismatches and is now absent from the scan. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/multi_fit.c | 1 - 1 file changed, 1 deletion(-) diff --git a/storm_analysis/sa_library/multi_fit.c b/storm_analysis/sa_library/multi_fit.c index ca81f46c..3e617501 100644 --- a/storm_analysis/sa_library/multi_fit.c +++ b/storm_analysis/sa_library/multi_fit.c @@ -62,7 +62,6 @@ void mFitAddPeak(fitData *fit_data) * Calculate the Anscombe transform. * * x - Signal / function value. - * var - (Gaussian) variance. */ double mFitAnscombe(double x) { From 54e252627302396188fdf6a0dcc8f08f376bfa47 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 13:38:09 -0400 Subject: [PATCH 05/13] Name the right parameters in three more comment blocks The scan that found the mFitAnscombe defect in b38b1f4c turned up three others, all the same shape but from renames rather than a parameter that never existed: decon_storm/mlem_sparse.c backwardCompressed() and backwardCompressedFixedBg() both document "cmpr", but take "compression". cmpr is a local inside each, 1.0/(1.0 + compression), so the comment describes the derived value rather than the argument. sa_library/matched_filter.c cleanup() documents "filter", but takes "flt". Worse than a plain mismatch, since filter is also the name of the struct type in its own signature, void cleanup(filter *flt). Corrects the names only. Both mlem_sparse functions also leave their first parameter, dec, undocumented, but writing documentation that was never there is a different job from fixing documentation that is wrong. Comment only. Verified by compiling each file at the parent commit and at this one and comparing the objects, which are byte identical: matched_filter 0242ef0177f1e28e78e3288fd3ef06c9, mlem_sparse 65c6e99e4335d80848df5e5b9a6830f9. Both have to be compiled from the same filename for this to mean anything, since gcc embeds the source name. Neither has a standing object file to compare against, which is worth knowing. matched_filter.c is compiled straight into libmatched_filter.so with no intermediate .o kept, and decon_storm is not in SConstruct at all, so mlem_sparse.c is never built as part of the package. That directory is largely MATLAB. Note that matched_filter.c has CRLF line endings, unlike almost everything else here, so it was edited in binary. Its CRLF count is unchanged at 239. Co-Authored-By: Claude Opus 5 --- storm_analysis/decon_storm/mlem_sparse.c | 4 ++-- storm_analysis/sa_library/matched_filter.c | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/storm_analysis/decon_storm/mlem_sparse.c b/storm_analysis/decon_storm/mlem_sparse.c index 46398bb2..d580b809 100644 --- a/storm_analysis/decon_storm/mlem_sparse.c +++ b/storm_analysis/decon_storm/mlem_sparse.c @@ -509,7 +509,7 @@ void backward(decon *dec) * Calculate ej and xj w/ compression parameter * following Lingenfelter L1 norm. * - * cmpr - compression parameter + * compression - compression parameter */ void backwardCompressed(decon *dec, double compression) { @@ -539,7 +539,7 @@ void backwardCompressed(decon *dec, double compression) * following Lingenfelter L1 norm, but do not * update the background. * - * cmpr - compression parameter + * compression - compression parameter */ void backwardCompressedFixedBg(decon *dec, double compression) { diff --git a/storm_analysis/sa_library/matched_filter.c b/storm_analysis/sa_library/matched_filter.c index f00c335b..d4c827fc 100644 --- a/storm_analysis/sa_library/matched_filter.c +++ b/storm_analysis/sa_library/matched_filter.c @@ -61,7 +61,7 @@ filter *initialize(double *, double, int, int, int); /* * cleanup() * - * filter - A pointer to a filter structure. + * flt - A pointer to a filter structure. */ void cleanup(filter *flt) { From b19b1b040ed1cdbf731ee99d12eaba8bbf7593f6 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 13:43:01 -0400 Subject: [PATCH 06/13] Describe mFitCalcErrFWLS as the fit weighted version Its comment read "The data weighted least squares version of the error function", copied verbatim from mFitCalcErrDWLS immediately above it. The distinction those words name is the entire reason both functions exist: DWLS err_i[k] = di*di/fit_data->x_data[k]; the measurement FWLS err_i[k] = di*di/fi; the fit Also says what follows from that. FWLS divides by the fit, so it checks for a negative fit and returns 1 when it finds one, incrementing n_neg_fi. DWLS divides by the measurement and needs no such guard. A reader given the DWLS description has no way to understand why only one of the two can fail. The other three variants, ALS, LS and DWLS, all describe themselves correctly, so FWLS was the only one. Comment only: multi_fit.o is byte identical before and after, md5 4b8058df50f0d1b692237950820df4aa. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/multi_fit.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/storm_analysis/sa_library/multi_fit.c b/storm_analysis/sa_library/multi_fit.c index 3e617501..d1c5af46 100644 --- a/storm_analysis/sa_library/multi_fit.c +++ b/storm_analysis/sa_library/multi_fit.c @@ -386,7 +386,9 @@ int mFitCalcErrDWLS(fitData *fit_data) /* * mFitCalcErrFWLS() * - * The data weighted least squares version of the error function. + * The fit weighted least squares version of the error function. This + * weights by the fit value where mFitCalcErrDWLS() weights by the + * measurement, so unlike DWLS it fails if the fit goes negative. * * fit_data - pointer to a fitData structure. * From 77eef16fa03b4efcb647851d42e82670b8d3c58c Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 14:28:42 -0400 Subject: [PATCH 07/13] Fix the ZCENTER delta in the convergence comment Every threshold in mFitDeltaConvergence() is annotated with the delta it represents at the default tolerance of 1e-6. Six of the seven agree with their multiplier, and they agree exactly: HEIGHT 10000 * tol 0.01 comment 0.01 BACKGROUND 10000 * tol 0.01 comment 0.01 XCENTER 1000 * tol 0.001 comment 0.001 YCENTER 1000 * tol 0.001 comment 0.001 ZCENTER 100 * tol 0.0001 comment 0.01 <- XWIDTH 1000 * tol 0.001 comment 0.001 YWIDTH 1000 * tol 0.001 comment 0.001 So 100 means 0.0001 and the comment is the thing that is wrong, not the multiplier. The units say the same: x and y converge at 0.001 pixels, which is 0.1nm at a 100nm pixel, and z at 0.0001 microns is also 0.1nm. Believing the comment would make the z criterion a hundred times looser than the lateral ones for no stated reason. Changing the multiplier instead would be a change to when fitting stops, so it would need measuring rather than reasoning, and does not belong in a documentation branch. Note that none of this runs today in any case: DELTA_CONVERGENCE defaults to 0, so the whole function is unreachable. Comment only: multi_fit.o is byte identical before and after, md5 4b8058df50f0d1b692237950820df4aa. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/multi_fit.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/storm_analysis/sa_library/multi_fit.c b/storm_analysis/sa_library/multi_fit.c index d1c5af46..61faaa19 100644 --- a/storm_analysis/sa_library/multi_fit.c +++ b/storm_analysis/sa_library/multi_fit.c @@ -612,7 +612,7 @@ int mFitDeltaConvergence(fitData *fit_data, int index) return 0; } - /* 0.01 delta */ + /* 0.0001 delta */ if(fabs(old_params[ZCENTER] - cur_params[ZCENTER]) > (100.0 * fit_data->tolerance)){ return 0; } From 275812c88233dc26f065821b6ebdcb5053838e5a Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 14:30:47 -0400 Subject: [PATCH 08/13] Say which of z_range and z_step divides the other The comment above ZScaler's assertion had the relationship backwards: # Assert that the z_step size is a multiple of the z_range. assert ((int(z_range*1.0e+3) % int(z_step*1.0e+3)) == 0), \ "The z range must be a multiple of the z step." z_range % z_step tests whether the range is a multiple of the step, which is what the assertion message says, and the comment says the reverse. Exercised directly: range=0.6 step=0.2 accepted 600 % 200 = 0 range=0.5 step=0.3 rejected 500 % 300 = 200 range=0.6 step=0.4 rejected 600 % 400 = 200 so a step larger than the range cannot be what is meant, and in fact the assertion two lines above already rejects that case outright. The code and the failure message were both right, only the comment was wrong, which is the awkward case: someone debugging a failure reads the comment first and it tells them to change the wrong number. Co-Authored-By: Claude Opus 5 --- storm_analysis/spliner/measure_psf_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/storm_analysis/spliner/measure_psf_utils.py b/storm_analysis/spliner/measure_psf_utils.py index 91f87c0a..95e02fe5 100644 --- a/storm_analysis/spliner/measure_psf_utils.py +++ b/storm_analysis/spliner/measure_psf_utils.py @@ -25,7 +25,7 @@ def __init__(self, z_range, z_step): assert(z_step > 0.0), "The z step must be positive." assert(z_range >= z_step), "The z range must be greater than or equal to the step size." - # Assert that the z_step size is a multiple of the z_range. + # Assert that the z_range is a multiple of the z_step size. assert ((int(z_range*1.0e+3) % int(z_step*1.0e+3)) == 0), "The z range must be a multiple of the z step." self.z_mid = int(round(z_range/z_step)) From 47e67fb17e22e06c5a23561b951fa021e2586c7a Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 14:32:57 -0400 Subject: [PATCH 09/13] State the condition under which the spline library is thread safe The header said: * 3. This library is thread safe.. Pretty sure.. Which is the least useful form the claim could take. It is unqualified, so it reads as a guarantee, and it is hedged, so it invites someone to test it. The truth is conditional and the condition matters. computeDelta2D() and computeDelta3D() write their results into the delta_f, delta_dxf, delta_dyf and delta_dzf arrays, which live in the splineData struct rather than on the stack. The accessors then read them back. So the evaluation of a point is a write followed by a read through shared state, and two threads evaluating the same splineData will overwrite each other's deltas in between. Per splineData it is safe, and initSpline2D() and initSpline3D() copy the coefficients into their own allocation rather than aliasing the caller's, so two splineData share nothing at all. Nothing is broken today, storm-analysis scales with multiprocessing rather than threads. Also fixes note 1 in the same block, "but this whether or not this works correctly has not been tested", which lost a word somewhere. Comment only: cubic_spline.o is byte identical before and after, md5 b0b37b59a7cd4270e7014183ed032699. Co-Authored-By: Claude Opus 5 --- storm_analysis/spliner/cubic_spline.c | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/storm_analysis/spliner/cubic_spline.c b/storm_analysis/spliner/cubic_spline.c index 60962d52..0b5ac484 100644 --- a/storm_analysis/spliner/cubic_spline.c +++ b/storm_analysis/spliner/cubic_spline.c @@ -3,13 +3,21 @@ * * Notes: * - * 1. Nominally xsize, ysize and zsize could be different, but this whether - * or not this works correctly has not been tested. + * 1. Nominally xsize, ysize and zsize could be different, but whether or + * not this works correctly has not been tested. * * 2. Unlike the Python version, this will not return the correct value if * x is exactly at the maximum value that the spline covers. * - * 3. This library is thread safe.. Pretty sure.. + * 3. This library is thread safe provided that each thread has its own + * splineData. Sharing one between threads is not safe: computeDelta2D() + * and computeDelta3D() write their results into the delta_f, delta_dxf + * and delta_dyf / delta_dzf arrays that live in the struct, and the + * f/dxf/dyf/dzf At2D() and At3D() accessors read them back, so two + * threads evaluating the same splineData will overwrite each other's + * deltas between the write and the read. initSpline2D() and + * initSpline3D() copy the coefficients into the struct, so separate + * splineData share nothing and need no locking. * * Hazen 11/16 * From 5d0f3db5dfa99da778ead3ee101b7ec951d86784 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 14:35:36 -0400 Subject: [PATCH 10/13] Say that the Z model uses five of the seven width terms getWidthParams() returns seven values, wo, c, d, A, B, C and D, and z_calibration.doFit() will fit all four of the higher order terms (assert(n_additional < 5)). The consumers then disagree: sa_utilities/fitz.c 7 terms, uses params[5] and params[6] sa_utilities/fitz_c.py 5 terms, calcSxSy() stops at params[4] sa_library/dao_fit.h 5 terms, double wx_z_params[5] daoInitializeZ() copies with for(i=0;i<5;i++), so a calibration with C and D terms is applied in full when z is calculated after fitting for the "3d" model, and silently truncated by the Z fitting model and by calcSxSy(). Nothing warns. Either could be the right design, the Z model may well not want the extra freedom, so this does not change any of it. What it changes is that you no longer have to read three files to find out. The note goes in all three places a reader can land: - parameters.py, on the shared doc string for all seven parameters, which is where they are actually documented. - dao_fit.h, on the [5] arrays, which is where the truncation happens. - fitz_c.py, on calcSxSy(), the other five term consumer. Worth recording that the default is safe: z_calibration takes --fit_order 2, which produces exactly five terms, so this only bites someone who deliberately fits a higher order curve. Comment only: dao_fit.o is byte identical before and after, md5 b2dc5b13076ae6d61574c86d6345c5ad. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/dao_fit.h | 13 +++++++++++-- storm_analysis/sa_library/parameters.py | 10 ++++++++++ storm_analysis/sa_utilities/fitz_c.py | 4 ++++ 3 files changed, 25 insertions(+), 2 deletions(-) diff --git a/storm_analysis/sa_library/dao_fit.h b/storm_analysis/sa_library/dao_fit.h index 04dd3340..52181c9f 100644 --- a/storm_analysis/sa_library/dao_fit.h +++ b/storm_analysis/sa_library/dao_fit.h @@ -32,8 +32,17 @@ typedef struct daoFit double width_min; /* Minimum allowed width for Gaussian fitting. */ double width_max; /* Maximum allowed width for Gaussian fitting. */ - double wx_z_params[5]; /* x width versus z parameters. */ - double wy_z_params[5]; /* y width versus z parameters. */ + /* + * x and y width versus z parameters, wo, c, d, A and B. + * + * Note that the calibration can have seven terms, and that + * sa_utilities/fitz.c uses all seven. The Z fitting model uses only + * these five, so daoInitializeZ() copies five and any C and D terms in + * the calibration are dropped without warning. See the wx/wy vs z + * parameter documentation in sa_library/parameters.py. + */ + double wx_z_params[5]; + double wy_z_params[5]; } daoFit; diff --git a/storm_analysis/sa_library/parameters.py b/storm_analysis/sa_library/parameters.py index 1007f423..4481693e 100644 --- a/storm_analysis/sa_library/parameters.py +++ b/storm_analysis/sa_library/parameters.py @@ -492,6 +492,16 @@ def __init__(self, **kwds): ww_doc_string = """wx/wy vs z parameters. Units are nanometers or dimensionless. See Huang, Science 2008 for a more detailed explanation. + +Note that not everything uses all seven terms. sa_utilities/fitz.c, which +calculates z from wx and wy after fitting for the "3d" model, uses all of +them. The "Z" fitting model and sa_utilities.fitz_c.calcSxSy() use only +the first five, wo, c, d, A and B, and silently ignore C and D. + +daostorm_3d.z_calibration defaults to --fit_order 2, which produces +exactly those five, so this only matters if you fit a calibration with +fit_order 3 or 4. Such a calibration is followed in full by the "3d" +model's z calculation and truncated everywhere else. """ diff --git a/storm_analysis/sa_utilities/fitz_c.py b/storm_analysis/sa_utilities/fitz_c.py index 879feb4c..8dc6b00f 100644 --- a/storm_analysis/sa_utilities/fitz_c.py +++ b/storm_analysis/sa_utilities/fitz_c.py @@ -44,6 +44,10 @@ def calcSxSy(wx_params, wy_params, z): """ Return sigma x and sigma y given the z calibration parameters. + + This uses the first five terms only, wo, c, d, A and B. A calibration + fitted with fit_order 3 or 4 also has C and D terms, and fitz.c does + use them, so the two do not agree on such a calibration. """ zx = (z - wx_params[1])/wx_params[2] sx = 0.5 * wx_params[0] * numpy.sqrt(1.0 + zx*zx + wx_params[3]*zx*zx*zx + wx_params[4]*zx*zx*zx*zx) From 4deff6215f5c62d1494417eb7789b8a3beb9bb21 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 14:38:03 -0400 Subject: [PATCH 11/13] Say that DriftFromFile takes z in microns The doc string said nanometers. drift() adds the third column straight onto h5_data["z"], which the PSF classes document as microns, so the file has to be microns too. The callers settle it. diagnostics/spliner_configure_common.py and diagnostics/measure_psf both build dz = numpy.arange(-spline_z_range, spline_z_range + 0.001, 0.01) and write that same array into two files: column 2 of drift.txt, and column 1 of z_offset.txt. The z offset file is documented, in measure_psf.py's own help text, as "the z offset of the frame relative to the focal plane in microns". One array cannot be in two different units, so drift.txt is microns. The values agree. spline_z_range is 0.75, which as microns is a plus or minus 750nm stack in 10nm steps, and as nanometers would be plus or minus 0.75nm in steps of ten picometers. Following the doc string would put every z offset out by 1000x. Also says what the file's three columns are, which was not written down anywhere. Co-Authored-By: Claude Opus 5 --- storm_analysis/simulator/drift.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/storm_analysis/simulator/drift.py b/storm_analysis/simulator/drift.py index e4ca8268..a2c3b73b 100644 --- a/storm_analysis/simulator/drift.py +++ b/storm_analysis/simulator/drift.py @@ -22,8 +22,9 @@ def __init__(self, sim_fp, x_size, y_size, h5_data): class DriftFromFile(Drift): """ - Add drift from a file. X and Y are in units of pixels, Z is in - nanometers. + Add drift from a file. The file has one row per frame and three + columns, x, y and z. X and Y are in units of pixels and Z is in + microns, the same units as the localizations it is added to. """ def __init__(self, sim_fp, x_size, y_size, h5_data, drift_file): super(DriftFromFile, self).__init__(sim_fp, x_size, y_size, h5_data) From 5fd112143ed4419045c123ad86d0a814317c142f Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 14:39:46 -0400 Subject: [PATCH 12/13] Name the right axis in the maxima finder's assertions findMaxima() sets the sizes correctly and then complains about the wrong axis four lines later: self.flm_data.xsize = int(images[0].shape[1]) self.flm_data.ysize = int(images[0].shape[0]) ... assert (image.shape[0] == ...), "Unexpected image x size!" assert (image.shape[1] == ...), "Unexpected image y size!" The assignments are right, and match the C, which indexes as images[zi][yi*xsize + xi], so x is the fast axis and xsize is shape[1]. The messages say the opposite. Message text only. The assertions compare an image against taken, so they are checking the right things and would fire correctly. But anyone who trips one is told the opposite of which axis is wrong, in the module that defines the x/y convention for the whole peak finder, and this package has two conventions in play: the simulator works in [x, y] order internally and transposes on write, see the comment in simulator/camera.py. Note that admm/admm_math.py also labels shape[0] as X, but it is self-consistent about it, mx is shape[0] and my is shape[1] throughout, and it describes matrix cells rather than images, so it is left alone. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/ia_utilities_c.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/storm_analysis/sa_library/ia_utilities_c.py b/storm_analysis/sa_library/ia_utilities_c.py index cc3c4228..307df2d3 100644 --- a/storm_analysis/sa_library/ia_utilities_c.py +++ b/storm_analysis/sa_library/ia_utilities_c.py @@ -197,8 +197,8 @@ def findMaxima(self, images, want_height = False): # Verify images are the correct size and C contiguous. for image in images: - assert (image.shape[0] == self.taken[0].shape[0]), "Unexpected image x size!" - assert (image.shape[1] == self.taken[0].shape[1]), "Unexpected image y size!" + assert (image.shape[0] == self.taken[0].shape[0]), "Unexpected image y size!" + assert (image.shape[1] == self.taken[0].shape[1]), "Unexpected image x size!" assert (image.flags['C_CONTIGUOUS']), "Image is not C contiguous!" # Create pointer array to the images. From 95e55d4756c8b8895ccca5164afc8deea143e240 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 14:41:29 -0400 Subject: [PATCH 13/13] Name the right derivative in the spline accessor comments Each d?fAt?D() comment block opens with a summary and closes with a "Return" line naming which derivative it gives. Four of those were wrong, in three of the five functions: dxfAt3D Return - The derivative in y. should be x dyfAt3D "the derivative of the in spline y" word order dzfAt3D "the derivative of the in spline z" word order dzfAt3D Return - The derivative in y. should be z so of the five accessors only dxfAt2D and dyfAt2D described themselves correctly, and two of the three 3D ones claimed to return the y derivative. These are the accessors underneath the Spliner Jacobian, so someone checking the derivative terms against the comments finds two of them apparently returning the wrong quantity, which is exactly the kind of doubt that costs an afternoon. The parallel d?fSpline?D() family, which takes double coordinates rather than integer ones, is correct throughout and is left alone. Comment only: cubic_spline.o is byte identical before and after, md5 b0b37b59a7cd4270e7014183ed032699. Co-Authored-By: Claude Opus 5 --- storm_analysis/spliner/cubic_spline.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/storm_analysis/spliner/cubic_spline.c b/storm_analysis/spliner/cubic_spline.c index 0b5ac484..443284d8 100644 --- a/storm_analysis/spliner/cubic_spline.c +++ b/storm_analysis/spliner/cubic_spline.c @@ -208,7 +208,7 @@ double dxfAt2D(splineData *spline_data, int yc, int xc) * yc - The y coordinate (integer). * xc - The x coordinate (integer). * - * Return - The derivative in y. + * Return - The derivative in x. */ double dxfAt3D(splineData *spline_data, int zc, int yc, int xc) { @@ -322,7 +322,7 @@ double dyfAt2D(splineData *spline_data, int yc, int xc) /* * dyfAt3D() * - * Compute the derivative of the in spline y at x,y,z coordinate (integer). + * Compute the derivative of the spline in y at x,y,z coordinate (integer). * In order for this to work correctly computeDelta3D should have already * been called. * @@ -418,7 +418,7 @@ double dyfSpline3D(splineData *spline_data, double z, double y, double x) /* * dzfAt3D() * - * Compute the derivative of the in spline z at x,y,z coordinate (integer). + * Compute the derivative of the spline in z at x,y,z coordinate (integer). * In order for this to work correctly computeDelta3D should have already * been called. * @@ -427,7 +427,7 @@ double dyfSpline3D(splineData *spline_data, double z, double y, double x) * yc - The y coordinate (integer). * xc - The x coordinate (integer). * - * Return - The derivative in y. + * Return - The derivative in z. */ double dzfAt3D(splineData *spline_data, int zc, int yc, int xc) {