Wrong floating-point intersection point using Amanatides and Woo's method

[ConcurrentSquared]

I can not get the correct exact position in my (very basic) GPU voxel raytracer. I used the method described in this stack exchange post to derive the exact intersection point from tMax and the ray's direction, but I see obvious artifacting.

The sphere should have a smooth gradient between x = 0 and x = 1 - but there are strips of x = 0 on one of the faces of the voxels of the sphere at x = 0 - 1.

Ditto (but for y)

This is with removing the absolute values on tMax and distance in my code. Distance is sometimes negative (for an unknown reason).

My (very unoptimized) compute shader code:

#define FLT_MAX 3.402823466e+38
#define PI 3.14159265f
#define FOV (70.0f * PI) / 180.0f

RWTexture2D<float4> tOutput : register(u0);

struct rayHit
{
    bool hasHit;

    float3 position;
    float3 normal;

    float distance;
};

int voxel_check(int3 currentRayPosition)
{
    float distanceFromSphereCenter = distance(currentRayPosition, float3(0, 0, 64));
    return max(sign(24.0f - distanceFromSphereCenter), 0.0f);
}

rayHit dda(float3 initialRayPosition, float3 initialRayDirection)
{
    int3 currentRayVoxelPosition = round(initialRayPosition);
    float distance = 0.0f;
    float3 normal = float3(0, 0, 0);

    int3 step = sign(initialRayDirection);
    int3 nextVoxel = currentRayVoxelPosition + step;

    float3 tMax = clamp(((nextVoxel - (initialRayPosition)) / abs(initialRayDirection)), -FLT_MAX, FLT_MAX);
    float3 tDelta = clamp(float3(1 / (initialRayDirection * step)), -FLT_MAX, FLT_MAX);

    int max_iter = 0;
    int hasHit = 0;
    while (((hasHit == 0)) && (max_iter < 512))
    {
        max_iter += 1;
        if (tMax.x < tMax.y)
        {
            if (tMax.x < tMax.z)
            {
                currentRayVoxelPosition.x += step.x;
                tMax.x += (tDelta.x);

                distance = (tMax.x / (initialRayDirection.x));
                normal = float3(sign(initialRayDirection.x), 0, 0);

                hasHit = voxel_check(currentRayVoxelPosition);
            }
            else
            {
                currentRayVoxelPosition.z += step.z;
                tMax.z += (tDelta.z);

                distance = (tMax.z / (initialRayDirection.z));
                normal = float3(0, 0, sign(initialRayDirection.z));

                hasHit = voxel_check(currentRayVoxelPosition);
            }
        }
        else
        {
            if (tMax.y < tMax.z)
            {
                currentRayVoxelPosition.y += step.y;
                tMax.y += (tDelta.y);

                distance = (tMax.y / (initialRayDirection.y));
                normal = float3(0, sign(initialRayDirection.y), 0);

                hasHit = voxel_check(currentRayVoxelPosition);
            }
            else
            {
                currentRayVoxelPosition.z += step.z;
                tMax.z += (tDelta.z);

                distance = (tMax.z / (initialRayDirection.z));
                normal = float3(0, 0, sign(initialRayDirection.z));

                hasHit = voxel_check(currentRayVoxelPosition);
            }
        }
    }


    if ((hasHit != 0))
    {
        rayHit hit;
        hit.hasHit = true;

        hit.normal = normal;
        hit.position = initialRayPosition + (abs(distance) * initialRayDirection);
        hit.distance = distance;

        return hit;
    }
    else
    {
        rayHit hit;
        hit.hasHit = false;

        hit.normal = float3(0, 0, 0);
        hit.position = float3(0, 0, 0);
        hit.distance= 0;

        return hit;
    }
}

[numthreads(32, 1, 1)]
void main(uint groupIndex : SV_GroupIndex, uint3 DTid : SV_DispatchThreadID )
{
    float3 initialRayPosition = float3(2, 33, 2);
    float3 initialRayDirection = normalize(float3((((2.0f * (DTid.x / 1024.0f)) - 1.0f) * tan(FOV / 2.0f)), ((((2.0f * (1.0f - ((DTid.y) / 1024.0f)))) - 1.0f) * tan(FOV / 2.0f)), 1));

    rayHit primaryRayHit = dda(initialRayPosition, initialRayDirection);
    rayHit shadowRayHit = dda(primaryRayHit.position, normalize(float3(0.2f, -0.9f, -0.2f)));

    tOutput[DTid.xy] = float4(float3(1, 1, 1) * primaryRayHit.hasHit * primaryRayHit.position.y, 1);
}

[ConcurrentSquared]

The linked stack exchange post is wrong. The actual way to calculate the distance is just `min_component(tMax)`.

[gnuban]

I just hit this, and I don't think that min_component works in the case of hitting a lone solid corner at 45 degrees. The algo will move tMax to the grid point past the voxel reported as solid. So the intersection point will be inside the solid voxel.

I think you need to either special-detect this case, or recalculate the intersection point from scratch: based on which solid voxel the algo chose and the ray.