摘 要： 將CUDA技術(shù)應(yīng)用于數(shù)字圖像拼接領(lǐng)域，闡述了圖像拼接的基本理論及其關(guān)鍵技術(shù)、多分辨率圖像融合的關(guān)鍵算法以及CUDA技術(shù)的基本原理和開發(fā)方法，并編寫了軟件以實現(xiàn)圖像快速拼接。采用對于尺度具有魯棒性的SIFT 算法進(jìn)行特征點的提取與匹配，使用穩(wěn)健的RANSAC 算法求出圖像間變換矩陣的值,并將圖像映射到拼接平面，最后使用基于CUDA的SIFT算法實現(xiàn)了圖像的無縫拼接。該方法提高了圖像拼接的效率，克服了傳統(tǒng)圖像拼接方法因計算量大而等待時間長的缺點。實驗結(jié)果表明,CUDA在數(shù)字圖像處理的實際應(yīng)用中卓有成效，有廣闊的應(yīng)用前景。
關(guān)鍵詞： CUDA; 圖像拼接; SIFT; 多分辨率融合
圖像拼接是計算機(jī)視覺領(lǐng)域的一個重要分支,它是一種將多幅相關(guān)的部分重疊圖像進(jìn)行無縫拼接從而獲得寬視角圖像的技術(shù)。利用計算機(jī)進(jìn)行匹配，將多幅具有重疊關(guān)系的圖像拼合成為一幅具有更大視野范圍的圖像，這就是圖像拼接的目的。CUDA是英偉達(dá)(NVIDIA)公司傾力開發(fā)和推廣的并行計算架構(gòu)，該架構(gòu)通過利用圖形處理器(GPU)的處理能力，能夠大幅提升計算性能。隨著微軟Windows 7與蘋果Snow Leopard操作系統(tǒng)的問世，GPU計算必將成為主流。本文基于SIFT算法，使用最新的CUDA并行計算技術(shù)重編算法并編制軟件，不僅可以克服傳統(tǒng)圖像拼接技術(shù)中的局限性(如光照、尺度變化的影響等)，實現(xiàn)光照和尺度變化條件下的多視角無縫圖像拼接，而且將提高圖像拼接的速度和效率。
1 數(shù)字圖像拼接的基本理論和方法
圖像拼接的基本流程如圖1所示，主要分為圖像預(yù)處理、圖像配準(zhǔn)和圖像融合與邊界平滑3個步驟。圖像預(yù)處理主要指對圖像進(jìn)行幾何畸變校正和噪聲點的抑制等,使參考圖像和待拼接圖像不存在明顯的幾何畸變。圖像預(yù)處理主要是為下一步圖像配準(zhǔn)做準(zhǔn)備，使圖像質(zhì)量能夠滿足圖像配準(zhǔn)的要求。圖像配準(zhǔn)主要指對參考圖像和待拼接圖像中的匹配信息進(jìn)行提取，在提取出的信息中尋找最佳的匹配，完成圖像間的對齊。圖像拼接的成功與否主要是圖像的配準(zhǔn)。待拼接的圖像之間可能存在平移、旋轉(zhuǎn)和縮放等多種變換或者大面積的同色區(qū)域等很難匹配的情況，一個好的圖像配準(zhǔn)算法應(yīng)該能夠在各種情況下準(zhǔn)確找到圖像間的對應(yīng)信息，將圖像對齊。圖像融合指在完成圖像匹配以后對圖像進(jìn)行縫合，并對縫合的邊界進(jìn)行平滑處理，使縫合自然過渡。由于任何兩幅相鄰圖像在采集條件上都不可能做到完全相同，因此，對于一些本應(yīng)該相同的圖像特性(如圖像的光照特性等)，在兩幅圖像中就不會表現(xiàn)得完全一樣。圖像拼接縫隙就是從一幅圖像的圖像區(qū)域過渡到另一幅圖像的圖像區(qū)域時，由于圖像中的某些相關(guān)特性發(fā)生了躍變而產(chǎn)生的。圖像融合就是使圖像間的拼接縫隙不明顯，拼接更自然。
2 CUDA技術(shù)的基本原理和開發(fā)方法
圖像處理的本質(zhì)是大規(guī)模矩陣運(yùn)算，特別適合并行處理。但CPU通用計算很難利用該特性。與此相反，GPU在并行數(shù)據(jù)運(yùn)算上具有強(qiáng)大的計算能力，特別適合作運(yùn)算符相同而運(yùn)算數(shù)據(jù)不同的運(yùn)算，當(dāng)執(zhí)行具有高運(yùn)算密度的多數(shù)據(jù)元素時，內(nèi)存訪問的延遲可以被忽略。
CUDA編程模型將CPU作為主機(jī)（Host），GPU作為協(xié)處理器(Coprocessor)或設(shè)備(Device)，一個系統(tǒng)中可以存在多個設(shè)備。在這個模型中，CPU與GPU共同工作，CPU負(fù)責(zé)邏輯性強(qiáng)的事務(wù)處理和串行計算，GPU則專注于執(zhí)行高度線程化的并行處理任務(wù)。圖 2顯示了CUDA的流程架構(gòu)。
CUDA對內(nèi)存的操作與一般的C程序基本相同，操作顯存則需要調(diào)用CUDA API中的存儲器管理函數(shù)。一旦確定好程序中的并行部分，就可以將這部分計算交給GPU。運(yùn)行在GPU上的CUDA并行計算函數(shù)稱為Kernel，即內(nèi)核函數(shù)，它并不是一個完整的程序，而是CUDA程序中可以被并行執(zhí)行的步驟。內(nèi)核函數(shù)必須通過_global_函數(shù)類型限定符定義，且只能在主機(jī)端代碼中調(diào)用。
CUDA線程結(jié)構(gòu)如圖3所示。Kernel以線程網(wǎng)格(Grid)為組織形式，每個Grid由若干個線程塊(Block)組成，每個Block又由若干個線程(Thread)組成。在程序運(yùn)行過程中，Kernel是以Block為單位執(zhí)行的，Grid只是一系列可以被執(zhí)行的Block的集合。不同Block是并行執(zhí)行的，沒有執(zhí)行的先后順序，而且相互無法通信。
CUDA軟件體系由CUDA Library、CUDA Runtime API和CUDA Driver API構(gòu)成，如圖 4所示。CUDA的核心是CUDA C語言，需要通過nvcc編譯器進(jìn)行編譯。編譯得到的僅是GPU端的代碼，要在GPU上分配顯存并啟動Kernel就需要借助CUDA Runtime API或者CUDA Driver API來實現(xiàn)。CUDA Runtime API和CUDA Driver API 提供了實現(xiàn)設(shè)備管理、上下文管理、存儲器管理、代碼塊管理和執(zhí)行控制等操作的應(yīng)用程序接口。CUDA Runtime API在CUDA Driver API的基礎(chǔ)上進(jìn)行了封裝，使得編程方便代碼簡潔, 因此通常采用CUDA Runtime API進(jìn)行項目開發(fā)。
3 SIFT特征匹配算法
SIFT算法首先在尺度空間進(jìn)行特征檢測，并確定關(guān)鍵點(Key Points)的位置和關(guān)鍵點所處的尺度，然后使用關(guān)鍵點鄰域梯度的主方向作為該點的方向特征，以實現(xiàn)算子對尺度和方向的無關(guān)性。
SIFT特征匹配算法包括兩個階段：(1)SIFT特征的生成，即從多幅待匹配圖像中提取出對尺度縮放、旋轉(zhuǎn)和亮度變化無關(guān)的特征向量；(2)SIFT特征向量的匹配。
在實際的尺度不變特征點提取中，SIFT算法將圖像金字塔引入了尺度空間。首先采用不同尺度因子的高斯核對圖像進(jìn)行卷積以得到圖像的不同尺度空間，將這一組圖像作為金字塔圖像的第一階。接著對其中的2倍尺度圖像(相對于該階第一幅圖像的2倍尺度)以2倍像素距離進(jìn)行下采樣來得到金字塔圖像第二階的第一幅圖像，對該圖像采用不同尺度因子的高斯核進(jìn)行卷積，以獲得金字塔圖像第二階的一組圖像。以此類推，獲得高斯金字塔圖像。每一階相鄰的高斯圖像相減，就得到了高斯差分圖像，即DoG圖像。通過擬和三維二次函數(shù)以精確確定關(guān)鍵點的位置和尺度，同時去除低對比度的關(guān)鍵點和不穩(wěn)定的邊緣響應(yīng)點(因為DoG算子會產(chǎn)生較強(qiáng)的邊緣響應(yīng))，以增強(qiáng)匹配穩(wěn)定性、提高抗噪聲能力。利用關(guān)鍵點鄰域像素的梯度方向分布特性為每個關(guān)鍵點指定方向參數(shù)，使算子具備旋轉(zhuǎn)不變性，生成SIFT特征向量。
接下來以關(guān)鍵點為中心取8×8的窗口。圖5(a)的中央黑點為當(dāng)前關(guān)鍵點的位置，每個小格代表關(guān)鍵點鄰域所在尺度空間的一個像素，箭頭方向代表該像素的梯度方向，箭頭長度代表梯度模值，圖中圈代表高斯加權(quán)的范圍（越靠近關(guān)鍵點的像素,其梯度方向信息貢獻(xiàn)越大）。然后在每4×4的小塊上計算8個方向的梯度方向直方圖，繪制每個梯度方向的累加值，即可形成一個種子點，如圖5(b)所示。圖5(b)中,一個關(guān)鍵點由2×2共4個種子點組成，每個種子點有8個方向向量信息。這種鄰域方向性信息聯(lián)合的思想增強(qiáng)了算法抗噪聲的能力，同時對于含有定位誤差的特征匹配也提供了較好的容錯性。
4 基于CUDA的圖像拼接軟件的設(shè)計
4.1 Host端實現(xiàn)
程序的Host端由C++編寫，負(fù)責(zé)控制整個程序的執(zhí)行流程、所有供CPU 和GPU 所用的數(shù)據(jù)的分配管理以及Device 端模塊的調(diào)用。界面使用最基本的Windows SDK編寫。
在數(shù)據(jù)初始化階段，包含所有之后處理步驟中所需要的圖像數(shù)據(jù)對象的生成，將輸入圖像作為高斯金字塔底層，通過系統(tǒng)中的PYRAMID_LEVEL宏指定金字塔的層數(shù)，在輸入圖像的尺寸基礎(chǔ)上循環(huán)計算各層金字塔圖像的分辨率，并對圖像進(jìn)行初始化。由于所有的圖像數(shù)據(jù)需要在設(shè)備端處理，使用cudaMallocPitch函數(shù)分配數(shù)據(jù)地址空間，數(shù)據(jù)結(jié)構(gòu)不再是OpenCV中的IplImage，而是GPU可以識別的uchar數(shù)組類型。
使用cudaMemcpy2D函數(shù)將IplImage結(jié)構(gòu)中的原始數(shù)據(jù)復(fù)制到相應(yīng)高斯金字塔的最底層，也就是uchar數(shù)組的第一個元素，供Device端函數(shù)使用。隨后進(jìn)行Kernel函數(shù)調(diào)用，對于每一個需要處理的金字塔層，Host端發(fā)起一次Kernel調(diào)用。例如:
reduce<<<(int)ceil((float)(imageSize[l]/THREAD_NUM)), HREAD_NUM>>>(lGaussianData[l+1],rGaussianData[l+1],lGaussianData[l], rGaussianData[l], lLaplacianData[l], rLaplacianData[l], stride[l+1], stride[l], width[l+1], height[l+1], width[l], height[l]);
4.2 Device端實現(xiàn)
主要Device函數(shù)如下：
(1) reduce（）函數(shù)對左、右圖和掩碼圖像各完成一次reduce操作，生成下一層高斯金字塔圖像。reduce變換按照前文所述的方法對目標(biāo)層的金字塔圖像進(jìn)行逐像素處理，每一個目標(biāo)像素的顏色值按一定的權(quán)重值對原始圖像中的一個5×5子塊進(jìn)行計算求得。
_global_ void reduce(uchar* lGaussianDataSrc, uchar* rGaussianDataSrc, uchar* mGaussianDataSrc,uchar* lGaussianDataDst, uchar* rGaussianDataDst, uchar* mGaussianDataDst, size_t strideSrc, size_t strideDst, int srcWidth, int srcHeight, int gauWidth, int gauHeight)；
(2) expand_and_minus（）函數(shù)對左、右圖像各完成一次expand操作和減法操作，生成下一層拉普拉斯金字塔圖像。expand變換相當(dāng)于reduce變換的逆過程，它對目標(biāo)層的金字塔圖像進(jìn)行逐像素處理，每一個目標(biāo)像素的顏色值也是按reduce變換中所使用的權(quán)重值對原始圖像中的一個5×5子塊進(jìn)行計算求得的。
_global_ void expand_and_minus(uchar*lGaussianDataSrcH,uchar* rGaussianDataSrcH,uchar*lGaussianDataSrcL,uchar*rGau-
ssianDataSrcL,uchar* lLaplacianDataDst,uchar*rLaplacianDataDst,size_t strideSrc,size_t strideDst,int srcWidth, int srcHeight, int expWidth, int expHeight)；
(3) blend( )函數(shù)根據(jù)掩碼圖像的高斯金字塔以及左、右圖像的拉普拉斯金字塔合成當(dāng)前層的目標(biāo)圖像的拉普拉斯金字塔,所有像素值均以掩碼圖像的高斯金字塔為權(quán)重而求得。
_global_ void blend(uchar* mGaussianDataSrc, uchar* lLaplacianDataSrc, uchar* rLaplacianDataSrc,uchar* sLaplacianDataDst, size_t stride,int lapWidth, int lapHeight)
(4) collapse（）函數(shù)對圖像的拉普拉斯金字塔分別完成一次expand操作和累加操作,本質(zhì)上等同于expand操作，兩者的基本算法是相同的。不同點在于expand模塊用于各層高斯金字塔的expand操作，從而生成各層拉普拉斯金字塔，而collapse函數(shù)則用于整個融合過程最后的圖像重構(gòu)步驟，將各層已經(jīng)求得的拉普拉斯金字塔作擴(kuò)展和累加操作，生成最后的拼接圖像。
_global_ void collapse(uchar* sLaplacianDataSrc, uchar* sExpandDataSrc, uchar* sExpandDataDst, size_t strideSrc, size_t strideDst, int srcWidth, int srcHeight, int expWidth, int expHeight)
本文借助于SIFT特征對于旋轉(zhuǎn)和尺度的不變性以及對于噪聲干擾良好的魯棒性進(jìn)行圖像拼接與匹配，使用CUDA技術(shù)簡單地對多分辨率融合算法進(jìn)行了優(yōu)化，提高了其執(zhí)行效率和速度。編寫了界面化的Demo程序，實現(xiàn)了基本的圖像拼接功能。
參考文獻(xiàn)
[1] 譚康.圖像拼接技術(shù)與實現(xiàn)[D].南京:南京理工大學(xué), 2006.
[2] HARRIS C, STEPHENS M. A combined corner and edge detector[C]. Proceedings of the 4th Alvey Vision Conference, 1988:147-151.
[3] 張小洪，李博，楊丹．一種新的Harris多尺度角點檢測[J]. 電子與信息學(xué)報，2007（7）：1735-1738.
[4] LOWE D G. Object recognition from local scale-invariant features[C]. The Proceedings of the Seventh IEEE International Conference on Computer Vision, 1999(2):1150-1157.
[5] 騫森，朱劍英. 基于改進(jìn)的SIFT特征的圖像雙向匹配算法[J]. 機(jī)械科學(xué)與技術(shù),2007(9):1179-1182.
[6] Peng Xiaoming, Ding Mingyue, Zhou Chengping, et al. Improved approach for object location under affine transformation using the Hausdorff distance[J].Optical Engineering, 2003,42(10):2794-2795.
[7] 張毓晉. 圖像工程(上冊)圖像處理(第2版)[M]. 北京：清華大學(xué)出版社，2006.
[8] LINDEBERG T． Detecting salient blob. like image structures and their scales with a scale-space primal sketch[J]．International Journal of Computer Vision,1993，11(3):283-318.
[9] NVIDIA. NVIDIA CUDA Programming Guide[Z].
[10] 張舒, 褚艷利. GPU高性能運(yùn)算之CUDA[M]. 北京：中國水利水電出版社，2009.