Compiler Explorer

Source code

fn filter_paeth_autovec<const N: usize>(a: [u8; N], b: [u8; N], c: [u8; N]) -> [u8; N] {
    // This is an optimized version of the paeth filter from the PNG specification, proposed by
    // Luca Versari for [FPNGE](https://www.lucaversari.it/FJXL_and_FPNGE.pdf). It operates
    // entirely on unsigned 8-bit quantities, making it more conducive to vectorization.
    //
    //     p = a + b - c
    //     pa = |p - a| = |a + b - c - a| = |b - c| = max(b, c) - min(b, c)
    //     pb = |p - b| = |a + b - c - b| = |a - c| = max(a, c) - min(a, c)
    //     pc = |p - c| = |a + b - c - c| = |(b - c) + (a - c)| = ...
    //
    // Further optimizing the calculation of `pc` a bit tricker. However, notice that:
    //
    //        a > c && b > c
    //    ==> (a - c) > 0 && (b - c) > 0
    //    ==> pc > (a - c) && pc > (b - c)
    //    ==> pc > |a - c| && pc > |b - c|
    //    ==> pc > pb && pc > pa
    //
    // Meaning that if `c` is smaller than `a` and `b`, the value of `pc` is irrelevant. Similar
    // reasoning applies if `c` is larger than the other two inputs. Assuming that `c >= b` and
    // `c <= b` or vice versa:
    //
    //     pc = ||b - c| - |a - c|| =  |pa - pb| = max(pa, pb) - min(pa, pb)
    
    
    let pa = simd_sub(simd_max(b,c), simd_min(c,b));
    let pb = simd_sub(simd_max(a,c), simd_min(c,a));

let mut pc = simd_sub(simd_max(pa, pb), simd_min(pa, pb));
    for (i, o) in pc.iter_mut().enumerate() {
        if (a[i] < c[i]) != (c[i] < b[i]) {
            *o = 255;
        }
    }

let mut crit_1 = [false; N];
    for (i, o) in crit_1.iter_mut().enumerate() {
        if pa[i] <= pb[i] && pa[i] <= pc[i] {
            *o = true;
        }
    }

let mut crit_2 = [false; N];
    for (i, o) in crit_1.iter_mut().enumerate() {
        if pb[i] <= pc[i] {
            *o = true;
        }
    }

let mut out = [0; N];
    for (i, o) in out.iter_mut().enumerate() {
        if crit_1[i] {
            *o = a[i];
        } else if crit_2[i] {
            *o = b[i];
        } else {
            *o = c[i];
        }
    }

out
}

pub fn instantiate(a: [u8; 16], b: [u8; 16], c: [u8; 16]) -> [u8; 16] {
    filter_paeth_autovec(a, b, c)
}

#[inline]
fn simd_max<const N: usize>(a: [u8; N], b: [u8; N]) -> [u8; N] {
    let mut out = [0; N];
    for (i, o) in out.iter_mut().enumerate() {
        *o = a[i].max(b[i])
    }
    out
}

#[inline]
fn simd_min<const N: usize>(a: [u8; N], b: [u8; N]) -> [u8; N] {
    let mut out = [0; N];
    for (i, o) in out.iter_mut().enumerate() {
        *o = a[i].min(b[i])
    }
    out
}

#[inline]
fn simd_sub<const N: usize>(a: [u8; N], b: [u8; N]) -> [u8; N] {
    let mut out = [0; N];
    for (i, o) in out.iter_mut().enumerate() {
        *o = a[i] - b[i]
    }
    out
}