src/px/neuclrtab.c

/* Copyright (c) 1994 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 * Neural-Net quantization algorithm based on work of Anthony Dekker
 */

#include "standard.h"

#include "color.h"

#include "random.h"

#ifdef COMPAT_MODE
#define neu_init        new_histo
#define neu_pixel       cnt_pixel
#define neu_colrs       cnt_colrs
#define neu_clrtab      new_clrtab
#define neu_map_pixel   map_pixel
#define neu_map_colrs   map_colrs
#define neu_dith_colrs  dith_colrs
#endif
                                /* our color table (global) */
extern BYTE     clrtab[256][3];
static int      clrtabsiz;

#ifndef DEFSMPFAC
#ifdef SPEED
#define DEFSMPFAC       (240/SPEED+3)
#else
#define DEFSMPFAC       30
#endif
#endif

int     samplefac = DEFSMPFAC;  /* sampling factor */

                /* Samples array starts off holding spacing between adjacent
                 * samples, and ends up holding actual BGR sample values.
                 */
static BYTE     *thesamples;
static int      nsamples;
static BYTE     *cursamp;
static long     skipcount;

#define MAXSKIP         (1<<24-1)

#define nskip(sp)       ((long)(sp)[0]<<16|(long)(sp)[1]<<8|(long)(sp)[2])

#define setskip(sp,n)   ((sp)[0]=(n)>>16,(sp)[1]=((n)>>8)&255,(sp)[2]=(n)&255)


neu_init(npixels)               /* initialize our sample array */
long    npixels;
{
        register int    nsleft;
        register long   sv;
        double  rval, cumprob;
        long    npleft;

        nsamples = npixels/samplefac;
        if (nsamples < 600)
                return(-1);
        thesamples = (BYTE *)malloc(nsamples*3);
        if (thesamples == NULL)
                return(-1);
        cursamp = thesamples;
        npleft = npixels;
        nsleft = nsamples;
        while (nsleft) {
                rval = frandom();       /* random distance to next sample */
                sv = 0;
                cumprob = 0.;
                while ((cumprob += (1.-cumprob)*nsleft/(npleft-sv)) < rval)
                        sv++;
                if (nsleft == nsamples)
                        skipcount = sv;
                else {
                        setskip(cursamp, sv);
                        cursamp += 3;
                }
                npleft -= sv+1;
                nsleft--;
        }
        setskip(cursamp, npleft);       /* tag on end to skip the rest */
        cursamp = thesamples;
        return(0);
}


neu_pixel(col)                  /* add pixel to our samples */
register BYTE   col[];
{
        if (!skipcount--) {
                skipcount = nskip(cursamp);
                cursamp[0] = col[BLU];
                cursamp[1] = col[GRN];
                cursamp[2] = col[RED];
                cursamp += 3;
        }
}


neu_colrs(cs, n)                /* add a scanline to our samples */
register COLR   *cs;
register int    n;
{
        while (n > skipcount) {
                cs += skipcount;
                n -= skipcount+1;
                skipcount = nskip(cursamp);
                cursamp[0] = cs[0][BLU];
                cursamp[1] = cs[0][GRN];
                cursamp[2] = cs[0][RED];
                cs++;
                cursamp += 3;
        }
        skipcount -= n;
}


neu_clrtab(ncolors)             /* make new color table using ncolors */
int     ncolors;
{
        clrtabsiz = ncolors;
        if (clrtabsiz > 256) clrtabsiz = 256;
        initnet();
        learn();
        unbiasnet();
        cpyclrtab();
        inxbuild();
                                /* we're done with our samples */
        free((char *)thesamples);
                                /* reset dithering function */
        neu_dith_colrs((BYTE *)NULL, (COLR *)NULL, 0);
                                /* return new color table size */
        return(clrtabsiz);
}


int
neu_map_pixel(col)              /* get pixel for color */
register BYTE   col[];
{
        return(inxsearch(col[BLU],col[GRN],col[RED]));
}


neu_map_colrs(bs, cs, n)        /* convert a scanline to color index values */
register BYTE   *bs;
register COLR   *cs;
register int    n;
{
        while (n-- > 0) {
                *bs++ = inxsearch(cs[0][BLU],cs[0][GRN],cs[0][RED]);
                cs++;
        }
}


