/*
 * pipes.c
 *
 * v0.1  1998-08-04  First working version
 * v0.2  1998-08-10  Allow curved secions, Radiance output
 *
 * Copyright 1998, Mark J. Stock
 *
 * compile with "cc -o pipes pipes.c"
 *
 * Pipes was originally a Java program to create modelling code for
 * a series of non-intersecting, regular-grid pipes. It was slow
 * and I don't care to support it. Plus, I was a poor programmer
 * back then. I'm not too much better, but this is C, right?
 *
 * Input is in 2 groups: the non-dimensional quantities, such as
 * voxel size of compute volume, active sides, number of pipes,
 * minimum pipe length, etc; and the dimensional side: voxel size
 * in world units and ratio of pipe diameter to voxel size
 *
 * TODO:
 *   force pipe to end at valid sides, not just when it can't go
 *     anywhere
 *   read in and write out connection points along the sides
 *   write radiance statements, either as xform, instance, or the
 *     actual geometry statements themselves
 *   add option for input of random seed
 *
 */


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define TRUE 1
#define FALSE 0
#define BOOL unsigned char

#define MAXDIM 100
#define MAXDEPTH 20

typedef struct vector_record {
   double x;
   double y;
   double z;
} VEC;

/* for grid sizes larger than 256x256x256, you'll need to change
 * the #define's, above, too. This is not reccommended */
typedef struct cell_record *cell_ptr;
typedef struct cell_record {
   unsigned char x;
   unsigned char y;
   unsigned char z;
   cell_ptr previous;
} CELL;
CELL size = {0, 0, 0, NULL};
cell_ptr pipe_tail;

typedef struct side_record {
   int x_low;
   int x_high;
   int y_low;
   int y_high;
   int z_low;
   int z_high;
} SIDE;

/* Validity of each side for pipe start */
SIDE valid_start = {FALSE, FALSE, FALSE, FALSE, TRUE, FALSE};

/* Area of each valid starting side */
SIDE valid_area = {0, 0, 0, 0, 0, 0};

/* The cell record, note, the primary dimensions are x and y */
/* FALSE means the cell is taken, TRUE means it is free */
BOOL v[MAXDIM][MAXDIM][MAXDEPTH];

/* the world size of the objects */
double cs = 1.0;	/* cell size */
VEC world = {1.0, 1.0, 1.0};	/* cell size, world dimensions */
double r1 = 0.5;		/* major radius */
double r2 = 0.2;		/* minor radius */

