ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.7
Committed: Tue Nov 22 12:19:21 1994 UTC (29 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.6: +19 -19 lines
Log Message:
bug fix

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