Carry the track length and start frame through to Insight3

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:

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: <bin_file> <smlc_file> <pix_to_nm>")
-    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: <bin_file> <smlc_file> <pix_to_nm>")
+        exit()
 
-print("Saved", index, "molecules.")
-smlc_file_fp.close()
+    binToLMChallenge(sys.argv[1], sys.argv[2], float(sys.argv[3]))
 
 #
 # The MIT License

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("<settings/>")
+        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()