Implement global motion parameter computation

This computes global motion parameters between 2 frames by
matching corresponding points using FAST feature and then
fitting a model using RANSAC.

Change-Id: Ib6664df44090e8cfa4db9f2f9e0556931ccfe5c8

third_party/fastfeat/LICENSE

+Copyright (c) 2006, 2008 Edward Rosten
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+
+	*Redistributions of source code must retain the above copyright
+	 notice, this list of conditions and the following disclaimer.
+
+	*Redistributions in binary form must reproduce the above copyright
+	 notice, this list of conditions and the following disclaimer in the
+	 documentation and/or other materials provided with the distribution.
+
+	*Neither the name of the University of Cambridge nor the names of 
+	 its contributors may be used to endorse or promote products derived 
+	 from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

third_party/fastfeat/README

+This code was taken from http://www.edwardrosten.com/work/fast.html. Fast 10, 11, and 12
+have been deleted.
+
+FAST feature detectors in C Version 2.0
+---------------------------------------
+
+The files are valid C and C++ code, and have no special requirements for
+compiling, and they do not depend on any libraries. Just compile them along with
+the rest of your project.
+
+To use the functions, #include "fast.h"
+
+The corner detectors have the following prototype (where X is 9, 10, 11 or 12):
+
+xy* fastX_detect_nonmax(const unsigned char * data, int xsize, int ysize, int stride, int threshold, int* numcorners)
+
+Where xy is the following simple struct typedef:
+
+typedef struct
+{
+	int x, y;
+} xy;
+
+The image is passed in as a block of data and dimensions, and the list of
+corners is returned as an array of xy structs, and an integer (numcorners)
+with the number of corners returned.  The data can be deallocated with free().
+Nonmaximal suppression is performed on the corners. Note that the stride
+is the number of bytes between rows. If your image has no padding, then this
+is the same as xsize.
+
+The detection, scoring and nonmaximal suppression are available as individual
+functions.  To see how to use the individual functions, see fast.c
+
+

third_party/fastfeat/fast.c

+// clang-format off
+#include <stdlib.h>
+#include "fast.h"
+
+
+xy* fast9_detect_nonmax(const byte* im, int xsize, int ysize, int stride, int b, int* ret_num_corners)
+{
+	xy* corners;
+	int num_corners;
+	int* scores;
+	xy* nonmax;
+
+	corners = fast9_detect(im, xsize, ysize, stride, b, &num_corners);
+	scores = fast9_score(im, stride, corners, num_corners, b);
+	nonmax = nonmax_suppression(corners, scores, num_corners, ret_num_corners);
+
+	free(corners);
+	free(scores);
+
+	return nonmax;
+}
+// clang-format on

third_party/fastfeat/fast.h

+// clang-format off
+#ifndef FAST_H
+#define FAST_H
+
+typedef struct { int x, y; } xy;
+typedef unsigned char byte;
+
+int fast9_corner_score(const byte* p, const int pixel[], int bstart);
+
+xy* fast9_detect(const byte* im, int xsize, int ysize, int stride, int b, int* ret_num_corners);
+
+int* fast9_score(const byte* i, int stride, xy* corners, int num_corners, int b);
+
+xy* fast9_detect_nonmax(const byte* im, int xsize, int ysize, int stride, int b, int* ret_num_corners);
+
+xy* nonmax_suppression(const xy* corners, const int* scores, int num_corners, int* ret_num_nonmax);
+
+
+#endif
+// clang-format on

third_party/fastfeat/nonmax.c

+// clang-format off
+#include <stdlib.h>
+#include "fast.h"
+
+
+#define Compare(X, Y) ((X)>=(Y))
+
+xy* nonmax_suppression(const xy* corners, const int* scores, int num_corners, int* ret_num_nonmax)
+{
+	int num_nonmax=0;
+	int last_row;
+	int* row_start;
+	int i, j;
+	xy* ret_nonmax;
+	const int sz = (int)num_corners;
+
+	/*Point above points (roughly) to the pixel above the one of interest, if there
+    is a feature there.*/
+	int point_above = 0;
+	int point_below = 0;
+
+
+	if(num_corners < 1)
+	{
+		*ret_num_nonmax = 0;
+		return 0;
+	}
+
+	ret_nonmax = (xy*)malloc(num_corners * sizeof(xy));
+
+	/* Find where each row begins
+	   (the corners are output in raster scan order). A beginning of -1 signifies
+	   that there are no corners on that row. */
+	last_row = corners[num_corners-1].y;
+	row_start = (int*)malloc((last_row+1)*sizeof(int));
+
+	for(i=0; i < last_row+1; i++)
+		row_start[i] = -1;
+
+	{
+		int prev_row = -1;
+		for(i=0; i< num_corners; i++)
+			if(corners[i].y != prev_row)
+			{
+				row_start[corners[i].y] = i;
+				prev_row = corners[i].y;
+			}
+	}
+
+
+
+	for(i=0; i < sz; i++)
+	{
+		int score = scores[i];
+		xy pos = corners[i];
+
+		/*Check left */
+		if(i > 0)
+			if(corners[i-1].x == pos.x-1 && corners[i-1].y == pos.y && Compare(scores[i-1], score))
+				continue;
+
+		/*Check right*/
+		if(i < (sz - 1))
+			if(corners[i+1].x == pos.x+1 && corners[i+1].y == pos.y && Compare(scores[i+1], score))
+				continue;
+
+		/*Check above (if there is a valid row above)*/
+		if(pos.y != 0 && row_start[pos.y - 1] != -1)
+		{
+			/*Make sure that current point_above is one
+			  row above.*/
+			if(corners[point_above].y < pos.y - 1)
+				point_above = row_start[pos.y-1];
+
+			/*Make point_above point to the first of the pixels above the current point,
+			  if it exists.*/
+			for(; corners[point_above].y < pos.y && corners[point_above].x < pos.x - 1; point_above++)
+			{}
+
+
+			for(j=point_above; corners[j].y < pos.y && corners[j].x <= pos.x + 1; j++)
+			{
+				int x = corners[j].x;
+				if( (x == pos.x - 1 || x ==pos.x || x == pos.x+1) && Compare(scores[j], score))
+					goto cont;
+			}
+
+		}
+
+		/*Check below (if there is anything below)*/
+		if(pos.y != last_row && row_start[pos.y + 1] != -1 && point_below < sz) /*Nothing below*/
+		{
+			if(corners[point_below].y < pos.y + 1)
+				point_below = row_start[pos.y+1];
+
+			/* Make point below point to one of the pixels belowthe current point, if it
+			   exists.*/
+			for(; point_below < sz && corners[point_below].y == pos.y+1 && corners[point_below].x < pos.x - 1; point_below++)
+			{}
+
+			for(j=point_below; j < sz && corners[j].y == pos.y+1 && corners[j].x <= pos.x + 1; j++)
+			{
+				int x = corners[j].x;
+				if( (x == pos.x - 1 || x ==pos.x || x == pos.x+1) && Compare(scores[j],score))
+					goto cont;
+			}
+		}
+
+		ret_nonmax[num_nonmax++] = corners[i];
+		cont:
+			;
+	}
+
+	free(row_start);
+	*ret_num_nonmax = num_nonmax;
+	return ret_nonmax;
+}
+
+// clang-format on

vp10/common/mv.h

