Merge pull request #167 from BrillouinMicroscopy/feature/noid/4-peak-fit

Implement 4-peak fitting

bmlab/controllers.py

@@ -703,24 +703,49 @@ def create_bounds(self, brillouin_regions, rayleigh_peaks):
         if bounds is None:
             return None
 
-        # We need a Rayleigh peak position for
-        # every brillouin_region
-        if rayleigh_peaks.shape[0] != (len(brillouin_regions)):
+        # We need two Rayleigh peak positions to determine the FSR
+        # and decide whether a peak is Stokes or Anti-Stokes
+        if rayleigh_peaks.shape[0] != 2:
             return None
 
         w0_bounds = []
         # We have to create a separate bound for every region
         for region_idx, region in enumerate(brillouin_regions):
             local_time = []
-            for f_rayleigh in rayleigh_peaks[region_idx, :]:
-                # In case this is an Anti-Stokes peak, we find the peak
-                # with the higher frequency on the left hand side and
-                # have to flip the bounds
-                is_anti_stokes = np.mean(region) <\
-                                 f_rayleigh
+            for rayleigh_idx in range(rayleigh_peaks.shape[1]):
+                anti_stokes_limit = np.nanmean(rayleigh_peaks[:, rayleigh_idx])
+                # We need to figure out, whether a peak is a
+                # Stokes or Anti-Stokes peak in order to decide
+                # to which Rayleigh peak we need to relate to.
+                tmp = region >= anti_stokes_limit
+                is_pure_region = (tmp == tmp[0]).all()
+
+                is_anti_stokes_region = np.mean(region) >\
+                    anti_stokes_limit
 
                 local_bound = []
                 for bound in bounds:
+                    parsed_bound = []
+                    for limit in bound:
+                        try:
+                            parsed_bound.append(float(limit))
+                        except ValueError:
+                            parsed_bound.append(np.NaN)
+                    # We don't treat Inf as a value
+                    with warnings.catch_warnings():
+                        warnings.filterwarnings(
+                            action='ignore',
+                            message='Mean of empty slice'
+                        )
+                        is_anti_stokes_peak =\
+                            np.nanmean(
+                                np.array(parsed_bound)[
+                                    np.isfinite(parsed_bound)]
+                            ) < 0
+
+                    is_anti_stokes = is_anti_stokes_region if\
+                        is_pure_region else is_anti_stokes_peak
+
                     local_limit = []
                     for limit in bound:
                         if limit.lower() == 'min':
@@ -739,10 +764,11 @@ def create_bounds(self, brillouin_regions, rayleigh_peaks):
                             # a value in pixel depending on the time
                             try:
                                 val = ((-1) ** is_anti_stokes)\
-                                          * 1e9 * float(limit) + f_rayleigh
+                                          * 1e9 * abs(float(limit))\
+                                          + rayleigh_peaks[
+                                          int(is_anti_stokes), rayleigh_idx]
                             except BaseException:
                                 val = np.Inf
-                        # print(val)
                         local_limit.append(val)
                     # Check that the bounds are sorted ascendingly
                     # (for anti-stokes, they might not).
@@ -890,55 +916,43 @@ def get_indices_from_key(resolution, key):
 
 def calculate_derived_values():
     """
-    We calculate the derived parameters here:
-    - Brillouin shift [GHz]
+    We calculate the Brillouin shift in GHz here
     """
     session = Session.get_instance()
     evm = session.evaluation_model()
     if not evm:
         return
 
-    if len(evm.results['time']) == 0:
+    if evm.results['brillouin_peak_position_f'].size == 0:
+        return
+
+    if evm.results['rayleigh_peak_position_f'].size == 0:
         return
 
     shape_brillouin = evm.results['brillouin_peak_position_f'].shape
     shape_rayleigh = evm.results['rayleigh_peak_position_f'].shape
