diff --git a/sys/include/geometry.h b/sys/include/geometry.h index 715469b7c..2f45fa54d 100644 --- a/sys/include/geometry.h +++ b/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); diff --git a/sys/man/2/geometry b/sys/man/2/geometry index 757c0fb5d..f5efe08ca 100644 --- a/sys/man/2/geometry +++ b/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 diff --git a/sys/src/libgeometry/point.c b/sys/src/libgeometry/point.c index 72d54b9fc..4d091ab02 100644 --- a/sys/src/libgeometry/point.c +++ b/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; +}