ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.5
Committed: Tue Aug 2 13:22:08 1994 UTC (29 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +3 -3 lines
Log Message:
bug fixes from Tony

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 /* */
213 /* NeuQuant */
214 /* -------- */
215 /* */
216 /* Copyright: Anthony Dekker, June 1994 */
217 /* */
218 /* This program performs colour quantization of graphics images (SUN */
219 /* raster files). It uses a Kohonen Neural Network. It produces */
220 /* better results than existing methods and runs faster, using minimal */
221 /* space (8kB plus the image itself). The algorithm is described in */
222 /* the paper "Kohonen Neural Networks for Optimal Colour Quantization" */
223 /* to appear in the journal "Network: Computation in Neural Systems". */
224 /* It is a significant improvement of an earlier algorithm. */
225 /* */
226 /* This program is distributed free for academic use or for evaluation */
227 /* by commercial organizations. */
228 /* */
229 /* Usage: NeuQuant -n inputfile > outputfile */
230 /* */
231 /* where n is a sampling factor for neural learning. */
232 /* */
233 /* Program performance compared with other methods is as follows: */
234 /* */
235 /* Algorithm | Av. CPU Time | Quantization Error */
236 /* ------------------------------------------------------------- */
237 /* NeuQuant -3 | 314 | 5.55 */
238 /* NeuQuant -10 | 119 | 5.97 */
239 /* NeuQuant -30 | 65 | 6.53 */
240 /* Oct-Trees | 141 | 8.96 */
241 /* Median Cut (XV -best) | 420 | 9.28 */
242 /* Median Cut (XV -slow) | 72 | 12.15 */
243 /* */
244 /* Author's address: Dept of ISCS, National University of Singapore */
245 /* Kent Ridge, Singapore 0511 */
246 /* Email: [email protected] */
247 /*----------------------------------------------------------------------*/
248
249 #define bool int
250 #define false 0
251 #define true 1
252
253 #define initrad 32
254 #define radiusdec 30
255 #define alphadec 30
256
257 /* defs for freq and bias */
258 #define gammashift 10
259 #define betashift gammashift
260 #define intbiasshift 16
261 #define intbias (1<<intbiasshift)
262 #define gamma (1<<gammashift)
263 #define beta (intbias>>betashift)
264 #define betagamma (intbias<<(gammashift-betashift))
265 #define gammaphi (intbias<<(gammashift-8))
266
267 /* defs for rad and alpha */
268 #define maxrad (initrad+1)
269 #define radiusbiasshift 6
270 #define radiusbias (1<<radiusbiasshift)
271 #define initradius ((int) (initrad*radiusbias))
272 #define alphabiasshift 10
273 #define initalpha (1<<alphabiasshift)
274 #define radbiasshift 8
275 #define radbias (1<<radbiasshift)
276 #define alpharadbshift (alphabiasshift+radbiasshift)
277 #define alpharadbias (1<<alpharadbshift)
278
279 /* other defs */
280 #define netbiasshift 4
281 #define funnyshift (intbiasshift-netbiasshift)
282 #define maxnetval ((256<<netbiasshift)-1)
283 #define ncycles 100
284 #define jump1 499 /* prime */
285 #define jump2 491 /* prime */
286 #define jump3 487 /* any pic whose size was divisible by all */
287 #define jump4 503 /* four primes would be simply enormous */
288
289 /* cheater definitions (GW) */
290 #define thepicture thesamples
291 #define lengthcount (nsamples*3)
292 #define samplefac 1
293
294 typedef int pixel[4]; /* BGRc */
295
296 static pixel network[256];
297
298 static int netindex[256];
299
300 static int bias [256];
301 static int freq [256];
302 static int radpower[256]; /* actually need only go up to maxrad */
303
304 /* fixed space overhead 256*4+256+256+256+256 words = 256*8 = 8kB */
305
306
307 static
308 initnet()
309 {
310 register int i;
311 register int *p;
312
313 for (i=0; i<clrtabsiz; i++) {
314 p = network[i];
315 p[0] =
316 p[1] =
317 p[2] = (i<<(netbiasshift+8))/clrtabsiz;
318 freq[i] = intbias/clrtabsiz; /* 1/256 */
319 bias[i] = 0;
320 }
321 }
322
323
324 static
325 inxbuild()
326 {
327 register int i,j,smallpos,smallval;
328 register int *p,*q;
329 int start,previous;
330
331 previous = 0;
332 start = 0;
333 for (i=0; i<clrtabsiz; i++) {
334 p = network[i];
335 smallpos = i;
336 smallval = p[1]; /* index on g */
337 /* find smallest in i+1..clrtabsiz-1 */
338 for (j=i+1; j<clrtabsiz; j++) {
339 q = network[j];
340 if (q[1] < smallval) { /* index on g */
341 smallpos = j;
342 smallval = q[1]; /* index on g */
343 }
344 }
345 q = network[smallpos];
346 if (i != smallpos) {
347 j = q[0]; q[0] = p[0]; p[0] = j;
348 j = q[1]; q[1] = p[1]; p[1] = j;
349 j = q[2]; q[2] = p[2]; p[2] = j;
350 j = q[3]; q[3] = p[3]; p[3] = j;
351 }
352 /* smallval entry is now in position i */
353 if (smallval != previous) {
354 netindex[previous] = (start+i)>>1;
355 for (j=previous+1; j<smallval; j++) netindex[j] = i;
356 previous = smallval;
357 start = i;
358 }
359 }
360 netindex[previous] = (start+255)>>1;
361 for (j=previous+1; j<256; j++) netindex[j] = 255;
362 }
363
364
365 static int
366 inxsearch(b,g,r) /* accepts real BGR values after net is unbiased */
367 register int b,g,r;
368 {
369 register int i,j,best,x,y,bestd;
370 register int *p;
371
372 bestd = 1000; /* biggest possible dist is 256*3 */
373 best = -1;
374 i = netindex[g]; /* index on g */
375 j = i-1;
376
377 while ((i<clrtabsiz) || (j>=0)) {
378 if (i<clrtabsiz) {
379 p = network[i];
380 x = p[1] - g; /* inx key */
381 if (x >= bestd) i = clrtabsiz; /* stop iter */
382 else {
383 i++;
384 if (x<0) x = -x;
385 y = p[0] - b;
386 if (y<0) y = -y;
387 x += y;
388 if (x<bestd) {
389 y = p[2] - r;
390 if (y<0) y = -y;
391 x += y; /* x holds distance */
392 if (x<bestd) {bestd=x; best=p[3];}
393 }
394 }
395 }
396 if (j>=0) {
397 p = network[j];
398 x = g - p[1]; /* inx key - reverse dif */
399 if (x >= bestd) j = -1; /* stop iter */
400 else {
401 j--;
402 if (x<0) x = -x;
403 y = p[0] - b;
404 if (y<0) y = -y;
405 x += y;
406 if (x<bestd) {
407 y = p[2] - r;
408 if (y<0) y = -y;
409 x += y; /* x holds distance */
410 if (x<bestd) {bestd=x; best=p[3];}
411 }
412 }
413 }
414 }
415 return(best);
416 }
417
418
419 static int
420 contest(b,g,r) /* accepts biased BGR values */
421 register int b,g,r;
422 {
423 register int i,best,bestbias,x,y,bestd,bestbiasd;
424 register int *p,*q, *pp;
425
426 bestd = ~(1<<31);
427 bestbiasd = bestd;
428 best = -1;
429 bestbias = best;
430 q = bias;
431 p = freq;
432 for (i=0; i<clrtabsiz; i++) {
433 pp = network[i];
434 x = pp[0] - b;
435 if (x<0) x = -x;
436 y = pp[1] - g;
437 if (y<0) y = -y;
438 x += y;
439 y = pp[2] - r;
440 if (y<0) y = -y;
441 x += y; /* x holds distance */
442 /* >> netbiasshift not needed if funnyshift used */
443 if (x<bestd) {bestd=x; best=i;}
444 y = x - ((*q)>>funnyshift); /* y holds biasd */
445 if (y<bestbiasd) {bestbiasd=y; bestbias=i;}
446 y = (*p >> betashift); /* y holds beta*freq */
447 *p -= y;
448 *q += (y<<gammashift);
449 p++;
450 q++;
451 }
452 freq[best] += beta;
453 bias[best] -= betagamma;
454 return(bestbias);
455 }
456
457
458 static
459 alterneigh(rad,i,b,g,r) /* accepts biased BGR values */
460 int rad,i;
461 register int b,g,r;
462 {
463 register int j,k,lo,hi,a;
464 register int *p, *q;
465
466 lo = i-rad;
467 if (lo<-1) lo= -1;
468 hi = i+rad;
469 if (hi>clrtabsiz) hi=clrtabsiz;
470
471 j = i+1;
472 k = i-1;
473 q = radpower;
474 while ((j<hi) || (k>lo)) {
475 a = (*(++q));
476 if (j<hi) {
477 p = network[j];
478 *p -= (a*(*p - b)) / alpharadbias;
479 p++;
480 *p -= (a*(*p - g)) / alpharadbias;
481 p++;
482 *p -= (a*(*p - r)) / alpharadbias;
483 j++;
484 }
485 if (k>lo) {
486 p = network[k];
487 *p -= (a*(*p - b)) / alpharadbias;
488 p++;
489 *p -= (a*(*p - g)) / alpharadbias;
490 p++;
491 *p -= (a*(*p - r)) / alpharadbias;
492 k--;
493 }
494 }
495 }
496
497
498 static
499 altersingle(alpha,j,b,g,r) /* accepts biased BGR values */
500 register int alpha,j,b,g,r;
501 {
502 register int *q;
503
504 q = network[j]; /* alter hit neuron */
505 *q -= (alpha*(*q - b)) / initalpha;
506 q++;
507 *q -= (alpha*(*q - g)) / initalpha;
508 q++;
509 *q -= (alpha*(*q - r)) / initalpha;
510 }
511
512
513 static
514 learn()
515 {
516 register int i,j,b,g,r;
517 int radius,rad,alpha,step,delta,upto;
518 register unsigned char *p;
519 unsigned char *lim;
520
521 upto = lengthcount/(3*samplefac);
522 delta = upto/ncycles;
523 lim = thepicture + lengthcount;
524 p = thepicture;
525 alpha = initalpha;
526 radius = initradius;
527 rad = radius >> radiusbiasshift;
528 if (rad <= 1) rad = 0;
529 for (i=0; i<rad; i++)
530 radpower[i] = alpha*(((rad*rad - i*i)*radbias)/(rad*rad));
531
532 if ((lengthcount%jump1) != 0) step = 3*jump1;
533 else {
534 if ((lengthcount%jump2) !=0) step = 3*jump2;
535 else {
536 if ((lengthcount%jump3) !=0) step = 3*jump3;
537 else step = 3*jump4;
538 }
539 }
540 i = 0;
541 while (i < upto) {
542 b = p[0] << netbiasshift;
543 g = p[1] << netbiasshift;
544 r = p[2] << netbiasshift;
545 j = contest(b,g,r);
546
547 altersingle(alpha,j,b,g,r);
548 if (rad) alterneigh(rad,j,b,g,r);
549 /* alter neighbours */
550
551 p += step;
552 if (p >= lim) p -= lengthcount;
553
554 i++;
555 if (i%delta == 0) {
556 alpha -= alpha / alphadec;
557 radius -= radius / radiusdec;
558 rad = radius >> radiusbiasshift;
559 if (rad <= 1) rad = 0;
560 for (j=0; j<rad; j++)
561 radpower[j] = alpha*(((rad*rad - j*j)*radbias)/(rad*rad));
562 }
563 }
564 }
565
566 static
567 unbiasnet()
568 {
569 int i,j;
570
571 for (i=0; i<clrtabsiz; i++) {
572 for (j=0; j<3; j++)
573 network[i][j] >>= netbiasshift;
574 network[i][3] = i; /* record colour no */
575 }
576 }
577
578 /* Don't do this until the network has been unbiased */
579
580 static
581 cpyclrtab()
582 {
583 register int i,j,k;
584
585 for (j=0; j<clrtabsiz; j++) {
586 k = network[j][3];
587 for (i = 0; i < 3; i++)
588 clrtab[k][i] = network[j][2-i];
589 }
590 }