home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / fbm / src / fbquant.c < prev    next >
C/C++ Source or Header  |  1990-06-24  |  29KB  |  1,007 lines
  1. /****************************************************************
  2.  * fbquant.c: FBM Release 1.0 25-Feb-90 Michael Mauldin
  3.  *
  4.  * Copyright (C) 1989,1990 by Michael Mauldin.  Permission is granted
  5.  * to use this file in whole or in part for any purpose, educational,
  6.  * recreational or commercial, provided that this copyright notice
  7.  * is retained unchanged.  This software is available to all free of
  8.  * charge by anonymous FTP and in the UUNET archives.
  9.  *
  10.  * fbquant: Convert an RGB color image to mapped color format (color
  11.  *        quantization step).  Floyd-Steinberg dithering is used
  12.  *        to reduce color banding.  The quantization used is a
  13.  *        modification of Heckbert's median cut.
  14.  *
  15.  * USAGE
  16.  *    % fbquant [ -c<numcolors> ] [ -r<bits> ] [ -m<mapimage> ]
  17.  *          [ -n ] [ -<type> ] < rgb > mapped"
  18.  *
  19.  * EDITLOG
  20.  *    LastEditDate = Mon Jun 25 00:04:18 1990 - Michael Mauldin
  21.  *    LastFileName = /usr2/mlm/src/misc/fbm/fbquant.c
  22.  *
  23.  * HISTORY
  24.  * 25-Jun-90  Michael Mauldin (mlm@cs.cmu.edu) Carnegie Mellon
  25.  *    Package for Release 1.0
  26.  *
  27.  * 04-Sep-89  Michael Mauldin (mlm) at Carnegie Mellon University
  28.  *    Add option for low-res colors, default assumes each RGB value
  29.  *    can vary from 0 to 255.  The default is 8 bits for more than 32
  30.  *    color images, and 4 bits for images with 32 or fewer colors.
  31.  *    For example, the Amiga has 4bit DACs and can only display
  32.  *     4^3 == 4096 colors.
  33.  *
  34.  * 03-May-89  Michael Mauldin (mlm) at Carnegie Mellon University
  35.  *    Clear FBM pointers before allocate
  36.  *
  37.  * 07-Apr-89  Michael Mauldin (mlm) at Carnegie Mellon University
  38.  *    Add -m<image> flag for using colormap from another image.
  39.  *
  40.  * 07-Mar-89  Michael Mauldin (mlm) at Carnegie Mellon University
  41.  *    Beta release (version 0.93) mlm@cs.cmu.edu
  42.  *
  43.  * 26-Feb-89  Michael Mauldin (mlm) at Carnegie Mellon University
  44.  *    Changes for small color maps.  Fixed bug with unsigned
  45.  *    arithmetic that ruined dithering for images with small
  46.  *    colormaps.  Added error limiting in the Floyd-Steinberg
  47.  *    code to prevent color "shadowing" that occurs with small
  48.  *    numbers of colors.  Also change to use colors 0..n-1 instead
  49.  *    of reserving colors 0 and n-1 for Sun foreground/background
  50.  *    colors, unless output type is SUN.
  51.  *
  52.  * 11-Nov-88  Michael Mauldin (mlm) at Carnegie Mellon University
  53.  *    Created.
  54.  *
  55.  * References: Uses a variant of Heckbert's adaptive partitioning
  56.  *           algorithm.  See Computer Graphics, v16n3 July 1982
  57.  *
  58.  ****************************************************************/
  59.  
  60. # include <stdio.h>
  61. # include "fbm.h"
  62.  
  63. int cmp_red(), cmp_grn(), cmp_blu(), cmp_cmap(), cmp_int();
  64.  
  65. # define RD 0
  66. # define GR 1
  67. # define BL 2
  68. # define REDMASK 0076000
  69. # define REDSHFT 10
  70. # define GRNMASK 0001740
  71. # define GRNSHFT 5
  72. # define BLUMASK 0000037
  73. # define BLUSHFT 0
  74. # define CUBITS  5
  75. # define CUBIGN  (8-CUBITS)
  76. # define CUBSID  32
  77. # define CUBSIZ  32768
  78. # define MAXSHRT 32767
  79. # define MAXERR     32
  80.  
  81. # define GETR(X) (((X) & REDMASK) >> REDSHFT)
  82. # define GETG(X) (((X) & GRNMASK) >> GRNSHFT)
  83. # define GETB(X)  ((X) & BLUMASK)
  84.  
  85. # define CLRINDEX(R,G,B)            \
  86.     (((R) << REDSHFT) & REDMASK |        \
  87.      ((G) << GRNSHFT) & GRNMASK |        \
  88.      ((B)  & BLUMASK))
  89.  
  90. # define CLRINDEX8(R,G,B)            \
  91.     (((R) << (REDSHFT-CUBIGN)) & REDMASK |    \
  92.      ((G) << (GRNSHFT-CUBIGN)) & GRNMASK |    \
  93.      ((B) >> (CUBIGN))  & BLUMASK)
  94.  
  95. # define GETR8(X) (((X) & REDMASK) >> (REDSHFT-CUBIGN))
  96. # define GETG8(X) (((X) & GRNMASK) >> (GRNSHFT-CUBIGN))
  97. # define GETB8(X) (((X) & BLUMASK) << CUBIGN)
  98.  
  99. # define abs(X) (((X)<0)?-(X):(X))
  100.  
  101. typedef struct cstruct {
  102.     unsigned char rd, gr, bl, indx;
  103. } COLOR;
  104.  
  105. COLOR *cmap = NULL;
  106.  
  107. typedef struct pix_struct {
  108.     short cnt;
  109.     short color;
  110. } PIXEL;
  111.  
  112. int debug=0, verbose=0, colors=256, showcolor=0, dither=1;
  113. int resbits = -1, resmask = 0xff;
  114. double flesh = 1.0, atof();
  115. int fr = 192, fg = 96, fb = 80;
  116.  
  117. /****************************************************************
  118.  * main
  119.  ****************************************************************/
  120.  
  121. # define USAGE \
  122. "Usage:  fbquant [ -c<numcolors> ] [ -r<bits> ] [ -m<mapimage> ]\n\
  123.         [ -f<flesh> ] [ -n ] [ -<type> ] < rgb > mapped"
  124.  
  125. #ifndef lint
  126. static char *fbmid = 
  127. "$FBM fbquant.c <1.0> 25-Jun-90  (C) 1989,1990 by Michael Mauldin, source \
  128. code available free from MLM@CS.CMU.EDU and from UUNET archives$";
  129. #endif
  130.  
  131. main (argc, argv)
  132. char *argv[];
  133. { FBM input, output, mapimage;    /* Images */
  134.   int hist[CUBSIZ];        /* Color cube 32x32x32 for histogram */
  135.   int outtype = DEF_8BIT;    /* Output format desired */
  136.   char *mapfile = NULL;        /* Name of file to copy map from */
  137.   register int c;
  138.  
  139.   /* Clear pointers */
  140.   input.bm = input.cm = (unsigned char *) NULL;
  141.   output.bm = output.cm = (unsigned char *) NULL;
  142.   mapimage.bm = mapimage.cm = (unsigned char *) NULL;
  143.  
  144.   /* Get the options */
  145.   while (--argc > 0 && (*++argv)[0] == '-')
  146.   { while (*++(*argv))
  147.     { switch (**argv)
  148.       { case 'c':    colors = atoi (*argv+1); SKIPARG; break;
  149.         case 'f':    flesh = atof (*argv+1); SKIPARG; break;
  150.         case 'r':    resbits = atoi (*argv+1); SKIPARG; break;
  151.         case 'm':    mapfile = *argv+1; SKIPARG; break;
  152.     case 'n':    dither = 0; break;
  153.     case 'd':    debug++; break;
  154.     case 'D':    showcolor++; break;
  155.     case 'v':    verbose++; break;
  156.     case 'A':    outtype = FMT_ATK; break;
  157.     case 'B':    outtype = FMT_FACE; break;
  158.     case 'F':    outtype = FMT_FBM; break;
  159.     case 'G':    outtype = FMT_GIF; break;
  160.     case 'I':    outtype = FMT_IFF; break;
  161.     case 'L':    outtype = FMT_LEAF; break;
  162.     case 'M':    outtype = FMT_MCP; break;
  163.     case 'P':    outtype = FMT_PBM; break;
  164.     case 'R':    outtype = FMT_RLE; break;
  165.     case 'S':    outtype = FMT_SUN; break;
  166.     case 'T':    outtype = FMT_TIFF; break;
  167.     case 'X':    outtype = FMT_X11; break;
  168.     case 'Z':    outtype = FMT_PCX; break;
  169.     default:    fprintf (stderr, "%s\n", USAGE);
  170.             exit (1);
  171.       }
  172.     }
  173.   }
  174.   
  175.   /* Check arguments */
  176.   if (colors > 256 || colors < 4)
  177.   { fprintf (stderr, "fbquant can only handle 4..256 colors\n");
  178.     exit (1);
  179.   }
  180.  
  181.   /* Choose default resbits: if using few colors, assume low-res */
  182.   if (resbits < 0)
  183.   { if (colors <= 32)    resbits = 4;
  184.     else        resbits = 8;
  185.   }
  186.  
  187.   if (resbits > 8 || resbits < 1)
  188.   { fprintf (stderr, "color resolution (-r) must be 1..8 bits\n");
  189.     exit (1);
  190.   }
  191.  
  192.   resmask = (0xff - (0xff >> resbits));
  193.  
  194.   /* Open file if name given as argument */
  195.   if (! read_bitmap (&input, (argc > 0) ? argv[0] : (char *) NULL))
  196.   { exit (1); }
  197.  
  198.   /* Read colormap from map image (if specified) */
  199.   if (mapfile != NULL)
  200.   { if (!read_bitmap (&mapimage, mapfile))
  201.     { perror (mapfile); exit (1); }
  202.     
  203.     if (mapimage.hdr.clrlen == 0)
  204.     { fprintf (stderr, "fbquant: %s: no map\n", mapfile); exit (1); }
  205.     
  206.     /* Dont care about the map image, just its colormap */
  207.     free (mapimage.bm); mapimage.bm = NULL;
  208.     
  209.     colors = mapimage.hdr.clrlen / 3;
  210.     cmap = (COLOR *) malloc ((unsigned) colors * sizeof (COLOR));
  211.     
  212.     for (c=0; c<colors; c++)
  213.     { cmap[c].rd = mapimage.cm[c];
  214.       cmap[c].gr = mapimage.cm[c+colors];
  215.       cmap[c].bl = mapimage.cm[c+colors*2];
  216.       cmap[c].indx = c;
  217.     }
  218.     
  219.     fprintf (stderr, "Quantizing \"%s\" [%dx%d] with %d colors from %s\n",
  220.          input.hdr.title, input.hdr.cols, input.hdr.rows, colors, mapfile);
  221.   }
  222.   else
  223.   { fprintf (stderr,
  224.          "Quantizing \"%s\" [%dx%d] with %d colors, %d bits/clr\n",
  225.          input.hdr.title, input.hdr.cols, input.hdr.rows, colors, resbits);
  226.   }
  227.  
  228.   if (input.hdr.planes != 3 || input.hdr.bits != 8)
  229.   { fprintf (stderr, "fbquant can only handle 24bit RGB inputs\n");
  230.     exit (1);
  231.   }
  232.  
  233.   /* Now build header for output bit map */
  234.   output.hdr = input.hdr;
  235.   output.hdr.planes = 1;
  236.   output.hdr.clrlen = 3 * colors;
  237.   
  238.   /* Allocate space for output */
  239.   alloc_fbm (&output);
  240.  
  241.   /* Build colormap by sampling and mcut if not specified */
  242.   if (mapfile == NULL)
  243.   {
  244.     cmap = (COLOR *) malloc ((unsigned) colors * sizeof (COLOR));
  245.     sample_image (&input, hist);
  246.     build_colormap (hist, cmap, colors);
  247.   }
  248.   
  249.   /* Use Floyd-Steinberg error diffusion to quantize using the new cmap */
  250.   clr_quantize (&input, &output, cmap, colors);
  251.   
  252.   /* Write out the result */
  253.   if (write_bitmap (&output, stdout, outtype))
  254.   { exit (0); }
  255.  
  256.   exit (1);
  257. }
  258.  
  259. /****************************************************************
  260.  * sample_image:
  261.  ****************************************************************/
  262.  
  263. sample_image (image, hist)
  264. FBM *image;
  265. int *hist;
  266. { register int i, j;
  267.   register unsigned char *rp, *gp, *bp;
  268.   int width = image->hdr.cols, height = image->hdr.rows;
  269.   int rowlen = image->hdr.rowlen, plnlen = image->hdr.plnlen;
  270.   int used=0;
  271.   
  272.   /* Clear the histogram */
  273.   for (i=0; i<CUBSIZ; i++) hist[i] = 0;
  274.   
  275.   /* Now count colors */
  276.   for (j=0; j<height; j++)
  277.   { rp = &(image->bm[j*rowlen]);
  278.     gp = rp+plnlen;
  279.     bp = gp+plnlen;
  280.  
  281.     /* If resolution is greater than cube size, just count */
  282.     if (resbits >= CUBITS)
  283.     { for (i=0; i<width; i++)
  284.       { if (++hist[ CLRINDEX8 (*rp++, *gp++, *bp++) ] == 1) used++; }
  285.     }
  286.  
  287.     /* If resolution is less than cube size, ignore extra bits */
  288.     else
  289.     { for (i=0; i<width; i++)
  290.       { if (++hist[ CLRINDEX8 (*rp++ & resmask,
  291.                    *gp++ & resmask,
  292.                    *bp++ & resmask) ] == 1)
  293.     { used++; }
  294.       }
  295.     }
  296.   }
  297.   
  298.   /* Bonus for flesh color */
  299.   if (flesh > 1.01 || flesh < 0.99)
  300.   { register int diff = 0;
  301.     int maxdiff = 0;
  302.     
  303.     maxdiff += (fr > 127) ? fr : 256 - fr;
  304.     maxdiff += (fg > 127) ? fr : 256 - fr;
  305.     maxdiff += (fb > 127) ? fr : 256 - fr;
  306.     
  307.     if (debug)
  308.     { fprintf (stderr, "\nFlesh bonus %1.4lf, max diff %d\n",
  309.            flesh, maxdiff);
  310.     }
  311.  
  312.     for (i=0; i<CUBSIZ; i++)
  313.     { if (hist[i] > 0)
  314.       { j = GETR8(i) - fr; j = abs (j); diff += j;
  315.         j = GETG8(i) - fr; j = abs (j); diff += j;
  316.     j = GETB8(i) - fr; j = abs (j); diff += j;
  317.     
  318.     if (maxdiff < 64)
  319.         { hist[i] *= (flesh - ((flesh * diff) / maxdiff)); }
  320.       }
  321.     }
  322.   }
  323.  
  324.   if (debug)
  325.   { fprintf (stderr, "Done counting, used %d colors for %d pixels\n",
  326.          used, width*height);
  327.   }
  328. }
  329.  
  330. /****************************************************************
  331.  * build_colormap:
  332.  ****************************************************************/
  333.  
  334. # define SWAP(A,B) (t=A,A=B,B=t)
  335.  
  336. build_colormap (hist, cmap, colors)
  337. int *hist;
  338. COLOR *cmap;
  339. int colors;
  340. { register int i, k;
  341.   PIXEL box[CUBSIZ];
  342.   register PIXEL *b;
  343.   int used=0, t;
  344.  
  345.   /* Build the first box, encompassing all pixels */  
  346.   for (b=box, i=0; i<CUBSIZ; i++)
  347.   { b->color = i;
  348.     k = hist[i];
  349.     b->cnt = (k > MAXSHRT) ? MAXSHRT : k;
  350.     b++;
  351.   }
  352.   
  353.   /* Move all non-zero count pixels to the front of the list */
  354.   for (i=0, used=0; i<CUBSIZ; i++)
  355.   { if (box[i].cnt > 0) box[used++] = box[i]; }
  356.  
  357.   /*-------- Special case if we didnt need all colors --------*/
  358.   if (used <= colors)
  359.   {
  360.     /* Copy the colors actually found */
  361.     for (i=0; i<used; i++)
  362.     { cmap[i].rd = GETR8 (box[i].color);
  363.       cmap[i].gr = GETG8 (box[i].color);
  364.       cmap[i].bl = GETB8 (box[i].color);
  365.     }
  366.  
  367.     /* Set the rest to WHITE */
  368.     for (; i<colors; i++)
  369.     { cmap[i].rd = 255;
  370.       cmap[i].gr = 255;
  371.       cmap[i].bl = 255;
  372.     }
  373.   }
  374.   
  375.   else                /* Build sets all colors */
  376.   { split_box (box, used, 0, colors, cmap); }
  377.   
  378.   /*----------------------------------------------------------------
  379.    * Now arrange the colors in the desired order.  Sun convention is that
  380.    * color 0 is white and color n-1 is black, so we sort by increasing
  381.    * intensity (why not?) and then swap the lightest and darkest colors
  382.    */
  383.  
  384.   /* Now sort 0..n-1 according to increasing intensity */
  385.   qsort (cmap, colors, sizeof (* cmap), cmp_int);
  386.  
  387.   /* Make first color lightest and last color darkest */
  388.   SWAP (cmap[0].rd, cmap[colors-1].rd);
  389.   SWAP (cmap[0].gr, cmap[colors-1].gr);
  390.   SWAP (cmap[0].bl, cmap[colors-1].bl);
  391.   
  392.   /* Set the output indices */
  393.   for (i=0; i<colors; i++) { cmap[i].indx = i; }
  394. }
  395.  
  396. /****************************************************************
  397.  * split_box: Basic recursive part of Heckberts adaptive partitioning
  398.  *          algorithm.
  399.  ****************************************************************/
  400.  
  401. split_box (box, boxlen, clr, numclr, cmap)
  402. PIXEL *box;
  403. int boxlen, clr, numclr;
  404. COLOR *cmap;
  405. { int maxv[3], minv[3], numv[3];
  406.   int pcnt[3][CUBSID];
  407.   int sbox, snum, split, half, maxdif, dif;
  408.   register PIXEL *top, *bot;
  409.   int topw, botw;
  410.   int red, grn, blu;
  411.   register int i, c;
  412.  
  413.   /* If numclr exceeds boxlen, we are in trouble */
  414.   if (numclr > boxlen)
  415.   { fprintf (stderr, "boxlen %d is less numclr %d, panic!\n than",
  416.          boxlen, numclr);
  417.     fflush (stderr);
  418.     abort ();
  419.   }
  420.  
  421.   /* Base case: only one color, assign the average for this cell */
  422.   if (numclr <= 1)
  423.   { red = box_avg_red (box, boxlen);
  424.     grn = box_avg_grn (box, boxlen);
  425.     blu = box_avg_blu (box, boxlen);
  426.     
  427.     /* Map x to x+4, because the histogram maps values to multiples of 8 */
  428.     cmap[clr].rd = red + 4;
  429.     cmap[clr].gr = grn + 4;
  430.     cmap[clr].bl = blu + 4;
  431.     
  432.     /* Mask off unused bits for color resolution */
  433.     cmap[clr].rd &= resmask;
  434.     cmap[clr].gr &= resmask;
  435.     cmap[clr].bl &= resmask;
  436.     
  437.     if (debug)
  438.     { fprintf (stderr, "\t\tassigning color %d  <%d,%d,%d>  (%d)\n",
  439.            clr, cmap[clr].rd, cmap[clr].gr, cmap[clr].bl,
  440.            box_weight (box, boxlen));
  441.     }
  442.     
  443.     return;
  444.   }
  445.  
  446.   /* Gather statistics about the boxes contents */
  447.   minv[RD] = minv[GR] = minv[BL] = CUBSID;
  448.   maxv[RD] = maxv[GR] = maxv[BL] = 0;
  449.   numv[RD] = numv[GR] = numv[BL] = 0;
  450.   for (i=0; i<CUBSID; i++) { pcnt[RD][i] = pcnt[GR][i] = pcnt[BL][i] = 0; }
  451.   
  452.   for (i=0; i<boxlen; i++)
  453.   { c = box[i].color;
  454.     red = GETR(c); grn = GETG(c); blu = GETB(c);
  455.     
  456.     if (red < minv[RD]) minv[RD] = red;
  457.     if (red > maxv[RD]) maxv[RD] = red;
  458.     if (grn < minv[GR]) minv[GR] = grn;
  459.     if (grn > maxv[GR]) maxv[GR] = grn;
  460.     if (blu < minv[BL]) minv[BL] = blu;
  461.     if (blu > maxv[BL]) maxv[BL] = blu;
  462.     
  463.     if (++pcnt[RD][red] == 1) numv[RD]++;
  464.     if (++pcnt[GR][grn] == 1) numv[GR]++;
  465.     if (++pcnt[BL][blu] == 1) numv[BL]++;
  466.   }
  467.  
  468.   /* Special case, boxlen = numclr, just assign each box one color */
  469.   if (boxlen == numclr)
  470.   { for (i=0; i<boxlen; i++)
  471.     { split_box (box+i, 1, clr+i, 1, cmap); }
  472.     return;
  473.   }
  474.  
  475.   /*
  476.    * Pick a dimension to split.  Currently splits longest dimension
  477.    * in RGB space.  Heckbert and other experts suggest splitting
  478.    * dimension with greatest variance may work better.
  479.    */
  480.  
  481.   split = -1; maxdif = -1;
  482.  
  483.   if ((dif = (maxv[RD] - minv[RD])) > maxdif) { maxdif = dif; split = RD; }
  484.   if ((dif = (maxv[GR] - minv[GR])) > maxdif) { maxdif = dif; split = GR; }
  485.   if ((dif = (maxv[BL] - minv[BL])) > maxdif) { maxdif = dif; split = BL; }
  486.   
  487.   /* Sort along the chosen dimension */
  488.   switch (split)
  489.   { case RD:    qsort (box, boxlen, sizeof (*box), cmp_red); break;
  490.     case GR:    qsort (box, boxlen, sizeof (*box), cmp_grn); break;
  491.     case BL:    qsort (box, boxlen, sizeof (*box), cmp_blu); break;
  492.     default:    fprintf (stderr, "panic in split_box, split %d\n", split);
  493.         fflush (stderr); fflush (stdout); abort ();
  494.   }
  495.   
  496.   /* 
  497.    * Split at the median, but make sure there are at least numclr/2
  498.    * different colors on each side of the split, to avoid wasting
  499.    * colors.
  500.    *
  501.    * Note: need to keep in mind that when the box is large, dont split
  502.    *       too close to one edge.
  503.    */
  504.    
  505.   half = numclr / 2;
  506.   top = box;        bot = box + (boxlen-1);
  507.   topw = top->cnt;    botw = bot->cnt;
  508.   
  509.   /* Set top and bot to point to min/max feasible splits */
  510.   while ((top-box)+1 < half)        { top++; topw += top->cnt; }
  511.   while ((boxlen-(bot-box)) < half)    { bot--; botw += bot->cnt; }
  512.  
  513.   /* Move top and bottom towards each other 1/8 of remaining distance */
  514.   c = (bot-top) / 8;
  515.   for (i=0; i<c; i++)            { top++; topw += top->cnt; }
  516.   for (i=0; i<c; i++)            { bot--; botw += bot->cnt; }
  517.  
  518.   /* Now search for median */
  519.   while (top < bot)
  520.   { if (topw < botw)            { top++; topw += top->cnt; }
  521.     else                { bot--; botw += bot->cnt; }
  522.   }
  523.  
  524.   /* Decide which half gets the midpoint */
  525.   if (topw > botw)            /* Median wants to go with top */
  526.   { sbox = (top-box) + 1;
  527.     snum = numclr - half;
  528.   }
  529.   else                    /* Median wants to go with bottom */
  530.   { sbox = (top - box);
  531.     snum = half;
  532.   }
  533.   
  534.   /* Handle boundary conditions with number of colors vs box size */
  535.   if (sbox == 0) sbox++;
  536.   else if (sbox == boxlen) sbox--;
  537.  
  538.   while (snum > sbox) snum--;
  539.   while (numclr-snum > boxlen-sbox) snum++;
  540.  
  541. # ifndef OPTIMIZE
  542.   /* Check for boundary condition errors */
  543.   if (snum <= 0 || snum >= numclr)
  544.   { fprintf (stderr, "panic, using zero colors for box\n");
  545.     fflush (stderr);
  546.     abort ();
  547.   }
  548.  
  549.   if (boxlen-sbox < numclr-snum)
  550.   { fprintf (stderr, "panic, about to used %d boxes for %d colors\n",
  551.          boxlen-sbox, numclr-snum);
  552.     fflush (stderr);
  553.     abort ();
  554.   }
  555.  
  556.   if (sbox < snum)
  557.   { fprintf (stderr, "panic, about to used %d boxes for %d colors\n",
  558.          sbox, snum);
  559.     fflush (stderr);
  560.     abort ();
  561.   }
  562. # endif
  563.  
  564.   if (debug)
  565.   { int count = numclr, depth = 8;
  566.     while (count > 0) { depth--; count /= 2; }
  567.     for (i=0; i<depth; i++) fprintf (stderr, "  ");
  568.     
  569.     fprintf (stderr, "box [%d..%d|%d] r(%d,%d,%d) g(%d,%d,%d) b(%d,%d,%d) =>",
  570.          0, boxlen-1, numclr,
  571.          minv[RD], maxv[RD], numv[RD],
  572.          minv[GR], maxv[GR], numv[GR],
  573.          minv[BL], maxv[BL], numv[BL]);
  574.     fprintf (stderr, " %c [%d..%d|%d] [%d..%d|%d]\n",
  575.          "RGB"[split], 0, sbox-1, snum, sbox, boxlen-1, numclr-snum);
  576.   }
  577.  
  578.   /* Now recurse and split each sub-box */
  579.   split_box (box,      sbox,          clr,      snum,        cmap);
  580.   split_box (box+sbox, boxlen - sbox, clr+snum, numclr-snum, cmap);
  581. }
  582.  
  583. /****************************************************************
  584.  * box_weight: Determine the total count of a box
  585.  ****************************************************************/
  586.  
  587. box_weight (box, boxlen)
  588. register PIXEL *box;
  589. register int boxlen;
  590. { register int sum = 0, i;
  591.  
  592.   for (i=0; i<boxlen; i++)
  593.   { sum += box[i].cnt; }
  594.   
  595.   return (sum);
  596. }
  597.  
  598. /****************************************************************
  599.  * box_avg_red: Determine the average red value [0..255] of a box
  600.  ****************************************************************/
  601.  
  602. box_avg_red (box, boxlen)
  603. register PIXEL *box;
  604. register int boxlen;
  605. { register int sum = 0, n = 0, i;
  606.  
  607.   for (i=0; i<boxlen; i++)
  608.   { sum += GETR8(box[i].color); n++; }
  609.   
  610.   return (sum / n);
  611. }
  612.  
  613. /****************************************************************
  614.  * box_avg_grn: Determine the average grn value [0..255] of a box
  615.  ****************************************************************/
  616.  
  617. box_avg_grn (box, boxlen)
  618. register PIXEL *box;
  619. register int boxlen;
  620. { register int sum = 0, n = 0, i;
  621.  
  622.   for (i=0; i<boxlen; i++)
  623.   { sum += GETG8(box[i].color); n++; }
  624.   
  625.   return (sum / n);
  626. }
  627.  
  628. /****************************************************************
  629.  * box_avg_blu: Determine the average blu value [0..255] of a box
  630.  ****************************************************************/
  631.  
  632. box_avg_blu (box, boxlen)
  633. register PIXEL *box;
  634. register int boxlen;
  635. { register int sum = 0, n = 0, i;
  636.  
  637.   for (i=0; i<boxlen; i++)
  638.   { sum += GETB8(box[i].color); n++; }
  639.   
  640.   return (sum / n);
  641. }
  642.  
  643.  
  644. /****************************************************************
  645.  * sort by increasing red ( then green and blue )
  646.  ****************************************************************/
  647.  
  648. cmp_red (a, b)
  649. PIXEL *a, *b;
  650. { register r;
  651.  
  652.   if (r = GETR(a->color) - GETR(b->color))
  653.   { return (r); }
  654.   
  655.   if (r = GETG(a->color) - GETG(b->color))
  656.   { return (r); }
  657.  
  658.   if (r = GETB(a->color) - GETB(b->color))
  659.   { return (r); }
  660.   
  661.   return (0);
  662. }
  663.  
  664. /****************************************************************
  665.  * sort by increasing green ( then blue and red )
  666.  ****************************************************************/
  667.  
  668. cmp_grn (a, b)
  669. PIXEL *a, *b;
  670. { register r;
  671.  
  672.   if (r = GETG(a->color) - GETG(b->color))
  673.   { return (r); }
  674.  
  675.   if (r = GETB(a->color) - GETB(b->color))
  676.   { return (r); }
  677.   
  678.   if (r = GETR(a->color) - GETR(b->color))
  679.   { return (r); }
  680.   
  681.   return (0);
  682. }
  683.  
  684. /****************************************************************
  685.  * sort by increasing blue ( then red and green )
  686.  ****************************************************************/
  687.  
  688. cmp_blu (a, b)
  689. PIXEL *a, *b;
  690. { register r;
  691.  
  692.   if (r = GETB(a->color) - GETB(b->color))
  693.   { return (r); }
  694.   
  695.   if (r = GETR(a->color) - GETR(b->color))
  696.   { return (r); }
  697.   
  698.   if (r = GETG(a->color) - GETG(b->color))
  699.   { return (r); }
  700.  
  701.   return (0);
  702. }
  703.  
  704. /****************************************************************
  705.  * clr_quantize: Do Floyd Steinberg quantizing on the image
  706.  ****************************************************************/
  707.  
  708. clr_quantize (input, output, cmap, colors)
  709. FBM *input, *output;
  710. COLOR *cmap;
  711. int colors;
  712. { int **cerr, **lerr, **terr;
  713.   int width = input->hdr.cols, height = input->hdr.rows;
  714.   int rowlen = input->hdr.rowlen, plnlen = input->hdr.plnlen;
  715.   register int i, j;
  716.   register int rd, gr, bl;
  717.   int rderr, grerr, blerr, clr;
  718.   unsigned char *rp, *gp, *bp, *obm;
  719.  
  720.   /*-------- Copy colormap into output bitmap --------*/
  721.   for (i=0; i<colors; i++)
  722.   { output->cm[i]               = cmap[i].rd;
  723.     output->cm[i+colors]        = cmap[i].gr;
  724.     output->cm[i+colors+colors] = cmap[i].bl;
  725.   }
  726.  
  727.   /*
  728.    * Precompute necessary nearest neighbor query info.  Note that we wait
  729.    * until after copying the colormap, since init reorders it
  730.    */
  731.  
  732.   init_nearest (cmap, colors);
  733.  
  734.   /*-------- Do halftoning --------*/
  735.   cerr = (int **) malloc (3 * sizeof (int *));
  736.   lerr = (int **) malloc (3 * sizeof (int *));
  737.   cerr[RD] = (int *) malloc (width * sizeof (int));
  738.   cerr[GR] = (int *) malloc (width * sizeof (int));
  739.   cerr[BL] = (int *) malloc (width * sizeof (int));
  740.   lerr[RD] = (int *) malloc (width * sizeof (int));
  741.   lerr[GR] = (int *) malloc (width * sizeof (int));
  742.   lerr[BL] = (int *) malloc (width * sizeof (int));
  743.   
  744.   /*-------- Just use the nearest color everywhere --------*/
  745.   if (!dither)
  746.   {
  747.     for (j=0; j<height; j++)
  748.     { rp = &(input->bm[j*rowlen]);
  749.       gp = rp+plnlen;
  750.       bp = gp+plnlen;
  751.       obm = &(output->bm[j*rowlen]);
  752.  
  753.       for (i=0; i<width; i++)
  754.       { rd = *rp++; gr = *gp++; bl = *bp++;
  755.  
  756.         obm[i] = cmap[nearest (rd, gr, bl, cmap, colors)].indx;
  757.       }
  758.     }
  759.     
  760.     return;
  761.   }
  762.  
  763.   /*------ Just use the nearest color around the left, right, and top ------*/
  764.  
  765.   /* Top border */
  766.   rp = input->bm; gp = rp+plnlen; bp = gp+plnlen;
  767.   obm = output->bm;
  768.   for (i=0; i<width; i++)
  769.   { obm[i] = cmap[nearest (rp[i], gp[i], bp[i], cmap, colors)].indx; }
  770.  
  771.   /* Left border */
  772.   rp = input->bm; gp = rp+plnlen; bp = gp+plnlen;
  773.   obm = output->bm;
  774.   for (j=0; j<height; j++)
  775.   { obm[j*rowlen] = cmap[nearest (rp[j*rowlen], gp[j*rowlen],
  776.              bp[j*rowlen], cmap, colors)].indx; }
  777.  
  778.   /* Right border */
  779.   rp = &(input->bm[width-1]); gp = rp + plnlen; bp = gp + plnlen;
  780.   obm = &(output->bm[width-1]);
  781.   for (j=0; j<height; j++)
  782.   { obm[j*rowlen] = cmap[nearest (rp[j*rowlen], gp[j*rowlen],
  783.                   bp[j*rowlen], cmap, colors)].indx; }
  784.  
  785.   /*-------- Clear error vectors --------*/
  786.   for (i=0; i<width; i++)
  787.   { cerr[RD][i] = cerr[GR][i] = cerr[BL][i] = 0;
  788.     lerr[RD][i] = lerr[GR][i] = lerr[BL][i] = 0;
  789.   }
  790.  
  791.   /*-------- Major loop for Floyd-Steinberg --------*/
  792.   for (j=1; j<height; j++)
  793.   { rp = &(input->bm[j*rowlen+1]);
  794.     gp = rp+plnlen;
  795.     bp = gp+plnlen;
  796.     obm = &(output->bm[j*rowlen+1]);
  797.  
  798.     for (i=1; i<width-1; i++)
  799.     { rd = *rp++; gr = *gp++; bl = *bp++;
  800.  
  801.       /* Sum up errors using Floyd-Steinberg weights */
  802.       rderr= cerr[RD][i-1] + 7*lerr[RD][i-1] + 5*lerr[RD][i] + 3*lerr[RD][i+1];
  803.       grerr= cerr[GR][i-1] + 7*lerr[GR][i-1] + 5*lerr[GR][i] + 3*lerr[GR][i+1];
  804.       blerr= cerr[BL][i-1] + 7*lerr[BL][i-1] + 5*lerr[BL][i] + 3*lerr[BL][i+1];
  805.  
  806.       rderr >>= 4;    /* Divide by 16 */
  807.       grerr >>= 4;    /* Divide by 16 */
  808.       blerr >>= 4;    /* Divide by 16 */
  809.  
  810.       /* Choose nearest color to adjusted RGB value */
  811.       rd += rderr; gr += grerr; bl += blerr;
  812.  
  813.       clr = nearest (rd, gr, bl, cmap, colors);
  814.       *obm++ = cmap[clr].indx;
  815.  
  816.       /* Compute accumulated error for this pixel */      
  817.       rd -= (int) cmap[clr].rd;
  818.       gr -= (int) cmap[clr].gr;
  819.       bl -= (int) cmap[clr].bl;
  820.       
  821.       /* Limit error (avoids color shadows) to 75% if over MAXERR */
  822.       if (rd < -MAXERR || rd > MAXERR)  rd = (rd * 3) >> 2;
  823.       if (gr < -MAXERR || gr > MAXERR)  gr = (gr * 3) >> 2;
  824.       if (bl < -MAXERR || bl > MAXERR)  bl = (bl * 3) >> 2;
  825.  
  826.       /* Store errors in error vectors */
  827.       cerr[RD][i] = rd;
  828.       cerr[GR][i] = gr;
  829.       cerr[BL][i] = bl;
  830.     }
  831.     
  832.     /* Swap error vectors */
  833.     terr = lerr; lerr = cerr; cerr = terr;
  834.   }
  835. }
  836.  
  837. /****************************************************************
  838.  * nearest: Choose nearest color
  839.  ****************************************************************/
  840.  
  841. short cache[CUBSIZ];
  842.  
  843. init_nearest (cmap, colors)
  844. COLOR *cmap;
  845. int colors;
  846. { register int i;
  847.  
  848.   /* Initialize the cache */
  849.   for (i=0; i<CUBSIZ; i++) { cache[i] = -1; }
  850.   
  851.   /* Sort the colormap by decreasing red, then green, then blue */
  852.   qsort (cmap, colors, sizeof (COLOR), cmp_cmap);
  853.  
  854.   if (debug || showcolor)
  855.   { fprintf (stderr, "\nFinal colormap:\n\n");
  856.     for (i=0; i<colors; i++)
  857.     { fprintf (stderr, "%3d:  <%3d,%3d,%3d> [%d]\n", 
  858.            i, cmap[i].rd, cmap[i].gr, cmap[i].bl, cmap[i].indx);
  859.     }
  860.   }
  861. }
  862.  
  863. /* Fast square macro, uses local variable to avoid mulitple eval of X */
  864. # define SQR(X) (sqrtmp = (X), sqrtmp*sqrtmp)
  865.  
  866. /* Fast distance macro */
  867. # define CDIST(CPTR,CR,CG,CB)        \
  868. (sumtmp =  SQR (((int) ((CPTR)->rd)) - CR),        \
  869.  sumtmp += SQR (((int) ((CPTR)->gr)) - CG),        \
  870.  sumtmp += SQR (((int) ((CPTR)->bl)) - CB),        \
  871.  sumtmp)
  872.  
  873. # define restrict(X) ((X) < 0 ? 0 : (X) > 255 ? 255 : (X))
  874.  
  875. nearest (rd, gr, bl, cmap, colors)
  876. int rd, gr, bl;
  877. COLOR *cmap;
  878. int colors;
  879. { register int clr, sqrtmp, sumtmp;
  880.   register COLOR *a, *b, *c, *best, *ctail;
  881.   int cindx, bestd, dist;
  882.  
  883.   rd = restrict (rd);
  884.   gr = restrict (gr);
  885.   bl = restrict (bl);
  886.  
  887.   /* Find array index in cache */
  888.   cindx = CLRINDEX8 (rd, gr, bl);
  889.  
  890.   /* Check cache for color value */  
  891.   if ((clr = cache[cindx]) >= 0)
  892.   { return (clr); }
  893.   
  894.   /*
  895.    * Search through colormap for nearest neighbor:
  896.    * 1) Find nearest red value by binary search
  897.    * 2) Search up and down, keeping best point
  898.    * 3) Terminate search when red difference is greater than best pt
  899.    */
  900.   
  901.   /* Binary search for nearest red value */
  902.   ctail = &cmap[colors];
  903.   for (a=cmap, b= ctail-1;  a<b;  )  
  904.   { c = a + (b-a)/2;    /* Find midpoint */
  905.  
  906.     if (debug && verbose)
  907.     { fprintf (stderr, "Searching for %d, a=%d (%d), b=%d (%d), c=%d (%d)\n",
  908.         rd, a-cmap, a->rd, b-cmap, b->rd, c-cmap, c->rd);
  909.     }
  910.  
  911.     if (c->rd == rd)
  912.     { break; }
  913.     else if (c->rd < rd)
  914.     { a = ++c; }
  915.     else
  916.     { b = c; }
  917.   }
  918.  
  919.   /*
  920.    * c now points to closest red value, search up and down for closer
  921.    * Euclidean distance, and stop search when the red distance alone
  922.    * exceeds the best found.
  923.    */
  924.  
  925.   /* Set best and bestd to best red value */
  926.   best = c;
  927.   bestd = CDIST (c, rd, gr, bl);
  928.  
  929.   if (debug && verbose)
  930.     fprintf (stderr, "Found c=%d (%d)  dist = %d\n", c-cmap, c->rd, bestd);
  931.  
  932.   /* a and b are search pointers up and down */
  933.   a = c-1;
  934.   b = c+1;
  935.  
  936.   while (bestd > 0 && (a >= cmap || b < ctail))
  937.   {
  938.     if (debug && verbose)
  939.     { fprintf (stderr, "  search: bestd %d, best %d|%d, a %d, b %d\n",
  940.         bestd, best-cmap, best->indx, a-cmap, b-cmap);
  941.     }
  942.  
  943.     if (a >= cmap)
  944.     { if ((dist = CDIST (a, rd, gr, bl)) < bestd)
  945.       { best = a; bestd = dist; }
  946.       
  947.       if (SQR ((int) a->rd - rd) >= bestd) a = cmap-1;
  948.       else                   a--;
  949.     }
  950.  
  951.     if (b < ctail)
  952.     { if ((dist = CDIST (b, rd, gr, bl)) < bestd)
  953.       { best = b; bestd = dist; }
  954.       
  955.       if (SQR ((int) b->rd - rd) >= bestd) b = ctail;
  956.       else                   b++;
  957.     }
  958.   }
  959.   
  960.   if (best < cmap || best >= ctail)
  961.   { perror ("returning invalid color\n");
  962.     abort ();
  963.   }
  964.  
  965.   clr = (best - cmap);
  966.   
  967.   if (debug && verbose)
  968.   { fprintf (stderr, "Nearest(%3d,%3d,%3d) => <%3d,%3d,%3d> [%d]\n",
  969.         rd, gr, bl, best->rd, best->gr, best->bl, best->indx);
  970.   }
  971.  
  972.   return ((cache[cindx] = clr));  
  973. }
  974.  
  975. /****************************************************************
  976.  * sort colormap by increasing red ( then green and blue )
  977.  ****************************************************************/
  978.  
  979. cmp_cmap (a, b)
  980. register COLOR *a, *b;
  981. { register int r;
  982.  
  983.   if (r = (a->rd - b->rd)) { return (r); }
  984.   if (r = (a->gr - b->gr)) { return (r); }
  985.   if (r = (a->bl - b->bl)) { return (r); }
  986.   
  987.   return (0);
  988. }
  989.  
  990. /****************************************************************
  991.  * sort colormap by increasing intensity (Use NTSC weights)
  992.  ****************************************************************/
  993.  
  994. # define RW 299
  995. # define GW 587
  996. # define BW 114
  997.  
  998. cmp_int (a, b)
  999. register COLOR *a, *b;
  1000. { register int ai, bi;
  1001.  
  1002.   ai = RW * a->rd + GW * a->gr + BW * a->bl;
  1003.   bi = RW * b->rd + GW * b->gr + BW * b->bl;
  1004.  
  1005.   return (ai - bi);
  1006. }
  1007.