全方向M型心動圖中非功能運動的分離和修正

楊秀芝 王衛(wèi)星 林 強(qiáng)

（福州大學(xué)物理與信息工程學(xué)院，福州 350108）

引言

心臟是人體的重要器官，它的運動與其功能緊密相關(guān)。因此在心臟疾病的診斷中，對心臟的運動分析是十分必要的。由于超聲診斷具有實時、直觀、無損傷、及時獲取等特點，為醫(yī)生對心臟的解剖結(jié)構(gòu)和動態(tài)過程的了解提供了全面的信息。這種方法在診斷心臟疾病時具有易于操作、靈活、成本低等特點，所以目前它仍是心臟疾病的重要檢查手段之一。

傳統(tǒng)的M型超聲心動圖是通過顯示采樣線上心臟結(jié)構(gòu)的運動曲線，測量室壁運動的各種參數(shù)（如室壁厚度、心室內(nèi)徑、室壁短軸縮短率、射血分?jǐn)?shù)、心血管運動時間等），對心肌功能進(jìn)行分析。這種檢測最大弱點是超聲的入射角不一定與運動方向一致，所獲得的速度誤差較大。解剖 M型（amatomic M-mode，AMM）超聲心動圖是對原始二維圖像的信息進(jìn)行后處理，從二維圖像上可任意調(diào)節(jié)采樣線的位置，獲取任意室壁階段或切面上感興趣的局部心肌運動曲線。但AMM超聲心動圖在同一個心動周期的二維圖像上只能設(shè)置三條取樣線，無法獲取同一心動周期內(nèi)所有室壁心肌運動的曲線。福州大學(xué)生物醫(yī)學(xué)工程研究所發(fā)明的全方向M型心動圖系統(tǒng)改進(jìn)了AMM超聲心動圖，可以在同一心動周期對不同部位的心肌運動進(jìn)行任意多條運動曲線的同步成像，能同時記錄下心臟內(nèi)任意結(jié)構(gòu)多個心動周期的 M型超聲圖像［1－2］。然而由于心臟受到呼吸、血液涌動以及各結(jié)構(gòu)的互相牽動造成的非功能性運動的影響，心臟室壁并不總是沿某條采樣線直線運動，初始時刻在這條采樣線上的心臟的某個結(jié)構(gòu)，在運動的過程中會脫離這條方向線，或者由于肌肉牽拉等在這條采樣線多（少）移動了一段距離，導(dǎo)致重建得到的灰度-時間波形圖含有其他結(jié)構(gòu)的運動信息。顯然，這樣所得到的灰度-時間波形圖不能準(zhǔn)確反映心臟結(jié)構(gòu)某個部位的功能運動情況。

修正心臟非功能性運動常用運動跟蹤、圖像配準(zhǔn)等方法［3－6］。圖像配準(zhǔn)常用于匹配一段時間內(nèi)來自不同或同一患者的醫(yī)學(xué)圖像，檢測的是圖像的相似度，對于含有心室舒縮運動的心臟序列圖像效果并不很理想?；谔卣鞯膱D像跟蹤是跟蹤圖像的整體或部分運動，達(dá)到圖像匹配的目的。文獻(xiàn)［7］采用遺傳算法跟蹤特征點的運動修正心臟非功能運動，但在心臟B超圖像中常常難以找到只含有非功能運動的特征點，因此它的應(yīng)用范圍有限。本研究通過獲取B超平面內(nèi)心臟的平移和心室的偏轉(zhuǎn)修正心臟的非功能性運動，改善獲取的心動波形圖的準(zhǔn)確性。因為B超圖像的局限性，僅能提供其平面內(nèi)的運動信息，因此為了減小心臟偏離B超平面的運動，在采集圖像時應(yīng)要患者盡量屏住呼吸，并且通過對超聲探頭位置和取向的準(zhǔn)確把握減小呼吸造成的心臟抖動和旋轉(zhuǎn)運動，進(jìn)而提高心臟運動參數(shù)測量的精準(zhǔn)性。針對心臟短軸B超圖像，采用基于心室輪廓線的運動跟蹤方法修正心臟非功能性運動。

1 改進(jìn)的Canny算法

