在 Matlab 中使用中心切片定理实现滤波反投影算法
Implementing a filtered backprojection algorithm using the central slice theorem in 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）。这是傅立叶 space.
从单角度取一行对应线积分信息的正弦图，对其应用一维傅立叶变换
根据中心切片定理，这个线积分的傅里叶变换是一条通过傅里叶域的线，它以与投影的角度相对应的角度穿过原点。因此，为了模拟这一点，我采用该线积分的 FT，并将其放在我创建的 diagDim x diagDim 矩阵的中心行
接下来我将我创建的一维斜坡滤波器的 FT 与线积分的 FT 相乘。傅立叶域中的乘法等同于空间域中的卷积，因此这对线积分与滤波器进行了卷积。
现在我按照投影的角度旋转整个矩阵。这应该给我一个 diagDim x diagDim 矩阵，其中一行信息以一定角度穿过中心。 matlab在旋转的时候增加了矩阵的大小但是由于一开始正弦图被填充了，所以没有信息丢失，矩阵仍然可以相加
如果将所有这些具有一条穿过中心的单线的空矩阵加在一起，应该可以得到图像的完整 2D FT。只需要做逆二维傅立叶变换，原图就是结果。
如果我 运行 遇到的问题是概念性的，如果有人能指出我搞砸的地方，我将不胜感激。相反，如果这是 Matlab 的东西（我对 Matlab 还是个新手），我会很高兴了解我错过了什么。
您发布的代码是过滤反投影 (FBP) 的一个很好的示例，我相信对于想要学习 FBP 基础的人可能会有用。可以使用 MATLAB 中的函数 iradon(...)（参见 here）来使用各种滤波器执行 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;
请注意，Ramp 滤波器在频域中是线性的（感谢@Cris Luengo 纠正了这个错误）。因此，您应该在 fft(rampFilter_1d) 中删除 fft，因为此滤波器应用于频域（记住 fft(x) 分解 x 的域，例如时间、space 等到它的频率内容）。
现在一个完整的例子来展示它是如何使用 modified 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');），这里是相同的图像：
注意重建图像中使用和不使用零填充的对象大小不同。当正弦图被零填充时对象收缩，因为径向内容被压缩。
我正在研究一种使用中心切片定理的过滤反投影算法来完成家庭作业，虽然我理解纸上的理论，但我 运行 遇到了在 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）。这是傅立叶 space.
从单角度取一行对应线积分信息的正弦图，对其应用一维傅立叶变换
根据中心切片定理，这个线积分的傅里叶变换是一条通过傅里叶域的线，它以与投影的角度相对应的角度穿过原点。因此，为了模拟这一点，我采用该线积分的 FT，并将其放在我创建的 diagDim x diagDim 矩阵的中心行
接下来我将我创建的一维斜坡滤波器的 FT 与线积分的 FT 相乘。傅立叶域中的乘法等同于空间域中的卷积，因此这对线积分与滤波器进行了卷积。
现在我按照投影的角度旋转整个矩阵。这应该给我一个 diagDim x diagDim 矩阵，其中一行信息以一定角度穿过中心。 matlab在旋转的时候增加了矩阵的大小但是由于一开始正弦图被填充了，所以没有信息丢失，矩阵仍然可以相加
如果将所有这些具有一条穿过中心的单线的空矩阵加在一起，应该可以得到图像的完整 2D FT。只需要做逆二维傅立叶变换，原图就是结果。
如果我 运行 遇到的问题是概念性的，如果有人能指出我搞砸的地方，我将不胜感激。相反，如果这是 Matlab 的东西（我对 Matlab 还是个新手），我会很高兴了解我错过了什么。
您发布的代码是过滤反投影 (FBP) 的一个很好的示例，我相信对于想要学习 FBP 基础的人可能会有用。可以使用 MATLAB 中的函数 iradon(...)（参见 here）来使用各种滤波器执行 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;
请注意，Ramp 滤波器在频域中是线性的（感谢@Cris Luengo 纠正了这个错误）。因此，您应该在 fft(rampFilter_1d) 中删除 fft，因为此滤波器应用于频域（记住 fft(x) 分解 x 的域，例如时间、space 等到它的频率内容）。
现在一个完整的例子来展示它是如何使用 modified 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');），这里是相同的图像：
注意重建图像中使用和不使用零填充的对象大小不同。当正弦图被零填充时对象收缩，因为径向内容被压缩。