袁桂霞

摘 要： 在醫(yī)學圖像處理領(lǐng)域中，醫(yī)學圖像配準技術(shù)極其重要，其價值體現(xiàn)在臨床醫(yī)學中對圖像處理技術(shù)的應用。在解決多模態(tài)圖像配準的相關(guān)問題時，基于互信息方法的應用最廣泛，但在某些特定的應用中該方法受到的約束仍然較多。針對這一情況，提出一種新的醫(yī)學圖像配準算法，模態(tài)變換的引入作為此研究算法的基礎(chǔ)，之后新的馬爾可夫能量函數(shù)則根據(jù)兩幅通過模態(tài)變換后的圖像矩陣以及原配準圖像得以構(gòu)建。同時，為了優(yōu)化能量函數(shù)引入了一種改進的梯度下降算法，從而得到配準結(jié)果。最后，運用不同的醫(yī)學圖像進行配準實驗來驗證該算法，通過實驗證明該配準算法具有良好的有效性及抗噪性能。

關(guān)鍵詞： 圖像配準； 模態(tài)變換； 馬爾可夫隨機場； 梯度下降算法； 非剛體配準； 圖像矩陣

中圖分類號： TN911.73?34； TP311.17 文獻標識碼： A 文章編號： 1004?373X（2018）01?0057?05

Abstract： In the field of medical image processing， the medical image registration technology is extremely important. Its value is embodied in the application of image processing technology in clinical medicine. For solving the relevant issues of multimodal image registration， the method based on mutual information is most?widely used， but constrained in some specific applications. In view of this situation， a new medical image registration algorithm is proposed. The modal transformation is introduced as the basis of the studied algorithm. After that， the new Markov energy function is constructed according to the two image matrixes after modal transformation and original registration image. In order to optimize the energy function， an improved gradient descent algorithm is introduced to get the registration results. The algorithm was verified with registration experiments of different medical images. The experimental results show that the registration algorithm has perfect effectiveness and high anti?noise performance.

Keywords： image registration； modal transformation； Markov random field； gradient descent algorithm； non?rigid registration； image matrix

0 引 言

圖像配準是將圖像信息進行匯總并對齊圖像的過程，多模態(tài)圖像配準即為一種將多模態(tài)圖像進行空間幾何關(guān)系對齊的技術(shù)。近年來，將臨床醫(yī)學和圖像處理技術(shù)相結(jié)合，更加促進了醫(yī)學圖像配準技術(shù)的發(fā)展。臨床上的大多數(shù)需求可以被目前已有的算法滿足，但是考慮到圖像個體之間內(nèi)容的差異性比較大，而且包含的信息也比較多，所以現(xiàn)今的各種配準算法中都有一定的限制因素存在，一般只能在特定情況下滿足配準?？紤]到在臨床醫(yī)學方面圖像配準具有廣泛的應用，而每一種不同的應用均需要不同的形變場，因此關(guān)于醫(yī)學圖像配準的問題并沒有一種具有普適性的算法解決。針對這一情況，本文提出了一種新的醫(yī)學圖像配準算法。

1 MRF隨機場的組成

在圖像分析中，馬爾可夫隨機場的問題可以視為一個建模問題，具體到圖像配準中即構(gòu)造一個有待優(yōu)化的能量函數(shù)問題。基于馬爾可夫隨機場構(gòu)造的能量函數(shù)（或稱為目標函數(shù)）不僅和其相鄰的像素點有關(guān)，而且還需要像素點之間的其他信息，這些信息可以視為一個約束項。一個能量函數(shù)可以從以下步驟推導得出。

定義[I（x）]和[J（x）]是二維空間中兩幅大小相同的圖像，[T（·）]是每一個像素點[x]從固定圖像[I（x）]到[J（x）]的一個變換模型，對于兩幅給定的圖像來說，定義一個觀測模型如下所示：

[I（x）=JT（x）+n（x）] （1）

式中：[n（x）]是一個高斯白噪聲模型，其對應的噪聲模型數(shù)學期望為[x～（0，σ2）]，[σ]是噪聲模型的標準差。

將兩幅圖像和變換模型視為隨機變量，[s]和[t]表示像素位置，則可以通過圖[?=（v，ε）]表示馬爾可夫隨機場模型中離散標記的問題。定義[?=（1，2，…，L）]為圖像的一個標號集合，則所有的標號集合可以定義為：

