Ray intersection with GLSL

I’m looking for to simulate the reflection effect using ray tracing on GLSL, however I could not find good references, examples or tutorial related to this topic. When I got some interesting data, the method is limited for specific object’s surfaces (e.g. sphere, box, cylinder…) and it needs to know object’s position and radius; it is not my case. I also know the GLSL does not support recursive functions, but as far as I know the ray tracing can be done iteratively.

My goal is to simulate the reverberation process for an acoustic sensor as follows: primary reflections by rasterization; and secondary reflections by ray tracing. When a ray hits the object’s surface, the distance and normal values are measured.

My scene with all objects is modeled using OpenSceneGraph, and I have used GLSL to handle the normal and position data from the viewpoint on GPU.

Follows below my current GLSL code. At this moment, I am able to calculate the ray parameters (world position and direction vector values, for each pixel), however I do not know how to calculate the data when a ray hits a surface.

Thanks in advance. Any help is very much welcome.

Vertex shader:

#version 130

uniform mat4 osg_ViewMatrixInverse;

out vec3 positionEyeSpace;
out vec3 normalEyeSpace;
uniform vec3 cameraPos;

// ray definition, with an origin point and a direction vector
struct Ray {
    vec3 origin;
    vec3 direction;
};

void main() {
    gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;

    // world space
    mat4 modelWorld = osg_ViewMatrixInverse * gl_ModelViewMatrix;
    vec3 positionWorldSpace = vec3(modelWorld * gl_Vertex);
    vec3 normalWorldSpace = mat3(modelWorld) * gl_Normal;

    // eye space
    positionEyeSpace = vec3(gl_ModelViewMatrix * gl_Vertex);
    normalEyeSpace = gl_NormalMatrix * gl_Normal;

    // calculate the reflection direction for an incident vector
    vec3 I = normalize(positionWorldSpace - cameraPos);
    vec3 N = normalize(normalWorldSpace);
    vec3 reflectedDirection = normalize(reflect(I, N));
}

Fragment shader:

#version 130

in vec3 positionEyeSpace;
in vec3 normalEyeSpace;

uniform float farPlane;
uniform bool drawNormal;
uniform bool drawDepth;

out vec4 out_data;

void main() {
    vec3 nNormalEyeSpace = normalize(normalEyeSpace);

    vec3 nPositionEyeSpace = normalize(-positionEyeSpace);

    float linearDepth = sqrt(positionEyeSpace.x * positionEyeSpace.x +
                             positionEyeSpace.y * positionEyeSpace.y +
                             positionEyeSpace.z * positionEyeSpace.z);

    linearDepth = linearDepth / farPlane;

    // output the normal and depth data as matrix
    out_data = vec4(0, 0, 0, 1);
    if (linearDepth <= 1) {
        if (drawNormal) out_data.z = abs(dot(nPositionEyeSpace, nNormalEyeSpace));
        if (drawDepth)  out_data.y = linearDepth;
    }

    gl_FragDepth = linearDepth;
}

The first surface is easy; as you note, you can just rasterise the surface, calculating the incident vector by subtracting the eye position from the fragment position, then the reflected vector from the incident vector and surface normal. At that point, you have a ray defined by starting position and direction.

That’s where it gets complicated: given a ray, you need to find the first surface which that ray intersects. This problem is central to any form of raytracing, whether on the CPU or GPU. If the number of surfaces is small, you can just test each one. Otherwise, you need some form of spatial index (3D array, octree, BSP tree, bounding volume hierarchy etc). Most techniques can be converted to use from GLSL, provided that they don’t require specific libraries.

So I would suggest just looking at articles on raytracing in general, finding one which fits your particular scene structure and use case, then figuring out how to implement that in GLSL.

Hi GClements,

thanks for the quick reply. I agree the first surface is easy by rasterization, and the ray properties can be calculated on reflected surface (starting position and direction). My idea is to simulate a single ray for each reflected point, this approach will save computational time. A priori, I have no information about the objects that compose the 3D scene (position, shape…).

This spatial index of each object present in the scene is a good suggestion. I will also look in this direction.

By the way, I still have a doubt: If I use ray tracing for the first surface, and I know a priori the objects’ positions, how can I calculate the intersection between the ray and any object surface?

