Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Fit circle to 3 points

Subject: Fit circle to 3 points

From: Peter Bone

Date: 6 Mar, 2008 10:59:02

Message: 1 of 9

I've written a function to calculate the circle that passes
through any non-collinear set of 3 points. There are a
number of examples on the file exchange that fit a circle to
any number of points using optimization techniques but none
that do just the 3 point problem analytically.
The problem is that my code (below) is quite long because of
a number of checks required to prevent divide by zeros. This
is because it uses the gradients of line segments in
Cartesian coordinates. Does anyone have a cleverer method of
doing this? If not then I'll upload my method to the file
exchange.

Peter


p1 = rand(1,2);
p2 = rand(1,2);
p3 = rand(1,2);

[c r] = calc_circle(p1, p2, p3);

figure(1)
cla
axis equal
hold on
if r == -1
    disp('COLLINEAR')
else
    rectangle('Position', [c(1)-r,c(2)-r,2*r,2*r],
'Curvature', [1,1], 'EdgeColor', 'g')
end
plot(p1(1), p1(2), '*')
plot(p2(1), p2(2), '*')
plot(p3(1), p3(2), '*')


function [centre radius] = calc_circle(pt1, pt2, pt3)

delta_a = pt2 - pt1;
delta_b = pt3 - pt2;

ax_is_0 = abs(delta_a(1)) <= 0.000000001;
bx_is_0 = abs(delta_b(1)) <= 0.000000001;

% check whether both lines are vertical - collinear
if (ax_is_0 && bx_is_0)
    centre = [0 0];
    radius = -1;
    return
end

% make sure delta gradients are not vertical
% re-arrange points to change deltas
if (ax_is_0)
    [centre radius] = calc_circle(pt1, pt3, pt2);
    return
end
if (bx_is_0)
    [centre radius] = calc_circle(pt2, pt1, pt3);
    return
end

grad_a = delta_a(2) / delta_a(1);
grad_b = delta_b(2) / delta_b(1);

% check whether the given points are collinear
if (abs(grad_a-grad_b) <= 0.000000001)
    centre = [0 0];
    radius = -1;
    return
end

% swap grads and points if grad_a is 0
if abs(grad_a) <= 0.000000001
    tmp = grad_a;
    grad_a = grad_b;
    grad_b = tmp;
    tmp = pt1;
    pt1 = pt3;
    pt3 = tmp;
end

% calculate centre - where the lines perpendicular to the
centre of
% segments a and b intersect.
centre(1) = ( grad_a*grad_b*(pt1(2)-pt3(2)) +
grad_b*(pt1(1)+pt2(1)) - grad_a*(pt2(1)+pt3(1)) ) /
(2*(grad_b-grad_a));
centre(2) = ((pt1(1)+pt2(1))/2 - centre(1)) / grad_a +
(pt1(2)+pt2(2))/2;

% calculate radius
radius = norm(centre - pt1);

Subject: Fit circle to 3 points

From: Jos

Date: 6 Mar, 2008 11:44:03

Message: 2 of 9

"Peter Bone" <peterbone@hotmail.com> wrote in message
<fqoipl$blh$1@fred.mathworks.com>...
> I've written a function to calculate the circle that
passes ...

Another solution which I found in my own archive:

function [xy,r] = points2circle(A,B,C)
% POINTS2CIRCLE - find circle given three points
%
% [XY,R] = POINTS2CIRCLE(A,B,C) finds the center point XY
and the radius R
% of the circle passing through the three points A, B, and C
%
% POINTS2CIRCLE(P) in which P = [A ; B ; C] can also be
used. P is a 3-by-2
% matrix.
%
% Example:
% POINTS2CIRCLE([1 0],[-4 3],[3,2])

if nargin==1,
    P = A ;
    if size(P,1) ~=3 | size(P,2) ~=2,
        error ('A single input should be a 3-by-2 matrix') ;
    end
elseif ~isequal(numel(A),numel(B),numel(C),2),
    error('The three points should all have 2 elements.')
else
    P = [A(:) B(:) C(:)] .' ;
end

