基于局部Radon變換和卡爾曼濾波的超聲圖像肌束方向自動(dòng)跟蹤方法

溫慧瑩 楊曉娟 郭燕榮 張 帥 張志國(guó) 汪天富 陳思平 陳 昕*

1(深圳大學(xué)醫(yī)學(xué)院生物醫(yī)學(xué)工程系 醫(yī)學(xué)超聲關(guān)鍵技術(shù)國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，廣東省生物醫(yī)學(xué)信息檢測(cè)與超聲成像重點(diǎn)實(shí)驗(yàn)室，深圳 518060)2(香港大學(xué)電氣與電子工程系，香港 999077)

基于局部Radon變換和卡爾曼濾波的超聲圖像肌束方向自動(dòng)跟蹤方法

溫慧瑩1楊曉娟1郭燕榮1張 帥2張志國(guó)2汪天富1陳思平1陳 昕1*

1(深圳大學(xué)醫(yī)學(xué)院生物醫(yī)學(xué)工程系 醫(yī)學(xué)超聲關(guān)鍵技術(shù)國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，廣東省生物醫(yī)學(xué)信息檢測(cè)與超聲成像重點(diǎn)實(shí)驗(yàn)室，深圳 518060)2(香港大學(xué)電氣與電子工程系，香港 999077)

超聲成像是一種無(wú)創(chuàng)、實(shí)時(shí)和便捷的成像方法，已被廣泛用于人體肌肉運(yùn)動(dòng)分析。臨床上目前主要采用手工處理的方法來(lái)提取超聲圖像中的形態(tài)結(jié)構(gòu)，不僅依賴于操作者的主觀經(jīng)驗(yàn)，且耗時(shí)嚴(yán)重、重復(fù)性低。提出一種自動(dòng)計(jì)算與跟蹤超聲圖像序列中肌束方向的新方法，利用基于約束互信息的自由變換算法自動(dòng)提取感興趣區(qū)域，然后基于局部Radon變換對(duì)超聲圖像進(jìn)行投影，并利用相鄰幀間肌束的連貫性，采用卡爾曼濾波的圖像處理方法。對(duì)6名試驗(yàn)對(duì)象進(jìn)行步行實(shí)驗(yàn)，共采集小腿內(nèi)側(cè)腓腸肌的動(dòng)態(tài)超聲影像1 080幀，并用所提出的方法計(jì)算肌束方向。實(shí)驗(yàn)結(jié)果表明，所提出的自動(dòng)檢測(cè)方法有較好的魯棒性。與人工檢測(cè)結(jié)果相比，兩者的平均相關(guān)系數(shù)為0.92±0.02，平均誤差為0.30°±0.62°。該方法可以一定程度地替代人工測(cè)量，在處理大量超聲數(shù)據(jù)時(shí)具有較大的優(yōu)勢(shì)。

超聲成像；肌束方向；Radon變換；卡爾曼濾波

引言

肌肉是構(gòu)成人體的重要組織，其主要功能為產(chǎn)生收縮和控制身體的運(yùn)動(dòng)。肌電圖學(xué)(electromyography, EMG)是最早用于分析人體肌肉運(yùn)動(dòng)的方法，其原理是記錄神經(jīng)肌肉收縮和靜止?fàn)顟B(tài)時(shí)的生物電信號(hào)。其中，表面肌電圖因其實(shí)時(shí)、無(wú)創(chuàng)等優(yōu)點(diǎn)，被廣泛用于肌肉運(yùn)動(dòng)研究中[1-2]。但是，表面肌電圖估計(jì)也存在一些不足，如對(duì)體內(nèi)深層肌肉難以區(qū)別、對(duì)干擾比較敏感等，這制約了肌電信號(hào)在肌肉運(yùn)動(dòng)分析中的應(yīng)用。另外，從上世紀(jì)90年代起，有學(xué)者開(kāi)始利用超聲評(píng)估肌肉的功能狀態(tài)，并將分析結(jié)果應(yīng)用在生物力學(xué)的研究領(lǐng)域[2-8]。超聲成像是一種無(wú)創(chuàng)、實(shí)時(shí)和便捷的成像方法。二維超聲圖像可提供肌肉收縮時(shí)肌束的實(shí)時(shí)動(dòng)態(tài)圖像，已被廣泛用于研究肌肉運(yùn)動(dòng)中形態(tài)結(jié)構(gòu)的改變，為臨床診斷和康復(fù)評(píng)估提供了途徑。