[X=xss∈V] （2）

式中：[xs]是[?]中的一個元素。給定一個包含平移、翻轉(zhuǎn)、旋轉(zhuǎn)和在[x，y]方向上的尺度變換的向量[θ，]可以得到標號集合的馬爾可夫隨機場模型的能量函數(shù)：

[Emrf（Xθ）=Edata+Esmooth] （3）endprint

1.1 數(shù)據(jù)項

式（3）反映了兩個相對應的像素點及其鄰域像素點之間的關(guān)系。通常來說，在整個定義域[Ω]中，構(gòu)造一個解決配準問題的模型是定義一個數(shù)據(jù)項。這個數(shù)據(jù)項能夠表示在固定圖像和浮動圖像之間的距離函數(shù)。根據(jù)式（1）和式（2），圖像[I（x）]的噪聲模型的期望是[JT（x）]，方差是[σ2，]因此馬爾可夫隨機場模型的數(shù)據(jù)項可以表示為：

[Edata=s∈Vθs（xs）=12σ2s∈VI（xs）-JT（x）2] （4）

參數(shù)[σ]是作為測度項的調(diào)節(jié)參數(shù)引入的，它的具體值取決于人為選定的噪聲模型以及實際的應用。盡管只采用數(shù)據(jù)項進行優(yōu)化的方法在某些特定應用中是可行的，但是這并非一種十分通用的方法。這一方法的局限性體現(xiàn)在：某些應用中方程的解不存在；某些情況下對該問題的解不惟一；某些問題中，方程的解對原始數(shù)據(jù)的依賴不具有連續(xù)性。

1.2 平滑項

上述提到的所有限制都可以歸結(jié)成一個不適定問題（ill?posed problem），這一問題最直接的解決方法是采用一個平滑項做約束。對于低層級的馬爾可夫隨機場問題，這一項可以視為一個歸一化項（regularization）。它假定了在相同的鄰域系統(tǒng)中相鄰像素點的物理特性是不會突然發(fā)生變化的。因此，對于一個鄰域系統(tǒng)[Ni]而言，其平滑項為：

[Esmooth=λs，t∈Niθstxs，xt=λs∈Vt∈Niθstxs，xt] （5）

式中：[λ]表示一個權(quán)重值，其值通常根據(jù)不同的配準形式以及配準所采用的不同函數(shù)[θst]來選取不同的常數(shù)值。函數(shù)[θst]的選擇有很多種，比如分段函數(shù)、多項式函數(shù)、對數(shù)函數(shù)等。鑒于這些函數(shù)在某些情況下會加大計算的復雜度，本文采用一種線性函數(shù)來簡化優(yōu)化過程。采用一種間斷自適應（Discontinuity Adaptive， DA）模型來構(gòu)造函數(shù)的平滑項，而且函數(shù)[θ（x）]需要滿足以下方程式：

[limx→∞θ（x）=limx→∞2xf（x）=C] （6）

式中：[C]是一個正整數(shù)；函數(shù)[f（x）]是一個在定義域中的偶函數(shù)。一般來說，馬爾可夫隨機場模型通常是指二階的隨機場模型。對于一個給定的概率模型，馬爾可夫隨機場問題可以看成是一個解決最大后驗概率的問題。在上述討論條件不變的情況下將固定圖像、浮動圖像和變換模型視為隨機變量，則根據(jù)三者之間的統(tǒng)計相關(guān)性可得：

[pI（xs）T，J（xt）=pIJT（xt）] （7）

但是對于一個隨機變量[X，]上述的分布往往難以求得，此時，可以考慮運用Hammersley?Clifford 定理將馬爾可夫隨機場模型轉(zhuǎn)換為吉布斯分布：

[p（Xθ）∝e-E（Xθ）] （8）

根據(jù)式（8）可得，將馬爾可夫的最大后驗概率的相關(guān)約束問題換一種思路進行考慮，即用最小化能量函數(shù)[minXE（Xθ）]的相關(guān)約束問題來代替則更容易解決。

2 需求分析和系統(tǒng)設計

為了更大限度地提高醫(yī)學圖像配準的精度，本文在對目標函數(shù)進行優(yōu)化之前引入了一個基于直方圖計算的模態(tài)變換方法。

