宋 驥,周松濤

(武漢大學(xué) 國際軟件學(xué)院，湖北 武漢 430079)

圖像匹配在數(shù)字?jǐn)z影測量、計算機視覺、遙感影像處理、醫(yī)學(xué)圖像、文件夾和數(shù)據(jù)庫的圖像搜索等領(lǐng)域有廣泛應(yīng)用.在現(xiàn)在文字檢索已然成熟的時代，以“圖”搜“圖”已經(jīng)成為下一步的需求.在軍事領(lǐng)域搜尋目標(biāo)時，甚至需要一張圖中迅速找到眾多目標(biāo)的精確位置.

目前用于圖像匹配的方法主要有基于圖像灰度相關(guān)方法、基于圖像特征方法、基于神經(jīng)網(wǎng)絡(luò)和遺傳算法等人工智能方法.基于圖像灰度的匹配算法簡單，匹配準(zhǔn)確率高，適合硬件實現(xiàn)，但是計算量大，不利于實時處理，對灰度變化、旋轉(zhuǎn)、形變以及遮擋等比較敏感；基于圖像特征的方法計算量相對較小，對灰度變化、形變及遮擋有較好的適應(yīng)性，但是匹配精度不高；基于神經(jīng)網(wǎng)絡(luò)、遺傳算法等人工智能的方法發(fā)展較晚，算法還不成熟[1].最近幾年在GPU基礎(chǔ)上發(fā)展起來的CUDA并行計算技術(shù)，在高性能計算方面得到了廣泛應(yīng)用，在圖像匹配上也用于圖像灰度相關(guān)算法的加速[2]，雖得到了一定的加速比，但是目前這種算法設(shè)計的效率還不夠高，當(dāng)數(shù)據(jù)規(guī)模增大時對問題的適應(yīng)性下降.本文基于圖像灰度相關(guān)算法，提出一種更加優(yōu)秀、高效的CUDA并行算法，使得在原來的基礎(chǔ)上進一步加速，并且對問題規(guī)模的繼續(xù)加大有很好的支持和適應(yīng).

圖1 圖像匹配過程Fig.1 Process of image matching

圖2 相鄰矩陣重復(fù)情況Fig.2 The repeat in neighborhood matrices

1 問題描述

假設(shè)需要在待搜索圖像中匹配給定目標(biāo)圖像.一般來說，目標(biāo)圖像較小(設(shè)為m×n)，待搜索圖像尺寸較大(設(shè)為(M+m)×(N+n)).圖像規(guī)格是0-255的8位灰度圖像，如果是彩色圖像，則先灰度化.

如圖1，在待搜索圖像中枚舉矩陣左上角(x,y)，逐一檢查待匹配的m×n矩陣，將其與目標(biāo)圖像進行匹配運算，若相關(guān)系數(shù)大于設(shè)定的閾值，則認(rèn)為該矩陣匹配成功.

記小圖像為g，大圖像為S，用Sx, y表示S中以(x,y)為左上角點與g大小相同的子塊.將兩個相同尺寸的矩陣的相關(guān)系數(shù)記為ρ(x,y).

將協(xié)方差和方差公式代入并化簡[1]可得公式1：

(1)

2 算法分析與CPU優(yōu)化

通過優(yōu)化，只需要一步運算(O(1)的時間)便可以完成相鄰像素點的遞推，式(1)中求平均值也可一步運算解決.如果搜索整個圖像，則式(1)的分母部分，以及分子中平均值相乘部分，基于CPU的串行計算的整體算法復(fù)雜度在O(M×N)，這已經(jīng)完全可以被CPU勝任，不再成為算法的瓶頸.

而式(1)的分子中，S與g相乘并求和的部分的算法復(fù)雜度依然是O(M×N×m×n)，這將導(dǎo)致該算法在CPU上的計算很慢，該部分已經(jīng)成為整個算法繼續(xù)加速和優(yōu)化的瓶頸所在.

為描述簡便，將Sx, y與g相乘并求和的結(jié)果記為Sum(x,y)，用來表示以(x,y)為左上角的m×n大小的像素矩陣與目標(biāo)矩陣對應(yīng)像素乘積之和.Sum(x,y)的計算成為問題的關(guān)鍵，因此下面的GPU算法主要針對Sum(x,y)進行計算.

(2)

3 利用GPU的初步并行算法

