feat: new POC method with polynomial fit and more verbose POC methods…

… with ret_details

CHANGELOG

@@ -1,3 +1,9 @@
+3.1.0
+ - feat: new contact point estimation method "fit_constant_polynomial"
+   which applies a piece-wise constant and polynomial fit
+ - feat: contact point estimation methods now return detailed
+   information about the procedure (currently plottable data to
+   understand the process)
 3.0.0
  - BREAKING CHANGE: The contact point estimation method "scheme_2020"
    has been removed, although it has been the default for some time.

docs/nanite.bib

@@ -73,4 +73,16 @@ @Article{Mueller19nanite
   publisher = {Springer Science and Business Media {LLC}},
 }
 
+@Article{Rusaczonek19,
+  author    = {M. Rusaczonek and B. Zapotoczny and M. Szymonski and J. Konior},
+  title     = {Application of a layered model for determination of the elasticity of biological systems},
+  journal   = {Micron},
+  year      = {2019},
+  volume    = {124},
+  pages     = {102705},
+  month     = {sep},
+  doi       = {10.1016/j.micron.2019.102705},
+  publisher = {Elsevier {BV}},
+}
+
 @Comment{jabref-meta: databaseType:bibtex;}

nanite/indent.py

@@ -25,6 +25,8 @@ def __init__(self, data, metadata, diskcache=None):
         self.preprocessing = []
         #: Preprocessing options
         self.preprocessing_options = {}
+        # preprocessing details (for user convenience)
+        self._preprocessing_details = {}
         # protected fit properties
         self._fit_properties = FitProperties()
 
@@ -53,7 +55,8 @@ def fit_properties(self):
     def fit_properties(self, fp):
         self._fit_properties.update(fp)
 
-    def apply_preprocessing(self, preprocessing=None, options=None):
+    def apply_preprocessing(self, preprocessing=None, options=None,
+                            ret_details=False):
         """Perform curve preprocessing steps
 
         Parameters
@@ -66,6 +69,8 @@ def apply_preprocessing(self, preprocessing=None, options=None):
         options: dict of dict
             Dictionary of keyword arguments for each preprocessing
             step (if applicable)
+        ret_details:
+            Return preprocessing details dictionary
         """
         if preprocessing is None:
             preprocessing = self.preprocessing
@@ -80,28 +85,34 @@ def apply_preprocessing(self, preprocessing=None, options=None):
         else:
             preproc_past = []
 
-        if preproc_past != [preprocessing, options]:
+        if ((preproc_past != [preprocessing, options])
+                or (not self._preprocessing_details and ret_details)):
             # Remember initial fit parameters for user convenience
             fp = self.fit_properties
             fp["preprocessing"] = preprocessing
             fp["preprocessing_options"] = options
             # Reset all data
             self.reset()
             # Apply preprocessing
-            IndentationPreprocessor.apply(apret=self,
-                                          identifiers=preprocessing,
-                                          options=options)
+            details = IndentationPreprocessor.apply(apret=self,
+                                                    identifiers=preprocessing,
+                                                    options=options,
+                                                    ret_details=ret_details)
+            self._preprocessing_details = details
             # Check availability of axes
             for ax in ["x_axis", "y_axis"]:
                 # make sure the fitting axes are defined
                 if ax in fp and not fp[ax] in self:
                     fp.pop(ax)
             # Set new fit properties
             self.fit_properties = fp
+
         # remember preprocessing
         self.preprocessing = preprocessing
         self.preprocessing_options = copy.deepcopy(options)
 
+        return self._preprocessing_details
+
     def compute_emodulus_mindelta(self, callback=None):
         """Elastic modulus in dependency of maximum indentation
 

nanite/poc.py

@@ -22,24 +22,44 @@ def compute_preproc_clip_approach(force):
     return fg
 
 
-def compute_poc(force, method="deviation_from_baseline"):
+def compute_poc(force, method="deviation_from_baseline", ret_details=False):
     """Compute the contact point from force data
 
+    Parameters
+    ----------
+    force: 1d ndarray
+        Force data
+    method: str
+        Name of the method for computing the POC (see :const:`POC_METHODS`)
+    ret_details: bool
+        Whether or not to return a dictionary with details alongside the
+        POC estimate.
+
+    Notes
+    -----
     If the POC method returns np.nan, then the center of the
