- Improved speed of MulDiv implementation for U16x2 images.

- Added optimisation for processing `U16x2` images by `MulDiv` with helps of `NEON SIMD` instructions.

CHANGELOG.md

@@ -2,7 +2,9 @@
 
 ### Crate
 
-- Improved speed of `MulDiv` implementation for `U8x2`, `U8x4` and `U16x4` images.
+- Improved speed of `MulDiv` implementation for `U8x2`, `U8x4`, `U16x2` and `U16x4` images.
+- Added optimisation for processing `U16x2` images by `MulDiv` with
+  helps of `NEON SIMD` instructions.
 - Excluded possibility of unnecessary operations during resize 
   of cropped image by convolution algorithm.
 

src/alpha/u16x2/avx2.rs

@@ -1,6 +1,7 @@
 use std::arch::x86_64::*;
 
 use crate::pixels::U16x2;
+use crate::utils::foreach_with_pre_reading;
 use crate::{ImageView, ImageViewMut};
 
 use super::sse4;
@@ -20,15 +21,63 @@ pub(crate) unsafe fn multiply_alpha(
 
 #[target_feature(enable = "avx2")]
 pub(crate) unsafe fn multiply_alpha_inplace(image: &mut ImageViewMut<U16x2>) {
-    for dst_row in image.iter_rows_mut() {
-        let src_row = std::slice::from_raw_parts(dst_row.as_ptr(), dst_row.len());
-        multiply_alpha_row(src_row, dst_row);
+    for row in image.iter_rows_mut() {
+        multiply_alpha_row_inplace(row);
     }
 }
 
 #[inline]
 #[target_feature(enable = "avx2")]
 pub(crate) unsafe fn multiply_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2]) {
+    let src_chunks = src_row.chunks_exact(8);
+    let src_remainder = src_chunks.remainder();
+    let mut dst_chunks = dst_row.chunks_exact_mut(8);
+    let src_dst = src_chunks.zip(&mut dst_chunks);
+    foreach_with_pre_reading(
+        src_dst,
+        |(src, dst)| {
+            let pixels = _mm256_loadu_si256(src.as_ptr() as *const __m256i);
+            let dst_ptr = dst.as_mut_ptr() as *mut __m256i;
+            (pixels, dst_ptr)
+        },
+        |(mut pixels, dst_ptr)| {
+            pixels = multiply_alpha_8_pixels(pixels);
+            _mm256_storeu_si256(dst_ptr, pixels);
+        },
+    );
+
+    if !src_remainder.is_empty() {
+        let dst_reminder = dst_chunks.into_remainder();
+        sse4::multiply_alpha_row(src_remainder, dst_reminder);
+    }
+}
+
+#[inline]
+#[target_feature(enable = "avx2")]
+pub(crate) unsafe fn multiply_alpha_row_inplace(row: &mut [U16x2]) {
+    let mut chunks = row.chunks_exact_mut(8);
+    foreach_with_pre_reading(
+        &mut chunks,
+        |chunk| {
+            let pixels = _mm256_loadu_si256(chunk.as_ptr() as *const __m256i);
+            let dst_ptr = chunk.as_mut_ptr() as *mut __m256i;
+            (pixels, dst_ptr)
+        },
+        |(mut pixels, dst_ptr)| {
+            pixels = multiply_alpha_8_pixels(pixels);
+            _mm256_storeu_si256(dst_ptr, pixels);
+        },
+    );
+
+    let reminder = chunks.into_remainder();
+    if !reminder.is_empty() {
+        sse4::multiply_alpha_row_inplace(reminder);
+    }
+}
+
+#[inline]
+#[target_feature(enable = "avx2")]
+unsafe fn multiply_alpha_8_pixels(pixels: __m256i) -> __m256i {
     let zero = _mm256_setzero_si256();
     let half = _mm256_set1_epi32(0x8000);
 
@@ -39,42 +88,27 @@ pub(crate) unsafe fn multiply_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2]
        |0001 0203| |0405 0607| |0809 1011| |1213 1415|
     */
     #[rustfmt::skip]
-    let factor_mask = _mm256_set_epi8(
+        let factor_mask = _mm256_set_epi8(
         15, 14, 15, 14, 11, 10, 11, 10, 7, 6, 7, 6, 3, 2, 3, 2,
         15, 14, 15, 14, 11, 10, 11, 10, 7, 6, 7, 6, 3, 2, 3, 2
     );
 
-    let src_chunks = src_row.chunks_exact(8);
-    let src_remainder = src_chunks.remainder();
-    let mut dst_chunks = dst_row.chunks_exact_mut(8);
+    let factor_pixels = _mm256_shuffle_epi8(pixels, factor_mask);
+    let factor_pixels = _mm256_or_si256(factor_pixels, max_alpha);
 
-    for (src, dst) in src_chunks.zip(&mut dst_chunks) {
-        let src_pixels = _mm256_loadu_si256(src.as_ptr() as *const __m256i);
+    let src_i32_lo = _mm256_unpacklo_epi16(pixels, zero);
+    let factors = _mm256_unpacklo_epi16(factor_pixels, zero);
+    let src_i32_lo = _mm256_add_epi32(_mm256_mullo_epi32(src_i32_lo, factors), half);
+    let dst_i32_lo = _mm256_add_epi32(src_i32_lo, _mm256_srli_epi32::<16>(src_i32_lo));
+    let dst_i32_lo = _mm256_srli_epi32::<16>(dst_i32_lo);
 
-        let factor_pixels = _mm256_shuffle_epi8(src_pixels, factor_mask);
-        let factor_pixels = _mm256_or_si256(factor_pixels, max_alpha);
+    let src_i32_hi = _mm256_unpackhi_epi16(pixels, zero);
+    let factors = _mm256_unpackhi_epi16(factor_pixels, zero);
+    let src_i32_hi = _mm256_add_epi32(_mm256_mullo_epi32(src_i32_hi, factors), half);
+    let dst_i32_hi = _mm256_add_epi32(src_i32_hi, _mm256_srli_epi32::<16>(src_i32_hi));
+    let dst_i32_hi = _mm256_srli_epi32::<16>(dst_i32_hi);
 
-        let src_i32_lo = _mm256_unpacklo_epi16(src_pixels, zero);
-        let factors = _mm256_unpacklo_epi16(factor_pixels, zero);
-        let src_i32_lo = _mm256_add_epi32(_mm256_mullo_epi32(src_i32_lo, factors), half);
-        let dst_i32_lo = _mm256_add_epi32(src_i32_lo, _mm256_srli_epi32::<16>(src_i32_lo));
-        let dst_i32_lo = _mm256_srli_epi32::<16>(dst_i32_lo);
-
-        let src_i32_hi = _mm256_unpackhi_epi16(src_pixels, zero);
-        let factors = _mm256_unpackhi_epi16(factor_pixels, zero);
-        let src_i32_hi = _mm256_add_epi32(_mm256_mullo_epi32(src_i32_hi, factors), half);
-        let dst_i32_hi = _mm256_add_epi32(src_i32_hi, _mm256_srli_epi32::<16>(src_i32_hi));
-        let dst_i32_hi = _mm256_srli_epi32::<16>(dst_i32_hi);
-
-        let dst_pixels = _mm256_packus_epi32(dst_i32_lo, dst_i32_hi);
-
-        _mm256_storeu_si256(dst.as_mut_ptr() as *mut __m256i, dst_pixels);
-    }
-
-    if !src_remainder.is_empty() {
-        let dst_reminder = dst_chunks.into_remainder();
-        sse4::multiply_alpha_row(src_remainder, dst_reminder);
-    }
+    _mm256_packus_epi32(dst_i32_lo, dst_i32_hi)
 }
 
 // Divide
