Intersection problems with ray-sphere intersection

Question

I'm writing a simple ray tracer and to keep it simple for now I've decided to just have spheres in my scene. I am at a stage now where I merely want to confirm that my rays are intersecting a sphere in the scene properly, nothing else. I've created a Ray and Sphere class and then a function in my main file which goes through each pixel to see if there's an intersection (relevant code will be posted below). The problem is that the whole intersection with the sphere is acting rather strangely. If I create a sphere with center (0, 0, -20) and a radius of 1 then I get only one intersection which is always at the very first pixel of what would be my image (upper-left corner). Once I reach a radius of 15 I suddenly get three intersections in the upper-left region. A radius of 18 gives me six intersections and once I reach a radius of 20+ I suddenly get an intersection for EACH pixel so something is acting as it's not supposed to do.

I was suspicious that my ray-sphere intersection code might be at fault here but having looked through it and looked through the net for more information most solutions describe the very same approach I use so I assume it shouldn't(!) be at fault here. So...I am not exactly sure what I am doing wrong, it could be my intersection code or it could be something else causing the problems. I just can't seem to find it. Could it be that I am thinking wrong when giving values for the sphere and rays? Below is relevant code

Sphere class:

Sphere::Sphere(glm::vec3 center, float radius) 
: m_center(center), m_radius(radius), m_radiusSquared(radius*radius)
{
}

//Sphere-ray intersection. Equation: (P-C)^2 - R^2 = 0, P = o+t*d
//(P-C)^2 - R^2 => (o+t*d-C)^2-R^2 => o^2+(td)^2+C^2+2td(o-C)-2oC-R^2
//=> at^2+bt+c, a = d*d, b = 2d(o-C), c = (o-C)^2-R^2
//o = ray origin, d = ray direction, C = sphere center, R = sphere radius
bool Sphere::intersection(Ray& ray) const
{
    //Squared distance between ray origin and sphere center
    float squaredDist = glm::dot(ray.origin()-m_center, ray.origin()-m_center);

    //If the distance is less than the squared radius of the sphere...
    if(squaredDist <= m_radiusSquared)
    {
        //Point is in sphere, consider as no intersection existing
        //std::cout << "Point inside sphere..." << std::endl;
        return false;
    }

    //Will hold solution to quadratic equation
    float t0, t1;

    //Calculating the coefficients of the quadratic equation
    float a = glm::dot(ray.direction(),ray.direction()); // a = d*d
    float b = 2.0f*glm::dot(ray.direction(),ray.origin()-m_center); // b = 2d(o-C)
    float c = glm::dot(ray.origin()-m_center, ray.origin()-m_center) - m_radiusSquared; // c = (o-C)^2-R^2

    //Calculate discriminant
    float disc = (b*b)-(4.0f*a*c);

    if(disc < 0) //If discriminant is negative no intersection happens
    {
        //std::cout << "No intersection with sphere..." << std::endl;
        return false;
    }
    else //If discriminant is positive one or two intersections (two solutions) exists
    {
        float sqrt_disc = glm::sqrt(disc);
        t0 = (-b - sqrt_disc) / (2.0f * a);
        t1 = (-b + sqrt_disc) / (2.0f * a);
    }

    //If the second intersection has a negative value then the intersections
    //happen behind the ray origin which is not considered. Otherwise t0 is
    //the intersection to be considered
    if(t1<0)
    {
        //std::cout << "No intersection with sphere..." << std::endl;
        return false;
    }
    else
    {
        //std::cout << "Intersection with sphere..." << std::endl;
        return true;
    }
}

Program:

#include "Sphere.h"
#include "Ray.h"

void renderScene(const Sphere& s);

const int imageWidth = 400;
const int imageHeight = 400;

int main()
{
    //Create sphere with center in (0, 0, -20) and with radius 10
    Sphere testSphere(glm::vec3(0.0f, 0.0f, -20.0f), 10.0f);

    renderScene(testSphere);

    return 0;
}

//Shoots rays through each pixel and check if there's an intersection with
//a given sphere. If an intersection exists then the counter is increased.
void renderScene(const Sphere& s)
{
    //Ray r(origin, direction)
    Ray r(glm::vec3(0.0f), glm::vec3(0.0f));

    //Will hold the total amount of intersections
    int counter = 0;

    //Loops through each pixel...
    for(int y=0; y<imageHeight; y++)
    {
        for(int x=0; x<imageWidth; x++)
        {
            //Change ray direction for each pixel being processed
            r.setDirection(glm::vec3(((x-imageWidth/2)/(float)imageWidth), ((imageHeight/2-y)/(float)imageHeight), -1.0f));

            //If current ray intersects sphere...
            if(s.intersection(r))
            {
                //Increase counter
                counter++;
            }
        }
    }

    std::cout << counter << std::endl;
}

Gareth Rees Gareth Rees · Accepted Answer · 2012-10-01T17:54:13

Your second solution (t1) to the quadratic equation is wrong in the case disc > 0, where you need something like:

float sqrt_disc = glm::sqrt(disc);
t0 = (-b - sqrt_disc) / (2 * a);
t1 = (-b + sqrt_disc) / (2 * a);

I think it's best to write out the equation in this form rather than turning the division by 2 into a multiplication by 0.5, because the more the code resembles the mathematics, the easier it is to check.

A few other minor comments:

It seemed confusing to re-use the name disc for sqrt(disc), so I used a new variable name above.
You don't need to test for t0 > t1, since you know that both a and sqrt_disc are positive, and so t1 is always greater than t0.
If the ray origin is inside the sphere, it's possible for t0 to be negative and t1 to be positive. You don't seem to handle this case.
You don't need a special case for disc == 0, as the general case computes the same values as the special case. (And the fewer special cases you have, the easier it is to check your code.)

Intersection problems with ray-sphere intersection

2 Answers