Back to index

radiance  4R0+20100331
clrtab.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: clrtab.c,v 2.18 2005/09/19 02:23:58 greg Exp $";
00003 #endif
00004 /*
00005  * Simple median-cut color quantization based on colortab.c
00006  */
00007 
00008 #include "copyright.h"
00009 
00010 #include <string.h>
00011 
00012 #include "standard.h"
00013 #include "color.h"
00014 #include "clrtab.h"
00015 
00016                             /* histogram resolution */
00017 #define NRED         36
00018 #define NGRN         48
00019 #define NBLU         24
00020 #define HMAX         NGRN
00021                             /* minimum box count for adaptive partition */
00022 #define MINSAMP             7
00023                             /* color partition */
00024 #define set_branch(p,c) ((c)<<2|(p))
00025 #define part(cn)     ((cn)>>2)
00026 #define prim(cn)     ((cn)&3)
00027                             /* our color table (global) */
00028 extern BYTE   clrtab[256][3];
00029                             /* histogram of colors / color assignments */
00030 static unsigned      histo[NRED][NGRN][NBLU];
00031 #define cndx(c)             histo[((c)[RED]*NRED)>>8][((c)[GRN]*NGRN)>>8][((c)[BLU]*NBLU)>>8]
00032                             /* initial color cube boundary */
00033 static int    CLRCUBE[3][2] = {{0,NRED},{0,NGRN},{0,NBLU}};
00034                             /* maximum propagated error during dithering */
00035 #define MAXERR              20
00036                             /* define CLOSEST to get closest colors */
00037 #ifndef CLOSEST
00038 #define CLOSEST             1      /* this step takes a little longer */
00039 #endif
00040 
00041 
00042 #ifdef CLOSEST
00043 static void closest(int n);
00044 static void setclosest(BYTE *nl[], int r, int g, int b);
00045 static void addneigh(BYTE *nl[], int i, int j);
00046 static unsigned int dist(BYTE col[3], int r, int g, int b);
00047 #endif
00048 static void cut(int box[3][2], int c0, int c1);
00049 static int split(int box[3][2]);
00050 static void mktabent(int p, int box[3][2]);
00051 
00052 
00053 extern int
00054 new_histo(           /* clear our histogram */
00055        int    n
00056 )
00057 {
00058        memset((void *)histo, '\0', sizeof(histo));
00059        return(0);
00060 }
00061 
00062 
00063 extern void
00064 cnt_pixel(           /* add pixel to our histogram */
00065        register BYTE col[]
00066 )
00067 {
00068        cndx(col)++;
00069 }
00070 
00071 
00072 extern void
00073 cnt_colrs(           /* add a scanline to our histogram */
00074        register COLR *cs,
00075        register int  n
00076 )
00077 {
00078        while (n-- > 0) {
00079               cndx(cs[0])++;
00080               cs++;
00081        }
00082 }
00083 
00084 
00085 extern int
00086 new_clrtab(          /* make new color table using ncolors */
00087        int    ncolors
00088 )
00089 {
00090        if (ncolors < 1)
00091               return(0);
00092        if (ncolors > 256)
00093               ncolors = 256;
00094                             /* partition color space */
00095        cut(CLRCUBE, 0, ncolors);
00096 #ifdef CLOSEST
00097        closest(ncolors);    /* ensure colors picked are closest */
00098 #endif
00099                             /* reset dithering function */
00100        dith_colrs((BYTE *)NULL, (COLR *)NULL, 0);
00101                             /* return new color table size */
00102        return(ncolors);
00103 }
00104 
00105 
00106 extern int
00107 map_pixel(                  /* get pixel for color */
00108        register BYTE col[]
00109 )
00110 {
00111     return(cndx(col));
00112 }
00113 
00114 
00115 extern void
00116 map_colrs(           /* convert a scanline to color index values */
00117        register BYTE *bs,
00118        register COLR *cs,
00119        register int  n
00120 )
00121 {
00122        while (n-- > 0) {
00123               *bs++ = cndx(cs[0]);
00124               cs++;
00125        }
00126 }
00127 
00128 
00129 extern void
00130 dith_colrs(          /* convert scanline to dithered index values */
00131        register BYTE *bs,
00132        register COLR *cs,
00133        int    n
00134 )
00135 {
00136        static short  (*cerr)[3] = NULL;
00137        static int    N = 0;
00138        int    err[3], errp[3];
00139        register int  x, i;
00140 
00141        if (n != N) {        /* get error propogation array */
00142               if (N) {
00143                      free((void *)cerr);
00144                      cerr = NULL;
00145               }
00146               if (n)
00147                      cerr = (short (*)[3])malloc(3*n*sizeof(short));
00148               if (cerr == NULL) {
00149                      N = 0;
00150                      map_colrs(bs, cs, n);
00151                      return;
00152               }
00153               N = n;
00154               memset((void *)cerr, '\0', 3*N*sizeof(short));
00155        }
00156        err[0] = err[1] = err[2] = 0;
00157        for (x = 0; x < n; x++) {
00158               for (i = 0; i < 3; i++) {   /* dither value */
00159                      errp[i] = err[i];
00160                      err[i] += cerr[x][i];
00161 #ifdef MAXERR
00162                      if (err[i] > MAXERR) err[i] = MAXERR;
00163                      else if (err[i] < -MAXERR) err[i] = -MAXERR;
00164 #endif
00165                      err[i] += cs[x][i];
00166                      if (err[i] < 0) err[i] = 0;
00167                      else if (err[i] > 255) err[i] = 255;
00168               }
00169               bs[x] = cndx(err);
00170               for (i = 0; i < 3; i++) {   /* propagate error */
00171                      err[i] -= clrtab[bs[x]][i];
00172                      err[i] /= 3;
00173                      cerr[x][i] = err[i] + errp[i];
00174               }
00175        }
00176 }
00177 
00178 
00179 static void
00180 cut(                 /* partition color space */
00181        register int  box[3][2],
00182        int    c0,
00183        int    c1
00184 )
00185 {
00186        register int  branch;
00187        int    kb[3][2];
00188        
00189        if (c1-c0 <= 1) {           /* assign pixel */
00190               mktabent(c0, box);
00191               return;
00192        }
00193                                    /* split box */
00194        branch = split(box);
00195        memcpy((void *)kb, (void *)box, sizeof(kb));
00196                                           /* do left (lesser) branch */
00197        kb[prim(branch)][1] = part(branch);
00198        cut(kb, c0, (c0+c1)>>1);
00199                                           /* do right branch */
00200        kb[prim(branch)][0] = part(branch);
00201        kb[prim(branch)][1] = box[prim(branch)][1];
00202        cut(kb, (c0+c1)>>1, c1);
00203 }
00204 
00205 
00206 static int
00207 split(                      /* find median cut for box */
00208        register int  box[3][2]
00209 )
00210 {
00211 #define c0    r
00212        register int  r, g, b;
00213        int    pri;
00214        long   t[HMAX], med;
00215                                    /* find dominant axis */
00216        pri = RED;
00217        if (box[GRN][1]-box[GRN][0] > box[pri][1]-box[pri][0])
00218               pri = GRN;
00219        if (box[BLU][1]-box[BLU][0] > box[pri][1]-box[pri][0])
00220               pri = BLU;
00221                                    /* sum histogram over box */
00222        med = 0;
00223        switch (pri) {
00224        case RED:
00225               for (r = box[RED][0]; r < box[RED][1]; r++) {
00226                      t[r] = 0;
00227                      for (g = box[GRN][0]; g < box[GRN][1]; g++)
00228                             for (b = box[BLU][0]; b < box[BLU][1]; b++)
00229                                    t[r] += histo[r][g][b];
00230                      med += t[r];
00231               }
00232               break;
00233        case GRN:
00234               for (g = box[GRN][0]; g < box[GRN][1]; g++) {
00235                      t[g] = 0;
00236                      for (b = box[BLU][0]; b < box[BLU][1]; b++)
00237                             for (r = box[RED][0]; r < box[RED][1]; r++)
00238                                    t[g] += histo[r][g][b];
00239                      med += t[g];
00240               }
00241               break;
00242        case BLU:
00243               for (b = box[BLU][0]; b < box[BLU][1]; b++) {
00244                      t[b] = 0;
00245                      for (r = box[RED][0]; r < box[RED][1]; r++)
00246                             for (g = box[GRN][0]; g < box[GRN][1]; g++)
00247                                    t[b] += histo[r][g][b];
00248                      med += t[b];
00249               }
00250               break;
00251        }
00252        if (med < MINSAMP)          /* if too sparse, split at midpoint */
00253               return(set_branch(pri,(box[pri][0]+box[pri][1])>>1));
00254                                    /* find median position */
00255        med >>= 1;
00256        for (c0 = box[pri][0]; med > 0; c0++)
00257               med -= t[c0];
00258        if (c0 > (box[pri][0]+box[pri][1])>>1)    /* if past the midpoint */
00259               c0--;                       /* part left of median */
00260        return(set_branch(pri,c0));
00261 #undef c0
00262 }
00263 
00264 
00265 static void
00266 mktabent(     /* compute average color for box and assign */
00267        int    p,
00268        register int  box[3][2]
00269 )
00270 {
00271        unsigned long sum[3];
00272        unsigned      r, g;
00273        unsigned long n;
00274        register unsigned    b, c;
00275                                           /* sum pixels in box */
00276        n = 0;
00277        sum[RED] = sum[GRN] = sum[BLU] = 0;
00278        for (r = box[RED][0]; r < box[RED][1]; r++)
00279            for (g = box[GRN][0]; g < box[GRN][1]; g++)
00280               for (b = box[BLU][0]; b < box[BLU][1]; b++) {
00281                   if ( (c = histo[r][g][b]) ) {
00282                      n += c;
00283                      sum[RED] += (long)c*r;
00284                      sum[GRN] += (long)c*g;
00285                      sum[BLU] += (long)c*b;
00286                   }
00287                   histo[r][g][b] = p;            /* assign pixel */
00288               }
00289        if (n >= (1L<<23)/HMAX) {          /* avoid overflow */
00290               sum[RED] /= n;
00291               sum[GRN] /= n;
00292               sum[BLU] /= n;
00293               n = 1;
00294        }
00295        if (n) {                           /* compute average */
00296               clrtab[p][RED] = sum[RED]*256/NRED/n;
00297               clrtab[p][GRN] = sum[GRN]*256/NGRN/n;
00298               clrtab[p][BLU] = sum[BLU]*256/NBLU/n;
00299        } else {                           /* empty box -- use midpoint */
00300               clrtab[p][RED] = (box[RED][0]+box[RED][1])*256/NRED/2;
00301               clrtab[p][GRN] = (box[GRN][0]+box[GRN][1])*256/NGRN/2;
00302               clrtab[p][BLU] = (box[BLU][0]+box[BLU][1])*256/NBLU/2;
00303        }
00304 }
00305 
00306 
00307 #ifdef CLOSEST
00308 #define NBSIZ        32
00309 static void
00310 closest(                    /* make sure we have the closest colors */
00311        int    n
00312 )
00313 {
00314        BYTE   *neigh[256];
00315        register int  r, g, b;
00316 #define i r
00317                                    /* get space for neighbor lists */
00318        for (i = 0; i < n; i++) {
00319               if ((neigh[i] = (BYTE *)malloc(NBSIZ)) == NULL) {
00320                      while (i--)
00321                             free(neigh[i]);
00322                      return;                     /* ENOMEM -- abandon effort */
00323               }
00324               neigh[i][0] = i;            /* identity is terminator */
00325        }
00326                                    /* make neighbor lists */
00327        for (r = 0; r < NRED; r++)
00328            for (g = 0; g < NGRN; g++)
00329               for (b = 0; b < NBLU; b++) {
00330                   if (r < NRED-1 && histo[r][g][b] != histo[r+1][g][b])
00331                      addneigh(neigh, histo[r][g][b], histo[r+1][g][b]);
00332                   if (g < NGRN-1 && histo[r][g][b] != histo[r][g+1][b])
00333                      addneigh(neigh, histo[r][g][b], histo[r][g+1][b]);
00334                   if (b < NBLU-1 && histo[r][g][b] != histo[r][g][b+1])
00335                      addneigh(neigh, histo[r][g][b], histo[r][g][b+1]);
00336               }
00337                                    /* assign closest values */
00338        for (r = 0; r < NRED; r++)
00339            for (g = 0; g < NGRN; g++)
00340               for (b = 0; b < NBLU; b++)
00341                   setclosest(neigh, r, g, b);
00342                                    /* free neighbor lists */
00343        for (i = 0; i < n; i++)
00344               free(neigh[i]);
00345 #undef i
00346 }
00347 
00348 
00349 static void
00350 addneigh(            /* i and j are neighbors; add them to list */
00351        register BYTE *nl[],
00352        register int  i,
00353        int    j
00354 )
00355 {
00356        int    nc;
00357        char   *nnl;
00358        register int  t;
00359        
00360        for (nc = 0; nc < 2; nc++) {              /* do both neighbors */
00361               for (t = 0; nl[i][t] != i; t++)
00362                      if (nl[i][t] == j)
00363                             break;        /* in list already */
00364               if (nl[i][t] == i) {        /* add to list */
00365                      nl[i][t++] = j;
00366                      if (t % NBSIZ == 0) {       /* enlarge list */
00367                             if ((nnl = realloc((void *)nl[i],
00368                                           t+NBSIZ)) == NULL)
00369                                    t--;
00370                             else
00371                                    nl[i] = (BYTE *)nnl;
00372                      }
00373                      nl[i][t] = i;        /* terminator */
00374               }
00375               t = i; i = j; j = t;        /* swap and do it again */
00376        }
00377 }
00378 
00379 
00380 static unsigned int
00381 dist(         /* find distance from clrtab entry to r,g,b */
00382        register BYTE col[3],
00383        int    r,
00384        int    g,
00385        int    b
00386 )
00387 {
00388        register int  tmp;
00389        register unsigned    sum;
00390        
00391        tmp = col[RED]*NRED/256 - r;
00392        sum = tmp*tmp;
00393        tmp = col[GRN]*NGRN/256 - g;
00394        sum += tmp*tmp;
00395        tmp = col[BLU]*NBLU/256 - b;
00396        sum += tmp*tmp;
00397        return(sum);
00398 }
00399 
00400 
00401 static void
00402 setclosest(          /* find index closest to color and assign */
00403        BYTE   *nl[],
00404        int    r,
00405        int    g,
00406        int    b
00407 )
00408 {
00409        int    ident;
00410        unsigned      min;
00411        register unsigned    d;
00412        register BYTE *p;
00413                                    /* get starting value */
00414        min = dist(clrtab[ident=histo[r][g][b]], r, g, b);
00415                                    /* find minimum */
00416        for (p = nl[ident]; *p != ident; p++)
00417               if ((d = dist(clrtab[*p], r, g, b)) < min) {
00418                      min = d;
00419                      histo[r][g][b] = *p;
00420               }
00421 }
00422 #endif