用python製作3D blob的更快方法？

Question

有沒有更好的方法來製作3D密度功能？

def make_spot_3d(bright, spread, x0,y0,z0): 
    # Create x and y indices 
    x = np.linspace(-50, 50, 200) 
    y = np.linspace(-50, 50, 200) 
    z = np.linspace(-50, 50, 200) 

    X, Y, Z = np.meshgrid(x, y, z) 
    Intensity = np.uint16(bright*np.exp(-((X-x0)/spread)**2 
             -((Y-y0)/spread)**2 
             -((Z-z0)/spread)**2)) 

    return Intensity 

該函數可以生成3D numpy的陣列，其可以與編寫Mayavi

然而，當函數用於生成斑點的集羣（〜100）被繪製成如下：

Spots = np.asarray([make_spot_3d(100,2, *loc) for loc in locations]) 
cluster = np.sum(Spots, axis=0) 

得到例如：

的執行時間是約1分鐘（CPU的i5）;我打賭第可能會更快。

Answer 1

一個明顯的改進是使用廣播來評估你的強度功能在「疏」網，而不是一個完整的meshgrid，如：

X, Y, Z = np.meshgrid(x, y, z, sparse=True) 

這由約4倍的因子來減小運行時我機：

%timeit make_spot_3d(1., 1., 0, 0, 0) 
1 loops, best of 3: 1.56 s per loop 

%timeit make_spot_3d_ogrid(1., 1., 0, 0, 0) 
1 loops, best of 3: 359 ms per loop 

你可以擺脫通過在位置，利差和亮度矢量化計算參與列表理解的開銷，例如：

def make_spots(bright, spread, x0, y0, z0): 

    # Create x and y indices 
    x = np.linspace(-50, 50, 200) 
    y = np.linspace(-50, 50, 200) 
    z = np.linspace(-50, 50, 200) 

    # this will broadcast out to an (nblobs, ny, nx, nz) array 
    dx = x[None, None, :, None] - x0[:, None, None, None] 
    dy = y[None, :, None, None] - y0[:, None, None, None] 
    dz = z[None, None, None, :] - z0[:, None, None, None] 
    spread = spread[:, None, None, None] 
    bright = bright[:, None, None, None] 

    # we can save time by performing the exponentiation over 2D arrays 
    # before broadcasting out to 4D, since exp(a + b) == exp(a) * exp(b) 
    s2 = spread * spread 
    a = np.exp(-(dx * dx)/s2) 
    b = np.exp(-(dy * dy)/s2) 
    c = np.exp(-(dz * dz)/s2) 

    intensity = bright * a * b * c 

    return intensity.astype(np.uint16) 

其中bright，spread，x0，y0和z0是1D向量。這將生成一個(nblobs, ny, nx, nz)數組，然後您可以總結第一個軸。根據您生成的blob數量以及您對它們進行評估的網格大小，創建此中間數組可能會在內存方面變得非常昂貴。

另一種選擇是初始化單個(ny, nx, nz)輸出陣列，並計算就地的總和：

def sum_spots_inplace(bright, spread, x0, y0, z0): 

    # Create x and y indices 
    x = np.linspace(-50, 50, 200) 
    y = np.linspace(-50, 50, 200) 
    z = np.linspace(-50, 50, 200) 

    dx = x[None, None, :, None] - x0[:, None, None, None] 
    dy = y[None, :, None, None] - y0[:, None, None, None] 
    dz = z[None, None, None, :] - z0[:, None, None, None] 
    spread = spread[:, None, None, None] 
    bright = bright[:, None, None, None] 

    s2 = spread * spread 
    a = np.exp(-(dx * dx)/s2) 
    b = np.exp(-(dy * dy)/s2) 
    c = np.exp(-(dz * dz)/s2) 

    out = np.zeros((200, 200, 200), dtype=np.uint16) 

    for ii in xrange(bright.shape[0]): 
     out += bright[ii] * a[ii] * b[ii] * c[ii] 

    return out 

這將需要較少的存儲器，但潛在的缺點是，它需要在Python循環。

爲了讓您的相對性能的一些想法：

def sum_spots_listcomp(bright, spread, x0, y0, z0): 
    return np.sum([make_spot_3d(bright[ii], spread[ii], x0[ii], y0[ii], z0[ii]) 
        for ii in xrange(len(bright))], axis=0) 

def sum_spots_vec(bright, spread, x0, y0, z0): 
    return make_spots(bright, spread, x0, y0, z0).sum(0) 

# some fake data 
bright = np.random.rand(10) * 100 
spread = np.random.rand(10) * 100 
x0 = (np.random.rand(10) - 0.5) * 50 
y0 = (np.random.rand(10) - 0.5) * 50 
z0 = (np.random.rand(10) - 0.5) * 50 

%timeit sum_spots_listcomp(bright, spread, x0, y0, z0) 
# 1 loops, best of 3: 16.6 s per loop 

%timeit sum_spots_vec(bright, spread, x0, y0, z0) 
# 1 loops, best of 3: 1.03 s per loop 

%timeit sum_spots_inplace(bright, spread, x0, y0, z0) 
# 1 loops, best of 3: 330 ms per loop 

Answer 2

所以，你800萬（= 200 * 200 * 200）的時間在做你的任期每個操作;首先，通過計算if和mirror的八分之一，可以將其降低到100萬（如果球恰好位於網格的中心）。鏡像並不是免費的，但仍然比exp便宜得多。

另外，很有可能你應該在強度值下降到0後停止計算。使用一個小對數魔法，你可以想出一個可能比200 * 200小得多的感興趣區域* 200格。

Answer 3

既然你有一個i5處理器，並且這些點是相互獨立的，那麼實現多線程將會很不錯。由於許多Numpy操作釋放GIL，因此您不一定需要多個進程。額外的代碼可以非常簡單：

from multiprocessing.dummy import Pool 

if __name__ == '__main__': 
    wrap = lambda pos: make_spot_3d(100, 2, *pos) 
    cluster = sum(Pool().imap_unordered(wrap, positions)) 

---

更新

我在工作電腦上的一些測試中，我必須承認，上面的代碼是太天真和低效之後。在8核心上，相對於單數性能而言，加速只有〜1.5倍。

我仍然認爲多線程將是一個好主意，但成功取決於實現。