利用超聲可測(cè)量肌肉厚度(muscle thickness)[3, 6, 8-9]、羽狀角(pennation angle)(或稱肌束方向)[3, 5-8]、肌束長(zhǎng)度(fascicle length)[3, 5-8]和肌肉橫截面積(cross-section area)[4]等結(jié)構(gòu)性參數(shù)來(lái)描述肌肉的狀態(tài)變化，這些參數(shù)和肌肉組織的力學(xué)特性直接相關(guān)。其中，肌束方向描述肌纖維的排列，是決定肌肉力學(xué)性能的主要特征之一。在肌肉收縮過(guò)程中，羽狀角越大，更多的平行排列且可收縮的肌纖維束將能量傳遞到肌腱上，使肌肉產(chǎn)生更大的力量。肌肉處于最大收縮狀態(tài)時(shí)，羽狀角的變化能達(dá)到靜息時(shí)的120%～170%[10-11]。最近的研究還表明，肌肉肌束特性的改變與神經(jīng)肌肉疾病相關(guān)[12-14]。目前，臨床上主要是通過(guò)手工處理的方法來(lái)提取肌束方向信息。這不僅依賴于操作者的主觀經(jīng)驗(yàn)，在處理實(shí)時(shí)超聲視頻數(shù)據(jù)時(shí)耗時(shí)嚴(yán)重，而且需通過(guò)人工多次測(cè)量來(lái)提高重復(fù)性。

多種自動(dòng)提取肌束方向的算法已被提出，根據(jù)方法大致可分為兩類：Hough變換和Radon變換。Zhou等[15]首先提出了用迭代Hough變換(revoting hough transform, RVHT)，計(jì)算圖像中的肌束方向。該方法通過(guò)將超聲的邊緣圖像映射到霍夫空間，并將局部極值對(duì)應(yīng)于該圖像的直線位置，在檢測(cè)線設(shè)定線寬范圍內(nèi)的所有特征點(diǎn)將被移除，重復(fù)上述過(guò)程，多條直線被依次提取，直至低于自定義閾值。然而，該方法極大地依賴邊緣檢測(cè)子的性能，容易受斑點(diǎn)噪聲的干擾。此外，RVHT的檢測(cè)結(jié)果高度依賴于所涉及迭代過(guò)程中的線寬參數(shù)，稍微不同的參數(shù)大小可能會(huì)導(dǎo)致完全不同的結(jié)果。臨床中的肌肉超聲圖像，通常很難找到一個(gè)合適的線寬代表所有的肌肉纖維。Zhao[16]采用局部Radon變換(localized radon transform, LRT)，自動(dòng)計(jì)算和跟蹤肌束方向，不需要邊緣檢測(cè)子，受斑點(diǎn)噪聲干擾小。然而，相同的迭代策略用于依次提取直線，仍存在相同的問(wèn)題，如果初始檢測(cè)不正確，后續(xù)難以恢復(fù)正確的肌束線結(jié)構(gòu)[17]。最近文獻(xiàn)[18]提出，利用局部梯度圖來(lái)估計(jì)主要的肌束方向，而不是單獨(dú)檢測(cè)幾根肌束方向，因此迭代策略是不必要的，并且更具穩(wěn)健性，但是對(duì)局部方向上的每個(gè)像素的計(jì)算復(fù)雜度比[15-17]高很多。此外，上述方法存在兩個(gè)缺點(diǎn)：一是需要在超聲圖像序列的第一幀圖像上設(shè)置感興趣區(qū)域(reading of interest, ROI)，對(duì)所有圖像的肌束方向跟蹤都是在該ROI范圍中進(jìn)行。然而，在真實(shí)的肌肉收縮過(guò)程中，深淺層筋膜區(qū)域是動(dòng)態(tài)變化的，可能超過(guò)了固定ROI的范圍，這將大幅度降低自動(dòng)跟蹤肌束方向的精度；二是在追蹤過(guò)程中沒(méi)考慮運(yùn)動(dòng)模型，因此容易受干擾，得到的結(jié)果曲線常出現(xiàn)毛刺。以往的研究都是針對(duì)靜態(tài)或者準(zhǔn)靜態(tài)運(yùn)動(dòng)，干擾相對(duì)比較小。對(duì)步行運(yùn)動(dòng)，超聲探頭相對(duì)滑動(dòng)等干擾因素比較嚴(yán)重，必須設(shè)法消除。

根據(jù)對(duì)現(xiàn)有方法的分析，本研究提出一種新的自動(dòng)化計(jì)算與跟蹤肌肉肌束方向的方法。與傳統(tǒng)檢測(cè)肌束直線的方法相比，該方法自動(dòng)檢測(cè)與跟蹤ROI，利用Gabor濾波來(lái)自適應(yīng)地增強(qiáng)超聲圖像，并應(yīng)用局部Radon變換和卡爾曼濾波(Kalman filter, KF)來(lái)自動(dòng)計(jì)算和跟蹤肌束方向，提高了算法的魯棒性和準(zhǔn)確性。

1 方法

1.1 方法流程

研究方法的流程如圖1所示。對(duì)連續(xù)的超聲影像，首先基于約束互信息的自由變換算法(constrained mutual information-based free-form deforamtion, C-MI-FFD)，實(shí)現(xiàn)ROI區(qū)域的自動(dòng)化選取，采用Gabor濾波抑制圖像斑點(diǎn)噪聲，并縱向增強(qiáng)圖像肌束區(qū)域。此外，使用遲滯閾值法(hysteresis thresholding, HT)和圖像細(xì)化(image thinning, IT)進(jìn)一步增強(qiáng)圖像，從而解決迭代策略中的線寬問(wèn)題。利用文獻(xiàn)[16]所述的方法來(lái)捕獲肌束方向，并為KF提供校正測(cè)量。由于兩個(gè)超聲圖像相鄰幀之間的肌束具有連貫性，可取前一幀KF預(yù)測(cè)的肌束方向變化作為先驗(yàn)知識(shí)來(lái)追蹤當(dāng)前圖像內(nèi)的肌束。

圖1 方法流程Fig.1 The diagram of the proposed method

1.2 基于約束互信息的自由變換算法

文獻(xiàn)[4]首次提出基于約束互信息的自由變換算法(C-MI-FFD)，即在連續(xù)的超聲圖像序列中自動(dòng)提取肌肉的橫截面積。本研究擴(kuò)展了該算法，自動(dòng)選取連續(xù)超聲圖像序列深淺層筋膜間的感興趣區(qū)域。通過(guò)最小化以互信息為基礎(chǔ)的目標(biāo)函數(shù)，確定描述連續(xù)兩幀圖像的轉(zhuǎn)換函數(shù)。為了精確地匹配感興趣區(qū)域，C-MI-FFD算法分為兩個(gè)步驟。

步驟1:全局變換。通過(guò)全局匹配兩幀圖像確定全局變換參數(shù)，對(duì)它們相對(duì)比例、平移和旋轉(zhuǎn)建模。

步驟2：局部變換。局部變換參數(shù)由二維樣條函數(shù)定義，并被細(xì)化為規(guī)則網(wǎng)絡(luò)下進(jìn)行插值樣條函數(shù)的位移值。圖2(a)、(b)分別為典型的原始和自動(dòng)獲取ROI后的超聲圖像。

1.3 Gabor濾波器

由于肌肉超聲圖像具有連貫性，采用Gabor濾波可以抑制斑點(diǎn)噪聲，并在縱向上增強(qiáng)肌束圖像[19]。Gabor濾波器(Gabor filter, GF)為帶通濾波器，目前廣泛應(yīng)用于圖像特征提取、邊緣檢測(cè)、圖像紋理分析等。二維Gabor函數(shù)的實(shí)數(shù)部分定義如下：

(1)

式中,λ為波長(zhǎng)，θ為Gabor函數(shù)并行條紋的方向，φ為相位偏移，σ為高斯因子標(biāo)準(zhǔn)差，γ為空間縱橫比，(x,y)為圖像坐標(biāo)。

圖2(b)、(c)分別表示Gabor濾波前后的超聲圖像?？梢钥闯?，使用Gabor濾波后的圖像斑點(diǎn)噪聲被抑制，肌束區(qū)域圖像被增強(qiáng)。值得注意的是，GF中的θ是自適應(yīng)更新的、由前一幀KF預(yù)測(cè)的方向作為處理當(dāng)前幀的并行條紋的方向。

