lavfi/nlmeans: add AArch64 SIMD for compute_safe_ssd_integral_image

ssd_integral_image_c: 49204.6
ssd_integral_image_neon: 28346.8

6 files changed, 175 insertions, 8 deletions

diff --git a/libavfilter/aarch64/Makefile b/libavfilter/aarch64/Makefile
new file mode 100644
index 0000000..b58daa3
--- /dev/null
+++ b/libavfilter/aarch64/Makefile
@@ -0,0 +1,3 @@
+OBJS-$(CONFIG_NLMEANS_FILTER)                += aarch64/vf_nlmeans_init.o
+
+NEON-OBJS-$(CONFIG_NLMEANS_FILTER)           += aarch64/vf_nlmeans_neon.o
diff --git a/libavfilter/aarch64/vf_nlmeans_init.c b/libavfilter/aarch64/vf_nlmeans_init.c
new file mode 100644
index 0000000..a1edefb
--- /dev/null
+++ b/libavfilter/aarch64/vf_nlmeans_init.c
@@ -0,0 +1,33 @@
+/*
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "libavutil/aarch64/cpu.h"
+#include "libavfilter/vf_nlmeans.h"
+
+void ff_compute_safe_ssd_integral_image_neon(uint32_t *dst, ptrdiff_t dst_linesize_32,
+                                             const uint8_t *s1, ptrdiff_t linesize1,
+                                             const uint8_t *s2, ptrdiff_t linesize2,
+                                             int w, int h);
+
+av_cold void ff_nlmeans_init_aarch64(NLMeansDSPContext *dsp)
+{
+    int cpu_flags = av_get_cpu_flags();
+
+    if (have_neon(cpu_flags))
+        dsp->compute_safe_ssd_integral_image = ff_compute_safe_ssd_integral_image_neon;
+}
diff --git a/libavfilter/aarch64/vf_nlmeans_neon.S b/libavfilter/aarch64/vf_nlmeans_neon.S
new file mode 100644
index 0000000..6308a42
--- /dev/null
+++ b/libavfilter/aarch64/vf_nlmeans_neon.S
@@ -0,0 +1,80 @@
+/*
+ * Copyright (c) 2018 Clément Bœsch <u pkh me>
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "libavutil/aarch64/asm.S"
+
+// acc_sum_store(ABCD) = {X+A, X+A+B, X+A+B+C, X+A+B+C+D}
+.macro acc_sum_store x, xb
+        dup             v24.4S, v24.4S[3]                               // ...X -> XXXX
+        ext             v25.16B, v26.16B, \xb, #12                      // ext(0000,ABCD,12)=0ABC
+        add             v24.4S, v24.4S, \x                              // XXXX+ABCD={X+A,X+B,X+C,X+D}
+        add             v24.4S, v24.4S, v25.4S                          // {X+A,X+B+A,X+C+B,X+D+C}       (+0ABC)
+        ext             v25.16B, v26.16B, v25.16B, #12                  // ext(0000,0ABC,12)=00AB
+        add             v24.4S, v24.4S, v25.4S                          // {X+A,X+B+A,X+C+B+A,X+D+C+B}   (+00AB)
+        ext             v25.16B, v26.16B, v25.16B, #12                  // ext(0000,00AB,12)=000A
+        add             v24.4S, v24.4S, v25.4S                          // {X+A,X+B+A,X+C+B+A,X+D+C+B+A} (+000A)
+        st1             {v24.4S}, [x0], #16                             // write 4x32-bit final values
+.endm
+
+function ff_compute_safe_ssd_integral_image_neon, export=1
+        movi            v26.4S, #0                                      // used as zero for the "rotations" in acc_sum_store
+        sub             x3, x3, w6, UXTW                                // s1 padding (s1_linesize - w)
+        sub             x5, x5, w6, UXTW                                // s2 padding (s2_linesize - w)
+        sub             x9, x0, x1, UXTW #2                             // dst_top
+        sub             x1, x1, w6, UXTW                                // dst padding (dst_linesize_32 - w)
+        lsl             x1, x1, #2                                      // dst padding expressed in bytes
+1:      mov             w10, w6                                         // width copy for each line
+        sub             x0, x0, #16                                     // beginning of the dst line minus 4 sums
+        sub             x8, x9, #4                                      // dst_top-1
+        ld1             {v24.4S}, [x0], #16                             // load ...X (contextual last sums)
+2:      ld1             {v0.16B}, [x2], #16                             // s1[x + 0..15]
+        ld1             {v1.16B}, [x4], #16                             // s2[x + 0..15]
+        ld1             {v16.4S,v17.4S}, [x8], #32                      // dst_top[x + 0..7 - 1]
+        usubl           v2.8H, v0.8B,  v1.8B                            // d[x + 0..7]  = s1[x + 0..7]  - s2[x + 0..7]
+        usubl2          v3.8H, v0.16B, v1.16B                           // d[x + 8..15] = s1[x + 8..15] - s2[x + 8..15]
+        ld1             {v18.4S,v19.4S}, [x8], #32                      // dst_top[x + 8..15 - 1]
+        smull           v4.4S, v2.4H, v2.4H                             // d[x + 0..3]^2
+        smull2          v5.4S, v2.8H, v2.8H                             // d[x + 4..7]^2
+        ld1             {v20.4S,v21.4S}, [x9], #32                      // dst_top[x + 0..7]
+        smull           v6.4S, v3.4H, v3.4H                             // d[x + 8..11]^2
+        smull2          v7.4S, v3.8H, v3.8H                             // d[x + 12..15]^2
+        ld1             {v22.4S,v23.4S}, [x9], #32                      // dst_top[x + 8..15]
+        sub             v0.4S, v20.4S, v16.4S                           // dst_top[x + 0..3] - dst_top[x + 0..3 - 1]
+        sub             v1.4S, v21.4S, v17.4S                           // dst_top[x + 4..7] - dst_top[x + 4..7 - 1]
+        add             v0.4S, v0.4S, v4.4S                             // + d[x + 0..3]^2
+        add             v1.4S, v1.4S, v5.4S                             // + d[x + 4..7]^2
+        sub             v2.4S, v22.4S, v18.4S                           // dst_top[x +  8..11] - dst_top[x +  8..11 - 1]
+        sub             v3.4S, v23.4S, v19.4S                           // dst_top[x + 12..15] - dst_top[x + 12..15 - 1]
+        add             v2.4S, v2.4S, v6.4S                             // + d[x +  8..11]^2
+        add             v3.4S, v3.4S, v7.4S                             // + d[x + 12..15]^2
+        acc_sum_store   v0.4S, v0.16B                                   // accumulate and store dst[ 0..3]
+        acc_sum_store   v1.4S, v1.16B                                   // accumulate and store dst[ 4..7]
+        acc_sum_store   v2.4S, v2.16B                                   // accumulate and store dst[ 8..11]
+        acc_sum_store   v3.4S, v3.16B                                   // accumulate and store dst[12..15]
+        subs            w10, w10, #16                                   // width dec
+        b.ne            2b                                              // loop til next line
+        add             x2, x2, x3                                      // skip to next line (s1)
+        add             x4, x4, x5                                      // skip to next line (s2)
+        add             x0, x0, x1                                      // skip to next line (dst)
+        add             x9, x9, x1                                      // skip to next line (dst_top)
+        subs            w7, w7, #1                                      // height dec
+        b.ne            1b
+        ret
+endfunc
diff --git a/libavfilter/tests/integral.c b/libavfilter/tests/integral.c
index 049fefa..2a8e8ff 100644
--- a/libavfilter/tests/integral.c
+++ b/libavfilter/tests/integral.c
@@ -57,6 +57,10 @@ int main(void)
     uint32_t *ii_start  = ii  + ii_lz_32 + 1; // skip top 0-line and left 0-column
     uint32_t *ii_start2 = ii2 + ii_lz_32 + 1; // skip top 0-line and left 0-column
 
+    NLMeansDSPContext dsp = {0};
+
+    ff_nlmeans_init(&dsp);
+
     if (!ii || !ii2)
         return -1;
 
@@ -64,7 +68,7 @@ int main(void)
         for (xoff = -e; xoff <= e; xoff++) {
             printf("xoff=%d yoff=%d\n", xoff, yoff);
 
-            compute_ssd_integral_image(ii_start, ii_lz_32,
+            compute_ssd_integral_image(&dsp, ii_start, ii_lz_32,
                                        src, lz, xoff, yoff, e, w, h);
             display_integral(ii_start, ii_w, ii_h, ii_lz_32);
 
diff --git a/libavfilter/vf_nlmeans.c b/libavfilter/vf_nlmeans.c
index b081a4e..c30e444 100644
--- a/libavfilter/vf_nlmeans.c
+++ b/libavfilter/vf_nlmeans.c
@@ -20,7 +20,6 @@
 
 /**
  * @todo
- * - SIMD for compute_safe_ssd_integral_image
  * - SIMD for final weighted averaging
  * - better automatic defaults? see "Parameters" @ http://www.ipol.im/pub/art/2011/bcm_nlm/
  * - temporal support (probably doesn't need any displacement according to
@@ -37,6 +36,7 @@
 #include "avfilter.h"
 #include "formats.h"
 #include "internal.h"
+#include "vf_nlmeans.h"
 #include "video.h"
 
 struct weighted_avg {
@@ -66,6 +66,7 @@ typedef struct NLMeansContext {
     double weight_lut[WEIGHT_LUT_SIZE];         // lookup table mapping (scaled) patch differences to their associated weights
     double pdiff_lut_scale;                     // scale factor for patch differences before looking into the LUT
     int max_meaningful_diff;                    // maximum difference considered (if the patch difference is too high we ignore the pixel)
+    NLMeansDSPContext dsp;
 } NLMeansContext;
 
 #define OFFSET(x) offsetof(NLMeansContext, x)
@@ -240,7 +241,8 @@ static inline void compute_unsafe_ssd_integral_image(uint32_t *dst, ptrdiff_t ds
  * @param h                 source height
  * @param e                 research padding edge
  */
