Implement Delaunay propagation

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

2 files changed, 21 insertions, 3 deletions

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