假設有兩幅不同模態(tài)的圖像[I（x）]和[J（x）]已經(jīng)歸一化，同時[I（x）]，[J（x）]∈[0，1]，[x]為像素點位置，通過對每一個像素點進行一個遍歷，每遍歷一個點就更新一次聯(lián)合直方圖，直到最后一個點遍歷完，得到圖像[I]和[J]的聯(lián)合直方圖[H（I，J）]。

[HI（x）N，J（x）N=HI（x）N，J（x）N+1] （9）

式中：[N]是聯(lián)合直方圖的樣條數(shù)；[x]是當前像素點的位置；[c]表示小于[c]的最大正整數(shù)值；聯(lián)合直方圖[H（i，j）]是一個二維矩陣；灰度值[i]在圖[I]以及灰度值[j]在圖像[J]中相關(guān)點具體的數(shù)目用[H（i，j）]來表示。

運用式（10）和式（11）能夠?qū)D像[I]和圖像[J]的灰度值表示形式的圖像矩陣轉(zhuǎn)換得大致相似。通過迭代過程，搜索圖像[I]在圖像[J]中重疊次數(shù)最多的一個像素點。在遍歷了所有像素點之后，得到一個新的圖像矩陣，也即經(jīng)過模態(tài)變換后的圖像矩陣。同理，可以由圖像[J]得到[JT。]

[IT（x）=argmaxjHI（x）N，jN] （10）

[JT（x）=argmaxiHiNN，J（x）N] （11）

圖1給出了利用MRI?T1序列腦部圖像和CT圖像進行模態(tài)變換的結(jié)果圖。圖1c）和圖1d）是進行模態(tài)變換前后的聯(lián)合直方圖，從圖中可以明顯看到在進行模態(tài)變換前，聯(lián)合直方圖上的點分布比較松散，變換之后的點分布緊湊，且向圖中的斜對角線靠攏，這表明模態(tài)變換的方法減少了圖像簇的分散性，在某種程度上可以提高圖像的配準精度。

3 基于MRF的模態(tài)變換多模態(tài)圖像配準

3.1 MRF能量函數(shù)

考慮到圖像[I]以及圖像[J，]在圖像的定義域[Ω=1，m×1，n]已經(jīng)給定的情況下，將得到的4個圖像矩陣結(jié)合可以得到如下數(shù)據(jù)項：

[Edata=12σ21s∈VIT（xs）-JT（x）2+12σ22s∈VI（xs）-JTT（x）2] （12）

式中：[σ1]和[σ2]是高斯白噪聲模型的標準差。在實驗中，為了簡化計算過程，取[σ1=σ2=σ]。通常來說，數(shù)據(jù)項來自于對距離函數(shù)的定義，它通常處理的是大多數(shù)的線性問題。為了提高配準精度以及配準過程的穩(wěn)定性，本文引入了另外一項表達式表示數(shù)據(jù)項，更多的細節(jié)信息的獲得來自于之前轉(zhuǎn)換后的圖像矩陣[IT]和[JT。]

能量函數(shù)只和局部形變有關(guān)，但是在優(yōu)化過程中仍然有可能會產(chǎn)生局部極值的情況，因此本文選擇了一個基本函數(shù)作為平滑項的勢能函數(shù)：

[θ（x）=ln（1+x2）] （13）endprint

[f（x）=11+x2] （14）

式中：函數(shù)[θ（x）]和[f（x）]均滿足表達式（6）的條件。由于在形變過程中，平滑項更多的是和圖像像素的鄰域系統(tǒng)有關(guān)，同時結(jié)合表達式（5）、（12）和（13），得到平滑項的相關(guān)計算表達式如下：

[Esmooth=λs∈Vt∈Niln1+xs-xt2] （15）

平滑項的引入能夠處理非線性問題，同時平滑形變場，從而可以在優(yōu)化過程中減少極值對配準結(jié)果的影響。

在得到數(shù)據(jù)項和平滑項之后，根據(jù)馬爾可夫隨機場模型理論可以得到一個新的能量函數(shù)模型：

[Emrf=Edata+Esmooth=12σ21s∈VIT（xs）-JT（x）2+ 12σ22s∈VI（xs）-JTT（x）2+λs∈Vt∈Niln1+xs-xt2] （16）

