/* File: voronoi.c
 *
 * Constructs 2D Voronoi diagram.  Original code by Bob McMichael.
 * Minor mods by Mike Donahue (mjd), June 2008.
 * Additional mods to interface by mjd, Sept 2008.
 *
 * Usage: See function usage() below.
 *
 * To build:
 *
 *      gcc -O voronoi.c -lm -o voronoi
 *
 * Note: If the grain count is large compared to the relative area,
 *       then some nucleation points can coincide and the actual
 *       number of grains may be slightly less than requested.
 *
 */

#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>


/* This program finds nearest nucleus match to each pixel by grouping
 * nuclei into locales, and searching only nearby locales.  Uncomment
 * the following line to use a slower search through the entire nuclei
 * list.  This is intended primarily for debugging and testing.
 */
#ifndef SEARCH_ALL
# define SEARCH_ALL 0
#endif

#define VORONOI_VERSION "1.1 20101101"

/* Set PGM_OUTPUT macro to 1 to get PGM grayscale output
 * instead of PPM P6 output.
 */
#ifndef PGM_OUTPUT
# define PGM_OUTPUT 0
#endif

/* On non-unix machines, stdout might not be opened in binary mode.  */
/* Use this macro to force the issue.                                */
#ifndef FORCE_BINARY_STDOUT
# if defined(_WIN32) || defined(WIN32) || defined(_WINNT) || defined(WINNT)
#  include <fcntl.h>
#  include <io.h>
#  define FORCE_BINARY_STDOUT _setmode(_fileno(stdout),_O_BINARY)
# else
#  define FORCE_BINARY_STDOUT
# endif
#endif

#define DEF_Nx 512
#define DEF_Ny 128
#define DEF_grains 2000
#if !PGM_OUTPUT
# define DEF_grain_bdry_red   0xFF
# define DEF_grain_bdry_green 0xFF
# define DEF_grain_bdry_blue  0xFF
#else
# define DEF_grain_bdry_red   0xEF
# define DEF_grain_bdry_green 0xEF
# define DEF_grain_bdry_blue  0xEF
#endif

#define DEF_basename "rect"

#ifndef SQRT1_2
# define SQRT1_2 0.70710678118654757
#endif

void usage(void)
{
  fprintf(stderr,"Usage: voronoi [-x Nx] [-y Ny] [-p edges] [-g ngrains] \\\n"
     "         [-gbdrywidth bdrywidth] [-gbdrycolor bdrycolor] [-s seed] \\\n"
#if PGM_OUTPUT
          "         [-o name] [-version]"
#else
          "         [-oldoutput] [-o name] [-version]"
#endif
          " [-nored] [-nogreen] [-noblue]\n"
          " where  Nx = number of pixels in x direction (default = %d)\n"
          "        Ny = number of pixels in y direction (default = %d\n"
          "     edges = periodic edges; one of x, y, xy, or none"
          "  (default = none)\n"
          "   ngrains = number of grains (default = %d)\n"
          " bdrywidth = width of grain boundaries, in pixels"
          " (default = 0.0)\n"
          " bdrycolor = color of grain boundaries, as rgb or rrggbb,\n"
          "             where r, g, b are hex digits"
          " (default = %02x%02x%02x)\n"
          "      seed = random number generator seed (default from clock)\n"
          "      name = basename for output images"
          "  (default = %s; use \"-\" for stdout)\n",
          DEF_Nx,DEF_Ny,DEF_grains,
          DEF_grain_bdry_red,DEF_grain_bdry_green,DEF_grain_bdry_blue,
          DEF_basename);
  fprintf(stderr,"The -nored, -noblue, and -nogreen options restrict"
          " the output color space\n   in the indicated manner.\n");
#if PGM_OUTPUT
  fprintf(stderr,"The output is a grain map in PGM P5 (binary) image"
          " format, \"name\".pgm,\n"
          "   where \"name\" is as specified by the -o option.\n");
#else
  fprintf(stderr,"The output is a grain map in PPM P6 (binary) image"
          " format, \"name\".ppm,\n"
          "   where \"name\" is as specified by the -o option.\n");
  fprintf(stderr,"If the -oldoutput flag is present, then a second"
          " reduced-contrast output\n"
          "   file is written to name-mif.ppm.\n");
#endif
  fprintf(stderr,"Return code is 0 on success, > 0 on error.\n");
  exit(1);
}

/*** IMAGE OUTPUT ROUTINES ********************************************/
/*
 * Support code for image output.  Original code by Bob McMichael.
 * DumpToPPM function added by Mike Donahue (June 2008).
 *
 */

#define ALLSCALE 1

#if PGM_OUTPUT

