How to handle incorrect index calculation for discretized ray tracing?

Question

The situation si as follows. I am trying to implement a linear voxel search in a glsl shader for efficient voxel ray tracing. In toehr words, I have a 3D texture and I am ray tracing on it but I am trying to ray trace such that I only ever check voxels intersected by the ray once.

To this effect I have written a program with the following results:

Not efficient but correct:

The above image was obtained by adding a small epsilon ray multiple times and sampling from the texture on each iteration. Which produces the correct results but it's very inefficient.

That would look like:

loop{
     start += direction*0.01;
     sample(start);
}

To make it efficient I decided to instead implement the following lookup function:

float bound(float val)
{
    if(val >= 0)
        return voxel_size;
    return 0;
}

float planeIntersection(vec3 ray, vec3 origin, vec3 n, vec3 q)
{
    n = normalize(n);
    if(dot(ray,n)!=0)
        return (dot(q,n)-dot(n,origin))/dot(ray,n);

    return -1;
}

vec3 get_voxel(vec3 start, vec3 direction)
{
    direction = normalize(direction);

    vec3 discretized_pos = ivec3((start*1.f/(voxel_size))) * voxel_size;

    vec3 n_x = vec3(sign(direction.x), 0,0);
    vec3 n_y = vec3(0, sign(direction.y),0);    
    vec3 n_z = vec3(0, 0,sign(direction.z));

    float bound_x, bound_y, bound_z;

    bound_x = bound(direction.x);
    bound_y = bound(direction.y);
    bound_z = bound(direction.z);

    float t_x, t_y, t_z;

    t_x = planeIntersection(direction, start, n_x, 
        discretized_pos+vec3(bound_x,0,0));

    t_y = planeIntersection(direction, start, n_y, 
        discretized_pos+vec3(0,bound_y,0));

    t_z = planeIntersection(direction, start, n_z, 
        discretized_pos+vec3(0,0,bound_z));

    if(t_x < 0)
        t_x = 1.f/0.f;
    if(t_y < 0)
        t_y = 1.f/0.f;
    if(t_z < 0)
        t_z = 1.f/0.f;

    float t = min(t_x, t_y);
    t = min(t, t_z);

    return start + direction*t;
}

Which produces the following result:

Notice the triangle aliasing on the left side of some surfaces.

It seems this aliasing occurs because some coordinates are not being set to their correct voxel.

For example modifying the truncation part as follows:

vec3 discretized_pos = ivec3((start*1.f/(voxel_size)) - vec3(0.1)) * voxel_size;

Creates:

So it has fixed the issue for some surfaces and caused it for others.

I wanted to know if there is a way in which I can correct this truncation so that this error does not happen.

Update:

I have narrowed down the issue a bit. Observe the following image:

The numbers represent the order in which I expect the boxes to be visited.

As you can see for some of the points the sampling of the fifth box seems to be ommitted.

The following is the sampling code:

vec4 grabVoxel(vec3 pos)
{

    pos *= 1.f/base_voxel_size;

    pos.x /= (width-1);
    pos.y /= (depth-1);
    pos.z /= (height-1);
    vec4 voxelVal = texture(voxel_map, pos);

    return voxelVal;
}

Answer 1

yep that was the +/- rounding I was talking about in my comments somewhere in your previous questions related to this. What you need to do is having step equal to grid size in one of the axises (and test 3 times once for |dx|=1 then for |dy|=1 and lastly |dz|=1).

Also you should create a debug draw 2D slice through your map to actually see where the hits for a single specific test ray occurred. Now based on direction of ray in each axis you set the rounding rules separately. Without this you are just blindly patching one case and corrupting other two ...

Now actually Look at this (I linked it to your before but you clearly did not):

• Wolf and Doom ray casting techniques

especially pay attention to:

On the right It shows you how to compute the ray step (your epsilon). You simply scale the ray direction so one of the coordinate is +/-1. For simplicity start with 2D slice through your map. The red dot is ray start position. Green is ray step vector for vertical grid lines hits and red is for horizontal grid lines hits (z will be analogically the same).

Now you should add the 2D overview of your map through some height slice that is visible (like on the image on the left) add a dot or marker to each intersection detected but distinguish between x,y and z hits by color. Do this for single ray only (I use the center of view ray). Fist handle view when you look at X+ directions than X- and when done move to Y,Z ...

In my GLSL volumetric 3D back raytracer I also linked you before look at these lines:

if (dir.x<0.0) { p+=dir*(((floor(p.x*n)-_zero)*_n)-ray_pos.x)/dir.x; nnor=vec3(+1.0,0.0,0.0); }
if (dir.x>0.0) { p+=dir*((( ceil(p.x*n)+_zero)*_n)-ray_pos.x)/dir.x; nnor=vec3(-1.0,0.0,0.0); }

if (dir.y<0.0) { p+=dir*(((floor(p.y*n)-_zero)*_n)-ray_pos.y)/dir.y; nnor=vec3(0.0,+1.0,0.0); }
if (dir.y>0.0) { p+=dir*((( ceil(p.y*n)+_zero)*_n)-ray_pos.y)/dir.y; nnor=vec3(0.0,-1.0,0.0); }

if (dir.z<0.0) { p+=dir*(((floor(p.z*n)-_zero)*_n)-ray_pos.z)/dir.z; nnor=vec3(0.0,0.0,+1.0); }
if (dir.z>0.0) { p+=dir*((( ceil(p.z*n)+_zero)*_n)-ray_pos.z)/dir.z; nnor=vec3(0.0,0.0,-1.0); }

they are how I did this. As you can see I use different rounding/flooring rule for each of the 6 cases. This way you handle case without corrupting the other. The rounding rule depends on a lot of stuff like how is your coordinate system offseted to (0,0,0) and more so it might be different in your code but the if conditions should be the same. Also as you can see I am handling this by offsetting the ray start position a bit instead of having these conditions inside the ray traversal loop castray.

That macro cast ray and look for intersections with grid and on top of that actually zsorts the intersections and use the first valid one (that is what l,ll are for and no other conditions or combination of ray results are needed). So my way of deal with this is cast ray for each type of intersection (x,y,z) starting on the first intersection with the grid for the same axis. You need to take into account the starting offset so the l,ll resembles the intersection distance to real start of ray not to offseted one ...

Also a good idea is to do this on CPU side first and when 100% working port to GLSL as in GLSL is very hard to debug things like this.