1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ The point the tangent line and circle meet is pre-determined by you? $\endgroup$ – AnilB May 22 '13 at 16:12
  • $\begingroup$ @AnilBaseski: I know only the set of points to be fitted and the line where the center of the circle lies. The optimal point where the two curves meet must be found. $\endgroup$ – digital-Ink May 22 '13 at 16:17
  • $\begingroup$ "..a circle whose center lies on...a half-line tangent to one of the endpoints of the arc circle" Are you sure this is what you mean? If the center lies on some (half-)line $\ell$, then $\ell$ cannot be tangent to the circle anywhere. $\endgroup$ – alex.jordan May 25 '13 at 6:36
  • $\begingroup$ Ah I see. You are asking for a "hook" shaped curve, the union of a circular arc and a line segment extending smoothly from the arc. Your current wording made me think the center was supposed to be at the intersection of some like and some half-line. $\endgroup$ – alex.jordan May 25 '13 at 6:38
  • $\begingroup$ @alex.jordan: You are right. I made an edit to the op such that the question is much clearer. $\endgroup$ – digital-Ink May 25 '13 at 6:46
1
+50
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ Can you also provide the Matlab code for testing? $\endgroup$ – digital-Ink May 26 '13 at 16:20
  • $\begingroup$ Also, I see that the objective function to be minimized has been constructed using only the vertical ($y$) displacements of the point w.r.t. the curve, and that you assumed that the curve "starts" (w.r.t. $x$) with the circle and then continues with the tangent, just like if it were the graph of a function. Please note the curve may very well not be the graph of a function. $\endgroup$ – digital-Ink May 26 '13 at 16:33
  • $\begingroup$ Maybe I can send m files with mail? Or do you prefer that I post it? Do you have a set of sample points? $\endgroup$ – AnilB May 26 '13 at 16:49
  • $\begingroup$ I guess you may integrate it with additional function and another break-point? $\endgroup$ – AnilB May 26 '13 at 16:50
  • $\begingroup$ I do have a set of points, but it is quite large. It is your choice how you share the code. $\endgroup$ – digital-Ink May 26 '13 at 16:52
1
$\begingroup$

Such curves have $5$ real parameters:

  • $(s,t)$, the point where the arc meets the line
  • $r$, the radius of the arc
  • $\theta$, the angle of the arc
  • $d$, the length of the line segment

And one discrete parameter that chooses the orientation of the hook.

You don't lose anything by letting the line extend forever and the arc complete its circle. It seems to me that doing so only increases your ability to obtain a smaller distance from the data to the curve. So then you have three real parameters and one discrete parameter.

Theoretically, if you sum the squared distances from your data points to the curve with these parameters, and use calculus to minimize, you will have your result.

But there are questions: what will be your measurement of distance? In least-squares regression, the distance is measured as the vertical separation. Do you instead want to minimize direct distance from data points to the curve?

$\endgroup$
  • $\begingroup$ I believe that the sign of $\theta$ is enough to give the orientation of the arc. Also, I don't think that $d$ (the lenght of the lien segment) is so important (it can be easily found or estimated from the other parameters and the given points, at the last stage). Also, I don't think that the center of the circle is uniquely determined from $(s,t)$ and $r$ (they may be no point, one point or two points at the distance $r$ from $(s,t)$ that lies on a given line. The considered distance clearly should be a direct distance between the curve and the given points, not the vertical separation. $\endgroup$ – digital-Ink May 25 '13 at 7:14
  • $\begingroup$ Yes, $d$ does not matter - that's what I mean by extending the line forever. And $\theta$ is also removed from consideration by letting the arc become a full circle. You are left with $(s,t)$, $r$, and the orientation to try calculus on, and as you rightly point out, another discrete parameter. Of course, there is a domain restriction in that the distance between $(s,t)$ and $\ell$ must be at least $r$. Then there is another discrete (two-option) parameter for a choice of center, which can be ordered based on a counterclockwise rotation from $(s,t)$. Can you give a sample data set? $\endgroup$ – alex.jordan May 25 '13 at 15:54
  • $\begingroup$ How should I give the sample data? As a picture or as coordinates set? $\endgroup$ – digital-Ink May 26 '13 at 16:36
  • $\begingroup$ With coordinates. With four points in general position, any three form a circle, and then there are two ways to connect the line tangentially. So give at least five points but not too many. $\endgroup$ – alex.jordan May 26 '13 at 16:46
  • $\begingroup$ Unfortunatelly, I have quite a large set of sample points (from a section to a laser-scanned object) $\endgroup$ – digital-Ink May 26 '13 at 16:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.