"""Inside-outside point queries: FCP straight from the descriptor (ring tables at the layer centres with y-bins, CPU numba and CUDA) against standard: trimesh.contains (Embree ray parity), Open3D RaycastingScene.compute_occupancy (Embree ray parity) SOTA: libigl fast_winding_number (Barill et al. 2018, CPU), NVIDIA Warp BVH closest point + winding number (GPU) on the original meshes. Two point sets per object: uniform in the box padded by 10 %, and near the surface (within 2 % of the box diagonal). Accuracy against the original mesh (Open3D occupancy with 3-ray voting) and against the decoded descriptor (fcp.decode mesh). Writes ../report/points.csv. env -u PYTHONPATH .venv/bin/python scripts/points_bench.py """ import csv import sys import time import numpy as np sys.path.insert(0, str(__import__("pathlib").Path(__file__).parent)) from common import NT, OBJ, RULE, REPORT, TAG, cuda_sync, descriptor, fcp, load_mesh, timeit # noqa: E402 import cupy as cp # noqa: E402 import igl # noqa: E402 import open3d as o3d # noqa: E402 import warp as wp # noqa: E402 from fcpq import FCPQuery # noqa: E402 wp.config.quiet = True wp.init() SIZES = [10_000, 1_000_000] @wp.kernel def _wp_inside(mesh: wp.uint64, P: wp.array(dtype=wp.vec3), out: wp.array(dtype=wp.int32)): i = wp.tid() q = wp.mesh_query_point_sign_winding_number(mesh, P[i], 1.0e6, 2.0, 0.5) out[i] = wp.where(q.result and q.sign < 0.0, 1, 0) def points(m, n, kind, rng): lo, hi = m.bounds if kind == "uniform": pad = 0.1 * (hi - lo) return rng.uniform(lo - pad, hi + pad, (n, 3)) diag = np.linalg.norm(hi - lo) S, fi = m.sample(n, return_index=True) return S + m.face_normals[fi] * rng.uniform(-0.02, 0.02, (n, 1)) * diag def main(): rows = [] for name in OBJ: m = load_mesh(name) d = descriptor(name) dm = fcp.decode(d, kind="mesh", resolution=min(d.layers, 256)) # the poly mode (exact polynomial per query) is # left out: 50-200x slower than the table at N 16 (report/a2048), and its CUDA kernel's shared Fourier table does not fit at N 64 V, Fc = np.asarray(m.vertices, np.float64), np.asarray(m.faces, np.int64) # builds FCPQuery(d, nt=NT).gpu(); cuda_sync() # warm-up (numba compile, CuPy context) outside the timed build t0 = time.perf_counter(); q = FCPQuery(d, nt=NT); q.gpu(); cuda_sync(); b_fcp = time.perf_counter() - t0 t0 = time.perf_counter(); scene = o3d.t.geometry.RaycastingScene(); scene.add_triangles(o3d.core.Tensor(V.astype(np.float32)), o3d.core.Tensor(Fc.astype(np.uint32))); b_o3d = time.perf_counter() - t0 # libigl's mesh overload builds its tree inside every call: its "build" is a one-point call (also in query_ms) b_igl = timeit(lambda: igl.fast_winding_number(V, Fc, V[:1]), runs=3) t0 = time.perf_counter(); wmesh = wp.Mesh(points=wp.array(V.astype(np.float32), dtype=wp.vec3), indices=wp.array(Fc.astype(np.int32).ravel(), dtype=wp.int32)); wp.synchronize(); b_wp = time.perf_counter() - t0 t0 = time.perf_counter(); m2 = m.copy(); _ = m2.ray.intersects_any(np.zeros((1, 3)), np.ones((1, 3))); b_tm = time.perf_counter() - t0 rng = np.random.default_rng(0) for kind in ("uniform", "near"): for n in SIZES: P = points(m, n, kind, rng) P32 = P.astype(np.float32) gt = scene.compute_occupancy(o3d.core.Tensor(P32), nsamples=3).numpy().astype(bool) gtd = dm.contains(P) if n <= 100_000 else None Pc = d.pose.to_canonical(P).astype(np.float32) Pc_g = cp.asarray(Pc) P_wp = wp.array(P32, dtype=wp.vec3) out_wp = wp.zeros(n, dtype=wp.int32) runs = 3 if n >= 1_000_000 else 7 methods = { "FCP table (CUDA)": (b_fcp, lambda: q.inside_gpu(Pc_g, RULE), cuda_sync, "gpu"), "FCP table (CPU)": (b_fcp, lambda: q.inside_cpu(Pc, RULE), None, "cpu"), "trimesh contains (Embree)": (b_tm, lambda: m2.contains(P), None, "cpu"), "Open3D occupancy (Embree)": (b_o3d, lambda: scene.compute_occupancy(o3d.core.Tensor(P32)), None, "cpu"), "libigl fast winding number": (b_igl, (lambda: igl.fast_winding_number(V, Fc, P)), None, "cpu"), "Warp BVH winding number (CUDA)": (b_wp, lambda: wp.launch(_wp_inside, dim=n, inputs=[wmesh.id, P_wp, out_wp]), wp.synchronize, "gpu"), } res = {} for meth, (build, fn, sync, dev) in methods.items(): if "trimesh" in meth and n > 100_000: continue # ~ minutes at 1M; timed at 10k sec = timeit(fn, runs=runs, sync=sync) r = fn() if sync: sync() if meth.startswith("FCP"): ins = r.get() if dev == "gpu" else r elif "Warp" in meth: ins = out_wp.numpy().astype(bool) elif "libigl" in meth: ins = np.asarray(r) > 0.5 elif "Open3D" in meth: ins = r.numpy().astype(bool) else: ins = np.asarray(r) res[meth] = ins row = dict(object=name, faces=len(Fc), parts=q.P, points=kind, n=n, method=meth, device=dev, build_ms=build * 1e3, query_ms=sec * 1e3, ns_per_point=sec / n * 1e9, err_mesh=float((ins != gt).mean()), err_decoded=float((ins != gtd).mean()) if gtd is not None else "") rows.append(row) print(name, kind, n, f"{meth:32s} {sec*1e3:9.2f} ms {sec/n*1e9:8.1f} ns/pt err {row['err_mesh']:.4f}", flush=True) if n == SIZES[-1] and kind == "near": np.savez_compressed(REPORT / f"pts{TAG}_{name}.npz", P=P[:60000], gt=gt[:60000], fcp=res["FCP table (CUDA)"][:60000]) with (REPORT / f"points{TAG}.csv").open("w", newline="") as f: w = csv.DictWriter(f, list(rows[0])) w.writeheader() w.writerows(rows) if __name__ == "__main__": main()