Added support for optimization with help of SSE4.1 and AVX2 for the F32 pixel type.

CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ### Added
 
+- Added support for optimization with help of `SSE4.1` and `AVX2` for
+  the `F32` pixel type.
 - Added support for the new pixel type `PixelType::F32x2` with
   optimizations for SSE4.1 and AVX2 (#30).
 - Added basic support for the new pixel type `PixelType::F32x3` (#30).

Cargo.toml

@@ -120,6 +120,11 @@ name = "bench_compare_la16"
 harness = false
 
 
+[[bench]]
+name = "bench_compare_l32f"
+harness = false
+
+
 [[bench]]
 name = "bench_compare_la32f"
 harness = false

README.md

@@ -8,7 +8,7 @@ Rust library for fast image resizing with using of SIMD instructions.
 
 [CHANGELOG](https://github.com/Cykooz/fast_image_resize/blob/main/CHANGELOG.md)
 
-Supported pixel formats and available optimisations:
+Supported pixel formats and available optimizations:
 
 | Format | Description                                                   | SSE4.1 | AVX2 | Neon | Wasm32 SIMD128 |
 |:------:|:--------------------------------------------------------------|:------:|:----:|:----:|:--------------:|
@@ -21,9 +21,10 @@ Supported pixel formats and available optimisations:
 | U16x3  | Three `u16` components per pixel (e.g. RGB16)                 |   +    |  +   |  +   |       +        |
 | U16x4  | Four `u16` components per pixel (e.g. RGBA16, RGBx16, CMYK16) |   +    |  +   |  +   |       +        |
 |  I32   | One `i32` component per pixel (e.g. L)                        |   -    |  -   |  -   |       -        |
-|  F32   | One `f32` component per pixel (e.g. L)                        |   -    |  -   |  -   |       -        |
+|  F32   | One `f32` component per pixel (e.g. L)                        |   +    |  +   |  -   |       -        |
 | F32x2  | Two `f32` components per pixel (e.g. LA32F)                   |   +    |  +   |  -   |       -        |
 | F32x3  | Three `f32` components per pixel (e.g. RGB32F)                |   -    |  -   |  -   |       -        |
+| F32x4  | Four `f32` components per pixel (e.g. RGBA32F)                |   -    |  -   |  -   |       -        |
 
 ## Colorspace
 

benches/bench_compare_l32f.rs

@@ -0,0 +1,27 @@
+use resize::Pixel::GrayF32;
+use rgb::FromSlice;
+
+use fast_image_resize::pixels::F32;
+use testing::PixelTestingExt;
+
+mod utils;
+
+pub fn bench_downscale_l32f(bench_group: &mut utils::BenchGroup) {
+    type P = F32;
+    let src_image = P::load_big_image();
+    utils::image_resize(bench_group, &src_image);
+    utils::resize_resize(
+        bench_group,
+        GrayF32,
+        src_image.as_raw().as_gray(),
+        src_image.width(),
+        src_image.height(),
+    );
+    utils::libvips_resize::<P>(bench_group, false);
+    utils::fir_resize::<P>(bench_group, false);
+}
+
+fn main() {
+    let res = utils::run_bench(bench_downscale_l32f, "Compare resize of F32 image");
+    utils::print_and_write_compare_result(&res);
+}

benches/bench_resize.rs

@@ -84,6 +84,10 @@ pub fn resize_in_one_dimension_bench(bench_group: &mut utils::BenchGroup) {
         PixelType::U16x2,
         PixelType::U16x3,
         PixelType::U16x4,
+        PixelType::F32,
+        PixelType::F32x2,
+        PixelType::F32x3,
+        PixelType::F32x4,
     ];
     let mut cpu_extensions = vec![CpuExtensions::None];
     #[cfg(target_arch = "x86_64")]
@@ -115,7 +119,10 @@ pub fn resize_in_one_dimension_bench(bench_group: &mut utils::BenchGroup) {
                 PixelType::U16x3 => U16x3::load_big_square_src_image(),
                 PixelType::U16x4 => U16x4::load_big_square_src_image(),
                 PixelType::I32 => I32::load_big_square_src_image(),
+                PixelType::F32 => F32::load_big_square_src_image(),
                 PixelType::F32x2 => F32x2::load_big_square_src_image(),
+                PixelType::F32x3 => F32x3::load_big_square_src_image(),
+                PixelType::F32x4 => F32x4::load_big_square_src_image(),
                 _ => unreachable!(),
             };
             downscale_bench(
@@ -151,7 +158,10 @@ pub fn resize_bench(bench_group: &mut utils::BenchGroup) {
         PixelType::U16x3,
         PixelType::U16x4,
         PixelType::I32,
+        PixelType::F32,
         PixelType::F32x2,
+        PixelType::F32x3,
+        PixelType::F32x4,
     ];
     let mut cpu_extensions = vec![CpuExtensions::None];
     #[cfg(target_arch = "x86_64")]
@@ -183,7 +193,10 @@ pub fn resize_bench(bench_group: &mut utils::BenchGroup) {
                 PixelType::U16x3 => U16x3::load_big_square_src_image(),
                 PixelType::U16x4 => U16x4::load_big_square_src_image(),
                 PixelType::I32 => I32::load_big_square_src_image(),
+                PixelType::F32 => F32::load_big_square_src_image(),
                 PixelType::F32x2 => F32x2::load_big_square_src_image(),
+                PixelType::F32x3 => F32x3::load_big_square_src_image(),
+                PixelType::F32x4 => F32x4::load_big_square_src_image(),
                 _ => unreachable!(),
             };
             downscale_bench(

benches/templates/bench_compare_l32f.md.tera

@@ -0,0 +1,11 @@
+### Resize L32F image (F32) 4928x3279 => 852x567
+
+Pipeline:
+
+`src_image => resize => dst_image`
+
+- Source image [nasa-4928x3279.png](https://github.com/Cykooz/fast_image_resize/blob/main/data/nasa-4928x3279.png)
+  has converted into grayscale image with two bytes per pixel.
+- Numbers in the table mean a duration of image resizing in milliseconds.
+
+{{ compare_results -}}

benchmarks-x86_64.md

@@ -212,6 +212,29 @@ Pipeline:
 
 <!-- bench_compare_la16 end -->
 
+<!-- bench_compare_l32f start -->
+
+### Resize L32F image (F32) 4928x3279 => 852x567
+
+Pipeline:
+
+`src_image => resize => dst_image`
+
+- Source image [nasa-4928x3279.png](https://github.com/Cykooz/fast_image_resize/blob/main/data/nasa-4928x3279.png)
+  has converted into grayscale image with two bytes per pixel.
+- Numbers in the table mean a duration of image resizing in milliseconds.
+
+|            | Nearest |  Box  | Bilinear | Bicubic | Lanczos3 |
+|------------|:-------:|:-----:|:--------:|:-------:|:--------:|
+| image      |  24.16  |   -   |  52.33   |  82.88  |  111.10  |
+| resize     |  4.99   | 8.85  |  13.38   |  30.00  |  45.71   |
+| libvips    |  7.37   | 25.94 |  20.30   |  40.97  |  70.76   |
+| fir rust   |  0.18   | 9.56  |  15.05   |  29.59  |  52.03   |
+| fir sse4.1 |  0.18   | 5.03  |   7.30   |  11.55  |  16.90   |
+| fir avx2   |  0.18   | 4.65  |   5.41   |  7.14   |  10.78   |
+
+<!-- bench_compare_l32f end -->
+
 <!-- bench_compare_la32f start -->
 
 ### Resize LA32F (luma with alpha channel) image (F32x2) 4928x3279 => 852x567
@@ -232,7 +255,7 @@ Pipeline:
 | libvips    |  11.85  | 70.31 |  101.80  | 177.59  |  254.22  |
 | fir rust   |  0.38   | 21.26 |  28.82   |  47.50  |  70.34   |
 | fir sse4.1 |  0.38   | 16.23 |  20.91   |  30.35  |  40.12   |
-| fir avx2   |  0.39   | 15.05 |  17.18   |  22.47  |  27.89   |
+| fir avx2   |  0.38   | 15.05 |  17.18   |  22.47  |  27.89   |
 
 <!-- bench_compare_la32f end -->
 

src/convolution/f32x1/avx2.rs

@@ -0,0 +1,127 @@
+use std::arch::x86_64::*;
+
+use crate::convolution::{Coefficients, CoefficientsChunk};
+use crate::pixels::F32;
+use crate::{simd_utils, ImageView, ImageViewMut};
+
+#[inline]
+pub(crate) fn horiz_convolution(
+    src_view: &impl ImageView<Pixel = F32>,
+    dst_view: &mut impl ImageViewMut<Pixel = F32>,
+    offset: u32,
+    coeffs: Coefficients,
+) {
+    let coefficients_chunks = coeffs.get_chunks();
+    let dst_height = dst_view.height();
+
+    let src_iter = src_view.iter_4_rows(offset, dst_height + offset);
+    let dst_iter = dst_view.iter_4_rows_mut();
+    for (src_rows, dst_rows) in src_iter.zip(dst_iter) {
+        unsafe {
+            horiz_convolution_rows(src_rows, dst_rows, &coefficients_chunks);
+        }
+    }
+
+    let yy = dst_height - dst_height % 4;
+    let src_rows = src_view.iter_rows(yy + offset);
+    let dst_rows = dst_view.iter_rows_mut(yy);
+    for (src_row, dst_row) in src_rows.zip(dst_rows) {
+        unsafe {
+            horiz_convolution_rows([src_row], [dst_row], &coefficients_chunks);
+        }
+    }
+}
+
+/// For safety, it is necessary to ensure the following conditions:
+/// - length of all rows in src_rows must be equal
+/// - length of all rows in dst_rows must be equal
+/// - coefficients_chunks.len() == dst_rows.0.len()
+/// - max(chunk.start + chunk.values.len() for chunk in coefficients_chunks) <= src_row.0.len()
+/// - precision <= MAX_COEFS_PRECISION
+#[target_feature(enable = "avx2")]
+unsafe fn horiz_convolution_rows<const ROWS_COUNT: usize>(
+    src_rows: [&[F32]; ROWS_COUNT],
+    dst_rows: [&mut [F32]; ROWS_COUNT],
+    coefficients_chunks: &[CoefficientsChunk],
+) {
+    let mut ll_buf = [0f64; 2];
+
+    for (dst_x, coeffs_chunk) in coefficients_chunks.iter().enumerate() {
+        let mut x: usize = coeffs_chunk.start as usize;
+        let mut sums = [_mm256_set1_pd(0.); ROWS_COUNT];
+
+        let mut coeffs = coeffs_chunk.values;
+
+        let coeffs_by_8 = coeffs.chunks_exact(8);
+        coeffs = coeffs_by_8.remainder();
+
+        for k in coeffs_by_8 {
+            let coeff03_f64x4 = simd_utils::loadu_pd256(k, 0);
+            let coeff47_f64x4 = simd_utils::loadu_pd256(k, 4);
+
+            for i in 0..ROWS_COUNT {
+                let mut sum = sums[i];
+                let source = simd_utils::loadu_ps256(src_rows[i], x);
+
+                let pixels03_f64x4 = _mm256_cvtps_pd(_mm256_extractf128_ps::<0>(source));
+                sum = _mm256_add_pd(sum, _mm256_mul_pd(pixels03_f64x4, coeff03_f64x4));
+
+                let pixels47_f64x4 = _mm256_cvtps_pd(_mm256_extractf128_ps::<1>(source));
+                sum = _mm256_add_pd(sum, _mm256_mul_pd(pixels47_f64x4, coeff47_f64x4));
+
+                sums[i] = sum;
+            }
+            x += 8;
+        }
+
+        let coeffs_by_4 = coeffs.chunks_exact(4);
+        coeffs = coeffs_by_4.remainder();
+
+        for k in coeffs_by_4 {
+            let coeff03_f64x4 = simd_utils::loadu_pd256(k, 0);
+
+            for i in 0..ROWS_COUNT {
+                let mut sum = sums[i];
+                let source = simd_utils::loadu_ps(src_rows[i], x);
+                let pixels03_f64x4 = _mm256_cvtps_pd(source);
+                sum = _mm256_add_pd(sum, _mm256_mul_pd(pixels03_f64x4, coeff03_f64x4));
+                sums[i] = sum;
+            }
+            x += 4;
+        }
+
+        let coeffs_by_2 = coeffs.chunks_exact(2);
+        coeffs = coeffs_by_2.remainder();
+        for k in coeffs_by_2 {
+            let coeff01_f64x4 = _mm256_set_pd(0., 0., k[1], k[0]);
+
+            for i in 0..ROWS_COUNT {
+                let pixel0 = src_rows[i].get_unchecked(x).0;
+                let pixel1 = src_rows[i].get_unchecked(x + 1).0;
+                let pixel01_f64x4 = _mm256_set_pd(0., 0., pixel1 as f64, pixel0 as f64);
+                sums[i] = _mm256_add_pd(sums[i], _mm256_mul_pd(pixel01_f64x4, coeff01_f64x4));
+            }
+            x += 2;
+        }
+
+        if let Some(&k) = coeffs.first() {
+            let coeff0_f64x4 = _mm256_set_pd(0., 0., 0., k);
+
+            for i in 0..ROWS_COUNT {
+                let pixel0 = src_rows[i].get_unchecked(x).0;
+                let pixel0_f64x4 = _mm256_set_pd(0., 0., 0., pixel0 as f64);
+                sums[i] = _mm256_add_pd(sums[i], _mm256_mul_pd(pixel0_f64x4, coeff0_f64x4));
+            }
+        }
+
+        for i in 0..ROWS_COUNT {
+            let sum_f64x2 = _mm_add_pd(
+                _mm256_extractf128_pd::<0>(sums[i]),
+                _mm256_extractf128_pd::<1>(sums[i]),
+            );
+            _mm_storeu_pd(ll_buf.as_mut_ptr(), sum_f64x2);
+            let dst_pixel = dst_rows[i].get_unchecked_mut(dst_x);
+            dst_pixel.0 = (ll_buf[0] + ll_buf[1]) as f32;
+        }
+    }
+}

src/convolution/f32x1/mod.rs

@@ -5,7 +5,15 @@ use crate::{ImageView, ImageViewMut};
 
 use super::{Coefficients, Convolution};
 
+#[cfg(target_arch = "x86_64")]
+mod avx2;
 mod native;
+// #[cfg(target_arch = "aarch64")]
+// mod neon;
+#[cfg(target_arch = "x86_64")]
+mod sse4;
+// #[cfg(target_arch = "wasm32")]
+// mod wasm32;
 
 impl Convolution for F32 {
     fn horiz_convolution(
@@ -13,9 +21,19 @@ impl Convolution for F32 {
         dst_view: &mut impl ImageViewMut<Pixel = Self>,
         offset: u32,
         coeffs: Coefficients,
-        _cpu_extensions: CpuExtensions,
+        cpu_extensions: CpuExtensions,
     ) {
-        native::horiz_convolution(src_view, dst_view, offset, coeffs);
+        match cpu_extensions {
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Avx2 => avx2::horiz_convolution(src_view, dst_view, offset, coeffs),
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Sse4_1 => sse4::horiz_convolution(src_view, dst_view, offset, coeffs),
+            // #[cfg(target_arch = "aarch64")]
+            // CpuExtensions::Neon => neon::horiz_convolution(src_view, dst_view, offset, coeffs),
+            // #[cfg(target_arch = "wasm32")]
+            // CpuExtensions::Simd128 => wasm32::horiz_convolution(src_view, dst_view, offset, coeffs),
+            _ => native::horiz_convolution(src_view, dst_view, offset, coeffs),
+        }
     }
 
     fn vert_convolution(

src/convolution/f32x1/sse4.rs

@@ -0,0 +1,107 @@
+use std::arch::x86_64::*;
+
+use crate::convolution::{Coefficients, CoefficientsChunk};
+use crate::pixels::F32;
+use crate::{simd_utils, ImageView, ImageViewMut};
+
+#[inline]
+pub(crate) fn horiz_convolution(
+    src_view: &impl ImageView<Pixel = F32>,
+    dst_view: &mut impl ImageViewMut<Pixel = F32>,
+    offset: u32,
+    coeffs: Coefficients,
+) {
+    let coefficients_chunks = coeffs.get_chunks();
+    let dst_height = dst_view.height();
+
+    let src_iter = src_view.iter_4_rows(offset, dst_height + offset);
+    let dst_iter = dst_view.iter_4_rows_mut();
+    for (src_rows, dst_rows) in src_iter.zip(dst_iter) {
+        unsafe {
+            horiz_convolution_rows(src_rows, dst_rows, &coefficients_chunks);
+        }
+    }
+
+    let yy = dst_height - dst_height % 4;
+    let src_rows = src_view.iter_rows(yy + offset);
+    let dst_rows = dst_view.iter_rows_mut(yy);
+    for (src_row, dst_row) in src_rows.zip(dst_rows) {
+        unsafe {
+            horiz_convolution_rows([src_row], [dst_row], &coefficients_chunks);
+        }
+    }
+}
+
+/// For safety, it is necessary to ensure the following conditions:
+/// - length of all rows in src_rows must be equal
+/// - length of all rows in dst_rows must be equal
+/// - coefficients_chunks.len() == dst_rows.0.len()
+/// - max(chunk.start + chunk.values.len() for chunk in coefficients_chunks) <= src_row.0.len()
+/// - precision <= MAX_COEFS_PRECISION
+#[target_feature(enable = "sse4.1")]
+unsafe fn horiz_convolution_rows<const ROWS_COUNT: usize>(
+    src_rows: [&[F32]; ROWS_COUNT],
+    dst_rows: [&mut [F32]; ROWS_COUNT],
+    coefficients_chunks: &[CoefficientsChunk],
+) {
+    let mut ll_buf = [0f64; 2];
+
+    for (dst_x, coeffs_chunk) in coefficients_chunks.iter().enumerate() {
+        let mut x: usize = coeffs_chunk.start as usize;
+        let mut sums = [_mm_set1_pd(0.); ROWS_COUNT];
+
+        let mut coeffs = coeffs_chunk.values;
+
+        let coeffs_by_4 = coeffs.chunks_exact(4);
+        coeffs = coeffs_by_4.remainder();
+
+        for k in coeffs_by_4 {
+            let coeff01_f64x2 = simd_utils::loadu_pd(k, 0);
+            let coeff23_f64x2 = simd_utils::loadu_pd(k, 2);
+
+            for i in 0..ROWS_COUNT {
+                let mut sum = sums[i];
+                let source = simd_utils::loadu_ps(src_rows[i], x);
+
+                let pixel01_f64 = _mm_cvtps_pd(source);
+                sum = _mm_add_pd(sum, _mm_mul_pd(pixel01_f64, coeff01_f64x2));
+
+                let pixel23_f64 = _mm_cvtps_pd(_mm_movehl_ps(source, source));
+                sum = _mm_add_pd(sum, _mm_mul_pd(pixel23_f64, coeff23_f64x2));
+
+                sums[i] = sum;
+            }
+            x += 4;
+        }
+
+        let coeffs_by_2 = coeffs.chunks_exact(2);
+        coeffs = coeffs_by_2.remainder();
+        for k in coeffs_by_2 {
+            let coeff01_f64x2 = simd_utils::loadu_pd(k, 0);
+
+            for i in 0..ROWS_COUNT {
+                let pixel0 = src_rows[i].get_unchecked(x).0;
+                let pixel1 = src_rows[i].get_unchecked(x + 1).0;
+                let pixel01_f64 = _mm_set_pd(pixel1 as f64, pixel0 as f64);
+                sums[i] = _mm_add_pd(sums[i], _mm_mul_pd(pixel01_f64, coeff01_f64x2));
+            }
+            x += 2;
+        }
+
+        if let Some(&k) = coeffs.first() {
+            let coeff0_f64x2 = _mm_set1_pd(k);
+
+            for i in 0..ROWS_COUNT {
+                let pixel0 = src_rows[i].get_unchecked(x).0;
+                let pixel0_f64 = _mm_set_pd(0., pixel0 as f64);
+                sums[i] = _mm_add_pd(sums[i], _mm_mul_pd(pixel0_f64, coeff0_f64x2));
+            }
+        }
+
+        for i in 0..ROWS_COUNT {
+            _mm_storeu_pd(ll_buf.as_mut_ptr(), sums[i]);
+            let dst_pixel = dst_rows[i].get_unchecked_mut(dst_x);
+            dst_pixel.0 = (ll_buf[0] + ll_buf[1]) as f32;
+        }
+    }
+}

src/simd_utils.rs

@@ -40,6 +40,16 @@ pub unsafe fn loadu_ps256<T>(buf: &[T], index: usize) -> __m256 {
     _mm256_loadu_ps(buf.get_unchecked(index..).as_ptr() as *const f32)
 }
 
+#[inline(always)]
+pub unsafe fn loadu_pd<T>(buf: &[T], index: usize) -> __m128d {
+    _mm_loadu_pd(buf.get_unchecked(index..).as_ptr() as *const f64)
+}
+
+#[inline(always)]
+pub unsafe fn loadu_pd256<T>(buf: &[T], index: usize) -> __m256d {
+    _mm256_loadu_pd(buf.get_unchecked(index..).as_ptr() as *const f64)
+}
+
 #[inline(always)]
 pub unsafe fn mm_cvtepu8_epi32(buf: &[U8x4], index: usize) -> __m128i {
     let v: i32 = transmute(buf.get_unchecked(index).0);

testing/src/lib.rs

@@ -318,10 +318,6 @@ impl PixelTestingExt for F32 {
     type ImagePixel = image::Luma<f32>;
     type Container = Vec<f32>;
 
-    fn cpu_extensions() -> Vec<CpuExtensions> {
-        vec![CpuExtensions::None]
-    }
-
     fn load_image_buffer(
         img_reader: Reader<BufReader<File>>,
     ) -> ImageBuffer<Self::ImagePixel, Self::Container> {