現(xiàn)代的顯示芯片已經(jīng)具有高度的可程序化能力、相當(dāng)高的內(nèi)存帶寬、以及大量的執(zhí)行單元，因此可利用顯示芯片來幫助進行一些計算工作，即 GPGPU[3].雖然GPU并沒有CPU所具備的強大的復(fù)雜邏輯指令的處理能力，卻非常適合勝任大量、簡單、重復(fù)的運算工作，以及高度并行化的并行計算.圖形硬件有強大的并行性，支撐高密度的運算，具有很高的浮點運算速度，因此采用圖形硬件主要是為了加速運算.

根據(jù)GPU的特點可以看出，式(2)是一個簡單、重復(fù)而大量的運算工作，十分適合并行處理.基于CUDA用GPU的并行計算來解決圖像匹配問題，已經(jīng)得到了初步的應(yīng)用[2].

3.1 CUDA模型簡介

CUDA是NVIDIA的GPGPU模型，它使用C語言為基礎(chǔ).在CUDA架構(gòu)下，顯示芯片執(zhí)行時的最小單位是thread，數(shù)個thread可以組成一個block[3].

從CUDA程序邏輯上看，采用兩層并行，每個block內(nèi)的thread細粒度并行處理，并且用共享寄存器進行通信，而各個block之間粗粒度并行，并行度視硬件情況而定.兩層并行的好處之一，在于實現(xiàn)了可擴展的程序模型.當(dāng)硬件的運算內(nèi)核數(shù)量增加時，CUDA可以自適應(yīng)的方式，自動增大粗粒度并行度(即同時運算的block數(shù)量)[4].好處之二在于細粒度并行中，不同線程(thread)之間可以通過共享寄存器進行通信，達到優(yōu)化數(shù)據(jù)傳輸速度和加速計算的目的.

圖3 算法2的分塊方式及一個block內(nèi)部線程分配Fig.3 Designs of blocks and threads in algorithm 2

圖4 使用共享存儲器的初步設(shè)想Fig.4 Abstract design of shared memory

3.2 初步的并行策略

通過對CUDA的并行機制的分析，做出以下分塊策略[2].首先，將一個Sum(x,y)的計算作為一個block，并使block(x,y)計算Sum(x,y)，此時共有M×N個block，根據(jù)硬件條件若干個block并行計算，直到M×N個Sum(x,y)計算完畢.而后將每一block分配m×n個線程，每個線程在S和g中各取一個像素點的灰度并將其相乘累加，計算一個Sum(x,y).

初步的并行方法的速度已經(jīng)較CPU有很大的提升，但是該并行方法存在以下諸多缺點：

1)數(shù)據(jù)重復(fù)調(diào)度：在GPU中計算，必須先把內(nèi)存中的數(shù)據(jù)拷貝到顯存中，然后訪問顯存進行計算，而訪問顯存需要消耗大量的時鐘周期，所以設(shè)計時應(yīng)該盡量減少對顯存的訪問次數(shù).CUDA在顯卡端(設(shè)備端)提供了一些可以高速緩存的寄存器，包括register和共享寄存器(shared memory)，可以將顯存中的數(shù)據(jù)放到共享寄存器中，以此減少數(shù)據(jù)傳輸時間.

在上述分塊方法中，每一個像素對應(yīng)的block之間粗粒度并行，并不能用共享寄存器共享數(shù)據(jù)，所以此時每一個block都必須把m×n的矩陣數(shù)據(jù)從顯存中讀取出來.然而相鄰像素之間(也即相鄰block之間)訪問的顯存數(shù)據(jù)大量重復(fù)(目標(biāo)矩陣重復(fù)m×n，即完全重復(fù)，待搜索矩陣重復(fù)(m-1)×n).

2)數(shù)據(jù)尺寸適應(yīng)性不強：在目前的主流GPU中，每一個block所能同時管理的線程數(shù)量是有上限的(512)，當(dāng)目標(biāo)矩陣的尺寸m×n大于線程上限時，線程數(shù)將不夠，無法使得每一個像素對應(yīng)一個線程.此時必須修改，讓每一個線程計算多個像素.而同時，當(dāng)目標(biāo)矩陣的尺寸變大時，在缺點1)中所闡述的數(shù)據(jù)的重復(fù)量將變大，目標(biāo)矩陣尺寸越大，數(shù)據(jù)重復(fù)問題越嚴(yán)重.在這樣的情況下，該算法對于數(shù)據(jù)規(guī)模的適應(yīng)性大打折扣.

4 高效的GPU并行算法

4.1 改變分塊方式

針對上述缺點對分塊方式進行改進，將連續(xù)的BLOCKSIZE*BLOCKSIZE(16×16)數(shù)量個Sum(x,y)的計算放到一個block中(形成一個方陣)，同時使一個block包含16×16個線程(參見圖3)，從而使每一個線程負(fù)責(zé)計算一個Sum(x,y).此時整個待搜索矩陣分塊為(M/16)×(N/16)個block，每個block計算16*16個Sum(x,y)(無法整除情況可填補像素).

