1// tQuantizeNeu.cpp 
2// 
3// This module implements the NeuQuant (Neural-Net) Quantization Algorithm written by Anthony Dekker. Only the parts 
4// modified by me are under the ISC licence. The majority is under what looks like the MIT licence, The original 
5// author's licence can be found below. Modifications include: 
6// * Placing it in a namespace. 
7// * Consolidating the state parameters so that it is threadsafe (no global state). 
8// * Bridging to a standardized Tacent interface. 
9// * Replacing the inxsearch and inxbuild with red-mean perceptual colour distance metric to choose best colours. 
10// * Support for an arbitrary number of colours between 2 and 256. 
11// 
12// The algrithm works well for larger numbers of colours (generally 128 to 256 or 255) but it can handle values as 
13// low as 2. 
14// 
15// Modifications Copyright (c) 2022-2024 Tristan Grimmer. 
16// Permission to use, copy, modify, and/or distribute this software for any purpose with or without fee is hereby 
17// granted, provided that the above copyright notice and this permission notice appear in all copies. 
18// 
19// THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL 
20// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, 
21// INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
22// AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR 
23// PERFORMANCE OF THIS SOFTWARE. 
24// 
25// Here is the original copyright: 
26// 
27// Copyright (c) 1994 Anthony Dekker 
28// 
29// NEUQUANT Neural-Net quantization algorithm by Anthony Dekker, 1994. See "Kohonen neural networks for optimal colour 
30// quantization" in "Network: Computation in Neural Systems" Vol. 5 (1994) pp351-367. for a discussion of the algorithm. 
31// See also http://members.ozemail.com.au/~dekker/NEUQUANT.HTML. 
32// 
33// Any party obtaining a copy of these files from the author, directly or indirectly, is granted, free of charge, a full 
34// and unrestricted irrevocable, world-wide, paid up, royalty-free, nonexclusive right and license to deal in this 
35// software and documentation files (the "Software"), including without limitation the rights to use, copy, modify, 
36// merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons who receive copies 
37// from any such party to do so, with the only requirement being that this copyright notice remain intact. 
38 
39#include <cstdint> 
40#include <Math/tColour.h> 
41#include "Image/tQuantize.h" 
42namespace tImage { 
43 
44 
45namespace tQuantizeNeu 
46{ 
47 // For 256 colours, fixed arrays need 8kb, plus space for the image. 
48 const int maxnetsize = 256; // Maximum network size (number of colours). 
49 
50 // Four primes near 500 - assume no image has a length so large that it is divisible by all four primes. 
51 const int prime1 = 499; 
52 const int prime2 = 491; 
53 const int prime3 = 487; 
54 const int prime4 = 503; 
55 
56 const int minpicturebytes = (3 * prime4); // Minimum size for input image. 
57 
58 // Network Definitions. 
59 const int maxnetpos = maxnetsize - 1; 
60 const int netbiasshift = 4; // Bias for colour values. 
61 const int ncycles = 100; // Number of learning cycles. 
62 
63 // Defs for freq and bias. 
64 const int intbiasshift = 16; // Bias for fractions. 
65 const int intbias = 65536; // 1 << intbiasshift. 
66 const int gammashift = 10; // Gamma = 1024. 
67 const int gamma = 1024; // 1 << gammashift. 
68 const int betashift = 10; 
69 const int beta = 64; // intbias >> betashift beta = 1/1024. 
70 const int betagamma = 65536; // intbias << (gammashift - betashift). 
71 
72 // Defs for decreasing radius factor. 
73 const int initrad = 32; // netsize >> 3 for 256 cols, radius starts 
74 const int radiusbiasshift = 6; // at 32.0 biased by 6 bits. 
75 const int radiusbias = 64; // 1 << radiusbiasshift 
76 const int initradius = 2048; // initrad * radiusbias and decreases by a 
77 const int radiusdec = 30; // factor of 1/30 each cycle. 
78 
79 // Defs for decreasing alpha factor. 
80 const int alphabiasshift = 10; // Alpha starts at 1.0 
81 const int initalpha = 1024; // 1 << alphabiasshift 
82 
83 // Radbias and alpharadbias used for radpower calculation. 
84 const int radbiasshift = 8; 
85 const int radbias = 256; // 1 << radbiasshift 
86 const int alpharadbshift = 18; // alphabiasshift + radbiasshift 
87 const int alpharadbias = 262144; 
88 
89 struct State 
90 { 
91 int netsize = maxnetsize; // This defaults to maxnetsize but can be reduced to as low as 2 colours. 
92 int alphadec; // Biased by 10 bits. 
93 unsigned char *thepicture; // The input image itself. 
94 int lengthcount; // lengthcount = H*W*3. 
95 int samplefac; // Sampling factor 1..30. 
96 typedef int pixel_bgr[4]; // BGRc 
97 pixel_bgr network[maxnetsize]; // The network itself. 
98 int netindex[maxnetsize]; // For network lookup - really 256. 
99 int bias[maxnetsize]; // Bias and freq arrays for learning. 
100 int freq[maxnetsize]; // Frequency array for learning. 
101 int radpower[initrad]; // radpower for precomputation. 
102 }; 
103 
104 // Initialise network in range (0,0,0) to (255,255,255) and set parameters. 
105 int getNetwork(State&, int i, int j); 
106 void initnet(State&, unsigned char *thepic, int len, int sample); 
107 
108 // Unbias network to give byte values 0..255 and record position i to prepare for sort. 
109 void unbiasnet(State&); 
110 
111 // This gets the palette in the out variable. 
112 int getColourMap(State&, tColour3b* out); 
113 
114 // Insertion sort of network and building of netindex[0..255] (to do after unbias). 
115 // We don't call this function since we do an exhaustive red-mean distance check. 
116 void inxbuild(State&); 
117 
118 // Search for BGR values 0..255 (after net is unbiased) and return colour index. 
119 // We don't call this function since we do an exhaustive red-mean distance check. 
120 int inxsearch(State&, int b, int g, int r); 
121 
122 // Main Learning Loop. 
123 void learn(State&); 
124 
125 int contest(State&, int b, int g, int r); 
126 void altersingle(State&, int alpha, int i, int b, int g, int r); 
127 void alterneigh(State&, int rad, int i, int b, int g, int r); 
128 
129 int FindIndexOfClosestColour_Redmean(const tColour3b* searchSpace, int searchSize, const tColour3b& searchColour); 
130} 
131 
132 
133int tQuantizeNeu::getNetwork(State& state, int i, int j) 
134{ 
135 return state.network[i][j]; 
136} 
137 
138 
139void tQuantizeNeu::initnet(State& state, unsigned char *thepic, int len, int sample) 
140{ 
141 // Initialise network in range (0,0,0) to (255,255,255) and set parameters. 
142 int i; 
143 int *p; 
144 
145 state.thepicture = thepic; 
146 state.lengthcount = len; 
147 state.samplefac = sample; 
148 
149 for (i = 0; i < state.netsize; i++) 
150 { 
151 p = state.network[i]; 
152 p[0] = p[1] = p[2] = (i << (netbiasshift + 8)) / state.netsize; 
153 state.freq[i] = intbias / state.netsize; /* 1/netsize */ 
154 state.bias[i] = 0; 
155 } 
156} 
157 
158 
159void tQuantizeNeu::unbiasnet(State& state) 
160{ 
161 // Unbias network to give byte values 0..255 and record position i to prepare for sort. 
162 int i, j, temp; 
163 
164 for (i = 0; i < state.netsize; i++) 
165 { 
166 for (j = 0; j < 3; j++) 
167 { 
168 // OLD CODE: network[i][j] >>= netbiasshift; 
169 // Fix based on bug report by Juergen Weigert jw@suse.de. 
170 temp = (state.network[i][j] + (1 << (netbiasshift - 1))) >> netbiasshift; 
171 if (temp > 255) temp = 255; 
172 state.network[i][j] = temp; 
173 } 
174 
175 // Record colour no. 
176 state.network[i][3] = i; 
177 } 
178} 
179 
180 
181int tQuantizeNeu::getColourMap(State& state, tColour3b* out) 
182{ 
183 // Output colour map. The palette. 
184 int index[maxnetsize]; 
185 for (int i = 0; i < state.netsize; i++) 
186 index[state.network[i][3]] = i; 
187 
188 for (int j = 0; j < state.netsize; j++) 
189 { 
190 out[j].R = uint8(state.network[j][0]); 
191 out[j].G = uint8(state.network[j][1]); 
192 out[j].B = uint8(state.network[j][2]); 
193 } 
194 return state.netsize; 
195} 
196 
197 
198void tQuantizeNeu::inxbuild(State& state) 
199{ 
200 // Insertion sort of network and building of netindex[0..255] (to do after unbias). 
201 int i, j, smallpos, smallval; 
202 int *p, *q; 
203 int previouscol, startpos; 
204 
205 previouscol = 0; 
206 startpos = 0; 
207 for (i = 0; i < state.netsize; i++) 
208 { 
209 p = state.network[i]; 
210 smallpos = i; 
211 smallval = p[1]; // Index on g. 
212 
213 // Find smallest in i..netsize-1. 
214 for (j = i + 1; j < state.netsize; j++) 
215 { 
216 q = state.network[j]; 
217 if (q[1] < smallval) // Index on g. 
218 { 
219 smallpos = j; 
220 smallval = q[1]; // Index on g. 
221 } 
222 } 
223 q = state.network[smallpos]; 
224 
225 // Swap p (i) and q (smallpos) entries. 
226 if (i != smallpos) 
227 { 
228 j = q[0]; 
229 q[0] = p[0]; 
230 p[0] = j; 
231 j = q[1]; 
232 q[1] = p[1]; 
233 p[1] = j; 
234 j = q[2]; 
235 q[2] = p[2]; 
236 p[2] = j; 
237 j = q[3]; 
238 q[3] = p[3]; 
239 p[3] = j; 
240 } 
241 
242 // Smallval entry is now in position i. 
243 if (smallval != previouscol) 
244 { 
245 state.netindex[previouscol] = (startpos + i) >> 1; 
246 for (j = previouscol + 1; j < smallval; j++) state.netindex[j] = i; 
247 previouscol = smallval; 
248 startpos = i; 
249 } 
250 } 
251 state.netindex[previouscol] = (startpos + maxnetpos) >> 1; 
252 for (j = previouscol + 1; j < 256; j++) 
253 state.netindex[j] = maxnetpos; 
254} 
255 
256 
257int tQuantizeNeu::inxsearch(State& state, int b, int g, int r) 
258{ 
259 // Search for BGR values 0..255 (after net is unbiased) and return colour index. 
260 int i, j, dist, a, bestd; 
261 int *p; 
262 int best; 
263 
264 bestd = 1000; // Biggest possible dist is 256*3. 
265 best = -1; 
266 i = state.netindex[g]; // Index on g. 
267 j = i - 1; // Start at netindex[g] and work outwards. 
268 
269 while ((i < state.netsize) || (j >= 0)) 
270 { 
271 if (i < state.netsize) 
272 { 
273 p = state.network[i]; 
274 dist = p[1] - g; // Inx key. 
275 if (dist >= bestd) 
276 i = state.netsize; // Stop iter. 
277 else 
278 { 
279 i++; 
280 if (dist < 0) dist = -dist; 
281 a = p[0] - b; 
282 if (a < 0) a = -a; 
283 dist += a; 
284 if (dist < bestd) 
285 { 
286 a = p[2] - r; 
287 if (a < 0) a = -a; 
288 dist += a; 
289 if (dist < bestd) 
290 { 
291 bestd = dist; 
292 best = p[3]; 
293 } 
294 } 
295 } 
296 } 
297 if (j >= 0) 
298 { 
299 p = state.network[j]; 
300 dist = g - p[1]; // Inx key - reverse dif. 
301 if (dist >= bestd) 
302 j = -1; // Stop iter. 
303 else 
304 { 
305 j--; 
306 if (dist < 0) dist = -dist; 
307 a = p[0] - b; 
308 if (a < 0) a = -a; 
309 dist += a; 
310 if (dist < bestd) 
311 { 
312 a = p[2] - r; 
313 if (a < 0) a = -a; 
314 dist += a; 
315 if (dist < bestd) 
316 { 
317 bestd = dist; 
318 best = p[3]; 
319 } 
320 } 
321 } 
322 } 
323 } 
324 return best; 
325} 
326 
327 
328int tQuantizeNeu::contest(State& state, int b, int g, int r) 
329{ 
330 // Search for biased BGR values. Finds closest neuron (min dist) and updates freq finds best neuron (min dist-bias) 
331 // and returns position for frequently chosen neurons, freq[i] is high and bias[i] is negative. 
332 // bias[i] = gamma*((1/netsize)-freq[i]) 
333 int i, dist, a, biasdist, betafreq; 
334 int bestpos, bestbiaspos, bestd, bestbiasd; 
335 int *p, *f, *n; 
336 
337 bestd = ~(((int) 1) << 31); 
338 bestbiasd = bestd; 
339 bestpos = -1; 
340 bestbiaspos = bestpos; 
341 p = state.bias; 
342 f = state.freq; 
343 
344 for (i = 0; i < state.netsize; i++) 
345 { 
346 n = state.network[i]; 
347 dist = n[0] - b; 
348 if (dist < 0) 
349 dist = -dist; 
350 
351 a = n[1] - g; 
352 if (a < 0) 
353 a = -a; 
354 
355 dist += a; 
356 a = n[2] - r; 
357 if (a < 0) 
358 a = -a; 
359 
360 dist += a; 
361 if (dist < bestd) 
362 { 
363 bestd = dist; 
364 bestpos = i; 
365 } 
366 biasdist = dist - ((*p) >> (intbiasshift - netbiasshift)); 
367 if (biasdist < bestbiasd) 
368 { 
369 bestbiasd = biasdist; 
370 bestbiaspos = i; 
371 } 
372 betafreq = (*f >> betashift); 
373 *f++ -= betafreq; 
374 *p++ += (betafreq << gammashift); 
375 } 
376 state.freq[bestpos] += beta; 
377 state.bias[bestpos] -= betagamma; 
378 return (bestbiaspos); 
379} 
380 
381 
382void tQuantizeNeu::altersingle(State& state, int alpha, int i, int b, int g, int r) 
383{ 
384 // Move neuron i towards biased (b,g,r) by factor alpha. 
385 int *n; 
386 
387 n = state.network[i]; // Alter hit neuron. 
388 *n -= (alpha * (*n - b)) / initalpha; 
389 n++; 
390 *n -= (alpha * (*n - g)) / initalpha; 
391 n++; 
392 *n -= (alpha * (*n - r)) / initalpha; 
393} 
394 
395 
396void tQuantizeNeu::alterneigh(State& state, int rad, int i, int b, int g, int r) 
397{ 
398 // Move adjacent neurons by precomputed alpha*(1-((i-j)^2/[r]^2)) in radpower[|i-j|]. 
399 int j, k, lo, hi, a; 
400 int *p, *q; 
401 
402 lo = i - rad; 
403 if (lo < -1) 
404 lo = -1; 
405 
406 hi = i + rad; 
407 if (hi > state.netsize) 
408 hi = state.netsize; 
409 
410 j = i + 1; 
411 k = i - 1; 
412 q = state.radpower; 
413 while ((j < hi) || (k > lo)) 
414 { 
415 a = (*(++q)); 
416 if (j < hi) 
417 { 
418 p = state.network[j]; 
419 *p -= (a * (*p - b)) / alpharadbias; 
420 p++; 
421 *p -= (a * (*p - g)) / alpharadbias; 
422 p++; 
423 *p -= (a * (*p - r)) / alpharadbias; 
424 j++; 
425 } 
426 if (k > lo) 
427 { 
428 p = state.network[k]; 
429 *p -= (a * (*p - b)) / alpharadbias; 
430 p++; 
431 *p -= (a * (*p - g)) / alpharadbias; 
432 p++; 
433 *p -= (a * (*p - r)) / alpharadbias; 
434 k--; 
435 } 
436 } 
437} 
438 
439 
440void tQuantizeNeu::learn(State& state) 
441{ 
442 // Main Learning Loop. 
443 int i, j, b, g, r; 
444 int radius, rad, alpha, step, delta, samplepixels; 
445 unsigned char *p; 
446 unsigned char *lim; 
447 
448 state.alphadec = 30 + ((state.samplefac - 1) / 3); 
449 p = state.thepicture; 
450 lim = state.thepicture + state.lengthcount; 
451 samplepixels = state.lengthcount / (3 * state.samplefac); 
452 delta = samplepixels / ncycles; 
453 alpha = initalpha; 
454 radius = initradius; 
455 
456 rad = radius >> radiusbiasshift; 
457 if (rad <= 1) 
458 rad = 0; 
459 
460 for (i = 0; i < rad; i++) 
461 state.radpower[i] = alpha * (((rad * rad - i * i) * radbias) / (rad * rad)); 
462 
463 if ((state.lengthcount % prime1) != 0) 
464 { 
465 step = 3 * prime1; 
466 } 
467 else 
468 { 
469 if ((state.lengthcount % prime2) != 0) 
470 { 
471 step = 3 * prime2; 
472 } 
473 else 
474 { 
475 if ((state.lengthcount % prime3) != 0) 
476 step = 3 * prime3; 
477 else 
478 step = 3 * prime4; 
479 } 
480 } 
481 
482 i = 0; 
483 while (i < samplepixels) 
484 { 
485 b = p[0] << netbiasshift; 
486 g = p[1] << netbiasshift; 
487 r = p[2] << netbiasshift; 
488 j = contest(state, b, g, r); 
489 
490 altersingle(state, alpha, i: j, b, g, r); 
491 if (rad) 
492 alterneigh(state, rad, i: j, b, g, r); // Alter neighbours. 
493 
494 p += step; 
495 if (p >= lim) 
496 p -= state.lengthcount; 
497 
498 i++; 
499 if (i % delta == 0) 
500 { 
501 alpha -= alpha / state.alphadec; 
502 radius -= radius / radiusdec; 
503 rad = radius >> radiusbiasshift; 
504 if (rad <= 1) 
505 rad = 0; 
506 
507 for (j = 0; j < rad; j++) 
508 state.radpower[j] = alpha * (((rad * rad - j * j) * radbias) / (rad * rad)); 
509 } 
510 } 
511} 
512 
513 
514int tQuantizeNeu::FindIndexOfClosestColour_Redmean(const tColour3b* searchSpace, int searchSize, const tColour3b& colour) 
515{ 
516 float closest = 1000.0f; 
517 int closestIndex = -1; 
518 
519 for (int i = 0; i < searchSize; i++) 
520 { 
521 float diff = tMath::tColourDiffRedmean(a: colour, b: searchSpace[i]); 
522 if (diff < closest) 
523 { 
524 closest = diff; 
525 closestIndex = i; 
526 } 
527 } 
528 return closestIndex; 
529} 
530 
531 
532// 
533// The functions below make up the external interface. 
534// 
535 
536 
537bool tQuantizeNeu::QuantizeImage 
538( 
539 int numColours, int width, int height, const tPixel3b* pixels, tColour3b* destPalette, uint8* destIndices, 
540 bool checkExact, int sampleFactor 
541) 
542{ 
543 if ((numColours < 2) || (numColours > 256) || (width <= 0) || (height <= 0) || !pixels || !destPalette || !destIndices) 
544 return false; 
545 
546 if ((sampleFactor < 1) || (sampleFactor > 30)) 
547 return false; 
548 
549 if (checkExact) 
550 { 
551 bool success = tQuantize::QuantizeImageExact(numColours, width, height, pixels, destPalette, destIndices); 
552 if (success) 
553 return true; 
554 } 
555 
556 State state; 
557 state.netsize = numColours; 
558 
559 initnet(state, thepic: (uint8*)pixels, len: width*height*3, sample: sampleFactor); 
560 learn(state); 
561 unbiasnet(state); 
562 int resultNumColours = getColourMap(state, out: destPalette); 
563 
564 // Exhaustive redmean is better. 
565 // inxbuild(state); 
566 for (int y = 0; y < height; y++) 
567 { 
568 for (int x = 0; x < width; x++) 
569 { 
570 const tPixel3b& pixel = pixels[x + y*width]; 
571 destIndices[x + y*width] = FindIndexOfClosestColour_Redmean(searchSpace: destPalette, searchSize: numColours, colour: pixel); 
572 
573 // Exhaustive redmean is better. 
574 // destIndices[x + y*width] = inxsearch(state, pixel.R, pixel.G, pixel.B); 
575 } 
576 } 
577 
578 return (resultNumColours > 0); 
579} 
580 
581 
582bool tQuantizeNeu::QuantizeImage 
583( 
584 int numColours, int width, int height, const tPixel4b* pixels, tColour3b* destPalette, uint8* destIndices, 
585 bool checkExact, int sampleFactor 
586) 
587{ 
588 if ((numColours < 2) || (numColours > 256) || (width <= 0) || (height <= 0) || !pixels || !destPalette || !destIndices) 
589 return false; 
590 
591 tPixel3b* pixels3 = new tPixel3b[width*height]; 
592 for (int p = 0; p < width*height; p++) 
593 pixels3[p].Set( r: pixels[p].R, g: pixels[p].G, b: pixels[p].B ); 
594 
595 bool success = QuantizeImage(numColours, width, height, pixels: pixels3, destPalette, destIndices, checkExact, sampleFactor); 
596 delete[] pixels3; 
597 return success; 
598} 
599 
600 
601} 
602