%?卷積反投影重建;莊天戈《CT原理與算法》/卷積反投影重建/f3.27;筆束平移旋轉掃描
%%?
clear;
clc;
close?all;
tic;
%%?initial
I=phantom（256）;
subplot（221）
imshow（I[]）;
title（‘256*256原始圖像‘）;
[NN]=size（I）;
z=2*ceil（norm（size（I）-floor（（size（I）-1）/2）-1））+3;%?radon變換默認平移點數/角度
Nt=360;%?角度采樣點數
Nd=N;%?平移數
x=pi/180;%?角度增量
d=N/Nd;%?平移步長
theta?=?1:Nt;
a=zeros（N）;
%%
%?產生無噪聲投影數據
[Rxp]?=?radon（Itheta）;%?產生I投影默認z點/角度即使指定N點也是z點.
????????????????????????%?所以為避免重建圖像放大或縮小下面計算取投影時需補償補償量e
????????????????????????%?如對256的圖像補償為55即pm的第55個點作為計算用的第一個投影
e=floor（（z-Nd）/2）+2;
R=R（e:（Nd+e-1）:）;
R1=reshape（R256360）;
%?添加噪聲
[mmnn]=size（R1）;
di=lognrnd（00.15mmnn）;
R1=?10*（R1-min（R1（:）））/（?max（R1（:））-min（R1（:）））;
I0?=?1.5e5;?%?incident?photons;?decrease?this?for?simulating?“low?dose“?scans
rand（‘state‘?0）?randn（‘state‘?0）;
yi=?poissrnd（I0?*?di.*exp（-R1））+3*randn（size（R1））;

if?any（yi（:）?==?0）
??warn（‘%d?of?%d?values?are?0?in?sinogram!‘?...
???????sum（yi（:）==0）?length（yi（:）））;
end
R1?=?log（I0?./?max（yi0.01））;?%?noisy?sinogram
R1=max（R10）;?%?R1?加噪的投影數據

%?load?R3.MAT
ff=2;
uu=22000;
v=ff*exp（R1/uu）;
subplot（222）
imagesc（R1）;
title（‘256*360有噪聲平行投影‘）;
colormap（gray）
colorbar
Q=reshape（R1256360）;
%%??

%?designing?RL?filter?
g=-（Nd/2-1）:（Nd/2）;
for?i=1:256
????if?g（i）==0
????????hl（i）=1/（4*d^2）;
????else?if??mod（g（i）2）==0
????????????hl（i）=0;
????????else
????????????hl（i）=（-1）/（pi^2*d^2*（g（i）^2））;
????????end
????end
end
k=Nd/2:（3*Nd/2-1）;%?取卷積結果時用
%%?
%?重建過程
for?m=1:Nt
????%?reading?projection
????pm=Q（:m）;%?讀取投影數據
????u=conv（hlpm）;%?卷積