469 lines
10 KiB
C++
469 lines
10 KiB
C++
/*
|
|
quantizer.cpp - image quantizer. based on Antony Dekker original code
|
|
Copyright (C) 2015 Uncle Mike
|
|
|
|
This program is free software: you can redistribute it and/or modify
|
|
it under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
GNU General Public License for more details.
|
|
*/
|
|
|
|
#include <windows.h>
|
|
#include "cmdlib.h"
|
|
#include "mathlib.h"
|
|
#include "stringlib.h"
|
|
#include "scriplib.h"
|
|
#include "filesystem.h"
|
|
#include "imagelib.h"
|
|
|
|
#define netsize 256 // number of colours used
|
|
#define prime1 499
|
|
#define prime2 491
|
|
#define prime3 487
|
|
#define prime4 503
|
|
|
|
#define minpicturebytes (3*prime4) // minimum size for input image
|
|
|
|
#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 (1<<intbiasshift)
|
|
#define gammashift 10 // gamma = 1024
|
|
#define gamma (1<<gammashift)
|
|
#define betashift 10
|
|
#define beta (intbias>>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 (1<<radiusbiasshift)
|
|
#define initradius (initrad * radiusbias) // and decreases by a
|
|
#define radiusdec 30 // factor of 1/30 each cycle
|
|
|
|
// defs for decreasing alpha factor
|
|
#define alphabiasshift 10 // alpha starts at 1.0
|
|
#define initalpha (1<<alphabiasshift)
|
|
int alphadec; // biased by 10 bits
|
|
|
|
// radbias and alpharadbias used for radpower calculation
|
|
#define radbiasshift 8
|
|
#define radbias (1<<radbiasshift)
|
|
#define alpharadbshift (alphabiasshift+radbiasshift)
|
|
#define alpharadbias (1<<alpharadbshift)
|
|
|
|
// types and global variables
|
|
static byte *thepicture; // the input image itself
|
|
static int lengthcount; // lengthcount = H*W*3
|
|
static int samplefac; // sampling factor 1..30
|
|
static int network[netsize][4]; // the network itself
|
|
static int netindex[256]; // for network lookup - really 256
|
|
static int bias[netsize]; // bias and freq arrays for learning
|
|
static int freq[netsize];
|
|
static int radpower[initrad]; // radpower for precomputation
|
|
|
|
void initnet( byte *thepic, int len, int sample )
|
|
{
|
|
register int i, *p;
|
|
|
|
thepicture = thepic;
|
|
lengthcount = len;
|
|
samplefac = sample;
|
|
|
|
for( i = 0; i < netsize; i++ )
|
|
{
|
|
p = network[i];
|
|
p[0] = p[1] = p[2] = (i << (netbiasshift + 8)) / netsize;
|
|
freq[i] = intbias / netsize; // 1 / netsize
|
|
bias[i] = 0;
|
|
}
|
|
}
|
|
|
|
// unbias network to give byte values 0..255 and record position i to prepare for sort
|
|
void unbiasnet( void )
|
|
{
|
|
for( int i = 0; i < netsize; i++ )
|
|
{
|
|
for( int j = 0; j < 3; j++ )
|
|
{
|
|
int temp = (network[i][j] + (1 << (netbiasshift - 1))) >> netbiasshift;
|
|
if( temp > 255 ) temp = 255;
|
|
network[i][j] = temp;
|
|
}
|
|
network[i][3] = i; // record colour num
|
|
}
|
|
}
|
|
|
|
// insertion sort of network and building of netindex[0..255] (to do after unbias)
|
|
void inxbuild( void )
|
|
{
|
|
register int *p, *q;
|
|
register int i, j, smallpos, smallval;
|
|
int previouscol, startpos;
|
|
|
|
previouscol = 0;
|
|
startpos = 0;
|
|
|
|
for( i = 0; i < netsize; i++ )
|
|
{
|
|
p = network[i];
|
|
smallpos = i;
|
|
smallval = p[1]; // index on g
|
|
|
|
// find smallest in i..netsize-1
|
|
for( j = i + 1; j < netsize; j++ )
|
|
{
|
|
q = network[j];
|
|
|
|
if( q[1] < smallval )
|
|
{
|
|
// index on g
|
|
smallpos = j;
|
|
smallval = q[1]; // index on g
|
|
}
|
|
}
|
|
|
|
q = network[smallpos];
|
|
|
|
// swap p (i) and q (smallpos) entries
|
|
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 != previouscol )
|
|
{
|
|
netindex[previouscol] = (startpos+i) >> 1;
|
|
|
|
for( j = previouscol + 1; j < smallval; j++ )
|
|
netindex[j] = i;
|
|
|
|
previouscol = smallval;
|
|
startpos = i;
|
|
}
|
|
}
|
|
|
|
netindex[previouscol] = (startpos + maxnetpos)>>1;
|
|
|
|
for( j = previouscol + 1; j < 256; j++ )
|
|
netindex[j] = maxnetpos; // really 256
|
|
}
|
|
|
|
// search for BGR values 0..255 (after net is unbiased) and return colour index
|
|
int inxsearch( int r, int g, int b )
|
|
{
|
|
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 < netsize ) || ( j >= 0 ))
|
|
{
|
|
if( i < netsize )
|
|
{
|
|
p = network[i];
|
|
dist = p[1] - g; // inx key
|
|
|
|
if( dist >= bestd )
|
|
{
|
|
i = netsize; // stop iter
|
|
}
|
|
else
|
|
{
|
|
i++;
|
|
if( dist < 0 ) dist = -dist;
|
|
a = p[2] - b;
|
|
if( a < 0 ) a = -a;
|
|
dist += a;
|
|
|
|
if( dist < bestd )
|
|
{
|
|
a = p[0] - r;
|
|
if( a < 0 ) a = -a;
|
|
dist += a;
|
|
|
|
if( dist < bestd )
|
|
{
|
|
bestd = dist;
|
|
best = p[3];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if( j >= 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[2] - b;
|
|
if( a < 0 ) a = -a;
|
|
dist += a;
|
|
|
|
if( dist < bestd )
|
|
{
|
|
a = p[0] - r;
|
|
if( a < 0 ) a = -a;
|
|
dist += a;
|
|
if( dist < bestd )
|
|
{
|
|
bestd = dist;
|
|
best = p[3];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return best;
|
|
}
|
|
|
|
// search for biased BGR values
|
|
int contest( int r, int g, int b )
|
|
{
|
|
// 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])
|
|
|
|
register int *p, *f, *n;
|
|
register int i, dist, a, biasdist, betafreq;
|
|
int bestpos, bestbiaspos, bestd, bestbiasd;
|
|
|
|
bestd = ~(1<<31);
|
|
bestbiasd = bestd;
|
|
bestpos = -1;
|
|
bestbiaspos = bestpos;
|
|
p = bias;
|
|
f = freq;
|
|
|
|
for( i = 0; i < netsize; i++ )
|
|
{
|
|
n = network[i];
|
|
dist = n[2] - b;
|
|
if( dist < 0 ) dist = -dist;
|
|
a = n[1] - g;
|
|
if( a < 0 ) a = -a;
|
|
dist += a;
|
|
a = n[0] - r;
|
|
if( a < 0 ) a = -a;
|
|
dist += a;
|
|
|
|
if( dist < bestd )
|
|
{
|
|
bestd = dist;
|
|
bestpos = i;
|
|
}
|
|
|
|
biasdist = dist - ((*p) >> (intbiasshift - netbiasshift));
|
|
|
|
if( biasdist < bestbiasd )
|
|
{
|
|
bestbiasd = biasdist;
|
|
bestbiaspos = i;
|
|
}
|
|
|
|
betafreq = (*f >> betashift);
|
|
*f++ -= betafreq;
|
|
*p++ += (betafreq << gammashift);
|
|
}
|
|
|
|
freq[bestpos] += beta;
|
|
bias[bestpos] -= betagamma;
|
|
|
|
return bestbiaspos;
|
|
}
|
|
|
|
// move neuron i towards biased (b,g,r) by factor alpha
|
|
void altersingle( int alpha, int i, int r, int g, int b )
|
|
{
|
|
register int *n;
|
|
|
|
n = network[i]; // alter hit neuron
|
|
*n -= (alpha * (*n - r)) / initalpha;
|
|
n++;
|
|
*n -= (alpha * (*n - g)) / initalpha;
|
|
n++;
|
|
*n -= (alpha * (*n - b)) / initalpha;
|
|
}
|
|
|
|
// move adjacent neurons by precomputed alpha*(1-((i-j)^2/[r]^2)) in radpower[|i-j|]
|
|
void alterneigh( int rad, int i, int r, int g, int b )
|
|
{
|
|
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(( j < hi ) || ( k > lo ))
|
|
{
|
|
a = (*(++q));
|
|
|
|
if( j < hi )
|
|
{
|
|
p = network[j];
|
|
*p -= (a * (*p - r)) / alpharadbias;
|
|
p++;
|
|
*p -= (a * (*p - g)) / alpharadbias;
|
|
p++;
|
|
*p -= (a * (*p - b)) / alpharadbias;
|
|
j++;
|
|
}
|
|
|
|
if( k > lo )
|
|
{
|
|
p = network[k];
|
|
*p -= (a * (*p - r)) / alpharadbias;
|
|
p++;
|
|
*p -= (a * (*p - g)) / alpharadbias;
|
|
p++;
|
|
*p -= (a * (*p - b)) / alpharadbias;
|
|
k--;
|
|
}
|
|
}
|
|
}
|
|
|
|
// main Learning Loop
|
|
void learn( void )
|
|
{
|
|
register byte *p;
|
|
register int i, j, r, g, b;
|
|
int radius, rad, alpha, step;
|
|
int delta, samplepixels;
|
|
byte *lim;
|
|
|
|
alphadec = 30 + ((samplefac - 1) / 3);
|
|
p = thepicture;
|
|
lim = thepicture + lengthcount;
|
|
samplepixels = lengthcount / (samplefac * 4); // RGBA
|
|
delta = samplepixels / ncycles;
|
|
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( delta <= 0 ) return;
|
|
|
|
if(( lengthcount % prime1 ) != 0 )
|
|
{
|
|
step = prime1 * 4; // RGBA
|
|
}
|
|
else if(( lengthcount % prime2 ) != 0 )
|
|
{
|
|
step = prime2 * 4; // RGBA
|
|
}
|
|
else if(( lengthcount % prime3 ) != 0 )
|
|
{
|
|
step = prime3 * 4; // RGBA
|
|
}
|
|
else
|
|
{
|
|
step = prime4 * 4; // RGBA;
|
|
}
|
|
|
|
i = 0;
|
|
|
|
while( i < samplepixels )
|
|
{
|
|
r = p[0] << netbiasshift;
|
|
g = p[1] << netbiasshift;
|
|
b = p[2] << netbiasshift;
|
|
j = contest( r, g, b );
|
|
|
|
altersingle( alpha, j, r, g, b );
|
|
if( rad ) alterneigh( rad, j, r, g, b ); // 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 ));
|
|
}
|
|
}
|
|
}
|
|
|
|
// returns the actual number of palette entries.
|
|
rgbdata_t *Image_Quantize( rgbdata_t *pic )
|
|
{
|
|
rgbdata_t *out;
|
|
int i, samples;
|
|
|
|
if( !pic ) return NULL;
|
|
|
|
// already quantized?
|
|
if( FBitSet( pic->flags, IMAGE_QUANTIZED ))
|
|
return pic;
|
|
|
|
out = Image_Alloc( pic->width, pic->height, true );
|
|
out->flags = pic->flags;
|
|
|
|
if( FBitSet( pic->flags, IMAGE_HAS_COLOR ))
|
|
samples = 1; // maximum quality
|
|
else samples = 10; // fast mode
|
|
|
|
initnet( pic->buffer, pic->size, samples );
|
|
learn();
|
|
unbiasnet();
|
|
|
|
for( i = 0; i < netsize; i++ )
|
|
{
|
|
out->palette[i*3+0] = network[i][0]; // red
|
|
out->palette[i*3+1] = network[i][1]; // green
|
|
out->palette[i*3+2] = network[i][2]; // blue
|
|
}
|
|
|
|
inxbuild();
|
|
|
|
for( i = 0; i < pic->width * pic->height; i++ )
|
|
out->buffer[i] = inxsearch( pic->buffer[i*4+0], pic->buffer[i*4+1], pic->buffer[i*4+2] );
|
|
Mem_Free( pic ); // release RGBA image
|
|
|
|
return out;
|
|
} |