Отфильтрованная обратная проекция в MATLAB и проектирование фильтра

Question

Я пытаюсь написать свой собственный код MATLAB для вычисления обратного преобразования Радона (iradon), и до сих пор мне удалось успешно восстановить изображение, используя линейный фильтр, окно Хэмминга, а также используя свертку 1D проекций в пространственном домен с окном h в моем коде на основе учебника как и шеки. Тем не менее, я думаю, что если я возьму БПФ окна h и умножу его на БПФ 1D-проекций, я смогу получить такую ​​же реконструкцию. К сожалению, я получаю беспорядок.

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).

Answer 1

Задача решена. Проблема заключалась в том, что я не брал абсолютное значение преобразования Фурье окна h, чтобы получить окно частот. Для тех, кто находит это, hfft = abs(fft(h,N)) должен заменить hfft = fft(h, N).