refacto function with test

Author: mdupaysign

lidar_for_fuel/pretreatment/filter_deviation_day.py

@@ -12,6 +12,66 @@
 logger = logging.getLogger(__name__)
 
 _SECONDS_PER_DAY = 86_400.0
+_EPSILON = 1e-3  # 1 ms — smaller than any realistic GpsTime resolution
+
+
+def gpstime_to_day_index(gpstime: np.ndarray, gpstime_ref_unix: float) -> np.ndarray:
+    """Convert a GpsTime array (relative seconds) to integer UTC day indices.
+
+    Each value is floored to midnight UTC, expressed as *days since the
+    UNIX epoch* (1970-01-01).
+
+    Args:
+        gpstime (np.ndarray): 1-D array of relative GPS timestamps (seconds).
+        gpstime_ref_unix (float): UNIX timestamp of the GPS reference epoch.
+
+    Returns:
+        np.ndarray: Integer array of shape ``(n,)`` with UNIX day indices.
+
+    """
+    unix_time = gpstime + gpstime_ref_unix
+    return np.floor(unix_time / _SECONDS_PER_DAY).astype(np.int64)
+
+
+def modal_day(day_index: np.ndarray) -> int:
+    """Return the UNIX day index of the most densely sampled calendar day.
+
+    Args:
+        day_index (np.ndarray): Integer array of UNIX day indices (one per point).
+
+    Returns:
+        int: UNIX day index of the mode.
+
+    """
+    unique_days, counts = np.unique(day_index, return_counts=True)
+    return int(unique_days[counts.argmax()])
+
+
+def compute_gpstime_window(
+    main_day: int,
+    deviation_day: int | float,
+    gpstime_ref_unix: float,
+) -> tuple[float, float]:
+    """Compute the GpsTime [t_min, t_max] filter window.
+
+    The window covers all calendar days in
+    ``[main_day - deviation_day, main_day + deviation_day]``.
+    ``t_max`` is nudged inward by ``_EPSILON`` so that a point landing
+    exactly on midnight of ``day_hi + 1`` is excluded.
+
+    Args:
+        main_day (int): UNIX day index of the modal acquisition day.
+        deviation_day (int | float): Half-width of the window in days.
+        gpstime_ref_unix (float): UNIX timestamp of the GPS reference epoch.
+
+    Returns:
+        tuple[float, float]: ``(t_min, t_max)`` in GpsTime space
+    """
+    day_lo = main_day - int(deviation_day)
+    day_hi = main_day + int(deviation_day)
+    t_min = max(day_lo, 0) * _SECONDS_PER_DAY - gpstime_ref_unix
+    t_max = (day_hi + 1) * _SECONDS_PER_DAY - gpstime_ref_unix - _EPSILON
+    return t_min, t_max
 
 
 def filter_points_by_date(
@@ -56,29 +116,13 @@ def filter_points_by_date(
 
     n_total = len(points)
 
-    # Convert relative GPS time to absolute UNIX timestamps
-    unix_time = points["GpsTime"] + gpstime_ref_unix
-
-    # Floor to calendar day (midnight UTC)
-    day_index = np.floor(unix_time / _SECONDS_PER_DAY).astype(np.int64)  # days since epoch
-
-    # Find the most densely sampled day (mode)
-    unique_days, counts = np.unique(day_index, return_counts=True)
-    id_max_count = counts.argmax()  # index of the modal day
-    main_day = unique_days[id_max_count]  # day index of the modal day
-
-    # Build the window of valid day indices
-    day_lo = main_day - int(deviation_day)
-    day_hi = main_day + int(deviation_day)
-
-    # Borne haute inclusive: include the *entire* last day
-    # PDAL filters.range uses a closed interval [t_min, t_max], so we subtract
-    # a small epsilon to exclude any point falling exactly on midnight of day_hi+1.
-    _EPSILON = 1e-3  # 1 ms — smaller than any realistic GpsTime resolution
-    t_min = max(day_lo, 0) * _SECONDS_PER_DAY - gpstime_ref_unix  # back to GpsTime space
-    t_max = (day_hi + 1) * _SECONDS_PER_DAY - gpstime_ref_unix - _EPSILON  # strictly excludes day_hi+1
+    day_index = gpstime_to_day_index(points["GpsTime"], gpstime_ref_unix)
+    main_day_idx = modal_day(day_index)
+    t_min, t_max = compute_gpstime_window(main_day_idx, deviation_day, gpstime_ref_unix)
 
-    # Warning: percentage of points removed
+    unix_time = points["GpsTime"] + gpstime_ref_unix
+    day_lo = main_day_idx - int(deviation_day)
+    day_hi = main_day_idx + int(deviation_day)
     n_retained = int(
         np.sum((unix_time >= max(day_lo, 0) * _SECONDS_PER_DAY) & (unix_time < (day_hi + 1) * _SECONDS_PER_DAY))
     )
@@ -92,18 +136,16 @@ def filter_points_by_date(
             stacklevel=2,
         )
 
+    unique_days = np.unique(day_index)
     logger.debug(
         "Main day index: %d | GpsTime window [%.1f, %.1f) | %d/%d days retained | %.1f %% points removed",
-        main_day,
+        main_day_idx,
         t_min,
         t_max,
         int(np.sum((unique_days >= day_lo) & (unique_days <= day_hi))),
         len(unique_days),
         pct_removed,
     )
 
-    # Create PDAL pipeline using range filter for timestamps
     filter_pipeline_json = {"pipeline": [{"type": "filters.range", "limits": f"GpsTime[{t_min}:{t_max}]"}]}
-
-    # Return a PDAL Pipeline object
     return pdal.Pipeline(json.dumps(filter_pipeline_json), arrays=[points])