matlab alpha-shape平面点云边界提取

clc;clear; %% ------------------------------读取点云-------------------------------- ptCloud = pcread('E://data//proo.pcd'); figure; pcshow(ptCloud); title('原始点云'); hold on; loc = ptCloud.Location; P = [loc(:,1),loc(:,2),loc(:,3)]; %% ------------------------------参数设置---------------------------------- alpha = 0.2; % %% ------------------------alpha-shape实现过程----------------------------- boundaryIdx =[]; for i = 1:length(loc) point_p = P(i,:); % 利用KD树根据搜索半径查找每个点的邻域 [indices,dist] = findNeighborsInRadius(ptCloud,point_p,2*alpha,'sort',true); % 查找到每个点的邻域，构造边界 for j = 2:length(indices) point_q = P(indices(j),:); x = point_p(1); y = point_p(2); x1 = point_q(1); y1 = point_q(2); % 计算s的平方 s_Square = dot(point_p(1:2) - point_q(1:2),point_p(1:2) - point_q(1:2)); % 计算H H = sqrt((power(alpha,2) / s_Square) - 0.25); % 计算过点p p1的两个圆的圆心位置 x2 = x + 0.5 * (x1 - x) - H * (y1 - y); y2 = y + 0.5 * (y1 - y) - H * (x - x1); x3 = x + 0.5 * (x1 - x) + H * (y1 - y); y3 = y + 0.5 * (y1 - y) + H * (x - x1); tem_dis_one = 1000000; tem_dis_two = 1000000; loop_num = 0; % 搜索剩余邻域点，计算剩余邻域点到过p,p1两点的圆的圆心的距离，判断是否大于搜索半径，是否为边界点 for k = 2:length(indices) if k == j continue; % 去除p p1两点 end % 计算邻域中剩余点到两个圆心的距离 dis_one = sqrt(power((P(indices(k),1) - x2), 2) + power((P(indices(k),2) - y2), 2)); dis_two = sqrt(power((P(indices(k),1) - x3), 2) + power((P(indices(k),2) - y3), 2)); if dis_one < alpha && dis_two < alpha break; end if tem_dis_one > dis_one tem_dis_one = dis_one; end if tem_dis_two > dis_two tem_dis_two = dis_two; end loop_num=loop_num+1; if loop_num == length(indices) - 2 if tem_dis_one > alpha || tem_dis_two > alpha boundaryIdx = [boundaryIdx;i]; % 获取边界点索引 break; end end end end end ptCloudB = select(ptCloud,boundaryIdx);% 根据索引提取边界点 pcwrite(ptCloudB,'boundaryPointCloud.pcd','Encoding','binary'); figure; pcshow(ptCloudB); title('边界点云');