home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_6.2 / xvora / xvora.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  33.1 KB  |  1,256 lines
  1. /* 
  2.  * xvora.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * XVOR(onoi)A(nalysis)
  11.  *
  12.  * routine to extract statistics of a given Voronoi diagram from
  13.  * the representation stored in a file of type .vdt; a histogram
  14.  * may be contructed from various quantities
  15.  *
  16.  * note: occasionally, very small V polygon edge lengths occur which
  17.  *    may not show on the display: currently, these are flagged in 
  18.  *    site_gem() (vora.c), but are kept for statistics; this may have
  19.  *    to be dealt with more carefully;
  20.  *
  21.  */
  22.  
  23. #include <stdio.h>
  24. #include <stdlib.h>
  25. #include <malloc.h>
  26. #include <math.h>
  27. #include <string.h>
  28. #include <stdarg.h>
  29. #include "xvora.h"
  30.  
  31. #define    NX        511
  32. #define    NY        511
  33. #define    NXNY        NX*NY             /* total area */
  34.  
  35. #define    WHITE        255
  36.  
  37. #define    LINE_LEN    80             /* max length of line in data file */
  38. #define    TITLE_LEN    30            /* max length of line title */
  39. #define    FN_LEN        10              /* length of data file name */
  40.  
  41. #define    BIN_NUMBER    10           /* no bins in histogram */
  42.  
  43. #define    EDGE_NUMBER    0
  44. #define    POLY_AREA    1
  45. #define    NND        2
  46. #define    EDGE_HIST    "\nhistogram for polygon edges"
  47. #define    AREA_HIST    "\nhistogram for polygon areas"
  48. #define    NND_HIST    "\nhistogram for nn-distances"
  49.  
  50. #define    ND        20
  51.  
  52.  
  53. #define    ON         1
  54. #define    OFF         0
  55. #define    DISPL_ALL    ON            /* when set, display all sites */
  56.  
  57. /*
  58.  * #define DISC_ONLY 1
  59.  * #define HEXA_ONLY 2
  60.  */
  61.  
  62. #define    FOUR        4
  63. #define    FIVE        5                 /* defect site coordination */
  64. #define    SIX        6
  65. #define    SEVEN        7
  66. #define    EIGHT        8
  67.  
  68. #define    NSQ        2                  /* limit distrib to NSQ*std_dev */
  69.  
  70. #undef    SHOW_HDATA
  71. #undef    SHOW_HSORT
  72. #define    SHOW_Q
  73. #define    SHOW_P
  74.  
  75. #undef    FRMS_DEBUG
  76.  
  77. /* global variables */
  78. extern char *optarg;
  79. extern int optind, opterr;
  80.  
  81. int V_INP = 0;                  /* input file type flags */
  82. int D_INP = 0;
  83.  
  84. int ACTIVE_AOI = 1;             /* enforce ``active'' AOI bounds */
  85.      /* the ``active'' AOI represents */
  86.      /* a subset of the orig displ AOI */
  87.  
  88.  
  89. char *HIST_HEADER;
  90. int HIST_TYPE = -1;
  91. int N_BINS = BIN_NUMBER;
  92. int MIN_SET = 0;
  93. int MAX_SET = 0;
  94. int SHOW_INPUT = 0;
  95. int EVAL_KNN = 0;
  96. int SHOW_ALL = 1;
  97. int COMPRESS = 1;               /* compr data set before eval histo */
  98.  
  99. int LABEL_MODE = -1;
  100. int NBRS = 0;
  101. int CLRP = 1;
  102. int BORDER = 0;
  103. int MARK_AOI = 0;
  104.  
  105. int P_STATS = 0;
  106. int Q_STATS = 0;
  107.  
  108. int WRITE_FILE = 0;
  109.  
  110. int SAVE_DATA = 0;
  111. int DISPL_SDATA = 0;
  112. int DISPL_PDATA = 0;
  113. int DUMP_VC = 1;                /* dump poly vertex coords to file */
  114. int CPV = 1;                    /* if set, coll poly vertex coords */
  115.  
  116. int KNN_MAX = 2;                /* max nof NN shells to be constr */
  117.  
  118.  
  119.  
  120. /*
  121.  * error message
  122.  */
  123. void
  124. exitmess (prompt, status)
  125.      char *prompt;
  126.      int status;
  127. {
  128.   printf (prompt);
  129.   printf ("\n");
  130.  
  131.   exit (status);
  132. }
  133.  
  134.  
  135. void
  136. fail_alloc (char *str, int code)
  137. {
  138.   printf ("\n...memory alloc for %s failed\n", str);
  139.   exit (code);
  140. }
  141.  
  142.  
  143.  
  144. /*
  145.  * gprintf() functions:
  146.  * write to std output and, if option WRITE_FILE set, to file wbuf (stream)
  147.  */
  148. void
  149. gprintf (FILE * fpOut, char *fmt,...)
  150. {
  151.   va_list arg_ptr;
  152.  
  153.   va_start (arg_ptr, fmt);
  154.   vprintf (fmt, arg_ptr);
  155.   if (WRITE_FILE == 1)
  156.     vfprintf (fpOut, fmt, arg_ptr);
  157.   va_end (arg_ptr);
  158.  
  159. }
  160. /*
  161.  * usage of routine
  162.  */
  163. void
  164. usage (char *progname)
  165. {
  166.   progname = last_bs (progname);
  167.   printf ("USAGE: %s infile outimg [-d] [-k S] [-p a] [-b] [-e] [-a] [-l]\n", progname);
  168.   printf ("                   [-n nbins] [-i f] [-f f] [-s f] [-w f] [-L]\n");
  169.   printf ("\n%s extracts statistics of a given Voronoi diagram from\n", progname);
  170.   printf ("the representation stored in a file of type .vdt; a histogram\n");
  171.   printf ("may be contructed from various quantities\n\n");
  172.   printf ("ARGUMENTS:\n");
  173.   printf ("        infile: input filename (ASCII) produced by xvor;\n");
  174.   printf ("                .vdt (Voronoi data)\n");
  175.   printf ("             or .ddt (Delaunay data)\n");
  176.   printf ("        outimg: output image filename (TIF)\n\n");
  177.   printf ("OPTIONS:\n");
  178.   printf ("    -d: show input data\n");
  179.   printf ("  -k S: construct k-NN shells, k<=S (default=%d)\n", KNN_MAX);
  180.   printf ("  -p a: enable labeling; a is a string, as follows:\n");
  181.   printf ("        a=nall: enable numbering of all sites\n");
  182.   printf ("        a=n57: enable numbering of 5-fold, 7-fold sites\n");
  183.   printf ("        a=n57a: ->in addition, mark active AOI\n");
  184.   printf ("    -b: no border strip\n\n");
  185.   printf ("        construct histogram for:\n");
  186.   printf ("    -e: edge number (per Voronoi polygon)\n");
  187.   printf ("    -a: polygon area\n");
  188.   printf ("    -l: nn-distance (edge bisector length)\n");
  189.   printf ("  -n n: set number of bins to n (default:=%d)\n", N_BINS);
  190.   printf ("  -i f: set initial value to f(float)\n");
  191.   printf ("  -f f: set final value to f(float)\n");
  192.   printf ("    -s: evaluate nn-distance statistics function Q = Q(t)\n\n");
  193.   printf ("  -w f: write data to file\n");
  194.   printf ("        f=h: hist data --> file of type .hdt\n");
  195.   printf ("        f=s: site topol data --> file of type .std\n");
  196.   printf ("        f=p: poly topol data --> file of type .ptd\n");
  197.   printf ("    -L: print Software License for this module\n");
  198.   exit (1);
  199. }
  200.  
  201.  
  202. /* 
  203.  * sort sites on y
  204.  */
  205. int
  206. sycomp (s1, s2)
  207.      struct vSite *s1, *s2;
  208. {
  209.   return ((int) (SIGN (s1->coord.y - s2->coord.y)));
  210. }
  211.  
  212.  
  213.  
  214. /* 
  215.  * sort sites on x
  216.  */
  217. int
  218. sxcomp (s1, s2)
  219.      struct vSite *s1, *s2;
  220. {
  221.   return ((int) (SIGN (s1->coord.x - s2->coord.x)));
  222. }
  223.  
  224.  
  225.  
  226. /*
  227.  * comparison function for qsort():
  228.  * sort array of tuples in order of increassing x-values
  229.  */
  230. int
  231. compare (t1, t2)
  232.      const void *t1, *t2;
  233. {
  234.   float f1, f2;
  235.  
  236.   f1 = *((float *) t1);
  237.   f2 = *((float *) t2);
  238.   return ((int) SIGN (f1 - f2));
  239. }
  240.  
  241. /*
  242.  * label site at (x, y) with numeral
  243.  */
  244. void
  245. nbr_site (x, y, n, imgIO, value)
  246.      double x, y;
  247.      int n;
  248.      Image *imgIO;
  249.      int value;
  250. {
  251.   static char *num[20] =
  252.   {" 0", " 1", " 2", " 3", " 4",
  253.    " 5", " 6", " 7", " 8", " 9",
  254.    "10", "11", "12", "13", "14",
  255.    "15", "16", "17", "18", "19"};
  256.   int ix, iy;
  257.  
  258.   ix = (int) x - 16;
  259.   iy = (int) y + 12;
  260.   if (DISC_ONLY == ON) {
  261.     if ((strcmp (*(num + n), " 5") == 0) ||
  262.         (strcmp (*(num + n), " 7") == 0))
  263.       gdImageString (gdFontSmall, (int) x, (int) y, *(num + n), imgIO, value);
  264.   }
  265.   else
  266.     gdImageString (gdFontSmall, (int) x, (int) y, *(num + n), imgIO, value);
  267. }
  268.  
  269.  
  270. void
  271. main (int argc, char **argv)
  272. {
  273.   int i, ich;
  274.  
  275.   int i_arg;
  276.   static char whb[32];
  277.   static char wsb[32];
  278.   static char wpb[32];
  279.  
  280.   static char *inp_suffix =
  281.   {".xxx"};                     /* default inp file suffix */
  282.   static char *vdt_type =
  283.   {".vdt"};                     /* suffix for .vdt inp file */
  284.   static char *ddt_type =
  285.   {".ddt"};                     /* suffix for .ddt inp file */
  286.   static char *whsuffix =
  287.   {".hdt"};                     /* suffix for output file name */
  288.   static char *wssuffix =
  289.   {".std"};                     /* suffix for output file name */
  290.   static char *wpsuffix =
  291.   {".ptd"};                     /* suffix for output file name */
  292.   char *rbuf;
  293.  
  294.   char ln_buf[LINE_LEN];
  295.  
  296.   int ic, is, ie;
  297.   int k;
  298.   int inn, nnn_max;
  299.   int ns, ne, nv;
  300.   int nas, ncl = FOUR, ncu = EIGHT;
  301.   int *nnna;
  302.   double c_6, c_ic, mu_2;
  303.   char rec;
  304.  
  305.   struct vSite *s;              /* array of Voronoi sites */
  306.   struct vSite *v;              /* array of Voronoi vertices */
  307.   struct vEdge *e;              /* array of Voronoi edges */
  308.   struct vPoly *p;              /* array of Voronoi polygons */
  309.   struct Tuple *pem;            /* array of ptrs to Tuple: */
  310.   /* midpoints of poly edges */
  311.  
  312.   unsigned int **a_n;           /* ar of pntrs to area of n-fold coord bubs */
  313.   unsigned int *a_nbar;         /* mean area for bubs of coord no n */
  314.   int nel = FOUR, neu = EIGHT;
  315.   int ipn, *npna;               /* array cont nof polygons with n edges */
  316.   double ma_vp;
  317.  
  318.   int xmin, xmax, ymin, ymax;
  319.  
  320.   int ih, nh;
  321.   float bw, min = (float) -1.0, max = (float) -1.0;
  322.   float *hd, *ha;
  323.   double del, mean, sdev;
  324.   struct Histo histo, *h = &histo;
  325.  
  326.   float *foo;
  327.   int np = 3;
  328.  
  329.   double r0;
  330.   struct Tuple *pf, *qf, *dstats;
  331.  
  332.   /* list structs */
  333.   struct kNNShell *sk;          /* ary of pntrs to lists of kNN shells */
  334.   struct linklist *sha;         /* array of kNN shells (lists) */
  335.   double frms;
  336.   double *tctck;
  337.  
  338.   FILE *fpIn, *fpOut;
  339.  
  340.   Image *imgOut;
  341.  
  342. /* 
  343.  * cmd line options:
  344.  */
  345.   static char *optstring = "dbp:ealn:i:f:sk:w:o:L";
  346.  
  347.  
  348.   if (argc < 3)
  349.     usage (argv[0]);
  350.  
  351. /* get the filename for input */
  352.   printf ("read data from file %s\n", rbuf = argv[1]);
  353.   if ((fpIn = fopen (rbuf, "r")) == NULL) {
  354.     printf ("\n...cannot open input file: %s\n", rbuf);
  355.     exit (1);
  356.   }
  357.  
  358.   /* construct output file name */
  359.   ich = 0;
  360.   while ((*(whb + ich) = *(wsb + ich) = *(wpb + ich) = *(rbuf + ich)) != '.')
  361.     ich++;
  362.   /* strip input file suffix */
  363.   for (is = 0; is < 4; is++)
  364.     *(inp_suffix + is) = *(rbuf + ich + is);
  365.  
  366.   if (strcmp (inp_suffix, vdt_type) == 0)
  367.     V_INP = 1;
  368.   else if (strcmp (inp_suffix, ddt_type) == 0)
  369.     D_INP = 1;
  370.   else {
  371.     printf ("\n...unknown input file type encountered:");
  372.     printf (" legal types: %s %s\n", vdt_type, ddt_type);
  373.     fclose (fpIn);
  374.     exit (1);
  375.   }
  376.  
  377.   /* append suffix to output file names */
  378.   for (is = 0; is < 4; is++) {
  379.     *(whb + ich + is) = *(whsuffix + is);
  380.     *(wsb + ich + is) = *(wssuffix + is);
  381.     *(wpb + ich + is) = *(wpsuffix + is);
  382.   }
  383.  
  384. /*
  385.  * parse command line
  386.  */
  387.   optind = 3;
  388.   opterr = ON;                  /* give error messages */
  389.  
  390.  
  391.   while ((i_arg = getopt (argc, argv, optstring)) != EOF) {
  392.     switch (i_arg) {
  393.     case 'd':
  394.       printf ("...option %c: ", i_arg);
  395.       printf ("display input data\n");
  396.       SHOW_INPUT = ON;
  397.       break;
  398.     case 'k':
  399.       printf ("...option %c: ", i_arg);
  400.       printf ("construct %d-NN shells for all Sites\n", KNN_MAX = atoi (optarg));
  401.       EVAL_KNN = ON;
  402.       break;
  403.     case 'p':
  404.       printf ("...option %c: ", i_arg);
  405.       if (strcmp (optarg, "nall") == 0) {
  406.         printf ("number all V. sites\n");
  407.         NBRS = ON;
  408.         LABEL_MODE = DISC_ONLY;
  409.         MARK_AOI = OFF;
  410.       }
  411.       else if (strcmp (optarg, "n57") == 0) {
  412.         printf ("number V. disclination sites\n");
  413.         NBRS = ON;
  414.         LABEL_MODE = DISC_ONLY;
  415.         MARK_AOI = OFF;
  416.       }
  417.       else if (strcmp (optarg, "n57a") == 0) {
  418.         printf ("number V. disclination sites\n");
  419.         NBRS = ON;
  420.         LABEL_MODE = DISC_ONLY;
  421.         MARK_AOI = ON;
  422.       }
  423.       else if (strcmp (optarg, "c57") == 0) {
  424.         printf ("color V. penta- and heptagons\n");
  425.         CLRP = ON;
  426.         LABEL_MODE = DISC_ONLY;
  427.         BORDER = ON;
  428.         MARK_AOI = OFF;
  429.       }
  430.       else if (strcmp (optarg, "c57a") == 0) {
  431.         printf ("color V. penta- and heptagons\n");
  432.         CLRP = ON;
  433.         LABEL_MODE = DISC_ONLY;
  434.         BORDER = ON;
  435.         MARK_AOI = ON;
  436.       }
  437.       else
  438.         exitmess ("\n...unknown p option", 1);
  439.       break;
  440.     case 'b':
  441.       printf ("...option %c: ", i_arg);
  442.       printf ("no zero border strip\n");
  443.       if (BORDER == ON)
  444.         BORDER = OFF;
  445.       break;
  446.     case 'e':
  447.       printf ("...option %c: ", i_arg);
  448.       printf ("construct histogram for poly edges\n");
  449.       HIST_TYPE = EDGE_NUMBER;
  450.       HIST_HEADER = EDGE_HIST;
  451.       break;
  452.     case 'a':
  453.       printf ("...option %c: ", i_arg);
  454.       printf ("construct histogram for polygon area\n");
  455.       HIST_TYPE = POLY_AREA;
  456.       HIST_HEADER = AREA_HIST;
  457.       break;
  458.     case 'l':
  459.       printf ("...option %c: ", i_arg);
  460.       printf ("construct histogram for vertex nn-dist\n");
  461.       HIST_TYPE = NND;
  462.       HIST_HEADER = NND_HIST;
  463.       break;
  464.     case 'n':
  465.       printf ("...option %c: ", i_arg);
  466.       printf ("set number of bins to %d\n",
  467.               N_BINS = atoi (optarg));
  468.       break;
  469.     case 'i':
  470.       printf ("...option %c: ", i_arg);
  471.       printf ("set min (init value) to %f\n",
  472.               min = (float) atof (optarg));
  473.       MIN_SET = 1;
  474.       break;
  475.     case 'f':
  476.       printf ("...option %c: ", i_arg);
  477.       printf ("set max (final value) to %f\n",
  478.               max = (float) atof (optarg));
  479.       MAX_SET = 1;
  480.       break;
  481.     case 's':
  482.       printf ("...option %c: ", i_arg);
  483.       printf ("distance statistics for q(t)\n");
  484.       Q_STATS = 1;
  485.       break;
  486.     case 'w':
  487.       printf ("...option %c: ", i_arg);
  488.       if (strcmp (optarg, "h") == 0) {
  489.         if ((fpOut = fopen (whb, "w")) == NULL) {
  490.           printf ("\n...cannot open %s for writing\n",
  491.                   whb);
  492.           exit (1);
  493.         }
  494.         printf ("write hist data to file %s\n", whb);
  495.         WRITE_FILE = 1;
  496.       }
  497.       else if (strcmp (optarg, "s") == 0) {
  498.         if ((fpOut = fopen (wsb, "w")) == NULL) {
  499.           printf ("\n...cannot open %s for writing\n",
  500.                   wsb);
  501.           exit (1);
  502.         }
  503.         printf ("log site topol data in file: %s\n", wsb);
  504.         DISPL_SDATA = 1;
  505.         WRITE_FILE = 1;
  506.       }
  507.       else if (strcmp (optarg, "p") == 0) {
  508.         if ((fpOut = fopen (wpb, "w")) == NULL) {
  509.           printf ("\n...cannot open %s for writing\n",
  510.                   wpb);
  511.           exit (1);
  512.         }
  513.         printf ("log polygon topol data in file: %s\n", wpb);
  514.         DISPL_PDATA = 1;
  515.         WRITE_FILE = 1;
  516.       }
  517.       else
  518.         exitmess ("\n...unknown -w option", 1);
  519.  
  520.       break;
  521.     case 'L':
  522.       print_sos_lic ();
  523.       exit (0);
  524.     default:
  525.       printf ("\n...unknown cmd line argument\n");
  526.       exit (1);
  527.       break;
  528.     }
  529.   }
  530.  
  531.   if ((Q_STATS == 1) && (HIST_TYPE != NND)) {
  532.     printf ("\n...before evaluating q-function:\n");
  533.     printf ("-->evaluate NND histogram (opt -l) with %d bins (opt -nND)\n", ND + 1);
  534.     exit (1);
  535.   }
  536.   if (D_INP == 1) {
  537.     printf ("\n...Delaunay input (.ddt) currently not handled\n");
  538.     exit (1);
  539.   }
  540.  
  541. /*
  542.  * first pass through .vdt file:
  543.  * -->determine number of sites (identical to the number of Voronoi polygons) 
  544.  * -->determine number of edges
  545.  */
  546.   ns = ne = nv = 0;
  547.   do {
  548.     while ((fgets (ln_buf, LINE_LEN, fpIn)) != (char *) NULL) {
  549.       sscanf (ln_buf, "%c", &rec);
  550.       switch (rec) {
  551.       case 's':
  552.         ns++;
  553.         break;
  554.       case 'v':
  555.         nv++;
  556.         break;
  557.       case 'e':
  558.         ne++;
  559.         break;
  560.       default:
  561.         break;
  562.       }
  563.     }
  564.   } while (feof (fpIn) == 0);
  565.   fclose (fpIn);
  566.   printf ("\n...first pass through data found:\n");
  567.   printf ("   %d sites, %d vertices, %d edges\n", ns, nv, ne);
  568.  
  569. /*
  570.  * allocate memory for input data
  571.  */
  572.   if ((s = (struct vSite *) calloc (ns, sizeof (struct vSite))) == NULL)
  573.       fail_alloc ("struct vSite *s", 1);
  574.   if ((v = (struct vSite *) calloc (nv, sizeof (struct vSite))) == NULL)
  575.       fail_alloc ("struct vSite *v", 1);
  576.   if ((e = (struct vEdge *) calloc (ne, sizeof (struct vEdge))) == NULL)
  577.       fail_alloc ("struct vEdge *e", 1);
  578.   if ((p = (struct vPoly *) calloc (ns, sizeof (struct vPoly))) == NULL)
  579.       fail_alloc ("struct vPoly *p", 1);
  580.  
  581. /*
  582.  * scan input file
  583.  */
  584.   if ((fpIn = fopen (rbuf, "r")) == NULL)
  585.     exitmess ("can't open input file", 1);
  586.  
  587.   printf ("\n...read and classify input data\n");
  588.   ns = nv = ne = 0;
  589.   do {
  590.     while ((fgets (ln_buf, LINE_LEN, fpIn)) != (char *) NULL) {
  591.       sscanf (ln_buf, "%c", &rec);
  592.       switch (rec) {
  593.       case 's':
  594.         sscanf (ln_buf, "%c %f %f",
  595.                 &rec, &((s + ns)->coord.x), &((s + ns)->coord.y));
  596.         (s + ns)->sitenbr = (p + ns)->sitenbr = ns;
  597.         ns++;
  598.         break;
  599.       case 'v':
  600.         sscanf (ln_buf, "%c %f %f",
  601.                 &rec, &((v + nv)->coord.x), &((v + nv)->coord.y));
  602.         (v + nv)->sitenbr = nv;
  603.         nv++;
  604.         break;
  605.       case 'e':
  606.         sscanf (ln_buf, "%c %d %d %d %d %d",
  607.                 &rec, &((e + ne)->edgenbr),
  608.                 &((e + ne)->vO), &((e + ne)->vF),
  609.                 &((e + ne)->sO), &((e + ne)->sF));
  610.         ne++;
  611.         break;
  612.       default:
  613.         break;
  614.       }
  615.     }
  616.   } while (feof (fpIn) == 0);
  617.   fclose (fpIn);
  618.  
  619.   if (SHOW_INPUT == ON) {
  620.     printf ("\n...second pass through data:\n");
  621.  
  622.     printf ("...sites:\n");
  623.     for (i = 0; i < ns; i++)
  624.       printf ("  %3d:  x: %f   y: %f\n",
  625.               (s + i)->sitenbr, (s + i)->coord.x, (s + i)->coord.y);
  626.  
  627.     printf ("...vertices:\n");
  628.     for (i = 0; i < nv; i++)
  629.       printf ("  %3d:  x: %f   y: %f\n",
  630.               (v + i)->sitenbr, (v + i)->coord.x, (v + i)->coord.y);
  631.  
  632.     printf ("...edges:\n");
  633.     for (i = 0; i < ne; i++)
  634.       printf ("  e %3d:  vO: %3d vF: %3d  sO: %3d sF: %3d\n",
  635.       (e + i)->edgenbr, (e + i)->vO, (e + i)->vF, (e + i)->sO, (e + i)->sF);
  636.   }
  637.  
  638.  
  639. /*
  640.  * check whether Site falls outside of the ``active'' AOI
  641.  * as defined in vora.h: if out of bounds, set aoiFlag Active;
  642.  */
  643.   mark_eSites (s, ns);
  644.  
  645. /*
  646.  * enforce ``active'' AOI bounds: mark all edges which cross 
  647.  * AOI boundaries and flag corresponding sites which lie out-of-bounds
  648.  */
  649.   if (ACTIVE_AOI == 1)
  650.     mark_eEdges (v, e, ne);
  651.  
  652. /*
  653.  * determine site coordination, 
  654.  * given by the number of edges of the associated Voronoi polygon
  655.  *
  656.  * if any of the polygon edges associated with a given Site 
  657.  * end out-of-bounds, mark that Site by setting eFlag: this
  658.  * flag indicates enforcement of stronger condition than that
  659.  * involved in setting aoiFlag (in mark_eSite())
  660.  *
  661.  * (the offending vertex should already have been marked, 
  662.  * by setting aoiFlag, in mark_eEdges())
  663.  */
  664.  
  665. /* 
  666.  * initialize structure member 
  667.  */
  668.   for (is = 0; is < ns; is++) {
  669.     (s + is)->nnn = 0;
  670.     (s + is)->eFlag = UnSet;
  671.   }
  672.   site_coordination (s, v, e, ns, ne);
  673.  
  674.  
  675. /*
  676.  * for each site, of known coord. number:
  677.  * -->determine the set of nearest neighbor distances;
  678.  * -->construct Voronoi polygons associated with given site by grouping edges;
  679.  * 
  680.  * flag ``edge'' polygons, by checking the status of 
  681.  * 
  682.  * set eFlag if the polygon does not lie entirely within the ``active'' AOI; 
  683.  * that is, enforce the criterion, more stringent than the one employed in 
  684.  * function mark_eSites(), that all polygon edges associated with Sites 
  685.  * considered to be part of the ``active'' AOI must lie entirely within that 
  686.  * (new) AOI; 
  687.  *
  688.  * in addition, if any polygon edge length is smaller than a preset limit
  689.  * (see EPS in vora.c), set lFlag;
  690.  */
  691.   for (i = 0; i < ns; i++) {
  692.     if (((s + i)->nnd = (float *) calloc ((s + i)->nnn, sizeof (float))) == NULL)
  693.         fail_alloc ("(s+i)->nnd", 1);
  694.  
  695.     if (((s + i)->nnsi = (int *) calloc ((s + i)->nnn, sizeof (int))) == NULL)
  696.         fail_alloc ("(s+i)->nnsi", 1);
  697.   }
  698.  
  699.   for (i = 0; i < ns; i++) {
  700.     if (((p + i)->pei = (int *) calloc ((s + i)->nnn, sizeof (int))) == NULL)
  701.         fail_alloc ("(p+i)->pei", 1);
  702.  
  703.     if (((p + i)->vel = (float *) calloc ((s + i)->nnn, sizeof (float))) == NULL)
  704.         fail_alloc ("(s+i)->nnd", 1);
  705.  
  706.     if (((p + i)->sphi = (float *) calloc ((s + i)->nnn, sizeof (float))) == NULL)
  707.         fail_alloc ("(p+i)->sphi", 1);
  708.  
  709.     if (CPV == 1) {
  710.       if (((p + i)->pvi = (int *) calloc ((s + i)->nnn, sizeof (int))) == NULL)
  711.           fail_alloc ("(p+i)->pvi", 1);
  712.     }
  713.   }
  714.   printf ("...evaluate properties of sites\n");
  715.   nnn_max = site_geom (s, v, e, p, ns, ne);
  716.  
  717. /*
  718.  * allocate memory for auxiliary array of edge midpoint coordinates,
  719.  * employed in sorting of polygon edges
  720.  */
  721.   if ((pem = (struct Tuple *) calloc (nnn_max, sizeof (struct Tuple))) == NULL)
  722.       fail_alloc ("struct Tuple *pem", 1);
  723.  
  724. /* 
  725.  * during construction of the V. diagram (xvor.c) edge endpoints 
  726.  * which lie out-of-bounds (with respect to the full display AOI) are
  727.  * flagged by setting the corresponding vertex array index (vO or vF)
  728.  * to -1; that is, in the array e, indices (e+is)->vO, (e+is)->vF) of 
  729.  * vertices delimiting edges are set to -1 if corresponding coordinates 
  730.  * are out-of-bounds with respect to the original AOI
  731.  */
  732.   printf ("...evaluate properties of polygons\n");
  733.   poly_geom (s, v, e, p, pem, ns, ne, CPV);
  734.  
  735.   if (DISPL_PDATA == 1) {
  736.     log_pt_data (fpOut, wpb, ns, p, e, v, CPV);
  737.     if (DUMP_VC == 1) {
  738.       printf ("\n...dump polygon vertices to file\n");
  739.       dump_poly_vc (fpOut, ns, p, v);
  740.     }
  741.   }
  742.   else {
  743.     printf ("\n...to display and log polygon data: option -wp\n");
  744.   }
  745.  
  746.  
  747.  
  748. /* 
  749.  * examine nof NNs as a function of V polygon domain area (->Lewis' law)
  750.  * (see also: xsgt.c)
  751.  */
  752.   npna = pe_dist (p, ns, nel, neu);
  753.   gprintf (fpOut, "\n...distribution of V poly edges:\n");
  754.   for (ie = 0; ie <= nnn_max; ie++) {
  755.     gprintf (fpOut, "\t%2d-sided V polygons : %3d\n", ie, *(npna + ie));
  756.   }
  757.  
  758.   if ((a_n = (unsigned int **) calloc (nnn_max + 1, sizeof (unsigned int *))) == NULL)
  759.       fail_alloc ("a_n", 1);
  760.  
  761.   for (inn = 0; inn <= nnn_max; inn++) {
  762.     if ((a_n[inn] = (unsigned int *) calloc (npna[inn] + 1, sizeof (unsigned int))) == NULL)
  763.         fail_alloc ("a_n[inn]", 1);
  764.  
  765.   }
  766.  
  767.   if ((a_nbar = (unsigned int *) calloc (nnn_max + 1, sizeof (unsigned int))) == NULL)
  768.       fail_alloc ("a_nbar", 1);
  769.  
  770.   for (ie = 0; ie <= nnn_max; ie++)
  771.     *(npna + ie) = 0;
  772.   ma_vp = p_of_nVA (npna, a_n, ns, p);
  773.   //for(inn=0; inn<nnn_max+10; inn++)
  774.   //  free(a_n[inn]);
  775.   //free(a_n);
  776.   //exit(0);
  777.  
  778.   gprintf (fpOut, "\nmean area of (active) Vor poly: %lf\n", ma_vp);
  779.   gprintf (fpOut, "V poly grouped by edge no n=2(1)%d:\n", nnn_max);
  780.   for (inn = 2; inn <= nnn_max; inn++) {
  781.     for (ipn = 0; ipn < npna[inn]; ipn++) {
  782.       if (SHOW_ALL == 1) {
  783.         gprintf (fpOut, "...%5u ", *(a_n[inn] + ipn));
  784.         if ((ipn + 1) % 8 == 0)
  785.           gprintf (fpOut, "\n");
  786.       }
  787.       *(a_nbar + inn) += *(a_n[inn] + ipn);
  788.     }
  789.     if (ipn > 0)
  790.       *(a_nbar + inn) /= ipn;
  791.  
  792.     gprintf (fpOut, "\n nof NN: %2d -- <A_n> (mean of %d): %6u ",
  793.              inn, ipn, *(a_nbar + inn));
  794.     gprintf (fpOut, " -- <A_n>/<A>: %lf\n\n",
  795.              (double) (*(a_nbar + inn)) / ma_vp);
  796.   }
  797.  
  798.  
  799.  
  800. /* copy vPoly area into vSite array (to ease handling of output) */
  801.   for (is = 0; is < ns; is++)
  802.     (s + is)->v_area = (unsigned int) (p + is)->area;
  803.  
  804.   if (DISPL_SDATA == 1)
  805.     log_st_data (fpOut, wsb, ns, s);
  806.  
  807.   else {
  808.     printf ("\n...to display and log site data: option -ws\n");
  809.   }
  810.  
  811. /*
  812.  * determine the number of topological defects:
  813.  * this is given by the number of sites not coordinated with 6 NN; 
  814.  * restrict count to an ``active'' AOI excluding the ``edge'' of 
  815.  * the display AOI;
  816.  */
  817.   printf ("   defect sites in active AOI (%3d, %3d), (%3d, %3d)\n",
  818.           AULX, AULY, ALRX, ALRY);
  819.  
  820.   nnna = count_defects (s, ns, ncl, ncu);
  821.   nas = 0;
  822.   for (ic = ncl; ic <= ncu; ic++)
  823.     nas += *(nnna + ic);
  824.  
  825.   gprintf (fpOut, "\n...distribution of NN:\n");
  826.   mu_2 = 0.0;
  827.   for (ic = ncl; ic <= ncu; ic++) {
  828.     c_ic = (double) (*(nnna + ic)) / (double) nas;
  829.     mu_2 += (ic - 6) * (ic - 6) * c_ic;
  830.  
  831.     gprintf (fpOut, "\t%2d-fold sites: %3d -- f_%2d: %lf\n",
  832.              ic, *(nnna + ic), ic, c_ic);
  833.   }
  834.   c_6 = (double) (*(nnna + SIX)) / (double) nas;
  835.   gprintf (fpOut, "...conc of 6-fold sites: ");
  836.   gprintf (fpOut, "\tc_6 = (n_6 = %d)/%d = %lf\n", *(nnna + SIX), nas, c_6);
  837.   gprintf (fpOut, "...second moment of NN distribution: %lf\n", mu_2);
  838.   gprintf (fpOut, "\n");
  839.  
  840.  
  841. /* 
  842.  * kNN shells and fractal measure
  843.  */
  844.   if (EVAL_KNN == 1) {
  845.     if ((sk = (struct kNNShell *) calloc (ns, sizeof (struct kNNShell))) == NULL)
  846.         fail_alloc ("sk", 1);
  847.  
  848.  
  849.     for (is = 0; is < ns; is++) {
  850.       if (((sk + is)->aka = (double *) calloc (KNN_MAX, sizeof (double))) == NULL)
  851.           fail_alloc ("sk->aka", 1);
  852.     }
  853.     for (is = 0; is < ns; is++) {
  854.       if (((sk + is)->tctc = (double *) calloc (KNN_MAX, sizeof (double))) == NULL)
  855.           fail_alloc ("sk->tctc", 1);
  856.     }
  857.  
  858.  
  859.     if ((sha = (struct linklist *) calloc (KNN_MAX, sizeof (struct linklist))) == NULL)
  860.         fail_alloc ("sha", 1);
  861.  
  862.  
  863.     gprintf (fpOut, "...constr kNN shells up to k=%d\n", KNN_MAX);
  864.     frms = construct_kNNs (ns, s, sha, sk);
  865.     gprintf (fpOut, "...avrg fractal meas of set: %lf\n", frms);
  866.  
  867.  
  868. /*
  869.  * display arrays sk->tctc, then average values over sites
  870.  */
  871.     if ((tctck = (double *) calloc (KNN_MAX + 1, sizeof (double))) == NULL)
  872.         fail_alloc ("tctck", 1);
  873.  
  874.  
  875.     gprintf (fpOut, "\n...charge-charge correlations (active Sites only)\n");
  876.  
  877.     nas = 0;
  878.     for (k = 0; k < KNN_MAX; k++)
  879.       *(tctck + k) = 0.0;
  880.     for (is = 0; is < ns; is++) {
  881.       gprintf (fpOut, "...Site %3d ", (sk + is)->sitenbr);
  882.       gprintf (fpOut, "charge: %d ", (sk + is)->tc);
  883.  
  884.       if (((s + is)->eFlag == UnSet) &&
  885.           ((s + is)->aoiFlag == UnSet)) {
  886.         nas++;
  887.         for (k = 0; k < (sk + is)->kNN; k++) {
  888.           gprintf (fpOut, "...%6.2lf ",
  889.                    (sk + is)->tctc[k]);
  890.  
  891.           *(tctck + k) += (sk + is)->tctc[k];
  892.         }
  893.         if ((k + 1) % 8 == 0)
  894.           gprintf (fpOut, "\n");
  895.       }
  896.       gprintf (fpOut, "\n");
  897.     }
  898.     gprintf (fpOut, "\n...charge-charge correlation function:\n");
  899.     gprintf (fpOut, "\taveraged over %d active Sites\n", nas);
  900.     for (k = 0; k < KNN_MAX; k++) {
  901.       gprintf (fpOut, "...%6.2lf", *(tctck + k) / (double) nas);
  902.  
  903.       if ((k + 1) % 8 == 0)
  904.         gprintf (fpOut, "\n");
  905.     }
  906.     gprintf (fpOut, "\n");
  907.  
  908. #ifdef FRMS_DEBUG
  909.     KNN_MAX = 0;
  910.     for (is = 0; is < ns; is++) {
  911.       printf ("Site %3d: ", (s + is)->sitenbr);
  912.       if (((s + is)->eFlag == UnSet) &&
  913.           ((s + is)->aoiFlag == UnSet)) {
  914.         printf ("kNN: %3d, frms: %6.3lf\n",
  915.                 (sk + is)->kNN, (sk + is)->frms);
  916.  
  917.         if ((sk + is)->kNN > KNN_MAX)
  918.           KNN_MAX = (sk + is)->kNN;
  919.       }
  920.       else {                    /* Site flagged out-of-bounds */
  921.         if ((sk + is)->kNN == 0)
  922.           printf (" out of bounds: -->ignored\n");
  923.         else
  924.           printf (" something wrong\n");
  925.       }
  926.       printf ("\tactual KNN_MAX: %3d\n", KNN_MAX);
  927.     }
  928. #endif
  929.   }
  930.  
  931.  
  932.  
  933. /*
  934.  * histogram evaluation
  935.  */
  936.   if (HIST_TYPE != -1) {
  937.     printf ("\n...initiate histogram evaluation\n");
  938.     switch (HIST_TYPE) {
  939.     case EDGE_NUMBER:          /* site coordination */
  940.       if ((hd = (float *) calloc (ns, sizeof (float))) == NULL)
  941.           fail_alloc ("hd", 1);
  942.       nh = 0;
  943.       for (is = 0; is < ns; is++) {
  944.         if (((s + is)->eFlag == UnSet) &&
  945.             ((s + is)->aoiFlag == UnSet)) {
  946.           *(hd + is) = (float) (s + is)->nnn;
  947.           nh++;
  948.         }
  949.       }
  950.       if ((hd = (float *) realloc (hd, nh * sizeof (float))) == NULL)
  951.           fail_alloc ("hd (realloc)", 1);
  952.       COMPRESS == 0;
  953.       break;
  954.     case POLY_AREA:
  955.       nh = ns;
  956.       if ((hd = (float *) calloc (nh, sizeof (float))) == NULL)
  957.           fail_alloc ("hd", 1);
  958.       for (i = 0; i < nh; i++) {
  959.         if ((*(hd + i) = (float) (p + i)->area) == -1.0)
  960.           *(hd + i) = (float) 0.0;
  961.       }
  962.       break;
  963.     case NND:                  /* note: each nn bond counted twice !! */
  964.       nh = 0;
  965.       for (is = 0; is < ns; is++)
  966.         nh += (s + is)->nnn;
  967.       if ((hd = (float *) calloc (nh, sizeof (float))) == NULL)
  968.           fail_alloc ("hd", 1);
  969.       i = 0;
  970.       for (is = 0; is < ns; is++) {
  971.         for (ie = 0; ie < (s + is)->nnn; ie++) {
  972.           *(hd + i) = (s + is)->nnd[ie];
  973.           i++;
  974.         }
  975.       }
  976.       break;
  977.     }
  978.  
  979. #ifdef SHOW_HDATA
  980.     printf ("\n...data to construct %s:\n", HIST_HEADER);
  981.     for (i = 0; i < nh; i++)
  982.       printf (" data[%d] = %f\n", i, *(hd + i));
  983.     printf ("\n");
  984. #endif
  985.  
  986. /*
  987.  * evaluate statistics
  988.  */
  989.     mean = find_mean (hd, nh);
  990.     sdev = find_sigma (hd, nh, mean);
  991.     printf ("...(raw) mean: %lf, std dev: %lf\n", mean, sdev);
  992.  
  993. /*
  994.  * limit values to be included in histogram to those within a spread of 
  995.  * NSQ*sdev of (raw) mean to eliminate artefacts
  996.  */
  997.     if (COMPRESS == 1) {
  998.       ih = 0;
  999.       for (i = 0; i < nh; i++) {
  1000.         del = (double) (*(hd + i)) - mean;
  1001.         if (fabs (del) <= NSQ * sdev) {
  1002.           *(hd + ih) = *(hd + i);
  1003.           ih++;
  1004.         }
  1005.       }
  1006.       nh = ih;
  1007.       printf ("...%d entries within %d*std dev of mean:\n", nh, NSQ);
  1008.  
  1009. /* 
  1010.  * realloc, to reduce memory requirement
  1011.  */
  1012.       hd = (float *) realloc (hd, nh * sizeof (float));
  1013.     }
  1014.  
  1015. /*
  1016.  * prepare to evaluate histogram
  1017.  */
  1018.     qsort (hd, nh, sizeof (float), compare);
  1019.  
  1020. #ifdef SHOW_HSORT
  1021.     printf ("\n...sorted data:\n");
  1022.     for (i = 0; i < nh; i++)
  1023.       printf (" data[%d] = %f\n", i, *(hd + i));
  1024. #endif
  1025.  
  1026.     if (MIN_SET == 0)
  1027.       min = *(hd + 0);
  1028.     if (MAX_SET == 0)
  1029.       max = *(hd + nh - 1);
  1030.     if ((-1.0e-10 < (min - max)) && ((min - max) < 1.0e-10))
  1031.       min = max;
  1032.  
  1033.     if (HIST_TYPE == EDGE_NUMBER) {  /* discrete variable */
  1034.       min = (float) 1.0;
  1035.       max = *(hd + nh - 1);
  1036.       N_BINS = (int) (max - min) + 1;
  1037.     }
  1038.     bw = (max - min) / ((float) (N_BINS - 1));
  1039.  
  1040.     if ((ha = (float *) calloc ((size_t) N_BINS, sizeof (float))) == NULL)
  1041.         fail_alloc ("hist\n", 1);
  1042.  
  1043.     if (bw != 0.0)
  1044.       construct_shist (nh, hd, ha, bw, min);
  1045.  
  1046.  
  1047.     printf ("%s\n", HIST_HEADER);
  1048.     printf ("......mean: %lf\n", mean);
  1049.     printf ("...std_dev: %lf\n", sdev);
  1050.     printf ("...min: %f, bw: %f, max: %f\n", min, bw, max);
  1051.     for (i = 0; i < N_BINS; i++) {
  1052.       printf (" %6.2f", *(ha + i));
  1053.       if ((i + 1) % 10 == 0)
  1054.         printf ("\n");
  1055.     }
  1056.   }
  1057.  
  1058. #if defined(LINUX)
  1059.   qsort (s, (size_t) ns, sizeof (struct vSite), (__compar_fn_t) sycomp);
  1060. #else
  1061.   qsort (s, (size_t) ns, sizeof (struct vSite), sycomp);
  1062. #endif
  1063.   ymin = (int) (s + 0)->coord.y;
  1064.   ymax = (int) (s + ns - 1)->coord.y;
  1065.  
  1066. #if defined(LINUX)
  1067.   qsort (s, ns, sizeof (struct vSite), (__compar_fn_t) sxcomp);
  1068. #else
  1069.   qsort (s, ns, sizeof (struct vSite), sxcomp);
  1070. #endif
  1071.   xmin = (int) (s + 0)->coord.x;
  1072.   xmax = (int) (s + ns - 1)->coord.x;
  1073.  
  1074.   printf ("...min: (%3d, %3d), max: (%3d, %3d)\n",
  1075.           xmin, ymin, xmax, ymax);
  1076. /* allocate default image buffer for image output */
  1077.   imgOut = ImageAlloc ((long) ymax, (long) xmax, BPS);
  1078.  
  1079. /*
  1080.  * label Voronoi diagram
  1081.  */
  1082.   if (NBRS == ON) {
  1083.     if (MARK_AOI == ON)
  1084.       draw_rect (AULX, AULY, ALRX, ALRY, imgOut, WHITE);
  1085.     for (is = 0; is < ns; is++) {
  1086.       if (DISPL_ALL == ON)
  1087.         nbr_site ((s + is)->coord.x, (s + is)->coord.y, (s + is)->nnn, imgOut, WHITE);
  1088.       else {
  1089.         if ((s + is)->eFlag != -1)
  1090.           nbr_site ((s + is)->coord.x, (s + is)->coord.y, (s + is)->nnn, imgOut, WHITE);
  1091.       }
  1092.     }
  1093.   }
  1094.   else if (CLRP == ON) {
  1095.     if (MARK_AOI == ON)
  1096.       draw_rect (AULX, AULY, ALRX, ALRY, imgOut, WHITE);
  1097.     for (is = 0; is < ns; is++) {
  1098.       if (DISPL_ALL == ON)
  1099.         //fill_poly( (s+is)->coord.x, (s+is)->coord.y, (s+is)->nnn, LABEL_MODE );            
  1100.         fill_Voronoi ((s + is)->coord.x, (s + is)->coord.y, (s + is)->nnn, imgOut, LABEL_MODE);
  1101.       else {
  1102.         if ((s + is)->eFlag != -1);
  1103.         //fill_poly( (s+is)->coord.x, (s+is)->coord.y, (s+is)->nnn, LABEL_MODE );
  1104.         fill_Voronoi ((s + is)->coord.x, (s + is)->coord.y, (s + is)->nnn, imgOut, LABEL_MODE);
  1105.       }
  1106.     }
  1107.  
  1108. /* 
  1109.  * draw border in active frame buffer 
  1110.  */
  1111.     if (BORDER == ON) {
  1112.       printf ("\n...generate border...\n");
  1113.  
  1114.       draw_border (xmin, ymin, xmax, ymax, imgOut, WHITE);
  1115.     }
  1116.   }
  1117.   else {
  1118.     if (MARK_AOI == ON)
  1119.       draw_rect (AULX, AULY, ALRX, ALRY, imgOut, WHITE);
  1120.  
  1121.   }
  1122.  
  1123.  
  1124.  
  1125. /*
  1126.  * distance statistics:
  1127.  * --> single-sites: p-function
  1128.  * -->nn site-pairs: q-function
  1129.  */
  1130.  
  1131. /* 
  1132.  * p-function: evaluate ``numerically'', employing graphics fill/flood fcts
  1133.  */
  1134.   if (P_STATS == ON) {
  1135.     printf ("\n...distance statistics: -->p-function\n");
  1136.  
  1137.     if ((pf = (struct Tuple *) calloc (ND + 1, sizeof (struct Tuple))) == NULL)
  1138.         fail_alloc ("pf", 1);
  1139.  
  1140.     r0 = sqrt ((float) NXNY / ((float) ns));
  1141.     eval_p (s, ns, pf, ND, (double) r0, (double) NXNY);
  1142.  
  1143. #ifdef SHOW_P
  1144.     printf ("\n...single-site distance statistics: p-function\n");
  1145.     for (i = 0; i <= ND; i++) {
  1146.       printf (" %6.2f", (pf + i)->y);
  1147.       if ((i + 1) % 10 == 0)
  1148.         printf ("\n");
  1149.     }
  1150. #endif
  1151.   }
  1152.  
  1153. /* 
  1154.  * q-function: evaluate by computing partial sums of nnd histogram
  1155.  */
  1156.   if (Q_STATS == ON) {
  1157.     printf ("\n...distance statistics: -->q-function\n");
  1158.  
  1159.     if ((qf = (struct Tuple *) calloc (ND + 1, sizeof (struct Tuple))) == NULL)
  1160.         fail_alloc ("qf", 1);
  1161.  
  1162.     r0 = sqrt ((float) NXNY / ((float) ns));
  1163.     eval_q (ha, qf, ND, (double) r0);
  1164.  
  1165. #ifdef SHOW_Q
  1166.     printf ("\n...nn distance statistics: q-function\n");
  1167.     for (i = 0; i <= ND; i++) {
  1168.       printf (" %6.2f", (qf + i)->y);
  1169.       if ((i + 1) % 10 == 0)
  1170.         printf ("\n");
  1171.     }
  1172. #endif
  1173.   }
  1174.  
  1175.  
  1176. /*
  1177.  * write output file of type .hdt (see also: .adt, .pdt)
  1178.  */
  1179.   if (WRITE_FILE == 1) {
  1180.  
  1181.     if ((foo = (float *) calloc (np, sizeof (float))) == NULL)
  1182.         fail_alloc ("foo", 1);
  1183.  
  1184.     for (i = 0; i < np; i++)
  1185.       *(foo + i) = (float) 1.0;
  1186.  
  1187.     if (HIST_TYPE != -1) {
  1188.       h->nh = nh;
  1189.       h->nb = N_BINS;
  1190.       h->bw = bw;
  1191.       h->amin = min;
  1192.       h->amax = max;
  1193.       h->mean = mean;
  1194.       h->sdev = sdev;
  1195.       h->phd = &hd[0];
  1196.       h->ph = &ha[0];
  1197.  
  1198.       if (P_STATS == ON)
  1199.         dstats = pf;
  1200.       if (Q_STATS == ON)
  1201.         dstats = qf;
  1202.       else
  1203.         dstats = (struct Tuple *) NULL;
  1204.  
  1205.       write_hdt_data (fpOut, np, foo, h, dstats, ND);
  1206.     }
  1207.     fclose (fpOut);
  1208.   }
  1209.  
  1210.   printf ("\n...writing output file %s\n", argv[2]);
  1211.   ImageOut (argv[2], imgOut);
  1212.  
  1213.   if (WRITE_FILE == 1)
  1214.     free (foo);
  1215.   if (HIST_TYPE != -1) {
  1216.     free (hd);
  1217.     free (ha);
  1218.   }
  1219.   free (nnna);                  /* alloc in vora.c */
  1220.   for (i = 0; i < ns; i++) {
  1221.     free ((s + i)->nnd);
  1222.     free ((s + i)->nnsi);
  1223.   }
  1224.   free (s);
  1225.   free (v);
  1226.   free (e);
  1227.   for (i = 0; i < ns; i++) {
  1228.     free ((p + i)->pei);
  1229.     free ((p + i)->vel);
  1230.     free ((p + i)->sphi);
  1231.     if (CPV == 1)
  1232.       free ((p + i)->pvi);
  1233.   }
  1234.   free (p);
  1235.   for (inn = 0; inn < nnn_max; inn++)
  1236.     free (a_n[inn]);
  1237.   free (a_n);
  1238.   free (a_nbar);
  1239.   free (npna);
  1240.  
  1241.   if (EVAL_KNN == 1) {
  1242.     for (is = 0; is <= ns; is++)
  1243.       free ((sk + is)->aka);
  1244.     for (is = 0; is <= ns; is++)
  1245.       free ((sk + is)->tctc);
  1246.     free (sk);
  1247.     free (sha);
  1248.     free (tctck);
  1249.   }
  1250.   free (pem);
  1251.   if (P_STATS == ON)
  1252.     free (pf);
  1253.   if (Q_STATS == ON)
  1254.     free (qf);
  1255. }
  1256.