ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
Revision: 2.14
Committed: Fri May 20 02:06:39 2011 UTC (12 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, rad5R2, rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R1, rad4R2P1, rad5R3, HEAD
Changes since 2.13: +10 -10 lines
Log Message:
Changed every instance of BYTE to uby8 to avoid conflicts

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: neuclrtab.c,v 2.13 2007/09/08 19:17:52 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 #include "clrtab.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 uby8 clrtab[256][3];
28 static int clrtabsiz;
29
30 #ifndef DEFSMPFAC
31 #define DEFSMPFAC 3
32 #endif
33
34 int samplefac = DEFSMPFAC; /* sampling factor */
35
36 /* Samples array starts off holding spacing between adjacent
37 * samples, and ends up holding actual BGR sample values.
38 */
39 static uby8 *thesamples;
40 static int nsamples;
41 static uby8 *cursamp;
42 static long skipcount;
43
44 #define MAXSKIP (1<<24-1)
45
46 #define nskip(sp) ((long)(sp)[0]<<16|(long)(sp)[1]<<8|(long)(sp)[2])
47
48 #define setskip(sp,n) ((sp)[0]=(n)>>16,(sp)[1]=((n)>>8)&255,(sp)[2]=(n)&255)
49
50 static void initnet(void);
51 static void inxbuild(void);
52 static int inxsearch(int b, int g, int r);
53 static int contest(int b, int g, int r);
54 static void altersingle(int alpha, int i, int b, int g, int r);
55 static void alterneigh(int rad, int i, int b, int g, int r);
56 static void learn(void);
57 static void unbiasnet(void);
58 static void cpyclrtab(void);
59
60
61 extern int
62 neu_init( /* initialize our sample array */
63 long npixels
64 )
65 {
66 register int nsleft;
67 register long sv;
68 double rval, cumprob;
69 long npleft;
70
71 nsamples = npixels/samplefac;
72 if (nsamples < 600)
73 return(-1);
74 thesamples = (uby8 *)malloc(nsamples*3);
75 if (thesamples == NULL)
76 return(-1);
77 cursamp = thesamples;
78 npleft = npixels;
79 nsleft = nsamples;
80 while (nsleft) {
81 rval = frandom(); /* random distance to next sample */
82 sv = 0;
83 cumprob = 0.;
84 while ((cumprob += (1.-cumprob)*nsleft/(npleft-sv)) < rval)
85 sv++;
86 if (nsleft == nsamples)
87 skipcount = sv;
88 else {
89 setskip(cursamp, sv);
90 cursamp += 3;
91 }
92 npleft -= sv+1;
93 nsleft--;
94 }
95 setskip(cursamp, npleft); /* tag on end to skip the rest */
96 cursamp = thesamples;
97 return(0);
98 }
99
100
101 extern void
102 neu_pixel( /* add pixel to our samples */
103 register uby8 col[]
104 )
105 {
106 if (!skipcount--) {
107 skipcount = nskip(cursamp);
108 cursamp[0] = col[BLU];
109 cursamp[1] = col[GRN];
110 cursamp[2] = col[RED];
111 cursamp += 3;
112 }
113 }
114
115
116 extern void
117 neu_colrs( /* add a scanline to our samples */
118 register COLR *cs,
119 register int n
120 )
121 {
122 while (n > skipcount) {
123 cs += skipcount;
124 n -= skipcount+1;
125 skipcount = nskip(cursamp);
126 cursamp[0] = cs[0][BLU];
127 cursamp[1] = cs[0][GRN];
128 cursamp[2] = cs[0][RED];
129 cs++;
130 cursamp += 3;
131 }
132 skipcount -= n;
133 }
134
135
136 extern int
137 neu_clrtab( /* make new color table using ncolors */
138 int ncolors
139 )
140 {
141 clrtabsiz = ncolors;
142 if (clrtabsiz > 256) clrtabsiz = 256;
143 initnet();
144 learn();
145 unbiasnet();
146 cpyclrtab();
147 inxbuild();
148 /* we're done with our samples */
149 free((void *)thesamples);
150 /* reset dithering function */
151 neu_dith_colrs((uby8 *)NULL, (COLR *)NULL, 0);
152 /* return new color table size */
153 return(clrtabsiz);
154 }
155
156
157 extern int
158 neu_map_pixel( /* get pixel for color */
159 register uby8 col[]
160 )
161 {
162 return(inxsearch(col[BLU],col[GRN],col[RED]));
163 }
164
165
166 extern void
167 neu_map_colrs( /* convert a scanline to color index values */
168 register uby8 *bs,
169 register COLR *cs,
170 register int n
171 )
172 {
173 while (n-- > 0) {
174 *bs++ = inxsearch(cs[0][BLU],cs[0][GRN],cs[0][RED]);
175 cs++;
176 }
177 }
178
179
180 extern void
181 neu_dith_colrs( /* convert scanline to dithered index values */
182 register uby8 *bs,
183 register COLR *cs,
184 int n
185 )
186 {
187 static short (*cerr)[3] = NULL;
188 static int N = 0;
189 int err[3], errp[3];
190 register int x, i;
191
192 if (n != N) { /* get error propogation array */
193 if (N) {
194 free((void *)cerr);
195 cerr = NULL;
196 }
197 if (n)
198 cerr = (short (*)[3])malloc(3*n*sizeof(short));
199 if (cerr == NULL) {
200 N = 0;
201 map_colrs(bs, cs, n);
202 return;
203 }
204 N = n;
205 memset((char *)cerr, '\0', 3*N*sizeof(short));
206 }
207 err[0] = err[1] = err[2] = 0;
208 for (x = 0; x < n; x++) {
209 for (i = 0; i < 3; i++) { /* dither value */
210 errp[i] = err[i];
211 err[i] += cerr[x][i];
212 #ifdef MAXERR
213 if (err[i] > MAXERR) err[i] = MAXERR;
214 else if (err[i] < -MAXERR) err[i] = -MAXERR;
215 #endif
216 err[i] += cs[x][i];
217 if (err[i] < 0) err[i] = 0;
218 else if (err[i] > 255) err[i] = 255;
219 }
220 bs[x] = inxsearch(err[BLU],err[GRN],err[RED]);
221 for (i = 0; i < 3; i++) { /* propagate error */
222 err[i] -= clrtab[bs[x]][i];
223 err[i] /= 3;
224 cerr[x][i] = err[i] + errp[i];
225 }
226 }
227 }
228
229 /* The following was adapted and modified from the original (GW) */
230
231 /* cheater definitions (GW) */
232 #define thepicture thesamples
233 #define lengthcount (nsamples*3)
234 #define samplefac 1
235
236 /* NeuQuant Neural-Net Quantization Algorithm Interface
237 * ----------------------------------------------------
238 *
239 * Copyright (c) 1994 Anthony Dekker
240 *
241 * NEUQUANT Neural-Net quantization algorithm by Anthony Dekker, 1994.
242 * See "Kohonen neural networks for optimal colour quantization"
243 * in "Network: Computation in Neural Systems" Vol. 5 (1994) pp 351-367.
244 * for a discussion of the algorithm.
245 * See also http://members.ozemail.com.au/~dekker/NEUQUANT.HTML
246 *
247 * Any party obtaining a copy of these files from the author, directly or
248 * indirectly, is granted, free of charge, a full and unrestricted irrevocable,
249 * world-wide, paid up, royalty-free, nonexclusive right and license to deal
250 * in this software and documentation files (the "Software"), including without
251 * limitation the rights to use, copy, modify, merge, publish, distribute, sublicense,
252 * and/or sell copies of the Software, and to permit persons who receive
253 * copies from any such party to do so, with the only requirement being
254 * that this copyright notice remain intact.
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 static void
314 initnet(void)
315 {
316 register int i;
317 register int *p;
318
319 for (i=0; i<netsize; i++) {
320 p = network[i];
321 p[0] = p[1] = p[2] = (i << (netbiasshift+8))/netsize;
322 freq[i] = intbias/netsize; /* 1/netsize */
323 bias[i] = 0;
324 }
325 }
326
327
328 /* do after unbias - insertion sort of network and build netindex[0..255] */
329
330 static void
331 inxbuild(void)
332 {
333 register int i,j,smallpos,smallval;
334 register int *p,*q;
335 int previouscol,startpos;
336
337 previouscol = 0;
338 startpos = 0;
339 for (i=0; i<netsize; i++) {
340 p = network[i];
341 smallpos = i;
342 smallval = p[1]; /* index on g */
343 /* find smallest in i..netsize-1 */
344 for (j=i+1; j<netsize; j++) {
345 q = network[j];
346 if (q[1] < smallval) { /* index on g */
347 smallpos = j;
348 smallval = q[1]; /* index on g */
349 }
350 }
351 q = network[smallpos];
352 /* swap p (i) and q (smallpos) entries */
353 if (i != smallpos) {
354 j = q[0]; q[0] = p[0]; p[0] = j;
355 j = q[1]; q[1] = p[1]; p[1] = j;
356 j = q[2]; q[2] = p[2]; p[2] = j;
357 j = q[3]; q[3] = p[3]; p[3] = j;
358 }
359 /* smallval entry is now in position i */
360 if (smallval != previouscol) {
361 netindex[previouscol] = (startpos+i)>>1;
362 for (j=previouscol+1; j<smallval; j++) netindex[j] = i;
363 previouscol = smallval;
364 startpos = i;
365 }
366 }
367 netindex[previouscol] = (startpos+maxnetpos)>>1;
368 for (j=previouscol+1; j<256; j++) netindex[j] = maxnetpos; /* really 256 */
369 }
370
371
372 static int
373 inxsearch( /* accepts real BGR values after net is unbiased */
374 register int b,
375 register int g,
376 register int r
377 )
378 {
379 register int i,j,dist,a,bestd;
380 register int *p;
381 int best;
382
383 bestd = 1000; /* biggest possible dist is 256*3 */
384 best = -1;
385 i = netindex[g]; /* index on g */
386 j = i-1; /* start at netindex[g] and work outwards */
387
388 while ((i<netsize) || (j>=0)) {
389 if (i<netsize) {
390 p = network[i];
391 dist = p[1] - g; /* inx key */
392 if (dist >= bestd) i = netsize; /* stop iter */
393 else {
394 i++;
395 if (dist<0) dist = -dist;
396 a = p[0] - b; if (a<0) a = -a;
397 dist += a;
398 if (dist<bestd) {
399 a = p[2] - r; if (a<0) a = -a;
400 dist += a;
401 if (dist<bestd) {bestd=dist; best=p[3];}
402 }
403 }
404 }
405 if (j>=0) {
406 p = network[j];
407 dist = g - p[1]; /* inx key - reverse dif */
408 if (dist >= bestd) j = -1; /* stop iter */
409 else {
410 j--;
411 if (dist<0) dist = -dist;
412 a = p[0] - b; if (a<0) a = -a;
413 dist += a;
414 if (dist<bestd) {
415 a = p[2] - r; if (a<0) a = -a;
416 dist += a;
417 if (dist<bestd) {bestd=dist; best=p[3];}
418 }
419 }
420 }
421 }
422 return(best);
423 }
424
425
426 /* finds closest neuron (min dist) and updates freq */
427 /* finds best neuron (min dist-bias) and returns position */
428 /* for frequently chosen neurons, freq[i] is high and bias[i] is negative */
429 /* bias[i] = gamma*((1/netsize)-freq[i]) */
430
431 static int
432 contest( /* accepts biased BGR values */
433 register int b,
434 register int g,
435 register int r
436 )
437 {
438 register int i,dist,a,biasdist,betafreq;
439 int bestpos,bestbiaspos,bestd,bestbiasd;
440 register int *p,*f, *n;
441
442 bestd = ~(((int) 1)<<31);
443 bestbiasd = bestd;
444 bestpos = -1;
445 bestbiaspos = bestpos;
446 p = bias;
447 f = freq;
448
449 for (i=0; i<netsize; i++) {
450 n = network[i];
451 dist = n[0] - b; if (dist<0) dist = -dist;
452 a = n[1] - g; if (a<0) a = -a;
453 dist += a;
454 a = n[2] - r; if (a<0) a = -a;
455 dist += a;
456 if (dist<bestd) {bestd=dist; bestpos=i;}
457 biasdist = dist - ((*p)>>(intbiasshift-netbiasshift));
458 if (biasdist<bestbiasd) {bestbiasd=biasdist; bestbiaspos=i;}
459 betafreq = (*f >> betashift);
460 *f++ -= betafreq;
461 *p++ += (betafreq<<gammashift);
462 }
463 freq[bestpos] += beta;
464 bias[bestpos] -= betagamma;
465 return(bestbiaspos);
466 }
467
468
469 /* move neuron i towards (b,g,r) by factor alpha */
470
471 static void
472 altersingle( /* accepts biased BGR values */
473 register int alpha,
474 register int i,
475 register int b,
476 register int g,
477 register int r
478 )
479 {
480 register int *n;
481
482 n = network[i]; /* alter hit neuron */
483 *n -= (alpha*(*n - b)) / initalpha;
484 n++;
485 *n -= (alpha*(*n - g)) / initalpha;
486 n++;
487 *n -= (alpha*(*n - r)) / initalpha;
488 }
489
490
491 /* move neurons adjacent to i towards (b,g,r) by factor */
492 /* alpha*(1-((i-j)^2/[r]^2)) precomputed as radpower[|i-j|]*/
493
494 static void
495 alterneigh( /* accents biased BGR values */
496 int rad,
497 int i,
498 register int b,
499 register int g,
500 register int r
501 )
502 {
503 register int j,k,lo,hi,a;
504 register int *p, *q;
505
506 lo = i-rad; if (lo<-1) lo= -1;
507 hi = i+rad; if (hi>netsize) hi=netsize;
508
509 j = i+1;
510 k = i-1;
511 q = radpower;
512 while ((j<hi) || (k>lo)) {
513 a = (*(++q));
514 if (j<hi) {
515 p = network[j];
516 *p -= (a*(*p - b)) / alpharadbias;
517 p++;
518 *p -= (a*(*p - g)) / alpharadbias;
519 p++;
520 *p -= (a*(*p - r)) / alpharadbias;
521 j++;
522 }
523 if (k>lo) {
524 p = network[k];
525 *p -= (a*(*p - b)) / alpharadbias;
526 p++;
527 *p -= (a*(*p - g)) / alpharadbias;
528 p++;
529 *p -= (a*(*p - r)) / alpharadbias;
530 k--;
531 }
532 }
533 }
534
535
536 static void
537 learn(void)
538 {
539 register int i,j,b,g,r;
540 int radius,rad,alpha,step,delta,samplepixels;
541 register unsigned char *p;
542 unsigned char *lim;
543
544 alphadec = 30 + ((samplefac-1)/3);
545 p = thepicture;
546 lim = thepicture + lengthcount;
547 samplepixels = lengthcount/(3*samplefac);
548 delta = samplepixels/ncycles;
549 alpha = initalpha;
550 radius = initradius;
551
552 rad = radius >> radiusbiasshift;
553 if (rad <= 1) rad = 0;
554 for (i=0; i<rad; i++)
555 radpower[i] = alpha*(((rad*rad - i*i)*radbias)/(rad*rad));
556
557 if ((lengthcount%prime1) != 0) step = 3*prime1;
558 else {
559 if ((lengthcount%prime2) !=0) step = 3*prime2;
560 else {
561 if ((lengthcount%prime3) !=0) step = 3*prime3;
562 else step = 3*prime4;
563 }
564 }
565
566 i = 0;
567 while (i < samplepixels) {
568 b = p[0] << netbiasshift;
569 g = p[1] << netbiasshift;
570 r = p[2] << netbiasshift;
571 j = contest(b,g,r);
572
573 altersingle(alpha,j,b,g,r);
574 if (rad) alterneigh(rad,j,b,g,r); /* alter neighbours */
575
576 p += step;
577 if (p >= lim) p -= lengthcount;
578
579 i++;
580 if (i%delta == 0) {
581 alpha -= alpha / alphadec;
582 radius -= radius / radiusdec;
583 rad = radius >> radiusbiasshift;
584 if (rad <= 1) rad = 0;
585 for (j=0; j<rad; j++)
586 radpower[j] = alpha*(((rad*rad - j*j)*radbias)/(rad*rad));
587 }
588 }
589 }
590
591 /* unbias network to give 0..255 entries */
592 /* which can then be used for colour map */
593 /* and record position i to prepare for sort */
594
595 static void
596 unbiasnet(void)
597 {
598 int i,j;
599
600 for (i=0; i<netsize; i++) {
601 for (j=0; j<3; j++)
602 network[i][j] >>= netbiasshift;
603 network[i][3] = i; /* record colour no */
604 }
605 }
606
607
608 /* Don't do this until the network has been unbiased (GW) */
609
610 static void
611 cpyclrtab(void)
612 {
613 register int i,j,k;
614
615 for (j=0; j<netsize; j++) {
616 k = network[j][3];
617 for (i = 0; i < 3; i++)
618 clrtab[k][i] = network[j][2-i];
619 }
620 }