Sync a GPX to the GPS clock of a timelapse, even one too noisy to publish

mapillary_tools/geotag/video_extractors/gpx.py

@@ -8,6 +8,7 @@
 import dataclasses
 import enum
 import logging
+import statistics
 import sys
 import typing as T
 from pathlib import Path
@@ -32,6 +33,11 @@
 # orders of magnitude.
 _IMPLAUSIBLE_OFFSET_SECONDS = 24 * 3600
 
+# A GPS clock that runs this much faster than the video clock is a timelapse.
+# Normal videos stay well below it, although over a few seconds of GPS
+# timestamped once per block the estimate can be 20% off
+_MIN_TIMELAPSE_RATE = 1.5
+
 
 class SyncMode(enum.Enum):
     # Sync by video GPS timestamps if found, otherwise rebase
@@ -42,6 +48,33 @@ class SyncMode(enum.Enum):
     REBASE = "rebase"
 
 
+@dataclasses.dataclass(frozen=True)
+class VideoClock:
+    """
+    Maps between the video time and the Unix time of the GPS clock in the video
+
+    >>> clock = VideoClock(video_time=2.0, unix_time=1000.0)
+    >>> clock.to_video_time(1030.0)
+    32.0
+    >>> timelapse = VideoClock(video_time=2.0, unix_time=1000.0, rate=15.0)
+    >>> timelapse.to_video_time(1030.0)
+    4.0
+    >>> timelapse.to_unix_time(4.0)
+    1030.0
+    """
+
+    video_time: float
+    unix_time: float
+    # GPS seconds per video second: 1 unless the video is a timelapse
+    rate: float = 1.0
+
+    def to_video_time(self, unix_time: float) -> float:
+        return self.video_time + (unix_time - self.unix_time) / self.rate
+
+    def to_unix_time(self, video_time: float) -> float:
+        return self.unix_time + (video_time - self.video_time) * self.rate
+
+
 class GPXVideoExtractor(BaseVideoExtractor):
     def __init__(
         self, video_path: Path, gpx_path: Path, sync_mode: SyncMode = SyncMode.SYNC
@@ -67,10 +100,8 @@ def extract(self) -> types.VideoMetadata:
 
         gpx_points: T.Sequence[geo.Point] = sum(gpx_tracks, [])
 
-        native_extractor = NativeVideoExtractor(self.video_path)
-
         try:
-            native_video_metadata = native_extractor.extract()
+            native_video_metadata, noisy = self._extract_native_video_metadata()
         except exceptions.MapillaryVideoGPSUnusableError as ex:
             # The GPX replaces the video's own GPS, so GPS that is missing,
             # empty or too noisy must not reject the video
@@ -84,23 +115,66 @@ def extract(self) -> types.VideoMetadata:
                 points=gpx_points,
             )
 
-        if self.sync_mode is SyncMode.REBASE:
+        clock: VideoClock | None = None
+        if self.sync_mode is not SyncMode.REBASE:
+            clock = self._read_video_clock(native_video_metadata.points)
+
+        if clock is not None and gpx_points:
+            offset = gpx_points[0].time - clock.to_unix_time(0)
+            if abs(offset) > _IMPLAUSIBLE_OFFSET_SECONDS:
+                if noisy and self.sync_mode is SyncMode.SYNC:
+                    # GPS rejected as noise may not have GPS time yet: GoPro
+                    # stamps a placeholder date until it acquires it. That GPS
+                    # is what the GPX replaces, so do not let its clock place
+                    # the GPX days away from the video
+                    LOG.warning(
+                        "The GPS clock of %s, whose GPS was rejected as noise, is %.1f days away from %s. "
+                        "Syncing the GPX from the start of the video instead",
+                        self.video_path,
+                        abs(offset) / 86400,
+                        self.gpx_path,
+                    )
+                    clock = None
+                else:
+                    LOG.warning(
+                        "Syncing %s against %s requires an offset of %.0f seconds (%.1f days). "
+                        "The GPX file probably does not belong to this video",
+                        self.video_path,
+                        self.gpx_path,
+                        offset,
+                        offset / 86400,
+                    )
+
+        if clock is None:
             self._rebase_times(gpx_points)
         else:
