2

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?

Matthew Cotton
  • 65
  • 5
  • Description of your "cone" is doubtful – MBo Jun 02 '20 at 08:24
  • 1
    [Cone is here](https://godotengine.org/qa/20875/how-to-create-a-cone-collisionshape-3d). Perhaps you mean union of cone and spherical cap - [spherical sector](https://en.wikipedia.org/wiki/Spherical_sector) (just found this term ;)) – MBo Jun 02 '20 at 11:10
  • Ahh spherical cap and sector is what I'm after! This is valuable information. Thanks! – Matthew Cotton Jun 02 '20 at 22:29
  • Does this answer your question? [Detect if a cube and a cone intersect each other?](https://stackoverflow.com/questions/22023977/detect-if-a-cube-and-a-cone-intersect-each-other) – Spektre Jun 03 '20 at 17:57
  • I did see that question previously and I think it's a partial solution. Because my problem involves the rounded base of the cone. The "circular sector" as I have now learned the term. I was hoping for a generic solution but I'm guessing the real solution is a collection of tests involving the circle and the cone/triangle. – Matthew Cotton Jun 03 '20 at 22:28
  • @MatthewCotton I added an answer with alternative solution using lower geometric primitives. Was too lazy to do the analytical solution (but that could be faster and more precise) – Spektre Jun 08 '20 at 08:55

1 Answers1

4

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:

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:

preview

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;
    }
Spektre
  • 49,595
  • 11
  • 110
  • 380
  • 1
    Accepting this as an answer since it more closely answers my original question than what I ultimately ended up doing. I ended up going with a 2D solution since my original problem could be easily reduced down to 2D and it greatly simplified the solution. For points I just made sure it was within a radius and within the angles making up the sides of the arc. And for more complex shapes I broke it down into testing corners, solving intersections and a little bit of logic. – Matthew Cotton Jun 11 '20 at 22:28
  • 1
    I think I found one of the "unhandeled edge cases" in `lin closest (line l0, line l1)`. As I understand it, it computes the closest points of two segments, say AB and CD. If you check the following input `A = {1,1,0}, B = {2,3,0}, C = {2,1,0}, D = {2,-2,0}` it says the closest points are `{2,3,0}` and `{2,1,0}`. See [this](https://godbolt.org/z/4G5sc7e87) example. I'd expected the fist point to be `{1.2, 1.4, 0}`. Shame on me in case I screwed up the implementation of `vec3` and helper functions (I am not using OpenGL). – StefanKssmr Apr 11 '21 at 21:03
  • 1
    @StefanKssmr looks like you're right nice catch +1 will investigate it is most likely some rounding relate error as the direction of resulting line is the same as one of the operands – Spektre Apr 12 '21 at 06:22
  • @Spektre I came to your post from [this](https://stackoverflow.com/q/66979936/8359552) thread. Maybe it will give you some clues. – StefanKssmr Apr 12 '21 at 06:42
  • 1
    @StefanKssmr I already found out what is the problem ... I made this from axis vs axis test so if the closest point is perpendicular distance and intersection between the 2 lines exist if they would be infinite lines then the perpendicular distance is ignored ... trying to repair it now with set of few insane checks ... – Spektre Apr 12 '21 at 06:44
  • @Spektre great, nice answer, though. +1 – StefanKssmr Apr 12 '21 at 06:49
  • @StefanKssmr I think its fixed (at least the visual check look OK to me) however one thing bothers me ... I need to divide the dot product by `2.0f*l1.l` instead of just by `l1.l` (which equations demand) not sure why... the result looks OK even if I change the lines or order of points/lines – Spektre Apr 12 '21 at 07:45
  • 1
    @Spektre I think that instead dividing by `2.0f*l1.l` you should divide by `l1.l*l1.l`, i.e. by the length squared. [Here](https://godbolt.org/z/Eqrqbsso8)'s an example. – StefanKssmr Apr 12 '21 at 08:22
  • @StefanKssmr yes ... `2.0f*` is used for reflect and it was just coincidence it worked ... once I change the angles of tested lines a lot the `2.0*` was wrong and `sqr()` was OK ... I forgot that the `t?` is multiplied by `l?.dp` latter on that is why I must use `sqr(size)` – Spektre Apr 12 '21 at 10:01
  • 1
    @StefanKssmr I finished repairing my 2D and 3D geometry lib ... however I can not post the upgraded code as I made many changes additions in the code since I post this post and the 3D code is ~48KByte which does not fit the SO 30K limit. The 2D version is ~28KByte but still does not fit as there is also text and code formatting would add a lot of too ... – Spektre Apr 12 '21 at 11:37
  • 1
    Thank you so much for this answer, just saved me tons of pain trying out different approaches with spherical cone vs my primitives (OBB, capsule, sphere). Simply calculating closest point between the sphere's center and the other primitive is fairly straightforward, and the overall solution thereafter is verrrry straightforward. Thanks again! – DragonJawad Sep 25 '21 at 07:43