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