void DumpToPGM(int *arr, int xsize, int ysize, char *outfile,
               int bdryval) {
  int i, j;
  int ij;
  int maxval, minval;
  unsigned char scaledval;
  FILE *out;
  size_t size;

  if(xsize < 1 || ysize < 1) return;

  /* determine the max and min values for scaling */

  maxval = 0;
  minval = INT_MAX;

  for(i = 0; i < xsize; i++){
    for(j = 0; j < ysize; j++){
      ij= i*ysize+j;
      if(arr[ij]>=0) {
        // arr[ij]==-1 is special case indicating grain boundary
        if ( arr[ij] > maxval ) maxval = arr[ij];
        if ( arr[ij] < minval ) minval = arr[ij];
      }
    }
  }

  if(strcmp("-",outfile)==0) {
    FORCE_BINARY_STDOUT; /* Make sure stdout is in binary mode */
    out = stdout;
  } else {
    out = fopen(outfile, "wb");
  }
  fprintf(out, "P5\n%d %d\n255\n", xsize, ysize);

  size = sizeof(unsigned char);
  /* scale the data and write it out */
  /* top row (large y) first */
  if(maxval > 1 || minval < -1 || ALLSCALE){
    for(j = ysize -1; j > -1; j--){
      for(i = 0; i < xsize; i++){
        int tval;
        ij= i*ysize+j;
        tval = arr[ij];
        if(tval == -1) {
          scaledval = bdryval;
        } else {
          if(maxval == minval) {
            scaledval = 127;
          } else {
            scaledval = (unsigned char)(255.*(tval-minval)
                                        /(double)(maxval-minval));
          }
        }
        fwrite(&scaledval, size, 1, out);
      }
    }
  } else {
    for(j = ysize -1; j > -1; j--){
      for(i = 0; i < xsize; i++){
        int tval;
        ij= i*ysize+j;
        tval = arr[ij];
        if(tval == -1) {
          scaledval = bdryval;
        } else {
          scaledval = (unsigned char)(127*(tval+1));
        }
        fwrite(&scaledval, size, 1, out);
      }
    }
  }
  if(strcmp("-",outfile)!=0) {
    fclose(out);
  }
}

#else /* PGM_OUTPUT */

void DumpToPPM(int *arr, int xsize, int ysize,
               int maxsize, char *outfile,
               int flag_nored,int flag_nogreen, int flag_noblue,
               int bdry_red, int bdry_green, int bdry_blue) {
  int i, j;
  FILE *out;

  unsigned int maxcomp;  /* Maximum component value */
  unsigned int compscale;  /* Component scaling (equivalently, stepping) */

  /* Determine component scaling */
  if(maxsize==0) {
    /* Special case: use colorstep = 1 scaling */
    maxcomp = 256;
  } else {
    int component_count = 3;
    if(flag_nored) --component_count;
    if(flag_nogreen) --component_count;
    if(flag_noblue) --component_count;
    if(component_count<1) {
      maxcomp = maxsize;
    } else {
      maxcomp = (unsigned int)ceil(pow((double)maxsize,
                                       1./(double)component_count));
    }
  }
  compscale = 256/maxcomp;  /* Integer division (with truncation) */
  if(compscale<1) compscale=1; /* Round-off protection */

  /* Dump data */
  if(strcmp("-",outfile)==0) {
    FORCE_BINARY_STDOUT; /* Make sure stdout is in binary mode */
    out = stdout;
  } else {
    out = fopen(outfile, "wb");
  }

  fprintf(out, "P6\n%d %d\n255\n", xsize, ysize);

  for(j = ysize -1; j > -1; j--){
    for(i = 0; i < xsize; i++){
      int ij= i*ysize+j;
      unsigned char pix[3];
      int tval = arr[ij];
      if(tval == -1) {
        /* Special code for boundary */
        pix[0] = bdry_red;
        pix[1] = bdry_green;
        pix[2] = bdry_blue;
      } else {
        unsigned int val = (unsigned int) tval;
        if(flag_nored && flag_nogreen && flag_noblue) {
          int tmp = val * compscale;
          if(tmp>255) tmp = 255;
          pix[2] = pix[1] = pix[0] = tmp;
        } else {
          if(flag_noblue) {
            pix[2] = 0;
          } else {
            int tmp = (val % maxcomp)*compscale;
            if(tmp>255) tmp = 255;
            pix[2] = tmp;
            val /= maxcomp;
          }
          if(flag_nogreen) {
            pix[1] = 0;
          } else {
            int tmp = (val % maxcomp)*compscale;
            if(tmp>255) tmp = 255;
            pix[1] = tmp;
            val /= maxcomp;
          }
          if(flag_nored) {
            pix[0] = 0;
          } else {
            int tmp = val * compscale;
            if(tmp>255) tmp = 255;
            pix[0] = tmp;
          }
        }
      }
      fwrite(pix,sizeof(unsigned char), 3, out);
    }
  }
  if(strcmp("-",outfile)!=0) {
    fclose(out);
  }
}
#endif /* PGM_OUTPUT */

#undef ALLSCALE

/**********************************************************************/

struct nucleus {
  int x;
  int y;
  int val;
  int li;
};

struct locale {
  int nucoff;  /* Index to first elt of nuc in this locale. */
  int ncnt;    /* Number of nuc in this locale.             */
               /* NB: The elts of nuc are sorted by locale. */
  int min_i;   /* Required search neighborhood */
  int max_i;
  int min_j;
  int max_j;
};

/* for sorting nuclei by locale */
int loccomp(const void *thing1,const void *thing2) {
  int li1,li2,x1,x2,y1,y2;
  li1 = ((const struct nucleus *)thing1)->li;
  li2 = ((const struct nucleus *)thing2)->li;
  if( li1 > li2 ) {
    return 1 ;
  } else if(li1<li2) {
    return -1 ;
  }
  /* li1 == li2 */
  x1 = ((const struct nucleus *)thing1)->x;
  x2 = ((const struct nucleus *)thing2)->x;
  if( x1 > x2 ) {
    return 1 ;
  } else if(x1<x2) {
    return -1 ;
  }
  /* x1 == x2 */
  y1 = ((const struct nucleus *)thing1)->y;
  y2 = ((const struct nucleus *)thing2)->y;
  if( y1 > y2 ) {
    return 1 ;
  } else if(y1<y2) {
    return -1 ;
  }
  return 0; /* Identical positions */
}

