Speed up pixel compositing with SIMD vectorization

Alpha blending of graphics protocol images and animation frames is
2-3.5x faster and glyph alpha masks are composited onto canvases using
the same vectorized primitives

docs/changelog.rst

@@ -37,6 +37,9 @@ should see a 15-35% improvement depending on workload. Some details:
 
 #. Improve throughput when processing large amounts of plain text by another ~10% by finding runs of printable ASCII chars with :term:`SIMD`, filling cells using wide stores and skipping unnecessary work in the scrolling hot path when there are no images
 
+#. Speed up pixel compositing with :term:`SIMD` vectorization: alpha blending of graphics protocol images and animation frames is 2-3.5x faster and glyph alpha masks are composited onto canvases using the same vectorized primitives.
+
+
 
 
 Vertical tabs [0.48]

kitty/fonts.c

@@ -14,6 +14,7 @@
 #include "decorations.h"
 #include "glyph-cache.h"
 #include "print-graphics.h"
+#include "simd-string.h"
 
 #define MISSING_GLYPH 1
 #define MAX_NUM_EXTRA_GLYPHS_PUA 4u
@@ -870,15 +871,12 @@ static PyObject *descriptor_for_idx = NULL;
 void
 render_alpha_mask(
     const uint8_t *alpha_mask, pixel *dest, const Region *src_rect, const Region *dest_rect, size_t src_stride, size_t dest_stride, pixel color_rgb) {
-    pixel col = (color_rgb << 8) & 0xffffff00;
+    const size_t src_width = src_rect->right > src_rect->left ? src_rect->right - src_rect->left : 0;
+    const size_t dest_width = dest_rect->right > dest_rect->left ? dest_rect->right - dest_rect->left : 0;
+    const size_t width = MIN(src_width, dest_width);
+    if (!width) return;
     for (size_t sr = src_rect->top, dr = dest_rect->top; sr < src_rect->bottom && dr < dest_rect->bottom; sr++, dr++) {
-        pixel *d = dest + dest_stride * dr;
-        const uint8_t *s = alpha_mask + src_stride * sr;
-        for (size_t sc = src_rect->left, dc = dest_rect->left; sc < src_rect->right && dc < dest_rect->right; sc++, dc++) {
-            uint8_t src_alpha = d[dc] & 0xff;
-            uint8_t alpha = s[sc];
-            d[dc] = col | MAX(alpha, src_alpha);
-        }
+        composite_alpha_mask(dest + dest_stride * dr + dest_rect->left, alpha_mask + src_stride * sr + src_rect->left, width, color_rgb);
     }
 }
 
@@ -1055,9 +1053,8 @@ render_scaled_decoration(FontCellMetrics unscaled_metrics, FontCellMetrics scale
     unsigned src_limit = MIN(scaled_metrics.cell_height, src.bottom), dest_limit = MIN(unscaled_metrics.cell_height, dest.bottom);
     unsigned cell_width = MIN(scaled_metrics.cell_width, unscaled_metrics.cell_width);
     for (unsigned srcy = src.top, desty = dest.top; srcy < src_limit && desty < dest_limit; srcy++, desty++) {
-        uint8_t *srcp = alpha_mask + cell_width * srcy;
-        pixel *destp = output + cell_width * desty;
-        for (unsigned x = 0; x < cell_width; x++) destp[x] = 0xffffff00 | srcp[x];
+        // the output was zeroed above so this sets destp[x] = 0xffffff00 | srcp[x]
+        composite_alpha_mask(output + cell_width * desty, alpha_mask + cell_width * srcy, cell_width, 0xffffff);
     }
 }
 
