Added full optimisation for convolution of U8 images with helps of SSE4.1 instructions.

CHANGELOG.md

@@ -7,6 +7,8 @@
 - Added support of optimisation with helps of `NEON SIMD` for convolution of `U16x4` images.
 - Added optimisation for processing `U16x4` images by `MulDiv` with
   helps of `NEON SIMD` instructions.
+- Added full optimisation for convolution of `U8` images with helps of
+  `SSE4.1` instructions.
 - Fixed link to documentation page in `README.md` file.
 - Fixed error in implementation of `MulDiv::divide_alpha()` and `MulDiv::divide_alpha_inplace()`
   for `U16x4` pixels with optimisation with helps of `SSE4.1` and `AVX2`.

README.md

@@ -10,18 +10,18 @@ Rust library for fast image resizing with using of SIMD instructions.
 
 Supported pixel formats and available optimisations:
 
-| Format | Description                                                   | Native Rust | SSE4.1  | AVX2 | Neon |
-|:------:|:--------------------------------------------------------------|:-----------:|:-------:|:----:|:----:|
-|   U8   | One `u8` component per pixel (e.g. L)                         |      +      | partial |  +   |  -   |
-|  U8x2  | Two `u8` components per pixel (e.g. LA)                       |      +      |    +    |  +   |  -   |
-|  U8x3  | Three `u8` components per pixel (e.g. RGB)                    |      +      |    +    |  +   |  -   |
-|  U8x4  | Four `u8` components per pixel (e.g. RGBA, RGBx, CMYK)        |      +      |    +    |  +   |  +   |
-|  U16   | One `u16` components per pixel (e.g. L16)                     |      +      |    +    |  +   |  -   |
-| U16x2  | Two `u16` components per pixel (e.g. LA16)                    |      +      |    +    |  +   |  -   |
-| 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                                 |      +      |    -    |  -   |  -   |
-|  F32   | One `f32` component per pixel                                 |      +      |    -    |  -   |  -   |
+| Format | Description                                                   | Native Rust | SSE4.1 | AVX2 | Neon |
+|:------:|:--------------------------------------------------------------|:-----------:|:------:|:----:|:----:|
+|   U8   | One `u8` component per pixel (e.g. L)                         |      +      |   +    |  +   |  -   |
+|  U8x2  | Two `u8` components per pixel (e.g. LA)                       |      +      |   +    |  +   |  -   |
+|  U8x3  | Three `u8` components per pixel (e.g. RGB)                    |      +      |   +    |  +   |  -   |
+|  U8x4  | Four `u8` components per pixel (e.g. RGBA, RGBx, CMYK)        |      +      |   +    |  +   |  +   |
+|  U16   | One `u16` components per pixel (e.g. L16)                     |      +      |   +    |  +   |  -   |
+| U16x2  | Two `u16` components per pixel (e.g. LA16)                    |      +      |   +    |  +   |  -   |
+| 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                                 |      +      |   -    |  -   |  -   |
+|  F32   | One `f32` component per pixel                                 |      +      |   -    |  -   |  -   |
 
 ## Colorspace
 

src/convolution/u8x1/avx2.rs