/* Global variables */
int Nx=DEF_Nx;
int Ny=DEF_Ny;
int periodic_x = 0;
int periodic_y = 0;
int grain_count=DEF_grains;
int nlx = 0; /* x-dimension locale count */
int nly = 0; /* y-dimension locale count */
struct nucleus *nuc;
struct locale *loc;

double grain_bdry_width = 0.0;
unsigned int grain_bdry_color[3] = {DEF_grain_bdry_red,
                                    DEF_grain_bdry_green,
                                    DEF_grain_bdry_blue};
int search_all_closest(int x,int y,int* nuc_index,double* nuc_dist)
{
  int index,min_index,xdiff,ydiff,testdist;
  double distsq,min_distsq;

  if(grain_count<1) return -1;

  xdiff = abs(x - nuc[0].x);
  if(periodic_x && 2*xdiff > Nx) xdiff = Nx - xdiff;
  ydiff = abs(y - nuc[0].y);
  if(periodic_y && 2*ydiff > Ny) ydiff = Ny - ydiff;
  testdist = xdiff + ydiff;
  min_distsq = (double)xdiff*(double)xdiff
             + (double)ydiff*(double)ydiff;
  min_index = 0;

  for(index=1;index<grain_count;++index) {
    xdiff = abs(x - nuc[index].x);
    if(periodic_x && 2*xdiff > Nx) xdiff = Nx - xdiff;
    if(xdiff >= testdist) continue;

    ydiff = abs(y - nuc[index].y);
    if(periodic_y && 2*ydiff > Ny) ydiff = Ny - ydiff;
    if(ydiff >= testdist) continue;

    distsq = (double)xdiff*(double)xdiff
           + (double)ydiff*(double)ydiff;

    if(distsq<min_distsq) {
      min_distsq = distsq;
      testdist = xdiff + ydiff;
      min_index = index;
    }
  }
#if 0
  fprintf(stderr,"x=%d, y=%d; closest: %d,%d (val=%d)\n",
         x,y,nuc[min_index].x,nuc[min_index].y,nuc[min_index].val);
#endif
  if(nuc_index != NULL) *nuc_index = min_index;
  if(nuc_dist  != NULL) *nuc_dist = sqrt(min_distsq);
  return nuc[min_index].val;
}

int search_all_closest2(int x,int y,int* nuc_index,double* nuc_dist)
{
  int index,min_index,xdiff,ydiff,testdist;
  double distsq,min_distsq;

  int min2_index,testdist2,min3_index,testdist3;
  double min2_distsq,min3_distsq;

  int inside_grain_boundary;
  double gbwsq;

  if(grain_count<1) return -1;

  gbwsq = grain_bdry_width*grain_bdry_width;

  min3_index = min2_index = -1;
  testdist3 = testdist2 = Nx + Ny;
  min3_distsq = min2_distsq = testdist2*testdist2;

  xdiff = abs(x - nuc[0].x);
  if(periodic_x && 2*xdiff > Nx) xdiff = Nx - xdiff;
  ydiff = abs(y - nuc[0].y);
  if(periodic_y && 2*ydiff > Ny) ydiff = Ny - ydiff;
  testdist = xdiff + ydiff;
  min_distsq = (double)xdiff*(double)xdiff
             + (double)ydiff*(double)ydiff;
  min_index = 0;

  for(index=1;index<grain_count;++index) {
    xdiff = abs(x - nuc[index].x);
    if(periodic_x && 2*xdiff > Nx) xdiff = Nx - xdiff;
    if(xdiff >= testdist3) continue;

    ydiff = abs(y - nuc[index].y);
    if(periodic_y && 2*ydiff > Ny) ydiff = Ny - ydiff;
    if(ydiff >= testdist3) continue;

    distsq = (double)xdiff*(double)xdiff
           + (double)ydiff*(double)ydiff;

    if(distsq<min_distsq) {
      min3_distsq = min2_distsq;
      min3_index  = min2_index;
      testdist3   = testdist2;
      min2_distsq = min_distsq;
      min2_index  = min_index;
      testdist2   = testdist;
      min_distsq = distsq;
      min_index = index;
      testdist = xdiff + ydiff;
    } else if(distsq<min2_distsq) {
      min3_distsq = min2_distsq;
      min3_index  = min2_index;
      testdist3   = testdist2;
      min2_distsq = distsq;
      min2_index = index;
      testdist2 = xdiff + ydiff;
    } else if(distsq<min3_distsq) {
      min3_distsq = distsq;
      min3_index = index;
      testdist3 = xdiff + ydiff;
    }
  }
#if 0
  fprintf(stderr,"x=%d, y=%d; closest: %d,%d (val=%d)\n",
         x,y,nuc[min_index].x,nuc[min_index].y,nuc[min_index].val);
#endif

  inside_grain_boundary = 0;
  if(gbwsq>0.0 && min2_index>=0) {
    /* Do draw grain boundaries */
    /* Use vector arithmetic to test grain boundary extents */
    double ax,ay,bx,by,cx,cy;
    double csq,test;
    ax = nuc[min_index].x - x;
    ay = nuc[min_index].y - y;
    bx = nuc[min2_index].x - x;
    by = nuc[min2_index].y - y;

    cx = nuc[min2_index].x - nuc[min_index].x;
    cy = nuc[min2_index].y - nuc[min_index].y;
    csq = cx*cx + cy*cy;

    test = (ax + bx)*cx + (ay + by)*cy;

    if(test*test <= gbwsq*csq) {
      inside_grain_boundary = 1;
    } else if(min3_index>=0) {
      /* Pixel probably not inside grain boundary, but check
       * next bisector too.
       */
      bx = nuc[min3_index].x - x;
      by = nuc[min3_index].y - y;
      cx = nuc[min3_index].x - nuc[min_index].x;
      cy = nuc[min3_index].y - nuc[min_index].y;
      csq = cx*cx + cy*cy;
      test = (ax + bx)*cx + (ay + by)*cy;
      if(test*test <= gbwsq*csq) {
        inside_grain_boundary = 1;
      }
      /* Otherwise, pixel even more likely not inside grain boundary.
       * Go with it.
       */
    }
  }

  if(inside_grain_boundary) {
    if(nuc_index != NULL) *nuc_index = -1; /* Code denoting grain bdry */
    if(nuc_dist  != NULL) *nuc_dist  = -1; /* Ditto */
    return -1;
  }

  if(nuc_index != NULL) *nuc_index = min_index;
  if(nuc_dist  != NULL) *nuc_dist = sqrt(min_distsq);
  return nuc[min_index].val;
}

