1 // NeuQuant Neural-Net Quantization Algorithm
2 // ------------------------------------------
3 //
4 // Copyright (c) 1994 Anthony Dekker
5 //
6 // NEUQUANT Neural-Net quantization algorithm by Anthony Dekker, 1994.
7 // See "Kohonen neural networks for optimal colour quantization"
8 // in "Network: Computation in Neural Systems" Vol. 5 (1994) pp 351-367.
9 // for a discussion of the algorithm.
10 //
11 // Any party obtaining a copy of these files from the author, directly or
12 // indirectly, is granted, free of charge, a full and unrestricted irrevocable,
13 // world-wide, paid up, royalty-free, nonexclusive right and license to deal
14 // in this software and documentation files (the "Software"), including without
15 // limitation the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 // and/or sell copies of the Software, and to permit persons who receive
17 // copies from any such party to do so, with the only requirement being
18 // that this copyright notice remain intact.
19
20 ///////////////////////////////////////////////////////////////////////
21 // History
22 // -------
23 // January 2001: Adaptation of the Neural-Net Quantization Algorithm
24 // for the FreeImage 2 library
25 // Author: Herv� Drolon (drolon@infonie.fr)
26 // March 2004: Adaptation for the FreeImage 3 library (port to big endian processors)
27 // Author: Herv� Drolon (drolon@infonie.fr)
28 // April 2004: Algorithm rewritten as a C++ class.
29 // Fixed a bug in the algorithm with handling of 4-byte boundary alignment.
30 // Author: Herv� Drolon (drolon@infonie.fr)
31 ///////////////////////////////////////////////////////////////////////
32
33 #include "Quantizers.h"
34 #include "FreeImage.h"
35 #include "Utilities.h"
36
37
38 // Four primes near 500 - assume no image has a length so large
39 // that it is divisible by all four primes
40 // ==========================================================
41
42 #define prime1 499
43 #define prime2 491
44 #define prime3 487
45 #define prime4 503
46
47 // ----------------------------------------------------------------
48
NNQuantizer(int PaletteSize)49 NNQuantizer::NNQuantizer(int PaletteSize)
50 {
51 netsize = PaletteSize;
52 maxnetpos = netsize - 1;
53 initrad = netsize < 8 ? 1 : (netsize >> 3);
54 initradius = (initrad * radiusbias);
55
56 network = NULL;
57
58 network = (pixel *)malloc(netsize * sizeof(pixel));
59 bias = (int *)malloc(netsize * sizeof(int));
60 freq = (int *)malloc(netsize * sizeof(int));
61 radpower = (int *)malloc(initrad * sizeof(int));
62
63 if( !network || !bias || !freq || !radpower ) {
64 if(network) free(network);
65 if(bias) free(bias);
66 if(freq) free(freq);
67 if(radpower) free(radpower);
68 throw FI_MSG_ERROR_MEMORY;
69 }
70 }
71
~NNQuantizer()72 NNQuantizer::~NNQuantizer()
73 {
74 if(network) free(network);
75 if(bias) free(bias);
76 if(freq) free(freq);
77 if(radpower) free(radpower);
78 }
79
80 ///////////////////////////////////////////////////////////////////////////
81 // Initialise network in range (0,0,0) to (255,255,255) and set parameters
82 // -----------------------------------------------------------------------
83
initnet()84 void NNQuantizer::initnet() {
85 int i, *p;
86
87 for (i = 0; i < netsize; i++) {
88 p = network[i];
89 p[FI_RGBA_BLUE] = p[FI_RGBA_GREEN] = p[FI_RGBA_RED] = (i << (netbiasshift+8))/netsize;
90 freq[i] = intbias/netsize; /* 1/netsize */
91 bias[i] = 0;
92 }
93 }
94
95 ///////////////////////////////////////////////////////////////////////////////////////
96 // Unbias network to give byte values 0..255 and record position i to prepare for sort
97 // ------------------------------------------------------------------------------------
98
unbiasnet()99 void NNQuantizer::unbiasnet() {
100 int i, j, temp;
101
102 for (i = 0; i < netsize; i++) {
103 for (j = 0; j < 3; j++) {
104 // OLD CODE: network[i][j] >>= netbiasshift;
105 // Fix based on bug report by Juergen Weigert jw@suse.de
106 temp = (network[i][j] + (1 << (netbiasshift - 1))) >> netbiasshift;
107 if (temp > 255) temp = 255;
108 network[i][j] = temp;
109 }
110 network[i][3] = i; // record colour no
111 }
112 }
113
114 //////////////////////////////////////////////////////////////////////////////////
115 // Insertion sort of network and building of netindex[0..255] (to do after unbias)
116 // -------------------------------------------------------------------------------
117
inxbuild()118 void NNQuantizer::inxbuild() {
119 int i,j,smallpos,smallval;
120 int *p,*q;
121 int previouscol,startpos;
122
123 previouscol = 0;
124 startpos = 0;
125 for (i = 0; i < netsize; i++) {
126 p = network[i];
127 smallpos = i;
128 smallval = p[FI_RGBA_GREEN]; // index on g
129 // find smallest in i..netsize-1
130 for (j = i+1; j < netsize; j++) {
131 q = network[j];
132 if (q[FI_RGBA_GREEN] < smallval) { // index on g
133 smallpos = j;
134 smallval = q[FI_RGBA_GREEN]; // index on g
135 }
136 }
137 q = network[smallpos];
138 // swap p (i) and q (smallpos) entries
139 if (i != smallpos) {
140 j = q[FI_RGBA_BLUE]; q[FI_RGBA_BLUE] = p[FI_RGBA_BLUE]; p[FI_RGBA_BLUE] = j;
141 j = q[FI_RGBA_GREEN]; q[FI_RGBA_GREEN] = p[FI_RGBA_GREEN]; p[FI_RGBA_GREEN] = j;
142 j = q[FI_RGBA_RED]; q[FI_RGBA_RED] = p[FI_RGBA_RED]; p[FI_RGBA_RED] = j;
143 j = q[3]; q[3] = p[3]; p[3] = j;
144 }
145 // smallval entry is now in position i
146 if (smallval != previouscol) {
147 netindex[previouscol] = (startpos+i)>>1;
148 for (j = previouscol+1; j < smallval; j++)
149 netindex[j] = i;
150 previouscol = smallval;
151 startpos = i;
152 }
153 }
154 netindex[previouscol] = (startpos+maxnetpos)>>1;
155 for (j = previouscol+1; j < 256; j++)
156 netindex[j] = maxnetpos; // really 256
157 }
158
159 ///////////////////////////////////////////////////////////////////////////////
160 // Search for BGR values 0..255 (after net is unbiased) and return colour index
161 // ----------------------------------------------------------------------------
162
inxsearch(int b,int g,int r)163 int NNQuantizer::inxsearch(int b, int g, int r) {
164 int i, j, dist, a, bestd;
165 int *p;
166 int best;
167
168 bestd = 1000; // biggest possible dist is 256*3
169 best = -1;
170 i = netindex[g]; // index on g
171 j = i-1; // start at netindex[g] and work outwards
172
173 while ((i < netsize) || (j >= 0)) {
174 if (i < netsize) {
175 p = network[i];
176 dist = p[FI_RGBA_GREEN] - g; // inx key
177 if (dist >= bestd)
178 i = netsize; // stop iter
179 else {
180 i++;
181 if (dist < 0)
182 dist = -dist;
183 a = p[FI_RGBA_BLUE] - b;
184 if (a < 0)
185 a = -a;
186 dist += a;
187 if (dist < bestd) {
188 a = p[FI_RGBA_RED] - r;
189 if (a<0)
190 a = -a;
191 dist += a;
192 if (dist < bestd) {
193 bestd = dist;
194 best = p[3];
195 }
196 }
197 }
198 }
199 if (j >= 0) {
200 p = network[j];
201 dist = g - p[FI_RGBA_GREEN]; // inx key - reverse dif
202 if (dist >= bestd)
203 j = -1; // stop iter
204 else {
205 j--;
206 if (dist < 0)
207 dist = -dist;
208 a = p[FI_RGBA_BLUE] - b;
209 if (a<0)
210 a = -a;
211 dist += a;
212 if (dist < bestd) {
213 a = p[FI_RGBA_RED] - r;
214 if (a<0)
215 a = -a;
216 dist += a;
217 if (dist < bestd) {
218 bestd = dist;
219 best = p[3];
220 }
221 }
222 }
223 }
224 }
225 return best;
226 }
227
228 ///////////////////////////////
229 // Search for biased BGR values
230 // ----------------------------
231
contest(int b,int g,int r)232 int NNQuantizer::contest(int b, int g, int r) {
233 // finds closest neuron (min dist) and updates freq
234 // finds best neuron (min dist-bias) and returns position
235 // for frequently chosen neurons, freq[i] is high and bias[i] is negative
236 // bias[i] = gamma*((1/netsize)-freq[i])
237
238 int i,dist,a,biasdist,betafreq;
239 int bestpos,bestbiaspos,bestd,bestbiasd;
240 int *p,*f, *n;
241
242 bestd = ~(((int) 1)<<31);
243 bestbiasd = bestd;
244 bestpos = -1;
245 bestbiaspos = bestpos;
246 p = bias;
247 f = freq;
248
249 for (i = 0; i < netsize; i++) {
250 n = network[i];
251 dist = n[FI_RGBA_BLUE] - b;
252 if (dist < 0)
253 dist = -dist;
254 a = n[FI_RGBA_GREEN] - g;
255 if (a < 0)
256 a = -a;
257 dist += a;
258 a = n[FI_RGBA_RED] - r;
259 if (a < 0)
260 a = -a;
261 dist += a;
262 if (dist < bestd) {
263 bestd = dist;
264 bestpos = i;
265 }
266 biasdist = dist - ((*p)>>(intbiasshift-netbiasshift));
267 if (biasdist < bestbiasd) {
268 bestbiasd = biasdist;
269 bestbiaspos = i;
270 }
271 betafreq = (*f >> betashift);
272 *f++ -= betafreq;
273 *p++ += (betafreq << gammashift);
274 }
275 freq[bestpos] += beta;
276 bias[bestpos] -= betagamma;
277 return bestbiaspos;
278 }
279
280 ///////////////////////////////////////////////////////
281 // Move neuron i towards biased (b,g,r) by factor alpha
282 // ----------------------------------------------------
283
altersingle(int alpha,int i,int b,int g,int r)284 void NNQuantizer::altersingle(int alpha, int i, int b, int g, int r) {
285 int *n;
286
287 n = network[i]; // alter hit neuron
288 n[FI_RGBA_BLUE] -= (alpha * (n[FI_RGBA_BLUE] - b)) / initalpha;
289 n[FI_RGBA_GREEN] -= (alpha * (n[FI_RGBA_GREEN] - g)) / initalpha;
290 n[FI_RGBA_RED] -= (alpha * (n[FI_RGBA_RED] - r)) / initalpha;
291 }
292
293 ////////////////////////////////////////////////////////////////////////////////////
294 // Move adjacent neurons by precomputed alpha*(1-((i-j)^2/[r]^2)) in radpower[|i-j|]
295 // ---------------------------------------------------------------------------------
296
alterneigh(int rad,int i,int b,int g,int r)297 void NNQuantizer::alterneigh(int rad, int i, int b, int g, int r) {
298 int j, k, lo, hi, a;
299 int *p, *q;
300
301 lo = i - rad; if (lo < -1) lo = -1;
302 hi = i + rad; if (hi > netsize) hi = netsize;
303
304 j = i+1;
305 k = i-1;
306 q = radpower;
307 while ((j < hi) || (k > lo)) {
308 a = (*(++q));
309 if (j < hi) {
310 p = network[j];
311 p[FI_RGBA_BLUE] -= (a * (p[FI_RGBA_BLUE] - b)) / alpharadbias;
312 p[FI_RGBA_GREEN] -= (a * (p[FI_RGBA_GREEN] - g)) / alpharadbias;
313 p[FI_RGBA_RED] -= (a * (p[FI_RGBA_RED] - r)) / alpharadbias;
314 j++;
315 }
316 if (k > lo) {
317 p = network[k];
318 p[FI_RGBA_BLUE] -= (a * (p[FI_RGBA_BLUE] - b)) / alpharadbias;
319 p[FI_RGBA_GREEN] -= (a * (p[FI_RGBA_GREEN] - g)) / alpharadbias;
320 p[FI_RGBA_RED] -= (a * (p[FI_RGBA_RED] - r)) / alpharadbias;
321 k--;
322 }
323 }
324 }
325
326 /////////////////////
327 // Main Learning Loop
328 // ------------------
329
330 /**
331 Get a pixel sample at position pos. Handle 4-byte boundary alignment.
332 @param pos pixel position in a WxHx3 pixel buffer
333 @param b blue pixel component
334 @param g green pixel component
335 @param r red pixel component
336 */
getSample(long pos,int * b,int * g,int * r)337 void NNQuantizer::getSample(long pos, int *b, int *g, int *r) {
338 // get equivalent pixel coordinates
339 // - assume it's a 24-bit image -
340 int x = pos % img_line;
341 int y = pos / img_line;
342
343 BYTE *bits = FreeImage_GetScanLine(dib_ptr, y) + x;
344
345 *b = bits[FI_RGBA_BLUE] << netbiasshift;
346 *g = bits[FI_RGBA_GREEN] << netbiasshift;
347 *r = bits[FI_RGBA_RED] << netbiasshift;
348 }
349
learn(int sampling_factor)350 void NNQuantizer::learn(int sampling_factor) {
351 int i, j, b, g, r;
352 int radius, rad, alpha, step, delta, samplepixels;
353 int alphadec; // biased by 10 bits
354 long pos, lengthcount;
355
356 // image size as viewed by the scan algorithm
357 lengthcount = img_width * img_height * 3;
358
359 // number of samples used for the learning phase
360 samplepixels = lengthcount / (3 * sampling_factor);
361
362 // decrease learning rate after delta pixel presentations
363 delta = samplepixels / ncycles;
364 if(delta == 0) {
365 // avoid a 'divide by zero' error with very small images
366 delta = 1;
367 }
368
369 // initialize learning parameters
370 alphadec = 30 + ((sampling_factor - 1) / 3);
371 alpha = initalpha;
372 radius = initradius;
373
374 rad = radius >> radiusbiasshift;
375 if (rad <= 1) rad = 0;
376 for (i = 0; i < rad; i++)
377 radpower[i] = alpha*(((rad*rad - i*i)*radbias)/(rad*rad));
378
379 // initialize pseudo-random scan
380 if ((lengthcount % prime1) != 0)
381 step = 3*prime1;
382 else {
383 if ((lengthcount % prime2) != 0)
384 step = 3*prime2;
385 else {
386 if ((lengthcount % prime3) != 0)
387 step = 3*prime3;
388 else
389 step = 3*prime4;
390 }
391 }
392
393 i = 0; // iteration
394 pos = 0; // pixel position
395
396 while (i < samplepixels) {
397 // get next learning sample
398 getSample(pos, &b, &g, &r);
399
400 // find winning neuron
401 j = contest(b, g, r);
402
403 // alter winner
404 altersingle(alpha, j, b, g, r);
405
406 // alter neighbours
407 if (rad) alterneigh(rad, j, b, g, r);
408
409 // next sample
410 pos += step;
411 while (pos >= lengthcount) pos -= lengthcount;
412
413 i++;
414 if (i % delta == 0) {
415 // decrease learning rate and also the neighborhood
416 alpha -= alpha / alphadec;
417 radius -= radius / radiusdec;
418 rad = radius >> radiusbiasshift;
419 if (rad <= 1) rad = 0;
420 for (j = 0; j < rad; j++)
421 radpower[j] = alpha * (((rad*rad - j*j) * radbias) / (rad*rad));
422 }
423 }
424
425 }
426
427 //////////////
428 // Quantizer
429 // -----------
430
Quantize(FIBITMAP * dib,int ReserveSize,RGBQUAD * ReservePalette,int sampling)431 FIBITMAP* NNQuantizer::Quantize(FIBITMAP *dib, int ReserveSize, RGBQUAD *ReservePalette, int sampling) {
432
433 if ((!dib) || (FreeImage_GetBPP(dib) != 24)) {
434 return NULL;
435 }
436
437 // 1) Select a sampling factor in range 1..30 (input parameter 'sampling')
438 // 1 => slower, 30 => faster. Default value is 1
439
440
441 // 2) Get DIB parameters
442
443 dib_ptr = dib;
444
445 img_width = FreeImage_GetWidth(dib); // DIB width
446 img_height = FreeImage_GetHeight(dib); // DIB height
447 img_line = FreeImage_GetLine(dib); // DIB line length in bytes (should be equal to 3 x W)
448
449 // For small images, adjust the sampling factor to avoid a 'divide by zero' error later
450 // (see delta in learn() routine)
451 int adjust = (img_width * img_height) / ncycles;
452 if(sampling >= adjust)
453 sampling = 1;
454
455
456 // 3) Initialize the network and apply the learning algorithm
457
458 if( netsize > ReserveSize ) {
459 netsize -= ReserveSize;
460 initnet();
461 learn(sampling);
462 unbiasnet();
463 netsize += ReserveSize;
464 }
465
466 // 3.5) Overwrite the last few palette entries with the reserved ones
467 for (int i = 0; i < ReserveSize; i++) {
468 network[netsize - ReserveSize + i][FI_RGBA_BLUE] = ReservePalette[i].rgbBlue;
469 network[netsize - ReserveSize + i][FI_RGBA_GREEN] = ReservePalette[i].rgbGreen;
470 network[netsize - ReserveSize + i][FI_RGBA_RED] = ReservePalette[i].rgbRed;
471 network[netsize - ReserveSize + i][3] = netsize - ReserveSize + i;
472 }
473
474 // 4) Allocate a new 8-bit DIB
475
476 FIBITMAP *new_dib = FreeImage_Allocate(img_width, img_height, 8);
477
478 if (new_dib == NULL)
479 return NULL;
480
481 // 5) Write the quantized palette
482
483 RGBQUAD *new_pal = FreeImage_GetPalette(new_dib);
484
485 for (int j = 0; j < netsize; j++) {
486 new_pal[j].rgbBlue = (BYTE)network[j][FI_RGBA_BLUE];
487 new_pal[j].rgbGreen = (BYTE)network[j][FI_RGBA_GREEN];
488 new_pal[j].rgbRed = (BYTE)network[j][FI_RGBA_RED];
489 }
490
491 inxbuild();
492
493 // 6) Write output image using inxsearch(b,g,r)
494
495 for (WORD rows = 0; rows < img_height; rows++) {
496 BYTE *new_bits = FreeImage_GetScanLine(new_dib, rows);
497 BYTE *bits = FreeImage_GetScanLine(dib_ptr, rows);
498
499 for (WORD cols = 0; cols < img_width; cols++) {
500 new_bits[cols] = (BYTE)inxsearch(bits[FI_RGBA_BLUE], bits[FI_RGBA_GREEN], bits[FI_RGBA_RED]);
501
502 bits += 3;
503 }
504 }
505
506 return (FIBITMAP*) new_dib;
507 }
508