卷積反投影重建(二維)

1、計算機(jī)圖像處理技術(shù)結(jié)課(報告)卷積反投影圖像重建1反投影重建基本介紹A設(shè)待重建圖像為a(x, y)，它的二維傅氏變換為 A( .1, .2)=A(幾)。根據(jù)中心切片定A理，A( ：?,)可通過a(x,y)在不同視角 '下的投影p (Xr)的一維傅氏變換求得。即：A( 1，2)= A(，)二 F1 | P (xr )二 P ()= P()待建圖像：_1a(rj)二 a(x,y)二 F?一A(匕)厶 i i A(2)ei( 1X '2y)d24 二=j" A(P)ei2用rcos(日P|d Pd©(1.1)0=r P(pf)ei2“PrcosZ)| p|d Pd

2、©0 -=fd 叮 °°|p P(Pf)ei2jlFrcosdd pJ0因為xrcosh -)，所以有：x r：'2y =2 J(xcos：；亠 ysin ) =2二x, =2：rcosf )同時：d創(chuàng)dco2 = J d Pd©2 cos2二 sin '-2 兀 Ps in©2 兀 Pcos©(1.2)先來看該式的第二個積分：門護(hù)(匕怩2叫叫戸=彷|時力嚴(yán)lecosg dP二皿片)以片，)lx嚴(yán) cos-J二g(x,)" cosd-)=grcos( -),式中:g(x , ) h (x ) p (x ,(1

