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.6 by greg, Mon Nov 21 15:56:20 1994 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines