/******************************************************************************/
/*                                                                            */
/*  COMPUTE_3DAG - All routines for computing 3D axial generation             */
/*                                                                            */
/******************************************************************************/

#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 "const.h"
#include "classes.h"
#include "funcdefs.h"
#include "lin.rh"
#include "extern.h"

#define DEBUG 0

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

   Generate a uniform in [0, 1).

   This version is low but reasonable quality and thread-safe
   It must NEVER be called with iparam=0 and will never return 0

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

static int irand = 12345 ;

#define IA 16807
#define IM 2147483647
#define IQ 127773
#define IR 2836

double fast_unif ( int *iparam )
{
   long k ;

   k = *iparam / IQ ;
   *iparam = IA * (*iparam - k * IQ) - IR * k ;
   if (*iparam < 0)
      *iparam += IM ;
   return *iparam / (double) IM ;
}
/*
--------------------------------------------------------------------------------

   Local routine counts how many destination items will be needed

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

static int count_advance (
   AG3D_Params *params , 
   int n ,         // Number of tokens in source
   int *source     // Source string (predecessor)
   )
{
   int i, j, k, maxlen, id ;

   j = 0 ;

   for (i=0 ; i<n ; i++) {
      if (source[i] == P_PLUS)           // Just copy turns, brackets, and shrinks
         ++j ;
      else if (source[i] == P_MINUS)
         ++j ;
      else if (source[i] == P_PUSH)
         ++j ;
      else if (source[i] == P_POP)
         ++j ;
      else if (source[i] == P_START_BRACE)
         ++j ;
      else if (source[i] == P_END_BRACE)
         ++j ;
      else if (source[i] == P_SHRINK)
         ++j ;
      else if (source[i] == P_STEP)
         ++j ;
      else if (source[i] == P_PITCH_DOWN)
         ++j ;
      else if (source[i] == P_PITCH_UP)
         ++j ;
      else if (source[i] == P_ROLL_LEFT)
         ++j ;
      else if (source[i] == P_ROLL_RIGHT)
         ++j ;
      else if (source[i] == P_TURN_AROUND)
         ++j ;
      else if (source[i] == P_COLOR)
         j += 2 ;
      else if (source[i] == P_FORWARD) {
         if (params->erp == NULL)
            ++j ;
         else
            j += params->n_erp ;
         }
      else if (source[i] == P_NAME) {
         id = source[++i] ;
         maxlen = 0 ;
         for (k=0 ; k<params->prod_nrand[id] ; k++) {
            if (params->prod_nops[id][k] > maxlen)
               maxlen = params->prod_nops[id][k] ;
            }
         j += maxlen ;
         }
      }

   return j ;
}

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

   Local routine does one iteration of string expansion

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

static int advance (
   AG3D_Params *params , 
   int n ,         // Number of tokens in source
   int *source ,   // Source string (predecessor)
   int *dest       // Destination string (successor)
   )
{
   int i, j, k, m, id ;
   double ff ;

   j = 0 ;

   for (i=0 ; i<n ; i++) {
      if (source[i] == P_PLUS)           // Just copy turns, brackets, and shrinks
         dest[j++] = P_PLUS ;
      else if (source[i] == P_MINUS)
         dest[j++] = P_MINUS ;
      else if (source[i] == P_PUSH)
         dest[j++] = P_PUSH ;
      else if (source[i] == P_POP)
         dest[j++] = P_POP ;
      else if (source[i] == P_START_BRACE)
         dest[j++] = P_START_BRACE ;
      else if (source[i] == P_END_BRACE)
         dest[j++] = P_END_BRACE ;
      else if (source[i] == P_SHRINK)
         dest[j++] = P_SHRINK ;
      else if (source[i] == P_STEP)
         dest[j++] = P_STEP ;
      else if (source[i] == P_PITCH_DOWN)
         dest[j++] = P_PITCH_DOWN ;
      else if (source[i] == P_PITCH_UP)
         dest[j++] = P_PITCH_UP ;
      else if (source[i] == P_ROLL_LEFT)
         dest[j++] = P_ROLL_LEFT ;
      else if (source[i] == P_ROLL_RIGHT)
         dest[j++] = P_ROLL_RIGHT ;
      else if (source[i] == P_TURN_AROUND)
         dest[j++] = P_TURN_AROUND ;
      else if (source[i] == P_COLOR) {
         dest[j++] = P_COLOR ;
         dest[j++] = source[++i] ;  // The selected color, 0-9
         }
      else if (source[i] == P_FORWARD) {
         if (params->erp == NULL)
            dest[j++] = P_FORWARD ;
         else {
            for (k=0 ; k<params->n_erp ; k++)
               dest[j++] = params->erp[k] ;
            }
         }
      else if (source[i] == P_NAME) {
         id = source[++i] ;
         ff = 0.999999 * fast_unif ( &irand ) ;
         for (m=0 ; m<params->prod_nrand[id] ; m++) {
            if (ff <= params->prod_prob[id][m])
               break ;
            }
         assert ( m < params->prod_nrand[id] ) ;
         for (k=0 ; k<params->prod_nops[id][m] ; k++)
            dest[j++] = params->prod_ops[id][m][k] ;
         }
      }

   return j ;
}


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

   Local routine counts how many nodes will be produced by turtle()

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


static int count_turtle (
   int n ,          // Number of operations in 'opstr'
   int *opstr       // The operation string
   )
{
   int i, node ;

   node = 1 ;           // Index into x and y arrays

   for (i=0 ; i<n ; i++) {
      if (opstr[i] == P_FORWARD)      // Take a normal step forward
         ++node ;
      if (opstr[i] == P_STEP)         // Take a small step for outline of filled plane
         ++node ;
      else if (opstr[i] == P_POP)     // Pop state from stack
         ++node ;
      else if (opstr[i] == P_END_BRACE)  // Pop state from stack
         ++node ;
      }

   return node ;
}


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

   Local routine converts an operation string to turtle-move node coordinates.
   This will leave x, y, and z in the range 0 - 1,
   with all three reaching a minimum of 0, and at least one of them
   reaching a maximum of 1.

   It assumes a starting position of (0,0,0), a heading of (0,0,1), pointing
   straight up the positive Z axis, and the turtle's 'up' (0,-1,0), pointing
   toward the observer out the -Y axis.
   The default color is color 0, which will be black if n_colors=0.

   It returns the number of (X,Y,Z) nodes.

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

static int turtle (
   int n ,             // Number of operations in 'opstr'
   int *opstr ,        // The operation string
   AG3D_Params *params ,
   double *x ,         // Output of x coordinates of nodes
   double *y ,         // And y coordinates
   double *z ,         // And z coordinates
   double *thickness , // Thickness of drawn line
   int *draw_line ,    // Color of line from prior to this point
                       //   If normal (F) line, color index+1
                       //   If step (S) outline, -(color index+1)
                       //   If ] or } no line; just move
   double *stack_x ,   // Pushed location work area
   double *stack_y ,
   double *stack_z ,
   double *stack_hx ,  // Pushed Heading work area
   double *stack_hy ,
   double *stack_hz ,
   double *stack_ux ,  // Pushed 'up' work area
   double *stack_uy ,
   double *stack_uz ,
   double *stack_t ,   // Pushed line thickness work area
   int *stack_c        // Current color work area
   )
{
   int i, node, kstack, current_color ;
   double current_x, current_y, current_z, angle ;
   double current_hx, current_hy, current_hz, current_ux, current_uy, current_uz ;
   double xmin, xmax, ymin, ymax, zmin, zmax, scale, tx, ty, tz ;
   double cos_a, sin_a, left_x, left_y, left_z, current_thickness ;

   angle = params->angle * PI / 180.0 ;  // Convert degrees to radians
   cos_a = cos ( angle ) ;
   sin_a = sin ( angle ) ;
   current_thickness = params->width ;
   current_color = 0 ;
   current_x = current_y = current_z = 0.0 ;
   current_hx = current_hy = 0.0 ;
   current_hz = 1.0 ;
   current_ux = current_uz = 0.0 ;
   current_uy = -1.0 ;

   xmin = ymin = zmin = 1.e30 ;  // Boundary of the figure
   xmax = ymax = zmax = -1.e30 ;

   x[0] = y[0] = z[0] = 0.0 ;  // Path starts here
   node = 1 ;                  // Index into node arrays
   kstack = 0 ;

   for (i=0 ; i<n ; i++) {
      left_x = current_uy * current_hz - current_uz * current_hy ; // We only need these for turn/roll/pitch
      left_y = current_uz * current_hx - current_ux * current_hz ; // But they are fast to compute and this
      left_z = current_ux * current_hy - current_uy * current_hx ; // avoids lots of repetition

      if (opstr[i] == P_FORWARD) {        // Take a normal step forward
#if DEBUG
         audit ( "F" ) ;
#endif
         current_x += current_hx ;
         current_y += current_hy ;
         current_z += current_hz ;
         if (current_x < xmin)
            xmin = current_x ;
         if (current_x > xmax)
            xmax = current_x ;
         if (current_y < ymin)
            ymin = current_y ;
         if (current_y > ymax)
            ymax = current_y ;
         if (current_z < zmin)
            zmin = current_z ;
         if (current_z > zmax)
            zmax = current_z ;
         x[node] = current_x ;
         y[node] = current_y ;
         z[node] = current_z ;
         thickness[node] = current_thickness ;
         draw_line[node++] = current_color + 1 ;
         }

      else if (opstr[i] == P_STEP) { // Take a small step forward for filled outline
#if DEBUG
         audit ( "S" ) ;
#endif
         current_x += params->stepsize * current_hx ;
         current_y += params->stepsize * current_hy ;
         current_z += params->stepsize * current_hz ;
         if (current_x < xmin)
            xmin = current_x ;
         if (current_x > xmax)
            xmax = current_x ;
         if (current_y < ymin)
            ymin = current_y ;
         if (current_y > ymax)
            ymax = current_y ;
         if (current_z < zmin)
            zmin = current_z ;
         if (current_z > zmax)
            zmax = current_z ;
         x[node] = current_x ;
         y[node] = current_y ;
         z[node] = current_z ;
         thickness[node] = current_thickness ; // Not needed
         draw_line[node++] = -(current_color + 1) ;
         }
      else if (opstr[i] == P_PLUS) {        // Turn CCW
#if DEBUG
         audit ( "+" ) ;
#endif
         current_hx = current_hx * cos_a + left_x * sin_a ;
         current_hy = current_hy * cos_a + left_y * sin_a ;
         current_hz = current_hz * cos_a + left_z * sin_a ;
         }
      else if (opstr[i] == P_MINUS) {       // Turn CW
#if DEBUG
         audit ( "-" ) ;
#endif
         current_hx = current_hx * cos_a - left_x * sin_a ;
         current_hy = current_hy * cos_a - left_y * sin_a ;
         current_hz = current_hz * cos_a - left_z * sin_a ;
         }
      else if (opstr[i] == P_ROLL_LEFT) {   // Roll left
#if DEBUG
         audit ( "\\" ) ;
#endif
         current_ux = current_ux * cos_a + left_x * sin_a ;
         current_uy = current_uy * cos_a + left_y * sin_a ;
         current_uz = current_uz * cos_a + left_z * sin_a ;
         }
      else if (opstr[i] == P_ROLL_RIGHT) {  // Roll right
#if DEBUG
         audit ( "/" ) ;
#endif
         current_ux = current_ux * cos_a - left_x * sin_a ;
         current_uy = current_uy * cos_a - left_y * sin_a ;
         current_uz = current_uz * cos_a - left_z * sin_a ;
         }
      else if (opstr[i] == P_PITCH_UP) {    // Pitch up
#if DEBUG
         audit ( "^" ) ;
#endif
         tx = current_hx ;
         ty = current_hy ;
         tz = current_hz ;
         current_hx = current_hx * cos_a + current_ux * sin_a ;
         current_hy = current_hy * cos_a + current_uy * sin_a ;
         current_hz = current_hz * cos_a + current_uz * sin_a ;
         current_ux = current_ux * cos_a - tx * sin_a ;
         current_uy = current_uy * cos_a - ty * sin_a ;
         current_uz = current_uz * cos_a - tz * sin_a ;
         }
      else if (opstr[i] == P_PITCH_DOWN) {    // Pitch down
#if DEBUG
         audit ( "&" ) ;
#endif
         tx = current_hx ;
         ty = current_hy ;
         tz = current_hz ;
         current_hx = current_hx * cos_a - current_ux * sin_a ;
         current_hy = current_hy * cos_a - current_uy * sin_a ;
         current_hz = current_hz * cos_a - current_uz * sin_a ;
         current_ux = current_ux * cos_a + tx * sin_a ;
         current_uy = current_uy * cos_a + ty * sin_a ;
         current_uz = current_uz * cos_a + tz * sin_a ;
         }
      else if (opstr[i] == P_TURN_AROUND) {        // Turn around (heading 180)
#if DEBUG
         audit ( "|" ) ;
#endif
         current_hx = -current_hx ;
         current_hy = -current_hy ;
         current_hz = -current_hz ;
         }
      else if (opstr[i] == P_SHRINK) {  // Shrink line thickness
#if DEBUG
         audit ( "V" ) ;
#endif
         current_thickness *= params->shrink ;
         if (current_thickness < 1.0)
            current_thickness = 1.0 ;
         }
      else if (opstr[i] == P_COLOR) {         // Set line color
         current_color = opstr[++i] ;
         assert ( current_color >= 0  &&  current_color <= 9 ) ;
         }
      else if (opstr[i] == P_PUSH  ||  opstr[i] == P_START_BRACE) {  // Push state onto stack
#if DEBUG
         audit ( "Push" ) ;
#endif
         stack_x[kstack] = current_x ;
         stack_y[kstack] = current_y ;
         stack_z[kstack] = current_z ;
         stack_hx[kstack] = current_hx ;
         stack_hy[kstack] = current_hy ;
         stack_hz[kstack] = current_hz ;
         stack_ux[kstack] = current_ux ;
         stack_uy[kstack] = current_uy ;
         stack_uz[kstack] = current_uz ;
         stack_t[kstack] = current_thickness ;
         stack_c[kstack] = current_color ;
         ++kstack ;
         }
      else if (opstr[i] == P_POP  ||  opstr[i] == P_END_BRACE) {   // Pop state from stack
#if DEBUG
         audit ( "Pop" ) ;
#endif
         --kstack ;
         current_x = stack_x[kstack] ;
         current_y = stack_y[kstack] ;
         current_z = stack_z[kstack] ;
         current_hx = stack_hx[kstack] ;
         current_hy = stack_hy[kstack] ;
         current_hz = stack_hz[kstack] ;
         current_ux = stack_ux[kstack] ;
         current_uy = stack_uy[kstack] ;
         current_uz = stack_uz[kstack] ;
         current_thickness = stack_t[kstack] ;
         current_color = stack_c[kstack] ;
         x[node] = current_x ;
         y[node] = current_y ;
         z[node] = current_z ;
         // No need to put thickness in node because no line drawn; just moved
         thickness[node] = 1.e30 ;  // Deliberately create chaos if we have a bug
         draw_line[node++] = 0 ;    // Flag that no line is to be drawn; just move
         }
      else if (opstr[i] == P_NAME) {
#if DEBUG
         audit ( "Name" ) ;
#endif
         ++i ;   // Must bypass index of production
         }
      }

/*
   The nodes are computed.
   Shift and scale them so that the figure touches the bottom and left
   sides, as well as either the top or right side or both.
*/

   if (xmax-xmin >= ymax-ymin  &&  xmax-xmin >= zmax-zmin)
      scale = 1.0 / (xmax - xmin + 1.e-30) ;
   else if (ymax-ymin >= xmax-xmin  &&  ymax-ymin >= zmax-zmin)
      scale = 1.0 / (ymax - ymin + 1.e-30) ;
   else if (zmax-zmin >= xmax-xmin  &&  zmax-zmin >= ymax-ymin)
      scale = 1.0 / (zmax - zmin + 1.e-30) ;
   else {
      scale = 1.0;
      assert(1 == 2) ;
      }

   for (i=0 ; i<node ; i++) {
      x[i] = scale * (x[i] - xmin) ;
      y[i] = scale * (y[i] - ymin) ;
      z[i] = scale * (z[i] - zmin) ;
      }

