基于曲線投影模型的電子斷層三維重構(gòu)并行算法

張靜蓉萬(wàn)曉華張 法

1(中國(guó)科學(xué)院計(jì)算技術(shù)研究所 北京 100190)

2(中國(guó)科學(xué)院大學(xué) 北京 100049)

基于曲線投影模型的電子斷層三維重構(gòu)并行算法

張靜蓉1,2萬(wàn)曉華1張 法1

1(中國(guó)科學(xué)院計(jì)算技術(shù)研究所 北京 100190)

2(中國(guó)科學(xué)院大學(xué) 北京 100049)

大尺度高精度的電子斷層三維重構(gòu)可以獲得在更大視角下的生物大分子三維結(jié)構(gòu)的細(xì)節(jié)信息。但研究尺度的增大給獲取高精度重構(gòu)結(jié)果和縮短數(shù)據(jù)處理時(shí)間帶來(lái)了巨大的挑戰(zhàn)。TxBR 提出的曲線模型顯著提高了重構(gòu)的精度，但其計(jì)算比直線模型更復(fù)雜耗時(shí)，且在曲線模型下，之前的并行策略不再可行。針對(duì)這一問(wèn)題，提出了一種在 GPU 平臺(tái)上實(shí)現(xiàn)的針對(duì)曲線模型的分塊迭代并行算法。通過(guò)對(duì)曲線模型的研究發(fā)現(xiàn)，曲線模型具有一定的空間局域性，利用這種性質(zhì)提出了一種縱向的分塊方式。在算法的實(shí)現(xiàn)階段，提出一個(gè)基于頁(yè)的數(shù)據(jù)傳輸策略，從而能夠去除冗余的數(shù)據(jù)傳輸，減少數(shù)據(jù)傳輸帶來(lái)的時(shí)間消耗。實(shí)驗(yàn)結(jié)果顯示，本算法可接近 40 倍的加速比。

電子斷層；迭代重構(gòu)算法；分塊；曲線投影模型；GPU 并行

1 引 言

電子斷層三維重構(gòu)(Electron Tomography， ET)是利用生物樣品的電子顯微投影圖像重構(gòu)出生物體三維結(jié)構(gòu)的技術(shù)。其主要策略是將制備好的生物樣品放在透射電子顯微鏡中，讓樣品沿著一個(gè)與電子束垂直的軸旋轉(zhuǎn)，每旋轉(zhuǎn)到一定角度，拍攝得到樣品的一張二維投影圖片[1]。然后對(duì)這一系列二維投影圖像進(jìn)行對(duì)位，利用三維重構(gòu)算法得到樣品的三維結(jié)構(gòu)。由于透射顯微鏡的物理限制，樣品的旋轉(zhuǎn)角度范圍在(－60°，＋60°)或(－70°，＋70°)之間，同時(shí)由于樣品可以承受的電子劑量有限，投影圖像往往噪聲很大。與直接重構(gòu)算法相比，迭代重構(gòu)算法可以更好地處理ET 中有信息缺失和高噪聲的圖像。

隨著投影圖像尺度的增大，大尺度高精度的電子斷層三維重構(gòu)研究越來(lái)越受到關(guān)注。目前大多數(shù)的重構(gòu)算法都采用直線投影模型。然而受到磁場(chǎng)的影響，電子束的軌跡是曲線，所以直線投影模型并不準(zhǔn)確，由此產(chǎn)生的誤差會(huì)隨著投影圖像尺度的增加變得更加明顯。為了減少這一誤差，研究人員提出了一種全局的非線性投影模型—TxBR[2]。與直線模型相比，這一模型可以提高重構(gòu)的精度，但其計(jì)算更加復(fù)雜耗時(shí)。

在大尺度高精度的電子斷層重構(gòu)中，我們常常使用高性能計(jì)算來(lái)幫助解決計(jì)算量大，耗時(shí)長(zhǎng)的問(wèn)題。GPU(Graphics Processing Units)由于其出色的加速效果和較低的硬件成本常常被大家采納。目前已有許多 GPU 加速的框架。TxBR 也采用了一種并行策略為曲線模型下的直接重構(gòu)算法進(jìn)行加速[3]。在此并行策略中，三維的重構(gòu)數(shù)據(jù)沿著z軸被分割為幾個(gè)部分(z-section)。每一個(gè)z-section 會(huì)被單獨(dú)地傳輸入 GPU 進(jìn)行重構(gòu)。然而這一并行策略并不適用于迭代算法的重構(gòu)。因?yàn)樵诿恳惠喌?，算法都需要將所有z-section 在GPU 與 CPU 傳輸一次。對(duì)于 GPU 來(lái)說(shuō)，這樣的數(shù)據(jù)傳輸十分耗時(shí)。目前的投影圖像的大小已經(jīng)達(dá)到了 8192×8192 像素，甚至更大[3,4]，相應(yīng)的重構(gòu)結(jié)果會(huì)達(dá)到幾個(gè) GByte[5]。由于 GPU 的顯存有限，這將進(jìn)一步使 GPU 并行變得困難。

