## 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:

 (x1, y1, z1), (x2, y2, z2), (x3, y3, z3), (x4, y4, z4)

It can be found by solving the following determinant equation:

 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

 r 2 − x M12 + y M13 − z M14 + M15 = 0 M11 M11 M11 M11

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.

### 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

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
```