home *** CD-ROM | disk | FTP | other *** search
/ Resource Library: Multimedia / Resource Library: Multimedia.iso / maestro / source / displytl / xv24to8.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-06-15  |  42.0 KB  |  1,719 lines
  1. /*
  2.  * xv24to8.c  -  contains the 24-to-8-bit Conv24to8 procedure
  3.  *
  4.  * The Conv24to8 procedure takes a pointer to a 24-bit image (loaded
  5.  * previously).  The image will be a w * h * 3 byte array of
  6.  * bytes.  The image will be arranged with 3 bytes per pixel (in order
  7.  * R, G, and B), pixel 0 at the top left corner.  (As normal.)
  8.  * The procedure also takes a maximum number of colors to use (numcols)
  9.  *
  10.  * The Conv24to8 procedure will set up the following:  it will allocate and
  11.  * create 'pic', a pWIDE*pHIGH 8-bit picture.  it will set up pWIDE and pHIGH.
  12.  * it will load up the r[], g[], and b[] colormap arrays.  it will NOT 
  13.  * calculate numcols, since the cmap sort procedure has to be called anyway
  14.  *
  15.  * Conv24to8 returns '0' if successful, '1' if it failed (presumably on a 
  16.  * malloc())
  17.  *
  18.  * The 'slow' code is based on Heckbert's Median Cut algorithm.
  19.  *
  20.  * contains:
  21.  *   Cont24to8()
  22.  *   Init24to8()
  23.  */
  24.  
  25. /*
  26.  * Copyright 1989, 1990, 1991, 1992 by John Bradley and
  27.  *                       The University of Pennsylvania
  28.  *
  29.  * Permission to use, copy, and distribute for non-commercial purposes,
  30.  * is hereby granted without fee, providing that the above copyright
  31.  * notice appear in all copies and that both the copyright notice and this
  32.  * permission notice appear in supporting documentation.
  33.  *
  34.  * The software may be modified for your own purposes, but modified versions
  35.  * may not be distributed.
  36.  *
  37.  * This software is provided "as is" without any expressed or implied warranty.
  38.  *
  39.  * The author may be contacted via:
  40.  *    US Mail:   John Bradley
  41.  *               GRASP Lab, Room 301C
  42.  *               3401 Walnut St.
  43.  *               Philadelphia, PA  19104
  44.  *
  45.  *    Phone:     (215) 898-8813
  46.  *    EMail:     bradley@cis.upenn.edu
  47.  */
  48.  
  49.  
  50. /*
  51.  * Portions Copyright (C) 1989, 1991 by Jef Poskanzer.  See copyright notice
  52.  * at beginning of relevant code.
  53.  */
  54.  
  55.  
  56. #include "xvimage.h"
  57.  
  58. #define    MAX_CMAP_SIZE    256
  59. #define    COLOR_DEPTH    8
  60. #define    MAX_COLOR    256
  61. #define    B_DEPTH        5        /* # bits/pixel to use */
  62. #define    B_LEN        (1<<B_DEPTH)
  63. #define    C_DEPTH        3
  64. #define    C_LEN        (1<<C_DEPTH)    /* # cells/color to use */
  65.  
  66. #define R2FACT 20 /* .300 * .300 * 256 = 23 */
  67. #define G2FACT 39 /* .586 * .586 * 256 = 88 */
  68. #define B2FACT 8  /* .114 * .114 * 256 =  3 */
  69.  
  70.  
  71. typedef    struct colorbox {
  72.   struct colorbox *next, *prev;
  73.   int              rmin,rmax, gmin,gmax, bmin,bmax;
  74.   int              total;
  75. } CBOX;
  76.  
  77. typedef struct {
  78.   int num_ents;
  79.   int entries[MAX_CMAP_SIZE][2];
  80. } CCELL;
  81.  
  82. static byte *pic24;
  83. static int   num_colors, WIDE, HIGH;
  84. static int   histogram[B_LEN][B_LEN][B_LEN];
  85.  
  86. CBOX   *freeboxes, *usedboxes;
  87. CCELL **ColorCells;
  88.  
  89. static void   get_histogram();
  90. static CBOX  *largest_box();
  91. static void   splitbox();
  92. static void   shrinkbox();
  93. static void   assign_color();
  94. static CCELL *create_colorcell();
  95. static void   map_colortable();
  96. static int    quant_fsdither();
  97. static int    Quick24to8();
  98. static int    QuickCheck();
  99. static int    ppmquant();
  100.  
  101. static int    tbl1[512],     /* tables used in F-S Dithering */
  102.               tbl3[512],     /* contain i/16, 3i/16, 5i/16, 7i/16, */
  103.               tbl5[512],     /* (i=-256..255) respectively */
  104.               tbl7[512];
  105.  
  106.  
  107. /****************************/
  108. void Init24to8()
  109. /****************************/
  110. {
  111.   /* initialize Floyd-Steinberg division tables */
  112.   /* make sure rounding is done correctly for negative values! */
  113.  
  114.   int i;
  115.  
  116.   for (i = -256; i < 0; i++) {
  117.     tbl1[i+256] = -((8 -i  )/16);
  118.     tbl3[i+256] = -((8 -3*i)/16);
  119.     tbl5[i+256] = -((8 -5*i)/16);
  120.     tbl7[i+256] = -((8 -7*i)/16);
  121.   }
  122.  
  123.   for (i = 0; i < 256; i++) {
  124.     tbl1[i+256] = (i  +8)/16;
  125.     tbl3[i+256] = (3*i+8)/16;
  126.     tbl5[i+256] = (5*i+8)/16;
  127.     tbl7[i+256] = (7*i+8)/16;
  128.   }
  129.  
  130. #ifdef FOO 
  131.   for (i=0; i<512; i++) {
  132.     printf("%3d:  tbl1=%3d, tbl3=%3d, tbl5=%3d, tbl7=%3d\n",
  133.        i-256, tbl1[i], tbl3[i], tbl5[i], tbl7[i]);
  134.   }
  135. #endif
  136. }
  137.  
  138.  
  139.  
  140. /****************************/
  141. int Conv24to8(p,w,h,nc)
  142. /****************************/
  143. byte *p;
  144. int   w,h,nc;
  145. {
  146.   int   i, j;
  147.   CBOX *box_list, *ptr;
  148.  
  149.  
  150.   if (nc<=0) nc = 255;  /* 'nc == 0' breaks code */
  151.  
  152.   /* copy arguments to local-global variables */
  153.   pic24 = p;  pWIDE = WIDE = w;  pHIGH = HIGH = h;  num_colors = nc;
  154.  
  155.   /* allocate pic immediately, so that if we can't allocate it, we don't
  156.      waste time running this algorithm */
  157.  
  158.   pic = (byte *) malloc(WIDE * HIGH);
  159.   if (pic == NULL) {
  160.     fprintf(stderr,"%s: Conv24to8() - failed to allocate 'pic'\n",cmd);
  161.     return(1);
  162.   }
  163.  
  164.   if (!noqcheck && QuickCheck(pic24,w,h,nc)) { 
  165.     /* see if it's a <256 color RGB pic */
  166.     printf("No color compression was necessary.\n");
  167.     return 0;   
  168.   }
  169.  
  170.   else if (conv24 == CONV24_FAST) {
  171.     printf("Doing 'quick' 24-bit to 8-bit conversion.\n");
  172.     return(Quick24to8(pic24,w,h));
  173.   }
  174.  
  175.   else if (conv24 == CONV24_BEST) {
  176.     printf("Doing 'best' 24-bit to 8-bit conversion.\n");
  177.     return(ppmquant(pic24,w,h,nc));
  178.   }
  179.  
  180.   else 
  181.     printf("Doing 'slow' 24-bit to 8-bit conversion.\n");
  182.  
  183.  
  184.  
  185.   /**** STEP 1:  create empty boxes ****/
  186.  
  187.   usedboxes = NULL;
  188.   box_list = freeboxes = (CBOX *) malloc(num_colors * sizeof(CBOX));
  189.  
  190.   if (box_list == NULL) {
  191.     fprintf(stderr,"%s: Conv24to8() - failed to allocate 'freeboxes'\n",cmd);
  192.     return(1);
  193.   }
  194.  
  195.   for (i=0; i<num_colors; i++) {
  196.     freeboxes[i].next = &freeboxes[i+1];
  197.     freeboxes[i].prev = &freeboxes[i-1];
  198.   }
  199.   freeboxes[0].prev = NULL;
  200.   freeboxes[num_colors-1].next = NULL;
  201.  
  202.  
  203.  
  204.  
  205.   /**** STEP 2: get histogram, initialize first box ****/
  206.  
  207.   ptr = freeboxes;
  208.   freeboxes = ptr->next;
  209.   if (freeboxes) freeboxes->prev = NULL;
  210.  
  211.   ptr->next = usedboxes;
  212.   usedboxes = ptr;
  213.   if (ptr->next) ptr->next->prev = ptr;
  214.     
  215.   get_histogram(ptr);
  216.  
  217.  
  218.  
  219.  
  220.   /**** STEP 3: continually subdivide boxes until no more free boxes remain */
  221.  
  222.   while (freeboxes) {
  223.     ptr = largest_box();
  224.     if (ptr) splitbox(ptr);
  225.     else break;
  226.   }
  227.  
  228.  
  229.  
  230.   
  231.   /**** STEP 4: assign colors to all boxes ****/
  232.  
  233.   for (i=0, ptr=usedboxes; i<num_colors && ptr; i++, ptr=ptr->next) {
  234.     assign_color(ptr, &r[i], &g[i], &b[i]);
  235.   }
  236.  
  237.   /* We're done with the boxes now */
  238.   num_colors = i;
  239.   free(box_list);
  240.   box_list = freeboxes = usedboxes = NULL;
  241.  
  242.  
  243.  
  244.  
  245.   /**** STEP 5: scan histogram and map all values to closest color */
  246.  
  247.   /* 5a: create cell list as described in Heckbert[2] */
  248.  
  249.   ColorCells = (CCELL **) calloc(C_LEN*C_LEN*C_LEN, sizeof(CCELL *));
  250.  
  251.  
  252.   /* 5b: create mapping from truncated pixel space to color table entries */
  253.  
  254.   map_colortable();
  255.  
  256.  
  257.  
  258.  
  259.   /**** STEP 6: scan image, match input values to table entries */
  260.  
  261.   i = quant_fsdither();
  262.  
  263.   /* free everything that can be freed */
  264.   for (j=0 ; j < C_LEN*C_LEN*C_LEN ; j++) {
  265.     if (ColorCells[j] != NULL) free (ColorCells[j]);
  266.   }
  267.   free(ColorCells);
  268.  
  269.   return i;
  270. }
  271.  
  272.  
  273. /****************************/
  274. static void get_histogram(box)
  275.      CBOX *box;
  276. /****************************/
  277. {
  278.   int   i,j,r,g,b,*ptr;
  279.   byte *p;
  280.  
  281.   box->rmin = box->gmin = box->bmin = 999;
  282.   box->rmax = box->gmax = box->bmax = -1;
  283.   box->total = WIDE * HIGH;
  284.  
  285.   /* zero out histogram */
  286.   ptr = &histogram[0][0][0];
  287.   for (i=B_LEN*B_LEN*B_LEN; i>0; i-- )  *ptr++ = 0;
  288.  
  289.   /* calculate histogram */
  290.   p = pic24;
  291.   for (i=0; i<HIGH; i++)
  292.     for (j=0; j<WIDE; j++) {
  293.       r = (*p++) >> (COLOR_DEPTH-B_DEPTH);
  294.       g = (*p++) >> (COLOR_DEPTH-B_DEPTH);
  295.       b = (*p++) >> (COLOR_DEPTH-B_DEPTH);
  296.  
  297.       if (r < box->rmin) box->rmin = r;
  298.       if (r > box->rmax) box->rmax = r;
  299.  
  300.       if (g < box->gmin) box->gmin = g;
  301.       if (g > box->gmax) box->gmax = g;
  302.  
  303.       if (b < box->bmin) box->bmin = b;
  304.       if (b > box->bmax) box->bmax = b;
  305.       histogram[r][g][b]++;
  306.     }
  307. }
  308.  
  309.  
  310.  
  311. /******************************/
  312. static CBOX *largest_box()
  313. /******************************/
  314. {
  315.   CBOX *tmp, *ptr;
  316.   int   size = -1;
  317.  
  318.   tmp = usedboxes;
  319.   ptr = NULL;
  320.  
  321.   while (tmp) {
  322.     if ( (tmp->rmax > tmp->rmin  ||
  323.       tmp->gmax > tmp->gmin  ||
  324.       tmp->bmax > tmp->bmin)  &&  tmp->total > size ) {
  325.       ptr = tmp;
  326.       size = tmp->total;
  327.     }
  328.     tmp = tmp->next;
  329.   }
  330.   return(ptr);
  331. }
  332.  
  333.  
  334.  
  335. /******************************/
  336. static void splitbox(ptr)
  337.      CBOX *ptr;
  338. /******************************/
  339. {
  340.   int   hist2[B_LEN], first, last, i, rdel, gdel, bdel;
  341.   CBOX *new;
  342.   int  *iptr, *histp, ir, ig, ib;
  343.   int  rmin,rmax,gmin,gmax,bmin,bmax;
  344.   enum {RED,GREEN,BLUE} which;
  345.  
  346.   /*
  347.    * see which axis is the largest, do a histogram along that
  348.    * axis.  Split at median point.  Contract both new boxes to
  349.    * fit points and return
  350.    */
  351.  
  352.   first = last = 0;   /* shut RT hcc compiler up */
  353.  
  354.   rmin = ptr->rmin;  rmax = ptr->rmax;
  355.   gmin = ptr->gmin;  gmax = ptr->gmax;
  356.   bmin = ptr->bmin;  bmax = ptr->bmax;
  357.  
  358.   rdel = rmax - rmin;
  359.   gdel = gmax - gmin;
  360.   bdel = bmax - bmin;
  361.  
  362.   if      (rdel>=gdel && rdel>=bdel) which = RED;
  363.   else if (gdel>=bdel)               which = GREEN;
  364.   else                               which = BLUE;
  365.  
  366.   /* get histogram along longest axis */
  367.   switch (which) {
  368.  
  369.   case RED:
  370.     histp = &hist2[rmin];
  371.     for (ir=rmin; ir<=rmax; ir++) {
  372.       *histp = 0;
  373.       for (ig=gmin; ig<=gmax; ig++) {
  374.     iptr = &histogram[ir][ig][bmin];
  375.     for (ib=bmin; ib<=bmax; ib++) {
  376.       *histp += *iptr;
  377.       ++iptr;
  378.     }
  379.       }
  380.       ++histp;
  381.     }
  382.     first = rmin;  last = rmax;
  383.     break;
  384.  
  385.   case GREEN:
  386.     histp = &hist2[gmin];
  387.     for (ig=gmin; ig<=gmax; ig++) {
  388.       *histp = 0;
  389.       for (ir=rmin; ir<=rmax; ir++) {
  390.     iptr = &histogram[ir][ig][bmin];
  391.     for (ib=bmin; ib<=bmax; ib++) {
  392.       *histp += *iptr;
  393.       ++iptr;
  394.     }
  395.       }
  396.       ++histp;
  397.     }
  398.     first = gmin;  last = gmax;
  399.     break;
  400.  
  401.   case BLUE:
  402.     histp = &hist2[bmin];
  403.     for (ib=bmin; ib<=bmax; ib++) {
  404.       *histp = 0;
  405.       for (ir=rmin; ir<=rmax; ir++) {
  406.     iptr = &histogram[ir][gmin][ib];
  407.     for (ig=gmin; ig<=gmax; ig++) {
  408.       *histp += *iptr;
  409.       iptr += B_LEN;
  410.     }
  411.       }
  412.       ++histp;
  413.     }
  414.     first = bmin;  last = bmax;
  415.     break;
  416.   }
  417.  
  418.  
  419.   /* find median point */
  420.   {
  421.     int sum, sum2;
  422.  
  423.     histp = &hist2[first];
  424.  
  425.     sum2 = ptr->total/2;
  426.     histp = &hist2[first];
  427.     sum = 0;
  428.         
  429.     for (i=first; i<=last && (sum += *histp++)<sum2; i++);
  430.     if (i==first) i++;
  431.   }
  432.  
  433.  
  434.   /* Create new box, re-allocate points */
  435.   
  436.   new = freeboxes;
  437.   freeboxes = new->next;
  438.   if (freeboxes) freeboxes->prev = NULL;
  439.  
  440.   if (usedboxes) usedboxes->prev = new;
  441.   new->next = usedboxes;
  442.   usedboxes = new;
  443.  
  444.   {
  445.     int sum1,sum2,j;
  446.     
  447.     histp = &hist2[first];
  448.     sum1 = 0;
  449.     for (j = first; j < i; ++j) sum1 += *histp++;
  450.     sum2 = 0;
  451.     for (j = i; j <= last; ++j) sum2 += *histp++;
  452.     new->total = sum1;
  453.     ptr->total = sum2;
  454.   }
  455.  
  456.  
  457.   new->rmin = rmin;  new->rmax = rmax;
  458.   new->gmin = gmin;  new->gmax = gmax;
  459.   new->bmin = bmin;  new->bmax = bmax;
  460.  
  461.   switch (which) {
  462.   case RED:    new->rmax = i-1;  ptr->rmin = i;  break;
  463.   case GREEN:  new->gmax = i-1;  ptr->gmin = i;  break;
  464.   case BLUE:   new->bmax = i-1;  ptr->bmin = i;  break;
  465.   }
  466.  
  467.   shrinkbox(new);
  468.   shrinkbox(ptr);
  469. }
  470.  
  471.  
  472. /****************************/
  473. static void shrinkbox(box)
  474.      CBOX *box;
  475. /****************************/
  476. {
  477.   int *histp,ir,ig,ib;
  478.   int  rmin,rmax,gmin,gmax,bmin,bmax;
  479.  
  480.   rmin = box->rmin;  rmax = box->rmax;
  481.   gmin = box->gmin;  gmax = box->gmax;
  482.   bmin = box->bmin;  bmax = box->bmax;
  483.  
  484.   if (rmax>rmin) {
  485.     for (ir=rmin; ir<=rmax; ir++)
  486.       for (ig=gmin; ig<=gmax; ig++) {
  487.     histp = &histogram[ir][ig][bmin];
  488.     for (ib=bmin; ib<=bmax; ib++)
  489.       if (*histp++ != 0) {
  490.         box->rmin = rmin = ir;
  491.         goto have_rmin;
  492.       }
  493.       }
  494.  
  495.   have_rmin:
  496.     if (rmax>rmin)
  497.       for (ir=rmax; ir>=rmin; --ir)
  498.     for (ig=gmin; ig<=gmax; ig++) {
  499.       histp = &histogram[ir][ig][bmin];
  500.       for (ib=bmin; ib<=bmax; ib++)
  501.         if (*histp++ != 0) {
  502.           box->rmax = rmax = ir;
  503.           goto have_rmax;
  504.         }
  505.     }
  506.   }
  507.  
  508.  
  509.  have_rmax:
  510.  
  511.   if (gmax>gmin) {
  512.     for (ig=gmin; ig<=gmax; ig++)
  513.       for (ir=rmin; ir<=rmax; ir++) {
  514.     histp = &histogram[ir][ig][bmin];
  515.     for (ib=bmin; ib<=bmax; ib++)
  516.       if (*histp++ != 0) {
  517.         box->gmin = gmin = ig;
  518.         goto have_gmin;
  519.       }
  520.       }
  521.   have_gmin:
  522.     if (gmax>gmin)
  523.       for (ig=gmax; ig>=gmin; --ig)
  524.     for (ir=rmin; ir<=rmax; ir++) {
  525.       histp = &histogram[ir][ig][bmin];
  526.       for (ib=bmin; ib<=bmax; ib++)
  527.         if (*histp++ != 0) {
  528.           box->gmax = gmax = ig;
  529.           goto have_gmax;
  530.         }
  531.     }
  532.   }
  533.  
  534.  
  535.  have_gmax:
  536.  
  537.   if (bmax>bmin) {
  538.     for (ib=bmin; ib<=bmax; ib++)
  539.       for (ir=rmin; ir<=rmax; ir++) {
  540.     histp = &histogram[ir][gmin][ib];
  541.     for (ig=gmin; ig<=gmax; ig++) {
  542.       if (*histp != 0) {
  543.         box->bmin = bmin = ib;
  544.         goto have_bmin;
  545.       }
  546.       histp += B_LEN;
  547.     }
  548.       }
  549.   have_bmin:
  550.     if (bmax>bmin)
  551.       for (ib=bmax; ib>=bmin; --ib)
  552.     for (ir=rmin; ir<=rmax; ir++) {
  553.       histp = &histogram[ir][gmin][ib];
  554.       for (ig=gmin; ig<=gmax; ig++) {
  555.         if (*histp != 0) {
  556.           bmax = ib;
  557.           goto have_bmax;
  558.         }
  559.         histp += B_LEN;
  560.       }
  561.     }
  562.   }
  563.  
  564.  have_bmax: return;
  565. }
  566.  
  567.  
  568.  
  569. /*******************************/
  570. static void assign_color(ptr,rp,gp,bp)
  571.      CBOX *ptr;
  572.      byte *rp,*gp,*bp;
  573. /*******************************/
  574. {
  575.   /* +1 ensures that color represents the middle of the box */
  576.   int r,g,b;
  577.  
  578.   r = ((ptr->rmin + ptr->rmax) << (COLOR_DEPTH - B_DEPTH)) / 2;
  579.   g = ((ptr->gmin + ptr->gmax) << (COLOR_DEPTH - B_DEPTH)) / 2;
  580.   b = ((ptr->bmin + ptr->bmax) << (COLOR_DEPTH - B_DEPTH)) / 2;
  581.  
  582.   *rp = (byte) r;  *gp = (byte) g;  *bp = (byte) b;
  583. }
  584.  
  585.  
  586.  
  587. /*******************************/
  588. static CCELL *create_colorcell(r1,g1,b1)
  589.      int r1,g1,b1;
  590. /*******************************/
  591. {
  592.   register int    i;
  593.   register CCELL *ptr;
  594.   register byte  *rp,*gp,*bp;
  595.   int             ir,ig,ib;
  596.   long            dist, mindist, tmp;
  597.  
  598.   ir = r1 >> (COLOR_DEPTH-C_DEPTH);
  599.   ig = g1 >> (COLOR_DEPTH-C_DEPTH);
  600.   ib = b1 >> (COLOR_DEPTH-C_DEPTH);
  601.  
  602.   r1 &= ~1 << (COLOR_DEPTH-C_DEPTH);
  603.   g1 &= ~1 << (COLOR_DEPTH-C_DEPTH);
  604.   b1 &= ~1 << (COLOR_DEPTH-C_DEPTH);
  605.  
  606.   ptr = (CCELL *) malloc(sizeof(CCELL));
  607.   *(ColorCells + ir*C_LEN*C_LEN + ig*C_LEN + ib) = ptr;
  608.   ptr->num_ents = 0;
  609.  
  610.   /* step 1: find all colors inside this cell, while we're at
  611.      it, find distance of centermost point to furthest
  612.      corner */
  613.  
  614.   mindist = 2000000000;
  615.  
  616.   rp=r;  gp=g;  bp=b;
  617.   for (i=0; i<num_colors; i++,rp++,gp++,bp++)
  618.     if( *rp>>(COLOR_DEPTH-C_DEPTH) == ir  &&
  619.         *gp>>(COLOR_DEPTH-C_DEPTH) == ig  &&
  620.         *bp>>(COLOR_DEPTH-C_DEPTH) == ib) {
  621.  
  622.       ptr->entries[ptr->num_ents][0] = i;
  623.       ptr->entries[ptr->num_ents][1] = 0;
  624.       ++ptr->num_ents;
  625.  
  626.       tmp = *rp - r1;
  627.       if (tmp < (MAX_COLOR/C_LEN/2)) tmp = MAX_COLOR/C_LEN-1 - tmp;
  628.       dist = (tmp*tmp * R2FACT);
  629.  
  630.       tmp = *gp - g1;
  631.       if (tmp < (MAX_COLOR/C_LEN/2)) tmp = MAX_COLOR/C_LEN-1 - tmp;
  632.       dist += (tmp*tmp * G2FACT);
  633.  
  634.       tmp = *bp - b1;
  635.       if (tmp < (MAX_COLOR/C_LEN/2)) tmp = MAX_COLOR/C_LEN-1 - tmp;
  636.       dist += (tmp*tmp * B2FACT);
  637.  
  638.       if (dist < mindist) mindist = dist;
  639.     }
  640.  
  641.  
  642.   /* step 3: find all points within that distance to box */
  643.  
  644.   rp=r;  gp=g;  bp=b;
  645.   for (i=0; i<num_colors; i++,rp++,gp++,bp++)
  646.     if (*rp >> (COLOR_DEPTH-C_DEPTH) != ir  ||
  647.     *gp >> (COLOR_DEPTH-C_DEPTH) != ig  ||
  648.     *bp >> (COLOR_DEPTH-C_DEPTH) != ib) {
  649.  
  650.       dist = 0;
  651.  
  652.       if ((tmp = r1 - *rp)>0 || (tmp = *rp - (r1 + MAX_COLOR/C_LEN-1)) > 0 )
  653.     dist += (tmp*tmp * R2FACT);
  654.  
  655.       if( (tmp = g1 - *gp)>0 || (tmp = *gp - (g1 + MAX_COLOR/C_LEN-1)) > 0 )
  656.     dist += (tmp*tmp * G2FACT);
  657.  
  658.       if( (tmp = b1 - *bp)>0 || (tmp = *bp - (b1 + MAX_COLOR/C_LEN-1)) > 0 )
  659.     dist += (tmp*tmp * B2FACT);
  660.  
  661.       if( dist < mindist ) {
  662.     ptr->entries[ptr->num_ents][0] = i;
  663.     ptr->entries[ptr->num_ents][1] = dist;
  664.     ++ptr->num_ents;
  665.       }
  666.     }
  667.  
  668.  
  669.   /* sort color cells by distance, use cheap exchange sort */
  670.   {
  671.     int n, next_n;
  672.  
  673.     n = ptr->num_ents - 1;
  674.     while (n>0) {
  675.       next_n = 0;
  676.       for (i=0; i<n; ++i) {
  677.     if (ptr->entries[i][1] > ptr->entries[i+1][1]) {
  678.       tmp = ptr->entries[i][0];
  679.       ptr->entries[i][0] = ptr->entries[i+1][0];
  680.       ptr->entries[i+1][0] = tmp;
  681.       tmp = ptr->entries[i][1];
  682.       ptr->entries[i][1] = ptr->entries[i+1][1];
  683.       ptr->entries[i+1][1] = tmp;
  684.       next_n = i;
  685.     }
  686.       }
  687.       n = next_n;
  688.     }
  689.   }
  690.   return (ptr);
  691. }
  692.  
  693.  
  694.  
  695.  
  696. /***************************/
  697. static void map_colortable()
  698. /***************************/
  699. {
  700.   int    ir,ig,ib, *histp;
  701.   CCELL *cell;
  702.  
  703.   histp  = &histogram[0][0][0];
  704.   for (ir=0; ir<B_LEN; ir++)
  705.     for (ig=0; ig<B_LEN; ig++)
  706.       for (ib=0; ib<B_LEN; ib++) {
  707.     if (*histp==0) *histp = -1;
  708.     else {
  709.       int  i, j;
  710.       long dist, d2, tmp;
  711.       
  712.       cell = *(ColorCells +
  713.            ( ((ir>>(B_DEPTH-C_DEPTH)) << C_DEPTH*2)
  714.            + ((ig>>(B_DEPTH-C_DEPTH)) << C_DEPTH)
  715.            +  (ib>>(B_DEPTH-C_DEPTH)) ) );
  716.         
  717.       if (cell==NULL)
  718.         cell = create_colorcell(ir<<(COLOR_DEPTH-B_DEPTH),
  719.                     ig<<(COLOR_DEPTH-B_DEPTH),
  720.                     ib<<(COLOR_DEPTH-B_DEPTH));
  721.  
  722.       dist = 2000000000;
  723.       for (i=0; i<cell->num_ents && dist>cell->entries[i][1]; i++) {
  724.         j = cell->entries[i][0];
  725.         d2 = r[j] - (ir << (COLOR_DEPTH-B_DEPTH));
  726.         d2 = (d2 * d2 * R2FACT);
  727.         tmp = g[j] - (ig << (COLOR_DEPTH-B_DEPTH));
  728.         d2 += (tmp*tmp * G2FACT);
  729.         tmp = b[j] - (ib << (COLOR_DEPTH-B_DEPTH));
  730.         d2 += (tmp*tmp * B2FACT);
  731.         if( d2 < dist ) { dist = d2;  *histp = j; }
  732.       }
  733.     }
  734.     histp++;
  735.       }
  736. }
  737.  
  738.  
  739.  
  740. /*****************************/
  741. static int quant_fsdither()
  742. /*****************************/
  743. {
  744.   register int  *thisptr, *nextptr;
  745.   int           *thisline, *nextline, *tmpptr;
  746.   int            r1, g1, b1, r2, g2, b2, r3, g3, b3;
  747.   int            i, j, k, d, imax, jmax, oval;
  748.   byte          *inptr, *outptr;
  749.   int            lastline, lastpixel;
  750.  
  751.   imax = HIGH - 1;
  752.   jmax = WIDE - 1;
  753.   
  754.   thisline = (int *) malloc(WIDE * 3 * sizeof(int));
  755.   nextline = (int *) malloc(WIDE * 3 * sizeof(int));
  756.  
  757.   if (thisline == NULL || nextline == NULL) {
  758.     fprintf(stderr,"%s: unable to allocate stuff for the 'dither' routine\n",
  759.         cmd);
  760.     return 1;
  761.   }
  762.  
  763.  
  764.   inptr  = (byte *) pic24;
  765.   outptr = (byte *) pic;
  766.  
  767.   /* get first line of picture */
  768.   for (j=WIDE * 3, tmpptr=nextline; j; j--)
  769.     *tmpptr++ = (int) *inptr++;
  770.  
  771.   for (i=0; i<HIGH; i++) {
  772.     /* swap thisline and nextline */
  773.     tmpptr = thisline;  thisline = nextline;  nextline = tmpptr;
  774.     lastline = (i==imax);
  775.  
  776.     /* read in next line */
  777.     if (!lastline)
  778.       for (j=WIDE * 3, tmpptr=nextline; j; j--)
  779.     *tmpptr++ = (int) *inptr++;
  780.  
  781.     /* dither this line and put it into the output picture */
  782.     thisptr = thisline;  nextptr = nextline;
  783.  
  784.     for (j=0; j<WIDE; j++) {
  785.       lastpixel = (j==jmax);
  786.  
  787.       r2 = *thisptr++;  g2 = *thisptr++;  b2 = *thisptr++;
  788.  
  789.       RANGE(r2, 0, MAX_COLOR-1);
  790.       RANGE(g2, 0, MAX_COLOR-1);
  791.       RANGE(b2, 0, MAX_COLOR-1);
  792.  
  793.       r1 = r2;  g1 = g2;  b1 = b2;
  794.  
  795.       r2 >>= (COLOR_DEPTH-B_DEPTH);
  796.       g2 >>= (COLOR_DEPTH-B_DEPTH);
  797.       b2 >>= (COLOR_DEPTH-B_DEPTH);
  798.  
  799.       if ( (oval=histogram[r2][g2][b2]) == -1) {
  800.     int ci, cj;
  801.     long dist, d2, tmp;
  802.     CCELL *cell;
  803.  
  804.     cell = *( ColorCells + 
  805.         ( ((r2>>(B_DEPTH-C_DEPTH)) << C_DEPTH*2)
  806.             + ((g2>>(B_DEPTH-C_DEPTH)) << C_DEPTH )
  807.         +  (b2>>(B_DEPTH-C_DEPTH)) ) );
  808.           
  809.     if (cell==NULL) cell = create_colorcell(r1,g1,b1);
  810.  
  811.     dist = 2000000000;
  812.     for (ci=0; ci<cell->num_ents && dist>cell->entries[ci][1]; ci++) {
  813.       cj = cell->entries[ci][0];
  814.       d2 = (r[cj] >> (COLOR_DEPTH-B_DEPTH)) - r2;
  815.       d2 = (d2*d2 * R2FACT);
  816.       tmp = (g[cj] >> (COLOR_DEPTH-B_DEPTH)) - g2;
  817.       d2 += (tmp*tmp * G2FACT);
  818.       tmp = (b[cj] >> (COLOR_DEPTH-B_DEPTH)) - b2;
  819.       d2 += (tmp*tmp * B2FACT);
  820.       if (d2<dist) { dist = d2;  oval = cj; }
  821.     }
  822.     histogram[r2][g2][b2] = oval;
  823.       }
  824.  
  825.       *outptr++ = oval;
  826.  
  827.       r1 -= r[oval];  g1 -= g[oval];  b1 -= b[oval];
  828.     
  829.       /* don't use tables, because r1,g1,b1 could go negative */
  830.       if (!lastpixel) {
  831.     thisptr[0] += (r1<0) ? (r1*7-8)/16 : (r1*7+8)/16;
  832.     thisptr[1] += (g1<0) ? (g1*7-8)/16 : (g1*7+8)/16;
  833.     thisptr[2] += (b1<0) ? (b1*7-8)/16 : (b1*7+8)/16;
  834.       }
  835.     
  836.       if (!lastline) {
  837.     if (j) {
  838.       nextptr[-3] += (r1<0) ? (r1*3-8)/16 : (r1*3+8)/16;
  839.       nextptr[-2] += (g1<0) ? (g1*3-8)/16 : (g1*3+8)/16;
  840.       nextptr[-1] += (b1<0) ? (b1*3-8)/16 : (b1*3+8)/16;
  841.     }
  842.  
  843.     nextptr[0] += (r1<0) ? (r1*5-8)/16 : (r1*5+8)/16;
  844.     nextptr[1] += (g1<0) ? (g1*5-8)/16 : (g1*5+8)/16;
  845.     nextptr[2] += (b1<0) ? (b1*5-8)/16 : (b1*5+8)/16;
  846.     
  847.     if (!lastpixel) {
  848.       nextptr[3] += (r1<0) ? (r1-8)/16 : (r1+8)/16;
  849.       nextptr[4] += (g1<0) ? (g1-8)/16 : (g1+8)/16;
  850.       nextptr[5] += (b1<0) ? (b1-8)/16 : (b1+8)/16;
  851.     }
  852.     nextptr += 3;
  853.       }
  854.     }
  855.   }
  856.  
  857.   free(thisline);  free(nextline);
  858.   return 0;
  859. }
  860.  
  861.  
  862.  
  863.  
  864.  
  865.  
  866.  
  867.  
  868.  
  869.  
  870. /************************************/
  871. static int Quick24to8(p24,w,h)
  872. byte *p24;
  873. int   w,h;
  874. {
  875.  
  876.   /* floyd-steinberg dithering.
  877.    *
  878.    * ----   x    7/16
  879.    * 3/16  5/16  1/16
  880.    *
  881.    */
  882.  
  883.   /* called after 'pic' has been alloced, pWIDE,pHIGH set up, mono/1-bit
  884.      checked already */
  885.  
  886.   byte *pp;
  887.   int  r1, g1, b1;
  888.   int  *thisline, *nextline, *thisptr, *nextptr, *tmpptr;
  889.   int  i, j, val, pwide3;
  890.   int  imax, jmax;
  891.  
  892.   pp = pic;  pwide3 = w * 3;  imax = h-1;  jmax = w-1;
  893.  
  894.   /* up to 256 colors:     3 bits R, 3 bits G, 2 bits B  (RRRGGGBB) */
  895. #define RMASK 0xe0
  896. #define R_SHIFT        0
  897. #define GMASK 0xe0
  898. #define G_SHIFT        3
  899. #define BMASK 0xc0
  900. #define B_SHIFT        6
  901.  
  902.   /* load up colormap */
  903.   /* note that 0 and 255 of each color are always in the map; */
  904.   /* intermediate values are evenly spaced. */
  905.   for (i=0; i<256; i++) {
  906.     r[i] = (((i<<R_SHIFT) & RMASK) * 255 + RMASK/2) / RMASK;
  907.     g[i] = (((i<<G_SHIFT) & GMASK) * 255 + GMASK/2) / GMASK;
  908.     b[i] = (((i<<B_SHIFT) & BMASK) * 255 + BMASK/2) / BMASK;
  909.   }
  910.  
  911.  
  912.   thisline = (int *) malloc(pwide3 * sizeof(int));
  913.   nextline = (int *) malloc(pwide3 * sizeof(int));
  914.   if (!thisline || !nextline) {
  915.     fprintf(stderr,"%s: unable to allocate memory in Quick24to8()\n", cmd);
  916.     return(1);
  917.     }
  918.  
  919.   /* get first line of picture */
  920.   for (j=pwide3, tmpptr=nextline; j; j--) *tmpptr++ = (int) *p24++;
  921.  
  922.   for (i=0; i<h; i++) {
  923.     tmpptr = thisline;  thisline = nextline;  nextline = tmpptr;   /* swap */
  924.  
  925.     if (i!=imax)   /* get next line */
  926.       for (j=pwide3, tmpptr=nextline; j; j--)
  927.     *tmpptr++ = (int) *p24++;
  928.  
  929.     for (j=0, thisptr=thisline, nextptr=nextline; j<w; j++,pp++) {
  930.       r1 = *thisptr++;  g1 = *thisptr++;  b1 = *thisptr++;
  931.       RANGE(r1,0,255);  RANGE(g1,0,255);  RANGE(b1,0,255);  
  932.  
  933.       /* choose actual pixel value */
  934.       val = (((r1&RMASK)>>R_SHIFT) | ((g1&GMASK)>>G_SHIFT) | 
  935.          ((b1&BMASK)>>B_SHIFT));
  936.       *pp = val;
  937.  
  938.       /* compute color errors */
  939.       r1 -= r[val];
  940.       g1 -= g[val];
  941.       b1 -= b[val];
  942.  
  943.       /* Add fractions of errors to adjacent pixels */
  944.       r1 += 256;              /* make positive for table indexing */
  945.       g1 += 256;
  946.       b1 += 256;
  947.  
  948.       if (j!=jmax) {  /* adjust RIGHT pixel */
  949.     thisptr[0] += tbl7[r1];
  950.     thisptr[1] += tbl7[g1];
  951.     thisptr[2] += tbl7[b1];
  952.       }
  953.       
  954.       if (i!=imax) {    /* do BOTTOM pixel */
  955.     nextptr[0] += tbl5[r1];
  956.     nextptr[1] += tbl5[g1];
  957.     nextptr[2] += tbl5[b1];
  958.  
  959.     if (j>0) {  /* do BOTTOM LEFT pixel */
  960.       nextptr[-3] += tbl3[r1];
  961.       nextptr[-2] += tbl3[g1];
  962.       nextptr[-1] += tbl3[b1];
  963.     }
  964.  
  965.     if (j!=jmax) {  /* do BOTTOM RIGHT pixel */
  966.       nextptr[3] += tbl1[r1];
  967.       nextptr[4] += tbl1[g1];
  968.       nextptr[5] += tbl1[b1];
  969.     }
  970.     nextptr += 3;
  971.       }
  972.     }
  973.   }
  974.  
  975.   return 0;
  976. }
  977.       
  978.  
  979.  
  980. /****************************/
  981. static int QuickCheck(pic24,w,h,maxcol)
  982. byte *pic24;
  983. int   w,h,maxcol;
  984. {
  985.   /* scans picture until it finds more than 'maxcol' different colors.  If it
  986.      finds more than 'maxcol' colors, it returns '0'.  If it DOESN'T, it does
  987.      the 24-to-8 conversion by simply sticking the colors it found into
  988.      a colormap, and changing instances of a color in pic24 into colormap
  989.      indicies (in pic) */
  990.  
  991.   unsigned long colors[256],col;
  992.   int           i, nc, low, high, mid;
  993.   byte         *p, *pix;
  994.  
  995.   if (maxcol>256) maxcol = 256;
  996.  
  997.   /* put the first color in the table by hand */
  998.   nc = 0;  mid = 0;  
  999.  
  1000.   for (i=w*h,p=pic24; i; i--) {
  1001.     col  = (((u_long) *p++) << 16);  
  1002.     col += (((u_long) *p++) << 8);
  1003.     col +=  *p++;
  1004.  
  1005.     /* binary search the 'colors' array to see if it's in there */
  1006.     low = 0;  high = nc-1;
  1007.     while (low <= high) {
  1008.       mid = (low+high)/2;
  1009.       if      (col < colors[mid]) high = mid - 1;
  1010.       else if (col > colors[mid]) low  = mid + 1;
  1011.       else break;
  1012.     }
  1013.  
  1014.     if (high < low) { /* didn't find color in list, add it. */
  1015.       /* WARNING: this is an overlapped memory copy.  memcpy doesn't do
  1016.      it correctly, hence 'bcopy', which claims to */
  1017.       if (nc>=maxcol) return 0;
  1018.       bcopy(&colors[low], &colors[low+1], (nc - low) * sizeof(unsigned long));
  1019.       colors[low] = col;
  1020.       nc++;
  1021.     }
  1022.   }
  1023.  
  1024.  
  1025.   /* run through the data a second time, this time mapping pixel values in
  1026.      pic24 into colormap offsets into 'colors' */
  1027.  
  1028.   for (i=w*h,p=pic24, pix=pic; i; i--,pix++) {
  1029.     col  = (((u_long) *p++) << 16);  
  1030.     col += (((u_long) *p++) << 8);
  1031.     col +=  *p++;
  1032.  
  1033.     /* binary search the 'colors' array.  It *IS* in there */
  1034.     low = 0;  high = nc-1;
  1035.     while (low <= high) {
  1036.       mid = (low+high)/2;
  1037.       if      (col < colors[mid]) high = mid - 1;
  1038.       else if (col > colors[mid]) low  = mid + 1;
  1039.       else break;
  1040.     }
  1041.  
  1042.     if (high < low) {
  1043.       fprintf(stderr,"QuickCheck:  impossible situation!\n");
  1044.       exit(1);
  1045.     }
  1046.     *pix = mid;
  1047.   }
  1048.  
  1049.   /* and load up the 'desired colormap' */
  1050.   for (i=0; i<nc; i++) {
  1051.     r[i] =  colors[i]>>16;  
  1052.     g[i] = (colors[i]>>8) & 0xff;
  1053.     b[i] =  colors[i]     & 0xff;
  1054.   }
  1055.  
  1056.   return 1;
  1057. }
  1058.  
  1059.  
  1060.  
  1061.  
  1062.  
  1063.  
  1064.  
  1065.  
  1066.  
  1067.  
  1068.  
  1069.  
  1070. /***************************************************************/
  1071. /* The following code based on code from the 'pbmplus' package */
  1072. /* written by Jef Poskanzer                                    */
  1073. /***************************************************************/
  1074.  
  1075.  
  1076. /* ppmquant.c - quantize the colors in a pixmap down to a specified number
  1077. **
  1078. ** Copyright (C) 1989, 1991 by Jef Poskanzer.
  1079. **
  1080. ** Permission to use, copy, modify, and distribute this software and its
  1081. ** documentation for any purpose and without fee is hereby granted, provided
  1082. ** that the above copyright notice appear in all copies and that both that
  1083. ** copyright notice and this permission notice appear in supporting
  1084. ** documentation.  This software is provided "as is" without express or
  1085. ** implied warranty.
  1086. */
  1087.  
  1088.  
  1089. #undef max
  1090. #define max(a,b) ((a) > (b) ? (a) : (b))
  1091. #undef min
  1092. #define min(a,b) ((a) < (b) ? (a) : (b))
  1093. #undef abs
  1094. #define abs(a) ((a) >= 0 ? (a) : -(a))
  1095. #undef odd
  1096. #define odd(n) ((n) & 1)
  1097.  
  1098. typedef unsigned char pixval;
  1099.  
  1100. #define PPM_MAXMAXVAL 255
  1101. typedef struct { pixval r, g, b; } pixel;
  1102.  
  1103. #define PPM_GETR(p) ((p).r)
  1104. #define PPM_GETG(p) ((p).g)
  1105. #define PPM_GETB(p) ((p).b)
  1106.  
  1107. #define PPM_ASSIGN(p,red,grn,blu) \
  1108.   { (p).r = (red); (p).g = (grn); (p).b = (blu); }
  1109.  
  1110. #define PPM_EQUAL(p,q) ( (p).r == (q).r && (p).g == (q).g && (p).b == (q).b )
  1111.  
  1112.  
  1113. /* Color scaling macro -- to make writing ppmtowhatever easier. */
  1114.  
  1115. #define PPM_DEPTH(newp,p,oldmaxval,newmaxval) \
  1116.     PPM_ASSIGN( (newp), \
  1117.             (int) PPM_GETR(p) * (newmaxval) / (oldmaxval), \
  1118.             (int) PPM_GETG(p) * (newmaxval) / (oldmaxval), \
  1119.             (int) PPM_GETB(p) * (newmaxval) / (oldmaxval) )
  1120.  
  1121.  
  1122. /* Luminance macro. */
  1123.  
  1124. /* 
  1125.  * #define PPM_LUMIN(p) \
  1126.  *   ( 0.299 * PPM_GETR(p) + 0.587 * PPM_GETG(p) + 0.114 * PPM_GETB(p) )
  1127.  */
  1128.  
  1129. /* Luminance macro, using only integer ops.  Returns an int (*256)  JHB */
  1130. #define PPM_LUMIN(p) \
  1131.   ( 77 * PPM_GETR(p) + 150 * PPM_GETG(p) + 29 * PPM_GETB(p) )
  1132.  
  1133. /* Color histogram stuff. */
  1134.  
  1135. typedef struct colorhist_item* colorhist_vector;
  1136. struct colorhist_item { pixel color;
  1137.             int value;
  1138.               };
  1139.  
  1140. typedef struct colorhist_list_item* colorhist_list;
  1141. struct colorhist_list_item { struct colorhist_item ch;
  1142.                  colorhist_list next;
  1143.                };
  1144.  
  1145. typedef colorhist_list* colorhash_table;
  1146.  
  1147.  
  1148. #define MAXCOLORS 32767
  1149. #define CLUSTER_MAXVAL 63
  1150.  
  1151. #define LARGE_LUM
  1152. #define REP_AVERAGE_PIXELS
  1153.  
  1154. #define FS_SCALE 1024
  1155.  
  1156. #define HASH_SIZE 6553
  1157.  
  1158. #define ppm_hashpixel(p) ((((int) PPM_GETR(p) * 33023 +    \
  1159.                 (int) PPM_GETG(p) * 30013 +    \
  1160.                 (int) PPM_GETB(p) * 27011) & 0x7fffffff)   \
  1161.               % HASH_SIZE)
  1162.  
  1163.  
  1164.  
  1165. /*** function defs ***/
  1166. static colorhist_vector mediancut();
  1167. static int              redcompare(), greencompare(), bluecompare();
  1168. static int              sumcompare();
  1169. static colorhist_vector ppm_computecolorhist(), ppm_colorhashtocolorhist();
  1170. static colorhash_table  ppm_computecolorhash(), ppm_alloccolorhash();
  1171. static void             ppm_addtocolorhash();
  1172. static void             ppm_freecolorhist(), ppm_freecolorhash();
  1173. static int              ppm_lookupcolor();
  1174.  
  1175.  
  1176. /****************************************************************************/
  1177. static int ppmquant(pic24, cols, rows, newcolors)
  1178. byte *pic24;
  1179. int  cols, rows, newcolors;
  1180. {
  1181.   pixel**          pixels;
  1182.   register pixel*  pP;
  1183.   int              argn, row;
  1184.   register int     col, limitcol;
  1185.   pixval           maxval, newmaxval;
  1186.   int              colors;
  1187.   register int     index;
  1188.   colorhist_vector chv, colormap;
  1189.   colorhash_table  cht;
  1190.   int              verbose, floyd;
  1191.   long             *thisrerr, *nextrerr, *thisgerr, *nextgerr;
  1192.   long             *thisberr, *nextberr, *temperr;
  1193.   register long    sr, sg, sb, err;
  1194.   int              fs_direction, i;
  1195.   unsigned char    *picptr;
  1196.   static char      *fn = "ppmquant()";
  1197.  
  1198.   maxval = 255;
  1199.  
  1200.   /*
  1201.    *  reformat 24-bit pic24 image (3 bytes per pixel) into 2-dimensional
  1202.    *  array of pixel structures
  1203.    */
  1204.  
  1205.   if (DEBUG) fprintf(stderr,"%s: remapping to ppm-style internal fmt\n", fn);
  1206.  
  1207.   pixels = (pixel **) malloc(rows * sizeof(pixel *));
  1208.   if (!pixels) FatalError("couldn't allocate 'pixels' array");
  1209.   for (row=0; row<rows; row++) {
  1210.     pixels[row] = (pixel *) malloc(cols * sizeof(pixel));
  1211.     if (!pixels[row]) FatalError("couldn't allocate a row of pixels array");
  1212.  
  1213.     for (col=0, pP=pixels[row]; col<cols; col++, pP++) {
  1214.       pP->r = *pic24++;
  1215.       pP->g = *pic24++;
  1216.       pP->b = *pic24++;
  1217.     }
  1218.   }
  1219.   if (DEBUG) fprintf(stderr,"%s: done format remapping\n", fn);
  1220.  
  1221.  
  1222.     
  1223.  
  1224.   /*
  1225.    *  attempt to make a histogram of the colors, unclustered.
  1226.    *  If at first we don't succeed, lower maxval to increase color
  1227.    *  coherence and try again.  This will eventually terminate, with
  1228.    *  maxval at worst 15, since 32^3 is approximately MAXCOLORS.
  1229.    */
  1230.  
  1231.   for ( ; ; ) {
  1232.     if (DEBUG) fprintf(stderr, "%s:  making histogram\n", fn);
  1233.  
  1234.     chv = ppm_computecolorhist(pixels, cols, rows, MAXCOLORS, &colors);
  1235.     if (chv != (colorhist_vector) 0) break;
  1236.     
  1237.     if (DEBUG) fprintf(stderr, "%s: too many colors!\n", fn);
  1238.     newmaxval = maxval / 2;
  1239.     if (DEBUG) fprintf(stderr, "%s: rescaling colors (maxval=%d) %s\n",
  1240.                fn, newmaxval, "to improve clustering");
  1241.  
  1242.     for (row=0; row<rows; ++row)
  1243.       for (col=0, pP=pixels[row]; col<cols; ++col, ++pP)
  1244.     PPM_DEPTH( *pP, *pP, maxval, newmaxval );
  1245.     maxval = newmaxval;
  1246.   }
  1247.  
  1248.   if (DEBUG) fprintf(stderr,"%s: %d colors found\n", fn, colors);
  1249.  
  1250.  
  1251.  
  1252.   /*
  1253.    * Step 3: apply median-cut to histogram, making the new colormap.
  1254.    */
  1255.  
  1256.   if (DEBUG) fprintf(stderr, "%s: choosing %d colors\n", fn, newcolors);
  1257.   colormap = mediancut(chv, colors, rows * cols, maxval, newcolors);
  1258.   ppm_freecolorhist(chv);
  1259.  
  1260.  
  1261.  
  1262.   /*
  1263.    *  Step 4: map the colors in the image to their closest match in the
  1264.    *  new colormap, and write 'em out.
  1265.    */
  1266.  
  1267.   if (DEBUG) fprintf(stderr,"%s: mapping image to new colors\n", fn);
  1268.   cht = ppm_alloccolorhash();
  1269.  
  1270.   picptr = pic;
  1271.   for (row = 0;  row < rows;  ++row) {
  1272.     col = 0;  limitcol = cols;  pP = pixels[row];
  1273.  
  1274.     do {
  1275.       int hash;
  1276.       colorhist_list chl;
  1277.  
  1278.       /* Check hash table to see if we have already matched this color. */
  1279.       /* index = ppm_lookupcolor( cht, pP ); */
  1280.  
  1281.       hash = ppm_hashpixel(*pP);
  1282.       for (chl = cht[hash];  chl;  chl = chl->next)
  1283.     if (PPM_EQUAL(chl->ch.color, *pP)) {index = chl->ch.value; break;}
  1284.  
  1285.       if (!chl /*index = -1*/) {/* No; search colormap for closest match. */
  1286.     register int i, r1, g1, b1, r2, g2, b2;
  1287.     register long dist, newdist;
  1288.  
  1289.     r1 = PPM_GETR( *pP );
  1290.     g1 = PPM_GETG( *pP );
  1291.     b1 = PPM_GETB( *pP );
  1292.     dist = 2000000000;
  1293.  
  1294.     for (i=0; i<newcolors; i++) {
  1295.       r2 = PPM_GETR( colormap[i].color );
  1296.       g2 = PPM_GETG( colormap[i].color );
  1297.       b2 = PPM_GETB( colormap[i].color );
  1298.  
  1299.       newdist = ( r1 - r2 ) * ( r1 - r2 ) +
  1300.                 ( g1 - g2 ) * ( g1 - g2 ) +
  1301.                 ( b1 - b2 ) * ( b1 - b2 );
  1302.  
  1303.       if (newdist<dist) { index = i;  dist = newdist; }
  1304.     }
  1305.  
  1306.     /* ppm_addtocolorhash(cht, pP, index); */
  1307.     hash = ppm_hashpixel(*pP);
  1308.     chl = (colorhist_list) malloc(sizeof(struct colorhist_list_item));
  1309.     if (!chl) FatalError("ran out of memory adding to hash table");
  1310.  
  1311.     chl->ch.color = *pP;
  1312.     chl->ch.value = index;
  1313.     chl->next = cht[hash];
  1314.     cht[hash] = chl;
  1315.       }
  1316.  
  1317.       *picptr++ = index;
  1318.  
  1319.       ++col;
  1320.       ++pP;
  1321.     }
  1322.     while (col != limitcol);
  1323.   }
  1324.  
  1325.   /* rescale the colormap and load the XV colormap */
  1326.   for (i=0; i<newcolors; i++) {
  1327.     PPM_DEPTH(colormap[i].color, colormap[i].color, maxval, 255);
  1328.     r[i] = PPM_GETR( colormap[i].color );
  1329.     g[i] = PPM_GETG( colormap[i].color );
  1330.     b[i] = PPM_GETB( colormap[i].color );
  1331.   }
  1332.  
  1333.   /* free the pixels array */
  1334.   for (i=0; i<rows; i++) free(pixels[i]);
  1335.   free(pixels);
  1336.  
  1337.   /* free cht and colormap */
  1338.   ppm_freecolorhist(colormap);
  1339.   ppm_freecolorhash(cht);
  1340.  
  1341.   return 0;
  1342. }
  1343.  
  1344.  
  1345.  
  1346. /*
  1347. ** Here is the fun part, the median-cut colormap generator.  This is based
  1348. ** on Paul Heckbert's paper "Color Image Quantization for Frame Buffer
  1349. ** Display", SIGGRAPH '82 Proceedings, page 297.
  1350. */
  1351.  
  1352. typedef struct box* box_vector;
  1353. struct box {
  1354.   int index;
  1355.   int colors;
  1356.   int sum;
  1357. };
  1358.  
  1359.  
  1360. /****************************************************************************/
  1361. static colorhist_vector mediancut( chv, colors, sum, maxval, newcolors )
  1362.      colorhist_vector chv;
  1363.      int colors, sum, newcolors;
  1364.      pixval maxval;
  1365. {
  1366.   colorhist_vector colormap;
  1367.   box_vector bv;
  1368.   register int bi, i;
  1369.   int boxes;
  1370.  
  1371.   bv = (box_vector) malloc(sizeof(struct box) * newcolors);
  1372.   colormap = (colorhist_vector) 
  1373.              malloc(sizeof(struct colorhist_item) * newcolors );
  1374.  
  1375.   if (!bv || !colormap) FatalError("unable to malloc in mediancut()");
  1376.  
  1377.   for (i=0; i<newcolors; i++)
  1378.     PPM_ASSIGN(colormap[i].color, 0, 0, 0);
  1379.  
  1380.   /*
  1381.    *  Set up the initial box.
  1382.    */
  1383.   bv[0].index = 0;
  1384.   bv[0].colors = colors;
  1385.   bv[0].sum = sum;
  1386.   boxes = 1;
  1387.  
  1388.  
  1389.   /*
  1390.    ** Main loop: split boxes until we have enough.
  1391.    */
  1392.  
  1393.   while ( boxes < newcolors ) {
  1394.     register int indx, clrs;
  1395.     int sm;
  1396.     register int minr, maxr, ming, maxg, minb, maxb, v;
  1397.     int halfsum, lowersum;
  1398.  
  1399.     /*
  1400.      ** Find the first splittable box.
  1401.      */
  1402.     for (bi=0; bv[bi].colors<2 && bi<boxes; bi++) ;
  1403.     if (bi == boxes) break;    /* ran out of colors! */
  1404.  
  1405.     indx = bv[bi].index;
  1406.     clrs = bv[bi].colors;
  1407.     sm = bv[bi].sum;
  1408.  
  1409.     /*
  1410.      ** Go through the box finding the minimum and maximum of each
  1411.      ** component - the boundaries of the box.
  1412.      */
  1413.     minr = maxr = PPM_GETR( chv[indx].color );
  1414.     ming = maxg = PPM_GETG( chv[indx].color );
  1415.     minb = maxb = PPM_GETB( chv[indx].color );
  1416.  
  1417.     for (i=1; i<clrs; i++) {
  1418.       v = PPM_GETR( chv[indx + i].color );
  1419.       if (v < minr) minr = v;
  1420.       if (v > maxr) maxr = v;
  1421.  
  1422.       v = PPM_GETG( chv[indx + i].color );
  1423.       if (v < ming) ming = v;
  1424.       if (v > maxg) maxg = v;
  1425.  
  1426.       v = PPM_GETB( chv[indx + i].color );
  1427.       if (v < minb) minb = v;
  1428.       if (v > maxb) maxb = v;
  1429.     }
  1430.  
  1431.     /*
  1432.      ** Find the largest dimension, and sort by that component.  I have
  1433.      ** included two methods for determining the "largest" dimension;
  1434.      ** first by simply comparing the range in RGB space, and second
  1435.      ** by transforming into luminosities before the comparison.  You
  1436.      ** can switch which method is used by switching the commenting on
  1437.      ** the LARGE_ defines at the beginning of this source file.
  1438.      */
  1439.     {
  1440.       /* LARGE_LUM version */
  1441.  
  1442.       pixel p;
  1443.       int rl, gl, bl;
  1444.  
  1445.       PPM_ASSIGN(p, maxr - minr, 0, 0);
  1446.       rl = PPM_LUMIN(p);
  1447.  
  1448.       PPM_ASSIGN(p, 0, maxg - ming, 0);
  1449.       gl = PPM_LUMIN(p);
  1450.  
  1451.       PPM_ASSIGN(p, 0, 0, maxb - minb);
  1452.       bl = PPM_LUMIN(p);
  1453.  
  1454.       if (rl >= gl && rl >= bl)
  1455.     qsort( (char*) &(chv[indx]), clrs, sizeof(struct colorhist_item),
  1456.           redcompare );
  1457.       else if (gl >= bl)
  1458.     qsort( (char*) &(chv[indx]), clrs, sizeof(struct colorhist_item),
  1459.           greencompare );
  1460.       else qsort( (char*) &(chv[indx]), clrs, sizeof(struct colorhist_item),
  1461.          bluecompare );
  1462.     }
  1463.  
  1464.     /*
  1465.      ** Now find the median based on the counts, so that about half the
  1466.      ** pixels (not colors, pixels) are in each subdivision.
  1467.      */
  1468.     lowersum = chv[indx].value;
  1469.     halfsum = sm / 2;
  1470.     for (i=1; i<clrs-1; i++) {
  1471.       if (lowersum >= halfsum) break;
  1472.       lowersum += chv[indx + i].value;
  1473.     }
  1474.  
  1475.     /*
  1476.      ** Split the box, and sort to bring the biggest boxes to the top.
  1477.      */
  1478.     bv[bi].colors = i;
  1479.     bv[bi].sum = lowersum;
  1480.     bv[boxes].index = indx + i;
  1481.     bv[boxes].colors = clrs - i;
  1482.     bv[boxes].sum = sm - lowersum;
  1483.     ++boxes;
  1484.     qsort( (char*) bv, boxes, sizeof(struct box), sumcompare );
  1485.   }  /* while (boxes ... */
  1486.  
  1487.   /*
  1488.    ** Ok, we've got enough boxes.  Now choose a representative color for
  1489.    ** each box.  There are a number of possible ways to make this choice.
  1490.    ** One would be to choose the center of the box; this ignores any structure
  1491.    ** within the boxes.  Another method would be to average all the colors in
  1492.    ** the box - this is the method specified in Heckbert's paper.  A third
  1493.    ** method is to average all the pixels in the box.  You can switch which
  1494.    ** method is used by switching the commenting on the REP_ defines at
  1495.    ** the beginning of this source file.
  1496.    */
  1497.   
  1498.   for (bi=0; bi<boxes; bi++) {
  1499.     /* REP_AVERAGE_PIXELS version */
  1500.     register int indx = bv[bi].index;
  1501.     register int clrs = bv[bi].colors;
  1502.     register long r = 0, g = 0, b = 0, sum = 0;
  1503.  
  1504.     for (i=0; i<clrs; i++) {
  1505.       r += PPM_GETR( chv[indx + i].color ) * chv[indx + i].value;
  1506.       g += PPM_GETG( chv[indx + i].color ) * chv[indx + i].value;
  1507.       b += PPM_GETB( chv[indx + i].color ) * chv[indx + i].value;
  1508.       sum += chv[indx + i].value;
  1509.     }
  1510.  
  1511.     r = r / sum;  if (r>maxval) r = maxval;    /* avoid math errors */
  1512.     g = g / sum;  if (g>maxval) g = maxval;
  1513.     b = b / sum;  if (b>maxval) b = maxval;
  1514.  
  1515.     PPM_ASSIGN( colormap[bi].color, r, g, b );
  1516.   }
  1517.  
  1518.   free(bv);
  1519.   return colormap;
  1520. }
  1521.  
  1522.  
  1523. /**********************************/
  1524. static int redcompare(ch1, ch2)
  1525. colorhist_vector ch1, ch2;
  1526. {
  1527.   return (int) PPM_GETR( ch1->color ) - (int) PPM_GETR( ch2->color );
  1528. }
  1529.  
  1530. /**********************************/
  1531. static int greencompare(ch1, ch2)
  1532. colorhist_vector ch1, ch2;
  1533. {
  1534.   return (int) PPM_GETG( ch1->color ) - (int) PPM_GETG( ch2->color );
  1535. }
  1536.  
  1537. /**********************************/
  1538. static int bluecompare( ch1, ch2 )
  1539. colorhist_vector ch1, ch2;
  1540. {
  1541.   return (int) PPM_GETB( ch1->color ) - (int) PPM_GETB( ch2->color );
  1542. }
  1543.  
  1544. /**********************************/
  1545. static int sumcompare( b1, b2 )
  1546. box_vector b1, b2;
  1547. {
  1548.   return b2->sum - b1->sum;
  1549. }
  1550.  
  1551.  
  1552.  
  1553. /****************************************************************************/
  1554. static colorhist_vector  ppm_computecolorhist(pixels, cols, rows, 
  1555.                           maxcolors, colorsP)
  1556.      pixel** pixels;
  1557.      int cols, rows;
  1558.      int* colorsP;
  1559. {
  1560.   colorhash_table cht;
  1561.   colorhist_vector chv;
  1562.  
  1563.   cht = ppm_computecolorhash(pixels, cols, rows, maxcolors, colorsP);
  1564.   if (!cht) return (colorhist_vector) 0;
  1565.  
  1566.   chv = ppm_colorhashtocolorhist(cht, maxcolors);
  1567.   ppm_freecolorhash(cht);
  1568.   return chv;
  1569. }
  1570.  
  1571.  
  1572. /****************************************************************************/
  1573. static colorhash_table ppm_computecolorhash(pixels, cols, rows, 
  1574.                         maxcolors, colorsP )
  1575.      pixel** pixels;
  1576.      int cols, rows;
  1577.      int* colorsP;
  1578. {
  1579.   colorhash_table cht;
  1580.   register pixel* pP;
  1581.   colorhist_list chl;
  1582.   int col, row, hash;
  1583.  
  1584.   cht = ppm_alloccolorhash( );
  1585.   *colorsP = 0;
  1586.  
  1587.   /* Go through the entire image, building a hash table of colors. */
  1588.   for (row=0; row<rows; row++)
  1589.     for (col=0, pP=pixels[row];  col<cols;  col++, pP++) {
  1590.       hash = ppm_hashpixel(*pP);
  1591.  
  1592.       for (chl = cht[hash]; chl != (colorhist_list) 0; chl = chl->next)
  1593.     if (PPM_EQUAL(chl->ch.color, *pP)) break;
  1594.       
  1595.       if (chl != (colorhist_list) 0) ++(chl->ch.value);
  1596.       else {
  1597.     if ((*colorsP)++ > maxcolors) {
  1598.       ppm_freecolorhash(cht);
  1599.       return (colorhash_table) 0;
  1600.     }
  1601.     
  1602.     chl = (colorhist_list) malloc(sizeof(struct colorhist_list_item));
  1603.     if (!chl) FatalError("ran out of memory computing hash table");
  1604.  
  1605.     chl->ch.color = *pP;
  1606.     chl->ch.value = 1;
  1607.     chl->next = cht[hash];
  1608.     cht[hash] = chl;
  1609.       }
  1610.     }
  1611.   
  1612.   return cht;
  1613. }
  1614.  
  1615.  
  1616. /****************************************************************************/
  1617. static colorhash_table ppm_alloccolorhash()
  1618. {
  1619.   colorhash_table cht;
  1620.   int i;
  1621.  
  1622.   cht = (colorhash_table) malloc( HASH_SIZE * sizeof(colorhist_list) );
  1623.   if (!cht) FatalError("ran out of memory allocating hash table");
  1624.  
  1625.   for (i=0; i<HASH_SIZE; i++ )
  1626.     cht[i] = (colorhist_list) 0;
  1627.  
  1628.   return cht;
  1629. }
  1630.  
  1631.  
  1632. /****************************************************************************/
  1633. static void ppm_addtocolorhash( cht, colorP, value )
  1634.      colorhash_table cht;
  1635.      pixel* colorP;
  1636.      int value;
  1637. {
  1638.   int hash;
  1639.   colorhist_list chl;
  1640.  
  1641.   hash = ppm_hashpixel( *colorP );
  1642.   chl = (colorhist_list) malloc( sizeof(struct colorhist_list_item) );
  1643.   if (!chl) FatalError("ran out of memory adding to hash table");
  1644.  
  1645.   chl->ch.color = *colorP;
  1646.   chl->ch.value = value;
  1647.   chl->next = cht[hash];
  1648.   cht[hash] = chl;
  1649. }
  1650.  
  1651.  
  1652. /****************************************************************************/
  1653. static colorhist_vector ppm_colorhashtocolorhist( cht, maxcolors )
  1654.      colorhash_table cht;
  1655.      int maxcolors;
  1656. {
  1657.   colorhist_vector chv;
  1658.   colorhist_list chl;
  1659.   int i, j;
  1660.  
  1661.   /* Now collate the hash table into a simple colorhist array. */
  1662.   chv = (colorhist_vector) malloc( maxcolors * sizeof(struct colorhist_item) );
  1663.  
  1664.   /* (Leave room for expansion by caller.) */
  1665.   if (!chv) FatalError("ran out of memory generating histogram");
  1666.  
  1667.   /* Loop through the hash table. */
  1668.   j = 0;
  1669.   for (i=0; i<HASH_SIZE; i++)
  1670.     for (chl = cht[i];  chl != (colorhist_list) 0;  chl = chl->next) {
  1671.       /* Add the new entry. */
  1672.       chv[j] = chl->ch;
  1673.       ++j;
  1674.     }
  1675.  
  1676.   return chv;
  1677. }
  1678.  
  1679.  
  1680. /****************************************************************************/
  1681. static int ppm_lookupcolor( cht, colorP )
  1682.      colorhash_table cht;
  1683.      pixel* colorP;
  1684. {
  1685.   int hash;
  1686.   colorhist_list chl;
  1687.  
  1688.   hash = ppm_hashpixel( *colorP );
  1689.   for (chl = cht[hash];  chl != (colorhist_list) 0;  chl = chl->next)
  1690.     if (PPM_EQUAL(chl->ch.color, *colorP)) return chl->ch.value;
  1691.  
  1692.   return -1;
  1693. }
  1694.  
  1695.  
  1696. /****************************************************************************/
  1697. static void ppm_freecolorhist( chv )
  1698.      colorhist_vector chv;
  1699. {
  1700.   free( (char*) chv );
  1701. }
  1702.  
  1703.  
  1704. /****************************************************************************/
  1705. static void ppm_freecolorhash( cht )
  1706.      colorhash_table cht;
  1707. {
  1708.   int i;
  1709.   colorhist_list chl, chlnext;
  1710.  
  1711.   for (i=0; i<HASH_SIZE; i++)
  1712.     for (chl = cht[i];  chl != (colorhist_list) 0; chl = chlnext) {
  1713.       chlnext = chl->next;
  1714.       free( (char*) chl );
  1715.     }
  1716.  
  1717.   free( (char*) cht );
  1718. }
  1719.