Speed up RGBA-stage image resampling - #32275

PR summary

Closes #32240.

Redrawing an imshow image got about 2x slower between 3.10.9 and 3.11. @ayshih traced this to #29776, which moved resampling into premultiplied alpha space. The premultiply and demultiply are both full-array NumPy passes, but in their current form, they can't use NumPy's fast paths.

Demultiply is the more expensive one:

np.divide(res[..., :3], res[..., 3:], out=res[..., :3], where=res[..., 3:] != 0)

This takes ~16 ms per draw in the code I used to reproduce the issue in #32240. Two reasons it's slow: out=res[..., :3] is a stride-4 view (so contiguous inner loops can't be used), and passing an array to where= is itself expensive. Running it on a 462x1560x4 float64 array takes 16.3 ms as written, 8.3 ms without where=, and only 2.8 ms operating on the whole contiguous array. An all-True mask array still costs 2.3-2.8x depending on whether the destination is strided, so the cost just comes from taking the masked path itself.

For both the demultiply and premultiply steps, this PR implements contiguous whole-array operations that save the alpha channel and then restore it.

• Demultiply: save the alpha values, divide the whole array by a version of alpha where zeros have been replaced by ones, then put the original alpha values back. This is equivalent to the where= form, and while dividing alpha by itself involves more operations, the fact that it runs contiguously means it's much faster.
• Premultiply: same idea. Making the multiply contiguous takes it from 13.4 ms to 4.6 ms per draw in my reproduction code. On top of that, when the alpha channel is all 1 (as with colormapped or opaque RGBA inputs), it can be skipped entirely for the cost of a 0.6 ms check, taking the line from 13.4 ms to 0.6 ms.

One positive side effect is that the line A = A.copy() is no longer needed. It existed to stop A[..., :3] *= A[..., 3:] from writing into the caller's array, but now A is not modified in place. When alpha is all 1 it's left untouched, and otherwise the multiplication just builds a new array.

Timings

Median of 15 repeats, same code from #32240 (macOS arm64, NumPy 2.5.2)

Verification

A new test, test_rgba_stage_demultiplies_resampled_alpha, checks that RGB values survive the premultiply/demultiply round trip for several alpha values and that a fully transparent image never divides by zero.

To double-check that rendered output is unchanged, I rendered a large set of configurations and compared the RGBA canvas buffers against main and against the 3.11.1 wheel. No differences were observed.

import sys, numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

rng = np.random.default_rng(1234)
alpha_arr = rng.random((7, 11))
arrays = {
    "scalar_f":    rng.random((7, 11)),
    "scalar_u8":   (rng.random((7, 11)) * 255).astype(np.uint8),
    "scalar_mask": np.ma.masked_array(rng.random((7, 11)), rng.random((7, 11)) < .3),
    "rgb_f":       rng.random((7, 11, 3)),
    "rgba_f":      rng.random((7, 11, 4)),
    "rgba_opaque": np.dstack([rng.random((7, 11, 3)), np.ones((7, 11))]),
    "rgb_u8":      (rng.random((7, 11, 3)) * 255).astype(np.uint8),
    "rgba_u8":     (rng.random((7, 11, 4)) * 255).astype(np.uint8),
    "wide":        rng.random((60, 900)),
}

out = {}
for name, arr in arrays.items():
    for interp in ["none", "nearest", "bilinear", "bicubic", "lanczos", "antialiased"]:
        for alpha in [None, 0.4, "array"]:
            for stage in ["auto", "data", "rgba"]:
                for origin in ["upper", "lower"]:
                    if alpha == "array" and arr.shape[:2] != alpha_arr.shape:
                        continue
                    fig = plt.figure(figsize=(2.37, 1.63), dpi=100)
                    ax = fig.add_axes((0.13, 0.11, 0.7, 0.72))
                    ax.imshow(arr, interpolation=interp, origin=origin, aspect="auto",
                              alpha=alpha_arr if alpha == "array" else alpha,
                              interpolation_stage=stage)
                    fig.canvas.draw()
                    out[f"{name}|{interp}|{alpha}|{stage}|{origin}"] = \
                        np.asarray(fig.canvas.buffer_rgba()).copy()
                    plt.close(fig)

np.savez_compressed(sys.argv[1], **out)
print(f"{len(out)} configurations -> {sys.argv[1]}")

a, b = np.load("before.npz"), np.load("after.npz")
print(sum(np.array_equal(a[k], b[k]) for k in a.files), "of", len(a.files), "identical")

AI Disclosure

I used Claude Code in several places and thoroughly reviewed its output. It profiled 3.10.9 against 3.11.1 to locate the two affected lines, benchmarked the approach implemented here, wrote the patch and the test, and ran the verification. I wrote or revised all prose, reviewed all code, and re-ran the verification script.

PR quality check

• Use an expressive title
• New and changed code is tested
• [N/A] Plotting related features are demonstrated in an example
• [N/A] New features and API changes have release notes
• [N/A] Documentation complies with general and docstring guidelines