src/rt/colortab.c

/* Copyright (c) 1989 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 * colortab.c - allocate and control dynamic color table.
 *
 *      We start off with a uniform partition of color space.
 *      As pixels are sent to the frame buffer, a histogram is built.
 *      When a new color table is requested, the histogram is used
 *      to make a pseudo-optimal partition, after which the
 *      histogram is cleared.  This algorithm
 *      performs only as well as the next drawing's color
 *      distribution is correlated to the last.
 */

#include "color.h"

#define NULL            0
                                /* histogram resolution */
#define NRED            24
#define NGRN            32
#define NBLU            16
#define HMAX            NGRN
                                /* minimum box count for adaptive partition */
#define MINSAMP         7
                                /* maximum distance^2 before color reassign */
#define MAXDST2         5
                                /* maximum frame buffer depth */
#define FBDEPTH         8
                                /* color map resolution */
#define MAPSIZ          256
                                /* map a color */
#define map_col(c,p)    clrmap[ colval(c,p)<1. ? \
                                (int)(colval(c,p)*MAPSIZ) : MAPSIZ-1 ]
                                /* color partition tree */
#define CNODE           short
#define set_branch(p,c) ((c)<<2|(p))
#define set_pval(pv)    ((pv)<<2|3)
#define is_pval(cn)     (((cn)&3)==3)
#define part(cn)        ((cn)>>2)
#define prim(cn)        ((cn)&3)
#define pval(cn)        ((cn)>>2)
                                /* our color table */
struct tabent {
        long    sum[3];         /* sum of colors using this entry */
        long    n;              /* number of colors */
        short   ent[3];         /* current table value */
}       clrtab[1<<FBDEPTH];
                                /* our color correction map */
static BYTE     clrmap[MAPSIZ];
                                /* histogram of colors used */
static unsigned histo[NRED][NGRN][NBLU];
                                /* initial color cube boundaries */
static int      CLRCUBE[3][2] = {0,NRED,0,NGRN,0,NBLU};
                                /* color cube partition */
static CNODE    ctree[1<<(FBDEPTH+1)];
                                /* callback for pixel assignment */
static int      (*set_pixel)();


int
new_ctab(ncolors, cset)         /* start new color table with max ncolors */
int     ncolors;
int     (*cset)();
{
        if (ncolors < 1 || ncolors > 1<<FBDEPTH || cset == NULL)
                return(0);
                                /* assign pixel callback routine */
        set_pixel = cset;
                                /* clear color table */
        bzero(clrtab, sizeof(clrtab));
                                /* partition color space */
        cut(ctree, FBDEPTH, CLRCUBE, 0, ncolors);
                                /* clear histogram */
        bzero(histo, sizeof(histo));
                                /* return number of colors used */
        return(ncolors);
}


int
get_pixel(col)                  /* get pixel for color */
COLOR   col;
{
        int     r, g, b;
        int     cv[3];
        register union { CNODE *t; struct tabent *e; }  p;
        register int    h;
                                                /* map color */
        r = map_col(col,RED);
        g = map_col(col,GRN);
        b = map_col(col,BLU);
                                                /* reduce resolution */
        cv[RED] = (r*NRED)>>8;
        cv[GRN] = (g*NGRN)>>8;
        cv[BLU] = (b*NBLU)>>8;
                                                /* add to histogram */
        histo[cv[RED]][cv[GRN]][cv[BLU]]++;
                                                /* find pixel in tree */
        p.t = ctree;
        for (h = FBDEPTH; h > 0; h--) {
                if (is_pval(*p.t))
                        break;
                if (cv[prim(*p.t)] < part(*p.t))
                        p.t++;          /* left branch */
                else
                        p.t += 1<<h;    /* right branch */
        }
        h = pval(*p.t);
                                                /* add to color table */
        p.e = clrtab + h;
                                        /* add to sum */
        p.e->sum[RED] += r;
        p.e->sum[GRN] += g;
        p.e->sum[BLU] += b;
        p.e->n++;
                                        /* recompute average */
        r = p.e->sum[RED] / p.e->n;
        g = p.e->sum[GRN] / p.e->n;
        b = p.e->sum[BLU] / p.e->n;
                                        /* check for movement */
        if (p.e->n == 1 ||
                        (r-p.e->ent[RED])*(r-p.e->ent[RED]) +
                        (g-p.e->ent[GRN])*(g-p.e->ent[GRN]) +
                        (b-p.e->ent[BLU])*(b-p.e->ent[BLU]) > MAXDST2) {
                p.e->ent[RED] = r;
                p.e->ent[GRN] = g;      /* reassign pixel */
                p.e->ent[BLU] = b;
#ifdef notdef
                printf("pixel %d = (%d,%d,%d) (%d refs)\n",
                                h, r, g, b, p.e->n);
#endif
                (*set_pixel)(h, r, g, b);
        }
        return(h);                              /* return pixel value */
}