@@ -94,9 +128,8 @@ pub(crate) unsafe fn divide_alpha(
 
 #[target_feature(enable = "avx2")]
 pub(crate) unsafe fn divide_alpha_inplace(image: &mut ImageViewMut<U16x2>) {
-    for dst_row in image.iter_rows_mut() {
-        let src_row = std::slice::from_raw_parts(dst_row.as_ptr(), dst_row.len());
-        divide_alpha_row(src_row, dst_row);
+    for row in image.iter_rows_mut() {
+        divide_alpha_row_inplace(row);
     }
 }
 
@@ -105,10 +138,19 @@ pub(crate) unsafe fn divide_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2])
     let src_chunks = src_row.chunks_exact(8);
     let src_remainder = src_chunks.remainder();
     let mut dst_chunks = dst_row.chunks_exact_mut(8);
-
-    for (src, dst) in src_chunks.zip(&mut dst_chunks) {
-        divide_alpha_eight_pixels(src.as_ptr(), dst.as_mut_ptr());
-    }
+    let src_dst = src_chunks.zip(&mut dst_chunks);
+    foreach_with_pre_reading(
+        src_dst,
+        |(src, dst)| {
+            let pixels = _mm256_loadu_si256(src.as_ptr() as *const __m256i);
+            let dst_ptr = dst.as_mut_ptr() as *mut __m256i;
+            (pixels, dst_ptr)
+        },
+        |(mut pixels, dst_ptr)| {
+            pixels = divide_alpha_8_pixels(pixels);
+            _mm256_storeu_si256(dst_ptr, pixels);
+        },
+    );
 
     if !src_remainder.is_empty() {
         let dst_reminder = dst_chunks.into_remainder();
@@ -119,7 +161,9 @@ pub(crate) unsafe fn divide_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2])
             .for_each(|(d, s)| *d = *s);
 
         let mut dst_pixels = [U16x2::new([0, 0]); 8];
-        divide_alpha_eight_pixels(src_pixels.as_ptr(), dst_pixels.as_mut_ptr());
+        let mut pixels = _mm256_loadu_si256(src_pixels.as_ptr() as *const __m256i);
+        pixels = divide_alpha_8_pixels(pixels);
+        _mm256_storeu_si256(dst_pixels.as_mut_ptr() as *mut __m256i, pixels);
 
         dst_pixels
             .iter()
@@ -128,9 +172,36 @@ pub(crate) unsafe fn divide_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2])
     }
 }
 
+#[target_feature(enable = "avx2")]
+pub(crate) unsafe fn divide_alpha_row_inplace(row: &mut [U16x2]) {
+    let mut chunks = row.chunks_exact_mut(8);
+    // Using a simple for-loop in this case is faster than implementation with pre-reading
+    for chunk in &mut chunks {
+        let mut pixels = _mm256_loadu_si256(chunk.as_ptr() as *const __m256i);
+        pixels = divide_alpha_8_pixels(pixels);
+        _mm256_storeu_si256(chunk.as_mut_ptr() as *mut __m256i, pixels);
+    }
+
+    let reminder = chunks.into_remainder();
+    if !reminder.is_empty() {
+        let mut src_pixels = [U16x2::new([0, 0]); 8];
+        src_pixels
+            .iter_mut()
+            .zip(reminder.iter())
+            .for_each(|(d, s)| *d = *s);
+
+        let mut dst_pixels = [U16x2::new([0, 0]); 8];
+        let mut pixels = _mm256_loadu_si256(src_pixels.as_ptr() as *const __m256i);
+        pixels = divide_alpha_8_pixels(pixels);
+        _mm256_storeu_si256(dst_pixels.as_mut_ptr() as *mut __m256i, pixels);
+
+        dst_pixels.iter().zip(reminder).for_each(|(s, d)| *d = *s);
+    }
+}
+
 #[inline]
 #[target_feature(enable = "avx2")]
