diff --git a/README.md b/README.md index 736756b..bb5200a 100644 --- a/README.md +++ b/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, diff --git a/scripts/scared_overview.py b/scripts/scared_overview.py new file mode 100755 index 0000000..509c9d5 --- /dev/null +++ b/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() diff --git a/src/cloud_point/scared_dataset_benchmark.cpp b/src/cloud_point/scared_dataset_benchmark.cpp index 5c2d3ba..81d61af 100644 --- a/src/cloud_point/scared_dataset_benchmark.cpp +++ b/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 #include +#include +#include #include +#include #include #include +#include #include #include +#include #include #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 z; + z.reserve(static_cast(cloud.width) * + static_cast(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(z.size()) / + (static_cast(cloud.width) * static_cast(cloud.height)); + if (z.empty()) + return stats; + std::sort(z.begin(), z.end()); + const auto at = [&](double q) { + return static_cast( + z[std::min(z.size() - 1, static_cast( + q * static_cast(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(r) * + static_cast(cloud.width) + + static_cast(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(r, c) = static_cast(255.0 * (1.0 - t)); + valid.at(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] - << " [num_disparities]\n"; + << " [num_disparities] [depth_png|-] " + "[min_depth_m max_depth_m]\n" + << " Accuracy metrics are reported when " + "/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 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(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';