Added support of new type of pixels PixelType::U16x4.

2026-10-08 01:11:09 +00:00 · 2022-06-26 19:16:22 +03:00
parent 091f78fac9
commit 891312a4bb
26 changed files with 1575 additions and 23 deletions
@@ -1,3 +1,7 @@
+## [Unreleased] - ReleaseDate
+
+- Added support of new type of pixels `PixelType::U16x4`.
+
 ## [0.9.5] - 2022-06-22

 - Fixed README.md
@@ -59,6 +59,11 @@ name = "bench_compare_rgba"
 harness = false


+[[bench]]
+name = "bench_compare_rgba16"
+harness = false
+
+
 [[bench]]
 name = "bench_compare_l"
 harness = false
@@ -30,6 +30,7 @@ fn multiplies_alpha(bench: &mut Bench, pixel_type: PixelType, cpu_extensions: Cp
        PixelType::U8x4 => &[255, 128, 0, 128],
        PixelType::U8x2 => &[255, 128],
        PixelType::U16x2 => &[255, 255, 0, 128],
+        PixelType::U16x4 => &[0, 255, 0, 128, 0, 0, 0, 128],
        _ => unreachable!(),
    };
    let src_data = get_src_image(width, height, pixel_type, pixel);
@@ -60,6 +61,7 @@ fn divides_alpha(bench: &mut Bench, pixel_type: PixelType, cpu_extensions: CpuEx
        PixelType::U8x4 => &[128, 64, 0, 128],
        PixelType::U8x2 => &[128, 128],
        PixelType::U16x2 => &[0, 128, 0, 128],
+        PixelType::U16x4 => &[0, 128, 0, 64, 0, 0, 0, 128],
        _ => unreachable!(),
    };
    let src_data = get_src_image(width, height, pixel_type, pixel);