-unsafe fn divide_alpha_eight_pixels(src: *const U16x2, dst: *mut U16x2) {
+unsafe fn divide_alpha_8_pixels(pixels: __m256i) -> __m256i {
     let alpha_mask = _mm256_set1_epi32(0xffff0000u32 as i32);
     let luma_mask = _mm256_set1_epi32(0xffff);
     let alpha_max = _mm256_set1_ps(65535.0);
@@ -144,15 +215,13 @@ unsafe fn divide_alpha_eight_pixels(src: *const U16x2, dst: *mut U16x2) {
         -1, -1, 15, 14, -1, -1, 11, 10, -1, -1, 7, 6, -1, -1, 3, 2,
     );
 
-    let src_pixels = _mm256_loadu_si256(src as *const __m256i);
-    let alpha_f32x8 = _mm256_cvtepi32_ps(_mm256_shuffle_epi8(src_pixels, alpha32_sh));
-    let luma_i32x8 = _mm256_and_si256(src_pixels, luma_mask);
+    let alpha_f32x8 = _mm256_cvtepi32_ps(_mm256_shuffle_epi8(pixels, alpha32_sh));
+    let luma_i32x8 = _mm256_and_si256(pixels, luma_mask);
     let luma_f32x8 = _mm256_cvtepi32_ps(luma_i32x8);
     let scaled_luma_f32x8 = _mm256_mul_ps(luma_f32x8, alpha_max);
     let divided_luma_f32x8 = _mm256_div_ps(scaled_luma_f32x8, alpha_f32x8);
     let divided_luma_i32x8 = _mm256_cvtps_epi32(divided_luma_f32x8);
 
-    let alpha = _mm256_and_si256(src_pixels, alpha_mask);
-    let dst_pixels = _mm256_blendv_epi8(divided_luma_i32x8, alpha, alpha_mask);
-    _mm256_storeu_si256(dst as *mut __m256i, dst_pixels);
+    let alpha = _mm256_and_si256(pixels, alpha_mask);
+    _mm256_blendv_epi8(divided_luma_i32x8, alpha, alpha_mask)
 }

src/alpha/u16x2/mod.rs

@@ -7,6 +7,8 @@ use super::AlphaMulDiv;
 #[cfg(target_arch = "x86_64")]
 mod avx2;
 mod native;
+#[cfg(target_arch = "aarch64")]
+mod neon;
 #[cfg(target_arch = "x86_64")]
 mod sse4;
 

src/alpha/u16x2/native.rs

@@ -12,9 +12,8 @@ pub(crate) fn multiply_alpha(src_image: &ImageView<U16x2>, dst_image: &mut Image
 }
 
 pub(crate) fn multiply_alpha_inplace(image: &mut ImageViewMut<U16x2>) {
-    for dst_row in image.iter_rows_mut() {
-        let src_row = unsafe { std::slice::from_raw_parts(dst_row.as_ptr(), dst_row.len()) };
-        multiply_alpha_row(src_row, dst_row);
+    for row in image.iter_rows_mut() {
+        multiply_alpha_row_inplace(row);
     }
 }
 
@@ -27,6 +26,15 @@ pub(crate) fn multiply_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2]) {
     }
 }
 
+#[inline(always)]
+pub(crate) fn multiply_alpha_row_inplace(row: &mut [U16x2]) {
+    for pixel in row {
+        let components: [u16; 2] = pixel.0;
+        let alpha = components[1];
+        pixel.0 = [mul_div_65535(components[0], alpha), alpha];
+    }
+}
+
 // Divide
 
 #[inline]
@@ -41,9 +49,8 @@ pub(crate) fn divide_alpha(src_image: &ImageView<U16x2>, dst_image: &mut ImageVi
 
 #[inline]
 pub(crate) fn divide_alpha_inplace(image: &mut ImageViewMut<U16x2>) {
-    for dst_row in image.iter_rows_mut() {
-        let src_row = unsafe { std::slice::from_raw_parts(dst_row.as_ptr(), dst_row.len()) };
-        divide_alpha_row(src_row, dst_row);
+    for row in image.iter_rows_mut() {
+        divide_alpha_row_inplace(row);
     }
 }
 
@@ -59,3 +66,13 @@ pub(crate) fn divide_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2]) {
             dst_pixel.0 = [div_and_clip16(components[0], recip_alpha), alpha];
         });
 }
+
+#[inline(always)]
+pub(crate) fn divide_alpha_row_inplace(row: &mut [U16x2]) {
+    for pixel in row {
+        let components: [u16; 2] = pixel.0;
+        let alpha = components[1];
+        let recip_alpha = RECIP_ALPHA16[alpha as usize];
+        pixel.0 = [div_and_clip16(components[0], recip_alpha), alpha];
+    }
+}

src/alpha/u16x2/neon.rs