neu_dith_colrs(bs, cs, n)       /* convert scanline to dithered index values */
register BYTE   *bs;
register COLR   *cs;
int     n;
{
        static short    (*cerr)[3] = NULL;
        static int      N = 0;
        int     err[3], errp[3];
        register int    x, i;

        if (n != N) {           /* get error propogation array */
                if (N) {
                        free((char *)cerr);
                        cerr = NULL;
                }
                if (n)
                        cerr = (short (*)[3])malloc(3*n*sizeof(short));
                if (cerr == NULL) {
                        N = 0;
                        map_colrs(bs, cs, n);
                        return;
                }
                N = n;
                bzero((char *)cerr, 3*N*sizeof(short));
        }
        err[0] = err[1] = err[2] = 0;
        for (x = 0; x < n; x++) {
                for (i = 0; i < 3; i++) {       /* dither value */
                        errp[i] = err[i];
                        err[i] += cerr[x][i];
#ifdef MAXERR
                        if (err[i] > MAXERR) err[i] = MAXERR;
                        else if (err[i] < -MAXERR) err[i] = -MAXERR;
#endif
                        err[i] += cs[x][i];
                        if (err[i] < 0) err[i] = 0;
                        else if (err[i] > 255) err[i] = 255;
                }
                bs[x] = inxsearch(err[BLU],err[GRN],err[RED]);
                for (i = 0; i < 3; i++) {       /* propagate error */
                        err[i] -= clrtab[bs[x]][i];
                        err[i] /= 3;
                        cerr[x][i] = err[i] + errp[i];
                }
        }
}

/* The following was adapted and modified from the original (GW)        */
/*----------------------------------------------------------------------*/
/*                                                                      */
/*                              NeuQuant                                */
/*                              --------                                */
/*                                                                      */
/*              Copyright: Anthony Dekker, June 1994                    */
/*                                                                      */
/* This program performs colour quantization of graphics images (SUN    */
/* raster files).  It uses a Kohonen Neural Network.  It produces       */
/* better results than existing methods and runs faster, using minimal  */
/* space (8kB plus the image itself).  The algorithm is described in    */
/* the paper "Kohonen Neural Networks for Optimal Colour Quantization"  */
/* to appear in the journal "Network: Computation in Neural Systems".   */
/* It is a significant improvement of an earlier algorithm.             */
/*                                                                      */
/* This program is distributed free for academic use or for evaluation  */
/* by commercial organizations.                                         */
/*                                                                      */
/*      Usage:  NeuQuant -n inputfile > outputfile                      */
/*                                                                      */
/* where n is a sampling factor for neural learning.                    */
/*                                                                      */
/* Program performance compared with other methods is as follows:       */
/*                                                                      */
/*      Algorithm               |  Av. CPU Time |  Quantization Error   */
/*      -------------------------------------------------------------   */
/*      NeuQuant -3             |  314          |  5.55                 */
/*      NeuQuant -10            |  119          |  5.97                 */
/*      NeuQuant -30            |  65           |  6.53                 */
/*      Oct-Trees               |  141          |  8.96                 */
/*      Median Cut (XV -best)   |  420          |  9.28                 */
/*      Median Cut (XV -slow)   |  72           |  12.15                */
/*                                                                      */
/* Author's address:    Dept of ISCS, National University of Singapore  */
/*                      Kent Ridge, Singapore 0511                      */
/* Email:       [email protected]                                     */
/*----------------------------------------------------------------------*/

#define bool    int
#define false   0
#define true    1

#define initrad                 32
#define radiusdec               30
#define alphadec;               30

/* defs for freq and bias */
#define gammashift      10
#define betashift       gammashift
#define intbiasshift    16
#define intbias         (1<<intbiasshift)
#define gamma           (1<<gammashift)
#define beta            (intbias>>betashift)
#define betagamma       (intbias<<(gammashift-betashift))
#define gammaphi        (intbias<<(gammashift-8))

/* defs for rad and alpha */
#define maxrad          (initrad+1)
#define radiusbiasshift 6
#define radiusbias      (1<<radiusbiasshift)
#define initradius      ((int) (initrad*radiusbias))
#define alphabiasshift  10
#define initalpha       (1<<alphabiasshift)
#define radbiasshift    8
#define radbias         (1<<radbiasshift)
#define alpharadbshift  (alphabiasshift+radbiasshift)
#define alpharadbias    (1<<alpharadbshift)

