C++ Mouse picking for ray to plane intersection

Question

I'm making a C++/OpenGL application (decidedly without GLM), and I'd like to understand how to create a ray from the camera eye to the mouse pointer, as well as detect an intersection point from the resulting ray onto a flat plane.

I have: - Camera position, camera pitch/yaw, camera's Fwd vector - mouse coordinates in pixel measurements (which I can convert to -1:+1 screen space) - A world system where XY are the horizontal and Z is up - A flat plane where Z=0 - My own matrix class which is basically just a float[16] and accessed (Get()) as a float* - Basic matrix multiply function, cross, dot product functions.

I'm still trying to learn 3D math (this project is somewhat of a reason for that), so I'm still somewhat hazy on how to accomplish this. Most tutorials I can find are Unity tutorials, or involve using an Unproject function from glu or GLM.

I'd hopefully like to understand how to accomplish this manually as I'd like to understand what's happening. Thank you for any wisdom or guidance.

Answer 1

Ok: Math hat on:

In computer graphics/games, you typically work in a couple of frames of reference:

1. local/object space
2. world space
3. view space
4. screen space

Local space is represented by relative to the origin of whatever 3d mesh you want to draw.

World space is relative to the game world, i.e. the mesh translated to a position in the world, and rotated to an orientation within the world.

View space is relative to the camera, i.e is it behind, in front of, above, below etc the camera, and by how much. Think of it this way: I Tell you to meet me at the entrance of the Empire State building in New York. That is the world position. From my perspective (view space), it is simply two metres in front of me.

Screen space is when everything is transposed onto the 2d screen and scaled down for display.

So, to get from local space to screen space we have to apply what are known as affine transformations, using 4x4 matrices.

The basic construction of a transform matrix is as follows:

[ R    U     F    T]

[ XR   XU    XF   X]
[ YR   YU    YF   Y]
[ ZR   ZU    ZF   Z]
[ SX   SY    SZ   0]

Ok, so some explanation of the columns is necessary: R stands for "Right", and represents one principle axis of the frame of reference. U stands for "Up", F for "Forward", and T for "Translate". Never mind about S, because it will just muddy the waters at this stage.

Right, up and forward are perpendicular axes, but also represent a change in orientation which will occur when applied to a point, or vector (important later!). Translation just means how much to move the point or vector.

So, now that we have this information, we can now start moving stuff around.

local space is usually defined thusly:

[ 1   0   0  0]
[ 0   1   0  0]
[ 0   0   1  0]
[ 1   1   1  0]

And is known as the "identity matrix"

Basically, if we multiply a vector by this matrix, it will be unchanged.

So, we construct a world matrix, based on the desired position and orientation in the world. If you want to know how to do this, I would recommend reading up on GLM, or whatever math library you want to use for specifics.

Then we construct a "view matrix", which represents the position, direction, and orientation of the camera, or eye. Again, math library specifics can vary, so choose a library and go with the API.

Finally, you need a projection matrix. This is usually comprised of resolution based information etc. Again, API's vary.

Now, to get a point from local space into screen space, we multiply the vector by each of the matrices in turn:

transformed_point = world_transform * point;
transformed_point = view_transform * transformed_point;
transformed_point = projection * transformed_point;

Or: transformed_point = projection * view_transform * world_transform * point;

No, to your question:

in order to transform a 2d point, into a 3d point, you must perform the transform in reverse, using the inverse of the matrices:

world_point = inverse(view) * inverse(projection) * mouse_point;

Then, you can draw a line between the camera position, and the normalised(length of one) mouse_point vector * ray length.

Some resources:

http://antongerdelan.net/opengl/raycasting.html

http://schabby.de/picking-opengl-ray-tracing/

http://www.opengl-tutorial.org/

I recommend the third link, and in particular, read the maths stuff. You'll need it.

Answer 2

I'll just write real code. It works for me. Modern OpenGL.

create a ray from the camera eye to the mouse pointer

