Given 4 points (non coplanar), how does one find the sphere, that is, center and radius, exactly fitting those points?
Call the four points A, B, C and D. Any three of them must be non-collinear (otherwise all three could not lie on the surface of a sphere) and all four must not be coplaner (otherwise either they cannot all lie on a sphere or they define an infinity of them).
A, B and C define a circle. The perpendicular bisectors of AB, BC and CA meet in a point (P, say) which is the centre of this circle. This circle must lie on the surface of the desired sphere.
Consider the normal to the plane ABC passing through P. All points on this normal are equidistant from A, B and C and its circle (in fact it is a diameter of the desired sphere). Take the plane containing this normal and D (if D lies on the normal any plane containing the normal will do); this plane is at right angles to the ABC one.
Let E be the point (there are normally two of them) on the circumference of the ABC circle which lies in this plane. We need a point Q on the normal such that EQ = DQ. But the intersection of the perpendicular bisector of ED and the normal is such a point (and it exists since D is not in the plane ABC, and so ED is not at right angles to the normal).
Algorithm:
Is the sphere well defined?
(1) Check that A and B are not coincident (=> failure).
(2) Find the line AB and check that C does not lie on it (=> failure).
(3) Find the plane ABC and check that D does not lie in it (=> failure).
Yes. Find its centre.
(1) Find the perpendicular bisectors of AB and AC.
(2) Find their point of intersection (P).
(3) Find the normal to the plane ABC passing through P (line N).
(4) Find the plane containing N and D; find the point E on the
ABC circle in this plane (if D lies on N, take E as A).
(4) Find the perpendicular bisector of ED (line L)
(5) Find the point of intersection of N and L (Q).
Q is the centre of the desired sphere
Pictures:
(1) In the plane ABC
A
P
B C
(2) At right-angles to ABC, in the plane containing N and D
E
D
line N
--------------------P-------------Q---------------------------
Numerically:
If ED << EP then Q will be very close to P (relative to the radius of the ABC circle) and subject to error. It's best to choose D so that the least of AD, BD and CD is larger than for any other choice.
Off the top of my head, I might try:
Given: p_1, p_2, p_3, p_4
Find: p_c (center of sphere determined by p_1, ..., p_4), dist(p_c, p_i) (radius)
p_c is the same distance from our four points, so dist(p_c,p_1) = dist(p_c,p_2) = dist(p_c,p_3) = dist(p_c,p_4)
Of course, we can square the whole thing to get rid of square roots: distsq(p_c,p_1) = distsq(p_c,p_2) = distsq(p_c,p_3) = distsq(p_c,p_4)
Plug in the variables into the distance formula, simplify, and the x^2_c, y^2_c, and z^2_c terms cancel out, leaving you with three linearly independent equations and three unknowns (x_c, y_c, z_c). Solve using your favorite method. :)
There is another useful method based on Least Sqyares Estimation of the sphere equation parameters.
The points (x,y,z) on a spherical surface with radius R and center (a,b,c) can be written as
(x-a)^2 + (y-b)^2 + (z-c)^2 = R^2
This equation can be rewritten into the following form:
2ax + 2by + 2cz + R^2 - a^2 - b^2 -c^2 = x^2 + y^2 + z^2
Approximate the left hand part by F(x,y,z) = p1.x + p2.x + p3.z + p4.1
For all datapoints, i.c. 4, determine the 4 parameters p1..p4 which minimise the average error |F(x,y,z) - x^2 - y^2 - z^2|^2.
In 'Numerical Recipes in C' can be found algorithms to solve these parameters.
The best fitting sphere will have
So, at last, will this solve you sphere estination problem, at least for the most situations I think ?.
I don't remember who the original poster of this problem was, but I hope the following information help if (s)he is reading:
"How to find the equation of a sphere given 4 surface points"
From analytic geometry, we know there is a unique sphere, say
c1(x^2+y^2+z^2) + c2(x) + c3(y) + c4(z) + c5 = 0
that passes through the four noncoplanar points
(x1,y1),(x2,y2),(x3,y3),(x4,y4)
and is given by the following determinant equation:
| ( x^2 + y^2 + z^2) x y z 1 |
| |
| (x1^2 + y1^2 + z1^2) x1 y1 z1 1 |
| |
| (x2^2 + y2^2 + z2^2) x2 y2 z2 1 | = 0
| |
| (x3^2 + y3^2 + z3^2) x3 y3 z3 1 |
| |
| (x4^2 + y4^2 + z4^2) x4 y4 z4 1 |
Example: Suppose, you have the four 3D points
(0,3,2), (1,-1,1), (2,1,0), and (5,1,3)
The equation of the sphere that passes through these points is
| (x^2 + y^2 + z^2) x y z 1 |
| |
| 13 0 3 2 1 |
| |
| 3 1 -1 1 1 | = 0
| |
| 5 2 1 0 1 |
| |
| 35 5 1 3 1 |
This reduces to
x^2 + y^2 + z^2 - 4x - 2y - 6z + 5 = 0
which in standard form is
(x-2)^2 + (y-1)^2 + (z-3)2 = 9
Q.E.D.
To answer the original poster's question, the center of the spere is given by (2,1,3).
Three points can never determine a unique sphere. If they are collinear and distinct, they do not lie on any circle (and hence not on any sphere either). If they are not collinear, they determine a unique circle, which is contained by an infinity of spheres: think of a ring with balls of varying sizes balanced on it.
Four points determine a unique sphere if, and only if, they are not coplanar. It is not necessary to add "and if no three are collinear," because in such a case the four will necessarily be coplanar. If they are coplanar, either there are no spheres through the 4 points, or an infinity of them (if the 4 points lie on some circle).
Let the 4 points be (x1, y1, z1) etc. They are coplanar if and only if the following determinant is zero:
The equation of the sphere in this case has been given in a previous post, but it is worth repeating; it is given by setting the following determinant to zero:| x1 y1 z1 1 | | x2 y2 z2 1 | | x3 y3 z3 1 | | x4 y4 z4 1 |
It is easy to see that this is in fact the equation of a sphere (when the points are not coplanar), and equally easy to see that the four points lie on the sphere thus determined. Hence it is the equation of the (unique) sphere that contains the points.| x^2 + y^2 + z^2 x y z 1 | | x1^2 + y1^2 + z1^2 x1 y1 z1 1 | | x2^2 + y2^2 + z2^2 x2 y2 z2 1 | = 0. | x3^2 + y3^2 + z3^2 x3 y3 z3 1 | | x4^2 + y4^2 + z4^2 x4 y4 z4 1 |