/*
* Copyright (c) 2020, Alliance for Open Media. All rights reserved.
*
* This source code is subject to the terms of the BSD 2 Clause License and
* the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
* was not distributed with this source code in the LICENSE file, you can
* obtain it at www.aomedia.org/license/software. If the Alliance for Open
* Media Patent License 1.0 was not distributed with this source code in the
* PATENTS file, you can obtain it at www.aomedia.org/license/patent.
*/
#include <arm_neon.h>
#include "config/aom_config.h"
#include "config/av1_rtcd.h"
#include "aom_dsp/arm/mem_neon.h"
#include "aom_dsp/arm/sum_neon.h"
#include "aom_dsp/arm/transpose_neon.h"
#include "av1/common/restoration.h"
#include "av1/encoder/arm/pickrst_neon.h"
#include "av1/encoder/pickrst.h"
int64_t av1_lowbd_pixel_proj_error_neon(
const uint8_t *src,
int width,
int height,
int src_stride,
const uint8_t *dat,
int dat_stride, int32_t *flt0,
int flt0_stride,
int32_t *flt1,
int flt1_stride,
int xq[2],
const sgr_params_type *params) {
int64_t sse = 0;
int64x2_t sse_s64 = vdupq_n_s64(0);
if (params->r[0] > 0 && params->r[1] > 0) {
int32x2_t xq_v = vld1_s32(xq);
int32x2_t xq_sum_v = vshl_n_s32(vpadd_s32(xq_v, xq_v), SGRPROJ_RST_BITS);
do {
int j = 0;
int32x4_t sse_s32 = vdupq_n_s32(0);
do {
const uint8x8_t d = vld1_u8(&dat[j]);
const uint8x8_t s = vld1_u8(&src[j]);
int32x4_t flt0_0 = vld1q_s32(&flt0[j]);
int32x4_t flt0_1 = vld1q_s32(&flt0[j + 4]);
int32x4_t flt1_0 = vld1q_s32(&flt1[j]);
int32x4_t flt1_1 = vld1q_s32(&flt1[j + 4]);
int32x4_t offset =
vdupq_n_s32(1 << (SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS - 1));
int32x4_t v0 = vmlaq_lane_s32(offset, flt0_0, xq_v, 0);
int32x4_t v1 = vmlaq_lane_s32(offset, flt0_1, xq_v, 0);
v0 = vmlaq_lane_s32(v0, flt1_0, xq_v, 1);
v1 = vmlaq_lane_s32(v1, flt1_1, xq_v, 1);
int16x8_t d_s16 = vreinterpretq_s16_u16(vmovl_u8(d));
v0 = vmlsl_lane_s16(v0, vget_low_s16(d_s16),
vreinterpret_s16_s32(xq_sum_v), 0);
v1 = vmlsl_lane_s16(v1, vget_high_s16(d_s16),
vreinterpret_s16_s32(xq_sum_v), 0);
int16x4_t vr0 = vshrn_n_s32(v0, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS);
int16x4_t vr1 = vshrn_n_s32(v1, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS);
int16x8_t diff = vreinterpretq_s16_u16(vsubl_u8(d, s));
int16x8_t e = vaddq_s16(vcombine_s16(vr0, vr1), diff);
int16x4_t e_lo = vget_low_s16(e);
int16x4_t e_hi = vget_high_s16(e);
sse_s32 = vmlal_s16(sse_s32, e_lo, e_lo);
sse_s32 = vmlal_s16(sse_s32, e_hi, e_hi);
j += 8;
}
while (j <= width - 8);
for (
int k = j; k < width; ++k) {
int32_t u = (dat[k] << SGRPROJ_RST_BITS);
int32_t v = (1 << (SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS - 1)) +
xq[0] * flt0[k] + xq[1] * flt1[k] - u * (xq[0] + xq[1]);
int32_t e =
(v >> (SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS)) + dat[k] - src[k];
sse += e * e;
}
sse_s64 = vpadalq_s32(sse_s64, sse_s32);
dat += dat_stride;
src += src_stride;
flt0 += flt0_stride;
flt1 += flt1_stride;
}
while (--height != 0);
}
else if (params->r[0] > 0 || params->r[1] > 0) {
int xq_active = (params->r[0] > 0) ? xq[0] : xq[1];
int32_t *flt = (params->r[0] > 0) ? flt0 : flt1;
int flt_stride = (params->r[0] > 0) ? flt0_stride : flt1_stride;
int32x2_t xq_v = vdup_n_s32(xq_active);
do {
int32x4_t sse_s32 = vdupq_n_s32(0);
int j = 0;
do {
const uint8x8_t d = vld1_u8(&dat[j]);
const uint8x8_t s = vld1_u8(&src[j]);
int32x4_t flt_0 = vld1q_s32(&flt[j]);
int32x4_t flt_1 = vld1q_s32(&flt[j + 4]);
int16x8_t d_s16 =
vreinterpretq_s16_u16(vshll_n_u8(d, SGRPROJ_RST_BITS));
int32x4_t sub_0 = vsubw_s16(flt_0, vget_low_s16(d_s16));
int32x4_t sub_1 = vsubw_s16(flt_1, vget_high_s16(d_s16));
int32x4_t offset =
vdupq_n_s32(1 << (SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS - 1));
int32x4_t v0 = vmlaq_lane_s32(offset, sub_0, xq_v, 0);
int32x4_t v1 = vmlaq_lane_s32(offset, sub_1, xq_v, 0);
int16x4_t vr0 = vshrn_n_s32(v0, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS);
int16x4_t vr1 = vshrn_n_s32(v1, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS);
int16x8_t diff = vreinterpretq_s16_u16(vsubl_u8(d, s));
int16x8_t e = vaddq_s16(vcombine_s16(vr0, vr1), diff);
int16x4_t e_lo = vget_low_s16(e);
int16x4_t e_hi = vget_high_s16(e);
sse_s32 = vmlal_s16(sse_s32, e_lo, e_lo);
sse_s32 = vmlal_s16(sse_s32, e_hi, e_hi);
j += 8;
}
while (j <= width - 8);
for (
int k = j; k < width; ++k) {
int32_t u = dat[k] << SGRPROJ_RST_BITS;
int32_t v = xq_active * (flt[k] - u);
int32_t e = ROUND_POWER_OF_TWO(v, SGRPROJ_RST_BITS + SGRPROJ_PRJ_BITS) +
dat[k] - src[k];
sse += e * e;
}
sse_s64 = vpadalq_s32(sse_s64, sse_s32);
dat += dat_stride;
src += src_stride;
flt += flt_stride;
}
while (--height != 0);
}
else {
uint32x4_t sse_s32 = vdupq_n_u32(0);
do {
int j = 0;
do {
const uint8x16_t d = vld1q_u8(&dat[j]);
const uint8x16_t s = vld1q_u8(&src[j]);
uint8x16_t diff = vabdq_u8(d, s);
uint8x8_t diff_lo = vget_low_u8(diff);
uint8x8_t diff_hi = vget_high_u8(diff);
sse_s32 = vpadalq_u16(sse_s32, vmull_u8(diff_lo, diff_lo));
sse_s32 = vpadalq_u16(sse_s32, vmull_u8(diff_hi, diff_hi));
j += 16;
}
while (j <= width - 16);
for (
int k = j; k < width; ++k) {
int32_t e = dat[k] - src[k];
sse += e * e;
}
dat += dat_stride;
src += src_stride;
}
while (--height != 0);
sse_s64 = vreinterpretq_s64_u64(vpaddlq_u32(sse_s32));
}
sse += horizontal_add_s64x2(sse_s64);
return sse;
}
// We can accumulate up to 32768 8-bit multiplication results in a signed
// 32-bit integer. We are processing 2 pixels at a time, so the accumulator max
// can be as high as 16384 for the compute stats.
#define STAT_ACCUMULATOR_MAX 16384
static inline uint8x8_t tbl2(uint8x16_t a, uint8x16_t b, uint8x8_t idx) {
#if AOM_ARCH_AARCH64
uint8x16x2_t table = { { a, b } };
return vqtbl2_u8(table, idx);
#else
uint8x8x4_t table = { { vget_low_u8(a), vget_high_u8(a), vget_low_u8(b),
vget_high_u8(b) } };
return vtbl4_u8(table, idx);
#endif
}
static inline uint8x16_t tbl2q(uint8x16_t a, uint8x16_t b, uint8x16_t idx) {
#if AOM_ARCH_AARCH64
uint8x16x2_t table = { { a, b } };
return vqtbl2q_u8(table, idx);
#else
uint8x8x4_t table = { { vget_low_u8(a), vget_high_u8(a), vget_low_u8(b),
vget_high_u8(b) } };
return vcombine_u8(vtbl4_u8(table, vget_low_u8(idx)),
vtbl4_u8(table, vget_high_u8(idx)));
#endif
}
// The M matrix is accumulated in STAT_ACCUMULATOR_MAX steps to speed-up the
// computation. This function computes the final M from the accumulated
// (src_s64) and the residual parts (src_s32). It also transposes the result as
// the output needs to be column-major.
static inline void acc_transpose_M(int64_t *dst,
const int64_t *src_s64,
const int32_t *src_s32,
const int wiener_win,
int scale) {
for (
int i = 0; i < wiener_win; ++i) {
for (
int j = 0; j < wiener_win; ++j) {
int tr_idx = j * wiener_win + i;
*dst++ += (int64_t)(src_s64[tr_idx] + src_s32[tr_idx]) * scale;
}
}
}
// The resulting H is a column-major matrix accumulated from the transposed
// (column-major) samples of the filter kernel (5x5 or 7x7) viewed as a single
// vector. For the 7x7 filter case: H(49x49) = [49 x 1] x [1 x 49]. This
// function transforms back to the originally expected format (double
// transpose). The H matrix is accumulated in STAT_ACCUMULATOR_MAX steps to
// speed-up the computation. This function computes the final H from the
// accumulated (src_s64) and the residual parts (src_s32). The computed H is
// only an upper triangle matrix, this function also fills the lower triangle of
// the resulting matrix.
static void update_H(int64_t *dst,
const int64_t *src_s64,
const int32_t *src_s32,
const int wiener_win,
int stride,
int scale) {
// For a simplified theoretical 3x3 case where `wiener_win` is 3 and
// `wiener_win2` is 9, the M matrix is 3x3:
// 0, 3, 6
// 1, 4, 7
// 2, 5, 8
//
// This is viewed as a vector to compute H (9x9) by vector outer product:
// 0, 3, 6, 1, 4, 7, 2, 5, 8
//
// Double transpose and upper triangle remapping for 3x3 -> 9x9 case:
// 0, 3, 6, 1, 4, 7, 2, 5, 8,
// 3, 30, 33, 12, 31, 34, 21, 32, 35,
// 6, 33, 60, 15, 42, 61, 24, 51, 62,
// 1, 12, 15, 10, 13, 16, 11, 14, 17,
// 4, 31, 42, 13, 40, 43, 22, 41, 44,
// 7, 34, 61, 16, 43, 70, 25, 52, 71,
// 2, 21, 24, 11, 22, 25, 20, 23, 26,
// 5, 32, 51, 14, 41, 52, 23, 50, 53,
// 8, 35, 62, 17, 44, 71, 26, 53, 80,
const int wiener_win2 = wiener_win * wiener_win;
// Loop through the indices according to the remapping above, along the
// columns:
// 0, wiener_win, 2 * wiener_win, ..., 1, 1 + 2 * wiener_win, ...,
// wiener_win - 1, wiener_win - 1 + wiener_win, ...
// For the 3x3 case `j` will be: 0, 3, 6, 1, 4, 7, 2, 5, 8.
for (
int i = 0; i < wiener_win; ++i) {
for (
int j = i; j < wiener_win2; j += wiener_win) {
// These two inner loops are the same as the two outer loops, but running
// along rows instead of columns. For the 3x3 case `l` will be:
// 0, 3, 6, 1, 4, 7, 2, 5, 8.
for (
int k = 0; k < wiener_win; ++k) {
for (
int l = k; l < wiener_win2; l += wiener_win) {
// The nominal double transpose indexing would be:
// int idx = stride * j + l;
// However we need the upper-triangle indices, it is easy with some
// min/max operations.
int tr_idx = stride * AOMMIN(j, l) + AOMMAX(j, l);
// Resulting matrix is filled by combining the 64-bit and the residual
// 32-bit matrices together with scaling.
*dst++ += (int64_t)(src_s64[tr_idx] + src_s32[tr_idx]) * scale;
}
}
}
}
}
// Load 7x7 matrix into 3 and a half 128-bit vectors from consecutive rows, the
// last load address is offset to prevent out-of-bounds access.
static inline void load_and_pack_u8_8x7(uint8x16_t dst[4],
const uint8_t *src,
ptrdiff_t stride) {
dst[0] = vcombine_u8(vld1_u8(src), vld1_u8(src + stride));
src += 2 * stride;
dst[1] = vcombine_u8(vld1_u8(src), vld1_u8(src + stride));
src += 2 * stride;
dst[2] = vcombine_u8(vld1_u8(src), vld1_u8(src + stride));
src += 2 * stride;
dst[3] = vcombine_u8(vld1_u8(src - 1), vdup_n_u8(0));
}
static inline void compute_stats_win7_downsampled_neon(
const uint8_t *dgd,
const uint8_t *src,
int width,
int height,
int dgd_stride,
int src_stride,
int avg, int64_t *M, int64_t *H,
int downsample_factor) {
// Matrix names are capitalized to help readability.
DECLARE_ALIGNED(64, int16_t, DGD_AVG0[WIENER_WIN2_ALIGN3]);
DECLARE_ALIGNED(64, int16_t, DGD_AVG1[WIENER_WIN2_ALIGN3]);
DECLARE_ALIGNED(64, int32_t, M_s32[WIENER_WIN2_ALIGN3]);
DECLARE_ALIGNED(64, int64_t, M_s64[WIENER_WIN2_ALIGN3]);
DECLARE_ALIGNED(64, int32_t, H_s32[WIENER_WIN2 * WIENER_WIN2_ALIGN2]);
DECLARE_ALIGNED(64, int64_t, H_s64[WIENER_WIN2 * WIENER_WIN2_ALIGN2]);
memset(M_s32, 0,
sizeof(M_s32));
memset(M_s64, 0,
sizeof(M_s64));
memset(H_s32, 0,
sizeof(H_s32));
memset(H_s64, 0,
sizeof(H_s64));
// Look-up tables to create 8x6 matrix with consecutive elements from two 7x7
// matrices.
// clang-format off
DECLARE_ALIGNED(16,
static const uint8_t, shuffle_stats7[96]) = {
0, 1, 2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 16, 17,
2, 3, 4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 16, 17, 18, 19,
4, 5, 6, 8, 9, 10, 11, 12, 13, 14, 17, 18, 19, 20, 21, 22,
1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 17, 18,
3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20,
5, 6, 7, 9, 10, 11, 12, 13, 14, 15, 18, 19, 20, 21, 22, 23,
};
// clang-format on
const uint8x16_t lut0 = vld1q_u8(shuffle_stats7 + 0);
const uint8x16_t lut1 = vld1q_u8(shuffle_stats7 + 16);
const uint8x16_t lut2 = vld1q_u8(shuffle_stats7 + 32);
const uint8x16_t lut3 = vld1q_u8(shuffle_stats7 + 48);
const uint8x16_t lut4 = vld1q_u8(shuffle_stats7 + 64);
const uint8x16_t lut5 = vld1q_u8(shuffle_stats7 + 80);
int acc_cnt = STAT_ACCUMULATOR_MAX;
const int src_next = downsample_factor * src_stride - width;
const int dgd_next = downsample_factor * dgd_stride - width;
const uint8x8_t avg_u8 = vdup_n_u8(avg);
do {
int j = width;
while (j >= 2) {
// Load two adjacent, overlapping 7x7 matrices: a 8x7 matrix with the
// middle 6x7 elements being shared.
uint8x16_t dgd_rows[4];
load_and_pack_u8_8x7(dgd_rows, dgd, dgd_stride);
const uint8_t *dgd_ptr = dgd + dgd_stride * 6;
dgd += 2;
// Re-arrange (and widen) the combined 8x7 matrix to have the 2 whole 7x7
// matrices (1 for each of the 2 pixels) separated into distinct
// int16x8_t[6] arrays. These arrays contain 48 elements of the 49 (7x7).
// Compute `dgd - avg` for both buffers. Each DGD_AVG buffer contains 49
// consecutive elements.
int16x8_t dgd_avg0[6];
int16x8_t dgd_avg1[6];
uint8x16_t dgd_shuf0 = tbl2q(dgd_rows[0], dgd_rows[1], lut0);
uint8x16_t dgd_shuf3 = tbl2q(dgd_rows[0], dgd_rows[1], lut3);
dgd_avg0[0] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf0), avg_u8));
dgd_avg0[1] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf0), avg_u8));
dgd_avg1[0] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf3), avg_u8));
dgd_avg1[1] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf3), avg_u8));
vst1q_s16(DGD_AVG0, dgd_avg0[0]);
vst1q_s16(DGD_AVG0 + 8, dgd_avg0[1]);
vst1q_s16(DGD_AVG1, dgd_avg1[0]);
vst1q_s16(DGD_AVG1 + 8, dgd_avg1[1]);
uint8x16_t dgd_shuf1 = tbl2q(dgd_rows[1], dgd_rows[2], lut1);
uint8x16_t dgd_shuf4 = tbl2q(dgd_rows[1], dgd_rows[2], lut4);
dgd_avg0[2] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf1), avg_u8));
dgd_avg0[3] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf1), avg_u8));
dgd_avg1[2] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf4), avg_u8));
dgd_avg1[3] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf4), avg_u8));
vst1q_s16(DGD_AVG0 + 16, dgd_avg0[2]);
vst1q_s16(DGD_AVG0 + 24, dgd_avg0[3]);
vst1q_s16(DGD_AVG1 + 16, dgd_avg1[2]);
vst1q_s16(DGD_AVG1 + 24, dgd_avg1[3]);
uint8x16_t dgd_shuf2 = tbl2q(dgd_rows[2], dgd_rows[3], lut2);
uint8x16_t dgd_shuf5 = tbl2q(dgd_rows[2], dgd_rows[3], lut5);
dgd_avg0[4] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf2), avg_u8));
dgd_avg0[5] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf2), avg_u8));
dgd_avg1[4] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf5), avg_u8));
dgd_avg1[5] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf5), avg_u8));
vst1q_s16(DGD_AVG0 + 32, dgd_avg0[4]);
vst1q_s16(DGD_AVG0 + 40, dgd_avg0[5]);
vst1q_s16(DGD_AVG1 + 32, dgd_avg1[4]);
vst1q_s16(DGD_AVG1 + 40, dgd_avg1[5]);
// The remaining last (49th) elements of `dgd - avg`.
DGD_AVG0[48] = dgd_ptr[6] - avg;
DGD_AVG1[48] = dgd_ptr[7] - avg;
// Accumulate into row-major variant of matrix M (cross-correlation) for 2
// output pixels at a time. M is of size 7 * 7. It needs to be filled such
// that multiplying one element from src with each element of a row of the
// wiener window will fill one column of M. However this is not very
// convenient in terms of memory access, as it means we do contiguous
// loads of dgd but strided stores to M. As a result, we use an
// intermediate matrix M_s32 which is instead filled such that one row of
// the wiener window gives one row of M_s32. Once fully computed, M_s32 is
// then transposed to return M.
int src_avg0 = *src++ - avg;
int src_avg1 = *src++ - avg;
int16x4_t src_avg0_s16 = vdup_n_s16(src_avg0);
int16x4_t src_avg1_s16 = vdup_n_s16(src_avg1);
update_M_2pixels(M_s32 + 0, src_avg0_s16, src_avg1_s16, dgd_avg0[0],
dgd_avg1[0]);
update_M_2pixels(M_s32 + 8, src_avg0_s16, src_avg1_s16, dgd_avg0[1],
dgd_avg1[1]);
update_M_2pixels(M_s32 + 16, src_avg0_s16, src_avg1_s16, dgd_avg0[2],
dgd_avg1[2]);
update_M_2pixels(M_s32 + 24, src_avg0_s16, src_avg1_s16, dgd_avg0[3],
dgd_avg1[3]);
update_M_2pixels(M_s32 + 32, src_avg0_s16, src_avg1_s16, dgd_avg0[4],
dgd_avg1[4]);
update_M_2pixels(M_s32 + 40, src_avg0_s16, src_avg1_s16, dgd_avg0[5],
dgd_avg1[5]);
// Last (49th) element of M_s32 can be computed as scalar more efficiently
// for 2 output pixels.
M_s32[48] += DGD_AVG0[48] * src_avg0 + DGD_AVG1[48] * src_avg1;
// Start accumulating into row-major version of matrix H
// (auto-covariance), it expects the DGD_AVG[01] matrices to also be
// row-major. H is of size 49 * 49. It is filled by multiplying every pair
// of elements of the wiener window together (vector outer product). Since
// it is a symmetric matrix, we only compute the upper-right triangle, and
// then copy it down to the lower-left later. The upper triangle is
// covered by 4x4 tiles. The original algorithm assumes the M matrix is
// column-major and the resulting H matrix is also expected to be
// column-major. It is not efficient to work with column-major matrices,
// so we accumulate into a row-major matrix H_s32. At the end of the
// algorithm a double transpose transformation will convert H_s32 back to
// the expected output layout.
update_H_7x7_2pixels(H_s32, DGD_AVG0, DGD_AVG1);
// The last element of the triangle of H_s32 matrix can be computed as a
// scalar more efficiently.
H_s32[48 * WIENER_WIN2_ALIGN2 + 48] +=
DGD_AVG0[48] * DGD_AVG0[48] + DGD_AVG1[48] * DGD_AVG1[48];
// Accumulate into 64-bit after STAT_ACCUMULATOR_MAX iterations to prevent
// overflow.
if (--acc_cnt == 0) {
acc_cnt = STAT_ACCUMULATOR_MAX;
accumulate_and_clear(M_s64, M_s32, WIENER_WIN2_ALIGN2);
// The widening accumulation is only needed for the upper triangle part
// of the matrix.
int64_t *lh = H_s64;
int32_t *lh32 = H_s32;
for (
int k = 0; k < WIENER_WIN2; ++k) {
// The widening accumulation is only run for the relevant parts
// (upper-right triangle) in a row 4-element aligned.
int k4 = k / 4 * 4;
accumulate_and_clear(lh + k4, lh32 + k4, 48 - k4);
// Last element of the row is computed separately.
lh[48] += lh32[48];
lh32[48] = 0;
lh += WIENER_WIN2_ALIGN2;
lh32 += WIENER_WIN2_ALIGN2;
}
}
j -= 2;
}
// Computations for odd pixel in the row.
if (width & 1) {
// Load two adjacent, overlapping 7x7 matrices: a 8x7 matrix with the
// middle 6x7 elements being shared.
uint8x16_t dgd_rows[4];
load_and_pack_u8_8x7(dgd_rows, dgd, dgd_stride);
const uint8_t *dgd_ptr = dgd + dgd_stride * 6;
++dgd;
// Re-arrange (and widen) the combined 8x7 matrix to have a whole 7x7
// matrix tightly packed into a int16x8_t[6] array. This array contains
// 48 elements of the 49 (7x7). Compute `dgd - avg` for the whole buffer.
// The DGD_AVG buffer contains 49 consecutive elements.
int16x8_t dgd_avg0[6];
uint8x16_t dgd_shuf0 = tbl2q(dgd_rows[0], dgd_rows[1], lut0);
dgd_avg0[0] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf0), avg_u8));
dgd_avg0[1] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf0), avg_u8));
vst1q_s16(DGD_AVG0, dgd_avg0[0]);
vst1q_s16(DGD_AVG0 + 8, dgd_avg0[1]);
uint8x16_t dgd_shuf1 = tbl2q(dgd_rows[1], dgd_rows[2], lut1);
dgd_avg0[2] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf1), avg_u8));
dgd_avg0[3] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf1), avg_u8));
vst1q_s16(DGD_AVG0 + 16, dgd_avg0[2]);
vst1q_s16(DGD_AVG0 + 24, dgd_avg0[3]);
uint8x16_t dgd_shuf2 = tbl2q(dgd_rows[2], dgd_rows[3], lut2);
dgd_avg0[4] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf2), avg_u8));
dgd_avg0[5] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf2), avg_u8));
vst1q_s16(DGD_AVG0 + 32, dgd_avg0[4]);
vst1q_s16(DGD_AVG0 + 40, dgd_avg0[5]);
// The remaining last (49th) element of `dgd - avg`.
DGD_AVG0[48] = dgd_ptr[6] - avg;
// Accumulate into row-major order variant of matrix M (cross-correlation)
// for 1 output pixel at a time. M is of size 7 * 7. It needs to be filled
// such that multiplying one element from src with each element of a row
// of the wiener window will fill one column of M. However this is not
// very convenient in terms of memory access, as it means we do
// contiguous loads of dgd but strided stores to M. As a result, we use an
// intermediate matrix M_s32 which is instead filled such that one row of
// the wiener window gives one row of M_s32. Once fully computed, M_s32 is
// then transposed to return M.
int src_avg0 = *src++ - avg;
int16x4_t src_avg0_s16 = vdup_n_s16(src_avg0);
update_M_1pixel(M_s32 + 0, src_avg0_s16, dgd_avg0[0]);
update_M_1pixel(M_s32 + 8, src_avg0_s16, dgd_avg0[1]);
update_M_1pixel(M_s32 + 16, src_avg0_s16, dgd_avg0[2]);
update_M_1pixel(M_s32 + 24, src_avg0_s16, dgd_avg0[3]);
update_M_1pixel(M_s32 + 32, src_avg0_s16, dgd_avg0[4]);
update_M_1pixel(M_s32 + 40, src_avg0_s16, dgd_avg0[5]);
// Last (49th) element of M_s32 can be computed as scalar more efficiently
// for 1 output pixel.
M_s32[48] += DGD_AVG0[48] * src_avg0;
// Start accumulating into row-major order version of matrix H
// (auto-covariance), it expects the DGD_AVG0 matrix to also be row-major.
// H is of size 49 * 49. It is filled by multiplying every pair of
// elements of the wiener window together (vector outer product). Since it
// is a symmetric matrix, we only compute the upper-right triangle, and
// then copy it down to the lower-left later. The upper triangle is
// covered by 4x4 tiles. The original algorithm assumes the M matrix is
// column-major and the resulting H matrix is also expected to be
// column-major. It is not efficient to work column-major matrices, so we
// accumulate into a row-major matrix H_s32. At the end of the algorithm a
// double transpose transformation will convert H_s32 back to the expected
// output layout.
update_H_1pixel(H_s32, DGD_AVG0, WIENER_WIN2_ALIGN2, 48);
// The last element of the triangle of H_s32 matrix can be computed as
// scalar more efficiently.
H_s32[48 * WIENER_WIN2_ALIGN2 + 48] += DGD_AVG0[48] * DGD_AVG0[48];
}
src += src_next;
dgd += dgd_next;
}
while (--height != 0);
acc_transpose_M(M, M_s64, M_s32, WIENER_WIN, downsample_factor);
update_H(H, H_s64, H_s32, WIENER_WIN, WIENER_WIN2_ALIGN2, downsample_factor);
}
// Load 5x5 matrix into 2 and a half 128-bit vectors from consecutive rows, the
// last load address is offset to prevent out-of-bounds access.
static inline void load_and_pack_u8_6x5(uint8x16_t dst[3],
const uint8_t *src,
ptrdiff_t stride) {
dst[0] = vcombine_u8(vld1_u8(src), vld1_u8(src + stride));
src += 2 * stride;
dst[1] = vcombine_u8(vld1_u8(src), vld1_u8(src + stride));
src += 2 * stride;
dst[2] = vcombine_u8(vld1_u8(src - 3), vdup_n_u8(0));
}
static inline void compute_stats_win5_downsampled_neon(
const uint8_t *dgd,
const uint8_t *src,
int width,
int height,
int dgd_stride,
int src_stride,
int avg, int64_t *M, int64_t *H,
int downsample_factor) {
// Matrix names are capitalized to help readability.
DECLARE_ALIGNED(64, int16_t, DGD_AVG0[WIENER_WIN2_REDUCED_ALIGN3]);
DECLARE_ALIGNED(64, int16_t, DGD_AVG1[WIENER_WIN2_REDUCED_ALIGN3]);
DECLARE_ALIGNED(64, int32_t, M_s32[WIENER_WIN2_REDUCED_ALIGN3]);
DECLARE_ALIGNED(64, int64_t, M_s64[WIENER_WIN2_REDUCED_ALIGN3]);
DECLARE_ALIGNED(64, int32_t,
H_s32[WIENER_WIN2_REDUCED * WIENER_WIN2_REDUCED_ALIGN2]);
DECLARE_ALIGNED(64, int64_t,
H_s64[WIENER_WIN2_REDUCED * WIENER_WIN2_REDUCED_ALIGN2]);
memset(M_s32, 0,
sizeof(M_s32));
memset(M_s64, 0,
sizeof(M_s64));
memset(H_s32, 0,
sizeof(H_s32));
memset(H_s64, 0,
sizeof(H_s64));
// Look-up tables to create 8x3 matrix with consecutive elements from two 5x5
// matrices.
// clang-format off
DECLARE_ALIGNED(16,
static const uint8_t, shuffle_stats5[48]) = {
0, 1, 2, 3, 4, 8, 9, 10, 11, 12, 16, 17, 18, 19, 20, 24,
1, 2, 3, 4, 5, 9, 10, 11, 12, 13, 17, 18, 19, 20, 21, 25,
9, 10, 11, 12, 19, 20, 21, 22, 10, 11, 12, 13, 20, 21, 22, 23,
};
// clang-format on
const uint8x16_t lut0 = vld1q_u8(shuffle_stats5 + 0);
const uint8x16_t lut1 = vld1q_u8(shuffle_stats5 + 16);
const uint8x16_t lut2 = vld1q_u8(shuffle_stats5 + 32);
int acc_cnt = STAT_ACCUMULATOR_MAX;
const int src_next = downsample_factor * src_stride - width;
const int dgd_next = downsample_factor * dgd_stride - width;
const uint8x8_t avg_u8 = vdup_n_u8(avg);
do {
int j = width;
while (j >= 2) {
// Load two adjacent, overlapping 5x5 matrices: a 6x5 matrix with the
// middle 4x5 elements being shared.
uint8x16_t dgd_rows[3];
load_and_pack_u8_6x5(dgd_rows, dgd, dgd_stride);
const uint8_t *dgd_ptr = dgd + dgd_stride * 4;
dgd += 2;
// Re-arrange (and widen) the combined 6x5 matrix to have the 2 whole 5x5
// matrices (1 for each of the 2 pixels) separated into distinct
// int16x8_t[3] arrays. These arrays contain 24 elements of the 25 (5x5).
// Compute `dgd - avg` for both buffers. Each DGD_AVG buffer contains 25
// consecutive elements.
int16x8_t dgd_avg0[3];
int16x8_t dgd_avg1[3];
uint8x16_t dgd_shuf0 = tbl2q(dgd_rows[0], dgd_rows[1], lut0);
uint8x16_t dgd_shuf1 = tbl2q(dgd_rows[0], dgd_rows[1], lut1);
uint8x16_t dgd_shuf2 = tbl2q(dgd_rows[1], dgd_rows[2], lut2);
dgd_avg0[0] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf0), avg_u8));
dgd_avg0[1] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf0), avg_u8));
dgd_avg0[2] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf2), avg_u8));
dgd_avg1[0] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf1), avg_u8));
dgd_avg1[1] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf1), avg_u8));
dgd_avg1[2] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf2), avg_u8));
vst1q_s16(DGD_AVG0 + 0, dgd_avg0[0]);
vst1q_s16(DGD_AVG0 + 8, dgd_avg0[1]);
vst1q_s16(DGD_AVG0 + 16, dgd_avg0[2]);
vst1q_s16(DGD_AVG1 + 0, dgd_avg1[0]);
vst1q_s16(DGD_AVG1 + 8, dgd_avg1[1]);
vst1q_s16(DGD_AVG1 + 16, dgd_avg1[2]);
// The remaining last (25th) elements of `dgd - avg`.
DGD_AVG0[24] = dgd_ptr[4] - avg;
DGD_AVG1[24] = dgd_ptr[5] - avg;
// Accumulate into row-major variant of matrix M (cross-correlation) for 2
// output pixels at a time. M is of size 5 * 5. It needs to be filled such
// that multiplying one element from src with each element of a row of the
// wiener window will fill one column of M. However this is not very
// convenient in terms of memory access, as it means we do contiguous
// loads of dgd but strided stores to M. As a result, we use an
// intermediate matrix M_s32 which is instead filled such that one row of
// the wiener window gives one row of M_s32. Once fully computed, M_s32 is
// then transposed to return M.
int src_avg0 = *src++ - avg;
int src_avg1 = *src++ - avg;
int16x4_t src_avg0_s16 = vdup_n_s16(src_avg0);
int16x4_t src_avg1_s16 = vdup_n_s16(src_avg1);
update_M_2pixels(M_s32 + 0, src_avg0_s16, src_avg1_s16, dgd_avg0[0],
dgd_avg1[0]);
update_M_2pixels(M_s32 + 8, src_avg0_s16, src_avg1_s16, dgd_avg0[1],
dgd_avg1[1]);
update_M_2pixels(M_s32 + 16, src_avg0_s16, src_avg1_s16, dgd_avg0[2],
dgd_avg1[2]);
// Last (25th) element of M_s32 can be computed as scalar more efficiently
// for 2 output pixels.
M_s32[24] += DGD_AVG0[24] * src_avg0 + DGD_AVG1[24] * src_avg1;
// Start accumulating into row-major version of matrix H
// (auto-covariance), it expects the DGD_AVG[01] matrices to also be
// row-major. H is of size 25 * 25. It is filled by multiplying every pair
// of elements of the wiener window together (vector outer product). Since
// it is a symmetric matrix, we only compute the upper-right triangle, and
// then copy it down to the lower-left later. The upper triangle is
// covered by 4x4 tiles. The original algorithm assumes the M matrix is
// column-major and the resulting H matrix is also expected to be
// column-major. It is not efficient to work with column-major matrices,
// so we accumulate into a row-major matrix H_s32. At the end of the
// algorithm a double transpose transformation will convert H_s32 back to
// the expected output layout.
update_H_5x5_2pixels(H_s32, DGD_AVG0, DGD_AVG1);
// The last element of the triangle of H_s32 matrix can be computed as a
// scalar more efficiently.
H_s32[24 * WIENER_WIN2_REDUCED_ALIGN2 + 24] +=
DGD_AVG0[24] * DGD_AVG0[24] + DGD_AVG1[24] * DGD_AVG1[24];
// Accumulate into 64-bit after STAT_ACCUMULATOR_MAX iterations to prevent
// overflow.
if (--acc_cnt == 0) {
acc_cnt = STAT_ACCUMULATOR_MAX;
accumulate_and_clear(M_s64, M_s32, WIENER_WIN2_REDUCED_ALIGN2);
// The widening accumulation is only needed for the upper triangle part
// of the matrix.
int64_t *lh = H_s64;
int32_t *lh32 = H_s32;
for (
int k = 0; k < WIENER_WIN2_REDUCED; ++k) {
// The widening accumulation is only run for the relevant parts
// (upper-right triangle) in a row 4-element aligned.
int k4 = k / 4 * 4;
accumulate_and_clear(lh + k4, lh32 + k4, 24 - k4);
// Last element of the row is computed separately.
lh[24] += lh32[24];
lh32[24] = 0;
lh += WIENER_WIN2_REDUCED_ALIGN2;
lh32 += WIENER_WIN2_REDUCED_ALIGN2;
}
}
j -= 2;
}
// Computations for odd pixel in the row.
if (width & 1) {
// Load two adjacent, overlapping 5x5 matrices: a 6x5 matrix with the
// middle 4x5 elements being shared.
uint8x16_t dgd_rows[3];
load_and_pack_u8_6x5(dgd_rows, dgd, dgd_stride);
const uint8_t *dgd_ptr = dgd + dgd_stride * 4;
++dgd;
// Re-arrange (and widen) the combined 6x5 matrix to have a whole 5x5
// matrix tightly packed into a int16x8_t[3] array. This array contains
// 24 elements of the 25 (5x5). Compute `dgd - avg` for the whole buffer.
// The DGD_AVG buffer contains 25 consecutive elements.
int16x8_t dgd_avg0[3];
uint8x16_t dgd_shuf0 = tbl2q(dgd_rows[0], dgd_rows[1], lut0);
uint8x8_t dgd_shuf1 = tbl2(dgd_rows[1], dgd_rows[2], vget_low_u8(lut2));
dgd_avg0[0] =
vreinterpretq_s16_u16(vsubl_u8(vget_low_u8(dgd_shuf0), avg_u8));
dgd_avg0[1] =
vreinterpretq_s16_u16(vsubl_u8(vget_high_u8(dgd_shuf0), avg_u8));
dgd_avg0[2] = vreinterpretq_s16_u16(vsubl_u8(dgd_shuf1, avg_u8));
vst1q_s16(DGD_AVG0 + 0, dgd_avg0[0]);
vst1q_s16(DGD_AVG0 + 8, dgd_avg0[1]);
vst1q_s16(DGD_AVG0 + 16, dgd_avg0[2]);
// The remaining last (25th) element of `dgd - avg`.
DGD_AVG0[24] = dgd_ptr[4] - avg;
// Accumulate into row-major order variant of matrix M (cross-correlation)
// for 1 output pixel at a time. M is of size 5 * 5. It needs to be filled
// such that multiplying one element from src with each element of a row
// of the wiener window will fill one column of M. However this is not
// very convenient in terms of memory access, as it means we do
// contiguous loads of dgd but strided stores to M. As a result, we use an
// intermediate matrix M_s32 which is instead filled such that one row of
// the wiener window gives one row of M_s32. Once fully computed, M_s32 is
// then transposed to return M.
int src_avg0 = *src++ - avg;
int16x4_t src_avg0_s16 = vdup_n_s16(src_avg0);
update_M_1pixel(M_s32 + 0, src_avg0_s16, dgd_avg0[0]);
update_M_1pixel(M_s32 + 8, src_avg0_s16, dgd_avg0[1]);
update_M_1pixel(M_s32 + 16, src_avg0_s16, dgd_avg0[2]);
// Last (25th) element of M_s32 can be computed as scalar more efficiently
// for 1 output pixel.
M_s32[24] += DGD_AVG0[24] * src_avg0;
// Start accumulating into row-major order version of matrix H
// (auto-covariance), it expects the DGD_AVG0 matrix to also be row-major.
// H is of size 25 * 25. It is filled by multiplying every pair of
// elements of the wiener window together (vector outer product). Since it
// is a symmetric matrix, we only compute the upper-right triangle, and
// then copy it down to the lower-left later. The upper triangle is
// covered by 4x4 tiles. The original algorithm assumes the M matrix is
// column-major and the resulting H matrix is also expected to be
// column-major. It is not efficient to work column-major matrices, so we
// accumulate into a row-major matrix H_s32. At the end of the algorithm a
// double transpose transformation will convert H_s32 back to the expected
// output layout.
update_H_1pixel(H_s32, DGD_AVG0, WIENER_WIN2_REDUCED_ALIGN2, 24);
// The last element of the triangle of H_s32 matrix can be computed as a
// scalar more efficiently.
H_s32[24 * WIENER_WIN2_REDUCED_ALIGN2 + 24] +=
DGD_AVG0[24] * DGD_AVG0[24];
}
src += src_next;
dgd += dgd_next;
}
while (--height != 0);
acc_transpose_M(M, M_s64, M_s32, WIENER_WIN_REDUCED, downsample_factor);
update_H(H, H_s64, H_s32, WIENER_WIN_REDUCED, WIENER_WIN2_REDUCED_ALIGN2,
downsample_factor);
}
static inline void hadd_update_6_stats_neon(
const int64_t *
const src,
const int32x4_t *deltas,
int64_t *
const dst) {
int32x4_t delta01 = horizontal_add_2d_s32(deltas[0], deltas[1]);
int32x4_t delta23 = horizontal_add_2d_s32(deltas[2], deltas[3]);
int32x4_t delta45 = horizontal_add_2d_s32(deltas[4], deltas[5]);
int64x2_t delta01_s64 = vpaddlq_s32(delta01);
int64x2_t delta23_s64 = vpaddlq_s32(delta23);
int64x2_t delta45_s64 = vpaddlq_s32(delta45);
int64x2_t src0 = vld1q_s64(src);
int64x2_t src1 = vld1q_s64(src + 2);
int64x2_t src2 = vld1q_s64(src + 4);
vst1q_s64(dst, vaddq_s64(src0, delta01_s64));
vst1q_s64(dst + 2, vaddq_s64(src1, delta23_s64));
vst1q_s64(dst + 4, vaddq_s64(src2, delta45_s64));
}
static inline void hadd_update_4_stats_neon(
const int64_t *
const src,
const int32x4_t *deltas,
int64_t *
const dst) {
int32x4_t delta01 = horizontal_add_2d_s32(deltas[0], deltas[1]);
int32x4_t delta23 = horizontal_add_2d_s32(deltas[2], deltas[3]);
int64x2_t delta01_s64 = vpaddlq_s32(delta01);
int64x2_t delta23_s64 = vpaddlq_s32(delta23);
int64x2_t src0 = vld1q_s64(src);
int64x2_t src1 = vld1q_s64(src + 2);
vst1q_s64(dst, vaddq_s64(src0, delta01_s64));
vst1q_s64(dst + 2, vaddq_s64(src1, delta23_s64));
}
static inline void compute_stats_win5_neon(
const int16_t *
const d,
const int32_t d_stride,
const int16_t *
const s,
const int32_t s_stride,
const int32_t width,
const int32_t height,
int64_t *
const M, int64_t *
const H) {
const int32_t wiener_win = WIENER_WIN_CHROMA;
const int32_t wiener_win2 = wiener_win * wiener_win;
const int32_t w16 = width & ~15;
const int32_t h8 = height & ~7;
int16x8_t mask[2];
mask[0] = vld1q_s16(&(mask_16bit[16]) - width % 16);
mask[1] = vld1q_s16(&(mask_16bit[16]) - width % 16 + 8);
const int bit_depth = 8;
int32_t i, j, x, y;
const int32_t num_bit_left =
32 - 1
/* sign */ - 2 * bit_depth /* energy */ + 2 /* SIMD */;
const int32_t h_allowed =
(1 << num_bit_left) / (w16 + ((w16 != width) ? 16 : 0));
// Step 1: Calculate the top edge of the whole matrix, i.e., the top
// edge of each triangle and square on the top row.
j = 0;
do {
const int16_t *s_t = s;
const int16_t *d_t = d;
int32_t height_t = 0;
int64x2_t sum_m[WIENER_WIN_CHROMA] = { vdupq_n_s64(0) };
int64x2_t sum_h[WIENER_WIN_CHROMA] = { vdupq_n_s64(0) };
int16x8_t src[2], dgd[2];
do {
const int32_t h_t =
((height - height_t) < h_allowed) ? (height - height_t) : h_allowed;
int32x4_t row_m[WIENER_WIN_CHROMA] = { vdupq_n_s32(0) };
int32x4_t row_h[WIENER_WIN_CHROMA] = { vdupq_n_s32(0) };
y = h_t;
do {
x = 0;
while (x < w16) {
src[0] = vld1q_s16(s_t + x + 0);
src[1] = vld1q_s16(s_t + x + 8);
dgd[0] = vld1q_s16(d_t + x + 0);
dgd[1] = vld1q_s16(d_t + x + 8);
stats_top_win5_neon(src, dgd, d_t + j + x, d_stride, row_m, row_h);
x += 16;
}
if (w16 != width) {
src[0] = vld1q_s16(s_t + w16 + 0);
src[1] = vld1q_s16(s_t + w16 + 8);
dgd[0] = vld1q_s16(d_t + w16 + 0);
dgd[1] = vld1q_s16(d_t + w16 + 8);
src[0] = vandq_s16(src[0], mask[0]);
src[1] = vandq_s16(src[1], mask[1]);
dgd[0] = vandq_s16(dgd[0], mask[0]);
dgd[1] = vandq_s16(dgd[1], mask[1]);
stats_top_win5_neon(src, dgd, d_t + j + w16, d_stride, row_m, row_h);
}
s_t += s_stride;
d_t += d_stride;
}
while (--y);
sum_m[0] = vpadalq_s32(sum_m[0], row_m[0]);
sum_m[1] = vpadalq_s32(sum_m[1], row_m[1]);
sum_m[2] = vpadalq_s32(sum_m[2], row_m[2]);
sum_m[3] = vpadalq_s32(sum_m[3], row_m[3]);
sum_m[4] = vpadalq_s32(sum_m[4], row_m[4]);
sum_h[0] = vpadalq_s32(sum_h[0], row_h[0]);
sum_h[1] = vpadalq_s32(sum_h[1], row_h[1]);
sum_h[2] = vpadalq_s32(sum_h[2], row_h[2]);
sum_h[3] = vpadalq_s32(sum_h[3], row_h[3]);
sum_h[4] = vpadalq_s32(sum_h[4], row_h[4]);
height_t += h_t;
}
while (height_t < height);
#if AOM_ARCH_AARCH64
int64x2_t sum_m0 = vpaddq_s64(sum_m[0], sum_m[1]);
int64x2_t sum_m2 = vpaddq_s64(sum_m[2], sum_m[3]);
vst1q_s64(&M[wiener_win * j + 0], sum_m0);
vst1q_s64(&M[wiener_win * j + 2], sum_m2);
M[wiener_win * j + 4] = vaddvq_s64(sum_m[4]);
int64x2_t sum_h0 = vpaddq_s64(sum_h[0], sum_h[1]);
int64x2_t sum_h2 = vpaddq_s64(sum_h[2], sum_h[3]);
vst1q_s64(&H[wiener_win * j + 0], sum_h0);
vst1q_s64(&H[wiener_win * j + 2], sum_h2);
H[wiener_win * j + 4] = vaddvq_s64(sum_h[4]);
#else
M[wiener_win * j + 0] = horizontal_add_s64x2(sum_m[0]);
M[wiener_win * j + 1] = horizontal_add_s64x2(sum_m[1]);
M[wiener_win * j + 2] = horizontal_add_s64x2(sum_m[2]);
M[wiener_win * j + 3] = horizontal_add_s64x2(sum_m[3]);
M[wiener_win * j + 4] = horizontal_add_s64x2(sum_m[4]);
H[wiener_win * j + 0] = horizontal_add_s64x2(sum_h[0]);
H[wiener_win * j + 1] = horizontal_add_s64x2(sum_h[1]);
H[wiener_win * j + 2] = horizontal_add_s64x2(sum_h[2]);
H[wiener_win * j + 3] = horizontal_add_s64x2(sum_h[3]);
H[wiener_win * j + 4] = horizontal_add_s64x2(sum_h[4]);
#endif // AOM_ARCH_AARCH64
}
while (++j < wiener_win);
// Step 2: Calculate the left edge of each square on the top row.
j = 1;
do {
const int16_t *d_t = d;
int32_t height_t = 0;
int64x2_t sum_h[WIENER_WIN_CHROMA - 1] = { vdupq_n_s64(0) };
int16x8_t dgd[2];
do {
const int32_t h_t =
((height - height_t) < h_allowed) ? (height - height_t) : h_allowed;
int32x4_t row_h[WIENER_WIN_CHROMA - 1] = { vdupq_n_s32(0) };
y = h_t;
do {
x = 0;
while (x < w16) {
dgd[0] = vld1q_s16(d_t + j + x + 0);
dgd[1] = vld1q_s16(d_t + j + x + 8);
stats_left_win5_neon(dgd, d_t + x, d_stride, row_h);
x += 16;
}
if (w16 != width) {
dgd[0] = vld1q_s16(d_t + j + x + 0);
dgd[1] = vld1q_s16(d_t + j + x + 8);
dgd[0] = vandq_s16(dgd[0], mask[0]);
dgd[1] = vandq_s16(dgd[1], mask[1]);
stats_left_win5_neon(dgd, d_t + x, d_stride, row_h);
}
d_t += d_stride;
}
while (--y);
sum_h[0] = vpadalq_s32(sum_h[0], row_h[0]);
sum_h[1] = vpadalq_s32(sum_h[1], row_h[1]);
sum_h[2] = vpadalq_s32(sum_h[2], row_h[2]);
sum_h[3] = vpadalq_s32(sum_h[3], row_h[3]);
height_t += h_t;
}
while (height_t < height);
#if AOM_ARCH_AARCH64
int64x2_t sum_h0 = vpaddq_s64(sum_h[0], sum_h[1]);
int64x2_t sum_h1 = vpaddq_s64(sum_h[2], sum_h[3]);
vst1_s64(&H[1 * wiener_win2 + j * wiener_win], vget_low_s64(sum_h0));
vst1_s64(&H[2 * wiener_win2 + j * wiener_win], vget_high_s64(sum_h0));
vst1_s64(&H[3 * wiener_win2 + j * wiener_win], vget_low_s64(sum_h1));
vst1_s64(&H[4 * wiener_win2 + j * wiener_win], vget_high_s64(sum_h1));
#else
H[1 * wiener_win2 + j * wiener_win] = horizontal_add_s64x2(sum_h[0]);
H[2 * wiener_win2 + j * wiener_win] = horizontal_add_s64x2(sum_h[1]);
H[3 * wiener_win2 + j * wiener_win] = horizontal_add_s64x2(sum_h[2]);
H[4 * wiener_win2 + j * wiener_win] = horizontal_add_s64x2(sum_h[3]);
#endif // AOM_ARCH_AARCH64
}
while (++j < wiener_win);
// Step 3: Derive the top edge of each triangle along the diagonal. No
// triangle in top row.
{
const int16_t *d_t = d;
if (height % 2) {
int32x4_t deltas[(WIENER_WIN + 1) * 2] = { vdupq_n_s32(0) };
int32x4_t deltas_tr[(WIENER_WIN + 1) * 2] = { vdupq_n_s32(0) };
int16x8_t ds[WIENER_WIN * 2];
load_s16_8x4(d_t, d_stride, &ds[0], &ds[2], &ds[4], &ds[6]);
load_s16_8x4(d_t + width, d_stride, &ds[1], &ds[3], &ds[5], &ds[7]);
d_t += 4 * d_stride;
step3_win5_oneline_neon(&d_t, d_stride, width, height, ds, deltas);
transpose_arrays_s32_8x8(deltas, deltas_tr);
update_5_stats_neon(H + 0 * wiener_win * wiener_win2 + 0 * wiener_win,
deltas_tr[0], vgetq_lane_s32(deltas_tr[4], 0),
H + 1 * wiener_win * wiener_win2 + 1 * wiener_win);
update_5_stats_neon(H + 1 * wiener_win * wiener_win2 + 1 * wiener_win,
deltas_tr[1], vgetq_lane_s32(deltas_tr[5], 0),
H + 2 * wiener_win * wiener_win2 + 2 * wiener_win);
update_5_stats_neon(H + 2 * wiener_win * wiener_win2 + 2 * wiener_win,
deltas_tr[2], vgetq_lane_s32(deltas_tr[6], 0),
H + 3 * wiener_win * wiener_win2 + 3 * wiener_win);
update_5_stats_neon(H + 3 * wiener_win * wiener_win2 + 3 * wiener_win,
deltas_tr[3], vgetq_lane_s32(deltas_tr[7], 0),
H + 4 * wiener_win * wiener_win2 + 4 * wiener_win);
}
else {
int32x4_t deltas[WIENER_WIN_CHROMA * 2] = { vdupq_n_s32(0) };
int16x8_t ds[WIENER_WIN_CHROMA * 2];
ds[0] = load_unaligned_s16_4x2(d_t + 0 * d_stride, width);
ds[1] = load_unaligned_s16_4x2(d_t + 1 * d_stride, width);
ds[2] = load_unaligned_s16_4x2(d_t + 2 * d_stride, width);
ds[3] = load_unaligned_s16_4x2(d_t + 3 * d_stride, width);
step3_win5_neon(d_t + 4 * d_stride, d_stride, width, height, ds, deltas);
transpose_elems_inplace_s32_4x4(&deltas[0], &deltas[1], &deltas[2],
&deltas[3]);
update_5_stats_neon(H + 0 * wiener_win * wiener_win2 + 0 * wiener_win,
deltas[0], vgetq_lane_s32(deltas[4], 0),
H + 1 * wiener_win * wiener_win2 + 1 * wiener_win);
update_5_stats_neon(H + 1 * wiener_win * wiener_win2 + 1 * wiener_win,
deltas[1], vgetq_lane_s32(deltas[4], 1),
H + 2 * wiener_win * wiener_win2 + 2 * wiener_win);
update_5_stats_neon(H + 2 * wiener_win * wiener_win2 + 2 * wiener_win,
deltas[2], vgetq_lane_s32(deltas[4], 2),
H + 3 * wiener_win * wiener_win2 + 3 * wiener_win);
update_5_stats_neon(H + 3 * wiener_win * wiener_win2 + 3 * wiener_win,
deltas[3], vgetq_lane_s32(deltas[4], 3),
H + 4 * wiener_win * wiener_win2 + 4 * wiener_win);
}
}
// Step 4: Derive the top and left edge of each square. No square in top and
// bottom row.
{
y = h8;
int16x4_t d_s[12];
int16x4_t d_e[12];
const int16_t *d_t = d;
int16x4_t zeros = vdup_n_s16(0);
load_s16_4x4(d_t, d_stride, &d_s[0], &d_s[1], &d_s[2], &d_s[3]);
load_s16_4x4(d_t + width, d_stride, &d_e[0], &d_e[1], &d_e[2], &d_e[3]);
int32x4_t deltas[6][18] = { { vdupq_n_s32(0) }, { vdupq_n_s32(0) } };
while (y >= 8) {
load_s16_4x8(d_t + 4 * d_stride, d_stride, &d_s[4], &d_s[5], &d_s[6],
&d_s[7], &d_s[8], &d_s[9], &d_s[10], &d_s[11]);
load_s16_4x8(d_t + width + 4 * d_stride, d_stride, &d_e[4], &d_e[5],
&d_e[6], &d_e[7], &d_e[8], &d_e[9], &d_e[10], &d_e[11]);
int16x8_t s_tr[8], e_tr[8];
transpose_elems_s16_4x8(d_s[0], d_s[1], d_s[2], d_s[3], d_s[4], d_s[5],
d_s[6], d_s[7], &s_tr[0], &s_tr[1], &s_tr[2],
&s_tr[3]);
transpose_elems_s16_4x8(d_s[8], d_s[9], d_s[10], d_s[11], zeros, zeros,
zeros, zeros, &s_tr[4], &s_tr[5], &s_tr[6],
&s_tr[7]);
transpose_elems_s16_4x8(d_e[0], d_e[1], d_e[2], d_e[3], d_e[4], d_e[5],
d_e[6], d_e[7], &e_tr[0], &e_tr[1], &e_tr[2],
&e_tr[3]);
transpose_elems_s16_4x8(d_e[8], d_e[9], d_e[10], d_e[11], zeros, zeros,
zeros, zeros, &e_tr[4], &e_tr[5], &e_tr[6],
&e_tr[7]);
int16x8_t start_col0[5], start_col1[5], start_col2[5], start_col3[5];
start_col0[0] = s_tr[0];
start_col0[1] = vextq_s16(s_tr[0], s_tr[4], 1);
start_col0[2] = vextq_s16(s_tr[0], s_tr[4], 2);
start_col0[3] = vextq_s16(s_tr[0], s_tr[4], 3);
start_col0[4] = vextq_s16(s_tr[0], s_tr[4], 4);
start_col1[0] = s_tr[1];
start_col1[1] = vextq_s16(s_tr[1], s_tr[5], 1);
start_col1[2] = vextq_s16(s_tr[1], s_tr[5], 2);
start_col1[3] = vextq_s16(s_tr[1], s_tr[5], 3);
start_col1[4] = vextq_s16(s_tr[1], s_tr[5], 4);
start_col2[0] = s_tr[2];
start_col2[1] = vextq_s16(s_tr[2], s_tr[6], 1);
start_col2[2] = vextq_s16(s_tr[2], s_tr[6], 2);
start_col2[3] = vextq_s16(s_tr[2], s_tr[6], 3);
start_col2[4] = vextq_s16(s_tr[2], s_tr[6], 4);
start_col3[0] = s_tr[3];
start_col3[1] = vextq_s16(s_tr[3], s_tr[7], 1);
start_col3[2] = vextq_s16(s_tr[3], s_tr[7], 2);
start_col3[3] = vextq_s16(s_tr[3], s_tr[7], 3);
start_col3[4] = vextq_s16(s_tr[3], s_tr[7], 4);
// i = 1, j = 2;
sub_deltas_step4(start_col0, start_col1, deltas[0]);
// i = 1, j = 3;
sub_deltas_step4(start_col0, start_col2, deltas[1]);
// i = 1, j = 4
sub_deltas_step4(start_col0, start_col3, deltas[2]);
// i = 2, j =3
sub_deltas_step4(start_col1, start_col2, deltas[3]);
// i = 2, j = 4
sub_deltas_step4(start_col1, start_col3, deltas[4]);
// i = 3, j = 4
sub_deltas_step4(start_col2, start_col3, deltas[5]);
int16x8_t end_col0[5], end_col1[5], end_col2[5], end_col3[5];
end_col0[0] = e_tr[0];
end_col0[1] = vextq_s16(e_tr[0], e_tr[4], 1);
end_col0[2] = vextq_s16(e_tr[0], e_tr[4], 2);
end_col0[3] = vextq_s16(e_tr[0], e_tr[4], 3);
end_col0[4] = vextq_s16(e_tr[0], e_tr[4], 4);
end_col1[0] = e_tr[1];
end_col1[1] = vextq_s16(e_tr[1], e_tr[5], 1);
end_col1[2] = vextq_s16(e_tr[1], e_tr[5], 2);
end_col1[3] = vextq_s16(e_tr[1], e_tr[5], 3);
end_col1[4] = vextq_s16(e_tr[1], e_tr[5], 4);
end_col2[0] = e_tr[2];
end_col2[1] = vextq_s16(e_tr[2], e_tr[6], 1);
end_col2[2] = vextq_s16(e_tr[2], e_tr[6], 2);
end_col2[3] = vextq_s16(e_tr[2], e_tr[6], 3);
end_col2[4] = vextq_s16(e_tr[2], e_tr[6], 4);
end_col3[0] = e_tr[3];
end_col3[1] = vextq_s16(e_tr[3], e_tr[7], 1);
end_col3[2] = vextq_s16(e_tr[3], e_tr[7], 2);
end_col3[3] = vextq_s16(e_tr[3], e_tr[7], 3);
end_col3[4] = vextq_s16(e_tr[3], e_tr[7], 4);
// i = 1, j = 2;
add_deltas_step4(end_col0, end_col1, deltas[0]);
// i = 1, j = 3;
add_deltas_step4(end_col0, end_col2, deltas[1]);
// i = 1, j = 4
add_deltas_step4(end_col0, end_col3, deltas[2]);
// i = 2, j =3
add_deltas_step4(end_col1, end_col2, deltas[3]);
// i = 2, j = 4
add_deltas_step4(end_col1, end_col3, deltas[4]);
// i = 3, j = 4
add_deltas_step4(end_col2, end_col3, deltas[5]);
d_s[0] = d_s[8];
d_s[1] = d_s[9];
d_s[2] = d_s[10];
d_s[3] = d_s[11];
d_e[0] = d_e[8];
d_e[1] = d_e[9];
d_e[2] = d_e[10];
d_e[3] = d_e[11];
d_t += 8 * d_stride;
y -= 8;
}
if (h8 != height) {
const int16x8_t mask_h = vld1q_s16(&mask_16bit[16] - (height % 8));
load_s16_4x8(d_t + 4 * d_stride, d_stride, &d_s[4], &d_s[5], &d_s[6],
&d_s[7], &d_s[8], &d_s[9], &d_s[10], &d_s[11]);
load_s16_4x8(d_t + width + 4 * d_stride, d_stride, &d_e[4], &d_e[5],
&d_e[6], &d_e[7], &d_e[8], &d_e[9], &d_e[10], &d_e[11]);
int16x8_t s_tr[8], e_tr[8];
transpose_elems_s16_4x8(d_s[0], d_s[1], d_s[2], d_s[3], d_s[4], d_s[5],
d_s[6], d_s[7], &s_tr[0], &s_tr[1], &s_tr[2],
&s_tr[3]);
transpose_elems_s16_4x8(d_s[8], d_s[9], d_s[10], d_s[11], zeros, zeros,
zeros, zeros, &s_tr[4], &s_tr[5], &s_tr[6],
&s_tr[7]);
transpose_elems_s16_4x8(d_e[0], d_e[1], d_e[2], d_e[3], d_e[4], d_e[5],
d_e[6], d_e[7], &e_tr[0], &e_tr[1], &e_tr[2],
&e_tr[3]);
transpose_elems_s16_4x8(d_e[8], d_e[9], d_e[10], d_e[11], zeros, zeros,
zeros, zeros, &e_tr[4], &e_tr[5], &e_tr[6],
&e_tr[7]);
int16x8_t start_col0[5], start_col1[5], start_col2[5], start_col3[5];
start_col0[0] = vandq_s16(s_tr[0], mask_h);
start_col0[1] = vandq_s16(vextq_s16(s_tr[0], s_tr[4], 1), mask_h);
start_col0[2] = vandq_s16(vextq_s16(s_tr[0], s_tr[4], 2), mask_h);
start_col0[3] = vandq_s16(vextq_s16(s_tr[0], s_tr[4], 3), mask_h);
start_col0[4] = vandq_s16(vextq_s16(s_tr[0], s_tr[4], 4), mask_h);
start_col1[0] = vandq_s16(s_tr[1], mask_h);
start_col1[1] = vandq_s16(vextq_s16(s_tr[1], s_tr[5], 1), mask_h);
start_col1[2] = vandq_s16(vextq_s16(s_tr[1], s_tr[5], 2), mask_h);
start_col1[3] = vandq_s16(vextq_s16(s_tr[1], s_tr[5], 3), mask_h);
start_col1[4] = vandq_s16(vextq_s16(s_tr[1], s_tr[5], 4), mask_h);
start_col2[0] = vandq_s16(s_tr[2], mask_h);
start_col2[1] = vandq_s16(vextq_s16(s_tr[2], s_tr[6], 1), mask_h);
start_col2[2] = vandq_s16(vextq_s16(s_tr[2], s_tr[6], 2), mask_h);
start_col2[3] = vandq_s16(vextq_s16(s_tr[2], s_tr[6], 3), mask_h);
start_col2[4] = vandq_s16(vextq_s16(s_tr[2], s_tr[6], 4), mask_h);
start_col3[0] = vandq_s16(s_tr[3], mask_h);
start_col3[1] = vandq_s16(vextq_s16(s_tr[3], s_tr[7], 1), mask_h);
start_col3[2] = vandq_s16(vextq_s16(s_tr[3], s_tr[7], 2), mask_h);
start_col3[3] = vandq_s16(vextq_s16(s_tr[3], s_tr[7], 3), mask_h);
start_col3[4] = vandq_s16(vextq_s16(s_tr[3], s_tr[7], 4), mask_h);
// i = 1, j = 2;
sub_deltas_step4(start_col0, start_col1, deltas[0]);
// i = 1, j = 3;
sub_deltas_step4(start_col0, start_col2, deltas[1]);
// i = 1, j = 4
sub_deltas_step4(start_col0, start_col3, deltas[2]);
// i = 2, j = 3
sub_deltas_step4(start_col1, start_col2, deltas[3]);
// i = 2, j = 4
sub_deltas_step4(start_col1, start_col3, deltas[4]);
// i = 3, j = 4
sub_deltas_step4(start_col2, start_col3, deltas[5]);
int16x8_t end_col0[5], end_col1[5], end_col2[5], end_col3[5];
end_col0[0] = vandq_s16(e_tr[0], mask_h);
end_col0[1] = vandq_s16(vextq_s16(e_tr[0], e_tr[4], 1), mask_h);
end_col0[2] = vandq_s16(vextq_s16(e_tr[0], e_tr[4], 2), mask_h);
end_col0[3] = vandq_s16(vextq_s16(e_tr[0], e_tr[4], 3), mask_h);
end_col0[4] = vandq_s16(vextq_s16(e_tr[0], e_tr[4], 4), mask_h);
end_col1[0] = vandq_s16(e_tr[1], mask_h);
end_col1[1] = vandq_s16(vextq_s16(e_tr[1], e_tr[5], 1), mask_h);
end_col1[2] = vandq_s16(vextq_s16(e_tr[1], e_tr[5], 2), mask_h);
end_col1[3] = vandq_s16(vextq_s16(e_tr[1], e_tr[5], 3), mask_h);
end_col1[4] = vandq_s16(vextq_s16(e_tr[1], e_tr[5], 4), mask_h);
end_col2[0] = vandq_s16(e_tr[2], mask_h);
end_col2[1] = vandq_s16(vextq_s16(e_tr[2], e_tr[6], 1), mask_h);
end_col2[2] = vandq_s16(vextq_s16(e_tr[2], e_tr[6], 2), mask_h);
end_col2[3] = vandq_s16(vextq_s16(e_tr[2], e_tr[6], 3), mask_h);
end_col2[4] = vandq_s16(vextq_s16(e_tr[2], e_tr[6], 4), mask_h);
end_col3[0] = vandq_s16(e_tr[3], mask_h);
end_col3[1] = vandq_s16(vextq_s16(e_tr[3], e_tr[7], 1), mask_h);
end_col3[2] = vandq_s16(vextq_s16(e_tr[3], e_tr[7], 2), mask_h);
end_col3[3] = vandq_s16(vextq_s16(e_tr[3], e_tr[7], 3), mask_h);
end_col3[4] = vandq_s16(vextq_s16(e_tr[3], e_tr[7], 4), mask_h);
// i = 1, j = 2;
add_deltas_step4(end_col0, end_col1, deltas[0]);
// i = 1, j = 3;
add_deltas_step4(end_col0, end_col2, deltas[1]);
// i = 1, j = 4
add_deltas_step4(end_col0, end_col3, deltas[2]);
// i = 2, j =3
add_deltas_step4(end_col1, end_col2, deltas[3]);
// i = 2, j = 4
add_deltas_step4(end_col1, end_col3, deltas[4]);
// i = 3, j = 4
add_deltas_step4(end_col2, end_col3, deltas[5]);
}
int32x4_t delta[6][2];
int32_t single_delta[6];
delta[0][0] = horizontal_add_4d_s32x4(&deltas[0][0]);
delta[1][0] = horizontal_add_4d_s32x4(&deltas[1][0]);
delta[2][0] = horizontal_add_4d_s32x4(&deltas[2][0]);
delta[3][0] = horizontal_add_4d_s32x4(&deltas[3][0]);
delta[4][0] = horizontal_add_4d_s32x4(&deltas[4][0]);
delta[5][0] = horizontal_add_4d_s32x4(&deltas[5][0]);
delta[0][1] = horizontal_add_4d_s32x4(&deltas[0][5]);
delta[1][1] = horizontal_add_4d_s32x4(&deltas[1][5]);
delta[2][1] = horizontal_add_4d_s32x4(&deltas[2][5]);
delta[3][1] = horizontal_add_4d_s32x4(&deltas[3][5]);
delta[4][1] = horizontal_add_4d_s32x4(&deltas[4][5]);
delta[5][1] = horizontal_add_4d_s32x4(&deltas[5][5]);
single_delta[0] = horizontal_add_s32x4(deltas[0][4]);
single_delta[1] = horizontal_add_s32x4(deltas[1][4]);
single_delta[2] = horizontal_add_s32x4(deltas[2][4]);
single_delta[3] = horizontal_add_s32x4(deltas[3][4]);
single_delta[4] = horizontal_add_s32x4(deltas[4][4]);
single_delta[5] = horizontal_add_s32x4(deltas[5][4]);
int idx = 0;
for (i = 1; i < wiener_win - 1; i++) {
for (j = i + 1; j < wiener_win; j++) {
update_4_stats_neon(
H + (i - 1) * wiener_win * wiener_win2 + (j - 1) * wiener_win,
delta[idx][0], H + i * wiener_win * wiener_win2 + j * wiener_win);
H[i * wiener_win * wiener_win2 + j * wiener_win + 4] =
H[(i - 1) * wiener_win * wiener_win2 + (j - 1) * wiener_win + 4] +
single_delta[idx];
H[(i * wiener_win + 1) * wiener_win2 + j * wiener_win] =
H[((i - 1) * wiener_win + 1) * wiener_win2 + (j - 1) * wiener_win] +
vgetq_lane_s32(delta[idx][1], 0);
H[(i * wiener_win + 2) * wiener_win2 + j * wiener_win] =
H[((i - 1) * wiener_win + 2) * wiener_win2 + (j - 1) * wiener_win] +
vgetq_lane_s32(delta[idx][1], 1);
H[(i * wiener_win + 3) * wiener_win2 + j * wiener_win] =
H[((i - 1) * wiener_win + 3) * wiener_win2 + (j - 1) * wiener_win] +
vgetq_lane_s32(delta[idx][1], 2);
H[(i * wiener_win + 4) * wiener_win2 + j * wiener_win] =
H[((i - 1) * wiener_win + 4) * wiener_win2 + (j - 1) * wiener_win] +
vgetq_lane_s32(delta[idx][1], 3);
idx++;
}
}
}
// Step 5: Derive other points of each square. No square in bottom row.
i = 0;
do {
const int16_t *
const di = d + i;
j = i + 1;
do {
const int16_t *
const dj = d + j;
int32x4_t deltas[WIENER_WIN_CHROMA - 1][WIENER_WIN_CHROMA - 1] = {
{ vdupq_n_s32(0) }, { vdupq_n_s32(0) }
};
int16x8_t d_is[WIN_CHROMA], d_ie[WIN_CHROMA];
int16x8_t d_js[WIN_CHROMA], d_je[WIN_CHROMA];
x = 0;
while (x < w16) {
load_square_win5_neon(di + x, dj + x, d_stride, height, d_is, d_ie,
d_js, d_je);
derive_square_win5_neon(d_is, d_ie, d_js, d_je, deltas);
x += 16;
}
if (w16 != width) {
load_square_win5_neon(di + x, dj + x, d_stride, height, d_is, d_ie,
d_js, d_je);
d_is[0] = vandq_s16(d_is[0], mask[0]);
d_is[1] = vandq_s16(d_is[1], mask[1]);
d_is[2] = vandq_s16(d_is[2], mask[0]);
d_is[3] = vandq_s16(d_is[3], mask[1]);
d_is[4] = vandq_s16(d_is[4], mask[0]);
d_is[5] = vandq_s16(d_is[5], mask[1]);
d_is[6] = vandq_s16(d_is[6], mask[0]);
d_is[7] = vandq_s16(d_is[7], mask[1]);
d_ie[0] = vandq_s16(d_ie[0], mask[0]);
d_ie[1] = vandq_s16(d_ie[1], mask[1]);
d_ie[2] = vandq_s16(d_ie[2], mask[0]);
d_ie[3] = vandq_s16(d_ie[3], mask[1]);
d_ie[4] = vandq_s16(d_ie[4], mask[0]);
d_ie[5] = vandq_s16(d_ie[5], mask[1]);
d_ie[6] = vandq_s16(d_ie[6], mask[0]);
d_ie[7] = vandq_s16(d_ie[7], mask[1]);
derive_square_win5_neon(d_is, d_ie, d_js, d_je, deltas);
}
hadd_update_4_stats_neon(
H + (i * wiener_win + 0) * wiener_win2 + j * wiener_win, deltas[0],
H + (i * wiener_win + 1) * wiener_win2 + j * wiener_win + 1);
hadd_update_4_stats_neon(
H + (i * wiener_win + 1) * wiener_win2 + j * wiener_win, deltas[1],
H + (i * wiener_win + 2) * wiener_win2 + j * wiener_win + 1);
hadd_update_4_stats_neon(
H + (i * wiener_win + 2) * wiener_win2 + j * wiener_win, deltas[2],
H + (i * wiener_win + 3) * wiener_win2 + j * wiener_win + 1);
hadd_update_4_stats_neon(
H + (i * wiener_win + 3) * wiener_win2 + j * wiener_win, deltas[3],
H + (i * wiener_win + 4) * wiener_win2 + j * wiener_win + 1);
}
while (++j < wiener_win);
}
while (++i < wiener_win - 1);
// Step 6: Derive other points of each upper triangle along the diagonal.
i = 0;
do {
const int16_t *
const di = d + i;
int32x4_t deltas[WIENER_WIN_CHROMA * 2 + 1] = { vdupq_n_s32(0) };
int16x8_t d_is[WIN_CHROMA], d_ie[WIN_CHROMA];
x = 0;
while (x < w16) {
load_triangle_win5_neon(di + x, d_stride, height, d_is, d_ie);
derive_triangle_win5_neon(d_is, d_ie, deltas);
x += 16;
}
if (w16 != width) {
load_triangle_win5_neon(di + x, d_stride, height, d_is, d_ie);
d_is[0] = vandq_s16(d_is[0], mask[0]);
d_is[1] = vandq_s16(d_is[1], mask[1]);
d_is[2] = vandq_s16(d_is[2], mask[0]);
d_is[3] = vandq_s16(d_is[3], mask[1]);
d_is[4] = vandq_s16(d_is[4], mask[0]);
d_is[5] = vandq_s16(d_is[5], mask[1]);
d_is[6] = vandq_s16(d_is[6], mask[0]);
d_is[7] = vandq_s16(d_is[7], mask[1]);
d_ie[0] = vandq_s16(d_ie[0], mask[0]);
d_ie[1] = vandq_s16(d_ie[1], mask[1]);
d_ie[2] = vandq_s16(d_ie[2], mask[0]);
d_ie[3] = vandq_s16(d_ie[3], mask[1]);
d_ie[4] = vandq_s16(d_ie[4], mask[0]);
d_ie[5] = vandq_s16(d_ie[5], mask[1]);
d_ie[6] = vandq_s16(d_ie[6], mask[0]);
d_ie[7] = vandq_s16(d_ie[7], mask[1]);
derive_triangle_win5_neon(d_is, d_ie, deltas);
}
// Row 1: 4 points
hadd_update_4_stats_neon(
H + (i * wiener_win + 0) * wiener_win2 + i * wiener_win, deltas,
H + (i * wiener_win + 1) * wiener_win2 + i * wiener_win + 1);
// Row 2: 3 points
int32x4_t deltas45 = horizontal_add_2d_s32(deltas[4], deltas[5]);
int32x4_t deltas78 = horizontal_add_2d_s32(deltas[7], deltas[8]);
int64x2_t deltas45_s64 = vpaddlq_s32(deltas45);
int64x2_t deltas78_s64 = vpaddlq_s32(deltas78);
int64x2_t src =
vld1q_s64(H + (i * wiener_win + 1) * wiener_win2 + i * wiener_win + 1);
int64x2_t dst = vaddq_s64(src, deltas45_s64);
vst1q_s64(H + (i * wiener_win + 2) * wiener_win2 + i * wiener_win + 2, dst);
int32x4_t delta69 = horizontal_add_2d_s32(deltas[6], deltas[9]);
int64x2_t delta69_s64 = vpaddlq_s32(delta69);
H[(i * wiener_win + 2) * wiener_win2 + i * wiener_win + 4] =
H[(i * wiener_win + 1) * wiener_win2 + i * wiener_win + 3] +
vgetq_lane_s64(delta69_s64, 0);
// Row 3: 2 points
vst1q_s64(H + (i * wiener_win + 3) * wiener_win2 + i * wiener_win + 3,
vaddq_s64(dst, deltas78_s64));
// Row 4: 1 point
H[(i * wiener_win + 4) * wiener_win2 + i * wiener_win + 4] =
H[(i * wiener_win + 3) * wiener_win2 + i * wiener_win + 3] +
vgetq_lane_s64(delta69_s64, 1);
}
while (++i < wiener_win);
}
static inline void compute_stats_win7_neon(
const int16_t *
const d,
const int32_t d_stride,
const int16_t *
const s,
const int32_t s_stride,
const int32_t width,
const int32_t height,
int64_t *
const M, int64_t *
const H) {
const int32_t wiener_win = WIENER_WIN;
const int32_t wiener_win2 = wiener_win * wiener_win;
const int32_t w16 = width & ~15;
const int32_t h8 = height & ~7;
int16x8_t mask[2];
mask[0] = vld1q_s16(&(mask_16bit[16]) - width % 16);
mask[1] = vld1q_s16(&(mask_16bit[16]) - width % 16 + 8);
const int bit_depth = 8;
int32_t i, j, x, y;
const int32_t num_bit_left =
32 - 1
/* sign */ - 2 * bit_depth /* energy */ + 2 /* SIMD */;
const int32_t h_allowed =
(1 << num_bit_left) / (w16 + ((w16 != width) ? 16 : 0));
// Step 1: Calculate the top edge of the whole matrix, i.e., the top
// edge of each triangle and square on the top row.
j = 0;
do {
const int16_t *s_t = s;
const int16_t *d_t = d;
int32_t height_t = 0;
int64x2_t sum_m[WIENER_WIN] = { vdupq_n_s64(0) };
int64x2_t sum_h[WIENER_WIN] = { vdupq_n_s64(0) };
int16x8_t src[2], dgd[2];
do {
const int32_t h_t =
((height - height_t) < h_allowed) ? (height - height_t) : h_allowed;
int32x4_t row_m[WIENER_WIN * 2] = { vdupq_n_s32(0) };
int32x4_t row_h[WIENER_WIN * 2] = { vdupq_n_s32(0) };
y = h_t;
do {
x = 0;
while (x < w16) {
src[0] = vld1q_s16(s_t + x);
src[1] = vld1q_s16(s_t + x + 8);
dgd[0] = vld1q_s16(d_t + x);
dgd[1] = vld1q_s16(d_t + x + 8);
stats_top_win7_neon(src, dgd, d_t + j + x, d_stride, row_m, row_h);
x += 16;
}
if (w16 != width) {
src[0] = vld1q_s16(s_t + w16);
--> --------------------
--> maximum size reached
--> --------------------
