Mojette 變換層析技術中的投影角度布局方法*

吳慎將 劉榮明 王佳 李黨娟 程軍霞

(西安工業大學光電工程學院,西安 710021)

Mojette 變換是一種最小冗余采樣的離散Radon 變換,能夠用較少角度的投影數據進行精確的計算層析(computed tomography,CT)重建,為少量投影角度CT 技術的實現提供了一種新思路.投影角度的空間布局決定了層析重建最少所需投影的數量.為了獲得Mojette 變換層析技術中的最優投影空間角度布局方案,本文對三維Mojette 變換數學模型及其精確重建條件進行了研究.以此為基礎,在考慮實際探測器像素數目受限的條件下,提出了確定最優投影角度的方法.研究結果表明: 所有探測器圍繞被測物體在同一水平面內進行平行投影采集是最優的投影角度布局方案,此時投影模型為二維Mojette 變換,所需的投影角度和探測器像素數最少,投影角度范圍最小; 若在實際的測量中該投影條件無法滿足,則投影矢量中|pi|和|qi|的值越小越好.該研究可為實際層析系統的建立提供理論基礎.

1 引 言

計算層析(computed tomography,CT)技術,是一種由低維投影數據重建高維目標的技術,已廣泛應用在Terahertz 波檢測[1]、量子態[2]、醫學三維成像[3]、風洞[4]、地質探測、激光打靶和燃燒場三維成像和檢測等領域中[5?7].它利用探測器采集測試目標在多個角度的二維投影,并結合層析理論進行三維重建.在實際應用中,由于探測對象或測量環境的影響,經常會遇到投影采集角度受限制的問題[8].如何利用有限角度投影進行精確層析重建,對于層析成像技術的發展和應用,具有非常重要的意義.

傳統的基于Radon 變換的層析技術很難在極少角度采樣的情況下獲得較好的重建結果.Mojette 變換是一種最小冗余采樣的離散Radon 變換,可以根據多個投影之間的相互獨立特性對投影個數、投影角度等進行變化,通過改變不同投影矢量下的采樣率來控制冗余度的大小,因此可以在最大程度上避免投影信息的重復和冗余采樣[9].Mojette變換的理論基礎由Katz[10]提出的離散角度概念以及Herman[11]提出的迭代算子共同構筑而成,該變換利用滿足Katz 引理的稀疏角度即可被精確重建,其重建所需的數據采集量遠小于Radon 變換所需的數據量[12?14].基于Mojette 變換的層析重建能夠顯著減少所需的投影角度和投影射線條數,其重建所需的數據采集量遠小于Radon 變換所需的數據量,在稀疏角度下具有良好的重建性能.并且實際的Radon 變換投影可以轉換為Mojette 投影,為基于Mojette 變換的實際投影層析重建提供可能[15,16].

Mojette 變換層析理論中的可精確重建條件以及最少投影角度布局對于實際層析系統的建立以及提高重建精度具有非常重要的指導意義[17].在傳統的基于二維Radon 變換[18,19]和二維Mojette變換的層析重建技術[20,21]中,探測器放置在被測物體周圍同一水平面內,進行平行投影的采集.當實際測量環境中水平面內投影角度受限時,可以在三維空間中進行投影采集,此時層析投影模型為三維Radon 變換或三維Mojette 變換.Cai 等[22]利用數值計算的方法討論了基于三維Radon 變換層析技術中投影角度的三維空間分布對重建精度的影響.目前,從理論上分析和解釋層析系統中投影角度的最優布局方案,特別是對三維Mojette 變換以及相應的層析重建理論的相關研究仍較少.

為了在理論上獲得Mojette 變換層析技術中的最優投影角度空間布局方案,本文將建立三維Mojette 變換數學模型,并且利用基于角的重建(corner based inversion,CBI)算法對精確重建條件進行研究.以此為基礎,結合實際探測器像素數目受限條件,提出并確定最優投影角度方案.

2 三維Mojette 變換數學模型

