#define PROG_NAME "ppmoquant"
#define PROG_DESC "quantize two PPM files into a partially transparent PPM"
#define PROG_VERS "1.0"

/* Copyright © 1989, 1991 by Jef Poskanzer.
** See the copyright, authorship, and warranty notice at end of file.
** Last edited on 2008-07-06 22:53:23 by stolfi
*/

/* TO DO: !!! Currently has syntax errors -- Update calls and fix inconsistencies. */

#define PROG_HELP \
  PROG_NAME " \\\n" \
  "  [ -fs | -nofs ] \\\n" \
  "  -transparent COLOR  -map PPMFILE \\\n" \
  "  PPMFILE0 PPMFILE1 > PPMFILE"

/* Similar to "ppmquant" but takes two pixmaps of an object over 
  black and white backgrounds, respectively, to deduce transparency. */

#include <stdlib.h>

#include <jspnm.h>
#include <jspnm_image.h>
#include <jsppm_xhist.h>
#include <jsppm_medcutx.h>

#include <ppmoquant.h>

int main(int argc, char* argv[])
  {
    int floyd;
    matchfn match;
    char *bimgname, *wimgname;
    char *cmapname;
    rgb_pixel_t transp;
    int transpmaxval;

    pnm_image_t *cmap, *bimg, *wimg;
    ppm_xhist_vector cm;
    int ncolors;
    int mapmaxval;

    parse_options(
      argc, argv, 
      &floyd, &match,
      &transp, &transpmaxval, 
      &cmapname,
      &bimgname, &wimgname
    );

    /* Read the black-bg and white-bg images: */
    
    bimg = pnm_image_read((bimgname == NULL ? "-" : bimgname));
    wimg = pnm_image_read((wimgname == NULL ? "-" : wimgname));
    
    if ((bimg->rows != wimg->rows) || (bimg->cols != wimg->cols))
      pm_error("images have incompatible dimensions");
    if (bimg->maxval != wimg->maxval)
      pm_error("images have different maxvals");
      
    /* Read the color map: */
    cmap = pnm_image_read(cmapname);
    mapmaxval = cmap->maxval;
    cm = collect_colors(cmap, &ncolors);
    
    /* Fix the "transparent" color to have the same maxval as the colormap: */
    transp = adjust_transp_color(transp, transpmaxval, mapmaxval);
    if (find_color(transp, cm, ncolors) != NULL)
      pm_error("\"-transparent\" color occurs in colormap");
    
    if (floyd)
      { floyd_quantize_pnm_images(bimg, wimg, cm, ncolors, transp, mapmaxval, match); }
    else                                                   
      { plain_quantize_pnm_images(bimg, wimg, cm, ncolors, transp, mapmaxval, match); }

    pnm_image_write("-", bimg, 0);
    exit(0);
  }

static void 
parse_options(
    int *argc, char **argv, 
    int *floydp,
    matchfn *matchp,
    rgb_pixel_t *transpp,
    int *transpmaxvalp,
    char **cmapnamep,
    char **bimgnamep,
    char **wimgnamep
  )
  {
    int argn;
    int transpgiven;
    char* usage = 
      PROG_NAME " version " PROG_VERS ", usage:\n" 
      PROG_HELP;

    ppm_init(argc, argv);
    argn = 1;
    (*floydp) = 1;
    (*matchp) = rgb_match_range;
    transpgiven = 0;
    (*transpmaxvalp) = 0;
    (*cmapnamep) = NULL;
    (*bimgnamep) = NULL;
    (*wimgnamep) = NULL;
    
    while (argn < (*argc) && argv[argn][0] == '-' && argv[argn][1] != '\0')
      { char *opt = argv[argn];
        if (
            pm_keymatch(opt, "-fs", 2) ||
	    pm_keymatch(opt, "-floyd", 2)
          )
          (*floydp) = 1;
        else if (
            pm_keymatch(opt, "-nofs", 2) ||
	    pm_keymatch(opt, "-nofloyd", 2)
          )
          (*floydp) = 0;
        else if (pm_keymatch(opt, "-map", 2))
          { ++argn;
            if (argn >= (*argc)) pm_usage(usage);
            (*cmapnamep) = argv[argn]; 
          }
        else if (pm_keymatch(opt, "-transparent", 2))
          { ++argn;
 	    if (argn >= (*argc)) pm_usage(usage);
            transpgiven = 1;
            (*transpp) = my_parsecolor(argv[argn], transpmaxvalp);
          }
        else pm_usage(usage);
        ++argn;
      }

    /* Colormapis mandatory: */
    if ((*cmapnamep) == NULL) pm_usage(usage);
    
    /* Transparent color is mandatory: */
    if (transpgiven == 0) pm_usage(usage);
    
    /* Black-background and white-background images are mandatory: */
    if (argn >= (*argc) - 1) pm_usage(usage);
    (*bimgnamep) = argv[argn]; argn++;
    (*wimgnamep) = argv[argn]; argn++;
    
    if (argn != (*argc)) pm_usage(usage);
  }

