ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.9
Committed: Sat Feb 22 02:07:27 2003 UTC (21 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R5
Changes since 2.8: +3 -6 lines
Log Message:
Changes and check-in for 3.5 release
Includes new source files and modifications not recorded for many years
See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release

File Contents

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