function??[a]=dps（P）
tic;
P=phantom（256）;
%P=imread（‘C:\0劉金豆\醫(yī)學圖像處理實驗重建\lenna.bmp‘）;
%P=rgb2gray（P）;
%P=?imresize（P[256?256]）;

[NN]=size（P）;%獲取圖像P的size，存入N
subplot（231）;imshow（P）;title（[int2str（N）‘*‘int2str（N）‘原始圖像‘]）;

%先進行自定義radon變換------------------------------------------------------------
thm=45;????????????????%45度時會出現(xiàn)最大尺寸
pre?=?imrotate（Pthm）;???%將圖像P（圖像的數(shù)據(jù)矩陣）繞圖像的中心點旋轉(zhuǎn)thm度，
?????????????????????????%正數(shù)表示逆時針旋轉(zhuǎn)，?負數(shù)表示順時針旋轉(zhuǎn)。返回旋轉(zhuǎn)后的圖像矩陣pre。
[mmaxnmax]?=?size（pre）;
s=1;????%定義步長為1

%創(chuàng)建一個180*nmax的空白圖片，用以存儲投影后的線狀圖片
Final?=?zeros（180/snmax）;%這里180代表180角度，每個角度投影成為一條線
t?=?1;
for?theta?=?1:s:179
%對原圖旋轉(zhuǎn)一個角度，求和（線積分）
Protate?=?imrotate（Ptheta）;?%對圖像P旋轉(zhuǎn)后的圖Protate???
Pf?=?sum（Protate1）;%對Protate矩陣列求和
[mrealnreal]=size（Pf）;???????????%計算實際尺寸
%確定起始點
if?（nmax?-?nreal）/2-floor（（nmax?-?nreal）/2）?==?0?????????%y?=?floor（x）?函數(shù)將x中元素向下取整
????From?=?floor（（nmax?-?nreal）/2?+?1）;%總點數(shù)為偶數(shù)時
???else
????From?=?floor（（nmax?-?nreal）/2）?+?1;%總點數(shù)為奇數(shù)時
????end
%確定結(jié)束點
End?=?floor（（nmax-nreal）/2）?+?nreal;
%將一個角度Radon變換后的線狀圖存入結(jié)果圖像的某一行
Final（180/s-tFrom:End）?=?Pf;?%從最底下一行開始存起
%上移一行
t?=?t?+?1;
end
%再逆時針旋轉(zhuǎn)
R=imrotate（Final90）;
subplot（232）;imshow（R[]）;title（‘自定義投影后圖像‘）;

z=2*ceil（norm（size（P）-floor（（size（P）-1）/2）-1））+3;%?radon變換默認平移點數(shù)/角度
e=floor（（z-N）/2）+2;???
R1=R（e:（N+e-1）:）;
[mmnn]=size（R1）;
subplot（233）;
imagesc（R1）;
title（[int2str（mm）‘*‘int2str（nn）‘正弦圖‘]）;
colormap（gray）;
colorbar;

%濾波函數(shù)構(gòu)造------------------------------------------------------------
g=1-N:N-1;?
%??d=1;
%??R-L濾波函數(shù)
%??for?i=1:2*N-1
%?????????if?g（i）==0
%?????????h（i）=1/（4*（d^2））;
%?????else?if?mod（g（i）2）==0
%?????????????h（i）=0;
%?????????else?
%?????????????h（i）=（-1）/（pi^2*d^2*（g（i）^2））;
%?????????end
%?????end
%??end

%S-L濾波函數(shù)
%?d=1;
%??for?i=1:2*N-1
%????????h（i）=2/（pi^2*d^2*（1-4*g（i）^2））;
%??end

%Parzen濾波函數(shù)
?for?i=1:2*N-1
?????h（i）=（24*pi*g（i）*cos（pi*g（i））-96*sin（pi*g（i））-48*pi*g（i）*cos（pi*g（i）/2）+384*sin（pi*g（i）/2）-2*（pi^3）*（g（i）^3）-72*pi*g（i））/（4*（pi^5）*（g（i）^5））;
?end
?h（N）=0.0435;
?
%將投影與濾波函數(shù)卷積----------------------------------------------------
G=zeros（N180）;
for?m=1:180??????
????for?i=1:N
????????b=0;
????????for?n=1:N
????????????b=b+h（N+n-i）*R1（nm）;
????????????G（im）=b;
????????end
????end
end

%投影濾波后線性內(nèi)插進行圖像重建----------------------------------------------
a=zeros（N）;???????%重建圖像初始化，每個像素點像素值為0
ns=（N+1）/2;
for?m=1:180?????????%遍歷每個投影角度
????r=pi/180;???????%將角度換為弧度
????for?i=1:N??????
????????for?j=1:N????%遍歷原圖的每一個像素點
?????????????Xrm=（i-ns）*cos（m*r）+（j-ns）*sin（m*r）+ns-1;??%坐標轉(zhuǎn)換
????????????if?Xrm<1
????????????????n=1;???????????????????????????%內(nèi)插運算整數(shù)值
????????????????t=abs（Xrm）-floor（abs（Xrm））;?????%內(nèi)插運算小數(shù)值
????????????else
????????????????n=floor（Xrm）;
????????????????t=Xrm-floor（Xrm）;
????????????end
????????????if?n>N-1
????????????????n=N-1;
????????????end
????????????c=（1-t）*G（nm）+t*G（n+1m）;???????%內(nèi)插后的濾波投影值
????????????a（N+1-ji）=a（N+1-ji）+c;?????????%像

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????4305??2018-12-05?11:18??醫(yī)學圖像重建作業(yè)\FBP_GUI.m
?????文件????????3786??2019-01-05?10:17??醫(yī)學圖像重建作業(yè)\dps.m
?????目錄???????????0??2019-01-17?14:14??醫(yī)學圖像重建作業(yè)\