From d4c8c90d1820eef06a01e56588807101565d6287 Mon Sep 17 00:00:00 2001 From: Hazen Babcock Date: Tue, 8 Sep 2026 14:50:00 -0400 Subject: [PATCH] Carry the track length and start frame through to Insight3 bin_to_lmchallenge_format.py expands a molecule that was on for several frames back into one row per frame: track_length = i3_block['tl'][i] for j in range(track_length): fr = i3_block['fr'][i] + j That worked against an Insight3 era file. tracker.c wrote object_data_int[TLEN] = cur->track_length, avemlist.c flagged TLEN as NOAVERAGE so it survived averaging, and insight.h defines FRAME as "frame the peak first occured in". hdf5ToBin(), which replaced that path in Dec 2017, wrote neither. createFromMultiFit() never touched 'tl', so it kept its default of 1, and every track was written into frame 1. The expansion was silently a no-op: in the HDF5 exported track_length / tl 1, 3, 10 1, 1, 1 frame_number / fr 0, 4, 20 1, 1, 1 rows out of lmchallenge 14 3 Map both fields in createFromMultiFit(), which is where every other HDF5 field is translated. Insight3 frame numbers start at 1 and HDF5 frame numbers start at 0, hence the +1. Localizations carry neither field so that path is unchanged. Restoring 'tl' on its own would have made things worse, which is why 2e43489b deliberately left it alone: each expanded row was given the whole of 'a', so a ten frame track would have reported its full photon count ten times. 'a' is a track total in our files, and in Insight3 era ones too since avemlist flagged AREA as TOTAL, so each row now gets its share. bin_to_lmchallenge_format.py was executable-only, all module level code reading sys.argv at import, so nothing could test it. Wrapped in binToLMChallenge() with a __main__ guard, like every other utility here, and it returns the row count. Note that a track with a gap in it is still expanded into consecutive frames. The .bin format records only where a track started and how long it was, so the information needed to do better is not there. Co-Authored-By: Claude Opus 5 --- storm_analysis/sa_library/i3dtype.py | 10 +++ .../sa_utilities/bin_to_lmchallenge_format.py | 79 +++++++++++------- storm_analysis/test/test_hdf5_to_bin.py | 81 +++++++++++++++++++ 3 files changed, 143 insertions(+), 27 deletions(-) diff --git a/storm_analysis/sa_library/i3dtype.py b/storm_analysis/sa_library/i3dtype.py index a7eefdca..c052bb6a 100644 --- a/storm_analysis/sa_library/i3dtype.py +++ b/storm_analysis/sa_library/i3dtype.py @@ -78,6 +78,14 @@ def createFromMultiFit(peaks, frame, nm_per_pixel): setI3Field(i3data, 'bg', peaks["background"]) if "error" in peaks: setI3Field(i3data, 'i', peaks["error"]) + if "frame_number" in peaks: + + # + # Tracks know which frame they started in, so use that rather than the + # 'frame' argument. HDF5 frame numbers start at 0 and Insight3 frame + # numbers start at 1. + # + setI3Field(i3data, 'fr', peaks["frame_number"] + 1) if "height" in peaks: setI3Field(i3data, 'h', peaks["height"]) # if "iterations" in peaks: @@ -88,6 +96,8 @@ def createFromMultiFit(peaks, frame, nm_per_pixel): setI3Field(i3data, 'fi', peaks["status"]) if "sum" in peaks: setI3Field(i3data, 'a', peaks["sum"]) + if "track_length" in peaks: + setI3Field(i3data, 'tl', peaks["track_length"]) if "x" in peaks: posSet(i3data, 'x', peaks["x"] + 1.0) if "xsigma" in peaks: diff --git a/storm_analysis/sa_utilities/bin_to_lmchallenge_format.py b/storm_analysis/sa_utilities/bin_to_lmchallenge_format.py index b4cbdde4..7a1d70e2 100644 --- a/storm_analysis/sa_utilities/bin_to_lmchallenge_format.py +++ b/storm_analysis/sa_utilities/bin_to_lmchallenge_format.py @@ -16,42 +16,67 @@ import storm_analysis.sa_library.readinsight3 as readinsight3 -if (len(sys.argv)!=4): - print("usage: ") - exit() -i3_reader = readinsight3.I3Reader(sys.argv[1]) -i3_block = i3_reader.nextBlock(block_size = 1000, good_only = False) +def binToLMChallenge(bin_name, smlc_name, pix_to_nm, verbose = True): + """ + A molecule that was on for several frames is stored as a single track, + so each one is expanded back into one row per frame. + + 'a' is the total number of photons over the whole track, both in the + files this package writes and in the Insight3 era ones, where avemlist + flagged AREA as a total. So each row gets its share of it rather than + all of it. + + Note that a track with a gap in it is expanded into consecutive frames + anyway, since the .bin format records only where the track started and + how long it was. + """ + i3_reader = readinsight3.I3Reader(bin_name) + i3_block = i3_reader.nextBlock(block_size = 1000, good_only = False) -smlc_file_fp = open(sys.argv[2], "w") -smlc_file_fp.write("index, frame, xnano, ynano, znano, intensity\n") + smlc_file_fp = open(smlc_name, "w") + smlc_file_fp.write("index, frame, xnano, ynano, znano, intensity\n") -pix_to_nm = float(sys.argv[3]) + if verbose: + print("Saving Localizations") + localization_number = 0 + index = 0 + while (type(i3_block) != type(False)): -print("Saving Localizations") -localization_number = 0 -index = 0 -while (type(i3_block) != type(False)): + if verbose: + print(" saving localization", localization_number) - print(" saving localization", localization_number) + for i in range(len(i3_block)): + track_length = int(i3_block['tl'][i]) + if (track_length < 1): + track_length = 1 + for j in range(track_length): + fr = i3_block['fr'][i] + j + xp = i3_block['xc'][i] * pix_to_nm + yp = i3_block['yc'][i] * pix_to_nm + zp = i3_block['zc'][i] + intensity = i3_block['a'][i] / float(track_length) - for i in range(len(i3_block)): - track_length = i3_block['tl'][i] - for j in range(track_length): - fr = i3_block['fr'][i] + j - xp = i3_block['xc'][i] * pix_to_nm - yp = i3_block['yc'][i] * pix_to_nm - zp = i3_block['zc'][i] - intensity = i3_block['a'][i] + index += 1 + smlc_file_fp.write("{0:d}, {1:d}, {2:.3f}, {3:.3f}, {4:.3f}, {5:.3f}\n".format(index, fr, xp, yp, zp, intensity)) - index += 1 - smlc_file_fp.write("{0:d}, {1:d}, {2:.3f}, {3:.3f}, {4:.3f}, {5:.3f}\n".format(index, fr, xp, yp, zp, intensity)) + localization_number += len(i3_block) + i3_block = i3_reader.nextBlock(block_size = 1000, good_only = False) - localization_number += len(i3_block) - i3_block = i3_reader.nextBlock(block_size = 1000, good_only = False) + if verbose: + print("Saved", index, "molecules.") + smlc_file_fp.close() + + return index + + +if (__name__ == "__main__"): + + if (len(sys.argv)!=4): + print("usage: ") + exit() -print("Saved", index, "molecules.") -smlc_file_fp.close() + binToLMChallenge(sys.argv[1], sys.argv[2], float(sys.argv[3])) # # The MIT License diff --git a/storm_analysis/test/test_hdf5_to_bin.py b/storm_analysis/test/test_hdf5_to_bin.py index a290aa2f..44039802 100644 --- a/storm_analysis/test/test_hdf5_to_bin.py +++ b/storm_analysis/test/test_hdf5_to_bin.py @@ -130,8 +130,89 @@ def test_hdf5_to_bin_track_normalization(): assert(numpy.allclose(i3_data['h'], lengths*height)) +def test_hdf5_to_bin_track_expansion(): + """ + Tracks carry their start frame and length through to the Insight3 file, + and bin_to_lmchallenge_format expands them back into one row per frame. + + hdf5ToBin() wrote every track into frame 1 and left 'tl' at its default + of 1, so the expansion was silently a no-op and a movie's worth of + tracks came out as one row each, all in frame 1. + """ + import storm_analysis.sa_utilities.bin_to_lmchallenge_format as binToLMC + + photons = 6000.0 + pixel_size = 100.0 + + lengths = numpy.array([1, 3, 10], dtype = numpy.int32) + starts = numpy.array([0, 4, 20], dtype = numpy.int32) + n = lengths.size + + tracks = {"x" : numpy.arange(n, dtype = numpy.float64) + 10.0, + "y" : numpy.arange(n, dtype = numpy.float64) + 20.0, + "z" : numpy.zeros(n), + "category" : numpy.zeros(n, dtype = numpy.int32), + "frame_number" : starts, + "track_id" : numpy.arange(n, dtype = numpy.int64), + "track_length" : lengths, + "xsigma" : lengths*1.5, + "ysigma" : lengths*1.5, + "height" : lengths*1000.0, + "background" : lengths*20.0, + "sum" : lengths*photons} + + h5_name = storm_analysis.getPathOutputTest("test_h5_to_bin_expand.hdf5") + storm_analysis.removeFile(h5_name) + + with saH5Py.SAH5Py(h5_name, is_existing = False) as h5: + h5.addMetadata("") + h5.setMovieInformation(256, 256, 40, "XYZZY") + h5.setPixelSize(pixel_size) + h5.addLocalizations({"x" : numpy.zeros(1), "y" : numpy.zeros(1)}, 0) + h5.addTracks(tracks) + + i3_name = storm_analysis.getPathOutputTest("test_h5_to_bin_expand.bin") + storm_analysis.removeFile(i3_name) + hdf5ToBin.hdf5ToBin(h5_name, i3_name) + + i3_data = readinsight3.loadI3File(i3_name, verbose = False) + + # The track length survives, and the frame is where the track started. + # Insight3 frame numbers start at 1. + assert(numpy.allclose(i3_data['tl'], lengths)) + assert(numpy.allclose(i3_data['fr'], starts + 1)) + + # 'a' is still the whole track's photons. + assert(numpy.allclose(i3_data['a'], lengths*photons)) + + ## Now expand it. + txt_name = storm_analysis.getPathOutputTest("test_h5_to_bin_expand.txt") + storm_analysis.removeFile(txt_name) + + n_rows = binToLMC.binToLMChallenge(i3_name, txt_name, pixel_size, verbose = False) + + assert(n_rows == numpy.sum(lengths)) + + with open(txt_name) as fp: + rows = [line.split(",") for line in fp.readlines()[1:]] + + assert(len(rows) == numpy.sum(lengths)) + + # Each track occupies consecutive frames from where it started, and each + # row carries its own share of the photons rather than the whole track's. + at = 0 + for i in range(n): + for j in range(lengths[i]): + frame = int(rows[at][1]) + intensity = float(rows[at][5]) + assert(frame == starts[i] + 1 + j) + assert(abs(intensity - photons) < 1.0e-3) + at += 1 + + if (__name__ == "__main__"): test_hdf5_to_bin_1() test_hdf5_to_bin_2() test_hdf5_to_bin_track_normalization() + test_hdf5_to_bin_track_expansion()