/******************************************************************************/
/*                                                                            */
/*  QUATERN - Handle complete program flow for Quaternion display             */
/*                                                                            */
/******************************************************************************/

#define STRICT
#include <windows.h>
#include <commctrl.h>
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <ctype.h>
#include <malloc.h>
#include <new.h>
#include <float.h>

#include "const.h"
#include "classes.h"
#include "extern.h"
#include "funcdefs.h"
#include "julia.rh"

#define DEBUG_TRACQUAT 0
#define DEBUG_DISPRT 0

#define MISS -10.  /* Flags that view ray missed object */

/*
   Flow for quaternion display is as follows

   execute
      For all rows
         Set parameters
         Point to this row in dwell map
         tracquat ( this row ) - Computes distances
            msetdist
      dispRT
         Set up lights
         for each row
            for each column
               given the distances computed by tracquat, compute this displayed pixel
*/


/*
--------------------------------------------------------------------------------

   msetdist_q - For any given point, this computes the approximate distance
                (a fairly tight lower bound) between that point and the
                nearest point to the object's boundary.
                This is the quaternion version.  The regular version,
                used for mandelbrot and Julia sets, is in MSETDIST.CPP

   This algorithm is from Peitgen and Saupe, "The Science of Fractal Images".


--------------------------------------------------------------------------------
*/

static double msetdist_q (
   QUATERNION *start ,   // Z0
   QUATERNION *param ,   // c 
   int max_dwell ,       // Dwell limit, should be as large as practical
   int julia ,           // Julia? (vs. Mandelbrot)
   QUATERNION *orbit     // Work area max_dwell long
   )
{
   int i, iter ;
   double real2, imag2, jmag2, kmag2, real, imag, jmag, kmag ;
   double rd, id, jd, kd, rtemp, itemp, jtemp, ktemp ;

   iter = 0 ;
   real = start->r ;
   imag = start->i ;
   jmag = start->j ;
   kmag = start->k ;

   do {
      orbit[iter].r = real ;
      orbit[iter].i = imag ;
      orbit[iter].j = jmag ;
      orbit[iter].k = kmag ;
      real2 = real * real ;
      imag2 = imag * imag ;
      jmag2 = jmag * jmag ;
      kmag2 = kmag * kmag ;
      if (real2 + imag2 + jmag2 + kmag2 > 100000.0)
         break ;
      rtemp = 2.0 * real ;
      real = real2 - imag2 - jmag2 - kmag2  +  param->r ;
      imag = rtemp * imag + param->i ;
      jmag = rtemp * jmag + param->j ;
      kmag = rtemp * kmag + param->k ;
      } while (++iter < max_dwell) ;

   if (real2 + imag2 + jmag2 + kmag2 > 100000.) {
      rd = 1.0 ;
      id = 0.0 ;
      jd = 0.0 ;
      kd = 0.0 ;
      for (i=0 ; i<iter ; i++) {
         rtemp = 2.0 * (orbit[i].r * rd  -  orbit[i].i * id  -
                        orbit[i].j * jd  -  orbit[i].k * kd) ;
         itemp = 2.0 * (orbit[i].r * id  +  orbit[i].i * rd  +
                        orbit[i].j * kd  -  orbit[i].k * jd) ;
         jtemp = 2.0 * (orbit[i].r * jd  -  orbit[i].i * kd  +
                        orbit[i].j * rd  +  orbit[i].k * id) ;
         ktemp = 2.0 * (orbit[i].r * kd  +  orbit[i].i * jd  -
                        orbit[i].j * id  +  orbit[i].k * rd) ;
         if (! julia)
            rtemp += 1.0 ;
         if ((fabs ( rtemp ) > 1.e30)  ||  (fabs ( itemp ) > 1.e30)  ||
             (fabs ( jtemp ) > 1.e30)  ||  (fabs ( ktemp ) > 1.e30))
            return 0.0 ;
         rd = rtemp ;
         id = itemp ;
         jd = jtemp ;
         kd = ktemp ;
         }
      return 0.25 * log  ( real2 + imag2 + jmag2 + kmag2 ) *     // 0.25 makes this very conservative
                    sqrt ( real2 + imag2 + jmag2 + kmag2 ) /
                    sqrt ( rd * rd  +  id * id  +  jd * jd  +  kd * kd ) ;
      }

   return 0.0 ;
}


/*
--------------------------------------------------------------------------------

   tracquat - Ray tracing part 1 of quaternion display

   For each pixel in the view plane this finds the 4d point (if one exists)
   at which the view ray intersects the object, finding the closest-to-eye
   such point if there are more than one.

--------------------------------------------------------------------------------
*/

void tracquat (
   double Rcent ,        // Object position at display center
   double Icent ,
   double Jcent ,
   double Kcent ,
   double magnif ,       // AMY magnification: 1 / radius.
   int collapse ,        // Dimension to collapse, 0=real thru 3=K
   QUATERNION *jparam ,  // additive constant 'c' for Julia sets
   double elev ,         // Elevation and azimuth (degrees CCW from +X) of view
   double az ,
   double rowpos ,       // .5 to -.5 for top to bottom view (0 in center)
   int ncols ,           // Number of left-to-right columns to do
   double *out ,         // Ncols of t parameter (-10=miss, else -2 (far) to 2)
   QUATERNION *orbit     // Work area maximum iterations long
   )
{
   int col ;
   double Rleft, Ileft, Jleft, Kleft, Rright, Iright, Jright, Kright ;
   double Rcolinc, Icolinc, Jcolinc, Kcolinc, Rrayinc, Irayinc, Jrayinc,Krayinc;
   double Rpt, Ipt, Jpt, Kpt, t, dist, eps ;
   double SinEl, CosEl, SinAz, CosAz, Xtemp, Ytemp, Ztemp ;
   double deg_rad = 3.141592653589793 / 180. ;
   QUATERNION start ;      /* Z0 */

   if (magnif < .5)
      magnif = .5 ;   // Safety factor

   start.r = start.i = start.j = start.k = 0.0 ;  // Not needed; shuts up LINT

/*
   Setup the coordinates of the left and right ends of this line (row)
   as if the elevation and azimuth were 0 (ie viewing from +X).

   Also do the view ray increments, which define the direction of
   moving toward the object from the eye.

   Assume for the moment that collapse=3, so (X,Y,Z) is (R,I,J).
   Fix that later by reassigning axes if necessary.
*/

   Rleft = Rright = 0. ;     // Looking in from +X
   Ileft = - 1. / magnif ;   // I goes across columns of this row
   Iright = 1. / magnif ;
   Jleft = Jright = 2. * rowpos / magnif ;  // J (vertical, Z) corresponds to the row
   Kleft = Kright = 0. ;     // Assume for now K collapsed out per user

   Rrayinc = -1. ;           // Looking in from +X
   Irayinc = Jrayinc = Krayinc = 0. ;

/*
   Now rotate per the view direction.
*/

   SinEl = sin ( deg_rad * elev ) ;
   CosEl = cos ( deg_rad * elev ) ;
   SinAz = sin ( deg_rad * az ) ;
   CosAz = cos ( deg_rad * az ) ;

   // Left end of line (row)
   Xtemp = - Jleft * SinEl ;                    // First we elevate (Assume x=0)
   Ytemp = Ileft ;
   Ztemp = Jleft * CosEl ;
   Rleft = Xtemp * CosAz  -  Ytemp * SinAz ;    // Then we rotate about Z
   Ileft = Xtemp * SinAz  +  Ytemp * CosAz ;
   Jleft = Ztemp ;

   // Right end of line (row)
   Xtemp = - Jright * SinEl ;                   // First we elevate (Assume x=0)
   Ytemp = Iright ;
   Ztemp = Jright * CosEl ;
   Rright = Xtemp * CosAz  -  Ytemp * SinAz ;   // Then we rotate about Z
   Iright = Xtemp * SinAz  +  Ytemp * CosAz ;
   Jright = Ztemp ;

   // View raw direction from eye to object
   Xtemp = Rrayinc * CosEl - Jrayinc * SinEl ;  // First we elevate
   Ytemp = Irayinc ;
   Ztemp = Jrayinc * CosEl + Rrayinc * SinEl ;
   Rrayinc = Xtemp * CosAz  -  Ytemp * SinAz ;  // Then we rotate about Z
   Irayinc = Xtemp * SinAz  +  Ytemp * CosAz ;
   Jrayinc = Ztemp ;

/*
   If collapse is not 3, fix that now by simply swapping axes
*/

   if (collapse == 0) {
      Kleft = Jleft ;
      Jleft = Ileft ;
      Ileft = Rleft ;
      Rleft = 0. ;
      Kright = Jright ;
      Jright = Iright ;
      Iright = Rright ;
      Rright = 0. ;
      Krayinc = Jrayinc ;
      Jrayinc = Irayinc ;
      Irayinc = Rrayinc ;
      Rrayinc = 0. ;
      }
   else if (collapse == 1) {
      Kleft = Jleft ;
      Jleft = Ileft ;
      Ileft = 0. ;
      Kright = Jright ;
      Jright = Iright ;
      Iright = 0. ;
      Krayinc = Jrayinc ;
      Jrayinc = Irayinc ;
      Irayinc = 0. ;
      }
   else if (collapse == 2) {
      Kleft = Jleft ;
      Jleft = 0. ;
      Kright = Jright ;
      Jright = 0. ;
      Krayinc = Jrayinc ;
      Jrayinc = 0. ;
      }

/*
   Position the center of the object to be the user's chosen center of view.
*/

   Rleft += Rcent ;
   Ileft += Icent ;
   Jleft += Jcent ;
   Kleft += Kcent ;
   Rright += Rcent ;
   Iright += Icent ;
   Jright += Jcent ;
   Kright += Kcent ;

/*
   Compute the tiny vector that will be added to the orthogonal view point
   each time we advance one column across the row.
   Moving from the first image column to the last will move us from the left
   end of the row to the right end.
*/

   Rcolinc = (Rright - Rleft) / (double) (ncols-1) ;
   Icolinc = (Iright - Ileft) / (double) (ncols-1) ;
   Jcolinc = (Jright - Jleft) / (double) (ncols-1) ;
   Kcolinc = (Kright - Kleft) / (double) (ncols-1) ;

#if DEBUG_TRACQUAT
   char msg[256] ;
   sprintf_s ( msg, "cent=(%.4lf %.4lf %.4lf %.4lf)", Rcent, Icent, Jcent, Kcent ) ;
   MEMTEXT ( msg ) ;
   sprintf_s ( msg, "Left=(%.4lf %.4lf %.4lf %.4lf)", Rleft, Ileft, Jleft, Kleft ) ;
   MEMTEXT ( msg ) ;
   sprintf_s ( msg, "Right=(%.4lf %.4lf %.4lf %.4lf)", Rright, Iright, Jright,Kright);
   MEMTEXT ( msg ) ;
   sprintf_s ( msg, "Inc=(%.4lf %.4lf %.4lf %.4lf)", Rcolinc,Icolinc,Jcolinc,Kcolinc);
   MEMTEXT ( msg ) ;
   sprintf_s ( msg, "Ray=(%.4lf %.4lf %.4lf %.4lf)", Rrayinc,Irayinc,Jrayinc,Krayinc);
   MEMTEXT ( msg ) ;
#endif

/*
--------------------------------------------------------------------------------

   This is the main loop that does each of the ncols points.

--------------------------------------------------------------------------------
*/

   eps = 0.25 * 2. / (magnif * (double) ncols) ; // Each pixel is 2 / (magnif*nc)
                                                 // Ileft=-1/magnif, Iright=1/magnif
                                                 // We specify accuracy to 1/4 pixel

   for (col=0 ; col<ncols ; col++) { // Move across the image row
      Rpt = Rleft + col * Rcolinc ;  // This is the 4D coordinate
      Ipt = Ileft + col * Icolinc ;  // corresponding to this view pixel
      Jpt = Jleft + col * Jcolinc ;
      Kpt = Kleft + col * Kcolinc ;

#if DEBUG_TRACQUAT
      sprintf_s ( msg, "ViewPt=(%.4lf %.4lf %.4lf %.4lf)", Rpt, Ipt, Jpt, Kpt) ;
      MEMTEXT ( msg ) ;
#endif

/*
   We have a view ray from the eye which passes through the viewpoint,
   and we seek the intersection of this ray with a boundary of the object,
   because the first such intersection is what the eye would see for
   this viewplane pixel.

   We do an orthographic projection, so the eye is at infinity and all rays
   from the eye through the viewpoint (Rpt, Ipt, Jpt, Kpt) are parallel.
   We start the search at the extreme in the direction opposite rayinc
   and move toward the object, thus guaranteeing that if there are multiple
   intersections of the view ray with the object, we will find the one
   closest to the eye, the visible one.  Others would be hidden.

   It's easy to know a decent place to start, because we know from
   mathematics that the entire object is confined within a sphere of
   radius 2 centered at the origin.  So all we have to do is start the
   search on the view ray just outside that sphere and on the 'eye' side
   of the object.  The code below finds such a starting point.

   There is an explicit formula for an intersection point that could be used
   to find the starting point, but the simple-minded approach used here
   is probably not much slower.
   Anyhow, this is not the CPU burning part!  That comes later with msetdist_q().

   We move along the view ray by TestPoint = ViewPoint - t * RayIncrement.
   We need not consider values of the parameter t outside (-4, 4),
   although here I set the limits at 4.1 for a margin of safety.

   Where does the number 4 come from?  It would be senseless for the user
   to specify a view center more than 2 units away from the origin,
   because the object is confined within a sphere of radius 2 around
   the origin.  If the center were at that extreme, it would take
   4 units of t (recall that the length of _rayinc is 1) to get to
   the opposite side of the containing sphere: 2 to reach the origin,
   and 2 more to reach the opposite side of the enclosing sphere.
*/

      t = 4.1 ;          // Start on the eye side of the sphere
      for (;;) {         // Until we reach the enclosing sphere
         if (t < -4.1) { // Passed right by the sphere?
#if DEBUG_TRACQUAT
            MEMTEXT ( "Outside guard" ) ;
#endif
            break ;
            }
         start.r = Rpt - t * Rrayinc ;
         start.i = Ipt - t * Irayinc ;
         start.j = Jpt - t * Jrayinc ;
         start.k = Kpt - t * Krayinc ;
         if (start.r * start.r + start.i * start.i + start.j * start.j +
             start.k * start.k > 4.4)  // Still outside the sphere?
            t -= .5 ;                  // Move a little closer to object
         else {
            t += .5 ;  // Crossed to inside.  Go back to prior t and quit.
            start.r = Rpt - t * Rrayinc ;
            start.i = Ipt - t * Irayinc ;
            start.j = Jpt - t * Jrayinc ;
            start.k = Kpt - t * Krayinc ;
            break ;
            }

         } // Until near but outside sphere

      // For problematic magnification and center we might miss the sphere.
      if (t < -4.1) {
         *out++ = MISS ;
         continue ;
         }

/*
   We now have a point outside but near the containing sphere, with the object
   in the 'rayinc' direction from here and the eye infinitely far behind.

   The tool we will use for closing in on the object's boundary along this ray
   is the msetdist_q() routine seen earlier in this file.  This routine returns
   an approximation (a lower bound due to shrinking) of the distance between a
   given point and the nearest point on the boundary.
   So we just march along the view ray toward the object, keeping track of
   how close we are, and using the return value of msetdist_q() to tell us
   how far to move.

   It's dangerous to make any use of change in the distance (such as decreasing
   and then increasing) because numerical issues may cause changes that
   are not meaningful.  Remember, this is just an approximation.

   But we can count on this return value being a lower bound.
   So we keep marching onward according to the return value.
   If the view ray intersects the boundary, the return value will eventually
   become tiny.

   If, however, the ray does not intersect the boundary (the object is
   usually full of holes), we will just shoot on past.  In fact, once we pass
   the nearest boundary our march will pick up speed fast; the further
   we get from the boundary, the larger will be the return value, and hence
   the larger will be our next jump.  Et cetera.

   So we just pick a large value for the distance of the point from the origin
   and use it as a threshold.  If we pass this threshold, all hope of hitting
   the boundary is lost.  I use 10 as the squared distance here.
*/

#if DEBUG_TRACQUAT
      sprintf_s( msg, "Starting at (%.4lf %.4lf %.4lf %.4lf)",
                 start.r, start.i, start.j, start.k);
      MEMTEXT ( msg ) ;
#endif

      dist = 0.0 ;  // Not needed; shuts up LINT
      for (;;) {
         if (start.r * start.r + start.i * start.i + start.j * start.j +
             start.k * start.k > 10.) {  /* Passed to other side without hit */
            *out++ = MISS ;
#if DEBUG_TRACQUAT
            MEMTEXT ( " Passed through" ) ;
#endif
            break ;
            }
         if (ObjectType == 4)        // Julia quaternion
            dist = msetdist_q ( &start , jparam , DwellLimit , 1 , orbit ) ;
         else if (ObjectType == 3)   // Mandelbrot quaternion
            dist = msetdist_q ( &start , &start , DwellLimit , 0 , orbit ) ;
         else
            assert ( 1 == 2 ) ;
#if DEBUG_TRACQUAT
         sprintf_s ( msg, "(%.4lf %.4lf %.4lf %.4lf)=%.6lf at %.4lf",
                     start.r, start.i, start.j, start.k, dist, t ) ;
         MEMTEXT ( msg ) ;
#endif
         if (dist < eps) {        // This point on the ray is effectively on the object
//            if (dist < 1.e-30)    // Should never happen!  Dangerously close.
//               printf ( "\a" ) ;  // Recall that msetdist_q UNDER-estimates distance
            *out++ = t ;          // We need store only t because we can recreate 4D point from t
#if DEBUG_TRACQUAT
            MEMTEXT ( " HIT" ) ;
#endif
            break ;
            }
         else {
            t -= dist ;           // Move a little closer to the object
            start.r = Rpt - t * Rrayinc ;
            start.i = Ipt - t * Irayinc ;
            start.j = Jpt - t * Jrayinc ;
            start.k = Kpt - t * Krayinc ;
            }
         } /* Until close enough */

      } /* For all columns */
}


