ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.10
Committed: Mon Jun 30 14:59:12 2003 UTC (20 years, 10 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 2.9: +6 -4 lines
Log Message:
Replaced most outdated BSD function calls with their posix equivalents, and cleaned up a few other platform dependencies.

File Contents

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