在三維直角坐標系 (x,y,z) 中,三維Mojette變換的投影方向用三維離散向量ξi=(pi,qi,ri) 來表示,其中pi ∈Z,qi ∈Z+和ri ∈Z分別表示投影向量在x,y和z軸的分量,投影方向限制在y軸正向,并且i=1,2,··· ,N表示投影角度數.如圖1所示,投影 向量ξi=(pi,qi,ri) 對應 的投影角度由方位角φi和天頂角θi確定,投影探測平面垂直于投影向量.將被測三維物體f(x,y,z) 均勻劃分為分辨率為P ×Q×R的離散網格f(k,l,m) .

圖1 三維Mojette 變換示意圖Fig.1.Schematic diagram of three-dimensional Mojette transform.

當pi=0 ,或qi=0 ,或ri=0 時,三維Mojette變換等效為相應方向上的二維Mojette 變換.以ri=0為例,三維Mojette 變換等效為m組水平方向上的二維Mojette 變換,其投影方程為

其中po1 為一個修正值,使投影像素的序號從1開 始,當pi >0 時,po1=1 ,當pi <0 時,po1=?(Q ?1)·pi+1 .每行投影像素數目B1和相鄰像素間隔h1分別為

當pi ′=0,qi ′=0 ,并且ri ′=0 時,如圖1 所示,三維Mojette 變換可以看成一組在平行于平面A的平面內的二維Mojette 變換結果.探測平面上每行像素的Mojette 變換投影值為在該行像素所對應的平面A內投影射線經過中心的所有網格數值的積分,投影的行數由平行于平面A的平面數決定.因此,三維Mojette 變換可以分解為兩次二維Mojette 變換:

步驟1矢量 (pi,qi) 確定了 (x,y) 平面內的投影方向,沿該方向的二維Mojette 變換決定了最終三維Mojette 變換投影的行數,行數和相鄰行之間的間隔由(2)式確定.根據二維Mojette 變換的要求,pi和qi互質,即 G CD(pi,qi)=1 .

步驟2在矢量 (pi,qi) 和ri構成的二維平面A內,以ξi為投影方向進行二維Mojette 變換,其結果決定了三維Mojette 變換投影的列數和最終的投影值.

由于第1 次二維Mojette 變換投影矢量(pi,qi)的長度為,并且其投影間隔為h1=,因 此 可 得 (pi,qi) 方 向 的 投 影 整 數 為p2i+qi2.平面A內水平方向的網格數為 (x,y) 平面內的離散網格沿 (?qi,pi) 方向進行二維Mojette變換后的投影數,可以表示為b′=(P ?1)|pi|+(Q ?1)|qi|+1,垂直方向的網格數為R.在平面A內以投影方向 (p2i+qi2,ri) 進行二維Mojette 變換,要求GCD(p2i+qi2,ri)=1 .

三維Mojette 變換投影的列數B2和相鄰投影間隔h2分別為

根據以上分析,三維Mojette 變換可表示為

其 中n=?qi(l ?1)+pi(k ?1)+po2 ,po1,po2 和po3 為投影像素序號的修正值,取值分別為

利用(4)式和(5)式計算一個P=Q=R=10的全1 矩陣在投影矢量 (pi,qi,ri)=(1,1,3) 和(pi,qi,ri)=(1,2,3)時 的 三 維Mojette 變 換,投影歸一化結果如圖2(a)所示.圖2(b)為全1 矩陣在該投影角度下的視覺投影顯示結果.可以看出,兩種投影結果在像素數目、投影分布上完全相同,驗證了該三維Mojette 變換數學模型的正確性.

3 精確重建條件

1978 年Katz[10]給出一個約束投影角度數量上界的公式,即著名的Katz 引理,該定理指出對于一簇互質的投影矢量對 (pi,qi) ,如果重建圖像的分辨率P×Q滿足關系

則最少可通過N個投影角度即可完成精確重建.(6)式為二維Mojette 變換的精確重建條件.根據不同重建圖像分辨率P×Q和投影角度數的要求,可以選擇合適的投影矢量(pi,qi) 完成精確重建.

從幾何角度出發,二維Mojette 變換的精確重建條件可理解為: 若所有投影矢量的絕對和超出被測區域,則該圖像可被精確重建.將該結論推廣至三維情況,可得三維Mojette 變換的精確重建條件為