void find_midpoint(CELL *, CELL *, VEC *);
void find_center(CELL *, VEC *);
void write_curve(FILE *,VEC *,VEC *,VEC *,int);

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

   int i, j, k, count, temp_i;
   char progname[80];
   int num_tries = 0;
   int pipe_count = 0;
   int seg_count = 0;
   int min_length = 10;
   int max_length = 10000;
   int max_pipes = 1;
   char outfile[80];
   int use_pgm = FALSE;
   double straight_frac = 0.7;
   BOOL temp_b;
   VEC temp_v;
   CELL temp_l;
   cell_ptr current,next_seg;
   FILE *outfp;

   void set_areas();
   void mark_off(CELL *);
   int choose_next_cell(CELL *,CELL *,double);
   void write_pipe(CELL *,FILE *,int);
   void write_pgm(int);
   void free_mem(CELL *);

   /* the solution space, voxel units, default size */
   (void) strcpy(outfile,"d");
   size.x = 20;
   size.y = 20;
   size.z = 5;

   /* Parse command-line args */
   (void) strcpy(progname,argv[0]);
   if (argc < 3) (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], "-outfile", 2) == 0)
         (void) strcpy(outfile,argv[++i]); 
      else if (strncmp(argv[i], "-s", 2) == 0) {
         size.x = atoi(argv[++i]); 
         size.y = atoi(argv[++i]); 
         size.z = atoi(argv[++i]); }
      else if (strncmp(argv[i], "-n", 2) == 0)
         max_pipes = atoi(argv[++i]); 
      else if (strncmp(argv[i], "-l", 2) == 0) {
         min_length = atoi(argv[++i]); 
         max_length = atoi(argv[++i]); }
      else if (strncmp(argv[i], "-p", 2) == 0)
         use_pgm = TRUE;
      else if (strncmp(argv[i], "-f", 2) == 0)
         straight_frac = atof(argv[++i]);
      else
         (void) Usage(progname,0);
   }

   /* Check input and act on user-defined options */
   if (size.x > MAXDIM || size.y > MAXDIM || size.z > MAXDEPTH) {
      fprintf(stderr,"Maximum dimensions are %d %d %d (in x,y,z)\n\n",MAXDIM,MAXDIM,MAXDEPTH);
      exit(0);
   }

   /* the solution space, voxel units */
   (void) set_areas();

   /* clean the solution space */
   for (i=0; i<size.x; i++) {
      for (j=0; j<size.y; j++) {
         for (k=0; k<size.z; k++) {
            v[i][j][k] = TRUE;
   }  }  }


   /* Open the output file */
   outfp = fopen(outfile,"w");
   if (outfp==NULL) {
      fprintf(stderr,"Could not open output file %s\n",outfile);
      exit(0);
   }


   /* Loop until there is no space left, or all pipes have been drawn */
   while (pipe_count < max_pipes) {

      /* determine starting point */
      current = (CELL *)malloc(sizeof(CELL));
      /* fprintf(stderr,"\nStarting new pipe\n"); */
      if (find_start_cell(current) == 0) {
         fprintf(stderr,"Could not find valid starting point, stopping after %d pipes.\n",pipe_count);
         free(current);
         break;
      } else 
         (void) mark_off(current);

      /* loop while pipe does not end */
      seg_count = 1;
      while (seg_count < max_length) {

         /* determine which way the segment goes */
         next_seg = (CELL *)malloc(sizeof(CELL));
         if (choose_next_cell(current,next_seg,straight_frac) == FALSE) {
            free(next_seg);
            break;
         }

         /* add the cell to the linked list for the pipe */
         next_seg->previous = current;

         /* mark the new current cell as used */
         mark_off(next_seg);

         /* reset current */
         current = next_seg;

         seg_count++;
      }

      /* if its big enough, write it, start at current and work back */
      if (seg_count >= min_length) {
         fprintf(stderr,"Pipe %d is %d segments long, writing...\n",pipe_count,seg_count);
         write_pipe(current,outfp,pipe_count);
         free_mem(current);
         num_tries = 0;
      } else {
         free_mem(current);
         if (num_tries++ > 100) break;
         continue;
      }

      pipe_count++;
   }


   /* Write stats and exit */
   if (use_pgm == TRUE) write_pgm(0);
   exit(0);
}


/*
 * Calculate the areas of each valid starting side
 */
void set_areas() {

   if (valid_start.x_low)
      valid_area.x_low = size.y*size.z;
   if (valid_start.x_high)
      valid_area.x_high = size.y*size.z;
   if (valid_start.y_low)
      valid_area.y_low = size.x*size.z;
   if (valid_start.y_high)
      valid_area.y_high = size.x*size.z;
   if (valid_start.z_low)
      valid_area.z_low = size.y*size.x;
   if (valid_start.z_high)
      valid_area.z_high = size.y*size.x;

   /* fprintf(stderr,"Areas are %d %d %d %d %d %d\n",valid_area.x_low,valid_area.x_high,
         valid_area.y_low,valid_area.y_high,valid_area.z_low,valid_area.z_high); */
}


/*
 * Determine an acceptable starting cell: must be on a valid
 * boundary (per valid_start), and must not be taken
 */