這樣的好處在于可以同時解決數(shù)據(jù)重復(fù)調(diào)度和尺寸適應(yīng)性問題.

首先，連續(xù)的16×16個線程處于同一個block中，他們可以通過共享寄存器共享數(shù)據(jù).在原來的方法中，16×16個Sum(x,y)的計算需要調(diào)度顯存共(m×n×2)×(16×16)次，而通過共享數(shù)據(jù)的方式，在所有線程開始階段，可以將他們共同使用的待搜索圖像的(16+m)×(16+n)灰度數(shù)據(jù)以及目標(biāo)圖像(m×n)讀入，柵欄同步(barrier synchronization)后，各個線程開始計算自己負(fù)責(zé)的m×n的矩陣，算出一個Sum(x,y)，但是計算過程中可以共享地訪問由前面讀入到共享寄存器內(nèi)的數(shù)據(jù).數(shù)據(jù)調(diào)度效率大約是原來的16×16倍.(參考圖4)

其次，每一個線程處理一個m×n的矩陣的計算，得到一個Sum(x,y)，當(dāng)m×n增大時，不會受到線程數(shù)量限制的影響.

4.2 共享寄存器的處理方法

雖然理論上可以通過共享寄存器一次性讀入數(shù)據(jù)，但是由于目前硬件限制，每一個SM中的共享寄存器空間被限制為一個上限(16 KB)，所以不能一次性讀入和計算.該限制可以用人工分塊的方法(blocking)解決.通過分塊的方法，分塊、分次讀入，同時也分塊、分次計算.

假設(shè)目標(biāo)圖像的尺寸m×n= 48×48，則將其分解成16×16的小塊，一共3×3=9個小塊.在一個block內(nèi)同步的所有線程中，枚舉這9個小塊，處理完所有的小塊，則累加計算完成，Sum(x,y)也就計算完畢了.處理一個小塊的方法和1中所述的類似，對于一個16×16的小塊，先將目標(biāo)矩陣g對應(yīng)的16×16個像素灰度值讀入共享寄存器，而對于待搜索圖像的數(shù)據(jù)此時則需要讀入32×32個像素的數(shù)據(jù)(如圖(5)，對于16×16個線程，Sum(x,y)對應(yīng)的16×16的小塊是從(x+i,y+j)到(x+i+15,y+j+15)，Sum(x+15,y+15)對應(yīng)的16×16的小塊是從(x+i+15,y+j+15)到(x+i+30,y+j+30)，其他Sum所用數(shù)據(jù)均在這31×31之中，但是為了滿足合并訪問，此處讀入32×32).柵欄同步(barrier synchronization)后，每個線程并行開始各自的16×16次的相乘計算，并且累加到用于加和的變量中(該變量是每個線程私有的，定義為register寄存器類型).計算完畢后再次柵欄同步，而后進入下一個16×16的小塊的計算，重復(fù)上述過程(參見圖5).

在改進算法的人工分塊(blocking)方法中，相鄰共享寄存器覆蓋的灰度數(shù)據(jù)還會出現(xiàn)重復(fù)的情況(如圖6).所以在共享寄存器的讀取數(shù)據(jù)時，可以采取將上一次共享寄存器的右邊16×32個灰度讀出來，直接移動到這一次共享寄存器的左邊16×32的位置中，而不用重新調(diào)度顯存中的這16×32個數(shù)據(jù)，只需要調(diào)度右邊16×32個位置(首次從顯存中讀入)對應(yīng)的灰度(如圖6).這樣進一步提高共享寄存器的利用率，也進一步減少對顯存的重復(fù)訪問.

圖5 人工分塊中共享存儲器的處理流程與方法 圖6 重復(fù)訪問的共享存儲器的處理Fig.5 The method of blocking using shared memory Fig.6 The method to solve repeat problem

4.3 改進后的并行算法的優(yōu)點

1)充分利用共享寄存器的作用，降低顯存訪問次數(shù)，速度大幅提升.

2)不會受到目標(biāo)圖像尺寸的影響，目標(biāo)圖像尺寸增大時，不受限制，適應(yīng)性好.

3)不會降低粗粒度層面(block)的并行度，雖然block數(shù)量下降(下降16×16倍，但是由于待搜索圖像的尺寸巨大和海量的特點，block數(shù)量依然足夠，而且不妨礙多GPU并行.

4)優(yōu)化計算指令，使用CUDA內(nèi)置整數(shù)乘法函數(shù)，達到最高運算速度.

