| 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 | 
} |