摘  要： 將CUDA技術(shù)應(yīng)用于數(shù)字圖像拼接領(lǐng)域，闡述了圖像拼接的基本理論及其關(guān)鍵技術(shù)、多分辨率圖像融合的關(guān)鍵算法以及CUDA技術(shù)的基本原理和開發(fā)方法，并編寫了軟件以實現(xiàn)圖像快速拼接。采用對于尺度具有魯棒性的SIFT 算法進行特征點的提取與匹配，使用穩(wěn)健的RANSAC 算法求出圖像間變換矩陣的值,并將圖像映射到拼接平面，最后使用基于CUDA的SIFT算法實現(xiàn)了圖像的無縫拼接。該方法提高了圖像拼接的效率，克服了傳統(tǒng)圖像拼接方法因計算量大而等待時間長的缺點。實驗結(jié)果表明,CUDA在數(shù)字圖像處理的實際應(yīng)用中卓有成效，有廣闊的應(yīng)用前景。
關(guān)鍵詞： CUDA; 圖像拼接; SIFT; 多分辨率融合

    圖像拼接是計算機視覺領(lǐng)域的一個重要分支,它是一種將多幅相關(guān)的部分重疊圖像進行無縫拼接從而獲得寬視角圖像的技術(shù)。利用計算機進行匹配，將多幅具有重疊關(guān)系的圖像拼合成為一幅具有更大視野范圍的圖像，這就是圖像拼接的目的。CUDA是英偉達(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ù)處理主要指對圖像進行幾何畸變校正和噪聲點的抑制等,使參考圖像和待拼接圖像不存在明顯的幾何畸變。圖像預(yù)處理主要是為下一步圖像配準(zhǔn)做準(zhǔn)備，使圖像質(zhì)量能夠滿足圖像配準(zhǔn)的要求。圖像配準(zhǔn)主要指對參考圖像和待拼接圖像中的匹配信息進行提取，在提取出的信息中尋找最佳的匹配，完成圖像間的對齊。圖像拼接的成功與否主要是圖像的配準(zhǔn)。待拼接的圖像之間可能存在平移、旋轉(zhuǎn)和縮放等多種變換或者大面積的同色區(qū)域等很難匹配的情況，一個好的圖像配準(zhǔn)算法應(yīng)該能夠在各種情況下準(zhǔn)確找到圖像間的對應(yīng)信息，將圖像對齊。圖像融合指在完成圖像匹配以后對圖像進行縫合，并對縫合的邊界進行平滑處理，使縫合自然過渡。由于任何兩幅相鄰圖像在采集條件上都不可能做到完全相同，因此，對于一些本應(yīng)該相同的圖像特性(如圖像的光照特性等)，在兩幅圖像中就不會表現(xiàn)得完全一樣。圖像拼接縫隙就是從一幅圖像的圖像區(qū)域過渡到另一幅圖像的圖像區(qū)域時，由于圖像中的某些相關(guān)特性發(fā)生了躍變而產(chǎn)生的。圖像融合就是使圖像間的拼接縫隙不明顯，拼接更自然。

