| /* |
| * Copyright (c) 2016, Alliance for Open Media. All rights reserved. |
| * |
| * This source code is subject to the terms of the BSD 2 Clause License and |
| * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License |
| * was not distributed with this source code in the LICENSE file, you can |
| * obtain it at www.aomedia.org/license/software. If the Alliance for Open |
| * Media Patent License 1.0 was not distributed with this source code in the |
| * PATENTS file, you can obtain it at www.aomedia.org/license/patent. |
| */ |
| |
| #include <stdlib.h> |
| #include <memory.h> |
| #include <math.h> |
| |
| #include "config/aom_dsp_rtcd.h" |
| |
| #include "aom_dsp/flow_estimation/corner_detect.h" |
| #include "aom_dsp/flow_estimation/corner_match.h" |
| #include "aom_dsp/flow_estimation/disflow.h" |
| #include "aom_dsp/flow_estimation/flow_estimation.h" |
| #include "aom_dsp/flow_estimation/ransac.h" |
| #include "aom_dsp/pyramid.h" |
| #include "aom_scale/yv12config.h" |
| |
| #define THRESHOLD_NCC 0.75 |
| |
| /* Compute mean and standard deviation of pixels in a window of size |
| MATCH_SZ by MATCH_SZ centered at (x, y). |
| Store results into *mean and *one_over_stddev |
| |
| Note: The output of this function is scaled by MATCH_SZ, as in |
| *mean = MATCH_SZ * <true mean> and |
| *one_over_stddev = 1 / (MATCH_SZ * <true stddev>) |
| |
| Combined with the fact that we return 1/stddev rather than the standard |
| deviation itself, this allows us to completely avoid divisions in |
| aom_compute_correlation, which is much hotter than this function is. |
| |
| Returns true if this feature point is usable, false otherwise. |
| */ |
| bool aom_compute_mean_stddev_c(const unsigned char *frame, int stride, int x, |
| int y, double *mean, double *one_over_stddev) { |
| int sum = 0; |
| int sumsq = 0; |
| for (int i = 0; i < MATCH_SZ; ++i) { |
| for (int j = 0; j < MATCH_SZ; ++j) { |
| sum += frame[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)]; |
| sumsq += frame[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)] * |
| frame[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)]; |
| } |
| } |
| *mean = (double)sum / MATCH_SZ; |
| const double variance = sumsq - (*mean) * (*mean); |
| if (variance < MIN_FEATURE_VARIANCE) { |
| *one_over_stddev = 0.0; |
| return false; |
| } |
| *one_over_stddev = 1.0 / sqrt(variance); |
| return true; |
| } |
| |
| /* Compute corr(frame1, frame2) over a window of size MATCH_SZ by MATCH_SZ. |
| To save on computation, the mean and (1 divided by the) standard deviation |
| of the window in each frame are precomputed and passed into this function |
| as arguments. |
| */ |
| double aom_compute_correlation_c(const unsigned char *frame1, int stride1, |
| int x1, int y1, double mean1, |
| double one_over_stddev1, |
| const unsigned char *frame2, int stride2, |
| int x2, int y2, double mean2, |
| double one_over_stddev2) { |
| int v1, v2; |
| int cross = 0; |
| for (int i = 0; i < MATCH_SZ; ++i) { |
| for (int j = 0; j < MATCH_SZ; ++j) { |
| v1 = frame1[(i + y1 - MATCH_SZ_BY2) * stride1 + (j + x1 - MATCH_SZ_BY2)]; |
| v2 = frame2[(i + y2 - MATCH_SZ_BY2) * stride2 + (j + x2 - MATCH_SZ_BY2)]; |
| cross += v1 * v2; |
| } |
| } |
| |
| // Note: In theory, the calculations here "should" be |
| // covariance = cross / N^2 - mean1 * mean2 |
| // correlation = covariance / (stddev1 * stddev2). |
| // |
| // However, because of the scaling in aom_compute_mean_stddev, the |
| // lines below actually calculate |
| // covariance * N^2 = cross - (mean1 * N) * (mean2 * N) |
| // correlation = (covariance * N^2) / ((stddev1 * N) * (stddev2 * N)) |
| // |
| // ie. we have removed the need for a division, and still end up with the |
| // correct unscaled correlation (ie, in the range [-1, +1]) |
| double covariance = cross - mean1 * mean2; |
| double correlation = covariance * (one_over_stddev1 * one_over_stddev2); |
| return correlation; |
| } |
| |
| 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; |
| } |
| |
| typedef struct { |
| int x; |
| int y; |
| double mean; |
| double one_over_stddev; |
| int best_match_idx; |
| double best_match_corr; |
| } PointInfo; |
| |
| static int determine_correspondence(const unsigned char *src, |
| const int *src_corners, int num_src_corners, |
| const unsigned char *ref, |
| const int *ref_corners, int num_ref_corners, |
| int width, int height, int src_stride, |
| int ref_stride, |
| Correspondence *correspondences) { |
| PointInfo *src_point_info = NULL; |
| PointInfo *ref_point_info = NULL; |
| int num_correspondences = 0; |
| |
| src_point_info = |
| (PointInfo *)aom_calloc(num_src_corners, sizeof(*src_point_info)); |
| if (!src_point_info) { |
| goto finished; |
| } |
| |
| ref_point_info = |
| (PointInfo *)aom_calloc(num_ref_corners, sizeof(*ref_point_info)); |
| if (!ref_point_info) { |
| goto finished; |
| } |
| |
| // First pass (linear): |
| // Filter corner lists and compute per-patch means and standard deviations, |
| // for the src and ref frames independently |
| int src_point_count = 0; |
| for (int i = 0; i < num_src_corners; i++) { |
| int src_x = src_corners[2 * i]; |
| int src_y = src_corners[2 * i + 1]; |
| if (!is_eligible_point(src_x, src_y, width, height)) continue; |
| |
| PointInfo *point = &src_point_info[src_point_count]; |
| point->x = src_x; |
| point->y = src_y; |
| point->best_match_corr = THRESHOLD_NCC; |
| if (!aom_compute_mean_stddev(src, src_stride, src_x, src_y, &point->mean, |
| &point->one_over_stddev)) |
| continue; |
| src_point_count++; |
| } |
| if (src_point_count == 0) { |
| goto finished; |
| } |
| |
| int ref_point_count = 0; |
| for (int j = 0; j < num_ref_corners; j++) { |
| int ref_x = ref_corners[2 * j]; |
| int ref_y = ref_corners[2 * j + 1]; |
| if (!is_eligible_point(ref_x, ref_y, width, height)) continue; |
| |
| PointInfo *point = &ref_point_info[ref_point_count]; |
| point->x = ref_x; |
| point->y = ref_y; |
| point->best_match_corr = THRESHOLD_NCC; |
| if (!aom_compute_mean_stddev(ref, ref_stride, ref_x, ref_y, &point->mean, |
| &point->one_over_stddev)) |
| continue; |
| ref_point_count++; |
| } |
| if (ref_point_count == 0) { |
| goto finished; |
| } |
| |
| // Second pass (quadratic): |
| // For each pair of points, compute correlation, and use this to determine |
| // the best match of each corner, in both directions |
| for (int i = 0; i < src_point_count; ++i) { |
| PointInfo *src_point = &src_point_info[i]; |
| for (int j = 0; j < ref_point_count; ++j) { |
| PointInfo *ref_point = &ref_point_info[j]; |
| if (!is_eligible_distance(src_point->x, src_point->y, ref_point->x, |
| ref_point->y, width, height)) |
| continue; |
| |
| double corr = aom_compute_correlation( |
| src, src_stride, src_point->x, src_point->y, src_point->mean, |
| src_point->one_over_stddev, ref, ref_stride, ref_point->x, |
| ref_point->y, ref_point->mean, ref_point->one_over_stddev); |
| |
| if (corr > src_point->best_match_corr) { |
| src_point->best_match_idx = j; |
| src_point->best_match_corr = corr; |
| } |
| if (corr > ref_point->best_match_corr) { |
| ref_point->best_match_idx = i; |
| ref_point->best_match_corr = corr; |
| } |
| } |
| } |
| |
| // Third pass (linear): |
| // Scan through source corners, generating a correspondence for each corner |
| // iff ref_best_match[src_best_match[i]] == i |
| // Then refine the generated correspondences using optical flow |
| for (int i = 0; i < src_point_count; i++) { |
| PointInfo *point = &src_point_info[i]; |
| |
| // Skip corners which were not matched, or which didn't find |
| // a good enough match |
| if (point->best_match_corr < THRESHOLD_NCC) continue; |
| |
| PointInfo *match_point = &ref_point_info[point->best_match_idx]; |
| if (match_point->best_match_idx == i) { |
| // Refine match using optical flow and store |
| const int sx = point->x; |
| const int sy = point->y; |
| const int rx = match_point->x; |
| const int ry = match_point->y; |
| double u = (double)(rx - sx); |
| double v = (double)(ry - sy); |
| |
| const int patch_tl_x = sx - DISFLOW_PATCH_CENTER; |
| const int patch_tl_y = sy - DISFLOW_PATCH_CENTER; |
| |
| aom_compute_flow_at_point(src, ref, patch_tl_x, patch_tl_y, width, height, |
| src_stride, &u, &v); |
| |
| Correspondence *correspondence = &correspondences[num_correspondences]; |
| correspondence->x = (double)sx; |
| correspondence->y = (double)sy; |
| correspondence->rx = (double)sx + u; |
| correspondence->ry = (double)sy + v; |
| num_correspondences++; |
| } |
| } |
| |
| finished: |
| aom_free(src_point_info); |
| aom_free(ref_point_info); |
| return num_correspondences; |
| } |
| |
| bool av1_compute_global_motion_feature_match( |
| TransformationType type, YV12_BUFFER_CONFIG *src, YV12_BUFFER_CONFIG *ref, |
| int bit_depth, int downsample_level, MotionModel *motion_models, |
| int num_motion_models, bool *mem_alloc_failed) { |
| int num_correspondences; |
| Correspondence *correspondences; |
| ImagePyramid *src_pyramid = src->y_pyramid; |
| CornerList *src_corners = src->corners; |
| ImagePyramid *ref_pyramid = ref->y_pyramid; |
| CornerList *ref_corners = ref->corners; |
| |
| // Precompute information we will need about each frame |
| if (aom_compute_pyramid(src, bit_depth, 1, src_pyramid) < 0) { |
| *mem_alloc_failed = true; |
| return false; |
| } |
| if (!av1_compute_corner_list(src, bit_depth, downsample_level, src_corners)) { |
| *mem_alloc_failed = true; |
| return false; |
| } |
| if (aom_compute_pyramid(ref, bit_depth, 1, ref_pyramid) < 0) { |
| *mem_alloc_failed = true; |
| return false; |
| } |
| if (!av1_compute_corner_list(src, bit_depth, downsample_level, ref_corners)) { |
| *mem_alloc_failed = true; |
| return false; |
| } |
| |
| const uint8_t *src_buffer = src_pyramid->layers[0].buffer; |
| const int src_width = src_pyramid->layers[0].width; |
| const int src_height = src_pyramid->layers[0].height; |
| const int src_stride = src_pyramid->layers[0].stride; |
| |
| const uint8_t *ref_buffer = ref_pyramid->layers[0].buffer; |
| assert(ref_pyramid->layers[0].width == src_width); |
| assert(ref_pyramid->layers[0].height == src_height); |
| const int ref_stride = ref_pyramid->layers[0].stride; |
| |
| // find correspondences between the two images |
| correspondences = (Correspondence *)aom_malloc(src_corners->num_corners * |
| sizeof(*correspondences)); |
| if (!correspondences) { |
| *mem_alloc_failed = true; |
| return false; |
| } |
| num_correspondences = determine_correspondence( |
| src_buffer, src_corners->corners, src_corners->num_corners, ref_buffer, |
| ref_corners->corners, ref_corners->num_corners, src_width, src_height, |
| src_stride, ref_stride, correspondences); |
| |
| bool result = ransac(correspondences, num_correspondences, type, |
| motion_models, num_motion_models, mem_alloc_failed); |
| |
| aom_free(correspondences); |
| return result; |
| } |