針對(duì)這一問(wèn)題，我們提出了一種改進(jìn)的分塊迭代重構(gòu)算法，并且在 GPU 平臺(tái)上實(shí)現(xiàn)了算法的并行。我們的貢獻(xiàn)主要在兩個(gè)方面：首先，分析了曲線投影模型的局部性特征，并根據(jù)此特征將數(shù)據(jù)縱向劃分，從而提出了一種新的分塊迭代重構(gòu)算法；其次，針對(duì)大規(guī)模 ET 重構(gòu)數(shù)據(jù)，提出了一種基于頁(yè)的 GPU 數(shù)據(jù)傳輸策略，以減少數(shù)據(jù)傳輸時(shí)間，提高并行效率。

2 相關(guān)研究

2.1 迭代重構(gòu)算法

電子斷層重構(gòu)通過(guò)透射投影圖像獲得樣品內(nèi)部結(jié)構(gòu)。在實(shí)域空間中，迭代算法通過(guò)將此問(wèn)題形式化為線性方程組求解問(wèn)題。以“立體像素”為最小的三維空間的單位，我們可以將三維重構(gòu)結(jié)果表示為N個(gè)立體像素點(diǎn)(N＝nwidth×nlength×nheight，nlength、nwidth和nheight分別是樣品的長(zhǎng)、寬和高)。同時(shí)假設(shè)二維投影像素點(diǎn)的數(shù)目為M(M＝npro_width×npro_length×nang，其中npro_length、npro_width和nang分別是投影序列的長(zhǎng)、寬和角度數(shù))。投影的過(guò)程可以用公式(1)來(lái)表示：

其中，pi是第i個(gè)像素點(diǎn)的值；sj是第j個(gè)立體像素點(diǎn)的值。在矩陣A中，aij表示第j個(gè)立體像素點(diǎn)對(duì)第i個(gè)二維投影像素點(diǎn)的貢獻(xiàn)。矩陣A中各元素的值是由投影模型來(lái)決定的。 我們通過(guò)迭代算法來(lái)求解向量的值。

迭代重構(gòu)算法大致可分為順序迭代算法、分塊迭代算法和同步迭代算法三大類(lèi)[6]。在本質(zhì)上，順序迭代算法和同步迭代算法都是分塊迭代算法的特殊情況[7]。假設(shè)將所有線性方程組中的方程分為B組，每組方程的數(shù)目為T(mén)。我們可以獲得一種通用的迭代重構(gòu)算法公式：其中，b＝kmodB是分塊的序號(hào)；i是方程在線性方程組中的編號(hào)； 是權(quán)重因子。松弛因子λk子對(duì)于迭代算法的收斂速度有較大的影響[8]。

對(duì)于順序迭代算法(B＝M，T＝1)，每次更新數(shù)據(jù)只考慮使用一個(gè)約束。同步迭代算法(B＝1)每次更新需要考慮方程組中的所有約束，例如經(jīng)典的同步迭代算法 SIRT[9](Simultaneous Iterative Reconstruction Technique)。分塊迭代算法每次更新數(shù)據(jù)的時(shí)候只考慮一組方程，一輪迭代會(huì)更新數(shù)據(jù)B次。傳統(tǒng)的分塊迭代算法如BICAV[10]和 ASART[11]，都是采取逐角度更新的，每一個(gè)分塊中方程的數(shù)目為T(mén)＝npro_width×npro_length，分塊的數(shù)目為B＝nang。

2.2 曲線投影模型

在投影數(shù)據(jù)獲取的過(guò)程中，有很多因素會(huì)影響重構(gòu)結(jié)果的精度。首先，由于地球磁場(chǎng)的影響，電子的運(yùn)動(dòng)軌跡是曲線而不是通常認(rèn)為的直線；其次，在電子光學(xué)設(shè)備中，球面像差很明顯，使得投影中較小的細(xì)節(jié)模糊不清；再次，樣品在旋轉(zhuǎn)的過(guò)程中，不同部分在磁場(chǎng)的位置不同，這會(huì)導(dǎo)致樣品各部分成像的放大倍數(shù)不一致。同時(shí)樣品的扭曲與損失也會(huì)使重構(gòu)的結(jié)果變差。

