ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/neuclrtab.c
(Generate patch)

Comparing ray/src/px/neuclrtab.c (file contents):
Revision 2.5 by greg, Tue Aug 2 13:22:08 1994 UTC vs.
Revision 2.10 by schorsch, Mon Jun 30 14:59:12 2003 UTC

# Line 1 | Line 1
1 /* Copyright (c) 1994 Regents of the University of California */
2
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6
4   /*
5   * Neural-Net quantization algorithm based on work of Anthony Dekker
6   */
7  
8 < #include "standard.h"
8 > #include "copyright.h"
9  
10 < #include "color.h"
10 > #include <string.h>
11  
12 + #include "standard.h"
13 + #include "color.h"
14   #include "random.h"
15  
16   #ifdef COMPAT_MODE
# Line 51 | Line 50 | static long    skipcount;
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   {
# Line 132 | Line 133 | int    ncolors;
133          cpyclrtab();
134          inxbuild();
135                                  /* we're done with our samples */
136 <        free((char *)thesamples);
136 >        free((void *)thesamples);
137                                  /* reset dithering function */
138          neu_dith_colrs((BYTE *)NULL, (COLR *)NULL, 0);
139                                  /* return new color table size */
# Line 172 | Line 173 | int    n;
173  
174          if (n != N) {           /* get error propogation array */
175                  if (N) {
176 <                        free((char *)cerr);
176 >                        free((void *)cerr);
177                          cerr = NULL;
178                  }
179                  if (n)
# Line 183 | Line 184 | int    n;
184                          return;
185                  }
186                  N = n;
187 <                bzero((char *)cerr, 3*N*sizeof(short));
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++) {
# Line 208 | Line 209 | int    n;
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, June 1994                    */
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       */
# Line 246 | Line 253 | int    n;
253   /* Email:       [email protected]                                     */
254   /*----------------------------------------------------------------------*/
255  
256 < #define bool    int
257 < #define false   0
258 < #define true    1
256 > #define bool            int
257 > #define false           0
258 > #define true            1
259  
260 < #define initrad                 32
261 < #define radiusdec               30
262 < #define alphadec                30
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 gammashift      10
268 < #define betashift       gammashift
269 < #define intbiasshift    16
270 < #define intbias         (1<<intbiasshift)
271 < #define gamma           (1<<gammashift)
272 < #define beta            (intbias>>betashift)
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))
265 #define gammaphi        (intbias<<(gammashift-8))
274  
275 < /* defs for rad and alpha */
276 < #define maxrad          (initrad+1)
277 < #define radiusbiasshift 6
278 < #define radiusbias      (1<<radiusbiasshift)
279 < #define initradius      ((int) (initrad*radiusbias))
280 < #define alphabiasshift  10
281 < #define initalpha       (1<<alphabiasshift)
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         (1<<radbiasshift)
289 > #define radbias         (((int) 1)<<radbiasshift)
290   #define alpharadbshift  (alphabiasshift+radbiasshift)
291 < #define alpharadbias    (1<<alpharadbshift)
291 > #define alpharadbias    (((int) 1)<<alpharadbshift)
292  
293 < /* other defs */
294 < #define netbiasshift    4
295 < #define funnyshift      (intbiasshift-netbiasshift)
296 < #define maxnetval       ((256<<netbiasshift)-1)
297 < #define ncycles         100
298 < #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 */
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  
289 /* cheater definitions (GW) */
290 #define thepicture      thesamples
291 #define lengthcount     (nsamples*3)
292 #define samplefac       1
293
300   typedef int pixel[4];  /* BGRc */
301 + pixel network[256];
302  
303 < static pixel network[256];
303 > int netindex[256];      /* for network lookup - really 256 */
304  
305 < static int netindex[256];
305 > int bias [256];         /* bias and freq arrays for learning */
306 > int freq [256];
307 > int radpower[initrad];  /* radpower for precomputation */
308  
300 static int bias [256];
301 static int freq [256];
302 static int radpower[256];       /* actually need only go up to maxrad */
309  
310 < /* fixed space overhead 256*4+256+256+256+256 words = 256*8 = 8kB */
310 > /* initialise network in range (0,0,0) to (255,255,255) */
311  
312 <
307 < static
308 < initnet()
312 > initnet()      
313   {
314          register int i;
315          register int *p;
316          
317 <        for (i=0; i<clrtabsiz; i++) {
317 >        for (i=0; i<netsize; i++) {
318                  p = network[i];
319 <                p[0] =
320 <                p[1] =
317 <                p[2] = (i<<(netbiasshift+8))/clrtabsiz;
318 <                freq[i] = intbias/clrtabsiz;  /* 1/256 */
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 < static
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 start,previous;
332 >        int previouscol,startpos;
333  
334 <        previous = 0;
335 <        start = 0;
336 <        for (i=0; i<clrtabsiz; i++) {
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+1..clrtabsiz-1 */
341 <                for (j=i+1; j<clrtabsiz; j++) {
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;
# Line 343 | Line 346 | inxbuild()
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;
# Line 350 | Line 354 | inxbuild()
354                          j = q[3];   q[3] = p[3];   p[3] = j;
355                  }
356                  /* smallval entry is now in position i */
357 <                if (smallval != previous) {
358 <                        netindex[previous] = (start+i)>>1;
359 <                        for (j=previous+1; j<smallval; j++) netindex[j] = i;
360 <                        previous = smallval;
361 <                        start = 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[previous] = (start+255)>>1;
365 <        for (j=previous+1; j<256; j++) netindex[j] = 255;
364 >        netindex[previouscol] = (startpos+maxnetpos)>>1;
365 >        for (j=previouscol+1; j<256; j++) netindex[j] = maxnetpos; /* really 256 */
366   }
367  
368  
369 < static int
366 < inxsearch(b,g,r)  /* accepts real BGR values after net is unbiased */
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,best,x,y,bestd;
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;
379 >        j = i-1;         /* start at netindex[g] and work outwards */
380  
381 <        while ((i<clrtabsiz) || (j>=0)) {
382 <                if (i<clrtabsiz) {
381 >        while ((i<netsize) || (j>=0)) {
382 >                if (i<netsize) {
383                          p = network[i];
384 <                        x = p[1] - g;   /* inx key */
385 <                        if (x >= bestd) i = clrtabsiz; /* stop iter */
384 >                        dist = p[1] - g;        /* inx key */
385 >                        if (dist >= bestd) i = netsize; /* stop iter */
386                          else {
387                                  i++;
388 <                                if (x<0) x = -x;
389 <                                y = p[0] - b;
390 <                                if (y<0) y = -y;
391 <                                x += y;
392 <                                if (x<bestd) {
393 <                                        y = p[2] - r;  
394 <                                        if (y<0) y = -y;
391 <                                        x += y; /* x holds distance */
392 <                                        if (x<bestd) {bestd=x; best=p[3];}
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 <                        x = g - p[1]; /* inx key - reverse dif */
401 <                        if (x >= bestd) j = -1; /* stop iter */
400 >                        dist = g - p[1]; /* inx key - reverse dif */
401 >                        if (dist >= bestd) j = -1; /* stop iter */
402                          else {
403                                  j--;
404 <                                if (x<0) x = -x;
405 <                                y = p[0] - b;
406 <                                if (y<0) y = -y;
407 <                                x += y;
408 <                                if (x<bestd) {
409 <                                        y = p[2] - r;  
410 <                                        if (y<0) y = -y;
409 <                                        x += y; /* x holds distance */
410 <                                        if (x<bestd) {bestd=x; best=p[3];}
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                  }
# Line 416 | Line 416 | register int b,g,r;
416   }
417  
418  
419 < static int
420 < contest(b,g,r)  /* accepts biased BGR values */
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,best,bestbias,x,y,bestd,bestbiasd;
428 <        register int *p,*q, *pp;
427 >        register int i,dist,a,biasdist,betafreq;
428 >        int bestpos,bestbiaspos,bestd,bestbiasd;
429 >        register int *p,*f, *n;
430  
431 <        bestd = ~(1<<31);
431 >        bestd = ~(((int) 1)<<31);
432          bestbiasd = bestd;
433 <        best = -1;
434 <        bestbias = best;
435 <        q = bias;
436 <        p = freq;
437 <        for (i=0; i<clrtabsiz; i++) {
438 <                pp = network[i];
439 <                x = pp[0] - b;
440 <                if (x<0) x = -x;
441 <                y = pp[1] - g;
442 <                if (y<0) y = -y;
443 <                x += y;
444 <                y = pp[2] - r;  
445 <                if (y<0) y = -y;
446 <                x += y; /* x holds distance */
447 <                        /* >> netbiasshift not needed if funnyshift used */
448 <                if (x<bestd) {bestd=x; best=i;}
449 <                y = x - ((*q)>>funnyshift);  /* y holds biasd */
450 <                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++;
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[best] += beta;
453 <        bias[best] -= betagamma;
454 <        return(bestbias);
452 >        freq[bestpos] += beta;
453 >        bias[bestpos] -= betagamma;
454 >        return(bestbiaspos);
455   }
456  
457  
458 < static
459 < alterneigh(rad,i,b,g,r) /* accepts biased BGR values */
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;
485 <        if (lo<-1) lo= -1;
468 <        hi = i+rad;
469 <        if (hi>clrtabsiz) hi=clrtabsiz;
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;
# Line 495 | Line 511 | register int b,g,r;
511   }
512  
513  
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;
517 >        int radius,rad,alpha,step,delta,samplepixels;
518          register unsigned char *p;
519          unsigned char *lim;
520  
521 <        upto = lengthcount/(3*samplefac);
522 <        delta = upto/ncycles;
523 <        lim = thepicture + lengthcount;
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%jump1) != 0) step = 3*jump1;
533 >        
534 >        if ((lengthcount%prime1) != 0) step = 3*prime1;
535          else {
536 <                if ((lengthcount%jump2) !=0) step = 3*jump2;
536 >                if ((lengthcount%prime2) !=0) step = 3*prime2;
537                  else {
538 <                        if ((lengthcount%jump3) !=0) step = 3*jump3;
539 <                        else step = 3*jump4;
538 >                        if ((lengthcount%prime3) !=0) step = 3*prime3;
539 >                        else step = 3*prime4;
540                  }
541          }
542 +        
543          i = 0;
544 <        while (i < upto) {
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);
549 <                                                /* alter neighbours */
551 >                if (rad) alterneigh(rad,j,b,g,r);   /* alter neighbours */
552  
553                  p += step;
554                  if (p >= lim) p -= lengthcount;
# Line 563 | Line 565 | learn()
565          }
566   }
567          
568 < static
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<clrtabsiz; i++) {
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 < /* Don't do this until the network has been unbiased */
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<clrtabsiz; j++) {
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];

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines