Overlay two numpy arrays treating fourth plane as alpha level [duplicate]

Question

I have two numpy arrays of shape (256, 256, 4). I would like to treat the fourth 256 x 256 plane as an alpha level, and export an image where these arrays have been overlayed.

Code example:

import numpy as np
from skimage import io

fg = np.ndarray((256, 256, 4), dtype=np.uint8)
one_plane = np.random.standard_normal((256, 256)) * 100 + 128
fg[:,:,0:3] = np.tile(one_plane, 3).reshape((256, 256, 3), order='F')
fg[:, :, 3] = np.zeros((256, 256), dtype=np.uint8)
fg[0:128, 0:128, 3] = np.ones((128, 128), dtype=np.uint8) * 255
fg[128:256, 128:256, 3] = np.ones((128, 128), dtype=np.uint8) * 128

bg = np.ndarray((256, 256, 4), dtype=np.uint8)
bg[:,:,0:3] = np.random.standard_normal((256, 256, 3)) * 100 + 128
bg[:, :, 3] = np.ones((256, 256), dtype=np.uint8) * 255

io.imsave('test_fg.png', fg)
io.imsave('test_bg.png', bg)

This creates two images, fg:

and bg:

:

I would like to be able to overlay the fg onto the bg. That is, the final image should have grey in the top left (because there the alpha of fg is 1), a blend of grey and colour noise in the bottom right, and pure colour noise in the other quandrants. I am looking for something like an add function that gives me a new np array.

Note that I don't think this is the same as this answer, which uses matplotlib.pyplot.plt to overlay the images and fiddles around with colour maps. I don't think I should need to fiddle with colour maps here, but maybe the answer is I do.

The reason I would like a new np.array returned by the operation is because I want to do this iteratively with many images, overlayed in order.

Answer 1

Alpha blending is usually done using the Porter & Duff equations:

where src and dst would correspond to your foreground and background images, and the A and RGB pixel values are assumed to be floating point, in the range [0, 1].

For your specific example:

src_rgb = fg[..., :3].astype(np.float32) / 255.0
src_a = fg[..., 3].astype(np.float32) / 255.0
dst_rgb = bg[..., :3].astype(np.float32) / 255.0
dst_a = bg[..., 3].astype(np.float32) / 255.0

out_a = src_a + dst_a*(1.0-src_a)
out_rgb = (src_rgb*src_a[..., None]
           + dst_rgb*dst_a[..., None]*(1.0-src_a[..., None])) / out_a[..., None]

out = np.zeros_like(bg)
out[..., :3] = out_rgb * 255
out[..., 3] = out_a * 255

Output: