腦腫瘤手術規劃及術中，術前磁共振（MR）圖像與術中超聲（US）圖像的配準甚為關鍵。考慮到兩種模態圖像具有不同密度范圍及分辨率，且超聲圖像存在較多的斑點噪聲干擾，采用一種基于局部鄰域信息的自相似性上下文（SSC）描述子定義圖像之間的相似性測度。將超聲圖像作為參考，使用三維微分運算提取其中角點作為關鍵點，并采用密集位移采樣離散優化算法實施配準。整個配準過程分為仿射配準和彈性配準兩個階段，在仿射配準階段，對圖像進行多分辨率分解，在彈性配準階段，采取最小卷積和均值場推理策略對關鍵點的位移向量進行正則化處理。對22名患者的術前MR和術中US圖像進行配準實驗，仿射配準后的誤差為（1.57 ± 0.30）mm，每對圖像配準平均耗時1.36 s；彈性配準后的誤差為（1.40 ± 0.28）mm，平均用時1.53 s。實驗結果證明本文采用的方法具有良好的配準精度和速度。
引用本文： 譚振霖, 郭圣文. 基于關鍵點的超聲圖像與磁共振圖像多分辨率離散優化配準方法. 生物醫學工程學雜志, 2023, 40(2): 202-207, 216. doi: 10.7507/1001-5515.202211022 復制
0 引言
以術前MR圖像為基礎，使用計算機外科神經導航是規劃和切除腦部腫瘤的重要方法[1]。雖然神經導航為外科醫生實施手術提供了一種高效的工具，但在手術過程中，腦組織往往發生移位和變形，使得基于術前MR圖像的神經導航系統失去準確性。術中磁共振成像（intraoperative magnetic resonance imaging，iMRI）可以在一定程度上解決此問題，但iMRI成本高且非常耗時，從而制約了其臨床應用。而術中超聲（intraoperative ultrasound，iUS）因具有實時性強、操作方便、成本低、對手術工作流程干擾少等優勢，成為iMRI的有力替代方案，可用于術中實時監測組織移位，更新神經導航系統中的組織位置變化信息[2]。因此，將iUS與術前MR進行配準可以對神經導航系統因為腦移位發生的誤差進行修正，幫助醫生進行更精準的手術。
根據臨床的需求，評價配準任務好壞最重要的兩個指標是配準時間和配準精度。根據配準算法使用的特征類型可分為基于幾何特征的圖像配準和基于灰度的圖像配準。其中基于幾何特征的配準算法，首先要檢測兩幅圖像中的關鍵幾何特征，如腫瘤邊界、腦溝或血管等。Reinertsen等[3]提出基于血管的配準方法進行腦移位校正，Coupé等[4]使用了一種基于大腦高回聲結構匹配的新概率函數，Farnia等[5]提出了一種基于MR圖像中的腦溝和腫瘤邊界等回聲結構與US圖像進行匹配的新型混合方法。該類方法直接利用了空間位置信息，配準過程簡單直接，但其難點在于特征的定位較為困難，并且容易出錯。如一幅圖像中存在的幾何特征，并不總是出現在另一幅圖像中，尤其是多模態圖像的情形，難以可靠地檢測到足夠的對應特征。因此，這種配準方法一般適宜于特定數據。
基于灰度的配準算法，利用兩幅圖像的體素灰度信息來尋求不同圖像之間的對應關系，其關鍵是體素灰度信息的相似性度量和尋找對應關系的優化過程。Wein等[6]將基于線性組合的線性相關性（linear correlation of linear combination，
）的相似性度量公式應用于MR和US的配準當中，Heinrich等[7]將量化自相似上下文（self-similarity context，SSC）描述符用于定義相似性度量，Rivaz等[8]提出一種基于塊的相關比率（correlation ratio，CR）的度量配準US和MR，Machado等[9]基于由量化匹配可靠性獲取的最佳Gabor屬性來定義相似性度量，采用DRAMMS[10]可變形配準算法對US和MR圖像進行配準。基于灰度的配準方法對不同模態圖像配準的穩定性較好，但該類方法的難點在于找到穩健的相似性度量函數和高效優化算法。這些方法雖然很好地解決了US與MR圖像的多模態相似性度量問題，但對大尺寸、高分辨率圖像的配準非常耗時，時間高達幾十秒甚至幾分鐘。
為了解決配準的時間問題，本文提出了一種基于關鍵點的多分辨率離散優化配準方法。首先，以超聲圖像作為目標，通過計算灰度梯度值獲取角點，以此作為關鍵點；然后，依據超聲圖像中的關鍵點，采用SSC與密集位移采樣離散優化算法，實施仿射配準。此外，為了進一步提高計算效率，在仿射配準階段，對圖像進行多分辨率分解，以逐步完成由粗到精的配準過程；最后，采用樣條變換進行彈性配準。為了防止發生非真實的彈性形變，采用最小卷積和均值場推理操作，對關鍵點的位移向量進行正則化。
1 配準方法
1.1 配準流程
本文將三維US圖像作為參考目標，以F表示，三維MR圖像作為移動對象，用M表示，配準的關鍵是在三維空間尋找二者之間相關點區域的匹配關系。本文采用的配準方法包括四步：首先，從F中提取角點，將其作為關鍵點 
，K的坐標用 
表示。然后進行仿射配準，此階段主要是消除剛性運動引起的組織位置差異，故只需采用部分的關鍵點進行匹配，F中關鍵點 
為中心，選擇包括其鄰域體素的小立方體（塊）VF，大小為m × n × d，在M中選擇具有相同坐標位置的中心點，選擇具有同樣大小的塊VM，經位移向量 
變化后的對應塊為VMM，計算VF與VMM之間的相似性度量，即SSC，通過密集位移采樣離散優化算法來尋找關鍵點的對應位移向量。為了提高計算效率，此階段對圖像進行多分辨率分解，采用由粗到精的多階段配準策略。最后，采用樣條變換進行彈性配準，此階段采用全部關鍵點，以實施更加精確的空間配準。為了提高配準精度，采用最小卷積和均值場推理方法，對關鍵點的位移向量進行正則化。本文方法命名為基于關鍵點的多分辨率離散優化配準方法（multi-resolution discrete optimization registration based on key points，MRDOR-KP），配準流程如圖1所示。
圖1
配準流程
Figure1.
Process of registration
在仿射配準階段，由于只考慮均勻的全局剛性運動，所以只使用部分關鍵點，得到關鍵點的變換參數后，使用最小二乘法[11]計算仿射變換矩陣；在彈性配準階段，則需考慮局部的非均勻運動，使用了全部的關鍵點，采用薄板樣條（thin plate spline，TPS）[12]函數計算彈性形變場。
1.2 關鍵點檢測
均勻劃分關鍵點的方法存在有效信息利用不充分及計算代價高等不足，故本研究采用三維角點檢測算法[13]，該算法的計算公式如下：
 |
