閆 利，王紫琦，葉志云

武漢大學(xué)測繪學(xué)院，湖北 武漢 430079

多模態(tài)遙感影像能夠彌補單一影像源的不足，提高影像的信息量。有效利用多模態(tài)遙感影像的關(guān)鍵是影像配準技術(shù)。影像配準技術(shù)是多模態(tài)遙感影像融合、目標識別、變化檢測等應(yīng)用必需的前期工作，是各類問題提高精度和有效性的瓶頸[1]。由于不同傳感器的成像機理不同，獲取影像的時間、角度、環(huán)境也不同，多模態(tài)影像的配準研究領(lǐng)域中還存在著很多的困難，特別是圖像或數(shù)據(jù)類別差異大(如光學(xué)與SAR影像)、波段差異大(如可見光與中波紅外影像)等情況下的影像高精度、自動配準技術(shù)的實現(xiàn)，更是存在較大的困難。

對于多模態(tài)遙感影像，灰度特性會發(fā)生非線性甚至非函數(shù)的變化，影像之間的灰度關(guān)系存在統(tǒng)計相關(guān)性及幾何結(jié)構(gòu)的相似性。傳統(tǒng)的遙感影像配準工作量大、自動化程度低。一般的基于灰度的配準方法(如相關(guān)函數(shù)[2])無法成功實現(xiàn)配準。特征描述符在多模態(tài)遙感影像中不具有相同的穩(wěn)健性[3]，基于特征的配準方法(如輪廓特征[4-5]、SUSAN算子[6]、SIFT[7]特征、基于相位一致性的特征[8-9]、擴展SURF描述符[10])也面臨著特征重復(fù)性探測和模態(tài)不變特征描述兩方面的挑戰(zhàn)。而基于統(tǒng)計依賴性的區(qū)域配準方法適用于多模態(tài)影像配準。最常用的互信息[11-13]主要利用影像之間灰度特性的統(tǒng)計相關(guān)性，最早用于醫(yī)學(xué)影像配準，后引入到多模態(tài)圖像配準中。歸一化互信息[14]改進了互信息對于影像重疊率變化敏感的缺點。高維互信息[15]、區(qū)域互信息[16]、特征互信息[17]存在數(shù)據(jù)維度高，計算復(fù)雜等缺陷。交叉累積剩余熵(CCRE)[18]將累積剩余熵引入互信息，具有較強的抗噪性和穩(wěn)定性。除此之外，基于互信息的配準方法還引入了小波變換[19]，實現(xiàn)紅外與光學(xué)遙感影像逐級由粗到精的配準。但這些方法的配準精度無法達到像素級，且均未顧及影像的梯度信息。多模態(tài)影像的地物邊緣屬于幾何特征，具有一定的相似性和不變性，梯度信息能很好地表示影像的邊緣特征，從而使影像配準結(jié)果更加準確。

因此，本文研究了一種顧及灰度和梯度信息的多模態(tài)影像自動配準算法，且滿足空間精度的要求。充分利用多模態(tài)影像的圖像信息進行相似性測量，并采用基于馬爾科夫隨機場(MRF)的非參數(shù)化配準模型，配準變換滿足微分同胚保持影像的拓撲結(jié)構(gòu)[20]。模型上同時考慮了灰度及梯度統(tǒng)計信息，求解方法上對值空間進行離散化，提高收斂速度。試驗表明該算法能較好地完成多模態(tài)遙感影像的自動配準。

1 顧及灰度和梯度信息的多模態(tài)影像配準

本文研究的顧及灰度和梯度信息的多模態(tài)影像配準算法是馬爾科夫隨機場(MRF)和B樣條插值策略的綜合。算法的整體框架是馬爾科夫隨機場(MRF)模型，通過最小化能量函數(shù)來估計最佳變換。首先，將影像進行網(wǎng)格劃分形成均勻控制點網(wǎng)格，考慮影像的灰度和梯度信息進行代價計算，根據(jù)馬爾科夫隨機模型通過最優(yōu)化解算(Fast PD算法[21-22])得到控制網(wǎng)格節(jié)點的位移量，再利用B樣條插值策略計算每個像素點的位移控制圖像變形實施配準。

