import numpy as np
# dummy data
nx, ny = 512, 512
i, j = np.indices((nx, ny))
red = np.zeros((nx, ny, 4), dtype=np.uint8)
red[..., 0] = 255
red[..., -1] = 255 * np.sin(2*i/nx*np.pi) * np.cos(3*j/ny*np.pi)
red_normed = red / 255 # shape (nx, ny, 4), dtype float
alpha = red_normed[..., -1:] # shape (nx, ny, 1) for broadcasting
red_normed_rgb = red_normed[..., :-1] # shape (nx, ny, 3)
#bg_normed = np.zeros_like(red_normed_rgb) # shape (nx, ny, 3) <-- black background
bg_normed = np.ones_like(red_normed_rgb) # shape (nx, ny, 3)
(see full text)