王興家 董利娜 李傳富 范 亞 馮煥清*（中國(guó)科學(xué)技術(shù)大學(xué)電子科學(xué)與技術(shù)系，合肥 3007）

2（安徽中醫(yī)學(xué)院第一附屬醫(yī)院影像中心，合肥 230031）

引言

隨著多層螺旋CT（multi-slice spiral CT，MSCT）技術(shù)的快速發(fā)展，64層乃至更高層數(shù)（如128、256）的MSCT裝備進(jìn)了許多醫(yī)院，320層的機(jī)器也已問(wèn)世。它們通過(guò)多排探測(cè)器技術(shù)來(lái)提供高掃描速度、高時(shí)間分辨率的體數(shù)據(jù)，同時(shí)采用心電門控技術(shù)來(lái)同步掃描，做到每個(gè)心動(dòng)周期掃描一周，能顯著減少心跳的偽影；通過(guò)指定不同的同步位置，還能獲得心臟不同搏動(dòng)時(shí)相的體數(shù)據(jù)。例如，以0.1s間隔，取R峰及前后共10個(gè)位置為掃描起始點(diǎn)，各心動(dòng)周期成像一次，就可獲得表征一個(gè)心動(dòng)周期不同搏動(dòng)時(shí)相的10個(gè)CT數(shù)據(jù)集，實(shí)現(xiàn)所謂的“四維CT成像”。對(duì)該4D-CT數(shù)據(jù)集做進(jìn)一步的處理和分析，就可能形成心、肺、冠狀動(dòng)脈等器官的動(dòng)態(tài)圖像，進(jìn)行心、肺功能的更深入研究。

四維CT數(shù)據(jù)處理的第一步是針對(duì)各個(gè)時(shí)相的數(shù)據(jù)集進(jìn)行心臟或肺的分割，工作量很大，必須實(shí)現(xiàn)自動(dòng)分割，而且算法要能保證較高的分割精度，這是后續(xù)研究順利進(jìn)行的根本保證。由于左心室在心臟功能研究中的重要地位，因此本研究以它為例，介紹從4D-CT心臟數(shù)據(jù)中分割出左心室腔和心肌的方法。在MRI心臟成像領(lǐng)域已有不少左心室分割的研究工作，但針對(duì)MSCT的左心室分割的報(bào)道較少，已有的研究基于主動(dòng)形狀模型或主動(dòng)表現(xiàn)模型，需要大量訓(xùn)練數(shù)據(jù)集，且不適合病變個(gè)體［1-2］。采用嵌入先驗(yàn)知識(shí)的改進(jìn)耦合水平集（improved couple level set，ICLS）模型，從MSCT心臟數(shù)據(jù)中精確提取左心室腔和心肌。下面在介紹水平集模型的基礎(chǔ)上，重點(diǎn)討論如何針對(duì)左心室分割的特點(diǎn)和先驗(yàn)知識(shí)，分別對(duì)水平集耦合因子、水平集起始輪廓、邊緣檢測(cè)因子進(jìn)行改進(jìn)的思路，并給出采用特定心相MSCT數(shù)據(jù)集進(jìn)行驗(yàn)證的結(jié)果。

1 原理和算法

1.1 水平集模型

水平集（level set）方法首先由Osher和Sethian提出，并用于曲線演化模型［3］。Caselles與Malladi等將該方法引入圖像處理領(lǐng)域，提出了幾何主動(dòng)輪廓模型，使初始曲線向能量泛函的極小值演化［4-5］，利用水平集方法，使模型能自適應(yīng)被檢物體的拓?fù)渥兓遗c參數(shù)無(wú)關(guān)。近年來(lái)，出現(xiàn)了許多基于水平集的圖像分割算法［6-9］。

在水平集模型中，用C表示演化曲線，C（t）=｛（x，y）|φ（t，x，y）=0｝是水平集函數(shù)φ（t，x，y）的零水平集，其中初始化輪廓為｛（x，y）|φ0（x，y）=0｝。假設(shè)I為待分割圖像，為使外部能量驅(qū)動(dòng)零水平集函數(shù)向目標(biāo)邊緣移動(dòng)，構(gòu)造圖像梯度的遞減函數(shù)為邊界檢測(cè)因子g，使得g在圖像邊緣的值近似為0。通常將g定義為

式中，Gσ是標(biāo)準(zhǔn)差為σ的高斯函數(shù)，*表示卷積運(yùn)算。

