Hardware Acceleration: the Same Algorithm, Ten Times Faster

#ifndef COMMON_H
#define COMMON_H

#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

/* Test image: 1920 x 1080, three bytes per pixel. */
#define WIDTH  1920
#define HEIGHT 1080
#define PIXELS ((size_t)(WIDTH) * (HEIGHT))

/* Monotonic clock: does not jump backward when the system clock is synced. */
static inline double now(void) {
    struct timespec t;
    clock_gettime(CLOCK_MONOTONIC, &t);
    return t.tv_sec + t.tv_nsec * 1e-9;
}

/* Deterministic image: the same data on every run, hence comparable. */
static void fill(uint8_t* rgb, size_t pixels) {
    for (size_t i = 0; i < pixels; i++) {
        rgb[i * 3 + 0] = (uint8_t)((i * 7) & 0xFF);
        rgb[i * 3 + 1] = (uint8_t)((i * 13) & 0xFF);
        rgb[i * 3 + 2] = (uint8_t)((i * 29) & 0xFF);
    }
}

/* Correctness check: a fast but wrong variant is worth nothing. */
static int compare(const uint8_t* a, const uint8_t* b, size_t n, const char* name) {
    for (size_t i = 0; i < n; i++) {
        if (a[i] != b[i]) {
            printf("  ERROR in %s at pixel %zu: %u != %u\n", name, i, a[i], b[i]);
            return 0;
        }
    }
    return 1;
}

/* Measurement with warm-up and repetitions: returns the MEDIAN time, not the mean.
   The median is not skewed by a single run disturbed by the system. */
static int cmp_double(const void* x, const void* y) {
    double a = *(const double*)x, b = *(const double*)y;
    return (a > b) - (a < b);
}

#define REPETITIONS 11

static double measure(void (*function)(const uint8_t*, uint8_t*, size_t),
                      const uint8_t* input, uint8_t* output, size_t pixels) {
    /* Warm-up: brings the data into cache and stabilizes the frequency. */
    for (int i = 0; i < 3; i++) function(input, output, pixels);

    double times[REPETITIONS];
    for (int i = 0; i < REPETITIONS; i++) {
        double t0 = now();
        function(input, output, pixels);
        times[i] = now() - t0;
    }
    qsort(times, REPETITIONS, sizeof(double), cmp_double);
    return times[REPETITIONS / 2];
}

#endif

#include "common.h"

void gray_scalar(const uint8_t* rgb, uint8_t* gray, size_t pixels) {
    for (size_t i = 0; i < pixels; i++) {
        uint32_t r = rgb[i * 3 + 0];
        uint32_t g = rgb[i * 3 + 1];
        uint32_t b = rgb[i * 3 + 2];
        gray[i] = (uint8_t)((77 * r + 150 * g + 29 * b) >> 8);
    }
}

#include "common.h"
#include <arm_neon.h>       /* the NEON intrinsic functions */

void gray_neon(const uint8_t* rgb, uint8_t* gray, size_t pixels) {
    /* The weights, replicated across all 8 lanes of the register. */
    const uint8x8_t w_r = vdup_n_u8(77);
    const uint8x8_t w_g = vdup_n_u8(150);
    const uint8x8_t w_b = vdup_n_u8(29);

    size_t i = 0;
    for (; i + 8 <= pixels; i += 8) {
        /* vld3_u8 loads 24 bytes and SPLITS them automatically across the
           three channels: 8 R values, 8 G values, 8 B values. Exactly what
           we need for interleaved data - an instruction a compiler
           rarely picks on its own. */
        uint8x8x3_t px = vld3_u8(rgb + i * 3);

        /* Widening multiply: 8 bits x 8 bits -> 16 bits, so nothing overflows.
           77 * 255 = 19635, so 8 bits would not have been enough. */
        uint16x8_t sum = vmull_u8(px.val[0], w_r);

        /* Widening multiply AND accumulate, in a single instruction. */
        sum = vmlal_u8(sum, px.val[1], w_g);
        sum = vmlal_u8(sum, px.val[2], w_b);

        /* Shift right by 8 AND narrow back to 8 bits,
           also in a single instruction. */
        vst1_u8(gray + i, vshrn_n_u16(sum, 8));
    }

    /* Tail: the remaining pixels, when their count is not a multiple of 8.
       This part is very easy to overlook and produces a wrong strip in the image. */
    for (; i < pixels; i++) {
        uint32_t r = rgb[i * 3 + 0], g = rgb[i * 3 + 1], b = rgb[i * 3 + 2];
        gray[i] = (uint8_t)((77 * r + 150 * g + 29 * b) >> 8);
    }
}

#include "common.h"
#include <omp.h>

void gray_neon(const uint8_t* rgb, uint8_t* gray, size_t pixels);   /* variant 3 */

void gray_parallel(const uint8_t* rgb, uint8_t* gray, size_t pixels) {
    int threads = omp_get_max_threads();

    /* Split into chunks that are a multiple of 8, so each thread can use
       the full NEON path, without its own tail. */
    size_t chunk = ((pixels / threads) / 8) * 8;

    #pragma omp parallel for schedule(static)
    for (int t = 0; t < threads; t++) {
        size_t start = t * chunk;
        size_t count = (t == threads - 1) ? (pixels - start) : chunk;
        gray_neon(rgb + start * 3, gray + start, count);
    }
}

#include "common.h"

void gray_scalar(const uint8_t*, uint8_t*, size_t);
void gray_neon(const uint8_t*, uint8_t*, size_t);
void gray_parallel(const uint8_t*, uint8_t*, size_t);

typedef struct {
    const char* name;
    void (*function)(const uint8_t*, uint8_t*, size_t);
} Variant;

int main(void) {
    uint8_t* rgb = malloc(PIXELS * 3);
    uint8_t* output = malloc(PIXELS);
    uint8_t* reference = malloc(PIXELS);
    if (!rgb || !output || !reference) return 1;

    fill(rgb, PIXELS);

    /* Correctness reference, computed once. */
    gray_scalar(rgb, reference, PIXELS);

    Variant variants[] = {
        { "scalar",         gray_scalar   },
        { "NEON",           gray_neon     },
        { "NEON + 4 threads", gray_parallel },
    };

    printf("image %d x %d = %.1f million pixels\n\n",
           WIDTH, HEIGHT, PIXELS / 1e6);
    printf("%-16s %10s %12s %10s %8s\n",
           "variant", "time", "Mpixels/s", "speedup", "correct");
    printf("%s\n", "------------------------------------------------------------");

    double t_reference = 0;
    for (size_t v = 0; v < sizeof(variants) / sizeof(variants[0]); v++) {
        memset(output, 0, PIXELS);
        double t = measure(variants[v].function, rgb, output, PIXELS);
        if (v == 0) t_reference = t;

        int correct = compare(output, reference, PIXELS, variants[v].name);

        printf("%-16s %8.2f ms %12.1f %9.2fx %8s\n",
               variants[v].name, t * 1000.0, PIXELS / t / 1e6,
               t_reference / t, correct ? "yes" : "NO");
    }

    /* How many frames per second does this mean, in practice? */
    printf("\nat 30 frames per second, the budget is %.1f ms per frame\n",
           1000.0 / 30.0);

    free(rgb); free(output); free(reference);
    return 0;
}