Sphere equation given 4 points

Question

Find the equation of the Sphere give the 4 points (3,2,1), (1,-2,-3), (2,1,3) and (-1,1,2).

The *failed* solution I tried is kinda straigh forward:

We need to find the center of the sphere.

Having the points:

$$ p_{1}(3,2,1),\, p_{2}(1,-2,-3),\, p_{3}(2,1,3),\, p_{4}(-1,1,2) $$

2 Triangles can be created using these points, let's call $A$ our triangle formed by the points $p_{1},\,p_{2}\, and\,\, p_{3}$; And $B$ our triangle formed by the points $p_{1},\, p_{3}\, and \,\,p_{4}$.

Calculate the centroids of each triangle: $$ CA = (2,1/3,1/3)\\ CB = (4/3,4/3,2) $$ And also, a normal vector for each triangle: $$ \overrightarrow{NA} = \overrightarrow{p_{1}p_{2}} \times \overrightarrow{p_{1}p_{3}}\\ \overrightarrow{NB} = \overrightarrow{p_{1}p_{3}} \times \overrightarrow{p_{1}p_{4}} $$

$$ \overrightarrow{p_{1}p_{2}} = <-2,-4,-4>\\ \overrightarrow{p_{1}p_{3}} = <-1,-1,2>\\ \overrightarrow{p_{1}p_{4}} = <-4,-1,1>\\ \:\\ \overrightarrow{NA} = <-12, 8, -2>\\ \overrightarrow{NB} = <1, -7, -3>\\ $$

With the centroids and normals of triangles $A$ and $B$, we can build two parametric equations for a line, the first one intersects the centroid of $A$ and the other one the centroid $B$. $$ Line \enspace A\\ x = 2-12t \quad y = 1/3+8t \quad z = 1/3-2t\\ \:\\ Line \enspace B\\ x = 4/3 + s \quad y = 4/3 - 7s \quad z = 2 - 3s $$

The point where this lines intersect should be the center of the sphere, unfortunately this system of equations is not linearly dependent, that means that they do not intersect each other. What could be the problem here?

Answer 1

I would cite the beautiful method from W.H.Beyer to find the center and radius of the sphere $(x-a)^2+(y-b^2)+(y-c)^2=R^2$

Answer 2

Another method is to start with the equation of the sphere:

$$(x-u)^2+(y-v)^2+(z-w)^2=r^2$$

$(u,v,w)$ are the coordinates of the center of the sphere and $r$ is the radius. Plug into the given points $p_1, p_2, p_3, p_4$ for $x,y,z$ in the equation and you get four equations with variables $u,v,w,r$. But these equations contain quadratic terms. Subtract one equation from each of the other three. You get three linear equations and the variables $u, v, w$. Such a system can be solved in the usual way to find the center $(u,v,w)$, e.g. Gaussian elimination. Now plug the calculated $u,v,w$ into one of the initial equations to find the radius $r$.

Note:

The three linear equations found are the equation of the planes described by @Jyrki Lahtonen

The problem can be solved with this Maxima program:

/* define the four points p[1],...,p[4] */
p[1]:[3,2,1]; p[2]:[1,-2,-3]; p[3]:[2,1,3]; p[4]:[-1,1,2];

/* ceq is the equation of the circle
(u,v,w) is the center, r is the radius */
ceq:(x-u)^2+(y-v)^2+(z-w)^2=r^2;

/* plug in the points in the circle equation 
to get 4 equation eq1[1],...,eq1[4] */
for i:1 thru 4 do 
 eq1[i]:ev(ceq,map("=",[x,y,z],p[i]));

/* display this 4 equations eq1*/
listarray(eq1);

/* subtract  the fourth equation for each of the first 
three equation to get three linear 
equations eq2[1],eq2[2],eq2[3] in u,v,w */

for i:1 thru 3 do 
 eq2[i]:ev(eq1[i]-eq1[4],expand);