在圖像配準的具體情況下，MRF模型定義如下[23]：節(jié)點p∈V(具有空間位置xp)對應(yīng)于均勻B樣條網(wǎng)格中的控制點，對于每個節(jié)點，存在一組離散標簽；離散標簽的集合L對應(yīng)于解空間的量化，表示允許的離散位移量，L≡d。標簽集L={l1,…l4n+1}的基數(shù)為4n+1，由沿x、y軸和對角線均勻采樣n個標簽加上原點構(gòu)建；隨機變量對應(yīng)于控制點的位移量；優(yōu)化的能量函數(shù)由兩個項組成：數(shù)據(jù)代價項Dp(lp)，采用塊匹配對源影像及其梯度影像進行相似性測度；正則代價項Ppq(lp,lq)，引入變形場平滑約束的先驗知識，λ是一個標量值，用于權(quán)衡正則項的影響。圖像配準的目標是為每個控制網(wǎng)格節(jié)點分配一個最佳標簽l，使得以下能量函數(shù)最小化

(1)

與此同時，本文算法結(jié)合影像高斯金字塔與變形網(wǎng)格模型的多尺度，高斯金字塔能夠降低計算成本，變形網(wǎng)格模型的多尺度通過將控制點的間距減半來增加變形網(wǎng)格的分辨率，因此，可以恢復(fù)較大的位移，逐步改善配準結(jié)果。對于各層影像金字塔下的各級網(wǎng)格，使用迭代方案來提高離散標簽的效率。保持標簽空間的基數(shù)合理，并在每次迭代時優(yōu)化標簽，并利用上一次迭代的結(jié)果來改進配準影像以捕獲較小的位移。配準流程概述如圖1所示。

1.1 均勻控制點網(wǎng)格劃分與代價計算

本文研究的配準模型是基于MRF的非參數(shù)化配準模型。將影像進行網(wǎng)格劃分形成均勻控制點網(wǎng)格G：[1，K]×[1，L](K和L顯著小于圖像尺寸；控制點之間的距離為δ)，利用多模態(tài)影像的灰度及梯度信息根據(jù)MRF模型進行均勻控制網(wǎng)格節(jié)點的代價計算。其中能量函數(shù)由兩部分組成。第1項為數(shù)據(jù)代價項，用于評估將所有允許的位移應(yīng)用于每個隨機變量的最大似然。對于數(shù)據(jù)代價項，采用塊匹配[24]進行相似性測度，即以每個控制網(wǎng)格節(jié)點為中心的分塊及其梯度影像來進行相似性度量

(2)

(4)

本文采用形態(tài)學(xué)梯度提取多模態(tài)影像的梯度信息。將形態(tài)學(xué)的腐蝕、膨脹、開運算以及閉運算等基本運算用于圖像處理，可構(gòu)成形態(tài)學(xué)梯度算子。常見的形態(tài)學(xué)梯度算子如下

G1[f(x,y)]=(f⊕B)(x,y)-f(x,y)

(5)

G2[f(x,y)]=f(x,y)-(fΘB)(x,y)

(6)

G3[f(x,y)]=(f⊕B)(x,y)-(fΘB)(x,y)

(7)

形態(tài)學(xué)梯度是膨脹圖與腐蝕圖之差，可以檢測圖像邊緣，提取特征信息，在配準時明顯減少相似性測度的計算，同時改善局部極值問題。

相似性測度ρ可以選用對光照變化不敏感的互信息[25]。兩張影像的互信息是通過它們的熵以及聯(lián)合熵來反映它們之間信息的相互包含程度。hI1、hI2、hI1,I2可以用于塊匹配策略的相似性測度計算

miI1,I2(i,k)=hI1(i)+hI2(k)-hI1,I2(i,k)

(8)

式中，hI1,I2可由相應(yīng)灰度的聯(lián)合概率分布PI1,I2計算。相應(yīng)像素的數(shù)量為n。與2D高斯(由?g(i)表示)的卷積進行Parzen估計[26]

(9)

相應(yīng)的，計算hI1,hI2可類似于hI1,I2

(10)

對于正則代價項，引入變形場平滑約束的先驗知識，即假設(shè)控制點的位移場是平滑的，緊密相關(guān)的變量p和q應(yīng)該分配給類似的標簽。采用標簽之間的向量差異的簡單策略，并通過網(wǎng)格距離δ進行歸一化

(11)

