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

    圖像拼接是計算機視覺領域的一個重要分支,它是一種將多幅相關的部分重疊圖像進行無縫拼接從而獲得寬視角圖像的技術。利用計算機進行匹配，將多幅具有重疊關系的圖像拼合成為一幅具有更大視野范圍的圖像，這就是圖像拼接的目的。CUDA是英偉達(NVIDIA)公司傾力開發和推廣的并行計算架構，該架構通過利用圖形處理器(GPU)的處理能力，能夠大幅提升計算性能。隨著微軟Windows 7與蘋果Snow Leopard操作系統的問世，GPU計算必將成為主流。本文基于SIFT算法，使用最新的CUDA并行計算技術重編算法并編制軟件，不僅可以克服傳統圖像拼接技術中的局限性(如光照、尺度變化的影響等)，實現光照和尺度變化條件下的多視角無縫圖像拼接，而且將提高圖像拼接的速度和效率。
1 數字圖像拼接的基本理論和方法
    圖像拼接的基本流程如圖1所示，主要分為圖像預處理、圖像配準和圖像融合與邊界平滑3個步驟。圖像預處理主要指對圖像進行幾何畸變校正和噪聲點的抑制等,使參考圖像和待拼接圖像不存在明顯的幾何畸變。圖像預處理主要是為下一步圖像配準做準備，使圖像質量能夠滿足圖像配準的要求。圖像配準主要指對參考圖像和待拼接圖像中的匹配信息進行提取，在提取出的信息中尋找最佳的匹配，完成圖像間的對齊。圖像拼接的成功與否主要是圖像的配準。待拼接的圖像之間可能存在平移、旋轉和縮放等多種變換或者大面積的同色區域等很難匹配的情況，一個好的圖像配準算法應該能夠在各種情況下準確找到圖像間的對應信息，將圖像對齊。圖像融合指在完成圖像匹配以后對圖像進行縫合，并對縫合的邊界進行平滑處理，使縫合自然過渡。由于任何兩幅相鄰圖像在采集條件上都不可能做到完全相同，因此，對于一些本應該相同的圖像特性(如圖像的光照特性等)，在兩幅圖像中就不會表現得完全一樣。圖像拼接縫隙就是從一幅圖像的圖像區域過渡到另一幅圖像的圖像區域時，由于圖像中的某些相關特性發生了躍變而產生的。圖像融合就是使圖像間的拼接縫隙不明顯，拼接更自然。

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

   CUDA對內存的操作與一般的C程序基本相同，操作顯存則需要調用CUDA API中的存儲器管理函數。一旦確定好程序中的并行部分，就可以將這部分計算交給GPU。運行在GPU上的CUDA并行計算函數稱為Kernel，即內核函數，它并不是一個完整的程序，而是CUDA程序中可以被并行執行的步驟。內核函數必須通過_global_函數類型限定符定義，且只能在主機端代碼中調用。
    CUDA線程結構如圖3所示。Kernel以線程網格(Grid)為組織形式，每個Grid由若干個線程塊(Block)組成，每個Block又由若干個線程(Thread)組成。在程序運行過程中，Kernel是以Block為單位執行的，Grid只是一系列可以被執行的Block的集合。不同Block是并行執行的，沒有執行的先后順序，而且相互無法通信。

    CUDA軟件體系由CUDA Library、CUDA Runtime API和CUDA Driver API構成，如圖 4所示。CUDA的核心是CUDA C語言，需要通過nvcc編譯器進行編譯。編譯得到的僅是GPU端的代碼，要在GPU上分配顯存并啟動Kernel就需要借助CUDA Runtime API或者CUDA Driver API來實現。CUDA Runtime API和CUDA Driver API 提供了實現設備管理、上下文管理、存儲器管理、代碼塊管理和執行控制等操作的應用程序接口。CUDA Runtime API在CUDA Driver API的基礎上進行了封裝，使得編程方便代碼簡潔, 因此通常采用CUDA Runtime API進行項目開發。

 3 SIFT特征匹配算法
    SIFT算法首先在尺度空間進行特征檢測，并確定關鍵點(Key Points)的位置和關鍵點所處的尺度，然后使用關鍵點鄰域梯度的主方向作為該點的方向特征，以實現算子對尺度和方向的無關性。
    SIFT特征匹配算法包括兩個階段：(1)SIFT特征的生成，即從多幅待匹配圖像中提取出對尺度縮放、旋轉和亮度變化無關的特征向量；(2)SIFT特征向量的匹配。
    在實際的尺度不變特征點提取中，SIFT算法將圖像金字塔引入了尺度空間。首先采用不同尺度因子的高斯核對圖像進行卷積以得到圖像的不同尺度空間，將這一組圖像作為金字塔圖像的第一階。接著對其中的2倍尺度圖像(相對于該階第一幅圖像的2倍尺度)以2倍像素距離進行下采樣來得到金字塔圖像第二階的第一幅圖像，對該圖像采用不同尺度因子的高斯核進行卷積，以獲得金字塔圖像第二階的一組圖像。以此類推，獲得高斯金字塔圖像。每一階相鄰的高斯圖像相減，就得到了高斯差分圖像，即DoG圖像。通過擬和三維二次函數以精確確定關鍵點的位置和尺度，同時去除低對比度的關鍵點和不穩定的邊緣響應點(因為DoG算子會產生較強的邊緣響應)，以增強匹配穩定性、提高抗噪聲能力。利用關鍵點鄰域像素的梯度方向分布特性為每個關鍵點指定方向參數，使算子具備旋轉不變性，生成SIFT特征向量。
    接下來以關鍵點為中心取8×8的窗口。圖5(a)的中央黑點為當前關鍵點的位置，每個小格代表關鍵點鄰域所在尺度空間的一個像素，箭頭方向代表該像素的梯度方向，箭頭長度代表梯度模值，圖中圈代表高斯加權的范圍（越靠近關鍵點的像素,其梯度方向信息貢獻越大）。然后在每4×4的小塊上計算8個方向的梯度方向直方圖，繪制每個梯度方向的累加值，即可形成一個種子點，如圖5(b)所示。圖5(b)中,一個關鍵點由2×2共4個種子點組成，每個種子點有8個方向向量信息。這種鄰域方向性信息聯合的思想增強了算法抗噪聲的能力，同時對于含有定位誤差的特征匹配也提供了較好的容錯性。

