feat: qualify K1 local surface over time

src/k1link/compute/__init__.py

@@ -111,6 +111,7 @@ from .lidar_local_surface import (
     K1_LOCAL_SURFACE_FRAME_SCHEMA,
     K1_LOCAL_SURFACE_REPORT_SCHEMA,
     K1_LOCAL_SURFACE_SCHEMA,
+    K1_LOCAL_SURFACE_TIMELINE_SCHEMA,
     K1LocalSurfaceProfile,
     K1LocalSurfaceV1,
     build_k1_local_surface,
@@ -209,6 +210,7 @@ __all__ = [
     "K1_LOCAL_SURFACE_FRAME_SCHEMA",
     "K1_LOCAL_SURFACE_REPORT_SCHEMA",
     "K1_LOCAL_SURFACE_SCHEMA",
+    "K1_LOCAL_SURFACE_TIMELINE_SCHEMA",
     "LIDAR_FIELD_REVIEW_REPORT_SCHEMA",
     "LIDAR_FIELD_REVIEW_SCHEMA",
     "LIDAR_FIELD_REVIEW_WINDOW_SCHEMA",

src/k1link/compute/lidar_local_surface.py

@@ -28,6 +28,7 @@ from .lidar_field_review import (
 K1_LOCAL_SURFACE_SCHEMA: Final = "missioncore.k1-local-surface/v1"
 K1_LOCAL_SURFACE_REPORT_SCHEMA: Final = "missioncore.k1-local-surface-report/v1"
 K1_LOCAL_SURFACE_FRAME_SCHEMA: Final = "missioncore.k1-local-surface-frame/v1"
+K1_LOCAL_SURFACE_TIMELINE_SCHEMA: Final = "missioncore.k1-local-surface-timeline/v1"
 K1_LOCAL_SURFACE_ARRAYS_NAME: Final = "local-surface.npz"
 K1_LOCAL_SURFACE_REPORT_NAME: Final = "local-surface.json"
 K1_LOCAL_SURFACE_MANIFEST_NAME: Final = "manifest.json"
@@ -67,6 +68,12 @@ class K1LocalSurfaceProfile:
     obstacle_max_height_m: float = 3.5
     maximum_pose_binding_ms: float = 100.0
     maximum_slope_deg: float = 40.0
+    step_min_height_m: float = 0.07
+    step_max_height_m: float = 0.32
+    step_max_plane_residual_m: float = 0.45
+    temporal_height_jump_m: float = 0.03
+    temporal_slope_jump_deg: float = 0.5
+    temporal_roughness_jump_m: float = 0.015
 
     def __post_init__(self) -> None:
         numeric = (
@@ -82,6 +89,12 @@ class K1LocalSurfaceProfile:
             self.obstacle_max_height_m,
             self.maximum_pose_binding_ms,
             self.maximum_slope_deg,
+            self.step_min_height_m,
+            self.step_max_height_m,
+            self.step_max_plane_residual_m,
+            self.temporal_height_jump_m,
+            self.temporal_slope_jump_deg,
+            self.temporal_roughness_jump_m,
         )
         if (
             not self.profile_id.strip()
@@ -101,6 +114,11 @@ class K1LocalSurfaceProfile:
             or not self.obstacle_min_height_m < self.obstacle_max_height_m <= 20.0
             or not 1.0 <= self.maximum_pose_binding_ms <= 10_000.0
             or not 1.0 <= self.maximum_slope_deg < 90.0
+            or not 0.02 <= self.step_min_height_m < self.step_max_height_m
+            or not self.step_max_height_m <= self.step_max_plane_residual_m <= 2.0
+            or not 0.02 <= self.temporal_height_jump_m <= 2.0
+            or not 0.1 <= self.temporal_slope_jump_deg <= 45.0
+            or not 0.005 <= self.temporal_roughness_jump_m <= 1.0
         ):
             raise LidarGroundError("K1 local-surface profile is invalid")
 
@@ -125,6 +143,9 @@ class K1LocalSurfaceProfile:
                 "robust_mad_scale": self.robust_mad_scale,
                 "minimum_inlier_band_m": self.minimum_inlier_band_m,
                 "maximum_slope_deg": self.maximum_slope_deg,
+                "step_min_height_m": self.step_min_height_m,
+                "step_max_height_m": self.step_max_height_m,
+                "step_max_plane_residual_m": self.step_max_plane_residual_m,
             },
             "classification": {
                 "surface_band_m": self.surface_band_m,
@@ -137,6 +158,13 @@ class K1LocalSurfaceProfile:
                 "basis": "recorded-nearest-host-monotonic-arrival",
                 "maximum_age_ms": self.maximum_pose_binding_ms,
             },
+            "temporal_qualification": {
+                "prediction_input": "previous-ttl-window-only",
+                "current_frame_excluded_from_prediction": True,
+                "height_jump_m": self.temporal_height_jump_m,
+                "slope_jump_deg": self.temporal_slope_jump_deg,
+                "roughness_jump_m": self.temporal_roughness_jump_m,
+            },
             "authority": {
                 "commands_enabled": False,
                 "navigation_or_safety_accepted": False,
@@ -186,7 +214,7 @@ class K1LocalSurfaceV1:
     def _validate(self) -> None:
         frame_count = _nonnegative_int(self.identity.get("frame_count"), "frame count")
         point_count = _nonnegative_int(self.identity.get("point_count"), "point count")
-        required = {
+        baseline = {
             "frame_valid",
             "frame_failure_code",
             "plane_coefficients_map",
@@ -205,8 +233,25 @@ class K1LocalSurfaceV1:
             "point_class",
             "point_height_m",
         }
-        if set(self.arrays.files) != required:
+        qualification = {
+            "prediction_available",
+            "prediction_cell_count",
+            "prediction_residual_p50_m",
+            "prediction_residual_p95_m",
+            "prediction_inlier_fraction",
+            "height_delta_m",
+            "slope_delta_deg",
+            "roughness_delta_m",
+            "temporal_compared",
+            "temporal_jump",
+            "step_candidate_cell_count",
+            "step_candidate_point_count",
+            "point_step_candidate",
+        }
+        files = set(self.arrays.files)
+        if files not in (baseline, baseline | qualification):
             raise LidarGroundError("K1 local-surface arrays are incomplete")
+        self.has_temporal_qualification = qualification <= files
         vector_f64 = (
             "sensor_height_m",
             "slope_deg",
@@ -261,6 +306,63 @@ class K1LocalSurfaceV1:
                 or np.any(value < 0)
             ):
                 raise LidarGroundError(f"K1 local-surface {name} is invalid")
+        if self.has_temporal_qualification:
+            qualification_f64 = (
+                "prediction_residual_p50_m",
+                "prediction_residual_p95_m",
+                "prediction_inlier_fraction",
+                "height_delta_m",
+                "slope_delta_deg",
+                "roughness_delta_m",
+            )
+            qualification_i64 = (
+                "prediction_cell_count",
+                "step_candidate_cell_count",
+                "step_candidate_point_count",
+            )
+            for name in qualification_f64:
+                value = self.arrays[name]
+                if (
+                    value.shape != (frame_count,)
+                    or value.dtype != np.dtype("<f8")
+                    or not np.isfinite(value).all()
+                ):
+                    raise LidarGroundError(
+                        f"K1 local-surface qualification {name} is invalid"
+                    )
+            for name in qualification_i64:
+                value = self.arrays[name]
+                if (
+                    value.shape != (frame_count,)
+                    or value.dtype != np.dtype("<i8")
+                    or np.any(value < 0)
+                ):
+                    raise LidarGroundError(
+                        f"K1 local-surface qualification {name} is invalid"
+                    )
+            for name in (
+                "prediction_available",
+                "temporal_compared",
+                "temporal_jump",
+            ):
+                value = self.arrays[name]
+                if value.shape != (frame_count,) or value.dtype != np.dtype("?"):
+                    raise LidarGroundError(
+                        f"K1 local-surface qualification {name} is invalid"
+                    )
+            step_candidate = self.arrays["point_step_candidate"]
+            if (
+                step_candidate.shape != (point_count,)
+                or step_candidate.dtype != np.dtype("u1")
+                or np.any(step_candidate > 1)
+                or np.any(
+                    self.arrays["prediction_inlier_fraction"][
+                        self.arrays["prediction_available"]
+                    ]
+                    > 1
+                )
+            ):
+                raise LidarGroundError("K1 local-surface step candidates are invalid")
         valid = self.arrays["frame_valid"]
         if (
             np.any(self.arrays["frame_failure_code"][valid] != FRAME_VALID)
@@ -270,6 +372,40 @@ class K1LocalSurfaceV1:
             or self.identity.get("valid_frame_count") != int(np.count_nonzero(valid))
         ):
             raise LidarGroundError("K1 local-surface frame validity is inconsistent")
+        if self.has_temporal_qualification:
+            metrics = _object(
+                self.report.get("metrics"),
+                "K1 local-surface metrics",
+            )
+            qualification_report = _object(
+                metrics.get("temporal_qualification"),
+                "K1 local-surface temporal qualification",
+            )
+            prediction_report = _object(
+                qualification_report.get("prediction"),
+                "K1 local-surface prediction report",
+            )
+            stability_report = _object(
+                qualification_report.get("stability"),
+                "K1 local-surface stability report",
+            )
+            step_report = _object(
+                qualification_report.get("step_candidates"),
+                "K1 local-surface step report",
+            )
+            if (
+                prediction_report.get("current_frame_excluded") is not True
+                or prediction_report.get("sample_count")
+                != int(np.count_nonzero(self.arrays["prediction_available"]))
+                or stability_report.get("sample_count")
+                != int(np.count_nonzero(self.arrays["temporal_compared"]))
+                or stability_report.get("jump_count")
+                != int(np.count_nonzero(self.arrays["temporal_jump"]))
+                or step_report.get("is_ground_truth") is not False
+            ):
+                raise LidarGroundError(
+                    "K1 local-surface temporal report is inconsistent"
+                )
         authority = _object(self.report.get("authority"), "K1 local-surface authority")
         policy = _object(self.report.get("occupancy_policy"), "K1 local-surface policy")
         if (
@@ -292,17 +428,27 @@ class K1LocalSurfaceV1:
         start = int(offsets[frame_index])
         end = int(offsets[frame_index + 1])
         point_class = self.arrays["point_class"][start:end]
+        if self.has_temporal_qualification:
+            step_candidate = self.arrays["point_step_candidate"][start:end]
+        else:
+            step_candidate = np.zeros(end - start, dtype=np.uint8)
         counts = {
             "classified": int(np.count_nonzero(point_class)),
             "surface": int(np.count_nonzero(point_class == POINT_SURFACE)),
             "occupied": int(np.count_nonzero(point_class == POINT_OCCUPIED)),
             "below_surface": int(np.count_nonzero(point_class == POINT_BELOW_SURFACE)),
+            "step_candidate": int(np.count_nonzero(step_candidate)),
         }
         expected = {
             "classified": int(self.arrays["classified_point_count"][frame_index]),
             "surface": int(self.arrays["surface_point_count"][frame_index]),
             "occupied": int(self.arrays["occupied_point_count"][frame_index]),
             "below_surface": int(self.arrays["below_surface_point_count"][frame_index]),
+            "step_candidate": (
+                int(self.arrays["step_candidate_point_count"][frame_index])
+                if self.has_temporal_qualification
+                else 0
+            ),
         }
         if counts != expected:
             raise LidarGroundError("K1 local-surface frame counts are inconsistent")
@@ -328,6 +474,7 @@ class K1LocalSurfaceV1:
             "point_height_m": self.arrays["point_height_m"][start:end]
             .astype(np.float64)
             .tolist(),
+            "point_step_candidate": step_candidate.astype(np.int64).tolist(),
             "pose": {
                 "position_xyz_m": source.arrays["pose_positions_map"][frame_index]
                 .astype(np.float64)
@@ -356,6 +503,66 @@ class K1LocalSurfaceV1:
                 ),
             },
             "counts": counts,
+            "prediction": {
+                "available": (
+                    bool(self.arrays["prediction_available"][frame_index])
+                    if self.has_temporal_qualification
+                    else False
+                ),
+                "current_frame_excluded": True,
+                "cell_count": (
+                    int(self.arrays["prediction_cell_count"][frame_index])
+                    if self.has_temporal_qualification
+                    else 0
+                ),
+                "residual_p50_m": (
+                    float(self.arrays["prediction_residual_p50_m"][frame_index])
+                    if self.has_temporal_qualification
+                    else 0.0
+                ),
+                "residual_p95_m": (
+                    float(self.arrays["prediction_residual_p95_m"][frame_index])
+                    if self.has_temporal_qualification
+                    else 0.0
+                ),
+                "inlier_fraction": (
+                    float(self.arrays["prediction_inlier_fraction"][frame_index])
+                    if self.has_temporal_qualification
+                    else 0.0
+                ),
+            },
+            "temporal": {
+                "compared": (
+                    bool(self.arrays["temporal_compared"][frame_index])
+                    if self.has_temporal_qualification
+                    else False
+                ),
+                "height_delta_m": (
+                    float(self.arrays["height_delta_m"][frame_index])
+                    if self.has_temporal_qualification
+                    else 0.0
+                ),
+                "slope_delta_deg": (
+                    float(self.arrays["slope_delta_deg"][frame_index])
+                    if self.has_temporal_qualification
+                    else 0.0
+                ),
+                "roughness_delta_m": (
+                    float(self.arrays["roughness_delta_m"][frame_index])
+                    if self.has_temporal_qualification
+                    else 0.0
+                ),
+                "jump": (
+                    bool(self.arrays["temporal_jump"][frame_index])
+                    if self.has_temporal_qualification
+                    else False
+                ),
+                "step_candidate_cell_count": (
+                    int(self.arrays["step_candidate_cell_count"][frame_index])
+                    if self.has_temporal_qualification
+                    else 0
+                ),
+            },
             "classes": {
                 "0": "unclassified-or-outside-local-radius",
                 "1": "observed-surface",
@@ -368,6 +575,62 @@ class K1LocalSurfaceV1:
             "authority": self.report["authority"],
         }
 
+    def timeline_detail(self, source: E10LidarFieldSource) -> dict[str, object]:
+        _validate_source_binding(self, source)
+        frame_count = source.frame_count
+        if self.has_temporal_qualification:
+            prediction_available = self.arrays["prediction_available"]
+            temporal_compared = self.arrays["temporal_compared"]
+            temporal_jump = self.arrays["temporal_jump"]
+            prediction_p50 = self.arrays["prediction_residual_p50_m"]
+            prediction_p95 = self.arrays["prediction_residual_p95_m"]
+            prediction_inlier = self.arrays["prediction_inlier_fraction"]
+            height_delta = self.arrays["height_delta_m"]
+            slope_delta = self.arrays["slope_delta_deg"]
+            roughness_delta = self.arrays["roughness_delta_m"]
+            step_points = self.arrays["step_candidate_point_count"]
+        else:
+            prediction_available = np.zeros(frame_count, dtype=np.bool_)
+            temporal_compared = np.zeros(frame_count, dtype=np.bool_)
+            temporal_jump = np.zeros(frame_count, dtype=np.bool_)
+            prediction_p50 = np.zeros(frame_count, dtype=np.float64)
+            prediction_p95 = np.zeros(frame_count, dtype=np.float64)
+            prediction_inlier = np.zeros(frame_count, dtype=np.float64)
+            height_delta = np.zeros(frame_count, dtype=np.float64)
+            slope_delta = np.zeros(frame_count, dtype=np.float64)
+            roughness_delta = np.zeros(frame_count, dtype=np.float64)
+            step_points = np.zeros(frame_count, dtype=np.int64)
+        return {
+            "schema_version": K1_LOCAL_SURFACE_TIMELINE_SCHEMA,
+            "model_id": self.model_id,
+            "source_pack_id": source.pack_id,
+            "session_id": source.identity["session_id"],
+            "frame_count": frame_count,
+            "source_frame_index": source.arrays["source_frame_indices"]
+            .astype(np.int64)
+            .tolist(),
+            "session_seconds": source.arrays["session_seconds"].astype(np.float64).tolist(),
+            "source_available": source.arrays["sample_available"].astype(np.int64).tolist(),
+            "valid": self.arrays["frame_valid"].astype(np.int64).tolist(),
+            "prediction_available": prediction_available.astype(np.int64).tolist(),
+            "prediction_residual_p50_m": prediction_p50.astype(np.float64).tolist(),
+            "prediction_residual_p95_m": prediction_p95.astype(np.float64).tolist(),
+            "prediction_inlier_fraction": prediction_inlier.astype(np.float64).tolist(),
+            "sensor_height_m": self.arrays["sensor_height_m"].astype(np.float64).tolist(),
+            "slope_deg": self.arrays["slope_deg"].astype(np.float64).tolist(),
+            "roughness_m": self.arrays["roughness_m"].astype(np.float64).tolist(),
+            "confidence": self.arrays["confidence"].astype(np.float64).tolist(),
+            "temporal_compared": temporal_compared.astype(np.int64).tolist(),
+            "height_delta_m": height_delta.astype(np.float64).tolist(),
+            "slope_delta_deg": slope_delta.astype(np.float64).tolist(),
+            "roughness_delta_m": roughness_delta.astype(np.float64).tolist(),
+            "temporal_jump": temporal_jump.astype(np.int64).tolist(),
+            "step_candidate_point_count": step_points.astype(np.int64).tolist(),
+            "ground_truth": False,
+            "access": "read-only",
+            "authority": self.report["authority"],
+        }
+
 
 def build_k1_local_surface(
     source: E10LidarFieldSource,
@@ -392,6 +655,7 @@ def build_k1_local_surface(
     times = source_arrays["session_seconds"]
     pose_delta = np.abs(source_arrays["pose_point_delta_ms"])
     cache: dict[tuple[int, int], tuple[float, float]] = {}
+    previous_surface: tuple[float, float, float, float] | None = None
 
     for frame_index in range(frame_count):
         start = int(offsets[frame_index])
@@ -416,8 +680,27 @@ def build_k1_local_surface(
         )
         local_cloud = cloud[local]
         _expire_cache(cache, session_seconds, position, profile)
+        _, prior_cell_points, _ = _local_cache_records(cache, position, profile)
+        _, current_cell_points = _cloud_cell_observations(local_cloud, profile)
+        prediction = _prediction_metrics(
+            prior_cell_points,
+            current_cell_points,
+            position,
+            profile,
+        )
+        if prediction is not None:
+            residual_p50, residual_p95, inlier_fraction, prediction_cells = prediction
+            arrays["prediction_available"][frame_index] = True
+            arrays["prediction_cell_count"][frame_index] = prediction_cells
+            arrays["prediction_residual_p50_m"][frame_index] = residual_p50
+            arrays["prediction_residual_p95_m"][frame_index] = residual_p95
+            arrays["prediction_inlier_fraction"][frame_index] = inlier_fraction
         _update_cache(cache, local_cloud, session_seconds, profile)
-        cell_points, cell_times = _local_cache_points(cache, position, profile)
+        cell_keys, cell_points, cell_times = _local_cache_records(
+            cache,
+            position,
+            profile,
+        )
         arrays["surface_cell_count"][frame_index] = cell_points.shape[0]
         if cell_points.shape[0] < profile.minimum_surface_cells:
             arrays["frame_failure_code"][frame_index] = FRAME_INSUFFICIENT_SURFACE
@@ -442,6 +725,14 @@ def build_k1_local_surface(
             & (heights <= profile.obstacle_max_height_m)
         ] = POINT_OCCUPIED
         local_classes[local & (heights < -profile.surface_band_m)] = POINT_BELOW_SURFACE
+        step_keys = _step_candidate_keys(cell_keys, cell_points, plane, profile)
+        step_candidate = _point_step_candidates(
+            cloud,
+            local,
+            heights,
+            step_keys,
+            profile,
+        )
         point_heights = np.zeros(cloud.shape[0], dtype=np.float32)
         point_heights[local] = heights[local].astype(np.float32)
         sensor_height = float(_height_above_plane(position.reshape(1, 3), plane)[0])
@@ -453,6 +744,23 @@ def build_k1_local_surface(
             1.0 - float(pose_delta[frame_index]) / profile.maximum_pose_binding_ms,
         )
         confidence = float(np.clip(coverage * roughness_confidence * pose_confidence, 0.0, 1.0))
+        if (
+            previous_surface is not None
+            and session_seconds - previous_surface[0] <= profile.surface_ttl_s
+        ):
+            height_delta = abs(sensor_height - previous_surface[1])
+            slope_delta = abs(slope_deg - previous_surface[2])
+            roughness_delta = abs(roughness - previous_surface[3])
+            arrays["height_delta_m"][frame_index] = height_delta
+            arrays["slope_delta_deg"][frame_index] = slope_delta
+            arrays["roughness_delta_m"][frame_index] = roughness_delta
+            arrays["temporal_compared"][frame_index] = True
+            arrays["temporal_jump"][frame_index] = (
+                height_delta > profile.temporal_height_jump_m
+                or slope_delta > profile.temporal_slope_jump_deg
+                or roughness_delta > profile.temporal_roughness_jump_m
+            )
+        previous_surface = (session_seconds, sensor_height, slope_deg, roughness)
         arrays["frame_valid"][frame_index] = True
         arrays["frame_failure_code"][frame_index] = FRAME_VALID
         arrays["plane_coefficients_map"][frame_index] = plane
@@ -467,6 +775,7 @@ def build_k1_local_surface(
         arrays["surface_inlier_cell_count"][frame_index] = int(np.count_nonzero(inliers))
         arrays["point_class"][start:end] = local_classes
         arrays["point_height_m"][start:end] = point_heights
+        arrays["point_step_candidate"][start:end] = step_candidate
         arrays["classified_point_count"][frame_index] = int(
             np.count_nonzero(local_classes)
         )
@@ -479,6 +788,10 @@ def build_k1_local_surface(
         arrays["below_surface_point_count"][frame_index] = int(
             np.count_nonzero(local_classes == POINT_BELOW_SURFACE)
         )
+        arrays["step_candidate_cell_count"][frame_index] = len(step_keys)
+        arrays["step_candidate_point_count"][frame_index] = int(
+            np.count_nonzero(step_candidate)
+        )
 
     logical_content_sha256 = _logical_sha256(arrays)
     valid = arrays["frame_valid"]
@@ -579,6 +892,58 @@ def build_k1_local_surface(
             "surface_max_age_ms": _valid_distribution(
                 arrays["surface_max_age_ms"], valid
             ),
+            "temporal_qualification": {
+                "prediction": {
+                    "current_frame_excluded": True,
+                    "sample_count": int(
+                        np.count_nonzero(arrays["prediction_available"])
+                    ),
+                    "residual_p50_m": _valid_distribution(
+                        arrays["prediction_residual_p50_m"],
+                        arrays["prediction_available"],
+                    ),
+                    "residual_p95_m": _valid_distribution(
+                        arrays["prediction_residual_p95_m"],
+                        arrays["prediction_available"],
+                    ),
+                    "inlier_fraction": _valid_distribution(
+                        arrays["prediction_inlier_fraction"],
+                        arrays["prediction_available"],
+                    ),
+                },
+                "stability": {
+                    "sample_count": int(
+                        np.count_nonzero(arrays["temporal_compared"])
+                    ),
+                    "height_delta_m": _valid_distribution(
+                        arrays["height_delta_m"],
+                        arrays["temporal_compared"],
+                    ),
+                    "slope_delta_deg": _valid_distribution(
+                        arrays["slope_delta_deg"],
+                        arrays["temporal_compared"],
+                    ),
+                    "roughness_delta_m": _valid_distribution(
+                        arrays["roughness_delta_m"],
+                        arrays["temporal_compared"],
+                    ),
+                    "jump_count": int(np.count_nonzero(arrays["temporal_jump"])),
+                },
+                "step_candidates": {
+                    "is_ground_truth": False,
+                    "frames_with_candidates": int(
+                        np.count_nonzero(arrays["step_candidate_cell_count"] > 0)
+                    ),
+                    "cell_count": _valid_distribution(
+                        arrays["step_candidate_cell_count"].astype(np.float64),
+                        valid,
+                    ),
+                    "point_count": _valid_distribution(
+                        arrays["step_candidate_point_count"].astype(np.float64),
+                        valid,
+                    ),
+                },
+            },
             "build_elapsed_ms": (time.perf_counter() - started) * 1_000.0,
         },
         "anchors": _anchors(source, valid),
@@ -666,8 +1031,21 @@ def _empty_arrays(
         "surface_point_count": np.zeros(frame_count, dtype="<i8"),
         "occupied_point_count": np.zeros(frame_count, dtype="<i8"),
         "below_surface_point_count": np.zeros(frame_count, dtype="<i8"),
+        "prediction_available": np.zeros(frame_count, dtype="?"),
+        "prediction_cell_count": np.zeros(frame_count, dtype="<i8"),
+        "prediction_residual_p50_m": np.zeros(frame_count, dtype="<f8"),
+        "prediction_residual_p95_m": np.zeros(frame_count, dtype="<f8"),
+        "prediction_inlier_fraction": np.zeros(frame_count, dtype="<f8"),
+        "height_delta_m": np.zeros(frame_count, dtype="<f8"),
+        "slope_delta_deg": np.zeros(frame_count, dtype="<f8"),
+        "roughness_delta_m": np.zeros(frame_count, dtype="<f8"),
+        "temporal_compared": np.zeros(frame_count, dtype="?"),
+        "temporal_jump": np.zeros(frame_count, dtype="?"),
+        "step_candidate_cell_count": np.zeros(frame_count, dtype="<i8"),
+        "step_candidate_point_count": np.zeros(frame_count, dtype="<i8"),
         "point_class": np.zeros(point_count, dtype="u1"),
         "point_height_m": np.zeros(point_count, dtype="<f4"),
+        "point_step_candidate": np.zeros(point_count, dtype="u1"),
     }
 
 
@@ -677,8 +1055,17 @@ def _update_cache(
     session_seconds: float,
     profile: K1LocalSurfaceProfile,
 ) -> None:
+    keys, points = _cloud_cell_observations(cloud, profile)
+    for key, point in zip(keys, points, strict=True):
+        cache[(int(key[0]), int(key[1]))] = (float(point[2]), session_seconds)
+
+
+def _cloud_cell_observations(
+    cloud: npt.NDArray[np.float64],
+    profile: K1LocalSurfaceProfile,
+) -> tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]]:
     if cloud.shape[0] == 0:
-        return
+        return np.empty((0, 2), dtype=np.int64), np.empty((0, 3), dtype=np.float64)
     cells = np.floor(cloud[:, :2] / profile.cell_size_m).astype(np.int64)
     order = np.lexsort((cells[:, 1], cells[:, 0]))
     sorted_cells = cells[order]
@@ -686,10 +1073,17 @@ def _update_cache(
     changes = np.flatnonzero(np.any(np.diff(sorted_cells, axis=0) != 0, axis=1)) + 1
     starts = np.concatenate((np.asarray([0]), changes))
     ends = np.concatenate((changes, np.asarray([cloud.shape[0]])))
-    for start, end in zip(starts, ends, strict=True):
-        key = (int(sorted_cells[start, 0]), int(sorted_cells[start, 1]))
-        z = float(np.percentile(sorted_z[start:end], profile.cell_lower_percentile))
-        cache[key] = (z, session_seconds)
+    keys = np.empty((starts.shape[0], 2), dtype=np.int64)
+    points = np.empty((starts.shape[0], 3), dtype=np.float64)
+    half_cell = profile.cell_size_m * 0.5
+    for index, (start, end) in enumerate(zip(starts, ends, strict=True)):
+        keys[index] = sorted_cells[start]
+        points[index] = (
+            float(sorted_cells[start, 0]) * profile.cell_size_m + half_cell,
+            float(sorted_cells[start, 1]) * profile.cell_size_m + half_cell,
+            float(np.percentile(sorted_z[start:end], profile.cell_lower_percentile)),
+        )
+    return keys, points
 
 
 def _expire_cache(
@@ -717,11 +1111,15 @@ def _expire_cache(
         del cache[key]
 
 
-def _local_cache_points(
+def _local_cache_records(
     cache: Mapping[tuple[int, int], tuple[float, float]],
     position: npt.NDArray[np.float64],
     profile: K1LocalSurfaceProfile,
-) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
+) -> tuple[
+    npt.NDArray[np.int64],
+    npt.NDArray[np.float64],
+    npt.NDArray[np.float64],
+]:
     values = [
         (
             (key[0] + 0.5) * profile.cell_size_m,
@@ -737,9 +1135,14 @@ def _local_cache_points(
         )
     ]
     if not values:
-        return np.empty((0, 3), dtype=np.float64), np.empty(0, dtype=np.float64)
+        return (
+            np.empty((0, 2), dtype=np.int64),
+            np.empty((0, 3), dtype=np.float64),
+            np.empty(0, dtype=np.float64),
+        )
     array = np.asarray(values, dtype=np.float64)
-    return array[:, :3], array[:, 3]
+    keys = np.floor(array[:, :2] / profile.cell_size_m).astype(np.int64)
+    return keys, array[:, :3], array[:, 3]
 
 
 def _fit_surface(
@@ -809,6 +1212,92 @@ def _fit_surface(
     return plane.astype("<f8"), inliers, residuals
 
 
+def _prediction_metrics(
+    prior_cell_points: npt.NDArray[np.float64],
+    current_cell_points: npt.NDArray[np.float64],
+    position: npt.NDArray[np.float64],
+    profile: K1LocalSurfaceProfile,
+) -> tuple[float, float, float, int] | None:
+    """Score current lower-cell evidence against a plane built without that frame."""
+
+    if (
+        prior_cell_points.shape[0] < profile.minimum_surface_cells
+        or current_cell_points.shape[0] < profile.minimum_surface_cells
+    ):
+        return None
+    prior_fit = _fit_surface(prior_cell_points, position, profile)
+    if prior_fit is None:
+        return None
+    prior_plane, _, _ = prior_fit
+    cutoff = float(
+        np.quantile(current_cell_points[:, 2], profile.initial_lower_fraction)
+    )
+    evaluation = current_cell_points[:, 2] <= cutoff
+    if int(np.count_nonzero(evaluation)) < profile.minimum_surface_cells:
+        return None
+    residual = np.abs(
+        _height_above_plane(current_cell_points[evaluation], prior_plane)
+    )
+    if residual.size == 0 or not np.isfinite(residual).all():
+        return None
+    return (
+        float(np.percentile(residual, 50)),
+        float(np.percentile(residual, 95)),
+        float(np.mean(residual <= profile.surface_band_m)),
+        int(residual.shape[0]),
+    )
+
+
+def _step_candidate_keys(
+    cell_keys: npt.NDArray[np.int64],
+    cell_points: npt.NDArray[np.float64],
+    plane: npt.NDArray[np.float64],
+    profile: K1LocalSurfaceProfile,
+) -> set[tuple[int, int]]:
+    """Find local discontinuities; the result remains an unverified candidate."""
+
+    if cell_keys.shape[0] != cell_points.shape[0]:
+        raise LidarGroundError("K1 local-surface cell alignment is invalid")
+    residual = _height_above_plane(cell_points, plane)
+    lookup = {
+        (int(key[0]), int(key[1])): float(value)
+        for key, value in zip(cell_keys, residual, strict=True)
+        if abs(float(value)) <= profile.step_max_plane_residual_m
+    }
+    candidates: set[tuple[int, int]] = set()
+    for key, value in lookup.items():
+        for neighbor in ((key[0] + 1, key[1]), (key[0], key[1] + 1)):
+            neighbor_value = lookup.get(neighbor)
+            if neighbor_value is None:
+                continue
+            delta = abs(value - neighbor_value)
+            if profile.step_min_height_m <= delta <= profile.step_max_height_m:
+                candidates.add(key)
+                candidates.add(neighbor)
+    return candidates
+
+
+def _point_step_candidates(
+    cloud: npt.NDArray[np.float64],
+    local: npt.NDArray[np.bool_],
+    heights: npt.NDArray[np.float64],
+    candidate_keys: set[tuple[int, int]],
+    profile: K1LocalSurfaceProfile,
+) -> npt.NDArray[np.uint8]:
+    result = np.zeros(cloud.shape[0], dtype=np.uint8)
+    if not candidate_keys:
+        return result
+    cells = np.floor(cloud[:, :2] / profile.cell_size_m).astype(np.int64)
+    for index in np.flatnonzero(local):
+        key = (int(cells[index, 0]), int(cells[index, 1]))
+        if (
+            key in candidate_keys
+            and abs(float(heights[index])) <= profile.step_max_plane_residual_m
+        ):
+            result[index] = 1
+    return result
+
+
 def _height_above_plane(
     points: npt.NDArray[np.float64],
     plane: npt.NDArray[np.float64],