ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.12
Committed: Mon Sep 19 02:23:58 2005 UTC (18 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R8
Changes since 2.11: +2 -6 lines
Log Message:
Eliminated SPEED macro from makeall and source tree

File Contents

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