锥体碰撞 [英] Cone to box collision

查看:24
本文介绍了锥体碰撞的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我希望在一个圆锥体(带有一个圆底.所以它基本上是一个球体的一部分)和一个盒子之间实现碰撞检测.我对它是 AABB 或 OBB 并不太大​​惊小怪,因为转换应该足够简单.我找到的每个解决方案都使用三角锥,但我的锥更像是一个具有角度和距离的圆弧".

是否有一个简单的解决方案来进行这种碰撞检测?或者是做几种类型的测试的情况?IE.类似于在球体上获取交点,r 是我的圆锥距离,然后测试它们是否在某个角度内相交?

解决方案

我很好奇并计划用 GLSL 数学风格做一些需要的东西.所以这里有一个不同的方法.让我们考虑一下你的锥体的这个定义:

  1. 创建一组基本几何图元

您需要支持点、线、三角形、凸三角网格、扇形(圆锥).

  1. 在点和三角形、网格、圆锥之间进行内部测试

对于 triangle 任何边和点之间交叉的结果 - 边原点应该指向三角形的同一边(像正常一样).如果不是点在外面.

对于凸面网格点-面原点和面法线指向之间的点积对于所有面应<=0.

对于cone,点应该在球体半径内,圆锥轴和点圆锥原点之间的角度应该是<= ang.再次点积可以用于此.

  1. 在基本原语之间实现最接近的线

这就像在形成一条线的每个基元上寻找最近​​的点.它类似于垂直距离.

point-point 很简单,因为它们是最近的线.

point-line 可以使用点投影到线上(点积)来完成.但是,您需要对结果进行限制,使其位于线内而不是外推.

point-triangle 可以得到所有圆周线对点组合和到表面的垂直距离(三角形法线的点积)的最小值.

所有其他原语组合都可以从这些基本组合中构建.

  1. 网格和圆锥之间的最近线

只需使用圆锥球体中心和网格之间的最近线.如果线位于圆锥内,则将其缩短球体半径 R.这将考虑所有帽相互作用.

然后在锥体表面测试线,沿其圆周采样,从锥体球心开始,到最外圈(锥体和帽之间的边缘)结束.如果您需要更高的精度,您还可以测试三角形.

  1. 网格与圆锥的交点

这很容易,只需计算网格和锥体之间的最近连线即可.然后测试它在网格侧的点是否在圆锥内.

检查

 `bool intersect(convex_mesh m0,spherical_sector s0);`

在下面的代码中实现.

这里是 C++/OpenGL 小例子(使用 GLSL 风格的数学):