-    # If we have the same number of Rayleigh and Brillouin regions,
-    # we can simply subtract the two arrays (regions are always
-    # sorted by center in the peak selection model, so corresponding
-    # regions should be at the same array index)
-    if shape_brillouin[4] == shape_rayleigh[4]:
-        evm.results['brillouin_shift_f'] = abs(
+
+    # We calculate every possible combination of
+    # Brillouin peak and Rayleigh peak position difference
+    # and then use the smallest absolute value.
+    # That saves us from sorting Rayleigh peaks to Brillouin peaks,
+    # because a Brillouin peak always belongs to the Rayleigh peak nearest.
+    brillouin_shift_f = np.NaN * np.ones((*shape_brillouin, shape_rayleigh[4]))
+    for idx in range(shape_rayleigh[4]):
+        brillouin_shift_f[:, :, :, :, :, :, idx] = abs(
             evm.results['brillouin_peak_position_f'] -
-            evm.results['rayleigh_peak_position_f']
+            np.tile(
+                evm.results['rayleigh_peak_position_f'][:, :, :, :, [idx], :],
+                (1, 1, 1, 1, 1, shape_brillouin[5])
+            )
         )
 
-    # Having a different number of Rayleigh and Brillouin regions
-    # doesn't really make sense. But in case I am missing something
-    # here, we assign each Brillouin region the nearest (by center)
-    # Rayleigh region.
-    else:
-        psm = session.peak_selection_model()
-        if not psm:
-            return
-        brillouin_centers = list(map(np.mean, psm.get_brillouin_regions()))
-        rayleigh_centers = list(map(np.mean, psm.get_rayleigh_regions()))
-
-        for idx in range(len(brillouin_centers)):
-            # Find corresponding (nearest) Rayleigh region
-            d = list(
-                map(
-                    lambda x: abs(x - brillouin_centers[idx]),
-                    rayleigh_centers
-                )
-            )
-            idx_r = d.index(min(d))
-            evm.results['brillouin_shift_f'][:, :, :, :, idx, :] = abs(
-                evm.results[
-                    'brillouin_peak_position_f'][:, :, :, :, idx, :] -
-                evm.results[
-                    'rayleigh_peak_position_f'][:, :, :, :, idx_r, :]
-            )
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            action='ignore',
+            message='All-NaN slice encountered'
+        )
+        evm.results['brillouin_shift_f'] = np.nanmin(brillouin_shift_f, 6)
 
 
 class Controller(object):

bmlab/fits.py

@@ -141,6 +141,122 @@ def error(params, xdata, ydata):
     return w0s, fwhms, intens, offset
 
 