@@ -84,7 +86,12 @@ fn divides_alpha(bench: &mut Bench, pixel_type: PixelType, cpu_extensions: CpuEx
 }

 fn bench_alpha(bench: &mut Bench) {
-    let pixel_types = [PixelType::U8x4, PixelType::U8x2, PixelType::U16x2];
+    let pixel_types = [
+        PixelType::U8x4,
+        PixelType::U8x2,
+        PixelType::U16x2,
+        PixelType::U16x4,
+    ];
    let mut cpu_extensions = vec![CpuExtensions::None];
    #[cfg(target_arch = "x86_64")]
    {
@@ -0,0 +1,131 @@
+use std::num::NonZeroU32;
+
+use glassbench::*;
+use image::imageops;
+use resize::px::RGBA;
+use resize::Pixel::RGBA16;
+use rgb::FromSlice;
+
+use fast_image_resize::pixels::U16x4;
+use fast_image_resize::{CpuExtensions, FilterType, Image, MulDiv, ResizeAlg, Resizer};
+use testing::PixelExt;
+
+mod utils;
+
+pub fn bench_downscale_rgba16(bench: &mut Bench) {
+    let src_image = U16x4::load_big_image().to_rgba16();
+    let new_width = NonZeroU32::new(852).unwrap();
+    let new_height = NonZeroU32::new(567).unwrap();
+
+    let alg_names = ["Nearest", "Bilinear", "CatmullRom", "Lanczos3"];
+
+    // image crate
+    // https://crates.io/crates/image
+    for alg_name in alg_names {
+        let filter = match alg_name {
+            "Nearest" => imageops::Nearest,
+            "Bilinear" => imageops::Triangle,
+            "CatmullRom" => imageops::CatmullRom,
+            "Lanczos3" => imageops::Lanczos3,
+            _ => continue,
+        };
+        bench.task(format!("image - {}", alg_name), |task| {
+            task.iter(|| {
+                imageops::resize(&src_image, new_width.get(), new_height.get(), filter);
+            })
+        });
+    }
+
+    // resize crate
+    // https://crates.io/crates/resize
+    for alg_name in alg_names {
+        let resize_src_image = src_image.as_raw().as_rgba();
+        let mut dst =
+            vec![RGBA::new(0u16, 0u16, 0u16, 0u16); (new_width.get() * new_height.get()) as usize];
+        bench.task(format!("resize - {}", alg_name), |task| {
+            let filter = match alg_name {
+                "Nearest" => {
+                    // resizer doesn't support "nearest" algorithm
+                    task.iter(|| {});
+                    return;
+                }
+                "Bilinear" => resize::Type::Triangle,
+                "CatmullRom" => resize::Type::Catrom,
+                "Lanczos3" => resize::Type::Lanczos3,
+                _ => return,
+            };
+            let mut resize = resize::new(
+                src_image.width() as usize,
+                src_image.height() as usize,
+                new_width.get() as usize,
+                new_height.get() as usize,
+                RGBA16,
+                filter,
+            )
+            .unwrap();
+            task.iter(|| {
+                resize.resize(resize_src_image, &mut dst).unwrap();
+            })
+        });
+    }
+
+    // fast_image_resize crate;
+    let src_image_data = U16x4::load_big_src_image();
+    let mut cpu_ext_and_name = vec![(CpuExtensions::None, "rust")];
+    #[cfg(target_arch = "x86_64")]
+    {
+        cpu_ext_and_name.push((CpuExtensions::Sse4_1, "sse4.1"));
+        cpu_ext_and_name.push((CpuExtensions::Avx2, "avx2"));
+    }
+    for (cpu_ext, ext_name) in cpu_ext_and_name {
+        for alg_name in alg_names {
+            let resize_alg = match alg_name {
+                "Nearest" => ResizeAlg::Nearest,
+                "Bilinear" => ResizeAlg::Convolution(FilterType::Bilinear),
+                "CatmullRom" => ResizeAlg::Convolution(FilterType::CatmullRom),
+                "Lanczos3" => ResizeAlg::Convolution(FilterType::Lanczos3),
+                _ => return,
+            };
+            let src_view = src_image_data.view();
+            let mut premultiplied_src_image = Image::new(
+                NonZeroU32::new(src_image.width()).unwrap(),
+                NonZeroU32::new(src_image.height()).unwrap(),
+                src_view.pixel_type(),
+            );
+            let mut dst_image = Image::new(new_width, new_height, src_view.pixel_type());
+            let mut dst_view = dst_image.view_mut();
+            let mut mul_div = MulDiv::default();
+
+            let mut fast_resizer = Resizer::new(resize_alg);
+
+            unsafe {
+                fast_resizer.reset_internal_buffers();
+                fast_resizer.set_cpu_extensions(cpu_ext);
+                mul_div.set_cpu_extensions(cpu_ext);
+            }
+
+            bench.task(format!("fir {} - {}", ext_name, alg_name), |task| {
+                task.iter(|| match resize_alg {
+                    ResizeAlg::Nearest => {
+                        fast_resizer
+                            .resize(&premultiplied_src_image.view(), &mut dst_view)
+                            .unwrap();
+                    }
+                    _ => {
+                        mul_div
+                            .multiply_alpha(&src_view, &mut premultiplied_src_image.view_mut())
+                            .unwrap();
+                        fast_resizer
+                            .resize(&premultiplied_src_image.view(), &mut dst_view)
+                            .unwrap();
+                        mul_div.divide_alpha_inplace(&mut dst_view).unwrap();
+                    }
+                })
+            });
+        }
+    }
+
+    utils::print_md_table(bench);
+}
+
+bench_main!("Compare resize of RGBA16 image", bench_downscale_rgba16,);
@@ -102,6 +102,7 @@ pub fn main() {
            PixelType::U16,
            PixelType::U16x2,
            PixelType::U16x3,
+            PixelType::U16x4,
            PixelType::I32,
        ];
        let mut cpu_extensions = vec![CpuExtensions::None];
@@ -120,6 +121,7 @@ pub fn main() {
                    PixelType::U16 => U16::load_big_src_image(),
                    PixelType::U16x2 => U16x2::load_big_src_image(),
                    PixelType::U16x3 => U16x3::load_big_src_image(),
+                    PixelType::U16x4 => U16x4::load_big_src_image(),
                    PixelType::I32 => I32::load_big_src_image(),
                    _ => unreachable!(),
                };
@@ -7,6 +7,7 @@ use crate::CpuExtensions;
 mod common;
 pub(crate) mod errors;
 mod u16x2;
+mod u16x4;
 mod u8x2;
 mod u8x4;

@@ -138,11 +138,9 @@ unsafe fn divide_alpha_four_pixels(src: *const U16x2, dst: *mut U16x2) {

    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_i32x4 = _mm_and_si128(src_pixels, luma_mask);
-    let luma_f32x4 = _mm_cvtepi32_ps(luma_i32x4);
+    let luma_f32x4 = _mm_cvtepi32_ps(_mm_and_si128(src_pixels, luma_mask));
    let scaled_luma_f32x4 = _mm_mul_ps(luma_f32x4, alpha_max);
-    let divided_luma_f32x4 = _mm_div_ps(scaled_luma_f32x4, alpha_f32x4);
-    let divided_luma_i32x4 = _mm_cvtps_epi32(divided_luma_f32x4);
+    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);
@@ -0,0 +1,172 @@
+use std::arch::x86_64::*;
+
+use crate::image_view::{TypedImageView, TypedImageViewMut};
+use crate::pixels::U16x4;
+
+use super::sse4;
+
+#[target_feature(enable = "avx2")]
+pub(crate) unsafe fn multiply_alpha(
+    src_image: TypedImageView<U16x4>,
+    mut dst_image: TypedImageViewMut<U16x4>,
+) {
+    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 = "avx2")]
+pub(crate) unsafe fn multiply_alpha_inplace(mut image: TypedImageViewMut<U16x4>) {
+    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);
+    }
+}
+
+#[inline]
+#[target_feature(enable = "avx2")]
+pub(crate) unsafe fn multiply_alpha_row(src_row: &[U16x4], dst_row: &mut [U16x4]) {
+    let zero = _mm256_setzero_si256();
+    let half = _mm256_set1_epi32(0x8000);
+
+    const MAX_A: i64 = 0xffff000000000000u64 as i64;
+    let max_alpha = _mm256_set1_epi64x(MAX_A);
+    /*
+       |R0   G0   B0   A0  | |R1   G1   B1   A1  | |R0   G0   B0   A0  | |R1   G1   B1   A1 |
+       |0001 0203 0405 0607| |0809 1011 1213 1415| |0001 0203 0405 0607| |0809 1011 1213 1415|
+    */
+    let factor_mask = _mm256_set_m128i(
+        _mm_set_epi8(15, 14, 15, 14, 15, 14, 15, 14, 7, 6, 7, 6, 7, 6, 7, 6),
+        _mm_set_epi8(15, 14, 15, 14, 15, 14, 15, 14, 7, 6, 7, 6, 7, 6, 7, 6),
+    );
+
+    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) {
+        let src_pixels = _mm256_loadu_si256(src.as_ptr() as *const __m256i);
+
+        let factor_pixels = _mm256_shuffle_epi8(src_pixels, factor_mask);
+        let factor_pixels = _mm256_or_si256(factor_pixels, max_alpha);
+
+        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);
+    }
+}
+
+// Divide
+
+#[target_feature(enable = "avx2")]
+pub(crate) unsafe fn divide_alpha(
+    src_image: TypedImageView<U16x4>,
+    mut dst_image: TypedImageViewMut<U16x4>,
+) {
+    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 = "avx2")]
+pub(crate) unsafe fn divide_alpha_inplace(mut image: TypedImageViewMut<U16x4>) {
+    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);
+    }
+}
+
+#[target_feature(enable = "avx2")]
+pub(crate) unsafe fn divide_alpha_row(src_row: &[U16x4], dst_row: &mut [U16x4]) {
+    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());
+    }
+
+    if !src_remainder.is_empty() {
+        let dst_reminder = dst_chunks.into_remainder();
+        let mut src_pixels = [U16x4([0, 0, 0, 0]); 4];
+        src_pixels
+            .iter_mut()
+            .zip(src_remainder)
+            .for_each(|(d, s)| *d = *s);
+
+        let mut dst_pixels = [U16x4([0, 0, 0, 0]); 4];
+        divide_alpha_four_pixels(src_pixels.as_ptr(), dst_pixels.as_mut_ptr());
+
+        dst_pixels
+            .iter()
+            .zip(dst_reminder)
+            .for_each(|(s, d)| *d = *s);
+    }
+}
+
+#[inline]
+#[target_feature(enable = "avx2")]
+unsafe fn divide_alpha_four_pixels(src: *const U16x4, dst: *mut U16x4) {
+    let zero = _mm256_setzero_si256();
+    let alpha_mask = _mm256_set1_epi64x(0xffff000000000000u64 as i64);
+    let alpha_max = _mm256_set1_ps(65535.0);
+    /*
+       |R0   G0   B0   A0  | |R1   G1   B1   A1  | |R0   G0   B0   A0  | |R1   G1   B1   A1 |
+       |0001 0203 0405 0607| |0809 1011 1213 1415| |0001 0203 0405 0607| |0809 1011 1213 1415|
+    */
+    let alpha32_sh0 = _mm256_set_m128i(
+        _mm_set_epi8(-1, -1, 7, 6, -1, -1, 7, 6, -1, -1, 7, 6, -1, -1, 7, 6),
+        _mm_set_epi8(-1, -1, 7, 6, -1, -1, 7, 6, -1, -1, 7, 6, -1, -1, 7, 6),
+    );
+    let alpha32_sh1 = _mm256_set_m128i(
+        _mm_set_epi8(
+            -1, -1, 15, 14, -1, -1, 15, 14, -1, -1, 15, 14, -1, -1, 15, 14,
+        ),
+        _mm_set_epi8(
+            -1, -1, 15, 14, -1, -1, 15, 14, -1, -1, 15, 14, -1, -1, 15, 14,
+        ),
+    );
+
+    let src_pixels = _mm256_loadu_si256(src as *const __m256i);
+
+    let alpha0_f32x8 = _mm256_cvtepi32_ps(_mm256_shuffle_epi8(src_pixels, alpha32_sh0));
+    let alpha1_f32x8 = _mm256_cvtepi32_ps(_mm256_shuffle_epi8(src_pixels, alpha32_sh1));
+
+    let pix0_f32x8 = _mm256_cvtepi32_ps(_mm256_unpacklo_epi16(src_pixels, zero));
+    let pix1_f32x8 = _mm256_cvtepi32_ps(_mm256_unpacklo_epi16(src_pixels, zero));
+
+    let scaled_pix0_f32x8 = _mm256_mul_ps(pix0_f32x8, alpha_max);
+    let scaled_pix1_f32x8 = _mm256_mul_ps(pix1_f32x8, alpha_max);
+
+    let divided_pix0_i32x8 = _mm256_cvtps_epi32(_mm256_div_ps(scaled_pix0_f32x8, alpha0_f32x8));
+    let divided_pix1_i32x8 = _mm256_cvtps_epi32(_mm256_div_ps(scaled_pix1_f32x8, alpha1_f32x8));
+
+    let two_pixels_i16x16 = _mm256_packus_epi32(divided_pix0_i32x8, divided_pix1_i32x8);
+    let alpha = _mm256_and_si256(src_pixels, alpha_mask);
+    let dst_pixels = _mm256_blendv_epi8(two_pixels_i16x16, alpha, alpha_mask);
+
+    _mm256_storeu_si256(dst as *mut __m256i, dst_pixels);
+}
@@ -0,0 +1,61 @@
+use crate::image_view::{TypedImageView, TypedImageViewMut};
+use crate::pixels::U16x4;
+use crate::CpuExtensions;
+
+use super::AlphaMulDiv;
+
+#[cfg(target_arch = "x86_64")]
+mod avx2;
+mod native;
+#[cfg(target_arch = "x86_64")]
+mod sse4;
+
+impl AlphaMulDiv for U16x4 {
+    fn multiply_alpha(
+        src_image: TypedImageView<Self>,
+        dst_image: TypedImageViewMut<Self>,
+        cpu_extensions: CpuExtensions,
+    ) {
+        match cpu_extensions {
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Avx2 => unsafe { avx2::multiply_alpha(src_image, dst_image) },
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Sse4_1 => unsafe { sse4::multiply_alpha(src_image, dst_image) },
+            _ => native::multiply_alpha(src_image, dst_image),
+        }
+    }
+
+    fn multiply_alpha_inplace(image: TypedImageViewMut<Self>, cpu_extensions: CpuExtensions) {
+        match cpu_extensions {
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Avx2 => unsafe { avx2::multiply_alpha_inplace(image) },
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Sse4_1 => unsafe { sse4::multiply_alpha_inplace(image) },
+            _ => native::multiply_alpha_inplace(image),
+        }
+    }
+
+    fn divide_alpha(
+        src_image: TypedImageView<Self>,
+        dst_image: TypedImageViewMut<Self>,
+        cpu_extensions: CpuExtensions,
+    ) {
+        match cpu_extensions {
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Avx2 => unsafe { avx2::divide_alpha(src_image, dst_image) },
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Sse4_1 => unsafe { sse4::divide_alpha(src_image, dst_image) },
+            _ => native::divide_alpha(src_image, dst_image),
+        }
+    }
+
+    fn divide_alpha_inplace(image: TypedImageViewMut<Self>, cpu_extensions: CpuExtensions) {
+        match cpu_extensions {
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Avx2 => unsafe { avx2::divide_alpha_inplace(image) },
+            #[cfg(target_arch = "x86_64")]
+            CpuExtensions::Sse4_1 => unsafe { sse4::divide_alpha_inplace(image) },
+            _ => native::divide_alpha_inplace(image),
+        }
+    }
+}
@@ -0,0 +1,77 @@
+use crate::alpha::common::{div_and_clip16, mul_div_65535, RECIP_ALPHA16};
+use crate::image_view::{TypedImageView, TypedImageViewMut};
+use crate::pixels::U16x4;
+
+pub(crate) fn multiply_alpha(
+    src_image: TypedImageView<U16x4>,
+    mut dst_image: TypedImageViewMut<U16x4>,
+) {
+    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);
+    }
+}
+
+pub(crate) fn multiply_alpha_inplace(mut image: TypedImageViewMut<U16x4>) {
+    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);
+    }
+}
+
+#[inline(always)]
+pub(crate) fn multiply_alpha_row(src_row: &[U16x4], dst_row: &mut [U16x4]) {
+    for (src_pixel, dst_pixel) in src_row.iter().zip(dst_row) {
+        let components: [u16; 4] = src_pixel.0;
+        let alpha = components[3];
+        dst_pixel.0 = [
+            mul_div_65535(components[0], alpha),
+            mul_div_65535(components[1], alpha),
+            mul_div_65535(components[2], alpha),
+            alpha,
+        ];
+    }
+}
+
+// Divide
+
+#[inline]
+pub(crate) fn divide_alpha(
+    src_image: TypedImageView<U16x4>,
+    mut dst_image: TypedImageViewMut<U16x4>,
+) {
+    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);
+    }
+}
+
+#[inline]
+pub(crate) fn divide_alpha_inplace(mut image: TypedImageViewMut<U16x4>) {
+    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);
+    }
+}
+
+#[inline(always)]
+pub(crate) fn divide_alpha_row(src_row: &[U16x4], dst_row: &mut [U16x4]) {
+    src_row
+        .iter()
+        .zip(dst_row)
+        .for_each(|(src_pixel, dst_pixel)| {
+            let components: [u16; 4] = src_pixel.0;
+            let alpha = components[3];
+            let recip_alpha = RECIP_ALPHA16[alpha as usize];
+            dst_pixel.0 = [
+                div_and_clip16(components[0], recip_alpha),
+                div_and_clip16(components[1], recip_alpha),
+                div_and_clip16(components[2], recip_alpha),
+                alpha,
+            ];
+        });
+}
@@ -0,0 +1,155 @@
+use std::arch::x86_64::*;
+
+use crate::image_view::{TypedImageView, TypedImageViewMut};
+use crate::pixels::U16x4;
+
+use super::native;
+
+#[target_feature(enable = "sse4.1")]
+pub(crate) unsafe fn multiply_alpha(
+    src_image: TypedImageView<U16x4>,
+    mut dst_image: TypedImageViewMut<U16x4>,
+) {
+    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 = "sse4.1")]
+pub(crate) unsafe fn multiply_alpha_inplace(mut image: TypedImageViewMut<U16x4>) {
+    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);
+    }
+}
+
+#[inline]
+#[target_feature(enable = "sse4.1")]
+pub(crate) unsafe fn multiply_alpha_row(src_row: &[U16x4], dst_row: &mut [U16x4]) {
+    let zero = _mm_setzero_si128();
+    let half = _mm_set1_epi32(0x8000);
+
+    const MAX_A: i64 = 0xffff000000000000u64 as i64;
+    let max_alpha = _mm_set1_epi64x(MAX_A);
+    /*
+       |R0   G0   B0   A0  | |R1   G1   B1   A1  |
+       |0001 0203 0405 0607| |0809 1011 1213 1415|
+    */
+    let factor_mask = _mm_set_epi8(15, 14, 15, 14, 15, 14, 15, 14, 7, 6, 7, 6, 7, 6, 7, 6);
+
+    let src_chunks = src_row.chunks_exact(2);
+    let src_remainder = src_chunks.remainder();
+    let mut dst_chunks = dst_row.chunks_exact_mut(2);
+
+    for (src, dst) in src_chunks.zip(&mut dst_chunks) {
+        let src_pixels = _mm_loadu_si128(src.as_ptr() as *const __m128i);
+
+        let factor_pixels = _mm_shuffle_epi8(src_pixels, factor_mask);
+        let factor_pixels = _mm_or_si128(factor_pixels, max_alpha);
+
+        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);
+    }
+}
+
+// Divide
+
+#[target_feature(enable = "sse4.1")]
+pub(crate) unsafe fn divide_alpha(
+    src_image: TypedImageView<U16x4>,
+    mut dst_image: TypedImageViewMut<U16x4>,
+) {
+    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 = "sse4.1")]
+pub(crate) unsafe fn divide_alpha_inplace(mut image: TypedImageViewMut<U16x4>) {
+    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);
+    }
+}
+
+#[target_feature(enable = "sse4.1")]
+pub(crate) unsafe fn divide_alpha_row(src_row: &[U16x4], dst_row: &mut [U16x4]) {
+    let src_chunks = src_row.chunks_exact(2);
+    let src_remainder = src_chunks.remainder();
+    let mut dst_chunks = dst_row.chunks_exact_mut(2);
+
+    for (src, dst) in src_chunks.zip(&mut dst_chunks) {
+        divide_alpha_two_pixels(src.as_ptr(), dst.as_mut_ptr());
+    }
+
+    if let Some(src) = src_remainder.get(0) {
+        let src_pixels = [*src, U16x4([0, 0, 0, 0])];
+        let mut dst_pixels = [U16x4([0, 0, 0, 0]); 2];
+        divide_alpha_two_pixels(src_pixels.as_ptr(), dst_pixels.as_mut_ptr());
+
+        let dst_reminder = dst_chunks.into_remainder();
+        if let Some(dst) = dst_reminder.get_mut(0) {
+            *dst = dst_pixels[0];
+        }
+    }
+}
+
+#[inline]
+#[target_feature(enable = "sse4.1")]
+unsafe fn divide_alpha_two_pixels(src: *const U16x4, dst: *mut U16x4) {
+    let zero = _mm_setzero_si128();
+    let alpha_mask = _mm_set1_epi64x(0xffff000000000000u64 as i64);
+    let alpha_max = _mm_set1_ps(65535.0);
+    /*
+       |R0   G0   B0   A0  | |R1   G1   B1   A1  |
+       |0001 0203 0405 0607| |0809 1011 1213 1415|
+    */
+    let alpha32_sh0 = _mm_set_epi8(-1, -1, 7, 6, -1, -1, 7, 6, -1, -1, 7, 6, -1, -1, 7, 6);
+    let alpha32_sh1 = _mm_set_epi8(
+        -1, -1, 15, 14, -1, -1, 15, 14, -1, -1, 15, 14, -1, -1, 15, 14,
+    );
+
+    let src_pixels = _mm_loadu_si128(src as *const __m128i);
+
+    let alpha0_f32x4 = _mm_cvtepi32_ps(_mm_shuffle_epi8(src_pixels, alpha32_sh0));
+    let alpha1_f32x4 = _mm_cvtepi32_ps(_mm_shuffle_epi8(src_pixels, alpha32_sh1));
+
+    let pix0_f32x4 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(src_pixels, zero));
+    let pix1_f32x4 = _mm_cvtepi32_ps(_mm_unpacklo_epi16(src_pixels, zero));
+
+    let scaled_pix0_f32x4 = _mm_mul_ps(pix0_f32x4, alpha_max);
+    let scaled_pix1_f32x4 = _mm_mul_ps(pix1_f32x4, alpha_max);
+
+    let divided_pix0_i32x4 = _mm_cvtps_epi32(_mm_div_ps(scaled_pix0_f32x4, alpha0_f32x4));
+    let divided_pix1_i32x4 = _mm_cvtps_epi32(_mm_div_ps(scaled_pix1_f32x4, alpha1_f32x4));
+
+    let two_pixels_i16x8 = _mm_packus_epi32(divided_pix0_i32x4, divided_pix1_i32x4);
+    let alpha = _mm_and_si128(src_pixels, alpha_mask);
+    let dst_pixels = _mm_blendv_epi8(two_pixels_i16x8, alpha, alpha_mask);
+
+    _mm_storeu_si128(dst as *mut __m128i, dst_pixels);
+}
@@ -15,6 +15,7 @@ mod optimisations;
 mod u16x1;
 mod u16x2;
 mod u16x3;
