利用MATLAB實(shí)現(xiàn)CT斷層圖像的三維重建

1、文章編號(hào):1004-4140(200402-0024-06利用MATLAB實(shí)現(xiàn)CT斷層圖像的三維重建曾箏1,董芳華1,陳曉2,周宏2,周建中1(1 華中科技大學(xué)水電與數(shù)字化工程學(xué)院,武漢430074; 2 總后軍需裝備研究所,北京100088摘要:醫(yī)學(xué)圖像三維重建為人體結(jié)構(gòu)提供了真實(shí)、直觀的反映,便于醫(yī)學(xué)人員對(duì)病灶的觀察及手術(shù)的進(jìn)行。但圖像三維重建編程實(shí)現(xiàn)困難,不易被非計(jì)算機(jī)專業(yè)人士所掌握。本文運(yùn)用MA TLAB6.5軟件中的圖像處理工具箱實(shí)現(xiàn)了CT斷層圖像的三維表面重建及體重建,原理簡(jiǎn)單,編程實(shí)現(xiàn)方便。在對(duì)頭部CT圖片進(jìn)行的三維表面重建及體重建實(shí)驗(yàn)中,重建速度快,顯示效果良好,便于各類非計(jì)算

2、機(jī)專業(yè)人士推廣應(yīng)用。關(guān)鍵詞:CT圖像;表面重建;體重建;MA TLAB中圖分類號(hào):TP391 文獻(xiàn)標(biāo)識(shí)碼:AThree Dimensions Reconstruction of CT Image by MATLAB ZENG Zheng1,DONG Fang-hua1,CHEN Xiao2,ZHOU Hong2,ZHOU Jian-zhong1(1. College of Hydroelectricity & Digital Engineering, HUST, Wuhan 430074, China;2. The Quartermaster Research Institute of

3、 the General Logistics Department of CPLA, Beijing 100088,ChinaAbstract: Three dimensions reconstruction of CT image gives a intuitionistic image of body. It is convenient to observation and operation. But it is difficult to those people not major in computer. In order to solve this problem, how to

4、reconstruct of CT image using image processing toolbox of MATLAB6.5 is discussed. MABLAB gives us a convenient method to reconstruct and acquires a good result. It can be mastered by anybody.Key words: CT image; surface reconstruction; volume reconstruction; MABLAB1引言隨著計(jì)算機(jī)科學(xué)技術(shù)的發(fā)展及其與其他各學(xué)科間的相互交叉滲透,特別是

5、計(jì)算機(jī)斷層掃描(CT、核磁共振掃描(MIT等技術(shù)逐漸應(yīng)用于臨床醫(yī)學(xué),使得獲取各種二維醫(yī)學(xué)圖像的能力得到了空前的提高1。如何將這些二維數(shù)據(jù)轉(zhuǎn)化為更容易被接受和處理的信息,正是醫(yī)學(xué)數(shù)據(jù)三維可視化所要解決的問(wèn)題。三維重建技術(shù)的實(shí)現(xiàn)方法包括兩種2:一種是通過(guò)幾何單元拼接擬合物體表面來(lái)描述物體的三維結(jié)構(gòu),稱為表面重建;另一種是直接將體像素以一定的顏色和透明度投影到顯示平面的方法,稱為體重建。比較而言,表面重建運(yùn)算量小,表面顯示清晰,但對(duì)邊緣檢測(cè)的要求比較高;而體重建直接基于體數(shù)據(jù)進(jìn)行顯示,避免了重建過(guò)程中所造成的偽像痕跡,但運(yùn)算量較大3。無(wú)論采用哪種重建方法,對(duì)于非計(jì)算機(jī)專業(yè)的人士,編程實(shí)現(xiàn)都是很困難的

6、,2期曾箏等:利用MA TLAB實(shí)現(xiàn)CT斷層圖像的三維重建 25 因此,本文介紹了基于MATLAB進(jìn)行三維表面重建及體重建的方法,實(shí)現(xiàn)起來(lái)靈活、方便,便于掌握。三維重建步驟如圖1所示。 圖1 三維重建流程圖為了有利于從CT圖像中準(zhǔn)確地提取出有用的信息,需要對(duì)原始圖像進(jìn)行預(yù)處理,以突出有效的圖像信息,消除或減少噪聲的干擾。(1 CT圖像格式的轉(zhuǎn)換與讀寫(xiě)。原始的CT圖像是采用DICOM3.0標(biāo)準(zhǔn)進(jìn)行存儲(chǔ)的,不能被MATLAB所識(shí)別,因此必須進(jìn)行圖像格式的轉(zhuǎn)換。在這里,使用Visual C+6.0自行開(kāi)發(fā)轉(zhuǎn)換程序,在正確讀取DICOM圖像后,通過(guò)選擇合適的窗寬、窗位,將窗寬范圍內(nèi)的值通過(guò)線性或非線性

7、變換轉(zhuǎn)換為小于256的值,將CT圖像轉(zhuǎn)換為256色BMP圖像4。使用MATLAB中的圖像讀入函數(shù)imread( ,可以讀取BMP圖像,同時(shí)可以使用圖像寫(xiě)出函數(shù)imwrite( 及圖像顯示函數(shù)image( 、imshow( 對(duì)圖像進(jìn)行寫(xiě)出和顯示5。(2 圖像增強(qiáng)。圖像增強(qiáng)就是根據(jù)某種應(yīng)用的需要,人為地突出輸入圖像中的某些信息,從而抑制或消除另一些信息的處理過(guò)程6。使輸入圖像具有更好的圖像質(zhì)量,有利于分析及識(shí)別。圖像增強(qiáng)技術(shù)主要包括直方圖修改、圖像平滑、圖像邊緣銳化及偽彩色增強(qiáng)等。分別采用以下函數(shù)實(shí)現(xiàn):灰度直方圖均衡化。均勻量化的自然圖像的灰度直方圖通常在低灰度區(qū)間上頻率較大,使得圖像中較暗區(qū)域中

8、的細(xì)節(jié)看不清楚,采用直方圖修整可使原圖像灰度集中的區(qū)域拉開(kāi)或使灰度分布均勻,從而增大反差,使圖像的細(xì)節(jié)清晰。直方圖均衡化在MATLAB 中使用histeq( 函數(shù)實(shí)現(xiàn)。灰度變換法。照片或電子方法得到的圖像,常表現(xiàn)出低對(duì)比度即整個(gè)圖像偏亮或偏暗,為此需要對(duì)圖像中的每一個(gè)像素的灰度級(jí)進(jìn)行標(biāo)度變換,擴(kuò)大圖像灰度范圍,以達(dá)到改善圖像質(zhì)量的目的。這一灰度調(diào)整過(guò)程可以用imadjust( 函數(shù)實(shí)現(xiàn)。平滑與銳化濾波。平滑技術(shù)用于平滑圖像中的噪聲,基本采用在空間域上的求平均值或中值,或在頻域上采取低通濾波。在MATLAB中,各種濾波方法都是在空間域中通過(guò)不同的卷積模板即濾波算子實(shí)現(xiàn),可用fspecial( 函

9、數(shù)創(chuàng)建預(yù)定義的濾波算子,然后用filter2( 或26 CT 理論與應(yīng)用研究 13卷conv2( 函數(shù)在實(shí)現(xiàn)卷積運(yùn)算的基礎(chǔ)上進(jìn)行濾波。中值濾波是一種基于排序統(tǒng)計(jì)理論的抑制噪聲的非線性信號(hào)處理技術(shù),其在除去圖像中的孤立點(diǎn)、線的噪聲的同時(shí),很好地保護(hù)了圖像的邊緣信息,適用于一些線性濾波器無(wú)法勝任地場(chǎng)合。在MATLAB 中使用medfilt( 函數(shù)實(shí)現(xiàn)中值濾波7。適當(dāng)運(yùn)用上述方法對(duì)原始圖像進(jìn)行處理,將使原始圖像變得較清晰,能夠較真實(shí)地反映圖像的結(jié)構(gòu)特征,便于三維重建的處理及顯示。2CT 圖像三維表面重建計(jì)算機(jī)三維表面重建是指首先運(yùn)用圖像技術(shù)從二維圖像中分割出興趣區(qū)的輪廓曲線,然后經(jīng)圖形處理,得到其三

10、維結(jié)構(gòu),從而再現(xiàn)原物體的空間結(jié)構(gòu)。因此,對(duì)于三維表面重建而言,邊界輪廓的提取尤為重要8。為了便于面部邊界的提取,先對(duì)各CT 圖片進(jìn)行顏色處理,去掉非有效區(qū),如頭發(fā)、支架等部分,并使其色素盡量減少,圖2a 為經(jīng)過(guò)格式轉(zhuǎn)換的頭部斷層圖像。在提取邊界時(shí),首先采用逐行掃描圖片的辦法,通過(guò)比較相鄰點(diǎn)的像素值,找到圖片邊界上的一個(gè)點(diǎn),作為切片邊界的起點(diǎn)。然后從邊界起點(diǎn)開(kāi)始,逐點(diǎn)判斷與之相鄰的八個(gè)點(diǎn),如果某點(diǎn)為圖片的邊界點(diǎn)則記錄下,并開(kāi)始下一步判斷,直到獲得所有的邊界點(diǎn)。圖2b 通過(guò)上述方法得到的面部邊界輪廓曲線。通過(guò)這種方法得到的輪廓曲線是以大量的點(diǎn)坐標(biāo)形式存儲(chǔ),在保證擬合精度的情況下,希望精簡(jiǎn)數(shù)據(jù)點(diǎn),以

11、減少存儲(chǔ)空間和提高計(jì)算機(jī)處理速度。因此將圖2b 中的面部邊緣輪廓曲線以中心點(diǎn)為極點(diǎn)轉(zhuǎn)換到極坐標(biāo)系中。然后,將圖像畫(huà)入以極角0,2ð為橫坐標(biāo),極徑為縱坐標(biāo)的直角坐標(biāo)系中,如圖3a 所示。采用12階傅立葉級(jí)數(shù)對(duì)該圖像進(jìn)行擬合,得到擬合后的圖像如圖3b 所示。 圖2面部邊界提取(a原始CT 圖像(b面部邊界輪廓y0 50 100 150 200 250 300 350 400 75808590 95100 105 110 115(b 擬合后的邊界曲線 (a 擬合前的邊界曲線2期 曾 箏等:利用MA TLAB 實(shí)現(xiàn)CT 斷層圖像的三維重建 27在這里,面部模型可以采用以下數(shù)學(xué)表示:=+=12

12、1,0,sin cos (2(n n i n i i i M n b n a a 式中n i a ,n i b ,為各層邊界曲線傅立葉級(jí)數(shù)擬合的系數(shù):d n f a i n i =cos (1,d n f b i n i sin (1,=。這樣,某一層斷層面部邊界曲線僅使用25個(gè)數(shù)據(jù)就可以表示,所占用的存儲(chǔ)空間很小。 下面,用這些表征面部邊界輪廓曲線的數(shù)據(jù)進(jìn)行三維表面重建。 (1 重建數(shù)據(jù)的采集。運(yùn)用上述傅立葉級(jí)數(shù)的系數(shù),求出邊界上若干個(gè)點(diǎn)x ,y 向坐標(biāo)值,并為其加上適當(dāng)?shù)膠 坐標(biāo)值。xo=0: pi/180:2*pi; % x 的值在0,2中選取yo=yo+a(i*cos(i-1*xo+b

13、(i*sin(i-1*xo; % 通過(guò)傅立葉系數(shù)求y 值,其中yo 初始值為a 0 consx=consx;yo.*cos(xo; %將x ,y 值從極坐標(biāo)系轉(zhuǎn)換到直角坐標(biāo)系 consy=consy;yo.*sin(xo;consz=consz;ones(1,length(xo*iLayer*(-4.0; %為每一切片層賦予z 坐標(biāo)值,iLayer 為層數(shù)(2 邊界輪廓曲線表面繪制。surf(consx,consy,consz; %利用surf( 函數(shù)進(jìn)行三維表面繪制。(3 設(shè)置圖像的顏色及陰影效果。colormap(gray; %利用colormap( 函數(shù)為圖像定義顏色集 shading

14、flat; %利用shading 定義顯示圖像的顏色陰影 (4 設(shè)置圖像光照效果。light('Position',-80,-262,-200,'style','infinite' %利用light( 函數(shù)為圖像設(shè)置光照效果 light('Position',-500,-0,-4500,'style','infinite' light('Position',5000,100,-300,'style','infinite'(5 設(shè)置圖像的顯示效果view

15、(-144,20; %利用view( 函數(shù)定義觀察者視角 lighting gouraud; %利用lighting 定義顯示圖像的光線陰影 axis equal; %利用axis 定義顯示圖像的軸圖4 頭部CT 圖像三維表面重建(a view(-144,20(b view(3(c view(-70,3028 CT理論與應(yīng)用研究 13卷當(dāng)函數(shù)view(AZ,EL取不同的值時(shí),可以得到圖像不同的視角,其中AZ表示圖像水平方向旋轉(zhuǎn)的角度,EL表示圖像垂直方向的高度,MATLAB6.5自行對(duì)view(2和view(3的方位進(jìn)行了定義。圖4a為view(-144,20得到的圖像;圖4b為view(3得

16、到的圖像;圖4c為view(-70,30得到的圖像。采用view( 函數(shù)可以從不同視角對(duì)三維重構(gòu)圖像進(jìn)行觀察。根據(jù)上述表面重建步驟,對(duì)頭部切片進(jìn)行表面重建,結(jié)果如圖4所示。由于采用上述的傅立葉級(jí)數(shù)擬合含有耳朵的CT圖像層時(shí),耳朵的輪廓信息無(wú)法很好地表達(dá),所以對(duì)這些層進(jìn)行處理時(shí),首先去除耳朵。從圖4中可以看出,運(yùn)用MATLAB程序在進(jìn)行CT圖像邊界輪廓提取的基礎(chǔ)上得到三維表面重建圖像。重建速度快、效果好;但是面繪制的缺點(diǎn)是信息的丟失比較大,運(yùn)算量與景物和物體形狀有關(guān)。3CT圖像三維體重建體繪制通過(guò)計(jì)算所有體素對(duì)光線的作用得到二維投影圖像,基于體繪制的三維體重建方法計(jì)算量不依賴于景物的復(fù)雜程度和物

17、體形狀的復(fù)雜程度,也不需要對(duì)切片的邊界輪廓進(jìn)行提取,其計(jì)算過(guò)程不依賴于視點(diǎn),處理三維采樣信號(hào)方便,便于顯示物體的內(nèi)部結(jié)構(gòu)。但是,三維體重建所需數(shù)據(jù)量大,運(yùn)算速度較慢。下面,介紹利用MATLAB進(jìn)行頭部CT圖像三維體重建的過(guò)程。(1 重建數(shù)據(jù)的采集。對(duì)現(xiàn)有的n幅頭部CT圖像數(shù)據(jù)進(jìn)行三維數(shù)據(jù)集D的構(gòu)造,得到的數(shù)據(jù)集D為一個(gè)x×y×n的矩陣。image1 = imread('01.bmp' %使用imread( 函數(shù)讀入現(xiàn)有的n幅圖像image2 = imread('02.bmp'imagen = imread('n.bmp'D=c

18、at(3,image1,image2,image3,imagen; %使用cat( 函數(shù)創(chuàng)建三維矩陣D(2 重建數(shù)據(jù)預(yù)處理采用上述方法構(gòu)造的三維數(shù)據(jù)集D,數(shù)據(jù)量大,在體重建中速度慢,并且可能在計(jì)算中超出內(nèi)存。因而,可以根據(jù)實(shí)際情況,對(duì)數(shù)據(jù)集D進(jìn)行預(yù)處理,減少數(shù)據(jù)量。x y z D = reducevolume(D,a b c; %使用reducevolume( 函數(shù)減少數(shù)據(jù)量,其中a,b,c為x,y和z軸數(shù)據(jù)抽取的比例,根據(jù)數(shù)據(jù)情況自行定義。D = smooth3(D; %使用smooth( 函數(shù)對(duì)數(shù)據(jù)進(jìn)行平滑處理(3 計(jì)算數(shù)據(jù)集在顯示平面累計(jì)投影。fv = isosurface(x,y,z,

19、D,isovalue; %使用isosurface( 函數(shù)計(jì)算數(shù)據(jù)集在顯示平面累計(jì)投影,isovalue根據(jù)實(shí)際情況自行定義(4 構(gòu)造三維體重建碎片p = patch(fv, FaceColor', 'yellow', 'EdgeColor', 'none' %使用patch( 函數(shù)對(duì)碎片進(jìn)行構(gòu)造,并對(duì)圖像的顏色,光線進(jìn)行定義,其中fv是第(3步中得到的。(5 設(shè)置圖像的顏色、陰影及顯示效果。2期 曾 箏等：利用 MATLAB 實(shí)現(xiàn) CT 斷層圖像的三維重建 colormap(gray; 29 %利用 colormap( 函數(shù)為圖像定義顏

20、色集 view(3; %利用 view( 函數(shù)定義觀察者視角 lighting gouraud; %利用 lighting 定義顯示圖像的光線陰影 axis equal; %利用 axis 定義顯示圖像的軸 daspect(x y z; %使用 daspect( 定義 x、y、z 軸的顯示比例 根據(jù)上述體重建步驟， 對(duì)頭部切片 進(jìn)行體重建，結(jié)果如圖 5 所示。 從圖 5 中可以看出，運(yùn)用 MATLAB 程序在構(gòu)造體重建碎片的基礎(chǔ)上，實(shí)現(xiàn)體三維重建。 重建速度較表面重建慢，但是體重建盡量保護(hù)了更多的表面信息，耳朵部分也不需要任何特 殊處理，就能得到很好的重建及顯示效果。 圖5 頭部 CT 圖像三維體重建 結(jié)論 本文介紹了運(yùn)用 MATLAB6.5 圖像處理工具箱進(jìn)行 CT 斷層圖像的三維重建的原理及實(shí) 現(xiàn)方法。并結(jié)合實(shí)際應(yīng)用，實(shí)現(xiàn)了頭部 CT 圖像的三維表面重建及體重建，編程實(shí)現(xiàn)簡(jiǎn)單， 顯示效果理想，大大提高了實(shí)驗(yàn)效率，