This problem can be solved using the method of least squares. Suppose we know
N
N
N sliced circles, the th
i
i
The radius of i sliced circles is
r
i
r_i
ri?, the coordinates of the center of the circle are
(
x
i
,
the y
i
,
z
i
)
(x_i,y_i,z_i)
(xi?, yi?, zi?). The coordinates of the center of the sphere we want to fit are
(
x
0
,
the y
0
,
z
0
)
(x_0,y_0,z_0)
(x0?,y0?,z0?), with a radius of
r
r
r.
for the first
i
i
i sliced circles, we can write the equation of the plane in which they lie:
(
x
?
x
i
)
2
+
(
the y
?
the y
i
)
2
+
(
z
?
z
i
)
2
=
r
i
2
(x-x_i)^2 + (y-y_i)^2 + (z-z_i)^2=r_i^2
(x?xi?)2 + (y?yi?)2 + (z?zi?)2=ri2?
Simplified:
x
2
+
the y
2
+
z
2
?
2
x
x
i
?
2
the y
the y
i
?
2
z
z
i
+
(
x
i
2
+
the y
i
2
+
z
i
2
?
r
i
2
)
=
0
x^2 + y^2 + z^2-2xx_i-2yy_i-2zz_i + (x_i^2 + y_i^2 + z_i^2-r_i^2)=0
x2 + y2 + z2?2xxi2yyi2zzi? + (xi2? + yi2? + zi2ri2?)=0
We can write the above equation in vector form:
[
2
x
i
2
the y
i
2
z
i
?
1
]
[
x
the y
z
?
1
]
x
i
2
+
the y
i
2
+
z
i
2
?
r
i
2
\begin{bmatrix} 2x_i & amp; 2y_i & amp; 2z_i & amp; -1 \end{bmatrix} \begin{bmatrix} x \ y \ z \ -1 \end{bmatrix} x_i ^2 + y_i^2 + z_i^2-r_i^2
[2xi2yi2zi?1?][x y z ?1?]xi2? + yi2? + zi2ri2?
we can combine all
N
N
N plane equations are written in matrix form:
[
2
x
1
2
the y
1
2
z
1
?
1
2
x
2
2
the y
2
2
z
2
?
1
?
?
?
?
2
x
N
2
the y
N
2
z
N
?
1
]
[
x
the y
z
?
1
]
[
x
1
2
+
the y
1
2
+
z
1
2
?
r
1
2
x
2
2
+
the y
2
2
+
z
2
2
?
r
2
2
?
x
N
2
+
the y
N
2
+
z
N
2
?
r
N
2
]
\begin{bmatrix} 2x_1 & amp; 2y_1 & amp; 2z_1 & amp; -1 \ 2x_2 & amp; 2y_2 & amp; 2z_2 & amp; -1 \ \vdots & amp; \vdots & amp; \vdots & amp; \vdots \ 2x_N & amp; 2y_N & amp; 2z_N & amp; -1 \ \end{bmatrix} \begin{bmatrix} x \ y \ z \ -1 \end{bmatrix} \begin{bmatrix} x_1^2 + y_1^2 + z_1^2-r_1^2 \ x_2^2 + y_2^2 + z_2^2-r_2^2 \ \vdots \ \ x_N^2 + y_N^2 + z_N^2-r_N^2 \ \end{bmatrix}
[2x12y12z1?1 2x22y22z2?1 ? 2xN2yN2zN?1 ?][x y z ?1?][ + y12? + z12r12? x22? + y22? + z22r22? ? xN2? + yN2? + zN2rN2? ]
Let the above matrix be
A
A
A, the vector is
b
b
b, the unknown vector is
x
x
x, then the above equation can be written as
A
x
=
b
Ax=b
of the form Ax=b. Our goal is to minimize the residual by
∣
A
x
?
b
∣
|Ax-b|
∣Ax?b∣ to get
x
x
The optimal solution for x, which is the coordinates of the center of the sphere
(
x
0
,
the y
0
,
z
0
)
(x_0,y_0,z_0)
(x0?,y0?,z0?) and radius
r
r
r.
The solution of the least squares method is relatively simple, directly solving
x
=
(
A
T
A
)
?
1
A
T
b
x=(A^TA)^{-1}A^Tb
x=(ATA)?1ATb is enough. Will
x
x
x is split into
(
x
0
,
the y
0
,
z
0
,
r
)
(x_0,y_0,z_0,r)
(x0?,y0?,z0?,r), the coordinates and radius of the center of the sphere can be obtained.
|The center and radius of the sphere can be fitted using the method of least squares. The specific implementation steps are as follows:
1. Define a sphere structure, including the center coordinates and radius:
struct Sphere {<!-- -->
double x, y, z; // coordinates of the center of the sphere
double r; // radius
};
2. Define a function, the input is the radius and center of multiple sliced circles, and the output is the fitted sphere structure:
Sphere fitSphere(vector<pair<double, tuple<double, double, double>>> circles) {<!-- -->
// Initialize matrix A and vector b
Eigen::MatrixXd A(circles. size(), 4);
Eigen::VectorXd b(circles. size());
for (int i = 0; i < circles. size(); i ++ ) {<!-- -->
double r = circles[i].first;
double x = get<0>(circles[i].second);
double y = get<1>(circles[i].second);
double z = get<2>(circles[i].second);
A(i, 0) = 2 * x;
A(i, 1) = 2 * y;
A(i, 2) = 2 * z;
A(i, 3) = -1;
b(i) = x * x + y * y + z * z - r * r;
}
// Solve using the least squares method
Eigen::VectorXd x = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
// Convert the solution result to a sphere structure
Sphere sphere;
sphere.x = x(0);
sphere.y = x(1);
sphere.z = x(2);
sphere.r = std::sqrt(x(0) * x(0) + x(1) * x(1) + x(2) * x(2) - x(3));
return sphere;
}
3. Call the above function for fitting:
vector<pair<double, tuple<double, double, double>>> circles; // Radius and center of slice circle
// add sliced circle
circles.push_back(make_pair(1.0, make_tuple(0.0, 0.0, 0.0)));
circles.push_back(make_pair(2.0, make_tuple(0.0, 0.0, 1.0)));
circles.push_back(make_pair(3.0, make_tuple(0.0, 0.0, 2.0)));
// fit the sphere
Sphere sphere = fitSphere(circles);
// output result
cout << "Sphere center: (" << sphere.x << ", " << sphere.y << ", " << sphere.z << ")" << endl ;
cout << "Sphere radius: " << sphere.r << endl;