FFmpeg
elbg.c
Go to the documentation of this file.
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/avassert.h"
29 #include "libavutil/common.h"
30 #include "libavutil/lfg.h"
31 #include "elbg.h"
32 #include "avcodec.h"
33 
34 #define DELTA_ERR_MAX 0.1 ///< Precision of the ELBG algorithm (as percentage error)
35 
36 /**
37  * In the ELBG jargon, a cell is the set of points that are closest to a
38  * codebook entry. Not to be confused with a RoQ Video cell. */
39 typedef struct cell_s {
40  int index;
41  struct cell_s *next;
42 } cell;
43 
44 /**
45  * ELBG internal data
46  */
47 typedef struct elbg_data {
48  int error;
49  int dim;
50  int numCB;
51  int *codebook;
52  cell **cells;
53  int *utility;
54  int64_t *utility_inc;
55  int *nearest_cb;
56  int *points;
58  int *scratchbuf;
59 } elbg_data;
60 
61 static inline int distance_limited(int *a, int *b, int dim, int limit)
62 {
63  int i, dist=0;
64  for (i=0; i<dim; i++) {
65  dist += (a[i] - b[i])*(a[i] - b[i]);
66  if (dist > limit)
67  return INT_MAX;
68  }
69 
70  return dist;
71 }
72 
73 static inline void vect_division(int *res, int *vect, int div, int dim)
74 {
75  int i;
76  if (div > 1)
77  for (i=0; i<dim; i++)
78  res[i] = ROUNDED_DIV(vect[i],div);
79  else if (res != vect)
80  memcpy(res, vect, dim*sizeof(int));
81 
82 }
83 
84 static int eval_error_cell(elbg_data *elbg, int *centroid, cell *cells)
85 {
86  int error=0;
87  for (; cells; cells=cells->next)
88  error += distance_limited(centroid, elbg->points + cells->index*elbg->dim, elbg->dim, INT_MAX);
89 
90  return error;
91 }
92 
93 static int get_closest_codebook(elbg_data *elbg, int index)
94 {
95  int i, pick=0, diff, diff_min = INT_MAX;
96  for (i=0; i<elbg->numCB; i++)
97  if (i != index) {
98  diff = distance_limited(elbg->codebook + i*elbg->dim, elbg->codebook + index*elbg->dim, elbg->dim, diff_min);
99  if (diff < diff_min) {
100  pick = i;
101  diff_min = diff;
102  }
103  }
104  return pick;
105 }
106 
108 {
109  int i=0;
110  /* Using linear search, do binary if it ever turns to be speed critical */
111  uint64_t r;
112 
113  if (elbg->utility_inc[elbg->numCB-1] < INT_MAX) {
114  r = av_lfg_get(elbg->rand_state) % (unsigned int)elbg->utility_inc[elbg->numCB-1] + 1;
115  } else {
116  r = av_lfg_get(elbg->rand_state);
117  r = (av_lfg_get(elbg->rand_state) + (r<<32)) % elbg->utility_inc[elbg->numCB-1] + 1;
118  }
119 
120  while (elbg->utility_inc[i] < r) {
121  i++;
122  }
123 
124  av_assert2(elbg->cells[i]);
125 
126  return i;
127 }
128 
129 /**
130  * Implementation of the simple LBG algorithm for just two codebooks
131  */
132 static int simple_lbg(elbg_data *elbg,
133  int dim,
134  int *centroid[3],
135  int newutility[3],
136  int *points,
137  cell *cells)
138 {
139  int i, idx;
140  int numpoints[2] = {0,0};
141  int *newcentroid[2] = {
142  elbg->scratchbuf + 3*dim,
143  elbg->scratchbuf + 4*dim
144  };
145  cell *tempcell;
146 
147  memset(newcentroid[0], 0, 2 * dim * sizeof(*newcentroid[0]));
148 
149  newutility[0] =
150  newutility[1] = 0;
151 
152  for (tempcell = cells; tempcell; tempcell=tempcell->next) {
153  idx = distance_limited(centroid[0], points + tempcell->index*dim, dim, INT_MAX)>=
154  distance_limited(centroid[1], points + tempcell->index*dim, dim, INT_MAX);
155  numpoints[idx]++;
156  for (i=0; i<dim; i++)
157  newcentroid[idx][i] += points[tempcell->index*dim + i];
158  }
159 
160  vect_division(centroid[0], newcentroid[0], numpoints[0], dim);
161  vect_division(centroid[1], newcentroid[1], numpoints[1], dim);
162 
163  for (tempcell = cells; tempcell; tempcell=tempcell->next) {
164  int dist[2] = {distance_limited(centroid[0], points + tempcell->index*dim, dim, INT_MAX),
165  distance_limited(centroid[1], points + tempcell->index*dim, dim, INT_MAX)};
166  int idx = dist[0] > dist[1];
167  newutility[idx] += dist[idx];
168  }
169 
170  return newutility[0] + newutility[1];
171 }
172 
173 static void get_new_centroids(elbg_data *elbg, int huc, int *newcentroid_i,
174  int *newcentroid_p)
175 {
176  cell *tempcell;
177  int *min = newcentroid_i;
178  int *max = newcentroid_p;
179  int i;
180 
181  for (i=0; i< elbg->dim; i++) {
182  min[i]=INT_MAX;
183  max[i]=0;
184  }
185 
186  for (tempcell = elbg->cells[huc]; tempcell; tempcell = tempcell->next)
187  for(i=0; i<elbg->dim; i++) {
188  min[i]=FFMIN(min[i], elbg->points[tempcell->index*elbg->dim + i]);
189  max[i]=FFMAX(max[i], elbg->points[tempcell->index*elbg->dim + i]);
190  }
191 
192  for (i=0; i<elbg->dim; i++) {
193  int ni = min[i] + (max[i] - min[i])/3;
194  int np = min[i] + (2*(max[i] - min[i]))/3;
195  newcentroid_i[i] = ni;
196  newcentroid_p[i] = np;
197  }
198 }
199 
200 /**
201  * Add the points in the low utility cell to its closest cell. Split the high
202  * utility cell, putting the separated points in the (now empty) low utility
203  * cell.
204  *
205  * @param elbg Internal elbg data
206  * @param indexes {luc, huc, cluc}
207  * @param newcentroid A vector with the position of the new centroids
208  */
209 static void shift_codebook(elbg_data *elbg, int *indexes,
210  int *newcentroid[3])
211 {
212  cell *tempdata;
213  cell **pp = &elbg->cells[indexes[2]];
214 
215  while(*pp)
216  pp= &(*pp)->next;
217 
218  *pp = elbg->cells[indexes[0]];
219 
220  elbg->cells[indexes[0]] = NULL;
221  tempdata = elbg->cells[indexes[1]];
222  elbg->cells[indexes[1]] = NULL;
223 
224  while(tempdata) {
225  cell *tempcell2 = tempdata->next;
226  int idx = distance_limited(elbg->points + tempdata->index*elbg->dim,
227  newcentroid[0], elbg->dim, INT_MAX) >
228  distance_limited(elbg->points + tempdata->index*elbg->dim,
229  newcentroid[1], elbg->dim, INT_MAX);
230 
231  tempdata->next = elbg->cells[indexes[idx]];
232  elbg->cells[indexes[idx]] = tempdata;
233  tempdata = tempcell2;
234  }
235 }
236 
237 static void evaluate_utility_inc(elbg_data *elbg)
238 {
239  int i;
240  int64_t inc=0;
241 
242  for (i=0; i < elbg->numCB; i++) {
243  if (elbg->numCB*elbg->utility[i] > elbg->error)
244  inc += elbg->utility[i];
245  elbg->utility_inc[i] = inc;
246  }
247 }
248 
249 
250 static void update_utility_and_n_cb(elbg_data *elbg, int idx, int newutility)
251 {
252  cell *tempcell;
253 
254  elbg->utility[idx] = newutility;
255  for (tempcell=elbg->cells[idx]; tempcell; tempcell=tempcell->next)
256  elbg->nearest_cb[tempcell->index] = idx;
257 }
258 
259 /**
260  * Evaluate if a shift lower the error. If it does, call shift_codebooks
261  * and update elbg->error, elbg->utility and elbg->nearest_cb.
262  *
263  * @param elbg Internal elbg data
264  * @param idx {luc (low utility cell, huc (high utility cell), cluc (closest cell to low utility cell)}
265  */
266 static void try_shift_candidate(elbg_data *elbg, int idx[3])
267 {
268  int j, k, olderror=0, newerror, cont=0;
269  int newutility[3];
270  int *newcentroid[3] = {
271  elbg->scratchbuf,
272  elbg->scratchbuf + elbg->dim,
273  elbg->scratchbuf + 2*elbg->dim
274  };
275  cell *tempcell;
276 
277  for (j=0; j<3; j++)
278  olderror += elbg->utility[idx[j]];
279 
280  memset(newcentroid[2], 0, elbg->dim*sizeof(int));
281 
282  for (k=0; k<2; k++)
283  for (tempcell=elbg->cells[idx[2*k]]; tempcell; tempcell=tempcell->next) {
284  cont++;
285  for (j=0; j<elbg->dim; j++)
286  newcentroid[2][j] += elbg->points[tempcell->index*elbg->dim + j];
287  }
288 
289  vect_division(newcentroid[2], newcentroid[2], cont, elbg->dim);
290 
291  get_new_centroids(elbg, idx[1], newcentroid[0], newcentroid[1]);
292 
293  newutility[2] = eval_error_cell(elbg, newcentroid[2], elbg->cells[idx[0]]);
294  newutility[2] += eval_error_cell(elbg, newcentroid[2], elbg->cells[idx[2]]);
295 
296  newerror = newutility[2];
297 
298  newerror += simple_lbg(elbg, elbg->dim, newcentroid, newutility, elbg->points,
299  elbg->cells[idx[1]]);
300 
301  if (olderror > newerror) {
302  shift_codebook(elbg, idx, newcentroid);
303 
304  elbg->error += newerror - olderror;
305 
306  for (j=0; j<3; j++)
307  update_utility_and_n_cb(elbg, idx[j], newutility[j]);
308 
309  evaluate_utility_inc(elbg);
310  }
311  }
312 
313 /**
314  * Implementation of the ELBG block
315  */
316 static void do_shiftings(elbg_data *elbg)
317 {
318  int idx[3];
319 
320  evaluate_utility_inc(elbg);
321 
322  for (idx[0]=0; idx[0] < elbg->numCB; idx[0]++)
323  if (elbg->numCB*elbg->utility[idx[0]] < elbg->error) {
324  if (elbg->utility_inc[elbg->numCB-1] == 0)
325  return;
326 
327  idx[1] = get_high_utility_cell(elbg);
328  idx[2] = get_closest_codebook(elbg, idx[0]);
329 
330  if (idx[1] != idx[0] && idx[1] != idx[2])
331  try_shift_candidate(elbg, idx);
332  }
333 }
334 
335 #define BIG_PRIME 433494437LL
336 
337 int avpriv_init_elbg(int *points, int dim, int numpoints, int *codebook,
338  int numCB, int max_steps, int *closest_cb,
339  AVLFG *rand_state)
340 {
341  int i, k, ret = 0;
342 
343  if (numpoints > 24*numCB) {
344  /* ELBG is very costly for a big number of points. So if we have a lot
345  of them, get a good initial codebook to save on iterations */
346  int *temp_points = av_malloc_array(dim, (numpoints/8)*sizeof(int));
347  if (!temp_points)
348  return AVERROR(ENOMEM);
349  for (i=0; i<numpoints/8; i++) {
350  k = (i*BIG_PRIME) % numpoints;
351  memcpy(temp_points + i*dim, points + k*dim, dim*sizeof(int));
352  }
353 
354  ret = avpriv_init_elbg(temp_points, dim, numpoints / 8, codebook,
355  numCB, 2 * max_steps, closest_cb, rand_state);
356  if (ret < 0) {
357  av_freep(&temp_points);
358  return ret;
359  }
360  ret = avpriv_do_elbg(temp_points, dim, numpoints / 8, codebook,
361  numCB, 2 * max_steps, closest_cb, rand_state);
362  av_free(temp_points);
363 
364  } else // If not, initialize the codebook with random positions
365  for (i=0; i < numCB; i++)
366  memcpy(codebook + i*dim, points + ((i*BIG_PRIME)%numpoints)*dim,
367  dim*sizeof(int));
368  return ret;
369 }
370 
371 int avpriv_do_elbg(int *points, int dim, int numpoints, int *codebook,
372  int numCB, int max_steps, int *closest_cb,
373  AVLFG *rand_state)
374 {
375  int dist;
376  elbg_data elbg_d;
377  elbg_data *elbg = &elbg_d;
378  int i, j, k, last_error, steps = 0, ret = 0;
379  int *dist_cb = av_malloc_array(numpoints, sizeof(int));
380  int *size_part = av_malloc_array(numCB, sizeof(int));
381  cell *list_buffer = av_malloc_array(numpoints, sizeof(cell));
382  cell *free_cells;
383  int best_dist, best_idx = 0;
384 
385  elbg->error = INT_MAX;
386  elbg->dim = dim;
387  elbg->numCB = numCB;
388  elbg->codebook = codebook;
389  elbg->cells = av_malloc_array(numCB, sizeof(cell *));
390  elbg->utility = av_malloc_array(numCB, sizeof(int));
391  elbg->nearest_cb = closest_cb;
392  elbg->points = points;
393  elbg->utility_inc = av_malloc_array(numCB, sizeof(*elbg->utility_inc));
394  elbg->scratchbuf = av_malloc_array(5*dim, sizeof(int));
395 
396  if (!dist_cb || !size_part || !list_buffer || !elbg->cells ||
397  !elbg->utility || !elbg->utility_inc || !elbg->scratchbuf) {
398  ret = AVERROR(ENOMEM);
399  goto out;
400  }
401 
402  elbg->rand_state = rand_state;
403 
404  do {
405  free_cells = list_buffer;
406  last_error = elbg->error;
407  steps++;
408  memset(elbg->utility, 0, numCB*sizeof(int));
409  memset(elbg->cells, 0, numCB*sizeof(cell *));
410 
411  elbg->error = 0;
412 
413  /* This loop evaluate the actual Voronoi partition. It is the most
414  costly part of the algorithm. */
415  for (i=0; i < numpoints; i++) {
416  best_dist = distance_limited(elbg->points + i*elbg->dim, elbg->codebook + best_idx*elbg->dim, dim, INT_MAX);
417  for (k=0; k < elbg->numCB; k++) {
418  dist = distance_limited(elbg->points + i*elbg->dim, elbg->codebook + k*elbg->dim, dim, best_dist);
419  if (dist < best_dist) {
420  best_dist = dist;
421  best_idx = k;
422  }
423  }
424  elbg->nearest_cb[i] = best_idx;
425  dist_cb[i] = best_dist;
426  elbg->error += dist_cb[i];
427  elbg->utility[elbg->nearest_cb[i]] += dist_cb[i];
428  free_cells->index = i;
429  free_cells->next = elbg->cells[elbg->nearest_cb[i]];
430  elbg->cells[elbg->nearest_cb[i]] = free_cells;
431  free_cells++;
432  }
433 
434  do_shiftings(elbg);
435 
436  memset(size_part, 0, numCB*sizeof(int));
437 
438  memset(elbg->codebook, 0, elbg->numCB*dim*sizeof(int));
439 
440  for (i=0; i < numpoints; i++) {
441  size_part[elbg->nearest_cb[i]]++;
442  for (j=0; j < elbg->dim; j++)
443  elbg->codebook[elbg->nearest_cb[i]*elbg->dim + j] +=
444  elbg->points[i*elbg->dim + j];
445  }
446 
447  for (i=0; i < elbg->numCB; i++)
448  vect_division(elbg->codebook + i*elbg->dim,
449  elbg->codebook + i*elbg->dim, size_part[i], elbg->dim);
450 
451  } while(((last_error - elbg->error) > DELTA_ERR_MAX*elbg->error) &&
452  (steps < max_steps));
453 
454 out:
455  av_free(dist_cb);
456  av_free(size_part);
457  av_free(elbg->utility);
458  av_free(list_buffer);
459  av_free(elbg->cells);
460  av_free(elbg->utility_inc);
461  av_free(elbg->scratchbuf);
462  return ret;
463 }
error
static void error(const char *err)
Definition: target_bsf_fuzzer.c:29
vect_division
static void vect_division(int *res, int *vect, int div, int dim)
Definition: elbg.c:73
r
const char * r
Definition: vf_curves.c:114
AVERROR
Filter the word “frame” indicates either a video frame or a group of audio as stored in an AVFrame structure Format for each input and each output the list of supported formats For video that means pixel format For audio that means channel sample they are references to shared objects When the negotiation mechanism computes the intersection of the formats supported at each end of a all references to both lists are replaced with a reference to the intersection And when a single format is eventually chosen for a link amongst the remaining all references to the list are updated That means that if a filter requires that its input and output have the same format amongst a supported all it has to do is use a reference to the same list of formats query_formats can leave some formats unset and return AVERROR(EAGAIN) to cause the negotiation mechanism toagain later. That can be used by filters with complex requirements to use the format negotiated on one link to set the formats supported on another. Frame references ownership and permissions
elbg_data::cells
cell ** cells
Definition: elbg.c:52
BIG_PRIME
#define BIG_PRIME
Definition: elbg.c:335
out
FILE * out
Definition: movenc.c:54
update_utility_and_n_cb
static void update_utility_and_n_cb(elbg_data *elbg, int idx, int newutility)
Definition: elbg.c:250
cell_s
In the ELBG jargon, a cell is the set of points that are closest to a codebook entry.
Definition: elbg.c:39
elbg_data::rand_state
AVLFG * rand_state
Definition: elbg.c:57
b
#define b
Definition: input.c:41
max
#define max(a, b)
Definition: cuda_runtime.h:33
eval_error_cell
static int eval_error_cell(elbg_data *elbg, int *centroid, cell *cells)
Definition: elbg.c:84
avassert.h
elbg_data::codebook
int * codebook
Definition: elbg.c:51
elbg_data::dim
int dim
Definition: elbg.c:49
av_lfg_get
static unsigned int av_lfg_get(AVLFG *c)
Get the next random unsigned 32-bit number using an ALFG.
Definition: lfg.h:53
lfg.h
distance_limited
static int distance_limited(int *a, int *b, int dim, int limit)
Definition: elbg.c:61
DELTA_ERR_MAX
#define DELTA_ERR_MAX
Precision of the ELBG algorithm (as percentage error)
Definition: elbg.c:34
elbg_data::points
int * points
Definition: elbg.c:56
elbg.h
evaluate_utility_inc
static void evaluate_utility_inc(elbg_data *elbg)
Definition: elbg.c:237
elbg_data::nearest_cb
int * nearest_cb
Definition: elbg.c:55
NULL
#define NULL
Definition: coverity.c:32
avpriv_init_elbg
int avpriv_init_elbg(int *points, int dim, int numpoints, int *codebook, int numCB, int max_steps, int *closest_cb, AVLFG *rand_state)
Initialize the **codebook vector for the elbg algorithm.
Definition: elbg.c:337
ROUNDED_DIV
#define ROUNDED_DIV(a, b)
Definition: common.h:56
get_new_centroids
static void get_new_centroids(elbg_data *elbg, int huc, int *newcentroid_i, int *newcentroid_p)
Definition: elbg.c:173
index
int index
Definition: gxfenc.c:89
cell_s::index
int index
Definition: elbg.c:40
cell_s::next
struct cell_s * next
Definition: elbg.c:41
AVLFG
Context structure for the Lagged Fibonacci PRNG.
Definition: lfg.h:33
FFMAX
#define FFMAX(a, b)
Definition: common.h:94
elbg_data::scratchbuf
int * scratchbuf
Definition: elbg.c:58
FFMIN
#define FFMIN(a, b)
Definition: common.h:96
a
The reader does not expect b to be semantically here and if the code is changed by maybe adding a a division or other the signedness will almost certainly be mistaken To avoid this confusion a new type was SUINT is the C unsigned type but it holds a signed int to use the same example SUINT a
Definition: undefined.txt:41
try_shift_candidate
static void try_shift_candidate(elbg_data *elbg, int idx[3])
Evaluate if a shift lower the error.
Definition: elbg.c:266
av_assert2
#define av_assert2(cond)
assert() equivalent, that does lie in speed critical code.
Definition: avassert.h:64
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:269
av_malloc_array
#define av_malloc_array(a, b)
Definition: tableprint_vlc.h:32
common.h
elbg_data::error
int error
Definition: elbg.c:48
elbg_data
ELBG internal data.
Definition: elbg.c:47
elbg_data::numCB
int numCB
Definition: elbg.c:50
avcodec.h
dim
int dim
Definition: vorbis_enc_data.h:451
get_closest_codebook
static int get_closest_codebook(elbg_data *elbg, int index)
Definition: elbg.c:93
ret
ret
Definition: filter_design.txt:187
get_high_utility_cell
static int get_high_utility_cell(elbg_data *elbg)
Definition: elbg.c:107
elbg_data::utility_inc
int64_t * utility_inc
Definition: elbg.c:54
simple_lbg
static int simple_lbg(elbg_data *elbg, int dim, int *centroid[3], int newutility[3], int *points, cell *cells)
Implementation of the simple LBG algorithm for just two codebooks.
Definition: elbg.c:132
diff
static av_always_inline int diff(const uint32_t a, const uint32_t b)
Definition: vf_palettegen.c:136
av_free
#define av_free(p)
Definition: tableprint_vlc.h:34
elbg_data::utility
int * utility
Definition: elbg.c:53
av_freep
#define av_freep(p)
Definition: tableprint_vlc.h:35
do_shiftings
static void do_shiftings(elbg_data *elbg)
Implementation of the ELBG block.
Definition: elbg.c:316
int
int
Definition: ffmpeg_filter.c:192
avpriv_do_elbg
int avpriv_do_elbg(int *points, int dim, int numpoints, int *codebook, int numCB, int max_steps, int *closest_cb, AVLFG *rand_state)
Implementation of the Enhanced LBG Algorithm Based on the paper "Neural Networks 14:1219-1237" that c...
Definition: elbg.c:371
shift_codebook
static void shift_codebook(elbg_data *elbg, int *indexes, int *newcentroid[3])
Add the points in the low utility cell to its closest cell.
Definition: elbg.c:209
min
float min
Definition: vorbis_enc_data.h:456