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

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

 

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103

%该函数是在其他网站看到的，以此共享。有两种方法（函数）实现。 %第一种比较费时： 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 = [];   %最主要的就是后两个函数了，第一个只是测试。