/src/aom/av1/encoder/global_motion.c
Line | Count | Source (jump to first uncovered line) |
1 | | /* |
2 | | * Copyright (c) 2016, Alliance for Open Media. All rights reserved |
3 | | * |
4 | | * This source code is subject to the terms of the BSD 2 Clause License and |
5 | | * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License |
6 | | * was not distributed with this source code in the LICENSE file, you can |
7 | | * obtain it at www.aomedia.org/license/software. If the Alliance for Open |
8 | | * Media Patent License 1.0 was not distributed with this source code in the |
9 | | * PATENTS file, you can obtain it at www.aomedia.org/license/patent. |
10 | | */ |
11 | | |
12 | | #include <stdio.h> |
13 | | #include <stdlib.h> |
14 | | #include <memory.h> |
15 | | #include <math.h> |
16 | | #include <assert.h> |
17 | | |
18 | | #include "config/aom_dsp_rtcd.h" |
19 | | |
20 | | #include "av1/encoder/global_motion.h" |
21 | | |
22 | | #include "av1/common/convolve.h" |
23 | | #include "av1/common/resize.h" |
24 | | #include "av1/common/warped_motion.h" |
25 | | |
26 | | #include "av1/encoder/segmentation.h" |
27 | | #include "av1/encoder/corner_detect.h" |
28 | | #include "av1/encoder/corner_match.h" |
29 | | #include "av1/encoder/ransac.h" |
30 | | |
31 | 0 | #define MIN_INLIER_PROB 0.1 |
32 | | |
33 | 0 | #define MIN_TRANS_THRESH (1 * GM_TRANS_DECODE_FACTOR) |
34 | | |
35 | | // Border over which to compute the global motion |
36 | 0 | #define ERRORADV_BORDER 0 |
37 | | |
38 | | // Number of pyramid levels in disflow computation |
39 | 0 | #define N_LEVELS 2 |
40 | | // Size of square patches in the disflow dense grid |
41 | 0 | #define PATCH_SIZE 8 |
42 | | // Center point of square patch |
43 | 0 | #define PATCH_CENTER ((PATCH_SIZE + 1) >> 1) |
44 | | // Step size between patches, lower value means greater patch overlap |
45 | 0 | #define PATCH_STEP 1 |
46 | | // Minimum size of border padding for disflow |
47 | | #define MIN_PAD 7 |
48 | | // Warp error convergence threshold for disflow |
49 | 0 | #define DISFLOW_ERROR_TR 0.01 |
50 | | // Max number of iterations if warp convergence is not found |
51 | 0 | #define DISFLOW_MAX_ITR 10 |
52 | | |
53 | | // Struct for an image pyramid |
54 | | typedef struct { |
55 | | int n_levels; |
56 | | int pad_size; |
57 | | int has_gradient; |
58 | | int widths[N_LEVELS]; |
59 | | int heights[N_LEVELS]; |
60 | | int strides[N_LEVELS]; |
61 | | int level_loc[N_LEVELS]; |
62 | | unsigned char *level_buffer; |
63 | | double *level_dx_buffer; |
64 | | double *level_dy_buffer; |
65 | | } ImagePyramid; |
66 | | |
67 | 0 | int av1_is_enough_erroradvantage(double best_erroradvantage, int params_cost) { |
68 | 0 | return best_erroradvantage < erroradv_tr && |
69 | 0 | best_erroradvantage * params_cost < erroradv_prod_tr; |
70 | 0 | } |
71 | | |
72 | 0 | static void convert_to_params(const double *params, int32_t *model) { |
73 | 0 | int i; |
74 | 0 | int alpha_present = 0; |
75 | 0 | model[0] = (int32_t)floor(params[0] * (1 << GM_TRANS_PREC_BITS) + 0.5); |
76 | 0 | model[1] = (int32_t)floor(params[1] * (1 << GM_TRANS_PREC_BITS) + 0.5); |
77 | 0 | model[0] = (int32_t)clamp(model[0], GM_TRANS_MIN, GM_TRANS_MAX) * |
78 | 0 | GM_TRANS_DECODE_FACTOR; |
79 | 0 | model[1] = (int32_t)clamp(model[1], GM_TRANS_MIN, GM_TRANS_MAX) * |
80 | 0 | GM_TRANS_DECODE_FACTOR; |
81 | |
|
82 | 0 | for (i = 2; i < 6; ++i) { |
83 | 0 | const int diag_value = ((i == 2 || i == 5) ? (1 << GM_ALPHA_PREC_BITS) : 0); |
84 | 0 | model[i] = (int32_t)floor(params[i] * (1 << GM_ALPHA_PREC_BITS) + 0.5); |
85 | 0 | model[i] = |
86 | 0 | (int32_t)clamp(model[i] - diag_value, GM_ALPHA_MIN, GM_ALPHA_MAX); |
87 | 0 | alpha_present |= (model[i] != 0); |
88 | 0 | model[i] = (model[i] + diag_value) * GM_ALPHA_DECODE_FACTOR; |
89 | 0 | } |
90 | 0 | for (; i < 8; ++i) { |
91 | 0 | model[i] = (int32_t)floor(params[i] * (1 << GM_ROW3HOMO_PREC_BITS) + 0.5); |
92 | 0 | model[i] = (int32_t)clamp(model[i], GM_ROW3HOMO_MIN, GM_ROW3HOMO_MAX) * |
93 | 0 | GM_ROW3HOMO_DECODE_FACTOR; |
94 | 0 | alpha_present |= (model[i] != 0); |
95 | 0 | } |
96 | |
|
97 | 0 | if (!alpha_present) { |
98 | 0 | if (abs(model[0]) < MIN_TRANS_THRESH && abs(model[1]) < MIN_TRANS_THRESH) { |
99 | 0 | model[0] = 0; |
100 | 0 | model[1] = 0; |
101 | 0 | } |
102 | 0 | } |
103 | 0 | } |
104 | | |
105 | | void av1_convert_model_to_params(const double *params, |
106 | 0 | WarpedMotionParams *model) { |
107 | 0 | convert_to_params(params, model->wmmat); |
108 | 0 | model->wmtype = get_wmtype(model); |
109 | 0 | model->invalid = 0; |
110 | 0 | } |
111 | | |
112 | | // Adds some offset to a global motion parameter and handles |
113 | | // all of the necessary precision shifts, clamping, and |
114 | | // zero-centering. |
115 | | static int32_t add_param_offset(int param_index, int32_t param_value, |
116 | 0 | int32_t offset) { |
117 | 0 | const int scale_vals[3] = { GM_TRANS_PREC_DIFF, GM_ALPHA_PREC_DIFF, |
118 | 0 | GM_ROW3HOMO_PREC_DIFF }; |
119 | 0 | const int clamp_vals[3] = { GM_TRANS_MAX, GM_ALPHA_MAX, GM_ROW3HOMO_MAX }; |
120 | | // type of param: 0 - translation, 1 - affine, 2 - homography |
121 | 0 | const int param_type = (param_index < 2 ? 0 : (param_index < 6 ? 1 : 2)); |
122 | 0 | const int is_one_centered = (param_index == 2 || param_index == 5); |
123 | | |
124 | | // Make parameter zero-centered and offset the shift that was done to make |
125 | | // it compatible with the warped model |
126 | 0 | param_value = (param_value - (is_one_centered << WARPEDMODEL_PREC_BITS)) >> |
127 | 0 | scale_vals[param_type]; |
128 | | // Add desired offset to the rescaled/zero-centered parameter |
129 | 0 | param_value += offset; |
130 | | // Clamp the parameter so it does not overflow the number of bits allotted |
131 | | // to it in the bitstream |
132 | 0 | param_value = (int32_t)clamp(param_value, -clamp_vals[param_type], |
133 | 0 | clamp_vals[param_type]); |
134 | | // Rescale the parameter to WARPEDMODEL_PRECISION_BITS so it is compatible |
135 | | // with the warped motion library |
136 | 0 | param_value *= (1 << scale_vals[param_type]); |
137 | | |
138 | | // Undo the zero-centering step if necessary |
139 | 0 | return param_value + (is_one_centered << WARPEDMODEL_PREC_BITS); |
140 | 0 | } |
141 | | |
142 | 0 | static void force_wmtype(WarpedMotionParams *wm, TransformationType wmtype) { |
143 | 0 | switch (wmtype) { |
144 | 0 | case IDENTITY: |
145 | 0 | wm->wmmat[0] = 0; |
146 | 0 | wm->wmmat[1] = 0; |
147 | 0 | AOM_FALLTHROUGH_INTENDED; |
148 | 0 | case TRANSLATION: |
149 | 0 | wm->wmmat[2] = 1 << WARPEDMODEL_PREC_BITS; |
150 | 0 | wm->wmmat[3] = 0; |
151 | 0 | AOM_FALLTHROUGH_INTENDED; |
152 | 0 | case ROTZOOM: |
153 | 0 | wm->wmmat[4] = -wm->wmmat[3]; |
154 | 0 | wm->wmmat[5] = wm->wmmat[2]; |
155 | 0 | AOM_FALLTHROUGH_INTENDED; |
156 | 0 | case AFFINE: break; |
157 | 0 | default: assert(0); |
158 | 0 | } |
159 | 0 | wm->wmtype = wmtype; |
160 | 0 | } |
161 | | |
162 | | #if CONFIG_AV1_HIGHBITDEPTH |
163 | | static int64_t highbd_warp_error( |
164 | | WarpedMotionParams *wm, const uint16_t *const ref, int width, int height, |
165 | | int stride, const uint16_t *const dst, int p_col, int p_row, int p_width, |
166 | | int p_height, int p_stride, int subsampling_x, int subsampling_y, int bd, |
167 | 0 | int64_t best_error, uint8_t *segment_map, int segment_map_stride) { |
168 | 0 | int64_t gm_sumerr = 0; |
169 | 0 | const int error_bsize_w = AOMMIN(p_width, WARP_ERROR_BLOCK); |
170 | 0 | const int error_bsize_h = AOMMIN(p_height, WARP_ERROR_BLOCK); |
171 | 0 | uint16_t tmp[WARP_ERROR_BLOCK * WARP_ERROR_BLOCK]; |
172 | |
|
173 | 0 | ConvolveParams conv_params = get_conv_params(0, 0, bd); |
174 | 0 | conv_params.use_dist_wtd_comp_avg = 0; |
175 | 0 | for (int i = p_row; i < p_row + p_height; i += WARP_ERROR_BLOCK) { |
176 | 0 | for (int j = p_col; j < p_col + p_width; j += WARP_ERROR_BLOCK) { |
177 | 0 | int seg_x = j >> WARP_ERROR_BLOCK_LOG; |
178 | 0 | int seg_y = i >> WARP_ERROR_BLOCK_LOG; |
179 | | // Only compute the error if this block contains inliers from the motion |
180 | | // model |
181 | 0 | if (!segment_map[seg_y * segment_map_stride + seg_x]) continue; |
182 | | // avoid warping extra 8x8 blocks in the padded region of the frame |
183 | | // when p_width and p_height are not multiples of WARP_ERROR_BLOCK |
184 | 0 | const int warp_w = AOMMIN(error_bsize_w, p_col + p_width - j); |
185 | 0 | const int warp_h = AOMMIN(error_bsize_h, p_row + p_height - i); |
186 | 0 | highbd_warp_plane(wm, ref, width, height, stride, tmp, j, i, warp_w, |
187 | 0 | warp_h, WARP_ERROR_BLOCK, subsampling_x, subsampling_y, |
188 | 0 | bd, &conv_params); |
189 | 0 | gm_sumerr += av1_calc_highbd_frame_error(tmp, WARP_ERROR_BLOCK, |
190 | 0 | dst + j + i * p_stride, warp_w, |
191 | 0 | warp_h, p_stride, bd); |
192 | 0 | if (gm_sumerr > best_error) return INT64_MAX; |
193 | 0 | } |
194 | 0 | } |
195 | 0 | return gm_sumerr; |
196 | 0 | } |
197 | | #endif |
198 | | |
199 | | static int64_t warp_error(WarpedMotionParams *wm, const uint8_t *const ref, |
200 | | int width, int height, int stride, |
201 | | const uint8_t *const dst, int p_col, int p_row, |
202 | | int p_width, int p_height, int p_stride, |
203 | | int subsampling_x, int subsampling_y, |
204 | | int64_t best_error, uint8_t *segment_map, |
205 | 0 | int segment_map_stride) { |
206 | 0 | int64_t gm_sumerr = 0; |
207 | 0 | int warp_w, warp_h; |
208 | 0 | const int error_bsize_w = AOMMIN(p_width, WARP_ERROR_BLOCK); |
209 | 0 | const int error_bsize_h = AOMMIN(p_height, WARP_ERROR_BLOCK); |
210 | 0 | DECLARE_ALIGNED(16, uint8_t, tmp[WARP_ERROR_BLOCK * WARP_ERROR_BLOCK]); |
211 | 0 | ConvolveParams conv_params = get_conv_params(0, 0, 8); |
212 | 0 | conv_params.use_dist_wtd_comp_avg = 0; |
213 | |
|
214 | 0 | for (int i = p_row; i < p_row + p_height; i += WARP_ERROR_BLOCK) { |
215 | 0 | for (int j = p_col; j < p_col + p_width; j += WARP_ERROR_BLOCK) { |
216 | 0 | int seg_x = j >> WARP_ERROR_BLOCK_LOG; |
217 | 0 | int seg_y = i >> WARP_ERROR_BLOCK_LOG; |
218 | | // Only compute the error if this block contains inliers from the motion |
219 | | // model |
220 | 0 | if (!segment_map[seg_y * segment_map_stride + seg_x]) continue; |
221 | | // avoid warping extra 8x8 blocks in the padded region of the frame |
222 | | // when p_width and p_height are not multiples of WARP_ERROR_BLOCK |
223 | 0 | warp_w = AOMMIN(error_bsize_w, p_col + p_width - j); |
224 | 0 | warp_h = AOMMIN(error_bsize_h, p_row + p_height - i); |
225 | 0 | warp_plane(wm, ref, width, height, stride, tmp, j, i, warp_w, warp_h, |
226 | 0 | WARP_ERROR_BLOCK, subsampling_x, subsampling_y, &conv_params); |
227 | |
|
228 | 0 | gm_sumerr += |
229 | 0 | av1_calc_frame_error(tmp, WARP_ERROR_BLOCK, dst + j + i * p_stride, |
230 | 0 | warp_w, warp_h, p_stride); |
231 | 0 | if (gm_sumerr > best_error) return INT64_MAX; |
232 | 0 | } |
233 | 0 | } |
234 | 0 | return gm_sumerr; |
235 | 0 | } |
236 | | |
237 | | int64_t av1_warp_error(WarpedMotionParams *wm, int use_hbd, int bd, |
238 | | const uint8_t *ref, int width, int height, int stride, |
239 | | uint8_t *dst, int p_col, int p_row, int p_width, |
240 | | int p_height, int p_stride, int subsampling_x, |
241 | | int subsampling_y, int64_t best_error, |
242 | 0 | uint8_t *segment_map, int segment_map_stride) { |
243 | 0 | if (wm->wmtype <= AFFINE) |
244 | 0 | if (!av1_get_shear_params(wm)) return INT64_MAX; |
245 | 0 | #if CONFIG_AV1_HIGHBITDEPTH |
246 | 0 | if (use_hbd) |
247 | 0 | return highbd_warp_error(wm, CONVERT_TO_SHORTPTR(ref), width, height, |
248 | 0 | stride, CONVERT_TO_SHORTPTR(dst), p_col, p_row, |
249 | 0 | p_width, p_height, p_stride, subsampling_x, |
250 | 0 | subsampling_y, bd, best_error, segment_map, |
251 | 0 | segment_map_stride); |
252 | 0 | #endif |
253 | 0 | (void)use_hbd; |
254 | 0 | (void)bd; |
255 | 0 | return warp_error(wm, ref, width, height, stride, dst, p_col, p_row, p_width, |
256 | 0 | p_height, p_stride, subsampling_x, subsampling_y, |
257 | 0 | best_error, segment_map, segment_map_stride); |
258 | 0 | } |
259 | | |
260 | | // Factors used to calculate the thresholds for av1_warp_error |
261 | | static double thresh_factors[GM_REFINEMENT_COUNT] = { 1.25, 1.20, 1.15, 1.10, |
262 | | 1.05 }; |
263 | | |
264 | | static INLINE int64_t calc_approx_erroradv_threshold( |
265 | 0 | double scaling_factor, int64_t erroradv_threshold) { |
266 | 0 | return erroradv_threshold < |
267 | 0 | (int64_t)(((double)INT64_MAX / scaling_factor) + 0.5) |
268 | 0 | ? (int64_t)(scaling_factor * erroradv_threshold + 0.5) |
269 | 0 | : INT64_MAX; |
270 | 0 | } |
271 | | |
272 | | int64_t av1_refine_integerized_param( |
273 | | WarpedMotionParams *wm, TransformationType wmtype, int use_hbd, int bd, |
274 | | uint8_t *ref, int r_width, int r_height, int r_stride, uint8_t *dst, |
275 | | int d_width, int d_height, int d_stride, int n_refinements, |
276 | | int64_t best_frame_error, uint8_t *segment_map, int segment_map_stride, |
277 | 0 | int64_t erroradv_threshold) { |
278 | 0 | static const int max_trans_model_params[TRANS_TYPES] = { 0, 2, 4, 6 }; |
279 | 0 | const int border = ERRORADV_BORDER; |
280 | 0 | int i = 0, p; |
281 | 0 | int n_params = max_trans_model_params[wmtype]; |
282 | 0 | int32_t *param_mat = wm->wmmat; |
283 | 0 | int64_t step_error, best_error; |
284 | 0 | int32_t step; |
285 | 0 | int32_t *param; |
286 | 0 | int32_t curr_param; |
287 | 0 | int32_t best_param; |
288 | |
|
289 | 0 | force_wmtype(wm, wmtype); |
290 | 0 | best_error = |
291 | 0 | av1_warp_error(wm, use_hbd, bd, ref, r_width, r_height, r_stride, |
292 | 0 | dst + border * d_stride + border, border, border, |
293 | 0 | d_width - 2 * border, d_height - 2 * border, d_stride, 0, |
294 | 0 | 0, best_frame_error, segment_map, segment_map_stride); |
295 | 0 | best_error = AOMMIN(best_error, best_frame_error); |
296 | 0 | step = 1 << (n_refinements - 1); |
297 | 0 | for (i = 0; i < n_refinements; i++, step >>= 1) { |
298 | 0 | int64_t error_adv_thresh = |
299 | 0 | calc_approx_erroradv_threshold(thresh_factors[i], erroradv_threshold); |
300 | 0 | for (p = 0; p < n_params; ++p) { |
301 | 0 | int step_dir = 0; |
302 | | // Skip searches for parameters that are forced to be 0 |
303 | 0 | param = param_mat + p; |
304 | 0 | curr_param = *param; |
305 | 0 | best_param = curr_param; |
306 | | // look to the left |
307 | 0 | *param = add_param_offset(p, curr_param, -step); |
308 | 0 | step_error = |
309 | 0 | av1_warp_error(wm, use_hbd, bd, ref, r_width, r_height, r_stride, |
310 | 0 | dst + border * d_stride + border, border, border, |
311 | 0 | d_width - 2 * border, d_height - 2 * border, d_stride, |
312 | 0 | 0, 0, AOMMIN(best_error, error_adv_thresh), |
313 | 0 | segment_map, segment_map_stride); |
314 | 0 | if (step_error < best_error) { |
315 | 0 | best_error = step_error; |
316 | 0 | best_param = *param; |
317 | 0 | step_dir = -1; |
318 | 0 | } |
319 | | |
320 | | // look to the right |
321 | 0 | *param = add_param_offset(p, curr_param, step); |
322 | 0 | step_error = |
323 | 0 | av1_warp_error(wm, use_hbd, bd, ref, r_width, r_height, r_stride, |
324 | 0 | dst + border * d_stride + border, border, border, |
325 | 0 | d_width - 2 * border, d_height - 2 * border, d_stride, |
326 | 0 | 0, 0, AOMMIN(best_error, error_adv_thresh), |
327 | 0 | segment_map, segment_map_stride); |
328 | 0 | if (step_error < best_error) { |
329 | 0 | best_error = step_error; |
330 | 0 | best_param = *param; |
331 | 0 | step_dir = 1; |
332 | 0 | } |
333 | 0 | *param = best_param; |
334 | | |
335 | | // look to the direction chosen above repeatedly until error increases |
336 | | // for the biggest step size |
337 | 0 | while (step_dir) { |
338 | 0 | *param = add_param_offset(p, best_param, step * step_dir); |
339 | 0 | step_error = |
340 | 0 | av1_warp_error(wm, use_hbd, bd, ref, r_width, r_height, r_stride, |
341 | 0 | dst + border * d_stride + border, border, border, |
342 | 0 | d_width - 2 * border, d_height - 2 * border, |
343 | 0 | d_stride, 0, 0, AOMMIN(best_error, error_adv_thresh), |
344 | 0 | segment_map, segment_map_stride); |
345 | 0 | if (step_error < best_error) { |
346 | 0 | best_error = step_error; |
347 | 0 | best_param = *param; |
348 | 0 | } else { |
349 | 0 | *param = best_param; |
350 | 0 | step_dir = 0; |
351 | 0 | } |
352 | 0 | } |
353 | 0 | } |
354 | 0 | } |
355 | 0 | force_wmtype(wm, wmtype); |
356 | 0 | wm->wmtype = get_wmtype(wm); |
357 | 0 | return best_error; |
358 | 0 | } |
359 | | |
360 | 0 | unsigned char *av1_downconvert_frame(YV12_BUFFER_CONFIG *frm, int bit_depth) { |
361 | 0 | int i, j; |
362 | 0 | uint16_t *orig_buf = CONVERT_TO_SHORTPTR(frm->y_buffer); |
363 | 0 | uint8_t *buf_8bit = frm->y_buffer_8bit; |
364 | 0 | assert(buf_8bit); |
365 | 0 | if (!frm->buf_8bit_valid) { |
366 | 0 | for (i = 0; i < frm->y_height; ++i) { |
367 | 0 | for (j = 0; j < frm->y_width; ++j) { |
368 | 0 | buf_8bit[i * frm->y_stride + j] = |
369 | 0 | orig_buf[i * frm->y_stride + j] >> (bit_depth - 8); |
370 | 0 | } |
371 | 0 | } |
372 | 0 | frm->buf_8bit_valid = 1; |
373 | 0 | } |
374 | 0 | return buf_8bit; |
375 | 0 | } |
376 | | |
377 | | static void get_inliers_from_indices(MotionModel *params, |
378 | 0 | int *correspondences) { |
379 | 0 | int *inliers_tmp = (int *)aom_malloc(2 * MAX_CORNERS * sizeof(*inliers_tmp)); |
380 | 0 | memset(inliers_tmp, 0, 2 * MAX_CORNERS * sizeof(*inliers_tmp)); |
381 | |
|
382 | 0 | for (int i = 0; i < params->num_inliers; i++) { |
383 | 0 | int index = params->inliers[i]; |
384 | 0 | inliers_tmp[2 * i] = correspondences[4 * index]; |
385 | 0 | inliers_tmp[2 * i + 1] = correspondences[4 * index + 1]; |
386 | 0 | } |
387 | 0 | memcpy(params->inliers, inliers_tmp, sizeof(*inliers_tmp) * 2 * MAX_CORNERS); |
388 | 0 | aom_free(inliers_tmp); |
389 | 0 | } |
390 | | |
391 | 0 | #define FEAT_COUNT_TR 3 |
392 | 0 | #define SEG_COUNT_TR 0.40 |
393 | | void av1_compute_feature_segmentation_map(uint8_t *segment_map, int width, |
394 | | int height, int *inliers, |
395 | 0 | int num_inliers) { |
396 | 0 | int seg_count = 0; |
397 | 0 | memset(segment_map, 0, sizeof(*segment_map) * width * height); |
398 | |
|
399 | 0 | for (int i = 0; i < num_inliers; i++) { |
400 | 0 | int x = inliers[i * 2]; |
401 | 0 | int y = inliers[i * 2 + 1]; |
402 | 0 | int seg_x = x >> WARP_ERROR_BLOCK_LOG; |
403 | 0 | int seg_y = y >> WARP_ERROR_BLOCK_LOG; |
404 | 0 | segment_map[seg_y * width + seg_x] += 1; |
405 | 0 | } |
406 | |
|
407 | 0 | for (int i = 0; i < height; i++) { |
408 | 0 | for (int j = 0; j < width; j++) { |
409 | 0 | uint8_t feat_count = segment_map[i * width + j]; |
410 | 0 | segment_map[i * width + j] = (feat_count >= FEAT_COUNT_TR); |
411 | 0 | seg_count += (segment_map[i * width + j]); |
412 | 0 | } |
413 | 0 | } |
414 | | |
415 | | // If this motion does not make up a large enough portion of the frame, |
416 | | // use the unsegmented version of the error metric |
417 | 0 | if (seg_count < (width * height * SEG_COUNT_TR)) |
418 | 0 | memset(segment_map, 1, width * height * sizeof(*segment_map)); |
419 | 0 | } |
420 | | |
421 | | static int compute_global_motion_feature_based( |
422 | | TransformationType type, unsigned char *src_buffer, int src_width, |
423 | | int src_height, int src_stride, int *src_corners, int num_src_corners, |
424 | | YV12_BUFFER_CONFIG *ref, int bit_depth, int *num_inliers_by_motion, |
425 | 0 | MotionModel *params_by_motion, int num_motions) { |
426 | 0 | int i; |
427 | 0 | int num_ref_corners; |
428 | 0 | int num_correspondences; |
429 | 0 | int *correspondences; |
430 | 0 | int ref_corners[2 * MAX_CORNERS]; |
431 | 0 | unsigned char *ref_buffer = ref->y_buffer; |
432 | 0 | RansacFunc ransac = av1_get_ransac_type(type); |
433 | |
|
434 | 0 | if (ref->flags & YV12_FLAG_HIGHBITDEPTH) { |
435 | 0 | ref_buffer = av1_downconvert_frame(ref, bit_depth); |
436 | 0 | } |
437 | |
|
438 | 0 | num_ref_corners = |
439 | 0 | av1_fast_corner_detect(ref_buffer, ref->y_width, ref->y_height, |
440 | 0 | ref->y_stride, ref_corners, MAX_CORNERS); |
441 | | |
442 | | // find correspondences between the two images |
443 | 0 | correspondences = |
444 | 0 | (int *)malloc(num_src_corners * 4 * sizeof(*correspondences)); |
445 | 0 | num_correspondences = av1_determine_correspondence( |
446 | 0 | src_buffer, (int *)src_corners, num_src_corners, ref_buffer, |
447 | 0 | (int *)ref_corners, num_ref_corners, src_width, src_height, src_stride, |
448 | 0 | ref->y_stride, correspondences); |
449 | |
|
450 | 0 | ransac(correspondences, num_correspondences, num_inliers_by_motion, |
451 | 0 | params_by_motion, num_motions); |
452 | | |
453 | | // Set num_inliers = 0 for motions with too few inliers so they are ignored. |
454 | 0 | for (i = 0; i < num_motions; ++i) { |
455 | 0 | if (num_inliers_by_motion[i] < MIN_INLIER_PROB * num_correspondences || |
456 | 0 | num_correspondences == 0) { |
457 | 0 | num_inliers_by_motion[i] = 0; |
458 | 0 | } else { |
459 | 0 | get_inliers_from_indices(¶ms_by_motion[i], correspondences); |
460 | 0 | } |
461 | 0 | } |
462 | |
|
463 | 0 | free(correspondences); |
464 | | |
465 | | // Return true if any one of the motions has inliers. |
466 | 0 | for (i = 0; i < num_motions; ++i) { |
467 | 0 | if (num_inliers_by_motion[i] > 0) return 1; |
468 | 0 | } |
469 | 0 | return 0; |
470 | 0 | } |
471 | | |
472 | | // Don't use points around the frame border since they are less reliable |
473 | 0 | static INLINE int valid_point(int x, int y, int width, int height) { |
474 | 0 | return (x > (PATCH_SIZE + PATCH_CENTER)) && |
475 | 0 | (x < (width - PATCH_SIZE - PATCH_CENTER)) && |
476 | 0 | (y > (PATCH_SIZE + PATCH_CENTER)) && |
477 | 0 | (y < (height - PATCH_SIZE - PATCH_CENTER)); |
478 | 0 | } |
479 | | |
480 | | static int determine_disflow_correspondence(int *frm_corners, |
481 | | int num_frm_corners, double *flow_u, |
482 | | double *flow_v, int width, |
483 | | int height, int stride, |
484 | 0 | double *correspondences) { |
485 | 0 | int num_correspondences = 0; |
486 | 0 | int x, y; |
487 | 0 | for (int i = 0; i < num_frm_corners; ++i) { |
488 | 0 | x = frm_corners[2 * i]; |
489 | 0 | y = frm_corners[2 * i + 1]; |
490 | 0 | if (valid_point(x, y, width, height)) { |
491 | 0 | correspondences[4 * num_correspondences] = x; |
492 | 0 | correspondences[4 * num_correspondences + 1] = y; |
493 | 0 | correspondences[4 * num_correspondences + 2] = x + flow_u[y * stride + x]; |
494 | 0 | correspondences[4 * num_correspondences + 3] = y + flow_v[y * stride + x]; |
495 | 0 | num_correspondences++; |
496 | 0 | } |
497 | 0 | } |
498 | 0 | return num_correspondences; |
499 | 0 | } |
500 | | |
501 | 0 | static double getCubicValue(double p[4], double x) { |
502 | 0 | return p[1] + 0.5 * x * |
503 | 0 | (p[2] - p[0] + |
504 | 0 | x * (2.0 * p[0] - 5.0 * p[1] + 4.0 * p[2] - p[3] + |
505 | 0 | x * (3.0 * (p[1] - p[2]) + p[3] - p[0]))); |
506 | 0 | } |
507 | | |
508 | | static void get_subcolumn(unsigned char *ref, double col[4], int stride, int x, |
509 | 0 | int y_start) { |
510 | 0 | int i; |
511 | 0 | for (i = 0; i < 4; ++i) { |
512 | 0 | col[i] = ref[(i + y_start) * stride + x]; |
513 | 0 | } |
514 | 0 | } |
515 | | |
516 | 0 | static double bicubic(unsigned char *ref, double x, double y, int stride) { |
517 | 0 | double arr[4]; |
518 | 0 | int k; |
519 | 0 | int i = (int)x; |
520 | 0 | int j = (int)y; |
521 | 0 | for (k = 0; k < 4; ++k) { |
522 | 0 | double arr_temp[4]; |
523 | 0 | get_subcolumn(ref, arr_temp, stride, i + k - 1, j - 1); |
524 | 0 | arr[k] = getCubicValue(arr_temp, y - j); |
525 | 0 | } |
526 | 0 | return getCubicValue(arr, x - i); |
527 | 0 | } |
528 | | |
529 | | // Interpolate a warped block using bicubic interpolation when possible |
530 | | static unsigned char interpolate(unsigned char *ref, double x, double y, |
531 | 0 | int width, int height, int stride) { |
532 | 0 | if (x < 0 && y < 0) |
533 | 0 | return ref[0]; |
534 | 0 | else if (x < 0 && y > height - 1) |
535 | 0 | return ref[(height - 1) * stride]; |
536 | 0 | else if (x > width - 1 && y < 0) |
537 | 0 | return ref[width - 1]; |
538 | 0 | else if (x > width - 1 && y > height - 1) |
539 | 0 | return ref[(height - 1) * stride + (width - 1)]; |
540 | 0 | else if (x < 0) { |
541 | 0 | int v; |
542 | 0 | int i = (int)y; |
543 | 0 | double a = y - i; |
544 | 0 | if (y > 1 && y < height - 2) { |
545 | 0 | double arr[4]; |
546 | 0 | get_subcolumn(ref, arr, stride, 0, i - 1); |
547 | 0 | return clamp((int)(getCubicValue(arr, a) + 0.5), 0, 255); |
548 | 0 | } |
549 | 0 | v = (int)(ref[i * stride] * (1 - a) + ref[(i + 1) * stride] * a + 0.5); |
550 | 0 | return clamp(v, 0, 255); |
551 | 0 | } else if (y < 0) { |
552 | 0 | int v; |
553 | 0 | int j = (int)x; |
554 | 0 | double b = x - j; |
555 | 0 | if (x > 1 && x < width - 2) { |
556 | 0 | double arr[4] = { ref[j - 1], ref[j], ref[j + 1], ref[j + 2] }; |
557 | 0 | return clamp((int)(getCubicValue(arr, b) + 0.5), 0, 255); |
558 | 0 | } |
559 | 0 | v = (int)(ref[j] * (1 - b) + ref[j + 1] * b + 0.5); |
560 | 0 | return clamp(v, 0, 255); |
561 | 0 | } else if (x > width - 1) { |
562 | 0 | int v; |
563 | 0 | int i = (int)y; |
564 | 0 | double a = y - i; |
565 | 0 | if (y > 1 && y < height - 2) { |
566 | 0 | double arr[4]; |
567 | 0 | get_subcolumn(ref, arr, stride, width - 1, i - 1); |
568 | 0 | return clamp((int)(getCubicValue(arr, a) + 0.5), 0, 255); |
569 | 0 | } |
570 | 0 | v = (int)(ref[i * stride + width - 1] * (1 - a) + |
571 | 0 | ref[(i + 1) * stride + width - 1] * a + 0.5); |
572 | 0 | return clamp(v, 0, 255); |
573 | 0 | } else if (y > height - 1) { |
574 | 0 | int v; |
575 | 0 | int j = (int)x; |
576 | 0 | double b = x - j; |
577 | 0 | if (x > 1 && x < width - 2) { |
578 | 0 | int row = (height - 1) * stride; |
579 | 0 | double arr[4] = { ref[row + j - 1], ref[row + j], ref[row + j + 1], |
580 | 0 | ref[row + j + 2] }; |
581 | 0 | return clamp((int)(getCubicValue(arr, b) + 0.5), 0, 255); |
582 | 0 | } |
583 | 0 | v = (int)(ref[(height - 1) * stride + j] * (1 - b) + |
584 | 0 | ref[(height - 1) * stride + j + 1] * b + 0.5); |
585 | 0 | return clamp(v, 0, 255); |
586 | 0 | } else if (x > 1 && y > 1 && x < width - 2 && y < height - 2) { |
587 | 0 | return clamp((int)(bicubic(ref, x, y, stride) + 0.5), 0, 255); |
588 | 0 | } else { |
589 | 0 | int i = (int)y; |
590 | 0 | int j = (int)x; |
591 | 0 | double a = y - i; |
592 | 0 | double b = x - j; |
593 | 0 | int v = (int)(ref[i * stride + j] * (1 - a) * (1 - b) + |
594 | 0 | ref[i * stride + j + 1] * (1 - a) * b + |
595 | 0 | ref[(i + 1) * stride + j] * a * (1 - b) + |
596 | 0 | ref[(i + 1) * stride + j + 1] * a * b); |
597 | 0 | return clamp(v, 0, 255); |
598 | 0 | } |
599 | 0 | } |
600 | | |
601 | | // Warps a block using flow vector [u, v] and computes the mse |
602 | | static double compute_warp_and_error(unsigned char *ref, unsigned char *frm, |
603 | | int width, int height, int stride, int x, |
604 | 0 | int y, double u, double v, int16_t *dt) { |
605 | 0 | int i, j; |
606 | 0 | unsigned char warped; |
607 | 0 | double x_w, y_w; |
608 | 0 | double mse = 0; |
609 | 0 | int16_t err = 0; |
610 | 0 | for (i = y; i < y + PATCH_SIZE; ++i) |
611 | 0 | for (j = x; j < x + PATCH_SIZE; ++j) { |
612 | 0 | x_w = (double)j + u; |
613 | 0 | y_w = (double)i + v; |
614 | 0 | warped = interpolate(ref, x_w, y_w, width, height, stride); |
615 | 0 | err = warped - frm[j + i * stride]; |
616 | 0 | mse += err * err; |
617 | 0 | dt[(i - y) * PATCH_SIZE + (j - x)] = err; |
618 | 0 | } |
619 | |
|
620 | 0 | mse /= (PATCH_SIZE * PATCH_SIZE); |
621 | 0 | return mse; |
622 | 0 | } |
623 | | |
624 | | // Computes the components of the system of equations used to solve for |
625 | | // a flow vector. This includes: |
626 | | // 1.) The hessian matrix for optical flow. This matrix is in the |
627 | | // form of: |
628 | | // |
629 | | // M = |sum(dx * dx) sum(dx * dy)| |
630 | | // |sum(dx * dy) sum(dy * dy)| |
631 | | // |
632 | | // 2.) b = |sum(dx * dt)| |
633 | | // |sum(dy * dt)| |
634 | | // Where the sums are computed over a square window of PATCH_SIZE. |
635 | | static INLINE void compute_flow_system(const double *dx, int dx_stride, |
636 | | const double *dy, int dy_stride, |
637 | | const int16_t *dt, int dt_stride, |
638 | 0 | double *M, double *b) { |
639 | 0 | for (int i = 0; i < PATCH_SIZE; i++) { |
640 | 0 | for (int j = 0; j < PATCH_SIZE; j++) { |
641 | 0 | M[0] += dx[i * dx_stride + j] * dx[i * dx_stride + j]; |
642 | 0 | M[1] += dx[i * dx_stride + j] * dy[i * dy_stride + j]; |
643 | 0 | M[3] += dy[i * dy_stride + j] * dy[i * dy_stride + j]; |
644 | |
|
645 | 0 | b[0] += dx[i * dx_stride + j] * dt[i * dt_stride + j]; |
646 | 0 | b[1] += dy[i * dy_stride + j] * dt[i * dt_stride + j]; |
647 | 0 | } |
648 | 0 | } |
649 | |
|
650 | 0 | M[2] = M[1]; |
651 | 0 | } |
652 | | |
653 | | // Solves a general Mx = b where M is a 2x2 matrix and b is a 2x1 matrix |
654 | | static INLINE void solve_2x2_system(const double *M, const double *b, |
655 | 0 | double *output_vec) { |
656 | 0 | double M_0 = M[0]; |
657 | 0 | double M_3 = M[3]; |
658 | 0 | double det = (M_0 * M_3) - (M[1] * M[2]); |
659 | 0 | if (det < 1e-5) { |
660 | | // Handle singular matrix |
661 | | // TODO(sarahparker) compare results using pseudo inverse instead |
662 | 0 | M_0 += 1e-10; |
663 | 0 | M_3 += 1e-10; |
664 | 0 | det = (M_0 * M_3) - (M[1] * M[2]); |
665 | 0 | } |
666 | 0 | const double det_inv = 1 / det; |
667 | 0 | const double mult_b0 = det_inv * b[0]; |
668 | 0 | const double mult_b1 = det_inv * b[1]; |
669 | 0 | output_vec[0] = M_3 * mult_b0 - M[1] * mult_b1; |
670 | 0 | output_vec[1] = -M[2] * mult_b0 + M_0 * mult_b1; |
671 | 0 | } |
672 | | |
673 | | /* |
674 | | static INLINE void image_difference(const uint8_t *src, int src_stride, |
675 | | const uint8_t *ref, int ref_stride, |
676 | | int16_t *dst, int dst_stride, int height, |
677 | | int width) { |
678 | | const int block_unit = 8; |
679 | | // Take difference in 8x8 blocks to make use of optimized diff function |
680 | | for (int i = 0; i < height; i += block_unit) { |
681 | | for (int j = 0; j < width; j += block_unit) { |
682 | | aom_subtract_block(block_unit, block_unit, dst + i * dst_stride + j, |
683 | | dst_stride, src + i * src_stride + j, src_stride, |
684 | | ref + i * ref_stride + j, ref_stride); |
685 | | } |
686 | | } |
687 | | } |
688 | | */ |
689 | | |
690 | | // Compute an image gradient using a sobel filter. |
691 | | // If dir == 1, compute the x gradient. If dir == 0, compute y. This function |
692 | | // assumes the images have been padded so that they can be processed in units |
693 | | // of 8. |
694 | | static INLINE void sobel_xy_image_gradient(const uint8_t *src, int src_stride, |
695 | | double *dst, int dst_stride, |
696 | 0 | int height, int width, int dir) { |
697 | 0 | double norm = 1.0; |
698 | | // TODO(sarahparker) experiment with doing this over larger block sizes |
699 | 0 | const int block_unit = 8; |
700 | | // Filter in 8x8 blocks to eventually make use of optimized convolve function |
701 | 0 | for (int i = 0; i < height; i += block_unit) { |
702 | 0 | for (int j = 0; j < width; j += block_unit) { |
703 | 0 | av1_convolve_2d_sobel_y_c(src + i * src_stride + j, src_stride, |
704 | 0 | dst + i * dst_stride + j, dst_stride, |
705 | 0 | block_unit, block_unit, dir, norm); |
706 | 0 | } |
707 | 0 | } |
708 | 0 | } |
709 | | |
710 | | static ImagePyramid *alloc_pyramid(int width, int height, int pad_size, |
711 | 0 | int compute_gradient) { |
712 | 0 | ImagePyramid *pyr = aom_malloc(sizeof(*pyr)); |
713 | 0 | pyr->has_gradient = compute_gradient; |
714 | | // 2 * width * height is the upper bound for a buffer that fits |
715 | | // all pyramid levels + padding for each level |
716 | 0 | const int buffer_size = sizeof(*pyr->level_buffer) * 2 * width * height + |
717 | 0 | (width + 2 * pad_size) * 2 * pad_size * N_LEVELS; |
718 | 0 | pyr->level_buffer = aom_malloc(buffer_size); |
719 | 0 | memset(pyr->level_buffer, 0, buffer_size); |
720 | |
|
721 | 0 | if (compute_gradient) { |
722 | 0 | const int gradient_size = |
723 | 0 | sizeof(*pyr->level_dx_buffer) * 2 * width * height + |
724 | 0 | (width + 2 * pad_size) * 2 * pad_size * N_LEVELS; |
725 | 0 | pyr->level_dx_buffer = aom_malloc(gradient_size); |
726 | 0 | pyr->level_dy_buffer = aom_malloc(gradient_size); |
727 | 0 | memset(pyr->level_dx_buffer, 0, gradient_size); |
728 | 0 | memset(pyr->level_dy_buffer, 0, gradient_size); |
729 | 0 | } |
730 | 0 | return pyr; |
731 | 0 | } |
732 | | |
733 | 0 | static void free_pyramid(ImagePyramid *pyr) { |
734 | 0 | aom_free(pyr->level_buffer); |
735 | 0 | if (pyr->has_gradient) { |
736 | 0 | aom_free(pyr->level_dx_buffer); |
737 | 0 | aom_free(pyr->level_dy_buffer); |
738 | 0 | } |
739 | 0 | aom_free(pyr); |
740 | 0 | } |
741 | | |
742 | 0 | static INLINE void update_level_dims(ImagePyramid *frm_pyr, int level) { |
743 | 0 | frm_pyr->widths[level] = frm_pyr->widths[level - 1] >> 1; |
744 | 0 | frm_pyr->heights[level] = frm_pyr->heights[level - 1] >> 1; |
745 | 0 | frm_pyr->strides[level] = frm_pyr->widths[level] + 2 * frm_pyr->pad_size; |
746 | | // Point the beginning of the next level buffer to the correct location inside |
747 | | // the padded border |
748 | 0 | frm_pyr->level_loc[level] = |
749 | 0 | frm_pyr->level_loc[level - 1] + |
750 | 0 | frm_pyr->strides[level - 1] * |
751 | 0 | (2 * frm_pyr->pad_size + frm_pyr->heights[level - 1]); |
752 | 0 | } |
753 | | |
754 | | // Compute coarse to fine pyramids for a frame |
755 | | static void compute_flow_pyramids(unsigned char *frm, const int frm_width, |
756 | | const int frm_height, const int frm_stride, |
757 | | int n_levels, int pad_size, int compute_grad, |
758 | 0 | ImagePyramid *frm_pyr) { |
759 | 0 | int cur_width, cur_height, cur_stride, cur_loc; |
760 | 0 | assert((frm_width >> n_levels) > 0); |
761 | 0 | assert((frm_height >> n_levels) > 0); |
762 | | |
763 | | // Initialize first level |
764 | 0 | frm_pyr->n_levels = n_levels; |
765 | 0 | frm_pyr->pad_size = pad_size; |
766 | 0 | frm_pyr->widths[0] = frm_width; |
767 | 0 | frm_pyr->heights[0] = frm_height; |
768 | 0 | frm_pyr->strides[0] = frm_width + 2 * frm_pyr->pad_size; |
769 | | // Point the beginning of the level buffer to the location inside |
770 | | // the padded border |
771 | 0 | frm_pyr->level_loc[0] = |
772 | 0 | frm_pyr->strides[0] * frm_pyr->pad_size + frm_pyr->pad_size; |
773 | | // This essentially copies the original buffer into the pyramid buffer |
774 | | // without the original padding |
775 | 0 | av1_resize_plane(frm, frm_height, frm_width, frm_stride, |
776 | 0 | frm_pyr->level_buffer + frm_pyr->level_loc[0], |
777 | 0 | frm_pyr->heights[0], frm_pyr->widths[0], |
778 | 0 | frm_pyr->strides[0]); |
779 | |
|
780 | 0 | if (compute_grad) { |
781 | 0 | cur_width = frm_pyr->widths[0]; |
782 | 0 | cur_height = frm_pyr->heights[0]; |
783 | 0 | cur_stride = frm_pyr->strides[0]; |
784 | 0 | cur_loc = frm_pyr->level_loc[0]; |
785 | 0 | assert(frm_pyr->has_gradient && frm_pyr->level_dx_buffer != NULL && |
786 | 0 | frm_pyr->level_dy_buffer != NULL); |
787 | | // Computation x gradient |
788 | 0 | sobel_xy_image_gradient(frm_pyr->level_buffer + cur_loc, cur_stride, |
789 | 0 | frm_pyr->level_dx_buffer + cur_loc, cur_stride, |
790 | 0 | cur_height, cur_width, 1); |
791 | | |
792 | | // Computation y gradient |
793 | 0 | sobel_xy_image_gradient(frm_pyr->level_buffer + cur_loc, cur_stride, |
794 | 0 | frm_pyr->level_dy_buffer + cur_loc, cur_stride, |
795 | 0 | cur_height, cur_width, 0); |
796 | 0 | } |
797 | | |
798 | | // Start at the finest level and resize down to the coarsest level |
799 | 0 | for (int level = 1; level < n_levels; ++level) { |
800 | 0 | update_level_dims(frm_pyr, level); |
801 | 0 | cur_width = frm_pyr->widths[level]; |
802 | 0 | cur_height = frm_pyr->heights[level]; |
803 | 0 | cur_stride = frm_pyr->strides[level]; |
804 | 0 | cur_loc = frm_pyr->level_loc[level]; |
805 | |
|
806 | 0 | av1_resize_plane(frm_pyr->level_buffer + frm_pyr->level_loc[level - 1], |
807 | 0 | frm_pyr->heights[level - 1], frm_pyr->widths[level - 1], |
808 | 0 | frm_pyr->strides[level - 1], |
809 | 0 | frm_pyr->level_buffer + cur_loc, cur_height, cur_width, |
810 | 0 | cur_stride); |
811 | |
|
812 | 0 | if (compute_grad) { |
813 | 0 | assert(frm_pyr->has_gradient && frm_pyr->level_dx_buffer != NULL && |
814 | 0 | frm_pyr->level_dy_buffer != NULL); |
815 | | // Computation x gradient |
816 | 0 | sobel_xy_image_gradient(frm_pyr->level_buffer + cur_loc, cur_stride, |
817 | 0 | frm_pyr->level_dx_buffer + cur_loc, cur_stride, |
818 | 0 | cur_height, cur_width, 1); |
819 | | |
820 | | // Computation y gradient |
821 | 0 | sobel_xy_image_gradient(frm_pyr->level_buffer + cur_loc, cur_stride, |
822 | 0 | frm_pyr->level_dy_buffer + cur_loc, cur_stride, |
823 | 0 | cur_height, cur_width, 0); |
824 | 0 | } |
825 | 0 | } |
826 | 0 | } |
827 | | |
828 | | static INLINE void compute_flow_at_point(unsigned char *frm, unsigned char *ref, |
829 | | double *dx, double *dy, int x, int y, |
830 | | int width, int height, int stride, |
831 | 0 | double *u, double *v) { |
832 | 0 | double M[4] = { 0 }; |
833 | 0 | double b[2] = { 0 }; |
834 | 0 | double tmp_output_vec[2] = { 0 }; |
835 | 0 | double error = 0; |
836 | 0 | int16_t dt[PATCH_SIZE * PATCH_SIZE]; |
837 | 0 | double o_u = *u; |
838 | 0 | double o_v = *v; |
839 | |
|
840 | 0 | for (int itr = 0; itr < DISFLOW_MAX_ITR; itr++) { |
841 | 0 | error = compute_warp_and_error(ref, frm, width, height, stride, x, y, *u, |
842 | 0 | *v, dt); |
843 | 0 | if (error <= DISFLOW_ERROR_TR) break; |
844 | 0 | compute_flow_system(dx, stride, dy, stride, dt, PATCH_SIZE, M, b); |
845 | 0 | solve_2x2_system(M, b, tmp_output_vec); |
846 | 0 | *u += tmp_output_vec[0]; |
847 | 0 | *v += tmp_output_vec[1]; |
848 | 0 | } |
849 | 0 | if (fabs(*u - o_u) > PATCH_SIZE || fabs(*v - o_u) > PATCH_SIZE) { |
850 | 0 | *u = o_u; |
851 | 0 | *v = o_v; |
852 | 0 | } |
853 | 0 | } |
854 | | |
855 | | // make sure flow_u and flow_v start at 0 |
856 | | static void compute_flow_field(ImagePyramid *frm_pyr, ImagePyramid *ref_pyr, |
857 | 0 | double *flow_u, double *flow_v) { |
858 | 0 | int cur_width, cur_height, cur_stride, cur_loc, patch_loc, patch_center; |
859 | 0 | double *u_upscale = |
860 | 0 | aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_u)); |
861 | 0 | double *v_upscale = |
862 | 0 | aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_v)); |
863 | |
|
864 | 0 | assert(frm_pyr->n_levels == ref_pyr->n_levels); |
865 | | |
866 | | // Compute flow field from coarsest to finest level of the pyramid |
867 | 0 | for (int level = frm_pyr->n_levels - 1; level >= 0; --level) { |
868 | 0 | cur_width = frm_pyr->widths[level]; |
869 | 0 | cur_height = frm_pyr->heights[level]; |
870 | 0 | cur_stride = frm_pyr->strides[level]; |
871 | 0 | cur_loc = frm_pyr->level_loc[level]; |
872 | |
|
873 | 0 | for (int i = PATCH_SIZE; i < cur_height - PATCH_SIZE; i += PATCH_STEP) { |
874 | 0 | for (int j = PATCH_SIZE; j < cur_width - PATCH_SIZE; j += PATCH_STEP) { |
875 | 0 | patch_loc = i * cur_stride + j; |
876 | 0 | patch_center = patch_loc + PATCH_CENTER * cur_stride + PATCH_CENTER; |
877 | 0 | compute_flow_at_point(frm_pyr->level_buffer + cur_loc, |
878 | 0 | ref_pyr->level_buffer + cur_loc, |
879 | 0 | frm_pyr->level_dx_buffer + cur_loc + patch_loc, |
880 | 0 | frm_pyr->level_dy_buffer + cur_loc + patch_loc, j, |
881 | 0 | i, cur_width, cur_height, cur_stride, |
882 | 0 | flow_u + patch_center, flow_v + patch_center); |
883 | 0 | } |
884 | 0 | } |
885 | | // TODO(sarahparker) Replace this with upscale function in resize.c |
886 | 0 | if (level > 0) { |
887 | 0 | int h_upscale = frm_pyr->heights[level - 1]; |
888 | 0 | int w_upscale = frm_pyr->widths[level - 1]; |
889 | 0 | int s_upscale = frm_pyr->strides[level - 1]; |
890 | 0 | for (int i = 0; i < h_upscale; ++i) { |
891 | 0 | for (int j = 0; j < w_upscale; ++j) { |
892 | 0 | u_upscale[j + i * s_upscale] = |
893 | 0 | flow_u[(int)(j >> 1) + (int)(i >> 1) * cur_stride]; |
894 | 0 | v_upscale[j + i * s_upscale] = |
895 | 0 | flow_v[(int)(j >> 1) + (int)(i >> 1) * cur_stride]; |
896 | 0 | } |
897 | 0 | } |
898 | 0 | memcpy(flow_u, u_upscale, |
899 | 0 | frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_u)); |
900 | 0 | memcpy(flow_v, v_upscale, |
901 | 0 | frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_v)); |
902 | 0 | } |
903 | 0 | } |
904 | 0 | aom_free(u_upscale); |
905 | 0 | aom_free(v_upscale); |
906 | 0 | } |
907 | | |
908 | | static int compute_global_motion_disflow_based( |
909 | | TransformationType type, unsigned char *frm_buffer, int frm_width, |
910 | | int frm_height, int frm_stride, int *frm_corners, int num_frm_corners, |
911 | | YV12_BUFFER_CONFIG *ref, int bit_depth, int *num_inliers_by_motion, |
912 | 0 | MotionModel *params_by_motion, int num_motions) { |
913 | 0 | unsigned char *ref_buffer = ref->y_buffer; |
914 | 0 | const int ref_width = ref->y_width; |
915 | 0 | const int ref_height = ref->y_height; |
916 | 0 | const int pad_size = AOMMAX(PATCH_SIZE, MIN_PAD); |
917 | 0 | int num_correspondences; |
918 | 0 | double *correspondences; |
919 | 0 | RansacFuncDouble ransac = av1_get_ransac_double_prec_type(type); |
920 | 0 | assert(frm_width == ref_width); |
921 | 0 | assert(frm_height == ref_height); |
922 | | |
923 | | // Ensure the number of pyramid levels will work with the frame resolution |
924 | 0 | const int msb = |
925 | 0 | frm_width < frm_height ? get_msb(frm_width) : get_msb(frm_height); |
926 | 0 | const int n_levels = AOMMIN(msb, N_LEVELS); |
927 | |
|
928 | 0 | if (ref->flags & YV12_FLAG_HIGHBITDEPTH) { |
929 | 0 | ref_buffer = av1_downconvert_frame(ref, bit_depth); |
930 | 0 | } |
931 | | |
932 | | // TODO(sarahparker) We will want to do the source pyramid computation |
933 | | // outside of this function so it doesn't get recomputed for every |
934 | | // reference. We also don't need to compute every pyramid level for the |
935 | | // reference in advance, since lower levels can be overwritten once their |
936 | | // flow field is computed and upscaled. I'll add these optimizations |
937 | | // once the full implementation is working. |
938 | | // Allocate frm image pyramids |
939 | 0 | int compute_gradient = 1; |
940 | 0 | ImagePyramid *frm_pyr = |
941 | 0 | alloc_pyramid(frm_width, frm_height, pad_size, compute_gradient); |
942 | 0 | compute_flow_pyramids(frm_buffer, frm_width, frm_height, frm_stride, n_levels, |
943 | 0 | pad_size, compute_gradient, frm_pyr); |
944 | | // Allocate ref image pyramids |
945 | 0 | compute_gradient = 0; |
946 | 0 | ImagePyramid *ref_pyr = |
947 | 0 | alloc_pyramid(ref_width, ref_height, pad_size, compute_gradient); |
948 | 0 | compute_flow_pyramids(ref_buffer, ref_width, ref_height, ref->y_stride, |
949 | 0 | n_levels, pad_size, compute_gradient, ref_pyr); |
950 | |
|
951 | 0 | double *flow_u = |
952 | 0 | aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_u)); |
953 | 0 | double *flow_v = |
954 | 0 | aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_v)); |
955 | |
|
956 | 0 | memset(flow_u, 0, |
957 | 0 | frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_u)); |
958 | 0 | memset(flow_v, 0, |
959 | 0 | frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_v)); |
960 | |
|
961 | 0 | compute_flow_field(frm_pyr, ref_pyr, flow_u, flow_v); |
962 | | |
963 | | // find correspondences between the two images using the flow field |
964 | 0 | correspondences = aom_malloc(num_frm_corners * 4 * sizeof(*correspondences)); |
965 | 0 | num_correspondences = determine_disflow_correspondence( |
966 | 0 | frm_corners, num_frm_corners, flow_u, flow_v, frm_width, frm_height, |
967 | 0 | frm_pyr->strides[0], correspondences); |
968 | 0 | ransac(correspondences, num_correspondences, num_inliers_by_motion, |
969 | 0 | params_by_motion, num_motions); |
970 | |
|
971 | 0 | free_pyramid(frm_pyr); |
972 | 0 | free_pyramid(ref_pyr); |
973 | 0 | aom_free(correspondences); |
974 | 0 | aom_free(flow_u); |
975 | 0 | aom_free(flow_v); |
976 | | // Set num_inliers = 0 for motions with too few inliers so they are ignored. |
977 | 0 | for (int i = 0; i < num_motions; ++i) { |
978 | 0 | if (num_inliers_by_motion[i] < MIN_INLIER_PROB * num_correspondences) { |
979 | 0 | num_inliers_by_motion[i] = 0; |
980 | 0 | } |
981 | 0 | } |
982 | | |
983 | | // Return true if any one of the motions has inliers. |
984 | 0 | for (int i = 0; i < num_motions; ++i) { |
985 | 0 | if (num_inliers_by_motion[i] > 0) return 1; |
986 | 0 | } |
987 | 0 | return 0; |
988 | 0 | } |
989 | | |
990 | | int av1_compute_global_motion(TransformationType type, |
991 | | unsigned char *src_buffer, int src_width, |
992 | | int src_height, int src_stride, int *src_corners, |
993 | | int num_src_corners, YV12_BUFFER_CONFIG *ref, |
994 | | int bit_depth, |
995 | | GlobalMotionEstimationType gm_estimation_type, |
996 | | int *num_inliers_by_motion, |
997 | 0 | MotionModel *params_by_motion, int num_motions) { |
998 | 0 | switch (gm_estimation_type) { |
999 | 0 | case GLOBAL_MOTION_FEATURE_BASED: |
1000 | 0 | return compute_global_motion_feature_based( |
1001 | 0 | type, src_buffer, src_width, src_height, src_stride, src_corners, |
1002 | 0 | num_src_corners, ref, bit_depth, num_inliers_by_motion, |
1003 | 0 | params_by_motion, num_motions); |
1004 | 0 | case GLOBAL_MOTION_DISFLOW_BASED: |
1005 | 0 | return compute_global_motion_disflow_based( |
1006 | 0 | type, src_buffer, src_width, src_height, src_stride, src_corners, |
1007 | 0 | num_src_corners, ref, bit_depth, num_inliers_by_motion, |
1008 | 0 | params_by_motion, num_motions); |
1009 | 0 | default: assert(0 && "Unknown global motion estimation type"); |
1010 | 0 | } |
1011 | 0 | return 0; |
1012 | 0 | } |