/* display these three equations 
eq2[1],eq2[2],eq2[3] */
listarray(eq2);

/* solve this system oflinear euations to get u,v,w */
ss:solve(listarray(eq2),[u,v,w]);

/* plugin the solutions in the fourth 
equation eq1[4] and calculate r^2 */
solve(eq1[4],r^2),ss[1];

Actually the problem can be solved using only one Maxima solve command:

(%i1) display2d:false$
(%i2) solve(makelist(ev((z-w)^2+(y-v)^2+(x-u)^2=r^2,
map("=",[x,y,z],p)),p,[[3,2,1],[1,-2,-3],[2,1,3],[-1,1,2]]),
[u,v,w,r]);
(%o2) [[u = 24/19,v = -16/19,w = 4/19,r = -3*sqrt(470)/19],
       [u = 24/19,v = -16/19,w = 4/19,r = 3*sqrt(470)/19]]

Answer 3

The centroid of a triangle is usually not at the same distance from the vertices. In the plane the center of the outer circle (= the intersection of the normals bisecting the sides) has this property.

But, because this is in 3D you might as well find the point of intersection of the planes the have $\vec{p_1p_2}$, $\vec{p_2p_3}$ and $\vec{p_3p_4}$ as their normals, and pass through the midpoints of those line segments. The points on the first plane are at the same distance from $p_1$ and $p_2$ et cetera. Thus the intersection of those three planes is what you want.

Answer 4

Hint

Using a totally different approach, you can also look at your problem using what JJacquelin (participant to MSE) proposed in his book

http://fr.scribd.com/doc/14819165/Regressions-coniques-quadriques-circulaire-spherique

Pages 17 and 18 give the full approach for a spherical regression. It is quite simple and reduces to a linear system of four equations for four unknowns from which are deduced the coordinates of the center and the radius of the sphere.

Since I do not want to spoil him, I shall not give more. I applied his method to your case and it works splendid.

Answer 5

For the ones calculating by hand, the error-prone numbers are:

the center of the sphere is at $p=(24/19,-16/19,4/19)$

the squared radius is $r^2=4230/361$.

Here is my small GNU Maxima script

---

display2d : false;

/*
* purpose:
*   given points x_k,
*   calculate circumsphere center p and squared radius r^2
**/

my_A(d,k,x) := block(

    [res,i,j],

    res : zeromatrix(d,d),
    for j : 0 thru k-1 do block(
        for i : 1 thru d do block(
            res[i,j+1] : x[j+1][i]-x[k+1][i]
        )
    ),
    for j : k+1 thru d do block(
        for i : 1 thru d do block(
            res[i,j] : x[j+1][i]-x[k+1][i]
        )
    ),

    return(res)

);

my_norm2(d,x) := apply("+",makelist(x[k]^2,k,1,d));
my_dist2(d,x,y) := apply("+",makelist((x[k]-y[k])^2,k,1,d));

x : [[3,2,1],[1,-2,-3],[2,1,3],[-1,1,2]];

d : length(x)-1;
A : makelist(my_A(d,k,x),k,0,d);
AinvT : makelist(transpose(invert(A[k+1])),k,0,d);
xnorm2 : makelist(my_norm2(d,x[k+1]),k,0,d);

z1 : transpose(makelist(1,k,1,d));
p : -1/2*apply("+",makelist(xnorm2[k+1]*(AinvT[k+1].z1),k,0,d));
r2 : makelist(my_dist2(d,transpose(p)[1],x[k+1]),k,0,d);

Answer 6

The center of the circle will lie on the plane going through the middle of each pair of points which is perpendicular to the line connecting them.

This means that for each point pair $P_i$, $P_j$ you can create a plane on which the center will lie:

$$(P_i-P_j)\times(\frac{P_i+P_j}{2}-C) =0 $$

Take 3 pairs (using all 4 non-planar points) and solve for C.