ct_fang.m

医学成像技术滤波反投影 内含R-L

function reconstruction
clear;
close all;
N=256;%N＊N大小图像
delta=3; %相当于Φ
a = 100;
b = 50;
d = 1;
%生成椭圆
for m =1:180/delta;
    L=a^2*cos(m*delta)^2+b^2*sin(m*delta)^2;
    for t=-N/2+1:N/2;
        if (t^2 <= L)
            proj(t+N/2,m)=2*a*b*sqrt(L-t^2)/L;
        else
            proj(t+N/2,m)=0;
        end
    end
end
%%%%各种滤波函数
for t = -N+1:N-1;
    h1(t+N)=-2/(pi^2*d^2*(4*t^2-1));  % S-L
end;
for t=-N+1:N-1;
    if t==0
        h2(t+N)=1/(4*d^2);   % R-L
    elseif mod(t,2) == 1
        h2(t+N)=-1/(t^2*pi^2*d^2);
    else
        h2(t+N)=0;
    end;
end;
% Lewitt
esp=0.5;
for t=-N+1:N-1;
    if t==0
        h3(t+N)=(1-2*esp/3)/(4*d^2);
    elseif mod(t,2)==1
        h3(t+N)=-(1-esp)/(t^2*pi^2*d^2);
    else
        h3(t+N)=-esp/(t^2*pi^2*d^2);
    end;
end;
%存储各种反投影结果
rProj1=zeros(N,N);
rProj2=zeros(N,N);
rProj3=zeros(N,N);
rProj4=zeros(N,N);
nProj=[];

%%调用重建函数
produce('ori');
produce('sl');
produce('rl');
produce('lew');

    function produce(filter)
        for m=1:180/delta;
            creat_nProj(m);%补零，用于卷积
            
            c_h1=conv(nProj,h1);%求卷积
            c_h2=conv(nProj,h2);
            c_h3=conv(nProj,h3);
            c=0.5*(N-1)*(1-cos(m*delta)-sin(m*delta));
            for i=1:N;
                for j=1:N;
                    L(i,j)=c+(i-1)*cos(m*delta)+(j-1)*sin(m*delta); %内插
                    n = fix(L(i,j));
                    cL=L(i,j)-n;
                    if strcmp(filter,'ori')
                       if (n>0)&(n<255)
                            rProj1(i,j)=rProj1(i,j)+(1-cL)*nProj(n+N)+cL*nProj(n+1+N);
                        elseif n==255
                            rProj1(i,j) =rProj1(i,j)+nProj(n+N);
                        elseif n==0
                            rProj1(i,j) =rProj1(i,j)+nProj(n+N+1);
                       end
                    elseif strcmp(filter,'sl')
                        if (n>0)&(n<255)
                            rProj2(i,j)=rProj2(i,j)+(1 - cL)*c_h1(n+N+N-1)+cL*c_h1(n+N+N);
                        elseif n==255
                            rProj2(i,j) =rProj2(i,j)+c_h1(n+N+N-1);
                        elseif n==0
                            rProj2(i,j)=rProj2(i,j)+c_h1(n+N+N);
                       end
                    elseif strcmp(filter,'rl')
                       if (n>0)&(n<255)
                            rProj3(i,j)=rProj3(i,j)+(1 - cL)*c_h2(n+N+N-1)+cL*c_h2(n+N+N);
                        elseif n==255
                            rProj3(i,j)=rProj3(i,j)+c_h2(n+N+N-1);
                        elseif n==0
                            rProj3(i,j) =rProj3(i,j)+c_h2(n+N+N);
                       end
                    elseif strcmp(filter,'lew')
                        if (n>0)&(n<255)
                            rProj4(i,j)=rProj4(i,j)+(1 - cL)*c_h3(n+N+N-1)+cL*c_h3(n+N+N);
                        elseif n==255
                            rProj4(i,j) =rProj4(i,j)+c_h3(n+N+N-1);
                        elseif n==0
                            rProj4(i,j) =rProj4(i,j)+c_h3(n+N+N);
                        end
                    end
                end
            end
        end
    end

    function creat_nProj(m)
        for t=1:((N-1) * 3 + 1);
            if (t<N)
                nProj(t)=(proj(1,m)+proj(2,m))/2;
            elseif (t >= 2*N)
                nProj(t)=(proj(N-1,m)+proj(N,m))/2;
            else
                nProj(t)=proj(t-255,m);
            end
        end
    end

%%显示输出结果
subplot(2,2,1);
imshow(rProj1, [min(min(rProj1)), max(max(rProj1))]);
title('直接反投影');
subplot(2,2,2);
imshow(rProj2, [min(min(rProj2)), max(max(rProj2))]);
title('S-L滤波');
subplot(2,2,3);
imshow(rProj3, [min(min(rProj3)), max(max(rProj3))]);
title('R-L滤波');
subplot(2,2,4);
imshow(rProj4, [min(min(rProj4)), max(max(rProj4))]);
title('Lewitt滤波');
end