package org.gif.encoder;

import snow.platform.native.utils.ByteArray;

/**
 * This class handles the Neural-Net quantization algorithm.
 * Based on https://code.google.com/p/as3gif/
 * 
 * NeuQuant Neural-Net Quantization Algorithm
 * ------------------------------------------
 * 
 * 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 from the author, directly or
 * indirectly, 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.
 * 
 * @author Steve Richey (Haxe/OpenFL version)
 * @author Thibault Imbert (AS3 version - bytearray.org)
 * @author Kevin Weiner (original Java version - kweiner@fmsware.com)
 * @author Stuart Coyle (Modified to quantize 32bit RGBA images for the pngnq program, modified to accept a number of colors argument.)
 * @author Kornel Lesinski (Euclidean distance, color matching dependent on alpha channel and with gamma correction. code refreshed for modern compilers/architectures: ANSI C, floats, removed pointer tricks and used arrays and structs.)
 * 
 * @version 0.1 Haxe implementation
 */

import haxe.io.Bytes;
import haxe.io.BytesBuffer;

class LearningStatus
{
	public var i:Int = 0;
	public var al:Int = 0;
	public var b:Int = 0;
	public var g:Int = 0;
	public var r:Int = 0;
	public var p:Int = 0;
	public var rad:Int = 0;
	public var step:Int = 0;
	public var delta:Int = 0;
	public var samplepixels:Int = 0;
	public var radius:Float = 0;
	public var alpha:Float = 0;
	public var started:Bool = false;
	public var finished:Bool = false;
	public function new() { }
}

class QuantizationStatus
{
	public var finished:Bool = false;
	public var percentage:Int = 0;
	public var workChunkTimeSize:Float = 0;
	public var chunkStartedTime:Float = 0;
	public function new() {	}
}

typedef NqPixel = {
	var al:Float;
	var b:Float;
	var g:Float;
	var r:Float;
}

typedef NqColormap = {
	var r:Int;
	var g:Int;
	var b:Int;
	var al:Int;
}

typedef NqLabColormap = {
	var l:Float;
	var a:Float;
	var b:Float;
	var al:Int;
}

class NeuQuant
{
	/**
	 * Maximum number of colours that can be used. actual number is now passed to initcolors.
	 */
	private static inline var MAX_NET_SIZE:Int = 256;
	/**
	 * Four primes near 500 - assume no image has a length so large that it is divisible by all four primes
	 */
	private static inline var PRIME_1:Int = 499;
    private static inline var PRIME_2:Int = 491;
    private static inline var PRIME_3:Int = 487;
    private static inline var PRIME_4:Int = 503;
    private static inline var MIN_PICTURE_BYTES:Int = 4 * PRIME_4;
	/**
	 * Network Definitions
	 */
	private static inline var MAX_NET_POS:Int = MAX_NET_SIZE - 1;
	/**
	 * Number of learning cycles
	 */
	private static inline var N_CYCLES:Int = 100;
	/**
	 * Defs for freq and bias
	 */
	private static inline var GAMMA_SHIFT:Int = 10; /* gamma = 1024 */
	private static inline var GAMMA:Float = 1 << GAMMA_SHIFT;
	private static inline var BETA_SHIFT:Int = 10;
	private static inline var BETA:Float = 1.0 / (1 << BETA_SHIFT); /* beta = 1/1024 */
	private static inline var BETA_GAMMA:Float = 1 << (GAMMA_SHIFT - BETA_SHIFT);
	/**
	 * defs for decreasing radius factor
	 */
	private static inline var INIT_RAD:Int = MAX_NET_SIZE >> 3; /* for 256 cols, radius starts */
	private static inline var INIT_RADIUS:Float = INIT_RAD * 1.0; /* and decreases by a */
	private static inline var RADIUS_DEC:Float = 30; /* factor of 1/30 each cycle */ 
	/**
	 * defs for decreasing alpha factor
	 */
	private static inline var ALPHA_BIAS_SHIFT:Int = 10; /* alpha starts at 1.0 */
	private static inline var INIT_ALPHA:Float = 1 << ALPHA_BIAS_SHIFT;
	private static var alphadec:Float = 0; /* biased by 10 bits, non-constant */
	/**
	 * radbias and alpharadbias used for radpower calculation
	 */
	private static inline var RAD_BIAS_SHIFT:Int = 8;
	private static inline var RAD_BIAS:Int = 1 << RAD_BIAS_SHIFT;
	private static inline var ALPHA_RAD_B_SHIFT:Int = ALPHA_BIAS_SHIFT + RAD_BIAS_SHIFT;
	private static inline var ALPHA_RAD_BIAS:Float = 1 << ALPHA_RAD_B_SHIFT;
	/**
	 * The input image itself
	 */
	private var thepicture:Bytes;
	/**
	 * ARGB (e.g. Flash) or RGBA (e.g. PNG)
	 */
	private var isARGB:Bool;
	/**
	 * Lengthcount = H*W*3
	 */
	private var lengthcount:Int;
	/**
	 * The network itself
	 */
	private var network: Array<NqPixel>;
	/**
	 * For network lookup - really 256
	 */
	private var netindex: Array<Int>;
	/**
	 * Bias and freq arrays for learning
	 */
	private var bias:Array<Float>;
	private var freq:Array<Float>;
	/**
	 * Radpower for precomputation
	 */
	private var radpower:Array<Float>;
	/**
	 * Number of colours to use. Made a global instead of #define
	 */
	private var netsize:Int = 0;
	/**
	 * 1.0/2.2 usually
	 */
	private var gamma_correction:Float = 1.0;
	
