Refraction in Raytracing?

I've been working on my raytracer again. I added reflection and multithreading support. Currently I am working on adding refractions, but its only half working.

As you can see, there is a center sphere(without specular highlight), a reflecting sphere(to the right) and a refracting sphere(left). I'm pretty happy about reflections, it does look very good. For refractions its kinda working...the light is refracted and all shadows of the spheres are visible in the sphere(refraction index 1.4), but there is an outer black ring.

EDIT: Apparently the black ring gets bigger, and therefore the sphere smaller, when I increase the refraction index of the sphere. On the contrary, when decreasing the index of refraction, the Sphere gets larger and the black ring smaller...until, with index of refraction set to one, the ring totally disappears. IOR = 1.9 IOR = 1.1 IOR = 1.00001 And interestingly enough at IOR = 1 the sphere loses its transparency and becomes white.

I think I covered total internal reflection and it is not the issue here.

Now the code: I'm using the operator | for dot product, so (vec|vec) is a dot product and the operator ~ to invert vectors. The objects, both ligths and spheres are stored in Object **objects;. Raytrace function

Colour raytrace(const Ray &r, const int &depth)
{
    //first find the nearest intersection of a ray with an object
    Colour finalColour = skyBlue *(r.getDirection()|Vector(0,0,-1)) * SKY_FACTOR;
    double t, t_min = INFINITY;
    int index_nearObj = -1;
    for(int i = 0; i < objSize; i++)
    {
        if(!dynamic_cast<Light *>(objects[i]))//skip light src
        {
            t = objects[i]->findParam(r);
            if(t > 0 && t < t_min)
            {
                t_min = t;
                index_nearObj = i;
            }
        }
    }
    //no intersection
    if(index_nearObj < 0)
        return finalColour;

    Vector intersect = r.getOrigin() + r.getDirection()*t_min;
    Vector normal = objects[index_nearObj]->NormalAtIntersect(intersect);
    Colour objectColor = objects[index_nearObj]->getColor();
    Ray rRefl, rRefr; //reflected and refracted Ray
    Colour refl = finalColour, refr = finalColour; //reflected and refracted colours
    double reflectance = 0, transmittance = 0;

    if(objects[index_nearObj]->isReflective() && depth < MAX_TRACE_DEPTH)
    {
        //handle reflection
        rRefl = objects[index_nearObj]->calcReflectingRay(r, intersect, normal);
        refl = raytrace(rRefl, depth + 1);
        reflectance = 1;
    }

    if(objects[index_nearObj]->isRefractive() && depth < MAX_TRACE_DEPTH)
    {
        //handle transmission
        rRefr = objects[index_nearObj]->calcRefractingRay(r, intersect, normal, reflectance, transmittance);
        refr = raytrace(rRefr, depth + 1);
    }

    Ray rShadow; //shadow ray
    bool shadowed;
    double t_light = -1;

    Colour localColour;
    Vector tmpv;

    //get material properties
    double ka = 0.2; //ambient coefficient
    double kd; //diffuse coefficient
    double ks; //specular coefficient

    Colour ambient = ka * objectColor; //ambient component
    Colour diffuse, specular;
    double brightness;
    localColour = ambient;
    //look if the object is in shadow or light
    //do this by casting a ray from the obj and
    // check if there is an intersection with another obj
    for(int i = 0; i < objSize; i++)
    {
        if(dynamic_cast<Light *>(objects[i])) //if object is a light
        {
            //for each light
            shadowed = false;
            //create Ray to light
            tmpv = objects[i]->getPosition() - intersect;
            rShadow = Ray(intersect  + (!tmpv) * BIAS, tmpv);
            t_light = objects[i]->findParam(rShadow);

            if(t_light < 0) //no imtersect, which is quite impossible
                continue;

            //then we check if that Ray intersects one object that is not a light
            for(int j = 0; j < objSize; j++)
            {
                    if(!dynamic_cast<Light *>(objects[j]) && j != index_nearObj)//if obj is not a light
                    {
                        t = objects[j]->findParam(rShadow);
                        //if it is smaller we know the light is behind the object
                        //--> shadowed by this light
                        if (t >= 0 && t < t_light)
                        {
                            // Set the flag and stop the cycle
                            shadowed = true;
                            break;
                        }
                    }
            }

            if(!shadowed)
            {
                rRefl = objects[index_nearObj]->calcReflectingRay(rShadow, intersect, normal);
                //reflected ray from ligh src, for ks
                kd = maximum(0.0, (normal|rShadow.getDirection()));
                if(objects[index_nearObj]->getShiny() <= 0)
                    ks = 0;
                else
                    ks = pow(maximum(0.0, (r.getDirection()|rRefl.getDirection())), objects[index_nearObj]->getShiny());
                diffuse = kd * objectColor;// * objects[i]->getColour();
                specular = ks * objects[i]->getColor();
                brightness = 1 /(1 + t_light * DISTANCE_DEPENDENCY_LIGHT);
                localColour += brightness * (diffuse + specular);
            }
        }
    }
    finalColour = localColour + (transmittance * refr + reflectance * refl);
    return finalColour;
}

