python-microscopy · barentine · Mar 7, 2021 · Mar 7, 2021 · Mar 8, 2021 · Mar 8, 2021
diff --git a/PYME/Acquire/Scripts/init_htsms.py b/PYME/Acquire/Scripts/init_htsms.py
@@ -256,6 +256,23 @@ def chained_analysis(main_frame, scope):
     defaults['htsms-flow'] = RuleChain([get_rule_tile(SpoolLocalLocalizationRuleFactory)(analysisMetadata=mdh)])
     defaults['htsms-staggered'] = RuleChain([get_rule_tile(SpoolLocalLocalizationRuleFactory)(analysisMetadata=mdh)])
 
+    # calibrations
+    mdh = DictMDHandler(mdh)
+    mdh['Analysis.subtractBackground'] = False
+    mdh['Analysis.GPUPCTBackground'] = False
+    mdh['Analysis.BGRange'] = [-0, 0]
+    mdh['Analysis.PCTBackground'] = 0
+    mdh['Analysis.StartAt'] = 0
+    reg_loc_tile = get_rule_tile(SpoolLocalLocalizationRuleFactory)(analysisMetadata=mdh)
+
+    reg_rec = os.path.join(rec_dir, 'shiftmaps-from-tiled-beads.yaml')
+    with open(reg_rec) as f:
+        reg_rec = f.read()
+    reg_rec_tile = get_rule_tile(RecipeRuleFactory)(recipe=reg_rec)
+    defaults['htsms-cal-registration'] = RuleChain([reg_loc_tile,
+                                                    reg_rec_tile])
+
+
     SMLMChainedAnalysisPanel.plug(main_frame, scope, defaults)
 
 @init_gui('Tiling')

