2015-11-01 143 views
0

我在python中實現的第一個項目是Monte Carlo模擬棒狀滲透。代碼不斷增長。第一部分是棒滲透的可視化。在寬度*長度的區域中,用一定長度的直杆的定義密度(棒/面積)用隨機開始座標和方向繪製。由於我經常使用gnuplot,因此我將生成的(x,y)開始和結束座標寫入文本文件,以便在之後進行gnuplot繪圖。直接「繪製」線段到numpy陣列

然後,我發現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快嗎?

回答

0

在數組中繪製線段是任何圖形庫的基本功能。最簡單的方法可能是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> 

下面是生成的情節:

sticks plot

+0

謝謝你指點我到這個方向。在我首次在純python中實現它之前,沒有使用過cython,這已經有了改進。然而,我現在正在努力讓你的代碼運行。下載完整的軟件包並保存在「cython」目錄中。運行時 – MrCyclophil

+0

最後在Linux下運行。不知道從哪裏獲得Visual Studio 2010 Express以使其在Windows下運行。 – MrCyclophil

+0

僅供參考,如果有人在Win8.1 64bit上安裝C編譯器時也有問題。對於Python 3.5,它變得更容易。只需下載並安裝免費的Visual Studio 2015社區版。 [本博客](http://blog.ionelmc.ro/2014/12/21/compiling-python-extensions-on-windows/)爲Python 2.7,3.4和3.5提供了一個很好的總結。 – MrCyclophil