+def fit_quadruple_lorentz(x, y, bounds_w0=None):
+    offset_guess = (y[0] + y[-1]) / 2.
+    # gam_guess = (w0_guess ** 2 * (np.max(y) - offset_guess)) ** -0.5
+    fwhm_guess = 10 * (x[-1] - x[0]) / x.shape[0]
+    intensity_guess = np.max(y) - offset_guess
+
+    # We run peak finding to get a good guess of the peak positions
+    # and use the four peaks with the highest prominence.
+    peaks, properties = find_peaks(y, prominence=1)
+    idx = np.argsort(properties['prominences'])[::-1]
+
+    # Return if we didn't find four peaks
+    if len(idx) < 4:
+        return
+
+    idx_sort = np.sort(peaks[idx[0:4]])
+    w0_guess = list(x[idx_sort])
+
+    def error(params, xdata, ydata):
+        return (ydata
+                - lorentz(xdata, *params[0:3])
+                - lorentz(xdata, *params[3:6])
+                - lorentz(xdata, *params[6:9])
+                - lorentz(xdata, *params[9:12])
+                - params[12]) ** 2
+
+    # Create the bounds array
+    if bounds_w0 is not None:
+        # Lower limits
+        bounds_lower = -np.Inf * np.ones(13)
+
+        # 1st peak:
+        # central position
+        bounds_lower[0] = bounds_w0[0][0]
+        # full-width-half-maximum
+        # The VIPA spectrometer has an instrument width of
+        # approx. 750 MHz for the FOB setup
+        # and 180 MHz for the 780 nm setup.
+        # This is far higher than the step size,
+        # so we limit it to this.
+        fwhm_lower_bound = (x[-1] - x[0]) / x.shape[0]
+
+        bounds_lower[1] = fwhm_lower_bound
+        # intensity
+        bounds_lower[2] = 0
+
+        # 2nd peak:
+        # central position
+        bounds_lower[3] = bounds_w0[1][0]
+        # full-width-half-maximum
+        bounds_lower[4] = fwhm_lower_bound
+        # intensity
+        bounds_lower[5] = 0
+
+        # 3rd peak:
+        # central position
+        bounds_lower[6] = bounds_w0[2][0]
+        # full-width-half-maximum
+        bounds_lower[7] = fwhm_lower_bound
+        # intensity
+        bounds_lower[8] = 0
+
+        # 4th peak:
+        # central position
+        bounds_lower[9] = bounds_w0[3][0]
+        # full-width-half-maximum
+        bounds_lower[10] = fwhm_lower_bound
+        # intensity
+        bounds_lower[11] = 0
+
+        # offset
+        bounds_lower[12] = 0
+
+        # Upper limits
+        bounds_upper = np.Inf * np.ones(13)
+
+        bounds_upper[0] = bounds_w0[0][1]
+        bounds_upper[3] = bounds_w0[1][1]
+        bounds_upper[6] = bounds_w0[2][1]
+        bounds_upper[9] = bounds_w0[3][1]
+        bounds = (bounds_lower, bounds_upper)
+
+        # Sort the guesses to the bounds
+        bounds_w0_center = [np.mean(
+            np.clip(bound, *x[::len(x) - 1])) for bound in bounds_w0]
+        w0_guess.sort(reverse=(bounds_w0_center[0] > bounds_w0_center[1]))
+
+        # Check that the initial guesses are within the bounds
+        w0_guess = [np.clip(
+            guess, *bounds_w0[idx]) for idx, guess in enumerate(w0_guess)]
+    else:
+        bounds = (-np.Inf, np.Inf)
+
+    opt_result = least_squares(
+        error,
+        x0=(w0_guess[0], fwhm_guess, intensity_guess,
+            w0_guess[1], fwhm_guess, intensity_guess,
+            w0_guess[2], fwhm_guess, intensity_guess,
+            w0_guess[3], fwhm_guess, intensity_guess,
+            offset_guess
+            ),
+        args=(x, y),
+        bounds=bounds
+    )
+
+    if not opt_result.success:
+        raise FitError('Lorentz fit failed.')
+
+    res = opt_result.x
+    w0s, fwhms, intens = (res[0], res[3], res[6], res[9]),\
+                         (res[1], res[4], res[7], res[10]),\
+                         (res[2], res[5], res[8], res[11])
+    offset = res[12]
+    return w0s, fwhms, intens, offset
+
+
 def fit_circle(points):
     """
     Fits a circle to a given set of points. Returnes the center and the radius
@@ -224,13 +340,17 @@ def fit_lorentz_region(region, xdata, ydata, nr_peaks=1, bounds_w0=None):
         idx_r = np.nanargmin(np.abs(xdata - region[1]))
         x = xdata[idx_l:idx_r]
         y = ydata[idx_l:idx_r]
-        if nr_peaks == 2:
+        if nr_peaks == 1:
+            w0s, fwhms, intensities, offset = fit_lorentz(
+                x, y)
+        elif nr_peaks == 2:
             w0s, fwhms, intensities, offset = fit_double_lorentz(
                 x, y,
                 bounds_w0=bounds_w0)
-        elif nr_peaks == 1:
-            w0s, fwhms, intensities, offset = fit_lorentz(
-                x, y)
+        elif nr_peaks == 4:
+            w0s, fwhms, intensities, offset = fit_quadruple_lorentz(
+                x, y,
+                bounds_w0=bounds_w0)
         else:
             return
     except Exception: