// 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));

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