static rgb_pixel_t
my_parsecolor(char *spec, int *maxvalp)
  { int r, g, b;
    rgb_pixel_t p;
    if ( (*spec) == '(' )
      { /* Absolute RGB, decimal: */
        if (sscanf(spec, "( %d, %d, %d )", &r, &g, &b) != 3)
          pm_error("bad colorspec - expecting \"(R,G,B)\"");
        if ((r < 0) || (g < 0) || (b < 0)) 
          pm_error("bad \"-transparent\" color value");
        (*maxvalp) = 0;
        PPM_ASSIGN(p, r, g, b);
      }
    else
      { p = ppm_parsecolor(spec, PPM_MAXMAXVAL);
        (*maxvalp) = PPM_MAXMAXVAL;
      }
    return(p);
  }

static rgb_pixel_t
adjust_transp_color(rgb_pixel_t p, int oldmaxval, int newmaxval)
  { 
    int r, g, b;
    rgb_pixel_t q;
    if ((oldmaxval != 0) && (oldmaxval != newmaxval))
      {  
        float s = ((float)newmaxval)/((float)oldmaxval);
        r = s*PPM_GETR(p) + 0.5;
        g = s*PPM_GETG(p) + 0.5;
        b = s*PPM_GETB(p) + 0.5;
        pm_message(
          "\"-transparent\" rgb_pixel_t scaled to (%d,%d,%d) to match colormap",
          r, g, b
        );
      }
    else
      { r = PPM_GETR(p);
        g = PPM_GETG(p);
        b = PPM_GETB(p);
      }
    if ((r > newmaxval) || (g > newmaxval) || (b > newmaxval)) 
      { pm_message("\"-transparent\" rgb_pixel_t = (%d,%d,%d), newmaxval = %d",
          r, g, b, newmaxval
        );
        pm_error("bad \"-transparent\" color value");
      }
    PPM_ASSIGN(q, r, g, b);
    return(q);
  }

static ppm_xhist_vector
collect_colors(pnm_image_t *map, int *ncolorsp)
  { ppm_xhist_vector cm;
    if (map->cols == 0 || map->rows == 0)
      { pm_error("null colormap??"); }
    cm = ppm_computecolorhist(
      map->pix, map->cols, map->rows, MAXCOLORS, ncolorsp
    );
    if (cm == (ppm_xhist_vector) NULL)
      pm_error("too many colors in colormap!");
    ppm_freearray(map->pix, map->rows);
      pm_message("%d colors found in colormap", (*ncolorsp));
    return(cm);
  }

static rgb_pixel_t* 
find_color(
    rgb_pixel_t p,
    ppm_xhist_vector cm, 
    int ncolors
  )
  {
    register int i;
    register int pr = PPM_GETR(p);
    register int pg = PPM_GETG(p);
    register int pb = PPM_GETB(p);
    for (i = 0; i < ncolors; ++i)
      { rgb_pixel_t *qq = &(cm[i].color);
        int r = PPM_GETR(*qq);
        int g = PPM_GETG(*qq);
        int b = PPM_GETB(*qq);
        if ((r == pr) && (g == pg) && (b == pb)) return (qq);
      }
    return(NULL);
  }

static void
floyd_quantize_pnm_images(
    pnm_image_t *bimg, 
    pnm_image_t *wimg, 
    ppm_xhist_vector cm, 
    int ncolors,
    rgb_pixel_t transp,
    int mapmaxval,
    matchfn match
  )
  {
#define FS_SCALE 1024
#define FS_WTA 7
#define FS_WTB 3
#define FS_WTC 5
#define FS_WTD 1
    fs_errors be, we;
    
    int row;
    register int col, limitcol;
    register rgb_pixel_t *bp, *wp, *qq;
    int maxval = bimg->maxval;
    long maxcor = 2*maxval;
    long mincor = -1*maxval;
    int fs_direction = 1;
    /* int addtohash = 1; */
    /* colorhash_table cht = ppm_alloccolorhash(); */
    register long br, bg, bb, wr, wg, wb;
    float coef = FS_SCALE*((float)maxval)/((float)mapmaxval);

    floyd_alloc_error_vectors(&be, bimg->cols);
    floyd_alloc_error_vectors(&we, wimg->cols);

    /* Initialize Floyd-Steinberg error vectors with randoms in [-1..+1]. */
    srand((int) 46157);
    for (col = 0; col < bimg->cols + 2; ++col)
      {
        be.thisr[col] = (rand() % (FS_SCALE * 2)) - FS_SCALE;
        be.thisg[col] = (rand() % (FS_SCALE * 2)) - FS_SCALE;
        be.thisb[col] = (rand() % (FS_SCALE * 2)) - FS_SCALE;

        we.thisr[col] = be.thisr[col];
        we.thisg[col] = be.thisg[col];
        we.thisb[col] = be.thisb[col];

      }
    for (row = 0; row < bimg->rows; ++row)
      { for (col = 0; col < bimg->cols + 2; ++col)
          { be.nextr[col] = be.nextg[col] = be.nextb[col] = 0;
            we.nextr[col] = we.nextg[col] = we.nextb[col] = 0;
          }
        if (fs_direction)
          { col = 0; limitcol = bimg->cols;
            bp = &(bimg->pix[row][0]);
            wp = &(wimg->pix[row][0]);
          }
        else
          { col = bimg->cols - 1; limitcol = -1;
            bp = &(bimg->pix[row][col]);
            wp = &(wimg->pix[row][col]);
          }
        do
          { /* Use Floyd-Steinberg errors to adjust actual color. */
            br = PPM_GETR(*bp) + be.thisr[col + 1] / FS_SCALE;
            bg = PPM_GETG(*bp) + be.thisg[col + 1] / FS_SCALE;
            bb = PPM_GETB(*bp) + be.thisb[col + 1] / FS_SCALE;
            
            if (br < mincor) br = mincor; else if (br > maxcor) br = maxcor;
            if (bg < mincor) bg = mincor; else if (bg > maxcor) bg = maxcor;
            if (bb < mincor) bb = mincor; else if (bb > maxcor) bb = maxcor;

            wr = PPM_GETR(*wp) + we.thisr[col + 1] / FS_SCALE;
            wg = PPM_GETG(*wp) + we.thisg[col + 1] / FS_SCALE;
            wb = PPM_GETB(*wp) + we.thisb[col + 1] / FS_SCALE;

            if (wr < mincor) wr = mincor; else if (wr > maxcor) wr = maxcor;
            if (wg < mincor) wg = mincor; else if (wg > maxcor) wg = maxcor;
            if (wb < mincor) wb = mincor; else if (wb > maxcor) wb = maxcor;

            /* Choose replacement color */
            qq = choose_color(
              br, bg, bb, wr, wg, wb, maxval,
              cm, ncolors, mapmaxval,
              match
            );
            /* Replace rgb_pixel_t in black image, and propagate error: */
            if (qq == NULL)
              { /* Best match is transparent: */
                (*bp) = transp;
                floyd_propagate_error_transp(
                  br, bg, bb,  wr, wg, wb,  
                  maxval,
                  &be, &we, col, fs_direction
                );
              }
            else
              { /* Best match is opaque: */
                (*bp) = (*qq);
                floyd_propagate_error_opaque(
                  br, bg, bb,  wr, wg, wb,
                  qq, coef,
                  &be, &we, col, fs_direction
                );
              }

            if (fs_direction)
              { ++col; ++bp; ++wp; }
            else
              { --col; --bp; --wp; }
            }
        while (col != limitcol);
        floyd_swap_error_vectors(&be);
        floyd_swap_error_vectors(&we);
        
        fs_direction = ! fs_direction;
      }
    bimg->maxval = mapmaxval;
  }

