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) composite_normed = (1 - alpha) * bg_normed + alpha * red_normed_rgb composite = (composite_normed * 255).round().astype(np.uint8)