Difference between revisions of "Least-squares circle fitting"

From JSXGraph Wiki
Jump to navigationJump to search
 
(21 intermediate revisions by the same user not shown)
Line 1: Line 1:
 +
This is an implementation of the linear least-squares algorithm by Coope (1993) for circle fitting.
 
<jsxgraph width="600" height="600">
 
<jsxgraph width="600" height="600">
 
var brd = JXG.JSXGraph.initBoard('jxgbox',{boundingbox:[-5,5,5,-5], keepaspectratio:true});
 
var brd = JXG.JSXGraph.initBoard('jxgbox',{boundingbox:[-5,5,5,-5], keepaspectratio:true});
 +
var i, p = [], angle, co, si, delta = 0.8;
  
 +
// Random points are constructed which lie roughly on a circle
 +
// of radius 4 having the origin as center.
 +
// delta*0.5 is the maximal distance in x- and y- direction of the random
 +
// points from the circle line.
 +
brd.suspendUpdate();
 +
for (i=0;i<100;i++) {
 +
  angle = Math.random()*2*Math.PI;
 +
 +
  co = 4*Math.cos(angle)+delta*(Math.random()-0.5);
 +
  si = 4*Math.sin(angle)+delta*(Math.random()-0.5);
 +
  p.push(brd.create('point',[co, si], {withLabel:false}));
 +
}
 +
brd.unsuspendUpdate();
 +
 +
// Having constructed the points, we can fit a circle
 +
// through the point set, consisting of n points.
 +
// The (n times 3) matrix consists of
 +
//  x_1, y_1, 1
 +
//  x_2, y_2, 1
 +
//      ...
 +
//  x_n, y_n, 1
 +
// where x_i, y_i is the position of point p_i
 +
// The vector y of length n consists of
 +
//    x_i*x_i+y_i*y_i
 +
// for i=1,...n.
 +
var M = [], y = [], MT, B, c, z, n;
 +
n = p.length;
 +
for (i=0;i<n;i++) {
 +
    M.push([p[i].X(), p[i].Y(), 1.0]);
 +
    y.push(p[i].X()*p[i].X() + p[i].Y()*p[i].Y());
 +
}
 +
 +
// Now, the general linear least-square fitting problem
 +
//    min_z || M*z - y||_2^2
 +
// is solved by solving the system of linear equations
 +
//    (M^T*M) * z = (M^T*y)
 +
// with Gauss elimination.
 +
MT = JXG.Math.transpose(M);
 +
B = JXG.Math.matMatMult(MT, M);
 +
c = JXG.Math.matVecMult(MT, y);
 +
z = JXG.Math.Numerics.Gauss(B, c);
 +
 +
// Finally, we can read from the solution vector z the coordinates [xm, ym] of the center
 +
// and the radius r and draw the circle.
 +
var xm = z[0]*0.5;
 +
var ym = z[1]*0.5;
 +
var r = Math.sqrt(z[2]+xm*xm+ym*ym);
 +
 +
brd.create('circle',[ [xm,ym], r]);
 
</jsxgraph>
 
</jsxgraph>
 +
===The underlying JavaScript code===
 +
<source lang="javascript">
 +
var brd = JXG.JSXGraph.initBoard('jxgbox',{boundingbox:[-5,5,5,-5], keepaspectratio:true});
 +
var i, p = [], angle, co, si, delta = 0.8;
 +
 +
// Random points are constructed which lie roughly on a circle
 +
// of radius 4 having the origin as center.
 +
// delta*0.5 is the maximal distance in x- and y- direction of the random
 +
// points from the circle line.
 +
brd.suspendUpdate();
 +
for (i=0;i<100;i++) {
 +
  angle = Math.random()*2*Math.PI;
 +
 +
  co = 4*Math.cos(angle)+delta*(Math.random()-0.5);
 +
  si = 4*Math.sin(angle)+delta*(Math.random()-0.5);
 +
  p.push(brd.create('point',[co, si], {withLabel:false}));
 +
}
 +
brd.unsuspendUpdate();
 +
 +
// Having constructed the points, we can fit a circle
 +
// through the point set, consisting of n points.
 +
// The (n times 3) matrix consists of
 +
//  x_1, y_1, 1
 +
//  x_2, y_2, 1
 +