static void
floyd_alloc_error_vectors(fs_errors *ep, int cols)
  {
    ep->thisr = (long*) pm_allocrow(cols + 2, sizeof(long));
    ep->nextr = (long*) pm_allocrow(cols + 2, sizeof(long));
    ep->thisg = (long*) pm_allocrow(cols + 2, sizeof(long));
    ep->nextg = (long*) pm_allocrow(cols + 2, sizeof(long));
    ep->thisb = (long*) pm_allocrow(cols + 2, sizeof(long));
    ep->nextb = (long*) pm_allocrow(cols + 2, sizeof(long));
  }

static void
floyd_propagate_error_opaque(
    long br, long bg, long bb,
    long wr, long wg, long wb,
    rgb_pixel_t *newp, float coef,
    fs_errors *bep, fs_errors *wep,
    int col,
    int fs_direction
  )
  {
    long err, new;
    if (fs_direction)
      {
        new = (long)(coef*PPM_GETR((*newp))+0.5);
        err = (br * FS_SCALE - new);
        bep->thisr[col + 2] += (err * FS_WTA) / 16;
        bep->nextr[col    ] += (err * FS_WTB) / 16;
        bep->nextr[col + 1] += (err * FS_WTC) / 16;
        bep->nextr[col + 2] += (err * FS_WTD) / 16;
        err = (wr * FS_SCALE - new);
        wep->thisr[col + 2] += (err * FS_WTA) / 16;
        wep->nextr[col    ] += (err * FS_WTB) / 16;
        wep->nextr[col + 1] += (err * FS_WTC) / 16;
        wep->nextr[col + 2] += (err * FS_WTD) / 16;

        new = (long)(coef*PPM_GETG((*newp))+0.5);
        err = (bg * FS_SCALE - new);
        bep->thisg[col + 2] += (err * FS_WTA) / 16;
        bep->nextg[col    ] += (err * FS_WTB) / 16;
        bep->nextg[col + 1] += (err * FS_WTC) / 16;
        bep->nextg[col + 2] += (err * FS_WTD) / 16;
        err = (wg * FS_SCALE - new);
        wep->thisg[col + 2] += (err * FS_WTA) / 16;
        wep->nextg[col    ] += (err * FS_WTB) / 16;
        wep->nextg[col + 1] += (err * FS_WTC) / 16;
        wep->nextg[col + 2] += (err * FS_WTD) / 16;

        new = (long)(coef*PPM_GETB((*newp))+0.5);
        err = (bb * FS_SCALE - new);
        bep->thisb[col + 2] += (err * FS_WTA) / 16;
        bep->nextb[col    ] += (err * FS_WTB) / 16;
        bep->nextb[col + 1] += (err * FS_WTC) / 16;
        bep->nextb[col + 2] += (err * FS_WTD) / 16;
        err = (wb * FS_SCALE - new);
        wep->thisb[col + 2] += (err * FS_WTA) / 16;
        wep->nextb[col    ] += (err * FS_WTB) / 16;
        wep->nextb[col + 1] += (err * FS_WTC) / 16;
        wep->nextb[col + 2] += (err * FS_WTD) / 16;
      }
    else
      {
        new = (long)(coef*PPM_GETR((*newp))+0.5);
        err = (br * FS_SCALE - new);
        bep->thisr[col    ] += (err * FS_WTA) / 16;
        bep->nextr[col + 2] += (err * FS_WTB) / 16;
        bep->nextr[col + 1] += (err * FS_WTC) / 16;
        bep->nextr[col    ] += (err * FS_WTD) / 16;
        err = (wr * FS_SCALE - new);
        wep->thisr[col    ] += (err * FS_WTA) / 16;
        wep->nextr[col + 2] += (err * FS_WTB) / 16;
        wep->nextr[col + 1] += (err * FS_WTC) / 16;
        wep->nextr[col    ] += (err * FS_WTD) / 16;
        
        new = (long)(coef*PPM_GETG((*newp))+0.5);
        err = (bg * FS_SCALE - new);
        bep->thisg[col    ] += (err * FS_WTA) / 16;
        bep->nextg[col + 2] += (err * FS_WTB) / 16;
        bep->nextg[col + 1] += (err * FS_WTC) / 16;
        bep->nextg[col    ] += (err * FS_WTD) / 16;
        err = (wg * FS_SCALE - new);
        wep->thisg[col    ] += (err * FS_WTA) / 16;
        wep->nextg[col + 2] += (err * FS_WTB) / 16;
        wep->nextg[col + 1] += (err * FS_WTC) / 16;
        wep->nextg[col    ] += (err * FS_WTD) / 16;
        
        new = (long)(coef*PPM_GETB((*newp))+0.5);
        err = (bb * FS_SCALE - new);
        bep->thisb[col    ] += (err * FS_WTA) / 16;
        bep->nextb[col + 2] += (err * FS_WTB) / 16;
        bep->nextb[col + 1] += (err * FS_WTC) / 16;
        bep->nextb[col    ] += (err * FS_WTD) / 16;
        err = (wb * FS_SCALE - new);
        wep->thisb[col    ] += (err * FS_WTA) / 16;
        wep->nextb[col + 2] += (err * FS_WTB) / 16;
        wep->nextb[col + 1] += (err * FS_WTC) / 16;
        wep->nextb[col    ] += (err * FS_WTD) / 16;
      }
  }

