基于MATLAB的CT滤波反投影算法的实验研究_iradon滤波matlab

实验四 滤波反投影算法的实验研究
1、利用”iradon”函数采用S-L与R-L滤波器进行滤波实现滤波反投影算法。

I=phantom(256);
theta=0:1:179;
P=radon(I,theta);
rec=iradon(P,theta,'linear','None');
rec_RL=iradon(P,theta,'Ram-Lak');
rec_SL=iradon(P,theta,'linear','Shepp-Logan');
figure;
subplot(2,2,1);imshow(I,[]),title('原始图像');
subplot(2,2,2);imshow(rec,[]),title('直接反投影图像');
subplot(2,2,3);imshow(rec_RL,[]),title('RL滤波反投影图像');
subplot(2,2,4);imshow(rec_SL,[]),title('SL滤波反投影图像');

1
2
3
4
5
6
7
8
9
10
11
12

2、利用FBP算法重建图像
1）对椭圆平行束投影数据进行滤波反投影重建，分别采用S-L与R-L滤波器进行滤波

clear;
N=256;
I = phantom('Shepp-Logan',N);
delta=pi/180;
theta=0:1:179;
theta_num=length(theta);
d=1;
%%产生投影数据
P=radon(I,theta);
[mm,nn]=size(P);
e=floor((mm-N-1)/2+1)+1;
P=P(e:N+e-1,:);
P1=reshape(P,N,theta_num);
%%
fh_RL=medfuncRlfilterfunction(N,d);
fh_SL=medfuncSlfilterfunction(N,d);
%%
rec=zeros(N);
for m=1:theta_num
    pm=P1(:,m);
    Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
            n=floor(Xrm);
            t=Xrm-floor(Xrm);
            n=max(1,n);
            n=min(n,N-1);
            p=(1-t)*pm(n)+t*pm(n+1);
            rec(N+1-k1,k2)=rec(N+1-k1,k2)+p;
        end
    end
