# How to determine the center and radius of a circle given some points in 3D?

I was wondering if anyone could give me a hand with this problem I have.

I have six points on a plane, and I am trying to determine if they form a circle or not.

I know that any three points in 2D can determine a circle as long as the three points don't all lie on the same line. The way I figured that problem was taking two points, finding the line that connects them, finding the midpoint of that segment and then finding a perpendicular to the segment at that point. I do that with two segments and the intersection of the perpendicular lines that I find intersect at the center of the circle.

When I follow the steps for the points I have in 3D, Mathematica "gets stuck" when I ask it to find the equation of the line with the slope and the two points I give it (because the points have 3 coordinates and not two like in 2D).

Any help is greatly appreciated!

Thank you,

• Could you please clarify the sense in which it is consistent first to state your points are "on a plane" and later to say they are "in 3D"? If you mean they are known to lie exactly on some affine plane in Euclidean 3-space, why don't you just adopt coordinates for that plane and work therein, reducing this to a 2D problem? – whuber Dec 12 '12 at 20:46
• Some code would be really nice. – Vitaliy Kaurov Dec 12 '12 at 20:48
• What does "I have six points on a plane" have to do with problem? A circle is always determined by three points. If you have six points, find two circles determined by disjoint subsets of three points and determine if the two circles are equal. If there are, then all six points lie on the same circle – m_goldberg Dec 12 '12 at 21:22
• Orthogonalize can be your friend here. First average the points (that is, compute their centroid) and subtract it to make all coordinates relative to the centroid. Applying Orthogonalize to the matrix of coordinates should produce just two nonzero orthogonal basis vectors. (If not, your points are not coplanar.) Rewrite each shifted point as a linear combination of those basis vectors: the coefficients are your new coordinates. Details are given at mathematica.stackexchange.com/questions/10333 (and its comments: J.M. provided some useful pointers). – whuber Dec 12 '12 at 21:22
• Are you dealing with points that lie exactly (or to high precision) on a circle? Or possibly only approximately so? For the former case you have solid responses. For the latter you could use those as a start and try to find a best fit solution from there. – Daniel Lichtblau Dec 13 '12 at 15:02

Here are three points in space.

SeedRandom[2];
{p1, p2, p3} = {{x1, y1, z1}, {x2, y2, z2}, {x3, y3, z3}} =
RandomReal[{-3, 3}, {3, 3}];


If the center is $p=(x,y,z)$ and the radius is $r$, then the distance from $p$ to each $p_i$ must be exactly $r$. Thus, for each $i=1,2,3$, we have

$$(x - x_i)^2 + (y - y_i)^2 + (z - z_i)^2 = r^2.$$

Furthermore, the point $p$ must lie on the plane determined by $p_1$, $p_2$, and $p_3$. If $\vec{n}$ is a vector normal to that plane, then an equation for the plane is

$$\vec{n}\cdot(p-p_1)=0.$$

That gives us four equations in the four unknowns $x$, $y$, $z$, and $r$. Using the fact that $\vec{v}_1 \times \vec{v}_2$ is always perpendicular to $\vec{v}_1$ and $\vec{v}_2$, we can set up this system of equations as follows:

v1 = p2 - p1;
v2 = p3 - p1;
{v1, v2} = Orthogonalize[{v1, v2}];
n = Cross[v1, v2];
eqs = {
(x - x1)^2 + (y - y1)^2 + (z - z1)^2 == r^2,
(x - x2)^2 + (y - y2)^2 + (z - z2)^2 == r^2,
(x - x3)^2 + (y - y3)^2 + (z - z3)^2 == r^2,
n.({x, y, z} - p1) == 0
};


Strictly speaking, the Orthogonalize step is not necessary, but it will be convenient a bit later.

Here is the solution of the system. We want the positive $r$ value, of course.

solution = NSolve[eqs, {x, y, z, r}]

(* Out:
{{x -> -0.333316, y -> -0.850768, z -> 0.653921, r -> 2.38624},
{x -> -0.333316, y -> -0.850768, z -> 0.653921, r -> -2.38624}}
*)


