GLSL Sphere - ray intersection geometric solution

Question

I'm trying to implement sphere ray intersection in GLSL, both the geometric and analytical solution. I'm having trouble solving the geom one, it should have something to do with how I return true or false:

bool hitSphere(Ray ray, Sphere sphere, float t_min, float t_max, out float t_out) {
    // Geometric solution
    float R2 = sphere.radius * sphere.radius;
    vec3 L = sphere.position - ray.origin;
    float tca = dot(L, normalize(ray.direction));
    // if(tca < 0) return false;

    float D2 = dot(L, L) - tca * tca;
    if(D2 > R2) return false;
    float thc = sqrt(R2 - D2);
    float t0 = tca - thc;
    float t1 = tca + thc;

    if(t0 < t_max && t0 > t_min) {
        t_out = t0;
        return true;
    }

    if(t1 < t_max && t1 > t_min) {
        t_out = t1;
        return true;
    }
    
    return false; 
}

I think the problem is with how I deal with t0 and t1 for none, one or both intersection cases.

Edit: the analytic version that does work:

vec3 oc = ray.origin - sphere.position;
float a = dot(ray.direction, ray.direction);
float b = dot(oc, ray.direction);
float c = dot(oc, oc) - sphere.radius * sphere.radius;

float discriminant = b * b - a * c;

if (discriminant > 0.0f) {
    if(b > 0) 
        t_out = (-b + sqrt(discriminant)) / a;
    else 
        t_out = (-b - sqrt(discriminant)) / a;
    
    if(t_out < t_max && t_out > t_min) {
        return true;
    }
}

return false;

Answer 1

The issue is caused by t_out. The algorithm has to compute t_out in that way, that X is the intersected point of the ray and the surface of the sphere, for:

X = ray.origin + ray.direction * t_out; 

In the working algorithm t_out depends on the length of ray.direction. t_out becomes smaller, if the magnitude of the vector ray.direction is greater.
In the algorithm, which doesn't work, ray.direction is normalized.

float tca = dot(L, normalize(ray.direction));

Hence t_out is computed for a ray direction length of 1. Actually you compute a t_out' where t_out' = t_out * length(ray.direction).

Divide t0 respectively t1 by the length of ray.direction:

bool hitSphere_2(Ray ray, Sphere sphere, float t_min, float t_max, out float t_out)
{
    float R2 = sphere.radius * sphere.radius;
    vec3 L = sphere.position - ray.origin;
    float tca = dot(L, normalize(ray.direction));
    // if(tca < 0) return false;

    float D2 = dot(L, L) - tca * tca;
    if(D2 > R2) return false;
    float thc = sqrt(R2 - D2);
    float t0 = tca - thc;
    float t1 = tca + thc;

    if (t0 < t_max && t0 > t_min) {
        t_out = t0 / length(ray.direction); // <---
        return true;
    } 

    if (t1 < t_max && t1 > t_min) {
        t_out = t1 / length(ray.direction); // <---
        return true;
    }

    return false; 
}