如何有效地将随机垂直段添加到 numpy 数组中？

Question

我正在尝试使用 NumPy 模拟下雨，他们说图像超过一千个字，所以这里是一个超过两千个字的描述：

我已经写好了代码，但是我觉得我的实现效率很低，所以我想知道NumPy有没有内置函数可以加速这个过程：

import numpy as np
from PIL import Image
from random import random, randbytes

def rain(width, strikes=360, color=True, lw=3):
    assert not width % 16
    height = int(width * 9 / 16)
    img = np.zeros((height, width, 3), dtype=np.uint8)
    half = height / 2
    for i in range(strikes):
        x = round(random() * width)
        y = round(height - random() * half)
        x1 = min(x + lw, width - 1)
        if color:
            rgb = list(randbytes(3))
        else:
            rgb = (178, 255, 255)
        img[0:y, x:x1] = rgb
    
    return img

img1 = Image.fromarray(rain(1920))
img1.show()
img1.save('D:/rain.jpg', format='jpeg', quality=80, optimize=True)

img2 = Image.fromarray(rain(1920, color=False))
img2.show()
img2.save('D:/rain_1.jpg', format='jpeg', quality=80, optimize=True)

Answer 1

因此，使用 NumPy 加速代码的最简单方法是利用广播和 element-by-element 操作，这样可以避免效率较低的 for-loops。下面是我的算法 (rain2) 和 OP 的 (rain1) 之间的性能比较：

import numpy.random as npr
from random import random, randbytes
from PIL import Image
import profile

def rain1(width, strikes=360, color=True, lw=3):
    assert not width % 16
    height = int(width * 9 / 16)
    img = np.zeros((height, width, 3), dtype=np.uint8)
    half = height / 2
    for i in range(strikes):
        x = round(random() * width)
        y = round(height - random() * half)
        x1 = min(x + lw, width - 1)
        if color:
            rgb = list(randbytes(3))
        else:
            rgb = (178, 255, 255)
    img[0:y, x:x1] = rgb
    
    return img


def rain2(width, strikes=360, color=True, lw=3):
    assert not width % 16
    height = width*9//16
    [inds,] = np.indices((width,))
    img = np.zeros((height, width, 4), dtype=np.uint8)
    img[:,:,3] = 255
    half = height/2
    # randint from numpy.random lets you 
    # define a lower and upper bound, 
    # and number of points.
    x = list(set(npr.randint(0, width-lw-1, (strikes,))))
    x = np.sort(x)
    y = npr.randint(half, height, (len(x),))
    if color:
        rgb = npr.randint(0, 255, (len(x), 3), dtype=np.uint8)
    else:
        rgb = np.array([178, 255, 255], dtype=np.uint8)
    for offset in range(lw):
        img[:,x+offset,3] = 0
        img[:,x+offset,:3] = rgb
    for xi, yi in zip(x, y):
        img[0:yi,xi:xi+lw,3] = 255
    return img


def example_test_old(strikes, disp_im=True):
    img1 = Image.fromarray(rain1(1920, strikes=strikes))
    if disp_im: img1.show()
    img1.save('rain1.jpg', format='jpeg', quality=80, optimize=True)

    img2 = Image.fromarray(rain1(1920, strikes=strikes, color=False))
    if disp_im: img2.show()
    img2.save('rain1.jpg', format='jpeg', quality=80, optimize=True)


def example_test_new(strikes, disp_im=True):
    img1 = Image.fromarray(rain2(1920, strikes=strikes))
    if disp_im: img1.show()
    img1.save('rain2.png', format='png', quality=80, optimize=True)
    
    img2 = Image.fromarray(rain2(1920, strikes=strikes, color=False))
    if disp_im: img2.show()
    img2.save('rain2.png', format='png', quality=80, optimize=True)


if __name__ == "__main__":
    # Execute only if this module is not imported into another script
    example_test_old(360)
    example_test_new(360)
    
    profile.run('example_test_old(100000, disp_im=False)')
    profile.run('example_test_new(100000, disp_im=False)')

在我的 PC 上，它的速度提高了 14.5 倍！希望这有帮助。

Answer 2

我的进步速度提高了 2 到 4 倍。

由于雨滴不会停留在图像的上半部分，所以在所有打击结束后，上半部分可以从下半部分伸展开
由于广播元组比较慢，所以改用32位格式的颜色。

def rain(width=1920, strikes=360, color=True, lw=3):
    assert not width % 16
    height = int(width * 9 / 16)
    img = np.zeros((height, width), dtype=np.uint32)
    half = height / 2
    upper_bottom = int(half) - 1
    alpha = 255 << 24

    # Paint background.
    # The upper half will always be overwritten and can be skipped.
    img[upper_bottom:] = alpha

    for i in range(strikes):
        x = round(random() * width)
        y = round(height - random() * half)
        x1 = min(x + lw, width - 1)
        if color:
            # Pack color into int. See below for details.
            rgb = int.from_bytes(randbytes(3), 'big') + alpha
        else:
            # This is how to pack color into int.
            r, b, g = 178, 255, 255
            rgb = r + (g << 8) + (b << 16) + alpha

        # Only the lower half needs to be painted in this loop.
        img[upper_bottom:y, x:x1] = rgb

    # The upper half can simply be stretched from the lower half.
    img[:upper_bottom] = img[upper_bottom]

    # Unpack uint32 to uint8 x4 without deep copying.
    img = img.view(np.uint8).reshape((height, width, 4))

    return img

注：

字节顺序被忽略。可能无法在某些平台上运行。
如果图像宽度非常大，性能会大大降低。
如果您要将 img 转换为 PIL.Image，请比较它的性能，因为它也有所改进。

由于雨水相互重叠（这使得移除for-loop变得困难）并且因为罢工的次数不多（这使得改进的空间很小），我觉得很难进一步优化。希望这就足够了。

如何有效地将随机垂直段添加到 numpy 数组中？

How do I efficiently add random vertical segments into a numpy array?

python

numpy