/* Last edited on 2009-01-07 01:31:23 by stolfi */

    /* Fix the "transparent" color to have the same maxval as the colormap: */
    o->transpcolor = adjust_transp_color(o->transpcolor, o->transpmaxval, mapmaxval);

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

ppm_pixel_t adjust_transp_color(ppm_pixel_t p, int oldmaxval, int newmaxval);
  /* Rescales the pixel p from [0..oldmaxval] to [0..newmaxval],
    unless oldmaxval = 0 (when p itself is returned). */

ppm_pixel_t adjust_transp_color(ppm_pixel_t p, int oldmaxval, int newmaxval)
  { 
    int r, g, b;
    ppm_pixel_t q;
    if ((oldmaxval != 0) && (oldmaxval != newmaxval))
      { double s = ((double)newmaxval)/((double)oldmaxval);
        q.c[0] = (int)floor(s*p.c[0] + 0.5);
        q.c[1] = (int)floor(s*p.c[1] + 0.5);
        q.c[2] = (int)floor(s*p.c[2] + 0.5);
        pnm_message
         ( "\"-transparent\" ppm_pixel_t scaled to (%d,%d,%d) to match colormap", q.c[0], q.c[1], q.c[2]);
      }
    else
      { q = p; }
    if ((q.c[0] > newmaxval) || (q.c[1] > newmaxval) || (q.c[2] > newmaxval)) 
      { pnm_message("\"-transparent\" ppm_pixel_t = (%d,%d,%d), newmaxval = %d", q.c[0], q.c[1], q.c[2], newmaxval);
        pnm_error(pp, "bad \"-transparent\" color value");
      }
    return q;
  }

/* ---------------------------------------------------------------------- */