為解決LRT迭代策略中的線寬問(wèn)題和提升其檢測(cè)精度，采用HT、IT兩個(gè)技術(shù)進(jìn)一步增強(qiáng)GF濾波后的圖像信息。

步驟1：首先在HT中設(shè)有兩個(gè)閾值，像素響應(yīng)高于閾值上界的輸出為1，而低于閾值下界的輸出為0，像素響應(yīng)位于閾值上界與下界之間的，且如果它們通過(guò)一條邊緣鏈連接到其他像素響應(yīng)均大于下界閾值的輸出為1。

步驟2：圖像細(xì)化。通過(guò)抑制非極大值部分的輸出，可減少輸出線結(jié)構(gòu)的單像素邊緣寬度，從而得到非常小的線寬，且LRT中重新表決的線寬問(wèn)題得到有效解決。圖2(d)、(e)分別表示HT和IT后的效果。

圖2 方法過(guò)程追蹤。(a)原始超聲圖像；(b)提取ROI(綠色線表示)；(c)Gabor濾波效果；(d)HT效果；(e)IT效果LRT得到的主要肌束方向用藍(lán)色線表示Fig.2 Illustration of the proposed muscle fiber orientation tracking approach. (a) A raw ultrasound video frame; (b) Extracted ROI (green line surrounded area in (a)); (c) GF result; (d) HT result; (e) IT result. The dominant orientation obtained by using LRT is labeled by blue line

1.4 局部Radon變換和線提取

Radon變換已廣泛用于直線特征的提取，標(biāo)準(zhǔn)Radon變換在2D變換空間內(nèi)的定義為

(4)

式中,f(x,y)為在(x,y)位置的圖像灰度，δ為狄拉克δ函數(shù)，ρ為直線到圖像中心的距離，θ為x軸和線的夾角，D為圖像網(wǎng)格。

標(biāo)準(zhǔn)Radon變換表示整個(gè)圖像平面內(nèi)的灰度積分。在具體應(yīng)用中，常采用矩形來(lái)定義ROI，無(wú)法真實(shí)反映淺層和深層筋膜之間的區(qū)域的真實(shí)輪廓。如果只對(duì)預(yù)定義范圍內(nèi)部的位置和方向積分，標(biāo)準(zhǔn)Radon變換則可修改為L(zhǎng)RT。運(yùn)用先驗(yàn)知識(shí)如追蹤ROI和肌束方向角度，LRT可以顯著增強(qiáng)捕獲肌束的準(zhǔn)確度，降低邊緣效應(yīng)。在本研究中，采用C-MI-FFD算法可排除有強(qiáng)回聲強(qiáng)度的淺層與深層筋膜，且同時(shí)最大程度地跟蹤其余區(qū)域。當(dāng)檢測(cè)到指定的直線數(shù)量，或是低于峰值振幅與平均振幅的比值，表明這幀圖像中所有直線已檢測(cè)完畢。

1.5 卡爾曼濾波

卡爾曼濾波是一種高效率的遞歸濾波器，它能夠從一系列的不完全及包含噪聲的測(cè)量中估計(jì)動(dòng)態(tài)系統(tǒng)的狀態(tài)。目前廣泛用于工程領(lǐng)域，如雷達(dá)、計(jì)算機(jī)視覺(jué)，近來(lái)更被應(yīng)用于計(jì)算機(jī)圖像處理，如人臉識(shí)別、圖像分割、圖像邊緣檢測(cè)等[20]。

正如上文描述，方向參數(shù)是Gabor濾波器和后續(xù)LRT直線特征提取的關(guān)鍵參數(shù)。筆者提出用KF自適應(yīng)地更新方向的變化，即用LRT獲得當(dāng)前幀的肌束方向作為KF的測(cè)量值，KF輸出的最優(yōu)估計(jì)值用于Gabor濾波器處理下一幀圖像中的方向參數(shù)的反饋輸入，使達(dá)到更好的圖像處理效果。KF的離散時(shí)間線性狀態(tài)空間模型如下：

(5)

zk=Ckxk+Dkyk+vk

(6)

相對(duì)傳統(tǒng)的LRT，使用KF后，跟蹤肌束方向的精度得到提高，運(yùn)行時(shí)間顯著降低。

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

選取試驗(yàn)對(duì)象6名，均為男性。年齡(24.6±1.5)歲，體重(69.4±6.3)kg，身高(170±5.6)cm。所有試驗(yàn)對(duì)象身體健康，均沒(méi)有肌肉、神經(jīng)類疾病病史。受試者均簽署知情同意書(shū)。

使用邁瑞DC-6B超實(shí)時(shí)采集肌肉的動(dòng)態(tài)超聲圖像，配備線陣超聲探頭，中心頻率為10MHz。用一個(gè)定制的可調(diào)節(jié)的支架固定超聲探頭，使探頭的長(zhǎng)軸與試驗(yàn)對(duì)象小腿內(nèi)側(cè)腓腸肌的表面走向平行，并通過(guò)大量超聲耦合劑使探頭緊貼在目標(biāo)肌肉表面。實(shí)驗(yàn)開(kāi)始前，試驗(yàn)對(duì)象被要求在跑步機(jī)上做適應(yīng)性步行，一般為2～3min。實(shí)驗(yàn)中，當(dāng)試驗(yàn)對(duì)象準(zhǔn)備好后，開(kāi)啟跑步機(jī)并且調(diào)節(jié)到合適的速度1.0km/h，試驗(yàn)對(duì)象在跑步機(jī)上做正常的步行運(yùn)動(dòng)。超聲掃描儀的視頻輸出通過(guò)視頻采集卡(大恒公司，CG400)來(lái)采集，視頻幀率為25Hz。實(shí)驗(yàn)中持續(xù)采集30s的視頻數(shù)據(jù)，對(duì)起始和結(jié)尾時(shí)間段的視頻進(jìn)行剪切，去除運(yùn)動(dòng)以及不穩(wěn)定因素的影響，最終對(duì)每名試驗(yàn)對(duì)象提取180幀連續(xù)的超聲圖像序列。

3 結(jié)果

利用所提出的方法，對(duì)6名試驗(yàn)對(duì)象(共6×180=1 080幀)進(jìn)行檢測(cè)。以人工測(cè)量的結(jié)果作為參考標(biāo)準(zhǔn)，在對(duì)所有方法的結(jié)果進(jìn)行比較和分析后，進(jìn)一步證明本方法的正確性。所提方法在Matlab2013b版本下實(shí)現(xiàn)。

在本實(shí)驗(yàn)中，手動(dòng)過(guò)程由兩位臨床專家完成。人工測(cè)量時(shí)，每一位臨床專家在每幀圖像的深淺層筋膜的感興趣區(qū)域內(nèi)，選擇標(biāo)記三段顯著的肌束，取平均后得到最終的人工測(cè)量結(jié)果，如圖3所示。

圖3 手工測(cè)量示意圖Fig.3 The diagram of the manual method

在步行運(yùn)動(dòng)中，肌肉收縮引起肌肉快速運(yùn)動(dòng)，并產(chǎn)生了較大的運(yùn)動(dòng)模糊。此時(shí)，人工和自動(dòng)檢測(cè)都很難準(zhǔn)確地檢測(cè)肌束方向，因此數(shù)據(jù)分析中將先排除因運(yùn)動(dòng)模糊而無(wú)法檢測(cè)的超聲圖。跟蹤的典型結(jié)果如圖4所示，其中取兩名臨床專家的測(cè)量結(jié)果的平均值作為人工測(cè)量的最終結(jié)果，用黑色虛線表示，本方法用黑色實(shí)線表示?？傮w來(lái)看，本方法能捕捉正常人步行中肌束方向的變化趨勢(shì)，與人工測(cè)量結(jié)果相差較小。進(jìn)一步分析同一組人工檢測(cè)與自動(dòng)方法檢測(cè)結(jié)果的相關(guān)性，Pearson相關(guān)系數(shù)r=0.93(見(jiàn)圖5)，即本方法與人工測(cè)量結(jié)果有高度相關(guān)性。

圖4 肌束方向跟蹤結(jié)果Fig.4 Tracking result of the muscle fascicle orientation

圖5 人工與自動(dòng)檢測(cè)結(jié)果的相關(guān)性分析(實(shí)線是線性趨勢(shì)線)Fig.5 The correlation analysis of the manual and automated test result, the sold line suggests the linear trend line

表1 自動(dòng)測(cè)量方法和人工測(cè)量方法比較Tab.1 The comparison between automated and manual measurement

表1總結(jié)了6名試驗(yàn)對(duì)象的肌束方向測(cè)量結(jié)果，包括人工與筆者提出的自動(dòng)測(cè)量結(jié)果的相關(guān)系數(shù)與誤差。從整體上看，相關(guān)系數(shù)的范圍是0.88～0.94，在統(tǒng)計(jì)學(xué)意義上的平均相關(guān)系數(shù)為0.92±0.02，即人工與自動(dòng)檢測(cè)結(jié)果有很好的相關(guān)性,證明本方法的測(cè)量精度較高。此外，自動(dòng)和人工檢測(cè)方法的平均誤差為0.30°±0.62°，可以推斷該方法的準(zhǔn)確性能夠滿足臨床需求?？梢?jiàn)，本方法在連續(xù)檢測(cè)正常步行中肌肉肌束的方向具有很強(qiáng)的魯棒性。

4 討論

實(shí)驗(yàn)表明，所提出的自動(dòng)測(cè)量方法與人工測(cè)量的平均誤差達(dá)到0.30°±0.62°，Pearson相關(guān)系數(shù)為0.92±0.02，性能優(yōu)于或者接近其他的一些自動(dòng)測(cè)量方法。例如，Zhou等提出的迭代Hough變換算法，平均誤差是0.18°±2.41°，Pearson相關(guān)系數(shù)為0.98[15]。Zhao等提出的局部Radon變換算法，誤差約為1°[16]。而且，上述兩個(gè)研究均是在靜態(tài)或準(zhǔn)靜態(tài)的情況下測(cè)得的，深淺層筋膜區(qū)域被設(shè)為固定。本研究是針對(duì)正常步行運(yùn)動(dòng)進(jìn)行分析，運(yùn)動(dòng)模型和噪聲等問(wèn)題更加復(fù)雜，自動(dòng)化計(jì)算與跟蹤超聲圖像中肌束方向的難度增大，因此本方法的準(zhǔn)確性和魯棒性更高。此外，本研究對(duì)提出方法的驗(yàn)證更充分。在實(shí)驗(yàn)設(shè)計(jì)中，搜集了6名試驗(yàn)對(duì)象，對(duì)每名試驗(yàn)對(duì)象處理至少180幀連續(xù)超聲圖像，總計(jì)處理幀數(shù)為1 080幀，充分驗(yàn)證了測(cè)量肌束方向新方法的穩(wěn)定性和準(zhǔn)確性。其他的研究一般只采用了1～2名試驗(yàn)對(duì)象, 圖像處理幀數(shù)有限[15-16]。綜上所述，本方法在自動(dòng)化計(jì)算與跟蹤肌束方向上有較高精度，且可實(shí)現(xiàn)正常步行運(yùn)動(dòng)中對(duì)感興趣區(qū)域的自動(dòng)化選取，實(shí)現(xiàn)了肌束方向的全自動(dòng)測(cè)量。

所提方法的計(jì)算復(fù)雜度主要集中在局部Radon變換的方向估計(jì)。若采用M表示ROI尺寸的規(guī)模， 表示Radon方向的總數(shù)規(guī)模，則LRT方法的計(jì)算復(fù)雜度為O(θ·M2)，而應(yīng)用中ROI的尺寸常小于圖像整體尺寸，因此該方法的實(shí)際運(yùn)算效率較高。下一步工作將采用執(zhí)行效率更高的語(yǔ)言(如C語(yǔ)言)來(lái)實(shí)現(xiàn)本方法，有望實(shí)現(xiàn)肌束角度的實(shí)時(shí)在線測(cè)量。

在實(shí)驗(yàn)中，誤差產(chǎn)生的原因有多個(gè)方面。首先，較靜態(tài)或準(zhǔn)靜態(tài)下超聲在肌肉運(yùn)動(dòng)的評(píng)估，在正常步行運(yùn)動(dòng)分析中，運(yùn)動(dòng)模型和噪聲等問(wèn)題更加復(fù)雜，難以準(zhǔn)確、動(dòng)態(tài)地計(jì)算與跟蹤肌束方向。本研究引入了一階運(yùn)動(dòng)模型，今后可能考慮在步行中采用更高階的運(yùn)動(dòng)模型，進(jìn)一步提高測(cè)量精度。其次，步行運(yùn)動(dòng)中肌肉收縮引起肌束快速運(yùn)動(dòng)，產(chǎn)生較大的運(yùn)動(dòng)模糊，對(duì)自動(dòng)和人工檢測(cè)方法的準(zhǔn)確性檢測(cè)造成一定的難度。在以后的研究中，將采用更高的視頻采樣幀率和更準(zhǔn)確的讀取圖像方法，獲取更快步行速度和更清晰的超聲圖像序列。最后，實(shí)驗(yàn)中如何在步行運(yùn)動(dòng)中固定超聲探頭，防止探頭不規(guī)則滑動(dòng)影響超聲數(shù)據(jù)的采集，也是急需解決的問(wèn)題。要從根本上解決該問(wèn)題，需要開(kāi)發(fā)出更加輕便靈活的超聲探頭。

5 結(jié)論

綜上所述，本研究實(shí)現(xiàn)了對(duì)感興趣區(qū)域的自動(dòng)提取，基于LRT和卡爾曼濾波，自動(dòng)化計(jì)算與跟蹤超聲圖像序列中的肌束方向。實(shí)驗(yàn)證明，在正常步行運(yùn)動(dòng)中，本研究提出的方法對(duì)連續(xù)檢測(cè)肌束方向有較高的魯棒性。選用卡爾曼濾波器作為輔助跟蹤的運(yùn)動(dòng)模型，預(yù)測(cè)存在觀測(cè)噪聲情況下的運(yùn)動(dòng)估計(jì)，然后再對(duì)估計(jì)的方向參數(shù)進(jìn)行修正，并利用Gabor濾波自適應(yīng)地增強(qiáng)圖像肌束區(qū)域，提高了測(cè)量結(jié)果的準(zhǔn)確性。本方法與臨床人工檢測(cè)結(jié)果具有高度的一致性，可以一定程度地替代人工測(cè)量，在處理大量超聲數(shù)據(jù)時(shí)具有較大的優(yōu)勢(shì)。隨著超聲在肌肉檢測(cè)上的應(yīng)用，利用重復(fù)性好的自動(dòng)方法定量分析肌肉肌束方向參數(shù)，在臨床診斷和康復(fù)評(píng)估上具有巨大潛力。

[1] Disselhorst-Klug C, Schmitz-Rode T, Rau G, et al. Surface electromyography and muscle force: Limits in semg-force relationship and new approaches for applications [J]. Clin Biomech (Bristol, Avon), 2009, 24(3): 225-35.

[2] 李喬亮，易萬(wàn)貫，陳昕，等. 融合實(shí)時(shí)超聲影像的多模態(tài)運(yùn)動(dòng)特性研究 [J]. 中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào), 2012, 31(4): 518-525.

[3] Blazevich AJ, Gill ND, Zhou S, et al. Intra- and intermuscular variation in human quadriceps femoris architecture assessed in vivo [J]. J Anat, 2006, 209(3): 289-310.

[4] Chen Xin, Zheng Yongping, Guo JingYi, et al. Sonomyographic responses during voluntary isometric ramp contraction of the human rectus femoris muscle [J]. Eur J Appl Physiol, 2012, 112(7): 2603-2614.

[5] Fukunaga T, Kubo K, Kawakami Y, et al. In vivo behaviour of human muscle tendon during walking [J]. Proc Biol Sci, 2001, 268(1464): 229-233.

[6] Hodges PW, Pengel LH, Herbert RD, et al. Measurement of muscle contraction with ultrasound imaging [J]. Muscle Nerve, 2003, 27(6): 682-692.

[7] Kawakami Y, Ichinose Y, Fukunaga T, et al. Architectural and functional features of human triceps surae muscles during contraction [J]. J Appl Physiol (1985), 1998, 85(2): 398-404.

[8] Maganaris CN, Baltzopoulos V, Sargeant AJ, et al. In vivo measurements of the triceps surae complex architecture in man: Implications for muscle function [J]. J Physiol, 1998, 512 (Pt 2): 603-614.

[9] 李喬亮，任盼盼，張會(huì)生，等. 基于光流的超聲圖像肌肉厚度自動(dòng)測(cè)量方法 [J]. 中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào), 2013, 32(2): 149-153.

[10] Herbert RD, Gandevia SC. Changes in pennation with joint angle and muscle torque-in-vivo measurements in human brachialis muscle [J]. Journal of Physiology-London, 1995, 484(2): 523-532.

[11] Narici MV, Binzoni T, Hiltbrand E, et al. In vivo human gastrocnemius architecture with changing joint angle at rest and during graded isometric contraction [J]. Journal of Physiology-London, 1996, 496(1): 287-297.