int find_start_cell(CELL *startcell) {

   int area_sum, area_tag;
   int num_tries = 0;

   /* being the first cell, it has no previous cell */
   startcell->previous = NULL;

   while (num_tries < 100) {

      /* first, choose which side to start on */
      area_sum = valid_area.x_low + valid_area.x_high + valid_area.y_low + 
         valid_area.y_high + valid_area.z_low + valid_area.z_high;
      area_tag = (int)((double)area_sum * rand()/(RAND_MAX+1.0));
      /* fprintf(stderr,"starting cell...%d / %d\n",area_tag,area_sum); */

      /* then, choose a starting cell based on that side */
      if (area_tag < valid_area.x_low) {
         startcell->x = 0;
         startcell->y = area_tag/size.z;
         startcell->z = area_tag%size.z;
      }
      area_tag -= valid_area.x_low;
      if (area_tag < valid_area.x_high) {
         startcell->x = size.x-1;
         startcell->y = area_tag/size.z;
         startcell->z = area_tag%size.z;
      }
      area_tag -= valid_area.x_high;
      if (area_tag < valid_area.y_low) {
         startcell->x = area_tag/size.x;
         startcell->y = 0;
         startcell->z = area_tag%size.x;
      }
      area_tag -= valid_area.y_low;
      if (area_tag < valid_area.y_high) {
         startcell->x = area_tag/size.x;
         startcell->y = size.y-1;
         startcell->z = area_tag%size.x;
      }
      area_tag -= valid_area.y_high;
      if (area_tag < valid_area.z_low) {
         startcell->x = area_tag/size.x;
         startcell->y = area_tag%size.x;
         startcell->z = 0;
      }
      area_tag -= valid_area.y_low;
      if (area_tag < valid_area.z_high) {
         startcell->x = area_tag/size.x;
         startcell->y = area_tag%size.x;
         startcell->z = size.z-1;
      }
      /* fprintf(stderr,"area_tag is %d, loc is %d %d %d\n",area_tag,startcell->x,startcell->y,startcell->z); */

      /* return succesful */
      if (v[startcell->x][startcell->y][startcell->z] == TRUE) return(1);

      num_tries++;
   }

   /* return unsuccesful */
   return(0);
}


/*
 * Mark this cell as taken, no other pipe can go thru it
 */
void mark_off(CELL *cell) {

   v[cell->x][cell->y][cell->z] = FALSE;
   /* fprintf(stderr,"   marked off cell %d %d %d\n",cell->x,cell->y,cell->z); */

}


/*
 * This one is pretty complicated: cell is a pointer to the
 * current end segment. The location stored there is not needed
 * any more, so the location of the new segment will replace 
 * it. CELL is passed by reference, see?
 */
