在MATLAB中過濾反投影和設計濾波器

我想寫我自己的MATLAB代碼來計算逆氡變換（iradon），並且到目前爲止我已經成功地使用斜坡過濾器，海明窗口重建圖像，並且還基於Kak和Shakey的教科書，使用空間域中的一維投影與我的代碼中的窗口h的卷積。然而，我認爲如果我對窗口h進行FFT並將其乘以1D投影的FFT，我應該能夠獲得相同的重構。不幸的是，我所得到的是一團糟。在MATLAB中過濾反投影和設計濾波器

function [out] = myfbp4(arg2); 


ph = phantom(); 
sino = radon(phantom,0:0.1:179); 

rho2 = 183; % max rho 
rho1 = -183; % min rho; 
step = 1; 
npos = 367; 
dtheta = 0.1; 
angles = deg2rad(0:dtheta:179); 
nproj = length(angles); 



n1 = ceil(log(npos)/log(2)); 
N = 2^n1;     % smallest power of 2 greater than npos (for fft efficiency) 
N = 1024; % for zero padding 
fs = 1; % sampling frequency 
T = 1; % sample spacing 

% grid for reconstructed image 
x1 = rho1:step:rho2; 
y1 = rho1:step:rho2; 

[yyy, xxx] = ndgrid(-y1, x1); 


% build filter h according to Kak and Shakey 
h = -floor(npos/2):T:floor(npos/2); 

for i = 1:length(h) 
    if (mod(h(i),2) == 1) 
     h(i) = -1./(h(i)^2 * pi^2 * T^2); 
    else 
     h(i) = 0; 
    end  
    h(ceil(npos/2)) = 1/(4*T^2); 
end 

%%%%%%%%%%% RAMP FILTER %%%%%%%%%%%%%%%% 
%%%%%%%%%% this is un needed when using h %% 
%%%%%%%%%% see below... %%%%%%%%%%%%%%%%%%%% 
rampFilterNum = [0:1:N/2-1 -N/2:1:-1]; 
rampFilterAbs = abs(rampFilterNum); 
rampFilter = rampFilterAbs *fs/N; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


% positions made to correspond to N point fft 
drho = (rho2 - rho1)/(npos-1); 
rho = rho1 + (0:(npos-1)) * drho; 
positions = zeros(nproj, length(rho)); 
for i = 1:nproj, 
    positions(i, :) = rho; 
end 

if (strcmp(arg2,'filter')) 
    % compute FT of h and multiply by fft projections 
    fftProj = fft(sino, N); 
    hfft = fft(h,N); 
    fftProjFiltered = bsxfun(@times, hfft', fftProj); 
    ifftProj = real(ifft(fftProjFiltered)); 
    filteredProjections = ifftProj; 
end 

if (strcmp(arg2,'conv')) 
    % make image my convolution of projections with h 
    for iproj = 1:nproj 
     sino(:, iproj) = conv(sino(:,iproj), h, 'same'); 
    end 
    filteredProjections = sino; 
end 


% display the image through backprojection 
fdata = zeros(size(xxx)); 
for iproj = 1:nproj 
    theta = angles(iproj); 
    rho1 = xxx*cos(theta) + yyy*sin(theta); % rotate coordinate system by theta 
    %r = x1; 
    r = positions(iproj,:); 

    fdata1 = filteredProjections(1:npos,iproj); % filtered projections 
    %fdata1 interp1(
    fdata2 = interp1(r, fdata1, rho1, 'linear', 0); 
    fdata = fdata + deg2rad(dtheta) * fdata2; %theta*fdata2; 
end 

out = fdata; 
end

剛用完= myfbp4（'conv'）或myfbp4（'filter'）會顯示不同的結果。看起來卷積工作正常，但過濾方法沒有按照我所希望的那樣工作。

任何人都可以看到問題嗎？（道歉，如果有任何多餘的代碼，我試圖削減它的大部分...我也應該提到，這個代碼是從某處借來的，並修改了一下，但我不記得我在哪裏找到它）。

在此先感謝

編輯：問題解決。問題在於我沒有取窗口h的傅立葉變換的絕對值來獲得頻率窗口。對於那些發現這一點的人來說，hfft = abs（fft（h，N））應該代替hfft = fft（h，N）。

來源

2016-04-05 Billy

您應該自己回答並將答案標記爲已接受，以便大家可以看到問題已解決。 http://stackoverflow.com/help/self-answer – Tapio

問題已解決。問題在於我沒有取窗口h的傅立葉變換的絕對值來獲得頻率窗口。對於那些發現這一點的人來說，hfft = abs（fft（h，N））應該代替hfft = fft（h，N）。

來源

2017-01-26 20:36:35 Billy

在MATLAB中過濾反投影和設計濾波器

回答

相關問題