一种得到点云精确边界球的方法

算法来自1991年EMO WELZL发表的一篇文章：Smallest enclosing disks (balls and ellipsoids)

可以精确计算出三维点云的用一个闭球将其包含进去，同时也可以将其拓展为n维空间中的应用

算法过程

给定平面中的n个点的集合P，令md（P）表示包含P中的所有点的最小半径的封闭盘。我们允许当$P=\emptyset$时，$md（P）=\emptyset$和$P = {p}$时，$md（P）= p$。
下面的例子全部针对于二维情况（三维情况可以类比得出）。可以证明得到的最小圆盘是唯一的。还可以得到一个结论：确定最小圆最多只需要点集中的至多三个点。
我们以增量的方式来计算最小圆盘。将点一个一个加入，如果新增加的点在这个圆盘里，直接进入下一个点的计算。如果不在的话，则这个点一定在新确定的圆的边界上。我们需要调用一个函数b_minidisk(A, p)来计算新的最小圆，A包含虽有旧点的闭圆，p在新圆的边界上。
这里我们先假设b_minidisk已经存在，固定点p，算法可以以下面的伪代码实现

function procedure minidisk(P); comment: returns md(P)
$\ \ \ $if P = $\emptyset$ then
$\ \ \ \ \ \ $D:=$\emptyset$
$\ \ \ $else
$\ \ \ \ \ \ $choose p$\in$P;
$\ \ \ \ \ \ $D := minidisk(P - {p});
$\ \ \ \ \ \ $if p $\notin$ D then
$\ \ \ \ \ \ \ \ \ $D := b_minidisk(P - {p}, p);
$\ \ \ $return D;

至于b_minidisk函数，由下面的伪代码给出

function procedure B_MINIDISK(P,R); comment: returns b_md(P,R)
$\ \ \ $if P = $\emptyset$ [or |R|=3]then
$\ \ \ \ \ \ $D:=b_md($\emptyset$,R)
$\ \ \ $else
$\ \ \ \ \ \ $choose random p$\in$P;
$\ \ \ \ \ \ $D := B_MINIDISK(P - {p},R);
$\ \ \ \ \ \ $if [D defined and] p $\notin$ D then
$\ \ \ \ \ \ \ \ \ $D := B_MINIDISK(P-{p}, R$\cup${p});
$\ \ \ $return D;

因此初始问题可以由下面的伪代码解决

function procedure MINIDISK(P); comment: returns md(P)
$\ \ \ $return B_MINIDTSK(P, $\emptyset$) ;

Matlab代码

对于三维的而言ExactMinBoundSphere3D.m

function [R,C,Xb]=ExactMinBoundSphere3D(X)
% Compute exact minimum bounding sphere of a 3D point cloud (or a 
% triangular surface mesh) using Welzl's algorithm. 
%
%   - X     : M-by-3 list of point co-ordinates or a triangular surface 
%             mesh specified as a TriRep object.
%   - R     : radius of the sphere.
%   - C     : 1-by-3 vector specifying the centroid of the sphere.
%   - Xb    : subset of X, listing K-by-3 list of point coordinates from 
%             which R and C were computed. See function titled 
%             'FitSphere2Points' for more info.
%
% REREFERENCES:
% [1] Welzl, E. (1991), 'Smallest enclosing disks (balls and ellipsoids)',
%     Lecture Notes in Computer Science, Vol. 555, pp. 359-370
%
% AUTHOR: Anton Semechko (a.semechko@gmail.com)
% DATE: Dec.2014
%


if isobject(X), X=X.X; end
if size(X,2)~=3
    error('This function only works for 3D data')
end
if sum(isnan(X(:)) | isinf(X(:)))>0
    error('Point data contains NaN or Inf entries. Remove them and try again.')
end

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

% Randomly permute the point set
idx=randperm(size(X,1));
X=X(idx(:),:);

% Get the minimum bounding sphere
if size(X,1)<=4
    [R,C]=FitSphere2Points(X); 
    Xb=X;
    return
end

if size(X,1)<1E3
    try

        % Center and radius of the sphere
        [R,C]=B_MinSphere(X,[]);

        % Coordiantes of the points used to compute parameters of the 
        % minimum bounding sphere
        D=sum(bsxfun(@minus,X,C).^2,2);
        [D,idx]=sort(abs(D-R^2));
        Xb=X(idx(1:4),:);
        D=D(1:4);
        Xb=Xb(D<1E-6,:);
        [~,idx]=sort(Xb(:,1));
        Xb=Xb(idx,:);
        return
    catch
    end
end

% If we got to this point, then the recursion depth limit was reached. So 
% need to break-up the the data into smaller sets and then recombine the 
% results.
M=size(X,1);
dM=min(floor(M/4),300);
res=mod(M,dM);
n=ceil(M/dM);  
idx=dM*ones(1,n);
if res>0
    idx(end)=res;
end

if res<=0.25*dM 
    idx(n-1)=idx(n-1)+idx(n);
    idx(n)=[];
    n=n-1;
end

X=mat2cell(X,idx,3);
Xb=[];
for i=1:n

    % Center and radius of the sphere
    [R,C,Xi]=B_MinSphere([Xb;X{i}],[]);    

    % 40 points closest to the sphere
    if i<1
        D=abs(sum(bsxfun(@minus,Xi,C).^2,2)-R^2);
    else
        D=abs(sqrt(sum(bsxfun(@minus,Xi,C).^2,2))-R);
    end
    [D,idx]=sort(D);
    Xb=Xi(idx(1:40),:);

