1 def _compute_density_stats(values: np.ndarray) -> dict:
 2     """Compute summary statistics for a density array."""
 3     return {
 4         "mean": float(np.mean(values)),
 5         "std": float(np.std(values)),
 6         "min": float(np.min(values)),
 7         "max": float(np.max(values)),
 8         "median": float(np.median(values)),
 9         "percentile_5": float(np.percentile(values, 5)),
10         "percentile_25": float(np.percentile(values, 25)),
11         "percentile_75": float(np.percentile(values, 75)),
12         "percentile_95": float(np.percentile(values, 95)),
13     }
14 
15 
16 def _density_knn(pcd: o3d.geometry.PointCloud, k: int) -> np.ndarray:
17     """KNN-based density: k / volume_of_sphere(r_k)."""
18     kdtree = o3d.geometry.KDTreeFlann(pcd)
19     points = np.asarray(pcd.points)
20     n = len(points)
21     densities = np.zeros(n, dtype=np.float64)
22 
23     for i in range(n):
24         _, idx, dists_sq = kdtree.search_knn_vector_3d(pcd.points[i], k + 1)
25         neighbor_dists = np.sqrt(np.array(dists_sq[1:], dtype=np.float64))
26         r_k = neighbor_dists[-1]
27         if r_k > 1e-10:
28             volume = (4.0 / 3.0) * np.pi * (r_k**3)
29             densities[i] = k / volume
30 
31     return densities
32 
33 
34 def _density_radius(pcd: o3d.geometry.PointCloud, radius: float) -> np.ndarray:
35     """Radius-based density: count of neighbors within fixed radius."""
36     kdtree = o3d.geometry.KDTreeFlann(pcd)
37     n = len(pcd.points)
38     densities = np.zeros(n, dtype=np.float64)
39 
40     for i in range(n):
41         count, _, _ = kdtree.search_radius_vector_3d(pcd.points[i], radius)
42         densities[i] = count - 1
43 
44     return densities
45 
46 
47 def _density_voxel(pcd: o3d.geometry.PointCloud, voxel_size: float) -> np.ndarray:
48     """Voxel-based density: points in enclosing voxel / voxel_volume."""
49     points = np.asarray(pcd.points)
50     voxel_indices = np.floor(points / voxel_size).astype(np.int64)
51 
52     voxel_keys = (
53         voxel_indices[:, 0].astype(np.int64) * 73856093
54         ^ voxel_indices[:, 1].astype(np.int64) * 19349663
55         ^ voxel_indices[:, 2].astype(np.int64) * 83492791
56     )
57 
58     counts = Counter(voxel_keys.tolist())
59     voxel_volume = voxel_size**3
60     densities = np.array([counts[k] / voxel_volume for k in voxel_keys], dtype=np.float64)
61     return densities