Center and Radius of a Sphere from Four Points
by Stephen R. Schmitt
Given 4 points, how does one find the sphere, that is, center and radius, exactly fitting those points? Four points determine a unique sphere if, and only if, they are not on the same plane. If they are on the same plane, either there are no spheres through the 4 points, or an infinity of them if the 4 points are on a circle. From analytic geometry, we know that there is a unique sphere that passes through the four noncoplanar points:
(x_{1}, y_{1}, z_{1}), (x_{2}, y_{2}, z_{2}), (x_{3}, y_{3}, z_{3}), (x_{4}, y_{4}, z_{4}) 
It can be found by solving the following determinant equation:
x^{2} + y^{2} + z^{2}  x_{ }  y_{ }  z_{ }  1
 x_{1}^{2} + y_{1}^{2} + z_{1}^{2}  x_{1}  y_{1}  z_{1}  1
 x_{2}^{2} + y_{2}^{2} + z_{2}^{2}  x_{2}  y_{2}  z_{2}  1
 x_{3}^{2} + y_{3}^{2} + z_{3}^{2}  x_{3}  y_{3}  z_{3}  1
 x_{4}^{2} + y_{4}^{2} + z_{4}^{2}  x_{4}  y_{4}  z_{4}  1

 = 0

This can be solved by evaluating the cofactors for the first row of the determinant. The determinant can be written as an equation of these cofactors:
(x^{2} + y^{2} + z^{2}) M_{11} − x M_{12} + y M_{13} − z M_{14} + M_{15} = 0 
Since, (x^{2} + y^{2} + z^{2}) = r^{2} this can be simplified to
r^{ 2} − x 
M_{12} 
+ y 
M_{13} 
− z 
M_{14} 
+ 
M_{15} 
= 0 
M_{11} 
M_{11} 
M_{11} 
M_{11} 
The general equation of a sphere with radius r_{0} and center (x_{0}, y_{0}, z_{0}) is
(x − x_{0})^{2} + (y − y_{0})^{2} + (z − z_{0})^{2} − r_{0}^{2} = 0 
Expanding this gives,
(x^{2} − 2 x x_{0} + x_{0}^{2}) + (y^{2} − 2 y y_{0} + y_{0}^{2}) + (z^{2} − 2 z z_{0} + z_{0}^{2}) − r_{0}^{2} = 0 
Rearranging terms and substitution gives,
r^{2} − 2 x x_{0} − 2 y y_{0} − 2 z z_{0} + x_{0}^{2} + y_{0}^{2} + z_{0}^{2} − r_{0}^{2} = 0 
Equating the like terms from the determinant equation and the general equation for the sphere gives:
x_{0} = + 0.5 
M_{12} 
M_{11} 

y_{0} = − 0.5 
M_{13} 
M_{11} 

z_{0} = + 0.5 
M_{14} 
M_{11} 

r_{0}^{2} = x_{0}^{2} + y_{0}^{2} + z_{0}^{2} − 
M_{15} 
M_{11} 

Note that there is no solution when M_{11} is equal to zero. In this case, the points are not on a sphere; they may all be on a plane or three point may be on a straight line.
Zeno source code
type POINT : record
x,y,z : real
end record
type FOURPOINTS : array[4] of POINT
type matrix : array[4,4] of real
program
var r : real
var c : POINT
var p : FOURPOINTS
p[1].x := 0
p[1].y := 4
p[1].z := 2
p[2].x := 1
p[2].y := 0
p[2].z := 1
p[3].x := 2
p[3].y := 2
p[3].z := 0
p[4].x := 5
p[4].y := 2
p[4].z := 3
r := sphere( c, p )
if r > 0 then
put "Sphere: (", c.x, ",", c.y, ",", c.z, "), ", r
else
put "Not a sphere!"
end if
end program
%
% Calculate center and radius of
% sphere given four points
%
function sphere( var c : POINT, var p : FOURPOINTS ) : real
var i : int
var r, m11, m12, m13, m14, m15 : real
var a : matrix
for i := 1...4 do % find minor 11
a[i,1] := p[i].x
a[i,2] := p[i].y
a[i,3] := p[i].z
a[i,4] := 1
end for
m11 := det( a, 4 )
for i := 1...4 do % find minor 12
a[i,1] := p[i].x^2 + p[i].y^2 + p[i].z^2
a[i,2] := p[i].y
a[i,3] := p[i].z
a[i,4] := 1
end for
m12 := det( a, 4 )
for i := 1...4 do % find minor 13
a[i,1] := p[i].x^2 + p[i].y^2 + p[i].z^2
a[i,2] := p[i].x
a[i,3] := p[i].z
a[i,4] := 1
end for
m13 := det( a, 4 )
for i := 1...4 do % find minor 14
a[i,1] := p[i].x^2 + p[i].y^2 + p[i].z^2
a[i,2] := p[i].x
a[i,3] := p[i].y
a[i,4] := 1
end for
m14 := det( a, 4 )
for i := 1...4 do % find minor 15
a[i,1] := p[i].x^2 + p[i].y^2 + p[i].z^2
a[i,2] := p[i].x
a[i,3] := p[i].y
a[i,4] := p[i].z
end for
m15 := det( a, 4 )
if m11 = 0 then
r := 0
else
c.x := 0.5 * m12 / m11 % center of sphere
c.y := 0.5 * m13 / m11
c.z := 0.5 * m14 / m11
r := sqrt( c.x^2 + c.y^2 + c.z^2  m15/m11 )
end if
return r % the radius
end function
%
% Calculate determinate using recursive
% expansion by minors.
%
function det( var a : matrix, n : int ) : real
var i, j, j1, j2 : int
var d : real := 0
var m : matrix
assert n > 1
if n = 2 then
d := a[1,1]*a[2,2]  a[2,1]*a[1,2]
else
d := 0
for j1 := 1...n do
% create minor
for i := 2...n do
j2 := 1
for j := 1...n do
continue when j = j1
m[i1,j2] := a[i,j]
incr j2
end for
end for
% calculate determinant
d := d + ( 1.0 )^(1 + j1 ) * a[1,j1] * det( m, n1 )
end for
end if
return d
end function
Copyright © 2005, Stephen R. Schmitt