	private var compareInLabColorSpace:Bool;
	
	/**
	 * Just instantiates this NeuQuant instance.
	 */
	public function new() { }
	
	/**
	 * Quantizes an image. An image should be a byte array with 32 bit per pixel. It supports ARGB and RGBA.
	 * 
	 * @param	image Bytes
	 * @param	isARGB Bool
	 * @param	colorsAmount A result number of colors
	 * @param	gammaCorrection A gamma correction factor
	 * @param	sampleFactor A quality factor in [1..30]. 1 is the best, 30 is the worst.
	 * @param	verbose Bool
	 */
	public function quantize(image:Bytes, isARGB:Bool, colorsAmount:Int, gammaCorrection:Float, sampleFactor:Int, verbose:Bool)
	{
		quantizeAsync(image, isARGB, colorsAmount, gammaCorrection, sampleFactor, verbose, 31536000000 /* a year */);
	}
	
	var quantizationStatus:QuantizationStatus;
	
	/**
	 * It quantizes an image in asynchronous way. Whole work is separated to small chunks with specified maximum operating time. 
	 * When each chunk finishes the method returns current status (QuantizationStatus).
	 * 
	 * @param	image Bytes
	 * @param	isARGB Bool
	 * @param	colorsAmount
	 * @param	colorsAmount A result number of colors
	 * @param	gammaCorrection A gamma correction factor
	 * @param	sampleFactor A quality factor in [1..30]. 1 is the best, 30 is the worst.
	 * @param	verbose Bool
	 * @param	workChunkTimeSize Maximum operating time for a current chunk in milliseconds.
	 * @return  A currect status of quantization (QuantizeStatus)
	 */
	public function quantizeAsync(image:Bytes, isARGB:Bool, colorsAmount:Int, gammaCorrection:Float, sampleFactor:Int, verbose:Bool, workChunkTimeSize: Float, ?compareInLabColorSpace: Bool):QuantizationStatus
	{
		this.isARGB = isARGB;
		this.compareInLabColorSpace = compareInLabColorSpace;
		
		if (quantizationStatus == null || quantizationStatus.finished)
		{
			quantizationStatus = new QuantizationStatus();
			quantizationStatus.workChunkTimeSize = workChunkTimeSize;
			learningStatus = new LearningStatus();
			initnet(image, colorsAmount, gammaCorrection);
		}
		
		quantizationStatus.chunkStartedTime = Date.now().getTime();
		
		learn(sampleFactor, verbose);
		quantizationStatus.percentage = Std.int(learningStatus.i / learningStatus.samplepixels * 100);
		
		if (!learningStatus.finished)
		{
			return quantizationStatus; //Learning was interruped
		}
		
		inxbuild();
		
		if (compareInLabColorSpace)
		{
			calculateLabColors();
		}
		
		quantizationStatus.finished = true;
		
		return quantizationStatus;
	}
	
