Add triangle map and quadrilateral representation and common operations

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

2 files changed, 122 insertions, 0 deletions

diff --git a/src/morpher/quadrilateral.c b/src/morpher/quadrilateral.c
new file mode 100644
index 0000000..b9e740b
--- /dev/null
+++ b/src/morpher/quadrilateral.c
@@ -0,0 +1,58 @@
+#include "morpher/quadrilateral.h"
+#include "common/geom.h"
+#include "morpher/matrix.h"
+
+static inline IntVector p2(IntVector n) {
+  return n * n;
+}
+
+static inline bool in_circumcircle(TriangleMap *t, CartesianVector v) {
+  CartesianVector a = t->vertices[0].origin, b = t->vertices[1].origin, c = t->vertices[2].origin;
+  IntVector v2 = p2(v.x) + p2(v.y);
+  return matrix_int_det3(a.x - v.x, a.y - v.y, p2(a.x) + p2(a.y) - v2,
+                         b.x - v.x, b.y - v.y, p2(b.x) + p2(b.y) - v2,
+                         c.x - v.x, c.y - v.y, p2(c.x) + p2(c.y) - v2) > 0;
+}
+
+static inline void rotate_vertices(TriangleMap *t1, TriangleMap *t2, int e1, int e2) {
+  CartesianMapping quad[] = {
+   t1->vertices[(e1 + 1) % 3],
+   t1->vertices[(e1 + 2) % 3],
+   t2->vertices[(e2 + 1) % 3],
+   t2->vertices[(e2 + 2) % 3]
+  };
+
+  t1->vertices[(e1 + 1) % 3] = quad[1];
+  t1->vertices[(e1 + 2) % 3] = quad[2];
+  t1->vertices[(e1 + 3) % 3] = quad[3];
+  t2->vertices[(e2 + 1) % 3] = quad[3];
+  t2->vertices[(e2 + 2) % 3] = quad[0];
+  t2->vertices[(e2 + 3) % 3] = quad[1];
+}
+
+static inline void rotate_neighbors(TriangleMap *t1, TriangleMap *t2, int e1, int e2) {
+  TriangleMap *neighbors[] = {
+   t1->neighbors[(e1 + 1) % 3],
+   t1->neighbors[(e1 + 2) % 3],
+   t2->neighbors[(e2 + 1) % 3],
+   t2->neighbors[(e2 + 2) % 3]
+  };
+
+  t1->neighbors[(e1 + 1) % 3] = neighbors[1];
+  t1->neighbors[(e1 + 2) % 3] = neighbors[2];
+  t2->neighbors[(e2 + 1) % 3] = neighbors[3];
+  t2->neighbors[(e2 + 2) % 3] = neighbors[0];
+  trianglemap_replace_neighbor(t1->neighbors[(e1 + 2) % 3], t2, t1);
+  trianglemap_replace_neighbor(t2->neighbors[(e2 + 2) % 3], t1, t2);
+}
+
+void quadrilateral_flip_diagonal(TriangleMap *t1, TriangleMap *t2) {
+  int 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) {
+  return in_circumcircle(t1, t2->vertices[(trianglemap_find_common_edge(t2, t1) + 2) % 3].origin) &&
+         in_circumcircle(t2, t1->vertices[(trianglemap_find_common_edge(t1, t2) + 2) % 3].origin);
+}
diff --git a/src/morpher/trianglemap.c b/src/morpher/trianglemap.c
new file mode 100644
index 0000000..45f4ade
--- /dev/null
+++ b/src/morpher/trianglemap.c
@@ -0,0 +1,64 @@
+#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) {
+  TriangleMap *triangle = malloc_or_die(sizeof(TriangleMap));
+  triangle->vertices[0] = v1;
+  triangle->vertices[1] = v2;
+  triangle->vertices[2] = v3;
+  return triangle;
+}
+
+TriangleMap *trianglemap_destroy(TriangleMap *t) {
+  TriangleMap *next = t->next;
+  free(t);
+  return next;
+}
+
+TriangleMap *trianglemap_to(TriangleMap *t, CartesianVector v) {
+  int edge;
+
+  for (edge = 0; edge < 3; ++edge)
+    if (triangle_area(t->vertices[edge].origin, t->vertices[(edge + 1) % 3].origin, v) >= 0)
+      return t->neighbors[edge];
+
+  return t;
+}
+
+int trianglemap_find_common_edge(TriangleMap *t, TriangleMap *neighbor) {
+  int edge;
+  for (edge = 0; t->neighbors[edge] != neighbor && edge < 3; ++edge);
+  assert(t->neighbors[edge] == neighbor && neighbor != NULL);
+  return edge;
+}
+
+void trianglemap_set_neighbors(TriangleMap *t, TriangleMap *n1, TriangleMap *n2, TriangleMap *n3, TriangleMap *next) {
+  t->neighbors[0] = n1;
+  t->neighbors[1] = n2;
+  t->neighbors[2] = n3;
+  t->next = next;
+}
+
+void trianglemap_replace_neighbor(TriangleMap *t, TriangleMap *old, TriangleMap *new) {
+  if (t != NULL) t->neighbors[trianglemap_find_common_edge(t, old)] = new;
+}
+
+TriangleMap *trianglemap_split(TriangleMap *t, CartesianMapping v) {
+  assert(trianglemap_to(t, v.origin) == t);
+
+  TriangleMap *triangle2 = trianglemap_create(t->vertices[1], t->vertices[2], v);
+  TriangleMap *triangle3 = trianglemap_create(t->vertices[2], t->vertices[0], v);
+  t->vertices[2] = v;
+
+  trianglemap_replace_neighbor(t->neighbors[1], t, triangle2);
+  trianglemap_replace_neighbor(t->neighbors[2], t, triangle3);
+
+  trianglemap_set_neighbors(triangle3, t->neighbors[2], t, triangle2, t->next);
+  trianglemap_set_neighbors(triangle2, t->neighbors[1], triangle3, t, triangle3);
+  trianglemap_set_neighbors(t, t->neighbors[0], triangle2, triangle3, triangle2);
+
+  return t;
+}