Python实现区域生长和RANSAC聚类

        本节我们分享区域生长和RANSAC进行点云聚类的方法，相比于上节，本节更专注同一块点云之间的聚类。

区域生长和RANSAC（Random Sample Consensus）是两种最常用的点云聚类思路，前者“由点到面、逐步扩张”，后者“大胆假设、小心验证”。下面用更简洁、连贯的语言重新梳理两者的原理、步骤与应用场景。

------------------------------------------------

1、区域生长：让“种子点”发芽成区域  
1）核心思想  
选取一个或几个种子点，按照“相邻且相似”的规则，像水波纹一样向外扩张，直到再也找不到符合规则的邻居为止；随后从剩余点再选新种子，重复此过程，直至所有点被归类。

        2）算法流程  
(1) 种子选取：可随机，也可优先选曲率最小的点（更平坦、更稳定）。  
(2)区域扩张：  
‑ 对当前种子，检查其 k 邻域。  
‑ 若某邻点同时满足以下三项条件，则并入当前区域并成为新种子：  
‑ 距离条件：邻点与种子点的欧氏距离 < d_th  
‑ 法向夹角：两法向量夹角 < θ_th  
‑ 曲率差：|邻点曲率 – 种子点曲率| < c_th  
(3)重复直至无新点可加入，标记该区域为完成。  
(4)从剩余未访问点重新选种子，开始新一轮生长，直至全部分割完毕。

        3) 典型应用  
• 复杂物体按曲面片/平面片进行精细分割  
• 场景理解中把不同材质或不同几何属性的表面区分开  
• 与 RANSAC 互补：RANSAC 先粗提平面，区域生长再细分细节

2、RANSAC：用“少数服从多数”的思路找模型  
1) 核心思想  
在含有大量离群点的数据中，通过“随机采样+投票验证”的方式，反复试错，最终选出一个被最多数据点支持的模型。

        2) 算法流程  
(1) 设定：最大迭代次数 max_iter、内点距离阈值 d_th。  
(2) 随机采样：从全部 N 个点里抽取“恰好能确定模型”的最小样本（例如平面只需 3 点）。  
(3) 计算模型：用这 3 点拟合平面方程。  
(4) 投票检验：计算所有点到该平面的距离，距离 ≤ d_th 的点记为内点。  
(5) 更新最佳：若当前内点数大于历史最佳，则保存此模型及其内点。  
(6) 重复迭代：直到 max_iter 次或内点比例已足够高。

        3) 典型应用  
• 点云平面分割（地面、墙壁）。  
• 图像拼接中的特征匹配与单应矩阵估计。 
• 3D 物体识别中的基本形状提取。

一句话总结  
RANSAC 擅长“大海捞针”式地找出一条直线或一个平面；区域生长则擅长“顺藤摸瓜”式地把相连且相似的点聚成一片。两者常结合使用：先用 RANSAC 快速提取主要平面，再用区域生长对剩余点进行精细化分割。

        本期我们的数据依然，哦不是，使用的是一个长方体点云，生成部分在程序里，第一个显示的就是。哈哈哈

一、区域生长和RANSAC聚类程序

