MayaFlux 1.0.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
10#include <execution>
11
12#include <Eigen/Core>
13
15
16namespace {
17
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)
22 {
23 const auto half = static_cast<int32_t>(kernel.size() / 2);
24 const auto iw = static_cast<int32_t>(w);
25
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(),
29 [&](uint32_t row) {
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;
32
33 // left border: pixels [0, half) require clamped reads
34 for (int32_t px = 0; px < std::min(half, iw); ++px) {
35 float acc = 0.0F;
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)];
39 }
40 dst_row[px] = acc;
41 }
42
43 const int32_t interior_end = iw - half;
44
45#ifdef MAYAFLUX_ARCH_X64
46 // AVX2: 8 output pixels per iteration
47 // interior: all kernel taps land in [0, w) without clamping
48 int32_t px = half;
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);
55 }
56 _mm256_storeu_ps(dst_row + px, acc);
57 }
58 // scalar tail (fewer than 8 pixels remain in interior + right border)
59 for (; px < iw; ++px) {
60 float acc = 0.0F;
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)];
64 }
65 dst_row[px] = acc;
66 }
67
68#elif defined(MAYAFLUX_ARCH_ARM64)
69 // NEON: 4 output pixels per iteration
70 int32_t px = half;
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);
77 }
78 vst1q_f32(dst_row + px, acc);
79 }
80 for (; px < iw; ++px) {
81 float acc = 0.0F;
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)];
85 dst_row[px] = acc;
86 }
87
88#else
89 // scalar interior (no clamping needed)
90 for (int32_t px = half; px < interior_end; ++px) {
91 float acc = 0.0F;
92 for (int32_t k = -half; k <= half; ++k)
93 acc += src_row[px + k] * kernel[static_cast<size_t>(k + half)];
94 dst_row[px] = acc;
95 }
96 // right border
97 for (int32_t px = interior_end; px < iw; ++px) {
98 float acc = 0.0F;
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)];
102 dst_row[px] = acc;
103 }
104#endif
105 });
106 }
107
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)
112 {
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);
116
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(),
120 [&](uint32_t row) {
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);
124
125#ifdef MAYAFLUX_ARCH_X64
126 int32_t px = 0;
127 if (!border) {
128 // interior row: no clamping needed, 8 columns at a time
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);
138 }
139 _mm256_storeu_ps(dst_row + px, acc);
140 }
141 }
142 // scalar tail and all border rows
143 for (; px < iw; ++px) {
144 float acc = 0.0F;
145 for (int32_t k = -half; k <= half; ++k) {
146 const int32_t ny = border
147 ? std::clamp(py + k, 0, ih - 1)
148 : py + k;
149 acc += src[static_cast<size_t>(ny) * w + px]
150 * kernel[static_cast<size_t>(k + half)];
151 }
152 dst_row[px] = acc;
153 }
154
155#elif defined(MAYAFLUX_ARCH_ARM64)
156 int32_t px = 0;
157 if (!border) {
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);
167 }
168 vst1q_f32(dst_row + px, acc);
169 }
170 }
171 for (; px < iw; ++px) {
172 float acc = 0.0F;
173 for (int32_t k = -half; k <= half; ++k) {
174 const int32_t ny = border
175 ? std::clamp(py + k, 0, ih - 1)
176 : py + k;
177 acc += src[static_cast<size_t>(ny) * w + px]
178 * kernel[static_cast<size_t>(k + half)];
179 }
180 dst_row[px] = acc;
181 }
182
183#else
184 for (int32_t px = 0; px < iw; ++px) {
185 float acc = 0.0F;
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)];
190 }
191 dst_row[px] = acc;
192 }
193#endif
194 });
195 }
196
197} // namespace
198
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)
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 // left border
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 }
233 d_ptr[i][px] = acc;
234 }
235 }
236
237 const int32_t interior_end = iw - half;
238
239#ifdef MAYAFLUX_ARCH_X64
240 int32_t px = half;
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 // right border + scalar tail
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)
284 int32_t px = half;
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)];
329 d_ptr[i][px] = acc;
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)];
336 d_ptr[i][px] = acc;
337 }
338 }
339#endif
340 });
341}
342
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)
348{
349 constexpr size_t Unroll = 4;
350
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();
356
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(),
360 [&](uint32_t row) {
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;
364
365 const float* s_ptr[32];
366 float* d_ptr[32];
367 for (size_t i = 0; i < n; ++i) {
368 s_ptr[i] = src[i];
369 d_ptr[i] = dst[i] + row_off;
370 }
371
372#ifdef MAYAFLUX_ARCH_X64
373 int32_t px = 0;
374 if (!border) {
375 for (; px + 8 <= iw; px += 8) {
376 size_t i = 0;
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);
390 }
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);
395 }
396 for (; i < n; ++i) {
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)]),
404 acc);
405 }
406 _mm256_storeu_ps(d_ptr[i] + px, acc);
407 }
408 }
409 }
410 for (; px < iw; ++px) {
411 for (size_t i = 0; i < n; ++i) {
412 float acc = 0.0F;
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)];
417 }
418 d_ptr[i][px] = acc;
419 }
420 }
421
422#elif defined(MAYAFLUX_ARCH_ARM64)
423 int32_t px = 0;
424 if (!border) {
425 for (; px + 4 <= iw; px += 4) {
426 size_t i = 0;
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);
440 }
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);
445 }
446 for (; i < n; ++i) {
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);
451 acc = vmlaq_f32(acc,
452 vld1q_f32(s_ptr[i] + off),
453 vdupq_n_f32(kdata[static_cast<size_t>(k + half)]));
454 }
455 vst1q_f32(d_ptr[i] + px, acc);
456 }
457 }
458 }
459 for (; px < iw; ++px) {
460 for (size_t i = 0; i < n; ++i) {
461 float acc = 0.0F;
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)];
466 }
467 d_ptr[i][px] = acc;
468 }
469 }
470
471#else
472 for (int32_t px = 0; px < iw; ++px) {
473 for (size_t i = 0; i < n; ++i) {
474 float acc = 0.0F;
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)];
479 }
480 d_ptr[i][px] = acc;
481 }
482 }
483#endif
484 });
485}
486
487std::vector<float> filter_separable(
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)
490{
491 const size_t n = static_cast<size_t>(w) * h;
492 std::vector<float> tmp(n);
493 std::vector<float> out(n);
494 filter_separable(src, tmp, out, w, h, kernel_x, kernel_y);
495 return out;
496}
497
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)
505{
506 convolve_horizontal(src, tmp, w, h, kernel_x);
507 convolve_vertical(tmp, dst, w, h, kernel_y);
508}
509
510std::vector<float> gaussian_blur(
511 std::span<const float> src, uint32_t w, uint32_t h, float sigma)
512{
513 const size_t n = static_cast<size_t>(w) * h;
514 std::vector<float> tmp(n);
515 std::vector<float> out(n);
516 gaussian_blur(src, tmp, out, w, h, sigma);
517 return out;
518}
519
521 std::span<const float> src,
522 std::span<float> tmp,
523 std::span<float> dst,
524 uint32_t w, uint32_t h,
525 float sigma)
526{
527 const auto radius = static_cast<int32_t>(std::ceil(3.0F * sigma));
528 const int32_t size = 2 * radius + 1;
529
530 std::vector<float> kernel(static_cast<size_t>(size));
531 float sum = 0.0F;
532
533 const float inv_2s2 = 1.0F / (2.0F * sigma * sigma);
534 for (int32_t i = -radius; i <= radius; ++i) {
535 const float v = std::exp(-static_cast<float>(i * i) * inv_2s2);
536 kernel[static_cast<size_t>(i + radius)] = v;
537 sum += v;
538 }
539
540 auto kmap = Eigen::Map<Eigen::ArrayXf>(kernel.data(), static_cast<Eigen::Index>(kernel.size()));
541 kmap /= sum;
542 filter_separable(src, tmp, dst, w, h, kernel, kernel);
543}
544
545} // namespace MayaFlux::Kinesis::Vision
bool size
2D spatial filtering on normalised float image spans.
float radius
uint32_t h
Definition InkPress.cpp:29
uint32_t * dst
const uint32_t * src
std::vector< float > * out
uint32_t py
uint32_t px
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.