@@ -0,0 +1,210 @@
+use std::arch::aarch64::*;
+
+use crate::neon_utils;
+use crate::pixels::U16x2;
+use crate::{ImageView, ImageViewMut};
+
+use super::native;
+
+#[target_feature(enable = "neon")]
+pub(crate) unsafe fn multiply_alpha(
+    src_image: &ImageView<U16x2>,
+    dst_image: &mut ImageViewMut<U16x2>,
+) {
+    let src_rows = src_image.iter_rows(0);
+    let dst_rows = dst_image.iter_rows_mut();
+
+    for (src_row, dst_row) in src_rows.zip(dst_rows) {
+        multiply_alpha_row(src_row, dst_row);
+    }
+}
+
+#[target_feature(enable = "neon")]
+pub(crate) unsafe fn multiply_alpha_inplace(image: &mut ImageViewMut<U16x2>) {
+    for row in image.iter_rows_mut() {
+        multiply_alpha_row_inplace(row);
+    }
+}
+
+#[inline(always)]
+unsafe fn multiply_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2]) {
+    let src_chunks = src_row.chunks_exact(8);
+    let src_remainder = src_chunks.remainder();
+    let mut dst_chunks = dst_row.chunks_exact_mut(8);
+    // Using a simple for-loop in this case is faster than implementation with pre-reading
+    for (src, dst) in src_chunks.zip(&mut dst_chunks) {
+        let mut pixels = neon_utils::load_deintrel_u16x8x2(src, 0);
+        pixels.0 = neon_utils::multiply_color_to_alpha_u16x8(pixels.0, pixels.1);
+        let dst_ptr = dst.as_mut_ptr() as *mut u16;
+        vst2q_u16(dst_ptr, pixels);
+    }
+
+    if !src_remainder.is_empty() {
+        let src_chunks = src_remainder.chunks_exact(4);
+        let src_remainder = src_chunks.remainder();
+        let dst_reminder = dst_chunks.into_remainder();
+        let mut dst_chunks = dst_reminder.chunks_exact_mut(4);
+        let mut src_dst = src_chunks.zip(&mut dst_chunks);
+        if let Some((src, dst)) = src_dst.next() {
+            let mut pixels = neon_utils::load_deintrel_u16x4x2(src, 0);
+            pixels.0 = neon_utils::multiply_color_to_alpha_u16x4(pixels.0, pixels.1);
+            let dst_ptr = dst.as_mut_ptr() as *mut u16;
+            vst2_u16(dst_ptr, pixels);
+        }
+
+        if !src_remainder.is_empty() {
+            let dst_reminder = dst_chunks.into_remainder();
+            native::multiply_alpha_row(src_remainder, dst_reminder);
+        }
+    }
+}
+
+#[inline(always)]
+unsafe fn multiply_alpha_row_inplace(row: &mut [U16x2]) {
+    let mut chunks = row.chunks_exact_mut(8);
+    // Using a simple for-loop in this case is faster than implementation with pre-reading
+    for chunk in &mut chunks {
+        let mut pixels = neon_utils::load_deintrel_u16x8x2(chunk, 0);
+        pixels.0 = neon_utils::multiply_color_to_alpha_u16x8(pixels.0, pixels.1);
+        let dst_ptr = chunk.as_mut_ptr() as *mut u16;
+        vst2q_u16(dst_ptr, pixels);
+    }
+
+    let reminder = chunks.into_remainder();
+    if !reminder.is_empty() {
+        let mut chunks = reminder.chunks_exact_mut(4);
+        if let Some(chunk) = chunks.next() {
+            let mut pixels = neon_utils::load_deintrel_u16x4x2(chunk, 0);
+            pixels.0 = neon_utils::multiply_color_to_alpha_u16x4(pixels.0, pixels.1);
+            let dst_ptr = chunk.as_mut_ptr() as *mut u16;
+            vst2_u16(dst_ptr, pixels);
+        }
+
+        let reminder = chunks.into_remainder();
+        if !reminder.is_empty() {
+            native::multiply_alpha_row_inplace(reminder);
+        }
+    }
+}
+
+// Divide
+
+#[target_feature(enable = "neon")]
+pub(crate) unsafe fn divide_alpha(
+    src_image: &ImageView<U16x2>,
+    dst_image: &mut ImageViewMut<U16x2>,
+) {
+    let src_rows = src_image.iter_rows(0);
+    let dst_rows = dst_image.iter_rows_mut();
+
+    for (src_row, dst_row) in src_rows.zip(dst_rows) {
+        divide_alpha_row(src_row, dst_row);
+    }
+}
+
+#[target_feature(enable = "neon")]
+pub(crate) unsafe fn divide_alpha_inplace(image: &mut ImageViewMut<U16x2>) {
+    for row in image.iter_rows_mut() {
+        divide_alpha_row_inplace(row);
+    }
+}
+
+#[inline(always)]
+pub(crate) unsafe fn divide_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2]) {
+    let src_chunks = src_row.chunks_exact(8);
+    let src_remainder = src_chunks.remainder();
+    let mut dst_chunks = dst_row.chunks_exact_mut(8);
+    // Using a simple for-loop in this case is faster than implementation with pre-reading
+    for (src, dst) in src_chunks.zip(&mut dst_chunks) {
+        let mut pixels = neon_utils::load_deintrel_u16x8x2(src, 0);
+        pixels = divide_alpha_8_pixels(pixels);
+        let dst_ptr = dst.as_mut_ptr() as *mut u16;
+        vst2q_u16(dst_ptr, pixels);
+    }
+
+    if !src_remainder.is_empty() {
+        let dst_reminder = dst_chunks.into_remainder();
+        let mut src_pixels = [U16x2::new([0, 0]); 8];
+        src_pixels
+            .iter_mut()
+            .zip(src_remainder)
+            .for_each(|(d, s)| *d = *s);
+
+        let mut dst_pixels = [U16x2::new([0, 0]); 8];
+        let mut pixels = neon_utils::load_deintrel_u16x8x2(&src_pixels, 0);
+        pixels = divide_alpha_8_pixels(pixels);
+        let dst_ptr = dst_pixels.as_mut_ptr() as *mut u16;
+        vst2q_u16(dst_ptr, pixels);
+
+        dst_pixels
+            .iter()
+            .zip(dst_reminder)
+            .for_each(|(s, d)| *d = *s);
+    }
+}
+
+#[inline(always)]
+pub(crate) unsafe fn divide_alpha_row_inplace(row: &mut [U16x2]) {
+    let mut chunks = row.chunks_exact_mut(8);
+    // Using a simple for-loop in this case is faster than implementation with pre-reading
+    for chunk in &mut chunks {
+        let mut pixels = neon_utils::load_deintrel_u16x8x2(chunk, 0);
+        pixels = divide_alpha_8_pixels(pixels);
+        let dst_ptr = chunk.as_mut_ptr() as *mut u16;
+        vst2q_u16(dst_ptr, pixels);
+    }
+
+    let reminder = chunks.into_remainder();
+    if !reminder.is_empty() {
+        let mut src_pixels = [U16x2::new([0, 0]); 8];
+        src_pixels
+            .iter_mut()
+            .zip(reminder.iter())
+            .for_each(|(d, s)| *d = *s);
+
+        let mut dst_pixels = [U16x2::new([0, 0]); 8];
+        let mut pixels = neon_utils::load_deintrel_u16x8x2(&src_pixels, 0);
+        pixels = divide_alpha_8_pixels(pixels);
+        let dst_ptr = dst_pixels.as_mut_ptr() as *mut u16;
+        vst2q_u16(dst_ptr, pixels);
+
+        dst_pixels.iter().zip(reminder).for_each(|(s, d)| *d = *s);
+    }
+}
+
+#[inline(always)]
+unsafe fn divide_alpha_8_pixels(mut pixels: uint16x8x2_t) -> uint16x8x2_t {
+    let zero = vdupq_n_u16(0);
+    let alpha_scale = vdupq_n_f32(65535.0);
+    let nonzero_alpha_mask = vmvnq_u16(vceqzq_u16(pixels.1));
+
+    // Low
+    let alpha_scaled_u32 = vreinterpretq_u32_u16(vzip1q_u16(pixels.1, zero));
+    let alpha_scaled_f32 = vcvtq_f32_u32(alpha_scaled_u32);
+    let recip_alpha_lo_f32 = vdivq_f32(alpha_scale, alpha_scaled_f32);
+
+    // High
+    let alpha_scaled_u32 = vreinterpretq_u32_u16(vzip2q_u16(pixels.1, zero));
+    let alpha_scaled_f32 = vcvtq_f32_u32(alpha_scaled_u32);
+    let recip_alpha_hi_f32 = vdivq_f32(alpha_scale, alpha_scaled_f32);
+
+    pixels.0 = mul_color_recip_alpha(pixels.0, recip_alpha_lo_f32, recip_alpha_hi_f32, zero);
+    pixels.0 = vandq_u16(pixels.0, nonzero_alpha_mask);
+    pixels
+}
+
+#[inline(always)]
+unsafe fn mul_color_recip_alpha(
+    color: uint16x8_t,
+    recip_alpha_lo: float32x4_t,
+    recip_alpha_hi: float32x4_t,
+    zero: uint16x8_t,
+) -> uint16x8_t {
+    let color_lo_f32 = vcvtq_f32_u32(vreinterpretq_u32_u16(vzip1q_u16(color, zero)));
+    let res_lo_u32 = vcvtaq_u32_f32(vmulq_f32(color_lo_f32, recip_alpha_lo));
+
+    let color_hi_f32 = vcvtq_f32_u32(vreinterpretq_u32_u16(vzip2q_u16(color, zero)));
+    let res_hi_u32 = vcvtaq_u32_f32(vmulq_f32(color_hi_f32, recip_alpha_hi));
+
+    vcombine_u16(vmovn_u32(res_lo_u32), vmovn_u32(res_hi_u32))
+}

