15821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// Copyright 2011 Google Inc. All Rights Reserved. 25821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// 3eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch// Use of this source code is governed by a BSD-style license 4eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch// that can be found in the COPYING file in the root of the source 5eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch// tree. An additional intellectual property rights grant can be found 6eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch// in the file PATENTS. All contributing project authors may 7eb525c5499e34cc9c4b825d6d9e75bb07cc06aceBen Murdoch// be found in the AUTHORS file in the root of the source tree. 85821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// ----------------------------------------------------------------------------- 95821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// 105821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// Quantize levels for specified number of quantization-levels ([2, 256]). 115821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// Min and max values are preserved (usual 0 and 255 for alpha plane). 125821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// 135821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// Author: Skal (pascal.massimino@gmail.com) 145821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 155821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#include <assert.h> 165821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 175821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#include "./quant_levels.h" 185821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 195821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#if defined(__cplusplus) || defined(c_plusplus) 205821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)extern "C" { 215821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#endif 225821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 235821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#define NUM_SYMBOLS 256 245821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 255821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#define MAX_ITER 6 // Maximum number of convergence steps. 265821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#define ERROR_THRESHOLD 1e-4 // MSE stopping criterion. 275821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 285821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// ----------------------------------------------------------------------------- 295821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)// Quantize levels. 305821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 315821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)int QuantizeLevels(uint8_t* const data, int width, int height, 325821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) int num_levels, uint64_t* const sse) { 335821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) int freq[NUM_SYMBOLS] = { 0 }; 345821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) int q_level[NUM_SYMBOLS] = { 0 }; 355821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) double inv_q_level[NUM_SYMBOLS] = { 0 }; 365821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) int min_s = 255, max_s = 0; 375821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) const size_t data_size = height * width; 385821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) int i, num_levels_in, iter; 395821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) double last_err = 1.e38, err = 0.; 405821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) const double err_threshold = ERROR_THRESHOLD * data_size; 415821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 425821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (data == NULL) { 435821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) return 0; 445821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 455821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 465821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (width <= 0 || height <= 0) { 475821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) return 0; 485821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 495821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 505821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (num_levels < 2 || num_levels > 256) { 515821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) return 0; 525821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 535821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 545821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) { 555821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) size_t n; 565821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) num_levels_in = 0; 575821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (n = 0; n < data_size; ++n) { 585821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) num_levels_in += (freq[data[n]] == 0); 595821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (min_s > data[n]) min_s = data[n]; 605821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (max_s < data[n]) max_s = data[n]; 615821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) ++freq[data[n]]; 625821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 635821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 645821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 655821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (num_levels_in <= num_levels) goto End; // nothing to do! 665821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 675821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Start with uniformly spread centroids. 685821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (i = 0; i < num_levels; ++i) { 695821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) inv_q_level[i] = min_s + (double)(max_s - min_s) * i / (num_levels - 1); 705821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 715821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 725821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Fixed values. Won't be changed. 735821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) q_level[min_s] = 0; 745821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) q_level[max_s] = num_levels - 1; 755821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) assert(inv_q_level[0] == min_s); 765821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) assert(inv_q_level[num_levels - 1] == max_s); 775821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 785821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // k-Means iterations. 795821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (iter = 0; iter < MAX_ITER; ++iter) { 805821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) double q_sum[NUM_SYMBOLS] = { 0 }; 815821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) double q_count[NUM_SYMBOLS] = { 0 }; 825821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) int s, slot = 0; 835821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 845821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Assign classes to representatives. 855821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (s = min_s; s <= max_s; ++s) { 865821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Keep track of the nearest neighbour 'slot' 875821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) while (slot < num_levels - 1 && 885821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 2 * s > inv_q_level[slot] + inv_q_level[slot + 1]) { 895821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) ++slot; 905821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 915821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (freq[s] > 0) { 925821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) q_sum[slot] += s * freq[s]; 935821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) q_count[slot] += freq[s]; 945821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 955821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) q_level[s] = slot; 965821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 975821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 985821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Assign new representatives to classes. 995821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (num_levels > 2) { 1005821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (slot = 1; slot < num_levels - 1; ++slot) { 1015821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) const double count = q_count[slot]; 1025821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (count > 0.) { 1035821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) inv_q_level[slot] = q_sum[slot] / count; 1045821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1055821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1065821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1075821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1085821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Compute convergence error. 1095821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) err = 0.; 1105821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (s = min_s; s <= max_s; ++s) { 1115821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) const double error = s - inv_q_level[q_level[s]]; 1125821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) err += freq[s] * error * error; 1135821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1145821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1155821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Check for convergence: we stop as soon as the error is no 1165821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // longer improving. 1175821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (last_err - err < err_threshold) break; 1185821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) last_err = err; 1195821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1205821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1215821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Remap the alpha plane to quantized values. 1225821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) { 1235821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // double->int rounding operation can be costly, so we do it 1245821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // once for all before remapping. We also perform the data[] -> slot 1255821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // mapping, while at it (avoid one indirection in the final loop). 1265821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) uint8_t map[NUM_SYMBOLS]; 1275821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) int s; 1285821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) size_t n; 1295821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (s = min_s; s <= max_s; ++s) { 1305821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) const int slot = q_level[s]; 1315821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) map[s] = (uint8_t)(inv_q_level[slot] + .5); 1325821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1335821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Final pass. 1345821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (n = 0; n < data_size; ++n) { 1355821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) data[n] = map[data[n]]; 1365821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1375821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1385821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) End: 1395821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) // Store sum of squared error if needed. 1405821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (sse != NULL) *sse = (uint64_t)err; 1415821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1425821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) return 1; 1435821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)} 1445821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1455821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#if defined(__cplusplus) || defined(c_plusplus) 1465821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)} // extern "C" 1475821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#endif 148