Share | home | download page | send comment | send link | add bookmark |

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 + z12x1y1z11
x22 + y22 + z22x2y2z21
x32 + y32 + z32x3y3z31
x42 + y42 + z42x4y4z41
 = 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) M11x M12 + y M13z M14 + M15 = 0

Since, (x2 + y2 + z2) = r2 this can be simplified to

r 2 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

(xx0)2 + (yy0)2 + (zz0)2r02 = 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 + z02r02 = 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

    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

Share | home | download page | send comment | send link | add bookmark |


Copyright © 2005, Stephen R. Schmitt