為了驗證該精確重建條件的正確性,利用CBI 算法進行層析重建[10,12].CBI 算法是Mojette變換最基本的重建理論.與各種迭代類、變換類算法不同,該算法是一種精確的求解線性方程組的算法,二維圖像中所有像素點的值被依次精確重建出來.CBI 算法能夠準確驗證線性方程組的可解性,其重建結果可以準確說明Mojette 變換的精確重建條件是否成立,并且通過CBI 驗證的精確重建條件對共軛梯度算法、反投影算法等均適用[11,14].

利用CBI 算法進行三維Mojette 變換重建時,除了對待重建三維物體進行正常投影外,在相同投影條件下對與重建物體相同維度的三維全1 矩陣和索引矩陣進行投影.全1 矩陣投影的作用是: 通過投影可以從投影向量值中看出某投影矢量下的射線穿過的像素中心點的個數,當投影向量值中的分量為1 時,說明該投影射線穿過的路徑上只有一個像素點值的貢獻,可直接求解出變量值.當向量值中有多個1 存在時,可依次進行重建.索引矩陣投影的作用是: 索引矩陣中的像素點值從左上角的(1,1,1)點賦值為0 開始,從左向右、從前到后、從上至下像素點以此遞增,每次加1,遍歷至右下角最后一個像素點時,該點的值為P×Q×R ?1 .設置索引矩陣的目的在于,在全1 矩陣的投影向量值中找到投影值為1 的分量時,可以讓計算機理解該投影對應的像素點的位置.其逆變換的求解過程是一種串行求解模式,即每次迭代求出全1 矩陣投影值為1 所對應的一批像素點,再根據這些像素點求出的值來更新被測物體、全1 矩陣、索引矩陣投影,從而產生出投影值即為點值的像素點,重復這一步驟,直至所有的待求點與探測器像元上更新完畢后的一個投影值一一對應,迭代結束.

選擇經典的Shepp-Logan 模型進行數值模擬實驗以驗證精確重建條件的正確性.模擬三維物體的分辨率為P=Q=R= 48,所有水平二維平面為相同的48×48 的Shepp-Logan 模型分布,如圖3(a)所示.對表1 中列出的不同投影矢量條件下的三維Mojette 變換投影及重建進行數值仿真.

圖2 (a)三維Mojette 變換投影; (b)視覺投影Fig.2.(a) Three-dimensional Mojette transform projection; (b) visual projection.

圖3(b)—(d)分別為Case 1,Case 3 和Case 5三種滿足精確重建條件時的重建結果.可以看出,當選取的投影矢量滿足精確重建條件時,都可以進行精確的層析重建.在相同的投影角度數下,即使一個投影矢量發生變化使得精確重建條件不滿足,就得不到正確的重建結果.數值模擬實驗驗證了精確重建條件的正確性.

4 最優的投影角度布局

由上述三維Mojette 變換投影模型可知,投影向量的選擇決定了投影采集所需探測器的像素數目和像素大小.對于三維Mojette 變換在水平和垂直方向的投影間隔不同的問題,在實際的投影探測中可以使用定制的成像鏡頭,使水平和垂直方向的放大率不同,可實現行和列上相鄰兩個像素的間隔不同.因此,本文只討論探測器像素數目的影響.由于實際探測器的像素個數有限,則滿足精確重建條件的實際投影向量和投影角度數應符合以下條件:

圖3 滿足精確重建條件下不同投影矢量的重建結果 (a) 模擬物體; (b) Case 1; (c) Case 3; (d) Case 5Fig.3.Reconstruction results of different projection vectors under accurate reconstruction condition: (a) Simulative object; (b) Case 1;(c) Case 3; (d) Case 5.

表1 模擬仿真中采用的不同投影矢量Table 1.Projection vectors used in simulation.

其中 (B1,B2) 為實際探測器的像素數.

確定最優投影向量的具體步驟如下.

步驟1確定滿足(8a)式、(8b)式和(8e)式的所有投影矢量 (pi,qi,ri) .

步驟2根據(8c)式,要滿足精確重建條件有三種方案,即或與其對應的條件是選取的投影矢量中pi,qi或ri的絕對值越大,投影角度數越少.且要保證投影像素數最少,則在投影矢量的一個分量保持較大值時,其他兩個分量的值盡可能小.因此,將步驟1 所確定的所有投影矢量對pi,qi或ri的絕對值進行降序排序,從大到小選取投影矢量,直至滿足(8c)式,則可確定投影角度數和其對應的投影矢量.

