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.4 by greg, Fri Jul 1 15:03:37 1994 UTC vs.
Revision 2.9 by greg, Sat Feb 22 02:07:27 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   */
# Line 51 | Line 48 | static long    skipcount;
48  
49   #define setskip(sp,n)   ((sp)[0]=(n)>>16,(sp)[1]=((n)>>8)&255,(sp)[2]=(n)&255)
50  
51 + static  cpyclrtab();
52  
53 +
54   neu_init(npixels)               /* initialize our sample array */
55   long    npixels;
56   {
# Line 132 | Line 131 | int    ncolors;
131          cpyclrtab();
132          inxbuild();
133                                  /* we're done with our samples */
134 <        free((char *)thesamples);
134 >        free((void *)thesamples);
135                                  /* reset dithering function */
136          neu_dith_colrs((BYTE *)NULL, (COLR *)NULL, 0);
137                                  /* return new color table size */
# Line 172 | Line 171 | int    n;
171  
172          if (n != N) {           /* get error propogation array */
173                  if (N) {
174 <                        free((char *)cerr);
174 >                        free((void *)cerr);
175                          cerr = NULL;
176                  }
177                  if (n)
# Line 208 | Line 207 | int    n;
207   }
208  
209   /* The following was adapted and modified from the original (GW)        */
210 +
211 + /* cheater definitions (GW) */
212 + #define thepicture      thesamples
213 + #define lengthcount     (nsamples*3)
214 + #define samplefac       1
215 +
216   /*----------------------------------------------------------------------*/
217   /*                                                                      */
218   /*                              NeuQuant                                */
219   /*                              --------                                */
220   /*                                                                      */
221 < /*              Copyright: Anthony Dekker, June 1994                    */
221 > /*              Copyright: Anthony Dekker, November 1994                */
222   /*                                                                      */
223   /* This program performs colour quantization of graphics images (SUN    */
224   /* raster files).  It uses a Kohonen Neural Network.  It produces       */
# Line 246 | Line 251 | int    n;
251   /* Email:       [email protected]                                     */
252   /*----------------------------------------------------------------------*/
253  
254 < #define bool    int
255 < #define false   0
256 < #define true    1
254 > #define bool            int
255 > #define false           0
256 > #define true            1
257  
258 < #define initrad                 32
259 < #define radiusdec               30
260 < #define alphadec                30
258 > /* network defs */
259 > #define netsize         clrtabsiz               /* number of colours - can change this */
260 > #define maxnetpos       (netsize-1)
261 > #define netbiasshift    4                       /* bias for colour values */
262 > #define ncycles         100                     /* no. of learning cycles */
263  
264   /* defs for freq and bias */
265 < #define gammashift      10
266 < #define betashift       gammashift
267 < #define intbiasshift    16
268 < #define intbias         (1<<intbiasshift)
269 < #define gamma           (1<<gammashift)
270 < #define beta            (intbias>>betashift)
265 > #define intbiasshift    16                      /* bias for fractions */
266 > #define intbias         (((int) 1)<<intbiasshift)
267 > #define gammashift      10                      /* gamma = 1024 */
268 > #define gamma           (((int) 1)<<gammashift)
269 > #define betashift       10
270 > #define beta            (intbias>>betashift)    /* beta = 1/1024 */
271   #define betagamma       (intbias<<(gammashift-betashift))
265 #define gammaphi        (intbias<<(gammashift-8))
272  
273 < /* defs for rad and alpha */
274 < #define maxrad          (initrad+1)
275 < #define radiusbiasshift 6
276 < #define radiusbias      (1<<radiusbiasshift)
277 < #define initradius      ((int) (initrad*radiusbias))
278 < #define alphabiasshift  10
279 < #define initalpha       (1<<alphabiasshift)
273 > /* defs for decreasing radius factor */
274 > #define initrad         (256>>3)                /* for 256 cols, radius starts */
275 > #define radiusbiasshift 6                       /* at 32.0 biased by 6 bits */
276 > #define radiusbias      (((int) 1)<<radiusbiasshift)
277 > #define initradius      (initrad*radiusbias)    /* and decreases by a */
278 > #define radiusdec       30                      /* factor of 1/30 each cycle */
279 >
280 > /* defs for decreasing alpha factor */
281 > #define alphabiasshift  10                      /* alpha starts at 1.0 */
282 > #define initalpha       (((int) 1)<<alphabiasshift)
283 > int alphadec;                                   /* biased by 10 bits */
284 >
285 > /* radbias and alpharadbias used for radpower calculation */
286   #define radbiasshift    8
287 < #define radbias         (1<<radbiasshift)
287 > #define radbias         (((int) 1)<<radbiasshift)
288   #define alpharadbshift  (alphabiasshift+radbiasshift)
289 < #define alpharadbias    (1<<alpharadbshift)
289 > #define alpharadbias    (((int) 1)<<alpharadbshift)
290  
291 < /* other defs */
292 < #define netbiasshift    4
293 < #define funnyshift      (intbiasshift-netbiasshift)
294 < #define maxnetval       ((256<<netbiasshift)-1)
295 < #define ncycles         100
296 < #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 */
291 > /* four primes near 500 - assume no image has a length so large */
292 > /* that it is divisible by all four primes */
293 > #define prime1          499
294 > #define prime2          491
295 > #define prime3          487
296 > #define prime4          503
297  
289 /* cheater definitions (GW) */
290 #define thepicture      thesamples
291 #define lengthcount     (nsamples*3)
292 #define samplefac       1
293
298   typedef int pixel[4];  /* BGRc */
299 + pixel network[256];
300  
301 < static pixel network[256];
301 > int netindex[256];      /* for network lookup - really 256 */
302  
303 < static int netindex[256];
303 > int bias [256];         /* bias and freq arrays for learning */
304 > int freq [256];
305 > int radpower[initrad];  /* radpower for precomputation */
306  
300 static int bias [256];
301 static int freq [256];
302 static int radpower[256];       /* actually need only go up to maxrad */
307  
308 < /* fixed space overhead 256*4+256+256+256+256 words = 256*8 = 8kB */
308 > /* initialise network in range (0,0,0) to (255,255,255) */
309  
310 <
307 < static
308 < initnet()
310 > initnet()      
311   {
312          register int i;
313          register int *p;
314          
315 <        for (i=0; i<clrtabsiz; i++) {
315 >        for (i=0; i<netsize; i++) {
316                  p = network[i];
317 <                p[0] =
318 <                p[1] =
317 <                p[2] = (i<<8) / clrtabsiz;
318 <                freq[i] = intbias/clrtabsiz;  /* 1/256 */
317 >                p[0] = p[1] = p[2] = (i << (netbiasshift+8))/netsize;
318 >                freq[i] = intbias/netsize;  /* 1/netsize */
319                  bias[i] = 0;
320          }
321   }
322  
323  
324 < static
324 > /* do after unbias - insertion sort of network and build netindex[0..255] */
325 >
326   inxbuild()
327   {
328          register int i,j,smallpos,smallval;
329          register int *p,*q;
330 <        int start,previous;
330 >        int previouscol,startpos;
331  
332 <        previous = 0;
333 <        start = 0;
334 <        for (i=0; i<clrtabsiz; i++) {
332 >        previouscol = 0;
333 >        startpos = 0;
334 >        for (i=0; i<netsize; i++) {
335                  p = network[i];
336                  smallpos = i;
337                  smallval = p[1];        /* index on g */
338 <                /* find smallest in i+1..clrtabsiz-1 */
339 <                for (j=i+1; j<clrtabsiz; j++) {
338 >                /* find smallest in i..netsize-1 */
339 >                for (j=i+1; j<netsize; j++) {
340                          q = network[j];
341                          if (q[1] < smallval) {  /* index on g */
342                                  smallpos = j;
# Line 343 | Line 344 | inxbuild()
344                          }
345                  }
346                  q = network[smallpos];
347 +                /* swap p (i) and q (smallpos) entries */
348                  if (i != smallpos) {
349                          j = q[0];   q[0] = p[0];   p[0] = j;
350                          j = q[1];   q[1] = p[1];   p[1] = j;
# Line 350 | Line 352 | inxbuild()
352                          j = q[3];   q[3] = p[3];   p[3] = j;
353                  }
354                  /* smallval entry is now in position i */
355 <                if (smallval != previous) {
356 <                        netindex[previous] = (start+i)>>1;
357 <                        for (j=previous+1; j<smallval; j++) netindex[j] = i;
358 <                        previous = smallval;
359 <                        start = i;
355 >                if (smallval != previouscol) {
356 >                        netindex[previouscol] = (startpos+i)>>1;
357 >                        for (j=previouscol+1; j<smallval; j++) netindex[j] = i;
358 >                        previouscol = smallval;
359 >                        startpos = i;
360                  }
361          }
362 <        netindex[previous] = (start+clrtabsiz-1)>>1;
363 <        for (j=previous+1; j<clrtabsiz; j++) netindex[j] = clrtabsiz-1;
362 >        netindex[previouscol] = (startpos+maxnetpos)>>1;
363 >        for (j=previouscol+1; j<256; j++) netindex[j] = maxnetpos; /* really 256 */
364   }
365  
366  
367 < static int
366 < inxsearch(b,g,r)  /* accepts real BGR values after net is unbiased */
367 > int inxsearch(b,g,r)  /* accepts real BGR values after net is unbiased */
368   register int b,g,r;
369   {
370 <        register int i,j,best,x,y,bestd;
370 >        register int i,j,dist,a,bestd;
371          register int *p;
372 +        int best;
373  
374          bestd = 1000;   /* biggest possible dist is 256*3 */
375          best = -1;
376          i = netindex[g]; /* index on g */
377 <        j = i-1;
377 >        j = i-1;         /* start at netindex[g] and work outwards */
378  
379 <        while ((i<clrtabsiz) || (j>=0)) {
380 <                if (i<clrtabsiz) {
379 >        while ((i<netsize) || (j>=0)) {
380 >                if (i<netsize) {
381                          p = network[i];
382 <                        x = p[1] - g;   /* inx key */
383 <                        if (x >= bestd) i = clrtabsiz; /* stop iter */
382 >                        dist = p[1] - g;        /* inx key */
383 >                        if (dist >= bestd) i = netsize; /* stop iter */
384                          else {
385                                  i++;
386 <                                if (x<0) x = -x;
387 <                                y = p[0] - b;
388 <                                if (y<0) y = -y;
389 <                                x += y;
390 <                                if (x<bestd) {
391 <                                        y = p[2] - r;  
392 <                                        if (y<0) y = -y;
391 <                                        x += y; /* x holds distance */
392 <                                        if (x<bestd) {bestd=x; best=p[3];}
386 >                                if (dist<0) dist = -dist;
387 >                                a = p[0] - b;   if (a<0) a = -a;
388 >                                dist += a;
389 >                                if (dist<bestd) {
390 >                                        a = p[2] - r;   if (a<0) a = -a;
391 >                                        dist += a;
392 >                                        if (dist<bestd) {bestd=dist; 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 */
398 >                        dist = g - p[1]; /* inx key - reverse dif */
399 >                        if (dist >= 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];}
402 >                                if (dist<0) dist = -dist;
403 >                                a = p[0] - b;   if (a<0) a = -a;
404 >                                dist += a;
405 >                                if (dist<bestd) {
406 >                                        a = p[2] - r;   if (a<0) a = -a;
407 >                                        dist += a;
408 >                                        if (dist<bestd) {bestd=dist; best=p[3];}
409                                  }
410                          }
411                  }
# Line 416 | Line 414 | register int b,g,r;
414   }
415  
416  
417 < static int
418 < contest(b,g,r)  /* accepts biased BGR values */
417 > /* finds closest neuron (min dist) and updates freq */
418 > /* finds best neuron (min dist-bias) and returns position */
419 > /* for frequently chosen neurons, freq[i] is high and bias[i] is negative */
420 > /* bias[i] = gamma*((1/netsize)-freq[i]) */
421 >
422 > int contest(b,g,r)      /* accepts biased BGR values */
423   register int b,g,r;
424   {
425 <        register int i,best,bestbias,x,y,bestd,bestbiasd;
426 <        register int *p,*q, *pp;
425 >        register int i,dist,a,biasdist,betafreq;
426 >        int bestpos,bestbiaspos,bestd,bestbiasd;
427 >        register int *p,*f, *n;
428  
429 <        bestd = ~(1<<31);
429 >        bestd = ~(((int) 1)<<31);
430          bestbiasd = bestd;
431 <        best = -1;
432 <        bestbias = best;
433 <        q = bias;
434 <        p = freq;
435 <        for (i=0; i<clrtabsiz; i++) {
436 <                pp = network[i];
437 <                x = pp[0] - b;
438 <                if (x<0) x = -x;
439 <                y = pp[1] - g;
440 <                if (y<0) y = -y;
441 <                x += y;
442 <                y = pp[2] - r;  
443 <                if (y<0) y = -y;
444 <                x += y; /* x holds distance */
445 <                        /* >> netbiasshift not needed if funnyshift used */
446 <                if (x<bestd) {bestd=x; best=i;}
447 <                y = x - ((*q)>>funnyshift);  /* y holds biasd */
448 <                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++;
431 >        bestpos = -1;
432 >        bestbiaspos = bestpos;
433 >        p = bias;
434 >        f = freq;
435 >
436 >        for (i=0; i<netsize; i++) {
437 >                n = network[i];
438 >                dist = n[0] - b;   if (dist<0) dist = -dist;
439 >                a = n[1] - g;   if (a<0) a = -a;
440 >                dist += a;
441 >                a = n[2] - r;   if (a<0) a = -a;
442 >                dist += a;
443 >                if (dist<bestd) {bestd=dist; bestpos=i;}
444 >                biasdist = dist - ((*p)>>(intbiasshift-netbiasshift));
445 >                if (biasdist<bestbiasd) {bestbiasd=biasdist; bestbiaspos=i;}
446 >                betafreq = (*f >> betashift);
447 >                *f++ -= betafreq;
448 >                *p++ += (betafreq<<gammashift);
449          }
450 <        freq[best] += beta;
451 <        bias[best] -= betagamma;
452 <        return(bestbias);
450 >        freq[bestpos] += beta;
451 >        bias[bestpos] -= betagamma;
452 >        return(bestbiaspos);
453   }
454  
455  
456 < static
457 < alterneigh(rad,i,b,g,r) /* accepts biased BGR values */
456 > /* move neuron i towards (b,g,r) by factor alpha */
457 >
458 > altersingle(alpha,i,b,g,r)      /* accepts biased BGR values */
459 > register int alpha,i,b,g,r;
460 > {
461 >        register int *n;
462 >
463 >        n = network[i];         /* alter hit neuron */
464 >        *n -= (alpha*(*n - b)) / initalpha;
465 >        n++;
466 >        *n -= (alpha*(*n - g)) / initalpha;
467 >        n++;
468 >        *n -= (alpha*(*n - r)) / initalpha;
469 > }
470 >
471 >
472 > /* move neurons adjacent to i towards (b,g,r) by factor */
473 > /* alpha*(1-((i-j)^2/[r]^2)) precomputed as radpower[|i-j|]*/
474 >
475 > alterneigh(rad,i,b,g,r) /* accents biased BGR values */
476   int rad,i;
477   register int b,g,r;
478   {
479          register int j,k,lo,hi,a;
480          register int *p, *q;
481  
482 <        lo = i-rad;
483 <        if (lo<-1) lo= -1;
468 <        hi = i+rad;
469 <        if (hi>clrtabsiz) hi=clrtabsiz;
482 >        lo = i-rad;   if (lo<-1) lo= -1;
483 >        hi = i+rad;   if (hi>netsize) hi=netsize;
484  
485          j = i+1;
486          k = i-1;
# Line 495 | Line 509 | register int b,g,r;
509   }
510  
511  
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
512   learn()
513   {
514          register int i,j,b,g,r;
515 <        int radius,rad,alpha,step,delta,upto;
515 >        int radius,rad,alpha,step,delta,samplepixels;
516          register unsigned char *p;
517          unsigned char *lim;
518  
519 <        upto = lengthcount/(3*samplefac);
522 <        delta = upto/ncycles;
523 <        lim = thepicture + lengthcount;
519 >        alphadec = 30 + ((samplefac-1)/3);
520          p = thepicture;
521 +        lim = thepicture + lengthcount;
522 +        samplepixels = lengthcount/(3*samplefac);
523 +        delta = samplepixels/ncycles;
524          alpha = initalpha;
525          radius = initradius;
526 +        
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;
531 >        
532 >        if ((lengthcount%prime1) != 0) step = 3*prime1;
533          else {
534 <                if ((lengthcount%jump2) !=0) step = 3*jump2;
534 >                if ((lengthcount%prime2) !=0) step = 3*prime2;
535                  else {
536 <                        if ((lengthcount%jump3) !=0) step = 3*jump3;
537 <                        else step = 3*jump4;
536 >                        if ((lengthcount%prime3) !=0) step = 3*prime3;
537 >                        else step = 3*prime4;
538                  }
539          }
540 +        
541          i = 0;
542 <        while (i < upto) {
542 >        while (i < samplepixels) {
543                  b = p[0] << netbiasshift;
544                  g = p[1] << netbiasshift;
545                  r = p[2] << netbiasshift;
546                  j = contest(b,g,r);
547  
548                  altersingle(alpha,j,b,g,r);
549 <                if (rad) alterneigh(rad,j,b,g,r);
549 <                                                /* alter neighbours */
549 >                if (rad) alterneigh(rad,j,b,g,r);   /* alter neighbours */
550  
551                  p += step;
552                  if (p >= lim) p -= lengthcount;
# Line 563 | Line 563 | learn()
563          }
564   }
565          
566 < static
566 > /* unbias network to give 0..255 entries */
567 > /* which can then be used for colour map */
568 > /* and record position i to prepare for sort */
569 >
570   unbiasnet()
571   {
572          int i,j;
573  
574 <        for (i=0; i<clrtabsiz; i++) {
574 >        for (i=0; i<netsize; i++) {
575                  for (j=0; j<3; j++)
576                          network[i][j] >>= netbiasshift;
577                  network[i][3] = i; /* record colour no */
578          }
579   }
580  
581 < /* Don't do this until the network has been unbiased */
581 >
582 > /* Don't do this until the network has been unbiased (GW) */
583                  
584   static
585   cpyclrtab()
586   {
587          register int i,j,k;
588          
589 <        for (j=0; j<clrtabsiz; j++) {
589 >        for (j=0; j<netsize; j++) {
590                  k = network[j][3];
591                  for (i = 0; i < 3; i++)
592                          clrtab[k][i] = network[j][2-i];

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines