Add more necessary math functions

diff --git a/refine/math.c b/refine/math.c
new file mode 100644 (file)
index 0000000..1c26159
--- /dev/null
+++ b/refine/math.c
@@ -0,0 +1,206 @@
+#include <glib.h>\r
+#include "math.h"\r
+\r
+#define EDGE_LEN_SQ(X,Y)  ((X) * (X) + (Y) * (Y))\r
+#define VECTOR2_LEN_SQ(V) (EDGE_LEN_SQ((V)->x, (V)->y))\r
+#define VECTOR2_DOT(V1,V2) ((V1)->x * (V2)->x + (V1)->y * (V2)->y)\r
+\r
+gdouble\r
+p2tr_math_length_sq (gdouble x1, gdouble y1,\r
+                     gdouble x2, gdouble y2)\r
+{\r
+  return EDGE_LEN_SQ (x1 - x2, y1 - y2);\r
+}\r
+\r
+gdouble\r
+p2tr_math_length_sq2 (const P2trVector2 *pt1,\r
+                      const P2trVector2 *pt2)\r
+{\r
+  return p2tr_math_length_sq (pt1->x, pt1->y, pt2->x, pt2->y);\r
+}\r
+\r
+static inline gdouble\r
+p2tr_matrix_det2 (gdouble a00, gdouble a01,\r
+                  gdouble a10, gdouble a11)\r
+{\r
+    return a00 * a11 - a10 * a01;\r
+}\r
+\r
+static inline gdouble\r
+p2tr_matrix_det3 (gdouble a00, gdouble a01, gdouble a02,\r
+                  gdouble a10, gdouble a11, gdouble a12,\r
+                  gdouble a20, gdouble a21, gdouble a22)\r
+{\r
+    return\r
+        + a00 * (a11 * a22 - a21 * a12)\r
+        - a01 * (a10 * a22 - a20 * a12)\r
+        + a02 * (a10 * a21 - a20 * a11);\r
+}\r
+\r
+static inline gdouble\r
+p2tr_matrix_det4 (gdouble a00, gdouble a01, gdouble a02, gdouble a03,\r
+                  gdouble a10, gdouble a11, gdouble a12, gdouble a13,\r
+                  gdouble a20, gdouble a21, gdouble a22, gdouble a23,\r
+                  gdouble a30, gdouble a31, gdouble a32, gdouble a33)\r
+{\r
+    return\r
+        + a00 * p2tr_matrix_det3 (a11, a12, a13,\r
+                                 a21, a22, a23,\r
+                                 a31, a32, a33)\r
+        - a01 * p2tr_matrix_det3 (a10, a12, a13,\r
+                                  a20, a22, a23,\r
+                                  a30, a32, a33)\r
+        + a02 * p2tr_matrix_det3 (a10, a11, a13,\r
+                                  a20, a21, a23,\r
+                                  a30, a31, a33)\r
+        - a03 * p2tr_matrix_det3 (a10, a11, a12,\r
+                                  a20, a21, a22,\r
+                                  a30, a31, a32);\r
+}\r
+\r
+/* The point in triangle test which is implemented below is based on the\r
+ * algorithm which appears on:\r
+ *\r
+ *    http://www.blackpawn.com/texts/pointinpoly/default.html\r
+ */\r
+#define INTRIANGLE_EPSILON 0e-9\r
+\r
+P2trInTriangle\r
+p2tr_math_intriangle(const P2trVector2 *A,\r
+                     const P2trVector2 *B,\r
+                     const P2trVector2 *C,\r
+                     const P2trVector2 *P)\r
+{\r
+  /* Compute the vectors offsetted so that A is the origin */\r
+  P2trVector2 v0, v1, v2;\r
+  p2tr_vector2_sub(C, A, &v0);\r
+  p2tr_vector2_sub(B, A, &v1);\r
+  p2tr_vector2_sub(P, A, &v2);\r
+\r
+  /* Compute dot products */\r
+  double dot00 = VECTOR2_DOT(&v0, &v0);\r
+  double dot01 = VECTOR2_DOT(&v0, &v1);\r
+  double dot02 = VECTOR2_DOT(&v0, &v2);\r
+  double dot11 = VECTOR2_DOT(&v1, &v1);\r
+  double dot12 = VECTOR2_DOT(&v1, &v2);\r
+\r
+  /* Compute barycentric coordinates */\r
+  double invDenom = 1 / (dot00 * dot11 - dot01 * dot01);\r
+  double u = (dot11 * dot02 - dot01 * dot12) * invDenom;\r
+  double v = (dot00 * dot12 - dot01 * dot02) * invDenom;\r
+\r
+  /* Check if point is in triangle - i.e. whether (u + v) < 1 and both\r
+   * u and v are positive */\r
+  if ((u >= INTRIANGLE_EPSILON) && (v >= INTRIANGLE_EPSILON) && (u + v < 1 - INTRIANGLE_EPSILON))\r
+    return P2TR_INTRIANGLE_IN;\r
+  else if ((u >= -INTRIANGLE_EPSILON) && (v >= -INTRIANGLE_EPSILON) && (u + v <= 1 + INTRIANGLE_EPSILON))\r
+    return P2TR_INTRIANGLE_ON;\r
+  else\r
+    return P2TR_INTRIANGLE_OUT;\r
+}\r
+\r
+/* The point in triangle circumcircle test, and the 3-point orientation\r
+ * test, are both based on the work of Jonathan Richard Shewchuk. The\r
+ * technique used here is described in his paper "Adaptive Precision\r
+ * Floating-Point Arithmetic and Fast Robust Geometric Predicates"\r
+ */\r
+\r
+#define ORIENT2D_EPSILON 1e-9\r
+\r
+P2trOrientation p2tr_math_orient2d (const P2trVector2 *A,\r
+                                    const P2trVector2 *B,\r
+                                    const P2trVector2 *C)\r
+{\r
+  /* We are trying to compute this determinant:\r
+   * |Ax Ay 1|\r
+   * |Bx By 1|\r
+   * |Cx Cy 1|\r
+   */\r
+  gdouble result = p2tr_matrix_det3 (\r
+      A->x, A->y, 1,\r
+      B->x, B->y, 1,\r
+      C->x, C->y, 1\r
+  );\r
+\r
+  if (result > ORIENT2D_EPSILON)\r
+    return P2TR_ORIENTATION_CCW;\r
+  else if (result < ORIENT2D_EPSILON)\r
+    return P2TR_ORIENTATION_CW;\r
+  else\r
+    return P2TR_ORIENTATION_LINEAR;\r
+}\r
+\r
+#define INCIRCLE_EPSILON 1e-9\r
+\r
+P2trInCircle\r
+p2tr_math_incircle (const P2trVector2 *A,\r
+                    const P2trVector2 *B,\r
+                    const P2trVector2 *C,\r
+                    const P2trVector2 *D)\r
+{\r
+  /* We are trying to compute this determinant:\r
+   * |Ax Ay Ax^2+Ay^2 1|\r
+   * |Bx By Bx^2+By^2 1|\r
+   * |Cx Cy Cx^2+Cy^2 1|\r
+   * |Dx Dy Dx^2+Dy^2 1|\r
+   */\r
+  gdouble result = p2tr_matrix_det4 (\r
+      A->x, A->y, VECTOR2_LEN_SQ(A), 1,\r
+      B->x, B->y, VECTOR2_LEN_SQ(B), 1,\r
+      C->x, C->y, VECTOR2_LEN_SQ(C), 1,\r
+      D->x, D->y, VECTOR2_LEN_SQ(D), 1\r
+  );\r
+\r
+  if (result > INCIRCLE_EPSILON)\r
+    return P2TR_INCIRCLE_IN;\r
+  else if (result < INCIRCLE_EPSILON)\r
+    return P2TR_INCIRCLE_OUT;\r
+  else\r
+    return P2TR_INCIRCLE_ON;\r
+}\r
+\r
+/* The point inside diametral-circle test and the point inside diametral\r
+ * lens test, are both based on the work of Jonathan Richard Shewchuk.\r
+ * The techniques used here are partially described in his paper\r
+ * "Delaunay Refinement Algorithms for Triangular Mesh Generation".\r
+ *\r
+ * W is inside the diametral circle (lens) of the line XY if and only if\r
+ * the angle XWY is larger than 90 (120) degrees. We know how to compute\r
+ * the cosine of that angle very easily like so:\r
+ *\r
+ * cos XWY = (WX * WY) / (|WX| * |WY|)\r
+ *\r
+ * Since XWY is larger than 90 (120) degrees, cos XWY <= 0 (-0.5) so:\r
+ *\r
+ * Diametral Circle               | Diametral Lens\r
+ * -------------------------------+-----------------------------------\r
+ * 0 >= (WX * WY) / (|WX| * |WY|) | - 0.5 >= (WX * WY) / (|WX| * |WY|)\r
+ * 0 >= WX * WY                   | - 0.5 * |WX| * |WY| >= WX * WY\r
+ */\r
+\r
+gboolean\r
+p2tr_math_diametral_circle_contains (const P2trVector2 *X,\r
+                                     const P2trVector2 *Y,\r
+                                     const P2trVector2 *W)\r
+{\r
+  P2trVector2 WX, WY;\r
+\r
+  p2tr_vector2_sub(X, W, &WX);\r
+  p2tr_vector2_sub(Y, W, &WY);\r
+\r
+  return VECTOR2_DOT(&WX, &WY) <= 0;\r
+}\r
+\r
+gboolean\r
+p2tr_math_diametral_lens_contains (const P2trVector2 *X,\r
+                                   const P2trVector2 *Y,\r
+                                   const P2trVector2 *W)\r
+{\r
+  P2trVector2 WX, WY;\r
+\r
+  p2tr_vector2_sub(X, W, &WX);\r
+  p2tr_vector2_sub(Y, W, &WY);\r
+\r
+  return VECTOR2_DOT(&WX, &WY)\r
+      <= 0.5 * p2tr_vector2_norm(&WX) * p2tr_vector2_norm(&WY);\r
+}\r