void screenToClient(const vec2f& screen_coord, const vec4f& client_rect, vec2f& out )
{
    out.x = screen_coord.x - client_rect.x;
    out.y = screen_coord.y - client_rect.y;
}
void getRay( Ray& ray, const vec2f& coord, const vec4f& viewport_rect,  const Matrix4& ViewProjection_invert )
{
    vec2f point;
    screenToClient(coord,viewport_rect,point);
    
    vec2f rect_size;
    rect_size.x = viewport_rect.z - viewport_rect.x;
    rect_size.y = viewport_rect.w - viewport_rect.y;

    float pt_x = (point.x / rect_size.x) * 2.f - 1.f;
    float pt_y = - (point.y / rect_size.y) * 2.f + 1.f;

    //                                           0.f - for d3d
    ray.m_origin = math::mul(Vector4(pt_x,pt_y,-1.f,1.f), ViewProjection_invert);
    ray.m_end   = math::mul(Vector4(pt_x,pt_y,1.f,1.f), ViewProjection_invert);

    ray.m_origin.w = 1.0f / ray.m_origin.w;
    ray.m_origin.x *= ray.m_origin.w;
    ray.m_origin.y *= ray.m_origin.w;
    ray.m_origin.z *= ray.m_origin.w;

    ray.m_end.w = 1.0f / ray.m_end.w;
    ray.m_end.x *= ray.m_end.w;
    ray.m_end.y *= ray.m_end.w;
    ray.m_end.z *= ray.m_end.w;
}

Matrix inversion is taken from: https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/matrix-inverse

//Vector-matrix product ... not matrix-vector
inline Vector4 mul( const Vector4& vec, const Matrix4& mat )
{
    return Vector4
    (
        mat[ 0u ].x * vec.x + mat[ 1u ].x * vec.y + mat[ 2u ].x * vec.z + mat[ 3u ].x * vec.w,
        mat[ 0u ].y * vec.x + mat[ 1u ].y * vec.y + mat[ 2u ].y * vec.z + mat[ 3u ].y * vec.w,
        mat[ 0u ].z * vec.x + mat[ 1u ].z * vec.y + mat[ 2u ].z * vec.z + mat[ 3u ].z * vec.w,
        mat[ 0u ].w * vec.x + mat[ 1u ].w * vec.y + mat[ 2u ].w * vec.z + mat[ 3u ].w * vec.w
        );
}

Calculating ray direction and other things...

void Ray::update()
{
    m_direction.x = m_end.x - m_origin.x;
    m_direction.y = m_end.y - m_origin.y;
    m_direction.z = m_end.z - m_origin.z;
    m_direction.normalize();

    m_invDir.x = 1.f / m_direction.x;
    m_invDir.y = 1.f / m_direction.y;
    m_invDir.z = 1.f / m_direction.z;
    m_invDir.w = 1.f / m_direction.w;
...

Ray-triangle intersection. You can modify this just for ray-plane intersection if you need.

struct Triangle
{
    Vector4 v1; // vertex 1
    Vector4 v2;
    Vector4 v3;
    Vector4 e1; // edge
    Vector4 e2;

    void update()
    {
        e1 = Vector4( v2.x - v1.x,
            v2.y - v1.y,
            v2.z - v1.z,
            0.f);
        e2 = Vector4( v3.x - v1.x,
            v3.y - v1.y,
            v3.z - v1.z,
            0.f);
    }

    // Möller–Trumbore ray-triangle intersection algorithm
    bool rayTest_MT( const Ray& ray, bool withBackFace, float& len, float& U, float& V, float& W )
    {
        Vector4  pvec = ray.m_direction.cross(e2);
        float det  = e1.dot(pvec);
        
        if( withBackFace )
        {
// Epsilon is #define Epsilon std::numeric_limits<float>::epsilon()
// or 1.192092896e-07F
            if( std::fabs(det) < Epsilon )
                return false;
        }
        else
        {
            if( det < Epsilon && det > -Epsilon )
                return false;
        }

        Vector4 tvec(
            ray.m_origin.x - v1.x,
            ray.m_origin.y - v1.y,
            ray.m_origin.z - v1.z,
            0.f);

        float inv_det = 1.f / det;
        U = tvec.dot(pvec) * inv_det;

        if( U < 0.f || U > 1.f )
            return false;

        Vector4  qvec = tvec.cross(e1);
        V    = ray.m_direction.dot(qvec) * inv_det;

        if( V < 0.f || U + V > 1.f )
            return false;

        len = e2.dot(qvec) * inv_det;
        
        if( len < Epsilon ) return false;
        
        W = 1.f - U - V;
        return true;
    }