+mod u16x4;
 mod u8x1;
 mod u8x2;
 mod u8x3;
@@ -248,7 +248,7 @@ unsafe fn horiz_convolution_one_row(

    for (dst_x, coeffs_chunk) in coefficients_chunks.iter().enumerate() {
        let mut x: usize = coeffs_chunk.start as usize;
-        let mut ll_sum = _mm256_set1_epi64x(0);
+        let mut ll_sum = _mm256_setzero_si256();
        let mut coeffs = coeffs_chunk.values;

        let coeffs_by_8 = coeffs.chunks_exact(8);
@@ -331,6 +331,6 @@ unsafe fn horiz_convolution_one_row(
        dst_pixel.0 = [
            normalizer_guard.clip(ll_buf[0] + ll_buf[2] + half_error),
            normalizer_guard.clip(ll_buf[1] + ll_buf[3] + half_error),
-        ]
+        ];
    }
 }
@@ -149,9 +149,6 @@ unsafe fn horiz_convolution_four_rows(
            x += 2;
        }

-        let coeffs_by_2 = coeffs.chunks_exact(2);
-        coeffs = coeffs_by_2.remainder();
-
        if let Some(&k) = coeffs.get(0) {
            let coeff0_i64x2 = _mm_set1_epi64x(k as i64);
            for i in 0..4 {
@@ -167,7 +164,7 @@ unsafe fn horiz_convolution_four_rows(
            dst_pixel.0 = [
                normalizer_guard.clip(ll_buf[0]),
                normalizer_guard.clip(ll_buf[1]),
-            ]
+            ];
        }
    }
 }
@@ -275,6 +272,6 @@ unsafe fn horiz_convolution_one_row(
        dst_pixel.0 = [
            normalizer_guard.clip(ll_buf[0]),
            normalizer_guard.clip(ll_buf[1]),
-        ]
+        ];
    }
 }