面對(duì)上述這些問(wèn)題，為了獲得較好的重構(gòu)結(jié)果，TxBR[2]使用了一種捆綁調(diào)整算法求得一種非線性的投影模型。在本文中我們采用二次的曲線投影模型來(lái)描述投影圖像與三維樣品之間的非線性變換關(guān)系，公式如下：

該公式描述了三維數(shù)據(jù)到二維投影圖像的映射關(guān)系。(X,Y, Z)是立體像素的坐標(biāo)，其在角度下θ下對(duì)應(yīng)的投影坐標(biāo)是(x,y)。此公式中的系數(shù)都是由捆綁調(diào)整算法計(jì)算得到的，其細(xì)節(jié)在Lawrence 等[2]的研究中有詳細(xì)論述。

2.3 并行迭代算法

由于 GPU 的內(nèi)存有限，對(duì)于大尺度的電子斷層重構(gòu)問(wèn)題，我們無(wú)法將所有的數(shù)據(jù)直接放入 GPU中，因此需要考慮如何完成數(shù)據(jù)的分割。

線性投影模型(如圖 1(a)所示)假設(shè)電子束沿著直線傳播，所以電子軌跡是彼此平行的。在線性投影模型前提下，樣品中垂直于旋轉(zhuǎn)軸的一個(gè)切片總是被認(rèn)為會(huì)投影到一條直線上，從而三維重構(gòu)問(wèn)題可以分解為多個(gè)二維重構(gòu)問(wèn)題。使用前面所描述的算法，二維重構(gòu)問(wèn)題可以利用對(duì)應(yīng)的一維投影圖像求解。目前大多數(shù)的并行重構(gòu)算法都是采用這樣的策略[12-14]。但是這一策略并不適用于曲線投影模型，因?yàn)樵谇€投影模型(如圖1(b)所示)下，我們認(rèn)為電子束的運(yùn)動(dòng)軌跡是曲線，彼此之間并不平行，所以三維重構(gòu)問(wèn)題不能簡(jiǎn)單地分解為多個(gè)二維重構(gòu)問(wèn)題。

圖 1 直線投影模型與曲線投影模型Fig. 1 Straight-line model and curvilinear model

TxBR 提出了一種在曲線模型下針對(duì)直接重構(gòu)算法的 GPU 并行策略[3]。該策略將三維數(shù)據(jù)沿著z軸分割為幾個(gè)部分，每次獨(dú)立的重構(gòu)其中一個(gè)部分。然而這一并行策略需要大量的數(shù)據(jù)重復(fù)傳輸，并不適用于迭代重構(gòu)算法。到目前為止，還沒(méi)有針對(duì)曲線模型下迭代重構(gòu)算法的 GPU 并行策略。

3 技術(shù)細(xì)節(jié)

如前所述，針對(duì)曲線投影模型，目前沒(méi)有相關(guān)的 GPU 并行迭代重構(gòu)算法。相對(duì)于大尺度的重構(gòu)數(shù)據(jù)來(lái)說(shuō)，GPU 的內(nèi)存十分有限，所以找到一種適合的數(shù)據(jù)劃分策略是解決迭代重構(gòu)算法GPU 高效并行的關(guān)鍵問(wèn)題。為了解決這個(gè)問(wèn)題，我們首先研究了曲線模型的特點(diǎn)，并據(jù)此提出了一種數(shù)據(jù)分割的策略，之后相應(yīng)地修改了傳統(tǒng)的SIRT 迭代重構(gòu)算法。最后為了更加有效地減少數(shù)據(jù)傳輸?shù)南?，提出了一種基于頁(yè)策略的 GPU數(shù)據(jù)傳輸方式。

3.1 曲線投影模型的局部性

在直線投影模型下，我們認(rèn)為電子束的運(yùn)動(dòng)軌跡是直線，它們的軌跡是互相平行的，切片的投影始終是一條直線。但由于受到地球磁場(chǎng)等的影響，一個(gè)縱向切片的投影實(shí)際上是一條曲線。反過(guò)來(lái)看，投影到同一條直線上的立體像素點(diǎn)實(shí)際上是位于一個(gè)曲面上的。這個(gè)曲面會(huì)橫跨很多個(gè)縱向切片(如圖 1(b)所示)，而且曲面的形狀會(huì)隨著樣品旋轉(zhuǎn)角度的變化而變化。但是這個(gè)曲面是具有一定局部特性的。直觀地講，假設(shè)旋轉(zhuǎn)軸是y軸，y軸坐標(biāo)值小的投影點(diǎn)對(duì)應(yīng)的立體像素點(diǎn)的Y值也不會(huì)很大。

