| /* |
| * Copyright (c) 2016, Alliance for Open Media. All rights reserved. |
| * |
| * This source code is subject to the terms of the BSD 2 Clause License and |
| * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License |
| * was not distributed with this source code in the LICENSE file, you can |
| * obtain it at www.aomedia.org/license/software. If the Alliance for Open |
| * Media Patent License 1.0 was not distributed with this source code in the |
| * PATENTS file, you can obtain it at www.aomedia.org/license/patent. |
| */ |
| #include <math.h> |
| #include <limits.h> |
| |
| #include "config/aom_config.h" |
| |
| #include "aom_dsp/mathutils.h" |
| #include "aom_mem/aom_mem.h" |
| |
| #include "av1/common/av1_common_int.h" |
| #include "av1/encoder/encoder.h" |
| #include "av1/encoder/optical_flow.h" |
| #include "av1/encoder/sparse_linear_solver.h" |
| #include "av1/encoder/reconinter_enc.h" |
| |
| #if CONFIG_OPTICAL_FLOW_API |
| |
| void av1_init_opfl_params(OPFL_PARAMS *opfl_params) { |
| opfl_params->pyramid_levels = OPFL_PYRAMID_LEVELS; |
| opfl_params->warping_steps = OPFL_WARPING_STEPS; |
| opfl_params->lk_params = NULL; |
| } |
| |
| void av1_init_lk_params(LK_PARAMS *lk_params) { |
| lk_params->window_size = OPFL_WINDOW_SIZE; |
| } |
| |
| // Helper function to determine whether a frame is encoded with high bit-depth. |
| static inline int is_frame_high_bitdepth(const YV12_BUFFER_CONFIG *frame) { |
| return (frame->flags & YV12_FLAG_HIGHBITDEPTH) ? 1 : 0; |
| } |
| |
| // Helper function to determine whether optical flow method is sparse. |
| static inline int is_sparse(const OPFL_PARAMS *opfl_params) { |
| return (opfl_params->flags & OPFL_FLAG_SPARSE) ? 1 : 0; |
| } |
| |
| static void gradients_over_window(const YV12_BUFFER_CONFIG *frame, |
| const YV12_BUFFER_CONFIG *ref_frame, |
| const double x_coord, const double y_coord, |
| const int window_size, const int bit_depth, |
| double *ix, double *iy, double *it, |
| LOCALMV *mv); |
| |
| // coefficients for bilinear interpolation on unit square |
| static int pixel_interp(const double x, const double y, const double b00, |
| const double b01, const double b10, const double b11) { |
| const int xint = (int)x; |
| const int yint = (int)y; |
| const double xdec = x - xint; |
| const double ydec = y - yint; |
| const double a = (1 - xdec) * (1 - ydec); |
| const double b = xdec * (1 - ydec); |
| const double c = (1 - xdec) * ydec; |
| const double d = xdec * ydec; |
| // if x, y are already integers, this results to b00 |
| int interp = (int)round(a * b00 + b * b01 + c * b10 + d * b11); |
| return interp; |
| } |
| |
| // Scharr filter to compute spatial gradient |
| static void spatial_gradient(const YV12_BUFFER_CONFIG *frame, const int x_coord, |
| const int y_coord, const int direction, |
| double *derivative) { |
| double *filter; |
| // Scharr filters |
| double gx[9] = { -3, 0, 3, -10, 0, 10, -3, 0, 3 }; |
| double gy[9] = { -3, -10, -3, 0, 0, 0, 3, 10, 3 }; |
| if (direction == 0) { // x direction |
| filter = gx; |
| } else { // y direction |
| filter = gy; |
| } |
| int idx = 0; |
| double d = 0; |
| for (int yy = -1; yy <= 1; yy++) { |
| for (int xx = -1; xx <= 1; xx++) { |
| d += filter[idx] * |
| frame->y_buffer[(y_coord + yy) * frame->y_stride + (x_coord + xx)]; |
| idx++; |
| } |
| } |
| // normalization scaling factor for scharr |
| *derivative = d / 32.0; |
| } |
| |
| // Determine the spatial gradient at subpixel locations |
| // For example, when reducing images for pyramidal LK, |
| // corners found in original image may be at subpixel locations. |
| static void gradient_interp(double *fullpel_deriv, const double x_coord, |
| const double y_coord, const int w, const int h, |
| double *derivative) { |
| const int xint = (int)x_coord; |
| const int yint = (int)y_coord; |
| double interp; |
| if (xint + 1 > w - 1 || yint + 1 > h - 1) { |
| interp = fullpel_deriv[yint * w + xint]; |
| } else { |
| interp = pixel_interp(x_coord, y_coord, fullpel_deriv[yint * w + xint], |
| fullpel_deriv[yint * w + (xint + 1)], |
| fullpel_deriv[(yint + 1) * w + xint], |
| fullpel_deriv[(yint + 1) * w + (xint + 1)]); |
| } |
| |
| *derivative = interp; |
| } |
| |
| static void temporal_gradient(const YV12_BUFFER_CONFIG *frame, |
| const YV12_BUFFER_CONFIG *frame2, |
| const double x_coord, const double y_coord, |
| const int bit_depth, double *derivative, |
| LOCALMV *mv) { |
| const int w = 2; |
| const int h = 2; |
| uint8_t pred1[4]; |
| uint8_t pred2[4]; |
| |
| const int y = (int)y_coord; |
| const int x = (int)x_coord; |
| const double ydec = y_coord - y; |
| const double xdec = x_coord - x; |
| const int is_intrabc = 0; // Is intra-copied? |
| const int is_high_bitdepth = is_frame_high_bitdepth(frame2); |
| const int subsampling_x = 0, subsampling_y = 0; // for y-buffer |
| const int_interpfilters interp_filters = |
| av1_broadcast_interp_filter(MULTITAP_SHARP); |
| const int plane = 0; // y-plane |
| const struct buf_2d ref_buf2 = { NULL, frame2->y_buffer, frame2->y_crop_width, |
| frame2->y_crop_height, frame2->y_stride }; |
| struct scale_factors scale; |
| av1_setup_scale_factors_for_frame(&scale, frame->y_crop_width, |
| frame->y_crop_height, frame->y_crop_width, |
| frame->y_crop_height); |
| InterPredParams inter_pred_params; |
| av1_init_inter_params(&inter_pred_params, w, h, y, x, subsampling_x, |
| subsampling_y, bit_depth, is_high_bitdepth, is_intrabc, |
| &scale, &ref_buf2, interp_filters); |
| inter_pred_params.interp_filter_params[0] = |
| &av1_interp_filter_params_list[interp_filters.as_filters.x_filter]; |
| inter_pred_params.interp_filter_params[1] = |
| &av1_interp_filter_params_list[interp_filters.as_filters.y_filter]; |
| inter_pred_params.conv_params = get_conv_params(0, plane, bit_depth); |
| MV newmv = { .row = (int16_t)round((mv->row + xdec) * 8), |
| .col = (int16_t)round((mv->col + ydec) * 8) }; |
| av1_enc_build_one_inter_predictor(pred2, w, &newmv, &inter_pred_params); |
| const struct buf_2d ref_buf1 = { NULL, frame->y_buffer, frame->y_crop_width, |
| frame->y_crop_height, frame->y_stride }; |
| av1_init_inter_params(&inter_pred_params, w, h, y, x, subsampling_x, |
| subsampling_y, bit_depth, is_high_bitdepth, is_intrabc, |
| &scale, &ref_buf1, interp_filters); |
| inter_pred_params.interp_filter_params[0] = |
| &av1_interp_filter_params_list[interp_filters.as_filters.x_filter]; |
| inter_pred_params.interp_filter_params[1] = |
| &av1_interp_filter_params_list[interp_filters.as_filters.y_filter]; |
| inter_pred_params.conv_params = get_conv_params(0, plane, bit_depth); |
| MV zeroMV = { .row = (int16_t)round(xdec * 8), |
| .col = (int16_t)round(ydec * 8) }; |
| av1_enc_build_one_inter_predictor(pred1, w, &zeroMV, &inter_pred_params); |
| |
| *derivative = pred2[0] - pred1[0]; |
| } |
| |
| // Numerical differentiate over window_size x window_size surrounding (x,y) |
| // location. Alters ix, iy, it to contain numerical partial derivatives |
| static void gradients_over_window(const YV12_BUFFER_CONFIG *frame, |
| const YV12_BUFFER_CONFIG *ref_frame, |
| const double x_coord, const double y_coord, |
| const int window_size, const int bit_depth, |
| double *ix, double *iy, double *it, |
| LOCALMV *mv) { |
| const double left = x_coord - window_size / 2.0; |
| const double top = y_coord - window_size / 2.0; |
| // gradient operators need pixel before and after (start at 1) |
| const double x_start = AOMMAX(1, left); |
| const double y_start = AOMMAX(1, top); |
| const int frame_height = frame->y_crop_height; |
| const int frame_width = frame->y_crop_width; |
| double deriv_x; |
| double deriv_y; |
| double deriv_t; |
| |
| const double x_end = AOMMIN(x_coord + window_size / 2.0, frame_width - 2); |
| const double y_end = AOMMIN(y_coord + window_size / 2.0, frame_height - 2); |
| const int xs = (int)AOMMAX(1, x_start - 1); |
| const int ys = (int)AOMMAX(1, y_start - 1); |
| const int xe = (int)AOMMIN(x_end + 2, frame_width - 2); |
| const int ye = (int)AOMMIN(y_end + 2, frame_height - 2); |
| // with normalization, gradients may be double values |
| double *fullpel_dx = aom_malloc((ye - ys) * (xe - xs) * sizeof(deriv_x)); |
| double *fullpel_dy = aom_malloc((ye - ys) * (xe - xs) * sizeof(deriv_y)); |
| if (!fullpel_dx || !fullpel_dy) { |
| aom_free(fullpel_dx); |
| aom_free(fullpel_dy); |
| return; |
| } |
| |
| // TODO(any): This could be more efficient in the case that x_coord |
| // and y_coord are integers.. but it may look more messy. |
| |
| // calculate spatial gradients at full pixel locations |
| for (int j = ys; j < ye; j++) { |
| for (int i = xs; i < xe; i++) { |
| spatial_gradient(frame, i, j, 0, &deriv_x); |
| spatial_gradient(frame, i, j, 1, &deriv_y); |
| int idx = (j - ys) * (xe - xs) + (i - xs); |
| fullpel_dx[idx] = deriv_x; |
| fullpel_dy[idx] = deriv_y; |
| } |
| } |
| // compute numerical differentiation for every pixel in window |
| // (this potentially includes subpixels) |
| for (double j = y_start; j < y_end; j++) { |
| for (double i = x_start; i < x_end; i++) { |
| temporal_gradient(frame, ref_frame, i, j, bit_depth, &deriv_t, mv); |
| gradient_interp(fullpel_dx, i - xs, j - ys, xe - xs, ye - ys, &deriv_x); |
| gradient_interp(fullpel_dy, i - xs, j - ys, xe - xs, ye - ys, &deriv_y); |
| int idx = (int)(j - top) * window_size + (int)(i - left); |
| ix[idx] = deriv_x; |
| iy[idx] = deriv_y; |
| it[idx] = deriv_t; |
| } |
| } |
| // TODO(any): to avoid setting deriv arrays to zero for every iteration, |
| // could instead pass these two values back through function call |
| // int first_idx = (int)(y_start - top) * window_size + (int)(x_start - left); |
| // int width = window_size - ((int)(x_start - left) + (int)(left + window_size |
| // - x_end)); |
| |
| aom_free(fullpel_dx); |
| aom_free(fullpel_dy); |
| } |
| |
| // To compute eigenvalues of 2x2 matrix: Solve for lambda where |
| // Determinant(matrix - lambda*identity) == 0 |
| static void eigenvalues_2x2(const double *matrix, double *eig) { |
| const double a = 1; |
| const double b = -1 * matrix[0] - matrix[3]; |
| const double c = -1 * matrix[1] * matrix[2] + matrix[0] * matrix[3]; |
| // quadratic formula |
| const double discriminant = b * b - 4 * a * c; |
| eig[0] = (-b - sqrt(discriminant)) / (2.0 * a); |
| eig[1] = (-b + sqrt(discriminant)) / (2.0 * a); |
| // double check that eigenvalues are ordered by magnitude |
| if (fabs(eig[0]) > fabs(eig[1])) { |
| double tmp = eig[0]; |
| eig[0] = eig[1]; |
| eig[1] = tmp; |
| } |
| } |
| |
| // Shi-Tomasi corner detection criteria |
| static double corner_score(const YV12_BUFFER_CONFIG *frame_to_filter, |
| const YV12_BUFFER_CONFIG *ref_frame, const int x, |
| const int y, double *i_x, double *i_y, double *i_t, |
| const int n, const int bit_depth) { |
| double eig[2]; |
| LOCALMV mv = { .row = 0, .col = 0 }; |
| // TODO(any): technically, ref_frame and i_t are not used by corner score |
| // so these could be replaced by dummy variables, |
| // or change this to spatial gradient function over window only |
| gradients_over_window(frame_to_filter, ref_frame, x, y, n, bit_depth, i_x, |
| i_y, i_t, &mv); |
| double Mres1[1] = { 0 }, Mres2[1] = { 0 }, Mres3[1] = { 0 }; |
| multiply_mat(i_x, i_x, Mres1, 1, n * n, 1); |
| multiply_mat(i_x, i_y, Mres2, 1, n * n, 1); |
| multiply_mat(i_y, i_y, Mres3, 1, n * n, 1); |
| double M[4] = { Mres1[0], Mres2[0], Mres2[0], Mres3[0] }; |
| eigenvalues_2x2(M, eig); |
| return fabs(eig[0]); |
| } |
| |
| // Finds corners in frame_to_filter |
| // For less strict requirements (i.e. more corners), decrease threshold |
| static int detect_corners(const YV12_BUFFER_CONFIG *frame_to_filter, |
| const YV12_BUFFER_CONFIG *ref_frame, |
| const int maxcorners, int *ref_corners, |
| const int bit_depth) { |
| const int frame_height = frame_to_filter->y_crop_height; |
| const int frame_width = frame_to_filter->y_crop_width; |
| // TODO(any): currently if maxcorners is decreased, then it only means |
| // corners will be omited from bottom-right of image. if maxcorners |
| // is actually used, then this algorithm would need to re-iterate |
| // and choose threshold based on that |
| assert(maxcorners == frame_height * frame_width); |
| int countcorners = 0; |
| const double threshold = 0.1; |
| double score; |
| const int n = 3; |
| double i_x[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 }; |
| double i_y[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 }; |
| double i_t[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 }; |
| const int fromedge = n; |
| double max_score = corner_score(frame_to_filter, ref_frame, fromedge, |
| fromedge, i_x, i_y, i_t, n, bit_depth); |
| // rough estimate of max corner score in image |
| for (int x = fromedge; x < frame_width - fromedge; x += 1) { |
| for (int y = fromedge; y < frame_height - fromedge; y += frame_height / 5) { |
| for (int i = 0; i < n * n; i++) { |
| i_x[i] = 0; |
| i_y[i] = 0; |
| i_t[i] = 0; |
| } |
| score = corner_score(frame_to_filter, ref_frame, x, y, i_x, i_y, i_t, n, |
| bit_depth); |
| if (score > max_score) { |
| max_score = score; |
| } |
| } |
| } |
| // score all the points and choose corners over threshold |
| for (int x = fromedge; x < frame_width - fromedge; x += 1) { |
| for (int y = fromedge; |
| (y < frame_height - fromedge) && countcorners < maxcorners; y += 1) { |
| for (int i = 0; i < n * n; i++) { |
| i_x[i] = 0; |
| i_y[i] = 0; |
| i_t[i] = 0; |
| } |
| score = corner_score(frame_to_filter, ref_frame, x, y, i_x, i_y, i_t, n, |
| bit_depth); |
| if (score > threshold * max_score) { |
| ref_corners[countcorners * 2] = x; |
| ref_corners[countcorners * 2 + 1] = y; |
| countcorners++; |
| } |
| } |
| } |
| return countcorners; |
| } |
| |
| // weights is an nxn matrix. weights is filled with a gaussian function, |
| // with independent variable: distance from the center point. |
| static void gaussian(const double sigma, const int n, const int normalize, |
| double *weights) { |
| double total_weight = 0; |
| for (int j = 0; j < n; j++) { |
| for (int i = 0; i < n; i++) { |
| double distance = sqrt(pow(n / 2 - i, 2) + pow(n / 2 - j, 2)); |
| double weight = exp(-0.5 * pow(distance / sigma, 2)); |
| weights[j * n + i] = weight; |
| total_weight += weight; |
| } |
| } |
| if (normalize == 1) { |
| for (int j = 0; j < n; j++) { |
| weights[j] = weights[j] / total_weight; |
| } |
| } |
| } |
| |
| static double convolve(const double *filter, const int *img, const int size) { |
| double result = 0; |
| for (int i = 0; i < size; i++) { |
| result += filter[i] * img[i]; |
| } |
| return result; |
| } |
| |
| // Applies a Gaussian low-pass smoothing filter to produce |
| // a corresponding lower resolution image with halved dimensions |
| static void reduce(uint8_t *img, int height, int width, int stride, |
| uint8_t *reduced_img) { |
| const int new_width = width / 2; |
| const int window_size = 5; |
| const double gaussian_filter[25] = { |
| 1. / 256, 1.0 / 64, 3. / 128, 1. / 64, 1. / 256, 1. / 64, 1. / 16, |
| 3. / 32, 1. / 16, 1. / 64, 3. / 128, 3. / 32, 9. / 64, 3. / 32, |
| 3. / 128, 1. / 64, 1. / 16, 3. / 32, 1. / 16, 1. / 64, 1. / 256, |
| 1. / 64, 3. / 128, 1. / 64, 1. / 256 |
| }; |
| // filter is 5x5 so need prev and forward 2 pixels |
| int img_section[25]; |
| for (int y = 0; y < height - 1; y += 2) { |
| for (int x = 0; x < width - 1; x += 2) { |
| int i = 0; |
| for (int yy = y - window_size / 2; yy <= y + window_size / 2; yy++) { |
| for (int xx = x - window_size / 2; xx <= x + window_size / 2; xx++) { |
| int yvalue = yy; |
| int xvalue = xx; |
| // copied pixels outside the boundary |
| if (yvalue < 0) yvalue = 0; |
| if (xvalue < 0) xvalue = 0; |
| if (yvalue >= height) yvalue = height - 1; |
| if (xvalue >= width) xvalue = width - 1; |
| img_section[i++] = img[yvalue * stride + xvalue]; |
| } |
| } |
| reduced_img[(y / 2) * new_width + (x / 2)] = (uint8_t)convolve( |
| gaussian_filter, img_section, window_size * window_size); |
| } |
| } |
| } |
| |
| static int cmpfunc(const void *a, const void *b) { |
| return (*(int *)a - *(int *)b); |
| } |
| static void filter_mvs(const MV_FILTER_TYPE mv_filter, const int frame_height, |
| const int frame_width, LOCALMV *localmvs, MV *mvs) { |
| const int n = 5; // window size |
| // for smoothing filter |
| const double gaussian_filter[25] = { |
| 1. / 256, 1. / 64, 3. / 128, 1. / 64, 1. / 256, 1. / 64, 1. / 16, |
| 3. / 32, 1. / 16, 1. / 64, 3. / 128, 3. / 32, 9. / 64, 3. / 32, |
| 3. / 128, 1. / 64, 1. / 16, 3. / 32, 1. / 16, 1. / 64, 1. / 256, |
| 1. / 64, 3. / 128, 1. / 64, 1. / 256 |
| }; |
| // for median filter |
| int mvrows[25]; |
| int mvcols[25]; |
| if (mv_filter != MV_FILTER_NONE) { |
| for (int y = 0; y < frame_height; y++) { |
| for (int x = 0; x < frame_width; x++) { |
| int center_idx = y * frame_width + x; |
| int i = 0; |
| double filtered_row = 0; |
| double filtered_col = 0; |
| for (int yy = y - n / 2; yy <= y + n / 2; yy++) { |
| for (int xx = x - n / 2; xx <= x + n / 2; xx++) { |
| int yvalue = yy; |
| int xvalue = xx; |
| // copied pixels outside the boundary |
| if (yvalue < 0) yvalue = 0; |
| if (xvalue < 0) xvalue = 0; |
| if (yvalue >= frame_height) yvalue = frame_height - 1; |
| if (xvalue >= frame_width) xvalue = frame_width - 1; |
| int index = yvalue * frame_width + xvalue; |
| if (mv_filter == MV_FILTER_SMOOTH) { |
| filtered_row += mvs[index].row * gaussian_filter[i]; |
| filtered_col += mvs[index].col * gaussian_filter[i]; |
| } else if (mv_filter == MV_FILTER_MEDIAN) { |
| mvrows[i] = mvs[index].row; |
| mvcols[i] = mvs[index].col; |
| } |
| i++; |
| } |
| } |
| |
| MV mv = mvs[center_idx]; |
| if (mv_filter == MV_FILTER_SMOOTH) { |
| mv.row = (int16_t)filtered_row; |
| mv.col = (int16_t)filtered_col; |
| } else if (mv_filter == MV_FILTER_MEDIAN) { |
| qsort(mvrows, 25, sizeof(mv.row), cmpfunc); |
| qsort(mvcols, 25, sizeof(mv.col), cmpfunc); |
| mv.row = mvrows[25 / 2]; |
| mv.col = mvcols[25 / 2]; |
| } |
| LOCALMV localmv = { .row = ((double)mv.row) / 8, |
| .col = ((double)mv.row) / 8 }; |
| localmvs[y * frame_width + x] = localmv; |
| // if mvs array is immediately updated here, then the result may |
| // propagate to other pixels. |
| } |
| } |
| for (int i = 0; i < frame_height * frame_width; i++) { |
| MV mv = { .row = (int16_t)round(8 * localmvs[i].row), |
| .col = (int16_t)round(8 * localmvs[i].col) }; |
| mvs[i] = mv; |
| } |
| } |
| } |
| |
| // Computes optical flow at a single pyramid level, |
| // using Lucas-Kanade algorithm. |
| // Modifies mvs array. |
| static void lucas_kanade(const YV12_BUFFER_CONFIG *from_frame, |
| const YV12_BUFFER_CONFIG *to_frame, const int level, |
| const LK_PARAMS *lk_params, const int num_ref_corners, |
| int *ref_corners, const int mv_stride, |
| const int bit_depth, LOCALMV *mvs) { |
| assert(lk_params->window_size > 0 && lk_params->window_size % 2 == 0); |
| const int n = lk_params->window_size; |
| // algorithm is sensitive to window size |
| double *i_x = (double *)aom_malloc(n * n * sizeof(*i_x)); |
| double *i_y = (double *)aom_malloc(n * n * sizeof(*i_y)); |
| double *i_t = (double *)aom_malloc(n * n * sizeof(*i_t)); |
| double *weights = (double *)aom_malloc(n * n * sizeof(*weights)); |
| if (!i_x || !i_y || !i_t || !weights) goto free_lk_buf; |
| |
| const int expand_multiplier = (int)pow(2, level); |
| double sigma = 0.2 * n; |
| // normalizing doesn't really affect anything since it's applied |
| // to every component of M and b |
| gaussian(sigma, n, 0, weights); |
| for (int i = 0; i < num_ref_corners; i++) { |
| const double x_coord = 1.0 * ref_corners[i * 2] / expand_multiplier; |
| const double y_coord = 1.0 * ref_corners[i * 2 + 1] / expand_multiplier; |
| int highres_x = ref_corners[i * 2]; |
| int highres_y = ref_corners[i * 2 + 1]; |
| int mv_idx = highres_y * (mv_stride) + highres_x; |
| LOCALMV mv_old = mvs[mv_idx]; |
| mv_old.row = mv_old.row / expand_multiplier; |
| mv_old.col = mv_old.col / expand_multiplier; |
| // using this instead of memset, since it's not completely |
| // clear if zero memset works on double arrays |
| for (int j = 0; j < n * n; j++) { |
| i_x[j] = 0; |
| i_y[j] = 0; |
| i_t[j] = 0; |
| } |
| gradients_over_window(from_frame, to_frame, x_coord, y_coord, n, bit_depth, |
| i_x, i_y, i_t, &mv_old); |
| double Mres1[1] = { 0 }, Mres2[1] = { 0 }, Mres3[1] = { 0 }; |
| double bres1[1] = { 0 }, bres2[1] = { 0 }; |
| for (int j = 0; j < n * n; j++) { |
| Mres1[0] += weights[j] * i_x[j] * i_x[j]; |
| Mres2[0] += weights[j] * i_x[j] * i_y[j]; |
| Mres3[0] += weights[j] * i_y[j] * i_y[j]; |
| bres1[0] += weights[j] * i_x[j] * i_t[j]; |
| bres2[0] += weights[j] * i_y[j] * i_t[j]; |
| } |
| double M[4] = { Mres1[0], Mres2[0], Mres2[0], Mres3[0] }; |
| double b[2] = { -1 * bres1[0], -1 * bres2[0] }; |
| double eig[2] = { 1, 1 }; |
| eigenvalues_2x2(M, eig); |
| double threshold = 0.1; |
| if (fabs(eig[0]) > threshold) { |
| // if M is not invertible, then displacement |
| // will default to zeros |
| double u[2] = { 0, 0 }; |
| linsolve(2, M, 2, b, u); |
| int mult = 1; |
| if (level != 0) |
| mult = expand_multiplier; // mv doubles when resolution doubles |
| LOCALMV mv = { .row = (mult * (u[0] + mv_old.row)), |
| .col = (mult * (u[1] + mv_old.col)) }; |
| mvs[mv_idx] = mv; |
| mvs[mv_idx] = mv; |
| } |
| } |
| free_lk_buf: |
| aom_free(weights); |
| aom_free(i_t); |
| aom_free(i_x); |
| aom_free(i_y); |
| } |
| |
| // Warp the src_frame to warper_frame according to mvs. |
| // mvs point to src_frame |
| static void warp_back_frame(YV12_BUFFER_CONFIG *warped_frame, |
| const YV12_BUFFER_CONFIG *src_frame, |
| const LOCALMV *mvs, int mv_stride) { |
| int w, h; |
| const int fw = src_frame->y_crop_width; |
| const int fh = src_frame->y_crop_height; |
| const int src_fs = src_frame->y_stride, warped_fs = warped_frame->y_stride; |
| const uint8_t *src_buf = src_frame->y_buffer; |
| uint8_t *warped_buf = warped_frame->y_buffer; |
| double temp; |
| for (h = 0; h < fh; h++) { |
| for (w = 0; w < fw; w++) { |
| double cord_x = (double)w + mvs[h * mv_stride + w].col; |
| double cord_y = (double)h + mvs[h * mv_stride + w].row; |
| cord_x = fclamp(cord_x, 0, (double)(fw - 1)); |
| cord_y = fclamp(cord_y, 0, (double)(fh - 1)); |
| const int floorx = (int)floor(cord_x); |
| const int floory = (int)floor(cord_y); |
| const double fracx = cord_x - (double)floorx; |
| const double fracy = cord_y - (double)floory; |
| |
| temp = 0; |
| for (int hh = 0; hh < 2; hh++) { |
| const double weighth = hh ? (fracy) : (1 - fracy); |
| for (int ww = 0; ww < 2; ww++) { |
| const double weightw = ww ? (fracx) : (1 - fracx); |
| int y = floory + hh; |
| int x = floorx + ww; |
| y = clamp(y, 0, fh - 1); |
| x = clamp(x, 0, fw - 1); |
| temp += (double)src_buf[y * src_fs + x] * weightw * weighth; |
| } |
| } |
| warped_buf[h * warped_fs + w] = (uint8_t)round(temp); |
| } |
| } |
| } |
| |
| // Same as warp_back_frame, but using a better interpolation filter. |
| static void warp_back_frame_intp(YV12_BUFFER_CONFIG *warped_frame, |
| const YV12_BUFFER_CONFIG *src_frame, |
| const LOCALMV *mvs, int mv_stride) { |
| int w, h; |
| const int fw = src_frame->y_crop_width; |
| const int fh = src_frame->y_crop_height; |
| const int warped_fs = warped_frame->y_stride; |
| uint8_t *warped_buf = warped_frame->y_buffer; |
| const int blk = 2; |
| uint8_t temp_blk[4]; |
| |
| const int is_intrabc = 0; // Is intra-copied? |
| const int is_high_bitdepth = is_frame_high_bitdepth(src_frame); |
| const int subsampling_x = 0, subsampling_y = 0; // for y-buffer |
| const int_interpfilters interp_filters = |
| av1_broadcast_interp_filter(MULTITAP_SHARP2); |
| const int plane = 0; // y-plane |
| const struct buf_2d ref_buf2 = { NULL, src_frame->y_buffer, |
| src_frame->y_crop_width, |
| src_frame->y_crop_height, |
| src_frame->y_stride }; |
| const int bit_depth = src_frame->bit_depth; |
| struct scale_factors scale; |
| av1_setup_scale_factors_for_frame( |
| &scale, src_frame->y_crop_width, src_frame->y_crop_height, |
| src_frame->y_crop_width, src_frame->y_crop_height); |
| |
| for (h = 0; h < fh; h++) { |
| for (w = 0; w < fw; w++) { |
| InterPredParams inter_pred_params; |
| av1_init_inter_params(&inter_pred_params, blk, blk, h, w, subsampling_x, |
| subsampling_y, bit_depth, is_high_bitdepth, |
| is_intrabc, &scale, &ref_buf2, interp_filters); |
| inter_pred_params.interp_filter_params[0] = |
| &av1_interp_filter_params_list[interp_filters.as_filters.x_filter]; |
| inter_pred_params.interp_filter_params[1] = |
| &av1_interp_filter_params_list[interp_filters.as_filters.y_filter]; |
| inter_pred_params.conv_params = get_conv_params(0, plane, bit_depth); |
| MV newmv = { .row = (int16_t)round((mvs[h * mv_stride + w].row) * 8), |
| .col = (int16_t)round((mvs[h * mv_stride + w].col) * 8) }; |
| av1_enc_build_one_inter_predictor(temp_blk, blk, &newmv, |
| &inter_pred_params); |
| warped_buf[h * warped_fs + w] = temp_blk[0]; |
| } |
| } |
| } |
| |
| #define DERIVATIVE_FILTER_LENGTH 7 |
| double filter[DERIVATIVE_FILTER_LENGTH] = { -1.0 / 60, 9.0 / 60, -45.0 / 60, 0, |
| 45.0 / 60, -9.0 / 60, 1.0 / 60 }; |
| |
| // Get gradient of the whole frame |
| static void get_frame_gradients(const YV12_BUFFER_CONFIG *from_frame, |
| const YV12_BUFFER_CONFIG *to_frame, double *ix, |
| double *iy, double *it, int grad_stride) { |
| int w, h, k, idx; |
| const int fw = from_frame->y_crop_width; |
| const int fh = from_frame->y_crop_height; |
| const int from_fs = from_frame->y_stride, to_fs = to_frame->y_stride; |
| const uint8_t *from_buf = from_frame->y_buffer; |
| const uint8_t *to_buf = to_frame->y_buffer; |
| |
| const int lh = DERIVATIVE_FILTER_LENGTH; |
| const int hleft = (lh - 1) / 2; |
| |
| for (h = 0; h < fh; h++) { |
| for (w = 0; w < fw; w++) { |
| // x |
| ix[h * grad_stride + w] = 0; |
| for (k = 0; k < lh; k++) { |
| // if we want to make this block dependent, need to extend the |
| // boundaries using other initializations. |
| idx = w + k - hleft; |
| idx = clamp(idx, 0, fw - 1); |
| ix[h * grad_stride + w] += filter[k] * 0.5 * |
| ((double)from_buf[h * from_fs + idx] + |
| (double)to_buf[h * to_fs + idx]); |
| } |
| // y |
| iy[h * grad_stride + w] = 0; |
| for (k = 0; k < lh; k++) { |
| // if we want to make this block dependent, need to extend the |
| // boundaries using other initializations. |
| idx = h + k - hleft; |
| idx = clamp(idx, 0, fh - 1); |
| iy[h * grad_stride + w] += filter[k] * 0.5 * |
| ((double)from_buf[idx * from_fs + w] + |
| (double)to_buf[idx * to_fs + w]); |
| } |
| // t |
| it[h * grad_stride + w] = |
| (double)to_buf[h * to_fs + w] - (double)from_buf[h * from_fs + w]; |
| } |
| } |
| } |
| |
| // Solve for linear equations given by the H-S method |
| static void solve_horn_schunck(const double *ix, const double *iy, |
| const double *it, int grad_stride, int width, |
| int height, const LOCALMV *init_mvs, |
| int init_mv_stride, LOCALMV *mvs, |
| int mv_stride) { |
| // TODO(bohanli): May just need to allocate the buffers once per optical flow |
| // calculation |
| int *row_pos = aom_calloc(width * height * 28, sizeof(*row_pos)); |
| int *col_pos = aom_calloc(width * height * 28, sizeof(*col_pos)); |
| double *values = aom_calloc(width * height * 28, sizeof(*values)); |
| double *mv_vec = aom_calloc(width * height * 2, sizeof(*mv_vec)); |
| double *mv_init_vec = aom_calloc(width * height * 2, sizeof(*mv_init_vec)); |
| double *temp_b = aom_calloc(width * height * 2, sizeof(*temp_b)); |
| double *b = aom_calloc(width * height * 2, sizeof(*b)); |
| if (!row_pos || !col_pos || !values || !mv_vec || !mv_init_vec || !temp_b || |
| !b) { |
| goto free_hs_solver_buf; |
| } |
| |
| // the location idx for neighboring pixels, k < 4 are the 4 direct neighbors |
| const int check_locs_y[12] = { 0, 0, -1, 1, -1, -1, 1, 1, 0, 0, -2, 2 }; |
| const int check_locs_x[12] = { -1, 1, 0, 0, -1, 1, -1, 1, -2, 2, 0, 0 }; |
| |
| int h, w, checkh, checkw, k, ret; |
| const int offset = height * width; |
| SPARSE_MTX A; |
| int c = 0; |
| const double lambda = 100; |
| |
| for (w = 0; w < width; w++) { |
| for (h = 0; h < height; h++) { |
| mv_init_vec[w * height + h] = init_mvs[h * init_mv_stride + w].col; |
| mv_init_vec[w * height + h + offset] = |
| init_mvs[h * init_mv_stride + w].row; |
| } |
| } |
| |
| // get matrix A |
| for (w = 0; w < width; w++) { |
| for (h = 0; h < height; h++) { |
| int center_num_direct = 4; |
| const int center_idx = w * height + h; |
| if (w == 0 || w == width - 1) center_num_direct--; |
| if (h == 0 || h == height - 1) center_num_direct--; |
| // diagonal entry for this row from the center pixel |
| double cor_w = center_num_direct * center_num_direct + center_num_direct; |
| row_pos[c] = center_idx; |
| col_pos[c] = center_idx; |
| values[c] = lambda * cor_w; |
| c++; |
| row_pos[c] = center_idx + offset; |
| col_pos[c] = center_idx + offset; |
| values[c] = lambda * cor_w; |
| c++; |
| // other entries from direct neighbors |
| for (k = 0; k < 4; k++) { |
| checkh = h + check_locs_y[k]; |
| checkw = w + check_locs_x[k]; |
| if (checkh < 0 || checkh >= height || checkw < 0 || checkw >= width) { |
| continue; |
| } |
| int this_idx = checkw * height + checkh; |
| int this_num_direct = 4; |
| if (checkw == 0 || checkw == width - 1) this_num_direct--; |
| if (checkh == 0 || checkh == height - 1) this_num_direct--; |
| cor_w = -center_num_direct - this_num_direct; |
| row_pos[c] = center_idx; |
| col_pos[c] = this_idx; |
| values[c] = lambda * cor_w; |
| c++; |
| row_pos[c] = center_idx + offset; |
| col_pos[c] = this_idx + offset; |
| values[c] = lambda * cor_w; |
| c++; |
| } |
| // entries from neighbors on the diagonal corners |
| for (k = 4; k < 8; k++) { |
| checkh = h + check_locs_y[k]; |
| checkw = w + check_locs_x[k]; |
| if (checkh < 0 || checkh >= height || checkw < 0 || checkw >= width) { |
| continue; |
| } |
| int this_idx = checkw * height + checkh; |
| cor_w = 2; |
| row_pos[c] = center_idx; |
| col_pos[c] = this_idx; |
| values[c] = lambda * cor_w; |
| c++; |
| row_pos[c] = center_idx + offset; |
| col_pos[c] = this_idx + offset; |
| values[c] = lambda * cor_w; |
| c++; |
| } |
| // entries from neighbors with dist of 2 |
| for (k = 8; k < 12; k++) { |
| checkh = h + check_locs_y[k]; |
| checkw = w + check_locs_x[k]; |
| if (checkh < 0 || checkh >= height || checkw < 0 || checkw >= width) { |
| continue; |
| } |
| int this_idx = checkw * height + checkh; |
| cor_w = 1; |
| row_pos[c] = center_idx; |
| col_pos[c] = this_idx; |
| values[c] = lambda * cor_w; |
| c++; |
| row_pos[c] = center_idx + offset; |
| col_pos[c] = this_idx + offset; |
| values[c] = lambda * cor_w; |
| c++; |
| } |
| } |
| } |
| ret = av1_init_sparse_mtx(row_pos, col_pos, values, c, 2 * width * height, |
| 2 * width * height, &A); |
| if (ret < 0) goto free_hs_solver_buf; |
| // subtract init mv part from b |
| av1_mtx_vect_multi_left(&A, mv_init_vec, temp_b, 2 * width * height); |
| for (int i = 0; i < 2 * width * height; i++) { |
| b[i] = -temp_b[i]; |
| } |
| av1_free_sparse_mtx_elems(&A); |
| |
| // add cross terms to A and modify b with ExEt / EyEt |
| for (w = 0; w < width; w++) { |
| for (h = 0; h < height; h++) { |
| int curidx = w * height + h; |
| // modify b |
| b[curidx] += -ix[h * grad_stride + w] * it[h * grad_stride + w]; |
| b[curidx + offset] += -iy[h * grad_stride + w] * it[h * grad_stride + w]; |
| // add cross terms to A |
| row_pos[c] = curidx; |
| col_pos[c] = curidx + offset; |
| values[c] = ix[h * grad_stride + w] * iy[h * grad_stride + w]; |
| c++; |
| row_pos[c] = curidx + offset; |
| col_pos[c] = curidx; |
| values[c] = ix[h * grad_stride + w] * iy[h * grad_stride + w]; |
| c++; |
| } |
| } |
| // Add diagonal terms to A |
| for (int i = 0; i < c; i++) { |
| if (row_pos[i] == col_pos[i]) { |
| if (row_pos[i] < offset) { |
| w = row_pos[i] / height; |
| h = row_pos[i] % height; |
| values[i] += pow(ix[h * grad_stride + w], 2); |
| } else { |
| w = (row_pos[i] - offset) / height; |
| h = (row_pos[i] - offset) % height; |
| values[i] += pow(iy[h * grad_stride + w], 2); |
| } |
| } |
| } |
| |
| ret = av1_init_sparse_mtx(row_pos, col_pos, values, c, 2 * width * height, |
| 2 * width * height, &A); |
| if (ret < 0) goto free_hs_solver_buf; |
| |
| // solve for the mvs |
| ret = av1_conjugate_gradient_sparse(&A, b, 2 * width * height, mv_vec); |
| if (ret < 0) goto free_hs_solver_buf; |
| |
| // copy mvs |
| for (w = 0; w < width; w++) { |
| for (h = 0; h < height; h++) { |
| mvs[h * mv_stride + w].col = mv_vec[w * height + h]; |
| mvs[h * mv_stride + w].row = mv_vec[w * height + h + offset]; |
| } |
| } |
| free_hs_solver_buf: |
| aom_free(row_pos); |
| aom_free(col_pos); |
| aom_free(values); |
| aom_free(mv_vec); |
| aom_free(mv_init_vec); |
| aom_free(b); |
| aom_free(temp_b); |
| av1_free_sparse_mtx_elems(&A); |
| } |
| |
| // Calculate optical flow from from_frame to to_frame using the H-S method. |
| static void horn_schunck(const YV12_BUFFER_CONFIG *from_frame, |
| const YV12_BUFFER_CONFIG *to_frame, const int level, |
| const int mv_stride, const int mv_height, |
| const int mv_width, const OPFL_PARAMS *opfl_params, |
| LOCALMV *mvs) { |
| // mvs are always on level 0, here we define two new mv arrays that is of size |
| // of this level. |
| const int fw = from_frame->y_crop_width; |
| const int fh = from_frame->y_crop_height; |
| const int factor = (int)pow(2, level); |
| int w, h, k, init_mv_stride; |
| LOCALMV *init_mvs = NULL, *refine_mvs = NULL; |
| double *ix = NULL, *iy = NULL, *it = NULL; |
| YV12_BUFFER_CONFIG temp_frame; |
| temp_frame.y_buffer = NULL; |
| if (level == 0) { |
| init_mvs = mvs; |
| init_mv_stride = mv_stride; |
| } else { |
| init_mvs = aom_calloc(fw * fh, sizeof(*mvs)); |
| if (!init_mvs) goto free_hs_buf; |
| init_mv_stride = fw; |
| for (h = 0; h < fh; h++) { |
| for (w = 0; w < fw; w++) { |
| init_mvs[h * init_mv_stride + w].row = |
| mvs[h * factor * mv_stride + w * factor].row / (double)factor; |
| init_mvs[h * init_mv_stride + w].col = |
| mvs[h * factor * mv_stride + w * factor].col / (double)factor; |
| } |
| } |
| } |
| refine_mvs = aom_calloc(fw * fh, sizeof(*mvs)); |
| if (!refine_mvs) goto free_hs_buf; |
| // temp frame for warping |
| temp_frame.y_buffer = |
| (uint8_t *)aom_calloc(fh * fw, sizeof(*temp_frame.y_buffer)); |
| if (!temp_frame.y_buffer) goto free_hs_buf; |
| temp_frame.y_crop_height = fh; |
| temp_frame.y_crop_width = fw; |
| temp_frame.y_stride = fw; |
| // gradient buffers |
| ix = aom_calloc(fw * fh, sizeof(*ix)); |
| iy = aom_calloc(fw * fh, sizeof(*iy)); |
| it = aom_calloc(fw * fh, sizeof(*it)); |
| if (!ix || !iy || !it) goto free_hs_buf; |
| // For each warping step |
| for (k = 0; k < opfl_params->warping_steps; k++) { |
| // warp from_frame with init_mv |
| if (level == 0) { |
| warp_back_frame_intp(&temp_frame, to_frame, init_mvs, init_mv_stride); |
| } else { |
| warp_back_frame(&temp_frame, to_frame, init_mvs, init_mv_stride); |
| } |
| // calculate frame gradients |
| get_frame_gradients(from_frame, &temp_frame, ix, iy, it, fw); |
| // form linear equations and solve mvs |
| solve_horn_schunck(ix, iy, it, fw, fw, fh, init_mvs, init_mv_stride, |
| refine_mvs, fw); |
| // update init_mvs |
| for (h = 0; h < fh; h++) { |
| for (w = 0; w < fw; w++) { |
| init_mvs[h * init_mv_stride + w].col += refine_mvs[h * fw + w].col; |
| init_mvs[h * init_mv_stride + w].row += refine_mvs[h * fw + w].row; |
| } |
| } |
| } |
| // copy back the mvs if needed |
| if (level != 0) { |
| for (h = 0; h < mv_height; h++) { |
| for (w = 0; w < mv_width; w++) { |
| mvs[h * mv_stride + w].row = |
| init_mvs[h / factor * init_mv_stride + w / factor].row * |
| (double)factor; |
| mvs[h * mv_stride + w].col = |
| init_mvs[h / factor * init_mv_stride + w / factor].col * |
| (double)factor; |
| } |
| } |
| } |
| free_hs_buf: |
| if (level != 0) aom_free(init_mvs); |
| aom_free(refine_mvs); |
| aom_free(temp_frame.y_buffer); |
| aom_free(ix); |
| aom_free(iy); |
| aom_free(it); |
| } |
| |
| // Apply optical flow iteratively at each pyramid level |
| static void pyramid_optical_flow(const YV12_BUFFER_CONFIG *from_frame, |
| const YV12_BUFFER_CONFIG *to_frame, |
| const int bit_depth, |
| const OPFL_PARAMS *opfl_params, |
| const OPTFLOW_METHOD method, LOCALMV *mvs) { |
| assert(opfl_params->pyramid_levels > 0 && |
| opfl_params->pyramid_levels <= MAX_PYRAMID_LEVELS); |
| int levels = opfl_params->pyramid_levels; |
| const int frame_height = from_frame->y_crop_height; |
| const int frame_width = from_frame->y_crop_width; |
| if ((frame_height / pow(2.0, levels - 1) < 50 || |
| frame_height / pow(2.0, levels - 1) < 50) && |
| levels > 1) |
| levels = levels - 1; |
| uint8_t *images1[MAX_PYRAMID_LEVELS] = { NULL }; |
| uint8_t *images2[MAX_PYRAMID_LEVELS] = { NULL }; |
| int *ref_corners = NULL; |
| |
| images1[0] = from_frame->y_buffer; |
| images2[0] = to_frame->y_buffer; |
| YV12_BUFFER_CONFIG *buffers1 = aom_malloc(levels * sizeof(*buffers1)); |
| YV12_BUFFER_CONFIG *buffers2 = aom_malloc(levels * sizeof(*buffers2)); |
| if (!buffers1 || !buffers2) goto free_pyramid_buf; |
| buffers1[0] = *from_frame; |
| buffers2[0] = *to_frame; |
| int fw = frame_width; |
| int fh = frame_height; |
| for (int i = 1; i < levels; i++) { |
| // TODO(bohanli): may need to extend buffers for better interpolation SIMD |
| images1[i] = (uint8_t *)aom_calloc(fh / 2 * fw / 2, sizeof(*images1[i])); |
| images2[i] = (uint8_t *)aom_calloc(fh / 2 * fw / 2, sizeof(*images2[i])); |
| if (!images1[i] || !images2[i]) goto free_pyramid_buf; |
| int stride; |
| if (i == 1) |
| stride = from_frame->y_stride; |
| else |
| stride = fw; |
| reduce(images1[i - 1], fh, fw, stride, images1[i]); |
| reduce(images2[i - 1], fh, fw, stride, images2[i]); |
| fh /= 2; |
| fw /= 2; |
| YV12_BUFFER_CONFIG a = { .y_buffer = images1[i], |
| .y_crop_width = fw, |
| .y_crop_height = fh, |
| .y_stride = fw }; |
| YV12_BUFFER_CONFIG b = { .y_buffer = images2[i], |
| .y_crop_width = fw, |
| .y_crop_height = fh, |
| .y_stride = fw }; |
| buffers1[i] = a; |
| buffers2[i] = b; |
| } |
| // Compute corners for specific frame |
| int num_ref_corners = 0; |
| if (is_sparse(opfl_params)) { |
| int maxcorners = from_frame->y_crop_width * from_frame->y_crop_height; |
| ref_corners = aom_malloc(maxcorners * 2 * sizeof(*ref_corners)); |
| if (!ref_corners) goto free_pyramid_buf; |
| num_ref_corners = detect_corners(from_frame, to_frame, maxcorners, |
| ref_corners, bit_depth); |
| } |
| const int stop_level = 0; |
| for (int i = levels - 1; i >= stop_level; i--) { |
| if (method == LUCAS_KANADE) { |
| assert(is_sparse(opfl_params)); |
| lucas_kanade(&buffers1[i], &buffers2[i], i, opfl_params->lk_params, |
| num_ref_corners, ref_corners, buffers1[0].y_crop_width, |
| bit_depth, mvs); |
| } else if (method == HORN_SCHUNCK) { |
| assert(!is_sparse(opfl_params)); |
| horn_schunck(&buffers1[i], &buffers2[i], i, buffers1[0].y_crop_width, |
| buffers1[0].y_crop_height, buffers1[0].y_crop_width, |
| opfl_params, mvs); |
| } |
| } |
| free_pyramid_buf: |
| for (int i = 1; i < levels; i++) { |
| aom_free(images1[i]); |
| aom_free(images2[i]); |
| } |
| aom_free(ref_corners); |
| aom_free(buffers1); |
| aom_free(buffers2); |
| } |
| // Computes optical flow by applying algorithm at |
| // multiple pyramid levels of images (lower-resolution, smoothed images) |
| // This accounts for larger motions. |
| // Inputs: |
| // from_frame Frame buffer. |
| // to_frame: Frame buffer. MVs point from_frame -> to_frame. |
| // from_frame_idx: Index of from_frame. |
| // to_frame_idx: Index of to_frame. Return all zero MVs when idx are equal. |
| // bit_depth: |
| // opfl_params: contains algorithm-specific parameters. |
| // mv_filter: MV_FILTER_NONE, MV_FILTER_SMOOTH, or MV_FILTER_MEDIAN. |
| // method: LUCAS_KANADE, HORN_SCHUNCK |
| // mvs: pointer to MVs. Contains initialization, and modified |
| // based on optical flow. Must have |
| // dimensions = from_frame->y_crop_width * from_frame->y_crop_height |
| void av1_optical_flow(const YV12_BUFFER_CONFIG *from_frame, |
| const YV12_BUFFER_CONFIG *to_frame, |
| const int from_frame_idx, const int to_frame_idx, |
| const int bit_depth, const OPFL_PARAMS *opfl_params, |
| const MV_FILTER_TYPE mv_filter, |
| const OPTFLOW_METHOD method, MV *mvs) { |
| const int frame_height = from_frame->y_crop_height; |
| const int frame_width = from_frame->y_crop_width; |
| // TODO(any): deal with the case where frames are not of the same dimensions |
| assert(frame_height == to_frame->y_crop_height && |
| frame_width == to_frame->y_crop_width); |
| if (from_frame_idx == to_frame_idx) { |
| // immediately return all zero mvs when frame indices are equal |
| for (int yy = 0; yy < frame_height; yy++) { |
| for (int xx = 0; xx < frame_width; xx++) { |
| MV mv = { .row = 0, .col = 0 }; |
| mvs[yy * frame_width + xx] = mv; |
| } |
| } |
| return; |
| } |
| |
| // Initialize double mvs based on input parameter mvs array |
| LOCALMV *localmvs = |
| aom_malloc(frame_height * frame_width * sizeof(*localmvs)); |
| if (!localmvs) return; |
| |
| filter_mvs(MV_FILTER_SMOOTH, frame_height, frame_width, localmvs, mvs); |
| |
| for (int i = 0; i < frame_width * frame_height; i++) { |
| MV mv = mvs[i]; |
| LOCALMV localmv = { .row = ((double)mv.row) / 8, |
| .col = ((double)mv.col) / 8 }; |
| localmvs[i] = localmv; |
| } |
| // Apply optical flow algorithm |
| pyramid_optical_flow(from_frame, to_frame, bit_depth, opfl_params, method, |
| localmvs); |
| |
| // Update original mvs array |
| for (int j = 0; j < frame_height; j++) { |
| for (int i = 0; i < frame_width; i++) { |
| int idx = j * frame_width + i; |
| if (j + localmvs[idx].row < 0 || j + localmvs[idx].row >= frame_height || |
| i + localmvs[idx].col < 0 || i + localmvs[idx].col >= frame_width) { |
| continue; |
| } |
| MV mv = { .row = (int16_t)round(8 * localmvs[idx].row), |
| .col = (int16_t)round(8 * localmvs[idx].col) }; |
| mvs[idx] = mv; |
| } |
| } |
| |
| filter_mvs(mv_filter, frame_height, frame_width, localmvs, mvs); |
| |
| aom_free(localmvs); |
| } |
| #endif |