make_cmap(gam)                  /* make gamma correction map */
double  gam;
{
        extern double   pow();
        register int    i;
        
        for (i = 0; i < MAPSIZ; i++)
                clrmap[i] = 256.0 * pow(i/(double)MAPSIZ, 1.0/gam);
}


static
cut(tree, height, box, c0, c1)          /* partition color space */
register CNODE  *tree;
int     height;
register int    box[3][2];
int     c0, c1;
{
        int     kb[3][2];
        
        if (c1-c0 <= 1) {               /* assign pixel */
                *tree = set_pval(c0);
                return;
        }
                                        /* split box */
        *tree = split(box);
        bcopy(box, kb, sizeof(kb));
                                                /* do left (lesser) branch */
        kb[prim(*tree)][1] = part(*tree);
        cut(tree+1, height-1, kb, c0, (c0+c1)>>1);
                                                /* do right branch */
        kb[prim(*tree)][0] = part(*tree);
        kb[prim(*tree)][1] = box[prim(*tree)][1];
        cut(tree+(1<<height), height-1, kb, (c0+c1)>>1, c1);
}


static int
split(box)                              /* find median cut for box */
register int    box[3][2];
{
#define c0      r
        register int    r, g, b;
        int     pri;
        int     t[HMAX], med;
                                        /* find dominant axis */
        pri = RED;
        if (box[GRN][1]-box[GRN][0] > box[pri][1]-box[pri][0])
                pri = GRN;
        if (box[BLU][1]-box[BLU][0] > box[pri][1]-box[pri][0])
                pri = BLU;
                                        /* sum histogram over box */
        med = 0;
        switch (pri) {
        case RED:
                for (r = box[RED][0]; r < box[RED][1]; r++) {
                        t[r] = 0;
                        for (g = box[GRN][0]; g < box[GRN][1]; g++)
                                for (b = box[BLU][0]; b < box[BLU][1]; b++)
                                        t[r] += histo[r][g][b];
                        med += t[r];
                }
                break;
        case GRN:
                for (g = box[GRN][0]; g < box[GRN][1]; g++) {
                        t[g] = 0;
                        for (b = box[BLU][0]; b < box[BLU][1]; b++)
                                for (r = box[RED][0]; r < box[RED][1]; r++)
                                        t[g] += histo[r][g][b];
                        med += t[g];
                }
                break;
        case BLU:
                for (b = box[BLU][0]; b < box[BLU][1]; b++) {
                        t[b] = 0;
                        for (r = box[RED][0]; r < box[RED][1]; r++)
                                for (g = box[GRN][0]; g < box[GRN][1]; g++)
                                        t[b] += histo[r][g][b];
                        med += t[b];
                }
                break;
        }
        if (med < MINSAMP)              /* if too sparse, split at midpoint */
                return(set_branch(pri,(box[pri][0]+box[pri][1])>>1));
                                        /* find median position */
        med >>= 1;
        for (c0 = box[pri][0]; med > 0; c0++)
                med -= t[c0];
        if (c0 > (box[pri][0]+box[pri][1])>>1)  /* if past the midpoint */
                c0--;                           /* part left of median */
        return(set_branch(pri,c0));
#undef c0
}
Revision:	1.3
Committed:	Tue Oct 3 11:10:10 1989 UTC (34 years, 7 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+75 -50 lines
Log Message:	Improved color table allocation, second cut.
#	User	Rev	Content
1	greg	1.1	/* Copyright (c) 1989 Regents of the University of California */
2
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8			* colortab.c - allocate and control dynamic color table.
9			*
10			* We start off with a uniform partition of color space.
11			* As pixels are sent to the frame buffer, a histogram is built.
12			* When a new color table is requested, the histogram is used
13			* to make a pseudo-optimal partition, after which the
14			* histogram is cleared. This algorithm
15			* performs only as well as the next drawing's color
16			* distribution is correlated to the last.
17			*/
18
19			#include "color.h"
20
21			#define NULL 0
22			/* histogram resolution */
23	greg	1.3	#define NRED 24
24	greg	1.1	#define NGRN 32
25	greg	1.3	#define NBLU 16
26	greg	1.1	#define HMAX NGRN
27			/* minimum box count for adaptive partition */
28			#define MINSAMP 7
29	greg	1.3	/* maximum distance^2 before color reassign */
30			#define MAXDST2 5
31	greg	1.1	/* maximum frame buffer depth */
32			#define FBDEPTH 8
33			/* color map resolution */
34	greg	1.3	#define MAPSIZ 256
35	greg	1.1	/* map a color */
36	greg	1.3	#define map_col(c,p) clrmap[ colval(c,p)<1. ? \
37	greg	1.1	(int)(colval(c,p)*MAPSIZ) : MAPSIZ-1 ]
38			/* color partition tree */
39			#define CNODE short
40			#define set_branch(p,c) ((c)<<2\|(p))
41			#define set_pval(pv) ((pv)<<2\|3)
42			#define is_pval(cn) (((cn)&3)==3)
43			#define part(cn) ((cn)>>2)
44			#define prim(cn) ((cn)&3)
45			#define pval(cn) ((cn)>>2)
46			/* our color table */
47	greg	1.3	struct tabent {
48			long sum[3]; /* sum of colors using this entry */
49			long n; /* number of colors */
50			short ent[3]; /* current table value */
51			} clrtab[1<<FBDEPTH];
52	greg	1.1	/* our color correction map */
53	greg	1.3	static BYTE clrmap[MAPSIZ];
54	greg	1.1	/* histogram of colors used */
55			static unsigned histo[NRED][NGRN][NBLU];
56			/* initial color cube boundaries */
57			static int CLRCUBE[3][2] = {0,NRED,0,NGRN,0,NBLU};
58			/* color cube partition */
59			static CNODE ctree[1<<(FBDEPTH+1)];
60	greg	1.3	/* callback for pixel assignment */
61			static int (*set_pixel)();
62	greg	1.1
63
64	greg	1.3	int
65			new_ctab(ncolors, cset) /* start new color table with max ncolors */
66	greg	1.1	int ncolors;
67	greg	1.3	int (*cset)();
68	greg	1.1	{
69	greg	1.3	if (ncolors < 1 \|\| ncolors > 1<<FBDEPTH \|\| cset == NULL)
70			return(0);
71			/* assign pixel callback routine */
72			set_pixel = cset;
73			/* clear color table */
74			bzero(clrtab, sizeof(clrtab));
75			/* partition color space */
76	greg	1.1	cut(ctree, FBDEPTH, CLRCUBE, 0, ncolors);
77			/* clear histogram */
78			bzero(histo, sizeof(histo));
79	greg	1.3	/* return number of colors used */
80			return(ncolors);
81	greg	1.1	}
82
83
84			int
85			get_pixel(col) /* get pixel for color */
86			COLOR col;
87			{
88	greg	1.3	int r, g, b;
89	greg	1.1	int cv[3];
90	greg	1.3	register union { CNODE t; struct tabent e; } p;
91	greg	1.1	register int h;
92	greg	1.3	/* map color */
93			r = map_col(col,RED);
94			g = map_col(col,GRN);
95			b = map_col(col,BLU);
96			/* reduce resolution */
97			cv[RED] = (r*NRED)>>8;
98			cv[GRN] = (g*NGRN)>>8;
99			cv[BLU] = (b*NBLU)>>8;
100			/* add to histogram */
101	greg	1.1	histo[cv[RED]][cv[GRN]][cv[BLU]]++;
102	greg	1.3	/* find pixel in tree */
103			p.t = ctree;
104	greg	1.1	for (h = FBDEPTH; h > 0; h--) {
105	greg	1.3	if (is_pval(*p.t))
106	greg	1.1	break;
107	greg	1.3	if (cv[prim(p.t)] < part(p.t))
108			p.t++; /* left branch */
109	greg	1.1	else
110	greg	1.3	p.t += 1<<h; /* right branch */
111	greg	1.1	}
112	greg	1.3	h = pval(*p.t);
113			/* add to color table */
114			p.e = clrtab + h;
115			/* add to sum */
116			p.e->sum[RED] += r;
117			p.e->sum[GRN] += g;
118			p.e->sum[BLU] += b;
119			p.e->n++;
120			/* recompute average */
121			r = p.e->sum[RED] / p.e->n;
122			g = p.e->sum[GRN] / p.e->n;
123			b = p.e->sum[BLU] / p.e->n;
124			/* check for movement */
125			if (p.e->n == 1 \|\|
126			(r-p.e->ent[RED])*(r-p.e->ent[RED]) +
127			(g-p.e->ent[GRN])*(g-p.e->ent[GRN]) +
128			(b-p.e->ent[BLU])*(b-p.e->ent[BLU]) > MAXDST2) {
129			p.e->ent[RED] = r;
130			p.e->ent[GRN] = g; /* reassign pixel */
131			p.e->ent[BLU] = b;
132	greg	1.1	#ifdef notdef
133	greg	1.3	printf("pixel %d = (%d,%d,%d) (%d refs)\n",
134			h, r, g, b, p.e->n);
135	greg	1.1	#endif
136	greg	1.3	(*set_pixel)(h, r, g, b);
137			}
138			return(h); /* return pixel value */
139	greg	1.1	}
140
141
142			make_cmap(gam) /* make gamma correction map */
143			double gam;
144			{
145			extern double pow();
146			register int i;
147
148	greg	1.3	for (i = 0; i < MAPSIZ; i++)
149			clrmap[i] = 256.0 * pow(i/(double)MAPSIZ, 1.0/gam);
150	greg	1.1	}
151
152
153			static
154			cut(tree, height, box, c0, c1) /* partition color space */
155			register CNODE *tree;
156			int height;
157			register int box[3][2];
158			int c0, c1;
159			{
160			int kb[3][2];
161
162	greg	1.3	if (c1-c0 <= 1) { /* assign pixel */
163	greg	1.1	*tree = set_pval(c0);
164			return;
165			}
166			/* split box */
167			*tree = split(box);
168			bcopy(box, kb, sizeof(kb));
169	greg	1.3	/* do left (lesser) branch */
170	greg	1.1	kb[prim(tree)][1] = part(tree);
171	greg	1.3	cut(tree+1, height-1, kb, c0, (c0+c1)>>1);
172			/* do right branch */
173	greg	1.1	kb[prim(tree)][0] = part(tree);
174			kb[prim(tree)][1] = box[prim(tree)][1];
175	greg	1.3	cut(tree+(1<<height), height-1, kb, (c0+c1)>>1, c1);
176	greg	1.1	}
177
178
179			static int
180			split(box) /* find median cut for box */
181			register int box[3][2];
182			{
183			#define c0 r
184			register int r, g, b;
185			int pri;
186			int t[HMAX], med;
187			/* find dominant axis */
188			pri = RED;
189			if (box[GRN][1]-box[GRN][0] > box[pri][1]-box[pri][0])
190			pri = GRN;
191			if (box[BLU][1]-box[BLU][0] > box[pri][1]-box[pri][0])
192			pri = BLU;
193			/* sum histogram over box */
194			med = 0;
195			switch (pri) {
196			case RED:
197			for (r = box[RED][0]; r < box[RED][1]; r++) {
198			t[r] = 0;
199			for (g = box[GRN][0]; g < box[GRN][1]; g++)
200			for (b = box[BLU][0]; b < box[BLU][1]; b++)
201			t[r] += histo[r][g][b];
202			med += t[r];
203			}
204			break;
205			case GRN:
206			for (g = box[GRN][0]; g < box[GRN][1]; g++) {
207			t[g] = 0;
208			for (b = box[BLU][0]; b < box[BLU][1]; b++)
209			for (r = box[RED][0]; r < box[RED][1]; r++)
210			t[g] += histo[r][g][b];
211			med += t[g];
212			}
213			break;
214			case BLU:
215			for (b = box[BLU][0]; b < box[BLU][1]; b++) {
216			t[b] = 0;
217			for (r = box[RED][0]; r < box[RED][1]; r++)
218			for (g = box[GRN][0]; g < box[GRN][1]; g++)
219			t[b] += histo[r][g][b];
220			med += t[b];
221			}
222			break;
223			}
224			if (med < MINSAMP) /* if too sparse, split at midpoint */
225			return(set_branch(pri,(box[pri][0]+box[pri][1])>>1));
226			/* find median position */
227			med >>= 1;
228			for (c0 = box[pri][0]; med > 0; c0++)
229			med -= t[c0];
230			if (c0 > (box[pri][0]+box[pri][1])>>1) /* if past the midpoint */
231			c0--; /* part left of median */
232			return(set_branch(pri,c0));
233			#undef c0
234			}