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()