python - 直接 "plot"线段到numpy数组

Question

我在 python 中实现的第一个项目之一是对棒渗流进行蒙特卡罗模拟。代码不断增长。第一部分是棍子渗滤的可视化。在宽度*长度的区域中，使用随机起始坐标和方向绘制具有一定长度的直棒的定义密度(棒/面积)。因为我经常使用 gnuplot，所以我将生成的 (x, y) 开始和结束坐标写入一个文本文件，以便之后用 gnuplot 绘制它们。

然后我找到了here使用 scipy.ndimage.measurements 分析图像数据的好方法。在灰度中使用 ndimage.imread 读取图像。生成的 numpy 数组进一步简化为 bool 值，因为我只对不同棒之间的连接感兴趣。然后使用 ndimage.measurements 分析生成的簇。这使我能够确定是否存在从一侧连接到另一侧的路径。此处为最小化示例。

import random
import math
from scipy.ndimage import measurements
from scipy.ndimage import imread
import numpy as np
import matplotlib.pyplot as plt

#dimensions of plot
width = 10
length = 8
stick_length = 1
fig = plt.figure(frameon=False)
ax = fig.add_axes([0, 0, 1, 1])
fig.set_figwidth(width)
fig.set_figheight(length)
ax.axis('off')

file = open("coordinates.txt", "w")

for i in range (300):
    # randomly create (x,y) start coordinates in channel and direction
    xstart = width * random.random()        # xstart = 18
    ystart = length * random.random()        # ystart = 2
    # randomly generate direction of stick from start coordinates and convert from GRAD in RAD
    dirgrad = 360 * random.random()
    dirrad = math.radians(dirgrad)
    # calculate (x,y) end coordinates
    xend = xstart + (math.cos(dirrad) * stick_length)
    yend = ystart + (math.sin(dirrad) * stick_length)
    # write start and end coordinates into text file for gnuplot plotting
    file.write(str(i) + ":\t" + str(xstart) + "\t" + str(ystart) + "\t" + str(dirgrad) + ":\t" + str(xend) + "\t" + str(yend) + "\n")
    file.write(str(i) + ":\t" + str(xend) + "\t" + str(yend) + "\n\n")
    # or plot directly with matplotlib
    ax.plot([xstart,xend],[ystart,yend],"black", lw=1)
fig.savefig("testimage.png", dpi=100)

# now read just saved image and do analysis with scipy.ndimage
fig1, ax1 = plt.subplots(1,1)
img_input = imread("testimage.png", flatten = True) # read image to np.ndarray in grey scales
img_bw = img_input < 255 # convert grey scales to b/w (boolean)
labeled_array, num_clusters = measurements.label(img_bw) #labeled_array: labeled clusters in array, num_clusters: number of clusters
area = measurements.sum(img_bw, labeled_array, index=np.arange(labeled_array.max() + 1)) # area of each cluster
areaImg = area[labeled_array] # label each cluster with labelnumber=area
cax = ax1.imshow(areaImg, origin='upper', interpolation='nearest', cmap = 'rainbow')
cbar = fig1.colorbar(cax)
fig1.savefig("testimage_analyzed.png")

虽然这主要适用于大量不同棒密度的 1000 次迭代的蒙特卡罗模拟最终运行 8 小时或更长时间。这部分是由于这样一个事实，即创建的图像和阵列非常大，并且为了更高的密度绘制了数千根棒。原因是我想模拟一定范围的几何形状(例如，长度在 500 到 20000 像素之间)，同时最大限度地减少像素化造成的误差。

我想最好的方法是不使用图像数据并将其视为矢量问题，但我什至不知道如何开始算法。许多连接也可能导致大型数据阵列。

继续使用上述方法，很明显将数据写入文件并重新读取它并不是很有效。因此，我正在寻找加快速度的方法。作为第一步，我使用 matplotlib 来创建图像，但是至少在使用单独的绘图调用绘制每根棍子时，对于更多的棍子，这会慢 10 倍。在数组中创建棒坐标列表并使用一次绘图调用绘制完整列表可能会加快速度，但仍然存在写入和读取图像的瓶颈。

你能告诉我一个有效的方法来直接生成代表木棍黑白图像的 bool 类型 numpy 数组吗？也许绘制坐标列表并以某种方式将图形转换为数组？我也发现这很有趣 discussion其中线条被绘制成 PIL 图像。这可能比 matplotlib 更快吗？

Answer 1

在数组中绘制线段是任何图形库的基本功能。最简单的方法可能是 Bresenham's algorithm .该算法简单而快速——当用快速语言实现时，就是这样。我不建议在纯 python 中实现它。该算法最简单版本的一个缺点是它不能抗锯齿。行显示 "jaggies" .搜索“画线算法”以获得具有更好抗锯齿的更高级方法。

我有一个 Cython implementation of Bresenham's algorithm在我的 eyediagram package .函数 bres_segment_count 沿着从 (x0, y0) 到 (x1, y1) 的直线递增输入数组中的值。将数组值简单地设置为 1 的修改对该代码来说是微不足道的更改。

例如，

In [21]: dim = 250

In [22]: num_sticks = 300

sticks 的每一行包含 [x0, y0, x1, y1]，“stick”的端点:

In [23]: sticks = np.random.randint(0, dim, size=(num_sticks, 4)).astype(np.int32)

In [24]: img = np.zeros((dim, dim), dtype=np.int32)

bres_segments_count 使用 Bresenham 算法绘制每根棍子。请注意，不是简单地将行中的值设置为 1，img 中的值沿行递增。

In [25]: from eyediagram._brescount import bres_segments_count

In [26]: bres_segments_count(sticks, img)

In [27]: plt.imshow(img, interpolation='nearest', cmap=cm.hot)
Out[27]: <matplotlib.image.AxesImage at 0x10f94b110>

这是生成的图: