Add function for working with triangle circumcircles

diff --git a/refine/triangle.c b/refine/triangle.c
index ae5ae02a33cc654f97fd3d23d8067131005d9a24..cd2df894392d0596ec4acf2cba4b96b615ca8e5e 100644 (file)
--- a/refine/triangle.c
+++ b/refine/triangle.c
@@ -1,3 +1,4 @@
+#include <math.h>
 #include <glib.h>
 #include "point.h"
 #include "edge.h"
@@ -158,6 +159,9 @@ p2tr_triangle_get_opposite_edge (P2trTriangle *self,
   p2tr_exception_programmatic ("The point is not in the triangle!");
 }
 
+/**
+ * Angles return by this function are always in the range [0,180]
+ */
 gdouble
 p2tr_triangle_get_angle_at (P2trTriangle *self,
                             P2trPoint    *p)
@@ -198,3 +202,55 @@ p2tr_triangle_smallest_non_constrained_angle (P2trTriangle *self)
 
     return result;
 }
+
+void
+p2tr_triangle_get_circum_circle (P2trTriangle *self,
+                                 P2trCircle   *circle)
+{
+  //       | Ax Bx Cx |
+  // D = + | Ay By Cy | * 2
+  //       | +1 +1 +1 |
+  //
+  //       | Asq Bsq Csq |
+  // X = + | Ay  By  Cy  | / D
+  //       | 1   1   1   |
+  //
+  //       | Asq Bsq Csq |
+  // Y = - | Ax  Bx  Cx  | / D
+  //       | 1   1   1   |
+  P2trPoint *A = self->edges[0]->end;
+  P2trPoint *B = self->edges[1]->end;
+  P2trPoint *C = self->edges[2]->end;
+
+  gdouble Asq = P2TR_VECTOR2_LEN_SQ (&A->c);
+  gdouble Bsq = P2TR_VECTOR2_LEN_SQ (&B->c);
+  gdouble Csq = P2TR_VECTOR2_LEN_SQ (&C->c);
+
+  gdouble invD = 1 / (2 * p2tr_matrix_det3 (
+      A->c.x, B->c.x, C->c.x,
+      A->c.y, B->c.y, C->c.y,
+      1,      1,      1));
+
+  circle->center.x = + p2tr_matrix_det3 (
+      Asq,    Bsq,    Csq,
+      A->c.y, B->c.y, C->c.y,
+      1,      1,      1) * invD;
+
+  circle->center.y = - p2tr_matrix_det3 (
+      Asq,    Bsq,    Csq,
+      A->c.x, B->c.x, C->c.x,
+      1,      1,      1) * invD;
+
+  circle->radius = sqrt (P2TR_VECTOR2_DISTANCE_SQ (&A->c, &circle->center));
+}
+
+P2trInCircle
+p2tr_triangle_circumcircle_contains_point (P2trTriangle *self,
+                                           P2trVector2  *pt)
+{
+  return p2tr_math_orient2d (
+      &self->edges[0]->end.c,
+      &self->edges[1]->end.c,
+      &self->edges[2]->end.c,
+      pt);
+}

diff --git a/refine/triangle.h b/refine/triangle.h
index 65ad46f6ba0a359d1534bc2f41213f1e403ca86c..d88fe6552d6bdcf30c54b926f9993304c800916c 100644 (file)
--- a/refine/triangle.h
+++ b/refine/triangle.h
@@ -40,6 +40,11 @@ P2trEdge*   p2tr_triangle_get_opposite_edge  (P2trTriangle *self,
 gdouble     p2tr_triangle_get_angle_at       (P2trTriangle *self,
                                               P2trPoint    *p);
 
-gdouble    p2tr_triangle_smallest_non_constrained_angle (P2trTriangle *self);
+gdouble     p2tr_triangle_smallest_non_constrained_angle (P2trTriangle *self);
 
+void        p2tr_triangle_get_circum_circle (P2trTriangle *self,
+                                             P2trCircle   *circle);
+
+P2trInCircle p2tr_triangle_circumcircle_contains_point (P2trTriangle *self,
+                                                        P2trVector2  *pt);
 #endif
\ No newline at end of file