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

算法过程

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;

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代码

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
%
% 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

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)';

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)';

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));

08-24

08-15
06-30
09-17
05-23
05-28 4478
10-04 4933
02-27
09-17
04-18 3726
07-24 2771
11-06 5547
08-03 6807
07-10 3376
06-27 1124
08-21 1197
07-26 6881
12-19 5415
11-18 28