@@ -62,7 +62,7 @@ unsafe fn horiz_convolution_8u4x(
     for (dst_x, coeffs_chunk) in coefficients_chunks.iter().enumerate() {
         let coeffs = coeffs_chunk.values;
         let mut x = coeffs_chunk.start as usize;
-        let mut result_i32x8 = [initial, initial, initial, initial];
+        let mut result_i32x8x4 = [initial, initial, initial, initial];
 
         let coeffs_by_16 = coeffs.chunks_exact(16);
         let reminder16 = coeffs_by_16.remainder();
@@ -71,8 +71,8 @@ unsafe fn horiz_convolution_8u4x(
             for i in 0..4 {
                 let pixels_u8x16 = simd_utils::loadu_si128(s_rows[i], x);
                 let pixels_i16x16 = _mm256_cvtepu8_epi16(pixels_u8x16);
-                result_i32x8[i] = _mm256_add_epi32(
-                    result_i32x8[i],
+                result_i32x8x4[i] = _mm256_add_epi32(
+                    result_i32x8x4[i],
                     _mm256_madd_epi16(pixels_i16x16, coeffs_i16x16),
                 );
             }
@@ -86,15 +86,15 @@ unsafe fn horiz_convolution_8u4x(
             for i in 0..4 {
                 let pixels_u8x8 = simd_utils::loadl_epi64(s_rows[i], x);
                 let pixels_i16x8 = _mm_cvtepu8_epi16(pixels_u8x8);
-                result_i32x8[i] = _mm256_add_epi32(
-                    result_i32x8[i],
+                result_i32x8x4[i] = _mm256_add_epi32(
+                    result_i32x8x4[i],
                     _mm256_set_m128i(zero, _mm_madd_epi16(pixels_i16x8, coeffs_i16x8)),
                 );
             }
             x += 8;
         }
 
-        let mut result_i32x4 = result_i32x8.map(|v| hsum_i32x8_avx2(v));
+        let mut result_i32x4 = result_i32x8x4.map(|v| hsum_i32x8_avx2(v));
 
         for &coeff in reminder8 {
             let coeff_i32 = coeff as i32;
@@ -180,7 +180,8 @@ unsafe fn hsum_i32x8_avx2(v: __m256i) -> i32 {
     hsum_epi32_avx(sum128)
 }
 
-#[inline(always)]
+#[inline]
+#[target_feature(enable = "avx2")]
 unsafe fn hsum_epi32_avx(x: __m128i) -> i32 {
     // 3-operand non-destructive AVX lets us save a byte without needing a movdqa
     let hi64 = _mm_unpackhi_epi64(x, x);

src/convolution/u8x1/mod.rs

@@ -8,6 +8,8 @@ use super::{Coefficients, Convolution};
 #[cfg(target_arch = "x86_64")]
 mod avx2;
 mod native;
+#[cfg(target_arch = "x86_64")]
+mod sse4;
 
 impl Convolution for U8 {
     fn horiz_convolution(
@@ -20,6 +22,8 @@ impl Convolution for U8 {
         match cpu_extensions {
             #[cfg(target_arch = "x86_64")]
             CpuExtensions::Avx2 => avx2::horiz_convolution(src_image, dst_image, offset, coeffs),
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Sse4_1 => sse4::horiz_convolution(src_image, dst_image, offset, coeffs),
             _ => native::horiz_convolution(src_image, dst_image, offset, coeffs),
         }
     }

src/convolution/u8x1/sse4.rs

@@ -0,0 +1,166 @@
+use std::arch::x86_64::*;
+
+use crate::convolution::{optimisations, Coefficients};
+use crate::image_view::{FourRows, FourRowsMut};
+use crate::pixels::U8;
+use crate::simd_utils;
+use crate::{ImageView, ImageViewMut};
+
+#[inline]
+pub(crate) fn horiz_convolution(
+    src_image: &ImageView<U8>,
+    dst_image: &mut ImageViewMut<U8>,
+    offset: u32,
+    coeffs: Coefficients,
+) {
+    let normalizer = optimisations::Normalizer16::new(coeffs);
+    let coefficients_chunks = normalizer.normalized_chunks();
+    let dst_height = dst_image.height().get();
+
+    let src_iter = src_image.iter_4_rows(offset, dst_height + offset);
+    let dst_iter = dst_image.iter_4_rows_mut();
+    for (src_rows, dst_rows) in src_iter.zip(dst_iter) {
+        unsafe {
+            horiz_convolution_four_rows(src_rows, dst_rows, &coefficients_chunks, &normalizer);
+        }
+    }
+
+    let mut yy = dst_height - dst_height % 4;
+    while yy < dst_height {
+        unsafe {
+            horiz_convolution_row(
+                src_image.get_row(yy + offset).unwrap(),
+                dst_image.get_row_mut(yy).unwrap(),
+                &coefficients_chunks,
+                &normalizer,
+            );
+        }
+        yy += 1;
+    }
+}
+
+/// 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
+#[inline]
+#[target_feature(enable = "sse4.1")]
+unsafe fn horiz_convolution_four_rows(
+    src_rows: FourRows<U8>,
+    dst_rows: FourRowsMut<U8>,
+    coefficients_chunks: &[optimisations::CoefficientsI16Chunk],
+    normalizer: &optimisations::Normalizer16,
+) {
+    let s_rows = [src_rows.0, src_rows.1, src_rows.2, src_rows.3];
+    let d_rows = [dst_rows.0, dst_rows.1, dst_rows.2, dst_rows.3];
+    let zero = _mm_setzero_si128();
+    let initial = 1 << (normalizer.precision() - 1);
+    let mut buf = [0, 0, 0, 0, initial];
+
+    for (dst_x, coeffs_chunk) in coefficients_chunks.iter().enumerate() {
+        let coeffs = coeffs_chunk.values;
+        let mut x = coeffs_chunk.start as usize;
+        let mut result_i32x4 = [zero, zero, zero, zero];
+
+        let coeffs_by_8 = coeffs.chunks_exact(8);
+        let reminder8 = coeffs_by_8.remainder();
+        for k in coeffs_by_8 {
+            let coeffs_i16x8 = _mm_loadu_si128(k.as_ptr() as *const __m128i);
+            for i in 0..4 {
+                let pixels_u8x8 = simd_utils::loadl_epi64(s_rows[i], x);
+                let pixels_i16x8 = _mm_cvtepu8_epi16(pixels_u8x8);
+                result_i32x4[i] =
+                    _mm_add_epi32(result_i32x4[i], _mm_madd_epi16(pixels_i16x8, coeffs_i16x8));
+            }
+            x += 8;
+        }
+
+        let mut coeffs_by_4 = reminder8.chunks_exact(4);
+        let reminder4 = coeffs_by_4.remainder();
+        if let Some(k) = coeffs_by_4.next() {
+            let coeffs_i16x4 = simd_utils::loadl_epi64(k, 0);
+            for i in 0..4 {
+                let pixels_u8x4 = simd_utils::loadl_epi32(s_rows[i], x);
+                let pixels_i16x4 = _mm_cvtepu8_epi16(pixels_u8x4);
+                result_i32x4[i] =
+                    _mm_add_epi32(result_i32x4[i], _mm_madd_epi16(pixels_i16x4, coeffs_i16x4));
+            }
+            x += 4;
+        }
+
+        let mut result_i32x4 = result_i32x4.map(|v| {
+            _mm_storeu_si128(buf.as_mut_ptr() as *mut __m128i, v);
+            buf.iter().sum()
+        });
+
+        for &coeff in reminder4 {
+            let coeff_i32 = coeff as i32;
+            for i in 0..4 {
+                result_i32x4[i] += s_rows[i].get_unchecked(x).0.to_owned() as i32 * coeff_i32;
+            }
+            x += 1;
+        }
+
+        let result_u8x4 = result_i32x4.map(|v| normalizer.clip(v));
+        for i in 0..4 {
+            d_rows[i].get_unchecked_mut(dst_x).0 = result_u8x4[i];
+        }
+    }
+}
+
+/// For safety, it is necessary to ensure the following conditions:
+/// - bounds.len() == dst_row.len()
+/// - coeffs.len() == dst_rows.0.len() * window_size
+/// - max(bound.start + bound.size for bound in bounds) <= src_row.len()
+/// - precision <= MAX_COEFS_PRECISION
+#[inline]
+#[target_feature(enable = "sse4.1")]
+unsafe fn horiz_convolution_row(
+    src_row: &[U8],
+    dst_row: &mut [U8],
+    coefficients_chunks: &[optimisations::CoefficientsI16Chunk],
+    normalizer: &optimisations::Normalizer16,
+) {
+    let zero = _mm_setzero_si128();
+    let initial = 1 << (normalizer.precision() - 1);
+    let mut buf = [0, 0, 0, 0, initial];
+
+    for (dst_x, &coeffs_chunk) in coefficients_chunks.iter().enumerate() {
+        let coeffs = coeffs_chunk.values;
+        let mut x = coeffs_chunk.start as usize;
+        let mut result_i32x4 = zero;
+
+        let coeffs_by_8 = coeffs.chunks_exact(8);
+        let reminder8 = coeffs_by_8.remainder();
+        for k in coeffs_by_8 {
+            let coeffs_i16x8 = _mm_loadu_si128(k.as_ptr() as *const __m128i);
+            let pixels_u8x8 = simd_utils::loadl_epi64(src_row, x);
+            let pixels_i16x8 = _mm_cvtepu8_epi16(pixels_u8x8);
+            result_i32x4 = _mm_add_epi32(result_i32x4, _mm_madd_epi16(pixels_i16x8, coeffs_i16x8));
+            x += 8;
+        }
+
+        let mut coeffs_by_4 = reminder8.chunks_exact(4);
+        let reminder4 = coeffs_by_4.remainder();
+        if let Some(k) = coeffs_by_4.next() {
+            let coeffs_i16x4 = simd_utils::loadl_epi64(k, 0);
+            let pixels_u8x4 = simd_utils::loadl_epi32(src_row, x);
+            let pixels_i16x4 = _mm_cvtepu8_epi16(pixels_u8x4);
+            result_i32x4 = _mm_add_epi32(result_i32x4, _mm_madd_epi16(pixels_i16x4, coeffs_i16x4));
+            x += 4;
+        }
+
+        _mm_storeu_si128(buf.as_mut_ptr() as *mut __m128i, result_i32x4);
+        let mut result_i32 = buf.iter().sum();
+
+        for &coeff in reminder4 {
+            let coeff_i32 = coeff as i32;
+            result_i32 += src_row.get_unchecked(x).0 as i32 * coeff_i32;
+            x += 1;
+        }
+
+        dst_row.get_unchecked_mut(dst_x).0 = normalizer.clip(result_i32);
+    }
+}