//---------------------------------------------------------------------//--- GL几何体-----------------------------------------------------------//---------------------------------------------------------------------------#ifndef _gl_geometry_h#define _gl_geometry_h//---------------------------------------------------------------------------const float deg=M_PI/180.0;const float rad=180.0/M_PI;浮动除法(浮动a,浮动b){如果(fabs(b)<1e-10)返回0.0;否则返回 a/b;}双除法(double a,double b){ if (fabs(b)<1e-10) return 0.0;否则返回 a/b;}#include "GLSL_math.h";#include "List.h";//---------------------------------------------------------------------------班级点{上市://配置文件vec3 p0;观点()     {}点(点& a){ *this = a;}〜点(){}point* operator = (const point *a) { *this=*a;返回这个;}//point* operator = (const point &a) { ...copy... return this;}点(vec3 _p0){p0=_p0;计算();}无效计算(){};无效绘制(){glBegin(GL_POINTS);glVertex3fv(p0.dat);glEnd();}};//---------------------------------------------------------------------------阶级轴心{上市://配置文件vec3 p0,dp;轴(){}轴(轴& a){ *this = a;}~轴(){}轴* 运算符 = (const 轴 *a) { *this=*a;返回这个;}//axis* operator = (const axis &a) { ...copy... return this;}轴(vec3 _p0,vec3 _dp){p0=_p0;dp=_dp;计算();}无效计算(){dp=标准化(dp);}无效绘制(){vec3 p;p=p0+100.0*dp;glBegin(GL_LINES);glVertex3fv(p0.dat);glVertex3fv(p .dat);glEnd();}};//---------------------------------------------------------------------------班级线{上市://配置文件vec3 p0,p1;//计算浮动 l;vec3 dp;线()  {}线(线& a){ *this = a;}〜线(){}line* operator = (const line *a) { *this=*a;返回这个;}//line* operator = (const line &a) { ...copy... return this;}线(vec3 _p0,vec3 _p1){p0=_p0;p1=_p1;计算();}无效交换(){vec3 p=p0;p0=p1;p1=p;}无效计算(){dp=p1-p0;l=长度(dp);}无效绘制(){glBegin(GL_LINES);glVertex3fv(p0.dat);glVertex3fv(p1.dat);glEnd();}};//---------------------------------------------------------------------------班级三角{上市://配置文件vec3 p0,p1,p2;//计算vec3 n;三角形()  {}三角形(三角形& a){ *this = a;}~三角形(){}三角形 * 运算符 = (const 三角形 *a) { *this=*a;返回这个;}//triangle* operator = (const triangle &a) { ...copy... return this;}三角形(vec3 _p0,vec3 _p1,vec3 _p2){p0=_p0;p1=_p1;p2=_p2;计算();}无效交换(){vec3 p=p1;p1=p2;p2=p;n=-n;}无效计算(){n=标准化(交叉(p1-p0,p2-p1));}无效绘制(){glBegin(GL_TRIANGLES);glNormal3fv(n.dat);glVertex3fv(p0.dat);glVertex3fv(p1.dat);glVertex3fv(p2.dat);glEnd();}};//---------------------------------------------------------------------------类convex_mesh{上市://配置文件列表<三角形>三;//计算vec3 p0;//中央凸网格(){tri.num=0;}凸网(凸网& a){ *this = a;}~convex_mesh() {}凸网格 * 运算符 = (const 凸网格 *a) { *this=*a;返回这个;}//convex_mesh* operator = (constconvex_mesh &a) { ...copy... return this;}void init_box(vec3 _p0,vec3 _u,vec3 _v,vec3 _w)//中心,一半大小{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,};三数=0;for (int i=0;i<36;i+=3) tri.add(triangle(p[ix[i+0]],p[ix[i+1]],p[ix[i+2]]));计算();}无效计算(){int i,n;p0=vec3(0.0,0.0,0.0);如果(!tri.num)返回;对于 (i=0,n=0;i0.75) u=vec3(0.0,1.0,0.0);v=cross(u,w);u=cross(v,w);u=归一化(u);v=归一化(v);//计算表da=2.0*M_PI/float(N-1);db=ang/float(M-1);对于 (a=0.0,i=0;i1e-6) 返回假;浮动 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));如果 (d0*d1<-1e-6) 返回假;如果 (d0*d2<-1e-6) 返回假;如果 (d1*d2<-1e-6) 返回假;返回真;}bool里面(点p0,convex_mesh m0){for (int i=0;i0.0)返回假;返回真;}布尔里面(点 p0,球面 s0){浮动 t,l;vec3 u;u=p0.p0-s0.p0;l=长度(u);如果(l>s0.R)返回假;t=divide(dot(u,s0.dp),(l*s0.R));if (t1e-6)&&(inside(p,t0))){ ll=line(p0.p0,p.p0);如果(cl.l>ll.l)cl=ll;}ll=最近的(p0,line(t0.p0,t0.p1));如果(cl.l>ll.l)cl=ll;ll=最近的(p0,line(t0.p1,t0.p2));如果(cl.l>ll.l)cl=ll;ll=最近的(p0,line(t0.p2,t0.p0));如果(cl.l>ll.l)cl=ll;返回 cl;}线最近(点 p0,convex_mesh m0){国际我;线 cl,ll;cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0));cl.l=1e300;对于 (i=0;ill.l)cl=ll;}返回 cl;}//---------------------------------------------------------------------------线最近(轴 a0,点 p0){ 线 cl;cl=最近的(p0,a0);cl.swap();返回 cl;}线最近(轴 a0,轴 a1){vec3 u=a0.dp;vec3 v=a1.dp;vec3 w=a0.p0-a1.p0;浮动 a=dot(u,u);//总是 >= 0浮动 b=点(u,v);浮动 c=点(v,v);//总是 >= 0浮动 d=点(u,w);浮动 e=点(v,w);浮动 D=a*c-b*b;//总是 >= 0浮动 t0,t1;//计算两个最近点的线参数if (D<1e-6)//线几乎平行{t0=0.0;t1=(b>c?d/b:e/c);//使用最大的分母}别的{t0=(b*e-c*d)/D;t1=(a*e-b*d)/D;}回车线(a0.p0+(a0.dp*t0),a1.p0+(a1.dp*t1));}线最近(轴 a0,线 l1){vec3 u=a0.dp;vec3 v=l1.dp;vec3 w=a0.p0-l1.p0;浮动 a=dot(u,u);//总是 >= 0浮动 b=点(u,v);浮动 c=点(v,v);//总是 >= 0浮动 d=点(u,w);浮动 e=点(v,w);浮动 D=a*c-b*b;//总是 >= 0浮动 t0,t1;//计算两个最近点的线参数if (D<1e-6)//线几乎平行{t0=0.0;t1=(b>c?d/b:e/c);//使用最大的分母}别的{t0=(b*e-c*d)/D;t1=(a*e-b*d)/D;}如果 (t1<0.0) t1=0.0;如果 (t1>1.0) t1=1.0;回车线(a0.p0+(a0.dp*t0),l1.p0+(l1.dp*t1));}线最近(轴 a0,三角形 t0){线 cl,ll;cl=closest(a0,line(t0.p0,t0.p1));ll=最近的(a0,line(t0.p1,t0.p2));如果(cl.l>ll.l)cl=ll;ll=最近的(a0,line(t0.p2,t0.p0));如果(cl.l>ll.l)cl=ll;返回 cl;}线最近(轴 a0,convex_mesh m0){国际我;线 cl,ll;cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0));cl.l=1e300;对于 (i=0;ill.l)cl=ll;}返回 cl;}//---------------------------------------------------------------------------线最近(线 l0,点 p0){ 线 cl;cl=最近的(p0,l0);cl.swap();返回 cl;}线最近(线 l0,轴 a0){ 线 cl;cl=最近的(a0,l0);cl.swap();返回 cl;}线最近(线 l0,线 l1){vec3 u=l0.p1-l0.p0;vec3 v=l1.p1-l1.p0;vec3 w=l0.p0-l1.p0;浮动 a=dot(u,u);//总是 >= 0浮动 b=点(u,v);浮动 c=点(v,v);//总是 >= 0浮动 d=点(u,w);浮动 e=点(v,w);浮动 D=a*c-b*b;//总是 >= 0浮动 t0,t1;//计算两个最近点的线参数if (D<1e-6)//线几乎平行{t0=0.0;t1=(b>c?d/b:e/c);//使用最大的分母}别的{t0=(b*e-c*d)/D;t1=(a*e-b*d)/D;}如果 (t0<0.0) t0=0.0;如果(t0>1.0)t0=1.0;如果 (t1<0.0) t1=0.0;如果 (t1>1.0) t1=1.0;回车线(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));}线最近(线 l0,三角形 t0){浮动 t;p点;线 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);如果(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);如果(cl.l>ll.l)cl=ll;}ll=最近的(l0,line(t0.p0,t0.p1));如果(cl.l>ll.l)cl=ll;ll=最近的(l0,line(t0.p1,t0.p2));如果(cl.l>ll.l)cl=ll;ll=最近的(l0,line(t0.p2,t0.p0));如果(cl.l>ll.l)cl=ll;返回 cl;}线最近(线 l0,convex_mesh m0){国际我;线 cl,ll;cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0));cl.l=1e300;对于 (i=0;ill.l)cl=ll;}返回 cl;}//---------------------------------------------------------------------------线最近(三角形 t0,点 p0){ 线 cl;cl=最近的(p0,t0);cl.swap();返回 cl;}线最近(三角形 t0,轴 a0){ 线 cl;cl=最近的(a0,t0);cl.swap();返回 cl;}线最近(三角形 t0,线 l0){ 线 cl;cl=最近的(l0,t0);cl.swap();返回 cl;}线最近(三角形 t0,三角形 t1){浮动 t;p点;线 l0,l1,l2,l3,l4,l5,cl,ll;l0=线(t0.p0,t0.p1);l3=线(t1.p0,t1.p1);l1=线(t0.p1,t0.p2);l4=线(t1.p1,t1.p2);l2=线(t0.p2,t0.p0);l5=线(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);如果(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);如果(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);如果(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);如果(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);如果(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);如果(cl.l>ll.l)cl=ll;}ll=最近的(l0,l3);如果(cl.l>ll.l)cl=ll;ll=最近的(l0,l4);如果(cl.l>ll.l)cl=ll;ll=最近的(l0,l5);如果(cl.l>ll.l)cl=ll;ll=最近的(l1,l3);如果(cl.l>ll.l)cl=ll;ll=最近的(l1,l4);如果(cl.l>ll.l)cl=ll;ll=最近的(l1,l5);如果(cl.l>ll.l)cl=ll;ll=最近的(l2,l3);如果(cl.l>ll.l)cl=ll;ll=最近的(l2,l4);如果(cl.l>ll.l)cl=ll;ll=最近的(l2,l5);如果(cl.l>ll.l)cl=ll;返回 cl;}线最近(三角形 t0,convex_mesh m0){国际我;线 cl,ll;cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0));cl.l=1e300;对于 (i=0;ill.l)cl=ll;}返回 cl;}//---------------------------------------------------------------------------线最近(凸网格 m0,点 p0){ 线 cl;cl=最近的(p0,m0);cl.swap();返回 cl;}线最近(凸网格 m0,轴 a0){ 线 cl;cl=最近的(a0,m0);cl.swap();返回 cl;}线最近(凸网格 m0,线 l0){ 线 cl;cl=最近的(l0,m0);cl.swap();返回 cl;}线最近(凸网格 m0,三角形 t0){ 线 cl;cl=最近的(t0,m0);cl.swap();返回 cl;}线最近(凸网格 m0,凸网格 m1){int i0,i1;线 cl,ll;cl=line(vec3(0.0,0.0,0.0),vec3(0.0,0.0,0.0));cl.l=1e300;对于 (i0=0;i0ll.l)cl=ll;}返回 cl;}线最近(凸网格 m0,球面 s0){int i,N=18;浮动a,da,ca,sa,cb,sb;vec3 u,v,w,q;线 cl,ll;//上限ll=最近(m0,点(s0.p0));//球体if (dot(ll.dp,s0.dp)/(ll.l*s0.R)>=cos(s0.ang))//上限ll=line(ll.p0,ll.p1+(ll.dp*s0.R/ll.l));cl=ll;//锥体w=标准化(s0.dp);u=vec3(1.0,0.0,0.0);如果 (fabs(dot(u,w))>0.75) u=vec3(0.0,1.0,0.0);v=cross(u,w);u=cross(v,w);u=标准化(u)*s0.r;v=归一化(v)*s0.r;da=2.0*M_PI/float(N-1);cb=cos(s0.ang);sb=sin(s0.ang);对于 (a=0.0,i=0;ill.l)cl=ll;}返回 cl;}//---------------------------------------------------------------------------布尔相交(凸网格 m0,球面 s0){线cl;cl=最近的(m0,s0);if (cl.l<=1e-6) 返回真;if (inside(cl.p0,s0)) 返回真;返回假;}//---------------------------------------------------------------------------#万一//---------------------------------------------------------------------------