其中，
為高斯核，
為圖像F在x、y、z軸三個方向的灰度梯度 
、
，將灰度梯度向量及其轉置相乘后得到梯度的協方差矩陣，將梯度的協方差矩陣求逆后，與 
進行卷積運算，trace為矩陣的跡計算。角點特征的定義為圖像灰度梯度 
、
的值同時盡可能地大，而梯度協方差矩陣的逆的跡越小表示 
、
、
的值越大，以此來對圖像的角點進行檢測。其中高斯核卷積的作用為對圖像梯度進行平滑處理，可以使圖像的邊緣模糊，抑制US圖像中斑點噪聲的響應值，使得獲取的關鍵點更有意義，增加配準的穩定性。
值代表每個體素的特征響應值，值越高表示該點越符合角點的特征，而角點往往是重要的解剖結構。為了確保角點的有效性和提高計算效率，對之進行了非極大值抑制，即在每個邊長為 
的立方體區域，選擇響應值最大者作為關鍵點。
1.3 相似性測度
為了在F中搜索到和M中關鍵點在空間上最匹配的點，采用基于三維立方體區域（子塊）的自相似性上下文SSC，定義兩個子塊的相似性測度[7]，其中SSC計算公式如下：
 |
x為圖像對應坐標點，y定義為與x相鄰的6個點，對這6個點中相鄰的點對計算像素灰度差的平方和（sum of squared differences，SSD），得到12個值，以反映對應點x領域的自相似性信息，
為該12個SSD的標準差。根據SSC，兩個子塊的相似性測度計算公式如下：
 |
其中 
和 
分別對應F和M中子塊的SSC值，k與l分別為F和M的子塊中心，p為位移，
，P表示局部立方體搜索區域。式(2)中的相似性測度，實際上是搜索區域內兩個子塊SSC值之差的絕對值均值，它不僅適合量度不同模態圖像塊之間的相似性，而且受噪聲影響小。
1.4 離散優化算法
根據F中關鍵點的空間位置(kx, ky, kz)，以M中相同位置點為中心，在其鄰域內搜索與之相匹配的點，用 
表示該鄰域，位移（步長）分別用dx、dy、dz表示，則有：
 |