/* Was opq-parts.c */
static void
opq_quantize_image(
    image *img, 
    image *opq,
    colorhist_vector cm, 
    int newcolors,
    matchfn match
  )
  /*
    Replaces each pixel in the image by a pixel from the given colormap.
  */
  {
#define FS_SCALE 1024

    typedef struct 
      { long *thisc;
        long *nextc;
        long *thisa;
        long *nexta;
      } channel_errors;
      
    channel_errors err[3];
    long* temperr;
    int row, k;
    register int col, limitcol;
    register pixel *pc, *pa, *qq;
    int maxc = img->maxval;
    int maxa = opq->maxval;
    int fs_direction = 1;
    int addtohash = 1;
    colorhash_table cht = ppm_alloccolorhash();
    register long sr, sg, sb, err;
    pixel coal, dork;

    PPM_ASSIGN(coal, 0, 0, 0);
    PPM_ASSIGN(dork, 1, 1, 1);
    
    if (img->rows != opq->cols) pm_error("image and mask have different shapes");
    if (img->cols != opq->cols) pm_error("image and mask have different shapes");

    /* Initialize Floyd-Steinberg error vectors with randoms in [-1..+1]. */
    srandom((int) 46157);
    for (k=0; k<3; k++)
      { err[k].thisc = (long*) pm_allocrow(img->cols + 2, sizeof(long));
        err[k].nextc = (long*) pm_allocrow(img->cols + 2, sizeof(long));
        err[k].thisa = (long*) pm_allocrow(img->cols + 2, sizeof(long));
        err[k].nexta = (long*) pm_allocrow(img->cols + 2, sizeof(long));
        for (col = 0; col < img->cols + 2; ++col)
          {
            err[k].thisc[col] = (random() % (FS_SCALE * 2)) - FS_SCALE;
            err[k].thisa[col] = (random() % (FS_SCALE * 2)) - FS_SCALE;
          }
      }
    
    for (row = 0; row < img->rows; ++row)
      { for (col = 0; col < img->cols + 2; ++col)
          for (k=0; k<3; k++)
            { err[k].nextc[col] = err[k].nexta[col] = 0; }
        if (fs_direction)
          { col = 0; limitcol = img->cols;
            pc = img->pixels[row];
            pa = opq->pixels[row];
          }
        else
          { col = img->cols - 1; limitcol = -1;
            pc = &(img->pixels[row][col]);
            pa = &(opq->pixels[row][col]);
          }
        do
          { /* Use Floyd-Steinberg errors to adjust actual color. */
            for (k=0, k<3; k++)
              { 
            sr = PPM_GETR(*pc) + err[0].thisc[col + 1] / FS_SCALE;
            sg = PPM_GETG(*pc) + err[1].thisc[col + 1] / FS_SCALE;
            sb = PPM_GETB(*pc) + err[2].thisc[col + 1] / FS_SCALE;
            if (sr < 0) sr = 0; else if (sr > maxc) sr = maxc;
            if (sg < 0) sg = 0; else if (sg > maxc) sg = maxc;
            if (sb < 0) sb = 0; else if (sb > maxc) sb = maxc;
            PPM_ASSIGN(*pc, sr, sg, sb);
            /* Ditto for transparency */
            sr = PPM_GETR(*pa) + err[0].thisa[col + 1] / FS_SCALE;
            sg = PPM_GETG(*pa) + err[1].thisa[col + 1] / FS_SCALE;
            sb = PPM_GETB(*pa) + err[2].thisa[col + 1] / FS_SCALE;
            if (sr < 0) sr = 0; else if (sr > maxa) sr = maxa;
            if (sg < 0) sg = 0; else if (sg > maxa) sg = maxa;
            if (sb < 0) sb = 0; else if (sb > maxa) sb = maxa;
            PPM_ASSIGN(*pa, sr, sg, sb);
            /* Find best match */
            qq = opq_choose_color(pc, pa, cm, newcolors, match, cht, &addtohash);
            /* Propagate Floyd-Steinberg error terms. */
            if (fs_direction)
              {
                err = (sr - (long) PPM_GETR((*qq))) * FS_SCALE;
                thisrerr[col + 2] += (err * 7) / 16;
                nextrerr[col    ] += (err * 3) / 16;
                nextrerr[col + 1] += (err * 5) / 16;
                nextrerr[col + 2] += (err    ) / 16;
                err = (sg - (long) PPM_GETG((*qq))) * FS_SCALE;
                thisgerr[col + 2] += (err * 7) / 16;
                nextgerr[col    ] += (err * 3) / 16;
                nextgerr[col + 1] += (err * 5) / 16;
                nextgerr[col + 2] += (err    ) / 16;
                err = (sb - (long) PPM_GETB((*qq))) * FS_SCALE;
                thisberr[col + 2] += (err * 7) / 16;
                nextberr[col    ] += (err * 3) / 16;
                nextberr[col + 1] += (err * 5) / 16;
                nextberr[col + 2] += (err    ) / 16;
              }
            else
              {
                err = (sr - (long) PPM_GETR((*qq))) * FS_SCALE;
                thisrerr[col    ] += (err * 7) / 16;
                nextrerr[col + 2] += (err * 3) / 16;
                nextrerr[col + 1] += (err * 5) / 16;
                nextrerr[col    ] += (err    ) / 16;
                err = (sg - (long) PPM_GETG((*qq))) * FS_SCALE;
                thisgerr[col    ] += (err * 7) / 16;
                nextgerr[col + 2] += (err * 3) / 16;
                nextgerr[col + 1] += (err * 5) / 16;
                nextgerr[col    ] += (err    ) / 16;
                err = (sb - (long) PPM_GETB((*qq))) * FS_SCALE;
                thisberr[col    ] += (err * 7) / 16;
                nextberr[col + 2] += (err * 3) / 16;
                nextberr[col + 1] += (err * 5) / 16;
                nextberr[col    ] += (err    ) / 16;
              }
            /* Replace pixel in image: */
            (*pp) = (*qq);

            if (fs_direction)
              { ++col; ++pp; }
            else
              { --col; --pp; }
            }
        while (col != limitcol);

        temperr = thisrerr;
        thisrerr = nextrerr;
        nextrerr = temperr;
        temperr = thisgerr;
        thisgerr = nextgerr;
        nextgerr = temperr;
        temperr = thisberr;
        thisberr = nextberr;
        nextberr = temperr;
        fs_direction = ! fs_direction;
      }
  }

/* End of opq-parts.c */
======================================================================


/* This is a loop to enumerate all corners of a (possibly degenerate) range. */
/* Should be useful for YUV range match, not for RGB match. */

        /* 
          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;
        zr = fbr; zg = fbg; zb = fbb;
        while(1)
          { 
            d = tr - zr; if (d < 0) d = -d; 
            if (d > newdist) { if (d > dist) goto skipit; newdist = d; }
            d = tg - zg; if (d < 0) d = -d;
            if (d > newdist) { if (d > dist) goto skipit; newdist = d; }
            d = tb - zb; 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;
            
            /* Get next corner: */
            if (zb != fwb) 
              { zb = fwb; }
            else if (zg != fwg)
              { zg = fwg; zb = fbb; }
            else if (zr != fwr)
              { zr = fwr; zg = fbg; zb = fbb; }
            else
              { break; }
          }
      skipit:
        /* OK */

    /* 
      This implementation assumes that the extremum of the
      distance function occurs at the corners of the 
      RGB cube.  This assumptions still needs to be proved,
      and is likely to be false for negative-intensity 
      colors under the YUV distance.
    */


