필터 지우기
필터 지우기

Non linear least square fitting

조회 수: 42 (최근 30일)
Jenny Andersen
Jenny Andersen 2019년 12월 9일
편집: cui,xingxing 2024년 4월 27일
Hi,
I am really stuck trying to figure out how to fit a circle to some data points. I have watched plenty of videos on the topic but do not really understand. Could you please explain to me in the most simple way how to plot some random data points, for example x = [2; 3; -2; 1]; and y = [3; 2; 4; -1];, in matlab?
I have only come up with the following but do not know how to proceed:
x = [2; 3; -2; 1];
y = [3; 2; 4; -1];
N = length(x);
% circle's equation x^2+y^2 = 2xc1+2yc2+c3
X = [ones(N,3),x];
Y=y;
phi = inv(X'*X)*X'*Y
plot(x,y,'*')
Thanks alot in advance!
  댓글 수: 4
Matt J
Matt J 2019년 12월 9일
Jenny's answer converted to a comment:
I tried the following this time, but it does not seem to be working. Anyone knows why?
x = [2; 3; -2; 1];
y = [3; 2; 4; -1];
N = length(x);
plot (x,y,'*')
% circle's equation x^2+y^2 = 2xc1+2yc2+c3
x = c1+sqrt(-y^2+2ay+a^2+c
y = c1+sqrt(-y^2+2ax+a^2+c
X = [ones(N,3),x];
Y = y;
plot(X,Y)
Matt J
Matt J 2019년 12월 9일
편집: Matt J 2019년 12월 9일
I have only come up with the following but do not know how to proceed ... phi = inv(X'*X)*X'*Y
You have tried to set up linear equations X*phi=Y and find its Ordinary Least Squares (OLS) solution. This is not a good approach here, because X and Y both have noise in them. But even if you were to do it this way, the better Matlab algorithm for OLS would be,
phi=X\Y

댓글을 달려면 로그인하십시오.

채택된 답변

Adam Danz
Adam Danz 2019년 12월 9일
편집: Adam Danz 2019년 12월 9일
Create noisy circular data
This section creates the noisy data we're fitting.
% Create randome points along a circle
r = 22.52; % radius
cnt = [9.923, -2.42]; % (x,y) center
theta = rand(1,1000)*2*pi; % Theta
xData = r*cos(theta)+cnt(1); % x coordinates
yData = r*sin(theta)+cnt(2); % y coordinates
% Add noise to coordinates
xData = xData + (rand(size(xData))*10-5);
yData = yData + (rand(size(yData))*10-5);
Plot the input data
The only inputs are xData and yData. The rest of the variables above are unknowns.
% Plot the data with noise
clf()
plot(xData, yData, 'ko')
axis equal
grid on
Fit the xData,yData coordinates to a circle
This section uses nonlinear least squares fitting x = lsqnonlin(fun,x0). The first line defines the function to fit and is the equation for a circle. The second line are estimated starting points. See the link for more info on this function. The output circFit is a 1x3 vector defining the [x_center, y_center, radius] of the fitted circle.
f = @(a) (xData-a(1)).^2 + (yData-a(2)).^2 - a(3).^2;
a0 = [mean(xData),mean(yData),max(xData)-mean(xData)];
circFit = lsqnonlin(f,a0);
Add the fitted circle to the noisy data
% use the circFit parameters to create the fitted circle
theta = linspace(0,2*pi,100); % arbitrary spacing
xFit = circFit(3)*cos(theta) + circFit(1);
yFit = circFit(3)*sin(theta) + circFit(2);
hold on
plot(xFit,yFit,'r-','LineWidth',3)
Compare the fitted values to the actual (noisy) values
Of course these values may differ when this code is repeated due to random noise.
This creates a table comparing the known and fitted results.
T = array2table([[cnt,r];circFit],'VariableNames',{'xCnt','yCnt','radius'},...
'RowNames',{'Real values','FitValues'})
T =
2×3 table
xCnt yCnt radius
______ ______ ______
Real values 9.923 -2.42 22.52
FitValues 9.8052 -2.548 22.926
  댓글 수: 4
Jenny Andersen
Jenny Andersen 2019년 12월 9일
Thanks alot Adam btw.
Adam Danz
Adam Danz 2019년 12월 9일
편집: Adam Danz 2019년 12월 9일
Thanks, @Matt J. I meant to mention alternative approaches to the special case when fitting a circle but by the time I finished typing out the answer, I had forgotten.
@Jenny Andersen, there are several option on the file exchange.
But Matt J's answer is even better since you won't need to use a function external to what Matlab already provides.

댓글을 달려면 로그인하십시오.

추가 답변 (2개)

Matt J
Matt J 2019년 12월 9일
편집: Matt J 2019년 12월 9일
Ok. Could you possibly show me the analytical alternatives to Adam's solution?
See attached. Note that my method does not minimize the same fitting error metric as Adam's. If you prefer his, though (and maybe you should), it might be worth mentioning that I do observe somewhat better convergence when initialized with mine. Notice the smaller final "Norm of step" and "First-order optimality" when the iterations terminate.
f = @(a) (xData-a(1)).^2 + (yData-a(2)).^2 - a(3).^2;
opts=optimoptions('lsqnonlin','Display','iter');
a0 = [mean(xData),mean(yData),max(xData)-mean(xData)];
lsqnonlin(f,a0,[],[],opts);
Norm of First-order
Iteration Func-count f(x) step optimality
0 4 6.04793e+07 1.13e+07
1 8 1.69499e+07 3.881 6.82e+05
2 12 1.67336e+07 0.317084 4.6e+03
3 16 1.67336e+07 0.00219778 0.221
[c,R]=circfitAlgebraic([xData;yData]);
lsqnonlin(f,[c,R],[],[],opts);
Norm of First-order
Iteration Func-count f(x) step optimality
0 4 1.74094e+07 1.22e+06
1 8 1.67337e+07 0.555282 1.41e+04
2 12 1.67336e+07 0.00672023 2.07
3 16 1.67336e+07 9.87177e-07 0.00653

cui,xingxing
cui,xingxing 2023년 4월 29일
편집: cui,xingxing 2024년 4월 27일
I offer here a robust least squares fit for reference purposes only.
%% Create randome points along a circle
r = 1;
numPts = 100;
theta = linspace(0,2*pi,numPts);
xdata = r*cos(theta)+0.1*randn(1,numPts);
ydata = r*sin(theta)+0.1*randn(1,numPts);
% add noise
xdata(1) = 5;xdata(2) = 4;xdata(3) = 5;
ydata(1) = 4;ydata(2) = 5;ydata(3) = 5;
h1 = scatter(xdata,ydata); % sample data
% noise least square fit circle
fun = @(x)(xdata-x(1)).^2+(ydata-x(2)).^2-x(3).^2;
x = lsqnonlin(fun,[0,0,1]);
Local minimum possible. lsqnonlin stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.
h2 = viscircles(x(1:2),x(3),color="red");hold on; % noise least square fit
h3 = viscircles([0,0],1,color="blue");% ground truth
%% robust fit
fitFcn = @(points)robustLeastSquareCircle(points);
distFcn = @(model, points)((xdata-model(1)).^2+(ydata-model(2)).^2-model(3).^2).^2;
sampleSize = 3;
maxDistance = 0.5;
data = [xdata(:),ydata(:)];
[modelRansac,inlierIdx] = ransac(data,fitFcn,distFcn,sampleSize,maxDistance);
Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance. Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.
h4 = viscircles(modelRansac(1:2),modelRansac(3),color="black");
modelRobust = fitFcn(data(inlierIdx,:));
Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.
h5 = viscircles(modelRobust(1:2),modelRobust(3),color="green");% robust fit'
legend([h1,h2,h3,h4,h5],{'point samples','least square fit','ground truth','Ransac samples fit','robust fit'},Location="northwest")
grid on;axis equal
%% support function
function model = robustLeastSquareCircle(points)
% least square fit circle
xdata = points(:,1);
ydata = points(:,2);
fun = @(x)(xdata-x(1)).^2+(ydata-x(2)).^2-x(3).^2;
model = lsqnonlin(fun,[0,0,1]);
end
-------------------------Off-topic interlude, 2024-------------------------------
I am currently looking for a job in the field of CV algorithm development, based in Shenzhen, Guangdong, China,or a remote support position. I would be very grateful if anyone is willing to offer me a job or make a recommendation. My preliminary resume can be found at: https://cuixing158.github.io/about/ . Thank you!
Email: cuixingxing150@gmail.com
  댓글 수: 2
Matt J
Matt J 2023년 4월 30일
This looks pretty much the same as @Adam Danz's solution.
cui,xingxing
cui,xingxing 2023년 5월 8일
I have improved it slightly in the back

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Shifting and Sorting Matrices에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by