2

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

然后我在这里找到了一种使用 scipy.ndimage.measurements 分析图像数据的好方法。使用 ndimage.imread 以灰度方式读取图像。生成的 numpy 数组进一步减少为布尔值，因为我只对不同棒之间的连接感兴趣。然后使用 ndimage.measurements 分析生成的集群。这使我能够找出是否有从一侧连接到另一侧的路径。一个最小化的例子是here。

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 倍。在数组中创建棒坐标列表并使用一次绘图调用绘制完整列表可能会加快速度，但仍然会留下写入和读取图像的瓶颈。

您能否指出一种有效的方法来直接生成表示棒的黑白图像的布尔型 numpy 数组？也许绘制坐标列表并以某种方式将图形转换为数组？我还发现了这个有趣的讨论，其中线条被绘制到 PIL 图像上。这可能比 matplotlib 快吗？

4

1 回答 1

2

在数组中绘制线段是任何图形库的基本功能。最简单的方法可能是Bresenham 算法。该算法既简单又快速——也就是说，当以快速语言实现时。我不建议在纯 python 中实现它。该算法最简单版本的一个缺点是它没有抗锯齿。线条显示“锯齿状”。搜索“线条绘制算法”以获得具有更好抗锯齿的更高级方法。

我的眼图包中有 Bresenham 算法的Cython 实现。该函数将输入数组中的值沿从 (x0, y0) 到 (x1, y1) 的直线递增。简单地将数组值设置为 1 的修改将对该代码进行微不足道的更改。bres_segment_count

例如，

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>

这是生成的图：

棒图

于 2015-11-01T23:28:47.907 回答