一种快速得到点云近似边界球的方法

最新推荐文章于 2025-03-13 16:29:02 发布

拉风小宇

最新推荐文章于 2025-03-13 16:29:02 发布

阅读量1.6k

点赞数 1

分类专栏： MATLAB CAGD/CG 文章标签：网格

本文链接：https://blog.csdn.net/lafengxiaoyu/article/details/73823595

版权

MATLAB 同时被 2 个专栏收录

26 篇文章

订阅专栏

CAGD/CG

24 篇文章

订阅专栏

算法来自1990年J Ritter发表的一篇文章：An efficient bounding sphere

方法简介

这个方法的目的是对于3维空间中的有N个点的集合找到一个接近最优边界球。它是Order（N），并且非常快。该计算的球体比理想的最小半径大约5％。

算法步骤

该算法在执行两遍：第一遍找到两个接近最大间隔的点。这对点作为初始球体。第二遍将每个点与当前球体进行比较，如果点在外面，则放大球体。算法如下：

1.快速遍历一下所有的N个点并且找出如下六点：
最小x的点，最大x的点
最小y的点，最大y的点
最小z的点，最大z的点
这给了三对点。每对都有其最大跨度尺寸。选择具有最大跨度的一对点（其可以大于最大尺寸跨度）。计算初始球体，使用这对点作为直径。

2.再次遍历N个点：对于当前球体的每个外面的点，将当前球体更新为较大的球体，以便通过这个点，以及旧球面的背面。每个新的球体将（几乎就是）旧的球体，加上新的点。
所需的更新将是N的一小部分。对于当前球体在测试每个点时，将点与当前球体的中心距离进行平方与当前球体的半径的平方进行比较，以避免做一个开根号(sqrt)计算。

伪代码

以下伪代码将点（x,y,z）与当前球体进行比较[center =（cenx,ceny,cenz）,radius = rad ]。如果（x,y,z）是在当前球体之外，（cenx,ceny,cenz）和rad更新到反映新的球。当前球的半径平方保持在rad_sq:

given x, y, z, cenx, ceny, cenz,rad, and rad_sq

dx ← x – cenx;

dy ← y – ceny;

dz ← z – cenz;

old_to_p_sq ← dx*dx + dy*dy + dz*dz;

%do economical r**2 test before calc sqrt

if (old_to_p_sq > rad_sq)

then

%Point is outsidecurrent sphere. update.

old_to_p ← old_ to_ p_ sq ;

rad ← (rad + old_to_p)/2.0;

%update square ofradius for next compare

rad_sq ← rad*rad;

old_to_new ← old_to_p–rad;

cenx ← (rad*cenx + old_to_new*x)/old_to_p;

ceny ← (rad*ceny + old_to_new*y) / old_to_p;

cenz ← (rad*cenz + old_to_new*z) / old_to_p;

end

测试例子

以下两个测试在Sun 3/50工作站（68020，具有MC68881浮点数协处理器）。
举例1：
在一个以原点为中心，128为半径的球形体积中随机填满10,000个点。其中五个放置在球体的边缘。这意味着最佳球体应该以128为半径。结果：中心=（3,5,4）; 半径=133（4％>理想）;处理时间：1.8秒。

举例2
在一个以原点为中心，128为半边长度的正方体体积中随机填满10,000个点分，包括其八个角点。意即理想半径sqrt(3)*128=222。注意：这是接近最糟糕的情况，在这种算法下，因为正交对齐的框意味着在初始猜测阶段不会发现任何角落。（一个正方体旋转沿着任何轴线旋转任何角度都将允许最初找到角落）。
结果：center =（5,21,2）; 半径= 237（7％>理想）; 处理时间：1.8秒。

MATLAB代码

在MathWorks上有Anton Semechko分享的代码Exact minimum bounding spheres/circles，其中ApproxMinBoundSphereND.m就是

function [R,C]=ApproxMinBoundSphereND(X)
% Find a near-optimal bounding sphere for a set of M points in 
% N-dimensional space using Ritter's algorithm. In 3D, the resulting sphere
% is approximately 5% bigger than the ideal minimum radius sphere.
%
%   - X     : M-by-N array of point coordinates, where N>1 and M is the
%             number of point samples.
%   - R     : radius of the bounding sphere.
%   - C     : centroid of the bounding sphere.
%
% REFERENCES:
% [1] Ritter, J. (1990), 'An efficient bounding sphere', in Graphics Gems,
%     A. Glassner, Ed. Academic Press, pp.301�03
%
% AUTHOR: Anton Semechko (a.semechko@gmail.com)
% DATE: Jan.2014
%

% Basic error checking
if size(X,2)<2
    error('Invalid entry for 1st input argument')
end
if size(X,1)==1
    R=0; C=X;
    return
end

% Get the convex hull of the point set
F=convhulln(X);
F=unique(F(:));
X=X(F,:);

% Find points with the most extreme coordinates
[~,idx_min]=min(X,[],1);
[~,idx_max]=max(X,[],1);

Xmin=X(idx_min(:),:);
Xmax=X(idx_max(:),:);

% Compute distance between the bounding box points
X1=permute(Xmin,[1 3 2]);
X2=permute(Xmax,[3 1 2]);
D2=bsxfun(@minus,X2,X1);
D2=sum(D2.^2,3);

% Select the pair with the largest distance
[D2_max,idx]=max(D2(:));
[i,j]=ind2sub(size(D2),idx);
Xmin=Xmin(i,:);
Xmax=Xmax(j,:);

% Initial radius (squared) and centroid
R2=D2_max/4;
R=sqrt(R2);
C=(Xmin+Xmax)/2;

% Loop through the M points adjusting the position and radius of the sphere
M=size(X,1);
for i=1:M
    
    di2=sum((X(i,:)-C).^2);
    if di2>R2
        di=sqrt(di2);
        R=(R+di)/2;
        R2=R^2;
        dR=di-R;
        C=(R*C+dR*X(i,:))/di;
    end
    
end

输入点云，输出的是球体的半径和球心

我的测试例子

测试代码为

function test
%TEST Summary of this function goes here
%   Detailed explanation goes here
[vertices,faces]=obj__read('mannequin.obj');
[R,C]=ApproxMinBoundSphereND(vertices');
[x,y,z]=sphere;
mesh(R*x+C(1),R*y+C(2),R*z+C(3),'EdgeColor','k');alpha(0.1);axis equal;hold on;trimesh(faces',vertices(1,:)',vertices(2,:)',vertices(3,:)');
end

效果如下