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
32
33 #define DELTA_ERR_MAX 0.1 ///< Precision of the ELBG algorithm (as percentage error)
34
35 /**
36 * In the ELBG jargon, a cell is the set of points that are closest to a
37 * codebook entry. Not to be confused with a RoQ Video cell. */
41 } cell;
42
43 /**
44 * ELBG internal data
45 */
61
62 /* Sizes for the buffers above. Pointers without such a field
63 * are not allocated by us and only valid for the duration
64 * of a single call to avpriv_elbg_do(). */
73
75 {
80 return INT_MAX;
81 }
82
83 return dist;
84 }
85
87 {
89 if (div > 1)
92 else if (res != vect)
93 memcpy(res, vect,
dim*
sizeof(
int));
94
95 }
96
98 {
100 for (; cells; cells=cells->next)
102
104 }
105
107 {
108 int pick = 0;
109 for (
int i = 0, diff_min = INT_MAX;
i < elbg->
num_cb;
i++)
113 if (
diff < diff_min) {
116 }
117 }
118 return pick;
119 }
120
122 {
124 /* Using linear search, do binary if it ever turns to be speed critical */
126
129 } else {
132 }
133
136 }
137
139
141 }
142
143 /**
144 * Implementation of the simple LBG algorithm for just two codebooks
145 */
148 int *centroid[3],
149 int newutility[3],
150 int *points,
151 cell *cells)
152 {
154 int numpoints[2] = {0,0};
155 int *newcentroid[2] = {
158 };
159 cell *tempcell;
160
161 memset(newcentroid[0], 0, 2 *
dim *
sizeof(*newcentroid[0]));
162
163 newutility[0] =
164 newutility[1] = 0;
165
166 for (tempcell = cells; tempcell; tempcell=tempcell->next) {
169 numpoints[idx]++;
171 newcentroid[idx][
i] += points[tempcell->index*
dim +
i];
172 }
173
176
177 for (tempcell = cells; tempcell; tempcell=tempcell->next) {
180 int idx = dist[0] > dist[1];
181 newutility[idx] += dist[idx];
182 }
183
184 return newutility[0] + newutility[1];
185 }
186
188 int *newcentroid_p)
189 {
190 cell *tempcell;
191 int *
min = newcentroid_i;
192 int *
max = newcentroid_p;
194
195 for (
i=0;
i< elbg->
dim;
i++) {
198 }
199
200 for (tempcell = elbg->
cells[huc]; tempcell; tempcell = tempcell->next)
201 for(
i=0;
i<elbg->
dim;
i++) {
204 }
205
206 for (
i=0;
i<elbg->
dim;
i++) {
209 newcentroid_i[
i] = ni;
210 newcentroid_p[
i] = np;
211 }
212 }
213
214 /**
215 * Add the points in the low utility cell to its closest cell. Split the high
216 * utility cell, putting the separated points in the (now empty) low utility
217 * cell.
218 *
219 * @param elbg Internal elbg data
220 * @param indexes {luc, huc, cluc}
221 * @param newcentroid A vector with the position of the new centroids
222 */
224 int *newcentroid[3])
225 {
226 cell *tempdata;
227 cell **pp = &elbg->
cells[indexes[2]];
228
229 while(*pp)
230 pp= &(*pp)->next;
231
232 *pp = elbg->
cells[indexes[0]];
233
235 tempdata = elbg->
cells[indexes[1]];
237
238 while(tempdata) {
239 cell *tempcell2 = tempdata->next;
241 newcentroid[0], elbg->
dim, INT_MAX) >
243 newcentroid[1], elbg->
dim, INT_MAX);
244
245 tempdata->next = elbg->
cells[indexes[idx]];
246 elbg->
cells[indexes[idx]] = tempdata;
247 tempdata = tempcell2;
248 }
249 }
250
252 {
253 int64_t inc=0;
254
259 }
260 }
261
262
264 {
265 cell *tempcell;
266
267 elbg->
utility[idx] = newutility;
268 for (tempcell=elbg->
cells[idx]; tempcell; tempcell=tempcell->next)
270 }
271
272 /**
273 * Evaluate if a shift lower the error. If it does, call shift_codebooks
274 * and update elbg->error, elbg->utility and elbg->nearest_cb.
275 *
276 * @param elbg Internal elbg data
277 * @param idx {luc (low utility cell, huc (high utility cell), cluc (closest cell to low utility cell)}
278 */
280 {
281 int j, k, cont=0;
282 int64_t olderror=0, newerror;
283 int newutility[3];
284 int *newcentroid[3] = {
288 };
289 cell *tempcell;
290
291 for (j=0; j<3; j++)
292 olderror += elbg->
utility[idx[j]];
293
294 memset(newcentroid[2], 0, elbg->
dim*
sizeof(
int));
295
296 for (k=0; k<2; k++)
297 for (tempcell=elbg->
cells[idx[2*k]]; tempcell; tempcell=tempcell->next) {
298 cont++;
299 for (j=0; j<elbg->
dim; j++)
300 newcentroid[2][j] += elbg->
points[tempcell->index*elbg->
dim + j];
301 }
302
304
306
309
310 newerror = newutility[2];
311
313 elbg->
cells[idx[1]]);
314
315 if (olderror > newerror) {
317
318 elbg->
error += newerror - olderror;
319
320 for (j=0; j<3; j++)
322
324 }
325 }
326
327 /**
328 * Implementation of the ELBG block
329 */
331 {
332 int idx[3];
333
335
336 for (idx[0]=0; idx[0] < elbg->
num_cb; idx[0]++)
339 return;
340
343
344 if (idx[1] != idx[0] && idx[1] != idx[2])
346 }
347 }
348
350 int max_steps)
351 {
352 int *const size_part = elbg->size_part;
354 int best_idx = 0;
355 int64_t last_error;
356
357 elbg->error = INT64_MAX;
358 elbg->points = points;
359
360 do {
361 cell *free_cells = elbg->cell_buffer;
362 last_error = elbg->error;
363 steps++;
364 memset(elbg->utility, 0, elbg->num_cb * sizeof(*elbg->utility));
365 memset(elbg->cells, 0, elbg->num_cb * sizeof(*elbg->cells));
366
367 elbg->error = 0;
368
369 /* This loop evaluate the actual Voronoi partition. It is the most
370 costly part of the algorithm. */
371 for (
i=0;
i < numpoints;
i++) {
373 elbg->codebook + best_idx * elbg->dim,
374 elbg->dim, INT_MAX);
375 for (int k = 0; k < elbg->num_cb; k++) {
377 elbg->codebook + k * elbg->dim,
378 elbg->dim, best_dist);
379 if (dist < best_dist) {
380 best_dist = dist;
381 best_idx = k;
382 }
383 }
384 elbg->nearest_cb[
i] = best_idx;
385 elbg->error += best_dist;
386 elbg->utility[elbg->nearest_cb[
i]] += best_dist;
387 free_cells->index =
i;
388 free_cells->next = elbg->cells[elbg->nearest_cb[
i]];
389 elbg->cells[elbg->nearest_cb[
i]] = free_cells;
390 free_cells++;
391 }
392
394
395 memset(size_part, 0, elbg->num_cb * sizeof(*size_part));
396
397 memset(elbg->codebook, 0, elbg->num_cb * elbg->dim * sizeof(*elbg->codebook));
398
399 for (
i=0;
i < numpoints;
i++) {
400 size_part[elbg->nearest_cb[
i]]++;
401 for (j=0; j < elbg->dim; j++)
402 elbg->codebook[elbg->nearest_cb[
i]*elbg->dim + j] +=
403 elbg->points[
i*elbg->dim + j];
404 }
405
406 for (
int i = 0;
i < elbg->num_cb;
i++)
408 elbg->codebook +
i*elbg->dim, size_part[
i], elbg->dim);
409
410 }
while(((last_error - elbg->error) >
DELTA_ERR_MAX*elbg->error) &&
411 (steps < max_steps));
412 }
413
414 #define BIG_PRIME 433494437LL
415
416 /**
417 * Initialize the codebook vector for the elbg algorithm.
418 * If numpoints <= 24 * num_cb this function fills codebook with random numbers.
419 * If not, it calls do_elbg for a (smaller) random sample of the points in
420 * points.
421 */
423 int numpoints, int max_steps)
424 {
426
427 if (numpoints > 24LL * elbg->num_cb) {
428 /* ELBG is very costly for a big number of points. So if we have a lot
429 of them, get a good initial codebook to save on iterations */
430 for (
int i = 0;
i < numpoints / 8;
i++) {
432 memcpy(temp_points +
i*
dim, points + k*
dim,
dim *
sizeof(*temp_points));
433 }
434
435 /* If anything is changed in the recursion parameters,
436 * the allocated size of temp_points will also need to be updated. */
437 init_elbg(elbg, temp_points, temp_points + numpoints / 8 *
dim,
438 numpoints / 8, 2 * max_steps);
439 do_elbg(elbg, temp_points, numpoints / 8, 2 * max_steps);
440 } else // If not, initialize the codebook with random positions
441 for (
int i = 0;
i < elbg->num_cb;
i++)
443 dim *
sizeof(*elbg->codebook));
444 }
445
447 int *
codebook,
int num_cb,
int max_steps,
448 int *closest_cb,
AVLFG *rand_state, uintptr_t
flags)
449 {
451
452 if (!elbg)
454 *elbgp = elbg;
455
457 elbg->rand_state = rand_state;
459 elbg->num_cb = num_cb;
461
462 #define ALLOCATE_IF_NECESSARY(field, new_elements, multiplicator) \
463 if (elbg->field ## _allocated < new_elements) { \
464 av_freep(&elbg->field); \
465 elbg->field = av_malloc_array(new_elements, \
466 multiplicator * sizeof(*elbg->field)); \
467 if (!elbg->field) { \
468 elbg->field ## _allocated = 0; \
469 return AVERROR(ENOMEM); \
470 } \
471 elbg->field ## _allocated = new_elements; \
472 }
473 /* Allocating the buffers for do_elbg() here once relies
474 * on their size being always the same even when do_elbg()
475 * is called from init_elbg(). It also relies on do_elbg()
476 * never calling itself recursively. */
483 if (numpoints > 24LL * elbg->num_cb) {
484 /* The first step in the recursion in init_elbg() needs a buffer with
485 * (numpoints / 8) * dim elements; the next step needs numpoints / 8 / 8
486 * * dim elements etc. The geometric series leads to an upper bound of
487 * numpoints / 8 * 8 / 7 * dim elements. */
488 uint64_t prod =
dim * (uint64_t)(numpoints / 7
U);
489 if (prod > INT_MAX)
492 }
493
494 init_elbg(elbg, points, elbg->temp_points, numpoints, max_steps);
495 do_elbg (elbg, points, numpoints, max_steps);
496 return 0;
497 }
498
500 {
502 if (!elbg)
503 return;
504
512
514 }