5 實驗

實驗所使用的硬件：CPU為AMD Turion 64×2，主頻2.2GHz.顯卡GeForce 8600M GS，多處理器數(shù)目4，CUDA核心數(shù)目32，顯存512M[4].

實驗數(shù)據(jù)如表1所示，記初步GPU并行算法為算法1；改進后的高效算法為算法2；用時均指Kernel時間；加速比均相較CPU用時.

CPU優(yōu)化中，已經(jīng)將除開式(2)的其他部分做了優(yōu)化，該部分用時在實驗數(shù)據(jù)中均在10ms左右，已經(jīng)不對最后算法時間造成過大影響，所以不參與比較.該表格中CPU用時是指用CPU的方法運算式(2)的所用時間.

表1 三種算法的匹配時間比較

實驗數(shù)據(jù)分析：

初步的GPU并行計算速度大概是CPU的5倍左右，數(shù)據(jù)尺寸越大，由于并行度提高，加速比有所提升.在數(shù)據(jù)較小時，計算量不大，算法1的重復(fù)訪問內(nèi)存的影響凸顯，在數(shù)據(jù)增大時，數(shù)據(jù)延遲的不足雖然被大計算量和大并行度稍微掩蓋，但是依然有不少的影響.

由于初步的GPU算法受限于傳輸數(shù)據(jù)，所以沒有達到最佳效率.以第1個數(shù)據(jù)為例，根據(jù)硬件條件和計算公式[5]，理論帶寬為：Theoretical Bandwidth=(內(nèi)存時鐘×106×(內(nèi)存接口寬度/8) ×2)/109=12.8GB/s(GeForce 8600M GS內(nèi)存時鐘400 MHz，內(nèi)存接口寬度128 bit).實際寬度為：Effective bandwidth = ((Br+Bw)/109)/time經(jīng)過計算，算法1的實際帶寬為：(16×16×512×512×4)/10^9/time=0.67 GB/s，遠遠小于理論帶寬，需要改進.而用改進后的高效并行算法2的實際帶寬為：(16×16×512×512×4)/10^9/time=9.2GB/s，已經(jīng)幾乎達到理論帶寬，此時改進后的算法的瓶頸已經(jīng)不會出現(xiàn)在數(shù)據(jù)傳輸上，而是計算指令上.

從數(shù)據(jù)可以看出，在目前主流軟硬件條件下，改進的算法是原有算法速度的7～12倍，大大提高了效率.

6 結(jié)語

本文提出一種高效GPU并行算法，避免了原有算法過多使用低速存儲器所帶來的問題，利用高速共享寄存器大幅提高效率，進一步加快圖像匹配問題的求解速度.

現(xiàn)今，GPU并行計算的高速性為各個領(lǐng)域解決復(fù)雜而海量的計算問題提供了一種通用方法.這種方法門檻低、見效快、效果好.但是，在設(shè)計CUDA并行計算算法時，應(yīng)綜合考慮計算指令優(yōu)化和數(shù)據(jù)調(diào)度傳輸優(yōu)化，充分利用粗粒度并行帶來的擴展性和細粒度并行中的數(shù)據(jù)共享通訊特性，將數(shù)據(jù)調(diào)度傳輸問題解決，不讓其成為瓶頸.

現(xiàn)在若需要將問題進一步加速，關(guān)鍵則在于以下三個方面：新型的問題求解方法和計算公式，計算機硬件本身的發(fā)展和加速，以及運用多設(shè)備和設(shè)備集群等大規(guī)模聯(lián)合計算，甚至是云計算.

[1] 李卓，邱慧娟.基于相關(guān)系數(shù)的快速圖像匹配研究[J].北京理工大學(xué)學(xué)報，2007，27 (11)：998-1000.

[2] 肖漢，張祖勛.基于GPGPU的并行影像匹配算法[J].測繪學(xué)報，2010，39(1):46-51.

[3] Hotball.深入淺出談CUDA[OL].[2008-06-04].http://www.kimicat.com/cuda簡介.

[4] NVIDIA.CUDA C Programming Guide 3.2[OL].[2010-10-22].http://developer.download.nvidia.com/compute/cuda/3_2_prod/toolkit/docs/CUDA_C_Programming_Guide.pdf.

[5] NVIDIA.CUDA C Best Practices Guide 3.2[OL].[2010-8-20].http://developer.download.nvidia.com/compute/cuda/3_2_prod/toolkit/docs/CUDA_C_Best_Practices_Guide.pdf.