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