home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11.lha / ccs-lib / rlelib / inv_cmap.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-01-15  |  23.1 KB  |  883 lines
  1. /*
  2.  * This software is copyrighted as noted below.  It may be freely copied,
  3.  * modified, and redistributed, provided that the copyright notice is
  4.  * preserved on all copies.
  5.  *
  6.  * There is no warranty or other guarantee of fitness for this software,
  7.  * it is provided solely "as is".  Bug reports or fixes may be sent
  8.  * to the author, who may or may not act on them as he desires.
  9.  *
  10.  * You may not include this software in a program or other software product
  11.  * without supplying the source, or without informing the end-user that the
  12.  * source is available for no extra charge.
  13.  *
  14.  * If you modify this software, you should include a notice giving the
  15.  * name of the person performing the modification, the date of modification,
  16.  * and the reason for such modification.
  17.  */
  18. /*
  19.  * inv_cmap.c - Compute an inverse colormap.
  20.  *
  21.  * Author:    Spencer W. Thomas
  22.  *         EECS Dept.
  23.  *         University of Michigan
  24.  * Date:    Thu Sep 20 1990
  25.  * Copyright (c) 1990, University of Michigan
  26.  *
  27.  * $Id: inv_cmap.c,v 3.0.1.1 90/11/29 15:09:43 spencer Exp $
  28.  */
  29.  
  30. #include <math.h>
  31. #include <stdio.h>
  32.  
  33. #ifndef NO_INV_CMAP_TRACKING
  34.  
  35. #ifdef DEBUG
  36. static int cred, cgreen, cblue;
  37. static int red, green;
  38. static unsigned char *zrgbp;
  39. #endif DEBUG
  40. static int bcenter, gcenter, rcenter;
  41. static long gdist, rdist, cdist;
  42. static long cbinc, cginc, crinc;
  43. static unsigned long *gdp, *rdp, *cdp;
  44. static unsigned char *grgbp, *rrgbp, *crgbp;
  45. static gstride, rstride;
  46. static long x, xsqr, colormax;
  47. static int cindex;
  48. #ifdef INSTRUMENT_IT
  49. static long outercount = 0, innercount = 0;
  50. #endif
  51.  
  52. /*****************************************************************
  53.  * TAG( inv_cmap )
  54.  *
  55.  * Compute an inverse colormap efficiently.
  56.  * Inputs:
  57.  *     colors:        Number of colors in the forward colormap.
  58.  *     colormap:    The forward colormap.
  59.  *     bits:        Number of quantization bits.  The inverse
  60.  *             colormap will have (2^bits)^3 entries.
  61.  *     dist_buf:    An array of (2^bits)^3 long integers to be
  62.  *             used as scratch space.
  63.  * Outputs:
  64.  *     rgbmap:        The output inverse colormap.  The entry
  65.  *             rgbmap[(r<<(2*bits)) + (g<<bits) + b]
  66.  *             is the colormap entry that is closest to the
  67.  *             (quantized) color (r,g,b).
  68.  * Assumptions:
  69.  *     Quantization is performed by right shift (low order bits are
  70.  *     truncated).  Thus, the distance to a quantized color is
  71.  *     actually measured to the color at the center of the cell
  72.  *     (i.e., to r+.5, g+.5, b+.5, if (r,g,b) is a quantized color).
  73.  * Algorithm:
  74.  *     Uses a "distance buffer" algorithm:
  75.  *     The distance from each representative in the forward color map
  76.  *     to each point in the rgb space is computed.  If it is less
  77.  *     than the distance currently stored in dist_buf, then the
  78.  *     corresponding entry in rgbmap is replaced with the current
  79.  *     representative (and the dist_buf entry is replaced with the
  80.  *     new distance).
  81.  *
  82.  *     The distance computation uses an efficient incremental formulation.
  83.  *
  84.  *     Distances are computed "outward" from each color.  If the
  85.  *     colors are evenly distributed in color space, the expected
  86.  *     number of cells visited for color I is N^3/I.
  87.  *     Thus, the complexity of the algorithm is O(log(K) N^3),
  88.  *     where K = colors, and N = 2^bits.
  89.  */
  90.  
  91. /*
  92.  * Here's the idea:  scan from the "center" of each cell "out"
  93.  * until we hit the "edge" of the cell -- that is, the point
  94.  * at which some other color is closer -- and stop.  In 1-D,
  95.  * this is simple:
  96.  *     for i := here to max do
  97.  *         if closer then buffer[i] = this color
  98.  *         else break
  99.  *     repeat above loop with i := here-1 to min by -1
  100.  *
  101.  * In 2-D, it's trickier, because along a "scan-line", the
  102.  * region might start "after" the "center" point.  A picture
  103.  * might clarify:
  104.  *         |    ...
  105.  *               | ...    .
  106.  *              ...        .
  107.  *           ... |      .
  108.  *          .    +         .
  109.  *           .          .
  110.  *            .         .
  111.  *             .........
  112.  *
  113.  * The + marks the "center" of the above region.  On the top 2
  114.  * lines, the region "begins" to the right of the "center".
  115.  *
  116.  * Thus, we need a loop like this:
  117.  *     detect := false
  118.  *     for i := here to max do
  119.  *         if closer then
  120.  *             buffer[..., i] := this color
  121.  *             if !detect then
  122.  *                 here = i
  123.  *                 detect = true
  124.  *         else
  125.  *             if detect then
  126.  *                 break
  127.  *                 
  128.  * Repeat the above loop with i := here-1 to min by -1.  Note that
  129.  * the "detect" value should not be reinitialized.  If it was
  130.  * "true", and center is not inside the cell, then none of the
  131.  * cell lies to the left and this loop should exit
  132.  * immediately.
  133.  *
  134.  * The outer loops are similar, except that the "closer" test
  135.  * is replaced by a call to the "next in" loop; its "detect"
  136.  * value serves as the test.  (No assignment to the buffer is
  137.  * done, either.)
  138.  *
  139.  * Each time an outer loop starts, the "here", "min", and
  140.  * "max" values of the next inner loop should be
  141.  * re-initialized to the center of the cell, 0, and cube size,
  142.  * respectively.  Otherwise, these values will carry over from
  143.  * one "call" to the inner loop to the next.  This tracks the
  144.  * edges of the cell and minimizes the number of
  145.  * "unproductive" comparisons that must be made.
  146.  *
  147.  * Finally, the inner-most loop can have the "if !detect"
  148.  * optimized out of it by splitting it into two loops: one
  149.  * that finds the first color value on the scan line that is
  150.  * in this cell, and a second that fills the cell until
  151.  * another one is closer:
  152.  *      if !detect then        {needed for "down" loop}
  153.  *         for i := here to max do
  154.  *         if closer then
  155.  *             buffer[..., i] := this color
  156.  *             detect := true
  157.  *             break
  158.  *    for i := i+1 to max do
  159.  *        if closer then
  160.  *             buffer[..., i] := this color
  161.  *         else
  162.  *             break
  163.  *
  164.  * In this implementation, each level will require the
  165.  * following variables.  Variables labelled (l) are local to each
  166.  * procedure.  The ? should be replaced with r, g, or b:
  167.  *      cdist:            The distance at the starting point.
  168.  *     ?center:    The value of this component of the color
  169.  *      c?inc:            The initial increment at the ?center position.
  170.  *     ?stride:    The amount to add to the buffer
  171.  *             pointers (dp and rgbp) to get to the
  172.  *             "next row".
  173.  *     min(l):        The "low edge" of the cell, init to 0
  174.  *     max(l):        The "high edge" of the cell, init to
  175.  *             colormax-1
  176.  *     detect(l):        True if this row has changed some
  177.  *                     buffer entries.
  178.  *      i(l):             The index for this row.
  179.  *      ?xx:            The accumulated increment value.
  180.  *      
  181.  *      here(l):        The starting index for this color.  The
  182.  *                      following variables are associated with here,
  183.  *                      in the sense that they must be updated if here
  184.  *                      is changed.
  185.  *      ?dist:            The current distance for this level.  The
  186.  *                      value of dist from the previous level (g or r,
  187.  *                      for level b or g) initializes dist on this
  188.  *                      level.  Thus gdist is associated with here(b)).
  189.  *      ?inc:            The initial increment for the row.
  190.  *      ?dp:            Pointer into the distance buffer.  The value
  191.  *                      from the previous level initializes this level.
  192.  *      ?rgbp:            Pointer into the rgb buffer.  The value
  193.  *                      from the previous level initializes this level.
  194.  * 
  195.  * The blue and green levels modify 'here-associated' variables (dp,
  196.  * rgbp, dist) on the green and red levels, respectively, when here is
  197.  * changed.
  198.  */
  199.  
  200. void
  201. inv_cmap( colors, colormap, bits, dist_buf, rgbmap )
  202. int colors, bits;
  203. unsigned char *colormap[3], *rgbmap;
  204. unsigned long *dist_buf;
  205. {
  206.     int nbits = 8 - bits;
  207.  
  208.     colormax = 1 << bits;
  209.     x = 1 << nbits;
  210.     xsqr = 1 << (2 * nbits);
  211.  
  212.     /* Compute "strides" for accessing the arrays. */
  213.     gstride = colormax;
  214.     rstride = colormax * colormax;
  215.  
  216. #ifdef INSTRUMENT_IT
  217.     outercount = 0;
  218.     innercount = 0;
  219. #endif
  220.     maxfill( dist_buf, colormax );
  221.  
  222.     for ( cindex = 0; cindex < colors; cindex++ )
  223.     {
  224.     /*
  225.      * Distance formula is
  226.      * (red - map[0])^2 + (green - map[1])^2 + (blue - map[2])^2
  227.      *
  228.      * Because of quantization, we will measure from the center of
  229.      * each quantized "cube", so blue distance is
  230.      *     (blue + x/2 - map[2])^2,
  231.      * where x = 2^(8 - bits).
  232.      * The step size is x, so the blue increment is
  233.      *     2*x*blue - 2*x*map[2] + 2*x^2
  234.      *
  235.      * Now, b in the code below is actually blue/x, so our
  236.      * increment will be 2*(b*x^2 + x^2 - x*map[2]).  For
  237.      * efficiency, we will maintain this quantity in a separate variable
  238.      * that will be updated incrementally by adding 2*x^2 each time.
  239.      */
  240.     /* The initial position is the cell containing the colormap
  241.      * entry.  We get this by quantizing the colormap values.
  242.      */
  243.     rcenter = colormap[0][cindex] >> nbits;
  244.     gcenter = colormap[1][cindex] >> nbits;
  245.     bcenter = colormap[2][cindex] >> nbits;
  246. #ifdef DEBUG
  247.     cred = colormap[0][cindex];
  248.     cgreen = colormap[1][cindex];
  249.     cblue = colormap[2][cindex];
  250.     fprintf( stderr, "---Starting %d: %d,%d,%d -> %d,%d,%d\n",
  251.          cindex, cred, cgreen, cblue,
  252.          rcenter, gcenter, bcenter );
  253.     zrgbp = rgbmap;
  254. #endif
  255.  
  256.     rdist = colormap[0][cindex] - (rcenter * x + x/2);
  257.     gdist = colormap[1][cindex] - (gcenter * x + x/2);
  258.     cdist = colormap[2][cindex] - (bcenter * x + x/2);
  259.     cdist = rdist*rdist + gdist*gdist + cdist*cdist;
  260.  
  261.     crinc = 2 * ((rcenter + 1) * xsqr - (colormap[0][cindex] * x));
  262.     cginc = 2 * ((gcenter + 1) * xsqr - (colormap[1][cindex] * x));
  263.     cbinc = 2 * ((bcenter + 1) * xsqr - (colormap[2][cindex] * x));
  264.  
  265.     /* Array starting points. */
  266.     cdp = dist_buf + rcenter * rstride + gcenter * gstride + bcenter;
  267.     crgbp = rgbmap + rcenter * rstride + gcenter * gstride + bcenter;
  268.  
  269.     (void)redloop();
  270.     }
  271. #ifdef INSTRUMENT_IT
  272.     fprintf( stderr, "K = %d, N = %d, outer count = %ld, inner count = %ld\n",
  273.          colors, colormax, outercount, innercount );
  274. #endif
  275. }
  276.  
  277. /* redloop -- loop up and down from red center. */
  278. int
  279. redloop()
  280. {
  281.     int detect;
  282.     int r, i = cindex;
  283.     int first;
  284.     long txsqr = xsqr + xsqr;
  285.     static int here, min, max;
  286.     static long rxx;
  287.  
  288.     detect = 0;
  289.  
  290.     /* Basic loop up. */
  291.     for ( r = rcenter, rdist = cdist, rxx = crinc,
  292.       rdp = cdp, rrgbp = crgbp, first = 1;
  293.       r < colormax;
  294.       r++, rdp += rstride, rrgbp += rstride,
  295.       rdist += rxx, rxx += txsqr, first = 0 )
  296.     {
  297. #ifdef DEBUG
  298.     red = r;
  299. #endif
  300.     if ( greenloop( first ) )
  301.         detect = 1;
  302.     else if ( detect )
  303.         break;
  304.     }
  305.     
  306.     /* Basic loop down. */
  307.     for ( r = rcenter - 1, rxx = crinc - txsqr, rdist = cdist - rxx,
  308.       rdp = cdp - rstride, rrgbp = crgbp - rstride, first = 1;
  309.       r >= 0;
  310.       r--, rdp -= rstride, rrgbp -= rstride,
  311.       rxx -= txsqr, rdist -= rxx, first = 0 )
  312.     {
  313. #ifdef DEBUG
  314.     red = r;
  315. #endif
  316.     if ( greenloop( first ) )
  317.         detect = 1;
  318.     else if ( detect )
  319.         break;
  320.     }
  321.     
  322.     return detect;
  323. }
  324.  
  325. /* greenloop -- loop up and down from green center. */
  326. int
  327. greenloop( restart )
  328. {
  329.     int detect;
  330.     int g, i = cindex;
  331.     int first;
  332.     long txsqr = xsqr + xsqr;
  333.     static int here, min, max;
  334. #ifdef MINMAX_TRACK
  335.     static int prevmax, prevmin;
  336.     int thismax, thismin;
  337. #endif
  338.     static long ginc, gxx, gcdist;    /* "gc" variables maintain correct */
  339.     static unsigned long *gcdp;        /*  values for bcenter position, */
  340.     static unsigned char *gcrgbp;    /*  despite modifications by blueloop */
  341.                     /*  to gdist, gdp, grgbp. */
  342.  
  343.     if ( restart )
  344.     {
  345.     here = gcenter;
  346.     min = 0;
  347.     max = colormax - 1;
  348.     ginc = cginc;
  349. #ifdef MINMAX_TRACK
  350.     prevmax = 0;
  351.     prevmin = colormax;
  352. #endif
  353.     }
  354.  
  355. #ifdef MINMAX_TRACK
  356.     thismin = min;
  357.     thismax = max;
  358. #endif
  359.     detect = 0;
  360.  
  361.     /* Basic loop up. */
  362.     for ( g = here, gcdist = gdist = rdist, gxx = ginc,
  363.       gcdp = gdp = rdp, gcrgbp = grgbp = rrgbp, first = 1;
  364.       g <= max;
  365.       g++, gdp += gstride, gcdp += gstride, grgbp += gstride, gcrgbp += gstride,
  366.       gdist += gxx, gcdist += gxx, gxx += txsqr, first = 0 )
  367.     {
  368. #ifdef DEBUG
  369.     green = g;
  370. #endif
  371.     if ( blueloop( first ) )
  372.     {
  373.         if ( !detect )
  374.         {
  375.         /* Remember here and associated data! */
  376.         if ( g > here )
  377.         {
  378.             here = g;
  379.             rdp = gcdp;
  380.             rrgbp = gcrgbp;
  381.             rdist = gcdist;
  382.             ginc = gxx;
  383. #ifdef MINMAX_TRACK
  384.             thismin = here;
  385. #endif
  386. #ifdef DEBUG
  387.             fprintf( stderr, "===Adjusting green here up at %d,%d\n",
  388.                  red, here );
  389. #endif
  390.         }
  391.         detect = 1;
  392.         }
  393.     }
  394.     else if ( detect )
  395.     {
  396. #ifdef MINMAX_TRACK
  397.         thismax = g - 1;
  398. #endif
  399.         break;
  400.     }
  401.     }
  402.     
  403.     /* Basic loop down. */
  404.     for ( g = here - 1, gxx = ginc - txsqr, gcdist = gdist = rdist - gxx,
  405.       gcdp = gdp = rdp - gstride, gcrgbp = grgbp = rrgbp - gstride,
  406.       first = 1;
  407.       g >= min;
  408.       g--, gdp -= gstride, gcdp -= gstride, grgbp -= gstride, gcrgbp -= gstride,
  409.       gxx -= txsqr, gdist -= gxx, gcdist -= gxx, first = 0 )
  410.     {
  411. #ifdef DEBUG
  412.     green = g;
  413. #endif
  414.     if ( blueloop( first ) )
  415.     {
  416.         if ( !detect )
  417.         {
  418.         /* Remember here! */
  419.         here = g;
  420.         rdp = gcdp;
  421.         rrgbp = gcrgbp;
  422.         rdist = gcdist;
  423.         ginc = gxx;
  424. #ifdef MINMAX_TRACK
  425.         thismax = here;
  426. #endif
  427. #ifdef DEBUG
  428.         fprintf( stderr, "===Adjusting green here down at %d,%d\n",
  429.              red, here );
  430. #endif
  431.         detect = 1;
  432.         }
  433.     }
  434.     else if ( detect )
  435.     {
  436. #ifdef MINMAX_TRACK
  437.         thismin = g + 1;
  438. #endif
  439.         break;
  440.     }
  441.     }
  442.     
  443. #ifdef MINMAX_TRACK
  444.     /* If we saw something, update the edge trackers.  For now, only
  445.      * tracks edges that are "shrinking" (min increasing, max
  446.      * decreasing.
  447.      */
  448.     if ( detect )
  449.     {
  450.     if ( thismax < prevmax )
  451.         max = thismax;
  452.  
  453.     prevmax = thismax;
  454.  
  455.     if ( thismin > prevmin )
  456.         min = thismin;
  457.  
  458.     prevmin = thismin;
  459.     }
  460. #endif
  461.  
  462.     return detect;
  463. }
  464.  
  465. /* blueloop -- loop up and down from blue center. */
  466. int
  467. blueloop( restart )
  468. {
  469.     int detect;
  470.     register unsigned long *dp;
  471.     register unsigned char *rgbp;
  472.     register long bdist, bxx;
  473.     register int b, i = cindex;
  474.     register long txsqr = xsqr + xsqr;
  475.     register int lim;
  476.     static int here, min, max;
  477. #ifdef MINMAX_TRACK
  478.     static int prevmin, prevmax;
  479.     int thismin, thismax;
  480. #ifdef DMIN_DMAX_TRACK
  481.     static int dmin, dmax;
  482. #endif /* DMIN_DMAX_TRACK */
  483. #endif /* MINMAX_TRACK */
  484.     static long binc;
  485. #ifdef DEBUG
  486.     long dist, tdist;
  487. #endif
  488.  
  489.     if ( restart )
  490.     {
  491.     here = bcenter;
  492.     min = 0;
  493.     max = colormax - 1;
  494.     binc = cbinc;
  495. #ifdef MINMAX_TRACK
  496.     prevmin = colormax;
  497.     prevmax = 0;
  498. #ifdef DMIN_DMAX_TRACK
  499.     dmin = 0;
  500.     dmax = 0;
  501. #endif /* DMIN_DMAX_TRACK */
  502. #endif /* MINMAX_TRACK */
  503.     }
  504.  
  505.     detect = 0;
  506. #ifdef MINMAX_TRACK
  507.     thismin = min;
  508.     thismax = max;
  509. #endif
  510. #ifdef DEBUG
  511.     tdist = cred - red * x - x/2;
  512.     dist = tdist*tdist;
  513.     tdist = cgreen - green * x - x/2;
  514.     dist += tdist*tdist;
  515.     tdist = cblue - here * x - x/2;
  516.     dist += tdist*tdist;
  517.     if ( gdist != dist )
  518.     fprintf( stderr, "*** At %d,%d,%d; dist = %ld != gdist = %ld\n",
  519.          red, green, here, dist, gdist );
  520.     if ( grgbp != zrgbp + red*rstride + green*gstride + here )
  521.     fprintf( stderr, "*** At %d,%d,%d: buffer pointer is at %d,%d,%d\n",
  522.          red, green, here,
  523.          (grgbp - zrgbp) / rstride,
  524.          ((grgbp - zrgbp) % rstride) / gstride,
  525.          (grgbp - zrgbp) % gstride );
  526. #endif DEBUG
  527.  
  528.     /* Basic loop up. */
  529.     /* First loop just finds first applicable cell. */
  530.     for ( b = here, bdist = gdist, bxx = binc, dp = gdp, rgbp = grgbp, lim = max;
  531.       b <= lim;
  532.       b++, dp++, rgbp++,
  533.       bdist += bxx, bxx += txsqr )
  534.     {
  535. #ifdef INSTRUMENT_IT
  536.     outercount++;
  537. #endif
  538.     if ( *dp > bdist )
  539.     {
  540.         /* Remember new 'here' and associated data! */
  541.         if ( b > here )
  542.         {
  543.         here = b;
  544.         gdp = dp;
  545.         grgbp = rgbp;
  546.         gdist = bdist;
  547.         binc = bxx;
  548. #ifdef MINMAX_TRACK
  549.         thismin = here;
  550. #endif
  551. #ifdef DEBUG
  552.         fprintf( stderr, "===Adjusting blue here up at %d,%d,%d\n",
  553.              red, green, here );
  554.  
  555.         tdist = cred - red * x - x/2;
  556.         dist = tdist*tdist;
  557.         tdist = cgreen - green * x - x/2;
  558.         dist += tdist*tdist;
  559.         tdist = cblue - here * x - x/2;
  560.         dist += tdist*tdist;
  561.         if ( gdist != dist )
  562.             fprintf( stderr,
  563.          "*** Adjusting here up at %d,%d,%d; dist = %ld != gdist = %ld\n",
  564.                  red, green, here, dist, gdist );
  565. #endif DEBUG
  566.         }
  567.         detect = 1;
  568. #ifdef INSTRUMENT_IT
  569.         outercount--;
  570. #endif
  571.         break;
  572.     }
  573.     }
  574.     /* Second loop fills in a run of closer cells. */
  575.     for ( ;
  576.       b <= lim;
  577.       b++, dp++, rgbp++,
  578.       bdist += bxx, bxx += txsqr )
  579.     {
  580. #ifdef INSTRUMENT_IT
  581.     outercount++;
  582. #endif
  583.     if ( *dp > bdist )
  584.     {
  585. #ifdef INSTRUMENT_IT
  586.         innercount++;
  587. #endif
  588.         *dp = bdist;
  589.         *rgbp = i;
  590.     }
  591.     else
  592.     {
  593. #ifdef MINMAX_TRACK
  594.         thismax = b - 1;
  595. #endif
  596.         break;
  597.     }
  598.     }
  599. #ifdef DEBUG
  600.     tdist = cred - red * x - x/2;
  601.     dist = tdist*tdist;
  602.     tdist = cgreen - green * x - x/2;
  603.     dist += tdist*tdist;
  604.     tdist = cblue - b * x - x/2;
  605.     dist += tdist*tdist;
  606.     if ( bdist != dist )
  607.     fprintf( stderr,
  608.          "*** After up loop at %d,%d,%d; dist = %ld != bdist = %ld\n",
  609.          red, green, b, dist, bdist );
  610. #endif DEBUG
  611.     
  612.     /* Basic loop down. */
  613.     /* Do initializations here, since the 'find' loop might not get
  614.      * executed. 
  615.      */
  616.     lim = min;
  617.     b = here - 1;
  618.     bxx = binc - txsqr;
  619.     bdist = gdist - bxx;
  620.     dp = gdp - 1;
  621.     rgbp = grgbp - 1;
  622.     /* The 'find' loop is executed only if we didn't already find
  623.      * something.
  624.      */
  625.     if ( !detect )
  626.     for ( ;
  627.           b >= lim;
  628.           b--, dp--, rgbp--,
  629.           bxx -= txsqr, bdist -= bxx )
  630.     {
  631. #ifdef INSTRUMENT_IT
  632.         outercount++;
  633. #endif
  634.         if ( *dp > bdist )
  635.         {
  636.         /* Remember here! */
  637.         /* No test for b against here necessary because b <
  638.          * here by definition.
  639.          */
  640.         here = b;
  641.         gdp = dp;
  642.         grgbp = rgbp;
  643.         gdist = bdist;
  644.         binc = bxx;
  645. #ifdef MINMAX_TRACK
  646.         thismax = here;
  647. #endif
  648.         detect = 1;
  649. #ifdef DEBUG
  650.         fprintf( stderr, "===Adjusting blue here down at %d,%d,%d\n",
  651.              red, green, here );
  652.         tdist = cred - red * x - x/2;
  653.         dist = tdist*tdist;
  654.         tdist = cgreen - green * x - x/2;
  655.         dist += tdist*tdist;
  656.         tdist = cblue - here * x - x/2;
  657.         dist += tdist*tdist;
  658.         if ( gdist != dist )
  659.             fprintf( stderr,
  660.          "*** Adjusting here down at %d,%d,%d; dist = %ld != gdist = %ld\n",
  661.                  red, green, here, dist, gdist );
  662. #endif DEBUG
  663. #ifdef INSTRUMENT_IT
  664.         outercount--;
  665. #endif
  666.         break;
  667.         }
  668.     }
  669.     /* The 'update' loop. */
  670.     for ( ;
  671.       b >= lim;
  672.       b--, dp--, rgbp--,
  673.       bxx -= txsqr, bdist -= bxx )
  674.     {
  675. #ifdef INSTRUMENT_IT
  676.     outercount++;
  677. #endif
  678.     if ( *dp > bdist )
  679.     {
  680. #ifdef INSTRUMENT_IT
  681.         innercount++;
  682. #endif
  683.         *dp = bdist;
  684.         *rgbp = i;
  685.     }
  686.     else
  687.     {
  688. #ifdef MINMAX_TRACK
  689.         thismin = b + 1;
  690. #endif
  691.         break;
  692.     }
  693.     }
  694.  
  695. #ifdef DEBUG
  696.     tdist = cred - red * x - x/2;
  697.     dist = tdist*tdist;
  698.     tdist = cgreen - green * x - x/2;
  699.     dist += tdist*tdist;
  700.     tdist = cblue - b * x - x/2;
  701.     dist += tdist*tdist;
  702.     if ( bdist != dist )
  703.     fprintf( stderr,
  704.          "*** After down loop at %d,%d,%d; dist = %ld != bdist = %ld\n",
  705.          red, green, b, dist, bdist );
  706. #endif DEBUG
  707.  
  708.     /* If we saw something, update the edge trackers. */
  709. #ifdef MINMAX_TRACK
  710.     if ( detect )
  711.     {
  712. #ifdef DMIN_DMAX_TRACK
  713.     /* Predictively track.  Not clear this is a win. */
  714.     /* If there was a previous line, update the dmin/dmax values. */
  715.     if ( prevmax >= prevmin )
  716.     {
  717.         if ( thismin > 0 )
  718.         dmin = thismin - prevmin - 1;
  719.         else
  720.         dmin = 0;
  721.         if ( thismax < colormax - 1 )
  722.         dmax = thismax - prevmax + 1;
  723.         else
  724.         dmax = 0;
  725.         /* Update the min and max values by the differences. */
  726.         max = thismax + dmax;
  727.         if ( max >= colormax )
  728.         max = colormax - 1;
  729.         min = thismin + dmin;
  730.         if ( min < 0 )
  731.         min = 0;
  732.     }
  733. #else /* !DMIN_DMAX_TRACK */
  734.     /* Only tracks edges that are "shrinking" (min increasing, max
  735.      * decreasing.
  736.      */
  737.     if ( thismax < prevmax )
  738.         max = thismax;
  739.  
  740.     if ( thismin > prevmin )
  741.         min = thismin;
  742. #endif /* DMIN_DMAX_TRACK */
  743.     
  744.     /* Remember the min and max values. */
  745.     prevmax = thismax;
  746.     prevmin = thismin;
  747.     }
  748. #endif /* MINMAX_TRACK */
  749.  
  750.     return detect;
  751. }
  752.  
  753. maxfill( buffer, side )
  754. unsigned long *buffer;
  755. long side;
  756. {
  757.     register unsigned long maxv = ~0L;
  758.     register long i;
  759.     register unsigned long *bp;
  760.  
  761.     for ( i = colormax * colormax * colormax, bp = buffer;
  762.       i > 0;
  763.       i--, bp++ )
  764.     *bp = maxv;
  765. }
  766.  
  767. #else /* !NO_INV_CMAP_TRACKING */
  768. /*****************************************************************
  769.  * TAG( inv_cmap )
  770.  *
  771.  * Compute an inverse colormap efficiently.
  772.  * Inputs:
  773.  *     colors:        Number of colors in the forward colormap.
  774.  *     colormap:    The forward colormap.
  775.  *     bits:        Number of quantization bits.  The inverse
  776.  *             colormap will have (2^bits)^3 entries.
  777.  *     dist_buf:    An array of (2^bits)^3 long integers to be
  778.  *             used as scratch space.
  779.  * Outputs:
  780.  *     rgbmap:        The output inverse colormap.  The entry
  781.  *             rgbmap[(r<<(2*bits)) + (g<<bits) + b]
  782.  *             is the colormap entry that is closest to the
  783.  *             (quantized) color (r,g,b).
  784.  * Assumptions:
  785.  *     Quantization is performed by right shift (low order bits are
  786.  *     truncated).  Thus, the distance to a quantized color is
  787.  *     actually measured to the color at the center of the cell
  788.  *     (i.e., to r+.5, g+.5, b+.5, if (r,g,b) is a quantized color).
  789.  * Algorithm:
  790.  *     Uses a "distance buffer" algorithm:
  791.  *     The distance from each representative in the forward color map
  792.  *     to each point in the rgb space is computed.  If it is less
  793.  *     than the distance currently stored in dist_buf, then the
  794.  *     corresponding entry in rgbmap is replaced with the current
  795.  *     representative (and the dist_buf entry is replaced with the
  796.  *     new distance).
  797.  *
  798.  *     The distance computation uses an efficient incremental formulation.
  799.  *
  800.  *     Right now, distances are computed for all entries in the rgb
  801.  *     space.  Thus, the complexity of the algorithm is O(K N^3),
  802.  *     where K = colors, and N = 2^bits.
  803.  */
  804. void
  805. inv_cmap( colors, colormap, bits, dist_buf, rgbmap )
  806. int colors, bits;
  807. unsigned char *colormap[3], *rgbmap;
  808. unsigned long *dist_buf;
  809. {
  810.     register unsigned long *dp;
  811.     register unsigned char *rgbp;
  812.     register long bdist, bxx;
  813.     register int b, i;
  814.     int nbits = 8 - bits;
  815.     register int colormax = 1 << bits;
  816.     register long xsqr = 1 << (2 * nbits);
  817.     int x = 1 << nbits;
  818.     int rinc, ginc, binc, r, g;
  819.     long rdist, gdist, rxx, gxx;
  820. #ifdef INSTRUMENT_IT
  821.     long outercount = 0, innercount = 0;
  822. #endif
  823.  
  824.     for ( i = 0; i < colors; i++ )
  825.     {
  826.     /*
  827.      * Distance formula is
  828.      * (red - map[0])^2 + (green - map[1])^2 + (blue - map[2])^2
  829.      *
  830.      * Because of quantization, we will measure from the center of
  831.      * each quantized "cube", so blue distance is
  832.      *     (blue + x/2 - map[2])^2,
  833.      * where x = 2^(8 - bits).
  834.      * The step size is x, so the blue increment is
  835.      *     2*x*blue - 2*x*map[2] + 2*x^2
  836.      *
  837.      * Now, b in the code below is actually blue/x, so our
  838.      * increment will be 2*x*x*b + (2*x^2 - 2*x*map[2]).  For
  839.      * efficiency, we will maintain this quantity in a separate variable
  840.      * that will be updated incrementally by adding 2*x^2 each time.
  841.      */
  842.     rdist = colormap[0][i] - x/2;
  843.     gdist = colormap[1][i] - x/2;
  844.     bdist = colormap[2][i] - x/2;
  845.     rdist = rdist*rdist + gdist*gdist + bdist*bdist;
  846.  
  847.     rinc = 2 * (xsqr - (colormap[0][i] << nbits));
  848.     ginc = 2 * (xsqr - (colormap[1][i] << nbits));
  849.     binc = 2 * (xsqr - (colormap[2][i] << nbits));
  850.     dp = dist_buf;
  851.     rgbp = rgbmap;
  852.     for ( r = 0, rxx = rinc;
  853.           r < colormax;
  854.           rdist += rxx, r++, rxx += xsqr + xsqr )
  855.         for ( g = 0, gdist = rdist, gxx = ginc;
  856.           g < colormax;
  857.           gdist += gxx, g++, gxx += xsqr + xsqr )
  858.         for ( b = 0, bdist = gdist, bxx = binc;
  859.               b < colormax;
  860.               bdist += bxx, b++, dp++, rgbp++,
  861.               bxx += xsqr + xsqr )
  862.         {
  863. #ifdef INSTRUMENT_IT
  864.             outercount++;
  865. #endif
  866.             if ( i == 0 || *dp > bdist )
  867.             {
  868. #ifdef INSTRUMENT_IT
  869.             innercount++;
  870. #endif
  871.             *dp = bdist;
  872.             *rgbp = i;
  873.             }
  874.         }
  875.     }
  876. #ifdef INSTRUMENT_IT
  877.     fprintf( stderr, "K = %d, N = %d, outer count = %ld, inner count = %ld\n",
  878.          colors, colormax, outercount, innercount );
  879. #endif
  880. }
  881.  
  882. #endif /* NO_INV_CMAP_TRACKING */
  883.