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:
It can be found by solving the following determinant equation:(x1, y1, z1), (x2, y2, z2), (x3, y3, z3), (x4, y4, z4)
x2 + y2 + z2 x y z 1 x12 + y12 + z12 x1 y1 z1 1 x22 + y22 + z22 x2 y2 z2 1 x32 + y32 + z32 x3 y3 z3 1 x42 + y42 + z42 x4 y4 z4 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:
(x2 + y2 + z2) M11 - x M12 + y M13 - z M14 + M15 = 0
Since, (x2 + y2 + z2) = r2 this can be simplified to
r2 - x M12 / M11 + y M13 / M11 - z M14 / M11 + M15 / M11 = 0
The general equation of a sphere with radius r0 and center (x0, y0, z0) is
(x - x0)2 + (y - y0)2 + (z - z0)2 - r02 = 0
Expanding this gives,
(x2 - 2 x x0 + x02) + (y2 - 2 y y0 + y02) + (z2 - 2 z z0 + z02) - r02 = 0
Re-arranging terms and substitution gives,
r2 - 2 x x0 - 2 y y0 - 2 z z0 + x02 + y02 + z02 - r02 = 0
Equating the like terms from the determinant equation and the general equation for the sphere gives:
x0 = + 0.5 M12 / M11 y0 = - 0.5 M13 / M11 z0 = + 0.5 M14 / M11 r02 = x02 + y02 + z02 - M15 / M11
Note that there is no solution when M11 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.
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[i-1,j2] := a[i,j]
incr j2
end for
end for
% calculate determinant
d := d + ( -1.0 )^(1 + j1 ) * a[1,j1] * det( m, n-1 )
end for
end if
return d
end function
AbCd Classics - free on-line books