How to find the center and radius of an any-dimensional sphere given dims+1 points

Question

Given a vector of N-dimensional points. The vector will be of size N+1.

Is there a generalized algorithm to find the center and radius of the ND sphere using those points where the sphere intersects every single one of those points?

Phiphenomenon

Hi! Could you please add a hyphen between "N" and "dimensional"? Right now it's a bit confusing and can be interpreted as "given a vector of N points" which in the next sentence "will be of size N+1". – Retrospective

I’m voting to close this question because it is about math rather than programming – Atonsah

An answer can be found here: Does a set of n+1 points that affinely span R^n lie on a unique (n-1)-sphere?. The idea is simple: use the N+1 points to get a system of N+1 quadratic equations in the N+1 parameters of the sphere (center coordinates and radius). Use a variable change to transform this system into a linear system. You now have a linear system of N+1 equations in N+1 variables. Solve this system (by inverting the matrix). You now have the center coordinates and the radius. – Retrospective

Answer 1

The same question has been asked on the mathematics stackexchange and has received a constructive answer:

• Does a set of n+1 points that affinely span R^n lie on a unique (n-1)-sphere?

Here is an implementation in python/numpy of the algorithm described at that answer.

import numpy as np

def find_sphere_through_points(points):
    n_points, n_dim = points.shape
    if (n_points != n_dim + 1):
        raise ValueError('Number of points must be equal to 1 + dimension')
    a = np.concatenate((points, np.ones((n_points, 1))), axis=1)
    b = (points**2).sum(axis=1)
    x = np.linalg.solve(a, b)
    center = x[:-1] / 2
    radius = x[-1] + center@center
    return center, radius

To test this method, we can generate random points on the surface of a sphere, using the method described in this related question:

• Generate a random sample of points distributed on the surface of a unit spher

import numpy as np

def sample_spherical(npoints, ndim=3, center=None):
    vec = np.random.randn(npoints, ndim)
    vec /= np.linalg.norm(vec, axis=1).reshape(npoints,1)
    if center is None:
        return vec
    else:
        return vec + center

n = 5
center = np.random.rand(n)
points = sample_spherical(n+1, ndim=n, center=center)

guessed_center, guessed_radius = find_sphere_through_points(points)

print('True center:\n    ', center)
print('Calc center:\n    ', guessed_center)
print('True radius:\n    ', 1.0)
print('Calc radius:\n    ', guessed_radius)

# True center:
#      [0.18150032 0.94979547 0.07719378 0.26561175 0.37509931]
# Calc center:
#      [0.18150032 0.94979547 0.07719378 0.26561175 0.37509931]
# True radius:
#      1.0
# Calc radius:
#      0.9999999999999997

Retrospective

Answer 2

The center of a circle by three points, let (X, Y) and its radius R are found as follows:

(X - X0)² + (Y - Y0)² = R²
(X - X1)² + (Y - Y1)² = R²
(X - X2)² + (Y - Y2)² = R²

Then subtracting pair-wise to eliminate R,

(2X - X0 - X1)(X1 - X0) + (2Y - Y0 - Y1)(Y1 - Y0) = 0
(2X - X0 - X2)(X2 - X0) + (2Y - Y0 - Y2)(Y2 - Y0) = 0

This is a system of two linear equations in two unknowns, giving the coordinates of the center (in fact we construct the intersection of two bisectors). The radius follows from the first equation.

This immediately generalizes to D dimensions.

Pictor