#if DEBUG
   char msg[256] ;
   audit ( "" ) ;
   audit ( "Pre-rotation nodes..." ) ;
   for (i=0 ; i<node ; i++) {
      sprintf_s ( msg, "%2d  %7.4lf %7.4lf %7.4lf", draw_line[i], x[i], y[i], z[i] ) ;
      audit ( msg ) ;
      }
#endif

   return node ;
}


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

   Local routine called from build_image() draws a filled circle

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

static void circle (
   double diameter ,   // Diameter of the circle (sort of; see comments)
   int color ,         // Index (0-9) of color, ignored and black if n_colors=0
   int row ,           // Row location of center (relative to top)
   int col ,           // And column
   int nr ,            // Row dimension of img
   int nc ,            // And column
   unsigned char *img  // The image
   )
{
   int i, j, extra, horz ;
   double radius ;

   radius = ((int) diameter) / 2 + 0.001 ;  // Not really radius, but captures the idea
   extra = (int) radius ;           // Number of extra pixels

   // Draw a single horizontal line with (row,col) at center
   for (i=col-extra ; i<=col+extra ; i++) {
      if (i >= 0  &&  i < nc) {
         if (AG3D_params.n_colors == 0) {
            img[3*(row*nc+i)+0] = 0 ;
            img[3*(row*nc+i)+1] = 0 ;
            img[3*(row*nc+i)+2] = 0 ;
            }
         else {
            img[3*(row*nc+i)+0] = (unsigned char)AG3D_params.colors[color][0] ;
            img[3*(row*nc+i)+1] = (unsigned char)AG3D_params.colors[color][1] ;
            img[3*(row*nc+i)+2] = (unsigned char)AG3D_params.colors[color][2] ;
            }
         }
      }

   // Draw horizontal lines above and below
   for (j=1 ; j<=extra ; j++) {
      horz = (int) (sqrt ( radius * radius - j * j )) ;
      for (i=col-horz ; i<=col+horz ; i++) {
         if (i >= 0  &&  i < nc  &&  row-j >= 0) {
            if (AG3D_params.n_colors == 0) {
               img[3*((row-j)*nc+i)+0] = 0 ;
               img[3*((row-j)*nc+i)+1] = 0 ;
               img[3*((row-j)*nc+i)+2] = 0 ;
               }
            else {
               img[3*((row-j)*nc+i)+0] = (unsigned char)AG3D_params.colors[color][0] ;
               img[3*((row-j)*nc+i)+1] = (unsigned char)AG3D_params.colors[color][1] ;
               img[3*((row-j)*nc+i)+2] = (unsigned char)AG3D_params.colors[color][2] ;
               }
            }
         if (i >= 0  &&  i < nc  &&  row+j < nr) {
            if (AG3D_params.n_colors == 0) {
               img[3*((row+j)*nc+i)+0] = 0 ;
               img[3*((row+j)*nc+i)+1] = 0 ;
               img[3*((row+j)*nc+i)+2] = 0 ;
               }
            else {
               img[3*((row+j)*nc+i)+0] = (unsigned char)AG3D_params.colors[color][0] ;
               img[3*((row+j)*nc+i)+1] = (unsigned char)AG3D_params.colors[color][1] ;
               img[3*((row+j)*nc+i)+2] = (unsigned char)AG3D_params.colors[color][2] ;
               }
            }
         }
      }
}


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

   Local routine builds an image from an array of nodes

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