Now the function that calculates the refracted Ray, I used several different sites for resource, and each had similar algorithms. This is the best I could do so far. It may just be a tiny detail I'm not seeing...

Ray Sphere::calcRefractingRay(const Ray &r, const Vector &intersection,Vector &normal, double & refl, double &trans)const
{
    double n1, n2, n;
    double cosI = (r.getDirection()|normal);
    if(cosI > 0.0)
    {
        n1 = 1.0;
        n2 = getRefrIndex();
        normal = ~normal;//invert
    }
    else
    {
        n1 = getRefrIndex();
        n2 = 1.0;
        cosI = -cosI;
    }
    n = n1/n2;
    double sinT2 = n*n * (1.0 - cosI * cosI);
    double cosT = sqrt(1.0 - sinT2);
    //fresnel equations
    double rn = (n1 * cosI - n2 * cosT)/(n1 * cosI + n2 * cosT);
    double rt = (n2 * cosI - n1 * cosT)/(n2 * cosI + n2 * cosT);
    rn *= rn;
    rt *= rt;
    refl = (rn + rt)*0.5;
    trans = 1.0 - refl;
    if(n == 1.0)
        return r;
    if(cosT*cosT < 0.0)//tot inner refl
    {
        refl = 1;
        trans = 0;
        return calcReflectingRay(r, intersection, normal);
    }
    Vector dir = n * r.getDirection() + (n * cosI - cosT)*normal;
    return Ray(intersection + dir * BIAS, dir);
}

EDIT: I also changed the refraction index around.From

    if(cosI > 0.0)
    {
        n1 = 1.0;
        n2 = getRefrIndex();
        normal = ~normal;
    }
    else
    {
        n1 = getRefrIndex();
        n2 = 1.0;
        cosI = -cosI;
    }

if(cosI > 0.0)
{
    n1 = getRefrIndex();
    n2 = 1.0;
    normal = ~normal;
}
else
{
    n1 = 1.0;
    n2 = getRefrIndex();
    cosI = -cosI;
}

Then I get this, and almost the same(still upside down) with an index of refraction at 1! And the reflection calculation:

Ray Sphere::calcReflectingRay(const Ray &r, const Vector &intersection, const Vector &normal)const
{
    Vector rdir = r.getDirection();
    Vector dir = rdir - 2 * (rdir|normal) * normal;
    return Ray(intersection + dir*BIAS, dir);
    //the Ray constructor automatically normalizes directions
}

So my question is: How do I fix the outer black circle? Which version is correct?

Help is greatly appreciated :)

This is compiled on Linux using g++ 4.8.2.

回答1:

Warning: the following is a guess, not a certainty. I'd have to look at the code in more detail to be sure what's happening and why.

That said, it looks to me like your original code is basically simulating a concave lens instead of convex.

A convex lens is basically a magnifying lens, bringing light rays from a relatively small area into focus on a plane:

This also shows why the corrected code shows an upside-down image. The rays of light coming from the top on one side get projected to the bottom on the other (and vice versa).

Getting back to the concave lens though: a concave lens is a reducing lens that shows a wide angle of picture from in front of the lens:

If you look at the bottom right corner here, it shows what I suspect is the problem: especially with a high index of refraction, the rays of light trying to come into the lens intersect the edge of the lens itself. For all the angles wider than that, you're typically going to see a black ring, because the front edge of the lens is acting as a shade to prevent light from entering.

Increasing the index of refraction increases the width of that black ring, because the light is bent more, so a larger portion at the edges is intersecting the outer edge of the lens.

In case you care about how they avoid this with things like wide-angle camera lenses, the usual route is to use a meniscus lens, at least for the front element:

This isn't a panacea, but does at least prevent incoming light rays from intersecting the outer edge of the front lens element. Depending on exactly how wide an angle the lens needs to cover, it'll often be quite a bit less radical of a meniscus than this (and in some cases it'll be a plano-concave) but you get the general idea.

Final warning: of course, all of these are hand-drawn, and intended only to give general idea, not (for example) reflect the design of any particular lens, an element with any particular index of refraction, etc.

回答2:

I stumbled across this exact issue as well when working on a ray tracer. @lightxbulb's comment about normalizing the ray direction vector fixed this problem for me.

Firstly, keep your code that computes the refraction indices prior to your edit. In other words, you should be seeing those black rings in your renderings.

Then, in your calcRefractingRay function where you compute cosI, use the dot product of normalize(r.getDirection()) and normal. Currently you're taking the dot product of r.getDirection() and normal.

Secondly, when you compute the refracted ray direction dir, use normalize(r.getDirection()) instead of r.getDirection(). Again, you're currently using r.getDirection() in your calculation.

Also, there is an issue with the way you're checking for total internal reflection. You should check that the term you're taking the square root of (1.0 - sinT2) is non-negative before actually computing the square root.

Hope that helps!