static void 
floyd_propagate_error_transp(
    long br, long bg, long bb,
    long wr, long wg, long wb,
    int maxval,
    fs_errors *bep, fs_errors *wep,
    int col,
    int fs_direction
  )
  {
    long err;
    if (fs_direction)
      {
        err = br * FS_SCALE;
        bep->thisr[col + 2] += (err * FS_WTA) / 16;
        bep->nextr[col    ] += (err * FS_WTB) / 16;
        bep->nextr[col + 1] += (err * FS_WTC) / 16;
        bep->nextr[col + 2] += (err * FS_WTD) / 16;
        err = (wr - maxval) * FS_SCALE;
        wep->thisr[col + 2] += (err * FS_WTA) / 16;
        wep->nextr[col    ] += (err * FS_WTB) / 16;
        wep->nextr[col + 1] += (err * FS_WTC) / 16;
        wep->nextr[col + 2] += (err * FS_WTD) / 16;

        err = bg * FS_SCALE;
        bep->thisg[col + 2] += (err * FS_WTA) / 16;
        bep->nextg[col    ] += (err * FS_WTB) / 16;
        bep->nextg[col + 1] += (err * FS_WTC) / 16;
        bep->nextg[col + 2] += (err * FS_WTD) / 16;
        err = (wg - maxval) * FS_SCALE;
        wep->thisg[col + 2] += (err * FS_WTA) / 16;
        wep->nextg[col    ] += (err * FS_WTB) / 16;
        wep->nextg[col + 1] += (err * FS_WTC) / 16;
        wep->nextg[col + 2] += (err * FS_WTD) / 16;

        err = bb * FS_SCALE;
        bep->thisb[col + 2] += (err * FS_WTA) / 16;
        bep->nextb[col    ] += (err * FS_WTB) / 16;
        bep->nextb[col + 1] += (err * FS_WTC) / 16;
        bep->nextb[col + 2] += (err * FS_WTD) / 16;
        err = (wb - maxval) * FS_SCALE;
        wep->thisb[col + 2] += (err * FS_WTA) / 16;
        wep->nextb[col    ] += (err * FS_WTB) / 16;
        wep->nextb[col + 1] += (err * FS_WTC) / 16;
        wep->nextb[col + 2] += (err * FS_WTD) / 16;
      }
    else
      {
        err = br * FS_SCALE;
        bep->thisr[col    ] += (err * FS_WTA) / 16;
        bep->nextr[col + 2] += (err * FS_WTB) / 16;
        bep->nextr[col + 1] += (err * FS_WTC) / 16;
        bep->nextr[col    ] += (err * FS_WTD) / 16;
        err = (wr - maxval) * FS_SCALE;
        wep->thisr[col    ] += (err * FS_WTA) / 16;
        wep->nextr[col + 2] += (err * FS_WTB) / 16;
        wep->nextr[col + 1] += (err * FS_WTC) / 16;
        wep->nextr[col    ] += (err * FS_WTD) / 16;
        
        err = bg * FS_SCALE;
        bep->thisg[col    ] += (err * FS_WTA) / 16;
        bep->nextg[col + 2] += (err * FS_WTB) / 16;
        bep->nextg[col + 1] += (err * FS_WTC) / 16;
        bep->nextg[col    ] += (err * FS_WTD) / 16;
        err = (wg - maxval) * FS_SCALE;
        wep->thisg[col    ] += (err * FS_WTA) / 16;
        wep->nextg[col + 2] += (err * FS_WTB) / 16;
        wep->nextg[col + 1] += (err * FS_WTC) / 16;
        wep->nextg[col    ] += (err * FS_WTD) / 16;
        
        err = bb * FS_SCALE;
        bep->thisb[col    ] += (err * FS_WTA) / 16;
        bep->nextb[col + 2] += (err * FS_WTB) / 16;
        bep->nextb[col + 1] += (err * FS_WTC) / 16;
        bep->nextb[col    ] += (err * FS_WTD) / 16;
        err = (wb - maxval) * FS_SCALE;
        wep->thisb[col    ] += (err * FS_WTA) / 16;
        wep->nextb[col + 2] += (err * FS_WTB) / 16;
        wep->nextb[col + 1] += (err * FS_WTC) / 16;
        wep->nextb[col    ] += (err * FS_WTD) / 16;
      }
  }
  