	/**
	 * Initialise network in range (0,0,0,0) to (255,255,255,255) and set parameters
	 */
	var biasvalues:Array<Float>;
	
	function initnet(thepic: Bytes, colours: Int, gamma_c: Float):Void
	{
		var i: Int;
		
		gamma_correction = gamma_c;
		
		thepicture = thepic;
		lengthcount = thepic.length;
		netsize = colours; 
		
		biasvalues = new Array<Float>();
		
		for (i in 0...MAX_NET_SIZE)
		{
			var temp: Float;
			temp = Math.pow(i/255.0, 1.0/gamma_correction) * 255.0;
			temp = Math.round(temp);
			biasvalues.push(temp);
		}
		
		network = new Array<NqPixel>();
		freq = new Array<Float>();
		bias = new Array<Float>();
		
		for (i in 0 ... netsize)
		{
			var v: Float = biasvalue(Std.int(i * 256 / netsize));
			network.push({b:v, g:v, r:v, al: if (i < 16) (i*16)*1.0 else 255.0});
			freq.push(1.0/netsize);	/* 1/netsize */
			bias.push(0.0);
		}
	}

	function unbiasvalue(temp: Float):Int
	{
		if (temp < 0)
		{
			return 0;
		}
		
		temp = Math.pow(temp/255.0, gamma_correction) * 255.0;    
		temp = Math.floor((temp / 255.0 * 256.0));
		
		if (temp > 255)
		{
			return 255;
		}
		
		return Std.int(temp);
	}

	function round_biased(temp:Float):Int
	{    
		if (temp < 0)
		{
			return 0;
		}
		
		temp = Math.floor((temp / 255.0 * 256.0));
		
		if (temp > 255)
		{
			return 255;    
		}
		
		return Std.int(temp);
	}
	
	inline function biasvalue(temp:Int):Float
	{    
		return biasvalues[temp];
	}

	/* Output colormap to unsigned char ptr in RGBA format */
	//function getcolormap(map: BytesBuffer): Void {
		//var j: Int;
		//for(j in 0 ... netsize) {
			//map.addByte(unbiasvalue(network[j].r));
			//map.addByte(unbiasvalue(network[j].g));
			//map.addByte(unbiasvalue(network[j].b));
			//map.addByte(round_biased(network[j].al));
		//}
	//}
	
	/**
	 * Insertion sort of network and building of netindex[0..255] (to do after unbias)
	 */
	var colormap: Array<NqColormap>;
	var labColorMap: Array<NqLabColormap>;

	function inxbuild(): Void
	{
		var i: Int, j: Int, smallpos: Int, smallval: Int;
		var previouscol: Int, startpos: Int;
		
		netindex = new Array<Int>();
		
		for (i in 0...256)
		{
			netindex.push(0);
		}
		
		colormap = new Array<NqColormap>();
		
		for (i in 0...netsize)
		{
			colormap.push(
				{
					r: Std.int(biasvalue(unbiasvalue(network[i].r))),
					g: Std.int(biasvalue(unbiasvalue(network[i].g))),
					b: Std.int(biasvalue(unbiasvalue(network[i].b))),
					al: round_biased(network[i].al)
				});
		}
		
		previouscol = 0;
		startpos = 0;
		
		for (i in 0 ... netsize)
		{
			smallpos = i;
			smallval = (colormap[i].g); // index on g
			
			// find smallest in i..netsize-1
			
			for (j in (i + 1)...netsize)
			{
				if ((colormap[j].g) < smallval)
				{
					// index on g
					smallpos = j;
					smallval = (colormap[j].g);	// index on g
				}
			}
			
			// swap colormap[i] (i) and colormap[smallpos] (smallpos) entries
			
			if (i != smallpos)
			{
				var temp = network[smallpos];   network[smallpos] = network[i];   network[i] = temp;
				var tempc = colormap[smallpos];   colormap[smallpos] = colormap[i];   colormap[i] = tempc;
			}
			
			// smallval entry is now in position i
			
			if (smallval != previouscol)
			{
				netindex[previouscol] = (startpos+i)>>1;
				
				for (j in (previouscol + 1)...(smallval))
				{
					netindex[j] = i;
				}
				
				previouscol = smallval;
				startpos = i;
			}
		}
		
		netindex[previouscol] = (startpos + MAX_NET_POS) >> 1;
		
		for (j in (previouscol + 1)...256)
		{
			netindex[j] = MAX_NET_POS; /* really 256 */
		}
	}
	
