home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / gems / gemsii / quantize.c < prev    next >
C/C++ Source or Header  |  1991-09-22  |  9KB  |  362 lines
  1. /**********************************************************************
  2.         C Implementation of Wu's Color Quantizer
  3.         (see Graphics Gems vol. II, pp. 126-133)
  4.  
  5. Author:    Xiaolin Wu
  6.     Dept. of Computer Science
  7.     Univ. of Western Ontario
  8.     London, Ontario N6A 5B7
  9.     wu@csd.uwo.ca
  10.  
  11. Algorithm: Greedy orthogonal bipartition of RGB space for variance
  12.        minimization aided by inclusion-exclusion tricks.
  13.        For speed no nearest neighbor search is done. Slightly
  14.        better performance can be expected by more sophisticated
  15.        but more expensive versions.
  16.  
  17. Free to distribute, comments and suggestions are appreciated.
  18. **********************************************************************/    
  19.  
  20. #include<stdio.h>
  21.  
  22. #define MAXCOLOR    256
  23. #define    RED    2
  24. #define    GREEN    1
  25. #define BLUE    0
  26.  
  27. struct box {
  28.     int r0;  
  29.     int r1;  
  30.     int g0;  
  31.     int g1;  
  32.     int b0;  
  33.     int b1;
  34.     int vol;
  35. };
  36.  
  37. float        m2[33][33][33];
  38. long int    wt[33][33][33], mr[33][33][33],    mg[33][33][33],    mb[33][33][33];
  39. unsigned char   *Ir, *Ig, *Ib;
  40. int            size; /*image size*/
  41. int        K;    /*color look-up table size*/
  42. unsigned short int *Qadd;
  43.  
  44.  
  45. Hist3d(vwt, vmr, vmg, vmb, m2) /* build 3-D color density */
  46. long int *vwt, *vmr, *vmg, *vmb;
  47. float    *m2;
  48. {
  49. register int ind, r, g, b;
  50. int         inr, ing, inb, table[256];
  51. register long int i;
  52.         
  53.     for(i=0; i<256; ++i) table[i]=i*i;
  54.     Qadd = (unsigned short int *)malloc(sizeof(short int)*size);
  55.     if (Qadd==NULL) {printf("Not enough space\n"); exit(1);}
  56.     for(i=0; i<size; ++i){
  57.         r = Ir[i]; g = Ig[i]; b = Ib[i];
  58.         inr=(r>>3)+1; 
  59.         ing=(g>>3)+1; 
  60.         inb=(b>>3)+1; 
  61.         Qadd[i]=ind=(inr<<10)+(inr<<6)+inr+(ing<<5)+ing+inb;
  62.         /*[inr][ing][inb]*/
  63.         ++vwt[ind];
  64.         vmr[ind] += r;
  65.         vmg[ind] += g;
  66.         vmb[ind] += b;
  67.              m2[ind] += (float)(table[r]+table[g]+table[b]);
  68.     }
  69. }
  70.  
  71.  
  72. M3d(vwt, vmr, vmg, vmb, m2) /* compute cumulative moments. */
  73. long int *vwt, *vmr, *vmg, *vmb;
  74. float    *m2;
  75. {
  76. register unsigned short int ind1, ind2;
  77. register unsigned char i, r, g, b;
  78. long int line, line_r, line_g, line_b,
  79.      area[33], area_r[33], area_g[33], area_b[33];
  80. float    line2, area2[33];
  81.  
  82.     for(r=1; r<=32; ++r){
  83.     for(i=0; i<=32; ++i) 
  84.         area2[i]=area[i]=area_r[i]=area_g[i]=area_b[i]=0;
  85.     for(g=1; g<=32; ++g){
  86.         line2 = line = line_r = line_g = line_b = 0;
  87.         for(b=1; b<=32; ++b){
  88.         ind1 = (r<<10) + (r<<6) + r + (g<<5) + g + b; /* [r][g][b] */
  89.         line += vwt[ind1];
  90.         line_r += vmr[ind1]; 
  91.         line_g += vmg[ind1]; 
  92.         line_b += vmb[ind1];
  93.         line2 += m2[ind1];
  94.         area[b] += line;
  95.         area_r[b] += line_r;
  96.         area_g[b] += line_g;
  97.         area_b[b] += line_b;
  98.         area2[b] += line2;
  99.         ind2 = ind1 - 1089; /* [r-1][g][b] */
  100.         vwt[ind1] = vwt[ind2] + area[b];
  101.         vmr[ind1] = vmr[ind2] + area_r[b];
  102.         vmg[ind1] = vmg[ind2] + area_g[b];
  103.         vmb[ind1] = vmb[ind2] + area_b[b];
  104.         m2[ind1] = m2[ind2] + area2[b];
  105.         }
  106.     }
  107.     }
  108. }
  109.  
  110.  
  111. long int Vol(cube, mmt) 
  112. struct box *cube;
  113. long int mmt[33][33][33];
  114. {
  115.     return( mmt[cube->r1][cube->g1][cube->b1] 
  116.        -mmt[cube->r1][cube->g1][cube->b0]
  117.        -mmt[cube->r1][cube->g0][cube->b1]
  118.        +mmt[cube->r1][cube->g0][cube->b0]
  119.        -mmt[cube->r0][cube->g1][cube->b1]
  120.        +mmt[cube->r0][cube->g1][cube->b0]
  121.        +mmt[cube->r0][cube->g0][cube->b1]
  122.        -mmt[cube->r0][cube->g0][cube->b0] );
  123. }
  124.  
  125.  
  126. long int Bottom(cube, dir, mmt)
  127. struct box *cube;
  128. unsigned char dir;
  129. long int mmt[33][33][33];
  130. {
  131.     switch(dir){
  132.     case RED:
  133.         return( -mmt[cube->r0][cube->g1][cube->b1]
  134.             +mmt[cube->r0][cube->g1][cube->b0]
  135.             +mmt[cube->r0][cube->g0][cube->b1]
  136.             -mmt[cube->r0][cube->g0][cube->b0] );
  137.         break;
  138.     case GREEN:
  139.         return( -mmt[cube->r1][cube->g0][cube->b1]
  140.             +mmt[cube->r1][cube->g0][cube->b0]
  141.             +mmt[cube->r0][cube->g0][cube->b1]
  142.             -mmt[cube->r0][cube->g0][cube->b0] );
  143.         break;
  144.     case BLUE:
  145.         return( -mmt[cube->r1][cube->g1][cube->b0]
  146.             +mmt[cube->r1][cube->g0][cube->b0]
  147.             +mmt[cube->r0][cube->g1][cube->b0]
  148.             -mmt[cube->r0][cube->g0][cube->b0] );
  149.         break;
  150.     }
  151. }
  152.  
  153.  
  154. long int Top(cube, dir, pos, mmt)
  155. struct box *cube;
  156. unsigned char dir;
  157. int   pos;
  158. long int mmt[33][33][33];
  159. {
  160.     switch(dir){
  161.     case RED:
  162.         return( mmt[pos][cube->g1][cube->b1] 
  163.            -mmt[pos][cube->g1][cube->b0]
  164.            -mmt[pos][cube->g0][cube->b1]
  165.            +mmt[pos][cube->g0][cube->b0] );
  166.         break;
  167.     case GREEN:
  168.         return( mmt[cube->r1][pos][cube->b1] 
  169.            -mmt[cube->r1][pos][cube->b0]
  170.            -mmt[cube->r0][pos][cube->b1]
  171.            +mmt[cube->r0][pos][cube->b0] );
  172.         break;
  173.     case BLUE:
  174.         return( mmt[cube->r1][cube->g1][pos]
  175.            -mmt[cube->r1][cube->g0][pos]
  176.            -mmt[cube->r0][cube->g1][pos]
  177.            +mmt[cube->r0][cube->g0][pos] );
  178.         break;
  179.     }
  180. }
  181.  
  182.  
  183. float Var(cube)
  184. struct box *cube;
  185. {
  186. float dr, dg, db, xx;
  187.  
  188.     dr = Vol(cube, mr); 
  189.     dg = Vol(cube, mg); 
  190.     db = Vol(cube, mb);
  191.     xx =  m2[cube->r1][cube->g1][cube->b1] 
  192.      -m2[cube->r1][cube->g1][cube->b0]
  193.      -m2[cube->r1][cube->g0][cube->b1]
  194.      +m2[cube->r1][cube->g0][cube->b0]
  195.      -m2[cube->r0][cube->g1][cube->b1]
  196.      +m2[cube->r0][cube->g1][cube->b0]
  197.      +m2[cube->r0][cube->g0][cube->b1]
  198.      -m2[cube->r0][cube->g0][cube->b0];
  199.  
  200.     return( xx - (dr*dr+dg*dg+db*db)/(float)Vol(cube,wt) );    
  201. }
  202.  
  203.  
  204. float Maximize(cube, dir, first, last, cut,
  205.         whole_r, whole_g, whole_b, whole_w)
  206. struct box *cube;
  207. unsigned char dir;
  208. int first, last, *cut;
  209. long int whole_r, whole_g, whole_b, whole_w;
  210. {
  211. register long int half_r, half_g, half_b, half_w;
  212. long int base_r, base_g, base_b, base_w;
  213. register int i;
  214. register float temp, max;
  215.  
  216.     base_r = Bottom(cube, dir, mr);
  217.     base_g = Bottom(cube, dir, mg);
  218.     base_b = Bottom(cube, dir, mb);
  219.     base_w = Bottom(cube, dir, wt);
  220.     max = 0.0;
  221.     for(i=first; i<last; ++i){
  222.     half_r = base_r + Top(cube, dir, i, mr);
  223.     half_g = base_g + Top(cube, dir, i, mg);
  224.     half_b = base_b + Top(cube, dir, i, mb);
  225.     half_w = base_w + Top(cube, dir, i, wt);
  226.     temp = ((float)half_r*half_r + (float)half_g*half_g + 
  227.                 (float)half_b*half_b)/half_w;
  228.  
  229.     half_r = whole_r - half_r;
  230.     half_g = whole_g - half_g;
  231.     half_b = whole_b - half_b;
  232.     half_w = whole_w - half_w;
  233.     temp = ((float)half_r*half_r + (float)half_g*half_g + 
  234.                 (float)half_b*half_b)/half_w + temp;
  235.  
  236.     if (temp > max) {max=temp; *cut=i;}
  237.     }
  238.     return(max);
  239. }
  240.  
  241. Cut(set1, set2)
  242. struct box *set1, *set2;
  243. {
  244. unsigned char dir;
  245. int cutr, cutg, cutb;
  246. float maxr, maxg, maxb;
  247. long int whole_r, whole_g, whole_b, whole_w;
  248.  
  249.     whole_r = Vol(set1, mr);
  250.     whole_g = Vol(set1, mg);
  251.     whole_b = Vol(set1, mb);
  252.     whole_w = Vol(set1, wt);
  253.  
  254.     maxr = Maximize(set1, RED, set1->r0+1, set1->r1, &cutr,
  255.             whole_r, whole_g, whole_b, whole_w);
  256.     maxg = Maximize(set1, GREEN, set1->g0+1, set1->g1, &cutg,
  257.             whole_r, whole_g, whole_b, whole_w);
  258.     maxb = Maximize(set1, BLUE, set1->b0+1, set1->b1, &cutb,
  259.             whole_r, whole_g, whole_b, whole_w);
  260.  
  261.     if( (maxr>=maxg)&&(maxr>=maxb) ) 
  262.     dir = RED; 
  263.     else
  264.     if( (maxg>=maxr)&&(maxg>=maxb) ) 
  265.     dir = GREEN;
  266.     else
  267.     dir = BLUE; 
  268.  
  269.     set2->r1 = set1->r1;
  270.     set2->g1 = set1->g1;
  271.     set2->b1 = set1->b1;
  272.  
  273.     switch (dir){
  274.     case RED:
  275.         set2->r0 = set1->r1 = cutr;
  276.         set2->g0 = set1->g0;
  277.         set2->b0 = set1->b0;
  278.         break;
  279.     case GREEN:
  280.         set2->g0 = set1->g1 = cutg;
  281.         set2->r0 = set1->r0;
  282.         set2->b0 = set1->b0;
  283.         break;
  284.     case BLUE:
  285.         set2->b0 = set1->b1 = cutb;
  286.         set2->r0 = set1->r0;
  287.         set2->g0 = set1->g0;
  288.         break;
  289.     }
  290.     set1->vol=(set1->r1-set1->r0)*(set1->g1-set1->g0)*(set1->b1-set1->b0);
  291.     set2->vol=(set2->r1-set2->r0)*(set2->g1-set2->g0)*(set2->b1-set2->b0);
  292. }
  293.  
  294.  
  295. Mark(cube, label, tag)
  296. struct box *cube;
  297. int label;
  298. unsigned char *tag;
  299. {
  300. register int r, g, b;
  301.  
  302.     for(r=cube->r0+1; r<=cube->r1; ++r)
  303.        for(g=cube->g0+1; g<=cube->g1; ++g)
  304.       for(b=cube->b0+1; b<=cube->b1; ++b)
  305.         tag[(r<<10) + (r<<6) + r + (g<<5) + g + b] = label;
  306. }
  307.  
  308. main()
  309. {
  310. struct box    cube[MAXCOLOR];
  311. unsigned char    *tag;
  312. unsigned char    lut_r[MAXCOLOR], lut_g[MAXCOLOR], lut_b[MAXCOLOR];
  313. int        next;
  314. register long int    i, weight;
  315. register int    k;
  316. float        vv[MAXCOLOR], temp;
  317.  
  318.     /* input R,G,B components into Ir, Ig, Ib;
  319.        set size to width*height */
  320.  
  321.     printf("no. of colors:\n");
  322.     scanf("%d", &K);
  323.  
  324.     Hist3d(wt, mr, mg, mb, m2); printf("Histogram done\n");
  325.     free(Ig); free(Ib); free(Ir);
  326.  
  327.     M3d(wt, mr, mg, mb, m2);    printf("Moments done\n");
  328.  
  329.     cube[0].r0 = cube[0].g0 = cube[0].b0 = 0;
  330.     cube[0].r1 = cube[0].g1 = cube[0].b1 = 32;
  331.     next = 0;
  332.     for(i=1; i<K; ++i){
  333.         Cut(&cube[next], &cube[i]);
  334.         vv[next] = (cube[next].vol>2) ? Var(&cube[next]) : 0.0;
  335.         vv[i] = (cube[i].vol>2) ? Var(&cube[i]) : 0.0;
  336.         next = 0; temp = vv[0];
  337.         for(k=1; k<=i; ++k)
  338.         if (vv[k] > temp) {
  339.             temp = vv[k]; next = k;
  340.         }
  341.     }
  342.         printf("Partition done\n");
  343.  
  344.     /* the space for array m2 can be freed now */
  345.  
  346.     tag = (unsigned char *)malloc(33*33*33);
  347.     if (tag==NULL) {printf("Not enough space\n"); exit(1);}
  348.     for(k=0; k<K; ++k){
  349.         Mark(&cube[k], k, tag);
  350.         weight = Vol(&cube[k], wt);
  351.         lut_r[k] = Vol(&cube[k], mr) / weight;
  352.         lut_g[k] = Vol(&cube[k], mg) / weight;
  353.         lut_b[k] = Vol(&cube[k], mb) / weight;
  354.     }
  355.  
  356.     for(i=0; i<size; ++i) Qadd[i] = tag[Qadd[i]];
  357.     
  358.     /* output lut_r, lut_g, lut_b as color look-up table contents,
  359.        Qadd as the quantized image (array of table addresses). */
  360. }
  361.  
  362.