卷積反投影重建二維

1、卷積反投影圖像重建1 1 反投影重建基本介紹A設(shè)待重建圖像為 a(x,y)a(x,y), ,它的二維傅氏變換為A(SI,CO2)=A(PQ)。根據(jù)中心切片定A理，A(P*)可通過 a(x,y)a(x,y)在不到視角由下的投影p/x)的一維傅氏變換求得。即：A(,2)=A(：,)=p(x)=P(：)=P(：,)待建圖像：1a(r)=a(x,y)=F2A(1,2)二11A(1,.2)ei(X.2y)d.，1d.24 二2-：_：、一1.A二A(P*)ei2*rc0s(e)|P|dPd1,0-00=PJP(p*)ei2JTPrcos(8_*)|PdPd 巾，0-OO=/d廣|pp(P,*)ei2jT

2、pcos(0)dP03H因?yàn)閤r=rcos(8-*),所以有:xf2y=2”.(xcos二，ysin)=2二Pxr=2：rcos(-)同時(shí)：do1d2=JdPd4先來看該式的第二個(gè)積分：口歸匕性=口耳-3dp=h(*)P(xr,)|xr=rcosf一)=g(x,)cosh-)-grcos(-),1式中:g(x,戶h(x)p(x,(1.4)(1.4)式(3.10)(3.10)的物理意義是投影p(xr仲)經(jīng)過傳遞函數(shù)為|H=|H=F F1h(x)的濾波器后得到的修正后的投影g(x2)在滿足xr=rcos(8-旬時(shí)的值。將(3.11)(3.11)代入(3.8)(3.8),得到：Aa(r,B)=(gr

3、cos(日一4),4d(1.5)(1.5)稱為濾波反投影方程，其物理意義是經(jīng)過給定點(diǎn)(r,的所有濾波后的投影在=0元范圍內(nèi)的累加一反投影重建，得出(r,B)點(diǎn)的像素值。(1.1)2二cos2二sin-2nPsinO2nPcos二4二2P(1.2)(1.3)可見，濾波(卷積)反投影算法的具體包含三大步：把在固定視角下測得的投影p(x,鶴經(jīng)過濾波，得到濾波后的投影g(xr,*)；(2)(2)對(duì)每一個(gè)把g(xr*)反投影于滿足xr=rcos(8-4)的射線上的所有各點(diǎn)(方);(3)(3)將步驟(2)(2)中的反投影值對(duì)所有 006gn6gn 進(jìn)行累加(積分)，得到重建后的圖像。2 重建流程2.1首先

4、我們利用 phantomphantom()()函數(shù)產(chǎn)生一個(gè)頭部幻影圖像，用以檢測二維重建算法，代碼如下：I=phantom(256);subplot(2,2,1)imshow(I,);Utle(256*256 原始圖像);效果圖如圖 1 1 所示,為一個(gè)大橢圓和幾個(gè)小橢圓。2.2初始參數(shù)設(shè)置重建采用的是平移加旋轉(zhuǎn)的掃描方式，射線源在某一角度下水平移動(dòng)，將物體全部照射后旋轉(zhuǎn)一角度，如此重復(fù)，在這個(gè)過程中探測器相應(yīng)地運(yùn)動(dòng)以接收 X X 射線。根據(jù)此原理，將重建程序的初始參數(shù)設(shè)置如下：N,N=size(I);z=2*ceil(norm(size(I)-floor(size(I)-1)/2)-1)+3

5、;%radon 變換默認(rèn)平移點(diǎn)數(shù)/角度Nt=360;%角度采樣點(diǎn)數(shù)Nd=N;%平移數(shù)x=pi/180;%角度增量d=N/Nd;%平移步長theta=1:Nt;a=zeros(N);2.3產(chǎn)生無噪聲投影數(shù)據(jù)256256原始圖圖1R,xp=radon(I,theta);e=floor(z-Nd)+2;R=R(e:(Nd+e-1),:);R1=reshape(R,256,360);radon(I,theta)產(chǎn)生 I I 投影，默認(rèn) z z 點(diǎn)/角度，即使指定 N N 點(diǎn)也是 z z 點(diǎn).所以為避免重建圖像放大或縮小，下面計(jì)算取投影時(shí)需補(bǔ)償，補(bǔ)償量 e e 如對(duì) 256256 的圖像，補(bǔ)償為 55,

6、55,即 pmllpmll 勺第 5555 個(gè)點(diǎn)作為計(jì)算用的第一個(gè)投影. .2.4添加噪聲并將有噪聲平行投影進(jìn)行顯小mm,nn=size(R1);di=lognrnd(0,0.15,mm,nn);R1=10*(R1-min(R1(:)/(max(R1(:)-min(R1(:);I0=1.5e5;%incidentphotons;decreasethisforsimulatinglowdosescansrand(state,0),randn(state,0);yi=poissrnd(I0*di.*exp(-R1)+3*randn(size(R1);ifany(yi(:)=0)warn(%dof%

7、dvaluesare0insinogram!,.sum(yi(:)=0),length(yi(:);endR1=log(I0./max(yi,0.01);%noisysinogramR1=max(R1,0);%顯示ff=2;uu=22000;v=ff*exp(R1/uu);subplot(2,2,2)imagesc(RI);title(256*360 有噪聲平行投影,);colormap(gray)colorbarQ=reshape(R1,256,360);效果圖如圖-2:2.5濾波器的選擇與設(shè)計(jì)Bf2h(s)=J|5exp(/2;rSs)dS=最基本的從投影重建圖像的濾波器:示)。19711

8、971 年提出的 R-LR-L 重建濾波器(下圖中實(shí)線表空域表達(dá)式為: :圖-3R-L-3R-L 濾波器示意圖(2.1(2.1) )其中 B 為截至頻率。若2,d 為空間采樣間隔，可解256#360 有噪聲平行投影50100150200250圖-2isnr100200300151050出離散的濾波器空間脈沖響應(yīng):fm=04d2hd)=fg(.v,w離散化就S得r*2.A/Afg式泄)二pm：d,心)=Zg/QH用-k)-?nLMJr=其中從濾波因尸方采樣得到Aft(HI)=考慮g(JW)線性內(nèi)插為歐必)gG心)=治(川)+(；-叫心(川+1)-京(加)a其中mdsin+l)rf反投影枳分NT/

9、(.Y.1)=Zg(MSQ5十sinn-0在角坐標(biāo)離散位NT-Ag|_(/Arcos?A+/AvsinnA),nAPI=O/,j0,h2,1t-代碼設(shè)計(jì)：form=1:Ntpm=Q(:,m);u=conv(hl,pm);pm=u(k);Cm=(N-1)/2)*(1-cos(m-1)*x)-sin(m-1)*x);fori=1:Nforj=1:NXrm=Cm+(j-1)*cos(m-1)*x)+(i-1)*sin(m-1)*x);ifXrm(Nd-1)n=Nd-1;endp=(1-t)*pm(n)+t*pm(n+1);a(N+1-i,j)=a(N+1-i,j)+p;endendend重建后的圖像

10、如圖-3-3 所示：程序中還包含對(duì)重建結(jié)果的評(píng)價(jià)/歸一化均方距離判據(jù)/歸一化平均絕對(duì)距離判據(jù)以及程序的運(yùn)行時(shí)間等.不再進(jìn)行詳細(xì)介紹。程序的最終運(yùn)行效果如下：重建圖像圖-31.510.50050100第230ff的像素圖-4 程序最終運(yùn)行效果圖3 分析反投影重建方法包括卷積反投影重建的缺點(diǎn)是會(huì)產(chǎn)生星狀偽跡，原因分析如下：斷層平面中某一點(diǎn)的密度值可以看作是這一平面內(nèi)所有經(jīng)過該點(diǎn)的射線的投影值之和（的均值）。整幅重建圖像可以看作是所有方向下的投影累加而成。射線標(biāo)號(hào)示于圖 5 5 中，像素值（代表密度）分別 x x1,x,x2,X,X3,X,X4,賦值如下：x1=5,x2=0,x3=2,x4=18根據(jù)

11、投影的定義（某條射線投影值為該條射線穿過的所有的像素值之和） ，每條射線的投影p（i=1,2- -）為：Pl=xi*x2=5, ,p2=x3*x4=20, ,p3=xi+x3=7p4=x2x4=18,p5=x3=2, ,p6=x1x4=23P7=x2=0256*256原始圖像256*360有噪聲平行投影inr上6iflV:w|重建圖像圖-5 斷層像素值和射線根據(jù)反投影重建算法的物理意義， 重建圖像中各像素， 得到:xi=pi+p3+p6=35,x2=pi+p4+p7=23,x3=p2p3p5=29x4=p2p4P6=61,共 9 頁第 8 頁重建后的圖像如圖 6(b)6(b)所示，可以看出原圖

12、像中像素值不為零的點(diǎn)反投影重建后仍較突出，但原圖中像素值為零的點(diǎn)，經(jīng)反投影重建后不再為零，即有偽跡。有時(shí)為了使重建后圖像的像素值更接近于原圖的像素值，在求反投影時(shí)，把數(shù)據(jù)除以投影的數(shù)目(即射線數(shù))，如圖 6(c)6(c)所示。因此有：1XkPk,inpi=1該式可作為反投影重建算法的計(jì)算式。其中 X Xk表示像素 k k 的值，Pk,i表示經(jīng)過像素k的第i條射線投影，np表示圖像內(nèi)的射線條數(shù)。圖 7(a)7(a)表示空間中一個(gè)孤立點(diǎn)源 A,A,密度為 1 1。經(jīng)過 A A 點(diǎn)的三條射線也示于圖中。射線束理論上可以很多， 取三條示意。 不經(jīng)過 A A 點(diǎn)的射線投影為零， 經(jīng)過 A A 點(diǎn)的射線投影值均為 1,1,P1=P2=P3=1。5021835262961(3.1)(3.1)(a)原圖像像素值(b)反投影重建后圖像圖 6 反投影示例(c)求平均后圖像圖 7 孤立點(diǎn)源的反投影重建及星狀偽跡經(jīng)反投影重建后，得到 A A 點(diǎn)的像素值為fA=(R+P2+P3)/3=1。