feat(benchmark): summarize SCARED reconstruction depth

- report valid coverage and depth percentiles even when XYZ ground truth is absent

- accept depth bounds and optional color-mapped PNG output for visual inspection

- retain accuracy metrics when point_cloud.obj is available

- add a batch helper that renders Markdown tables and optional JSON and PNG artifacts

- document benchmark inputs, dataset limits, and visualization tradeoffs

README.md

@@ -175,6 +175,8 @@ cloud_point:
   num_disparities: 160  # fx ~1024 px, baseline ~4.35 mm -> up to ~160 px
   min_depth_m: 0.02     # endoscopic working range: 20 mm .. 300 mm
   max_depth_m: 0.30
+  ply_stride: 1         # option 5: 4 = 16x smaller, block-averaged mesh
+  wls_filter: false     # true = smoother mesh for viewers, less accurate
 ```
 
 Depth limits are a physical bound on the scene: with this rig depth is
@@ -198,16 +200,48 @@ slivers across the surface and makes a correct cloud look like a fan of
 rays. Pass `PlyOptions{false, 0.0f, false}` to `write_ply` for a points-only ASCII
 file.
 
+**Surface roughness vs accuracy.** SGBM's sub-pixel disparity noise
+(~0.2 px, correlated over several pixels) is ~1 mm of depth on this rig,
+far more than the 0.07 mm lateral pixel pitch, so a mesh built from the raw
+cloud is "hairy": face normals sit a median 45° off the camera axis on
+tissue that faces the camera. Colour depth maps hide this; shaded mesh
+viewers show it as fuzz. `wls_filter: true` applies OpenCV's edge-aware
+WLS disparity filter (needs `opencv_ximgproc`, doubles matching time) and
+brings the median normal to ~23° with neighbour depth jumps down from 0.15 mm
+to 0.06 mm, but it also costs accuracy on SCARED (keyframe 1: MAE 0.84 →
+0.92 mm, within 2 mm 80 → 78 %; dataset_3: MAE 1.75 → 2.20 mm). It is
+therefore off by default: use it for pictures, not for measurements. SGBM's
+own holes are never filled by the filter.
+
 When the keyframe contains `point_cloud.obj`, run the benchmark to compare the
 reconstruction with its pixel-aligned XYZ ground truth:
 
 ```bash
 ./build/src/cloud_point/scared_dataset_benchmark \
-    /path/to/test_dataset_8/keyframe_0 160
+    /path/to/dataset_1/keyframe_1 160 [depth.png|-] [min_depth_m max_depth_m]
 ```
 
-The benchmark reports coverage, component-wise and 3-D errors, threshold
-accuracy, and matching/reconstruction timings as JSON. OBJ coordinates are
+The benchmark always reports the valid-point fraction, depth percentiles
+and matching/reconstruction timings as JSON, and optionally writes a
+colour-mapped depth image (third argument, `-` to skip) for visual
+inspection. The optional depth range applies the same filter as the CLI's
+`cloud_point` section; without it the generic 0.01–10 m defaults are used,
+which lets a few residual mismatches at metres of depth inflate RMSE. When
+`point_cloud.obj` is present it also reports coverage, component-wise and
+3-D errors and threshold accuracy. Note that the `test_dataset_*` archives
+ship without `point_cloud.obj`; ground truth is only in the full
+`dataset_N.zip` archives (13–40 GB each). The zips are served with HTTP
+range support, so single keyframes can be extracted remotely with Python's
+`zipfile` over a seekable HTTP file object instead of downloading the whole
+archive.
+
+`scripts/scared_overview.py` runs the benchmark over many keyframes and
+prints a Markdown table:
+
+```bash
+scripts/scared_overview.py --png-dir out/depth --depth-range 0.02 0.30 \
+    datasets/scared/dataset_1/keyframe_* datasets/scared/test_dataset_8/keyframe_*
+``` OBJ coordinates are
 converted from millimetres to metres and rectified into the same left-camera
 frame as the reconstructed cloud before evaluation. Reference numbers for
 `dataset_1/keyframe_1` with the tuned SGBM configuration: coverage ≈ 0.85,

scripts/scared_overview.py

@@ -0,0 +1,83 @@
+#!/usr/bin/env python3
+"""Run scared_dataset_benchmark over many keyframes and print a Markdown table.
+
+Usage:
+    scripts/scared_overview.py [--bench PATH] [--disparities N] [--png-dir DIR]
+                               [--json-dir DIR] [--depth-range MIN_M MAX_M]
+                               KEYFRAME_DIR...
+
+Each KEYFRAME_DIR must hold Left_Image.png, Right_Image.png and
+endoscope_calibration.yaml. Accuracy columns are filled in only for keyframes
+that also contain point_cloud.obj.
+"""
+import argparse
+import json
+import os
+import subprocess
+import sys
+
+
+def run(bench, kf, disparities, png_dir, json_dir, depth_range):
+    label = "/".join(kf.rstrip("/").split("/")[-2:])
+    cmd = [bench, kf, str(disparities)]
+    if png_dir:
+        os.makedirs(png_dir, exist_ok=True)
+        cmd.append(os.path.join(png_dir, label.replace("/", "_") + ".png"))
+    else:
+        cmd.append("-")
+    if depth_range:
+        cmd += [str(depth_range[0]), str(depth_range[1])]
+    proc = subprocess.run(cmd, capture_output=True, text=True)
+    if proc.returncode != 0:
+        print(f"{label}: benchmark failed\n{proc.stderr}", file=sys.stderr)
+        return label, None
+    result = json.loads(proc.stdout)
+    if json_dir:
+        os.makedirs(json_dir, exist_ok=True)
+        with open(os.path.join(json_dir, label.replace("/", "_") + ".json"), "w") as f:
+            json.dump(result, f, indent=2)
+    return label, result
+
+
+def fmt_mm(v):
+    return "" if v is None else f"{v * 1000:.2f}"
+
+
+def fmt_pct(v):
+    return "" if v is None else f"{v * 100:.1f}"
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--bench", default="build-opencv/src/cloud_point/scared_dataset_benchmark")
+    ap.add_argument("--disparities", type=int, default=160)
+    ap.add_argument("--png-dir")
+    ap.add_argument("--json-dir")
+    ap.add_argument("--depth-range", nargs=2, type=float, metavar=("MIN_M", "MAX_M"),
+                    help="depth filter in metres (default: builder defaults 0.01..10)")
+    ap.add_argument("keyframes", nargs="+")
+    args = ap.parse_args()
+
+    rows = [run(args.bench, kf, args.disparities, args.png_dir, args.json_dir,
+                args.depth_range)
+            for kf in args.keyframes]
+
+    print("| keyframe | valid % | z p5 / median / p95 (mm) | match ms | GT | coverage % | MAE3D mm | RMSE3D mm | median mm | <1 mm % | <2 mm % | <5 mm % |")
+    print("|---|---|---|---|---|---|---|---|---|---|---|---|")
+    for label, r in rows:
+        if r is None:
+            print(f"| {label} | failed | | | | | | | | | | |")
+            continue
+        z = f"{r['z_p05_m']*1000:.0f} / {r['z_median_m']*1000:.0f} / {r['z_p95_m']*1000:.0f}"
+        gt = r.get("has_ground_truth", False)
+        print("| {} | {} | {} | {:.0f} | {} | {} | {} | {} | {} | {} | {} | {} |".format(
+            label, fmt_pct(r["valid_fraction"]), z, r["matching_ms"],
+            "yes" if gt else "no",
+            fmt_pct(r.get("coverage")), fmt_mm(r.get("mae_3d_m")),
+            fmt_mm(r.get("rmse_3d_m")), fmt_mm(r.get("median_3d_m")),
+            fmt_pct(r.get("within_1mm")), fmt_pct(r.get("within_2mm")),
+            fmt_pct(r.get("within_5mm"))))
+
+
+if __name__ == "__main__":
+    main()

src/cloud_point/scared_dataset_benchmark.cpp

@@ -1,13 +1,20 @@
 /// @file scared_dataset_benchmark.cpp
-/// @brief Evaluate the CPU stereo reconstruction against SCARED XYZ ground
-/// truth.
+/// @brief Reconstruct a SCARED keyframe with the CPU stereo pipeline, report
+/// depth statistics and, when point_cloud.obj is present, accuracy against
+/// the XYZ ground truth.
+#include <algorithm>
 #include <chrono>
+#include <cmath>
+#include <filesystem>
 #include <iostream>
+#include <optional>
 #include <stdexcept>
 #include <string>
+#include <vector>
 
 #include <glog/logging.h>
 #include <nlohmann/json.hpp>
+#include <opencv2/imgcodecs.hpp>
 #include <opencv2/imgproc.hpp>
 
 #include "cloud_point/imageFactory.h"
@@ -37,34 +44,114 @@ cv::Mat to_gray(const score::ImageRPC &rpc) {
     return gray;
 }
 
+/// Depth percentiles (metres) over the valid points of a cloud.
+struct DepthStats {
+    std::size_t valid_points{0};
+    double valid_fraction{0.0};
+    double z_min_m{0.0}, z_p05_m{0.0}, z_median_m{0.0}, z_p95_m{0.0},
+        z_max_m{0.0};
+};
+
+DepthStats depth_stats(const score::PointCloud &cloud) {
+    std::vector<float> z;
+    z.reserve(static_cast<std::size_t>(cloud.width) *
+              static_cast<std::size_t>(cloud.height));
+    for (const auto &pt : cloud.valid_points())
+        z.push_back(pt[2]);
+    DepthStats stats;
+    stats.valid_points = z.size();
+    stats.valid_fraction =
+        static_cast<double>(z.size()) /
+        (static_cast<double>(cloud.width) * static_cast<double>(cloud.height));
+    if (z.empty())
+        return stats;
+    std::sort(z.begin(), z.end());
+    const auto at = [&](double q) {
+        return static_cast<double>(
+            z[std::min(z.size() - 1, static_cast<std::size_t>(
+                                         q * static_cast<double>(z.size())))]);
+    };
+    stats.z_min_m    = z.front();
+    stats.z_p05_m    = at(0.05);
+    stats.z_median_m = at(0.5);
+    stats.z_p95_m    = at(0.95);
+    stats.z_max_m    = z.back();
+    return stats;
+}
+
+/// Colour-mapped depth image (invalid pixels black) for visual inspection.
+void write_depth_png(const score::PointCloud &cloud, const std::string &path,
+                     double z_lo, double z_hi) {
+    cv::Mat depth8(cloud.height, cloud.width, CV_8UC1, cv::Scalar(0));
+    cv::Mat valid(cloud.height, cloud.width, CV_8UC1, cv::Scalar(0));
+    const double span = std::max(z_hi - z_lo, 1e-6);
+    for (int r = 0; r < cloud.height; ++r) {
+        for (int c = 0; c < cloud.width; ++c) {
+            const std::size_t idx = (static_cast<std::size_t>(r) *
+                                         static_cast<std::size_t>(cloud.width) +
+                                     static_cast<std::size_t>(c)) *
+                                    3u;
+            const float z         = cloud.data[idx + 2];
+            if (std::isnan(z))
+                continue;
+            const double t         = std::clamp((z - z_lo) / span, 0.0, 1.0);
+            depth8.at<uchar>(r, c) = static_cast<uchar>(255.0 * (1.0 - t));
+            valid.at<uchar>(r, c)  = 255;
+        }
+    }
+    cv::Mat colour;
+    cv::applyColorMap(depth8, colour, cv::COLORMAP_TURBO);
+    colour.setTo(cv::Scalar(0, 0, 0), valid == 0);
+    if (!cv::imwrite(path, colour))
+        throw std::runtime_error("failed to write " + path);
+}
+
 } // namespace
 
 int main(int argc, char *argv[]) {
     google::InitGoogleLogging(argv[0]);
     FLAGS_alsologtostderr = true;
 
-    if (argc < 2 || argc > 3) {
+    if (argc < 2 || argc > 6 || argc == 5) {
         std::cerr << "Usage: " << argv[0]
-                  << " <keyframe_dir> [num_disparities]\n";
+                  << " <keyframe_dir> [num_disparities] [depth_png|-] "
+                     "[min_depth_m max_depth_m]\n"
+                  << "  Accuracy metrics are reported when "
+                     "<keyframe_dir>/point_cloud.obj exists;\n"
+                  << "  depth statistics are always reported.\n";
         return 1;
     }
 
     try {
         const std::string keyframe_dir = argv[1];
-        const int num_disparities      = argc == 3 ? std::stoi(argv[2]) : 160;
+        const int num_disparities      = argc >= 3 ? std::stoi(argv[2]) : 160;
+        std::string depth_png          = argc >= 4 ? argv[3] : "";
+        if (depth_png == "-")
+            depth_png.clear();
+        score::PointCloudBuilder::Options depth_range;
+        if (argc == 6) {
+            depth_range = score::PointCloudBuilder::Options{std::stof(argv[4]),
+                                                            std::stof(argv[5])};
+        }
 
         score::ScaredDatasetLoader dataset(keyframe_dir);
         const auto &calibration_rpc = dataset.calibration();
-        score::ScaredGroundTruthLoader ground_truth(
-            keyframe_dir,
-            cv::Size(calibration_rpc.width, calibration_rpc.height));
+        const cv::Size image_size(calibration_rpc.width,
+                                  calibration_rpc.height);
+        std::optional<score::ScaredGroundTruthLoader> ground_truth;
+        if (std::filesystem::exists(keyframe_dir + "/point_cloud.obj")) {
+            ground_truth.emplace(keyframe_dir, image_size);
+        } else {
+            LOG(WARNING) << "No point_cloud.obj in " << keyframe_dir
+                         << "; reporting depth statistics only";
+        }
 
         const auto calibration =
             score::StereoRectifier::Calibration::from_rpc(calibration_rpc);
         score::StereoRectifier rectifier(calibration);
         auto matcher = score::StereoMatcherFactory::create(
             score::StereoAlgorithmType::CPU, num_disparities);
-        score::PointCloudBuilder builder(rectifier.q());
+        score::PointCloudBuilder builder(rectifier.q(), depth_range);
 
         const auto pair          = dataset.image_pair(0);
         const cv::Mat left_gray  = to_gray(pair.left);
@@ -79,10 +166,10 @@ int main(int argc, char *argv[]) {
         const score::PointCloud cloud = builder.build(disparity);
         const auto reconstruction_end = std::chrono::steady_clock::now();
 
-        const cv::Mat rectified_ground_truth =
-            rectifier.rectify_left_point_map(ground_truth.point_map());
-        const auto metrics =
-            score::PointCloudEvaluator::evaluate(cloud, rectified_ground_truth);
+        const auto stats              = depth_stats(cloud);
+        if (!depth_png.empty()) {
+            write_depth_png(cloud, depth_png, stats.z_p05_m, stats.z_p95_m);
+        }
 
         const auto matching_ms =
             std::chrono::duration<double, std::milli>(matching_end - start)
@@ -92,25 +179,46 @@ int main(int argc, char *argv[]) {
                                                       matching_end)
                 .count();
 
-        const nlohmann::json output = {
+        nlohmann::json output = {
             {"keyframe_dir", keyframe_dir},
             {"algorithm", "StereoSGBM"},
             {"num_disparities", num_disparities},
-            {"ground_truth_points", metrics.ground_truth_points},
-            {"matched_points", metrics.matched_points},
-            {"coverage", metrics.coverage},
-            {"mae_x_m", metrics.mae_x_m},
-            {"mae_y_m", metrics.mae_y_m},
-            {"mae_z_m", metrics.mae_z_m},
-            {"mae_3d_m", metrics.mae_3d_m},
-            {"rmse_3d_m", metrics.rmse_3d_m},
-            {"median_3d_m", metrics.median_3d_m},
-            {"within_1mm", metrics.within_1mm},
-            {"within_2mm", metrics.within_2mm},
-            {"within_5mm", metrics.within_5mm},
+            {"min_depth_m", depth_range.min_depth_m},
+            {"max_depth_m", depth_range.max_depth_m},
+            {"image_width", cloud.width},
+            {"image_height", cloud.height},
+            {"valid_points", stats.valid_points},
+            {"valid_fraction", stats.valid_fraction},
+            {"z_min_m", stats.z_min_m},
+            {"z_p05_m", stats.z_p05_m},
+            {"z_median_m", stats.z_median_m},
+            {"z_p95_m", stats.z_p95_m},
+            {"z_max_m", stats.z_max_m},
             {"matching_ms", matching_ms},
             {"reconstruction_ms", reconstruction_ms},
+            {"has_ground_truth", ground_truth.has_value()},
         };
+        if (!depth_png.empty())
+            output["depth_png"] = depth_png;
+
+        if (ground_truth) {
+            const cv::Mat rectified_ground_truth =
+                rectifier.rectify_left_point_map(ground_truth->point_map());
+            const auto metrics = score::PointCloudEvaluator::evaluate(
+                cloud, rectified_ground_truth);
+            output["ground_truth_points"] = metrics.ground_truth_points;
+            output["matched_points"]      = metrics.matched_points;
+            output["coverage"]            = metrics.coverage;
+            output["mae_x_m"]             = metrics.mae_x_m;
+            output["mae_y_m"]             = metrics.mae_y_m;
+            output["mae_z_m"]             = metrics.mae_z_m;
+            output["mae_3d_m"]            = metrics.mae_3d_m;
+            output["rmse_3d_m"]           = metrics.rmse_3d_m;
+            output["median_3d_m"]         = metrics.median_3d_m;
+            output["within_1mm"]          = metrics.within_1mm;
+            output["within_2mm"]          = metrics.within_2mm;
+            output["within_5mm"]          = metrics.within_5mm;
+        }
         std::cout << output.dump(2) << '\n';
     } catch (const std::exception &error) {
         std::cerr << "Benchmark failed: " << error.what() << '\n';