ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.1
Committed: Fri Jun 10 12:33:35 1994 UTC (29 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial revision

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