GPU加速電子剪切散斑干涉圖像處理

邵 珩, 周 勇, 祁俊峰, 聶中原

(北京衛星制造廠有限公司, 北京 100094)

1 引 言

電子剪切散斑干涉測量技術是一種全場、非接觸、高精度和高靈敏度的光學測量方法,測量系統光路簡單，調節方便，對環境要求不高，特別適合于無損檢測實時測量[1-3]。散斑圖像的處理需要大量的運算。例如采用時間相移法時，首先需要將若干散斑圖進行相減、相除和反正切計算，得到散斑干涉相位差圖，然后通過濾波處理去掉噪聲，獲得散斑條紋圖，最后通過解包裹算法得到物體變形量的空間導數。隨著激光散斑干涉測量技術的應用，相關濾波算法也得到較多研究[4-8]，但這些研究一般都關注于改善濾波效果，對濾波時間基本不做考慮，無法應用于電子剪切散斑干涉連續測量的實時處理。在電子剪切散斑干涉連續測量中，需要計算機在短時間內對硬件獲得的大量圖像數據進行處理，因此計算機的運算能力成為能否實時處理測量數據的關鍵。傳統的CPU計算靈活性高，適合處理各種復雜情況，但CPU單核運算能力有限，隨著圖像處理規模的不斷擴大，CPU的處理能力越來越成為測量結果實時處理顯示的瓶頸；若采用多核CPU計算，一則計算任務的分割、多核CPU間的通信處理較為繁瑣，二則隨著核心數量的增多，CPU平臺的成本也迅速提高，不利于無損檢測實時測量的商業應用。

GPU傳統上用于游戲圖形顯示，近年來因其運算能力、存儲帶寬遠超CPU而被用于加速計算。GPU計算核心數量眾多而單核處理能力較弱，尤其適用于大規模并行計算任務，被廣泛應用于流體動力學模擬[9-10]，分子動力學模擬[11-12]，神經網絡與深度學習算法等[13-15]。由于GPU具有強大的計算能力，易于實現數據的實時處理[16-17]，因此，將GPU用于圖像處理有望大幅提升計算機的運算能力，實現電子剪切散斑干涉的實時測量；同時，GPU的運算能力和能源效率遠超同價位的CPU，有利于相關商業平臺的開發。本研究擬在電子剪切散斑干涉連續測量的數據處理中采用GPU加速計算技術，以提高圖像處理效率，提高測量效率。

2 電子剪切散斑干涉原理

激光照射到粗糙物體表面時發生漫反射，生成無數子波。這些子波在空間中發生干涉，形成明暗不一的斑點，稱為激光散斑。由于激光散斑源于激光在粗糙表面上發生的漫反射，攜帶了粗糙表面的形貌信息，當粗糙表面的形貌發生變化時，激光散斑的分布也隨之發生變化，因而可以利用激光散斑的這一特性對物體表面的形變進行測量。

圖1 基于邁克爾遜剪切裝置的剪切散斑干涉光路圖Fig.1 Shearing speckle interferometry optical system based on Michelson shearing device

圖1為剪切干涉光路示意圖。相干光經粗糙被測物S漫反射，部分光線被分光鏡反射后經反射鏡M1再次反射，穿過分光鏡后進入CCD成實像S1；另一部分光線穿過分光鏡，被反射鏡M2反射后經分光鏡反射進入CCD成實像S2。如果M1和M2相互垂直并且到分光鏡的距離相等時，S1與S2大小相等且互相重合；當反射鏡M2的角度發生偏轉時，S2的位置將發生移動，與S1形成剪切。將S1與S2之間距離的大小和方向映射到被測物上即得到剪切量和剪切方向。此時，S1和S2將在像面上互相干涉，生成散斑圖像。若記錄被測物發生形變前后的散斑干涉圖像，經過處理可以得到表示物體沿剪切方向的位移偏導數的條紋圖案。

3 GPU圖像處理方法

GPU加速計算將計算任務分派給大量的線程同時處理，利用遠高于CPU的線程數量獲得遠超出CPU的性能。GPU本身單核計算能力很弱，因此只有將計算任務分割為互相獨立的小任務，使每個線程都可以不依賴其他線程的計算結果而完成本身的計算任務時，才能發揮GPU線程數量多的優勢。

采用CUDA平臺處理計算任務時，CPU作為控制器，負責讀取數據，并將數據傳輸給GPU；GPU作為協處理器，接收CPU傳遞的數據，在顯存中開辟相應存儲空間，存儲傳入的數據和中間處理數據，將這些數據分配給GPU上的多個線程進行處理；處理完成后，將數據傳回CPU，由CPU進行后續處理。

電子剪切散斑干涉實時測量過程中，圖像處理包括兩部分：解相位差和濾波。解相位差采用四步相移法[18]：首先在加載前對被測物采集4幅散斑圖，每次采集后對實像S1引入π/2的相移。4幅圖的光強可以表示為：

(1)

式中：Ii為第i步測得的光強分布，i=1,2,3,4，a為背景光強，b為調制度，α為未知的隨機相位。

加載后，再次采集4幅散斑圖，光強為：

(2)

式中：β為加載引入的相位差。

根據Ii和Ii′，可以解出每個點在加載前后位置的相位差：

(3)

(4)

θ=θ2-θ1.

(5)

為濾除相位差圖的噪聲，采用3種速度較快的濾波方法，包括正余弦均值濾波、復數均值濾波、Butterworth低通濾波方法，分別比較它們的精度、計算效率以及GPU加速效果。

(6)

(7)

對正余弦均值濾波方法采用GPU加速時，首先將圖像分割為32×32的小塊，將每個小塊及其外圍一圈像素分配給一個線程塊，即采用3×3窗口時，每個線程塊分配34×34的小塊，其中超出圖像邊界的部分采用鏡像處理。對小塊內像素的相位值分別求正弦和余弦值，然后采用式(6)、(7)求得正余弦均值并通過反正切函數還原為相位。為提高計算效率，每個線程塊內采用共享存儲存取正余弦值。

復數均值濾波中，首先將相位值轉換為復數：

cx,y=exp(θx,yi)，

(8)

然后按照類似式(6)的方式對其求取均值：

(9)

多次濾波完成之后再通過反正切函數將復數平均值還原為相位值。

比較正余弦均值濾波和復數均值濾波可以發現，在進行一次濾波時，兩者是等價的。在進行多次濾波時，復數均值濾波減少了正余弦計算和反正切計算次數，降低了計算量。

對復數均值濾波進行GPU加速時，將相位與復數的轉換與還原、復數求均值分別進行GPU加速計算，其中求均值部分采用類似正余弦濾波圖像分割的方法采用共享存儲加速計算。

Butterworth低通濾波中，對相位差圖任意點(x,y)的相位值θ，將其轉換為復數c。對cx,y進行傅里葉變換，得到Cu,v，然后進行低通濾波：

(10)

式中，Du,v是到原點的距離，D0為截止距離，n為濾波器的階數。

對Butterworth低通濾波采用GPU加速時，可以將轉換與還原、正反傅里葉變換和低通濾波每個步驟均采用相應的GPU計算。

4 實驗結果

4.1 相位差圖生成

4.2.1 CPU解相位差

采用CPU解相位差時，首先讀取8幅圖片的數據，然后掃描圖片，解出每個點的相位差值。解相位差的全部時長約500 ms，其中25 ms用于讀取圖片，465 ms用于解相位差計算，10 ms用于保存圖片，解相位差的時間占總時間的90%以上。采用的CPU為i7-6500U@2.5 GHz，圖片尺寸為2 592×2 048，下同。圖2為計算得到的相位差圖。

圖2 四步相移法得到的相位差圖Fig.2 Phase-difference map obtained by four-step phase shift method

4.2.2 GPU解相位差

從式(3)～(5)可知，每個點的相位差求解過程都與周圍點無關，因此適用于GPU并行處理。本研究中采用兩個不同的GPU進行圖像處理，其中GPU1型號為GeForce 940MX@1004 MHz，GPU2型號為Quadro P3000@1088 MHz，下同。

求解結果表明，采用GPU解相位差的耗時遠少于CPU解相位差，如表1所示。采用GPU1時，全部求解時長約77 ms，其中26 ms用于讀取圖片，43 ms用于解相位差計算，8 ms用于保存圖片，解相位差的時間下降到總時間的約56%；采用GPU2時，解相位差時間下降到18 ms，總耗時52 ms。

表1 GPU加速解相位差與CPU解相位差耗時對比

