一、問題描述
最近在開發(fā)過程中遇到了這樣的問題：
在醫(yī)學(xué)圖像開發(fā)過程中，我們將醫(yī)學(xué)圖像通過深度學(xué)習(xí)算法進(jìn)行分割，現(xiàn)在想要通過這一套二維圖像進(jìn)行三維重構(gòu)。
以下是分割結(jié)果：
圖一：前列腺核磁圖像分割結(jié)果 圖一：前列腺核磁圖像分割結(jié)果 圖一：前列腺核磁圖像分割結(jié)果
以下是讀取的遮罩mask：
圖二：圖像分割遮罩
圖二：圖像分割遮罩
圖二：圖像分割遮罩
如何將這些二維圖像進(jìn)行三維重建，是個(gè)棘手問題，筆者通過vtk進(jìn)行建模操作。
二、解決方案
0. 寫在前面
醫(yī)學(xué)圖像的三維重建本身就是熱點(diǎn)技術(shù)，這項(xiàng)技術(shù)也并非新鮮技術(shù)，筆者調(diào)研多份前者的博客與其余資料，整理出了自己的解決方案，旨在與大家共同交流，如果您有更好的建模方案，歡迎隨時(shí)與我交流！
1. 準(zhǔn)備工作
進(jìn)行醫(yī)學(xué)圖像的三維重建，首先需要提供清晰可見的輪廓與遮罩。（如圖二所示）
所用到的庫：
- vtk，您可以通過
pip install vtk直接安裝
2. 文件結(jié)構(gòu)
-
mask文件夾 （用于存放分割結(jié)果遮罩，圖片名為
mask_0.png,mask_1.png,mask_2.png等20張圖片） - vtk_gaussian.py （python腳本，用于執(zhí)行并進(jìn)行三維重建）
如圖所示：
圖三：項(xiàng)目采用的文件結(jié)構(gòu) 圖三：項(xiàng)目采用的文件結(jié)構(gòu) 圖三：項(xiàng)目采用的文件結(jié)構(gòu)
3. 代碼講解
3.0 完整代碼
我知道有的朋友比較急，這里先給出完整代碼：
import vtk
# 定義渲染窗口、交互模式
aRender = vtk.vtkRenderer()
Renwin = vtk.vtkRenderWindow()
Renwin.AddRenderer(aRender)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(Renwin)
# 定義個(gè)圖片讀取接口
# 讀取PNG圖片就換成PNG_Reader = vtk.vtkPNGReader()
PNG_Reader = vtk.vtkPNGReader()
PNG_Reader.SetNumberOfScalarComponents(1)
PNG_Reader.SetFileDimensionality(2) # 說明圖像是三維的
# 定義圖像大小，本行表示圖像大小為（512*512*240）
PNG_Reader.SetDataExtent(0, 256, 0, 256, 0, 19)
# 設(shè)置圖像的存放位置
name_prefix = ['mask/mask_']
PNG_Reader.SetFilePrefix(name_prefix[0])
# 設(shè)置圖像前綴名字
# 表示圖像前綴為數(shù)字（如：0.jpg）
PNG_Reader.SetFilePattern("%s%d.png")
PNG_Reader.Update()
PNG_Reader.SetDataByteOrderToLittleEndian()
spacing = [1.0, 1.0, 2.5] # x, y 方向上的間距為 2 像素，z 方向上的間距為 2.5 像素
PNG_Reader.GetOutput().SetSpacing(spacing)
# 高斯平滑
gauss = vtk.vtkImageGaussianSmooth()
gauss.SetInputConnection(PNG_Reader.GetOutputPort())
gauss.SetStandardDeviations(1.0, 1.0, 1.0)
gauss.SetRadiusFactors(1.0, 1.0, 1.0)
gauss.Update()
# 計(jì)算輪廓的方法
contour = vtk.vtkMarchingCubes()
gauss.GetOutput().SetSpacing(spacing)
contour.SetInputConnection(gauss.GetOutputPort())
contour.ComputeNormalsOn()
contour.SetValue(0, 100)
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(contour.GetOutputPort())
mapper.ScalarVisibilityOff()
actor = vtk.vtkActor()
actor.SetMapper(mapper)
renderer = vtk.vtkRenderer()
renderer.SetBackground([1.0, 1.0, 1.0])
renderer.AddActor(actor)
window = vtk.vtkRenderWindow()
window.SetSize(512, 512)
window.AddRenderer(renderer)
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(window)
# 開始顯示
if __name__ == '__main__':
window.Render()
interactor.Initialize()
interactor.Start()
3.1 定義渲染窗口、交互模式
import vtk
# 定義渲染窗口、交互模式
aRender = vtk.vtkRenderer()
Renwin = vtk.vtkRenderWindow()
Renwin.AddRenderer(aRender)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(Renwin)
本人的所有三維重建腳本中幾乎都包含這一塊內(nèi)容，同時(shí)這也是進(jìn)行vtk交互窗口初始化的部分，更多信息您可以查閱vtk官方文檔或者其他技術(shù)博客。
3.2 讀取二維圖像序列 - 定義讀取接口
要將mask_0.png到mask_19.png的圖像全部讀取，需要先定義個(gè)圖片讀取接口 (vtk.vtkXxxReader)。
筆者的圖像為.png格式，因此使用vtk.vtkPNGReader() 進(jìn)行圖像讀取。
PNG_Reader = vtk.vtkPNGReader()
PNG_Reader.SetNumberOfScalarComponents(1)
PNG_Reader.SetFileDimensionality(2) # 說明圖像是二維的
在這里，可以根據(jù)不同的圖片格式選取不同的vtkReader，如果是 .jpg 格式的圖像，可以有如下更改：
JPG_Reader = vtk.vtkJPEGReader()
# your code here...
3.3 讀取二維圖像序列 - 前置設(shè)置 & 圖像讀取
# 定義圖像大小，本人表示圖像大小為（256*256）
# 后兩個(gè)參數(shù)是圖片的數(shù)目，本人這里所用的圖像共20張，所以就輸入0, 19
# 在后續(xù)讀取的時(shí)候，會(huì)根據(jù)這個(gè)序列進(jìn)行讀取
PNG_Reader.SetDataExtent(0, 256, 0, 256, 0, 19)
# 設(shè)置圖像的存放位置
name_prefix = ['mask/mask_']
PNG_Reader.SetFilePrefix(name_prefix[0])
# 表示圖像前綴為數(shù)字（如：0.jpg）
PNG_Reader.SetFilePattern("%s%d.png")
PNG_Reader.Update()
PNG_Reader.SetDataByteOrderToLittleEndian()
這段代碼是進(jìn)行圖像讀取的一些前置設(shè)置。SetDataExtent()函數(shù)的參數(shù)設(shè)置會(huì)對(duì)后續(xù)圖像的處理會(huì)有一定影響，請(qǐng)正確填寫您所使用的圖片的大小和數(shù)目！
SetFilePrefix() 函數(shù)會(huì)根據(jù)傳入的字符串進(jìn)行鎖定。在這里一定要特別注意：
- 本人的圖像存放在
mask文件夾下，每張圖片的名字為：mask_0.png,mask_1.png等 - 在這里設(shè)置Prefix時(shí)，就要輸入
mask/mask_ - 后面的
SetFilePattern()函數(shù)會(huì)自動(dòng)讀取數(shù)字，因?yàn)榍熬Y已經(jīng)設(shè)置好，不需要再在此處進(jìn)行一些正則運(yùn)算符操作
3.4 圖像數(shù)據(jù)量少，三維建模的結(jié)果很扁平
解決方案：您可以增大圖像之間的間距來解決這個(gè)問題，緊接著上面的代碼：
PNG_Reader.SetDataByteOrderToLittleEndian()
spacing = [1.0, 1.0, 2.5] # x, y 方向上的間距為 2 像素，z 方向上的間距為 2.5 像素
PNG_Reader.GetOutput().SetSpacing(spacing)
在讀取好圖片之后，可以設(shè)置圖像之間的 spacing 這個(gè)列表，分別代表x, y, z 三個(gè)維度的間距，其中我將z維度的間距增大為2.5， 這樣的操作在后面的建模中有顯著效果。具體內(nèi)容如下：
圖四：兩種不同間距的建模結(jié)果，左圖為 z = 1.0 ，右圖為 z = 2.5 圖四：兩種不同間距的建模結(jié)果，左圖為z=1.0，右圖為z=2.5 圖四：兩種不同間距的建模結(jié)果，左圖為z=1.0，右圖為z=2.5
為了使建模結(jié)果更接近與器官的形狀，我建議設(shè)置好二維圖像之間的間距。
3.5 三維重建結(jié)果的平滑—高斯平滑
原本建模的結(jié)果“層次分明”，并不是特別美觀。筆者采用高斯平滑的方案對(duì)圖像進(jìn)行平滑。
如果您有更好的平滑方案，歡迎您與我交流！
# 高斯平滑
gauss = vtk.vtkImageGaussianSmooth()
gauss.SetInputConnection(PNG_Reader.GetOutputPort())
gauss.SetStandardDeviations(1.0, 1.0, 1.0)
gauss.SetRadiusFactors(1.0, 1.0, 1.0)
gauss.Update()
圖五：兩種不同間距的平滑結(jié)果，左圖為未平滑，右圖為使用高斯平滑 圖五：兩種不同間距的平滑結(jié)果，左圖為未平滑，右圖為使用高斯平滑 圖五：兩種不同間距的平滑結(jié)果，左圖為未平滑，右圖為使用高斯平滑
3.6 計(jì)算輪廓與邊緣提取
在進(jìn)行高斯平滑之后，進(jìn)行邊緣提取。
# 計(jì)算輪廓的方法
contour = vtk.vtkMarchingCubes()
gauss.GetOutput().SetSpacing(spacing)
contour.SetInputConnection(gauss.GetOutputPort())
contour.ComputeNormalsOn()
contour.SetValue(0, 100)
3.7 管道操作與可視化展示
后面這段代碼是vtk的顯示部分，筆者一般不去動(dòng)它，也是每一份腳本中的固有內(nèi)容，如果您對(duì)該部分感興趣，您應(yīng)該查閱vtk官方文檔。
mapper = vtk.vtkPolyDataMapper()
mapper.SetInputConnection(contour.GetOutputPort())
mapper.ScalarVisibilityOff()
actor = vtk.vtkActor()
actor.SetMapper(mapper)
renderer = vtk.vtkRenderer()
renderer.SetBackground([1.0, 1.0, 1.0])
renderer.AddActor(actor)
window = vtk.vtkRenderWindow()
window.SetSize(512, 512)
window.AddRenderer(renderer)
interactor = vtk.vtkRenderWindowInteractor()
interactor.SetRenderWindow(window)
# 開始顯示
if __name__ == '__main__':
window.Render()
interactor.Initialize()
interactor.Start()
三、一些vtk庫的使用經(jīng)驗(yàn)分享
1. python vtk
個(gè)人感覺python vtk的開發(fā)并沒有pythonic風(fēng)格，開發(fā)者有一些將cpp的開發(fā)思路帶入python庫/接口的設(shè)計(jì)，讓我寫起來味如嚼蠟，從上面的代碼亦可以看出，具有濃厚的cpp風(fēng)格。
不過代碼能跑就行，不得不說vtk仍然是很強(qiáng)大的三維重建工具！
2. 數(shù)據(jù)傳入的兩種方式
2.1 .GetOutputPort() 和 .SetInputConnection()
您會(huì)看見諸如 contour.SetInputConnection(gauss.GetOutputPort()) 的句子。
這兩個(gè)函數(shù)一般是成對(duì)出現(xiàn)，上下傳遞的。
2.1 .GetOutput() 和 .SetInputData()
筆者在參閱其他博主的博客時(shí)，同樣看見這樣的寫法，例如:
contour.SetInputData(gauss.GetOutput())
這兩個(gè)函數(shù)一般是成對(duì)出現(xiàn)的，進(jìn)行上下傳遞處理結(jié)果。
3. vtkImageData 和 vtkPolyData
在開發(fā)過程中，您可能會(huì)遇到不少報(bào)錯(cuò)，其中肯定會(huì)有vtk數(shù)據(jù)類型報(bào)錯(cuò)的問題。我整理了一份表格：
| Name | Input Type | Return Type | Variable |
|---|---|---|---|
| vtk.vtkPNGReader() | ? | vtkImageData | PNG_Reader |
| vtk.vtkImageGaussianSmooth() | vtkImageData | vtkImageData | gauss |
| vtk.vtkMarchingCubes() | vtkImageData | vtkPolyData | contour |
| vtk.vtkPolyDataNormals() | vtkImageData | vtkPolyData | normfilter |
希望能夠幫助您解決開發(fā)過程中的一些疑惑。
后記
醫(yī)學(xué)圖像的三維重建工作部分博客較少，筆者希望提供一些星星之火，大家共同進(jìn)步！
筆者采用的重建代碼已經(jīng)打包至百度云盤，您可以通過下面的鏈接下載：
下載鏈接文章來源：http://www.zghlxwxcb.cn/news/detail-428536.html
提取碼：jpt3文章來源地址http://www.zghlxwxcb.cn/news/detail-428536.html
到了這里，關(guān)于【Python VTK】讀取二維序列醫(yī)學(xué)圖像分割結(jié)果并進(jìn)行三維重建的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
