1、基于快速匹配算法的頸動脈超聲圖像的運動分析                     作者：徐明才，陶進緒連娟，李昕，吳學詠，汪華【摘要】  頸動脈超聲圖像中的運動信息能夠間接地反應頸動脈彈性等狀況，結合頸動脈內中膜厚度(CIMT)能夠為心腦血管疾病診斷提供定性與定量的依據。我們將改進后的金字塔快速匹配算法(modified block sum pyramid，MBSP)應用于頸動脈超聲波

2、圖像斑點跟蹤獲得運動信息。理論和實驗結果都表明，改進后的金字塔塊匹配算法能有效地減少運動跟蹤的運算量，并且有著和改進前的金字塔塊匹配算法相同的準確度。運動跟蹤結果能夠為醫生診斷心腦血管疾病起到一定的輔助作用。 【關鍵詞】  頸動脈;斑點跟蹤;金字塔算法;運動分析;平均絕對差函數Abstract：The motion information in the carotid ultrasound image can reflect the condition of carotid whether it has hardened, expanded and contracted in the

3、 right way. Then it can help doctor to diagnose cardio cerebral vascular diseases together with the carotid intima-medial thickness(CIMT).We used the modified block sum pyramid (MBSP) algorithm to analyze the motion of carotid ultrasound images. The result shows that the modified algorithm effective

4、ly reduced the computation while speckle tracking. Besides that, the MBSP algorithm has the same accuracy with the block sum pyramid. The tracking result can assist the doctor to evaluate the disk of cardio cerebral vascular diseases to a certain degree.Key words：Carotid; Speckle tracking; Block sum

5、 pyramid; Motion analysis; Mean absolute difference1 引 言心腦血管疾病與頸動脈內中膜厚度及其運動特征等有著密切的相關性1。對于超聲波圖像而言，可以通過斑點跟蹤的方法獲得運動信息。斑點作為噪聲不可避免地存在于超聲波圖像中。雖然斑點噪聲影響超聲波圖像質量，但是根據斑點噪聲的形成原理，可以通過跟蹤斑點來獲取圖像中相關目標的運動信息。為研究組織的運動特征和彈性信息，進而輔助評估心腦血管疾病提供了一種可能。斑點跟蹤算法一般包括兩個方面：匹配和搜索。匹配的準則有許多種，常用的準則有歸一化相關函數(NCCF),平均均方誤差函數(MSD)以及平均絕對差函數

6、(MAD)等。因MAD不涉及乘法等復雜運算，所以較為常用，本研究也采用該匹配準則。搜索方法也有好多種，有全搜索法2，十字搜索法3，菱形搜索法4等。此外，在這些算法基礎上出了許多快速算法，如金字塔算法5,連續淘汰法6,改進的自適應2BIT變換法7等。臨床實踐中醫生除了用肉眼定性地去觀察頸動脈的運動特征，有時需要定量測量出感興趣的頸動脈某一目標的運動信息(如：短軸切面輪廓在一個心動周期內收縮和擴張的運動狀況)。這樣就能夠更加可靠地評估心腦血管疾病和風險。由于輔助診斷對運動跟蹤的實時性要求不是很高，而對運動跟蹤的準確性要求較高。故我們選擇能達到全局最優的全搜索法和金字塔快速匹配算法的組合進行斑點運動

7、跟蹤，并且對金字塔快速匹配算法進行了改進，將其應用于頸動脈運動分析，高效準確地跟蹤出目標興趣點運動的相關信息。2 原理與方法2.1 金字塔塊匹配算法設有塊X和Y,塊大小均為：2M×2M，對塊X和Y均建立 M+1層的金字塔。從塔頂往下依次為第0，1，M層。X,Y的金字塔的第m 層分別記為Xm,Ym。第m-1層和第m層關系如方程(1)所示：Xm-1(i,j)=Xm(2i-1,2j-1)+Xm(2i-1，2j)+Xm(2i,2j-1)+Xm(2i,2j)(1)令：MADm(X,Y)=2mi=12mj=1|Xm(i,j)-Ym(i,j)|(2)可以得出結論(3)：MAD0(X，Y)MAD1(

8、X,Y)MAD2(X,Y)MADm(X,Y)MADm(X,Y)(3)兩個塊在生成了各自的金字塔之后，從頂層往下逐層MADm(X,Y)。如果發現其大于或者等于最新的MADmin時，說明本候選塊一定不是最優塊，舍棄該候選塊繼續搜索下一個候選塊，如果比較到最底層,仍然比最新的MADmin還小時，將MADmin更新為當前候選塊與塊的MADm(X,Y)。直到所有的候選塊均搜索完后，最小的MADmin的候選塊相對于參考塊的偏移量即為所要求解的運動向量。 金字塔匹配算法就是通過盡快地舍棄那些不是最優的候選塊而減少運算量的。雖然建立金字塔會增加額外的計算量，但是由于大部分的塊在比較時就被舍棄了，所以總的運算量

