update geometry lib

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

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 <[email protected]>
+ */
 
 #include <stdbool.h>
 #include <complex.h>

src/geom/logo.c

@@ -1,4 +1,4 @@
-
+ 
 #include "logo.h"
 
 const double bart_logo[77][2][4] = { // cubic hermite spline (really just floats)

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 <[email protected]>
  */
 
-
 #include "polygon.h"
 
 

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 <[email protected]>
+ */
+
+#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;
+}
+
+
+

src/geom/polyhedron.h

@@ -0,0 +1,4 @@
+
+extern double polyhedron_vol(int N, const double tri[N][3][3]);
+
+

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 <[email protected]>
+ */
+
+#include <math.h>
+#include <stdbool.h>
+
+#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.);
+}
+
+

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]);
+

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 <[email protected]>
+ */
+
+#include <math.h>
+
+#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);
+
+
+