[12] Gao F, Zhang LQ. Altered contractile properties of the gastrocnemius muscle poststroke [J]. Journal of Applied Physiology, 2008, 105(6): 1802-1808.

[13] Mohagheghi AA, Khan T, Meadows TH, et al. In vivo gastrocnemius muscle fascicle length in children with and without diplegic cerebral palsy [J]. Developmental Medicine and Child Neurology, 2008, 50(1): 44-50.

[14] 劉芳，蒲傳強(qiáng)，時(shí)宵冰. 肌炎特異性自身抗體在多發(fā)性肌炎/皮肌炎及其他神經(jīng)肌肉疾病的表達(dá) [J]. 臨床神經(jīng)病學(xué)雜志, 2009 (3): 219-220.

[15] Zhou Yongjin, Zheng Yongping. Estimation of muscle fiber orientation in ultrasound images using revoting hough transform (RVHT) [J]. Ultrasound In Medicine And Biology, 2008, 34(9): 1474-1481.

[16] Zhao Heng, Zhang Liqun. Automatic tracking of muscle fascicles in ultrasound images using localized radon transform [J]. IEEE Transactions on Biomedical Engineering, 2011, 58(7): 2094-2101.

[17] Rana M, Hamarneh G, Wakeling JM. Automated tracking of muscle fascicle orientation in b-mode ultrasound images [J]. Journal of Biomechanics, 2009, 42(13): 2068-2073.

[18] Zhou Yongjin, Li Jizhou, Zhou Guangquan, et al. Dynamic measurement of pennation angle of gastrocnemius muscles during contractions based on ultrasound imaging [J]. Biomedical Engineering Online, 2012, 11(17):63.

[19] Zhou Yongjin, Zheng Yongpin. Enhancement of muscle fibers in ultrasound images using Gabor filters[C]// 2009 IEEE International Ultrasonics Symposium (IUS). Roma:IEEE, 2009: 2296-2299.

[20] 牟鍇鈺，韋明，郭建平. 基于小波和卡爾曼平滑的事件相關(guān)電位單次提取 [J]. 中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào), 2012, 31(2): 167-174.

Automatic Tracking of Muscle Fascicle Orientation in Ultrasound Images: A Novel Approach Based on Local Radon Transform and Kalman Filter

Wen Huiying1Yang Xiaojuan1Guo Yanrong1Zhang Shuai2Zhang Zhiguo2Wang Tianfu1Chen Siping1Chen Xin1*

1(DepartmentofBiomedicalEngineering,SchoolofMedicine,ShenzhenUniversity,National-RegionalKeyTechnologyEngineeringLaboratoryforMedicalUltrasound,GuangdongKeyLaboratoryforBiomedicalMeasurementandUltrasoundImaging,Shenzhen518060,Guangdong,China)2(DepartmentofElectricalandElectronicEngineering,TheUniversityofHongKong,HongKong999077,China)

Ultrasound imaging is a noninvasive, real-time and convenient imaging method, which has been widely applied in analysis of human muscle contraction. So far manual operation is commonly used in clinical applications to extract morphological parameters from ultrasound image, however, may bring outcomes such as subjective, time-consuming and low repetitive. In this paper, an automatic method is proposed to calculate and track the fascicle orientation in sequences of ultrasound images. We constrained mutual information-based free-form deformation algorithm to automatically detect reading of interest. Then the detected results were projected based on local Radon transform. Then the Kalman filter image processing technique was applied to predict and track the next frame change of fascicle orientation. Six subjects participated in the walking experiment. We collected the dynamic ultrasound images from the medial gastrocnemius muscle during the process and calculated the muscle fascicle orientation using the proposed method. Experimental results showed that the proposed automatic detection method achieved better robustness. Compared with the manual measurement, the average correlation coefficient was 0.92±0.02, and the average error was 0.30°±0.62°. To some extent, the proposed method can replace the manual measurement, presenting a significant advantage especially when dealing with a large number of ultrasound data.

ultrasonography; fascicle orientation; radon transform; Kalman filter

10.3969/j.issn.0258-8021. 2016. 02.003

2015-12-26， 錄用日期:2015-11-26

國(guó)家自然科學(xué)基金(81000637)

R318

A

0258-8021(2016) 02-0141-07

*通信作者(Corresponding author)， E-mail: chenxin@szu.edu.cn