本發(fā)明涉及醫(yī)療圖像處理領(lǐng)域，尤其是一種能夠提高病灶或感興趣部位的配準(zhǔn)精度，有助于臨床診斷、放射治療計(jì)劃的制定和評(píng)價(jià)的基于SURF算法的多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)方法。
背景技術(shù)：
目前，隨著計(jì)算機(jī)科學(xué)技術(shù)和醫(yī)療影響工程學(xué)的快速發(fā)展，世界上出現(xiàn)了許多先進(jìn)的醫(yī)療成像設(shè)備，為臨床醫(yī)學(xué)診斷提供了多種模態(tài)的醫(yī)學(xué)圖像，這些圖像從不同方面反映了人體結(jié)構(gòu)、臟器和病變組織的不同信息。
比如CT(ComputedTomography)圖像具有較強(qiáng)的空間分辨率和幾何特性，對(duì)骨骼的成像非常清晰，它可對(duì)病灶定位提供較好的參照，但對(duì)軟組織的對(duì)比度則相對(duì)較低。MR(MagneticResonance)圖像可以清晰地反映軟組織、器官、血管等的解剖結(jié)構(gòu)，有利于確定病灶范圍，但是MR圖像對(duì)鈣化點(diǎn)不敏感，且受到磁干擾會(huì)發(fā)生幾何失真。SPEC、PET圖像能得到人體任意角度斷層面的放射性濃度分布，可以反映組織器官的代謝水平和血流狀況，對(duì)腫瘤病變呈現(xiàn)“熱點(diǎn)”，提供人體的功能信息，但是它們的分辨率差，很難得到精確的解剖結(jié)構(gòu)，也不易分辨組織、器官的邊界。由此可見，不同成像技術(shù)有著自身的優(yōu)勢(shì)也同時(shí)擁有一些局限性，這些圖像對(duì)人體同一解剖結(jié)構(gòu)所得到的形態(tài)和功能信息是互為差異、互為補(bǔ)充的。
在臨床診斷中，單一模態(tài)的圖像往往不能提供醫(yī)生所需要的足夠信息，因此，如果能將不同模態(tài)的醫(yī)學(xué)圖像進(jìn)行適當(dāng)?shù)娜诤?，使解剖信息和功能信息有機(jī)地結(jié)合起來，在一幅圖像上同時(shí)綜合地表達(dá)來自多種成像源的信息，以便醫(yī)生了解病變組織或器官的綜合情況，并做出更加準(zhǔn)確的診斷或制定出更加科學(xué)優(yōu)化的治療方案，這必將推動(dòng)現(xiàn)代醫(yī)學(xué)臨床技術(shù)的巨大進(jìn)步。
由于從醫(yī)學(xué)圖像角度來看MR與CT圖像的圖像內(nèi)容是不同的，因此我們?cè)诖税l(fā)明中引入了基于改良的SURF算法的圖像配準(zhǔn)方法來介紹如何進(jìn)行精準(zhǔn)的MR與CT圖像的配準(zhǔn)。
技術(shù)實(shí)現(xiàn)要素：
本發(fā)明的目的是針對(duì)多源圖像間的配準(zhǔn)問題，提出一種基于SURF算法的多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)方法。此種基于SURF算法的圖像配準(zhǔn)方法可以實(shí)現(xiàn)多模態(tài)醫(yī)學(xué)圖像融合過程中至關(guān)重要的一步——圖像配準(zhǔn)。此方法可實(shí)現(xiàn)精準(zhǔn)圖像配準(zhǔn)，同時(shí)相比于類似的SIFT算法，此方法使用SURF算法可以快速完成圖像配準(zhǔn)過程。
本發(fā)明的技術(shù)方案是：
一種基于SURF算法的多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)方法，它包括以下步驟：
S1、利用SURF算法提取不同模態(tài)醫(yī)學(xué)圖像中的特征點(diǎn)；
S2、確定各模態(tài)醫(yī)學(xué)圖像中特征點(diǎn)的主方向，并構(gòu)造對(duì)應(yīng)的SURF特征描述子；
S3、對(duì)特征點(diǎn)的SURF特征描述子進(jìn)行矩陣運(yùn)算，獲取轉(zhuǎn)換矩陣；
S4、使用轉(zhuǎn)換矩陣完成不同模態(tài)醫(yī)學(xué)圖像配準(zhǔn)。
本發(fā)明的步驟S1中，利用SURF算法提取的不同模態(tài)醫(yī)學(xué)圖像中的特征點(diǎn)具體為：
S1.1、根據(jù)Hessian矩陣，計(jì)算特征值α，具體步驟如下：
S1.1-1、利用下述Hessian矩陣公式計(jì)算出不同模態(tài)醫(yī)學(xué)圖像中各個(gè)像素點(diǎn)的Hessian矩陣：
其中，x表示各模態(tài)醫(yī)學(xué)圖像中各個(gè)像素點(diǎn)的灰度值，Lxx(x,σ),Lxy(x,σ),Lyy(x,σ)表示當(dāng)前像素在x，y方向上的二階偏導(dǎo)數(shù)，也即二階標(biāo)準(zhǔn)高斯函數(shù)對(duì)圖像的卷積；
S1.1-2、對(duì)于各模態(tài)醫(yī)學(xué)圖像，分別采用下屬公式計(jì)算各個(gè)像素點(diǎn)的Hessian矩陣行列式的近似值作為對(duì)應(yīng)像素點(diǎn)的特征值α：
α＝LxxLyy-(0.9Lxy)2
S1.2、根據(jù)SURF算法特性構(gòu)造高斯金字塔，根據(jù)像素點(diǎn)的特征值α是否為鄰域極大值判斷特征點(diǎn)，具體步驟為：
S1.2-1、構(gòu)造SURF高斯金字塔，前述金字塔分為若干層，每一層均作為一個(gè)頻度范圍Octave，每個(gè)Octave尺度不同的圖片，在SURF算法中圖片大小即尺寸始終保持不變；
S1.2-2、在構(gòu)造好的SURF高斯金字塔中，將經(jīng)過步驟S1.1-1中Hessian矩陣處理過的每個(gè)像素點(diǎn)與其在SURF高斯金字塔中三維領(lǐng)域的26個(gè)點(diǎn)進(jìn)行大小數(shù)值比較，若該像素點(diǎn)為這26個(gè)點(diǎn)中的最大值或最小值，則將像素點(diǎn)保留下來，作為初步的特征點(diǎn)，否則，采用三維線性插值算法得到亞像素級(jí)特征點(diǎn)。
本發(fā)明的步驟S1.2-1中，每一層Octave的獲取方式是：通過對(duì)原始圖片進(jìn)行不同尺度的高斯模糊得到的，同一Octave中的各個(gè)圖片也是通過不同高斯模糊尺度的模糊得到的。
本發(fā)明的步驟S2中，確定特征點(diǎn)主方向并構(gòu)造SURF特征描述子，提取不同模態(tài)醫(yī)學(xué)圖像的特征點(diǎn)具體為：
S2.1、根據(jù)SURF算法特性統(tǒng)計(jì)特征點(diǎn)鄰域內(nèi)的Harr小波特征，選擇最長(zhǎng)矢量的方向?yàn)樵撎卣鼽c(diǎn)的主方向，具體步驟為：
S2.1-1、統(tǒng)計(jì)任一特征點(diǎn)領(lǐng)域內(nèi)的Harr小波特征，即以特征點(diǎn)為中心，計(jì)算半徑為6s的鄰域內(nèi)，s為特征點(diǎn)所在的尺度值，統(tǒng)計(jì)60°扇形內(nèi)所有點(diǎn)在水平x和垂直y方向的Harr小波響應(yīng)總和(Harr小波邊長(zhǎng)取4s，s為特征點(diǎn)所在的尺度值，并給這些Harr小波響應(yīng)值賦上高斯權(quán)重系數(shù)，高斯權(quán)重系數(shù)由高斯模型在不同的角度和距離上確定，權(quán)重系數(shù)的和為1)；然后將60°范圍內(nèi)的響應(yīng)相加以形成特征矢量(特征值加上方向信息即形成特征矢量)，遍歷整個(gè)圓形區(qū)域，覆蓋整個(gè)360°，選擇最長(zhǎng)矢量的方向?yàn)樵撎卣鼽c(diǎn)的主方向。
遍歷所有的特征點(diǎn)，得到每個(gè)特征點(diǎn)的主方向。
S2.2、根據(jù)S2.1步驟中得到的方向構(gòu)造對(duì)應(yīng)特征點(diǎn)的SURF特征描述子，具體方法為：
S2.2-1、在特征點(diǎn)鄰域范圍內(nèi)取一個(gè)正方形框，框的邊長(zhǎng)為20s(s是該特征點(diǎn)所在的尺度值)，該框的方向即為步驟S2.1得到的主方向；
S2.2-2、把該框分為16個(gè)子區(qū)域，每個(gè)子區(qū)域統(tǒng)計(jì)25個(gè)像素的水平方向和垂直方向的Harr小波特征，前述水平和垂直方向均為相對(duì)特征點(diǎn)的主方向而言。該Harr小波特征為：水平方向灰度值之和∑dx，水平方向灰度值絕對(duì)值之和∑|dx|，垂直方向灰度值之和∑dy，以及垂直方向灰度值絕對(duì)值之和∑|dy|，對(duì)于每個(gè)特征點(diǎn)，建立64維向量作為該特征點(diǎn)的SURF特征描述子，其中列向量對(duì)應(yīng)16個(gè)子區(qū)域，行向量對(duì)應(yīng)各子區(qū)域的4個(gè)Harr小波特征參數(shù)值。
本發(fā)明的步驟S3具體為：
S3.1對(duì)于各模態(tài)醫(yī)學(xué)圖像中相應(yīng)位置的圖像以遍歷的方式計(jì)算兩張圖像所有特征點(diǎn)描述子的內(nèi)積，
即對(duì)于來自兩張圖像的對(duì)應(yīng)的兩個(gè)特征點(diǎn)，計(jì)算64維特征向量的內(nèi)積，按照數(shù)值從大到小進(jìn)行排序，得到排序后的特征點(diǎn)序列，其中數(shù)值最大者為最匹配的點(diǎn)；
S3.2、對(duì)前述排序后的特征點(diǎn)序列進(jìn)行矩陣運(yùn)算，得到配準(zhǔn)需要的3X3轉(zhuǎn)換矩陣。
本發(fā)明的步驟S3還包括以下步驟：設(shè)定閾值，在排序后的特征點(diǎn)序列中，選擇大于前述閾值的配對(duì)特征點(diǎn)進(jìn)行矩陣運(yùn)算。
本發(fā)明的閾值的設(shè)置方法為：通過對(duì)大量不同模態(tài)的醫(yī)學(xué)圖像的機(jī)器學(xué)習(xí)，設(shè)置特征點(diǎn)的評(píng)判閾值。
本發(fā)明的配準(zhǔn)過程如下：
S4.1使用S3.3計(jì)算得到的轉(zhuǎn)換矩陣對(duì)待配準(zhǔn)圖像進(jìn)行卷積，得到每個(gè)像素點(diǎn)新的坐標(biāo)，并使用雙三次插值方法計(jì)算配準(zhǔn)后的像素值。從而完成兩幅不同模態(tài)圖像的配準(zhǔn)。
本發(fā)明的有益效果：
本發(fā)明的SURF算法相比于SIFT算法，計(jì)算速率大大提高了。利用GPU等并行化技術(shù)可以使得大分辨率圖像的線性配準(zhǔn)達(dá)到實(shí)時(shí)。
具體實(shí)施方式
下面結(jié)合實(shí)施例對(duì)本發(fā)明作進(jìn)一步的說明。
一種基于SURF算法的多模態(tài)醫(yī)學(xué)圖像配準(zhǔn)方法，它包括以下步驟：
S1、利用SURF算法提取不同模態(tài)醫(yī)學(xué)圖像中的特征點(diǎn)；步驟S1中，利用SURF算法提取的不同模態(tài)醫(yī)學(xué)圖像中的特征點(diǎn)具體為：
S1.1、根據(jù)Hessian矩陣，計(jì)算特征值α，具體步驟如下：
S1.1-1、利用下述Hessian矩陣公式計(jì)算出不同模態(tài)醫(yī)學(xué)圖像中各個(gè)像素點(diǎn)的Hessian矩陣：
其中，x表示各模態(tài)醫(yī)學(xué)圖像中各個(gè)像素點(diǎn)的灰度值，Lxx(x,σ),Lxy(x,σ),Lyy(x,σ)表示當(dāng)前像素在x，y方向上的二階偏導(dǎo)數(shù)，也即二階標(biāo)準(zhǔn)高斯函數(shù)對(duì)圖像的卷積；
S1.1-2、對(duì)于各模態(tài)醫(yī)學(xué)圖像，分別采用下屬公式計(jì)算各個(gè)像素點(diǎn)的Hessian矩陣行列式的近似值作為對(duì)應(yīng)像素點(diǎn)的特征值α：
α＝LxxLyy-(0.9Lxy)2
S1.2、根據(jù)SURF算法特性構(gòu)造高斯金字塔，根據(jù)像素點(diǎn)的特征值α是否為鄰域極大值判斷特征點(diǎn)，具體步驟為：
S1.2-1、構(gòu)造SURF高斯金字塔，前述金字塔分為若干層，每一層均作為一個(gè)頻度范圍Octave，每個(gè)Octave尺度不同的圖片，在SURF算法中圖片大小即尺寸始終保持不變；(其中：每一層Octave的獲取方式是：通過對(duì)原始圖片進(jìn)行不同尺度的高斯模糊得到的，同一Octave中的各個(gè)圖片也是通過不同高斯模糊尺度的模糊得到的)
S1.2-2、在構(gòu)造好的SURF高斯金字塔中，將經(jīng)過步驟S1.1-1中Hessian矩陣處理過的每個(gè)像素點(diǎn)與其在SURF高斯金字塔中三維領(lǐng)域的26個(gè)點(diǎn)進(jìn)行大小數(shù)值比較，若該像素點(diǎn)為這26個(gè)點(diǎn)中的最大值或最小值，則將像素點(diǎn)保留下來，作為初步的特征點(diǎn)，否則，采用三維線性插值算法得到亞像素級(jí)特征點(diǎn)。
S2、確定各模態(tài)醫(yī)學(xué)圖像中特征點(diǎn)的主方向，并構(gòu)造對(duì)應(yīng)的SURF特征描述子，S2中，確定特征點(diǎn)主方向并構(gòu)造SURF特征描述子，提取不同模態(tài)醫(yī)學(xué)圖像的特征點(diǎn)具體為：
S2.1、根據(jù)SURF算法特性統(tǒng)計(jì)特征點(diǎn)鄰域內(nèi)的Harr小波特征，選擇最長(zhǎng)矢量的方向?yàn)樵撎卣鼽c(diǎn)的主方向，具體步驟為：
S2.1-1、統(tǒng)計(jì)任一特征點(diǎn)領(lǐng)域內(nèi)的Harr小波特征，即以特征點(diǎn)為中心，計(jì)算半徑為6s的鄰域內(nèi)，s為特征點(diǎn)所在的尺度值，統(tǒng)計(jì)60°扇形內(nèi)所有點(diǎn)在水平x和垂直y方向的Harr小波響應(yīng)總和(Harr小波邊長(zhǎng)取4s，s為特征點(diǎn)所在的尺度值，并給這些Harr小波響應(yīng)值賦上高斯權(quán)重系數(shù)，高斯權(quán)重系數(shù)由高斯模型在不同的角度和距離上確定，權(quán)重系數(shù)的和為1)；然后將60°范圍內(nèi)的響應(yīng)相加以形成特征矢量(特征值加上方向信息即形成特征矢量)，遍歷整個(gè)圓形區(qū)域，覆蓋整個(gè)360°，選擇最長(zhǎng)矢量的方向?yàn)樵撎卣鼽c(diǎn)的主方向。
遍歷所有的特征點(diǎn)，得到每個(gè)特征點(diǎn)的主方向。
S2.2、根據(jù)S2.1步驟中得到的方向構(gòu)造對(duì)應(yīng)特征點(diǎn)的SURF特征描述子，具體方法為：
S2.2-1、在特征點(diǎn)鄰域范圍內(nèi)取一個(gè)正方形框，框的邊長(zhǎng)為20s(s是該特征點(diǎn)所在的尺度值)，該框的方向即為步驟S2.1得到的主方向；
S2.2-2、把該框分為16個(gè)子區(qū)域，每個(gè)子區(qū)域統(tǒng)計(jì)25個(gè)像素的水平方向和垂直方向的Harr小波特征，前述水平和垂直方向均為相對(duì)特征點(diǎn)的主方向而言。該Harr小波特征為：水平方向灰度值之和∑dx，水平方向灰度值絕對(duì)值之和∑|dx|，垂直方向灰度值之和∑dy，以及垂直方向灰度值絕對(duì)值之和∑|dy|，對(duì)于每個(gè)特征點(diǎn)，建立64維向量作為該特征點(diǎn)的SURF特征描述子，其中列向量對(duì)應(yīng)16個(gè)子區(qū)域，行向量對(duì)應(yīng)各子區(qū)域的4個(gè)Harr小波特征參數(shù)值。
S3、對(duì)特征點(diǎn)的SURF特征描述子進(jìn)行矩陣運(yùn)算，獲取轉(zhuǎn)換矩陣，步驟S3具體為：
S3.1對(duì)于各模態(tài)醫(yī)學(xué)圖像中相應(yīng)位置的圖像以遍歷的方式計(jì)算兩張圖像所有特征點(diǎn)描述子的內(nèi)積，
即對(duì)于來自兩張圖像的對(duì)應(yīng)的兩個(gè)特征點(diǎn)，計(jì)算64維特征向量的內(nèi)積，按照數(shù)值從大到小進(jìn)行排序，得到排序后的特征點(diǎn)序列，其中數(shù)值最大者為最匹配的點(diǎn)；
S3.2、對(duì)前述排序后的特征點(diǎn)序列進(jìn)行矩陣運(yùn)算，得到配準(zhǔn)需要的3X3轉(zhuǎn)換矩陣。
步驟S3還包括以下步驟：設(shè)定閾值，在排序后的特征點(diǎn)序列中，選擇大于前述閾值的配對(duì)特征點(diǎn)進(jìn)行矩陣運(yùn)算；閾值的設(shè)置方法為：通過對(duì)大量不同模態(tài)的醫(yī)學(xué)圖像的機(jī)器學(xué)習(xí)，設(shè)置特征點(diǎn)的評(píng)判閾值。
S4、使用轉(zhuǎn)換矩陣完成不同模態(tài)醫(yī)學(xué)圖像配準(zhǔn)，配準(zhǔn)過程如下：
S4.1使用S3.3計(jì)算得到的轉(zhuǎn)換矩陣對(duì)待配準(zhǔn)圖像進(jìn)行卷積，得到每個(gè)像素點(diǎn)新的坐標(biāo)，并使用雙三次插值方法計(jì)算配準(zhǔn)后的像素值。從而完成兩幅不同模態(tài)圖像的配準(zhǔn)。
本發(fā)明的SURF算法相比于SIFT算法，計(jì)算速率大大提高了。利用GPU等并行化技術(shù)可以使得大分辨率圖像的線性配準(zhǔn)達(dá)到實(shí)時(shí)。
本發(fā)明未涉及部分均與現(xiàn)有技術(shù)相同或可采用現(xiàn)有技術(shù)加以實(shí)現(xiàn)。