Clean up RANSAC code
* Merge integer and double precision paths
Only the integer path was being used in practice, which meant that
the double precision path had developed subtle incompatibilities
over time.
Having one unified path prevents this reoccurring.
* Consistently represent success/failure results as booleans, with
true = success and false = failure
* Simplify the calling interface of ransac(), eliminating the need to
request a function pointer and then call that pointer
* Remove redundant num_inliers_by_motion array
* Rename params_by_motion and num_models to motion_models and
num_motion_models, respectively
* Pack correspondence data into an array of N structs, rather than
an array of 4*N integers
* Eliminate intermediate matrix when doing least squares, which could
potentially be quite large (4096 * 6 * 8 bytes)
Change-Id: I4ea57764a2035b92ddf801b12e42619e5fa5f7ef
diff --git a/aom_dsp/flow_estimation/corner_match.c b/aom_dsp/flow_estimation/corner_match.c
index a5ea528..8331a06 100644
--- a/aom_dsp/flow_estimation/corner_match.c
+++ b/aom_dsp/flow_estimation/corner_match.c
@@ -95,19 +95,21 @@
for (i = 0; i < num_correspondences; ++i) {
int x, y, best_x = 0, best_y = 0;
double best_match_ncc = 0.0;
+ // For this algorithm, all points have integer coordinates.
+ // It's a little more efficient to convert them to ints once,
+ // before the inner loops
+ 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(correspondences[i].rx + x,
- correspondences[i].ry + y, width, height))
+ if (!is_eligible_point(rx0 + x, ry0 + y, width, height)) continue;
+ if (!is_eligible_distance(x0, y0, rx0 + x, ry0 + 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);
+ match_ncc = av1_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;
@@ -121,19 +123,18 @@
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(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))
+ 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 = av1_compute_cross_correlation(
- ref, ref_stride, correspondences[i].rx, correspondences[i].ry, frm,
- frm_stride, correspondences[i].x + x, correspondences[i].y + y);
+ 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;
@@ -150,10 +151,10 @@
const unsigned char *ref,
const int *ref_corners, int num_ref_corners,
int width, int height, int src_stride,
- int ref_stride, int *correspondence_pts) {
+ int ref_stride,
+ Correspondence *correspondences) {
// 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;
@@ -198,29 +199,12 @@
return num_correspondences;
}
-static bool get_inliers_from_indices(MotionModel *params,
- int *correspondences) {
- int *inliers_tmp = (int *)aom_calloc(2 * MAX_CORNERS, sizeof(*inliers_tmp));
- if (!inliers_tmp) return false;
-
- 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);
- return true;
-}
-
int av1_compute_global_motion_feature_based(
TransformationType type, YV12_BUFFER_CONFIG *src, YV12_BUFFER_CONFIG *ref,
- int bit_depth, int *num_inliers_by_motion, MotionModel *params_by_motion,
- int num_motions) {
+ int bit_depth, MotionModel *motion_models, int num_motion_models) {
int i;
int num_correspondences;
- int *correspondences;
- RansacFunc ransac = av1_get_ransac_type(type);
+ Correspondence *correspondences;
ImagePyramid *src_pyramid = src->y_pyramid;
CornerList *src_corners = src->corners;
ImagePyramid *ref_pyramid = ref->y_pyramid;
@@ -243,34 +227,30 @@
const int ref_stride = ref_pyramid->layers[0].stride;
// find correspondences between the two images
- correspondences = (int *)aom_malloc(src_corners->num_corners * 4 *
- sizeof(*correspondences));
+ correspondences = (Correspondence *)aom_malloc(src_corners->num_corners *
+ sizeof(*correspondences));
if (!correspondences) return 0;
num_correspondences = aom_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);
- ransac(correspondences, num_correspondences, num_inliers_by_motion,
- params_by_motion, num_motions);
+ ransac(correspondences, num_correspondences, type, motion_models,
+ num_motion_models);
// 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 ||
+ for (i = 0; i < num_motion_models; ++i) {
+ if (motion_models[i].num_inliers < MIN_INLIER_PROB * num_correspondences ||
num_correspondences == 0) {
- num_inliers_by_motion[i] = 0;
- } else if (!get_inliers_from_indices(¶ms_by_motion[i],
- correspondences)) {
- aom_free(correspondences);
- return 0;
+ motion_models[i].num_inliers = 0;
}
}
aom_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;
+ for (i = 0; i < num_motion_models; ++i) {
+ if (motion_models[i].num_inliers > 0) return 1;
}
return 0;
}
diff --git a/aom_dsp/flow_estimation/corner_match.h b/aom_dsp/flow_estimation/corner_match.h
index 2d175d9..77cc501 100644
--- a/aom_dsp/flow_estimation/corner_match.h
+++ b/aom_dsp/flow_estimation/corner_match.h
@@ -28,22 +28,17 @@
#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 aom_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, int *correspondence_pts);
+ int ref_stride,
+ Correspondence *correspondences);
int av1_compute_global_motion_feature_based(
TransformationType type, YV12_BUFFER_CONFIG *src, YV12_BUFFER_CONFIG *ref,
- int bit_depth, int *num_inliers_by_motion, MotionModel *params_by_motion,
- int num_motions);
+ int bit_depth, MotionModel *motion_models, int num_motion_models);
#ifdef __cplusplus
}
diff --git a/aom_dsp/flow_estimation/disflow.c b/aom_dsp/flow_estimation/disflow.c
index 04795a6..107823f 100644
--- a/aom_dsp/flow_estimation/disflow.c
+++ b/aom_dsp/flow_estimation/disflow.c
@@ -63,17 +63,17 @@
int num_frm_corners, double *flow_u,
double *flow_v, int width,
int height, int stride,
- double *correspondences) {
+ 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[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];
+ 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++;
}
}
@@ -541,12 +541,10 @@
int av1_compute_global_motion_disflow_based(
TransformationType type, YV12_BUFFER_CONFIG *frm, YV12_BUFFER_CONFIG *ref,
- int bit_depth, int *num_inliers_by_motion, MotionModel *params_by_motion,
- int num_motions) {
+ int bit_depth, MotionModel *motion_models, int num_motion_models) {
const int pad_size = AOMMAX(PATCH_SIZE, MIN_PAD);
int num_correspondences;
- double *correspondences;
- RansacFuncDouble ransac = av1_get_ransac_double_prec_type(type);
+ Correspondence *correspondences;
ImagePyramid *frm_pyramid = frm->y_pyramid;
CornerList *frm_corners = frm->corners;
ImagePyramid *ref_pyramid = ref->y_pyramid;
@@ -613,25 +611,25 @@
if (!compute_flow_field(frm_pyr, ref_pyr, flow_u, flow_v)) goto Error;
// find correspondences between the two images using the flow field
- correspondences =
- aom_malloc(frm_corners->num_corners * 4 * sizeof(*correspondences));
+ correspondences = (Correspondence *)aom_malloc(frm_corners->num_corners *
+ sizeof(*correspondences));
if (!correspondences) goto Error;
num_correspondences = determine_disflow_correspondence(
frm_corners->corners, frm_corners->num_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);
+ ransac(correspondences, num_correspondences, type, motion_models,
+ num_motion_models);
// 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;
+ for (int i = 0; i < num_motion_models; ++i) {
+ if (motion_models[i].num_inliers < MIN_INLIER_PROB * num_correspondences) {
+ motion_models[i].num_inliers = 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) {
+ for (int i = 0; i < num_motion_models; ++i) {
+ if (motion_models[i].num_inliers > 0) {
ret = 1;
break;
}
diff --git a/aom_dsp/flow_estimation/disflow.h b/aom_dsp/flow_estimation/disflow.h
index 738b21a..771d3f9 100644
--- a/aom_dsp/flow_estimation/disflow.h
+++ b/aom_dsp/flow_estimation/disflow.h
@@ -21,8 +21,7 @@
int av1_compute_global_motion_disflow_based(
TransformationType type, YV12_BUFFER_CONFIG *frm, YV12_BUFFER_CONFIG *ref,
- int bit_depth, int *num_inliers_by_motion, MotionModel *params_by_motion,
- int num_motions);
+ int bit_depth, MotionModel *motion_models, int num_motion_models);
#ifdef __cplusplus
}
diff --git a/aom_dsp/flow_estimation/flow_estimation.c b/aom_dsp/flow_estimation/flow_estimation.c
index f7684a5..ea8ceec 100644
--- a/aom_dsp/flow_estimation/flow_estimation.c
+++ b/aom_dsp/flow_estimation/flow_estimation.c
@@ -21,17 +21,15 @@
int aom_compute_global_motion(TransformationType type, YV12_BUFFER_CONFIG *src,
YV12_BUFFER_CONFIG *ref, int bit_depth,
GlobalMotionEstimationType gm_estimation_type,
- int *num_inliers_by_motion,
- MotionModel *params_by_motion, int num_motions) {
+ MotionModel *motion_models,
+ int num_motion_models) {
switch (gm_estimation_type) {
case GLOBAL_MOTION_FEATURE_BASED:
return av1_compute_global_motion_feature_based(
- type, src, ref, bit_depth, num_inliers_by_motion, params_by_motion,
- num_motions);
+ type, src, ref, bit_depth, motion_models, num_motion_models);
case GLOBAL_MOTION_DISFLOW_BASED:
return av1_compute_global_motion_disflow_based(
- type, src, ref, bit_depth, num_inliers_by_motion, params_by_motion,
- num_motions);
+ type, src, ref, bit_depth, motion_models, num_motion_models);
default: assert(0 && "Unknown global motion estimation type");
}
return 0;
diff --git a/aom_dsp/flow_estimation/flow_estimation.h b/aom_dsp/flow_estimation/flow_estimation.h
index 0b30110..6d84129 100644
--- a/aom_dsp/flow_estimation/flow_estimation.h
+++ b/aom_dsp/flow_estimation/flow_estimation.h
@@ -49,11 +49,16 @@
int num_inliers;
} MotionModel;
+typedef struct {
+ double x, y;
+ double rx, ry;
+} Correspondence;
+
int aom_compute_global_motion(TransformationType type, YV12_BUFFER_CONFIG *src,
YV12_BUFFER_CONFIG *ref, int bit_depth,
GlobalMotionEstimationType gm_estimation_type,
- int *num_inliers_by_motion,
- MotionModel *params_by_motion, int num_motions);
+ MotionModel *motion_models,
+ int num_motion_models);
#ifdef __cplusplus
}
diff --git a/aom_dsp/flow_estimation/ransac.c b/aom_dsp/flow_estimation/ransac.c
index 8ffc30d..d153385 100644
--- a/aom_dsp/flow_estimation/ransac.c
+++ b/aom_dsp/flow_estimation/ransac.c
@@ -13,11 +13,12 @@
#include <math.h>
#include <time.h>
#include <stdio.h>
-#include <stdlib.h>
+#include <stdbool.h>
#include <assert.h>
#include "aom_dsp/flow_estimation/ransac.h"
#include "aom_dsp/mathutils.h"
+#include "aom_mem/aom_mem.h"
// TODO(rachelbarker): Remove dependence on code in av1/encoder/
#include "av1/encoder/random.h"
@@ -29,21 +30,36 @@
#define INLIER_THRESHOLD 1.25
#define MIN_TRIALS 20
+// Flag to enable functions for finding TRANSLATION type models.
+//
+// These modes are not considered currently due to a spec bug (see comments
+// in gm_get_motion_vector() in av1/common/mv.h). Thus we don't need to compile
+// the corresponding search functions, but it is nice to keep the source around
+// but disabled, for completeness.
+#define ALLOW_TRANSLATION_MODELS 0
+
////////////////////////////////////////////////////////////////////////////////
// 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);
+typedef bool (*IsDegenerateFunc)(double *p);
+typedef bool (*FindTransformationFunc)(int points, const double *points1,
+ const double *points2, double *params);
+typedef void (*ProjectPointsFunc)(const double *mat, const 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) {
+// vtable-like structure which stores all of the information needed by RANSAC
+// for a particular model type
+typedef struct {
+ IsDegenerateFunc is_degenerate;
+ FindTransformationFunc find_transformation;
+ ProjectPointsFunc project_points;
+ int minpts;
+} RansacModelInfo;
+
+#if ALLOW_TRANSLATION_MODELS
+static void project_points_translation(const double *mat, const 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++);
@@ -53,23 +69,11 @@
proj += stride_proj - 2;
}
}
+#endif // ALLOW_TRANSLATION_MODELS
-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) {
+static void project_points_affine(const double *mat, const 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++);
@@ -80,241 +84,137 @@
}
}
-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;
+#if ALLOW_TRANSLATION_MODELS
+static bool find_translation(int np, const double *pts1, const double *pts2,
+ double *params) {
+ double sumx = 0;
+ double sumy = 0;
- 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++);
+ for (int i = 0; i < np; ++i) {
+ double dx = *(pts2++);
+ double dy = *(pts2++);
+ double sx = *(pts1++);
+ double sy = *(pts1++);
sumx += dx - sx;
sumy += dy - sy;
}
- mat[0] = sumx / np;
- mat[1] = sumy / np;
- denormalize_translation_reorder(mat, T1, T2);
- return 0;
+
+ params[0] = sumx / np;
+ params[1] = sumy / np;
+ params[2] = 1;
+ params[3] = 0;
+ params[4] = 0;
+ params[5] = 1;
+ return true;
+}
+#endif // ALLOW_TRANSLATION_MODELS
+
+static bool find_rotzoom(int np, const double *pts1, const double *pts2,
+ double *params) {
+ const int n = 4; // Size of least-squares problem
+ double mat[4 * 4]; // Accumulator for A'A
+ double y[4]; // Accumulator for A'b
+ double x[4]; // Output vector
+ double a[4]; // Single row of A
+ double b; // Single element of b
+
+ least_squares_init(mat, y, n);
+ for (int i = 0; i < np; ++i) {
+ double dx = *(pts2++);
+ double dy = *(pts2++);
+ double sx = *(pts1++);
+ double sy = *(pts1++);
+
+ a[0] = sx;
+ a[1] = sy;
+ a[2] = 1;
+ a[3] = 0;
+ b = dx;
+ least_squares_accumulate(mat, y, a, b, n);
+
+ a[0] = sy;
+ a[1] = -sx;
+ a[2] = 0;
+ a[3] = 1;
+ b = dy;
+ least_squares_accumulate(mat, y, a, b, n);
+ }
+
+ if (!least_squares_solve(mat, y, x, n)) {
+ return false;
+ }
+
+ // Rearrange least squares result to form output model
+ params[0] = x[2];
+ params[1] = x[3];
+ params[2] = x[0];
+ params[3] = x[1];
+ params[4] = -params[3];
+ params[5] = params[2];
+
+ return true;
}
-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));
- if (a == NULL) return 1;
- double *b = a + np2 * 4;
- double *temp = b + np2;
- int i;
- double sx, sy, dx, dy;
+// TODO(rachelbarker): As the x and y equations are decoupled in find_affine(),
+// the least-squares problem can be split this into two 3-dimensional problems,
+// which should be faster to solve.
+static bool find_affine(int np, const double *pts1, const double *pts2,
+ double *params) {
+ const int n = 6; // Size of least-squares problem
+ double mat[6 * 6]; // Accumulator for A'A
+ double y[6]; // Accumulator for A'b
+ double x[6]; // Output vector
+ double a[6]; // Single row of A
+ double b; // Single element of b
- double T1[9], T2[9];
- normalize_homography(pts1, np, T1);
- normalize_homography(pts2, np, T2);
+ least_squares_init(mat, y, n);
+ for (int i = 0; i < np; ++i) {
+ double dx = *(pts2++);
+ double dy = *(pts2++);
+ double sx = *(pts1++);
+ double sy = *(pts1++);
- for (i = 0; i < np; ++i) {
- dx = *(pts2++);
- dy = *(pts2++);
- sx = *(pts1++);
- sy = *(pts1++);
+ a[0] = sx;
+ a[1] = sy;
+ a[2] = 0;
+ a[3] = 0;
+ a[4] = 1;
+ a[5] = 0;
+ b = dx;
+ least_squares_accumulate(mat, y, a, b, n);
- 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;
+ a[0] = 0;
+ a[1] = 0;
+ a[2] = sx;
+ a[3] = sy;
+ a[4] = 0;
+ a[5] = 1;
+ b = dy;
+ least_squares_accumulate(mat, y, a, b, n);
}
- if (!least_squares(4, a, np2, 4, b, temp, mat)) {
- aom_free(a);
- return 1;
+
+ if (!least_squares_solve(mat, y, x, n)) {
+ return false;
}
- denormalize_rotzoom_reorder(mat, T1, T2);
- aom_free(a);
- return 0;
+
+ // Rearrange least squares result to form output model
+ params[0] = mat[4];
+ params[1] = mat[5];
+ params[2] = mat[0];
+ params[3] = mat[1];
+ params[4] = mat[2];
+ params[5] = mat[3];
+
+ return true;
}
-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) {
+// 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 0;
+ if (minpts > npoints) return false;
indices[0] = ptr;
ptr = (ptr == npoints - 1 ? 0 : ptr + 1);
i = 1;
@@ -329,7 +229,7 @@
}
indices[i++] = ptr;
}
- return 1;
+ return true;
}
typedef struct {
@@ -350,8 +250,8 @@
return 0;
}
-static int is_better_motion(const RANSAC_MOTION *motion_a,
- const RANSAC_MOTION *motion_b) {
+static bool is_better_motion(const RANSAC_MOTION *motion_a,
+ const RANSAC_MOTION *motion_b) {
return compare_motions(motion_a, motion_b) < 0;
}
@@ -373,15 +273,14 @@
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) {
+// Returns true on success, false on error
+static bool ransac_internal(const Correspondence *matched_points, int npoints,
+ MotionModel *motion_models, int num_desired_motions,
+ const RansacModelInfo *model_info) {
int trial_count = 0;
int i = 0;
- int ret_val = 0;
+ int minpts = model_info->minpts;
+ bool ret_val = true;
unsigned int seed = (unsigned int)npoints;
@@ -401,13 +300,11 @@
// 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;
+ motion_models[i].num_inliers = 0;
}
if (npoints < minpts * MINPTS_MULTIPLIER || npoints == 0) {
- return 1;
+ return false;
}
points1 = (double *)aom_malloc(sizeof(*points1) * npoints * 2);
@@ -415,36 +312,31 @@
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_calloc(num_desired_motions, sizeof(RANSAC_MOTION));
- current_motion.inlier_indices =
- (int *)aom_malloc(sizeof(*current_motion.inlier_indices) * npoints);
- if (!(points1 && points2 && corners1 && corners2 && image1_coord && motions &&
- current_motion.inlier_indices)) {
- ret_val = 1;
- goto finish_ransac;
- }
+ 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);
- if (!motions[i].inlier_indices) {
- ret_val = 1;
- goto finish_ransac;
- }
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;
- cnp1 = corners1;
- cnp2 = corners2;
+ if (!(points1 && points2 && corners1 && corners2 && image1_coord && motions &&
+ current_motion.inlier_indices)) {
+ ret_val = false;
+ goto finish_ransac;
+ }
+
for (i = 0; i < npoints; ++i) {
- *(cnp1++) = *(matched_points++);
- *(cnp1++) = *(matched_points++);
- *(cnp2++) = *(matched_points++);
- *(cnp2++) = *(matched_points++);
+ corners1[2 * i + 0] = matched_points[i].x;
+ corners1[2 * i + 1] = matched_points[i].y;
+ corners2[2 * i + 0] = matched_points[i].rx;
+ corners2[2 * i + 1] = matched_points[i].ry;
}
while (MIN_TRIALS > trial_count) {
@@ -459,26 +351,28 @@
while (degenerate) {
num_degenerate_iter++;
if (!get_rand_indices(npoints, minpts, indices, &seed)) {
- ret_val = 1;
+ 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);
+ degenerate = model_info->is_degenerate(points1);
if (num_degenerate_iter > MAX_DEGENERATE_ITER) {
- ret_val = 1;
+ ret_val = false;
goto finish_ransac;
}
}
- if (find_transformation(minpts, points1, points2, params_this_motion)) {
+ if (!model_info->find_transformation(minpts, points1, points2,
+ params_this_motion)) {
trial_count++;
continue;
}
- projectpoints(params_this_motion, corners1, image1_coord, npoints, 2, 2);
+ model_info->project_points(params_this_motion, corners1, image1_coord,
+ npoints, 2, 2);
for (i = 0; i < npoints; ++i) {
double dx = image1_coord[i * 2] - corners2[i * 2];
@@ -526,19 +420,24 @@
// 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,
- motions[i].num_inliers);
+ num_inliers);
copy_points_at_indices(points2, corners2, motions[i].inlier_indices,
- motions[i].num_inliers);
+ num_inliers);
- find_transformation(motions[i].num_inliers, points1, points2,
- params_by_motion[i].params);
+ model_info->find_transformation(num_inliers, points1, points2,
+ motion_models[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;
+ // Populate inliers array
+ for (int j = 0; j < num_inliers; j++) {
+ int index = motions[i].inlier_indices[j];
+ const Correspondence *corr = &matched_points[index];
+ motion_models[i].inliers[2 * j + 0] = (int)rint(corr->x);
+ motion_models[i].inliers[2 * j + 1] = (int)rint(corr->y);
+ }
}
+ motion_models[i].num_inliers = motions[i].num_inliers;
}
finish_ransac:
@@ -558,277 +457,49 @@
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_calloc(num_desired_motions, sizeof(RANSAC_MOTION));
- current_motion.inlier_indices =
- (int *)aom_malloc(sizeof(*current_motion.inlier_indices) * npoints);
- if (!(points1 && points2 && corners1 && corners2 && image1_coord && motions &&
- current_motion.inlier_indices)) {
- ret_val = 1;
- goto finish_ransac;
- }
-
- for (i = 0; i < num_desired_motions; ++i) {
- motions[i].inlier_indices =
- (int *)aom_malloc(sizeof(*motions->inlier_indices) * npoints);
- if (!motions[i].inlier_indices) {
- ret_val = 1;
- goto finish_ransac;
- }
- clear_motion(motions + i, npoints);
- }
- clear_motion(¤t_motion, npoints);
-
- worst_kept_motion = motions;
-
- 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);
- if (motions) {
- 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 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 int is_degenerate_translation(double *p) {
+#if ALLOW_TRANSLATION_MODELS
+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;
}
+#endif // ALLOW_TRANSLATION_MODELS
-static int is_degenerate_affine(double *p) {
+static bool 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 const RansacModelInfo ransac_model_info[TRANS_TYPES] = {
+ // IDENTITY
+ { NULL, NULL, NULL, 0 },
+// TRANSLATION
+#if ALLOW_TRANSLATION_MODELS
+ { is_degenerate_translation, find_translation, project_points_translation,
+ 3 },
+#else
+ { NULL, NULL, NULL, 0 },
+#endif
+ // ROTZOOM
+ { is_degenerate_affine, find_rotzoom, project_points_affine, 3 },
+ // AFFINE
+ { is_degenerate_affine, find_affine, project_points_affine, 3 },
+};
-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);
-}
+// Returns true on success, false on error
+bool ransac(const Correspondence *matched_points, int npoints,
+ TransformationType type, MotionModel *motion_models,
+ int num_desired_motions) {
+#if ALLOW_TRANSLATION_MODELS
+ assert(type > IDENTITY && type < TRANS_TYPES);
+#else
+ assert(type > TRANSLATION && type < TRANS_TYPES);
+#endif // ALLOW_TRANSLATION_MODELS
-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;
- }
+ return ransac_internal(matched_points, npoints, motion_models,
+ num_desired_motions, &ransac_model_info[type]);
}
diff --git a/aom_dsp/flow_estimation/ransac.h b/aom_dsp/flow_estimation/ransac.h
index aa3a243..6047580 100644
--- a/aom_dsp/flow_estimation/ransac.h
+++ b/aom_dsp/flow_estimation/ransac.h
@@ -16,6 +16,7 @@
#include <stdlib.h>
#include <math.h>
#include <memory.h>
+#include <stdbool.h>
#include "aom_dsp/flow_estimation/flow_estimation.h"
@@ -23,14 +24,9 @@
extern "C" {
#endif
-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);
+bool ransac(const Correspondence *matched_points, int npoints,
+ TransformationType type, MotionModel *motion_models,
+ int num_desired_motions);
#ifdef __cplusplus
}
diff --git a/aom_dsp/mathutils.h b/aom_dsp/mathutils.h
index 22b0202..72572eb 100644
--- a/aom_dsp/mathutils.h
+++ b/aom_dsp/mathutils.h
@@ -63,32 +63,45 @@
// 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 i, j, k;
- double *scratch_ = NULL;
- double *AtA, *Atb;
- if (!scratch) {
- scratch_ = (double *)aom_malloc(sizeof(*scratch) * n * (n + 1));
- if (!scratch_) return 0;
- scratch = scratch_;
- }
- AtA = scratch;
- Atb = scratch + n * n;
+//
+// This process is split into three steps in order to avoid needing to
+// explicitly allocate the A matrix, which may be very large if there
+// are many equations to solve.
+//
+// The process for using this is (in pseudocode):
+//
+// Allocate mat (size n*n), y (size n), a (size n), x (size n)
+// least_squares_init(mat, y, n)
+// for each equation a . x = b {
+// least_squares_accumulate(mat, y, a, b, n)
+// }
+// least_squares_solve(mat, y, x, n)
+//
+// where:
+// * mat, y are accumulators for the values A'A and A'b respectively,
+// * a, b are the coefficients of each individual equation,
+// * x is the result vector
+// * and n is the problem size
+static INLINE void least_squares_init(double *mat, double *y, int n) {
+ memset(mat, 0, n * n * sizeof(double));
+ memset(y, 0, n * sizeof(double));
+}
- for (i = 0; i < n; ++i) {
- for (j = i; j < n; ++j) {
- AtA[i * n + j] = 0.0;
- for (k = 0; k < rows; ++k)
- AtA[i * n + j] += A[k * stride + i] * A[k * stride + j];
- AtA[j * n + i] = AtA[i * n + j];
+static INLINE void least_squares_accumulate(double *mat, double *y,
+ const double *a, double b, int n) {
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ mat[i * n + j] += a[i] * a[j];
}
- Atb[i] = 0;
- for (k = 0; k < rows; ++k) Atb[i] += A[k * stride + i] * b[k];
}
- int ret = linsolve(n, AtA, n, Atb, x);
- aom_free(scratch_);
- return ret;
+ for (int i = 0; i < n; i++) {
+ y[i] += a[i] * b;
+ }
+}
+
+static INLINE int least_squares_solve(double *mat, double *y, double *x,
+ int n) {
+ return linsolve(n, mat, n, y, x);
}
// Matrix multiply
diff --git a/av1/encoder/ethread.c b/av1/encoder/ethread.c
index 5f2f6d9..f94e675 100644
--- a/av1/encoder/ethread.c
+++ b/av1/encoder/ethread.c
@@ -2220,7 +2220,7 @@
static AOM_INLINE void init_gm_thread_data(
const GlobalMotionInfo *gm_info, GlobalMotionThreadData *thread_data) {
for (int m = 0; m < RANSAC_NUM_MOTIONS; m++) {
- MotionModel motion_params = thread_data->params_by_motion[m];
+ MotionModel motion_params = thread_data->motion_models[m];
av1_zero(motion_params.params);
motion_params.num_inliers = 0;
}
@@ -2277,7 +2277,7 @@
// 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_thread_data->params_by_motion,
+ cpi, gm_info->ref_buf, ref_buf_idx, gm_thread_data->motion_models,
gm_thread_data->segment_map, gm_info->segment_map_w,
gm_info->segment_map_h);
@@ -2357,7 +2357,7 @@
aom_free(thread_data->segment_map);
for (int m = 0; m < RANSAC_NUM_MOTIONS; m++)
- aom_free(thread_data->params_by_motion[m].inliers);
+ aom_free(thread_data->motion_models[m].inliers);
}
aom_free(gm_sync_data->thread_data);
}
@@ -2386,8 +2386,8 @@
for (int m = 0; m < RANSAC_NUM_MOTIONS; m++) {
CHECK_MEM_ERROR(
- cm, thread_data->params_by_motion[m].inliers,
- aom_malloc(sizeof(*thread_data->params_by_motion[m].inliers) * 2 *
+ cm, thread_data->motion_models[m].inliers,
+ aom_malloc(sizeof(*thread_data->motion_models[m].inliers) * 2 *
MAX_CORNERS));
}
}
diff --git a/av1/encoder/global_motion.h b/av1/encoder/global_motion.h
index 4fa3253..478be8e 100644
--- a/av1/encoder/global_motion.h
+++ b/av1/encoder/global_motion.h
@@ -34,9 +34,9 @@
typedef struct {
// Array of structure which holds the global motion parameters for a given
- // motion model. params_by_motion[i] holds the parameters for a given motion
+ // motion model. motion_models[i] holds the parameters for a given motion
// model for the ith ransac motion.
- MotionModel params_by_motion[RANSAC_NUM_MOTIONS];
+ MotionModel motion_models[RANSAC_NUM_MOTIONS];
// Pointer to hold inliers from motion model.
uint8_t *segment_map;
diff --git a/av1/encoder/global_motion_facade.c b/av1/encoder/global_motion_facade.c
index 9c9f9f3..ee3b525 100644
--- a/av1/encoder/global_motion_facade.c
+++ b/av1/encoder/global_motion_facade.c
@@ -82,9 +82,8 @@
// different motion models and finds the best.
static AOM_INLINE void compute_global_motion_for_ref_frame(
AV1_COMP *cpi, YV12_BUFFER_CONFIG *ref_buf[REF_FRAMES], int frame,
- MotionModel *params_by_motion, uint8_t *segment_map,
- const int segment_map_w, const int segment_map_h,
- const WarpedMotionParams *ref_params) {
+ MotionModel *motion_models, 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;
@@ -100,7 +99,6 @@
// 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;
int bit_depth = cpi->common.seq_params->bit_depth;
@@ -120,19 +118,19 @@
int64_t best_warp_error = INT64_MAX;
// Initially set all params to identity.
for (i = 0; i < RANSAC_NUM_MOTIONS; ++i) {
- memcpy(params_by_motion[i].params, kIdentityParams,
- (MAX_PARAMDIM - 1) * sizeof(*(params_by_motion[i].params)));
- params_by_motion[i].num_inliers = 0;
+ memcpy(motion_models[i].params, kIdentityParams,
+ (MAX_PARAMDIM - 1) * sizeof(*(motion_models[i].params)));
+ motion_models[i].num_inliers = 0;
}
aom_compute_global_motion(model, cpi->source, ref_buf[frame], bit_depth,
- gm_estimation_type, inliers_by_motion,
- params_by_motion, RANSAC_NUM_MOTIONS);
+ gm_estimation_type, motion_models,
+ 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 (motion_models[i].num_inliers == 0) continue;
- params_this_motion = params_by_motion[i].params;
+ params_this_motion = motion_models[i].params;
av1_convert_model_to_params(params_this_motion, &tmp_wm_params);
// Work around a bug in the AV1 specification
@@ -147,8 +145,8 @@
if (tmp_wm_params.wmtype != IDENTITY) {
av1_compute_feature_segmentation_map(
- segment_map, segment_map_w, segment_map_h,
- params_by_motion[i].inliers, params_by_motion[i].num_inliers);
+ segment_map, segment_map_w, segment_map_h, motion_models[i].inliers,
+ motion_models[i].num_inliers);
ref_frame_error = av1_segmented_frame_error(
is_cur_buf_hbd(xd), xd->bd, ref_buf[frame]->y_buffer,
@@ -220,14 +218,14 @@
// Computes global motion for the given reference frame.
void av1_compute_gm_for_valid_ref_frames(
AV1_COMP *cpi, YV12_BUFFER_CONFIG *ref_buf[REF_FRAMES], int frame,
- MotionModel *params_by_motion, uint8_t *segment_map, int segment_map_w,
+ MotionModel *motion_models, uint8_t *segment_map, int segment_map_w,
int segment_map_h) {
AV1_COMMON *const cm = &cpi->common;
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, params_by_motion,
+ compute_global_motion_for_ref_frame(cpi, ref_buf, frame, motion_models,
segment_map, segment_map_w, segment_map_h,
ref_params);
}
@@ -236,16 +234,16 @@
static AOM_INLINE void compute_global_motion_for_references(
AV1_COMP *cpi, YV12_BUFFER_CONFIG *ref_buf[REF_FRAMES],
FrameDistPair reference_frame[REF_FRAMES - 1], int num_ref_frames,
- MotionModel *params_by_motion, uint8_t *segment_map,
- const int segment_map_w, const int segment_map_h) {
+ MotionModel *motion_models, 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,
- params_by_motion, segment_map,
- segment_map_w, segment_map_h);
+ av1_compute_gm_for_valid_ref_frames(cpi, ref_buf, ref_frame, motion_models,
+ 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
@@ -372,26 +370,26 @@
}
// Deallocates segment_map and inliers.
-static AOM_INLINE void dealloc_global_motion_data(MotionModel *params_by_motion,
+static AOM_INLINE void dealloc_global_motion_data(MotionModel *motion_models,
uint8_t *segment_map) {
aom_free(segment_map);
for (int m = 0; m < RANSAC_NUM_MOTIONS; m++) {
- aom_free(params_by_motion[m].inliers);
+ aom_free(motion_models[m].inliers);
}
}
// Allocates and initializes memory for segment_map and MotionModel.
-static AOM_INLINE bool alloc_global_motion_data(MotionModel *params_by_motion,
+static AOM_INLINE bool alloc_global_motion_data(MotionModel *motion_models,
uint8_t **segment_map,
const int segment_map_w,
const int segment_map_h) {
- av1_zero_array(params_by_motion, RANSAC_NUM_MOTIONS);
+ av1_zero_array(motion_models, RANSAC_NUM_MOTIONS);
for (int m = 0; m < RANSAC_NUM_MOTIONS; m++) {
- params_by_motion[m].inliers =
- aom_malloc(sizeof(*(params_by_motion[m].inliers)) * 2 * MAX_CORNERS);
- if (!params_by_motion[m].inliers) {
- dealloc_global_motion_data(params_by_motion, NULL);
+ motion_models[m].inliers =
+ aom_malloc(sizeof(*(motion_models[m].inliers)) * 2 * MAX_CORNERS);
+ if (!motion_models[m].inliers) {
+ dealloc_global_motion_data(motion_models, NULL);
return false;
}
}
@@ -399,7 +397,7 @@
*segment_map = (uint8_t *)aom_calloc(segment_map_w * segment_map_h,
sizeof(*segment_map));
if (!*segment_map) {
- dealloc_global_motion_data(params_by_motion, NULL);
+ dealloc_global_motion_data(motion_models, NULL);
return false;
}
return true;
@@ -437,11 +435,11 @@
// Computes global motion w.r.t. valid reference frames.
static AOM_INLINE void global_motion_estimation(AV1_COMP *cpi) {
GlobalMotionInfo *const gm_info = &cpi->gm_info;
- MotionModel params_by_motion[RANSAC_NUM_MOTIONS];
+ MotionModel motion_models[RANSAC_NUM_MOTIONS];
uint8_t *segment_map = NULL;
- alloc_global_motion_data(params_by_motion, &segment_map,
- gm_info->segment_map_w, gm_info->segment_map_h);
+ alloc_global_motion_data(motion_models, &segment_map, gm_info->segment_map_w,
+ gm_info->segment_map_h);
// Compute global motion w.r.t. past reference frames and future reference
// frames
@@ -449,11 +447,11 @@
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], params_by_motion, segment_map,
+ gm_info->num_ref_frames[dir], motion_models, segment_map,
gm_info->segment_map_w, gm_info->segment_map_h);
}
- dealloc_global_motion_data(params_by_motion, segment_map);
+ dealloc_global_motion_data(motion_models, segment_map);
}
// Global motion estimation for the current frame is computed.This computation
diff --git a/av1/encoder/global_motion_facade.h b/av1/encoder/global_motion_facade.h
index 9d33153..dfdedf7 100644
--- a/av1/encoder/global_motion_facade.h
+++ b/av1/encoder/global_motion_facade.h
@@ -20,7 +20,7 @@
void av1_compute_gm_for_valid_ref_frames(
AV1_COMP *cpi, YV12_BUFFER_CONFIG *ref_buf[REF_FRAMES], int frame,
- MotionModel *params_by_motion, uint8_t *segment_map, int segment_map_w,
+ MotionModel *motion_models, uint8_t *segment_map, int segment_map_w,
int segment_map_h);
void av1_compute_global_motion_facade(struct AV1_COMP *cpi);
#ifdef __cplusplus