diff --git a/refine/math.h b/refine/math.h
new file mode 100644 (file)
index 0000000..266aaa9
--- /dev/null
+++ b/refine/math.h
@@ -0,0 +1,50 @@
+#ifndef __P2TC_REFINE_MATH_H__\r
+#define __P2TC_REFINE_MATH_H__\r
+\r
+#include <glib.h>\r
+#include "vector2.h"
+\r
+gdouble   p2tr_math_length_sq  (gdouble x1,
+                                gdouble y1,
+                                gdouble x2,
+                                gdouble y2);
+\r
+gdouble   p2tr_math_length_sq2 (const P2trVector2 *pt1,
+                                const P2trVector2 *pt2);\r
+
+typedef enum\r
+{\r
+  P2TR_INTRIANGLE_OUT = -1,\r
+  P2TR_INTRIANGLE_ON = 0,\r
+  P2TR_INTRIANGLE_IN = 1\r
+} P2trInTriangle;
+
+P2trInTriangle p2tr_math_intriangle (const P2trVector2 *A,\r
+                                     const P2trVector2 *B,\r
+                                     const P2trVector2 *C,\r
+                                     const P2trVector2 *P);\r
+
+typedef enum\r
+{\r
+  P2TR_ORIENTATION_CW = -1,\r
+  P2TR_ORIENTATION_LINEAR = 0,\r
+  P2TR_ORIENTATION_CCW = 1\r
+} P2trOrientation;\r
+
+P2trOrientation p2tr_math_orient2d (const P2trVector2 *A,
+                                    const P2trVector2 *B,
+                                    const P2trVector2 *C);
+
+typedef enum
+{
+  P2TR_INCIRCLE_IN,
+  P2TR_INCIRCLE_ON,
+  P2TR_INCIRCLE_OUT
+} P2trInCircle;
+
+P2trInCircle p2tr_math_incircle (const P2trVector2 *A,\r
+                                 const P2trVector2 *B,\r
+                                 const P2trVector2 *C,\r
+                                 const P2trVector2 *D);\r
+\r
+#endif
\ No newline at end of file