% matrix
M = [...
        1 1 1 1 ; ...
        (P(1,1).^2 + P(1,2).^2) P(1,1) P(1,2) 1 ; ...
        (P(2,1).^2 + P(2,2).^2) P(2,1) P(2,2) 1 ; ...
        (P(3,1).^2 + P(3,2).^2) P(3,1) P(3,2) 1 ...
    ] ;

M11 = local_minordet(M,1,1) ;
if M11==0,
    xy = [] ;
    r = [] ;
    warning('No solution! Points may be on a straight line.') ;
else
    xy(1) = 0.5 * (local_minordet(M,1,2) ./ M11) ;
    xy(2) = -0.5 * (local_minordet(M,1,3) ./ M11) ;
    r = sqrt(xy(1).^2 + xy(2).^2 + (local_minordet(M,1,4) ./
M11)) ;
end


function md = local_minordet(M,i,j) ;
% minor determinant
M(i,:)=[] ;
M(:,j)=[] ;
md = det(M) ;

% END OF FUNCTION

beware of line wraps ...

hth
Jos

Subject: Fit circle to 3 points

From: Peter Bone

Date: 6 Mar, 2008 12:02:02

Message: 3 of 9

"Jos " <DELjos@jasenDEL.nl> wrote in message
<fqole3$edg$1@fred.mathworks.com>...
> "Peter Bone" <peterbone@hotmail.com> wrote in message
> <fqoipl$blh$1@fred.mathworks.com>...
> > I've written a function to calculate the circle that
> passes ...
>
> Another solution which I found in my own archive:

Thanks. That method looks interesting but I'll stick with
mine because it converts more easily to C and is almost 4
times faster on my PC.

Peter

Subject: Fit circle to 3 points

From: Yumnam Kirani Singh

Date: 6 Mar, 2008 12:14:37

Message: 4 of 9

You can do it just within 4/5 lines. There is a direct formula to find the centre of circle that passes through three points. You can visit my webpage http://www.geocities.com/kiranisingh/center.html and use the formula 2.

Subject: Fit circle to 3 points

From: Jos

Date: 6 Mar, 2008 12:37:02

Message: 5 of 9

"Peter Bone" <peterbone@hotmail.com> wrote in message
<fqomfq$6ub$1@fred.mathworks.com>...
> "Jos " <DELjos@jasenDEL.nl> wrote in message
> <fqole3$edg$1@fred.mathworks.com>...
> > "Peter Bone" <peterbone@hotmail.com> wrote in message
> > <fqoipl$blh$1@fred.mathworks.com>...
> > > I've written a function to calculate the circle that
> > passes ...
> >
> > Another solution which I found in my own archive:
>
> Thanks. That method looks interesting but I'll stick with
> mine because it converts more easily to C and is almost 4
> times faster on my PC.
>
> Peter

Fair enough. I think I will put my function of the FEX as
well. What were the timings you btw?

Best,
Jos

Subject: Fit circle to 3 points

From: Roger Stafford

Date: 6 Mar, 2008 21:59:02

Message: 6 of 9

"Peter Bone" <peterbone@hotmail.com> wrote in message <fqoipl$blh
$1@fred.mathworks.com>...
> I've written a function to calculate the circle that passes
> through any non-collinear set of 3 points. There are a
> number of examples on the file exchange that fit a circle to
> any number of points using optimization techniques but none
> that do just the 3 point problem analytically.
> The problem is that my code (below) is quite long because of
> a number of checks required to prevent divide by zeros. This
> is because it uses the gradients of line segments in
> Cartesian coordinates. Does anyone have a cleverer method of
> doing this? If not then I'll upload my method to the file
> exchange.
>
> Peter
> ...........

  Just in case it is of interest to you, Peter, here is a solution to your problem
in three-dimensional space.

  Let p1 = (x1,y1,z1), p2 = (x2,y2,z2), and p3 = (x3,y3,z3) be three points in