模型上同時考慮了灰度以及梯度統(tǒng)計信息，同時采用適用于多模態(tài)遙感數(shù)據(jù)集的塊匹配策略，可以局部評估配準圖像之間的差異，提高計算效率。求解方法上，對值空間進行離散化，提高收斂速度。代價計算是利用多模態(tài)影像的灰度及梯度信息計算基于互信息(MI)的控制網(wǎng)格節(jié)點的匹配成本，其流程圖如圖2所示。其中，代價計算的中間數(shù)值可線性縮放實現(xiàn)可視化，越深的點表示越大的值。

圖1 多模態(tài)影像配準流程概述Fig.1 Overview of presented approach for multimodel image registration

1.2 B樣條的彈性變形模型

通過優(yōu)化MRF模型可以解算控制網(wǎng)格節(jié)點的標簽即位移量。采用基于網(wǎng)格的B樣條的彈性變形模型，可以提供一對一和可逆轉(zhuǎn)換。變形模型的基本思想是通過疊加均勻變形網(wǎng)格G：[1，K]×[1，L]，可以通過控制網(wǎng)格的節(jié)點的位移量并使用插值策略來計算圖像域的其余像素點的位移量來使底層圖像變形

(12)

式中，T(x)為目標影像；d表示一個控制點的位移；η對應(yīng)于控制網(wǎng)格節(jié)點ij對圖像點x的影響的內(nèi)插或加權(quán)函數(shù)。圖像點越近，控制點的影響越大?？刂泣c在圖像上均勻分布行為網(wǎng)格，網(wǎng)格間距為δ。

插值策略采用均勻三次B樣條函數(shù)。令F2為(m+3)×(n+3)控制格網(wǎng)定義的2D均勻三次B樣條函數(shù)。函數(shù)F2由m×n個2D面片組成，每個面片由R2中的4×4控制點決定。通過面片f2來表示二維均勻三次B樣條函數(shù)[27]

(13)

式中，0≤u,v≤1；B0、B1、B2和B3為均勻立方B樣條基函數(shù)；φij為控制點的位移量。

2 試驗與分析

2.1 試驗數(shù)據(jù)

本文共設(shè)置3組多模態(tài)影像配準試驗，具體試驗數(shù)據(jù)如表1所示。配準前，將影像疊加在一起，以棋盤格的方式交替分塊顯示。圖3顯示了3組多模態(tài)影像疊加時的棋盤格，觀察棋盤格的每個分塊邊緣，將未重合的區(qū)域用紅色圓圈突出表示，并從中選擇3個未重合區(qū)域作局部放大顯示。

表1 試驗數(shù)據(jù)介紹

2.2 多模態(tài)影像配準

為了驗證本文算法的性能優(yōu)劣，以3組多模態(tài)影像為研究對象，共進行3組對比試驗：①人工選取10對同名點利用多項式模型進行配準；②只考慮灰度信息的基于MRF的非參數(shù)化配準[23]；③顧及灰度和梯度信息的基于MRF的非參數(shù)化配準。

本文使用了以下試驗配置：三級高斯金字塔，三級變形網(wǎng)格，初始控制點網(wǎng)格間距δ為128像素，其他兩個層級的距離分別設(shè)置為64和32像素。對于各層金字塔下各級網(wǎng)絡(luò)，使用5次迭代方案，采用歸一化互信息作為相似性測度，形態(tài)學(xué)梯度保留物體的邊緣輪廓。標簽集的基數(shù)為41(沿x、y軸和對角線均勻采樣10個標簽加上原點構(gòu)建標簽集)。在第1次迭代中，采樣的最大位移等于0.4×δ，滿足微分同胚要求從而保證影像的拓撲結(jié)構(gòu)不發(fā)生改變[20]。在以下迭代中，最大采樣位移對應(yīng)于前一次迭代的最大位移的0.67。ω的初始值為1。λ為經(jīng)驗選擇，取值在10左右是穩(wěn)健的。面片尺寸等于2δ×2δ。本文研究算法對于中波紅外影像與可見光正射影像配準計算耗時控制在30 min之內(nèi)；對于Terra SAR影像與光學(xué)影像配準計算耗時在12 min左右；對于Worldview-3近紅外與多光譜影像配準計算耗時在9 min左右。與只考慮灰度信息的基于MRF的非參數(shù)化配準方法相比，計算耗時相差不多于1～2 min。

