# Difference between revisions of "Least-squares circle fitting"

From JSXGraph Wiki

Jump to navigationJump to searchA WASSERMANN (talk | contribs) |
A WASSERMANN (talk | contribs) |
||

(5 intermediate revisions by the same user not shown) | |||

Line 4: | Line 4: | ||

var i, p = [], angle, co, si, delta = 0.8; | 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. | // of radius 4 having the origin as center. | ||

// delta*0.5 is the maximal distance in x- and y- direction of the random | // delta*0.5 is the maximal distance in x- and y- direction of the random | ||

Line 25: | Line 25: | ||

// ... | // ... | ||

// x_n, y_n, 1 | // x_n, y_n, 1 | ||

− | // where x_i, y_i is the | + | // where x_i, y_i is the position of point p_i |

// The vector y of length n consists of | // The vector y of length n consists of | ||

// x_i*x_i+y_i*y_i | // x_i*x_i+y_i*y_i | ||

Line 36: | Line 36: | ||

} | } | ||

− | // Now, the general | + | // Now, the general linear least-square fitting problem |

// min_z || M*z - y||_2^2 | // min_z || M*z - y||_2^2 | ||

− | // is solved by solving the system of linear equations | + | // is solved by solving the system of linear equations |

// (M^T*M) * z = (M^T*y) | // (M^T*M) * z = (M^T*y) | ||

// with Gauss elimination. | // with Gauss elimination. | ||

Line 59: | Line 59: | ||

var i, p = [], angle, co, si, delta = 0.8; | 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. | // of radius 4 having the origin as center. | ||

// delta*0.5 is the maximal distance in x- and y- direction of the random | // delta*0.5 is the maximal distance in x- and y- direction of the random | ||

Line 80: | Line 80: | ||

// ... | // ... | ||

// x_n, y_n, 1 | // x_n, y_n, 1 | ||

− | // where x_i, y_i is the | + | // where x_i, y_i is the position of point p_i |

// The vector y of length n consists of | // The vector y of length n consists of | ||

// x_i*x_i+y_i*y_i | // x_i*x_i+y_i*y_i | ||

Line 91: | Line 91: | ||

} | } | ||

− | // Now, the general | + | // Now, the general linear least-square fitting problem |

// min_z || M*z - y||_2^2 | // min_z || M*z - y||_2^2 | ||

− | // is solved by solving the system of linear equations | + | // is solved by solving the system of linear equations |

// (M^T*M) * z = (M^T*y) | // (M^T*M) * z = (M^T*y) | ||

// with Gauss elimination. | // with Gauss elimination. | ||

Line 101: | Line 101: | ||

z = JXG.Math.Numerics.Gauss(B, c); | z = JXG.Math.Numerics.Gauss(B, c); | ||

− | // Finally, we can read | + | // Finally, we can read from the solution vector z the coordinates [xm, ym] of the center |

− | // and the radius r | + | // and the radius r and draw the circle. |

var xm = z[0]*0.5; | var xm = z[0]*0.5; | ||

var ym = z[1]*0.5; | var ym = z[1]*0.5; |

## 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