static void build_image (
   int n ,              // Number of nodes
   double *x ,          // They are here
   double *y ,          // Must range from 0 to slightly less than 1
   double *z ,          // Must range from 0 to slightly less than 1
   double *thickness ,  // Line thickness from prior to this point
   int *draw_line ,     // Color of line from prior to this point
                        //   If normal (F) line, color index+1
                        //   If step (S) outline, -(color index+1)
                        //   If ] or } no line; just move
   int nr ,             // Image dimensions
   int nc ,
   unsigned char *img , // Image is output here
   double *z_buffer ,   // nr * nc work area
   double *fill_buffer  // nr * nc work area
   )
{
   int i, j, k, node, first, last, color, xa, xb, ya, yb, dx, dy, row, col, final, g ;
   int inc_nochg, inc_chg, pos_slope ;
   double rate, offset, distance, first_dist, last_dist ;

   memset ( img , 255 , 3 * nr * nc ) ; // Initialize to all white
   for (i=0 ; i<nr ; i++) {             // Initialize Z buffer to infinitely distant
      for (j=0 ; j<nc ; j++) {
         z_buffer[i*nc+j] = 1.e60 ;
         fill_buffer[i*nc+j] = 1.e60 ;
         }
      }

   // We are viewing from -Y, so +X is right (x here) and +Z is up (y here)
   xa = (int) (x[0] * nc + 0.5) ;
   if (xa >= nc)
      xa = nc - 1 ;
   ya = (int) (z[0] * nr + 0.5) ;
   if (ya >= nr)
      ya = nr - 1 ;

   for (node=1 ; node<n ; node++) {

      if (node % 10 == 0) {
         setpos_wait_message ( (double) node / (double) n ) ;
         user_pressed_escape () ;
         }

      xb = (int) (x[node] * nc + 0.5) ;
      if (xb >= nc)
         xb = nc - 1 ;
      yb = (int) (z[node] * nr + 0.5) ;
      if (yb >= nr)
         yb = nr - 1 ;

      if (draw_line[node] == 0) {  // Just move without drawing?
         first_dist = last_dist = 0.0;  // Not needed.  Shuts up LINT.
         // If we just drew a filled object, this is the POP due to closing curly brace.
         // Fill in the object and then zero the fill buffer to get ready for the next.
         if (draw_line[node-1] < 0) {  // Prior node was perimeter of a filled object?
            color = -(draw_line[node-1] + 1) ;
            for (i=0 ; i<nr ; i++) {
               k = 0 ;  // Will count perimeter flags in this row
               first = last = -1 ;
               for (j=0 ; j<nc ; j++) {
                  if (fill_buffer[i*nc+j] < 1.e59) {
                     ++k ;           // Count perimeter flags
                     if (first < 0) {
                        first = j ;
                        first_dist = fill_buffer[i*nc+j] ;
                        }
                     last = j ;
                     last_dist = fill_buffer[i*nc+j] ;
                     }
                  }
               if (k < 2)
                  continue ;
               for (j=first ; j<=last ; j++) {
                  distance = first_dist + (double) (j - first) / (last - first) * (last_dist - first_dist) ;
                  if (distance < z_buffer[i*nc+j]) {
                     z_buffer[i*nc+j] = distance ;
                     if (AG3D_params.n_colors == 0) {
                        img[3*((nr-i-1)*nc+j)+0] = 0 ;
                        img[3*((nr-i-1)*nc+j)+1] = 0 ;
                        img[3*((nr-i-1)*nc+j)+2] = 0 ;
                        }
                     else {
                        img[3*((nr-i-1)*nc+j)+0] = (unsigned char) AG3D_params.colors[color][0] ;
                        img[3*((nr-i-1)*nc+j)+1] = (unsigned char) AG3D_params.colors[color][1] ;
                        img[3*((nr-i-1)*nc+j)+2] = (unsigned char) AG3D_params.colors[color][2] ;
                        }
                     } // If distance < z buffer (so visible, not hidden)
                  } // For each column from first through last
               } // Checking each row
            // Zero fill buffer to prepare for next object
            for (i=0 ; i<nr ; i++) {
               for (j=0 ; j<nc ; j++)
                  fill_buffer[i*nc+j] = 1.e60 ;
               }
            } // Prior node was perimeter of a filled object

         xa = xb ;
         ya = yb ;
         continue ;
         } // Just move due to POP

      dx = xb - xa ;
      dy = yb - ya ;

      pos_slope = (dx > 0) ;
      if (dy < 0)
         pos_slope = ! pos_slope ;

      if (abs(dx) > abs(dy)) {  // Gentle slope case

         rate = (y[node] - y[node-1]) / (dx + 1.e-20) ;  // These are for Z buffer, not line
         offset = x[node-1] ;                            // Lets us compute distance along line

         if (dx > 0) {          // We want to always use positive increment
            col = xa ;
            row = ya ;
            final = xb ;        // We are incrementing x
            } // dx > 0
         else {
            col = xb ;
            row = yb ;
            final = xa ;
            } // dx <= 0

         inc_chg = 2 * (abs(dy) - abs(dx)) ;
         inc_nochg = 2 * abs(dy) ;
         g = inc_nochg - abs(dx) ;

         if (pos_slope) {
            while (col < final) {
               distance = offset + (col - xa) * rate ;  // Relative distance from eye
               if (draw_line[node] < 0)  // Drawing a filled object with 'S'?  Flag perimeter
                  fill_buffer[row*nc+col] = distance ;   // We'll check z_buffer when object is filled
               if (distance < z_buffer[row*nc+col]) {
                  z_buffer[row*nc+col] = distance ;
                  circle ( thickness[node] , abs(draw_line[node])-1 , nr-1-row , col , nr , nc , img ) ;
                  }
               ++col ;
               if (g >= 0) {
                  ++row ;
                  g += inc_chg ;
                  } // if g >= 0
               else
                  g += inc_nochg ;
               } // while col < final
            } // pos_slope

         else {  // NOT pos_slope
            while (col <= final) {
               distance = offset + (col - xa) * rate ;  // Relative distance from eye
               if (draw_line[node] < 0)  // Drawing a filled object with 'S'?  Flag perimeter
                  fill_buffer[row*nc+col] = distance ;   // We'll check z_buffer when object is filled
               if (distance < z_buffer[row*nc+col]) {
                  z_buffer[row*nc+col] = distance ;
                  circle ( thickness[node] , abs(draw_line[node])-1 , nr-1-row , col , nr , nc , img ) ;
                  }
               ++col ;
               if (g > 0) {
                  --row ;
                  g += inc_chg ;
                  } // if g > 0
               else
                  g += inc_nochg ;
               } // while col < final
            } // NOT pos_slope
         } // (abs(dx) > abs(dy)) Gentle slope case

      else {   // Steep slope case -----------------

         rate = (y[node] - y[node-1]) / (dy + 1.e-20) ;  // These are for Z buffer, not line
         offset = y[node-1] ;                            // Lets us compute distance along line

         if (dy > 0) {          // We want to always use positive increment
            col = xa ;
            row = ya ;
            final = yb ;        // We are incrementing y
            } // dx > 0
         else {
            col = xb ;
            row = yb ;
            final = ya ;
            } // dx <= 0

         inc_chg = 2 * (abs(dx) - abs(dy)) ;
         inc_nochg = 2 * abs(dx) ;
         g = inc_nochg - abs(dy) ;

         if (pos_slope) {
            while (row < final) {
               distance = offset + (row - ya) * rate ;  // Relative distance from eye
               if (draw_line[node] < 0)  // Drawing a filled object with 'S'?  Flag perimeter
                  fill_buffer[row*nc+col] = distance ;   // We'll check z_buffer when object is filled
               if (distance < z_buffer[row*nc+col]) {
                  z_buffer[row*nc+col] = distance ;
                  circle ( thickness[node] , abs(draw_line[node])-1 , nr-1-row , col , nr , nc , img ) ;
                  }
               ++row ;
               if (g >= 0) {
                  ++col ;
                  g += inc_chg ;
                  } // if g >= 0
               else
                  g += inc_nochg ;
               } // while row < final
            } // pos_slope

         else {  // NOT pos_slope
            while (row <= final) {
               distance = offset + (row - ya) * rate ;  // Relative distance from eye
               if (draw_line[node] < 0)  // Drawing a filled object with 'S'?  Flag perimeter
                  fill_buffer[row*nc+col] = distance ;   // We'll check z_buffer when object is filled
               if (distance < z_buffer[row*nc+col]) {
                  z_buffer[row*nc+col] = distance ;
                  circle ( thickness[node] , abs(draw_line[node])-1 , nr-1-row , col , nr , nc , img ) ;
                  }
               ++row ;
               if (g > 0) {
                  --col ;
                  g += inc_chg ;
                  } // if g > 0
               else
                  g += inc_nochg ;
               } // while row < final
            } // NOT pos_slope
         } // Steep slope case

      xa = xb ;
      ya = yb ;
      } // For all nodes
}


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

   Local routine rotates the object to change the view position

   We assume +X is to the right, +Y is straight ahead, and +Z is up.
   We are looking at the object from infinitely far out the -Y axis.
   Azimuth is CCW rotation 0-360 degrees of eye around Z axis.
   Elevation is elevation 0-90 degrees of eye above the XY plane

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