Thanks a lot.

GClements,

do you have any example with RayTracing and GLSL?

That depends upon the nature of the surface. For a triangle, you’d treat the vertex positions as a 3x3 matrix which transforms barycentric coordinates to spatial coordinates:


[x]   [x1 x2 x3] [a]
[y] = [y1 y2 y3].[b]
[z]   [z1 z2 z3] [c]

Inverting the matrix gives a matrix M which transforms spatial coordinates to barycentric coordinates.


[a]     [x]
[b] = M.[y]
[c]     [z]

Note that this only needs to be done once, when the geometry is defined, not every frame.

For each ray p=s+td (in world space), transform the start position and direction of the ray into barycentric coordinates:
b=M.s+t
M.d

Finding the intersection with the triangle’s plane is just solving a+b+c=1 for t, i.e. u+tv=1 => t=(1-u)/v where u and v are just the sums of the coordinates of M.s and M.d respectively. Note that if v=0 then the ray is parallel to the plane of the triangle and there is no intersection. Substituting t back into [a b c]T=M.s+tM.d gives the barycentric coordinates of the intersection point; if all three are positive, the intersection point lies inside the triangle, otherwise it’s outside.

Only this example of ray-sphere intersection. But it’s only the first surface; it doesn’t handle reflection. And it’s only one object; the spatial index (so you aren’t testing every ray against every surface) tends to be the hard part.

Hi GClements,

I’m going to implement this way to store all objects as triangles and pass this information to shader. I will update here with the progresses.

Hi GClements,

I am able to store all vertices of 3D models as triangles and pass them to shader as texture, and then perform the ray-triangle intersection using classic algorithms (e.g. Möller–Trumbore ray-triangle intersection).

My idea was to distribute these triangles into a k-d tree, and I got excellent query performance results using C++ (it is a really good approach). However, my problem comes when I tried to pass the Nearest Neighbor Searching to GLSL, once shaders doesn’t allow recursive calls. Do you know how can I figure out this point?

My current code in C++ is:


// Calculate the nearest node for a specific value (using node)
void nearest(KDnode *root, KDnode *nd, int i, int dim, KDnode **best, double *best_dist, int *visited)
{
    if (!root)
        return;

    double d = dist(root, nd);
    double dx = root->data[3][i] - nd->data[3][i];
    double dx2 = dx * dx;

    (*visited)++;

    if (!*best || d < *best_dist)
    {
        *best_dist = d;
        *best = root;
    }

    // If chance of exact match is high
    if (!*best_dist)
        return;

    i = ++i % dim;

    nearest(dx > 0 ? root->left : root->right, nd, i, dim, best, best_dist, visited);

    if (dx2 < *best_dist) {
        nearest1(dx > 0 ? root->right : root->left, nd, i, dim, best, best_dist, visited);
    }
}


Thanks in advance,

Rômulo.

The usual alternative to recursion in situations where it can’t be used is a trampoline. Essentially, you call the function repeatedly in a loop. Rather than taking parameters, the function retrieves its parameters from the top of a stack. Instead of making recursive calls, the function pushes the parameters onto the stack and returns, so the next iteration of the loop makes the recursive call(s). The loop terminates when the stack is empty.

The main thing to bear in mind is that a GPU gets much of its performance from vector parallelism (SIMD), so you need to take this into account to get decent performance. One thing which is bad for parallelism is variable-length loops, as the number of iterations required is the maximum of all threads in the group. If you can’t figure out how to make the algorithm SIMD-friendly, you may be better off running it on the CPU.

@GClements, thanks for all your support until now!

Could you kindly help me with the reflection part on GLSL? For each fragment, I would like to compute its world position and world surface normal, and reflect it over the surface normal. I will use these information to simulate the ray from each fragment with the ray’s origin at the world position, and the direction I calculated earlier using surface normal.

Thanks in advance.

Vertex:

#version 130
 
uniform mat4 osg_ViewMatrixInverse;
 
out vec3 positionEyeSpace;
out vec3 normalEyeSpace;
uniform vec3 cameraPos;
 
