I was curious and planned to do stuff needed for this in GLSL math style anyway. So here a different approach. Let consider this definition of your cone:
- create a set of basic geometry primitives
You need to support points,lines,triangles,convex triangulated mesh,spherical sector (cone).
- implement inside test between point and triangle,mesh,cone
for triangle the results of cross between any side and point – side origin should point on the same side of triangle (like normal). If not point is outside.
for convex mesh dot product between point-face origin and face normal pointing out should be <=0 for all faces.
for cone the point should be inside sphere radius and angle between cone axis and point-cone origin should be <= ang. again dot product can be used for this.
- implement closest line between basic primitives
this is like finding closest points on each primitive that forms a line. Its similar to perpendicular distance.
point-point its easy as they are the closest line.
point-line can be done using projection of point onto the line (dot product). However you need to bound the result so it is inside line and not extrapolated beond it.
point-triangle can be obtained as minimum of all the circumference lines vs point combinations and perpendicular distance to surface (dot product with triangle normal).
All the other combinations of primitives can be build from these basic ones.
- closest line between mesh and cone
simply use closest line between cone sphere center and mesh. If the line lies inside cone shorten it by sphere radius R. This will account all cap interactions.
Then test lines on surface of the cone so sample along its circumference starting on cone sphere center and ending on the outermost circle (edge between cone and cap). You van test also triangles instead if you need better precision.
- intersection between mesh and cone
this one is easy just compute closest lien between mesh and cone. And then test if its point on mesh side is inside cone or not.
check the
`bool intersect(convex_mesh m0,spherical_sector s0);`
implementation in the code below.
Here small C++/OpenGL example (using GLSL style math):
//---------------------------------------------------------------------------
//--- GL geometry -----------------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _gl_geometry_h
#define _gl_geometry_h
//---------------------------------------------------------------------------
const float deg=M_PI/180.0;
const float rad=180.0/M_PI;
float divide(float a,float b){ if (fabs(b)<1e-10) return 0.0; else return a/b; }
double divide(double a,double b){ if (fabs(b)<1e-10) return 0.0; else return a/b; }
#include "GLSL_math.h"
#include "List.h"
//---------------------------------------------------------------------------
class point
{
public:
// cfg
vec3 p0;
point() {}
point(point& a) { *this=a; }
~point() {}
point* operator = (const point *a) { *this=*a; return this; }
//point* operator = (const point &a) { ...copy... return this; }
point(vec3 _p0)
{
p0=_p0;
compute();
}
void compute(){};
void draw()
{
glBegin(GL_POINTS);
glVertex3fv(p0.dat);
glEnd();
}
};
//---------------------------------------------------------------------------
class axis
{
public:
// cfg
vec3 p0,dp;
axis() {}
axis(axis& a) { *this=a; }
~axis() {}
axis* operator = (const axis *a) { *this=*a; return this; }
//axis* operator = (const axis &a) { ...copy... return this; }
axis(vec3 _p0,vec3 _dp)
{
p0=_p0;
dp=_dp;
compute();
}
void compute()
{
dp=normalize(dp);
}
void draw()
{
vec3 p; p=p0+100.0*dp;
glBegin(GL_LINES);
glVertex3fv(p0.dat);
glVertex3fv(p .dat);
glEnd();
}
};
//---------------------------------------------------------------------------
class line
{
public:
// cfg
vec3 p0,p1;
// computed
float l;
vec3 dp;
line() {}
line(line& a) { *this=a; }
~line() {}
line* operator = (const line *a) { *this=*a; return this; }
//line* operator = (const line &a) { ...copy... return this; }
line(vec3 _p0,vec3 _p1)
{
p0=_p0;
p1=_p1;
compute();
}
void swap()
{
vec3 p=p0; p0=p1; p1=p;
}
void compute()
{
dp=p1-p0;
l=length(dp);
}
void draw()
{
glBegin(GL_LINES);
glVertex3fv(p0.dat);
glVertex3fv(p1.dat);
glEnd();
}
};
//---------------------------------------------------------------------------
class triangle
{
public:
// cfg
vec3 p0,p1,p2;
// computed
vec3 n;
triangle() {}
triangle(triangle& a) { *this=a; }
~triangle() {}
triangle* operator = (const triangle *a) { *this=*a; return this; }
//triangle* operator = (const triangle &a) { ...copy... return this; }
triangle(vec3 _p0,vec3 _p1,vec3 _p2)
{
p0=_p0;
p1=_p1;
p2=_p2;
compute();
}
void swap()
{
vec3 p=p1; p1=p2; p2=p;
n=-n;
}
void compute()
{
n=normalize(cross(p1-p0,p2-p1));
}
void draw()
{
glBegin(GL_TRIANGLES);
glNormal3fv(n.dat);
glVertex3fv(p0.dat);
glVertex3fv(p1.dat);
glVertex3fv(p2.dat);
glEnd();
}
};
//---------------------------------------------------------------------------
class convex_mesh
{
public:
// cfg
List<triangle> tri;
// computed
vec3 p0; // center
convex_mesh() { tri.num=0; }
convex_mesh(convex_mesh& a) { *this=a; }
~convex_mesh() {}
convex_mesh* operator = (const convex_mesh *a) { *this=*a; return this; }
//convex_mesh* operator = (const convex_mesh &a) { ...copy... return this; }
void init_box(vec3 _p0,vec3 _u,vec3 _v,vec3 _w) // center, half sizes
{
const vec3 p[8]=
{
_p0-_u+_v-_w,
_p0+_u+_v-_w,
_p0+_u-_v-_w,
_p0-_u-_v-_w,
_p0-_u-_v+_w,
_p0+_u-_v+_w,
_p0+_u+_v+_w,
_p0-_u+_v+_w,
};
const int ix[36]=
{
0,1,2,0,2,3,
4,5,6,4,6,7,
3,2,5,3,5,4,
2,1,6,2,6,5,
1,0,7,1,7,6,
0,3,4,0,4,7,
};
tri.num=0;
for (int i=0;i<36;i+=3) tri.add(triangle(p[ix[i+0]],p[ix[i+1]],p[ix[i+2]]));
compute();
}
void compute()
{
int i,n;
p0=vec3(0.0,0.0,0.0);
if (!tri.num) return;
for (i=0,n=0;i<tri.num;i++,n+=3)
{
p0+=tri.dat[i].p0;
p0+=tri.dat[i].p1;
p0+=tri.dat[i].p2;
} p0/=float(n);
for (i=0;i<tri.num;i++)
if (dot(tri.dat[i].p0-p0,tri.dat[i].n)<0.0)
tri.dat[i].swap();
}
void draw()
{
int i;
glBegin(GL_TRIANGLES);
for (i=0;i<tri.num;i++) tri.dat[i].draw();
glEnd();
}
};
//---------------------------------------------------------------------------
class spherical_sector
{
public:
// cfg
vec3 p0,p1;
float ang;
// computed
vec3 dp;
float r,R;
spherical_sector() {}
spherical_sector(spherical_sector& a) { *this=a; }
~spherical_sector() {}
spherical_sector* operator = (const spherical_sector *a) { *this=*a; return this; }
//spherical_sector* operator = (const spherical_sector &a) { ...copy... return this; }
spherical_sector(vec3 _p0,vec3 _p1,float _ang)
{
p0=_p0;
p1=_p1;
ang=_ang;
compute();
}
void compute()
{
dp=p1-p0;
R=length(dp);
r=R*tan(ang);
}
void draw()
{
const int N=32;
const int M=16;
vec3 pnt[M][N]; // points
vec3 n0[N]; // normals for cine
vec3 n1[M][N]; // normals for cap
int i,j;
float a,b,da,db,ca,sa,cb,sb;
vec3 q,u,v,w;
// basis vectors
w=normalize(dp); u=vec3(1.0,0.0,0.0);
if (fabs(dot(u,w))>0.75) u=vec3(0.0,1.0,0.0);
v=cross(u,w);
u=cross(v,w);
u=normalize(u);
v=normalize(v);
// compute tables
da=2.0*M_PI/float(N-1);
db=ang/float(M-1);
for (a=0.0,i=0;i<N;i++,a+=da)
{
ca=cos(a);
sa=sin(a);
n0[i]=u*ca+v*sa;
for (b=0.0,j=0;j<M;j++,b+=db)
{
cb=cos(b);
sb=sin(b);
q=vec3(ca*sb,sa*sb,cb);
pnt[j][i]=p0+((q.x*u+q.y*v+q.z*w)*R);
n1[j][i]=normalize(pnt[j][i]);
}
}
// render
glBegin(GL_TRIANGLES);
for (i=1,j=M-1;i<N;i++)
{
glNormal3fv(n0[i].dat); // p0 should have average 0.5*(n0[i]+n0[i-1]) as nomal
glVertex3fv(p0.dat);
glVertex3fv(pnt[j][i+0].dat);
glNormal3fv(n0[i-1].dat);
glVertex3fv(pnt[j][i-1].dat);
glNormal3fv( n1[0][0].dat);
glVertex3fv(pnt[0][0].dat);
glNormal3fv( n1[1][i-1].dat);
glVertex3fv(pnt[1][i-1].dat);
glNormal3fv( n1[1][i+0].dat);
glVertex3fv(pnt[1][i+0].dat);
}
glEnd();
glBegin(GL_QUADS);
for (i=0;i<N;i++)
for (j=2;j<M;j++)
{
glNormal3fv( n1[j-1][i+0].dat);
glVertex3fv(pnt[j-1][i+0].dat);
glNormal3fv( n1[j-1][i-1].dat);
glVertex3fv(pnt[j-1][i-1].dat);
glNormal3fv( n1[j+0][i-1].dat);
glVertex3fv(pnt[j+0][i-1].dat);
glNormal3fv( n1[j+0][i+0].dat);
glVertex3fv(pnt[j+0][i+0].dat);
}
glEnd();
}
};
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
bool inside(point p0,triangle t0);
bool inside(point p0,convex_mesh m0);
bool inside(point p0,spherical_sector s0);
//---------------------------------------------------------------------------
line closest(point p0,axis a0);
line closest(point p0,line l0);
line closest(point p0,triangle t0);
line closest(point p0,convex_mesh m0);
//---------------------------------------------------------------------------
line closest(axis a0,point p0);
line closest(axis a0,axis a1);
line closest(axis a0,line l1);
line closest(axis a0,triangle t0);
line closest(axis a0,convex_mesh m0);
//---------------------------------------------------------------------------
line closest(line l0,point p0);
line closest(line l0,axis a0);
line closest(line l0,line l1);
line closest(line l0,triangle t0);
line closest(line l0,convex_mesh m0);
//---------------------------------------------------------------------------
line closest(triangle t0,point p0);
line closest(triangle t0,axis a0);
line closest(triangle t0,line l0);
line closest(triangle t0,triangle t1);
line closest(triangle t0,convex_mesh m0);
//---------------------------------------------------------------------------
line closest(convex_mesh m0,point p0);
line closest(convex_mesh m0,axis a0);
line closest(convex_mesh m0,line l0);
line closest(convex_mesh m0,triangle t0);
line closest(convex_mesh m0,spherical_sector s0);
//---------------------------------------------------------------------------
bool intersect(convex_mesh m0,spherical_sector s0);
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
bool inside(point p0,triangle t0)
{
if (fabs(dot(p0.p0-t0.p0,t0.n))>1e-6) return false;
float d0,d1,d2;
d0=dot(t0.n,cross(p0.p0-t0.p0,t0.p1-t0.p0));
d1=dot(t0.n,cross(p0.p0-t0.p1,t0.p2-t0.p1));
d2=dot(t0.n,cross(p0.p0-t0.p2,t0.p0-t0.p2));
if (d0*d1<-1e-6) return false;
if (d0*d2<-1e-6) return false;
if (d1*d2<-1e-6) return false;
return true;
}
bool inside(point p0,convex_mesh m0)
{
for (int i=0;i<m0.tri.num;i++)
if (dot(p0.p0-m0.tri.dat[i].p0,m0.tri.dat[i].n)>0.0)
return false;
return true;
}
bool inside(point p0,spherical_sector s0)
{
float t,l;
vec3 u;
u=p0.p0-s0.p0;
l=length(u);
if (l>s0.R) return false;
t=divide(dot(u,s0.dp),(l*s0.R));
if (t<cos(s0.ang)) return false;
return true;
}
//---------------------------------------------------------------------------
line closest(point p0,axis a0){ return line(p0.p0,a0.p0+(a0.dp*dot(p0.p0-a0.p0,a0.dp))); }
line closest(point p0,line l0)
{
float t=dot(p0.p0-l0.p0,l0.dp);
if (t<0.0) t=0.0;
if (t>1.0) t=1.0;
return line(p0.p0,l0.p0+(l0.dp*t));
}
line closest(point p0,triangle t0)
{
float t;
point p;
line cl,ll;
cl.l=1e300;
t=dot(p0.p0-t0.p0,t0.n); p=p0.p0-t*t0.n; if ((fabs(t)>1e-6)&&(inside(p,t0))){ ll=line(p0.p0,p.p0); if (cl.l>ll.l) cl=ll; }
ll=closest(p0,line(t0.p0,t0.p1)); if (cl.l>ll.l) cl=ll;
ll=closest(p0,line(t0.p1,t0.p2)); if (cl.l>ll.l) cl=ll;
ll=closest(p0,line(t0.p2,t0.p0)); if (cl.l>ll.l) cl=ll;
return cl;
}
line closest(point p0,convex_mesh m0)
{
int i;
line cl,ll;
cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0)); cl.l=1e300;
for (i=0;i<m0.tri.num;i++)
{
ll=closest(p0,m0.tri.dat[i]);
if (cl.l>ll.l) cl=ll;
}
return cl;
}
//---------------------------------------------------------------------------
line closest(axis a0,point p0){ line cl; cl=closest(p0,a0); cl.swap(); return cl; }
line closest(axis a0,axis a1)
{
vec3 u=a0.dp;
vec3 v=a1.dp;
vec3 w=a0.p0-a1.p0;
float a=dot(u,u); // always >= 0
float b=dot(u,v);
float c=dot(v,v); // always >= 0
float d=dot(u,w);
float e=dot(v,w);
float D=a*c-b*b; // always >= 0
float t0,t1;
// compute the line parameters of the two closest points
if (D<1e-6) // the lines are almost parallel
{
t0=0.0;
t1=(b>c ? d/b : e/c); // use the largest denominator
}
else{
t0=(b*e-c*d)/D;
t1=(a*e-b*d)/D;
}
return line(a0.p0+(a0.dp*t0),a1.p0+(a1.dp*t1));
}
line closest(axis a0,line l1)
{
vec3 u=a0.dp;
vec3 v=l1.dp;
vec3 w=a0.p0-l1.p0;
float a=dot(u,u); // always >= 0
float b=dot(u,v);
float c=dot(v,v); // always >= 0
float d=dot(u,w);
float e=dot(v,w);
float D=a*c-b*b; // always >= 0
float t0,t1;
// compute the line parameters of the two closest points
if (D<1e-6) // the lines are almost parallel
{
t0=0.0;
t1=(b>c ? d/b : e/c); // use the largest denominator
}
else{
t0=(b*e-c*d)/D;
t1=(a*e-b*d)/D;
}
if (t1<0.0) t1=0.0;
if (t1>1.0) t1=1.0;
return line(a0.p0+(a0.dp*t0),l1.p0+(l1.dp*t1));
}
line closest(axis a0,triangle t0)
{
line cl,ll;
cl=closest(a0,line(t0.p0,t0.p1));
ll=closest(a0,line(t0.p1,t0.p2)); if (cl.l>ll.l) cl=ll;
ll=closest(a0,line(t0.p2,t0.p0)); if (cl.l>ll.l) cl=ll;
return cl;
}
line closest(axis a0,convex_mesh m0)
{
int i;
line cl,ll;
cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0)); cl.l=1e300;
for (i=0;i<m0.tri.num;i++)
{
ll=closest(a0,m0.tri.dat[i]);
if (cl.l>ll.l) cl=ll;
}
return cl;
}
//---------------------------------------------------------------------------
line closest(line l0,point p0){ line cl; cl=closest(p0,l0); cl.swap(); return cl; }
line closest(line l0,axis a0) { line cl; cl=closest(a0,l0); cl.swap(); return cl; }
line closest(line l0,line l1)
{
vec3 u=l0.p1-l0.p0;
vec3 v=l1.p1-l1.p0;
vec3 w=l0.p0-l1.p0;
float a=dot(u,u); // always >= 0
float b=dot(u,v);
float c=dot(v,v); // always >= 0
float d=dot(u,w);
float e=dot(v,w);
float D=a*c-b*b; // always >= 0
float t0,t1;
// compute the line parameters of the two closest points
if (D<1e-6) // the lines are almost parallel
{
t0=0.0;
t1=(b>c ? d/b : e/c); // use the largest denominator
}
else{
t0=(b*e-c*d)/D;
t1=(a*e-b*d)/D;
}
if (t0<0.0) t0=0.0;
if (t0>1.0) t0=1.0;
if (t1<0.0) t1=0.0;
if (t1>1.0) t1=1.0;
return line(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));
}
line closest(line l0,triangle t0)
{
float t;
point p;
line cl,ll;
cl.l=1e300;
t=dot(l0.p0-t0.p0,t0.n); p=l0.p0-t*t0.n; if ((fabs(t)>1e-6)&&(inside(p,t0))){ ll=line(l0.p0,p.p0); if (cl.l>ll.l) cl=ll; }
t=dot(l0.p1-t0.p0,t0.n); p=l0.p1-t*t0.n; if ((fabs(t)>1e-6)&&(inside(p,t0))){ ll=line(l0.p1,p.p0); if (cl.l>ll.l) cl=ll; }
ll=closest(l0,line(t0.p0,t0.p1)); if (cl.l>ll.l) cl=ll;
ll=closest(l0,line(t0.p1,t0.p2)); if (cl.l>ll.l) cl=ll;
ll=closest(l0,line(t0.p2,t0.p0)); if (cl.l>ll.l) cl=ll;
return cl;
}
line closest(line l0,convex_mesh m0)
{
int i;
line cl,ll;
cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0)); cl.l=1e300;
for (i=0;i<m0.tri.num;i++)
{
ll=closest(l0,m0.tri.dat[i]);
if (cl.l>ll.l) cl=ll;
}
return cl;
}
//---------------------------------------------------------------------------
line closest(triangle t0,point p0){ line cl; cl=closest(p0,t0); cl.swap(); return cl; }
line closest(triangle t0,axis a0) { line cl; cl=closest(a0,t0); cl.swap(); return cl; }
line closest(triangle t0,line l0) { line cl; cl=closest(l0,t0); cl.swap(); return cl; }
line closest(triangle t0,triangle t1)
{
float t;
point p;
line l0,l1,l2,l3,l4,l5,cl,ll;
l0=line(t0.p0,t0.p1); l3=line(t1.p0,t1.p1);
l1=line(t0.p1,t0.p2); l4=line(t1.p1,t1.p2);
l2=line(t0.p2,t0.p0); l5=line(t1.p2,t1.p0);
cl.l=1e300;
t=dot(t0.p0-t1.p0,t1.n); p=t0.p0-t*t1.n; if ((fabs(t)>1e-6)&&(inside(p,t1))){ ll=line(t0.p0,p.p0); if (cl.l>ll.l) cl=ll; }
t=dot(t0.p1-t1.p0,t1.n); p=t0.p1-t*t1.n; if ((fabs(t)>1e-6)&&(inside(p,t1))){ ll=line(t0.p1,p.p0); if (cl.l>ll.l) cl=ll; }
t=dot(t0.p2-t1.p0,t1.n); p=t0.p2-t*t1.n; if ((fabs(t)>1e-6)&&(inside(p,t1))){ ll=line(t0.p2,p.p0); if (cl.l>ll.l) cl=ll; }
t=dot(t1.p0-t0.p0,t0.n); p=t1.p0-t*t0.n; if ((fabs(t)>1e-6)&&(inside(p,t0))){ ll=line(p.p0,t1.p0); if (cl.l>ll.l) cl=ll; }
t=dot(t1.p1-t0.p0,t0.n); p=t1.p1-t*t0.n; if ((fabs(t)>1e-6)&&(inside(p,t0))){ ll=line(p.p0,t1.p1); if (cl.l>ll.l) cl=ll; }
t=dot(t1.p2-t0.p0,t0.n); p=t1.p2-t*t0.n; if ((fabs(t)>1e-6)&&(inside(p,t0))){ ll=line(p.p0,t1.p2); if (cl.l>ll.l) cl=ll; }
ll=closest(l0,l3); if (cl.l>ll.l) cl=ll;
ll=closest(l0,l4); if (cl.l>ll.l) cl=ll;
ll=closest(l0,l5); if (cl.l>ll.l) cl=ll;
ll=closest(l1,l3); if (cl.l>ll.l) cl=ll;
ll=closest(l1,l4); if (cl.l>ll.l) cl=ll;
ll=closest(l1,l5); if (cl.l>ll.l) cl=ll;
ll=closest(l2,l3); if (cl.l>ll.l) cl=ll;
ll=closest(l2,l4); if (cl.l>ll.l) cl=ll;
ll=closest(l2,l5); if (cl.l>ll.l) cl=ll;
return cl;
}
line closest(triangle t0,convex_mesh m0)
{
int i;
line cl,ll;
cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0)); cl.l=1e300;
for (i=0;i<m0.tri.num;i++)
{
ll=closest(m0.tri.dat[i],t0);
if (cl.l>ll.l) cl=ll;
}
return cl;
}
//---------------------------------------------------------------------------
line closest(convex_mesh m0,point p0) { line cl; cl=closest(p0,m0); cl.swap(); return cl; }
line closest(convex_mesh m0,axis a0) { line cl; cl=closest(a0,m0); cl.swap(); return cl; }
line closest(convex_mesh m0,line l0) { line cl; cl=closest(l0,m0); cl.swap(); return cl; }
line closest(convex_mesh m0,triangle t0){ line cl; cl=closest(t0,m0); cl.swap(); return cl; }
line closest(convex_mesh m0,convex_mesh m1)
{
int i0,i1;
line cl,ll;
cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0)); cl.l=1e300;
for (i0=0;i0<m0.tri.num;i0++)
for (i1=0;i1<m1.tri.num;i1++)
{
ll=closest(m0.tri.dat[i0],m1.tri.dat[i1]);
if (cl.l>ll.l) cl=ll;
}
return cl;
}
line closest(convex_mesh m0,spherical_sector s0)
{
int i,N=18;
float a,da,ca,sa,cb,sb;
vec3 u,v,w,q;
line cl,ll;
// cap
ll=closest(m0,point(s0.p0)); // sphere
if (dot(ll.dp,s0.dp)/(ll.l*s0.R)>=cos(s0.ang)) // cap
ll=line(ll.p0,ll.p1+(ll.dp*s0.R/ll.l));
cl=ll;
// cone
w=normalize(s0.dp); u=vec3(1.0,0.0,0.0);
if (fabs(dot(u,w))>0.75) u=vec3(0.0,1.0,0.0);
v=cross(u,w);
u=cross(v,w);
u=normalize(u)*s0.r;
v=normalize(v)*s0.r;
da=2.0*M_PI/float(N-1);
cb=cos(s0.ang);
sb=sin(s0.ang);
for (a=0.0,i=0;i<N;i++)
{
ca=cos(a);
sa=sin(a);
q=vec3(ca*sb,sa*sb,cb);
q=s0.p0+((q.x*u+q.y*v+q.z*w)*s0.R);
ll=line(s0.p0,q);
ll=closest(m0,ll);
if (cl.l>ll.l) cl=ll;
}
return cl;
}
//---------------------------------------------------------------------------
bool intersect(convex_mesh m0,spherical_sector s0)
{
line cl;
cl=closest(m0,s0);
if (cl.l<=1e-6) return true;
if (inside(cl.p0,s0)) return true;
return false;
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
The GLSL math can be created by this or use GLM or whatever else instead.
I also used mine dynamic list template (just to store the list of triangles in mesh) so:
List<double> xxx; is the same as double xxx[];
xxx.add(5); adds 5 to end of the list
xxx[7] access array element (safe)
xxx.dat[7] access array element (unsafe but fast direct access)
xxx.num is the actual used size of the array
xxx.reset() clears the array and set xxx.num=0
xxx.allocate(100) preallocate space for 100 items
You can use whatever list you have at disposal.
And here test preview testing correctness of this:
The cone is rotating and changing color according to result of intersect test. The yellow line is the closest line result.
I busted this for fun during this weekend so its not yet extensively tested and there still might be unhandled edge cases.
I wanted the code to be as readable as I could so its not optimized at all. Also I did not comment much (as the low level primitives and basic vector math should be obvious enough if not you should learn first before implementing stuff like this)
[edit1]
Looks like I mess some closest test ignoring some edge cases.. It need major rework (apply the fix for all related functions) which I do not have enough time for right now (will update the code once finished) so for now here a quick fix for line vs. line closest test only:
line3D closest(line3D l0,line3D l1)
{
vec3 u=l0.p1-l0.p0;
vec3 v=l1.p1-l1.p0;
vec3 w=l0.p0-l1.p0;
float a=dot(u,u); // always >= 0
float b=dot(u,v);
float c=dot(v,v); // always >= 0
float d=dot(u,w);
float e=dot(v,w);
float D=a*c-b*b; // always >= 0
float t0,t1;
point3D p;
line3D r,rr;
int f; // check distance perpendicular to: 1: l0, 2: l1
f=0; r.l=-1.0;
// compute the line3D parameters of the two closest points
if (D<acc_dot) f=3; // the lines are almost parallel
else{
t0=(b*e-c*d)/D;
t1=(a*e-b*d)/D;
if (t0<0.0){ f|=1; t0=0.0; }
if (t0>1.0){ f|=1; t0=1.0; }
if (t1<0.0){ f|=2; t1=0.0; }
if (t1>1.0){ f|=2; t1=1.0; }
r=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));
}
// check perpendicular distance from each endpoint marked in f
if (int(f&1))
{
t0=0.0;
t1=divide(dot(l0.p0-l1.p0,l1.dp),l1.l*l1.l);
if (t1<0.0) t1=0.0;
if (t1>1.0) t1=1.0;
rr=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));
if ((r.l<0.0)||(r.l>rr.l)) r=rr;
t0=1.0;
t1=divide(dot(l0.p1-l1.p0,l1.dp),l1.l*l1.l);
if (t1<0.0) t1=0.0;
if (t1>1.0) t1=1.0;
rr=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));
if ((r.l<0.0)||(r.l>rr.l)) r=rr;
}
if (int(f&2))
{
t1=0.0;
t0=divide(dot(l1.p0-l0.p0,l0.dp),l0.l*l0.l);
if (t0<0.0) t0=0.0;
if (t0>1.0) t0=1.0;
rr=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));
if ((r.l<0.0)||(r.l>rr.l)) r=rr;
t1=1.0;
t0=divide(dot(l1.p1-l0.p0,l0.dp),l0.l*l0.l);
if (t0<0.0) t0=0.0;
if (t0>1.0) t0=1.0;
rr=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));
if ((r.l<0.0)||(r.l>rr.l)) r=rr;
}
return r;
}