Libav
|
00001 /* 00002 * Copyright (C) 2007 Vitor Sessak <vitor1001@gmail.com> 00003 * 00004 * This file is part of FFmpeg. 00005 * 00006 * FFmpeg is free software; you can redistribute it and/or 00007 * modify it under the terms of the GNU Lesser General Public 00008 * License as published by the Free Software Foundation; either 00009 * version 2.1 of the License, or (at your option) any later version. 00010 * 00011 * FFmpeg is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00014 * Lesser General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU Lesser General Public 00017 * License along with FFmpeg; if not, write to the Free Software 00018 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00019 */ 00020 00026 #include <string.h> 00027 00028 #include "libavutil/lfg.h" 00029 #include "elbg.h" 00030 #include "avcodec.h" 00031 00032 #define DELTA_ERR_MAX 0.1 ///< Precision of the ELBG algorithm (as percentual error) 00033 00037 typedef struct cell_s { 00038 int index; 00039 struct cell_s *next; 00040 } cell; 00041 00045 typedef struct{ 00046 int error; 00047 int dim; 00048 int numCB; 00049 int *codebook; 00050 cell **cells; 00051 int *utility; 00052 int *utility_inc; 00053 int *nearest_cb; 00054 int *points; 00055 AVLFG *rand_state; 00056 } elbg_data; 00057 00058 static inline int distance_limited(int *a, int *b, int dim, int limit) 00059 { 00060 int i, dist=0; 00061 for (i=0; i<dim; i++) { 00062 dist += (a[i] - b[i])*(a[i] - b[i]); 00063 if (dist > limit) 00064 return INT_MAX; 00065 } 00066 00067 return dist; 00068 } 00069 00070 static inline void vect_division(int *res, int *vect, int div, int dim) 00071 { 00072 int i; 00073 if (div > 1) 00074 for (i=0; i<dim; i++) 00075 res[i] = ROUNDED_DIV(vect[i],div); 00076 else if (res != vect) 00077 memcpy(res, vect, dim*sizeof(int)); 00078 00079 } 00080 00081 static int eval_error_cell(elbg_data *elbg, int *centroid, cell *cells) 00082 { 00083 int error=0; 00084 for (; cells; cells=cells->next) 00085 error += distance_limited(centroid, elbg->points + cells->index*elbg->dim, elbg->dim, INT_MAX); 00086 00087 return error; 00088 } 00089 00090 static int get_closest_codebook(elbg_data *elbg, int index) 00091 { 00092 int i, pick=0, diff, diff_min = INT_MAX; 00093 for (i=0; i<elbg->numCB; i++) 00094 if (i != index) { 00095 diff = distance_limited(elbg->codebook + i*elbg->dim, elbg->codebook + index*elbg->dim, elbg->dim, diff_min); 00096 if (diff < diff_min) { 00097 pick = i; 00098 diff_min = diff; 00099 } 00100 } 00101 return pick; 00102 } 00103 00104 static int get_high_utility_cell(elbg_data *elbg) 00105 { 00106 int i=0; 00107 /* Using linear search, do binary if it ever turns to be speed critical */ 00108 int r = av_lfg_get(elbg->rand_state)%elbg->utility_inc[elbg->numCB-1] + 1; 00109 while (elbg->utility_inc[i] < r) 00110 i++; 00111 00112 assert(!elbg->cells[i]); 00113 00114 return i; 00115 } 00116 00120 static int simple_lbg(int dim, 00121 int *centroid[3], 00122 int newutility[3], 00123 int *points, 00124 cell *cells) 00125 { 00126 int i, idx; 00127 int numpoints[2] = {0,0}; 00128 int newcentroid[2][dim]; 00129 cell *tempcell; 00130 00131 memset(newcentroid, 0, sizeof(newcentroid)); 00132 00133 newutility[0] = 00134 newutility[1] = 0; 00135 00136 for (tempcell = cells; tempcell; tempcell=tempcell->next) { 00137 idx = distance_limited(centroid[0], points + tempcell->index*dim, dim, INT_MAX)>= 00138 distance_limited(centroid[1], points + tempcell->index*dim, dim, INT_MAX); 00139 numpoints[idx]++; 00140 for (i=0; i<dim; i++) 00141 newcentroid[idx][i] += points[tempcell->index*dim + i]; 00142 } 00143 00144 vect_division(centroid[0], newcentroid[0], numpoints[0], dim); 00145 vect_division(centroid[1], newcentroid[1], numpoints[1], dim); 00146 00147 for (tempcell = cells; tempcell; tempcell=tempcell->next) { 00148 int dist[2] = {distance_limited(centroid[0], points + tempcell->index*dim, dim, INT_MAX), 00149 distance_limited(centroid[1], points + tempcell->index*dim, dim, INT_MAX)}; 00150 int idx = dist[0] > dist[1]; 00151 newutility[idx] += dist[idx]; 00152 } 00153 00154 return newutility[0] + newutility[1]; 00155 } 00156 00157 static void get_new_centroids(elbg_data *elbg, int huc, int *newcentroid_i, 00158 int *newcentroid_p) 00159 { 00160 cell *tempcell; 00161 int min[elbg->dim]; 00162 int max[elbg->dim]; 00163 int i; 00164 00165 for (i=0; i< elbg->dim; i++) { 00166 min[i]=INT_MAX; 00167 max[i]=0; 00168 } 00169 00170 for (tempcell = elbg->cells[huc]; tempcell; tempcell = tempcell->next) 00171 for(i=0; i<elbg->dim; i++) { 00172 min[i]=FFMIN(min[i], elbg->points[tempcell->index*elbg->dim + i]); 00173 max[i]=FFMAX(max[i], elbg->points[tempcell->index*elbg->dim + i]); 00174 } 00175 00176 for (i=0; i<elbg->dim; i++) { 00177 newcentroid_i[i] = min[i] + (max[i] - min[i])/3; 00178 newcentroid_p[i] = min[i] + (2*(max[i] - min[i]))/3; 00179 } 00180 } 00181 00191 static void shift_codebook(elbg_data *elbg, int *indexes, 00192 int *newcentroid[3]) 00193 { 00194 cell *tempdata; 00195 cell **pp = &elbg->cells[indexes[2]]; 00196 00197 while(*pp) 00198 pp= &(*pp)->next; 00199 00200 *pp = elbg->cells[indexes[0]]; 00201 00202 elbg->cells[indexes[0]] = NULL; 00203 tempdata = elbg->cells[indexes[1]]; 00204 elbg->cells[indexes[1]] = NULL; 00205 00206 while(tempdata) { 00207 cell *tempcell2 = tempdata->next; 00208 int idx = distance_limited(elbg->points + tempdata->index*elbg->dim, 00209 newcentroid[0], elbg->dim, INT_MAX) > 00210 distance_limited(elbg->points + tempdata->index*elbg->dim, 00211 newcentroid[1], elbg->dim, INT_MAX); 00212 00213 tempdata->next = elbg->cells[indexes[idx]]; 00214 elbg->cells[indexes[idx]] = tempdata; 00215 tempdata = tempcell2; 00216 } 00217 } 00218 00219 static void evaluate_utility_inc(elbg_data *elbg) 00220 { 00221 int i, inc=0; 00222 00223 for (i=0; i < elbg->numCB; i++) { 00224 if (elbg->numCB*elbg->utility[i] > elbg->error) 00225 inc += elbg->utility[i]; 00226 elbg->utility_inc[i] = inc; 00227 } 00228 } 00229 00230 00231 static void update_utility_and_n_cb(elbg_data *elbg, int idx, int newutility) 00232 { 00233 cell *tempcell; 00234 00235 elbg->utility[idx] = newutility; 00236 for (tempcell=elbg->cells[idx]; tempcell; tempcell=tempcell->next) 00237 elbg->nearest_cb[tempcell->index] = idx; 00238 } 00239 00247 static void try_shift_candidate(elbg_data *elbg, int idx[3]) 00248 { 00249 int j, k, olderror=0, newerror, cont=0; 00250 int newutility[3]; 00251 int newcentroid[3][elbg->dim]; 00252 int *newcentroid_ptrs[3]; 00253 cell *tempcell; 00254 00255 newcentroid_ptrs[0] = newcentroid[0]; 00256 newcentroid_ptrs[1] = newcentroid[1]; 00257 newcentroid_ptrs[2] = newcentroid[2]; 00258 00259 for (j=0; j<3; j++) 00260 olderror += elbg->utility[idx[j]]; 00261 00262 memset(newcentroid[2], 0, elbg->dim*sizeof(int)); 00263 00264 for (k=0; k<2; k++) 00265 for (tempcell=elbg->cells[idx[2*k]]; tempcell; tempcell=tempcell->next) { 00266 cont++; 00267 for (j=0; j<elbg->dim; j++) 00268 newcentroid[2][j] += elbg->points[tempcell->index*elbg->dim + j]; 00269 } 00270 00271 vect_division(newcentroid[2], newcentroid[2], cont, elbg->dim); 00272 00273 get_new_centroids(elbg, idx[1], newcentroid[0], newcentroid[1]); 00274 00275 newutility[2] = eval_error_cell(elbg, newcentroid[2], elbg->cells[idx[0]]); 00276 newutility[2] += eval_error_cell(elbg, newcentroid[2], elbg->cells[idx[2]]); 00277 00278 newerror = newutility[2]; 00279 00280 newerror += simple_lbg(elbg->dim, newcentroid_ptrs, newutility, elbg->points, 00281 elbg->cells[idx[1]]); 00282 00283 if (olderror > newerror) { 00284 shift_codebook(elbg, idx, newcentroid_ptrs); 00285 00286 elbg->error += newerror - olderror; 00287 00288 for (j=0; j<3; j++) 00289 update_utility_and_n_cb(elbg, idx[j], newutility[j]); 00290 00291 evaluate_utility_inc(elbg); 00292 } 00293 } 00294 00298 static void do_shiftings(elbg_data *elbg) 00299 { 00300 int idx[3]; 00301 00302 evaluate_utility_inc(elbg); 00303 00304 for (idx[0]=0; idx[0] < elbg->numCB; idx[0]++) 00305 if (elbg->numCB*elbg->utility[idx[0]] < elbg->error) { 00306 if (elbg->utility_inc[elbg->numCB-1] == 0) 00307 return; 00308 00309 idx[1] = get_high_utility_cell(elbg); 00310 idx[2] = get_closest_codebook(elbg, idx[0]); 00311 00312 if (idx[1] != idx[0] && idx[1] != idx[2]) 00313 try_shift_candidate(elbg, idx); 00314 } 00315 } 00316 00317 #define BIG_PRIME 433494437LL 00318 00319 void ff_init_elbg(int *points, int dim, int numpoints, int *codebook, 00320 int numCB, int max_steps, int *closest_cb, 00321 AVLFG *rand_state) 00322 { 00323 int i, k; 00324 00325 if (numpoints > 24*numCB) { 00326 /* ELBG is very costly for a big number of points. So if we have a lot 00327 of them, get a good initial codebook to save on iterations */ 00328 int *temp_points = av_malloc(dim*(numpoints/8)*sizeof(int)); 00329 for (i=0; i<numpoints/8; i++) { 00330 k = (i*BIG_PRIME) % numpoints; 00331 memcpy(temp_points + i*dim, points + k*dim, dim*sizeof(int)); 00332 } 00333 00334 ff_init_elbg(temp_points, dim, numpoints/8, codebook, numCB, 2*max_steps, closest_cb, rand_state); 00335 ff_do_elbg(temp_points, dim, numpoints/8, codebook, numCB, 2*max_steps, closest_cb, rand_state); 00336 00337 av_free(temp_points); 00338 00339 } else // If not, initialize the codebook with random positions 00340 for (i=0; i < numCB; i++) 00341 memcpy(codebook + i*dim, points + ((i*BIG_PRIME)%numpoints)*dim, 00342 dim*sizeof(int)); 00343 00344 } 00345 00346 void ff_do_elbg(int *points, int dim, int numpoints, int *codebook, 00347 int numCB, int max_steps, int *closest_cb, 00348 AVLFG *rand_state) 00349 { 00350 int dist; 00351 elbg_data elbg_d; 00352 elbg_data *elbg = &elbg_d; 00353 int i, j, k, last_error, steps=0; 00354 int *dist_cb = av_malloc(numpoints*sizeof(int)); 00355 int *size_part = av_malloc(numCB*sizeof(int)); 00356 cell *list_buffer = av_malloc(numpoints*sizeof(cell)); 00357 cell *free_cells; 00358 int best_dist, best_idx = 0; 00359 00360 elbg->error = INT_MAX; 00361 elbg->dim = dim; 00362 elbg->numCB = numCB; 00363 elbg->codebook = codebook; 00364 elbg->cells = av_malloc(numCB*sizeof(cell *)); 00365 elbg->utility = av_malloc(numCB*sizeof(int)); 00366 elbg->nearest_cb = closest_cb; 00367 elbg->points = points; 00368 elbg->utility_inc = av_malloc(numCB*sizeof(int)); 00369 00370 elbg->rand_state = rand_state; 00371 00372 do { 00373 free_cells = list_buffer; 00374 last_error = elbg->error; 00375 steps++; 00376 memset(elbg->utility, 0, numCB*sizeof(int)); 00377 memset(elbg->cells, 0, numCB*sizeof(cell *)); 00378 00379 elbg->error = 0; 00380 00381 /* This loop evaluate the actual Voronoi partition. It is the most 00382 costly part of the algorithm. */ 00383 for (i=0; i < numpoints; i++) { 00384 best_dist = distance_limited(elbg->points + i*elbg->dim, elbg->codebook + best_idx*elbg->dim, dim, INT_MAX); 00385 for (k=0; k < elbg->numCB; k++) { 00386 dist = distance_limited(elbg->points + i*elbg->dim, elbg->codebook + k*elbg->dim, dim, best_dist); 00387 if (dist < best_dist) { 00388 best_dist = dist; 00389 best_idx = k; 00390 } 00391 } 00392 elbg->nearest_cb[i] = best_idx; 00393 dist_cb[i] = best_dist; 00394 elbg->error += dist_cb[i]; 00395 elbg->utility[elbg->nearest_cb[i]] += dist_cb[i]; 00396 free_cells->index = i; 00397 free_cells->next = elbg->cells[elbg->nearest_cb[i]]; 00398 elbg->cells[elbg->nearest_cb[i]] = free_cells; 00399 free_cells++; 00400 } 00401 00402 do_shiftings(elbg); 00403 00404 memset(size_part, 0, numCB*sizeof(int)); 00405 00406 memset(elbg->codebook, 0, elbg->numCB*dim*sizeof(int)); 00407 00408 for (i=0; i < numpoints; i++) { 00409 size_part[elbg->nearest_cb[i]]++; 00410 for (j=0; j < elbg->dim; j++) 00411 elbg->codebook[elbg->nearest_cb[i]*elbg->dim + j] += 00412 elbg->points[i*elbg->dim + j]; 00413 } 00414 00415 for (i=0; i < elbg->numCB; i++) 00416 vect_division(elbg->codebook + i*elbg->dim, 00417 elbg->codebook + i*elbg->dim, size_part[i], elbg->dim); 00418 00419 } while(((last_error - elbg->error) > DELTA_ERR_MAX*elbg->error) && 00420 (steps < max_steps)); 00421 00422 av_free(dist_cb); 00423 av_free(size_part); 00424 av_free(elbg->utility); 00425 av_free(list_buffer); 00426 av_free(elbg->cells); 00427 av_free(elbg->utility_inc); 00428 }