4 基于CUDA的圖像拼接軟件的設計
4.1 Host端實現
    程序的Host端由C++編寫，負責控制整個程序的執行流程、所有供CPU 和GPU 所用的數據的分配管理以及Device 端模塊的調用。界面使用最基本的Windows SDK編寫。
   在數據初始化階段，包含所有之后處理步驟中所需要的圖像數據對象的生成，將輸入圖像作為高斯金字塔底層，通過系統中的PYRAMID_LEVEL宏指定金字塔的層數，在輸入圖像的尺寸基礎上循環計算各層金字塔圖像的分辨率，并對圖像進行初始化。由于所有的圖像數據需要在設備端處理，使用cudaMallocPitch函數分配數據地址空間，數據結構不再是OpenCV中的IplImage，而是GPU可以識別的uchar數組類型。
    使用cudaMemcpy2D函數將IplImage結構中的原始數據復制到相應高斯金字塔的最底層，也就是uchar數組的第一個元素，供Device端函數使用。隨后進行Kernel函數調用，對于每一個需要處理的金字塔層，Host端發起一次Kernel調用。例如:
    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端實現
   主要Device函數如下：
    (1) reduce（）函數對左、右圖和掩碼圖像各完成一次reduce操作，生成下一層高斯金字塔圖像。reduce變換按照前文所述的方法對目標層的金字塔圖像進行逐像素處理，每一個目標像素的顏色值按一定的權重值對原始圖像中的一個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（）函數對左、右圖像各完成一次expand操作和減法操作，生成下一層拉普拉斯金字塔圖像。expand變換相當于reduce變換的逆過程，它對目標層的金字塔圖像進行逐像素處理，每一個目標像素的顏色值也是按reduce變換中所使用的權重值對原始圖像中的一個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( )函數根據掩碼圖像的高斯金字塔以及左、右圖像的拉普拉斯金字塔合成當前層的目標圖像的拉普拉斯金字塔,所有像素值均以掩碼圖像的高斯金字塔為權重而求得。
    _global_ void blend(uchar* mGaussianDataSrc, uchar* lLaplacianDataSrc, uchar* rLaplacianDataSrc,uchar* sLaplacianDataDst, size_t stride,int lapWidth, int lapHeight)
    (4) collapse（）函數對圖像的拉普拉斯金字塔分別完成一次expand操作和累加操作,本質上等同于expand操作，兩者的基本算法是相同的。不同點在于expand模塊用于各層高斯金字塔的expand操作，從而生成各層拉普拉斯金字塔，而collapse函數則用于整個融合過程最后的圖像重構步驟，將各層已經求得的拉普拉斯金字塔作擴展和累加操作，生成最后的拼接圖像。
    _global_ void collapse(uchar* sLaplacianDataSrc, uchar* sExpandDataSrc, uchar* sExpandDataDst, size_t strideSrc, size_t strideDst, int srcWidth, int srcHeight, int expWidth, int expHeight)
    本文借助于SIFT特征對于旋轉和尺度的不變性以及對于噪聲干擾良好的魯棒性進行圖像拼接與匹配，使用CUDA技術簡單地對多分辨率融合算法進行了優化，提高了其執行效率和速度。編寫了界面化的Demo程序，實現了基本的圖像拼接功能。
參考文獻
[1] 譚康.圖像拼接技術與實現[D].南京:南京理工大學, 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]. 電子與信息學報，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]. 機械科學與技術,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]. 北京：清華大學出版社，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.