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

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?

c++mathtriangulation


Solution 1:[1]

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

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:

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

Solution 2:[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.

Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution Source
Solution 1 Stef
Solution 2

Related Questions