void main() {
    gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;
 
    // world space
    mat4 modelWorld = osg_ViewMatrixInverse * gl_ModelViewMatrix;
    vec3 positionWorldSpace = vec3(modelWorld * gl_Vertex);
    vec3 normalWorldSpace = mat3(modelWorld) * gl_Normal;
 
    // eye space
    positionEyeSpace = vec3(gl_ModelViewMatrix * gl_Vertex);
    normalEyeSpace = gl_NormalMatrix * gl_Normal;
 
    // calculate the reflection direction for an incident vector
    vec3 I = normalize(positionWorldSpace - cameraPos);
    vec3 N = normalize(normalWorldSpace);
    vec3 reflectedDirection = normalize(reflect(I, N));
}

Just make positionWorldSpace and normalWorldSpace vertex shader outputs, i.e.


out vec3 positionWorldSpace;
out vec3 normalWorldSpace;

The fragment shader needs inputs with the same names, which will be the interpolated values. You can then move the reflection computation (the last 3 lines of the vertex shader) to the fragment shader.

Note that the ray direction doesn’t need to be normalised. reflect() requires the normal to be normalised, but not the incidence vector; the reflected vector will have the same magnitude as the incidence vector. Intersection calculations shouldn’t require a normalised direction. Also, if a 3x3 matrix is orthonormal (only rotation, no scaling), it will preserve lengths, so you don’t need to re-normalise after transformation.

You probably don’t need anything in eye space. Eye space is used in the fixed-function pipeline and conventional forward-rendering because lighting calculations work in any coordinate system that’s affine to world space, it avoids the need to subtract the view position (which is always the origin in eye space), and you need to use the view transformation at some point to get to clip space. But if you’re doing raytracing with a spatial index that uses world space, you may as well just use clip space for gl_Position and world space for everything else, and forget about eye space.

Hi GClements,

thank you so much. I have followed your suggestion and moved the reflection computation to fragment, as well as I converted the position and normal in the world space as inputs.

Since I’ve tried to simulate the first reflections by rasterization and second ones by ray-tracing, so I need some data in eye space for the first process. In this case, can I avoid the camera position?

Follows my current code. Could you have a look again, please?

Vertex:

#version 130

uniform mat4 osg_ViewMatrixInverse;

out vec3 positionEyeSpace;
out vec3 normalEyeSpace;
out vec3 positionWorldSpace;
out vec3 normalWorldSpace;

void main() {
    gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;

    // eye space
    positionEyeSpace = vec3(gl_ModelViewMatrix * gl_Vertex);
    normalEyeSpace = gl_NormalMatrix * gl_Normal;

    // world space
    mat4 modelWorld = osg_ViewMatrixInverse * gl_ModelViewMatrix;
    positionWorldSpace = vec3(modelWorld * gl_Vertex);
    normalWorldSpace = mat3(modelWorld) * gl_Normal;

    // Texture for normal mapping (irregularities surfaces)
    gl_TexCoord[0] = gl_TextureMatrix[0] * gl_MultiTexCoord0;
}

Fragment:

#version 130

in vec3 positionEyeSpace;
in vec3 normalEyeSpace;
in vec3 positionWorldSpace;
in vec3 normalWorldSpace;

uniform float farPlane;

// ray definition
struct Ray {
    vec3 origin;        // starting point
    vec3 direction;     // ray direction
};


// primary reflections: rasterization
vec4 primaryReflections() {
    vec3 nNormalEyeSpace = normalize(normalEyeSpace);

    // Depth calculation
    vec3 nPositionEyeSpace = normalize(-positionEyeSpace);
    float linearDepth = length(positionEyeSpace);

    // Normalize depth using range value (farPlane)
    linearDepth = linearDepth / farPlane;
    gl_FragDepth = linearDepth;

    // presents the normal and depth data as matrix
    vec4 output = vec4(0, 0, 0, 1);
    if (linearDepth <= 1) {
        output.y = linearDepth;
        output.z = abs(dot(nPositionEyeSpace, nNormalEyeSpace));
    }

    return output;
}

// secondary reflections: ray-triangle intersection
vec4 secondaryReflections(vec4 firstR) {
    // calculate the reflection direction for an incident vector
    vec3 nNormalWorldSpace = normalize(normalWorldSpace);
    vec3 reflectedDir = reflect(positionWorldSpace, nNormalWorldSpace);

    Ray ray;
    ray.origin = positionWorld;
    ray.direction = reflectedDir;

    // TODO: perform the ray-triangle intersection
    vec4 output = vec4(0,0,0,1);

    return output;
}