Tab.1 Time consumption of phase difference calculation by CPU and GPU

ms

進一步分析表明，對于GPU1，在43 ms的解相位差時間中，GPU初始化和分配顯存耗時10 ms，從CPU到GPU的數據傳輸耗時約23 ms，計算耗時約6 ms，數據從GPU傳回至CPU耗時約4 ms；對于GPU2，初始化和分配顯存耗時約7 ms，從CPU到GPU的數據傳輸耗時約7 ms，計算和數據回傳耗時均約2 ms。若僅考慮計算時間，GPU1的計算速度是CPU單核的約80倍，GPU2的速度則達到200倍以上。由此可見，在解相位差計算中，GPU初始化和顯存分配、數據傳輸占據了絕大部分時間，真正用于計算的時間非常短。

計算時間對比表明， CPU解相位差的瓶頸在于求解過程，而采用GPU加速計算時，求解過程耗時占比極小，絕大部分時間用于數據存取與傳輸，其次為GPU的初始化過程，計算總時間降到原來的約10.4%～15.4%，大幅加速了求相位差圖過程。

4.2 相位差圖濾波處理

4.2.1 CPU串行算法濾波

為便于比較不同濾波算法的效果，首先采用模擬散斑圖進行測試。模擬相位分布為：

d1=sqrt[(x-1696)2+(y-1024)2]，

(11)

d2=sqrt[(x-896)2+(y-1024)2]，

(12)

θx,y=3000×[exp(d1/1000)-exp(d2/1000)].

(13)

然后給模擬相位附加均值為0、標準差為60的高斯噪聲，得到含噪聲散斑圖。無噪聲相位圖和含噪聲散斑圖如圖3所示。

圖3 模擬相位圖。(a)無噪聲； (b)含噪聲。Fig.3 Simulated phase diagram. (a)Without noise； (b)With noise.

對圖3(b)采用不同的算法進行濾波，并比較其標準峰值信噪比 PSNR的變化。其中PSNR的定義為：

(14)

式中:KMSE為濾波后圖像與不含噪聲的原始相位圖之間的均方誤差。

由定義可知，PSNR越大，濾波后圖像與原始相位圖差別越小，濾波失真越小。圖4為正余弦均值濾波和復數均值濾波的PSNR隨濾波次數和濾波窗口尺寸變化的曲線及不同階次Butterworth低通濾波的PSNR隨截止距離的變化。

圖4 正余弦均值濾波(SC)和復數均值濾波(Mean)的PSNR隨濾波次數和濾波窗口大小變化(a)， 與不同階次Butterworth低通濾波的PSNR隨截止距離的變化(b)。Fig.4 Influence of the filtering times and window size on PSNR of sin-cos filtering (SC) and complex mean filtering (Mean) (a)， and influence of the order and cut-off distance on PSNR of Butterworth filtering (b).

圖4表明，對于正余弦均值濾波和復數均值濾波而言，濾波窗口越大，PSNR上升越快，且在相同的濾波窗口下正余弦濾波效果較復數濾波更好。對于Butterworth低通濾波，階次越低，達到最佳效果的截止距離越短?？傮w而言，正余弦均值濾波和復數均值濾波的最好效果好于Butterworth低通濾波。由于Butterworth濾波的計算時間比較穩定，不隨階次和截止距離而變化，以其達到的最佳效果為標準(n=2、D0=37時，PSNR為88.84)，選取PSNR≥88.84的正余弦均值濾波和復數均值濾波不同濾波窗口下最短計算時間(“3×3”窗口下濾波100次以內PSNR值均未達88.84，因此不再考慮)，結果見表2。

表2表明，對于正余弦均值濾波，每次濾波時重新計算正余弦值和相位值耗時較大，因此濾波次數對濾波時間的影響比濾波窗口影響更大；而對于復數均值濾波，濾波次數與濾波窗口對總計算時間的影響相當，因此達到相同效果時采用不同濾波窗口的計算時間較為接近。盡管達到同樣效果時復數均值濾波所需次數大于正余弦均值濾波，但由于不必每次濾波都重新計算正余弦值和相位值，其計算速度顯著快于正余弦均值濾波。

表2 達到Butterworth濾波最佳效果時正余弦均值濾波和復數均值濾波不同濾波窗口下最短計算時間

