/* 
 * FIG : Facility for Interactive Generation of figures 
 * Copyright (c) 1994 Anthony Dekker 
 * 
 * [NEUQUANT Neural-Net quantization algorithm by Anthony Dekker, 1994. 
 * See "Kohonen neural networks for optimal colour quantization" 
 * in "Network: Computation in Neural Systems" Vol. 5 (1994) pp 351-367. 
 * for a discussion of the algorithm.] 
 * 
 * Any party obtaining a copy of these files is granted, free of charge, a 
 * full and unrestricted irrevocable, world-wide, paid up, royalty-free, 
 * nonexclusive right and license to deal in this software and 
 * documentation files (the "Software"), including without limitation the 
 * rights to use, copy, modify, merge, publish, distribute, sublicense, 
 * and/or sell copies of the Software, and to permit persons who receive 
 * copies from any such party to do so, with the only requirement being 
 * that this copyright notice remain intact. 
 * 
 */ 
 
 
/* 
 * Neural-Net quantization algorithm based on work of Anthony Dekker 
 */ 
 
#include  
#include "fig.h" 
#include "f_neuclrtab.h" 
#include "f_util.h" 
 
static	void initnet(); 
static	void learn(); 
static	void unbiasnet(); 
static	void cpyclrtab(); 
static	void inxbuild(); 
static	int  contest(); 
static	void alterneigh(); 
static	void altersingle(); 
static	int  inxsearch(); 
 
#define MAXNETSIZE	256 
 
/* our color table (global) */

 
BYTE	clrtab[256][3]; 


 
static int	netsize; 
 
#ifndef DEFSMPFAC 
#ifdef SPEED 
#define DEFSMPFAC	(240/SPEED+3) 
#else 
#define DEFSMPFAC	30	/* only sample every 30th pixel */ 
#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) 
 
/* npixels is the total number of pixels defined */ 
 
int

 
neu_init(npixels)		/* initialize our sample array */ 
long	npixels; 
{ 
	if (npixels < MIN_NEU_SAMPLES) 
		samplefac = 1; 
	else 
		samplefac = DEFSMPFAC; 
	return neu_init2(npixels); 
} 
 