將關鍵點的搜索位移空間定義為密集采樣的三維集合d∈L=
，其中，lmax為采樣點數，q為搜索步長。通過對
的值進行設置，就可以覆蓋兩幅圖像之間潛在的所有位移運動。當進行大范圍的搜索時，采用較多的采樣點數和較大的搜索步長，并隨著搜索范圍的一步步減少，采樣點數和搜索步長也相應減少，該方法稱為密集位移采樣離散優化算法。
本研究將關鍵點采樣和塊相似性離散優化策略相結合，如在大范圍搜索時，采用少量的關鍵點及較大圖像塊P，從而獲取大的感受野范圍來保證配準精度，而較少的關鍵點可提高計算速度；反之，在小范圍搜索時使用較多的關鍵點及較小圖像塊P，可獲取較小的感受野范圍，從而更精準地捕捉局部微小形變，此時使用較多的關鍵點，可提升配準精度。在進行大范圍搜索時，圖像塊感受野較大，配準的計算復雜度也隨之加大。為了解決該問題，采用多分辨率配準的方法，在大范圍搜索時使用分辨率較低的圖像，得到粗略的位置偏移；隨著范圍搜索的減少，圖像分辨率提高，以感知精細的局部形變，確保配準精度。
1.5 位移向量的正則化
在圖像經局部相似性搜索和離散優化后，所有關鍵點都獲得了密集位移采樣的相似性損失矩陣。為了進一步消除配準誤差，采用最小卷積和均值場推理的位移正則化方法[14]。
該方法采用了最小池化和平均池化兩個操作來對點的相似性損失進行正則化，首先使用最小池化找到局部相似性損失D的近似最小值，然后使用兩個平均池化操作進行二次平滑來進行正則化。為了實現上述池化操作，本文將整個相似性損失矩陣排列成了一個6維的矩陣，其中1～3維按照關鍵點
采樣位置排列，4～6維按照密集位移采樣d對應位置進行排列。首先對4～6維進行最小池化（最小卷積）和平均池化操作，來找到搜索空間上的近似局部最小值和平滑，接著對1～3維進行平均池化（均值場推理），來對相鄰關鍵點的位移進行約束。由于池化操作可以在GPU上進行并行運算，故可大大提升計算速度。
2 實驗數據與配準結果
2.1 實驗數據
數據來源于公共數據集RESECT[15]，該數據集使用1.5T Siemens Magnetom Avanto對22名患者進行T1w和T2-FLAIR兩種序列MRI掃描，體素大小為1 mm × 1 mm × 1 mm，分辨率為256 × 256 × 192。使用12FLA-L線性超聲探頭分別對患者腫瘤切除前、切除中和切除后進行掃描，體素大小為0.14 mm × 0.14 mm × 0.14 mm～0.24 mm × 0.24 mm × 0.24 mm。本研究選擇手術前T2-FLAIR MR圖像和術中US圖像進行配準實驗，在每對MR圖像和US圖像之間，均事先由人工確定了15～16個對應的標記點，以供各種配準算法進行性能評估和比較。
2.2 關鍵點檢測
角點是圖像亮度變化劇烈的點或圖像邊緣曲線上曲率極大值的點。這些點在保留圖像圖形重要特征的同時，可以有效地減少信息的數據量，使其信息的含量很高，有效地提高了計算速度，有利于圖像的可靠匹配。圖2顯示了US圖像的角點檢測的響應圖、用局部最高響應值篩選后的關鍵點，以及配準后MR圖像中對應的關鍵點。
圖2
關鍵點可視化
Figure2.
Visualization of key points
從圖中可以看出，在US圖像上檢測到的關鍵點，大部分位于重要解剖結構上，如腦脊液、灰質和腫瘤邊緣等。
2.3 實驗結果與分析
2.3.1 配準實驗環境及參數設置
實驗使用Dell precision tower 7910工作站，CPU為Intel Xeon E5-2690 v4，配置NVIDIA GeForce GTX 1080 Ti顯卡的GPU。整個配準流程的計算工作在GPU上實現。采用對應標記點之間的平均歐氏距離，作為平均配準誤差（mean target registration error of landmarks，mTRE），對配準結果進行評估。
配準參數選取如下：在仿射配準進行了3個階段的從粗到細的配準，分別將圖像的分辨率大小采樣為原圖的0.5、0.75和1，關鍵點的非極大值抑制的參數 
。該階段只使用了角點特征響應值最高的部分關鍵點，關鍵點數量的參數為64、128、216，密集搜索空間的搜索點數 
，搜索步長 
。在彈性配準時，圖像已經進行了仿射的預配準，圖像之間只有微小的誤差，因此只進行了單個階段的精細配準，關鍵點非極大化抑制的參數 
。該階段使用了所有的關鍵點，并對關鍵點進行了閾值化處理，閾值化處理采用的是Otsu算法[16]，密集搜索空間的搜索點數 
，搜索步長 
。
2.3.2 配準結果和對比分析
為證明本文所提出配準方法MRDOR-KP的配準效果，我們在RESECT數據集的22個病例中進行了驗證。該數據集的初始配準誤差為（5.41 ± 4.19）mm，使用MRDOR-KP方法經過仿射配準后數據的配準誤差為（1.57 ± 0.30）mm，平均配準時間為1.36 s，彈性配準后的配準誤差為（1.40 ± 0.28）mm，平均配準時間為1.53 s。對于仿射配準，MRDOR-KP方法能夠有效地對不同大小形變的圖像進行配準，即使圖像存在比較大的形變，如病例12有19.74 mm的初始誤差，該方法依舊有效。通過彈性形變進行擬合后，絕大部分患者腦圖像的配準誤差在2 mm以內，最大配準誤差為2.02 mm，表明對于所有案例的配準結果都比較好。
為了驗證多分辨率策略和正則化方法的有效性，表1列出了使用/不使用多分辨率分解和正則化策略的配準誤差及運算時間：
表1
不同模塊性能比較
Table1.
Performance comparison of different modules
從表1可以看出，仿射配準階段使用多分辨率策略來進行配準，不僅可以顯著提高配準的速度，并且能夠提升配準精度；彈性配準階段，正則化操作能有效降配準誤差，因采用了并行運算，故正則化計算時間僅增加了0.19 s。
為了評估本文方法的性能，我們將本文的配準結果與CuRIOUS2018挑戰賽[17]排名前3的算法ImFusion[18]、DeedsSSC[19]和cDRAMMS[20]進行了比較，結果如表2所示。
表2
不同方法性能比較
Table2.
Performance comparison of different methods
從表2的配準結果比較可以看出，本文方法成功地將配準時間降低到了2.89 s，比其他方法大幅縮短，并且取得了最低的配準誤差。結果表明，我們的方法可以自動精準地對US與MR圖像進行配準，只需要幾秒就能實現，能夠滿足醫生對于臨床手術精度和時間的需求。
2.3.3 結果可視化
圖3展示了RESECT數據集中三名不同患者（case#12、case#17和case#23）的MR到US的配準結果，為了進一步評估配準結果的好壞，具有臨床經驗的神經外科醫生對配準前后的MR圖像和US圖像對關鍵區域（腦溝和腫瘤邊界）進行了標記。其中白色和藍色箭頭分別指出US和MR的腫瘤邊界，黃色和紅色箭頭分別指出US和MR的腦溝。從配準前后箭頭的對齊效果可以看出經過我們的方法配準后，圖像關鍵區域得到了明顯的對齊。
圖3
部分配準結果
白色和黃色箭頭分別指向US圖像的腫瘤邊界和腦溝，藍色和紅色箭頭分別指向MR圖像的腫瘤邊界和腦溝
Figure3. Partial registration resultsthe white and yellow arrows point to the tumor boundary and sulcus of the US image respectively, and the blue and red arrows point to the tumor boundary and sulcus of the MR image respectively
圖4展示了患者#25在本文算法的配準結果與CuRIOUS2018挑戰賽給出的配準結果的對比。該病例的初始配準誤差為10.06 mm，ImFusion、cDRAMMS和DeedsSSC算法的配準誤差分別為2.12、2.76、1.45 mm，而本文算法的配準誤差僅為1.37 mm。圖中使用箭頭指向了腦溝區域的對齊情況，可以看出我們的算法取得了更好的對齊效果和更小的配準誤差。
圖4
多算法配準結果
白色箭頭指向US圖像的腦溝，右下角為不同算法的mTRE
Figure4. Registration results of multiple algorithmsthe white arrow points to the sulcus of the US image, and the lower right corner is the mTRE of different algorithms
3 結論
本文提出了一種基于關鍵點的多分辨率離散優化兩階段圖像配準算法，采用自相似性上下文描述子定義圖像間相似性測度，并對關鍵點的位移向量進行正則化處理，以減少配準誤差。實驗證明，本文提出的配準方法，具有良好的配準精度和速度，能有效實現MR圖像和US圖像的實時配準。本文的創新點包括：① 采用一種基于角點作為關鍵特征點的策略，根據運動特性不同，分別在仿射配準和彈性配準階段采用不同數量的關鍵點，以充分利用有效的解剖結構信息，并加速配準；② 在仿射配準階段，對圖像進行多分辨率分解，實施由粗到精的配準，以提高配準精度；③ 在彈性配準階段，采取最小卷積和均值場推理策略對位移向量進行正則化處理，有效地減少配準誤差。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：譚振霖負責算法研究、實驗、分析及論文撰寫，郭圣文負責總體設計、方法與論文撰寫指導和審校。
0 引言
以術前MR圖像為基礎，使用計算機外科神經導航是規劃和切除腦部腫瘤的重要方法[1]。雖然神經導航為外科醫生實施手術提供了一種高效的工具，但在手術過程中，腦組織往往發生移位和變形，使得基于術前MR圖像的神經導航系統失去準確性。術中磁共振成像（intraoperative magnetic resonance imaging，iMRI）可以在一定程度上解決此問題，但iMRI成本高且非常耗時，從而制約了其臨床應用。而術中超聲（intraoperative ultrasound，iUS）因具有實時性強、操作方便、成本低、對手術工作流程干擾少等優勢，成為iMRI的有力替代方案，可用于術中實時監測組織移位，更新神經導航系統中的組織位置變化信息[2]。因此，將iUS與術前MR進行配準可以對神經導航系統因為腦移位發生的誤差進行修正，幫助醫生進行更精準的手術。
根據臨床的需求，評價配準任務好壞最重要的兩個指標是配準時間和配準精度。根據配準算法使用的特征類型可分為基于幾何特征的圖像配準和基于灰度的圖像配準。其中基于幾何特征的配準算法，首先要檢測兩幅圖像中的關鍵幾何特征，如腫瘤邊界、腦溝或血管等。Reinertsen等[3]提出基于血管的配準方法進行腦移位校正，Coupé等[4]使用了一種基于大腦高回聲結構匹配的新概率函數，Farnia等[5]提出了一種基于MR圖像中的腦溝和腫瘤邊界等回聲結構與US圖像進行匹配的新型混合方法。該類方法直接利用了空間位置信息，配準過程簡單直接，但其難點在于特征的定位較為困難，并且容易出錯。如一幅圖像中存在的幾何特征，并不總是出現在另一幅圖像中，尤其是多模態圖像的情形，難以可靠地檢測到足夠的對應特征。因此，這種配準方法一般適宜于特定數據。
基于灰度的配準算法，利用兩幅圖像的體素灰度信息來尋求不同圖像之間的對應關系，其關鍵是體素灰度信息的相似性度量和尋找對應關系的優化過程。Wein等[6]將基于線性組合的線性相關性（linear correlation of linear combination，）的相似性度量公式應用于MR和US的配準當中，Heinrich等[7]將量化自相似上下文（self-similarity context，SSC）描述符用于定義相似性度量，Rivaz等[8]提出一種基于塊的相關比率（correlation ratio，CR）的度量配準US和MR，Machado等[9]基于由量化匹配可靠性獲取的最佳Gabor屬性來定義相似性度量，采用DRAMMS[10]可變形配準算法對US和MR圖像進行配準。基于灰度的配準方法對不同模態圖像配準的穩定性較好，但該類方法的難點在于找到穩健的相似性度量函數和高效優化算法。這些方法雖然很好地解決了US與MR圖像的多模態相似性度量問題，但對大尺寸、高分辨率圖像的配準非常耗時，時間高達幾十秒甚至幾分鐘。
為了解決配準的時間問題，本文提出了一種基于關鍵點的多分辨率離散優化配準方法。首先，以超聲圖像作為目標，通過計算灰度梯度值獲取角點，以此作為關鍵點；然后，依據超聲圖像中的關鍵點，采用SSC與密集位移采樣離散優化算法，實施仿射配準。此外，為了進一步提高計算效率，在仿射配準階段，對圖像進行多分辨率分解，以逐步完成由粗到精的配準過程；最后，采用樣條變換進行彈性配準。為了防止發生非真實的彈性形變，采用最小卷積和均值場推理操作，對關鍵點的位移向量進行正則化。
1 配準方法
1.1 配準流程
本文將三維US圖像作為參考目標，以F表示，三維MR圖像作為移動對象，用M表示，配準的關鍵是在三維空間尋找二者之間相關點區域的匹配關系。本文采用的配準方法包括四步：首先，從F中提取角點，將其作為關鍵點 ，K的坐標用 表示。然后進行仿射配準，此階段主要是消除剛性運動引起的組織位置差異，故只需采用部分的關鍵點進行匹配，F中關鍵點 為中心，選擇包括其鄰域體素的小立方體（塊）VF，大小為m × n × d，在M中選擇具有相同坐標位置的中心點，選擇具有同樣大小的塊VM，經位移向量 變化后的對應塊為VMM，計算VF與VMM之間的相似性度量，即SSC，通過密集位移采樣離散優化算法來尋找關鍵點的對應位移向量。為了提高計算效率，此階段對圖像進行多分辨率分解，采用由粗到精的多階段配準策略。最后，采用樣條變換進行彈性配準，此階段采用全部關鍵點，以實施更加精確的空間配準。為了提高配準精度，采用最小卷積和均值場推理方法，對關鍵點的位移向量進行正則化。本文方法命名為基于關鍵點的多分辨率離散優化配準方法（multi-resolution discrete optimization registration based on key points，MRDOR-KP），配準流程如圖1所示。
圖1
配準流程
Figure1.
Process of registration
在仿射配準階段，由于只考慮均勻的全局剛性運動，所以只使用部分關鍵點，得到關鍵點的變換參數后，使用最小二乘法[11]計算仿射變換矩陣；在彈性配準階段，則需考慮局部的非均勻運動，使用了全部的關鍵點，采用薄板樣條（thin plate spline，TPS）[12]函數計算彈性形變場。
1.2 關鍵點檢測
均勻劃分關鍵點的方法存在有效信息利用不充分及計算代價高等不足，故本研究采用三維角點檢測算法[13]，該算法的計算公式如下：
|
其中， 為高斯核， 為圖像F在x、y、z軸三個方向的灰度梯度 、，將灰度梯度向量及其轉置相乘后得到梯度的協方差矩陣，將梯度的協方差矩陣求逆后，與 進行卷積運算，trace為矩陣的跡計算。角點特征的定義為圖像灰度梯度 、 的值同時盡可能地大，而梯度協方差矩陣的逆的跡越小表示 、、 的值越大，以此來對圖像的角點進行檢測。其中高斯核卷積的作用為對圖像梯度進行平滑處理，可以使圖像的邊緣模糊，抑制US圖像中斑點噪聲的響應值，使得獲取的關鍵點更有意義，增加配準的穩定性。 值代表每個體素的特征響應值，值越高表示該點越符合角點的特征，而角點往往是重要的解剖結構。為了確保角點的有效性和提高計算效率，對之進行了非極大值抑制，即在每個邊長為 的立方體區域，選擇響應值最大者作為關鍵點。
1.3 相似性測度
為了在F中搜索到和M中關鍵點在空間上最匹配的點，采用基于三維立方體區域（子塊）的自相似性上下文SSC，定義兩個子塊的相似性測度[7]，其中SSC計算公式如下：
|
x為圖像對應坐標點，y定義為與x相鄰的6個點，對這6個點中相鄰的點對計算像素灰度差的平方和（sum of squared differences，SSD），得到12個值，以反映對應點x領域的自相似性信息， 為該12個SSD的標準差。根據SSC，兩個子塊的相似性測度計算公式如下：
|
其中 和 分別對應F和M中子塊的SSC值，k與l分別為F和M的子塊中心，p為位移，，P表示局部立方體搜索區域。式(2)中的相似性測度，實際上是搜索區域內兩個子塊SSC值之差的絕對值均值，它不僅適合量度不同模態圖像塊之間的相似性，而且受噪聲影響小。
1.4 離散優化算法
根據F中關鍵點的空間位置(kx, ky, kz)，以M中相同位置點為中心，在其鄰域內搜索與之相匹配的點，用 表示該鄰域，位移（步長）分別用dx、dy、dz表示，則有：
|
將關鍵點的搜索位移空間定義為密集采樣的三維集合d∈L=，其中，lmax為采樣點數，q為搜索步長。通過對的值進行設置，就可以覆蓋兩幅圖像之間潛在的所有位移運動。當進行大范圍的搜索時，采用較多的采樣點數和較大的搜索步長，并隨著搜索范圍的一步步減少，采樣點數和搜索步長也相應減少，該方法稱為密集位移采樣離散優化算法。
本研究將關鍵點采樣和塊相似性離散優化策略相結合，如在大范圍搜索時，采用少量的關鍵點及較大圖像塊P，從而獲取大的感受野范圍來保證配準精度，而較少的關鍵點可提高計算速度；反之，在小范圍搜索時使用較多的關鍵點及較小圖像塊P，可獲取較小的感受野范圍，從而更精準地捕捉局部微小形變，此時使用較多的關鍵點，可提升配準精度。在進行大范圍搜索時，圖像塊感受野較大，配準的計算復雜度也隨之加大。為了解決該問題，采用多分辨率配準的方法，在大范圍搜索時使用分辨率較低的圖像，得到粗略的位置偏移；隨著范圍搜索的減少，圖像分辨率提高，以感知精細的局部形變，確保配準精度。
1.5 位移向量的正則化
在圖像經局部相似性搜索和離散優化后，所有關鍵點都獲得了密集位移采樣的相似性損失矩陣。為了進一步消除配準誤差，采用最小卷積和均值場推理的位移正則化方法[14]。
該方法采用了最小池化和平均池化兩個操作來對點的相似性損失進行正則化，首先使用最小池化找到局部相似性損失D的近似最小值，然后使用兩個平均池化操作進行二次平滑來進行正則化。為了實現上述池化操作，本文將整個相似性損失矩陣排列成了一個6維的矩陣，其中1～3維按照關鍵點采樣位置排列，4～6維按照密集位移采樣d對應位置進行排列。首先對4～6維進行最小池化（最小卷積）和平均池化操作，來找到搜索空間上的近似局部最小值和平滑，接著對1～3維進行平均池化（均值場推理），來對相鄰關鍵點的位移進行約束。由于池化操作可以在GPU上進行并行運算，故可大大提升計算速度。
2 實驗數據與配準結果
2.1 實驗數據
數據來源于公共數據集RESECT[15]，該數據集使用1.5T Siemens Magnetom Avanto對22名患者進行T1w和T2-FLAIR兩種序列MRI掃描，體素大小為1 mm × 1 mm × 1 mm，分辨率為256 × 256 × 192。使用12FLA-L線性超聲探頭分別對患者腫瘤切除前、切除中和切除后進行掃描，體素大小為0.14 mm × 0.14 mm × 0.14 mm～0.24 mm × 0.24 mm × 0.24 mm。本研究選擇手術前T2-FLAIR MR圖像和術中US圖像進行配準實驗，在每對MR圖像和US圖像之間，均事先由人工確定了15～16個對應的標記點，以供各種配準算法進行性能評估和比較。
2.2 關鍵點檢測
角點是圖像亮度變化劇烈的點或圖像邊緣曲線上曲率極大值的點。這些點在保留圖像圖形重要特征的同時，可以有效地減少信息的數據量，使其信息的含量很高，有效地提高了計算速度，有利于圖像的可靠匹配。圖2顯示了US圖像的角點檢測的響應圖、用局部最高響應值篩選后的關鍵點，以及配準后MR圖像中對應的關鍵點。
圖2
關鍵點可視化
Figure2.
Visualization of key points
從圖中可以看出，在US圖像上檢測到的關鍵點，大部分位于重要解剖結構上，如腦脊液、灰質和腫瘤邊緣等。
2.3 實驗結果與分析
2.3.1 配準實驗環境及參數設置
實驗使用Dell precision tower 7910工作站，CPU為Intel Xeon E5-2690 v4，配置NVIDIA GeForce GTX 1080 Ti顯卡的GPU。整個配準流程的計算工作在GPU上實現。采用對應標記點之間的平均歐氏距離，作為平均配準誤差（mean target registration error of landmarks，mTRE），對配準結果進行評估。
配準參數選取如下：在仿射配準進行了3個階段的從粗到細的配準，分別將圖像的分辨率大小采樣為原圖的0.5、0.75和1，關鍵點的非極大值抑制的參數 。該階段只使用了角點特征響應值最高的部分關鍵點，關鍵點數量的參數為64、128、216，密集搜索空間的搜索點數 ，搜索步長 。在彈性配準時，圖像已經進行了仿射的預配準，圖像之間只有微小的誤差，因此只進行了單個階段的精細配準，關鍵點非極大化抑制的參數 。該階段使用了所有的關鍵點，并對關鍵點進行了閾值化處理，閾值化處理采用的是Otsu算法[16]，密集搜索空間的搜索點數 ，搜索步長 。
2.3.2 配準結果和對比分析
為證明本文所提出配準方法MRDOR-KP的配準效果，我們在RESECT數據集的22個病例中進行了驗證。該數據集的初始配準誤差為（5.41 ± 4.19）mm，使用MRDOR-KP方法經過仿射配準后數據的配準誤差為（1.57 ± 0.30）mm，平均配準時間為1.36 s，彈性配準后的配準誤差為（1.40 ± 0.28）mm，平均配準時間為1.53 s。對于仿射配準，MRDOR-KP方法能夠有效地對不同大小形變的圖像進行配準，即使圖像存在比較大的形變，如病例12有19.74 mm的初始誤差，該方法依舊有效。通過彈性形變進行擬合后，絕大部分患者腦圖像的配準誤差在2 mm以內，最大配準誤差為2.02 mm，表明對于所有案例的配準結果都比較好。
為了驗證多分辨率策略和正則化方法的有效性，表1列出了使用/不使用多分辨率分解和正則化策略的配準誤差及運算時間：
表1
不同模塊性能比較
Table1.
Performance comparison of different modules
從表1可以看出，仿射配準階段使用多分辨率策略來進行配準，不僅可以顯著提高配準的速度，并且能夠提升配準精度；彈性配準階段，正則化操作能有效降配準誤差，因采用了并行運算，故正則化計算時間僅增加了0.19 s。
為了評估本文方法的性能，我們將本文的配準結果與CuRIOUS2018挑戰賽[17]排名前3的算法ImFusion[18]、DeedsSSC[19]和cDRAMMS[20]進行了比較，結果如表2所示。
表2
不同方法性能比較
Table2.
Performance comparison of different methods
從表2的配準結果比較可以看出，本文方法成功地將配準時間降低到了2.89 s，比其他方法大幅縮短，并且取得了最低的配準誤差。結果表明，我們的方法可以自動精準地對US與MR圖像進行配準，只需要幾秒就能實現，能夠滿足醫生對于臨床手術精度和時間的需求。
2.3.3 結果可視化
圖3展示了RESECT數據集中三名不同患者（case#12、case#17和case#23）的MR到US的配準結果，為了進一步評估配準結果的好壞，具有臨床經驗的神經外科醫生對配準前后的MR圖像和US圖像對關鍵區域（腦溝和腫瘤邊界）進行了標記。其中白色和藍色箭頭分別指出US和MR的腫瘤邊界，黃色和紅色箭頭分別指出US和MR的腦溝。從配準前后箭頭的對齊效果可以看出經過我們的方法配準后，圖像關鍵區域得到了明顯的對齊。
圖3
部分配準結果
白色和黃色箭頭分別指向US圖像的腫瘤邊界和腦溝，藍色和紅色箭頭分別指向MR圖像的腫瘤邊界和腦溝
Figure3. Partial registration resultsthe white and yellow arrows point to the tumor boundary and sulcus of the US image respectively, and the blue and red arrows point to the tumor boundary and sulcus of the MR image respectively
圖4展示了患者#25在本文算法的配準結果與CuRIOUS2018挑戰賽給出的配準結果的對比。該病例的初始配準誤差為10.06 mm，ImFusion、cDRAMMS和DeedsSSC算法的配準誤差分別為2.12、2.76、1.45 mm，而本文算法的配準誤差僅為1.37 mm。圖中使用箭頭指向了腦溝區域的對齊情況，可以看出我們的算法取得了更好的對齊效果和更小的配準誤差。
圖4
多算法配準結果
白色箭頭指向US圖像的腦溝，右下角為不同算法的mTRE
Figure4. Registration results of multiple algorithmsthe white arrow points to the sulcus of the US image, and the lower right corner is the mTRE of different algorithms
3 結論
本文提出了一種基于關鍵點的多分辨率離散優化兩階段圖像配準算法，采用自相似性上下文描述子定義圖像間相似性測度，并對關鍵點的位移向量進行正則化處理，以減少配準誤差。實驗證明，本文提出的配準方法，具有良好的配準精度和速度，能有效實現MR圖像和US圖像的實時配準。本文的創新點包括：① 采用一種基于角點作為關鍵特征點的策略，根據運動特性不同，分別在仿射配準和彈性配準階段采用不同數量的關鍵點，以充分利用有效的解剖結構信息，并加速配準；② 在仿射配準階段，對圖像進行多分辨率分解，實施由粗到精的配準，以提高配準精度；③ 在彈性配準階段，采取最小卷積和均值場推理策略對位移向量進行正則化處理，有效地減少配準誤差。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：譚振霖負責算法研究、實驗、分析及論文撰寫，郭圣文負責總體設計、方法與論文撰寫指導和審校。
