From 68569dee1416593955c1570d638b3d9250b33012 Mon Sep 17 00:00:00 2001 From: trav90 Date: Mon, 15 Oct 2018 21:45:30 -0500 Subject: Import aom library This is the reference implementation for the Alliance for Open Media's av1 video code. The commit used was 4d668d7feb1f8abd809d1bca0418570a7f142a36. --- third_party/aom/av1/encoder/ransac.c | 1210 ++++++++++++++++++++++++++++++++++ 1 file changed, 1210 insertions(+) create mode 100644 third_party/aom/av1/encoder/ransac.c (limited to 'third_party/aom/av1/encoder/ransac.c') diff --git a/third_party/aom/av1/encoder/ransac.c b/third_party/aom/av1/encoder/ransac.c new file mode 100644 index 000000000..5d5dd7572 --- /dev/null +++ b/third_party/aom/av1/encoder/ransac.c @@ -0,0 +1,1210 @@ +/* + * 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. + */ +#define _POSIX_C_SOURCE 200112L // rand_r() +#include +#include +#include +#include +#include +#include + +#include "av1/encoder/ransac.h" + +#define MAX_MINPTS 4 +#define MAX_DEGENERATE_ITER 10 +#define MINPTS_MULTIPLIER 5 + +#define INLIER_THRESHOLD 1.0 +#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, const int n, + const int stride_points, + const int stride_proj); + +static void project_points_double_translation(double *mat, double *points, + double *proj, const int n, + const int stride_points, + const 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, const int n, + const int stride_points, + const 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, const int n, + const int stride_points, + const 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_double_hortrapezoid(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[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[5] * y + mat[1]) * Z; + points += stride_points - 2; + proj += stride_proj - 2; + } +} + +static void project_points_double_vertrapezoid(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 + 1; + assert(fabs(Z_inv) > 0.000001); + Z = 1. / Z_inv; + *(proj++) = (mat[2] * x + mat[0]) * Z; + *(proj++) = (mat[4] * x + mat[5] * y + mat[1]) * Z; + points += stride_points - 2; + proj += stride_proj - 2; + } +} + +static void project_points_double_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; + } +} + +/////////////////////////////////////////////////////////////////////////////// +// svdcmp +// Adopted from Numerical Recipes in C + +static const double TINY_NEAR_ZERO = 1.0E-12; + +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 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 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; + double anorm, c, f, g, h, s, scale, x, y, z; + double *rv1 = (double *)aom_malloc(sizeof(*rv1) * (n + 1)); + g = scale = anorm = 0.0; + for (i = 0; i < n; i++) { + l = i + 1; + rv1[i] = scale * g; + g = s = scale = 0.0; + if (i < m) { + for (k = i; k < m; k++) scale += fabs(u[k][i]); + if (scale != 0.) { + for (k = i; k < m; k++) { + u[k][i] /= scale; + s += u[k][i] * u[k][i]; + } + f = u[i][i]; + g = -sign(sqrt(s), f); + h = f * g - s; + u[i][i] = f - g; + for (j = l; j < n; j++) { + for (s = 0.0, k = i; k < m; k++) s += u[k][i] * u[k][j]; + f = s / h; + for (k = i; k < m; k++) u[k][j] += f * u[k][i]; + } + for (k = i; k < m; k++) u[k][i] *= scale; + } + } + w[i] = scale * g; + g = s = scale = 0.0; + if (i < m && i != n - 1) { + for (k = l; k < n; k++) scale += fabs(u[i][k]); + if (scale != 0.) { + for (k = l; k < n; k++) { + u[i][k] /= scale; + s += u[i][k] * u[i][k]; + } + f = u[i][l]; + g = -sign(sqrt(s), f); + h = f * g - s; + u[i][l] = f - g; + for (k = l; k < n; k++) rv1[k] = u[i][k] / h; + for (j = l; j < m; j++) { + for (s = 0.0, k = l; k < n; k++) s += u[j][k] * u[i][k]; + for (k = l; k < n; k++) u[j][k] += s * rv1[k]; + } + for (k = l; k < n; k++) u[i][k] *= scale; + } + } + anorm = fmax(anorm, (fabs(w[i]) + fabs(rv1[i]))); + } + + for (i = n - 1; i >= 0; i--) { + if (i < n - 1) { + if (g != 0.) { + for (j = l; j < n; j++) v[j][i] = (u[i][j] / u[i][l]) / g; + for (j = l; j < n; j++) { + for (s = 0.0, k = l; k < n; k++) s += u[i][k] * v[k][j]; + for (k = l; k < n; k++) v[k][j] += s * v[k][i]; + } + } + for (j = l; j < n; j++) v[i][j] = v[j][i] = 0.0; + } + v[i][i] = 1.0; + g = rv1[i]; + l = i; + } + for (i = AOMMIN(m, n) - 1; i >= 0; i--) { + l = i + 1; + g = w[i]; + for (j = l; j < n; j++) u[i][j] = 0.0; + if (g != 0.) { + g = 1.0 / g; + for (j = l; j < n; j++) { + for (s = 0.0, k = l; k < m; k++) s += u[k][i] * u[k][j]; + f = (s / u[i][i]) * g; + for (k = i; k < m; k++) u[k][j] += f * u[k][i]; + } + for (j = i; j < m; j++) u[j][i] *= g; + } else { + for (j = i; j < m; j++) u[j][i] = 0.0; + } + ++u[i][i]; + } + for (k = n - 1; k >= 0; k--) { + for (its = 0; its < max_its; its++) { + flag = 1; + for (l = k; l >= 0; l--) { + nm = l - 1; + if ((double)(fabs(rv1[l]) + anorm) == anorm || nm < 0) { + flag = 0; + break; + } + if ((double)(fabs(w[nm]) + anorm) == anorm) break; + } + if (flag) { + c = 0.0; + s = 1.0; + for (i = l; i <= k; i++) { + f = s * rv1[i]; + rv1[i] = c * rv1[i]; + if ((double)(fabs(f) + anorm) == anorm) break; + g = w[i]; + h = pythag(f, g); + w[i] = h; + h = 1.0 / h; + c = g * h; + s = -f * h; + for (j = 0; j < m; j++) { + y = u[j][nm]; + z = u[j][i]; + u[j][nm] = y * c + z * s; + u[j][i] = z * c - y * s; + } + } + } + z = w[k]; + if (l == k) { + if (z < 0.0) { + w[k] = -z; + for (j = 0; j < n; j++) v[j][k] = -v[j][k]; + } + break; + } + if (its == max_its - 1) { + aom_free(rv1); + return 1; + } + assert(k > 0); + x = w[l]; + nm = k - 1; + y = w[nm]; + g = rv1[nm]; + h = rv1[k]; + f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y); + g = pythag(f, 1.0); + f = ((x - z) * (x + z) + h * ((y / (f + sign(g, f))) - h)) / x; + c = s = 1.0; + for (j = l; j <= nm; j++) { + i = j + 1; + g = rv1[i]; + y = w[i]; + h = s * g; + g = c * g; + z = pythag(f, h); + rv1[j] = z; + c = f / z; + s = h / z; + f = x * c + g * s; + g = g * c - x * s; + h = y * s; + y *= c; + for (jj = 0; jj < n; jj++) { + x = v[jj][j]; + z = v[jj][i]; + v[jj][j] = x * c + z * s; + v[jj][i] = z * c - x * s; + } + z = pythag(f, h); + w[j] = z; + if (z != 0.) { + z = 1.0 / z; + c = f * z; + s = h * z; + } + f = c * g + s * y; + x = c * y - s * g; + for (jj = 0; jj < m; jj++) { + y = u[jj][j]; + z = u[jj][i]; + u[jj][j] = y * c + z * s; + u[jj][i] = z * c - y * s; + } + } + rv1[l] = 0.0; + rv1[k] = f; + w[k] = x; + } + } + aom_free(rv1); + return 0; +} + +static 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)); + int problem, i; + + problem = !(nrU && nrV); + if (!problem) { + for (i = 0; i < M; i++) { + nrU[i] = &U[i * N]; + } + for (i = 0; i < N; i++) { + nrV[i] = &V[i * N]; + } + } else { + if (nrU) aom_free(nrU); + if (nrV) aom_free(nrV); + return 1; + } + + /* copy from given matx into nrU */ + for (i = 0; i < M; i++) { + memcpy(&(nrU[i][0]), matx + N * i, N * sizeof(*matx)); + } + + /* HERE IT IS: do SVD */ + if (svdcmp(nrU, M, N, W, nrV)) { + aom_free(nrU); + aom_free(nrV); + return 1; + } + + /* aom_free Numerical Recipes arrays */ + aom_free(nrU); + aom_free(nrV); + + return 0; +} + +int pseudo_inverse(double *inv, double *matx, const int M, const int N) { + double ans; + int i, j, k; + double *const U = (double *)aom_malloc(M * N * sizeof(*matx)); + double *const W = (double *)aom_malloc(N * sizeof(*matx)); + double *const V = (double *)aom_malloc(N * N * sizeof(*matx)); + + if (!(U && W && V)) { + return 1; + } + if (SVD(U, W, V, matx, M, N)) { + aom_free(U); + aom_free(W); + aom_free(V); + return 1; + } + for (i = 0; i < N; i++) { + if (fabs(W[i]) < TINY_NEAR_ZERO) { + aom_free(U); + aom_free(W); + aom_free(V); + return 1; + } + } + + for (i = 0; i < N; i++) { + for (j = 0; j < M; j++) { + ans = 0; + for (k = 0; k < N; k++) { + ans += V[k + N * i] * U[k + N * j] / W[k]; + } + inv[j + M * i] = ans; + } + } + aom_free(U); + aom_free(W); + aom_free(V); + return 0; +} + +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; + 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 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 * 9); + 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 (pseudo_inverse(temp, a, np2, 4)) { + aom_free(a); + return 1; + } + multiply_mat(temp, b, mat, 4, np2, 1); + denormalize_rotzoom_reorder(mat, T1, T2); + aom_free(a); + return 0; +} + +static int find_affine(int np, double *pts1, double *pts2, double *mat) { + const int np2 = np * 2; + double *a = (double *)aom_malloc(sizeof(*a) * np2 * 13); + 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 (pseudo_inverse(temp, a, np2, 6)) { + aom_free(a); + return 1; + } + multiply_mat(temp, b, mat, 6, np2, 1); + denormalize_affine_reorder(mat, T1, T2); + aom_free(a); + return 0; +} + +static int find_vertrapezoid(int np, double *pts1, double *pts2, double *mat) { + const int np3 = np * 3; + double *a = (double *)aom_malloc(sizeof(*a) * np3 * 14); + double *U = a + np3 * 7; + double S[7], V[7 * 7], 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 * 7 + 0] = a[i * 3 * 7 + 1] = 0; + a[i * 3 * 7 + 2] = -sx; + a[i * 3 * 7 + 3] = -sy; + a[i * 3 * 7 + 4] = -1; + a[i * 3 * 7 + 5] = dy * sx; + a[i * 3 * 7 + 6] = dy; + + a[(i * 3 + 1) * 7 + 0] = sx; + a[(i * 3 + 1) * 7 + 1] = 1; + a[(i * 3 + 1) * 7 + 2] = a[(i * 3 + 1) * 7 + 3] = a[(i * 3 + 1) * 7 + 4] = + 0; + a[(i * 3 + 1) * 7 + 5] = -dx * sx; + a[(i * 3 + 1) * 7 + 6] = -dx; + + a[(i * 3 + 2) * 7 + 0] = -dy * sx; + a[(i * 3 + 2) * 7 + 1] = -dy; + a[(i * 3 + 2) * 7 + 2] = dx * sx; + a[(i * 3 + 2) * 7 + 3] = dx * sy; + a[(i * 3 + 2) * 7 + 4] = dx; + a[(i * 3 + 2) * 7 + 5] = a[(i * 3 + 2) * 7 + 6] = 0; + } + if (SVD(U, S, V, a, np3, 7)) { + aom_free(a); + return 1; + } else { + double minS = 1e12; + mini = -1; + for (i = 0; i < 7; ++i) { + if (S[i] < minS) { + minS = S[i]; + mini = i; + } + } + } + H[1] = H[7] = 0; + for (i = 0; i < 1; i++) H[i] = V[i * 7 + mini]; + for (; i < 6; i++) H[i + 1] = V[i * 7 + mini]; + for (; i < 7; i++) H[i + 2] = V[i * 7 + mini]; + + denormalize_homography_reorder(H, T1, T2); + aom_free(a); + if (H[8] == 0.0) { + return 1; + } else { + // normalize + double f = 1.0 / H[8]; + for (i = 0; i < 8; i++) mat[i] = f * H[i]; + } + return 0; +} + +static int find_hortrapezoid(int np, double *pts1, double *pts2, double *mat) { + const int np3 = np * 3; + double *a = (double *)aom_malloc(sizeof(*a) * np3 * 14); + double *U = a + np3 * 7; + double S[7], V[7 * 7], 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 * 7 + 0] = a[i * 3 * 7 + 1] = a[i * 3 * 7 + 2] = 0; + a[i * 3 * 7 + 3] = -sy; + a[i * 3 * 7 + 4] = -1; + a[i * 3 * 7 + 5] = dy * sy; + a[i * 3 * 7 + 6] = dy; + + a[(i * 3 + 1) * 7 + 0] = sx; + a[(i * 3 + 1) * 7 + 1] = sy; + a[(i * 3 + 1) * 7 + 2] = 1; + a[(i * 3 + 1) * 7 + 3] = a[(i * 3 + 1) * 7 + 4] = 0; + a[(i * 3 + 1) * 7 + 5] = -dx * sy; + a[(i * 3 + 1) * 7 + 6] = -dx; + + a[(i * 3 + 2) * 7 + 0] = -dy * sx; + a[(i * 3 + 2) * 7 + 1] = -dy * sy; + a[(i * 3 + 2) * 7 + 2] = -dy; + a[(i * 3 + 2) * 7 + 3] = dx * sy; + a[(i * 3 + 2) * 7 + 4] = dx; + a[(i * 3 + 2) * 7 + 5] = a[(i * 3 + 2) * 7 + 6] = 0; + } + + if (SVD(U, S, V, a, np3, 7)) { + aom_free(a); + return 1; + } else { + double minS = 1e12; + mini = -1; + for (i = 0; i < 7; ++i) { + if (S[i] < minS) { + minS = S[i]; + mini = i; + } + } + } + H[3] = H[6] = 0; + for (i = 0; i < 3; i++) H[i] = V[i * 7 + mini]; + for (; i < 5; i++) H[i + 1] = V[i * 7 + mini]; + for (; i < 7; i++) H[i + 2] = V[i * 7 + mini]; + + denormalize_homography_reorder(H, T1, T2); + aom_free(a); + if (H[8] == 0.0) { + return 1; + } else { + // normalize + double f = 1.0 / H[8]; + for (i = 0; i < 8; i++) mat[i] = f * H[i]; + } + return 0; +} + +static int 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 1; + } 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 1; + } else { + // normalize + double f = 1.0 / H[8]; + for (i = 0; i < 8; i++) mat[i] = f * H[i]; + } + return 0; +} + +static int get_rand_indices(int npoints, int minpts, int *indices, + unsigned int *seed) { + int i, j; + int ptr = rand_r(seed) % npoints; + if (minpts > npoints) return 0; + indices[0] = ptr; + ptr = (ptr == npoints - 1 ? 0 : ptr + 1); + i = 1; + while (i < minpts) { + int index = rand_r(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, double *params_by_motion, + int num_desired_motions, const int minpts, + IsDegenerateFunc is_degenerate, + FindTransformationFunc find_transformation, + ProjectPointsDoubleFunc projectpoints) { + static const double PROBABILITY_REQUIRED = 0.9; + static const double EPS = 1e-12; + + int N = 10000, 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; + + 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 (N > 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) { + int temp; + double fracinliers, pNoOutliers, 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); + fracinliers = (double)current_motion.num_inliers / (double)npoints; + pNoOutliers = 1 - pow(fracinliers, minpts); + pNoOutliers = fmax(EPS, pNoOutliers); + pNoOutliers = fmin(1 - EPS, pNoOutliers); + temp = (int)(log(1.0 - PROBABILITY_REQUIRED) / log(pNoOutliers)); + + if (temp > 0 && temp < N) { + N = AOMMAX(temp, MIN_TRIALS); + } + + // 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) { + 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 + (MAX_PARAMDIM - 1) * i); + 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 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); +} + +int ransac_translation(int *matched_points, int npoints, + int *num_inliers_by_motion, double *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); +} + +int ransac_rotzoom(int *matched_points, int npoints, int *num_inliers_by_motion, + double *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); +} + +int ransac_affine(int *matched_points, int npoints, int *num_inliers_by_motion, + double *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); +} + +int ransac_homography(int *matched_points, int npoints, + int *num_inliers_by_motion, double *params_by_motion, + int num_desired_motions) { + return ransac(matched_points, npoints, num_inliers_by_motion, + params_by_motion, num_desired_motions, 4, + is_degenerate_homography, find_homography, + project_points_double_homography); +} + +int ransac_hortrapezoid(int *matched_points, int npoints, + int *num_inliers_by_motion, double *params_by_motion, + int num_desired_motions) { + return ransac(matched_points, npoints, num_inliers_by_motion, + params_by_motion, num_desired_motions, 4, + is_degenerate_homography, find_hortrapezoid, + project_points_double_hortrapezoid); +} + +int ransac_vertrapezoid(int *matched_points, int npoints, + int *num_inliers_by_motion, double *params_by_motion, + int num_desired_motions) { + return ransac(matched_points, npoints, num_inliers_by_motion, + params_by_motion, num_desired_motions, 4, + is_degenerate_homography, find_vertrapezoid, + project_points_double_vertrapezoid); +} -- cgit v1.2.3 From df9477dfa60ebb5d31bc142e58ce46535c17abce Mon Sep 17 00:00:00 2001 From: trav90 Date: Wed, 17 Oct 2018 05:59:08 -0500 Subject: Update aom to slightly newer commit ID --- third_party/aom/av1/encoder/ransac.c | 335 +++-------------------------------- 1 file changed, 21 insertions(+), 314 deletions(-) (limited to 'third_party/aom/av1/encoder/ransac.c') diff --git a/third_party/aom/av1/encoder/ransac.c b/third_party/aom/av1/encoder/ransac.c index 5d5dd7572..bbd2d179c 100644 --- a/third_party/aom/av1/encoder/ransac.c +++ b/third_party/aom/av1/encoder/ransac.c @@ -8,7 +8,6 @@ * 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. */ -#define _POSIX_C_SOURCE 200112L // rand_r() #include #include #include @@ -17,6 +16,7 @@ #include #include "av1/encoder/ransac.h" +#include "av1/encoder/mathutils.h" #define MAX_MINPTS 4 #define MAX_DEGENERATE_ITER 10 @@ -133,309 +133,6 @@ static void project_points_double_homography(double *mat, double *points, } } -/////////////////////////////////////////////////////////////////////////////// -// svdcmp -// Adopted from Numerical Recipes in C - -static const double TINY_NEAR_ZERO = 1.0E-12; - -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 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 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; - double anorm, c, f, g, h, s, scale, x, y, z; - double *rv1 = (double *)aom_malloc(sizeof(*rv1) * (n + 1)); - g = scale = anorm = 0.0; - for (i = 0; i < n; i++) { - l = i + 1; - rv1[i] = scale * g; - g = s = scale = 0.0; - if (i < m) { - for (k = i; k < m; k++) scale += fabs(u[k][i]); - if (scale != 0.) { - for (k = i; k < m; k++) { - u[k][i] /= scale; - s += u[k][i] * u[k][i]; - } - f = u[i][i]; - g = -sign(sqrt(s), f); - h = f * g - s; - u[i][i] = f - g; - for (j = l; j < n; j++) { - for (s = 0.0, k = i; k < m; k++) s += u[k][i] * u[k][j]; - f = s / h; - for (k = i; k < m; k++) u[k][j] += f * u[k][i]; - } - for (k = i; k < m; k++) u[k][i] *= scale; - } - } - w[i] = scale * g; - g = s = scale = 0.0; - if (i < m && i != n - 1) { - for (k = l; k < n; k++) scale += fabs(u[i][k]); - if (scale != 0.) { - for (k = l; k < n; k++) { - u[i][k] /= scale; - s += u[i][k] * u[i][k]; - } - f = u[i][l]; - g = -sign(sqrt(s), f); - h = f * g - s; - u[i][l] = f - g; - for (k = l; k < n; k++) rv1[k] = u[i][k] / h; - for (j = l; j < m; j++) { - for (s = 0.0, k = l; k < n; k++) s += u[j][k] * u[i][k]; - for (k = l; k < n; k++) u[j][k] += s * rv1[k]; - } - for (k = l; k < n; k++) u[i][k] *= scale; - } - } - anorm = fmax(anorm, (fabs(w[i]) + fabs(rv1[i]))); - } - - for (i = n - 1; i >= 0; i--) { - if (i < n - 1) { - if (g != 0.) { - for (j = l; j < n; j++) v[j][i] = (u[i][j] / u[i][l]) / g; - for (j = l; j < n; j++) { - for (s = 0.0, k = l; k < n; k++) s += u[i][k] * v[k][j]; - for (k = l; k < n; k++) v[k][j] += s * v[k][i]; - } - } - for (j = l; j < n; j++) v[i][j] = v[j][i] = 0.0; - } - v[i][i] = 1.0; - g = rv1[i]; - l = i; - } - for (i = AOMMIN(m, n) - 1; i >= 0; i--) { - l = i + 1; - g = w[i]; - for (j = l; j < n; j++) u[i][j] = 0.0; - if (g != 0.) { - g = 1.0 / g; - for (j = l; j < n; j++) { - for (s = 0.0, k = l; k < m; k++) s += u[k][i] * u[k][j]; - f = (s / u[i][i]) * g; - for (k = i; k < m; k++) u[k][j] += f * u[k][i]; - } - for (j = i; j < m; j++) u[j][i] *= g; - } else { - for (j = i; j < m; j++) u[j][i] = 0.0; - } - ++u[i][i]; - } - for (k = n - 1; k >= 0; k--) { - for (its = 0; its < max_its; its++) { - flag = 1; - for (l = k; l >= 0; l--) { - nm = l - 1; - if ((double)(fabs(rv1[l]) + anorm) == anorm || nm < 0) { - flag = 0; - break; - } - if ((double)(fabs(w[nm]) + anorm) == anorm) break; - } - if (flag) { - c = 0.0; - s = 1.0; - for (i = l; i <= k; i++) { - f = s * rv1[i]; - rv1[i] = c * rv1[i]; - if ((double)(fabs(f) + anorm) == anorm) break; - g = w[i]; - h = pythag(f, g); - w[i] = h; - h = 1.0 / h; - c = g * h; - s = -f * h; - for (j = 0; j < m; j++) { - y = u[j][nm]; - z = u[j][i]; - u[j][nm] = y * c + z * s; - u[j][i] = z * c - y * s; - } - } - } - z = w[k]; - if (l == k) { - if (z < 0.0) { - w[k] = -z; - for (j = 0; j < n; j++) v[j][k] = -v[j][k]; - } - break; - } - if (its == max_its - 1) { - aom_free(rv1); - return 1; - } - assert(k > 0); - x = w[l]; - nm = k - 1; - y = w[nm]; - g = rv1[nm]; - h = rv1[k]; - f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y); - g = pythag(f, 1.0); - f = ((x - z) * (x + z) + h * ((y / (f + sign(g, f))) - h)) / x; - c = s = 1.0; - for (j = l; j <= nm; j++) { - i = j + 1; - g = rv1[i]; - y = w[i]; - h = s * g; - g = c * g; - z = pythag(f, h); - rv1[j] = z; - c = f / z; - s = h / z; - f = x * c + g * s; - g = g * c - x * s; - h = y * s; - y *= c; - for (jj = 0; jj < n; jj++) { - x = v[jj][j]; - z = v[jj][i]; - v[jj][j] = x * c + z * s; - v[jj][i] = z * c - x * s; - } - z = pythag(f, h); - w[j] = z; - if (z != 0.) { - z = 1.0 / z; - c = f * z; - s = h * z; - } - f = c * g + s * y; - x = c * y - s * g; - for (jj = 0; jj < m; jj++) { - y = u[jj][j]; - z = u[jj][i]; - u[jj][j] = y * c + z * s; - u[jj][i] = z * c - y * s; - } - } - rv1[l] = 0.0; - rv1[k] = f; - w[k] = x; - } - } - aom_free(rv1); - return 0; -} - -static 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)); - int problem, i; - - problem = !(nrU && nrV); - if (!problem) { - for (i = 0; i < M; i++) { - nrU[i] = &U[i * N]; - } - for (i = 0; i < N; i++) { - nrV[i] = &V[i * N]; - } - } else { - if (nrU) aom_free(nrU); - if (nrV) aom_free(nrV); - return 1; - } - - /* copy from given matx into nrU */ - for (i = 0; i < M; i++) { - memcpy(&(nrU[i][0]), matx + N * i, N * sizeof(*matx)); - } - - /* HERE IT IS: do SVD */ - if (svdcmp(nrU, M, N, W, nrV)) { - aom_free(nrU); - aom_free(nrV); - return 1; - } - - /* aom_free Numerical Recipes arrays */ - aom_free(nrU); - aom_free(nrV); - - return 0; -} - -int pseudo_inverse(double *inv, double *matx, const int M, const int N) { - double ans; - int i, j, k; - double *const U = (double *)aom_malloc(M * N * sizeof(*matx)); - double *const W = (double *)aom_malloc(N * sizeof(*matx)); - double *const V = (double *)aom_malloc(N * N * sizeof(*matx)); - - if (!(U && W && V)) { - return 1; - } - if (SVD(U, W, V, matx, M, N)) { - aom_free(U); - aom_free(W); - aom_free(V); - return 1; - } - for (i = 0; i < N; i++) { - if (fabs(W[i]) < TINY_NEAR_ZERO) { - aom_free(U); - aom_free(W); - aom_free(V); - return 1; - } - } - - for (i = 0; i < N; i++) { - for (j = 0; j < M; j++) { - ans = 0; - for (k = 0; k < N; k++) { - ans += V[k + N * i] * U[k + N * j] / W[k]; - } - inv[j + M * i] = ans; - } - } - aom_free(U); - aom_free(W); - aom_free(V); - return 0; -} - static void normalize_homography(double *pts, int n, double *T) { double *p = pts; double mean[2] = { 0, 0 }; @@ -597,7 +294,7 @@ static int find_translation(int np, double *pts1, double *pts2, double *mat) { 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 * 9); + double *a = (double *)aom_malloc(sizeof(*a) * (np2 * 5 + 20)); double *b = a + np2 * 4; double *temp = b + np2; int i; @@ -625,11 +322,10 @@ static int find_rotzoom(int np, double *pts1, double *pts2, double *mat) { b[2 * i] = dx; b[2 * i + 1] = dy; } - if (pseudo_inverse(temp, a, np2, 4)) { + if (!least_squares(4, a, np2, 4, b, temp, mat)) { aom_free(a); return 1; } - multiply_mat(temp, b, mat, 4, np2, 1); denormalize_rotzoom_reorder(mat, T1, T2); aom_free(a); return 0; @@ -637,7 +333,7 @@ static int find_rotzoom(int np, double *pts1, double *pts2, double *mat) { static int find_affine(int np, double *pts1, double *pts2, double *mat) { const int np2 = np * 2; - double *a = (double *)aom_malloc(sizeof(*a) * np2 * 13); + double *a = (double *)aom_malloc(sizeof(*a) * (np2 * 7 + 42)); double *b = a + np2 * 6; double *temp = b + np2; int i; @@ -669,11 +365,10 @@ static int find_affine(int np, double *pts1, double *pts2, double *mat) { b[2 * i] = dx; b[2 * i + 1] = dy; } - if (pseudo_inverse(temp, a, np2, 6)) { + if (!least_squares(6, a, np2, 6, b, temp, mat)) { aom_free(a); return 1; } - multiply_mat(temp, b, mat, 6, np2, 1); denormalize_affine_reorder(mat, T1, T2); aom_free(a); return 0; @@ -890,16 +585,22 @@ static int find_homography(int np, double *pts1, double *pts2, double *mat) { return 0; } +// Generate a random number in the range [0, 32768). +static unsigned int lcg_rand16(unsigned int *state) { + *state = (unsigned int)(*state * 1103515245ULL + 12345); + return *state / 65536 % 32768; +} + static int get_rand_indices(int npoints, int minpts, int *indices, unsigned int *seed) { int i, j; - int ptr = rand_r(seed) % npoints; + 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 = rand_r(seed) % npoints; + int index = lcg_rand16(seed) % npoints; while (index) { ptr = (ptr == npoints - 1 ? 0 : ptr + 1); for (j = 0; j < i; ++j) { @@ -986,6 +687,9 @@ static int ransac(const int *matched_points, int npoints, 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; } @@ -1072,7 +776,7 @@ static int ransac(const int *matched_points, int npoints, if (current_motion.num_inliers >= worst_kept_motion->num_inliers && current_motion.num_inliers > 1) { int temp; - double fracinliers, pNoOutliers, mean_distance; + double fracinliers, pNoOutliers, mean_distance, dtemp; mean_distance = sum_distance / ((double)current_motion.num_inliers); current_motion.variance = sum_distance_squared / ((double)current_motion.num_inliers - 1.0) - @@ -1092,7 +796,10 @@ static int ransac(const int *matched_points, int npoints, pNoOutliers = 1 - pow(fracinliers, minpts); pNoOutliers = fmax(EPS, pNoOutliers); pNoOutliers = fmin(1 - EPS, pNoOutliers); - temp = (int)(log(1.0 - PROBABILITY_REQUIRED) / log(pNoOutliers)); + dtemp = log(1.0 - PROBABILITY_REQUIRED) / log(pNoOutliers); + temp = (dtemp > (double)INT32_MAX) + ? INT32_MAX + : dtemp < (double)INT32_MIN ? INT32_MIN : (int)dtemp; if (temp > 0 && temp < N) { N = AOMMAX(temp, MIN_TRIALS); -- cgit v1.2.3 From 7369c7d7a5eed32963d8af37658286617919f91c Mon Sep 17 00:00:00 2001 From: trav90 Date: Thu, 18 Oct 2018 06:04:57 -0500 Subject: Update aom to commit id f5bdeac22930ff4c6b219be49c843db35970b918 --- third_party/aom/av1/encoder/ransac.c | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) (limited to 'third_party/aom/av1/encoder/ransac.c') diff --git a/third_party/aom/av1/encoder/ransac.c b/third_party/aom/av1/encoder/ransac.c index bbd2d179c..c6e3675be 100644 --- a/third_party/aom/av1/encoder/ransac.c +++ b/third_party/aom/av1/encoder/ransac.c @@ -139,6 +139,8 @@ static void normalize_homography(double *pts, int n, double *T) { 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]; @@ -821,13 +823,15 @@ static int ransac(const int *matched_points, int npoints, // Recompute the motions using only the inliers. for (i = 0; i < num_desired_motions; ++i) { - 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 + (MAX_PARAMDIM - 1) * 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 + (MAX_PARAMDIM - 1) * i); + } num_inliers_by_motion[i] = motions[i].num_inliers; } -- cgit v1.2.3 From ec910d81405c736a4490383a250299a7837c2e64 Mon Sep 17 00:00:00 2001 From: trav90 Date: Thu, 18 Oct 2018 21:53:44 -0500 Subject: Update aom to commit id e87fb2378f01103d5d6e477a4ef6892dc714e614 --- third_party/aom/av1/encoder/ransac.c | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) (limited to 'third_party/aom/av1/encoder/ransac.c') diff --git a/third_party/aom/av1/encoder/ransac.c b/third_party/aom/av1/encoder/ransac.c index c6e3675be..6d2eb4183 100644 --- a/third_party/aom/av1/encoder/ransac.c +++ b/third_party/aom/av1/encoder/ransac.c @@ -17,6 +17,7 @@ #include "av1/encoder/ransac.h" #include "av1/encoder/mathutils.h" +#include "av1/encoder/random.h" #define MAX_MINPTS 4 #define MAX_DEGENERATE_ITER 10 @@ -587,12 +588,6 @@ static int find_homography(int np, double *pts1, double *pts2, double *mat) { return 0; } -// Generate a random number in the range [0, 32768). -static unsigned int lcg_rand16(unsigned int *state) { - *state = (unsigned int)(*state * 1103515245ULL + 12345); - return *state / 65536 % 32768; -} - static int get_rand_indices(int npoints, int minpts, int *indices, unsigned int *seed) { int i, j; -- cgit v1.2.3 From bbcc64772580c8a979288791afa02d30bc476d2e Mon Sep 17 00:00:00 2001 From: trav90 Date: Fri, 19 Oct 2018 21:52:15 -0500 Subject: Update aom to v1.0.0 Update aom to commit id d14c5bb4f336ef1842046089849dee4a301fbbf0. --- third_party/aom/av1/encoder/ransac.c | 313 ----------------------------------- 1 file changed, 313 deletions(-) (limited to 'third_party/aom/av1/encoder/ransac.c') diff --git a/third_party/aom/av1/encoder/ransac.c b/third_party/aom/av1/encoder/ransac.c index 6d2eb4183..781f528eb 100644 --- a/third_party/aom/av1/encoder/ransac.c +++ b/third_party/aom/av1/encoder/ransac.c @@ -80,60 +80,6 @@ static void project_points_double_affine(double *mat, double *points, } } -static void project_points_double_hortrapezoid(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[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[5] * y + mat[1]) * Z; - points += stride_points - 2; - proj += stride_proj - 2; - } -} - -static void project_points_double_vertrapezoid(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 + 1; - assert(fabs(Z_inv) > 0.000001); - Z = 1. / Z_inv; - *(proj++) = (mat[2] * x + mat[0]) * Z; - *(proj++) = (mat[4] * x + mat[5] * y + mat[1]) * Z; - points += stride_points - 2; - proj += stride_proj - 2; - } -} - -static void project_points_double_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 }; @@ -193,22 +139,6 @@ static void denormalize_homography(double *params, double *T1, double *T2) { 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]; @@ -377,217 +307,6 @@ static int find_affine(int np, double *pts1, double *pts2, double *mat) { return 0; } -static int find_vertrapezoid(int np, double *pts1, double *pts2, double *mat) { - const int np3 = np * 3; - double *a = (double *)aom_malloc(sizeof(*a) * np3 * 14); - double *U = a + np3 * 7; - double S[7], V[7 * 7], 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 * 7 + 0] = a[i * 3 * 7 + 1] = 0; - a[i * 3 * 7 + 2] = -sx; - a[i * 3 * 7 + 3] = -sy; - a[i * 3 * 7 + 4] = -1; - a[i * 3 * 7 + 5] = dy * sx; - a[i * 3 * 7 + 6] = dy; - - a[(i * 3 + 1) * 7 + 0] = sx; - a[(i * 3 + 1) * 7 + 1] = 1; - a[(i * 3 + 1) * 7 + 2] = a[(i * 3 + 1) * 7 + 3] = a[(i * 3 + 1) * 7 + 4] = - 0; - a[(i * 3 + 1) * 7 + 5] = -dx * sx; - a[(i * 3 + 1) * 7 + 6] = -dx; - - a[(i * 3 + 2) * 7 + 0] = -dy * sx; - a[(i * 3 + 2) * 7 + 1] = -dy; - a[(i * 3 + 2) * 7 + 2] = dx * sx; - a[(i * 3 + 2) * 7 + 3] = dx * sy; - a[(i * 3 + 2) * 7 + 4] = dx; - a[(i * 3 + 2) * 7 + 5] = a[(i * 3 + 2) * 7 + 6] = 0; - } - if (SVD(U, S, V, a, np3, 7)) { - aom_free(a); - return 1; - } else { - double minS = 1e12; - mini = -1; - for (i = 0; i < 7; ++i) { - if (S[i] < minS) { - minS = S[i]; - mini = i; - } - } - } - H[1] = H[7] = 0; - for (i = 0; i < 1; i++) H[i] = V[i * 7 + mini]; - for (; i < 6; i++) H[i + 1] = V[i * 7 + mini]; - for (; i < 7; i++) H[i + 2] = V[i * 7 + mini]; - - denormalize_homography_reorder(H, T1, T2); - aom_free(a); - if (H[8] == 0.0) { - return 1; - } else { - // normalize - double f = 1.0 / H[8]; - for (i = 0; i < 8; i++) mat[i] = f * H[i]; - } - return 0; -} - -static int find_hortrapezoid(int np, double *pts1, double *pts2, double *mat) { - const int np3 = np * 3; - double *a = (double *)aom_malloc(sizeof(*a) * np3 * 14); - double *U = a + np3 * 7; - double S[7], V[7 * 7], 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 * 7 + 0] = a[i * 3 * 7 + 1] = a[i * 3 * 7 + 2] = 0; - a[i * 3 * 7 + 3] = -sy; - a[i * 3 * 7 + 4] = -1; - a[i * 3 * 7 + 5] = dy * sy; - a[i * 3 * 7 + 6] = dy; - - a[(i * 3 + 1) * 7 + 0] = sx; - a[(i * 3 + 1) * 7 + 1] = sy; - a[(i * 3 + 1) * 7 + 2] = 1; - a[(i * 3 + 1) * 7 + 3] = a[(i * 3 + 1) * 7 + 4] = 0; - a[(i * 3 + 1) * 7 + 5] = -dx * sy; - a[(i * 3 + 1) * 7 + 6] = -dx; - - a[(i * 3 + 2) * 7 + 0] = -dy * sx; - a[(i * 3 + 2) * 7 + 1] = -dy * sy; - a[(i * 3 + 2) * 7 + 2] = -dy; - a[(i * 3 + 2) * 7 + 3] = dx * sy; - a[(i * 3 + 2) * 7 + 4] = dx; - a[(i * 3 + 2) * 7 + 5] = a[(i * 3 + 2) * 7 + 6] = 0; - } - - if (SVD(U, S, V, a, np3, 7)) { - aom_free(a); - return 1; - } else { - double minS = 1e12; - mini = -1; - for (i = 0; i < 7; ++i) { - if (S[i] < minS) { - minS = S[i]; - mini = i; - } - } - } - H[3] = H[6] = 0; - for (i = 0; i < 3; i++) H[i] = V[i * 7 + mini]; - for (; i < 5; i++) H[i + 1] = V[i * 7 + mini]; - for (; i < 7; i++) H[i + 2] = V[i * 7 + mini]; - - denormalize_homography_reorder(H, T1, T2); - aom_free(a); - if (H[8] == 0.0) { - return 1; - } else { - // normalize - double f = 1.0 / H[8]; - for (i = 0; i < 8; i++) mat[i] = f * H[i]; - } - return 0; -} - -static int 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 1; - } 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 1; - } else { - // normalize - double f = 1.0 / H[8]; - for (i = 0; i < 8; i++) mat[i] = f * H[i]; - } - return 0; -} - static int get_rand_indices(int npoints, int minpts, int *indices, unsigned int *seed) { int i, j; @@ -860,11 +579,6 @@ static int is_degenerate_affine(double *p) { return is_collinear3(p, p + 2, p + 4); } -static int 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); -} - int ransac_translation(int *matched_points, int npoints, int *num_inliers_by_motion, double *params_by_motion, int num_desired_motions) { @@ -887,30 +601,3 @@ int ransac_affine(int *matched_points, int npoints, int *num_inliers_by_motion, params_by_motion, num_desired_motions, 3, is_degenerate_affine, find_affine, project_points_double_affine); } - -int ransac_homography(int *matched_points, int npoints, - int *num_inliers_by_motion, double *params_by_motion, - int num_desired_motions) { - return ransac(matched_points, npoints, num_inliers_by_motion, - params_by_motion, num_desired_motions, 4, - is_degenerate_homography, find_homography, - project_points_double_homography); -} - -int ransac_hortrapezoid(int *matched_points, int npoints, - int *num_inliers_by_motion, double *params_by_motion, - int num_desired_motions) { - return ransac(matched_points, npoints, num_inliers_by_motion, - params_by_motion, num_desired_motions, 4, - is_degenerate_homography, find_hortrapezoid, - project_points_double_hortrapezoid); -} - -int ransac_vertrapezoid(int *matched_points, int npoints, - int *num_inliers_by_motion, double *params_by_motion, - int num_desired_motions) { - return ransac(matched_points, npoints, num_inliers_by_motion, - params_by_motion, num_desired_motions, 4, - is_degenerate_homography, find_vertrapezoid, - project_points_double_vertrapezoid); -} -- cgit v1.2.3