Back to index

radiance  4R0+20100331
colortab.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: colortab.c,v 2.9 2004/03/30 16:13:01 schorsch Exp $";
00003 #endif
00004 /*
00005  * colortab.c - allocate and control dynamic color table.
00006  *
00007  *     We start off with a uniform partition of color space.
00008  *     As pixels are sent to the frame buffer, a histogram is built.
00009  *     When a new color table is requested, the histogram is used
00010  *     to make a pseudo-optimal partition, after which the
00011  *     histogram is cleared.  This algorithm
00012  *     performs only as well as the next drawing's color
00013  *     distribution is correlated to the last.
00014  *
00015  *  External symbols declared in drvier.h
00016  */
00017 
00018 #include "copyright.h"
00019 
00020 #include <string.h>
00021 
00022 #include "standard.h"
00023 #include "color.h"
00024 #include "driver.h"
00025 
00026                             /* histogram resolution */
00027 #define NRED         24
00028 #define NGRN         32
00029 #define NBLU         16
00030 #define HMAX         NGRN
00031                             /* minimum box count for adaptive partition */
00032 #define MINSAMP             7
00033                             /* maximum distance^2 before color reassign */
00034 #define MAXDST2             12
00035                             /* map a color */
00036 #define map_col(c,p) clrmap[p][ colval(c,p)<1. ? \
00037                             (int)(colval(c,p)*256.) : 255 ]
00038                             /* color partition tree */
00039 #define CNODE        short
00040 #define set_branch(p,c) ((c)<<2|(p))
00041 #define set_pval(pv) ((pv)<<2|3)
00042 #define is_branch(cn)       (((cn)&3)!=3)
00043 #define is_pval(cn)  (((cn)&3)==3)
00044 #define part(cn)     ((cn)>>2)
00045 #define prim(cn)     ((cn)&3)
00046 #define pval(cn)     ((cn)>>2)
00047                             /* our color table */
00048 static struct tabent {
00049        long   sum[3];              /* sum of colors using this entry */
00050        int    n;            /* number of colors */
00051        BYTE   ent[3];              /* current table value */
00052 }      *clrtab = NULL;
00053                             /* color cube partition */
00054 static CNODE  *ctree = NULL;
00055                             /* our color correction map */
00056 static BYTE   clrmap[3][256];
00057                             /* histogram of colors used */
00058 static unsigned short       histo[NRED][NGRN][NBLU];
00059                             /* initial color cube boundary */
00060 static int    CLRCUBE[3][2] = {{0,NRED},{0,NGRN},{0,NBLU}};
00061 
00062 static int    split();
00063 static void   cut();
00064 
00065 
00066 int
00067 new_ctab(ncolors)           /* start new color table with max ncolors */
00068 int    ncolors;
00069 {
00070        int    treesize;
00071 
00072        if (ncolors < 1)
00073               return(0);
00074                             /* free old tables */
00075        if (clrtab != NULL)
00076               free((void *)clrtab);
00077        if (ctree != NULL)
00078               free((void *)ctree);
00079                             /* get new tables */
00080        for (treesize = 1; treesize < ncolors; treesize <<= 1)
00081               ;
00082        treesize <<= 1;
00083        clrtab = (struct tabent *)calloc(ncolors, sizeof(struct tabent));
00084        ctree = (CNODE *)malloc(treesize*sizeof(CNODE));
00085        if (clrtab == NULL || ctree == NULL)
00086               return(0);
00087                             /* partition color space */
00088        cut(ctree, 0, CLRCUBE, 0, ncolors);
00089                             /* clear histogram */
00090        memset((void *)histo, '\0', sizeof(histo));
00091                             /* return number of colors used */
00092        return(ncolors);
00093 }
00094 
00095 
00096 extern int
00097 get_pixel(    /* get pixel for color */
00098        COLOR  col,
00099        dr_newcolrf_t *newcolr
00100 )
00101 {
00102        int    r, g, b;
00103        int    cv[3];
00104        register CNODE       *tp;
00105        register int  h;
00106                                           /* map color */
00107        r = map_col(col,RED);
00108        g = map_col(col,GRN);
00109        b = map_col(col,BLU);
00110                                           /* reduce resolution */
00111        cv[RED] = (r*NRED)>>8;
00112        cv[GRN] = (g*NGRN)>>8;
00113        cv[BLU] = (b*NBLU)>>8;
00114                                           /* add to histogram */
00115        histo[cv[RED]][cv[GRN]][cv[BLU]]++;
00116                                           /* find pixel in tree */
00117        for (tp = ctree, h = 0; is_branch(*tp); h++)
00118               if (cv[prim(*tp)] < part(*tp))
00119                      tp += 1<<h;          /* left branch */
00120               else
00121                      tp += 1<<(h+1);             /* right branch */
00122        h = pval(*tp);
00123                                           /* add to color table */
00124        clrtab[h].sum[RED] += r;
00125        clrtab[h].sum[GRN] += g;
00126        clrtab[h].sum[BLU] += b;
00127        clrtab[h].n++;
00128                                    /* recompute average */
00129        r = clrtab[h].sum[RED] / clrtab[h].n;
00130        g = clrtab[h].sum[GRN] / clrtab[h].n;
00131        b = clrtab[h].sum[BLU] / clrtab[h].n;
00132                                    /* check for movement */
00133        if (clrtab[h].n == 1 ||
00134                      (r-clrtab[h].ent[RED])*(r-clrtab[h].ent[RED]) +
00135                      (g-clrtab[h].ent[GRN])*(g-clrtab[h].ent[GRN]) +
00136                      (b-clrtab[h].ent[BLU])*(b-clrtab[h].ent[BLU]) > MAXDST2) {
00137               clrtab[h].ent[RED] = r;
00138               clrtab[h].ent[GRN] = g; /* reassign pixel */
00139               clrtab[h].ent[BLU] = b;
00140 #ifdef DEBUG
00141               sprintf(errmsg, "pixel %d = (%d,%d,%d) (%d refs)\n",
00142                             h, r, g, b, clrtab[h].n);
00143               eputs(errmsg);
00144 #endif
00145               (*newcolr)(h, r, g, b);
00146        }
00147        return(h);                         /* return pixel value */
00148 }
00149 
00150 
00151 void
00152 make_gmap(gam)                     /* make gamma correction map */
00153 double gam;
00154 {
00155        register int  i;
00156        
00157        for (i = 0; i < 256; i++)
00158               clrmap[RED][i] =
00159               clrmap[GRN][i] =
00160               clrmap[BLU][i] = 256.0 * pow((i+0.5)/256.0, 1.0/gam);
00161 }
00162 
00163 
00164 void
00165 set_cmap(rmap, gmap, bmap)  /* set custom color correction map */
00166 BYTE   *rmap, *gmap, *bmap;
00167 {
00168        memcpy((void *)clrmap[RED], (void *)rmap, 256);
00169        memcpy((void *)clrmap[GRN], (void *)gmap, 256);
00170        memcpy((void *)clrmap[BLU], (void *)bmap, 256);
00171 }
00172 
00173 
00174 void
00175 map_color(rgb, col)         /* map a color to a byte triplet */
00176 BYTE   rgb[3];
00177 COLOR  col;
00178 {
00179        rgb[RED] = map_col(col,RED);
00180        rgb[GRN] = map_col(col,GRN);
00181        rgb[BLU] = map_col(col,BLU);
00182 }
00183 
00184 
00185 static void
00186 cut(tree, level, box, c0, c1)             /* partition color space */
00187 register CNODE       *tree;
00188 int    level;
00189 register int  box[3][2];
00190 int    c0, c1;
00191 {
00192        int    kb[3][2];
00193        
00194        if (c1-c0 <= 1) {           /* assign pixel */
00195               *tree = set_pval(c0);
00196               return;
00197        }
00198                                    /* split box */
00199        *tree = split(box);
00200        memcpy((void *)kb, (void *)box, sizeof(kb));
00201                                           /* do left (lesser) branch */
00202        kb[prim(*tree)][1] = part(*tree);
00203        cut(tree+(1<<level), level+1, kb, c0, (c0+c1)>>1);
00204                                           /* do right branch */
00205        kb[prim(*tree)][0] = part(*tree);
00206        kb[prim(*tree)][1] = box[prim(*tree)][1];
00207        cut(tree+(1<<(level+1)), level+1, kb, (c0+c1)>>1, c1);
00208 }
00209 
00210 
00211 static int
00212 split(box)                         /* find median cut for box */
00213 register int  box[3][2];
00214 {
00215 #define c0    r
00216        register int  r, g, b;
00217        int    pri;
00218        long   t[HMAX], med;
00219                                    /* find dominant axis */
00220        pri = RED;
00221        if (box[GRN][1]-box[GRN][0] > box[pri][1]-box[pri][0])
00222               pri = GRN;
00223        if (box[BLU][1]-box[BLU][0] > box[pri][1]-box[pri][0])
00224               pri = BLU;
00225                                    /* sum histogram over box */
00226        med = 0;
00227        switch (pri) {
00228        case RED:
00229               for (r = box[RED][0]; r < box[RED][1]; r++) {
00230                      t[r] = 0;
00231                      for (g = box[GRN][0]; g < box[GRN][1]; g++)
00232                             for (b = box[BLU][0]; b < box[BLU][1]; b++)
00233                                    t[r] += histo[r][g][b];
00234                      med += t[r];
00235               }
00236               break;
00237        case GRN:
00238               for (g = box[GRN][0]; g < box[GRN][1]; g++) {
00239                      t[g] = 0;
00240                      for (b = box[BLU][0]; b < box[BLU][1]; b++)
00241                             for (r = box[RED][0]; r < box[RED][1]; r++)
00242                                    t[g] += histo[r][g][b];
00243                      med += t[g];
00244               }
00245               break;
00246        case BLU:
00247               for (b = box[BLU][0]; b < box[BLU][1]; b++) {
00248                      t[b] = 0;
00249                      for (r = box[RED][0]; r < box[RED][1]; r++)
00250                             for (g = box[GRN][0]; g < box[GRN][1]; g++)
00251                                    t[b] += histo[r][g][b];
00252                      med += t[b];
00253               }
00254               break;
00255        }
00256        if (med < MINSAMP)          /* if too sparse, split at midpoint */
00257               return(set_branch(pri,(box[pri][0]+box[pri][1])>>1));
00258                                    /* find median position */
00259        med >>= 1;
00260        for (c0 = box[pri][0]; med > 0; c0++)
00261               med -= t[c0];
00262        if (c0 > (box[pri][0]+box[pri][1])>>1)    /* if past the midpoint */
00263               c0--;                       /* part left of median */
00264        return(set_branch(pri,c0));
00265 #undef c0
00266 }