int choose_next_cell(CELL *cell, CELL *new_cell, double frac) {

   int ct = 0;		/* number of turn cells within bounds and unoccupied, max=4 */
   int choose;		/* which turn cell is next */
   int temp_x, temp_y, temp_z;
   CELL straight, turn[4];
   int go_straight = FALSE;
   int go_turn[4] = {FALSE, FALSE, FALSE, FALSE};

   /* fprintf(stderr,"Find next cell\n"); */

   /* find all valid next cells, start with the straight option */
   if (cell->previous == NULL) {
      /* we're finding the 2nd cell, straight is away from the wall */
      /* fprintf(stderr,"   find 2nd valid straight cell\n"); */
      if (cell->x == 0) {
         straight.x = 1;
         straight.y = cell->y;
         straight.z = cell->z;
      } else if (cell->x == size.x-1) {
         straight.x = size.x-2;
         straight.y = cell->y;
         straight.z = cell->z;
      } else if (cell->y == 0) {
         straight.x = cell->x;
         straight.y = 1;
         straight.z = cell->z;
      } else if (cell->y == size.y-1) {
         straight.x = cell->x;
         straight.y = size.y-2;
         straight.z = cell->z;
      } else if (cell->z == 0) {
         straight.x = cell->x;
         straight.y = cell->y;
         straight.z = 1;
      } else if (cell->z == size.z-1) {
         straight.x = cell->x;
         straight.y = cell->y;
         straight.z = size.z-2;
      }
      go_straight = is_cell_in_bounds(straight.x,straight.y,straight.z);
      if (go_straight) go_straight = v[straight.x][straight.y][straight.z];
   } else {
      /* we're finding the next cell, straight is away from previous */
      /* fprintf(stderr,"   find next valid straight cell, this= %d %d %d, last= %d %d %d\n", */
         /* cell->x,cell->y,cell->z,cell->previous->x,cell->previous->y,cell->previous->z); */
      temp_x = 2*cell->x - cell->previous->x;
      temp_y = 2*cell->y - cell->previous->y;
      temp_z = 2*cell->z - cell->previous->z;
      go_straight = FALSE;
      if (is_cell_in_bounds(temp_x,temp_y,temp_z) == TRUE) {
         straight.x = temp_x;
         straight.y = temp_y;
         straight.z = temp_z;
         go_straight = v[straight.x][straight.y][straight.z];
      }
   }

   /* now, fill the turn arrays */
   if (cell->previous == NULL) {
      /* do the 4 possible cells around the first cell */
      /* fprintf(stderr,"   find 4 turn cells off first cell\n"); */
      if (cell->z > 0 && cell->z < size.z-1) {
         /* 2 possible cells are +-1 in z-direction */
         if (is_cell_in_bounds(cell->x,cell->y,cell->z + 1) == TRUE) {
            turn[ct].x = cell->x;
            turn[ct].y = cell->y;
            turn[ct].z = cell->z + 1;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
         if (is_cell_in_bounds(cell->x,cell->y,cell->z - 1) == TRUE) {
            turn[ct].x = cell->x;
            turn[ct].y = cell->y;
            turn[ct].z = cell->z - 1;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
      }
      if (cell->y > 0 && cell->y < size.y-1) {
         /* 2 possible cells are +-1 in y-direction */
         if (is_cell_in_bounds(cell->x,cell->y + 1,cell->z) == TRUE) {
            turn[ct].x = cell->x;
            turn[ct].y = cell->y + 1;
            turn[ct].z = cell->z;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
         if (is_cell_in_bounds(cell->x,cell->y - 1,cell->z) == TRUE) {
            turn[ct].x = cell->x;
            turn[ct].y = cell->y - 1;
            turn[ct].z = cell->z;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
      }
      if (cell->x > 0 && cell->x < size.x-1) {
         /* 2 possible cells are +-1 in x-direction */
         if (is_cell_in_bounds(cell->x + 1,cell->y,cell->z) == TRUE) {
            turn[ct].x = cell->x + 1;
            turn[ct].y = cell->y;
            turn[ct].z = cell->z;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
         if (is_cell_in_bounds(cell->x - 1,cell->y,cell->z) == TRUE) {
            turn[ct].x = cell->x - 1;
            turn[ct].y = cell->y;
            turn[ct].z = cell->z;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
      }
      /* fprintf(stderr,"   %d turn cells in bounds and unoccupied\n",ct); */
   } else {
      /* find the 4 arbitrary possible cells */
      /* fprintf(stderr,"   find 4 cells off current cell\n"); */
      if (cell->x - cell->previous->x == 0) {
         /* 2 possible cells are +-1 in x-direction */
         if (is_cell_in_bounds(cell->x + 1,cell->y,cell->z) == TRUE) {
            turn[ct].x = cell->x + 1;
            turn[ct].y = cell->y;
            turn[ct].z = cell->z;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
         if (is_cell_in_bounds(cell->x - 1,cell->y,cell->z) == TRUE) {
            turn[ct].x = cell->x - 1;
            turn[ct].y = cell->y;
            turn[ct].z = cell->z;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
      }
      if (cell->y - cell->previous->y == 0) {
         /* 2 possible cells are +-1 in y-direction */
         if (is_cell_in_bounds(cell->x,cell->y + 1,cell->z) == TRUE) {
            turn[ct].x = cell->x;
            turn[ct].y = cell->y + 1;
            turn[ct].z = cell->z;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
         if (is_cell_in_bounds(cell->x,cell->y - 1,cell->z) == TRUE) {
            turn[ct].x = cell->x;
            turn[ct].y = cell->y - 1;
            turn[ct].z = cell->z;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
      }
      if (cell->z - cell->previous->z == 0) {
         /* 2 possible cells are +-1 in z-direction */
         if (is_cell_in_bounds(cell->x,cell->y,cell->z + 1) == TRUE) {
            turn[ct].x = cell->x;
            turn[ct].y = cell->y;
            turn[ct].z = cell->z + 1;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
         if (is_cell_in_bounds(cell->x,cell->y,cell->z - 1) == TRUE) {
            turn[ct].x = cell->x;
            turn[ct].y = cell->y;
            turn[ct].z = cell->z - 1;
            if (v[turn[ct].x][turn[ct].y][turn[ct].z] == TRUE) ct++;
         }
      }
      /* fprintf(stderr,"   %d turn cells in bounds and unoccupied\n",ct); */
   }

   /* if there are no available cells, return zero */
   if (go_straight == FALSE && ct == 0) return(0);

   /* select to go straight or turn */
   if (go_straight == TRUE && ((double)rand()/(RAND_MAX+1.0) < frac || ct == 0)) {
      /* fprintf(stderr,"   chose straight direction\n"); */
      /* fprintf(stderr,"s"); */
      new_cell->x = straight.x;      
      new_cell->y = straight.y;      
      new_cell->z = straight.z;      
   } else {
      /* choose one of the valid turn cells */
      if (ct == 1) choose = 0;
      else choose = (int)((double)ct * (double)rand()/(RAND_MAX+1.0));
      /* fprintf(stderr,"."); */
      /* fprintf(stderr,"   chose cell %d\n",choose); */
      new_cell->x = turn[choose].x;
      new_cell->y = turn[choose].y;
      new_cell->z = turn[choose].z;
   }

   return(TRUE);
}


/*
 * Returns TRUE (1) if the cell indexes point to a cell within
 * the bounds defined by the user
 */
int is_cell_in_bounds(int x, int y, int z) {

   if (x>-1 && x<size.x && y>-1 && y<size.y && z>-1 && z<size.z) return(TRUE);
   else return(FALSE);
}


/*
 * Write one pipe, start at its tail and work thru to the
 * first segment, try to detect long straight segments, in
 * order to reduce the load on the renderer (nah.)
 */
void write_pipe(CELL *end, FILE *fp, int pnum) {

   int seg_ct = 0;
   cell_ptr curr, next, last;
   VEC *p1, *p2, *p3;
   int curve_detail = 5;
   /*double r1;		/* major radius (of turns) */
   /*double r2;		/* minor radius (of pipe segments) */

   curr= end;
   next= end->previous;
   last= NULL;

   p1 = (VEC *)malloc(sizeof(VEC));
   p2 = (VEC *)malloc(sizeof(VEC));
   p3 = (VEC *)malloc(sizeof(VEC));

   fprintf(fp,"\n# starting pipe number %d\n\n",pnum);
   fflush(fp);

   /* there are 12 different permutaions of a turn, 3 of a straight; 
    * the 2 end segments need to run right, flat up against the outer
    * boundary; but, pipes that end in a dead end need to be capped
    * with a sphere, or something special. Plus, we should try to 
    * simplify long straight sections with single cylinders (later) */


   /* while end is still a cell */
   while (curr) {

      /* is this segment one of the ends? */
      if (next== NULL) {
         /* this segment is the head, first cell created */
         /* it is always at an edge */
         /* fprintf(fp,"# last cell written\n"); */
      } else if (last== NULL) {
         /* this segment is the tail, last cell created */
         /* find if the dead end is internal or on an edge */
         /* fprintf(fp,"# first cell written\n"); */
         (void) find_center(curr,p1);
         (void) find_midpoint(curr,next,p2);
         fprintf(fp,"c%d sphere p%ds%d\n0 0 4 %lf %lf %lf  %lf\n",
            pnum,pnum,seg_ct,p1->x,p1->y,p1->z,r2);
         if (curve_detail == 1) fprintf(fp,"c%d sphere p%ds%ds\n0 0 4 %lf %lf %lf  %lf\n",
            pnum,pnum,seg_ct,p2->x,p2->y,p2->z,r2);
         fprintf(fp,"c%d cylinder p%ds%dc\n0 0 7 %lf %lf %lf  %lf %lf %lf  %lf\n",
            pnum,pnum,seg_ct,p1->x,p1->y,p1->z,p2->x,p2->y,p2->z,r2);
         if (curr->z == 0) {
         } else if (curr->z == size.z-1) {
         } else if (curr->y == 0) {
         } else if (curr->y == size.y-1) {
         } else if (curr->x == 0) {
         } else if (curr->x == size.x-1) {
         } else {
            /* last cell created is internal */
         }
      } else {
         /* this is one of the middle segments, we need to know
          * what two cells it is between */
         /* it's either a straight segment within this cell,
          * or it's a bend (then it's tough) */
         if ((next->x == last->x && next->y == last->y) ||
             (next->x == last->x && next->z == last->z) ||
             (next->z == last->z && next->y == last->y)) {
            /* curr is a straight segment */
            (void) find_midpoint(last,curr,p1);
            (void) find_midpoint(curr,next,p2);
            /* fprintf(fp,"# middle cell written - straight\n"); */
            fprintf(fp,"c%d cylinder p%ds%d\n0 0 7 %lf %lf %lf  %lf %lf %lf  %lf\n",
               pnum,pnum,seg_ct,p1->x,p1->y,p1->z,p2->x,p2->y,p2->z,r2);
            if (curve_detail == 1) fprintf(fp,"c%d sphere p%ds%ds\n0 0 4 %lf %lf %lf  %lf\n",
               pnum,pnum,seg_ct,p2->x,p2->y,p2->z,r2);
         } else {
            /* curr is a curved segment */
            (void) find_midpoint(last,curr,p1);
            (void) find_midpoint(curr,next,p3);
            /* fprintf(fp,"# middle cell written - curved\n"); */
            if (curve_detail == 1) {
               /* use a straight segment to denote the curve */
               fprintf(fp,"c%d cylinder p%ds%d\n0 0 7 %lf %lf %lf  %lf %lf %lf  %lf\n",
                  pnum,pnum,seg_ct,p1->x,p1->y,p1->z,p3->x,p3->y,p3->z,r2);
               fprintf(fp,"c%d sphere p%ds%ds\n0 0 4 %lf %lf %lf  %lf\n",
                  pnum,pnum,seg_ct,p3->x,p3->y,p3->z,r2);
            } else {
               /* use multiple segments/spheres for the curve */
               (void) find_center(curr,p2);
               fprintf(fp,"!genworm c%d p%ds%d",pnum,pnum,seg_ct);
               (void) write_curve(fp,p1,p2,p3,curve_detail);
            }
         }
      }

      fprintf(fp,"# segment %d, cell %d %d %d\n",seg_ct,curr->x,curr->y,curr->z);
      fflush(fp);
      seg_ct++;
      last= curr;
      curr= next;
      if (curr) next= next->previous;
   }
}


/*
 * Write a simple PGM description of the top view of the pipes
 */
void write_pgm(int pnum) {

   int i,j,k,max_val,raw_val;
   int xloc, yloc;
   char a[MAXDIM][MAXDIM];
   int height,width;
   char filename[80];
   FILE *fp;

   (void) strcpy(filename,"preview.pgm");
   height = size.y;
   width = size.x;

   fprintf(stderr,"Writing %s...\n\n",filename);

   /*
    * Make the output array
    */

   for (i=height-1; i>-1; i--) {
      for (j=0; j<width; j++) {
         raw_val = -1;
         for (k=0; k<size.z; k++) {
            if (v[i][j][k] == FALSE) raw_val = k;
         }
         a[i][j] = raw_val++;
      }
   }


   /*
    * Write the PGM
    */

   /* find the maximum value of the array */
   /* for (i=0; i<height; i++) {
      for (j=0; j<width; j++) {
         if (a[i][j] > max_val) max_val = a[i][j];
      }
   } */
   max_val = size.z;

   /* open the .pgm file for writing */
   fp = fopen(filename,"w");
   if (fp==NULL) {
      fprintf(stderr,"Could not open output file %s\n",filename);
      exit(0);
   }

   /* write a header */
   fprintf(fp,"P5\n");
   fprintf(fp,"%d %d %d\n",width, height, max_val);

   /* write each data point as a char to the file */
   for (i=0; i<height; i++) {
      for (j=0; j<width; j++) {
         raw_val = a[i][j];
         if (raw_val > 255) raw_val = 255;
         putc((char)raw_val,fp);
      }
   }

   /* close the file before returning to main */
   (void) fclose(fp);
}


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

   /* Usage for hfcolor */
   static char **cpp, *help_message[] =
   {
       "where options include:",
       "  -help          writes this help information",
       "  -out name      name of Radiance output file",
       "  -s x y z       integer size of available space, default is 20 20 5",
       "  -n val         number of pipes, def= 1",
       "  -l min max     length of pipes, minimum and maximum, def= 10 10000",
       "  -f val         probability that a segment will go straight if it can, def= 0.7, range 0-1",
       "  -p             write PGM preview also",
       " ",
       "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);
}


/*
 * Free the memory associated with the last pipe
 */
void free_mem(CELL *end) {

   CELL *prev;

   /* while end is still a cell */
   while (end) {
      prev = end->previous;
      free(end);
      end = prev;
   }
}


/*
 * Return the midpoint of the common face of two cells
 */
void find_midpoint(CELL *start, CELL *end, VEC *point) {

   double csd2 = cs/2.0;
   point->x = csd2 * (1 + (int)start->x + (int)end->x);
   point->y = csd2 * (1 + (int)start->y + (int)end->y);
   point->z = csd2 * (1 + (int)start->z + (int)end->z);
}


/*
 * Return the center of the cell
 */
void find_center(CELL *cell, VEC *point) {

   double csd2 = cs/2.0;
   point->x = csd2 * (1 + 2*(int)cell->x);
   point->y = csd2 * (1 + 2*(int)cell->y);
   point->z = csd2 * (1 + 2*(int)cell->z);
}


/*
 * Draw a curved section of pipe in the given cell
 */
void write_curve(FILE *fp, VEC *st, VEC *ce, VEC *en, int segs) {

   double po2 = 3.14159265358979323846/2.0;

   /* put them together */
   fprintf(fp," '%lf + %lf*(1-cos(t*%lf)) + %lf*sin(t*%lf)'",st->x,en->x-ce->x,po2,ce->x-st->x,po2);
   fprintf(fp," '%lf + %lf*(1-cos(t*%lf)) + %lf*sin(t*%lf)'",st->y,en->y-ce->y,po2,ce->y-st->y,po2);
   fprintf(fp," '%lf + %lf*(1-cos(t*%lf)) + %lf*sin(t*%lf)'",st->z,en->z-ce->z,po2,ce->z-st->z,po2);
   fprintf(fp," '%lf' %d\n",r2,segs);

}