src/alpha/u16x2/sse4.rs

@@ -1,6 +1,7 @@
 use std::arch::x86_64::*;
 
 use crate::pixels::U16x2;
+use crate::utils::foreach_with_pre_reading;
 use crate::{ImageView, ImageViewMut};
 
 use super::native;
@@ -20,15 +21,63 @@ pub(crate) unsafe fn multiply_alpha(
 
 #[target_feature(enable = "sse4.1")]
 pub(crate) unsafe fn multiply_alpha_inplace(image: &mut ImageViewMut<U16x2>) {
-    for dst_row in image.iter_rows_mut() {
-        let src_row = std::slice::from_raw_parts(dst_row.as_ptr(), dst_row.len());
-        multiply_alpha_row(src_row, dst_row);
+    for row in image.iter_rows_mut() {
+        multiply_alpha_row_inplace(row);
     }
 }
 
 #[inline]
 #[target_feature(enable = "sse4.1")]
 pub(crate) unsafe fn multiply_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2]) {
+    let src_chunks = src_row.chunks_exact(4);
+    let src_remainder = src_chunks.remainder();
+    let mut dst_chunks = dst_row.chunks_exact_mut(4);
+    let src_dst = src_chunks.zip(&mut dst_chunks);
+    foreach_with_pre_reading(
+        src_dst,
+        |(src, dst)| {
+            let pixels = _mm_loadu_si128(src.as_ptr() as *const __m128i);
+            let dst_ptr = dst.as_mut_ptr() as *mut __m128i;
+            (pixels, dst_ptr)
+        },
+        |(mut pixels, dst_ptr)| {
+            pixels = multiplies_alpha_4_pixels(pixels);
+            _mm_storeu_si128(dst_ptr, pixels);
+        },
+    );
+
+    if !src_remainder.is_empty() {
+        let dst_reminder = dst_chunks.into_remainder();
+        native::multiply_alpha_row(src_remainder, dst_reminder);
+    }
+}
+
+#[inline]
+#[target_feature(enable = "sse4.1")]
+pub(crate) unsafe fn multiply_alpha_row_inplace(row: &mut [U16x2]) {
+    let mut chunks = row.chunks_exact_mut(4);
+    foreach_with_pre_reading(
+        &mut chunks,
+        |chunk| {
+            let pixels = _mm_loadu_si128(chunk.as_ptr() as *const __m128i);
+            let dst_ptr = chunk.as_mut_ptr() as *mut __m128i;
+            (pixels, dst_ptr)
+        },
+        |(mut pixels, dst_ptr)| {
+            pixels = multiplies_alpha_4_pixels(pixels);
+            _mm_storeu_si128(dst_ptr, pixels);
+        },
+    );
+
+    let reminder = chunks.into_remainder();
+    if !reminder.is_empty() {
+        native::multiply_alpha_row_inplace(reminder);
+    }
+}
+
+#[inline]
+#[target_feature(enable = "sse4.1")]
+unsafe fn multiplies_alpha_4_pixels(pixels: __m128i) -> __m128i {
     let zero = _mm_setzero_si128();
     let half = _mm_set1_epi32(0x8000);
 
@@ -40,37 +89,22 @@ pub(crate) unsafe fn multiply_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2]
     */
     let factor_mask = _mm_set_epi8(15, 14, 15, 14, 11, 10, 11, 10, 7, 6, 7, 6, 3, 2, 3, 2);
 
-    let src_chunks = src_row.chunks_exact(4);
-    let src_remainder = src_chunks.remainder();
-    let mut dst_chunks = dst_row.chunks_exact_mut(4);
+    let factor_pixels = _mm_shuffle_epi8(pixels, factor_mask);
+    let factor_pixels = _mm_or_si128(factor_pixels, max_alpha);
 