static void rotate_image ( int n , double *x , double *y , double *z ,
                           double azimuth , double elevation )
{
   int i ;
   double sine_azimuth, cosine_azimuth, sine_elevation, cosine_elevation, tempval ;
   double xmin, xmax, ymin, ymax, zmin, zmax, scale ;

   azimuth *= PI / 180.0 ;    // Convert from degrees to radians
   elevation *= PI / 180.0 ;

   sine_azimuth = sin ( azimuth ) ;
   cosine_azimuth = cos ( azimuth ) ;
   sine_elevation = sin ( elevation ) ;
   cosine_elevation = cos ( elevation ) ;

   for (i=0 ; i<n ; i++) {

      // First, rotate clockwise about Z, because the eye rotates CCW
      tempval = x[i] ;
      x[i] = cosine_azimuth * x[i] + sine_azimuth * y[i] ;
      y[i] = cosine_azimuth * y[i] - sine_azimuth * tempval ;

      // Second, rotate downward about X, because the eye rises above XY plane
      tempval = y[i] ;
      y[i] = cosine_elevation * y[i] - sine_elevation * z[i] ;
      z[i] = cosine_elevation * z[i] + sine_elevation * tempval ;
      }


   xmin = ymin = zmin = 1.e30 ;  // Boundary of the figure
   xmax = ymax = zmax = -1.e30 ;

   for (i=0 ; i<n ; i++) {
      if (x[i] < xmin)
         xmin = x[i] ;
      if (x[i] > xmax)
         xmax = x[i] ;
      if (y[i] < ymin)
         ymin = y[i] ;
      if (y[i] > ymax)
         ymax = y[i] ;
      if (z[i] < zmin)
         zmin = z[i] ;
      if (z[i] > zmax)
         zmax = z[i] ;
      }

/*
   The nodes are rotated.
   Shift and scale them so that the figure touches the bottom and left
   sides, as well as either the top or right side or both.
   Recall that the image will be displayed with +X pointing right,
   +Z pointing up, and Y being the distance from the eye.
*/

   if (xmax-xmin >= zmax-zmin)
      scale = 0.98 / (xmax - xmin + 1.e-30) ;
   else
      scale = 0.98 / (zmax - zmin + 1.e-30) ;

   for (i=0 ; i<n ; i++) {
      x[i] = scale * (x[i] - xmin) + 0.01 ;
      y[i] = scale * (y[i] - ymin) + 0.01 ;
      z[i] = scale * (z[i] - zmin) + 0.01 ;
      }
}


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

   Local routine prints an operation string to log file

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