值得注意的一點是，當圖像完全配準之后，馬爾可夫隨機場的能量函數(shù)值[Emrf]將會最小。通常來說，一個馬爾可夫隨機場模型指的是一個二階模型，在計算過程中需要求像素點相鄰的八個像素點（一階系統(tǒng)中，斜對角線上的像素點不考慮在內(nèi)）。由于二階鄰域系統(tǒng)計算了更多的圖像信息，從而可以得到更精確的配準結(jié)果，本文中的實驗均采用二階鄰域系統(tǒng)。

3.2 優(yōu)化算法

對于醫(yī)學圖像配準來說，優(yōu)化方法主要分為兩大類：參數(shù)化和非參數(shù)化配準。因為非參數(shù)化配準廣泛應用于非剛體醫(yī)學圖像配準中，本文主要討論這一種優(yōu)化方法。優(yōu)化方法的通用模型如下：

[xi+1=xi+ai?di，i=0，1，2，…] （17）

式中：[xi]和[xi+1]分別為第[i]次和第[i+1]次迭代過程中的優(yōu)化值；[di]表示當前優(yōu)化過程的搜索方向，搜索方向既可以為正值也可以為負值。沿著搜索方向，[ai]是一個在每次迭代過程中用來控制步長的標量尺度系數(shù)。

一個最基本的優(yōu)化方法是梯度下降方法，它是眾多無約束最優(yōu)化方法中的一種。這種方法迭代簡單，而且相對應的配準結(jié)果也較為準確。但是這種算法在優(yōu)化時并不是逐次向最優(yōu)項靠近，而且由于迭代步長較小，所以迭代時間相對較長。因此，通過用圖形的子集計算互信息的差分方法來提高多模態(tài)圖像配準的計算時間。這一方法用預定義的一個衰減函數(shù)取代了常數(shù)值作為每次迭代的增益。其中衰減函數(shù)的表達式如下：

[ai=ai+Aα] （18）

式中：[a>0，A≥1，0<α≤1。]這一函數(shù)能夠更好地選擇步長值，同時也減少了優(yōu)化過程中的迭代次數(shù)。在此基礎(chǔ)上，用自適應隨機梯度算法（Adaptive Stochastic Gradient Descent，ASGD）對該算法進行改進。這一方法通過計算[di+1]和[di]內(nèi)積值來調(diào)整步長[ai]的大?。?/p>

[xi+1=xi-ai?di，i=0，1，2，…] （19）

[ai+1=ai?f-dTi?di-1+] （20）

式中：[[x]+]表示[x]和0進行比較后取較大的數(shù)，即[max（x，0）；]非線性sigmoid函數(shù)用函數(shù)[f]進行表示，該函數(shù)的一般表達式為[f（x）=11+e-xw]。在迭代過程中，若目標函數(shù)的當前值小于最后一次迭代值，則保持當前步長不變，否則當前步長將被式（20）的步長值取代。

通過將上述兩種優(yōu)化方法進行結(jié)合的方法改進優(yōu)化算法。在同一個配準實驗中，這一方法提供了一個更加靈活的選擇步長的方式，同時使得在優(yōu)化過程中減少了迭代次數(shù)，且配準輸出的結(jié)果和文獻中的算法具有可比性。算法中的步長值[ai+1]定義為：

[ai+1=ai?f-dTi?di-1+，Δ≥εα?ai，Δ<ε] （21）

式中：[α]是一個在定義域為[0.9，1]的常數(shù)；[Δ]是[ai]和[ai+1]差值的絕對值[ai-ai+1]；[ε]是一個根據(jù)經(jīng)驗而選取的值，通常它是一個小于1的常數(shù)。將梯度[di]由[ΔEmrf]代替，即得到最終的優(yōu)化算法：

[xi+1=xi-ai?ΔEmrf（xi）] （22）

由式（22）可以看到，對于能量函數(shù)以及能量函數(shù)的導數(shù)，通常需要逐點進行求解，這樣對于較大尺寸的圖像進行配準時，會相當耗時。因此，文中所做的實驗均采用一種多分辨率（multiresolution）的配準方法。這一方法的引入可以避免配準過程中遇到的局部極值，從而提高配準精度。將這一方法應用到研究中，首先對一個16×16的子網(wǎng)格像素點進行配準，然后將得到的形變場和原始圖像調(diào)整到32×32的子網(wǎng)格再進行配準，以此類推，直到配準過程中這一網(wǎng)格達到待配準圖像的原始大小，配準結(jié)束。圖2中，三條不同的曲線分別表示在每一個層級優(yōu)化過程的迭代情況。

