CT 滤波反投影重建 直接反投影重建

%% Andrew ID:jiag1---------------------------guojia514@gmail.com %% %% ---------filtered-backprojection and backprojection---------- %% function [I_recon_lin, I_recon_nn,I_recon_lin1, I_recon_nn1]=BP(IMG,M) close all; I=IMG; [rows cols]=size(I); theta=creat_THETA(M); [R,xp] = radon(I,theta); % Compute (forward ) projections [N, K] = size(R); % K is number of angles % N is number of detectors xp_offset = abs(min(xp) + 1); % Fourier Transform of the projections fft_N = 512 ; % Nearest power of 2 corresponding to 367 is 512 R_fft = fft(R, fft_N); % Define the ramp filter in frequency ramlak = [2 * [0:(fft_N/2-1), fft_N/2:-1:1 ]' / fft_N ]; plot(ramlak); % Filtering of the projections in frequency [m,n]=size(ramlak); for i = 1:K R_filtered(:,i) = R_fft(:,i) .* ramlak; % R_unfiltered(:,i) = R_fft(:,i)./m; R_unfiltered(:,i) = R_fft(:,i); end % Inverse Fourier Transform of the filtered projections R_ifft1 = real(ifft(R_filtered)); % Inverse Fourier Transform of the unfiltered projections R_ifft2 = real(ifft(R_unfiltered)); % the filtered backprojections I_recon_lin = zeros(rows,cols); I_recon_nn = zeros(rows,cols); for(i = 1:K) Q = R_ifft1(:,i); theta_rad = theta(i) * pi / 180; for x = 1 : cols for y = 1 : rows t_temp = (x-128) * cos(theta_rad) - (y-128) * sin(theta_rad) + xp_offset ; % Nearest Neighbour Interpolation t = round(t_temp) ; wt_nn = Q(t); I_recon_nn(y,x) = I_recon_nn(y,x) + wt_nn; % Linear Interpolation t_next = ceil(t_temp); t_prev = floor(t_temp); if ( t_prev == 0 || t_next == 0 ) wt_lin = Q(t_next); else wt_lin = ( Q(t_next) * (t_temp - t_prev) ) + ( Q(t_prev) * (t_next - t_temp)); end I_recon_lin(y,x) = I_recon_lin(y,x) + wt_lin; end end end I_recon_lin = (pi / K) * I_recon_lin ; I_recon_nn = (pi / K) * I_recon_nn ; % the unfiltered backprojections I_recon_lin1 = zeros(rows,cols); I_recon_nn1 = zeros(rows,cols); for(i = 1:K) Q = R_ifft2(:,i); theta_rad = theta(i) * pi / 180; for x = 1 : cols for y = 1 : rows t_temp = (x-128) * cos(theta_rad) - (y-128) * sin(theta_rad) + xp_offset ; % Nearest Neighbour Interpolation t = round(t_temp) ; wt_nn = Q(t); I_recon_nn1(y,x) = I_recon_nn1(y,x) + wt_nn; % Linear Interpolation t_next = ceil(t_temp); t_prev = floor(t_temp); if ( t_prev == 0 || t_next == 0 ) wt_lin = Q(t_next); else wt_lin = ( Q(t_next) * (t_temp - t_prev) ) + ( Q(t_prev) * (t_next - t_temp)); end I_recon_lin1(y,x) = I_recon_lin1(y,x) + wt_lin; end end end I_recon_lin1 = (pi / K) * I_recon_lin1 ; I_recon_nn1 = (pi / K) * I_recon_nn1 ;