Libav 0.7.1
|
00001 /* 00002 * Copyright (C) 2007 Vitor Sessak <vitor1001@gmail.com> 00003 * 00004 * This file is part of Libav. 00005 * 00006 * Libav 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 * Libav 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 Libav; 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 int *scratchbuf; 00057 } elbg_data; 00058 00059 static inline int distance_limited(int *a, int *b, int dim, int limit) 00060 { 00061 int i, dist=0; 00062 for (i=0; i<dim; i++) { 00063 dist += (a[i] - b[i])*(a[i] - b[i]); 00064 if (dist > limit) 00065 return INT_MAX; 00066 } 00067 00068 return dist; 00069 } 00070 00071 static inline void vect_division(int *res, int *vect, int div, int dim) 00072 { 00073 int i; 00074 if (div > 1) 00075 for (i=0; i<dim; i++) 00076 res[i] = ROUNDED_DIV(vect[i],div); 00077 else if (res != vect) 00078 memcpy(res, vect, dim*sizeof(int)); 00079 00080 } 00081 00082 static int eval_error_cell(elbg_data *elbg, int *centroid, cell *cells) 00083 { 00084 int error=0; 00085 for (; cells; cells=cells->next) 00086 error += distance_limited(centroid, elbg->points + cells->index*elbg->dim, elbg->dim, INT_MAX); 00087 00088 return error; 00089 } 00090 00091 static int get_closest_codebook(elbg_data *elbg, int index) 00092 { 00093 int i, pick=0, diff, diff_min = INT_MAX; 00094 for (i=0; i<elbg->numCB; i++) 00095 if (i != index) { 00096 diff = distance_limited(elbg->codebook + i*elbg->dim, elbg->codebook + index*elbg->dim, elbg->dim, diff_min); 00097 if (diff < diff_min) { 00098 pick = i; 00099 diff_min = diff; 00100 } 00101 } 00102 return pick; 00103 } 00104 00105 static int get_high_utility_cell(elbg_data *elbg) 00106 { 00107 int i=0; 00108 /* Using linear search, do binary if it ever turns to be speed critical */ 00109 int r = av_lfg_get(elbg->rand_state)%elbg->utility_inc[elbg->numCB-1] + 1; 00110 while (elbg->utility_inc[i] < r) 00111 i++; 00112 00113 assert(!elbg->cells[i]); 00114 00115 return i; 00116 } 00117 00121 static int simple_lbg(elbg_data *elbg, 00122 int dim, 00123 int *centroid[3], 00124 int newutility[3], 00125 int *points, 00126 cell *cells) 00127 { 00128 int i, idx; 00129 int numpoints[2] = {0,0}; 00130 int *newcentroid[2] = { 00131 elbg->scratchbuf + 3*dim, 00132 elbg->scratchbuf + 4*dim 00133 }; 00134 cell *tempcell; 00135 00136 memset(newcentroid[0], 0, 2 * dim * sizeof(*newcentroid[0])); 00137 00138 newutility[0] = 00139 newutility[1] = 0; 00140 00141 for (tempcell = cells; tempcell; tempcell=tempcell->next) { 00142 idx = distance_limited(centroid[0], points + tempcell->index*dim, dim, INT_MAX)>= 00143 distance_limited(centroid[1], points + tempcell->index*dim, dim, INT_MAX); 00144 numpoints[idx]++; 00145 for (i=0; i<dim; i++) 00146 newcentroid[idx][i] += points[tempcell->index*dim + i]; 00147 } 00148 00149 vect_division(centroid[0], newcentroid[0], numpoints[0], dim); 00150 vect_division(centroid[1], newcentroid[1], numpoints[1], dim); 00151 00152 for (tempcell = cells; tempcell; tempcell=tempcell->next) { 00153 int dist[2] = {distance_limited(centroid[0], points + tempcell->index*dim, dim, INT_MAX), 00154 distance_limited(centroid[1], points + tempcell->index*dim, dim, INT_MAX)}; 00155 int idx = dist[0] > dist[1]; 00156 newutility[idx] += dist[idx]; 00157 } 00158 00159 return newutility[0] + newutility[1]; 00160 } 00161 00162 static void get_new_centroids(elbg_data *elbg, int huc, int *newcentroid_i, 00163 int *newcentroid_p) 00164 { 00165 cell *tempcell; 00166 int *min = newcentroid_i; 00167 int *max = newcentroid_p; 00168 int i; 00169 00170 for (i=0; i< elbg->dim; i++) { 00171 min[i]=INT_MAX; 00172 max[i]=0; 00173 } 00174 00175 for (tempcell = elbg->cells[huc]; tempcell; tempcell = tempcell->next) 00176 for(i=0; i<elbg->dim; i++) { 00177 min[i]=FFMIN(min[i], elbg->points[tempcell->index*elbg->dim + i]); 00178 max[i]=FFMAX(max[i], elbg->points[tempcell->index*elbg->dim + i]); 00179 } 00180 00181 for (i=0; i<elbg->dim; i++) { 00182 int ni = min[i] + (max[i] - min[i])/3; 00183 int np = min[i] + (2*(max[i] - min[i]))/3; 00184 newcentroid_i[i] = ni; 00185 newcentroid_p[i] = np; 00186 } 00187 } 00188 00198 static void shift_codebook(elbg_data *elbg, int *indexes, 00199 int *newcentroid[3]) 00200 { 00201 cell *tempdata; 00202 cell **pp = &elbg->cells[indexes[2]]; 00203 00204 while(*pp) 00205 pp= &(*pp)->next; 00206 00207 *pp = elbg->cells[indexes[0]]; 00208 00209 elbg->cells[indexes[0]] = NULL; 00210 tempdata = elbg->cells[indexes[1]]; 00211 elbg->cells[indexes[1]] = NULL; 00212 00213 while(tempdata) { 00214 cell *tempcell2 = tempdata->next; 00215 int idx = distance_limited(elbg->points + tempdata->index*elbg->dim, 00216 newcentroid[0], elbg->dim, INT_MAX) > 00217 distance_limited(elbg->points + tempdata->index*elbg->dim, 00218 newcentroid[1], elbg->dim, INT_MAX); 00219 00220 tempdata->next = elbg->cells[indexes[idx]]; 00221 elbg->cells[indexes[idx]] = tempdata; 00222 tempdata = tempcell2; 00223 } 00224 } 00225 00226 static void evaluate_utility_inc(elbg_data *elbg) 00227 { 00228 int i, inc=0; 00229 00230 for (i=0; i < elbg->numCB; i++) { 00231 if (elbg->numCB*elbg->utility[i] > elbg->error) 00232 inc += elbg->utility[i]; 00233 elbg->utility_inc[i] = inc; 00234 } 00235 } 00236 00237 00238 static void update_utility_and_n_cb(elbg_data *elbg, int idx, int newutility) 00239 { 00240 cell *tempcell; 00241 00242 elbg->utility[idx] = newutility; 00243 for (tempcell=elbg->cells[idx]; tempcell; tempcell=tempcell->next) 00244 elbg->nearest_cb[tempcell->index] = idx; 00245 } 00246 00254 static void try_shift_candidate(elbg_data *elbg, int idx[3]) 00255 { 00256 int j, k, olderror=0, newerror, cont=0; 00257 int newutility[3]; 00258 int *newcentroid[3] = { 00259 elbg->scratchbuf, 00260 elbg->scratchbuf + elbg->dim, 00261 elbg->scratchbuf + 2*elbg->dim 00262 }; 00263 cell *tempcell; 00264 00265 for (j=0; j<3; j++) 00266 olderror += elbg->utility[idx[j]]; 00267 00268 memset(newcentroid[2], 0, elbg->dim*sizeof(int)); 00269 00270 for (k=0; k<2; k++) 00271 for (tempcell=elbg->cells[idx[2*k]]; tempcell; tempcell=tempcell->next) { 00272 cont++; 00273 for (j=0; j<elbg->dim; j++) 00274 newcentroid[2][j] += elbg->points[tempcell->index*elbg->dim + j]; 00275 } 00276 00277 vect_division(newcentroid[2], newcentroid[2], cont, elbg->dim); 00278 00279 get_new_centroids(elbg, idx[1], newcentroid[0], newcentroid[1]); 00280 00281 newutility[2] = eval_error_cell(elbg, newcentroid[2], elbg->cells[idx[0]]); 00282 newutility[2] += eval_error_cell(elbg, newcentroid[2], elbg->cells[idx[2]]); 00283 00284 newerror = newutility[2]; 00285 00286 newerror += simple_lbg(elbg, elbg->dim, newcentroid, newutility, elbg->points, 00287 elbg->cells[idx[1]]); 00288 00289 if (olderror > newerror) { 00290 shift_codebook(elbg, idx, newcentroid); 00291 00292 elbg->error += newerror - olderror; 00293 00294 for (j=0; j<3; j++) 00295 update_utility_and_n_cb(elbg, idx[j], newutility[j]); 00296 00297 evaluate_utility_inc(elbg); 00298 } 00299 } 00300 00304 static void do_shiftings(elbg_data *elbg) 00305 { 00306 int idx[3]; 00307 00308 evaluate_utility_inc(elbg); 00309 00310 for (idx[0]=0; idx[0] < elbg->numCB; idx[0]++) 00311 if (elbg->numCB*elbg->utility[idx[0]] < elbg->error) { 00312 if (elbg->utility_inc[elbg->numCB-1] == 0) 00313 return; 00314 00315 idx[1] = get_high_utility_cell(elbg); 00316 idx[2] = get_closest_codebook(elbg, idx[0]); 00317 00318 if (idx[1] != idx[0] && idx[1] != idx[2]) 00319 try_shift_candidate(elbg, idx); 00320 } 00321 } 00322 00323 #define BIG_PRIME 433494437LL 00324 00325 void ff_init_elbg(int *points, int dim, int numpoints, int *codebook, 00326 int numCB, int max_steps, int *closest_cb, 00327 AVLFG *rand_state) 00328 { 00329 int i, k; 00330 00331 if (numpoints > 24*numCB) { 00332 /* ELBG is very costly for a big number of points. So if we have a lot 00333 of them, get a good initial codebook to save on iterations */ 00334 int *temp_points = av_malloc(dim*(numpoints/8)*sizeof(int)); 00335 for (i=0; i<numpoints/8; i++) { 00336 k = (i*BIG_PRIME) % numpoints; 00337 memcpy(temp_points + i*dim, points + k*dim, dim*sizeof(int)); 00338 } 00339 00340 ff_init_elbg(temp_points, dim, numpoints/8, codebook, numCB, 2*max_steps, closest_cb, rand_state); 00341 ff_do_elbg(temp_points, dim, numpoints/8, codebook, numCB, 2*max_steps, closest_cb, rand_state); 00342 00343 av_free(temp_points); 00344 00345 } else // If not, initialize the codebook with random positions 00346 for (i=0; i < numCB; i++) 00347 memcpy(codebook + i*dim, points + ((i*BIG_PRIME)%numpoints)*dim, 00348 dim*sizeof(int)); 00349 00350 } 00351 00352 void ff_do_elbg(int *points, int dim, int numpoints, int *codebook, 00353 int numCB, int max_steps, int *closest_cb, 00354 AVLFG *rand_state) 00355 { 00356 int dist; 00357 elbg_data elbg_d; 00358 elbg_data *elbg = &elbg_d; 00359 int i, j, k, last_error, steps=0; 00360 int *dist_cb = av_malloc(numpoints*sizeof(int)); 00361 int *size_part = av_malloc(numCB*sizeof(int)); 00362 cell *list_buffer = av_malloc(numpoints*sizeof(cell)); 00363 cell *free_cells; 00364 int best_dist, best_idx = 0; 00365 00366 elbg->error = INT_MAX; 00367 elbg->dim = dim; 00368 elbg->numCB = numCB; 00369 elbg->codebook = codebook; 00370 elbg->cells = av_malloc(numCB*sizeof(cell *)); 00371 elbg->utility = av_malloc(numCB*sizeof(int)); 00372 elbg->nearest_cb = closest_cb; 00373 elbg->points = points; 00374 elbg->utility_inc = av_malloc(numCB*sizeof(int)); 00375 elbg->scratchbuf = av_malloc(5*dim*sizeof(int)); 00376 00377 elbg->rand_state = rand_state; 00378 00379 do { 00380 free_cells = list_buffer; 00381 last_error = elbg->error; 00382 steps++; 00383 memset(elbg->utility, 0, numCB*sizeof(int)); 00384 memset(elbg->cells, 0, numCB*sizeof(cell *)); 00385 00386 elbg->error = 0; 00387 00388 /* This loop evaluate the actual Voronoi partition. It is the most 00389 costly part of the algorithm. */ 00390 for (i=0; i < numpoints; i++) { 00391 best_dist = distance_limited(elbg->points + i*elbg->dim, elbg->codebook + best_idx*elbg->dim, dim, INT_MAX); 00392 for (k=0; k < elbg->numCB; k++) { 00393 dist = distance_limited(elbg->points + i*elbg->dim, elbg->codebook + k*elbg->dim, dim, best_dist); 00394 if (dist < best_dist) { 00395 best_dist = dist; 00396 best_idx = k; 00397 } 00398 } 00399 elbg->nearest_cb[i] = best_idx; 00400 dist_cb[i] = best_dist; 00401 elbg->error += dist_cb[i]; 00402 elbg->utility[elbg->nearest_cb[i]] += dist_cb[i]; 00403 free_cells->index = i; 00404 free_cells->next = elbg->cells[elbg->nearest_cb[i]]; 00405 elbg->cells[elbg->nearest_cb[i]] = free_cells; 00406 free_cells++; 00407 } 00408 00409 do_shiftings(elbg); 00410 00411 memset(size_part, 0, numCB*sizeof(int)); 00412 00413 memset(elbg->codebook, 0, elbg->numCB*dim*sizeof(int)); 00414 00415 for (i=0; i < numpoints; i++) { 00416 size_part[elbg->nearest_cb[i]]++; 00417 for (j=0; j < elbg->dim; j++) 00418 elbg->codebook[elbg->nearest_cb[i]*elbg->dim + j] += 00419 elbg->points[i*elbg->dim + j]; 00420 } 00421 00422 for (i=0; i < elbg->numCB; i++) 00423 vect_division(elbg->codebook + i*elbg->dim, 00424 elbg->codebook + i*elbg->dim, size_part[i], elbg->dim); 00425 00426 } while(((last_error - elbg->error) > DELTA_ERR_MAX*elbg->error) && 00427 (steps < max_steps)); 00428 00429 av_free(dist_cb); 00430 av_free(size_part); 00431 av_free(elbg->utility); 00432 av_free(list_buffer); 00433 av_free(elbg->cells); 00434 av_free(elbg->utility_inc); 00435 av_free(elbg->scratchbuf); 00436 }