Tab.2 The shortest times for sin-cos filtering and complex mean filtering with different filtering window sizes to reach the best of Butterworth filtering

算法PSNR計算參數耗時/msButterworth濾波88.84n=2,D0=371416正余弦均值濾波88.92“7×7”,24次11995正余弦均值濾波89.00“5×5”,39次15427復數均值濾波89.08“7×7”,26次7023復數均值濾波88.88“5×5”,41次6719

圖5為對模擬圖進行不同方法濾波所得結果。從圖中可以看出，Butterworth濾波對細節的平滑最為徹底，條紋最為平滑，剩余微小起伏最小，但其缺點在于邊緣效應明顯，在圖像左右邊緣出現了平行于邊緣的多余條紋，甚至產生了殘差點。正余弦均值與復數均值濾波結果基本相同，條紋較不平滑，殘余微小起伏較大，邊緣效應不明顯，濾波圖中無殘差點。

圖5 模擬圖濾波結果。(a)Butterworth濾波，n=2,D0=37;(b)正余弦均值濾波“7×7”,24次；(c)復數均值濾波，“7×7”，26次。Fig.5 Filtered images of simulated diagram with different filtering methods. (a)Butterworth filtering, n=2, D0=37； (b) Sin-cos filtering, “7×7”, 24 times； (c) Complex mean filtering, “7×7”, 26 times.

圖5表明，Butterworth濾波對圖像內部的濾波效果較好，但邊緣部分結果較差；正余弦均值與復數均值濾波對圖像內部濾波效果略差，但邊緣部分結果較為穩定，濾波效果與圖像內部相當。

總體而言，正余弦均值濾波和復數均值濾波達成與之相當的效果所需的計算時間遠遠超出Butterworth濾波。雖然可以通過減少濾波次數縮短計算耗時，但將時間縮短至Butterworth濾波耗時將嚴重惡化濾波效果。因此，在CPU平臺上，Butterworth濾波效率遠高于正余弦均值濾波和復數均值濾波。

4.2.2 GPU并行算法濾波

由于GPU濾波加速算法基于CPU濾波算法，兩者結果理論上一致，因此仍按表2中的計算參數對各種算法重新測定計算時間。結果如表3所示。

表3表明，雖然仍然是Butterworth濾波速度最快，但其他算法的速度劣勢已經大幅減小，可以犧牲一定的精度換取更短的計算時間。例如在GPU1上，采用“7×7”窗口進行濾波時，若將濾波次數減少到15次，正余弦均值濾波和復數均值濾波需要的時間將縮短為145 ms和161 ms，時間短于Butterworth濾波，相應PSNR值為85.39和85.15，與Butterworth濾波精度最好時相差不大。就GPU加速而言，Butterworth濾波加速比最低，僅有正余弦均值濾波的1/9～1/6，復數均值濾波的1/3左右。究其原因，主要是因為Butterworth濾波需要進行正反快速傅里葉變換，計算中內存讀寫不連續，控制較為復雜，而正余弦均值濾波和復數均值濾波計算中內存讀寫連續，有利于提高內存讀寫效率。GPU加速正余弦均值濾波快于復數均值濾波表明，在這兩種GPU算法中，主要瓶頸已經由不再是正余弦和反正切計算，而是內存的讀寫，包括線程塊內共享內存的讀寫和線程塊與全局顯存間的讀寫。復數均值濾波的全局內存讀寫量為正余弦均值濾波的2倍，因此在消除正余弦和反正切重復計算之后速度仍然慢于正余弦均值濾波。

表3 采用GPU加速計算時，達到Butterworth濾波最佳效果時正余弦均值濾波和復數均值濾波不同濾波窗口下最短計算時間

Tab.3 The shortest times for sin-cos filtering and complex mean filtering with different filtering window sizes to reach the best of Butterworth filtering adopting GPU accelerated computing

算法計算參數耗時/msGPU1GPU2加速比GPU1GPU2Butterworth濾波n=2,D0=37170508.328.3正余弦均值濾波“7×7”,24次2067158169正余弦均值濾波“5×5”,39次2157272214復數均值濾波“7×7”,26次2468028.688復數均值濾波“5×5”,41次2968222.782

4.2.3 解相位差與濾波合并計算

