/* * 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 #include #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