/*
--------------------------------------------------------------------------------

   smooth_quat - Smooth the quaternion t matrix

   We allocated twice the required space so we can use the second half
   as scratch storage

--------------------------------------------------------------------------------
*/

void smooth_quat (
   int nrows ,             // Number of rows
   int ncols ,             // And columns
   double *quat_dist       // nrows by ncols matrix of t values from tracquat()
   )
{
   int irow, icol, n ;
   double *tptr, *aptr, *bptr, *dest_ptr, sum, smooth_eps ;

   smooth_eps = 40.0 / (nrows + ncols) ; // See the comment on 'factor' in next routine
                                        // The 4.0 is somewhat arbitrary and may be changed

   for (irow=1 ; irow<nrows-1 ; irow++) {
      for (icol=1 ; icol<ncols-1 ; icol++) {
         tptr = quat_dist + irow * ncols + icol ; // This pixel (source)
         aptr = tptr - ncols ;                    // Pixel above
         bptr = tptr + ncols ;                    // Pixel below
         dest_ptr = tptr + nrows * ncols ;        // Destination
         if (*tptr <= MISS)
            *dest_ptr = MISS ;
         else {
            sum = *tptr ;
            n = 1 ;
            if (*(aptr-1) > MISS  &&  fabs(*(aptr-1) - *tptr) < smooth_eps) {
               sum += *(aptr-1) ;
               ++n ;
               }
            if (*aptr > MISS  &&  fabs(*aptr - *tptr) < smooth_eps) {
               sum += *aptr ;
               ++n ;
               }
            if (*(aptr+1) > MISS  &&  fabs(*(aptr+1) - *tptr) < smooth_eps) {
               sum += *(aptr+1) ;
               ++n ;
               }
            if (*(tptr-1) > MISS  &&  fabs(*(tptr-1) - *tptr) < smooth_eps) {
               sum += *(tptr-1) ;
               ++n ;
               }
            if (*(tptr+1) > MISS  &&  fabs(*(tptr+1) - *tptr) < smooth_eps) {
               sum += *(tptr+1) ;
               ++n ;
               }
            if (*(bptr-1) > MISS  &&  fabs(*(bptr-1) - *tptr) < smooth_eps) {
               sum += *(bptr-1) ;
               ++n ;
               }
            if (*bptr > MISS  &&  fabs(*bptr - *tptr) < smooth_eps) {
               sum += *bptr ;
               ++n ;
               }
            if (*(bptr+1) > MISS  &&  fabs(*(bptr+1) - *tptr) < smooth_eps) {
               sum += *(bptr+1) ;
               ++n ;
               }
            *dest_ptr = sum / n ;
            }
         }
      }

   for (irow=1 ; irow<nrows-1 ; irow++) {
      for (icol=1 ; icol<ncols-1 ; icol++) {
         dest_ptr = quat_dist + irow * ncols + icol ; // Destination
         tptr = dest_ptr + nrows * ncols ;            // This pixel (source)
         aptr = tptr - ncols ;                        // Pixel above
         bptr = tptr + ncols ;                        // Pixel below
         if (*tptr <= MISS)
            *dest_ptr = MISS ;
         else {
            sum = *tptr ;
            n = 1 ;
            if (*(aptr-1) > MISS  &&  fabs(*(aptr-1) - *tptr) < smooth_eps) {
               sum += *(aptr-1) ;
               ++n ;
               }
            if (*aptr > MISS  &&  fabs(*aptr - *tptr) < smooth_eps) {
               sum += *aptr ;
               ++n ;
               }
            if (*(aptr+1) > MISS  &&  fabs(*(aptr+1) - *tptr) < smooth_eps) {
               sum += *(aptr+1) ;
               ++n ;
               }
            if (*(tptr-1) > MISS  &&  fabs(*(tptr-1) - *tptr) < smooth_eps) {
               sum += *(tptr-1) ;
               ++n ;
               }
            if (*(tptr+1) > MISS  &&  fabs(*(tptr+1) - *tptr) < smooth_eps) {
               sum += *(tptr+1) ;
               ++n ;
               }
            if (*(bptr-1) > MISS  &&  fabs(*(bptr-1) - *tptr) < smooth_eps) {
               sum += *(bptr-1) ;
               ++n ;
               }
            if (*bptr > MISS  &&  fabs(*bptr - *tptr) < smooth_eps) {
               sum += *bptr ;
               ++n ;
               }
            if (*(bptr+1) > MISS  &&  fabs(*(bptr+1) - *tptr) < smooth_eps) {
               sum += *(bptr+1) ;
               ++n ;
               }
            *dest_ptr = sum / n ;
            }
         }
      }
}


