blob: 81c5f2c6256e9114f13297d2ff6e75c0f5cd6cb2 [file] [log] [blame]
Sarah Parker4dc0f1b2016-08-09 17:40:53 -07001/*
Wan-Teh Chang0289a9f2022-12-14 09:33:45 -08002 * Copyright (c) 2016, Alliance for Open Media. All rights reserved
Sarah Parker4dc0f1b2016-08-09 17:40:53 -07003 *
Wan-Teh Chang0289a9f2022-12-14 09:33:45 -08004 * This source code is subject to the terms of the BSD 2 Clause License and
5 * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
6 * was not distributed with this source code in the LICENSE file, you can
7 * obtain it at www.aomedia.org/license/software. If the Alliance for Open
8 * Media Patent License 1.0 was not distributed with this source code in the
9 * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
Sarah Parker4dc0f1b2016-08-09 17:40:53 -070010 */
Rachel Barker53e3a692022-11-30 15:45:11 +000011
Sarah Parker4dc0f1b2016-08-09 17:40:53 -070012#include <memory.h>
13#include <math.h>
14#include <time.h>
15#include <stdio.h>
Rachel Barker738c4f42023-01-11 17:36:18 +000016#include <stdbool.h>
Rachel Barkera02796a2023-03-02 20:06:47 +000017#include <string.h>
Sarah Parker4dc0f1b2016-08-09 17:40:53 -070018#include <assert.h>
19
Rachel Barker53e3a692022-11-30 15:45:11 +000020#include "aom_dsp/flow_estimation/ransac.h"
Yaowu Xu9a0d1b72021-07-13 09:56:50 -070021#include "aom_dsp/mathutils.h"
Rachel Barker738c4f42023-01-11 17:36:18 +000022#include "aom_mem/aom_mem.h"
Rachel Barker53e3a692022-11-30 15:45:11 +000023
24// TODO(rachelbarker): Remove dependence on code in av1/encoder/
Alex Converse9d068c12017-08-03 11:48:19 -070025#include "av1/encoder/random.h"
Sarah Parker4dc0f1b2016-08-09 17:40:53 -070026
Sarah Parker4dc0f1b2016-08-09 17:40:53 -070027#define MAX_MINPTS 4
Sarah Parker4dc0f1b2016-08-09 17:40:53 -070028#define MINPTS_MULTIPLIER 5
29
Debargha Mukherjee3b16e002018-09-19 13:44:05 -070030#define INLIER_THRESHOLD 1.25
Rachel Barker36042692023-01-11 18:09:43 +000031#define INLIER_THRESHOLD_SQUARED (INLIER_THRESHOLD * INLIER_THRESHOLD)
Rachel Barker41f1bb02023-02-13 16:22:38 +000032#define NUM_TRIALS 20
Debargha Mukherjee5dfa9302017-02-10 05:00:08 -080033
Rachel Barker738c4f42023-01-11 17:36:18 +000034// Flag to enable functions for finding TRANSLATION type models.
35//
36// These modes are not considered currently due to a spec bug (see comments
37// in gm_get_motion_vector() in av1/common/mv.h). Thus we don't need to compile
38// the corresponding search functions, but it is nice to keep the source around
39// but disabled, for completeness.
40#define ALLOW_TRANSLATION_MODELS 0
41
Sarah Parker4dc0f1b2016-08-09 17:40:53 -070042////////////////////////////////////////////////////////////////////////////////
43// ransac
Rachel Barker738c4f42023-01-11 17:36:18 +000044typedef bool (*IsDegenerateFunc)(double *p);
45typedef bool (*FindTransformationFunc)(int points, const double *points1,
46 const double *points2, double *params);
47typedef void (*ProjectPointsFunc)(const double *mat, const double *points,
48 double *proj, int n, int stride_points,
49 int stride_proj);
Sarah Parker97fa6da2016-09-23 11:17:27 -070050
Rachel Barker738c4f42023-01-11 17:36:18 +000051// vtable-like structure which stores all of the information needed by RANSAC
52// for a particular model type
53typedef struct {
54 IsDegenerateFunc is_degenerate;
55 FindTransformationFunc find_transformation;
56 ProjectPointsFunc project_points;
57 int minpts;
58} RansacModelInfo;
59
60#if ALLOW_TRANSLATION_MODELS
61static void project_points_translation(const double *mat, const double *points,
62 double *proj, int n, int stride_points,
63 int stride_proj) {
Sarah Parker97fa6da2016-09-23 11:17:27 -070064 int i;
65 for (i = 0; i < n; ++i) {
66 const double x = *(points++), y = *(points++);
Debargha Mukherjee8db4c772016-11-07 12:54:21 -080067 *(proj++) = x + mat[0];
68 *(proj++) = y + mat[1];
Sarah Parker97fa6da2016-09-23 11:17:27 -070069 points += stride_points - 2;
70 proj += stride_proj - 2;
71 }
72}
Rachel Barker738c4f42023-01-11 17:36:18 +000073#endif // ALLOW_TRANSLATION_MODELS
Sarah Parker97fa6da2016-09-23 11:17:27 -070074
Rachel Barker738c4f42023-01-11 17:36:18 +000075static void project_points_affine(const double *mat, const double *points,
76 double *proj, int n, int stride_points,
77 int stride_proj) {
Sarah Parker97fa6da2016-09-23 11:17:27 -070078 int i;
79 for (i = 0; i < n; ++i) {
80 const double x = *(points++), y = *(points++);
Debargha Mukherjee8db4c772016-11-07 12:54:21 -080081 *(proj++) = mat[2] * x + mat[3] * y + mat[0];
82 *(proj++) = mat[4] * x + mat[5] * y + mat[1];
Sarah Parker97fa6da2016-09-23 11:17:27 -070083 points += stride_points - 2;
84 proj += stride_proj - 2;
85 }
86}
87
Rachel Barker738c4f42023-01-11 17:36:18 +000088#if ALLOW_TRANSLATION_MODELS
89static bool find_translation(int np, const double *pts1, const double *pts2,
90 double *params) {
91 double sumx = 0;
92 double sumy = 0;
Debargha Mukherjeeafe7c5f2017-06-22 13:47:23 -070093
Rachel Barker738c4f42023-01-11 17:36:18 +000094 for (int i = 0; i < np; ++i) {
95 double dx = *(pts2++);
96 double dy = *(pts2++);
97 double sx = *(pts1++);
98 double sy = *(pts1++);
Debargha Mukherjeee6eb3b52017-02-26 08:50:56 -080099
100 sumx += dx - sx;
101 sumy += dy - sy;
102 }
Rachel Barker738c4f42023-01-11 17:36:18 +0000103
104 params[0] = sumx / np;
105 params[1] = sumy / np;
106 params[2] = 1;
107 params[3] = 0;
108 params[4] = 0;
109 params[5] = 1;
110 return true;
111}
112#endif // ALLOW_TRANSLATION_MODELS
113
114static bool find_rotzoom(int np, const double *pts1, const double *pts2,
115 double *params) {
116 const int n = 4; // Size of least-squares problem
117 double mat[4 * 4]; // Accumulator for A'A
118 double y[4]; // Accumulator for A'b
Rachel Barker738c4f42023-01-11 17:36:18 +0000119 double a[4]; // Single row of A
120 double b; // Single element of b
121
122 least_squares_init(mat, y, n);
123 for (int i = 0; i < np; ++i) {
124 double dx = *(pts2++);
125 double dy = *(pts2++);
126 double sx = *(pts1++);
127 double sy = *(pts1++);
128
Rachel Barker24ec16a2023-02-09 13:47:35 +0000129 a[0] = 1;
130 a[1] = 0;
131 a[2] = sx;
132 a[3] = sy;
Rachel Barker738c4f42023-01-11 17:36:18 +0000133 b = dx;
134 least_squares_accumulate(mat, y, a, b, n);
135
Rachel Barker24ec16a2023-02-09 13:47:35 +0000136 a[0] = 0;
137 a[1] = 1;
138 a[2] = sy;
139 a[3] = -sx;
Rachel Barker738c4f42023-01-11 17:36:18 +0000140 b = dy;
141 least_squares_accumulate(mat, y, a, b, n);
142 }
143
Rachel Barker24ec16a2023-02-09 13:47:35 +0000144 // Fill in params[0] .. params[3] with output model
145 if (!least_squares_solve(mat, y, params, n)) {
Rachel Barker738c4f42023-01-11 17:36:18 +0000146 return false;
147 }
148
Rachel Barker24ec16a2023-02-09 13:47:35 +0000149 // Fill in remaining parameters
Rachel Barker738c4f42023-01-11 17:36:18 +0000150 params[4] = -params[3];
151 params[5] = params[2];
152
153 return true;
Debargha Mukherjeee6eb3b52017-02-26 08:50:56 -0800154}
155
Rachel Barker738c4f42023-01-11 17:36:18 +0000156static bool find_affine(int np, const double *pts1, const double *pts2,
157 double *params) {
Rachel Barker24ec16a2023-02-09 13:47:35 +0000158 // Note: The least squares problem for affine models is 6-dimensional,
159 // but it splits into two independent 3-dimensional subproblems.
160 // Solving these two subproblems separately and recombining at the end
161 // results in less total computation than solving the 6-dimensional
162 // problem directly.
163 //
164 // The two subproblems correspond to all the parameters which contribute
165 // to the x output of the model, and all the parameters which contribute
166 // to the y output, respectively.
Debargha Mukherjeee6eb3b52017-02-26 08:50:56 -0800167
Rachel Barker24ec16a2023-02-09 13:47:35 +0000168 const int n = 3; // Size of each least-squares problem
169 double mat[2][3 * 3]; // Accumulator for A'A
170 double y[2][3]; // Accumulator for A'b
171 double x[2][3]; // Output vector
172 double a[2][3]; // Single row of A
173 double b[2]; // Single element of b
174
175 least_squares_init(mat[0], y[0], n);
176 least_squares_init(mat[1], y[1], n);
Rachel Barker738c4f42023-01-11 17:36:18 +0000177 for (int i = 0; i < np; ++i) {
178 double dx = *(pts2++);
179 double dy = *(pts2++);
180 double sx = *(pts1++);
181 double sy = *(pts1++);
Debargha Mukherjeee6eb3b52017-02-26 08:50:56 -0800182
Rachel Barker24ec16a2023-02-09 13:47:35 +0000183 a[0][0] = 1;
184 a[0][1] = sx;
185 a[0][2] = sy;
186 b[0] = dx;
187 least_squares_accumulate(mat[0], y[0], a[0], b[0], n);
Debargha Mukherjeee6eb3b52017-02-26 08:50:56 -0800188
Rachel Barker24ec16a2023-02-09 13:47:35 +0000189 a[1][0] = 1;
190 a[1][1] = sx;
191 a[1][2] = sy;
192 b[1] = dy;
193 least_squares_accumulate(mat[1], y[1], a[1], b[1], n);
Debargha Mukherjeee6eb3b52017-02-26 08:50:56 -0800194 }
Rachel Barker738c4f42023-01-11 17:36:18 +0000195
Rachel Barker24ec16a2023-02-09 13:47:35 +0000196 if (!least_squares_solve(mat[0], y[0], x[0], n)) {
197 return false;
198 }
199 if (!least_squares_solve(mat[1], y[1], x[1], n)) {
Rachel Barker738c4f42023-01-11 17:36:18 +0000200 return false;
Debargha Mukherjeee6eb3b52017-02-26 08:50:56 -0800201 }
Rachel Barker738c4f42023-01-11 17:36:18 +0000202
203 // Rearrange least squares result to form output model
Rachel Barker24ec16a2023-02-09 13:47:35 +0000204 params[0] = x[0][0];
205 params[1] = x[1][0];
206 params[2] = x[0][1];
207 params[3] = x[0][2];
208 params[4] = x[1][1];
209 params[5] = x[1][2];
Rachel Barker738c4f42023-01-11 17:36:18 +0000210
211 return true;
Debargha Mukherjeee6eb3b52017-02-26 08:50:56 -0800212}
213
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500214typedef struct {
215 int num_inliers;
Rachel Barker36042692023-01-11 18:09:43 +0000216 double sse; // Sum of squared errors of inliers
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500217 int *inlier_indices;
218} RANSAC_MOTION;
219
220// Return -1 if 'a' is a better motion, 1 if 'b' is better, 0 otherwise.
221static int compare_motions(const void *arg_a, const void *arg_b) {
222 const RANSAC_MOTION *motion_a = (RANSAC_MOTION *)arg_a;
223 const RANSAC_MOTION *motion_b = (RANSAC_MOTION *)arg_b;
224
225 if (motion_a->num_inliers > motion_b->num_inliers) return -1;
226 if (motion_a->num_inliers < motion_b->num_inliers) return 1;
Rachel Barker36042692023-01-11 18:09:43 +0000227 if (motion_a->sse < motion_b->sse) return -1;
228 if (motion_a->sse > motion_b->sse) return 1;
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500229 return 0;
230}
231
Rachel Barker738c4f42023-01-11 17:36:18 +0000232static bool is_better_motion(const RANSAC_MOTION *motion_a,
233 const RANSAC_MOTION *motion_b) {
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500234 return compare_motions(motion_a, motion_b) < 0;
235}
236
237static void copy_points_at_indices(double *dest, const double *src,
238 const int *indices, int num_points) {
239 for (int i = 0; i < num_points; ++i) {
240 const int index = indices[i];
241 dest[i * 2] = src[index * 2];
242 dest[i * 2 + 1] = src[index * 2 + 1];
243 }
244}
245
Rachel Barker738c4f42023-01-11 17:36:18 +0000246// Returns true on success, false on error
247static bool ransac_internal(const Correspondence *matched_points, int npoints,
248 MotionModel *motion_models, int num_desired_motions,
249 const RansacModelInfo *model_info) {
Rachel Barkera02796a2023-03-02 20:06:47 +0000250 assert(npoints >= 0);
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500251 int i = 0;
Rachel Barker738c4f42023-01-11 17:36:18 +0000252 int minpts = model_info->minpts;
253 bool ret_val = true;
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500254
Sarah Parkerefa65822016-10-11 12:29:07 -0700255 unsigned int seed = (unsigned int)npoints;
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700256
Sarah Parkerefa65822016-10-11 12:29:07 -0700257 int indices[MAX_MINPTS] = { 0 };
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700258
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500259 double *points1, *points2;
260 double *corners1, *corners2;
Rachel Barker8e997222023-03-02 15:12:41 +0000261 double *projected_corners;
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700262
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500263 // Store information for the num_desired_motions best transformations found
264 // and the worst motion among them, as well as the motion currently under
265 // consideration.
266 RANSAC_MOTION *motions, *worst_kept_motion = NULL;
267 RANSAC_MOTION current_motion;
268
269 // Store the parameters and the indices of the inlier points for the motion
270 // currently under consideration.
271 double params_this_motion[MAX_PARAMDIM];
272
Sarah Parkerfa75ae02016-10-26 12:48:01 -0700273 if (npoints < minpts * MINPTS_MULTIPLIER || npoints == 0) {
Rachel Barker738c4f42023-01-11 17:36:18 +0000274 return false;
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700275 }
276
Rachel Barkera02796a2023-03-02 20:06:47 +0000277 int min_inliers = AOMMAX((int)(MIN_INLIER_PROB * npoints), minpts);
278
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500279 points1 = (double *)aom_malloc(sizeof(*points1) * npoints * 2);
280 points2 = (double *)aom_malloc(sizeof(*points2) * npoints * 2);
Sarah Parkerf9a961c2016-09-06 11:25:04 -0700281 corners1 = (double *)aom_malloc(sizeof(*corners1) * npoints * 2);
Sarah Parkerf9a961c2016-09-06 11:25:04 -0700282 corners2 = (double *)aom_malloc(sizeof(*corners2) * npoints * 2);
Rachel Barker8e997222023-03-02 15:12:41 +0000283 projected_corners =
284 (double *)aom_malloc(sizeof(*projected_corners) * npoints * 2);
Rachel Barker738c4f42023-01-11 17:36:18 +0000285 motions =
Rachel Barkerdfd80462023-03-02 20:49:41 +0000286 (RANSAC_MOTION *)aom_calloc(num_desired_motions, sizeof(RANSAC_MOTION));
James Zern9a0dd4b2022-05-17 18:20:49 -0700287
Rachel Barker79fc28c2023-03-02 19:41:08 +0000288 // Allocate one large buffer which will be carved up to store the inlier
289 // indices for the current motion plus the num_desired_motions many
290 // output models
291 // This allows us to keep the allocation/deallocation logic simple, without
292 // having to (for example) check that `motions` is non-null before allocating
293 // the inlier arrays
294 int *inlier_buffer = (int *)aom_malloc(sizeof(*inlier_buffer) * npoints *
295 (num_desired_motions + 1));
James Zern9a0dd4b2022-05-17 18:20:49 -0700296
Rachel Barker8e997222023-03-02 15:12:41 +0000297 if (!(points1 && points2 && corners1 && corners2 && projected_corners &&
Rachel Barker79fc28c2023-03-02 19:41:08 +0000298 motions && inlier_buffer)) {
Rachel Barker738c4f42023-01-11 17:36:18 +0000299 ret_val = false;
300 goto finish_ransac;
301 }
302
Rachel Barker79fc28c2023-03-02 19:41:08 +0000303 // Once all our allocations are known-good, we can fill in our structures
304 worst_kept_motion = motions;
305
306 for (i = 0; i < num_desired_motions; ++i) {
307 motions[i].inlier_indices = inlier_buffer + i * npoints;
Rachel Barker79fc28c2023-03-02 19:41:08 +0000308 }
Rachel Barkerdfd80462023-03-02 20:49:41 +0000309 memset(&current_motion, 0, sizeof(current_motion));
Rachel Barker79fc28c2023-03-02 19:41:08 +0000310 current_motion.inlier_indices = inlier_buffer + num_desired_motions * npoints;
Rachel Barker79fc28c2023-03-02 19:41:08 +0000311
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500312 for (i = 0; i < npoints; ++i) {
Rachel Barker738c4f42023-01-11 17:36:18 +0000313 corners1[2 * i + 0] = matched_points[i].x;
314 corners1[2 * i + 1] = matched_points[i].y;
315 corners2[2 * i + 0] = matched_points[i].rx;
316 corners2[2 * i + 1] = matched_points[i].ry;
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700317 }
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700318
Rachel Barker41f1bb02023-02-13 16:22:38 +0000319 for (int trial_count = 0; trial_count < NUM_TRIALS; trial_count++) {
Rachel Barkerfebe1a02023-03-07 15:10:17 +0000320 lcg_pick(npoints, minpts, indices, &seed);
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500321
Rachel Barker41f1bb02023-02-13 16:22:38 +0000322 copy_points_at_indices(points1, corners1, indices, minpts);
323 copy_points_at_indices(points2, corners2, indices, minpts);
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500324
Rachel Barker41f1bb02023-02-13 16:22:38 +0000325 if (model_info->is_degenerate(points1)) {
326 continue;
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700327 }
328
Rachel Barker738c4f42023-01-11 17:36:18 +0000329 if (!model_info->find_transformation(minpts, points1, points2,
330 params_this_motion)) {
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700331 continue;
332 }
333
Rachel Barker8e997222023-03-02 15:12:41 +0000334 model_info->project_points(params_this_motion, corners1, projected_corners,
Rachel Barker738c4f42023-01-11 17:36:18 +0000335 npoints, 2, 2);
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700336
Rachel Barker41f1bb02023-02-13 16:22:38 +0000337 current_motion.num_inliers = 0;
Rachel Barker36042692023-01-11 18:09:43 +0000338 double sse = 0.0;
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700339 for (i = 0; i < npoints; ++i) {
Rachel Barker8e997222023-03-02 15:12:41 +0000340 double dx = projected_corners[i * 2] - corners2[i * 2];
341 double dy = projected_corners[i * 2 + 1] - corners2[i * 2 + 1];
Rachel Barker36042692023-01-11 18:09:43 +0000342 double squared_error = dx * dx + dy * dy;
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700343
Rachel Barker36042692023-01-11 18:09:43 +0000344 if (squared_error < INLIER_THRESHOLD_SQUARED) {
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500345 current_motion.inlier_indices[current_motion.num_inliers++] = i;
Rachel Barker36042692023-01-11 18:09:43 +0000346 sse += squared_error;
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700347 }
348 }
349
Rachel Barkera02796a2023-03-02 20:06:47 +0000350 if (current_motion.num_inliers < min_inliers) {
Rachel Barker41f1bb02023-02-13 16:22:38 +0000351 // Reject models with too few inliers
352 continue;
353 }
354
355 current_motion.sse = sse;
356 if (is_better_motion(&current_motion, worst_kept_motion)) {
357 // This motion is better than the worst currently kept motion. Remember
358 // the inlier points and sse. The parameters for each kept motion
359 // will be recomputed later using only the inliers.
360 worst_kept_motion->num_inliers = current_motion.num_inliers;
361 worst_kept_motion->sse = current_motion.sse;
Rachel Barkerdfd80462023-03-02 20:49:41 +0000362
363 // Rather than copying the (potentially many) inlier indices from
364 // current_motion.inlier_indices to worst_kept_motion->inlier_indices,
365 // we can swap the underlying pointers.
366 //
367 // This is okay because the next time current_motion.inlier_indices
368 // is used will be in the next trial, where we ignore its previous
369 // contents anyway. And both arrays will be deallocated together at the
370 // end of this function, so there are no lifetime issues.
371 int *tmp = worst_kept_motion->inlier_indices;
372 worst_kept_motion->inlier_indices = current_motion.inlier_indices;
373 current_motion.inlier_indices = tmp;
374
Rachel Barker41f1bb02023-02-13 16:22:38 +0000375 // Determine the new worst kept motion and its num_inliers and sse.
376 for (i = 0; i < num_desired_motions; ++i) {
377 if (is_better_motion(worst_kept_motion, &motions[i])) {
378 worst_kept_motion = &motions[i];
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500379 }
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700380 }
381 }
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700382 }
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500383
384 // Sort the motions, best first.
385 qsort(motions, num_desired_motions, sizeof(RANSAC_MOTION), compare_motions);
386
387 // Recompute the motions using only the inliers.
388 for (i = 0; i < num_desired_motions; ++i) {
Rachel Barker41f1bb02023-02-13 16:22:38 +0000389 int num_inliers = motions[i].num_inliers;
390 if (num_inliers > 0) {
391 assert(num_inliers >= minpts);
392
Debargha Mukherjeeafe7c5f2017-06-22 13:47:23 -0700393 copy_points_at_indices(points1, corners1, motions[i].inlier_indices,
Rachel Barker738c4f42023-01-11 17:36:18 +0000394 num_inliers);
Debargha Mukherjeeafe7c5f2017-06-22 13:47:23 -0700395 copy_points_at_indices(points2, corners2, motions[i].inlier_indices,
Rachel Barker738c4f42023-01-11 17:36:18 +0000396 num_inliers);
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500397
Rachel Barkera02796a2023-03-02 20:06:47 +0000398 if (!model_info->find_transformation(num_inliers, points1, points2,
399 motion_models[i].params)) {
400 // In the unlikely event that this model fitting fails,
401 // we don't have a good fallback. So just clear the output
402 // model and move on
403 memcpy(motion_models[i].params, kIdentityParams,
Rachel Barker0d9df422023-03-07 20:51:55 +0000404 MAX_PARAMDIM * sizeof(*(motion_models[i].params)));
Rachel Barkera02796a2023-03-02 20:06:47 +0000405 motion_models[i].num_inliers = 0;
406 continue;
407 }
Sarah Parker0c74c852019-05-09 15:15:32 -0700408
Rachel Barker738c4f42023-01-11 17:36:18 +0000409 // Populate inliers array
410 for (int j = 0; j < num_inliers; j++) {
411 int index = motions[i].inlier_indices[j];
412 const Correspondence *corr = &matched_points[index];
413 motion_models[i].inliers[2 * j + 0] = (int)rint(corr->x);
414 motion_models[i].inliers[2 * j + 1] = (int)rint(corr->y);
415 }
Rachel Barkera02796a2023-03-02 20:06:47 +0000416 motion_models[i].num_inliers = num_inliers;
417 } else {
418 memcpy(motion_models[i].params, kIdentityParams,
Rachel Barker0d9df422023-03-07 20:51:55 +0000419 MAX_PARAMDIM * sizeof(*(motion_models[i].params)));
Rachel Barkera02796a2023-03-02 20:06:47 +0000420 motion_models[i].num_inliers = 0;
Debargha Mukherjeeafe7c5f2017-06-22 13:47:23 -0700421 }
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500422 }
423
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700424finish_ransac:
Rachel Barker79fc28c2023-03-02 19:41:08 +0000425 aom_free(inlier_buffer);
426 aom_free(motions);
Rachel Barker8e997222023-03-02 15:12:41 +0000427 aom_free(projected_corners);
Rachel Barker79fc28c2023-03-02 19:41:08 +0000428 aom_free(corners2);
429 aom_free(corners1);
430 aom_free(points2);
431 aom_free(points1);
emilkeyder@google.comf3477632017-03-01 16:29:14 -0500432
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700433 return ret_val;
434}
435
Rachel Barker738c4f42023-01-11 17:36:18 +0000436static bool is_collinear3(double *p1, double *p2, double *p3) {
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700437 static const double collinear_eps = 1e-3;
438 const double v =
439 (p2[0] - p1[0]) * (p3[1] - p1[1]) - (p2[1] - p1[1]) * (p3[0] - p1[0]);
440 return fabs(v) < collinear_eps;
441}
442
Rachel Barker738c4f42023-01-11 17:36:18 +0000443#if ALLOW_TRANSLATION_MODELS
444static bool is_degenerate_translation(double *p) {
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700445 return (p[0] - p[2]) * (p[0] - p[2]) + (p[1] - p[3]) * (p[1] - p[3]) <= 2;
446}
Rachel Barker738c4f42023-01-11 17:36:18 +0000447#endif // ALLOW_TRANSLATION_MODELS
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700448
Rachel Barker738c4f42023-01-11 17:36:18 +0000449static bool is_degenerate_affine(double *p) {
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700450 return is_collinear3(p, p + 2, p + 4);
451}
452
Rachel Barker738c4f42023-01-11 17:36:18 +0000453static const RansacModelInfo ransac_model_info[TRANS_TYPES] = {
454 // IDENTITY
455 { NULL, NULL, NULL, 0 },
456// TRANSLATION
457#if ALLOW_TRANSLATION_MODELS
458 { is_degenerate_translation, find_translation, project_points_translation,
459 3 },
460#else
461 { NULL, NULL, NULL, 0 },
462#endif
463 // ROTZOOM
464 { is_degenerate_affine, find_rotzoom, project_points_affine, 3 },
465 // AFFINE
466 { is_degenerate_affine, find_affine, project_points_affine, 3 },
467};
Sarah Parker4dc0f1b2016-08-09 17:40:53 -0700468
Rachel Barker738c4f42023-01-11 17:36:18 +0000469// Returns true on success, false on error
470bool ransac(const Correspondence *matched_points, int npoints,
471 TransformationType type, MotionModel *motion_models,
472 int num_desired_motions) {
473#if ALLOW_TRANSLATION_MODELS
474 assert(type > IDENTITY && type < TRANS_TYPES);
475#else
476 assert(type > TRANSLATION && type < TRANS_TYPES);
477#endif // ALLOW_TRANSLATION_MODELS
Sarah Parkerf9a961c2016-09-06 11:25:04 -0700478
Rachel Barker738c4f42023-01-11 17:36:18 +0000479 return ransac_internal(matched_points, npoints, motion_models,
480 num_desired_motions, &ransac_model_info[type]);
Yaowu Xu0e4509c2019-05-02 15:59:17 -0700481}