void main() {
    // output: primary reflections by rasterization
    vec4 firstR = primaryReflections();
    vec4 secndR = secondaryReflections(firstR);

    // gl_FragData[0] = firstR;
    gl_FragData[0] = secndR;
}

This is the distance from the viewpoint, which isn’t the same thing as depth.

Why are you writing out linear depth (or distance)? Writing to gl_FragCoord can be relatively expensive, as it disables early depth tests.

Note that using world-space vectors would compute the same value (but note that nPositionEyeSpace is actually a direction rather than a position, so you’d need to subtract from the view position rather than simply negating; but that’s a trivial operation).

The incident vector should be cameraPos-positionWorldSpace. In eye space, the camera position is (0,0,0), so the calculation simplifies to -positionEyeSpace. But that doesn’t apply to world space.

In short, if you don’t need the eye-space position or normal, you shouldn’t calculate them. The more outputs the vertex shader has, the fewer vertices will fit in the cache. Fixed-function lighting uses eye space instead of world space, but raytracing is (presumably) going to need to use world-space rays for the intersection calculations (either that, or you’d need to transform your scene description into eye space each frame). If possible, you should only use world space

Hi GClement,

thanks for your kindly explanation and support. I followed your tips and now all calculations are done on world-space.

My current question: can I associate incident vector to Ray.origin?

Thanks in advance.

Vertex:

#version 130

uniform mat4 osg_ViewMatrixInverse;

out vec3 posWorldSpace;
out vec3 normalWorldSpace;
out vec3 cameraPos;

void main() {
    gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;

    // position/normal in world space
    mat4 modelWorld = osg_ViewMatrixInverse * gl_ModelViewMatrix;
    posWorldSpace = vec3(modelWorld * gl_Vertex);
    normalWorldSpace = mat3(modelWorld) * gl_Normal;

    // incident vector
    cameraPos = osg_ViewMatrixInverse[3].xyz;

    // Texture for normal mapping (irregularities surfaces)
    gl_TexCoord[0] = gl_TextureMatrix[0] * gl_MultiTexCoord0;
}

Fragment:

#version 130

in vec3 posWorldSpace;
in vec3 normalWorldSpace;
in vec3 cameraPos;

uniform float farPlane;

// ray definition
struct Ray {
    vec3 origin;        // starting point
    vec3 direction;     // ray direction
};

// primary reflections: rasterization
vec4 primaryReflections() {
    vec3 nPosWorldSpace = normalize(cameraPos - posWorldSpace);
    vec3 nNormalWorldSpace = normalize(normalWorldSpace);

    // Distance calculation
    float viewDistance = length(cameraPos - posWorldSpace);

    // Normalize distance using range value (farPlane)
    float nViewDistance = viewDistance / farPlane;

    // presents the normal and depth data as matrix
    vec4 output = vec4(0, 0, 0, 1);
    if (nViewDistance <= 1) {
        output.y = nViewDistance;
        output.z = abs(dot(nPosWorldSpace, nNormalWorldSpace));
    }

    return output;
}

// secondary reflections: ray-triangle intersection
vec4 secondaryReflections(vec4 firstR) {

    // calculate the reflection direction for an incident vector
    vec3 incidentVec = cameraPos - posWorldSpace;
    vec3 nNormalWorldSpace = normalize(normalWorldSpace);
    vec3 reflectedDir = reflect(incidentVec, nNormalWorldSpace);

    // setup ray
    Ray ray;
    ray.origin = incidentVec;
    ray.direction = reflectedDir;

    // TODO: perform ray-triangle intersection
    vec4 output = vec4(0,0,0,1);
    return output;
}

void main() {
    // output: primary reflections by rasterization
    vec4 firstR = primaryReflections();

    // output: secondary reflections by ray-tracing
    vec4 secndR = secondaryReflections(firstR);

    // gl_FragData[0] = firstR;
    gl_FragData[0] = secndR;
}

The first ray has origin [var]cameraPos[/var], direction [var]posWorldSpace-cameraPos[/var]. The second ray has origin [var]posWorldSpace[/var], direction [var]reflect(posWorldSpace-cameraPos, normalWorldSpace)[/var].