9、還是大幅地減少。 2.2 改進的金字塔匹配算法金字塔塊匹配算法從頂層開始往下逐層比較時，層與層間的MAD值相差很大。當層數較大時，最后一層與倒數第二層的MAD值相差更大，因此，比較時仍有不少的計算冗余。為此，希望在現在的金字塔基礎上不改變運算量，在層與層之間新增一些閾值，能夠盡早地舍棄那些不是最優的塊。為了敘述方便，以下討論均以在最后一層(第M層)與倒數第二層(第M-1層)之間增加閾值為例。兩個要比較的塊分別記為X,Y,塊X和Y的金字塔的第M-1層分別記為：XM-1，YM-1。設SM-1為第M-1層塔的所有像素點空間。將XM-1，YM-1按照相同的方式均分割為N個區域：S1M-1,S2M-1，

10、S3M-1，SkM-1，,SNM-1，1kN，見圖1。具體的分割方法還需要根據實際情況進行合理的選擇。分割出的N個區域要滿足以下兩個約束關系：Nk=1SM-1k=SM-1(4)SM-1kSM-1l=，1klN(5)圖1 像素空間分割方法Fig 1 Pixels partition method圖2 實驗所采取的分割方法Fig 2 Experiments method同時對塊X，Y的金字塔的第M層也劃分為N個區域SM1,SM2,SM3,SMk，SMN，1kN。也要滿足約束關系：Nk=1SMk=SM和SMkSMl=，1klN。SM-1k區域中的每一個點對應于SMk區域中的4個點，此對應關系和構建金

11、字塔時從M層到M-1層的4個點累加成1個點的對應關系一致。構造如下函數表達式：(0nN)MADM-1n(X，Y)=MADM-1(X，Y)-nk=1(i,j)SM-1k|XM-1(i,j)YM-1(i,j)|+nk=1(i,j)SMk |XM(i,j)YM(i,j)|(6)根據層間區域關系式(1)，上式又可以寫成(7)式：MADM-1n(X，Y)=MADM-1(X，Y)-f(n)(7)其中f(n)可寫成(8)式:f(n)=nk=1(i,j)SM-1k|XM-1(i,j)YM-1(i,j)|-1         |XM

12、(2i-1,2j-1)YM(2i-1,2j-1)|+|XM(2i-1,2j)YM(2i-1,2j)|+|XM(2i,2j-1)YM(2i,2j-1)|+|XM(2i,2j)YM(2i,2j)|(8)又據三角不等式及上下層點的對應關系知：f(n)0，且f(n)隨著n的增大單調遞減,據此得出(9)、(10)、(11)的3個結論：MADM-1(X，Y)=MADM-10(X,Y)(9)MADM-1N(X，Y)=MADM(X,Y)(10)MADM-10(X，Y)MADM-11(X，Y)MADM-1n(X，Y)MADM-1N(X，Y)(11)如果將SM-1像素空間一列一列依次從左往右分割時，所得到的有關結

13、論則和8是一致的。文獻8中所提到的逐行更新算法可看作是本研究改進算法的特例。現將改進的金字塔匹配算法(以在兩幀中跟蹤一個興趣點為例)歸納如下：記前后兩幀圖像分別為fi(x,y),fi+1(x,y)，被跟蹤點坐標為(a,b)，前一幀圖像點(a,b)所對應的塊稱為塊，搜索范圍為A。(1)相對于參考塊運動位移為(0,0)的候選塊與參考塊兩者的MAD值作為MADmin，并且確定后一幀圖像中的起始候選塊。MAD計算公式如(12)式所示：MADmin=2M-1-1x=-2M-1 2M-1-1y=-2M-1|fi(a+x,b+y)-fi+1(a+x,b+y)|(12)(2)依據金字塔算法構建當前候選塊和參考

