基于已知点云数据的最小外接圆matlab函数

基于已知点云数据的最小外接圆matlab函数 – MATLAB中文论坛 (ilovematlab.cn)

 

%该函数是在其他网站看到的，以此共享。有两种方法（函数）实现。
%第一种比较费时：
function [xc,yc,r] = smallestcircle(x,y)

% This finds the circle of smallest area containing all
% the points of vectors x and y.  The center is (xc,yc)
% and the radius is r.   RAS - 2/28/05

if any(size(x)~=size(y)) | min(size(x))>1
error('Args must be vectors of equal length.')
end
n = length(x);
if n == 1, xc = x(1); yc = y(1); r = 0;
else
  bestr2 = inf;
  for i = 1:n-1
    x1 = x(i); y1 = y(i);
    for j = i+1:n
      x2 = x(j); y2 = y(j); x0 = (x1+x2)/2; y0 = (y1+y2)/2;
      w = sqrt((x2-x1)^2+(y2-y1)^2);
      if w==0, break, end
      A = (x2-x1)/w; B = (y2-y1)/w; C = (x1*y2-y1*x2)/w;
      h = w/2; h2 = h^2;
      pmax = -inf; pmin = inf;
      for k = [1:i-1,i+1:j-1,j+1:n]
        x3 = x(k); y3 = y(k);
        o = A*y3-B*x3+C;
        if o==0, o = 2^(-1074); end
        p = ((x3-x0)^2+(y3-y0)^2-h2)/(2*o);
        if o > 0, if p > pmax, pmax = p; uk = k; end
        else if p < pmin, pmin = p; lk = k; end
        end
        if pmax > pmin, break, end
      end
      if pmax <= pmin
        if pmax > 0, p = pmax; k = uk;
        elseif pmin < 0, p = pmin; k = lk;
        else xc = x0; yc = y0; r = h; return
        end
        r2 = p^2+h2;
        if r2 < bestr2, bestr2 = r2; si = i; sj = j; sk = k; end
      end
    end
  end
  a = x(si); b = x(sj); c = x(sk);
  d = y(si); e = y(sj); f = y(sk);
  t = 2*(a*e+b*f+c*d-b*d-c*e-a*f);
  xc = ((a^2+d^2)*(e-f)+(b^2+e^2)*(f-d)+(c^2+f^2)*(d-e))/t;
  yc = ((a^2+d^2)*(c-b)+(b^2+e^2)*(a-c)+(c^2+f^2)*(b-a))/t;
  r = sqrt(((b-a)^2+(e-d)^2)*((a-c)^2+(d-f)^2)*((c-b)^2+(f-e)^2))/abs(t);
end

第二种比较优化：
function testcircumcircle
%% Test the CIRCUMCIRCLE function
% Generate random data
x = randn(20, 1);
y = randn(20, 1);
% Call CIRCUMCIRCLE
[c, r] = circumcircle(x, y);
% Plot the data as red + symbols
plot(x, y, 'r+');
% Wait for the next plot ...
hold on;
% Create a vector of theta values
dtheta = 0.1;
theta = 0:dtheta:(2*pi+dtheta);
% Plot the circle whose center is [c(1), c(2)] and whose radius is r
plot(r*cos(theta)+c(1), r*sin(theta)+c(2), 'g-')


function [c, r] = circumcircle(x, y)
%% The CIRCUMCIRCLE function
% Generate a feasible random guess
% The last elemeent of p0 is the maximum distance from the origin to a
% point in our data set. All points must be within that distance from the
% origin, so we're guaranteed that the nonlinear constraint in TestIfInside
% will be satisfied
p0 = [0; 0; sqrt(max(x.^2+y.^2))];
% Creating an options structure
options = optimset('Display', 'iter');
% Call FMINCON
P = fmincon(@(p) p(3), p0, [], [], [], [], [-Inf;-Inf;0], [], ...
   @(p) TestIfInside(p, x, y), options);
% Note that for versions of MATLAB prior to MATLAB 7.0 (R14), you would
% need to replace the anonymous functions @(p) p(3) and
% @(p) TestIfInside(p, x, y) with inline functions or regular function
% handles
% Extracting the elements of the P vector; the first two are the center of
% the circle, the third is the radius
c = P(1:2);
r = P(3);


function [c, ceq] = TestIfInside(p, x, y)
%% Nonlinear constraint function for CIRCUMCIRCLE
% The squared distance between our points and the center of the circle is
% less than the radius squared
c = (p(1)-x).^2 + (p(2)-y).^2 - p(3).^2;
% No nonlinear equality constraint
ceq = [];

%最主要的就是后两个函数了，第一个只是测试。