bring in the edge function and ptinpoly algorithm.

2 files changed, 42 insertions, 0 deletions

diff --git a/geometry.h b/geometry.h
index 3d99b41..7b4402d 100644
--- a/geometry.h
+++ b/geometry.h
@@ -7,6 +7,7 @@ typedef double Matrix3[4][4];
 typedef struct Quaternion Quaternion;
 typedef struct RFrame RFrame;
 typedef struct RFrame3 RFrame3;
+typedef struct Triangle2 Triangle2;
 typedef struct Triangle3 Triangle3;
 
 struct Point2 {
@@ -31,6 +32,11 @@ struct RFrame3 {
 	Point3 bx, by, bz;
 };
 
+struct Triangle2
+{
+	Point2 p0, p1, p2;
+};
+
 struct Triangle3 {
 	Point3 p0, p1, p2;
 };
diff --git a/point.c b/point.c
index 507bb7f..e891b02 100644
--- a/point.c
+++ b/point.c
@@ -75,6 +75,42 @@ normvec2(Point2 v)
 	return Pt2(v.x/len, v.y/len, 0);
 }
 
+/*
+ * the edge function, from:
+ *
+ * Juan Pineda, A Parallel Algorithm for Polygon Rasterization,
+ * Computer Graphics, Vol. 22, No. 8, August 1988
+ */
+int
+edgeptcmp(Point2 e0, Point2 e1, Point2 p)
+{
+	Point3 e0p, e01, r;
+
+	p = subpt2(p, e0);
+	e1 = subpt2(e1, e0);
+	e0p = Vec3(p.x,p.y,0);
+	e01 = Vec3(e1.x,e1.y,0);
+	r = crossvec3(e0p, e01);
+
+	/* clamp to avoid overflow */
+	return fclamp(r.z, -1, 1); /* e0.x*e1.y - e0.y*e1.x */
+}
+
+/*
+ * algorithm by W. Randolph Franklin
+ */
+int
+ptinpoly(Point2 p, Point2 *pts, ulong npts)
+{
+	int i, j, c;
+
+	for(i = c = 0, j = npts-1; i < npts; j = i++)
+		if(p.y < pts[i].y != p.y < pts[j].y &&
+			p.x < (pts[j].x - pts[i].x) * (p.y - pts[i].y)/(pts[j].y - pts[i].y) + pts[i].x)
+		c ^= 1;
+	return c;
+}
+
 /* 3D */
 
 Point3