	inline function colorimportance(al:Float):Float
	{
		var transparency: Float = 1.0 - al/255.0;
		return (1.0 - transparency*transparency);
	}

	function calculateLabColors()
	{
		labColorMap = [];
		
		for (i in 0...netsize)
		{
			var xyz = rgb2xyz(colormap[i].r, colormap[i].g, colormap[i].b);
			var lab = xyz2cielab(xyz.x, xyz.y, xyz.z );
			labColorMap.unshift({ al: colormap[i].al, l: lab.l, a: lab.a, b: lab.b });
		}
	}
	
	function getLabDistance(i:Int, al:Int, l:Float, a:Float, b:Float):Float
	{
		var	x: Float = labColorMap[i].l - l;
		var dist: Float = x * x;
		
		x = labColorMap[i].a - a;
		dist += x*x;
		
		x = labColorMap[i].b - b;
		dist += x*x;
		
		x = labColorMap[i].al - al;
		dist += x * x;
		
		return dist;
	}
	
	function getRgbDistance(i:Int, al:Int, b:Int, g:Int, r:Int)
	{
		var	a: Float = colormap[i].r - r;
		var dist: Float = a*a;
		
		a = colormap[i].g - g;
		dist += a*a;
		
		a = colormap[i].b - b;
		dist += a*a;
		
		a = colormap[i].al - al;
		dist += a * a;
		return dist;
		
	}
	
	/**
	 * Search for ABGR values 0..255 (after net is unbiased) and return color index
	 * 
	 * @param	al
	 * @param	b
	 * @param	g
	 * @param	r
	 * @return
	 */
	function slowinxsearch(al: Int, b: Int, g: Int, r: Int): Int
	{
		var i:Int = 0;
		var best:Int = 0;
		var bestd:Float = -1;
		var dist:Float = 0;
		
		r = Std.int(biasvalue(r));
		g = Std.int(biasvalue(g));
		b = Std.int(biasvalue(b));
		
		var lab:Dynamic = null;
		
		if (compareInLabColorSpace)
		{
			var xyz = rgb2xyz(r, g, b);
			lab = xyz2cielab(xyz.x, xyz.y, xyz.z );
		}
		
		for (i in 0...netsize)
		{
			if (compareInLabColorSpace)
			{
				dist = getLabDistance(i, al, lab.l, lab.a, lab.b);
			}
			else
			{
				dist = getRgbDistance(i, al, b, g, r);
			}
			
			if (dist < bestd || bestd == -1)
			{
				bestd = dist; best = i;
			}
		}
		
		return best;
	}

    /**
     * Getting an index of a specified color in a resulted palette.
     * @param	color 32-bit integer
     */
	public function lookup (a:Int,r:Int, g:Int, b:Int)
	{
	    return slowinxsearch(a, b, g, r);
	}

	/**
	 * Getting a 32-bit ARGB integer for a index of color using a current palette.
	 * @param	index Int
	 * @return
	 */
	public function getColor(index:Int):Int
	{
		return (colormap[index].al << 24) | (colormap[index].r << 16) | (colormap[index].g << 8) | colormap[index].b;
	}
	
	/**
	 * Returns the current palette as a ByteArray.
	 */
	public function getColorMap():ByteArray
	{
		var map:ByteArray = new ByteArray(netsize*3);
		var index:Array<Int> = [for (i in 0...netsize) 0];
		
		for (i in 0...netsize)
		{
			index[Std.int(network[i].al * 255)] = i;
		}
		
		var k:Int = 0;
		var j:Int = 0;
		for (i in 0...netsize)
		{
			j = i;// index[i];
			map[k++] = Std.int(network[j].r * 255);
			map[k++] = Std.int(network[j].g * 255);
			map[k++] = Std.int(network[j].b * 255);
		}
		
		return map;
	}
	