void find_closest(int index,int x,int y,int* nuc_index,double* nuc_dist)
{
  int i,j,k,min_i,max_i,min_j,max_j,min_index,locoff;
  int xdiff,ydiff,testdist;
  double distsq, min_distsq;
  struct nucleus *nucptr;
  min_i = loc[index].min_i;
  max_i = loc[index].max_i;
  min_j = loc[index].min_j;
  max_j = loc[index].max_j;

  testdist = Nx+Ny;
  min_distsq = (double)(testdist); /* Upper bound */
  min_distsq *= min_distsq;
  min_index = -1;

  for(i=min_i;i<=max_i;++i) {
    int iw = i;
    if(periodic_x) {
      if(iw<0)         iw += nlx;
      else if(iw>=nlx) iw -= nlx;
    }
    iw *= nly;
    for(j=min_j;j<=max_j;++j) {
      int jw = j;
      if(periodic_y) {
        if(jw<0)         jw += nly;
        else if(jw>=nly) jw -= nly;
      }
      locoff = iw + jw;
      nucptr = nuc + loc[locoff].nucoff;
      for(k=0;k<loc[locoff].ncnt;++k) {
        xdiff = abs(nucptr[k].x - x);
        if(periodic_x && 2*xdiff > Nx) xdiff = Nx - xdiff;
        if(xdiff>testdist) continue;
        ydiff = abs(nucptr[k].y - y);
        if(periodic_y && 2*ydiff > Ny) ydiff = Ny - ydiff;
        if(ydiff>testdist) continue;
        distsq = (double)xdiff*(double)xdiff
               + (double)ydiff*(double)ydiff;
        if(distsq<min_distsq
           || (distsq == min_distsq && loc[locoff].nucoff + k < min_index) ) {
          min_distsq = distsq;
          testdist = xdiff + ydiff;
          min_index = loc[locoff].nucoff + k;
        }
      }
    }
  }

  if(min_index>=0) {
    /* Match found */
    if(nuc_index != NULL) *nuc_index = min_index;
    if(nuc_dist  != NULL) *nuc_dist  = sqrt(min_distsq);
  } else {
    /* No points in specified neighborhood.  Fall back to full search. */
    search_all_closest(x,y,nuc_index,nuc_dist);
  }

}

