FITCIRCLE3D Fit a 3D circle to a set of points.
[FITTEDCIRCLE, CIRCLENORMAL, RESIDUALS] = fitCircle3d(PTS)
Example
% points on a 2d circle with noise
nop = randi([5 50],1,1);
radius = randi([5 25],1,1);
points2d = circleToPolygon([0 0 radius], nop);
points2d(1,:) = [];
points2d = points2d + rand(size(points2d));
points2d(:,3)=rand(length(nop),1);
% apply random rotation and translation
[theta, phi] = randomAngle3d;
theta = rad2deg(theta);
phi = rad2deg(phi);
tfm = eulerAnglesToRotation3d(phi, theta, 0);
trans = randi([-250 250],3,1);
tfm(1:3,4)=trans;
points3d = awgn(transformPoint3d(points2d, tfm),1);
% fit 3d circle
[fittedCircle, circleNormal] = fitCircle3d(points3d);
% plot 3d points and 3d circle
figure('Color','w'); hold on; axis equal tight; view(3);
xlabel('X');ylabel('Y');zlabel('Z');
drawPoint3d(points3d)
drawCircle3d(fittedCircle, 'k')
drawVector3d(fittedCircle(1:3), circleNormal*fittedCircle(4))
See also
circle3dOrigin, circle3dPosition, circle3dPoint, intersectPlaneSphere
drawCircle3d, drawCircleArc3d, drawEllipse3d0001 function [fittedCircle, circleNormal, residuals] = fitCircle3d(pts, varargin) 0002 %FITCIRCLE3D Fit a 3D circle to a set of points. 0003 % 0004 % [FITTEDCIRCLE, CIRCLENORMAL, RESIDUALS] = fitCircle3d(PTS) 0005 % 0006 % Example 0007 % % points on a 2d circle with noise 0008 % nop = randi([5 50],1,1); 0009 % radius = randi([5 25],1,1); 0010 % points2d = circleToPolygon([0 0 radius], nop); 0011 % points2d(1,:) = []; 0012 % points2d = points2d + rand(size(points2d)); 0013 % points2d(:,3)=rand(length(nop),1); 0014 % % apply random rotation and translation 0015 % [theta, phi] = randomAngle3d; 0016 % theta = rad2deg(theta); 0017 % phi = rad2deg(phi); 0018 % tfm = eulerAnglesToRotation3d(phi, theta, 0); 0019 % trans = randi([-250 250],3,1); 0020 % tfm(1:3,4)=trans; 0021 % points3d = awgn(transformPoint3d(points2d, tfm),1); 0022 % % fit 3d circle 0023 % [fittedCircle, circleNormal] = fitCircle3d(points3d); 0024 % % plot 3d points and 3d circle 0025 % figure('Color','w'); hold on; axis equal tight; view(3); 0026 % xlabel('X');ylabel('Y');zlabel('Z'); 0027 % drawPoint3d(points3d) 0028 % drawCircle3d(fittedCircle, 'k') 0029 % drawVector3d(fittedCircle(1:3), circleNormal*fittedCircle(4)) 0030 % 0031 % See also 0032 % circle3dOrigin, circle3dPosition, circle3dPoint, intersectPlaneSphere 0033 % drawCircle3d, drawCircleArc3d, drawEllipse3d 0034 0035 % ------ 0036 % Authors: oqilipo 0037 % E-mail: N/A 0038 % Created: 2017-05-09 0039 % Copyright 2017-2024 0040 0041 parser = inputParser; 0042 addRequired(parser, 'pts', @(x) validateattributes(x, {'numeric'},... 0043 {'ncols',3,'real','finite','nonnan'})); 0044 addOptional(parser,'verbose',true,@islogical); 0045 parse(parser,pts,varargin{:}); 0046 pts = parser.Results.pts; 0047 verbose = parser.Results.verbose; 0048 0049 % Mean of all points 0050 meanPoint = mean(pts,1); 0051 0052 % Center points by subtracting the meanPoint 0053 centeredPoints = pts - repmat(meanPoint,size(pts,1),1); 0054 0055 % Project 3D data to a plane 0056 [~,~,V]=svd(centeredPoints); 0057 tfmPoints = transformPoint3d(centeredPoints, V'); 0058 0059 % Fit a circle to the points in the xy-plane 0060 circleParamter = CircleFitByTaubin(tfmPoints(:,1:2), verbose); 0061 center2d = circleParamter(1:2); 0062 radius=circleParamter(3); 0063 center3d = transformPoint3d([center2d, 0], [inv(V'), meanPoint']); 0064 circleNormal = V(:,3)'; 0065 [theta, phi, ~] = cart2sph2(circleNormal); 0066 fittedCircle = [center3d radius rad2deg(theta) rad2deg(phi) 0]; 0067 0068 % Residuals 0069 residuals = radius - distancePoints(tfmPoints(:,1:2),center2d); 0070 0071 end 0072 0073 % Circle Fit (Taubin method) 0074 % version 1.0 (2.24 KB) by Nikolai Chernov 0075 % http://www.mathworks.com/matlabcentral/fileexchange/22678 0076 function Par = CircleFitByTaubin(XY, varargin) 0077 %__________________________________________________________________________ 0078 % 0079 % Circle fit by Taubin 0080 % G. Taubin, "Estimation Of Planar Curves, Surfaces And Nonplanar 0081 % Space Curves Defined By Implicit Equations, With 0082 % Applications To Edge And Range Image Segmentation", 0083 % IEEE Trans. PAMI, Vol. 13, pages 1115-1138, (1991) 0084 % 0085 % Input: XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2) 0086 % 0087 % Output: Par = [a b R] is the fitting circle: 0088 % center (a,b) and radius R 0089 % 0090 % Note: this fit does not use built-in matrix functions (except "mean"), 0091 % so it can be easily programmed in any programming language 0092 % 0093 %__________________________________________________________________________ 0094 0095 parser = inputParser; 0096 addRequired(parser, 'XY', @(x) validateattributes(x, {'numeric'},... 0097 {'ncols',2,'real','finite','nonnan'})); 0098 addOptional(parser,'verbose',true,@islogical); 0099 parse(parser,XY,varargin{:}); 0100 XY=parser.Results.XY; 0101 verbose = parser.Results.verbose; 0102 0103 n = size(XY,1); % number of data points 0104 0105 centroid = mean(XY); % the centroid of the data set 0106 0107 % computing moments (note: all moments will be normed, i.e. divided by n) 0108 Mxx = 0; Myy = 0; Mxy = 0; Mxz = 0; Myz = 0; Mzz = 0; 0109 0110 for i=1:n 0111 Xi = XY(i,1) - centroid(1); % centering data 0112 Yi = XY(i,2) - centroid(2); % centering data 0113 Zi = Xi*Xi + Yi*Yi; 0114 Mxy = Mxy + Xi*Yi; 0115 Mxx = Mxx + Xi*Xi; 0116 Myy = Myy + Yi*Yi; 0117 Mxz = Mxz + Xi*Zi; 0118 Myz = Myz + Yi*Zi; 0119 Mzz = Mzz + Zi*Zi; 0120 end 0121 0122 Mxx = Mxx/n; 0123 Myy = Myy/n; 0124 Mxy = Mxy/n; 0125 Mxz = Mxz/n; 0126 Myz = Myz/n; 0127 Mzz = Mzz/n; 0128 0129 % computing the coefficients of the characteristic polynomial 0130 Mz = Mxx + Myy; 0131 Cov_xy = Mxx*Myy - Mxy*Mxy; 0132 A3 = 4*Mz; 0133 A2 = -3*Mz*Mz - Mzz; 0134 A1 = Mzz*Mz + 4*Cov_xy*Mz - Mxz*Mxz - Myz*Myz - Mz*Mz*Mz; 0135 A0 = Mxz*Mxz*Myy + Myz*Myz*Mxx - Mzz*Cov_xy - 2*Mxz*Myz*Mxy + Mz*Mz*Cov_xy; 0136 A22 = A2 + A2; 0137 A33 = A3 + A3 + A3; 0138 0139 xnew = 0; 0140 ynew = 1e+20; 0141 epsilon = 1e-12; 0142 IterMax = 20; 0143 0144 % Newton's method starting at x=0 0145 for iter=1:IterMax 0146 yold = ynew; 0147 ynew = A0 + xnew*(A1 + xnew*(A2 + xnew*A3)); 0148 if abs(ynew) > abs(yold) 0149 if verbose 0150 disp('Newton-Taubin goes wrong direction: |ynew| > |yold|'); 0151 end 0152 xnew = 0; 0153 break; 0154 end 0155 Dy = A1 + xnew*(A22 + xnew*A33); 0156 xold = xnew; 0157 xnew = xold - ynew/Dy; 0158 if (abs((xnew-xold)/xnew) < epsilon), break, end 0159 if (iter >= IterMax) 0160 if verbose 0161 disp('Newton-Taubin will not converge'); 0162 end 0163 xnew = 0; 0164 end 0165 if (xnew<0.) 0166 if verbose 0167 fprintf(1,'Newton-Taubin negative root: x=%f\n',xnew); 0168 end 0169 xnew = 0; 0170 end 0171 end 0172 0173 % computing the circle parameters 0174 DET = xnew*xnew - xnew*Mz + Cov_xy; 0175 Center = [Mxz*(Myy-xnew)-Myz*Mxy , Myz*(Mxx-xnew)-Mxz*Mxy]/DET/2; 0176 0177 Par = [Center+centroid , sqrt(Center*Center'+Mz)]; 0178 0179 end