@@ -52,7 +52,7 @@ typedef struct mv32 {
 // for each parameter. In other words, after a parameter is integerized
 // it is clamped between -(1 << ABS_XXX_BITS) and (1 << ABS_XXX_BITS).
 //
-// XXX_PREC_DIFF, XXX_ENCODE_FACTOR and XXX_DECODE_FACTOR
+// XXX_PREC_DIFF and XXX_DECODE_FACTOR
 // are computed once here to prevent repetitive
 // computation on the decoder side. These are
 // to allow the global motion parameters to be encoded in a lower
@@ -64,12 +64,10 @@ typedef struct mv32 {
 #define GM_TRANS_PREC_BITS 5
 #define GM_TRANS_PREC_DIFF (WARPEDMODEL_PREC_BITS - GM_TRANS_PREC_BITS)
 #define GM_TRANS_DECODE_FACTOR (1 << GM_TRANS_PREC_DIFF)
-#define GM_TRANS_ENCODE_FACTOR (1 / (GM_TRANS_DECODE_FACTOR))
 
 #define GM_ALPHA_PREC_BITS 5
 #define GM_ALPHA_PREC_DIFF (WARPEDMODEL_PREC_BITS - GM_ALPHA_PREC_BITS)
 #define GM_ALPHA_DECODE_FACTOR (1 << GM_ALPHA_PREC_DIFF)
-#define GM_ALPHA_ENCODE_FACTOR (1 / (GM_ALPHA_DECODE_FACTOR))
 
 #define GM_ABS_ALPHA_BITS 8
 #define GM_ABS_TRANS_BITS 8

vp10/common/warped_motion.h

@@ -37,6 +37,24 @@
 
 #define WARPEDDIFF_PREC_BITS (WARPEDMODEL_PREC_BITS - WARPEDPIXEL_PREC_BITS)
 
+typedef void (*ProjectPointsType)(int *mat, int *points, int *proj, const int n,
+                                  const int stride_points,
+                                  const int stride_proj,
+                                  const int subsampling_x,
+                                  const int subsampling_y);
+void projectPointsHomography(int *mat, int *points, int *proj, const int n,
+                             const int stride_points, const int stride_proj,
+                             const int subsampling_x, const int subsampling_y);
+void projectPointsAffine(int *mat, int *points, int *proj, const int n,
+                         const int stride_points, const int stride_proj,
+                         const int subsampling_x, const int subsampling_y);
+void projectPointsRotZoom(int *mat, int *points, int *proj, const int n,
+                          const int stride_points, const int stride_proj,
+                          const int subsampling_x, const int subsampling_y);
+void projectPointsTranslation(int *mat, int *points, int *proj, const int n,
+                              const int stride_points, const int stride_proj,
+                              const int subsampling_x, const int subsampling_y);
+
 typedef enum {
   UNKNOWN_TRANSFORM = -1,
   HOMOGRAPHY,   // homography, 8-parameter
@@ -54,6 +72,15 @@ typedef struct {
   int wmmat[8];  // For homography wmmat[9] is assumed to be 1
 } WarpedMotionParams;
 
+double vp10_warp_erroradv(WarpedMotionParams *wm,
+#if CONFIG_VP9_HIGHBITDEPTH
+                          int use_hbd, int bd,
+#endif  // CONFIG_VP9_HIGHBITDEPTH
+                          uint8_t *ref, int width, int height, int stride,
+                          uint8_t *dst, int p_col, int p_row, int p_width,
+                          int p_height, int p_stride, int subsampling_x,
+                          int subsampling_y, int x_scale, int y_scale);
+
 void vp10_warp_plane(WarpedMotionParams *wm,
 #if CONFIG_VP9_HIGHBITDEPTH
                      int use_hbd, int bd,

vp10/encoder/bitstream.c

@@ -3203,28 +3203,28 @@ static void write_global_motion_params(Global_Motion_Params *params,
     case GLOBAL_ZERO: break;
     case GLOBAL_AFFINE:
       vp10_write_primitive_symmetric(
-          w, params->motion_params.wmmat[4] * GM_ALPHA_ENCODE_FACTOR,
+          w, params->motion_params.wmmat[4] >> GM_ALPHA_PREC_DIFF,
           GM_ABS_ALPHA_BITS);
       vp10_write_primitive_symmetric(
-          w, (params->motion_params.wmmat[5] * GM_ALPHA_ENCODE_FACTOR) -
+          w, (params->motion_params.wmmat[5] >> GM_ALPHA_PREC_DIFF) -
                  (1 << GM_ALPHA_PREC_BITS),
           GM_ABS_ALPHA_BITS);
     // fallthrough intended
     case GLOBAL_ROTZOOM:
       vp10_write_primitive_symmetric(
-          w, (params->motion_params.wmmat[2] * GM_ALPHA_ENCODE_FACTOR) -
+          w, (params->motion_params.wmmat[2] >> GM_ALPHA_PREC_DIFF) -
                  (1 << GM_ALPHA_PREC_BITS),
           GM_ABS_ALPHA_BITS);
       vp10_write_primitive_symmetric(
-          w, params->motion_params.wmmat[3] * GM_ALPHA_ENCODE_FACTOR,
+          w, params->motion_params.wmmat[3] >> GM_ALPHA_PREC_DIFF,
           GM_ABS_ALPHA_BITS);
     // fallthrough intended
     case GLOBAL_TRANSLATION:
       vp10_write_primitive_symmetric(
-          w, params->motion_params.wmmat[0] * GM_TRANS_ENCODE_FACTOR,
+          w, params->motion_params.wmmat[0] >> GM_TRANS_PREC_DIFF,
           GM_ABS_TRANS_BITS);
       vp10_write_primitive_symmetric(
-          w, params->motion_params.wmmat[1] * GM_TRANS_ENCODE_FACTOR,
+          w, params->motion_params.wmmat[1] >> GM_TRANS_PREC_DIFF,
           GM_ABS_TRANS_BITS);
       break;
     default: assert(0);

vp10/encoder/corner_detect.c

+/*
+ *  Copyright (c) 2010 The WebM project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <memory.h>
+#include <math.h>
+#include <assert.h>
+
+#include "vp10/encoder/corner_detect.h"
+#include "third_party/fastfeat/fast.h"
+
+// Fast_9 wrapper
+#define FAST_BARRIER 40
+int FastCornerDetect(unsigned char *buf, int width, int height, int stride,
+                     int *points, int max_points) {
+  int num_points;
+  xy *frm_corners_xy = fast9_detect_nonmax(buf, width, height, stride,
+                                           FAST_BARRIER, &num_points);
+  num_points = (num_points <= max_points ? num_points : max_points);
+  if (num_points > 0 && frm_corners_xy) {
+    memcpy(points, frm_corners_xy, sizeof(xy) * num_points);
+    free(frm_corners_xy);
+    return num_points;
+  } else {
+    return 0;
+  }
+}

vp10/encoder/corner_detect.h

+/*
+ *  Copyright (c) 2010 The WebM project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#ifndef VP10_ENCODER_CORNER_DETECT_H_
+#define VP10_ENCODER_CORNER_DETECT_H_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <memory.h>
+
+int FastCornerDetect(unsigned char *buf, int width, int height, int stride,
+                     int *points, int max_points);
+
+#endif  // VP10_ENCODER_CORNER_DETECT_H

vp10/encoder/corner_match.c

+/*
+ *  Copyright (c) 2010 The WebM project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <memory.h>
+#include <math.h>
+
+#include "vp10/encoder/corner_match.h"
+
+#define MATCH_SZ 15
+#define MATCH_SZ_BY2 ((MATCH_SZ - 1) / 2)
+#define MATCH_SZ_SQ (MATCH_SZ * MATCH_SZ)
+#define SEARCH_SZ 9
+#define SEARCH_SZ_BY2 ((SEARCH_SZ - 1) / 2)
+
+#define THRESHOLD_NCC 0.80
+
+static double compute_variance(unsigned char *im, int stride, int x, int y,
+                               double *mean) {
+  double sum = 0.0;
+  double sumsq = 0.0;
+  double var;
+  int i, j;
+  for (i = 0; i < MATCH_SZ; ++i)
+    for (j = 0; j < MATCH_SZ; ++j) {
+      sum += im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)];
+      sumsq += im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)] *
+               im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)];
+    }
+  var = (sumsq * MATCH_SZ_SQ - sum * sum) / (MATCH_SZ_SQ * MATCH_SZ_SQ);
+  if (mean) *mean = sum / MATCH_SZ_SQ;
+  return var;
+}
+
+static double compute_cross_correlation(unsigned char *im1, int stride1, int x1,
+                                        int y1, unsigned char *im2, int stride2,
+                                        int x2, int y2) {
+  double sum1 = 0;
+  double sum2 = 0;
+  double cross = 0;
+  double corr;
+  int i, j;
+  for (i = 0; i < MATCH_SZ; ++i)
+    for (j = 0; j < MATCH_SZ; ++j) {
+      sum1 += im1[(i + y1 - MATCH_SZ_BY2) * stride1 + (j + x1 - MATCH_SZ_BY2)];
+      sum2 += im2[(i + y2 - MATCH_SZ_BY2) * stride2 + (j + x2 - MATCH_SZ_BY2)];
+      cross +=
+          im1[(i + y1 - MATCH_SZ_BY2) * stride1 + (j + x1 - MATCH_SZ_BY2)] *
+          im2[(i + y2 - MATCH_SZ_BY2) * stride2 + (j + x2 - MATCH_SZ_BY2)];
+    }
+  corr = (cross * MATCH_SZ_SQ - sum1 * sum2) / (MATCH_SZ_SQ * MATCH_SZ_SQ);
+  return corr;
+}
+
+static int is_eligible_point(double pointx, double pointy, int width,
+                             int height) {
+  return (pointx >= MATCH_SZ_BY2 && pointy >= MATCH_SZ_BY2 &&
+          pointx + MATCH_SZ_BY2 < width && pointy + MATCH_SZ_BY2 < height);
+}
+
+static int is_eligible_distance(double point1x, double point1y, double point2x,
+                                double point2y, int width, int height) {
+  const int thresh = (width < height ? height : width) >> 4;
+  return ((point1x - point2x) * (point1x - point2x) +
+          (point1y - point2y) * (point1y - point2y)) <= thresh * thresh;
+}
+
+static void improve_correspondence(unsigned char *frm, unsigned char *ref,
+                                   int width, int height, int frm_stride,
+                                   int ref_stride,
+                                   correspondence *correspondences,
+                                   int num_correspondences) {
+  int i;
+  for (i = 0; i < num_correspondences; ++i) {
+    double template_norm =
+        compute_variance(frm, frm_stride, (int)correspondences[i].x,
+                         (int)correspondences[i].y, NULL);
+    int x, y, best_x = 0, best_y = 0;
+    double best_match_ncc = 0.0;
+    for (y = -SEARCH_SZ_BY2; y <= SEARCH_SZ_BY2; ++y) {
+      for (x = -SEARCH_SZ_BY2; x <= SEARCH_SZ_BY2; ++x) {
+        double match_ncc;
+        double subimage_norm;
+        if (!is_eligible_point((int)correspondences[i].rx + x,
+                               (int)correspondences[i].ry + y, width, height))
+          continue;
+        if (!is_eligible_distance(
+                (int)correspondences[i].x, (int)correspondences[i].y,
+                (int)correspondences[i].rx + x, (int)correspondences[i].ry + y,
+                width, height))
+          continue;
+        subimage_norm =
+            compute_variance(ref, ref_stride, (int)correspondences[i].rx + x,
+                             (int)correspondences[i].ry + y, NULL);
+        match_ncc = compute_cross_correlation(
+                        frm, frm_stride, (int)correspondences[i].x,
+                        (int)correspondences[i].y, ref, ref_stride,
+                        (int)correspondences[i].rx + x,
+                        (int)correspondences[i].ry + y) /
+                    sqrt(template_norm * subimage_norm);
+        if (match_ncc > best_match_ncc) {
+          best_match_ncc = match_ncc;
+          best_y = y;
+          best_x = x;
+        }
+      }
+    }
+    correspondences[i].rx += (double)best_x;
+    correspondences[i].ry += (double)best_y;
+  }
+  for (i = 0; i < num_correspondences; ++i) {
+    double template_norm =
+        compute_variance(ref, ref_stride, (int)correspondences[i].rx,
+                         (int)correspondences[i].ry, NULL);
+    int x, y, best_x = 0, best_y = 0;
+    double best_match_ncc = 0.0;
+    for (y = -SEARCH_SZ_BY2; y <= SEARCH_SZ_BY2; ++y)
+      for (x = -SEARCH_SZ_BY2; x <= SEARCH_SZ_BY2; ++x) {
+        double match_ncc;
+        double subimage_norm;
+        if (!is_eligible_point((int)correspondences[i].x + x,
+                               (int)correspondences[i].y + y, width, height))
+          continue;
+        if (!is_eligible_distance((int)correspondences[i].x + x,
+                                  (int)correspondences[i].y + y,
+                                  (int)correspondences[i].rx,
+                                  (int)correspondences[i].ry, width, height))
+          continue;
+        subimage_norm =
+            compute_variance(frm, frm_stride, (int)correspondences[i].x + x,
+                             (int)correspondences[i].y + y, NULL);
+        match_ncc =
+            compute_cross_correlation(
+                frm, frm_stride, (int)correspondences[i].x + x,
+                (int)correspondences[i].y + y, ref, ref_stride,
+                (int)correspondences[i].rx, (int)correspondences[i].ry) /
+            sqrt(template_norm * subimage_norm);
+        if (match_ncc > best_match_ncc) {
+          best_match_ncc = match_ncc;
+          best_y = y;
+          best_x = x;
+        }
+      }
+    correspondences[i].x += best_x;
+    correspondences[i].y += best_y;
+  }
+}
+
+int determine_correspondence(unsigned char *frm, int *frm_corners,
+                             int num_frm_corners, unsigned char *ref,
+                             int *ref_corners, int num_ref_corners, int width,
+                             int height, int frm_stride, int ref_stride,
+                             double *correspondence_pts) {
+  // TODO(sarahparker) Improve this to include 2-way match
+  int i, j;
+  correspondence *correspondences = (correspondence *)correspondence_pts;
+  int num_correspondences = 0;
+  for (i = 0; i < num_frm_corners; ++i) {
+    double best_match_ncc = 0.0;
+    double template_norm;
+    int best_match_j = -1;
+    if (!is_eligible_point(frm_corners[2 * i], frm_corners[2 * i + 1], width,
+                           height))
+      continue;
+    template_norm = compute_variance(frm, frm_stride, frm_corners[2 * i],
+                                     frm_corners[2 * i + 1], NULL);
+    for (j = 0; j < num_ref_corners; ++j) {
+      double match_ncc;
+      double subimage_norm;
+      if (!is_eligible_point(ref_corners[2 * j], ref_corners[2 * j + 1], width,
+                             height))
+        continue;
+      if (!is_eligible_distance(frm_corners[2 * i], frm_corners[2 * i + 1],
+                                ref_corners[2 * j], ref_corners[2 * j + 1],
+                                width, height))
+        continue;
+      subimage_norm = compute_variance(ref, ref_stride, ref_corners[2 * j],
+                                       ref_corners[2 * j + 1], NULL);
+      match_ncc = compute_cross_correlation(frm, frm_stride, frm_corners[2 * i],
+                                            frm_corners[2 * i + 1], ref,
+                                            ref_stride, ref_corners[2 * j],
+                                            ref_corners[2 * j + 1]) /
+                  sqrt(template_norm * subimage_norm);
+      if (match_ncc > best_match_ncc) {
+        best_match_ncc = match_ncc;
+        best_match_j = j;
+      }
+    }
+    if (best_match_ncc > THRESHOLD_NCC) {
+      correspondences[num_correspondences].x = (double)frm_corners[2 * i];
+      correspondences[num_correspondences].y = (double)frm_corners[2 * i + 1];
+      correspondences[num_correspondences].rx =
+          (double)ref_corners[2 * best_match_j];
+      correspondences[num_correspondences].ry =
+          (double)ref_corners[2 * best_match_j + 1];
+      num_correspondences++;
+    }
+  }
+  improve_correspondence(frm, ref, width, height, frm_stride, ref_stride,
+                         correspondences, num_correspondences);
+  return num_correspondences;
+}

vp10/encoder/corner_match.h

+/*
+ *  Copyright (c) 2010 The WebM project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#ifndef VP10_ENCODER_CORNER_MATCH_H_
+#define VP10_ENCODER_CORNER_MATCH_H_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <memory.h>
+
+typedef struct {
+  double x, y;
+  double rx, ry;