為了驗證本文優(yōu)化方法的有效性，本文將上述方法和前人所用的研究方法進行了對比。實驗中除了優(yōu)化方法以外，所有的參數(shù)都設定為相同值。在式（18）和式（21）中，參數(shù)[A=10，a=1，α=0.6，ε=0.1。]同時，多分辨率方法設定層次為三層。實驗用腦部圖像參見圖1。迭代過程中兩種優(yōu)化方法在迭代過程中的曲線如圖2所示。從圖2a）中可以看出，Myronenko方法在前20～40次迭代時陷入了局部極值，且這一情況在優(yōu)化的最后一個層級仍然存在。而本文提出的優(yōu)化算法不僅在迭代次數(shù)上少于Myronenko方法（本文算法的迭代次數(shù)少于100次），而且在整個迭代過程中曲線較為平滑，較大程度上解決了圖2a）出現(xiàn)的局部極值問題。

4 算法驗證

為了驗證本文所提算法的可行性以及配準效率，本研究做了一系列圖像配準實驗。當配準過程開始時，需要應用一個全局性的仿射變換以增加配準的準確性。初始的仿射變換參數(shù)[x=][0，0，0，100，100，0，0]分別代表在[x，y]方向上的平移、旋轉(zhuǎn)、尺度變換以及裁剪（shearing）。之后，一個三層的多分辨率方法同時應用于固定圖像和浮動圖像，所有圖像的灰度值均歸一化為[0，1]的區(qū)間。在每一個層級的迭代過程中，都應用模態(tài)變化的方法。之后，優(yōu)化方法對已經(jīng)構(gòu)建好的能量函數(shù)進行優(yōu)化，能量函數(shù)的參數(shù)設定為[σ1=σ2=σ=1]。所有最大迭代次數(shù)都設定為200次。endprint

在第一個剛體實驗中，選取了MRI?T1序列腦部圖像和MRI?PD序列腦部圖像分別作為固定圖像和浮動圖像，可參見圖1a）和圖1b）。浮動圖像是在固定圖像的基礎(chǔ)上人為旋轉(zhuǎn)了10°，同時在[x]和[y]方向上分別有13 mm和17 mm的位移。圖像大小均為221×257。圖3給出了兩幅實驗圖像的配準后圖像以及配準前后的棋盤式疊加圖。通過棋盤圖可以看到圖像邊緣有了很好的配準效果。

實驗2中，驗證了3種不同的非剛體形變情況下，本文算法的配準效果。如圖4所示為本實驗的配準用圖及配準結(jié)果。圖中的每一列分別為固定圖像、浮動圖像、配準后圖像以及配準前后的棋盤格疊加圖。從配準前后疊加圖的效果來看，本文提出的算法輸出了準確的配準結(jié)果，圖中從上到下的圖像模態(tài)分別為MRI?T2/CT、MRIT2/MRIT1和CT/MRI?T2。

5 結(jié) 語

本文在結(jié)合MRF隨機場和模態(tài)變換的基礎(chǔ)上，構(gòu)造了一個新的待優(yōu)化能量函數(shù)。此能量函數(shù)由數(shù)據(jù)項和平滑項構(gòu)成。數(shù)據(jù)項包含了4個圖像矩陣，分別為固定圖像和浮動圖像，以及由這兩幅圖像進行模態(tài)變換后得到的2個圖像矩陣。增加兩幅圖像矩陣可以增加圖像的細節(jié)，通過獲得更多的圖像信息從而提高圖像配準的精度，平滑項選取了自然對數(shù)函數(shù)。

在優(yōu)化方法的選取上，本文改進了梯度下降算法，將兩種優(yōu)化方法進行結(jié)合，使優(yōu)化過程中的迭代步長選擇更加靈活。通過實驗驗證，此方法能夠有效地減少圖像迭代次數(shù)，大大降低了優(yōu)化過程陷入局部極值的情況。為了驗證本文算法的配準精度以及抗噪性能，本文用MRI不同序列腦部圖像和CT圖像做了一系列的配準實驗。實驗結(jié)果表明，本文提出的配準方法配準效果更為準確，同時抗噪性能更好。

參考文獻

[1] 張冉，王雷，夏威，等.2D/3D圖像配準中的相似性測度和優(yōu)化算法[J].激光與紅外，2014（1）：98?102.

