在python中計算體積或表面積的好算法

我想要計算3D numpy數組的體積（或表面積）。在很多情況下，體素是各向異性的，我在每個方向都有像素到釐米的轉換因子。在python中計算體積或表面積的好算法

有沒有人知道一個好的地方找到一個工具包或包做上述?? ??

現在，我有一些內部代碼，但我期待升級到更精確的工業實力。

編輯1：這是一些（差）樣品data。這比典型的球體小得多。我可以在生成它時添加更好的數據！它在（自身）腫瘤腦腫瘤。

2012-01-18 tylerthemiler

你需要定義你的意思有點更清晰。如果你的數據實際上是一個3D數組，那麼整個網格所佔據的體積是'nx * dx * ny * dy * nz * dz'，但我敢肯定你並不是這個意思......你想要體積內的等值面？ – 2012-01-18 00:56:49

我認爲你是對的。它是一個X×Y×Z的二進制數組，我想要包含在其二進制部分周邊的所有內容。它通常（但不總是）球形。 – tylerthemiler 2012-01-18 01:01:25

這聽起來很有趣，你可以發佈一些示例數據的鏈接嗎？用'pickle.dump'保存numpy數組 – wim 2012-01-18 01:14:56

一個選項是使用VTK。（我將在這裏使用tvtk python綁定...）

至少在某些情況下，獲取等值面內的區域會更精確一些。

另外，就表面積而言，tvtk.MassProperties也計算表面積。它是mass.surface_area（下面的代碼中有mass對象）。

import numpy as np 
from tvtk.api import tvtk 

def main(): 
    # Generate some data with anisotropic cells... 
    # x,y,and z will range from -2 to 2, but with a 
    # different (20, 15, and 5 for x, y, and z) number of steps 
    x,y,z = np.mgrid[-2:2:20j, -2:2:15j, -2:2:5j] 
    r = np.sqrt(x**2 + y**2 + z**2) 

    dx, dy, dz = [np.diff(it, axis=a)[0,0,0] for it, a in zip((x,y,z),(0,1,2))] 

    # Your actual data is a binary (logical) array 
    max_radius = 1.5 
    data = (r <= max_radius).astype(np.int8) 

    ideal_volume = 4.0/3 * max_radius**3 * np.pi 
    coarse_volume = data.sum() * dx * dy * dz 
    est_volume = vtk_volume(data, (dx, dy, dz), (x.min(), y.min(), z.min())) 

    coarse_error = 100 * (coarse_volume - ideal_volume)/ideal_volume 
    vtk_error = 100 * (est_volume - ideal_volume)/ideal_volume 

    print 'Ideal volume', ideal_volume 
    print 'Coarse approximation', coarse_volume, 'Error', coarse_error, '%' 
    print 'VTK approximation', est_volume, 'Error', vtk_error, '%' 

def vtk_volume(data, spacing=(1,1,1), origin=(0,0,0)): 
    data[data == 0] = -1 
    grid = tvtk.ImageData(spacing=spacing, origin=origin) 
    grid.point_data.scalars = data.T.ravel() # It wants fortran order??? 
    grid.point_data.scalars.name = 'scalars' 
    grid.dimensions = data.shape 

    iso = tvtk.ImageMarchingCubes(input=grid) 
    mass = tvtk.MassProperties(input=iso.output) 
    return mass.volume 

main()

這產生了：

Ideal volume 14.1371669412 
Coarse approximation 14.7969924812 Error 4.66731094565 % 
VTK approximation 14.1954890878 Error 0.412544796894 %

來源

2012-01-18 06:42:43

這看起來很棒！當我的老闆回來時我會試一試（我甚至沒有管理員權限，所以我不能安裝vtk :(）。 – tylerthemiler 2012-01-19 00:08:41

體素將是相當簡單，規則的多面體，不是嗎？計算每一個的數量並對它們進行求和。

來源

2012-01-18 00:52:23

主要問題是z方向太粗糙。我們採用一種貝殼方法來平均切片，但這不夠好。 – tylerthemiler 2012-01-18 00:54:52

在python中計算體積或表面積的好算法

回答

相關問題