Зернистая сфера в моем трассировщике лучей

Я пытаюсь написать простой трассировщик лучей. Окончательное изображение должно выглядеть следующим образом: Я читал об этом, и ниже показано, что я делаю. :

 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;

Теперь данные сцены такие: серая сфера радиуса 1 устанавливается в (0,0,-3). Он устанавливает источник белого света в начале координат.

 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

Мой выглядит очень странно:

//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;
}

}

Ниже приведена функция trace():

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;
                        } 
                }
        }
    }           
}

Код, вызывающий trace(), выглядит так:

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();
 }
}

Может быть много-много проблем с кодом/пониманием.