ZHANG Ran， WANG Lei， XIA Wei， et al. Comparison of similarity measurement and optimization [J]. Laser & infrared， 2014（1）： 98?102.

[2] 王雷，黃晨雪，王高波.改進的MRF彩色圖像分割算法[J].計算機工程與應用，2016，52（13）：191?194.

WANG Lei， HUANG Chenxue， WANG Gaobo. Improved MRF color image segmentation algorithm [J]. Computer engineering and applications， 2016， 52（13）： 191?194.

[3] 張靜亞，王加俊.一種改進的非剛性醫(yī)學圖像配準算法[J].計算機應用研究，2015，32（4）：1261?1264.

ZHANG Jingya， WANG Jiajun. Improved non?rigid medical image registration algorithm [J]. Application research of computers， 2015， 32（4）： 1261?1264.

[4] 金斌，周偉，叢瑜，等.一種基于局部不變特征的SAR圖像配準新算法[J].哈爾濱工業(yè)大學學報，2014，46（11）：112?118.

JIN Bin， ZHOU Wei， CONG Yu， et al. A novel and efficient algorithm using local invariant feature for SAR image registration [J]. Journal of Harbin Institute of Technology， 2014， 46（11）： 112?118.

[5] 趙雪梅，李玉，趙泉華.結(jié)合高斯回歸模型和隱馬爾可夫隨機場的模糊聚類圖像分割[J].電子與信息學報，2014，36（11）：2730?2736.

ZHAO Xuemei， LI Yu， ZHAO Quanhua. Image segmentation by fuzzy clustering algorithm combining hidden Markov random field and Gaussian regression model [J]. Journal of electronics & information technology， 2014， 36（11）： 2730?2736.

[6] 徐勝軍，韓九強，何波，等.融合邊緣特征的馬爾可夫隨機場模型及分割算法[J].西安交通大學學報，2014，48（2）：14?19.

XU Shengjun， HAN Jiuqiang， HE Bo， et al. A region Markov random field model with integrated edge feature and image segmentation algorithm [J]. Journal of Xian Jiaotong University， 2014， 48（2）： 14?19.

[7] 余麗玲，徐彬鋒，金浩宇，等.基于馬爾科夫隨機場的乳腺DCE?MRI圖像序列配準[J].計算機與現(xiàn)代化，2015（11）：74?78.

YU Liling， XU Binfeng， JIN Haoyu， et al. Images registration of breast DCE?MRI by Markov random field [J]. Computer and modernization， 2015（11）： 74?78.endprint

[8] 林江，戴齊，歐陽婷雪，等.一種邊界和馬爾可夫隨機場相結(jié)合的腦MRI醫(yī)學圖像分割方法[J].中國醫(yī)學物理學雜志，2015，32（5）：717?720.

LIN Jiang， DAI Qi， OUYANG Tingxue， et al. Brain magnetic resonance image segmentation combined boundary and Markov random field [J]. Chinese journal of medical physics， 2015， 32（5）： 717?720.

[9] 郝培博，陳震，江少鋒，等.基于Demons算法的2D、3D多模態(tài)醫(yī)學圖像非剛性配準研究[J].生物醫(yī)學工程學雜志，2014（1）：161?165.

HAO Peibo， CHEN Zhen， JIANG Shaofeng， et al. Research on non?regid registration of multi?model medical image based on demons algorithm [J]. Journal of biomedical engineering， 2014（1）： 161?165.

[10] 劉曉慧，楊新鋒.基于歸一化互信息和金字塔分解優(yōu)化的圖像配準算法研究[J].微型電腦應用，2016，32（11）：19?22.

LIU Xiaohui， YANG Xinfeng. Research on image registration algorithm based on normalized mutual information and pyramid decomposition optimization [J]. Microcomputer applications， 2016， 32 （11）： 19?22.

[11] 秦緒佳，肖佳吉，陳珊，等.局部更新的分層B樣條醫(yī)學圖像非剛性配準算法[J].小型微型計算機系統(tǒng)，2016，37（10）：2338?2342.

QIN Xujia， XIAO Jiaji， CHEN Shan， et al. Local updating algorithm of hierarchical B?spline based non?rigid registration for medical images [J]. Journal of Chinese computer systems， 2016， 37（10）： 2338?2342.endprint