void find_closest2(int index,int x,int y,int* nuc_index,double* nuc_dist)
{ /* Modified version of find_closest that checks two closest nuclei
   * so grain boundaries can be identified.
   */
  int i,j,k,min_i,max_i,min_j,max_j,min_index,min2_index,min3_index,locoff;
  int xdiff,ydiff,testdist,testdist2,testdist3;
  double distsq,min_distsq,min2_distsq,min3_distsq,gbwsq;
  min_i = loc[index].min_i;
  max_i = loc[index].max_i;
  min_j = loc[index].min_j;
  max_j = loc[index].max_j;


  testdist = Nx+Ny;
  min_distsq = (double)(testdist); /* Upper bound */
  min_distsq *= min_distsq;
  min_index = -1;

  min3_distsq = min2_distsq = min_distsq;
  min3_index  = min2_index  = min_index;
  testdist3   = testdist2   = testdist;

  gbwsq = grain_bdry_width*grain_bdry_width;

  for(i=min_i;i<=max_i;++i) {
    int iw = i;
    if(periodic_x) {
      if(iw<0)         iw += nlx;
      else if(iw>=nlx) iw -= nlx;
    }
    iw *= nly;
    for(j=min_j;j<=max_j;++j) {
      struct nucleus *nucptr;
      int jw = j;
      if(periodic_y) {
        if(jw<0)         jw += nly;
        else if(jw>=nly) jw -= nly;
      }
      locoff = iw + jw;
      nucptr = nuc + loc[locoff].nucoff;
      for(k=0;k<loc[locoff].ncnt;++k) {
        xdiff = abs(nucptr[k].x - x);
        if(periodic_x && 2*xdiff > Nx) xdiff = Nx - xdiff;
        if(xdiff>testdist3) continue;
        ydiff = abs(nucptr[k].y - y);
        if(periodic_y && 2*ydiff > Ny) ydiff = Ny - ydiff;
        if(ydiff>testdist3) continue;
        distsq = (double)xdiff*(double)xdiff
               + (double)ydiff*(double)ydiff;

        if(distsq<min_distsq
           || (distsq == min_distsq && loc[locoff].nucoff + k < min_index) ) {
          min3_distsq = min2_distsq;
          min3_index  = min2_index;
          testdist3   = testdist2;
          min2_distsq = min_distsq;
          min2_index  = min_index;
          testdist2   = testdist;
          min_distsq = distsq;
          min_index = loc[locoff].nucoff + k;
          testdist = xdiff + ydiff;
        } else if(distsq<min2_distsq
           || (distsq == min2_distsq && loc[locoff].nucoff + k < min2_index) ) {
          min3_distsq = min2_distsq;
          min3_index  = min2_index;
          testdist3   = testdist2;
          min2_distsq = distsq;
          min2_index = loc[locoff].nucoff + k;
          testdist2 = xdiff + ydiff;
        } else if(distsq<min3_distsq
           || (distsq == min3_distsq && loc[locoff].nucoff + k < min3_index) ) {
          min3_distsq = distsq;
          min3_index = loc[locoff].nucoff + k;
          testdist3 = xdiff + ydiff;
        }
      }
    }
  }

  if(gbwsq==0.0) { /* gbwsq := grain boundary width squared */
    /* Don't draw grain boundaries */
    if(min_index>=0) {
      /* Match found */
      if(nuc_index != NULL) *nuc_index = min_index;
      if(nuc_dist  != NULL) *nuc_dist  = sqrt(min_distsq);
    } else {
      /* No points in specified neighborhood.  Fall back to full search. */
      search_all_closest(x,y,nuc_index,nuc_dist);
    }
  } else {
    /* Do draw grain boundaries */
    if(min_index>=0 && min2_index>=0) {
      /* Use vector arithmetic to test grain boundary extents */
      double ax,ay,bx,by,cx,cy;
      double csq,test;
      int inside_grain_boundary = 0;
      ay = nuc[min_index].y - y;
      bx = nuc[min2_index].x - x;
      by = nuc[min2_index].y - y;
        
      ax = nuc[min_index].x - x;
      ay = nuc[min_index].y - y;
      bx = nuc[min2_index].x - x;
      by = nuc[min2_index].y - y;

      cx = nuc[min2_index].x - nuc[min_index].x;
      cy = nuc[min2_index].y - nuc[min_index].y;
      csq = cx*cx + cy*cy;

      test = (ax + bx)*cx + (ay + by)*cy;

      if(test*test <= gbwsq*csq) {
        inside_grain_boundary = 1;
      } else if(min3_index>=0) {
        /* Pixel probably not inside grain boundary, but check
         * next bisector too.
         */
        bx = nuc[min3_index].x - x;
        by = nuc[min3_index].y - y;
        cx = nuc[min3_index].x - nuc[min_index].x;
        cy = nuc[min3_index].y - nuc[min_index].y;
        csq = cx*cx + cy*cy;
        test = (ax + bx)*cx + (ay + by)*cy;
        if(test*test <= gbwsq*csq) {
          inside_grain_boundary = 1;
        }
        /* Otherwise, pixel even more likely not inside grain boundary.
         * Go with it.
         */
      }
      if(inside_grain_boundary) {
        if(nuc_index != NULL) *nuc_index = -1; /* Code denoting grain bdry */
        if(nuc_dist  != NULL) *nuc_dist  = -1; /* Ditto */
      } else {
        if(nuc_index != NULL) *nuc_index = min_index;
        if(nuc_dist  != NULL) *nuc_dist  = sqrt(min_distsq);
      }
    } else {
      /* Didn't find two points in specified neighborhood.
       * Fall back to full search.
       */
      search_all_closest2(x,y,nuc_index,nuc_dist);
    }
  }
}