diff --git a/PYME/Analysis/points/multiview.py b/PYME/Analysis/points/multiview.py
@@ -308,6 +308,114 @@ def pair_molecules(t_index, x0, y0, which_chan, delta_x=[None], appear_in=np.ara
     else:
         return assigned
 
+def raw_shifts_by_frame(x0, y0, t0, x0err, y0err, x1, y1, t1, x1err, y1err,
+                        clump_distance=1000):
+    """Calculate raw shifts between points and uncertainty appropriate for
+    creating a shiftmap between two lists of points. Points are grouped only
+    if they share the same `t` and are within `clump_distance`
+
+    Parameters
+    ----------
+    x0 : ndarray
+        x positions in channel 0
+    y0 : ndarray
+        y positions in channel 0
+    t0 : ndarray
+        frame index for points in channel 0
+    x0err : ndarray
+        uncertainties for x positions in channel 0
+    y0err : ndarray
+        uncertainties for y positions in channel 0
+    x1 : ndarray
+        x positions in channel 1
+    y1 : ndarray
+        y positions in channel 1
+    t1 : ndarray
+        frame index for points in channel 1
+    x1err : ndarray
+        uncertainties for x positions in channel 1
+    y1err : ndarray
+        uncertainties for y positions in channel 1
+    clump_distance : int, ndarray, optional
+        distance to search, by default 1000
+
+    Returns
+    -------
+    xc : ndarray
+        x position (in channel 0) of grouped pairs
+    yc : ndarray
+        y position (in channel 0) of grouped pairs
+    dx : ndarray
+        shift in x, per pair of grouped points
+    dy : ndarray
+        shift in y, per pair of grouped points
+    dxerr : ndarray
+        uncertainty in `dx`
+    dyerr : ndarray
+        uncertainty in `dy`
+    """
+    from PYME.Analysis.points.DeClump import findClumps
+    # take out any large linear shifts for the sake of easier pairing
+    # x, y = correlative_shift(x0, y0, which_chan, clump_distance / 10)
+
+    x = np.concatenate([x0, x1])
+    y = np.concatenate([y0, y1])
+    t = np.concatenate([t0, t1])
+    xerr = np.concatenate([x0err, x1err])
+    yerr = np.concatenate([y0err, y1err])
+    channel = np.concatenate([np.zeros(len(t0)), np.ones(len(t1))])
+
+    # sort by time for pyDeclump.findClumps
+    I = np.argsort(t)
+    x = x[I]
+    y = y[I]
+    t = t[I]
+    xerr = xerr[I]
+    yerr = yerr[I]
+    channel = channel[I]
+
+    if np.isscalar(clump_distance):
+        clump_distance = clump_distance*np.ones_like(x)
+
+    assigned = findClumps(t.astype(np.int32), x, y, clump_distance, -1, True)
+    ids, count = np.unique(assigned, return_counts=True)
+
+    # reorder by id (is this already done actually?)
+    I = np.argsort(ids)
+    x = x[I]
+    y = y[I]
+    t = t[I]
+    xerr = xerr[I]
+    yerr = yerr[I]
+    channel = channel[I]
+
+    # iterate through and calculate shifts for clumps with both channels
+    xc, yc, dx, dy, dxerr, dyerr = [], [], [], [], [], []
+    start = 0
+    for ind, clump in enumerate(ids):
+        end = start + count[ind]
+        chan = channel[start:end]
+        # only look at pairings with exactly one localization from each channel
+        if (len(chan) == 2) and (0 in chan) and (1 in chan):
+            # grab the 0 and 1 channel indices for the full arrays
+            chan0 = start + np.argmin(chan)
+            chan1 = start + np.argmax(chan)
+
+            xc.append(x[chan0])
+            yc.append(y[chan0])
+
+            # weirdly looks like we're supposed to do 0 - 1 instead of 1 - 0?
+            dx.append(x[chan0] - x[chan1])
+            dy.append(y[chan0] - y[chan1])
+            # add error in quadrature
+            dxerr.append(np.sqrt(xerr[chan0] ** 2 + xerr[chan1] ** 2))
+            dyerr.append(np.sqrt(yerr[chan0] ** 2 + yerr[chan1] ** 2))
+
+        start = end
+
+    return np.array(xc), np.array(yc), np.array(dx), np.array(dy), \
+           np.array(dxerr), np.array(dyerr)
+
 def calc_shifts_for_points(datasource, shiftWallet):
     import importlib
     model = shiftWallet['shiftModel'].split('.')[-1]

diff --git a/PYME/recipes/multiview.py b/PYME/recipes/multiview.py
@@ -1,5 +1,5 @@
 from .base import register_module, ModuleBase, Filter
-from .traits import Input, Output, Float, CStr, Bool, Int, FileOrURI
+from .traits import Input, Output, Float, CStr, Bool, Int, FileOrURI, Enum
 import numpy as np
 from PYME.IO import tabular
 from PYME.IO import MetaDataHandler
@@ -108,7 +108,7 @@ class FindClumps(ModuleBase):
         localizations as PYME.IO.Tabular types
     time_gap_tolerance : traits.Int
         Number of frames which a localizations is allowed to be missing and still be considered the same molecule if it
-        reappears
+        reappears. Currently this is ~off by 2, so time_gap_tolerance=1 means link only on the same frame
     radius_scale : traits.Float
         Factor by which the localization precision is multiplied to determine the search radius for clustering. The
         default of 2 sigma means that we link ~95% of the points which should be linked (if Gaussian statistics hold)
@@ -271,6 +271,13 @@ class CalibrateShifts(ModuleBase):
         localizations as PYME.IO.Tabular types
     search_radius_nm : traits.Float
         radius within which bead localizations should be clumped if the bead appears in all channels. Units of nm.
+    mode : str
+        `average first` averages all points found within `search_radius_nm` by
+        channel before generating shifts between clumps containing all
+        multiview channels. `per-frame` pairs localizations between points on
+        the 0th channel and other channels on the same frame only, and does not
+        require all channels to be present to generate a valid shift at a given
+        location.
 
     Returns
     -------
@@ -284,97 +291,114 @@ class CalibrateShifts(ModuleBase):
     input_name = Input('folded')
 
     search_radius_nm = Float(250.)
+    mode = Enum(['average first', 'per-frame'])
 
     output_name = Output('shiftmap')
 
     def execute(self, namespace):
         from PYME.Analysis.points import twoColour
         from PYME.Analysis.points import multiview
-        from PYME.IO.MetaDataHandler import NestedClassMDHandler
 
         inp = namespace[self.input_name]
 
         try:  # make sure we're looking at multiview data
-            n_chan = inp.mdh['Multiview.NumROIs']
+            n_chan = inp.mdh['Multiview.NumROIs']  # TODO- use Multiview.ActiveViews + index appropriately
         except AttributeError:
             raise AttributeError('multiview metadata is missing or incomplete')
 
-        # sort in frame order
-        I = inp['tIndex'].argsort()
-        x_sort, y_sort = inp['x'][I], inp['y'][I]
-        chan_sort = inp['multiviewChannel'][I]
-
-        clump_id, keep = multiview.pair_molecules(inp['tIndex'][I], x_sort, y_sort, chan_sort,
-                                                  self.search_radius_nm * np.ones_like(x_sort),
-                                                  appear_in=np.arange(n_chan), n_frame_sep=inp['tIndex'].max(),
-                                                  pix_size_nm=inp.mdh.voxelsize_nm.x)
-
-        # only look at the clumps which showed up in all channels
-        x = x_sort[keep]
-        y = y_sort[keep]
-        chan = chan_sort[keep]
-        clump_id = clump_id[keep]
-
-        # Generate raw shift vectors (map of displacements between channels) for each channel
-        mol_list = np.unique(clump_id)
-        n_mols = len(mol_list)
-        if n_mols < 3:
-            raise ValueError('Need at 3 clusters containing points from each channel - try increasing search radius')
-
-        dx = np.zeros((n_chan - 1, n_mols))
-        dy = np.zeros_like(dx)
-        dx_err = np.zeros_like(dx)
-        dy_err = np.zeros_like(dx)
-        x_clump, y_clump, x_std, y_std, x_shifted, y_shifted = [], [], [], [], [], []
-
         shift_map_dtype = [('mx', '<f4'), ('mx2', '<f4'), ('mx3', '<f4'),  # x terms
                            ('my', '<f4'), ('my2', '<f4'), ('my3', '<f4'),  # y terms
                            ('mxy', '<f4'), ('mx2y', '<f4'), ('mxy2', '<f4'),  # cross terms
                            ('x0', '<f4')]  # 0th order shift
 
         shift_maps = np.zeros(2*(n_chan - 1), dtype=shift_map_dtype)
-        mdh = NestedClassMDHandler(inp.mdh)
+        mdh = MetaDataHandler.DictMDHandler(inp.mdh)
         mdh['Multiview.shift_map.legend'] = {}
-
-        for ii in range(n_chan):
-            chan_mask = (chan == ii)
-            x_chan = np.zeros(n_mols)
-            y_chan = np.zeros(n_mols)
-            x_chan_std = np.zeros(n_mols)
-            y_chan_std = np.zeros(n_mols)
-
-            for ind in range(n_mols):
-                # merge clumps within channels
-                clump_mask = np.where(np.logical_and(chan_mask, clump_id == mol_list[ind]))
-                x_chan[ind] = x[clump_mask].mean()
-                y_chan[ind] = y[clump_mask].mean()
-                x_chan_std[ind] = x[clump_mask].std()
-                y_chan_std[ind] = y[clump_mask].std()
-
-            x_clump.append(x_chan)
-            y_clump.append(y_chan)
-            x_std.append(x_chan_std)
-            y_std.append(y_chan_std)
-
-            if ii > 0:
-                dx[ii - 1, :] = x_clump[0] - x_clump[ii]
-                dy[ii - 1, :] = y_clump[0] - y_clump[ii]
-                dx_err[ii - 1, :] = np.sqrt(x_std[ii] ** 2 + x_std[0] ** 2)
-                dy_err[ii - 1, :] = np.sqrt(y_std[ii] ** 2 + y_std[0] ** 2)
-                # generate shiftmap between ii-th channel and the 0th channel
-                dxx, dyy, spx, spy, good = twoColour.genShiftVectorFieldQ(x_clump[0], y_clump[0], dx[ii - 1, :],
-                                                                          dy[ii - 1, :], dx_err[ii - 1, :],
-                                                                          dy_err[ii - 1, :])
-                # store shiftmaps in structured array
-                mdh['Multiview.shift_map.legend']['Chan0%s.X' % ii] = 2*(ii - 1)
-                mdh['Multiview.shift_map.legend']['Chan0%s.Y' % ii] = 2*(ii - 1) + 1
+
+        if self.mode == 'average first':
+            # sort in frame order
+            I = inp['tIndex'].argsort()
+            x_sort, y_sort = inp['x'][I], inp['y'][I]
+            chan_sort = inp['multiviewChannel'][I]
+
+            clump_id, keep = multiview.pair_molecules(inp['tIndex'][I], x_sort, y_sort, chan_sort,
+                                                    self.search_radius_nm * np.ones_like(x_sort),
+                                                    appear_in=np.arange(n_chan), n_frame_sep=inp['tIndex'].max(),
+                                                    pix_size_nm=inp.mdh.voxelsize_nm.x)
+
+            # only look at the clumps which showed up in all channels
+            x = x_sort[keep]
+            y = y_sort[keep]
+            chan = chan_sort[keep]
+            clump_id = clump_id[keep]
+
+            # Generate raw shift vectors (map of displacements between channels) for each channel
+            mol_list = np.unique(clump_id)
+            n_mols = len(mol_list)
+            if n_mols < 3:
+                raise ValueError('Need at 3 clusters containing points from each channel - try increasing search radius')
+
+            dx = np.zeros((n_chan - 1, n_mols))
+            dy = np.zeros_like(dx)
+            dx_err = np.zeros_like(dx)
+            dy_err = np.zeros_like(dx)
+            x_clump, y_clump, x_std, y_std, x_shifted, y_shifted = [], [], [], [], [], []
+
+            for ii in range(n_chan):
+                chan_mask = (chan == ii)
+                x_chan = np.zeros(n_mols)
+                y_chan = np.zeros(n_mols)
+                x_chan_std = np.zeros(n_mols)
+                y_chan_std = np.zeros(n_mols)
+
+                for ind in range(n_mols):
+                    # merge clumps within channels
+                    clump_mask = np.where(np.logical_and(chan_mask, clump_id == mol_list[ind]))
+                    x_chan[ind] = x[clump_mask].mean()
+                    y_chan[ind] = y[clump_mask].mean()
+                    x_chan_std[ind] = x[clump_mask].std()
+                    y_chan_std[ind] = y[clump_mask].std()
+
+                x_clump.append(x_chan)
+                y_clump.append(y_chan)
+                x_std.append(x_chan_std)
+                y_std.append(y_chan_std)
+
+                if ii > 0:
+                    dx[ii - 1, :] = x_clump[0] - x_clump[ii]
+                    dy[ii - 1, :] = y_clump[0] - y_clump[ii]
+                    dx_err[ii - 1, :] = np.sqrt(x_std[ii] ** 2 + x_std[0] ** 2)
+                    dy_err[ii - 1, :] = np.sqrt(y_std[ii] ** 2 + y_std[0] ** 2)
+                    # generate shiftmap between ii-th channel and the 0th channel
+                    dxx, dyy, spx, spy, good = twoColour.genShiftVectorFieldQ(x_clump[0], y_clump[0], dx[ii - 1, :],
+                                                                            dy[ii - 1, :], dx_err[ii - 1, :],
+                                                                            dy_err[ii - 1, :])
+                    # store shiftmaps in structured array
+                    mdh['Multiview.shift_map.legend']['Chan0%s.X' % ii] = 2*(ii - 1)
+                    mdh['Multiview.shift_map.legend']['Chan0%s.Y' % ii] = 2*(ii - 1) + 1
+                    for ki in range(len(shift_map_dtype)):
+                        k = shift_map_dtype[ki][0]
+                        shift_maps[2*(ii - 1)][k] = spx.__getattribute__(k)
+                        shift_maps[2*(ii - 1) + 1][k] = spy.__getattribute__(k)
+        else:
+            # calculate shifts from each frame
+            from PYME.IO.tabular import SelectionFilter
+            chan0 = SelectionFilter(inp, inp['multiviewChannel'] == 0)
+            for ind in range(1, n_chan):
+                other = SelectionFilter(inp, inp['multiviewChannel'] == ind)
+                xc, yc, dx, dy, dxerr, dyerr = multiview.raw_shifts_by_frame(chan0['x'], chan0['y'], chan0['tIndex'], chan0['error_x'], chan0['error_y'],
+                                                                                other['x'], other['y'], other['tIndex'], other['error_x'], other['error_y'],
+                                                                                self.search_radius_nm)
+                dxx, dyy, spx, spy, good = twoColour.genShiftVectorFieldQ(xc, yc, dx, dy, dxerr, dyerr)
+
+                # add into structured array, marking index in metadata
+                mdh['Multiview.shift_map.legend']['Chan0%s.X' % ind] = 2*(ind - 1)
+                mdh['Multiview.shift_map.legend']['Chan0%s.Y' % ind] = 2*(ind - 1) + 1
                 for ki in range(len(shift_map_dtype)):
                     k = shift_map_dtype[ki][0]
-                    shift_maps[2*(ii - 1)][k] = spx.__getattribute__(k)
-                    shift_maps[2*(ii - 1) + 1][k] = spy.__getattribute__(k)
-
-
-                # shift_maps['Chan0%s.X' % ii], shift_maps['Chan0%s.Y' % ii] = spx.__dict__, spy.__dict__
+                    shift_maps[2*(ind - 1)][k] = spx.__getattribute__(k)
+                    shift_maps[2*(ind - 1) + 1][k] = spy.__getattribute__(k)
+
 
         mdh['Multiview.shift_map.model'] = '.'.join([spx.__class__.__module__, spx.__class__.__name__])