Linear least squares circle fit in C or C++

Question

0.00/5 (No votes)

See more:

Anyone aware of C/C++ code to fit a circle to a set of 2d datapoints using the linear least squares method as described by I.D.Coope?

Ron H.

Posted 19-Jul-11 0:03am

Ron H.

Add a Solution

3 solutions

Solution 1

Unsurprisingly, Google has lots of suggestions[^].

Posted 19-Jul-11 0:45am

Richard MacCutchan

Comments

Ron H. 19-Jul-11 7:17am

Also unsurprisingly Richard, most are irrelevant! While I didn't follow-up on all 475,000 hits ( now of course my question here is number 1) most address either "linear least squares fit" or "least squares circle fitting" ( non-linear ). Coope's method applies linear methods to the fitting of a circle to 2d data points. His solution seems to be much more tolerant of outliers thus providing a better "geometric" fit than most algebraic fits. Yes, Google is your friend and I am not sure what we did without it but sometimes a good forum of knowledgeable people provides a better result!

Ron H.

Richard MacCutchan 19-Jul-11 7:24am

I hear what you say, but this is a technical Q&A forum and less suitable for general research questions such as yours. I appreciate that you are one of the exceptions to the rule, but we do get a lot of questions here of the form "how can I ...", where the poster has not done a basic Google search of the subject.

Add a Solution

Add your solution here

Treat my content as plain text, not as HTML

Preview 0

…

Existing Members

Sign in to your account

...or Join us

Download, Vote, Comment, Publish.

Your Email
Password
Forgot your password?

Your Email
This email is in use. Do you need your password?
Optional Password

I have read and agree to the Terms of Service and Privacy Policy
Please subscribe me to the CodeProject newsletters

When answering a question please:

Read the question carefully.
Understand that English isn't everyone's first language so be lenient of bad spelling and grammar.
If a question is poorly phrased then either ask for clarification, ignore it, or edit the question and fix the problem. Insults are not welcome.
Don't tell someone to read the manual. Chances are they have and don't get it. Provide an answer or move on to the next question.

Let's work to help developers, not make them feel stupid.

This content, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

YvesDaoust · Accepted Answer · 2011-07-19T04:58:00

Solution 3

You can do-it-yourself from this:

- use a 3x3 array A and a 3x1 vector B and accumulate the normal least-squares equations, using the formulas in the paper (this is an easy task);

- solve the 3x3 system using this function (this is the Cholesky method for a 3x3 symmetric positive defined matrix):

C#

bool Cholesky3x3(double* A, double* B)
{
#define A(i,j) A[i + j * 3]
    double Sum;
    double Diagonal[3];
    Sum= A(0,0);
    if (Sum <= 0.f) return false;
    Diagonal[0]= Sqrt(Sum);
    Sum= A(0,1);
    A(1,0)= Sum / Diagonal[0];
    Sum= A(0,2);
    A(2,0)= Sum / Diagonal[0];
    Sum= A(1,1) - A(1,0) * A(1,0);
    if (Sum <= 0.f) return false;
    Diagonal[1]= Sqrt(Sum);
    Sum= A(1,2) - A(1,0) * A(2,0);
    A(2,1)= Sum / Diagonal[1];
    Sum= A(2,2) - A(2,1) * A(2,1) - A(2,0) * A(2,0);
    if (Sum <= 0.f) return false;
    Diagonal[2]= Sqrt(Sum);
    Sum= B[0];
    B[0]= Sum / Diagonal[0];
    Sum= B[1] - A(1,0) * B[0];
    B[1]= Sum / Diagonal[1];
    Sum= B[2] - A(2,1) * B[1] - A(2,0) * B[0];
    B[2]= Sum / Diagonal[2];
    Sum= B[2];
    B[2]= Sum / Diagonal[2];
    Sum= B[1] - A(2,1) * B[2];
    B[1]= Sum / Diagonal[1];
    Sum= B[0] - A(1,0) * B[1] - A(2,0) * B[2];
    B[0]= Sum / Diagonal[0];
    return true;
}

Solutions in B. Compute Xc, Yc, R, you are done. Requires no library.

Posted 19-Jul-11 4:58am

YvesDaoust

Updated 19-Jul-11 5:07am

v7

Comments

Ron H. 19-Jul-11 11:08am

Thanks! I appreciate your response!

Ron H.

YvesDaoust 19-Jul-11 16:26pm

You can rate it then :)

YvesDaoust 20-Jul-11 4:31am

Thx

Swatieson 1-Mar-15 18:12pm

What are the formulas for the arrays? Apparently they are not in the Internet.

YvesDaoust 2-Mar-15 2:47am

Lookup http://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Derivation_of_the_normal_equations

Member 11511083 9-Mar-15 16:34pm

Thanks for your solution YvesDaoust. I'm having trouble putting together all the pieces. I follow the 3x3 Cholesky decomposition. But I don't understand where this 3x3 matrix comes from. I have read the link you've posted several times.

Considering I have N points, I see instead that I have a 3xN matrix, rather than 3x3.

How dissimilar is your solution to the one posted in Solution 2.

Thank you.

Member 11511083 9-Mar-15 16:43pm

I think I'm starting to get there, that my system of points X (3xN) times it's transpose must always be a square matrix 3x3... been a while since I've done linear algebra.

YvesDaoust 9-Mar-15 16:44pm

Read this: http://mathworld.wolfram.com/NormalEquation.html

Member 11511083 9-Mar-15 17:50pm

Thank you. What I was missing was that if a matrix X is Nx3, then it's transpose times itself (and not the other way around) results in a 3x3 matrix.

Thanks again.

Member 12328488 18-May-16 9:15am

Hi I am trying this code and ı wrote this
int main(int argc, char** argv) {
double matris[9]={25,30,50,40,20,71,43,78,12};
double matris2[3]={40,25,50};
//double matris=25,matris2=50;
Cholesky3x3(matris,matris2);

return 0;
}
But promram no write circle equation

YvesDaoust 18-May-16 9:32am

Don't worry, this is normal.

Member 12328488 18-May-16 9:35am

O thanks :) but what should ı do for this

YvesDaoust 18-May-16 9:44am

Write a complete program. If you don't have the skills, ask someone around you.

Member 12328488 18-May-16 9:46am

Ok Thank you

Member 12328488 18-May-16 14:51pm

Hi, How I can write complete program for this algorithm Can you help me?

YvesDaoust 18-May-16 15:08pm

Sorry, no, I am busy.

Stefan_Lang · Accepted Answer · 2011-07-19T04:08:00

Solution 2

Here's a JavaScript implementation: Least-squares circle fitting

It shouldn't be too hard to translate.

Posted 19-Jul-11 4:08am

Stefan_Lang

Comments

Ron H. 19-Jul-11 11:07am

Thanks! I am not a Java programmer but this looks pretty straight forward!

Ron H.