Add ransac interfaces for floating point correspondences

DISflow flow vectors will be floating point and so all correspondences
passed to ransac will also be floating point. Ransac already uses
floating point internally, so the only change needed is to the function
interface.

Change-Id: If9f1c6fbeceaa31f625e6fd38a55848dfca20261
diff --git a/av1/encoder/global_motion.c b/av1/encoder/global_motion.c
index 3a66c38..e35a208 100644
--- a/av1/encoder/global_motion.c
+++ b/av1/encoder/global_motion.c
@@ -324,6 +324,16 @@
   return 0;
 }
 #else
+static INLINE RansacFuncDouble
+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;
+  }
+}
+
 double getCubicValue(double p[4], double x) {
   return p[1] + 0.5 * x *
                     (p[2] - p[0] +
diff --git a/av1/encoder/ransac.c b/av1/encoder/ransac.c
index 821ac41..6a8854c 100644
--- a/av1/encoder/ransac.c
+++ b/av1/encoder/ransac.c
@@ -546,6 +546,182 @@
   return ret_val;
 }
 
+static int ransac_double_prec(const double *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) {
+  int trial_count = 0;
+  int i = 0;
+  int ret_val = 0;
+
+  unsigned int seed = (unsigned int)npoints;
+
+  int indices[MAX_MINPTS] = { 0 };
+
+  double *points1, *points2;
+  double *corners1, *corners2;
+  double *image1_coord;
+
+  // Store information for the num_desired_motions best transformations found
+  // and the worst motion among them, as well as the motion currently under
+  // consideration.
+  RANSAC_MOTION *motions, *worst_kept_motion = NULL;
+  RANSAC_MOTION current_motion;
+
+  // Store the parameters and the indices of the inlier points for the motion
+  // currently under consideration.
+  double params_this_motion[MAX_PARAMDIM];
+
+  double *cnp1, *cnp2;
+
+  for (i = 0; i < num_desired_motions; ++i) {
+    num_inliers_by_motion[i] = 0;
+  }
+  if (npoints < minpts * MINPTS_MULTIPLIER || npoints == 0) {
+    return 1;
+  }
+
+  points1 = (double *)aom_malloc(sizeof(*points1) * npoints * 2);
+  points2 = (double *)aom_malloc(sizeof(*points2) * npoints * 2);
+  corners1 = (double *)aom_malloc(sizeof(*corners1) * npoints * 2);
+  corners2 = (double *)aom_malloc(sizeof(*corners2) * npoints * 2);
+  image1_coord = (double *)aom_malloc(sizeof(*image1_coord) * npoints * 2);
+
+  motions =
+      (RANSAC_MOTION *)aom_malloc(sizeof(RANSAC_MOTION) * num_desired_motions);
+  for (i = 0; i < num_desired_motions; ++i) {
+    motions[i].inlier_indices =
+        (int *)aom_malloc(sizeof(*motions->inlier_indices) * npoints);
+    clear_motion(motions + i, npoints);
+  }
+  current_motion.inlier_indices =
+      (int *)aom_malloc(sizeof(*current_motion.inlier_indices) * npoints);
+  clear_motion(&current_motion, npoints);
+
+  worst_kept_motion = motions;
+
+  if (!(points1 && points2 && corners1 && corners2 && image1_coord && motions &&
+        current_motion.inlier_indices)) {
+    ret_val = 1;
+    goto finish_ransac;
+  }
+
+  cnp1 = corners1;
+  cnp2 = corners2;
+  for (i = 0; i < npoints; ++i) {
+    *(cnp1++) = *(matched_points++);
+    *(cnp1++) = *(matched_points++);
+    *(cnp2++) = *(matched_points++);
+    *(cnp2++) = *(matched_points++);
+  }
+
+  while (MIN_TRIALS > trial_count) {
+    double sum_distance = 0.0;
+    double sum_distance_squared = 0.0;
+
+    clear_motion(&current_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(&current_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 + (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 =
@@ -583,3 +759,33 @@
                 params_by_motion, num_desired_motions, 3, is_degenerate_affine,
                 find_affine, project_points_double_affine);
 }
+
+int ransac_translation_double_prec(double *matched_points, int npoints,
+                                   int *num_inliers_by_motion,
+                                   double *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);
+}
+
+int ransac_rotzoom_double_prec(double *matched_points, int npoints,
+                               int *num_inliers_by_motion,
+                               double *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);
+}
+
+int ransac_affine_double_prec(double *matched_points, int npoints,
+                              int *num_inliers_by_motion,
+                              double *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);
+}
diff --git a/av1/encoder/ransac.h b/av1/encoder/ransac.h
index c429f2c..b754bac 100644
--- a/av1/encoder/ransac.h
+++ b/av1/encoder/ransac.h
@@ -22,6 +22,9 @@
 typedef int (*RansacFunc)(int *matched_points, int npoints,
                           int *num_inliers_by_motion, double *params_by_motion,
                           int num_motions);
+typedef int (*RansacFuncDouble)(double *matched_points, int npoints,
+                                int *num_inliers_by_motion,
+                                double *params_by_motion, int num_motions);
 
 /* Each of these functions fits a motion model from a set of
    corresponding points in 2 frames using RANSAC. */
@@ -32,4 +35,16 @@
 int ransac_translation(int *matched_points, int npoints,
                        int *num_inliers_by_motion, double *params_by_motion,
                        int num_motions);
+int ransac_translation_double_prec(double *matched_points, int npoints,
+                                   int *num_inliers_by_motion,
+                                   double *params_by_motion,
+                                   int num_desired_motions);
+int ransac_rotzoom_double_prec(double *matched_points, int npoints,
+                               int *num_inliers_by_motion,
+                               double *params_by_motion,
+                               int num_desired_motions);
+int ransac_affine_double_prec(double *matched_points, int npoints,
+                              int *num_inliers_by_motion,
+                              double *params_by_motion,
+                              int num_desired_motions);
 #endif  // AOM_AV1_ENCODER_RANSAC_H_