	function inxsearch(al:Int, b:Int, g:Int, r:Int):Int
	{
		var i:Int; 
		var j:Int; 
		var dist: Float, a: Float, bestd: Float;
		var best: Int;
		
		bestd = 1<<30;		/* biggest possible dist */
		best = 0;
	 
		if (al > 0)	{       
			r = Std.int(biasvalue(r));
			g = Std.int(biasvalue(g));
			b = Std.int(biasvalue(b));
		} else {
			r=g=b=0;
		}
		
		i = netindex[(g)];	/* index on g */
		j = i-1;		/* start at netindex[g] and work outwards */
		
		var colimp: Float = colorimportance(al);
		
		while ((i<netsize) || (j>=0)) {
			if (i<netsize) {
				a = colormap[i].g - g;		/* inx key */
				dist = a*a * colimp;
				if (dist > bestd) break;	/* stop iter */
				else {								
					a = colormap[i].r - r;
					dist += a*a * colimp;				
					if (dist<bestd) {
						a = colormap[i].b - b;
						dist += a*a * colimp;
						if(dist<bestd) {
							a = colormap[i].al - al;
							dist += a*a;
							if (dist<bestd) {bestd=dist; best=i;}
						}                    
					}								
					i++;
				}            
			}
			if (j>=0) {
				a = colormap[j].g - g; /* inx key - reverse dif */
				dist = a*a * colimp;
				if (dist > bestd) break; /* stop iter */
				else {				
					a = colormap[j].b - b;
					dist += a*a * colimp;				
					if (dist<bestd) {
						a = colormap[j].r - r;
						dist += a*a * colimp;                    
						if(dist<bestd) {
							a = colormap[j].al - al;
							dist += a*a;
							if (dist<bestd) {bestd=dist; best=j;}
						}			                    
					}				
					j--;
				}            
			}
		}
		
		return(best);
	}
	
	/**
	 * Search for biased ABGR values.
	 * 
	 * @param	al
	 * @param	b
	 * @param	g
	 * @param	r
	 * @return
	 */
	 function contest(al: Float, b: Float, g: Float, r: Float): Int {
		/* finds closest neuron (min dist) and updates freq */
		/* finds best neuron (min dist-bias) and returns position */
		/* for frequently chosen neurons, freq[i] is high and bias[i] is negative */
		/* bias[i] = gamma*((1/netsize)-freq[i]) */

		var i: Int; 
		var dist: Float, a: Float, betafreq: Float;
		var bestpos: Int, bestbiaspos: Int;
		var bestd: Float, bestbiasd: Float;
		
		bestd = 1<<30;
		bestbiasd = bestd;
		bestpos = 0;
		bestbiaspos = bestpos;

		var colimp: Float = colorimportance(al);
		
		for (i in 0...netsize) {
			var bestbiasd_biased: Float = bestbiasd + bias[i];
			
			a = network[i].b - b;
			dist = Math.abs(a) * colimp;
			a = network[i].r - r;
			dist += Math.abs(a) * colimp;
			
			if (dist < bestd || dist < bestbiasd_biased)
			{                 
				a = network[i].g - g;
				dist += Math.abs(a) * colimp;
				a = network[i].al - al; 
				dist += Math.abs(a);    
				
				if (dist<bestd) {bestd=dist; bestpos=i;}
				if (dist<bestbiasd_biased) {bestbiasd=dist - bias[i]; bestbiaspos=i;}
			}
			betafreq = freq[i] / (1 << BETA_SHIFT);
			freq[i] -= betafreq;
			bias[i] += betafreq * (1 << GAMMA_SHIFT);
		}
		freq[bestpos] += BETA;
		bias[bestpos] -= BETA_GAMMA;
		
		return bestbiaspos;
	}
	