-static void compute_ssd_integral_image(uint32_t *ii, ptrdiff_t ii_linesize_32,
+static void compute_ssd_integral_image(const NLMeansDSPContext *dsp,
+                                       uint32_t *ii, ptrdiff_t ii_linesize_32,
                                        const uint8_t *src, ptrdiff_t linesize, int offx, int offy,
                                        int e, int w, int h)
 {
@@ -291,10 +293,10 @@ static void compute_ssd_integral_image(uint32_t *ii, ptrdiff_t ii_linesize_32,
     av_assert1(startx_safe - s2x >= 0); av_assert1(startx_safe - s2x < w);
     av_assert1(starty_safe - s2y >= 0); av_assert1(starty_safe - s2y < h);
     if (safe_pw && safe_ph)
-        compute_safe_ssd_integral_image_c(ii + starty_safe*ii_linesize_32 + startx_safe, ii_linesize_32,
-                                          src + (starty_safe - s1y) * linesize + (startx_safe - s1x), linesize,
-                                          src + (starty_safe - s2y) * linesize + (startx_safe - s2x), linesize,
-                                          safe_pw, safe_ph);
+        dsp->compute_safe_ssd_integral_image(ii + starty_safe*ii_linesize_32 + startx_safe, ii_linesize_32,
+                                             src + (starty_safe - s1y) * linesize + (startx_safe - s1x), linesize,
+                                             src + (starty_safe - s2y) * linesize + (startx_safe - s2x), linesize,
+                                             safe_pw, safe_ph);
 
     // right part of the integral
     compute_unsafe_ssd_integral_image(ii, ii_linesize_32,
@@ -431,7 +433,7 @@ static int nlmeans_plane(AVFilterContext *ctx, int w, int h, int p, int r,
                     .p            = p,
                 };
 
-                compute_ssd_integral_image(s->ii, s->ii_lz_32,
+                compute_ssd_integral_image(&s->dsp, s->ii, s->ii_lz_32,
                                            src, src_linesize,
                                            offx, offy, e, w, h);
                 ctx->internal->execute(ctx, nlmeans_slice, &td, NULL,
@@ -489,6 +491,14 @@ static int filter_frame(AVFilterLink *inlink, AVFrame *in)
     }                                                           \
 } while (0)
 
+void ff_nlmeans_init(NLMeansDSPContext *dsp)
+{
+    dsp->compute_safe_ssd_integral_image = compute_safe_ssd_integral_image_c;
+
+    if (ARCH_AARCH64)
+        ff_nlmeans_init_aarch64(dsp);
+}
+
 static av_cold int init(AVFilterContext *ctx)
 {
     int i;
@@ -520,6 +530,8 @@ static av_cold int init(AVFilterContext *ctx)
            s->research_size, s->research_size, s->research_size_uv, s->research_size_uv,
            s->patch_size,    s->patch_size,    s->patch_size_uv,    s->patch_size_uv);
 
+    ff_nlmeans_init(&s->dsp);
+
     return 0;
 }
 
diff --git a/libavfilter/vf_nlmeans.h b/libavfilter/vf_nlmeans.h
new file mode 100644
index 0000000..0a9aab2
--- /dev/null
+++ b/libavfilter/vf_nlmeans.h
@@ -0,0 +1,35 @@
+/*
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifndef AVFILTER_NLMEANS_H
+#define AVFILTER_NLMEANS_H
+
+#include <stddef.h>
+#include <stdint.h>
+
+typedef struct NLMeansDSPContext {
+    void (*compute_safe_ssd_integral_image)(uint32_t *dst, ptrdiff_t dst_linesize_32,
+                                            const uint8_t *s1, ptrdiff_t linesize1,
+                                            const uint8_t *s2, ptrdiff_t linesize2,
+                                            int w, int h);
+} NLMeansDSPContext;
+
+void ff_nlmeans_init(NLMeansDSPContext *dsp);
+void ff_nlmeans_init_aarch64(NLMeansDSPContext *dsp);
+
+#endif /* AVFILTER_NLMEANS_H */