將公式(3)和(4)的左邊的值x和y固定，則這兩個(gè)式子描述了一條二次的三維空間曲線；如果只是固定x值(x＝c)，則式子(3)描述的是一個(gè)三維空間中的曲面(如圖 1(b)中展示的曲面)。根據(jù)此表達(dá)式，我們可以進(jìn)一步利用數(shù)學(xué)中關(guān)于曲面極值的相關(guān)知識(shí)來(lái)定量研究曲面的局部性。

將方程(3)看作一個(gè)隱函數(shù)，其中X是因變量而Y和Z是自變量。要定量地研究曲面的局部性，首先要將該問(wèn)題形式化地表示，即在一個(gè)有限的范圍內(nèi)求解曲面的極值點(diǎn)。該有限的范圍是樣品的空間范圍?；谒懻摰木唧w問(wèn)題，這個(gè)曲面是一個(gè)連續(xù)的曲面，其局部最大值和局部最小值可能是曲面的極值點(diǎn)，也可能是在邊緣曲線上。因此首先要求解曲面在給定范圍內(nèi)的極值點(diǎn)，再分析邊緣曲線的極值。這些點(diǎn)中的最大值和最小值就是我們需要的結(jié)果。

對(duì)于多元可微分函數(shù)，極值點(diǎn)首先是偏導(dǎo)數(shù)等于零的點(diǎn)，所以我們需要加入對(duì)于偏導(dǎo)數(shù)為零的約束條件，得到的方程如下：

接下來(lái)考慮邊緣數(shù)據(jù)的局部最大值和最小值，空間曲面與樣品邊界相交會(huì)出現(xiàn)這些邊界線。很明顯曲面一定會(huì)和樣品的上下兩個(gè)面相交，但并不會(huì)與樣品的所有面相交，隨著角度的變化，不同的曲面與樣品的交界并不一定相同。考慮所有情況，如求解曲面與平面Z＝Zmax交線的極值點(diǎn)需要用到公式：

邊界曲線的端點(diǎn)也可能是極值點(diǎn)，例如：

通過(guò)這幾步，我們可以獲得某一個(gè)曲面的局部最大值和最小值，亦即可得到投影點(diǎn)對(duì)應(yīng)的空間三維點(diǎn)所在的范圍。假設(shè)結(jié)果是，則顯然該結(jié)果會(huì)隨著角度的變化而變化，也會(huì)隨著c值的變化而變化。

3.2 縱向分塊迭代重構(gòu)算法

與同步迭代算法重構(gòu)不同，分塊迭代重構(gòu)算法每次更新數(shù)據(jù)只需要一部分約束條件，整個(gè)迭代過(guò)程是通過(guò)每次選取不同的約束方程子集來(lái)進(jìn)行的。目前將約束方程分塊的方法是基于角度的，例如 SIRT 是按照投影點(diǎn)的順序?qū)?duì)應(yīng)方程分塊。每一個(gè)塊的大小相同，是一個(gè)角度下所有投影點(diǎn)對(duì)應(yīng)的方程，塊的數(shù)目與投影角度數(shù)一致。

如前所述，由于 GPU 的內(nèi)存有限，所以對(duì)于大尺度三維重構(gòu)，需要考慮如何完成數(shù)據(jù)的分割。在直線模型下，三維重構(gòu)問(wèn)題可以分解為多個(gè)二維重構(gòu)問(wèn)題。這樣的數(shù)據(jù)分割不會(huì)引入多余的通信消耗。但是在曲線投影模型下，每次迭代都需要將所有的數(shù)據(jù)傳入 GPU 并對(duì)其進(jìn)行更新，這對(duì)于全局內(nèi)存小于重構(gòu)數(shù)據(jù)的情況就意味著，需要在一次迭代過(guò)程中將所有的數(shù)據(jù)分批次地在 GPU 與內(nèi)存之間傳輸，這將導(dǎo)致計(jì)算性能的急劇下降。在之前直線投影模型分割方式的啟發(fā)下，我們提出了一種新的基于 GPU 的分塊迭代重構(gòu)并行算法，將三維重構(gòu)結(jié)果利用曲線投影模型的局部性分割為橫跨幾個(gè)切片的縱向切塊，同時(shí)利用分塊迭代重構(gòu)算法每次只更新需要約束方程子集特點(diǎn)的數(shù)據(jù)。