static void
floyd_swap_error_vectors(fs_errors *ep)
  {
    long* t;
    t = ep->thisr; ep->thisr = ep->nextr; ep->nextr = t;
    t = ep->thisg; ep->thisg = ep->nextg; ep->nextg = t;
    t = ep->thisb; ep->thisb = ep->nextb; ep->nextb = t;
  }

static void
plain_quantize_pnm_images(
    pnm_image_t *bimg, 
    pnm_image_t *wimg, 
    ppm_xhist_vector cm, 
    int ncolors,
    rgb_pixel_t transp,
    int mapmaxval,
    matchfn match
  )
  {
    int row;
    register int col, limitcol;
    register rgb_pixel_t *bp, *wp, *qq;
    long br, bg, bb, wr, wg, wb;
    int maxval = bimg->maxval;
    
    for (row = 0; row < bimg->rows; ++row)
      {
        col = 0;
        limitcol = bimg->cols;
        bp = bimg->pix[row];
        wp = wimg->pix[row];

        do
          {
            br = PPM_GETR(*bp);
            bg = PPM_GETG(*bp);
            bb = PPM_GETB(*bp);

            wr = PPM_GETR(*wp);
            wg = PPM_GETG(*wp);
            wb = PPM_GETB(*wp);
                
            /* Choose replacement color* */
            qq = choose_color(
              br, bg, bb, wr, wg, wb, maxval,
              cm, ncolors, mapmaxval,
              match
            );
            /* Replace rgb_pixel_t in black image, and propagate error: */
            if (qq == NULL)
              { (*bp) = transp; }
            else
              { (*bp) = (*qq); }
            ++col; ++bp; ++wp;
          }
        while (col != limitcol);
      }
    bimg->maxval = mapmaxval;
  }

