ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/clrtab.c
Revision: 2.10
Committed: Fri Jun 10 12:51:16 1994 UTC (29 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.9: +4 -2 lines
Log Message:
added Anthony Dekker's neural net color quantization code

File Contents

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