直接“绘制”线段到numpy阵列

Question

我在python中实现的第一个项目是Monte Carlo模拟棒状渗透。代码不断增长。第一部分是棒渗透的可视化。在宽度*长度的区域中，用一定长度的直杆的定义密度（棒/面积）用随机开始坐标和方向绘制。由于我经常使用gnuplot，因此我将生成的（x，y）开始和结束坐标写入文本文件，以便在之后进行gnuplot绘图。

然后，我发现here一个很好的方式来分析图像数据使用scipy.ndimage.measurements。通过使用灰度图中的ndimage.imread来读取图像。生成的numpy数组进一步减少为布尔值，因为我只对不同的棒之间的连接感兴趣。然后用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") 

虽然这个工程主要就好了Monte Carlo模拟与1000次迭代不同棒密度较大数量的最终运行8小时以上。这部分是由于这样一个事实，即创建的图像&阵列非常大，并且为更高的密度绘制了数千个棒。原因是我想模拟一系列几何图形（例如长度在500和20000像素之间），同时最大限度地减少由于像素化造成的误差。

我想最好的方法是不使用图像数据并将其视为矢量问题，但我不知道如何启动算法。许多连接也可能导致大数据数组。保持上述方法，很明显，将数据写入文件并重新读取文件并不是非常有效。因此，我正在寻找方法加快这一进程。作为第一步，我使用matplotlib来创建图像，但至少当用单独的绘图函数绘制每个棒时，对于大量的棒来说，这比速度慢10倍。在数组中创建棍子坐标列表并通过一次调用绘制完整列表可能会加快速度，但仍会留下写入和读取图像的瓶颈。

你能指点我一个有效的方法来直接生成布尔类型numpy数组，代表棒的黑白图像吗？也许绘制坐标列表并以某种方式将数字转换为数组？我还发现这个有趣的discussion，其中线绘制为PIL图像。这可能比matplotlib快吗？

Answer 1

在数组中绘制线段是任何图形库的基本功能。最简单的方法可能是Bresenham's algorithm。该算法简单快捷 - 以快速语言实现时，即。我不会推荐在纯Python中实现它。最简单版本的算法的缺点是它不是反锯齿。线条显示"jaggies"。搜索“线条绘制算法”以获得更好的抗锯齿更高级的方法。我的eyediagram package中有一个Cython implementation of Bresenham's algorithm。函数bres_segment_count将输入数组中的值沿着从（x0，y0）到（x1，y1）的直线递增。只需将数组值设置为1的修改对该代码而言是微不足道的变化。

例如，

In [21]: dim = 250 

In [22]: num_sticks = 300 

的sticks每一行包含[X0，Y0，X1，Y1]中， “粘” 的端点：

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> 

下面是生成的情节：