Fitting Geometric Primitives to Points Using Nonlinear Least Squares
Tuesday, August 9th, 2011We were recently contacted by a customer looking for help on how to use NMath to fit geometric primitives to clouds of 2D points. The solution is to cast the problem as a minimization. In the space of the parameters which define the geometric object, minimize the residuals with respect to the 2D points. NMath provides .NET classes for solving nonlinear least squares problems such as this, using the Trust-Region or Levenberg-Marquardt methods.
For instance, let’s try fitting a circle to a set of points. A circle is defined by three parameters: a center (a,b), and a radius r. The circle is all points such that (x-a)^2 + (y – b)^2 = r^2. To setup the minimization problem, we first need to define a function which given a set of circle parameters, returns the residuals with respect to a cloud of 2D points. There are several methods for doing this in NMath. In the following C# code, we extend the abstract base class DoubleMultiVariableFunction, and override the Evaluate() method:
public class CircleFitFunction : DoubleMultiVariableFunction { public CircleFitFunction( DoubleVector x, DoubleVector y ) : base( 3, x.Length ) { if( x.Length != y.Length ) throw new Exception( "Unequal number of x,y values." ); X = x; Y = y; } public DoubleVector X { get; internal set; } public DoubleVector Y { get; internal set; } public override void Evaluate( DoubleVector parameters, ref DoubleVector residuals ) { // parameters of circle with center (a,b) and radius r double a = parameters[0]; double b = parameters[1]; double r = parameters[2]; for( int i = 0; i < X.Length; i++ ) { // distance of point from circle center double d = Math.Sqrt( Math.Pow( X[i] - a, 2.0 ) + Math.Pow( Y[i] - b, 2.0 ) ); residuals[i] = d - r; } } }
The constructor takes the set of x,y points to fit, which are stored on the class. Note that we call the base constructor with the dimensions of the domain and range of the function:
: base( 3, x.Length )
In this case, the function maps the 3-dimensional space of the circle parameters to the x.Length-dimension space of the residuals.
Next we override the Evaluate() method, which takes an input vector and output vector passed by reference. Our implementation computes the residual for each point with respect to a circle defined by the given parameters. We calculate the distance, d, of each point from the circle center; the residual is then equal to d – r. The nonlinear least squares method will minimize the L2 norm of this function.
To demonstrate the fitting process, let’s first start with a set of points generated from a circle of known parameters. For example, this C# code creates 20 x,y points evenly spaced around a circle with center (0.5, 0.25) and radius 2, plus some added noise:
int n = 20; DoubleVector x = new DoubleVector( n ); DoubleVector y = new DoubleVector( n ); double a = 0.5; double b = 0.25; double r = 2; RandGenUniform noise = new RandGenUniform( -.75, .75 ); for( int i = 0; i < n; i++ ) { double t = i * 2 * Math.PI / n; x[i] = a + r * Math.Cos( t ) + noise.Next(); y[i] = b + r * Math.Sin( t ) + noise.Next(); }
Now that we have our function defined and some points to fit, performing the minimization takes only a few lines of code. This C# code finds the minimum and prints the solution and final residual:
CircleFitFunction f = new CircleFitFunction( x, y ); TrustRegionMinimizer minimizer = new TrustRegionMinimizer(); DoubleVector start = new DoubleVector( "0.1 0.1 0.1" ); DoubleVector solution = minimizer.Minimize( f, start ); Console.WriteLine( "solution = " + solution ); Console.WriteLine( "final residual = " + minimizer.FinalResidual );
Sample output:
solution = [ 0.446122611523468 0.175483486509012 1.93511389538286 ] final residual = 1.62414259527024
Starting from point (0.1, 0.1, 0.1) in the parameter space of the circle, we minimized the sum of the squared residuals with respect to the (x,y) points. In the run above, the minimum found was center (0.45, 0.18) and radius 1.9, close to the actual parameters which generated our (noisy) data. To visually inspect the fit, we can use NMath with the Microsoft Chart Controls for .NET.
Now that we’re sure the minimization is working as desired, let’s try it again with some random points:
int n = 10; RandGenUniform rng = new RandGenUniform( -1, 2 ); DoubleVector x = new DoubleVector( n, rng ); DoubleVector y = new DoubleVector( n, rng ); CircleFitFunction f = new CircleFitFunction( x, y ); TrustRegionMinimizer minimizer = new TrustRegionMinimizer(); DoubleVector start = new DoubleVector( "0.1 0.1 0.1" ); DoubleVector solution = minimizer.Minimize( f, start );
Again, plotting the solution:
In some runs, the fit may seem counter-intuitive. For example:
But if you think about it, this sort of fit makes perfect sense, given how we defined our fit function. We’re asking the minimizer to minimize the residuals with respect to the circle perimeter, and allowing it vary the circle center and radius however it can to achieve that goal. In an extreme case, if our points fall approximately along a line, the best fit will be a circle with a very large radius, so that the circle perimeter is nearly linear as it passes through the points.
In this case, the fitted circle has center (428.0, -757.8) and radius 870.6. If we wish to avoid solutions such as these, we could define our fit function differently. For example, we might constrain the circle center to be the center of mass of the points, and vary only the radius.
Obviously, a similar technique can be used to fit other geometric primitives, though as the complexity of the shape increases, so does the complexity of the fit function.
Ken
