Add a function to get the barycentric coords of a point.

It "works" for lines, triangles and tetrahedrons. For lines and triangles,
it gives the barycentric coordinates of the perpendicular projection of the
point onto to features. Only tetrahedrons are guaranteed to reproduce the
original point.
This commit is contained in:
Bill Currie 2013-03-12 14:10:32 +09:00
parent fe55bb678e
commit 0cd6d93030
4 changed files with 174 additions and 2 deletions

View file

@ -213,6 +213,9 @@ R_CullSphere (const vec3_t origin, const float radius)
return false;
}
void BarycentricCoords (const vec_t **points, int num_points, const vec3_t p,
vec_t *lambda);
//@}
#endif // __mathlib_h

View file

@ -1151,3 +1151,52 @@ Mat4Decompose (const mat4_t mat, quat_t rot, vec3_t shear, vec3_t scale,
Mat4toMat3 (mat, m3);
return Mat3Decompose (m3, rot, shear, scale);
}
void
BarycentricCoords (const vec_t **points, int num_points, const vec3_t p,
vec_t *lambda)
{
vec3_t a, b, c, x, ab, bc, ca, n;
vec_t div;
if (num_points > 4)
Sys_Error ("Don't know how to compute the barycentric coordinates "
"for %d points", num_points);
switch (num_points) {
case 1:
lambda[0] = 1;
return;
case 2:
VectorSubtract (p, points[0], x);
VectorSubtract (points[1], points[0], a);
lambda[1] = DotProduct (x, a) / DotProduct (a, a);
lambda[0] = 1 - lambda[1];
return;
case 3:
VectorSubtract (p, points[0], x);
VectorSubtract (points[1], points[0], a);
VectorSubtract (points[2], points[0], b);
CrossProduct (a, b, ab);
div = DotProduct (ab, ab);
CrossProduct (x, b, n);
lambda[1] = DotProduct (n, ab) / div;
CrossProduct (a, x, n);
lambda[2] = DotProduct (n, ab) / div;
lambda[0] = 1 - lambda[1] - lambda[2];
return;
case 4:
VectorSubtract (p, points[0], x);
VectorSubtract (points[1], points[0], a);
VectorSubtract (points[2], points[0], b);
VectorSubtract (points[3], points[0], c);
CrossProduct (a, b, ab);
CrossProduct (b, c, bc);
CrossProduct (c, a, ca);
div = DotProduct (a, bc);
lambda[1] = DotProduct (x, bc) / div;
lambda[2] = DotProduct (x, ca) / div;
lambda[3] = DotProduct (x, ab) / div;
lambda[0] = 1 - lambda[1] - lambda[2] - lambda[3];
return;
}
Sys_Error ("Not enough points to project or enclose the point");
}

View file

@ -3,8 +3,12 @@ AUTOMAKE_OPTIONS= foreign
INCLUDES= -I$(top_srcdir)/include
check_PROGRAMS= \
test-dq test-half test-mat3 test-mat4 test-plist test-qfs test-quat \
test-set test-vrect
test-bary test-dq test-half test-mat3 test-mat4 test-plist test-qfs \
test-quat test-set test-vrect
test_bary_SOURCES=test-bary.c
test_bary_LDADD=$(top_builddir)/libs/util/libQFutil.la
test_bary_DEPENDENCIES=$(top_builddir)/libs/util/libQFutil.la
test_dq_SOURCES=test-dq.c
test_dq_LDADD=$(top_builddir)/libs/util/libQFutil.la

116
libs/util/test/test-bary.c Normal file
View file

@ -0,0 +1,116 @@
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include "QF/mathlib.h"
#include "QF/mersenne.h"
#include "QF/sys.h"
vec3_t points[] = {
{-1, -1, 1},
{ 1, 1, 1},
{-1, 1, -1},
{ 1, -1, -1},
{-1, -1, -1},
{ 1, 1, -1},
{-1, 1, 1},
{ 1, -1, 1},
{ 0, 0, 0},
};
const vec_t *line[] = {points[0], points[1]};
const vec_t *tri[] = {points[0], points[1], points[2]};
const vec_t *tetra[] = {points[0], points[1], points[2], points[3]};
struct {
const vec_t **points;
int num_points;
vec_t *x;
vec_t expect[4];
} tests[] = {
{line, 2, points[0], {1, 0}},
{line, 2, points[1], {0, 1}},
{line, 2, points[2], {0.5, 0.5}},
{line, 2, points[3], {0.5, 0.5}},
{line, 2, points[8], {0.5, 0.5}},
{tri, 3, points[0], {1, 0, 0}},
{tri, 3, points[1], {0, 1, 0}},
{tri, 3, points[2], {0, 0, 1}},
{tri, 3, points[3], {0.333333284, 0.333333333, 0.333333333}},//rounding :P
{tri, 3, points[8], {0.333333284, 0.333333333, 0.333333333}},//rounding :P
{tetra, 4, points[0], {1, 0, 0, 0}},
{tetra, 4, points[1], {0, 1, 0, 0}},
{tetra, 4, points[2], {0, 0, 1, 0}},
{tetra, 4, points[3], {0, 0, 0, 1}},
{tetra, 4, points[4], { 0.5, -0.5, 0.5, 0.5}},
{tetra, 4, points[5], {-0.5, 0.5, 0.5, 0.5}},
{tetra, 4, points[6], { 0.5, 0.5, 0.5, -0.5}},
{tetra, 4, points[7], { 0.5, 0.5, -0.5, 0.5}},
{tetra, 4, points[8], {0.25, 0.25, 0.25, 0.25}},
};
#define num_tests (sizeof (tests) / sizeof (tests[0]))
static inline float
rnd (mtstate_t *mt)
{
union {
uint32_t u;
float f;
} uf;
do {
uf.u = mtwist_rand (mt) & 0x3fffffff;
} while (!uf.u);
return uf.f - 1.0;
}
int
main (int argc, const char **argv)
{
int res = 0;
size_t i;
int j;
vec_t lambda[4];
mtstate_t mt;
double start, end;
for (i = 0; i < num_tests; i ++) {
BarycentricCoords (tests[i].points, tests[i].num_points, tests[i].x,
lambda);
for (j = 0; j < tests[i].num_points; j++) {
if (tests[i].expect[j] != lambda[j])
break;
}
if (j != tests[i].num_points) {
res = 1;
printf ("test %d failed\n", (int) i);
printf ("expect:");
for (j = 0; j < tests[i].num_points; j++)
printf (" %.9g", tests[i].expect[j]);
printf ("\ngot :");
for (j = 0; j < tests[i].num_points; j++)
printf (" %.9g", lambda[j]);
printf ("\n");
}
}
mtwist_seed (&mt, 0);
start = Sys_DoubleTime ();
for (i = 0; i < 1000000; i++) {
vec3_t p, x;
VectorSet (rnd (&mt), rnd (&mt), rnd (&mt), p);
BarycentricCoords (tetra, 4, p, lambda);
VectorZero (x);
for (j = 0; j < 4; j++)
VectorMultAdd (x, lambda[j], tetra[j], x);
if (VectorDistance_fast (x, p) > 1e-4) {
res = 1;
printf ("[%.9g %.9g %.9g] != [%.9g %.9g %.9g]: [%g %g %g %g]\n",
VectorExpand (x), VectorExpand (p), QuatExpand (lambda));
break;
}
}
end = Sys_DoubleTime ();
printf ("%d itterations in %gs: %g itters/second\n", (int) i, end - start,
i / (end - start));
return res;
}