3#ifdef MAYAFLUX_ARCH_X64
6#ifdef MAYAFLUX_ARCH_ARM64
18 void convolve_horizontal(
19 std::span<const float>
src, std::span<float>
dst,
20 uint32_t w, uint32_t
h,
21 std::span<const float> kernel)
23 const auto half =
static_cast<int32_t
>(kernel.size() / 2);
24 const auto iw =
static_cast<int32_t
>(w);
26 std::for_each(std::execution::par_unseq,
27 std::views::iota(uint32_t { 0 },
h).begin(),
28 std::views::iota(uint32_t { 0 },
h).end(),
30 const float* src_row =
src.data() +
static_cast<size_t>(row) * w;
31 float* dst_row =
dst.data() +
static_cast<size_t>(row) * w;
34 for (int32_t
px = 0;
px < std::min(half, iw); ++
px) {
36 for (int32_t k = -half; k <= half; ++k) {
37 acc += src_row[std::clamp(
px + k, 0, iw - 1)]
38 * kernel[
static_cast<size_t>(k + half)];
43 const int32_t interior_end = iw - half;
45#ifdef MAYAFLUX_ARCH_X64
49 for (;
px + 8 <= interior_end;
px += 8) {
50 __m256 acc = _mm256_setzero_ps();
51 for (int32_t k = -half; k <= half; ++k) {
52 const __m256 sv = _mm256_loadu_ps(src_row +
px + k);
53 const __m256 kv = _mm256_set1_ps(kernel[
static_cast<size_t>(k + half)]);
54 acc = _mm256_fmadd_ps(sv, kv, acc);
56 _mm256_storeu_ps(dst_row +
px, acc);
59 for (;
px < iw; ++
px) {
61 for (int32_t k = -half; k <= half; ++k) {
62 acc += src_row[std::clamp(
px + k, 0, iw - 1)]
63 * kernel[
static_cast<size_t>(k + half)];
68#elif defined(MAYAFLUX_ARCH_ARM64)
71 for (;
px + 4 <= interior_end;
px += 4) {
72 float32x4_t acc = vdupq_n_f32(0.0F);
73 for (int32_t k = -half; k <= half; ++k) {
74 const float32x4_t sv = vld1q_f32(src_row +
px + k);
75 const float32x4_t kv = vdupq_n_f32(kernel[
static_cast<size_t>(k + half)]);
76 acc = vmlaq_f32(acc, sv, kv);
78 vst1q_f32(dst_row +
px, acc);
80 for (;
px < iw; ++
px) {
82 for (int32_t k = -half; k <= half; ++k)
83 acc += src_row[std::clamp(
px + k, 0, iw - 1)]
84 * kernel[
static_cast<size_t>(k + half)];
90 for (int32_t
px = half;
px < interior_end; ++
px) {
92 for (int32_t k = -half; k <= half; ++k)
93 acc += src_row[
px + k] * kernel[
static_cast<size_t>(k + half)];
97 for (int32_t
px = interior_end;
px < iw; ++
px) {
99 for (int32_t k = -half; k <= half; ++k)
100 acc += src_row[std::clamp(
px + k, 0, iw - 1)]
101 * kernel[
static_cast<size_t>(k + half)];
108 void convolve_vertical(
109 std::span<const float>
src, std::span<float>
dst,
110 uint32_t w, uint32_t
h,
111 std::span<const float> kernel)
113 const auto half =
static_cast<int32_t
>(kernel.size() / 2);
114 const auto ih =
static_cast<int32_t
>(
h);
115 const auto iw =
static_cast<int32_t
>(w);
117 std::for_each(std::execution::par_unseq,
118 std::views::iota(uint32_t { 0 },
h).begin(),
119 std::views::iota(uint32_t { 0 },
h).end(),
121 const auto py =
static_cast<int32_t
>(row);
122 float* dst_row =
dst.data() +
static_cast<size_t>(row) * w;
123 const bool border = (py < half || py >= ih - half);
125#ifdef MAYAFLUX_ARCH_X64
129 for (;
px + 8 <= iw;
px += 8) {
130 __m256 acc = _mm256_setzero_ps();
131 for (int32_t k = -half; k <= half; ++k) {
132 const float* src_row =
src.data()
133 +
static_cast<size_t>(
py + k) * w +
px;
134 const __m256 sv = _mm256_loadu_ps(src_row);
135 const __m256 kv = _mm256_set1_ps(
136 kernel[
static_cast<size_t>(k + half)]);
137 acc = _mm256_fmadd_ps(sv, kv, acc);
139 _mm256_storeu_ps(dst_row +
px, acc);
143 for (;
px < iw; ++
px) {
145 for (int32_t k = -half; k <= half; ++k) {
146 const int32_t ny = border
147 ? std::clamp(
py + k, 0, ih - 1)
149 acc +=
src[
static_cast<size_t>(ny) * w +
px]
150 * kernel[
static_cast<size_t>(k + half)];
155#elif defined(MAYAFLUX_ARCH_ARM64)
158 for (;
px + 4 <= iw;
px += 4) {
159 float32x4_t acc = vdupq_n_f32(0.0F);
160 for (int32_t k = -half; k <= half; ++k) {
161 const float* src_row =
src.data()
162 +
static_cast<size_t>(
py + k) * w +
px;
163 const float32x4_t sv = vld1q_f32(src_row);
164 const float32x4_t kv = vdupq_n_f32(
165 kernel[
static_cast<size_t>(k + half)]);
166 acc = vmlaq_f32(acc, sv, kv);
168 vst1q_f32(dst_row +
px, acc);
171 for (;
px < iw; ++
px) {
173 for (int32_t k = -half; k <= half; ++k) {
174 const int32_t ny = border
175 ? std::clamp(
py + k, 0, ih - 1)
177 acc +=
src[
static_cast<size_t>(ny) * w +
px]
178 * kernel[
static_cast<size_t>(k + half)];
184 for (int32_t
px = 0;
px < iw; ++
px) {
186 for (int32_t k = -half; k <= half; ++k) {
187 const int32_t ny = std::clamp(
py + k, 0, ih - 1);
188 acc +=
src[
static_cast<size_t>(ny) * w +
px]
189 * kernel[
static_cast<size_t>(k + half)];
200 std::span<const float* const>
src,
201 std::span<float* const>
dst,
202 uint32_t w, uint32_t
h,
203 std::span<const float> kernel)
205 constexpr size_t Unroll = 4;
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();
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(),
216 const size_t row_off =
static_cast<size_t>(row) * w;
218 const float* s_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;
226 for (
size_t i = 0; i < n; ++i) {
227 for (int32_t
px = 0;
px < std::min(half, iw); ++
px) {
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)];
237 const int32_t interior_end = iw - half;
239#ifdef MAYAFLUX_ARCH_X64
241 for (;
px + 8 <= interior_end;
px += 8) {
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);
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);
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)]),
268 _mm256_storeu_ps(d_ptr[i] +
px, acc);
272 for (
size_t i = 0; i < n; ++i) {
273 for (int32_t px2 =
px; px2 < iw; ++px2) {
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)];
283#elif defined(MAYAFLUX_ARCH_ARM64)
285 for (;
px + 4 <= interior_end;
px += 4) {
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);
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);
305 float32x4_t acc = vdupq_n_f32(0.0F);
306 for (int32_t k = -half; k <= half; ++k)
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);
313 for (
size_t i = 0; i < n; ++i) {
314 for (int32_t px2 =
px; px2 < iw; ++px2) {
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)];
324 for (
size_t i = 0; i < n; ++i) {
325 for (int32_t
px = half;
px < interior_end; ++
px) {
327 for (int32_t k = -half; k <= half; ++k)
328 acc += s_ptr[i][
px + k] * kdata[
static_cast<size_t>(k + half)];
331 for (int32_t
px = interior_end;
px < iw; ++
px) {
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)];
344 std::span<const float* const>
src,
345 std::span<float* const>
dst,
346 uint32_t w, uint32_t
h,
347 std::span<const float> kernel)
349 constexpr size_t Unroll = 4;
351 const auto half =
static_cast<int32_t
>(kernel.size() / 2);
352 const auto ih =
static_cast<int32_t
>(
h);
353 const auto iw =
static_cast<int32_t
>(w);
354 const float* kdata = kernel.data();
355 const size_t n =
src.size();
357 std::for_each(std::execution::par_unseq,
358 std::views::iota(uint32_t { 0 },
h).begin(),
359 std::views::iota(uint32_t { 0 },
h).end(),
361 const auto py =
static_cast<int32_t
>(row);
362 const bool border = (py < half || py >= ih - half);
363 const size_t row_off =
static_cast<size_t>(row) * w;
365 const float* s_ptr[32];
367 for (
size_t i = 0; i < n; ++i) {
369 d_ptr[i] =
dst[i] + row_off;
372#ifdef MAYAFLUX_ARCH_X64
375 for (;
px + 8 <= iw;
px += 8) {
377 for (; i + Unroll - 1 < n; i += Unroll) {
378 __m256 acc0 = _mm256_setzero_ps();
379 __m256 acc1 = _mm256_setzero_ps();
380 __m256 acc2 = _mm256_setzero_ps();
381 __m256 acc3 = _mm256_setzero_ps();
382 for (int32_t k = -half; k <= half; ++k) {
383 const size_t off =
static_cast<size_t>(
py + k) * w
384 +
static_cast<size_t>(
px);
385 const __m256 kv = _mm256_set1_ps(kdata[
static_cast<size_t>(k + half)]);
386 acc0 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 0] + off), kv, acc0);
387 acc1 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 1] + off), kv, acc1);
388 acc2 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 2] + off), kv, acc2);
389 acc3 = _mm256_fmadd_ps(_mm256_loadu_ps(s_ptr[i + 3] + off), kv, acc3);
391 _mm256_storeu_ps(d_ptr[i + 0] +
px, acc0);
392 _mm256_storeu_ps(d_ptr[i + 1] +
px, acc1);
393 _mm256_storeu_ps(d_ptr[i + 2] +
px, acc2);
394 _mm256_storeu_ps(d_ptr[i + 3] +
px, acc3);
397 __m256 acc = _mm256_setzero_ps();
398 for (int32_t k = -half; k <= half; ++k) {
399 const size_t off =
static_cast<size_t>(
py + k) * w
400 +
static_cast<size_t>(
px);
401 acc = _mm256_fmadd_ps(
402 _mm256_loadu_ps(s_ptr[i] + off),
403 _mm256_set1_ps(kdata[
static_cast<size_t>(k + half)]),
406 _mm256_storeu_ps(d_ptr[i] +
px, acc);
410 for (;
px < iw; ++
px) {
411 for (
size_t i = 0; i < n; ++i) {
413 for (int32_t k = -half; k <= half; ++k) {
414 const int32_t ny = border ? std::clamp(
py + k, 0, ih - 1) :
py + k;
415 const size_t off =
static_cast<size_t>(ny) * w +
static_cast<size_t>(
px);
416 acc += s_ptr[i][off] * kdata[
static_cast<size_t>(k + half)];
422#elif defined(MAYAFLUX_ARCH_ARM64)
425 for (;
px + 4 <= iw;
px += 4) {
427 for (; i + Unroll - 1 < n; i += Unroll) {
428 float32x4_t acc0 = vdupq_n_f32(0.0F);
429 float32x4_t acc1 = vdupq_n_f32(0.0F);
430 float32x4_t acc2 = vdupq_n_f32(0.0F);
431 float32x4_t acc3 = vdupq_n_f32(0.0F);
432 for (int32_t k = -half; k <= half; ++k) {
433 const size_t off =
static_cast<size_t>(
py + k) * w
434 +
static_cast<size_t>(
px);
435 const float32x4_t kv = vdupq_n_f32(kdata[
static_cast<size_t>(k + half)]);
436 acc0 = vmlaq_f32(acc0, vld1q_f32(s_ptr[i + 0] + off), kv);
437 acc1 = vmlaq_f32(acc1, vld1q_f32(s_ptr[i + 1] + off), kv);
438 acc2 = vmlaq_f32(acc2, vld1q_f32(s_ptr[i + 2] + off), kv);
439 acc3 = vmlaq_f32(acc3, vld1q_f32(s_ptr[i + 3] + off), kv);
441 vst1q_f32(d_ptr[i + 0] +
px, acc0);
442 vst1q_f32(d_ptr[i + 1] +
px, acc1);
443 vst1q_f32(d_ptr[i + 2] +
px, acc2);
444 vst1q_f32(d_ptr[i + 3] +
px, acc3);
447 float32x4_t acc = vdupq_n_f32(0.0F);
448 for (int32_t k = -half; k <= half; ++k) {
449 const size_t off =
static_cast<size_t>(
py + k) * w
450 +
static_cast<size_t>(
px);
452 vld1q_f32(s_ptr[i] + off),
453 vdupq_n_f32(kdata[
static_cast<size_t>(k + half)]));
455 vst1q_f32(d_ptr[i] +
px, acc);
459 for (;
px < iw; ++
px) {
460 for (
size_t i = 0; i < n; ++i) {
462 for (int32_t k = -half; k <= half; ++k) {
463 const int32_t ny = border ? std::clamp(
py + k, 0, ih - 1) :
py + k;
464 const size_t off =
static_cast<size_t>(ny) * w +
static_cast<size_t>(
px);
465 acc += s_ptr[i][off] * kdata[
static_cast<size_t>(k + half)];
472 for (int32_t
px = 0;
px < iw; ++
px) {
473 for (
size_t i = 0; i < n; ++i) {
475 for (int32_t k = -half; k <= half; ++k) {
476 const int32_t ny = std::clamp(
py + k, 0, ih - 1);
477 const size_t off =
static_cast<size_t>(ny) * w +
static_cast<size_t>(
px);
478 acc += s_ptr[i][off] * kdata[
static_cast<size_t>(k + half)];
488 std::span<const float>
src, uint32_t w, uint32_t
h,
489 std::span<const float> kernel_x, std::span<const float> kernel_y)
491 const size_t n =
static_cast<size_t>(w) *
h;
492 std::vector<float> tmp(n);
493 std::vector<float>
out(n);
499 std::span<const float>
src,
500 std::span<float> tmp,
501 std::span<float>
dst,
502 uint32_t w, uint32_t
h,
503 std::span<const float> kernel_x,
504 std::span<const float> kernel_y)
506 convolve_horizontal(
src, tmp, w,
h, kernel_x);
507 convolve_vertical(tmp,
dst, w,
h, kernel_y);
511 std::span<const float>
src, uint32_t w, uint32_t
h,
float sigma)
513 const size_t n =
static_cast<size_t>(w) *
h;
514 std::vector<float> tmp(n);
515 std::vector<float>
out(n);
521 std::span<const float>
src,
522 std::span<float> tmp,
523 std::span<float>
dst,
524 uint32_t w, uint32_t
h,
527 const auto radius =
static_cast<int32_t
>(std::ceil(3.0F * sigma));
530 std::vector<float> kernel(
static_cast<size_t>(
size));
533 const float inv_2s2 = 1.0F / (2.0F * sigma * sigma);
535 const float v = std::exp(-
static_cast<float>(i * i) * inv_2s2);
536 kernel[
static_cast<size_t>(i +
radius)] = v;
540 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 > * out
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.