直接“绘制”线段到numpy数组

Question

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

然后我发现here是一种很好的方法，可以使用scipy.ndimage.测量. 使用图像.imread在灰阶上。生成的numpy数组进一步简化为布尔值，因为我只对不同棒之间的连接感兴趣。然后用n图像测量. 这让我可以找出是否有路径从一边连接到另一边。这里是一个最小化的例子。在

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")

虽然这基本上是很好的montecarlo模拟和1000次迭代，以更大数量的不同棍子密度结束运行8小时或更长时间。这部分是由于这样一个事实，即所创建的图像和阵列相当大，并且绘制了数千个木条以获得更高的密度。原因是我想模拟一系列的几何图形（例如长度在500到20000像素之间），同时尽量减少像素化带来的误差。在

我想最好的办法是不使用图像数据，把它当作一个向量问题来处理，但我甚至不知道如何开始一个算法。而且许多连接也可能导致大数据阵列。在

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

你能告诉我一个有效的方法来直接生成代表木棒黑白图像的布尔型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]，一个“棒”的端点：

^{pr2}$

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>

以下是生成的图：