diff --git a/refine/utils.h b/refine/utils.h
index a8461558c5c43df3e5394ec0d66a1771db7c5e6b..b98931b2bdbe9da3c2a352bdb4a390d6b52df2b8 100755 (executable)
--- a/refine/utils.h
+++ b/refine/utils.h
@@ -33,10 +33,10 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
-#ifndef __POLY2TRI_C_REFINE_UTILS_H__
-#define        __POLY2TRI_C_REFINE_UTILS_H__
+#ifndef __P2TC_REFINE_UTILS_H__
+#define __P2TC_REFINE_UTILS_H__
 
-#ifdef __cplusplus
+#ifdef __cplusplus
 extern "C"
 {
 #endif
@@ -47,7 +47,7 @@ extern "C"
    * http://developer.gnome.org/glib/2.29/glib-Hash-Tables.html
    */
 
-  typedef GHashTable     P2tRHashSet;
+  typedef GHashTable     P2trHashSet;
   typedef GHashTableIter P2trHashSetIter;
 
 #define p2tr_hash_set_set_new(hash_func, equal_func, destroy) g_hash_table_new_full ((hash_func), (equal_func), (destroy),NULL)
@@ -64,8 +64,8 @@ extern "C"
 
 #define foreach(iter,list) for ((iter) = (list); (iter) != NULL; (iter) = (iter)->next)
 
-#ifdef __cplusplus
+#ifdef __cplusplus
 }
 #endif
 