-    force data is used.
+    force data is returned (to allow fitting algorithms to proceed).
     """
     # compute POC according to method chosen
     for mfunc in POC_METHODS:
         if mfunc.identifier == method:
             if "clip_approach" in mfunc.preprocessing:
                 force = compute_preproc_clip_approach(force)
-            cp = mfunc(force)
+            data = mfunc(force, ret_details=ret_details)
+            if ret_details:
+                cp, details = data
+                details["method"] = method
+            else:
+                cp, details = data, None
             break
     else:
         raise ValueError(f"Undefined POC method '{method}'!")
     if np.isnan(cp):
         cp = force.size // 2
-    return cp
+    if ret_details:
+        return cp, details
+    else:
+        return cp
 
 
 def poc(identifier, name, preprocessing):
@@ -78,7 +98,7 @@ def attribute_setter(func):
 @poc(identifier="deviation_from_baseline",
      name="Deviation from baseline",
      preprocessing=["clip_approach"])
-def poc_deviation_from_baseline(force):
+def poc_deviation_from_baseline(force, ret_details=False):
     """Deviation from baseline
 
     1. Obtain the baseline (initial 10% of the gradient curve)
@@ -87,6 +107,7 @@ def poc_deviation_from_baseline(force):
        twice of the maximum deviation
     """
     cp = np.nan
+    details = {}
     # Crop the slow approach trace (10% of the curve)
     baseline = force[:int(force.size * .1)]
     if baseline.size:
@@ -97,33 +118,55 @@ def poc_deviation_from_baseline(force):
         maxid = np.argmax(bl_dev)
         if bl_dev[maxid]:
             cp = maxid
-    return cp
+        if ret_details:
+            x = [0, force.size-1]
+            details["plot force"] = [np.arange(force.size), force]
+            details["plot baseline mean"] = [x, [bl_avg, bl_avg]]
+            details["plot baseline threshold"] = [x, [bl_avg + bl_rng,
+                                                      bl_avg + bl_rng]]
+            details["plot poc"] = [[cp, cp],
+                                   [force.min(), force.max()]]
+
+    if ret_details:
+        return cp, details
+    else:
+        return cp
 
 
 @poc(identifier="fit_constant_line",
      name="Piecewise fit with constant and line",
      preprocessing=["clip_approach"])
-def poc_fit_constant_line(force):
-    """Piecewise fit with constant and line
+def poc_fit_constant_line(force, ret_details=False):
+    r"""Piecewise fit with constant and line
 
     Fit a piecewise function (constant+linear) to the baseline
-    and indentation part.
+    and indentation part:
+
+    .. math::
+
+       F = \text{max}(m\delta, a)
 
     The point of contact is the intersection of a horizontal line
-    (constant) for the baseline and a linear function (constant slope)
+    at :math:`a` (baseline) and a linear function with slope :math:`m`
     for the indentation part.
+
+    The point of contact is defined as :math:`\delta=0` (It's another
+    fitting parameter).
     """
-    def residual(params, x, data):
+    def model(params, x):
         off = params["off"]
         x0 = params["x0"]
         m = params["m"]
         one = off
         two = m * (x - x0) + off
-        curve = np.maximum(one, two)
-        return data - curve
+        return np.maximum(one, two)
+
+    def residual(params, x, data):
+        return data - model(params, x)
 
     y = np.array(force, copy=True)
     cp = np.nan
+    details = {}
     if y.size > 4:  # 3 fit parameters
         ymin, ymax = np.min(y), np.max(y)
         y = (y - ymin) / (ymax - ymin)
@@ -139,24 +182,125 @@ def residual(params, x, data):
         out = lmfit.minimize(residual, params, args=(x, y))
         if out.success:
             cp = int(out.params["x0"])