import open3d as o3d
import numpy as np
from collections import deque
from typing import List, Tuple# -----------------------------------------------------------------------------
# 1. 区域生长分割
# -----------------------------------------------------------------------------
class RegionGrowing:"""区域生长分割（Region Growing Segmentation）仅依赖 Open3D 与 Numpy，无外部依赖。"""def __init__(self,cloud: o3d.geometry.PointCloud,*,min_pts_per_cluster: int = 1,max_pts_per_cluster: int = int(1e6),angle_threshold_deg: float = 30.0,curvature_threshold: float = 0.05,k_neighbours: int = 30):self.pcd = cloudself.min_pts = min_pts_per_clusterself.max_pts = max_pts_per_clusterself.cos_thresh = np.cos(np.deg2rad(angle_threshold_deg))self.curv_thresh = curvature_thresholdself.k = k_neighboursself.curvature: np.ndarray | None = Noneself.labels: np.ndarray | None = Noneself.clusters: List[List[int]] = []# ---------------------------------------------------------------------def _pre_check(self) -> bool:if not self.pcd.has_points():return Falseif self.k <= 0:return Falseif not self.pcd.has_normals():self.pcd.estimate_normals(o3d.geometry.KDTreeSearchParamKNN(self.k))return True# ---------------------------------------------------------------------def _build_kdtree(self) -> np.ndarray:kdtree = o3d.geometry.KDTreeFlann(self.pcd)n_pts = len(self.pcd.points)neighbours = np.empty((n_pts, self.k), dtype=int)for i in range(n_pts):_, idx, _ = kdtree.search_knn_vector_3d(self.pcd.points[i], self.k)neighbours[i] = idxreturn neighbours# ---------------------------------------------------------------------def _compute_curvature(self) -> np.ndarray:self.pcd.estimate_covariances(o3d.geometry.KDTreeSearchParamKNN(self.k))covs = np.asarray(self.pcd.covariances)curv = np.empty(len(covs))for i, cov in enumerate(covs):eigval = np.linalg.eigvals(cov)eigval.sort()curv[i] = eigval[0] / eigval.sum()return curv# ---------------------------------------------------------------------def _validate_point(self, seed: int, cand: int) -> Tuple[bool, bool]:normal_seed = self.pcd.normals[seed]normal_cand = self.pcd.normals[cand]if np.fabs(np.dot(normal_seed, normal_cand)) < self.cos_thresh:return False, Falseis_seed = self.curvature[cand] <= self.curv_threshreturn True, is_seed# ---------------------------------------------------------------------def _grow_region(self, seed: int, seg_id: int) -> int:seeds = deque([seed])self.labels[seed] = seg_idnum = 1while seeds:current = seeds.popleft()for nb in self.neighbours[current]:if self.labels[nb] != -1:continueok, add_seed = self._validate_point(current, nb)if not ok:continueself.labels[nb] = seg_idnum += 1if add_seed:seeds.append(nb)return num# ---------------------------------------------------------------------def extract(self) -> List[List[int]]:if not self._pre_check():raise RuntimeError("RegionGrowing pre-check failed.")self.curvature = self._compute_curvature()self.neighbours = self._build_kdtree()n_pts = len(self.pcd.points)self.labels = -np.ones(n_pts, dtype=int)order = np.argsort(self.curvature)clusters: List[List[int]] = []seg_id = 0for idx in order:if self.labels[idx] != -1:continuecnt = self._grow_region(idx, seg_id)if self.min_pts <= cnt <= self.max_pts:clusters.append(np.where(self.labels == seg_id)[0].tolist())else:self.labels[self.labels == seg_id] = -1seg_id += 1self.clusters = clustersreturn clusters# -----------------------------------------------------------------------------
# 2. RANSAC 多平面分割
# -----------------------------------------------------------------------------
class MultiPlaneRANSAC:"""使用 RANSAC 连续提取多个平面"""def __init__(self,cloud: o3d.geometry.PointCloud,*,distance_threshold: float = 0.01,ransac_n: int = 3,max_iter: int = 1000):self.pcd = cloudself.dist_thresh = distance_thresholdself.ransac_n = ransac_nself.max_iter = max_iter# ---------------------------------------------------------------------def extract(self, max_planes: int = 10):planes, inlier_clouds = [], []cloud = self.pcd  # 保留原始点云for _ in range(max_planes):# 关键：提前退出if len(cloud.points) < self.ransac_n:breakplane_model, idx = cloud.segment_plane(distance_threshold=self.dist_thresh,ransac_n=self.ransac_n,num_iterations=self.max_iter)if len(idx) < self.ransac_n:  # 保险起见再判一次breakinlier = cloud.select_by_index(idx)planes.append(np.asarray(plane_model))inlier_clouds.append(inlier)# 从剩余点云中剔除已提取的平面cloud = cloud.select_by_index(idx, invert=True)return planes, inlier_clouds, cloud# 3. 长方体点云生成
# -----------------------------------------------------------------------------
def create_box_pcd(center=(0, 0, 0),extent=(5, 2.5, 2.5),resolution=0.02,noise=0.002) -> o3d.geometry.PointCloud:"""在 center 附近生成一个 extent 大小的长方体点云（表面采样）"""cx, cy, cz = centerdx, dy, dz = extentxs = np.linspace(cx - dx / 2, cx + dx / 2, int(dx / resolution))ys = np.linspace(cy - dy / 2, cy + dy / 2, int(dy / resolution))zs = np.linspace(cz - dz / 2, cz + dz / 2, int(dz / resolution))# 六个面采样points = []for x in xs:for y in ys:points.append([x, y, cz - dz / 2])points.append([x, y, cz + dz / 2])for x in xs:for z in zs:points.append([x, cy - dy / 2, z])points.append([x, cy + dy / 2, z])for y in ys:for z in zs:points.append([cx - dx / 2, y, z])points.append([cx + dx / 2, y, z])points = np.array(points)points += np.random.normal(0, noise, points.shape)pcd = o3d.geometry.PointCloud(o3d.utility.Vector3dVector(points))pcd.estimate_normals()return pcd# -----------------------------------------------------------------------------
# 4. 可视化辅助
# -----------------------------------------------------------------------------
def colorize_clouds(clouds: List[o3d.geometry.PointCloud]) -> None:for c in clouds:c.paint_uniform_color(np.random.rand(3))# -----------------------------------------------------------------------------
# 5. main
# -----------------------------------------------------------------------------
def main():# ----------------- 读取点云 -----------------pcd = create_box_pcd()o3d.visualization.draw_geometries([pcd],window_name="原始点云",width=1200, height=800,left=50, top=50)# ----------------- 区域生长 -----------------rg = RegionGrowing(pcd,min_pts_per_cluster=50,max_pts_per_cluster=int(1e6),angle_threshold_deg=30,curvature_threshold=0.005,k_neighbours=30)clusters = rg.extract()print("区域生长聚类:", len(clusters))rg_clouds = [pcd.select_by_index(c) for c in clusters]colorize_clouds(rg_clouds)o3d.visualization.draw_geometries(rg_clouds,window_name="区域生长聚类",width=1200, height=800,left=50, top=50)# ----------------- RANSAC 多平面 -----------------ransac = MultiPlaneRANSAC(pcd,distance_threshold=0.01,ransac_n=3,max_iter=100000)planes, plane_clouds, rest = ransac.extract(max_planes=10)print("RANSAC聚类:", len(planes))colorize_clouds(plane_clouds)rest.paint_uniform_color([0.5, 0.5, 0.5])# o3d.visualization.draw_geometries(plane_clouds + [rest], # 这块也显示剩下的部分#                                   window_name="区域生长聚类",#                                       width=1200, height=800,#                                       left=50, top=50)o3d.visualization.draw_geometries(plane_clouds,window_name="区域生长聚类",width=1200, height=800,left=50, top=50)if __name__ == "__main__":main()

二、区域生长和RANSAC聚类结果

        可以看到，第二个RANSAC聚类的是非常的快(其实整体加速了，CSDN上面限制大小了，太大的动图上传不了)，这也说明如果能用RANSAC聚类，就别用别的了，man。同时这里注意到，本次使用的数据是一个整体，不是上一节的好几个离的老远的点云。其实更严格的来说，可以叫做分割。不过，随便吧，一个名字鹅以。
就酱，下次见^-^