ref: 7059938e9449f161805c79a7129da5ede7a3dafe
parent: a584d29458e3fec3e1e3162068fc5ab516a2d00e
author: rodri <rgl@antares-labs.eu>
date: Mon Aug 19 17:30:14 EDT 2024
libgeometry: add lineXsphere function
--- a/sys/include/geometry.h
+++ b/sys/include/geometry.h
@@ -77,6 +77,7 @@
Point3 crossvec3(Point3, Point3);
double vec3len(Point3);
Point3 normvec3(Point3);
+int lineXsphere(Point3*, Point3, Point3, Point3, double, int);
/* Matrix */
void identity(Matrix);
--- a/sys/man/2/geometry
+++ b/sys/man/2/geometry
@@ -88,6 +88,8 @@
Point3 crossvec3(Point3 a, Point3 b);
double vec3len(Point3 v);
Point3 normvec3(Point3 v);
+int lineXsphere(Point3 *rp, Point3 p0, Point3 p1,
+ Point3 c, double rad, int isaray);
.PB
/* Matrix */
void identity(Matrix m);
@@ -366,6 +368,25 @@
Normalizes the vector
.I v
and returns a new 3D point.
+.TP
+.B lineXsphere(\fIrp\fP,\fIp0\fP,\fIp1\fP,\fIc\fP,\fIrad\fP,\fIisaray\fP)
+Finds the intersection between the line defined by
+.I p0
+and
+.I p1
+and the circle with center at
+.I c
+and a radius of
+.IR rad .
+If
+.I isaray
+is not zero, the line will be interpreted as a ray directed from
+.I p0
+to
+.IR p1 .
+The function returns the number of intersections (up to two) and, if
+.I rp
+is not nil, it is filled with the resulting points.
.SS Matrices
.TP
Name
--- a/sys/src/libgeometry/point.c
+++ b/sys/src/libgeometry/point.c
@@ -216,3 +216,35 @@
return Pt3(0,0,0,0);
return Pt3(v.x/len, v.y/len, v.z/len, 0);
}
+
+int
+lineXsphere(Point3 *rp, Point3 p0, Point3 p1, Point3 c, double r, int isaray)
+{+ Point3 u, dp;
+ double u·dp, Δ, d;
+
+ u = normvec3(subpt3(p1, p0));
+ dp = subpt3(p1, c);
+ u·dp = dotvec3(u, dp);
+ if(isaray && u·dp > 0) /* ignore what's behind */
+ return 0;
+
+ Δ = u·dp*u·dp - dotvec3(dp, dp) + r*r;
+ if(Δ < 0) /* no intersection */
+ return 0;
+ else if(Δ == 0){ /* tangent */+ if(rp != nil){+ d = -u·dp;
+ rp[0] = addpt3(p0, mulpt3(u, d));
+ }
+ return 1;
+ }else{ /* secant */+ if(rp != nil){+ d = -u·dp + sqrt(Δ);
+ rp[0] = addpt3(p0, mulpt3(u, d));
+ d = -u·dp - sqrt(Δ);
+ rp[1] = addpt3(p0, mulpt3(u, d));
+ }
+ return 2;
+ }
+}
--
⑨