Grainy looking sphere in my ray tracer

Question

I am trying to write a simple ray tracer. The final image should like this: I have read stuff about it and below is what I am doing:

 create an empty image (to fill each pixel, via ray tracing)
 for each pixel [for each row, each column]
 create the equation of the ray emanating from our pixel
 trace() ray:
 if ray intersects SPHERE
 compute local shading (including shadow determination)
 return color;

Now, the scene data is like: It sets a gray sphere of radius 1 at (0,0,-3). It sets a white light source at the origin.

 2
 amb: 0.3 0.3 0.3
 sphere
 pos: 0.0 0.0 -3.0
 rad: 1
 dif: 0.3 0.3 0.3
 spe: 0.5 0.5 0.5
 shi: 1
 light
 pos: 0 0 0
 col: 1 1 1

Mine looks very weird :

//check ray intersection with the sphere
boolean intersectsWithSphere(struct point rayPosition, struct point rayDirection,    Sphere sp,float* t){

//float a = (rayDirection.x * rayDirection.x) + (rayDirection.y * rayDirection.y) +(rayDirection.z * rayDirection.z);
// value for a is 1 since rayDirection vector is normalized
double radius = sp.radius;
double xc = sp.position[0];
double yc =sp.position[1];
double zc =sp.position[2];

double xo = rayPosition.x;
double yo = rayPosition.y;
double zo = rayPosition.z;

double xd = rayDirection.x;
double yd = rayDirection.y;
double zd = rayDirection.z;

double b = 2 * ((xd*(xo-xc))+(yd*(yo-yc))+(zd*(zo-zc)));
double c = (xo-xc)*(xo-xc) + (yo-yc)*(yo-yc) + (zo-zc)*(zo-zc) - (radius * radius);
float D = b*b + (-4.0f)*c;

//ray does not intersect the sphere
if(D < 0 ){
    return false;
}

D = sqrt(D);
float t0 = (-b - D)/2 ;
float t1 = (-b + D)/2;

//printf("D=%f",D);
//printf(" t0=%f",t0);
//printf(" t1=%f\n",t1);

if((t0 > 0) && (t1 > 0)){
    *t = min(t0,t1);
    return true;
}
else {
    *t = 0;
    return false;
}

}

Below is the trace() function:

unsigned char* trace(struct point rayPosition, struct point rayDirection, Sphere * totalspheres) {

struct point tempRayPosition = rayPosition;
struct point tempRayDirection = rayDirection;
float  f=0;
float tnear = INFINITY;
boolean sphereIntersectionFound = false;
int sphereIndex = -1;
for(int i=0; i < num_spheres ; i++){
    float t = INFINITY;
    if(intersectsWithSphere(tempRayPosition,tempRayDirection,totalspheres[i],&t)){
        if(t < tnear){
            tnear = t;
            sphereIntersectionFound = true;
            sphereIndex = i;
        }
    }
}

if(sphereIndex < 0){
    //printf("No interesection found\n");
    mycolor[0] = 1;
    mycolor[1] = 1;
    mycolor[2] = 1;
    return mycolor;
}
else {
       Sphere sp = totalspheres[sphereIndex];
        //intersection point
        hitPoint[0].x = tempRayPosition.x + tempRayDirection.x * tnear;
        hitPoint[0].y = tempRayPosition.y + tempRayDirection.y * tnear;
        hitPoint[0].z = tempRayPosition.z + tempRayDirection.z * tnear;

        //normal at the intersection point
        normalAtHitPoint[0].x = (hitPoint[0].x - totalspheres[sphereIndex].position[0])/ totalspheres[sphereIndex].radius;
        normalAtHitPoint[0].y = (hitPoint[0].y - totalspheres[sphereIndex].position[1])/ totalspheres[sphereIndex].radius;
        normalAtHitPoint[0].z = (hitPoint[0].z - totalspheres[sphereIndex].position[2])/ totalspheres[sphereIndex].radius;
        normalizedNormalAtHitPoint[0] = normalize(normalAtHitPoint[0]);

        for(int j=0; j < num_lights ; j++) {

            for(int k=0; k < num_spheres ; k++){

                shadowRay[0].x = lights[j].position[0] - hitPoint[0].x;
                shadowRay[0].y = lights[j].position[1] - hitPoint[0].y;
                shadowRay[0].z = lights[j].position[2] - hitPoint[0].z;
                normalizedShadowRay[0]  = normalize(shadowRay[0]);

                //R = 2 * ( N dot L) * N - L
                reflectionRay[0].x = - 2 * dot(normalizedShadowRay[0],normalizedNormalAtHitPoint[0]) * normalizedNormalAtHitPoint[0].x +normalizedShadowRay[0].x;
                reflectionRay[0].y = - 2 * dot(normalizedShadowRay[0],normalizedNormalAtHitPoint[0]) * normalizedNormalAtHitPoint[0].y +normalizedShadowRay[0].y;
                reflectionRay[0].z = - 2 * dot(normalizedShadowRay[0],normalizedNormalAtHitPoint[0]) * normalizedNormalAtHitPoint[0].z +normalizedShadowRay[0].z;
                normalizeReflectionRay[0] = normalize(reflectionRay[0]);

                        struct point temp;
                        temp.x = hitPoint[0].x + (shadowRay[0].x * 0.0001 );
                        temp.y = hitPoint[0].y + (shadowRay[0].y * 0.0001);
                        temp.z = hitPoint[0].z + (shadowRay[0].z * 0.0001);

                        struct point ntemp = normalize(temp);
                        float f=0;
                        struct point tempHitPoint;
                        tempHitPoint.x = hitPoint[0].x + 0.001;
                        tempHitPoint.y = hitPoint[0].y + 0.001;
                        tempHitPoint.z = hitPoint[0].z + 0.001;

                        if(intersectsWithSphere(hitPoint[0],ntemp,totalspheres[k],&f)){
                        //  if(intersectsWithSphere(tempHitPoint,ntemp,totalspheres[k],&f)){
                            printf("In shadow\n");
                            float r = lights[j].color[0];                       
                            float g = lights[j].color[1];
                            float b = lights[j].color[2];
                            mycolor[0] = ambient_light[0] + r;
                            mycolor[1] = ambient_light[1] + g;
                            mycolor[2] = ambient_light[2] + b;
                            return mycolor;

                    } else {

                        // point is not is shadow , use Phong shading to determine the color of the point.      
                        //I = lightColor * (kd * (L dot N) + ks * (R dot V) ^ sh)       
                        //(for each color channel separately; note that if L dot N < 0, you should clamp L dot N to zero; same for R dot V)

                       float x = dot(normalizedShadowRay[0],normalizedNormalAtHitPoint[0]);
                       if(x < 0)
                           x = 0;

                       V[0].x = - rayDirection.x;
                       V[0].x = - rayDirection.y;
                       V[0].x = - rayDirection.z;
                       normalizedV[0] = normalize(V[0]);
                       float y = dot(normalizeReflectionRay[0],normalizedV[0]);
                       if(y < 0)
                           y = 0;
                       float ar = totalspheres[sphereIndex].color_diffuse[0] * x;
                       float br =  totalspheres[sphereIndex].color_specular[0] * pow(y,totalspheres[sphereIndex].shininess);
                       float r = lights[j].color[0] * (ar+br);
                    //----------------------------------------------------------------------------------
                       float bg =  totalspheres[sphereIndex].color_specular[1] * pow(y,totalspheres[sphereIndex].shininess);
                       float ag = totalspheres[sphereIndex].color_diffuse[1] * x;
                       float g = lights[j].color[1] * (ag+bg);
                    //----------------------------------------------------------------------------------
                       float bb =  totalspheres[sphereIndex].color_specular[2] * pow(y,totalspheres[sphereIndex].shininess);
                       float ab = totalspheres[sphereIndex].color_diffuse[2] * x;
                       float b = lights[j].color[2] * (ab+bb);
                        mycolor[0] =  r + ambient_light[0];
                        mycolor[1] =  g + ambient_light[1];
                        mycolor[2] =  b+ ambient_light[2];
                        return mycolor;
                        } 
                }
        }
    }           
}

