Fused vertical separable filter over N planes.
Same contract as filter_horizontal_planes but operates vertically.
350{
351 constexpr size_t Unroll = 4;
352
353 const auto half = static_cast<int32_t>(kernel.size() / 2);
354 const auto ih =
static_cast<int32_t
>(
h);
355 const auto iw = static_cast<int32_t>(w);
356 const float* kdata = kernel.data();
357 const size_t n = src.size();
358
359 P::for_each(P::par_unseq,
360 std::views::iota(uint32_t { 0 },
h).begin(),
361 std::views::iota(uint32_t { 0 },
h).end(),
362 [&](uint32_t row) {
363 const auto py = static_cast<int32_t>(row);
364 const bool border = (py < half || py >= ih - half);
365 const size_t row_off = static_cast<size_t>(row) * w;
366
367 const float* s_ptr[32];
368 float* d_ptr[32];
369 for (size_t i = 0; i < n; ++i) {
370 s_ptr[i] = src[i];
371 d_ptr[i] = dst[i] + row_off;
372 }
373
374#ifdef MAYAFLUX_ARCH_X64
375 int32_t px = 0;
376 if (!border) {
377 for (; px + 8 <= iw; px += 8) {
378 size_t i = 0;
379 for (; i + Unroll - 1 < n; i += Unroll) {
380 __m256 acc0 = _mm256_setzero_ps();
381 __m256 acc1 = _mm256_setzero_ps();
382 __m256 acc2 = _mm256_setzero_ps();
383 __m256 acc3 = _mm256_setzero_ps();
384 for (int32_t
k = -half;
k <= half; ++
k) {
385 const size_t off =
static_cast<size_t>(py +
k) * w
386 + static_cast<size_t>(px);
387 const __m256 kv = _mm256_set1_ps(kdata[
static_cast<size_t>(
k + half)]);
388 acc0 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 0] + off), kv, acc0);
389 acc1 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 1] + off), kv, acc1);
390 acc2 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 2] + off), kv, acc2);
391 acc3 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 3] + off), kv, acc3);
392 }
393 _mm256_storeu_ps(d_ptr[i + 0] + px, acc0);
394 _mm256_storeu_ps(d_ptr[i + 1] + px, acc1);
395 _mm256_storeu_ps(d_ptr[i + 2] + px, acc2);
396 _mm256_storeu_ps(d_ptr[i + 3] + px, acc3);
397 }
398 for (; i < n; ++i) {
399 __m256 acc = _mm256_setzero_ps();
400 for (int32_t
k = -half;
k <= half; ++
k) {
401 const size_t off =
static_cast<size_t>(py +
k) * w
402 + static_cast<size_t>(px);
403 acc = _mm256_fmadd_ps(
404 _mm256_loadu_ps(s_ptr[i] + off),
405 _mm256_set1_ps(kdata[
static_cast<size_t>(
k + half)]),
406 acc);
407 }
408 _mm256_storeu_ps(d_ptr[i] + px, acc);
409 }
410 }
411 }
412 for (; px < iw; ++px) {
413 for (size_t i = 0; i < n; ++i) {
414 float acc = 0.0F;
415 for (int32_t
k = -half;
k <= half; ++
k) {
416 const int32_t ny = border ? std::clamp(py +
k, 0, ih - 1) : py +
k;
417 const size_t off = static_cast<size_t>(ny) * w + static_cast<size_t>(px);
418 acc += s_ptr[i][off] * kdata[
static_cast<size_t>(
k + half)];
419 }
420 d_ptr[i][px] = acc;
421 }
422 }
423
424#elif defined(MAYAFLUX_ARCH_ARM64)
425 int32_t px = 0;
426 if (!border) {
427 for (; px + 4 <= iw; px += 4) {
428 size_t i = 0;
429 for (; i + Unroll - 1 < n; i += Unroll) {
430 float32x4_t acc0 = vdupq_n_f32(0.0F);
431 float32x4_t acc1 = vdupq_n_f32(0.0F);
432 float32x4_t acc2 = vdupq_n_f32(0.0F);
433 float32x4_t acc3 = vdupq_n_f32(0.0F);
434 for (int32_t
k = -half;
k <= half; ++
k) {
435 const size_t off =
static_cast<size_t>(py +
k) * w
436 + static_cast<size_t>(px);
437 const float32x4_t kv = vdupq_n_f32(kdata[
static_cast<size_t>(
k + half)]);
438 acc0 = vmlaq_f32(acc0, vld1q_f32(s_ptr[i + 0] + off), kv);
439 acc1 = vmlaq_f32(acc1, vld1q_f32(s_ptr[i + 1] + off), kv);
440 acc2 = vmlaq_f32(acc2, vld1q_f32(s_ptr[i + 2] + off), kv);
441 acc3 = vmlaq_f32(acc3, vld1q_f32(s_ptr[i + 3] + off), kv);
442 }
443 vst1q_f32(d_ptr[i + 0] + px, acc0);
444 vst1q_f32(d_ptr[i + 1] + px, acc1);
445 vst1q_f32(d_ptr[i + 2] + px, acc2);
446 vst1q_f32(d_ptr[i + 3] + px, acc3);
447 }
448 for (; i < n; ++i) {
449 float32x4_t acc = vdupq_n_f32(0.0F);
450 for (int32_t
k = -half;
k <= half; ++
k) {
451 const size_t off =
static_cast<size_t>(py +
k) * w
452 + static_cast<size_t>(px);
453 acc = vmlaq_f32(acc,
454 vld1q_f32(s_ptr[i] + off),
455 vdupq_n_f32(kdata[
static_cast<size_t>(
k + half)]));
456 }
457 vst1q_f32(d_ptr[i] + px, acc);
458 }
459 }
460 }
461 for (; px < iw; ++px) {
462 for (size_t i = 0; i < n; ++i) {
463 float acc = 0.0F;
464 for (int32_t
k = -half;
k <= half; ++
k) {
465 const int32_t ny = border ? std::clamp(py +
k, 0, ih - 1) : py +
k;
466 const size_t off = static_cast<size_t>(ny) * w + static_cast<size_t>(px);
467 acc += s_ptr[i][off] * kdata[
static_cast<size_t>(
k + half)];
468 }
469 d_ptr[i][px] = acc;
470 }
471 }
472
473#else
474 for (int32_t px = 0; px < iw; ++px) {
475 for (size_t i = 0; i < n; ++i) {
476 float acc = 0.0F;
477 for (int32_t
k = -half;
k <= half; ++
k) {
478 const int32_t ny = std::clamp(py +
k, 0, ih - 1);
479 const size_t off = static_cast<size_t>(ny) * w + static_cast<size_t>(px);
480 acc += s_ptr[i][off] * kdata[
static_cast<size_t>(
k + half)];
481 }
482 d_ptr[i][px] = acc;
483 }
484 }
485#endif
486 });
487}