home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / radsrc22 / src / px / clrtab.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-13  |  8.5 KB  |  370 lines
  1. /* Copyright (c) 1992 Regents of the University of California */
  2.  
  3. #ifndef lint
  4. static char SCCSid[] = "@(#)clrtab.c 2.3 10/13/92 LBL";
  5. #endif
  6.  
  7. /*
  8.  * Simple median-cut color quantization based on colortab.c
  9.  */
  10.  
  11. #include "standard.h"
  12.  
  13. #include "color.h"
  14.                 /* histogram resolution */
  15. #define NRED        36
  16. #define NGRN        48
  17. #define NBLU        24
  18. #define HMAX        NGRN
  19.                 /* minimum box count for adaptive partition */
  20. #define MINSAMP        7
  21.                 /* color partition */
  22. #define set_branch(p,c) ((c)<<2|(p))
  23. #define part(cn)    ((cn)>>2)
  24. #define prim(cn)    ((cn)&3)
  25.                 /* our color table (global) */
  26. BYTE    clrtab[256][3];
  27.                 /* histogram of colors / color assignments */
  28. static unsigned    histo[NRED][NGRN][NBLU];
  29. #define cndx(c)        histo[((c)[RED]*NRED)>>8][((c)[GRN]*NGRN)>>8][((c)[BLU]*NBLU)>>8]
  30.                 /* initial color cube boundary */
  31. static int    CLRCUBE[3][2] = {0,NRED,0,NGRN,0,NBLU};
  32.                 /* maximum propagated error during dithering */
  33. #define MAXERR        20
  34.                 /* define CLOSEST to get closest colors */
  35. #define CLOSEST        1
  36.  
  37.  
  38. new_histo()        /* clear our histogram */
  39. {
  40.     bzero((char *)histo, sizeof(histo));
  41. }
  42.  
  43.  
  44. cnt_pixel(col)        /* add pixel to our histogram */
  45. register BYTE    col[];
  46. {
  47.     cndx(col)++;
  48. }
  49.  
  50.  
  51. cnt_colrs(cs, n)        /* add a scanline to our histogram */
  52. register COLR    *cs;
  53. register int    n;
  54. {
  55.     while (n-- > 0) {
  56.         cndx(cs[0])++;
  57.         cs++;
  58.     }
  59. }
  60.  
  61.  
  62. new_clrtab(ncolors)        /* make new color table using ncolors */
  63. int    ncolors;
  64. {
  65.     if (ncolors < 1)
  66.         return(0);
  67.     if (ncolors > 256)
  68.         ncolors = 256;
  69.                 /* partition color space */
  70.     cut(CLRCUBE, 0, ncolors);
  71. #ifdef CLOSEST
  72.     closest(ncolors);    /* ensure colors picked are closest */
  73. #endif
  74.                 /* return new color table size */
  75.     return(ncolors);
  76. }
  77.  
  78.  
  79. int
  80. map_pixel(col)            /* get pixel for color */
  81. register BYTE    col[];
  82. {
  83.     return(cndx(col));
  84. }
  85.  
  86.  
  87. map_colrs(bs, cs, n)        /* convert a scanline to color index values */
  88. register BYTE    *bs;
  89. register COLR    *cs;
  90. register int    n;
  91. {
  92.     while (n-- > 0) {
  93.         *bs++ = cndx(cs[0]);
  94.         cs++;
  95.     }
  96. }
  97.  
  98.  
  99. dith_colrs(bs, cs, n)        /* convert scanline to dithered index values */
  100. register BYTE    *bs;
  101. register COLR    *cs;
  102. int    n;
  103. {
  104.     static short    (*cerr)[3];
  105.     static int    N = 0;
  106.     int    err[3], errp[3];
  107.     register int    x, i;
  108.  
  109.     if (n != N) {        /* get error propogation array */
  110.         if (N)
  111.             cerr = (short (*)[3])realloc((char *)cerr,
  112.                     3*n*sizeof(short));
  113.         else
  114.             cerr = (short (*)[3])malloc(3*n*sizeof(short));
  115.         if (cerr == NULL) {
  116.             N = 0;
  117.             map_colrs(bs, cs, n);
  118.             return;
  119.         }
  120.         N = n;
  121.         bzero((char *)cerr, 3*N*sizeof(short));
  122.     }
  123.     err[0] = err[1] = err[2] = 0;
  124.     for (x = 0; x < n; x++) {
  125.         for (i = 0; i < 3; i++) {    /* dither value */
  126.             errp[i] = err[i];
  127.             err[i] += cerr[x][i];
  128. #ifdef MAXERR
  129.             if (err[i] > MAXERR) err[i] = MAXERR;
  130.             else if (err[i] < -MAXERR) err[i] = -MAXERR;
  131. #endif
  132.             err[i] += cs[x][i];
  133.             if (err[i] < 0) err[i] = 0;
  134.             else if (err[i] > 255) err[i] = 255;
  135.         }
  136.         bs[x] = cndx(err);
  137.         for (i = 0; i < 3; i++) {    /* propagate error */
  138.             err[i] -= clrtab[bs[x]][i];
  139.             err[i] /= 3;
  140.             cerr[x][i] = err[i] + errp[i];
  141.         }
  142.     }
  143. }
  144.  
  145.  
  146. static
  147. cut(box, c0, c1)            /* partition color space */
  148. register int    box[3][2];
  149. int    c0, c1;
  150. {
  151.     register int    branch;
  152.     int    kb[3][2];
  153.     
  154.     if (c1-c0 <= 1) {        /* assign pixel */
  155.         mktabent(c0, box);
  156.         return;
  157.     }
  158.                     /* split box */
  159.     branch = split(box);
  160.     bcopy((char *)box, (char *)kb, sizeof(kb));
  161.                         /* do left (lesser) branch */
  162.     kb[prim(branch)][1] = part(branch);
  163.     cut(kb, c0, (c0+c1)>>1);
  164.                         /* do right branch */
  165.     kb[prim(branch)][0] = part(branch);
  166.     kb[prim(branch)][1] = box[prim(branch)][1];
  167.     cut(kb, (c0+c1)>>1, c1);
  168. }
  169.  
  170.  
  171. static int
  172. split(box)                /* find median cut for box */
  173. register int    box[3][2];
  174. {
  175. #define c0    r
  176.     register int    r, g, b;
  177.     int    pri;
  178.     int    t[HMAX], med;
  179.                     /* find dominant axis */
  180.     pri = RED;
  181.     if (box[GRN][1]-box[GRN][0] > box[pri][1]-box[pri][0])
  182.         pri = GRN;
  183.     if (box[BLU][1]-box[BLU][0] > box[pri][1]-box[pri][0])
  184.         pri = BLU;
  185.                     /* sum histogram over box */
  186.     med = 0;
  187.     switch (pri) {
  188.     case RED:
  189.         for (r = box[RED][0]; r < box[RED][1]; r++) {
  190.             t[r] = 0;
  191.             for (g = box[GRN][0]; g < box[GRN][1]; g++)
  192.                 for (b = box[BLU][0]; b < box[BLU][1]; b++)
  193.                     t[r] += histo[r][g][b];
  194.             med += t[r];
  195.         }
  196.         break;
  197.     case GRN:
  198.         for (g = box[GRN][0]; g < box[GRN][1]; g++) {
  199.             t[g] = 0;
  200.             for (b = box[BLU][0]; b < box[BLU][1]; b++)
  201.                 for (r = box[RED][0]; r < box[RED][1]; r++)
  202.                     t[g] += histo[r][g][b];
  203.             med += t[g];
  204.         }
  205.         break;
  206.     case BLU:
  207.         for (b = box[BLU][0]; b < box[BLU][1]; b++) {
  208.             t[b] = 0;
  209.             for (r = box[RED][0]; r < box[RED][1]; r++)
  210.                 for (g = box[GRN][0]; g < box[GRN][1]; g++)
  211.                     t[b] += histo[r][g][b];
  212.             med += t[b];
  213.         }
  214.         break;
  215.     }
  216.     if (med < MINSAMP)        /* if too sparse, split at midpoint */
  217.         return(set_branch(pri,(box[pri][0]+box[pri][1])>>1));
  218.                     /* find median position */
  219.     med >>= 1;
  220.     for (c0 = box[pri][0]; med > 0; c0++)
  221.         med -= t[c0];
  222.     if (c0 > (box[pri][0]+box[pri][1])>>1)    /* if past the midpoint */
  223.         c0--;                /* part left of median */
  224.     return(set_branch(pri,c0));
  225. #undef c0
  226. }
  227.  
  228.  
  229. static
  230. mktabent(p, box)    /* compute average color for box and assign */
  231. int    p;
  232. register int    box[3][2];
  233. {
  234.     long    sum[3];
  235.     int    r, g, n;
  236.     register int    b, c;
  237.                         /* sum pixels in box */
  238.     n = 0;
  239.     sum[RED] = sum[GRN] = sum[BLU] = 0;
  240.     for (r = box[RED][0]; r < box[RED][1]; r++)
  241.         for (g = box[GRN][0]; g < box[GRN][1]; g++)
  242.             for (b = box[BLU][0]; b < box[BLU][1]; b++) {
  243.                 if (c = histo[r][g][b]) {
  244.             n += c;
  245.             sum[RED] += (long)c*r;
  246.             sum[GRN] += (long)c*g;
  247.             sum[BLU] += (long)c*b;
  248.             }
  249.             histo[r][g][b] = p;        /* assign pixel */
  250.         }
  251.     if (n) {                /* compute average */
  252.         clrtab[p][RED] = sum[RED]*256/NRED/n;
  253.         clrtab[p][GRN] = sum[GRN]*256/NGRN/n;
  254.         clrtab[p][BLU] = sum[BLU]*256/NBLU/n;
  255.     } else {                /* empty box -- use midpoint */
  256.         clrtab[p][RED] = (box[RED][0]+box[RED][1])*256/NRED/2;
  257.         clrtab[p][GRN] = (box[GRN][0]+box[GRN][1])*256/NGRN/2;
  258.         clrtab[p][BLU] = (box[BLU][0]+box[BLU][1])*256/NBLU/2;
  259.     }
  260. }
  261.  
  262.  
  263. #ifdef CLOSEST
  264. #define NBSIZ        32
  265. static
  266. closest(n)            /* make sure we have the closest colors */
  267. int    n;
  268. {
  269.     BYTE    *neigh[256];
  270.     register int    r, g, b;
  271. #define i r
  272.                     /* get space for neighbor lists */
  273.     for (i = 0; i < n; i++) {
  274.         if ((neigh[i] = (BYTE *)malloc(NBSIZ)) == NULL) {
  275.             while (i--)
  276.                 free(neigh[i]);
  277.             return;            /* ENOMEM -- abandon effort */
  278.         }
  279.         neigh[i][0] = i;        /* identity is terminator */
  280.     }
  281.                     /* make neighbor lists */
  282.     for (r = 0; r < NRED; r++)
  283.         for (g = 0; g < NGRN; g++)
  284.         for (b = 0; b < NBLU; b++) {
  285.             if (r < NRED-1 && histo[r][g][b] != histo[r+1][g][b])
  286.             addneigh(neigh, histo[r][g][b], histo[r+1][g][b]);
  287.             if (g < NGRN-1 && histo[r][g][b] != histo[r][g+1][b])
  288.             addneigh(neigh, histo[r][g][b], histo[r][g+1][b]);
  289.             if (b < NBLU-1 && histo[r][g][b] != histo[r][g][b+1])
  290.             addneigh(neigh, histo[r][g][b], histo[r][g][b+1]);
  291.         }
  292.                     /* assign closest values */
  293.     for (r = 0; r < NRED; r++)
  294.         for (g = 0; g < NGRN; g++)
  295.         for (b = 0; b < NBLU; b++)
  296.             setclosest(neigh, r, g, b);
  297.                     /* free neighbor lists */
  298.     for (i = 0; i < n; i++)
  299.         free(neigh[i]);
  300. #undef i
  301. }
  302.  
  303.  
  304. static
  305. addneigh(nl, i, j)        /* i and j are neighbors; add them to list */
  306. register BYTE    *nl[];
  307. register int    i;
  308. int    j;
  309. {
  310.     int    nc;
  311.     char    *nnl;
  312.     register int    t;
  313.     
  314.     for (nc = 0; nc < 2; nc++) {        /* do both neighbors */
  315.         for (t = 0; nl[i][t] != i; t++)
  316.             if (nl[i][t] == j)
  317.                 break;        /* in list already */
  318.         if (nl[i][t] == i) {        /* add to list */
  319.             nl[i][t++] = j;
  320.             if (t % NBSIZ == 0) {    /* enlarge list */
  321.                 if ((nnl = realloc(nl[i], t+NBSIZ)) == NULL)
  322.                     t--;
  323.                 else
  324.                     nl[i] = (BYTE *)nnl;
  325.             }
  326.             nl[i][t] = i;        /* terminator */
  327.         }
  328.         t = i; i = j; j = t;        /* swap and do it again */
  329.     }
  330. }
  331.  
  332.  
  333. static unsigned
  334. dist(col, r, g, b)        /* find distance from clrtab entry to r,g,b */
  335. register BYTE    col[3];
  336. int    r, g, b;
  337. {
  338.     register unsigned    tmp;
  339.     register unsigned    sum;
  340.     
  341.     tmp = col[RED]*NRED/256 - r;
  342.     sum = tmp*tmp;
  343.     tmp = col[GRN]*NGRN/256 - g;
  344.     sum += tmp*tmp;
  345.     tmp = col[BLU]*NBLU/256 - b;
  346.     sum += tmp*tmp;
  347.     return(sum);
  348. }
  349.  
  350.  
  351. static
  352. setclosest(nl, r, g, b)        /* find index closest to color and assign */
  353. BYTE    *nl[];
  354. int    r, g, b;
  355. {
  356.     int    ident;
  357.     unsigned    min;
  358.     register unsigned    d;
  359.     register BYTE    *p;
  360.                     /* get starting value */
  361.     min = dist(clrtab[ident=histo[r][g][b]], r, g, b);
  362.                     /* find minimum */
  363.     for (p = nl[ident]; *p != ident; p++)
  364.         if ((d = dist(clrtab[*p], r, g, b)) < min) {
  365.             min = d;
  366.             histo[r][g][b] = *p;
  367.         }
  368. }
  369. #endif
  370.