-    for (src, dst) in src_chunks.zip(&mut dst_chunks) {
-        let src_pixels = _mm_loadu_si128(src.as_ptr() as *const __m128i);
+    let src_i32_lo = _mm_unpacklo_epi16(pixels, zero);
+    let factors = _mm_unpacklo_epi16(factor_pixels, zero);
+    let src_i32_lo = _mm_add_epi32(_mm_mullo_epi32(src_i32_lo, factors), half);
+    let dst_i32_lo = _mm_add_epi32(src_i32_lo, _mm_srli_epi32::<16>(src_i32_lo));
+    let dst_i32_lo = _mm_srli_epi32::<16>(dst_i32_lo);
 
-        let factor_pixels = _mm_shuffle_epi8(src_pixels, factor_mask);
-        let factor_pixels = _mm_or_si128(factor_pixels, max_alpha);
+    let src_i32_hi = _mm_unpackhi_epi16(pixels, zero);
+    let factors = _mm_unpackhi_epi16(factor_pixels, zero);
+    let src_i32_hi = _mm_add_epi32(_mm_mullo_epi32(src_i32_hi, factors), half);
+    let dst_i32_hi = _mm_add_epi32(src_i32_hi, _mm_srli_epi32::<16>(src_i32_hi));
+    let dst_i32_hi = _mm_srli_epi32::<16>(dst_i32_hi);
 
-        let src_i32_lo = _mm_unpacklo_epi16(src_pixels, zero);
-        let factors = _mm_unpacklo_epi16(factor_pixels, zero);
-        let src_i32_lo = _mm_add_epi32(_mm_mullo_epi32(src_i32_lo, factors), half);
-        let dst_i32_lo = _mm_add_epi32(src_i32_lo, _mm_srli_epi32::<16>(src_i32_lo));
-        let dst_i32_lo = _mm_srli_epi32::<16>(dst_i32_lo);
-
-        let src_i32_hi = _mm_unpackhi_epi16(src_pixels, zero);
-        let factors = _mm_unpackhi_epi16(factor_pixels, zero);
-        let src_i32_hi = _mm_add_epi32(_mm_mullo_epi32(src_i32_hi, factors), half);
-        let dst_i32_hi = _mm_add_epi32(src_i32_hi, _mm_srli_epi32::<16>(src_i32_hi));
-        let dst_i32_hi = _mm_srli_epi32::<16>(dst_i32_hi);
-
-        let dst_pixels = _mm_packus_epi32(dst_i32_lo, dst_i32_hi);
-
-        _mm_storeu_si128(dst.as_mut_ptr() as *mut __m128i, dst_pixels);
-    }
-
-    if !src_remainder.is_empty() {
-        let dst_reminder = dst_chunks.into_remainder();
-        native::multiply_alpha_row(src_remainder, dst_reminder);
-    }
+    _mm_packus_epi32(dst_i32_lo, dst_i32_hi)
 }
 
 // Divide
@@ -90,9 +124,8 @@ pub(crate) unsafe fn divide_alpha(
 
 #[target_feature(enable = "sse4.1")]
 pub(crate) unsafe fn divide_alpha_inplace(image: &mut ImageViewMut<U16x2>) {
-    for dst_row in image.iter_rows_mut() {
-        let src_row = std::slice::from_raw_parts(dst_row.as_ptr(), dst_row.len());
-        divide_alpha_row(src_row, dst_row);
+    for row in image.iter_rows_mut() {
+        divide_alpha_row_inplace(row);
     }
 }
 
@@ -101,10 +134,19 @@ pub(crate) unsafe fn divide_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2])
     let src_chunks = src_row.chunks_exact(4);
     let src_remainder = src_chunks.remainder();
     let mut dst_chunks = dst_row.chunks_exact_mut(4);
-
-    for (src, dst) in src_chunks.zip(&mut dst_chunks) {
-        divide_alpha_four_pixels(src.as_ptr(), dst.as_mut_ptr());
-    }
+    let src_dst = src_chunks.zip(&mut dst_chunks);
+    foreach_with_pre_reading(
+        src_dst,
+        |(src, dst)| {
+            let pixels = _mm_loadu_si128(src.as_ptr() as *const __m128i);
+            let dst_ptr = dst.as_mut_ptr() as *mut __m128i;
+            (pixels, dst_ptr)
+        },
+        |(mut pixels, dst_ptr)| {
+            pixels = divide_alpha_4_pixels(pixels);
+            _mm_storeu_si128(dst_ptr, pixels);
+        },
+    );
 
     if !src_remainder.is_empty() {
         let dst_reminder = dst_chunks.into_remainder();
@@ -115,7 +157,9 @@ pub(crate) unsafe fn divide_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2])
             .for_each(|(d, s)| *d = *s);
 
         let mut dst_pixels = [U16x2::new([0, 0]); 4];
-        divide_alpha_four_pixels(src_pixels.as_ptr(), dst_pixels.as_mut_ptr());
+        let mut pixels = _mm_loadu_si128(src_pixels.as_ptr() as *const __m128i);
+        pixels = divide_alpha_4_pixels(pixels);
+        _mm_storeu_si128(dst_pixels.as_mut_ptr() as *mut __m128i, pixels);
 
         dst_pixels
             .iter()
@@ -124,9 +168,36 @@ pub(crate) unsafe fn divide_alpha_row(src_row: &[U16x2], dst_row: &mut [U16x2])
     }
 }
 
+#[target_feature(enable = "sse4.1")]
+pub(crate) unsafe fn divide_alpha_row_inplace(row: &mut [U16x2]) {
+    let mut chunks = row.chunks_exact_mut(4);
+    // Using a simple for-loop in this case is faster than implementation with pre-reading
+    for chunk in &mut chunks {
+        let mut pixels = _mm_loadu_si128(chunk.as_ptr() as *const __m128i);
+        pixels = divide_alpha_4_pixels(pixels);
+        _mm_storeu_si128(chunk.as_mut_ptr() as *mut __m128i, pixels);
+    }
+
+    let reminder = chunks.into_remainder();
+    if !reminder.is_empty() {
+        let mut src_pixels = [U16x2::new([0, 0]); 4];
+        src_pixels
+            .iter_mut()
+            .zip(reminder.iter())
+            .for_each(|(d, s)| *d = *s);
+
+        let mut dst_pixels = [U16x2::new([0, 0]); 4];
+        let mut pixels = _mm_loadu_si128(src_pixels.as_ptr() as *const __m128i);
+        pixels = divide_alpha_4_pixels(pixels);
+        _mm_storeu_si128(dst_pixels.as_mut_ptr() as *mut __m128i, pixels);
+
+        dst_pixels.iter().zip(reminder).for_each(|(s, d)| *d = *s);
+    }
+}
+
 #[inline]
 #[target_feature(enable = "sse4.1")]
