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;
      }
  }