/* other defs */
#define netbiasshift    4
#define funnyshift      (intbiasshift-netbiasshift)
#define maxnetval       ((256<<netbiasshift)-1)
#define ncycles         100
#define jump1           499     /* prime */
#define jump2           491     /* prime */
#define jump3           487     /* any pic whose size was divisible by all */
#define jump4           503     /* four primes would be simply enormous */

/* cheater definitions (GW) */
#define thepicture      thesamples
#define lengthcount     (nsamples*3)
#define samplefac       1

typedef int pixel[4];  /* BGRc */

static pixel network[256];

static int netindex[256];

static int bias [256];
static int freq [256];
static int radpower[256];       /* actually need only go up to maxrad */

/* fixed space overhead 256*4+256+256+256+256 words = 256*8 = 8kB */


static
initnet()
{
        register int i;
        register int *p;
        
        for (i=0; i<clrtabsiz; i++) {
                p = network[i];
                p[0] = i << netbiasshift;
                p[1] = i << netbiasshift;
                p[2] = i << netbiasshift;
                freq[i] = intbias >> 8;  /* 1/256 */
                bias[i] = 0;
        }
}


static
inxbuild()
{
        register int i,j,smallpos,smallval;
        register int *p,*q;
        int start,previous;

        previous = 0;
        start = 0;
        for (i=0; i<clrtabsiz; i++) {
                p = network[i];
                smallpos = i;
                smallval = p[1];        /* index on g */
                /* find smallest in i+1..clrtabsiz-1 */
                for (j=i+1; j<clrtabsiz; j++) {
                        q = network[j];
                        if (q[1] < smallval) {  /* index on g */
                                smallpos = j;
                                smallval = q[1]; /* index on g */
                        }
                }
                q = network[smallpos];
                if (i != smallpos) {
                        j = q[0];   q[0] = p[0];   p[0] = j;
                        j = q[1];   q[1] = p[1];   p[1] = j;
                        j = q[2];   q[2] = p[2];   p[2] = j;
                        j = q[3];   q[3] = p[3];   p[3] = j;
                }
                /* smallval entry is now in position i */
                if (smallval != previous) {
                        netindex[previous] = (start+i)>>1;
                        for (j=previous+1; j<smallval; j++) netindex[j] = i;
                        previous = smallval;
                        start = i;
                }
        }
        netindex[previous] = (start+clrtabsiz-1)>>1;
        for (j=previous+1; j<clrtabsiz; j++) netindex[j] = clrtabsiz-1;
}


static int
inxsearch(b,g,r)  /* accepts real BGR values after net is unbiased */
register int b,g,r;
{
        register int i,j,best,x,y,bestd;
        register int *p;

        bestd = 1000;   /* biggest possible dist is 256*3 */
        best = -1;
        i = netindex[g]; /* index on g */
        j = i-1;

        while ((i<clrtabsiz) || (j>=0)) {
                if (i<clrtabsiz) {
                        p = network[i];
                        x = p[1] - g;   /* inx key */
                        if (x >= bestd) i = clrtabsiz; /* stop iter */
                        else {
                                i++;
                                if (x<0) x = -x;
                                y = p[0] - b;
                                if (y<0) y = -y;
                                x += y;
                                if (x<bestd) {
                                        y = p[2] - r;   
                                        if (y<0) y = -y;
                                        x += y; /* x holds distance */
                                        if (x<bestd) {bestd=x; best=p[3];}
                                }
                        }
                }
                if (j>=0) {
                        p = network[j];
                        x = g - p[1]; /* inx key - reverse dif */
                        if (x >= bestd) j = -1; /* stop iter */
                        else {
                                j--;
                                if (x<0) x = -x;
                                y = p[0] - b;
                                if (y<0) y = -y;
                                x += y;
                                if (x<bestd) {
                                        y = p[2] - r;   
                                        if (y<0) y = -y;
                                        x += y; /* x holds distance */
                                        if (x<bestd) {bestd=x; best=p[3];}
                                }
                        }
                }
        }
        return(best);
}