14、塊各自的金字塔。金字塔總共假設有M+1層。沿金字塔至上而下逐層進行比較，因為這里只考慮在M-1層和M層間新增閾值，所以前M-1層的比較和原金字塔塊匹配算法一樣。當M-1層比較完后，發現MADM-1(X，Y)仍然小于最新的MADmin。此時則需要按照構造函數表達式生成一個新的中間閾值，再進行比較，其間發現新生成的閾值大于或者等于MADmin時，舍棄，直接進入步驟(3)進行下一個候選塊的匹配。否則一直進行到最后一層為止。如果其結果仍然小于MADmin時，更新MADmin。(3)按照選取的搜索算法(如全搜索法)，確認下一點被候選點及相應的候選塊。重復步驟(2)，直到A所有需要搜索的點搜索完為止。(4

15、)最后，對應于MADmin值最小的候選塊即為最優候選塊，其相對于參考塊的偏移量，即為所要跟蹤的運動矢量。2.3 利用MBSP算法分析頸動脈圖像(1)取頸動脈短軸切面一個心動周期的圖像序列。在第一幀圖像中，沿頸動脈短軸輪廓，手動大致均勻地選取出若干個興趣點，興趣點應盡量選取在整個運動序列圖像中亮度較為平穩且肉眼觀察較為明顯的斑點噪聲顆粒。(2)利用改進的金字塔跟蹤算法(MBSP)逐幀地跟蹤所有興趣點，得到各個興趣點在各幀中的運動矢量，再將各幀上各點的運動矢量按照適合人眼觀看的比例顯示在各幀圖像上。最后以動態的方式播放附有運動矢量的頸動脈圖像序列，醫生可以根據播放視頻中的運動矢量定性地掌握興趣點運

16、動大小、方向、扭轉、是否出現異常等情況，做出一些如是否可能出現斑塊硬化等初步診斷。(3)最后給出各個興趣點的收縮擴張曲線，每條曲線給出每個點隨時間的運動速度變化。設第j幀圖像中跟蹤到的N個點的分別為坐標(X1，Y1),(Xi,Yi),(XN,YN)。其均值為(X,Y)。記Vi=(Xi,Yi)-(X,Y),13 實驗與結果3.1 實驗結果實驗數據由合肥105超聲科提供，為一個心動周期的頸動脈短軸切片序列圖像，共25幀。第1、11、22幀圖像分別見圖3、4、5。圖3中小白框即為我們要跟蹤的6個興趣點，左上角為1號點，按順時針依次排列。為了表述方便，這里只選取了6個興趣點為例。分別用BSP和MBSP

17、算法對該頸動脈序列圖像進行斑點運動跟蹤，其中MBSP算法在第3層和第4層增加了3個新的閾值。將第3層的像素空間從左到右依次分割成4個相同的區域，見圖2。實驗中塊大小取為16×16，搜索范圍為8×8。圖4附上了第11幀圖像上各興趣點的運動矢量圖，兩個算法跟蹤出來的結果一樣。圖6給出了各個興趣點(從左往右，從上至下依次為1到6)隨時間變化的擴張收縮運動速度圖。從圖中我們可以看出，興趣點2、3、4的運動表現出較為一致的運動，在一個心動周期內經歷一個明顯的擴張和之后的收縮過程。但是觀察興趣點1、5、6的運動速度曲線，發現其運動起伏相對來說沒有這些規律，從而可以認為1、5、6附近區域

18、可能出現一些病變。這和醫生的診斷結果：本試驗者頸動脈左下部區域內膜形成了較為明顯的斑塊，而斑塊區域收縮運動常出現紊亂是相互吻合的。3.2 性能分析(1)準確性。分析改進后的算法和改進前的金字塔匹配算法(BSP)，容易得出前后算法運動跟蹤結果應該是一致的，差異只在于計算量上。實際跟蹤的結果也證明了這一點。圖7、8給出了前后兩種算法的跟蹤興趣點1的運動軌跡圖，圖中橫軸、縱軸分別對應于圖像的X軸和Y軸坐標軸。(2)效率。表1給出了改進前后算法在上述實驗過程中需要計算的各MAD值統計次數。計算第m層的MADm值所需要運算量(含加減法及絕對值)可以公式化9為：3.22m-1,計算MADmn運算量約為：5

19、n·(3·22m+2-1)/16。以第4層為例，如果直接采用表1 金字塔算法改進前后各MAD次數原來的金字塔算法，需要計算總次數為821 457;但采用改進后的算法需要的計算量僅為721 747。可見因為第3層與第4層3個新增閾值的引入，第4層判決過程的總體運算就減少12.11%，可以考慮將區域進行更細劃分，但是要權衡因為閾值增加而造成的比較運算量。如果在第3層等層也引入同樣機制的話，運算量將進一步減少。4 本研究將改進后的金字塔匹配算法應用于斑點運動跟蹤，實驗結果表明該算法能準確而更加有效地跟蹤出運動矢量。并且將該方法應用于頸動脈運動分析，結合頸動脈的自身特點，給出有利于輔助醫生診斷的擴張與收縮曲線等依據，為頸動脈的進一步分析奠定了一定基礎。【參考文獻】  1胡瑋,張軍,孔國喜，等. 超聲檢測頸內動脈內-中膜厚度及其形態結構與相關心腦血管疾病的關系J. 第四軍醫大學學報 ,2007,28(15):1431-1433.2Brünig M, Niehsen W. Fast full-search block matchingJ. IEEE Trans Circuits Syst Video