1// tQuantizeSpatial.cpp 
2// 
3// This module implements a spatial quantize algorithm written by Derrick Coetzee. Only the parts modified by me are 
4// under the ISC licence. The majority is under what looks like the MIT licence, The original author's licence can be 
5// found below. Modifications include placing it in a namespace, putting what used to be in 'main' into a function, 
6// modifying it so that all random numbers and operators are deterministic, indenting it for readability, and making 
7// sure there is no global state so calls are threadsafe. 
8// 
9// The algrithm works well for smaller numbers of colours (generally 32 or fewer) but it can handle from 2 to 256 
10// colours. The official name of this algorith is 'scolorq'. Running on colours more than 32 takes a LONG time. 
11// See https://github.com/samhocevar/scolorq/blob/master/spatial_color_quant.cpp 
12// 
13// Modifications Copyright (c) 2022-2024 Tristan Grimmer. 
14// Permission to use, copy, modify, and/or distribute this software for any purpose with or without fee is hereby 
15// granted, provided that the above copyright notice and this permission notice appear in all copies. 
16// 
17// THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL 
18// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, 
19// INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
20// AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR 
21// PERFORMANCE OF THIS SOFTWARE. 
22// 
23// Here is the original copyright: 
24// 
25// Copyright (c) 2006 Derrick Coetzee 
26// 
27// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated 
28// documentation files (the "Software"), to deal in the Software without restriction, including without limitation the 
29// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to 
30// permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright 
31// notice and this permission notice shall be included in all copies or substantial portions of the Software. 
32// 
33// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE 
34// WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR 
35// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR 
36// OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 
37 
38#include <deque> 
39#include <random> 
40#include <cmath> 
41#include <algorithm> 
42#include <Math/tColour.h> 
43#include <Math/tRandom.h> 
44#include "Image/tQuantize.h" 
45namespace tImage { 
46  
47  
48namespace tQuantizeSpatial 
49{ 
50 template <typename T, int length> class vector_fixed 
51 { 
52 public: 
53 vector_fixed() { for (int i=0; i<length; i++) data[i] = 0; } 
54 vector_fixed(const vector_fixed<T, length>& rhs) { for (int i=0; i<length; i++) data[i] = rhs.data[i]; } 
55 vector_fixed(const std::vector<T>& rhs) { for (int i=0; i<length; i++) data[i] = rhs[i]; } 
56 
57 T& operator()(int i) { return data[i]; } 
58 int get_length() { return length; } 
59 T norm_squared() { T result = 0; for (int i=0; i<length; i++) result += (*this)(i) * (*this)(i); return result; } 
60 
61 vector_fixed<T, length>& operator=(const vector_fixed<T, length> rhs) { for (int i=0; i<length; i++) data[i] = rhs.data[i]; return *this; } 
62 vector_fixed<T, length> direct_product(vector_fixed<T, length>& rhs) { vector_fixed<T, length> result; for (int i=0; i<length; i++) result(i) = (*this)(i) * rhs(i); return result; } 
63 
64 double dot_product(vector_fixed<T, length> rhs) { T result = 0; for (int i=0; i<length; i++) result += (*this)(i) * rhs(i); return result; } 
65 vector_fixed<T, length>& operator+=(vector_fixed<T, length> rhs) { for (int i=0; i<length; i++) data[i] += rhs(i); return *this; } 
66 vector_fixed<T, length> operator+(vector_fixed<T, length> rhs) { vector_fixed<T, length> result(*this); result += rhs; return result; } 
67 vector_fixed<T, length>& operator-=(vector_fixed<T, length> rhs) { for (int i=0; i<length; i++) data[i] -= rhs(i); return *this; } 
68 vector_fixed<T, length> operator-(vector_fixed<T, length> rhs) { vector_fixed<T, length> result(*this); result -= rhs; return result; } 
69 vector_fixed<T, length>& operator*=(T scalar) { for (int i=0; i<length; i++) data[i] *= scalar; return *this; } 
70 vector_fixed<T, length> operator*(T scalar) { vector_fixed<T, length> result(*this); result *= scalar; return result; } 
71 
72 private: 
73 T data[length]; 
74 }; 
75 
76 template <typename T, int length> vector_fixed<T, length> operator*(T scalar, vector_fixed<T, length> vec) { return vec*scalar; } 
77 
78 template <typename T> class array2d 
79 { 
80 public: 
81 // @tacent Added the () to initialize when new called. 
82 array2d(int w, int h) { width = w; height = h; data = new T[width * height](); } 
83 
84 array2d(const array2d<T>& rhs) { width = rhs.width; height = rhs.height; data = new T[width * height]; for (int i=0; i<width; i++) for (int j=0; j<height; j++) (*this)(i,j) = rhs.data[j*width + i]; } 
85 ~array2d() { delete [] data; } 
86 
87 T& operator()(int col, int row) { return data[row*width + col]; } 
88 int get_width() { return width; } 
89 int get_height() { return height; } 
90 
91 array2d<T>& operator*=(T scalar) { for (int i=0; i<width; i++) for (int j=0; j<height; j++) (*this)(i,j) *= scalar; return *this; } 
92 array2d<T> operator*(T scalar) { array2d<T> result(*this); result *= scalar; return result; } 
93 std::vector<T> operator*(std::vector<T> vec) { std::vector<T> result; T sum; for(int row=0; row<get_height(); row++) { sum = 0; for (int col=0; col<get_width(); col++) sum += (*this)(col,row) * vec[col]; result.push_back(sum); } return result; } 
94 array2d<T>& multiply_row_scalar(int row, double mult) { for (int i=0; i<get_width(); i++) { (*this)(i,row) *= mult; } return *this; } 
95 array2d<T>& add_row_multiple(int from_row, int to_row, double mult) { for (int i=0; i<get_width(); i++) { (*this)(i,to_row) += mult*(*this)(i,from_row); } return *this; } 
96 
97 // We use simple Gaussian elimination - perf doesn't matter since the matrices will be K x K, where K = number of palette entries. 
98 array2d<T> matrix_inverse() 
99 { 
100 array2d<T> result(get_width(), get_height()); 
101 array2d<T>& a = *this; 
102 
103 // Set result to identity matrix 
104 result *= 0; 
105 for (int i=0; i<get_width(); i++) 
106 result(i,i) = 1; 
107 
108 // Reduce to echelon form, mirroring in result 
109 for (int i=0; i<get_width(); i++) 
110 { 
111 result.multiply_row_scalar(i, 1/a(i,i)); 
112 multiply_row_scalar(row: i, mult: 1/a(i,i)); 
113 for (int j=i+1; j<get_height(); j++) 
114 { 
115 result.add_row_multiple(i, j, -a(i,j)); 
116 add_row_multiple(from_row: i, to_row: j, mult: -a(i,j)); 
117 } 
118 } 
119 // Back substitute, mirroring in result 
120 for (int i=get_width()-1; i>=0; i--) 
121 { 
122 for (int j=i-1; j>=0; j--) 
123 { 
124 result.add_row_multiple(i, j, -a(i,j)); 
125 add_row_multiple(from_row: i, to_row: j, mult: -a(i,j)); 
126 } 
127 } 
128 
129 // result is now the inverse 
130 return result; 
131 } 
132 
133 private: 
134 T* data; 
135 int width, height; 
136 }; 
137 
138 template <typename T> array2d<T> operator*(T scalar, array2d<T> a) { return a*scalar; } 
139 
140 template <typename T> class array3d 
141 { 
142 public: 
143 // @tacent Added the () to initialize when new called. 
144 array3d(int w, int h, int d) { width = w; height = h; depth = d; data = new T[width * height * depth](); } 
145 
146 array3d(const array3d<T>& rhs) { width = rhs.width; height = rhs.height; depth = rhs.depth; data = new T[width * height * depth]; for (int i=0; i<width; i++) for (int j=0; j<height; j++) for (int k=0; k<depth; k++) (*this)(i,j,k) = rhs.data[j*width*depth + i*depth + k]; } 
147 ~array3d() { delete [] data; } 
148 
149 T& operator()(int col, int row, int layer) { return data[row*width*depth + col*depth + layer]; } 
150 
151 int get_width() { return width; } 
152 int get_height() { return height; } 
153 int get_depth() { return depth; } 
154 
155 private: 
156 T* data; 
157 int width, height, depth; 
158 }; 
159 
160 int compute_max_coarse_level(int width, int height); 
161 void fill_random(array3d<double>& a, tMath::tRandom::tGenerator& gen); 
162 double get_initial_temperature(); 
163 double get_final_temperature(); 
164 void random_permutation(int count, std::vector<int>& result, std::default_random_engine& randEng); 
165 void random_permutation_2d(int width, int height, std::deque< std::pair<int, int> >& result, std::default_random_engine& randEng); 
166 void compute_b_array(array2d< vector_fixed<double, 3> >& filter_weights, array2d< vector_fixed<double, 3> >& b); 
167 vector_fixed<double, 3> b_value(array2d< vector_fixed<double, 3> >& b, int i_x, int i_y, int j_x, int j_y); 
168 void compute_a_image(array2d< vector_fixed<double, 3> >& image, array2d< vector_fixed<double, 3> >& b, array2d< vector_fixed<double, 3> >& a); 
169 void sum_coarsen(array2d< vector_fixed<double, 3> >& fine, array2d< vector_fixed<double, 3> >& coarse); 
170 template <typename T, int length> array2d<T> extract_vector_layer_2d(array2d< vector_fixed<T, length> > s, int k); 
171 template <typename T, int length> std::vector<T> extract_vector_layer_1d(std::vector< vector_fixed<T, length> > s, int k); 
172 int best_match_color(array3d<double>& vars, int i_x, int i_y, std::vector< vector_fixed<double, 3> >& palette); 
173 void zoom_double(array3d<double>& small, array3d<double>& big); 
174 void compute_initial_s(array2d< vector_fixed<double,3> >& s, array3d<double>& coarse_variables, array2d< vector_fixed<double, 3> >& b); 
175 void update_s(array2d< vector_fixed<double,3> >& s, array3d<double>& coarse_variables, array2d< vector_fixed<double, 3> >& b, int j_x, int j_y, int alpha, double delta); 
176 void refine_palette(array2d< vector_fixed<double,3> >& s, array3d<double>& coarse_variables, array2d< vector_fixed<double, 3> >& a, std::vector< vector_fixed<double, 3> >& palette); 
177 void compute_initial_j_palette_sum(array2d< vector_fixed<double, 3> >& j_palette_sum, array3d<double>& coarse_variables, std::vector< vector_fixed<double, 3> >& palette); 
178 bool spatial_color_quant 
179 ( 
180 array2d< vector_fixed<double, 3> >& image, array2d< vector_fixed<double, 3> >& filter_weights, array2d< int >& quantized_image, 
181 std::vector< vector_fixed<double, 3> >& palette, array3d<double>*& p_coarse_variables, 
182 double initial_temperature, double final_temperature, int temps_per_level, int repeats_per_temp, 
183 tMath::tRandom::tGenerator& randGen, std::default_random_engine& randEng 
184 ); 
185} 
186 
187 
188int tQuantizeSpatial::compute_max_coarse_level(int width, int height) 
189{ 
190 // We want the coarsest layer to have at most MAX_PIXELS pixels 
191 const int MAX_PIXELS = 4000; 
192 int result = 0; 
193 while (width * height > MAX_PIXELS) 
194 { 
195 width >>= 1; 
196 height >>= 1; 
197 result++; 
198 } 
199 return result; 
200} 
201 
202 
203void tQuantizeSpatial::fill_random(array3d<double>& a, tMath::tRandom::tGenerator& gen) 
204{ 
205 for (int i=0; i<a.get_width(); i++) 
206 for (int j=0; j<a.get_height(); j++) 
207 for (int k=0; k<a.get_depth(); k++) 
208 //a(i,j,k) = ((double)rand())/RAND_MAX; 
209 a(i,j,k) = tMath::tRandom::tGetDouble(gen); 
210} 
211 
212 
213double tQuantizeSpatial::get_initial_temperature() 
214{ 
215 return 2.0; // TODO: Figure out what to make this 
216} 
217 
218 
219double tQuantizeSpatial::get_final_temperature() 
220{ 
221 return 0.02; // TODO: Figure out what to make this 
222} 
223 
224 
225void tQuantizeSpatial::random_permutation(int count, std::vector<int>& result, std::default_random_engine& randEng) 
226{ 
227 result.clear(); 
228 for(int i=0; i<count; i++) 
229 result.push_back(x: i); 
230 
231 std::shuffle(first: result.begin(), last: result.end(), g&: randEng); 
232} 
233 
234 
235void tQuantizeSpatial::random_permutation_2d(int width, int height, std::deque< std::pair<int, int> >& result, std::default_random_engine& randEng) 
236{ 
237 std::vector<int> perm1d; 
238 random_permutation(count: width*height, result&: perm1d, randEng); 
239 while(!perm1d.empty()) 
240 { 
241 int idx = perm1d.back(); 
242 perm1d.pop_back(); 
243 result.push_back(x: std::pair<int,int>(idx % width, idx / width)); 
244 } 
245} 
246 
247 
248void tQuantizeSpatial::compute_b_array(array2d< vector_fixed<double, 3> >& filter_weights, array2d< vector_fixed<double, 3> >& b) 
249{ 
250 // Assume that the pixel i is always located at the center of b, and vary pixel j's location through each location in b. 
251 int radius_width = (filter_weights.get_width() - 1)/2, radius_height = (filter_weights.get_height() - 1)/2; 
252 int offset_x = (b.get_width() - 1)/2 - radius_width; 
253 int offset_y = (b.get_height() - 1)/2 - radius_height; 
254 for (int j_y = 0; j_y < b.get_height(); j_y++) 
255 { 
256 for (int j_x = 0; j_x < b.get_width(); j_x++) 
257 { 
258 for (int k_y = 0; k_y < filter_weights.get_height(); k_y++) 
259 { 
260 for (int k_x = 0; k_x < filter_weights.get_width(); k_x++) 
261 { 
262 if (k_x+offset_x >= j_x - radius_width && 
263 k_x+offset_x <= j_x + radius_width && 
264 k_y+offset_y >= j_y - radius_width && 
265 k_y+offset_y <= j_y + radius_width) 
266 { 
267 b(j_x,j_y) += filter_weights(k_x,k_y).direct_product(rhs&: filter_weights(k_x+offset_x-j_x+radius_width,k_y+offset_y-j_y+radius_height)); 
268 } 
269 } 
270 }  
271 } 
272 } 
273} 
274 
275 
276tQuantizeSpatial::vector_fixed<double, 3> tQuantizeSpatial::b_value(array2d< vector_fixed<double, 3> >& b, int i_x, int i_y, int j_x, int j_y) 
277{ 
278 int radius_width = (b.get_width() - 1)/2, radius_height = (b.get_height() - 1)/2; 
279 int k_x = j_x - i_x + radius_width; 
280 int k_y = j_y - i_y + radius_height; 
281 if (k_x >= 0 && k_y >= 0 && k_x < b.get_width() && k_y < b.get_height()) 
282 return b(k_x, k_y); 
283 else 
284 return vector_fixed<double, 3>(); 
285} 
286 
287 
288void tQuantizeSpatial::compute_a_image(array2d< vector_fixed<double, 3> >& image, array2d< vector_fixed<double, 3> >& b, array2d< vector_fixed<double, 3> >& a) 
289{ 
290 int radius_width = (b.get_width() - 1)/2, radius_height = (b.get_height() - 1)/2; 
291 for (int i_y = 0; i_y < a.get_height(); i_y++) 
292 { 
293 for (int i_x = 0; i_x < a.get_width(); i_x++) 
294 { 
295 for (int j_y = i_y - radius_height; j_y <= i_y + radius_height; j_y++) 
296 { 
297 if (j_y < 0) j_y = 0; 
298 if (j_y >= a.get_height()) break; 
299 
300 for (int j_x = i_x - radius_width; j_x <= i_x + radius_width; j_x++) 
301 { 
302 if (j_x < 0) j_x = 0; 
303 if (j_x >= a.get_width()) break; 
304 
305 a(i_x,i_y) += b_value(b, i_x, i_y, j_x, j_y).direct_product(rhs&: image(j_x,j_y)); 
306 } 
307 } 
308 a(i_x, i_y) *= -2.0; 
309 } 
310 } 
311} 
312 
313 
314void tQuantizeSpatial::sum_coarsen(array2d< vector_fixed<double, 3> >& fine, array2d< vector_fixed<double, 3> >& coarse) 
315{ 
316 for (int y=0; y<coarse.get_height(); y++) 
317 { 
318 for (int x=0; x<coarse.get_width(); x++) 
319 { 
320 double divisor = 1.0; 
321 vector_fixed<double, 3> val = fine(x*2, y*2); 
322 if (x*2 + 1 < fine.get_width()) 
323 { 
324 divisor += 1; val += fine(x*2 + 1, y*2); 
325 } 
326 if (y*2 + 1 < fine.get_height()) 
327 { 
328 divisor += 1; val += fine(x*2, y*2 + 1); 
329 } 
330 if (x*2 + 1 < fine.get_width() && y*2 + 1 < fine.get_height()) 
331 { 
332 divisor += 1; val += fine(x*2 + 1, y*2 + 1); 
333 } 
334 coarse(x, y) = /*(1/divisor)**/val; 
335 } 
336 } 
337} 
338 
339 
340template <typename T, int length> tQuantizeSpatial::array2d<T> tQuantizeSpatial::extract_vector_layer_2d(array2d< vector_fixed<T, length> > s, int k) 
341{ 
342 array2d<T> result(s.get_width(), s.get_height()); 
343 for (int i=0; i < s.get_width(); i++) 
344 { 
345 for (int j=0; j < s.get_height(); j++) 
346 { 
347 result(i,j) = s(i,j)(k); 
348 } 
349 } 
350 return result; 
351} 
352 
353 
354template <typename T, int length> std::vector<T> tQuantizeSpatial::extract_vector_layer_1d(std::vector< vector_fixed<T, length> > s, int k) 
355{ 
356 std::vector<T> result; 
357 for (unsigned int i=0; i < s.size(); i++) 
358 { 
359 result.push_back(s[i](k)); 
360 } 
361 return result; 
362} 
363 
364 
365int tQuantizeSpatial::best_match_color(array3d<double>& vars, int i_x, int i_y, std::vector< vector_fixed<double, 3> >& palette) 
366{ 
367 int max_v = 0; 
368 double max_weight = vars(i_x, i_y, 0); 
369 for (unsigned int v=1; v < palette.size(); v++) 
370 { 
371 if (vars(i_x, i_y, v) > max_weight) 
372 { 
373 max_v = v; 
374 max_weight = vars(i_x, i_y, v); 
375 } 
376 } 
377 return max_v; 
378} 
379 
380 
381void tQuantizeSpatial::zoom_double(array3d<double>& small, array3d<double>& big) 
382{ 
383 // Simple scaling of the weights array based on mixing the four pixels falling under each fine pixel, weighted by area. 
384 // To mix the pixels a little, we assume each fine pixel is 1.2 fine pixels wide and high. 
385 for (int y=0; y<big.get_height()/2*2; y++) 
386 { 
387 for (int x=0; x<big.get_width()/2*2; x++) 
388 { 
389 double left = std::max(a: 0.0, b: (x-0.1)/2.0), right = std::min(a: small.get_width()-0.001, b: (x+1.1)/2.0); 
390 double top = std::max(a: 0.0, b: (y-0.1)/2.0), bottom = std::min(a: small.get_height()-0.001, b: (y+1.1)/2.0); 
391 int x_left = (int)floor(x: left), x_right = (int)floor(x: right); 
392 int y_top = (int)floor(x: top), y_bottom = (int)floor(x: bottom); 
393 double area = (right-left)*(bottom-top); 
394 double top_left_weight = (ceil(x: left) - left)*(ceil(x: top) - top)/area; 
395 double top_right_weight = (right - floor(x: right))*(ceil(x: top) - top)/area; 
396 double bottom_left_weight = (ceil(x: left) - left)*(bottom - floor(x: bottom))/area; 
397 double bottom_right_weight = (right - floor(x: right))*(bottom - floor(x: bottom))/area; 
398 double top_weight = (right-left)*(ceil(x: top) - top)/area; 
399 double bottom_weight = (right-left)*(bottom - floor(x: bottom))/area; 
400 double left_weight = (bottom-top)*(ceil(x: left) - left)/area; 
401 double right_weight = (bottom-top)*(right - floor(x: right))/area; 
402 for(int z=0; z<big.get_depth(); z++) 
403 { 
404 if (x_left == x_right && y_top == y_bottom) 
405 { 
406 big(x, y, z) = small(x_left, y_top, z); 
407 } 
408 else if (x_left == x_right) 
409 { 
410 big(x, y, z) = top_weight*small(x_left, y_top, z) + bottom_weight*small(x_left, y_bottom, z); 
411 } 
412 else if (y_top == y_bottom) 
413 { 
414 big(x, y, z) = left_weight*small(x_left, y_top, z) + right_weight*small(x_right, y_top, z); 
415 } 
416 else 
417 { 
418 big(x, y, z) = top_left_weight*small(x_left, y_top, z) + top_right_weight*small(x_right, y_top, z) + bottom_left_weight*small(x_left, y_bottom, z) + bottom_right_weight*small(x_right, y_bottom, z); 
419 } 
420 } 
421 } 
422 } 
423} 
424 
425 
426void tQuantizeSpatial::compute_initial_s(array2d< vector_fixed<double,3> >& s, array3d<double>& coarse_variables, array2d< vector_fixed<double, 3> >& b) 
427{ 
428 int palette_size = s.get_width(); 
429 int coarse_width = coarse_variables.get_width(); 
430 int coarse_height = coarse_variables.get_height(); 
431 int center_x = (b.get_width()-1)/2, center_y = (b.get_height()-1)/2; 
432 vector_fixed<double,3> center_b = b_value(b,i_x: 0,i_y: 0,j_x: 0,j_y: 0); 
433 vector_fixed<double,3> zero_vector; 
434 for (int v=0; v<palette_size; v++) 
435 { 
436 for (int alpha=v; alpha<palette_size; alpha++) 
437 { 
438 s(v,alpha) = zero_vector; 
439 } 
440 } 
441 for (int i_y=0; i_y<coarse_height; i_y++) 
442 { 
443 for (int i_x=0; i_x<coarse_width; i_x++) 
444 { 
445 int max_j_x = std::min(a: coarse_width, b: i_x - center_x + b.get_width()); 
446 int max_j_y = std::min(a: coarse_height, b: i_y - center_y + b.get_height()); 
447 for (int j_y=std::max(a: 0, b: i_y - center_y); j_y<max_j_y; j_y++) 
448 { 
449 for (int j_x=std::max(a: 0, b: i_x - center_x); j_x<max_j_x; j_x++) 
450 { 
451 if (i_x == j_x && i_y == j_y) continue; 
452 vector_fixed<double,3> b_ij = b_value(b,i_x,i_y,j_x,j_y); 
453 for (int v=0; v<palette_size; v++) 
454 { 
455 for (int alpha=v; alpha<palette_size; alpha++) 
456 { 
457 double mult = coarse_variables(i_x,i_y,v)*coarse_variables(j_x,j_y,alpha); 
458 s(v,alpha)(0) += mult * b_ij(0); 
459 s(v,alpha)(1) += mult * b_ij(1); 
460 s(v,alpha)(2) += mult * b_ij(2); 
461 } 
462 } 
463 } 
464 }  
465 for (int v=0; v<palette_size; v++) 
466 { 
467 s(v,v) += coarse_variables(i_x,i_y,v)*center_b; 
468 } 
469 } 
470 } 
471} 
472 
473 
474void tQuantizeSpatial::update_s(array2d< vector_fixed<double,3> >& s, array3d<double>& coarse_variables, array2d< vector_fixed<double, 3> >& b, int j_x, int j_y, int alpha, double delta) 
475{ 
476 int palette_size = s.get_width(); 
477 int coarse_width = coarse_variables.get_width(); 
478 int coarse_height = coarse_variables.get_height(); 
479 int center_x = (b.get_width()-1)/2, center_y = (b.get_height()-1)/2; 
480 int max_i_x = std::min(a: coarse_width, b: j_x + center_x + 1); 
481 int max_i_y = std::min(a: coarse_height, b: j_y + center_y + 1); 
482 for (int i_y=std::max(a: 0, b: j_y - center_y); i_y<max_i_y; i_y++) 
483 { 
484 for (int i_x=std::max(a: 0, b: j_x - center_x); i_x<max_i_x; i_x++) 
485 { 
486 vector_fixed<double,3> delta_b_ij = delta*b_value(b,i_x,i_y,j_x,j_y); 
487 if (i_x == j_x && i_y == j_y) continue; 
488 for (int v=0; v <= alpha; v++) 
489 { 
490 double mult = coarse_variables(i_x,i_y,v); 
491 s(v,alpha)(0) += mult * delta_b_ij(0); 
492 s(v,alpha)(1) += mult * delta_b_ij(1); 
493 s(v,alpha)(2) += mult * delta_b_ij(2); 
494 } 
495 for (int v=alpha; v<palette_size; v++) 
496 { 
497 double mult = coarse_variables(i_x,i_y,v); 
498 s(alpha,v)(0) += mult * delta_b_ij(0); 
499 s(alpha,v)(1) += mult * delta_b_ij(1); 
500 s(alpha,v)(2) += mult * delta_b_ij(2); 
501 } 
502 } 
503 } 
504 s(alpha,alpha) += delta*b_value(b,i_x: 0,i_y: 0,j_x: 0,j_y: 0); 
505} 
506 
507 
508void tQuantizeSpatial::refine_palette(array2d< vector_fixed<double,3> >& s, array3d<double>& coarse_variables, array2d< vector_fixed<double, 3> >& a, std::vector< vector_fixed<double, 3> >& palette) 
509{ 
510 // We only computed the half of S above the diagonal - reflect it 
511 for (int v=0; v<s.get_width(); v++) 
512 { 
513 for (int alpha=0; alpha<v; alpha++) 
514 { 
515 s(v,alpha) = s(alpha,v); 
516 } 
517 } 
518 
519 std::vector< vector_fixed<double,3> > r(palette.size()); 
520 for (unsigned int v=0; v<palette.size(); v++) 
521 { 
522 for (int i_y=0; i_y<coarse_variables.get_height(); i_y++) 
523 { 
524 for (int i_x=0; i_x<coarse_variables.get_width(); i_x++) 
525 { 
526 r[v] += coarse_variables(i_x,i_y,v)*a(i_x,i_y); 
527 } 
528 } 
529 } 
530 
531 for (unsigned int k=0; k<3; k++) 
532 { 
533 array2d<double> S_k = extract_vector_layer_2d(s, k); 
534 std::vector<double> R_k = extract_vector_layer_1d(s: r, k); 
535 std::vector<double> palette_channel = -1.0*((2.0*S_k).matrix_inverse())*R_k; 
536 for (unsigned int v=0; v<palette.size(); v++) 
537 { 
538 double val = palette_channel[v]; 
539 if (val < 0) val = 0; 
540 if (val > 1) val = 1; 
541 palette[v](k) = val; 
542 }  
543 } 
544} 
545 
546 
547void tQuantizeSpatial::compute_initial_j_palette_sum(array2d< vector_fixed<double, 3> >& j_palette_sum, array3d<double>& coarse_variables, std::vector< vector_fixed<double, 3> >& palette) 
548{ 
549 for (int j_y=0; j_y<coarse_variables.get_height(); j_y++) 
550 { 
551 for (int j_x=0; j_x<coarse_variables.get_width(); j_x++) 
552 { 
553 vector_fixed<double, 3> palette_sum = vector_fixed<double, 3>(); 
554 for (unsigned int alpha=0; alpha < palette.size(); alpha++) 
555 { 
556 palette_sum += coarse_variables(j_x,j_y,alpha)*palette[alpha]; 
557 } 
558 j_palette_sum(j_x, j_y) = palette_sum; 
559 } 
560 } 
561} 
562 
563 
564bool tQuantizeSpatial::spatial_color_quant 
565( 
566 array2d< vector_fixed<double, 3> >& image, array2d< vector_fixed<double, 3> >& filter_weights, array2d< int >& quantized_image, 
567 std::vector< vector_fixed<double, 3> >& palette, array3d<double>*& p_coarse_variables, 
568 double initial_temperature, double final_temperature, int temps_per_level, int repeats_per_temp, 
569 tMath::tRandom::tGenerator& randGen, std::default_random_engine& randEng 
570) 
571{ 
572 int max_coarse_level = compute_max_coarse_level(width: image.get_width(), height: image.get_height()); //1; 
573 p_coarse_variables = new array3d<double>(image.get_width() >> max_coarse_level, image.get_height() >> max_coarse_level, palette.size()); 
574 // For syntactic convenience 
575 array3d<double>& coarse_variables = *p_coarse_variables; 
576 fill_random(a&: coarse_variables, gen&: randGen); 
577 
578 double temperature = initial_temperature; 
579 
580 // Compute a_i, b_{ij} according to (11) 
581 int extended_neighborhood_width = filter_weights.get_width()*2 - 1; 
582 int extended_neighborhood_height = filter_weights.get_height()*2 - 1; 
583 array2d< vector_fixed<double, 3> > b0(extended_neighborhood_width, extended_neighborhood_height); 
584 compute_b_array(filter_weights, b&: b0); 
585 
586 array2d< vector_fixed<double, 3> > a0(image.get_width(), image.get_height()); 
587 compute_a_image(image, b&: b0, a&: a0); 
588 
589 // Compute a_I^l, b_{IJ}^l according to (18) 
590 std::vector< array2d< vector_fixed<double, 3> > > a_vec, b_vec; 
591 a_vec.push_back(x: a0); 
592 b_vec.push_back(x: b0); 
593 
594 int coarse_level; 
595 for (coarse_level=1; coarse_level <= max_coarse_level; coarse_level++) 
596 { 
597 int radius_width = (filter_weights.get_width() - 1)/2, radius_height = (filter_weights.get_height() - 1)/2; 
598 array2d< vector_fixed<double, 3> > 
599 bi(std::max(a: 3, b: b_vec.back().get_width()-2), std::max(a: 3, b: b_vec.back().get_height()-2)); 
600 for (int J_y=0; J_y<bi.get_height(); J_y++) 
601 { 
602 for (int J_x=0; J_x<bi.get_width(); J_x++) 
603 { 
604 for (int i_y=radius_height*2; i_y<radius_height*2+2; i_y++) 
605 { 
606 for (int i_x=radius_width*2; i_x<radius_width*2+2; i_x++) 
607 { 
608 for (int j_y=J_y*2; j_y<J_y*2+2; j_y++) 
609 { 
610 for (int j_x=J_x*2; j_x<J_x*2+2; j_x++) 
611 { 
612 bi(J_x,J_y) += b_value(b&: b_vec.back(), i_x, i_y, j_x, j_y); 
613 } 
614 } 
615 } 
616 } 
617 } 
618 } 
619 b_vec.push_back(x: bi); 
620 array2d< vector_fixed<double, 3> > ai(image.get_width() >> coarse_level, image.get_height() >> coarse_level); 
621 sum_coarsen(fine&: a_vec.back(), coarse&: ai); 
622 a_vec.push_back(x: ai); 
623 } 
624 
625 // Multiscale annealing 
626 coarse_level = max_coarse_level; 
627 const int iters_per_level = temps_per_level; 
628 double temperature_multiplier = std::pow(x: final_temperature/initial_temperature, y: 1.0/(std::max(a: 3, b: max_coarse_level*iters_per_level))); 
629 
630 int iters_at_current_level = 0; 
631 bool skip_palette_maintenance = false; 
632 array2d< vector_fixed<double,3> > s(palette.size(), palette.size()); 
633 compute_initial_s(s, coarse_variables, b&: b_vec[coarse_level]); 
634 array2d< vector_fixed<double, 3> >* j_palette_sum = 
635 new array2d< vector_fixed<double, 3> >(coarse_variables.get_width(), coarse_variables.get_height()); 
636 compute_initial_j_palette_sum(j_palette_sum&: *j_palette_sum, coarse_variables, palette); 
637 while (coarse_level >= 0 || temperature > final_temperature) 
638 { 
639 // Need to reseat this reference in case we changed p_coarse_variables 
640 array3d<double>& coarse_variables = *p_coarse_variables; 
641 array2d< vector_fixed<double, 3> >& a = a_vec[coarse_level]; 
642 array2d< vector_fixed<double, 3> >& b = b_vec[coarse_level]; 
643 vector_fixed<double,3> middle_b = b_value(b,i_x: 0,i_y: 0,j_x: 0,j_y: 0); 
644 int center_x = (b.get_width()-1)/2, center_y = (b.get_height()-1)/2; 
645 int step_counter = 0; 
646 for (int repeat=0; repeat<repeats_per_temp; repeat++) 
647 { 
648 int pixels_changed = 0, pixels_visited = 0; 
649 std::deque< std::pair<int, int> > visit_queue; 
650 random_permutation_2d(width: coarse_variables.get_width(), height: coarse_variables.get_height(), result&: visit_queue, randEng); 
651 
652 // Compute 2*sum(j in extended neighborhood of i, j != i) b_ij 
653 while (!visit_queue.empty()) 
654 { 
655 // If we get to 10% above initial size, just revisit them all 
656 if ((int)visit_queue.size() > coarse_variables.get_width()*coarse_variables.get_height()*11/10) 
657 { 
658 visit_queue.clear(); 
659 random_permutation_2d(width: coarse_variables.get_width(), height: coarse_variables.get_height(), result&: visit_queue, randEng); 
660 } 
661 
662 int i_x = visit_queue.front().first, i_y = visit_queue.front().second; 
663 visit_queue.pop_front(); 
664 
665 // Compute (25) 
666 vector_fixed<double,3> p_i; 
667 for (int y=0; y<b.get_height(); y++) 
668 { 
669 for (int x=0; x<b.get_width(); x++) 
670 { 
671 int j_x = x - center_x + i_x, j_y = y - center_y + i_y; 
672 if (i_x == j_x && i_y == j_y) continue; 
673 if (j_x < 0 || j_y < 0 || j_x >= coarse_variables.get_width() || j_y >= coarse_variables.get_height()) continue; 
674 vector_fixed<double,3> b_ij = b_value(b, i_x, i_y, j_x, j_y); 
675 vector_fixed<double,3> j_pal = (*j_palette_sum)(j_x,j_y); 
676 p_i(0) += b_ij(0)*j_pal(0); 
677 p_i(1) += b_ij(1)*j_pal(1); 
678 p_i(2) += b_ij(2)*j_pal(2); 
679 } 
680 } 
681 p_i *= 2.0; 
682 p_i += a(i_x, i_y); 
683 
684 std::vector<double> meanfield_logs, meanfields; 
685 double max_meanfield_log = -std::numeric_limits<double>::infinity(); 
686 double meanfield_sum = 0.0; 
687 for (unsigned int v=0; v < palette.size(); v++) 
688 { 
689 // Update m_{pi(i)v}^I according to (23). We can subtract an arbitrary factor to prevent overflow, 
690 // since only the weight relative to the sum matters, so we will choose a value that makes the maximum e^100. 
691 meanfield_logs.push_back(x: -(palette[v].dot_product( 
692 rhs: p_i + middle_b.direct_product(rhs&: palette[v])))/temperature); 
693 if (meanfield_logs.back() > max_meanfield_log) 
694 { 
695 max_meanfield_log = meanfield_logs.back(); 
696 } 
697 } 
698 
699 for (unsigned int v=0; v < palette.size(); v++) 
700 { 
701 meanfields.push_back(x: exp(x: meanfield_logs[v]-max_meanfield_log+100)); 
702 meanfield_sum += meanfields.back(); 
703 } 
704 
705 
706 if (meanfield_sum == 0.0) 
707 return false; 
708 int old_max_v = best_match_color(vars&: coarse_variables, i_x, i_y, palette); 
709 vector_fixed<double,3> & j_pal = (*j_palette_sum)(i_x,i_y); 
710 for (unsigned int v=0; v < palette.size(); v++) 
711 { 
712 double new_val = meanfields[v]/meanfield_sum; 
713 // Prevent the matrix S from becoming singular 
714 if (new_val <= 0) new_val = 1e-10; 
715 if (new_val >= 1) new_val = 1 - 1e-10; 
716 double delta_m_iv = new_val - coarse_variables(i_x,i_y,v); 
717 coarse_variables(i_x,i_y,v) = new_val; 
718 j_pal(0) += delta_m_iv*palette[v](0); 
719 j_pal(1) += delta_m_iv*palette[v](1); 
720 j_pal(2) += delta_m_iv*palette[v](2); 
721 if (abs(x: delta_m_iv) > 0.001 && !skip_palette_maintenance) 
722 update_s(s, coarse_variables, b, j_x: i_x, j_y: i_y, alpha: v, delta: delta_m_iv); 
723 } 
724 int max_v = best_match_color(vars&: coarse_variables, i_x, i_y, palette); 
725 // Only consider it a change if the colors are different enough 
726 if ((palette[max_v]-palette[old_max_v]).norm_squared() >= 1.0/(255.0*255.0)) 
727 { 
728 pixels_changed++; 
729 // We don't add the outer layer of pixels , because there isn't much weight there, and if it does need 
730 // to be visited, it'll probably be added when we visit neighboring pixels. The commented out loops are 
731 // faster but cause a little bit of distortion 
732 //for (int y=center_y-1; y<center_y+1; y++) { 
733 // for (int x=center_x-1; x<center_x+1; x++) { 
734 for (int y=std::min(a: 1,b: center_y-1); y<std::max(a: b.get_height()-1,b: center_y+1); y++) 
735 { 
736 for (int x=std::min(a: 1,b: center_x-1); x<std::max(a: b.get_width()-1,b: center_x+1); x++) 
737 { 
738 int j_x = x - center_x + i_x, j_y = y - center_y + i_y; 
739 if (j_x < 0 || j_y < 0 || j_x >= coarse_variables.get_width() || j_y >= coarse_variables.get_height()) continue; 
740 visit_queue.push_back(x: std::pair<int,int>(j_x,j_y)); 
741 } 
742 } 
743 } 
744 pixels_visited++; 
745 
746 // Show progress with dots - in a graphical interface, we'd show progressive refinements of the image instead, and maybe a palette preview. 
747 #if 0 
748 step_counter++; 
749 if ((step_counter % 10000) == 0) 
750 { 
751 cout << "."; 
752 cout.flush(); 
753 } 
754 #endif 
755 } 
756 
757 if (skip_palette_maintenance) 
758 { 
759 compute_initial_s(s, coarse_variables&: *p_coarse_variables, b&: b_vec[coarse_level]); 
760 } 
761 refine_palette(s, coarse_variables, a, palette); 
762 compute_initial_j_palette_sum(j_palette_sum&: *j_palette_sum, coarse_variables, palette); 
763 } 
764 
765 iters_at_current_level++; 
766 skip_palette_maintenance = false; 
767 if ((temperature <= final_temperature || coarse_level > 0) && iters_at_current_level >= iters_per_level) 
768 { 
769 coarse_level--; 
770 if (coarse_level < 0) break; 
771 array3d<double>* p_new_coarse_variables = new array3d<double>(image.get_width() >> coarse_level, image.get_height() >> coarse_level, palette.size()); 
772 zoom_double(small&: coarse_variables, big&: *p_new_coarse_variables); 
773 delete p_coarse_variables; 
774 p_coarse_variables = p_new_coarse_variables; 
775 iters_at_current_level = 0; 
776 delete j_palette_sum; 
777 j_palette_sum = new array2d< vector_fixed<double, 3> >((*p_coarse_variables).get_width(), (*p_coarse_variables).get_height()); 
778 compute_initial_j_palette_sum(j_palette_sum&: *j_palette_sum, coarse_variables&: *p_coarse_variables, palette); 
779 skip_palette_maintenance = true; 
780 } 
781 if (temperature > final_temperature) 
782 { 
783 temperature *= temperature_multiplier; 
784 } 
785 } 
786 
787 // This is normally not used, but is handy sometimes for debugging 
788 while (coarse_level > 0) 
789 { 
790 coarse_level--; 
791 array3d<double>* p_new_coarse_variables = new array3d<double>(image.get_width() >> coarse_level, image.get_height() >> coarse_level, palette.size()); 
792 zoom_double(small&: *p_coarse_variables, big&: *p_new_coarse_variables); 
793 delete p_coarse_variables; 
794 p_coarse_variables = p_new_coarse_variables; 
795 } 
796 
797 { 
798 // Need to reseat this reference in case we changed p_coarse_variables 
799 array3d<double>& coarse_variables = *p_coarse_variables; 
800 
801 for (int i_x = 0; i_x < image.get_width(); i_x++) 
802 { 
803 for(int i_y = 0; i_y < image.get_height(); i_y++) 
804 { 
805 quantized_image(i_x,i_y) = 
806 best_match_color(vars&: coarse_variables, i_x, i_y, palette); 
807 } 
808 } 
809 for (unsigned int v=0; v<palette.size(); v++) 
810 { 
811 for (unsigned int k=0; k<3; k++) 
812 { 
813 if (palette[v](k) > 1.0) palette[v](k) = 1.0; 
814 if (palette[v](k) < 0.0) palette[v](k) = 0.0; 
815 } 
816 } 
817 } 
818 
819 return true; 
820} 
821 
822 
823// 
824// The functions below make up the external interface. 
825// 
826 
827 
828double tQuantizeSpatial::ComputeBaseDither(int width, int height, int numColours) 
829{  
830 double ditherLevel = 0.09*log(x: double(width*height)) - 0.04*log(x: double(numColours)) + 0.001; 
831 if (ditherLevel < 0.0) 
832 ditherLevel = 0.001; 
833 
834 return ditherLevel; 
835} 
836 
837 
838bool tQuantizeSpatial::QuantizeImage 
839( 
840 int numColours, int width, int height, const tPixel3b* pixels, tColour3b* destPalette, uint8* destIndices, 
841 bool checkExact, double ditherLevel, int filterSize 
842) 
843{ 
844 if ((numColours < 2) || (numColours > 256) || (width <= 0) || (height <= 0) || !pixels || !destPalette || !destIndices) 
845 return false; 
846 
847 if ((filterSize != 1) && (filterSize != 3) && (filterSize != 5)) 
848 return false; 
849 
850 if (ditherLevel <= 0.0) 
851 ditherLevel = tQuantizeSpatial::ComputeBaseDither(width, height, numColours); 
852 
853 if (checkExact) 
854 { 
855 bool success = tQuantize::QuantizeImageExact(numColours, width, height, pixels, destPalette, destIndices); 
856 if (success) 
857 return true; 
858 } 
859 
860 // Seeding the generator with the same value every time guarantees repeatability. 
861 tMath::tRandom::tGeneratorMersenneTwister randGen(uint32(147)); 
862 std::default_random_engine randEng(137); 
863 
864 array2d< vector_fixed<double, 3> > image(width, height); 
865 array2d< vector_fixed<double, 3> > filter1_weights(1, 1); 
866 array2d< vector_fixed<double, 3> > filter3_weights(3, 3); 
867 array2d< vector_fixed<double, 3> > filter5_weights(5, 5); 
868 array2d< int > quantized_image(width, height); 
869 std::vector< vector_fixed<double, 3> > palette; 
870 
871 for(int k=0; k<3; k++) 
872 filter1_weights(0,0)(k) = 1.0; 
873 
874 for (int i = 0; i < numColours; i++) 
875 { 
876 vector_fixed<double, 3> v; 
877 v(0) = tMath::tRandom::tGetDouble(randGen); // ((double)rand())/RAND_MAX; 
878 v(1) = tMath::tRandom::tGetDouble(randGen); // ((double)rand())/RAND_MAX; 
879 v(2) = tMath::tRandom::tGetDouble(randGen); // ((double)rand())/RAND_MAX; 
880 palette.push_back(x: v); 
881 } 
882 
883 for (int y = 0; y < height; y++) 
884 { 
885 for (int x = 0; x < width; x++) 
886 { 
887 image(x,y)(0) = (pixels[x + y*width].R)/255.0; 
888 image(x,y)(1) = (pixels[x + y*width].G)/255.0; 
889 image(x,y)(2) = (pixels[x + y*width].B)/255.0; 
890 } 
891 }  
892 
893 array3d<double>* coarse_variables; 
894 double stddev = ditherLevel; 
895 double sum = 0.0; 
896 for(int i = 0; i < 3; i++) 
897 for(int j = 0; j < 3; j++) 
898 for(int k = 0; k < 3; k++) 
899 sum += filter3_weights(i,j)(k) = exp(x: -sqrt(x: (double)((i-1)*(i-1) + (j-1)*(j-1)))/(stddev*stddev)); 
900 
901 sum /= 3; 
902 for(int i = 0; i < 3; i++) 
903 for(int j = 0; j < 3; j++) 
904 for(int k = 0; k < 3; k++) 
905 filter3_weights(i,j)(k) /= sum; 
906 
907 sum = 0.0; 
908 for(int i = 0; i < 5; i++) 
909 for(int j = 0; j < 5; j++) 
910 for(int k = 0; k < 3; k++) 
911 sum += filter5_weights(i,j)(k) = exp(x: -sqrt(x: (double)((i-2)*(i-2) + (j-2)*(j-2)))/(stddev*stddev)); 
912 
913 sum /= 3; 
914 for(int i = 0; i < 5; i++) 
915 for(int j = 0; j < 5; j++) 
916 for(int k = 0; k < 3; k++) 
917 filter5_weights(i,j)(k) /= sum; 
918 
919 array2d< vector_fixed<double, 3> >* filters[] = {NULL, &filter1_weights, NULL, &filter3_weights, NULL, &filter5_weights}; 
920 bool success = spatial_color_quant(image, filter_weights&: *filters[filterSize], quantized_image, palette, p_coarse_variables&: coarse_variables, initial_temperature: 1.0, final_temperature: 0.001, temps_per_level: 3, repeats_per_temp: 1, randGen, randEng); 
921 // bool success = spatial_color_quant(image, filter3_weights, quantized_image, palette, coarse_variables, 0.05, 0.02); 
922 if (!success) 
923 return false; 
924 
925 for (int y = 0; y < height; y++) 
926 for (int x = 0; x < width; x++) 
927 destIndices[x + y*width] = uint8(quantized_image(x,y)); 
928 
929 for (int c = 0; c < numColours; c++) 
930 { 
931 destPalette[c].R = uint8(255.0 * palette[c](0)); 
932 destPalette[c].G = uint8(255.0 * palette[c](1)); 
933 destPalette[c].B = uint8(255.0 * palette[c](2)); 
934 } 
935 
936 return true; 
937} 
938 
939 
940bool tQuantizeSpatial::QuantizeImage 
941( 
942 int numColours, int width, int height, const tPixel4b* pixels, tColour3b* destPalette, uint8* destIndices, 
943 bool checkExact, double ditherLevel, int filterSize 
944) 
945{ 
946 if ((numColours < 2) || (numColours > 256) || (width <= 0) || (height <= 0) || !pixels || !destPalette || !destIndices) 
947 return false; 
948 
949 tPixel3b* pixels3 = new tPixel3b[width*height]; 
950 for (int p = 0; p < width*height; p++) 
951 pixels3[p].Set( r: pixels[p].R, g: pixels[p].G, b: pixels[p].B ); 
952 
953 bool success = QuantizeImage(numColours, width, height, pixels: pixels3, destPalette, destIndices, checkExact, ditherLevel, filterSize); 
954 delete[] pixels3; 
955 return success; 
956} 
957 
958 
959} 
960