ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.8
Committed: Mon Dec 12 12:19:04 1994 UTC (29 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.7: +2 -0 lines
Log Message:
added static function definitions

File Contents

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