Listing 2: The lambert function

#include "cips.h"

#define SOURCE         100
#define AMBIENT        200
#define DEGREESRADIANS 57.29577951
#define NINETYDEGREES  1.570796327
#define ONEEIGHTY      3.141592654

lambert(the_image, out_image,
        k_diffuse, k_specular, eta,
        L,
        rows, cols)

   float  eta, k_diffuse, k_specular, *L;
   short  **the_image,
          **out_image;
   long   cols, rows;
{
   float  theta1, theta2, theta3, theta4;
   float  diffuse_term, specular_term;
   float  nl, vr;
   float  N[3], R[3], V[3], v1[3], v2[3];
   float  l, n, r, v;
   int    i, j;


      /******************************************
      *
      *   V, the viewer, is straight up
      *
      ******************************************/

   V[0] =  0.0;
   V[1] =  0.0;
   V[2] = -1.0;


   for(i=1; i<rows-1; i++){
      for(j=1; j<cols-1; j++){

            /*************************************
            *
            *   v1 is the vector down one row
            *   v2 is the vector right one column
            *   v1 and v2 form the surface plane
            *   Use them to calculate the vector N
            *
            *************************************/

         v1[0] = 0.0;
         v1[1] = 1.0;
         v1[2] = (float)(the_image[i][j] 
                 - the_image[i+1][j]);
         v2[0] = 1.0;
         v2[1] = 0.0;
         v2[2] = (float)(the_image[i][j] 
                 - the_image[i][j+1]);

         cross_product(v1, v2, N);

            /*************************************
            *
            *   theta1 = angle between L and N
            *   theta2 = angle between N and V
            *   theta3 = angle between L and V
            *   theta4 = angle between V and R
            *
            *************************************/

         angle_between(L, N, &theta1);
         theta1 = ONEEIGHTY - theta1;
         angle_between(N, V, &theta2);
         angle_between(L, V, &theta3);
         theta3 = ONEEIGHTY - theta3;
         if(theta3 >= theta1)
            theta4 = theta1 - theta2;
         else
            theta4 = theta1 + theta2;

         vr = cos(theta4);
         if(vr <= 0.0)
            vr = 0.0;
         vr = pow(vr, eta);
         specular_term = k_specular*vr;

         magnitude_of(N, &n);
         magnitude_of(L, &l);
         nl           = cos(theta1);
         diffuse_term = k_diffuse*nl;

            /*************************************
            *
            *   If the surface is in the shadow
            *   of the light, the diffuse and
            *   specular terms are zero. 
            *
            *************************************/

         if(theta1 >= NINETYDEGREES){
            diffuse_term  = 0.0;
            specular_term = 0.0;
         }  /* ends if */


         out_image[i][j] = k_diffuse*AMBIENT + 
              SOURCE*(specular_term+diffuse_term);
         if(out_image[i][j] < 0) 
            out_image[i][j] = 0;
         if(out_image[i][j] > GRAY_LEVELS-1)
            out_image[i][j] = GRAY_LEVELS-1;

      }  /* ends loop over j */
   }  /* ends loop over i */

   .
   .
   .

}  /* ends lambert */