1/* 2 * Copyright (C) 2007 Vitor Sessak <vitor1001@gmail.com> 3 * 4 * This file is part of FFmpeg. 5 * 6 * FFmpeg is free software; you can redistribute it and/or 7 * modify it under the terms of the GNU Lesser General Public 8 * License as published by the Free Software Foundation; either 9 * version 2.1 of the License, or (at your option) any later version. 10 * 11 * FFmpeg is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 14 * Lesser General Public License for more details. 15 * 16 * You should have received a copy of the GNU Lesser General Public 17 * License along with FFmpeg; if not, write to the Free Software 18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 19 */ 20 21/** 22 * @file 23 * Codebook Generator using the ELBG algorithm 24 */ 25 26#include <string.h> 27 28#include "libavutil/lfg.h" 29#include "elbg.h" 30#include "avcodec.h" 31 32#define DELTA_ERR_MAX 0.1 ///< Precision of the ELBG algorithm (as percentual error) 33 34/** 35 * In the ELBG jargon, a cell is the set of points that are closest to a 36 * codebook entry. Not to be confused with a RoQ Video cell. */ 37typedef struct cell_s { 38 int index; 39 struct cell_s *next; 40} cell; 41 42/** 43 * ELBG internal data 44 */ 45typedef struct{ 46 int error; 47 int dim; 48 int numCB; 49 int *codebook; 50 cell **cells; 51 int *utility; 52 int *utility_inc; 53 int *nearest_cb; 54 int *points; 55 AVLFG *rand_state; 56} elbg_data; 57 58static inline int distance_limited(int *a, int *b, int dim, int limit) 59{ 60 int i, dist=0; 61 for (i=0; i<dim; i++) { 62 dist += (a[i] - b[i])*(a[i] - b[i]); 63 if (dist > limit) 64 return INT_MAX; 65 } 66 67 return dist; 68} 69 70static inline void vect_division(int *res, int *vect, int div, int dim) 71{ 72 int i; 73 if (div > 1) 74 for (i=0; i<dim; i++) 75 res[i] = ROUNDED_DIV(vect[i],div); 76 else if (res != vect) 77 memcpy(res, vect, dim*sizeof(int)); 78 79} 80 81static int eval_error_cell(elbg_data *elbg, int *centroid, cell *cells) 82{ 83 int error=0; 84 for (; cells; cells=cells->next) 85 error += distance_limited(centroid, elbg->points + cells->index*elbg->dim, elbg->dim, INT_MAX); 86 87 return error; 88} 89 90static int get_closest_codebook(elbg_data *elbg, int index) 91{ 92 int i, pick=0, diff, diff_min = INT_MAX; 93 for (i=0; i<elbg->numCB; i++) 94 if (i != index) { 95 diff = distance_limited(elbg->codebook + i*elbg->dim, elbg->codebook + index*elbg->dim, elbg->dim, diff_min); 96 if (diff < diff_min) { 97 pick = i; 98 diff_min = diff; 99 } 100 } 101 return pick; 102} 103 104static int get_high_utility_cell(elbg_data *elbg) 105{ 106 int i=0; 107 /* Using linear search, do binary if it ever turns to be speed critical */ 108 int r = av_lfg_get(elbg->rand_state)%elbg->utility_inc[elbg->numCB-1] + 1; 109 while (elbg->utility_inc[i] < r) 110 i++; 111 112 assert(!elbg->cells[i]); 113 114 return i; 115} 116 117/** 118 * Implementation of the simple LBG algorithm for just two codebooks 119 */ 120static int simple_lbg(int dim, 121 int *centroid[3], 122 int newutility[3], 123 int *points, 124 cell *cells) 125{ 126 int i, idx; 127 int numpoints[2] = {0,0}; 128 int newcentroid[2][dim]; 129 cell *tempcell; 130 131 memset(newcentroid, 0, sizeof(newcentroid)); 132 133 newutility[0] = 134 newutility[1] = 0; 135 136 for (tempcell = cells; tempcell; tempcell=tempcell->next) { 137 idx = distance_limited(centroid[0], points + tempcell->index*dim, dim, INT_MAX)>= 138 distance_limited(centroid[1], points + tempcell->index*dim, dim, INT_MAX); 139 numpoints[idx]++; 140 for (i=0; i<dim; i++) 141 newcentroid[idx][i] += points[tempcell->index*dim + i]; 142 } 143 144 vect_division(centroid[0], newcentroid[0], numpoints[0], dim); 145 vect_division(centroid[1], newcentroid[1], numpoints[1], dim); 146 147 for (tempcell = cells; tempcell; tempcell=tempcell->next) { 148 int dist[2] = {distance_limited(centroid[0], points + tempcell->index*dim, dim, INT_MAX), 149 distance_limited(centroid[1], points + tempcell->index*dim, dim, INT_MAX)}; 150 int idx = dist[0] > dist[1]; 151 newutility[idx] += dist[idx]; 152 } 153 154 return newutility[0] + newutility[1]; 155} 156 157static void get_new_centroids(elbg_data *elbg, int huc, int *newcentroid_i, 158 int *newcentroid_p) 159{ 160 cell *tempcell; 161 int min[elbg->dim]; 162 int max[elbg->dim]; 163 int i; 164 165 for (i=0; i< elbg->dim; i++) { 166 min[i]=INT_MAX; 167 max[i]=0; 168 } 169 170 for (tempcell = elbg->cells[huc]; tempcell; tempcell = tempcell->next) 171 for(i=0; i<elbg->dim; i++) { 172 min[i]=FFMIN(min[i], elbg->points[tempcell->index*elbg->dim + i]); 173 max[i]=FFMAX(max[i], elbg->points[tempcell->index*elbg->dim + i]); 174 } 175 176 for (i=0; i<elbg->dim; i++) { 177 newcentroid_i[i] = min[i] + (max[i] - min[i])/3; 178 newcentroid_p[i] = min[i] + (2*(max[i] - min[i]))/3; 179 } 180} 181 182/** 183 * Add the points in the low utility cell to its closest cell. Split the high 184 * utility cell, putting the separed points in the (now empty) low utility 185 * cell. 186 * 187 * @param elbg Internal elbg data 188 * @param indexes {luc, huc, cluc} 189 * @param newcentroid A vector with the position of the new centroids 190 */ 191static void shift_codebook(elbg_data *elbg, int *indexes, 192 int *newcentroid[3]) 193{ 194 cell *tempdata; 195 cell **pp = &elbg->cells[indexes[2]]; 196 197 while(*pp) 198 pp= &(*pp)->next; 199 200 *pp = elbg->cells[indexes[0]]; 201 202 elbg->cells[indexes[0]] = NULL; 203 tempdata = elbg->cells[indexes[1]]; 204 elbg->cells[indexes[1]] = NULL; 205 206 while(tempdata) { 207 cell *tempcell2 = tempdata->next; 208 int idx = distance_limited(elbg->points + tempdata->index*elbg->dim, 209 newcentroid[0], elbg->dim, INT_MAX) > 210 distance_limited(elbg->points + tempdata->index*elbg->dim, 211 newcentroid[1], elbg->dim, INT_MAX); 212 213 tempdata->next = elbg->cells[indexes[idx]]; 214 elbg->cells[indexes[idx]] = tempdata; 215 tempdata = tempcell2; 216 } 217} 218 219static void evaluate_utility_inc(elbg_data *elbg) 220{ 221 int i, inc=0; 222 223 for (i=0; i < elbg->numCB; i++) { 224 if (elbg->numCB*elbg->utility[i] > elbg->error) 225 inc += elbg->utility[i]; 226 elbg->utility_inc[i] = inc; 227 } 228} 229 230 231static void update_utility_and_n_cb(elbg_data *elbg, int idx, int newutility) 232{ 233 cell *tempcell; 234 235 elbg->utility[idx] = newutility; 236 for (tempcell=elbg->cells[idx]; tempcell; tempcell=tempcell->next) 237 elbg->nearest_cb[tempcell->index] = idx; 238} 239 240/** 241 * Evaluate if a shift lower the error. If it does, call shift_codebooks 242 * and update elbg->error, elbg->utility and elbg->nearest_cb. 243 * 244 * @param elbg Internal elbg data 245 * @param indexes {luc (low utility cell, huc (high utility cell), cluc (closest cell to low utility cell)} 246 */ 247static void try_shift_candidate(elbg_data *elbg, int idx[3]) 248{ 249 int j, k, olderror=0, newerror, cont=0; 250 int newutility[3]; 251 int newcentroid[3][elbg->dim]; 252 int *newcentroid_ptrs[3]; 253 cell *tempcell; 254 255 newcentroid_ptrs[0] = newcentroid[0]; 256 newcentroid_ptrs[1] = newcentroid[1]; 257 newcentroid_ptrs[2] = newcentroid[2]; 258 259 for (j=0; j<3; j++) 260 olderror += elbg->utility[idx[j]]; 261 262 memset(newcentroid[2], 0, elbg->dim*sizeof(int)); 263 264 for (k=0; k<2; k++) 265 for (tempcell=elbg->cells[idx[2*k]]; tempcell; tempcell=tempcell->next) { 266 cont++; 267 for (j=0; j<elbg->dim; j++) 268 newcentroid[2][j] += elbg->points[tempcell->index*elbg->dim + j]; 269 } 270 271 vect_division(newcentroid[2], newcentroid[2], cont, elbg->dim); 272 273 get_new_centroids(elbg, idx[1], newcentroid[0], newcentroid[1]); 274 275 newutility[2] = eval_error_cell(elbg, newcentroid[2], elbg->cells[idx[0]]); 276 newutility[2] += eval_error_cell(elbg, newcentroid[2], elbg->cells[idx[2]]); 277 278 newerror = newutility[2]; 279 280 newerror += simple_lbg(elbg->dim, newcentroid_ptrs, newutility, elbg->points, 281 elbg->cells[idx[1]]); 282 283 if (olderror > newerror) { 284 shift_codebook(elbg, idx, newcentroid_ptrs); 285 286 elbg->error += newerror - olderror; 287 288 for (j=0; j<3; j++) 289 update_utility_and_n_cb(elbg, idx[j], newutility[j]); 290 291 evaluate_utility_inc(elbg); 292 } 293 } 294 295/** 296 * Implementation of the ELBG block 297 */ 298static void do_shiftings(elbg_data *elbg) 299{ 300 int idx[3]; 301 302 evaluate_utility_inc(elbg); 303 304 for (idx[0]=0; idx[0] < elbg->numCB; idx[0]++) 305 if (elbg->numCB*elbg->utility[idx[0]] < elbg->error) { 306 if (elbg->utility_inc[elbg->numCB-1] == 0) 307 return; 308 309 idx[1] = get_high_utility_cell(elbg); 310 idx[2] = get_closest_codebook(elbg, idx[0]); 311 312 if (idx[1] != idx[0] && idx[1] != idx[2]) 313 try_shift_candidate(elbg, idx); 314 } 315} 316 317#define BIG_PRIME 433494437LL 318 319void ff_init_elbg(int *points, int dim, int numpoints, int *codebook, 320 int numCB, int max_steps, int *closest_cb, 321 AVLFG *rand_state) 322{ 323 int i, k; 324 325 if (numpoints > 24*numCB) { 326 /* ELBG is very costly for a big number of points. So if we have a lot 327 of them, get a good initial codebook to save on iterations */ 328 int *temp_points = av_malloc(dim*(numpoints/8)*sizeof(int)); 329 for (i=0; i<numpoints/8; i++) { 330 k = (i*BIG_PRIME) % numpoints; 331 memcpy(temp_points + i*dim, points + k*dim, dim*sizeof(int)); 332 } 333 334 ff_init_elbg(temp_points, dim, numpoints/8, codebook, numCB, 2*max_steps, closest_cb, rand_state); 335 ff_do_elbg(temp_points, dim, numpoints/8, codebook, numCB, 2*max_steps, closest_cb, rand_state); 336 337 av_free(temp_points); 338 339 } else // If not, initialize the codebook with random positions 340 for (i=0; i < numCB; i++) 341 memcpy(codebook + i*dim, points + ((i*BIG_PRIME)%numpoints)*dim, 342 dim*sizeof(int)); 343 344} 345 346void ff_do_elbg(int *points, int dim, int numpoints, int *codebook, 347 int numCB, int max_steps, int *closest_cb, 348 AVLFG *rand_state) 349{ 350 int dist; 351 elbg_data elbg_d; 352 elbg_data *elbg = &elbg_d; 353 int i, j, k, last_error, steps=0; 354 int *dist_cb = av_malloc(numpoints*sizeof(int)); 355 int *size_part = av_malloc(numCB*sizeof(int)); 356 cell *list_buffer = av_malloc(numpoints*sizeof(cell)); 357 cell *free_cells; 358 int best_dist, best_idx = 0; 359 360 elbg->error = INT_MAX; 361 elbg->dim = dim; 362 elbg->numCB = numCB; 363 elbg->codebook = codebook; 364 elbg->cells = av_malloc(numCB*sizeof(cell *)); 365 elbg->utility = av_malloc(numCB*sizeof(int)); 366 elbg->nearest_cb = closest_cb; 367 elbg->points = points; 368 elbg->utility_inc = av_malloc(numCB*sizeof(int)); 369 370 elbg->rand_state = rand_state; 371 372 do { 373 free_cells = list_buffer; 374 last_error = elbg->error; 375 steps++; 376 memset(elbg->utility, 0, numCB*sizeof(int)); 377 memset(elbg->cells, 0, numCB*sizeof(cell *)); 378 379 elbg->error = 0; 380 381 /* This loop evaluate the actual Voronoi partition. It is the most 382 costly part of the algorithm. */ 383 for (i=0; i < numpoints; i++) { 384 best_dist = distance_limited(elbg->points + i*elbg->dim, elbg->codebook + best_idx*elbg->dim, dim, INT_MAX); 385 for (k=0; k < elbg->numCB; k++) { 386 dist = distance_limited(elbg->points + i*elbg->dim, elbg->codebook + k*elbg->dim, dim, best_dist); 387 if (dist < best_dist) { 388 best_dist = dist; 389 best_idx = k; 390 } 391 } 392 elbg->nearest_cb[i] = best_idx; 393 dist_cb[i] = best_dist; 394 elbg->error += dist_cb[i]; 395 elbg->utility[elbg->nearest_cb[i]] += dist_cb[i]; 396 free_cells->index = i; 397 free_cells->next = elbg->cells[elbg->nearest_cb[i]]; 398 elbg->cells[elbg->nearest_cb[i]] = free_cells; 399 free_cells++; 400 } 401 402 do_shiftings(elbg); 403 404 memset(size_part, 0, numCB*sizeof(int)); 405 406 memset(elbg->codebook, 0, elbg->numCB*dim*sizeof(int)); 407 408 for (i=0; i < numpoints; i++) { 409 size_part[elbg->nearest_cb[i]]++; 410 for (j=0; j < elbg->dim; j++) 411 elbg->codebook[elbg->nearest_cb[i]*elbg->dim + j] += 412 elbg->points[i*elbg->dim + j]; 413 } 414 415 for (i=0; i < elbg->numCB; i++) 416 vect_division(elbg->codebook + i*elbg->dim, 417 elbg->codebook + i*elbg->dim, size_part[i], elbg->dim); 418 419 } while(((last_error - elbg->error) > DELTA_ERR_MAX*elbg->error) && 420 (steps < max_steps)); 421 422 av_free(dist_cb); 423 av_free(size_part); 424 av_free(elbg->utility); 425 av_free(list_buffer); 426 av_free(elbg->cells); 427 av_free(elbg->utility_inc); 428} 429