Fitting circle to points with lsfitnonlinear

mustrhaz · **Joined:** Thu May 13, 2010 2:14 pm **Posts:** 2

Hello,

I'm working on a project and i've to find the center of circles.
I use OpenCv to find approximatively the center, then i use lsfitnonlineariteration function on the edges points (by canny), but the program crashes.
Thanks u.

Sources:

Code:

// the result of Hough algorithm on the picture is (161, 121) with a radius of 50 px.
   int centre_x = 161, centre_y = 121, r = 50;
   Mat image=imread("image.bmp", -1), edges;
   Canny(image, edges, 50, 100);
   vector<Point> pointsContour; // pointsContour contains the edges points (the superior part of the circle)
   for(int i=0; i<centre_y; i++)
   {
      uchar *ptr = edges.ptr<uchar>(i);
      for(int j=0; j<edges.cols; j++)
      {
         if(ptr[j])
         {
            pointsContour.push_back(Point(j, i));
         }
      }
   }
   int n = pointsContour.size(), m = 1, k = 2;
   ap::real_1d_array y, c;
   ap::real_2d_array x;
   lsfitreport rep;
   lsfitstate state;
   x.setlength(n, m);
   y.setlength(n);
   c.setlength(k);
   for(int i=0; i<n; i++)
   {
      x(i, 0) = pointsContour[i].x;
      y(i) = pointsContour[i].y;
   }
   c(0) = centre_x;
   c(1) = centre_y;

   lsfitnonlinearfgh(x, y, c, n, m, k, state);
   lsfitnonlinearsetcond(state, 0.1, 0.1, 0);
   while(lsfitnonlineariteration(state))
   {
      /*
         F = sqrt(r? - (x - xo)?) + yo
      */
      state.f = sqrt(carre(r) - carre(state.x(0) - state.c(0))) + state.c(1);
      if(state.needfg || state.needfgh)
      {
         /*
            DF/Dxo = (-2*xo + x) / (2 * sqrt(r? - (x - xo)?))
            DF/Dyo = 1;
         */
         state.g(0) = (- 2.0 * state.c(0) + state.x(0)) / 
            (2.0 * sqrt(carre(r) - carre(state.x(0) - state.c(0))));
         state.g(1) = 1.0;
      }
      if(state.needfgh)
      {
         /*
            D?F/Dxo? = -(-2*xo + x) / (4 * (r? - (x - xo)?)^(3/2))
         */
         state.h(0, 0) = - (-2.0 * state.c(0) + state.x(0)) / 
            (4.0 * (carre(r) - carre(state.x(0) - state.c(0))) * sqrt((carre(r) - carre(state.x(0) - state.c(0)))));
         state.h(0, 1) = 0.0;
         state.h(1, 0) = 0.0;
         state.h(1, 1) = 0.0;
      }
   }
   int info;
   lsfitnonlinearresults(state, info, c, rep);
   printf("%f\t%f\n", c(0), c(1));

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

You didn't provide additional info, but I suspect that your model tries to calculate sqrt() of negative number. That's because you choose incorrect model. You try to minimize vertical distance from point to the circle. But it is much better to minimize difference between (distance(point,circle_center)^2 and circle_radius^2. It is called "total least squares". Actually, it is the only way to solve your problem.

Read http://www.math.sunysb.edu/~scott/Book3 ... ircle.html - it contains some info on your problem. You can use MinLM unit to minimize E(a,b,r) function which were mentioned there.

The only drawback is that you have to write slightly more complex code to calculate F. LSFit unit does this work for you (LSFitNonlinear is just a wrapper around MinLM unit), but it can't be used to solve total least squares problems.

mustrhaz · **Joined:** Thu May 13, 2010 2:14 pm **Posts:** 2

Yes, the crash is due to sqrt function. I will use minlmcreatefgh function.
Thanks you very much.

alexander · **Joined:** Sun May 30, 2010 12:33 pm **Posts:** 2

I have a very similar problem. What I have is a list of 2D points to that I want to fit a circle or ellipse with unknown center and radius. The problem is not the mathematical model for, let's say, the circle, but how to apply it in alglib.
As far as I understand lsfitnonlinear is not capabale of solving this minimization problem. So I try to use the minlm unit, but I don't see how to pass all my points to minlmcreatefgh().

Let's say we want to minimize the function E(a,b,r) = SUM[ ((xi-a)^2 + (yi-b)^2 - r^2)^2 ] according to http://www.math.sunysb.edu/~scott/Book3 ... ircle.html. Could you give me a short example on how to apply this with alglib?

Sergey.Bochkanov · **Joined:** Fri May 07, 2010 7:06 am **Posts:** 927

I don't know what language you use, so I'll post it in pseudocode:

Code:

var
    State:  MinLMState;
    Rep:    MinLMReport;
    S:      array of Real;
    X, Y:   Real;
begin
    SetLength(S, 3);
    S[0]:=0; // A
    S[1]:=0; // B
    S[2]:=0; // R
    MinLMCreateFJ(3, N, S, State);
    MinLMSetCond(State, 0.0, 0.0, 0.001, 0);
    while MinLMIteration(State) do
    begin
        if State.NeedF then
            State.F:=SUM[ ((x[]-State.X[0])^2 + (y[]-State.X[1])^2 - State.X[2]^2)^2 ];
        if State.NeedFiJ then
        begin
            State.Fi[]:=(x[]-State.X[0])^2 + (y[]-State.X[1])^2 - State.X[2]^2;
            State.J[*,0]:=-2*(x[]-State.X[0]);
            State.J[*,1]:=-2*(y[]-State.X[1]);
            State.J[*,2]:=-2*State.X[2];
        end;
   end;
   MinLMResults(State, S, Rep);

I hope that I typed correct code for calculation of F=SUM(Fi), Fi and J={dFi/dA, dFi/dB, dFi/dR}. Please, recheck it before actually fitting something.

alexander · **Joined:** Sun May 30, 2010 12:33 pm **Posts:** 2

Thanks for your quick answer. The derivatives are correct ;-)
What I found is, that I have to initialize S[3] with a non-zero value, otherwise R will always be zero after the minimization, A is correct and B is wrong. If you keep that in mind and remember that r should be non-negative, everything works fine! Thank you!

forum.alglib.net

Forum rules

Fitting circle to points with lsfitnonlinear

Who is online