static int
contest(b,g,r)  /* accepts biased BGR values */
register int b,g,r;
{
        register int i,best,bestbias,x,y,bestd,bestbiasd;
        register int *p,*q, *pp;

        bestd = ~(1<<31);
        bestbiasd = bestd;
        best = -1;
        bestbias = best;
        q = bias;
        p = freq;
        for (i=0; i<clrtabsiz; i++) {
                pp = network[i];
                x = pp[0] - b;
                if (x<0) x = -x;
                y = pp[1] - g;
                if (y<0) y = -y;
                x += y;
                y = pp[2] - r;  
                if (y<0) y = -y;
                x += y; /* x holds distance */
                        /* >> netbiasshift not needed if funnyshift used */
                if (x<bestd) {bestd=x; best=i;}
                y = x - ((*q)>>funnyshift);  /* y holds biasd */
                if (y<bestbiasd) {bestbiasd=y; bestbias=i;}
                y = (*p >> betashift);       /* y holds beta*freq */
                *p -= y;
                *q += (y<<gammashift);
                p++;
                q++;
        }
        freq[best] += beta;
        bias[best] -= betagamma;
        return(bestbias);
}


static
alterneigh(rad,i,b,g,r) /* accepts biased BGR values */
int rad,i;
register int b,g,r;
{
        register int j,k,lo,hi,a;
        register int *p, *q;

        lo = i-rad;
        if (lo<-1) lo= -1;
        hi = i+rad;
        if (hi>clrtabsiz) hi=clrtabsiz;

        j = i+1;
        k = i-1;
        q = radpower;
        while ((j<hi) || (k>lo)) {
                a = (*(++q));
                if (j<hi) {
                        p = network[j];
                        *p -= (a*(*p - b)) / alpharadbias;
                        p++;
                        *p -= (a*(*p - g)) / alpharadbias;
                        p++;
                        *p -= (a*(*p - r)) / alpharadbias;
                        j++;
                }
                if (k>lo) {
                        p = network[k];
                        *p -= (a*(*p - b)) / alpharadbias;
                        p++;
                        *p -= (a*(*p - g)) / alpharadbias;
                        p++;
                        *p -= (a*(*p - r)) / alpharadbias;
                        k--;
                }
        }
}


static
altersingle(alpha,j,b,g,r)      /* accepts biased BGR values */
register int alpha,j,b,g,r;
{
        register int *q;

        q = network[j];         /* alter hit neuron */
        *q -= (alpha*(*q - b)) / initalpha;
        q++;
        *q -= (alpha*(*q - g)) / initalpha;
        q++;
        *q -= (alpha*(*q - r)) / initalpha;
}


static
learn()
{
        register int i,j,b,g,r;
        int radius,rad,alpha,step,delta,upto;
        register unsigned char *p;
        unsigned char *lim;

        upto = lengthcount/(3*samplefac);
        delta = upto/ncycles;
        lim = thepicture + lengthcount;
        p = thepicture;
        alpha = initalpha;
        radius = initradius;
        rad = radius >> radiusbiasshift;
        if (rad <= 1) rad = 0;
        for (i=0; i<rad; i++) 
                radpower[i] = alpha*(((rad*rad - i*i)*radbias)/(rad*rad));

        if ((lengthcount%jump1) != 0) step = 3*jump1;
        else {
                if ((lengthcount%jump2) !=0) step = 3*jump2;
                else {
                        if ((lengthcount%jump3) !=0) step = 3*jump3;
                        else step = 3*jump4;
                }
        }
        i = 0;
        while (i < upto) {
                b = p[0] << netbiasshift;
                g = p[1] << netbiasshift;
                r = p[2] << netbiasshift;
                j = contest(b,g,r);

                altersingle(alpha,j,b,g,r);
                if (rad) alterneigh(rad,j,b,g,r);
                                                /* alter neighbours */

                p += step;
                if (p >= lim) p -= lengthcount;
        
                i++;
                if (i%delta == 0) {     
                        alpha -= alpha / alphadec;
                        radius -= radius / radiusdec;
                        rad = radius >> radiusbiasshift;
                        if (rad <= 1) rad = 0;
                        for (j=0; j<rad; j++) 
                                radpower[j] = alpha*(((rad*rad - j*j)*radbias)/(rad*rad));
                }
        }
}
        
static
unbiasnet()
{
        int i,j;

        for (i=0; i<clrtabsiz; i++) {
                for (j=0; j<3; j++)
                        network[i][j] >>= netbiasshift;
                network[i][3] = i; /* record colour no */
        }
}

/* Don't do this until the network has been unbiased */
                