static void print_string ( int n , int *string )
{
   int i, k ;
   char msg[2048], msg2[2] ;

   if (n >= 2048)
      n = 2047 ;

   strcpy_s ( msg , "" ) ;

   for (i=0 ; i<n ; i++) {
      if (string[i] == P_PLUS)
         strcat_s ( msg , "+" ) ;
      else if (string[i] == P_MINUS)
         strcat_s ( msg , "-" ) ;
      else if (string[i] == P_PUSH)
         strcat_s ( msg , "[" ) ;
      else if (string[i] == P_POP)
         strcat_s ( msg , "]" ) ;
      else if (string[i] == P_START_BRACE)
         strcat_s ( msg , "{" ) ;
      else if (string[i] == P_END_BRACE)
         strcat_s ( msg , "}" ) ;
      else if (string[i] == P_SHRINK)
         strcat_s ( msg , "V" ) ;
      else if (string[i] == P_FORWARD)
         strcat_s ( msg , "F" ) ;
      else if (string[i] == P_STEP)
         strcat_s ( msg , "S" ) ;
      else if (string[i] == P_PITCH_DOWN)
         strcat_s ( msg , "&" ) ;
      else if (string[i] == P_PITCH_UP)
         strcat_s ( msg , "^" ) ;
      else if (string[i] == P_ROLL_LEFT)
         strcat_s ( msg , "\\" ) ;
      else if (string[i] == P_ROLL_RIGHT)
         strcat_s ( msg , "/" ) ;
      else if (string[i] == P_TURN_AROUND)
         strcat_s ( msg , "|" ) ;
      else if (string[i] == P_COLOR) {
         strcat_s ( msg , "C" ) ;
         sprintf_s ( msg2 , "%d", string[++i] ) ;
         strcat_s ( msg , msg2 ) ;
         }
      else if (string[i] == P_NAME) {
         strcat_s ( msg , "\"" ) ;
         k = string[++i] ;
         strcat_s ( msg , AG3D_params.prod_name[k] ) ;
         strcat_s ( msg , "\"" ) ;
         }
      else
         strcat_s ( msg , "ERROR" ) ; // Should never happen
      }

   audit ( msg ) ;
}


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

   This is the main routine, called externally

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


