本發(fā)明涉及醫(yī)學圖像、計算機視覺、圖像處理等領域，具體為一種基于icp算法和b樣條的錐形束ct和ct圖像配準方法。
背景技術：
：錐形束ct具有成像快、非侵入性、低輻射等優(yōu)點，并且它能夠?qū)崟r的顯示病人軟組織和骨骼信息，因此錐形束ct成像技術可用于實時手術導航中。然而在骨科手術導航中，由于ct成像對骨骼的清晰度較高，所以術中錐形束ct圖像常常要與術前ct圖像進行圖像配準以進行精準治療。但是骨科患者的病變部位可能產(chǎn)生變形，從而導致錐形束ct和ct進行配準時的精度成為難點，并且配準的速度也是手術導航中彈性配準的技術重點。目前錐形束ct和ct圖像的配準大致分為兩種：基于特征點的配準和基于體素的配準。這兩種配準方法的優(yōu)勢有所不同，基于特征的配準方法只使用局部特征作為配準要素，所以其配準的速度較快，但由于舍棄一部分圖像信息，因此其精確度有所欠缺。基于體素的配準方法由于使用了圖像的全部像素點信息，所以其配準準確度較高。但是計算量龐大導致基于體素的配準的速度較慢。近年來，將兩種方法結合的混合算法相繼出現(xiàn)，但仍未應用到錐形束ct和ct圖像配準中。技術實現(xiàn)要素：針對以上問題，本發(fā)明提供了一種基于icp算法和b樣條的錐形束ct和ct圖像配準方法，這種方法結合基于點集的配準(icp算法)和基于體素的配準(b樣條配準)，將基于點集的配準結果作為基于體素配準的初值可以快速進行配準，同時基于體素的配準也使本發(fā)明的方法具有很高的精度，因此本發(fā)明可以快速準確的配準圖像，滿足其在實時骨科手術導航中的應用。為實現(xiàn)上述目的，本發(fā)明提供如下技術方案：一種基于icp和b樣條的錐形束ct和ct圖像配準方法，步驟如下：為實現(xiàn)上述目的，本發(fā)明提供如下技術方案：一種基于icp和b樣條的錐形束ct和ct圖像配準方法，步驟如下：步驟1：獲取待配準圖像(ct圖像)i1和參考圖像(錐形束ct圖像)i2；步驟2：對所述待配準圖像和參考圖像進行圖像分割，將所述待配準圖像中的目標物分割出來并提取其坐標點，然后對坐標點進行采樣生成點集數(shù)據(jù)p＝{p1,p2,…,pt}，其中p1,p2,…,pt表示目標物的坐標點，t表示點集p中坐標點個數(shù)，將所述參考圖像中的目標物分割出來并提取其坐標點，然后對坐標點進行采樣生成點集數(shù)據(jù)q＝{q1,q2,…,qt},其中q1,q2,…,qt表示目標物的坐標點，點集q中坐標點個數(shù)同點集p中的個數(shù)一樣；步驟3：對所述待配準圖像點集數(shù)據(jù)p和參考圖像點集數(shù)據(jù)q進行基于點集的配準即icp仿射配準，其步驟主要為：首先對所述仿射變換陣maffine進行初始化；利用仿射變換陣maffine對所述待配準點集p進行變換得到點集p，即p，＝p*maffine；然后使用公式計算所述參考圖像點集q和p，之間的均方誤差，其中pi和qi分別表示第i個坐標點；使用lbfgs算法不斷優(yōu)化maffine對均方誤差求最小化，得到最優(yōu)解maffine；步驟4：根據(jù)所述仿射變換陣maffine對所述待配準圖像進行仿射配準得到icp配準圖像iicp；步驟5：將獲取的icp配準圖像iicp作為基于體素的配準初值，即對所述仿射配準結果iicp和所述參考圖像i2進行b樣條配準得到最終結果ifinal。優(yōu)選的，所述步驟2中的對所述待配準圖像和參考圖像進行圖像分割的步驟，包括：對所述待配準圖像和參考圖像的像素灰度進行歸一化為[0,255]；對所述歸一化處理后的待配準圖像和參考圖像進行閾值分割；對所述閾值分割后的待配準圖像和參考圖像進行像素坐標點提取和采樣并生成點集數(shù)據(jù)。優(yōu)選的，對所述待配準圖像和參考圖像的像素灰度進行歸一化時采用的公式為：其中g和g，分別表示歸一化后灰度值和原始圖像ct值，a和b分別用公式a＝wc-ww/2，b＝wc+ww/2，其中ww表示圖像上ct值的窗口范圍，wc表示圖像上ct值窗口中心。優(yōu)選的，對所述歸一化處理后的待配準圖像和參考圖像進行閾值分割的目標物為骨骼區(qū)域，其中骨骼的ct值范圍為[100,1000],使用歸一化公式計算其灰度值。優(yōu)選的，所述閾值分割后的待配準圖像和參考圖像進行像素坐標點提取和采樣并生成點集數(shù)據(jù)時首先提取出256灰度級中骨骼范圍的坐標點，然后對提取出的坐標點進行1/50等間距采樣得到點集圖像。優(yōu)選的，所述步驟5中，b樣條配準采用多分辨率b樣條彈性配準方法，步驟為：使用公式對輸入圖像進行高斯平滑處理，根據(jù)圖像大小進行降采樣實施多分辨率配準，先將圖像在低分辨進行b樣條粗配準。首先初始化控制點，然后對兩幅圖像的相似度測量，通過對相似度值的判斷修改控制點，將修改后的控制點代入b樣條變換函數(shù)對待配準圖像進行變換。使用優(yōu)化方法不斷修改控制點對相似度值進行最小化求解，從而得到最終配準結果ifinal。優(yōu)選的，所述多分辨率配準的分辨率層數(shù)分為1到3層，每層中xyz軸最小分辨率參數(shù)分別設為128×128×10,256×256×20,512×512×40。優(yōu)選的，所述步驟5中的b樣條配準所用變換模型為b樣條函數(shù)，用公式求取像素點在(x，y，z)位置處的形變場，其中其中表示向下取整，φi,j,k表示x，y，x軸控制點間距分別為δi,δj,δk，大小為nx×ny×nz的網(wǎng)格中序號為ijk的控制點的值，其中i，j，k為了求像素點在(x,y,z)處的位移值t(x,y,z)僅與相鄰的43個控制點有關；bl,bm,bn均為三次b樣條函數(shù)基函數(shù)，bl表示b樣條函數(shù)中第l個基函數(shù)，其計算公式為b0(u)＝(1-u)3/6，b1(u)＝(3u3-6u2+4)/6，b2(u)＝(-3u3+3u2+3u+1)/6，b3(u)＝u3/6,其中0≤u≤1bm和bn的計算同bl，這些函數(shù)起到權函數(shù)的作用，他們根據(jù)控制點到(x,y,z)的距離大小來加權每個控制點對t(x,y,z)的貢獻。優(yōu)選的，所述多分辨率配準的每層分辨率的b樣條網(wǎng)格間距分別設為：32×32×16,16×16×8,8×8×4。優(yōu)選的，所述相似度測量和優(yōu)化方法采用灰度均方差作為相似度測量，優(yōu)化方法選用lbfgs算法。本發(fā)明與現(xiàn)有技術相比的優(yōu)點在于：本發(fā)明基于icp和b樣條的錐形束ct和ct圖像配準方法，采用結合基于點的配準和基于體素的配準，在對錐形束ct和ct圖像配準時，相對于傳統(tǒng)方法具有速度快的同時保持配準精度高的特點，可以很好的應用于基于錐形束ct成像技術的骨科手術導航系統(tǒng)中。附圖說明圖1為基于icp和b樣條的配準方法流程圖；圖2為對圖像分割和點云提取采樣后的點集；其中a是從錐形束ct圖像中獲取的點集，b是從ct圖像中獲取的點集；圖3為對點集圖像進行icp配準的效果展示；其中a為icp配準前的點集圖像，b為icp配準后的點集圖像；圖4為本發(fā)明配準效果展示圖；其中a為初始錐形束ct圖像，b為初始ct圖像，c為icp配準結果，d為最終配準結果；圖5為仿真圖像生成原理展示圖；其中a為生成標準形變場過程示意圖，b為生成目標形變場過程示意圖；圖6為配準結果的差值對比；其中a為配準前初始錐形束ct和ct圖像的差值圖，b為icp配準結果與初始錐形束ct圖像的差值圖，c為最終配準結果與初始錐形束ct圖像的差值圖；圖7為配準誤差的直方圖顯示。具體實施方式由于錐形束ct和ct圖像取自不同成像裝置并且拍攝時間不同，因此術前規(guī)劃系統(tǒng)中的待配準圖像即ct圖像和術中圖像即錐形束ct圖像并不相同，例如兩者成像視角不同，每層切片的對應關系可能有出入，圖像會有偏移和形變，需要進行配準矯正。因此本發(fā)明采用如下的基于icp和b樣條的配準方法。如圖1所示，本發(fā)明一種基于icp和b樣條的配準方法流程圖，包括以下步驟：步驟s1,獲取待配準圖像(ct圖像)i1和參考圖像(錐形束ct圖像)i2。步驟s2：對所述待配準圖像和參考圖像進行圖像分割，將所述待配準圖像中的目標物分割出來并提取其坐標點，然后對坐標點進行采樣生成點集數(shù)據(jù)p＝{p1,p2,…,pt}，其中p1,p2,…,pt表示目標物的坐標點，t表示點集p中坐標點個數(shù)，將所述參考圖像中的目標物分割出來并提取其坐標點，然后對坐標點進行采樣生成點集數(shù)據(jù)q＝{q1,q2,…,qt},其中q1,q2,…,qt表示目標物的坐標點，點集q中坐標點個數(shù)同點集p中的個數(shù)一樣；在本實例中，步驟s2中的對待配準圖像和參考圖像進行圖像分割的步驟，包括：(1)對所述待配準圖像和參考圖像的像素灰度進行歸一化為[0,255]；(2)對所述歸一化處理后的待配準圖像和參考圖像進行閾值分割；(3)對所述閾值分割后的待配準圖像和參考圖像進行像素坐標點提取和采樣并生成點集數(shù)據(jù)。由于圖像是取自不同時間，不同設備或者不同視角的同一場景的圖像，因此ct圖像和錐形束ct圖像的分辨率會有所不同，所以將圖像的像素灰度值進行歸一化處理。歸一化處理所使用的公式為：其中g和g，分別表示歸一化后灰度值和原始圖像ct值，a和b分別用公式a＝wc-ww/2，b＝wc+ww/2，其中ww表示圖像上ct值的窗口范圍，wc表示圖像上ct值窗口中心。在本實例中，閾值分割的目標可以是：肺部、骨骼。其中肺部的ct值范圍為[-950,200],骨骼的ct值范圍為[100,1000]。優(yōu)選地，對骨骼進行分割。在本實例中，提出256灰度級中相應器官范圍的坐標點，然后對坐標點進行采樣，采樣方法的采用如下方式之一：隨機采樣、等間隔采樣。優(yōu)選地，采用1/50等間隔采樣。如圖2所示，為對待配準圖像和參考圖像進行閾值分割后提取坐標點并進行采樣后的實驗展示圖，其中a是從錐形束ct圖像中獲取的點集，b是從ct圖像中獲取的點集。s3：對所述待配準圖像點集數(shù)據(jù)p和參考圖像點集數(shù)據(jù)q進行基于點集的配準即icp仿射配準，其步驟主要為：首先對所述仿射變換陣maffine進行初始化；利用仿射變換陣maffine對所述待配準點集p進行變換得到點集p，即p，＝p*maffine；然后使用公式計算所述參考圖像點集q和p，之間的均方誤差，其中pi和qi分別表示第i個坐標點；使用lbfgs算法不斷優(yōu)化maffine對均方誤差求最小化，得到最優(yōu)解maffine。在本實例中，icp配準是基于點的配準，其中所使用的參數(shù)優(yōu)化方法為lbfgs算法，其搜索步長為0.1mm、梯度收斂公差為0.01，最大迭代次數(shù)為100次在本發(fā)明實例中，icp配準的變換模型采用如下方式之一：剛性變換(旋轉(zhuǎn)、平移)，仿射變換(旋轉(zhuǎn)、平移、縮放、剪切)。優(yōu)選地，采用仿射變換。如圖3所示，為對待配準圖像點集和參考圖像點集進行icp配準前(圖3中的a)和配準后(圖3中的b)的實驗對比圖。步驟s4：根據(jù)仿射變換陣maffine對待配準圖像進行仿射配準得到icp配準圖像iicp；在本實例中，對待配準圖像進行仿射配準時的插值方法可采用如下方式之一：線性插值、三次插值。優(yōu)選地，選用線性插值。步驟s5：將獲取的icp配準圖像iicp作為基于體素的配準初值，即對仿射配準結果iicp和參考圖像i2進行多分辨率b樣條彈性配準得到最終結果ifinal。將icp配準圖像和參考圖像作為b樣條彈性配準的初始圖像以提高配準精度和配準速度。同時b樣條配準前，對icp配準圖像和參考圖像進行高斯濾波處理用以提高配準速度。在本實例中，b樣條配準的主要步驟為：對輸入圖像進行平滑處理，根據(jù)圖像大小進行降采樣實施多分辨率配準，先將圖像在低分辨進行b樣條粗配準。首先初始化控制點，然后對兩幅圖像的相似度測量，通過對相似度值的判斷修改控制點，將修改后的控制點代入b樣條變換函數(shù)對待配準圖像進行變換。使用優(yōu)化方法不斷修改控制點對相似度值進行最小化求解，從而得到最終配準結果ifinal。在本實例中，b樣條彈性配準部分采用多分辨率配準，分辨率根據(jù)具體圖像大小可分為1到3層，其中每層分辨率的最小分辨率參數(shù)分別為128×128×10,256×256×20,512×512×40。例如當圖像大小等于或超過512×512×40(如512×512×40)時可以分為三層進行配準，相鄰層數(shù)之間為倍數(shù)關系，(即第一層128×128×10，第二層256×256×20，第三層512×512×40)。首先進行低分辨率配準，然后得出的結果作為高分辨率的輸入再進行高分辨率配準直至結束。在本實例中，每層分辨率的b樣條網(wǎng)格間距分別選為：32×32×16,16×16×8,8×8×4。在本實例中，相似度測量方法可選以下方式之一：互信息、互相關、灰度均方差。優(yōu)選地，選用灰度均方差作為相似度測量。參數(shù)優(yōu)化選用lbfgs算法。如圖4所示為整個配準實驗其中一個切片的結果展示，其中圖4中的(a)為參考圖像(錐形束ct)，圖4中的(b)為待配準圖像(ct圖像)，圖4中的(c)為icp配準結果，圖4中的(d)為最終配準結果?？梢钥闯鲈诮?jīng)過配準后，待配準圖像與參考圖像的接近程度逐漸增加。如表1示，為ct(512*512*53)和錐形束ct(410*410*42)圖像配準的結果，其中灰度均方差值作為相似度測量。配準輸出參數(shù)包括配準前參考圖像和待配準圖像的灰度均方差值，icp配準結果和b樣條配結果和參考圖像的灰度均方差值對比以及配準所需時間。從配準前后的灰度均方差值對比可以看出灰度均方差值提升較多，說明配準后待配準圖像和參考圖像的一致性有很大提高，此配準方法對于ct和錐形束ct圖像配準效果很好。表1配準前icp配準結果b樣條配準結果配準時間灰度均方差值25951208575630為對本發(fā)明的準確度進行測量，設計仿真實驗。如圖5為仿真圖像生成的流程展示，圖5中的(a)將錐形束ct和ct進行配準時產(chǎn)生的形變場作為標準值，圖5中的(b)中用標準形變場的逆變換對錐形束ct進行形變生成仿真ct圖像，然后對仿真ct圖像和錐形束ct圖像進行配準生成目標形變場作為目標值。圖6為對仿真圖像進行icp配準和b樣條配準后的圖像差值顯示，圖6中的(a)為配準前的差值效果展示，圖6中的(b)為icp配準結果和錐形束ct之間像素差值效果展示，圖6中的(c)為最終配準結果和錐形束ct之間像素差值效果展示。從圖中可以看出，icp配準后參考圖像和待配準圖像的空間位置大致重合，但仍有局部差異，在進行b樣條配準之后其配準結果和參考圖像基本重合。將標準值和目標值相減求出形變誤差，圖7畫出誤差的直方圖，從圖中可以看出誤差大多分布在較低值，通過求均值算出誤差值為0.51mm。準確度較高，滿足手術導航的精度要求。表2為比較本實例中配準方法與傳統(tǒng)基于體素配準方法的配準速度而設計的實驗結果。該實驗中傳統(tǒng)的配準方法也是使用放射變化和b樣條配準作為配準流程，同時相似度測量選用灰度均方差，優(yōu)化方法選用lbfgs算法。在此基礎上，兩種方法的收斂條件設為相同。傳統(tǒng)方法的配準時間和本發(fā)明的配準方法的配準時間如下表所示：表2傳統(tǒng)方法本發(fā)明配準方法時間(s)ct/錐形束ct1610630以上所述僅為本發(fā)明的較佳實施例而已，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內(nèi)所作的任何修改、等同替換和改進等，均應包含在本發(fā)明的保護范圍之內(nèi)。當前第1頁12