3、.4)式(3.10)的物理意義是投影p(xr,)經(jīng)過傳遞函數(shù)為Fdh(Xr)的濾波器后得到的修正后的投影g(x,)在滿足xr二rcos()時的值。將(3.11)代入(3.8)，得到:Aa(r,日)=J0 gr cos(日一©),©d©(1.5)稱為濾波反投影方程，其物理意義是經(jīng)過給定點(r,R的所有濾波后的投影在 -0-：范圍內(nèi)的累加一反投影重建，得出(rj)點的像素值。計算機(jī)圖像處理技術(shù)結(jié)課(報告)可見，濾波(卷積)反投影算法的具體包含三大步：(1)把在固定視角】下測得的投影P(Xr,)經(jīng)過濾波，得到濾波后的投影g(Xr,)； 對每一個1，把g(xr, i)反投

4、影于滿足xr rcos(J的射線上的所有各點(r,v)；(3)將步驟 中的反投影值對所有Ov乜二進(jìn)行累加(積分)，得到重建后的圖像。2重建流程256-2 圖像圖12.1首先我們利用phantom ()函數(shù)產(chǎn)生一個頭部幻影圖像，用以檢測二維重建算 法，代碼如下：I=pha ntom(256);subplot(2,2,1)imshow(l,);title('256*256 原始圖像');效果圖如圖1所示，為一個大橢圓和幾個小橢圓。2.2初始參數(shù)設(shè)置重建采用的是平移加旋轉(zhuǎn)的掃描方式，射線源在某一角度下水平移動，將物體全部 照射后旋轉(zhuǎn)一角度，如此重復(fù)，在這個過程中探測器相應(yīng)地運(yùn)動以接收

5、X射線。根據(jù)此原理,將重建程序的初始參數(shù)設(shè)置如下：N,N=size(l);z=2*ceil(norm(size(I)-floor(size(I)-1)/2)-1)+3;% radon 變換默認(rèn)平移點數(shù)/角度Nt=360;%角度采樣點數(shù)Nd=N; %平移數(shù)x二pi/180; %角度增量 d=N/Nd; %平移步長theta = 1:Nt;a=zeros(N);2.3產(chǎn)生無噪聲投影數(shù)據(jù)共9頁第4頁計算機(jī)圖像處理技術(shù)結(jié)課（報告）R,xp = rad on (l,theta);e二floor(z-Nd)/2)+2;R=R(e:(Nd+e-1),:);R仁 reshape(R,256,360);rado

6、n(I,theta)產(chǎn)生I投影，默認(rèn)z點/角度，即使指定N點也是z點.所以為避免重建 圖像放大或縮小,下面計算取投影時需補(bǔ)償,補(bǔ)償量e如對256的圖像,補(bǔ)償為55,即pm勺 第55個點作為計算用的第一個投影.2.4添加噪聲并將有噪聲平行投影進(jìn)行顯示mm, nn =size(R1);di=log nrn d(0,0.15,mm, nn);R1= 10*(R1-mi n(R1(:)/( max(R1(:)-mi n(R1(:);I0 = 1.5e5; % in cide nt phot ons; decrease this for simulati ng "low dose"

7、sca ns ran d('state', 0), randn ('state', 0);yi= poissrnd(I0 * di.*exp(-R1)+3*ra ndn (size(R1);if any (yi(:) = 0)warn ('%d of %d values are 0 in sinogram!', .sum(yi(:)=0), le ngth(yi(:);endR1 = log(I0 ./ max(yi,0.01); % n oisy si nogramR1=max(R1,0);%顯示ff=2;uu=22000;v=ff*exp(R

8、1/uu);共9頁第5頁計算機(jī)圖像處理技術(shù)結(jié)課（報告）15subplot(2,2,2)imagesc(RI);title('256*360有噪聲平行投影');colormap(gray)colorbarQ=reshape(R1,256,360);效果圖如圖-2:256*360有噪聲平行投影50100150200250100 200 300圖-2-102.5濾波器的選擇與設(shè)計最基本的從投影重建圖像的濾波器：1971年提出的R-L重建濾波器（下圖中實線表 示）。圖-3 R-L濾波器示意圖 空域表達(dá)式為：Bf2h(s) = J |ijexp(y2S5)</5 =*嚴(yán)、sm()a

9、*訂s、 i snr()jrs=JIS-()2* J_ d1(2.1)其中B為截至頻率。若, d為空間采樣間隔，可解共9頁第7頁計算機(jī)圖像處理技術(shù)結(jié)課（報告）出離散的濾波器空間脈沖響應(yīng):7m = 0=«0iv is eivn(2.2)丁丿=m is oddMATLAB弋碼如下：g=-(Nd/2-1):(Nd/2);for i=1:256if g(i)=0hl(i)=1/(4*dA2);else if mod(g(i),2)=0hl(i)=0;elsehl(i)=(-1)/(piA2*dA2*(g(i)A2);endendendk=Nd/2:(3*Nd/2-1);2.6重建圖像卷積 g

10、(斗 0) = gGE)"G)反投影=0離散化3)得-11<tn< gff(»?) = g(d,?a)= Y gK)h(m-k)買屮從濾波禺MG）采樣得到h(in) = dh(itnl)考慮g (加)線性內(nèi)插為£(/")g(兵 /?A) =乳(m) + (J - 加)衣(川 +1) - 或(山) a其中 wd <s< (m + l)rf反投影積分JV-1/ (.Y. 1 ) = A(ACQSf7 十 win M )，/? n-0il角坐標(biāo)r(X. V)離散值i9J 0»±L ±2< -代碼設(shè)計：

11、for m=1:Ntpm=Q(:,m);u=c on v(hl,pm);pm=u(k);Cm=(N-1)/2)*(1-cos(m-1)*x)-si n(m-1)*x);for i=1:Nfor j=1:NXrm二Cm+(j-1)*cos(m-1)*x)+(i-1)*si n(m-1)*x);if Xrm<1n=1;t=abs(Xrm)-floor(abs(Xrm);elsen=floor(Xrm); t=Xrm-floor(Xrm); end if n>(Nd-1)n=Nd-1;endp=(1-t)*pm( n)+t*pm( n+1); a(N+1-i,j)=a(N+1-i,j)+

12、p;end重建圖像圖-3endend重建后的圖像如圖-3所示:程序中還包含對重建結(jié)果的評價/歸一化均方距離判據(jù)/歸一化平均絕對距離判據(jù) 以及程序的運(yùn)行時間等不再進(jìn)行詳細(xì)介紹。程序的最終運(yùn)行效果如下：共9頁第10頁計算機(jī)圖像處理技術(shù)結(jié)課（報告）2於360有噪聲平行投影256*256原始圖像50100150200250100 200 300151050重建圖像1.50.550100第2301亍的像素圖-5斷層像素值和射線圖-4程序最終運(yùn)行效果圖3分析反投影重建方法包括卷積反投影重建的缺點是會產(chǎn)生星狀偽跡，原因分析如下： 斷層平面中某一點的密度值可以看作是這一平面內(nèi)所有經(jīng)過該點的射線的投影值之和 （

13、的均值）。整幅重建圖像可以看作是所有方向下的投影累加而成。射線標(biāo)號示于圖5中，像素值（代表密度）分別 , x2, x3, x4，賦值如下：x1 =5, x2 =0 , X3 =2 , x4 =18根據(jù)投影的定義（某條射線投影值為該條射線穿過的所 有的像素值之和），每條射線的投影Pi （ i 1,2）為：Pl = Xi x2 = 5 ,P = x3 x4 = 20 , p3 = xi x3 = 7p = x2 x4 = 18,P5 = X3 = 2 ,p6 =為 X4 = 23P7 = X2 二 0根據(jù)反投影重建算法的物理意義，重建圖像中各像素，得到:X P1 P3 P6 =35，X2 =P1

14、P4 P7 =23，X3 =P2 P3 P5 =29 X4 =P2 P4 P6 =61，共9頁第8頁計算機(jī)圖像處理技術(shù)結(jié)課(報告)50218(a)原圖像像素值35262961(b)反投影重建后圖像圖6反投影示例8.7(c)求平均后圖像共9頁第13頁重建后的圖像如圖6(b)所示，可以看出原圖像中像素值不為零的點反投影重建后 仍較突出，但原圖中像素值為零的點，經(jīng)反投影重建后不再為零，即有偽跡。有時為了 使重建后圖像的像素值更接近于原圖的像素值，在求反投影時，把數(shù)據(jù)除以投影的數(shù)目 (即射線數(shù))，如圖6(c)所示。因此有：nPXknpi 二1Pk,i(3.1)共9頁第#頁共9頁第#頁該式可作為反投影重建算法的計算式。其中Xk表示像素k的值，Pk,i表示經(jīng)過像素 k的第i條射線投影，np表示圖像內(nèi)的射線條數(shù)。圖7(a)表示空間中一個孤立點源 A，密度為1。經(jīng)過A點的三條射線也示于圖中。 射線束理論上可以很多，取三條示意。不經(jīng)過 A點的射線投影為零，經(jīng)過 A點的射線 投影值均為1， p1 p2 p3 =1。(a)孤立點源A及三條射線(b