package org.freehep.graphicsio.gif;
/* 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.
* See also http://www.acm.org/~dekker/NEUQUANT.HTML
*
* 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.
*/
public class NeuQuant {
public static final int ncycles = 100; // no. of learning cycles
public static final int netsize = 255; // number of colours used
public static final int specials = 3; // number of reserved colours used
public static final int bgColour = specials-1; // reserved background colour
public static final int cutnetsize = netsize - specials;
public static final int maxnetpos = netsize-1;
public static final int initrad = netsize/8; // for 256 cols, radius starts at 32
public static final int radiusbiasshift = 6;
public static final int radiusbias = 1 << radiusbiasshift;
public static final int initBiasRadius = initrad*radiusbias;
public static final int radiusdec = 30; // factor of 1/30 each cycle
public static final int alphabiasshift = 10; // alpha starts at 1
public static final int initalpha = 1<<alphabiasshift; // biased by 10 bits
public static final double gamma = 1024.0;
public static final double beta = 1.0/1024.0;
public static final double betagamma = beta * gamma;
private double [] [] network = new double [netsize] [3]; // the network itself
protected int [] [] colormap = new int [netsize] [4]; // the network itself
private int [] netindex = new int [256]; // for network lookup - really 256
private double [] bias = new double [netsize]; // bias and freq arrays for learning
private double [] freq = new double [netsize];
// four primes near 500 - assume no image has a length so large
// that it is divisible by all four primes
public static final int prime1 = 499;
public static final int prime2 = 491;
public static final int prime3 = 487;
public static final int prime4 = 503;
public static final int maxprime= prime4;
protected int [][] pixels = null;
private int samplefac = 0;
public NeuQuant (int sample, int[][] pixels) {
if (sample < 1) throw new RuntimeException ("Sample must be 1..30");
if (sample > 30) throw new RuntimeException ("Sample must be 1..30");
samplefac = sample;
setPixels (pixels);
setUpArrays ();
}
public int getColorCount () {
return netsize;
}
public int[] getColorMap() {
// keep entry 0 free for transparent color
int[] c = new int[netsize+1];
c[0] = 0x00000000;
for (int i=0; i<netsize; i++) {
c[i+1] = (colormap[i][0] ) |
(colormap[i][1] << 8) |
(colormap[i][2] << 16) |
(0xFF << 24);
}
return c;
}
protected void setUpArrays () {
network [0] [0] = 0.0; // black
network [0] [1] = 0.0;
network [0] [2] = 0.0;
network [1] [0] = 1.0; // white
network [1] [1] = 1.0;
network [1] [2] = 1.0;
// RESERVED bgColour // background
for (int i=0; i<specials; i++) {
freq[i] = 1.0 / netsize;
bias[i] = 0.0;
}
for (int i=specials; i<netsize; i++) {
double [] p = network [i];
p[0] = (256.0 * (i-specials)) / cutnetsize;
p[1] = (256.0 * (i-specials)) / cutnetsize;
p[2] = (256.0 * (i-specials)) / cutnetsize;
freq[i] = 1.0 / netsize;
bias[i] = 0.0;
}
}
private void setPixels (int[][] pixels) {
if (pixels.length*pixels[0].length < maxprime) throw new RuntimeException ("Image is too small");
this.pixels = pixels;
}
public void init () {
learn ();
fix ();
inxbuild ();
}
private void altersingle(double alpha, int i, double b, double g, double r) {
// Move neuron i towards biased (b,g,r) by factor alpha
double [] n = network[i]; // alter hit neuron
n[0] -= (alpha*(n[0] - b));
n[1] -= (alpha*(n[1] - g));
n[2] -= (alpha*(n[2] - r));
}
private void alterneigh(double alpha, int rad, int i, double b, double g, double r) {
int lo = i-rad; if (lo<specials-1) lo=specials-1;
int hi = i+rad; if (hi>netsize) hi=netsize;
int j = i+1;
int k = i-1;
int q = 0;
while ((j<hi) || (k>lo)) {
double a = (alpha * (rad*rad - q*q)) / (rad*rad);
q ++;
if (j<hi) {
double [] p = network[j];
p[0] -= (a*(p[0] - b));
p[1] -= (a*(p[1] - g));
p[2] -= (a*(p[2] - r));
j++;
}
if (k>lo) {
double [] p = network[k];
p[0] -= (a*(p[0] - b));
p[1] -= (a*(p[1] - g));
p[2] -= (a*(p[2] - r));
k--;
}
}
}
private int contest (double b, double g, double r) { // Search for biased BGR values
// 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])
double bestd = Float.MAX_VALUE;
double bestbiasd = bestd;
int bestpos = -1;
int bestbiaspos = bestpos;
for (int i=specials; i<netsize; i++) {
double [] n = network[i];
double dist = n[0] - b; if (dist<0) dist = -dist;
double a = n[1] - g; if (a<0) a = -a;
dist += a;
a = n[2] - r; if (a<0) a = -a;
dist += a;
if (dist<bestd) {bestd=dist; bestpos=i;}
double biasdist = dist - bias [i];
if (biasdist<bestbiasd) {bestbiasd=biasdist; bestbiaspos=i;}
freq [i] -= beta * freq [i];
bias [i] += betagamma * freq [i];
}
freq[bestpos] += beta;
bias[bestpos] -= betagamma;
return bestbiaspos;
}
private int specialFind (double b, double g, double r) {
for (int i=0; i<specials; i++) {
double [] n = network[i];
if (n[0] == b && n[1] == g && n[2] == r) return i;
}
return -1;
}
private void learn() {
int biasRadius = initBiasRadius;
int alphadec = 30 + ((samplefac-1)/3);
int lengthcount = pixels.length*pixels[0].length;
int samplepixels = lengthcount / samplefac;
int delta = samplepixels / ncycles;
int alpha = initalpha;
int i = 0;
int rad = biasRadius >> radiusbiasshift;
if (rad <= 1) rad = 0;
// System.err.println("beginning 1D learning: samplepixels=" + samplepixels + " rad=" + rad);
int step = 0;
int pos = 0;
if ((lengthcount%prime1) != 0) step = prime1;
else {
if ((lengthcount%prime2) !=0) step = prime2;
else {
if ((lengthcount%prime3) !=0) step = prime3;
else step = prime4;
}
}
i = 0;
while (i < samplepixels) {
int p = pixels [pos / pixels[0].length][pos % pixels[0].length];
int red = (p >> 16) & 0xff;
int green = (p >> 8) & 0xff;
int blue = (p ) & 0xff;
double b = blue;
double g = green;
double r = red;
if (i == 0) { // remember background colour
network [bgColour] [0] = b;
network [bgColour] [1] = g;
network [bgColour] [2] = r;
}
int j = specialFind (b, g, r);
j = j < 0 ? contest (b, g, r) : j;
if (j >= specials) { // don't learn for specials
double a = (1.0 * alpha) / initalpha;
altersingle (a, j, b, g, r);
if (rad > 0) alterneigh (a, rad, j, b, g, r); // alter neighbours
}
pos += step;
while (pos >= lengthcount) pos -= lengthcount;
i++;
if (i%delta == 0) {
alpha -= alpha / alphadec;
biasRadius -= biasRadius / radiusdec;
rad = biasRadius >> radiusbiasshift;
if (rad <= 1) rad = 0;
}
}
// System.err.println("finished 1D learning: final alpha=" + (1.0 * alpha)/initalpha + "!");
}
private void fix() {
for (int i=0; i<netsize; i++) {
for (int j=0; j<3; j++) {
int x = (int) (0.5 + network[i][j]);
if (x < 0) x = 0;
if (x > 255) x = 255;
colormap[i][j] = x;
}
colormap[i][3] = i;
}
}
private void inxbuild() {
// Insertion sort of network and building of netindex[0..255]
int previouscol = 0;
int startpos = 0;
for (int i=0; i<netsize; i++) {
int[] p = colormap[i];
int[] q = null;
int smallpos = i;
int smallval = p[1]; // index on g
// find smallest in i..netsize-1
for (int j=i+1; j<netsize; j++) {
q = colormap[j];
if (q[1] < smallval) { // index on g
smallpos = j;
smallval = q[1]; // index on g
}
}
q = colormap[smallpos];
// swap p (i) and q (smallpos) entries
if (i != smallpos) {
int 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 != previouscol) {
netindex[previouscol] = (startpos+i)>>1;
for (int j=previouscol+1; j<smallval; j++) netindex[j] = i;
previouscol = smallval;
startpos = i;
}
}
netindex[previouscol] = (startpos+maxnetpos)>>1;
for (int j=previouscol+1; j<256; j++) netindex[j] = maxnetpos; // really 256
}
public int convert (int pixel) {
int alfa = (pixel >> 24) & 0xff;
int r = (pixel >> 16) & 0xff;
int g = (pixel >> 8) & 0xff;
int b = (pixel ) & 0xff;
int i = inxsearch(b, g, r);
int bb = colormap[i][0];
int gg = colormap[i][1];
int rr = colormap[i][2];
return (alfa << 24) | (rr << 16) | (gg << 8) | (bb);
}
public int lookup (int pixel) {
int r = (pixel >> 16) & 0xff;
int g = (pixel >> 8) & 0xff;
int b = (pixel ) & 0xff;
int i = inxsearch(b, g, r);
// compensate for transparent color
return i+1;
}
private int not_used_slow_inxsearch(int b, int g, int r) {
// Search for BGR values 0..255 and return colour index
int bestd = 1000; // biggest possible dist is 256*3
int best = -1;
for (int i = 0; i<netsize; i++) {
int [] p = colormap[i];
int dist = p[1] - g;
if (dist<0) dist = -dist;
int a = p[0] - b; if (a<0) a = -a;
dist += a;
a = p[2] - r; if (a<0) a = -a;
dist += a;
if (dist<bestd) {bestd=dist; best=i;}
}
return best;
}
protected int inxsearch(int b, int g, int r) {
// Search for BGR values 0..255 and return colour index
int bestd = 1000; // biggest possible dist is 256*3
int best = -1;
int i = netindex[g]; // index on g
int j = i-1; // start at netindex[g] and work outwards
while ((i<netsize) || (j>=0)) {
if (i<netsize) {
int [] p = colormap[i];
int dist = p[1] - g; // inx key
if (dist >= bestd) i = netsize; // stop iter
else {
if (dist<0) dist = -dist;
int a = p[0] - b; if (a<0) a = -a;
dist += a;
if (dist<bestd) {
a = p[2] - r; if (a<0) a = -a;
dist += a;
if (dist<bestd) {bestd=dist; best=i;}
}
i++;
}
}
if (j>=0) {
int [] p = colormap[j];
int dist = g - p[1]; // inx key - reverse dif
if (dist >= bestd) j = -1; // stop iter
else {
if (dist<0) dist = -dist;
int a = p[0] - b; if (a<0) a = -a;
dist += a;
if (dist<bestd) {
a = p[2] - r; if (a<0) a = -a;
dist += a;
if (dist<bestd) {bestd=dist; best=j;}
}
j--;
}
}
}
return best;
}
}