2.2.1 配準結(jié)果定性評價

本文通過棋盤格疊加顯示及興趣區(qū)域的局部圖定性評價3種方法的配準效果。將配準后的多模態(tài)遙感影像疊加在一起，以棋盤格的方式交替分塊顯示。圖4、圖5、圖6分別顯示了中波紅外影像與可見光正射影像、Terra SAR影像與光學(xué)影像、Worldview-3近紅外與多光譜影像配準疊加時的棋盤格。與此同時，在棋盤格疊加顯示結(jié)果上選取3個相同的興趣區(qū)域(黃色矩形圈中顯示)作局部放大顯示，便于比較分析。其中，(A)代表基于人工刺點的多項式模型配準結(jié)果；(B)代表只考慮灰度信息的基于MRF的非參數(shù)化配準結(jié)果；(C)代表顧及灰度和梯度信息的基于MRF的非參數(shù)化配準結(jié)果；(D)代表3種方法的局部區(qū)域?qū)Ρ葓D，每一行表示相同區(qū)域不同方法的配準結(jié)果，每一列表示同一方法不同區(qū)域的配準結(jié)果。觀察多模態(tài)影像疊加顯示的棋盤格的分塊邊緣，可以發(fā)現(xiàn)：

(1) 對于中波紅外影像與可見光正射影像配準、Worldview-3近紅外與多光譜影像配準，基于人工刺點的多項式模型效果最差，而顧及灰度和梯度信息的基于MRF的非參數(shù)化配準效果與只考慮灰度信息的基于MRF的非參數(shù)化配準效果較好，且配準結(jié)果大致相當(dāng)。

(2) 對于Terra SAR影像與光學(xué)影像配準，基于人工刺點的多項式模型與只考慮灰度信息的基于MRF的非參數(shù)化配準方法均不能實現(xiàn)成功配準，而顧及灰度和梯度信息的基于MRF的非參數(shù)化配準結(jié)果的棋盤格分塊邊緣均整齊地重合，配準效果遠優(yōu)于其他兩種方法。

綜上所述，本文算法的棋盤可視化結(jié)果最優(yōu)。

2.2.2 配準結(jié)果定量評價

為準確地評價3種方法的配準效果，本文通過15組同名點對的像素坐標的差異定量評價3種方法的配準效果(表2)?？梢园l(fā)現(xiàn)：

(1) 從同名點對的匹配精度來看，基于人工刺點的多項式模型配準方法的精度最差，而顧及灰度及梯度信息的基于MRF的非參數(shù)化配準模型的精度最優(yōu)，3組多模態(tài)影像的x、y方向的平均配準誤差均小于1個像素，最大誤差小于2個像素。

(2) 本文算法對于3組多模態(tài)影像的配準精度均有明顯改進，尤其是Terra SAR影像與光學(xué)影像的配準。在其余兩種方法均不能實現(xiàn)成功配準的情況下，本文研究的算法針對Terra SAR影像與光學(xué)影像可達到像素級的配準精度。

綜上所述，對于多組多模態(tài)影像的配準，本文研究的算法最優(yōu)且滿足像素級配準精度要求。

表2 3種方法的配準結(jié)果及誤差

2.3 大形變影像的配準試驗

上述試驗中，多模態(tài)影像的初始最大視差約為50個像素，為了進一步驗證本文研究的算法的適應(yīng)性，即測試本文算法針對較大視差的多模態(tài)影像是否有效，將配準前的中波紅外影像進行人工平移，對較大初始視差的中波紅外數(shù)據(jù)進行了試驗。

根據(jù)位移離散空間的最大標簽值和金字塔影像級數(shù)，本文設(shè)置如圖7所示的3組試驗：對中波紅外影像進行x方向、y方向、對角線方向各200個像素的平移，將平移后的中波紅外影像與可見光正射影像進行配準。試驗配置與上述試驗相同，采用歸一化互信息作為相似性測度。配準

前的平移后的中波紅外影像及配準后的棋盤格疊加顯示如圖7所示，配準精度見表3。

