Fix nine latent defects

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

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;i<fit_data->fit_size_x;i++){

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):
         """

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
 

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;i<fit_data->nfit;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;i<fit_data->nfit;i++){
       values[i] = fit_data->fit[i].params[HEIGHT];

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 *);

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

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,23 +16,25 @@
 
 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])
 
 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
 
 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
 

storm_analysis/simulator/dhPSF.py

@@ -4,16 +4,20 @@
 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
 """
 
 import numpy
 
 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):

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

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

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

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