首先將投影圖像序列垂直于旋轉(zhuǎn)軸均勻地分割為B個(gè)部分。對(duì)于給定的分割數(shù)目B＝b，利用曲線投影模型的局部性可以估計(jì)出每個(gè)部分的投影點(diǎn)集合St：

其次，利用曲線投影模型的空間局部性，可以根據(jù)投影點(diǎn)的范圍計(jì)算得到對(duì)應(yīng)三維空間點(diǎn)的范圍，再利用坐標(biāo)信息通過(guò)代碼完成數(shù)據(jù)分割。這些投影點(diǎn)方程中的非零點(diǎn)(X,Y,Z)一定在此范圍內(nèi)。根據(jù)求出的X的范圍，可將三維重構(gòu)數(shù)據(jù)分割為幾個(gè)縱向的切塊。每次更新數(shù)據(jù)，只需要將一個(gè)切塊傳輸入 GPU 就可以。這樣極大地減少了數(shù)據(jù)傳輸對(duì)于 GPU 并行的影響。

圖 2 改進(jìn)的分塊迭代重構(gòu)算法Fig. 2 Modified block iterative method

綜上所述，首先對(duì)所有投影點(diǎn)(x,y,g)進(jìn)行整體分割，運(yùn)用之前所述的曲面求極值的方法求出它們所對(duì)應(yīng)的所有三維空間點(diǎn)的范圍。接著按照這一計(jì)算結(jié)果，對(duì)三維空間點(diǎn)(X,Y,Z)進(jìn)行縱向的分割，如圖 2 所示，每次更新過(guò)程中，只將一部分投影點(diǎn)和它對(duì)應(yīng)的三維空間點(diǎn)傳入 GPU 以完成算法的更新。在一次迭代結(jié)束之后，再將下一塊數(shù)據(jù)放入 GPU。

其中，分割的數(shù)目B受限于 GPU 全局內(nèi)存的大小，用sizeslab表示最大的縱向切塊需要的內(nèi)存大小，sizepro表示投影圖像的大小，sizeinter表示其他中間變量占的內(nèi)存大小。同一時(shí)刻需要的內(nèi)存大小可以表示為：

3.3 基于頁(yè)的數(shù)據(jù)傳輸策略

由于相鄰的縱向切塊之間有重復(fù)的部分，如果每次重構(gòu)一個(gè)縱向切塊后就將其傳輸出去，重復(fù)的數(shù)據(jù)會(huì)被重復(fù)傳輸。我們提出了一種基于分頁(yè)的數(shù)據(jù)傳輸策略消除冗余的數(shù)據(jù)傳輸。在該策略下，相鄰縱向切塊之間的重疊數(shù)據(jù)不會(huì)在前一個(gè)切塊重構(gòu)之后傳輸回 CPU，而是留在 GPU 中等待第二個(gè)縱向切塊的其他數(shù)據(jù)傳入后，一起完成第二個(gè)縱向切塊的重構(gòu)。這樣，每次只將不再需要的數(shù)據(jù)傳輸回 CPU。

這種策略的最基本思想是延遲拷貝，只有不需要的數(shù)據(jù)才會(huì)傳輸出去。但是如果每次只將不用的數(shù)據(jù)傳輸，然后將有用的數(shù)據(jù)傳入 GPU，將會(huì)導(dǎo)致傳出的數(shù)據(jù)比新傳入的數(shù)據(jù)量小，使得傳入的數(shù)據(jù)沒(méi)有足夠的空間進(jìn)行存儲(chǔ)。借鑒內(nèi)存管理的思想，我們將 GPU 的全局內(nèi)存分割成相同大小的頁(yè)，并使其作為每次數(shù)據(jù)傳輸?shù)幕締挝?。以?shí)際數(shù)據(jù)為例，使用投影大小為 512×512的數(shù)據(jù)，則整個(gè)投影圖像的序列有 121 張。其重構(gòu)的結(jié)構(gòu)的尺寸是 650×683×101，空間坐標(biāo)分布為{(X,Y,Z)|X∈[－87, 563],Y∈[－31, 652],Z∈[－62, 39]}。將投影圖像沿著垂直于旋轉(zhuǎn)軸的方向，把投影圖像分給為 8 份，對(duì)應(yīng)的樣品空間點(diǎn)也會(huì)沿著垂直于旋轉(zhuǎn)軸的方向分割為 8份，其對(duì)應(yīng)的縱向切塊的分布如表 1。

