在醫學圖像的處理及分析中, 常常需要對圖像進行插值運算。盡管高斯徑向基(GRBF)插值有插值精度高的優點, 但運算時間長的不足仍限制了它在圖像插值中的應用。因此, 本文提出采用基于計算統一設備架構(CUDA)的方法實現二維和三維醫學圖像的GRBF快速插值。根據CUDA單指令多線程(SIMT)的執行模型, 采用合并訪存、共享內存等各種合適的內存優化措施。并且在應用對數據空間進行二維分塊, 三維分體策略的過程中使用基于重疊區域的自然縫合算法來消除圖像插值連接邊界的失真現象。在保持較高圖像插值精度的基礎上, 二維和三維醫學圖像GRBF插值各基本計算步驟都得到了極大的加速。實驗結果表明:基于CUDA平臺的GRBF插值執行效率與傳統CPU運算相比明顯提高, 對其在圖像插值中的應用具有相當的參考價值。
引用本文： 陳浩, 陳兆學, 喻海中. 基于計算統一設備架構的高斯徑向基圖像插值快速實現方法研究. 生物醫學工程學雜志, 2014, 31(2): 237-244. doi: 10.7507/1001-5515.20140045 復制
引言
在醫學圖像處理和分析過程中，往往需要基于原始的具有較低分辨率的圖像數據得到具有更高分辨率的圖像數據。如對于二維醫學圖像而言，人們常常需要將原始圖像放大，從而更好地觀察細胞或者腫瘤的形狀、位置以及和周圍物質的相對關系等等。對于三維圖像序列，由于成像設備本身精度、成像時間以及存儲器容量的限制，采集過程中獲得的任意兩切片間的空間距離均大于片內各像素之間的距離，即存在不同方向上分辨率的差異。當用所獲得的圖像序列來重建三維結構時，可能會使重建結果產生比例變形。這往往都需要通過插值過程予以修正。因此，圖像插值在二維醫學圖像縮放和三維圖像重建過程中是一種常用的方法。一般說來，常用的醫學圖像插值方法可以分為基于圖像的插值方法、基于對象的插值方法和基于小波的插值方法[1]。基于圖像的插值方法包括最近鄰插值、線性插值、匹配插值[2]等；基于對象的插值方法主要是采用形態學中的膨脹與腐蝕；基于小波的插值方法是通過小波變換獲得相應小波系數后對斷層圖像進行強度和位置插值[3]。本文基于高斯徑向基(Gaussian radial basis function, GRBF)插值的方法是一種基于圖像的插值方法。目前較多地應用于曲線和曲面的插值逼近場合，繆報通等[4]基于徑向基神經網絡也將GRBF用于散亂數據插值； Fornberg等[5]采用RBF-QR方法消除形狀參數對GRBF插值穩定性計算的影響，Liao等[6]擴展了有理GRBF，結合形狀參數實現了曲線和曲面的插值。如果說最近鄰插值、線性插值、樣條插值是一種“局部”插值，即插值點只與周圍少量數據點相關，那么GRBF插值則是一種相對“全局”的插值，高精度的優點使GRBF插值圖像不存在邊緣鋸齒現象。但基于GRBF的插值過程因其數據計算量龐大，需要的運算時間比較長，這限制了其在圖像縮放和重建過程中的應用。
隨著圖形處理單元(GPU)并行運算能力的大幅提升，GPU通用計算方法逐漸成為研究的熱點。而計算統一設備架構(compute unified device architecture, CUDA)技術基于擴展的C語言進一步改進了編程模型和內存模型，使其具有更為廣闊的應用前景。基于CUDA本文提出了一種醫學圖像GRBF二維/三維插值的方法，可望大大節約運算時間，為醫學圖像處理過程中高精度圖像縮放和三維重建應用奠定了基礎。
1 二維圖像縮放插值
1.1 GRBF插值原理
目前普遍采用的圖像縮放方法是將輸入圖像用連續函數進行插值擬合，再根據圖像縮放的倍數要求完成從輸出圖像到輸入圖像的坐標逆映射，即在已經擬合好的函數中對逆映射點進行重采樣，完成目標圖像各像素點的像素值計算[7]。給定一個一元函數φ：R+→R 在定義域x∈Rd上, 所有形如φ(x-c)=φ‖x-c‖)及其線性組合張成的函數空間稱為由函數φ的徑向基函數空間。在一定的條件下, 人們只要取{xj}兩兩不同, {φ(x-xj)}是線性無關的, 從而形成徑向基函數空間中某子空間的一組基。當{xj}幾乎充滿Rd時, {φ(x-xj)}及其線性組合可以逼近幾乎任何函數[8]。
對于二維GRBF可以采用表達式(1)，其中‖·‖為歐幾里德范數。
| $\begin{align} & \varphi \left( \|x-{{x}_{i}}\|,\|y-{{y}_{j}}\| \right)= \\ & \text{exp(}-\frac{\|x-{{x}_{i}}{{\|}^{2}}+\|y-{{y}_{j}}{{\|}^{2}}}{2{{\sigma }^{2}}}) \\ \end{align}$ |
GRBF 插值就是基于平移操作構造基函數系并尋找連續插值函數S(x, y)，如公式(2)所示：
| $S\left( x,y \right)=\sum\limits_{j=1}^{m}{\sum\limits_{i=1}^{n}{{{a}_{i,j}}}}\centerdot \varphi \left( \|x-{{x}_{i}}\|,\|y-{{y}_{j}}\| \right)$ |
GRBF插值正向計算過程是使給定的原始離散圖像數據滿足S(xj，yj)=fi,j,。其中m為圖像的高度, n為圖像的寬度。根據這些數據值以及GRBF基函數，便可以求出控制系數ai,j。GRBF插值逆向計算過程是用求得的控制系數ai,j構造出連續的高斯插值曲面，然后根據圖像的縮放倍數將目標圖像每一點的坐標逆映射到高斯插值曲面的坐標系中，對圖像進行重采樣。
1.2 基于CUDA的GRBF二維圖像插值與優化
CUDA是NVIDIA公司在2007年為其GPU開發的一個并行計算統一設備架構，在大規模通用計算領域，它采用單指令多線程(single instruction multiple threads, SIMT)的執行模型，即一條指令控制多個線程同時執行，并且在計算能力l.x的GPU的運算上不需要借助傳統圖形API, 只需要遵守固定的流水管線。CUDA的核心機制包括三個模型：硬件模型、編程模型和內存模型。在GPU硬件上，包含有很多個流多處理器(stream multiprocessor, SM), 每個SM上又包含有8個流處理器(stream processor, SP), SP是硬件執行的最小單位。編程人員根據需要，將多個任務分成多個網格(Grid)，每個網格又劃分成多個線程塊(Block), 每個線程塊又劃分成多個線程(Thread)，每個線程塊以及線程都有自己唯一的標識ID，在程序運行時，Block映射到SM上，Thread映射到SP上。在內存模型上，全局內存、常量內存、紋理內存是對所有線程透明的，全局內存容量較大，但讀寫時間較長；常量內存、紋理內存雖然訪問速度較快，但它們是只讀的，且容量不大。共享內存對同一個線程塊里的所有線程透明，它是片上存儲器，讀寫速度跟寄存器相當。
根據原始圖像數據的像素值以及坐標，可以將公式(2)表示成矩陣乘法的形式：
| $F=\Phi \centerdot A,$ |
其中F=(F11, F21, ……Fn1……Fnm) 是圖像的像素值矩陣，Φ是一個寬度和高度均為n·m的GRBF基函數方陣，它是關于主對角線對稱的。
| $\Phi =\left\{ \begin{matrix} \varphi (\left\| {{x}_{1}}-{{x}_{1}} \right\|,\left\| {{y}_{1}}-{{y}_{1}} \right\|) & \varphi (\left\| {{x}_{2}}-{{x}_{1}} \right\|,\left\| {{y}_{1}}-{{y}_{1}} \right\|) & \cdots \cdots & \varphi (\left\| {{x}_{n}}-{{x}_{1}} \right\|,\left\| {{y}_{m}}-{{y}_{1}} \right\|) \\ \varphi \left( \left\| {{x}_{1}}-{{x}_{2}} \right\|,\left\| {{y}_{1}}-{{y}_{1}} \right\| \right) & \varphi \left( \left\| {{x}_{2}}-{{x}_{2}} \right\|,\left\| {{y}_{1}}-{{y}_{1}} \right\| \right) & \cdots \cdots & \varphi \left( \left\| {{x}_{n}}-{{x}_{2}} \right\|,\left\| {{y}_{m}}-{{y}_{1}} \right\| \right) \\ \vdots & \vdots & \ddots & \vdots \\ \varphi \left( \left\| {{x}_{1}}-{{x}_{n}} \right\|,\left\| {{y}_{1}}-{{y}_{m}} \right\| \right) & \varphi \left( \left\| {{x}_{2}}-{{x}_{n}} \right\|,\left\| {{y}_{1}}-{{y}_{m}} \right\| \right) & \cdots \cdots & \varphi \left( \left\| {{x}_{n}}-{{x}_{n}} \right\|,\left\| {{y}_{m}}-{{y}_{m}} \right\| \right) \\ \end{matrix} \right\}$ |
根據GRBF插值正向求控制系數矩陣A，以及逆向重采樣的兩個過程，可以將CUDA編程模型設計成兩個Grid。正向計算的過程在數學描述上就是解一個大規模的線性方程組。目前求解線性方程組的一般方法包括直接法和迭代法。迭代法求解時行與行方程之間的迭代都是獨立的，非常適用于并行計算，但是GRBF基函數方陣并不是嚴格主對角占優[9]，實驗結果證明，使用雅克比迭代法求出來的A是發散的。
楊梅等[10]提出一個grid-strip-group-block-thread五層CUDA編程模型結構來進行無回代過程的高斯-約當主元消去，能夠在一定程度上降低過多線程頻繁訪問顯存所消耗的時間，但此模型實現起來步驟較多，且在group組織結構消元過程中并沒有考慮到合并訪存的要求。為了充分降低線程訪問全局內存的時間，本文采用有回代過程的列主元高斯消去法。如果將GRBF基函數方陣Φ和像素值矩陣F表示成如圖 1(a)所示的增廣矩陣，以第i行(i∈(1, n·m-1))為主元進行消元時，自然考慮設計n·m-i個線程，每個線程分別負責i+1行到最后1行與第i行進行消元。以i=1為例，線程1~(n·m-1)分別負責φ12，……φ1, nm所在行與φ11所在主元行進行消元，結果使φ12，……φ1, nm都消為0。當i遍歷其他取值時，最終使增廣矩陣化為上三角形式。每次消元過程可以通過將每個線程映射到SP來并行執行，但是因為增廣矩陣在顯存中是按照行的順序來存儲，線程1~(n·m-i)訪問的地址空間是不連續的，也就是說每個線程都要獨立的進行一次訪存通信，而一次訪存需要400~600個時鐘周期，當矩陣的維數較大時，線程的并行性完全可以被訪存延時掩蓋。實驗表明，以行為主元來消元時，GPU運行時間比CPU運行時間僅僅稍微有所降低。
圖1
并行高斯消去法示意圖
(a)行主元;(b)列主元
Figure1. Diagram of parallel Gaussian elimination method(a) row principle; (b) column principle
將圖 1(a)所示增廣矩陣進行轉置，便得到圖 1(b)所示矩陣，將以行為主元進行消元的方法改為以列為主元進行消元。這次以i=2為例，線程1~(n·m-2)分別負責φ23，……φ2, nm所在列與φ22所在主元列進行消元，使φ23，……φ2, nm都消元為0。因為這(n*m-2)個線程所訪存的地址空間在每次消元時都是連續的，在線程執行時，同一個線程束(wrap)中的16個線程就只需要1次顯存讀取，這充分利用了合并訪存的優勢，大大節約了時間，并且因為GRBF基函數矩陣是對稱的，實際操作時并不需要進行轉置。求解控制系數矩陣A就是采用列主元高斯消去法。編程時將增廣矩陣載入到顯存之中，設置每個block中有t個線程，t小于512且最好是16的整數倍，CPU控制消元的迭代循環，一共有(n·m-1)次。每次消元所需要的線程總數為(n·m-i)，每次循環block的大小動態的設置為((n·m-i)+(t-1))/t, 經過(n·m-1)次列主元迭代消元后，增廣矩成為了下三角形，最后將消元后的增廣矩陣傳回CPU，通過回代過程求解出控制系數矩陣A。
將求得的控制系數矩陣A代入公式(2)便得到了連續的高斯插值曲面，在GRBF插值逆向重采樣過程中，假設原始圖像在水平和垂直方向的像素坐標為(x, y)∈(0~n-1, 0~m-1)，縮放倍數分別為z1、z2，則目標圖像的像素坐標是(i, j)∈(0~z1·n-1, 0~z2·m-1)。接著將(i/z1, j/z2)的坐標值代入連續高斯插值曲面進行重采樣。由于目標圖像每一點像素值重采樣過程都是獨立的，可以設置每一個線程負責一個像素點重采樣，在CUDA編程模型第二個grid中設置每一個block中有(t, t)的二維線程陣列，則block的維度是((z1·n+t-1)/t, (z2·m+t-1)/t)。當控制系數矩陣A所占大小小于16 KB時，可以通過block中線程的ID值并行地將A載入到共享內存里面。最后所有線程就可以并行地計算其ID所對應的采樣點的像素值。
1.3 插值結果及分析
為了驗證GRBF插值在CUDA上實現的高效性，本節以顯微鏡下二維尿沉渣的圖像為例，分別實現了同一尺寸圖像下不同縮放倍數、不同尺寸圖像下同一縮放倍數的情況，給出CPU、GPU運行時間以及插值精度的相關信息。本節最后簡略地分析了GRBF形狀參數對插值精度的影響。
實驗中所使用的平臺為CPU：AMD Sempron 140，主頻2.7 GHz。GPU：Geforce GTX 280, 處理器頻率1 296 MHz, SP數量為240個。其中表 1顯示的是對于尺寸大小為60×60的尿沉渣圖像，縮放倍率在水平和垂直方向均取為0.5、2、4、8時在CPU以及GPU上運行的結果。其中變量CPU是CPU實現方式程序運行的總時間，CPU1是采用高斯消去的時間，CPU2是逆采樣時間，MSEC是CPU上GRBF插值的均方誤差，GPU是CUDA程序運行的總時間，GPU1、GPU2分別是CPU1和CPU2的并行實現，MSEG是GPU上GRBF插值的均方誤差，Acc是加速比，也就等于CPU/GPU。其中，基函數的標準差σ取為1。并行消元過程block中設置1×256個線程，重采樣過程block中設置16×16個線程。均方誤差的計算是將原始圖像映射到插值圖像進行求解。
表1
60×60尿沉渣圖像GRBF插值時間比較
Table1.
Time comparison of (60×60) urinary sediment image GRBF interpolation
從表 1可以看出，對于同一尺寸的圖像隨著縮放倍率的增大，加速比逐漸增大。計算時間的消耗主要體現在高斯消去過程，不同縮放倍率下CPU1和GPU1時間基本保持不變，這是因為GRBF基函數方陣是固定大小的(此例中為3 600階)，且運算時間和加速比均優于文獻[10]中的結果。隨著縮放倍率的增大，CPU2時間急劇增長，而GPU2時間則是緩慢增長，兩者相比有超過1 500倍的加速。并且對比CPU和GPU上GRBF插值的均方誤差，兩者都很小且數值比較接近，其都采用單精度浮點運算，插值圖像質量比較接近。圖 2(a)~(d)分別顯示了原圖、縮放倍率為0.5、2、4時的插值圖像。
圖2
插值圖像
(a)原圖;(b)縮放0.5倍;(c) 縮放2倍;(d) 縮放4倍
Figure2. Interpolated images(a) the original image; (b) 0.5 times scaling; (c) 2 times scaling; (d) 4 times scaling
表 2顯示的是不同尺寸尿沉渣圖像放大2倍時的GRBF插值計算時間比較。從表中可以看出，隨著圖像尺寸的增大，運算時間尤其是CPU1和GPU1大幅增長。當圖像尺寸為80×80時，控制系數矩陣A的大小超過了16 KB，不能將它放在共享內存里面，實驗中可以將它放在具有高速緩存的常量內存里。比較80×80與120×120時的數據發現，后者在保持較高插值精度的同時，CPU1和GPU1時間要比前者低。這是因為在對大尺寸圖像運用GRBF插值時，為了降低運算量，考慮到一幅圖像距離近的像素點相關性大，距離遠的像素點相關性小，可以對圖像進行二維分塊，然后對每個小塊分別進行GRBF插值，最后把插值結果連接起來，如果把120×120劃分成4個60×60小塊，GRBF插值的結果如圖 3(b)所示，在圖像水平和垂直中線處出現細長的干擾線，即在塊連接處出現失真。解決辦法是，保持相鄰小塊有足夠的重疊區域，比如4個62×62小塊，連接時相鄰塊運用自然縫合算法[11]。結果如圖 3(c)所示。
表2
不同尺寸尿沉渣圖像GRBF插值時間比較
Table2.
Time comparison between different size urinary sediment image GRBF interpolation
圖3
120×120尿沉渣插值圖像
(a)120×120原圖;(b) 連接處失真圖像;(c)改進后圖像
Figure3. 120×120 Interpolated urinary sediment images(a) the 120×120 original image; (b) the edge distorted image; (c) the improved image
在以上兩個實驗中，式(1)中的σ均取為1，而當形狀參數1/(2σ2)取不同值時將得到不同的GRBF插值精度。圖 4顯示了對于60×60的尿沉渣圖像放大2倍時，σ取不同值時，在GPU上得出來的均方誤差關系圖。隨著σ的增大，形狀參數減小，使GRBF形狀趨于平坦，當σ大于等于1.6時，均方誤差變得較大且不穩定，插值圖像完全失真。
圖4
σ-MSEG關系圖
Figure4.
Diagram of σ-MSEG
2 三維斷層圖像序列插值
2.1 LSCM圖像提高軸向分辨率
激光掃描共焦顯微鏡(laser scanning confocal microscopy, LSCM)是一種用于分析細胞學的儀器。可以對活細胞進行分層掃描，進行三維重建。通過熒光標記可以對細胞內微細結構進行定性、定量的觀察。受載物臺沿軸向移動步距的影響，LSCM斷層圖像序列軸向的分辨率要低于水平方向的分辨率，在三維重建過程中，可能使欲觀察的結構發生比例變形。使用GRBF插值可以在斷層圖像之間插入新的圖像，從而提高軸向分辨率，消除變形的影響。
本節中采用的原始數據是Rubelab實驗室提供的小鼠耳蝸螺旋器LSCM斷層掃描圖像。切片尺寸為512×512，一共有192層，層間距離為0.4 μm。在對512×512×192的三維數據進行GRBF插值時，首先將其均勻地劃分成多個體空間，然后對每個體空間進行正向和逆向求解。在正向求解過程中，如果每個體空間都進行一次列主元高斯消元，對于擁有大量數據的三維插值將耗費比較多的時間。如果對GRBF基函數方陣進行一次求逆，然后用逆矩陣和每個體空間數據進行乘積，將大大節約時間。LSCM斷層掃描圖像GRBF插值流程圖如圖 5所示。
圖5
LSCM斷層掃描圖像GRBF插值流程圖
Figure5.
GRBF interpolation process chart of LSCM tomographic images
2.2 三維插值CUDA實現及結果
在小鼠耳蝸螺旋器LSCM斷層掃描圖像數據中，設定體空間大小為16×16×8，則原始數據被均勻劃分為24 576個子體空間。應注意各個子體空間大小的選擇要綜合考慮顯存大小、block維度上限以及運算復雜度等的影響。實際操作中為了消除子體空間在連接處的失真現象，可以保持上下體空間有一層重疊層，然后在插值圖像的重疊區域采用自然縫合算法。三維GRBF可以在式(1)的基礎上增加z分量, σ仍然取為1。在用CUDA實現GRBF方陣求逆時，將此方陣與同尺寸單位矩陣E上下組合，采用列主元高斯消去法，將GRBF方陣消為單位矩陣，E則變為GRBF方陣的逆矩陣。逆矩陣與子體空間數據相乘時，也要充分考慮合并訪存的影響。表 3顯示了CPU與GPU運算的對比。在CPU方面，時間主要消耗在重采樣過程，且采用雙精度浮點運算，MSEC等于0。在GPU方面，GPU1、GPU2、GPU3分別代表求逆、相乘和重采樣時間，CUDA程序大部分時間消耗在三維數據的讀寫上面。比較CPU和GPU上GRBF插值的均方誤差，后者要大于前者，表明GPU上圖像插值質量要稍低于CPU上圖像插值質量。圖 6顯示的是對小鼠耳蝸螺旋器LSCM斷層掃描圖像通過CUDA進行GRBF插值得到的結果，圖 6(a)從左至右依次顯示的是軸向分辨率提高2倍時的第一層圖像、插值圖像、第二層圖像。從圖 6(a)可以看出，插值圖像結合了上下兩層圖像的結構特征。圖 6(b)從左至右依次顯示的是LSCM原始圖像、軸向分辨率提高2倍、軸向分辨率提高3倍的三維重建插值圖像。圖 6(b)中原始三維圖像熒光標記結構形狀扁平，且結構之間較密集。經過GRBF插值提高軸向分辨率后，熒光標記結構形狀變得修長，各結構之間能較為清楚地區分開來，且三維插值圖像不存在子體空間連接處的失真現象。
圖6
LSCM斷層掃描圖像插值結果
(a)第一層圖像、插值圖像、第二層圖像;(b)原始三維圖像、軸向分辨率提高2倍和3倍插值圖像
Figure6. Results of LSCM tomographic images GRBF interpolation(a) the first layer image, interpolated image and second layer image; (b) the original 3D image, the axial resolution 2 times and 3 times increased interpolated image
表3
LSCM斷層掃描圖像GRBF插值時間比較
Table3.
Time comparison of LSCM tomographic images GRBF interpolation
3 結束語
本文在基于CUDA對二維、三維醫學圖像的GRBF插值快速實現過程中，充分考慮到合并訪存、共享內存對GPU運算時間的影響，采用列主元的高斯消去法來求解控制系數矩陣。為了進一步降低運算量并有效提高運算速度，提出對數據空間進行分塊、分體的思想，并且在實現過程中采用自然縫合算法消除連接處的失真現象。根據上文中實驗數據可以看出，無論對于兩維還是三維生物醫學圖像，GRBF插值在正向求解控制系數矩陣A和逆向重采樣過程都得到了明顯的加速。這表明基于CUDA來實現醫學圖像的GRBF快速插值是完全可行的，能夠彌補GRBF插值速度慢的固有缺陷，可有效拓展GRBF插值方法的應用范圍。
引言
在醫學圖像處理和分析過程中，往往需要基于原始的具有較低分辨率的圖像數據得到具有更高分辨率的圖像數據。如對于二維醫學圖像而言，人們常常需要將原始圖像放大，從而更好地觀察細胞或者腫瘤的形狀、位置以及和周圍物質的相對關系等等。對于三維圖像序列，由于成像設備本身精度、成像時間以及存儲器容量的限制，采集過程中獲得的任意兩切片間的空間距離均大于片內各像素之間的距離，即存在不同方向上分辨率的差異。當用所獲得的圖像序列來重建三維結構時，可能會使重建結果產生比例變形。這往往都需要通過插值過程予以修正。因此，圖像插值在二維醫學圖像縮放和三維圖像重建過程中是一種常用的方法。一般說來，常用的醫學圖像插值方法可以分為基于圖像的插值方法、基于對象的插值方法和基于小波的插值方法[1]。基于圖像的插值方法包括最近鄰插值、線性插值、匹配插值[2]等；基于對象的插值方法主要是采用形態學中的膨脹與腐蝕；基于小波的插值方法是通過小波變換獲得相應小波系數后對斷層圖像進行強度和位置插值[3]。本文基于高斯徑向基(Gaussian radial basis function, GRBF)插值的方法是一種基于圖像的插值方法。目前較多地應用于曲線和曲面的插值逼近場合，繆報通等[4]基于徑向基神經網絡也將GRBF用于散亂數據插值； Fornberg等[5]采用RBF-QR方法消除形狀參數對GRBF插值穩定性計算的影響，Liao等[6]擴展了有理GRBF，結合形狀參數實現了曲線和曲面的插值。如果說最近鄰插值、線性插值、樣條插值是一種“局部”插值，即插值點只與周圍少量數據點相關，那么GRBF插值則是一種相對“全局”的插值，高精度的優點使GRBF插值圖像不存在邊緣鋸齒現象。但基于GRBF的插值過程因其數據計算量龐大，需要的運算時間比較長，這限制了其在圖像縮放和重建過程中的應用。
隨著圖形處理單元(GPU)并行運算能力的大幅提升，GPU通用計算方法逐漸成為研究的熱點。而計算統一設備架構(compute unified device architecture, CUDA)技術基于擴展的C語言進一步改進了編程模型和內存模型，使其具有更為廣闊的應用前景。基于CUDA本文提出了一種醫學圖像GRBF二維/三維插值的方法，可望大大節約運算時間，為醫學圖像處理過程中高精度圖像縮放和三維重建應用奠定了基礎。
1 二維圖像縮放插值
1.1 GRBF插值原理
目前普遍采用的圖像縮放方法是將輸入圖像用連續函數進行插值擬合，再根據圖像縮放的倍數要求完成從輸出圖像到輸入圖像的坐標逆映射，即在已經擬合好的函數中對逆映射點進行重采樣，完成目標圖像各像素點的像素值計算[7]。給定一個一元函數φ：R+→R 在定義域x∈Rd上, 所有形如φ(x-c)=φ‖x-c‖)及其線性組合張成的函數空間稱為由函數φ的徑向基函數空間。在一定的條件下, 人們只要取{xj}兩兩不同, {φ(x-xj)}是線性無關的, 從而形成徑向基函數空間中某子空間的一組基。當{xj}幾乎充滿Rd時, {φ(x-xj)}及其線性組合可以逼近幾乎任何函數[8]。
對于二維GRBF可以采用表達式(1)，其中‖·‖為歐幾里德范數。
| $\begin{align} & \varphi \left( \|x-{{x}_{i}}\|,\|y-{{y}_{j}}\| \right)= \\ & \text{exp(}-\frac{\|x-{{x}_{i}}{{\|}^{2}}+\|y-{{y}_{j}}{{\|}^{2}}}{2{{\sigma }^{2}}}) \\ \end{align}$ |
GRBF 插值就是基于平移操作構造基函數系并尋找連續插值函數S(x, y)，如公式(2)所示：
| $S\left( x,y \right)=\sum\limits_{j=1}^{m}{\sum\limits_{i=1}^{n}{{{a}_{i,j}}}}\centerdot \varphi \left( \|x-{{x}_{i}}\|,\|y-{{y}_{j}}\| \right)$ |
GRBF插值正向計算過程是使給定的原始離散圖像數據滿足S(xj，yj)=fi,j,。其中m為圖像的高度, n為圖像的寬度。根據這些數據值以及GRBF基函數，便可以求出控制系數ai,j。GRBF插值逆向計算過程是用求得的控制系數ai,j構造出連續的高斯插值曲面，然后根據圖像的縮放倍數將目標圖像每一點的坐標逆映射到高斯插值曲面的坐標系中，對圖像進行重采樣。
1.2 基于CUDA的GRBF二維圖像插值與優化
CUDA是NVIDIA公司在2007年為其GPU開發的一個并行計算統一設備架構，在大規模通用計算領域，它采用單指令多線程(single instruction multiple threads, SIMT)的執行模型，即一條指令控制多個線程同時執行，并且在計算能力l.x的GPU的運算上不需要借助傳統圖形API, 只需要遵守固定的流水管線。CUDA的核心機制包括三個模型：硬件模型、編程模型和內存模型。在GPU硬件上，包含有很多個流多處理器(stream multiprocessor, SM), 每個SM上又包含有8個流處理器(stream processor, SP), SP是硬件執行的最小單位。編程人員根據需要，將多個任務分成多個網格(Grid)，每個網格又劃分成多個線程塊(Block), 每個線程塊又劃分成多個線程(Thread)，每個線程塊以及線程都有自己唯一的標識ID，在程序運行時，Block映射到SM上，Thread映射到SP上。在內存模型上，全局內存、常量內存、紋理內存是對所有線程透明的，全局內存容量較大，但讀寫時間較長；常量內存、紋理內存雖然訪問速度較快，但它們是只讀的，且容量不大。共享內存對同一個線程塊里的所有線程透明，它是片上存儲器，讀寫速度跟寄存器相當。
根據原始圖像數據的像素值以及坐標，可以將公式(2)表示成矩陣乘法的形式：
| $F=\Phi \centerdot A,$ |
其中F=(F11, F21, ……Fn1……Fnm) 是圖像的像素值矩陣，Φ是一個寬度和高度均為n·m的GRBF基函數方陣，它是關于主對角線對稱的。
| $\Phi =\left\{ \begin{matrix} \varphi (\left\| {{x}_{1}}-{{x}_{1}} \right\|,\left\| {{y}_{1}}-{{y}_{1}} \right\|) & \varphi (\left\| {{x}_{2}}-{{x}_{1}} \right\|,\left\| {{y}_{1}}-{{y}_{1}} \right\|) & \cdots \cdots & \varphi (\left\| {{x}_{n}}-{{x}_{1}} \right\|,\left\| {{y}_{m}}-{{y}_{1}} \right\|) \\ \varphi \left( \left\| {{x}_{1}}-{{x}_{2}} \right\|,\left\| {{y}_{1}}-{{y}_{1}} \right\| \right) & \varphi \left( \left\| {{x}_{2}}-{{x}_{2}} \right\|,\left\| {{y}_{1}}-{{y}_{1}} \right\| \right) & \cdots \cdots & \varphi \left( \left\| {{x}_{n}}-{{x}_{2}} \right\|,\left\| {{y}_{m}}-{{y}_{1}} \right\| \right) \\ \vdots & \vdots & \ddots & \vdots \\ \varphi \left( \left\| {{x}_{1}}-{{x}_{n}} \right\|,\left\| {{y}_{1}}-{{y}_{m}} \right\| \right) & \varphi \left( \left\| {{x}_{2}}-{{x}_{n}} \right\|,\left\| {{y}_{1}}-{{y}_{m}} \right\| \right) & \cdots \cdots & \varphi \left( \left\| {{x}_{n}}-{{x}_{n}} \right\|,\left\| {{y}_{m}}-{{y}_{m}} \right\| \right) \\ \end{matrix} \right\}$ |
根據GRBF插值正向求控制系數矩陣A，以及逆向重采樣的兩個過程，可以將CUDA編程模型設計成兩個Grid。正向計算的過程在數學描述上就是解一個大規模的線性方程組。目前求解線性方程組的一般方法包括直接法和迭代法。迭代法求解時行與行方程之間的迭代都是獨立的，非常適用于并行計算，但是GRBF基函數方陣并不是嚴格主對角占優[9]，實驗結果證明，使用雅克比迭代法求出來的A是發散的。
楊梅等[10]提出一個grid-strip-group-block-thread五層CUDA編程模型結構來進行無回代過程的高斯-約當主元消去，能夠在一定程度上降低過多線程頻繁訪問顯存所消耗的時間，但此模型實現起來步驟較多，且在group組織結構消元過程中并沒有考慮到合并訪存的要求。為了充分降低線程訪問全局內存的時間，本文采用有回代過程的列主元高斯消去法。如果將GRBF基函數方陣Φ和像素值矩陣F表示成如圖 1(a)所示的增廣矩陣，以第i行(i∈(1, n·m-1))為主元進行消元時，自然考慮設計n·m-i個線程，每個線程分別負責i+1行到最后1行與第i行進行消元。以i=1為例，線程1~(n·m-1)分別負責φ12，……φ1, nm所在行與φ11所在主元行進行消元，結果使φ12，……φ1, nm都消為0。當i遍歷其他取值時，最終使增廣矩陣化為上三角形式。每次消元過程可以通過將每個線程映射到SP來并行執行，但是因為增廣矩陣在顯存中是按照行的順序來存儲，線程1~(n·m-i)訪問的地址空間是不連續的，也就是說每個線程都要獨立的進行一次訪存通信，而一次訪存需要400~600個時鐘周期，當矩陣的維數較大時，線程的并行性完全可以被訪存延時掩蓋。實驗表明，以行為主元來消元時，GPU運行時間比CPU運行時間僅僅稍微有所降低。
圖1
并行高斯消去法示意圖
(a)行主元;(b)列主元
Figure1. Diagram of parallel Gaussian elimination method(a) row principle; (b) column principle
將圖 1(a)所示增廣矩陣進行轉置，便得到圖 1(b)所示矩陣，將以行為主元進行消元的方法改為以列為主元進行消元。這次以i=2為例，線程1~(n·m-2)分別負責φ23，……φ2, nm所在列與φ22所在主元列進行消元，使φ23，……φ2, nm都消元為0。因為這(n*m-2)個線程所訪存的地址空間在每次消元時都是連續的，在線程執行時，同一個線程束(wrap)中的16個線程就只需要1次顯存讀取，這充分利用了合并訪存的優勢，大大節約了時間，并且因為GRBF基函數矩陣是對稱的，實際操作時并不需要進行轉置。求解控制系數矩陣A就是采用列主元高斯消去法。編程時將增廣矩陣載入到顯存之中，設置每個block中有t個線程，t小于512且最好是16的整數倍，CPU控制消元的迭代循環，一共有(n·m-1)次。每次消元所需要的線程總數為(n·m-i)，每次循環block的大小動態的設置為((n·m-i)+(t-1))/t, 經過(n·m-1)次列主元迭代消元后，增廣矩成為了下三角形，最后將消元后的增廣矩陣傳回CPU，通過回代過程求解出控制系數矩陣A。
將求得的控制系數矩陣A代入公式(2)便得到了連續的高斯插值曲面，在GRBF插值逆向重采樣過程中，假設原始圖像在水平和垂直方向的像素坐標為(x, y)∈(0~n-1, 0~m-1)，縮放倍數分別為z1、z2，則目標圖像的像素坐標是(i, j)∈(0~z1·n-1, 0~z2·m-1)。接著將(i/z1, j/z2)的坐標值代入連續高斯插值曲面進行重采樣。由于目標圖像每一點像素值重采樣過程都是獨立的，可以設置每一個線程負責一個像素點重采樣，在CUDA編程模型第二個grid中設置每一個block中有(t, t)的二維線程陣列，則block的維度是((z1·n+t-1)/t, (z2·m+t-1)/t)。當控制系數矩陣A所占大小小于16 KB時，可以通過block中線程的ID值并行地將A載入到共享內存里面。最后所有線程就可以并行地計算其ID所對應的采樣點的像素值。
1.3 插值結果及分析
為了驗證GRBF插值在CUDA上實現的高效性，本節以顯微鏡下二維尿沉渣的圖像為例，分別實現了同一尺寸圖像下不同縮放倍數、不同尺寸圖像下同一縮放倍數的情況，給出CPU、GPU運行時間以及插值精度的相關信息。本節最后簡略地分析了GRBF形狀參數對插值精度的影響。
實驗中所使用的平臺為CPU：AMD Sempron 140，主頻2.7 GHz。GPU：Geforce GTX 280, 處理器頻率1 296 MHz, SP數量為240個。其中表 1顯示的是對于尺寸大小為60×60的尿沉渣圖像，縮放倍率在水平和垂直方向均取為0.5、2、4、8時在CPU以及GPU上運行的結果。其中變量CPU是CPU實現方式程序運行的總時間，CPU1是采用高斯消去的時間，CPU2是逆采樣時間，MSEC是CPU上GRBF插值的均方誤差，GPU是CUDA程序運行的總時間，GPU1、GPU2分別是CPU1和CPU2的并行實現，MSEG是GPU上GRBF插值的均方誤差，Acc是加速比，也就等于CPU/GPU。其中，基函數的標準差σ取為1。并行消元過程block中設置1×256個線程，重采樣過程block中設置16×16個線程。均方誤差的計算是將原始圖像映射到插值圖像進行求解。
表1
60×60尿沉渣圖像GRBF插值時間比較
Table1.
Time comparison of (60×60) urinary sediment image GRBF interpolation
從表 1可以看出，對于同一尺寸的圖像隨著縮放倍率的增大，加速比逐漸增大。計算時間的消耗主要體現在高斯消去過程，不同縮放倍率下CPU1和GPU1時間基本保持不變，這是因為GRBF基函數方陣是固定大小的(此例中為3 600階)，且運算時間和加速比均優于文獻[10]中的結果。隨著縮放倍率的增大，CPU2時間急劇增長，而GPU2時間則是緩慢增長，兩者相比有超過1 500倍的加速。并且對比CPU和GPU上GRBF插值的均方誤差，兩者都很小且數值比較接近，其都采用單精度浮點運算，插值圖像質量比較接近。圖 2(a)~(d)分別顯示了原圖、縮放倍率為0.5、2、4時的插值圖像。
圖2
插值圖像
(a)原圖;(b)縮放0.5倍;(c) 縮放2倍;(d) 縮放4倍
Figure2. Interpolated images(a) the original image; (b) 0.5 times scaling; (c) 2 times scaling; (d) 4 times scaling
表 2顯示的是不同尺寸尿沉渣圖像放大2倍時的GRBF插值計算時間比較。從表中可以看出，隨著圖像尺寸的增大，運算時間尤其是CPU1和GPU1大幅增長。當圖像尺寸為80×80時，控制系數矩陣A的大小超過了16 KB，不能將它放在共享內存里面，實驗中可以將它放在具有高速緩存的常量內存里。比較80×80與120×120時的數據發現，后者在保持較高插值精度的同時，CPU1和GPU1時間要比前者低。這是因為在對大尺寸圖像運用GRBF插值時，為了降低運算量，考慮到一幅圖像距離近的像素點相關性大，距離遠的像素點相關性小，可以對圖像進行二維分塊，然后對每個小塊分別進行GRBF插值，最后把插值結果連接起來，如果把120×120劃分成4個60×60小塊，GRBF插值的結果如圖 3(b)所示，在圖像水平和垂直中線處出現細長的干擾線，即在塊連接處出現失真。解決辦法是，保持相鄰小塊有足夠的重疊區域，比如4個62×62小塊，連接時相鄰塊運用自然縫合算法[11]。結果如圖 3(c)所示。
表2
不同尺寸尿沉渣圖像GRBF插值時間比較
Table2.
Time comparison between different size urinary sediment image GRBF interpolation
圖3
120×120尿沉渣插值圖像
(a)120×120原圖;(b) 連接處失真圖像;(c)改進后圖像
Figure3. 120×120 Interpolated urinary sediment images(a) the 120×120 original image; (b) the edge distorted image; (c) the improved image
在以上兩個實驗中，式(1)中的σ均取為1，而當形狀參數1/(2σ2)取不同值時將得到不同的GRBF插值精度。圖 4顯示了對于60×60的尿沉渣圖像放大2倍時，σ取不同值時，在GPU上得出來的均方誤差關系圖。隨著σ的增大，形狀參數減小，使GRBF形狀趨于平坦，當σ大于等于1.6時，均方誤差變得較大且不穩定，插值圖像完全失真。
圖4
σ-MSEG關系圖
Figure4.
Diagram of σ-MSEG
2 三維斷層圖像序列插值
2.1 LSCM圖像提高軸向分辨率
激光掃描共焦顯微鏡(laser scanning confocal microscopy, LSCM)是一種用于分析細胞學的儀器。可以對活細胞進行分層掃描，進行三維重建。通過熒光標記可以對細胞內微細結構進行定性、定量的觀察。受載物臺沿軸向移動步距的影響，LSCM斷層圖像序列軸向的分辨率要低于水平方向的分辨率，在三維重建過程中，可能使欲觀察的結構發生比例變形。使用GRBF插值可以在斷層圖像之間插入新的圖像，從而提高軸向分辨率，消除變形的影響。
本節中采用的原始數據是Rubelab實驗室提供的小鼠耳蝸螺旋器LSCM斷層掃描圖像。切片尺寸為512×512，一共有192層，層間距離為0.4 μm。在對512×512×192的三維數據進行GRBF插值時，首先將其均勻地劃分成多個體空間，然后對每個體空間進行正向和逆向求解。在正向求解過程中，如果每個體空間都進行一次列主元高斯消元，對于擁有大量數據的三維插值將耗費比較多的時間。如果對GRBF基函數方陣進行一次求逆，然后用逆矩陣和每個體空間數據進行乘積，將大大節約時間。LSCM斷層掃描圖像GRBF插值流程圖如圖 5所示。
圖5
LSCM斷層掃描圖像GRBF插值流程圖
Figure5.
GRBF interpolation process chart of LSCM tomographic images
2.2 三維插值CUDA實現及結果
在小鼠耳蝸螺旋器LSCM斷層掃描圖像數據中，設定體空間大小為16×16×8，則原始數據被均勻劃分為24 576個子體空間。應注意各個子體空間大小的選擇要綜合考慮顯存大小、block維度上限以及運算復雜度等的影響。實際操作中為了消除子體空間在連接處的失真現象，可以保持上下體空間有一層重疊層，然后在插值圖像的重疊區域采用自然縫合算法。三維GRBF可以在式(1)的基礎上增加z分量, σ仍然取為1。在用CUDA實現GRBF方陣求逆時，將此方陣與同尺寸單位矩陣E上下組合，采用列主元高斯消去法，將GRBF方陣消為單位矩陣，E則變為GRBF方陣的逆矩陣。逆矩陣與子體空間數據相乘時，也要充分考慮合并訪存的影響。表 3顯示了CPU與GPU運算的對比。在CPU方面，時間主要消耗在重采樣過程，且采用雙精度浮點運算，MSEC等于0。在GPU方面，GPU1、GPU2、GPU3分別代表求逆、相乘和重采樣時間，CUDA程序大部分時間消耗在三維數據的讀寫上面。比較CPU和GPU上GRBF插值的均方誤差，后者要大于前者，表明GPU上圖像插值質量要稍低于CPU上圖像插值質量。圖 6顯示的是對小鼠耳蝸螺旋器LSCM斷層掃描圖像通過CUDA進行GRBF插值得到的結果，圖 6(a)從左至右依次顯示的是軸向分辨率提高2倍時的第一層圖像、插值圖像、第二層圖像。從圖 6(a)可以看出，插值圖像結合了上下兩層圖像的結構特征。圖 6(b)從左至右依次顯示的是LSCM原始圖像、軸向分辨率提高2倍、軸向分辨率提高3倍的三維重建插值圖像。圖 6(b)中原始三維圖像熒光標記結構形狀扁平，且結構之間較密集。經過GRBF插值提高軸向分辨率后，熒光標記結構形狀變得修長，各結構之間能較為清楚地區分開來，且三維插值圖像不存在子體空間連接處的失真現象。
圖6
LSCM斷層掃描圖像插值結果
(a)第一層圖像、插值圖像、第二層圖像;(b)原始三維圖像、軸向分辨率提高2倍和3倍插值圖像
Figure6. Results of LSCM tomographic images GRBF interpolation(a) the first layer image, interpolated image and second layer image; (b) the original 3D image, the axial resolution 2 times and 3 times increased interpolated image
表3
LSCM斷層掃描圖像GRBF插值時間比較
Table3.
Time comparison of LSCM tomographic images GRBF interpolation
3 結束語
本文在基于CUDA對二維、三維醫學圖像的GRBF插值快速實現過程中，充分考慮到合并訪存、共享內存對GPU運算時間的影響，采用列主元的高斯消去法來求解控制系數矩陣。為了進一步降低運算量并有效提高運算速度，提出對數據空間進行分塊、分體的思想，并且在實現過程中采用自然縫合算法消除連接處的失真現象。根據上文中實驗數據可以看出，無論對于兩維還是三維生物醫學圖像，GRBF插值在正向求解控制系數矩陣A和逆向重采樣過程都得到了明顯的加速。這表明基于CUDA來實現醫學圖像的GRBF快速插值是完全可行的，能夠彌補GRBF插值速度慢的固有缺陷，可有效拓展GRBF插值方法的應用范圍。
