%%投影
I?=?phantom（‘Modified?Shepp-Logan‘256）;
nmax=ceil（length（I）*（sqrt（2）））;
I1=zeros（512）;
I1（128:383128:383）=I;
proj=zeros（512180）;
for?theta=1:180
????I2=imrotate（I1-theta‘bilinear‘）;%原圖旋轉theta度，得到I2
????temp=length（I2）;
????temp=round（temp/2-255）;
????I3=I2（temp:temp+511temp:temp+511）;
????b=sum（I3）;%b是投影向量
????proj（:theta）=b‘;%賦值投影矩陣
end
imshow（proj[]）
%%投影

%?width?=?2^nextpow2（size（I1））;
width=512;
H?=?2*[0:（width/2-1）?width/2:-1:1]‘/width;%R-L濾波
proj_filtered?=?zeros（length（proj）180）;?
%?plot（H）
proj_fft?=?fft（proj?512）;
for?i?=?1:180????
????proj_filtered（:i）?=?proj_fft（:i）.*H;??
end??

%?M=256;
%?fbp?=?zeros（M）;?%?假設初值為0
proj_ifft?=?real（ifft（proj_filtered））;?
%?imshow（proj_ifft[]）
fbp?=?zeros（256）;?

%?5.?back-projection?to?the?x-?and?y-?axis??

for?i?=?1:180??
????%?rad?is?the?angle?of?the?projection?line??not?projection?angle??
????rad?=?i*pi/180;??
????for?x?=?（-256/2+1）:256/2??
????????for?y?=?（-256/2+1）:256/2??
????????????t?=?round（x*cos（rad+pi/2）+y*sin（rad+pi/2））;??
????????????fbp（x+256/2y+256/2）=fbp（x+256/2y+256/2）+proj_ifft（t+256i）/180;??
????????????
????????end??
????end??
%?????imshow（fbp[]）
end??
fbp?=?fbp/360;??

subplot（1?2?1）?imshow（I）?title（‘Original‘）
subplot（1?2?2）?imshow（fbp[]）?title（‘FBP‘）

%?I=dicomread（‘C:\Users\56326\Desktop\dcm\100.dcm‘）;
%?imshow（I[]）
%?info=?dicominfo（‘C:\Users\56326\Desktop\dcm\100.dcm‘）;
%?I1=imread（‘C:\Users\56326\Desktop\sonoscape\12345\20180425_145827\20180425_150454_1.jpg‘）;
%?I2=imread（‘C:\Users\56326\Desktop\sonoscape\12345\20180425_145827\20180425_150545_2.jpg‘）;
%?I3=imread（‘C:\Users\56326\Desktop\sonoscape\12345\20180425_145827\20180425_150603_3.jpg‘）;
%?
%?figure
%%
theta?=?1:180;??
??P=phantom（128）;
%?1.?projection??using?radon?function??
[Rxp]?=?radon（Ptheta）;??
width?=?2^nextpow2（size（R1））;??%?set?width?for?fft?transformation??
??
%?2.?do?fft?to?the?projection??
proj_fft?=?fft（R?width）;??
??
%?3.?filter??
%?Ramp?filter?function??from?0?to?width?then?to??0??
filter?=?2*[0:（width/2-1）?width/2:-1:1]‘/width;??
proj_filtered?=?zeros（width180）;??
for?i?=?1:180??
????proj_filtered（:i）?=?proj_fft（:i）.*filter;??
end??
??
%?4.?do?ifft?to?the?filtered?projection??
proj_ifft?=?real（ifft（proj_filtered））;?%?get?the?real?part?of?the?result??
??
%?5.?back-projection?to?the?x-?and?y-?axis??
fbp?=?zeros（128）;?%?asign?the?original?value?0??
for?i?=?1:180??
????%?rad?is?the?angle?of?the?projection?line??not?projection?angle??
????rad?=?i*pi/180;??
????for?x?=?（-128/2+1）:128/2??
????????for?y?=?（-128/2+1）:128/2??
????????????t?=?round（x*cos（rad+pi/2）+y*sin（rad+pi/2））;??
????????????fbp（x+128/2y+128/2）=fbp（x+128/2y+128/2）+proj_ifft（t+round（size（R1）/2）i）;??
????????end??
????end??
????imshow（fbp[]）
end??
fbp?=?fbp/180;??
??
%?6.?show?the?result??
subplot（1?2?1）?imshow（P）?title（‘Original‘）??
subplot（1?2?2）?imshow（fbp）?title（‘FBP‘）??
%%

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????3022??2018-05-16?14:41??fbp-.m

-----------?---------??----------?-----??----

?????????????????3022????????????????????1