GLSL 数学可以通过

锥体根据相交测试的结果旋转并改变颜色.黄线是最近的线结果.

我在这个周末为了好玩而破坏了它,因此它尚未经过广泛测试,并且仍然可能存在未处理的边缘情况.

我希望代码尽可能具有可读性,因此根本没有对其进行优化.我也没有过多评论(因为低级原语和基本向量数学应该足够明显,如果不是你应该在实现这样的东西之前先学习)

[edit1]

看起来我搞砸了一些 closest 测试而忽略了一些边缘情况..它需要大量返工(对所有相关功能应用修复),我现在没有足够的时间(将更新代码一旦完成)所以现在这里只针对线路与线路最近的测试进行快速修复:

line3D最近的(line3D l0,line3D l1){vec3 u=l0.p1-l0.p0;vec3 v=l1.p1-l1.p0;vec3 w=l0.p0-l1.p0;浮动 a=dot(u,u);//总是 >= 0浮动 b=点(u,v);浮动 c=点(v,v);//总是 >= 0浮动 d=点(u,w);浮动 e=点(v,w);浮动 D=a*c-b*b;//总是 >= 0浮动 t0,t1;point3D p;line3D r,rr;国际f;//检查垂直于:1: l0, 2: l1 的距离f=0;r.l=-1.0;//计算两个最近点的line3D参数如果 (D1.0){ f|=1;t0=1.0;}如果 (t1<0.0){ f|=2;t1=0.0;}如果 (t1>1.0){ f|=2;t1=1.0;}r=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));}//检查与 f 中标记的每个端点的垂直距离如果(整数(f&1)){t0=0.0;t1=divide(dot(l0.p0-l1.p0,l1.dp),l1.l*l1.l);如果 (t1<0.0) t1=0.0;如果 (t1>1.0) t1=1.0;rr=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));如果 ((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);如果 (t1<0.0) t1=0.0;如果 (t1>1.0) t1=1.0;rr=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));如果 ((r.l<0.0)||(r.l>rr.l)) r=rr;}如果 (int(f&2)){t1=0.0;t0=divide(dot(l1.p0-l0.p0,l0.dp),l0.l*l0.l);如果 (t0<0.0) t0=0.0;如果(t0>1.0)t0=1.0;rr=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));如果 ((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);如果 (t0<0.0) t0=0.0;如果(t0>1.0)t0=1.0;rr=line3D(l0.p0+(l0.dp*t0),l1.p0+(l1.dp*t1));如果 ((r.l<0.0)||(r.l>rr.l)) r=rr;}返回 r;}

I'm looking to implement collision detection between a cone (With a round bottom. So it's basically a slice of a sphere) and a box. I'm not too fussed about it being AABB or OBB because transforming should be simple enough. Every solution I find uses a triangular cone but my cone is more of an "arc" that has an angle and distance.

Is there a simple solution to doing this collision detection? Or is it a case of doing several types of tests? ie. something like getting intersection points on a sphere with r being my cone distance then testing if they intersect within an angle or something?

解决方案

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:

  1. create a set of basic geometry primitives

You need to support points,lines,triangles,convex triangulated mesh,spherical sector (cone).

  1. 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.

  1. 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.

  1. 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.

  1. 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;
    }

这篇关于锥体碰撞的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