表 1 中，“Slab”表示切塊的序號(hào)；“Start”和“End”分別對(duì)應(yīng)于每一個(gè)分塊的起始位置和終止位置的X值；“Out”對(duì)應(yīng)的數(shù)據(jù)是在將本分塊更新完成后需要傳出數(shù)據(jù)的大小，它是以X值所涉及的范圍表示的，比如第一個(gè)分塊在迭代完后，需要傳出的實(shí)際數(shù)據(jù)量是(－12－(－78)＋1)×683×101＝67×683×101；“In”和“Range”都是直接用X的范圍來(lái)表示數(shù)據(jù)量的大小，其中“In”表示了本分塊在迭代開(kāi)始之前需要傳輸進(jìn)去的數(shù)據(jù)量大??；“Range”表示的是一個(gè)分塊所涉及的數(shù)據(jù)量的大小。

由表 1 可看出，相鄰分塊之間重疊的部分達(dá)到了整個(gè)分塊的 3/4，所以去除冗余的數(shù)據(jù)傳輸很有必要。此外，每次傳入和傳出的數(shù)據(jù)量相差不大，這就保證了該策略可以有效地利用全局內(nèi)存。

表 1 縱向分塊的情況Table 1 The distribution of slabs

對(duì)數(shù)據(jù)按照傳輸?shù)倪^(guò)程重新進(jìn)行分割，數(shù)據(jù)自動(dòng)分頁(yè)的算法如下：首先，以每個(gè)塊的開(kāi)始位置作為頁(yè)的起始處，保證每次傳輸出去的數(shù)據(jù)都是下一次迭代中不需要的。在保證頁(yè)的起始位置正確的前提下，盡可能多地傳入數(shù)據(jù)以達(dá)到充分利用全局內(nèi)存的作用。

這里的旋轉(zhuǎn)軸是x軸，將生物樣品的三維密度點(diǎn)按照求出的X值進(jìn)行分割，得到的數(shù)據(jù)分布如表 2。

表 2 數(shù)據(jù)的頁(yè)分布Table 2 The distribution of pages

在每一輪迭代的開(kāi)始，都先將第一次更新需要的數(shù)據(jù)(Page 1 to 5)全部傳入 GPU。在接下來(lái)的分塊迭代中，將不再需要的一頁(yè)數(shù)據(jù)(Page 1)傳出，并將新的一頁(yè)數(shù)據(jù)(Page 6)傳入。在樣品的最后一頁(yè)數(shù)據(jù)(Page 10)傳入 GPU 的同時(shí)停止數(shù)據(jù)的傳入和傳出。在本輪的迭代結(jié)束以后，可以開(kāi)始新的一輪的迭代。

4 實(shí) 驗(yàn)

在實(shí)驗(yàn)部分，我們首先對(duì)分塊迭代的 SIRT算法與同步迭代的 SIRT 算法的重構(gòu)結(jié)果進(jìn)行了比較；其次，討論了并行算法的時(shí)間消耗。

本實(shí)驗(yàn)使用的生物樣品厚度為 350 nm，使用放大倍數(shù)為 37 K 的 300 kV FEI Titan TEM透射顯微鏡拍攝透射投影圖像，投影圖像序列有 121 張，角度范圍為(－60°，＋60°)，投影圖像的大小為 512×512。重構(gòu)結(jié)果的大小為651×684×101。實(shí)驗(yàn)中的數(shù)據(jù)分割方法在“基于頁(yè)的數(shù)據(jù)傳輸策略”中已討論。

所有實(shí)驗(yàn)均在 64 位 Ubuntu12.04 操作系統(tǒng)上實(shí)現(xiàn)，使用的 GPU 顯卡是 NVIDIA GTX480，具有 480 SPs 和 1536 M 全局內(nèi)存。

4.1 重構(gòu)結(jié)果

所有實(shí)驗(yàn)均采用二次曲線投影模型。在本實(shí)驗(yàn)中，我們考慮了三種情況：(1)使用傳統(tǒng)的同步迭代重構(gòu)算法 SIRT；(2)將數(shù)據(jù)分割為 4 個(gè)切塊，使用分塊迭代版本的 SIRT 算法更新數(shù)據(jù)；(3)將數(shù)據(jù)分割為 8 個(gè)切塊，使用分塊迭代版本的 SIRT 算法更新數(shù)據(jù)。當(dāng)分塊的數(shù)目為 1 時(shí)，分塊迭代重構(gòu)算法就變成了同步迭代重構(gòu)算法。實(shí)驗(yàn)中主要考慮兩個(gè)參數(shù)：松弛因子和迭代次數(shù)。這里將所有投影點(diǎn)對(duì)應(yīng)的方程都定義為一輪迭代。因?yàn)楫?dāng)松弛因子 0＜λ＜2 時(shí)，SIRT 算法保證收斂[15]，本實(shí)驗(yàn)中松弛因子λ分別取值為1、0.5、0.2 和 0.1。