-    return cp
+            if ret_details:
+                details["plot force"] = [x, y]
+                details["plot fit"] = [np.arange(force.size),
+                                       model(out.params, x)]
+                details["plot poc"] = [[cp, cp],
+                                       [y.min(), y.max()]]
+
+    if ret_details:
+        return cp, details
+    else:
+        return cp
+
+
+@poc(identifier="fit_constant_polynomial",
+     name="Piecewise fit with constant and polynomial",
+     preprocessing=["clip_approach"])
+def poc_fit_constant_polynomial(force, ret_details=False):
+    r"""Piecewise fit with constant and line
+
+    Fit a piecewise function (constant + polynomial) to the baseline
+    and indentation part.
+
+    .. math::
+
+       F = \frac{\delta^3}{a\delta^2 + b\delta + c} + d
+
+    This function is defined for all :math:`\delta>0`. For all
+    :math:`\delta<0` the model evaluates to :math:`d` (baseline).
+
+    I'm not sure where this has been described initially, but it is
+    used e.g. in :cite:`Rusaczonek19`.
+
+    For small indentations, this function exhibits a cubic behavior:
+
+    .. math::
+
+       y \approx \delta^3/c
+
+    And for large indentations, this function is linear:
+
+    .. math::
+
+       y \approx \delta/a - b/a^2
+
+    The point of contact is defined as :math:`\delta=0` (It's another
+    fitting parameter).
+    """
+    def model(params, x):
+        off = params["off"].value
+        x0 = params["x0"].value
+        a = params["a"].value
+        b = params["b"].value
+        c = params["c"].value
+        x1 = x - x0
+        curve = x1**3 / (a*x1**2 + b*x1 + c) + off
+        curve[x1 <= 0] = off
+        return curve
+
+    def residual(params, x, data):
+        curve = model(params, x)
+        return data - curve
+
+    y = np.array(force, copy=True)
+    cp = np.nan
+    details = {}
+    if y.size > 4:  # 3 fit parameters
+        ymin, ymax = np.min(y), np.max(y)
+        y = (y - ymin) / (ymax - ymin)
+        x = np.arange(y.size)
+        x0 = poc_deviation_from_baseline(force)
+        if np.isnan(x0):
+            x0 = y.size // 2
+        params = lmfit.Parameters()
+        params.add('off', value=np.mean(y[:10]))
+        params.add('x0', value=x0)
+        # The polynomial fitting parameters are supposed to be
+        # greater than zero (source?). We set the minimum to 1e-3 so
+        # the fitting algorithm becomes more stable. Also, the initial
+        # values for b and c are more or less arbitrary (this is a heuristic
+        # approach).
+        # for larger x, a is something like an inverse slope. Since we
+        # normalized the y-values to 1, we just take the x-difference.
+        params.add('a', value=(y.size-x0), min=1e-3, max=100*(y.size-x0))
+        params.add('b', value=y.size, min=1e-3)
+        params.add('c', value=.5, min=1e-3)
+
+        out = lmfit.minimize(residual, params, args=(x, y), method="leastsq")
+
+        if out.success:
+            cp = int(out.params["x0"])
+            if ret_details:
+                details["plot force"] = [x, y]
+                details["plot fit"] = [np.arange(force.size),
+                                       model(out.params, x)]
+                details["plot poc"] = [[cp, cp],
+                                       [y.min(), y.max()]]
+
+    if ret_details:
+        return cp, details
+    else:
+        return cp
 
 
 @poc(identifier="gradient_zero_crossing",
      name="Gradient zero-crossing of indentation part",
      preprocessing=["clip_approach"])
-def poc_gradient_zero_crossing(force):
+def poc_gradient_zero_crossing(force, ret_details=False):
     """Gradient zero-crossing of indentation part
 
     1. Apply a moving average filter to the curve
     2. Compute the gradient
-    3. Cut off gradient at maxiumum with a 10 point reserve
+    3. Cut off gradient at maximum with a 10 point reserve
     4. Apply a moving average filter to the gradient
     5. The POC is the index of the averaged gradient curve where
        the values are below 1% of the gradient maximum, measured
-       from the indentation maxium (not from baseline).
+       from the indentation maximum (not from baseline).
     """
     cp = np.nan
+    details = {}
     # Perform a median filter to smooth the array
     filtsize = max(5, int(force.size*.01))
     y = uniform_filter1d(force, size=filtsize)
@@ -167,7 +311,8 @@ def poc_gradient_zero_crossing(force):
         if grad.size > 50:
             # Use the point where the gradient becomes small enough.
             gradn = uniform_filter1d(grad, size=filtsize)
-            gradpos = gradn <= 0.01 * np.max(gradn)
+            thresh = 0.01 * np.max(gradn)
+            gradpos = gradn <= thresh
             if np.sum(gradpos):
                 # The gradient contains positive values.
                 # Flip `gradpos`, because we want the first value from the
@@ -177,4 +322,16 @@ def poc_gradient_zero_crossing(force):
                 # center of the array (and two times, because we did two
                 # filter operations).
                 cp = y.size - np.where(gradpos[::-1])[0][0] - cutoff + filtsize
-    return cp
+
+                if ret_details:
+                    x = np.arange(gradn.size)
+                    details["plot force gradient"] = [x, gradn]
+                    details["plot threshold"] = [[x[0], x[-1]],
+                                                 [thresh, thresh]]
+                    details["plot poc"] = [[cp, cp],
+                                           [gradn.min(), gradn.max()]]
+
+    if ret_details:
+        return cp, details
+    else:
+        return cp