邊界檢測(cè)因子g在圖像分割中扮演極其重要的角色。從式（1）可看出，在圖像梯度較大區(qū)域，g值近似為0，在灰度一致區(qū)域?yàn)檎?。采用文獻(xiàn)［7］的WRLS模型為原始水平集模型，該模型包含了額外的能量項(xiàng)，用于控制水平集函數(shù)和符號(hào)距離函數(shù)的位置偏差，避免計(jì)算量大且復(fù)雜的重新初始化。能量函數(shù)由懲罰項(xiàng)作為內(nèi)部能量P（φ）和外部能量εg，λ，ν（φ）構(gòu)成，有

式中，λ（λ＞0）和ν為常數(shù)，δ為狄拉克函數(shù)，H為海維塞德函數(shù)。

外部能量項(xiàng)εg，λ，ν（φ）驅(qū)動(dòng)零水平集曲線向目標(biāo)邊界移動(dòng)，同時(shí)內(nèi)部能量項(xiàng)μP（φ）約束水平集函數(shù)趨近符號(hào)距離函數(shù)，從而保持水平集函數(shù)的平滑性。最小化能量函數(shù)，得到水平集演化方程，即

下面將針對(duì)MSCT圖像的特點(diǎn)以及左心室的先驗(yàn)知識(shí)改進(jìn)WRLS模型，使原始模型能適用于左心室的定位和精確分割。

1.2 水平集耦合因子

在結(jié)構(gòu)上，左心室可看成由類圓形的腔和類環(huán)形的心肌構(gòu)成，因此采用了不同的水平集函數(shù)φ1、φ2來(lái)表示左心室的腔和心肌。為避免兩水平集函數(shù)發(fā)生重疊，定義一個(gè)耦合水平集（coupled level set，CLS）因子，它的正負(fù)和權(quán)重控制著水平集曲線的方向和演化速度，即

式中，d表示從左心室重心到內(nèi)、外心膜的平均距離，w表示平均距離的容許范圍。

在實(shí)際應(yīng)用中，可根據(jù)左心室的結(jié)構(gòu)和大小對(duì)d、w取值。通過(guò)耦合水平集，水平集函數(shù)在內(nèi)、外心膜的對(duì)應(yīng)區(qū)域演化而不相互影響。

1.3 起始輪廓的初始化

1.3.1 左心室定位

由于分辨率很高，而且Z軸分辨率與層片內(nèi)的分辨率一致，即數(shù)據(jù)是各向同性的，所以MSCT數(shù)據(jù)具有結(jié)構(gòu)連續(xù)性的顯著特征，即相鄰層片間的結(jié)構(gòu)變化是連續(xù)的，結(jié)構(gòu)差異不大，灰度具有較大的相似性［10］。根據(jù)結(jié)構(gòu)連續(xù)特征，提出MSCT左心室的定位算法和粗輪廓提取算法，其中粗輪廓用于水平集起始函數(shù)的初始化。

設(shè)I是心臟某心相的左心室短軸MSCT的3D數(shù)據(jù)集，它包含N個(gè)層片In：Ω?R2→R+，n∈［1，2，…，N］。選取層片Is（s∈［1，2，…，N］）作為序列I的定位模板，其中Is需滿足條件：左心室腔面積在層片中較大（腔面積要比主動(dòng)脈等大血管面積大），左心室腔有較光滑的邊界，腔周圍有較小的組織噪聲。由于左心室腔的灰度值較大，邊緣較明顯，所以可根據(jù)閾值法和形態(tài)學(xué)操作獲得左心室腔的定位初始模板。

1）利用閾值法和形態(tài)學(xué)模型，去除胸腔的肺、骨頭、大血管等組織，初步獲取定位模板圖像中的左心室腔粗輪廓模板M，即

式中，T1和T2是去除左心室腔以外組織的低閾值和高閾值，*表示二值形態(tài)學(xué)操作，結(jié)構(gòu)元素SE1用于填充左心室腔中的小孔，Bwselect用于提取閾值分割結(jié)果，max用于提取出面積最大的左心室腔。

2）通過(guò)層片與模板相乘，提取出左心室腔的粗輪廓ΓtM為

式中，*表示矩陣對(duì)應(yīng)點(diǎn)的值相乘。

沿著左心室的短軸從腔由大變小的方向，使用初始定位模板來(lái)獲取鄰近層片左心室腔的粗輪廓，層片In的腔室掩板可由層片In-1得到，并利用圖1所示流程更新定位模板，從而獲取整個(gè)數(shù)據(jù)集中層片的左心室腔粗輪廓，完成左心室的定位。

圖1 左心室腔粗輪廓的獲取流程Fig.1 The flowchart to extract coarse contour of LV cavity

1.3.2 起始輪廓的初始化

水平集曲線演化之前必須確定初始化輪廓的位置、面積，不合適的起始位置、形狀和大小可能導(dǎo)致錯(cuò)誤的分割結(jié)果。在大多數(shù)水平集模型中，初始輪廓的形狀是圓形、矩形等規(guī)則形狀，與待分割目標(biāo)形狀無(wú)任何關(guān)系。本模型則按如下方法、利用左心室腔粗輪廓來(lái)定義起始形狀，有

式中，Γtn是左心室的粗輪廓，area（Γtn）表示左心室腔的區(qū)域，Centroid（ΓLVn-1）表示上一層片左心室的質(zhì)心，*表示二值形態(tài)學(xué)運(yùn)算，結(jié)構(gòu)因子SE3、SE4用于平滑區(qū)域輪廓。

利用結(jié)構(gòu)連續(xù)性和前一層片的分割結(jié)果，可以精確地定位腔室和心肌的起始輪廓，而且由于起始輪廓近似于待分割目標(biāo)邊緣，因此有效降低了水平集演化時(shí)間。

1.4 基于先驗(yàn)區(qū)域的邊緣檢測(cè)因子

邊緣檢測(cè)因子又稱為停止項(xiàng)，文獻(xiàn)［8］將其改進(jìn)為

式中，β是由高斯函數(shù)Gσ和引導(dǎo)因子λ卷積構(gòu)成的矩陣，*表示矩陣對(duì)應(yīng)點(diǎn)的值相乘。

引導(dǎo)因子定義為式中，Θ表示給定的感興趣區(qū)域。

引導(dǎo)因子可以把先驗(yàn)知識(shí)融合到曲線的演化當(dāng)中。由式（11）可知，當(dāng)曲線演化到感興趣區(qū)域時(shí)，引導(dǎo)因子退化為常數(shù)，g還原成式（1）的邊緣檢測(cè)因子，曲線保持正常的演化。當(dāng)引導(dǎo)因子為0時(shí)，曲線會(huì)越過(guò)不感興趣的區(qū)域。由于MSCT圖像層片間左心室的結(jié)構(gòu)和位置相對(duì)不變，假設(shè)Γn-1是上一層片圖像分割的左心室腔區(qū)域，則下一層片邊緣和層片n-1相近。本模型中設(shè)Γn-1為感興趣區(qū)域，則分割室腔的水平集引導(dǎo)函數(shù)可定義為

同樣，可以設(shè)置心肌的感興趣區(qū)域。通過(guò)重新設(shè)計(jì)引導(dǎo)函數(shù)，水平集演化方程在每個(gè)層片上都有對(duì)應(yīng)的檢測(cè)因子和停止函數(shù)，降低了曲線的演化泄露和曲線在目標(biāo)區(qū)域內(nèi)局部極值上的迭代。

1.5 改進(jìn)耦合水平集模型

把融合區(qū)域先驗(yàn)?zāi)Ｐ偷倪吘墮z測(cè)因子、水平集耦合因子加入水平集演化方程，便可得到ICLS模型的演化方程為

步驟1：通過(guò)式（7）和式（8），獲取層片In的左心室腔粗輪廓。

步驟3：通過(guò)式（10）和式（12），更新邊緣檢測(cè)因子gn。

步驟5：計(jì)算|φt+1n-φtn|＜θ，θ是自定義的小常數(shù)。如果對(duì)比水平集演化前后的面積差滿足該條件，便停止演化，否則繼續(xù)步驟4。

2 實(shí)驗(yàn)方法

在應(yīng)用ICLS模型分割前，必須對(duì)獲取的MSCT圖像數(shù)據(jù)進(jìn)行有效的預(yù)處理。首先對(duì)原始CT數(shù)據(jù)作非線性灰度變換，把一定范圍內(nèi)的灰度值映射到特定范圍內(nèi)，使得在該灰度范圍內(nèi)的圖像更加便于后續(xù)圖像處理［11］。此外，CT數(shù)據(jù)會(huì)受到數(shù)據(jù)采集系統(tǒng)的電子噪聲和重建噪聲等的影響。為提高分割質(zhì)量，必須去除其中的局部噪聲和斑點(diǎn)噪聲，同時(shí)要保留待分割目標(biāo)的邊緣。采用Paris等提出的快速雙邊濾波器［12］來(lái)去除CT圖像噪聲，與傳統(tǒng)的高斯濾波器相比，雙邊濾波器能較好地保持圖像邊緣，而高斯濾波后的左心室邊緣卻比較模糊（見圖2）。

圖2 濾波效果。（a）高斯濾波器；（b）雙邊濾波器Fig.2 Filtering results.（a）Gaussian filter；（b）bilateral filter

