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

# User Rev Content
1 greg 2.4 /* Copyright (c) 1993 Regents of the University of California */
2 greg 2.1
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 greg 2.10 extern BYTE clrtab[256][3];
27 greg 2.1 /* 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 greg 2.2 /* define CLOSEST to get closest colors */
35 greg 2.4 #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 greg 2.1
43    
44 greg 2.10 new_histo(n) /* clear our histogram */
45     int n;
46 greg 2.1 {
47     bzero((char *)histo, sizeof(histo));
48 greg 2.10 return(0);
49 greg 2.1 }
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 greg 2.2 #ifdef CLOSEST
80     closest(ncolors); /* ensure colors picked are closest */
81     #endif
82 greg 2.9 /* reset dithering function */
83     dith_colrs((BYTE *)NULL, (COLR *)NULL, 0);
84 greg 2.1 /* 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 greg 2.8 static short (*cerr)[3] = NULL;
115 greg 2.1 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 greg 2.8 if (N) {
121     free((char *)cerr);
122     cerr = NULL;
123     }
124     if (n)
125 greg 2.1 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 greg 2.6 long t[HMAX], med;
190 greg 2.1 /* 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 greg 2.5 unsigned long sum[3];
246 greg 2.7 unsigned r, g;
247     unsigned long n;
248 greg 2.5 register unsigned b, c;
249 greg 2.1 /* 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 greg 2.7 if (n >= (1L<<23)/HMAX) { /* avoid overflow */
264 greg 2.5 sum[RED] /= n;
265     sum[GRN] /= n;
266     sum[BLU] /= n;
267     n = 1;
268     }
269 greg 2.1 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 greg 2.2
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 greg 2.3 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 greg 2.2 addneigh(neigh, histo[r][g][b], histo[r+1][g][b]);
305 greg 2.3 if (g < NGRN-1 && histo[r][g][b] != histo[r][g+1][b])
306 greg 2.2 addneigh(neigh, histo[r][g][b], histo[r][g+1][b]);
307 greg 2.3 if (b < NBLU-1 && histo[r][g][b] != histo[r][g][b+1])
308 greg 2.2 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 greg 2.4 register int tmp;
357 greg 2.2 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