滤波反投影fbp算法matlab

%%投影 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') %%