home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / x / gui / x3d.lha / x3d / hull / hull.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-05  |  42.1 KB  |  1,518 lines
  1. /*
  2.   Copyright 1992 Mark Spychalla
  3.  
  4.   Permission to use, copy, modify, distribute, and sell this software and
  5.   its documentation for any purpose is hereby granted without fee,
  6.   provided that the above copyright notice appear in all copies and that
  7.   both that copyright notice and this permission notice appear in
  8.   supporting documentation, and that the name of Mark Spychalla not be used
  9.   in advertising or publicity pertaining to distribution of the software
  10.   without specific, written prior permission.  Mark Spychalla makes no
  11.   representations about the suitability of this software for any purpose.
  12.   It is provided "as is" without express or implied warranty.
  13.  
  14.   Mark Spychalla disclaims all warranties with regard to this software,
  15.   including all implied warranties of merchantability and fitness, in no
  16.   event shall Mark Spychalla be liable for any special, indirect or
  17.   consequential damages or any damages whatsoever resulting from loss of use,
  18.   data or profits, whether in an action of contract, negligence or other
  19.   tortious action, arising out of or in connection with the use or performance
  20.   of this software.
  21. */
  22.  
  23. #include <stdio.h>
  24. #include <math.h>
  25. #include <string.h>
  26.  
  27. #define NUMCOLORS      11 
  28. #define MAXLINE        8192 
  29. #define EPSILON        0.01
  30. #define TRUE           1
  31. #define FALSE          0
  32.  
  33. #define EOK            0
  34. #define SYSTEM       -1
  35.  
  36. #define    RED    2
  37. #define    GREEN    1
  38. #define BLUE    0
  39.  
  40. typedef struct POLY Poly;
  41.  
  42. typedef struct POINT{
  43.    double x, y, z;
  44.    int index;
  45.    int used;
  46. } Point;
  47.  
  48. typedef struct SEG{
  49.    Point *p, *q;
  50.    int index;
  51.    Poly *poly;
  52.    int lookupColor;
  53. } Seg;
  54.  
  55. struct POLY{
  56.    Point **points;
  57.    Seg   **segs;
  58.    int red;
  59.    int green;
  60.    int blue;
  61.    int numPoints;
  62.    int numSegs;
  63.    int index;
  64.    int used;
  65.    int lookupColor;
  66. };
  67.  
  68. typedef struct HASHENTRY{
  69.    Seg   *seg;
  70.    int used;
  71. } hashEntry;
  72.  
  73. struct box {
  74.     int r0;             /* min value, exclusive */
  75.     int r1;             /* max value, inclusive */
  76.     int g0;  
  77.     int g1;  
  78.     int b0;  
  79.     int b1;
  80.     int vol;
  81. };
  82.  
  83. /* Histogram is in elements 1..HISTSIZE along each axis,
  84.  * element 0 is for base or marginal value
  85.  * NB: these must start out 0!
  86.  */
  87.  
  88. float        m2[33][33][33];
  89. long int    wt[33][33][33], mr[33][33][33],    mg[33][33][33],    mb[33][33][33];
  90. unsigned char   *Ir, *Ig, *Ib;
  91. int            size; /*image size*/
  92. int        K;    /*color look-up table size*/
  93. unsigned int *Qadd;
  94.  
  95.  
  96. Point *dotCompareA, *dotCompareB;
  97. hashEntry *hashTable;
  98. int hashTableSize;
  99.  
  100.  
  101.  
  102. /**********************************************************************
  103.    Code for quantizer from Graphics Gems II starts here.
  104. **********************************************************************/
  105.  
  106.  
  107.  
  108. /**********************************************************************
  109.         C Implementation of Wu's Color Quantizer (v. 2)
  110.         (see Graphics Gems vol. II, pp. 126-133)
  111.  
  112. Author:    Xiaolin Wu
  113.     Dept. of Computer Science
  114.     Univ. of Western Ontario
  115.     London, Ontario N6A 5B7
  116.     wu@csd.uwo.ca
  117.  
  118. Algorithm: Greedy orthogonal bipartition of RGB space for variance
  119.        minimization aided by inclusion-exclusion tricks.
  120.        For speed no nearest neighbor search is done. Slightly
  121.        better performance can be expected by more sophisticated
  122.        but more expensive versions.
  123.  
  124. The author thanks Tom Lane at Tom_Lane@G.GP.CS.CMU.EDU for much of
  125. additional documentation and a cure to a previous bug.
  126.  
  127. Free to distribute, comments and suggestions are appreciated.
  128. **********************************************************************/    
  129.  
  130.  
  131. static void
  132. Hist3d(vwt, vmr, vmg, vmb, M2) 
  133. /* build 3-D color histogram of counts, r/g/b, c^2 */
  134. long int *vwt, *vmr, *vmg, *vmb;
  135. float *M2;
  136. {
  137. register int ind, r, g, b;
  138. int         inr, ing, inb, table[256];
  139. register long int i;
  140.         
  141.     for(i=0; i<256; ++i) table[i]=i*i;
  142.     Qadd = (unsigned int *)malloc(sizeof(int)*size);
  143.     if (Qadd==NULL) {fprintf(stderr, "Not enough space\n"); exit(1);}
  144.     for(i=0; i<size; ++i){
  145.         r = Ir[i]; g = Ig[i]; b = Ib[i];
  146.         inr=(r>>3)+1; 
  147.         ing=(g>>3)+1; 
  148.         inb=(b>>3)+1; 
  149.         Qadd[i]=ind=(inr<<10)+(inr<<6)+inr+(ing<<5)+ing+inb;
  150.         /*[inr][ing][inb]*/
  151.         ++vwt[ind];
  152.         vmr[ind] += r;
  153.         vmg[ind] += g;
  154.         vmb[ind] += b;
  155.              M2[ind] += (float)(table[r]+table[g]+table[b]);
  156.     }
  157. }
  158.  
  159. /* At conclusion of the histogram step, we can interpret
  160.  *   wt[r][g][b] = sum over voxel of P(c)
  161.  *   mr[r][g][b] = sum over voxel of r*P(c)  ,  similarly for mg, mb
  162.  *   m2[r][g][b] = sum over voxel of c^2*P(c)
  163.  * Actually each of these should be divided by 'size' to give the usual
  164.  * interpretation of P() as ranging from 0 to 1, but we needn't do that here.
  165.  */
  166.  
  167. /* We now convert histogram into moments so that we can rapidly calculate
  168.  * the sums of the above quantities over any desired box.
  169.  */
  170.  
  171.  
  172. static void
  173. M3d(vwt, vmr, vmg, vmb, M2) /* compute cumulative moments. */
  174. long int *vwt, *vmr, *vmg, *vmb;
  175. float    *M2;
  176. {
  177. register unsigned int ind1, ind2;
  178. register unsigned char i, r, g, b;
  179. long int line, line_r, line_g, line_b,
  180.      area[33], area_r[33], area_g[33], area_b[33];
  181. float    line2, area2[33];
  182.  
  183.     for(r=1; r<=32; ++r){
  184.     for(i=0; i<=32; ++i) 
  185.         area2[i]=area[i]=area_r[i]=area_g[i]=area_b[i]=0;
  186.     for(g=1; g<=32; ++g){
  187.         line2 = line = line_r = line_g = line_b = 0;
  188.         for(b=1; b<=32; ++b){
  189.         ind1 = (r<<10) + (r<<6) + r + (g<<5) + g + b; /* [r][g][b] */
  190.         line += vwt[ind1];
  191.         line_r += vmr[ind1]; 
  192.         line_g += vmg[ind1]; 
  193.         line_b += vmb[ind1];
  194.         line2 += M2[ind1];
  195.         area[b] += line;
  196.         area_r[b] += line_r;
  197.         area_g[b] += line_g;
  198.         area_b[b] += line_b;
  199.         area2[b] += line2;
  200.         ind2 = ind1 - 1089; /* [r-1][g][b] */
  201.         vwt[ind1] = vwt[ind2] + area[b];
  202.         vmr[ind1] = vmr[ind2] + area_r[b];
  203.         vmg[ind1] = vmg[ind2] + area_g[b];
  204.         vmb[ind1] = vmb[ind2] + area_b[b];
  205.         M2[ind1] = M2[ind2] + area2[b];
  206.         }
  207.     }
  208.     }
  209. }
  210.  
  211.  
  212. static long int Vol(cube, mmt) 
  213. /* Compute sum over a box of any given statistic */
  214. struct box *cube;
  215. long int mmt[33][33][33];
  216. {
  217.     return( mmt[cube->r1][cube->g1][cube->b1] 
  218.        -mmt[cube->r1][cube->g1][cube->b0]
  219.        -mmt[cube->r1][cube->g0][cube->b1]
  220.        +mmt[cube->r1][cube->g0][cube->b0]
  221.        -mmt[cube->r0][cube->g1][cube->b1]
  222.        +mmt[cube->r0][cube->g1][cube->b0]
  223.        +mmt[cube->r0][cube->g0][cube->b1]
  224.        -mmt[cube->r0][cube->g0][cube->b0] );
  225. }
  226.  
  227. /* The next two routines allow a slightly more efficient calculation
  228.  * of Vol() for a proposed subbox of a given box.  The sum of Top()
  229.  * and Bottom() is the Vol() of a subbox split in the given direction
  230.  * and with the specified new upper bound.
  231.  */
  232.  
  233. static long int Bottom(cube, dir, mmt)
  234. /* Compute part of Vol(cube, mmt) that doesn't depend on r1, g1, or b1 */
  235. /* (depending on dir) */
  236. struct box *cube;
  237. unsigned char dir;
  238. long int mmt[33][33][33];
  239. {
  240.     switch(dir){
  241.     case RED:
  242.         return( -mmt[cube->r0][cube->g1][cube->b1]
  243.             +mmt[cube->r0][cube->g1][cube->b0]
  244.             +mmt[cube->r0][cube->g0][cube->b1]
  245.             -mmt[cube->r0][cube->g0][cube->b0] );
  246.         break;
  247.     case GREEN:
  248.         return( -mmt[cube->r1][cube->g0][cube->b1]
  249.             +mmt[cube->r1][cube->g0][cube->b0]
  250.             +mmt[cube->r0][cube->g0][cube->b1]
  251.             -mmt[cube->r0][cube->g0][cube->b0] );
  252.         break;
  253.     case BLUE:
  254.         return( -mmt[cube->r1][cube->g1][cube->b0]
  255.             +mmt[cube->r1][cube->g0][cube->b0]
  256.             +mmt[cube->r0][cube->g1][cube->b0]
  257.             -mmt[cube->r0][cube->g0][cube->b0] );
  258.         break;
  259.     }
  260.    return 0;
  261. }
  262.  
  263.  
  264. static long int Top(cube, dir, pos, mmt)
  265. /* Compute remainder of Vol(cube, mmt), substituting pos for */
  266. /* r1, g1, or b1 (depending on dir) */
  267. struct box *cube;
  268. unsigned char dir;
  269. int   pos;
  270. long int mmt[33][33][33];
  271. {
  272.     switch(dir){
  273.     case RED:
  274.         return( mmt[pos][cube->g1][cube->b1] 
  275.            -mmt[pos][cube->g1][cube->b0]
  276.            -mmt[pos][cube->g0][cube->b1]
  277.            +mmt[pos][cube->g0][cube->b0] );
  278.         break;
  279.     case GREEN:
  280.         return( mmt[cube->r1][pos][cube->b1] 
  281.            -mmt[cube->r1][pos][cube->b0]
  282.            -mmt[cube->r0][pos][cube->b1]
  283.            +mmt[cube->r0][pos][cube->b0] );
  284.         break;
  285.     case BLUE:
  286.         return( mmt[cube->r1][cube->g1][pos]
  287.            -mmt[cube->r1][cube->g0][pos]
  288.            -mmt[cube->r0][cube->g1][pos]
  289.            +mmt[cube->r0][cube->g0][pos] );
  290.         break;
  291.     }
  292.    return 0;
  293. }
  294.  
  295.  
  296. static float Var(cube)
  297. /* Compute the weighted variance of a box */
  298. /* NB: as with the raw statistics, this is really the variance * size */
  299. struct box *cube;
  300. {
  301. float dr, dg, db, xx;
  302.  
  303.     dr = Vol(cube, mr); 
  304.     dg = Vol(cube, mg); 
  305.     db = Vol(cube, mb);
  306.     xx =  m2[cube->r1][cube->g1][cube->b1] 
  307.      -m2[cube->r1][cube->g1][cube->b0]
  308.      -m2[cube->r1][cube->g0][cube->b1]
  309.      +m2[cube->r1][cube->g0][cube->b0]
  310.      -m2[cube->r0][cube->g1][cube->b1]
  311.      +m2[cube->r0][cube->g1][cube->b0]
  312.      +m2[cube->r0][cube->g0][cube->b1]
  313.      -m2[cube->r0][cube->g0][cube->b0];
  314.  
  315.     return( xx - (dr*dr+dg*dg+db*db)/(float)Vol(cube,wt) );    
  316. }
  317.  
  318. /* We want to minimize the sum of the variances of two subboxes.
  319.  * The sum(c^2) terms can be ignored since their sum over both subboxes
  320.  * is the same (the sum for the whole box) no matter where we split.
  321.  * The remaining terms have a minus sign in the variance formula,
  322.  * so we drop the minus sign and MAXIMIZE the sum of the two terms.
  323.  */
  324.  
  325.  
  326. static float Maximize(cube, dir, first, last, cut,
  327.         whole_r, whole_g, whole_b, whole_w)
  328. struct box *cube;
  329. unsigned char dir;
  330. int first, last, *cut;
  331. long int whole_r, whole_g, whole_b, whole_w;
  332. {
  333. register long int half_r, half_g, half_b, half_w;
  334. long int base_r, base_g, base_b, base_w;
  335. register int i;
  336. register float temp, max;
  337.  
  338.     base_r = Bottom(cube, dir, mr);
  339.     base_g = Bottom(cube, dir, mg);
  340.     base_b = Bottom(cube, dir, mb);
  341.     base_w = Bottom(cube, dir, wt);
  342.     max = 0.0;
  343.     *cut = -1;
  344.     for(i=first; i<last; ++i){
  345.     half_r = base_r + Top(cube, dir, i, mr);
  346.     half_g = base_g + Top(cube, dir, i, mg);
  347.     half_b = base_b + Top(cube, dir, i, mb);
  348.     half_w = base_w + Top(cube, dir, i, wt);
  349.         /* now half_x is sum over lower half of box, if split at i */
  350.         if (half_w == 0) {      /* subbox could be empty of pixels! */
  351.           continue;             /* never split into an empty box */
  352.     } else
  353.         temp = ((float)half_r*half_r + (float)half_g*half_g +
  354.                 (float)half_b*half_b)/half_w;
  355.  
  356.     half_r = whole_r - half_r;
  357.     half_g = whole_g - half_g;
  358.     half_b = whole_b - half_b;
  359.     half_w = whole_w - half_w;
  360.         if (half_w == 0) {      /* subbox could be empty of pixels! */
  361.           continue;             /* never split into an empty box */
  362.     } else
  363.         temp += ((float)half_r*half_r + (float)half_g*half_g +
  364.                  (float)half_b*half_b)/half_w;
  365.  
  366.         if (temp > max) {max=temp; *cut=i;}
  367.     }
  368.     return(max);
  369. }
  370.  
  371. static int
  372. Cut(set1, set2)
  373. struct box *set1, *set2;
  374. {
  375. unsigned char dir;
  376. int cutr, cutg, cutb;
  377. float maxr, maxg, maxb;
  378. long int whole_r, whole_g, whole_b, whole_w;
  379.  
  380.     whole_r = Vol(set1, mr);
  381.     whole_g = Vol(set1, mg);
  382.     whole_b = Vol(set1, mb);
  383.     whole_w = Vol(set1, wt);
  384.  
  385.     maxr = Maximize(set1, RED, set1->r0+1, set1->r1, &cutr,
  386.             whole_r, whole_g, whole_b, whole_w);
  387.     maxg = Maximize(set1, GREEN, set1->g0+1, set1->g1, &cutg,
  388.             whole_r, whole_g, whole_b, whole_w);
  389.     maxb = Maximize(set1, BLUE, set1->b0+1, set1->b1, &cutb,
  390.             whole_r, whole_g, whole_b, whole_w);
  391.  
  392.     if( (maxr>=maxg)&&(maxr>=maxb) ) {
  393.     dir = RED;
  394.     if (cutr < 0) return 0; /* can't split the box */
  395.     }
  396.     else
  397.     if( (maxg>=maxr)&&(maxg>=maxb) ) 
  398.     dir = GREEN;
  399.     else
  400.     dir = BLUE; 
  401.  
  402.     set2->r1 = set1->r1;
  403.     set2->g1 = set1->g1;
  404.     set2->b1 = set1->b1;
  405.  
  406.     switch (dir){
  407.     case RED:
  408.         set2->r0 = set1->r1 = cutr;
  409.         set2->g0 = set1->g0;
  410.         set2->b0 = set1->b0;
  411.         break;
  412.     case GREEN:
  413.         set2->g0 = set1->g1 = cutg;
  414.         set2->r0 = set1->r0;
  415.         set2->b0 = set1->b0;
  416.         break;
  417.     case BLUE:
  418.         set2->b0 = set1->b1 = cutb;
  419.         set2->r0 = set1->r0;
  420.         set2->g0 = set1->g0;
  421.         break;
  422.     }
  423.     set1->vol=(set1->r1-set1->r0)*(set1->g1-set1->g0)*(set1->b1-set1->b0);
  424.     set2->vol=(set2->r1-set2->r0)*(set2->g1-set2->g0)*(set2->b1-set2->b0);
  425.     return 1;
  426. }
  427.  
  428.  
  429. static void Mark(cube, label, tag)
  430. struct box *cube;
  431. int label;
  432. unsigned char *tag;
  433. {
  434. register int r, g, b;
  435.  
  436.     for(r=cube->r0+1; r<=cube->r1; ++r)
  437.        for(g=cube->g0+1; g<=cube->g1; ++g)
  438.       for(b=cube->b0+1; b<=cube->b1; ++b)
  439.         tag[(r<<10) + (r<<6) + r + (g<<5) + g + b] = label;
  440. }
  441.  
  442. static void quantize(red, green, blue, numColors, quantColors, lookup)
  443. unsigned char *red, *green, *blue;
  444. int numColors, *quantColors;
  445. unsigned *lookup;
  446. {
  447. struct box    *cube;
  448. unsigned char    *tag;
  449. unsigned char    *lut_r, *lut_g, *lut_b;
  450. int        next;
  451. register long int    i, weight;
  452. register int    k;
  453. float        *vv, temp;
  454.  
  455.         if((cube = (struct box *)malloc(*quantColors * sizeof(struct box)))
  456.         == NULL){
  457.            fprintf(stderr, "Unable to allocate memory for cube\n");
  458.            exit(1);
  459.            }
  460.  
  461.         if((lut_r = (unsigned char *)malloc(*quantColors *
  462.         sizeof(unsigned char))) == NULL){
  463.            fprintf(stderr, "Unable to allocate memory for lut_r\n");
  464.            exit(1);
  465.            }
  466.  
  467.         if((lut_g = (unsigned char *)malloc(*quantColors *
  468.         sizeof(unsigned char))) == NULL){
  469.            fprintf(stderr, "Unable to allocate memory for lut_g\n");
  470.            exit(1);
  471.            }
  472.  
  473.         if((lut_b = (unsigned char *)malloc(*quantColors *
  474.         sizeof(unsigned char))) == NULL){
  475.            fprintf(stderr, "Unable to allocate memory for lut_b\n");
  476.            exit(1);
  477.            }
  478.  
  479.         if((vv = (float *)malloc(*quantColors *
  480.         sizeof(float))) == NULL){
  481.            fprintf(stderr, "Unable to allocate memory for lut_b\n");
  482.            exit(1);
  483.            }
  484.  
  485.     /* input R,G,B components into Ir, Ig, Ib;
  486.        set size to width*height */
  487.  
  488.         Ir = red;
  489.         Ig = green;
  490.         Ib = blue;
  491.       size = numColors;
  492.          K = *quantColors;
  493.  
  494.         if((Ig = (unsigned char *)malloc(numColors * sizeof(unsigned char)))
  495.         == NULL){
  496.            fprintf(stderr, "Unable to allocate memory for Ig\n");
  497.            exit(1);
  498.            }
  499.  
  500.         if((Ib = (unsigned char *)malloc(numColors * sizeof(unsigned char)))
  501.         == NULL){
  502.            fprintf(stderr, "Unable to allocate memory for Ig\n");
  503.            exit(1);
  504.            }
  505.  
  506.         if((Ir = (unsigned char *)malloc(numColors * sizeof(unsigned char)))
  507.         == NULL){
  508.            fprintf(stderr, "Unable to allocate memory for Ig\n");
  509.            exit(1);
  510.            }
  511.  
  512.         memcpy(Ir, red,   numColors * sizeof(unsigned char));
  513.         memcpy(Ig, green, numColors * sizeof(unsigned char));
  514.         memcpy(Ib, blue,  numColors * sizeof(unsigned char));
  515.  
  516. /*
  517.     printf("no. of colors:\n");
  518.     scanf("%d", &K);
  519. */
  520.  
  521.     Hist3d(wt, mr, mg, mb, m2); fprintf(stderr, "Histogram done\n");
  522.  
  523.     free(Ig); free(Ib); free(Ir);
  524.  
  525.     M3d(wt, mr, mg, mb, m2);    fprintf(stderr, "Moments done\n");
  526.  
  527.     cube[0].r0 = cube[0].g0 = cube[0].b0 = 0;
  528.     cube[0].r1 = cube[0].g1 = cube[0].b1 = 32;
  529.     next = 0;
  530.         for(i=1; i<K; ++i){
  531.             if (Cut(&cube[next], &cube[i])) {
  532.               /* volume test ensures we won't try to cut one-cell box */
  533.               vv[next] = (cube[next].vol>1) ? Var(&cube[next]) : 0.0;
  534.               vv[i] = (cube[i].vol>1) ? Var(&cube[i]) : 0.0;
  535.         } else {
  536.               vv[next] = 0.0;   /* don't try to split this box again */
  537.               i--;              /* didn't create box i */
  538.         }
  539.             next = 0; temp = vv[0];
  540.             for(k=1; k<=i; ++k)
  541.                 if (vv[k] > temp) {
  542.                     temp = vv[k]; next = k;
  543.         }
  544.             if (temp <= 0.0) {
  545.               K = i+1;
  546.               fprintf(stderr, "Only got %d boxes\n", K);
  547.               break;
  548.         }
  549.     }
  550.         fprintf(stderr, "Partition done\n");
  551.  
  552.     /* the space for array m2 can be freed now */
  553.  
  554.     tag = (unsigned char *)malloc(33*33*33);
  555.     if (tag==NULL) {fprintf(stderr, "Not enough space\n"); exit(1);}
  556.     for(k=0; k<K; ++k){
  557.         Mark(&cube[k], k, tag);
  558.         weight = Vol(&cube[k], wt);
  559.         if (weight) {
  560.         lut_r[k] = Vol(&cube[k], mr) / weight;
  561.         lut_g[k] = Vol(&cube[k], mg) / weight;
  562.         lut_b[k] = Vol(&cube[k], mb) / weight;
  563.         }
  564.         else{
  565.           fprintf(stderr, "bogus box %d\n", k);
  566.           lut_r[k] = lut_g[k] = lut_b[k] = 0;        
  567.         }
  568.     }
  569.  
  570.     for(i=0; i<size; ++i) Qadd[i] = tag[Qadd[i]];
  571.     
  572.     /* output lut_r, lut_g, lut_b as color look-up table contents,
  573.        Qadd as the quantized image (array of table addresses). */
  574.  
  575.         free(tag);
  576.         memcpy(red,   lut_r, *quantColors);
  577.         memcpy(green, lut_g, *quantColors);
  578.         memcpy(blue,  lut_b, *quantColors);
  579.         memcpy(lookup, Qadd, numColors * sizeof(int));
  580.         *quantColors = K;
  581.         free(cube);
  582.         free(lut_r);
  583.         free(lut_g);
  584.         free(lut_b);
  585.         free(vv);
  586. }
  587.  
  588.  
  589.  
  590.  
  591. /******************************************************************************
  592.    My own hull.c code starts here. 
  593. ******************************************************************************/
  594.  
  595.  
  596.  
  597.  
  598. static int fequal(x,y)
  599. double x, y;
  600. /******************************************************************************/
  601. /* Are x and y scaled equal within epsilon?                                   */
  602. /******************************************************************************/
  603. {
  604.    return (((fabs((x) - (y))) / (fabs(x) + 1.0)) < EPSILON);
  605. }
  606.  
  607.  
  608.  
  609. static int pointEqual(p,q)
  610. Point *p, *q;
  611. /******************************************************************************/
  612. /* Are p and q equal?                                                         */
  613. /******************************************************************************/
  614. {
  615.    return(fequal(p->x,q->x) && fequal(p->y,q->y) && fequal(p->z,q->z));
  616. }
  617.  
  618.  
  619.  
  620. static void hash(p,q,segs,numSegs,seg,maxSegs,poly)
  621. Point *p, *q;
  622. Seg   *segs;
  623. int   *numSegs;
  624. Seg   **seg;
  625. int   maxSegs;
  626. Poly  *poly;
  627. /******************************************************************************/
  628. /* Compute hash function for segment                                          */
  629. /******************************************************************************/
  630. {
  631. int value, index;
  632.  
  633.    value = ((int)(p->x * q->x) + (int)(p->y * q->y) +
  634.             (int)(p->z * q->z));
  635.  
  636.    if(value < 0){
  637.       value = (-value) % hashTableSize;
  638.    }else{
  639.       value = value % hashTableSize;
  640.       }
  641.  
  642.    index = value;
  643.  
  644.    while((hashTable[index].used == TRUE) && (
  645.      !((p->x == hashTable[index].seg->p->x) &&
  646.        (p->y == hashTable[index].seg->p->y) &&
  647.        (p->z == hashTable[index].seg->p->z) &&
  648.        (q->x == hashTable[index].seg->q->x) &&
  649.        (q->y == hashTable[index].seg->q->y) &&
  650.        (q->z == hashTable[index].seg->q->z))   )&&(
  651.      !((q->x == hashTable[index].seg->p->x) &&
  652.        (q->y == hashTable[index].seg->p->y) &&
  653.        (q->z == hashTable[index].seg->p->z) &&
  654.        (p->x == hashTable[index].seg->q->x) &&
  655.        (p->y == hashTable[index].seg->q->y) &&
  656.        (p->z == hashTable[index].seg->q->z)))) {
  657.           index = (index + 1) % hashTableSize; 
  658.           }
  659.   
  660.    if(hashTable[index].used == TRUE){
  661.       *seg = hashTable[index].seg;
  662.    }else{
  663.  
  664.       if(*numSegs == maxSegs){
  665.          fprintf(stderr, "Too many Segments\n");
  666.          exit(1);
  667.          }
  668.  
  669.       hashTable[index].seg  = &(segs[*numSegs]);
  670.       segs[*numSegs].p = p;
  671.       segs[*numSegs].q = q;
  672.       segs[*numSegs].poly = poly;
  673.       p->used = TRUE;
  674.       q->used = TRUE;
  675.       segs[*numSegs].index = *numSegs;
  676.       (*numSegs)++;
  677.  
  678.       hashTable[index].used = TRUE;
  679.       *seg = hashTable[index].seg;
  680.       } 
  681. }
  682.  
  683.  
  684.  
  685. static int comparePoints(p,q)
  686. Point *p, *q;
  687. /******************************************************************************/
  688. /* Compare p and q.                                                           */
  689. /******************************************************************************/
  690. {
  691. /* Sorted by X, Y, then Z */
  692.  
  693. /* Lower X comes before higher X */
  694.  
  695.    if(p->x < q->x)
  696.       return -1;
  697.  
  698.    if(p->x > q->x)
  699.       return 1;
  700.  
  701. /* Lower Y comes before higher Y */
  702.  
  703.    if(p->y < q->y)
  704.       return -1;
  705.  
  706.    if(p->y > q->y)
  707.       return 1;
  708.  
  709. /* Lower Z comes before higher Z */
  710.  
  711.    if(p->z < q->z)
  712.       return -1;
  713.  
  714.    if(p->z > q->z)
  715.       return 1;
  716.  
  717.    return 0;
  718. }
  719.  
  720.  
  721.  
  722. static int pComparePoints(p,q)
  723. Point **p, **q;
  724. /******************************************************************************/
  725. /* Compare p and q.                                                           */
  726. /******************************************************************************/
  727. {
  728.    return comparePoints(*p,*q);
  729. }
  730.  
  731.  
  732.  
  733. static int comparePolys(p,q)
  734. Poly *p, *q;
  735. /******************************************************************************/
  736. /* Compare p and q.                                                           */
  737. /******************************************************************************/
  738. {
  739. int index;
  740.  
  741.    if(p->numPoints < q->numPoints){
  742.       return -1;
  743.       }
  744.    
  745.    if(p->numPoints > q->numPoints){
  746.       return 1;
  747.       }
  748.  
  749.    index = 0;
  750.  
  751.    while(index < p->numPoints){
  752.  
  753.       if(p->points[index] < q->points[index])
  754.          return -1;
  755.    
  756.       if(p->points[index] > q->points[index])
  757.          return 1;
  758.    
  759.       index++;
  760.       }
  761.  
  762.    return 0;
  763. }
  764.  
  765.  
  766.  
  767. static int colorComparePolys(p,q)
  768. Poly *p, *q;
  769. /******************************************************************************/
  770. /* Compare p and q by color.                                                  */
  771. /******************************************************************************/
  772. {
  773.    if(p->lookupColor < q->lookupColor){
  774.       return -1;
  775.       }
  776.  
  777.    if(p->lookupColor > q->lookupColor){
  778.       return 1;
  779.       }
  780.  
  781.    return 0;
  782. }
  783.  
  784.  
  785.  
  786. static int colorCompareSegs(p,q)
  787. Seg *p, *q;
  788. /******************************************************************************/
  789. /* Compare p and q by color.                                                  */
  790. /******************************************************************************/
  791. {
  792.    if(p->lookupColor < q->lookupColor){
  793.       return -1;
  794.       }
  795.  
  796.    if(p->lookupColor > q->lookupColor){
  797.       return 1;
  798.       }
  799.  
  800.    return 0;
  801. }
  802.  
  803.  
  804.  
  805. static double dot(a, b, c)
  806. Point *a, *b, *c;
  807. /******************************************************************************/
  808. /* Compute  ab (dot) ac.                                                      */
  809. /******************************************************************************/
  810. {
  811. double X1, Y1, Z1, X2, Y2, Z2, len1, len2;
  812.  
  813.    if((comparePoints(a, b) == 0) || (comparePoints(a, c) == 0)){
  814.       return 2.0;
  815.       }
  816.  
  817.    X1 = b->x - a->x; Y1 = b->y - a->y; Z1 = b->z - a->z;
  818.    X2 = c->x - a->x; Y2 = c->y - a->y; Z2 = c->z - a->z;
  819.  
  820.    len1 = sqrt(X1 * X1 + Y1 * Y1 + Z1 * Z1);
  821.    len2 = sqrt(X2 * X2 + Y2 * Y2 + Z2 * Z2);
  822.  
  823.    X1 /= len1; Y1 /= len1; Z1 /= len1;
  824.    X2 /= len2; Y2 /= len2; Z2 /= len2;
  825.  
  826.    return(X1 * X2 + Y1 * Y2 + Z1 * Z2);
  827. }
  828.  
  829.  
  830.  
  831. static int dotComparePoints(p,q)
  832. Point **p, **q;
  833. /******************************************************************************/
  834. /* Compare p and q.                                                           */
  835. /******************************************************************************/
  836. {
  837. double dotP, dotQ;
  838.  
  839.    dotP = dot(dotCompareA, dotCompareB, *p);
  840.    dotQ = dot(dotCompareA, dotCompareB, *q);
  841.  
  842.    if(dotP < dotQ)
  843.       return -1;
  844.  
  845.    if(dotP > dotQ)
  846.       return 1;
  847.  
  848.    return 0;
  849. }
  850.  
  851.  
  852.  
  853. static int readPoints(filename, points, numPoints, polys, numPolys, maxPoints,
  854. maxPolys)
  855. char  *filename;
  856. Point **points;
  857. int   *numPoints; 
  858. Poly  **polys;
  859. int   *numPolys; 
  860. int   *maxPoints;
  861. int   *maxPolys;
  862. /******************************************************************************/
  863. /* Read in a list of polygons from a file into the points, lines, and poly    */
  864. /* arrays.                                                                    */
  865. /******************************************************************************/
  866. {
  867. float X, Y, Z;
  868. double x, y, z;
  869. int count = 0, lineNo = 0, needColor = TRUE;
  870. int index, index2;
  871. char line[MAXLINE];
  872. FILE *fd;
  873.  
  874.    *numPoints = 0;
  875.    *numPolys  = 0;
  876.  
  877. /* Try to open the file */
  878.  
  879.    if(strcmp(filename, "-") == 0){
  880.       fd = stdin;
  881.    }else{
  882.       if((fd = fopen(filename, "r")) == NULL){
  883.          fprintf(stderr, "Unable to open file %s\n", filename);
  884.          exit(1);
  885.          }
  886.       }
  887.  
  888.    while(fgets(line, MAXLINE, fd) != NULL){
  889.  
  890. /* Warn about long lines */
  891.  
  892.       if(strlen(line) > (MAXLINE - 5)){
  893.          fprintf(stderr, "Line too long on line %d, line may be trucated.\n",
  894.          lineNo);
  895.          exit(1);
  896.          }
  897.  
  898. /* ignore comments */
  899.  
  900.       lineNo++;
  901.  
  902.       count = sscanf(line, "%f %f %f\n", &X, &Y, &Z);
  903.  
  904.       x = (double)X; y = (double)Y; z = (double)Z;
  905.  
  906.       if((strcmp(line, "\n") == 0) || (count == 0)){
  907.  
  908.          if((*polys)[*numPolys].numPoints != 0){
  909.             (*numPolys)++;
  910.             if(*numPolys == *maxPolys){
  911.                if((*polys = (Poly *)realloc(*polys, ((*numPolys) + 100) *
  912.                sizeof(Poly))) == NULL){
  913.                   fprintf(stderr,
  914.                   "Unable to allocate memory for polygons on line %d\n",
  915.                   lineNo);
  916.                   exit(1);
  917.                   }
  918.                bzero(&((*polys)[*maxPolys]), 100 * sizeof(Poly));
  919.                (*maxPolys) += 100;
  920.                }
  921.             needColor = TRUE;
  922.             }
  923.       }else{ 
  924.  
  925.          if(count != 3){
  926.             fprintf(stderr, "Mismatch on line %d, %d items matched\n",
  927.             lineNo, count);
  928.             exit(1);
  929.             }
  930.      
  931.          if(needColor == FALSE){
  932.  
  933.             if(*numPoints == *maxPoints){
  934.                if((*points = (Point *)realloc(*points, (*maxPoints + 100)
  935.                * sizeof(Point))) == NULL){
  936.                   fprintf(stderr,
  937.                   "Unable to allocate memory for points on line %d\n", lineNo);
  938.                   exit(1);
  939.                   }
  940.                bzero(&((*points)[*maxPoints]), 100 * sizeof(Point));
  941.                (*maxPoints) += 100;
  942.                }
  943.  
  944.             if((*polys)[*numPolys].numPoints == 0){
  945.                if(((*polys)[*numPolys].points = (Point **)malloc(
  946.                sizeof(Point *))) == NULL){
  947.                   fprintf(stderr,
  948.                   "Unable to allocate memory for polygon points on line %d\n",
  949.                   lineNo);
  950.                   exit(1);
  951.                   }
  952.             }else{
  953.                if((((*polys)[*numPolys].points = (Point **)realloc(
  954.               (*polys)[*numPolys].points, ((*polys)[*numPolys].numPoints
  955.               + 1) * sizeof(Point *)))) == NULL){
  956.                   fprintf(stderr,
  957.                   "Unable to allocate memory for polygon points on line %d\n",
  958.                   lineNo);
  959.                   exit(1);
  960.                   }
  961.                }
  962.              
  963. /* BEWARE: Evil cast just below!!! */
  964.   
  965.             (*polys)[*numPolys].points[(*polys)[*numPolys].numPoints] = 
  966.             (Point *)(&((*points)[*numPoints]) - (*points));
  967.  
  968.             ((*polys)[*numPolys].numPoints)++;
  969.  
  970.             (*points)[*numPoints].x = x; 
  971.             (*points)[*numPoints].y = y; 
  972.             (*points)[*numPoints].z = z;
  973.  
  974.             (*numPoints)++;
  975.  
  976.          }else{
  977.  
  978.             if((x < 0.0) || (x > 255.0) ||
  979.             (y < 0.0) || (y > 255.0) ||
  980.             (z < 0.0) || (z > 255.0)){
  981.                fprintf(stderr,
  982.                "Color out of range ( %f %f %f ) on line %d\n",
  983.                 x, y, z, lineNo);
  984.                exit(1);
  985.                }
  986.  
  987.             (*polys)[*numPolys].red   = (int)x;
  988.             (*polys)[*numPolys].green = (int)y;
  989.             (*polys)[*numPolys].blue  = (int)z;
  990.             needColor = FALSE;
  991.             }
  992.          }
  993.       }
  994.  
  995.    if(feof(fd) == FALSE)
  996.       return SYSTEM;
  997.    
  998.    if((*polys)[*numPolys].numPoints != 0){
  999.       (*numPolys)++;
  1000.       }
  1001.     
  1002.    fclose(fd);
  1003.  
  1004. /* Undo the Evilness of the Evil cast above */
  1005.  
  1006.    for(index = 0; index < *numPolys; index++){
  1007.       for(index2 = 0; index2 < (*polys)[index].numPoints; index2++){
  1008.          (*polys)[index].points[index2] = 
  1009.          &((*points)[(int)(*polys)[index].points[index2]]);
  1010.          }
  1011.       }
  1012.  
  1013.    return(EOK);
  1014. }
  1015.  
  1016.  
  1017.  
  1018. static void eliminateDuplicatePoints(points, numPoints, polys, numPolys,
  1019. maxPoints)
  1020. Point **points;
  1021. int   *numPoints;
  1022. Poly  *polys;
  1023. int   numPolys;
  1024. int   *maxPoints;
  1025. /******************************************************************************/
  1026. /* Eliminate duplicate points.                                                */
  1027. /******************************************************************************/
  1028. {
  1029. Point *pointsCopy;
  1030. Point *outPoints;
  1031. int index, numOutPoints, index2, oldIndex;
  1032.  
  1033. /* Sort the point list by coordinates */
  1034.  
  1035.    for(index = 0; index < *numPoints; index++){
  1036.       (*points)[index].index = index;
  1037.       }
  1038.  
  1039.    if((pointsCopy = (Point *)malloc(*maxPoints * sizeof(Point))) == NULL){
  1040.       fprintf(stderr, "Unable to allocate memory for points\n");
  1041.       exit(1);
  1042.       }
  1043.  
  1044.    if((outPoints = (Point *)malloc(*maxPoints * sizeof(Point))) == NULL){
  1045.       fprintf(stderr, "Unable to allocate memory for points\n");
  1046.       exit(1);
  1047.       }
  1048.  
  1049.    memcpy(pointsCopy, *points, *numPoints * sizeof(Point));
  1050.    qsort(pointsCopy, *numPoints, sizeof(Point), comparePoints);
  1051.  
  1052. /* Eliminate duplicate points */
  1053.  
  1054.    index = 0;
  1055.    numOutPoints = 0;
  1056.  
  1057.    oldIndex = pointsCopy[index].index;
  1058.    outPoints[numOutPoints] = pointsCopy[index];
  1059.    pointsCopy[index].index = numOutPoints;
  1060.    (*points)[oldIndex].index = pointsCopy[index].index;
  1061.    numOutPoints++;
  1062.  
  1063.    for(index = 1; index < *numPoints; index++){
  1064.       oldIndex = pointsCopy[index].index;
  1065.       if(pointEqual(&(pointsCopy[index]), &(pointsCopy[index - 1]))){
  1066.          pointsCopy[index].index = pointsCopy[index - 1].index;
  1067.       }else{
  1068.          outPoints[numOutPoints] = pointsCopy[index];
  1069.          pointsCopy[index].index = numOutPoints;
  1070.          numOutPoints++;
  1071.          }
  1072.       (*points)[oldIndex].index = pointsCopy[index].index;
  1073.       }
  1074.  
  1075.    if((outPoints = (Point *)realloc(outPoints, numOutPoints * sizeof(Point)))
  1076.    == NULL){
  1077.       fprintf(stderr, "Unable to allocate memory for points\n");
  1078.       exit(1);
  1079.       }
  1080.  
  1081.    for(index = 0; index < numPolys; index++){
  1082.  
  1083.       if(polys[index].numPoints < 3){
  1084.          fprintf(stderr, "Too few points (%d) in polygon %d\n", 
  1085.          polys[index].numPoints, index);
  1086.          exit(1);
  1087.          }
  1088.      
  1089.       for(index2 = 0; index2 < polys[index].numPoints; index2++){
  1090.          polys[index].points[index2] =
  1091.          &(outPoints[polys[index].points[index2]->index]);
  1092.          }
  1093.       }
  1094.  
  1095.    free(pointsCopy);
  1096.    free(*points);
  1097.  
  1098.    *numPoints = numOutPoints;
  1099.    *points = outPoints;
  1100. }
  1101.  
  1102.  
  1103.  
  1104. static int formPolygonsSegments(points, numPoints, polys, numPolys, segs,
  1105. numSegs, numPointsUsed, numPolysUsed, maxSegs, dupPolys, invalidPolys)
  1106. Point *points;
  1107. int   *numPoints;
  1108. Poly  *polys;
  1109. int   *numPolys;
  1110. Seg   *segs;
  1111. int   *numSegs;
  1112. int   *numPointsUsed;
  1113. int   *numPolysUsed;
  1114. int   maxSegs;
  1115. int   *dupPolys;
  1116. int   *invalidPolys;
  1117. /******************************************************************************/
  1118. /* Generate the list of line segments for the polygons.                       */
  1119. /******************************************************************************/
  1120. {
  1121. int index, index2, index3;
  1122. Point **polyPoints;
  1123. Seg *seg;
  1124.  
  1125.    *dupPolys     = 0;
  1126.    *invalidPolys = 0;
  1127.    *numSegs = 0;
  1128.  
  1129.    for(index = 0; index < *numPoints; index++){
  1130.       points[index].index = index;
  1131.       points[index].used = FALSE;
  1132.       }
  1133.  
  1134.    for(index = 0; index < *numPolys; index++){
  1135.  
  1136. /* sort the poly's points by coordinates for finding duplicate polys */
  1137.  
  1138.       qsort(polys[index].points, polys[index].numPoints, sizeof(Point *),
  1139.       pComparePoints);
  1140.  
  1141.       polys[index].used = TRUE; 
  1142.       }
  1143.  
  1144. /* Sort the polygons */
  1145.  
  1146.    qsort(polys, *numPolys, sizeof(Poly), comparePolys);
  1147.  
  1148. /* Tag duplicate polygons */
  1149.  
  1150.    if(*numPolys > 0){
  1151.       *numPolysUsed = 1;
  1152.       }
  1153.  
  1154.    for(index = 1; index < *numPolys; index++){
  1155.  
  1156.       if((comparePolys((&(polys[index - 1])), (&(polys[index])))) == 0){
  1157.          polys[index].used = FALSE;
  1158.          (*dupPolys)++;
  1159.       }else{
  1160.          (*numPolysUsed)++;
  1161.          }
  1162.       }
  1163.  
  1164. /* Assemble valid polygons */
  1165.  
  1166.    for(index = 0; index < *numPolys; index++){
  1167.  
  1168.       if(polys[index].used == TRUE){
  1169.          polyPoints = polys[index].points;
  1170.          polys[index].numSegs = 0;
  1171.  
  1172. /* Sort copy of the poly's points by angle to get another point on the */
  1173. /* convex hull                                                         */
  1174.  
  1175.          dotCompareA = polyPoints[0]; 
  1176.          dotCompareB = polyPoints[1]; 
  1177.  
  1178.          qsort(polyPoints, polys[index].numPoints, sizeof(Point *),
  1179.          dotComparePoints);
  1180.  
  1181. /* Sort copy of the poly's points by angle with respect to the segment on the */
  1182. /* convex hull                                                                */
  1183.  
  1184.          dotCompareB = dotCompareA;
  1185.          dotCompareA = polyPoints[0];
  1186.  
  1187.          qsort(polyPoints, polys[index].numPoints, sizeof(Point *),
  1188.          dotComparePoints);
  1189.  
  1190. /* Set line segments */
  1191.  
  1192.          dotCompareA = polyPoints[0];
  1193.          dotCompareB = polyPoints[polys[index].numPoints - 1];
  1194.  
  1195.          index2 = 1;
  1196.  
  1197.          while(index2 < (polys[index].numPoints - 1)){
  1198.    
  1199.             if(dotComparePoints(&(polyPoints[index2]),
  1200.             &(polyPoints[index2 + 1])) >= 0){
  1201.      
  1202. /* if point would make polygon non-convex delete it */
  1203.  
  1204.                polys[index].numPoints--;
  1205.                for(index3 = index2; index3 < polys[index].numPoints; index3++){
  1206.                   polyPoints[index3] = polyPoints[index3 + 1];
  1207.                   }
  1208.             }else{
  1209.  
  1210. /* keep the point and add the line segment */
  1211.  
  1212.                if(polys[index].numSegs == 0){
  1213.                   if((polys[index].segs = (Seg **)malloc(
  1214.                   sizeof(Seg *))) == NULL){
  1215.                     fprintf(stderr,
  1216.                      "Unable to allocate memory for polygon segments\n");
  1217.                      exit(1);
  1218.                      }
  1219.                }else{
  1220.                   if((polys[index].segs = (Seg **)realloc(polys[index].segs,
  1221.                   (polys[index].numSegs + 1) * sizeof(Seg *))) == NULL){
  1222.                     fprintf(stderr,
  1223.                      "Unable to allocate memory for polygon segments\n");
  1224.                      exit(1);
  1225.                      }
  1226.                   }
  1227.  
  1228.                hash(dotCompareA,polyPoints[index2],segs,numSegs,&seg,
  1229.                maxSegs, &(polys[index]));
  1230.  
  1231.                polys[index].segs[polys[index].numSegs] = seg;
  1232.                (polys[index].numSegs)++;
  1233.  
  1234.                dotCompareB = dotCompareA;
  1235.                dotCompareA = polyPoints[index2];
  1236.  
  1237.                index2++;
  1238.                }               
  1239.             }
  1240.  
  1241. /* Close off the polygon */
  1242.  
  1243.          if(polys[index].numSegs == 0){
  1244.             polys[index].used = FALSE;
  1245.             (*numPolysUsed)--;
  1246.             (*invalidPolys)++;
  1247.          }else{
  1248.     
  1249.             if((polys[index].segs = (Seg **)realloc(polys[index].segs,
  1250.             (polys[index].numSegs + 1) * sizeof(Seg *))) == NULL){
  1251.               fprintf(stderr,
  1252.                "Unable to allocate memory for polygon segments\n");
  1253.                exit(1);
  1254.                }
  1255.  
  1256.             hash(polyPoints[polys[index].numPoints - 2],
  1257.             polyPoints[polys[index].numPoints - 1], segs, numSegs, &seg,
  1258.             maxSegs, &(polys[index]));
  1259.  
  1260.             polys[index].segs[polys[index].numSegs] = seg;
  1261.             (polys[index].numSegs)++;
  1262.  
  1263.             if((polys[index].segs = (Seg **)realloc(polys[index].segs,
  1264.             (polys[index].numSegs + 1) * sizeof(Seg *))) == NULL){
  1265.               fprintf(stderr,
  1266.                "Unable to allocate memory for polygon segments\n");
  1267.                exit(1);
  1268.                }
  1269.  
  1270.             hash(polyPoints[polys[index].numPoints - 1],
  1271.             polyPoints[0], segs, numSegs, &seg,
  1272.             maxSegs, &(polys[index]));
  1273.  
  1274.             polys[index].segs[polys[index].numSegs] = seg;
  1275.             (polys[index].numSegs)++;
  1276.             }
  1277.          }
  1278.       }
  1279.  
  1280.    index2 = 0;
  1281.  
  1282.    for(index = 0; index < *numPoints; index++){
  1283.       if(points[index].used == TRUE){
  1284.          points[index].index = index2++;
  1285.       }else{
  1286.          points[index].index = -1;
  1287.          }
  1288.       }
  1289.  
  1290.    *numPointsUsed = index2;
  1291.  
  1292.    return EOK;
  1293. }
  1294.  
  1295.  
  1296.  
  1297. void main(argc, argv)
  1298. int    argc;
  1299. char **argv;
  1300. {
  1301. FILE *fd;
  1302. Point  *points;
  1303. Poly   *polys;
  1304. Seg    *segs, *segCopy;
  1305. int numPoints, numPolys, numSegs, numPointsUsed, numPolysUsed;
  1306. int maxPoints, maxPolys, maxSegs, numColors;
  1307. int index, index2, dupPolys, invalidPolys;
  1308. unsigned char *red, *green, *blue;
  1309. unsigned int *lookup;
  1310. char line[MAXLINE];
  1311.  
  1312.    if(argc == 4){
  1313.       numColors = atol(argv[3]);
  1314.       (void)sprintf(line, "%d", numColors);
  1315.  
  1316.       if((strcmp(line, argv[3])) || (numColors < 1)){
  1317.          fprintf(stderr, "\n\nIllegal number of colors specified\n");
  1318.          fprintf(stderr, "\nUsage: %s <infile> <outfile> [colors]\n", argv[0]);
  1319.          fprintf(stderr, "Use - for standard input or output\n");
  1320.          fprintf(stderr, "<infile>  File to read polygons from\n");
  1321.          fprintf(stderr, "<outfile> File to write object to\n");
  1322.          fprintf(stderr, "[colors]  Number of colors to quantize to\n");
  1323.          fprintf(stderr, "          Default = %d\n", NUMCOLORS);
  1324.          fprintf(stderr, "\n");
  1325.          exit(1);
  1326.          }
  1327.    }else{
  1328.       numColors = NUMCOLORS;
  1329.       }
  1330.  
  1331.    if((argc != 3) && (argc != 4)){
  1332.       fprintf(stderr, "\nUsage: %s <infile> <outfile> [colors]\n", argv[0]);
  1333.       fprintf(stderr, "Use - for standard input or output\n");
  1334.       fprintf(stderr, "<infile>  File to read polygons from\n");
  1335.       fprintf(stderr, "<outfile> File to write object to\n");
  1336.       fprintf(stderr, "[colors]  Number of colors to quantize to\n");
  1337.       fprintf(stderr, "          Default = %d\n", NUMCOLORS);
  1338.       fprintf(stderr, "\n");
  1339.       exit(1);
  1340.       }
  1341.  
  1342.    maxPoints = 100;
  1343.    maxPolys  = 100;
  1344.  
  1345.    if((points = (Point *)calloc(1, maxPoints * sizeof(Point))) == NULL){
  1346.       fprintf(stderr, "Unable to allocate memory for points\n");
  1347.       exit(1);
  1348.       }
  1349.  
  1350.    if((polys = (Poly *)calloc(1, maxPolys * sizeof(Poly))) == NULL){
  1351.       fprintf(stderr, "Unable to allocate memory for polygons\n");
  1352.       exit(1);
  1353.       }
  1354.  
  1355.    readPoints(argv[1], &points, &numPoints, &polys, &numPolys, &maxPoints,
  1356.    &maxPolys);
  1357.  
  1358.    if(numPolys <= 0){
  1359.       fprintf(stderr, "No polygons found from which to generate object.\n");
  1360.       exit(1);
  1361.       }
  1362.  
  1363.    eliminateDuplicatePoints(&points, &numPoints, polys, numPolys, &maxPoints);
  1364.  
  1365.    maxSegs = (maxPoints * 3) / 2;
  1366.  
  1367.    if((segs = (Seg *)calloc(1, maxSegs * sizeof(Seg))) == NULL){
  1368.       fprintf(stderr, "Unable to allocate memory for segments\n");
  1369.       exit(1);
  1370.       }
  1371.  
  1372.    hashTableSize = maxSegs + 1000;
  1373.  
  1374.    if((hashTable = (hashEntry *)calloc(1, hashTableSize *  sizeof(hashEntry)))
  1375.    == NULL){
  1376.       fprintf(stderr, "Unable to allocate memory for hashTable\n");
  1377.       exit(1);
  1378.       }
  1379.  
  1380.    formPolygonsSegments(points, &numPoints, polys, &numPolys, segs, &numSegs,
  1381.    &numPointsUsed, &numPolysUsed, maxSegs, &dupPolys, &invalidPolys);
  1382.  
  1383.    if(strcmp(argv[2], "-") == 0){
  1384.       fd = stdout;
  1385.    }else{
  1386.       if((fd = fopen(argv[2], "w")) == NULL){
  1387.          fprintf(stderr, "Unable to open file %s\n", argv[2]);
  1388.          exit(1);
  1389.          }
  1390.       }
  1391.  
  1392.    if((red = (unsigned char *)calloc(1, numPolys * sizeof(unsigned char)))
  1393.     == NULL){
  1394.       fprintf(stderr, "Unable to allocate memory for red\n");
  1395.       exit(1);
  1396.       }
  1397.  
  1398.    if((green = (unsigned char *)calloc(1, numPolys * sizeof(unsigned char)))
  1399.     == NULL){
  1400.       fprintf(stderr, "Unable to allocate memory for green\n");
  1401.       exit(1);
  1402.       }
  1403.  
  1404.    if((blue = (unsigned char *)calloc(1, numPolys * sizeof(unsigned char)))
  1405.     == NULL){
  1406.       fprintf(stderr, "Unable to allocate memory for blue\n");
  1407.       exit(1);
  1408.       }
  1409.   
  1410.    if((lookup = (unsigned int*)calloc(1, numPolys *
  1411.    sizeof(unsigned int))) == NULL){
  1412.       fprintf(stderr, "Unable to allocate memory for lookup\n");
  1413.       exit(1);
  1414.       }
  1415.  
  1416.    for(index = 0; index < numPolys; index++){
  1417.       red[index]   = polys[index].red; 
  1418.       green[index] = polys[index].green; 
  1419.       blue[index]  = polys[index].blue;
  1420.       }
  1421.  
  1422.    quantize(red, green, blue, numPolys, &numColors, lookup);
  1423.  
  1424.    fprintf(fd, "# x3d 2.0 object file generated by %s using %s\n\n", 
  1425.    argv[0], argv[1]);
  1426.    fprintf(fd, "# number of colors used in object\n\n");
  1427.    fprintf(fd, "%d\n\n", numColors);
  1428.    fprintf(fd, "# colors used in object (color number red green blue)\n\n");
  1429.  
  1430.    for(index = 0; index < numColors; index++){
  1431.       fprintf(fd, "%d   %d %d %d\n", index, red[index],
  1432.       green[index], blue[index]);
  1433.       }
  1434.  
  1435.    free(red);
  1436.    free(green);
  1437.    free(blue);
  1438.  
  1439.    fprintf(fd, "\n");
  1440.    fprintf(fd, "# number of points used in object\n\n");
  1441.    fprintf(fd, "%d\n\n", numPointsUsed);
  1442.    fprintf(fd, "# points used in object (point number x y z)\n\n");
  1443.  
  1444.    for(index = 0; index < numPoints; index++){
  1445.       if(points[index].used == TRUE){
  1446.          fprintf(fd, "%d   %f %f %f\n", points[index].index, 
  1447.         (float)points[index].x, (float)points[index].y, (float)points[index].z);
  1448.          }
  1449.       }
  1450.  
  1451. /* Sort segments and polygons by color in hope of reducing foreground */
  1452. /* color changes in x3d */
  1453.  
  1454.    for(index = 0; index < numSegs; index++){
  1455.       segs[index].lookupColor = lookup[segs[index].poly - polys];
  1456.       } 
  1457.  
  1458.    if((segCopy = (Seg *)calloc(1, numSegs * sizeof(Seg))) == NULL){
  1459.       fprintf(stderr, "Unable to allocate memory for segCopy\n");
  1460.       exit(1);
  1461.       }
  1462.  
  1463.    memcpy(segCopy, segs, numSegs * sizeof(Seg));
  1464.  
  1465.    qsort(segCopy, numSegs, sizeof(Seg), colorCompareSegs);
  1466.  
  1467.    for(index = 0; index < numSegs; index++){
  1468.       segs[segCopy[index].index].index = index;
  1469.       segCopy[index].index = index;
  1470.       } 
  1471.  
  1472.    for(index = 0; index < numPolys; index++){
  1473.       polys[index].lookupColor = lookup[index];
  1474.       } 
  1475.  
  1476.    qsort(polys, numPolys, sizeof(Poly), colorComparePolys);
  1477.  
  1478.    fprintf(fd, "\n");
  1479.    fprintf(fd, "# number segments used in object\n\n");
  1480.    fprintf(fd, "%d\n\n", numSegs);
  1481.    fprintf(fd, "# segments used in object (segment number color p q)\n\n");
  1482.  
  1483.    for(index = 0; index < numSegs; index++){
  1484.       fprintf(fd, "%d   %d   %d %d\n", segCopy[index].index, 
  1485.       segCopy[index].lookupColor, segCopy[index].p->index,
  1486.       segCopy[index].q->index);
  1487.       }
  1488.  
  1489.    fprintf(fd, "\n");
  1490.    fprintf(fd, "# number of polygons used in object\n\n");
  1491.    fprintf(fd, "%d\n\n", numPolysUsed);
  1492.    fprintf(fd, "# polygons used in object (polygon number color\n");
  1493.    fprintf(fd, "# number of segments in polygon s0 s1 s2 ... sn)\n\n");
  1494.  
  1495.    index2 = 0;
  1496.  
  1497.    for(index = 0; index < numPolys; index++){
  1498.       if(polys[index].used == TRUE){
  1499.          fprintf(fd, "%d   %d   %d  ", index, polys[index].lookupColor,
  1500.          polys[index].numSegs);
  1501.          index2++;
  1502.          for(index2 = 0; index2 < polys[index].numSegs; index2++){
  1503.             fprintf(fd, " %d", polys[index].segs[index2]->index);
  1504.             } 
  1505.          fprintf(fd, "\n");
  1506.          } 
  1507.       } 
  1508.  
  1509.    fprintf(stderr, "\nConversion successful:\n\n");
  1510.    fprintf(stderr, "%8d colors used in object.\n",       numColors);
  1511.    fprintf(stderr, "%8d points used in object\n",        numPointsUsed);
  1512.    fprintf(stderr, "%8d segments used in object\n",      numSegs);
  1513.    fprintf(stderr, "%8d polygons used in object\n",      numPolysUsed);
  1514.    fprintf(stderr, "%8d duplicate polygons removed.\n",  dupPolys);
  1515.    fprintf(stderr, "%8d invalid polygons removed.\n",    invalidPolys);
  1516.    fprintf(stderr, "\n");
  1517. }
  1518.