libgeometry: add ptincylinder and ptincone

sys/include/geometry.h

@@ -78,6 +78,8 @@ Point3 crossvec3(Point3, Point3);
 double vec3len(Point3);
 Point3 normvec3(Point3);
 int lineXsphere(Point3*, Point3, Point3, Point3, double, int);
+int ptincylinder(Point3, Point3, Point3, double);
+int ptincone(Point3, Point3, Point3, double);
 
 /* Matrix */
 void identity(Matrix);

sys/man/2/geometry

@@ -90,6 +90,8 @@ double vec3len(Point3 v)
 Point3 normvec3(Point3 v)
 int lineXsphere(Point3 *rp, Point3 p0, Point3 p1,
 		Point3 c, double rad, int isaray)
+int ptincylinder(Point3 p, Point3 p0, Point3 p1, double r)
+int ptincone(Point3 p, Point3 p0, Point3 p1, double br)
 .PB
 /* Matrix */
 void identity(Matrix m)
@@ -389,6 +391,24 @@ to
 The function returns the number of intersections (up to two) and, if
 .I rp
 is not nil, it is filled with the resulting points.
+.TP
+.B ptincylinder(\fIp\fP,\fIp0\fP,\fIp1\fP,\fIr\fP)
+Tests if the point
+.I p
+is inside a cylinder with an axis defined by the line between
+.I p0
+and
+.IR p1 ,
+and a radius r.
+.TP
+.B ptincone(\fIp\fP,\fIp0\fP,\fIp1\fP,\fIbr\fP)
+Tests if the point
+.I p
+is inside a cone with its apex at
+.I p0
+and its base centered at
+.IR p1 ,
+with a radius of br.
 .SS Matrices
 .TP
 Name

sys/src/libgeometry/point.c

@@ -248,3 +248,51 @@ lineXsphere(Point3 *rp, Point3 p0, Point3 p1, Point3 c, double r, int isaray)
 		return 2;
 	}
 }
+
+/*
+ * p is the point to test
+ * p0 and p1 are the centers of the circles at each end of the cylinder
+ * r is the radius of these circles
+ */
+int
+ptincylinder(Point3 p, Point3 p0, Point3 p1, double r)
+{
+	Point3 p01, p0p, p1p;
+	double h;
+
+	p01 = subpt3(p1, p0);
+	p0p = subpt3(p, p0);
+	p1p = subpt3(p, p1);
+	h = vec3len(p01);
+
+	if(h == 0)
+		return 0;
+
+	return dotvec3(p0p, p01) >= 0 &&
+		dotvec3(p1p, p01) <= 0 &&
+		vec3len(crossvec3(p0p, p01))/h <= r;
+}
+
+/*
+ * p is the point to test
+ * p0 is the apex
+ * p1 is the center of the base
+ * br is the radius of the base
+ */
+int
+ptincone(Point3 p, Point3 p0, Point3 p1, double br)
+{
+	Point3 p01, p0p;
+	double h, d, r;
+
+	p01 = subpt3(p1, p0);
+	p0p = subpt3(p, p0);
+	h = vec3len(p01);
+	d = dotvec3(p0p, normvec3(p01));
+
+	if(h == 0 || d < 0 || d > h)
+		return 0;
+
+	r = d/h * br;
+	return vec3len(crossvec3(p0p, p01))/h <= r;
+}