We now have enough information to parametrize the circle and plot it together with the points, so we can see how this all looks. Note that in the Mathematica code, I'm using {x0,y0,z0} to denote the center and r0 to denote the radius. Also, this depends on the fact that v1 and v2 are perpendicular unit vectors, which is one reason we used the Orthogonalize command.

{x0, y0, z0} = {x, y, z} /. First[solution];
r0 = r /. First[solution];
Show[{
ParametricPlot3D[{x0, y0, z0} + r0*v1*Cos[t] + r0*v2*Sin[t], {t, 0, 2 Pi}],
Graphics3D[{PointSize[Large], Point[{p1, p2, p3}]}]
}]


The tight integration of symbolics and graphics is one of the great strengths of Mathematica and this example shows why it's so nice - it's super simple to see that we've done things correctly!

Now, in order to check if your six points all lie on the same circle, we just need a way to check if any particular point lies on that circle. Thus, suppose we have another point that we suspect must lie on this circle, say:

{x4, y4, z4} = {x0, y0, z0} + r0*v1*Cos[t] + r0*v2*Sin[t] /.
t -> RandomReal[{0, 2 Pi}];


It's distance to {x0,y0,z0} should be r0, so the following should be zero:

Chop[(x4 - x0)^2 + (y4 - y0)^2 + (z4 - z0)^2 - r0^2]

(* Out: 0 *)


and it should also lie on the plane:

Chop[n.({x4, y4, z4} - p1)]

(* Out: 0 *)


## Mapping to a plane

I don't really see a major advantage to working in a plane, but it can certainly be done. Consider the following $3\times3$ matrix $M$ defined using column vectors $\vec{v}_1$, $\vec{v}_2$, and $\vec{v}_3$:

$$M = \left( \begin{array}{ccc} \vec{v}_1 & \vec{v}_2 & \vec{v}_3 \end{array} \right).$$

If ${\mathbb R}^3$ has the standard unit basis vectors $\vec{e}_1$, $\vec{e}_2$, and $\vec{e}_3$, then direct computation shows that

$$M \, \vec{e}_i = \vec{v}_i.$$

As a result, $M^{-1}$ does the opposite:

$$M^{-1} \, \vec{v}_i = \vec{e}_i.$$

In particular, if we take $\vec{v}_1$ and $\vec{v}_2$ to be the vectors defined above and then set the third vector to be our normal vector, $\vec{v}_3 = \vec{n}$, then $M^{-1}$ should map our three original points to some horizontal plane.

This is an absolutely standard trick in linear algebra that helps us understand the basics behind changes of bases. It's also a second point where it is nice that v1 and v2 are perpendicular, unit vectors, for this implies that $M$ and $M^{-1}$ are orthogonal matrices and, therefore, preserve the shape of the original configuration. Without this assumption, $M$ could be a shear transformation so the image of the points under $M^{-1}$ wouldn't necessarily lie on a circle, though they would lie on a plane.