static long*
new_long_vector(int n)
  { long *lv = (long*) malloc(n*sizeof(long));
    if (lv == (long*)NULL) pm_error("out of memory");
    return(lv);
  }

    /* Read or build the opacity map: */
    
     if (omapname == NULL)
      { om = choose_opacmap(wimg, &newopacs); }
    else
      { omap = pgm_image_read(omapname);
        om = adjust_opacmap(omap, wimg, &newopacs);
      }
    
        else if (strcmp(opt, "-opacmap") == 0)
          { ++argn;
            if (argn >= (*argc)) pm_usage(usage);
            (*wimgnamep) = argv[argn]; 
          }
        else if (strcmp(opt, "-nopacs") == 0)
          { ++argn;
            if (argn >= (*argc)) pm_usage(usage);
            if (sscanf(argv[argn], "%d", newopacsp) != 1) pm_usage(usage);
            if ((*newopacsp) <= 1)
              { pm_error("number of opacs must be > 1"); }
          }


static pgm_image*
pgm_image_read (char *name)
  { FILE *f = NULL;
    pgm_image *bimg;
    if (name == NULL) pm_error("no file name given");
    else if (strlen(name) == 0) pm_error("empty file name");
    else { f = pm_openr(name); }
    bimg = pgm_image_new();
    bimg->pixels = pgm_readpgm( f, &(bimg->cols), &(bimg->rows), &(bimg->maxval));
    if (f != stdin)
      { pm_close(f); }
    return(bimg);
  }

static grayhist_vector
choose_opacmap(pgm_image *wimg, int *newopacsp)
  /*
    Chooses a opacmap with at most *newopacs entries for *wimg.
    On exit *newopacs may have been reduced.
  */
  { 
    int opacs;
    grayhist_vector gh, om;

    /*
      Attempt to make a histogram of the grays, unclustered.
      If at first we don't succeed, lower wimg->maxval 
      and try again.
    */

    for ( ; ; )
      { gh = pgm_grayhist(
          wimg->pixels,
          wimg->cols, wimg->rows, MAXOPACS, 
          &opacs
        );
        if (gh != (grayhist_vector)NULL) 
          { break; }
        else
          { int newmaxval = wimg->maxval / 2;
            int row;
            register int col;
            gray *pp;
            pm_message( 
              "too many opacs, reducing maxval from %d to %d...",
              wimg->maxval, newmaxval
            );
            for (row = 0; row < wimg->rows; ++row)
              for (col = 0, pp = wimg->pixels[row]; col < wimg->cols; ++col, ++pp)
                PGM_DEPTH((*pp), (*pp), wimg->maxval, newmaxval);
            wimg->maxval = newmaxval;
          }
      }
    
    if (opacs < (*newopacsp)) (*newopacsp) = opacs;
    pm_message("found %d opacs, choosing %d...", opacs, (*newopacsp));

    /* Now apply median-cut to histogram, making the new opacmap. */

    om = median_cut(gh, opacs, bimg->rows*bimg->cols, bimg->maxval, (*newopacsp));

    ppm_freegrayhist(gh);
    return (om);
  }