2 * Copyright (C) 2007 Vitor Sessak <vitor1001@gmail.com>
4 * This file is part of FFmpeg.
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.
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.
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
23 * Codebook Generator using the ELBG algorithm
28 #include "libavutil/avassert.h"
29 #include "libavutil/common.h"
30 #include "libavutil/lfg.h"
34 #define DELTA_ERR_MAX 0.1 ///< Precision of the ELBG algorithm (as percentual error)
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
{
61 static inline int distance_limited(int *a
, int *b
, int dim
, int limit
)
64 for (i
=0; i
<dim
; i
++) {
65 dist
+= (a
[i
] - b
[i
])*(a
[i
] - b
[i
]);
73 static inline void vect_division(int *res
, int *vect
, int div
, int dim
)
78 res
[i
] = ROUNDED_DIV(vect
[i
],div
);
80 memcpy(res
, vect
, dim
*sizeof(int));
84 static int eval_error_cell(elbg_data
*elbg
, int *centroid
, cell
*cells
)
87 for (; cells
; cells
=cells
->next
)
88 error
+= distance_limited(centroid
, elbg
->points
+ cells
->index
*elbg
->dim
, elbg
->dim
, INT_MAX
);
93 static int get_closest_codebook(elbg_data
*elbg
, int index
)
95 int i
, pick
=0, diff
, diff_min
= INT_MAX
;
96 for (i
=0; i
<elbg
->numCB
; i
++)
98 diff
= distance_limited(elbg
->codebook
+ i
*elbg
->dim
, elbg
->codebook
+ index
*elbg
->dim
, elbg
->dim
, diff_min
);
99 if (diff
< diff_min
) {
107 static int get_high_utility_cell(elbg_data
*elbg
)
110 /* Using linear search, do binary if it ever turns to be speed critical */
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;
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;
120 while (elbg
->utility_inc
[i
] < r
) {
124 av_assert2(elbg
->cells
[i
]);
130 * Implementation of the simple LBG algorithm for just two codebooks
132 static int simple_lbg(elbg_data
*elbg
,
140 int numpoints
[2] = {0,0};
141 int *newcentroid
[2] = {
142 elbg
->scratchbuf
+ 3*dim
,
143 elbg
->scratchbuf
+ 4*dim
147 memset(newcentroid
[0], 0, 2 * dim
* sizeof(*newcentroid
[0]));
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
);
156 for (i
=0; i
<dim
; i
++)
157 newcentroid
[idx
][i
] += points
[tempcell
->index
*dim
+ i
];
160 vect_division(centroid
[0], newcentroid
[0], numpoints
[0], dim
);
161 vect_division(centroid
[1], newcentroid
[1], numpoints
[1], dim
);
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
];
170 return newutility
[0] + newutility
[1];
173 static void get_new_centroids(elbg_data
*elbg
, int huc
, int *newcentroid_i
,
177 int *min
= newcentroid_i
;
178 int *max
= newcentroid_p
;
181 for (i
=0; i
< elbg
->dim
; i
++) {
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
]);
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
;
201 * Add the points in the low utility cell to its closest cell. Split the high
202 * utility cell, putting the separate points in the (now empty) low utility
205 * @param elbg Internal elbg data
206 * @param indexes {luc, huc, cluc}
207 * @param newcentroid A vector with the position of the new centroids
209 static void shift_codebook(elbg_data
*elbg
, int *indexes
,
213 cell
**pp
= &elbg
->cells
[indexes
[2]];
218 *pp
= elbg
->cells
[indexes
[0]];
220 elbg
->cells
[indexes
[0]] = NULL
;
221 tempdata
= elbg
->cells
[indexes
[1]];
222 elbg
->cells
[indexes
[1]] = NULL
;
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
);
231 tempdata
->next
= elbg
->cells
[indexes
[idx
]];
232 elbg
->cells
[indexes
[idx
]] = tempdata
;
233 tempdata
= tempcell2
;
237 static void evaluate_utility_inc(elbg_data
*elbg
)
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
;
250 static void update_utility_and_n_cb(elbg_data
*elbg
, int idx
, int newutility
)
254 elbg
->utility
[idx
] = newutility
;
255 for (tempcell
=elbg
->cells
[idx
]; tempcell
; tempcell
=tempcell
->next
)
256 elbg
->nearest_cb
[tempcell
->index
] = idx
;
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.
263 * @param elbg Internal elbg data
264 * @param idx {luc (low utility cell, huc (high utility cell), cluc (closest cell to low utility cell)}
266 static void try_shift_candidate(elbg_data
*elbg
, int idx
[3])
268 int j
, k
, olderror
=0, newerror
, cont
=0;
270 int *newcentroid
[3] = {
272 elbg
->scratchbuf
+ elbg
->dim
,
273 elbg
->scratchbuf
+ 2*elbg
->dim
278 olderror
+= elbg
->utility
[idx
[j
]];
280 memset(newcentroid
[2], 0, elbg
->dim
*sizeof(int));
283 for (tempcell
=elbg
->cells
[idx
[2*k
]]; tempcell
; tempcell
=tempcell
->next
) {
285 for (j
=0; j
<elbg
->dim
; j
++)
286 newcentroid
[2][j
] += elbg
->points
[tempcell
->index
*elbg
->dim
+ j
];
289 vect_division(newcentroid
[2], newcentroid
[2], cont
, elbg
->dim
);
291 get_new_centroids(elbg
, idx
[1], newcentroid
[0], newcentroid
[1]);
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]]);
296 newerror
= newutility
[2];
298 newerror
+= simple_lbg(elbg
, elbg
->dim
, newcentroid
, newutility
, elbg
->points
,
299 elbg
->cells
[idx
[1]]);
301 if (olderror
> newerror
) {
302 shift_codebook(elbg
, idx
, newcentroid
);
304 elbg
->error
+= newerror
- olderror
;
307 update_utility_and_n_cb(elbg
, idx
[j
], newutility
[j
]);
309 evaluate_utility_inc(elbg
);
314 * Implementation of the ELBG block
316 static void do_shiftings(elbg_data
*elbg
)
320 evaluate_utility_inc(elbg
);
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)
327 idx
[1] = get_high_utility_cell(elbg
);
328 idx
[2] = get_closest_codebook(elbg
, idx
[0]);
330 if (idx
[1] != idx
[0] && idx
[1] != idx
[2])
331 try_shift_candidate(elbg
, idx
);
335 #define BIG_PRIME 433494437LL
337 void avpriv_init_elbg(int *points
, int dim
, int numpoints
, int *codebook
,
338 int numCB
, int max_steps
, int *closest_cb
,
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 for (i
=0; i
<numpoints
/8; i
++) {
348 k
= (i
*BIG_PRIME
) % numpoints
;
349 memcpy(temp_points
+ i
*dim
, points
+ k
*dim
, dim
*sizeof(int));
352 avpriv_init_elbg(temp_points
, dim
, numpoints
/8, codebook
, numCB
, 2*max_steps
, closest_cb
, rand_state
);
353 avpriv_do_elbg(temp_points
, dim
, numpoints
/8, codebook
, numCB
, 2*max_steps
, closest_cb
, rand_state
);
355 av_free(temp_points
);
357 } else // If not, initialize the codebook with random positions
358 for (i
=0; i
< numCB
; i
++)
359 memcpy(codebook
+ i
*dim
, points
+ ((i
*BIG_PRIME
)%numpoints
)*dim
,
364 void avpriv_do_elbg(int *points
, int dim
, int numpoints
, int *codebook
,
365 int numCB
, int max_steps
, int *closest_cb
,
370 elbg_data
*elbg
= &elbg_d
;
371 int i
, j
, k
, last_error
, steps
=0;
372 int *dist_cb
= av_malloc_array(numpoints
, sizeof(int));
373 int *size_part
= av_malloc_array(numCB
, sizeof(int));
374 cell
*list_buffer
= av_malloc_array(numpoints
, sizeof(cell
));
376 int best_dist
, best_idx
= 0;
378 elbg
->error
= INT_MAX
;
381 elbg
->codebook
= codebook
;
382 elbg
->cells
= av_malloc_array(numCB
, sizeof(cell
*));
383 elbg
->utility
= av_malloc_array(numCB
, sizeof(int));
384 elbg
->nearest_cb
= closest_cb
;
385 elbg
->points
= points
;
386 elbg
->utility_inc
= av_malloc_array(numCB
, sizeof(*elbg
->utility_inc
));
387 elbg
->scratchbuf
= av_malloc_array(5*dim
, sizeof(int));
389 elbg
->rand_state
= rand_state
;
392 free_cells
= list_buffer
;
393 last_error
= elbg
->error
;
395 memset(elbg
->utility
, 0, numCB
*sizeof(int));
396 memset(elbg
->cells
, 0, numCB
*sizeof(cell
*));
400 /* This loop evaluate the actual Voronoi partition. It is the most
401 costly part of the algorithm. */
402 for (i
=0; i
< numpoints
; i
++) {
403 best_dist
= distance_limited(elbg
->points
+ i
*elbg
->dim
, elbg
->codebook
+ best_idx
*elbg
->dim
, dim
, INT_MAX
);
404 for (k
=0; k
< elbg
->numCB
; k
++) {
405 dist
= distance_limited(elbg
->points
+ i
*elbg
->dim
, elbg
->codebook
+ k
*elbg
->dim
, dim
, best_dist
);
406 if (dist
< best_dist
) {
411 elbg
->nearest_cb
[i
] = best_idx
;
412 dist_cb
[i
] = best_dist
;
413 elbg
->error
+= dist_cb
[i
];
414 elbg
->utility
[elbg
->nearest_cb
[i
]] += dist_cb
[i
];
415 free_cells
->index
= i
;
416 free_cells
->next
= elbg
->cells
[elbg
->nearest_cb
[i
]];
417 elbg
->cells
[elbg
->nearest_cb
[i
]] = free_cells
;
423 memset(size_part
, 0, numCB
*sizeof(int));
425 memset(elbg
->codebook
, 0, elbg
->numCB
*dim
*sizeof(int));
427 for (i
=0; i
< numpoints
; i
++) {
428 size_part
[elbg
->nearest_cb
[i
]]++;
429 for (j
=0; j
< elbg
->dim
; j
++)
430 elbg
->codebook
[elbg
->nearest_cb
[i
]*elbg
->dim
+ j
] +=
431 elbg
->points
[i
*elbg
->dim
+ j
];
434 for (i
=0; i
< elbg
->numCB
; i
++)
435 vect_division(elbg
->codebook
+ i
*elbg
->dim
,
436 elbg
->codebook
+ i
*elbg
->dim
, size_part
[i
], elbg
->dim
);
438 } while(((last_error
- elbg
->error
) > DELTA_ERR_MAX
*elbg
->error
) &&
439 (steps
< max_steps
));
443 av_free(elbg
->utility
);
444 av_free(list_buffer
);
445 av_free(elbg
->cells
);
446 av_free(elbg
->utility_inc
);
447 av_free(elbg
->scratchbuf
);