2

我正在研究使用中心切片定理进行家庭作业的过滤反投影算法，虽然我理解了纸上的理论，但我在 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 并将其放在我创建的 diagDim x diagDim 矩阵的中心行

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

  • 现在我将整个矩阵旋转投影的角度。这应该给我一个 diagDim x diagDim 矩阵，其中单行信息以一定角度穿过中心。Matlab在旋转时增加了矩阵的大小，但由于正弦图在开始时被填充，没有信息丢失，矩阵仍然可以添加

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

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

4

1 回答 1

1

您发布的代码是过滤反投影 (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 纠正此错误）。因此，您应该将其fft放入，fft(rampFilter_1d)因为此过滤器应用于频域（请记住fft(x)将 x 的域，例如时间、空间等分解为其频率内容）。

现在一个完整的例子来展示它是如何使用修改后的 Shepp-Logan phantom工作的：

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');），这是相同的图像：

没有零填充正弦图的重建图像

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

于 2020-04-17T16:33:18.650 回答