Let's try it - I'm using rr to denote the original matrix $M$ (I'm thinking of it as being rotation, which it almost might be) and I'm using RR to denote the inverse.

rr = Transpose[{v1, v2, n}];
RR = Inverse[rr];
RR.# & /@ {p1, p2, p3}

(* Out: {
{-2.46374, -1.07249, -0.282502},
{0.770168, -1.07249, -0.282502},
{0.729536, 2.47387, -0.282502}
}
*)


After application of the matrix RR, all the points have the same $z$ value. We could simply map Most onto this result if we want 2D points.

• If you check first three points and next three other points, they could lie on two separate circles. So you have provided necessary but insufficient conditions. – whuber Dec 12 '12 at 21:07
• Fair enough: My comment was intended as a (perhaps too subtle) suggestion that your answer would be better for making that point explicit and showing just how the check would be carried out. You will also find that your solution fails in special circumstances, such as where points are coincident or collinear or where coordinates are subject to floating point imprecision. – whuber Dec 12 '12 at 21:18
• @whuber Is there any way you could help me get my original problem back to a 2D problem? – RedPotatoe Dec 12 '12 at 21:22
• @RedPotatoe I made some pretty major changes. Does that help? – Mark McClure Dec 13 '12 at 5:15
• @whuber Is that any better? – Mark McClure Dec 13 '12 at 5:15

This question of determining whether points form a circle or not can be resolved by solving two sub-problems:

1. Express the points with just two coordinates rather than $n=3$.

2. Test whether a collection of points in a plane is concyclic.

We might as well solve the first sub-problem for dimensions $n\ge 3$, because it all works the same: we find an orthonormal basis of $n$-space where the first two basis elements (appropriately translated) span the plane containing the points. Expressing the points in that basis gives the desired 2D coordinates.

An elegant way to identify concyclic planar points is to view them as complex numbers and check them in groups of four each to see whether they have a real cross ratio: if so, each group lies on a common generalized circle (that is, they are either collinear or concyclic). We may pick any three distinct points and, adjoining each of the remaining points in turn, check each of these groups of four. If they are all found to have real cross-ratios, then they all must lie on the (unique) generalized circle determined by those three distinct points.

This plan is executed by concyclicQ, which includes the checks for degenerate situations (such as all points collinear). Some care is needed in the use of Orthogonalize: just a little bit of noise in the data orthogonal to the plane can cause it to conclude the points are not coplanar. To force the points to be considered coplanar, just remove the check If[Length[f] > 2.... (Notice the uses of Chop to reduce the effects of floating point imprecision.)

crossRatio[z1_, z2_, z3_, z4_] := (z1 - z3) (z2 - z4) / ((z2 - z3) (z1 - z4));
concyclicQ[pts_List] /; Length[pts] >= 3 := Module[{points, k, f, g, q, p, test, tol=10^-6},
(* Step 1: project into the Complex plane *)
points = Union[pts];                             (* Eliminate any duplicates *)
q = # - Mean[points] & /@ points;                (* Center the points *)
f = Select[Chop[Orthogonalize[q, Tolerance -> tol]], Norm[#] > 0 &]; (* Find a basis *)
If[Length[f] == 0, Return[True]];                (* All are coincident *)
If[Length[f] == 1, Return[Length[points] <= 2]]; (* Collinear *)
If[Length[f] > 2, Return[False]];                (* Not coplanar *)
f = Select[Chop[Orthogonalize[f~Join~IdentityMatrix[Dimensions[points][[2]]]]],
Norm[#] > 0 &];                                (* Find an adapted basis *)
g = ( Inverse[f] )[[All, 1 ;; 2]];               (* Change-of-basis matrix *)
p = Complex @@@ (q . g);                         (* Change the basis *)

(* Step 2: test for concyclicality *)
test = Append[p[[1 ;; 3]], #] & /@ p[[4 ;;]];    (* Test last n-3 against the first 3 *)
Min[Boole[# \[Element] Reals] & /@ Chop[crossRatio @@@ test, tol]] == 1
];
concyclicQ[points_List] := True (* Three points or fewer *)


This code works in greater than three dimensions, too.

Examples:

With[{n = 600, e = Orthogonalize[RandomReal[NormalDistribution[0, 1], {3, 3}]]},
points = Join[#, ConstantArray[1,  3 - 2]] & /@ (Through[{Cos, Sin}[
RandomReal[{0, 2 \[Pi]}, n]]]\[Transpose]) . e;
]; (* Many random concyclic points *)
concyclicQ[points]


True

q = Append[points, {0, 0, 0}]; (* Throw in a "bad" point *)
concyclicQ[q]


False

concyclicQ[RandomReal[{0, 1}, {4, 3}]] (* Four points in general position *)


False

concyclicQ[{RandomReal[{0, 1}, 4]}\[Transpose] . {{1, 1, 1}}] (* Collinear *)


False

concyclicQ[ConstantArray[ {0, 0, 0}, 4]] (* Coincident *)


True

### Edit

In comments to the question, Daniel Lichtblau pointed out that a solution that handles a little "noise" or positional error would be advantageous. Here is one using Singular Value Decomposition to project the points to a plane, find a least-squares fit to a circle there, and test the accuracy of the fit relative to the fitted radius. This one makes no effort to look for exceptional cases, because they shouldn't exist: when the points are collinear, the fit ought to find an extremely large radius and conclude that--up to the specified tolerance--the points are indeed concyclic!

Clear[concyclicQ];
concyclicQ[points_List, tol_: 10^-8] :=
Module[{u, w, v, s2, solution, fit, x, y, r, objective},
fit[p_] := FindMinimum[objective[{x, y}, r, p],
{{x, 0}, {y, 0}, {r,  Mean[Max[#] - Min[#] & /@ Transpose[p]]}}];
objective[center_, radius_, q_] := Total[(Norm[#] - radius)^2 & /@ (# - center & /@ q)];
{u, w, v} = SingularValueDecomposition[q = # - Mean[points] & /@ points];
{s2, solution} = fit[u[[All, 1 ;; 2]]. w[[1 ;; 2, 1 ;; 2]]];
(Sqrt[s2]/r /. solution) <= tol
]


The first two lines define a fitting function fit which estimates reasonable starting values of the center and radius of the circle and uses objective to achieve the least-squares solution. The next line computes the SVD of the points, which is then limited to the eigenspaces of the two largest eigenvalues in the call to fit. The root mean square returned by fit is divided by the fitted radius r. This would be an acceptable value to return, but to make this version compatible with the previous one, I then compare this relative residual error to a user-supplied tolerance: when the error is sufficiently small, it is fair to conclude the points are concyclic.

Examples:

concyclicQ[points]


True

concyclicQ[RandomReal[NormalDistribution[0, 1], {6, 3}]]


False

• It might be useful to note that SingularValueDecomposition[] allows the output of just the first few singular values/vectors (e.g. SingularValueDecomposition[mat, 2] to take the system corresponding to the two largest singular values). – J. M.'s technical difficulties Apr 13 '13 at 16:21

Here is an approach to finding a best fit circle given a set of points in R^3. I'll start with points that are a perturbation of cocircular points centerend at the origin in the x-y plane. I then rotate and translate them to get our example set.

SeedRandom[11111];
pts = RandomReal[{0, 2*Pi}, 20];
circpts =
Map[4*{Sin[#], Cos[#], 0.} &, pts] +
RandomReal[{-.1, .1}, Dimensions[pts]];

tf = TranslationTransform[{4, -2, 3}];
rf = RotationTransform[3*Pi/5, {-3, 1, 2}];

circpts3D = tf[rf[circpts]];


Here they are.

ptplot = Graphics3D[{PointSize[Large], Red, Point /@ circpts3D},
BoxRatios -> {1, 1, 1}]


First I'll find the best fit to a plane containing these points. Tha's a very simple least squares computation. (All "best" fits herein will be in terms of a least squares measure. Of something.)

planedir = LeastSquares[circpts3D, ConstantArray[1, Length[circpts3D]]];

(* {0.104003, -0.325088, -0.0227852} *)


We now have to find a center and radius. This involves a certain amount of magic, below.

sphermat = ArrayPad[2*circpts3D, {{0, 0}, {0, 1}}, 1];
spherrhs = Map[#.# &, circpts3D];
spherparams = LeastSquares[sphermat, spherrhs];
cen = Most[spherparams]

(* {3.96641, -1.90807, 3.01841}
4.00409 *)


Notice that these are indeed quite close to the parameters of the unperturbed circle.

Here we'll see points, plance, and sphere all playing nicely with one another.

planeplot =
ContourPlot3D[
planedir.{x, y, z} == 1, {x, -4, 10}, {y, -4, 10}, {z, -4, 10},
Mesh -> False];

spherplot =
Graphics3D[{Opacity[.5], Sphere[cen, rad]}, BoxRatios -> {1, 1, 1}];

Show[planeplot, ptplot, spherplot,
PlotRange -> {{-1, 10}, {-7, 4}, {-2, 9}}, BoxRatios -> {1, 1, 1},
AspectRatio -> 1]


The "magic" I alluded to is that computing the center and radius as I did gives no guarantee that the center will not be elsewhere on the line through the true center and orthogonal to the plane. I believe we get the "correct" solution from LeastSquares using a pseudoinverse type of method, because those give minimal norm solutions. Were this not the case we'd need to also use our plane direction vector to put, in the least squares setup, an approximate orthogonality condition between that and the vectors formed by subtracting the center coordinates from each point.

I think this can all be wrapped into a single least squares computation but I already made that claim once and got lost trying to put the pieces into place. (Might not be all that difficult, some days I just get lost easily.)

For fitting a circle orthogonally in three dimensions, I use a procedure adapted from this paper. I'll show it here in steps; bundling all this into a routine is left up to you.

BlockRandom[SeedRandom[42, Method -> "MersenneTwister"];
rnd = AffineTransform[{Orthogonalize[RandomVariate[NormalDistribution[], {3, 3}]],
RandomReal[{-1, 1}, 3]}];
pts = rnd[Append[(2 + RandomVariate[NormalDistribution[]]/10) Normalize[#],0.]] & /@
RandomVariate[NormalDistribution[], {20, 2}]];

Graphics3D[{AbsolutePointSize[4], Point[pts]}]


First, find parameters for the Hessian normal form of the best fit plane:

cent = Mean[pts]; (* data mean, point on the best-fit plane *)
(* normal to best-fit plane *)
nrm = Flatten[Last[SingularValueDecomposition[TranslationTransform[-cent] /@ pts, -1]]]


Now, make a few transformations so that the data lies in the $x$-$y$ plane and is centered at the origin:

toPlane = Composition[RotationTransform[{nrm, {0, 0, 1}}], TranslationTransform[-cent]];
newpts = Composition[Most, toPlane] /@ pts;


We make an initial circle fit on the transformed points:

{a, b1, b2, c} =
Flatten[Last[SingularValueDecomposition[Flatten[{Norm[#]^2, #, 1}] & /@ newpts, -1]]];
{cen, rad} = {-{b1, b2}/(2 a), Sqrt[(b1^2 + b2^2)/(4 a^2) - c/a]}


Check:

Graphics[{Circle[cen, rad], {AbsolutePointSize[4], Point[newpts]}}]


Now, display the initial 3D circle:

c3 = cent + RotationTransform[{{0, 0, 1}, nrm}][Append[cen, 0]];
Show[Graphics3D[{AbsolutePointSize[4], Point[pts], Arrow[{cent, cent + nrm}]}],
ParametricPlot3D[Evaluate @
c3 + RotationTransform[{{0, 0, 1}, nrm}][Append[rad {Cos[u], Sin[u]}, 0]],
{u, -π, π}]]


This is not yet a genuine orthogonal fit. We need to do one more polishing for that purpose:

{c3, nrm, {rad}} = Partition[FindArgMin[
Total[With[{pd = ({x1, x2, x3} - #).{n1, n2, n3}},
Norm[{pd, Norm[({x1, x2, x3} - #) - pd {n1, n2, n3}] - r}]] & /@ pts],
Append[Transpose[{{x1, x2, x3}, c3}] ~Join~ Transpose[{{n1, n2, n3}, nrm}], {r, rad}]
// Evaluate], 3, 3, 1, {}]
{{-0.39890712440524145, -0.5627585871305216, -0.011641064393032342},
{-0.7373064372198214, 0.5460388072274082, 0.39654401814726603},
{2.0349303367150706}}


We now show the final orthogonally-fitted circle:

circle3D[cent_?VectorQ, r_?NumericQ, no_?VectorQ] := BSplineCurve[
Map[AffineTransform[{RotationMatrix[{{0, 0, 1}, no}], cent}],
r {{1, 0, 0}, {1, 1, 0}, {-1, 1, 0}, {-1, 0, 0}, {-1, -1, 0}, {1, -1, 0}, {1, 0, 0}}],
SplineDegree -> 2, SplineKnots -> {0, 0, 0, 1/4, 1/2, 1/2, 3/4, 1, 1, 1},
SplineWeights -> {1, 1/2, 1/2, 1, 1/2, 1/2, 1}]

Graphics3D[{{AbsolutePointSize[4], Point[pts]}, Arrow[{c3, c3 + nrm}],


Here is another way to test whether points in space lie on a circle. It is composed of two parts. First, points that lie on a sphere satisfy an equation of the form $$A(x^2+y^2+z^2) + Bx + Cy + Dz + E = 0 \,,$$ for some unknowns $A$, $B$, $C$, $D$, $E$. Second, points that lie in a plane satisfy an equation of the form $$Fx + Gy + Hz + K = 0 \,,$$ for some unknowns $F$, $G$, $H$, $K$. For a given system of points, these can viewed as two systems of linear equations in the unknowns, one for a sphere and one for a plane (in nondegenerate configurations), the intersection of which would be a circle.

On the one hand, the rank of the first will always be at least the rank of the second, since a solution of the second is with $A=0$ also a solution of the first. So if the rank of the first system is at most 3, then the points satisfy equations of both kinds. On the other hand, there are degeneracies that one might wish to exclude, such as when the points lie on a line or are coincident. To exclude these, the ranks need to be at least 3. So the test for a circle is that both ranks should be exactly 3.

Two matrices represent the systems. MatrixRank has a Tolerance option, which can be passed through circle3DQ, and controls how close to a circle the points have to be. (The default tolerance is, I think, the same as Chop, 10^-10.)

sphereMatrix[pts_?MatrixQ] := ArrayFlatten[{{List /@ Total /@ (pts^2), pts, 1}}];
planeMatrix[pts_?MatrixQ] := ArrayFlatten[{{pts, 1}}];
circle3DQ[pts_?(MatrixQ[#, NumericQ] &), opts___] :=
MatrixRank[planeMatrix[pts], opts] == 3 && MatrixRank[sphereMatrix[pts], opts] == 3;


Some test points, and a function tweak that nudges the points off the circle.

radii = Orthogonalize[RandomReal[{-1, 1}, {2, 3}]~Join~IdentityMatrix[3]][[1 ;; 2]];
center = RandomReal[{-1, 1}, 3];
pts = center + {Cos[t], Sin[t]}.radii /. List /@ Thread[t -> RandomReal[{0, 2π}, 6]];
tweak[tolerance_] := # + RandomReal[{-tolerance, tolerance}, 3] &;


Tests:

circle3DQ[pts]
(* True *)

circle3DQ[tweak[10^-6] /@ pts]
(* False *)

circle3DQ[tweak[10^-6] /@ pts, Tolerance -> 10^-6]
(* True *)


Tweaking the points up to about $\pm 10^{-14}$ doesn't change the result of circle3DQ. The following counts the number of tweaked sets of points identified as non-circles under the default Tolerance.

Count[Table[circle3DQ[tweak[10^-14] /@ pts], {100}], False]
(* 0 *)

Count[Table[circle3DQ[tweak[5 * 10^-14] /@ pts], {100}], False]
(* 52 *)


One can adjust the tolerance so that if the points are tweaked by 10^-13, circle3DQ still tests True.

Count[Table[circle3DQ[tweak[5 * 10^-14] /@ pts, Tolerance -> 10^-9], {100}], False]
(* 0 *)


For Tolerance -> 0, the coordinates need to be exact.

circle3DQ[pts, Tolerance -> 0]
(* False *)


With exact coordinates:

radii1 = Orthogonalize[RandomInteger[{-100, 100}, {2, 3}]~Join~IdentityMatrix[3]][[1 ;; 2]];
center1 = RandomInteger[{-100, 100}, 3];
pts1 = center1 + {Cos[t], Sin[t]}.radii1 /.
List /@ Thread[t -> Round[RandomReal[{0, 2π}, 6], 10^-2]];

circle3DQ[pts1, Tolerance -> 0]
(* True *)


Some exceptional test configuration:

coincident = ConstantArray[{1, 2, 3}, 6];
colinear = Array[# {1, 2, 3} &, 6];
coplanar = RandomReal[{-1, 1}, {6, 3}];
coplanar[[All, 3]] = 1;
spherical = Normalize /@ RandomReal[{-1, 1}, {6, 3}];
random = RandomReal[{-1, 1}, {6, 3}];


The table entries can be computed with

{circle3DQ@#, MatrixRank@planeMatrix@#, MatrixRank@sphereMatrix@#} & /@
{coincident, colinear, coplanar, spherical, random}


Anyone interested in an analytic solution from Yumnam Kirani Singh may check this link.

He explains some formulae for finding the centre of a circle or sphere, or a hypersphere in general, which will pass through some given number of points.

• Analytic in the sense that it is the determinantal form of a solution to a system of linear equations. Okay for low dimensions, I guess. – Daniel Lichtblau Dec 15 '12 at 23:10
• Mainly yes. In the 4th formula says"If the number of points in n-dimensions is less than or equal to n, then we can find infinite number of n-dimensional hypersphere" – s.s.o Dec 15 '12 at 23:54