在心臟短軸B超圖像中，心腔輪廓的徑向變化體現(xiàn)了心臟的功能性舒縮，而輪廓的整體平移和在二維B超平面中的旋轉(zhuǎn)體現(xiàn)了心臟的非功能性運動。目前邊界提取算法的改進(jìn)使得獲取心臟心室輪廓成為可能，一個正常人的心動周期大約為0.75～1 s。DICOM標(biāo)準(zhǔn)圖像為75幀/s，每一幀的時間間隔約為13 ms，這相對于心臟的運動周期要小得多，因此每幀圖像之間具有較強(qiáng)的空間相關(guān)性。為此，將室壁輪廓線的提取分為兩個階段：（1）首幀圖像室壁輪廓的提?。唬?）后續(xù)幀室壁輪廓的提取。經(jīng)過對多種邊界提取算法的實驗比較，首幀輪廓的提取采用改進(jìn)的 Canny算法［8－10］并輔以改進(jìn)的基于徑向灰度搜索算法［11］，后續(xù)幀采用改進(jìn)的 Canny算法并輔以相鄰幀匹配的方法獲得左心室短軸輪廓。下面介紹改進(jìn)的Canny算法。

在Canny算法中，由于二維高斯濾波器可以由兩個一維高斯濾波器疊加而成，所以討論一維高斯濾波器。其數(shù)學(xué)表達(dá)式為

參數(shù)σ2決定平滑圖像的程度，如果太小，則不能完全消除噪聲，而太大則會產(chǎn)生假邊緣。解決該矛盾的方法有兩種：一種是采用不同的濾波器和不同的平滑參數(shù)對圖像進(jìn)行預(yù)處理，再進(jìn)行圖像的邊緣檢測，最后通過分析這一系列結(jié)果來合成圖像的邊緣。另一種是采用Hodson提出的方法。該方法基于圖像的不同部分的噪聲和邊緣類型不同而采用不同的平滑方法。具體如下：

信號F（x）的高斯濾波可表示為。省略高階導(dǎo)數(shù)可近似表達(dá)為

定義預(yù)處理的誤差

所以有σ2（x）=2ε/|F（x）|（5）

采用后一種方法并進(jìn)行了改進(jìn)：假設(shè)可獲得的信號S（x）等于原始信號F（x）加噪聲N（x），則有

根據(jù)人的視覺系統(tǒng)對平緩部分的噪聲和邊緣比較敏感，而邊緣部分的二階導(dǎo)數(shù)比較大，所以 σ2必須盡量小；平緩部分的二階導(dǎo)數(shù)比較小，所以 σ2必須大，用式（7）代替式（5），有

式中，σs（x）是可獲得信號的平滑參數(shù)，σf（x）是原始信號的平滑參數(shù)，σn（x）是噪聲的平滑參數(shù)。再根據(jù)上述結(jié)論有

式中，k為比例系數(shù)。當(dāng)x在圖像邊緣時，F(xiàn)（x）變化很快，所以：（x）?，即：σ2（x）? k/（x）。當(dāng)x在圖像平緩部分：（x）?，σ2（x）? k，可以看出σ2是隨x的變化而變化的，因此在二維圖像中σ2是隨x，y的不同而不同的。根據(jù)這一特點，在圖像的不同位置采用不同的σ2，從而達(dá)到自適應(yīng)的效果。圖1是算法的部分結(jié)果，兩種方法結(jié)合提取心室輪廓可以獲得較好的效果。圖1中圖像來自一位男性心臟左心室短軸B超序列圖像，該段B超圖像描述了病人3個心動周期的心臟運動情況，周期長度2.5 min，共200幀。所有分析數(shù)據(jù)都來自這組B超圖像。

2 心室輪廓非功能性運動參數(shù)的獲取

正常的成年人心率約75次/min，平均每—心動周期約占0.8 s。一個心動周期里，心房收縮期占0.1 s，舒張期占0.7 s；而心室收縮期占0.3 s，舒張期占0.5 s。在心臟舒縮的各個期間，心室受到的血液沖擊、肌肉牽拉的方向和力量大小都不同，使每一幀的整體位移和在B超照射平面內(nèi)的偏轉(zhuǎn)也不相同。從輪廓圖中可以看出，由于心室的運動、乳頭肌的存在以及B超圖像的噪聲，所得到的輪廓形狀不規(guī)則，輪廓曲線含有高次諧波分量，稱輪廓曲線為類橢圓。隨著心臟的舒縮，這個類橢圓的面積大小不斷地變化。類橢圓的軸徑方向的變化表示了心臟的功能性運動—舒張和收縮，類橢圓的整體位移和偏轉(zhuǎn)則表達(dá)了心臟的非功能性運動。傅里葉描述子和質(zhì)心主軸法都可描述物體的整體運動而忽略細(xì)節(jié)部分（高頻分量），因此采用傅里葉描述子和經(jīng)典的質(zhì)心主軸法［12－14］對左心室的輪廓進(jìn)行了分析和計算，獲取其整體運動數(shù)據(jù)。

圖1 心臟左心室室壁輪廓圖Fig.1 Contours of left ventricular wall

設(shè)圖像的橫坐標(biāo)為實軸，縱坐標(biāo)為虛軸，心室輪廓由N點組成。從心室輪廓上的任一點開始繞輪廓曲線一周，可得到一個一維的復(fù)數(shù)序列

s（k）的離散傅里葉變換為

式中，u為頻率，S（u）為s（k）的頻譜系數(shù)。由傅里葉變換性質(zhì)可知，傅里葉變換的高頻分量對應(yīng)一些細(xì)節(jié)，低頻分量對應(yīng)輪廓的總體形狀，直流分量對應(yīng)輪廓所包圍區(qū)域的質(zhì)心位置

通過傅里葉逆變換，可以由傅里葉系數(shù) S（u）重建輪廓曲線s（k）：

因為傅里葉系數(shù)的高頻分量對應(yīng)輪廓細(xì)節(jié)，低頻分量對應(yīng)輪廓的總體形狀，可用對應(yīng)低頻分量的前M個傅里葉系數(shù)來近似地重建邊界，得到封閉邊界的一個近似輪廓（k）

根據(jù)輪廓重建實驗得知，當(dāng) M取40時，^s（k）可以保證既不丟失邊界的特征，同時忽略邊界的細(xì)節(jié)和微小波動，這樣也恰好減小了心室邊緣乳頭肌的影響以及輪廓獲取過程中的噪聲影響和邊界誤差。

為了消除心室邊緣噪聲的影響，心室輪廓的質(zhì)心用重建邊界來計算。設(shè)前后兩幀中的心室輪廓分別為si（k）和si+1（k），則第i幀的質(zhì)心位置為若輪廓在平面內(nèi)平移量為（Δx，Δy），則取坐標(biāo)原點在質(zhì)心位置，當(dāng)心臟的舒縮引起心室輪廓大小變化時，輪廓曲線可表示為式中，C為邊界放大（縮?。┑谋稊?shù)。它的傅里葉變換：

即傅里葉系數(shù)也放大C倍，因此由心臟舒縮引起的心室大小變化的影響可以通過傅里葉系數(shù)的歸一化消除。當(dāng)選取輪廓起點變化時，設(shè)k0為選取不同起點時所引起的序列原點平移量，則：sp（k）的傅里葉變換為：

取坐標(biāo)原點在心室輪廓的質(zhì)心上，輪廓曲線在空域旋轉(zhuǎn)一個角度θ后可表示為

其傅里葉變換為：

旋轉(zhuǎn)后的傅里葉系數(shù)等于原傅里葉系數(shù)乘以exp（jθ）。即不同頻率成分的傅里葉系數(shù)的相位都增加了θ角度，因此可通過傅里葉系數(shù)相位的變化規(guī)律來檢測邊界的轉(zhuǎn)動。

去掉直流分量，由式（21）可得前后兩幀圖像中對應(yīng)輪廓的前M個傅里葉系數(shù)的比值為

上式中，C為傅立葉系數(shù)的幅度比值，反映邊界的尺度縮放，而Δφ（u）為不同頻率下傅立葉系數(shù)的相位變化。由式（19）～（22）可得出前后兩幀圖像邊界的傅立葉系數(shù)的相位變化

由式（24）中的直線方程，已知不同頻率下傅立葉系數(shù)的相位變化，則可以通過Hough變換來確定直線的截距θ和斜率 q。由于u，Δφ（u）空間的一點對應(yīng)著θ，z參數(shù)空間上的一條直線［14］

在u，Δφ（u）空間共線的點所對應(yīng) θ，z參數(shù)空間的直線交于一點，該交點 θ，z值就是所要確定的式（25）中的截距和斜率。求出式（24）中的截距，則可以確定邊界曲線的二維平面內(nèi)的旋轉(zhuǎn)角度［14］。

利用心室輪廓的傅里葉變換的直流分量變化檢測輪廓在平面內(nèi)的平移，利用傅里葉系數(shù)的相位變化，再通過Hough變換來檢測輪廓的轉(zhuǎn)動。

質(zhì)心主軸法［15］是根據(jù)經(jīng)典力學(xué)物體質(zhì)量分布的原理，利用檢測出的輪廓曲線數(shù)據(jù)獲取心室的質(zhì)心位置和在二維平面中的主軸方向，通過比對質(zhì)心的坐標(biāo)和主軸角度，計算出輪廓的整體運動參數(shù)。

3 實驗結(jié)果及分析

比較數(shù)據(jù)對應(yīng)的是一段左心室短軸B超圖像的第 11 ～20、22、24、26、28、30、32、38、40、42、44、46、48～59幀，該段B超圖像的幀率為80幀/s，患者心率為74次/min，一個心動周期大約對應(yīng)65幀圖像，第11～38幀處于心室收縮期階段，40～59幀處于心室的舒張期階段。對于質(zhì)心位置比較，按照上述方法得到的心臟B超左心室質(zhì)心位置變化情況見圖2。

在圖2（b）中，在心室收縮到最小時（約38～44幀），輪廓中心的 y軸坐標(biāo)明顯減小，見圖中標(biāo)志，這是因為此時的心室輪廓下部收縮得很窄（見圖1），導(dǎo)致采用直接求輪廓質(zhì)心的y軸坐標(biāo)減?。ㄙ|(zhì)心位置上移）。采用傅里葉描述子重建邊界曲線后，減弱了乳頭肌等心內(nèi)結(jié)構(gòu)的影響，獲取的心室輪廓更接近于實際情況，得到的輪廓中心坐標(biāo)見圖2（a）。從圖2可看出，在心室的收縮和舒張過程中，心室輪廓的質(zhì)心做小幅的平移，即心室壁也做相同的平移，雖然每一幀的平移量不大，一般不超過2個像素，這里每一像素等效于0.44 mm，但多幀的累積位移仍會對采樣波形圖產(chǎn)生明顯的影響。因此對心臟平移運動的修正仍是必要的。

關(guān)于主軸偏轉(zhuǎn)角度的比較，計算了心室輪廓主軸方向的變化，如表1所示。表1中ya為按照傅里葉描述子方法計算的主軸角度，yb為按照質(zhì)心主軸法計算的主軸角度。從表1中可以看出，心室收縮和舒張期輪廓主軸的方向略有變化，收縮期主軸角度增大，舒張期主軸角度減小，在心室收縮到最小時，主軸近乎垂直。兩種方法計算得到的主軸變化規(guī)律和數(shù)值大小基本一致。

圖2 心臟左心室輪廓中心坐標(biāo)圖。（a）用傅里葉描述子計算；（b）用經(jīng)典力學(xué)原理計算Fig.2 The central coordinates of left ventricular contours.（a）By Fourier descriptors；（b）By principle of classical mechanics

表1 心室輪廓主軸角度列表Tab.1 Angle list of main axis in chamber contour

從上述結(jié)果可看出，兩種方法都可以分離心室的功能和非功能性運動，得到心室輪廓的非功能性運動參數(shù)，基于傅里葉描述子的計算結(jié)果更接近于心臟實際運動情況，但主軸質(zhì)心法計算簡單，需要更少的計算量。

4 全方向M型心動圖的修正

利用傅里葉描述子算法得到的每幀輪廓曲線的位移和在B超照射平面內(nèi)的偏轉(zhuǎn)，建立非功能性運動參考模型，指導(dǎo)采樣線的位置和方向做等效運動，使所要分析的室壁結(jié)構(gòu)點始終在該方向線上上下運動。

利用這種方法嘗試重做了全方向M型心動圖。圖3顯示了采用固定方向線（未去除非功能性運動的方向線）和跟蹤采樣線（去除了非功能性運動的方向線）時采集的全方向M型心動圖。圖中給出了4個采樣位置時的結(jié)果比較，每組圖的左側(cè)為采樣線的位置，右側(cè)為采樣結(jié)果。從圖3中可以看出經(jīng)過非功能運動的修正，采樣線基本跟蹤住了特定結(jié)構(gòu)，如圖3（a）所示，能夠獲得該結(jié)構(gòu)在整個心動周期中的運動情況；獲取的全方向M型心動圖周期更加明顯，波形具有了很強(qiáng)的重復(fù)性，如圖3中（a）、（b）、（d）所示；獲取的邊界輪廓清晰，基本不會出現(xiàn)特定點失落和跑偏情況，如圖3（b）所示；心動圖波形完整，不會丟失采樣結(jié)構(gòu)或出現(xiàn)其他組織結(jié)構(gòu)情況，如圖3中（b）和（c）所示；心動圖波峰比較尖銳，能夠反映出心臟受到?jīng)_力時的運動情況，而不被整體的大運動所平滑。

5 結(jié)論

圖3 消除非功能性運動前后的M型心動圖比較。（a）左上方采樣線位置；（b）上方采樣線位置；（c）左下方采樣線位置；（d）右下方采樣線位置Fig.3 Comparison of M-mode echocardiography before and after eliminating non-functional exercise.（a）The site of sampling on the right and upward side；（b）The site of sampling on the upward side；（c）The site of sampling on the left and downward side；（d）The site of sampling on the right and downward side