在本實(shí)驗(yàn)中，我們使用投影誤差作為評(píng)判結(jié)果好壞的標(biāo)準(zhǔn)。投影誤差反應(yīng)了再投影與原始投影圖像之間的誤差。這里使用絕對(duì)平均誤差(MAE)作為投影誤差的依據(jù)。

圖 3 絕對(duì)平均誤差Fig. 3 Mean absolute projection error

其中，pi是實(shí)驗(yàn)獲取的投影原始圖像；是利用重構(gòu)結(jié)果求出的再投影結(jié)果。圖 3 中的曲線顯示了絕對(duì)平均誤差的結(jié)果?？梢钥闯?，對(duì)于不同的松弛因子，分塊迭代算法都比同步迭代算法在前20 輪收斂的速度更快。這說(shuō)明分塊迭代算法的收斂速度更快。將數(shù)據(jù)分割為 8 個(gè)切塊的分塊迭代算法比分為 4 塊的分塊迭代算法的投影誤差更小，結(jié)果精度更高。通過(guò)松弛因子λ的不同取值可以發(fā)現(xiàn)，本測(cè)試數(shù)據(jù)，比較大的松弛因子可以獲得更好的重構(gòu)結(jié)果。

圖 4 展示了 50 輪迭代之后同步迭代算法與分塊迭代算法的重構(gòu)結(jié)果。圖 4(a)為同步迭代算法 SIRT 的重構(gòu)結(jié)果，圖 4(b)是分塊迭代版本的SIRT 算法的重構(gòu)結(jié)果。兩者的結(jié)果差距不大，但同步迭代算法在一些細(xì)節(jié)上更加清晰。

圖 4 重構(gòu)結(jié)果Fig. 4 Reconstruction result

4.2 并行算法的效果

我們?cè)?GTX480 上測(cè)試了分塊迭代重構(gòu)算法的運(yùn)行時(shí)間。具體的運(yùn)行時(shí)間和加速比如表 3 所示。表 3 中，B表示分塊的數(shù)目(B＝4 表示將所有投影點(diǎn)分為 4 組，相應(yīng)地，三維結(jié)果被分為了4 個(gè)縱向切塊)。

可以看到運(yùn)行時(shí)間基本與迭代的次數(shù)成正比?？傮w來(lái)說(shuō)，分塊數(shù)目為 4 的迭代算法的運(yùn)行時(shí)間少于分塊為 8 的迭代算法的運(yùn)行時(shí)間。這是因?yàn)榉謮K數(shù)目越多，并行處理的數(shù)據(jù)越少，算法的并行度會(huì)降低。所以在硬件允許的前提下，應(yīng)該盡量減少分塊的數(shù)量。

表 3 算法執(zhí)行時(shí)間與加速比Table 3 Running time and speedup

5 結(jié)論與討論

曲線投影模型能夠提高大尺度電子斷層重構(gòu)的精度，但它的計(jì)算復(fù)雜高，限制了曲線投影模型的廣泛應(yīng)用。在電子斷層三維重構(gòu)領(lǐng)域，有許多關(guān)于三維重構(gòu)算法的并行加速工作。絕大部分的工作是在直線模型下對(duì)迭代重構(gòu)算法采用 MPI或 GPU 并行。 目前針對(duì)曲線投影模型的工作并不多。已有的并行加速工作主要是針對(duì)直接背投影算法。我們采用 GPU 對(duì)于更加耗時(shí)的迭代重構(gòu)算法進(jìn)行加速， 使得在曲線模型下的迭代算法也可以快速地獲得結(jié)果。

本文提出了曲線模型下，針對(duì)迭代重構(gòu)算法的 GPU 并行算法。在我們的工作中，我們通過(guò)對(duì)曲線模型的研究發(fā)現(xiàn)曲線模型具有一定的空間局域性，利用這種性質(zhì)我們提出了一種縱向的分塊方式。在算法的實(shí)現(xiàn)階段，我們提出一個(gè)基于頁(yè)的數(shù)據(jù)傳輸策略，從而能夠去除冗余的數(shù)據(jù)傳輸。實(shí)驗(yàn)結(jié)果顯示，本算法可以接近 40 倍的加速比。在以后的工作中，我們會(huì)將本文提出的并行策略應(yīng)用到其他的并行平臺(tái)上，并采用更多的實(shí)驗(yàn)數(shù)據(jù)來(lái)驗(yàn)證我們的算法。

