c++opengl glm-math perspectivecamera ray-picking

Computing&Plotting 3D Points on Surface of a Mesh from mouse position

I need to plot a set of 3D points (could be curve/fonts) on surface of a mesh such that they can be retrieved later, so rendering points on a texture and attaching the texture does not work. The points will be drawed/edited by mouse clicks. This is for CAD purposes so precision is important.

How can I get 3D vertices in mesh local coordinates from 2D mouse position?

I am using OpenGL for the rendering part with perspective projection created using glm::perspective() with:
```
FOV = 45.0f
aspect ratio = 16 : 9
zNear = 0.1f
zFar = 100.0f
triangulated mesh
```
Is it possible to do Ray-Triangle Intersection calculations in the object space?
Does the camera Position (Origin) have to be World Space?

Solution

what you need is (Assuming old api OpenGL default notations):

perspective projection: FOVx,FOVy,znear
inverse of Model*View matrix
mouse position converted to <-1,+1> range
list of triangles your mesh is composed of

You can do this like this:

cast ray from camera origin
convert it to mesh local coordinates
compute closest ray/triangle intersection

Yes you can compute this in mesh local coordinates and yes in such case you need Ray in the same coordinates.

Here simple C++/old api OpenGL/VCL example:

//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#include "gl_simple.h"
#include "GLSL_math.h"
#pragma hdrstop
#include "Unit1.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
float mx=0.0,my=0.0;    // mouse position
//---------------------------------------------------------------------------
// Icosahedron
#define icoX .525731112119133606
#define icoZ .850650808352039932
const GLfloat vdata[12][3] =
    {
    {-icoX,0.0,icoZ}, {icoX,0.0,icoZ}, {-icoX,0.0,-icoZ}, {icoX,0.0,-icoZ},
    {0.0,icoZ,icoX}, {0.0,icoZ,-icoX}, {0.0,-icoZ,icoX}, {0.0,-icoZ,-icoX},
    {icoZ,icoX,0.0}, {-icoZ,icoX,0.0}, {icoZ,-icoX,0.0}, {-icoZ,-icoX,0.0},
    };
const int tindices=20;
const GLuint tindice[tindices][3] =
    {
    {0,4,1}, {0,9,4}, {9,5,4}, {4,5,8}, {4,8,1},
    {8,10,1}, {8,3,10}, {5,3,8}, {5,2,3}, {2,7,3},
    {7,10,3}, {7,6,10}, {7,11,6}, {11,0,6}, {0,1,6},
    {6,1,10}, {9,0,11}, {9,11,2}, {9,2,5}, {7,2,11}
    };
//---------------------------------------------------------------------------
void icosahedron_draw() // renders mesh using old api
    {
    int i;
    GLfloat nx,ny,nz;
    glEnable(GL_CULL_FACE);
    glFrontFace(GL_CW);
    glBegin(GL_TRIANGLES);
    for (i=0;i<tindices;i++)
        {
        nx =vdata[tindice[i][0]][0];
        ny =vdata[tindice[i][0]][1];
        nz =vdata[tindice[i][0]][2];
        nx+=vdata[tindice[i][1]][0];
        ny+=vdata[tindice[i][1]][1];
        nz+=vdata[tindice[i][1]][2];
        nx+=vdata[tindice[i][2]][0]; nx/=3.0;
        ny+=vdata[tindice[i][2]][1]; ny/=3.0;
        nz+=vdata[tindice[i][2]][2]; nz/=3.0;
        glNormal3f(nx,ny,nz);
        glVertex3fv(vdata[tindice[i][0]]);
        glVertex3fv(vdata[tindice[i][1]]);
        glVertex3fv(vdata[tindice[i][2]]);
        }
    glEnd();
    }
//---------------------------------------------------------------------------
vec3 ray_pick(float mx,float my,mat4 _mv)   // return closest intersection using mouse mx,my <-1,+1> position and inverse of ModelView _mv
    {
    // Perspective settings
    const float deg=M_PI/180.0;
    const float _zero=1e-6;
    float znear=0.1;
    float FOVy=45.0*deg;
    float FOVx=FOVy*xs/ys;  // use aspect ratio if you do not know screen resolution
    // Ray endpoints in camera local coordinates
    vec3 pos=vec3(mx*tan(0.5*FOVx)*znear,my*tan(0.5*FOVy)*znear,-znear);
    vec3 dir=vec3(0.0,0.0,0.0);
    // Transform to mesh local coordinates
    pos=(_mv*vec4(pos,1.0)).xyz;
    dir=(_mv*vec4(dir,1.0)).xyz;
    // convert endpoint to direction
    dir=normalize(pos-dir);
    // needed variables
    vec3 pnt=vec3(0.0,0.0,0.0);
    vec3 v0,v1,v2,e1,e2,n,p,q,r;
    int i,ii=1;
    float t=-1.0,tt=-1.0,u,v,det,idet;
    // loop through all triangles
    for (int i=0;i<tindices;i++)
        {
        // load v0,v1,v2 with actual triangle
        v0.x=vdata[tindice[i][0]][0];
        v0.y=vdata[tindice[i][0]][1];
        v0.z=vdata[tindice[i][0]][2];
        v1.x=vdata[tindice[i][1]][0];
        v1.y=vdata[tindice[i][1]][1];
        v1.z=vdata[tindice[i][1]][2];
        v2.x=vdata[tindice[i][2]][0];
        v2.y=vdata[tindice[i][2]][1];
        v2.z=vdata[tindice[i][2]][2];
        //compute ray(pos,dir) triangle(v0,v1,v2) intersection
        e1=v1-v0;
        e2=v2-v0;
        // Calculate planes normal vector
        p=cross(dir,e2);
        det=dot(e1,p);
        // Ray is parallel to plane
        if (abs(det)<1e-8) continue;
        idet=1.0/det;
        r=pos-v0;
        u=dot(r,p)*idet;
        if ((u<0.0)||(u>1.0)) continue;
        q=cross(r,e1);
        v=dot(dir,q)*idet;
        if ((v<0.0)||(u+v>1.0)) continue;
        t=dot(e2,q)*idet;
        // remember closest intersection to camera
        if ((t>_zero)&&((t<=tt)||(ii!=0)))
            {
            ii=0; tt=t;
            // barycentric interpolate position
            t=1.0-u-v;
            pnt=(v0*t)+(v1*u)+(v2*v);
            }
        }
    return pnt; // if (ii==1) no intersection found
    }