int

 
neu_init2(npixels)		/* initialize our sample array */ 
long	npixels; 
{ 
	register int	nsleft; 
	register long	sv; 
	double	rval, cumprob; 
	long	npleft; 
 
	nsamples = npixels/samplefac; 
	if (nsamples < MIN_NEU_SAMPLES) 
		return((int)-MIN_NEU_SAMPLES/nsamples-1); /* THIS IS DIFFERENT FROM THE ORIGINAL -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[N_BLU]; 
		cursamp[1] = col[N_GRN]; 
		cursamp[2] = col[N_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][N_BLU]; 
		cursamp[1] = cs[0][N_GRN]; 
		cursamp[2] = cs[0][N_RED]; 
		cs++; 
		cursamp += 3; 
	} 
	skipcount -= n; 
} 
 


 
neu_clrtab(ncolors)		/* make new color table using ncolors */ 
int	ncolors; 
{ 
	netsize = ncolors; 
	if (netsize > MAXNETSIZE) netsize = MAXNETSIZE; 
	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(netsize); 
} 
 
 
int

 
neu_map_pixel(col)		/* get pixel for color */ 
register BYTE	col[]; 
{ 
	return(inxsearch(col[N_BLU],col[N_GRN],col[N_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][N_BLU],cs[0][N_GRN],cs[0][N_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; 
			neu_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[N_BLU],err[N_GRN],err[N_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 slightly from the original */ 
 
#define bool    	int 
#define false   	0 
#define true    	1 
 
/* network defs */ 
#define maxnetpos	(netsize-1) 
#define netbiasshift	4			/* bias for colour values */ 
#define ncycles		100			/* no. of learning cycles */ 
 
/* defs for freq and bias */ 
#define intbiasshift    16			/* bias for fractions */ 
#define intbias		(((int) 1)<>betashift)	/* beta = 1/1024 */ 
#define betagamma	(intbias<<(gammashift-betashift)) 
 
/* defs for decreasing radius factor */ 
#define initrad		(netsize>>3)		/* for 256 cols, radius starts */ 
#define radiusbiasshift	6			/* at 32.0 biased by 6 bits */ 
#define radiusbias	(((int) 1)< 
int alphadec;					/* biased by 10 bits */ 
 
/* radbias and alpharadbias used for radpower calculation */ 
#define radbiasshift	8 
#define radbias		(((int) 1)< 
pixel network[MAXNETSIZE]; 


 
int netindex[256];	/* for network lookup - really 256 */ 


 
int bias [MAXNETSIZE];	/* bias and freq arrays for learning */

 
int freq [MAXNETSIZE];

 
int radpower[MAXNETSIZE>>3];	/* radpower for precomputation */ 
 
 
/* initialise network in range (0,0,0) to (255,255,255) */ 
 
static void

 
initnet()	 
{ 
	register int i; 
	register int *p; 
	 
	for (i=0; i 
inxbuild() 
{ 
	register int i,j,smallpos,smallval; 
	register int *p,*q; 
	int previouscol,startpos; 
 
	previouscol = 0; 
	startpos = 0; 
	for (i=0; i>1; 
			for (j=previouscol+1; j>1; 
	for (j=previouscol+1; j<256; j++) netindex[j] = maxnetpos; /* really 256 */ 
} 
 
static int

 
inxsearch(b,g,r)  /* accepts real BGR values after net is unbiased */ 
register int b,g,r; 
{ 
	register int i,j,dist,a,bestd; 
	register int *p; 
	int best; 
 
	bestd = 1000;	/* biggest possible dist is 256*3 */ 
	best = -1; 
	i = netindex[g]; /* index on g */ 
	j = i-1;	 /* start at netindex[g] and work outwards */ 
 
	while ((i=0)) { 
		if (i= bestd) i = netsize; /* stop iter */ 
			else { 
				i++; 
				if (dist<0) dist = -dist; 
				a = p[0] - b;   if (a<0) a = -a; 
				dist += a; 
				if (dist=0) { 
			p = network[j]; 
			dist = g - p[1]; /* inx key - reverse dif */ 
			if (dist >= bestd) j = -1; /* stop iter */ 
			else { 
				j--; 
				if (dist<0) dist = -dist; 
				a = p[0] - b;   if (a<0) a = -a; 
				dist += a; 
				if (dist 
contest(b,g,r)	/* accepts biased BGR values */ 
register int b,g,r; 
{ 
	register int i,dist,a,biasdist,betafreq; 
	int bestpos,bestbiaspos,bestd,bestbiasd; 
	register int *p,*f, *n; 
 
	bestd = INT_MAX; 
	bestbiasd = bestd; 
	bestpos = -1; 
	bestbiaspos = bestpos; 
	p = bias; 
	f = freq; 
 
	for (i=0; i>(intbiasshift-netbiasshift)); 
		if (biasdist> betashift); 
		*f++ -= betafreq; 
		*p++ += (betafreq< 
altersingle(alpha,i,b,g,r)	/* accepts biased BGR values */ 
register int alpha,i,b,g,r; 
{ 
	register int *n; 
 
	n = network[i];		/* alter hit neuron */ 
	*n -= (alpha*(*n - b)) / initalpha; 
	n++; 
	*n -= (alpha*(*n - g)) / initalpha; 
	n++; 
	*n -= (alpha*(*n - r)) / initalpha; 
} 
 
 
/* move neurons adjacent to i towards (b,g,r) by factor */ 
/* alpha*(1-((i-j)^2/[r]^2)) precomputed as radpower[|i-j|]*/ 
 
static void

 
alterneigh(rad,i,b,g,r)	/* accents 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>netsize) hi=netsize; 
 
	j = i+1; 
	k = i-1; 
	q = radpower; 
	while ((jlo)) { 
		a = (*(++q)); 
		if (jlo) { 
			p = network[k]; 
			*p -= (a*(*p - b)) / alpharadbias; 
			p++; 
			*p -= (a*(*p - g)) / alpharadbias; 
			p++; 
			*p -= (a*(*p - r)) / alpharadbias; 
			k--; 
		} 
	} 
} 
 
 
static void

 
learn() 
{ 
	register int i,j,b,g,r; 
	int radius,rad,alpha,step,delta,samplepixels; 
	register unsigned char *p; 
	unsigned char *lim; 
 
	alphadec = 30 + ((samplefac-1)/3); 
	p = thepicture; 
	lim = thepicture + lengthcount; 
	samplepixels = lengthcount/(3*samplefac); 
	delta = samplepixels/ncycles; 
	alpha = initalpha; 
	radius = initradius; 
	 
	rad = radius >> radiusbiasshift; 
	if (rad <= 1) rad = 0; 
	for (i=0; i= 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 
unbiasnet() 
{ 
	int i,j; 
 
	for (i=0; i>= netbiasshift; 
		network[i][3] = i; /* record colour no */ 
	} 
} 
 
/* Don't do this until the network has been unbiased */ 
		 
static void

 
cpyclrtab() 
{ 
	register int i,j,k; 
	 
	for (j=0; j