2012-01-18 73 views
3

3D numpy配列のボリューム(または表面積)を計算しようとしています。ボクセルは多くの場合異方性であり、各方向にピクセルからcmへの変換係数があります。Pythonでボリュームまたは表面積を計算するためのアルゴリズムが良い

誰もが、上記を行うためのツールキットまたはパッケージを見つけるのに良い場所を知っていますか?

現時点では社内コードがいくつかありますが、正確性の面でより強固なものにアップグレードすることを検討しています。

Edit1:ここにいくつかの(悪い)サンプルdataがあります。これは典型的な球よりはるかに小さい。私はそれを生成することができます私はより良いデータを追加します!それは(自己)tumorBrain.tumorにあります。

+0

あなたは、もう少し明確に何を意味するかを定義する必要があります。データが実際に3D配列であれば、グリッド全体が占める音量は 'nx * dx * ny * dy * nz * dz'ですが、あなたはそれを意味していないと確信しています。 isosurfaceの中のボリューム? –

+0

私はあなたが正しいと思う。これはX X Y Y Zバイナリ配列であり、バイナリ部分の周囲に含まれるすべてのボリュームが必要です。それは典型的に(常にではないが)球形である。 – tylerthemiler

+0

これは楽しいように聞こえますが、いくつかのサンプルデータへのリンクを投稿できますか? 'pickle.dump'でnumpy配列を保存してください – wim

答えて

4

1つのオプションは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 % 
+0

これは素晴らしいですね!私は上司が戻ってくるときに試してみます(私はここでも管理者権限を持っていないので、vtk :()をインストールすることはできません)。 – tylerthemiler

1

ボクセルはかなりシンプルで、規則的な多面体ではないでしょうか?それぞれの音量を計算し、合計します。

+0

これの主な問題は、z方向が粗すぎることです。私たちはスライス間の平均をとるための一種のシェルメソッドを実行しますが、これでは十分ではありません。 – tylerthemiler

関連する問題