實(shí)驗(yàn)數(shù)據(jù)采自GE公司的64排螺旋CT，采用心電門控技術(shù)獲取的3D心臟MSCT數(shù)據(jù)，分辨率512像素×512像素×297像素；以及來(lái)自Philips公司的256排MSCT的8例3D心臟數(shù)據(jù)集，分辨率為512像素×512像素×331像素。首先選用64排3D心臟數(shù)據(jù)集，初步驗(yàn)證ICLS模型的分割效果。ICLS模型的參數(shù)設(shè)置為μ=0.04，λ=5，ν=3，ε=1.5，θ=10，曲線演化步長(zhǎng)Δt=5。為了進(jìn)一步驗(yàn)證算法對(duì)三維數(shù)據(jù)分割的適用性和有效性，保持ICLS模型的參數(shù)不變，對(duì)8例256排MSCT三維數(shù)據(jù)集進(jìn)行左心室心肌和血腔的分割與提取。

為了精確比較并評(píng)價(jià)自動(dòng)分割結(jié)果與手工分割結(jié)果的相似度，用ICLS模型自動(dòng)提取該數(shù)據(jù)集的短軸位［13］左心室腔和心肌，然后將自動(dòng)分割的結(jié)果與手工勾勒的結(jié)果進(jìn)行比較。采用如下方法［14］進(jìn)行相似度（Similarity，SM）對(duì)比：設(shè)CT層片經(jīng)計(jì)算機(jī)分割處理后，得到目標(biāo)面積S1、手工勾勒輪廓包圍的面積S2和共同包圍的面積S，相似度的計(jì)算公式為

SM的范圍為0～1。當(dāng)S1與S2完全重合時(shí)值為1，完全不重合時(shí)值為0；當(dāng)相似度數(shù)值越接近于1，表明算法的分割結(jié)果越接近于手工分割，算法的準(zhǔn)確性就越高。

3 結(jié)果

圖3是原水平集模型（WRLS）和改進(jìn)耦合水平集模型（ICLS）起始輪廓的對(duì)比，后者起始輪廓已近似于待分割輪廓，并能保證左心室的精確定位。圖4是兩種模型在相同參數(shù)下左心室的分割結(jié)果，可見在灰度近似區(qū)域ICLS模型減少了邊緣泄露。

圖3 起始輪廓。（a）WRLS模型；（b）ICLS模型Fig.3 The initial contour.（a）WRLS model；（b）ICLS model

圖4 分割結(jié)果。（a）WRLS模型；（b）ICLS模型Fig.4 Segmentation results.（a）WRLS model；（b）ICLS model

圖5是相似度分析后的結(jié)果，其中左心室腔的平均相似度0.950 3，標(biāo)準(zhǔn)偏差0.043 5；左心室心肌的平均相似度0.920 6，標(biāo)準(zhǔn)偏差為0.026 8。在腔和包圍腔的心肌之間，灰度對(duì)比度較為明顯，邊緣較為清晰，腔的定位和分割效果較好。心肌和鄰近組織的灰度相似度很高，心肌邊緣不是特別明顯。不過(guò)，除極個(gè)別層片外，其他層片的分割情況都較好，自動(dòng)分割與手工分割的結(jié)果基本吻合。圖6顯示的是對(duì)分割結(jié)果進(jìn)行3D表面重建后的結(jié)果，可以看出分割結(jié)果在左心房短軸位具有很好的同一性和完整性。

表1是ICLS模型與分割對(duì)心肌的分割結(jié)果，用式（14）計(jì)算獲得平均相似度?？傻肐CLS模型的血腔分割結(jié)果的平均相似度為0.959 2，平均標(biāo)準(zhǔn)偏差為0.019 0；心肌分割結(jié)果的平均相似度為0.906 3，平均標(biāo)準(zhǔn)偏差為0.025 6?？梢钥闯?，對(duì)于不同來(lái)源的三維左心室數(shù)據(jù)，ICLS模型同樣具有較高的分割精度，對(duì)心肌的平均分割準(zhǔn)確率可以達(dá)到95％以上，對(duì)心肌的分割準(zhǔn)確率可以達(dá)到90％以上，說(shuō)明ICLS算法具有普遍的適用性。因此，ICLS模型可以用于提取整個(gè)3D+t的MSCT心臟數(shù)據(jù)，這為后續(xù)多心相的心功能分析研究打下良好的基礎(chǔ)。

圖5 分割相似度。（a）左心室腔；（b）心肌Fig.5 Similaritiy of Segmentation results.（a）LV cavity；（b）myocardium