由于濾波讀取的數據源自解相位差得到的數據，將相位差數據從GPU傳到CPU，保存至硬盤，再從硬盤中讀取數據，傳到GPU進行濾波，可以看出中間傳遞數據的部分實際上做了無用功。因此，將解相位差與濾波合并計算，即在GPU中解相位差后直接進行濾波可進一步提高效率。表4為合并解相位差和濾波計算后的總計算時間，其較解相位差全過程時間和單純濾波時間之和減少約10～20 ms，這主要是因為去除了冗余的數據傳送過程。

因此，采用GPU加速計算后，求解相位差和濾波的總時間可以從CPU計算時的1 900 ms降低到239 ms或86 ms。采用性能較弱的GPU時，可以通過降低濾波次數將正余弦均值濾波時間進一步縮短，從而保證濾波的實時性，如正余弦均值濾波采用“7×7”濾波10次時，總計算時間為190 ms，代價是濾波不充分，PSNR下降至81.5；采用性能較強的GPU時，可以直接采用Butterworth低通濾波，同時獲得較好的濾波效果和濾波實時性。

表4 解相位差合并不同濾波算法的計算總耗時Tab.4 Total procssing time of difforent filting method combined with phase calculation

圖6為對原始圖進行解相位差和濾波所得結果。與模擬圖最大的不同在于，因為激光照明質量和遮擋等原因，實際散斑圖各處相干程度相差很大。從圖2(相位差圖)可以看出，圖像中心部位條紋明顯，相干程度強，外圍條紋暗弱，相干程度弱，圖像四角和左側則無條紋，為非相干區。圖6表明，Butterworth濾波對細節的平滑最為徹底，不僅外圍弱相干和非相干區的殘差點最少，中心強相干區的微小起伏也最小，但其缺點在于邊緣效應明顯，在圖像邊緣平行于邊緣的多余條紋。正余弦均值濾波在弱相干區留下的殘差點最多，強相干區的微小起伏也最大。復數均值濾波介于兩者之間，對弱相干區殘差點的消除更接近Butterworth濾波，強相干區的微小起伏程度更接近正余弦均值濾波，且和正余弦均值濾波一樣邊緣效應不明顯。

圖6 解相位差和濾波合并計算結果。(a)Butterworth濾波，n=2,D0=37;(b)正余弦均值濾波“7×7”,24次;(c)復數均值濾波，“7×7”，26次。Fig.6 Results of combined phase difference and filtering calculation with different filtering methods:(a)Butterworth filtering, n=2, D0=37; (b) Sin-cos filtering, “7×7”, 24 times; (c) Complex mean filtering, “7×7”, 26 times.

5 結 論

(1)采用GPU計算可以大幅提高運算速度。相對于i7-6500U單線程，GeForce 940MX的濾波計算速度可以提高8.3倍(Butterworth低通濾波)、58倍(正余弦均值濾波，“7×7” 窗口)和28.6倍(復數均值濾波，“7×7”窗口)，Quadro P3000可提高28.3倍(Butterworth低通濾波)、169倍(正余弦均值濾波，“7×7” 窗口)和88倍(復數均值濾波，“7×7” 窗口)。

(2)不同濾波方法有不同的優缺點，Butterworth低通濾波對強相干區和弱相干區濾波效果均較好，濾波結果平滑，但在圖像邊緣有嚴重邊緣效應，易出現異常條紋；正余弦濾波平滑程度低于Butterworth低通濾波，在弱相干區容易殘留大量殘差點，而邊緣效應不明顯；復數均值濾波平滑程度與正余弦濾波相似，在弱相干區效果與Butterworth低通濾波相似，邊緣效應不明顯。

(3)采用GPU計算時，速度瓶頸很大程度上出現在數據傳輸中，包括系統內存與GPU顯存之間的傳輸和GPU核心與顯存間的傳輸。數據傳輸量越少、顯存讀寫連續性越好，越有利于發揮GPU計算速度快的優點。

(4)對于性能不同的GPU，可以采取不同的算法在濾波速度和濾波精度間取得平衡。對性能較低的GPU，可以適當降低精度要求，采用正余弦均值濾波通過減少濾波次數的方式提高濾波速度，保證實時性；對于性能較強的GPU，可以采用Butterworth低通濾波同時獲得較好的濾波效果和濾波實時性；如果GPU性能更強，可以采用正余弦均值濾波增加濾波次數以獲得更好的濾波效果。