end
D=D(1:4);
Xb=Xb(D/R*100<1E-3,:);
[~,idx]=sort(Xb(:,1));
Xb=Xb(idx,:);


    function [R,C,P]=B_MinSphere(P,B)

    if size(B,1)==4 || isempty(P)
        [R,C]=FitSphere2Points(B); % fit sphere to boundary points
        return
    end

    % Remove the last (i.e., end) point, p, from the list
    P_new=P;
    P_new(end,:)=[];
    p=P(end,:);

    % Check if p is on or inside the bounding sphere. If not, it must be
    % part of the new boundary.
    [R,C,P_new]=B_MinSphere(P_new,B); 
    if isnan(R) || isinf(R) || R<=eps
        chk=true;
    else
        chk=norm(p-C)>(R+eps);
    end

    if chk
        B=[p;B];
        [R,C]=B_MinSphere(P_new,B);
        P=[p;P_new];
    end

    end


end

以及FitSphere2Points.m

function [R,C]=FitSphere2Points(X)
% Fit a sphere to a set of 2, 3, or at most 4 points in 3D space. Note that
% point configurations with 3 collinear or 4 coplanar points do not have 
% well-defined solutions (i.e., they lie on spheres with infinite radius).
%
%   - X     : M-by-3 array of point coordinates, where M<=4.
%   - R     : radius of the sphere. R=Inf when the sphere is undefined, as 
%             specified above.
%   - C     : 1-by-3 vector specifying the centroid of the sphere. 
%             C=nan(1,3) when the sphere is undefined, as specified above.
%
% AUTHOR: Anton Semechko (a.semechko@gmail.com)
% DATE: Dec.2014
%


N=size(X,1);
if N>4 
    error('Input must a N-by-3 array of point coordinates, with N<=4')
end

% Empty set
if isempty(X)
    C=nan(1,3);
    R=nan; 
    return
end

% A single point
if N==1
    C=X;
    R=0;
    return
end

% Line segment
if N==2
    C=mean(X,1);
    R=norm(X(2,:)-X(1,:))/2;
    return
end

% Remove duplicate vertices, if there are any
D=bsxfun(@minus,permute(X,[1 3 2]),permute(X,[3 1 2]));
D=sqrt(sum(D.^2,3));
D(1:(N+1):end)=Inf;
chk=D<=1E-12;
if sum(chk(:))>0
    for i=1:N
        if size(X,1)<=i, break; end
        idx=chk(i,:);
        idx(1:i)=false;
        idx=find(idx);
        chk(idx,:)=[];
        chk(:,idx)=[];
        X(idx,:)=[];
    end
    [R,C]=FitSphere2Points(X);
    return
end


% Three unique, though possibly collinear points
tol=1E-2; % collinearity/co-planarity threshold (in degrees)
if N==3

    % Check for collinearity
    D12=X(2,:)-X(1,:); D12=D12/norm(D12);
    D13=X(3,:)-X(1,:); D13=D13/norm(D13);

    chk=abs(D12*D13(:));
    chk(chk>1)=1;
    if acos(chk)/pi*180<tol
        R=inf;
        C=nan(1,3);
        return
    end

    % Make plane formed by the points parallel with the xy-plane  
    n=cross(D13,D12);
    n=n/norm(n);
    r=cross(n,[0 0 1]); r=acos(n(3))*r/norm(r); % Euler rotation vector
    Rmat=expm([0 -r(3) r(2); r(3) 0 -r(1); -r(2) r(1) 0]);
    Xr=(Rmat*X')'; 

    % Circle centroid
    x=Xr(:,1:2);
    A=2*bsxfun(@minus,x(2:3,:),x(1,:));
    b=sum(bsxfun(@minus,x(2:3,:).^2,x(1,:).^2),2);
    C=(A\b)';

    % Circle radius
    R=sqrt(sum((x(1,:)-C).^2,2));     

    % Rotate centroid back into the original frame of reference
    C(3)=mean(Xr(:,3));
    C=(Rmat'*C(:))';
    return
end


% If we got to this point then we have 4 unique, though possibly co-linear
% or co-planar points. 

% Check if the the points are co-linear
D12=X(2,:)-X(1,:); D12=D12/norm(D12);
D13=X(3,:)-X(1,:); D13=D13/norm(D13);
D14=X(4,:)-X(1,:); D14=D14/norm(D14);

chk1=abs(D12*D13(:));
chk1(chk1>1)=1;
chk2=abs(D12*D14(:));
chk2(chk2>1)=1;

if acos(chk1)/pi*180<tol || acos(chk2)/pi*180<tol
    R=inf;
    C=nan(1,3);
    return
end

% Check if the the points are co-planar
n1=cross(D12,D13); n1=norm(n1);
n2=cross(D12,D14); n2=norm(n2);

chk=abs(n1*n2(:));
chk(chk>1)=1;
if acos(chk)/pi*180<tol
    R=inf;
    C=nan(1,3);
    return
end

% Centroid of the sphere
A=2*bsxfun(@minus,X(2:end,:),X(1,:));
b=sum(bsxfun(@minus,X(2:end,:).^2,X(1,:).^2),2);
C=(A\b)';

% Radius of the sphere
R=sqrt(sum((X(1,:)-C).^2,2));

距离测试

再次拿上次的人头像举例

[vertex,face]=obj__read('mannequin.obj');
X=vertex';
[R,C,Xb]=ExactMinBoundSphere3D(X);
[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(face',X(:,1),X(:,2),X(:,3));