ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.11
Committed: Sun Mar 28 20:33:14 2004 UTC (20 years ago) by schorsch
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad3R6, rad3R6P1
Changes since 2.10: +82 -37 lines
Log Message:
Continued ANSIfication, and other fixes and clarifications.

File Contents

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