-unsafe fn divide_alpha_four_pixels(src: *const U16x2, dst: *mut U16x2) {
+unsafe fn divide_alpha_4_pixels(pixels: __m128i) -> __m128i {
     let alpha_mask = _mm_set1_epi32(0xffff0000u32 as i32);
     let luma_mask = _mm_set1_epi32(0xffff);
     let alpha_max = _mm_set1_ps(65535.0);
@@ -136,13 +207,11 @@ unsafe fn divide_alpha_four_pixels(src: *const U16x2, dst: *mut U16x2) {
     */
     let alpha32_sh = _mm_set_epi8(-1, -1, 15, 14, -1, -1, 11, 10, -1, -1, 7, 6, -1, -1, 3, 2);
 
-    let src_pixels = _mm_loadu_si128(src as *const __m128i);
-    let alpha_f32x4 = _mm_cvtepi32_ps(_mm_shuffle_epi8(src_pixels, alpha32_sh));
-    let luma_f32x4 = _mm_cvtepi32_ps(_mm_and_si128(src_pixels, luma_mask));
+    let alpha_f32x4 = _mm_cvtepi32_ps(_mm_shuffle_epi8(pixels, alpha32_sh));
+    let luma_f32x4 = _mm_cvtepi32_ps(_mm_and_si128(pixels, luma_mask));
     let scaled_luma_f32x4 = _mm_mul_ps(luma_f32x4, alpha_max);
     let divided_luma_i32x4 = _mm_cvtps_epi32(_mm_div_ps(scaled_luma_f32x4, alpha_f32x4));
 
-    let alpha = _mm_and_si128(src_pixels, alpha_mask);
-    let dst_pixels = _mm_blendv_epi8(divided_luma_i32x4, alpha, alpha_mask);
-    _mm_storeu_si128(dst as *mut __m128i, dst_pixels);
+    let alpha = _mm_and_si128(pixels, alpha_mask);
+    _mm_blendv_epi8(divided_luma_i32x4, alpha, alpha_mask)
 }

src/alpha/u16x4/neon.rs

@@ -41,9 +41,9 @@ unsafe fn multiply_alpha_row(src_row: &[U16x4], dst_row: &mut [U16x4]) {
             (pixels, dst_ptr)
         },
         |(mut pixels, dst_ptr)| {
-            pixels.0 = multiply_color_to_alpha_u16x8(pixels.0, pixels.3);
-            pixels.1 = multiply_color_to_alpha_u16x8(pixels.1, pixels.3);
-            pixels.2 = multiply_color_to_alpha_u16x8(pixels.2, pixels.3);
+            pixels.0 = neon_utils::multiply_color_to_alpha_u16x8(pixels.0, pixels.3);
+            pixels.1 = neon_utils::multiply_color_to_alpha_u16x8(pixels.1, pixels.3);
+            pixels.2 = neon_utils::multiply_color_to_alpha_u16x8(pixels.2, pixels.3);
             vst4q_u16(dst_ptr, pixels);
         },
     );
@@ -55,9 +55,9 @@ unsafe fn multiply_alpha_row(src_row: &[U16x4], dst_row: &mut [U16x4]) {
     let mut src_dst = src_chunks.zip(&mut dst_chunks);
     if let Some((src, dst)) = src_dst.next() {
         let mut pixels = neon_utils::load_deintrel_u16x4x4(src, 0);
-        pixels.0 = multiply_color_to_alpha_u16x4(pixels.0, pixels.3);
-        pixels.1 = multiply_color_to_alpha_u16x4(pixels.1, pixels.3);
-        pixels.2 = multiply_color_to_alpha_u16x4(pixels.2, pixels.3);
+        pixels.0 = neon_utils::multiply_color_to_alpha_u16x4(pixels.0, pixels.3);
+        pixels.1 = neon_utils::multiply_color_to_alpha_u16x4(pixels.1, pixels.3);
+        pixels.2 = neon_utils::multiply_color_to_alpha_u16x4(pixels.2, pixels.3);
         let dst_ptr = dst.as_mut_ptr() as *mut u16;
         vst4_u16(dst_ptr, pixels);
     }
@@ -79,9 +79,9 @@ unsafe fn multiply_alpha_row_inplace(row: &mut [U16x4]) {
             (pixels, dst_ptr)
         },
         |(mut pixels, dst_ptr)| {
-            pixels.0 = multiply_color_to_alpha_u16x8(pixels.0, pixels.3);
-            pixels.1 = multiply_color_to_alpha_u16x8(pixels.1, pixels.3);
-            pixels.2 = multiply_color_to_alpha_u16x8(pixels.2, pixels.3);
+            pixels.0 = neon_utils::multiply_color_to_alpha_u16x8(pixels.0, pixels.3);
+            pixels.1 = neon_utils::multiply_color_to_alpha_u16x8(pixels.1, pixels.3);
+            pixels.2 = neon_utils::multiply_color_to_alpha_u16x8(pixels.2, pixels.3);
             vst4q_u16(dst_ptr, pixels);
         },
     );
@@ -90,9 +90,9 @@ unsafe fn multiply_alpha_row_inplace(row: &mut [U16x4]) {
     let mut chunks = reminder.chunks_exact_mut(4);
     if let Some(chunk) = chunks.next() {
         let mut pixels = neon_utils::load_deintrel_u16x4x4(chunk, 0);
-        pixels.0 = multiply_color_to_alpha_u16x4(pixels.0, pixels.3);
-        pixels.1 = multiply_color_to_alpha_u16x4(pixels.1, pixels.3);
-        pixels.2 = multiply_color_to_alpha_u16x4(pixels.2, pixels.3);
+        pixels.0 = neon_utils::multiply_color_to_alpha_u16x4(pixels.0, pixels.3);
+        pixels.1 = neon_utils::multiply_color_to_alpha_u16x4(pixels.1, pixels.3);
+        pixels.2 = neon_utils::multiply_color_to_alpha_u16x4(pixels.2, pixels.3);
         let dst_ptr = chunk.as_mut_ptr() as *mut u16;
         vst4_u16(dst_ptr, pixels);
     }
@@ -103,23 +103,6 @@ unsafe fn multiply_alpha_row_inplace(row: &mut [U16x4]) {
     }
 }
 
-#[inline(always)]
-unsafe fn multiply_color_to_alpha_u16x8(color: uint16x8_t, alpha: uint16x8_t) -> uint16x8_t {
-    let rounder = vdupq_n_u32(0x8000);
-    let color_lo_u32 = vmlal_u16(rounder, vget_low_u16(color), vget_low_u16(alpha));
-    let color_hi_u32 = vmlal_high_u16(rounder, color, alpha);
-    let color_lo_u32 = vaddhn_u32(color_lo_u32, vshrq_n_u32::<16>(color_lo_u32));
-    let color_hi_u32 = vaddhn_u32(color_hi_u32, vshrq_n_u32::<16>(color_hi_u32));
-    vcombine_u16(color_lo_u32, color_hi_u32)
-}
-
-#[inline(always)]
-unsafe fn multiply_color_to_alpha_u16x4(color: uint16x4_t, alpha: uint16x4_t) -> uint16x4_t {
-    let rounder = vdupq_n_u32(0x8000);
-    let color_u32 = vmlal_u16(rounder, color, alpha);
-    vaddhn_u32(color_u32, vshrq_n_u32::<16>(color_u32))
-}
-
 // Divide
 
 #[target_feature(enable = "neon")]

src/neon_utils.rs

@@ -128,6 +128,16 @@ pub unsafe fn load_deintrel_u16x4x4<T>(buf: &[T], index: usize) -> uint16x4x4_t
     vld4_u16(buf.get_unchecked(index..).as_ptr() as *const u16)
 }
 