[1] Phan S, Lawrence A, Molina T，et al. Txbr montage reconstruction for large field electron tomography [J]. Journal of Structural Biology, 2012, 180: 154-164.

[2] Lawrence A, Bouwer JC, Perkins G, et al. Transform-based backprojection for volume reconstruction of large format electron microscope tilt series [J]. Journal of Structural Biology, 2006, 154(2): 144-167.

[3] Lawrence A, Phan S, Singh R. Parallel processing and large-field electron microscope tomography [C] // World Congress on Computer Science and Information Engineering, 2009: 339-343.

[4] Wan X, Phan S, Lawrence A, et al. Iterative methods in large field electron microscope tomography [J]. SIAM Journal on Scientific Computing, 2013, 35(5): S402-S419.

[5] Agulleiro JI, Fernández JJ. Evaluation of a multicore optimized implementation for tomographic reconstruction [J]. PloS One, 2012, 7(11): e48261.

[6] Censor Y. Parallel optimization: theory, algorithms and applications [D]. Oxford: Oxford University Press, 1997: 127-190.

[7] Aharoni R, Censor Y. Block-iterative projection methods for parallel computation of solutions to convex feasibility problems [J]. Linear Algebra and Its Applications, 1989, 120: 165-175.

[8] Herman GT, Meyer LB. Algebraic reconstruction techniques can be made computationally efficient [J]. IEEE Transactions on Medical Imaging, 1993, 12(3): 600-609.

[9] Gilbert P. Iterative methods for the threedimensional reconstruction of an object from projections [J]. Journal of Theoretical Biology, 1972, 36(1): 105-117.

[10] Censor Y, Gordon D, Gordon R. BICAV: A blockiterative parallel algorithm for sparse systems with pixel-related weighting [J]. IEEE Transactions on Medical Imaging, 2001, 20(10): 1050-1060.

[11] Wan X, Zhang F, Chu Q, et al. Three dimensional reconstruction using an adaptive simultaneous algebraic reconstruction technique in electron tomography [J]. Journal of Structural Biology, 2011, 175(3): 277-287.

[12] Wan X, Zhang F, Chu Q, et al. High-performance blob-based iterative reconstruction of electron tomography on multi-gpus [C] // Bioinformatics Research and Applications, the 7th International Symposiu, ISBRA 2011, 2011: 61-72.

[13] Fernández JJ. High performance computing in structural determination by electron cryomicroscopy [J]. Journal of Structural Biology, 2008, 164(1): 1-6.

[14] Casta?o Díez D, Mueller H, Frangakis AS. Implementation and performance evaluation of reconstruction algorithms on graphics processors [J]. Journal of Structural Biology, 2007, 157(1): 288-295.

[15] van der Sluis A, van der Vorst HA. SIRT-and CG-type methods for the iterative solution of sparse linear least-squares problems [J]. Linear Algebra and its Applications, 1990, 130: 257-303.

A Parallel Reconstruction Algorithm Based on Curvilinear Projection Model in Tomography

ZHANG Jingrong1,2WAN Xiaohua1ZHANG Fa1

1(Institute of Computing Technology,Chinese Academy of Sciences,Beijing100190,China)
2(University of Chinese Academy of Sciences,Beijing100049,China)

Large-field high-resolution electron tomography enables visualizing detailed mechanisms under global structure. As field enlarges, the distortions of reconstruction and processing time become more critical. TxBR has proposed a curvilinear projection model, which can dramatically improve the quality of reconstruction. But its computation is more complex and time-consuming. Furthermore, previous parallel strategies are not suitable for curvilinear projection model. In this work, a block iterative parallel algorithm using curvilinear projection model on GPU platform was proposed. By studying the locality of curvilinear projection model, we proposed a vertical data decomposition method. We also adopt a page-based data transfer scheme to reduce the processing time. Experimental results show that our method can yield speedups of approximate 40 times.

electron tomography; iterative reconstruction methods; block; curvilinear projection model; GPU parallel

Q 617

A

2014-01-15

：2015-03-17

張靜蓉，博士研究生，研究方向?yàn)樯镄畔⑴c GPU 并行；萬(wàn)曉華，博士，助理研究員，研究方向?yàn)樯镄畔?、生物醫(yī)學(xué)圖像處理與高性能計(jì)算；張法(通訊作者)，博士，副研究員，研究方向?yàn)樯镄畔⑴c GPU 并行，E-mail：zhangfa@ict.ac.cn。