From 3acb33f014f51370d3e971762976d7eb97f9f205 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 3 Sep 2026 11:13:58 -0400 Subject: [PATCH 1/9] Use yoff, not xoff, for the y peak shape offset daoCalcPeakShape() computed the y offset of the Gaussian within the ROI from fit_data->xoff: yo = -(peak->params[YCENTER] - (double)peak->yi) - fit_data->xoff; daoInitialize() assigns both offsets 0.5*roi_size - 1.0, in both of its branches, so xoff == yoff identically and nothing observable changes today. That is a property of this one code path rather than of the codebase: spliner/cubic_fit.c derives xoff from sx and yoff from sy, which differ for a non-square spline. Working through it, the fixed line gives yo = yi - Y for any yoff, where Y is the position in the image, while the original leaves an extra (yoff - xoff) behind. Built the library four ways and fit a single emitter at x = 25.3, y = 13.7 to confirm: recovered x recovered y original, yoff == xoff 25.3000 13.7000 fixed, yoff == xoff 25.3000 13.7000 original, yoff = xoff + 1 25.3000 14.7000 fixed, yoff = xoff + 1 25.3000 13.7000 The fix changes nothing as the code stands and makes the y position independent of the offset convention, which is what the x line already does. The perturbed yoff was a temporary local patch for the measurement and is not part of this commit. No regression test. With xoff == yoff there is no input that tells the two versions apart, and nothing in the Python API can set them to different values. A test here could only pin the coincidence that makes the bug invisible. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/dao_fit.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/storm_analysis/sa_library/dao_fit.c b/storm_analysis/sa_library/dao_fit.c index 0be46b66..6c671d2b 100644 --- a/storm_analysis/sa_library/dao_fit.c +++ b/storm_analysis/sa_library/dao_fit.c @@ -983,7 +983,7 @@ void daoCalcPeakShape(fitData *fit_data) dao_peak = (daoPeak *)peak->peak_model; xo = -(peak->params[XCENTER] - (double)peak->xi) - fit_data->xoff; - yo = -(peak->params[YCENTER] - (double)peak->yi) - fit_data->xoff; + yo = -(peak->params[YCENTER] - (double)peak->yi) - fit_data->yoff; /* Calculate peak shape. */ for(i=0;ifit_size_x;i++){ From 0847a07d6896c68f02c6716cbd4d495c768be178 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 3 Sep 2026 11:17:48 -0400 Subject: [PATCH 2/9] Finish removing the sensitivity corrected foreground sum mFitGetPeakPropertyDouble() answered a request for "fg_sum_sc" with the uncorrected sum, the intended call sitting commented out beside it: else if (!strcmp(what, "fg_sum_sc")){ for(i=0;infit;i++){ values[i] = mFitPeakFgSum(fit_data, &fit_data->fit[i]); // values[i] = mFitPeakFgSumSensitivityCorrected(...); } } and mFitPeakFgSumSensitivityCorrected() was declared in multi_fit.h with no definition anywhere in the tree. This is not an open question, it is an unfinished removal. 0fe8987f added sensitivity correction; a1520d61 took it out again, saying "Remove 'sensitivity_correction' parameter, decided this option wasn't worth the complexity", and left the branch and the declaration behind. Delete both. Note that the property was already unreachable from Python. The peak_properties dict in dao_fit_c.py lists fg_sum but not fg_sum_sc, and getPeakProperty() raises MultiFitterException for anything not in it. A C caller now falls through to the function's own else clause, which prints "Unrecognized parameter 'fg_sum_sc'!" and leaves the values alone, rather than quietly answering a different question. No regression test. The branch was unreachable from Python, so nothing can observe its removal, and a test asserting that getPeakProperty() rejects the name would pass on the parent commit too. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/multi_fit.c | 6 ------ storm_analysis/sa_library/multi_fit.h | 1 - 2 files changed, 7 deletions(-) diff --git a/storm_analysis/sa_library/multi_fit.c b/storm_analysis/sa_library/multi_fit.c index 2aeb697c..54c94be5 100644 --- a/storm_analysis/sa_library/multi_fit.c +++ b/storm_analysis/sa_library/multi_fit.c @@ -777,12 +777,6 @@ void mFitGetPeakPropertyDouble(fitData *fit_data, double *values, char *what) values[i] = mFitPeakFgSum(fit_data, &fit_data->fit[i]); } } - else if (!strcmp(what, "fg_sum_sc")){ - for(i=0;infit;i++){ - values[i] = mFitPeakFgSum(fit_data, &fit_data->fit[i]); - // values[i] = mFitPeakFgSumSensitivityCorrected(fit_data, &fit_data->fit[i]); - } - } else if (!strcmp(what, "height")){ for(i=0;infit;i++){ values[i] = fit_data->fit[i].params[HEIGHT]; diff --git a/storm_analysis/sa_library/multi_fit.h b/storm_analysis/sa_library/multi_fit.h index 29afacf6..2c26fc91 100644 --- a/storm_analysis/sa_library/multi_fit.h +++ b/storm_analysis/sa_library/multi_fit.h @@ -178,7 +178,6 @@ void mFitNewImage(fitData *, double *); void mFitNewPeaks(fitData *, int); double mFitPeakBgSum(fitData *, peakData *); double mFitPeakFgSum(fitData *, peakData *); -double mFitPeakFgSumSensitivityCorrected(fitData *, peakData *); double mFitPeakSum(fitData *, peakData *); void mFitRecenterPeaks(fitData *); void mFitRemoveErrorPeaks(fitData *); From 6302c5e05d4f32e1b619937e3eb23cbd4cb81301 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 3 Sep 2026 11:26:02 -0400 Subject: [PATCH 3/9] Mask the last row and column of a rectangular AOI peakMask() ended both of its stop slices at -1: peak_mask[:,parameters.getAttr("x_stop")+margin:-1] = 0.0 peak_mask[parameters.getAttr("y_stop")+margin:-1,:] = 0.0 which stops one element short, so the final column and the final row kept the value 1.0 rather than being masked out. The two start slices, and the circular AOI branch above them, are all written correctly. Nothing changes in practice. 'shape' here is the padded image shape and the stops are already offset by 'margin', so the row and column in question are padding, and the maxima finder excludes them independently through its own margin. The full test suite is unchanged either way. The slice was still wrong, and it is the kind of wrong that only stays harmless while two unrelated margins happen to agree. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/fitting.py | 4 ++-- storm_analysis/test/test_fitting.py | 28 ++++++++++++++++++++++++++++ 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/storm_analysis/sa_library/fitting.py b/storm_analysis/sa_library/fitting.py index ae1d6e82..c30cf27a 100644 --- a/storm_analysis/sa_library/fitting.py +++ b/storm_analysis/sa_library/fitting.py @@ -210,11 +210,11 @@ def peakMask(shape, parameters, margin): if parameters.hasAttr("x_start"): peak_mask[:,0:parameters.getAttr("x_start")+margin] = 0.0 if parameters.hasAttr("x_stop"): - peak_mask[:,parameters.getAttr("x_stop")+margin:-1] = 0.0 + peak_mask[:,parameters.getAttr("x_stop")+margin:] = 0.0 if parameters.hasAttr("y_start"): peak_mask[0:parameters.getAttr("y_start")+margin,:] = 0.0 if parameters.hasAttr("y_stop"): - peak_mask[parameters.getAttr("y_stop")+margin:-1,:] = 0.0 + peak_mask[parameters.getAttr("y_stop")+margin:,:] = 0.0 return peak_mask diff --git a/storm_analysis/test/test_fitting.py b/storm_analysis/test/test_fitting.py index 93f8300e..da4b1a15 100644 --- a/storm_analysis/test/test_fitting.py +++ b/storm_analysis/test/test_fitting.py @@ -84,6 +84,33 @@ def test_fitting_peak_mask_square_aoi(): assert (numpy.count_nonzero(mask[:10+margin,:]) == 0) +def test_fitting_peak_mask_stop(): + """ + The x_stop and y_stop edges of a rectangular AOI. + + Both slices ended at -1, which stops one element short, so the last + column and the last row kept the value 1.0 instead of being masked. + """ + margin = 5 + [x_stop, y_stop] = [70, 30] + + p = params.ParametersDAO() + p.setAttr("x_stop", "int", x_stop) + p.setAttr("y_stop", "int", y_stop) + + # More columns than rows, so the two axes cannot be confused. + shape = (60, 120) + mask = fitting.peakMask(shape, p, margin) + + # Nothing at or beyond the stops survives, including the final + # column and the final row. + assert (numpy.count_nonzero(mask[:,x_stop+margin:]) == 0) + assert (numpy.count_nonzero(mask[y_stop+margin:,:]) == 0) + + # And the region inside both stops is untouched. + assert (numpy.all(mask[:y_stop+margin,:x_stop+margin] == 1.0)) + + def test_fitting_check_mode(): """ check_mode writes diagnostic images from the arbitrary PSF peak finder, @@ -125,4 +152,5 @@ def test_fitting_check_mode(): test_fitting_no_peak_finder() test_fitting_peak_mask_circular() test_fitting_peak_mask_square_aoi() + test_fitting_peak_mask_stop() test_fitting_check_mode() From 071124673fb79504d80f6b022e2d1127a5bb5858 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 3 Sep 2026 11:32:50 -0400 Subject: [PATCH 4/9] Stop CSpline2D.py_f from transposing its arguments Spline2D.f() takes (y, x), and so does CSpline2D.py_f(), but it passed them on the other way round: def py_f(self, y, x): return self.py_spline.f(x, y) so it evaluated the transposed point. CSpline3D.py_f() hands (z, y, x) straight through and has always been right. py_f() exists to evaluate in Python the same point that f() evaluates in C, which is how you check one implementation against the other. A PSF spline is not symmetric in x and y, so the two answers genuinely differ: on test_spliner_psf_2d.spline at y = 5.0, x = 6.5, f() returns 0.003209 and py_f() returned 0.002391. Nothing calls py_f(), which is why this has never mattered. It is also why it is worth fixing: the one place it would be used is validating the C library, and there it would report a discrepancy that is not there. The new tests make py_f() do the job it was written for, comparing it against f() over a thousand random points. The 2D test fails on the parent commit; the 3D one is the control and passes on both. With the fix, py_f() and f() agree to 3.5e-18. Co-Authored-By: Claude Opus 5 --- storm_analysis/spliner/cubic_spline_c.py | 2 +- storm_analysis/test/test_spline.py | 50 ++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/storm_analysis/spliner/cubic_spline_c.py b/storm_analysis/spliner/cubic_spline_c.py index c6daad51..e87ef557 100644 --- a/storm_analysis/spliner/cubic_spline_c.py +++ b/storm_analysis/spliner/cubic_spline_c.py @@ -143,7 +143,7 @@ def getPSF(self, dy, dx): return psf def py_f(self, y, x): - return self.py_spline.f(x, y) + return self.py_spline.f(y, x) class CSpline3D(CSpline): diff --git a/storm_analysis/test/test_spline.py b/storm_analysis/test/test_spline.py index 23a192a4..172bf7b5 100644 --- a/storm_analysis/test/test_spline.py +++ b/storm_analysis/test/test_spline.py @@ -171,6 +171,54 @@ def test_psf_3D_dz(): assert (abs(py_spline.dzf(x, y, z) - c_spline.dzf(x, y, z)) < 1.0e-6) +def test_psf_2D_py_f(): + """ + py_f() is the Python evaluation of the same point that f() evaluates in + C, so the two have to agree. CSpline2D.py_f() passed its arguments on in + the opposite order, which evaluates the transposed point. A PSF spline is + not symmetric in x and y, so the answers differ. + """ + if (sys.version_info < (3, 0)): + return + + spline_filename = storm_analysis.getData("test/data/test_spliner_psf_2d.spline") + with open(spline_filename, "rb") as fp: + spline_data = pickle.load(fp) + + py_spline = spline2D.Spline2D(spline_data["spline"], spline_data["coeff"]) + c_spline = cubicSplineC.CSpline2D(py_spline) + + size = py_spline.getSize() - 1.0e-6 + + for i in range(reps): + x = random.uniform(1.0e-6, size) + y = random.uniform(1.0e-6, size) + assert (abs(c_spline.py_f(y, x) - c_spline.f(y, x)) < 1.0e-6) + +def test_psf_3D_py_f(): + """ + The same check on CSpline3D, which has always passed its arguments + through in the right order. + """ + if (sys.version_info < (3, 0)): + return + + spline_filename = storm_analysis.getData("test/data/test_spliner_psf.spline") + with open(spline_filename, "rb") as fp: + spline_data = pickle.load(fp) + + py_spline = spline3D.Spline3D(spline_data["spline"], spline_data["coeff"]) + c_spline = cubicSplineC.CSpline3D(py_spline) + + size = py_spline.getSize() - 1.0e-6 + + for i in range(reps): + x = random.uniform(1.0e-6, size) + y = random.uniform(1.0e-6, size) + z = random.uniform(1.0e-6, size) + assert (abs(c_spline.py_f(z, y, x) - c_spline.f(z, y, x)) < 1.0e-6) + + if (__name__ == "__main__"): test_psf_2D_f() test_psf_2D_dx() @@ -179,3 +227,5 @@ def test_psf_3D_dz(): test_psf_3D_dx() test_psf_3D_dy() test_psf_3D_dz() + test_psf_2D_py_f() + test_psf_3D_py_f() From c919815009055b2527ff2c56e1bc102a270ad38f Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 3 Sep 2026 11:39:19 -0400 Subject: [PATCH 5/9] Fail clearly on an incomplete .inf or a missing .dax DaxReader gives the frame dimensions a documented fallback, 256x256 with a printed warning, but two other things it reads from the .inf file had none, and neither did the movie itself: complete .inf and .dax ok .inf with no frame count AttributeError: 'DaxReader' object has no attribute 'number_frames' .inf with no endian line AttributeError: ... no attribute 'bigendian', from inside loadAFrame() .dax missing constructs, then AttributeError: 'NoneType' object has no attribute 'seek' None of these put correctness at stake, but a truncated .inf or a moved .dax surfaced as an AttributeError deep in the call stack rather than a message naming the file. The missing movie was the worst of the three: the warning was gated behind 'verbose', so by default the reader constructed silently into an unusable state. Initialize bigendian and number_frames alongside the image size, then: - endianness defaults to little, with a warning. Everything this package writes is little endian. - a missing .dax raises IOError naming the file, which is what inferReader() already does for an unsupported extension. - the frame count comes from the file size. A dax file is a bare block of 16 bit pixels with no header, so this is exact rather than a guess, and raising above guarantees the file is there to measure. Note that raising is a behaviour change: code that constructed a reader for a file that does not exist, and never read from it, now fails. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/datareader.py | 26 ++++++++--- storm_analysis/test/test_dataio.py | 61 ++++++++++++++++++++++++- 2 files changed, 80 insertions(+), 7 deletions(-) diff --git a/storm_analysis/sa_library/datareader.py b/storm_analysis/sa_library/datareader.py index a2132e99..e55d83b7 100644 --- a/storm_analysis/sa_library/datareader.py +++ b/storm_analysis/sa_library/datareader.py @@ -167,8 +167,10 @@ def __init__(self, filename, verbose = False): self.inf_filename = dirname + os.path.splitext(os.path.basename(filename))[0] + ".inf" # defaults + self.bigendian = None self.image_height = None self.image_width = None + self.number_frames = None # extract the movie information from the associated inf file size_re = re.compile(r'frame dimensions = ([\d]+) x ([\d]+)') @@ -217,17 +219,29 @@ def __init__(self, filename, verbose = False): # set defaults, probably correct, but warn the user # that they couldn't be determined from the inf file. - if not self.image_height: + if self.image_height is None: print("Could not determine image size, assuming 256x256.") self.image_height = 256 self.image_width = 256 + # Not every .inf file records the endianness. Everything this package + # writes is little endian, so that is the better guess. + if self.bigendian is None: + print("Could not determine endianness of", self.inf_filename, "assuming little endian.") + self.bigendian = 0 + # open the dax file - if os.path.exists(filename): - self.fileptr = open(filename, "rb") - else: - if self.verbose: - print("dax data not found", filename) + if not os.path.exists(filename): + raise IOError("dax data not found '" + filename + "'") + self.fileptr = open(filename, "rb") + + # A dax file is a bare block of 16 bit pixels with no header, so if the + # .inf file did not say how many frames there are the file size does. + if self.number_frames is None: + frame_size = 2 * self.image_height * self.image_width + self.number_frames = os.path.getsize(filename)//frame_size + print("Could not determine the number of frames in", self.inf_filename + ",", + "the file size gives", self.number_frames) def loadAFrame(self, frame_number): """ diff --git a/storm_analysis/test/test_dataio.py b/storm_analysis/test/test_dataio.py index fa8f1d52..f829545d 100644 --- a/storm_analysis/test/test_dataio.py +++ b/storm_analysis/test/test_dataio.py @@ -211,10 +211,69 @@ def test_io_6(): assert(ml == 1) assert(numpy.allclose(data, rd.loadAFrame(0))) - +def test_io_7(): + """ + DaxReader on an incomplete .inf file, and on a missing .dax. + + 'number of frames' and the endianness had no fallback, so a .inf file + without them gave AttributeError from inside loadAFrame(). A missing + .dax constructed a reader that then failed with AttributeError on a + None file pointer. + """ + import pytest + + movie_h = 8 + movie_w = 16 + movie_l = 5 + + data = numpy.random.randint(0, 60000, (movie_h, movie_w)).astype(numpy.uint16) + + movie_name = storm_analysis.getPathOutputTest("test_dataio_inf.dax") + inf_name = storm_analysis.getPathOutputTest("test_dataio_inf.inf") + + wr = datawriter.inferWriter(movie_name) + for i in range(movie_l): + wr.addFrame(data) + wr.close() + + with open(inf_name) as fp: + inf_lines = fp.readlines() + + def writeInf(drop): + with open(inf_name, "w") as fp: + for line in inf_lines: + if (drop is None) or (not drop in line): + fp.write(line) + + ## No 'number of frames', the file size gives it instead. + writeInf("number of frames") + rd = datareader.inferReader(movie_name) + assert(rd.filmSize() == [movie_w, movie_h, movie_l]) + assert(numpy.allclose(data, rd.loadAFrame(0))) + rd.close() + + ## No endianness, little endian is assumed. + writeInf("data type") + rd = datareader.inferReader(movie_name) + assert(numpy.allclose(data, rd.loadAFrame(0))) + rd.close() + + ## A complete .inf but no .dax, and the error names the file. + missing_name = storm_analysis.getPathOutputTest("test_dataio_missing.dax") + storm_analysis.removeFile(missing_name) + with open(storm_analysis.getPathOutputTest("test_dataio_missing.inf"), "w") as fp: + for line in inf_lines: + fp.write(line) + + with pytest.raises(IOError) as err: + datareader.inferReader(missing_name) + assert("test_dataio_missing.dax" in str(err.value)) + + if (__name__ == "__main__"): test_io_1() test_io_2() test_io_3() test_io_6() + test_io_7() From 50e2c56ad0b1dabe1a49014dd842f01e8dc6b806 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 3 Sep 2026 11:48:15 -0400 Subject: [PATCH 6/9] Return a list from the v0 calibration branch loadCMOSCalibration() says it returns "the CMOS calibration data as a list [offset, variance, gain, rqe]". The v1 and v2 branches do. The v0 branch returned return map(numpy.transpose, data) which in Python 3 is a lazy iterator, not a list. It cannot be indexed, and it is exhausted by a single pass. Every caller in the tree tuple unpacks the result, and that works on an iterator, which is why this has gone unnoticed. It breaks for any caller that indexes the result or iterates it twice, and only for users on legacy v0 calibration files, so it would look like a problem with their data rather than with the reader. test_cal_v0 already covered the v0 path but only by unpacking, which an iterator satisfies. It now checks what the doc string actually promises. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/analysis_io.py | 2 +- storm_analysis/test/test_analysis_io.py | 12 +++++++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/storm_analysis/sa_library/analysis_io.py b/storm_analysis/sa_library/analysis_io.py index 2de99593..b8c36526 100644 --- a/storm_analysis/sa_library/analysis_io.py +++ b/storm_analysis/sa_library/analysis_io.py @@ -54,7 +54,7 @@ def loadCMOSCalibration(filename, verbose = False): rqe = numpy.ones_like(data[0]) data = list(data) + [rqe] - return map(numpy.transpose, data) + return list(map(numpy.transpose, data)) # Check for v1 format. elif (len(data) == 4): diff --git a/storm_analysis/test/test_analysis_io.py b/storm_analysis/test/test_analysis_io.py index 87fb9cc5..7842a119 100644 --- a/storm_analysis/test/test_analysis_io.py +++ b/storm_analysis/test/test_analysis_io.py @@ -28,7 +28,17 @@ def test_cal_v0(): numpy.transpose(variance), numpy.transpose(gain)], fp) - [o, v, g, r] = analysisIO.loadCMOSCalibration(cal_file) + cal = analysisIO.loadCMOSCalibration(cal_file) + + # The doc string promises a list, and the v1 and v2 branches return one. + # v0 returned map(), which is not indexable and is exhausted by a single + # pass. + assert(isinstance(cal, list)) + assert(cal[0].shape == cal_size) + assert(len(list(cal)) == 4) + assert(len(list(cal)) == 4) + + [o, v, g, r] = cal assert(numpy.allclose(offset, o)) assert(numpy.allclose(variance, v)) assert(numpy.allclose(gain, g)) From 5919fe607c3531672527f9ca7a006920ad16d4f6 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 3 Sep 2026 11:58:36 -0400 Subject: [PATCH 7/9] Document the drift_corrected defaults The three localization read paths disagree about drift correction, and none of them said so: getLocalizations(drift_corrected = False) no doc string getLocalizationsInFrame(drift_corrected = False) no doc string getLocalizationsInFrameRange(drift_corrected = False) doc string, silent localizationsIterator(drift_corrected = True) doc string, silent So switching between "load them all" and "iterate over them" changes whether dx, dy and dz have been applied, with nothing to warn you. This is the mechanism behind the marker bug fixed in 0377f488: zCheck() read corrected z and compared it against a marker written in raw z. The defaults are not wrong. Each suits its method's primary use: the getters hand back what is stored, and the iterator is described as the easiest way to process a whole file, where corrected coordinates are what you want. So this is a documentation change, not a behaviour one. State the default and the reason in all four, cross referencing so that someone reading one sees the other. Also note on tracksIterator() why it has no such argument: tracks are built from corrected localizations, so their coordinates are already corrected and the raw values are gone. test_sa_h5py_20 pins the defaults. It cannot fail on the parent commit, since nothing behavioural changes here, but it does fail if the defaults are ever harmonized: flipping getLocalizations() to True fails it. Given that documentation is the whole of this fix, that test is the only part of it that resists being undone. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/sa_h5py.py | 29 ++++++++++++++++ storm_analysis/test/test_sa_h5py.py | 49 ++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+) diff --git a/storm_analysis/sa_library/sa_h5py.py b/storm_analysis/sa_library/sa_h5py.py index 85a68cb4..4de1fefb 100644 --- a/storm_analysis/sa_library/sa_h5py.py +++ b/storm_analysis/sa_library/sa_h5py.py @@ -338,13 +338,28 @@ def getGroupName(self, frame_number): return "/fr_" + str(frame_number) def getLocalizations(self, drift_corrected = False, fields = None): + """ + Return all of the localizations in the file as a single dictionary + of numpy arrays. + + drift_corrected defaults to False here, so x, y and z come back + exactly as they are stored. Note that localizationsIterator() + defaults to True, so the two ways of reading a file do not give you + the same coordinates unless you say which you want. + """ return self.getLocalizationsInFrameRange(0, self.hdf5.attrs['movie_l'], drift_corrected = drift_corrected, fields = fields) def getLocalizationsInFrame(self, frame_number, drift_corrected = False, fields = None): + """ + Return the localizations in a single frame. + This is where the drift correction is actually applied, by adding + the frame's dx, dy and dz attributes. drift_corrected defaults to + False, so the stored values are returned unchanged. + """ locs = {} grp = self.getGroup(frame_number) @@ -364,6 +379,9 @@ def getLocalizationsInFrame(self, frame_number, drift_corrected = False, fields def getLocalizationsInFrameRange(self, start, stop, drift_corrected = False, fields = None): """ Return the localizations in the range start <= frame number < stop. + + As with getLocalizations(), drift_corrected defaults to False and + the stored coordinates are returned unchanged. """ assert(stop > start) locs = {} @@ -528,6 +546,13 @@ def localizationsIterator(self, drift_corrected = True, fields = None, skip_empt If skip_empty is false you will get empty locs dictionaries for frames that have no localizations. + + Note that drift_corrected defaults to True here, unlike + getLocalizations() and getLocalizationsInFrame(), which default to + False. Corrected coordinates are usually what you want when you are + just walking a finished file, but it does mean that switching + between loading and iterating silently changes whether dx, dy and dz + have been applied. """ # This should be a zero length generator if there are no localizations. for i in range(self.getMovieLength()): @@ -605,6 +630,10 @@ def tracksIterator(self, fields = None): The track center in x,y,z are the 'x', 'y' and 'z' fields. The other fields may need to be normalized by the track length. + + There is no drift_corrected argument because tracks are built from + drift corrected localizations, so their coordinates are already + corrected and there is no way to get back the raw values. """ # This should be a zero length generator if there are no tracks. if (not self.hasTracks()): diff --git a/storm_analysis/test/test_sa_h5py.py b/storm_analysis/test/test_sa_h5py.py index e29e0b1b..ed039340 100644 --- a/storm_analysis/test/test_sa_h5py.py +++ b/storm_analysis/test/test_sa_h5py.py @@ -572,6 +572,54 @@ def test_sa_h5py_19(): assert not elt in locs +def test_sa_h5py_20(): + """ + Pin the drift_corrected defaults. + + getLocalizations() and getLocalizationsInFrame() default to False and + return the stored coordinates. localizationsIterator() defaults to True + and returns corrected ones. The defaults suit each method's primary use, + but they disagree, so they are worth holding still: harmonizing them + would silently change what every existing caller gets back. + """ + [dx, dy, dz] = [1.0, -2.0, 0.5] + + peaks = {"x" : numpy.zeros(3), + "y" : numpy.zeros(3), + "z" : numpy.zeros(3)} + + h5_name = storm_analysis.getPathOutputTest("test_sa_hdf5_drift.hdf5") + storm_analysis.removeFile(h5_name) + + with saH5Py.SAH5Py(h5_name, is_existing = False, overwrite = True) as h5: + h5.setMovieInformation(100, 100, 1, "") + h5.addLocalizations(peaks, 0) + h5.setDriftCorrection(0, dx = dx, dy = dy, dz = dz) + + with saH5Py.SAH5Py(h5_name) as h5: + + # The getters hand back what is stored. + for locs in [h5.getLocalizations(), + h5.getLocalizationsInFrame(0), + h5.getLocalizationsInFrameRange(0, 1)]: + assert(numpy.allclose(locs["x"], numpy.zeros(3))) + assert(numpy.allclose(locs["y"], numpy.zeros(3))) + assert(numpy.allclose(locs["z"], numpy.zeros(3))) + + # The iterator corrects. + for fnum, locs in h5.localizationsIterator(): + assert(numpy.allclose(locs["x"], dx*numpy.ones(3))) + assert(numpy.allclose(locs["y"], dy*numpy.ones(3))) + assert(numpy.allclose(locs["z"], dz*numpy.ones(3))) + + # And both are explicitly overridable. + locs = h5.getLocalizations(drift_corrected = True) + assert(numpy.allclose(locs["x"], dx*numpy.ones(3))) + + for fnum, locs in h5.localizationsIterator(drift_corrected = False): + assert(numpy.allclose(locs["x"], numpy.zeros(3))) + + if (__name__ == "__main__"): test_sa_h5py_1() test_sa_h5py_2() @@ -592,4 +640,5 @@ def test_sa_h5py_19(): test_sa_h5py_17() test_sa_h5py_18() test_sa_h5py_19() + test_sa_h5py_20() From cc6f128bca86626cf7dae074fa8bc44e17e4d2f0 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 3 Sep 2026 12:01:52 -0400 Subject: [PATCH 8/9] Take z in microns in astigmaticPSF astigmaticPSF.PSF() and PSFIntegral() scaled z by 0.001 before calling calcSxSy(), which says the incoming z is nanometers. Every other PSF module in the simulator takes microns: PupilFunctionBase.getPSFs() and Spline.getPSFs() both document "x, y in pixels and z in microns", and that is what h5_data['z'] holds. The parameters settle which side is wrong. calcSxSy() compares z against wx_params[1] and [2], which are 0.150 and 0.40. As microns that is a 150nm offset and a 400nm defocusing scale, ordinary astigmatism. As nanometers it would be 0.15nm and 0.4nm, which is meaningless. So the parameters are microns and the conversion was backwards. Fed the simulator's own convention, every z landed within a nanometer of focus: z (um) with z*0.001 z used directly -0.40 sx=1.068 sy=1.068 sx=1.700 sy=1.179 -0.20 sx=1.068 sy=1.068 sx=1.329 sy=1.008 0.00 sx=1.068 sy=1.068 sx=1.068 sy=1.068 +0.20 sx=1.068 sy=1.068 sx=1.008 sy=1.329 +0.40 sx=1.068 sy=1.068 sx=1.179 sy=1.700 A round, z-independent PSF at every z. Nothing calls this module, which is the only reason it has been dormant: wired into the simulator the way its siblings are it would silently produce data with no astigmatism and no z dependence, and fail nothing. That is exactly the degenerate case that made the first Z fitter measurement in this audit meaningless, round fixed-width spots that a Z model test cannot tell from correct output. Checked that nothing depends on the old convention. No test imports storm_analysis.diagnostics, the diagnostics scripts build their PSFs from psf.PupilFunction, psf.GaussianPSF or a spline file, and across the whole repository the only references to astigmaticPSF are the new test and one line of simulator/README.txt. calcSxSy() itself is unchanged. Note that the sibling dhPSF.py has z_min and z_max of -500 and 500, which are nanometers, so the two hand written PSF modules now disagree with each other. dhPSF is equally uncalled and is left alone here. Co-Authored-By: Claude Opus 5 --- storm_analysis/simulator/astigmaticPSF.py | 11 ++++-- storm_analysis/test/test_simulator.py | 41 +++++++++++++++++++++++ 2 files changed, 49 insertions(+), 3 deletions(-) diff --git a/storm_analysis/simulator/astigmaticPSF.py b/storm_analysis/simulator/astigmaticPSF.py index fca5343e..9fc7a77d 100644 --- a/storm_analysis/simulator/astigmaticPSF.py +++ b/storm_analysis/simulator/astigmaticPSF.py @@ -4,6 +4,9 @@ of objects that emulates an astigmatic PSF in a format compatible with drawgaussians. +x and y are in pixels and z is in microns, which is the convention the +rest of the simulator uses and what h5_data['z'] holds. + Hazen 01/16 """ @@ -13,7 +16,9 @@ psf_type = "astigmatic" -# Parameters for astigmatic PSF. +# Parameters for astigmatic PSF. calcSxSy() compares z against elements 1 +# and 2, so those are an offset of 150nm and a defocusing scale of 400nm, +# both in microns like z. wx_params = numpy.array([2.0, 0.150, 0.40, 0.0, 0.0]) wy_params = numpy.array([2.0, -0.150, 0.40, 0.0, 0.0]) @@ -21,7 +26,7 @@ def PSF(x, y, z, h): num_objects = x.size objects = numpy.zeros((num_objects, 5)) for i in range(num_objects): - [sx, sy] = fitzC.calcSxSy(wx_params, wy_params, z[i] * 0.001) + [sx, sy] = fitzC.calcSxSy(wx_params, wy_params, z[i]) objects[i,:] = [x[i], y[i], h[i], sx, sy] return objects @@ -29,7 +34,7 @@ def PSF(x, y, z, h): def PSFIntegral(z, h): integral = numpy.zeros(z.size) for i in range(z.size): - [sx, sy] = fitzC.calcSxSy(wx_params, wy_params, z[i] * 0.001) + [sx, sy] = fitzC.calcSxSy(wx_params, wy_params, z[i]) integral[i] = 2.0 * numpy.pi * h[i] * sx * sy return integral diff --git a/storm_analysis/test/test_simulator.py b/storm_analysis/test/test_simulator.py index 8d474cd0..4c7743a5 100644 --- a/storm_analysis/test/test_simulator.py +++ b/storm_analysis/test/test_simulator.py @@ -7,6 +7,7 @@ import storm_analysis.spliner.spline_to_psf as splineToPSF +import storm_analysis.simulator.astigmaticPSF as astigmaticPSF import storm_analysis.simulator.background as background import storm_analysis.simulator.camera as camera import storm_analysis.simulator.photophysics as photophysics @@ -149,6 +150,45 @@ def test_simulate_3(): sim.simulate(dax_name, bin_name, 5) +def test_astigmatic_psf(): + """ + astigmaticPSF takes z in microns, like everything else in the simulator. + + It used to multiply z by 0.001 before calling calcSxSy(), which says the + incoming z is nanometers. Fed the simulator's own convention that put + every z within a nanometer of focus, so it returned a round, + z-independent PSF, sx = sy = 1.068 px at every z. + """ + z = numpy.array([-0.4, -0.2, 0.0, 0.2, 0.4]) + x = numpy.ones(z.size) + y = numpy.ones(z.size) + h = numpy.ones(z.size)*100.0 + + objects = astigmaticPSF.PSF(x, y, z, h) + sx = objects[:,3] + sy = objects[:,4] + + # Round at focus. + assert(abs(sx[2] - sy[2]) < 1.0e-9) + + # Elongated in x below focus, in y above it, and it is a real effect and + # not rounding. wx_params puts the x waist at +150nm and the y waist at + # -150nm, with a 400nm defocusing scale. + assert(sx[0] > 1.4*sy[0]) + assert(sy[4] > 1.4*sx[4]) + + # Symmetric about focus. + assert(abs(sx[0] - sy[4]) < 1.0e-9) + assert(abs(sy[0] - sx[4]) < 1.0e-9) + + # And the width actually varies with z, which is what a Z fitter needs. + assert(sx[0] > 1.5*sx[2]) + + # PSFIntegral uses the same widths. + integral = astigmaticPSF.PSFIntegral(z, h) + assert(numpy.allclose(integral, 2.0*numpy.pi*h*sx*sy)) + + if (__name__ == "__main__"): test_psf_spline2D_1() test_psf_spline2D_2() @@ -157,3 +197,4 @@ def test_simulate_3(): test_simulate_1() test_simulate_2() test_simulate_3() + test_astigmatic_psf() From b98571c7a9895791bb71d6a1dd9732f42e606dfa Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Thu, 3 Sep 2026 12:13:29 -0400 Subject: [PATCH 9/9] Take z in microns in dhPSF The other hand written PSF module had the same unit problem as astigmaticPSF, in its own idiom. dhPSF has no explicit conversion, but its range constants z_min = -500.0 z_max = 500.0 are nanometers, and PSF() interpolates the lobe angle linearly between them. Fed the simulator's own z in microns, every value sits in the middle of a range a thousand times too wide, so the angle never moves: z (um) dx dy angle -0.50 +1.415 +1.413 45.0 -0.25 +1.415 +1.414 45.0 0.00 +1.414 +1.414 45.0 +0.25 +1.414 +1.415 45.0 +0.50 +1.413 +1.415 45.0 A double helix PSF that does not rotate carries no z information at all, which is the same degeneracy cc6f128b fixed in astigmaticPSF. Fixing only one of the pair would have left them disagreeing with each other. Give z_min and z_max the same physical range in microns, so the lobes sweep 0 to 90 degrees across it as intended, and say the convention in the module doc string. Left alone, both pre-existing: the commented out alternative angle formula, and the question above it about warning when z falls outside the range. The range is still unenforced, it is just in the right units now. Nothing calls this module either. Across the repository the only references to dhPSF are the new test and one line of simulator/README.txt. Co-Authored-By: Claude Opus 5 --- storm_analysis/simulator/dhPSF.py | 10 ++++-- storm_analysis/test/test_simulator.py | 44 +++++++++++++++++++++++++++ 2 files changed, 51 insertions(+), 3 deletions(-) diff --git a/storm_analysis/simulator/dhPSF.py b/storm_analysis/simulator/dhPSF.py index 69822d3d..7c3daa9d 100644 --- a/storm_analysis/simulator/dhPSF.py +++ b/storm_analysis/simulator/dhPSF.py @@ -4,6 +4,9 @@ of objects that emulates a double-helical PSF in a format compatible with drawgaussians. +x and y are in pixels and z is in microns, which is the convention the +rest of the simulator uses and what h5_data['z'] holds. + Hazen 01/16 """ @@ -11,9 +14,10 @@ psf_type = "double-helix" -# Double helix PSF range. -z_min = -500.0 -z_max = 500.0 +# Double helix PSF range, in microns like z. The lobes rotate through 90 +# degrees between z_min and z_max. +z_min = -0.5 +z_max = 0.5 sigma = 1.0 def PSF(x, y, z, h): diff --git a/storm_analysis/test/test_simulator.py b/storm_analysis/test/test_simulator.py index 4c7743a5..9f3e15ba 100644 --- a/storm_analysis/test/test_simulator.py +++ b/storm_analysis/test/test_simulator.py @@ -10,6 +10,7 @@ import storm_analysis.simulator.astigmaticPSF as astigmaticPSF import storm_analysis.simulator.background as background import storm_analysis.simulator.camera as camera +import storm_analysis.simulator.dhPSF as dhPSF import storm_analysis.simulator.photophysics as photophysics import storm_analysis.simulator.psf as psf import storm_analysis.simulator.simulate as simulate @@ -189,6 +190,48 @@ def test_astigmatic_psf(): assert(numpy.allclose(integral, 2.0*numpy.pi*h*sx*sy)) +def test_dh_psf(): + """ + dhPSF takes z in microns too. + + z_min and z_max were -500 and 500, which are nanometers, so feeding it + the simulator's own z froze the helix angle at 45 degrees across the + whole range. A double helix PSF that does not rotate carries no z + information at all, the same degeneracy astigmaticPSF had. + """ + [xc, yc] = [10.0, 20.0] + z = numpy.array([-0.5, 0.0, 0.5]) + + objects = dhPSF.PSF(numpy.ones(z.size)*xc, + numpy.ones(z.size)*yc, + z, + numpy.ones(z.size)*100.0) + + # Two lobes per emitter, ordered [emitter0 +, emitter0 -, emitter1 +, ..]. + assert(objects.shape[0] == 2*z.size) + + lobe = objects[::2,:] + dx = lobe[:,0] - xc + dy = lobe[:,1] - yc + + r = 2.0*dhPSF.sigma + + # At z_min the lobes lie along x, at z_max along y, and at focus the + # angle is halfway between. + assert(numpy.allclose([dx[0], dy[0]], [r, 0.0], atol = 1.0e-9)) + assert(numpy.allclose([dx[2], dy[2]], [0.0, r], atol = 1.0e-9)) + assert(abs(dx[1] - dy[1]) < 1.0e-9) + + # The angle has to actually sweep, which is the whole point. + angles = numpy.degrees(numpy.arctan2(dy, dx)) + assert(numpy.allclose(angles, [0.0, 45.0, 90.0], atol = 1.0e-6)) + + # The second lobe of each pair is opposite the first. + other = objects[1::2,:] + assert(numpy.allclose(other[:,0] - xc, -dx)) + assert(numpy.allclose(other[:,1] - yc, -dy)) + + if (__name__ == "__main__"): test_psf_spline2D_1() test_psf_spline2D_2() @@ -198,3 +241,4 @@ def test_astigmatic_psf(): test_simulate_2() test_simulate_3() test_astigmatic_psf() + test_dh_psf()