%%投影
I = phantom('Modified Shepp-Logan',256);
nmax=ceil(length(I)*(sqrt(2)));
I1=zeros(512);
I1(128:383,128:383)=I;
proj=zeros(512,180);
for theta=1:180
I2=imrotate(I1,-theta,'bilinear');%原图旋转theta度，得到I2
temp=length(I2);
temp=round(temp/2-255);
I3=I2(temp:temp+511,temp:temp+511);
b=sum(I3);%b是投影向量
proj(:,theta)=b';%赋值投影矩阵
end
imshow(proj,[])
%%投影
% width = 2^nextpow2(size(I,1));
width=512;
H = 2*[0:(width/2-1), width/2:-1:1]'/width;%R-L滤波
proj_filtered = zeros(length(proj),180);
% plot(H)
proj_fft = fft(proj, 512);
for i = 1:180
proj_filtered(:,i) = proj_fft(:,i).*H;
end
% M=256;
% fbp = zeros(M); % 假设初值为0
proj_ifft = real(ifft(proj_filtered));
% imshow(proj_ifft,[])
fbp = zeros(256);
% 5. back-projection to the x- and y- axis
for i = 1:180
% rad is the angle of the projection line , not projection angle
rad = i*pi/180;
for x = (-256/2+1):256/2
for y = (-256/2+1):256/2
t = round(x*cos(rad+pi/2)+y*sin(rad+pi/2));
fbp(x+256/2,y+256/2)=fbp(x+256/2,y+256/2)+proj_ifft(t+256,i)/180;
end
end
% imshow(fbp,[])
end
fbp = fbp/360;
subplot(1, 2, 1), imshow(I), title('Original')
subplot(1, 2, 2), imshow(fbp,[]), title('FBP')
% I=dicomread('C:\Users\56326\Desktop\dcm\100.dcm');
% imshow(I,[])
% info= dicominfo('C:\Users\56326\Desktop\dcm\100.dcm');
% I1=imread('C:\Users\56326\Desktop\sonoscape\12345\20180425_145827\20180425_150454_1.jpg');
% I2=imread('C:\Users\56326\Desktop\sonoscape\12345\20180425_145827\20180425_150545_2.jpg');
% I3=imread('C:\Users\56326\Desktop\sonoscape\12345\20180425_145827\20180425_150603_3.jpg');
%
% figure,
%%
theta = 1:180;
P=phantom(128);
% 1. projection using radon function
[R,xp] = radon(P,theta);
width = 2^nextpow2(size(R,1)); % set width for fft transformation
% 2. do fft to the projection
proj_fft = fft(R, width);
% 3. filter
% Ramp filter function from 0 to width then to 0
filter = 2*[0:(width/2-1), width/2:-1:1]'/width;
proj_filtered = zeros(width,180);
for i = 1:180
proj_filtered(:,i) = proj_fft(:,i).*filter;
end
% 4. do ifft to the filtered projection
proj_ifft = real(ifft(proj_filtered)); % get the real part of the result
% 5. back-projection to the x- and y- axis
fbp = zeros(128); % asign the original value 0
for i = 1:180
% rad is the angle of the projection line , not projection angle
rad = i*pi/180;
for x = (-128/2+1):128/2
for y = (-128/2+1):128/2
t = round(x*cos(rad+pi/2)+y*sin(rad+pi/2));
fbp(x+128/2,y+128/2)=fbp(x+128/2,y+128/2)+proj_ifft(t+round(size(R,1)/2),i);
end
end
imshow(fbp,[])
end
fbp = fbp/180;
% 6. show the result
subplot(1, 2, 1), imshow(P), title('Original')
subplot(1, 2, 2), imshow(fbp), title('FBP')
%%
