Last update 09-Aug-2012
3次元(3D)画像データをどのように観るか?
2次元(2D)データの場合は画面に映し出すだけでよく、いろいろなソフトウェアが利用出来ます。特にImageJはいろいろなフォーマットの画像データファイルに対応しており、広く用いられています。ImageJの特徴は、本体がJavaで書かれておりJavaで比較的容易にプラグインが作成できることです。実際、いろいろなプラグインが利用可能です。しかしImageJはもともと2Dイメージ用に開発されたものであり、3Dイメージデータの可視化や解析にはあまり向いていません。
3D画像データの扱いは、多くの人たち、とりわけCT画像やMRI画像を扱う医療関係者にとって大きな問題であり、この分野でソフトウェアの開発が行われて来ています。汎用性のある基礎的なプログラムを集めたライブラリとしては、OpenCV (Computer Vision)、VTK (Visualization Toolkit)、ITK (Insight Segmentation and Registration Toolkit)などがあります。ここではまず、3D画像データの可視化を主な目的としているVTKの利用を目指します。
VTKは、画像処理用のいろいろなツールを集めた巨大なライブラリです。それらをプログラムから呼んで利用します。VTK自体はC/C++で書かれていますが、呼ぶ側のプログラム言語としては、C/C++、Python、Tclの利用が可能です。手軽さを重視すると、現時点では呼ぶプログラム側の言語にPythonを使用するのが妥当と思われます。
VTKは膨大なライブラリであり、その全体像を把握することは容易ではありません。利用するにはexample programを走らせコードを理解し、自らの目的にあうパーツを集めてプログラムを作成していく必要があるようです。そのためにはかなりの経験(時間の浪費)が必要と予想されます。
とりあえずこの説明では、VTKを利用するためのシステムをWindowsで準備し、簡単なプログラムを走らせてみることにします。3D可視化のもう少し理論的なことは別の機会に。
ある程度のメモリ量とグラフィック機能が必要なようです。
MacBook Air 11インチ(2 GBメモリ、NVIDIA GeForce 320M)でもそこそこ快適。
DELL Optiplex 760 (4 GBメモリ、Intel Core 2 Duo, E7400, on-chip graphics)では、かなりしんどい。
Windows機でグラフィックカードを試してみたところ、
NVIDIA GeForce 210ではまずまずの快適さ、
NVIDIA GeForce GT440を期待して入れたが、速さは一応OKだが、3D画像を回転した時に残像が残る、
ATI Radeon HD 2600 XTあるいはNVIDIA Quadro 600では問題のない快適、
となりました。
VTKのグラフィックはOpenGLを利用するので、同じNVIDIAの製品でもDirect X向きではなくOpenGL向きのボードの方が適しているのではないかと思います(NVIDIAではGeForceよりQuadro系列)。
CMakeのダウンロードサイトより, Windows用のインストールファイル(cmake-2.8.8-win32-x86.exe)をダウンロード。
インストールするにはファイルをクリック。オプションはデフォールトのままでOK。
cmakeを起動するには、Windowsのスタート→すべてのプログラム→CMake 2.8→CMake (cmake-gui)をクリックします。
PythonにはVersion 2の系列とVersion 3の系列がありますが、Version 3に後方互換性はありません。現在のVTKはVersion 3に未対応であるため、Python version 2.7.3を使用。
Pythonの公式サイトのダウンロードサイトより、Windows用のインストーラをダウンロード。32ビットの64ビットのどちらを選ぶかは、使用しているWindowsが32ビット版か64ビット版であるかに依存。もしWindowsが64ビット版であれば64ビットを選ぶ方がよい。それ以外は32ビットを選ぶ。
ダウンロードした.msiファイルをクリック。オプションはデフォールトでOK。PATHを通しておく必要があります。
numpyとscipyは、コンパイルしてライブラリを作成する方法もりますが、ここでは簡便に作成されたライブラリをダウンロードします。
http://www.lfd.uci.edu/~gohlke/pythonlibs/より、numpy、scipyのインストールファイルをダウンロード。
Pythonが32ビット版の場合はnumpy-MKL-1.6.1-win32-py2.7.exe、Pythonが64ビット版の場合はnumpy-MKL-1.6.1-amd64-py2.7.exe。同じようにnumpyをダウンロード。numpy、scipyの順序でインストール。いずれも順番にクリックしていくだけでインストール出来るはずです。
これらが問題なくインストール出来ているかどうかの確認は、
MacOSもしくはLinuxの場合は、cmakeを利用してMakefileを作成し、makeでライブラリをbuildします。そして管理者権限を有する状態で、make installを行います。MacOSの場合は、cmakeで.xcodeprojを作成し、XCodeを用いてbuildすることも出来ます。cmakeでの設定条件は、基本的にはWindowsの場合と同じです。
レーザー走査顕微鏡で得たデータとしては、モノクロ16ビットデータが基本的なものだと思われます。モノクロ16ビットデータのTIFFファイルをVTKで読むことは出来ますが、ファイル内に複数の画像が含まれる形式(multi-page format)はプログラムのバグのために読めない様です(下記参照)。
このバグを回避するためには、multi-pageファイルを単ページファイルに変換し、それを取り込むクラスvtkTIFFReaderがあります。ファイルの変換は、ImageJ(またはImageJ64)を用いて行うのが便利です。変換の時に、Voxelサイズ(x, y, z)値を正確に設定しておくことが大切です(特にz)。
vtkTIFFReaderはTIFFファイルを読み込むためのクラスです。単一画像のファイルは問題なく読めます。モノクロ16ビットでmulti-pageのファイルはエラーで中断してしまいます。ソースファイルを調べたところ、multi-pageファイルを読む場合は、ReadVolumeという名前の関数が用いられますが、その関数をよぶ前にOutputExtent変数が設定されていないためにエラーが起きるという事がわかりました。その応急的な修正は下記の通りです。この修正が他に影響を及ぼす可能性については検討していません。もともとこの応急的修正はvtk-5.8用のものですが、vtk-5.10でも働きます。
別の方法として、vtkTIFFReader.cxxの570行目あたりの
self->ReadVolume( outPtr );
の前に、その20行程下で用いられている関数
vtkTIFFReaderUpdate2(self, reader, outPtr2, outExtent);
を加えるだけでもよいのかもしれませんが、十分テストしていません。
バグ修正後のライブラリで動作します。
作業フォルダに、このpythonプログラムファイルとデータファイルが置かれていることを想定しています。コマンドプロンプトを起動し、cdにより作業フォルダに移動し、
python vol05.py
でプログラムを動かします。
#!/usr/local/bin/python
import os
import vtk
from vtk.util.colors import *
def readFile(filename):
print filename
reader = vtk.vtkTIFFReader()
reader.SetFileName(filename)
reader.GetOutput().SetOrigin(0.0,0.0,0.0)
reader.SetDataExtent(0,511,0,511,0,40)
reader.SetDataSpacing(1.0, 1.0, 1.0)
reader.Update()
print reader.GetOutput().GetDimensions()
print reader.GetOutput().GetScalarSize()
return reader.GetOutput()
def cellSurface(cell):
med = vtk.vtkImageMedian3D()
med.SetInput(cell)
med.SetKernelSize(3,3,2)
med.ReleaseDataFlagOff()
surface = vtk.vtkMarchingCubes()
surface.SetInputConnection(med.GetOutputPort())
surface.ComputeGradientsOn()
surface.ComputeScalarsOff()
surface.SetValue(0, 380)
return surface.GetOutput()
def main():
neuron = readFile('./Neuron2.tif')
glia = readFile('./Microglia2.tif')
neuronSurface = cellSurface(neuron)
gliaSurface = cellSurface(glia)
# mapper
mapper1 = vtk.vtkDataSetMapper()
mapper1.SetInput(neuronSurface)
mapper2 = vtk.vtkDataSetMapper()
mapper2.SetInput(gliaSurface)
# actor
actor1 = vtk.vtkActor()
actor1.SetMapper(mapper1)
actor1.GetProperty().SetColor(green)
actor2 = vtk.vtkActor()
actor2.SetMapper(mapper2)
actor2.GetProperty().SetColor(red)
# renderer
ren = vtk.vtkRenderer()
ren.AddActor(actor1)
ren.AddActor(actor2)
ren.SetBackground(0.3,0.4,0.2)
# render window and interactor
renWin = vtk.vtkRenderWindow()
renWin.AddRenderer(ren)
renWin.SetSize(800,800)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
iren.Initialize()
iren.Start()
if __name__ == '__main__':
main()
mpegファイルに出力する事が可能です。しかし著作権の問題を避けるために、VTKのパッケージにはmepgのエンコーダーは含まれていないので、mpeg2encoderを別にダウンロードし、ライブラリlibMPEG2Encodeを作らなくてはなりません。