由配準結(jié)果可知，棋盤格分塊的邊緣均整齊地重合。從同名檢查點的匹配精度來看，x、y方向的平均配準精度均小于1像素，最大誤差小于2個像素。相比與較小形變的多模態(tài)影像配準，同名點對的匹配精度基本相當(dāng)。因此針對較大形變的多模態(tài)影像，本文研究的算法也具有較好的配準效果，配準精度滿足需求。由此說明本文研究的算法針對較大形變的多模態(tài)影像配準也具有一定的適用性，利用高精度的POS定向數(shù)據(jù)與快速拼接能夠滿足大部分多模態(tài)影像的配準要求。

圖2 代價計算流程圖Fig.2 Overview of cost calculation visualization

圖3 配準前多模態(tài)影像的棋盤格疊加顯示及其局部圖Fig. 3 Chessboard overlay and local graphs of multimodel images before registration

圖4 中波紅外影像與可見光正射影像配準后棋盤格疊加顯示及其局部圖Fig.4 Chessboard overlay and local graph of mid-infrared image and visible orthophoto after registration

圖5 Terra SAR影像與光學(xué)影像配準后棋盤格疊加顯示及其局部圖Fig.5 Chessboard overlay and local graph of Terra SAR image and optical image after registration

圖6 Worldview-3近紅外與多光譜影像配準后棋盤格疊加顯示及其局部圖Fig.6 Chessboard overlay and local graph of worldview-3 near infrared and multi-spectral images after registration

3 結(jié) 論

本文研究的基于MRF的非參數(shù)化配準模型同時考慮了灰度以及梯度統(tǒng)計信息，能夠?qū)崿F(xiàn)多種多模態(tài)圖像的自動配準，具有較好的穩(wěn)定性和配準精度。文中進行了多組多模態(tài)影像配準試驗，試驗結(jié)果證明：本文算法的平均配準誤差小于1個像素，最大配準誤差小于2個像素。該算法利用互信息進行相似性測度時需要計算多模態(tài)影像的聯(lián)合概率分布，因此待配準影像的初始視差范圍是有限的，對于具有顯著旋轉(zhuǎn)、尺度變化和視差的多模態(tài)數(shù)據(jù)的配準，有待進一步研究。

[1] 倪國強, 劉瓊. 多源圖像配準技術(shù)分析與展望[J]. 光電工程, 2004, 31(9): 1-6. NI Guoqiang, LIU Qiong. Analysis and Prospect of Multi-Source Image Registration Techniques[J]. Opto-Electronic Engineering, 2004, 31(9): 1-6.

[2] BROWN L G. A Survey of Image Registration Techniques[J]. ACM Computing Surveys, 1992, 24(4): 325-376.

[3] FANBin, HUO Chunlei, PAN Chunhong, et al. Registration of Optical and SAR Satellite Images by Exploring the Spatial Relationship of the Improved SIFT[J]. IEEE Geoscience and Remote Sensing Letters, 2013,10(4):657-661.

[4] LI Hui, MANJUNATH B S,MITRA S K.A Contour-Based Approach to Multisensor Image Registration[J]. IEEE Transactions on Image Processing, 1995,4(3):320-334.

[5] LI HH,ZHOU Yitong.Automatic Visual/IR Image Registration[J]. Optical Engineering, 1996,35(2):391-400.

[6] 張遷, 劉政凱, 龐彥偉, 等. 基于SUSAN算法的航空影像的自動配準[J]. 測繪學(xué)報, 2003, 32(3): 245-250. ZHANG Qian, LIU Zhengkai, PANG Yanwei, et al. Automatic Registration of Aerophotos Based on SUSAN Operator[J]. Acta Geodaetica et Cartographica Sinica, 2003, 32(3): 245-250.

[7] KELMANA, SOFKA M,STEWART C V.Keypoint Descriptors for Matching Across Multiple Image Modalities and Non-Linear Intensity Variations[C]∥Proceedings of 2007 IEEE Conference on Computer Vision and Pattern Recognition.Minneapolis, MN, USA: IEEE,2007.

[8] KOVESI P.Image Featuresfrom Phase Congruency[J]. Videre: Journal of Computer Vision Research, 1999,1(3): 1-26.

[9] WONG A,CLAUSI D A.AISIR: Automated Inter-Sensor/Inter-Band Satellite Image Registration Using Robust Complex Wavelet Feature Representations[J]. Pattern Recognition Letters, 2010,31(10):1160-1167.