	/**
	 * Move neuron i towards biased (a,b,g,r) by factor alpha.
	 * 
	 * @param	alpha
	 * @param	i
	 * @param	al
	 * @param	b
	 * @param	g
	 * @param	r
	 */
	function altersingle(alpha:Float, i:Int, al:Float, b:Float, g:Float, r:Float):Void
	{    
		var colorimp: Float = 1.0;//0.5;// + 0.7*colorimportance(al);
		
		alpha /= INIT_ALPHA;
		
		/* alter hit neuron */
		network[i].al -= alpha*(network[i].al - al);
		network[i].b -= colorimp*alpha*(network[i].b - b);
		network[i].g -= colorimp*alpha*(network[i].g - g);
		network[i].r -= colorimp*alpha*(network[i].r - r);
	}
	
	/**
	 * Move adjacent neurons by precomputed alpha*(1-((i-j)^2/[r]^2)) in radpower[|i-j|].
	 * 
	 * @param	rad
	 * @param	i
	 * @param	al
	 * @param	b
	 * @param	g
	 * @param	r
	 */
	function alterneigh(rad:Int, i:Int, al:Float, b:Float, g:Float, r:Float):Void
	{
		var j: Int, hi: Int;
		var k: Int, lo: Int, q: Int;
		var a: Float;
		
		lo = i-rad;   if (lo<0) lo=0;
		hi = i+rad;   if (hi>netsize-1) hi=netsize-1;
		
		j = i+1;
		k = i-1;
		q = 0;
		
		while ((j<=hi) || (k>=lo)) {
			a = (radpower[q++]) / ALPHA_RAD_BIAS;
			if (j<=hi) {
				network[j].al -= a*(network[j].al - al);
				network[j].b  -= a*(network[j].b  - b) ;
				network[j].g  -= a*(network[j].g  - g) ;
				network[j].r  -= a*(network[j].r  - r) ;
				j++;
			}
			if (k>=lo) {
				network[k].al -= a*(network[k].al - al);
				network[k].b  -= a*(network[k].b  - b) ;
				network[k].g  -= a*(network[k].g  - g) ;
				network[k].r  -= a*(network[k].r  - r) ;
				k--;
			}
		}
	}
	
