Fused horizontal separable filter over N planes.
No aliasing between any src pointer and any dst pointer.
204{
205 constexpr size_t Unroll = 4;
206
207 const auto half = static_cast<int32_t>(kernel.size() / 2);
208 const auto iw = static_cast<int32_t>(w);
209 const float* kdata = kernel.data();
210 const size_t n =
src.size();
211
212 std::for_each(std::execution::par_unseq,
213 std::views::iota(uint32_t { 0 },
h).begin(),
214 std::views::iota(uint32_t { 0 },
h).end(),
215 [&](uint32_t row) {
216 const size_t row_off = static_cast<size_t>(row) * w;
217
218 const float* s_ptr[32];
219 float* d_ptr[32];
220 for (size_t i = 0; i < n; ++i) {
221 s_ptr[i] =
src[i] + row_off;
222 d_ptr[i] =
dst[i] + row_off;
223 }
224
225
226 for (size_t i = 0; i < n; ++i) {
227 for (int32_t
px = 0;
px < std::min(half, iw); ++
px) {
228 float acc = 0.0F;
229 for (int32_t k = -half; k <= half; ++k) {
230 acc += s_ptr[i][std::clamp(
px + k, 0, iw - 1)]
231 * kdata[static_cast<size_t>(k + half)];
232 }
234 }
235 }
236
237 const int32_t interior_end = iw - half;
238
239#ifdef MAYAFLUX_ARCH_X64
241 for (;
px + 8 <= interior_end;
px += 8) {
242 size_t i = 0;
243 for (; i + Unroll - 1 < n; i += Unroll) {
244 __m256 acc0 = _mm256_setzero_ps();
245 __m256 acc1 = _mm256_setzero_ps();
246 __m256 acc2 = _mm256_setzero_ps();
247 __m256 acc3 = _mm256_setzero_ps();
248 for (int32_t k = -half; k <= half; ++k) {
249 const __m256 kv = _mm256_set1_ps(kdata[static_cast<size_t>(k + half)]);
250 acc0 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 0] +
px + k), kv, acc0);
251 acc1 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 1] +
px + k), kv, acc1);
252 acc2 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 2] +
px + k), kv, acc2);
253 acc3 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 3] +
px + k), kv, acc3);
254 }
255 _mm256_storeu_ps(d_ptr[i + 0] +
px, acc0);
256 _mm256_storeu_ps(d_ptr[i + 1] +
px, acc1);
257 _mm256_storeu_ps(d_ptr[i + 2] +
px, acc2);
258 _mm256_storeu_ps(d_ptr[i + 3] +
px, acc3);
259 }
260 for (; i < n; ++i) {
261 __m256 acc = _mm256_setzero_ps();
262 for (int32_t k = -half; k <= half; ++k) {
263 acc = _mm256_fmadd_ps(
264 _mm256_loadu_ps(s_ptr[i] +
px + k),
265 _mm256_set1_ps(kdata[static_cast<size_t>(k + half)]),
266 acc);
267 }
268 _mm256_storeu_ps(d_ptr[i] +
px, acc);
269 }
270 }
271
272 for (size_t i = 0; i < n; ++i) {
273 for (int32_t px2 =
px; px2 < iw; ++px2) {
274 float acc = 0.0F;
275 for (int32_t k = -half; k <= half; ++k) {
276 acc += s_ptr[i][std::clamp(px2 + k, 0, iw - 1)]
277 * kdata[static_cast<size_t>(k + half)];
278 }
279 d_ptr[i][px2] = acc;
280 }
281 }
282
283#elif defined(MAYAFLUX_ARCH_ARM64)
285 for (;
px + 4 <= interior_end;
px += 4) {
286 size_t i = 0;
287 for (; i + Unroll - 1 < n; i += Unroll) {
288 float32x4_t acc0 = vdupq_n_f32(0.0F);
289 float32x4_t acc1 = vdupq_n_f32(0.0F);
290 float32x4_t acc2 = vdupq_n_f32(0.0F);
291 float32x4_t acc3 = vdupq_n_f32(0.0F);
292 for (int32_t k = -half; k <= half; ++k) {
293 const float32x4_t kv = vdupq_n_f32(kdata[static_cast<size_t>(k + half)]);
294 acc0 = vmlaq_f32(acc0, vld1q_f32(s_ptr[i + 0] +
px + k), kv);
295 acc1 = vmlaq_f32(acc1, vld1q_f32(s_ptr[i + 1] +
px + k), kv);
296 acc2 = vmlaq_f32(acc2, vld1q_f32(s_ptr[i + 2] +
px + k), kv);
297 acc3 = vmlaq_f32(acc3, vld1q_f32(s_ptr[i + 3] +
px + k), kv);
298 }
299 vst1q_f32(d_ptr[i + 0] +
px, acc0);
300 vst1q_f32(d_ptr[i + 1] +
px, acc1);
301 vst1q_f32(d_ptr[i + 2] +
px, acc2);
302 vst1q_f32(d_ptr[i + 3] +
px, acc3);
303 }
304 for (; i < n; ++i) {
305 float32x4_t acc = vdupq_n_f32(0.0F);
306 for (int32_t k = -half; k <= half; ++k)
307 acc = vmlaq_f32(acc,
308 vld1q_f32(s_ptr[i] +
px + k),
309 vdupq_n_f32(kdata[static_cast<size_t>(k + half)]));
310 vst1q_f32(d_ptr[i] +
px, acc);
311 }
312 }
313 for (size_t i = 0; i < n; ++i) {
314 for (int32_t px2 =
px; px2 < iw; ++px2) {
315 float acc = 0.0F;
316 for (int32_t k = -half; k <= half; ++k)
317 acc += s_ptr[i][std::clamp(px2 + k, 0, iw - 1)]
318 * kdata[static_cast<size_t>(k + half)];
319 d_ptr[i][px2] = acc;
320 }
321 }
322
323#else
324 for (size_t i = 0; i < n; ++i) {
325 for (int32_t
px = half;
px < interior_end; ++
px) {
326 float acc = 0.0F;
327 for (int32_t k = -half; k <= half; ++k)
328 acc += s_ptr[i][
px + k] * kdata[
static_cast<size_t>(k + half)];
330 }
331 for (int32_t
px = interior_end;
px < iw; ++
px) {
332 float acc = 0.0F;
333 for (int32_t k = -half; k <= half; ++k)
334 acc += s_ptr[i][std::clamp(
px + k, 0, iw - 1)]
335 * kdata[static_cast<size_t>(k + half)];
337 }
338 }
339#endif
340 });
341}