The code calling trace() looks like :

void draw_scene()
{
//Aspect Ratio
double a = WIDTH / HEIGHT;
double angel = tan(M_PI * 0.5 * fov/ 180);
ray[0].x = 0.0;
ray[0].y = 0.0;
ray[0].z = 0.0;
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
unsigned int x,y;
float sx, sy;
for(x=0;x < WIDTH;x++)
{
    glPointSize(2.0); 
    glBegin(GL_POINTS);
    for(y=0;y < HEIGHT;y++)
    {
        sx = (((x + 0.5) / WIDTH) * 2.0 ) - 1;
        sy = (((y + 0.5) / HEIGHT) * 2.0 ) - 1;;
        sx  = sx * angel * a;
        sy = sy * angel;
        //set ray direction
        ray[1].x = sx;
        ray[1].y = sy;
        ray[1].z = -1;
        normalizedRayDirection[0] = normalize(ray[1]);
        unsigned char* color = trace(ray[0],normalizedRayDirection[0],spheres);
        unsigned char  x1 = color[0] * 255;
        unsigned char  y1 = color[1] * 255;
        unsigned char  z1 = color[2] * 255;
        plot_pixel(x,y,x1 %256,y1%256,z1%256);
    }
   glEnd();
   glFlush();
 }
}

There could be many, many problems with the code/understanding.

Answer 1

What possible reason do you have for doing this (in the non-shadow branch of trace (...)):

V[0].x = - rayDirection.x;
V[0].x = - rayDirection.y;
V[0].x = - rayDirection.z;

You might as well comment out the first two computations since you write the results of each to the same component. I think you probably meant to do this instead:

V[0].x = - rayDirection.x;
V[0].y = - rayDirection.y;
V[0].z = - rayDirection.z;

That said, you should also avoid using GL_POINT primitives to cover a 2x2 pixel quad. Point primitives are not guaranteed to be square, and OpenGL implementations are not required to support any size other than 1.0. In practice, most support 1.0 - ~64.0 but glDrawPixels (...) is a much better way of writing 2x2 pixels, since it skips primitive assembly and the above mentioned limitations. You are using immediate mode in this example anyway, so glRasterPos (...) and glDrawPixels (...) are still a valid approach.

Answer 2

It seems you are implementing the formula here, but you deviate at the end from the direction the article takes.

First the article warns that D & b can be very close in value, so that -b + D gets you a very limited number. They suggest an alternative.

Also, you are testing that both t0 & t1 > 0. This doesn't have to be true for you to hit the sphere, you could be inside of it (though you obviously should not be in your test scene).

Finally, I would add a test at the beginning to confirm that the direction vector is normalized. I've messed that up more than once in my renderers.

Answer 3

I haven't taken the time to understand all your code, and I'm definitely not a graphics expert, but I believe the problem you have is called "surface acne". In this case it's probably happening because your shadow rays are intersecting with the object itself. What I did in my code to fix this is add epsilon * hitPoint.normal to the shadow ray origin. This effectively moves the ray away from your object a bit, so they don't intersect.

The value I'm using for epsilon is the square root of 1.19209290 * 10^-7, as that is the square root of a constant called EPSILON that is defined in the particular language I'm using.