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

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

(西安工業(yè)大學(xué)光電工程學(xué)院,西安 710021)

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

1 引 言

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

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

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

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

2 三維Mojette 變換數(shù)學(xué)模型

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

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

當(dāng)pi=0 ,或qi=0 ,或ri=0 時(shí),三維Mojette變換等效為相應(yīng)方向上的二維Mojette 變換.以ri=0為例,三維Mojette 變換等效為m組水平方向上的二維Mojette 變換,其投影方程為

其中po1 為一個(gè)修正值,使投影像素的序號(hào)從1開 始,當(dāng)pi >0 時(shí),po1=1 ,當(dāng)pi <0 時(shí),po1=?(Q ?1)·pi+1 .每行投影像素?cái)?shù)目B1和相鄰像素間隔h1分別為

當(dāng)pi ′=0,qi ′=0 ,并且ri ′=0 時(shí),如圖1 所示,三維Mojette 變換可以看成一組在平行于平面A的平面內(nèi)的二維Mojette 變換結(jié)果.探測(cè)平面上每行像素的Mojette 變換投影值為在該行像素所對(duì)應(yīng)的平面A內(nèi)投影射線經(jīng)過中心的所有網(wǎng)格數(shù)值的積分,投影的行數(shù)由平行于平面A的平面數(shù)決定.因此,三維Mojette 變換可以分解為兩次二維Mojette 變換:

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

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

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

三維Mojette 變換投影的列數(shù)B2和相鄰?fù)队伴g隔h2分別為

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

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

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

3 精確重建條件

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

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

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

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

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

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

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

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

4 最優(yōu)的投影角度布局

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

圖3 滿足精確重建條件下不同投影矢量的重建結(jié)果 (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) 為實(shí)際探測(cè)器的像素?cái)?shù).

確定最優(yōu)投影向量的具體步驟如下.

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

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

步驟3選擇步驟2 中確定的三種方案投影角度數(shù)的最小值,其投影矢量即為最優(yōu)的投影角度布局方案.或者在實(shí)際的測(cè)量系統(tǒng)中,根據(jù)測(cè)試條件對(duì)投影角度的限制,選擇三種方案中最好實(shí)現(xiàn)的一種為最優(yōu)投影角度布局方案.

圖4 不同投影矢量對(duì)應(yīng)的投影角度空間分布 (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 ,像 素 數(shù) 為1024×2048的探測(cè)器為例對(duì)最優(yōu)投影角度的布局進(jìn)行說明.當(dāng)|ri|=0時(shí),最少需要5 個(gè)投影角度可實(shí)現(xiàn)精確重建.若要求投影像素?cái)?shù)最少,則選取的投影矢量為[(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)所示,該條件下探測(cè)器限制在較小的角度范圍內(nèi).該投影矢量方向的二維Mojette變換投影結(jié)果如圖5(a)和圖5(b)所示,投影像素?cái)?shù)分別為64×1009,64×1009,64×946,64 ×946 和64×883.若探測(cè)器角度范圍不受限制,選取投影矢量為[(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)所示,投影像素?cái)?shù)分別為64×1009,64×1072,64×1198,64×1450和64×1513.

若 選取 的投 影矢 量中 | pi|>1 或 | qi|>1 ,此時(shí)對(duì)應(yīng)的 | ri| 變小,則所需最少投影角度數(shù)會(huì)增多.例如滿足精確重建條件的另一組投影矢量為[(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)所示.對(duì)應(yīng)的投影結(jié)果如圖5(f)所示.

圖5 不同投影矢量對(duì)應(yīng)的投影結(jié)果 (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 的結(jié)果綜合表明: 在探測(cè)器像素?cái)?shù)受限的條件下,最優(yōu)的投影角度布局方案為水平面投影,即基于二維Mojette 變換的層析重建,此時(shí)所需投影角的數(shù)目最少.并且不論投影角度是否受限,其所需的投影探測(cè)器的像素?cái)?shù)都比三維布局方案的少.這種方式與傳統(tǒng)層析技術(shù)中所選取的探測(cè)器的空間布局方案完全一致.當(dāng)投影矢量無法滿足水平面投影時(shí),要選擇為精確重建條件,并且選擇|pi|和|qi|的值越小,所需的投影角度數(shù)和探測(cè)器的像素?cái)?shù)越少.

5 結(jié) 論

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