From eba5010c4d8fb7708b2859c7af0a6019560c30af Mon Sep 17 00:00:00 2001 From: Martin Uecker Date: Sun, 4 Oct 2020 10:11:19 +0200 Subject: [PATCH] update geometry lib --- Makefile | 4 +++ src/geom/draw.c | 7 ++++ src/geom/logo.c | 2 +- src/geom/polygon.c | 4 +-- src/geom/polyhedron.c | 53 +++++++++++++++++++++++++++ src/geom/polyhedron.h | 4 +++ src/geom/triangle.c | 70 ++++++++++++++++++++++++++++++++++++ src/geom/triangle.h | 7 ++++ utests/test_geom.c | 83 +++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 231 insertions(+), 3 deletions(-) create mode 100644 src/geom/polyhedron.c create mode 100644 src/geom/polyhedron.h create mode 100644 src/geom/triangle.c create mode 100644 src/geom/triangle.h create mode 100644 utests/test_geom.c diff --git a/Makefile b/Makefile index b53fb953d..de200fed3 100644 --- a/Makefile +++ b/Makefile @@ -523,6 +523,10 @@ UTARGETS += test_ode_bloch test_biot_savart MODULES_test_ode_bloch += -lsimu MODULES_test_biot_savart += -lsimu +# lib geom +UTARGETS += test_geom +MODULES_test_geom += -lgeom + # lib iter UTARGETS += test_iter test_prox MODULES_test_iter += -liter -lnlops -llinops diff --git a/src/geom/draw.c b/src/geom/draw.c index ad35c96f1..99b6b4468 100644 --- a/src/geom/draw.c +++ b/src/geom/draw.c @@ -1,3 +1,10 @@ +/* Copyright 2019-2020. Martin Uecker. + * All rights reserved. Use of this source code is governed by + * a BSD-style license which can be found in the LICENSE file. + * + * Authors: + * 2019-2020 Martin Uecker + */ #include #include diff --git a/src/geom/logo.c b/src/geom/logo.c index d0adac92f..fe1438994 100644 --- a/src/geom/logo.c +++ b/src/geom/logo.c @@ -1,4 +1,4 @@ - + #include "logo.h" const double bart_logo[77][2][4] = { // cubic hermite spline (really just floats) diff --git a/src/geom/polygon.c b/src/geom/polygon.c index f3382374c..2a2672f8f 100644 --- a/src/geom/polygon.c +++ b/src/geom/polygon.c @@ -2,10 +2,10 @@ * All rights reserved. Use of this source code is governed by * a BSD-style license which can be found in the LICENSE file. * - * 2019-2020 Martin Uecker + * Authors: + * 2019-2020 Martin Uecker */ - #include "polygon.h" diff --git a/src/geom/polyhedron.c b/src/geom/polyhedron.c new file mode 100644 index 000000000..a0d128193 --- /dev/null +++ b/src/geom/polyhedron.c @@ -0,0 +1,53 @@ +/* Copyright 2020. Martin Uecker. + * All rights reserved. Use of this source code is governed by + * a BSD-style license which can be found in the LICENSE file. + * + * Authors: + * 2019-2020 Martin Uecker + */ + +#include "polyhedron.h" + +static void vec3d_sub(double n[3], const double a[3], const double b[3]) +{ + for (int j = 0; j < 3; j++) + n[j] = a[j] - b[j]; +} + +static void vec3d_rot(double n[3], const double a[3], const double b[3]) +{ + n[0] = a[1] * b[2] - a[2] * b[1]; + n[1] = a[2] * b[0] - a[0] * b[2]; + n[2] = a[0] * b[1] - a[1] * b[0]; +} + +static double vec3d_dot(const double a[3], const double b[3]) +{ + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + + +double polyhedron_vol(int N, const double tri[N][3][3]) +{ + double sum = 0.; + + for (int i = 0; i < N; i++) { + + double a[3]; + double b[3]; + + vec3d_sub(a, tri[i][0], tri[i][1]); + vec3d_sub(b, tri[i][0], tri[i][2]); + + double n[3]; + vec3d_rot(n, a, b); + double r = vec3d_dot(n, tri[i][0]); + + sum += r / 6.; + } + + return sum; +} + + + diff --git a/src/geom/polyhedron.h b/src/geom/polyhedron.h new file mode 100644 index 000000000..65f7ff2e7 --- /dev/null +++ b/src/geom/polyhedron.h @@ -0,0 +1,4 @@ + +extern double polyhedron_vol(int N, const double tri[N][3][3]); + + diff --git a/src/geom/triangle.c b/src/geom/triangle.c new file mode 100644 index 000000000..e2bba46bb --- /dev/null +++ b/src/geom/triangle.c @@ -0,0 +1,70 @@ +/* Copyright 2019-2020. Martin Uecker. + * All rights reserved. Use of this source code is governed by + * a BSD-style license which can be found in the LICENSE file. + * + * Authors: + * 2019-2020 Martin Uecker + */ + +#include +#include + +#include "num/vec3.h" + +#include "triangle.h" + + +// Moeller-Trumbore + +float (triangle_intersect)(float uv[2], const vec3_t o, const vec3_t d, const vec3_t tri[3]) +{ + vec3_t e1, e2; + vec3_sub(e1, tri[1], tri[0]); + vec3_sub(e2, tri[2], tri[0]); + + vec3_t pvec; + vec3_rot(pvec, d, e2); + float det = vec3_sdot(e1, pvec); + + if (fabsf(det) < 1.E-8) // parallel + return 0.; + + float idet = 1. / det; + + vec3_t tvec; + vec3_sub(tvec, o, tri[0]); + float u = vec3_sdot(tvec, pvec) * idet; + + if ((u < 0.) || (u > 1.)) + return 0.; + + vec3_t qvec; + vec3_rot(qvec, tvec, e1); + float v = vec3_sdot(d, qvec) * idet; + + if ((v < 0.) || (u + v > 1.)) + return 0.; + + uv[0] = u; + uv[1] = v; + + return vec3_sdot(e2, qvec) * idet; +} + + +static float det2d(const float a[2], const float b[2]) +{ + return a[0] * b[1] - a[1] * b[0]; +} + +bool (triangle2d_inside)(const float tri[3][2], const float p[2]) +{ + float ab[2] = { tri[1][0] - tri[0][0], tri[1][1] - tri[0][1] }; + float ac[2] = { tri[2][0] - tri[0][0], tri[2][1] - tri[0][1] }; + float d = det2d(ab, ac); + float u = (det2d(p, ac) - det2d(tri[0], ac)) / d; + float v = (det2d(p, ab) - det2d(tri[0], ab)) / d; + return (0. < u) && (0. < v) && (u + v < 1.); +} + + diff --git a/src/geom/triangle.h b/src/geom/triangle.h new file mode 100644 index 000000000..d4127ca65 --- /dev/null +++ b/src/geom/triangle.h @@ -0,0 +1,7 @@ + + +typedef float vec3_t[3]; + +extern float triangle_intersect(float uv[2], const vec3_t o, const vec3_t d, const vec3_t tri[3]); +extern _Bool triangle2d_inside(const float tri[3][2], const float p[2]); + diff --git a/utests/test_geom.c b/utests/test_geom.c new file mode 100644 index 000000000..562e794fd --- /dev/null +++ b/utests/test_geom.c @@ -0,0 +1,83 @@ +/* Copyright 2020. Martin Uecker. + * All rights reserved. Use of this source code is governed by + * a BSD-style license which can be found in the LICENSE file. + * + * Authors: + * 2019-2020 Martin Uecker + */ + +#include + +#include "geom/polygon.h" +#include "geom/polyhedron.h" +#include "geom/triangle.h" + +#include "utest.h" + + +static bool test_winding_number(void) +{ + double p[2] = { 0., 0. }; + + double pg[4][2] = { + { 1., -1. }, + { 1., 1. }, + { -1., 1. }, + { -1., -1. }, + }; + + return (1 == polygon_winding_number(4, pg, p)); +} + +UT_REGISTER_TEST(test_winding_number); + + +static bool test_triangle_intersect(void) +{ + float o[3] = { 0., 0., 0. }; + float d[3] = { 1., 0., 0. }; + float t[3][3] = { + + { 1., -1., -1. }, + { 1., -1., 1. }, + { 1., 1., 0. }, + }; + + float uv[2]; + + return (1. == triangle_intersect(uv, o, d, t)); +} + + +UT_REGISTER_TEST(test_triangle_intersect); + + + +static bool test_polygon_area(void) +{ + double pg[4][2] = { { 0., 0. }, { 1., 0. }, { 1., 1. }, { 0., 1. } }; + + return (1. == polygon_area(4, pg)); +} + +UT_REGISTER_TEST(test_polygon_area); + + + + +static bool test_polyhedron_volume(void) +{ + double ph[4][3][3] = { + { { +0.5, +0.5, +0.5 }, { +0.5, -0.5, +0.5 }, { +0.5, +0.5, -0.5 } }, + { { +0.5, +0.5, +0.5 }, { +0.5, +0.5, -0.5 }, { -0.5, +0.5, +0.5 } }, + { { +0.5, -0.5, +0.5 }, { +0.5, +0.5, +0.5 }, { -0.5, +0.5, +0.5 } }, + { { +0.5, -0.5, +0.5 }, { -0.5, +0.5, +0.5 }, { +0.5, +0.5, -0.5 } }, + }; + + return fabs(1. / 6. - polyhedron_vol(4, ph)) < 1.E-16; +} + +UT_REGISTER_TEST(test_polyhedron_volume); + + +