int compute_3Dag ( AG3D_Params *params , unsigned char *img )
{
   int nr, nc, i, k, n, ret_val ;
   int *source, *dest, *final, *draw_line, * stack_c;
   double *x, *y, *z, *thickness, *z_buffer ;
   double *stack_x, *stack_y, *stack_z, *stack_hx, *stack_hy, *stack_hz ;
   double *stack_ux, *stack_uy, *stack_uz, *stack_t, *fill_buffer ;
   char msg[256];

   begin_wait_message ( "Creating 3D axial generation image" , 1 ) ;
   setpos_wait_message ( 0.0 ) ;
   user_pressed_escape () ;  // Spin the message loop to display message

   source = NULL ;
   dest = NULL ;
   x = NULL ;
   y = NULL ;
   z = NULL ;
   thickness = NULL ;
   draw_line = NULL ;
   stack_x = NULL ;
   stack_y = NULL ;
   stack_z = NULL ;
   stack_hx = NULL ;
   stack_hy = NULL ;
   stack_hz = NULL ;
   stack_ux = NULL ;
   stack_uy = NULL ;
   stack_uz = NULL ;
   stack_t = NULL ;
   stack_c = NULL ;
   z_buffer = NULL ;
   fill_buffer = NULL ;

   ret_val = 0 ;
         
   nr = AG3D_params.nrows ;
   nc = AG3D_params.ncols ;
   memset ( img , 255 , 3 * nr * nc ) ;


/*
   Allocate source, z buffer, fill buffer, and stacks.
   Other allocation will be done later, on the fly as needed.
*/

   source = (int *) malloc ( params->n_axiom * sizeof(int) ) ;

   stack_x = (double *) malloc ( 1000 * sizeof(double) ) ;  // Ridiculously conservative
   stack_y = (double *) malloc ( 1000 * sizeof(double) ) ;
   stack_z = (double *) malloc ( 1000 * sizeof(double) ) ;
   stack_hx = (double *) malloc ( 1000 * sizeof(double) ) ;
   stack_hy = (double *) malloc ( 1000 * sizeof(double) ) ;
   stack_hz = (double *) malloc ( 1000 * sizeof(double) ) ;
   stack_ux = (double *) malloc ( 1000 * sizeof(double) ) ;
   stack_uy = (double *) malloc ( 1000 * sizeof(double) ) ;
   stack_uz = (double *) malloc ( 1000 * sizeof(double) ) ;
   stack_t = (double *) malloc ( 1000 * sizeof(double) ) ;
   stack_c = (int *) malloc ( 1000 * sizeof(int) ) ;

   z_buffer = (double *) malloc ( nr * nc * sizeof(double) ) ;
   fill_buffer = (double *) malloc ( nr * nc * sizeof(double) ) ;

/*
   Check for insufficient memory
*/

   if (source == NULL  ||  stack_x == NULL  ||  stack_y == NULL  ||  stack_z == NULL  ||
       stack_hx == NULL  ||  stack_hy == NULL  ||  stack_hz == NULL  ||
       stack_ux == NULL  ||  stack_uy == NULL  ||  stack_uz == NULL  ||
       stack_t == NULL  ||  stack_c == NULL  ||  z_buffer == NULL  ||  fill_buffer == NULL) {
      if (source != NULL) {
         free ( source ) ;
         source = NULL ;
         }
      if (stack_x != NULL) {
         free ( stack_x ) ;
         stack_x = NULL ;
         }
      if (stack_y != NULL) {
         free ( stack_y ) ;
         stack_y = NULL ;
         }
      if (stack_z != NULL) {
         free ( stack_z ) ;
         stack_z = NULL ;
         }
      if (stack_hx != NULL) {
         free ( stack_hx ) ;
         stack_hx = NULL ;
         }
      if (stack_hy != NULL) {
         free ( stack_hy ) ;
         stack_hy = NULL ;
         }
      if (stack_hz != NULL) {
         free ( stack_hz ) ;
         stack_hz = NULL ;
         }
      if (stack_ux != NULL) {
         free ( stack_ux ) ;
         stack_ux = NULL ;
         }
      if (stack_uy != NULL) {
         free ( stack_uy ) ;
         stack_uy = NULL ;
         }
      if (stack_uz != NULL) {
         free ( stack_uz ) ;
         stack_uz = NULL ;
         }
      if (stack_t != NULL) {
         free ( stack_t ) ;
         stack_t = NULL ;
         }
      if (stack_c != NULL) {
         free ( stack_c ) ;
         stack_c = NULL ;
         }
      if (z_buffer != NULL) {
         free ( z_buffer ) ;
         z_buffer = NULL ;
         }
      if (fill_buffer != NULL) {
         free ( fill_buffer ) ;
         fill_buffer = NULL ;
         }
      audit ( "ERROR... Insufficient memory" ) ;
      ret_val = ERROR_INSUFFICIENT_MEMORY ;
      goto FINISH ;
      }

         
/*
   Do all required productions.
   When the depth loop is finished, the final string will be in 'final', n long.

   For the first iteration (i=0) the axiom will be in 'source' and the
   output string will be in 'dest'.  On the next iteration (i=1) the roles will
   be reversed, with 'dest' the source and 'source' the destination.
   They bounce back and forth, reversing roles each time.
*/

   for (i=0 ; i<params->n_axiom ; i++)
      source[i] = params->axiom[i] ;
   final = source ;  // User might have depth=0, so be ready for that
   n = params->n_axiom ;

   for (i=0 ; i<params->depth ; i++) {
      if (i % 2 == 0) {
         k = count_advance ( params , n , source ) ;
         dest = (int *) realloc ( dest , k * sizeof(int) ) ;
         if (dest == NULL) {
            audit ( "ERROR... Insufficient memory" ) ;
            ret_val = ERROR_INSUFFICIENT_MEMORY ;
            goto FINISH ;
            }
         n = advance ( params , n , source , dest ) ;
         final = dest ;
         }
      else {
         k = count_advance ( params , n , dest ) ;
         source = (int *) realloc ( source , k * sizeof(int) ) ;
         if (source == NULL) {
            audit ( "ERROR... Insufficient memory" ) ;
            ret_val = ERROR_INSUFFICIENT_MEMORY ;
            goto FINISH ;
            }
         n = advance ( params , n , dest , source ) ;
         final = source ;
         }
      sprintf_s ( msg, "Iteration %d string length = %d", i+1, n ) ;
      audit ( msg ) ;
      if (n < 300)
         print_string ( n , final ) ;

      setpos_wait_message ( (i + 0.5) / params->depth ) ;
      if (user_pressed_escape()) {
         audit ( "User preseed ESCape" ) ;
         goto FINISH ;
         }
      }
         
/*
   Convert the operation string to turtle-move node coordinates.
   Then build the image using the turtle moves.
*/

   k = count_turtle ( n , final ) ;
   x = (double *) malloc ( k * sizeof(double) ) ;
   y = (double *) malloc ( k * sizeof(double) ) ;
   z = (double *) malloc ( k * sizeof(double) ) ;
   thickness = (double *) malloc ( k * sizeof(double) ) ;
   draw_line = (int *) malloc ( k * sizeof(int) ) ;
   if (x == NULL  ||  y == NULL  ||  z == NULL  ||  thickness == NULL  ||  draw_line == NULL) {
      audit ( "ERROR... Insufficient memory" ) ;
      ret_val = ERROR_INSUFFICIENT_MEMORY ;
      goto FINISH ;
      }
   n = turtle ( n , final , params , x , y , z , thickness , draw_line ,
                stack_x , stack_y , stack_z , stack_hx , stack_hy , stack_hz , stack_ux , stack_uy , stack_uz ,
                stack_t , stack_c ) ;

   rotate_image ( n , x , y , z , params->azimuth , params->elevation ) ; // This also rescales, so MUST be called

   build_image ( n , x , y , z , thickness , draw_line , nr , nc , img , z_buffer , fill_buffer ) ;


FINISH:
   end_wait_message () ;

   audit ( "" ) ;
   audit ( "-------------------------> Finished project <-------------------------" ) ;
   audit ( "" ) ;

   if (source != NULL)
      free ( source ) ;
   if (dest != NULL)
      free ( dest ) ;
   if (x != NULL)
      free ( x ) ;
   if (y != NULL)
      free ( y ) ;
   if (z != NULL)
      free ( z ) ;
   if (thickness != NULL)
      free ( thickness ) ;
   if (draw_line != NULL)
      free ( draw_line ) ;
   if (stack_x != NULL)
      free ( stack_x ) ;
   if (stack_y != NULL)
      free ( stack_y ) ;
   if (stack_z != NULL)
      free ( stack_z ) ;
   if (stack_hx != NULL)
      free ( stack_hx ) ;
   if (stack_hy != NULL)
      free ( stack_hy ) ;
   if (stack_hz != NULL)
      free ( stack_hz ) ;
   if (stack_ux != NULL)
      free ( stack_ux ) ;
   if (stack_uy != NULL)
      free ( stack_uy ) ;
   if (stack_uz != NULL)
      free ( stack_uz ) ;
   if (stack_t != NULL)
      free ( stack_t ) ;
   if (stack_c != NULL)
      free ( stack_c ) ;
   if (z_buffer != NULL)
      free ( z_buffer ) ;
   if (fill_buffer != NULL)
      free ( fill_buffer ) ;
   return ret_val ;
}