cljradon.m

医学仪器课程做的CT成像原理——radon变换

P = zeros(100,100);
for x=[1:1:100]
    for y=[1:1:100]
        if((x-50)^2+(y-50)^2<=1225)
          P(x,y) = 0.5;
        end
    end
end
for x=[1:1:100]
    for y=[1:1:100]
        if((x-35)^2+(y-35)^2<=100)
           P(x,y) = 1;
        end
    end
end
theta = 0:179;
[R,xp] = radon(P,theta);
subplot(2,3,1);mesh(P);title('initial image');
subplot(2,3,2);plot(R(:,1));title('\theta=0(degrees)');
subplot(2,3,3);plot(R(:,45));title('\theta=45(degrees)');
subplot(2,3,4);plot(R(:,90));title('\theta=90(degrees)');
subplot(2,3,5);plot(R(:,135));title('\theta=135(degrees)');
subplot(2,3,6);plot(R(:,179));title('\theta=180(degrees)');

figure;subplot(2,3,1);imshow(P);title('initial image');
theta1 = 0:10:170; [R1,xp] = radon(P,theta1);
theta2 = 0:5:175;  [R2,xp] = radon(P,theta2);
theta3 = 0:2:178;  [R3,xp] = radon(P,theta3);
theta4 = 0:1:179;  [R4,xp] = radon(P,theta4);
subplot(2,3,2); imagesc(theta4,xp,R4); colormap(hot); colorbar
xlabel('\theta'); ylabel('x\prime');title('projection');
I1 = iradon(R1,10);
subplot(2,3,3); imshow(I1);title('iradon D-theta=10');
I2 = iradon(R2,5);
subplot(2,3,4);imshow(I2);title('iradon D-theta=5');
I3 = iradon(R3,2);
subplot(2,3,5); imshow(I3);title('iradon D-theta=2');
I4 = iradon(R4,1);
subplot(2,3,6);imshow(I4);title('iradon D-theta=1');