OGS
MinimalBoundingSphere.cpp
Go to the documentation of this file.
1
16
17#include <Eigen/Geometry>
18#include <ctime>
19
21#include "MathLib/MathTools.h"
22#include "MathLib/Point3d.h"
23
24namespace GeoLib
25{
27
29 double radius)
30 : _radius(radius), _center(p)
31{
32}
33
35 MathLib::Point3d const& q)
36 : _radius(std::numeric_limits<double>::epsilon()), _center(p)
37{
38 Eigen::Vector3d const a = q.asEigenVector3d() - p.asEigenVector3d();
39
40 Eigen::Vector3d o = a / 2;
41 _radius = o.norm() + std::numeric_limits<double>::epsilon();
42 o += p.asEigenVector3d();
43 _center = MathLib::Point3d{{o[0], o[1], o[2]}};
44}
45
47 MathLib::Point3d const& q,
48 MathLib::Point3d const& r)
49{
50 auto const& vp = p.asEigenVector3d();
51 Eigen::Vector3d const a = r.asEigenVector3d() - vp;
52 Eigen::Vector3d const b = q.asEigenVector3d() - vp;
53 Eigen::Vector3d const axb = a.cross(b);
54
55 if (axb.squaredNorm() > 0)
56 {
57 double const denom = 2.0 * axb.dot(axb);
58 Eigen::Vector3d o =
59 (b.dot(b) * axb.cross(a) + a.dot(a) * b.cross(axb)) / denom;
60 _radius = o.norm() + std::numeric_limits<double>::epsilon();
61 o += vp;
62 _center = MathLib::Point3d{{o[0], o[1], o[2]}};
63 }
64 else
65 {
66 MinimalBoundingSphere two_pnts_sphere;
67 if (a.squaredNorm() > b.squaredNorm())
68 {
69 two_pnts_sphere = MinimalBoundingSphere(p, r);
70 }
71 else
72 {
73 two_pnts_sphere = MinimalBoundingSphere(p, q);
74 }
75 _radius = two_pnts_sphere.getRadius();
76 _center = two_pnts_sphere.getCenter();
77 }
78}
79
81 MathLib::Point3d const& q,
82 MathLib::Point3d const& r,
83 MathLib::Point3d const& s)
84{
85 Eigen::Vector3d const va = q.asEigenVector3d() - p.asEigenVector3d();
86 Eigen::Vector3d const vb = r.asEigenVector3d() - p.asEigenVector3d();
87 Eigen::Vector3d const vc = s.asEigenVector3d() - p.asEigenVector3d();
88
89 if (!MathLib::isCoplanar(p, q, r, s))
90 {
91 double const denom = 2.0 * va.cross(vb).dot(vc);
92 Eigen::Vector3d o =
93 (vc.dot(vc) * va.cross(vb) + vb.dot(vb) * vc.cross(va) +
94 va.dot(va) * vb.cross(vc)) /
95 denom;
96
97 _radius = o.norm() + std::numeric_limits<double>::epsilon();
98 o += p.asEigenVector3d();
99 _center = MathLib::Point3d({o[0], o[1], o[2]});
100 }
101 else
102 {
103 MinimalBoundingSphere const pqr(p, q, r);
104 MinimalBoundingSphere const pqs(p, q, s);
105 MinimalBoundingSphere const prs(p, r, s);
106 MinimalBoundingSphere const qrs(q, r, s);
107 _radius = pqr.getRadius();
108 _center = pqr.getCenter();
109 if (_radius < pqs.getRadius())
110 {
111 _radius = pqs.getRadius();
112 _center = pqs.getCenter();
113 }
114 if (_radius < prs.getRadius())
115 {
116 _radius = prs.getRadius();
117 _center = prs.getCenter();
118 }
119 if (_radius < qrs.getRadius())
120 {
121 _radius = qrs.getRadius();
122 _center = qrs.getCenter();
123 }
124 }
125}
126
128 std::vector<MathLib::Point3d*> const& points)
129 : _radius(-1), _center({0, 0, 0})
130{
131 const std::vector<MathLib::Point3d*>& sphere_points(points);
132 MinimalBoundingSphere const bounding_sphere =
133 recurseCalculation(sphere_points, 0, sphere_points.size(), 0);
134 _center = bounding_sphere.getCenter();
135 _radius = bounding_sphere.getRadius();
136}
137
139 std::vector<MathLib::Point3d*> sphere_points,
140 std::size_t start_idx,
141 std::size_t length,
142 std::size_t n_boundary_points)
143{
145 switch (n_boundary_points)
146 {
147 case 0:
148 sphere = MinimalBoundingSphere();
149 break;
150 case 1:
151 sphere = MinimalBoundingSphere(*sphere_points[start_idx - 1]);
152 break;
153 case 2:
154 sphere = MinimalBoundingSphere(*sphere_points[start_idx - 1],
155 *sphere_points[start_idx - 2]);
156 break;
157 case 3:
158 sphere = MinimalBoundingSphere(*sphere_points[start_idx - 1],
159 *sphere_points[start_idx - 2],
160 *sphere_points[start_idx - 3]);
161 break;
162 case 4:
163 sphere = MinimalBoundingSphere(
164 *sphere_points[start_idx - 1], *sphere_points[start_idx - 2],
165 *sphere_points[start_idx - 3], *sphere_points[start_idx - 4]);
166 return sphere;
167 }
168
169 for (std::size_t i = 0; i < length; ++i)
170 {
171 // current point is located outside of sphere
172 if (sphere.pointDistanceSquared(*sphere_points[start_idx + i]) > 0)
173 {
174 if (i > start_idx)
175 {
176 using DiffType =
177 std::vector<MathLib::Point3d*>::iterator::difference_type;
178 std::vector<MathLib::Point3d*> const tmp_ps(
179 sphere_points.cbegin() + static_cast<DiffType>(start_idx),
180 sphere_points.cbegin() +
181 static_cast<DiffType>(start_idx + i + 1));
182 std::copy(tmp_ps.cbegin(), --tmp_ps.cend(),
183 sphere_points.begin() +
184 static_cast<DiffType>(start_idx + 1));
185 sphere_points[start_idx] = tmp_ps.back();
186 }
187 sphere = recurseCalculation(sphere_points, start_idx + 1, i,
188 n_boundary_points + 1);
189 }
190 }
191 return sphere;
192}
193
195 MathLib::Point3d const& pnt) const
196{
197 return MathLib::sqrDist(_center, pnt) - (_radius * _radius);
198}
199
200} // namespace GeoLib
Calculation of a minimum bounding sphere for a vector of points.
Definition of the Point3d class.
static MinimalBoundingSphere recurseCalculation(std::vector< MathLib::Point3d * > sphere_points, std::size_t start_idx, std::size_t length, std::size_t n_boundary_points)
MinimalBoundingSphere()
Constructor using no points.
double getRadius() const
Returns the radius of the sphere.
MathLib::Point3d getCenter() const
Returns the center point of the sphere.
double pointDistanceSquared(MathLib::Point3d const &pnt) const
Eigen::Vector3d const & asEigenVector3d() const
Definition Point3d.h:64
bool isCoplanar(const MathLib::Point3d &a, const MathLib::Point3d &b, const MathLib::Point3d &c, const MathLib::Point3d &d)
Checks if the four given points are located on a plane.
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition Point3d.cpp:26