• Main Page
  • Related Pages
  • Modules
  • Data Structures
  • Files
  • File List
  • Globals

libavcodec/elbg.c

Go to the documentation of this file.
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 }

Generated on Fri Sep 16 2011 17:17:36 for FFmpeg by  doxygen 1.7.1