Add triangle map and quadrilateral representation and common operations

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

6 files changed, 426 insertions, 0 deletions

diff --git a/include/morpher/quadrilateral.h b/include/morpher/quadrilateral.h
new file mode 100644
index 0000000..c8cf3a1
--- /dev/null
+++ b/include/morpher/quadrilateral.h
@@ -0,0 +1,39 @@
+#ifndef UPEM_MORPHING_QUADRILATERAL
+#define UPEM_MORPHING_QUADRILATERAL
+
+/**
+ * File: quadrilateral.h
+ * Operations on quadrilaterals formed by adjacent triangles pairs.
+ */
+
+#include <stdbool.h>
+#include "common/geom.h"
+#include "morpher/trianglemap.h"
+
+/**
+ * Function: quadrilateral_flip_diagonal
+ * Flips the diagonal of a quadrilateral formed by two triangles sharing a common edge,
+ * using a positive rotation-like transform inverting both references.
+ * This transforms keeps the positive orientation of the vertices.
+ * Pointers to surrounding triangles are updated accordingly.
+ *
+ * Parameters:
+ *   *t1 - the first triangle
+ *   *t2 - the second triangle
+ */
+void quadrilateral_flip_diagonal(TriangleMap *t1, TriangleMap *t2);
+
+/**
+ * Function: quadrilateral_is_delaunay
+ * Tells whether the quadrilateral formed by the two supplied adjacent triangle forms a Delaunay triangulation.
+ *
+ * Parameters:
+ *   *t1 - first triangle
+ *   *t2 - second triangle adjacent to the first one
+ *
+ * Returns:
+ *   T(t1 and t2 are a Delaunay triangulation in the quadrilateral formed by the twos)
+ */
+bool quadrilateral_is_delaunay(TriangleMap *t1, TriangleMap *t2);
+
+#endif
diff --git a/include/morpher/trianglemap.h b/include/morpher/trianglemap.h
new file mode 100644
index 0000000..05e7b87
--- /dev/null
+++ b/include/morpher/trianglemap.h
@@ -0,0 +1,123 @@
+#ifndef UPEM_MORPHING_TRIANGLEMAP
+#define UPEM_MORPHING_TRIANGLEMAP
+
+/**
+ * File: trianglemap.h
+ * Represents a triangle map in a triangulation, in a plane oriented left to right and top to bottom.
+ */
+
+#include <stdbool.h>
+#include "common/geom.h"
+
+/**
+ * Struct: TriangleMap
+ * Represents a triangle mapping.
+ *
+ * Fields:
+ *   vertices[]   - array of vertices
+ *   *neighbors[] - array of neighboring triangles ordered by the edges spawned by the vertices
+ *   *next        - pointer to another triangle for linear browsing, or NULL
+ */
+typedef struct _TriangleMap {
+  CartesianMapping vertices[3];
+  struct _TriangleMap *neighbors[3];
+  struct _TriangleMap *next;
+} TriangleMap;
+
+/**
+ * Function: trianglemap_create
+ * Creates a TriangleMap, instantiating it on the heap.
+ *
+ * Parameters:
+ *   vertex1 - first vertex
+ *   vertex2 - second vertex
+ *   vertex3 - third vertex
+ *
+ * Returns:
+ *   A pointer to the newly created triangle
+ */
+TriangleMap *trianglemap_create(CartesianMapping v1, CartesianMapping v2, CartesianMapping v3);
+
+/**
+ * Function: trianglemap_destroy
+ * Destroys a triangle, freeing its dynamically allocated resources.
+ * Also returns the next triangle for convenience.
+ * Does not update surrounding triangles.
+ *
+ * Parameters:
+ *   *t - pointer to the triangle to destroy
+ *
+ * Returns:
+ *   A pointer to the next linear triangle
+ */
+TriangleMap *trianglemap_destroy(TriangleMap *t);
+
+/**
+ * Function: trianglemap_to
+ * Returns a pointer to the current or the closest adjacent neighbour TriangleMap
+ * minimizing the distance to the targeted position.
+ *
+ * Parameters:
+ *   *t - the origin triangle
+ *   v  - the target position vector
+ *
+ * Returns:
+ *   A pointer to the current or an immediate neighbour TriangleMap
+ */
+TriangleMap *trianglemap_to(TriangleMap *t, CartesianVector v);
+
+/**
+ * Function: trianglemap_find_common_edge
+ * Finds the index of the commonly shared edge in the neighbourhood of t.
+ *
+ * Parameters:
+ *   *t        - the base triangle
+ *   *neighbor - the neighbour to search for
+ *
+ * Returns:
+ *   The index of the common neighbour in the listing of t.
+ */
+int trianglemap_find_common_edge(TriangleMap *t, TriangleMap *neighbor);
+
+/**
+ * Function: trianglemap_set_neighbors
+ * Updates a triangle neighbors.
+ * Vertices must be given in a positively oriented (trigonometric) order.
+ * Neighbours must be given in the same order as the vertices.
+ *
+ * Parameters:
+ *   *t    - the triangle to modify
+ *   *n1   - first neighbour
+ *   *n2   - second neighbour
+ *   *n3   - third neighbour
+ *   *next - linear neighbour
+ */
+void trianglemap_set_neighbors(TriangleMap *t, TriangleMap *n1, TriangleMap *n2, TriangleMap *n3, TriangleMap *next);
+
+/**
+ * Function: triangle_replace_neighbor
+ * Substitutes a given neighbour in a neighbourhood.
+ *
+ * Parameters:
+ *   *t   - the base triangle
+ *   *old - the neighbour to replace
+ *   *new - the new neighbour
+ */
+void trianglemap_replace_neighbor(TriangleMap *t, TriangleMap *old, TriangleMap *new);
+
+/**
+ * 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.
+ * Those generated triangles each have one of the original surrounding neighbour as first element in their listings.
+ *
+ * Parameters:
+ *   *t - the triangle to split
+ *   v  - the new vertex to add
+ *
+ * Returns:
+ *   The first generated TriangleMap
+ */
+TriangleMap *trianglemap_split(TriangleMap *t, CartesianMapping v);
+
+#endif
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));