[util] Make a number of improvements to SEB

Attempting to vis ad_tears drags a few lurking bugs out of
SmallestEnclosingBall_vf: poor calculation of 2-point affine space, poor
handling of duplicate points and dropped support points, poor
calculation of the new center (related to duplicate points), and
insufficient iterations for large point sets. qfvis (modified for
cluster spheres) now loads ad_tears.
This commit is contained in:
Bill Currie 2021-07-30 13:50:12 +09:00
parent 9461779ba7
commit 9d819254d4
4 changed files with 399 additions and 62 deletions

View file

@ -35,6 +35,8 @@
# include <strings.h>
#endif
#include "qfalloca.h"
#include <math.h>
#define IMPLEMENT_R_Cull
@ -42,6 +44,7 @@
#include "QF/mathlib.h"
#include "QF/qtypes.h"
#include "QF/set.h"
#include "QF/sys.h"
VISIBLE int nanmask = 255 << 23;
@ -1348,12 +1351,13 @@ test_support_points(const vec_t **points, int *num_points, const vec3_t center)
break;
case 2:
VectorSubtract (points[1], points[0], v);
VectorSubtract (center, points[0], d);
VectorAdd (points[0], points[1], d);
VectorScale (d, 0.5, d);
VectorSubtract (center, d, d);
CrossProduct (v, d, n);
nn = DotProduct (n, n);
dd = DotProduct (d, d);
vv = DotProduct (v, v);
in_affine = nn < 1e-5 * vv * dd;
in_affine = nn < 1e-5 * vv * vv;
break;
case 3:
VectorSubtract (points[1], points[0], a);
@ -1403,12 +1407,14 @@ test_support_points(const vec_t **points, int *num_points, const vec3_t center)
sphere_t
SmallestEnclosingBall (const vec3_t points[], int num_points)
{
set_t was_support = SET_STATIC_INIT (num_points, alloca);
sphere_t sphere;
const vec_t *best;
const vec_t *support[4];
int num_support;
vec_t dist, best_dist;
int i;
int best_i = 0;
int iters = 0;
if (num_points < 1) {
@ -1419,6 +1425,7 @@ SmallestEnclosingBall (const vec3_t points[], int num_points)
for (i = 0; i < 4; i++)
support[i] = 0;
set_empty (&was_support);
VectorCopy (points[0], sphere.center);
best_dist = dist = 0;
@ -1427,55 +1434,51 @@ SmallestEnclosingBall (const vec3_t points[], int num_points)
dist = VectorDistance_fast (points[i], sphere.center);
if (dist > best_dist) {
best_dist = dist;
best_i = i;
best = points[i];
}
}
num_support = 1;
support[0] = best;
sphere.radius = best_dist; // note: radius squared until the end
set_add (&was_support, best_i);
while (!test_support_points (support, &num_support, sphere.center)) {
vec3_t affine;
vec3_t center_to_affine, center_to_point;
vec_t affine_dist, point_proj, point_dist, bound;
vec_t scale = 1;
vec3_t affine, v, p, r;
vec_t x, best_x = 0, rr, pv, pp;
int i;
if (iters++ > 10)
if (iters++ > 2 * num_points)
Sys_Error ("stuck SEB");
closest_affine_point (support, num_support, sphere.center, affine);
VectorSubtract (affine, sphere.center, center_to_affine);
affine_dist = DotProduct (center_to_affine, center_to_affine);
if (affine_dist < sphere.radius * 1e-5) {
// It's possible test_support_points failed due to precision
// issues
break;
}
VectorSubtract (support[0], affine, r);
rr = DotProduct (r, r);
VectorSubtract (sphere.center, affine, v);
best = 0;
for (i = 0; i < num_points; i++) {
if (points[i] == support[0] || points[i] == support[1]
|| points[i] == support[2])
if (SET_TEST_MEMBER (&was_support, i)) {
continue;
VectorSubtract (points[i], sphere.center, center_to_point);
point_proj = DotProduct (center_to_affine, center_to_point);
if (affine_dist - point_proj <= 0
|| ((affine_dist - point_proj) * (affine_dist - point_proj)
< 1e-6 * sphere.radius * affine_dist))
}
VectorSubtract (points[i], affine, p);
pp = DotProduct (p, p);
pv = DotProduct (p, v);
if (pp <= rr || pv <= 0 || pv * pv < 1e-6 * rr) {
continue;
point_dist = DotProduct (center_to_point, center_to_point);
bound = sphere.radius - point_dist;
bound /= 2 * (affine_dist - point_proj);
if (bound < scale) {
}
x = (pp - rr) / (2 * pv);
if (x > best_x) {
best = points[i];
scale = bound;
best_i = i;
best_x = x;
}
}
VectorMultAdd (sphere.center, scale, center_to_affine, sphere.center);
VectorMultAdd (affine, best_x, v, sphere.center);
sphere.radius = VectorDistance_fast (sphere.center, support[0]);
if (best) {
support[num_support++] = best;
set_add (&was_support, best_i);
}
}
best_dist = 0;

View file

@ -28,6 +28,8 @@
# include "config.h"
#endif
#include "qfalloca.h"
#include <math.h>
#define IMPLEMENT_VEC4F_Funcs
@ -37,6 +39,7 @@
#include "QF/simd/vec4d.h"
#include "QF/simd/vec4f.h"
#include "QF/simd/mat4f.h"
#include "QF/set.h"
#include "QF/sys.h"
vec4f_t
@ -192,12 +195,11 @@ test_support_points(const vec4f_t **points, int *num_points, vec4f_t center)
break;
case 2:
v = *points[1] - *points[0];
d = center - *points[0];
d = center - (*points[0] + *points[1]) / 2;
n = crossf (v, d);
nn = dotf (n, n)[0];
dd = dotf (d, d)[0];
vv = dotf (v, v)[0];
in_affine = nn < 1e-5 * vv * dd;
in_affine = nn <= 1e-5 * vv * vv;
break;
case 3:
a = *points[1] - *points[0];
@ -247,6 +249,7 @@ test_support_points(const vec4f_t **points, int *num_points, vec4f_t center)
vspheref_t
SmallestEnclosingBall_vf (const vec4f_t *points, int num_points)
{
set_t was_support = SET_STATIC_INIT (num_points, alloca);
vspheref_t sphere = {};
vec4f_t center = {};
const vec4f_t *best;
@ -262,6 +265,7 @@ SmallestEnclosingBall_vf (const vec4f_t *points, int num_points)
for (i = 0; i < 4; i++) {
support[i] = 0;
}
set_empty (&was_support);
vec4f_t dist = {};
float best_dist = 0;
@ -278,50 +282,51 @@ SmallestEnclosingBall_vf (const vec4f_t *points, int num_points)
num_support = 1;
support[0] = best;
sphere.radius = best_dist; // note: radius squared until the end
set_add (&was_support, best - points);
while (!test_support_points (support, &num_support, center)) {
vec4f_t affine;
vec4f_t center_to_affine, center_to_point;
float affine_dist, point_proj, point_dist, bound;
float scale = 1;
vec4f_t affine, r, rr;
vec4f_t v, p, pv, pp, x;
vec4f_t best_x = { };
int i;
if (iters++ > 10)
//Sys_Printf ("%d: "VEC4F_FMT", %.9g, %d\n", iters, VEC4_EXP (center), sqrt(sphere.radius), num_support);
if (iters++ > 2 * num_points) {
//for (i = 0; i < num_points; i++) {
// Sys_Printf (VEC4F_FMT",\n", VEC4_EXP (points[i]));
//}
Sys_Error ("stuck SEB");
}
affine = closest_affine_point (support, num_support, center);
center_to_affine = affine - center;
affine_dist = dotf (center_to_affine, center_to_affine)[0];
if (affine_dist < sphere.radius * 1e-5) {
// It's possible test_support_points failed due to precision
// issues
break;
}
r = *support[0] - affine; //FIXME better choice
rr = dotf (r, r);
v = center - affine;
best = 0;
for (i = 0; i < num_points; i++) {
if (&points[i] == support[0] || &points[i] == support[1]
|| &points[i] == support[2])
if (SET_TEST_MEMBER (&was_support, i)) {
continue;
center_to_point = points[i] - center;
point_proj = dotf (center_to_affine, center_to_point)[0];
if (affine_dist - point_proj <= 0
|| ((affine_dist - point_proj) * (affine_dist - point_proj)
< 1e-6 * sphere.radius * affine_dist))
}
p = points[i] - affine;
pp = dotf (p, p);
pv = dotf (p, v);
if (pp[0] <= rr[0] || pv[0] <= 0
|| (pv[0] * pv[0]) < 1e-6 * rr[0]) {
continue;
point_dist = dotf (center_to_point, center_to_point)[0];
bound = sphere.radius - point_dist;
bound /= 2 * (affine_dist - point_proj);
if (bound < scale) {
}
x = (pp - rr) / (2 * pv);
if (x[0] > best_x[0]) {
best = &points[i];
scale = bound;
best_x = x;
}
}
center = center + scale * center_to_affine;
center = affine + best_x * v;
dist = center - *support[0];
sphere.radius = dotf (dist, dist)[0];
if (best) {
support[num_support++] = best;
set_add (&was_support, best - points);
}
}
best_dist = 0;

View file

@ -6,6 +6,8 @@
#include "QF/mersenne.h"
#include "QF/sys.h"
#define SIZEOF(x) (sizeof(x) / sizeof(x[0]))
const vec3_t points[] = {
{-1, -1, 1},
{ 1, 1, 1},
@ -24,9 +26,155 @@ const vec3_t points[] = {
// due to qfbsp not culling the portal for being too small.
const vec3_t tears_triangle[] = {
{2201.82007, -1262, -713.450012},
{2201.8501, -1262, -713.593994},
{2201.8501, -1262, -713.593994},
{2201.84009, -1262, -713.445007},
};
const vec3_t tears_triangle2[] = {
{-519.989014, 1111.94995, -2592.05005},
{-520.002991, 1112, -2592.02002},
{-520.005005, 1112, -2592.04004},
};
const vec3_t tears_triangle3[] = {
{284, 624, 2959.76001},
{284, 624, 2960},
{284.234985, 624.471008, 2959.53003},
};
const vec3_t tears_quad[] = {
{1792, 1155.53003, 240},
{1792, 1030.32996, 240},
{1792.01001, 1030.32996, 240},
{1792.01001, 1155.53003, 240},
};
const vec3_t tears_quad2[] = {
{304, 3248, 3208},
{-61.2307014, 3248, 3208},
{-61.2307014, 3248, 3192},
{304, 3248, 3192},
};
const vec3_t tears_cluster[] = {
{551.638, -439.277008, 2804.63989},
{551, -438, 2806.02002},
{555.999023, -447.998993, 2805.1499},
{555.999023, -447.998993, 2804.32007},
{551.638, -439.277008, 2804.63989},
{555.999023, -447.998993, 2804.32007},
{555.999023, -447.998993, 2795.19995},
{550.97699, -437.954987, 2806.02002},
{551, -438, 2806.02002},
{555.999023, -447.998993, 2795.19995},
{555.999023, -447.998993, 2794.6499},
{550.978027, -437.954987, 2806.02002},
{545.945984, -442.945007, 2804.78003},
{546, -442.998993, 2804.77002},
{546, -443, 2804.69995},
{546, -443, 2804.78003},
{545.945007, -442.945007, 2804.78003},
{556, -448, 2794.6499},
{556, -448, 2805.1499},
{546, -443, 2804.78003},
{546, -443, 2804.69995},
{550.97699, -437.954987, 2806.02002},
{556, -448, 2794.6499},
{546.000977, -443, 2804.69995},
{545.945007, -442.945007, 2804.78003},
};
const vec3_t tears_cluster2[] = {
{-2160, -192, -704.242004},
{-2160, -192, -662},
{-2160, -196, -658},
{-2160, -220, -646},
{-2160, -222.332001, -645.416992},
{-2160, -320, -752},
{-2160, -192, -752},
{-2160, -192, -704.242004},
{-2160, -222.332001, -645.416992},
{-2160, -244, -640},
{-2160, -268, -640},
{-2160, -292, -646},
{-2160, -316, -658},
{-2160, -320, -662},
{-2160, -196, -814},
{-2160, -192, -810},
{-2160, -192, -752},
{-2160, -320, -752},
{-2160, -320, -810},
{-2160, -316, -814},
{-2160, -312, -816},
{-2160, -200, -816},
{-2174, -334, -676},
{-2186, -346, -700},
{-2192, -352, -724},
{-2192, -352, -748},
{-2186, -346, -772},
{-2174, -334, -796},
{-2160, -319.998993, -810},
{-2160, -319.998993, -661.999023},
{-2174, -178, -796},
{-2186, -166, -772},
{-2192, -160, -748},
{-2192, -160, -724},
{-2186, -166, -700},
{-2174, -178, -676},
{-2160, -192, -662},
{-2160, -192, -810},
{-2224, -320, -810},
{-2224, -330, -800},
{-2224, -320, -800},
{-2224, -200, -816},
{-2224, -312, -816},
{-2224, -316, -814},
{-2224, -320, -810},
{-2224, -320, -800},
{-2224, -192, -800},
{-2224, -192, -810},
{-2224, -196, -814},
{-2224, -192, -800},
{-2224, -182, -800},
{-2224, -192, -810},
{-2224, -330, -800},
{-2224, -334, -796},
{-2224, -178, -796},
{-2224, -182, -800},
{-2224, -334, -796},
{-2224, -346, -772},
{-2224, -351, -752},
{-2224, -161, -752},
{-2224, -166, -772},
{-2224, -178, -796},
{-2224, -351, -752},
{-2224, -352, -748},
{-2224, -352, -724},
{-2224, -351, -720},
{-2224, -161, -720},
{-2224, -160, -724},
{-2224, -160, -748},
{-2224, -161, -752},
{-2224, -351, -720},
{-2224, -346, -700},
{-2224, -334, -676},
{-2224, -316, -658},
{-2224, -292, -646},
{-2224, -268, -640},
{-2224, -244, -640},
{-2224, -220, -646},
{-2224, -196, -658},
{-2224, -178, -676},
{-2224, -166, -700},
{-2224, -161, -720},
};
// random numbers dug this up. needed more than 6 iterations
const vec3_t cloud_points[] = {
{1.70695281, 2.60889673, 2.86233163},
{2.11825109, 2.59883952, 1.48193693},
{2.76728868, 1.57646465, 2.39769602},
{1.88795376, 2.73980999, 1.58031297},
{2.79736757, 1.89582968, 1.80514407},
};
struct {
const vec3_t *points;
@ -38,9 +186,24 @@ struct {
{points, 2, {{ 0, 0, 1}, 1.41421356}},
{points, 3, {{-0.333333343, 0.333333343, 0.333333343}, 1.63299322}},
{points, 4, {{0, 0, 0}, 1.73205081}},
{tears_triangle, 3, {{2201.84521, -1262, -713.519531}, 0.0747000724}},
{tears_triangle, SIZEOF (tears_triangle),
{{2201.84521, -1262, -713.519531}, 0.0747000724}},//FIXME numbers?
{tears_cluster, SIZEOF (tears_cluster),
{{552.678101, -443.460876, 2800.40405}, 8.04672813}},//FIXME numbers?
{tears_quad, SIZEOF (tears_quad),
{{1792.005, 1092.92999, 240}, 62.6000302}},
{tears_triangle2, SIZEOF (tears_triangle2),
{{-519.995972, 1111.97498, -2592.03516}, 0.0300767105}},//FIXME numbers?
{tears_quad2, SIZEOF (tears_quad2),
{{121.384659, 3248, 3200}, 182.790497}},
{tears_triangle3, SIZEOF (tears_triangle3),
{{284.117493, 624.235535, 2959.76489}, 0.352925777}},
{tears_cluster2, SIZEOF (tears_cluster2),
{{-2192, -256.000031, -736}, 103.479492}},//FIXME numbers?
{cloud_points, SIZEOF (cloud_points),
{{2.20056152, 2.23369908, 2.2332375}, 0.88327992}},
};
#define num_tests (sizeof (tests) / sizeof (tests[0]))
#define num_tests SIZEOF (tests)
static inline float
rnd (mtstate_t *mt)

View file

@ -7,6 +7,8 @@
#include "QF/sys.h"
#include "QF/simd/vec4f.h"
#define SIZEOF(x) (sizeof(x) / sizeof(x[0]))
const vec4f_t points[] = {
{-1, -1, 1, 1},
{ 1, 1, 1, 1},
@ -25,23 +27,186 @@ const vec4f_t points[] = {
// due to qfbsp not culling the portal for being too small.
const vec4f_t tears_triangle[] = {
{2201.82007, -1262, -713.450012, 1},
{2201.8501, -1262, -713.593994, 1},
{2201.8501, -1262, -713.593994, 1},
{2201.84009, -1262, -713.445007, 1},
};
const vec4f_t tears_triangle2[] = {
{-519.989014, 1111.94995, -2592.05005, 1},
{-520.002991, 1112, -2592.02002, 1},
{-520.005005, 1112, -2592.04004, 1},
};
const vec4f_t tears_triangle3[] = {
{284, 624, 2959.76001, 1},
{284, 624, 2960, 1},
{284.234985, 624.471008, 2959.53003, 1},
};
const vec4f_t tears_quad[] = {
{1792, 1155.53003, 240, 1},
{1792, 1030.32996, 240, 1},
{1792.01001, 1030.32996, 240, 1},
{1792.01001, 1155.53003, 240, 1},
};
const vec4f_t tears_quad2[] = {
{304, 3248, 3208, 1},
{-61.2307014, 3248, 3208, 1},
{-61.2307014, 3248, 3192, 1},
{304, 3248, 3192, 1},
};
const vec4f_t tears_cluster[] = {
{551.638, -439.277008, 2804.63989, 1},
{551, -438, 2806.02002, 1},
{555.999023, -447.998993, 2805.1499, 1},
{555.999023, -447.998993, 2804.32007, 1},
{551.638, -439.277008, 2804.63989, 1},
{555.999023, -447.998993, 2804.32007, 1},
{555.999023, -447.998993, 2795.19995, 1},
{550.97699, -437.954987, 2806.02002, 1},
{551, -438, 2806.02002, 1},
{555.999023, -447.998993, 2795.19995, 1},
{555.999023, -447.998993, 2794.6499, 1},
{550.978027, -437.954987, 2806.02002, 1},
{545.945984, -442.945007, 2804.78003, 1},
{546, -442.998993, 2804.77002, 1},
{546, -443, 2804.69995, 1},
{546, -443, 2804.78003, 1},
{545.945007, -442.945007, 2804.78003, 1},
{556, -448, 2794.6499, 1},
{556, -448, 2805.1499, 1},
{546, -443, 2804.78003, 1},
{546, -443, 2804.69995, 1},
{550.97699, -437.954987, 2806.02002, 1},
{556, -448, 2794.6499, 1},
{546.000977, -443, 2804.69995, 1},
{545.945007, -442.945007, 2804.78003, 1},
};
const vec4f_t tears_cluster2[] = {
{-2160, -192, -704.242004, 1},
{-2160, -192, -662, 1},
{-2160, -196, -658, 1},
{-2160, -220, -646, 1},
{-2160, -222.332001, -645.416992, 1},
{-2160, -320, -752, 1},
{-2160, -192, -752, 1},
{-2160, -192, -704.242004, 1},
{-2160, -222.332001, -645.416992, 1},
{-2160, -244, -640, 1},
{-2160, -268, -640, 1},
{-2160, -292, -646, 1},
{-2160, -316, -658, 1},
{-2160, -320, -662, 1},
{-2160, -196, -814, 1},
{-2160, -192, -810, 1},
{-2160, -192, -752, 1},
{-2160, -320, -752, 1},
{-2160, -320, -810, 1},
{-2160, -316, -814, 1},
{-2160, -312, -816, 1},
{-2160, -200, -816, 1},
{-2174, -334, -676, 1},
{-2186, -346, -700, 1},
{-2192, -352, -724, 1},
{-2192, -352, -748, 1},
{-2186, -346, -772, 1},
{-2174, -334, -796, 1},
{-2160, -319.998993, -810, 1},
{-2160, -319.998993, -661.999023, 1},
{-2174, -178, -796, 1},
{-2186, -166, -772, 1},
{-2192, -160, -748, 1},
{-2192, -160, -724, 1},
{-2186, -166, -700, 1},
{-2174, -178, -676, 1},
{-2160, -192, -662, 1},
{-2160, -192, -810, 1},
{-2224, -320, -810, 1},
{-2224, -330, -800, 1},
{-2224, -320, -800, 1},
{-2224, -200, -816, 1},
{-2224, -312, -816, 1},
{-2224, -316, -814, 1},
{-2224, -320, -810, 1},
{-2224, -320, -800, 1},
{-2224, -192, -800, 1},
{-2224, -192, -810, 1},
{-2224, -196, -814, 1},
{-2224, -192, -800, 1},
{-2224, -182, -800, 1},
{-2224, -192, -810, 1},
{-2224, -330, -800, 1},
{-2224, -334, -796, 1},
{-2224, -178, -796, 1},
{-2224, -182, -800, 1},
{-2224, -334, -796, 1},
{-2224, -346, -772, 1},
{-2224, -351, -752, 1},
{-2224, -161, -752, 1},
{-2224, -166, -772, 1},
{-2224, -178, -796, 1},
{-2224, -351, -752, 1},
{-2224, -352, -748, 1},
{-2224, -352, -724, 1},
{-2224, -351, -720, 1},
{-2224, -161, -720, 1},
{-2224, -160, -724, 1},
{-2224, -160, -748, 1},
{-2224, -161, -752, 1},
{-2224, -351, -720, 1},
{-2224, -346, -700, 1},
{-2224, -334, -676, 1},
{-2224, -316, -658, 1},
{-2224, -292, -646, 1},
{-2224, -268, -640, 1},
{-2224, -244, -640, 1},
{-2224, -220, -646, 1},
{-2224, -196, -658, 1},
{-2224, -178, -676, 1},
{-2224, -166, -700, 1},
{-2224, -161, -720, 1},
};
// random numbers dug this up. needed more than 6 iterations
const vec4f_t cloud_points[] = {
{1.70695281, 2.60889673, 2.86233163, 1},
{2.11825109, 2.59883952, 1.48193693, 1},
{2.76728868, 1.57646465, 2.39769602, 1},
{1.88795376, 2.73980999, 1.58031297, 1},
{2.79736757, 1.89582968, 1.80514407, 1},
};
struct {
const vec4f_t *points;
int num_points;
vspheref_t expect;
} tests[] = {
#if 0
{0, 0, {{ 0, 0, 0, 1}, 0}},
{points, 1, {{-1, -1, 1, 1}, 0}},
{points, 2, {{ 0, 0, 1, 1}, 1.41421356}},
{points, 3, {{-0.333333343, 0.333333343, 0.333333343, 1}, 1.63299322}},
{points, 4, {{0, 0, 0, 1}, 1.73205081}},
{tears_triangle, 3, {{2201.84521, -1262, -713.519531}, 0.0747000724}},
{tears_triangle, SIZEOF (tears_triangle),
{{2201.84521, -1262, -713.519531}, 0.0747000724}},//FIXME numbers?
{tears_cluster, SIZEOF (tears_cluster),
{{552.678101, -443.460876, 2800.40405}, 8.04672813}},//FIXME numbers?
{tears_quad, SIZEOF (tears_quad),
{{1792.005, 1092.92999, 240}, 62.6000302}},
{tears_triangle2, SIZEOF (tears_triangle2),
{{-519.995972, 1111.97498, -2592.03516}, 0.0300767105}},//FIXME numbers?
{tears_quad2, SIZEOF (tears_quad2),
{{121.384659, 3248, 3200}, 182.790497}},
{tears_triangle3, SIZEOF (tears_triangle3),
{{284.117493, 624.235535, 2959.76489}, 0.352925777}},
{tears_cluster2, SIZEOF (tears_cluster2),
{{-2192, -256.000031, -736}, 103.479492}},//FIXME numbers?
#endif
{cloud_points, SIZEOF (cloud_points),
{{2.20056152, 2.23369908, 2.2332375}, 0.88327992}},
};
#define num_tests (sizeof (tests) / sizeof (tests[0]))
#define num_tests SIZEOF (tests)
static inline float
rnd (mtstate_t *mt)
@ -92,6 +257,7 @@ main (int argc, const char **argv)
for (j = 0; j < 5; j++) {
VectorSet (rnd (&mt), rnd (&mt), rnd (&mt), cloud[j]);
cloud[j][3] = 1;
}
seb = SmallestEnclosingBall_vf (cloud, 5);
r2 = seb.radius * seb.radius;