@@ -2494,7 +2491,9 @@ concat_cells(PyObject UNUSED *self, PyObject *args) {
             void *s = ((uint8_t *)PyBytes_AS_STRING(PyTuple_GET_ITEM(cells, c)));
             if (is_32_bit) {
                 pixel *src = (pixel *)s + cell_width * r;
-                for (i = 0; i < cell_width; i++, dest++) dest[0] = alpha_blend(src[0], bgcolor);
+                for (i = 0; i < cell_width; i++) dest[i] = (pixel)bgcolor;
+                blend_over_opaque((uint8_t *)dest, 4, (const uint8_t *)src, cell_width);
+                dest += cell_width;
             } else {
                 uint8_t *src = (uint8_t *)s + cell_width * r;
                 for (i = 0; i < cell_width; i++, dest++) dest[0] = alpha_blend(0x00ffffff | ((src[i] & 0xff) << 24), bgcolor);

kitty/graphics.c

@@ -11,6 +11,7 @@
 #include "disk-cache.h"
 #include "iqsort.h"
 #include "safe-wrappers.h"
+#include "simd-string.h"
 
 #include <sys/types.h>
 #include <sys/stat.h>
@@ -1542,29 +1543,6 @@ typedef struct {
     bool is_4byte_aligned, is_opaque, transient;
 } CoalescedFrameData;
 
-static void
-blend_on_opaque(uint8_t *under_px, const uint8_t *over_px) {
-    const float alpha = (float)over_px[3] / 255.f;
-    const float alpha_op = 1.f - alpha;
-    for (unsigned i = 0; i < 3; i++) under_px[i] = (uint8_t)(over_px[i] * alpha + under_px[i] * alpha_op);
-}
-
-static void
-alpha_blend(uint8_t *dest_px, const uint8_t *src_px) {
-    if (src_px[3]) {
-        const float dest_a = (float)dest_px[3] / 255.f, src_a = (float)src_px[3] / 255.f;
-        const float alpha = src_a + dest_a * (1.f - src_a);
-        dest_px[3] = (uint8_t)(255 * alpha);
-        if (!dest_px[3]) {
-            dest_px[0] = 0;
-            dest_px[1] = 0;
-            dest_px[2] = 0;
-            return;
-        }
-        for (unsigned i = 0; i < 3; i++) dest_px[i] = (uint8_t)((src_px[i] * src_a + dest_px[i] * dest_a * (1.f - src_a)) / alpha);
-    }
-}
-
 typedef struct {
     bool needs_blending;
     uint32_t over_px_sz, under_px_sz;
@@ -1572,93 +1550,63 @@ typedef struct {
     uint32_t stride;
 } ComposeData;
 
-#define COPY_RGB              \
-    under_px[0] = over_px[0]; \
-    under_px[1] = over_px[1]; \
-    under_px[2] = over_px[2];
-#define COPY_PIXELS                                                                         \
-    if (d.needs_blending) {                                                                 \
-        if (d.under_px_sz == 3) { ROW_ITER PIX_ITER blend_on_opaque(under_px, over_px); }   \
-    }                                                                                       \
-    }                                                                                       \
-    else { ROW_ITER PIX_ITER alpha_blend(under_px, over_px); }                              \
-    }                                                                                       \
-    }                                                                                       \
-    }                                                                                       \
-    else {                                                                                  \
-        if (d.under_px_sz == 4) {                                                           \
-            if (d.over_px_sz == 4) { ROW_ITER PIX_ITER COPY_RGB under_px[3] = over_px[3]; } \
-        }                                                                                   \
-    }                                                                                       \
-    else { ROW_ITER PIX_ITER COPY_RGB under_px[3] = 255; }                                  \
-    }                                                                                       \
-    }                                                                                       \
-    }                                                                                       \
-    else { ROW_ITER PIX_ITER COPY_RGB }                                                     \
-    }                                                                                       \
-    }                                                                                       \
+static void
+copy_or_blend_row(
+    uint8_t *under_row, const uint8_t *over_row, const unsigned num_px, const uint32_t under_px_sz, const uint32_t over_px_sz, const bool needs_blending) {
+    if (needs_blending && over_px_sz == 4) {
+        // over is straight alpha RGBA
+        if (under_px_sz == 4) blend_over_straight(under_row, over_row, num_px);
+        else blend_over_opaque(under_row, 3, over_row, num_px);
+    } else {
+        for (unsigned x = 0; x < num_px; x++) {
+            uint8_t *under_px = under_row + (size_t)under_px_sz * x;
+            const uint8_t *over_px = over_row + (size_t)over_px_sz * x;
+            under_px[0] = over_px[0];
+            under_px[1] = over_px[1];
+            under_px[2] = over_px[2];
+            if (under_px_sz == 4) under_px[3] = over_px_sz == 4 ? over_px[3] : 255;
+        }
     }
-
+}
 
 static void
 compose_rectangles(const ComposeData d, uint8_t *under_data, const size_t under_data_sz, const uint8_t *over_data, const size_t over_data_sz) {
     // compose two equal sized, non-overlapping rectangles at different offsets.
-    // The per-row guard in ROW_ITER ensures that inconsistent geometry (offsets
-    // or dimensions larger than the actual buffers) can never cause an out of
+    // The per-row guard ensures that inconsistent geometry (offsets or
+    // dimensions larger than the actual buffers) can never cause an out of
     // bounds read or write. All offset arithmetic is done in size_t to avoid
     // 32-bit overflow.
     const bool can_copy_rows = !d.needs_blending && d.over_px_sz == d.under_px_sz;
     const unsigned min_width = MIN(d.under_width, d.over_width);
-#define ROW_ITER                                                                                                                                \
-    for (unsigned y = 0; y < d.under_height && y < d.over_height; y++) {                                                                        \
-        const size_t under_off = (size_t)(y + d.under_offset_y) * d.under_px_sz * d.stride + (size_t)d.under_offset_x * d.under_px_sz;          \
-        const size_t over_off = (size_t)(y + d.over_offset_y) * d.over_px_sz * d.stride + (size_t)d.over_offset_x * d.over_px_sz;               \
-        if (under_off + (size_t)d.under_px_sz * min_width > under_data_sz || over_off + (size_t)d.over_px_sz * min_width > over_data_sz) break; \
-        uint8_t *under_row = under_data + under_off;                                                                                            \
+    for (unsigned y = 0; y < d.under_height && y < d.over_height; y++) {
+        const size_t under_off = (size_t)(y + d.under_offset_y) * d.under_px_sz * d.stride + (size_t)d.under_offset_x * d.under_px_sz;
+        const size_t over_off = (size_t)(y + d.over_offset_y) * d.over_px_sz * d.stride + (size_t)d.over_offset_x * d.over_px_sz;
+        if (under_off + (size_t)d.under_px_sz * min_width > under_data_sz || over_off + (size_t)d.over_px_sz * min_width > over_data_sz) break;
+        uint8_t *under_row = under_data + under_off;
         const uint8_t *over_row = over_data + over_off;
-    if (can_copy_rows) { ROW_ITER memcpy(under_row, over_row, (size_t)d.over_px_sz * min_width); }
-    return;
-}
-#define PIX_ITER                                             \
-    for (unsigned x = 0; x < min_width; x++) {               \
-        uint8_t *under_px = under_row + (d.under_px_sz * x); \
-        const uint8_t *over_px = over_row + (d.over_px_sz * x);
-COPY_PIXELS
-#undef PIX_ITER
-#undef ROW_ITER
+        if (can_copy_rows) memcpy(under_row, over_row, (size_t)d.over_px_sz * min_width);
+        else copy_or_blend_row(under_row, over_row, min_width, d.under_px_sz, d.over_px_sz, d.needs_blending);
+    }
 }
 
 static void
 compose(const ComposeData d, uint8_t *under_data, const size_t under_data_sz, const uint8_t *over_data, const size_t over_data_sz) {
-    // The per-row guard in ROW_ITER ensures that inconsistent geometry (an
-    // overlay declaring more data than it actually has, or an oversized base)
-    // can never cause an out of bounds read or write. All offset arithmetic is
+    // The per-row guard ensures that inconsistent geometry (an overlay
+    // declaring more data than it actually has, or an oversized base) can
+    // never cause an out of bounds read or write. All offset arithmetic is
     // done in size_t to avoid 32-bit overflow.
     const bool can_copy_rows = !d.needs_blending && d.over_px_sz == d.under_px_sz;
     unsigned min_row_sz = d.over_offset_x < d.under_width ? d.under_width - d.over_offset_x : 0;
     min_row_sz = MIN(min_row_sz, d.over_width);
-#define ROW_ITER                                                                                                                                  \
-    for (unsigned y = 0; y + d.over_offset_y < d.under_height && y < d.over_height; y++) {                                                        \
-        const size_t under_off = (size_t)(y + d.over_offset_y) * d.under_px_sz * d.under_width + (size_t)d.under_px_sz * d.over_offset_x;         \
-        const size_t over_off = (size_t)y * d.over_px_sz * d.over_width;                                                                          \
-        if (under_off + (size_t)d.under_px_sz * min_row_sz > under_data_sz || over_off + (size_t)d.over_px_sz * min_row_sz > over_data_sz) break; \
-        uint8_t *under_row = under_data + under_off;                                                                                              \
+    for (unsigned y = 0; y + d.over_offset_y < d.under_height && y < d.over_height; y++) {
+        const size_t under_off = (size_t)(y + d.over_offset_y) * d.under_px_sz * d.under_width + (size_t)d.under_px_sz * d.over_offset_x;
+        const size_t over_off = (size_t)y * d.over_px_sz * d.over_width;
+        if (under_off + (size_t)d.under_px_sz * min_row_sz > under_data_sz || over_off + (size_t)d.over_px_sz * min_row_sz > over_data_sz) break;
+        uint8_t *under_row = under_data + under_off;
         const uint8_t *over_row = over_data + over_off;
-#define END_ITER }
-    if (can_copy_rows) {
-        ROW_ITER memcpy(under_row, over_row, (size_t)d.over_px_sz * min_row_sz);
-        END_ITER
-        return;
+        if (can_copy_rows) memcpy(under_row, over_row, (size_t)d.over_px_sz * min_row_sz);
+        else copy_or_blend_row(under_row, over_row, min_row_sz, d.under_px_sz, d.over_px_sz, d.needs_blending);
     }
-#define PIX_ITER                                             \
-    for (unsigned x = 0; x < min_row_sz; x++) {              \
-        uint8_t *under_px = under_row + (d.under_px_sz * x); \
-        const uint8_t *over_px = over_row + (d.over_px_sz * x);
-    COPY_PIXELS
-#undef COPY_RGB
-#undef PIX_ITER
-#undef ROW_ITER
-#undef END_ITER
 }
 
 static CoalescedFrameData

kitty/simd-string-512.c

@@ -5,7 +5,8 @@
  * Distributed under terms of the GPL3 license.
  *
  * AVX-512 (F, BW, VL, VBMI2) implementations of utf8_decode_to_esc,
- * find_either_of_two_bytes, xor_data64 and printable_ascii_run_length. Unlike
+ * find_either_of_two_bytes, xor_data64, printable_ascii_run_length and the
+ * pixel compositing functions. Unlike
  * the 128/256 bit implementations in simd-string-impl.h these are x86-64 only
  * and use native intrinsics, since the algorithms are built around mask
  * registers and masked loads/stores (and vpcompressb for UTF-8 decoding),
@@ -79,6 +80,93 @@ printable_ascii_run_length_512(const uint32_t *chars, const size_t sz) {
     return ans;
 }
 
+// Pixel compositing {{{
+
+void
+composite_alpha_mask_512(uint32_t *dst, const uint8_t *mask, const size_t num_pixels, const uint32_t color_rgb) {
+    const __m512i col = _mm512_set1_epi32((int32_t)((color_rgb << 8) & 0xffffff00)), low_byte = _mm512_set1_epi32(0xff);
+    for (size_t i = 0; i < num_pixels; i += 16) {
+        const __mmask16 k = num_pixels - i >= 16 ? (__mmask16)0xffff : (__mmask16)((1u << (num_pixels - i)) - 1);
+        const __m512i m = _mm512_cvtepu8_epi32(_mm_maskz_loadu_epi8(k, mask + i));
+        const __m512i d = _mm512_maskz_loadu_epi32(k, dst + i);
+        _mm512_mask_storeu_epi32(dst + i, k, _mm512_or_si512(col, _mm512_max_epu32(m, _mm512_and_si512(d, low_byte))));
+    }
+    _mm256_zeroupper();
+}
+
+static inline __m512i
+div255_epu16_512(const __m512i x) {
+    // rounding division of 16-bit lanes by 255, exact for values <= 65407, matches div255_round()
+    const __m512i y = _mm512_add_epi16(x, _mm512_set1_epi16(128));
+    return _mm512_srli_epi16(_mm512_add_epi16(y, _mm512_srli_epi16(y, 8)), 8);
+}
+
+// Blend the 4-byte RGBA pixels in over onto the pixels in under, with under considered fully
+// opaque: out_c = round((over_c * alpha + under_c * (255 - alpha)) / 255) for each of the first
+// three channels, with the alpha bytes set to 255
+static inline __m512i
+blend_opaque_pixels_512(const __m512i under, const __m512i over) {
+    const __m512i zero = _mm512_setzero_si512();
+    const __m512i alpha = _mm512_shuffle_epi8(over, _mm512_broadcast_i32x4(_mm_set_epi8(15, 15, 15, 15, 11, 11, 11, 11, 7, 7, 7, 7, 3, 3, 3, 3)));
+    const __m512i inv_alpha = _mm512_xor_si512(alpha, _mm512_set1_epi64(-1)); // 255 - alpha in every byte
+    const __m512i lo = div255_epu16_512(_mm512_add_epi16(
+        _mm512_mullo_epi16(_mm512_unpacklo_epi8(over, zero), _mm512_unpacklo_epi8(alpha, zero)),
+        _mm512_mullo_epi16(_mm512_unpacklo_epi8(under, zero), _mm512_unpacklo_epi8(inv_alpha, zero))));
+    const __m512i hi = div255_epu16_512(_mm512_add_epi16(
+        _mm512_mullo_epi16(_mm512_unpackhi_epi8(over, zero), _mm512_unpackhi_epi8(alpha, zero)),
+        _mm512_mullo_epi16(_mm512_unpackhi_epi8(under, zero), _mm512_unpackhi_epi8(inv_alpha, zero))));
+    // the per 128-bit lane unpacks and pack are symmetric so byte order is preserved
+    return _mm512_or_si512(_mm512_packus_epi16(lo, hi), _mm512_set1_epi32((int32_t)0xff000000));
+}
+
+void
+blend_over_opaque_512(uint8_t *dst, const unsigned dst_bpp, const uint8_t *src, const size_t num_pixels) {
+    if (dst_bpp == 4) {
+        for (size_t i = 0; i < num_pixels; i += 16) {
+            const size_t px = MIN(num_pixels - i, (size_t)16u);
+            const __mmask64 k = px == 16 ? ~0ull : ((1ull << (4 * px)) - 1);
+            const __m512i s = _mm512_maskz_loadu_epi8(k, src + 4 * i), d = _mm512_maskz_loadu_epi8(k, dst + 4 * i);
+            _mm512_mask_storeu_epi8(dst + 4 * i, k, blend_opaque_pixels_512(d, s));
+        }
+    } else {
+        // 3-byte destination pixels: the byte expand and compress instructions this would need
+        // (VBMI2) are slower than the 128-bit shuffle based kernel, at least on the Zen CPUs
+        // tested, so just use that.
+        blend_over_opaque_128(dst, dst_bpp, src, num_pixels);
+        return;
+    }
+    _mm256_zeroupper();
+}
+
+void
+blend_over_straight_512(uint8_t *dst, const uint8_t *src, const size_t num_pixels) {
+    // Four pixels per iteration unpacked into 32-bit lanes, four consecutive lanes per pixel.
+    // Matches the arithmetic of blend_pixel_over_straight() bit for bit: exact integer numerator
+    // and denominator, IEEE single precision division, round to nearest.
+    const __m512i c255 = _mm512_set1_epi32(255), c128 = _mm512_set1_epi32(128), zero = _mm512_setzero_si512();
+    for (size_t i = 0; i < num_pixels; i += 4) {
+        const size_t px = MIN(num_pixels - i, (size_t)4u);
+        const __mmask16 k = px == 4 ? (__mmask16)0xffff : (__mmask16)((1u << (4 * px)) - 1);
+        const __m512i s = _mm512_cvtepu8_epi32(_mm_maskz_loadu_epi8(k, src + 4 * i));
+        const __m512i d = _mm512_cvtepu8_epi32(_mm_maskz_loadu_epi8(k, dst + 4 * i));
+        const __m512i alpha = _mm512_shuffle_epi32(s, _MM_PERM_DDDD), dst_alpha = _mm512_shuffle_epi32(d, _MM_PERM_DDDD);
+        const __m512i inv_alpha = _mm512_sub_epi32(c255, alpha);
+        const __m512i denom = _mm512_add_epi32(_mm512_mullo_epi32(alpha, c255), _mm512_mullo_epi32(dst_alpha, inv_alpha));
+        const __m512i num =
+            _mm512_add_epi32(_mm512_mullo_epi32(_mm512_mullo_epi32(s, alpha), c255), _mm512_mullo_epi32(_mm512_mullo_epi32(d, dst_alpha), inv_alpha));
+        __m512i r = _mm512_cvtps_epi32(_mm512_div_ps(_mm512_cvtepi32_ps(num), _mm512_cvtepi32_ps(denom)));
+        const __m512i y = _mm512_add_epi32(denom, c128); // out alpha = round(denom / 255)
+        const __m512i out_alpha = _mm512_srli_epi32(_mm512_add_epi32(y, _mm512_srli_epi32(y, 8)), 8);
+        r = _mm512_mask_mov_epi32(r, 0x8888, out_alpha);
+        // both src and dst fully transparent => leave dst unchanged (this also discards the NaN from the 0/0 above)
+        r = _mm512_mask_mov_epi32(r, _mm512_cmpeq_epi32_mask(denom, zero), d);
+        _mm_mask_storeu_epi8(dst + 4 * i, k, _mm512_cvtepi32_epi8(r));
+    }
+    _mm256_zeroupper();
+}
+
+// }}}
+
 #define do_one_byte                                                                   \
     const uint8_t ch = src[pos++];                                                    \
     switch (decode_utf8(&d->state.cur, &d->state.codep, ch)) {                        \
@@ -376,4 +464,19 @@ printable_ascii_run_length_512(const uint32_t *chars UNUSED, const size_t sz UNU
     fatal("No AVX-512 implementation for this platform");
 }
 
+void
+blend_over_straight_512(uint8_t *dst UNUSED, const uint8_t *src UNUSED, size_t num_pixels UNUSED) {
+    fatal("No AVX-512 implementation for this platform");
+}
+
+void
+blend_over_opaque_512(uint8_t *dst UNUSED, unsigned dst_bpp UNUSED, const uint8_t *src UNUSED, size_t num_pixels UNUSED) {
+    fatal("No AVX-512 implementation for this platform");
+}
+
+void
+composite_alpha_mask_512(uint32_t *dst UNUSED, const uint8_t *mask UNUSED, size_t num_pixels UNUSED, uint32_t color_rgb UNUSED) {
+    fatal("No AVX-512 implementation for this platform");
+}
+
 #endif // x86-64

kitty/simd-string-impl.h

@@ -24,6 +24,9 @@ FUNC(utf8_decode_to_esc)(UTF8Decoder *d UNUSED, const uint8_t *src UNUSED, size_
     const uint8_t *haystack UNUSED, const size_t sz UNUSED, const uint8_t a UNUSED, const uint8_t b UNUSED) NOSIMD
     void FUNC(xor_data64)(const uint8_t key[64] UNUSED, uint8_t *data UNUSED, const size_t data_sz UNUSED) NOSIMD size_t
     FUNC(printable_ascii_run_length)(const uint32_t *chars UNUSED, const size_t sz UNUSED) NOSIMD
+    void FUNC(blend_over_straight)(uint8_t *dst UNUSED, const uint8_t *src UNUSED, size_t num_pixels UNUSED) NOSIMD
+    void FUNC(blend_over_opaque)(uint8_t *dst UNUSED, unsigned dst_bpp UNUSED, const uint8_t *src UNUSED, size_t num_pixels UNUSED) NOSIMD
+    void FUNC(composite_alpha_mask)(uint32_t *dst UNUSED, const uint8_t *mask UNUSED, size_t num_pixels UNUSED, uint32_t color_rgb UNUSED) NOSIMD
 #undef NOSIMD
 #else
 
@@ -96,6 +99,31 @@ _Pragma("clang diagnostic pop")
 #define create_zero_integer simde_mm_setzero_si128
 #define create_all_ones_integer() simde_mm_set1_epi64x(-1)
 #define zero_upper()
+#define set1_epi16(x) simde_mm_set1_epi16((short)(x))
+#define add_epi16 simde_mm_add_epi16
+#define mullo_epi16 simde_mm_mullo_epi16
+#define shift_right_by_bits16 simde_mm_srli_epi16
+#define unpacklo_epi8 simde_mm_unpacklo_epi8
+#define unpackhi_epi8 simde_mm_unpackhi_epi8
+#define packus_epi16 simde_mm_packus_epi16
+#define packus_epi32 simde_mm_packus_epi32
+#define add_epi32 simde_mm_add_epi32
+#define sub_epi32 simde_mm_sub_epi32
+#define mullo_epi32 simde_mm_mullo_epi32
+#define cmpeq_epi32 simde_mm_cmpeq_epi32
+#define max_epu32 simde_mm_max_epu32
+#define shuffle_epi32 simde_mm_shuffle_epi32
+#define cvtepi32_ps simde_mm_cvtepi32_ps
+#define cvtps_epi32 simde_mm_cvtps_epi32
+#define div_ps simde_mm_div_ps
+// set the same four 32-bit values in every 128-bit lane, a is the highest lane
+#define set_epi32_in_lanes(a, b, c, d) simde_mm_set_epi32(a, b, c, d)
+// replicate the alpha byte of each 4-byte pixel into all four of its bytes
+#define alpha_broadcast_pattern() set_epi8(15, 15, 15, 15, 11, 11, 11, 11, 7, 7, 7, 7, 3, 3, 3, 3)
+// widen the qth quarter (sizeof(integer_t)/4 bytes) of A into 32-bit lanes, q must be a literal
+#define widen_quarter(A, q) simde_mm_cvtepu8_epi32(simde_mm_srli_si128(A, 4 * (q)))
+// restore byte order after packus_epi32+packus_epi16 of four widened quarters
+#define fixup_packed_dword_order(v) (v)
 
     static inline int FUNC(is_zero)(const integer_t a) {
     return simde_mm_testz_si128(a, a);
@@ -170,6 +198,31 @@ w(right, one_byte, 1) w(right, two_bytes, 2) w(right, four_bytes, 4) w(right, ei
 #define shift_right_by_bits32 simde_mm256_srli_epi32
 #define create_zero_integer simde_mm256_setzero_si256
 #define create_all_ones_integer() simde_mm256_set1_epi64x(-1)
+#define set1_epi16(x) simde_mm256_set1_epi16((short)(x))
+#define add_epi16 simde_mm256_add_epi16
+#define mullo_epi16 simde_mm256_mullo_epi16
+#define shift_right_by_bits16 simde_mm256_srli_epi16
+#define unpacklo_epi8 simde_mm256_unpacklo_epi8
+#define unpackhi_epi8 simde_mm256_unpackhi_epi8
+#define packus_epi16 simde_mm256_packus_epi16
+#define packus_epi32 simde_mm256_packus_epi32
+#define add_epi32 simde_mm256_add_epi32
+#define sub_epi32 simde_mm256_sub_epi32
+#define mullo_epi32 simde_mm256_mullo_epi32
+#define cmpeq_epi32 simde_mm256_cmpeq_epi32
+#define max_epu32 simde_mm256_max_epu32
+#define shuffle_epi32 simde_mm256_shuffle_epi32
+#define cvtepi32_ps simde_mm256_cvtepi32_ps
+#define cvtps_epi32 simde_mm256_cvtps_epi32
+#define div_ps simde_mm256_div_ps
+// set the same four 32-bit values in every 128-bit lane, a is the highest lane
+#define set_epi32_in_lanes(a, b, c, d) simde_mm256_set_epi32(a, b, c, d, a, b, c, d)
+// replicate the alpha byte of each 4-byte pixel into all four of its bytes (in-lane indices)
+#define alpha_broadcast_pattern() set_epi8(15, 15, 15, 15, 11, 11, 11, 11, 7, 7, 7, 7, 3, 3, 3, 3, 15, 15, 15, 15, 11, 11, 11, 11, 7, 7, 7, 7, 3, 3, 3, 3)
+// widen the qth quarter (sizeof(integer_t)/4 bytes) of A into 32-bit lanes, q must be a literal
+#define widen_quarter(A, q) simde_mm256_cvtepu8_epi32(simde_mm_srli_si128(simde_mm256_extracti128_si256(A, (q) / 2), 8 * ((q) % 2)))
+// restore byte order after the per 128-bit lane packus_epi32+packus_epi16 of four widened quarters
+#define fixup_packed_dword_order(v) simde_mm256_permutevar8x32_epi32(v, simde_mm256_setr_epi32(0, 4, 1, 5, 2, 6, 3, 7))
 #define numbered_bytes() set_epi8(31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0)
 #define reverse_numbered_bytes() \
     simde_mm256_setr_epi8(31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0)
@@ -954,6 +1007,153 @@ FUNC(utf8_decode_to_esc)(UTF8Decoder *d, const uint8_t *src_data, size_t src_len
 #undef handle_trailing_bytes
 }
 
+// Pixel compositing {{{
+
+static inline integer_t
+FUNC(load_bytes_as_dwords)(const uint8_t *p) {
+#if KITTY_SIMD_LEVEL == 128
+    uint32_t v;
+    memcpy(&v, p, sizeof(v));
+    return simde_mm_cvtepu8_epi32(simde_mm_cvtsi32_si128((int32_t)v));
+#else
+    return simde_mm256_cvtepu8_epi32(simde_mm_loadl_epi64((const simde__m128i *)p));
+#endif
+}
+
+void
+FUNC(composite_alpha_mask)(uint32_t *dst, const uint8_t *mask, const size_t num_pixels, const uint32_t color_rgb) {
+    const uint32_t col = (color_rgb << 8) & 0xffffff00;
+    const integer_t col_vec = set1_epi32((int32_t)col), low_byte = set1_epi32(0xff);
+    const size_t px_per_iter = sizeof(integer_t) / 4;
+    size_t i = 0;
+    for (; i + px_per_iter <= num_pixels; i += px_per_iter) {
+        const integer_t m = FUNC(load_bytes_as_dwords)(mask + i);
+        const integer_t d = load_unaligned((const integer_t *)(dst + i));
+        store_unaligned((integer_t *)(dst + i), or_si(col_vec, max_epu32(m, and_si(d, low_byte))));
+    }
+    for (; i < num_pixels; i++) {
+        const uint32_t dst_alpha = dst[i] & 0xff, mask_alpha = mask[i];
+        dst[i] = col | MAX(mask_alpha, dst_alpha);
+    }
+    zero_upper();
+}
+
+static inline integer_t
+FUNC(div255_epu16)(const integer_t x) {
+    // rounding division of 16-bit lanes by 255, exact for values <= 65407, matches div255_round()
+    const integer_t y = add_epi16(x, set1_epi16(128));
+    return shift_right_by_bits16(add_epi16(y, shift_right_by_bits16(y, 8)), 8);
+}
+
+// Blend the 4-byte RGBA pixels in over onto the pixels in under, with under considered fully
+// opaque: out_c = round((over_c * alpha + under_c * (255 - alpha)) / 255) for each of the first
+// three channels, with the alpha bytes set to 255
+static inline integer_t
+FUNC(blend_opaque_pixels)(const integer_t under, const integer_t over) {
+    const integer_t zero = create_zero_integer();
+    const integer_t alpha = shuffle_epi8(over, alpha_broadcast_pattern());
+    const integer_t inv_alpha = xor_si(alpha, create_all_ones_integer()); // 255 - alpha in every byte
+    const integer_t lo = FUNC(div255_epu16)(
+        add_epi16(mullo_epi16(unpacklo_epi8(over, zero), unpacklo_epi8(alpha, zero)), mullo_epi16(unpacklo_epi8(under, zero), unpacklo_epi8(inv_alpha, zero))));
+    const integer_t hi = FUNC(div255_epu16)(
+        add_epi16(mullo_epi16(unpackhi_epi8(over, zero), unpackhi_epi8(alpha, zero)), mullo_epi16(unpackhi_epi8(under, zero), unpackhi_epi8(inv_alpha, zero))));
+    // the per 128-bit lane unpacks and pack are symmetric so byte order is preserved
+    return or_si(packus_epi16(lo, hi), set1_epi32((int32_t)0xff000000));
+}
+
+#if KITTY_SIMD_LEVEL == 128
+#define blend_opaque_pixels_128bit FUNC(blend_opaque_pixels)
+#else
+// 128-bit version of blend_opaque_pixels for the 3 bytes per pixel destination case, whose
+// expand/compact shuffles cannot cross 128-bit lane boundaries
+static inline simde__m128i
+blend_opaque_pixels_128bit(const simde__m128i under, const simde__m128i over) {
+    const simde__m128i zero = simde_mm_setzero_si128(), c128 = simde_mm_set1_epi16(128);
+    const simde__m128i alpha = simde_mm_shuffle_epi8(over, simde_mm_set_epi8(15, 15, 15, 15, 11, 11, 11, 11, 7, 7, 7, 7, 3, 3, 3, 3));
+    const simde__m128i inv_alpha = simde_mm_xor_si128(alpha, simde_mm_set1_epi64x(-1));
+    simde__m128i t, y;
+    t = simde_mm_add_epi16(
+        simde_mm_mullo_epi16(simde_mm_unpacklo_epi8(over, zero), simde_mm_unpacklo_epi8(alpha, zero)),
+        simde_mm_mullo_epi16(simde_mm_unpacklo_epi8(under, zero), simde_mm_unpacklo_epi8(inv_alpha, zero)));
+    y = simde_mm_add_epi16(t, c128);
+    const simde__m128i lo = simde_mm_srli_epi16(simde_mm_add_epi16(y, simde_mm_srli_epi16(y, 8)), 8);
+    t = simde_mm_add_epi16(
+        simde_mm_mullo_epi16(simde_mm_unpackhi_epi8(over, zero), simde_mm_unpackhi_epi8(alpha, zero)),
+        simde_mm_mullo_epi16(simde_mm_unpackhi_epi8(under, zero), simde_mm_unpackhi_epi8(inv_alpha, zero)));
+    y = simde_mm_add_epi16(t, c128);
+    const simde__m128i hi = simde_mm_srli_epi16(simde_mm_add_epi16(y, simde_mm_srli_epi16(y, 8)), 8);
+    return simde_mm_or_si128(simde_mm_packus_epi16(lo, hi), simde_mm_set1_epi32((int32_t)0xff000000));
+}
+#endif
+
+void
+FUNC(blend_over_opaque)(uint8_t *dst, const unsigned dst_bpp, const uint8_t *src, const size_t num_pixels) {
+    size_t i = 0;
+    if (dst_bpp == 4) {
+        const size_t px_per_iter = sizeof(integer_t) / 4;
+        for (; i + px_per_iter <= num_pixels; i += px_per_iter) {
+            const integer_t s = load_unaligned((const integer_t *)(src + 4 * i)), d = load_unaligned((const integer_t *)(dst + 4 * i));
+            store_unaligned((integer_t *)(dst + 4 * i), FUNC(blend_opaque_pixels)(d, s));
+        }
+    } else {
+        // Expand each group of four 3-byte pixels to 4-byte pixels with a byte shuffle, blend, and
+        // compact back with another shuffle, using 128-bit registers at every level since the
+        // shuffles cannot cross 128-bit lane boundaries. The 16-byte destination loads cover 4
+        // bytes beyond the four pixels being processed, which stay inside dst as long as at least
+        // six pixels remain. Only the 12 blended bytes are stored, both to leave the extra bytes
+        // untouched and because a 16-byte store would overlap the next iteration's load, stalling
+        // on failed store to load forwarding.
+        const simde__m128i expand = simde_mm_set_epi8(-1, 11, 10, 9, -1, 8, 7, 6, -1, 5, 4, 3, -1, 2, 1, 0);
+        const simde__m128i compact = simde_mm_set_epi8(-1, -1, -1, -1, 14, 13, 12, 10, 9, 8, 6, 5, 4, 2, 1, 0);
+        for (; i + 6 <= num_pixels; i += 4) {
+            const simde__m128i s = simde_mm_loadu_si128((const simde__m128i *)(src + 4 * i));
+            const simde__m128i d = simde_mm_loadu_si128((const simde__m128i *)(dst + 3 * i));
+            const simde__m128i blended = simde_mm_shuffle_epi8(blend_opaque_pixels_128bit(simde_mm_shuffle_epi8(d, expand), s), compact);
+            simde_mm_storel_epi64((simde__m128i *)(dst + 3 * i), blended);
+            const uint32_t last = (uint32_t)simde_mm_extract_epi32(blended, 2);
+            memcpy(dst + 3 * i + 8, &last, sizeof(last));
+        }
+    }
+    for (; i < num_pixels; i++) blend_pixel_over_opaque(dst + dst_bpp * i, src + 4 * i, dst_bpp);
+    zero_upper();
+}
+
+// Straight alpha over blend of pixels unpacked into 32-bit lanes, four consecutive lanes per
+// pixel. Matches the arithmetic of blend_pixel_over_straight() bit for bit: exact integer
+// numerator and denominator, IEEE single precision division, round to nearest.
+static inline integer_t
+FUNC(blend_straight_dwords)(const integer_t d, const integer_t s) {
+    const integer_t c255 = set1_epi32(255);
+    const integer_t alpha = shuffle_epi32(s, _MM_SHUFFLE(3, 3, 3, 3)), dst_alpha = shuffle_epi32(d, _MM_SHUFFLE(3, 3, 3, 3));
+    const integer_t inv_alpha = sub_epi32(c255, alpha);
+    const integer_t denom = add_epi32(mullo_epi32(alpha, c255), mullo_epi32(dst_alpha, inv_alpha));
+    const integer_t num = add_epi32(mullo_epi32(mullo_epi32(s, alpha), c255), mullo_epi32(mullo_epi32(d, dst_alpha), inv_alpha));
+    integer_t r = cvtps_epi32(div_ps(cvtepi32_ps(num), cvtepi32_ps(denom)));
+    const integer_t y = add_epi32(denom, set1_epi32(128)); // out alpha = round(denom / 255)
+    const integer_t out_alpha = shift_right_by_bits32(add_epi32(y, shift_right_by_bits32(y, 8)), 8);
+    r = blendv_epi8(r, out_alpha, set_epi32_in_lanes(-1, 0, 0, 0));
+    // both src and dst fully transparent => leave dst unchanged (this also discards the NaN from the 0/0 above)
+    return blendv_epi8(r, d, cmpeq_epi32(denom, create_zero_integer()));
+}
+
+void
+FUNC(blend_over_straight)(uint8_t *dst, const uint8_t *src, const size_t num_pixels) {
+    const size_t px_per_iter = sizeof(integer_t) / 4;
+    size_t i = 0;
+    for (; i + px_per_iter <= num_pixels; i += px_per_iter) {
+        const integer_t s = load_unaligned((const integer_t *)(src + 4 * i)), d = load_unaligned((const integer_t *)(dst + 4 * i));
+        const integer_t r0 = FUNC(blend_straight_dwords)(widen_quarter(d, 0), widen_quarter(s, 0));
+        const integer_t r1 = FUNC(blend_straight_dwords)(widen_quarter(d, 1), widen_quarter(s, 1));
+        const integer_t r2 = FUNC(blend_straight_dwords)(widen_quarter(d, 2), widen_quarter(s, 2));
+        const integer_t r3 = FUNC(blend_straight_dwords)(widen_quarter(d, 3), widen_quarter(s, 3));
+        store_unaligned((integer_t *)(dst + 4 * i), fixup_packed_dword_order(packus_epi16(packus_epi32(r0, r1), packus_epi32(r2, r3))));
+    }
+    for (; i < num_pixels; i++) blend_pixel_over_straight(dst + 4 * i, src + 4 * i);
+    zero_upper();
+}
+
+#undef blend_opaque_pixels_128bit
+// }}}
 
 #undef FUNC
 #undef integer_t
@@ -1012,4 +1212,25 @@ FUNC(utf8_decode_to_esc)(UTF8Decoder *d, const uint8_t *src_data, size_t src_len
 #undef is_zero
 #undef zero_upper
 #undef print_register_as_bytes
+#undef set1_epi16
+#undef add_epi16
+#undef mullo_epi16
+#undef shift_right_by_bits16
+#undef unpacklo_epi8
+#undef unpackhi_epi8
+#undef packus_epi16
+#undef packus_epi32
+#undef add_epi32
+#undef sub_epi32
+#undef mullo_epi32
+#undef cmpeq_epi32
+#undef max_epu32
+#undef shuffle_epi32
+#undef cvtepi32_ps
+#undef cvtps_epi32
+#undef div_ps
+#undef set_epi32_in_lanes
+#undef alpha_broadcast_pattern
+#undef widen_quarter
+#undef fixup_packed_dword_order
 #endif // KITTY_NO_SIMD

kitty/simd-string.c

@@ -39,6 +39,42 @@ find_either_of_two_bytes(const uint8_t *haystack, const size_t sz, const uint8_t
 }
 // }}}
 
+// pixel compositing {{{
+static void
+blend_over_straight_scalar(uint8_t *dst, const uint8_t *src, const size_t num_pixels) {
+    for (size_t i = 0; i < num_pixels; i++) blend_pixel_over_straight(dst + 4 * i, src + 4 * i);
+}
+static void (*blend_over_straight_impl)(uint8_t *, const uint8_t *, size_t) = blend_over_straight_scalar;
+void
+blend_over_straight(uint8_t *dst, const uint8_t *src, size_t num_pixels) {
+    blend_over_straight_impl(dst, src, num_pixels);
+}
+
+static void
+blend_over_opaque_scalar(uint8_t *dst, const unsigned dst_bpp, const uint8_t *src, const size_t num_pixels) {
+    for (size_t i = 0; i < num_pixels; i++) blend_pixel_over_opaque(dst + dst_bpp * i, src + 4 * i, dst_bpp);
+}
+static void (*blend_over_opaque_impl)(uint8_t *, unsigned, const uint8_t *, size_t) = blend_over_opaque_scalar;
+void
+blend_over_opaque(uint8_t *dst, unsigned dst_bpp, const uint8_t *src, size_t num_pixels) {
+    blend_over_opaque_impl(dst, dst_bpp, src, num_pixels);
+}
+
+static void
+composite_alpha_mask_scalar(uint32_t *dst, const uint8_t *mask, const size_t num_pixels, const uint32_t color_rgb) {
+    const uint32_t col = (color_rgb << 8) & 0xffffff00;
+    for (size_t i = 0; i < num_pixels; i++) {
+        const uint32_t dst_alpha = dst[i] & 0xff, mask_alpha = mask[i];
+        dst[i] = col | MAX(mask_alpha, dst_alpha);
+    }
+}
+static void (*composite_alpha_mask_impl)(uint32_t *, const uint8_t *, size_t, uint32_t) = composite_alpha_mask_scalar;
+void
+composite_alpha_mask(uint32_t *dst, const uint8_t *mask, size_t num_pixels, uint32_t color_rgb) {
+    composite_alpha_mask_impl(dst, mask, num_pixels, color_rgb);
+}
+// }}}
+
 // printable_ascii_run_length {{{
 static size_t
 printable_ascii_run_length_scalar(const uint32_t *chars, const size_t sz) {
@@ -212,6 +248,85 @@ test_xor64(PyObject *self UNUSED, PyObject *args) {
     return ans;
 }
 
+static PyObject *
+test_blend_pixels(PyObject *self UNUSED, PyObject *args) {
+    const char *kind;
+    RAII_PY_BUFFER(dst);
+    RAII_PY_BUFFER(src);
+    int which_function = 0, align_offset = 0;
+    unsigned long color = 0xffffff;
+    if (!PyArg_ParseTuple(args, "ss*s*|iki", &kind, &dst, &src, &which_function, &color, &align_offset)) return NULL;
+    align_offset &= 63;
+    void (*straight)(uint8_t *, const uint8_t *, size_t) = NULL;
+    void (*opaque)(uint8_t *, unsigned, const uint8_t *, size_t) = NULL;
+    void (*mask)(uint32_t *, const uint8_t *, size_t, uint32_t) = NULL;
+    unsigned dst_bpp = 4;
+    size_t num_pixels;
+    if (strcmp(kind, "straight") == 0) {
+        num_pixels = (size_t)src.len / 4;
+        switch (which_function) {
+            case 0: straight = blend_over_straight; break;
+            case 1: straight = blend_over_straight_scalar; break;
+            case 2: straight = blend_over_straight_128; break;
+            case 3: straight = blend_over_straight_256; break;
+            case 4: straight = blend_over_straight_512; break;
+        }
+    } else if (strcmp(kind, "opaque") == 0 || strcmp(kind, "opaque3") == 0) {
+        if (strcmp(kind, "opaque3") == 0) dst_bpp = 3;
+        num_pixels = (size_t)src.len / 4;
+        switch (which_function) {
+            case 0: opaque = blend_over_opaque; break;
+            case 1: opaque = blend_over_opaque_scalar; break;
+            case 2: opaque = blend_over_opaque_128; break;
+            case 3: opaque = blend_over_opaque_256; break;
+            case 4: opaque = blend_over_opaque_512; break;
+        }
+    } else if (strcmp(kind, "mask") == 0) {
+        num_pixels = (size_t)src.len;
+        switch (which_function) {
+            case 0: mask = composite_alpha_mask; break;
+            case 1: mask = composite_alpha_mask_scalar; break;
+            case 2: mask = composite_alpha_mask_128; break;
+            case 3: mask = composite_alpha_mask_256; break;
+            case 4: mask = composite_alpha_mask_512; break;
+        }
+    } else {
+        PyErr_Format(PyExc_KeyError, "Unknown kind: %s", kind);
+        return NULL;
+    }
+    if (!straight && !opaque && !mask) {
+        PyErr_SetString(PyExc_ValueError, "Unknown which_function");
+        return NULL;
+    }
+    const size_t expected_dst_sz = num_pixels * (mask ? 4 : dst_bpp);
+    if ((size_t)dst.len != expected_dst_sz) {
+        PyErr_Format(PyExc_ValueError, "dst must be %zu bytes not %zd", expected_dst_sz, dst.len);
+        return NULL;
+    }
+    uint8_t *abuf;
+    if (posix_memalign((void **)&abuf, 64, 192 + dst.len) != 0) return PyErr_NoMemory();
+    uint8_t *p = abuf + 64 + align_offset;
+    memset(abuf, '<', 64 + align_offset);
+    memcpy(p, dst.buf, dst.len);
+    memset(p + dst.len, '>', 64);
+    if (straight) straight(p, src.buf, num_pixels);
+    else if (opaque) opaque(p, dst_bpp, src.buf, num_pixels);
+    else mask((uint32_t *)p, src.buf, num_pixels, (uint32_t)color);
+    PyObject *ans = NULL;
+    for (int i = 0; i < 64 + align_offset; i++)
+        if (abuf[i] != '<') {
+            PyErr_SetString(PyExc_SystemError, "blend wrote before start of data region");
+            break;
+        }
+    for (int i = 0; i < 64; i++)
+        if (p[i + dst.len] != '>') {
+            PyErr_SetString(PyExc_SystemError, "blend wrote after end of data region");
+            break;
+        }
+    if (!PyErr_Occurred()) ans = PyBytes_FromStringAndSize((const char *)p, dst.len);
+    free(abuf);
+    return ans;
+}
 
 // }}}
 
@@ -220,6 +335,7 @@ static PyMethodDef module_methods[] = {
     METHODB(test_find_either_of_two_bytes, METH_VARARGS),
     METHODB(test_printable_ascii_run_length, METH_VARARGS),
     METHODB(test_xor64, METH_VARARGS),
+    METHODB(test_blend_pixels, METH_VARARGS),
     {NULL, NULL, 0, NULL} /* Sentinel */
 };
 
@@ -280,6 +396,9 @@ init_simd(void *x) {
         find_either_of_two_bytes_impl = find_either_of_two_bytes_512;
         xor_data64_impl = xor_data64_512;
         printable_ascii_run_length_impl = printable_ascii_run_length_512;
+        blend_over_straight_impl = blend_over_straight_512;
+        blend_over_opaque_impl = blend_over_opaque_512;
+        composite_alpha_mask_impl = composite_alpha_mask_512;
     } else {
         A(has_avx512, False);
     }
@@ -289,6 +408,9 @@ init_simd(void *x) {
         if (utf8_decode_to_esc_impl == utf8_decode_to_esc_scalar) utf8_decode_to_esc_impl = utf8_decode_to_esc_256;
         if (xor_data64_impl == xor_data64_scalar) xor_data64_impl = xor_data64_256;
         if (printable_ascii_run_length_impl == printable_ascii_run_length_scalar) printable_ascii_run_length_impl = printable_ascii_run_length_256;
+        if (blend_over_straight_impl == blend_over_straight_scalar) blend_over_straight_impl = blend_over_straight_256;
+        if (blend_over_opaque_impl == blend_over_opaque_scalar) blend_over_opaque_impl = blend_over_opaque_256;
+        if (composite_alpha_mask_impl == composite_alpha_mask_scalar) composite_alpha_mask_impl = composite_alpha_mask_256;
     } else {
         A(has_avx2, False);
     }
@@ -298,6 +420,9 @@ init_simd(void *x) {
         if (utf8_decode_to_esc_impl == utf8_decode_to_esc_scalar) utf8_decode_to_esc_impl = utf8_decode_to_esc_128;
         if (xor_data64_impl == xor_data64_scalar) xor_data64_impl = xor_data64_128;
         if (printable_ascii_run_length_impl == printable_ascii_run_length_scalar) printable_ascii_run_length_impl = printable_ascii_run_length_128;
+        if (blend_over_straight_impl == blend_over_straight_scalar) blend_over_straight_impl = blend_over_straight_128;
+        if (blend_over_opaque_impl == blend_over_opaque_scalar) blend_over_opaque_impl = blend_over_opaque_128;
+        if (composite_alpha_mask_impl == composite_alpha_mask_scalar) composite_alpha_mask_impl = composite_alpha_mask_128;
     } else {
         A(has_sse4_2, False);
     }

kitty/simd-string.h

@@ -9,6 +9,7 @@
 #include "data-types.h"
 #include <stddef.h>
 #include <stdalign.h>
+#include <math.h>
 
 typedef void (*control_byte_callback)(void *data, uint8_t ch);
 typedef void (*output_chars_callback)(void *data, const uint32_t *chars, unsigned count);
@@ -62,6 +63,52 @@ void xor_data64(const uint8_t key[64], uint8_t *data, const size_t data_sz);
 // Returns the length of the prefix of chars that contains only printable ASCII codepoints (32 <= ch <= 126)
 size_t printable_ascii_run_length(const uint32_t *chars, const size_t sz);
 
+// Composite src over dst in place (Porter-Duff over with straight, aka non-premultiplied, alpha).
+// Both are arrays of num_pixels 4-byte RGBA pixels, alpha in the fourth byte of each pixel.
+void blend_over_straight(uint8_t *dst, const uint8_t *src, size_t num_pixels);
+
+// As above, but dst is treated as fully opaque so no tracking of dst alpha is needed. dst_bpp must
+// be 4 (RGBA pixels whose alpha bytes are set to 255) or 3 (RGB pixels).
+void blend_over_opaque(uint8_t *dst, unsigned dst_bpp, const uint8_t *src, size_t num_pixels);
+
+// Composite an 8-bit alpha mask onto 32-bit pixels of the form (r << 24) | (g << 16) | (b << 8) | alpha:
+// dst[i] = ((color_rgb << 8) & 0xffffff00) | MAX(mask[i], dst[i] & 0xff)
+void composite_alpha_mask(uint32_t *dst, const uint8_t *mask, size_t num_pixels, uint32_t color_rgb);
+
+// Scalar per-pixel helpers shared by the scalar implementations and the tail handling of the SIMD
+// implementations, which must produce bit identical results to them.
+
+static inline uint32_t
+div255_round(const uint32_t x) {
+    // rounding division by 255, exact for x <= 65407 which covers x <= 255 * 255
+    const uint32_t y = x + 128;
+    return (y + (y >> 8)) >> 8;
+}
+
+static inline void
+blend_pixel_over_opaque(uint8_t *dst, const uint8_t *src, const unsigned dst_bpp) {
+    const uint32_t alpha = src[3], inv_alpha = 255 - alpha;
+    for (unsigned c = 0; c < 3; c++) dst[c] = (uint8_t)div255_round(src[c] * alpha + dst[c] * inv_alpha);
+    if (dst_bpp == 4) dst[3] = 255;
+}
+
+static inline void
+blend_pixel_over_straight(uint8_t *dst, const uint8_t *src) {
+    const uint32_t alpha = src[3], dst_alpha = dst[3], inv_alpha = 255 - alpha;
+    // The over operator on straight alpha in 8-bit integer arithmetic. Multiplying out the
+    // normalizations by 255, the output alpha and color channels are:
+    // out_alpha = (alpha * 255 + dst_alpha * (255 - alpha)) / 255
+    // out_c = (src_c * alpha * 255 + dst_c * dst_alpha * (255 - alpha)) / (alpha * 255 + dst_alpha * (255 - alpha))
+    const uint32_t denom = alpha * 255 + dst_alpha * inv_alpha;
+    if (!denom) return; // both src and dst fully transparent, leave dst unchanged
+    for (unsigned c = 0; c < 3; c++) {
+        const uint32_t num = src[c] * alpha * 255 + dst[c] * dst_alpha * inv_alpha;
+        // float division and round to nearest, matching the SIMD implementations bit for bit
+        dst[c] = (uint8_t)lrintf((float)num / (float)denom);
+    }
+    dst[3] = (uint8_t)div255_round(denom);
+}
+
 // SIMD implementations, internal use
 bool utf8_decode_to_esc_128(UTF8Decoder *d, const uint8_t *src, size_t src_sz);
 bool utf8_decode_to_esc_256(UTF8Decoder *d, const uint8_t *src, size_t src_sz);
@@ -75,3 +122,12 @@ void xor_data64_512(const uint8_t key[64], uint8_t *data, const size_t data_sz);
 size_t printable_ascii_run_length_128(const uint32_t *chars, const size_t sz);
 size_t printable_ascii_run_length_256(const uint32_t *chars, const size_t sz);
 size_t printable_ascii_run_length_512(const uint32_t *chars, const size_t sz);
+void blend_over_straight_128(uint8_t *dst, const uint8_t *src, size_t num_pixels);
+void blend_over_straight_256(uint8_t *dst, const uint8_t *src, size_t num_pixels);
+void blend_over_straight_512(uint8_t *dst, const uint8_t *src, size_t num_pixels);
+void blend_over_opaque_128(uint8_t *dst, unsigned dst_bpp, const uint8_t *src, size_t num_pixels);
+void blend_over_opaque_256(uint8_t *dst, unsigned dst_bpp, const uint8_t *src, size_t num_pixels);
+void blend_over_opaque_512(uint8_t *dst, unsigned dst_bpp, const uint8_t *src, size_t num_pixels);
+void composite_alpha_mask_128(uint32_t *dst, const uint8_t *mask, size_t num_pixels, uint32_t color_rgb);
+void composite_alpha_mask_256(uint32_t *dst, const uint8_t *mask, size_t num_pixels, uint32_t color_rgb);
+void composite_alpha_mask_512(uint32_t *dst, const uint8_t *mask, size_t num_pixels, uint32_t color_rgb);

kitty_tests/parser.py

@@ -13,6 +13,7 @@ from kitty.fast_data_types import (
     has_avx2,
     has_avx512,
     has_sse4_2,
+    test_blend_pixels,
     test_find_either_of_two_bytes,
     test_printable_ascii_run_length,
     test_utf8_decode_to_sentinel,
@@ -836,6 +837,42 @@ class TestParser(BaseTest):
                     data = base + base_data[:extra]
                     t(key, data, align_offset)
 
+    def test_blend_pixels(self):
+        rng = random.Random(0x1B7F5)
+        impls = []
+        if has_sse4_2:
+            impls.append(2)
+        if has_avx2:
+            impls.append(3)
+        if has_avx512:
+            impls.append(4)
+        impls.append(0)
+
+        def t(kind, dst, src, color=0xFFFFFF, align_offsets=(0, 1, 2, 3, 7, 13)):
+            expected = test_blend_pixels(kind, dst, src, 1, color, 0)
+            for which_function in impls:
+                for align_offset in align_offsets:
+                    actual = test_blend_pixels(kind, dst, src, which_function, color, align_offset)
+                    self.ae(expected, actual, f'{kind=} {which_function=} {align_offset=} num_pixels={len(src) // 4}')
+
+        for num_pixels in (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 15, 16, 17, 31, 32, 33, 63, 64, 65, 255, 256, 257):
+            src = bytearray(rng.getrandbits(8) for _ in range(4 * num_pixels))
+            # ensure fully transparent and fully opaque source pixels are well represented
+            for i in range(num_pixels):
+                if i % 5 == 0:
+                    src[4 * i + 3] = 0
+                elif i % 5 == 1:
+                    src[4 * i + 3] = 255
+            src = bytes(src)
+            dst4 = bytes(rng.getrandbits(8) for _ in range(4 * num_pixels))
+            dst3 = bytes(rng.getrandbits(8) for _ in range(3 * num_pixels))
+            mask = bytes(rng.getrandbits(8) for _ in range(num_pixels))
+            t('straight', dst4, src)
+            t('opaque', dst4, src)
+            t('opaque3', dst3, src)
+            # keep the destination pixels 4-byte aligned for the mask variant
+            t('mask', dst4, mask, color=0x123456, align_offsets=(0, 4, 12, 36))
+
     def test_esc_codes(self):
         s = self.create_screen()
         pb = partial(self.parse_bytes_dump, s)