//---------------------------------------------------------------------------
void gl_draw()
    {
    glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
    glDisable(GL_TEXTURE_2D);
    glEnable(GL_DEPTH_TEST);
    glEnable(GL_LIGHT0);
    glEnable(GL_CULL_FACE);
//  glDisable(GL_CULL_FACE);
    glFrontFace(GL_CCW);
    glEnable(GL_COLOR_MATERIAL);
/*
    glPolygonMode(GL_FRONT,GL_FILL);
    glPolygonMode(GL_BACK,GL_LINE);
    glDisable(GL_CULL_FACE);
*/
    // set projection
    glMatrixMode(GL_PROJECTION);        // operacie s projekcnou maticou
    glLoadIdentity();                   // jednotkova matica projekcie
    gluPerspective(45,float(xs)/float(ys),0.1,100.0); // matica=perspektiva,120 stupnov premieta z viewsize do 0.1

    // set view
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslatef(0.2,0.0,-5.0);
    static float ang=0.0;
    glRotatef(ang,0.2,0.7,0.2); ang+=5.0; if (ang>=360.0) ang-=360.0;

    // obtain actual modelview matrix (mv) and its inverse (_mv)
    mat4 mv,_mv;
    float m[16];
    glGetFloatv(GL_MODELVIEW_MATRIX,m);
    mv.set(m);
    _mv=inverse(mv);

    // render mesh
    glColor3f(0.5,0.5,0.5);
    glEnable(GL_LIGHTING);
    icosahedron_draw();
    glDisable(GL_LIGHTING);

    // get point mouse points to
    vec3 p=ray_pick(mx,my,_mv);

    // render it for visual check
    float r=0.1;
    glColor3f(1.0,1.0,0.0);
    glBegin(GL_LINES);
    glVertex3f(p.x-r,p.y,p.z); glVertex3f(p.x+r,p.y,p.z);
    glVertex3f(p.x,p.y-r,p.z); glVertex3f(p.x,p.y+r,p.z);
    glVertex3f(p.x,p.y,p.z-r); glVertex3f(p.x,p.y,p.z+r);
    glEnd();

//  glFlush();
    glFinish();
    SwapBuffers(hdc);
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    // Init of program
    gl_init(Handle);    // init OpenGL
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
    {
    // Exit of program
    gl_exit();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    // repaint
    gl_draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
    {
    // resize
    gl_resize(ClientWidth,ClientHeight);
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::tim_redrawTimer(TObject *Sender)
    {
    gl_draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormMouseMove(TObject *Sender, TShiftState Shift, int X, int Y)
    {
    // just event to obtain actual mouse position
    // and convert it from screen coordinates <0,xs),<0,ys) to <-1,+1> range
    mx=X; mx=(2.0*mx/float(xs-1))-1.0;
    my=Y; my=1.0-(2.0*my/float(ys-1)); // y is mirrored in OpenGL
    }
//---------------------------------------------------------------------------

The only imnportant thing is function ray_pick which returns your 3D point based on mouse 2D position and actual inverse of ModelView matrix...

Here preview:

I render yellow cross at the found 3D position as you can see its in direct contact to surface (as half of its line are below surface).

On top of usual stuff I used mine libs: gl_simple.h for the OpenGL context creation and GLSL_math.h instead of GLM for vector and matrix math. But the GL context can be created anyhow and you can use what you have for the math too (I think even the syntax is the same as GLM as they tried to mimic GLSL too...)

Looks like for aspects 1:1 this works perfectly, and in rectangular aspects its slightly imprecise the further away from screen center you get (most likely because I used raw gluPerspective which has some imprecise terms in it, or I missed some correction while ray creation but I doubt that)

In case you do not need high precision (not your case as CAD/CAM needs as high precision as you can get) You can get rid of the ray mesh intersection and directly pick depth buffer at mouse position and compute the resulting point from that (no need for mesh).

For more info see related QAs: