Implement Delaunay propagation

Signed-off-by: pacien <pacien.trangirard@pacien.net>

5 files changed, 67 insertions, 4 deletions

diff --git a/include/morpher/quadrilateral.h b/include/morpher/quadrilateral.h
index c8cf3a1..8691263 100644
--- a/include/morpher/quadrilateral.h
+++ b/include/morpher/quadrilateral.h
@@ -36,4 +36,15 @@ void quadrilateral_flip_diagonal(TriangleMap *t1, TriangleMap *t2);
  */
 bool quadrilateral_is_delaunay(TriangleMap *t1, TriangleMap *t2);
 
+/**
+ * Function: quadrilateral_propagate_delaunay
+ * Ensures that the quadrilateral spawned by the given triangles fulfills the Delaunay criterion,
+ * flipping the diagonal if necessary and propagating the changes to the neighbouring triangles.
+ *
+ * Parameters:
+ *   *start    - the starting triangle
+ *   *neighbor - a neighboring triangle
+ */
+void quadrilateral_propagate_delaunay(TriangleMap *start, TriangleMap *neighbor);
+
 #endif
diff --git a/include/morpher/trianglemap.h b/include/morpher/trianglemap.h
index 05e7b87..0bb6a1d 100644
--- a/include/morpher/trianglemap.h
+++ b/include/morpher/trianglemap.h
@@ -106,6 +106,16 @@ void trianglemap_set_neighbors(TriangleMap *t, TriangleMap *n1, TriangleMap *n2,
 void trianglemap_replace_neighbor(TriangleMap *t, TriangleMap *old, TriangleMap *new);
 
 /**
+ * Function: trianglemap_foreach_neighbor
+ * Executes the given function for each existing (non-NULL) neighbour of the supplied triangle.
+ *
+ * Parameters:
+ *   *t - the base triangle
+ *   *f - the function
+ */
+void trianglemap_foreach_neighbor(TriangleMap *t, void (*f)(TriangleMap *current, TriangleMap *neighbor));
+
+/**
  * Function: trianglemap_split
  * Splits a triangle into three sub-triangles at the supplied center vertex, updating the surrounding triangles.
  * The first triangle resulting from the split is returned, with the two others chained as linear neighbours.
diff --git a/src/morpher/quadrilateral.c b/src/morpher/quadrilateral.c
index 9671116..af77f0f 100644
--- a/src/morpher/quadrilateral.c
+++ b/src/morpher/quadrilateral.c
@@ -1,5 +1,6 @@
 #include "morpher/quadrilateral.h"
-#include "common/geom.h"
+#include <stdlib.h>
+#include <assert.h>
 #include "morpher/matrix.h"
 
 static inline IntVector p2(IntVector n) {
@@ -47,12 +48,24 @@ static inline void rotate_neighbors(TriangleMap *t1, TriangleMap *t2, int e1, in
 }
 
 void quadrilateral_flip_diagonal(TriangleMap *t1, TriangleMap *t2) {
-  int e1 = trianglemap_find_common_edge(t1, t2), e2 = trianglemap_find_common_edge(t2, t1);
+  int e1, e2;
+  assert(t1 != NULL && t2 != NULL);
+  e1 = trianglemap_find_common_edge(t1, t2);
+  e2 = trianglemap_find_common_edge(t2, t1);
   rotate_vertices(t1, t2, e1, e2);
   rotate_neighbors(t1, t2, e1, e2);
 }
 
 bool quadrilateral_is_delaunay(TriangleMap *t1, TriangleMap *t2) {
+  assert(t1 != NULL && t2 != NULL);
   return not_in_circumcircle(t1, t2->vertices[(trianglemap_find_common_edge(t2, t1) + 2) % 3].origin) &&
          not_in_circumcircle(t2, t1->vertices[(trianglemap_find_common_edge(t1, t2) + 2) % 3].origin);
 }
+
+void quadrilateral_propagate_delaunay(TriangleMap *start, TriangleMap *neighbor) {
+  assert(start != NULL && neighbor != NULL);
+  if (!quadrilateral_is_delaunay(start, neighbor)) {
+    quadrilateral_flip_diagonal(start, neighbor);
+    trianglemap_foreach_neighbor(neighbor, quadrilateral_propagate_delaunay);
+  }
+}
diff --git a/src/morpher/trianglemap.c b/src/morpher/trianglemap.c
index 45f4ade..01b66c9 100644
--- a/src/morpher/trianglemap.c
+++ b/src/morpher/trianglemap.c
@@ -1,7 +1,6 @@
 #include "morpher/trianglemap.h"
 #include <stdlib.h>
 #include <assert.h>
-#include "common/geom.h"
 #include "common/mem.h"
 
 TriangleMap *trianglemap_create(CartesianMapping v1, CartesianMapping v2, CartesianMapping v3) {
@@ -46,6 +45,12 @@ void trianglemap_replace_neighbor(TriangleMap *t, TriangleMap *old, TriangleMap
   if (t != NULL) t->neighbors[trianglemap_find_common_edge(t, old)] = new;
 }
 
+void trianglemap_foreach_neighbor(TriangleMap *t, void (*f)(TriangleMap *, TriangleMap *)) {
+  int cursor;
+  assert(t != NULL);
+  for (cursor = 0; cursor < 3; ++cursor) if (t->neighbors[cursor] != NULL) f(t, t->neighbors[cursor]);
+}
+
 TriangleMap *trianglemap_split(TriangleMap *t, CartesianMapping v) {
   assert(trianglemap_to(t, v.origin) == t);
 
diff --git a/test/morpher/quadrilateral.c b/test/morpher/quadrilateral.c
index 963d1d6..b278e2f 100644
--- a/test/morpher/quadrilateral.c
+++ b/test/morpher/quadrilateral.c
@@ -21,7 +21,6 @@ static inline bool vertices_equals(CartesianMapping maps[],
  */
 static void test_quadrilateral_flip_diagonal() {
   CartesianMapping A = m(50, 0), B = m(0, 100), C = m(100, 100), D = m(0, 200), E = m(100, 200);
-
   TriangleMap *t = trianglemap_create(A, B, C);
   TriangleMap *r = trianglemap_create(C, B, E);
   TriangleMap *l = trianglemap_create(B, D, E);
@@ -54,9 +53,34 @@ static void test_quadrilateral_is_delaunay(CartesianMapping A, CartesianMapping
   trianglemap_destroy(r);
 }
 
+static void test_quadrilateral_propagate_delaunay() {
+  CartesianMapping A = m(0, 0), B = m(0, 10), C = m(10, 10), D = m(10, 0), E = m(4, 6), F = m(8, 7);
+  TriangleMap *l = trianglemap_create(A, B, C);
+  TriangleMap *r = trianglemap_create(A, C, D);
+  trianglemap_set_neighbors(l, NULL, NULL, r, r);
+  trianglemap_set_neighbors(r, l, NULL, NULL, NULL);
+  trianglemap_split(l, E);
+  trianglemap_split(r, F);
+
+  trianglemap_foreach_neighbor(r, quadrilateral_propagate_delaunay);
+
+  {
+    TriangleMap *triangle;
+    int neighbor_index;
+    for (triangle = l; triangle != NULL; triangle = triangle->next)
+      for (neighbor_index = 0; neighbor_index < 3; ++neighbor_index)
+        if (triangle->neighbors[neighbor_index] != NULL)
+          assert(quadrilateral_is_delaunay(triangle, triangle->neighbors[neighbor_index]));
+  }
+
+  trianglemap_destroy(l);
+  trianglemap_destroy(r);
+}
+
 int main(int argc, char **argv) {
   test_quadrilateral_flip_diagonal();
   test_quadrilateral_is_delaunay(m(0, 0), m(0, 3), m(3, 3), m(2, 1), false);
   test_quadrilateral_is_delaunay(m(0, 0), m(0, 3), m(3, 3), m(4, -1), true);
+  test_quadrilateral_propagate_delaunay();
   return 0;
 }