How would I go about calculating the center of a sphere based on points from it's surface?

Question

I've recorded random Cartesian coordinates from the surface of a sphere. I'm trying to extrapolate where the center x,y,z coordinate of this sphere would be based on the points below.

I'm trying to create a script to solve this, but the math has me stuck.

Appreciate any help. Thanks!

Coordinates x,y,z:
-1.190185,0.7824033,0.1691585
-1.152931,0.7811859,0.1401751
-1.110813,0.775885, 0.1348239
-1.078485,0.772419, 0.1545465
-1.101839,0.7839018,0.2090827
-1.120642,0.7859643,0.2231803
-1.136541,0.7895643,0.2040433
-1.135001,0.7879777,0.1715891
-1.111366,0.7854245,0.1738369
-1.115193,0.78763 , 0.1931185
-1.130671,0.7892017,0.1985596
-1.132487,0.7872546,0.2240417
-1.161051,0.7876478,0.171136
-1.114549,0.7809206,0.1476536
-1.083838,0.7764122,0.1626304
-1.06429, 0.769419,0.2177327
-1.142459,0.7805369,0.2528498
-1.162153,0.7863876,0.1668259

Answer 1

The equation of a sphere with center $(x_0, y_0, z_0)$ and radius $r$ is

$ (x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2 $

Expanding, we get,

$ x^2 + y^2 + z^2 + A x + B y + C z + D = 0 $

where

$ A = -2 x_0 , \ B = - 2 y_0, \ C = - 2 z_0, \ D = x_0^2 + y_0^2 + z_0^2 - r^2 $

To identify the parameters $A,B,C,D$, we build a linear regression model, where we set

$ A x_i + B y_i + C z_i + D = -( x_i^2 + y_i^2 + z_i^2) $

If the above model equation is repeated for $i = 1, 2,.., N $, then we end up with the following regression model

$ M X = Y $

where the $i$-th row of $M$ is

$M_i = [x_i , y_i, z_i, 1]$

and $X$ is the vector of parameters:

$ X = [A, B, C, D]^T$

and $Y$ is the data vector, whose $i$-th entry is

$Y_i = -( x_i^2 + y_i^2 + z_i^2) $

The least squares estimate of $X$ is given by

$ \hat{X} = (M^T M)^{-1} M^T Y $

The computation of $\hat{X}$ involves the inversion of the $4 \times 4$ matrix $M^T M$.

Once $\hat{X}$ is determined, then we know the estimate of $A, B, C, D$, and hence,

$x_0 = - \dfrac{1}{2} A$

$y_0 = - \dfrac{1}{2} B $

$z_0 = -\dfrac{1}{2} C $

$ r^2 = x_0^2 + y_0^2 + z_0^2 - D$

Answer 2

The equation of a 3D sphere is $$(x-x_0)^2+(y-y_0)^2+(z-z_0)^2=r^2$$, where $(x_0, y_0, z_0)$ is the center of the sphere and $r$ is the radius. Plug in any four noncoplanar points $(x,y,z)$ to get a system of equations with four unknowns $x_0, y_0, z_0, r$ in order to solve for $z$. Of course, this system is quadratic, so it would be tricky to solve by hand, but it can be solved exactly either by means of an existing computer program or approximated through numerical means.

Answer 3

One way to do it is the following:

Take any three distinct points $A,B,C,$ and find the circumcentre of those $3$ points. To do this, just find the perpendicular bisector of $AB$ and then find the perpendicular bisector of $BC.$ The circumcentre of triangle $ABC$ is the intersection of these two perpendicular bisectors.

Then the centre of the sphere lies of the line perpendicular to the plane formed by the three points and that passes through the circumcentre of $ABC.$

Now do this again with three other points $DEF$: and where the line perpendicular to triangle $ABC$ and the line perpendicular to triangle $DEF$ intersect is the centre of the sphere.

Answer 4

Least Squared Spherical Regression

Port of a MATLAB solution in Mathematica.

Clear["Global`*"];
(*Given coordinates*)
data = {{-1.190185, 0.7824033, 0.1691585}, {-1.152931, 0.7811859, 
    0.1401751}, {-1.110813, 0.775885, 0.1348239}, {-1.078485, 
    0.772419, 0.1545465}, {-1.101839, 0.7839018, 
    0.2090827}, {-1.120642, 0.7859643, 0.2231803}, {-1.136541, 
    0.7895643, 0.2040433}, {-1.135001, 0.7879777, 
    0.1715891}, {-1.111366, 0.7854245, 0.1738369}, {-1.115193, 
    0.78763, 0.1931185}, {-1.130671, 0.7892017, 
    0.1985596}, {-1.132487, 0.7872546, 0.2240417}, {-1.161051, 
    0.7876478, 0.171136}, {-1.114549, 0.7809206, 
    0.1476536}, {-1.083838, 0.7764122, 0.1626304}, {-1.06429, 
    0.769419, 0.2177327}, {-1.142459, 0.7805369, 
    0.2528498}, {-1.162153, 0.7863876, 0.1668259}};


sphereFit[data_] := 
 Module[{A, B, Center, Radius, meanX, meanY, 
   meanZ},(*Compute the means*)meanX = Mean[data[[All, 1]]];
  meanY = Mean[data[[All, 2]]];
  meanZ = Mean[data[[All, 3]]];
  (*Compute matrix A*)
  A = {{Mean[(data[[All, 1]] - meanX)*data[[All, 1]]], 
     2*Mean[data[[All, 1]]*(data[[All, 2]] - meanY)], 
     2*Mean[data[[All, 1]]*(data[[All, 3]] - meanZ)]}, {0, 
     Mean[(data[[All, 2]] - meanY)*data[[All, 2]]], 
     2*Mean[data[[All, 2]]*(data[[All, 3]] - meanZ)]}, {0, 0, 
     Mean[(data[[All, 3]] - meanZ)*data[[All, 3]]]}};
  A = A + Transpose[A];
  (*Compute vector B*)
  B = {Mean[(data[[All, 1]]^2 + data[[All, 2]]^2 + 
        data[[All, 3]]^2)*(data[[All, 1]] - meanX)], 
    Mean[(data[[All, 1]]^2 + data[[All, 2]]^2 + 
        data[[All, 3]]^2)*(data[[All, 2]] - meanY)], 
    Mean[(data[[All, 1]]^2 + data[[All, 2]]^2 + 
        data[[All, 3]]^2)*(data[[All, 3]] - meanZ)]};
  (*Compute the center*)Center = LinearSolve[A, B];
  (*Compute the radius*)
  Radius = Sqrt[
    Mean[Total[(data - ConstantArray[Center, Length[data]])^2, {2}]]];
  (*Return center and radius*){Center, Radius}]

{center, radius} = sphereFit[data];
Print["The Least Squares(LR) Sphere: 
Center=", center, " Radius=", radius];

(*Visualize data and best fit sphere*)
Show[Graphics3D[{PointSize[0.02], Point[data], Opacity[0.3], 
   Sphere[center, radius]}], Axes -> True, 
 AxesLabel -> {"x", "y", "z"}]

Answer 5

You might be able to simplify the algorithm. Any three non-collinear points determines a circle. 4 non-coplanar points determines a sphere.

You can combine these ideas together. Pick three points, find the circle going through them. Pick another three points and do the same.

The normal to two circles made by planes intersecting the same sphere intersects at the center of the sphere. The sum of the square of distance from that intersection to the center of the circle and the square of the radius of that circle gives you the square of the radius of the sphere.

The perpendicular bisectors of the legs of a triangle intersect at the center of their circumscribed circle.

You can find the perpendicular bisectors using cross products. In fact that will also give you the normals to the planes you need.

So have a subroutine that can convert coordinates to vectors and another subroutine that can take the vector cross product of those vectors. Should allow you to calculate the center and the radius of the sphere.