static
cpyclrtab()
{
        register int i,j,k;
        
        for (j=0; j<clrtabsiz; j++) {
                k = network[j][3];
                for (i = 0; i < 3; i++)
                        clrtab[k][i] = network[j][2-i];
        }
}
Revision:	2.2
Committed:	Mon Jun 13 13:29:56 1994 UTC (31 years, 4 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.1:	+12 -9 lines
Log Message:	minor bug fixes
#	User	Rev	Content
1	greg	2.1	/* Copyright (c) 1994 Regents of the University of California */
2
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8			* Neural-Net quantization algorithm based on work of Anthony Dekker
9			*/
10
11			#include "standard.h"
12
13			#include "color.h"
14
15			#include "random.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 BYTE clrtab[256][3];
28			static int clrtabsiz;
29
30			#ifndef DEFSMPFAC
31			#ifdef SPEED
32			#define DEFSMPFAC (240/SPEED+3)
33			#else
34			#define DEFSMPFAC 30
35			#endif
36			#endif
37
38			int samplefac = DEFSMPFAC; /* sampling factor */
39
40			/* Samples array starts off holding spacing between adjacent
41			* samples, and ends up holding actual BGR sample values.
42			*/
43			static BYTE *thesamples;
44			static int nsamples;
45			static BYTE *cursamp;
46			static long skipcount;
47
48			#define MAXSKIP (1<<24-1)
49
50			#define nskip(sp) ((long)(sp)[0]<<16\|(long)(sp)[1]<<8\|(long)(sp)[2])
51
52			#define setskip(sp,n) ((sp)[0]=(n)>>16,(sp)[1]=((n)>>8)&255,(sp)[2]=(n)&255)
53
54
55			neu_init(npixels) /* initialize our sample array */
56			long npixels;
57			{
58			register int nsleft;
59			register long sv;
60			double rval, cumprob;
61			long npleft;
62
63			nsamples = npixels/samplefac;
64			if (nsamples < 600)
65			return(-1);
66	greg	2.2	thesamples = (BYTE )malloc(nsamples3);
67	greg	2.1	if (thesamples == NULL)
68			return(-1);
69			cursamp = thesamples;
70			npleft = npixels;
71			nsleft = nsamples;
72			while (nsleft) {
73			rval = frandom(); /* random distance to next sample */
74			sv = 0;
75			cumprob = 0.;
76			while ((cumprob += (1.-cumprob)*nsleft/(npleft-sv)) < rval)
77			sv++;
78	greg	2.2	if (nsleft == nsamples)
79			skipcount = sv;
80			else {
81			setskip(cursamp, sv);
82			cursamp += 3;
83			}
84			npleft -= sv+1;
85	greg	2.1	nsleft--;
86			}
87	greg	2.2	setskip(cursamp, npleft); /* tag on end to skip the rest */
88	greg	2.1	cursamp = thesamples;
89			return(0);
90			}
91
92
93			neu_pixel(col) /* add pixel to our samples */
94			register BYTE col[];
95			{
96			if (!skipcount--) {
97	greg	2.2	skipcount = nskip(cursamp);
98	greg	2.1	cursamp[0] = col[BLU];
99			cursamp[1] = col[GRN];
100			cursamp[2] = col[RED];
101			cursamp += 3;
102			}
103			}
104
105
106			neu_colrs(cs, n) /* add a scanline to our samples */
107			register COLR *cs;
108			register int n;
109			{
110			while (n > skipcount) {
111			cs += skipcount;
112	greg	2.2	n -= skipcount+1;
113			skipcount = nskip(cursamp);
114	greg	2.1	cursamp[0] = cs[0][BLU];
115			cursamp[1] = cs[0][GRN];
116			cursamp[2] = cs[0][RED];
117			cs++;
118			cursamp += 3;
119			}
120			skipcount -= n;
121			}
122
123
124			neu_clrtab(ncolors) /* make new color table using ncolors */
125			int ncolors;
126			{
127			clrtabsiz = ncolors;
128			if (clrtabsiz > 256) clrtabsiz = 256;
129			initnet();
130			learn();
131			unbiasnet();
132			cpyclrtab();
133			inxbuild();
134			/* we're done with our samples */
135			free((char *)thesamples);
136			/* reset dithering function */
137			neu_dith_colrs((BYTE )NULL, (COLR )NULL, 0);
138			/* return new color table size */
139			return(clrtabsiz);
140			}
141
142
143			int
144			neu_map_pixel(col) /* get pixel for color */
145			register BYTE col[];
146			{
147			return(inxsearch(col[BLU],col[GRN],col[RED]));
148			}
149
150
151			neu_map_colrs(bs, cs, n) /* convert a scanline to color index values */
152			register BYTE *bs;
153			register COLR *cs;
154			register int n;
155			{
156			while (n-- > 0) {
157			*bs++ = inxsearch(cs[0][BLU],cs[0][GRN],cs[0][RED]);
158			cs++;
159			}
160			}
161
162
163			neu_dith_colrs(bs, cs, n) /* convert scanline to dithered index values */
164			register BYTE *bs;
165			register COLR *cs;
166			int n;
167			{
168			static short (*cerr)[3] = NULL;
169			static int N = 0;
170			int err[3], errp[3];
171			register int x, i;
172
173			if (n != N) { /* get error propogation array */
174			if (N) {
175			free((char *)cerr);
176			cerr = NULL;
177			}
178			if (n)
179			cerr = (short ()[3])malloc(3n*sizeof(short));
180			if (cerr == NULL) {
181			N = 0;
182			map_colrs(bs, cs, n);
183			return;
184			}
185			N = n;
186			bzero((char )cerr, 3N*sizeof(short));
187			}
188			err[0] = err[1] = err[2] = 0;
189			for (x = 0; x < n; x++) {
190			for (i = 0; i < 3; i++) { /* dither value */
191			errp[i] = err[i];
192			err[i] += cerr[x][i];
193			#ifdef MAXERR
194			if (err[i] > MAXERR) err[i] = MAXERR;
195			else if (err[i] < -MAXERR) err[i] = -MAXERR;
196			#endif
197			err[i] += cs[x][i];
198			if (err[i] < 0) err[i] = 0;
199			else if (err[i] > 255) err[i] = 255;
200			}
201			bs[x] = inxsearch(err[BLU],err[GRN],err[RED]);
202			for (i = 0; i < 3; i++) { /* propagate error */
203			err[i] -= clrtab[bs[x]][i];
204			err[i] /= 3;
205			cerr[x][i] = err[i] + errp[i];
206			}
207			}
208			}
209
210			/* The following was adapted and modified from the original (GW) */
211			/----------------------------------------------------------------------/
212			/* */
213			/* NeuQuant */
214			/* -------- */
215			/* */
216			/* Copyright: Anthony Dekker, June 1994 */
217			/* */
218			/* This program performs colour quantization of graphics images (SUN */
219			/* raster files). It uses a Kohonen Neural Network. It produces */
220			/* better results than existing methods and runs faster, using minimal */
221			/* space (8kB plus the image itself). The algorithm is described in */
222			/* the paper "Kohonen Neural Networks for Optimal Colour Quantization" */
223			/* to appear in the journal "Network: Computation in Neural Systems". */
224			/* It is a significant improvement of an earlier algorithm. */
225			/* */
226			/* This program is distributed free for academic use or for evaluation */
227			/* by commercial organizations. */
228			/* */
229			/* Usage: NeuQuant -n inputfile > outputfile */
230			/* */
231			/* where n is a sampling factor for neural learning. */
232			/* */
233			/* Program performance compared with other methods is as follows: */
234			/* */
235			/* Algorithm \| Av. CPU Time \| Quantization Error */
236			/* ------------------------------------------------------------- */
237			/* NeuQuant -3 \| 314 \| 5.55 */
238			/* NeuQuant -10 \| 119 \| 5.97 */
239			/* NeuQuant -30 \| 65 \| 6.53 */
240			/* Oct-Trees \| 141 \| 8.96 */
241			/* Median Cut (XV -best) \| 420 \| 9.28 */
242			/* Median Cut (XV -slow) \| 72 \| 12.15 */
243			/* */
244			/* Author's address: Dept of ISCS, National University of Singapore */
245			/* Kent Ridge, Singapore 0511 */
246			/* Email: [email protected] */
247			/----------------------------------------------------------------------/
248
249			#define bool int
250			#define false 0
251			#define true 1
252
253			#define initrad 32
254			#define radiusdec 30
255			#define alphadec; 30
256
257			/* defs for freq and bias */
258			#define gammashift 10
259			#define betashift gammashift
260			#define intbiasshift 16
261			#define intbias (1<<intbiasshift)
262			#define gamma (1<<gammashift)
263			#define beta (intbias>>betashift)
264			#define betagamma (intbias<<(gammashift-betashift))
265			#define gammaphi (intbias<<(gammashift-8))
266
267			/* defs for rad and alpha */
268			#define maxrad (initrad+1)
269			#define radiusbiasshift 6
270			#define radiusbias (1<<radiusbiasshift)
271			#define initradius ((int) (initrad*radiusbias))
272			#define alphabiasshift 10
273			#define initalpha (1<<alphabiasshift)
274			#define radbiasshift 8
275			#define radbias (1<<radbiasshift)
276			#define alpharadbshift (alphabiasshift+radbiasshift)
277			#define alpharadbias (1<<alpharadbshift)
278
279			/* other defs */
280			#define netbiasshift 4
281			#define funnyshift (intbiasshift-netbiasshift)
282			#define maxnetval ((256<<netbiasshift)-1)
283			#define ncycles 100
284			#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 */
288
289			/* cheater definitions (GW) */
290			#define thepicture thesamples
291			#define lengthcount (nsamples*3)
292			#define samplefac 1
293
294			typedef int pixel[4]; /* BGRc */
295
296			static pixel network[256];
297
298			static int netindex[256];
299
300			static int bias [256];
301			static int freq [256];
302			static int radpower[256]; /* actually need only go up to maxrad */
303
304			/* fixed space overhead 2564+256+256+256+256 words = 2568 = 8kB */
305
306
307			static
308			initnet()
309			{
310			register int i;
311			register int *p;
312
313			for (i=0; i<clrtabsiz; i++) {
314			p = network[i];
315			p[0] = i << netbiasshift;
316			p[1] = i << netbiasshift;
317			p[2] = i << netbiasshift;
318			freq[i] = intbias >> 8; /* 1/256 */
319			bias[i] = 0;
320			}
321			}
322
323
324			static
325			inxbuild()
326			{
327			register int i,j,smallpos,smallval;
328			register int p,q;
329			int start,previous;
330
331			previous = 0;
332			start = 0;
333			for (i=0; i<clrtabsiz; i++) {
334			p = network[i];
335			smallpos = i;
336			smallval = p[1]; /* index on g */
337			/* find smallest in i+1..clrtabsiz-1 */
338			for (j=i+1; j<clrtabsiz; j++) {
339			q = network[j];
340			if (q[1] < smallval) { /* index on g */
341			smallpos = j;
342			smallval = q[1]; /* index on g */
343			}
344			}
345			q = network[smallpos];
346			if (i != smallpos) {
347			j = q[0]; q[0] = p[0]; p[0] = j;
348			j = q[1]; q[1] = p[1]; p[1] = j;
349			j = q[2]; q[2] = p[2]; p[2] = j;
350			j = q[3]; q[3] = p[3]; p[3] = j;
351			}
352			/* smallval entry is now in position i */
353			if (smallval != previous) {
354			netindex[previous] = (start+i)>>1;
355			for (j=previous+1; j<smallval; j++) netindex[j] = i;
356			previous = smallval;
357			start = i;
358			}
359			}
360			netindex[previous] = (start+clrtabsiz-1)>>1;
361			for (j=previous+1; j<clrtabsiz; j++) netindex[j] = clrtabsiz-1;
362			}
363
364
365			static int
366			inxsearch(b,g,r) /* accepts real BGR values after net is unbiased */
367			register int b,g,r;
368			{
369			register int i,j,best,x,y,bestd;
370			register int *p;
371
372			bestd = 1000; /* biggest possible dist is 2563 /
373			best = -1;
374			i = netindex[g]; /* index on g */
375			j = i-1;
376
377			while ((i<clrtabsiz) \|\| (j>=0)) {
378			if (i<clrtabsiz) {
379			p = network[i];
380			x = p[1] - g; /* inx key */
381			if (x >= bestd) i = clrtabsiz; /* stop iter */
382			else {
383			i++;
384			if (x<0) x = -x;
385			y = p[0] - b;
386			if (y<0) y = -y;
387			x += y;
388			if (x<bestd) {
389			y = p[2] - r;
390			if (y<0) y = -y;
391			x += y; /* x holds distance */
392			if (x<bestd) {bestd=x; 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 */
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];}
411			}
412			}
413			}
414			}
415			return(best);
416			}
417
418
419			static int
420			contest(b,g,r) /* accepts biased BGR values */
421			register int b,g,r;
422			{
423			register int i,best,bestbias,x,y,bestd,bestbiasd;
424			register int p,q, *pp;
425
426			bestd = ~(1<<31);
427			bestbiasd = bestd;
428			best = -1;
429			bestbias = best;
430			q = bias;
431			p = freq;
432			for (i=0; i<clrtabsiz; i++) {
433			pp = network[i];
434			x = pp[0] - b;
435			if (x<0) x = -x;
436			y = pp[1] - g;
437			if (y<0) y = -y;
438			x += y;
439			y = pp[2] - r;
440			if (y<0) y = -y;
441			x += y; /* x holds distance */
442			/* >> netbiasshift not needed if funnyshift used */
443			if (x<bestd) {bestd=x; best=i;}
444			y = x - ((q)>>funnyshift); / y holds biasd */
445			if (y<bestbiasd) {bestbiasd=y; bestbias=i;}
446			y = (p >> betashift); / y holds betafreq /
447			*p -= y;
448			*q += (y<<gammashift);
449			p++;
450			q++;
451			}
452			freq[best] += beta;
453			bias[best] -= betagamma;
454			return(bestbias);
455			}
456
457
458			static
459			alterneigh(rad,i,b,g,r) /* accepts biased BGR values */
460			int rad,i;
461			register int b,g,r;
462			{
463			register int j,k,lo,hi,a;
464			register int p, q;
465
466			lo = i-rad;
467			if (lo<-1) lo= -1;
468			hi = i+rad;
469			if (hi>clrtabsiz) hi=clrtabsiz;
470
471			j = i+1;
472			k = i-1;
473			q = radpower;
474			while ((j<hi) \|\| (k>lo)) {
475			a = (*(++q));
476			if (j<hi) {
477			p = network[j];
478			p -= (a(*p - b)) / alpharadbias;
479			p++;
480			p -= (a(*p - g)) / alpharadbias;
481			p++;
482			p -= (a(*p - r)) / alpharadbias;
483			j++;
484			}
485			if (k>lo) {
486			p = network[k];
487			p -= (a(*p - b)) / alpharadbias;
488			p++;
489			p -= (a(*p - g)) / alpharadbias;
490			p++;
491			p -= (a(*p - r)) / alpharadbias;
492			k--;
493			}
494			}
495			}
496
497
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;
518			register unsigned char *p;
519			unsigned char *lim;
520
521			upto = lengthcount/(3*samplefac);
522			delta = upto/ncycles;
523			lim = thepicture + lengthcount;
524			p = thepicture;
525			alpha = initalpha;
526			radius = initradius;
527			rad = radius >> radiusbiasshift;
528			if (rad <= 1) rad = 0;
529			for (i=0; i<rad; i++)
530			radpower[i] = alpha(((radrad - ii)radbias)/(rad*rad));
531
532			if ((lengthcount%jump1) != 0) step = 3*jump1;
533			else {
534			if ((lengthcount%jump2) !=0) step = 3*jump2;
535			else {
536			if ((lengthcount%jump3) !=0) step = 3*jump3;
537			else step = 3*jump4;
538			}
539			}
540			i = 0;
541			while (i < upto) {
542			b = p[0] << netbiasshift;
543			g = p[1] << netbiasshift;
544			r = p[2] << netbiasshift;
545			j = contest(b,g,r);
546
547			altersingle(alpha,j,b,g,r);
548			if (rad) alterneigh(rad,j,b,g,r);
549			/* alter neighbours */
550
551			p += step;
552			if (p >= lim) p -= lengthcount;
553
554			i++;
555			if (i%delta == 0) {
556			alpha -= alpha / alphadec;
557			radius -= radius / radiusdec;
558			rad = radius >> radiusbiasshift;
559			if (rad <= 1) rad = 0;
560			for (j=0; j<rad; j++)
561			radpower[j] = alpha(((radrad - jj)radbias)/(rad*rad));
562			}
563			}
564			}
565
566			static
567			unbiasnet()
568			{
569			int i,j;
570
571			for (i=0; i<clrtabsiz; i++) {
572			for (j=0; j<3; j++)
573			network[i][j] >>= netbiasshift;
574			network[i][3] = i; /* record colour no */
575			}
576			}
577
578			/* Don't do this until the network has been unbiased */
579
580			static
581			cpyclrtab()
582			{
583			register int i,j,k;
584
585			for (j=0; j<clrtabsiz; j++) {
586			k = network[j][3];
587			for (i = 0; i < 3; i++)
588			clrtab[k][i] = network[j][2-i];
589			}
590			}