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 <immintrin.h> // AVX2
13 #include "aom_dsp/x86/synonyms.h"
14 #include "aom_dsp/x86/synonyms_avx2.h"
15 #include "aom_dsp/x86/transpose_sse2.h"
16
17 #include "config/av1_rtcd.h"
18 #include "av1/common/restoration.h"
19 #include "av1/encoder/pickrst.h"
20
acc_stat_avx2(int32_t * dst,const uint8_t * src,const __m128i * shuffle,const __m256i * kl)21 static INLINE void acc_stat_avx2(int32_t *dst, const uint8_t *src,
22 const __m128i *shuffle, const __m256i *kl) {
23 const __m128i s = _mm_shuffle_epi8(xx_loadu_128(src), *shuffle);
24 const __m256i d0 = _mm256_madd_epi16(*kl, _mm256_cvtepu8_epi16(s));
25 const __m256i dst0 = yy_load_256(dst);
26 const __m256i r0 = _mm256_add_epi32(dst0, d0);
27 yy_store_256(dst, r0);
28 }
29
acc_stat_win7_one_line_avx2(const uint8_t * dgd,const uint8_t * src,int h_start,int h_end,int dgd_stride,const __m128i * shuffle,int32_t * sumX,int32_t sumY[WIENER_WIN][WIENER_WIN],int32_t M_int[WIENER_WIN][WIENER_WIN],int32_t H_int[WIENER_WIN2][WIENER_WIN * 8])30 static INLINE void acc_stat_win7_one_line_avx2(
31 const uint8_t *dgd, const uint8_t *src, int h_start, int h_end,
32 int dgd_stride, const __m128i *shuffle, int32_t *sumX,
33 int32_t sumY[WIENER_WIN][WIENER_WIN], int32_t M_int[WIENER_WIN][WIENER_WIN],
34 int32_t H_int[WIENER_WIN2][WIENER_WIN * 8]) {
35 int j, k, l;
36 const int wiener_win = WIENER_WIN;
37 for (j = h_start; j < h_end; j += 2) {
38 const uint8_t X1 = src[j];
39 const uint8_t X2 = src[j + 1];
40 *sumX += X1 + X2;
41 const uint8_t *dgd_ij = dgd + j;
42 for (k = 0; k < wiener_win; k++) {
43 const uint8_t *dgd_ijk = dgd_ij + k * dgd_stride;
44 for (l = 0; l < wiener_win; l++) {
45 int32_t *H_ = &H_int[(l * wiener_win + k)][0];
46 const uint8_t D1 = dgd_ijk[l];
47 const uint8_t D2 = dgd_ijk[l + 1];
48 sumY[k][l] += D1 + D2;
49 M_int[k][l] += D1 * X1 + D2 * X2;
50
51 const __m256i kl =
52 _mm256_cvtepu8_epi16(_mm_set1_epi16(*((uint16_t *)(dgd_ijk + l))));
53 acc_stat_avx2(H_ + 0 * 8, dgd_ij + 0 * dgd_stride, shuffle, &kl);
54 acc_stat_avx2(H_ + 1 * 8, dgd_ij + 1 * dgd_stride, shuffle, &kl);
55 acc_stat_avx2(H_ + 2 * 8, dgd_ij + 2 * dgd_stride, shuffle, &kl);
56 acc_stat_avx2(H_ + 3 * 8, dgd_ij + 3 * dgd_stride, shuffle, &kl);
57 acc_stat_avx2(H_ + 4 * 8, dgd_ij + 4 * dgd_stride, shuffle, &kl);
58 acc_stat_avx2(H_ + 5 * 8, dgd_ij + 5 * dgd_stride, shuffle, &kl);
59 acc_stat_avx2(H_ + 6 * 8, dgd_ij + 6 * dgd_stride, shuffle, &kl);
60 }
61 }
62 }
63 }
64
compute_stats_win7_opt_avx2(const uint8_t * dgd,const uint8_t * src,int h_start,int h_end,int v_start,int v_end,int dgd_stride,int src_stride,int64_t * M,int64_t * H)65 static INLINE void compute_stats_win7_opt_avx2(
66 const uint8_t *dgd, const uint8_t *src, int h_start, int h_end, int v_start,
67 int v_end, int dgd_stride, int src_stride, int64_t *M, int64_t *H) {
68 int i, j, k, l, m, n;
69 const int wiener_win = WIENER_WIN;
70 const int pixel_count = (h_end - h_start) * (v_end - v_start);
71 const int wiener_win2 = wiener_win * wiener_win;
72 const int wiener_halfwin = (wiener_win >> 1);
73 uint8_t avg = find_average(dgd, h_start, h_end, v_start, v_end, dgd_stride);
74
75 int32_t M_int32[WIENER_WIN][WIENER_WIN] = { { 0 } };
76 int64_t M_int64[WIENER_WIN][WIENER_WIN] = { { 0 } };
77
78 DECLARE_ALIGNED(32, int32_t,
79 H_int32[WIENER_WIN2][WIENER_WIN * 8]) = { { 0 } };
80 int64_t H_int64[WIENER_WIN2][WIENER_WIN * 8] = { { 0 } };
81 int32_t sumY[WIENER_WIN][WIENER_WIN] = { { 0 } };
82 int32_t sumX = 0;
83 const uint8_t *dgd_win = dgd - wiener_halfwin * dgd_stride - wiener_halfwin;
84
85 const __m128i shuffle = xx_loadu_128(g_shuffle_stats_data);
86 for (j = v_start; j < v_end; j += 64) {
87 const int vert_end = AOMMIN(64, v_end - j) + j;
88 for (i = j; i < vert_end; i++) {
89 acc_stat_win7_one_line_avx2(
90 dgd_win + i * dgd_stride, src + i * src_stride, h_start, h_end,
91 dgd_stride, &shuffle, &sumX, sumY, M_int32, H_int32);
92 }
93 for (k = 0; k < wiener_win; ++k) {
94 for (l = 0; l < wiener_win; ++l) {
95 M_int64[k][l] += M_int32[k][l];
96 M_int32[k][l] = 0;
97 }
98 }
99 for (k = 0; k < WIENER_WIN2; ++k) {
100 for (l = 0; l < WIENER_WIN * 8; ++l) {
101 H_int64[k][l] += H_int32[k][l];
102 H_int32[k][l] = 0;
103 }
104 }
105 }
106
107 const int64_t avg_square_sum = (int64_t)avg * (int64_t)avg * pixel_count;
108 for (k = 0; k < wiener_win; k++) {
109 for (l = 0; l < wiener_win; l++) {
110 const int32_t idx0 = l * wiener_win + k;
111 M[idx0] =
112 M_int64[k][l] + (avg_square_sum - (int64_t)avg * (sumX + sumY[k][l]));
113 int64_t *H_ = H + idx0 * wiener_win2;
114 int64_t *H_int_ = &H_int64[idx0][0];
115 for (m = 0; m < wiener_win; m++) {
116 for (n = 0; n < wiener_win; n++) {
117 H_[m * wiener_win + n] = H_int_[n * 8 + m] + avg_square_sum -
118 (int64_t)avg * (sumY[k][l] + sumY[n][m]);
119 }
120 }
121 }
122 }
123 }
124
125 #if CONFIG_AV1_HIGHBITDEPTH
acc_stat_highbd_avx2(int64_t * dst,const uint16_t * dgd,const __m256i * shuffle,const __m256i * dgd_ijkl)126 static INLINE void acc_stat_highbd_avx2(int64_t *dst, const uint16_t *dgd,
127 const __m256i *shuffle,
128 const __m256i *dgd_ijkl) {
129 // Load two 128-bit chunks from dgd
130 const __m256i s0 = _mm256_inserti128_si256(
131 _mm256_castsi128_si256(_mm_loadu_si128((__m128i *)dgd)),
132 _mm_loadu_si128((__m128i *)(dgd + 4)), 1);
133 // s0 = [11 10 9 8 7 6 5 4] [7 6 5 4 3 2 1 0] as u16 (values are dgd indices)
134 // The weird order is so the shuffle stays within 128-bit lanes
135
136 // Shuffle 16x u16 values within lanes according to the mask:
137 // [0 1 1 2 2 3 3 4] [0 1 1 2 2 3 3 4]
138 // (Actually we shuffle u8 values as there's no 16-bit shuffle)
139 const __m256i s1 = _mm256_shuffle_epi8(s0, *shuffle);
140 // s1 = [8 7 7 6 6 5 5 4] [4 3 3 2 2 1 1 0] as u16 (values are dgd indices)
141
142 // Multiply 16x 16-bit integers in dgd_ijkl and s1, resulting in 16x 32-bit
143 // integers then horizontally add pairs of these integers resulting in 8x
144 // 32-bit integers
145 const __m256i d0 = _mm256_madd_epi16(*dgd_ijkl, s1);
146 // d0 = [a b c d] [e f g h] as u32
147
148 // Take the lower-half of d0, extend to u64, add it on to dst (H)
149 const __m256i d0l = _mm256_cvtepu32_epi64(_mm256_extracti128_si256(d0, 0));
150 // d0l = [a b] [c d] as u64
151 const __m256i dst0 = yy_load_256(dst);
152 yy_store_256(dst, _mm256_add_epi64(d0l, dst0));
153
154 // Take the upper-half of d0, extend to u64, add it on to dst (H)
155 const __m256i d0h = _mm256_cvtepu32_epi64(_mm256_extracti128_si256(d0, 1));
156 // d0h = [e f] [g h] as u64
157 const __m256i dst1 = yy_load_256(dst + 4);
158 yy_store_256(dst + 4, _mm256_add_epi64(d0h, dst1));
159 }
160
acc_stat_highbd_win7_one_line_avx2(const uint16_t * dgd,const uint16_t * src,int h_start,int h_end,int dgd_stride,const __m256i * shuffle,int32_t * sumX,int32_t sumY[WIENER_WIN][WIENER_WIN],int64_t M_int[WIENER_WIN][WIENER_WIN],int64_t H_int[WIENER_WIN2][WIENER_WIN * 8])161 static INLINE void acc_stat_highbd_win7_one_line_avx2(
162 const uint16_t *dgd, const uint16_t *src, int h_start, int h_end,
163 int dgd_stride, const __m256i *shuffle, int32_t *sumX,
164 int32_t sumY[WIENER_WIN][WIENER_WIN], int64_t M_int[WIENER_WIN][WIENER_WIN],
165 int64_t H_int[WIENER_WIN2][WIENER_WIN * 8]) {
166 int j, k, l;
167 const int wiener_win = WIENER_WIN;
168 for (j = h_start; j < h_end; j += 2) {
169 const uint16_t X1 = src[j];
170 const uint16_t X2 = src[j + 1];
171 *sumX += X1 + X2;
172 const uint16_t *dgd_ij = dgd + j;
173 for (k = 0; k < wiener_win; k++) {
174 const uint16_t *dgd_ijk = dgd_ij + k * dgd_stride;
175 for (l = 0; l < wiener_win; l++) {
176 int64_t *H_ = &H_int[(l * wiener_win + k)][0];
177 const uint16_t D1 = dgd_ijk[l];
178 const uint16_t D2 = dgd_ijk[l + 1];
179 sumY[k][l] += D1 + D2;
180 M_int[k][l] += D1 * X1 + D2 * X2;
181
182 // Load two u16 values from dgd_ijkl combined as a u32,
183 // then broadcast to 8x u32 slots of a 256
184 const __m256i dgd_ijkl =
185 _mm256_set1_epi32(*((uint32_t *)(dgd_ijk + l)));
186 // dgd_ijkl = [y x y x y x y x] [y x y x y x y x] where each is a u16
187
188 acc_stat_highbd_avx2(H_ + 0 * 8, dgd_ij + 0 * dgd_stride, shuffle,
189 &dgd_ijkl);
190 acc_stat_highbd_avx2(H_ + 1 * 8, dgd_ij + 1 * dgd_stride, shuffle,
191 &dgd_ijkl);
192 acc_stat_highbd_avx2(H_ + 2 * 8, dgd_ij + 2 * dgd_stride, shuffle,
193 &dgd_ijkl);
194 acc_stat_highbd_avx2(H_ + 3 * 8, dgd_ij + 3 * dgd_stride, shuffle,
195 &dgd_ijkl);
196 acc_stat_highbd_avx2(H_ + 4 * 8, dgd_ij + 4 * dgd_stride, shuffle,
197 &dgd_ijkl);
198 acc_stat_highbd_avx2(H_ + 5 * 8, dgd_ij + 5 * dgd_stride, shuffle,
199 &dgd_ijkl);
200 acc_stat_highbd_avx2(H_ + 6 * 8, dgd_ij + 6 * dgd_stride, shuffle,
201 &dgd_ijkl);
202 }
203 }
204 }
205 }
206
compute_stats_highbd_win7_opt_avx2(const uint8_t * dgd8,const uint8_t * src8,int h_start,int h_end,int v_start,int v_end,int dgd_stride,int src_stride,int64_t * M,int64_t * H,aom_bit_depth_t bit_depth)207 static INLINE void compute_stats_highbd_win7_opt_avx2(
208 const uint8_t *dgd8, const uint8_t *src8, int h_start, int h_end,
209 int v_start, int v_end, int dgd_stride, int src_stride, int64_t *M,
210 int64_t *H, aom_bit_depth_t bit_depth) {
211 int i, j, k, l, m, n;
212 const int wiener_win = WIENER_WIN;
213 const int pixel_count = (h_end - h_start) * (v_end - v_start);
214 const int wiener_win2 = wiener_win * wiener_win;
215 const int wiener_halfwin = (wiener_win >> 1);
216 const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
217 const uint16_t *dgd = CONVERT_TO_SHORTPTR(dgd8);
218 const uint16_t avg =
219 find_average_highbd(dgd, h_start, h_end, v_start, v_end, dgd_stride);
220
221 int64_t M_int[WIENER_WIN][WIENER_WIN] = { { 0 } };
222 DECLARE_ALIGNED(32, int64_t, H_int[WIENER_WIN2][WIENER_WIN * 8]) = { { 0 } };
223 int32_t sumY[WIENER_WIN][WIENER_WIN] = { { 0 } };
224 int32_t sumX = 0;
225 const uint16_t *dgd_win = dgd - wiener_halfwin * dgd_stride - wiener_halfwin;
226
227 const __m256i shuffle = yy_loadu_256(g_shuffle_stats_highbd_data);
228 for (j = v_start; j < v_end; j += 64) {
229 const int vert_end = AOMMIN(64, v_end - j) + j;
230 for (i = j; i < vert_end; i++) {
231 acc_stat_highbd_win7_one_line_avx2(
232 dgd_win + i * dgd_stride, src + i * src_stride, h_start, h_end,
233 dgd_stride, &shuffle, &sumX, sumY, M_int, H_int);
234 }
235 }
236
237 uint8_t bit_depth_divider = 1;
238 if (bit_depth == AOM_BITS_12)
239 bit_depth_divider = 16;
240 else if (bit_depth == AOM_BITS_10)
241 bit_depth_divider = 4;
242
243 const int64_t avg_square_sum = (int64_t)avg * (int64_t)avg * pixel_count;
244 for (k = 0; k < wiener_win; k++) {
245 for (l = 0; l < wiener_win; l++) {
246 const int32_t idx0 = l * wiener_win + k;
247 M[idx0] = (M_int[k][l] +
248 (avg_square_sum - (int64_t)avg * (sumX + sumY[k][l]))) /
249 bit_depth_divider;
250 int64_t *H_ = H + idx0 * wiener_win2;
251 int64_t *H_int_ = &H_int[idx0][0];
252 for (m = 0; m < wiener_win; m++) {
253 for (n = 0; n < wiener_win; n++) {
254 H_[m * wiener_win + n] =
255 (H_int_[n * 8 + m] +
256 (avg_square_sum - (int64_t)avg * (sumY[k][l] + sumY[n][m]))) /
257 bit_depth_divider;
258 }
259 }
260 }
261 }
262 }
263
acc_stat_highbd_win5_one_line_avx2(const uint16_t * dgd,const uint16_t * src,int h_start,int h_end,int dgd_stride,const __m256i * shuffle,int32_t * sumX,int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],int64_t M_int[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],int64_t H_int[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8])264 static INLINE void acc_stat_highbd_win5_one_line_avx2(
265 const uint16_t *dgd, const uint16_t *src, int h_start, int h_end,
266 int dgd_stride, const __m256i *shuffle, int32_t *sumX,
267 int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],
268 int64_t M_int[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],
269 int64_t H_int[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8]) {
270 int j, k, l;
271 const int wiener_win = WIENER_WIN_CHROMA;
272 for (j = h_start; j < h_end; j += 2) {
273 const uint16_t X1 = src[j];
274 const uint16_t X2 = src[j + 1];
275 *sumX += X1 + X2;
276 const uint16_t *dgd_ij = dgd + j;
277 for (k = 0; k < wiener_win; k++) {
278 const uint16_t *dgd_ijk = dgd_ij + k * dgd_stride;
279 for (l = 0; l < wiener_win; l++) {
280 int64_t *H_ = &H_int[(l * wiener_win + k)][0];
281 const uint16_t D1 = dgd_ijk[l];
282 const uint16_t D2 = dgd_ijk[l + 1];
283 sumY[k][l] += D1 + D2;
284 M_int[k][l] += D1 * X1 + D2 * X2;
285
286 // Load two u16 values from dgd_ijkl combined as a u32,
287 // then broadcast to 8x u32 slots of a 256
288 const __m256i dgd_ijkl =
289 _mm256_set1_epi32(*((uint32_t *)(dgd_ijk + l)));
290 // dgd_ijkl = [x y x y x y x y] [x y x y x y x y] where each is a u16
291
292 acc_stat_highbd_avx2(H_ + 0 * 8, dgd_ij + 0 * dgd_stride, shuffle,
293 &dgd_ijkl);
294 acc_stat_highbd_avx2(H_ + 1 * 8, dgd_ij + 1 * dgd_stride, shuffle,
295 &dgd_ijkl);
296 acc_stat_highbd_avx2(H_ + 2 * 8, dgd_ij + 2 * dgd_stride, shuffle,
297 &dgd_ijkl);
298 acc_stat_highbd_avx2(H_ + 3 * 8, dgd_ij + 3 * dgd_stride, shuffle,
299 &dgd_ijkl);
300 acc_stat_highbd_avx2(H_ + 4 * 8, dgd_ij + 4 * dgd_stride, shuffle,
301 &dgd_ijkl);
302 }
303 }
304 }
305 }
306
compute_stats_highbd_win5_opt_avx2(const uint8_t * dgd8,const uint8_t * src8,int h_start,int h_end,int v_start,int v_end,int dgd_stride,int src_stride,int64_t * M,int64_t * H,aom_bit_depth_t bit_depth)307 static INLINE void compute_stats_highbd_win5_opt_avx2(
308 const uint8_t *dgd8, const uint8_t *src8, int h_start, int h_end,
309 int v_start, int v_end, int dgd_stride, int src_stride, int64_t *M,
310 int64_t *H, aom_bit_depth_t bit_depth) {
311 int i, j, k, l, m, n;
312 const int wiener_win = WIENER_WIN_CHROMA;
313 const int pixel_count = (h_end - h_start) * (v_end - v_start);
314 const int wiener_win2 = wiener_win * wiener_win;
315 const int wiener_halfwin = (wiener_win >> 1);
316 const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
317 const uint16_t *dgd = CONVERT_TO_SHORTPTR(dgd8);
318 const uint16_t avg =
319 find_average_highbd(dgd, h_start, h_end, v_start, v_end, dgd_stride);
320
321 int64_t M_int64[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
322 DECLARE_ALIGNED(
323 32, int64_t,
324 H_int64[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8]) = { { 0 } };
325 int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
326 int32_t sumX = 0;
327 const uint16_t *dgd_win = dgd - wiener_halfwin * dgd_stride - wiener_halfwin;
328
329 const __m256i shuffle = yy_loadu_256(g_shuffle_stats_highbd_data);
330 for (j = v_start; j < v_end; j += 64) {
331 const int vert_end = AOMMIN(64, v_end - j) + j;
332 for (i = j; i < vert_end; i++) {
333 acc_stat_highbd_win5_one_line_avx2(
334 dgd_win + i * dgd_stride, src + i * src_stride, h_start, h_end,
335 dgd_stride, &shuffle, &sumX, sumY, M_int64, H_int64);
336 }
337 }
338
339 uint8_t bit_depth_divider = 1;
340 if (bit_depth == AOM_BITS_12)
341 bit_depth_divider = 16;
342 else if (bit_depth == AOM_BITS_10)
343 bit_depth_divider = 4;
344
345 const int64_t avg_square_sum = (int64_t)avg * (int64_t)avg * pixel_count;
346 for (k = 0; k < wiener_win; k++) {
347 for (l = 0; l < wiener_win; l++) {
348 const int32_t idx0 = l * wiener_win + k;
349 M[idx0] = (M_int64[k][l] +
350 (avg_square_sum - (int64_t)avg * (sumX + sumY[k][l]))) /
351 bit_depth_divider;
352 int64_t *H_ = H + idx0 * wiener_win2;
353 int64_t *H_int_ = &H_int64[idx0][0];
354 for (m = 0; m < wiener_win; m++) {
355 for (n = 0; n < wiener_win; n++) {
356 H_[m * wiener_win + n] =
357 (H_int_[n * 8 + m] +
358 (avg_square_sum - (int64_t)avg * (sumY[k][l] + sumY[n][m]))) /
359 bit_depth_divider;
360 }
361 }
362 }
363 }
364 }
365
av1_compute_stats_highbd_avx2(int wiener_win,const uint8_t * dgd8,const uint8_t * src8,int h_start,int h_end,int v_start,int v_end,int dgd_stride,int src_stride,int64_t * M,int64_t * H,aom_bit_depth_t bit_depth)366 void av1_compute_stats_highbd_avx2(int wiener_win, const uint8_t *dgd8,
367 const uint8_t *src8, int h_start, int h_end,
368 int v_start, int v_end, int dgd_stride,
369 int src_stride, int64_t *M, int64_t *H,
370 aom_bit_depth_t bit_depth) {
371 if (wiener_win == WIENER_WIN) {
372 compute_stats_highbd_win7_opt_avx2(dgd8, src8, h_start, h_end, v_start,
373 v_end, dgd_stride, src_stride, M, H,
374 bit_depth);
375 } else if (wiener_win == WIENER_WIN_CHROMA) {
376 compute_stats_highbd_win5_opt_avx2(dgd8, src8, h_start, h_end, v_start,
377 v_end, dgd_stride, src_stride, M, H,
378 bit_depth);
379 } else {
380 av1_compute_stats_highbd_c(wiener_win, dgd8, src8, h_start, h_end, v_start,
381 v_end, dgd_stride, src_stride, M, H, bit_depth);
382 }
383 }
384 #endif // CONFIG_AV1_HIGHBITDEPTH
385
acc_stat_win5_one_line_avx2(const uint8_t * dgd,const uint8_t * src,int h_start,int h_end,int dgd_stride,const __m128i * shuffle,int32_t * sumX,int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],int32_t M_int[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],int32_t H_int[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8])386 static INLINE void acc_stat_win5_one_line_avx2(
387 const uint8_t *dgd, const uint8_t *src, int h_start, int h_end,
388 int dgd_stride, const __m128i *shuffle, int32_t *sumX,
389 int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],
390 int32_t M_int[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],
391 int32_t H_int[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8]) {
392 int j, k, l;
393 const int wiener_win = WIENER_WIN_CHROMA;
394 for (j = h_start; j < h_end; j += 2) {
395 const uint8_t X1 = src[j];
396 const uint8_t X2 = src[j + 1];
397 *sumX += X1 + X2;
398 const uint8_t *dgd_ij = dgd + j;
399 for (k = 0; k < wiener_win; k++) {
400 const uint8_t *dgd_ijk = dgd_ij + k * dgd_stride;
401 for (l = 0; l < wiener_win; l++) {
402 int32_t *H_ = &H_int[(l * wiener_win + k)][0];
403 const uint8_t D1 = dgd_ijk[l];
404 const uint8_t D2 = dgd_ijk[l + 1];
405 sumY[k][l] += D1 + D2;
406 M_int[k][l] += D1 * X1 + D2 * X2;
407
408 const __m256i kl =
409 _mm256_cvtepu8_epi16(_mm_set1_epi16(*((uint16_t *)(dgd_ijk + l))));
410 acc_stat_avx2(H_ + 0 * 8, dgd_ij + 0 * dgd_stride, shuffle, &kl);
411 acc_stat_avx2(H_ + 1 * 8, dgd_ij + 1 * dgd_stride, shuffle, &kl);
412 acc_stat_avx2(H_ + 2 * 8, dgd_ij + 2 * dgd_stride, shuffle, &kl);
413 acc_stat_avx2(H_ + 3 * 8, dgd_ij + 3 * dgd_stride, shuffle, &kl);
414 acc_stat_avx2(H_ + 4 * 8, dgd_ij + 4 * dgd_stride, shuffle, &kl);
415 }
416 }
417 }
418 }
419
compute_stats_win5_opt_avx2(const uint8_t * dgd,const uint8_t * src,int h_start,int h_end,int v_start,int v_end,int dgd_stride,int src_stride,int64_t * M,int64_t * H)420 static INLINE void compute_stats_win5_opt_avx2(
421 const uint8_t *dgd, const uint8_t *src, int h_start, int h_end, int v_start,
422 int v_end, int dgd_stride, int src_stride, int64_t *M, int64_t *H) {
423 int i, j, k, l, m, n;
424 const int wiener_win = WIENER_WIN_CHROMA;
425 const int pixel_count = (h_end - h_start) * (v_end - v_start);
426 const int wiener_win2 = wiener_win * wiener_win;
427 const int wiener_halfwin = (wiener_win >> 1);
428 uint8_t avg = find_average(dgd, h_start, h_end, v_start, v_end, dgd_stride);
429
430 int32_t M_int32[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
431 int64_t M_int64[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
432 DECLARE_ALIGNED(
433 32, int32_t,
434 H_int32[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8]) = { { 0 } };
435 int64_t H_int64[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8] = { { 0 } };
436 int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
437 int32_t sumX = 0;
438 const uint8_t *dgd_win = dgd - wiener_halfwin * dgd_stride - wiener_halfwin;
439
440 const __m128i shuffle = xx_loadu_128(g_shuffle_stats_data);
441 for (j = v_start; j < v_end; j += 64) {
442 const int vert_end = AOMMIN(64, v_end - j) + j;
443 for (i = j; i < vert_end; i++) {
444 acc_stat_win5_one_line_avx2(
445 dgd_win + i * dgd_stride, src + i * src_stride, h_start, h_end,
446 dgd_stride, &shuffle, &sumX, sumY, M_int32, H_int32);
447 }
448 for (k = 0; k < wiener_win; ++k) {
449 for (l = 0; l < wiener_win; ++l) {
450 M_int64[k][l] += M_int32[k][l];
451 M_int32[k][l] = 0;
452 }
453 }
454 for (k = 0; k < WIENER_WIN2_CHROMA; ++k) {
455 for (l = 0; l < WIENER_WIN_CHROMA * 8; ++l) {
456 H_int64[k][l] += H_int32[k][l];
457 H_int32[k][l] = 0;
458 }
459 }
460 }
461
462 const int64_t avg_square_sum = (int64_t)avg * (int64_t)avg * pixel_count;
463 for (k = 0; k < wiener_win; k++) {
464 for (l = 0; l < wiener_win; l++) {
465 const int32_t idx0 = l * wiener_win + k;
466 M[idx0] =
467 M_int64[k][l] + (avg_square_sum - (int64_t)avg * (sumX + sumY[k][l]));
468 int64_t *H_ = H + idx0 * wiener_win2;
469 int64_t *H_int_ = &H_int64[idx0][0];
470 for (m = 0; m < wiener_win; m++) {
471 for (n = 0; n < wiener_win; n++) {
472 H_[m * wiener_win + n] = H_int_[n * 8 + m] + avg_square_sum -
473 (int64_t)avg * (sumY[k][l] + sumY[n][m]);
474 }
475 }
476 }
477 }
478 }
479
av1_compute_stats_avx2(int wiener_win,const uint8_t * dgd,const uint8_t * src,int h_start,int h_end,int v_start,int v_end,int dgd_stride,int src_stride,int64_t * M,int64_t * H)480 void av1_compute_stats_avx2(int wiener_win, const uint8_t *dgd,
481 const uint8_t *src, int h_start, int h_end,
482 int v_start, int v_end, int dgd_stride,
483 int src_stride, int64_t *M, int64_t *H) {
484 if (wiener_win == WIENER_WIN) {
485 compute_stats_win7_opt_avx2(dgd, src, h_start, h_end, v_start, v_end,
486 dgd_stride, src_stride, M, H);
487 } else if (wiener_win == WIENER_WIN_CHROMA) {
488 compute_stats_win5_opt_avx2(dgd, src, h_start, h_end, v_start, v_end,
489 dgd_stride, src_stride, M, H);
490 } else {
491 av1_compute_stats_c(wiener_win, dgd, src, h_start, h_end, v_start, v_end,
492 dgd_stride, src_stride, M, H);
493 }
494 }
495
pair_set_epi16(int a,int b)496 static INLINE __m256i pair_set_epi16(int a, int b) {
497 return _mm256_set1_epi32(
498 (int32_t)(((uint16_t)(a)) | (((uint32_t)(b)) << 16)));
499 }
500
av1_lowbd_pixel_proj_error_avx2(const uint8_t * src8,int width,int height,int src_stride,const uint8_t * dat8,int dat_stride,int32_t * flt0,int flt0_stride,int32_t * flt1,int flt1_stride,int xq[2],const sgr_params_type * params)501 int64_t av1_lowbd_pixel_proj_error_avx2(
502 const uint8_t *src8, int width, int height, int src_stride,
503 const uint8_t *dat8, int dat_stride, int32_t *flt0, int flt0_stride,
504 int32_t *flt1, int flt1_stride, int xq[2], const sgr_params_type *params) {
505 int i, j, k;
506 const int32_t shift = SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS;
507 const __m256i rounding = _mm256_set1_epi32(1 << (shift - 1));
508 __m256i sum64 = _mm256_setzero_si256();
509 const uint8_t *src = src8;
510 const uint8_t *dat = dat8;
511 int64_t err = 0;
512 if (params->r[0] > 0 && params->r[1] > 0) {
513 __m256i xq_coeff = pair_set_epi16(xq[0], xq[1]);
514 for (i = 0; i < height; ++i) {
515 __m256i sum32 = _mm256_setzero_si256();
516 for (j = 0; j <= width - 16; j += 16) {
517 const __m256i d0 = _mm256_cvtepu8_epi16(xx_loadu_128(dat + j));
518 const __m256i s0 = _mm256_cvtepu8_epi16(xx_loadu_128(src + j));
519 const __m256i flt0_16b = _mm256_permute4x64_epi64(
520 _mm256_packs_epi32(yy_loadu_256(flt0 + j),
521 yy_loadu_256(flt0 + j + 8)),
522 0xd8);
523 const __m256i flt1_16b = _mm256_permute4x64_epi64(
524 _mm256_packs_epi32(yy_loadu_256(flt1 + j),
525 yy_loadu_256(flt1 + j + 8)),
526 0xd8);
527 const __m256i u0 = _mm256_slli_epi16(d0, SGRPROJ_RST_BITS);
528 const __m256i flt0_0_sub_u = _mm256_sub_epi16(flt0_16b, u0);
529 const __m256i flt1_0_sub_u = _mm256_sub_epi16(flt1_16b, u0);
530 const __m256i v0 = _mm256_madd_epi16(
531 xq_coeff, _mm256_unpacklo_epi16(flt0_0_sub_u, flt1_0_sub_u));
532 const __m256i v1 = _mm256_madd_epi16(
533 xq_coeff, _mm256_unpackhi_epi16(flt0_0_sub_u, flt1_0_sub_u));
534 const __m256i vr0 =
535 _mm256_srai_epi32(_mm256_add_epi32(v0, rounding), shift);
536 const __m256i vr1 =
537 _mm256_srai_epi32(_mm256_add_epi32(v1, rounding), shift);
538 const __m256i e0 = _mm256_sub_epi16(
539 _mm256_add_epi16(_mm256_packs_epi32(vr0, vr1), d0), s0);
540 const __m256i err0 = _mm256_madd_epi16(e0, e0);
541 sum32 = _mm256_add_epi32(sum32, err0);
542 }
543 for (k = j; k < width; ++k) {
544 const int32_t u = (int32_t)(dat[k] << SGRPROJ_RST_BITS);
545 int32_t v = xq[0] * (flt0[k] - u) + xq[1] * (flt1[k] - u);
546 const int32_t e = ROUND_POWER_OF_TWO(v, shift) + dat[k] - src[k];
547 err += ((int64_t)e * e);
548 }
549 dat += dat_stride;
550 src += src_stride;
551 flt0 += flt0_stride;
552 flt1 += flt1_stride;
553 const __m256i sum64_0 =
554 _mm256_cvtepi32_epi64(_mm256_castsi256_si128(sum32));
555 const __m256i sum64_1 =
556 _mm256_cvtepi32_epi64(_mm256_extracti128_si256(sum32, 1));
557 sum64 = _mm256_add_epi64(sum64, sum64_0);
558 sum64 = _mm256_add_epi64(sum64, sum64_1);
559 }
560 } else if (params->r[0] > 0 || params->r[1] > 0) {
561 const int xq_active = (params->r[0] > 0) ? xq[0] : xq[1];
562 const __m256i xq_coeff =
563 pair_set_epi16(xq_active, (-xq_active * (1 << SGRPROJ_RST_BITS)));
564 const int32_t *flt = (params->r[0] > 0) ? flt0 : flt1;
565 const int flt_stride = (params->r[0] > 0) ? flt0_stride : flt1_stride;
566 for (i = 0; i < height; ++i) {
567 __m256i sum32 = _mm256_setzero_si256();
568 for (j = 0; j <= width - 16; j += 16) {
569 const __m256i d0 = _mm256_cvtepu8_epi16(xx_loadu_128(dat + j));
570 const __m256i s0 = _mm256_cvtepu8_epi16(xx_loadu_128(src + j));
571 const __m256i flt_16b = _mm256_permute4x64_epi64(
572 _mm256_packs_epi32(yy_loadu_256(flt + j),
573 yy_loadu_256(flt + j + 8)),
574 0xd8);
575 const __m256i v0 =
576 _mm256_madd_epi16(xq_coeff, _mm256_unpacklo_epi16(flt_16b, d0));
577 const __m256i v1 =
578 _mm256_madd_epi16(xq_coeff, _mm256_unpackhi_epi16(flt_16b, d0));
579 const __m256i vr0 =
580 _mm256_srai_epi32(_mm256_add_epi32(v0, rounding), shift);
581 const __m256i vr1 =
582 _mm256_srai_epi32(_mm256_add_epi32(v1, rounding), shift);
583 const __m256i e0 = _mm256_sub_epi16(
584 _mm256_add_epi16(_mm256_packs_epi32(vr0, vr1), d0), s0);
585 const __m256i err0 = _mm256_madd_epi16(e0, e0);
586 sum32 = _mm256_add_epi32(sum32, err0);
587 }
588 for (k = j; k < width; ++k) {
589 const int32_t u = (int32_t)(dat[k] << SGRPROJ_RST_BITS);
590 int32_t v = xq_active * (flt[k] - u);
591 const int32_t e = ROUND_POWER_OF_TWO(v, shift) + dat[k] - src[k];
592 err += ((int64_t)e * e);
593 }
594 dat += dat_stride;
595 src += src_stride;
596 flt += flt_stride;
597 const __m256i sum64_0 =
598 _mm256_cvtepi32_epi64(_mm256_castsi256_si128(sum32));
599 const __m256i sum64_1 =
600 _mm256_cvtepi32_epi64(_mm256_extracti128_si256(sum32, 1));
601 sum64 = _mm256_add_epi64(sum64, sum64_0);
602 sum64 = _mm256_add_epi64(sum64, sum64_1);
603 }
604 } else {
605 __m256i sum32 = _mm256_setzero_si256();
606 for (i = 0; i < height; ++i) {
607 for (j = 0; j <= width - 16; j += 16) {
608 const __m256i d0 = _mm256_cvtepu8_epi16(xx_loadu_128(dat + j));
609 const __m256i s0 = _mm256_cvtepu8_epi16(xx_loadu_128(src + j));
610 const __m256i diff0 = _mm256_sub_epi16(d0, s0);
611 const __m256i err0 = _mm256_madd_epi16(diff0, diff0);
612 sum32 = _mm256_add_epi32(sum32, err0);
613 }
614 for (k = j; k < width; ++k) {
615 const int32_t e = (int32_t)(dat[k]) - src[k];
616 err += ((int64_t)e * e);
617 }
618 dat += dat_stride;
619 src += src_stride;
620 }
621 const __m256i sum64_0 =
622 _mm256_cvtepi32_epi64(_mm256_castsi256_si128(sum32));
623 const __m256i sum64_1 =
624 _mm256_cvtepi32_epi64(_mm256_extracti128_si256(sum32, 1));
625 sum64 = _mm256_add_epi64(sum64_0, sum64_1);
626 }
627 int64_t sum[4];
628 yy_storeu_256(sum, sum64);
629 err += sum[0] + sum[1] + sum[2] + sum[3];
630 return err;
631 }
632
633 // When params->r[0] > 0 and params->r[1] > 0. In this case all elements of
634 // C and H need to be computed.
calc_proj_params_r0_r1_avx2(const uint8_t * src8,int width,int height,int src_stride,const uint8_t * dat8,int dat_stride,int32_t * flt0,int flt0_stride,int32_t * flt1,int flt1_stride,int64_t H[2][2],int64_t C[2])635 static AOM_INLINE void calc_proj_params_r0_r1_avx2(
636 const uint8_t *src8, int width, int height, int src_stride,
637 const uint8_t *dat8, int dat_stride, int32_t *flt0, int flt0_stride,
638 int32_t *flt1, int flt1_stride, int64_t H[2][2], int64_t C[2]) {
639 const int size = width * height;
640 const uint8_t *src = src8;
641 const uint8_t *dat = dat8;
642 __m256i h00, h01, h11, c0, c1;
643 const __m256i zero = _mm256_setzero_si256();
644 h01 = h11 = c0 = c1 = h00 = zero;
645
646 for (int i = 0; i < height; ++i) {
647 for (int j = 0; j < width; j += 8) {
648 const __m256i u_load = _mm256_cvtepu8_epi32(
649 _mm_loadl_epi64((__m128i *)(dat + i * dat_stride + j)));
650 const __m256i s_load = _mm256_cvtepu8_epi32(
651 _mm_loadl_epi64((__m128i *)(src + i * src_stride + j)));
652 __m256i f1 = _mm256_loadu_si256((__m256i *)(flt0 + i * flt0_stride + j));
653 __m256i f2 = _mm256_loadu_si256((__m256i *)(flt1 + i * flt1_stride + j));
654 __m256i d = _mm256_slli_epi32(u_load, SGRPROJ_RST_BITS);
655 __m256i s = _mm256_slli_epi32(s_load, SGRPROJ_RST_BITS);
656 s = _mm256_sub_epi32(s, d);
657 f1 = _mm256_sub_epi32(f1, d);
658 f2 = _mm256_sub_epi32(f2, d);
659
660 const __m256i h00_even = _mm256_mul_epi32(f1, f1);
661 const __m256i h00_odd = _mm256_mul_epi32(_mm256_srli_epi64(f1, 32),
662 _mm256_srli_epi64(f1, 32));
663 h00 = _mm256_add_epi64(h00, h00_even);
664 h00 = _mm256_add_epi64(h00, h00_odd);
665
666 const __m256i h01_even = _mm256_mul_epi32(f1, f2);
667 const __m256i h01_odd = _mm256_mul_epi32(_mm256_srli_epi64(f1, 32),
668 _mm256_srli_epi64(f2, 32));
669 h01 = _mm256_add_epi64(h01, h01_even);
670 h01 = _mm256_add_epi64(h01, h01_odd);
671
672 const __m256i h11_even = _mm256_mul_epi32(f2, f2);
673 const __m256i h11_odd = _mm256_mul_epi32(_mm256_srli_epi64(f2, 32),
674 _mm256_srli_epi64(f2, 32));
675 h11 = _mm256_add_epi64(h11, h11_even);
676 h11 = _mm256_add_epi64(h11, h11_odd);
677
678 const __m256i c0_even = _mm256_mul_epi32(f1, s);
679 const __m256i c0_odd =
680 _mm256_mul_epi32(_mm256_srli_epi64(f1, 32), _mm256_srli_epi64(s, 32));
681 c0 = _mm256_add_epi64(c0, c0_even);
682 c0 = _mm256_add_epi64(c0, c0_odd);
683
684 const __m256i c1_even = _mm256_mul_epi32(f2, s);
685 const __m256i c1_odd =
686 _mm256_mul_epi32(_mm256_srli_epi64(f2, 32), _mm256_srli_epi64(s, 32));
687 c1 = _mm256_add_epi64(c1, c1_even);
688 c1 = _mm256_add_epi64(c1, c1_odd);
689 }
690 }
691
692 __m256i c_low = _mm256_unpacklo_epi64(c0, c1);
693 const __m256i c_high = _mm256_unpackhi_epi64(c0, c1);
694 c_low = _mm256_add_epi64(c_low, c_high);
695 const __m128i c_128bit = _mm_add_epi64(_mm256_extracti128_si256(c_low, 1),
696 _mm256_castsi256_si128(c_low));
697
698 __m256i h0x_low = _mm256_unpacklo_epi64(h00, h01);
699 const __m256i h0x_high = _mm256_unpackhi_epi64(h00, h01);
700 h0x_low = _mm256_add_epi64(h0x_low, h0x_high);
701 const __m128i h0x_128bit = _mm_add_epi64(_mm256_extracti128_si256(h0x_low, 1),
702 _mm256_castsi256_si128(h0x_low));
703
704 // Using the symmetric properties of H, calculations of H[1][0] are not
705 // needed.
706 __m256i h1x_low = _mm256_unpacklo_epi64(zero, h11);
707 const __m256i h1x_high = _mm256_unpackhi_epi64(zero, h11);
708 h1x_low = _mm256_add_epi64(h1x_low, h1x_high);
709 const __m128i h1x_128bit = _mm_add_epi64(_mm256_extracti128_si256(h1x_low, 1),
710 _mm256_castsi256_si128(h1x_low));
711
712 xx_storeu_128(C, c_128bit);
713 xx_storeu_128(H[0], h0x_128bit);
714 xx_storeu_128(H[1], h1x_128bit);
715
716 H[0][0] /= size;
717 H[0][1] /= size;
718 H[1][1] /= size;
719
720 // Since H is a symmetric matrix
721 H[1][0] = H[0][1];
722 C[0] /= size;
723 C[1] /= size;
724 }
725
726 // When only params->r[0] > 0. In this case only H[0][0] and C[0] are
727 // non-zero and need to be computed.
calc_proj_params_r0_avx2(const uint8_t * src8,int width,int height,int src_stride,const uint8_t * dat8,int dat_stride,int32_t * flt0,int flt0_stride,int64_t H[2][2],int64_t C[2])728 static AOM_INLINE void calc_proj_params_r0_avx2(const uint8_t *src8, int width,
729 int height, int src_stride,
730 const uint8_t *dat8,
731 int dat_stride, int32_t *flt0,
732 int flt0_stride,
733 int64_t H[2][2], int64_t C[2]) {
734 const int size = width * height;
735 const uint8_t *src = src8;
736 const uint8_t *dat = dat8;
737 __m256i h00, c0;
738 const __m256i zero = _mm256_setzero_si256();
739 c0 = h00 = zero;
740
741 for (int i = 0; i < height; ++i) {
742 for (int j = 0; j < width; j += 8) {
743 const __m256i u_load = _mm256_cvtepu8_epi32(
744 _mm_loadl_epi64((__m128i *)(dat + i * dat_stride + j)));
745 const __m256i s_load = _mm256_cvtepu8_epi32(
746 _mm_loadl_epi64((__m128i *)(src + i * src_stride + j)));
747 __m256i f1 = _mm256_loadu_si256((__m256i *)(flt0 + i * flt0_stride + j));
748 __m256i d = _mm256_slli_epi32(u_load, SGRPROJ_RST_BITS);
749 __m256i s = _mm256_slli_epi32(s_load, SGRPROJ_RST_BITS);
750 s = _mm256_sub_epi32(s, d);
751 f1 = _mm256_sub_epi32(f1, d);
752
753 const __m256i h00_even = _mm256_mul_epi32(f1, f1);
754 const __m256i h00_odd = _mm256_mul_epi32(_mm256_srli_epi64(f1, 32),
755 _mm256_srli_epi64(f1, 32));
756 h00 = _mm256_add_epi64(h00, h00_even);
757 h00 = _mm256_add_epi64(h00, h00_odd);
758
759 const __m256i c0_even = _mm256_mul_epi32(f1, s);
760 const __m256i c0_odd =
761 _mm256_mul_epi32(_mm256_srli_epi64(f1, 32), _mm256_srli_epi64(s, 32));
762 c0 = _mm256_add_epi64(c0, c0_even);
763 c0 = _mm256_add_epi64(c0, c0_odd);
764 }
765 }
766 const __m128i h00_128bit = _mm_add_epi64(_mm256_extracti128_si256(h00, 1),
767 _mm256_castsi256_si128(h00));
768 const __m128i h00_val =
769 _mm_add_epi64(h00_128bit, _mm_srli_si128(h00_128bit, 8));
770
771 const __m128i c0_128bit = _mm_add_epi64(_mm256_extracti128_si256(c0, 1),
772 _mm256_castsi256_si128(c0));
773 const __m128i c0_val = _mm_add_epi64(c0_128bit, _mm_srli_si128(c0_128bit, 8));
774
775 const __m128i c = _mm_unpacklo_epi64(c0_val, _mm256_castsi256_si128(zero));
776 const __m128i h0x = _mm_unpacklo_epi64(h00_val, _mm256_castsi256_si128(zero));
777
778 xx_storeu_128(C, c);
779 xx_storeu_128(H[0], h0x);
780
781 H[0][0] /= size;
782 C[0] /= size;
783 }
784
785 // When only params->r[1] > 0. In this case only H[1][1] and C[1] are
786 // non-zero and need to be computed.
calc_proj_params_r1_avx2(const uint8_t * src8,int width,int height,int src_stride,const uint8_t * dat8,int dat_stride,int32_t * flt1,int flt1_stride,int64_t H[2][2],int64_t C[2])787 static AOM_INLINE void calc_proj_params_r1_avx2(const uint8_t *src8, int width,
788 int height, int src_stride,
789 const uint8_t *dat8,
790 int dat_stride, int32_t *flt1,
791 int flt1_stride,
792 int64_t H[2][2], int64_t C[2]) {
793 const int size = width * height;
794 const uint8_t *src = src8;
795 const uint8_t *dat = dat8;
796 __m256i h11, c1;
797 const __m256i zero = _mm256_setzero_si256();
798 c1 = h11 = zero;
799
800 for (int i = 0; i < height; ++i) {
801 for (int j = 0; j < width; j += 8) {
802 const __m256i u_load = _mm256_cvtepu8_epi32(
803 _mm_loadl_epi64((__m128i *)(dat + i * dat_stride + j)));
804 const __m256i s_load = _mm256_cvtepu8_epi32(
805 _mm_loadl_epi64((__m128i *)(src + i * src_stride + j)));
806 __m256i f2 = _mm256_loadu_si256((__m256i *)(flt1 + i * flt1_stride + j));
807 __m256i d = _mm256_slli_epi32(u_load, SGRPROJ_RST_BITS);
808 __m256i s = _mm256_slli_epi32(s_load, SGRPROJ_RST_BITS);
809 s = _mm256_sub_epi32(s, d);
810 f2 = _mm256_sub_epi32(f2, d);
811
812 const __m256i h11_even = _mm256_mul_epi32(f2, f2);
813 const __m256i h11_odd = _mm256_mul_epi32(_mm256_srli_epi64(f2, 32),
814 _mm256_srli_epi64(f2, 32));
815 h11 = _mm256_add_epi64(h11, h11_even);
816 h11 = _mm256_add_epi64(h11, h11_odd);
817
818 const __m256i c1_even = _mm256_mul_epi32(f2, s);
819 const __m256i c1_odd =
820 _mm256_mul_epi32(_mm256_srli_epi64(f2, 32), _mm256_srli_epi64(s, 32));
821 c1 = _mm256_add_epi64(c1, c1_even);
822 c1 = _mm256_add_epi64(c1, c1_odd);
823 }
824 }
825
826 const __m128i h11_128bit = _mm_add_epi64(_mm256_extracti128_si256(h11, 1),
827 _mm256_castsi256_si128(h11));
828 const __m128i h11_val =
829 _mm_add_epi64(h11_128bit, _mm_srli_si128(h11_128bit, 8));
830
831 const __m128i c1_128bit = _mm_add_epi64(_mm256_extracti128_si256(c1, 1),
832 _mm256_castsi256_si128(c1));
833 const __m128i c1_val = _mm_add_epi64(c1_128bit, _mm_srli_si128(c1_128bit, 8));
834
835 const __m128i c = _mm_unpacklo_epi64(_mm256_castsi256_si128(zero), c1_val);
836 const __m128i h1x = _mm_unpacklo_epi64(_mm256_castsi256_si128(zero), h11_val);
837
838 xx_storeu_128(C, c);
839 xx_storeu_128(H[1], h1x);
840
841 H[1][1] /= size;
842 C[1] /= size;
843 }
844
845 // AVX2 variant of av1_calc_proj_params_c.
av1_calc_proj_params_avx2(const uint8_t * src8,int width,int height,int src_stride,const uint8_t * dat8,int dat_stride,int32_t * flt0,int flt0_stride,int32_t * flt1,int flt1_stride,int64_t H[2][2],int64_t C[2],const sgr_params_type * params)846 void av1_calc_proj_params_avx2(const uint8_t *src8, int width, int height,
847 int src_stride, const uint8_t *dat8,
848 int dat_stride, int32_t *flt0, int flt0_stride,
849 int32_t *flt1, int flt1_stride, int64_t H[2][2],
850 int64_t C[2], const sgr_params_type *params) {
851 if ((params->r[0] > 0) && (params->r[1] > 0)) {
852 calc_proj_params_r0_r1_avx2(src8, width, height, src_stride, dat8,
853 dat_stride, flt0, flt0_stride, flt1,
854 flt1_stride, H, C);
855 } else if (params->r[0] > 0) {
856 calc_proj_params_r0_avx2(src8, width, height, src_stride, dat8, dat_stride,
857 flt0, flt0_stride, H, C);
858 } else if (params->r[1] > 0) {
859 calc_proj_params_r1_avx2(src8, width, height, src_stride, dat8, dat_stride,
860 flt1, flt1_stride, H, C);
861 }
862 }
863
864 #if CONFIG_AV1_HIGHBITDEPTH
av1_highbd_pixel_proj_error_avx2(const uint8_t * src8,int width,int height,int src_stride,const uint8_t * dat8,int dat_stride,int32_t * flt0,int flt0_stride,int32_t * flt1,int flt1_stride,int xq[2],const sgr_params_type * params)865 int64_t av1_highbd_pixel_proj_error_avx2(
866 const uint8_t *src8, int width, int height, int src_stride,
867 const uint8_t *dat8, int dat_stride, int32_t *flt0, int flt0_stride,
868 int32_t *flt1, int flt1_stride, int xq[2], const sgr_params_type *params) {
869 int i, j, k;
870 const int32_t shift = SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS;
871 const __m256i rounding = _mm256_set1_epi32(1 << (shift - 1));
872 __m256i sum64 = _mm256_setzero_si256();
873 const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
874 const uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
875 int64_t err = 0;
876 if (params->r[0] > 0 && params->r[1] > 0) { // Both filters are enabled
877 const __m256i xq0 = _mm256_set1_epi32(xq[0]);
878 const __m256i xq1 = _mm256_set1_epi32(xq[1]);
879 for (i = 0; i < height; ++i) {
880 __m256i sum32 = _mm256_setzero_si256();
881 for (j = 0; j <= width - 16; j += 16) { // Process 16 pixels at a time
882 // Load 16 pixels each from source image and corrupted image
883 const __m256i s0 = yy_loadu_256(src + j);
884 const __m256i d0 = yy_loadu_256(dat + j);
885 // s0 = [15 14 13 12 11 10 9 8] [7 6 5 4 3 2 1 0] as u16 (indices)
886
887 // Shift-up each pixel to match filtered image scaling
888 const __m256i u0 = _mm256_slli_epi16(d0, SGRPROJ_RST_BITS);
889
890 // Split u0 into two halves and pad each from u16 to i32
891 const __m256i u0l = _mm256_cvtepu16_epi32(_mm256_castsi256_si128(u0));
892 const __m256i u0h =
893 _mm256_cvtepu16_epi32(_mm256_extracti128_si256(u0, 1));
894 // u0h, u0l = [15 14 13 12] [11 10 9 8], [7 6 5 4] [3 2 1 0] as u32
895
896 // Load 16 pixels from each filtered image
897 const __m256i flt0l = yy_loadu_256(flt0 + j);
898 const __m256i flt0h = yy_loadu_256(flt0 + j + 8);
899 const __m256i flt1l = yy_loadu_256(flt1 + j);
900 const __m256i flt1h = yy_loadu_256(flt1 + j + 8);
901 // flt?l, flt?h = [15 14 13 12] [11 10 9 8], [7 6 5 4] [3 2 1 0] as u32
902
903 // Subtract shifted corrupt image from each filtered image
904 const __m256i flt0l_subu = _mm256_sub_epi32(flt0l, u0l);
905 const __m256i flt0h_subu = _mm256_sub_epi32(flt0h, u0h);
906 const __m256i flt1l_subu = _mm256_sub_epi32(flt1l, u0l);
907 const __m256i flt1h_subu = _mm256_sub_epi32(flt1h, u0h);
908
909 // Multiply basis vectors by appropriate coefficients
910 const __m256i v0l = _mm256_mullo_epi32(flt0l_subu, xq0);
911 const __m256i v0h = _mm256_mullo_epi32(flt0h_subu, xq0);
912 const __m256i v1l = _mm256_mullo_epi32(flt1l_subu, xq1);
913 const __m256i v1h = _mm256_mullo_epi32(flt1h_subu, xq1);
914
915 // Add together the contributions from the two basis vectors
916 const __m256i vl = _mm256_add_epi32(v0l, v1l);
917 const __m256i vh = _mm256_add_epi32(v0h, v1h);
918
919 // Right-shift v with appropriate rounding
920 const __m256i vrl =
921 _mm256_srai_epi32(_mm256_add_epi32(vl, rounding), shift);
922 const __m256i vrh =
923 _mm256_srai_epi32(_mm256_add_epi32(vh, rounding), shift);
924 // vrh, vrl = [15 14 13 12] [11 10 9 8], [7 6 5 4] [3 2 1 0]
925
926 // Saturate each i32 to an i16 then combine both halves
927 // The permute (control=[3 1 2 0]) fixes weird ordering from AVX lanes
928 const __m256i vr =
929 _mm256_permute4x64_epi64(_mm256_packs_epi32(vrl, vrh), 0xd8);
930 // intermediate = [15 14 13 12 7 6 5 4] [11 10 9 8 3 2 1 0]
931 // vr = [15 14 13 12 11 10 9 8] [7 6 5 4 3 2 1 0]
932
933 // Add twin-subspace-sgr-filter to corrupt image then subtract source
934 const __m256i e0 = _mm256_sub_epi16(_mm256_add_epi16(vr, d0), s0);
935
936 // Calculate squared error and add adjacent values
937 const __m256i err0 = _mm256_madd_epi16(e0, e0);
938
939 sum32 = _mm256_add_epi32(sum32, err0);
940 }
941
942 const __m256i sum32l =
943 _mm256_cvtepu32_epi64(_mm256_castsi256_si128(sum32));
944 sum64 = _mm256_add_epi64(sum64, sum32l);
945 const __m256i sum32h =
946 _mm256_cvtepu32_epi64(_mm256_extracti128_si256(sum32, 1));
947 sum64 = _mm256_add_epi64(sum64, sum32h);
948
949 // Process remaining pixels in this row (modulo 16)
950 for (k = j; k < width; ++k) {
951 const int32_t u = (int32_t)(dat[k] << SGRPROJ_RST_BITS);
952 int32_t v = xq[0] * (flt0[k] - u) + xq[1] * (flt1[k] - u);
953 const int32_t e = ROUND_POWER_OF_TWO(v, shift) + dat[k] - src[k];
954 err += ((int64_t)e * e);
955 }
956 dat += dat_stride;
957 src += src_stride;
958 flt0 += flt0_stride;
959 flt1 += flt1_stride;
960 }
961 } else if (params->r[0] > 0 || params->r[1] > 0) { // Only one filter enabled
962 const int32_t xq_on = (params->r[0] > 0) ? xq[0] : xq[1];
963 const __m256i xq_active = _mm256_set1_epi32(xq_on);
964 const __m256i xq_inactive =
965 _mm256_set1_epi32(-xq_on * (1 << SGRPROJ_RST_BITS));
966 const int32_t *flt = (params->r[0] > 0) ? flt0 : flt1;
967 const int flt_stride = (params->r[0] > 0) ? flt0_stride : flt1_stride;
968 for (i = 0; i < height; ++i) {
969 __m256i sum32 = _mm256_setzero_si256();
970 for (j = 0; j <= width - 16; j += 16) {
971 // Load 16 pixels from source image
972 const __m256i s0 = yy_loadu_256(src + j);
973 // s0 = [15 14 13 12 11 10 9 8] [7 6 5 4 3 2 1 0] as u16
974
975 // Load 16 pixels from corrupted image and pad each u16 to i32
976 const __m256i d0 = yy_loadu_256(dat + j);
977 const __m256i d0h =
978 _mm256_cvtepu16_epi32(_mm256_extracti128_si256(d0, 1));
979 const __m256i d0l = _mm256_cvtepu16_epi32(_mm256_castsi256_si128(d0));
980 // d0 = [15 14 13 12 11 10 9 8] [7 6 5 4 3 2 1 0] as u16
981 // d0h, d0l = [15 14 13 12] [11 10 9 8], [7 6 5 4] [3 2 1 0] as i32
982
983 // Load 16 pixels from the filtered image
984 const __m256i flth = yy_loadu_256(flt + j + 8);
985 const __m256i fltl = yy_loadu_256(flt + j);
986 // flth, fltl = [15 14 13 12] [11 10 9 8], [7 6 5 4] [3 2 1 0] as i32
987
988 const __m256i flth_xq = _mm256_mullo_epi32(flth, xq_active);
989 const __m256i fltl_xq = _mm256_mullo_epi32(fltl, xq_active);
990 const __m256i d0h_xq = _mm256_mullo_epi32(d0h, xq_inactive);
991 const __m256i d0l_xq = _mm256_mullo_epi32(d0l, xq_inactive);
992
993 const __m256i vh = _mm256_add_epi32(flth_xq, d0h_xq);
994 const __m256i vl = _mm256_add_epi32(fltl_xq, d0l_xq);
995
996 // Shift this down with appropriate rounding
997 const __m256i vrh =
998 _mm256_srai_epi32(_mm256_add_epi32(vh, rounding), shift);
999 const __m256i vrl =
1000 _mm256_srai_epi32(_mm256_add_epi32(vl, rounding), shift);
1001 // vrh, vrl = [15 14 13 12] [11 10 9 8], [7 6 5 4] [3 2 1 0] as i32
1002
1003 // Saturate each i32 to an i16 then combine both halves
1004 // The permute (control=[3 1 2 0]) fixes weird ordering from AVX lanes
1005 const __m256i vr =
1006 _mm256_permute4x64_epi64(_mm256_packs_epi32(vrl, vrh), 0xd8);
1007 // intermediate = [15 14 13 12 7 6 5 4] [11 10 9 8 3 2 1 0] as u16
1008 // vr = [15 14 13 12 11 10 9 8] [7 6 5 4 3 2 1 0] as u16
1009
1010 // Subtract twin-subspace-sgr filtered from source image to get error
1011 const __m256i e0 = _mm256_sub_epi16(_mm256_add_epi16(vr, d0), s0);
1012
1013 // Calculate squared error and add adjacent values
1014 const __m256i err0 = _mm256_madd_epi16(e0, e0);
1015
1016 sum32 = _mm256_add_epi32(sum32, err0);
1017 }
1018
1019 const __m256i sum32l =
1020 _mm256_cvtepu32_epi64(_mm256_castsi256_si128(sum32));
1021 sum64 = _mm256_add_epi64(sum64, sum32l);
1022 const __m256i sum32h =
1023 _mm256_cvtepu32_epi64(_mm256_extracti128_si256(sum32, 1));
1024 sum64 = _mm256_add_epi64(sum64, sum32h);
1025
1026 // Process remaining pixels in this row (modulo 16)
1027 for (k = j; k < width; ++k) {
1028 const int32_t u = (int32_t)(dat[k] << SGRPROJ_RST_BITS);
1029 int32_t v = xq_on * (flt[k] - u);
1030 const int32_t e = ROUND_POWER_OF_TWO(v, shift) + dat[k] - src[k];
1031 err += ((int64_t)e * e);
1032 }
1033 dat += dat_stride;
1034 src += src_stride;
1035 flt += flt_stride;
1036 }
1037 } else { // Neither filter is enabled
1038 for (i = 0; i < height; ++i) {
1039 __m256i sum32 = _mm256_setzero_si256();
1040 for (j = 0; j <= width - 32; j += 32) {
1041 // Load 2x16 u16 from source image
1042 const __m256i s0l = yy_loadu_256(src + j);
1043 const __m256i s0h = yy_loadu_256(src + j + 16);
1044
1045 // Load 2x16 u16 from corrupted image
1046 const __m256i d0l = yy_loadu_256(dat + j);
1047 const __m256i d0h = yy_loadu_256(dat + j + 16);
1048
1049 // Subtract corrupted image from source image
1050 const __m256i diffl = _mm256_sub_epi16(d0l, s0l);
1051 const __m256i diffh = _mm256_sub_epi16(d0h, s0h);
1052
1053 // Square error and add adjacent values
1054 const __m256i err0l = _mm256_madd_epi16(diffl, diffl);
1055 const __m256i err0h = _mm256_madd_epi16(diffh, diffh);
1056
1057 sum32 = _mm256_add_epi32(sum32, err0l);
1058 sum32 = _mm256_add_epi32(sum32, err0h);
1059 }
1060
1061 const __m256i sum32l =
1062 _mm256_cvtepu32_epi64(_mm256_castsi256_si128(sum32));
1063 sum64 = _mm256_add_epi64(sum64, sum32l);
1064 const __m256i sum32h =
1065 _mm256_cvtepu32_epi64(_mm256_extracti128_si256(sum32, 1));
1066 sum64 = _mm256_add_epi64(sum64, sum32h);
1067
1068 // Process remaining pixels (modulu 16)
1069 for (k = j; k < width; ++k) {
1070 const int32_t e = (int32_t)(dat[k]) - src[k];
1071 err += ((int64_t)e * e);
1072 }
1073 dat += dat_stride;
1074 src += src_stride;
1075 }
1076 }
1077
1078 // Sum 4 values from sum64l and sum64h into err
1079 int64_t sum[4];
1080 yy_storeu_256(sum, sum64);
1081 err += sum[0] + sum[1] + sum[2] + sum[3];
1082 return err;
1083 }
1084 #endif // CONFIG_AV1_HIGHBITDEPTH
1085