int main(int argc,const char **argv)
{
  double max;
  int block_size,last_block_xsize,last_block_ysize;
  int i, j, index, ncnt;
  int offset;
  int count;
  int *array, **pa;
  int seed=1, seed_set=0;
  int flag_nored=0, flag_nogreen=0, flag_noblue=0, flag_version=0;
  int dup_count=0;
  int dup_removal_attempts = 100;
  int return_code = 0;
  int flag_oldoutput = 0;
  const char* basename = DEF_basename;
  size_t namebuf_size;
  char* namebuf;
  char* namebuf_mif;

  i=1;

  while(i<argc) {
    if(strcmp(argv[i],"-x")==0) {
      if(i+1>=argc) usage();
      Nx = atoi(argv[i+1]);
      if(Nx<1) usage();
      for(j=i+2;j<argc;++j) argv[j-2] = argv[j];
      argc -= 2;
    } else if(strcmp(argv[i],"-y")==0) {
      if(i+1>=argc) usage();
      Ny = atoi(argv[i+1]);
      if(Ny<1) usage();
      for(j=i+2;j<argc;++j) argv[j-2] = argv[j];
      argc -= 2;
    } else if(strcmp(argv[i],"-p")==0
              || strcmp(argv[i],"-periodic")==0
              || strcmp(argv[i],"--periodic")==0) {
      if(i+1>=argc) usage();
      if(strcmp(argv[i+1],"x")==0) {
        periodic_x = 1;
        periodic_y = 0;
      } else if(strcmp(argv[i+1],"y")==0) {
        periodic_x = 0;
        periodic_y = 1;
      } else if(strcmp(argv[i+1],"xy")==0 || strcmp(argv[i+1],"yx")==0) {
        periodic_x = 1;
        periodic_y = 1;
      } else if(strcmp(argv[i+1],"none")==0) {
        periodic_x = 0;
        periodic_y = 0;
      } else {
        usage();
      }
      for(j=i+2;j<argc;++j) argv[j-2] = argv[j];
      argc -= 2;
    } else if(strcmp(argv[i],"-g")==0) {
      if(i+1>=argc) usage();
      grain_count = atoi(argv[i+1]);
      if(grain_count<1) usage();
      for(j=i+2;j<argc;++j) argv[j-2] = argv[j];
      argc -= 2;

    } else if(strcmp(argv[i],"-gbdrywidth")==0) {
      if(i+1>=argc) usage();
      grain_bdry_width = atof(argv[i+1]);
      if(grain_bdry_width<0.0) usage();
      for(j=i+2;j<argc;++j) argv[j-2] = argv[j];
      argc -= 2;

    } else if(strcmp(argv[i],"-gbdrycolor")==0) {
      unsigned long gbc;
      const char *cptr;
      if(i+1>=argc) usage();
      cptr = argv[i+1];
      if(strlen(cptr)!=3 && strlen(cptr)!=6) {
        usage();
      }
      gbc = strtoul(cptr,NULL,16);
      if(strlen(cptr)==3) {
        /* #rgb    */
        grain_bdry_color[0] = 16*((gbc & 0xF00) >> 8);
        grain_bdry_color[1] = 16*((gbc & 0x0F0) >> 4);
        grain_bdry_color[2] = 16*((gbc & 0x00F));
      } else {
        /* #rrggbb */
        grain_bdry_color[0] = (gbc & 0xFF0000) >> 16;
        grain_bdry_color[1] = (gbc & 0x00FF00) >>  8;
        grain_bdry_color[2] = (gbc & 0x0000FF);
      }
      if(grain_bdry_color[0]>255
         || grain_bdry_color[1]>255 || grain_bdry_color[2]>255) {
        usage();
      }
      for(j=i+2;j<argc;++j) argv[j-2] = argv[j];
      argc -= 2;
    } else if(strcmp(argv[i],"-s")==0) {
      if(i+1>=argc) usage();
      seed = atoi(argv[i+1]);
      seed_set = 1;
      for(j=i+2;j<argc;++j) argv[j-2] = argv[j];
      argc -= 2;
    } else if(strcmp(argv[i],"-oldoutput")==0) {
      flag_oldoutput=1;
      for(j=i+1;j<argc;++j) argv[j-1] = argv[j];
      argc -= 1;
    } else if(strcmp(argv[i],"-o")==0) {
      if(i+1>=argc) usage();
      basename = argv[i+1];
      if(strlen(basename)<1) usage();
      for(j=i+2;j<argc;++j) argv[j-2] = argv[j];
      argc -= 2;
    } else if(strcmp(argv[i],"-nored")==0) {
      flag_nored=1;
      for(j=i+1;j<argc;++j) argv[j-1] = argv[j];
      argc -= 1;
    } else if(strcmp(argv[i],"-nogreen")==0) {
      flag_nogreen=1;
      for(j=i+1;j<argc;++j) argv[j-1] = argv[j];
      argc -= 1;
    } else if(strcmp(argv[i],"-noblue")==0) {
      flag_noblue=1;
      for(j=i+1;j<argc;++j) argv[j-1] = argv[j];
      argc -= 1;
    } else if(strcmp(argv[i],"--version")==0
              || strcmp(argv[i],"-version")==0) {
      flag_version=1;
      for(j=i+1;j<argc;++j) argv[j-1] = argv[j];
      argc -= 1;
    } else {
      ++i;
    }
  }
  if(argc != 1) usage();
  if(flag_version) {
    fprintf(stdout,"voronoi version %s\n",VORONOI_VERSION);
    exit(0);
  }
  if(grain_count > Nx * Ny) {
    fprintf(stderr,
            "\nERROR: More grains (%d) requested than pixels (%d).\n\n",
            grain_count, Nx * Ny);
    usage();
  }
  if(flag_nored && flag_nogreen && flag_noblue) {
    /* Gray-scale output request */
    if(grain_count>256) {
      fprintf(stderr,"\nERROR: More grains requested than"
              " available gray shades.\n\n");
      usage();
    }
  } else {
    int color_count = 1;
    if(!flag_nored)   color_count *= 256;
    if(!flag_nogreen) color_count *= 256;
    if(!flag_noblue)  color_count *= 256;
    if(grain_count>color_count) {
      fprintf(stderr,"\nERROR: More grains requested (%d) than"
              " available colors (%d).\n\n",grain_count,color_count);
      usage();
    }
  }

  if(seed_set) {
    srand(seed);
  } else {
    srand(time(NULL));
  }

  max=RAND_MAX;

  fprintf(stderr,"Nx=%d, Ny=%d, grain_count=%d\n", Nx, Ny, grain_count);

  /* allocate memory */
  array = (int *)malloc(Nx*Ny*sizeof(int));
  pa = (int **)malloc(Nx*sizeof(int *));
  for(i=0;i<Nx;i++){
    pa[i]=array+i*Ny;
  };
  for(i=0;i<Nx*Ny;i++) array[i] = -1;  /* initialize */

  block_size = (int)ceil(sqrt((10.*Nx*Ny)/grain_count));
  if (block_size > Nx) block_size = Nx;
  if (block_size > Ny) block_size = Ny;
  nlx = (Nx+block_size-1)/block_size;
  nly = (Ny+block_size-1)/block_size;
  loc = (struct locale *)malloc(nlx*nly*sizeof(struct locale));
  nuc = (struct nucleus *)malloc(grain_count*sizeof(struct nucleus));

  last_block_xsize = Nx - (nlx-1)*block_size;
  last_block_ysize = Ny - (nly-1)*block_size;

  /* load array of nuclei */
  for (i=0 ; i<grain_count ; ++i) {
    nuc[i].x   = (int)floor(Nx*((double)rand()/((double)max+1)));
    if(nuc[i].x > Nx-1) nuc[i].x = Nx - 1;  /* Safety */
    nuc[i].y   = (int)floor(Ny*((double)rand()/((double)max+1)));
    if(nuc[i].y > Ny-1) nuc[i].y = Ny - 1;  /* Safety */
    nuc[i].val =  i;
    nuc[i].li  = nly*(nuc[i].x/block_size) + nuc[i].y/block_size;
  }

  /* sort the nuclei by locale  */
  qsort(nuc, grain_count, sizeof(struct nucleus), loccomp);

  do {
    int x0,y0,x1,y1;

    /* Check for duplicate positions.  This depends upon proper
     * sorting in the previous step.
     */
    dup_count = 0;
    x0 = nuc[0].x;
    y0 = nuc[0].y;
    for(i=1;i<grain_count;++i) {
      x1 = nuc[i].x;
      y1 = nuc[i].y;
      if(x1 == x0 && y1 == y0) {
        /* Duplicate; try again. */
        nuc[i].x   = (int)floor(Nx*((double)rand()/((double)max+1)));
        if(nuc[i].x > Nx-1) nuc[i].x = Nx - 1;  /* Safety */
        nuc[i].y   = (int)floor(Ny*((double)rand()/((double)max+1)));
        if(nuc[i].y > Ny-1) nuc[i].y = Ny - 1;  /* Safety */
        nuc[i].li  = nly*(nuc[i].x/block_size) + nuc[i].y/block_size;
        ++dup_count;
      }
      x0 = x1;
      y0 = y1;
    }

    /* sort the nuclei by locale.  At this point the list is already
     * nearly sorted, so one might want to use something other than
     * qsort.
     */
    qsort(nuc, grain_count, sizeof(struct nucleus), loccomp);

    if(dup_count>0 && dup_removal_attempts<1) {
      /* Get fully accurate dup count */
      dup_count = 0;
      for(i=1;i<grain_count;++i) {
        x1 = nuc[i].x;
        y1 = nuc[i].y;
        if(x1 == x0 && y1 == y0) {
          ++dup_count;
        }
        x0 = x1;
        y0 = y1;
      }
      if(dup_count>0) {
        fflush(stdout);
        fprintf(stderr,"WARNING: Only %d distinct grains created.\n",
                grain_count - dup_count);
        fflush(stderr);
        return_code |= 0x02;
      }
    }

  } while(dup_count>0 && dup_removal_attempts-- > 0);

  /* Set up locale nucleus pointers to point into nuc array */
  for(index = 0;index<nlx*nly;++index) {
    loc[index].nucoff = -1;
    loc[index].ncnt = 0;
  }
  for(i=0;i<grain_count;++i) { /* Fill nuc info */
    int index = nuc[i].li;
    if(loc[index].nucoff == -1) {
      loc[index].nucoff = i;
    }
    ++loc[index].ncnt;
  }
  index = 0;
  for(i=0;i<nlx;++i) {
    int cix;
    int block_width = block_size;
    if(i==nlx-1) block_width = last_block_xsize;
    cix = i*block_size + block_width/2; /* x-center of block */
    for(j=0;j<nly;++j) {
      int nuc_index,itmp,jtmp;
      double nuc_dist; /* In pixels */
      int xsearch_dist,ysearch_dist; /* In pixels */
      int cjy;
      int block_height = block_size;
      if(j==nly-1) block_height = last_block_ysize;
      cjy = j*block_size + block_height/2; /* y-center of block */

      loc[index].min_i = i;
      loc[index].max_i = i;
      loc[index].min_j = j;
      loc[index].max_j = j;
      find_closest(index,cix,cjy,&nuc_index,&nuc_dist);
      xsearch_dist = ysearch_dist = (int)ceil(nuc_dist + block_size*SQRT1_2);
      if(periodic_x) {
        if(2*xsearch_dist>Nx) xsearch_dist = Nx/2;
      } else {
        if(xsearch_dist>=Nx) xsearch_dist = Nx-1;
      }
      if(periodic_y) {
        if(2*ysearch_dist>Ny) ysearch_dist = Ny/2;
      } else {
        if(ysearch_dist>=Ny) ysearch_dist = Ny-1;
      }

      itmp = cix - (block_width+1)/2 - xsearch_dist;
      if(itmp<0) {
        if(!periodic_x) {
          itmp = 0;
        } else {
          itmp = -itmp - last_block_xsize;
          if(itmp<0) itmp = -1;
          else       itmp = -1*(1 + (itmp+block_size-1)/block_size);
        }
      } else {
        itmp /= block_size;
      }
      loc[index].min_i = itmp;

      itmp = cix + (block_width+1)/2 + xsearch_dist;
      if(itmp>=Nx) {
        if(!periodic_x) {
          itmp = nlx-1;
        } else {
          itmp -= Nx;
          itmp = nlx + itmp/block_size;
        }
      } else {
        itmp /= block_size;
      }
      loc[index].max_i = itmp;

      jtmp = cjy - (block_height+1)/2 - ysearch_dist;
      if(jtmp<0) {
        if(!periodic_y) {
          jtmp = 0;
        } else {
          jtmp = -jtmp - last_block_ysize;
          if(jtmp<0) jtmp = -1;
          else       jtmp = -1*(1 + (jtmp+block_size-1)/block_size);
        }
      } else {
        jtmp /= block_size;
      }
      loc[index].min_j = jtmp;

      jtmp = cjy + (block_height+1)/2 + ysearch_dist;
      if(jtmp>=Ny) {
        if(!periodic_y) {
          jtmp = nly-1;
        } else {
          jtmp -= Ny;
          jtmp = nly + jtmp/block_size;
        }
      } else {
        jtmp /= block_size;
      }
      loc[index].max_j = jtmp;
      ++index;
    }
  }

  /* now, for every point on the map, we have a locale & neighbor locales,
     each with a list of the nuclei it contains */

  for (i = 0; i < Nx ; ++i) {
    int block_i = i/block_size;
    block_i *= nly;
    for (j = 0; j < Ny; ++j) {
      double dist;
#if !SEARCH_ALL
      int block_j = j/block_size;
      if(grain_bdry_width>0.0) {
        find_closest2(block_i + block_j,i,j,&index,&dist);
      } else {
        find_closest(block_i + block_j,i,j,&index,&dist);
      }
      if(index == -1) {
        pa[i][j] = -1; /* Code for grain boundary */
      } else {
        pa[i][j] = nuc[index].val;
      }
#else /* SEARCH_ALL */
      if(grain_bdry_width>0.0) {
        pa[i][j] = search_all_closest2(i,j,&index,&dist);
      } else {
        pa[i][j] = search_all_closest(i,j,&index,&dist);
      }
#endif
#if 0
      fprintf(stderr,"pa[%3d][%3d] = %3d  at %6d (%3d,%3d) dist = %14.8f\n",
             i,j,pa[i][j],index,nuc[index].x,nuc[index].y,dist);
#endif
    }
  }

  /* Output image file */
  namebuf_size = strlen(basename)+32;
  namebuf = malloc(namebuf_size);

#if PGM_OUTPUT
  fprintf(stderr,"Output file: ");
#else
  if(flag_oldoutput) {
    fprintf(stderr,"Output files: ");
  } else {
    fprintf(stderr,"Output file: ");
  }
#endif
  if(strcmp("-",basename)==0) {
    strncpy(namebuf,"-",namebuf_size);
  } else {
    if(strchr(basename,'.')) {
      /* basename has extension already affixed. */
      strncpy(namebuf,basename,namebuf_size);
    } else {
      /* Append extension */
#if PGM_OUTPUT
      sprintf(namebuf,"%s.pgm",basename);
#else
      sprintf(namebuf,"%s.ppm",basename);
#endif
    }
  }

#if PGM_OUTPUT
  DumpToPGM(array, Nx, Ny, namebuf,
            (grain_bdry_color[0]+grain_bdry_color[1]+grain_bdry_color[2]+1)/3);
#else
  DumpToPPM(array,Nx,Ny,grain_count,namebuf,
            flag_nored,flag_nogreen,flag_noblue,
            grain_bdry_color[0],grain_bdry_color[1],grain_bdry_color[2]);
#endif
  fprintf(stderr,"%s",namebuf);
  free(namebuf);

#if !PGM_OUTPUT
  if(flag_oldoutput) {
    namebuf_mif = malloc(namebuf_size+4);
    if(strcmp("-",basename)==0) {
      strncpy(namebuf_mif,"-",namebuf_size);
    } else {
      if(strchr(basename,'.')) {
        /* basename has extension already affixed. */
        sprintf(namebuf_mif,"%s-mif",basename);
      } else {
        sprintf(namebuf_mif,"%s-mif.ppm",basename);
      }
    }
    DumpToPPM(array,Nx,Ny,0,namebuf_mif,
              flag_nored,flag_nogreen,flag_noblue,
              grain_bdry_color[0],grain_bdry_color[1],grain_bdry_color[2]);
    fprintf(stderr,", %s",namebuf_mif);
    free(namebuf_mif);
  }
#endif /* PGM_OUTPUT */

  fprintf(stderr,"\n");

  free(nuc);
  free(loc);
  free(pa);
  free(array);

  return return_code;
}