# -*- coding: utf-8 -*- """ Evaluate a reconstructed dense point cloud against GT. - Reads ASCII / binary PLY point clouds (x,y,z + optional nx,ny,nz + red,green,blue). - Aligns the reconstruction into the GT frame using: (a) similarity transform from camera centres (COLMAP poses -> GT transforms.json), then (b) point-to-point ICP refinement. - Reports Chamfer distance, Precision / Recall / F1 at several thresholds, and saves a coloured 3D scatter rendering of the reconstruction. Usage: python eval_pointclouds.py --recon output/colmap/dense/fused.ply \ --gt mvs_project/gt/dense_pointcloud.ply \ --poses output/colmap/cameras_export/poses_c2w.npz \ --gt_transforms mvs_project/poses/transforms.json \ --outdir output/eval --label mvs """ import os, sys, argparse, json import numpy as np from scipy.spatial import cKDTree # ---------------------------------------------------------------- PLY reader def read_ply(path, max_points=None): import io with open(path, "rb") as f: raw = f.read() hdr_end = raw.find(b"end_header") hdr_end = raw.find(b"\n", hdr_end) + 1 header = raw[:hdr_end].decode("ascii", "ignore") body = raw[hdr_end:] fmt = "ascii" props = [] # (name, type) n_verts = 0 for line in header.splitlines(): line = line.strip() if line.startswith("format"): fmt = line.split()[1] elif line.startswith("element vertex"): n_verts = int(line.split()[2]) elif line.startswith("property"): parts = line.split() props.append((parts[2], parts[1])) dtype_map = {"float": "f4", "float32": "f4", "double": "f8", "int": "i4", "int32": "i4", "uchar": "u1", "uint8": "u1", "short": "i2", "ushort": "u2", "uint": "u4"} idx = {p[0]: i for i, p in enumerate(props)} if fmt == "ascii": if max_points is not None and max_points < n_verts: n_verts = max_points data = np.loadtxt(io.BytesIO(body), max_rows=n_verts) if data.ndim == 1: data = data[None, :] else: endian = "<" if "little" in fmt else ">" dt = np.dtype([(p[0], endian + dtype_map[p[1]]) for p in props]) verts = np.frombuffer(body, dtype=dt, count=n_verts) data = np.zeros((len(verts), len(props)), dtype=np.float32) for i, (name, _) in enumerate(props): data[:, i] = verts[name] xi, yi, zi = idx["x"], idx["y"], idx["z"] pts = data[:, [xi, yi, zi]].astype(np.float64) col = None if "red" in idx and "green" in idx and "blue" in idx: col = data[:, [idx["red"], idx["green"], idx["blue"]]].astype(np.uint8) return pts, col # ---------------------------------------------------------------- transforms def umeyama(src, dst): mu_s = src.mean(0); mu_d = dst.mean(0) s = src - mu_s; d = dst - mu_d cov = d.T @ s / s.shape[0] U, S, Vt = np.linalg.svd(cov) R = U @ Vt if np.linalg.det(R) < 0: U[:, -1] *= -1; R = U @ Vt scale = np.trace(np.diag(S)) / ((s ** 2).sum() / s.shape[0]) t = mu_d - scale * R @ mu_s return scale, R, t def apply_sim(P, s, R, t): return (s * R @ P.T).T + t def icp(src, dst, iters=40, subsample=50000): """point-to-point ICP, aligns src -> dst.""" cur = src.copy() for it in range(iters): if len(cur) > subsample: sel = np.random.choice(len(cur), subsample, replace=False) cur_s = cur[sel] else: cur_s = cur tree = cKDTree(dst) dd, ii = tree.query(cur_s, k=1) corr = dst[ii] # filter correspondences that are too far (outliers) thresh = np.quantile(dd, 0.9) mask = dd < max(thresh, 1.0) if mask.sum() < 50: break s, R, t = umeyama(cur_s[mask], corr[mask]) cur = apply_sim(cur, s, R, t) if it % 10 == 0: print(f" ICP iter {it}: mean dist {dd[mask].mean():.4f}") return cur def chamfer(P, Q, subsample=100000): """bidirectional chamfer + precision/recall/F1 at thresholds.""" if len(P) > subsample: P = P[np.random.choice(len(P), subsample, replace=False)] if len(Q) > subsample: Q = Q[np.random.choice(len(Q), subsample, replace=False)] tq = cKDTree(Q) tp = cKDTree(P) d_PQ, _ = tq.query(P, k=1) # each recon point -> GT d_QP, _ = tp.query(Q, k=1) # each GT point -> recon cd = 0.5 * (d_PQ.mean() + d_QP.mean()) res = {"chamfer": cd, "mean_recon_to_gt": d_PQ.mean(), "mean_gt_to_recon": d_QP.mean(), "median_recon_to_gt": np.median(d_PQ), "median_gt_to_recon": np.median(d_QP)} for th in (0.01, 0.02, 0.05, 0.10): prec = (d_PQ < th).mean() # fraction of recon points within th of GT rec = (d_QP < th).mean() # fraction of GT points within th of recon f1 = 2 * prec * rec / (prec + rec + 1e-12) res[f"th{int(th*100)}cm"] = {"precision": float(prec), "recall": float(rec), "f1": float(f1)} return res def render_scatter(pts, col, out_png, title, azim=60, elev=25): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(111, projection="3d") if len(pts) > 80000: sel = np.random.choice(len(pts), 80000, replace=False) pts, col = pts[sel], (col[sel] if col is not None else None) c = (col / 255.0) if col is not None else (pts[:, 2] - pts[:, 2].min()) ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], c=c, s=0.5, marker=".") ax.set_title(title, fontsize=12) ax.set_xlabel("x"); ax.set_ylabel("y"); ax.set_zlabel("z") ax.view_init(elev=elev, azim=azim) plt.tight_layout() plt.savefig(out_png, dpi=130) plt.close(fig) def main(): ap = argparse.ArgumentParser() ap.add_argument("--recon", required=True) ap.add_argument("--gt", required=True) ap.add_argument("--poses", default=None) ap.add_argument("--gt_transforms", default=None) ap.add_argument("--outdir", required=True) ap.add_argument("--label", default="recon") ap.add_argument("--no_icp", action="store_true") args = ap.parse_args() os.makedirs(args.outdir, exist_ok=True) print(f"loading recon: {args.recon}") P, Pc = read_ply(args.recon) print(f" recon points: {len(P)}") print(f"loading GT: {args.gt}") Q, Qc = read_ply(args.gt) print(f" GT points: {len(Q)}") # initial similarity alignment from camera centres (if provided) if args.poses and args.gt_transforms: d = np.load(args.poses, allow_pickle=True) names, c2w = d["names"], d["c2w"] gt_tr = json.load(open(args.gt_transforms)) gmap = {f["file_path"].split('/')[-1]: np.array(f["transform_matrix"]) for f in gt_tr["frames"]} A, B = [], [] for n, c in zip(names, c2w): if n in gmap: A.append(c[:3, 3]); B.append(gmap[n][:3, 3]) A, B = np.array(A), np.array(B) s, R, t = umeyama(A, B) P = apply_sim(P, s, R, t) print(f"camera-centre similarity: s={s:.4f} applied") if not args.no_icp: print("running ICP refinement ...") P = icp(P, Q) res = chamfer(P, Q) print("\n===== EVALUATION (%s) =====" % args.label) print(f"Chamfer distance: {res['chamfer']*100:.3f} cm") print(f"mean recon->GT: {res['mean_recon_to_gt']*100:.3f} cm | median {res['median_recon_to_gt']*100:.3f} cm") print(f"mean GT->recon: {res['mean_gt_to_recon']*100:.3f} cm | median {res['median_gt_to_recon']*100:.3f} cm") for k in ("th1cm", "th2cm", "th5cm", "th10cm"): r = res[k] print(f" @{k[2:]} : precision {r['precision']:.3f} | recall {r['recall']:.3f} | F1 {r['f1']:.3f}") # save aligned recon cloud out_ply = os.path.join(args.outdir, f"{args.label}_aligned.ply") write_ply(out_ply, P, Pc) print("saved aligned cloud:", out_ply) png = os.path.join(args.outdir, f"{args.label}_aligned.png") render_scatter(P, Pc, png, f"{args.label}: {len(P)} points") print("saved render:", png) with open(os.path.join(args.outdir, f"{args.label}_metrics.json"), "w") as f: json.dump(res, f, indent=2) print("saved metrics:", os.path.join(args.outdir, f"{args.label}_metrics.json")) def write_ply(path, pts, col): with open(path, "w") as f: f.write("ply\nformat ascii 1.0\n") f.write(f"element vertex {len(pts)}\n") f.write("property float x\nproperty float y\nproperty float z\n") f.write("property uchar red\nproperty uchar green\nproperty uchar blue\n") f.write("end_header\n") for i in range(len(pts)): r, g, b = (col[i] if col is not None else [128, 128, 128]) f.write(f"{pts[i,0]:.6f} {pts[i,1]:.6f} {pts[i,2]:.6f} {int(r)} {int(g)} {int(b)}\n") if __name__ == "__main__": main()