本發(fā)明涉及圖像融合和圖像序列監(jiān)測(cè)技術(shù)領(lǐng)域，具體涉及一種非特征的3d圖像快速剛性配準(zhǔn)方法。

背景技術(shù)：

3d圖像在應(yīng)用中越來越受歡迎。例如在mri(磁共振成像)中，人們習(xí)慣于從3d物體中獲取一組3d圖像來研究3d物體(如：患者的頭部)的生物機(jī)制。為了合理地去比較同一個(gè)對(duì)象在不同時(shí)刻兩個(gè)3d圖像的差異，需要事先對(duì)其進(jìn)行配準(zhǔn)(ir)。圖像配準(zhǔn)技術(shù)廣泛的應(yīng)用在醫(yī)療成像，目標(biāo)跟蹤，遙感技術(shù)，指紋和臉部識(shí)別，圖像壓縮和視頻增強(qiáng)等場(chǎng)合。

文獻(xiàn)中存在的大部分的2d圖像配準(zhǔn)方法大致可以分為基于特征和基于灰度方法的兩類。基于特征的方法首先要考慮選取兩張圖像的兩個(gè)特征集，然后去尋找能夠很好的匹配兩個(gè)特征集的幾何變換f。常用的特征包括角點(diǎn)、lbp、surf、質(zhì)心或模板，還有圖像灰度函數(shù)得到的退化像素。由于特征提取通常是一個(gè)耗時(shí)多、隨意性強(qiáng)且具有挑戰(zhàn)性的任務(wù)，最近ir研究更多關(guān)注直接對(duì)兩幅圖像灰度值的變換f，包括基于參數(shù)的轉(zhuǎn)換，和更靈活的非參數(shù)的轉(zhuǎn)換。