@@ -0,0 +1,313 @@
+use std::arch::x86_64::*;
+
+use crate::convolution::{optimisations, Coefficients};
+use crate::image_view::{FourRows, FourRowsMut, TypedImageView, TypedImageViewMut};
+use crate::pixels::U16x4;
+use crate::simd_utils;
+
+#[inline]
+pub(crate) fn horiz_convolution(
+    src_image: TypedImageView<U16x4>,
+    mut dst_image: TypedImageViewMut<U16x4>,
+    offset: u32,
+    coeffs: Coefficients,
+) {
+    let (values, window_size, bounds_per_pixel) =
+        (coeffs.values, coeffs.window_size, coeffs.bounds);
+
+    let normalizer_guard = optimisations::NormalizerGuard32::new(values);
+    let coefficients_chunks = normalizer_guard.normalized_chunks(window_size, &bounds_per_pixel);
+    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_guard,
+            );
+        }
+    }
+
+    let mut yy = dst_height - dst_height % 4;
+    while yy < dst_height {
+        unsafe {
+            horiz_convolution_one_row(
+                src_image.get_row(yy + offset).unwrap(),
+                dst_image.get_row_mut(yy).unwrap(),
+                &coefficients_chunks,
+                &normalizer_guard,
+            );
+        }
+        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
+#[target_feature(enable = "avx2")]
+unsafe fn horiz_convolution_four_rows(
+    src_rows: FourRows<U16x4>,
+    dst_rows: FourRowsMut<U16x4>,
+    coefficients_chunks: &[optimisations::CoefficientsI32Chunk],
+    normalizer_guard: &optimisations::NormalizerGuard32,
+) {
+    let (s_row0, s_row1, s_row2, s_row3) = src_rows;
+    let s_rows = [s_row0, s_row1, s_row2, s_row3];
+    let (d_row0, d_row1, d_row2, d_row3) = dst_rows;
+    let d_rows = [d_row0, d_row1, d_row2, d_row3];
+    let precision = normalizer_guard.precision();
+    let half_error = 1i64 << (precision - 1);
+    let mut rg_buf = [0i64; 4];
+    let mut ba_buf = [0i64; 4];
+
+    /*
+       |R0   G0   B0   A0  | |R1   G1   B1   A1  |
+       |0001 0203 0405 0607| |0809 1011 1213 1415|
+
+        Shuffle to extract R0 and G0 as i64:
+        -1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0
+
+        Shuffle to extract R1 and G1 as i64:
+        -1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8
+
+        Shuffle to extract B0 and A0 as i64:
+        -1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4
+
+        Shuffle to extract B1 and A1 as i64:
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12
+    */
+    #[rustfmt::skip]
+    let rg0_shuffle = _mm256_set_epi8(
+        -1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0,
+        -1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0,
+    );
+    #[rustfmt::skip]
+    let rg1_shuffle = _mm256_set_epi8(
+        -1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8,
+        -1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8,
+    );
+    #[rustfmt::skip]
+    let ba0_shuffle = _mm256_set_epi8(
+        -1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4,
+        -1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4,
+    );
+    #[rustfmt::skip]
+    let ba1_shuffle = _mm256_set_epi8(
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12,
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12,
+    );
+
+    for (dst_x, coeffs_chunk) in coefficients_chunks.iter().enumerate() {
+        let mut x: usize = coeffs_chunk.start as usize;
+        let mut rg_sum = [_mm256_set1_epi64x(half_error); 2];
+        let mut ba_sum = [_mm256_set1_epi64x(half_error); 2];
+
+        let mut coeffs = coeffs_chunk.values;
+
+        let coeffs_by_2 = coeffs.chunks_exact(2);
+        coeffs = coeffs_by_2.remainder();
+
+        for k in coeffs_by_2 {
+            let coeff0_i64x4 = _mm256_set1_epi64x(k[0] as i64);
+            let coeff1_i64x4 = _mm256_set1_epi64x(k[1] as i64);
+
+            for i in 0..2 {
+                let source = _mm256_set_m128i(
+                    simd_utils::loadu_si128(s_rows[i * 2 + 1], x),
+                    simd_utils::loadu_si128(s_rows[i * 2], x),
+                );
+
+                let mut sum = rg_sum[i];
+                let rg_i64x4 = _mm256_shuffle_epi8(source, rg0_shuffle);
+                sum = _mm256_add_epi64(sum, _mm256_mul_epi32(rg_i64x4, coeff0_i64x4));
+                let rg_i64x4 = _mm256_shuffle_epi8(source, rg1_shuffle);
+                sum = _mm256_add_epi64(sum, _mm256_mul_epi32(rg_i64x4, coeff1_i64x4));
+                rg_sum[i] = sum;
+
+                let mut sum = ba_sum[i];
+                let ba_i64x4 = _mm256_shuffle_epi8(source, ba0_shuffle);
+                sum = _mm256_add_epi64(sum, _mm256_mul_epi32(ba_i64x4, coeff0_i64x4));
+                let ba_i64x4 = _mm256_shuffle_epi8(source, ba1_shuffle);
+                sum = _mm256_add_epi64(sum, _mm256_mul_epi32(ba_i64x4, coeff1_i64x4));
+                ba_sum[i] = sum;
+            }
+            x += 2;
+        }
+
+        if let Some(&k) = coeffs.get(0) {
+            let coeff0_i64x4 = _mm256_set1_epi64x(k as i64);
+
+            for i in 0..2 {
+                let source = _mm256_set_m128i(
+                    simd_utils::loadl_epi64(s_rows[i * 2 + 1], x),
+                    simd_utils::loadl_epi64(s_rows[i * 2], x),
+                );
+
+                let mut sum = rg_sum[i];
+                let rg_i64x4 = _mm256_shuffle_epi8(source, rg0_shuffle);
+                sum = _mm256_add_epi64(sum, _mm256_mul_epi32(rg_i64x4, coeff0_i64x4));
+                rg_sum[i] = sum;
+
+                let mut sum = ba_sum[i];
+                let ba_i64x4 = _mm256_shuffle_epi8(source, ba0_shuffle);
+                sum = _mm256_add_epi64(sum, _mm256_mul_epi32(ba_i64x4, coeff0_i64x4));
+                ba_sum[i] = sum;
+            }
+        }
+
+        for i in 0..2 {
+            _mm256_storeu_si256((&mut rg_buf).as_mut_ptr() as *mut __m256i, rg_sum[i]);
+            _mm256_storeu_si256((&mut ba_buf).as_mut_ptr() as *mut __m256i, ba_sum[i]);
+
+            let dst_pixel = d_rows[i * 2].get_unchecked_mut(dst_x);
+            dst_pixel.0 = [
+                normalizer_guard.clip(rg_buf[0]),
+                normalizer_guard.clip(rg_buf[1]),
+                normalizer_guard.clip(ba_buf[0]),
+                normalizer_guard.clip(ba_buf[1]),
+            ];
+
+            let dst_pixel = d_rows[i * 2 + 1].get_unchecked_mut(dst_x);
+            dst_pixel.0 = [
+                normalizer_guard.clip(rg_buf[2]),
+                normalizer_guard.clip(rg_buf[3]),
+                normalizer_guard.clip(ba_buf[2]),
+                normalizer_guard.clip(ba_buf[3]),
+            ];
+        }
+    }
+}
+
+/// 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 = "avx2")]
+unsafe fn horiz_convolution_one_row(
+    src_row: &[U16x4],
+    dst_row: &mut [U16x4],
+    coefficients_chunks: &[optimisations::CoefficientsI32Chunk],
+    normalizer_guard: &optimisations::NormalizerGuard32,
+) {
+    let precision = normalizer_guard.precision();
+    let half_error = 1i64 << (precision - 1);
+    let mut rg_buf = [0i64; 4];
+    let mut ba_buf = [0i64; 4];
+
+    /*
+       |R0   G0   B0   A0  | |R1   G1   B1   A1  |
+       |0001 0203 0405 0607| |0809 1011 1213 1415|
+
+        Shuffle to extract R0 and G0 as i64:
+        -1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0
+
+        Shuffle to extract R1 and G1 as i64:
+        -1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8
+
+        Shuffle to extract B0 and A0 as i64:
+        -1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4
+
+        Shuffle to extract B1 and A1 as i64:
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12
+    */
+
+    #[rustfmt::skip]
+    let rg02_shuffle = _mm256_set_epi8(
+        -1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0,
+        -1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0,
+    );
+    #[rustfmt::skip]
+    let rg13_shuffle = _mm256_set_epi8(
+        -1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8,
+        -1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8,
+    );
+    #[rustfmt::skip]
+    let ba02_shuffle = _mm256_set_epi8(
+        -1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4,
+        -1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4,
+    );
+    #[rustfmt::skip]
+    let ba13_shuffle = _mm256_set_epi8(
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12,
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12,
+    );
+
+    for (dst_x, coeffs_chunk) in coefficients_chunks.iter().enumerate() {
+        let mut x: usize = coeffs_chunk.start as usize;
+        let mut coeffs = coeffs_chunk.values;
+        let mut rg_sum = _mm256_setzero_si256();
+        let mut ba_sum = _mm256_setzero_si256();
+
+        let coeffs_by_4 = coeffs.chunks_exact(4);
+        coeffs = coeffs_by_4.remainder();
+
+        for k in coeffs_by_4 {
+            let coeff02_i64x4 =
+                _mm256_set_epi64x(k[2] as i64, k[2] as i64, k[0] as i64, k[0] as i64);
+            let coeff13_i64x4 =
+                _mm256_set_epi64x(k[3] as i64, k[3] as i64, k[1] as i64, k[1] as i64);
+
+            let source = simd_utils::loadu_si256(src_row, x);
+
+            let rg_i64x4 = _mm256_shuffle_epi8(source, rg02_shuffle);
+            rg_sum = _mm256_add_epi64(rg_sum, _mm256_mul_epi32(rg_i64x4, coeff02_i64x4));
+            let rg_i64x4 = _mm256_shuffle_epi8(source, rg13_shuffle);
+            rg_sum = _mm256_add_epi64(rg_sum, _mm256_mul_epi32(rg_i64x4, coeff13_i64x4));
+
+            let ba_i64x4 = _mm256_shuffle_epi8(source, ba02_shuffle);
+            ba_sum = _mm256_add_epi64(ba_sum, _mm256_mul_epi32(ba_i64x4, coeff02_i64x4));
+            let ba_i64x4 = _mm256_shuffle_epi8(source, ba13_shuffle);
+            ba_sum = _mm256_add_epi64(ba_sum, _mm256_mul_epi32(ba_i64x4, coeff13_i64x4));
+            x += 4;
+        }
+
+        let coeffs_by_2 = coeffs.chunks_exact(2);
+        coeffs = coeffs_by_2.remainder();
+        for k in coeffs_by_2 {
+            let coeff01_i64x4 =
+                _mm256_set_epi64x(k[1] as i64, k[1] as i64, k[0] as i64, k[0] as i64);
+
+            let source = _mm256_set_m128i(
+                simd_utils::loadl_epi64(src_row, x + 1),
+                simd_utils::loadl_epi64(src_row, x),
+            );
+
+            let rg_i64x4 = _mm256_shuffle_epi8(source, rg02_shuffle);
+            rg_sum = _mm256_add_epi64(rg_sum, _mm256_mul_epi32(rg_i64x4, coeff01_i64x4));
+            let ba_i64x4 = _mm256_shuffle_epi8(source, ba02_shuffle);
+            ba_sum = _mm256_add_epi64(ba_sum, _mm256_mul_epi32(ba_i64x4, coeff01_i64x4));
+
+            x += 2;
+        }
+
+        if let Some(&k) = coeffs.get(0) {
+            let coeff_i64x4 = _mm256_set_epi64x(0, 0, k as i64, k as i64);
+            let source = _mm256_set_m128i(_mm_setzero_si128(), simd_utils::loadl_epi64(src_row, x));
+
+            let rg_i64x4 = _mm256_shuffle_epi8(source, rg02_shuffle);
+            rg_sum = _mm256_add_epi64(rg_sum, _mm256_mul_epi32(rg_i64x4, coeff_i64x4));
+            let ba_i64x4 = _mm256_shuffle_epi8(source, ba02_shuffle);
+            ba_sum = _mm256_add_epi64(ba_sum, _mm256_mul_epi32(ba_i64x4, coeff_i64x4));
+        }
+
+        _mm256_storeu_si256((&mut rg_buf).as_mut_ptr() as *mut __m256i, rg_sum);
+        _mm256_storeu_si256((&mut ba_buf).as_mut_ptr() as *mut __m256i, ba_sum);
+        let dst_pixel = dst_row.get_unchecked_mut(dst_x);
+        dst_pixel.0 = [
+            normalizer_guard.clip(rg_buf[0] + rg_buf[2] + half_error),
+            normalizer_guard.clip(rg_buf[1] + rg_buf[3] + half_error),
+            normalizer_guard.clip(ba_buf[0] + ba_buf[2] + half_error),
+            normalizer_guard.clip(ba_buf[1] + ba_buf[3] + half_error),
+        ];
+    }
+}
@@ -0,0 +1,38 @@
+use super::{Coefficients, Convolution};
+use crate::convolution::vertical_u16::vert_convolution_u16;
+use crate::image_view::{TypedImageView, TypedImageViewMut};
+use crate::pixels::U16x4;
+use crate::CpuExtensions;
+
+#[cfg(target_arch = "x86_64")]
+mod avx2;
+mod native;
+#[cfg(target_arch = "x86_64")]
+mod sse4;
+
+impl Convolution for U16x4 {
+    fn horiz_convolution(
+        src_image: TypedImageView<Self>,
+        dst_image: TypedImageViewMut<Self>,
+        offset: u32,
+        coeffs: Coefficients,
+        cpu_extensions: CpuExtensions,
+    ) {
+        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),
+        }
+    }
+
+    fn vert_convolution(
+        src_image: TypedImageView<Self>,
+        dst_image: TypedImageViewMut<Self>,
+        coeffs: Coefficients,
+        cpu_extensions: CpuExtensions,
+    ) {
+        vert_convolution_u16(src_image, dst_image, coeffs, cpu_extensions);
+    }
+}
@@ -0,0 +1,36 @@
+use crate::convolution::{optimisations, Coefficients};
+use crate::image_view::{TypedImageView, TypedImageViewMut};
+use crate::pixels::U16x4;
+
+#[inline(always)]
+pub(crate) fn horiz_convolution(
+    src_image: TypedImageView<U16x4>,
+    mut dst_image: TypedImageViewMut<U16x4>,
+    offset: u32,
+    coeffs: Coefficients,
+) {
+    let (values, window_size, bounds) = (coeffs.values, coeffs.window_size, coeffs.bounds);
+
+    let normalizer_guard = optimisations::NormalizerGuard32::new(values);
+    let precision = normalizer_guard.precision();
+    let coefficients_chunks = normalizer_guard.normalized_chunks(window_size, &bounds);
+    let initial: i64 = 1 << (precision - 1);
+
+    let src_rows = src_image.iter_rows(offset);
+    let dst_rows = dst_image.iter_rows_mut();
+    for (dst_row, src_row) in dst_rows.zip(src_rows) {
+        for (&coeffs_chunk, dst_pixel) in coefficients_chunks.iter().zip(dst_row.iter_mut()) {
+            let first_x_src = coeffs_chunk.start as usize;
+            let mut ss = [initial; 4];
+            let src_pixels = unsafe { src_row.get_unchecked(first_x_src..) };
+            for (&k, src_pixel) in coeffs_chunk.values.iter().zip(src_pixels) {
+                for (i, s) in ss.iter_mut().enumerate() {
+                    *s += src_pixel.0[i] as i64 * (k as i64);
+                }
+            }
+            for (i, s) in ss.iter().copied().enumerate() {
+                dst_pixel.0[i] = normalizer_guard.clip(s);
+            }
+        }
+    }
+}
@@ -0,0 +1,242 @@
+use std::arch::x86_64::*;
+
+use crate::convolution::{optimisations, Coefficients};
+use crate::image_view::{FourRows, FourRowsMut, TypedImageView, TypedImageViewMut};
+use crate::pixels::U16x4;
+use crate::simd_utils;
+
+#[inline]
+pub(crate) fn horiz_convolution(
+    src_image: TypedImageView<U16x4>,
+    mut dst_image: TypedImageViewMut<U16x4>,
+    offset: u32,
+    coeffs: Coefficients,
+) {
+    let (values, window_size, bounds_per_pixel) =
+        (coeffs.values, coeffs.window_size, coeffs.bounds);
+
+    let normalizer_guard = optimisations::NormalizerGuard32::new(values);
+    let coefficients_chunks = normalizer_guard.normalized_chunks(window_size, &bounds_per_pixel);
+    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_guard,
+            );
+        }
+    }
+
+    let mut yy = dst_height - dst_height % 4;
+    while yy < dst_height {
+        unsafe {
+            horiz_convolution_one_row(
+                src_image.get_row(yy + offset).unwrap(),
+                dst_image.get_row_mut(yy).unwrap(),
+                &coefficients_chunks,
+                &normalizer_guard,
+            );
+        }
+        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
+#[target_feature(enable = "sse4.1")]
+unsafe fn horiz_convolution_four_rows(
+    src_rows: FourRows<U16x4>,
+    dst_rows: FourRowsMut<U16x4>,
+    coefficients_chunks: &[optimisations::CoefficientsI32Chunk],
+    normalizer_guard: &optimisations::NormalizerGuard32,
+) {
+    let (s_row0, s_row1, s_row2, s_row3) = src_rows;
+    let s_rows = [s_row0, s_row1, s_row2, s_row3];
+    let (d_row0, d_row1, d_row2, d_row3) = dst_rows;
+    let d_rows = [d_row0, d_row1, d_row2, d_row3];
+    let precision = normalizer_guard.precision();
+    let half_error = 1i64 << (precision - 1);
+    let mut rg_buf = [0i64; 2];
+    let mut ba_buf = [0i64; 2];
+
+    /*
+       |R0   G0   B0   A0  | |R1   G1   B1   A1  |
+       |0001 0203 0405 0607| |0809 1011 1213 1415|
+
+        Shuffle to extract R0 and G0 as i64:
+        -1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0
+
+        Shuffle to extract R1 and G1 as i64:
+        -1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8
+
+        Shuffle to extract B0 and A0 as i64:
+        -1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4
+
+        Shuffle to extract B1 and A1 as i64:
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12
+    */
+
+    let rg0_shuffle = _mm_set_epi8(-1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0);
+    let rg1_shuffle = _mm_set_epi8(-1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8);
+    let ba0_shuffle = _mm_set_epi8(-1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4);
+    let ba1_shuffle = _mm_set_epi8(
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12,
+    );
+
+    for (dst_x, coeffs_chunk) in coefficients_chunks.iter().enumerate() {
+        let mut x: usize = coeffs_chunk.start as usize;
+        let mut rg_sum = [_mm_set1_epi64x(half_error); 4];
+        let mut ba_sum = [_mm_set1_epi64x(half_error); 4];
+
+        let mut coeffs = coeffs_chunk.values;
+
+        let coeffs_by_2 = coeffs.chunks_exact(2);
+        coeffs = coeffs_by_2.remainder();
+
+        for k in coeffs_by_2 {
+            let coeff0_i64x2 = _mm_set1_epi64x(k[0] as i64);
+            let coeff1_i64x2 = _mm_set1_epi64x(k[1] as i64);
+
+            for i in 0..4 {
+                let source = simd_utils::loadu_si128(s_rows[i], x);
+                let mut sum = rg_sum[i];
+                let rg_i64x2 = _mm_shuffle_epi8(source, rg0_shuffle);
+                sum = _mm_add_epi64(sum, _mm_mul_epi32(rg_i64x2, coeff0_i64x2));
+                let rg_i64x2 = _mm_shuffle_epi8(source, rg1_shuffle);
+                sum = _mm_add_epi64(sum, _mm_mul_epi32(rg_i64x2, coeff1_i64x2));
+                rg_sum[i] = sum;
+
+                let mut sum = ba_sum[i];
+                let ba_i64x2 = _mm_shuffle_epi8(source, ba0_shuffle);
+                sum = _mm_add_epi64(sum, _mm_mul_epi32(ba_i64x2, coeff0_i64x2));
+                let ba_i64x2 = _mm_shuffle_epi8(source, ba1_shuffle);
+                sum = _mm_add_epi64(sum, _mm_mul_epi32(ba_i64x2, coeff1_i64x2));
+                ba_sum[i] = sum;
+            }
+            x += 2;
+        }
+
+        if let Some(&k) = coeffs.get(0) {
+            let coeff0_i64x2 = _mm_set1_epi64x(k as i64);
+            for i in 0..4 {
+                let source = simd_utils::loadl_epi64(s_rows[i], x);
+                let rg_i64x2 = _mm_shuffle_epi8(source, rg0_shuffle);
+                rg_sum[i] = _mm_add_epi64(rg_sum[i], _mm_mul_epi32(rg_i64x2, coeff0_i64x2));
+                let ba_i64x2 = _mm_shuffle_epi8(source, ba0_shuffle);
+                ba_sum[i] = _mm_add_epi64(ba_sum[i], _mm_mul_epi32(ba_i64x2, coeff0_i64x2));
+            }
+        }
+
+        for i in 0..4 {
+            _mm_storeu_si128((&mut rg_buf).as_mut_ptr() as *mut __m128i, rg_sum[i]);
+            _mm_storeu_si128((&mut ba_buf).as_mut_ptr() as *mut __m128i, ba_sum[i]);
+            let dst_pixel = d_rows[i].get_unchecked_mut(dst_x);
+            dst_pixel.0 = [
+                normalizer_guard.clip(rg_buf[0]),
+                normalizer_guard.clip(rg_buf[1]),
+                normalizer_guard.clip(ba_buf[0]),
+                normalizer_guard.clip(ba_buf[1]),
+            ];
+        }
+    }
+}
+
+/// 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_one_row(
+    src_row: &[U16x4],
+    dst_row: &mut [U16x4],
+    coefficients_chunks: &[optimisations::CoefficientsI32Chunk],
+    normalizer_guard: &optimisations::NormalizerGuard32,
+) {
+    let precision = normalizer_guard.precision();
+    let half_error = 1i64 << (precision - 1);
+    let mut rg_buf = [0i64; 2];
+    let mut ba_buf = [0i64; 2];
+
+    /*
+       |R0   G0   B0   A0  | |R1   G1   B1   A1  |
+       |0001 0203 0405 0607| |0809 1011 1213 1415|
+
+        Shuffle to extract R0 and G0 as i64:
+        -1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0
+
+        Shuffle to extract R1 and G1 as i64:
+        -1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8
+
+        Shuffle to extract B0 and A0 as i64:
+        -1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4
+
+        Shuffle to extract B1 and A1 as i64:
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12
+    */
+
+    let rg0_shuffle = _mm_set_epi8(-1, -1, -1, -1, -1, -1, 3, 2, -1, -1, -1, -1, -1, -1, 1, 0);
+    let rg1_shuffle = _mm_set_epi8(-1, -1, -1, -1, -1, -1, 11, 10, -1, -1, -1, -1, -1, -1, 9, 8);
+    let ba0_shuffle = _mm_set_epi8(-1, -1, -1, -1, -1, -1, 7, 6, -1, -1, -1, -1, -1, -1, 5, 4);
+    let ba1_shuffle = _mm_set_epi8(
+        -1, -1, -1, -1, -1, -1, 15, 14, -1, -1, -1, -1, -1, -1, 13, 12,
+    );
+
+    for (dst_x, coeffs_chunk) in coefficients_chunks.iter().enumerate() {
+        let mut x: usize = coeffs_chunk.start as usize;
+        let mut coeffs = coeffs_chunk.values;
+        let mut rg_sum = _mm_set1_epi64x(half_error);
+        let mut ba_sum = _mm_set1_epi64x(half_error);
+
+        let coeffs_by_2 = coeffs.chunks_exact(2);
+        coeffs = coeffs_by_2.remainder();
+
+        for k in coeffs_by_2 {
+            let coeff0_i64x2 = _mm_set1_epi64x(k[0] as i64);
+            let coeff1_i64x2 = _mm_set1_epi64x(k[1] as i64);
+
+            let source = simd_utils::loadu_si128(src_row, x);
+
+            let rg_i64x2 = _mm_shuffle_epi8(source, rg0_shuffle);
+            rg_sum = _mm_add_epi64(rg_sum, _mm_mul_epi32(rg_i64x2, coeff0_i64x2));
+            let rg_i64x2 = _mm_shuffle_epi8(source, rg1_shuffle);
+            rg_sum = _mm_add_epi64(rg_sum, _mm_mul_epi32(rg_i64x2, coeff1_i64x2));
+
+            let ba_i64x2 = _mm_shuffle_epi8(source, ba0_shuffle);
+            ba_sum = _mm_add_epi64(ba_sum, _mm_mul_epi32(ba_i64x2, coeff0_i64x2));
+            let ba_i64x2 = _mm_shuffle_epi8(source, ba1_shuffle);
+            ba_sum = _mm_add_epi64(ba_sum, _mm_mul_epi32(ba_i64x2, coeff1_i64x2));
+
+            x += 2;
+        }
+
+        if let Some(&k) = coeffs.get(0) {
+            let coeff0_i64x2 = _mm_set1_epi64x(k as i64);
+            let source = simd_utils::loadl_epi64(src_row, x);
+            let rg_i64x2 = _mm_shuffle_epi8(source, rg0_shuffle);
+            rg_sum = _mm_add_epi64(rg_sum, _mm_mul_epi32(rg_i64x2, coeff0_i64x2));
+            let ba_i64x2 = _mm_shuffle_epi8(source, ba0_shuffle);
+            ba_sum = _mm_add_epi64(ba_sum, _mm_mul_epi32(ba_i64x2, coeff0_i64x2));
+        }
+
+        _mm_storeu_si128((&mut rg_buf).as_mut_ptr() as *mut __m128i, rg_sum);
+        _mm_storeu_si128((&mut ba_buf).as_mut_ptr() as *mut __m128i, ba_sum);
+        let dst_pixel = dst_row.get_unchecked_mut(dst_x);
+        dst_pixel.0 = [
+            normalizer_guard.clip(rg_buf[0]),
+            normalizer_guard.clip(rg_buf[1]),
+            normalizer_guard.clip(ba_buf[0]),
+            normalizer_guard.clip(ba_buf[1]),
+        ];
+    }
+}
@@ -1,7 +1,7 @@
 use std::num::NonZeroU32;

 use crate::image_view::{ImageRows, ImageRowsMut, TypedImageView, TypedImageViewMut};