//      ...
 +
//  x_n, y_n, 1
 +
// where x_i, y_i is the position of point p_i
 +
// The vector y of length n consists of
 +
//    x_i*x_i+y_i*y_i
 +
// for i=1,...n.
 +
var M = [], y = [], MT, B, c, z, n;
 +
n = p.length;
 +
for (i=0;i<n;i++) {
 +
    M.push([p[i].X(), p[i].Y(), 1.0]);
 +
    y.push(p[i].X()*p[i].X() + p[i].Y()*p[i].Y());
 +
}
 +
 +
// Now, the general linear least-square fitting problem
 +
//    min_z || M*z - y||_2^2
 +
// is solved by solving the system of linear equations
 +
//    (M^T*M) * z = (M^T*y)
 +
// with Gauss elimination.
 +
MT = JXG.Math.transpose(M);
 +
B = JXG.Math.matMatMult(MT, M);
 +
c = JXG.Math.matVecMult(MT, y);
 +
z = JXG.Math.Numerics.Gauss(B, c);
 +
 +
// Finally, we can read from the solution vector z the coordinates [xm, ym] of the center
 +
// and the radius r and draw the circle.
 +
var xm = z[0]*0.5;
 +
var ym = z[1]*0.5;
 +
var r = Math.sqrt(z[2]+xm*xm+ym*ym);
 +
 +
brd.create('circle',[ [xm,ym], r]);
 +
</source>
 +
 +
===References===
 +
* Coope, I.D., ''Circle fitting by linear and nonlinear least squares'', Journal of Optimization Theory and Applications Volume 76, Issue 2, New York: Plenum Press, February 1993
  
 
[[Category:Examples]]
 
[[Category:Examples]]
[[Catgegory:Statistics]]
+
[[Category:Statistics]]

Latest revision as of 22:45, 1 January 2011

This is an implementation of the linear least-squares algorithm by Coope (1993) for circle fitting.

The underlying JavaScript code

var brd = JXG.JSXGraph.initBoard('jxgbox',{boundingbox:[-5,5,5,-5], keepaspectratio:true});
var i, p = [], angle, co, si, delta = 0.8;

// Random points are constructed which lie roughly on a circle
// of radius 4 having the origin as center.
// delta*0.5 is the maximal distance in x- and y- direction of the random
// points from the circle line.
brd.suspendUpdate();
for (i=0;i<100;i++) {
  angle = Math.random()*2*Math.PI;

  co = 4*Math.cos(angle)+delta*(Math.random()-0.5);
  si = 4*Math.sin(angle)+delta*(Math.random()-0.5);
  p.push(brd.create('point',[co, si], {withLabel:false}));
}
brd.unsuspendUpdate();

// Having constructed the points, we can fit a circle 
// through the point set, consisting of n points.
// The (n times 3) matrix consists of
//   x_1, y_1, 1
//   x_2, y_2, 1
//      ...
//   x_n, y_n, 1
// where x_i, y_i is the position of point p_i
// The vector y of length n consists of
//    x_i*x_i+y_i*y_i 
// for i=1,...n.
var M = [], y = [], MT, B, c, z, n;
n = p.length;
for (i=0;i<n;i++) {
    M.push([p[i].X(), p[i].Y(), 1.0]);
    y.push(p[i].X()*p[i].X() + p[i].Y()*p[i].Y());
}

// Now, the general linear least-square fitting problem
//    min_z || M*z - y||_2^2
// is solved by solving the system of linear equations
//    (M^T*M) * z = (M^T*y)
// with Gauss elimination.
MT = JXG.Math.transpose(M);
B = JXG.Math.matMatMult(MT, M);
c = JXG.Math.matVecMult(MT, y);
z = JXG.Math.Numerics.Gauss(B, c);

// Finally, we can read from the solution vector z the coordinates [xm, ym] of the center
// and the radius r and draw the circle.
var xm = z[0]*0.5;
var ym = z[1]*0.5;
var r = Math.sqrt(z[2]+xm*xm+ym*ym);

brd.create('circle',[ [xm,ym], r]);

References

  • Coope, I.D., Circle fitting by linear and nonlinear least squares, Journal of Optimization Theory and Applications Volume 76, Issue 2, New York: Plenum Press, February 1993