[10] 羅楠, 孫權(quán)森, 耿蕾蕾, 等. 一種擴展SURF描述符及其在遙感圖像配準中的應(yīng)用[J]. 測繪學(xué)報, 2013, 42(3): 383-388. LUO Nan, SUN Quansen, GENG Leilei, et al. An Extended SURF Descriptor and Its Application in Remote Sensing Images Registration[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(3): 383-388.

[11] VIOLA P,WELLS III W M. Alignment by Maximization of Mutual Information[J]. International Journal of Computer Vision,1997, 24(2):137-154.

[12] COLLIGNONA, MAES F, DELAERE D, et al. Automated Multi-Modality Image Registration based on Information Theory[C]∥Proceedings of the International Conference on Information Processing in Medical Imaging.Ile de Berder, France:Kluwer Academic Publishers,1995.

[13] MAES F, COLLIGNON A, VANDERMEULEN D, et al. Multimodality Image Registration by Maximization of Mutual Information[J]. IEEE Transactions on Medical Imaging,1997, 16(2): 187-198.

[14] STUDHOLME C, HILL D L G, HAWKES D J. An Overlap Invariant Entropy Measure of 3D Medical Image Alignment[J]. Pattern Recognition,1999, 32(1):71-86.

[15] RUECKERT D,CLARKSONM J, HILL D L G, et al. Non-Rigid Registration Using Higher-Order Mutual Information[C]∥Proceedings of SPIE-Medical Imaging 2000: Image Processing.San Diego, CA, United States: SPIE, 2000: 438-447.

[16] RUSSAKOFF D B, TOMASI C, ROHLFING T, et al. Image Similarity Using Mutual Information of Regions[C]∥Proceedings of the 8th European Conference on Computer Vision.Prague, Czech Republic: Springer, 2004: 596-607.

[18] WANG Fei, VEMURI B C. Non-Rigid Multi-Modal Image Registration Using Cross-CumulativeResidual Entropy[J]. International Journal of Computer Vision,2007, 74(2):201-215.

[19] FAN Xiaofeng. Automatic Registration of Multi-Modal Airborne Imagery[D].Rochester:Rochester Institute of Technology, 2011.

[20] 閆德勤, 劉彩鳳, 劉勝藍, 等. 大形變微分同胚圖像配準快速算法[J]. 自動化學(xué)報, 2015, 41(8): 1461-1470. YAN Deqin, LIU Caifeng, LIU Shenglan, et al.A Fast Image Registration Algorithm for DiffeomorphicImage with Large Deformation[J]. ActaAutomaticaSinica, 2015, 41(8): 1461-1470.

[21] STREKALOVSKIY E, CREMERS D. Real-Time Minimization of the Piecewise Smooth Mumford-Shah Functional[C]∥Proceedings of the 13th European Conference on Computer Vision. Zurich, Switzerland:Springer, 2014: 127-141.

[22] KOMODAKIS N, TZIRITAS G, PARAGIOS N. Performance vs Computational Efficiency for Optimizing Single and Dynamic MRFs: Setting the State of the Art with Primal-Dual Strategies[J]. Computer Vision and Image Understanding,2008, 112(1): 14-29.

[23] KARANTZALOSK, SOTIRAS A,PARAGIOS N. Efficient and Automated Multimodal Satellite Data Registration through MRFs and Linear Programming[C]∥Proceedings of 2014 IEEE Conference on Computer Vision and Pattern Recognition Workshops.Columbus, OH, USA: IEEE, 2014.

[24] 張劍清, 潘勵, 王樹根. 攝影測量學(xué)[M]. 2版. 武漢: 武漢大學(xué)出版社, 2009. ZHANG Jianqing, PAN Li, WANG Shugen. Photogrammetry[M]. 2nd ed. Wuhan: Wuhan University Press, 2009.

[25] HIRSCHMULLERH.Stereo Processing by Semiglobal Matching and Mutual Information[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2008,30(2):328-341.

[26] KIM J, KOLMOGOROV V,ZABIH R. Visual Correspondence Using Energy Minimization and Mutual Information[C]∥Proceedings of the 9th IEEE International Conference on Computer Vision. Nice, France: IEEE, 2003.

[27] CHOIY,LEE S.Injectivity Conditions of 2D and 3D Uniform Cubic B-Spline Functions[J]. Graphical Models, 2000,62(6):411-427.