戴虹
(上海第二工業(yè)大學 電子與電氣工程學院,上海 201209)
醫(yī)學圖像可視化技術通過對二維斷層圖像序列,如CT、MRI圖像進行三維重建,從而得到人體組織的三維體數(shù)據(jù),這些體數(shù)據(jù)能為醫(yī)學診斷、手術規(guī)劃和治療中的監(jiān)控等提供更完整的信息,因而成為醫(yī)學圖像領域重要的研究方向[1]。美國國家醫(yī)學圖書館的可視人項目(Visible Human Project,VHP)提供了一男一女彩色照片截面圖像以及CT和MRI斷層掃描圖像,可用于構(gòu)造各種虛擬人體三維模型,在醫(yī)學教育和臨床中發(fā)揮重要作用[2]。目前常用VC++軟件,并使用OpenGL(Open Graphics Library,開放的圖形接口)進行醫(yī)學圖像可視化編程工作,然而復雜度仍然較高。美國Kitware公司推出了VTK (Visualization Toolkit,可視化工具箱),是一種基于OpenGL的用于3D圖形學、圖像處理及三維可視化開發(fā)包[3],可嵌入到VC++中,用簡單代碼實現(xiàn)斷層圖像的三維重建。
文中對可視人項目提供的男性頭部CT圖像在matlab中編程,采用最大類間方差閾值法分別求得其皮膚和骨骼的分割閾值,將其作為等值面抽取值;利用VTK編程分別實現(xiàn)頭部皮膚和骨骼的三維表面重建。
在對頭部進行三維表面繪制之前須分別求得皮膚和骨骼等值面抽取值,即皮膚和骨骼圖像分割閾值。由于憑經(jīng)驗手工設定閾值的時間較長,故采用最大類間方差自動閾值分割算法,也稱為OTSU算法,它通過搜索圖像的灰度范圍尋找最佳分類點,即各類圖像的分割閾值,使得類間方差最大。設一幅圖像的灰度范圍是1~L級,有M-1個分割閾值t1,t2,…,tM-1,將灰度 范圍分為 C1,C2,… ,CM,共 M 類:C1←[1,t1],C2←[t1+1,t2],…,…,CM←[tM-1,L],最佳閾值{t1*,t2*,…,t*M-1}使類間方差 σ2B最大[4],即:
本課題實驗數(shù)據(jù)來自可視人項目一位男性頭部的CT序列圖像,尺寸為512×512,共245幅,為dicom格式圖像,dicom標準是醫(yī)學數(shù)字成像與通信標準,其主要目的是在各種醫(yī)療成像產(chǎn)品之間提供一致性接口,使醫(yī)療設備之間實現(xiàn)互操作?;谧畲箢愰g方差法利用Matlab2007軟件編程實現(xiàn)頭部皮膚與骨骼圖像分割的主要程序與實驗結(jié)果如下:
1)設計頭部CT圖像讀取與信息顯示界面。
利用Matlab中的dicomread與dicominfo函數(shù),分別讀入頭部CT圖像序列并顯示圖像的信
息,并用圖形界面(GUI)工具設計界面,如圖1所示(以第125幅圖像為例):
圖1 頭部CT圖片與信息顯示界面Fig.1 Interface of the head CT image and information
2)設計頭部皮膚與骨骼閾值分割程序。
根據(jù)公式(1)設計的Matlab主要程序代碼如下(設f為原圖,仍以第125幅圖像為例):
[M,N]=size (f);%圖像的尺寸為 M*N; [T,C]=imhist(f);p=T/(M*N);%p 為圖像的歸一化直方圖 L=256; max1=0; %L為最大的灰度級;設類間方差σB2最大值的初始值max1=0
for t1=1:L%對皮膚閾值t1和骨骼閾值t2進行遍歷,求得使max1最大的t1和t2(設為t10和t20)
for t2=t1:L
w0=sum(p (1:t1));w1=sum (p(t1+1:t2));w2=sum(p(t2+1:L));%計算 ω0,ω1和 ω2
i1=1:t1;i2=t1+1:t2;i3=t2+1:L;
u0=sum(i1.*p(1:t1))/w0;u1=sum(i2.*p(t1+1:t2))/w1;u2=sum(i3.*p(t2+1:L))/w2;%計算 μ0,μ1和 μ2
uT=w0*u0+w1*u1+w2*u2; %計算μT
sigB=w0*(u0-uT)^2+w1*(u1-uT)^2+w2*(u2-uT)^2;%計算
max1=sigB;t10=t1;t20=t2;
end
end
end
g1=f.*(f>=t10&f<t20);g2=f.*(f>=t20);%g1 和 g2 分別是頭部皮膚與骨骼圖像
程序運行結(jié)果為:皮膚與骨骼的最佳閾值分別為t1*=51,t2*=145,頭部皮膚(含肌肉)與骨骼圖像如圖 2(a)和(b)所示,可見該算法分割效果較好。
圖2 基于最大類間方差閾值法的頭部皮膚與骨骼分割圖像Fig.2 The segmented images of head and bone based on maximum between-class variance threshold algorithm
如果將CT序列圖像數(shù)據(jù)看成是空間區(qū)域內(nèi)采樣點的集合,則該空間區(qū)域內(nèi)具有相同值的點集將定義一個曲面,稱之為等值面,不同組織具有不同的物理屬性,因此以適當值定義的等值面可表示各組織的外表面[4]。VTK支持一系列可視化算法和高級建模技術,許多算法已經(jīng)集成到VTK中。本課題利用其所包含的vtkContourFilter類提取可視人頭部皮膚和骨骼等值面實現(xiàn)表面重建,也稱“面繪制”。等值面抽取值即圖像分割閾值,則頭部皮膚和骨骼等值面抽取值分別是t1*=51,t2*=145。
VTK采用的是管道模型(Pipeline)機制,在這種模式下,首先建立一個數(shù)據(jù)流水線并選擇相應算法來處理數(shù)據(jù),其次建立適當?shù)哪繕藞D形來演示數(shù)據(jù)[5]。采用VTK類庫實現(xiàn)頭部面繪制的管道模型如圖3所示。
圖3 基于VTK的頭部面繪制管道模型Fig.3 The surface rendering pipeline model of the head based on VTK
根據(jù)上述管道模型,基于VTK采用VC++6.0編程實現(xiàn)頭部三維重建的步驟、環(huán)境設置和主要程序代碼如下:
1)程序環(huán)境設置
到Kitware公司的VTK網(wǎng)站上下載VTK5.0開發(fā)包,安裝于本地目錄下。打開VC++6.0:
(i)選擇菜單:工具---選項,加入VTK的頭文件和庫文件目錄。
(ii)選擇:工程---設置---Link,在 Object/Library Modules中加入VTK的庫文件,如:vtkIO.lib,vtkHybrid.lib,vtkRendering.lib,vtkFiltering.lib 等。
2)建立一個vtkDICOMImageReader對象[6],讀取存儲在D盤“Head”目錄下的dicom格式
可視人頭部CT序列圖像。
#include"vtkRenderer.h"
…… //此處省略了包含的頭文件
void main()
{vtkDICOMImageReader*v16=vtkDICOMImageReader::New();
v16 ->SetDataByteOrderToLittleEndian ();//設置為小端字節(jié)序
v16 ->SetDirectoryName ("D://Head");v16 ->SetDataSpacing (1, 1, 1);//設置圖像間隔
3)建立一個vtkContourFilter對象,并設定頭部皮膚或骨骼等值面抽取值。
vtkContourFilter*skinExtractor=vtkContourFilter::New();
skinExtractor->SetInputConnection (v16->GetOutputPort());
skinExtractor->SetValue (0,51); //皮膚等值面抽取值為51
//skinExtractor->SetValue(0,143); //若提取骨骼等值面則抽取值改為143
4)利用vtkPolyDataNormals類產(chǎn)生等值面上的法線,以使繪制過程中的表面光滑著色。
vtkPolyDataNormals*skinNormals=vtkPolyDataNormals::New();
skinNormals ->SetInputConnection (skinExtractor ->GetOutputPort());
skinNormals->SetFeatureAngle(60.0);
5)由vtkPolyDataMapper類將圖像數(shù)據(jù)映射為可被計算機繪制的圖元;由vtkActor類創(chuàng)建一個代表頭部皮膚或骨骼的“角色”,設置其顏色,并與一個映射器關聯(lián)。
vtkPolyDataMapper*skinMapper= vtkPolyDataMapper::New();
skinMapper ->SetInputConnection (smoother ->GetOutputPort());
skinMapper->ScalarVisibilityOff();
vtkActor*skin=vtkActor::New();
skin->SetMapper (skinMapper); skin->GetProperty ()->SetColor(1,0.7,0.4);//設置顏色
6)利用vtkRenderer類創(chuàng)建一個繪制器,添加被繪制的“角色”;用vtkRenderWindow在繪制窗口中顯示;通過vtkRenderWindowInteractor實現(xiàn)繪制結(jié)果與用戶的交互。
vtkRenderer*aRenderer=vtkRenderer::New();
aRenderer->AddActor(skin); aRenderer->SetBackground(1,1,1);//設置背景顏色為白色vtkRenderWindow*renWin=vtkRenderWindow::New();renWin->AddRenderer (aRenderer); renWin->SetSize(640, 480);//設置繪制窗口大小為 640*480
vtkRenderWindowInteractor *iren =vtkRenderWindowInteractor::New();
iren->SetRenderWindow(renWin);}
利用VTK進行頭部皮膚與骨骼表面重建結(jié)果分別如圖3(a)和圖 3(b)所示。 可見基于 vtkContourfilter類的面繪制算法產(chǎn)生了清晰的重建效果。在程序中采用vtkPolyDataNormals類產(chǎn)生等值面上的法線其作用是對重建表面進行平滑處理,從而使結(jié)果更為逼真;利用鼠標還能直接拖動三維模型進行旋轉(zhuǎn),增強了系統(tǒng)的交互能力;程序代碼僅20~30句,故編程十分簡單;程序運行時間為:8.40秒,因此三維繪制速度較快。
圖4 基于VTK的可視人頭部皮膚與骨骼三維表面重建Fig.4 Three dimensional surface reconstruction of the visible human’s skin and bone of head based on VTK
提出一種基于VTK[7]的可視人頭部皮膚與骨骼三維可視化方法:1)采用最大類間方差閾值法獲得皮膚與骨骼等值面抽取值。2)基于VTK編程實現(xiàn)頭部的三維表面重建。實驗結(jié)果表明,該方法重建結(jié)果比較理想,編程簡單且運行速度很快,也可推廣至對人體所有器官進行三維可視化工作。
[1]曹毅.基于體繪制技術的三維醫(yī)學圖像可視化平臺研究[D].成都:四川大學,2007.
[2]Michael J.Ackerman.The Visible Human Project Getting the Data [EB/OL].[2013-06-11].http://www.nlm.nih.gov/research/visible/
[3]WANG Hong-jian.3D Medical CT Images Reconstruction Based on VTK and Visual C++ [C]//The 3rd International Conference on Bioinformatics and Biomedical Engineering,Beijing,2009:1-4.
[4]謝術富.一種醫(yī)學圖像的3-D表面重建方法 [J].計算機工程與應用,2006,42(6):230-232.XIE Shu-fu.A Method of 3-D surface reconstruction of medical images[J].Computer Engineering and Application,2006,42(6):230-232.
[5]WANG Yan-an,WAN Wang-gen,WANG Rui.An improved interpolation algorithm using nearest neighbor from VTK[C]//2010 International Conference On Image Processing,Hong Kong,2010:1062-1065.
[6]賀銳鋼,杜紅.基于VTK的醫(yī)學圖像三維重建的研究與實現(xiàn)[J].電腦知識與技術,2013,9(1):165-167.HE Rui-gang,DU Hong.Research and implementation of 3D medical image reconstruction based on VTK[J].Computer Knowledge and Technology,2013,9(1):165-167.
[7]李永毅,鐘雪麗.基于掛片協(xié)議的醫(yī)學影像顯示系統(tǒng)架構(gòu)[J].電子科技,2012(1):81-84.LI Yong-yi,ZHONG Xue-li.Medical imaging display system architecture based on hanging protocols[J].Electronic Science and Technology,2012(1):81-84.