The “eye” vector used for lighting calculations is the negation of the incident vector, but it’s probably cheaper to negate the result of the dot product rather than the vector: dot(-a,b)=-dot(a,b).

Thanks, GClement.

The reflect requires the normalized normal vector, right? i.e. reflect(posWorldSpace - cameraPos, normalize(normalWorldSpace));

Yes.

Although if you aren’t going to be normalising for other reasons, it might be more efficient to use e.g.


vec3 reflect1(vec3 I, vec3 N)
{
    return I - (2.0 * dot(N, I) / dot(N,N)) * N;
}

which doesn’t require N to be normalised.

Thanks again, GClement. I will do some tests and update you here. :slight_smile:

Hi GClement,

I have gotten some progress with my challenge. At this moment, I am trying to get secondary reflection with ray-tracing without success - I didn’t get any intersection. Could you kindly have a look in my current fragment shader?

Given:
[ul]
[li] Use the primary reflections (with rasterization) as input of secondary reflections (with ray tracing);[/li][li] Triangles data (vertices, centroid and normal in world coordinates) are sorted as level order k-d tree and sent to shader as texture;[/li][li] The nearest function calculates the closest triangle to the current ray, performing a nearest neighbor search;[/li][li] Use of Möller–Trumbore ray-triangle intersection algorithm;[/li][/ul]

Thanks in advance.

#version 130

in vec3 worldPos;
in vec3 worldNormal;
in vec3 cameraPos;

uniform float farPlane;
uniform sampler2D trianglesTex;         // all triangles/meshes collected from the simulated scene
uniform vec2 trianglesTexSize;          // texture size of triangles

// ray definition
struct Ray {
    vec3 origin, direction;
};

// triangle definition
struct Triangle {
    vec3 v0;        // vertex A
    vec3 v1;        // vertex B
    vec3 v2;        // vertex C
    vec3 center;    // centroid
    vec3 normal;    // normal
};

...

float getTexData(int i, int j) {
    return texelFetch(trianglesTex, ivec2(i,j), 0).r;
}

Triangle getTriangleData(int idx) {
    Triangle triangle;
    triangle.v0     = vec3(getTexData(idx,0), getTexData(idx,1), getTexData(idx,2));
    triangle.v1     = vec3(getTexData(idx,3), getTexData(idx,4), getTexData(idx,5));
    triangle.v2     = vec3(getTexData(idx,6), getTexData(idx,7), getTexData(idx,8));
    triangle.center = vec3(getTexData(idx,9), getTexData(idx,10), getTexData(idx,11));
    triangle.normal = vec3(getTexData(idx,12), getTexData(idx,13), getTexData(idx,14));
    return triangle;
}

// Calculate the distance between two points
float dist(vec3 a, vec3 b) {
    vec3 c = a - b;
    return (c.x * c.x +
            c.y * c.y +
            c.z * c.z);
}

// Calculate the closest triangle to the ray
Triangle nearest(vec3 nd) {

    Queue queue = createQueue();
    enqueue(queue, 0);
    int i = 0;

    bool first = true;
    Triangle best;
    float best_dist;

    while (!isEmpty(queue)) {
        int idx = dequeue(queue);

        if (idx > (trianglesTexSize.x - 1))
            continue;

        Triangle root = getTriangleData(idx);
        float d = dist(root.center, nd);
        float dx = 0;
        switch (i) {
            case 0: dx = root.center.x - nd.x; break;
            case 1: dx = root.center.y - nd.y; break;
            case 2: dx = root.center.z - nd.z; break;
        }

        if (first || d < best_dist) {
            best_dist = d;
            best = root;
            first = false;
        }

        // If chance of exact match is high
        if (best_dist == 0)
            return best;

        i = (i + 1) % 3;

        int idx_left  = 2 * idx + 1;
        int idx_right = 2 * idx + 2;
        if (dx > 0)
            enqueue(queue, idx_left);
        else
            enqueue(queue, idx_right);
    }
    return best;
}