一些2d圖像配準(zhǔn)方法被引入到3dslicer軟件包(http://slicer.org/)中從而推廣到3d情況下。最近比較流行的icp算法是基于點(diǎn)云數(shù)據(jù)集，點(diǎn)云數(shù)據(jù)是3d掃描設(shè)備產(chǎn)生的，主要用來代表物體的外表面形狀，而現(xiàn)實(shí)中很多的3d圖像是存放整個(gè)圖像的各個(gè)位置的灰度信息(比如mri圖像)，基于點(diǎn)云的配準(zhǔn)算法不適用這些情況。

綜上可知，現(xiàn)有的一些配準(zhǔn)方法在用于3d圖像配準(zhǔn)時(shí)，往往存在配準(zhǔn)結(jié)果不準(zhǔn)確的問題。

技術(shù)實(shí)現(xiàn)要素：

在下文中給出了關(guān)于本發(fā)明的簡(jiǎn)要概述，以便提供關(guān)于本發(fā)明的某些方面的基本理解。應(yīng)當(dāng)理解，這個(gè)概述并不是關(guān)于本發(fā)明的窮舉性概述。它并不是意圖確定本發(fā)明的關(guān)鍵或重要部分，也不是意圖限定本發(fā)明的范圍。其目的僅僅是以簡(jiǎn)化的形式給出某些概念，以此作為稍后論述的更詳細(xì)描述的前序。

鑒于此，本發(fā)明提出一種非特征的3d圖像快速剛性配準(zhǔn)方法，以至少解決現(xiàn)有配準(zhǔn)方法在用于3d圖像配準(zhǔn)時(shí)存在配準(zhǔn)結(jié)果不準(zhǔn)確的問題。

根據(jù)本發(fā)明的一個(gè)方面，提供了一種非特征的3d圖像快速剛性配準(zhǔn)方法，所述非特征的3d圖像快速剛性配準(zhǔn)方法包括：步驟一、獲得實(shí)際參考圖像zr(xi,yj,zk)和實(shí)際移動(dòng)圖像zm(xi,yj,zk)；執(zhí)行步驟二；步驟二、估計(jì)下式中的和

其中，和分別表示移動(dòng)圖像fm(x,y,z)在(x,y,z)上的一階偏導(dǎo)數(shù)f′mx(x,y,z)、f′my(x,y,z)和f′mz(x,y,z)的估計(jì)函數(shù)；(xi,yj,zk)∈ω表示參考圖像與移動(dòng)圖像對(duì)應(yīng)位置的像素點(diǎn)，i,j,k指的是各像素點(diǎn)的位置下標(biāo)，zm(xi,yj,zk)表示的是實(shí)際移動(dòng)圖像，kh是以單位圓支撐的核函數(shù)，kh(x,y,z)＝(1-x²)(1-y²)(1-z²)，h＞0是一個(gè)帶寬參數(shù)；執(zhí)行步驟三；步驟三、由θ＝(x^tx)^-1x^ty獲取參數(shù)θ的當(dāng)前估計(jì)值其中，θ＝(α,β,γ,δx,δy,δz)^t，其中，x是n³×6設(shè)計(jì)矩陣，它的第[(i-1)n²+(j-1)n+k]行元素是由確定的；y是n³維向量，它的第[(i-1)n²+(j-1)n+k]元素由zr(xi,yj,zk)-zm(xi,yj,zk)確定；根據(jù)參數(shù)θ的當(dāng)前估計(jì)值確定幾何變換函數(shù)f的當(dāng)前估計(jì)值其中，執(zhí)行步驟四；步驟四、判定當(dāng)前是否滿足以及若是，執(zhí)行步驟六；否則，執(zhí)行步驟五；其中，ε1和ε2為0至1之間的實(shí)數(shù)；步驟五、根據(jù)幾何變換函數(shù)f的當(dāng)前估計(jì)值對(duì)當(dāng)前zm(xi,yj,zk)進(jìn)行變換，將zm(xi,yj,zk)變換后的值記為并用更新zm(xi,yj,zk)后返回執(zhí)行步驟二；步驟六、將參數(shù)θ的當(dāng)前估計(jì)值確定為計(jì)算結(jié)果，以完成圖像配準(zhǔn)。

進(jìn)一步地，ε1＝ε2＝0.1。

進(jìn)一步地，在定義幾何變換函數(shù)f(x,y,z)時(shí)，(x,y,z)是fm(x,y,z)的像素點(diǎn)；||f(x,y,z)-(x,y,z)||小于預(yù)設(shè)實(shí)數(shù)，以使fm(f1(x,y,z),f2(x,y,z),f3(x,y,z))的泰勒展開式有效。

現(xiàn)有技術(shù)中，大都基于特征點(diǎn)和針對(duì)特定的約束條件，帶來了特征選擇耗時(shí)多，隨機(jī)性強(qiáng)，而且約束條件使用不靈活等問題。相比之下，本發(fā)明的非特征的3d圖像快速剛性配準(zhǔn)方法針對(duì)這些問題，提出直接使用圖像灰度值的無特征3d剛性配準(zhǔn)方法，該方法使用泰勒展開式和最小二乘法直接計(jì)算待配準(zhǔn)圖像的變換參數(shù)，并且使用較少的數(shù)據(jù)點(diǎn)完成快速的配準(zhǔn)。本發(fā)明的算法能夠獲得較高的精度，并且在數(shù)據(jù)減少時(shí)仍可以有效計(jì)算；而且，該方法在需要使用數(shù)據(jù)壓縮的應(yīng)用里是非常實(shí)用的，提出的方法適合于應(yīng)用在當(dāng)數(shù)據(jù)量較大或者需要快速進(jìn)行圖像配準(zhǔn)場(chǎng)合。

通過以下結(jié)合附圖對(duì)本發(fā)明的最佳實(shí)施例的詳細(xì)說明，本發(fā)明的這些以及其他優(yōu)點(diǎn)將更加明顯。

附圖說明

本發(fā)明可以通過參考下文中結(jié)合附圖所給出的描述而得到更好的理解，其中在所有附圖中使用了相同或相似的附圖標(biāo)記來表示相同或者相似的部件。所述附圖連同下面的詳細(xì)說明一起包含在本說明書中并且形成本說明書的一部分，而且用來進(jìn)一步舉例說明本發(fā)明的優(yōu)選實(shí)施例和解釋本發(fā)明的原理和優(yōu)點(diǎn)。在附圖中：

圖1是示意性地示出本發(fā)明的非特征的3d圖像快速剛性配準(zhǔn)方法的一個(gè)示例的結(jié)構(gòu)圖。

本領(lǐng)域技術(shù)人員應(yīng)當(dāng)理解，附圖中的元件僅僅是為了簡(jiǎn)單和清楚起見而示出的，而且不一定是按比例繪制的。例如，附圖中某些元件的尺寸可能相對(duì)于其他元件放大了，以便有助于提高對(duì)本發(fā)明實(shí)施例的理解。

具體實(shí)施方式

在下文中將結(jié)合附圖對(duì)本發(fā)明的示范性實(shí)施例進(jìn)行描述。為了清楚和簡(jiǎn)明起見，在說明書中并未描述實(shí)際實(shí)施方式的所有特征。然而，應(yīng)該了解，在開發(fā)任何這種實(shí)際實(shí)施例的過程中必須做出很多特定于實(shí)施方式的決定，以便實(shí)現(xiàn)開發(fā)人員的具體目標(biāo)，例如，符合與系統(tǒng)及業(yè)務(wù)相關(guān)的那些限制條件，并且這些限制條件可能會(huì)隨著實(shí)施方式的不同而有所改變。此外，還應(yīng)該了解，雖然開發(fā)工作有可能是非常復(fù)雜和費(fèi)時(shí)的，但對(duì)得益于本公開內(nèi)容的本領(lǐng)域技術(shù)人員來說，這種開發(fā)工作僅僅是例行的任務(wù)。

在此，還需要說明的一點(diǎn)是，為了避免因不必要的細(xì)節(jié)而模糊了本發(fā)明，在附圖中僅僅示出了與根據(jù)本發(fā)明的方案密切相關(guān)的裝置結(jié)構(gòu)和/或處理步驟，而省略了與本發(fā)明關(guān)系不大的其他細(xì)節(jié)。

本發(fā)明的實(shí)施例提供了一種非特征的3d圖像快速剛性配準(zhǔn)方法，所述非特征的3d圖像快速剛性配準(zhǔn)方法包括：步驟一、獲得實(shí)際參考圖像zr(xi,yj,zk)和實(shí)際移動(dòng)圖像zm(xi,yj,zk)；執(zhí)行步驟二；步驟二、估計(jì)下式中的和

其中，和分別表示移動(dòng)圖像fm(x,y,z)在(x,y,z)上的一階偏導(dǎo)數(shù)f′mx(x,y,z)、f′my(x,y,z)和f′mz(x,y,z)的估計(jì)函數(shù)；(xi,yj,zk)∈ω表示參考圖像與移動(dòng)圖像對(duì)應(yīng)位置的像素點(diǎn)，i,j,k指的是各像素點(diǎn)的位置下標(biāo)，zm(xi,yj,zk)表示的是實(shí)際移動(dòng)圖像，kh是以單位圓支撐的核函數(shù)，kh(x,y,z)＝(1-x²)(1-y²)(1-z²)，h＞0是一個(gè)帶寬參數(shù)；執(zhí)行步驟三；步驟三、由θ＝(x^tx)^-1x^ty獲取參數(shù)θ的當(dāng)前估計(jì)值其中，θ＝(α,β,γ,δx,δy,δz)^t，其中，x是n³×6設(shè)計(jì)矩陣，它的第[(i-1)n²+(j-1)n+k]行元素是由確定的；y是n³維向量，它的第[(i-1)n²+(j-1)n+k]元素由zr(xi,yj,zk)-zm(xi,yj,zk)確定；根據(jù)參數(shù)θ的當(dāng)前估計(jì)值確定幾何變換函數(shù)f的當(dāng)前估計(jì)值其中，執(zhí)行步驟四；步驟四、判定當(dāng)前是否滿足以及若是，執(zhí)行步驟六；否則，執(zhí)行步驟五；其中，ε1和ε2為0至1之間的實(shí)數(shù)；步驟五、根據(jù)幾何變換函數(shù)f的當(dāng)前估計(jì)值對(duì)當(dāng)前zm(xi,yj,zk)進(jìn)行變換，將zm(xi,yj,zk)變換后的值記為并用更新zm(xi,yj,zk)后返回執(zhí)行步驟二；步驟六、將參數(shù)θ的當(dāng)前估計(jì)值確定為計(jì)算結(jié)果，以完成圖像配準(zhǔn)。

圖1給出了本發(fā)明的非特征的3d圖像快速剛性配準(zhǔn)方法的一個(gè)示例性處理。

如圖1所示，該方法開始后，首先執(zhí)行步驟一。

在步驟一中，獲得實(shí)際參考圖像zr(xi,yj,zk)和實(shí)際移動(dòng)圖像zm(xi,yj,zk)，然后，執(zhí)行步驟二。

在步驟二中，估計(jì)公式一一至公式一三中的和然后，執(zhí)行步驟三。

公式一一：

公式一二：

公式一三：

其中，和分別表示移動(dòng)圖像fm(x,y,z)在(x,y,z)上的一階偏導(dǎo)數(shù)f′mx(x,y,z)、f′my(x,y,z)和f′mz(x,y,z)的估計(jì)函數(shù)；(xi,yj,zk)∈ω表示參考圖像與移動(dòng)圖像對(duì)應(yīng)位置的像素點(diǎn)，i,j,k指的是各像素點(diǎn)的位置下標(biāo)，zm(xi,yj,zk)表示的是實(shí)際移動(dòng)圖像，kh是以單位圓支撐的核函數(shù)，kh(x,y,z)＝(1-x²)(1-y²)(1-z²)，h＞0是一個(gè)帶寬參數(shù)。然后，執(zhí)行步驟三。

在步驟三中，由公式二獲取參數(shù)θ的當(dāng)前估計(jì)值

公式二：θ＝(x^tx)^-1x^ty。

其中，θ＝(α，β，γ，δx，δy，δz)^t，θ是幾何變換函數(shù)f的參數(shù)(即幾何變換參數(shù))。

其中，在參數(shù)θ中，α、β和γ分別表示沿x軸、y軸和z軸的旋轉(zhuǎn)角度，而δx、δy和δz分別表示沿著x軸、y軸和z軸這三個(gè)軸的平移量。

其中，x是n³×6設(shè)計(jì)矩陣，它的第[(i-1)n²+(j-1)n+k]行元素是由公式三確定的。

公式三：

此外，y是n³維向量，它的第[(i-1)n²+(j-1)n+k]元素由zr(xi,yj,zk)-zm(xi,yj,zk)確定

這樣，可以根據(jù)參數(shù)θ的當(dāng)前估計(jì)值確定幾何變換函數(shù)f的當(dāng)前估計(jì)值幾何變換函數(shù)f可由公式四一至公式四四確定，相應(yīng)地，幾何變換函數(shù)f的估計(jì)值可由公式五一至公式五四確定。然后，執(zhí)行步驟四。

公式四一：f(x,y,z)＝(f1(x,y,z),f2(x,y,z),f3(x,y,z))

公式四二：f1(x,y,z)＝x(cosαcosβ)+y(sinαcosβ)-z(sinβ)+δx

公式四三：

f2(x,y,z)＝x(cosαsinβsinγ-sinαcosγ)+y(sinαsinβsinγ+cosαcosγ)+z(cosβsinγ)+δy

公式四四：

f3(x,y,z)＝x(cosαsinβcosγ+sinαsinγ)+y(sinαsinβcosγ-cosαsinγ)+z(cosβcosγ)+δz

公式五一：

公式五二：

公式五三：

公式五四：

在步驟四中，判定當(dāng)前是否滿足以及若是(即和同時(shí)滿足時(shí))，執(zhí)行步驟六；否則(即和中任一個(gè)未滿足，或二者都不滿足時(shí))，執(zhí)行步驟五。

其中，ε1和ε2為0至1之間的實(shí)數(shù)，例如，ε1＝ε2＝0.1。

在步驟五中，根據(jù)幾何變換函數(shù)f的當(dāng)前估計(jì)值對(duì)當(dāng)前zm(xi,yj,zk)進(jìn)行變換，將zm(xi,yj,zk)變換后的值記為并用更新zm(xi,yj,zk)后返回執(zhí)行步驟二。

步驟六、將參數(shù)θ的當(dāng)前估計(jì)值確定為計(jì)算結(jié)果，以完成圖像配準(zhǔn)。

其中，在定義幾何變換函數(shù)f(x,y,z)時(shí)，(x,y,z)是fm(x,y,z)的像素點(diǎn)，而不是fr(x,y,z)中的；同樣，這個(gè)討論同樣是基于||f(x,y,z)-(x,y,z)||小于預(yù)設(shè)實(shí)數(shù)，以使fm(f1(x,y,z),f2(x,y,z),f3(x,y,z))的泰勒展開式有效；在實(shí)際中，從fr(x,y,z)中獲取像素點(diǎn)(x,y,z)要更合理。

優(yōu)選實(shí)施例

該方法適用于一個(gè)圖像到另一圖像間的幾何剛性變換的情況，即同一幅圖中的任意兩個(gè)像素間的歐氏距離變換后不變。剛體變換在實(shí)踐中是最為常見的。本優(yōu)選實(shí)施例采用泰勒展開式和最小二乘法，通過這樣的方法，可以很好地得出參數(shù)θ＝(α,β,γ,δx,δy,δz)^t的估計(jì)值。這使得可以在3d圖像采集數(shù)據(jù)大幅減少的情況下進(jìn)行有效計(jì)算。

在數(shù)學(xué)上，3d圖像配準(zhǔn)問題描述如下：參考3d圖像設(shè)為fr(x,y,z)(即上文所述的參考圖像)，移動(dòng)3d圖像設(shè)為fm(x,y,z)(即上文所述的移動(dòng)圖像)。3d圖像配準(zhǔn)的主要目的是為了找到一個(gè)令fm(f(x,y,z))盡可能逼近fr(x,y,z)的幾何變換f(x,y,z)＝(f1(x,y,z),f2(x,y,z),f3(x,y,z))。對(duì)剛性變換來說，幾何變換核心參數(shù)包括3d旋轉(zhuǎn)、3d平移以及尺度縮放。在該優(yōu)選實(shí)施例中假設(shè)沒有尺度縮放，即尺度縮放系數(shù)為1，這也與大部分實(shí)際應(yīng)用相符合。

假設(shè)α、β和γ分別表示沿著x軸、y軸和z軸的旋轉(zhuǎn)角。(δx,δy,δz)是沿著三個(gè)軸(即x軸、y軸和z軸)的平移參數(shù)，則3d剛體變換可用下面的表達(dá)式(1)來描述：

實(shí)際上，真實(shí)的圖像fr(x,y,z)和fm(x,y,z)通常是不存在的。實(shí)際采集的圖像含有噪聲，由此可以建立模型如下:

在式(2)中，i＝1,2,…,n，j＝1,2,…,n，k＝1,2,…,n，其中n為正整數(shù)。例如，n可以為1000或10000，或者n也可以為無窮大整數(shù)。

其中的(xi,yj,zk)∈ω表示的是對(duì)應(yīng)位置的像素；zr(xi,yj,zk)表示實(shí)際的參考圖像，而zm(xi,yj,zk)表示實(shí)際的移動(dòng)圖像；εr(xi,yj,zk)表示zr(xi,yj,zk)中所含的噪聲，而εm(xi,yj,zk)表示zm(xi,yj,zk)中所含的噪聲。

因此，為了解決圖像的配準(zhǔn)問題，需從實(shí)際移動(dòng)圖像zm(xi,yj,zk)和實(shí)際參考圖像zr(xi,yj,zk)中估計(jì)出式(1)中的(α,β,γ)和(δx,δy,δz)這六個(gè)參數(shù)值，以確保zm(f(xi,yj,zk))盡可能地逼近zr(xi,yj,zk)，可以按照下面描述的過程來估計(jì)(α,β,γ)和(δx,δy,δz)。

首先，修改幾何變換f(x,y,z)＝(f1(x,y,z),f2(x,y,z),f3(x,y,z))為下式表達(dá)形式：

其中，b(x,y,z)、c(x,y,z)和d(x,y,z)分別表示在x軸、y軸和z軸上的變化量，b(x,y,z)＝f1(x,y,z)-x，c(x,y,z)＝f2(x,y,z)-y，d(x,y,z)＝f3(x,y,z)-z。因此，對(duì)f(x,y,z)估計(jì)，相當(dāng)于對(duì)(b(x,y,z),c(x,y,z),d(x,y,z))估計(jì)。假若(b(x,y,z),c(x,y,z),d(x,y,z))很小，而且fm(x,y,z)有在(x,y,z)上的一階偏導(dǎo)數(shù)。由泰勒展開式，可以得到:

上式中，f′mx(x,y,z)、f′my(x,y,z)和f′mz(x,y,z)分別表示fm(x,y,z)對(duì)x、y和z的一階偏導(dǎo)數(shù)，||·||表示歐幾里得范數(shù)。

為了計(jì)算于(b(x,y,z),c(x,y,z),d(x,y,z))的估計(jì)需令下式的計(jì)算誤差盡可能?。?/p>

fr(x,y,z)-[fm(x,y,z)+f′mx(x,y,z)b(x,y,z)+f′my(x,y,z)c(x,y,z)+fm'z(x,y,z)d(x,y,z)]。

然而，事實(shí)上，fr(x,y,z)、fm(x,y,z)、fm'x(x,y,z)、fm'y(x,y,z)和fm'z(x,y,z)都是不存在的。只有式(2)中定義的有噪聲的3d圖像{zr(xi,yj,zk)}和{zm(xi,yj,zk)}是可以得到的。

在本優(yōu)選實(shí)施例中，使用統(tǒng)計(jì)回歸分析中的最小二乘法(ls)對(duì)上述變換中的六個(gè)參數(shù)進(jìn)行估計(jì)，具體如下：首先，假設(shè)f′mx(x,y,z)、f′my(x,y,z)和f′mz(x,y,z)事先已經(jīng)由對(duì)應(yīng)的估計(jì)函數(shù)和估計(jì)出來，具體估計(jì)方法見公式(6)。于是，b(x,y,z),c(x,y,z),d(x,y,z)中的參數(shù)θ＝(α,β,γ,δx,δy,δz)^t便可由最小化問題的解決方案來估計(jì)近似：

其中，kh是一個(gè)以單位圓支撐的核函數(shù)，h＞0是一個(gè)帶寬參數(shù)。假設(shè)參數(shù)α、β和γ都很小，那么sinα≈α，sinβ≈β，sinγ≈γ，cosα≈1，cosβ≈1，cosγ≈1。在這種條件下，由式(4)得出的θ的最小二乘法估計(jì)值由下式計(jì)算：

θ＝(x^tx)^-1x^ty(5)

其中x是一個(gè)n³×6設(shè)計(jì)矩陣，它的第[(i-1)n²+(j-1)n+k]行元素是由下式確定的：

此外，y是一個(gè)n³維向量，它的第[(i-1)n²+(j-1)n+k]元素由zr(xi,yj,zk)-zm(xi,yj,zk)確定。

為了計(jì)算式(4)和式(5)所定義的參數(shù)估計(jì)值，可根據(jù)式(6)得到和的估計(jì)值：

核函數(shù)例如采用kh(x,y,z)＝(1-x²)(1-y²)(1-z²)。

只有在θ＝(α,β,γ,δx,δy,δz)^t中的所有參數(shù)值很小時(shí)，幾何變換f(x,y,z)才能被正確的估計(jì)。對(duì)于3d移動(dòng)圖像比參考圖像變換較大的時(shí)候，可以使用迭代運(yùn)算。由此可以通過上文結(jié)合圖1描述的步驟一至步驟六完成迭代計(jì)算，獲得幾何變換參數(shù)，進(jìn)而完成圖像配準(zhǔn)。

在該3d圖像快速剛性配準(zhǔn)方法中，通過使用3d圖像灰度值建立非參數(shù)變換統(tǒng)計(jì)模型(例如)，突破已有方法必須滿足特定的條件、需要預(yù)先設(shè)定參數(shù)等局限性。

盡管根據(jù)有限數(shù)量的實(shí)施例描述了本發(fā)明，但是受益于上面的描述，本技術(shù)領(lǐng)域內(nèi)的技術(shù)人員明白，在由此描述的本發(fā)明的范圍內(nèi)，可以設(shè)想其它實(shí)施例。此外，應(yīng)當(dāng)注意，本說明書中使用的語言主要是為了可讀性和教導(dǎo)的目的而選擇的，而不是為了解釋或者限定本發(fā)明的主題而選擇的。因此，在不偏離所附權(quán)利要求書的范圍和精神的情況下，對(duì)于本技術(shù)領(lǐng)域的普通技術(shù)人員來說許多修改和變更都是顯而易見的。對(duì)于本發(fā)明的范圍，對(duì)本發(fā)明所做的公開是說明性的，而非限制性的，本發(fā)明的范圍由所附權(quán)利要求書限定。