pub_soft/libcifpp
3
0
Fork 0
mirror of https://github.com/PDB-REDO/libcifpp.git synced 2026-06-04 13:54:25 +08:00
Code Issues Packages Projects Releases Wiki Activity

fixing smallest sphere

Browse Source
This commit is contained in:
Maarten L. Hekkelman
2025-11-05 13:18:58 +01:00
parent f94e9aece9
commit 251fb55d6a
1 changed files with 65 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files

131
src/point.cpp
View File

@@ -330,98 +330,97 @@ point nudge(point p, float offset)
// --------------------------------------------------------------------
std::tuple<point, float> smallest_sphere_around_points(point p1, point p2)
std::tuple<point, float> smallest_sphere_around_2_points(std::array<cif::point, 2> pts)
{
return { (p1 + p2) / 2, distance(p1, p2) / 2 };
return { (pts[0] + pts[1]) / 2, distance(pts[0], pts[1]) / 2 };
}
std::tuple<point, float> smallest_sphere_around_points(point p1, point p2, point p3)
std::tuple<point, float> smallest_sphere_around_3_points(std::array<cif::point, 3> pts)
{
// Find two bisectors
auto vz = cross_product(p2 - p1, p3 - p1);
auto vz = cross_product(pts[1] - pts[0], pts[2] - pts[0]);
auto bs1 = cross_product(vz, p2 - p1);
auto bs1 = cross_product(vz, pts[1] - pts[0]);
bs1.normalize();
auto v1 = (p2 - p1);
auto v1 = (pts[1] - pts[0]);
v1.normalize();
auto s1 = p1 + (distance(p2, p1) / 2) * v1;
auto s1 = pts[0] + (distance(pts[1], pts[0]) / 2) * v1;
auto bs2 = cross_product(vz, p3 - p1);
auto bs2 = cross_product(vz, pts[2] - pts[0]);
bs2.normalize();
auto v2 = (p3 - p1);
auto v2 = (pts[2] - pts[0]);
v2.normalize();
auto s2 = p1 + (distance(p3, p1) / 2) * v2;
auto s2 = pts[0] + (distance(pts[2], pts[0]) / 2) * v2;
auto c = line_line_intersection(s1, s1 + bs1, s2, s2 + bs2);
if (c)
return { *c, distance(*c, p1) };
return { *c, distance(*c, pts[0]) };
// Colinear points I guess, try something else
auto l1 = distance_squared(p1, p2);
auto l2 = distance_squared(p1, p3);
auto l3 = distance_squared(p2, p3);
auto l1 = distance_squared(pts[0], pts[1]);
auto l2 = distance_squared(pts[0], pts[2]);
auto l3 = distance_squared(pts[1], pts[2]);
if (l1 > l2 and l1 > l3)
return smallest_sphere_around_points(p1, p2);
return smallest_sphere_around_2_points({ pts[0], pts[1] });
else if (l2 > l1 and l2 > l3)
return smallest_sphere_around_points(p1, p3);
return smallest_sphere_around_2_points({ pts[0], pts[2] });
else
return smallest_sphere_around_points(p2, p3);
return smallest_sphere_around_2_points({ pts[1], pts[2] });
}
std::tuple<point, float> smallest_sphere_around_points(point p1, point p2, point p3, point p4)
std::tuple<point, float> smallest_sphere_around_4_points(std::array<cif::point, 4> pts)
{
auto t0 = -norm_squared(p1);
auto t1 = -norm_squared(p2);
auto t2 = -norm_squared(p3);
auto t3 = -norm_squared(p4);
auto t0 = -norm_squared(pts[0]);
auto t1 = -norm_squared(pts[1]);
auto t2 = -norm_squared(pts[2]);
auto t3 = -norm_squared(pts[3]);
// clang-format off
matrix4x4<float> Tm({
p1.m_x, p1.m_y, p1.m_z, 1,
p2.m_x, p2.m_y, p2.m_z, 1,
p3.m_x, p3.m_y, p3.m_z, 1,
p4.m_x, p4.m_y, p4.m_z, 1
pts[0].m_x, pts[0].m_y, pts[0].m_z, 1,
pts[1].m_x, pts[1].m_y, pts[1].m_z, 1,
pts[2].m_x, pts[2].m_y, pts[2].m_z, 1,
pts[3].m_x, pts[3].m_y, pts[3].m_z, 1
});
auto T = determinant(Tm);
if (T != 0)
{
matrix4x4<float> Dm({
t0, p1.m_y, p1.m_z, 1,
t1, p2.m_y, p2.m_z, 1,
t2, p3.m_y, p3.m_z, 1,
t3, p4.m_y, p4.m_z, 1
t0, pts[0].m_y, pts[0].m_z, 1,
t1, pts[1].m_y, pts[1].m_z, 1,
t2, pts[2].m_y, pts[2].m_z, 1,
t3, pts[3].m_y, pts[3].m_z, 1
});
auto D = determinant(Dm) / T;
matrix4x4<float> Em({
p1.m_x, t0, p1.m_z, 1,
p2.m_x, t1, p2.m_z, 1,
p3.m_x, t2, p3.m_z, 1,
p4.m_x, t3, p4.m_z, 1
pts[0].m_x, t0, pts[0].m_z, 1,
pts[1].m_x, t1, pts[1].m_z, 1,
pts[2].m_x, t2, pts[2].m_z, 1,
pts[3].m_x, t3, pts[3].m_z, 1
});
auto E = determinant(Em) / T;
matrix4x4<float> Fm({
p1.m_x, p1.m_y, t0, 1,
p2.m_x, p2.m_y, t1, 1,
p3.m_x, p3.m_y, t2, 1,
p4.m_x, p4.m_y, t3, 1
pts[0].m_x, pts[0].m_y, t0, 1,
pts[1].m_x, pts[1].m_y, t1, 1,
pts[2].m_x, pts[2].m_y, t2, 1,
pts[3].m_x, pts[3].m_y, t3, 1
});
auto F = determinant(Fm) / T;
matrix4x4<float> Gm({
p1.m_x, p1.m_y, p1.m_z, t0,
p2.m_x, p2.m_y, p2.m_z, t1,
p3.m_x, p3.m_y, p3.m_z, t2,
p4.m_x, p4.m_y, p4.m_z, t3
pts[0].m_x, pts[0].m_y, pts[0].m_z, t0,
pts[1].m_x, pts[1].m_y, pts[1].m_z, t1,
pts[2].m_x, pts[2].m_y, pts[2].m_z, t2,
pts[3].m_x, pts[3].m_y, pts[3].m_z, t3
});
auto G = determinant(Gm) / T;
@@ -435,19 +434,19 @@ std::tuple<point, float> smallest_sphere_around_points(point p1, point p2, point
// Perhaps some colinear points, try something else:
// for (auto ix : std::initializer_list<std::array<size_t, 4>>{
// { 1, 2, 3, 0 },
// { 0, 2, 3, 1 },
// { 0, 1, 3, 2 },
// { 0, 1, 2, 3 },
// })
// {
// auto [center, radius] =
// smallest_sphere_around_3_points(pts[ix[0]], pts[ix[1]], pts[ix[2]]);
for (auto ix : std::initializer_list<std::array<size_t, 4>>{
{ 1, 2, 3, 0 },
{ 0, 2, 3, 1 },
{ 0, 1, 3, 2 },
{ 0, 1, 2, 3 },
})
{
auto [center, radius] =
smallest_sphere_around_3_points({ pts[ix[0]], pts[ix[1]], pts[ix[2]] });
// if (distance(pts[ix[3]], center) <= radius)
// return { center, radius };
// }
if (distance(pts[ix[3]], center) <= radius)
return { center, radius };
}
assert(false);
exit(1);
@@ -463,13 +462,13 @@ std::tuple<point, float> smallest_sphere_around_all_points(std::vector<point> P,
return { R[0], 0 };
case 2:
return smallest_sphere_around_points(R[0], R[1]);
return smallest_sphere_around_2_points({ R[0], R[1] });
case 3:
return smallest_sphere_around_points(R[0], R[1], R[2]);
return smallest_sphere_around_3_points({ R[0], R[1], R[2] });
case 4:
return smallest_sphere_around_points(R[0], R[1], R[2], R[3]);
return smallest_sphere_around_4_points({ R[0], R[1], R[2], R[3] });
default:
assert(false);
@@ -500,19 +499,19 @@ bool point_in_circle(point p, std::vector<point> c)
case 2:
{
auto [center, radius] = smallest_sphere_around_points(c[0], c[1]);
auto [center, radius] = smallest_sphere_around_2_points({ c[0], c[1] });
return cif::distance_squared(p, center) <= radius * radius;
}
case 3:
{
auto [center, radius] = smallest_sphere_around_points(c[0], c[1], c[2]);
auto [center, radius] = smallest_sphere_around_3_points({ c[0], c[1], c[2] });
return cif::distance_squared(p, center) <= radius * radius;
}
case 4:
{
auto [center, radius] = smallest_sphere_around_points(c[0], c[1], c[2], c[3]);
auto [center, radius] = smallest_sphere_around_4_points({ c[0], c[1], c[2], c[3] });
return cif::distance_squared(p, center) <= radius * radius;
}
@@ -550,13 +549,13 @@ std::tuple<point, float> smallest_sphere_around_points(std::vector<point> pts)
else
{
cix.erase(std::remove_if(cix.begin(), cix.end(), [i](size_t j)
{ return j < i; }),
{ return j < i; }),
cix.end());
cix.push_back(i);
if (cix.size() < 4)
i = 0;
else
++i;
++i;
}
}
@@ -565,11 +564,11 @@ std::tuple<point, float> smallest_sphere_around_points(std::vector<point> pts)
case 1:
return { pts[cix[0]], 0 };
case 2:
return smallest_sphere_around_points(pts[cix[0]], pts[cix[1]]);
return smallest_sphere_around_2_points({ pts[cix[0]], pts[cix[1]] });
case 3:
return smallest_sphere_around_points(pts[cix[0]], pts[cix[1]], pts[cix[2]]);
return smallest_sphere_around_3_points({ pts[cix[0]], pts[cix[1]], pts[cix[2]] });
case 4:
return smallest_sphere_around_points(pts[cix[0]], pts[cix[1]], pts[cix[2]], pts[cix[3]]);
return smallest_sphere_around_4_points({ pts[cix[0]], pts[cix[1]], pts[cix[2]], pts[cix[3]] });
default:
assert(false);
throw std::runtime_error("Error finding smallest sphere");