圖6 左心室三維重建結(jié)果。（a）不同角度下的心?。唬╞）不同角度下的血腔Fig.6 3D surface reconstruction results of LV.（a）myocardium from different views；（b）cavity from different views

4 結(jié)論

針對(duì)MSCT圖像左心室定位與分割的特定應(yīng)用，本研究提出了融合左心室先驗(yàn)形狀和結(jié)構(gòu)連續(xù)性的耦合水平集模型——ICLS模型。它首先根據(jù)結(jié)構(gòu)連續(xù)性和閾值法定位，提取出左心室腔的粗輪廓；然后利用左心室先驗(yàn)知識(shí)以及結(jié)構(gòu)連續(xù)性，指導(dǎo)后續(xù)水平集曲線的演化。ICLS模型對(duì)MSCT數(shù)據(jù)集的實(shí)驗(yàn)表明，它能實(shí)現(xiàn)左心室腔和心肌的自動(dòng)分割，與其他基于單幅梯度邊緣、無(wú)耦合因子、無(wú)先驗(yàn)知識(shí)的水平集方法比較，結(jié)果更準(zhǔn)確和完整，還提高了曲線演化速度并減少了邊緣泄露，這為進(jìn)一步分析左心室運(yùn)動(dòng)提供良好的數(shù)據(jù)保證。

表1 ICLS模型分割結(jié)果與手工分割結(jié)果的相似度與標(biāo)準(zhǔn)偏差Tab.1 The similarity and standard deviation between ICLS model and manual segmentation results

［1］Fritz D，Kroll J，Dillmann R，et al.Automatic 4D-segmentation of the left ventricle in cardiac-CT-data［C］.//Pluim JPW，Reinhardt JM.，eds.Medical Image 2007：Imaging Proceeding.Bellingham：SPIE-INT Soc Optical Engineering，2007，6512：N5123.

［2］Garreau M，Simon A，Boulmie D，et al.Assessment of left ventricular function in cardiacmsct imaging by a 4D hierarchical surface-volume matching process［J］.International Journal of Biomedical Imaging，2006，1-10.

［3］OsherS，Sethian JA.Fronts propagating with curvature dependent speed：algorithms based on the Hamilton-Jacobi formulation［J］.Journal of Computational Physics，1988，79：12-49.

［4］Caselles V，Catte F，Coll T，et al.A geometric method for active contours［J］.Numeriche Meathematik，1993，66：1-31.

［5］Malladi R，Sethian JA，Vemuri BC.Shape modeling with front propagation：a level set approach［J］.IEEE Trans Pattern Analysis and Machine Intelligence，1995，17（2）：158-175.

［6］Chan T，Vese L.Active contours without edges［J］.IEEE Transactions on Image Processing，2001，10（2）：266-277.

［7］Li CM，Xu CY，Gui C，et al.Level set evolution without reinitialization：a new variational formulation［C］.//Schmid C，Soatto S，Tomasi C.，eds.IEEE International Conference on Computer Vision and Pattern Recognition.Los Alamitos：IEEE Computer Soc，2005，1：430-436.

［8］Zhang H，Morrow P，McClean S，et al.Incorporating feature based priors into the geodesic active contour model and its application in biomedical imagery［C］//McDonald J，Markham C，Ghent J，eds.International Machine Vision and Image Processing Conference.Los Alamitos：IEEE Computer Soc，2007：67-74.

［9］Lynch M，Ghita O，Whelam PF.Left-ventricle myocardium segmentation using a coupled level set with a priori knowledge［J］.Computerized Medical Imaging and Graphics，2006，30：255-262.

［10］Liu Junwei.，F(xiàn)eng Huanqing，Zhou Yingyue，et al.A novel automatic extraction method of lung texture tree from HRCT image［J］.Acta Automatica Sinica，2009，35（4）：345-349.

［11］Heinemann EG.Simultaneous brightness induction as a function of inducing and test field luminances［J］.Journal of Experimental Psychology，1955，50（2）：89-96.

［12］Paris S，Durand F.A fast approximation of the bilateral filter using a signal processing approach［C］//Bischof H，Leonardis A.，eds.European conference on Computer Vision.Berlin：Springer-Verlag，2006，81（1）：24-52.

［13］American Heart Association Writing Group on Myocardial Segmentation and Registration for Cardiac Image.Standardized myocardial segmentation and nomenclature fortomographic imaging of the heart［J］.Circulation，2002，105：539-542.

［14］周穎玥，馮煥清，李傳富，等.一種從胸HRCT圖像序列分割肺的自動(dòng)化方法［J］.北京生物醫(yī)學(xué)工程，2008，27（1）：6-10.