Introduction

This library is for determining the best-fitting 2D line, circle or rotated ellipse of a set of input points.

Take for example a set of 2D x,y points that closely but not accurately approximates a circle. Then there is a centre point and radius that represents the best circle that matches the points. This library will give you that point and radius. Similarly for a line: you get the gradient and y-axis intercept, and ellipse: the centre, major and minor axes and rotation.

Background

Best-fitting shapes have uses in engineering, surveying and computational geometry. There is plenty of literature out there on this, but very little reusable code seems to exist. This library offers a means to incorporate this ability into a client application. It does so using the least-squares method.

Using the code

The code is standard C++ with a dependency on the Boost ublas library for matrix operations. A Visual Studio 2008 solution file, and a Makefile are provided for building on Windows and Linux respectively. Boost needs to be pre-installed with the relevant directory paths setup in your build environment.

General overview

Use the factory class to make an instance of the BestFit-derived computation object.

Decide what std::ostream, if any, you want textual output of the computation to go to. e.g. Output to the console or to a std::stringstream, or none at all.

Create two instances of the BestFitIO struct, one for input parameters and one for output.

At the very least, two members need to be set for the input struct: numPoints is the number of points to be provided, and points being the starting address of a sequence of pairs of doubles being the x,y points.

The output struct can be left as default constructed, though see the Further input and output section below for additional options.

Call Compute() on the computation object, passing in the I/O structs.

If needed, get the results from the output BestFitIO instance, or rely simply on what was output to the ostream.

delete the factory-made object.

A simple example

Here is an example of computing the best-fitting circle of 10 test points that approximate a circle of unit radius centered at coordinate (100,200). The true centre of these points is then found to be (99.969,200.016) and the radius 0.977.

In addition, some quality statistics are output to the specified ostream, in this case std::cout. More or less output is set with the input BestFitIO::verbosity setting. No output is also possible with the one-parameter override of the factory constructor.

In addition to solving for the geometry parameters, the library can also calculate:

the adjusted points - for each input point, the equivalent point that lies on the best-fitting shape

the residuals - the corrections to the input points

In the simple case above, no adjusted points were output. However, if output.wantAdjustedObs=true, the output.points array will contain the points described in (1) above. Either allocate memory for output.points, or a trick to save memory is to set output.points = input.points. The input points will then be overwritten with the adjusted points. Also by setting output.wantResiduals=true, the output.residuals array is allocated and populated with the corrections described in (2) above. Both these predicates default to false.

The Demo application

As a testing tool and a graphical way of demonstrating the use of the library I wrote a demo application. This is dependent on MFC and the Dundas Ultimate Grid and therefore is Windows only. Note that by using the library it is then also dependent on Boost.

The demo application uses the best-fit code as a statically linked-in library and allows you to:

See the provisional and adjusted coordinate data, and a graphical view of the points;

Right-click to add, and left-click to remove points from the graphical view;

See the results of the adjustment dynamically change;

Open and save the points to file;

Generate test data of a random distribution of points approximating given line, circle or ellipse parameters.

Points of Interest

Command-line program

The library code can also be compiled as a stand-alone console program. This takes stream input from stdio. To build as such, all that you need to do is to undefine BEST_FIT_LIB in the preprocessor defines and set the build target as an .exe. For Linux the Makefile has this mode as the default target. When building with BEST_FIT_LIB defined there is an additional dependency on Boost program_options for parsing of command line switches. Command line usage is as follows:

Speed/performance metrics

Testing the speed of best-fit on a Core 2 Duo shows the following speed stats for increasing number of points.

No. of unknowns

100 pts

1 000 pts

10 000 pts

Line

2

<1ms

2ms

15ms

Circle

3

<1ms

4ms

60ms

Ellipse

5

1ms

10ms

124ms

There is a truly massive difference in matrix multiplication time between the Release and Debug builds of Boost uBlas. This is well documented on the Internet. I had unfortunately missed this in an earlier version of this article which showed drastically incorrect timings. This has been corrected and the above timings are with NDEBUG defined to disable the time consuming debug mode checks.

Algorithm

Linear least squares adjustment is the algorithm behind the code. Specifically the parametric case of the adjustment. Given provisional values for the unknowns, this least squares method solves for small adjustments that refine the provisionals which are then reused iteratively until the solution converges. Matrix operations are integral to this, unfortunately including matrix inversion, which means it will be slow for large matrices.

In the case of the best-fitting line, only y values are treated as observations. The x value is treated as invariable. The sum of the squares of the shortest vertical distance of each point to the line is minimised. This is the ordinary parametric case of the least squares adjustment.

In the case of the circle and the rotated ellipse, both the x and the y values are treated as observations. The sum of the squares of the shortest distance of each point to the shape is minimised. This is a special case of the general least squares adjustment known as the quasi-parametric case.

Both the parametric and quasi-parametric case are solved in the same way. The difference is in the formulation of the weight matrix and the meaning of the residuals. The residuals of the quasi-parametric case do not immediately equate to a useful geometric quantity as with the parametric case.