static rgb_pixel_t*
choose_color (
    long br, long bg, long bb,
    long wr, long wg, long wb,
    int maxval,
    ppm_xhist_vector cm, 
    int ncolors,
    int mapmaxval,
    matchfn match
  )
  {
    int ind = match(
      br, bg, bb, wr, wg, wb, maxval,
      cm, ncolors, mapmaxval
    );
    if (ind == -1) return(NULL);
    else return(&(cm[ind].color));
  }

static int
rgb_match_range (
    long br, long bg, long bb,
    long wr, long wg, long wb,
    int maxval,
    ppm_xhist_vector cm, 
    int ncolors,
    int mapmaxval
  )
  { 
#define RGBMATCH_SCALE 4
    register int i;
    register int ind = -1;
    register long tr, tg, tb;
    register long d, dist, newdist;
        
    /* Scale the rgb_pixel_t coordinates to match the map's maxval: */
    float s = RGBMATCH_SCALE * ((float)mapmaxval)/((float)maxval);
    
    long fbr = s*br + 0.5;
    long fbg = s*bg + 0.5;
    long fbb = s*bb + 0.5;

    long fwr = s*wr + 0.5;
    long fwg = s*wg + 0.5;
    long fwb = s*wb + 0.5;
    
    long fmx = s*maxval + 0.5;
    
    /* 
      Initialize dist with the distance between the transparent
      color range and the given color range:
    */
    dist = 0;
    d = fbr; if (d < 0) d = -d; if (d > dist) dist = d;
    d = fbg; if (d < 0) d = -d; if (d > dist) dist = d;
    d = fbb; if (d < 0) d = -d; if (d > dist) dist = d;
    
    d = fwr - fmx; if (d < 0) d = -d; if (d > dist) dist = d;
    d = fwg - fmx; if (d < 0) d = -d; if (d > dist) dist = d;
    d = fwb - fmx; if (d < 0) d = -d; if (d > dist) dist = d;

    /* Now see if an opaque color would work better: */
    for (i = 0; ((dist > 0) && (i < ncolors)); ++i)
      { rgb_pixel_t *qq = &(cm[i].color);
        tr = RGBMATCH_SCALE * PPM_GETR(*qq);
        tg = RGBMATCH_SCALE * PPM_GETG(*qq);
        tb = RGBMATCH_SCALE * PPM_GETB(*qq);
        /* 
          Find the maximum distance newdist between this color and the corners 
          of the given interval. If it is larger than dist, skip it.
        */
        newdist = 0;
        d = tr - fbr; if (d < 0) d = -d; 
        if (d > newdist) { if (d >= dist) goto skipit; newdist = d; }
        d = tg - fbg; if (d < 0) d = -d;
        if (d > newdist) { if (d >= dist) goto skipit; newdist = d; }
        d = tb - fbb; if (d < 0) d = -d;
        if (d > newdist) { if (d >= dist) goto skipit; newdist = d; }
            
        d = tr - fwr; if (d < 0) d = -d; 
        if (d > newdist) { if (d >= dist) goto skipit; newdist = d; }
        d = tg - fwg; if (d < 0) d = -d;
        if (d > newdist) { if (d >= dist) goto skipit; newdist = d; }
        d = tb - fwb; if (d < 0) d = -d;
        if (d > newdist) { if (d >= dist) goto skipit; newdist = d; }
            
        /* If we got here, newdist is smaller than dist: */
        ind = i; dist = newdist;
      skipit:
        /* OK */;
      }
    return (ind);
  }

/* Modification history:
**
**   18/may/1996: Opacity logic added by J. Stolfi
**   1989-1991: Original "ppmquant" by Jef Poskanzer.
*/

/* Copyright © 1989, 1991 by Jef Poskanzer.
**
** Permission to use, copy, modify, and distribute this software and its
** documentation for any purpose and without fee is hereby granted, provided
** that the above copyright notice appear in all copies and that both that
** copyright notice and this permission notice appear in supporting
** documentation.  This software is provided "as is" without express or
** implied warranty.
*/