-#endif /* UTILS_H */
+#endif
\ No newline at end of file

diff --git a/refine/vector2.c b/refine/vector2.c
index 8a19e8d9f8a8884047edef01948ed4ce6f50b095..fde20b6a07aaffa2d3053920bfe62fe5d6ed4f0f 100644 (file)
--- a/refine/vector2.c
+++ b/refine/vector2.c
@@ -1,17 +1,17 @@
 #include <math.h>\r
 #include <glib.h>\r
-#include "basic-geometry.h"\r
+#include "vector2.h"\r
 \r
 gdouble\r
-p2tr_vector2_dot (P2trVector2 *a,\r
-                  P2trVector2 *b)\r
+p2tr_vector2_dot (const P2trVector2 *a,\r
+                  const P2trVector2 *b)\r
 {\r
   return a->x * b->x + a->y * b->y;\r
 }\r
 \r
 gboolean\r
-p2tr_vector2_is_same (P2trVector2 *a,\r
-                      P2trVector2 *b)\r
+p2tr_vector2_is_same (const P2trVector2 *a,\r
+                      const P2trVector2 *b)\r
 {\r
   if (a == NULL || b == NULL)\r
     return ! ((a == NULL) ^ (b == NULL));\r
@@ -20,23 +20,23 @@ p2tr_vector2_is_same (P2trVector2 *a,
 }\r
 \r
 void\r
-p2tr_vector2_sub (P2trVector2 *a,\r
-                  P2trVector2 *b,\r
-                  P2trVector2 *dest)\r
+p2tr_vector2_sub (const P2trVector2 *a,\r
+                  const P2trVector2 *b,\r
+                  P2trVector2       *dest)\r
 {\r
   dest->x = a->x - b->x;\r
   dest->y = a->y - b->y;\r
 }\r
 \r
 gdouble\r
-p2tr_vector2_norm (P2trVector2 *v)\r
+p2tr_vector2_norm (const P2trVector2 *v)\r
 {\r
   return sqrt (v->x * v->x + v->y * v->y);\r
 }\r
 \r
 void\r
-p2tr_vector2_copy (P2trVector2 *dest,\r
-                   P2trVector2 *src)\r
+p2tr_vector2_copy (P2trVector2       *dest,\r
+                   const P2trVector2 *src)\r
 {\r
   dest->x = src->x;\r
   dest->y = src->y;\r

diff --git a/refine/vector2.h b/refine/vector2.h
index 9d6a58cd70b885105a301fc367b443094ad62038..36f39452699725607b36e90d27b7508a2262fc66 100644 (file)
--- a/refine/vector2.h
+++ b/refine/vector2.h
@@ -8,10 +8,10 @@ typedef struct {
   gdouble y;\r
 } P2trVector2;\r
 \r
-gdouble       p2tr_vector2_dot       (P2trVector2 *a, P2trVector2 *b);\r
-gboolean      p2tr_vector2_is_same   (P2trVector2 *a, P2trVector2 *b);\r
-void          p2tr_vector2_sub       (P2trVector2 *a, P2trVector2 *b, P2trVector2 *dest);\r
-gdouble       p2tr_vector2_norm      (P2trVector2 *v);\r
-void          p2tr_vector2_copy      (P2trVector2 *dest, P2trVector2 *src);\r
+gdouble       p2tr_vector2_dot       (const P2trVector2 *a, const P2trVector2 *b);\r
+gboolean      p2tr_vector2_is_same   (const P2trVector2 *a, const P2trVector2 *b);\r
+void          p2tr_vector2_sub       (const P2trVector2 *a, const P2trVector2 *b, P2trVector2 *dest);\r
+gdouble       p2tr_vector2_norm      (const P2trVector2 *v);\r
+void          p2tr_vector2_copy      (P2trVector2 *dest, const P2trVector2 *src);\r
 \r
 #endif
\ No newline at end of file