/*
--------------------------------------------------------------------------------

   dispRT - Display a ray-traced quaternion object.

   This must be given the nrows by ncols matrix of t values that were
   computed by tracquat().

--------------------------------------------------------------------------------
*/

void dispRT (
   int nrows ,             // Number of rows
   int ncols ,             // And columns
   double *quat_dist ,     // nrows by ncols matrix of t values from tracquat()
   unsigned char *image ,  // The nrows by ncols image computed here
   int lighting ,          // Type of lighting
   double magnif           // Magnification
   )
{
   int row, col ;
   double xeye, yeye, zeye, xnor, ynor, znor ;
   double CosEl, SinEl, CosAz, SinAz, CosEl_d, SinEl_d ;
   double x1, y1, z1, x2, y2, z2, x3, y3, z3, d, length, proj ;
   double x12d, y12d, z12d, x23d, y23d, z23d, x31d, y31d, z31d ;
   double x12s, y12s, z12s, x23s, y23s, z23s, x31s, y31s, z31s ;
   double RlightX, RlightY, RlightZ, GlightX, GlightY, GlightZ,
          BlightX, BlightY, BlightZ ;
   double factor ;
   double deg_rad = 3.141592653589793 / 180. ;
   double rhalf = sqrt ( 1. / 2. ) ;
   double rthird = sqrt ( 1. / 3. ) ;
   double *this_row, *row_above, *row_below ;
   unsigned char *iptr, *iptr1, *iptr2 ;

   RlightX = RlightY = RlightZ = GlightX = GlightY = GlightZ = BlightX = BlightY = BlightZ = 0.0 ;  // Not needed; shuts up LINT

   SinEl = sin ( deg_rad * Elevation ) ;
   CosEl = cos ( deg_rad * Elevation ) ;
   SinAz = sin ( deg_rad * Azimuth ) ;
   CosAz = cos ( deg_rad * Azimuth ) ;

/*
   Find the direction of the eye.  This will be needed later to determine
   if a light source is on the same side of the object as the eye.
   It also is needed for specular reflectance of lights, if done some day.
*/

   xeye = CosEl * CosAz ;
   yeye = CosEl * SinAz ;
   zeye = SinEl ;

/*
   Initialize lights
*/

   switch (lighting) {
      case 2:                  // RGB staggered around eye
         /* Red at (1,0,1) */
         RlightX = rhalf * CosAz * (CosEl - SinEl) ;
         RlightY = rhalf * SinAz * (CosEl - SinEl) ;
         RlightZ = rhalf * (CosEl + SinEl) ;
         /* Green at (1,-1,-1) */
         GlightX = rthird * (CosAz * SinEl + CosAz * CosEl + SinAz) ;
         GlightY = rthird * (SinAz * SinEl + SinAz * CosEl - CosAz) ;
         GlightZ = rthird * (SinEl - CosEl) ;
         /* Blue at (1,1,-1) */
         BlightX = rthird * (CosAz * SinEl + CosAz * CosEl - SinAz) ;
         BlightY = rthird * (SinAz * SinEl + SinAz * CosEl + CosAz) ;
         BlightZ = rthird * (SinEl - CosEl) ;
         break ;
      case 3:                // White above
         RlightX = 0. ;
         RlightY = 0. ;
         RlightZ = 1. ;
         break ;
      case 4:                // RGB staggered above
         RlightX = -rhalf ;
         RlightY = 0. ;
         RlightZ = rhalf ;
         GlightX = rthird ;
         GlightY = -rthird ;
         GlightZ = rthird ;
         BlightX = rthird ;
         BlightY = rthird ;
         BlightZ = rthird ;
         break ;
      }

   // This is a scaling factor that controls how changes in t
   // map to changes in the distance in the view ray direction.
   // To understand, assume that magnif=1 and nrows=ncols,
   // which we will call n.  Recall that traversing the image from
   // top to bottom or left to right is a distance of 2 (-1 to 1).
   // Thus, changing one row or column moves a distance of 2/n units.
   // Equivalently, a unit change is a move of n/2 rows or columns.
   // So 'factor' here multiplies t by n/2, keeping the scaling
   // commensurate.  The magnification affects all dimensions equally.
   // If nrows != ncols we lose this perfect scaling, unavoidably,
   // but we do the best we can by averaging the two dimensions.

   factor = magnif * (double) (nrows+ncols) / 4. ; // Scale t commensurate with row/col

#if DEBUG_DISPRT
   char msg[256] ;
   sprintf_s ( msg, "eye at (%.4lf %.4lf %.4lf)  factor=%.4lf", xeye, yeye, zeye, factor ) ;
   MEMTEXT ( msg ) ;
#endif


/*
   Main loop
   Because we need neighbors of each pixel to compute image,
   we cannot do the outer rows and columns
*/

   for (row=1 ; row<nrows-1 ; row++) {
      this_row = quat_dist + row * ncols ;  // Current row in t matrix
      row_above = this_row - ncols ;
      row_below = this_row + ncols ;

      for (col=1 ; col<ncols-1 ; col++) {
         iptr = image + (row * ncols + col) * 3 ;  // 24-bit RGB

         if (this_row[col] <= MISS) {       // This ray misses object?
            iptr[0] = iptr[1] = iptr[2] = 0 ;
            continue ;
            }
         if ((row_above[col] <= MISS) && (row_below[col-1] <= MISS) &&
             (row_below[col+1] <= MISS)) {  // All three neighbors miss?
            iptr[0] = iptr[1] = iptr[2] = 0 ;
            continue ;
            }
         if (! lighting) {    /* No lighting (just distance-based tone)? */
            iptr[0] = iptr[1] = iptr[2] = (unsigned char) (127.0 + 30.0 * this_row[col]) ;
            continue ;
            }

/*
   Get the coordinates of the pixels directly above, below left and below right
*/

         d = row_above[col] * factor ;
         CosEl_d = CosEl * d ;
         SinEl_d = SinEl * d ;
         x1 = CosAz * (CosEl_d - SinEl) ;
         y1 = SinAz * (CosEl_d - SinEl) ;
         z1 = CosEl + SinEl_d ;
         d = row_below[col-1] * factor ;
         CosEl_d = CosEl * d ;
         SinEl_d = SinEl * d ;
         x2 = CosAz * (SinEl + CosEl_d) + SinAz ;
         y2 = SinAz * (SinEl + CosEl_d) - CosAz ;
         z2 = SinEl_d - CosEl ;
         d = row_below[col+1] * factor ;
         CosEl_d = CosEl * d ;
         SinEl_d = SinEl * d ;
         x3 = CosAz * (SinEl + CosEl_d) - SinAz ;
         y3 = SinAz * (SinEl + CosEl_d) + CosAz ;
         z3 = SinEl_d - CosEl ;

/*
   Now compute the surface normal for these 3 points.  Point it toward eye.
*/

         x12d = x1 - x2 ;
         y12d = y1 - y2 ;
         z12d = z1 - z2 ;
         x23d = x2 - x3 ;
         y23d = y2 - y3 ;
         z23d = z2 - z3 ;
         x31d = x3 - x1 ;
         y31d = y3 - y1 ;
         z31d = z3 - z1 ;
         x12s = x1 + x2 ;
         y12s = y1 + y2 ;
         z12s = z1 + z2 ;
         x23s = x2 + x3 ;
         y23s = y2 + y3 ;
         z23s = z2 + z3 ;
         x31s = x3 + x1 ;
         y31s = y3 + y1 ;
         z31s = z3 + z1 ;

         xnor = y12d * z12s  +  y23d * z23s  +  y31d * z31s  ;
         ynor = z12d * x12s  +  z23d * x23s  +  z31d * x31s  ;
         znor = x12d * y12s  +  x23d * y23s  +  x31d * y31s  ;

         length = sqrt( 1.e-20 + xnor * xnor + ynor * ynor + znor * znor ) ;
         proj = xnor * xeye  +  ynor * yeye  +  znor * zeye ;

         if (proj < 0.)
            length = -length ;

         xnor /= length ;
         ynor /= length ;
         znor /= length ;

#if DEBUG_DISPRT > 1
         sprintf_s ( msg, "%d %d: (%.2lf %.2lf %.2lf) (%.2lf %.2lf %.2lf) (%.2lf %.2lf %.2lf)",
                  row, col, x1, y1, z1, x2, y2, z2, x3, y3, z3) ;
         MEMTEXT ( msg ) ;
         sprintf_s ( msg, "Normal=(%.4lf %.4lf %.4lf)", xnor, ynor, znor ) ;
         MEMTEXT ( msg ) ;
#endif


/*
   Compute lighting
*/

         if (lighting == 1)  // Single white light at eye
            iptr[0] = iptr[1] = iptr[2] = (unsigned char) (255.9 * proj / length) ;

         else if (lighting == 3) {  // Single white light above
            proj = xnor * RlightX  +  ynor * RlightY  +  znor * RlightZ ;
            if (proj > 0.)          // Light must be on same side as eye!
               iptr[0] = iptr[1] = iptr[2] = (unsigned char) (255.9 * proj) ;
            else
               iptr[0] = iptr[1] = iptr[2] = (unsigned char) 0 ;
            }

         else if ((lighting == 2)  ||  (lighting == 4)) { // 3 RGB staggered; note Windows RGB order
            proj = xnor * RlightX  +  ynor * RlightY  +  znor * RlightZ ;
            if (proj > 0.0)
               iptr[2] = (unsigned char) (255.9 * proj) ;
            else 
               iptr[2] = 0 ;
            proj = xnor * GlightX  +  ynor * GlightY  +  znor * GlightZ ;
            if (proj > 0.0)
               iptr[1] = (unsigned char) (255.9 * proj) ;
            else 
               iptr[1] = 0 ;
            proj = xnor * BlightX  +  ynor * BlightY  +  znor * BlightZ ;
            if (proj > 0.0)
               iptr[0] = (unsigned char) (255.9 * proj) ;
            else 
               iptr[0] = 0 ;
            } // RGB lighting
         } // For all columns
      } // For all rows

/*
   We just did the inside.  Copy to the outermost rows and columns, display.
*/

   // Copy second row to top row
   iptr1 = image + 3 * ncols ;   // 24-bit RGB; this is the second row
   for (col=1 ; col<ncols-1 ; col++) {
      image[3*col+0] = iptr1[3*col+0] ;
      image[3*col+1] = iptr1[3*col+1] ;
      image[3*col+2] = iptr1[3*col+2] ;
      }

   // Copy second-last row to bottom row
   iptr1 = image + 3 * (nrows-1) * ncols ;  // Bottom row
   iptr2 = iptr1 - 3 * ncols ;
   for (col=1 ; col<ncols-1 ; col++) {
      iptr1[3*col+0] = iptr2[3*col+0] ;
      iptr1[3*col+1] = iptr2[3*col+1] ;
      iptr1[3*col+2] = iptr2[3*col+2] ;
      }

   // Copy left and right columns from those just inside
   for (row=0 ; row<nrows ; row++) {
      iptr = image + 3 * ncols * row ;        // This row
      iptr[0] = iptr[3] ;
      iptr[1] = iptr[4] ;
      iptr[2] = iptr[5] ;
      iptr[3*(ncols-1)+0] = iptr[3*(ncols-2)+0] ;
      iptr[3*(ncols-1)+1] = iptr[3*(ncols-2)+1] ;
      iptr[3*(ncols-1)+2] = iptr[3*(ncols-2)+2] ;
      }
}

