MATLAB'da filtrelenmiş geri yansıtma ve filtre tasarlama

Question

Ters radon dönüşümü (iradon) hesaplamak için kendi MATLAB kodumu yazmaya çalışıyorum ve şu ana kadar bir rampa filtresi, bir görüntü penceresi kullanarak bir görüntüyü başarılı bir şekilde yeniden oluşturmayı başardım ve Ayrıca Kak ve Shakey tarafından ders kitabına dayanan kodumda bir pencere h ile uzamsal etki alanındaki 1D projeksiyonlarının evrişimini kullanarak. Ancak, sanırım pencerenin FFT'sini alıp 1D projeksiyonlarının FFT'si ile çarptığımda, aynı rekonstrüksiyonu elde edebilmem gerekir. Ne yazık ki, yaptığım şey bir karmaşa.

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 

Sadece farklı sonuçlar gösterecektir = myfbp4 ('dönüşüm') veya myfbp4 ('filtre') azalıyor. Konvolüsyon iyi çalışıyor gibi görünüyor, ama filtreleme yaklaşımı umduğum gibi çalışmıyor.

Problemi gören var mı? (Herhangi bir yedek kod varsa özür dilerim, çoğunu kesmeye çalıştım ... Ayrıca bu kodun bir yerden ödünç aldığını ve biraz değiştirildiğini de belirtmeliyim, ama nerede bulduğumu hatırlamıyorum). peşin

Teşekkür

DÜZENLEME: Sorun çözüldü. Sorun, frekans penceresini elde etmek için pencerenin (h) fourier dönüşümünün mutlak değerini almamamdı. Bunu bulanlar için, hfft = abs (fft (h, N)) hfft = fft (h, N) ile değiştirilmelidir.

Answer 1

Sorun çözüldü. Sorun, frekans penceresini elde etmek için pencerenin (h) fourier dönüşümünün mutlak değerini almamamdı. Bunu bulanlar için, hfft = abs (fft (h, N)) hfft = fft (h, N) ile değiştirilmelidir.