Merge pull request #80 from awxkee/f16

Better complex multiplication

Author: awxkee

.github/workflows/build_push.yml

@@ -134,6 +134,7 @@ jobs:
       - uses: dtolnay/rust-toolchain@nightly
       - run: cargo install cargo-fuzz
       - run: cargo fuzz run filter_complex --features neon,nightly_fcma -- -max_total_time=17
+      - run: cargo fuzz run motion --features neon,nightly_fcma -- -max_total_time=17
 
   fuzz_filters_x86:
     name: Fuzzing Filters 1D/2D x86

Cargo.toml

@@ -2,7 +2,7 @@ workspace = { members = ["fuzz", "app/accelerate", "app"] }
 
 [package]
 name = "libblur"
-version = "0.20.0"
+version = "0.20.1"
 edition = "2021"
 description = "Fast image blurring in pure Rust"
 readme = "./README.md"

app/benches/gauss/main.rs

@@ -60,7 +60,7 @@ pub fn criterion_benchmark(c: &mut Criterion) {
                 &src_image,
                 &mut dst_bytes,
                 GaussianBlurParams::new_from_kernel(3.),
-                EdgeMode2D::new(EdgeMode::Clamp.as_2d()),
+                EdgeMode::Clamp.as_2d(),
                 ThreadingPolicy::Adaptive,
                 ConvolutionMode::FixedPoint,
             )
@@ -107,7 +107,7 @@ pub fn criterion_benchmark(c: &mut Criterion) {
                 &src_image,
                 &mut dst_bytes,
                 GaussianBlurParams::new_from_kernel(13.),
-                EdgeMode2D::new(EdgeMode::Clamp.as_2d()),
+                EdgeMode::Clamp.as_2d(),
                 ThreadingPolicy::Adaptive,
                 ConvolutionMode::Exact,
             )
@@ -123,7 +123,7 @@ pub fn criterion_benchmark(c: &mut Criterion) {
                 &src_image,
                 &mut dst_bytes,
                 GaussianBlurParams::new_from_kernel(13.),
-                EdgeMode2D::new(EdgeMode::Clamp.as_2d()),
+                EdgeMode::Clamp.as_2d(),
                 ThreadingPolicy::Adaptive,
                 ConvolutionMode::FixedPoint,
             )

app/src/main.rs

@@ -36,7 +36,7 @@ use libblur::{
     bilateral_filter, complex_gaussian_kernel, fast_bilateral_filter, fast_bilateral_filter_u16,
     filter_1d_complex, filter_1d_complex_fixed_point, gaussian_blur, gaussian_kernel_1d,
     lens_kernel, sigma_size, AnisotropicRadius, BilateralBlurParams, BlurImage, BlurImageMut,
-    BoxBlurParameters, CLTParameters, ConvolutionMode, EdgeMode, FastBlurChannels,
+    BoxBlurParameters, CLTParameters, ConvolutionMode, EdgeMode, EdgeMode2D, FastBlurChannels,
     GaussianBlurParams, KernelShape, Scalar, ThreadingPolicy, TransferFunction,
 };
 use num_complex::Complex;