2 CUDA技術(shù)的基本原理和開發(fā)方法
    圖像處理的本質(zhì)是大規(guī)模矩陣運算，特別適合并行處理。但CPU通用計算很難利用該特性。與此相反，GPU在并行數(shù)據(jù)運算上具有強大的計算能力，特別適合作運算符相同而運算數(shù)據(jù)不同的運算，當(dāng)執(zhí)行具有高運算密度的多數(shù)據(jù)元素時，內(nèi)存訪問的延遲可以被忽略。
    CUDA編程模型將CPU作為主機（Host），GPU作為協(xié)處理器(Coprocessor)或設(shè)備(Device)，一個系統(tǒng)中可以存在多個設(shè)備。在這個模型中，CPU與GPU共同工作，CPU負(fù)責(zé)邏輯性強的事務(wù)處理和串行計算，GPU則專注于執(zhí)行高度線程化的并行處理任務(wù)。圖 2顯示了CUDA的流程架構(gòu)。

   CUDA對內(nèi)存的操作與一般的C程序基本相同，操作顯存則需要調(diào)用CUDA API中的存儲器管理函數(shù)。一旦確定好程序中的并行部分，就可以將這部分計算交給GPU。運行在GPU上的CUDA并行計算函數(shù)稱為Kernel，即內(nèi)核函數(shù)，它并不是一個完整的程序，而是CUDA程序中可以被并行執(zhí)行的步驟。內(nèi)核函數(shù)必須通過_global_函數(shù)類型限定符定義，且只能在主機端代碼中調(diào)用。
    CUDA線程結(jié)構(gòu)如圖3所示。Kernel以線程網(wǎng)格(Grid)為組織形式，每個Grid由若干個線程塊(Block)組成，每個Block又由若干個線程(Thread)組成。在程序運行過程中，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編譯器進行編譯。編譯得到的僅是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ǔ)上進行了封裝，使得編程方便代碼簡潔, 因此通常采用CUDA Runtime API進行項目開發(fā)。

 3 SIFT特征匹配算法
    SIFT算法首先在尺度空間進行特征檢測，并確定關(guān)鍵點(Key Points)的位置和關(guān)鍵點所處的尺度，然后使用關(guān)鍵點鄰域梯度的主方向作為該點的方向特征，以實現(xiàn)算子對尺度和方向的無關(guān)性。
    SIFT特征匹配算法包括兩個階段：(1)SIFT特征的生成，即從多幅待匹配圖像中提取出對尺度縮放、旋轉(zhuǎn)和亮度變化無關(guān)的特征向量；(2)SIFT特征向量的匹配。
    在實際的尺度不變特征點提取中，SIFT算法將圖像金字塔引入了尺度空間。首先采用不同尺度因子的高斯核對圖像進行卷積以得到圖像的不同尺度空間，將這一組圖像作為金字塔圖像的第一階。接著對其中的2倍尺度圖像(相對于該階第一幅圖像的2倍尺度)以2倍像素距離進行下采樣來得到金字塔圖像第二階的第一幅圖像，對該圖像采用不同尺度因子的高斯核進行卷積，以獲得金字塔圖像第二階的一組圖像。以此類推，獲得高斯金字塔圖像。每一階相鄰的高斯圖像相減，就得到了高斯差分圖像，即DoG圖像。通過擬和三維二次函數(shù)以精確確定關(guān)鍵點的位置和尺度，同時去除低對比度的關(guān)鍵點和不穩(wěn)定的邊緣響應(yīng)點(因為DoG算子會產(chǎn)生較強的邊緣響應(yī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)鍵點的像素,其梯度方向信息貢獻越大）。然后在每4×4的小塊上計算8個方向的梯度方向直方圖，繪制每個梯度方向的累加值，即可形成一個種子點，如圖5(b)所示。圖5(b)中,一個關(guān)鍵點由2×2共4個種子點組成，每個種子點有8個方向向量信息。這種鄰域方向性信息聯(lián)合的思想增強了算法抗噪聲的能力，同時對于含有定位誤差的特征匹配也提供了較好的容錯性。

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)計算各層金字塔圖像的分辨率，并對圖像進行初始化。由于所有的圖像數(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ù)使用。隨后進行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)層的金字塔圖像進行逐像素處理，每一個目標(biāo)像素的顏色值按一定的權(quán)重值對原始圖像中的一個5×5子塊進行計算求得。
       _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)層的金字塔圖像進行逐像素處理，每一個目標(biāo)像素的顏色值也是按reduce變換中所使用的權(quán)重值對原始圖像中的一個5×5子塊進行計算求得的。
    _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)求得的拉普拉斯金字塔作擴展和累加操作，生成最后的拼接圖像。
    _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)和尺度的不變性以及對于噪聲干擾良好的魯棒性進行圖像拼接與匹配，使用CUDA技術(shù)簡單地對多分辨率融合算法進行了優(yōu)化，提高了其執(zhí)行效率和速度。編寫了界面化的Demo程序，實現(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] 騫森，朱劍英. 基于改進的SIFT特征的圖像雙向匹配算法[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高性能運算之CUDA[M]. 北京：中國水利水電出版社，2009.