scale smoothing by signal-to-noise ratio

1 files changed, 53 insertions, 7 deletions

diff --git a/fietsboek/transformers/elevation.py b/fietsboek/transformers/elevation.py
index 6bdcbb4..15a5061 100644
--- a/fietsboek/transformers/elevation.py
+++ b/fietsboek/transformers/elevation.py
@@ -1,6 +1,7 @@
 """Transformers that deal with elevation changes in the track."""
 
 import logging
+import statistics
 from collections.abc import Callable, Iterable
 from itertools import islice, zip_longest
 
@@ -155,7 +156,39 @@ class FixElevationJumps(Transformer):
 
 
 class SmoothElevation(Transformer):
-    """A transformer that smooths the elevation."""
+    """A transformer that smooths the elevation.
+
+    **Technical description:**
+
+    This transformer has the goal of smoothing out noise in the GPS elevation
+    recording, while keeping the overall elevation profile true to the original.
+    We do this in two parts: First, we estimate the noisiness of the input data.
+    Then we apply a `Savitzky-Golay`_ filter, whose strength is scaled by the
+    noisiness estimated in the first step.
+
+    .. _Savitzky-Golay: https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter
+
+    **Step 1:** Noise estimation
+
+    To estimate the noisiness, we walk through the file in windows of 5 points.
+    We do a linear regression through those points, and compute the sum of
+    absolute errors.
+
+    The reason behind this is that we assume that the profile given by 5 points
+    can be approximated reasonably with a straight line, so deviations from this
+    line are likely caused by noise in the recording.
+
+    We average the sum-of-errors over all windows in the track to reach the
+    track's total noisiness, this prevents outlier windows from affecting the
+    noisiness too much.
+
+    **Step 2:** Filter
+
+    The actual smoothing is done with a pass of the Savitzky-Golay filter. The
+    window size :math:`w` is scaled linearly according to the noisiness
+    :math:`w` determined in step 1 as :math:`w = n \\cdot f`. The factor
+    :math:`f` is fixed and determined manually with some example tracks.
+    """
 
     @classmethod
     def identifier(cls) -> str:
@@ -191,13 +224,26 @@ class SmoothElevation(Transformer):
         if radius:
             self._savitzky_golay(path, radius)
 
+    def _snr(self, path: geo.Path) -> float:
+        window_size = 5
+        errors = []
+        for start in range(0, len(path.points) - window_size):
+            elevations = [point.elevation for point in path.points[start : start + window_size]]
+            xs = list(range(len(elevations)))
+
+            sl, intercept = statistics.linear_regression(xs, elevations)
+            errors.extend(
+                abs(elevation - (x * sl + intercept))
+                for x, elevation in enumerate(elevations)
+            )
+        return statistics.mean(errors)
+
     def _estimate_window_size(self, path: geo.Path) -> int:
-        length = 0.0
-        for i, (a, b) in enumerate(zip(path.points, path.points[1:])):
-            length += a.flat_distance(b)
-            if length >= ELEVATION_SMOOTHING_WINDOW_LENGTH:
-                return i
-        return 0
+        points_per_meter = len(path.points) / path.movement_data().length
+        snr = self._snr(path)
+        window_length = 6 * ELEVATION_SMOOTHING_WINDOW_LENGTH * snr
+        LOGGER.debug("signal-to-noise=%s, points_per_meter=%s", snr, points_per_meter)
+        return int(points_per_meter * window_length)
 
     def _savitzky_golay(self, path: geo.Path, radius: int):
         elevations = [0.0] * len(path.points)