通過以上分析和實驗可以看出，用改進(jìn)的Canny算法并輔以基于徑向灰度搜索算法提取短軸左心室輪廓是可行的，結(jié)果令人滿意。采用基于心室輪廓線的運動跟蹤方法修正心臟非功能性運動是可行的，通過對非功能性運動的修正，采樣方向線基本能夠跟蹤心臟同一結(jié)構(gòu)某個部位在不同時刻的運動，能夠較準(zhǔn)確地反映心臟結(jié)構(gòu)特定部位的功能性運動。B超平面內(nèi)的非功能性運動可以通過心室輪廓的整體位移和主軸偏轉(zhuǎn)來表示。采用傅里葉描述子可以較準(zhǔn)確地獲得一系列圖像輪廓曲線的平移和主軸偏轉(zhuǎn)，但計算量較大；采用傳統(tǒng)的質(zhì)心主軸方法也可以快速獲得輪廓曲線的平移和偏轉(zhuǎn)，除心室輪廓變形劇烈的圖像外（如38～44幀），同樣可以改善心動圖的準(zhǔn)確性。通過對心臟非功能性運動的修正，可以使全方向M型心動波形圖更加準(zhǔn)確，進(jìn)而提高心臟運動速度、加速度等功能性動態(tài)特征參數(shù)的測量準(zhǔn)確性。采用本方法對其他患者的心臟短軸B超圖像進(jìn)行的分析實驗，可以得到類似的結(jié)果。所進(jìn)行修正的基礎(chǔ)建立在心室輪廓的準(zhǔn)確提取上，因此高效、快速、準(zhǔn)確的提取心室輪廓還將是筆者進(jìn)一步努力的方向。

［1］Lin Qiang，Wu Wenji，Huang Liqin，et al.An omnidirectional M-mode echocardiography system andits clinical application［J］.Computerized Medical Imaging and Graphics，2006，30：333 －338.

［2］林強(qiáng)，石江宏.超聲心動圖的一種動態(tài)信息——全方向M型心動圖［J］.儀器儀表學(xué)報，2005，26（4）：437－440.

［3］Guidry DL，F(xiàn)arag AA，Usingactivecontours and fourier descriptors for motion tracking with applications in MRI［C］//Mita D，eds.Proceedings of Image Processing（Volume 2）.Kobe：IEEE Signal Processing Society，1999：177－181.

［4］徐牧，王雪松，肖順平.基于目標(biāo)輪廓特征的SAR圖像目標(biāo)識別［J］.系統(tǒng)工程與電子技術(shù)，2006，28（12）：1812－1815.

［5］鄢余武，劉進(jìn)忙.非剛性醫(yī)學(xué)圖像的博弈配準(zhǔn)方法［J］.儀器儀表學(xué)報，2010，31（9）：2049 －2055.

［6］吳鋒，錢宗才，杭洽時.基于輪廓的力矩主軸法在醫(yī)學(xué)圖像配準(zhǔn)中的應(yīng)用［J］.第四軍醫(yī)大學(xué)學(xué)報，2001，22（6）：567－569.

［7］黃立勤，林強(qiáng).全方向M型心動圖系統(tǒng)精確檢測研究［J］.儀器儀表學(xué)報，2009，30（5）：1105 －1109.

［8］張斌，賀賽先.基于Canny算子的邊緣提取改善方法［J］.紅外技術(shù)，2006，28（3）：165－169.

［9］邵曉芳，孫即祥，王亮亮，等.改進(jìn)的Canny算法［J］.電光與控制，2006，13（6）：53 －55.

［10］Canny J.A computational approach to edge detection［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8（6）：679 －698.

［11］王同，余建國，王威琪.提取心臟截面圖像上內(nèi)腔的邊緣曲線［J］.復(fù)旦學(xué)報（自然科學(xué)版），2002，41（2）：174 －176，181.

［12］Wang Q，Clarke RJ.Estimation of general 2D motion using Fourier descriptors［J］.Vision，Image and Signal Processing，1994，141（2）：115－121.

［13］Oirrak AE，Daoudi M，Aboutajdine D.Estimation of general 2D affine motion using Fourier descriptors［J］.Pattern Recognition，2002，35（1）：223－228.

［14］鄭楚君，，楊志勇，何惠玲，等.傅里葉描述子和 Hough變換檢測封閉邊界運動［J］.計算機(jī)工程與應(yīng)用，2005，28：68－69，80.

［15］吳一全，付曉莉.采用角點信息和慣性主軸的車牌傾斜檢測與校正方法［J］.工程圖學(xué)學(xué)報，2009，6：127－131.