3#ifdef MAYAFLUX_ARCH_X64
6#ifdef MAYAFLUX_ARCH_ARM64
20 void convolve_horizontal(
21 std::span<const float> src, std::span<float> dst,
22 uint32_t w, uint32_t
h,
23 std::span<const float> kernel)
25 const auto half =
static_cast<int32_t
>(kernel.size() / 2);
26 const auto iw =
static_cast<int32_t
>(w);
28 P::for_each(P::par_unseq,
29 std::views::iota(uint32_t { 0 },
h).begin(),
30 std::views::iota(uint32_t { 0 },
h).end(),
32 const float* src_row = src.data() +
static_cast<size_t>(row) * w;
33 float* dst_row = dst.data() +
static_cast<size_t>(row) * w;
36 for (int32_t px = 0; px < std::min(half, iw); ++px) {
38 for (int32_t
k = -half;
k <= half; ++
k) {
39 acc += src_row[std::clamp(px +
k, 0, iw - 1)]
40 * kernel[
static_cast<size_t>(
k + half)];
45 const int32_t interior_end = iw - half;
47#ifdef MAYAFLUX_ARCH_X64
51 for (; px + 8 <= interior_end; px += 8) {
52 __m256 acc = _mm256_setzero_ps();
53 for (int32_t
k = -half;
k <= half; ++
k) {
54 const __m256 sv = _mm256_loadu_ps(src_row + px +
k);
55 const __m256 kv = _mm256_set1_ps(kernel[
static_cast<size_t>(
k + half)]);
56 acc = _mm256_fmadd_ps(sv, kv, acc);
58 _mm256_storeu_ps(dst_row + px, acc);
61 for (; px < iw; ++px) {
63 for (int32_t
k = -half;
k <= half; ++
k) {
64 acc += src_row[std::clamp(px +
k, 0, iw - 1)]
65 * kernel[
static_cast<size_t>(
k + half)];
70#elif defined(MAYAFLUX_ARCH_ARM64)
73 for (; px + 4 <= interior_end; px += 4) {
74 float32x4_t acc = vdupq_n_f32(0.0F);
75 for (int32_t
k = -half;
k <= half; ++
k) {
76 const float32x4_t sv = vld1q_f32(src_row + px +
k);
77 const float32x4_t kv = vdupq_n_f32(kernel[
static_cast<size_t>(
k + half)]);
78 acc = vmlaq_f32(acc, sv, kv);
80 vst1q_f32(dst_row + px, acc);
82 for (; px < iw; ++px) {
84 for (int32_t
k = -half;
k <= half; ++
k)
85 acc += src_row[std::clamp(px +
k, 0, iw - 1)]
86 * kernel[
static_cast<size_t>(
k + half)];
92 for (int32_t px = half; px < interior_end; ++px) {
94 for (int32_t
k = -half;
k <= half; ++
k)
95 acc += src_row[px +
k] * kernel[
static_cast<size_t>(
k + half)];
99 for (int32_t px = interior_end; px < iw; ++px) {
101 for (int32_t
k = -half;
k <= half; ++
k)
102 acc += src_row[std::clamp(px +
k, 0, iw - 1)]
103 * kernel[
static_cast<size_t>(
k + half)];
110 void convolve_vertical(
111 std::span<const float> src, std::span<float> dst,
112 uint32_t w, uint32_t
h,
113 std::span<const float> kernel)
115 const auto half =
static_cast<int32_t
>(kernel.size() / 2);
116 const auto ih =
static_cast<int32_t
>(
h);
117 const auto iw =
static_cast<int32_t
>(w);
119 P::for_each(P::par_unseq,
120 std::views::iota(uint32_t { 0 },
h).begin(),
121 std::views::iota(uint32_t { 0 },
h).end(),
123 const auto py =
static_cast<int32_t
>(row);
124 float* dst_row = dst.data() +
static_cast<size_t>(row) * w;
125 const bool border = (py < half || py >= ih - half);
127#ifdef MAYAFLUX_ARCH_X64
131 for (; px + 8 <= iw; px += 8) {
132 __m256 acc = _mm256_setzero_ps();
133 for (int32_t
k = -half;
k <= half; ++
k) {
134 const float* src_row = src.data()
135 +
static_cast<size_t>(py +
k) * w + px;
136 const __m256 sv = _mm256_loadu_ps(src_row);
137 const __m256 kv = _mm256_set1_ps(
138 kernel[
static_cast<size_t>(
k + half)]);
139 acc = _mm256_fmadd_ps(sv, kv, acc);
141 _mm256_storeu_ps(dst_row + px, acc);
145 for (; px < iw; ++px) {
147 for (int32_t
k = -half;
k <= half; ++
k) {
148 const int32_t ny = border
149 ? std::clamp(py +
k, 0, ih - 1)
151 acc += src[
static_cast<size_t>(ny) * w + px]
152 * kernel[
static_cast<size_t>(
k + half)];
157#elif defined(MAYAFLUX_ARCH_ARM64)
160 for (; px + 4 <= iw; px += 4) {
161 float32x4_t acc = vdupq_n_f32(0.0F);
162 for (int32_t
k = -half;
k <= half; ++
k) {
163 const float* src_row = src.data()
164 +
static_cast<size_t>(py +
k) * w + px;
165 const float32x4_t sv = vld1q_f32(src_row);
166 const float32x4_t kv = vdupq_n_f32(
167 kernel[
static_cast<size_t>(
k + half)]);
168 acc = vmlaq_f32(acc, sv, kv);
170 vst1q_f32(dst_row + px, acc);
173 for (; px < iw; ++px) {
175 for (int32_t
k = -half;
k <= half; ++
k) {
176 const int32_t ny = border
177 ? std::clamp(py +
k, 0, ih - 1)
179 acc += src[
static_cast<size_t>(ny) * w + px]
180 * kernel[
static_cast<size_t>(
k + half)];
186 for (int32_t px = 0; px < iw; ++px) {
188 for (int32_t
k = -half;
k <= half; ++
k) {
189 const int32_t ny = std::clamp(py +
k, 0, ih - 1);
190 acc += src[
static_cast<size_t>(ny) * w + px]
191 * kernel[
static_cast<size_t>(
k + half)];
202 std::span<const float* const> src,
203 std::span<float* const> dst,
204 uint32_t w, uint32_t
h,
205 std::span<const float> kernel)
207 constexpr size_t Unroll = 4;
209 const auto half =
static_cast<int32_t
>(kernel.size() / 2);
210 const auto iw =
static_cast<int32_t
>(w);
211 const float* kdata = kernel.data();
212 const size_t n = src.size();
214 P::for_each(P::par_unseq,
215 std::views::iota(uint32_t { 0 },
h).begin(),
216 std::views::iota(uint32_t { 0 },
h).end(),
218 const size_t row_off =
static_cast<size_t>(row) * w;
220 const float* s_ptr[32];
222 for (
size_t i = 0; i < n; ++i) {
223 s_ptr[i] = src[i] + row_off;
224 d_ptr[i] = dst[i] + row_off;
228 for (
size_t i = 0; i < n; ++i) {
229 for (int32_t px = 0; px < std::min(half, iw); ++px) {
231 for (int32_t
k = -half;
k <= half; ++
k) {
232 acc += s_ptr[i][std::clamp(px +
k, 0, iw - 1)]
233 * kdata[
static_cast<size_t>(
k + half)];
239 const int32_t interior_end = iw - half;
241#ifdef MAYAFLUX_ARCH_X64
243 for (; px + 8 <= interior_end; px += 8) {
245 for (; i + Unroll - 1 < n; i += Unroll) {
246 __m256 acc0 = _mm256_setzero_ps();
247 __m256 acc1 = _mm256_setzero_ps();
248 __m256 acc2 = _mm256_setzero_ps();
249 __m256 acc3 = _mm256_setzero_ps();
250 for (int32_t
k = -half;
k <= half; ++
k) {
251 const __m256 kv = _mm256_set1_ps(kdata[
static_cast<size_t>(
k + half)]);
252 acc0 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 0] + px +
k), kv, acc0);
253 acc1 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 1] + px +
k), kv, acc1);
254 acc2 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 2] + px +
k), kv, acc2);
255 acc3 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 3] + px +
k), kv, acc3);
257 _mm256_storeu_ps(d_ptr[i + 0] + px, acc0);
258 _mm256_storeu_ps(d_ptr[i + 1] + px, acc1);
259 _mm256_storeu_ps(d_ptr[i + 2] + px, acc2);
260 _mm256_storeu_ps(d_ptr[i + 3] + px, acc3);
263 __m256 acc = _mm256_setzero_ps();
264 for (int32_t
k = -half;
k <= half; ++
k) {
265 acc = _mm256_fmadd_ps(
266 _mm256_loadu_ps(s_ptr[i] + px +
k),
267 _mm256_set1_ps(kdata[
static_cast<size_t>(
k + half)]),
270 _mm256_storeu_ps(d_ptr[i] + px, acc);
274 for (
size_t i = 0; i < n; ++i) {
275 for (int32_t px2 = px; px2 < iw; ++px2) {
277 for (int32_t
k = -half;
k <= half; ++
k) {
278 acc += s_ptr[i][std::clamp(px2 +
k, 0, iw - 1)]
279 * kdata[
static_cast<size_t>(
k + half)];
285#elif defined(MAYAFLUX_ARCH_ARM64)
287 for (; px + 4 <= interior_end; px += 4) {
289 for (; i + Unroll - 1 < n; i += Unroll) {
290 float32x4_t acc0 = vdupq_n_f32(0.0F);
291 float32x4_t acc1 = vdupq_n_f32(0.0F);
292 float32x4_t acc2 = vdupq_n_f32(0.0F);
293 float32x4_t acc3 = vdupq_n_f32(0.0F);
294 for (int32_t
k = -half;
k <= half; ++
k) {
295 const float32x4_t kv = vdupq_n_f32(kdata[
static_cast<size_t>(
k + half)]);
296 acc0 = vmlaq_f32(acc0, vld1q_f32(s_ptr[i + 0] + px +
k), kv);
297 acc1 = vmlaq_f32(acc1, vld1q_f32(s_ptr[i + 1] + px +
k), kv);
298 acc2 = vmlaq_f32(acc2, vld1q_f32(s_ptr[i + 2] + px +
k), kv);
299 acc3 = vmlaq_f32(acc3, vld1q_f32(s_ptr[i + 3] + px +
k), kv);
301 vst1q_f32(d_ptr[i + 0] + px, acc0);
302 vst1q_f32(d_ptr[i + 1] + px, acc1);
303 vst1q_f32(d_ptr[i + 2] + px, acc2);
304 vst1q_f32(d_ptr[i + 3] + px, acc3);
307 float32x4_t acc = vdupq_n_f32(0.0F);
308 for (int32_t
k = -half;
k <= half; ++
k)
310 vld1q_f32(s_ptr[i] + px +
k),
311 vdupq_n_f32(kdata[
static_cast<size_t>(
k + half)]));
312 vst1q_f32(d_ptr[i] + px, acc);
315 for (
size_t i = 0; i < n; ++i) {
316 for (int32_t px2 = px; px2 < iw; ++px2) {
318 for (int32_t
k = -half;
k <= half; ++
k)
319 acc += s_ptr[i][std::clamp(px2 +
k, 0, iw - 1)]
320 * kdata[
static_cast<size_t>(
k + half)];
326 for (
size_t i = 0; i < n; ++i) {
327 for (int32_t px = half; px < interior_end; ++px) {
329 for (int32_t
k = -half;
k <= half; ++
k)
330 acc += s_ptr[i][px +
k] * kdata[
static_cast<size_t>(
k + half)];
333 for (int32_t px = interior_end; px < iw; ++px) {
335 for (int32_t
k = -half;
k <= half; ++
k)
336 acc += s_ptr[i][std::clamp(px +
k, 0, iw - 1)]
337 * kdata[
static_cast<size_t>(
k + half)];
346 std::span<const float* const> src,
347 std::span<float* const> dst,
348 uint32_t w, uint32_t
h,
349 std::span<const float> kernel)
351 constexpr size_t Unroll = 4;
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();
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(),
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;
367 const float* s_ptr[32];
369 for (
size_t i = 0; i < n; ++i) {
371 d_ptr[i] = dst[i] + row_off;
374#ifdef MAYAFLUX_ARCH_X64
377 for (; px + 8 <= iw; px += 8) {
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);
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);
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)]),
408 _mm256_storeu_ps(d_ptr[i] + px, acc);
412 for (; px < iw; ++px) {
413 for (
size_t i = 0; i < n; ++i) {
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)];
424#elif defined(MAYAFLUX_ARCH_ARM64)
427 for (; px + 4 <= iw; px += 4) {
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);
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);
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);
454 vld1q_f32(s_ptr[i] + off),
455 vdupq_n_f32(kdata[
static_cast<size_t>(
k + half)]));
457 vst1q_f32(d_ptr[i] + px, acc);
461 for (; px < iw; ++px) {
462 for (
size_t i = 0; i < n; ++i) {
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)];
474 for (int32_t px = 0; px < iw; ++px) {
475 for (
size_t i = 0; i < n; ++i) {
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)];
490 std::span<const float> src, uint32_t w, uint32_t
h,
491 std::span<const float> kernel_x, std::span<const float> kernel_y)
493 const size_t n =
static_cast<size_t>(w) *
h;
494 std::vector<float> tmp(n);
495 std::vector<float> out(n);
501 std::span<const float> src,
502 std::span<float> tmp,
503 std::span<float> dst,
504 uint32_t w, uint32_t
h,
505 std::span<const float> kernel_x,
506 std::span<const float> kernel_y)
508 convolve_horizontal(src, tmp, w,
h, kernel_x);
509 convolve_vertical(tmp, dst, w,
h, kernel_y);
513 std::span<const float> src, uint32_t w, uint32_t
h,
float sigma)
515 const size_t n =
static_cast<size_t>(w) *
h;
516 std::vector<float> tmp(n);
517 std::vector<float> out(n);
523 std::span<const float> src,
524 std::span<float> tmp,
525 std::span<float> dst,
526 uint32_t w, uint32_t
h,
529 const auto radius =
static_cast<int32_t
>(std::ceil(3.0F *
sigma));
530 const int32_t size = 2 *
radius + 1;
532 std::vector<float> kernel(
static_cast<size_t>(size));
535 const float inv_2s2 = 1.0F / (2.0F *
sigma *
sigma);
537 const float v = std::exp(-
static_cast<float>(i * i) * inv_2s2);
538 kernel[
static_cast<size_t>(i +
radius)] = v;
542 auto kmap = Eigen::Map<Eigen::ArrayXf>(kernel.data(),
static_cast<Eigen::Index
>(kernel.size()));
2D spatial filtering on normalised float image spans.
std::vector< float > filter_separable(std::span< const float > src, uint32_t w, uint32_t h, std::span< const float > kernel_x, std::span< const float > kernel_y)
Apply a separable 2D filter via two 1D passes.
void filter_horizontal_planes(std::span< const float *const > src, std::span< float *const > dst, uint32_t w, uint32_t h, std::span< const float > kernel)
Fused horizontal separable filter over N planes.
std::vector< float > gaussian_blur(std::span< const float > src, uint32_t w, uint32_t h, float sigma)
Separable Gaussian blur.
void filter_vertical_planes(std::span< const float *const > src, std::span< float *const > dst, uint32_t w, uint32_t h, std::span< const float > kernel)
Fused vertical separable filter over N planes.