end
%%
rec_RL=zeros(N)
for m1=1:theta_num
    pm1=P1(:,m1);
    pm_RL=conv(fh_RL,pm1,'same');
    Cm1=(N/2)*(1-cos((m1-1)*delta)-sin((m1-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm1=Cm1+(k2-1)*cos((m1-1)*delta)+(k1-1)*sin((m1-1)*delta);
            n1=floor(Xrm1);
            t1=Xrm1-floor(Xrm1);
            n1=max(1,n1);
            n1=min(n1,N-1);
            p1=(1-t1)*pm_RL(n1)+t1*pm_RL(n1+1);
            rec_RL(N+1-k1,k2)=rec_RL(N+1-k1,k2)+p1;
        end
    end
end
%%
rec_SL=zeros(N)
for m2=1:theta_num
    pm2=P1(:,m2);
    pm_SL=conv(fh_SL,pm2,'same');
    Cm2=(N/2)*(1-cos((m2-1)*delta)-sin((m2-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm2=Cm2+(k2-1)*cos((m2-1)*delta)+(k1-1)*sin((m2-1)*delta);
            n2=floor(Xrm2);
            t2=Xrm2-floor(Xrm2);
            n2=max(1,n2);
            n2=min(n2,N-1);
            p2=(1-t2)*pm_SL(n2)+t2*pm_SL(n2+1);
            rec_SL(N+1-k1,k2)=rec_SL(N+1-k1,k2)+p2;
        end
    end
end
%%
figure;
subplot(2,2,1);imshow(I,[]),title('原始图像');
subplot(2,2,2);imshow(rec,[]),title('直接反投影图像');
subplot(2,2,3);imshow(rec_RL,[]),title('RL滤波反投影图像');
subplot(2,2,4);imshow(rec_SL,[]),title('SL滤波反投影图像');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

2）对Shepp-Logan图平行束投影数据进行滤波反投影重建，别采用S-L与R-L滤波器进行滤波

clear;
N=256;
I=phantom(N);
delta=pi/180;
theta=0:1:179;
theta_num=length(theta);
d=1;
%%产生投影数据
P=radon(I,theta);
[mm,nn]=size(P);
e=floor((mm-N-1)/2+1)+1;
P=P(e:N+e-1,:);
P1=reshape(P,N,theta_num);
%%
fh_RL=medfuncRlfilterfunction(N,d);
fh_SL=medfuncSlfilterfunction(N,d);
%%
rec=zeros(N);
for m=1:theta_num
    pm=P1(:,m);
    Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
            n=floor(Xrm);
            t=Xrm-floor(Xrm);
            n=max(1,n);
            n=min(n,N-1);
            p=(1-t)*pm(n)+t*pm(n+1);
            rec(N+1-k1,k2)=rec(N+1-k1,k2)+p;
        end
    end
end
%%
rec_RL=zeros(N)
for m1=1:theta_num
    pm1=P1(:,m1);
    pm_RL=conv(fh_RL,pm1,'same');
    Cm1=(N/2)*(1-cos((m1-1)*delta)-sin((m1-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm1=Cm1+(k2-1)*cos((m1-1)*delta)+(k1-1)*sin((m1-1)*delta);
            n1=floor(Xrm1);
            t1=Xrm1-floor(Xrm1);
            n1=max(1,n1);
            n1=min(n1,N-1);
            p1=(1-t1)*pm_RL(n1)+t1*pm_RL(n1+1);
            rec_RL(N+1-k1,k2)=rec_RL(N+1-k1,k2)+p1;
        end
    end
end
%%
rec_SL=zeros(N)
for m2=1:theta_num
    pm2=P1(:,m2);
    pm_SL=conv(fh_SL,pm2,'same');
    Cm2=(N/2)*(1-cos((m2-1)*delta)-sin((m2-1)*delta));
    for k1=1:N
        for k2=1:N
            Xrm2=Cm2+(k2-1)*cos((m2-1)*delta)+(k1-1)*sin((m2-1)*delta);
            n2=floor(Xrm2);
            t2=Xrm2-floor(Xrm2);
            n2=max(1,n2);
            n2=min(n2,N-1);
            p2=(1-t2)*pm_SL(n2)+t2*pm_SL(n2+1);
            rec_SL(N+1-k1,k2)=rec_SL(N+1-k1,k2)+p2;
        end
    end
end
%%
figure;
subplot(2,2,1);imshow(I,[]),title('原始图像');
subplot(2,2,2);imshow(rec,[]),title('直接反投影图像');
subplot(2,2,3);imshow(rec_RL,[]),title('RL滤波反投影图像');
subplot(2,2,4);imshow(rec_SL,[]),title('SL滤波反投影图像');
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

function P = medfuncParallelBeamForwardProjection(theta,N,N_d)
    shep=[   1   .69     .92     0       0       0
        -.8 .6624   08740   0       -.0184  0
        -.2 .1100   .3100   .22     0       -10
        -.2 .1600   .4100   -.22    0       15
        .1  .2100   .2500   0       .35     5
        .1  .0460   .0560   .06       .1      15
        .1  .0260   .0460   .13       -.1     20
        .1  .0460   .0130   -.08    -.605   -5
        .1  .0430   .0230   .03       -.606   3
        .1  .0230   .0460   .06     -.605   0];
    theta_num = length(theta);
    P = zeros(N_d,theta_num); % 存放投影数据
    rho = shep(:,1).';  % 圆对应密度
    be =0.5*N*shep(:,2).'; %椭圆短半轴
    ae = 0.5*N*shep(:,3).'; % 椭圆长半轴
    xe =0.5*N*shep(:,4).'; %圆中心x坐标
    ye = 0.5*N*shep(:,5).'; %圆中心y坐标
    alpha = shep(:,6).'; % 圆旋转角度
    alpha = alpha * pi/180;
    theta = theta * pi/180; % 角度换算成弧度
    TT =-(N_d-1)/2:(N_d-1)/2;% 探测器坐标
    for k1 = 1 : theta_num
        P_theta = zeros(1,N_d);% 存放每一角度投影数据
        for k2 = 1: max(size(xe))
            a =(ae(k2) * cos(theta(k1)-alpha(k2)))^2+(be(k2) * sin(theta(k1)-alpha(k2)))^2;
            temp = a- (TT-xe(k2) * cos(theta(k1))-ye(k2) * sin(theta(k1))).^2;
            ind = temp>0;% 根号内值需为非负
            P_theta(ind)= P_theta(ind) +rho(k2) * (2 * ae(k2) * be(k2) * sqrt(temp(ind)))./a;
        end
        P(:,k1) = P_theta.';
    end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

function fh_RL=medfuncRlfilterfunction(N,d)
fh_RL=zeros(1,N);
for k1=1:N
    fh_RL(k1)=-1/(pi*pi*((k1-N/2-1)*d)^2);
    if mod(k1-N/2-1,2)==0
        fh_RL(k1)=0;
    end
end
fh_RL(N/2+1)=1/(4*d^2);
1
2
3
4
5
6
7
8
9

function fh_SL=medfuncSlfilterfunction(N,d)
fh_PL=zeros(1,N);
for k1=1:N
    fh_SL(k1)=-2/(pi^2*d^2*(4*(k1-N/2-1)^2-1));
end
1
2
3
4
5