MayaFlux 0.5.0
Digital-First Multimedia Processing Framework
Loading...
Searching...
No Matches
ImageFilter.cpp
Go to the documentation of this file.
1#include "ImageFilter.hpp"
2
3#ifdef MAYAFLUX_ARCH_X64
4#include <immintrin.h>
5#endif
6#ifdef MAYAFLUX_ARCH_ARM64
7#include <arm_neon.h>
8#endif
9
11
12#include <Eigen/Core>
13
14namespace P = MayaFlux::Parallel;
15
17
18namespace {
19
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)
24 {
25 const auto half = static_cast<int32_t>(kernel.size() / 2);
26 const auto iw = static_cast<int32_t>(w);
27
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(),
31 [&](uint32_t row) {
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;
34
35 // left border: pixels [0, half) require clamped reads
36 for (int32_t px = 0; px < std::min(half, iw); ++px) {
37 float acc = 0.0F;
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)];
41 }
42 dst_row[px] = acc;
43 }
44
45 const int32_t interior_end = iw - half;
46
47#ifdef MAYAFLUX_ARCH_X64
48 // AVX2: 8 output pixels per iteration
49 // interior: all kernel taps land in [0, w) without clamping
50 int32_t px = half;
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);
57 }
58 _mm256_storeu_ps(dst_row + px, acc);
59 }
60 // scalar tail (fewer than 8 pixels remain in interior + right border)
61 for (; px < iw; ++px) {
62 float acc = 0.0F;
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)];
66 }
67 dst_row[px] = acc;
68 }
69
70#elif defined(MAYAFLUX_ARCH_ARM64)
71 // NEON: 4 output pixels per iteration
72 int32_t px = half;
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);
79 }
80 vst1q_f32(dst_row + px, acc);
81 }
82 for (; px < iw; ++px) {
83 float acc = 0.0F;
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)];
87 dst_row[px] = acc;
88 }
89
90#else
91 // scalar interior (no clamping needed)
92 for (int32_t px = half; px < interior_end; ++px) {
93 float acc = 0.0F;
94 for (int32_t k = -half; k <= half; ++k)
95 acc += src_row[px + k] * kernel[static_cast<size_t>(k + half)];
96 dst_row[px] = acc;
97 }
98 // right border
99 for (int32_t px = interior_end; px < iw; ++px) {
100 float acc = 0.0F;
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)];
104 dst_row[px] = acc;
105 }
106#endif
107 });
108 }
109
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)
114 {
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);
118
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(),
122 [&](uint32_t row) {
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);
126
127#ifdef MAYAFLUX_ARCH_X64
128 int32_t px = 0;
129 if (!border) {
130 // interior row: no clamping needed, 8 columns at a time
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);
140 }
141 _mm256_storeu_ps(dst_row + px, acc);
142 }
143 }
144 // scalar tail and all border rows
145 for (; px < iw; ++px) {
146 float acc = 0.0F;
147 for (int32_t k = -half; k <= half; ++k) {
148 const int32_t ny = border
149 ? std::clamp(py + k, 0, ih - 1)
150 : py + k;
151 acc += src[static_cast<size_t>(ny) * w + px]
152 * kernel[static_cast<size_t>(k + half)];
153 }
154 dst_row[px] = acc;
155 }
156
157#elif defined(MAYAFLUX_ARCH_ARM64)
158 int32_t px = 0;
159 if (!border) {
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);
169 }
170 vst1q_f32(dst_row + px, acc);
171 }
172 }
173 for (; px < iw; ++px) {
174 float acc = 0.0F;
175 for (int32_t k = -half; k <= half; ++k) {
176 const int32_t ny = border
177 ? std::clamp(py + k, 0, ih - 1)
178 : py + k;
179 acc += src[static_cast<size_t>(ny) * w + px]
180 * kernel[static_cast<size_t>(k + half)];
181 }
182 dst_row[px] = acc;
183 }
184
185#else
186 for (int32_t px = 0; px < iw; ++px) {
187 float acc = 0.0F;
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)];
192 }
193 dst_row[px] = acc;
194 }
195#endif
196 });
197 }
198
199} // namespace
200
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)
206{
207 constexpr size_t Unroll = 4;
208
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();
213
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(),
217 [&](uint32_t row) {
218 const size_t row_off = static_cast<size_t>(row) * w;
219
220 const float* s_ptr[32];
221 float* d_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;
225 }
226
227 // left border
228 for (size_t i = 0; i < n; ++i) {
229 for (int32_t px = 0; px < std::min(half, iw); ++px) {
230 float acc = 0.0F;
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)];
234 }
235 d_ptr[i][px] = acc;
236 }
237 }
238
239 const int32_t interior_end = iw - half;
240
241#ifdef MAYAFLUX_ARCH_X64
242 int32_t px = half;
243 for (; px + 8 <= interior_end; px += 8) {
244 size_t i = 0;
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);
256 }
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);
261 }
262 for (; i < n; ++i) {
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)]),
268 acc);
269 }
270 _mm256_storeu_ps(d_ptr[i] + px, acc);
271 }
272 }
273 // right border + scalar tail
274 for (size_t i = 0; i < n; ++i) {
275 for (int32_t px2 = px; px2 < iw; ++px2) {
276 float acc = 0.0F;
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)];
280 }
281 d_ptr[i][px2] = acc;
282 }
283 }
284
285#elif defined(MAYAFLUX_ARCH_ARM64)
286 int32_t px = half;
287 for (; px + 4 <= interior_end; px += 4) {
288 size_t i = 0;
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);
300 }
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);
305 }
306 for (; i < n; ++i) {
307 float32x4_t acc = vdupq_n_f32(0.0F);
308 for (int32_t k = -half; k <= half; ++k)
309 acc = vmlaq_f32(acc,
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);
313 }
314 }
315 for (size_t i = 0; i < n; ++i) {
316 for (int32_t px2 = px; px2 < iw; ++px2) {
317 float acc = 0.0F;
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)];
321 d_ptr[i][px2] = acc;
322 }
323 }
324
325#else
326 for (size_t i = 0; i < n; ++i) {
327 for (int32_t px = half; px < interior_end; ++px) {
328 float acc = 0.0F;
329 for (int32_t k = -half; k <= half; ++k)
330 acc += s_ptr[i][px + k] * kdata[static_cast<size_t>(k + half)];
331 d_ptr[i][px] = acc;
332 }
333 for (int32_t px = interior_end; px < iw; ++px) {
334 float acc = 0.0F;
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)];
338 d_ptr[i][px] = acc;
339 }
340 }
341#endif
342 });
343}
344
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)
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}
488
489std::vector<float> filter_separable(
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)
492{
493 const size_t n = static_cast<size_t>(w) * h;
494 std::vector<float> tmp(n);
495 std::vector<float> out(n);
496 filter_separable(src, tmp, out, w, h, kernel_x, kernel_y);
497 return out;
498}
499
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)
507{
508 convolve_horizontal(src, tmp, w, h, kernel_x);
509 convolve_vertical(tmp, dst, w, h, kernel_y);
510}
511
512std::vector<float> gaussian_blur(
513 std::span<const float> src, uint32_t w, uint32_t h, float sigma)
514{
515 const size_t n = static_cast<size_t>(w) * h;
516 std::vector<float> tmp(n);
517 std::vector<float> out(n);
518 gaussian_blur(src, tmp, out, w, h, sigma);
519 return out;
520}
521
523 std::span<const float> src,
524 std::span<float> tmp,
525 std::span<float> dst,
526 uint32_t w, uint32_t h,
527 float sigma)
528{
529 const auto radius = static_cast<int32_t>(std::ceil(3.0F * sigma));
530 const int32_t size = 2 * radius + 1;
531
532 std::vector<float> kernel(static_cast<size_t>(size));
533 float sum = 0.0F;
534
535 const float inv_2s2 = 1.0F / (2.0F * sigma * sigma);
536 for (int32_t i = -radius; i <= radius; ++i) {
537 const float v = std::exp(-static_cast<float>(i * i) * inv_2s2);
538 kernel[static_cast<size_t>(i + radius)] = v;
539 sum += v;
540 }
541
542 auto kmap = Eigen::Map<Eigen::ArrayXf>(kernel.data(), static_cast<Eigen::Index>(kernel.size()));
543 kmap /= sum;
544 filter_separable(src, tmp, dst, w, h, kernel, kernel);
545}
546
547} // namespace MayaFlux::Kinesis::Vision
2D spatial filtering on normalised float image spans.
uint32_t h
Definition InkPress.cpp:28
uint32_t radius
float k
float sigma
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.