home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / GRAPHICS / MISC / PVQUAN15.ZIP / QUANT.ZIP / HECKBERT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-03  |  16.4 KB  |  565 lines
  1. /************************************************************************
  2.  *                                                                      *
  3.  *                  Copyright (c) 1991, Frank van der Hulst             *
  4.  *                          All Rights Reserved                         *
  5.  *                                                                      *
  6.  * Authors:                                                             *
  7.  *        FvdH - Frank van der Hulst (Wellington, NZ)                   *
  8.  *                                                                      *
  9.  * Versions:                                                            *
  10.  *      V1.1 910626 FvdH - QUANT released for DBW_RENDER                *
  11.  *      V1.2 911021 FvdH - QUANT released for PoV Ray                   *
  12.  *      V1.4 920303 FvdH - Ported to GNU                                *
  13.  *                                                                      *
  14.  ************************************************************************/
  15. /*
  16.  * This software is copyrighted as noted below.  It may be freely copied,
  17.  * modified, and redistributed, provided that the copyright notice is
  18.  * preserved on all copies.
  19.  *
  20.  * There is no warranty or other guarantee of fitness for this software,
  21.  * it is provided solely "as is".  Bug reports or fixes may be sent
  22.  * to the author, who may or may not act on them as he desires.
  23.  *
  24.  * You may not include this software in a program or other software product
  25.  * without supplying the source, or without informing the end-user that the
  26.  * source is available for no extra charge.
  27.  *
  28.  * If you modify this software, you should include a notice giving the
  29.  * name of the person performing the modification, the date of modification,
  30.  * and the reason for such modification.
  31. */
  32. /*
  33.  * colorquant.c
  34.  *
  35.  * Perform variance-based color quantization on a "full color" image.
  36.  * Author:    Craig Kolb
  37.  *        Department of Mathematics
  38.  *        Yale University
  39.  *        kolb@yale.edu
  40.  * Date:    Tue Aug 22 1989
  41.  * Copyright (C) 1989 Craig E. Kolb
  42.  * $Id: colorquant.c,v 1.3 89/12/03 18:27:16 craig Exp $
  43.  *
  44.  * $Log:    colorquant.c,v $
  45.  *
  46.  * Revision 1.4  91/06/26  16:00:00  Frank van der Hulst
  47.  * Ported to Turbo C;
  48.  *       Call farmalloc rather than malloc
  49.  *       Virtual memory added to swap box data to/from disk
  50.  *       Rewritten in ANSI C
  51.  *       Removed call to QuantHistogram() from colorquant, to allow two
  52.  *       image files to create one palette
  53.  *       Changed QuantHistogram() to read from file, rather than from an
  54.  *       array
  55.  *          Changed format of palette to conform with the VGA palette
  56.  *
  57.  * Revision 1.3  89/12/03  18:27:16  craig
  58.  * Removed bogus integer casts in distance calculation in makenearest().
  59.  *
  60.  * Revision 1.2  89/12/03  18:13:12  craig
  61.  * FindCutpoint now returns FALSE if the given box cannot be cut.  This
  62.  * to avoid overflow problems in CutBox.
  63.  * "whichbox" in GreatestVariance() is now initialized to 0.
  64.  *
  65.  */
  66.  
  67. #ifdef __TURBOC__
  68. #include <mem.h>
  69. #define HUGE 1.79e308
  70. #endif
  71. #include <math.h>
  72.  
  73. #include "quant.h"
  74. #include "heckbert.h"
  75.  
  76. #define MAX(x,y)    ((x) > (y) ? (x) : (y))
  77.  
  78. static unsigned long        NPixels = 0L;            /* total # of pixels */
  79.  
  80. static int neighbours[MAXCOLORS];
  81.  
  82. /*
  83.  * Compute the histogram of the image as well as the projected frequency
  84.  * arrays for the first world-encompassing box.
  85.  * We compute both the histogram and the proj. frequencies of
  86.  * the first box at the same time to save a pass through the
  87.  * entire image.
  88.  * The projected frequency arrays of the largest box are zeroed out as
  89.  * as part of open_box_file(), called from main().
  90.  */
  91.  
  92. void QuantHistogram(Box *box)
  93. {
  94. unsigned long *rf, *gf, *bf;
  95. UCHAR pixel[3];
  96.  
  97.     rf = box->freq[RED];
  98.     gf = box->freq[GREEN];
  99.     bf = box->freq[BLUE];
  100.  
  101.     while (get_pixel(pixel)) {
  102.         rf[pixel[RED]]++;
  103.         gf[pixel[GREEN]]++;
  104.         bf[pixel[BLUE]]++;
  105.         Histogram[(((pixel[RED]<<INPUT_BITS)|pixel[GREEN])<<INPUT_BITS)|pixel[BLUE]]++;
  106.         NPixels++;
  107.     }
  108. }
  109.  
  110. /*
  111.  * Compute mean and weighted variance of the given box.
  112.  */
  113. void BoxStats(Box HUGE_PTR box)
  114. {
  115. int i, color;
  116. unsigned long *freq;
  117. double mean, var;
  118.  
  119.     if(box->weight == 0) {
  120.         box->weightedvar = 0;
  121.         return;
  122.     }
  123.  
  124.     box->weightedvar = 0.0;
  125.     for (color = 0; color < 3; color++) {
  126.         var = mean = 0;
  127.         i = box->low[color];
  128.         freq = &box->freq[color][i];
  129.         for (; i < box->high[color]; i++, freq++) {
  130.             mean += i * *freq;
  131.             var += i*i* *freq;
  132.         }
  133.         box->mean[color] = mean / box->weight;
  134.         box->weightedvar += var - box->mean[color]*box->mean[color]*box->weight;
  135.     }
  136.     box->weightedvar /= NPixels;
  137. }
  138.  
  139. /*
  140.  * Return the number of the box in 'boxes' with the greatest variance.
  141.  * Restrict the search to those boxes with indices between 0 and n-1.
  142.  */
  143. int GreatestVariance(int n)
  144. {
  145.     int i, whichbox = 0;
  146.     double max;
  147.     Box *box;
  148.  
  149.     max = -1;
  150.     for (i = 0; i < n; i++) {
  151.         box = get_box_tmp(i);
  152.         if (box->weightedvar > max) {
  153.             max = box->weightedvar;
  154.             whichbox = i;
  155.         }
  156.     }
  157.     return whichbox;
  158. }
  159.  
  160. /*
  161.  * Update projected frequency arrays for two boxes which used to be
  162.  * a single box.
  163.  */
  164. void UpdateFrequencies(Box HUGE_PTR box1, Box HUGE_PTR box2)
  165. {
  166. unsigned long myfreq, HUGE_PTR h;
  167. int b, g, r;
  168. int roff;
  169.  
  170.     memset(box1->freq[0], 0, IN_COLOURS * sizeof(unsigned long));
  171.     memset(box1->freq[1], 0, IN_COLOURS * sizeof(unsigned long));
  172.     memset(box1->freq[2], 0, IN_COLOURS * sizeof(unsigned long));
  173.  
  174.     for (r = box1->low[0]; r < box1->high[0]; r++) {
  175.         roff = r << INPUT_BITS;
  176.         for (g = box1->low[1];g < box1->high[1]; g++) {
  177.             b = box1->low[2];
  178.             h = Histogram + (((roff | g) << INPUT_BITS) | b);
  179.             for (; b < box1->high[2]; b++) {
  180.                 if ((myfreq = *h++) == 0)
  181.                     continue;
  182.                 box1->freq[0][r] += myfreq;
  183.                 box1->freq[1][g] += myfreq;
  184.                 box1->freq[2][b] += myfreq;
  185.                 box2->freq[0][r] -= myfreq;
  186.                 box2->freq[1][g] -= myfreq;
  187.                 box2->freq[2][b] -= myfreq;
  188.             }
  189.         }
  190.     }
  191. }
  192.  
  193. /*
  194.  * Compute the 'optimal' cutpoint for the given box along the axis
  195.  * indicated by 'color'.  Store the boxes which result from the cut
  196.  * in newbox1 and newbox2.
  197.  */
  198. int FindCutpoint(Box HUGE_PTR box, int color, Box HUGE_PTR newbox1, Box HUGE_PTR newbox2)
  199. {
  200.     double u, v, max;
  201.     int i, maxindex, minindex, cutpoint;
  202.     unsigned long optweight, curweight;
  203.  
  204.     if (box->low[color] + 1 == box->high[color])
  205.         return FALSE;    /* Cannot be cut. */
  206.     minindex = (int)((box->low[color] + box->mean[color]) * 0.5);
  207.     maxindex = (int)((box->mean[color] + box->high[color]) * 0.5);
  208.  
  209.     cutpoint = minindex;
  210.     optweight = box->weight;
  211.  
  212.     curweight = 0.;
  213.     for (i = box->low[color] ; i < minindex ; i++)
  214.         curweight += box->freq[color][i];
  215.     u = 0.;
  216.     max = -1;
  217.     for (i = minindex; i <= maxindex ; i++) {
  218.         curweight += box->freq[color][i];
  219.         if (curweight == box->weight)
  220.             break;
  221.         u += (double)(i * box->freq[color][i]) / box->weight;
  222.         v = ((double)curweight / (box->weight-curweight)) *
  223.                 (box->mean[color]-u)*(box->mean[color]-u);
  224.         if (v > max) {
  225.             max = v;
  226.             cutpoint = i;
  227.             optweight = curweight;
  228.         }
  229.     }
  230.     cutpoint++;
  231.     *newbox1 = *newbox2 = *box;
  232.     newbox1->weight = optweight;
  233.     newbox2->weight -= optweight;
  234.     newbox1->high[color] = cutpoint;
  235.     newbox2->low[color] = cutpoint;
  236.     UpdateFrequencies(newbox1, newbox2);
  237.     BoxStats(newbox1);
  238.     BoxStats(newbox2);
  239.  
  240.     return TRUE;    /* Found cutpoint. */
  241. }
  242.  
  243. /*
  244.  * Cut the given box.  Returns TRUE if the box could be cut, FALSE otherwise.
  245.  */
  246. int CutBox(Box HUGE_PTR box, Box HUGE_PTR newbox)
  247. {
  248.     int i;
  249.     double totalvar[3];
  250.     static Box newboxes[3][2];  /* Only used by CutBox, but don't want it on stack */
  251.  
  252.     if (box->weightedvar == 0. || box->weight == 0)
  253.         /*
  254.          * Can't cut this box.
  255.          */
  256.         return FALSE;
  257.  
  258.     /*
  259.      * Find 'optimal' cutpoint along each of the red, green and blue
  260.      * axes.  Sum the variances of the two boxes which would result
  261.      * by making each cut and store the resultant boxes for
  262.      * (possible) later use.
  263.      */
  264.     for (i = 0; i < 3; i++) {
  265.         if (FindCutpoint(box, i, &newboxes[i][0], &newboxes[i][1]))
  266.             totalvar[i] = newboxes[i][0].weightedvar +
  267.                 newboxes[i][1].weightedvar;
  268.         else
  269.             totalvar[i] = HUGE;
  270.     }
  271.  
  272.     /*
  273.      * Find which of the three cuts minimized the total variance
  274.      * and make that the 'real' cut.
  275.      */
  276.     if (totalvar[RED] <= totalvar[GREEN] &&
  277.         totalvar[RED] <= totalvar[BLUE]) {
  278.         *box = newboxes[RED][0];
  279.         *newbox = newboxes[RED][1];
  280.     } else if (totalvar[GREEN] <= totalvar[RED] &&
  281.          totalvar[GREEN] <= totalvar[BLUE]) {
  282.         *box = newboxes[GREEN][0];
  283.         *newbox = newboxes[GREEN][1];
  284.     } else {
  285.         *box = newboxes[BLUE][0];
  286.         *newbox = newboxes[BLUE][1];
  287.     }
  288.  
  289.     return TRUE;
  290. }
  291.  
  292. /*
  293.  * Iteratively cut the boxes.
  294.  */
  295. CutBoxes(int colors)
  296. {
  297. int curbox, varbox;
  298. Box *box = get_box(0);
  299.  
  300.     box->low[RED] = box->low[GREEN] = box->low[BLUE] = 0;
  301.     box->high[RED] = box->high[GREEN] = box->high[BLUE] = IN_COLOURS;
  302.     box->weight = NPixels;
  303.  
  304.     BoxStats(box);
  305.     free_box(0);
  306.  
  307.     printf("%d Boxes: cutting box: ", colors);
  308.     for (curbox = 1; curbox < colors; curbox++) {
  309.         printf("%3d\b\b\b", curbox);
  310.         varbox = GreatestVariance(curbox);
  311.         if (CutBox(get_box(varbox), get_box(curbox)) == FALSE)
  312.                 break;
  313.         free_box(curbox);
  314.         free_box(varbox);
  315.     }
  316.     printf("Done\n");
  317.  
  318.     return curbox;
  319. }
  320.  
  321. /*
  322.  * Make the centroid of "boxnum" serve as the representative for
  323.  * each color in the box.
  324.  */
  325. void SetRGBmap(int boxnum, Box *box, int bits)
  326. {
  327. int r, g, b;
  328.  
  329.     for (r = box->low[RED]; r < box->high[RED]; r++) {
  330.         for (g = box->low[GREEN]; g < box->high[GREEN]; g++) {
  331.             for (b = box->low[BLUE]; b < box->high[BLUE]; b++) {
  332.                 RGBmap[(((r<<bits)|g)<<bits)|b]=(char)boxnum;
  333.             }
  334.         }
  335.     }
  336. }
  337.  
  338. /*
  339.  * In order to minimize our search for 'best representative', we form the
  340.  * 'neighbors' array.  This array contains the number of the boxes whose
  341.  * centroids *might* be used as a representative for some color in the
  342.  * current box.  We need only consider those boxes whose centroids are closer
  343.  * to one or more of the current box's corners than is the centroid of the
  344.  * current box. 'Closeness' is measured by Euclidean distance.
  345.  */
  346. int getneighbors(int num, int colors)
  347. {
  348.     int i, j;
  349.     Box *bp;
  350.     double dist, LowR, LowG, LowB, HighR, HighG, HighB, ldiff, hdiff;
  351.  
  352.     bp = get_box_tmp(num);
  353.  
  354.     ldiff = bp->low[RED] - bp->mean[RED];
  355.     ldiff *= ldiff;
  356.     hdiff = bp->high[RED] - bp->mean[RED];
  357.     hdiff *= hdiff;
  358.     dist = MAX(ldiff, hdiff);
  359.  
  360.     ldiff = bp->low[GREEN] - bp->mean[GREEN];
  361.     ldiff *= ldiff;
  362.     hdiff = bp->high[GREEN] - bp->mean[GREEN];
  363.     hdiff *= hdiff;
  364.     dist += MAX(ldiff, hdiff);
  365.  
  366.     ldiff = bp->low[BLUE] - bp->mean[BLUE];
  367.     ldiff *= ldiff;
  368.     hdiff = bp->high[BLUE] - bp->mean[BLUE];
  369.     hdiff *= hdiff;
  370.     dist += MAX(ldiff, hdiff);
  371.  
  372.     dist = sqrt(dist);
  373.  
  374.     /*
  375.      * Loop over all colors in the colormap, the ith entry of which
  376.      * corresponds to the ith box.
  377.      *
  378.      * If the centroid of a box is as close to any corner of the
  379.      * current box as is the centroid of the current box, add that
  380.      * box to the list of "neighbors" of the current box.
  381.      */
  382.     HighR = (double)bp->high[RED] + dist;
  383.     HighG = (double)bp->high[GREEN] + dist;
  384.     HighB = (double)bp->high[BLUE] + dist;
  385.     LowR = (double)bp->low[RED] - dist;
  386.     LowG = (double)bp->low[GREEN] - dist;
  387.     LowB = (double)bp->low[BLUE] - dist;
  388.     for (i = j = 0; i < colors; i++) {
  389.         bp = get_box_tmp(i);
  390.         if (LowR <= bp->mean[RED] && HighR >= bp->mean[RED] &&
  391.             LowG <= bp->mean[GREEN] && HighG >= bp->mean[GREEN] &&
  392.             LowB <= bp->mean[BLUE] && HighB >= bp->mean[BLUE])
  393.             neighbours[j++] = i;
  394.     }
  395.  
  396.     return j;    /* Return the number of neighbors found. */
  397. }
  398.  
  399. /*
  400.  * Assign representative colors to every pixel in a given box through
  401.  * the construction of the NearestColor array.  For each color in the
  402.  * given box, we look at the list of neighbors passed to find the
  403.  * one whose centroid is closest to the current color.
  404.  */
  405. void makenearest(int boxnum, int nneighbors, int bits)
  406. {
  407.     int n, b, g, r;
  408.     double rdist, gdist, bdist, dist, mindist;
  409.     int which, *np;
  410.     Box *box, *bp;
  411.  
  412.     box = get_box(boxnum);
  413.  
  414.     for (r = box->low[RED]; r < box->high[RED]; r++) {
  415.         for (g = box->low[GREEN]; g < box->high[GREEN]; g++) {
  416.             for (b = box->low[BLUE]; b < box->high[BLUE]; b++) {
  417. /*
  418.  * The following two lines should be commented out if the RGBmap is going
  419.  * to be used for images other than the one given.
  420.  */
  421.                 if (Histogram[(((r<<bits)|g)<<bits)|b] == 0)
  422.                     continue;
  423.                 mindist = HUGE;
  424.                 /*
  425.                  * Find the colormap entry which is
  426.                  * closest to the current color.
  427.                  */
  428.                 np = neighbours;
  429.                 for (n = 0; n < nneighbors; n++, np++) {
  430.                     bp = get_box_tmp(*np);
  431.                     rdist = r - bp->mean[RED];
  432.                     gdist = g - bp->mean[GREEN];
  433.                     bdist = b - bp->mean[BLUE];
  434.                     dist = rdist*rdist + gdist*gdist + bdist*bdist;
  435.                     if (dist < mindist) {
  436.                         mindist = dist;
  437.                         which = *np;
  438.                     }
  439.                 }
  440.                 /*
  441.                  * The colormap entry closest to this
  442.                  * color is used as a representative.
  443.                  */
  444.                 RGBmap[(((r<<bits)|g)<<bits)|b] = which;
  445.             }
  446.         }
  447.     }
  448.     free_box(boxnum);
  449. }
  450. /*
  451.  * Form colormap and NearestColor arrays.
  452.  */
  453. void find_colors(int colors, int bits)
  454. {
  455. int i;
  456. int num;
  457.  
  458.     /*
  459.      * Form map of representative (nearest) colors.
  460.      */
  461.     printf("Mapping colours for box: ");
  462.     for (i = 0; i < colors; i++) {
  463.         printf("%3d\b\b\b", i);
  464.         /*
  465.          * Create list of candidate neighbors and
  466.          * find closest representative for each
  467.          * color in the box.
  468.          */
  469.         num = getneighbors(i, colors);
  470.         makenearest(i, num, bits);
  471.     }
  472.     printf("Done\n");
  473. }
  474.  
  475. /*
  476.  * Compute RGB to colormap index map.
  477.  */
  478. void ComputeRGBMap(int colors, int bits, int fast)
  479. {
  480. int i;
  481.  
  482.     if (fast) {
  483.         /*
  484.          * The centroid of each box serves as the representative
  485.          * for each color in the box.
  486.          */
  487.         for (i = 0; i < colors; i++)
  488.             SetRGBmap(i, get_box_tmp(i), bits);
  489.     } else
  490.         /*
  491.          * Find the 'nearest' representative for each pixel.
  492.          */
  493.         find_colors(colors, bits);
  494. }
  495.  
  496. /*
  497.  * Perform variance-based color quantization on a 24-bit image.
  498.  *
  499.  * Input consists of:
  500.  *    in_file    Pointer to file containing of red, green and blue pixel
  501.  *                intensities stored as unsigned characters.
  502.  *                The color of the ith pixel is given consecutive bytes, in the
  503.  *                order red, then green, then blue. Only the LS 4 bits are used.
  504.  *                0 indicates zero intensity, 15 full intensity.
  505.  *    pixels    The length of the red, green and blue arrays
  506.  *                in bytes, stored as an unsigned long int.
  507.  *    colormap    Points to the colormap.  The colormap
  508.  *                consists of red, green and blue arrays.
  509.  *                The red/green/blue values of the ith
  510.  *                colormap entry are given respectively by
  511.  *                colormap[0][i], colormap[1][i] and
  512.  *                colormap[2][i].  Each entry is an unsigned char.
  513.  *    colors    The number of colormap entries, stored
  514.  *                as an integer.
  515.  *    bits        The number of significant bits in each entry
  516.  *                of the red, green and blue arrays. An integer.
  517.  *    rgbmap    An array of unsigned chars of size (2^bits)^3.
  518.  *                This array is used to map from pixels to
  519.  *                colormap entries.  The 'prequantized' red,
  520.  *                green and blue components of a pixel
  521.  *                are used as an index into rgbmap to retrieve
  522.  *                the index which should be used into the colormap
  523.  *                to represent the pixel.  In short:
  524.  *                index = rgbmap[(((r << bits) | g) << bits) | b];
  525.  *     fast    If non-zero, the rgbmap will be constructed
  526.  *                quickly.  If zero, the rgbmap will be built
  527.  *                much slower, but more accurately.  In most
  528.  *                cases, fast should be non-zero, as the error
  529.  *                introduced by the approximation is usually
  530.  *                small.  'Fast' is stored as an integer.
  531.  *    Cfactor    Conversion factor.
  532.  *
  533.  * colorquant returns the number of colors to which the image was
  534.  * quantized.
  535.  */
  536.  
  537. int colorquant(int colors, int bits, int fast, double Cfactor)
  538. {
  539. int    i,                        /* Counter */
  540.         OutColors;            /* # of entries computed */
  541. Box    *box;
  542.  
  543.     OutColors = CutBoxes(colors);
  544.     /*
  545.      * We now know the set of representative colors.  We now
  546.      * must fill in the colormap and convert the representatives
  547.      * from their 'prequantized' range to 0-FULLINTENSITY.
  548.      */
  549.     for (i = 0; i < OutColors; i++) {
  550.         box = get_box_tmp(i);
  551.         palette[i][RED]   = (char)(box->mean[RED] * Cfactor + 0.5);
  552.         palette[i][GREEN] = (char)(box->mean[GREEN] * Cfactor + 0.5);
  553.         palette[i][BLUE]  = (char)(box->mean[BLUE] * Cfactor + 0.5);
  554.     }
  555.  
  556.     ComputeRGBMap(OutColors, bits, fast);
  557.  
  558.     return OutColors;        /* Return # of colormap entries */
  559. }
  560.  
  561. int pal_index(UCHAR *pixel)
  562. {
  563.     return RGBmap[(((pixel[RED]<<INPUT_BITS)|pixel[GREEN])<<INPUT_BITS)|pixel[BLUE]];
  564. }
  565.