@@ -80,9 +80,9 @@ fn main() {
 
     println!("{:?}", dyn_image.color());
 
-    let img = dyn_image.to_rgba8();
+    let img = dyn_image.to_rgb8();
     let mut src_bytes = img.as_bytes();
-    let components = 4;
+    let components = 3;
     let stride = dimensions.0 as usize * components;
     let mut bytes: Vec<u8> = src_bytes.to_vec();
     let mut dst_bytes: Vec<u8> = src_bytes.to_vec();
@@ -92,10 +92,10 @@ fn main() {
     let mut v_vec = src_bytes
         .to_vec()
         .iter()
-        // .map(|&x| x)
+        .map(|&x| x)
         // .map(|&x| (x as f32 / 255.))
-        .map(|&x| u16::from_ne_bytes([x, x]))
-        .collect::<Vec<u16>>();
+        // .map(|&x| u16::from_ne_bytes([x, x]))
+        .collect::<Vec<u8>>();
 
     // let mut dst_image = BlurImageMut::borrow(
     //     &mut v_vec,
@@ -109,7 +109,7 @@ fn main() {
         &v_vec,
         dyn_image.width(),
         dyn_image.height(),
-        FastBlurChannels::Channels4,
+        FastBlurChannels::Channels3,
     );
     // let vcvt = cvt.linearize(TransferFunction::Srgb, true).unwrap();
 
@@ -157,13 +157,11 @@ fn main() {
 
     // }
 
-    libblur::box_blur_u16(
+    libblur::sobel(
         &cvt,
         &mut dst_image,
-        BoxBlurParameters {
-            x_axis_kernel: 7,
-            y_axis_kernel: 7,
-        },
+        EdgeMode2D::default(),
+        Scalar::default(),
         ThreadingPolicy::Single,
     )
     .unwrap();
@@ -186,17 +184,17 @@ fn main() {
     // )
     // .unwrap();
 
-    let j_dag = dst_image.to_immutable_ref();
+    // let j_dag = dst_image.to_immutable_ref();
 
     // let gamma = j_dag.gamma8(TransferFunction::Srgb, true).unwrap();
 
     dst_bytes = dst_image
         .data
         .borrow_mut()
         .iter()
-        // .map(|&x| x)
+        .map(|&x| x)
         // .map(|&x| (x * 255f32).round() as u8)
-        .map(|&x| (x >> 8) as u8)
+        // .map(|&x| (x >> 8) as u8)
         .collect::<Vec<u8>>();
 
     // dst_bytes = dst_image.data.borrow().to_vec();
@@ -225,7 +223,7 @@ fn main() {
 
     if components == 3 {
         image::save_buffer(
-            "blurred_stack_next.jpg",
+            "blurred_stack_next1.jpg",
             bytes.as_bytes(),
             dimensions.0,
             dimensions.1,
@@ -234,7 +232,7 @@ fn main() {
         .unwrap();
     } else {
         image::save_buffer(
-            "blurred_stack_next_f.png",
+            "blurred_stack_next_f1.png",
             bytes.as_bytes(),
             dimensions.0,
             dimensions.1,

src/filter2d/avx/mod.rs

@@ -28,8 +28,8 @@
  */
 mod convolve_op_u8_i16_fp;
 #[cfg(feature = "fft")]
-mod mul_spectrum;
+mod mul_spectrum_f32;
 
 pub(crate) use convolve_op_u8_i16_fp::convolve_segment_sse_2d_u8_i16_fp;
 #[cfg(feature = "fft")]
-pub(crate) use mul_spectrum::avx_fma_mul_spectrum_in_place_f32;
+pub(crate) use mul_spectrum_f32::avx_fma_mul_spectrum_in_place_f32;

src/filter2d/avx/mul_spectrum.rs src/filter2d/avx/mul_spectrum_f32.rssrc/filter2d/avx/mul_spectrum.rs renamed to src/filter2d/avx/mul_spectrum_f32.rs

@@ -65,13 +65,13 @@ unsafe fn avx_interleave(a: __m256, b: __m256) -> (__m256, __m256) {
 }
 
 #[inline]
-#[target_feature(enable = "avx2")]
-unsafe fn sse_unpacklo_ps(a: __m128i) -> (__m128, __m128) {
-    let v2 = _mm_unpacklo_epi32(a, _mm_setzero_si128()); // a0 a2 b0 b2
-
-    let a = _mm_unpacklo_epi32(v2, _mm_setzero_si128()); // a0 a1 a2 a3
-    let b = _mm_unpackhi_epi32(v2, _mm_setzero_si128()); // b0 b1 ab b3
-    (_mm_castsi128_ps(a), _mm_castsi128_ps(b))
+#[target_feature(enable = "avx2", enable = "fma")]
+unsafe fn complex_mul_fma(a: __m128, b: __m128) -> __m128 {
+    let temp1 = _mm_shuffle_ps::<0xA0>(b, b);
+    let temp2 = _mm_shuffle_ps::<0xF5>(b, b);
+    let mul2 = _mm_mul_ps(a, temp2);
+    let mul2 = _mm_shuffle_ps::<0xB1>(mul2, mul2);
+    _mm_fmaddsub_ps(a, temp1, mul2)
 }
 
 #[target_feature(enable = "avx2", enable = "fma")]
@@ -89,15 +89,15 @@ unsafe fn mul_spectrum_in_place_f32_impl(
         let other = &other[..complex_size];
 
         for (dst, kernel) in value1.chunks_exact_mut(16).zip(other.chunks_exact(16)) {
-            let vd0 = _mm256_loadu_ps(dst.as_ptr() as *const f32);
-            let vd1 = _mm256_loadu_ps(dst.get_unchecked(4..).as_ptr() as *const f32);
-            let vd2 = _mm256_loadu_ps(dst.get_unchecked(8..).as_ptr() as *const f32);
-            let vd3 = _mm256_loadu_ps(dst.get_unchecked(12..).as_ptr() as *const f32);
+            let vd0 = _mm256_loadu_ps(dst.as_ptr().cast());
+            let vd1 = _mm256_loadu_ps(dst.get_unchecked(4..).as_ptr().cast());
+            let vd2 = _mm256_loadu_ps(dst.get_unchecked(8..).as_ptr().cast());
+            let vd3 = _mm256_loadu_ps(dst.get_unchecked(12..).as_ptr().cast());
 
-            let vk0 = _mm256_loadu_ps(kernel.as_ptr() as *const f32);
-            let vk1 = _mm256_loadu_ps(kernel.get_unchecked(4..).as_ptr() as *const f32);
-            let vk2 = _mm256_loadu_ps(kernel.get_unchecked(8..).as_ptr() as *const f32);
-            let vk3 = _mm256_loadu_ps(kernel.get_unchecked(12..).as_ptr() as *const f32);
+            let vk0 = _mm256_loadu_ps(kernel.as_ptr().cast());
+            let vk1 = _mm256_loadu_ps(kernel.get_unchecked(4..).as_ptr().cast());
+            let vk2 = _mm256_loadu_ps(kernel.get_unchecked(8..).as_ptr().cast());
+            let vk3 = _mm256_loadu_ps(kernel.get_unchecked(12..).as_ptr().cast());
 
             let (ar0, ai0) = avx_deinterleave(vd0, vd1);
             let (ar1, ai1) = avx_deinterleave(vd2, vd3);
@@ -123,18 +123,18 @@ unsafe fn mul_spectrum_in_place_f32_impl(
             let (d0, d1) = avx_interleave(prod_r0, prod_i0);
             let (d2, d3) = avx_interleave(prod_r1, prod_i1);
 
-            _mm256_storeu_ps(dst.as_mut_ptr() as *mut f32, d0);
-            _mm256_storeu_ps(dst.get_unchecked_mut(4..).as_mut_ptr() as *mut f32, d1);
-            _mm256_storeu_ps(dst.get_unchecked_mut(8..).as_mut_ptr() as *mut f32, d2);
-            _mm256_storeu_ps(dst.get_unchecked_mut(12..).as_mut_ptr() as *mut f32, d3);
+            _mm256_storeu_ps(dst.as_mut_ptr().cast(), d0);
+            _mm256_storeu_ps(dst.get_unchecked_mut(4..).as_mut_ptr().cast(), d1);
+            _mm256_storeu_ps(dst.get_unchecked_mut(8..).as_mut_ptr().cast(), d2);
+            _mm256_storeu_ps(dst.get_unchecked_mut(12..).as_mut_ptr().cast(), d3);
         }
 
         let dst_rem = value1.chunks_exact_mut(16).into_remainder();
         let src_rem = other.chunks_exact(16).remainder();
 
         for (dst, kernel) in dst_rem.chunks_exact_mut(4).zip(src_rem.chunks_exact(4)) {
-            let a0 = _mm256_loadu_ps(dst.as_ptr() as *const f32);
-            let b0 = _mm256_loadu_ps(kernel.as_ptr() as *const f32);
+            let a0 = _mm256_loadu_ps(dst.as_ptr().cast());
+            let b0 = _mm256_loadu_ps(kernel.as_ptr().cast());
 
             let (ar0, ai0) = avx_deinterleave(a0, _mm256_setzero_ps());
             let (br0, bi0) = avx_deinterleave(b0, _mm256_setzero_ps());
@@ -149,7 +149,7 @@ unsafe fn mul_spectrum_in_place_f32_impl(
 
             let (d0, _) = avx_interleave(prod_r0, prod_i0);
 
-            _mm256_storeu_ps(dst.as_mut_ptr() as *mut f32, d0);
+            _mm256_storeu_ps(dst.as_mut_ptr().cast(), d0);
         }
 
         let dst_rem = dst_rem.chunks_exact_mut(4).into_remainder();
@@ -159,18 +159,7 @@ unsafe fn mul_spectrum_in_place_f32_impl(
             let v0 = _mm_loadu_si64(dst as *const Complex<f32> as *const _);
             let v1 = _mm_loadu_si64(kernel as *const Complex<f32> as *const _);
 
-            let (ar0, ai0) = sse_unpacklo_ps(v0);
-            let (br0, bi0) = sse_unpacklo_ps(v1);
-
-            let mut prod_r0 = _mm_mul_ps(ar0, br0);
-            let mut prod_i0 = _mm_mul_ps(ar0, bi0);
-            prod_r0 = _mm_fnmadd_ps(ai0, bi0, prod_r0);
-            prod_i0 = _mm_fmadd_ps(ai0, br0, prod_i0);
-
-            prod_r0 = _mm_mul_ps(prod_r0, _mm256_castps256_ps128(v_norm_factor));
-            prod_i0 = _mm_mul_ps(prod_i0, _mm256_castps256_ps128(v_norm_factor));
-
-            let lo = _mm_unpacklo_ps(prod_r0, prod_i0);
+            let lo = complex_mul_fma(_mm_castsi128_ps(v0), _mm_castsi128_ps(v1));
 
             _mm_storeu_si64(dst as *mut Complex<f32> as *mut _, _mm_castps_si128(lo));
         }

src/filter2d/mul_spectrum.rs

@@ -72,8 +72,8 @@ impl SpectrumMultiplier<f32> for f32 {
         #[cfg(all(target_arch = "aarch64", feature = "nightly_fcma"))]
         {
             if std::arch::is_aarch64_feature_detected!("fcma") {
-                use crate::filter2d::neon::neon_mul_spectrum_in_place_f32;
-                return neon_mul_spectrum_in_place_f32(value1, other, width, height);
+                use crate::filter2d::neon::fcma_mul_spectrum_in_place_f32;
+                return fcma_mul_spectrum_in_place_f32(value1, other, width, height);
             }
         }
         #[cfg(all(target_arch = "x86_64", feature = "avx"))]
@@ -96,7 +96,15 @@ impl SpectrumMultiplier<f32> for f32 {
                 }
             }
         }
-        mul_spectrum_in_place_impl(value1, other, width, height);
+        #[cfg(all(target_arch = "aarch64", feature = "neon"))]
+        {
+            use crate::filter2d::neon::neon_mul_spectrum_in_place_f32;
+            neon_mul_spectrum_in_place_f32(value1, other, width, height);
+        }
+        #[cfg(not(all(target_arch = "aarch64", feature = "neon")))]
+        {
+            mul_spectrum_in_place_impl(value1, other, width, height);
+        }
     }
 }
 
@@ -111,6 +119,7 @@ impl SpectrumMultiplier<f64> for f64 {
     }
 }
 
+#[allow(dead_code)]
 #[inline(always)]
 fn mul_spectrum_in_place_impl<V: FftNum + Mul<V>>(
     value1: &mut [Complex<V>],

src/filter2d/neon/mod.rs

@@ -30,10 +30,14 @@ mod convolve_op_u8_f32;
 mod convolve_op_u8_i16;
 mod convolve_op_u8_i16_fp;
 #[cfg(all(feature = "nightly_fcma", feature = "fft"))]
-mod mul_spectrum;
+mod mul_spectrum_fcma;
+#[cfg(feature = "fft")]
+mod mul_spectrum_neon;
 
 pub(crate) use convolve_op_u8_f32::convolve_segment_neon_2d_u8_f32;
 pub(crate) use convolve_op_u8_i16::convolve_segment_neon_2d_u8_i16;
 pub(crate) use convolve_op_u8_i16_fp::convolve_segment_neon_2d_u8_i16_fp;
 #[cfg(all(feature = "nightly_fcma", feature = "fft"))]
-pub(crate) use mul_spectrum::neon_mul_spectrum_in_place_f32;
+pub(crate) use mul_spectrum_fcma::fcma_mul_spectrum_in_place_f32;
+#[cfg(feature = "fft")]
+pub(crate) use mul_spectrum_neon::neon_mul_spectrum_in_place_f32;

src/filter2d/neon/mul_spectrum.rs src/filter2d/neon/mul_spectrum_fcma.rssrc/filter2d/neon/mul_spectrum.rs renamed to src/filter2d/neon/mul_spectrum_fcma.rs

@@ -29,7 +29,7 @@
 use rustfft::num_complex::Complex;
 use std::arch::aarch64::*;
 
-pub(crate) fn neon_mul_spectrum_in_place_f32(
+pub(crate) fn fcma_mul_spectrum_in_place_f32(
     value1: &mut [Complex<f32>],
     other: &[Complex<f32>],
     width: usize,
@@ -56,15 +56,15 @@ unsafe fn mul_spectrum_in_place_f32_impl(
         let zero = vdupq_n_f32(0.);
 
         for (dst, kernel) in value1.chunks_exact_mut(8).zip(other.chunks_exact(8)) {
-            let vd0 = vld1q_f32(dst.as_ptr() as *const f32);
-            let vd1 = vld1q_f32(dst.as_ptr().add(2) as *const f32);
-            let vd2 = vld1q_f32(dst.as_ptr().add(4) as *const f32);
-            let vd3 = vld1q_f32(dst.as_ptr().add(6) as *const f32);
+            let vd0 = vld1q_f32(dst.as_ptr().cast());
+            let vd1 = vld1q_f32(dst.as_ptr().add(2).cast());
+            let vd2 = vld1q_f32(dst.as_ptr().add(4).cast());
+            let vd3 = vld1q_f32(dst.as_ptr().add(6).cast());
 
-            let vk0 = vld1q_f32(kernel.as_ptr() as *const f32);
-            let vk1 = vld1q_f32(kernel.as_ptr().add(2) as *const f32);
-            let vk2 = vld1q_f32(kernel.as_ptr().add(4) as *const f32);
-            let vk3 = vld1q_f32(kernel.as_ptr().add(6) as *const f32);
+            let vk0 = vld1q_f32(kernel.as_ptr().cast());
+            let vk1 = vld1q_f32(kernel.as_ptr().add(2).cast());
+            let vk2 = vld1q_f32(kernel.as_ptr().add(4).cast());
+            let vk3 = vld1q_f32(kernel.as_ptr().add(6).cast());
 
             let p0 = vmulq_f32(
                 vcmlaq_rot90_f32(vcmlaq_f32(zero, vd0, vk0), vd0, vk0),
@@ -83,21 +83,21 @@ unsafe fn mul_spectrum_in_place_f32_impl(
                 v_norm_factor,
             );
 
-            vst1q_f32(dst.as_mut_ptr() as *mut f32, p0);
-            vst1q_f32(dst.get_unchecked_mut(2..).as_mut_ptr() as *mut f32, p1);
-            vst1q_f32(dst.get_unchecked_mut(4..).as_mut_ptr() as *mut f32, p2);
-            vst1q_f32(dst.get_unchecked_mut(6..).as_mut_ptr() as *mut f32, p3);
+            vst1q_f32(dst.as_mut_ptr().cast(), p0);
+            vst1q_f32(dst.get_unchecked_mut(2..).as_mut_ptr().cast(), p1);
+            vst1q_f32(dst.get_unchecked_mut(4..).as_mut_ptr().cast(), p2);
+            vst1q_f32(dst.get_unchecked_mut(6..).as_mut_ptr().cast(), p3);
         }
 
         let dst_rem = value1.chunks_exact_mut(8).into_remainder();
         let src_rem = other.chunks_exact(8).remainder();
 
         for (dst, kernel) in dst_rem.chunks_exact_mut(2).zip(src_rem.chunks_exact(2)) {
-            let v0 = vld1q_f32(dst.as_ptr() as *const f32);
-            let v1 = vld1q_f32(kernel.as_ptr() as *const f32);
+            let v0 = vld1q_f32(dst.as_ptr().cast());
+            let v1 = vld1q_f32(kernel.as_ptr().cast());
             let p0 = vcmlaq_rot90_f32(vcmlaq_f32(zero, v0, v1), v0, v1);
             let p1 = vmulq_f32(p0, v_norm_factor);
-            vst1q_f32(dst.as_mut_ptr() as *mut f32, p1);
+            vst1q_f32(dst.as_mut_ptr().cast(), p1);
         }
 
         let dst_rem = dst_rem.chunks_exact_mut(2).into_remainder();