// Möller–Trumbore ray-triangle intersection algorithm
// source: http://bit.ly/2MxnPMG
bool rayIntersectsTriangle(Ray ray, Triangle triangle)
{
    float EPSILON = 0.0000001;

    vec3 edge1, edge2, h, s, q;
    float a, f, u, v;
    edge1 = triangle.v1 - triangle.v0;
    edge2 = triangle.v2 - triangle.v0;

    h = cross(ray.direction, edge2);
    a = dot(edge1, h);
    if (a > -EPSILON && a < EPSILON)
        return false;

    f = 1 / a;
    s = ray.origin - triangle.v0;
    u = f * dot(s, h);
    if (u < 0.0 || u > 1.0)
        return false;

    q = cross(s, edge1);
    v = f * dot(ray.direction, q);
    if (v < 0.0 || (u + v) > 1.0)
        return false;

    // At this stage we can compute t to find out where the intersection point is on the line.
    float t = f * dot(edge2, q);
    if (t <= EPSILON)   // this means that there is a line intersection but not a ray intersection.
        return false;

    return true;        // ray intersection
}


// ============================================================================================================================

// primary reflections: rasterization
vec4 primaryReflections() {
    vec3 worldIncident = cameraPos - worldPos;
    vec3 nWorldPos = normalize(worldIncident);
    vec3 nWorldNormal = normalize(worldNormal);

    // Distance calculation
    float viewDistance = length(worldIncident);

    // Normalize distance using range value (farPlane)
    float nViewDistance = viewDistance / farPlane;

    // presents the normal and depth data as matrix
    vec4 output = vec4(0, 0, 0, 1);
    if (nViewDistance <= 1) {
        output.y = nViewDistance;
        output.z = abs(dot(nWorldPos, nWorldNormal));
    }

    return output;
}

// ============================================================================================================================

// secondary reflections: ray-triangle intersection
vec4 secondaryReflections(vec4 firstR) {

    // calculate the reflection direction for an incident vector
    vec3 worldIncident = cameraPos - worldPos;
    vec3 nWorldNormal = normalize(worldNormal);
    vec3 reflectedDir = reflect(-worldIncident, nWorldNormal);

    // set current ray
    Ray ray = Ray(worldPos, reflectedDir);

    // perform ray-triangle intersection only for pixels with valid normal values
    vec4 output = vec4(0,0,0,1);
    if (firstR.z > 0) {
        // find the closest triangle to the origin point
        Triangle triangle = nearest(ray.origin);

        // ray-triangle intersection
        bool intersected = rayIntersectsTriangle(ray, triangle);
        if (intersected) {
            output = vec4(1,1,1,1);
        }
   }

    return output;
}

// ============================================================================================================================

void main() {
    // output: primary reflections by rasterization
    vec4 firstR = primaryReflections();

    // output: secondary reflections by ray-tracing
    vec4 secndR = secondaryReflections(firstR);

    // gl_FragData[0] = firstR;
    gl_FragData[0] = secndR;
}

0

The first thing I’d do is write out some information regarding the triangle which nearest() finds, and check whether it makes sense.

Hi GClements,

I have good news: I can advance a lot with my adventure with ray tracing! Thanks for all your help!

After I updated my shader code to work in world space coordinates, I faced with some problems with my normal mapping process. Once normal maps are built in tangent space, interpolating the vertex normal and a RGB texture, we need to create a TBN matrix to convert between these two spaces (tangent and world space). I believe my problem is due to work with tangent space.

1) In the vertex shader, I calculate the TBN matrix as follows:

    vec3 N = worldNormal;
    vec3 T = mat3(modelWorld) * gl_MultiTexCoord0.xyz;
    vec3 B = cross(N, T);
    TBN = transpose(mat3(T, B, N));

and I got the following sonar image result. Please look the normal mapping is applied on the center of sonar image, however the corners are black.

[ATTACH=CONFIG]1854[/ATTACH]

2) If I remove the tranpose operation on TBN matrix, I got the following result, with corners filled and the center of sonar image in black:

    TBN = mat3(T, B, N);

[ATTACH=CONFIG]1855[/ATTACH]

If I just modify the bitangent vector by sum of cross(N,T) and cross(T,N), I got a weird result.

    vec3 B = cross(N, T) + cross(T, N);
    TBN = transpose(mat3(T, B, N));

[ATTACH=CONFIG]1856[/ATTACH]

Do you have any tip how can I solve this problem?

Thanks in advance,