算法来自1991年EMO WELZL发表的一篇文章:Smallest enclosing disks (balls and ellipsoids)
可以精确计算出三维点云的用一个闭球将其包含进去,同时也可以将其拓展为n维空间中的应用
算法过程
给定平面中的n个点的集合P,令md(P)表示包含P中的所有点的最小半径的封闭盘。我们允许当
P=∅
时,
md(P)=∅
和
P=p
时,
md(P)=p
。
下面的例子全部针对于二维情况(三维情况可以类比得出)。可以证明得到的最小圆盘是唯一的。还可以得到一个结论:确定最小圆最多只需要点集中的至多三个点。
我们以增量的方式来计算最小圆盘。将点一个一个加入,如果新增加的点在这个圆盘里,直接进入下一个点的计算。如果不在的话,则这个点一定在新确定的圆的边界上。我们需要调用一个函数b_minidisk(A, p)来计算新的最小圆,A包含虽有旧点的闭圆,p在新圆的边界上。
这里我们先假设b_minidisk已经存在,固定点p,算法可以以下面的伪代码实现
function procedure minidisk(P); comment: returns md(P)
if P = ∅ then
D:= ∅
else
choose p ∈ P;
D := minidisk(P - {p});
if p ∉ 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 = ∅ [or |R|=3]then
D:=b_md( ∅ ,R)
else
choose random p ∈ P;
D := B_MINIDISK(P - {p},R);
if [D defined and] p ∉ D then
D := B_MINIDISK(P-{p}, R ∪ {p});
return D;
因此初始问题可以由下面的伪代码解决
function procedure MINIDISK(P); comment: returns md(P)
return B_MINIDTSK(P, ∅ ) ;
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));