ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.1
Committed: Fri Jun 10 12:33:35 1994 UTC (29 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial revision

File Contents

# User Rev Content
1 greg 2.1 /* Copyright (c) 1994 Regents of the University of California */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ LBL";
5     #endif
6    
7     /*
8     * Neural-Net quantization algorithm based on work of Anthony Dekker
9     */
10    
11     #include "standard.h"
12    
13     #include "color.h"
14    
15     #include "random.h"
16    
17     #ifdef COMPAT_MODE
18     #define neu_init new_histo
19     #define neu_pixel cnt_pixel
20     #define neu_colrs cnt_colrs
21     #define neu_clrtab new_clrtab
22     #define neu_map_pixel map_pixel
23     #define neu_map_colrs map_colrs
24     #define neu_dith_colrs dith_colrs
25     #endif
26     /* our color table (global) */
27     extern BYTE clrtab[256][3];
28     static int clrtabsiz;
29    
30     #ifndef DEFSMPFAC
31     #ifdef SPEED
32     #define DEFSMPFAC (240/SPEED+3)
33     #else
34     #define DEFSMPFAC 30
35     #endif
36     #endif
37    
38     int samplefac = DEFSMPFAC; /* sampling factor */
39    
40     /* Samples array starts off holding spacing between adjacent
41     * samples, and ends up holding actual BGR sample values.
42     */
43     static BYTE *thesamples;
44     static int nsamples;
45     static BYTE *cursamp;
46     static long skipcount;
47    
48     #define MAXSKIP (1<<24-1)
49    
50     #define nskip(sp) ((long)(sp)[0]<<16|(long)(sp)[1]<<8|(long)(sp)[2])
51    
52     #define setskip(sp,n) ((sp)[0]=(n)>>16,(sp)[1]=((n)>>8)&255,(sp)[2]=(n)&255)
53    
54    
55     neu_init(npixels) /* initialize our sample array */
56     long npixels;
57     {
58     register int nsleft;
59     register long sv;
60     double rval, cumprob;
61     long npleft;
62    
63     nsamples = npixels/samplefac;
64     if (nsamples < 600)
65     return(-1);
66     thesamples = (BYTE *)malloc((nsamples+1)*3);
67     if (thesamples == NULL)
68     return(-1);
69     cursamp = thesamples;
70     npleft = npixels;
71     nsleft = nsamples;
72     while (nsleft) {
73     rval = frandom(); /* random distance to next sample */
74     sv = 0;
75     cumprob = 0.;
76     while ((cumprob += (1.-cumprob)*nsleft/(npleft-sv)) < rval)
77     sv++;
78     setskip(cursamp, sv);
79     cursamp += 3;
80     npleft -= sv;
81     nsleft--;
82     }
83     setskip(cursamp, 0); /* dummy tagged onto end */
84     cursamp = thesamples;
85     skipcount = nskip(cursamp);
86     return(0);
87     }
88    
89    
90     neu_pixel(col) /* add pixel to our samples */
91     register BYTE col[];
92     {
93     if (!skipcount--) {
94     cursamp[0] = col[BLU];
95     cursamp[1] = col[GRN];
96     cursamp[2] = col[RED];
97     cursamp += 3;
98     skipcount = nskip(cursamp);
99     }
100     }
101    
102    
103     neu_colrs(cs, n) /* add a scanline to our samples */
104     register COLR *cs;
105     register int n;
106     {
107     while (n > skipcount) {
108     cs += skipcount;
109     cursamp[0] = cs[0][BLU];
110     cursamp[1] = cs[0][GRN];
111     cursamp[2] = cs[0][RED];
112     cs++;
113     n -= skipcount+1;
114     cursamp += 3;
115     skipcount = nskip(cursamp);
116     }
117     skipcount -= n;
118     }
119    
120    
121     neu_clrtab(ncolors) /* make new color table using ncolors */
122     int ncolors;
123     {
124     clrtabsiz = ncolors;
125     if (clrtabsiz > 256) clrtabsiz = 256;
126     initnet();
127     learn();
128     unbiasnet();
129     cpyclrtab();
130     inxbuild();
131     /* we're done with our samples */
132     free((char *)thesamples);
133     /* reset dithering function */
134     neu_dith_colrs((BYTE *)NULL, (COLR *)NULL, 0);
135     /* return new color table size */
136     return(clrtabsiz);
137     }
138    
139    
140     int
141     neu_map_pixel(col) /* get pixel for color */
142     register BYTE col[];
143     {
144     return(inxsearch(col[BLU],col[GRN],col[RED]));
145     }
146    
147    
148     neu_map_colrs(bs, cs, n) /* convert a scanline to color index values */
149     register BYTE *bs;
150     register COLR *cs;
151     register int n;
152     {
153     while (n-- > 0) {
154     *bs++ = inxsearch(cs[0][BLU],cs[0][GRN],cs[0][RED]);
155     cs++;
156     }
157     }
158    
159    
160     neu_dith_colrs(bs, cs, n) /* convert scanline to dithered index values */
161     register BYTE *bs;
162     register COLR *cs;
163     int n;
164     {
165     static short (*cerr)[3] = NULL;
166     static int N = 0;
167     int err[3], errp[3];
168     register int x, i;
169    
170     if (n != N) { /* get error propogation array */
171     if (N) {
172     free((char *)cerr);
173     cerr = NULL;
174     }
175     if (n)
176     cerr = (short (*)[3])malloc(3*n*sizeof(short));
177     if (cerr == NULL) {
178     N = 0;
179     map_colrs(bs, cs, n);
180     return;
181     }
182     N = n;
183     bzero((char *)cerr, 3*N*sizeof(short));
184     }
185     err[0] = err[1] = err[2] = 0;
186     for (x = 0; x < n; x++) {
187     for (i = 0; i < 3; i++) { /* dither value */
188     errp[i] = err[i];
189     err[i] += cerr[x][i];
190     #ifdef MAXERR
191     if (err[i] > MAXERR) err[i] = MAXERR;
192     else if (err[i] < -MAXERR) err[i] = -MAXERR;
193     #endif
194     err[i] += cs[x][i];
195     if (err[i] < 0) err[i] = 0;
196     else if (err[i] > 255) err[i] = 255;
197     }
198     bs[x] = inxsearch(err[BLU],err[GRN],err[RED]);
199     for (i = 0; i < 3; i++) { /* propagate error */
200     err[i] -= clrtab[bs[x]][i];
201     err[i] /= 3;
202     cerr[x][i] = err[i] + errp[i];
203     }
204     }
205     }
206    
207     /* The following was adapted and modified from the original (GW) */
208     /*----------------------------------------------------------------------*/
209     /* */
210     /* NeuQuant */
211     /* -------- */
212     /* */
213     /* Copyright: Anthony Dekker, June 1994 */
214     /* */
215     /* This program performs colour quantization of graphics images (SUN */
216     /* raster files). It uses a Kohonen Neural Network. It produces */
217     /* better results than existing methods and runs faster, using minimal */
218     /* space (8kB plus the image itself). The algorithm is described in */
219     /* the paper "Kohonen Neural Networks for Optimal Colour Quantization" */
220     /* to appear in the journal "Network: Computation in Neural Systems". */
221     /* It is a significant improvement of an earlier algorithm. */
222     /* */
223     /* This program is distributed free for academic use or for evaluation */
224     /* by commercial organizations. */
225     /* */
226     /* Usage: NeuQuant -n inputfile > outputfile */
227     /* */
228     /* where n is a sampling factor for neural learning. */
229     /* */
230     /* Program performance compared with other methods is as follows: */
231     /* */
232     /* Algorithm | Av. CPU Time | Quantization Error */
233     /* ------------------------------------------------------------- */
234     /* NeuQuant -3 | 314 | 5.55 */
235     /* NeuQuant -10 | 119 | 5.97 */
236     /* NeuQuant -30 | 65 | 6.53 */
237     /* Oct-Trees | 141 | 8.96 */
238     /* Median Cut (XV -best) | 420 | 9.28 */
239     /* Median Cut (XV -slow) | 72 | 12.15 */
240     /* */
241     /* Author's address: Dept of ISCS, National University of Singapore */
242     /* Kent Ridge, Singapore 0511 */
243     /* Email: [email protected] */
244     /*----------------------------------------------------------------------*/
245    
246     #define bool int
247     #define false 0
248     #define true 1
249    
250     #define initrad 32
251     #define radiusdec 30
252     #define alphadec; 30
253    
254     /* defs for freq and bias */
255     #define gammashift 10
256     #define betashift gammashift
257     #define intbiasshift 16
258     #define intbias (1<<intbiasshift)
259     #define gamma (1<<gammashift)
260     #define beta (intbias>>betashift)
261     #define betagamma (intbias<<(gammashift-betashift))
262     #define gammaphi (intbias<<(gammashift-8))
263    
264     /* defs for rad and alpha */
265     #define maxrad (initrad+1)
266     #define radiusbiasshift 6
267     #define radiusbias (1<<radiusbiasshift)
268     #define initradius ((int) (initrad*radiusbias))
269     #define alphabiasshift 10
270     #define initalpha (1<<alphabiasshift)
271     #define radbiasshift 8
272     #define radbias (1<<radbiasshift)
273     #define alpharadbshift (alphabiasshift+radbiasshift)
274     #define alpharadbias (1<<alpharadbshift)
275    
276     /* other defs */
277     #define netbiasshift 4
278     #define funnyshift (intbiasshift-netbiasshift)
279     #define maxnetval ((256<<netbiasshift)-1)
280     #define ncycles 100
281     #define jump1 499 /* prime */
282     #define jump2 491 /* prime */
283     #define jump3 487 /* any pic whose size was divisible by all */
284     #define jump4 503 /* four primes would be simply enormous */
285    
286     /* cheater definitions (GW) */
287     #define thepicture thesamples
288     #define lengthcount (nsamples*3)
289     #define samplefac 1
290    
291     typedef int pixel[4]; /* BGRc */
292    
293     static pixel network[256];
294    
295     static int netindex[256];
296    
297     static int bias [256];
298     static int freq [256];
299     static int radpower[256]; /* actually need only go up to maxrad */
300    
301     /* fixed space overhead 256*4+256+256+256+256 words = 256*8 = 8kB */
302    
303    
304     static
305     initnet()
306     {
307     register int i;
308     register int *p;
309    
310     for (i=0; i<clrtabsiz; i++) {
311     p = network[i];
312     p[0] = i << netbiasshift;
313     p[1] = i << netbiasshift;
314     p[2] = i << netbiasshift;
315     freq[i] = intbias >> 8; /* 1/256 */
316     bias[i] = 0;
317     }
318     }
319    
320    
321     static
322     inxbuild()
323     {
324     register int i,j,smallpos,smallval;
325     register int *p,*q;
326     int start,previous;
327    
328     previous = 0;
329     start = 0;
330     for (i=0; i<clrtabsiz; i++) {
331     p = network[i];
332     smallpos = i;
333     smallval = p[1]; /* index on g */
334     /* find smallest in i+1..clrtabsiz-1 */
335     for (j=i+1; j<clrtabsiz; j++) {
336     q = network[j];
337     if (q[1] < smallval) { /* index on g */
338     smallpos = j;
339     smallval = q[1]; /* index on g */
340     }
341     }
342     q = network[smallpos];
343     if (i != smallpos) {
344     j = q[0]; q[0] = p[0]; p[0] = j;
345     j = q[1]; q[1] = p[1]; p[1] = j;
346     j = q[2]; q[2] = p[2]; p[2] = j;
347     j = q[3]; q[3] = p[3]; p[3] = j;
348     }
349     /* smallval entry is now in position i */
350     if (smallval != previous) {
351     netindex[previous] = (start+i)>>1;
352     for (j=previous+1; j<smallval; j++) netindex[j] = i;
353     previous = smallval;
354     start = i;
355     }
356     }
357     netindex[previous] = (start+clrtabsiz-1)>>1;
358     for (j=previous+1; j<clrtabsiz; j++) netindex[j] = clrtabsiz-1;
359     }
360    
361    
362     static int
363     inxsearch(b,g,r) /* accepts real BGR values after net is unbiased */
364     register int b,g,r;
365     {
366     register int i,j,best,x,y,bestd;
367     register int *p;
368    
369     bestd = 1000; /* biggest possible dist is 256*3 */
370     best = -1;
371     i = netindex[g]; /* index on g */
372     j = i-1;
373    
374     while ((i<clrtabsiz) || (j>=0)) {
375     if (i<clrtabsiz) {
376     p = network[i];
377     x = p[1] - g; /* inx key */
378     if (x >= bestd) i = clrtabsiz; /* stop iter */
379     else {
380     i++;
381     if (x<0) x = -x;
382     y = p[0] - b;
383     if (y<0) y = -y;
384     x += y;
385     if (x<bestd) {
386     y = p[2] - r;
387     if (y<0) y = -y;
388     x += y; /* x holds distance */
389     if (x<bestd) {bestd=x; best=p[3];}
390     }
391     }
392     }
393     if (j>=0) {
394     p = network[j];
395     x = g - p[1]; /* inx key - reverse dif */
396     if (x >= bestd) j = -1; /* stop iter */
397     else {
398     j--;
399     if (x<0) x = -x;
400     y = p[0] - b;
401     if (y<0) y = -y;
402     x += y;
403     if (x<bestd) {
404     y = p[2] - r;
405     if (y<0) y = -y;
406     x += y; /* x holds distance */
407     if (x<bestd) {bestd=x; best=p[3];}
408     }
409     }
410     }
411     }
412     return(best);
413     }
414    
415    
416     static int
417     contest(b,g,r) /* accepts biased BGR values */
418     register int b,g,r;
419     {
420     register int i,best,bestbias,x,y,bestd,bestbiasd;
421     register int *p,*q, *pp;
422    
423     bestd = ~(1<<31);
424     bestbiasd = bestd;
425     best = -1;
426     bestbias = best;
427     q = bias;
428     p = freq;
429     for (i=0; i<clrtabsiz; i++) {
430     pp = network[i];
431     x = pp[0] - b;
432     if (x<0) x = -x;
433     y = pp[1] - g;
434     if (y<0) y = -y;
435     x += y;
436     y = pp[2] - r;
437     if (y<0) y = -y;
438     x += y; /* x holds distance */
439     /* >> netbiasshift not needed if funnyshift used */
440     if (x<bestd) {bestd=x; best=i;}
441     y = x - ((*q)>>funnyshift); /* y holds biasd */
442     if (y<bestbiasd) {bestbiasd=y; bestbias=i;}
443     y = (*p >> betashift); /* y holds beta*freq */
444     *p -= y;
445     *q += (y<<gammashift);
446     p++;
447     q++;
448     }
449     freq[best] += beta;
450     bias[best] -= betagamma;
451     return(bestbias);
452     }
453    
454    
455     static
456     alterneigh(rad,i,b,g,r) /* accepts biased BGR values */
457     int rad,i;
458     register int b,g,r;
459     {
460     register int j,k,lo,hi,a;
461     register int *p, *q;
462    
463     lo = i-rad;
464     if (lo<-1) lo= -1;
465     hi = i+rad;
466     if (hi>clrtabsiz) hi=clrtabsiz;
467    
468     j = i+1;
469     k = i-1;
470     q = radpower;
471     while ((j<hi) || (k>lo)) {
472     a = (*(++q));
473     if (j<hi) {
474     p = network[j];
475     *p -= (a*(*p - b)) / alpharadbias;
476     p++;
477     *p -= (a*(*p - g)) / alpharadbias;
478     p++;
479     *p -= (a*(*p - r)) / alpharadbias;
480     j++;
481     }
482     if (k>lo) {
483     p = network[k];
484     *p -= (a*(*p - b)) / alpharadbias;
485     p++;
486     *p -= (a*(*p - g)) / alpharadbias;
487     p++;
488     *p -= (a*(*p - r)) / alpharadbias;
489     k--;
490     }
491     }
492     }
493    
494    
495     static
496     altersingle(alpha,j,b,g,r) /* accepts biased BGR values */
497     register int alpha,j,b,g,r;
498     {
499     register int *q;
500    
501     q = network[j]; /* alter hit neuron */
502     *q -= (alpha*(*q - b)) / initalpha;
503     q++;
504     *q -= (alpha*(*q - g)) / initalpha;
505     q++;
506     *q -= (alpha*(*q - r)) / initalpha;
507     }
508    
509    
510     static
511     learn()
512     {
513     register int i,j,b,g,r;
514     int radius,rad,alpha,step,delta,upto;
515     register unsigned char *p;
516     unsigned char *lim;
517    
518     upto = lengthcount/(3*samplefac);
519     delta = upto/ncycles;
520     lim = thepicture + lengthcount;
521     p = thepicture;
522     alpha = initalpha;
523     radius = initradius;
524     rad = radius >> radiusbiasshift;
525     if (rad <= 1) rad = 0;
526     for (i=0; i<rad; i++)
527     radpower[i] = alpha*(((rad*rad - i*i)*radbias)/(rad*rad));
528    
529     if ((lengthcount%jump1) != 0) step = 3*jump1;
530     else {
531     if ((lengthcount%jump2) !=0) step = 3*jump2;
532     else {
533     if ((lengthcount%jump3) !=0) step = 3*jump3;
534     else step = 3*jump4;
535     }
536     }
537     i = 0;
538     while (i < upto) {
539     b = p[0] << netbiasshift;
540     g = p[1] << netbiasshift;
541     r = p[2] << netbiasshift;
542     j = contest(b,g,r);
543    
544     altersingle(alpha,j,b,g,r);
545     if (rad) alterneigh(rad,j,b,g,r);
546     /* alter neighbours */
547    
548     p += step;
549     if (p >= lim) p -= lengthcount;
550    
551     i++;
552     if (i%delta == 0) {
553     alpha -= alpha / alphadec;
554     radius -= radius / radiusdec;
555     rad = radius >> radiusbiasshift;
556     if (rad <= 1) rad = 0;
557     for (j=0; j<rad; j++)
558     radpower[j] = alpha*(((rad*rad - j*j)*radbias)/(rad*rad));
559     }
560     }
561     }
562    
563     static
564     unbiasnet()
565     {
566     int i,j;
567    
568     for (i=0; i<clrtabsiz; i++) {
569     for (j=0; j<3; j++)
570     network[i][j] >>= netbiasshift;
571     network[i][3] = i; /* record colour no */
572     }
573     }
574    
575     /* Don't do this until the network has been unbiased */
576    
577     static
578     cpyclrtab()
579     {
580     register int i,j,k;
581    
582     for (j=0; j<clrtabsiz; j++) {
583     k = network[j][3];
584     for (i = 0; i < 3; i++)
585     clrtab[k][i] = network[j][2-i];
586     }
587     }