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

我正在研究一种使用中心切片定理的滤波反投影算法作为家庭作业，虽然我理解纸上的理论，但在 Matlab 中实现它时遇到了问题。我得到了一个可以遵循的框架，但我认为我可能误解了一个步骤。这是我所拥有的：

function img = sampleFBP(sino,angs)

% 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');

% 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
    % 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(floor(diagDim/2), :) = fourierCurRow;
    imContrib = imContrib * 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)));

我输入的正弦图只是 Shepp-Logan 体模上 0 到 179 度的氡气函数的输出。现在运行代码会给我一个黑色图像。我认为我在将行的 FT 添加到图像的循环中遗漏了一些内容。根据我对中心切片定理的理解，我认为应该发生的是：

• 初始化一个与 2DFT 大小相同的数组（即 diagDim x diagDim）。这就是傅里叶空间。

• 取一行正弦图，对应于单个角度的线积分信息，并对其应用一维 FT

• 根据中心切片定理，该线积分的 FT 是一条穿过傅立叶域的线，该线以与投影的角度相对应的角度穿过原点。因此，为了模拟这一点，我采用该线积分的 FT 并将其放置在我创建的 diagDim x diagDim 矩阵的中心行中

• 接下来，我将创建的一维斜坡滤波器的 FT 与线积分的 FT 相乘。傅立叶域中的乘法相当于空间域中的卷积，因此这将线积分与滤波器进行卷积。

• 现在，我将整个矩阵旋转投影拍摄的角度。这应该给我一个 diagDim x diagDim 矩阵，其中一条信息以一定角度穿过中心。 Matlab在旋转时会增加矩阵的大小，但由于正弦图在开始时被填充，因此不会丢失任何信息，并且矩阵仍然可以添加

• 如果所有这些中心有一条线的空矩阵加在一起，它应该给我图像的完整 2D FT。所需要做的就是进行逆 2D FT，原始图像应该就是结果。

如果我遇到的问题是概念性的，如果有人能指出我搞砸的地方，我将不胜感激。相反，如果这是一个 Matlab 的事情（我对 Matlab 还是个新手），我会很高兴了解我错过了什么。

您发布的代码是滤波反投影（FBP）的一个很好的例子，我相信对于想要学习 FBP 基础的人来说可能很有用。可以使用该功能iradon(...)在 MATLAB 中（参见here https://www.mathworks.com/help/images/ref/iradon.html）使用各种滤波器执行 FBP。当然，就您而言，重点是学习中心切片定理的基础，因此寻找捷径并不是重点。通过回答你的问题我也学到了很多，刷新了我的知识！

现在您的代码已得到完美注释并描述了需要采取的步骤。有一些微妙的[编程]问题需要修复，以便代码能够正常工作。

首先，傅里叶域中的图像表示可能最终会由于以下原因而丢失数组：floor(diagDim/2)取决于正弦图的大小。我会把它改成round(diagDim/2)拥有完整的数据集fimg。请注意，如果处理不当，这可能会导致某些正弦图尺寸出现错误。我鼓励你想象fimg了解缺失的数组是什么以及它为何重要。

第二个问题是您的正弦图需要转置以与您的算法保持一致。因此添加sino = sino'。再次，我鼓励您尝试不使用此代码的代码，看看会发生什么！请注意，必须沿视图进行零填充，以避免由于锯齿而产生伪像。我将在这个答案中演示一个例子。

第三，也是最重要的一点，imContrib是一个数组的临时持有者fimg。因此，它必须保持与fimg, so

imContrib = imContrib * fft(rampFilter_1d);

应替换为

imContrib(floor(diagDim/2), :) = imContrib(floor(diagDim/2), :)' .* rampFilter_1d;

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

现在是一个完整的示例来展示它如何使用改进的谢普-洛根体模 https://www.mathworks.com/help/images/ref/phantom.html#d120e206667:

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');):

请注意带零填充和不带零填充的重建图像中不同的对象大小。当正弦图用零填充时，对象会收缩，因为径向内容被压缩。