/******************************************************************************/
/*                                                                            */
/*  COMPUTE_ER - All routines for computing edge rewriting                    */
/*                                                                            */
/******************************************************************************/

#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"


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

   Local routine counts how many destination items will be needed

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

static int count_replace (
   int n ,            // Number of tokens in source
   int *source ,      // Source string (predecessor)
   int n_primary ,    // Number of primary production tokens (must be positive)
   int n_secondary    // Number of secondary production tokens (may be zero)
   )
{
   int i, j ;

   j = 0 ;

   for (i=0 ; i<n ; i++) {
      if (source[i] == P_PLUS)           // Just copy turns
         ++j ;
      else if (source[i] == P_MINUS)
         ++j ;
      else if (source[i] == P_PRIMARY)   // Insert a production
         j += n_primary ;
      else if (source[i] == P_SECONDARY)
         j += n_secondary ;
      }

   return j ;
}

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

   Local routine does the edge replacement for one iteration

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

static int replace (
   int n ,            // Number of tokens in source
   int *source ,      // Source string (predecessor)
   int *dest ,        // Destination string (successor)
   int n_primary ,    // Number of primary production tokens (must be positive)
   int *primary ,     // They are here
   int n_secondary ,  // Number of secondary production tokens (may be zero)
   int *secondary     // They are here
   )
{
   int i, j, k ;

   j = 0 ;

   for (i=0 ; i<n ; i++) {
      if (source[i] == P_PLUS)           // Just copy turns
         dest[j++] = P_PLUS ;
      else if (source[i] == P_MINUS)
         dest[j++] = P_MINUS ;
      else if (source[i] == P_PRIMARY) { // Insert a production
         for (k=0 ; k<n_primary ; k++)
            dest[j++] = primary[k] ;
         }
      else if (source[i] == P_SECONDARY) {
         for (k=0 ; k<n_secondary ; k++)
            dest[j++] = secondary[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 ;

   for (i=0 ; i<n ; i++) {
      if (opstr[i] == P_PRIMARY  ||  opstr[i] == P_SECONDARY)   // Take a step
         ++node ;
      }

   return node ;
}


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

   Local routine converts an operation string to turtle-move node coordinates.
   This will leave x (column) and y (row) in the range 0.01 - 0.99,
   with both x and y reaching a minimum of 0.01, and at least one of them
   reaching a maximum of 0.99.

   It assumes a starting position of (0,0), the lower left corner,
   and an angle of 0, pointing to the right.

   It returns the number of nodes in x and y.

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


static int turtle (
   int n ,          // Number of operations in 'opstr'
   int *opstr ,     // The operation string
   double angle ,   // Angle (0-90 degrees) by which to turn; + means CCW
   double *x ,      // Output of x coordinates of nodes
   double *y        // And y coordinates
   )
{
   int i, node ;
   double current_x, current_y, current_angle, xmin, xmax, ymin, ymax, scale ;
   double cos_a, sin_a ;
   double pi = 3.141592653589793 ;
   double two_pi = 2.0 * 3.141592653589793 ;
   double half_pi = 0.5 * 3.141592653589793 ;

   angle *= 3.141592653589793 / 180.0 ;  // Convert degrees to radians
   current_x = current_y = current_angle = 0.0 ;
   xmin = xmax = ymin = ymax = 0.0 ;  // Boundary of the figure

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

   for (i=0 ; i<n ; i++) {
      if (opstr[i] == P_PLUS) {        // Turn CCW
         current_angle += angle ;
         if (current_angle >= two_pi)
            current_angle -= two_pi ;
         }
      else if (opstr[i] == P_MINUS) {  // Turn CW
         current_angle -= angle ;
         if (current_angle < 0.0)
            current_angle += two_pi ;
         }
      else if (opstr[i] == P_PRIMARY  ||  opstr[i] == P_SECONDARY) {  // Take a step
         // Periodically reset angle to exact value to prevent error cumulation
         if (fabs(current_angle) < 1.e-12)
            current_angle = 0.0 ;
         else if (fabs(current_angle - half_pi) < 1.e-12)
            current_angle = half_pi ;
         else if (fabs(current_angle - pi) < 1.e-12)
            current_angle = pi ;
         else if (fabs(current_angle - pi - half_pi) < 1.e-12)
            current_angle = pi + half_pi ;
         cos_a = cos ( current_angle ) ;
         sin_a = sin ( current_angle ) ;
         // Make sure vertical/horizontal moves are perfect
         if (fabs(cos_a) < 1.e-12)
            cos_a = 0.0 ;
         else if (fabs(cos_a - 1.0) < 1.e-12)
            cos_a = 1.0 ;
         else if (fabs(cos_a + 1.0) < 1.e-12)
            cos_a = -1.0 ;
         if (fabs(sin_a) < 1.e-12)
            sin_a = 0.0 ;
         else if (fabs(sin_a - 1.0) < 1.e-12)
            sin_a = 1.0 ;
         else if (fabs(sin_a + 1.0) < 1.e-12)
            sin_a = -1.0 ;
         current_x += cos_a ;
         current_y += sin_a ;
         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 ;
         x[node] = current_x ;
         y[node++] = current_y ;
         }
      }

/*
   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)
      scale = 0.98 / (xmax - xmin + 1.e-30) ;
   else
      scale = 0.98 / (ymax - ymin + 1.e-30) ;

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

   return node ;
}


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

   Local routine called from build_image() draws a black circle

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

static void circle (
   int diameter ,     // Diameter of the circle
   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 = 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)
         img[row*nc+i] = 0 ;
      }

   // 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)
            img[(row-j)*nc+i] = 0 ;
         if (i >= 0  &&  i < nc  &&  row+j < nr)
            img[(row+j)*nc+i] = 0 ;
         }
      }
}


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

   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
   int width ,          // Line width, at least 1, 3 often best
   int nr ,             // Image dimensions
   int nc ,
   unsigned char *img   // Image is output here
   )
{
   int node, xa, xb, ya, yb, dx, dy, row, col, final, g ;
   int inc_nochg, inc_chg, pos_slope ;

   memset ( img , 255 , nr * nc ) ; // Initialize to all white

   xa = (int) (x[0] * nc + 0.5) ;
   if (xa >= nc)
      xa = nc - 1 ;
   ya = (int) (y[0] * nr + 0.5) ;
   if (ya >= nr)
      ya = nr - 1 ;

   for (node=1 ; node<n ; node++) {
      xb = (int) (x[node] * nc + 0.5) ;
      if (xb >= nc)
         xb = nc - 1 ;
      yb = (int) (y[node] * nr + 0.5) ;
      if (yb >= nr)
         yb = nr - 1 ;

      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
         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) {
               circle ( width , 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) {
               circle ( width , 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 -----------------
         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) {
               circle ( width , 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 {  // pos_slope
            while (row <= final) {
               circle ( width , 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 prints an operation string to log file

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

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

   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_PRIMARY)
         strcat_s ( msg , "P" ) ;
      else if (string[i] == P_SECONDARY)
         strcat_s ( msg , "S" ) ;
      else
         strcat_s ( msg , "!" ) ;
      }

   audit ( msg ) ;
}


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

   This is the main routine, called externally

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


int compute_er ( ER_Params *params , unsigned char *img )
{
   int nr, nc, i, k, n, n_primary, n_secondary, ret_val ;
   int *primary_production, *secondary_production, *source, *dest, *final ;
   double *x, *y ;
   
   ret_val = 0 ;
   n_primary = (int) strlen ( params->primary ) ;
   n_secondary = (int) strlen ( params->secondary ) ;
         
   nr = ER_params.nrows ;
   nc = ER_params.ncols ;
   memset ( img , 255 , nr * nc ) ;

   primary_production = NULL ;
   secondary_production = NULL ;
   source = NULL ;
   dest = NULL ;
   x = NULL ;
   y = NULL ;


/*
   Allocate the production recoding vectors and the initial source vector.
   Other memory allocation will be done later, on the fly.
*/

   n = (int) strlen ( ER_params.axiom ) ;
   source = (int *) malloc ( n * sizeof(int) ) ;
   primary_production = (int *) malloc ( n_primary * sizeof(int) ) ;
   if (n_secondary)
      secondary_production = (int *) malloc ( n_secondary * sizeof(int) ) ;
   else
      secondary_production = NULL ;

   if (source == NULL  || primary_production == NULL  || (n_secondary && secondary_production == NULL)) {
      if (source != NULL) {
         free ( source ) ;
         source = NULL ;
         }
      if (primary_production != NULL) {
         free ( primary_production ) ;
         primary_production = NULL ;
         }
      if (secondary_production != NULL) {
         free ( secondary_production ) ;
         secondary_production = NULL ;
         }
      audit ( "ERROR... Insufficient memory" ) ;
      ret_val = ERROR_INSUFFICIENT_MEMORY ;
      goto FINISH ;
      }

   begin_wait_message ( "Creating edge rewriting image" , 1 ) ;
   user_pressed_escape () ;  // Spin the message loop to display message

/*
   Recode the productions
*/

   for (i=0 ; i<n_primary ; i++) {
      if (params->primary[i] == '+')
         primary_production[i] = P_PLUS ;
      else if (params->primary[i] == '-')
         primary_production[i] = P_MINUS ;
      else if (params->primary[i] == 'P')
         primary_production[i] = P_PRIMARY ;
      else if (params->primary[i] == 'S')
         primary_production[i] = P_SECONDARY ;
      else {
         audit ( "Invalid character in primary production" ) ;
         ret_val = ERROR_INVALID_PRODUCTION_OPERATION ;
         goto FINISH ;
         }
      }

   for (i=0 ; i<n_secondary ; i++) {
      if (params->secondary[i] == '+')
         secondary_production[i] = P_PLUS ;
      else if (params->secondary[i] == '-')
         secondary_production[i] = P_MINUS ;
      else if (params->secondary[i] == 'P')
         secondary_production[i] = P_PRIMARY ;
      else if (params->secondary[i] == 'S')
         secondary_production[i] = P_SECONDARY ;
      else {
         audit ( "Invalid character in secondary production" ) ;
         ret_val = ERROR_INVALID_PRODUCTION_OPERATION ;
         goto FINISH ;
         }
      }
         
/*
   Recode the axiom and place it in 'source'
*/

   n = (int) strlen(params->axiom) ;
   for (i=0 ; i<n ; i++) {
      if (params->axiom[i] == '+')
         source[i] = P_PLUS ;
      else if (params->axiom[i] == '-')
         source[i] = P_MINUS ;
      else if (params->axiom[i] == 'P')
         source[i] = P_PRIMARY ;
      else if (params->axiom[i] == 'S')
         source[i] = P_SECONDARY ;
      else {
         audit ( "Invalid character in axiom" ) ;
         ret_val = ERROR_INVALID_PRODUCTION_OPERATION ;
         goto FINISH ;
         }
      }
         
/*
   If a secondary production is mentioned, there had better be one!
*/

   if (n_secondary == 0) {
      n = (int) strlen(params->axiom) ;
      for (i=0 ; i<n ; i++) {
         if (source[i] == P_SECONDARY) {
            ret_val = ERROR_REFERENCE_NO_SECONDARY ;
            goto FINISH ;
            }
         }
      for (i=0 ; i<n_primary ; i++) {
         if (primary_production[i] == P_SECONDARY) {
            ret_val = ERROR_REFERENCE_NO_SECONDARY ;
            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.
*/

   final = source ;  // User might have depth=0, so be ready for that
   n = (int) strlen(params->axiom) ;
   for (i=0 ; i<params->depth ; i++) {
      if (i % 2 == 0) {
         k = count_replace ( n , source , n_primary , n_secondary ) ;
         dest = (int *) realloc ( dest , k * sizeof(int) ) ;
         if (dest == NULL) {
            audit ( "ERROR... Insufficient memory" ) ;
            ret_val = ERROR_INSUFFICIENT_MEMORY ;
            goto FINISH ;
            }
         n = replace ( n , source , dest , n_primary , primary_production , n_secondary , secondary_production ) ;
         final = dest ;
         }
      else {
         k = count_replace ( n , dest , n_primary , n_secondary ) ;
         source = (int *) realloc ( source , k * sizeof(int) ) ;
         if (source == NULL) {
            audit ( "ERROR... Insufficient memory" ) ;
            ret_val = ERROR_INSUFFICIENT_MEMORY ;
            goto FINISH ;
            }
         n = replace ( n , dest , source , n_primary , primary_production , n_secondary , secondary_production ) ;
         final = source ;
         }
      if (i == 0) {  // This is optional, purely for user's edification
         audit ( "String after first iteration" ) ;
         print_string ( n , final ) ;
         }
      if (i == 1  &&  n < 200) {  // This is optional, purely for user's edification
         audit ( "String after second iteration" ) ;
         print_string ( n , final ) ;
         }
      else if (i == 2  &&  n < 200) {  // This is optional, purely for user's edification
         audit ( "String after third iteration" ) ;
         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) ) ;
   if (x == NULL  ||  y == NULL) {
      audit ( "ERROR... Insufficient memory" ) ;
      ret_val = ERROR_INSUFFICIENT_MEMORY ;
      goto FINISH ;
      }
   n = turtle ( n , final , params->angle , x , y ) ;
   build_image ( n , x , y , params->width , nr , nc , img ) ;


FINISH:
   end_wait_message () ;

   free ( primary_production ) ;
   if (secondary_production != NULL)
      free ( secondary_production ) ;
   if (source != NULL)
      free ( source ) ;
   if (dest != NULL)
      free ( dest ) ;
   if (x != NULL)
      free ( x ) ;
   if (y != NULL)
      free ( y ) ;
   return ret_val ;
}
