image processing - 在Matlab中使用中心切片定理实现滤波反投影算法

您发布的代码是过滤反投影（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;

请注意，斜坡滤波器在频域中是线性的（感谢@Cris Luengo纠正了此错误）。因此，当在频域中应用此滤波器时，您应该插入fftin fft(rampFilter_1d)（请记住fft(x)将x的域分解为频率内容，例如时间，空间等）。

现在是一个完整的示例，以显示使用修改后的Shepp-Logan幻像的工作方式：

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');），则这是同一张图片：

请注意，在有和没有零填充的情况下，重建图像中的对象大小不同。当正弦图填充为零时，对象会收缩，因为会压缩径向内容。