Код, который вы опубликовали, является довольно хорошим примером отфильтрованной обратной проекции (FBP), и я считаю, что он может быть полезен для людей, которые хотят изучить основы FBP. Можно использовать функцию iradon(...) в MATLAB (см. здесь ) для выполнения FBP с использованием различных фильтров. В вашем случае, конечно, смысл состоит в том, чтобы изучить основы теоремы о центральном срезе, и поэтому поиск кратчайшего пути - не главное. Я также многому научился и обновил свои знания, ответив на ваш вопрос!
Теперь ваш код был отлично прокомментирован и описывает шаги, которые необходимо предпринять. Есть пара тонких проблем [программирования], которые нужно исправить, чтобы код работал нормально.
Во-первых, у вашего представления изображения в области Фурье может быть отсутствующий массив из-за floor(diagDim/2) в зависимости от размера синограммы. Я бы изменил это на round(diagDim/2), чтобы иметь полный набор данных в fimg. Имейте в виду, что это может привести к ошибке для определенных размеров синограммы, если не обрабатывается правильно. Я бы посоветовал вам визуализировать fimg, чтобы понять, что это за отсутствующий массив и почему он имеет значение.
Второй вопрос заключается в том, что ваша синограмма должна быть транспонирована для соответствия вашему алгоритму. Отсюда добавление sino = sino'. Опять же, я призываю вас попробовать код без этого, чтобы увидеть, что происходит! Обратите внимание, что заполнение нулями должно происходить вдоль видов, чтобы избежать появления артефактов из-за наложения. Я продемонстрирую пример этого в этом ответе.
В-третьих, и самое главное, imContrib является временным держателем для массива вдоль fimg. Следовательно, он должен поддерживать тот же размер, что и fimg, поэтому
imContrib = imContrib * fft(rampFilter_1d);
следует заменить на
imContrib(floor(diagDim/2), :) = imContrib(floor(diagDim/2), :)' .* rampFilter_1d;
Обратите внимание, что фильтр линейного изменения линейен в частотной области (благодаря @ Крис Луен go за исправление этой ошибки). Следовательно, вы должны сбросить fft в fft(rampFilter_1d), поскольку этот фильтр применяется в частотной области (помните, fft(x) разбивает область x, такую как время, пространство и т. Д. c, на ее частотное содержимое).
Теперь полный пример, демонстрирующий, как это работает с использованием модифицированного фантома Шепп-Логана :
angs = 0:359; % angles of rotation 0, 1, 2... 359
init_img = phantom('Modified Shepp-Logan', 100); % Initial image 2D [100 x 100]
sino = radon(init_img, angs); % Create a sinogram using radon transform
% Here is your function ....
% This step is necessary so that frequency information is preserved: we pad
% the sinogram with zeros so that this is ensured.
sino = padarray(sino, floor(size(sino,1)/2), 'both');
% Rotate the sinogram 90-degree to be compatible with your codes definition of view and radial positions
% dim 1 -> view
% dim 2 -> Radial position
sino = sino';
% diagDim should be the length of an individual row of the sinogram - don't
% hardcode this!
diagDim = size(sino, 2);
% The 2DFT (2D Fourier transform) of our image will start as a matrix of
% all zeros.
fimg = zeros(diagDim);
% Design your 1-d ramp filter.
rampFilter_1d = abs(linspace(-1, 1, diagDim))';
rowIndex = 1;
for nn = angs
% fprintf('rowIndex = %g => nn = %g\n', rowIndex, nn);
% Each contribution to the image's 2DFT will also begin as all zero.
imContrib = zeros(diagDim);
% Get the current row of the sinogram - use rowIndex.
curRow = sino(rowIndex,:);
% Take the 1D Fourier transform the current row - be careful, as it's
% necessary to perform ifftshift and fftshift as Matlab tends to
% place zero-frequency components of a spectrum at the edges.
fourierCurRow = fftshift(fft(ifftshift(curRow)));
% Place the Fourier-transformed sino row and place it at the center of
% the next image contribution. Add the ramp filter in Fourier domain.
imContrib(round(diagDim/2), :) = fourierCurRow;
imContrib(round(diagDim/2), :) = imContrib(round(diagDim/2), :)' .* rampFilter_1d; % <-- NOT fft(rampFilter_1d)
% Rotate the current image contribution to be at the correct angle on
% the 2D Fourier-space image.
imContrib = imrotate(imContrib, nn, 'crop');
% Add the current image contribution to the running representation of
% the image in Fourier space!
fimg = fimg + imContrib;
rowIndex = rowIndex + 1;
end
% Finally, just take the inverse 2D Fourier transform of the image! Don't
% forget - you may need an fftshift or ifftshift here.
rcon = fftshift(ifft2(ifftshift(fimg)));
Обратите внимание, что ваше изображение имеет сложное значение. Итак, я использую imshow(abs(rcon),[]), чтобы показать изображение. Пара полезных изображений (пища для размышлений) с окончательным восстановленным изображением rcon:
А вот то же изображение, если вы закомментируете шаг заполнения нулями (т. е. комментарий sino = padarray(sino, floor(size(sino,1)/2), 'both');):
Обратите внимание на различный размер объекта в восстановленных изображениях с заполнением нулями и без него. Объект сжимается, когда синограмма заполнена нулями, поскольку радиальное содержимое сжато.