	/**
	 * Main Learning Loop
	 */
	var learningStatus:LearningStatus;
	/**
	 * Sampling factor 1..30
	 * 
	 * @param	samplefac
	 * @param	verbose
	 */
	function learn(samplefac: Int, verbose: Bool)
	{ /* Stu: N.B. added parameter so that main() could control verbosity. */
		var s: LearningStatus = learningStatus;
		
		if (s.finished)
		{
			return; //All work was done
		}
		
		if (!s.started)
		{
			alphadec = 30 + ((samplefac - 1) / 3);
			s.samplepixels = Std.int(lengthcount / (4 * samplefac)); 
			s.delta = Std.int(s.samplepixels / N_CYCLES);  /* here's a problem with small images: samplepixels < ncycles => delta = 0 */
			if(s.delta==0) s.delta = 1;        /* kludge to fix */
			s.alpha = INIT_ALPHA;
			s.radius = INIT_RADIUS;
			
			s.p = 0;
			
			s.rad = Std.int(s.radius);
			if (s.rad <= 1) s.rad = 0;
			radpower = new Array<Float>();
			
			for (i in 0...(s.rad))
			{
				radpower.push(Math.floor(s.alpha * (((s.rad * s.rad - i * i) * RAD_BIAS) / (s.rad * s.rad))));
			}
			
			if (verbose)
			{
				trace("beginning 1D learning: initial radius=" + s.rad + "\n");
			}
			
			if ((lengthcount % PRIME_1) != 0)
			{
				s.step = 4 * PRIME_1;
			}
			else
			{
				if ((lengthcount % PRIME_2) != 0)
				{
					s.step = 4 * PRIME_2;
				}
				else
				{
					if ((lengthcount % PRIME_3) != 0)
					{
						s.step = 4 * PRIME_3;
					}
					else
					{
						s.step = PRIME_4;
					}
				}
			}
			
			s.i = 0;
			s.started = true;
		}
		
		while (s.i < s.samplepixels)
		{
			if (isARGB && thepicture.get(s.p) > 0 || !isARGB && thepicture.get(s.p + 3) > 0)
			{            
				if (isARGB)
				{
					s.al = thepicture.get(s.p);
					s.b = Std.int(biasvalue(thepicture.get(s.p + 3)));
					s.g = Std.int(biasvalue(thepicture.get(s.p + 2)));
					s.r = Std.int(biasvalue(thepicture.get(s.p + 1)));
				} 
				else
				{
					s.al = thepicture.get(s.p + 3);
					s.b = Std.int(biasvalue(thepicture.get(s.p + 2)));
					s.g = Std.int(biasvalue(thepicture.get(s.p + 1)));
					s.r = Std.int(biasvalue(thepicture.get(s.p)));
				}
			} 
			else
			{
				s.al=s.r=s.g=s.b=0;
			}
			
			var j = contest(s.al,s.b,s.g,s.r);
			altersingle(s.alpha, j, s.al, s.b, s.g, s.r);
			
			if (s.rad > 0)
			{
				alterneigh(s.rad,j,s.al,s.b,s.g,s.r);   /* alter neighbours */
			}
			
			s.p += s.step;
			
			while (s.p >= lengthcount)
			{
				s.p -= lengthcount;
			}
			
			s.i++;
			
			// FPE here if delta=0
			
			if (s.i % s.delta == 0)
			{
				s.alpha -= s.alpha / alphadec;
				s.radius -= s.radius / RADIUS_DEC;
				s.rad = Std.int(s.radius);
				if (s.rad <= 1)
				{
					s.rad = 0;
				}
				
				for (j in 0...(s.rad))
				{
					radpower[j] = Math.floor( s.alpha * (((s.rad * s.rad - j * j) * RAD_BIAS) / (s.rad * s.rad)) );
				}
			}
			
			if (Date.now().getTime() - quantizationStatus.chunkStartedTime > quantizationStatus.workChunkTimeSize)
			{
				return;
			}
		}
		
		if (verbose)
		{
			trace("finished 1D learning: final alpha=" + (1.0 * s.alpha / INIT_ALPHA) + " !\n");
		}
		
		s.finished = true;
	}
	
	private static function rgb2xyz(r:Int, g:Int, b: Int):Dynamic
	{
		var fR: Float = ( r / 255);        //R from 0 to 255
		var fG: Float = ( g / 255 );        //G from 0 to 255
		var fB: Float = ( b / 255 );        //B from 0 to 255
		
		fR = if ( fR > 0.04045 ) Math.pow(( fR + 0.055 ) / 1.055, 2.4)	else fR / 12.92;
		fG = if ( fG > 0.04045 ) Math.pow(( fG + 0.055 ) / 1.055, 2.4) else fG / 12.92;
		fB = if ( fB > 0.04045 ) Math.pow(( fB + 0.055 ) / 1.055, 2.4) else fB / 12.92;
		
		fR *= 100;
		fG *= 100;
		fB *= 100;
		
		//Observer. = 2°, Illuminant = D65
		
		return {
			x: fR * 0.4124 + fG * 0.3576 + fB * 0.1805,
			y: fR * 0.2126 + fG * 0.7152 + fB * 0.0722,
			z: fR * 0.0193 + fG * 0.1192 + fB * 0.9505	
		};
	}
	
	private static function xyz2cielab(x:Float, y: Float, z: Float):Dynamic
	{
		x /= 95.047;          //ref_X =  95.047   Observer= 2°, Illuminant= D65
		y /= 100.0;          //ref_Y = 100.000
		z /= 108.883;          //ref_Z = 108.883
		
		x = if ( x > 0.008856 ) Math.pow(x, 1/3) else ( 7.787 * x ) + ( 16 / 116 );
		y = if ( y > 0.008856 ) Math.pow(y, 1/3) else ( 7.787 * y ) + ( 16 / 116 );
		z = if ( z > 0.008856 ) Math.pow(z, 1/3) else ( 7.787 * z ) + ( 16 / 116 );
		
		return {
			l: ( 116 * y ) - 16,
			a: 500 * ( x - y ),
			b: 200 * ( y - z )
		};
	}
}