/***********************************************************
 *
 *  pgmblob.c - a grayscale morphogenesis tool
 *
 *  Mark J. Stock, March 27, 1998, version 0.1
 *
 ***********************************************************/


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "const.h"

#define MAX_STEPS 100

int inheight, inwidth;		/* size of the input file */

void main(int argc,char **argv) {

   int i,j,step;
   double delsqa,delsqb,absq;
   double F,k,Da,Db;
   int num_steps = 50;
   int readfile = 1;
   int randomize = 1;
   int use_subtract = 0;
   char infile[80];
   char subfile[80];
   char outfile[80];
   char progname[80];

   /* a is white, b is black */
   unsigned char a[MAX_SIZE][MAX_SIZE],b[MAX_SIZE][MAX_SIZE];
   int new_a[MAX_SIZE][MAX_SIZE],new_b[MAX_SIZE][MAX_SIZE];

   /* function declarations */
   unsigned char flatten(int val);

   F = 0.06;			/* constants of the morphogenetic algo */
   k = 0.055;
   Da = 0.1;
   Db = 0.1;
   inheight = 256;		/* default values */
   inwidth = 256;
   if (argc > 1) num_steps = atoi(argv[1]);

   /*
    * Parse command-line args
    */
   (void) strcpy(progname,argv[0]);
   if (argc < 4) (void) Usage(progname,0);
   for (i=1; i<argc; i++) {
      if (strncmp(argv[i], "-help", 2) == 0)
         (void) Usage(progname,0);
      else if (strncmp(argv[i], "-infile", 2) == 0) {
         (void) strcpy(infile,argv[++i]); 
         randomize = 0; }
      else if (strncmp(argv[i], "-outfile", 2) == 0)
         (void) strcpy(outfile,argv[++i]); 
      else if (strncmp(argv[i], "-subtract", 2) == 0) {
         (void) strcpy(subfile,argv[++i]);
         use_subtract = 1; }
      else if (strncmp(argv[i], "-n", 2) == 0)
         num_steps = atoi(argv[++i]); 
      else
         (void) Usage(progname,0);
   }


   /* initially, set all values to zero */
   for (i=0; i<MAX_SIZE; i++) {
      for (j=0; j<MAX_SIZE; j++) {
         a[i][j] = 255;
         b[i][j] = 200;
      }
   }

   /* read in the input file */
   if (readfile) (void)read_pgm_char(infile,a);
   for (i=0; i<inheight; i++) {
      for (j=0;j<inwidth;j++) {
         /* a[i][j] = a[i][j]*(0.3+0.7*rand()/(RAND_MAX+1.0)); */
      }
   }

   /* or, populate the arrays with random values */
   if (randomize) {
      for (i=0; i<inheight; i++) {
         for (j=0;j<inwidth;j++) {
            a[i][j] = a[i][j]*rand()/(RAND_MAX+1.0);
         }
      }
   }

   /* populate the enemy cells randomly */
   for (i=0; i<inheight; i++) {
      for (j=0;j<inwidth;j++) {
         b[i][j] = b[i][j]*rand()/(RAND_MAX+1.0);
      }
   }

   /* then, iterate through steps */
   for (step=0; step<num_steps; step++) {

      /* calculate new values from old values */
      for (i=0;i<inheight;i++) {
         for (j=0;j<inwidth;j++) {
            if (i==0) {
               if (j==0) {
                  delsqa = Da*(a[i][j+1]+a[i+1][j]-2*a[i][j]);
                  delsqb = Db*(b[i][j+1]+b[i+1][j]-2*b[i][j]);
               } else if (j==inwidth-1) {
                  delsqa = Da*(a[i][j-1]+a[i+1][j]-2*a[i][j]);
                  delsqb = Db*(b[i][j-1]+b[i+1][j]-2*b[i][j]);
               } else {
                  delsqa = Da*(a[i][j+1]+a[i][j-1]+a[i+1][j]-3*a[i][j]);
                  delsqb = Db*(b[i][j+1]+b[i][j-1]+b[i+1][j]-3*b[i][j]);
               }
            } else if (i==inheight-1) {
               if (j==0) {
                  delsqa = Da*(a[i][j+1]+a[i-1][j]-2*a[i][j]);
                  delsqb = Db*(b[i][j+1]+b[i-1][j]-2*b[i][j]);
               } else if (j==inwidth-1) {
                  delsqa = Da*(a[i][j-1]+a[i-1][j]-2*a[i][j]);
                  delsqb = Db*(b[i][j-1]+b[i-1][j]-2*b[i][j]);
               } else {
                  delsqa = Da*(a[i][j+1]+a[i][j-1]+a[i-1][j]-3*a[i][j]);
                  delsqb = Db*(b[i][j+1]+b[i][j-1]+b[i-1][j]-3*b[i][j]);
               }
            } else if (j==0) {
               delsqa = Da*(a[i][j+1]+a[i+1][j]+a[i-1][j]-3*a[i][j]);
               delsqb = Db*(b[i][j+1]+b[i+1][j]+b[i-1][j]-3*b[i][j]);
            } else if (j==inwidth-1) {
               delsqa = Da*(a[i][j-1]+a[i+1][j]+a[i-1][j]-3*a[i][j]);
               delsqb = Db*(b[i][j-1]+b[i+1][j]+b[i-1][j]-3*b[i][j]);
            } else {
               delsqa = Da*(a[i][j+1]+a[i][j-1]+a[i+1][j]+a[i-1][j]-4*a[i][j]);
               delsqb = Db*(b[i][j+1]+b[i][j-1]+b[i+1][j]+b[i-1][j]-4*b[i][j]);
            }
            /* competition */
            /* absq = a[i][j]*b[i][j]*b[i][j];
            new_a[i][j] = a[i][j] + (delsqa + F*(1.0-a[i][j]))*step/10.0;
            new_b[i][j] = b[i][j] + (delsqb - b[i][j]*(F+k))*step/10.0; */
            /* overcrowding */
            new_a[i][j] = a[i][j] + (delsqa + a[i][j]/10.0 - a[i][j]*b[i][j]/500.0);
            new_b[i][j] = b[i][j] + (delsqb + b[i][j]/10.0 - b[i][j]*a[i][j]/500.0);
         }
      }

      /* replace the old values with the new values */
      for (i=0;i<inheight;i++) {
         for (j=0;j<inwidth;j++) {
            a[i][j] = flatten(new_a[i][j]);
            b[i][j] = flatten(new_b[i][j]);
         }
      }
   }

   /* write out a frame at the end */

   /* use the subtract file to clean areas of white */
   if (use_subtract) {
      (void)read_pgm_char(subfile,b);
      for (i=0; i<inheight; i++) {
         for (j=0;j<inwidth;j++) {
            if (a[i][j] > b[i][j]) a[i][j] = a[i][j] - b[i][j];
            else a[i][j] = 0;
         }
      }
   }
   (void) write_pgm_char(outfile,a,inheight,inwidth);

}


/*
 * Flatten an int down to 0-255 for PGM writing
 */
unsigned char flatten(int val) {
   if (val>255) return 255;
   else if (val<0) return 0;
   else return (unsigned char)val;
}


/*
 * This function writes basic usage information to stderr,
 * and then quits. Too bad.
 */
int Usage(char progname[80],int status) {

   /* Usage for pgmblob */
   static char **cpp, *help_message[] =
   {
       "where options include:",
       "  -help          writes this help information",
       "  -in name       name of input PGM file, binary or raw, max value of 255",
       "  -sub name      name of subtract file, same range, 0-255, binary",
       "  -out name      name of output file, same range, 0-255, binary",
       "  -n val         number of cycles, def= 50",
       " ",
       "Options may be abbreviated to an unambiguous length.",
       NULL
   };

   fprintf(stderr, "usage:\n  %s [-options ...]\n\n", progname);
   for (cpp = help_message; *cpp; cpp++)
      fprintf(stderr, "%s\n", *cpp);
   exit(status);
   return(0);
}