-            offset = self._gpx_offset(gpx_points, native_video_metadata.points)
-            if abs(offset) > _IMPLAUSIBLE_OFFSET_SECONDS:
-                LOG.warning(
-                    "Syncing %s against %s requires an offset of %.0f seconds (%.1f days). "
-                    "The GPX file probably does not belong to this video",
+            if clock.rate != 1.0:
+                LOG.info(
+                    "The GPS clock in %s runs %.1f times faster than the video (timelapse). Syncing %s accordingly",
                     self.video_path,
+                    clock.rate,
                     self.gpx_path,
-                    offset,
-                    offset / 86400,
                 )
-            self._rebase_times(gpx_points, offset=offset)
+            for p in gpx_points:
+                p.time = clock.to_video_time(p.time)
 
         return dataclasses.replace(native_video_metadata, points=gpx_points)
 
+    def _extract_native_video_metadata(self) -> tuple[types.VideoMetadata, bool]:
+        """
+        Extract the video metadata, and whether its GPS was rejected as noise
+        """
+        try:
+            return NativeVideoExtractor(self.video_path).extract(), False
+        except exceptions.MapillaryGPSNoiseError:
+            # The GPX replaces the video's own GPS, which is then only needed
+            # for make/model and a clock to sync against. Noisy points are
+            # never published
+            metadata = NativeVideoExtractor(
+                self.video_path, filter_noisy_points=False
+            ).extract()
+            return metadata, True
+
     @classmethod
     def _rebase_times(cls, points: T.Sequence[geo.Point], offset: float = 0.0) -> None:
         """
@@ -112,34 +186,94 @@ def _rebase_times(cls, points: T.Sequence[geo.Point], offset: float = 0.0) -> No
                 p.time = (p.time - first_timestamp) + offset
 
     @classmethod
-    def _gpx_offset(
-        cls, gpx_points: T.Sequence[geo.Point], video_gps_points: T.Sequence[geo.Point]
-    ) -> float:
+    def _read_video_clock(
+        cls, video_gps_points: T.Sequence[geo.Point]
+    ) -> VideoClock | None:
+        """
+        Read the GPS clock from the video GPS points, or None if they have no timestamps
+
+        >>> GPXVideoExtractor._read_video_clock([]) is None
+        True
+        >>> point = geo.Point(time=5, lat=1, lon=1, alt=None, angle=None)
+        >>> GPXVideoExtractor._read_video_clock([point]) is None
+        True
         """
-        Calculate the offset that needs to be applied to the GPX points to sync with the video GPS points.
-
-        >>> gpx_points = [geo.Point(time=5, lat=1, lon=1, alt=None, angle=None)]
-        >>> GPXVideoExtractor._gpx_offset(gpx_points, gpx_points)
-        0.0
-        >>> GPXVideoExtractor._gpx_offset(gpx_points, [])
-        0.0
-        >>> GPXVideoExtractor._gpx_offset([], gpx_points)
-        0.0
+        # (video time, Unix time) pairs
+        pairs: list[tuple[float, float]] = []
+
+        for p in video_gps_points:
+            # Video GPS timestamps are stored in whatever epoch their container
+            # uses (CAMM records GPS time, GoPro records Unix time), so go
+            # through get_unix_time(), which also skips zero/invalid timestamps
+            unix_time = p.get_unix_time()
+            if unix_time is None:
+                continue
+            if pairs and unix_time == pairs[-1][1]:
+                # GoPro GPS5 stamps a whole block (~1s) of points with one
+                # timestamp. Keep the first point of the block
+                continue
+            pairs.append((p.time, unix_time))
+
+        if not pairs:
+            return None
+
+        # Anchored at its own video time, which is not 0 when the video's GPS
+        # starts late (e.g. GoPro drops noisy fixes at the start of a recording)
+        video_time, unix_time = pairs[0]
+
+        rate = cls._estimate_clock_rate(pairs)
+        # One timestamp far off (e.g. a placeholder date before GPS time is
+        # acquired) can pull the least-squares slope far off, so a timelapse
+        # must show in the median step between neighbouring pairs too. That is
+        # too coarse to measure the rate, but not to tell a timelapse from a
+        # normal video
+        step_rate = cls._median_step_rate(pairs)
+        if (
+            rate is not None
+            and _MIN_TIMELAPSE_RATE < rate
+            and step_rate is not None
+            and _MIN_TIMELAPSE_RATE < step_rate
+        ):
+            return VideoClock(video_time=video_time, unix_time=unix_time, rate=rate)
+
+        return VideoClock(video_time=video_time, unix_time=unix_time)
+
+    @classmethod
+    def _estimate_clock_rate(
+        cls, pairs: T.Sequence[tuple[float, float]]
+    ) -> float | None:
         """
-        offset: float = 0.0
+        Least-squares slope of the Unix time over the video time, or None if the
+        pairs span less than a second of video
 
-        if not gpx_points or not video_gps_points:
-            return offset
+        >>> GPXVideoExtractor._estimate_clock_rate([(0, 100), (1, 115), (2, 130)])
+        15.0
+        >>> GPXVideoExtractor._estimate_clock_rate([(0, 100), (0.5, 107.5)]) is None
+        True
+        """
+        if len(pairs) < 2 or pairs[-1][0] - pairs[0][0] < 1.0:
+            return None
 
-        # Both sides must be Unix time here. Video GPS timestamps are stored in
-        # whatever epoch their container uses (CAMM records GPS time, GoPro
-        # records Unix time), so go through get_unix_time() rather than reading
-        # the raw attributes -- that also skips zero/invalid timestamps.
-        video_unix_time = video_gps_points[0].get_unix_time()
+        n = len(pairs)
+        mean_video_time = sum(v for v, _ in pairs) / n
+        mean_unix_time = sum(u for _, u in pairs) / n
+        covariance = sum((v - mean_video_time) * (u - mean_unix_time) for v, u in pairs)
+        variance = sum((v - mean_video_time) ** 2 for v, _ in pairs)
+        return covariance / variance
 
-        if video_unix_time is not None:
-            # Sync to video time 0, not to the first video GPS point, which can
-            # come later (GoPro drops noisy fixes at the start of a recording)
-            offset = gpx_points[0].time - (video_unix_time - video_gps_points[0].time)
+    @classmethod
+    def _median_step_rate(cls, pairs: T.Sequence[tuple[float, float]]) -> float | None:
+        """
+        Median rate of the GPS clock between neighbouring pairs, or None if
+        there are none
 
-        return offset
+        >>> pairs = [(0, 100), (1, 101), (2, 102), (3, 140), (4, 141)]
+        >>> GPXVideoExtractor._median_step_rate(pairs)
+        1.0
+        """
+        steps = [
+            (u2 - u1) / (v2 - v1)
+            for (v1, u1), (v2, u2) in zip(pairs, pairs[1:])
+            if v1 < v2
+        ]
+        return statistics.median(steps) if steps else None

mapillary_tools/geotag/video_extractors/native.py

@@ -22,6 +22,13 @@
 
 
 class GoProVideoExtractor(BaseVideoExtractor):
+    def __init__(self, video_path: Path, filter_noisy_points: bool = True):
+        super().__init__(video_path)
+        # The noise filter is a quality gate on the track we are about to
+        # publish. A caller that replaces the track (e.g. with a GPX) and only
+        # needs the make/model and the GPS clock of the video passes False
+        self.filter_noisy_points = filter_noisy_points
+
     @override
     def extract(self) -> types.VideoMetadata:
         with self.video_path.open("rb") as fp:
@@ -37,11 +44,13 @@ def extract(self) -> types.VideoMetadata:
         if not gps_points:
             raise exceptions.MapillaryGPXEmptyError("Empty GPS data found")
 
-        gps_points = T.cast(
-            T.List[telemetry.GPSPoint], gpmf_gps_filter.remove_noisy_points(gps_points)
-        )
-        if not gps_points:
-            raise exceptions.MapillaryGPSNoiseError("GPS is too noisy")
+        if self.filter_noisy_points:
+            gps_points = T.cast(
+                T.List[telemetry.GPSPoint],
+                gpmf_gps_filter.remove_noisy_points(gps_points),
+            )
+            if not gps_points:
+                raise exceptions.MapillaryGPSNoiseError("GPS is too noisy")
 
         video_metadata = types.VideoMetadata(
             filename=self.video_path,
@@ -106,17 +115,25 @@ def extract(self) -> types.VideoMetadata:
 
 
 class NativeVideoExtractor(BaseVideoExtractor):
-    def __init__(self, video_path: Path, filetypes: set[types.FileType] | None = None):
+    def __init__(
+        self,
+        video_path: Path,
+        filetypes: set[types.FileType] | None = None,
+        filter_noisy_points: bool = True,
+    ):
         super().__init__(video_path)
         self.filetypes = filetypes
+        self.filter_noisy_points = filter_noisy_points
 
     @override
     def extract(self) -> types.VideoMetadata:
         ft = self.filetypes
         extractor: BaseVideoExtractor
 
         if ft is None or types.FileType.VIDEO in ft or types.FileType.GOPRO in ft:
-            extractor = GoProVideoExtractor(self.video_path)
+            extractor = GoProVideoExtractor(
+                self.video_path, filter_noisy_points=self.filter_noisy_points
+            )
             try:
                 return extractor.extract()
             except simple_mp4_parser.BoxNotFoundError as ex:

tests/unit/test_gps_epoch.py

@@ -129,50 +129,48 @@ def test_invalid_timestamps_are_ignored(self):
         )
 
 
-class TestGPXOffset:
+class TestGPXSync:
     """A GPX recorded alongside the video must sync to ~0, not to ~10 years."""
 
     def test_camm_video_syncs_to_zero(self):
-        gpx_points = [_camm_point(time=A_UNIX_TIME, epoch_time=A_UNIX_TIME)]
         video_points = [_camm_point(time=0.0, epoch_time=A_UNIX_TIME)]
-        assert GPXVideoExtractor._gpx_offset(gpx_points, video_points) == 0.0
+        clock = GPXVideoExtractor._read_video_clock(video_points)
+        assert clock is not None
+        assert clock.to_video_time(A_UNIX_TIME) == 0.0
 
     def test_gopro_video_syncs_to_zero(self):
-        gpx_points = [_camm_point(time=A_UNIX_TIME, epoch_time=A_UNIX_TIME)]
         video_points = [_gps_point(time=0.0, epoch_time=A_UNIX_TIME)]
-        assert GPXVideoExtractor._gpx_offset(gpx_points, video_points) == 0.0
+        clock = GPXVideoExtractor._read_video_clock(video_points)
+        assert clock is not None
+        assert clock.to_video_time(A_UNIX_TIME) == 0.0
 
     def test_real_offset_is_preserved(self):
-        gpx_points = [_camm_point(time=A_UNIX_TIME + 30, epoch_time=A_UNIX_TIME + 30)]
         video_points = [_camm_point(time=0.0, epoch_time=A_UNIX_TIME)]
-        assert GPXVideoExtractor._gpx_offset(gpx_points, video_points) == 30.0
+        clock = GPXVideoExtractor._read_video_clock(video_points)
+        assert clock is not None
+        assert clock.to_video_time(A_UNIX_TIME + 30) == 30.0
 
-    def test_missing_video_timestamp_yields_no_offset(self):
-        gpx_points = [_camm_point(time=A_UNIX_TIME, epoch_time=A_UNIX_TIME)]
+    def test_missing_video_timestamp_yields_no_clock(self):
         video_points = [_camm_point(time=0.0, epoch_time=0.0)]
-        assert GPXVideoExtractor._gpx_offset(gpx_points, video_points) == 0.0
+        assert GPXVideoExtractor._read_video_clock(video_points) is None
 
     @pytest.mark.parametrize("video_point", [_camm_point, _gps_point])
     def test_video_gps_starting_late_syncs_to_video_time(self, video_point):
         """
-        The offset is to video time 0, not to the first video GPS point, or
-        the GPX track would land early by that point's time.
+        The clock is anchored at the video time of the first video GPS point,
+        not at video time 0, or the GPX track would land early by that time.
         """
-        gpx_points = [
-            _camm_point(time=A_UNIX_TIME + t, epoch_time=A_UNIX_TIME + t)
-            for t in range(21)
-        ]
         # The video's GPS starts 5s into the video
         video_points = [
             video_point(time=5.0 + t, epoch_time=A_UNIX_TIME + 5 + t) for t in range(10)
         ]
 
-        offset = GPXVideoExtractor._gpx_offset(gpx_points, video_points)
-        GPXVideoExtractor._rebase_times(gpx_points, offset=offset)
+        clock = GPXVideoExtractor._read_video_clock(video_points)
 
         # Recorded at the same instant as the first video GPS point, so it
         # lands at the same video time
-        assert gpx_points[5].time == 5.0
+        assert clock is not None
+        assert clock.to_video_time(A_UNIX_TIME + 5) == 5.0
 
 
 class TestEditListOverflow: