1 /*
2 * Copyright (c) 2018, 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 <math.h>
13
14 #include <immintrin.h>
15 #include "config/aom_dsp_rtcd.h"
16
17 #include "aom_ports/mem.h"
18 #include "aom_dsp/flow_estimation/corner_match.h"
19
20 DECLARE_ALIGNED(32, static const uint16_t, ones_array[16]) = { 1, 1, 1, 1, 1, 1,
21 1, 1, 1, 1, 1, 1,
22 1, 1, 1, 1 };
23
24 #if MATCH_SZ != 16
25 #error "Need to apply pixel mask in corner_match_avx2.c if MATCH_SZ != 16"
26 #endif
27
28 /* Compute mean and standard deviation of pixels in a window of size
29 MATCH_SZ by MATCH_SZ centered at (x, y).
30 Store results into *mean and *one_over_stddev
31
32 Note: The output of this function is scaled by MATCH_SZ, as in
33 *mean = MATCH_SZ * <true mean> and
34 *one_over_stddev = 1 / (MATCH_SZ * <true stddev>)
35
36 Combined with the fact that we return 1/stddev rather than the standard
37 deviation itself, this allows us to completely avoid divisions in
38 aom_compute_correlation, which is much hotter than this function is.
39
40 Returns true if this feature point is usable, false otherwise.
41 */
aom_compute_mean_stddev_avx2(const unsigned char * frame,int stride,int x,int y,double * mean,double * one_over_stddev)42 bool aom_compute_mean_stddev_avx2(const unsigned char *frame, int stride, int x,
43 int y, double *mean,
44 double *one_over_stddev) {
45 __m256i sum_vec = _mm256_setzero_si256();
46 __m256i sumsq_vec = _mm256_setzero_si256();
47
48 frame += (y - MATCH_SZ_BY2) * stride + (x - MATCH_SZ_BY2);
49
50 for (int i = 0; i < MATCH_SZ; ++i) {
51 const __m256i v = _mm256_cvtepu8_epi16(_mm_loadu_si128((__m128i *)frame));
52
53 sum_vec = _mm256_add_epi16(sum_vec, v);
54 sumsq_vec = _mm256_add_epi32(sumsq_vec, _mm256_madd_epi16(v, v));
55
56 frame += stride;
57 }
58
59 // Reduce sum_vec and sumsq_vec into single values
60 // Start by reducing each vector to 8x32-bit values, hadd() to perform 8
61 // additions, sum vertically to do 4 more, then the last 2 in scalar code.
62 const __m256i ones = _mm256_load_si256((__m256i *)ones_array);
63 const __m256i partial_sum = _mm256_madd_epi16(sum_vec, ones);
64 const __m256i tmp_8x32 = _mm256_hadd_epi32(partial_sum, sumsq_vec);
65 const __m128i tmp_4x32 = _mm_add_epi32(_mm256_extracti128_si256(tmp_8x32, 0),
66 _mm256_extracti128_si256(tmp_8x32, 1));
67 const int sum =
68 _mm_extract_epi32(tmp_4x32, 0) + _mm_extract_epi32(tmp_4x32, 1);
69 const int sumsq =
70 _mm_extract_epi32(tmp_4x32, 2) + _mm_extract_epi32(tmp_4x32, 3);
71
72 *mean = (double)sum / MATCH_SZ;
73 const double variance = sumsq - (*mean) * (*mean);
74 if (variance < MIN_FEATURE_VARIANCE) {
75 *one_over_stddev = 0.0;
76 return false;
77 }
78 *one_over_stddev = 1.0 / sqrt(variance);
79 return true;
80 }
81
82 /* Compute corr(frame1, frame2) over a window of size MATCH_SZ by MATCH_SZ.
83 To save on computation, the mean and (1 divided by the) standard deviation
84 of the window in each frame are precomputed and passed into this function
85 as arguments.
86 */
aom_compute_correlation_avx2(const unsigned char * frame1,int stride1,int x1,int y1,double mean1,double one_over_stddev1,const unsigned char * frame2,int stride2,int x2,int y2,double mean2,double one_over_stddev2)87 double aom_compute_correlation_avx2(const unsigned char *frame1, int stride1,
88 int x1, int y1, double mean1,
89 double one_over_stddev1,
90 const unsigned char *frame2, int stride2,
91 int x2, int y2, double mean2,
92 double one_over_stddev2) {
93 __m256i cross_vec = _mm256_setzero_si256();
94
95 frame1 += (y1 - MATCH_SZ_BY2) * stride1 + (x1 - MATCH_SZ_BY2);
96 frame2 += (y2 - MATCH_SZ_BY2) * stride2 + (x2 - MATCH_SZ_BY2);
97
98 for (int i = 0; i < MATCH_SZ; ++i) {
99 const __m256i v1 = _mm256_cvtepu8_epi16(_mm_loadu_si128((__m128i *)frame1));
100 const __m256i v2 = _mm256_cvtepu8_epi16(_mm_loadu_si128((__m128i *)frame2));
101
102 cross_vec = _mm256_add_epi32(cross_vec, _mm256_madd_epi16(v1, v2));
103
104 frame1 += stride1;
105 frame2 += stride2;
106 }
107
108 // Sum cross_vec into a single value
109 const __m128i tmp = _mm_add_epi32(_mm256_extracti128_si256(cross_vec, 0),
110 _mm256_extracti128_si256(cross_vec, 1));
111 const int cross = _mm_extract_epi32(tmp, 0) + _mm_extract_epi32(tmp, 1) +
112 _mm_extract_epi32(tmp, 2) + _mm_extract_epi32(tmp, 3);
113
114 // Note: In theory, the calculations here "should" be
115 // covariance = cross / N^2 - mean1 * mean2
116 // correlation = covariance / (stddev1 * stddev2).
117 //
118 // However, because of the scaling in aom_compute_mean_stddev, the
119 // lines below actually calculate
120 // covariance * N^2 = cross - (mean1 * N) * (mean2 * N)
121 // correlation = (covariance * N^2) / ((stddev1 * N) * (stddev2 * N))
122 //
123 // ie. we have removed the need for a division, and still end up with the
124 // correct unscaled correlation (ie, in the range [-1, +1])
125 const double covariance = cross - mean1 * mean2;
126 const double correlation = covariance * (one_over_stddev1 * one_over_stddev2);
127 return correlation;
128 }
129