1 /*
2  * Copyright (c) 2016, Alliance for Open Media. All rights reserved
3  *
4  * This source code is subject to the terms of the BSD 2 Clause License and
5  * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
6  * was not distributed with this source code in the LICENSE file, you can
7  * obtain it at www.aomedia.org/license/software. If the Alliance for Open
8  * Media Patent License 1.0 was not distributed with this source code in the
9  * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
10  */
11 
12 #include <stdlib.h>
13 #include <memory.h>
14 #include <math.h>
15 
16 #include "config/av1_rtcd.h"
17 
18 #include "aom_ports/system_state.h"
19 #include "av1/encoder/corner_match.h"
20 
21 #define SEARCH_SZ 9
22 #define SEARCH_SZ_BY2 ((SEARCH_SZ - 1) / 2)
23 
24 #define THRESHOLD_NCC 0.75
25 
26 /* Compute var(im) * MATCH_SZ_SQ over a MATCH_SZ by MATCH_SZ window of im,
27    centered at (x, y).
28 */
compute_variance(unsigned char * im,int stride,int x,int y)29 static double compute_variance(unsigned char *im, int stride, int x, int y) {
30   int sum = 0;
31   int sumsq = 0;
32   int var;
33   int i, j;
34   for (i = 0; i < MATCH_SZ; ++i)
35     for (j = 0; j < MATCH_SZ; ++j) {
36       sum += im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)];
37       sumsq += im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)] *
38                im[(i + y - MATCH_SZ_BY2) * stride + (j + x - MATCH_SZ_BY2)];
39     }
40   var = sumsq * MATCH_SZ_SQ - sum * sum;
41   return (double)var;
42 }
43 
44 /* Compute corr(im1, im2) * MATCH_SZ * stddev(im1), where the
45    correlation/standard deviation are taken over MATCH_SZ by MATCH_SZ windows
46    of each image, centered at (x1, y1) and (x2, y2) respectively.
47 */
av1_compute_cross_correlation_c(unsigned char * im1,int stride1,int x1,int y1,unsigned char * im2,int stride2,int x2,int y2)48 double av1_compute_cross_correlation_c(unsigned char *im1, int stride1, int x1,
49                                        int y1, unsigned char *im2, int stride2,
50                                        int x2, int y2) {
51   int v1, v2;
52   int sum1 = 0;
53   int sum2 = 0;
54   int sumsq2 = 0;
55   int cross = 0;
56   int var2, cov;
57   int i, j;
58   for (i = 0; i < MATCH_SZ; ++i)
59     for (j = 0; j < MATCH_SZ; ++j) {
60       v1 = im1[(i + y1 - MATCH_SZ_BY2) * stride1 + (j + x1 - MATCH_SZ_BY2)];
61       v2 = im2[(i + y2 - MATCH_SZ_BY2) * stride2 + (j + x2 - MATCH_SZ_BY2)];
62       sum1 += v1;
63       sum2 += v2;
64       sumsq2 += v2 * v2;
65       cross += v1 * v2;
66     }
67   var2 = sumsq2 * MATCH_SZ_SQ - sum2 * sum2;
68   cov = cross * MATCH_SZ_SQ - sum1 * sum2;
69   aom_clear_system_state();
70   return cov / sqrt((double)var2);
71 }
72 
is_eligible_point(int pointx,int pointy,int width,int height)73 static int is_eligible_point(int pointx, int pointy, int width, int height) {
74   return (pointx >= MATCH_SZ_BY2 && pointy >= MATCH_SZ_BY2 &&
75           pointx + MATCH_SZ_BY2 < width && pointy + MATCH_SZ_BY2 < height);
76 }
77 
is_eligible_distance(int point1x,int point1y,int point2x,int point2y,int width,int height)78 static int is_eligible_distance(int point1x, int point1y, int point2x,
79                                 int point2y, int width, int height) {
80   const int thresh = (width < height ? height : width) >> 4;
81   return ((point1x - point2x) * (point1x - point2x) +
82           (point1y - point2y) * (point1y - point2y)) <= thresh * thresh;
83 }
84 
improve_correspondence(unsigned char * frm,unsigned char * ref,int width,int height,int frm_stride,int ref_stride,Correspondence * correspondences,int num_correspondences)85 static void improve_correspondence(unsigned char *frm, unsigned char *ref,
86                                    int width, int height, int frm_stride,
87                                    int ref_stride,
88                                    Correspondence *correspondences,
89                                    int num_correspondences) {
90   int i;
91   for (i = 0; i < num_correspondences; ++i) {
92     int x, y, best_x = 0, best_y = 0;
93     double best_match_ncc = 0.0;
94     for (y = -SEARCH_SZ_BY2; y <= SEARCH_SZ_BY2; ++y) {
95       for (x = -SEARCH_SZ_BY2; x <= SEARCH_SZ_BY2; ++x) {
96         double match_ncc;
97         if (!is_eligible_point(correspondences[i].rx + x,
98                                correspondences[i].ry + y, width, height))
99           continue;
100         if (!is_eligible_distance(correspondences[i].x, correspondences[i].y,
101                                   correspondences[i].rx + x,
102                                   correspondences[i].ry + y, width, height))
103           continue;
104         match_ncc = av1_compute_cross_correlation(
105             frm, frm_stride, correspondences[i].x, correspondences[i].y, ref,
106             ref_stride, correspondences[i].rx + x, correspondences[i].ry + y);
107         if (match_ncc > best_match_ncc) {
108           best_match_ncc = match_ncc;
109           best_y = y;
110           best_x = x;
111         }
112       }
113     }
114     correspondences[i].rx += best_x;
115     correspondences[i].ry += best_y;
116   }
117   for (i = 0; i < num_correspondences; ++i) {
118     int x, y, best_x = 0, best_y = 0;
119     double best_match_ncc = 0.0;
120     for (y = -SEARCH_SZ_BY2; y <= SEARCH_SZ_BY2; ++y)
121       for (x = -SEARCH_SZ_BY2; x <= SEARCH_SZ_BY2; ++x) {
122         double match_ncc;
123         if (!is_eligible_point(correspondences[i].x + x,
124                                correspondences[i].y + y, width, height))
125           continue;
126         if (!is_eligible_distance(
127                 correspondences[i].x + x, correspondences[i].y + y,
128                 correspondences[i].rx, correspondences[i].ry, width, height))
129           continue;
130         match_ncc = av1_compute_cross_correlation(
131             ref, ref_stride, correspondences[i].rx, correspondences[i].ry, frm,
132             frm_stride, correspondences[i].x + x, correspondences[i].y + y);
133         if (match_ncc > best_match_ncc) {
134           best_match_ncc = match_ncc;
135           best_y = y;
136           best_x = x;
137         }
138       }
139     correspondences[i].x += best_x;
140     correspondences[i].y += best_y;
141   }
142 }
143 
av1_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,int * correspondence_pts)144 int av1_determine_correspondence(unsigned char *frm, int *frm_corners,
145                                  int num_frm_corners, unsigned char *ref,
146                                  int *ref_corners, int num_ref_corners,
147                                  int width, int height, int frm_stride,
148                                  int ref_stride, int *correspondence_pts) {
149   // TODO(sarahparker) Improve this to include 2-way match
150   int i, j;
151   Correspondence *correspondences = (Correspondence *)correspondence_pts;
152   int num_correspondences = 0;
153   for (i = 0; i < num_frm_corners; ++i) {
154     double best_match_ncc = 0.0;
155     double template_norm;
156     int best_match_j = -1;
157     if (!is_eligible_point(frm_corners[2 * i], frm_corners[2 * i + 1], width,
158                            height))
159       continue;
160     for (j = 0; j < num_ref_corners; ++j) {
161       double match_ncc;
162       if (!is_eligible_point(ref_corners[2 * j], ref_corners[2 * j + 1], width,
163                              height))
164         continue;
165       if (!is_eligible_distance(frm_corners[2 * i], frm_corners[2 * i + 1],
166                                 ref_corners[2 * j], ref_corners[2 * j + 1],
167                                 width, height))
168         continue;
169       match_ncc = av1_compute_cross_correlation(
170           frm, frm_stride, frm_corners[2 * i], frm_corners[2 * i + 1], ref,
171           ref_stride, ref_corners[2 * j], ref_corners[2 * j + 1]);
172       if (match_ncc > best_match_ncc) {
173         best_match_ncc = match_ncc;
174         best_match_j = j;
175       }
176     }
177     // Note: We want to test if the best correlation is >= THRESHOLD_NCC,
178     // but need to account for the normalization in
179     // av1_compute_cross_correlation.
180     template_norm = compute_variance(frm, frm_stride, frm_corners[2 * i],
181                                      frm_corners[2 * i + 1]);
182     if (best_match_ncc > THRESHOLD_NCC * sqrt(template_norm)) {
183       correspondences[num_correspondences].x = frm_corners[2 * i];
184       correspondences[num_correspondences].y = frm_corners[2 * i + 1];
185       correspondences[num_correspondences].rx = ref_corners[2 * best_match_j];
186       correspondences[num_correspondences].ry =
187           ref_corners[2 * best_match_j + 1];
188       num_correspondences++;
189     }
190   }
191   improve_correspondence(frm, ref, width, height, frm_stride, ref_stride,
192                          correspondences, num_correspondences);
193   return num_correspondences;
194 }
195