3D space, each represented by a three-element (row or column) vector. Then
a circle through these three points has its center at the location given by
vector c with radius r as given in the following:

 t = p2-p1; u = p3-p1; v = p3-p2;
 w = cross(t,u);
 c = p1+(dot(t,t)*dot(u,v)*u-dot(u,u)*dot(t,v)*t)/(2*dot(w,w));
 r = 1/2*norm(t)*norm(u)*norm(v)/norm(w);

  However, the following is probably a bit more efficient computation-wise:

 t = p2-p1; u = p3-p1; v = p3-p2;
 w = cross(t,u);
 t2 = sum(t.^2); u2 = sum(u.^2); w2 = sum(w.^2);
 c = p1+(t2*sum(u.*v)*u-u2*sum(t.*v)*t)/(2*w2);
 r = 1/2*sqrt(t2*u2*sum(v.^2)/w2);

  In these computations the only thing you need to guard against is obtaining
a w2 value of zero (or very near it.) That would occur when triangle p1p2p3
has zero area and consequently an infinitely large radius for the circle.

  There is a corresponding two-dimensional version of the above, but I
suspect it would prove to be fairly similar to Yumnam Kirani Singh's formula
2.

Roger Stafford

Subject: Fit circle to 3 points

From: Peter Bone

Date: 11 Mar, 2008 16:39:03

Message: 7 of 9

> Just in case it is of interest to you, Peter, here is a
solution to your problem
> in three-dimensional space.
>
> Let p1 = (x1,y1,z1), p2 = (x2,y2,z2), and p3 =
(x3,y3,z3) be three points in
> 3D space, each represented by a three-element (row or
column) vector. Then
> a circle through these three points has its center at the
location given by
> vector c with radius r as given in the following:
>
> t = p2-p1; u = p3-p1; v = p3-p2;
> w = cross(t,u);
> c =
p1+(dot(t,t)*dot(u,v)*u-dot(u,u)*dot(t,v)*t)/(2*dot(w,w));
> r = 1/2*norm(t)*norm(u)*norm(v)/norm(w);
>
> However, the following is probably a bit more efficient
computation-wise:
>
> t = p2-p1; u = p3-p1; v = p3-p2;
> w = cross(t,u);
> t2 = sum(t.^2); u2 = sum(u.^2); w2 = sum(w.^2);
> c = p1+(t2*sum(u.*v)*u-u2*sum(t.*v)*t)/(2*w2);
> r = 1/2*sqrt(t2*u2*sum(v.^2)/w2);
>
> In these computations the only thing you need to guard
against is obtaining
> a w2 value of zero (or very near it.) That would occur
when triangle p1p2p3
> has zero area and consequently an infinitely large radius
for the circle.
>
> There is a corresponding two-dimensional version of the
above, but I
> suspect it would prove to be fairly similar to Yumnam
Kirani Singh's formula
> 2.
>
> Roger Stafford

Thanks Roger. How about the problem of fitting a hypersphere
to the minimum number of points that produce a unique
solution. I guess this would be n+1 points where n is the
dimension of the hypersphere. I don't need to do this. I'm
just interested if there's an analytical solution.

Pete

Subject: Fit circle to 3 points

From: Roger Stafford

Date: 11 Mar, 2008 23:37:01

Message: 8 of 9

"Peter Bone" <peterbone@hotmail.com> wrote in message <fr6cj7$555
$1@fred.mathworks.com>...
> Thanks Roger. How about the problem of fitting a hypersphere
> to the minimum number of points that produce a unique
> solution. I guess this would be n+1 points where n is the
> dimension of the hypersphere. I don't need to do this. I'm
> just interested if there's an analytical solution.
>
> Pete
---------
  Remarkably enough, the center of such a hypersphere can be found from
the solution to a certain set of n simultaneous linear equations in its n
coordinates. Let the given n+1 points have coordinates

 P0 = (p01,p02,p03,...,p0n)
 P1 = (p11,p12,p13,...,p1n)
 P2 = (p21,p22,p23,...,p2n)
 ...
 Pn = (pn1,pn2,pn3,...,pnn)

Call the coordinates of the hypersphere center C = (x1,x2,...,xn) and its
radius r.

  The following n+1 equations must be satisfied

 (p01-x1)^2+(p02-x2)^2+...+(p0n-xn)^2 = r^2
 (p11-x1)^2+(p12-x2)^2+...+(p1n-xn)^2 = r^2
 (p21-x1)^2+(p22-x2)^2+...+(p2n-xn)^2 = r^2
 ...
 (pn1-x1)^2+(pn2-x2)^2+...+(pnn-xn)^2 = r^2