+#[inline(always)]
+pub unsafe fn load_deintrel_u16x4x2<T>(buf: &[T], index: usize) -> uint16x4x2_t {
+    vld2_u16(buf.get_unchecked(index..).as_ptr() as *const u16)
+}
+
+#[inline(always)]
+pub unsafe fn load_deintrel_u16x8x2<T>(buf: &[T], index: usize) -> uint16x8x2_t {
+    vld2q_u16(buf.get_unchecked(index..).as_ptr() as *const u16)
+}
+
 #[inline(always)]
 pub unsafe fn load_deintrel_u16x8x3<T>(buf: &[T], index: usize) -> uint16x8x3_t {
     vld3q_u16(buf.get_unchecked(index..).as_ptr() as *const u16)
@@ -305,6 +315,23 @@ pub unsafe fn mul_color_to_alpha_u8x8(
     vqmovn_u16(res_u16)
 }
 
+#[inline(always)]
+pub unsafe fn multiply_color_to_alpha_u16x8(color: uint16x8_t, alpha: uint16x8_t) -> uint16x8_t {
+    let rounder = vdupq_n_u32(0x8000);
+    let color_lo_u32 = vmlal_u16(rounder, vget_low_u16(color), vget_low_u16(alpha));
+    let color_hi_u32 = vmlal_high_u16(rounder, color, alpha);
+    let color_lo_u16 = vaddhn_u32(color_lo_u32, vshrq_n_u32::<16>(color_lo_u32));
+    let color_hi_u16 = vaddhn_u32(color_hi_u32, vshrq_n_u32::<16>(color_hi_u32));
+    vcombine_u16(color_lo_u16, color_hi_u16)
+}
+
+#[inline(always)]
+pub unsafe fn multiply_color_to_alpha_u16x4(color: uint16x4_t, alpha: uint16x4_t) -> uint16x4_t {
+    let rounder = vdupq_n_u32(0x8000);
+    let color_u32 = vmlal_u16(rounder, color, alpha);
+    vaddhn_u32(color_u32, vshrq_n_u32::<16>(color_u32))
+}
+
 #[inline]
 #[target_feature(enable = "neon")]
 pub unsafe fn mul_color_recip_alpha_u8x16(
@@ -330,7 +357,6 @@ pub unsafe fn mul_color_recip_alpha_u8x8(
     let color_u16_lo = vreinterpret_u16_u8(vzip1_u8(zero, color));
     let color_u16_hi = vreinterpret_u16_u8(vzip2_u8(zero, color));
     let color_u16 = vcombine_u16(color_u16_lo, color_u16_hi);
-
     let res_u16 = mulhi_u16x8(color_u16, recip_alpha);
     vmovn_u16(res_u16)
 }

tests/alpha_tests.rs

@@ -252,6 +252,12 @@ mod multiply_alpha_u16x2 {
         mul_div_alpha_test(Oper::Mul, SRC_PIXELS, RES_PIXELS, CpuExtensions::Sse4_1);
     }
 
+    #[cfg(target_arch = "aarch64")]
+    #[test]
+    fn neon_test() {
+        mul_div_alpha_test(Oper::Mul, SRC_PIXELS, RES_PIXELS, CpuExtensions::Neon);
+    }
+
     #[test]
     fn native_test() {
         mul_div_alpha_test(Oper::Mul, SRC_PIXELS, RES_PIXELS, CpuExtensions::None);
@@ -439,6 +445,12 @@ mod divide_alpha_u16x2 {
         mul_div_alpha_test(OPER, SRC_PIXELS, SIMD_RES_PIXELS, CpuExtensions::Sse4_1);
     }
 
+    #[cfg(target_arch = "aarch64")]
+    #[test]
+    fn neon_test() {
+        mul_div_alpha_test(OPER, SRC_PIXELS, SIMD_RES_PIXELS, CpuExtensions::Neon);
+    }
+
     #[test]
     fn native_test() {
         mul_div_alpha_test(OPER, SRC_PIXELS, RES_PIXELS, CpuExtensions::None);