-use crate::pixels::{Pixel, PixelType, U16x2, U16x3, U8x2, U8x3, U8x4, F32, I32, U16, U8};
+use crate::pixels::{Pixel, PixelType, U16x2, U16x3, U16x4, U8x2, U8x3, U8x4, F32, I32, U16, U8};
 use crate::{ImageBufferError, ImageView, ImageViewMut};

 #[derive(Debug)]
@@ -174,6 +174,15 @@ impl<'a> Image<'a> {
                        .collect(),
                )
            }
+            PixelType::U16x4 => {
+                let pixels = unsafe { buffer.align_to::<U16x4>().1 };
+                ImageRows::U16x4(
+                    pixels
+                        .chunks_exact(self.width.get() as usize)
+                        .take(rows_count)
+                        .collect(),
+                )
+            }
            PixelType::I32 => {
                let pixels = unsafe { buffer.align_to::<I32>().1 };
                ImageRows::I32(
@@ -267,6 +276,15 @@ impl<'a> Image<'a> {
                        .collect(),
                )
            }
+            PixelType::U16x4 => {
+                let pixels = unsafe { buffer.align_to_mut::<U16x4>().1 };
+                ImageRowsMut::U16x4(
+                    pixels
+                        .chunks_exact_mut(width.get() as usize)
+                        .take(rows_count)
+                        .collect(),
+                )
+            }
            PixelType::I32 => {
                let pixels = unsafe { buffer.align_to_mut::<I32>().1 };
                ImageRowsMut::I32(
@@ -2,7 +2,7 @@ use std::num::NonZeroU32;
 use std::slice;

 use crate::errors::{CropBoxError, ImageBufferError, ImageRowsError};
-use crate::pixels::{Pixel, PixelType, U16x2, U16x3, U8x2, U8x3, U8x4, F32, I32, U16, U8};
+use crate::pixels::{Pixel, PixelType, U16x2, U16x3, U16x4, U8x2, U8x3, U8x4, F32, I32, U16, U8};

 pub(crate) type RowMut<'a, 'b, T> = &'a mut &'b mut [T];
 pub(crate) type TwoRows<'a, T> = (&'a [T], &'a [T]);
@@ -34,6 +34,7 @@ pub enum ImageRows<'a> {
    U16(Vec<&'a [U16]>),
    U16x2(Vec<&'a [U16x2]>),
    U16x3(Vec<&'a [U16x3]>),
+    U16x4(Vec<&'a [U16x4]>),
    I32(Vec<&'a [I32]>),
    F32(Vec<&'a [F32]>),
 }
@@ -52,6 +53,7 @@ impl<'a> ImageRows<'a> {
            ImageRows::U16(rows) => check_rows_count_and_size(width, height, rows),
            ImageRows::U16x2(rows) => check_rows_count_and_size(width, height, rows),
            ImageRows::U16x3(rows) => check_rows_count_and_size(width, height, rows),
+            ImageRows::U16x4(rows) => check_rows_count_and_size(width, height, rows),
            ImageRows::I32(rows) => check_rows_count_and_size(width, height, rows),
            ImageRows::F32(rows) => check_rows_count_and_size(width, height, rows),
        }
@@ -66,6 +68,7 @@ impl<'a> ImageRows<'a> {
            Self::U16(_) => PixelType::U16,
            Self::U16x2(_) => PixelType::U16x2,
            Self::U16x3(_) => PixelType::U16x3,
+            Self::U16x4(_) => PixelType::U16x4,
            Self::I32(_) => PixelType::I32,
            Self::F32(_) => PixelType::F32,
        }
@@ -89,6 +92,7 @@ image_rows_from!(U8x4, ImageRows::U8x4);
 image_rows_from!(U16, ImageRows::U16);
 image_rows_from!(U16x2, ImageRows::U16x2);
 image_rows_from!(U16x3, ImageRows::U16x3);
+image_rows_from!(U16x4, ImageRows::U16x4);
 image_rows_from!(I32, ImageRows::I32);
 image_rows_from!(F32, ImageRows::F32);

@@ -102,6 +106,7 @@ pub enum ImageRowsMut<'a> {
    U16(Vec<&'a mut [U16]>),
    U16x2(Vec<&'a mut [U16x2]>),
    U16x3(Vec<&'a mut [U16x3]>),
+    U16x4(Vec<&'a mut [U16x4]>),
    I32(Vec<&'a mut [I32]>),
    F32(Vec<&'a mut [F32]>),
    U8(Vec<&'a mut [U8]>),
@@ -120,6 +125,7 @@ impl<'a> ImageRowsMut<'a> {
            Self::U16(rows) => check_rows_count_and_size(width, height, rows),
            Self::U16x2(rows) => check_rows_count_and_size(width, height, rows),
            Self::U16x3(rows) => check_rows_count_and_size(width, height, rows),
+            Self::U16x4(rows) => check_rows_count_and_size(width, height, rows),
            Self::I32(rows) => check_rows_count_and_size(width, height, rows),
            Self::F32(rows) => check_rows_count_and_size(width, height, rows),
            Self::U8(rows) => check_rows_count_and_size(width, height, rows),
@@ -134,6 +140,7 @@ impl<'a> ImageRowsMut<'a> {
            Self::U16(_) => PixelType::U16,
            Self::U16x2(_) => PixelType::U16x2,
            Self::U16x3(_) => PixelType::U16x3,
+            Self::U16x4(_) => PixelType::U16x4,
            Self::I32(_) => PixelType::I32,
            Self::F32(_) => PixelType::F32,
            Self::U8(_) => PixelType::U8,
@@ -236,6 +243,15 @@ impl<'a> ImageView<'a> {
                        .collect(),
                )
            }
+            PixelType::U16x4 => {
+                let pixels = align_buffer_to(buffer)?;
+                ImageRows::U16x4(
+                    pixels
+                        .chunks_exact(width.get() as usize)
+                        .take(rows_count)
+                        .collect(),
+                )
+            }
            PixelType::I32 => {
                let pixels = align_buffer_to(buffer)?;
                ImageRows::I32(
@@ -450,6 +466,19 @@ impl<'a> ImageView<'a> {
        }
    }

+    pub(crate) fn u16x4_image(&self) -> Option<TypedImageView<U16x4>> {
+        if let ImageRows::U16x4(ref rows) = self.rows {
+            Some(TypedImageView {
+                width: self.width,
+                height: self.height,
+                crop_box: self.crop_box,
+                rows,
+            })
+        } else {
+            None
+        }
+    }
+
    pub(crate) fn i32_image(&self) -> Option<TypedImageView<I32>> {
        if let ImageRows::I32(ref rows) = self.rows {
            Some(TypedImageView {
@@ -682,6 +711,15 @@ impl<'a> ImageViewMut<'a> {
                        .collect(),
                )
            }
+            PixelType::U16x4 => {
+                let pixels = align_buffer_to_mut(buffer)?;
+                ImageRowsMut::U16x4(
+                    pixels
+                        .chunks_exact_mut(width.get() as usize)
+                        .take(rows_count)
+                        .collect(),
+                )
+            }
            PixelType::I32 => {
                let pixels = align_buffer_to_mut(buffer)?;
                ImageRowsMut::I32(
@@ -804,6 +842,18 @@ impl<'a> ImageViewMut<'a> {
        }
    }

+    pub(crate) fn u16x4_image<'s>(&'s mut self) -> Option<TypedImageViewMut<'s, 'a, U16x4>> {
+        if let ImageRowsMut::U16x4(rows) = &mut self.rows {
+            Some(TypedImageViewMut {
+                width: self.width,
+                height: self.height,
+                rows,
+            })
+        } else {
+            None
+        }
+    }
+
    pub(crate) fn i32_image<'s>(&'s mut self) -> Option<TypedImageViewMut<'s, 'a, I32>> {
        if let ImageRowsMut::I32(rows) = &mut self.rows {
            Some(TypedImageViewMut {
@@ -1,12 +1,12 @@
 use crate::alpha::AlphaMulDiv;
 use crate::image_view::{TypedImageView, TypedImageViewMut};
-use crate::pixels::{U16x2, U8x2, U8x4};
+use crate::pixels::{U16x2, U16x4, U8x2, U8x4};
 use crate::{
    CpuExtensions, ImageView, ImageViewMut, MulDivImageError, MulDivImagesError, PixelType,
 };

 /// Methods of this structure used to multiply or divide color-channels (RGB or Luma)
-/// by alpha-channel. Supported pixel types: U8x2, U8x4.
+/// by alpha-channel. Supported pixel types: U8x2, U8x4, U16x2 and U16x4.
 ///
 /// By default, instance of `MulDiv` created with best CPU-extensions provided by your CPU.
 /// You can change this by use method [MulDiv::set_cpu_extensions].
@@ -67,6 +67,12 @@ impl MulDiv {
                multiply_alpha(typed_src_image, typed_dst_image, self.cpu_extensions);
                Ok(())
            }
+            PixelType::U16x4 => {
+                let (typed_src_image, typed_dst_image) = assert_images_u16x4(src_image, dst_image)?;
+                multiply_alpha(typed_src_image, typed_dst_image, self.cpu_extensions);
+                Ok(())
+            }
+
            _ => Err(MulDivImagesError::UnsupportedPixelType),
        }
    }
@@ -89,6 +95,11 @@ impl MulDiv {
                multiply_alpha_inplace(typed_image, self.cpu_extensions);
                Ok(())
            }
+            PixelType::U16x4 => {
+                let typed_image = assert_image_u16x4(image)?;
+                multiply_alpha_inplace(typed_image, self.cpu_extensions);
+                Ok(())
+            }
            _ => Err(MulDivImageError::UnsupportedPixelType),
        }
    }
@@ -116,6 +127,11 @@ impl MulDiv {
                divide_alpha(typed_src_image, typed_dst_image, self.cpu_extensions);
                Ok(())
            }
+            PixelType::U16x4 => {
+                let (typed_src_image, typed_dst_image) = assert_images_u16x4(src_image, dst_image)?;
+                divide_alpha(typed_src_image, typed_dst_image, self.cpu_extensions);
+                Ok(())
+            }
            _ => Err(MulDivImagesError::UnsupportedPixelType),
        }
    }
@@ -138,6 +154,11 @@ impl MulDiv {
                divide_alpha_inplace(typed_image, self.cpu_extensions);
                Ok(())
            }
+            PixelType::U16x4 => {
+                let typed_image = assert_image_u16x4(image)?;
+                divide_alpha_inplace(typed_image, self.cpu_extensions);
+                Ok(())
+            }
            _ => Err(MulDivImageError::UnsupportedPixelType),
        }
    }
@@ -245,6 +266,40 @@ fn assert_image_u16x2<'a, 'b>(
        .ok_or(MulDivImageError::UnsupportedPixelType)
 }

+#[inline]
+fn assert_images_u16x4<'s, 'd, 'da>(
+    src_image: &'s ImageView<'s>,
+    dst_image: &'d mut ImageViewMut<'da>,
+) -> Result<
+    (
+        TypedImageView<'s, 's, U16x4>,
+        TypedImageViewMut<'d, 'da, U16x4>,
+    ),
+    MulDivImagesError,
+> {
+    let src_image_u16x4 = src_image
+        .u16x4_image()
+        .ok_or(MulDivImagesError::UnsupportedPixelType)?;
+    let dst_image_u16x4 = dst_image
+        .u16x4_image()
+        .ok_or(MulDivImagesError::UnsupportedPixelType)?;
+    if src_image_u16x4.width() != dst_image_u16x4.width()
+        || src_image_u16x4.height() != dst_image_u16x4.height()
+    {
+        return Err(MulDivImagesError::SizeIsDifferent);
+    }
+    Ok((src_image_u16x4, dst_image_u16x4))
+}
+
+#[inline]
+fn assert_image_u16x4<'a, 'b>(
+    image: &'a mut ImageViewMut<'b>,
+) -> Result<TypedImageViewMut<'a, 'b, U16x4>, MulDivImageError> {
+    image
+        .u16x4_image()
+        .ok_or(MulDivImageError::UnsupportedPixelType)
+}
+
 fn multiply_alpha<P>(
    src_image: TypedImageView<P>,
    dst_image: TypedImageViewMut<P>,
@@ -12,6 +12,7 @@ pub enum PixelType {
    U16,
    U16x2,
    U16x3,
+    U16x4,
    I32,
    F32,
    U8,
@@ -26,6 +27,7 @@ impl PixelType {
            Self::U16 => 2,
            Self::U16x2 => 4,
            Self::U16x3 => 6,
+            Self::U16x4 => 8,
            _ => 4,
        }
    }
@@ -40,6 +42,7 @@ impl PixelType {
            Self::U16 => unsafe { buffer.align_to::<U16>().0.is_empty() },
            Self::U16x2 => unsafe { buffer.align_to::<U16x2>().0.is_empty() },
            Self::U16x3 => unsafe { buffer.align_to::<U16x3>().0.is_empty() },
+            Self::U16x4 => unsafe { buffer.align_to::<U16x4>().0.is_empty() },
            Self::I32 => unsafe { buffer.align_to::<I32>().0.is_empty() },
            Self::F32 => unsafe { buffer.align_to::<F32>().0.is_empty() },
        }
@@ -108,14 +111,14 @@ macro_rules! pixel_struct {
    };
 }

-pixel_struct!(U8, u8, u8, 1, PixelType::U8, "One byte per pixel");
+pixel_struct!(U8, u8, u8, 1, PixelType::U8, "One byte per pixel (e.g. L8)");
 pixel_struct!(
    U8x2,
    u16,
    u8,
    2,
    PixelType::U8x2,
-    "Two bytes per pixel (e.g. LA)"
+    "Two bytes per pixel (e.g. LA8)"
 );
 pixel_struct!(
    U8x3,
@@ -123,7 +126,7 @@ pixel_struct!(
    u8,
    3,
    PixelType::U8x3,
-    "Three bytes per pixel (e.g. RGB)"
+    "Three bytes per pixel (e.g. RGB8)"
 );
 pixel_struct!(
    U8x4,
@@ -131,7 +134,7 @@ pixel_struct!(
    u8,
    4,
    PixelType::U8x4,
-    "Four bytes per pixel (RGBA, RGBx, CMYK and other)"
+    "Four bytes per pixel (RGBA8, RGBx8, CMYK8 and other)"
 );
 pixel_struct!(
    U16,
@@ -147,7 +150,7 @@ pixel_struct!(
    u16,
    2,
    PixelType::U16x2,
-    "Two `u16` components per pixel (e.g. LA)"
+    "Two `u16` components per pixel (e.g. LA16)"
 );
 pixel_struct!(
    U16x3,
@@ -155,7 +158,15 @@ pixel_struct!(
    u16,
    3,
    PixelType::U16x3,
-    "Three `u16` components per pixel (e.g. RGB)"
+    "Three `u16` components per pixel (e.g. RGB16)"
+);
+pixel_struct!(
+    U16x4,
+    [u16; 4],
+    u16,
+    4,
+    PixelType::U16x4,
+    "Four `u16` components per pixel (e.g. RGBA16)"
 );
 pixel_struct!(
    I32,
@@ -125,6 +125,13 @@ impl Resizer {
                    }
                }
            }
+            PixelType::U16x4 => {
+                if let Some(src_rows) = src_image.u16x4_image() {
+                    if let Some(dst_rows) = dst_image.u16x4_image() {
+                        self.resize_inner(src_rows, dst_rows);
+                    }
+                }
+            }
            PixelType::I32 => {
                if let Some(src_rows) = src_image.i32_image() {
                    if let Some(dst_rows) = dst_image.i32_image() {
@@ -33,6 +33,7 @@ pub trait PixelExt: Pixel {
            PixelType::U16 => "u16",
            PixelType::U16x2 => "u16x2",
            PixelType::U16x3 => "u16x3",
+            PixelType::U16x4 => "u16x4",
            PixelType::I32 => "i32",
            PixelType::F32 => "f32",
            _ => unreachable!(),
@@ -181,6 +182,30 @@ impl PixelExt for U16x3 {
    }
 }

+impl PixelExt for U16x4 {
+    fn load_big_image() -> DynamicImage {
+        ImageReader::open("./data/nasa-4928x3279-rgba.png")
+            .unwrap()
+            .decode()
+            .unwrap()
+    }
+
+    fn load_small_image() -> DynamicImage {
+        ImageReader::open("./data/nasa-852x567-rgba.png")
+            .unwrap()
+            .decode()
+            .unwrap()
+    }
+
+    fn img_into_bytes(img: DynamicImage) -> Vec<u8> {
+        img.to_rgba16()
+            .as_raw()
+            .iter()
+            .flat_map(|&c| c.to_le_bytes())
+            .collect()
+    }
+}
+
 impl PixelExt for I32 {
    fn img_into_bytes(img: DynamicImage) -> Vec<u8> {
        img.to_luma16()
@@ -217,6 +242,7 @@ pub fn save_result(image: &Image, name: &str) {
        PixelType::U16 => ColorType::L16,
        PixelType::U16x2 => ColorType::La16,
        PixelType::U16x3 => ColorType::Rgb16,
+        PixelType::U16x4 => ColorType::Rgba16,
        _ => panic!("Unsupported type of pixels"),
    };
    image::save_buffer(
@@ -1,6 +1,6 @@
 use std::num::NonZeroU32;

-use fast_image_resize::pixels::{Pixel, U16x2, U8x2, U8x4};
+use fast_image_resize::pixels::{Pixel, U16x2, U16x4, U8x2, U8x4};
 use fast_image_resize::{
    CpuExtensions, Image, ImageRows, ImageRowsMut, ImageView, ImageViewMut, MulDiv, PixelType,
 };
@@ -44,6 +44,16 @@ impl IntoImageRows for U16x2 {
    }
 }

+impl IntoImageRows for U16x4 {
+    fn into_image_rows(rows: Vec<&[Self]>) -> ImageRows {
+        ImageRows::U16x4(rows)
+    }
+
+    fn into_image_rows_mut(rows: Vec<&mut [Self]>) -> ImageRowsMut {
+        ImageRowsMut::U16x4(rows)
+    }
+}
+
 const fn p2(l: u8, a: u8) -> U8x2 {
    U8x2(u16::from_le_bytes([l, a]))
 }
@@ -272,6 +282,46 @@ mod multiply_alpha_u16x2 {
    }
 }

+#[cfg(test)]
+mod multiply_alpha_u16x4 {
+    use super::*;
+    use fast_image_resize::pixels::U16x4;
+
+    const SRC_PIXELS: [U16x4; 3] = [
+        U16x4([0xffff, 0x8000, 0, 0x8000]),
+        U16x4([0xffff, 0x8000, 0, 0xffff]),
+        U16x4([0xffff, 0x8000, 0, 0]),
+    ];
+    const RES_PIXELS: [U16x4; 3] = [
+        U16x4([0x8000, 0x4000, 0, 0x8000]),
+        U16x4([0xffff, 0x8000, 0, 0xffff]),
+        U16x4([0, 0, 0, 0]),
+    ];
+
+    #[cfg(target_arch = "x86_64")]
+    #[test]
+    fn avx2_test() {
+        for (s, r) in SRC_PIXELS.into_iter().zip(RES_PIXELS) {
+            mul_div_alpha_test(Oper::Mul, s, r, CpuExtensions::Avx2);
+        }
+    }
+
+    #[cfg(target_arch = "x86_64")]
+    #[test]
+    fn sse4_test() {
+        for (s, r) in SRC_PIXELS.into_iter().zip(RES_PIXELS) {
+            mul_div_alpha_test(Oper::Mul, s, r, CpuExtensions::Sse4_1);
+        }
+    }
+
+    #[test]
+    fn native_test() {
+        for (s, r) in SRC_PIXELS.into_iter().zip(RES_PIXELS) {
+            mul_div_alpha_test(Oper::Mul, s, r, CpuExtensions::None);
+        }
+    }
+}
+
 // Divides by alpha

 #[cfg(test)]
@@ -425,6 +475,51 @@ mod divide_alpha_u16x2 {
    }
 }

+#[cfg(test)]
+mod divide_alpha_u16x4 {
+    use super::*;
+
+    const OPER: Oper = Oper::Div;
+    const SRC_PIXELS: [U16x4; 3] = [
+        U16x4([0x8000, 0x4000, 0, 0x8000]),
+        U16x4([0xffff, 0x8000, 0, 0xffff]),
+        U16x4([0xffff, 0x8000, 0, 0]),
+    ];
+    const RES_PIXELS: [U16x4; 3] = [
+        U16x4([0xffff, 0x7fff, 0, 0x8000]),
+        U16x4([0xffff, 0x8000, 0, 0xffff]),
+        U16x4([0, 0, 0, 0]),
+    ];
+    const SIMD_RES_PIXELS: [U16x4; 3] = [
+        U16x4([0xffff, 0x8000, 0, 0x8000]),
+        U16x4([0xffff, 0x8000, 0, 0xffff]),
+        U16x4([0, 0, 0, 0]),
+    ];
+
+    #[cfg(target_arch = "x86_64")]
+    #[test]
+    fn avx2_test() {
+        for (s, r) in SRC_PIXELS.into_iter().zip(SIMD_RES_PIXELS) {
+            mul_div_alpha_test(OPER, s, r, CpuExtensions::Avx2);
+        }
+    }
+
+    #[cfg(target_arch = "x86_64")]
+    #[test]
+    fn sse4_test() {
+        for (s, r) in SRC_PIXELS.into_iter().zip(SIMD_RES_PIXELS) {
+            mul_div_alpha_test(OPER, s, r, CpuExtensions::Sse4_1);
+        }
+    }
+
+    #[test]
+    fn native_test() {
+        for (s, r) in SRC_PIXELS.into_iter().zip(RES_PIXELS) {
+            mul_div_alpha_test(OPER, s, r, CpuExtensions::None);
+        }
+    }
+}
+
 #[test]
 fn multiply_alpha_real_image_test() {
    let mut pixels = vec![0u8; 256 * 256 * 4];
@@ -449,3 +449,53 @@ fn upscale_u16x3() {
        );
    }
 }
+
+#[test]
+fn downscale_u16x4() {
+    type P = U16x4;
+    let buffer = downscale_test::<P>(ResizeAlg::Nearest, CpuExtensions::None);
+    assert_eq!(
+        testing::image_u16_checksum::<4>(&buffer),
+        [755050580, 756962660, 740848503, 1573303114]
+    );
+
+    let mut cpu_extensions_vec = vec![CpuExtensions::None];
+    #[cfg(target_arch = "x86_64")]
+    {
+        cpu_extensions_vec.push(CpuExtensions::Sse4_1);
+        cpu_extensions_vec.push(CpuExtensions::Avx2);
+    }
+    for cpu_extensions in cpu_extensions_vec {
+        let buffer =
+            downscale_test::<P>(ResizeAlg::Convolution(FilterType::Lanczos3), cpu_extensions);
+        assert_eq!(
+            testing::image_u16_checksum::<4>(&buffer),
+            [756269847, 757632467, 741478612, 1573563971]
+        );
+    }
+}
+
+#[test]
+fn upscale_u16x4() {
+    type P = U16x4;
+    let buffer = upscale_test::<P>(ResizeAlg::Nearest, CpuExtensions::None);
+    assert_eq!(
+        testing::image_u16_checksum::<4>(&buffer),
+        [296859917949, 296229709231, 288684470903, 607778112660]
+    );
+
+    let mut cpu_extensions_vec = vec![CpuExtensions::None];
+    #[cfg(target_arch = "x86_64")]
+    {
+        cpu_extensions_vec.push(CpuExtensions::Sse4_1);
+        cpu_extensions_vec.push(CpuExtensions::Avx2);
+    }
+    for cpu_extensions in cpu_extensions_vec {
+        let buffer =
+            upscale_test::<P>(ResizeAlg::Convolution(FilterType::Lanczos3), cpu_extensions);
+        assert_eq!(
+            testing::image_u16_checksum::<4>(&buffer),
+            [296888688348, 296243667797, 288698172180, 607776760273],
+        );
+    }
+}