Subtracting the first of these from each of the others and simplifying will give
us

 (p11-p01)*x1+(p12-p02)*x2+...+(p1n-p0n)*xn =
    ((p11^2+p12^2+...+p1n^2)-(p01^2+p02^2+...+p0n^2))/2

 (p21-p01)*x1+(p22-p02)*x2+...+(p2n-p0n)*xn =
    ((p21^2+p22^2+...+p2n^2)-(p01^2+p02^2+...+p0n^2))/2

 ...

 (pn1-p01)*x1+(pn2-p02)*x2+...+(pnn-p0n)*xn =
    ((pn1^2+pn2^2+...+pnn^2)-(p01^2+p02^2+...+p0n^2))/2

These are the n equations in n unknowns.

  Matlab can easily solve them to get the center coordinates and then r is
determined from the first equation above involving p01,p02,...p0n. As for an
analytic solution, it is of course possible using Cramer's rule, but the
resulting expressions would be rather cumbersome.

  One observation here. Finding a hypersphere and finding a circle are not
the same things for more than two dimensions. It still requires only three
points in n-dmensional space to determine a circle there. Note, by the way,
that in that previous article in the 3D case, I failed to identify the orientation
of the circle whose center and radius had been determined. The orthogonal
direction to the circle's plane would be along the w = cross(t,u) direction.

Roger Stafford

Subject: Fit circle to 3 points

From: rue mohr

Date: 17 Mar, 2008 01:08:03

Message: 9 of 9

here is some C code I just wrote for the same thing I just
went to a lot of work to register just to post this, so I
hope its of some use. it does need adjustments :)

rue_mohr#robotics@irc.linux.org

#include "2d.h"
#include <stdio.h>

typedef struct linef_s {
  float m;
  float b;
} linef_t;

typedef struct circlef_s {
  float x;
  float y;
  float r;
} circlef_t;


void printPoint( char * name, point2d_t * P) {
   printf("%s: %0.4f,%0.4f\n", name, P->x, P->y);
}

void printLine( char * name, linef_t * L) {
   printf("%s: y = %0.4fx + %0.4f\n", name, L->m, L->b);
}

inline void midpoint( point2d_t * A, point2d_t * B,
point2d_t * result ) {
   result->x = (A->x+B->x)/2.0;
   result->y = (A->y+B->y)/2.0;
   return;
}

inline void slope( point2d_t * A, point2d_t * B, linef_t *
result ) {
  result->m = (A->y-B->y)/(A->x-B->x);
  return;
}

inline void perpendicular( linef_t * result ) {
  result->m = (-1.0/result->m);
  return;
}

inline void pointSlopeToOff( point2d_t * A, linef_t * result) {
  result->b = A->y - (result->m * A->x);
  return;
}

void intersection ( linef_t * A, linef_t * B, point2d_t *
result) {
  result->x = (A->b-B->b)/(B->m-A->m);
  result->y = ((A->m*B->b) - (A->b*B->m))/(A->m-B->m);
  return;
}


// take 3 points and find the circle that goes thu em.
void circleThru( point2d_t * A, point2d_t * B, point2d_t *
C, point2d_t * result) {
  point2d_t Mab, Mbc;
  linef_t AB, BC;
  
  // get middle of segments
  midpoint( A, B, &Mab );
  midpoint( B, C, &Mbc );
  
  // find slopes
  slope( A, B, &AB);
  slope( B, C, &BC);
  
  // invert slopes
  perpendicular( &AB);
  perpendicular( &BC);
  
  // solve lines for points
  pointSlopeToOff( &Mab, &AB);
  pointSlopeToOff( &Mbc, &BC);
  
  // find line intersections
  intersection( &AB, &BC, result);
  
  return;
}



int main(void) {
  point2d_t A = {2,3};
  point2d_t B = {10,2};
  point2d_t C = {20,30};
  
  point2d_t answer;
  
  circleThru( &A, &B, &C, &answer);
    
  printf("Circle middle is: %f,%f\n", answer.x, answer.y);

  return 0;
}

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us