• Home
  • History
  • Annotate
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 <assert.h>
13 #include <emmintrin.h>
14 #include "aom_dsp/x86/synonyms.h"
15 
16 #include "config/av1_rtcd.h"
17 #include "av1/common/restoration.h"
18 #include "av1/encoder/pickrst.h"
19 
acc_stat_sse41(int32_t * dst,const uint8_t * src,const __m128i * shuffle,const __m128i * kl)20 static INLINE void acc_stat_sse41(int32_t *dst, const uint8_t *src,
21                                   const __m128i *shuffle, const __m128i *kl) {
22   const __m128i s = _mm_shuffle_epi8(xx_loadu_128(src), *shuffle);
23   const __m128i d0 = _mm_madd_epi16(*kl, _mm_cvtepu8_epi16(s));
24   const __m128i d1 =
25       _mm_madd_epi16(*kl, _mm_cvtepu8_epi16(_mm_srli_si128(s, 8)));
26   const __m128i dst0 = xx_loadu_128(dst);
27   const __m128i dst1 = xx_loadu_128(dst + 4);
28   const __m128i r0 = _mm_add_epi32(dst0, d0);
29   const __m128i r1 = _mm_add_epi32(dst1, d1);
30   xx_storeu_128(dst, r0);
31   xx_storeu_128(dst + 4, r1);
32 }
33 
acc_stat_win7_one_line_sse4_1(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])34 static INLINE void acc_stat_win7_one_line_sse4_1(
35     const uint8_t *dgd, const uint8_t *src, int h_start, int h_end,
36     int dgd_stride, const __m128i *shuffle, int32_t *sumX,
37     int32_t sumY[WIENER_WIN][WIENER_WIN], int32_t M_int[WIENER_WIN][WIENER_WIN],
38     int32_t H_int[WIENER_WIN2][WIENER_WIN * 8]) {
39   const int wiener_win = 7;
40   int j, k, l;
41   for (j = h_start; j < h_end; j += 2) {
42     const uint8_t *dgd_ij = dgd + j;
43     const uint8_t X1 = src[j];
44     const uint8_t X2 = src[j + 1];
45     *sumX += X1 + X2;
46     for (k = 0; k < wiener_win; k++) {
47       const uint8_t *dgd_ijk = dgd_ij + k * dgd_stride;
48       for (l = 0; l < wiener_win; l++) {
49         int32_t *H_ = &H_int[(l * wiener_win + k)][0];
50         const uint8_t D1 = dgd_ijk[l];
51         const uint8_t D2 = dgd_ijk[l + 1];
52         sumY[k][l] += D1 + D2;
53         M_int[k][l] += D1 * X1 + D2 * X2;
54 
55         const __m128i kl =
56             _mm_cvtepu8_epi16(_mm_set1_epi16(*((uint16_t *)(dgd_ijk + l))));
57         acc_stat_sse41(H_ + 0 * 8, dgd_ij + 0 * dgd_stride, shuffle, &kl);
58         acc_stat_sse41(H_ + 1 * 8, dgd_ij + 1 * dgd_stride, shuffle, &kl);
59         acc_stat_sse41(H_ + 2 * 8, dgd_ij + 2 * dgd_stride, shuffle, &kl);
60         acc_stat_sse41(H_ + 3 * 8, dgd_ij + 3 * dgd_stride, shuffle, &kl);
61         acc_stat_sse41(H_ + 4 * 8, dgd_ij + 4 * dgd_stride, shuffle, &kl);
62         acc_stat_sse41(H_ + 5 * 8, dgd_ij + 5 * dgd_stride, shuffle, &kl);
63         acc_stat_sse41(H_ + 6 * 8, dgd_ij + 6 * dgd_stride, shuffle, &kl);
64       }
65     }
66   }
67 }
68 
compute_stats_win7_opt_sse4_1(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)69 static INLINE void compute_stats_win7_opt_sse4_1(
70     const uint8_t *dgd, const uint8_t *src, int h_start, int h_end, int v_start,
71     int v_end, int dgd_stride, int src_stride, int64_t *M, int64_t *H) {
72   int i, j, k, l, m, n;
73   const int wiener_win = WIENER_WIN;
74   const int pixel_count = (h_end - h_start) * (v_end - v_start);
75   const int wiener_win2 = wiener_win * wiener_win;
76   const int wiener_halfwin = (wiener_win >> 1);
77   const uint8_t avg =
78       find_average(dgd, h_start, h_end, v_start, v_end, dgd_stride);
79 
80   int32_t M_int32[WIENER_WIN][WIENER_WIN] = { { 0 } };
81   int64_t M_int64[WIENER_WIN][WIENER_WIN] = { { 0 } };
82   int32_t H_int32[WIENER_WIN2][WIENER_WIN * 8] = { { 0 } };
83   int64_t H_int64[WIENER_WIN2][WIENER_WIN * 8] = { { 0 } };
84   int32_t sumY[WIENER_WIN][WIENER_WIN] = { { 0 } };
85   int32_t sumX = 0;
86   const uint8_t *dgd_win = dgd - wiener_halfwin * dgd_stride - wiener_halfwin;
87 
88   const __m128i shuffle = xx_loadu_128(g_shuffle_stats_data);
89   for (j = v_start; j < v_end; j += 64) {
90     const int vert_end = AOMMIN(64, v_end - j) + j;
91     for (i = j; i < vert_end; i++) {
92       acc_stat_win7_one_line_sse4_1(
93           dgd_win + i * dgd_stride, src + i * src_stride, h_start, h_end,
94           dgd_stride, &shuffle, &sumX, sumY, M_int32, H_int32);
95     }
96     for (k = 0; k < wiener_win; ++k) {
97       for (l = 0; l < wiener_win; ++l) {
98         M_int64[k][l] += M_int32[k][l];
99         M_int32[k][l] = 0;
100       }
101     }
102     for (k = 0; k < WIENER_WIN2; ++k) {
103       for (l = 0; l < WIENER_WIN * 8; ++l) {
104         H_int64[k][l] += H_int32[k][l];
105         H_int32[k][l] = 0;
106       }
107     }
108   }
109 
110   const int64_t avg_square_sum = (int64_t)avg * (int64_t)avg * pixel_count;
111   for (k = 0; k < wiener_win; k++) {
112     for (l = 0; l < wiener_win; l++) {
113       const int32_t idx0 = l * wiener_win + k;
114       M[idx0] =
115           M_int64[k][l] + (avg_square_sum - (int64_t)avg * (sumX + sumY[k][l]));
116       int64_t *H_ = H + idx0 * wiener_win2;
117       int64_t *H_int_ = &H_int64[idx0][0];
118       for (m = 0; m < wiener_win; m++) {
119         for (n = 0; n < wiener_win; n++) {
120           H_[m * wiener_win + n] = H_int_[n * 8 + m] + avg_square_sum -
121                                    (int64_t)avg * (sumY[k][l] + sumY[n][m]);
122         }
123       }
124     }
125   }
126 }
127 
acc_stat_highbd_sse41(int64_t * dst,const uint16_t * dgd,const __m128i * shuffle,const __m128i * dgd_ijkl)128 static INLINE void acc_stat_highbd_sse41(int64_t *dst, const uint16_t *dgd,
129                                          const __m128i *shuffle,
130                                          const __m128i *dgd_ijkl) {
131   // Load 256 bits from dgd in two chunks
132   const __m128i s0l = xx_loadu_128(dgd);
133   const __m128i s0h = xx_loadu_128(dgd + 4);
134   // s0l = [7 6 5 4 3 2 1 0] as u16 values (dgd indices)
135   // s0h = [11 10 9 8 7 6 5 4] as u16 values (dgd indices)
136   // (Slightly strange order so we can apply the same shuffle to both halves)
137 
138   // Shuffle the u16 values in each half (actually using 8-bit shuffle mask)
139   const __m128i s1l = _mm_shuffle_epi8(s0l, *shuffle);
140   const __m128i s1h = _mm_shuffle_epi8(s0h, *shuffle);
141   // s1l = [4 3 3 2 2 1 1 0] as u16 values (dgd indices)
142   // s1h = [8 7 7 6 6 5 5 4] as u16 values (dgd indices)
143 
144   // Multiply s1 by dgd_ijkl resulting in 8x u32 values
145   // Horizontally add pairs of u32 resulting in 4x u32
146   const __m128i dl = _mm_madd_epi16(*dgd_ijkl, s1l);
147   const __m128i dh = _mm_madd_epi16(*dgd_ijkl, s1h);
148   // dl = [d c b a] as u32 values
149   // dh = [h g f e] as u32 values
150 
151   // Add these 8x u32 results on to dst in four parts
152   const __m128i dll = _mm_cvtepu32_epi64(dl);
153   const __m128i dlh = _mm_cvtepu32_epi64(_mm_srli_si128(dl, 8));
154   const __m128i dhl = _mm_cvtepu32_epi64(dh);
155   const __m128i dhh = _mm_cvtepu32_epi64(_mm_srli_si128(dh, 8));
156   // dll = [b a] as u64 values, etc.
157 
158   const __m128i rll = _mm_add_epi64(xx_loadu_128(dst), dll);
159   xx_storeu_128(dst, rll);
160   const __m128i rlh = _mm_add_epi64(xx_loadu_128(dst + 2), dlh);
161   xx_storeu_128(dst + 2, rlh);
162   const __m128i rhl = _mm_add_epi64(xx_loadu_128(dst + 4), dhl);
163   xx_storeu_128(dst + 4, rhl);
164   const __m128i rhh = _mm_add_epi64(xx_loadu_128(dst + 6), dhh);
165   xx_storeu_128(dst + 6, rhh);
166 }
167 
acc_stat_highbd_win7_one_line_sse4_1(const uint16_t * dgd,const uint16_t * src,int h_start,int h_end,int dgd_stride,const __m128i * 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])168 static INLINE void acc_stat_highbd_win7_one_line_sse4_1(
169     const uint16_t *dgd, const uint16_t *src, int h_start, int h_end,
170     int dgd_stride, const __m128i *shuffle, int32_t *sumX,
171     int32_t sumY[WIENER_WIN][WIENER_WIN], int64_t M_int[WIENER_WIN][WIENER_WIN],
172     int64_t H_int[WIENER_WIN2][WIENER_WIN * 8]) {
173   int j, k, l;
174   const int wiener_win = WIENER_WIN;
175   for (j = h_start; j < h_end; j += 2) {
176     const uint16_t X1 = src[j];
177     const uint16_t X2 = src[j + 1];
178     *sumX += X1 + X2;
179     const uint16_t *dgd_ij = dgd + j;
180     for (k = 0; k < wiener_win; k++) {
181       const uint16_t *dgd_ijk = dgd_ij + k * dgd_stride;
182       for (l = 0; l < wiener_win; l++) {
183         int64_t *H_ = &H_int[(l * wiener_win + k)][0];
184         const uint16_t D1 = dgd_ijk[l];
185         const uint16_t D2 = dgd_ijk[l + 1];
186         sumY[k][l] += D1 + D2;
187         M_int[k][l] += D1 * X1 + D2 * X2;
188 
189         // Load two u16 values from dgd as a single u32
190         // Then broadcast to 4x u32 slots of a 128
191         const __m128i dgd_ijkl = _mm_set1_epi32(*((uint32_t *)(dgd_ijk + l)));
192         // dgd_ijkl = [y x y x y x y x] as u16
193 
194         acc_stat_highbd_sse41(H_ + 0 * 8, dgd_ij + 0 * dgd_stride, shuffle,
195                               &dgd_ijkl);
196         acc_stat_highbd_sse41(H_ + 1 * 8, dgd_ij + 1 * dgd_stride, shuffle,
197                               &dgd_ijkl);
198         acc_stat_highbd_sse41(H_ + 2 * 8, dgd_ij + 2 * dgd_stride, shuffle,
199                               &dgd_ijkl);
200         acc_stat_highbd_sse41(H_ + 3 * 8, dgd_ij + 3 * dgd_stride, shuffle,
201                               &dgd_ijkl);
202         acc_stat_highbd_sse41(H_ + 4 * 8, dgd_ij + 4 * dgd_stride, shuffle,
203                               &dgd_ijkl);
204         acc_stat_highbd_sse41(H_ + 5 * 8, dgd_ij + 5 * dgd_stride, shuffle,
205                               &dgd_ijkl);
206         acc_stat_highbd_sse41(H_ + 6 * 8, dgd_ij + 6 * dgd_stride, shuffle,
207                               &dgd_ijkl);
208       }
209     }
210   }
211 }
212 
compute_stats_highbd_win7_opt_sse4_1(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)213 static INLINE void compute_stats_highbd_win7_opt_sse4_1(
214     const uint8_t *dgd8, const uint8_t *src8, int h_start, int h_end,
215     int v_start, int v_end, int dgd_stride, int src_stride, int64_t *M,
216     int64_t *H, aom_bit_depth_t bit_depth) {
217   int i, j, k, l, m, n;
218   const int wiener_win = WIENER_WIN;
219   const int pixel_count = (h_end - h_start) * (v_end - v_start);
220   const int wiener_win2 = wiener_win * wiener_win;
221   const int wiener_halfwin = (wiener_win >> 1);
222   const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
223   const uint16_t *dgd = CONVERT_TO_SHORTPTR(dgd8);
224   const uint16_t avg =
225       find_average_highbd(dgd, h_start, h_end, v_start, v_end, dgd_stride);
226 
227   int64_t M_int[WIENER_WIN][WIENER_WIN] = { { 0 } };
228   int64_t H_int[WIENER_WIN2][WIENER_WIN * 8] = { { 0 } };
229   int32_t sumY[WIENER_WIN][WIENER_WIN] = { { 0 } };
230   int32_t sumX = 0;
231   const uint16_t *dgd_win = dgd - wiener_halfwin * dgd_stride - wiener_halfwin;
232 
233   // Load just half of the 256-bit shuffle control used for the AVX2 version
234   const __m128i shuffle = xx_loadu_128(g_shuffle_stats_highbd_data);
235   for (j = v_start; j < v_end; j += 64) {
236     const int vert_end = AOMMIN(64, v_end - j) + j;
237     for (i = j; i < vert_end; i++) {
238       acc_stat_highbd_win7_one_line_sse4_1(
239           dgd_win + i * dgd_stride, src + i * src_stride, h_start, h_end,
240           dgd_stride, &shuffle, &sumX, sumY, M_int, H_int);
241     }
242   }
243 
244   uint8_t bit_depth_divider = 1;
245   if (bit_depth == AOM_BITS_12)
246     bit_depth_divider = 16;
247   else if (bit_depth == AOM_BITS_10)
248     bit_depth_divider = 4;
249 
250   const int64_t avg_square_sum = (int64_t)avg * (int64_t)avg * pixel_count;
251   for (k = 0; k < wiener_win; k++) {
252     for (l = 0; l < wiener_win; l++) {
253       const int32_t idx0 = l * wiener_win + k;
254       M[idx0] = (M_int[k][l] +
255                  (avg_square_sum - (int64_t)avg * (sumX + sumY[k][l]))) /
256                 bit_depth_divider;
257       int64_t *H_ = H + idx0 * wiener_win2;
258       int64_t *H_int_ = &H_int[idx0][0];
259       for (m = 0; m < wiener_win; m++) {
260         for (n = 0; n < wiener_win; n++) {
261           H_[m * wiener_win + n] =
262               (H_int_[n * 8 + m] +
263                (avg_square_sum - (int64_t)avg * (sumY[k][l] + sumY[n][m]))) /
264               bit_depth_divider;
265         }
266       }
267     }
268   }
269 }
270 
acc_stat_highbd_win5_one_line_sse4_1(const uint16_t * dgd,const uint16_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],int64_t M_int[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],int64_t H_int[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8])271 static INLINE void acc_stat_highbd_win5_one_line_sse4_1(
272     const uint16_t *dgd, const uint16_t *src, int h_start, int h_end,
273     int dgd_stride, const __m128i *shuffle, int32_t *sumX,
274     int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],
275     int64_t M_int[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],
276     int64_t H_int[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8]) {
277   int j, k, l;
278   const int wiener_win = WIENER_WIN_CHROMA;
279   for (j = h_start; j < h_end; j += 2) {
280     const uint16_t X1 = src[j];
281     const uint16_t X2 = src[j + 1];
282     *sumX += X1 + X2;
283     const uint16_t *dgd_ij = dgd + j;
284     for (k = 0; k < wiener_win; k++) {
285       const uint16_t *dgd_ijk = dgd_ij + k * dgd_stride;
286       for (l = 0; l < wiener_win; l++) {
287         int64_t *H_ = &H_int[(l * wiener_win + k)][0];
288         const uint16_t D1 = dgd_ijk[l];
289         const uint16_t D2 = dgd_ijk[l + 1];
290         sumY[k][l] += D1 + D2;
291         M_int[k][l] += D1 * X1 + D2 * X2;
292 
293         // Load two u16 values from dgd as a single u32
294         // then broadcast to 4x u32 slots of a 128
295         const __m128i dgd_ijkl = _mm_set1_epi32(*((uint32_t *)(dgd_ijk + l)));
296         // dgd_ijkl = [y x y x y x y x] as u16
297 
298         acc_stat_highbd_sse41(H_ + 0 * 8, dgd_ij + 0 * dgd_stride, shuffle,
299                               &dgd_ijkl);
300         acc_stat_highbd_sse41(H_ + 1 * 8, dgd_ij + 1 * dgd_stride, shuffle,
301                               &dgd_ijkl);
302         acc_stat_highbd_sse41(H_ + 2 * 8, dgd_ij + 2 * dgd_stride, shuffle,
303                               &dgd_ijkl);
304         acc_stat_highbd_sse41(H_ + 3 * 8, dgd_ij + 3 * dgd_stride, shuffle,
305                               &dgd_ijkl);
306         acc_stat_highbd_sse41(H_ + 4 * 8, dgd_ij + 4 * dgd_stride, shuffle,
307                               &dgd_ijkl);
308       }
309     }
310   }
311 }
312 
compute_stats_highbd_win5_opt_sse4_1(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)313 static INLINE void compute_stats_highbd_win5_opt_sse4_1(
314     const uint8_t *dgd8, const uint8_t *src8, int h_start, int h_end,
315     int v_start, int v_end, int dgd_stride, int src_stride, int64_t *M,
316     int64_t *H, aom_bit_depth_t bit_depth) {
317   int i, j, k, l, m, n;
318   const int wiener_win = WIENER_WIN_CHROMA;
319   const int pixel_count = (h_end - h_start) * (v_end - v_start);
320   const int wiener_win2 = wiener_win * wiener_win;
321   const int wiener_halfwin = (wiener_win >> 1);
322   const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
323   const uint16_t *dgd = CONVERT_TO_SHORTPTR(dgd8);
324   const uint16_t avg =
325       find_average_highbd(dgd, h_start, h_end, v_start, v_end, dgd_stride);
326 
327   int64_t M_int[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
328   int64_t H_int[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8] = { { 0 } };
329   int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
330   int32_t sumX = 0;
331   const uint16_t *dgd_win = dgd - wiener_halfwin * dgd_stride - wiener_halfwin;
332 
333   // Load just half of the 256-bit shuffle control used for the AVX2 version
334   const __m128i shuffle = xx_loadu_128(g_shuffle_stats_highbd_data);
335   for (j = v_start; j < v_end; j += 64) {
336     const int vert_end = AOMMIN(64, v_end - j) + j;
337     for (i = j; i < vert_end; i++) {
338       acc_stat_highbd_win5_one_line_sse4_1(
339           dgd_win + i * dgd_stride, src + i * src_stride, h_start, h_end,
340           dgd_stride, &shuffle, &sumX, sumY, M_int, H_int);
341     }
342   }
343 
344   uint8_t bit_depth_divider = 1;
345   if (bit_depth == AOM_BITS_12)
346     bit_depth_divider = 16;
347   else if (bit_depth == AOM_BITS_10)
348     bit_depth_divider = 4;
349 
350   const int64_t avg_square_sum = (int64_t)avg * (int64_t)avg * pixel_count;
351   for (k = 0; k < wiener_win; k++) {
352     for (l = 0; l < wiener_win; l++) {
353       const int32_t idx0 = l * wiener_win + k;
354       M[idx0] = (M_int[k][l] +
355                  (avg_square_sum - (int64_t)avg * (sumX + sumY[k][l]))) /
356                 bit_depth_divider;
357       int64_t *H_ = H + idx0 * wiener_win2;
358       int64_t *H_int_ = &H_int[idx0][0];
359       for (m = 0; m < wiener_win; m++) {
360         for (n = 0; n < wiener_win; n++) {
361           H_[m * wiener_win + n] =
362               (H_int_[n * 8 + m] +
363                (avg_square_sum - (int64_t)avg * (sumY[k][l] + sumY[n][m]))) /
364               bit_depth_divider;
365         }
366       }
367     }
368   }
369 }
370 
av1_compute_stats_highbd_sse4_1(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)371 void av1_compute_stats_highbd_sse4_1(int wiener_win, const uint8_t *dgd8,
372                                      const uint8_t *src8, int h_start,
373                                      int h_end, int v_start, int v_end,
374                                      int dgd_stride, int src_stride, int64_t *M,
375                                      int64_t *H, aom_bit_depth_t bit_depth) {
376   if (wiener_win == WIENER_WIN) {
377     compute_stats_highbd_win7_opt_sse4_1(dgd8, src8, h_start, h_end, v_start,
378                                          v_end, dgd_stride, src_stride, M, H,
379                                          bit_depth);
380   } else if (wiener_win == WIENER_WIN_CHROMA) {
381     compute_stats_highbd_win5_opt_sse4_1(dgd8, src8, h_start, h_end, v_start,
382                                          v_end, dgd_stride, src_stride, M, H,
383                                          bit_depth);
384   } else {
385     av1_compute_stats_highbd_c(wiener_win, dgd8, src8, h_start, h_end, v_start,
386                                v_end, dgd_stride, src_stride, M, H, bit_depth);
387   }
388 }
389 
acc_stat_win5_one_line_sse4_1(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])390 static INLINE void acc_stat_win5_one_line_sse4_1(
391     const uint8_t *dgd, const uint8_t *src, int h_start, int h_end,
392     int dgd_stride, const __m128i *shuffle, int32_t *sumX,
393     int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],
394     int32_t M_int[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA],
395     int32_t H_int[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8]) {
396   const int wiener_win = WIENER_WIN_CHROMA;
397   int j, k, l;
398   for (j = h_start; j < h_end; j += 2) {
399     const uint8_t *dgd_ij = dgd + j;
400     const uint8_t X1 = src[j];
401     const uint8_t X2 = src[j + 1];
402     *sumX += X1 + X2;
403     for (k = 0; k < wiener_win; k++) {
404       const uint8_t *dgd_ijk = dgd_ij + k * dgd_stride;
405       for (l = 0; l < wiener_win; l++) {
406         int32_t *H_ = &H_int[(l * wiener_win + k)][0];
407         const uint8_t D1 = dgd_ijk[l];
408         const uint8_t D2 = dgd_ijk[l + 1];
409         sumY[k][l] += D1 + D2;
410         M_int[k][l] += D1 * X1 + D2 * X2;
411 
412         const __m128i kl =
413             _mm_cvtepu8_epi16(_mm_set1_epi16(*((uint16_t *)(dgd_ijk + l))));
414         acc_stat_sse41(H_ + 0 * 8, dgd_ij + 0 * dgd_stride, shuffle, &kl);
415         acc_stat_sse41(H_ + 1 * 8, dgd_ij + 1 * dgd_stride, shuffle, &kl);
416         acc_stat_sse41(H_ + 2 * 8, dgd_ij + 2 * dgd_stride, shuffle, &kl);
417         acc_stat_sse41(H_ + 3 * 8, dgd_ij + 3 * dgd_stride, shuffle, &kl);
418         acc_stat_sse41(H_ + 4 * 8, dgd_ij + 4 * dgd_stride, shuffle, &kl);
419       }
420     }
421   }
422 }
423 
compute_stats_win5_opt_sse4_1(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)424 static INLINE void compute_stats_win5_opt_sse4_1(
425     const uint8_t *dgd, const uint8_t *src, int h_start, int h_end, int v_start,
426     int v_end, int dgd_stride, int src_stride, int64_t *M, int64_t *H) {
427   int i, j, k, l, m, n;
428   const int wiener_win = WIENER_WIN_CHROMA;
429   const int pixel_count = (h_end - h_start) * (v_end - v_start);
430   const int wiener_win2 = wiener_win * wiener_win;
431   const int wiener_halfwin = (wiener_win >> 1);
432   const uint8_t avg =
433       find_average(dgd, h_start, h_end, v_start, v_end, dgd_stride);
434 
435   int32_t M_int32[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
436   int64_t M_int64[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
437   int32_t H_int32[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8] = { { 0 } };
438   int64_t H_int64[WIENER_WIN2_CHROMA][WIENER_WIN_CHROMA * 8] = { { 0 } };
439   int32_t sumY[WIENER_WIN_CHROMA][WIENER_WIN_CHROMA] = { { 0 } };
440   int32_t sumX = 0;
441   const uint8_t *dgd_win = dgd - wiener_halfwin * dgd_stride - wiener_halfwin;
442 
443   const __m128i shuffle = xx_loadu_128(g_shuffle_stats_data);
444   for (j = v_start; j < v_end; j += 64) {
445     const int vert_end = AOMMIN(64, v_end - j) + j;
446     for (i = j; i < vert_end; i++) {
447       acc_stat_win5_one_line_sse4_1(
448           dgd_win + i * dgd_stride, src + i * src_stride, h_start, h_end,
449           dgd_stride, &shuffle, &sumX, sumY, M_int32, H_int32);
450     }
451     for (k = 0; k < wiener_win; ++k) {
452       for (l = 0; l < wiener_win; ++l) {
453         M_int64[k][l] += M_int32[k][l];
454         M_int32[k][l] = 0;
455       }
456     }
457     for (k = 0; k < WIENER_WIN_CHROMA * WIENER_WIN_CHROMA; ++k) {
458       for (l = 0; l < WIENER_WIN_CHROMA * 8; ++l) {
459         H_int64[k][l] += H_int32[k][l];
460         H_int32[k][l] = 0;
461       }
462     }
463   }
464 
465   const int64_t avg_square_sum = (int64_t)avg * (int64_t)avg * pixel_count;
466   for (k = 0; k < wiener_win; k++) {
467     for (l = 0; l < wiener_win; l++) {
468       const int32_t idx0 = l * wiener_win + k;
469       M[idx0] =
470           M_int64[k][l] + (avg_square_sum - (int64_t)avg * (sumX + sumY[k][l]));
471       int64_t *H_ = H + idx0 * wiener_win2;
472       int64_t *H_int_ = &H_int64[idx0][0];
473       for (m = 0; m < wiener_win; m++) {
474         for (n = 0; n < wiener_win; n++) {
475           H_[m * wiener_win + n] = H_int_[n * 8 + m] + avg_square_sum -
476                                    (int64_t)avg * (sumY[k][l] + sumY[n][m]);
477         }
478       }
479     }
480   }
481 }
av1_compute_stats_sse4_1(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)482 void av1_compute_stats_sse4_1(int wiener_win, const uint8_t *dgd,
483                               const uint8_t *src, int h_start, int h_end,
484                               int v_start, int v_end, int dgd_stride,
485                               int src_stride, int64_t *M, int64_t *H) {
486   if (wiener_win == WIENER_WIN) {
487     compute_stats_win7_opt_sse4_1(dgd, src, h_start, h_end, v_start, v_end,
488                                   dgd_stride, src_stride, M, H);
489   } else if (wiener_win == WIENER_WIN_CHROMA) {
490     compute_stats_win5_opt_sse4_1(dgd, src, h_start, h_end, v_start, v_end,
491                                   dgd_stride, src_stride, M, H);
492   } else {
493     av1_compute_stats_c(wiener_win, dgd, src, h_start, h_end, v_start, v_end,
494                         dgd_stride, src_stride, M, H);
495   }
496 }
497 
pair_set_epi16(uint16_t a,uint16_t b)498 static INLINE __m128i pair_set_epi16(uint16_t a, uint16_t b) {
499   return _mm_set1_epi32((int32_t)(((uint16_t)(a)) | (((uint32_t)(b)) << 16)));
500 }
501 
av1_lowbd_pixel_proj_error_sse4_1(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)502 int64_t av1_lowbd_pixel_proj_error_sse4_1(
503     const uint8_t *src8, int width, int height, int src_stride,
504     const uint8_t *dat8, int dat_stride, int32_t *flt0, int flt0_stride,
505     int32_t *flt1, int flt1_stride, int xq[2], const sgr_params_type *params) {
506   int i, j, k;
507   const int32_t shift = SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS;
508   const __m128i rounding = _mm_set1_epi32(1 << (shift - 1));
509   __m128i sum64 = _mm_setzero_si128();
510   const uint8_t *src = src8;
511   const uint8_t *dat = dat8;
512   int64_t err = 0;
513   if (params->r[0] > 0 && params->r[1] > 0) {
514     __m128i xq_coeff = pair_set_epi16(xq[0], xq[1]);
515     for (i = 0; i < height; ++i) {
516       __m128i sum32 = _mm_setzero_si128();
517       for (j = 0; j <= width - 8; j += 8) {
518         const __m128i d0 = _mm_cvtepu8_epi16(xx_loadl_64(dat + j));
519         const __m128i s0 = _mm_cvtepu8_epi16(xx_loadl_64(src + j));
520         const __m128i flt0_16b =
521             _mm_packs_epi32(xx_loadu_128(flt0 + j), xx_loadu_128(flt0 + j + 4));
522         const __m128i flt1_16b =
523             _mm_packs_epi32(xx_loadu_128(flt1 + j), xx_loadu_128(flt1 + j + 4));
524         const __m128i u0 = _mm_slli_epi16(d0, SGRPROJ_RST_BITS);
525         const __m128i flt0_0_sub_u = _mm_sub_epi16(flt0_16b, u0);
526         const __m128i flt1_0_sub_u = _mm_sub_epi16(flt1_16b, u0);
527         const __m128i v0 = _mm_madd_epi16(
528             xq_coeff, _mm_unpacklo_epi16(flt0_0_sub_u, flt1_0_sub_u));
529         const __m128i v1 = _mm_madd_epi16(
530             xq_coeff, _mm_unpackhi_epi16(flt0_0_sub_u, flt1_0_sub_u));
531         const __m128i vr0 = _mm_srai_epi32(_mm_add_epi32(v0, rounding), shift);
532         const __m128i vr1 = _mm_srai_epi32(_mm_add_epi32(v1, rounding), shift);
533         const __m128i e0 =
534             _mm_sub_epi16(_mm_add_epi16(_mm_packs_epi32(vr0, vr1), d0), s0);
535         const __m128i err0 = _mm_madd_epi16(e0, e0);
536         sum32 = _mm_add_epi32(sum32, err0);
537       }
538       for (k = j; k < width; ++k) {
539         const int32_t u = (int32_t)(dat[k] << SGRPROJ_RST_BITS);
540         int32_t v = xq[0] * (flt0[k] - u) + xq[1] * (flt1[k] - u);
541         const int32_t e = ROUND_POWER_OF_TWO(v, shift) + dat[k] - src[k];
542         err += ((int64_t)e * e);
543       }
544       dat += dat_stride;
545       src += src_stride;
546       flt0 += flt0_stride;
547       flt1 += flt1_stride;
548       const __m128i sum64_0 = _mm_cvtepi32_epi64(sum32);
549       const __m128i sum64_1 = _mm_cvtepi32_epi64(_mm_srli_si128(sum32, 8));
550       sum64 = _mm_add_epi64(sum64, sum64_0);
551       sum64 = _mm_add_epi64(sum64, sum64_1);
552     }
553   } else if (params->r[0] > 0 || params->r[1] > 0) {
554     const int xq_active = (params->r[0] > 0) ? xq[0] : xq[1];
555     const __m128i xq_coeff =
556         pair_set_epi16(xq_active, -(xq_active << SGRPROJ_RST_BITS));
557     const int32_t *flt = (params->r[0] > 0) ? flt0 : flt1;
558     const int flt_stride = (params->r[0] > 0) ? flt0_stride : flt1_stride;
559     for (i = 0; i < height; ++i) {
560       __m128i sum32 = _mm_setzero_si128();
561       for (j = 0; j <= width - 8; j += 8) {
562         const __m128i d0 = _mm_cvtepu8_epi16(xx_loadl_64(dat + j));
563         const __m128i s0 = _mm_cvtepu8_epi16(xx_loadl_64(src + j));
564         const __m128i flt_16b =
565             _mm_packs_epi32(xx_loadu_128(flt + j), xx_loadu_128(flt + j + 4));
566         const __m128i v0 =
567             _mm_madd_epi16(xq_coeff, _mm_unpacklo_epi16(flt_16b, d0));
568         const __m128i v1 =
569             _mm_madd_epi16(xq_coeff, _mm_unpackhi_epi16(flt_16b, d0));
570         const __m128i vr0 = _mm_srai_epi32(_mm_add_epi32(v0, rounding), shift);
571         const __m128i vr1 = _mm_srai_epi32(_mm_add_epi32(v1, rounding), shift);
572         const __m128i e0 =
573             _mm_sub_epi16(_mm_add_epi16(_mm_packs_epi32(vr0, vr1), d0), s0);
574         const __m128i err0 = _mm_madd_epi16(e0, e0);
575         sum32 = _mm_add_epi32(sum32, err0);
576       }
577       for (k = j; k < width; ++k) {
578         const int32_t u = (int32_t)(dat[k] << SGRPROJ_RST_BITS);
579         int32_t v = xq_active * (flt[k] - u);
580         const int32_t e = ROUND_POWER_OF_TWO(v, shift) + dat[k] - src[k];
581         err += ((int64_t)e * e);
582       }
583       dat += dat_stride;
584       src += src_stride;
585       flt += flt_stride;
586       const __m128i sum64_0 = _mm_cvtepi32_epi64(sum32);
587       const __m128i sum64_1 = _mm_cvtepi32_epi64(_mm_srli_si128(sum32, 8));
588       sum64 = _mm_add_epi64(sum64, sum64_0);
589       sum64 = _mm_add_epi64(sum64, sum64_1);
590     }
591   } else {
592     __m128i sum32 = _mm_setzero_si128();
593     for (i = 0; i < height; ++i) {
594       for (j = 0; j <= width - 16; j += 16) {
595         const __m128i d = xx_loadu_128(dat + j);
596         const __m128i s = xx_loadu_128(src + j);
597         const __m128i d0 = _mm_cvtepu8_epi16(d);
598         const __m128i d1 = _mm_cvtepu8_epi16(_mm_srli_si128(d, 8));
599         const __m128i s0 = _mm_cvtepu8_epi16(s);
600         const __m128i s1 = _mm_cvtepu8_epi16(_mm_srli_si128(s, 8));
601         const __m128i diff0 = _mm_sub_epi16(d0, s0);
602         const __m128i diff1 = _mm_sub_epi16(d1, s1);
603         const __m128i err0 = _mm_madd_epi16(diff0, diff0);
604         const __m128i err1 = _mm_madd_epi16(diff1, diff1);
605         sum32 = _mm_add_epi32(sum32, err0);
606         sum32 = _mm_add_epi32(sum32, err1);
607       }
608       for (k = j; k < width; ++k) {
609         const int32_t e = (int32_t)(dat[k]) - src[k];
610         err += ((int64_t)e * e);
611       }
612       dat += dat_stride;
613       src += src_stride;
614     }
615     const __m128i sum64_0 = _mm_cvtepi32_epi64(sum32);
616     const __m128i sum64_1 = _mm_cvtepi32_epi64(_mm_srli_si128(sum32, 8));
617     sum64 = _mm_add_epi64(sum64_0, sum64_1);
618   }
619   int64_t sum[2];
620   xx_storeu_128(sum, sum64);
621   err += sum[0] + sum[1];
622   return err;
623 }
624 
av1_highbd_pixel_proj_error_sse4_1(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)625 int64_t av1_highbd_pixel_proj_error_sse4_1(
626     const uint8_t *src8, int width, int height, int src_stride,
627     const uint8_t *dat8, int dat_stride, int32_t *flt0, int flt0_stride,
628     int32_t *flt1, int flt1_stride, int xq[2], const sgr_params_type *params) {
629   int i, j, k;
630   const int32_t shift = SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS;
631   const __m128i rounding = _mm_set1_epi32(1 << (shift - 1));
632   __m128i sum64 = _mm_setzero_si128();
633   const uint16_t *src = CONVERT_TO_SHORTPTR(src8);
634   const uint16_t *dat = CONVERT_TO_SHORTPTR(dat8);
635   int64_t err = 0;
636   if (params->r[0] > 0 && params->r[1] > 0) {  // Both filters are enabled
637     const __m128i xq0 = _mm_set1_epi32(xq[0]);
638     const __m128i xq1 = _mm_set1_epi32(xq[1]);
639 
640     for (i = 0; i < height; ++i) {
641       __m128i sum32 = _mm_setzero_si128();
642       for (j = 0; j <= width - 8; j += 8) {
643         // Load 8x pixels from source image
644         const __m128i s0 = xx_loadu_128(src + j);
645         // s0 = [7 6 5 4 3 2 1 0] as i16 (indices of src[])
646 
647         // Load 8x pixels from corrupted image
648         const __m128i d0 = xx_loadu_128(dat + j);
649         // d0 = [7 6 5 4 3 2 1 0] as i16 (indices of dat[])
650 
651         // Shift each pixel value up by SGRPROJ_RST_BITS
652         const __m128i u0 = _mm_slli_epi16(d0, SGRPROJ_RST_BITS);
653 
654         // Split u0 into two halves and pad each from u16 to i32
655         const __m128i u0l = _mm_cvtepu16_epi32(u0);
656         const __m128i u0h = _mm_cvtepu16_epi32(_mm_srli_si128(u0, 8));
657         // u0h = [7 6 5 4] as i32, u0l = [3 2 1 0] as i32, all dat[] indices
658 
659         // Load 8 pixels from first and second filtered images
660         const __m128i flt0l = xx_loadu_128(flt0 + j);
661         const __m128i flt0h = xx_loadu_128(flt0 + j + 4);
662         const __m128i flt1l = xx_loadu_128(flt1 + j);
663         const __m128i flt1h = xx_loadu_128(flt1 + j + 4);
664         // flt0 = [7 6 5 4] [3 2 1 0] as i32 (indices of flt0+j)
665         // flt1 = [7 6 5 4] [3 2 1 0] as i32 (indices of flt1+j)
666 
667         // Subtract shifted corrupt image from each filtered image
668         // This gives our two basis vectors for the projection
669         const __m128i flt0l_subu = _mm_sub_epi32(flt0l, u0l);
670         const __m128i flt0h_subu = _mm_sub_epi32(flt0h, u0h);
671         const __m128i flt1l_subu = _mm_sub_epi32(flt1l, u0l);
672         const __m128i flt1h_subu = _mm_sub_epi32(flt1h, u0h);
673         // flt?h_subu = [ f[7]-u[7] f[6]-u[6] f[5]-u[5] f[4]-u[4] ] as i32
674         // flt?l_subu = [ f[3]-u[3] f[2]-u[2] f[1]-u[1] f[0]-u[0] ] as i32
675 
676         // Multiply each basis vector by the corresponding coefficient
677         const __m128i v0l = _mm_mullo_epi32(flt0l_subu, xq0);
678         const __m128i v0h = _mm_mullo_epi32(flt0h_subu, xq0);
679         const __m128i v1l = _mm_mullo_epi32(flt1l_subu, xq1);
680         const __m128i v1h = _mm_mullo_epi32(flt1h_subu, xq1);
681 
682         // Add together the contribution from each scaled basis vector
683         const __m128i vl = _mm_add_epi32(v0l, v1l);
684         const __m128i vh = _mm_add_epi32(v0h, v1h);
685 
686         // Right-shift v with appropriate rounding
687         const __m128i vrl = _mm_srai_epi32(_mm_add_epi32(vl, rounding), shift);
688         const __m128i vrh = _mm_srai_epi32(_mm_add_epi32(vh, rounding), shift);
689 
690         // Saturate each i32 value to i16 and combine lower and upper halves
691         const __m128i vr = _mm_packs_epi32(vrl, vrh);
692 
693         // Add twin-subspace-sgr-filter to corrupt image then subtract source
694         const __m128i e0 = _mm_sub_epi16(_mm_add_epi16(vr, d0), s0);
695 
696         // Calculate squared error and add adjacent values
697         const __m128i err0 = _mm_madd_epi16(e0, e0);
698 
699         sum32 = _mm_add_epi32(sum32, err0);
700       }
701 
702       const __m128i sum32l = _mm_cvtepu32_epi64(sum32);
703       sum64 = _mm_add_epi64(sum64, sum32l);
704       const __m128i sum32h = _mm_cvtepu32_epi64(_mm_srli_si128(sum32, 8));
705       sum64 = _mm_add_epi64(sum64, sum32h);
706 
707       // Process remaining pixels in this row (modulo 8)
708       for (k = j; k < width; ++k) {
709         const int32_t u = (int32_t)(dat[k] << SGRPROJ_RST_BITS);
710         int32_t v = xq[0] * (flt0[k] - u) + xq[1] * (flt1[k] - u);
711         const int32_t e = ROUND_POWER_OF_TWO(v, shift) + dat[k] - src[k];
712         err += ((int64_t)e * e);
713       }
714       dat += dat_stride;
715       src += src_stride;
716       flt0 += flt0_stride;
717       flt1 += flt1_stride;
718     }
719   } else if (params->r[0] > 0 || params->r[1] > 0) {  // Only one filter enabled
720     const int32_t xq_on = (params->r[0] > 0) ? xq[0] : xq[1];
721     const __m128i xq_active = _mm_set1_epi32(xq_on);
722     const __m128i xq_inactive =
723         _mm_set1_epi32(-xq_on * (1 << SGRPROJ_RST_BITS));
724     const int32_t *flt = (params->r[0] > 0) ? flt0 : flt1;
725     const int flt_stride = (params->r[0] > 0) ? flt0_stride : flt1_stride;
726     for (i = 0; i < height; ++i) {
727       __m128i sum32 = _mm_setzero_si128();
728       for (j = 0; j <= width - 8; j += 8) {
729         // Load 8x pixels from source image
730         const __m128i s0 = xx_loadu_128(src + j);
731         // s0 = [7 6 5 4 3 2 1 0] as u16 (indices of src[])
732 
733         // Load 8x pixels from corrupted image and pad each u16 to i32
734         const __m128i d0 = xx_loadu_128(dat + j);
735         const __m128i d0h = _mm_cvtepu16_epi32(_mm_srli_si128(d0, 8));
736         const __m128i d0l = _mm_cvtepu16_epi32(d0);
737         // d0h, d0l = [7 6 5 4], [3 2 1 0] as u32 (indices of dat[])
738 
739         // Load 8 pixels from the filtered image
740         const __m128i flth = xx_loadu_128(flt + j + 4);
741         const __m128i fltl = xx_loadu_128(flt + j);
742         // flth, fltl = [7 6 5 4], [3 2 1 0] as i32 (indices of flt+j)
743 
744         const __m128i flth_xq = _mm_mullo_epi32(flth, xq_active);
745         const __m128i fltl_xq = _mm_mullo_epi32(fltl, xq_active);
746         const __m128i d0h_xq = _mm_mullo_epi32(d0h, xq_inactive);
747         const __m128i d0l_xq = _mm_mullo_epi32(d0l, xq_inactive);
748 
749         const __m128i vh = _mm_add_epi32(flth_xq, d0h_xq);
750         const __m128i vl = _mm_add_epi32(fltl_xq, d0l_xq);
751         // vh = [ xq0(f[7]-d[7]) xq0(f[6]-d[6]) xq0(f[5]-d[5]) xq0(f[4]-d[4]) ]
752         // vl = [ xq0(f[3]-d[3]) xq0(f[2]-d[2]) xq0(f[1]-d[1]) xq0(f[0]-d[0]) ]
753 
754         // Shift this down with appropriate rounding
755         const __m128i vrh = _mm_srai_epi32(_mm_add_epi32(vh, rounding), shift);
756         const __m128i vrl = _mm_srai_epi32(_mm_add_epi32(vl, rounding), shift);
757 
758         // Saturate vr0 and vr1 from i32 to i16 then pack together
759         const __m128i vr = _mm_packs_epi32(vrl, vrh);
760 
761         // Subtract twin-subspace-sgr filtered from source image to get error
762         const __m128i e0 = _mm_sub_epi16(_mm_add_epi16(vr, d0), s0);
763 
764         // Calculate squared error and add adjacent values
765         const __m128i err0 = _mm_madd_epi16(e0, e0);
766 
767         sum32 = _mm_add_epi32(sum32, err0);
768       }
769 
770       const __m128i sum32l = _mm_cvtepu32_epi64(sum32);
771       sum64 = _mm_add_epi64(sum64, sum32l);
772       const __m128i sum32h = _mm_cvtepu32_epi64(_mm_srli_si128(sum32, 8));
773       sum64 = _mm_add_epi64(sum64, sum32h);
774 
775       // Process remaining pixels in this row (modulo 8)
776       for (k = j; k < width; ++k) {
777         const int32_t u = (int32_t)(dat[k] << SGRPROJ_RST_BITS);
778         int32_t v = xq_on * (flt[k] - u);
779         const int32_t e = ROUND_POWER_OF_TWO(v, shift) + dat[k] - src[k];
780         err += ((int64_t)e * e);
781       }
782       dat += dat_stride;
783       src += src_stride;
784       flt += flt_stride;
785     }
786   } else {  // Neither filter is enabled
787     for (i = 0; i < height; ++i) {
788       __m128i sum32 = _mm_setzero_si128();
789       for (j = 0; j <= width - 16; j += 16) {
790         // Load 2x8 u16 from source image
791         const __m128i s0 = xx_loadu_128(src + j);
792         const __m128i s1 = xx_loadu_128(src + j + 8);
793         // Load 2x8 u16 from corrupted image
794         const __m128i d0 = xx_loadu_128(dat + j);
795         const __m128i d1 = xx_loadu_128(dat + j + 8);
796 
797         // Subtract corrupted image from source image
798         const __m128i diff0 = _mm_sub_epi16(d0, s0);
799         const __m128i diff1 = _mm_sub_epi16(d1, s1);
800 
801         // Square error and add adjacent values
802         const __m128i err0 = _mm_madd_epi16(diff0, diff0);
803         const __m128i err1 = _mm_madd_epi16(diff1, diff1);
804 
805         sum32 = _mm_add_epi32(sum32, err0);
806         sum32 = _mm_add_epi32(sum32, err1);
807       }
808 
809       const __m128i sum32l = _mm_cvtepu32_epi64(sum32);
810       sum64 = _mm_add_epi64(sum64, sum32l);
811       const __m128i sum32h = _mm_cvtepu32_epi64(_mm_srli_si128(sum32, 8));
812       sum64 = _mm_add_epi64(sum64, sum32h);
813 
814       // Process remaining pixels (modulu 8)
815       for (k = j; k < width; ++k) {
816         const int32_t e = (int32_t)(dat[k]) - src[k];
817         err += ((int64_t)e * e);
818       }
819       dat += dat_stride;
820       src += src_stride;
821     }
822   }
823 
824   // Sum 4 values from sum64l and sum64h into err
825   int64_t sum[2];
826   xx_storeu_128(sum, sum64);
827   err += sum[0] + sum[1];
828   return err;
829 }
830