[Math] Fitting a set of points in a plane to a smooth curve obtained by joining a half-line and an arc of a circle

approximationplane-curves

I have a set of points in the plane and I want to find a curve that best fits these points (e.g., in a least squares manner, or using some other convenient "measure"). I want that the curve be either:

  • a line

  • a circle

  • a circle whose center lies on a given line

  • a curve that looks like a hook, that is obtained from:

    (a) an arc of a circle whose center lies on a given line

    (b) a half-line tangent to one of the endpoints of the arc circle (in this way I can get a curve that is smooth at the point where the arc of the circle and the half-line meet).

I managed to get an answer for all but the last problem. Any suggestions are welcomed.

Best Answer

Here is a numerical solution. The parameters are taken as $A$,$B$,$C$,$M$,$K$,$c$. The equation of the arc is given by $$(y-C)^2+(x-B)^2=A^2$$ and for the line $$y=M\,x+K$$ The break-point is taken as $c$. Therefore the objective function to be minimized becomes $$F=\sum_1^n\bigg((y_i-F_{1,i})^2+(y_i-F_{2,i})^2\bigg)$$ $$F_{1,i}=-\sqrt{A^2-(x_i-B)^2}+C\qquad \text{if }x_i\le c\text{ and }0\text{ elsewhere}$$ $$F_{2,i}=M\,x_i+K\qquad \text{if }x_i>c\text{ and }0\text{ elsewhere}$$ The first constraint function assures that the arc and line meet at the break-point $c$ $$-\sqrt{A^2-(c-B)^2}+C-M\,c-K=0$$ The second constraint that the line is tangent to the arc at $c$ $$\frac{c-B}{\sqrt{A^2-(x-B)^2}}-M=0$$ And the last constraint is that the center of arc lies on a line given by $y=m_1x+k$ $$C-m_1\,B+k_1=0$$ To test the algorithm I generated 120 points and added error terms of $N(0,0.1)$ as given in below graph. The black dots are for circle with $A=3$, $B=1$, $C=4$; and the blue dots for line with $M=0.8944$, $K=-0.9193$. The break-point is given as $c=3$ Data

In Matlab I used "fmincon" with "interior-point" algorithm. The initial guess is randomly generated as $[A=0.9045;B=0.8663;C=0.9225;M=0.5637;K=0.1863;c=0.6766]$. After 170 iterations the algorithm stopped convergent with objective function value of $1.16322$ and maximum constraint violation of $1.66\cdot 10^{-16}$. The parameters are found to be $$A=2.97\qquad B=0.99\qquad C=3.96\qquad M=0.88\qquad K=-0.88\qquad c=2.95$$ As you can see the calculated parameters are very close to the original parameters even though random error is added to the data. The final grap shows the data (red) and the calculated fitting function (blue) calculated

-------------------ADDITION OF CODE------------------------

Please note that $A\rightarrow x(1)$, $B\rightarrow x(2)$, $C\rightarrow x(3)$, $M\rightarrow x(4)$, $K\rightarrow x(5)$, $c\rightarrow x(6)$

First the function file for objective function (save it as a seperate file as objfun.m)

function f = objfun( x )
global n
global Y_er X
sum=0;
for i=1:2*n
    if X(i)<=x(6)
        f1=-sqrt(x(1)^2-(X(i)-x(2))^2)+x(3);
        sum=sum+(Y_er(i)-f1)^2;
    else
        f2=x(4)*X(i)+x(5);
        sum=sum+(Y_er(i)-f2)^2;
    end
end
f=sum;
end

As the second constraint equations (save it as a seperate file as constfun.m)

function [ c ceq ] = constfun( x )
c=[];
f1=-sqrt(x(1)^2-(x(6)-x(2))^2)+x(3);
f2=x(4)*x(6)+x(5);
g1=x(6)-x(2);
g2=x(4)*(sqrt(x(1)^2-(x(6)-x(2))^2));
ceq(1)=f1-f2;
ceq(2)=g1-g2;
ceq(3)=x(3)-3*x(2)-1;
end

And the final one is main file (save it as m.file with any name in the same folder with objfun.m and constfun.m)

clc
clear
global Y_er X n
global beta
%------------------- GENERATION OF SAMPLE POINTS -------------------------
% Generate 60 uniformly distributed random points between 0 and 3
n=60;
a=0;
b=3;
X1=a+(b-a)*rand(n,1);
% Generate y=-sqrt(A-(x-B)^2) + C function for random data
A=3;
B=1;
C=4;
for i=1:n
    Y1(i,1)=-sqrt(A^2-(X1(i)-B)^2)+C;
end
plot(X1,Y1,'k.');

% Generate errors with normal dist (mean 0 stdev 0.1)
Er=0.1*randn(n,1);
Y1_er=Y1+Er;
hold on;
plot(X1,Y1_er,'ro');

% Now generate 60 more samples between 3 and 6
clear a b
a=3;
b=6;
X2=a+(b-a)*rand(n,1);

% Generate y=m*x+k function for random data
m=(-B+a)/sqrt(A^2-(a-B)^2);
k=-sqrt(-(a-A-B)*(a+A-B))-a*m+C;
for i=1:n
    Y2(i,1)=m*X2(i)+k;
end
plot(X2,Y2,'b.');

% Generate errors with normal dist (mean 0 stdev 0.1)
Er=0.1*randn(n,1);
Y2_er=Y2+Er;
hold on;
plot(X2,Y2_er,'ro');



%Construct the final data set
X=[X1;X2];
Y_er=[Y1_er;Y2_er];
figure(2)
plot(X,Y_er,'b.');


%--------------------------------OPTIMIZATION---------------------------
x0=rand(6,1);
options = optimset('Display','iter-detailed','Algorithm','interior-point','TolFun',1e-12,'TolX',1e-12);
x = fmincon(@objfun,x0,[],[],[],[],[],[],@constfun,options);
figure(3)
hold on
for i=1:2*n
    if X(i)<=x(6)
        r=-sqrt(x(1)^2-(X(i)-x(2))^2)+x(3);
    else
        r=x(4)*X(i)+x(5);
    end
    plot(X(i),r,'b.')
end

PS: If you have your own data set put x-data in variable X and y-data in Y_er. In addition put n= your number of data points. Hope it works for you.

PS2: Please modify your constraint "center of arc on a given line" in "constfun.m" in ceq(3).

Related Question