步驟3選擇步驟2 中確定的三種方案投影角度數的最小值,其投影矢量即為最優的投影角度布局方案.或者在實際的測量系統中,根據測試條件對投影角度的限制,選擇三種方案中最好實現的一種為最優投影角度布局方案.

圖4 不同投影矢量對應的投影角度空間分布 (a)方案1; (b)方案2; (c)方案3; (d)方案4; (e)方案5; (f)方案6Fig.4.Spatial distribution of projection angles corresponding to different projection vectors: (a) Scheme 1; (b) Scheme 2;(c) Scheme 3; (d) Scheme 4; (e) Scheme 5; (f) Scheme 6.

以P=Q=R=64 ,像 素 數 為1024×2048的探測器為例對最優投影角度的布局進行說明.當|ri|=0時,最少需要5 個投影角度可實現精確重建.若要求投影像素數最少,則選取的投影矢量為[(15,1,0),(–15,1,0),(14,1,0),(–14,1,0),(13,1,0)](方案1)或[(1,15,0),(–1,15,0),(1,14,0),(–1,14,0),(1,13,0)] (方案2),其空間 分布如圖4(a)和圖4(b)所示,該條件下探測器限制在較小的角度范圍內.該投影矢量方向的二維Mojette變換投影結果如圖5(a)和圖5(b)所示,投影像素數分別為64×1009,64×1009,64×946,64 ×946 和64×883.若探測器角度范圍不受限制,選取投影矢量為[(15,1,0),(–15,2,0),(14,5,0),(–14,9,0),(13,11,0)] (方 案3)或[(1,15,0),(–2,15,0),(5,14,0),(–9,14,0),(11,13,0)] (方案4),空間分布如圖4(c)和圖4(d)所示.Mojette變換投影如圖5(c)和圖5(d)所示,投影像素數分別為64×1009,64×1072,64×1198,64×1450和64×1513.

若 選取 的投 影矢 量中 | pi|>1 或 | qi|>1 ,此時對應的 | ri| 變小,則所需最少投影角度數會增多.例如滿足精確重建條件的另一組投影矢量為[(1,2,9),(–1,2,9),(1,2,–9),(–1,2,–9),(2,1,9),(–2,1,9),(2,1,–9),(–2,1,–9)] (方案6),空間分布如圖4(f)所示.對應的投影結果如圖5(f)所示.

圖5 不同投影矢量對應的投影結果 (a)方案1; (b)方案2; (c)方案3; (d)方案4; (e)方案5; (f)方案6Fig.5.Projections corresponding to different projection vectors: (a) Scheme 1; (b) Scheme 2; (c) Scheme 3; (d) Scheme 4;(e) Scheme; (f) Scheme 6.

圖4和圖5 的結果綜合表明: 在探測器像素數受限的條件下,最優的投影角度布局方案為水平面投影,即基于二維Mojette 變換的層析重建,此時所需投影角的數目最少.并且不論投影角度是否受限,其所需的投影探測器的像素數都比三維布局方案的少.這種方式與傳統層析技術中所選取的探測器的空間布局方案完全一致.當投影矢量無法滿足水平面投影時,要選擇為精確重建條件,并且選擇|pi|和|qi|的值越小,所需的投影角度數和探測器的像素數越少.

5 結 論

為了解決計算層析技術中投影采集角度受限制的問題,利用有限角度的投影實現高精度的層析重建,本文在建立三維Mojette 變換數學模型及其精確重建條件的基礎上,對Mojette 變換層析技術中的最優投影空間角度布局方案進行了研究.在綜合考慮精確重建條件和實際探測器像素數目受限的條件下,提出了確定最優投影角度的方法.研究結果表明: 1)若要求層析采集系統中投影角度數和投影像素數盡可能少,則探測器要分布在被測物體周圍的同一水平面內進行平行投影的采集,此時層析模型為二維Mojette 變換及重建; 2)當投影條件受限,無法實現水平面投影采集時,則投影矢量中|pi|和|qi|的值越小越好.該投影角度布局方案與傳統的層析系統中探測器的空間布局方案完全一致,本文首次從理論上說明了這種布局方案的優越性.