The solution vector, (or corrections to the provisional unknowns) is given by:

x=(ATPA)-1ATPℓ

and the vector of residuals, (or corrections to the observations) is given by:

v=Ax-ℓ

where A is known as the design matrix of the observations, P the weight matrix and ℓ the observation vector. The inverse that exists in the equation for the solution vector is what necessitates a matrix inversion and negatively affects the performance.

To formulate the design matrix the equations of the geometry shapes need to be linear. The equation of the straight line is already linear, but the equations for a circle and a rotated ellipse need to be linearised first by Taylor expansion.

Provisional values for the unknowns are first determined by approximation. Each coordinate observation and the provisional unknown values are substituted into the linearised equation to form one row of the A matrix. The ℓ matrix is the observed less provisional.

The solution vector x is then determined, and added to the provisional values to refine them. If the change was sufficiently small then the solution has converged and the provisional values are final. Otherwise the new provisionals are used to reformulate the matrices A, P and ℓ before iterating again.

In conclusion, this is one approach, not the only approach, and not necessarily the best approach. However this approach doesn't need to concern itself with differences between lines, circles and ellipses, and this makes it convenient for concise implementation as a general purpose library.

Yes, that is expected. Not so much a bug, but rather impossible to model a vertical line with the line equation y=mx+c because there is no y-intercept. However all you need to do is swap the y and x, solve, then the solved y-intercept is actually your x constant.

Thanks for creating this. A suggested improvement would be to add capabilities for non-linear best fit, such as polynomial best-fit or a Savitsky-Golay (did I spell this correctly?) filter. I would do it myself but it's very complicated.

Hi Alasdair Craig,
Thanks for your article.
It is useful and interesting.
In the other comments I read some interesting improvements,
but that does not give them the right to downvote this article.
I am sure you worked hard to have this method programmed and tested.
I am glad you want to share this with the world.
Who needs 10000 or more points ?
Often a few 100 points are sufficient.
So I do not understand why the other comments are so negative.
I think it is a valuable contribution to the code-project comunity.
I give: vote = max value;
Keep up the good work !
Greetings,
Frans.

Hi, thank you. It turns out I incorrectly did my performance measuring on the DEBUG build of boost::ublas and the release timings are massively faster. An update on this, and some bug fixes are in the pipeline.

Reasons for up/down-voting:
- Lack of description of the actually implemented algorithm

- matrix inversion is impractical or downright unusable for large dimensions. I'd never use it for dimensions greater than 10: it's numerically unstable, at O(N^3) way too slow for large sized problems, and it doesn't solve the problem when for various reasons the inverse cannot be determined. Better: iterative numerical solutions to solve large linear equation systems. Best: Express target function (sum of square point distances) as function of the unknown parameters (up to 5 for an ellipsis), and find the root(s) for its gradient function. Admittedly the later method delivers a non-linear equation system and is mathematically and numerically challenging. But it's at most 5 equations rather than O(N)!

- for lines it's known that the best approximate for a given set of points will pass through the mean point. So you really only need to find the optimum orientation. I. e. you can easily calculate one unknown with O(N) and then approximate the other.

Bonus point for going to the length of testing with large number of points even though it reveals the algorithms weakness.

another Bonus point for recognizing the lack of (reusable) functionality in this area, and offering a solution, even if it's not the most efficient.

Thank you for your comments, but this is a least-squares implementation. You're suggesting a different solution entirely for speed - that is not my intention. Matrix inversion is integral to the least-squares method and cannot be avoided. The article will be updated in due course to explain this.

Sorry, I haven't been clear. I do understand what you did, I did some similar things myself over the years.

The misunderstanding is that the least squares method does not in fact prescribe the exact way of how to determine the solution. Text books that introduce the method often like to provide an easy example that can be solved directly, with linear algebra. But that doesn't mean you can't use a different solver! Nor is the direct solution (involving matrix inversion) considered an integral part of the least squares method.

GOTOs are a bit like wire coat hangers: they tend to breed in the darkness, such that where there once were few, eventually there are many, and the program's architecture collapses beneath them. (Fran Poretto)

Both of those are good reasons for down-voting the article, but 1 seems harsh when the range is 1 to 5. The article itself should be nearly unintelligible to deserve a 1. (Based on votes on my own articles, this doesn't seem to be a consensus. Just my opinion.)

Not interested... don't vote.
Not fantastic, clear, concise... loss of one.
Inefficient... loss of one.
Code itself is barely readable... loss of one.
If you know of clearer, more efficient algorithms... loss of one.
If I agreed with just your written assessment, I think I'd vote 3 or 4. Maybe with egregious inefficiency, a 2.

I don't disagree with you, just don't have the interest to look into it in more detail at this time, so I'm not going to vote at all.

The algorithm is the parametric case of the least-squares adjustment. In time I'll update the article with further information. The time taken is almost solely due to the matrix inversion (I should have explained that). I agree it's unacceptable. Cholesky inversion should massively speed it up and this is a planned future enhancement. Boost unfortunately doesn't natively support Cholesky.