Rework global motion estimation code
diff --git a/aom/exports_test b/aom/exports_test index 452a532..fc708c1 100644 --- a/aom/exports_test +++ b/aom/exports_test
@@ -1,4 +1,5 @@ text aom_copy_metadata_to_frame_buffer text aom_dsp_rtcd +text aom_fit_motion_model text aom_remove_metadata_from_frame_buffer text aom_scale_rtcd
diff --git a/aom_dsp/aom_dsp.cmake b/aom_dsp/aom_dsp.cmake index 4cc8f57..8eb939c 100644 --- a/aom_dsp/aom_dsp.cmake +++ b/aom_dsp/aom_dsp.cmake
@@ -185,6 +185,7 @@ "${AOM_ROOT}/aom_dsp/fwd_txfm.c" "${AOM_ROOT}/aom_dsp/grain_table.c" "${AOM_ROOT}/aom_dsp/grain_table.h" + "${AOM_ROOT}/aom_dsp/linalg.c" "${AOM_ROOT}/aom_dsp/noise_model.c" "${AOM_ROOT}/aom_dsp/noise_model.h" "${AOM_ROOT}/aom_dsp/noise_util.c" @@ -299,6 +300,23 @@ "${AOM_ROOT}/aom_dsp/mips/variance_msa.c" "${AOM_ROOT}/aom_dsp/mips/sub_pixel_variance_msa.c") + # Flow estimation library + list( + APPEND + AOM_DSP_ENCODER_SOURCES + "${AOM_ROOT}/aom_dsp/flow_estimation/corner_detect.c" + "${AOM_ROOT}/aom_dsp/flow_estimation/corner_match.c" + "${AOM_ROOT}/aom_dsp/flow_estimation/disflow.c" + "${AOM_ROOT}/aom_dsp/flow_estimation/flow_estimation.c" + "${AOM_ROOT}/aom_dsp/flow_estimation/pyramid.c" + "${AOM_ROOT}/aom_dsp/flow_estimation/ransac.c") + + list(APPEND AOM_DSP_ENCODER_INTRIN_SSE4_1 + "${AOM_ROOT}/aom_dsp/flow_estimation/x86/corner_match_sse4.c") + + list(APPEND AOM_DSP_ENCODER_INTRIN_AVX2 + "${AOM_ROOT}/aom_dsp/flow_estimation/x86/corner_match_avx2.c") + if(CONFIG_INTERNAL_STATS) list(APPEND AOM_DSP_ENCODER_SOURCES "${AOM_ROOT}/aom_dsp/fastssim.c" "${AOM_ROOT}/aom_dsp/psnrhvs.c" "${AOM_ROOT}/aom_dsp/ssim.c"
diff --git a/aom_dsp/aom_dsp_rtcd_defs.pl b/aom_dsp/aom_dsp_rtcd_defs.pl index 255fc54..e3a3fc5 100755 --- a/aom_dsp/aom_dsp_rtcd_defs.pl +++ b/aom_dsp/aom_dsp_rtcd_defs.pl
@@ -1851,4 +1851,10 @@ specialize qw/aom_highbd_comp_mask_pred sse2 avx2/; } # CONFIG_AV1_ENCODER +# Flow estimation library +if (aom_config("CONFIG_AV1_ENCODER") eq "yes") { + add_proto qw/double aom_compute_cross_correlation/, "unsigned char *im1, int stride1, int x1, int y1, unsigned char *im2, int stride2, int x2, int y2"; + specialize qw/aom_compute_cross_correlation sse4_1 avx2/; +} + 1;
diff --git a/aom_dsp/flow_estimation/corner_detect.c b/aom_dsp/flow_estimation/corner_detect.c new file mode 100644 index 0000000..5af942a --- /dev/null +++ b/aom_dsp/flow_estimation/corner_detect.c
@@ -0,0 +1,58 @@ +/* + * Copyright (c) 2021, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#include <stdlib.h> +#include <stdio.h> +#include <memory.h> +#include <math.h> +#include <assert.h> + +#include "third_party/fastfeat/fast.h" + +#include "aom_dsp/flow_estimation/corner_detect.h" +#include "aom_dsp/flow_estimation/flow_estimation.h" +#include "aom_dsp/flow_estimation/pyramid.h" +#include "aom_mem/aom_mem.h" + +// Fast_9 wrapper +#define FAST_BARRIER 18 +int aom_fast_corner_detect(unsigned char *buf, int width, int height, + int stride, int *points, int max_points) { + int num_points; + xy *const frm_corners_xy = aom_fast9_detect_nonmax(buf, width, height, stride, + FAST_BARRIER, &num_points); + num_points = (num_points <= max_points ? num_points : max_points); + if (num_points > 0 && frm_corners_xy) { + memcpy(points, frm_corners_xy, sizeof(*frm_corners_xy) * num_points); + free(frm_corners_xy); + return num_points; + } + free(frm_corners_xy); + return 0; +} + +void aom_find_corners_in_frame(YV12_BUFFER_CONFIG *frm, int bit_depth) { + if (!frm->y_pyramid) { + frm->y_pyramid = aom_compute_pyramid(frm, bit_depth, MAX_PYRAMID_LEVELS); + assert(frm->y_pyramid); + } + ImagePyramid *pyr = frm->y_pyramid; + + unsigned char *buffer = pyr->level_buffer + pyr->level_loc[0]; + int width = pyr->widths[0]; + int height = pyr->heights[0]; + int stride = pyr->strides[0]; + + frm->corners = aom_malloc(2 * MAX_CORNERS * sizeof(*frm->corners)); + frm->num_corners = aom_fast_corner_detect(buffer, width, height, stride, + frm->corners, MAX_CORNERS); +}
diff --git a/av1/encoder/corner_detect.h b/aom_dsp/flow_estimation/corner_detect.h similarity index 64% rename from av1/encoder/corner_detect.h rename to aom_dsp/flow_estimation/corner_detect.h index c4a5668..fd8fe28 100644 --- a/av1/encoder/corner_detect.h +++ b/aom_dsp/flow_estimation/corner_detect.h
@@ -10,14 +10,26 @@ * aomedia.org/license/patent-license/. */ -#ifndef AOM_AV1_ENCODER_CORNER_DETECT_H_ -#define AOM_AV1_ENCODER_CORNER_DETECT_H_ +#ifndef AOM_FLOW_ESTIMATION_CORNER_DETECT_H_ +#define AOM_FLOW_ESTIMATION_CORNER_DETECT_H_ + +#include "aom_scale/yv12config.h" #include <stdio.h> #include <stdlib.h> #include <memory.h> -int av1_fast_corner_detect(unsigned char *buf, int width, int height, +#ifdef __cplusplus +extern "C" { +#endif + +int aom_fast_corner_detect(unsigned char *buf, int width, int height, int stride, int *points, int max_points); -#endif // AOM_AV1_ENCODER_CORNER_DETECT_H_ +void aom_find_corners_in_frame(YV12_BUFFER_CONFIG *frm, int bit_depth); + +#ifdef __cplusplus +} +#endif + +#endif // AOM_FLOW_ESTIMATION_CORNER_DETECT_H_
diff --git a/aom_dsp/flow_estimation/corner_match.c b/aom_dsp/flow_estimation/corner_match.c new file mode 100644 index 0000000..c8a41fa --- /dev/null +++ b/aom_dsp/flow_estimation/corner_match.c
@@ -0,0 +1,315 @@ +/* + * Copyright (c) 2021, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#include <stdlib.h> +#include <memory.h> +#include <math.h> + +#include "config/aom_dsp_rtcd.h" + +#include "aom_ports/system_state.h" +#include "aom_dsp/flow_estimation/corner_detect.h" +#include "aom_dsp/flow_estimation/corner_match.h" +#include "aom_dsp/flow_estimation/ransac.h" +#include "aom_mem/aom_mem.h" + +#define SEARCH_SZ 9 +#define SEARCH_SZ_BY2 ((SEARCH_SZ - 1) / 2) + +#define THRESHOLD_NCC 0.75 + +/* Compute var(im) * MATCH_SZ_SQ over a MATCH_SZ by MATCH_SZ window of im, + centered at (x, y). +*/ +static double compute_variance(unsigned char *im, int stride, int x, int y) { + int sum = 0; + int sumsq = 0; + int var; + int i, j; + for (i = 0; i < MATCH_SZ; ++i) + for (j = 0; j < MATCH_SZ; ++j) { + sum += im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)]; + sumsq += im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)] * + im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)]; + } + var = sumsq * MATCH_SZ_SQ - sum * sum; + return (double)var; +} + +/* Compute corr(im1, im2) * MATCH_SZ * stddev(im1), where the + correlation/standard deviation are taken over MATCH_SZ by MATCH_SZ windows + of each image, centered at (x1, y1) and (x2, y2) respectively. +*/ +double aom_compute_cross_correlation_c(unsigned char *im1, int stride1, int x1, + int y1, unsigned char *im2, int stride2, + int x2, int y2) { + int v1, v2; + int sum1 = 0; + int sum2 = 0; + int sumsq2 = 0; + int cross = 0; + int var2, cov; + int i, j; + for (i = 0; i < MATCH_SZ; ++i) + for (j = 0; j < MATCH_SZ; ++j) { + v1 = im1[(i + y1 - MATCH_SZ_BY2) * stride1 + (j + x1 - MATCH_SZ_BY2)]; + v2 = im2[(i + y2 - MATCH_SZ_BY2) * stride2 + (j + x2 - MATCH_SZ_BY2)]; + sum1 += v1; + sum2 += v2; + sumsq2 += v2 * v2; + cross += v1 * v2; + } + var2 = sumsq2 * MATCH_SZ_SQ - sum2 * sum2; + cov = cross * MATCH_SZ_SQ - sum1 * sum2; + aom_clear_system_state(); + return cov / sqrt((double)var2); +} + +static int is_eligible_point(int pointx, int pointy, int width, int height) { + return (pointx >= MATCH_SZ_BY2 && pointy >= MATCH_SZ_BY2 && + pointx + MATCH_SZ_BY2 < width && pointy + MATCH_SZ_BY2 < height); +} + +static int is_eligible_distance(int point1x, int point1y, int point2x, + int point2y, int width, int height) { + const int thresh = (width < height ? height : width) >> 4; + return ((point1x - point2x) * (point1x - point2x) + + (point1y - point2y) * (point1y - point2y)) <= thresh * thresh; +} + +static void improve_correspondence(unsigned char *frm, unsigned char *ref, + int width, int height, int frm_stride, + int ref_stride, + Correspondence *correspondences, + int num_correspondences) { + int i; + for (i = 0; i < num_correspondences; ++i) { + int x, y, best_x = 0, best_y = 0; + double best_match_ncc = 0.0; + int x0 = (int)correspondences[i].x; + int y0 = (int)correspondences[i].y; + int rx0 = (int)correspondences[i].rx; + int ry0 = (int)correspondences[i].ry; + for (y = -SEARCH_SZ_BY2; y <= SEARCH_SZ_BY2; ++y) { + for (x = -SEARCH_SZ_BY2; x <= SEARCH_SZ_BY2; ++x) { + double match_ncc; + if (!is_eligible_point(rx0 + x, ry0 + y, width, height)) continue; + if (!is_eligible_distance(x0, y0, rx0 + x, ry0 + y, width, height)) + continue; + match_ncc = aom_compute_cross_correlation(frm, frm_stride, x0, y0, ref, + ref_stride, rx0 + x, ry0 + y); + if (match_ncc > best_match_ncc) { + best_match_ncc = match_ncc; + best_y = y; + best_x = x; + } + } + } + correspondences[i].rx += best_x; + correspondences[i].ry += best_y; + } + for (i = 0; i < num_correspondences; ++i) { + int x, y, best_x = 0, best_y = 0; + double best_match_ncc = 0.0; + int x0 = (int)correspondences[i].x; + int y0 = (int)correspondences[i].y; + int rx0 = (int)correspondences[i].rx; + int ry0 = (int)correspondences[i].ry; + for (y = -SEARCH_SZ_BY2; y <= SEARCH_SZ_BY2; ++y) + for (x = -SEARCH_SZ_BY2; x <= SEARCH_SZ_BY2; ++x) { + double match_ncc; + if (!is_eligible_point(x0 + x, y0 + y, width, height)) continue; + if (!is_eligible_distance(x0 + x, y0 + y, rx0, ry0, width, height)) + continue; + match_ncc = aom_compute_cross_correlation( + ref, ref_stride, rx0, ry0, frm, frm_stride, x0 + x, y0 + y); + if (match_ncc > best_match_ncc) { + best_match_ncc = match_ncc; + best_y = y; + best_x = x; + } + } + correspondences[i].x += best_x; + correspondences[i].y += best_y; + } +} + +static INLINE int determine_correspondence(unsigned char *src, int *src_corners, + int num_src_corners, + unsigned char *ref, int *ref_corners, + int num_ref_corners, int width, + int height, int src_stride, + int ref_stride, + Correspondence *correspondences) { + // TODO(sarahparker) Improve this to include 2-way match + int i, j; + int num_correspondences = 0; + for (i = 0; i < num_src_corners; ++i) { + double best_match_ncc = 0.0; + double template_norm; + int best_match_j = -1; + if (!is_eligible_point(src_corners[2 * i], src_corners[2 * i + 1], width, + height)) + continue; + for (j = 0; j < num_ref_corners; ++j) { + double match_ncc; + if (!is_eligible_point(ref_corners[2 * j], ref_corners[2 * j + 1], width, + height)) + continue; + if (!is_eligible_distance(src_corners[2 * i], src_corners[2 * i + 1], + ref_corners[2 * j], ref_corners[2 * j + 1], + width, height)) + continue; + match_ncc = aom_compute_cross_correlation( + src, src_stride, src_corners[2 * i], src_corners[2 * i + 1], ref, + ref_stride, ref_corners[2 * j], ref_corners[2 * j + 1]); + if (match_ncc > best_match_ncc) { + best_match_ncc = match_ncc; + best_match_j = j; + } + } + // Note: We want to test if the best correlation is >= THRESHOLD_NCC, + // but need to account for the normalization in + // aom_compute_cross_correlation. + template_norm = compute_variance(src, src_stride, src_corners[2 * i], + src_corners[2 * i + 1]); + if (best_match_ncc > THRESHOLD_NCC * sqrt(template_norm)) { + correspondences[num_correspondences].x = src_corners[2 * i]; + correspondences[num_correspondences].y = src_corners[2 * i + 1]; + correspondences[num_correspondences].rx = ref_corners[2 * best_match_j]; + correspondences[num_correspondences].ry = + ref_corners[2 * best_match_j + 1]; + num_correspondences++; + } + } + improve_correspondence(src, ref, width, height, src_stride, ref_stride, + correspondences, num_correspondences); + return num_correspondences; +} + +CorrespondenceList *aom_compute_corner_match(YV12_BUFFER_CONFIG *src, + YV12_BUFFER_CONFIG *ref, + int bit_depth) { + // Ensure that all relevant per-frame data is available + if (!src->y_pyramid) { + src->y_pyramid = aom_compute_pyramid(src, bit_depth, MAX_PYRAMID_LEVELS); + assert(src->y_pyramid); + } + if (!ref->y_pyramid) { + ref->y_pyramid = aom_compute_pyramid(ref, bit_depth, MAX_PYRAMID_LEVELS); + assert(ref->y_pyramid); + } + if (!src->corners) { + aom_find_corners_in_frame(src, bit_depth); + } + if (!ref->corners) { + aom_find_corners_in_frame(ref, bit_depth); + } + + ImagePyramid *src_pyr = src->y_pyramid; + + unsigned char *src_buffer = src_pyr->level_buffer + src_pyr->level_loc[0]; + int src_width = src_pyr->widths[0]; + int src_height = src_pyr->heights[0]; + int src_stride = src_pyr->strides[0]; + + ImagePyramid *ref_pyr = ref->y_pyramid; + + unsigned char *ref_buffer = ref_pyr->level_buffer + ref_pyr->level_loc[0]; + int ref_stride = ref_pyr->strides[0]; + assert(ref_pyr->widths[0] == src_width); + assert(ref_pyr->heights[0] == src_height); + + // Compute correspondences + CorrespondenceList *list = aom_malloc(sizeof(CorrespondenceList)); + list->correspondences = (Correspondence *)aom_malloc( + src->num_corners * sizeof(*list->correspondences)); + list->num_correspondences = determine_correspondence( + src_buffer, src->corners, src->num_corners, ref_buffer, ref->corners, + ref->num_corners, src_width, src_height, src_stride, ref_stride, + list->correspondences); + return list; +} + +bool aom_fit_global_model_to_correspondences(CorrespondenceList *corrs, + TransformationType type, + MotionModel *params_by_motion, + int num_motions) { + int num_correspondences = corrs->num_correspondences; + + ransac(corrs->correspondences, num_correspondences, type, params_by_motion, + num_motions); + + // Set num_inliers = 0 for motions with too few inliers so they are ignored. + for (int i = 0; i < num_motions; ++i) { + if (params_by_motion[i].num_inliers < + MIN_INLIER_PROB * num_correspondences || + num_correspondences == 0) { + params_by_motion[i].num_inliers = 0; + } + } + + // Return true if any one of the motions has inliers. + for (int i = 0; i < num_motions; ++i) { + if (params_by_motion[i].num_inliers > 0) return true; + } + return false; +} + +bool aom_fit_local_model_to_correspondences(CorrespondenceList *corrs, + PixelRect *rect, + TransformationType type, + double *mat) { + int width = rect_height(rect); + int height = rect_width(rect); + int num_points = width * height; + + // TODO(rachelbarker): Downsample if num_points is > some threshold? + double *pts1 = aom_malloc(num_points * 2 * sizeof(double)); + double *pts2 = aom_malloc(num_points * 2 * sizeof(double)); + int point_index = 0; + + for (int i = 0; i < corrs->num_correspondences; i++) { + Correspondence *corr = &corrs->correspondences[i]; + int x = (int)corr->x; + int y = (int)corr->y; + if (is_inside_rect(x, y, rect)) { + pts1[2 * point_index + 0] = corr->x; + pts1[2 * point_index + 1] = corr->y; + pts2[2 * point_index + 0] = corr->rx; + pts2[2 * point_index + 1] = corr->ry; + point_index++; + } + } + assert(point_index <= num_points); + + num_points = point_index; + + bool result; + if (num_points < 4) { + // Too few points to fit a model + result = false; + } else { + result = aom_fit_motion_model(type, num_points, pts1, pts2, mat); + } + + aom_free(pts1); + aom_free(pts2); + return result; +} + +void aom_free_correspondence_list(CorrespondenceList *list) { + if (list) { + aom_free(list->correspondences); + aom_free(list); + } +}
diff --git a/aom_dsp/flow_estimation/corner_match.h b/aom_dsp/flow_estimation/corner_match.h new file mode 100644 index 0000000..fb8f2a0 --- /dev/null +++ b/aom_dsp/flow_estimation/corner_match.h
@@ -0,0 +1,51 @@ +/* + * Copyright (c) 2021, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ +#ifndef AOM_FLOW_ESTIMATION_CORNER_MATCH_H_ +#define AOM_FLOW_ESTIMATION_CORNER_MATCH_H_ + +#include <stdio.h> +#include <stdlib.h> +#include <stdbool.h> +#include <memory.h> + +#include "aom_dsp/flow_estimation/flow_estimation.h" +#include "aom_scale/yv12config.h" + +#ifdef __cplusplus +extern "C" { +#endif + +#define MATCH_SZ 13 +#define MATCH_SZ_BY2 ((MATCH_SZ - 1) / 2) +#define MATCH_SZ_SQ (MATCH_SZ * MATCH_SZ) + +CorrespondenceList *aom_compute_corner_match(YV12_BUFFER_CONFIG *src, + YV12_BUFFER_CONFIG *ref, + int bit_depth); + +bool aom_fit_global_model_to_correspondences(CorrespondenceList *corrs, + TransformationType type, + MotionModel *params_by_motion, + int num_motions); + +bool aom_fit_local_model_to_correspondences(CorrespondenceList *corrs, + PixelRect *rect, + TransformationType type, + double *mat); + +void aom_free_correspondence_list(CorrespondenceList *list); + +#ifdef __cplusplus +} +#endif + +#endif // AOM_FLOW_ESTIMATION_CORNER_MATCH_H_
diff --git a/aom_dsp/flow_estimation/disflow.c b/aom_dsp/flow_estimation/disflow.c new file mode 100644 index 0000000..3e76b4f --- /dev/null +++ b/aom_dsp/flow_estimation/disflow.c
@@ -0,0 +1,497 @@ +/* + * Copyright (c) 2022, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#include "aom_dsp/aom_dsp_common.h" +#include "aom_dsp/flow_estimation/disflow.h" +#include "aom_dsp/flow_estimation/corner_detect.h" +#include "aom_dsp/flow_estimation/pyramid.h" +#include "aom_dsp/flow_estimation/ransac.h" +#include "aom_mem/aom_mem.h" + +#include "config/av1_rtcd.h" + +// TODO(rachelbarker): Move needed code from av1/ to aom_dsp/ +#include "av1/common/resize.h" + +#include <assert.h> + +// Size of square patches in the disflow dense grid +#define PATCH_SIZE 8 +// Center point of square patch +#define PATCH_CENTER ((PATCH_SIZE + 1) >> 1) +// Step size between patches, lower value means greater patch overlap +#define PATCH_STEP 1 +// Warp error convergence threshold for disflow +#define DISFLOW_ERROR_TR 0.01 +// Max number of iterations if warp convergence is not found +// TODO(rachelbarker): Experiment with different numbers of pyramid levels +// and numbers of refinement steps per level +#define DISFLOW_MAX_ITR 10 + +// Don't use points around the frame border since they are less reliable +static INLINE int valid_point(int x, int y, int width, int height) { + return (x > (PATCH_SIZE + PATCH_CENTER)) && + (x < (width - PATCH_SIZE - PATCH_CENTER)) && + (y > (PATCH_SIZE + PATCH_CENTER)) && + (y < (height - PATCH_SIZE - PATCH_CENTER)); +} + +static int determine_disflow_correspondence(int *frm_corners, + int num_frm_corners, double *flow_u, + double *flow_v, int width, + int height, int stride, + Correspondence *correspondences) { + int num_correspondences = 0; + int x, y; + for (int i = 0; i < num_frm_corners; ++i) { + x = frm_corners[2 * i]; + y = frm_corners[2 * i + 1]; + if (valid_point(x, y, width, height)) { + correspondences[num_correspondences].x = x; + correspondences[num_correspondences].y = y; + correspondences[num_correspondences].rx = x + flow_u[y * stride + x]; + correspondences[num_correspondences].ry = y + flow_v[y * stride + x]; + num_correspondences++; + } + } + return num_correspondences; +} + +static void getCubicKernel(double x, double *kernel) { + assert(0 <= x && x < 1); + double x2 = x * x; + double x3 = x2 * x; + kernel[0] = -0.5 * x + x2 - 0.5 * x3; + kernel[1] = 1.0 - 2.5 * x2 + 1.5 * x3; + kernel[2] = 0.5 * x + 2.0 * x2 - 1.5 * x3; + kernel[3] = -0.5 * x2 + 0.5 * x3; +} + +static double getCubicValue(double *p, double *kernel) { + return kernel[0] * p[0] + kernel[1] * p[1] + kernel[2] * p[2] + + kernel[3] * p[3]; +} + +// Warps a block using flow vector [u, v] and computes the mse +static double compute_warp_and_error(unsigned char *ref, unsigned char *frm, + int width, int height, int stride, int x, + int y, double u, double v, int16_t *dt) { + unsigned char warped; + int x_w, y_w; + double mse = 0; + int16_t err = 0; + + // Split offset into integer and fractional parts, and compute cubic + // interpolation kernels + int u_int = (int)floor(u); + int v_int = (int)floor(v); + double u_frac = u - u_int; + double v_frac = v - v_int; + + double h_kernel[4]; + double v_kernel[4]; + getCubicKernel(u_frac, h_kernel); + getCubicKernel(v_frac, v_kernel); + + // Storage for intermediate values between the two convolution directions + double tmp_[PATCH_SIZE * (PATCH_SIZE + 3)]; + double *tmp = tmp_ + PATCH_SIZE; // Offset by one row + + // Clamp coordinates so that all pixels we fetch will remain within the + // allocated border region, but allow them to go far enough out that + // the border pixels' values do not change. + // Since we are calculating an 8x8 block, the bottom-right pixel + // in the block has coordinates (x0 + 7, y0 + 7). Then, the cubic + // interpolation has 4 taps, meaning that the output of pixel + // (x_w, y_w) depends on the pixels in the range + // ([x_w - 1, x_w + 2], [y_w - 1, y_w + 2]). + // + // Thus the most extreme coordinates which will be fetched are + // (x0 - 1, y0 - 1) and (x0 + 9, y0 + 9). + int x0 = clamp(x + u_int, -9, width); + int y0 = clamp(y + v_int, -9, height); + + // Horizontal convolution + for (int i = -1; i < PATCH_SIZE + 2; ++i) { + y_w = y0 + i; + for (int j = 0; j < PATCH_SIZE; ++j) { + x_w = x0 + j; + double arr[4]; + + arr[0] = (double)ref[y_w * stride + (x_w - 1)]; + arr[1] = (double)ref[y_w * stride + (x_w + 0)]; + arr[2] = (double)ref[y_w * stride + (x_w + 1)]; + arr[3] = (double)ref[y_w * stride + (x_w + 2)]; + + tmp[i * PATCH_SIZE + j] = getCubicValue(arr, h_kernel); + } + } + + // Vertical convolution + for (int i = 0; i < PATCH_SIZE; ++i) { + for (int j = 0; j < PATCH_SIZE; ++j) { + double *p = &tmp[i * PATCH_SIZE + j]; + double arr[4] = { p[-PATCH_SIZE], p[0], p[PATCH_SIZE], + p[2 * PATCH_SIZE] }; + double result = getCubicValue(arr, v_kernel); + + warped = clamp((int)(result + 0.5), 0, 255); + err = warped - frm[(x + j) + (y + i) * stride]; + mse += err * err; + dt[i * PATCH_SIZE + j] = err; + } + } + + mse /= (PATCH_SIZE * PATCH_SIZE); + return mse; +} + +// Computes the components of the system of equations used to solve for +// a flow vector. This includes: +// 1.) The hessian matrix for optical flow. This matrix is in the +// form of: +// +// M = |sum(dx * dx) sum(dx * dy)| +// |sum(dx * dy) sum(dy * dy)| +// +// 2.) b = |sum(dx * dt)| +// |sum(dy * dt)| +// Where the sums are computed over a square window of PATCH_SIZE. +static INLINE void compute_flow_system(const double *dx, int dx_stride, + const double *dy, int dy_stride, + const int16_t *dt, int dt_stride, + double *M, double *b) { + for (int i = 0; i < PATCH_SIZE; i++) { + for (int j = 0; j < PATCH_SIZE; j++) { + M[0] += dx[i * dx_stride + j] * dx[i * dx_stride + j]; + M[1] += dx[i * dx_stride + j] * dy[i * dy_stride + j]; + M[3] += dy[i * dy_stride + j] * dy[i * dy_stride + j]; + + b[0] += dx[i * dx_stride + j] * dt[i * dt_stride + j]; + b[1] += dy[i * dy_stride + j] * dt[i * dt_stride + j]; + } + } + + M[2] = M[1]; +} + +// Solves a general Mx = b where M is a 2x2 matrix and b is a 2x1 matrix +static INLINE void solve_2x2_system(const double *M, const double *b, + double *output_vec) { + double M_0 = M[0]; + double M_3 = M[3]; + double det = (M_0 * M_3) - (M[1] * M[2]); + if (det < 1e-5) { + // Handle singular matrix + // TODO(sarahparker) compare results using pseudo inverse instead + M_0 += 1e-10; + M_3 += 1e-10; + det = (M_0 * M_3) - (M[1] * M[2]); + } + const double det_inv = 1 / det; + const double mult_b0 = det_inv * b[0]; + const double mult_b1 = det_inv * b[1]; + output_vec[0] = M_3 * mult_b0 - M[1] * mult_b1; + output_vec[1] = -M[2] * mult_b0 + M_0 * mult_b1; +} + +/* +static INLINE void image_difference(const uint8_t *src, int src_stride, + const uint8_t *ref, int ref_stride, + int16_t *dst, int dst_stride, int height, + int width) { + const int block_unit = 8; + // Take difference in 8x8 blocks to make use of optimized diff function + for (int i = 0; i < height; i += block_unit) { + for (int j = 0; j < width; j += block_unit) { + aom_subtract_block(block_unit, block_unit, dst + i * dst_stride + j, + dst_stride, src + i * src_stride + j, src_stride, + ref + i * ref_stride + j, ref_stride); + } + } +} +*/ + +static INLINE void compute_flow_at_point(unsigned char *frm, unsigned char *ref, + int x, int y, int width, int height, + int stride, double *u, double *v) { + double M[4] = { 0 }; + double b[2] = { 0 }; + double tmp_output_vec[2] = { 0 }; + double error = 0; + int16_t dt[PATCH_SIZE * PATCH_SIZE]; + double o_u = *u; + double o_v = *v; + + double dx_tmp[PATCH_SIZE * PATCH_SIZE]; + double dy_tmp[PATCH_SIZE * PATCH_SIZE]; + + // Compute gradients within this patch + unsigned char *frm_patch = &frm[y * stride + x]; + av1_convolve_2d_sobel_y_c(frm_patch, stride, dx_tmp, PATCH_SIZE, PATCH_SIZE, + PATCH_SIZE, 1, 1.0); + av1_convolve_2d_sobel_y_c(frm_patch, stride, dy_tmp, PATCH_SIZE, PATCH_SIZE, + PATCH_SIZE, 0, 1.0); + + for (int itr = 0; itr < DISFLOW_MAX_ITR; itr++) { + error = compute_warp_and_error(ref, frm, width, height, stride, x, y, *u, + *v, dt); + if (error <= DISFLOW_ERROR_TR) break; + compute_flow_system(dx_tmp, PATCH_SIZE, dy_tmp, PATCH_SIZE, dt, PATCH_SIZE, + M, b); + solve_2x2_system(M, b, tmp_output_vec); + *u += tmp_output_vec[0]; + *v += tmp_output_vec[1]; + } + if (fabs(*u - o_u) > PATCH_SIZE || fabs(*v - o_u) > PATCH_SIZE) { + *u = o_u; + *v = o_v; + } +} + +static void fill_flow_field_borders(double *flow, int width, int height, + int stride) { + // Calculate the bounds of the rectangle which was filled in by + // compute_flow_field() before calling this function. + // These indices are inclusive on both ends. + const int left_index = PATCH_CENTER; + const int right_index = (width - PATCH_SIZE - 1) + PATCH_CENTER; + const int top_index = PATCH_CENTER; + const int bottom_index = (height - PATCH_SIZE - 1) + PATCH_CENTER; + + // Left area + for (int i = top_index; i <= bottom_index; i += PATCH_STEP) { + double *row = flow + i * stride; + double left = row[left_index]; + for (int j = 0; j < left_index; j++) { + row[j] = left; + } + } + + // Right area + for (int i = top_index; i <= bottom_index; i += PATCH_STEP) { + double *row = flow + i * stride; + double right = row[right_index]; + for (int j = right_index + 1; j < width; j++) { + row[j] = right; + } + } + + // Top area + double *top_row = flow + top_index * stride; + for (int i = 0; i < top_index; i++) { + double *row = flow + i * stride; + memcpy(row, top_row, width * sizeof(double)); + } + + // Bottom area + double *bottom_row = flow + bottom_index * stride; + for (int i = bottom_index + 1; i < height; i++) { + double *row = flow + i * stride; + memcpy(row, bottom_row, width * sizeof(double)); + } +} + +// make sure flow_u and flow_v start at 0 +static void compute_flow_field(ImagePyramid *frm_pyr, ImagePyramid *ref_pyr, + double *flow_u, double *flow_v) { + int cur_width, cur_height, cur_stride, cur_loc, patch_loc, patch_center; + double *u_upscale = + aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_u)); + double *v_upscale = + aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_v)); + + assert(frm_pyr->n_levels == ref_pyr->n_levels); + + // Compute flow field from coarsest to finest level of the pyramid +#if PATCH_STEP != 1 + // TODO(rachelbarker): This function, as written, only works if PATCH_STEP + // == 1. For any other value, the border filling and interpolation code will + // need to be reworked to handle the fact that there will be rows and columns + // with no flow data +#error "compute_flow_field() needs updating for PATCH_STEP != 1" +#endif + + for (int level = frm_pyr->n_levels - 1; level >= 0; --level) { + cur_width = frm_pyr->widths[level]; + cur_height = frm_pyr->heights[level]; + cur_stride = frm_pyr->strides[level]; + cur_loc = frm_pyr->level_loc[level]; + + for (int i = 0; i < cur_height - PATCH_SIZE; i += PATCH_STEP) { + for (int j = 0; j < cur_width - PATCH_SIZE; j += PATCH_STEP) { + patch_loc = i * cur_stride + j; + patch_center = patch_loc + PATCH_CENTER * cur_stride + PATCH_CENTER; + compute_flow_at_point(frm_pyr->level_buffer + cur_loc, + ref_pyr->level_buffer + cur_loc, j, i, cur_width, + cur_height, cur_stride, flow_u + patch_center, + flow_v + patch_center); + } + } + + // Fill in the areas which we haven't explicitly computed, with copies + // of the outermost values which we did compute + fill_flow_field_borders(flow_u, cur_width, cur_height, cur_stride); + fill_flow_field_borders(flow_v, cur_width, cur_height, cur_stride); + + if (level > 0) { + int h_upscale = frm_pyr->heights[level - 1]; + int w_upscale = frm_pyr->widths[level - 1]; + int s_upscale = frm_pyr->strides[level - 1]; + av1_upscale_plane_double_prec(flow_u, cur_height, cur_width, cur_stride, + u_upscale, h_upscale, w_upscale, s_upscale); + av1_upscale_plane_double_prec(flow_v, cur_height, cur_width, cur_stride, + v_upscale, h_upscale, w_upscale, s_upscale); + + // Multiply all flow vectors by 2. + // When we move down a pyramid level, the image resolution doubles. + // Thus we need to double all vectors in order for them to represent + // the same translation at the next level down + for (int i = 0; i < h_upscale; i++) { + for (int j = 0; j < w_upscale; j++) { + int index = i * s_upscale + j; + flow_u[index] = u_upscale[index] * 2.0; + flow_v[index] = v_upscale[index] * 2.0; + } + } + } + } + aom_free(u_upscale); + aom_free(v_upscale); +} + +FlowField *aom_alloc_flow_field(int width, int height, int stride) { + FlowField *flow = (FlowField *)aom_malloc(sizeof(FlowField)); + if (flow == NULL) return NULL; + + flow->width = width; + flow->height = height; + flow->stride = stride; + + size_t flow_size = stride * (size_t)height; + flow->u = aom_calloc(flow_size, sizeof(double)); + flow->v = aom_calloc(flow_size, sizeof(double)); + + if (flow->u == NULL || flow->v == NULL) { + aom_free(flow->u); + aom_free(flow->v); + aom_free(flow); + return NULL; + } + + return flow; +} + +void aom_free_flow_field(FlowField *flow) { + aom_free(flow->u); + aom_free(flow->v); + aom_free(flow); +} + +FlowField *aom_compute_flow_field(YV12_BUFFER_CONFIG *frm, + YV12_BUFFER_CONFIG *ref, int bit_depth) { + const int frm_width = frm->y_width; + const int frm_height = frm->y_height; + assert(frm->y_width == ref->y_width); + assert(frm->y_height == ref->y_height); + + // Compute pyramids if necessary. + // These are cached alongside the framebuffer to avoid unnecessary + // recomputation. When the framebuffer is freed, or reused for a new frame, + // these pyramids will be automatically freed. + if (!frm->y_pyramid) { + frm->y_pyramid = aom_compute_pyramid(frm, bit_depth, MAX_PYRAMID_LEVELS); + assert(frm->y_pyramid); + } + if (!ref->y_pyramid) { + ref->y_pyramid = aom_compute_pyramid(ref, bit_depth, MAX_PYRAMID_LEVELS); + assert(ref->y_pyramid); + } + + ImagePyramid *frm_pyr = frm->y_pyramid; + ImagePyramid *ref_pyr = ref->y_pyramid; + + FlowField *flow = + aom_alloc_flow_field(frm_width, frm_height, frm_pyr->strides[0]); + + compute_flow_field(frm_pyr, ref_pyr, flow->u, flow->v); + + return flow; +} + +bool aom_fit_global_model_to_flow_field(FlowField *flow, + TransformationType type, + YV12_BUFFER_CONFIG *frm, int bit_depth, + MotionModel *params_by_motion, + int num_motions) { + int num_correspondences; + + if (!frm->corners) { + aom_find_corners_in_frame(frm, bit_depth); + } + + // find correspondences between the two images using the flow field + Correspondence *correspondences = + aom_malloc(frm->num_corners * sizeof(*correspondences)); + num_correspondences = determine_disflow_correspondence( + frm->corners, frm->num_corners, flow->u, flow->v, flow->width, + flow->height, flow->stride, correspondences); + ransac(correspondences, num_correspondences, type, params_by_motion, + num_motions); + + aom_free(correspondences); + // Set num_inliers = 0 for motions with too few inliers so they are ignored. + for (int i = 0; i < num_motions; ++i) { + if (params_by_motion[i].num_inliers < + MIN_INLIER_PROB * num_correspondences) { + params_by_motion[i].num_inliers = 0; + } + } + + // Return true if any one of the motions has inliers. + for (int i = 0; i < num_motions; ++i) { + if (params_by_motion[i].num_inliers > 0) return true; + } + return false; +} + +bool aom_fit_local_model_to_flow_field(const FlowField *flow, + const PixelRect *rect, + TransformationType type, double *mat) { + int width = rect_height(rect); + int height = rect_width(rect); + int num_points = width * height; + + // TODO(rachelbarker): Downsample if num_points is > some threshold? + double *pts1 = aom_malloc(num_points * 2 * sizeof(double)); + double *pts2 = aom_malloc(num_points * 2 * sizeof(double)); + int index = 0; + + for (int y = rect->top; y < rect->bottom; y++) { + for (int x = rect->left; x < rect->right; x++) { + pts1[2 * index + 0] = (double)x; + pts1[2 * index + 1] = (double)y; + pts2[2 * index + 0] = (double)x + flow->u[y * flow->stride + x]; + pts2[2 * index + 1] = (double)y + flow->v[y * flow->stride + x]; + index++; + } + } + + // Check that we filled the expected number of points + assert(index == num_points); + + bool result = aom_fit_motion_model(type, num_points, pts1, pts2, mat); + + aom_free(pts1); + aom_free(pts2); + return result; +}
diff --git a/aom_dsp/flow_estimation/disflow.h b/aom_dsp/flow_estimation/disflow.h new file mode 100644 index 0000000..5ebece5 --- /dev/null +++ b/aom_dsp/flow_estimation/disflow.h
@@ -0,0 +1,48 @@ +/* + * Copyright (c) 2022, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ +#ifndef AOM_FLOW_ESTIMATION_DISFLOW_BASED_H_ +#define AOM_FLOW_ESTIMATION_DISFLOW_BASED_H_ + +#include <stdbool.h> + +#include "aom_dsp/flow_estimation/flow_estimation.h" +#include "aom_dsp/rect.h" +#include "aom_scale/yv12config.h" + +#ifdef __cplusplus +extern "C" { +#endif + +// Number of pyramid levels in disflow computation +#define DISFLOW_PYRAMID_LEVELS 2 + +FlowField *aom_alloc_flow_field(int width, int height, int stride); +void aom_free_flow_field(FlowField *flow); + +FlowField *aom_compute_flow_field(YV12_BUFFER_CONFIG *frm, + YV12_BUFFER_CONFIG *ref, int bit_depth); + +bool aom_fit_global_model_to_flow_field(FlowField *flow, + TransformationType type, + YV12_BUFFER_CONFIG *frm, int bit_depth, + MotionModel *params_by_motion, + int num_motions); + +bool aom_fit_local_model_to_flow_field(const FlowField *flow, + const PixelRect *rect, + TransformationType type, double *mat); + +#ifdef __cplusplus +} +#endif + +#endif // AOM_FLOW_ESTIMATION_DISFLOW_BASED_H_
diff --git a/aom_dsp/flow_estimation/flow_estimation.c b/aom_dsp/flow_estimation/flow_estimation.c new file mode 100644 index 0000000..5b4b2e7 --- /dev/null +++ b/aom_dsp/flow_estimation/flow_estimation.c
@@ -0,0 +1,90 @@ +/* + * Copyright (c) 2022, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#include <assert.h> + +#include "aom_dsp/flow_estimation/corner_match.h" +#include "aom_dsp/flow_estimation/flow_estimation.h" +#include "aom_dsp/flow_estimation/disflow.h" +#include "aom_scale/yv12config.h" +#include "aom_mem/aom_mem.h" + +FlowData *aom_compute_flow_data(YV12_BUFFER_CONFIG *src, + YV12_BUFFER_CONFIG *ref, int bit_depth, + GlobalMotionEstimationType gm_estimation_type) { + FlowData *flow_data = aom_malloc(sizeof(*flow_data)); + if (!flow_data) { + return NULL; + } + + flow_data->method = gm_estimation_type; + + if (flow_data->method == GLOBAL_MOTION_FEATURE_BASED) { + flow_data->corrs = aom_compute_corner_match(src, ref, bit_depth); + } else if (flow_data->method == GLOBAL_MOTION_DISFLOW_BASED) { + flow_data->flow = aom_compute_flow_field(src, ref, bit_depth); + } else { + assert(0 && "Unknown global motion estimation type"); + aom_free(flow_data); + return NULL; + } + + return flow_data; +} + +// Fit one or several models of a given type to the specified flow data. +// This function fits models to the entire frame, using the RANSAC method +// to fit models in a noise-resilient way, and returns the list of inliers +// for each model found +bool aom_fit_global_motion_model(FlowData *flow_data, TransformationType type, + YV12_BUFFER_CONFIG *src, int bit_depth, + MotionModel *params_by_motion, + int num_motions) { + if (flow_data->method == GLOBAL_MOTION_FEATURE_BASED) { + return aom_fit_global_model_to_correspondences( + flow_data->corrs, type, params_by_motion, num_motions); + } else if (flow_data->method == GLOBAL_MOTION_DISFLOW_BASED) { + return aom_fit_global_model_to_flow_field( + flow_data->flow, type, src, bit_depth, params_by_motion, num_motions); + } else { + assert(0 && "Unknown global motion estimation type"); + return 0; + } +} + +// Fit a model of a given type to a subset of the specified flow data. +// This does not used the RANSAC method, so is more noise-sensitive than +// aom_fit_global_motion_model(), but in the context of fitting models +// to single blocks this is not an issue. +bool aom_fit_local_motion_model(FlowData *flow_data, PixelRect *rect, + TransformationType type, double *mat) { + if (flow_data->method == GLOBAL_MOTION_FEATURE_BASED) { + return aom_fit_local_model_to_correspondences(flow_data->corrs, rect, type, + mat); + } else if (flow_data->method == GLOBAL_MOTION_DISFLOW_BASED) { + return aom_fit_local_model_to_flow_field(flow_data->flow, rect, type, mat); + } else { + assert(0 && "Unknown global motion estimation type"); + return 0; + } +} + +void aom_free_flow_data(FlowData *flow_data) { + if (flow_data->method == GLOBAL_MOTION_FEATURE_BASED) { + aom_free_correspondence_list(flow_data->corrs); + } else if (flow_data->method == GLOBAL_MOTION_DISFLOW_BASED) { + aom_free_flow_field(flow_data->flow); + } else { + assert(0 && "Unknown global motion estimation type"); + } + aom_free(flow_data); +}
diff --git a/aom_dsp/flow_estimation/flow_estimation.h b/aom_dsp/flow_estimation/flow_estimation.h new file mode 100644 index 0000000..8876615 --- /dev/null +++ b/aom_dsp/flow_estimation/flow_estimation.h
@@ -0,0 +1,136 @@ +/* + * Copyright (c) 2022, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ +#ifndef AOM_FLOW_ESTIMATION_H_ +#define AOM_FLOW_ESTIMATION_H_ + +#include "aom_dsp/rect.h" +#include "aom_ports/mem.h" +#include "aom_scale/yv12config.h" + +#ifdef __cplusplus +extern "C" { +#endif + +#define MAX_PARAMDIM 9 +#define MAX_CORNERS 4096 +#define MIN_INLIER_PROB 0.1 + +/* clang-format off */ +enum { + IDENTITY = 0, // identity transformation, 0-parameter + TRANSLATION = 1, // translational motion 2-parameter + ROTATION = 2, // rotation about some point, 3-parameter + ZOOM = 3, // zoom in/out on some point, 3-parameter + VERTSHEAR = 4, // translation + vertical shear, 3-parameter + HORZSHEAR = 5, // translation + horizontal shear, 3-parameter + UZOOM = 6, // unequal zoom, 4-parameter + ROTZOOM = 7, // equal zoom, then rotate, 4-parameter + ROTUZOOM = 8, // unequal zoom, then rotate, 5-parameter + AFFINE = 9, // general affine, 6-parameter + VERTRAPEZOID = 10, // vertical-only perspective, 6-parameter + HORTRAPEZOID = 11, // horizontal-only perspective, 6-parameter + HOMOGRAPHY = 12, // general perspective transformation, 8-parameter + TRANS_TYPES, +} UENUM1BYTE(TransformationType); +/* clang-format on */ + +// number of parameters used by each transformation in TransformationTypes +static const int trans_model_params[TRANS_TYPES] = { 0, 2, 3, 3, 3, 3, 4, + 4, 5, 6, 6, 6, 8 }; + +typedef enum { + GLOBAL_MOTION_FEATURE_BASED, + GLOBAL_MOTION_DISFLOW_BASED, +} GlobalMotionEstimationType; + +typedef struct { + double params[MAX_PARAMDIM - 1]; + int *inliers; + int num_inliers; +} MotionModel; + +typedef struct { + double x, y; + double rx, ry; +} Correspondence; + +typedef struct { + int num_correspondences; + Correspondence *correspondences; +} CorrespondenceList; + +typedef struct { + // x and y directions of flow, per patch + double *u; + double *v; + + // Sizes of the above arrays + size_t width; + size_t height; + size_t stride; +} FlowField; + +// We want to present external code with a generic type, which holds whatever +// data is needed for the desired motion estimation method. +// As different methods use different data, we store this in a tagged union, +// with the selected motion estimation type as the tag. +typedef struct { + GlobalMotionEstimationType method; + union { + CorrespondenceList *corrs; + FlowField *flow; + }; +} FlowData; + +FlowData *aom_compute_flow_data(YV12_BUFFER_CONFIG *src, + YV12_BUFFER_CONFIG *ref, int bit_depth, + GlobalMotionEstimationType gm_estimation_type); + +/* + Computes "num_motions" candidate global motion parameters between two frames. + The array "params_by_motion" should be length 8 * "num_motions". The ordering + of each set of parameters is best described by the homography: + + [x' (m2 m3 m0 [x + z . y' = m4 m5 m1 * y + 1] m6 m7 1) 1] + + where m{i} represents the ith value in any given set of parameters. + + "num_inliers" should be length "num_motions", and will be populated with the + number of inlier feature points for each motion. Params for which the + num_inliers entry is 0 should be ignored by the caller. +*/ +// Fit one or several models of a given type to the specified flow data. +// This function fits models to the entire frame, using the RANSAC method +// to fit models in a noise-resilient way, and returns the list of inliers +// for each model found +// TODO: Cleanup this comment +bool aom_fit_global_motion_model(FlowData *flow_data, TransformationType type, + YV12_BUFFER_CONFIG *src, int bit_depth, + MotionModel *params_by_motion, + int num_motions); + +// Fit a model of a given type to a subset of the specified flow data. +// This does not used the RANSAC method, so is more noise-sensitive than +// aom_fit_global_motion_model(), but in the context of fitting models +// to single blocks this is not an issue. +bool aom_fit_local_motion_model(FlowData *flow_data, PixelRect *rect, + TransformationType type, double *mat); + +void aom_free_flow_data(FlowData *flow_data); + +#ifdef __cplusplus +} +#endif + +#endif // AOM_FLOW_ESTIMATION_H_
diff --git a/aom_dsp/flow_estimation/pyramid.c b/aom_dsp/flow_estimation/pyramid.c new file mode 100644 index 0000000..fd35614 --- /dev/null +++ b/aom_dsp/flow_estimation/pyramid.c
@@ -0,0 +1,174 @@ +/* + * Copyright (c) 2022, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#include "aom_dsp/flow_estimation/pyramid.h" +#include "aom_mem/aom_mem.h" +#include "aom_ports/bitops.h" + +// TODO(rachelbarker): Move needed code from av1/ to aom_dsp/ +#include "av1/common/resize.h" + +#include <assert.h> +#include <string.h> + +// TODO(rachelbarker): Check for allocations returning NULL +// TODO(rachelbarker): Align the first image pixel of each level to some +// convenient power of two (eg, to 16 bytes for SIMD) +static INLINE ImagePyramid *alloc_pyramid(int width, int height, int n_levels) { + assert(n_levels <= MAX_PYRAMID_LEVELS); + + // Limit number of levels on small frames + const int msb = get_msb(AOMMIN(width, height)); + const int max_levels = AOMMAX(msb - MIN_PYRAMID_SIZE_LOG2, 1); + n_levels = AOMMIN(n_levels, max_levels); + + ImagePyramid *pyr = aom_malloc(sizeof(*pyr)); + pyr->n_levels = n_levels; + + // Compute sizes and offsets for each pyramid level + size_t buffer_size = 0; + + for (int level = 0; level < n_levels; level++) { + int level_width = width >> level; + int level_height = height >> level; + int level_alloc_width = level_width + 2 * PYRAMID_PADDING; + int level_alloc_height = level_height + 2 * PYRAMID_PADDING; + + pyr->widths[level] = level_width; + pyr->heights[level] = level_height; + pyr->strides[level] = level_alloc_width; + + // Offset the level_loc table so that each element points to the first image + // pixel, not the first padding pixel + size_t level_alloc_start = buffer_size; + pyr->level_loc[level] = level_alloc_start + + PYRAMID_PADDING * level_alloc_width + + PYRAMID_PADDING; + + buffer_size += level_alloc_width * level_alloc_height; + } + + // TODO(rachelbarker): Do we need to zero this buffer? + pyr->level_buffer = aom_malloc(buffer_size * sizeof(*pyr->level_buffer)); + + return pyr; +} + +// Fill the border region of a pyramid frame. +// This must be called after the main image area is filled out. +// `img_buf` should point to the first pixel in the image area, +// ie. it should be pyr->level_buffer + pyr->level_loc[level]. +static INLINE void fill_border(unsigned char *img_buf, const int width, + const int height, const int stride) { + // Fill left and right areas + for (int row = 0; row < height; row++) { + unsigned char *row_start = &img_buf[row * stride]; + unsigned char left_pixel = row_start[0]; + memset(row_start - PYRAMID_PADDING, left_pixel, PYRAMID_PADDING); + unsigned char right_pixel = row_start[width - 1]; + memset(row_start + width, right_pixel, PYRAMID_PADDING); + } + + // Fill top area + for (int row = -PYRAMID_PADDING; row < 0; row++) { + unsigned char *row_start = &img_buf[row * stride]; + memcpy(row_start - PYRAMID_PADDING, img_buf - PYRAMID_PADDING, + width + 2 * PYRAMID_PADDING); + } + + // Fill bottom area + unsigned char *last_row_start = &img_buf[(height - 1) * stride]; + for (int row = height; row < height + PYRAMID_PADDING; row++) { + unsigned char *row_start = &img_buf[row * stride]; + memcpy(row_start - PYRAMID_PADDING, last_row_start - PYRAMID_PADDING, + width + 2 * PYRAMID_PADDING); + } +} + +// Compute coarse to fine pyramids for a frame +static INLINE void fill_pyramid(YV12_BUFFER_CONFIG *frm, int bit_depth, + ImagePyramid *frm_pyr) { + int n_levels = frm_pyr->n_levels; + const int frm_width = frm->y_width; + const int frm_height = frm->y_height; + const int frm_stride = frm->y_stride; + assert((frm_width >> n_levels) >= 0); + assert((frm_height >> n_levels) >= 0); + + int cur_width, cur_height, cur_stride, cur_loc; + cur_width = frm_pyr->widths[0]; + cur_height = frm_pyr->heights[0]; + cur_stride = frm_pyr->strides[0]; + cur_loc = frm_pyr->level_loc[0]; + + assert(frm_width == frm_pyr->widths[0]); + assert(frm_height == frm_pyr->heights[0]); + + // Fill out the initial pyramid level + if (frm->flags & YV12_FLAG_HIGHBITDEPTH) { + // For frames stored in 16-bit buffers, we need to downconvert to 8 bits + uint16_t *frm_buffer = CONVERT_TO_SHORTPTR(frm->y_buffer); + uint8_t *pyr_buffer = frm_pyr->level_buffer + cur_loc; + for (int y = 0; y < frm_height; y++) { + uint16_t *frm_row = frm_buffer + y * frm_stride; + uint8_t *pyr_row = pyr_buffer + y * cur_stride; + for (int x = 0; x < frm_width; x++) { + pyr_row[x] = frm_row[x] >> (bit_depth - 8); + } + } + } else { + // For frames stored in 8-bit buffers, we can simply copy the frame data + uint8_t *frm_buffer = frm->y_buffer; + uint8_t *pyr_buffer = frm_pyr->level_buffer + cur_loc; + for (int y = 0; y < frm_height; y++) { + uint8_t *frm_row = frm_buffer + y * frm_stride; + uint8_t *pyr_row = pyr_buffer + y * cur_stride; + memcpy(pyr_row, frm_row, frm_width); + } + } + + fill_border(frm_pyr->level_buffer + cur_loc, cur_width, cur_height, + cur_stride); + + // Start at the finest level and resize down to the coarsest level + for (int level = 1; level < n_levels; ++level) { + cur_width = frm_pyr->widths[level]; + cur_height = frm_pyr->heights[level]; + cur_stride = frm_pyr->strides[level]; + cur_loc = frm_pyr->level_loc[level]; + + // TODO(rachelbarker): Implement a special downsample-by-2 function + // to make this more efficient + av1_resize_plane(frm_pyr->level_buffer + frm_pyr->level_loc[level - 1], + frm_pyr->heights[level - 1], frm_pyr->widths[level - 1], + frm_pyr->strides[level - 1], + frm_pyr->level_buffer + cur_loc, cur_height, cur_width, + cur_stride); + fill_border(frm_pyr->level_buffer + cur_loc, cur_width, cur_height, + cur_stride); + } +} + +// Allocate and fill out a pyramid structure for a given frame +ImagePyramid *aom_compute_pyramid(YV12_BUFFER_CONFIG *frm, int bit_depth, + int n_levels) { + ImagePyramid *frm_pyr = alloc_pyramid(frm->y_width, frm->y_height, n_levels); + fill_pyramid(frm, bit_depth, frm_pyr); + return frm_pyr; +} + +void aom_free_pyramid(ImagePyramid *pyr) { + if (pyr) { + aom_free(pyr->level_buffer); + aom_free(pyr); + } +}
diff --git a/aom_dsp/flow_estimation/pyramid.h b/aom_dsp/flow_estimation/pyramid.h new file mode 100644 index 0000000..2082ab9 --- /dev/null +++ b/aom_dsp/flow_estimation/pyramid.h
@@ -0,0 +1,83 @@ +/* + * Copyright (c) 2022, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#ifndef AOM_FLOW_ESTIMATION_PYRAMID_H_ +#define AOM_FLOW_ESTIMATION_PYRAMID_H_ + +#ifdef __cplusplus +extern "C" { +#endif + +#include "config/aom_config.h" + +// Maximum number of pyramid levels +#if CONFIG_GM_USE_DISFLOW +// Disflow requires two pyramid levels +#define MAX_PYRAMID_LEVELS 2 +#else +// Feature based code only requires one pyramid level +#define MAX_PYRAMID_LEVELS 1 +#endif // CONFIG_GM_USE_DISFLOW + +// Minimum dimensions of a downsampled image +#define MIN_PYRAMID_SIZE_LOG2 3 +#define MIN_PYRAMID_SIZE (1 << MIN_PYRAMID_SIZE_LOG2) + +// Size of border around each pyramid image, in pixels +// Similarly to the border around regular image buffers, this border is filled +// with copies of the outermost pixels of the frame, to allow for more efficient +// convolution code +#define PYRAMID_PADDING 16 + +// Forward declare this struct rather than including "aom_scale/yv12config.h", +// so that that file can include this one without causing circular dependencies +struct yv12_buffer_config; + +// Struct for an image pyramid +typedef struct { + int n_levels; + int widths[MAX_PYRAMID_LEVELS]; + int heights[MAX_PYRAMID_LEVELS]; + int strides[MAX_PYRAMID_LEVELS]; + int level_loc[MAX_PYRAMID_LEVELS]; + unsigned char *level_buffer; +} ImagePyramid; + +// Allocate and fill out a downsampling pyramid for a given frame. +// +// The top level (index 0) will always be an 8-bit copy of the input frame, +// regardless of the input bit depth. Additional levels are then downscaled +// by powers of 2. +// +// Note on n_levels: +// * For feature-based global motion, n_levels need only be 1, +// which just constructs an 8-bit version of the input frame. +// * For disflow-based global motion, n_levels should equal +// DISFLOW_PYRAMID_LEVELS +// * In any case, n_levels must be <= MAX_PYRAMID_LEVELS +// +// For small input frames, the number of levels actually constructed +// will be limited so that the smallest image is at least MIN_PYRAMID_SIZE +// pixels along each side. +// +// However, if the input frame has a side of length < MIN_PYRAMID_SIZE, +// we will still construct the top level. +ImagePyramid *aom_compute_pyramid(struct yv12_buffer_config *frm, int bit_depth, + int n_levels); + +void aom_free_pyramid(ImagePyramid *pyr); + +#ifdef __cplusplus +} +#endif + +#endif // AOM_FLOW_ESTIMATION_PYRAMID_H_
diff --git a/aom_dsp/flow_estimation/ransac.c b/aom_dsp/flow_estimation/ransac.c new file mode 100644 index 0000000..6217899 --- /dev/null +++ b/aom_dsp/flow_estimation/ransac.c
@@ -0,0 +1,1773 @@ +/* + * Copyright (c) 2021, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ +#include <memory.h> +#include <math.h> +#include <time.h> +#include <stdio.h> +#include <stdbool.h> +#include <assert.h> + +#include "aom_dsp/linalg.h" +#include "aom_dsp/flow_estimation/ransac.h" +#include "aom_mem/aom_mem.h" + +// TODO(rachelbarker): Remove dependence on code in av1/encoder/ +#include "av1/encoder/random.h" + +#define MAX_MINPTS 4 +#define MAX_DEGENERATE_ITER 10 +#define MINPTS_MULTIPLIER 5 + +#define INLIER_THRESHOLD 1.25 +#define MIN_TRIALS 20 + +// Choose between three different algorithms for finding homographies. +// TODO(rachelbarker): Select one of these +// TODO(rachelbarker): See if these algorithms' stability can be improved +// by some kind of refinement method. eg, take the SVD result and do gradient +// descent from there +#define HORZTRAP_ALGORITHM 0 +#define VERTTRAP_ALGORITHM 0 +#define HOMOGRAPHY_ALGORITHM 0 + +//////////////////////////////////////////////////////////////////////////////// +// ransac +typedef bool (*IsDegenerateFunc)(double *p); +typedef bool (*FindTransformationFunc)(int points, double *points1, + double *points2, double *params); +typedef void (*ProjectPointsFunc)(double *mat, double *points, double *proj, + int n, int stride_points, int stride_proj); + +static void project_points_translation(double *mat, double *points, + double *proj, int n, int stride_points, + int stride_proj) { + int i; + for (i = 0; i < n; ++i) { + const double x = *(points++), y = *(points++); + *(proj++) = x + mat[0]; + *(proj++) = y + mat[1]; + points += stride_points - 2; + proj += stride_proj - 2; + } +} + +static void project_points_affine(double *mat, double *points, double *proj, + int n, int stride_points, int stride_proj) { + int i; + for (i = 0; i < n; ++i) { + const double x = *(points++), y = *(points++); + *(proj++) = mat[2] * x + mat[3] * y + mat[0]; + *(proj++) = mat[4] * x + mat[5] * y + mat[1]; + points += stride_points - 2; + proj += stride_proj - 2; + } +} + +static void project_points_homography(double *mat, double *points, double *proj, + const int n, const int stride_points, + const int stride_proj) { + int i; + double x, y, Z, Z_inv; + for (i = 0; i < n; ++i) { + x = *(points++), y = *(points++); + Z_inv = mat[6] * x + mat[7] * y + 1; + assert(fabs(Z_inv) > 0.000001); + Z = 1. / Z_inv; + *(proj++) = (mat[2] * x + mat[3] * y + mat[0]) * Z; + *(proj++) = (mat[4] * x + mat[5] * y + mat[1]) * Z; + points += stride_points - 2; + proj += stride_proj - 2; + } +} + +static void normalize_homography(double *pts, int n, double *T) { + double *p = pts; + double mean[2] = { 0, 0 }; + double msqe = 0; + double scale; + int i; + + assert(n > 0); + for (i = 0; i < n; ++i, p += 2) { + mean[0] += p[0]; + mean[1] += p[1]; + } + mean[0] /= n; + mean[1] /= n; + for (p = pts, i = 0; i < n; ++i, p += 2) { + p[0] -= mean[0]; + p[1] -= mean[1]; + msqe += sqrt(p[0] * p[0] + p[1] * p[1]); + } + msqe /= n; + scale = (msqe == 0 ? 1.0 : sqrt(2) / msqe); + T[0] = scale; + T[1] = 0; + T[2] = -scale * mean[0]; + T[3] = 0; + T[4] = scale; + T[5] = -scale * mean[1]; + T[6] = 0; + T[7] = 0; + T[8] = 1; + for (p = pts, i = 0; i < n; ++i, p += 2) { + p[0] *= scale; + p[1] *= scale; + } +} + +static void invnormalize_mat(double *T, double *iT) { + double is = 1.0 / T[0]; + double m0 = -T[2] * is; + double m1 = -T[5] * is; + iT[0] = is; + iT[1] = 0; + iT[2] = m0; + iT[3] = 0; + iT[4] = is; + iT[5] = m1; + iT[6] = 0; + iT[7] = 0; + iT[8] = 1; +} + +static void denormalize_homography(double *params, double *T1, double *T2) { + double iT2[9]; + double params2[9]; + invnormalize_mat(T2, iT2); + multiply_mat(params, T1, params2, 3, 3, 3); + multiply_mat(iT2, params2, params, 3, 3, 3); +} + +/* +static void denormalize_homography_reorder(double *params, double *T1, + double *T2) { + double params_denorm[MAX_PARAMDIM]; + memcpy(params_denorm, params, sizeof(*params) * 8); + params_denorm[8] = 1.0; + denormalize_homography(params_denorm, T1, T2); + params[0] = params_denorm[2]; + params[1] = params_denorm[5]; + params[2] = params_denorm[0]; + params[3] = params_denorm[1]; + params[4] = params_denorm[3]; + params[5] = params_denorm[4]; + params[6] = params_denorm[6]; + params[7] = params_denorm[7]; +} +*/ + +static void denormalize_affine_reorder(double *params, double *T1, double *T2) { + double params_denorm[MAX_PARAMDIM]; + params_denorm[0] = params[0]; + params_denorm[1] = params[1]; + params_denorm[2] = params[4]; + params_denorm[3] = params[2]; + params_denorm[4] = params[3]; + params_denorm[5] = params[5]; + params_denorm[6] = params_denorm[7] = 0; + params_denorm[8] = 1; + denormalize_homography(params_denorm, T1, T2); + params[0] = params_denorm[2]; + params[1] = params_denorm[5]; + params[2] = params_denorm[0]; + params[3] = params_denorm[1]; + params[4] = params_denorm[3]; + params[5] = params_denorm[4]; + params[6] = params[7] = 0; +} + +static void denormalize_rotzoom_reorder(double *params, double *T1, + double *T2) { + double params_denorm[MAX_PARAMDIM]; + params_denorm[0] = params[0]; + params_denorm[1] = params[1]; + params_denorm[2] = params[2]; + params_denorm[3] = -params[1]; + params_denorm[4] = params[0]; + params_denorm[5] = params[3]; + params_denorm[6] = params_denorm[7] = 0; + params_denorm[8] = 1; + denormalize_homography(params_denorm, T1, T2); + params[0] = params_denorm[2]; + params[1] = params_denorm[5]; + params[2] = params_denorm[0]; + params[3] = params_denorm[1]; + params[4] = -params[3]; + params[5] = params[2]; + params[6] = params[7] = 0; +} + +static void denormalize_translation_reorder(double *params, double *T1, + double *T2) { + double params_denorm[MAX_PARAMDIM]; + params_denorm[0] = 1; + params_denorm[1] = 0; + params_denorm[2] = params[0]; + params_denorm[3] = 0; + params_denorm[4] = 1; + params_denorm[5] = params[1]; + params_denorm[6] = params_denorm[7] = 0; + params_denorm[8] = 1; + denormalize_homography(params_denorm, T1, T2); + params[0] = params_denorm[2]; + params[1] = params_denorm[5]; + params[2] = params[5] = 1; + params[3] = params[4] = 0; + params[6] = params[7] = 0; +} + +/* +static void denormalize_zoom_reorder(double *params, double *T1, double *T2) { + double params_denorm[MAX_PARAMDIM]; + params_denorm[0] = params[0]; + params_denorm[1] = 0; + params_denorm[2] = params[1]; + params_denorm[3] = 0; + params_denorm[4] = params[0]; + params_denorm[5] = params[2]; + params_denorm[6] = params_denorm[7] = 0; + params_denorm[8] = 1; + denormalize_homography(params_denorm, T1, T2); + params[0] = params_denorm[2]; + params[1] = params_denorm[5]; + params[2] = params_denorm[0]; + params[3] = params_denorm[1]; + params[4] = -params[3]; + params[5] = params[2]; + params[6] = params[7] = 0; +} +*/ + +static double norm(double *x, int len) { + double normsq = 0.0; + for (int i = 0; i < len; ++i) normsq += x[i] * x[i]; + return sqrt(normsq); +} + +#if VERTTRAP_ALGORITHM == 0 +static bool find_vertrapezoid(int np, double *pts1, double *pts2, double *mat) { + // Implemented from Peter Kovesi's normalized implementation + const int nvar = 7; + const int np3 = np * 3; + double *a = (double *)aom_malloc(sizeof(*a) * np3 * nvar * 2); + double *U = a + np3 * nvar; + double S[7], V[7 * 7]; + int i, mini; + double sx, sy, dx, dy; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 3 * nvar + 0] = 0; + a[i * 3 * nvar + 1] = 0; + a[i * 3 * nvar + 2] = -sx; + a[i * 3 * nvar + 3] = -sy; + a[i * 3 * nvar + 4] = -1; + a[i * 3 * nvar + 5] = dy * sx; + a[i * 3 * nvar + 6] = dy; + + a[(i * 3 + 1) * nvar + 0] = sx; + a[(i * 3 + 1) * nvar + 1] = 1; + a[(i * 3 + 1) * nvar + 2] = 0; + a[(i * 3 + 1) * nvar + 3] = 0; + a[(i * 3 + 1) * nvar + 4] = 0; + a[(i * 3 + 1) * nvar + 5] = -dx * sx; + a[(i * 3 + 1) * nvar + 6] = -dx; + + a[(i * 3 + 2) * nvar + 0] = -dy * sx; + a[(i * 3 + 2) * nvar + 1] = -dy; + a[(i * 3 + 2) * nvar + 2] = dx * sx; + a[(i * 3 + 2) * nvar + 3] = dx * sy; + a[(i * 3 + 2) * nvar + 4] = dx; + a[(i * 3 + 2) * nvar + 5] = 0; + a[(i * 3 + 2) * nvar + 6] = 0; + } + + if (SVD(U, S, V, a, np3, nvar)) { + aom_free(a); + return false; + } else { + double minS = 1e12; + mini = -1; + for (i = 0; i < nvar; ++i) { + if (S[i] < minS) { + minS = S[i]; + mini = i; + } + } + } + double H[9]; + H[0] = V[0 * nvar + mini]; + H[1] = 0; + H[2] = V[1 * nvar + mini]; + H[3] = V[2 * nvar + mini]; + H[4] = V[3 * nvar + mini]; + H[5] = V[4 * nvar + mini]; + H[6] = V[5 * nvar + mini]; + H[7] = 0; + H[8] = V[6 * nvar + mini]; + // denormalize_homography_reorder(H, T1, T2); + aom_free(a); + if (H[8] == 0.0) { + return false; + } else { + // normalize + double f = 1.0 / H[8]; + // for (i = 0; i < 8; i++) mat[i] = f * H[i]; + mat[0] = f * H[2]; + mat[1] = f * H[5]; + mat[2] = f * H[0]; + mat[3] = f * H[1]; + mat[4] = f * H[3]; + mat[5] = f * H[4]; + mat[6] = f * H[6]; + mat[7] = f * H[7]; + } + return true; +} +#elif VERTTRAP_ALGORITHM == 1 +static bool find_vertrapezoid(int np, double *pts1, double *pts2, double *mat) { + // Implemented from Peter Kovesi's normalized implementation + const int nvar = 7; + const int np2 = np * 2; + double *a = (double *)aom_malloc(sizeof(*a) * np2 * nvar * 2); + double *U = a + np2 * nvar; + double S[7], V[7 * 7]; + int i, mini; + double sx, sy, dx, dy; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * nvar + 0] = 0; + a[i * 2 * nvar + 1] = 0; + a[i * 2 * nvar + 2] = -sx; + a[i * 2 * nvar + 3] = -sy; + a[i * 2 * nvar + 4] = -1; + a[i * 2 * nvar + 5] = dy * sx; + a[i * 2 * nvar + 6] = dy; + + a[(i * 2 + 1) * nvar + 0] = sx; + a[(i * 2 + 1) * nvar + 1] = 1; + a[(i * 2 + 1) * nvar + 2] = 0; + a[(i * 2 + 1) * nvar + 3] = 0; + a[(i * 2 + 1) * nvar + 4] = 0; + a[(i * 2 + 1) * nvar + 5] = -dx * sx; + a[(i * 2 + 1) * nvar + 6] = -dx; + } + + if (SVD(U, S, V, a, np2, nvar)) { + aom_free(a); + return false; + } else { + double minS = 1e12; + mini = -1; + for (i = 0; i < nvar; ++i) { + if (S[i] < minS) { + minS = S[i]; + mini = i; + } + } + } + double H[9]; + H[0] = V[0 * nvar + mini]; + H[1] = 0; + H[2] = V[1 * nvar + mini]; + H[3] = V[2 * nvar + mini]; + H[4] = V[3 * nvar + mini]; + H[5] = V[4 * nvar + mini]; + H[6] = V[5 * nvar + mini]; + H[7] = 0; + H[8] = V[6 * nvar + mini]; + // denormalize_homography_reorder(H, T1, T2); + aom_free(a); + if (H[8] == 0.0) { + return false; + } else { + // normalize + double f = 1.0 / H[8]; + // for (i = 0; i < 8; i++) mat[i] = f * H[i]; + mat[0] = f * H[2]; + mat[1] = f * H[5]; + mat[2] = f * H[0]; + mat[3] = f * H[1]; + mat[4] = f * H[3]; + mat[5] = f * H[4]; + mat[6] = f * H[6]; + mat[7] = f * H[7]; + } + return true; +} +#elif VERTTRAP_ALGORITHM == 2 +static bool find_vertrapezoid(int np, double *pts1, double *pts2, double *mat) { + // Based on straight Least-squares + const int np2 = np * 2; + const int nvar = 6; + double *a = + (double *)aom_malloc(sizeof(*a) * (np2 * (nvar + 1) + (nvar + 1) * nvar)); + if (a == NULL) return false; + double *b = a + np2 * nvar; + double *temp = b + np2; + int i; + double sx, sy, dx, dy; + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * nvar + 0] = sx; + a[i * 2 * nvar + 1] = 1; + a[i * 2 * nvar + 2] = 0; + a[i * 2 * nvar + 3] = 0; + a[i * 2 * nvar + 4] = 0; + a[i * 2 * nvar + 5] = -dx * sx; + + a[(i * 2 + 1) * nvar + 0] = 0; + a[(i * 2 + 1) * nvar + 1] = 0; + a[(i * 2 + 1) * nvar + 2] = sx; + a[(i * 2 + 1) * nvar + 3] = sy; + a[(i * 2 + 1) * nvar + 4] = 1; + a[(i * 2 + 1) * nvar + 5] = -dy * sx; + + b[2 * i] = dx; + b[2 * i + 1] = dy; + } + double sol[8]; + if (!least_squares(nvar, a, np2, nvar, b, temp, sol)) { + aom_free(a); + return false; + } + mat[0] = sol[1]; + mat[1] = sol[4]; + mat[2] = sol[0]; + mat[3] = 0; + mat[4] = sol[2]; + mat[5] = sol[3]; + mat[6] = sol[5]; + mat[7] = 0; + aom_free(a); + return true; +} +#else +#error "Invalid value of VERTTRAP_ALGORITHM" +#endif + +#if HORZTRAP_ALGORITHM == 0 +static bool find_hortrapezoid(int np, double *pts1, double *pts2, double *mat) { + // Implemented from Peter Kovesi's normalized implementation + const int nvar = 7; + const int np3 = np * 3; + double *a = (double *)aom_malloc(sizeof(*a) * np3 * nvar * 2); + double *U = a + np3 * nvar; + double S[7], V[7 * 7]; + int i, mini; + double sx, sy, dx, dy; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 3 * nvar + 0] = 0; + a[i * 3 * nvar + 1] = 0; + a[i * 3 * nvar + 2] = 0; + a[i * 3 * nvar + 3] = -sy; + a[i * 3 * nvar + 4] = -1; + a[i * 3 * nvar + 5] = dy * sy; + a[i * 3 * nvar + 6] = dy; + + a[(i * 3 + 1) * nvar + 0] = sx; + a[(i * 3 + 1) * nvar + 1] = sy; + a[(i * 3 + 1) * nvar + 2] = 1; + a[(i * 3 + 1) * nvar + 3] = 0; + a[(i * 3 + 1) * nvar + 4] = 0; + a[(i * 3 + 1) * nvar + 5] = -dx * sy; + a[(i * 3 + 1) * nvar + 6] = -dx; + + a[(i * 3 + 2) * nvar + 0] = -dy * sx; + a[(i * 3 + 2) * nvar + 1] = -dy * sy; + a[(i * 3 + 2) * nvar + 2] = -dy; + a[(i * 3 + 2) * nvar + 3] = dx * sy; + a[(i * 3 + 2) * nvar + 4] = dx; + a[(i * 3 + 2) * nvar + 5] = 0; + a[(i * 3 + 2) * nvar + 6] = 0; + } + + if (SVD(U, S, V, a, np3, nvar)) { + aom_free(a); + return false; + } else { + double minS = 1e12; + mini = -1; + for (i = 0; i < nvar; ++i) { + if (S[i] < minS) { + minS = S[i]; + mini = i; + } + } + } + double H[9]; + H[0] = V[0 * nvar + mini]; + H[1] = V[1 * nvar + mini]; + H[2] = V[2 * nvar + mini]; + H[3] = 0; + H[4] = V[3 * nvar + mini]; + H[5] = V[4 * nvar + mini]; + H[6] = 0; + H[7] = V[5 * nvar + mini]; + H[8] = V[6 * nvar + mini]; + // denormalize_homography_reorder(H, T1, T2); + aom_free(a); + if (H[8] == 0.0) { + return false; + } else { + // normalize + double f = 1.0 / H[8]; + // for (i = 0; i < 8; i++) mat[i] = f * H[i]; + mat[0] = f * H[2]; + mat[1] = f * H[5]; + mat[2] = f * H[0]; + mat[3] = f * H[1]; + mat[4] = f * H[3]; + mat[5] = f * H[4]; + mat[6] = f * H[6]; + mat[7] = f * H[7]; + } + return true; +} +#elif HORZTRAP_ALGORITHM == 1 +static bool find_hortrapezoid(int np, double *pts1, double *pts2, double *mat) { + // Based on SVD decomposition of homogeneous equation and using the right + // unitary vector corresponding to the smallest singular value + const int nvar = 7; + const int np2 = np * 2; + double *a = (double *)aom_malloc(sizeof(*a) * np2 * nvar * 2); + double *U = a + np2 * nvar; + double S[7], V[7 * 7]; + int i, mini; + double sx, sy, dx, dy; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * nvar + 0] = 0; + a[i * 2 * nvar + 1] = 0; + a[i * 2 * nvar + 2] = 0; + a[i * 2 * nvar + 3] = -sy; + a[i * 2 * nvar + 4] = -1; + a[i * 2 * nvar + 5] = dy * sy; + a[i * 2 * nvar + 6] = dy; + + a[(i * 2 + 1) * nvar + 0] = -sx; + a[(i * 2 + 1) * nvar + 1] = -sy; + a[(i * 2 + 1) * nvar + 2] = -1; + a[(i * 2 + 1) * nvar + 3] = 0; + a[(i * 2 + 1) * nvar + 4] = 0; + a[(i * 2 + 1) * nvar + 5] = dx * sy; + a[(i * 2 + 1) * nvar + 6] = dx; + } + + if (SVD(U, S, V, a, np2, nvar)) { + aom_free(a); + return false; + } else { + double minS = 1e12; + mini = -1; + for (i = 0; i < nvar; ++i) { + if (S[i] < minS) { + minS = S[i]; + mini = i; + } + } + } + + double H[9]; + H[0] = V[0 * nvar + mini]; + H[1] = V[1 * nvar + mini]; + H[2] = V[2 * nvar + mini]; + H[3] = 0; + H[4] = V[3 * nvar + mini]; + H[5] = V[4 * nvar + mini]; + H[6] = 0; + H[7] = V[5 * nvar + mini]; + H[8] = V[6 * nvar + mini]; + // denormalize_homography_reorder(H, T1, T2); + aom_free(a); + if (H[8] == 0.0) { + return false; + } else { + // normalize + double f = 1.0 / H[8]; + // for (i = 0; i < 8; i++) mat[i] = f * H[i]; + mat[0] = f * H[2]; + mat[1] = f * H[5]; + mat[2] = f * H[0]; + mat[3] = f * H[1]; + mat[4] = f * H[3]; + mat[5] = f * H[4]; + mat[6] = f * H[6]; + mat[7] = f * H[7]; + } + return true; +} +#elif HORZTRAP_ALGORITHM == 2 +static bool find_hortrapezoid(int np, double *pts1, double *pts2, double *mat) { + // Based on straight Least-squares + const int np2 = np * 2; + const int nvar = 8; + double *a = + (double *)aom_malloc(sizeof(*a) * (np2 * (nvar + 1) + (nvar + 1) * nvar)); + if (a == NULL) return false; + double *b = a + np2 * nvar; + double *temp = b + np2; + int i; + double sx, sy, dx, dy; + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * nvar + 0] = sx; + a[i * 2 * nvar + 1] = sy; + a[i * 2 * nvar + 2] = 1; + a[i * 2 * nvar + 3] = 0; + a[i * 2 * nvar + 4] = 0; + a[i * 2 * nvar + 5] = -dx * sy; + + a[(i * 2 + 1) * nvar + 0] = 0; + a[(i * 2 + 1) * nvar + 1] = 0; + a[(i * 2 + 1) * nvar + 2] = 0; + a[(i * 2 + 1) * nvar + 3] = sy; + a[(i * 2 + 1) * nvar + 4] = 1; + a[(i * 2 + 1) * nvar + 5] = -dy * sy; + + b[2 * i] = dx; + b[2 * i + 1] = dy; + } + double sol[8]; + if (!least_squares(nvar, a, np2, nvar, b, temp, sol)) { + aom_free(a); + return false; + } + mat[0] = sol[2]; + mat[1] = sol[4]; + mat[2] = sol[0]; + mat[3] = sol[1]; + mat[4] = 0.0; + mat[5] = sol[3]; + mat[6] = 0.0; + mat[7] = sol[5]; + aom_free(a); + return true; +} +#else +#error "Invalid value of HORZTRAP_ALGORITHM" +#endif + +#if HOMOGRAPHY_ALGORITHM == 0 +static bool find_homography(int np, double *pts1, double *pts2, double *mat) { + // Implemented from Peter Kovesi's normalized implementation + const int np3 = np * 3; + double *a = (double *)aom_malloc(sizeof(*a) * np3 * 18); + double *U = a + np3 * 9; + double S[9], V[9 * 9], H[9]; + int i, mini; + double sx, sy, dx, dy; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 3 * 9 + 0] = a[i * 3 * 9 + 1] = a[i * 3 * 9 + 2] = 0; + a[i * 3 * 9 + 3] = -sx; + a[i * 3 * 9 + 4] = -sy; + a[i * 3 * 9 + 5] = -1; + a[i * 3 * 9 + 6] = dy * sx; + a[i * 3 * 9 + 7] = dy * sy; + a[i * 3 * 9 + 8] = dy; + + a[(i * 3 + 1) * 9 + 0] = sx; + a[(i * 3 + 1) * 9 + 1] = sy; + a[(i * 3 + 1) * 9 + 2] = 1; + a[(i * 3 + 1) * 9 + 3] = a[(i * 3 + 1) * 9 + 4] = a[(i * 3 + 1) * 9 + 5] = + 0; + a[(i * 3 + 1) * 9 + 6] = -dx * sx; + a[(i * 3 + 1) * 9 + 7] = -dx * sy; + a[(i * 3 + 1) * 9 + 8] = -dx; + + a[(i * 3 + 2) * 9 + 0] = -dy * sx; + a[(i * 3 + 2) * 9 + 1] = -dy * sy; + a[(i * 3 + 2) * 9 + 2] = -dy; + a[(i * 3 + 2) * 9 + 3] = dx * sx; + a[(i * 3 + 2) * 9 + 4] = dx * sy; + a[(i * 3 + 2) * 9 + 5] = dx; + a[(i * 3 + 2) * 9 + 6] = a[(i * 3 + 2) * 9 + 7] = a[(i * 3 + 2) * 9 + 8] = + 0; + } + + if (SVD(U, S, V, a, np3, 9)) { + aom_free(a); + return false; + } else { + double minS = 1e12; + mini = -1; + for (i = 0; i < 9; ++i) { + if (S[i] < minS) { + minS = S[i]; + mini = i; + } + } + } + + for (i = 0; i < 9; i++) H[i] = V[i * 9 + mini]; + // denormalize_homography_reorder(H, T1, T2); + aom_free(a); + if (H[8] == 0.0) { + return false; + } else { + // normalize + double f = 1.0 / H[8]; + // for (i = 0; i < 8; i++) mat[i] = f * H[i]; + mat[0] = f * H[2]; + mat[1] = f * H[5]; + mat[2] = f * H[0]; + mat[3] = f * H[1]; + mat[4] = f * H[3]; + mat[5] = f * H[4]; + mat[6] = f * H[6]; + mat[7] = f * H[7]; + } + return true; +} +#elif HOMOGRAPHY_ALGORITHM == 1 +static bool find_homography(int np, double *pts1, double *pts2, double *mat) { + // Based on SVD decomposition of homogeneous equation and using the right + // unitary vector corresponding to the smallest singular value + const int np2 = np * 2; + double *a = (double *)aom_malloc(sizeof(*a) * np2 * 18); + double *U = a + np2 * 9; + double S[9], V[9 * 9], H[9]; + int i, mini; + double sx, sy, dx, dy; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * 9 + 0] = a[i * 2 * 9 + 1] = a[i * 2 * 9 + 2] = 0; + a[i * 2 * 9 + 3] = -sx; + a[i * 2 * 9 + 4] = -sy; + a[i * 2 * 9 + 5] = -1; + a[i * 2 * 9 + 6] = dy * sx; + a[i * 2 * 9 + 7] = dy * sy; + a[i * 2 * 9 + 8] = dy; + + a[(i * 2 + 1) * 9 + 0] = -sx; + a[(i * 2 + 1) * 9 + 1] = -sy; + a[(i * 2 + 1) * 9 + 2] = -1; + a[(i * 2 + 1) * 9 + 3] = a[(i * 2 + 1) * 9 + 4] = a[(i * 2 + 1) * 9 + 5] = + 0; + a[(i * 2 + 1) * 9 + 6] = dx * sx; + a[(i * 2 + 1) * 9 + 7] = dx * sy; + a[(i * 2 + 1) * 9 + 8] = dx; + } + + if (SVD(U, S, V, a, np2, 9)) { + aom_free(a); + return false; + } else { + double minS = 1e12; + mini = -1; + for (i = 0; i < 9; ++i) { + if (S[i] < minS) { + minS = S[i]; + mini = i; + } + } + } + + for (i = 0; i < 9; i++) H[i] = V[i * 9 + mini]; + // denormalize_homography_reorder(H, T1, T2); + aom_free(a); + if (H[8] == 0.0) { + return false; + } else { + // normalize + double f = 1.0 / H[8]; + // for (i = 0; i < 8; i++) mat[i] = f * H[i]; + mat[0] = f * H[2]; + mat[1] = f * H[5]; + mat[2] = f * H[0]; + mat[3] = f * H[1]; + mat[4] = f * H[3]; + mat[5] = f * H[4]; + mat[6] = f * H[6]; + mat[7] = f * H[7]; + } + return true; +} +#elif HOMOGRAPHY_ALGORITHM == 2 +static bool find_homography(int np, double *pts1, double *pts2, double *mat) { + // Based on straight Least-squares + const int np2 = np * 2; + const int nvar = 8; + double *a = + (double *)aom_malloc(sizeof(*a) * (np2 * (nvar + 1) + (nvar + 1) * nvar)); + if (a == NULL) return false; + double *b = a + np2 * nvar; + double *temp = b + np2; + int i; + double sx, sy, dx, dy; + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * nvar + 0] = sx; + a[i * 2 * nvar + 1] = sy; + a[i * 2 * nvar + 2] = 1; + a[i * 2 * nvar + 3] = 0; + a[i * 2 * nvar + 4] = 0; + a[i * 2 * nvar + 5] = 0; + a[i * 2 * nvar + 6] = -dx * sx; + a[i * 2 * nvar + 7] = -dx * sy; + + a[(i * 2 + 1) * nvar + 0] = 0; + a[(i * 2 + 1) * nvar + 1] = 0; + a[(i * 2 + 1) * nvar + 2] = 0; + a[(i * 2 + 1) * nvar + 3] = sx; + a[(i * 2 + 1) * nvar + 4] = sy; + a[(i * 2 + 1) * nvar + 5] = 1; + a[(i * 2 + 1) * nvar + 6] = -dy * sx; + a[(i * 2 + 1) * nvar + 7] = -dy * sy; + + b[2 * i] = dx; + b[2 * i + 1] = dy; + } + double sol[8]; + if (!least_squares(nvar, a, np2, nvar, b, temp, sol)) { + aom_free(a); + return false; + } + mat[0] = sol[2]; + mat[1] = sol[5]; + mat[2] = sol[0]; + mat[3] = sol[1]; + mat[4] = sol[3]; + mat[5] = sol[4]; + mat[6] = sol[6]; + mat[7] = sol[7]; + aom_free(a); + return true; +} +#else +#error "Invalid value of HOMOGRAPHY_ALGORITHM" +#endif // HOMOGRAPHY_ALGORITHM + +static bool find_translation(int np, double *pts1, double *pts2, double *mat) { + int i; + double sx, sy, dx, dy; + double sumx, sumy; + + double T1[9], T2[9]; + normalize_homography(pts1, np, T1); + normalize_homography(pts2, np, T2); + + sumx = 0; + sumy = 0; + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + sumx += dx - sx; + sumy += dy - sy; + } + mat[0] = sumx / np; + mat[1] = sumy / np; + denormalize_translation_reorder(mat, T1, T2); + return true; +} + +static bool find_rotzoom(int np, double *pts1, double *pts2, double *mat) { + const int np2 = np * 2; + double *a = (double *)aom_malloc(sizeof(*a) * (np2 * 5 + 20)); + double *b = a + np2 * 4; + double *temp = b + np2; + int i; + double sx, sy, dx, dy; + + double T1[9], T2[9]; + normalize_homography(pts1, np, T1); + normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * 4 + 0] = sx; + a[i * 2 * 4 + 1] = sy; + a[i * 2 * 4 + 2] = 1; + a[i * 2 * 4 + 3] = 0; + a[(i * 2 + 1) * 4 + 0] = sy; + a[(i * 2 + 1) * 4 + 1] = -sx; + a[(i * 2 + 1) * 4 + 2] = 0; + a[(i * 2 + 1) * 4 + 3] = 1; + + b[2 * i] = dx; + b[2 * i + 1] = dy; + } + if (!least_squares(4, a, np2, 4, b, temp, mat)) { + aom_free(a); + return false; + } + denormalize_rotzoom_reorder(mat, T1, T2); + aom_free(a); + return true; +} + +static bool find_affine(int np, double *pts1, double *pts2, double *mat) { + assert(np > 0); + const int np2 = np * 2; + double *a = (double *)aom_malloc(sizeof(*a) * (np2 * 7 + 42)); + if (a == NULL) return false; + double *b = a + np2 * 6; + double *temp = b + np2; + int i; + double sx, sy, dx, dy; + + double T1[9], T2[9]; + normalize_homography(pts1, np, T1); + normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * 6 + 0] = sx; + a[i * 2 * 6 + 1] = sy; + a[i * 2 * 6 + 2] = 0; + a[i * 2 * 6 + 3] = 0; + a[i * 2 * 6 + 4] = 1; + a[i * 2 * 6 + 5] = 0; + a[(i * 2 + 1) * 6 + 0] = 0; + a[(i * 2 + 1) * 6 + 1] = 0; + a[(i * 2 + 1) * 6 + 2] = sx; + a[(i * 2 + 1) * 6 + 3] = sy; + a[(i * 2 + 1) * 6 + 4] = 0; + a[(i * 2 + 1) * 6 + 5] = 1; + + b[2 * i] = dx; + b[2 * i + 1] = dy; + } + if (!least_squares(6, a, np2, 6, b, temp, mat)) { + aom_free(a); + return false; + } + denormalize_affine_reorder(mat, T1, T2); + aom_free(a); + return true; +} + +static bool find_rotation(int np, double *pts1, double *pts2, double *mat) { + // Note(rachelbarker): + // Unlike the other model types, a rotational model has a nonlinear + // constraint: The output model must satisfy + // mat[2] * mat[2] + mat[3] * mat[3] = 1 + // Thus we cannot use the same linear least-squares approach as the + // other model types. However, we can use an alternative algorithm + // called the Kabsch algorithm to solve this problem. + + double mean1[2] = { 0.0, 0.0 }; + double mean2[2] = { 0.0, 0.0 }; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + double *p, *q; + double inp = 1.0 / np; + int i; + for (i = 0, p = pts1; i < np; ++i, p += 2) { + mean1[0] += p[0]; + mean1[1] += p[1]; + } + mean1[0] *= inp; + mean1[1] *= inp; + for (i = 0, p = pts2; i < np; ++i, p += 2) { + mean2[0] += p[0]; + mean2[1] += p[1]; + } + mean2[0] *= inp; + mean2[1] *= inp; + double A[4] = { 0.0, 0.0, 0.0, 0.0 }; + for (p = pts1, q = pts2, i = 0; i < np; ++i, p += 2, q += 2) { + A[0] += (p[0] - mean1[0]) * (q[0] - mean2[0]); + A[1] += (p[0] - mean1[0]) * (q[1] - mean2[1]); + A[2] += (p[1] - mean1[1]) * (q[0] - mean2[0]); + A[3] += (p[1] - mean1[1]) * (q[1] - mean2[1]); + } + double V[4], S[2], W[4]; + if (SVD(V, S, W, A, 2, 2)) return false; + // printf("V: %f %f %f %f\n", V[0], V[1], V[2], V[3]); + // printf("S: %f %f\n", S[0], S[1]); + // printf("W: %f %f %f %f\n", W[0], W[1], W[2], W[3]); + double detA = A[0] * A[3] - A[1] * A[2]; + if (detA < 0) { + V[1] = -V[1]; + V[3] = -V[3]; + } + mat[2] = W[0] * V[0] + W[1] * V[1]; + mat[3] = W[0] * V[2] + W[1] * V[3]; + mat[4] = W[2] * V[0] + W[3] * V[1]; + mat[5] = W[2] * V[2] + W[3] * V[3]; + mat[6] = mat[7] = 0.0; + mat[0] = mean2[0] - mean1[0] * mat[2] - mean1[1] * mat[3]; + mat[1] = mean2[1] - mean1[0] * mat[4] - mean1[1] * mat[5]; + // denormalize_homography_general_reorder(mat, T1, T2); + return true; +} + +static bool find_zoom(int np, double *pts1, double *pts2, double *mat) { + const int np2 = np * 2; + double *a = (double *)aom_malloc(sizeof(*a) * (np2 * 4 + 12)); + double *b = a + np2 * 3; + double *temp = b + np2; + int i; + double sx, sy, dx, dy; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * 3 + 0] = sx; + a[i * 2 * 3 + 1] = 1; + a[i * 2 * 3 + 2] = 0; + a[(i * 2 + 1) * 3 + 0] = sy; + a[(i * 2 + 1) * 3 + 1] = 0; + a[(i * 2 + 1) * 3 + 2] = 1; + + b[2 * i] = dx; + b[2 * i + 1] = dy; + } + double sol[3]; + if (!least_squares(3, a, np2, 3, b, temp, sol)) { + aom_free(a); + return false; + } + // denormalize_zoom_reorder(mat, T1, T2); + mat[0] = sol[1]; + mat[1] = sol[2]; + mat[2] = mat[5] = sol[0]; + mat[3] = mat[4] = mat[6] = mat[7] = 0.0; + + aom_free(a); + return true; +} + +static bool find_uzoom(int np, double *pts1, double *pts2, double *mat) { + const int np2 = np * 2; + const int nvar = 4; + double *a = + (double *)aom_malloc(sizeof(*a) * (np2 * (nvar + 1) + (nvar + 1) * nvar)); + if (a == NULL) return false; + double *b = a + np2 * nvar; + double *temp = b + np2; + int i; + double sx, sy, dx, dy; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (i = 0; i < np; ++i) { + dx = *(pts2++); + dy = *(pts2++); + sx = *(pts1++); + sy = *(pts1++); + + a[i * 2 * nvar + 0] = sx; + a[i * 2 * nvar + 1] = 0; + a[i * 2 * nvar + 2] = 1; + a[i * 2 * nvar + 3] = 0; + a[(i * 2 + 1) * nvar + 0] = 0; + a[(i * 2 + 1) * nvar + 1] = sy; + a[(i * 2 + 1) * nvar + 2] = 0; + a[(i * 2 + 1) * nvar + 3] = 1; + + b[2 * i] = dx; + b[2 * i + 1] = dy; + } + double sol[4]; + if (!least_squares(nvar, a, np2, nvar, b, temp, sol)) { + aom_free(a); + return false; + } + // denormalize_rotzoom_reorder(mat, T1, T2); + mat[0] = sol[2]; + mat[1] = sol[3]; + mat[2] = sol[0]; + mat[3] = mat[4] = 0; + mat[5] = sol[1]; + mat[6] = mat[7] = 0.0; + aom_free(a); + return true; +} + +static bool find_rotuzoom(int np, double *pts1, double *pts2, double *mat) { + // The affine matrix is assumed to be the product of a rotation matrix by + // theta, and a zoom matrix of the form: ( zx 0 + // 0 zy ) + // So the resultant affine matrix is of the form: + // ( a bt + // -at b ) + // where a = zx * cos(theta), b = zy * cos(theta), t = tan(theta) + // We are required to find the best (a, b, t) values and the best motion + // vector (vx, vy) so that the error in projection of the points (x, y) to + // (x', y') following: + // ( x' ) = ( a bt ) * ( x ) + ( vx ) + // ( y' ) (-at b ) ( y ) ( vy ) + // is minimized. + // + // This optimizer uses a gradient descent algorithm in the (a, b, t) space. + // For a given (a, b, t) the optimal motion vector (vx, vy) can be computed + // by setting the derivatives of the projection error to 0. Therefore it + // is sufficient to run graduient descent in the (a, b, t) 3-parameter space. + // + double Sx = 0.0; // mean of source x + double Sy = 0.0; // mean of source y + double Px = 0.0; // mean of projected x + double Py = 0.0; // mean of projected y + double Sxx = 0.0; // mean of source x^2 + double Syy = 0.0; // mean of source y^2 + double Kxx = 0.0; // mean of source x * projected x + double Kxy = 0.0; // mean of source x * projected y + double Kyx = 0.0; // mean of source y * projected x + double Kyy = 0.0; // mean of source y * projected y + for (int i = 0; i < np; ++i) { + const double dx = *(pts2++); + const double dy = *(pts2++); + const double sx = *(pts1++); + const double sy = *(pts1++); + + Sx += sx; + Sy += sy; + Px += dx; + Py += dy; + Sxx += sx * sx; + Syy += sy * sy; + + Kxx += sx * dx; + Kxy += sx * dy; + Kyx += sy * dx; + Kyy += sy * dy; + } + Sx /= np; + Sy /= np; + Sxx /= np; + Syy /= np; + Px /= np; + Py /= np; + Kxx /= np; + Kxy /= np; + Kyx /= np; + Kyy /= np; + + // Step size + // + // By using a large initial step size, we can rapidly search the parameter + // space for a good model. However, gradient descent with a large step size + // can end up oscillating around the solution rather than converging. + // We detect that situation and reduce alpha when it occurs, so that we + // can converge in on the minimum which has been located. + double alpha = 1.0; + + const int iters_thresh = 1000; + // Threshold for deciding when we're at a minimum + const double termination_threshold = 1e-5; + // Threshold for detecting oscillatory behaviour + const double oscillation_threshold = -0.90; + + // Initialize z = (a, b, t) + double z[3] = { 1, 1, 0 }; + // Derivatives + double dz[3]; + double dz_prev[3] = { 0.0, 0.0, 0.0 }; + // Motion vector + double v[2]; + + int iters = 0; + while (1) { + const double a = z[0]; + const double b = z[1]; + const double t = z[2]; + // Optimal motion vector obtained by setting partial derivatives to 0 + v[0] = Px - a * Sx - b * t * Sy; + v[1] = Py + a * t * Sx - b * Sy; + // These are from partial derivatives of the projection error + dz[0] = + 2 * (a * (1 + t * t) * Sxx + (v[0] - v[1] * t) * Sx - Kxx + t * Kxy); + dz[1] = + 2 * (b * (1 + t * t) * Syy + (v[0] * t + v[1]) * Sy - Kyy - t * Kyx); + dz[2] = 2 * (t * (b * b * Syy + a * a * Sxx) + v[0] * b * Sy - + a * v[1] * Sx - b * Kyx + a * Kxy); + + // Test termination criteria + double dz_norm = norm(dz, 3); + if (iters >= iters_thresh) { + // Could not find a good enough model + return false; + } else if (dz_norm < termination_threshold) { + // At a local minimum or saddle point + break; + } + + // Normalize partial derivative vector + dz[0] /= dz_norm; + dz[1] /= dz_norm; + dz[2] /= dz_norm; + + // Decide when to reduce step size + // + // The gradient descent method with a fixed step size tends to oscillate + // around the solution, so we check for cases where the normalized gradient + // vector reverses between iterations. + // + // Since dz and dz_prev are both normalized, we have + // dot(dz, dz_prev) = cos(angle between dz and dz_prev) + // + // Then there are a few cases to think about: + // 1) When walking toward a minimum, dz and dz_prev will be in similar + // directions, so cos(angle) is positive + // 2) If we're spiralling in toward a minimum, then cos(angle) will be + // negative but small + // 3) If we're oscillating around a minimum, then cos(angle) will be + // close to -1 + // + // So our oscillation criterion is that dot(dz, dz_prev) is sufficiently + // close to -1. + double dot = dz[0] * dz_prev[0] + dz[1] * dz_prev[1] + dz[2] * dz_prev[2]; + if (dot < oscillation_threshold) { + alpha *= 0.5; + } + + // Gradient Descent Updates + z[0] -= alpha * dz[0]; + z[1] -= alpha * dz[1]; + z[2] -= alpha * dz[2]; + + // Prepare for next iteration + memcpy(dz_prev, dz, sizeof(dz)); + iters++; + } + + mat[0] = v[0]; + mat[1] = v[1]; + mat[2] = z[0]; + mat[3] = z[1] * z[2]; + mat[4] = -z[0] * z[2]; + mat[5] = z[1]; + mat[6] = mat[7] = 0.0; + return true; +} + +static bool find_vertshear(int np, double *pts1, double *pts2, double *mat) { + const int nvar = 3; + const int np2 = np * 2; + double *a = + (double *)aom_malloc(sizeof(*a) * (np2 * (nvar + 1) + (nvar + 1) * nvar)); + if (a == NULL) return false; + double *b = a + np2 * nvar; + double *temp = b + np2; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (int i = 0; i < np; ++i) { + const double dx = *(pts2++); + const double dy = *(pts2++); + const double sx = *(pts1++); + const double sy = *(pts1++); + + a[i * 2 * nvar + 0] = 0; + a[i * 2 * nvar + 1] = 1; + a[i * 2 * nvar + 2] = 0; + a[(i * 2 + 1) * nvar + 0] = sx; + a[(i * 2 + 1) * nvar + 1] = 0; + a[(i * 2 + 1) * nvar + 2] = 1; + + b[2 * i] = dx - sx; + b[2 * i + 1] = dy - sy; + } + double sol[3]; + if (!least_squares(nvar, a, np2, nvar, b, temp, sol)) { + aom_free(a); + return false; + } + // denormalize_zoom_reorder(mat, T1, T2); + mat[0] = sol[1]; + mat[1] = sol[2]; + mat[2] = 1.0; + mat[3] = 0; + mat[4] = sol[0]; + mat[5] = 1.0; + mat[6] = mat[7] = 0.0; + aom_free(a); + return true; +} + +static bool find_horzshear(int np, double *pts1, double *pts2, double *mat) { + const int nvar = 3; + const int np2 = np * 2; + double *a = + (double *)aom_malloc(sizeof(*a) * (np2 * (nvar + 1) + (nvar + 1) * nvar)); + if (a == NULL) return false; + double *b = a + np2 * nvar; + double *temp = b + np2; + + // double T1[9], T2[9]; + // normalize_homography(pts1, np, T1); + // normalize_homography(pts2, np, T2); + + for (int i = 0; i < np; ++i) { + const double dx = *(pts2++); + const double dy = *(pts2++); + const double sx = *(pts1++); + const double sy = *(pts1++); + + a[i * 2 * nvar + 0] = sy; + a[i * 2 * nvar + 1] = 1; + a[i * 2 * nvar + 2] = 0; + a[(i * 2 + 1) * nvar + 0] = 0; + a[(i * 2 + 1) * nvar + 1] = 0; + a[(i * 2 + 1) * nvar + 2] = 1; + + b[2 * i] = dx - sx; + b[2 * i + 1] = dy - sy; + } + double sol[3]; + if (!least_squares(nvar, a, np2, nvar, b, temp, sol)) { + aom_free(a); + return false; + } + // denormalize_zoom_reorder(mat, T1, T2); + mat[0] = sol[1]; + mat[1] = sol[2]; + mat[2] = 1.0; + mat[3] = sol[0]; + mat[4] = 0.0; + mat[5] = 1.0; + mat[6] = mat[7] = 0.0; + aom_free(a); + return true; +} + +// Returns true on success, false if not enough points provided +static bool get_rand_indices(int npoints, int minpts, int *indices, + unsigned int *seed) { + int i, j; + int ptr = lcg_rand16(seed) % npoints; + if (minpts > npoints) return false; + indices[0] = ptr; + ptr = (ptr == npoints - 1 ? 0 : ptr + 1); + i = 1; + while (i < minpts) { + int index = lcg_rand16(seed) % npoints; + while (index) { + ptr = (ptr == npoints - 1 ? 0 : ptr + 1); + for (j = 0; j < i; ++j) { + if (indices[j] == ptr) break; + } + if (j == i) index--; + } + indices[i++] = ptr; + } + return true; +} + +typedef struct { + int num_inliers; + double variance; + int *inlier_indices; +} RANSAC_MOTION; + +// Return -1 if 'a' is a better motion, 1 if 'b' is better, 0 otherwise. +static int compare_motions(const void *arg_a, const void *arg_b) { + const RANSAC_MOTION *motion_a = (RANSAC_MOTION *)arg_a; + const RANSAC_MOTION *motion_b = (RANSAC_MOTION *)arg_b; + + if (motion_a->num_inliers > motion_b->num_inliers) return -1; + if (motion_a->num_inliers < motion_b->num_inliers) return 1; + if (motion_a->variance < motion_b->variance) return -1; + if (motion_a->variance > motion_b->variance) return 1; + return 0; +} + +static bool is_better_motion(const RANSAC_MOTION *motion_a, + const RANSAC_MOTION *motion_b) { + return compare_motions(motion_a, motion_b) < 0; +} + +static void copy_points_at_indices(double *dest, const double *src, + const int *indices, int num_points) { + for (int i = 0; i < num_points; ++i) { + const int index = indices[i]; + dest[i * 2] = src[index * 2]; + dest[i * 2 + 1] = src[index * 2 + 1]; + } +} + +static const double kInfiniteVariance = 1e12; + +static void clear_motion(RANSAC_MOTION *motion, int num_points) { + motion->num_inliers = 0; + motion->variance = kInfiniteVariance; + memset(motion->inlier_indices, 0, + sizeof(*motion->inlier_indices) * num_points); +} + +// Returns true on success, false on error +static bool ransac_internal(const Correspondence *matched_points, int npoints, + MotionModel *params_by_motion, + int num_desired_motions, int minpts, + IsDegenerateFunc is_degenerate, + FindTransformationFunc find_transformation, + ProjectPointsFunc projectpoints) { + int trial_count = 0; + int i = 0; + bool ret_val = true; + + unsigned int seed = (unsigned int)npoints; + + int indices[MAX_MINPTS] = { 0 }; + + double *points1, *points2; + double *corners1, *corners2; + double *image1_coord; + + // Store information for the num_desired_motions best transformations found + // and the worst motion among them, as well as the motion currently under + // consideration. + RANSAC_MOTION *motions, *worst_kept_motion = NULL; + RANSAC_MOTION current_motion; + + // Store the parameters and the indices of the inlier points for the motion + // currently under consideration. + double params_this_motion[MAX_PARAMDIM]; + + double *cnp1, *cnp2; + + for (i = 0; i < num_desired_motions; ++i) { + params_by_motion[i].num_inliers = 0; + } + if (npoints < minpts * MINPTS_MULTIPLIER || npoints == 0) { + return 1; + } + + points1 = (double *)aom_malloc(sizeof(*points1) * npoints * 2); + points2 = (double *)aom_malloc(sizeof(*points2) * npoints * 2); + corners1 = (double *)aom_malloc(sizeof(*corners1) * npoints * 2); + corners2 = (double *)aom_malloc(sizeof(*corners2) * npoints * 2); + image1_coord = (double *)aom_malloc(sizeof(*image1_coord) * npoints * 2); + + motions = + (RANSAC_MOTION *)aom_malloc(sizeof(RANSAC_MOTION) * num_desired_motions); + for (i = 0; i < num_desired_motions; ++i) { + motions[i].inlier_indices = + (int *)aom_malloc(sizeof(*motions->inlier_indices) * npoints); + clear_motion(motions + i, npoints); + } + current_motion.inlier_indices = + (int *)aom_malloc(sizeof(*current_motion.inlier_indices) * npoints); + clear_motion(¤t_motion, npoints); + + worst_kept_motion = motions; + + if (!(points1 && points2 && corners1 && corners2 && image1_coord && motions && + current_motion.inlier_indices)) { + ret_val = false; + goto finish_ransac; + } + + cnp1 = corners1; + cnp2 = corners2; + for (i = 0; i < npoints; ++i) { + cnp1[2 * i + 0] = matched_points[i].x; + cnp1[2 * i + 1] = matched_points[i].y; + cnp2[2 * i + 0] = matched_points[i].rx; + cnp2[2 * i + 1] = matched_points[i].ry; + } + + while (MIN_TRIALS > trial_count) { + double sum_distance = 0.0; + double sum_distance_squared = 0.0; + + clear_motion(¤t_motion, npoints); + + int degenerate = 1; + int num_degenerate_iter = 0; + + while (degenerate) { + num_degenerate_iter++; + if (!get_rand_indices(npoints, minpts, indices, &seed)) { + ret_val = false; + goto finish_ransac; + } + + copy_points_at_indices(points1, corners1, indices, minpts); + copy_points_at_indices(points2, corners2, indices, minpts); + + degenerate = is_degenerate(points1); + if (num_degenerate_iter > MAX_DEGENERATE_ITER) { + ret_val = false; + goto finish_ransac; + } + } + + if (!find_transformation(minpts, points1, points2, params_this_motion)) { + trial_count++; + continue; + } + + projectpoints(params_this_motion, corners1, image1_coord, npoints, 2, 2); + + for (i = 0; i < npoints; ++i) { + double dx = image1_coord[i * 2] - corners2[i * 2]; + double dy = image1_coord[i * 2 + 1] - corners2[i * 2 + 1]; + double distance = sqrt(dx * dx + dy * dy); + + if (distance < INLIER_THRESHOLD) { + current_motion.inlier_indices[current_motion.num_inliers++] = i; + sum_distance += distance; + sum_distance_squared += distance * distance; + } + } + + if (current_motion.num_inliers >= worst_kept_motion->num_inliers && + current_motion.num_inliers > 1) { + double mean_distance; + mean_distance = sum_distance / ((double)current_motion.num_inliers); + current_motion.variance = + sum_distance_squared / ((double)current_motion.num_inliers - 1.0) - + mean_distance * mean_distance * ((double)current_motion.num_inliers) / + ((double)current_motion.num_inliers - 1.0); + if (is_better_motion(¤t_motion, worst_kept_motion)) { + // This motion is better than the worst currently kept motion. Remember + // the inlier points and variance. The parameters for each kept motion + // will be recomputed later using only the inliers. + worst_kept_motion->num_inliers = current_motion.num_inliers; + worst_kept_motion->variance = current_motion.variance; + memcpy(worst_kept_motion->inlier_indices, current_motion.inlier_indices, + sizeof(*current_motion.inlier_indices) * npoints); + assert(npoints > 0); + // Determine the new worst kept motion and its num_inliers and variance. + for (i = 0; i < num_desired_motions; ++i) { + if (is_better_motion(worst_kept_motion, &motions[i])) { + worst_kept_motion = &motions[i]; + } + } + } + } + trial_count++; + } + + // Sort the motions, best first. + qsort(motions, num_desired_motions, sizeof(RANSAC_MOTION), compare_motions); + + // Recompute the motions using only the inliers. + for (i = 0; i < num_desired_motions; ++i) { + if (motions[i].num_inliers >= minpts) { + int num_inliers = motions[i].num_inliers; + copy_points_at_indices(points1, corners1, motions[i].inlier_indices, + num_inliers); + copy_points_at_indices(points2, corners2, motions[i].inlier_indices, + num_inliers); + + find_transformation(num_inliers, points1, points2, + params_by_motion[i].params); + + // Populate inliers array + for (int j = 0; j < num_inliers; j++) { + int index = motions[i].inlier_indices[j]; + const Correspondence *corr = &matched_points[index]; + params_by_motion[i].inliers[2 * j + 0] = (int)rint(corr->x); + params_by_motion[i].inliers[2 * j + 1] = (int)rint(corr->y); + } + } + params_by_motion[i].num_inliers = motions[i].num_inliers; + } + +finish_ransac: + aom_free(points1); + aom_free(points2); + aom_free(corners1); + aom_free(corners2); + aom_free(image1_coord); + aom_free(current_motion.inlier_indices); + for (i = 0; i < num_desired_motions; ++i) { + aom_free(motions[i].inlier_indices); + } + aom_free(motions); + + return ret_val; +} + +static bool is_collinear3(double *p1, double *p2, double *p3) { + static const double collinear_eps = 1e-3; + const double v = + (p2[0] - p1[0]) * (p3[1] - p1[1]) - (p2[1] - p1[1]) * (p3[0] - p1[0]); + return fabs(v) < collinear_eps; +} + +static bool is_degenerate_homography(double *p) { + return is_collinear3(p, p + 2, p + 4) || is_collinear3(p, p + 2, p + 6) || + is_collinear3(p, p + 4, p + 6) || is_collinear3(p + 2, p + 4, p + 6); +} + +static bool is_degenerate_translation(double *p) { + return (p[0] - p[2]) * (p[0] - p[2]) + (p[1] - p[3]) * (p[1] - p[3]) <= 2; +} + +static bool is_degenerate_affine(double *p) { + return is_collinear3(p, p + 2, p + 4); +} + +static IsDegenerateFunc is_degenerate[TRANS_TYPES] = { + NULL, // IDENTITY + is_degenerate_translation, // TRANSLATION + is_degenerate_affine, // ROTATION + is_degenerate_affine, // ZOOM + is_degenerate_affine, // VERTSHEAR + is_degenerate_affine, // HORZSHEAR + is_degenerate_affine, // UZOOM + is_degenerate_affine, // ROTZOOM + is_degenerate_affine, // ROTUZOOM + is_degenerate_affine, // AFFINE + is_degenerate_homography, // VERTRAPEZOID + is_degenerate_homography, // HORTRAPEZOID + is_degenerate_homography // HOMOGRAPHY +}; + +static FindTransformationFunc find_transform[TRANS_TYPES] = { + NULL, // IDENTITY + find_translation, // TRANSLATION + find_rotation, // ROTATION + find_zoom, // ZOOM + find_vertshear, // VERTSHEAR + find_horzshear, // HORZSHEAR + find_uzoom, // UZOOM + find_rotzoom, // ROTZOOM + find_rotuzoom, // ROTUZOOM + find_affine, // AFFINE + find_vertrapezoid, // VERTRAPEZOID + find_hortrapezoid, // HORTRAPEZOID + find_homography, // HOMOGRAPHY +}; + +static ProjectPointsFunc project_points[TRANS_TYPES] = { + NULL, // IDENTITY + project_points_translation, // TRANSLATION + project_points_affine, // ROTATION + project_points_affine, // ZOOM + project_points_affine, // VERTSHEAR + project_points_affine, // HORZSHEAR + project_points_affine, // UZOOM + project_points_affine, // ROTZOOM + project_points_affine, // ROTUZOOM + project_points_affine, // AFFINE + project_points_homography, // VERTRAPEZOID + project_points_homography, // HORTRAPEZOID + project_points_homography // HOMOGRAPHY +}; + +// Returns true on success, false on error +bool ransac(Correspondence *matched_points, int npoints, + TransformationType type, MotionModel *params_by_motion, + int num_desired_motions) { + assert(type > IDENTITY && type < TRANS_TYPES); + + int minpts = 3; + + return ransac_internal(matched_points, npoints, params_by_motion, + num_desired_motions, minpts, is_degenerate[type], + find_transform[type], project_points[type]); +} + +// Fit a specified type of motion model to a set of correspondences. +// The input consists of `np` points, where pts1 stores the source position +// and pts2 stores the destination position for each correspondence. +// The resulting model is stored in `mat`. +// Returns true on success, false on error +// +// Note: The input points lists may be modified during processing +bool aom_fit_motion_model(TransformationType type, int np, double *pts1, + double *pts2, double *mat) { + assert(type > IDENTITY && type < TRANS_TYPES); + return find_transform[type](np, pts1, pts2, mat); +}
diff --git a/aom_dsp/flow_estimation/ransac.h b/aom_dsp/flow_estimation/ransac.h new file mode 100644 index 0000000..c2da385 --- /dev/null +++ b/aom_dsp/flow_estimation/ransac.h
@@ -0,0 +1,46 @@ +/* + * Copyright (c) 2021, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#ifndef AOM_FLOW_ESTIMATION_RANSAC_H_ +#define AOM_FLOW_ESTIMATION_RANSAC_H_ + +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include <memory.h> +#include <stdbool.h> + +#include "aom_dsp/flow_estimation/flow_estimation.h" + +#ifdef __cplusplus +extern "C" { +#endif + +bool ransac(Correspondence *matched_points, int npoints, + TransformationType type, MotionModel *params_by_motion, + int num_desired_motions); + +// Fit a specified type of motion model to a set of correspondences. +// The input consists of `np` points, where pts1 stores the source position +// and pts2 stores the destination position for each correspondence. +// The resulting model is stored in `mat` +// Returns true on success, false on error +// +// Note: The input points lists are modified during processing +bool aom_fit_motion_model(TransformationType type, int np, double *pts1, + double *pts2, double *mat); + +#ifdef __cplusplus +} +#endif + +#endif // AOM_FLOW_ESTIMATION_RANSAC_H_
diff --git a/av1/encoder/x86/corner_match_avx2.c b/aom_dsp/flow_estimation/x86/corner_match_avx2.c similarity index 96% rename from av1/encoder/x86/corner_match_avx2.c rename to aom_dsp/flow_estimation/x86/corner_match_avx2.c index 8ee997a..6aae6ad 100644 --- a/av1/encoder/x86/corner_match_avx2.c +++ b/aom_dsp/flow_estimation/x86/corner_match_avx2.c
@@ -17,7 +17,7 @@ #include "aom_ports/mem.h" #include "aom_ports/system_state.h" -#include "av1/encoder/corner_match.h" +#include "aom_dsp/flow_estimation/corner_match.h" DECLARE_ALIGNED(16, static const uint8_t, byte_mask[16]) = { 255, 255, 255, 255, 255, 255, 255, 255, @@ -30,7 +30,7 @@ correlation/standard deviation are taken over MATCH_SZ by MATCH_SZ windows of each image, centered at (x1, y1) and (x2, y2) respectively. */ -double av1_compute_cross_correlation_avx2(unsigned char *im1, int stride1, +double aom_compute_cross_correlation_avx2(unsigned char *im1, int stride1, int x1, int y1, unsigned char *im2, int stride2, int x2, int y2) { int i, stride1_i = 0, stride2_i = 0;
diff --git a/av1/encoder/x86/corner_match_sse4.c b/aom_dsp/flow_estimation/x86/corner_match_sse4.c similarity index 96% rename from av1/encoder/x86/corner_match_sse4.c rename to aom_dsp/flow_estimation/x86/corner_match_sse4.c index a718e85..f48201f 100644 --- a/av1/encoder/x86/corner_match_sse4.c +++ b/aom_dsp/flow_estimation/x86/corner_match_sse4.c
@@ -21,7 +21,7 @@ #include "aom_ports/mem.h" #include "aom_ports/system_state.h" -#include "av1/encoder/corner_match.h" +#include "aom_dsp/flow_estimation/corner_match.h" DECLARE_ALIGNED(16, static const uint8_t, byte_mask[16]) = { 255, 255, 255, 255, 255, 255, 255, 255, @@ -34,7 +34,7 @@ correlation/standard deviation are taken over MATCH_SZ by MATCH_SZ windows of each image, centered at (x1, y1) and (x2, y2) respectively. */ -double av1_compute_cross_correlation_sse4_1(unsigned char *im1, int stride1, +double aom_compute_cross_correlation_sse4_1(unsigned char *im1, int stride1, int x1, int y1, unsigned char *im2, int stride2, int x2, int y2) { int i;
diff --git a/av1/encoder/mathutils.h b/aom_dsp/linalg.c similarity index 84% rename from av1/encoder/mathutils.h rename to aom_dsp/linalg.c index 3319a5c..0d97def 100644 --- a/av1/encoder/mathutils.h +++ b/aom_dsp/linalg.c
@@ -10,19 +10,20 @@ * aomedia.org/license/patent-license/. */ -#ifndef AOM_AV1_ENCODER_MATHUTILS_H_ -#define AOM_AV1_ENCODER_MATHUTILS_H_ - #include <memory.h> #include <math.h> #include <stdio.h> #include <stdlib.h> #include <assert.h> +#include "aom_dsp/aom_dsp_common.h" +#include "aom_dsp/linalg.h" +#include "aom_mem/aom_mem.h" + static const double TINY_NEAR_ZERO = 1.0E-16; // Solves Ax = b, where x and b are column vectors of size nx1 and A is nxn -static INLINE int linsolve(int n, double *A, int stride, double *b, double *x) { +int linsolve(int n, double *A, int stride, double *b, double *x) { int i, j, k; double c; // Forward elimination @@ -63,8 +64,8 @@ // Solves for n-dim x in a least squares sense to minimize |Ax - b|^2 // The solution is simply x = (A'A)^-1 A'b or simply the solution for // the system: A'A x = A'b -static INLINE int least_squares(int n, double *A, int rows, int stride, - double *b, double *scratch, double *x) { +int least_squares(int n, double *A, int rows, int stride, double *b, + double *scratch, double *x) { int i, j, k; double *scratch_ = NULL; double *AtA, *Atb; @@ -90,49 +91,9 @@ return ret; } -// Matrix multiply -static INLINE void multiply_mat(const double *m1, const double *m2, double *res, - const int m1_rows, const int inner_dim, - const int m2_cols) { - double sum; - - int row, col, inner; - for (row = 0; row < m1_rows; ++row) { - for (col = 0; col < m2_cols; ++col) { - sum = 0; - for (inner = 0; inner < inner_dim; ++inner) - sum += m1[row * inner_dim + inner] * m2[inner * m2_cols + col]; - *(res++) = sum; - } - } -} - -// -// The functions below are needed only for homography computation -// Remove if the homography models are not used. -// /////////////////////////////////////////////////////////////////////////////// // svdcmp // Adopted from Numerical Recipes in C - -static INLINE double sign(double a, double b) { - return ((b) >= 0 ? fabs(a) : -fabs(a)); -} - -static INLINE double pythag(double a, double b) { - double ct; - const double absa = fabs(a); - const double absb = fabs(b); - - if (absa > absb) { - ct = absb / absa; - return absa * sqrt(1.0 + ct * ct); - } else { - ct = absa / absb; - return (absb == 0) ? 0 : absb * sqrt(1.0 + ct * ct); - } -} - static INLINE int svdcmp(double **u, int m, int n, double w[], double **v) { const int max_its = 30; int flag, i, its, j, jj, k, l, nm; @@ -317,8 +278,7 @@ return 0; } -static INLINE int SVD(double *U, double *W, double *V, double *matx, int M, - int N) { +int SVD(double *U, double *W, double *V, double *matx, int M, int N) { // Assumes allocation for U is MxN double **nrU = (double **)aom_malloc((M) * sizeof(*nrU)); double **nrV = (double **)aom_malloc((N) * sizeof(*nrV)); @@ -356,5 +316,3 @@ return 0; } - -#endif // AOM_AV1_ENCODER_MATHUTILS_H_
diff --git a/aom_dsp/linalg.h b/aom_dsp/linalg.h new file mode 100644 index 0000000..bb38209 --- /dev/null +++ b/aom_dsp/linalg.h
@@ -0,0 +1,76 @@ +/* + * Copyright (c) 2021, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#ifndef AOM_AOM_DSP_LINALG_H_ +#define AOM_AOM_DSP_LINALG_H_ + +#include <math.h> + +#include "config/aom_config.h" + +#ifdef __cplusplus +extern "C" { +#endif // __cplusplus + +// Solves Ax = b, where x and b are column vectors of size nx1 and A is nxn +int linsolve(int n, double *A, int stride, double *b, double *x); + +//////////////////////////////////////////////////////////////////////////////// +// Least-squares +// Solves for n-dim x in a least squares sense to minimize |Ax - b|^2 +// The solution is simply x = (A'A)^-1 A'b or simply the solution for +// the system: A'A x = A'b +int least_squares(int n, double *A, int rows, int stride, double *b, + double *scratch, double *x); + +// Matrix multiply +static INLINE void multiply_mat(const double *m1, const double *m2, double *res, + const int m1_rows, const int inner_dim, + const int m2_cols) { + double sum; + + int row, col, inner; + for (row = 0; row < m1_rows; ++row) { + for (col = 0; col < m2_cols; ++col) { + sum = 0; + for (inner = 0; inner < inner_dim; ++inner) + sum += m1[row * inner_dim + inner] * m2[inner * m2_cols + col]; + *(res++) = sum; + } + } +} + +static INLINE double sign(double a, double b) { + return ((b) >= 0 ? fabs(a) : -fabs(a)); +} + +static INLINE double pythag(double a, double b) { + double ct; + const double absa = fabs(a); + const double absb = fabs(b); + + if (absa > absb) { + ct = absb / absa; + return absa * sqrt(1.0 + ct * ct); + } else { + ct = absa / absb; + return (absb == 0) ? 0 : absb * sqrt(1.0 + ct * ct); + } +} + +int SVD(double *U, double *W, double *V, double *matx, int M, int N); + +#ifdef __cplusplus +} // extern "C" +#endif // __cplusplus + +#endif // AOM_AOM_DSP_LINALG_H_
diff --git a/aom_dsp/noise_model.c b/aom_dsp/noise_model.c index cc47019..1d22ca4 100644 --- a/aom_dsp/noise_model.c +++ b/aom_dsp/noise_model.c
@@ -16,11 +16,11 @@ #include <string.h> #include "aom_dsp/aom_dsp_common.h" +#include "aom_dsp/linalg.h" #include "aom_dsp/noise_model.h" #include "aom_dsp/noise_util.h" #include "aom_mem/aom_mem.h" #include "av1/common/common.h" -#include "av1/encoder/mathutils.h" #define kLowPolyNumParams 3
diff --git a/aom_dsp/rect.h b/aom_dsp/rect.h new file mode 100644 index 0000000..64c7298 --- /dev/null +++ b/aom_dsp/rect.h
@@ -0,0 +1,36 @@ +/* + * Copyright (c) 2022, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#ifndef AOM_AOM_DSP_RECT_H_ +#define AOM_AOM_DSP_RECT_H_ + +#include "config/aom_config.h" + +#include <stdbool.h> + +// Struct representing a rectangle of pixels. +// The axes are inclusive-exclusive, ie. the point (top, left) is included +// in the rectangle but (bottom, right) is not. +typedef struct { + int left, right, top, bottom; +} PixelRect; + +static INLINE int rect_width(const PixelRect *r) { return r->right - r->left; } + +static INLINE int rect_height(const PixelRect *r) { return r->bottom - r->top; } + +static INLINE bool is_inside_rect(const int x, const int y, + const PixelRect *r) { + return (r->left <= x && x < r->right) && (r->top <= y && y < r->bottom); +} + +#endif // AOM_AOM_DSP_RECT_H_
diff --git a/aom_scale/generic/yv12config.c b/aom_scale/generic/yv12config.c index dc3668f..0c1a811 100644 --- a/aom_scale/generic/yv12config.c +++ b/aom_scale/generic/yv12config.c
@@ -13,6 +13,7 @@ #include <assert.h> #include "aom/internal/aom_image_internal.h" +#include "aom_dsp/flow_estimation/pyramid.h" #include "aom_mem/aom_mem.h" #include "aom_ports/mem.h" #include "aom_scale/yv12config.h" @@ -32,7 +33,14 @@ if (ybf->buffer_alloc_sz > 0) { aom_free(ybf->buffer_alloc); } - if (ybf->y_buffer_8bit) aom_free(ybf->y_buffer_8bit); +#if CONFIG_AV1_ENCODER + if (ybf->y_pyramid) { + aom_free_pyramid(ybf->y_pyramid); + } + if (ybf->corners) { + aom_free(ybf->corners); + } +#endif // CONFIG_AV1_ENCODER aom_remove_metadata_from_frame_buffer(ybf); /* buffer_alloc isn't accessed by most functions. Rather y_buffer, u_buffer and v_buffer point to buffer_alloc and are used. Clear out @@ -44,6 +52,25 @@ return AOM_CODEC_MEM_ERROR; } +#if CONFIG_AV1_ENCODER +// Discard global motion data +// This should be called whenever a frame buffer is reused for a new frame, +// to avoid using stale data +void aom_invalidate_gm_data(YV12_BUFFER_CONFIG *ybf) { + // TODO(rachelbarker): Erase data but keep allocations + // This requires appropriate resizing logic + if (ybf->y_pyramid) { + aom_free_pyramid(ybf->y_pyramid); + ybf->y_pyramid = NULL; + } + if (ybf->corners) { + aom_free(ybf->corners); + ybf->corners = NULL; + } + ybf->num_corners = 0; +} +#endif // CONFIG_AV1_ENCODER + static int realloc_frame_buffer_aligned( YV12_BUFFER_CONFIG *ybf, int width, int height, int ss_x, int ss_y, int use_highbitdepth, int border, int byte_alignment, @@ -62,8 +89,7 @@ #if defined AOM_MAX_ALLOCABLE_MEMORY // The size of ybf->buffer_alloc. uint64_t alloc_size = frame_size; - // The size of ybf->y_buffer_8bit. - if (use_highbitdepth) alloc_size += yplane_size; + // TODO(rachelbarker): Add pyramid size here // The decoder may allocate REF_FRAMES frame buffers in the frame buffer // pool. Bound the total amount of allocated memory as if these REF_FRAMES // frame buffers were allocated in a single allocation. @@ -154,19 +180,11 @@ ybf->use_external_reference_buffers = 0; - if (use_highbitdepth) { - if (ybf->y_buffer_8bit) aom_free(ybf->y_buffer_8bit); - ybf->y_buffer_8bit = (uint8_t *)aom_memalign(32, (size_t)yplane_size); - if (!ybf->y_buffer_8bit) return AOM_CODEC_MEM_ERROR; - } else { - if (ybf->y_buffer_8bit) { - aom_free(ybf->y_buffer_8bit); - ybf->y_buffer_8bit = NULL; - } - } - // y_buffer_8bit may have been allocated above, but it has not been filled - // in yet. So, mark it as invalid. - ybf->buf_8bit_valid = 0; +#if CONFIG_AV1_ENCODER + // Discard global motion data, so that stale data is not used + aom_invalidate_gm_data(ybf); +#endif // CONFIG_AV1_ENCODER + ybf->corrupted = 0; /* assume not corrupted by errors */ return 0; }
diff --git a/aom_scale/yv12config.h b/aom_scale/yv12config.h index 6172c77..c271a53 100644 --- a/aom_scale/yv12config.h +++ b/aom_scale/yv12config.h
@@ -24,6 +24,8 @@ #include "aom/aom_integer.h" #include "aom/internal/aom_image_internal.h" +#include "aom_dsp/flow_estimation/pyramid.h" + /*!\cond */ #define AOMINNERBORDERINPIXELS 160 @@ -90,10 +92,12 @@ // external reference frame is no longer used. uint8_t *store_buf_adr[3]; - // If the frame is stored in a 16-bit buffer, this stores an 8-bit version - // for use in global motion detection. It is allocated on-demand. - uint8_t *y_buffer_8bit; - int buf_8bit_valid; +#if CONFIG_AV1_ENCODER + // Data needed for global motion estimation + ImagePyramid *y_pyramid; + int *corners; + int num_corners; +#endif // CONFIG_AV1_ENCODER uint8_t *buffer_alloc; size_t buffer_alloc_sz; @@ -125,6 +129,13 @@ int ss_x, int ss_y, int use_highbitdepth, int border, int byte_alignment); +#if CONFIG_AV1_ENCODER +// Discard global motion data +// This should be called whenever a frame buffer is reused for a new frame, +// to avoid using stale data +void aom_invalidate_gm_data(YV12_BUFFER_CONFIG *ybf); +#endif // CONFIG_AV1_ENCODER + // Updates the yv12 buffer config with the frame buffer. |byte_alignment| must // be a power of 2, from 32 to 1024. 0 sets legacy alignment. If cb is not // NULL, then libaom is using the frame buffer callbacks to handle memory.
diff --git a/av1/av1.cmake b/av1/av1.cmake index 0b94a40..ec17b6f 100644 --- a/av1/av1.cmake +++ b/av1/av1.cmake
@@ -154,10 +154,6 @@ "${AOM_ROOT}/av1/encoder/compound_type.h" "${AOM_ROOT}/av1/encoder/context_tree.c" "${AOM_ROOT}/av1/encoder/context_tree.h" - "${AOM_ROOT}/av1/encoder/corner_detect.c" - "${AOM_ROOT}/av1/encoder/corner_detect.h" - "${AOM_ROOT}/av1/encoder/corner_match.c" - "${AOM_ROOT}/av1/encoder/corner_match.h" "${AOM_ROOT}/av1/encoder/cost.c" "${AOM_ROOT}/av1/encoder/cost.h" "${AOM_ROOT}/av1/encoder/encodeframe.c" @@ -225,8 +221,6 @@ "${AOM_ROOT}/av1/encoder/picklpf.h" "${AOM_ROOT}/av1/encoder/pickrst.c" "${AOM_ROOT}/av1/encoder/pickrst.h" - "${AOM_ROOT}/av1/encoder/ransac.c" - "${AOM_ROOT}/av1/encoder/ransac.h" "${AOM_ROOT}/av1/encoder/ratectrl.c" "${AOM_ROOT}/av1/encoder/ratectrl.h" "${AOM_ROOT}/av1/encoder/rc_utils.h" @@ -369,7 +363,6 @@ "${AOM_ROOT}/av1/encoder/x86/av1_fwd_txfm1d_sse4.c" "${AOM_ROOT}/av1/encoder/x86/av1_fwd_txfm2d_sse4.c" "${AOM_ROOT}/av1/encoder/x86/av1_highbd_quantize_sse4.c" - "${AOM_ROOT}/av1/encoder/x86/corner_match_sse4.c" "${AOM_ROOT}/av1/encoder/x86/encodetxb_sse4.c" "${AOM_ROOT}/av1/encoder/x86/highbd_fwd_txfm_sse4.c" "${AOM_ROOT}/av1/encoder/x86/rdopt_sse4.c" @@ -380,7 +373,6 @@ AOM_AV1_ENCODER_INTRIN_AVX2 "${AOM_ROOT}/av1/encoder/x86/av1_quantize_avx2.c" "${AOM_ROOT}/av1/encoder/x86/av1_highbd_quantize_avx2.c" - "${AOM_ROOT}/av1/encoder/x86/corner_match_avx2.c" "${AOM_ROOT}/av1/encoder/x86/error_intrin_avx2.c" "${AOM_ROOT}/av1/encoder/x86/highbd_block_error_intrin_avx2.c" "${AOM_ROOT}/av1/encoder/x86/av1_fwd_txfm_avx2.h"
diff --git a/av1/common/av1_common_int.h b/av1/common/av1_common_int.h index 0ade4b0..e3f3035 100644 --- a/av1/common/av1_common_int.h +++ b/av1/common/av1_common_int.h
@@ -1625,7 +1625,9 @@ if (new_fb_idx == INVALID_IDX) return NULL; cm->cur_frame = &cm->buffer_pool->frame_bufs[new_fb_idx]; - cm->cur_frame->buf.buf_8bit_valid = 0; +#if CONFIG_AV1_ENCODER + aom_invalidate_gm_data(&cm->cur_frame->buf); +#endif // CONFIG_AV1_ENCODER av1_zero(cm->cur_frame->interp_filter_selected); return cm->cur_frame; }
diff --git a/av1/common/av1_rtcd_defs.pl b/av1/common/av1_rtcd_defs.pl index 4b4d704..5323778 100644 --- a/av1/common/av1_rtcd_defs.pl +++ b/av1/common/av1_rtcd_defs.pl
@@ -506,11 +506,6 @@ add_proto qw/int64_t av1_calc_frame_error/, "const uint8_t *const ref, int stride, const uint8_t *const dst, int p_width, int p_height, int p_stride"; specialize qw/av1_calc_frame_error sse2 avx2/; -if (aom_config("CONFIG_AV1_ENCODER") eq "yes") { - add_proto qw/double av1_compute_cross_correlation/, "unsigned char *im1, int stride1, int x1, int y1, unsigned char *im2, int stride2, int x2, int y2"; - specialize qw/av1_compute_cross_correlation sse4_1 avx2/; -} - # LOOP_RESTORATION functions add_proto qw/void av1_apply_selfguided_restoration/, "const uint8_t *dat, int width, int height, int stride, int eps, const int *xqd, uint8_t *dst, int dst_stride, int32_t *tmpbuf, int bit_depth, int highbd";
diff --git a/av1/common/mv.h b/av1/common/mv.h index a6b197e..4a91cba 100644 --- a/av1/common/mv.h +++ b/av1/common/mv.h
@@ -16,6 +16,7 @@ #include "av1/common/common.h" #include "av1/common/common_data.h" #include "aom_dsp/aom_filter.h" +#include "aom_dsp/flow_estimation/flow_estimation.h" #ifdef __cplusplus extern "C" { @@ -261,39 +262,18 @@ #define WARPEDDIFF_PREC_BITS (WARPEDMODEL_PREC_BITS - WARPEDPIXEL_PREC_BITS) -/* clang-format off */ -enum { - IDENTITY = 0, // identity transformation, 0-parameter - TRANSLATION = 1, // translational motion 2-parameter - ROTZOOM = 2, // simplified affine with rotation + zoom only, 4-parameter - AFFINE = 3, // affine, 6-parameter - TRANS_TYPES, -} UENUM1BYTE(TransformationType); -/* clang-format on */ - -// Number of types used for global motion (must be >= 3 and <= TRANS_TYPES) -// The following can be useful: -// GLOBAL_TRANS_TYPES 3 - up to rotation-zoom -// GLOBAL_TRANS_TYPES 4 - up to affine -// GLOBAL_TRANS_TYPES 6 - up to hor/ver trapezoids -// GLOBAL_TRANS_TYPES 7 - up to full homography -#define GLOBAL_TRANS_TYPES 4 - typedef struct { int global_warp_allowed; int local_warp_allowed; } WarpTypesAllowed; -// number of parameters used by each transformation in TransformationTypes -static const int trans_model_params[TRANS_TYPES] = { 0, 2, 4, 6 }; - // The order of values in the wmmat matrix below is best described // by the homography: // [x' (m2 m3 m0 [x // z . y' = m4 m5 m1 * y // 1] m6 m7 1) 1] typedef struct { - int32_t wmmat[8]; + int32_t wmmat[MAX_PARAMDIM - 1]; int16_t alpha, beta, gamma, delta; TransformationType wmtype; int8_t invalid;
diff --git a/av1/common/restoration.c b/av1/common/restoration.c index e66f9c6..1d2bb78 100644 --- a/av1/common/restoration.c +++ b/av1/common/restoration.c
@@ -40,8 +40,8 @@ { { 2, 0 }, { 56, -1 } }, { { 2, 0 }, { 22, -1 } }, }; -AV1PixelRect av1_whole_frame_rect(const AV1_COMMON *cm, int is_uv) { - AV1PixelRect rect; +PixelRect av1_whole_frame_rect(const AV1_COMMON *cm, int is_uv) { + PixelRect rect; int ss_x = is_uv && cm->seq_params.subsampling_x; int ss_y = is_uv && cm->seq_params.subsampling_y; @@ -71,7 +71,7 @@ // top-left and we can use av1_get_tile_rect(). With CONFIG_MAX_TILE, we have // to do the computation ourselves, iterating over the tiles and keeping // track of the largest width and height, then upscaling. - const AV1PixelRect tile_rect = av1_whole_frame_rect(cm, is_uv); + const PixelRect tile_rect = av1_whole_frame_rect(cm, is_uv); const int max_tile_w = tile_rect.right - tile_rect.left; const int max_tile_h = tile_rect.bottom - tile_rect.top; @@ -244,7 +244,7 @@ // av1_loop_restoration_save_boundary_lines() function, so here we just need // to decide if we're overwriting the above/below boundary pixels or not. static void get_stripe_boundary_info(const RestorationTileLimits *limits, - const AV1PixelRect *tile_rect, int ss_y, + const PixelRect *tile_rect, int ss_y, int *copy_above, int *copy_below) { *copy_above = 1; *copy_below = 1; @@ -1010,7 +1010,7 @@ void av1_loop_restoration_filter_unit( const RestorationTileLimits *limits, const RestorationUnitInfo *rui, const RestorationStripeBoundaries *rsb, RestorationLineBuffers *rlbs, - const AV1PixelRect *tile_rect, int tile_stripe0, int ss_x, int ss_y, + const PixelRect *tile_rect, int tile_stripe0, int ss_x, int ss_y, int highbd, int bit_depth, uint8_t *data8, int stride, uint8_t *dst8, int dst_stride, int32_t *tmpbuf, int optimized_lr) { RestorationType unit_rtype = rui->restoration_type; @@ -1076,8 +1076,8 @@ } static void filter_frame_on_unit(const RestorationTileLimits *limits, - const AV1PixelRect *tile_rect, - int rest_unit_idx, void *priv, int32_t *tmpbuf, + const PixelRect *tile_rect, int rest_unit_idx, + void *priv, int32_t *tmpbuf, RestorationLineBuffers *rlbs) { FilterFrameCtxt *ctxt = (FilterFrameCtxt *)priv; const RestorationInfo *rsi = ctxt->rsi; @@ -1152,7 +1152,7 @@ assert(num_planes <= 3); for (int plane = 0; plane < num_planes; ++plane) { if (cm->rst_info[plane].frame_restoration_type == RESTORE_NONE) continue; - AV1PixelRect tile_rect = loop_rest_ctxt->ctxt[plane].tile_rect; + PixelRect tile_rect = loop_rest_ctxt->ctxt[plane].tile_rect; copy_funs[plane](loop_rest_ctxt->dst, loop_rest_ctxt->frame, tile_rect.left, tile_rect.right, tile_rect.top, tile_rect.bottom); } @@ -1190,7 +1190,7 @@ } void av1_foreach_rest_unit_in_row( - RestorationTileLimits *limits, const AV1PixelRect *tile_rect, + RestorationTileLimits *limits, const PixelRect *tile_rect, rest_unit_visitor_t on_rest_unit, int row_number, int unit_size, int unit_idx0, int hunits_per_tile, int vunits_per_tile, int plane, void *priv, int32_t *tmpbuf, RestorationLineBuffers *rlbs, @@ -1245,7 +1245,7 @@ } static void foreach_rest_unit_in_tile( - const AV1PixelRect *tile_rect, int tile_row, int tile_col, int tile_cols, + const PixelRect *tile_rect, int tile_row, int tile_col, int tile_cols, int hunits_per_tile, int vunits_per_tile, int units_per_tile, int unit_size, int ss_y, int plane, rest_unit_visitor_t on_rest_unit, void *priv, int32_t *tmpbuf, RestorationLineBuffers *rlbs) { @@ -1281,7 +1281,7 @@ void av1_foreach_rest_unit_in_plane(const struct AV1Common *cm, int plane, rest_unit_visitor_t on_rest_unit, - void *priv, AV1PixelRect *tile_rect, + void *priv, PixelRect *tile_rect, int32_t *tmpbuf, RestorationLineBuffers *rlbs) { const int is_uv = plane > 0; @@ -1308,7 +1308,7 @@ const int is_uv = plane > 0; - const AV1PixelRect tile_rect = av1_whole_frame_rect(cm, is_uv); + const PixelRect tile_rect = av1_whole_frame_rect(cm, is_uv); const int tile_w = tile_rect.right - tile_rect.left; const int tile_h = tile_rect.bottom - tile_rect.top; @@ -1486,7 +1486,7 @@ // Get the tile rectangle, with height rounded up to the next multiple of 8 // luma pixels (only relevant for the bottom tile of the frame) - const AV1PixelRect tile_rect = av1_whole_frame_rect(cm, is_uv); + const PixelRect tile_rect = av1_whole_frame_rect(cm, is_uv); const int stripe0 = 0; RestorationStripeBoundaries *boundaries = &cm->rst_info[plane].boundaries;
diff --git a/av1/common/restoration.h b/av1/common/restoration.h index 198ea7f..4dd6ba2 100644 --- a/av1/common/restoration.h +++ b/av1/common/restoration.h
@@ -332,7 +332,7 @@ } RestorationTileLimits; typedef void (*rest_unit_visitor_t)(const RestorationTileLimits *limits, - const AV1PixelRect *tile_rect, + const PixelRect *tile_rect, int rest_unit_idx, void *priv, int32_t *tmpbuf, RestorationLineBuffers *rlbs); @@ -344,7 +344,7 @@ int highbd, bit_depth; uint8_t *data8, *dst8; int data_stride, dst_stride; - AV1PixelRect tile_rect; + PixelRect tile_rect; } FilterFrameCtxt; typedef struct AV1LrStruct { @@ -403,7 +403,7 @@ void av1_loop_restoration_filter_unit( const RestorationTileLimits *limits, const RestorationUnitInfo *rui, const RestorationStripeBoundaries *rsb, RestorationLineBuffers *rlbs, - const AV1PixelRect *tile_rect, int tile_stripe0, int ss_x, int ss_y, + const PixelRect *tile_rect, int tile_stripe0, int ss_x, int ss_y, int highbd, int bit_depth, uint8_t *data8, int stride, uint8_t *dst8, int dst_stride, int32_t *tmpbuf, int optimized_lr); @@ -438,7 +438,7 @@ // Call on_rest_unit for each loop restoration unit in the plane. void av1_foreach_rest_unit_in_plane(const struct AV1Common *cm, int plane, rest_unit_visitor_t on_rest_unit, - void *priv, AV1PixelRect *tile_rect, + void *priv, PixelRect *tile_rect, int32_t *tmpbuf, RestorationLineBuffers *rlbs); @@ -466,13 +466,13 @@ void av1_loop_restoration_copy_planes(AV1LrStruct *loop_rest_ctxt, struct AV1Common *cm, int num_planes); void av1_foreach_rest_unit_in_row( - RestorationTileLimits *limits, const AV1PixelRect *tile_rect, + RestorationTileLimits *limits, const PixelRect *tile_rect, rest_unit_visitor_t on_rest_unit, int row_number, int unit_size, int unit_idx0, int hunits_per_tile, int vunits_per_tile, int plane, void *priv, int32_t *tmpbuf, RestorationLineBuffers *rlbs, sync_read_fn_t on_sync_read, sync_write_fn_t on_sync_write, struct AV1LrSyncData *const lr_sync); -AV1PixelRect av1_whole_frame_rect(const struct AV1Common *cm, int is_uv); +PixelRect av1_whole_frame_rect(const struct AV1Common *cm, int is_uv); int av1_lr_count_units_in_tile(int unit_size, int tile_size); void av1_lr_sync_read_dummy(void *const lr_sync, int r, int c, int plane); void av1_lr_sync_write_dummy(void *const lr_sync, int r, int c,
diff --git a/av1/common/thread_common.c b/av1/common/thread_common.c index 92fdb20..403656c 100644 --- a/av1/common/thread_common.c +++ b/av1/common/thread_common.c
@@ -720,7 +720,7 @@ const int is_uv = plane > 0; const int ss_y = is_uv && cm->seq_params.subsampling_y; - AV1PixelRect tile_rect = ctxt[plane].tile_rect; + PixelRect tile_rect = ctxt[plane].tile_rect; const int unit_size = ctxt[plane].rsi->restoration_unit_size; const int tile_h = tile_rect.bottom - tile_rect.top; @@ -863,7 +863,7 @@ for (int plane = 0; plane < num_planes; plane++) { if (cm->rst_info[plane].frame_restoration_type == RESTORE_NONE) continue; - const AV1PixelRect tile_rect = ctxt[plane].tile_rect; + const PixelRect tile_rect = ctxt[plane].tile_rect; const int max_tile_h = tile_rect.bottom - tile_rect.top; const int unit_size = cm->rst_info[plane].restoration_unit_size;
diff --git a/av1/common/tile_common.c b/av1/common/tile_common.c index e7c6c8f..e7940eb 100644 --- a/av1/common/tile_common.c +++ b/av1/common/tile_common.c
@@ -168,9 +168,9 @@ return sb_cols; } -AV1PixelRect av1_get_tile_rect(const TileInfo *tile_info, const AV1_COMMON *cm, - int is_uv) { - AV1PixelRect r; +PixelRect av1_get_tile_rect(const TileInfo *tile_info, const AV1_COMMON *cm, + int is_uv) { + PixelRect r; // Calculate position in the Y plane r.left = tile_info->mi_col_start * MI_SIZE;
diff --git a/av1/common/tile_common.h b/av1/common/tile_common.h index 5bbea5f..30f657a 100644 --- a/av1/common/tile_common.h +++ b/av1/common/tile_common.h
@@ -18,6 +18,7 @@ #endif #include "config/aom_config.h" +#include "aom_dsp/rect.h" struct AV1Common; struct SequenceHeader; @@ -43,13 +44,9 @@ int av1_get_sb_rows_in_tile(struct AV1Common *cm, TileInfo tile); int av1_get_sb_cols_in_tile(struct AV1Common *cm, TileInfo tile); -typedef struct { - int left, top, right, bottom; -} AV1PixelRect; - // Return the pixel extents of the given tile -AV1PixelRect av1_get_tile_rect(const TileInfo *tile_info, - const struct AV1Common *cm, int is_uv); +PixelRect av1_get_tile_rect(const TileInfo *tile_info, + const struct AV1Common *cm, int is_uv); // Define tile maximum width and area // There is no maximum height since height is limited by area and width limits
diff --git a/av1/common/warped_motion.h b/av1/common/warped_motion.h index 7dd2550..b957151 100644 --- a/av1/common/warped_motion.h +++ b/av1/common/warped_motion.h
@@ -26,7 +26,6 @@ #include "av1/common/mv.h" #include "av1/common/convolve.h" -#define MAX_PARAMDIM 9 #define LEAST_SQUARES_SAMPLES_MAX_BITS 3 #define LEAST_SQUARES_SAMPLES_MAX (1 << LEAST_SQUARES_SAMPLES_MAX_BITS) #define SAMPLES_ARRAY_SIZE (LEAST_SQUARES_SAMPLES_MAX * 2)
diff --git a/av1/encoder/corner_detect.c b/av1/encoder/corner_detect.c deleted file mode 100644 index 1dd7705..0000000 --- a/av1/encoder/corner_detect.c +++ /dev/null
@@ -1,38 +0,0 @@ -/* - * Copyright (c) 2021, Alliance for Open Media. All rights reserved - * - * This source code is subject to the terms of the BSD 3-Clause Clear License - * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear - * License was not distributed with this source code in the LICENSE file, you - * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 - * aomedia.org/license/patent-license/. - */ - -#include <stdlib.h> -#include <stdio.h> -#include <memory.h> -#include <math.h> -#include <assert.h> - -#include "third_party/fastfeat/fast.h" - -#include "av1/encoder/corner_detect.h" - -// Fast_9 wrapper -#define FAST_BARRIER 18 -int av1_fast_corner_detect(unsigned char *buf, int width, int height, - int stride, int *points, int max_points) { - int num_points; - xy *const frm_corners_xy = aom_fast9_detect_nonmax(buf, width, height, stride, - FAST_BARRIER, &num_points); - num_points = (num_points <= max_points ? num_points : max_points); - if (num_points > 0 && frm_corners_xy) { - memcpy(points, frm_corners_xy, sizeof(*frm_corners_xy) * num_points); - free(frm_corners_xy); - return num_points; - } - free(frm_corners_xy); - return 0; -}
diff --git a/av1/encoder/corner_match.c b/av1/encoder/corner_match.c deleted file mode 100644 index 686d8ec..0000000 --- a/av1/encoder/corner_match.c +++ /dev/null
@@ -1,195 +0,0 @@ -/* - * Copyright (c) 2021, Alliance for Open Media. All rights reserved - * - * This source code is subject to the terms of the BSD 3-Clause Clear License - * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear - * License was not distributed with this source code in the LICENSE file, you - * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 - * aomedia.org/license/patent-license/. - */ - -#include <stdlib.h> -#include <memory.h> -#include <math.h> - -#include "config/av1_rtcd.h" - -#include "aom_ports/system_state.h" -#include "av1/encoder/corner_match.h" - -#define SEARCH_SZ 9 -#define SEARCH_SZ_BY2 ((SEARCH_SZ - 1) / 2) - -#define THRESHOLD_NCC 0.75 - -/* Compute var(im) * MATCH_SZ_SQ over a MATCH_SZ by MATCH_SZ window of im, - centered at (x, y). -*/ -static double compute_variance(unsigned char *im, int stride, int x, int y) { - int sum = 0; - int sumsq = 0; - int var; - int i, j; - for (i = 0; i < MATCH_SZ; ++i) - for (j = 0; j < MATCH_SZ; ++j) { - sum += im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)]; - sumsq += im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)] * - im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)]; - } - var = sumsq * MATCH_SZ_SQ - sum * sum; - return (double)var; -} - -/* Compute corr(im1, im2) * MATCH_SZ * stddev(im1), where the - correlation/standard deviation are taken over MATCH_SZ by MATCH_SZ windows - of each image, centered at (x1, y1) and (x2, y2) respectively. -*/ -double av1_compute_cross_correlation_c(unsigned char *im1, int stride1, int x1, - int y1, unsigned char *im2, int stride2, - int x2, int y2) { - int v1, v2; - int sum1 = 0; - int sum2 = 0; - int sumsq2 = 0; - int cross = 0; - int var2, cov; - int i, j; - for (i = 0; i < MATCH_SZ; ++i) - for (j = 0; j < MATCH_SZ; ++j) { - v1 = im1[(i + y1 - MATCH_SZ_BY2) * stride1 + (j + x1 - MATCH_SZ_BY2)]; - v2 = im2[(i + y2 - MATCH_SZ_BY2) * stride2 + (j + x2 - MATCH_SZ_BY2)]; - sum1 += v1; - sum2 += v2; - sumsq2 += v2 * v2; - cross += v1 * v2; - } - var2 = sumsq2 * MATCH_SZ_SQ - sum2 * sum2; - cov = cross * MATCH_SZ_SQ - sum1 * sum2; - aom_clear_system_state(); - return cov / sqrt((double)var2); -} - -static int is_eligible_point(int pointx, int pointy, int width, int height) { - return (pointx >= MATCH_SZ_BY2 && pointy >= MATCH_SZ_BY2 && - pointx + MATCH_SZ_BY2 < width && pointy + MATCH_SZ_BY2 < height); -} - -static int is_eligible_distance(int point1x, int point1y, int point2x, - int point2y, int width, int height) { - const int thresh = (width < height ? height : width) >> 4; - return ((point1x - point2x) * (point1x - point2x) + - (point1y - point2y) * (point1y - point2y)) <= thresh * thresh; -} - -static void improve_correspondence(unsigned char *frm, unsigned char *ref, - int width, int height, int frm_stride, - int ref_stride, - Correspondence *correspondences, - int num_correspondences) { - int i; - for (i = 0; i < num_correspondences; ++i) { - int x, y, best_x = 0, best_y = 0; - double best_match_ncc = 0.0; - for (y = -SEARCH_SZ_BY2; y <= SEARCH_SZ_BY2; ++y) { - for (x = -SEARCH_SZ_BY2; x <= SEARCH_SZ_BY2; ++x) { - double match_ncc; - if (!is_eligible_point(correspondences[i].rx + x, - correspondences[i].ry + y, width, height)) - continue; - if (!is_eligible_distance(correspondences[i].x, correspondences[i].y, - correspondences[i].rx + x, - correspondences[i].ry + y, width, height)) - continue; - match_ncc = av1_compute_cross_correlation( - frm, frm_stride, correspondences[i].x, correspondences[i].y, ref, - ref_stride, correspondences[i].rx + x, correspondences[i].ry + y); - if (match_ncc > best_match_ncc) { - best_match_ncc = match_ncc; - best_y = y; - best_x = x; - } - } - } - correspondences[i].rx += best_x; - correspondences[i].ry += best_y; - } - for (i = 0; i < num_correspondences; ++i) { - int x, y, best_x = 0, best_y = 0; - double best_match_ncc = 0.0; - for (y = -SEARCH_SZ_BY2; y <= SEARCH_SZ_BY2; ++y) - for (x = -SEARCH_SZ_BY2; x <= SEARCH_SZ_BY2; ++x) { - double match_ncc; - if (!is_eligible_point(correspondences[i].x + x, - correspondences[i].y + y, width, height)) - continue; - if (!is_eligible_distance( - correspondences[i].x + x, correspondences[i].y + y, - correspondences[i].rx, correspondences[i].ry, width, height)) - continue; - match_ncc = av1_compute_cross_correlation( - ref, ref_stride, correspondences[i].rx, correspondences[i].ry, frm, - frm_stride, correspondences[i].x + x, correspondences[i].y + y); - if (match_ncc > best_match_ncc) { - best_match_ncc = match_ncc; - best_y = y; - best_x = x; - } - } - correspondences[i].x += best_x; - correspondences[i].y += best_y; - } -} - -int av1_determine_correspondence(unsigned char *src, int *src_corners, - int num_src_corners, unsigned char *ref, - int *ref_corners, int num_ref_corners, - int width, int height, int src_stride, - int ref_stride, int *correspondence_pts) { - // TODO(sarahparker) Improve this to include 2-way match - int i, j; - Correspondence *correspondences = (Correspondence *)correspondence_pts; - int num_correspondences = 0; - for (i = 0; i < num_src_corners; ++i) { - double best_match_ncc = 0.0; - double template_norm; - int best_match_j = -1; - if (!is_eligible_point(src_corners[2 * i], src_corners[2 * i + 1], width, - height)) - continue; - for (j = 0; j < num_ref_corners; ++j) { - double match_ncc; - if (!is_eligible_point(ref_corners[2 * j], ref_corners[2 * j + 1], width, - height)) - continue; - if (!is_eligible_distance(src_corners[2 * i], src_corners[2 * i + 1], - ref_corners[2 * j], ref_corners[2 * j + 1], - width, height)) - continue; - match_ncc = av1_compute_cross_correlation( - src, src_stride, src_corners[2 * i], src_corners[2 * i + 1], ref, - ref_stride, ref_corners[2 * j], ref_corners[2 * j + 1]); - if (match_ncc > best_match_ncc) { - best_match_ncc = match_ncc; - best_match_j = j; - } - } - // Note: We want to test if the best correlation is >= THRESHOLD_NCC, - // but need to account for the normalization in - // av1_compute_cross_correlation. - template_norm = compute_variance(src, src_stride, src_corners[2 * i], - src_corners[2 * i + 1]); - if (best_match_ncc > THRESHOLD_NCC * sqrt(template_norm)) { - correspondences[num_correspondences].x = src_corners[2 * i]; - correspondences[num_correspondences].y = src_corners[2 * i + 1]; - correspondences[num_correspondences].rx = ref_corners[2 * best_match_j]; - correspondences[num_correspondences].ry = - ref_corners[2 * best_match_j + 1]; - num_correspondences++; - } - } - improve_correspondence(src, ref, width, height, src_stride, ref_stride, - correspondences, num_correspondences); - return num_correspondences; -}
diff --git a/av1/encoder/corner_match.h b/av1/encoder/corner_match.h deleted file mode 100644 index 3392630..0000000 --- a/av1/encoder/corner_match.h +++ /dev/null
@@ -1,34 +0,0 @@ -/* - * Copyright (c) 2021, Alliance for Open Media. All rights reserved - * - * This source code is subject to the terms of the BSD 3-Clause Clear License - * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear - * License was not distributed with this source code in the LICENSE file, you - * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 - * aomedia.org/license/patent-license/. - */ -#ifndef AOM_AV1_ENCODER_CORNER_MATCH_H_ -#define AOM_AV1_ENCODER_CORNER_MATCH_H_ - -#include <stdio.h> -#include <stdlib.h> -#include <memory.h> - -#define MATCH_SZ 13 -#define MATCH_SZ_BY2 ((MATCH_SZ - 1) / 2) -#define MATCH_SZ_SQ (MATCH_SZ * MATCH_SZ) - -typedef struct { - int x, y; - int rx, ry; -} Correspondence; - -int av1_determine_correspondence(unsigned char *src, int *src_corners, - int num_src_corners, unsigned char *ref, - int *ref_corners, int num_ref_corners, - int width, int height, int src_stride, - int ref_stride, int *correspondence_pts); - -#endif // AOM_AV1_ENCODER_CORNER_MATCH_H_
diff --git a/av1/encoder/encoder.c b/av1/encoder/encoder.c index 9936381..aaf765c 100644 --- a/av1/encoder/encoder.c +++ b/av1/encoder/encoder.c
@@ -61,6 +61,7 @@ #include "av1/encoder/encodetxb.h" #include "av1/encoder/ethread.h" #include "av1/encoder/firstpass.h" +#include "av1/encoder/global_motion_facade.h" #include "av1/encoder/hash_motion.h" #include "av1/encoder/intra_mode_search.h" #include "av1/encoder/mv_prec.h" @@ -2420,7 +2421,7 @@ int phase_scaler = 0; set_size_independent_vars(cpi); - cpi->source->buf_8bit_valid = 0; + aom_invalidate_gm_data(cpi->source); av1_setup_frame_size(cpi); av1_set_size_dependent_vars(cpi, &q, &bottom_index, &top_index); @@ -2595,7 +2596,7 @@ assert(IMPLIES(oxcf->rc_cfg.min_cr > 0, allow_recode)); set_size_independent_vars(cpi); - cpi->source->buf_8bit_valid = 0; + aom_invalidate_gm_data(cpi->source); av1_setup_frame_size(cpi); @@ -3527,6 +3528,37 @@ // for the purpose to verify no mismatch between encoder and decoder. if (cm->show_frame) cpi->last_show_frame_buf = cm->cur_frame; + av1_free_flow_fields(cpi); + +#if CONFIG_GM_USE_SRC_FRAMES + // Before storing the reconstructed frame (cm->cur_frame) into the reference + // buffers, transfer the original frame's (cpi->source's) global motion + // information to it. + // + // This means that, when this frame is used as a reference, all of the + // global motion estimation functions will use the pyramid and corner + // list which were constructed from the original frame, not from + // the reconstructed frame. + // + // Note: key/intra frames will not have computed the pyramid yet, + // so we need to do that before we discard the source frame. + // We do not need to compute the corner list, as this will be derived + // from the pyramid when needed. + if (!cpi->source->y_pyramid) { + cpi->source->y_pyramid = aom_compute_pyramid( + cpi->source, cm->seq_params.bit_depth, MAX_PYRAMID_LEVELS); + assert(cpi->source->y_pyramid); + } + + cm->cur_frame->buf.y_pyramid = cpi->source->y_pyramid; + cm->cur_frame->buf.corners = cpi->source->corners; + cm->cur_frame->buf.num_corners = cpi->source->num_corners; + + cpi->source->y_pyramid = NULL; + cpi->source->corners = NULL; + cpi->source->num_corners = 0; +#endif // CONFIG_GM_USE_SRC_FRAMES + refresh_reference_frames(cpi); #if CONFIG_ENTROPY_STATS
diff --git a/av1/encoder/encoder.h b/av1/encoder/encoder.h index 05772c7..fc12609 100644 --- a/av1/encoder/encoder.h +++ b/av1/encoder/encoder.h
@@ -61,6 +61,7 @@ #endif #include "aom/internal/aom_codec_internal.h" +#include "aom_dsp/flow_estimation/flow_estimation.h" #include "aom_util/aom_thread.h" #ifdef __cplusplus @@ -1833,12 +1834,6 @@ */ typedef struct { /*! - * Array to store the cost for signalling each global motion model. - * gmtype_cost[i] stores the cost of signalling the ith Global Motion model. - */ - int type_cost[TRANS_TYPES]; - - /*! * Array to store the cost for signalling a particular global motion model for * each reference frame. gmparams_cost[i] stores the cost of signalling global * motion for the ith reference frame. @@ -1866,11 +1861,6 @@ #endif // CONFIG_NEW_REF_SIGNALING /*! - * Pointer to the source frame buffer. - */ - unsigned char *src_buffer; - - /*! * Holds the number of valid reference frames in past and future directions * w.r.t. the current frame. num_ref_frames[i] stores the total number of * valid reference frames in 'i' direction. @@ -1889,6 +1879,13 @@ FrameDistPair reference_frames[MAX_DIRECTIONS][REF_FRAMES - 1]; #endif // CONFIG_NEW_REF_SIGNALING + /*! + * Array of structures which hold flow information per ref frame + * Exactly what information is stored depends on the motion estimation + * method - see aom_dsp/flow_estimation/flow_estimation.h for details + */ + FlowData *flow_data[REF_FRAMES]; + /** * \name Dimensions for which segment map is allocated. */ @@ -1896,18 +1893,6 @@ int segment_map_w; /*!< segment map width */ int segment_map_h; /*!< segment map height */ /**@}*/ - - /*! - * Holds the total number of corner points detected in the source frame. - */ - int num_src_corners; - - /*! - * Holds the x and y co-ordinates of the corner points detected in the source - * frame. src_corners[i] holds the x co-ordinate and src_corners[i+1] holds - * the y co-ordinate of the ith corner point detected. - */ - int src_corners[2 * MAX_CORNERS]; } GlobalMotionInfo; /*!
diff --git a/av1/encoder/ethread.c b/av1/encoder/ethread.c index 2fc2b17..af9b1db 100644 --- a/av1/encoder/ethread.c +++ b/av1/encoder/ethread.c
@@ -1391,10 +1391,9 @@ // Compute global motion for the given ref_buf_idx. av1_compute_gm_for_valid_ref_frames( - cpi, gm_info->ref_buf, ref_buf_idx, gm_info->num_src_corners, - gm_info->src_corners, gm_info->src_buffer, - gm_thread_data->params_by_motion, gm_thread_data->segment_map, - gm_info->segment_map_w, gm_info->segment_map_h); + cpi, gm_info->ref_buf, ref_buf_idx, gm_thread_data->params_by_motion, + gm_thread_data->segment_map, gm_info->segment_map_w, + gm_info->segment_map_h); #if CONFIG_MULTITHREAD pthread_mutex_lock(gm_mt_mutex_);
diff --git a/av1/encoder/global_motion.c b/av1/encoder/global_motion.c index 66f9bef..e442802 100644 --- a/av1/encoder/global_motion.c +++ b/av1/encoder/global_motion.c
@@ -21,50 +21,15 @@ #include "av1/encoder/global_motion.h" #include "av1/common/convolve.h" -#include "av1/common/resize.h" #include "av1/common/warped_motion.h" #include "av1/encoder/segmentation.h" -#include "av1/encoder/corner_detect.h" -#include "av1/encoder/corner_match.h" -#include "av1/encoder/ransac.h" - -#define MIN_INLIER_PROB 0.1 #define MIN_TRANS_THRESH (1 * GM_TRANS_DECODE_FACTOR) // Border over which to compute the global motion #define ERRORADV_BORDER 0 -// Number of pyramid levels in disflow computation -#define N_LEVELS 2 -// Size of square patches in the disflow dense grid -#define PATCH_SIZE 8 -// Center point of square patch -#define PATCH_CENTER ((PATCH_SIZE + 1) >> 1) -// Step size between patches, lower value means greater patch overlap -#define PATCH_STEP 1 -// Minimum size of border padding for disflow -#define MIN_PAD 7 -// Warp error convergence threshold for disflow -#define DISFLOW_ERROR_TR 0.01 -// Max number of iterations if warp convergence is not found -#define DISFLOW_MAX_ITR 10 - -// Struct for an image pyramid -typedef struct { - int n_levels; - int pad_size; - int has_gradient; - int widths[N_LEVELS]; - int heights[N_LEVELS]; - int strides[N_LEVELS]; - int level_loc[N_LEVELS]; - unsigned char *level_buffer; - double *level_dx_buffer; - double *level_dy_buffer; -} ImagePyramid; - int av1_is_enough_erroradvantage(double best_erroradvantage, int params_cost) { return best_erroradvantage < erroradv_tr && best_erroradvantage * params_cost < erroradv_prod_tr; @@ -270,10 +235,9 @@ int d_width, int d_height, int d_stride, int n_refinements, int64_t best_frame_error, uint8_t *segment_map, int segment_map_stride, int64_t erroradv_threshold) { - static const int max_trans_model_params[TRANS_TYPES] = { 0, 2, 4, 6 }; const int border = ERRORADV_BORDER; int i = 0, p; - int n_params = max_trans_model_params[wmtype]; + int n_params = trans_model_params[wmtype]; int32_t *param_mat = wm->wmmat; int64_t step_error, best_error; int32_t step; @@ -352,48 +316,6 @@ return best_error; } -unsigned char *av1_downconvert_frame(YV12_BUFFER_CONFIG *frm, int bit_depth) { - int i, j; - uint16_t *orig_buf = CONVERT_TO_SHORTPTR(frm->y_buffer); - uint8_t *buf_8bit = frm->y_buffer_8bit; - assert(buf_8bit); - if (!frm->buf_8bit_valid) { - for (i = 0; i < frm->y_height; ++i) { - for (j = 0; j < frm->y_width; ++j) { - buf_8bit[i * frm->y_stride + j] = - orig_buf[i * frm->y_stride + j] >> (bit_depth - 8); - } - } - frm->buf_8bit_valid = 1; - } -#if CONFIG_DEBUG - else { - // frm->buf_8bit_valid == 1. So, double check that 'buf_8bit' is correct. - for (i = 0; i < frm->y_height; ++i) { - for (j = 0; j < frm->y_width; ++j) { - assert(buf_8bit[i * frm->y_stride + j] == - (orig_buf[i * frm->y_stride + j] >> (bit_depth - 8))); - } - } - } -#endif // CONFIG_DEBUG - return buf_8bit; -} - -static void get_inliers_from_indices(MotionModel *params, - int *correspondences) { - int *inliers_tmp = (int *)aom_malloc(2 * MAX_CORNERS * sizeof(*inliers_tmp)); - memset(inliers_tmp, 0, 2 * MAX_CORNERS * sizeof(*inliers_tmp)); - - for (int i = 0; i < params->num_inliers; i++) { - int index = params->inliers[i]; - inliers_tmp[2 * i] = correspondences[4 * index]; - inliers_tmp[2 * i + 1] = correspondences[4 * index + 1]; - } - memcpy(params->inliers, inliers_tmp, sizeof(*inliers_tmp) * 2 * MAX_CORNERS); - aom_free(inliers_tmp); -} - #define FEAT_COUNT_TR 3 #define SEG_COUNT_TR 0.40 void av1_compute_feature_segmentation_map(uint8_t *segment_map, int width, @@ -423,596 +345,3 @@ if (seg_count < (width * height * SEG_COUNT_TR)) memset(segment_map, 1, width * height * sizeof(*segment_map)); } - -static int compute_global_motion_feature_based( - TransformationType type, unsigned char *src_buffer, int src_width, - int src_height, int src_stride, int *src_corners, int num_src_corners, - YV12_BUFFER_CONFIG *ref, int bit_depth, int *num_inliers_by_motion, - MotionModel *params_by_motion, int num_motions) { - int i; - int num_ref_corners; - int num_correspondences; - int *correspondences; - int ref_corners[2 * MAX_CORNERS]; - unsigned char *ref_buffer = ref->y_buffer; - RansacFunc ransac = av1_get_ransac_type(type); - - if (ref->flags & YV12_FLAG_HIGHBITDEPTH) { - ref_buffer = av1_downconvert_frame(ref, bit_depth); - } - - num_ref_corners = - av1_fast_corner_detect(ref_buffer, ref->y_width, ref->y_height, - ref->y_stride, ref_corners, MAX_CORNERS); - - // find correspondences between the two images - correspondences = - (int *)malloc(num_src_corners * 4 * sizeof(*correspondences)); - num_correspondences = av1_determine_correspondence( - src_buffer, (int *)src_corners, num_src_corners, ref_buffer, - (int *)ref_corners, num_ref_corners, src_width, src_height, src_stride, - ref->y_stride, correspondences); - - ransac(correspondences, num_correspondences, num_inliers_by_motion, - params_by_motion, num_motions); - - // Set num_inliers = 0 for motions with too few inliers so they are ignored. - for (i = 0; i < num_motions; ++i) { - if (num_inliers_by_motion[i] < MIN_INLIER_PROB * num_correspondences || - num_correspondences == 0) { - num_inliers_by_motion[i] = 0; - } else { - get_inliers_from_indices(¶ms_by_motion[i], correspondences); - } - } - - free(correspondences); - - // Return true if any one of the motions has inliers. - for (i = 0; i < num_motions; ++i) { - if (num_inliers_by_motion[i] > 0) return 1; - } - return 0; -} - -// Don't use points around the frame border since they are less reliable -static INLINE int valid_point(int x, int y, int width, int height) { - return (x > (PATCH_SIZE + PATCH_CENTER)) && - (x < (width - PATCH_SIZE - PATCH_CENTER)) && - (y > (PATCH_SIZE + PATCH_CENTER)) && - (y < (height - PATCH_SIZE - PATCH_CENTER)); -} - -static int determine_disflow_correspondence(int *frm_corners, - int num_frm_corners, double *flow_u, - double *flow_v, int width, - int height, int stride, - double *correspondences) { - int num_correspondences = 0; - int x, y; - for (int i = 0; i < num_frm_corners; ++i) { - x = frm_corners[2 * i]; - y = frm_corners[2 * i + 1]; - if (valid_point(x, y, width, height)) { - correspondences[4 * num_correspondences] = x; - correspondences[4 * num_correspondences + 1] = y; - correspondences[4 * num_correspondences + 2] = x + flow_u[y * stride + x]; - correspondences[4 * num_correspondences + 3] = y + flow_v[y * stride + x]; - num_correspondences++; - } - } - return num_correspondences; -} - -static double getCubicValue(double p[4], double x) { - return p[1] + 0.5 * x * - (p[2] - p[0] + - x * (2.0 * p[0] - 5.0 * p[1] + 4.0 * p[2] - p[3] + - x * (3.0 * (p[1] - p[2]) + p[3] - p[0]))); -} - -static void get_subcolumn(unsigned char *ref, double col[4], int stride, int x, - int y_start) { - int i; - for (i = 0; i < 4; ++i) { - col[i] = ref[(i + y_start) * stride + x]; - } -} - -static double bicubic(unsigned char *ref, double x, double y, int stride) { - double arr[4]; - int k; - int i = (int)x; - int j = (int)y; - for (k = 0; k < 4; ++k) { - double arr_temp[4]; - get_subcolumn(ref, arr_temp, stride, i + k - 1, j - 1); - arr[k] = getCubicValue(arr_temp, y - j); - } - return getCubicValue(arr, x - i); -} - -// Interpolate a warped block using bicubic interpolation when possible -static unsigned char interpolate(unsigned char *ref, double x, double y, - int width, int height, int stride) { - if (x < 0 && y < 0) - return ref[0]; - else if (x < 0 && y > height - 1) - return ref[(height - 1) * stride]; - else if (x > width - 1 && y < 0) - return ref[width - 1]; - else if (x > width - 1 && y > height - 1) - return ref[(height - 1) * stride + (width - 1)]; - else if (x < 0) { - int v; - int i = (int)y; - double a = y - i; - if (y > 1 && y < height - 2) { - double arr[4]; - get_subcolumn(ref, arr, stride, 0, i - 1); - return clamp((int)(getCubicValue(arr, a) + 0.5), 0, 255); - } - v = (int)(ref[i * stride] * (1 - a) + ref[(i + 1) * stride] * a + 0.5); - return clamp(v, 0, 255); - } else if (y < 0) { - int v; - int j = (int)x; - double b = x - j; - if (x > 1 && x < width - 2) { - double arr[4] = { ref[j - 1], ref[j], ref[j + 1], ref[j + 2] }; - return clamp((int)(getCubicValue(arr, b) + 0.5), 0, 255); - } - v = (int)(ref[j] * (1 - b) + ref[j + 1] * b + 0.5); - return clamp(v, 0, 255); - } else if (x > width - 1) { - int v; - int i = (int)y; - double a = y - i; - if (y > 1 && y < height - 2) { - double arr[4]; - get_subcolumn(ref, arr, stride, width - 1, i - 1); - return clamp((int)(getCubicValue(arr, a) + 0.5), 0, 255); - } - v = (int)(ref[i * stride + width - 1] * (1 - a) + - ref[(i + 1) * stride + width - 1] * a + 0.5); - return clamp(v, 0, 255); - } else if (y > height - 1) { - int v; - int j = (int)x; - double b = x - j; - if (x > 1 && x < width - 2) { - int row = (height - 1) * stride; - double arr[4] = { ref[row + j - 1], ref[row + j], ref[row + j + 1], - ref[row + j + 2] }; - return clamp((int)(getCubicValue(arr, b) + 0.5), 0, 255); - } - v = (int)(ref[(height - 1) * stride + j] * (1 - b) + - ref[(height - 1) * stride + j + 1] * b + 0.5); - return clamp(v, 0, 255); - } else if (x > 1 && y > 1 && x < width - 2 && y < height - 2) { - return clamp((int)(bicubic(ref, x, y, stride) + 0.5), 0, 255); - } else { - int i = (int)y; - int j = (int)x; - double a = y - i; - double b = x - j; - int v = (int)(ref[i * stride + j] * (1 - a) * (1 - b) + - ref[i * stride + j + 1] * (1 - a) * b + - ref[(i + 1) * stride + j] * a * (1 - b) + - ref[(i + 1) * stride + j + 1] * a * b); - return clamp(v, 0, 255); - } -} - -// Warps a block using flow vector [u, v] and computes the mse -static double compute_warp_and_error(unsigned char *ref, unsigned char *frm, - int width, int height, int stride, int x, - int y, double u, double v, int16_t *dt) { - int i, j; - unsigned char warped; - double x_w, y_w; - double mse = 0; - int16_t err = 0; - for (i = y; i < y + PATCH_SIZE; ++i) - for (j = x; j < x + PATCH_SIZE; ++j) { - x_w = (double)j + u; - y_w = (double)i + v; - warped = interpolate(ref, x_w, y_w, width, height, stride); - err = warped - frm[j + i * stride]; - mse += err * err; - dt[(i - y) * PATCH_SIZE + (j - x)] = err; - } - - mse /= (PATCH_SIZE * PATCH_SIZE); - return mse; -} - -// Computes the components of the system of equations used to solve for -// a flow vector. This includes: -// 1.) The hessian matrix for optical flow. This matrix is in the -// form of: -// -// M = |sum(dx * dx) sum(dx * dy)| -// |sum(dx * dy) sum(dy * dy)| -// -// 2.) b = |sum(dx * dt)| -// |sum(dy * dt)| -// Where the sums are computed over a square window of PATCH_SIZE. -static INLINE void compute_flow_system(const double *dx, int dx_stride, - const double *dy, int dy_stride, - const int16_t *dt, int dt_stride, - double *M, double *b) { - for (int i = 0; i < PATCH_SIZE; i++) { - for (int j = 0; j < PATCH_SIZE; j++) { - M[0] += dx[i * dx_stride + j] * dx[i * dx_stride + j]; - M[1] += dx[i * dx_stride + j] * dy[i * dy_stride + j]; - M[3] += dy[i * dy_stride + j] * dy[i * dy_stride + j]; - - b[0] += dx[i * dx_stride + j] * dt[i * dt_stride + j]; - b[1] += dy[i * dy_stride + j] * dt[i * dt_stride + j]; - } - } - - M[2] = M[1]; -} - -// Solves a general Mx = b where M is a 2x2 matrix and b is a 2x1 matrix -static INLINE void solve_2x2_system(const double *M, const double *b, - double *output_vec) { - double M_0 = M[0]; - double M_3 = M[3]; - double det = (M_0 * M_3) - (M[1] * M[2]); - if (det < 1e-5) { - // Handle singular matrix - // TODO(sarahparker) compare results using pseudo inverse instead - M_0 += 1e-10; - M_3 += 1e-10; - det = (M_0 * M_3) - (M[1] * M[2]); - } - const double det_inv = 1 / det; - const double mult_b0 = det_inv * b[0]; - const double mult_b1 = det_inv * b[1]; - output_vec[0] = M_3 * mult_b0 - M[1] * mult_b1; - output_vec[1] = -M[2] * mult_b0 + M_0 * mult_b1; -} - -/* -static INLINE void image_difference(const uint8_t *src, int src_stride, - const uint8_t *ref, int ref_stride, - int16_t *dst, int dst_stride, int height, - int width) { - const int block_unit = 8; - // Take difference in 8x8 blocks to make use of optimized diff function - for (int i = 0; i < height; i += block_unit) { - for (int j = 0; j < width; j += block_unit) { - aom_subtract_block(block_unit, block_unit, dst + i * dst_stride + j, - dst_stride, src + i * src_stride + j, src_stride, - ref + i * ref_stride + j, ref_stride); - } - } -} -*/ - -// Compute an image gradient using a sobel filter. -// If dir == 1, compute the x gradient. If dir == 0, compute y. This function -// assumes the images have been padded so that they can be processed in units -// of 8. -static INLINE void sobel_xy_image_gradient(const uint8_t *src, int src_stride, - double *dst, int dst_stride, - int height, int width, int dir) { - double norm = 1.0; - // TODO(sarahparker) experiment with doing this over larger block sizes - const int block_unit = 8; - // Filter in 8x8 blocks to eventually make use of optimized convolve function - for (int i = 0; i < height; i += block_unit) { - for (int j = 0; j < width; j += block_unit) { - av1_convolve_2d_sobel_y_c(src + i * src_stride + j, src_stride, - dst + i * dst_stride + j, dst_stride, - block_unit, block_unit, dir, norm); - } - } -} - -static ImagePyramid *alloc_pyramid(int width, int height, int pad_size, - int compute_gradient) { - ImagePyramid *pyr = aom_malloc(sizeof(*pyr)); - pyr->has_gradient = compute_gradient; - // 2 * width * height is the upper bound for a buffer that fits - // all pyramid levels + padding for each level - const int buffer_size = sizeof(*pyr->level_buffer) * 2 * width * height + - (width + 2 * pad_size) * 2 * pad_size * N_LEVELS; - pyr->level_buffer = aom_malloc(buffer_size); - memset(pyr->level_buffer, 0, buffer_size); - - if (compute_gradient) { - const int gradient_size = - sizeof(*pyr->level_dx_buffer) * 2 * width * height + - (width + 2 * pad_size) * 2 * pad_size * N_LEVELS; - pyr->level_dx_buffer = aom_malloc(gradient_size); - pyr->level_dy_buffer = aom_malloc(gradient_size); - memset(pyr->level_dx_buffer, 0, gradient_size); - memset(pyr->level_dy_buffer, 0, gradient_size); - } - return pyr; -} - -static void free_pyramid(ImagePyramid *pyr) { - aom_free(pyr->level_buffer); - if (pyr->has_gradient) { - aom_free(pyr->level_dx_buffer); - aom_free(pyr->level_dy_buffer); - } - aom_free(pyr); -} - -static INLINE void update_level_dims(ImagePyramid *frm_pyr, int level) { - frm_pyr->widths[level] = frm_pyr->widths[level - 1] >> 1; - frm_pyr->heights[level] = frm_pyr->heights[level - 1] >> 1; - frm_pyr->strides[level] = frm_pyr->widths[level] + 2 * frm_pyr->pad_size; - // Point the beginning of the next level buffer to the correct location inside - // the padded border - frm_pyr->level_loc[level] = - frm_pyr->level_loc[level - 1] + - frm_pyr->strides[level - 1] * - (2 * frm_pyr->pad_size + frm_pyr->heights[level - 1]); -} - -// Compute coarse to fine pyramids for a frame -static void compute_flow_pyramids(unsigned char *frm, const int frm_width, - const int frm_height, const int frm_stride, - int n_levels, int pad_size, int compute_grad, - ImagePyramid *frm_pyr) { - int cur_width, cur_height, cur_stride, cur_loc; - assert((frm_width >> n_levels) > 0); - assert((frm_height >> n_levels) > 0); - - // Initialize first level - frm_pyr->n_levels = n_levels; - frm_pyr->pad_size = pad_size; - frm_pyr->widths[0] = frm_width; - frm_pyr->heights[0] = frm_height; - frm_pyr->strides[0] = frm_width + 2 * frm_pyr->pad_size; - // Point the beginning of the level buffer to the location inside - // the padded border - frm_pyr->level_loc[0] = - frm_pyr->strides[0] * frm_pyr->pad_size + frm_pyr->pad_size; - // This essentially copies the original buffer into the pyramid buffer - // without the original padding - av1_resize_plane(frm, frm_height, frm_width, frm_stride, - frm_pyr->level_buffer + frm_pyr->level_loc[0], - frm_pyr->heights[0], frm_pyr->widths[0], - frm_pyr->strides[0]); - - if (compute_grad) { - cur_width = frm_pyr->widths[0]; - cur_height = frm_pyr->heights[0]; - cur_stride = frm_pyr->strides[0]; - cur_loc = frm_pyr->level_loc[0]; - assert(frm_pyr->has_gradient && frm_pyr->level_dx_buffer != NULL && - frm_pyr->level_dy_buffer != NULL); - // Computation x gradient - sobel_xy_image_gradient(frm_pyr->level_buffer + cur_loc, cur_stride, - frm_pyr->level_dx_buffer + cur_loc, cur_stride, - cur_height, cur_width, 1); - - // Computation y gradient - sobel_xy_image_gradient(frm_pyr->level_buffer + cur_loc, cur_stride, - frm_pyr->level_dy_buffer + cur_loc, cur_stride, - cur_height, cur_width, 0); - } - - // Start at the finest level and resize down to the coarsest level - for (int level = 1; level < n_levels; ++level) { - update_level_dims(frm_pyr, level); - cur_width = frm_pyr->widths[level]; - cur_height = frm_pyr->heights[level]; - cur_stride = frm_pyr->strides[level]; - cur_loc = frm_pyr->level_loc[level]; - - av1_resize_plane(frm_pyr->level_buffer + frm_pyr->level_loc[level - 1], - frm_pyr->heights[level - 1], frm_pyr->widths[level - 1], - frm_pyr->strides[level - 1], - frm_pyr->level_buffer + cur_loc, cur_height, cur_width, - cur_stride); - - if (compute_grad) { - assert(frm_pyr->has_gradient && frm_pyr->level_dx_buffer != NULL && - frm_pyr->level_dy_buffer != NULL); - // Computation x gradient - sobel_xy_image_gradient(frm_pyr->level_buffer + cur_loc, cur_stride, - frm_pyr->level_dx_buffer + cur_loc, cur_stride, - cur_height, cur_width, 1); - - // Computation y gradient - sobel_xy_image_gradient(frm_pyr->level_buffer + cur_loc, cur_stride, - frm_pyr->level_dy_buffer + cur_loc, cur_stride, - cur_height, cur_width, 0); - } - } -} - -static INLINE void compute_flow_at_point(unsigned char *frm, unsigned char *ref, - double *dx, double *dy, int x, int y, - int width, int height, int stride, - double *u, double *v) { - double M[4] = { 0 }; - double b[2] = { 0 }; - double tmp_output_vec[2] = { 0 }; - double error = 0; - int16_t dt[PATCH_SIZE * PATCH_SIZE]; - double o_u = *u; - double o_v = *v; - - for (int itr = 0; itr < DISFLOW_MAX_ITR; itr++) { - error = compute_warp_and_error(ref, frm, width, height, stride, x, y, *u, - *v, dt); - if (error <= DISFLOW_ERROR_TR) break; - compute_flow_system(dx, stride, dy, stride, dt, PATCH_SIZE, M, b); - solve_2x2_system(M, b, tmp_output_vec); - *u += tmp_output_vec[0]; - *v += tmp_output_vec[1]; - } - if (fabs(*u - o_u) > PATCH_SIZE || fabs(*v - o_u) > PATCH_SIZE) { - *u = o_u; - *v = o_v; - } -} - -// make sure flow_u and flow_v start at 0 -static void compute_flow_field(ImagePyramid *frm_pyr, ImagePyramid *ref_pyr, - double *flow_u, double *flow_v) { - int cur_width, cur_height, cur_stride, cur_loc, patch_loc, patch_center; - double *u_upscale = - aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_u)); - double *v_upscale = - aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_v)); - - assert(frm_pyr->n_levels == ref_pyr->n_levels); - - // Compute flow field from coarsest to finest level of the pyramid - for (int level = frm_pyr->n_levels - 1; level >= 0; --level) { - cur_width = frm_pyr->widths[level]; - cur_height = frm_pyr->heights[level]; - cur_stride = frm_pyr->strides[level]; - cur_loc = frm_pyr->level_loc[level]; - - for (int i = PATCH_SIZE; i < cur_height - PATCH_SIZE; i += PATCH_STEP) { - for (int j = PATCH_SIZE; j < cur_width - PATCH_SIZE; j += PATCH_STEP) { - patch_loc = i * cur_stride + j; - patch_center = patch_loc + PATCH_CENTER * cur_stride + PATCH_CENTER; - compute_flow_at_point(frm_pyr->level_buffer + cur_loc, - ref_pyr->level_buffer + cur_loc, - frm_pyr->level_dx_buffer + cur_loc + patch_loc, - frm_pyr->level_dy_buffer + cur_loc + patch_loc, j, - i, cur_width, cur_height, cur_stride, - flow_u + patch_center, flow_v + patch_center); - } - } - // TODO(sarahparker) Replace this with upscale function in resize.c - if (level > 0) { - int h_upscale = frm_pyr->heights[level - 1]; - int w_upscale = frm_pyr->widths[level - 1]; - int s_upscale = frm_pyr->strides[level - 1]; - for (int i = 0; i < h_upscale; ++i) { - for (int j = 0; j < w_upscale; ++j) { - u_upscale[j + i * s_upscale] = - flow_u[(int)(j >> 1) + (int)(i >> 1) * cur_stride]; - v_upscale[j + i * s_upscale] = - flow_v[(int)(j >> 1) + (int)(i >> 1) * cur_stride]; - } - } - memcpy(flow_u, u_upscale, - frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_u)); - memcpy(flow_v, v_upscale, - frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_v)); - } - } - aom_free(u_upscale); - aom_free(v_upscale); -} - -static int compute_global_motion_disflow_based( - TransformationType type, unsigned char *frm_buffer, int frm_width, - int frm_height, int frm_stride, int *frm_corners, int num_frm_corners, - YV12_BUFFER_CONFIG *ref, int bit_depth, int *num_inliers_by_motion, - MotionModel *params_by_motion, int num_motions) { - unsigned char *ref_buffer = ref->y_buffer; - const int ref_width = ref->y_width; - const int ref_height = ref->y_height; - const int pad_size = AOMMAX(PATCH_SIZE, MIN_PAD); - int num_correspondences; - double *correspondences; - RansacFuncDouble ransac = av1_get_ransac_double_prec_type(type); - assert(frm_width == ref_width); - assert(frm_height == ref_height); - - // Ensure the number of pyramid levels will work with the frame resolution - const int msb = - frm_width < frm_height ? get_msb(frm_width) : get_msb(frm_height); - const int n_levels = AOMMIN(msb, N_LEVELS); - - if (ref->flags & YV12_FLAG_HIGHBITDEPTH) { - ref_buffer = av1_downconvert_frame(ref, bit_depth); - } - - // TODO(sarahparker) We will want to do the source pyramid computation - // outside of this function so it doesn't get recomputed for every - // reference. We also don't need to compute every pyramid level for the - // reference in advance, since lower levels can be overwritten once their - // flow field is computed and upscaled. I'll add these optimizations - // once the full implementation is working. - // Allocate frm image pyramids - int compute_gradient = 1; - ImagePyramid *frm_pyr = - alloc_pyramid(frm_width, frm_height, pad_size, compute_gradient); - compute_flow_pyramids(frm_buffer, frm_width, frm_height, frm_stride, n_levels, - pad_size, compute_gradient, frm_pyr); - // Allocate ref image pyramids - compute_gradient = 0; - ImagePyramid *ref_pyr = - alloc_pyramid(ref_width, ref_height, pad_size, compute_gradient); - compute_flow_pyramids(ref_buffer, ref_width, ref_height, ref->y_stride, - n_levels, pad_size, compute_gradient, ref_pyr); - - double *flow_u = - aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_u)); - double *flow_v = - aom_malloc(frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_v)); - - memset(flow_u, 0, - frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_u)); - memset(flow_v, 0, - frm_pyr->strides[0] * frm_pyr->heights[0] * sizeof(*flow_v)); - - compute_flow_field(frm_pyr, ref_pyr, flow_u, flow_v); - - // find correspondences between the two images using the flow field - correspondences = aom_malloc(num_frm_corners * 4 * sizeof(*correspondences)); - num_correspondences = determine_disflow_correspondence( - frm_corners, num_frm_corners, flow_u, flow_v, frm_width, frm_height, - frm_pyr->strides[0], correspondences); - ransac(correspondences, num_correspondences, num_inliers_by_motion, - params_by_motion, num_motions); - - free_pyramid(frm_pyr); - free_pyramid(ref_pyr); - aom_free(correspondences); - aom_free(flow_u); - aom_free(flow_v); - // Set num_inliers = 0 for motions with too few inliers so they are ignored. - for (int i = 0; i < num_motions; ++i) { - if (num_inliers_by_motion[i] < MIN_INLIER_PROB * num_correspondences) { - num_inliers_by_motion[i] = 0; - } - } - - // Return true if any one of the motions has inliers. - for (int i = 0; i < num_motions; ++i) { - if (num_inliers_by_motion[i] > 0) return 1; - } - return 0; -} - -int av1_compute_global_motion(TransformationType type, - unsigned char *src_buffer, int src_width, - int src_height, int src_stride, int *src_corners, - int num_src_corners, YV12_BUFFER_CONFIG *ref, - int bit_depth, - GlobalMotionEstimationType gm_estimation_type, - int *num_inliers_by_motion, - MotionModel *params_by_motion, int num_motions) { - switch (gm_estimation_type) { - case GLOBAL_MOTION_FEATURE_BASED: - return compute_global_motion_feature_based( - type, src_buffer, src_width, src_height, src_stride, src_corners, - num_src_corners, ref, bit_depth, num_inliers_by_motion, - params_by_motion, num_motions); - case GLOBAL_MOTION_DISFLOW_BASED: - return compute_global_motion_disflow_based( - type, src_buffer, src_width, src_height, src_stride, src_corners, - num_src_corners, ref, bit_depth, num_inliers_by_motion, - params_by_motion, num_motions); - default: assert(0 && "Unknown global motion estimation type"); - } - return 0; -}
diff --git a/av1/encoder/global_motion.h b/av1/encoder/global_motion.h index bfc1fec..d3347ca 100644 --- a/av1/encoder/global_motion.h +++ b/av1/encoder/global_motion.h
@@ -14,34 +14,22 @@ #define AOM_AV1_ENCODER_GLOBAL_MOTION_H_ #include "aom/aom_integer.h" +#include "aom_dsp/flow_estimation/flow_estimation.h" #include "aom_scale/yv12config.h" #include "aom_util/aom_thread.h" #include "av1/common/mv.h" #include "av1/common/warped_motion.h" +#include "av1/encoder/cost.h" #ifdef __cplusplus extern "C" { #endif -#define MAX_CORNERS 4096 #define RANSAC_NUM_MOTIONS 1 #define GM_REFINEMENT_COUNT 5 #define MAX_DIRECTIONS 2 -typedef enum { - GLOBAL_MOTION_FEATURE_BASED, - GLOBAL_MOTION_DISFLOW_BASED, -} GlobalMotionEstimationType; - -unsigned char *av1_downconvert_frame(YV12_BUFFER_CONFIG *frm, int bit_depth); - -typedef struct { - double params[MAX_PARAMDIM - 1]; - int *inliers; - int num_inliers; -} MotionModel; - // The structure holds a valid reference frame type and its temporal distance // from the source frame. typedef struct { @@ -97,6 +85,30 @@ int8_t allocated_workers; } AV1GlobalMotionSync; +// Cost (in bits) to signal a given global motion type +// This is only the cost of signaling a given type, and does not +// take into account the cost of the associated parameters. +// +// Unusable types are assigned a cost of 1000 bits here so that, +// even if they somehow make it through the encoder, they are extremely +// unlikely to be chosen +#define GM_TYPE_UNUSABLE_COST (1000 << AV1_PROB_COST_SHIFT) +static const int gm_type_cost[TRANS_TYPES] = { + 1 << AV1_PROB_COST_SHIFT, // IDENTITY + 3 << AV1_PROB_COST_SHIFT, // TRANSLATION + GM_TYPE_UNUSABLE_COST, // ROTATION + GM_TYPE_UNUSABLE_COST, // ZOOM + GM_TYPE_UNUSABLE_COST, // VERTSHEAR + GM_TYPE_UNUSABLE_COST, // HORZSHEAR + GM_TYPE_UNUSABLE_COST, // UZOOM + 2 << AV1_PROB_COST_SHIFT, // ROTZOOM + GM_TYPE_UNUSABLE_COST, // ROTUZOOM + 3 << AV1_PROB_COST_SHIFT, // AFFINE + GM_TYPE_UNUSABLE_COST, // VERTRAPEZOID + GM_TYPE_UNUSABLE_COST, // HORTRAPEZOID + GM_TYPE_UNUSABLE_COST, // HOMOGRAPHY +}; + void av1_convert_model_to_params(const double *params, WarpedMotionParams *model); @@ -130,29 +142,6 @@ int64_t best_frame_error, uint8_t *segment_map, int segment_map_stride, int64_t erroradv_threshold); -/* - Computes "num_motions" candidate global motion parameters between two frames. - The array "params_by_motion" should be length 8 * "num_motions". The ordering - of each set of parameters is best described by the homography: - - [x' (m2 m3 m0 [x - z . y' = m4 m5 m1 * y - 1] m6 m7 1) 1] - - where m{i} represents the ith value in any given set of parameters. - - "num_inliers" should be length "num_motions", and will be populated with the - number of inlier feature points for each motion. Params for which the - num_inliers entry is 0 should be ignored by the caller. -*/ -int av1_compute_global_motion(TransformationType type, - unsigned char *src_buffer, int src_width, - int src_height, int src_stride, int *src_corners, - int num_src_corners, YV12_BUFFER_CONFIG *ref, - int bit_depth, - GlobalMotionEstimationType gm_estimation_type, - int *num_inliers_by_motion, - MotionModel *params_by_motion, int num_motions); #ifdef __cplusplus } // extern "C" #endif
diff --git a/av1/encoder/global_motion_facade.c b/av1/encoder/global_motion_facade.c index 1dda3c5..fa4a688 100644 --- a/av1/encoder/global_motion_facade.c +++ b/av1/encoder/global_motion_facade.c
@@ -13,16 +13,21 @@ #include "aom_dsp/binary_codes_writer.h" #include "aom_ports/system_state.h" +#include "aom_dsp/flow_estimation/flow_estimation.h" +#include "aom_dsp/flow_estimation/corner_detect.h" + #if CONFIG_FLEX_MVRES #include "av1/common/mv.h" #endif -#include "av1/encoder/corner_detect.h" #include "av1/encoder/encoder.h" #include "av1/encoder/ethread.h" #include "av1/encoder/rdopt.h" -// Highest motion model to search. -#define GLOBAL_TRANS_TYPES_ENC 3 +// Motion models to search +#define NUM_MODELS_TO_SEARCH 1 +static const TransformationType models_to_search[NUM_MODELS_TO_SEARCH] = { + ROTZOOM +}; // Computes the cost for the warp parameters. static int gm_get_params_cost(const WarpedMotionParams *gm, @@ -102,18 +107,19 @@ #else YV12_BUFFER_CONFIG *ref_buf[REF_FRAMES], #endif // CONFIG_NEW_REF_SIGNALING - int frame, int num_src_corners, int *src_corners, unsigned char *src_buffer, - MotionModel *params_by_motion, uint8_t *segment_map, + int frame, MotionModel *params_by_motion, uint8_t *segment_map, const int segment_map_w, const int segment_map_h, const WarpedMotionParams *ref_params) { ThreadData *const td = &cpi->td; MACROBLOCK *const x = &td->mb; AV1_COMMON *const cm = &cpi->common; MACROBLOCKD *const xd = &x->e_mbd; + GlobalMotionInfo *const gm_info = &cpi->gm_info; int i; int src_width = cpi->source->y_width; int src_height = cpi->source->y_height; int src_stride = cpi->source->y_stride; + int bit_depth = cm->seq_params.bit_depth; // clang-format off static const double kIdentityParams[MAX_PARAMDIM - 1] = { 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0 @@ -121,30 +127,27 @@ // clang-format on WarpedMotionParams tmp_wm_params; const double *params_this_motion; - int inliers_by_motion[RANSAC_NUM_MOTIONS]; assert(ref_buf[frame] != NULL); - TransformationType model; aom_clear_system_state(); - // TODO(sarahparker, debargha): Explore do_adaptive_gm_estimation = 1 - const int do_adaptive_gm_estimation = 0; - -#if CONFIG_NEW_REF_SIGNALING - const int ref_frame_dist = get_relative_dist( - &cm->seq_params.order_hint_info, cm->current_frame.order_hint, - cm->cur_frame->ref_order_hints[frame]); -#else - const int ref_frame_dist = get_relative_dist( - &cm->seq_params.order_hint_info, cm->current_frame.order_hint, - cm->cur_frame->ref_order_hints[frame - LAST_FRAME]); -#endif // CONFIG_NEW_REF_SIGNALING +#if CONFIG_GM_USE_DISFLOW + // TODO(rachelbarker): Test a hybrid method, where we use disflow for + // refs with order hint distance <= 2, and corner matching for the rest const GlobalMotionEstimationType gm_estimation_type = - cm->seq_params.order_hint_info.enable_order_hint && - abs(ref_frame_dist) <= 2 && do_adaptive_gm_estimation - ? GLOBAL_MOTION_DISFLOW_BASED - : GLOBAL_MOTION_FEATURE_BASED; - for (model = ROTZOOM; model < GLOBAL_TRANS_TYPES_ENC; ++model) { + GLOBAL_MOTION_DISFLOW_BASED; +#else + const GlobalMotionEstimationType gm_estimation_type = + GLOBAL_MOTION_FEATURE_BASED; +#endif // CONFIG_GM_USE_DISFLOW + + gm_info->flow_data[frame] = aom_compute_flow_data( + cpi->source, ref_buf[frame], bit_depth, gm_estimation_type); + FlowData *flow_data = gm_info->flow_data[frame]; + + for (int model_type_index = 0; model_type_index < NUM_MODELS_TO_SEARCH; + model_type_index++) { + TransformationType model_type = models_to_search[model_type_index]; int64_t best_warp_error = INT64_MAX; // Initially set all params to identity. for (i = 0; i < RANSAC_NUM_MOTIONS; ++i) { @@ -153,14 +156,12 @@ params_by_motion[i].num_inliers = 0; } - av1_compute_global_motion(model, src_buffer, src_width, src_height, - src_stride, src_corners, num_src_corners, - ref_buf[frame], cpi->common.seq_params.bit_depth, - gm_estimation_type, inliers_by_motion, - params_by_motion, RANSAC_NUM_MOTIONS); + aom_fit_global_motion_model(flow_data, model_type, cpi->source, bit_depth, + params_by_motion, RANSAC_NUM_MOTIONS); + int64_t ref_frame_error = 0; for (i = 0; i < RANSAC_NUM_MOTIONS; ++i) { - if (inliers_by_motion[i] == 0) continue; + if (params_by_motion[i].num_inliers == 0) continue; params_this_motion = params_by_motion[i].params; av1_convert_model_to_params(params_this_motion, &tmp_wm_params); @@ -250,18 +251,17 @@ #else YV12_BUFFER_CONFIG *ref_buf[REF_FRAMES], #endif // CONFIG_NEW_REF_SIGNALING - int frame, int num_src_corners, int *src_corners, unsigned char *src_buffer, - MotionModel *params_by_motion, uint8_t *segment_map, int segment_map_w, - int segment_map_h) { + int frame, MotionModel *params_by_motion, uint8_t *segment_map, + int segment_map_w, int segment_map_h) { AV1_COMMON *const cm = &cpi->common; GlobalMotionInfo *const gm_info = &cpi->gm_info; const WarpedMotionParams *ref_params = cm->prev_frame ? &cm->prev_frame->global_motion[frame] : &default_warp_params; - compute_global_motion_for_ref_frame( - cpi, ref_buf, frame, num_src_corners, src_corners, src_buffer, - params_by_motion, segment_map, segment_map_w, segment_map_h, ref_params); + compute_global_motion_for_ref_frame(cpi, ref_buf, frame, params_by_motion, + segment_map, segment_map_w, segment_map_h, + ref_params); gm_info->params_cost[frame] = gm_get_params_cost(&cm->global_motion[frame], ref_params, @@ -270,8 +270,7 @@ #else cm->features.fr_mv_precision) + #endif - gm_info->type_cost[cm->global_motion[frame].wmtype] - - gm_info->type_cost[IDENTITY]; + gm_type_cost[cm->global_motion[frame].wmtype] - gm_type_cost[IDENTITY]; } // Loops over valid reference frames and computes global motion estimation. @@ -284,19 +283,16 @@ YV12_BUFFER_CONFIG *ref_buf[REF_FRAMES], FrameDistPair reference_frame[REF_FRAMES - 1], #endif // CONFIG_NEW_REF_SIGNALING - int num_ref_frames, int num_src_corners, int *src_corners, - unsigned char *src_buffer, MotionModel *params_by_motion, - uint8_t *segment_map, const int segment_map_w, const int segment_map_h) { - // Computation of frame corners for the source frame will be done already. - assert(num_src_corners != -1); + int num_ref_frames, MotionModel *params_by_motion, uint8_t *segment_map, + const int segment_map_w, const int segment_map_h) { AV1_COMMON *const cm = &cpi->common; // Compute global motion w.r.t. reference frames starting from the nearest ref // frame in a given direction. for (int frame = 0; frame < num_ref_frames; frame++) { int ref_frame = reference_frame[frame].frame; - av1_compute_gm_for_valid_ref_frames( - cpi, ref_buf, ref_frame, num_src_corners, src_corners, src_buffer, - params_by_motion, segment_map, segment_map_w, segment_map_h); + av1_compute_gm_for_valid_ref_frames(cpi, ref_buf, ref_frame, + params_by_motion, segment_map, + segment_map_w, segment_map_h); // If global motion w.r.t. current ref frame is // INVALID/TRANSLATION/IDENTITY, skip the evaluation of global motion w.r.t // the remaining ref frames in that direction. The below exit is disabled @@ -499,14 +495,6 @@ GlobalMotionInfo *const gm_info = &cpi->gm_info; YV12_BUFFER_CONFIG *source = cpi->source; - gm_info->src_buffer = source->y_buffer; - if (source->flags & YV12_FLAG_HIGHBITDEPTH) { - // The source buffer is 16-bit, so we need to convert to 8 bits for the - // following code. We cache the result until the source frame is released. - gm_info->src_buffer = - av1_downconvert_frame(source, cpi->common.seq_params.bit_depth); - } - gm_info->segment_map_w = (source->y_width + WARP_ERROR_BLOCK) >> WARP_ERROR_BLOCK_LOG; gm_info->segment_map_h = @@ -535,15 +523,6 @@ sizeof(gm_info->reference_frames[0][0]), compare_distance); qsort(gm_info->reference_frames[1], gm_info->num_ref_frames[1], sizeof(gm_info->reference_frames[1][0]), compare_distance); - - gm_info->num_src_corners = -1; - // If atleast one valid reference frame exists in past/future directions, - // compute interest points of source frame using FAST features. - if (gm_info->num_ref_frames[0] > 0 || gm_info->num_ref_frames[1] > 0) { - gm_info->num_src_corners = av1_fast_corner_detect( - gm_info->src_buffer, source->y_width, source->y_height, - source->y_stride, gm_info->src_corners, MAX_CORNERS); - } } // Computes global motion w.r.t. valid reference frames. @@ -561,9 +540,8 @@ if (gm_info->num_ref_frames[dir] > 0) compute_global_motion_for_references( cpi, gm_info->ref_buf, gm_info->reference_frames[dir], - gm_info->num_ref_frames[dir], gm_info->num_src_corners, - gm_info->src_corners, gm_info->src_buffer, params_by_motion, - segment_map, gm_info->segment_map_w, gm_info->segment_map_h); + gm_info->num_ref_frames[dir], params_by_motion, segment_map, + gm_info->segment_map_w, gm_info->segment_map_h); } dealloc_global_motion_data(params_by_motion, segment_map); @@ -591,3 +569,15 @@ memcpy(cm->cur_frame->global_motion, cm->global_motion, sizeof(cm->cur_frame->global_motion)); } + +// After encoding each frame, this function should be called to free any +// flow fields which were allocated +void av1_free_flow_fields(AV1_COMP *cpi) { + GlobalMotionInfo *const gm_info = &cpi->gm_info; + for (int ref = 0; ref < REF_FRAMES; ref++) { + if (gm_info->flow_data[ref] != NULL) { + aom_free_flow_data(gm_info->flow_data[ref]); + gm_info->flow_data[ref] = NULL; + } + } +}
diff --git a/av1/encoder/global_motion_facade.h b/av1/encoder/global_motion_facade.h index 474cebb..65dba3d 100644 --- a/av1/encoder/global_motion_facade.h +++ b/av1/encoder/global_motion_facade.h
@@ -26,10 +26,14 @@ #else YV12_BUFFER_CONFIG *ref_buf[REF_FRAMES], #endif // CONFIG_NEW_REF_SIGNALING - int frame, int num_src_corners, int *src_corners, unsigned char *src_buffer, - MotionModel *params_by_motion, uint8_t *segment_map, int segment_map_w, - int segment_map_h); + int frame, MotionModel *params_by_motion, uint8_t *segment_map, + int segment_map_w, int segment_map_h); void av1_compute_global_motion_facade(struct AV1_COMP *cpi); + +// After encoding each frame, this function should be called to free any +// flow fields which were allocated +void av1_free_flow_fields(AV1_COMP *cpi); + #ifdef __cplusplus } // extern "C" #endif
diff --git a/av1/encoder/optical_flow.c b/av1/encoder/optical_flow.c index 5899fc2..54e0850 100644 --- a/av1/encoder/optical_flow.c +++ b/av1/encoder/optical_flow.c
@@ -13,9 +13,9 @@ #include <limits.h> #include "config/aom_config.h" +#include "aom_dsp/linalg.h" #include "av1/common/av1_common_int.h" #include "av1/encoder/encoder.h" -#include "av1/encoder/mathutils.h" #include "av1/encoder/optical_flow.h" #include "av1/encoder/reconinter_enc.h" #include "aom_mem/aom_mem.h"
diff --git a/av1/encoder/pickrst.c b/av1/encoder/pickrst.c index ae1bfd4..e5e1940 100644 --- a/av1/encoder/pickrst.c +++ b/av1/encoder/pickrst.c
@@ -20,6 +20,7 @@ #include "aom_dsp/aom_dsp_common.h" #include "aom_dsp/binary_codes_writer.h" +#include "aom_dsp/linalg.h" #include "aom_dsp/psnr.h" #include "aom_mem/aom_mem.h" #include "aom_ports/mem.h" @@ -30,7 +31,6 @@ #include "av1/encoder/av1_quantize.h" #include "av1/encoder/encoder.h" -#include "av1/encoder/mathutils.h" #include "av1/encoder/picklpf.h" #include "av1/encoder/pickrst.h" @@ -142,7 +142,7 @@ // tile in the frame. SgrprojInfo sgrproj; WienerInfo wiener; - AV1PixelRect tile_rect; + PixelRect tile_rect; } RestSearchCtxt; static AOM_INLINE void rsc_on_tile(void *priv) { @@ -185,7 +185,7 @@ static int64_t try_restoration_unit(const RestSearchCtxt *rsc, const RestorationTileLimits *limits, - const AV1PixelRect *tile_rect, + const PixelRect *tile_rect, const RestorationUnitInfo *rui) { const AV1_COMMON *const cm = rsc->cm; const int plane = rsc->plane; @@ -849,9 +849,8 @@ } static AOM_INLINE void search_sgrproj(const RestorationTileLimits *limits, - const AV1PixelRect *tile, - int rest_unit_idx, void *priv, - int32_t *tmpbuf, + const PixelRect *tile, int rest_unit_idx, + void *priv, int32_t *tmpbuf, RestorationLineBuffers *rlbs) { (void)rlbs; RestSearchCtxt *rsc = (RestSearchCtxt *)priv; @@ -1328,7 +1327,7 @@ #define USE_WIENER_REFINEMENT_SEARCH 1 static int64_t finer_tile_search_wiener(const RestSearchCtxt *rsc, const RestorationTileLimits *limits, - const AV1PixelRect *tile, + const PixelRect *tile, RestorationUnitInfo *rui, int wiener_win) { const int plane_off = (WIENER_WIN - wiener_win) >> 1; @@ -1434,7 +1433,7 @@ } static AOM_INLINE void search_wiener(const RestorationTileLimits *limits, - const AV1PixelRect *tile_rect, + const PixelRect *tile_rect, int rest_unit_idx, void *priv, int32_t *tmpbuf, RestorationLineBuffers *rlbs) { @@ -1569,7 +1568,7 @@ } static AOM_INLINE void search_norestore(const RestorationTileLimits *limits, - const AV1PixelRect *tile_rect, + const PixelRect *tile_rect, int rest_unit_idx, void *priv, int32_t *tmpbuf, RestorationLineBuffers *rlbs) { @@ -1588,7 +1587,7 @@ } static AOM_INLINE void search_switchable(const RestorationTileLimits *limits, - const AV1PixelRect *tile_rect, + const PixelRect *tile_rect, int rest_unit_idx, void *priv, int32_t *tmpbuf, RestorationLineBuffers *rlbs) {
diff --git a/av1/encoder/ransac.c b/av1/encoder/ransac.c deleted file mode 100644 index 53beb15..0000000 --- a/av1/encoder/ransac.c +++ /dev/null
@@ -1,821 +0,0 @@ -/* - * Copyright (c) 2021, Alliance for Open Media. All rights reserved - * - * This source code is subject to the terms of the BSD 3-Clause Clear License - * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear - * License was not distributed with this source code in the LICENSE file, you - * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 - * aomedia.org/license/patent-license/. - */ -#include <memory.h> -#include <math.h> -#include <time.h> -#include <stdio.h> -#include <stdlib.h> -#include <assert.h> - -#include "av1/encoder/ransac.h" -#include "av1/encoder/mathutils.h" -#include "av1/encoder/random.h" - -#define MAX_MINPTS 4 -#define MAX_DEGENERATE_ITER 10 -#define MINPTS_MULTIPLIER 5 - -#define INLIER_THRESHOLD 1.25 -#define MIN_TRIALS 20 - -//////////////////////////////////////////////////////////////////////////////// -// ransac -typedef int (*IsDegenerateFunc)(double *p); -typedef void (*NormalizeFunc)(double *p, int np, double *T); -typedef void (*DenormalizeFunc)(double *params, double *T1, double *T2); -typedef int (*FindTransformationFunc)(int points, double *points1, - double *points2, double *params); -typedef void (*ProjectPointsDoubleFunc)(double *mat, double *points, - double *proj, int n, int stride_points, - int stride_proj); - -static void project_points_double_translation(double *mat, double *points, - double *proj, int n, - int stride_points, - int stride_proj) { - int i; - for (i = 0; i < n; ++i) { - const double x = *(points++), y = *(points++); - *(proj++) = x + mat[0]; - *(proj++) = y + mat[1]; - points += stride_points - 2; - proj += stride_proj - 2; - } -} - -static void project_points_double_rotzoom(double *mat, double *points, - double *proj, int n, - int stride_points, int stride_proj) { - int i; - for (i = 0; i < n; ++i) { - const double x = *(points++), y = *(points++); - *(proj++) = mat[2] * x + mat[3] * y + mat[0]; - *(proj++) = -mat[3] * x + mat[2] * y + mat[1]; - points += stride_points - 2; - proj += stride_proj - 2; - } -} - -static void project_points_double_affine(double *mat, double *points, - double *proj, int n, int stride_points, - int stride_proj) { - int i; - for (i = 0; i < n; ++i) { - const double x = *(points++), y = *(points++); - *(proj++) = mat[2] * x + mat[3] * y + mat[0]; - *(proj++) = mat[4] * x + mat[5] * y + mat[1]; - points += stride_points - 2; - proj += stride_proj - 2; - } -} - -static void normalize_homography(double *pts, int n, double *T) { - double *p = pts; - double mean[2] = { 0, 0 }; - double msqe = 0; - double scale; - int i; - - assert(n > 0); - for (i = 0; i < n; ++i, p += 2) { - mean[0] += p[0]; - mean[1] += p[1]; - } - mean[0] /= n; - mean[1] /= n; - for (p = pts, i = 0; i < n; ++i, p += 2) { - p[0] -= mean[0]; - p[1] -= mean[1]; - msqe += sqrt(p[0] * p[0] + p[1] * p[1]); - } - msqe /= n; - scale = (msqe == 0 ? 1.0 : sqrt(2) / msqe); - T[0] = scale; - T[1] = 0; - T[2] = -scale * mean[0]; - T[3] = 0; - T[4] = scale; - T[5] = -scale * mean[1]; - T[6] = 0; - T[7] = 0; - T[8] = 1; - for (p = pts, i = 0; i < n; ++i, p += 2) { - p[0] *= scale; - p[1] *= scale; - } -} - -static void invnormalize_mat(double *T, double *iT) { - double is = 1.0 / T[0]; - double m0 = -T[2] * is; - double m1 = -T[5] * is; - iT[0] = is; - iT[1] = 0; - iT[2] = m0; - iT[3] = 0; - iT[4] = is; - iT[5] = m1; - iT[6] = 0; - iT[7] = 0; - iT[8] = 1; -} - -static void denormalize_homography(double *params, double *T1, double *T2) { - double iT2[9]; - double params2[9]; - invnormalize_mat(T2, iT2); - multiply_mat(params, T1, params2, 3, 3, 3); - multiply_mat(iT2, params2, params, 3, 3, 3); -} - -static void denormalize_affine_reorder(double *params, double *T1, double *T2) { - double params_denorm[MAX_PARAMDIM]; - params_denorm[0] = params[0]; - params_denorm[1] = params[1]; - params_denorm[2] = params[4]; - params_denorm[3] = params[2]; - params_denorm[4] = params[3]; - params_denorm[5] = params[5]; - params_denorm[6] = params_denorm[7] = 0; - params_denorm[8] = 1; - denormalize_homography(params_denorm, T1, T2); - params[0] = params_denorm[2]; - params[1] = params_denorm[5]; - params[2] = params_denorm[0]; - params[3] = params_denorm[1]; - params[4] = params_denorm[3]; - params[5] = params_denorm[4]; - params[6] = params[7] = 0; -} - -static void denormalize_rotzoom_reorder(double *params, double *T1, - double *T2) { - double params_denorm[MAX_PARAMDIM]; - params_denorm[0] = params[0]; - params_denorm[1] = params[1]; - params_denorm[2] = params[2]; - params_denorm[3] = -params[1]; - params_denorm[4] = params[0]; - params_denorm[5] = params[3]; - params_denorm[6] = params_denorm[7] = 0; - params_denorm[8] = 1; - denormalize_homography(params_denorm, T1, T2); - params[0] = params_denorm[2]; - params[1] = params_denorm[5]; - params[2] = params_denorm[0]; - params[3] = params_denorm[1]; - params[4] = -params[3]; - params[5] = params[2]; - params[6] = params[7] = 0; -} - -static void denormalize_translation_reorder(double *params, double *T1, - double *T2) { - double params_denorm[MAX_PARAMDIM]; - params_denorm[0] = 1; - params_denorm[1] = 0; - params_denorm[2] = params[0]; - params_denorm[3] = 0; - params_denorm[4] = 1; - params_denorm[5] = params[1]; - params_denorm[6] = params_denorm[7] = 0; - params_denorm[8] = 1; - denormalize_homography(params_denorm, T1, T2); - params[0] = params_denorm[2]; - params[1] = params_denorm[5]; - params[2] = params[5] = 1; - params[3] = params[4] = 0; - params[6] = params[7] = 0; -} - -static int find_translation(int np, double *pts1, double *pts2, double *mat) { - int i; - double sx, sy, dx, dy; - double sumx, sumy; - - double T1[9], T2[9]; - normalize_homography(pts1, np, T1); - normalize_homography(pts2, np, T2); - - sumx = 0; - sumy = 0; - for (i = 0; i < np; ++i) { - dx = *(pts2++); - dy = *(pts2++); - sx = *(pts1++); - sy = *(pts1++); - - sumx += dx - sx; - sumy += dy - sy; - } - mat[0] = sumx / np; - mat[1] = sumy / np; - denormalize_translation_reorder(mat, T1, T2); - return 0; -} - -static int find_rotzoom(int np, double *pts1, double *pts2, double *mat) { - const int np2 = np * 2; - double *a = (double *)aom_malloc(sizeof(*a) * (np2 * 5 + 20)); - double *b = a + np2 * 4; - double *temp = b + np2; - int i; - double sx, sy, dx, dy; - - double T1[9], T2[9]; - normalize_homography(pts1, np, T1); - normalize_homography(pts2, np, T2); - - for (i = 0; i < np; ++i) { - dx = *(pts2++); - dy = *(pts2++); - sx = *(pts1++); - sy = *(pts1++); - - a[i * 2 * 4 + 0] = sx; - a[i * 2 * 4 + 1] = sy; - a[i * 2 * 4 + 2] = 1; - a[i * 2 * 4 + 3] = 0; - a[(i * 2 + 1) * 4 + 0] = sy; - a[(i * 2 + 1) * 4 + 1] = -sx; - a[(i * 2 + 1) * 4 + 2] = 0; - a[(i * 2 + 1) * 4 + 3] = 1; - - b[2 * i] = dx; - b[2 * i + 1] = dy; - } - if (!least_squares(4, a, np2, 4, b, temp, mat)) { - aom_free(a); - return 1; - } - denormalize_rotzoom_reorder(mat, T1, T2); - aom_free(a); - return 0; -} - -static int find_affine(int np, double *pts1, double *pts2, double *mat) { - assert(np > 0); - const int np2 = np * 2; - double *a = (double *)aom_malloc(sizeof(*a) * (np2 * 7 + 42)); - if (a == NULL) return 1; - double *b = a + np2 * 6; - double *temp = b + np2; - int i; - double sx, sy, dx, dy; - - double T1[9], T2[9]; - normalize_homography(pts1, np, T1); - normalize_homography(pts2, np, T2); - - for (i = 0; i < np; ++i) { - dx = *(pts2++); - dy = *(pts2++); - sx = *(pts1++); - sy = *(pts1++); - - a[i * 2 * 6 + 0] = sx; - a[i * 2 * 6 + 1] = sy; - a[i * 2 * 6 + 2] = 0; - a[i * 2 * 6 + 3] = 0; - a[i * 2 * 6 + 4] = 1; - a[i * 2 * 6 + 5] = 0; - a[(i * 2 + 1) * 6 + 0] = 0; - a[(i * 2 + 1) * 6 + 1] = 0; - a[(i * 2 + 1) * 6 + 2] = sx; - a[(i * 2 + 1) * 6 + 3] = sy; - a[(i * 2 + 1) * 6 + 4] = 0; - a[(i * 2 + 1) * 6 + 5] = 1; - - b[2 * i] = dx; - b[2 * i + 1] = dy; - } - if (!least_squares(6, a, np2, 6, b, temp, mat)) { - aom_free(a); - return 1; - } - denormalize_affine_reorder(mat, T1, T2); - aom_free(a); - return 0; -} - -static int get_rand_indices(int npoints, int minpts, int *indices, - unsigned int *seed) { - int i, j; - int ptr = lcg_rand16(seed) % npoints; - if (minpts > npoints) return 0; - indices[0] = ptr; - ptr = (ptr == npoints - 1 ? 0 : ptr + 1); - i = 1; - while (i < minpts) { - int index = lcg_rand16(seed) % npoints; - while (index) { - ptr = (ptr == npoints - 1 ? 0 : ptr + 1); - for (j = 0; j < i; ++j) { - if (indices[j] == ptr) break; - } - if (j == i) index--; - } - indices[i++] = ptr; - } - return 1; -} - -typedef struct { - int num_inliers; - double variance; - int *inlier_indices; -} RANSAC_MOTION; - -// Return -1 if 'a' is a better motion, 1 if 'b' is better, 0 otherwise. -static int compare_motions(const void *arg_a, const void *arg_b) { - const RANSAC_MOTION *motion_a = (RANSAC_MOTION *)arg_a; - const RANSAC_MOTION *motion_b = (RANSAC_MOTION *)arg_b; - - if (motion_a->num_inliers > motion_b->num_inliers) return -1; - if (motion_a->num_inliers < motion_b->num_inliers) return 1; - if (motion_a->variance < motion_b->variance) return -1; - if (motion_a->variance > motion_b->variance) return 1; - return 0; -} - -static int is_better_motion(const RANSAC_MOTION *motion_a, - const RANSAC_MOTION *motion_b) { - return compare_motions(motion_a, motion_b) < 0; -} - -static void copy_points_at_indices(double *dest, const double *src, - const int *indices, int num_points) { - for (int i = 0; i < num_points; ++i) { - const int index = indices[i]; - dest[i * 2] = src[index * 2]; - dest[i * 2 + 1] = src[index * 2 + 1]; - } -} - -static const double kInfiniteVariance = 1e12; - -static void clear_motion(RANSAC_MOTION *motion, int num_points) { - motion->num_inliers = 0; - motion->variance = kInfiniteVariance; - memset(motion->inlier_indices, 0, - sizeof(*motion->inlier_indices) * num_points); -} - -static int ransac(const int *matched_points, int npoints, - int *num_inliers_by_motion, MotionModel *params_by_motion, - int num_desired_motions, int minpts, - IsDegenerateFunc is_degenerate, - FindTransformationFunc find_transformation, - ProjectPointsDoubleFunc projectpoints) { - int trial_count = 0; - int i = 0; - int ret_val = 0; - - unsigned int seed = (unsigned int)npoints; - - int indices[MAX_MINPTS] = { 0 }; - - double *points1, *points2; - double *corners1, *corners2; - double *image1_coord; - - // Store information for the num_desired_motions best transformations found - // and the worst motion among them, as well as the motion currently under - // consideration. - RANSAC_MOTION *motions, *worst_kept_motion = NULL; - RANSAC_MOTION current_motion; - - // Store the parameters and the indices of the inlier points for the motion - // currently under consideration. - double params_this_motion[MAX_PARAMDIM]; - - double *cnp1, *cnp2; - - for (i = 0; i < num_desired_motions; ++i) { - num_inliers_by_motion[i] = 0; - } - if (npoints < minpts * MINPTS_MULTIPLIER || npoints == 0) { - return 1; - } - - points1 = (double *)aom_malloc(sizeof(*points1) * npoints * 2); - points2 = (double *)aom_malloc(sizeof(*points2) * npoints * 2); - corners1 = (double *)aom_malloc(sizeof(*corners1) * npoints * 2); - corners2 = (double *)aom_malloc(sizeof(*corners2) * npoints * 2); - image1_coord = (double *)aom_malloc(sizeof(*image1_coord) * npoints * 2); - - motions = - (RANSAC_MOTION *)aom_malloc(sizeof(RANSAC_MOTION) * num_desired_motions); - for (i = 0; i < num_desired_motions; ++i) { - motions[i].inlier_indices = - (int *)aom_malloc(sizeof(*motions->inlier_indices) * npoints); - clear_motion(motions + i, npoints); - } - current_motion.inlier_indices = - (int *)aom_malloc(sizeof(*current_motion.inlier_indices) * npoints); - clear_motion(¤t_motion, npoints); - - worst_kept_motion = motions; - - if (!(points1 && points2 && corners1 && corners2 && image1_coord && motions && - current_motion.inlier_indices)) { - ret_val = 1; - goto finish_ransac; - } - - cnp1 = corners1; - cnp2 = corners2; - for (i = 0; i < npoints; ++i) { - *(cnp1++) = *(matched_points++); - *(cnp1++) = *(matched_points++); - *(cnp2++) = *(matched_points++); - *(cnp2++) = *(matched_points++); - } - - while (MIN_TRIALS > trial_count) { - double sum_distance = 0.0; - double sum_distance_squared = 0.0; - - clear_motion(¤t_motion, npoints); - - int degenerate = 1; - int num_degenerate_iter = 0; - - while (degenerate) { - num_degenerate_iter++; - if (!get_rand_indices(npoints, minpts, indices, &seed)) { - ret_val = 1; - goto finish_ransac; - } - - copy_points_at_indices(points1, corners1, indices, minpts); - copy_points_at_indices(points2, corners2, indices, minpts); - - degenerate = is_degenerate(points1); - if (num_degenerate_iter > MAX_DEGENERATE_ITER) { - ret_val = 1; - goto finish_ransac; - } - } - - if (find_transformation(minpts, points1, points2, params_this_motion)) { - trial_count++; - continue; - } - - projectpoints(params_this_motion, corners1, image1_coord, npoints, 2, 2); - - for (i = 0; i < npoints; ++i) { - double dx = image1_coord[i * 2] - corners2[i * 2]; - double dy = image1_coord[i * 2 + 1] - corners2[i * 2 + 1]; - double distance = sqrt(dx * dx + dy * dy); - - if (distance < INLIER_THRESHOLD) { - current_motion.inlier_indices[current_motion.num_inliers++] = i; - sum_distance += distance; - sum_distance_squared += distance * distance; - } - } - - if (current_motion.num_inliers >= worst_kept_motion->num_inliers && - current_motion.num_inliers > 1) { - double mean_distance; - mean_distance = sum_distance / ((double)current_motion.num_inliers); - current_motion.variance = - sum_distance_squared / ((double)current_motion.num_inliers - 1.0) - - mean_distance * mean_distance * ((double)current_motion.num_inliers) / - ((double)current_motion.num_inliers - 1.0); - if (is_better_motion(¤t_motion, worst_kept_motion)) { - // This motion is better than the worst currently kept motion. Remember - // the inlier points and variance. The parameters for each kept motion - // will be recomputed later using only the inliers. - worst_kept_motion->num_inliers = current_motion.num_inliers; - worst_kept_motion->variance = current_motion.variance; - memcpy(worst_kept_motion->inlier_indices, current_motion.inlier_indices, - sizeof(*current_motion.inlier_indices) * npoints); - assert(npoints > 0); - // Determine the new worst kept motion and its num_inliers and variance. - for (i = 0; i < num_desired_motions; ++i) { - if (is_better_motion(worst_kept_motion, &motions[i])) { - worst_kept_motion = &motions[i]; - } - } - } - } - trial_count++; - } - - // Sort the motions, best first. - qsort(motions, num_desired_motions, sizeof(RANSAC_MOTION), compare_motions); - - // Recompute the motions using only the inliers. - for (i = 0; i < num_desired_motions; ++i) { - if (motions[i].num_inliers >= minpts) { - copy_points_at_indices(points1, corners1, motions[i].inlier_indices, - motions[i].num_inliers); - copy_points_at_indices(points2, corners2, motions[i].inlier_indices, - motions[i].num_inliers); - - find_transformation(motions[i].num_inliers, points1, points2, - params_by_motion[i].params); - - params_by_motion[i].num_inliers = motions[i].num_inliers; - memcpy(params_by_motion[i].inliers, motions[i].inlier_indices, - sizeof(*motions[i].inlier_indices) * npoints); - num_inliers_by_motion[i] = motions[i].num_inliers; - } - } - -finish_ransac: - aom_free(points1); - aom_free(points2); - aom_free(corners1); - aom_free(corners2); - aom_free(image1_coord); - aom_free(current_motion.inlier_indices); - for (i = 0; i < num_desired_motions; ++i) { - aom_free(motions[i].inlier_indices); - } - aom_free(motions); - - return ret_val; -} - -static int ransac_double_prec(const double *matched_points, int npoints, - int *num_inliers_by_motion, - MotionModel *params_by_motion, - int num_desired_motions, int minpts, - IsDegenerateFunc is_degenerate, - FindTransformationFunc find_transformation, - ProjectPointsDoubleFunc projectpoints) { - int trial_count = 0; - int i = 0; - int ret_val = 0; - - unsigned int seed = (unsigned int)npoints; - - int indices[MAX_MINPTS] = { 0 }; - - double *points1, *points2; - double *corners1, *corners2; - double *image1_coord; - - // Store information for the num_desired_motions best transformations found - // and the worst motion among them, as well as the motion currently under - // consideration. - RANSAC_MOTION *motions, *worst_kept_motion = NULL; - RANSAC_MOTION current_motion; - - // Store the parameters and the indices of the inlier points for the motion - // currently under consideration. - double params_this_motion[MAX_PARAMDIM]; - - double *cnp1, *cnp2; - - for (i = 0; i < num_desired_motions; ++i) { - num_inliers_by_motion[i] = 0; - } - if (npoints < minpts * MINPTS_MULTIPLIER || npoints == 0) { - return 1; - } - - points1 = (double *)aom_malloc(sizeof(*points1) * npoints * 2); - points2 = (double *)aom_malloc(sizeof(*points2) * npoints * 2); - corners1 = (double *)aom_malloc(sizeof(*corners1) * npoints * 2); - corners2 = (double *)aom_malloc(sizeof(*corners2) * npoints * 2); - image1_coord = (double *)aom_malloc(sizeof(*image1_coord) * npoints * 2); - - motions = - (RANSAC_MOTION *)aom_malloc(sizeof(RANSAC_MOTION) * num_desired_motions); - for (i = 0; i < num_desired_motions; ++i) { - motions[i].inlier_indices = - (int *)aom_malloc(sizeof(*motions->inlier_indices) * npoints); - clear_motion(motions + i, npoints); - } - current_motion.inlier_indices = - (int *)aom_malloc(sizeof(*current_motion.inlier_indices) * npoints); - clear_motion(¤t_motion, npoints); - - worst_kept_motion = motions; - - if (!(points1 && points2 && corners1 && corners2 && image1_coord && motions && - current_motion.inlier_indices)) { - ret_val = 1; - goto finish_ransac; - } - - cnp1 = corners1; - cnp2 = corners2; - for (i = 0; i < npoints; ++i) { - *(cnp1++) = *(matched_points++); - *(cnp1++) = *(matched_points++); - *(cnp2++) = *(matched_points++); - *(cnp2++) = *(matched_points++); - } - - while (MIN_TRIALS > trial_count) { - double sum_distance = 0.0; - double sum_distance_squared = 0.0; - - clear_motion(¤t_motion, npoints); - - int degenerate = 1; - int num_degenerate_iter = 0; - - while (degenerate) { - num_degenerate_iter++; - if (!get_rand_indices(npoints, minpts, indices, &seed)) { - ret_val = 1; - goto finish_ransac; - } - - copy_points_at_indices(points1, corners1, indices, minpts); - copy_points_at_indices(points2, corners2, indices, minpts); - - degenerate = is_degenerate(points1); - if (num_degenerate_iter > MAX_DEGENERATE_ITER) { - ret_val = 1; - goto finish_ransac; - } - } - - if (find_transformation(minpts, points1, points2, params_this_motion)) { - trial_count++; - continue; - } - - projectpoints(params_this_motion, corners1, image1_coord, npoints, 2, 2); - - for (i = 0; i < npoints; ++i) { - double dx = image1_coord[i * 2] - corners2[i * 2]; - double dy = image1_coord[i * 2 + 1] - corners2[i * 2 + 1]; - double distance = sqrt(dx * dx + dy * dy); - - if (distance < INLIER_THRESHOLD) { - current_motion.inlier_indices[current_motion.num_inliers++] = i; - sum_distance += distance; - sum_distance_squared += distance * distance; - } - } - - if (current_motion.num_inliers >= worst_kept_motion->num_inliers && - current_motion.num_inliers > 1) { - double mean_distance; - mean_distance = sum_distance / ((double)current_motion.num_inliers); - current_motion.variance = - sum_distance_squared / ((double)current_motion.num_inliers - 1.0) - - mean_distance * mean_distance * ((double)current_motion.num_inliers) / - ((double)current_motion.num_inliers - 1.0); - if (is_better_motion(¤t_motion, worst_kept_motion)) { - // This motion is better than the worst currently kept motion. Remember - // the inlier points and variance. The parameters for each kept motion - // will be recomputed later using only the inliers. - worst_kept_motion->num_inliers = current_motion.num_inliers; - worst_kept_motion->variance = current_motion.variance; - memcpy(worst_kept_motion->inlier_indices, current_motion.inlier_indices, - sizeof(*current_motion.inlier_indices) * npoints); - assert(npoints > 0); - // Determine the new worst kept motion and its num_inliers and variance. - for (i = 0; i < num_desired_motions; ++i) { - if (is_better_motion(worst_kept_motion, &motions[i])) { - worst_kept_motion = &motions[i]; - } - } - } - } - trial_count++; - } - - // Sort the motions, best first. - qsort(motions, num_desired_motions, sizeof(RANSAC_MOTION), compare_motions); - - // Recompute the motions using only the inliers. - for (i = 0; i < num_desired_motions; ++i) { - if (motions[i].num_inliers >= minpts) { - copy_points_at_indices(points1, corners1, motions[i].inlier_indices, - motions[i].num_inliers); - copy_points_at_indices(points2, corners2, motions[i].inlier_indices, - motions[i].num_inliers); - - find_transformation(motions[i].num_inliers, points1, points2, - params_by_motion[i].params); - memcpy(params_by_motion[i].inliers, motions[i].inlier_indices, - sizeof(*motions[i].inlier_indices) * npoints); - } - num_inliers_by_motion[i] = motions[i].num_inliers; - } - -finish_ransac: - aom_free(points1); - aom_free(points2); - aom_free(corners1); - aom_free(corners2); - aom_free(image1_coord); - aom_free(current_motion.inlier_indices); - for (i = 0; i < num_desired_motions; ++i) { - aom_free(motions[i].inlier_indices); - } - aom_free(motions); - - return ret_val; -} - -static int is_collinear3(double *p1, double *p2, double *p3) { - static const double collinear_eps = 1e-3; - const double v = - (p2[0] - p1[0]) * (p3[1] - p1[1]) - (p2[1] - p1[1]) * (p3[0] - p1[0]); - return fabs(v) < collinear_eps; -} - -static int is_degenerate_translation(double *p) { - return (p[0] - p[2]) * (p[0] - p[2]) + (p[1] - p[3]) * (p[1] - p[3]) <= 2; -} - -static int is_degenerate_affine(double *p) { - return is_collinear3(p, p + 2, p + 4); -} - -static int ransac_translation(int *matched_points, int npoints, - int *num_inliers_by_motion, - MotionModel *params_by_motion, - int num_desired_motions) { - return ransac(matched_points, npoints, num_inliers_by_motion, - params_by_motion, num_desired_motions, 3, - is_degenerate_translation, find_translation, - project_points_double_translation); -} - -static int ransac_rotzoom(int *matched_points, int npoints, - int *num_inliers_by_motion, - MotionModel *params_by_motion, - int num_desired_motions) { - return ransac(matched_points, npoints, num_inliers_by_motion, - params_by_motion, num_desired_motions, 3, is_degenerate_affine, - find_rotzoom, project_points_double_rotzoom); -} - -static int ransac_affine(int *matched_points, int npoints, - int *num_inliers_by_motion, - MotionModel *params_by_motion, - int num_desired_motions) { - return ransac(matched_points, npoints, num_inliers_by_motion, - params_by_motion, num_desired_motions, 3, is_degenerate_affine, - find_affine, project_points_double_affine); -} - -RansacFunc av1_get_ransac_type(TransformationType type) { - switch (type) { - case AFFINE: return ransac_affine; - case ROTZOOM: return ransac_rotzoom; - case TRANSLATION: return ransac_translation; - default: assert(0); return NULL; - } -} - -static int ransac_translation_double_prec(double *matched_points, int npoints, - int *num_inliers_by_motion, - MotionModel *params_by_motion, - int num_desired_motions) { - return ransac_double_prec(matched_points, npoints, num_inliers_by_motion, - params_by_motion, num_desired_motions, 3, - is_degenerate_translation, find_translation, - project_points_double_translation); -} - -static int ransac_rotzoom_double_prec(double *matched_points, int npoints, - int *num_inliers_by_motion, - MotionModel *params_by_motion, - int num_desired_motions) { - return ransac_double_prec(matched_points, npoints, num_inliers_by_motion, - params_by_motion, num_desired_motions, 3, - is_degenerate_affine, find_rotzoom, - project_points_double_rotzoom); -} - -static int ransac_affine_double_prec(double *matched_points, int npoints, - int *num_inliers_by_motion, - MotionModel *params_by_motion, - int num_desired_motions) { - return ransac_double_prec(matched_points, npoints, num_inliers_by_motion, - params_by_motion, num_desired_motions, 3, - is_degenerate_affine, find_affine, - project_points_double_affine); -} - -RansacFuncDouble av1_get_ransac_double_prec_type(TransformationType type) { - switch (type) { - case AFFINE: return ransac_affine_double_prec; - case ROTZOOM: return ransac_rotzoom_double_prec; - case TRANSLATION: return ransac_translation_double_prec; - default: assert(0); return NULL; - } -}
diff --git a/av1/encoder/ransac.h b/av1/encoder/ransac.h deleted file mode 100644 index 18e5468..0000000 --- a/av1/encoder/ransac.h +++ /dev/null
@@ -1,32 +0,0 @@ -/* - * Copyright (c) 2021, Alliance for Open Media. All rights reserved - * - * This source code is subject to the terms of the BSD 3-Clause Clear License - * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear - * License was not distributed with this source code in the LICENSE file, you - * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 - * aomedia.org/license/patent-license/. - */ - -#ifndef AOM_AV1_ENCODER_RANSAC_H_ -#define AOM_AV1_ENCODER_RANSAC_H_ - -#include <stdio.h> -#include <stdlib.h> -#include <math.h> -#include <memory.h> - -#include "av1/common/warped_motion.h" -#include "av1/encoder/global_motion.h" - -typedef int (*RansacFunc)(int *matched_points, int npoints, - int *num_inliers_by_motion, - MotionModel *params_by_motion, int num_motions); -typedef int (*RansacFuncDouble)(double *matched_points, int npoints, - int *num_inliers_by_motion, - MotionModel *params_by_motion, int num_motions); -RansacFunc av1_get_ransac_type(TransformationType type); -RansacFuncDouble av1_get_ransac_double_prec_type(TransformationType type); -#endif // AOM_AV1_ENCODER_RANSAC_H_
diff --git a/av1/encoder/rd.c b/av1/encoder/rd.c index 2868724..8da584a 100644 --- a/av1/encoder/rd.c +++ b/av1/encoder/rd.c
@@ -970,16 +970,6 @@ dvcost, &cm->fc->ndvc, MV_SUBPEL_NONE); #endif } - - if (!is_stat_generation_stage(cpi)) { - for (int i = 0; i < TRANS_TYPES; ++i) - // IDENTITY: 1 bit - // TRANSLATION: 3 bits - // ROTZOOM: 2 bits - // AFFINE: 3 bits - cpi->gm_info.type_cost[i] = (1 + (i > 0 ? (i == ROTZOOM ? 1 : 2) : 0)) - << AV1_PROB_COST_SHIFT; - } } static void model_rd_norm(int xsq_q10, int *r_q10, int *d_q10) {
diff --git a/build/cmake/aom_config_defaults.cmake b/build/cmake/aom_config_defaults.cmake index 6cc5bbd..35019f3 100644 --- a/build/cmake/aom_config_defaults.cmake +++ b/build/cmake/aom_config_defaults.cmake
@@ -203,6 +203,16 @@ "AV2 flexible mv precision experiment flag") set_aom_config_var(CONFIG_DERIVED_MV 0 NUMBER "AV2 derived motion vector experiment flag") +set_aom_config_var( + CONFIG_GM_USE_DISFLOW + 0 + "Encoder only: Use disflow method for flow estimation, rather than corner matching" +) +set_aom_config_var( + CONFIG_GM_USE_SRC_FRAMES + 0 + "Encoder-only: Compute flow using source frames rather than reconstructed frames" +) # # Variables in this section control optional features of the build system.
diff --git a/test/corner_match_test.cc b/test/corner_match_test.cc index 39927e6..5c97723 100644 --- a/test/corner_match_test.cc +++ b/test/corner_match_test.cc
@@ -11,7 +11,7 @@ */ #include <tuple> -#include "config/av1_rtcd.h" +#include "config/aom_dsp_rtcd.h" #include "third_party/googletest/src/googletest/include/gtest/gtest.h" #include "test/acm_random.h" @@ -19,7 +19,7 @@ #include "test/clear_system_state.h" #include "test/register_state_check.h" -#include "av1/encoder/corner_match.h" +#include "aom_dsp/flow_estimation/corner_match.h" namespace test_libaom { @@ -92,13 +92,13 @@ int y2 = MATCH_SZ_BY2 + rnd_.PseudoUniform(h - 2 * MATCH_SZ_BY2); double res_c = - av1_compute_cross_correlation_c(input1, w, x1, y1, input2, w, x2, y2); + aom_compute_cross_correlation_c(input1, w, x1, y1, input2, w, x2, y2); double res_simd = target_func(input1, w, x1, y1, input2, w, x2, y2); if (run_times > 1) { aom_usec_timer_start(&ref_timer); for (j = 0; j < run_times; j++) { - av1_compute_cross_correlation_c(input1, w, x1, y1, input2, w, x2, y2); + aom_compute_cross_correlation_c(input1, w, x1, y1, input2, w, x2, y2); } aom_usec_timer_mark(&ref_timer); const int elapsed_time_c = @@ -131,15 +131,15 @@ #if HAVE_SSE4_1 INSTANTIATE_TEST_SUITE_P( SSE4_1, AV1CornerMatchTest, - ::testing::Values(make_tuple(0, &av1_compute_cross_correlation_sse4_1), - make_tuple(1, &av1_compute_cross_correlation_sse4_1))); + ::testing::Values(make_tuple(0, &aom_compute_cross_correlation_sse4_1), + make_tuple(1, &aom_compute_cross_correlation_sse4_1))); #endif #if HAVE_AVX2 INSTANTIATE_TEST_SUITE_P( AVX2, AV1CornerMatchTest, - ::testing::Values(make_tuple(0, &av1_compute_cross_correlation_avx2), - make_tuple(1, &av1_compute_cross_correlation_avx2))); + ::testing::Values(make_tuple(0, &aom_compute_cross_correlation_avx2), + make_tuple(1, &aom_compute_cross_correlation_avx2))); #endif } // namespace AV1CornerMatch
diff --git a/test/flow_estimation_test.cc b/test/flow_estimation_test.cc new file mode 100644 index 0000000..6274926 --- /dev/null +++ b/test/flow_estimation_test.cc
@@ -0,0 +1,356 @@ +/* + * Copyright (c) 2021, Alliance for Open Media. All rights reserved + * + * This source code is subject to the terms of the BSD 3-Clause Clear License + * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear + * License was not distributed with this source code in the LICENSE file, you + * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. 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 + * aomedia.org/license/patent-license/. + */ + +#include <tuple> +#include <stdlib.h> + +// Needed on Windows to define M_PI_4 (== pi/4) +// Source: +// https://docs.microsoft.com/en-us/cpp/c-runtime-library/math-constants?view=msvc-170 +#define _USE_MATH_DEFINES +#include <math.h> +#include <stdbool.h> + +#include "third_party/googletest/src/googletest/include/gtest/gtest.h" + +#include "aom_dsp/flow_estimation/flow_estimation.h" +#include "aom_dsp/flow_estimation/ransac.h" + +#include "test/acm_random.h" +#include "test/util.h" + +namespace { + +using libaom_test::ACMRandom; +using std::tuple; + +typedef tuple<TransformationType> TestParams; + +// Fixed test parameters +const int npoints = 100; +const double noise_level = 0.5; +const int test_iters = 100; + +static const char *model_type_names[] = { + "IDENTITY", "TRANSLATION", "ROTATION", "ZOOM", "VERTSHEAR", + "HORZSHEAR", "UZOOM", "ROTZOOM", "ROTUZOOM", "AFFINE", + "VERTRAPEZOID", "HORTRAPEZOID", "HOMOGRAPHY" +}; + +static double random_double(ACMRandom &rnd, double min_, double max_) { + return min_ + (max_ - min_) * (rnd.Rand31() / (double)(1LL << 31)); +} + +static const double default_model[MAX_PARAMDIM] = { 0.0, 0.0, 1.0, 0.0, + 0.0, 1.0, 0.0, 0.0 }; + +static void generate_model(const TransformationType type, ACMRandom &rnd, + double *model) { + memcpy(model, default_model, sizeof(default_model)); + + switch (type) { + case TRANSLATION: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + break; + + case ROTATION: { + double angle = random_double(rnd, -M_PI_4, M_PI_4); + double c = cos(angle); + double s = sin(angle); + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = c; + model[3] = s; + model[4] = -s; + model[5] = c; + } break; + + case ZOOM: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = 1.0 + random_double(rnd, -0.25, 0.25); + model[3] = 0; + model[4] = 0; + model[5] = model[2]; + break; + + case VERTSHEAR: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = 1.0; + model[3] = 0; + model[4] = random_double(rnd, -0.25, 0.25); + model[5] = 1.0; + break; + + case HORZSHEAR: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = 1.0; + model[3] = random_double(rnd, -0.25, 0.25); + model[4] = 0; + model[5] = 1.0; + break; + + case UZOOM: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = 1.0 + random_double(rnd, -0.25, 0.25); + model[3] = 0; + model[4] = 0; + model[5] = 1.0 + random_double(rnd, -0.25, 0.25); + break; + + case ROTZOOM: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = 1.0 + random_double(rnd, -0.25, 0.25); + model[3] = random_double(rnd, -0.25, 0.25); + model[4] = -model[3]; + model[5] = model[2]; + break; + + case ROTUZOOM: { + // ROTUZOOM models consist of a zoom followed by a rotation, + // which can be expressed as: + // + // ( c s) * (a 0) = ( a*c b*s) + // (-s c) (0 b) (-a*s b*c) + double zoom_x = 1.0 + random_double(rnd, -0.25, 0.25); + double zoom_y = 1.0 + random_double(rnd, -0.25, 0.25); + double angle = random_double(rnd, -M_PI_4, M_PI_4); + double c = cos(angle); + double s = sin(angle); + + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = zoom_x * c; + model[3] = zoom_y * s; + model[4] = -zoom_x * s; + model[5] = zoom_y * c; + } break; + + case AFFINE: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = 1.0 + random_double(rnd, -0.25, 0.25); + model[3] = random_double(rnd, -0.25, 0.25); + model[4] = random_double(rnd, -0.25, 0.25); + model[5] = 1.0 + random_double(rnd, -0.25, 0.25); + break; + + case VERTRAPEZOID: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = 1.0 + random_double(rnd, -0.25, 0.25); + model[3] = 0.0; + model[4] = random_double(rnd, -0.25, 0.25); + model[5] = 1.0 + random_double(rnd, -0.25, 0.25); + model[6] = random_double(rnd, -0.25, 0.25); + model[7] = 0.0; + break; + + case HORTRAPEZOID: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = 1.0 + random_double(rnd, -0.25, 0.25); + model[3] = random_double(rnd, -0.25, 0.25); + model[4] = 0.0; + model[5] = 1.0 + random_double(rnd, -0.25, 0.25); + model[6] = 0.0; + model[7] = random_double(rnd, -0.25, 0.25); + break; + + case HOMOGRAPHY: + model[0] = random_double(rnd, -128.0, 128.0); + model[1] = random_double(rnd, -128.0, 128.0); + model[2] = 1.0 + random_double(rnd, -0.25, 0.25); + model[3] = random_double(rnd, -0.25, 0.25); + model[4] = random_double(rnd, -0.25, 0.25); + model[5] = 1.0 + random_double(rnd, -0.25, 0.25); + model[6] = random_double(rnd, -0.25, 0.25); + model[7] = random_double(rnd, -0.25, 0.25); + break; + + default: assert(0); break; + } +} + +static void apply_model_plus_noise(const int npoints, const double *src_points, + const double *model, ACMRandom &rnd, + const double noise_level, + double *dst_points) { + for (int i = 0; i < npoints; i++) { + double src_x = src_points[2 * i + 0]; + double src_y = src_points[2 * i + 1]; + + double dst_sx = model[0] + src_x * model[2] + src_y * model[3]; + double dst_sy = model[1] + src_x * model[4] + src_y * model[5]; + double dst_s = 1.0 + src_x * model[6] + src_y * model[7]; + + double noise_x = random_double(rnd, -noise_level, noise_level); + double noise_y = random_double(rnd, -noise_level, noise_level); + + double dst_x = dst_sx / dst_s + noise_x; + double dst_y = dst_sy / dst_s + noise_y; + + dst_points[2 * i + 0] = dst_x; + dst_points[2 * i + 1] = dst_y; + } +} + +static double get_rms_err(const int npoints, const double *src_points, + const double *dst_points, const double *model) { + double sse = 0.0; + for (int i = 0; i < npoints; i++) { + double src_x = src_points[2 * i + 0]; + double src_y = src_points[2 * i + 1]; + double dst_x = dst_points[2 * i + 0]; + double dst_y = dst_points[2 * i + 1]; + + double proj_sx = model[0] + src_x * model[2] + src_y * model[3]; + double proj_sy = model[1] + src_x * model[4] + src_y * model[5]; + double proj_s = 1.0 + src_x * model[6] + src_y * model[7]; + + double proj_x = proj_sx / proj_s; + double proj_y = proj_sy / proj_s; + + sse += (proj_x - dst_x) * (proj_x - dst_x); + sse += (proj_y - dst_y) * (proj_y - dst_y); + } + return sqrt(sse / npoints); +} + +static void print_model(double *model) { + printf("{%f, %f, %f, %f, %f, %f, %f, %f}", model[0], model[1], model[2], + model[3], model[4], model[5], model[6], model[7]); +} + +class AomFlowEstimationTest : public ::testing::TestWithParam<TestParams> { + public: + virtual ~AomFlowEstimationTest() {} + virtual void SetUp() {} + virtual void TearDown() {} + + protected: + void RunTest() const { + // Outline: + // * Generate a set of input points + // * Generate a "ground truth" model of the relevant type + // * Apply ground truth model + noise + // * Fit model using aom_fit_motion_model() + // * Compare RMS error of ground truth model and fitted model + + ACMRandom rnd(ACMRandom::DeterministicSeed()); + + double *src_points = (double *)malloc(npoints * 2 * sizeof(*src_points)); + double *dst_points = (double *)malloc(npoints * 2 * sizeof(*dst_points)); + double *src_points2 = (double *)malloc(npoints * 2 * sizeof(*src_points2)); + double *dst_points2 = (double *)malloc(npoints * 2 * sizeof(*dst_points2)); + double ground_truth_model[MAX_PARAMDIM]; + double fitted_model[MAX_PARAMDIM]; + + TransformationType type = GET_PARAM(0); + + for (int iter = 0; iter < test_iters; iter++) { + // Simulate a dataset which could come from an 8K x 4K video frame + for (int i = 0; i < npoints; i++) { + double src_x = random_double(rnd, 0, 8192); + double src_y = random_double(rnd, 0, 4096); + + src_points[2 * i + 0] = src_x; + src_points[2 * i + 1] = src_y; + } + + generate_model(type, rnd, ground_truth_model); + + apply_model_plus_noise(npoints, src_points, ground_truth_model, rnd, + noise_level, dst_points); + + // Copy point arrays, as they will be modified by the fitting code + memcpy(src_points2, src_points, npoints * 2 * sizeof(*src_points)); + memcpy(dst_points2, dst_points, npoints * 2 * sizeof(*dst_points)); + bool result = aom_fit_motion_model(type, npoints, src_points2, + dst_points2, fitted_model); + ASSERT_EQ(result, true) + << "Model fitting failed for type = " << model_type_names[type] + << ", iter = " << iter; + + // Calculate projection errors + double ground_truth_rms = + get_rms_err(npoints, src_points, dst_points, ground_truth_model); + double fitted_rms = + get_rms_err(npoints, src_points, dst_points, fitted_model); + + // Code to aid with debugging +#if 0 + if (type == ... && iter == ...) { + printf("Model type: %s\n", model_type_names[type]); + + printf("Ground truth model: "); + print_model(ground_truth_model); + printf("\n"); + printf("Fitted model: "); + print_model(fitted_model); + printf("\n"); + printf("RMS error: Ground truth = %f, fitted = %f\n", ground_truth_rms, + fitted_rms); + } +#else + // Suppress unused variable warnings + (void)print_model; +#endif + + // In theory, since the models are fitted by a least-squares process, + // we should have fitted_rms <= ground_truth_rms. + // This is because the ground truth model is *a* valid model, and the + // fitted model should minimize the RMS error among *all* valid models. + // + // However, in practice, we want to allow a bit of leeway for numerical + // imprecision. + // + // Note: The trapezoid and homography models seem to have an overall + // error which is very high, and grows greater-than-linearly with the + // noise level, whereas the other models' errors grow remain close to + // optimal. This has not been fully investigated yet, but suggests that + // the condition number of these problems is very high. + double relative_threshold; + if (type <= AFFINE) { + relative_threshold = 1.25; + } else { + relative_threshold = 15.0; + } + ASSERT_LE(fitted_rms, relative_threshold * ground_truth_rms) + << "Fitted model for type = " << model_type_names[type] + << ", iter = " << iter << " is worse than ground truth model"; + } + + free(src_points); + free(dst_points); + free(src_points2); + free(dst_points2); + } +}; + +TEST_P(AomFlowEstimationTest, Test) { RunTest(); } + +INSTANTIATE_TEST_SUITE_P(C, AomFlowEstimationTest, + ::testing::Values(TRANSLATION, ROTATION, ZOOM, + VERTSHEAR, HORZSHEAR, UZOOM, ROTZOOM, + ROTUZOOM, AFFINE + // VERTRAPEZOID, + // HORTRAPEZOID, + // HOMOGRAPHY + )); + +} // namespace
diff --git a/test/test.cmake b/test/test.cmake index 1c7cd8e..ac29fc4 100644 --- a/test/test.cmake +++ b/test/test.cmake
@@ -71,6 +71,7 @@ "${AOM_ROOT}/test/encode_test_driver.cc" "${AOM_ROOT}/test/encode_test_driver.h" "${AOM_ROOT}/test/end_to_end_test.cc" + "${AOM_ROOT}/test/flow_estimation_test.cc" "${AOM_ROOT}/test/gf_pyr_height_test.cc" "${AOM_ROOT}/test/horz_superres_test.cc" "${AOM_ROOT}/test/i420_video_source.h"