CT断层成像系列04——平行束直接反投影和滤波反投影FBP(附C++|Matlab代码
时间地点_20260201上海光源
终于到第一个重头戏了, 本节我讲CT断层成像中最火最简单最好用的重建方法滤波反投影(Filtered Backprojection, FBP)了.
然后X光束线的分布是最简单的平行束, 虽然平行束它非常简单, 但是它是有应用场景的, 没错, 就是同步辐射光源, 光源能够出射平行分布的X射线, 随着距离几乎不发散, 这在CT界比较罕见.
开门见山, 先唯象的介绍直接反投影和滤波反投影究竟是怎么成像的.
重建一个点源的图像
我们以一个点源图像为例来拆解FBP成像原理, 它足够简单, 而一般的图像就是点(\(\delta\)函数)在不同位置的堆叠而已.
下图为图1, 这是一个点源物体图像投影数据的采集过程, 由于物体只有一个点, 所以采集的投影数据是一个立柱, 也可以说是\(\delta\)函数, 采集过程中我们不同角度采集五张投影, 就会得到图二中的图(a).
下图为图2, 按图1这样对点源旋转采集后, 我么就得到了五个投影sino图(如图a), 解析法的思路很简单直接——我拍到什么投影, 我就认为你长这样, 原图的射线这一条轨道上直接全部涂抹这个点的数值(这肯定不对, 导致了后面的问题), 如下图中(b)所示, 直接匀匀一条线直接反抹回去, 那把5个角度的点的投影都反抹回去, 重建出来的图像如图(c), 所有角度都叠加起来, 我们重建出来的图像就如图(d)所示, 与点源相比, 看着是不是有点区别?
没错, 除了点源位置外被正确重构外, 图像全领域都弥漫着一层"薄雾", 这是因为我们在图(b)直接反投影的过程中, 多加了很多不该存在的点, 本来反投影只该把点源处反抹投影值的, 但是我们怎么知道这个点在哪? 在这个角度我们确实不知道点源在哪. 只能沿着射线全给反抹了——这就导致抹多了, 每一个角度都多抹了点源前后射线轨迹的点, 旋转反抹一圈, 就是除了点位置, 全图范围内都被多反抹了, 而且, 对于原图的每一个点源, 不同角度的反抹线都围绕着他们反抹, 所以越是离原图的点源靠近, 被多反摸的情况就越剧烈, 就会如图d一般, 重建图像从点源处弥散了一个"纱帐".
要想解决这个问题, 依然要解决那个核心问题——在本角度下，我们不知道点源在哪， 只能沿直线全反投影， 导致了额外的噪声。我们是知道其他角度下点源的位置的， 可以在其他角度把本角度造成的多抹给消除掉， 我们给其他角度的除点投影处外都乘上负数卷积核，这样从其他角度除了点源这一条线是正的， 其他全图位置都是负的， 就把本角度除了点源位置外的多抹的多余射线反投影都正负相消掉了。之前说过，点源附近抹的更多，所以我们给\(\delta\)函数的卷积核改造改造, 给点投影附近负的多一点, 远处的负的少一点, 这样所有角度除了原图点处不变, 点附近强烈的正负相消, 得到的图就近似于原图了. 这个卷积核就是图(e)的负值"小翅膀", 正负相消后, 得到的原图就是图(f)了. 对投影图卷积就是滤波, 所以图(f)就是先滤波后反投影FBP, 直接反投影就是图(d). 这个例子是真的非常好, 直接把FBP的重建过程几何可视化了.
中心切片定理
说完了FBP成像的原理, 我们要开始讲FBP的前置定理中心切片定理了, 为啥前置条件后讲? 因为我觉得这东西有点问题.
先上图, 这是图3(二维成像的中心切片定理). 这个定理也叫投影切片定理和傅里叶中心切片定理. 内容如下:
在探测器旋转到某个角度\(\theta\)后:
(1) 物体原图\(f(x,y)\)的投影正弦图\(p(s)\), \(p(s)\)的一维傅里叶变换为\(P(\omega)\).
(2) 物体原图\(f(x,y)\)的二维傅里叶变换图\(F(\omega_x,\omega_y)\), 取图\(F(\omega_x,\omega_y)\)中过旋转中心且与探测器平行的一条线\(P(\omega)\).
中心定理说以上两个\(P(\omega)\)是同一个.
由中心切片定理我们可知, 我们拍一个角度\(\theta\)就能获得物体原图的二维傅里叶变换图沿这个角度的切线\(F(\omega_x,\omega_y)\), 180°转半圈, 就能得到完整的傅里叶频谱图\(F(\omega_x,\omega_y)\), 然后二维逆傅里叶变换回去, 我们就得到原图了. 没错这其实就是不滤波直接反投影, 中心切片定理在给FBP提供合法性.
但是一个问题呼之欲出, 如果中心切片定理真成立, 那我还滤波干什么, 直接反投影傅里叶逆小波走起, 就以及拿下了. 是的, 没错, 中心切片定理有问题(折磨我半年多悟出来了)——那就是定理中的(1)(2)并不相等.
道理很简单, 定理的条件(2), 逆傅里叶变换完, 是物体\(f(x,y)\)过中心旋转\(\theta\)的图像, 这没毛病. 问题在左边, 左边探测器得到的是\(\theta\)角度下, 垂直的射线们累积起来的灰度值. 投影图不是仅有过中心那行的投影值(如果只有他就是相等的), 而是平行的所有行的叠加, 左侧要比右侧大的多. 没错, 这里左侧也多加了,所以中心定理只是近似,所以也需要滤波处理, 这下就通顺了.与上面讲FBP原理是反投影多抹了有异曲同工之妙.
-+
从信号与系统的角度看为何要滤波以及如何求得过滤器
这里再推荐另一本黄力宇老师的书<<医学影像的数字处理>>
这蓝皮书的第192-197页讲的直接反投影，和黄的另一本白皮书一起看，这蓝书从另一个角度讲了为什么要需要滤波——直接反投影会产生星状伪影，然后用数学公式推导出了滤波算子的最优形态。
先上图，这是本文图6， 依然是对一个密度为1的单个体素（点源）进行投影。 经过四个方向扫描正向投影后， 我们再将其直接回抹来重建物体图像。原图（a）的正方形变成了图（f）的”星“状物，星星以外，全图还弥漫着多抹的伪影。
我们从信号与背景的角度来思考问题，先看下图图6，左侧图(a)是我们的原图一个点，经过180°直接反投影后留下一堆星状伪影(图(b))， 重构后得到的原图我们把它立起来看，所谓星状物它近似是一个\(\delta\)函数(实际上不是,如果真是\(\delta\)就是完美重建系统无需滤波了), 也可以是一个冲激函数. 一个点过了一个黑箱子系统变成了一个脉冲函数, 非常符合信号系统.
如下图图7, 我们把点源信号当作某个系统的输入, 而最终重建的图像看作该系统的输出时, 根据信号与系统理论, 输出实际上是该系统的冲激响应, 这个系统就是反投影重建的变换过程. 把对原图的扫描前向投影和反投影重建当作一个系统时, 原图就可以看作系统的输入\(\mu(x,y)\), 反投影后重建的图像可以看作系统的输出\(f(x,y)\), 如下图图7:
这个扫描然后反投影重建的过程可以视做对原图每一个点的滤波函数, 经过推导(书里跳了, 略过略过哈哈哈)我们可以得知滤波核为:
那么整个成像系统就可以写为:
**表示卷积, 直接滤波重建得到的图像\(f(x,y)\), 可以视作是原图\(\mu(x,y)\)对这个等效系统的卷积操作\(h(x,y)\)的响应.
最理想的重建就是原图一个点重建完还是一个点, 现在的这个系统肯定不是. 那我们可以再追加一个滤波卷积核\(q(x,y)\), 来把重建的输出图像矫正为\(\delta\)函数. 也就是要求:
没错这就是我们先滤波后反投影的那个滤波核, 现在我们来求他的理论最完美形态, 先把卷积变成傅里叶波相乘, 这样才能等式左右两边倒腾. 等式两边傅里叶变换, 得到:
即:
也就是说, 直接反投影后, 把图像二维傅里叶变换, 再乘上\(Q(\omega_1+\omega_2)\)这个二维滤波算子, 就能得到完美的原图. 蓝书说这是个理想的二维滤波器, 即使加上近似也不容易实现.
经我师兄提醒, 我可能理解这个不容易实现是啥意思了, 首先我们傅里叶变换后的二维频谱图, 肯定只截取中间部分, 无穷远处的极高频信号不要了, 这导致了\(Q(\omega_1+\omega_2)\)的第一个不好实现, 我成为频谱"不够远(大)", 第二个不好实现是这一看就是频域上的一个圆, 计算机是离散的, 频谱图像素只要不是无限多\(Q(\omega_1+\omega_2)\)终究不够圆, 我称之为频谱"不够密". 从信号与系统的角度看, CT重建系统, 与可见光聚焦凹凸镜是极为类似的, 因此透镜不能无限大, 把点光源无穷远处的光折射聚焦回来, 导致成像点变成了爱丽光斑. CT系统由于滤波频谱图不能无限大无限密, 也导致了原图的点源重建后变成了一个类似光斑的"星"状物.
直接反投影
终于到了正题了, 正题反而非常简单. 反投影比正向投影简单的多, 可以对比系列03的正投影来看, 这里直接加速讲了, 没有新东西. 首先, 对于某个角度\(\theta_m\)时的探测器垂直接收的平行束射线\(t_m\), 它在探测器通道的位置编号可以用横纵坐标\((x,y)\)求得, 有:
这样我们就知道我们要重建的图像上的每一个像素点\((x_i，y_j)\)应该反抹哪根射线的投影值\(p(n_0+\delta)\)了，但是我们求得的射线编号几乎都是小数(\(n_0是整数部分,\delta是小数部分\))，它总是有零有整介于第\(n_0\)根射线和第\(n_0+1\)根射线之间, 因此重建的图像上的像素点\((x_i，y_j)\)需要结合这左右两根射线的投影值, 也就是需要线性插值
没了, 就这么简单, 就两步: 1算重建图像每个点要插那根射线. 2.用线性插值算这根射线的投影值, 附上一张线性插值的示意图.
剩下的折磨人脑子的就是各个坐标系的转换了, 计算机语言中的图像都是用数组实现的, 因此都是从左上角往下往右从0到N-1的计数的, 而我们在推导出的数学公式都是在笛卡尔坐标系下, 零点在图像中心, 上正下负左正右负. 而书中给出的图像计数坐标, 为了最终公式美观直观, 是从左下角开始计数的. 写代码的时候各种坐标转换绕的狠, 现在搞明白了, 过一点时间忘了又混乱了. 不过影响不大, 因为哪怕错一位, 也只是重建图像偏一个像素而已, 图像坐标转换的错误无伤大雅. 最关键的是探测器旋转的中心位置的选定, 这个错了影响图像质量, 但由于实际实验中这玩意都是自己填的, 所以代码里旋转中心大大方方\(N/2\)附近就行了, 错了也不影响, 实战里肯定不是正中间.
下面开始推导最终公式, 就是把坐标平移带入公式(6), \(x_i,y_j\)相对于矩阵计数\(i,j\)往左|负平移\((\frac{N}{2}-1)\)(书上为啥多个-1我疑惑了很久也没懂), 探测器计数\(t\)也得往正平移\(\frac{N}{2}\), 以使得for循环的时候, t=0时的值映射探测器最左侧第一个像素. 先把横纵坐标平移至矩阵计数坐标, 于是有:
然后再把探测器坐标平移到for循环迭代能够使用:
算出投影编号, 线性插值就行了, 多说无益, 下面进入咱们最关键的代码环节:
先介绍一下我偶然意外发现的宝藏Github代码库, 我发现了黄力宇老师的书对应的代码, 不仅有书里的Matlab版代码, 居然还有C++版本的代码, 非常强大(C++的直接运行不起来, 让我研究研究). 这里给出链接和截图,
https://github.com/zhengpu-lab/CT-reconstruction
先上Matlab最简单的调包代码:
clc;
close all;
clear;
I=phantom(256);
theta=0:179;
P=radon(I,theta);
rec=iradon(P,theta,'None');//反拉东iradon就是直接反投影了
figure;
imshow(I,[]),title('原始图像');
figure;
imshow(rec,[]),title('直接反投影图像');
运行结果如下:
然后上Matlab手搓版本:
%%===主程序===%%
clc;
close all;
clear;
N=256;
I=phantom(N);
delta=pi/180;
theta=0:1:179;
theta_num=length(theta);
P=radon(I,theta);
[mm,nn]=size(P);
e=floor((mm-N-1)/2+1)+1;
P=P(e:N+e-1,:);
P1=reshape(P,N,theta_num);
rec=medfuncBackprojection(theta_num,N,P1,delta);
figure;
imshow(I,[]),title('原始图像');
figure;
imshow(rec,[]),title('直接反投影重建图像');
%%===子程序===%%%
function rec=medfuncBackprojection(theta_num,N,R1,delta)
rec=zeros(N);
for m=1:theta_num
pm=R1(:,m);
Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
for k1=1:N
for k2=1:N
Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
n=floor(Xrm);
t=Xrm-floor(Xrm);
n=max(1,n);n=min(n,N-1);
p=(1-t)*pm(n)+t*pm(n+1);
rec(N+1-k1,k2)=rec(N+1-k1,k2)+p;
end
end
end
end
运行结果如下:
C++版本由于所有的函数都需要自己手搓实现, 很长, 而且和03的代码一模一样, 在本文的组后一起贴出来, 滤波不滤波只需要改一行代码
滤波反投影FBP
先滤波后反投影算法是名扬CT界的最强大的算法, 简单好用又快. 仅需要在直接反投影算法之前对投影Sino图或者之后对重建出来的图像, 进行滤波|卷积处理即可. 而这么做仅需要在前面的直接反投影算法上, 加个滤波函数即可. 因此滤波反投影法的核心之一就是滤波, 频域滤波等于空间域卷积, 白书中使用的是卷积来实现.
这一节就简单讲讲滤波核, 其他的前面都讲过了, 滤波核就讲两个, Ram-Lak(R-L)滤波函数和Shepp-Logan(S-L)滤波函数 . 下图是R-L的频域和时域的分布. 分析滤波函数为啥长这样要结合上面讲的知识(即直接投影法的缺陷), 频域(左图(a))的比较直观, 去掉低频保留高频, 因为我们反抹多的那些全图弥漫的"雾"它们从点源往外变化不激烈的逐步变低, 都是低频信号, 滤掉低频就伪影就没了. 然后从卷积特行看(右图(b)), 我们卷积是为了不同角度的点源的投影经过卷积处理后非点源处多抹的地方正负相抵消, 因此卷积核点源处不变(中间是最高峰), 这样重建出的信息没有被影响, 然后两把快速下降且反复震荡, 这样不同角度的卷积后的投影正负相消再叠加这些噪声就消干净了.
R-L滤波器的核函数就是个\(|x|\), 简单好用直观, 我们代码里用它的卷积形式:
不直观, 直接给图, 下图就是这个曲线的离散形式, 也就是我们用的版本:
公式里的\(\delta\)我寻思应该就是\(d\)的意思, 也就是一个像素, 每隔一个像素的地方离散的取卷积函数的值. 实战的时候可以无视. 用R-L滤波核的缺陷是存在吉布斯现象, 表现为较为明显的振荡响应. 为了缓解振荡响应, Shepp-Logan压低了滤波函数在边缘处的值(这块我没了解, 大致说一嘴).
S-L的频域形式, 时域形式为:
S-L的离散形式, 线性内插为
滤波反投影FBP就是在直接反投影前, 对Sino图继续滤波处理, 就搞定了, 下面开始上FBP的代码.
先附上Matlab直接调包的版本:
clc;
close all;
clear;
N=256;
I=phantom(N);
theta=0:179;
P=radon(I,theta);
rec=iradon(P,theta,'linear','None');
rec_RL=iradon(P,theta,'linear','Ram-Lak');
rec_SL=iradon(P,theta,'linear','Shepp-Logan');
figure;
subplot(2,2,1),imshow(I),title('原始图像');
subplot(2,2,2),imshow(rec,[]),title('直接反投影重建图像');
subplot(2,2,3),imshow(rec_RL,[]),title('R-L反投影重建图像');
subplot(2,2,4),imshow(rec_SL,[]),title('S-L反投影重建图像');
运行结果如下:
然后是matlab解析法手搓的FBP版本:
%%===主程序===%%
clc;
close all;
clear;
N=256;
I=phantom(N);
delta=pi/180;
theta=0:179;
theta_num=length(theta);
d=1;
P=radon(I,theta);
[mm,nn]=size(P);
e=floor((mm-N-1)/2+1)+1;
P=P(e:N+e-1,:);
P1=reshape(P,N,theta_num);
fh_RL=medfuncRlfilterfunction(N,d);
fh_SL=medfuncSlfilterfunction(N,d);
rec=medfuncBackprojection(theta_num,N,P1,delta);
rec_RL=medfuncRLfilteredbackprojection(theta_num,N,P1,delta,fh_RL);
rec_SL=medfuncSLfilteredbackprojection(theta_num,N,P1,delta,fh_SL);
figure;
subplot(2,2,1),imshow(I),xlabel('（a）256×256头模型（原始图像）');
subplot(2,2,2),imshow(rec,[]),xlabel('（b）直接反投影重建图像');
subplot(2,2,3),imshow(rec_RL,[]),xlabel('（c）R-L反投影重建图像');
subplot(2,2,4),imshow(rec_SL,[]),xlabel('（d）S-L反投影重建图像');
%%===子程序1===%%
%R-L滤波函数
function fh_RL=medfuncRlfilterfunction(N,d)
fh_RL=zeros(1,N);
for k1=1:N
fh_RL(k1)=-1/(pi*pi*((k1-N/2-1)*d)^2);
if mod(k1-N/2-1,2)==0
fh_RL(k1)=0;
end
end
fh_RL(N/2+1)=1/(4*d^2);
end
%%===子程序2===%%
%S-L滤波函数
function fh_SL=medfuncSlfilterfunction(N,d)
fh_SL=zeros(1,N);
for k1=1:N
fh_SL(k1)=-2/(pi^2*d^2*(4*(k1-N/2-1)^2-1));
end
end
%%===子程序3===%%
%R-L反投影函数
function rec_RL=medfuncRLfilteredbackprojection(theta_num,N,R1,delta,fh_RL)
rec_RL=zeros(N);
for m=1:theta_num
pm=R1(:,m);
pm_RL=conv(fh_RL,pm,'same');
Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
for k1=1:N
for k2=1:N
Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
n=floor(Xrm);
t=Xrm-floor(Xrm);
n=max(1,n);n=min(n,N-1);
p_RL=(1-t)*pm_RL(n)+t*pm_RL(n+1);
rec_RL(N+1-k1,k2)=rec_RL(N+1-k1,k2)+p_RL;
end
end
end
end
%%===子程序4===%%
%S-L反投影函数
function rec_SL=medfuncSLfilteredbackprojection(theta_num,N,R1,delta,fh_SL)
rec_SL=zeros(N);
for m=1:theta_num
pm=R1(:,m);
pm_SL=conv(fh_SL,pm,'same');
Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
for k1=1:N
for k2=1:N
Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
n=floor(Xrm);
t=Xrm-floor(Xrm);
n=max(1,n);n=min(n,N-1);
p_SL=(1-t)*pm_SL(n)+t*pm_SL(n+1);
rec_SL(N+1-k1,k2)=rec_SL(N+1-k1,k2)+p_SL;
end
end
end
end
%%===子程序5===%%
%不滤波直接反投影
function rec=medfuncBackprojection(theta_num,N,R1,delta)
rec=zeros(N);
for m=1:theta_num
pm=R1(:,m);
Cm=(N/2)*(1-cos((m-1)*delta)-sin((m-1)*delta));
for k1=1:N
for k2=1:N
Xrm=Cm+(k2-1)*cos((m-1)*delta)+(k1-1)*sin((m-1)*delta);
n=floor(Xrm);
t=Xrm-floor(Xrm);
n=max(1,n);n=min(n,N-1);
p=(1-t)*pm(n)+t*pm(n+1);
rec(N+1-k1,k2)=rec(N+1-k1,k2)+p;
end
end
end
end
我把代码文件全和一块了, 这样方便你复制直接就能跑, 运行结果如下:
C++版本的代码就在系列03里, 滤波不滤波你直接把滤波那一行注释掉就行, 给出代码:
#include <opencv2/opencv.hpp>
#include <cmath>
using namespace cv;
using namespace std;
/* -------------------------------------------------
* 1. 生成 Shepp-Logan 头模型（256×256）
* ------------------------------------------------- */
Mat phantomSheppLogan(int size = 256)
{
Mat img(size, size, CV_32F, Scalar(0));
const int c = size / 2;
const float scale = size / 2.0f;
// 标准 Shepp-Logan 参数
//椭圆绘图的六个参数 A密度|振幅 a长轴 b短轴 x0椭圆心横坐标 y0椭圆形纵坐标 θ旋转角度
struct Ellipse { float A, a, b, x, y, theta; };
const Ellipse E[] = {
{1.0f, .69f, .92f, 0.f, 0.f, 0.f}, // 大椭圆
{-.8f, .6624f, .874f, 0.f, -.0184f, 0.f}, // 内椭圆
{-.2f, .11f, .31f, .22f, 0.f, -18.f * CV_PI / 180},
{-.2f, .16f, .41f, -.22f, 0.f, 18.f * CV_PI / 180},
{.1f, .21f, .25f, 0.f, .35f, 0.f},
{.1f, .046f, .046f, 0.f, .1f, 0.f},
{.1f, .046f, .046f, 0.f, -.1f, 0.f},
{.1f, .046f, .023f, -.08f, -.605f, 0.f},
{.1f, .023f, .023f, 0.f, -.605f, 0.f},
{.1f, .023f, .046f, .06f, -.605f, 0.f}
};
for (auto& e : E)
{
float ct = cos(e.theta), st = sin(e.theta);
for (int v = 0; v < size; ++v)
{
for (int u = 0; u < size; ++u)
{
float x = (u - c) / scale;
float y = (v - c) / scale;
float xp = x * ct + y * st;
float yp = -x * st + y * ct;
//float d2 = (xp - e.x) * (xp - e.x) / (e.a * e.a) + (yp - e.y) * (yp - e.y) / (e.b * e.b);
float d2 = ((x - e.x) * ct + (y - e.y) * st) * ((x - e.x) * ct + (y - e.y) * st)
/ (e.a * e.a) + (-(x - e.x) * st + (y - e.y) * ct)* (-(x - e.x) * st + (y - e.y) * ct) / (e.b * e.b);
if (d2 <= 1.0f)
img.at<float>(v, u) += e.A;
}
}
}
normalize(img, img, 0, 1, NORM_MINMAX);
Mat img_updownflip;
flip(img, img_updownflip, 0);
return img_updownflip;
}
/* -------------------------------------------------
* 2. 平行束前向投影（Radon）
* angles: 0~179°，共 nAngles 个
* ------------------------------------------------- */
/* -------------------------------------------------
* 重写：平行束 Radon 变换（射线驱动 + 双线性插值）
* 接口：Mat radonParallel(const Mat& img, int nAngles = 360)
* ------------------------------------------------- */
Mat radonParallel(const Mat& img, int nAngles)
{
int w = img.cols, h = img.rows;
int diag = cvCeil(std::sqrt(float(w * w + h * h))); // 对角线像素数
int rBins = diag; // 探测器通道数
Mat sino(nAngles, rBins, CV_32F, Scalar(0.0f));
// 图像中心坐标（像素单位）
const float cx = (w - 1) * 0.5f;
const float cy = (h - 1) * 0.5f;
const int nSamples = 2 * diag; // 每条射线采样点数
for (int ia = 0; ia < nAngles; ++ia)
{
float theta = ia * CV_PI / nAngles;
float ct = cos(theta), st = sin(theta);
for (int ir = 0; ir < rBins; ++ir)
{
// 探测器坐标 r（像素单位，范围 [-diag/2, diag/2]）
float r = ir - (rBins - 1) * 0.5f;
// 射线中心点（在像素坐标系中）
float x_center = cx - st * r;
float y_center = cy + ct * r;
float sum = 0.0f;
// 沿射线积分，采样范围 [-diag/2, diag/2] 像素
for (int k = 0; k < nSamples; ++k)
{
// 采样偏移量 t（像素单位）
float t = (k + 0.5f) / nSamples * diag - diag * 0.5f;
// 当前采样点坐标（像素坐标系）
float x = x_center + t * ct;
float y = y_center + t * st;
// 双线性插值
if (x >= 0 && x <= w - 1 && y >= 0 && y <= h - 1)
{
int x0 = int(x), y0 = int(y);
int x1 = std::min(x0 + 1, w - 1);
int y1 = std::min(y0 + 1, h - 1);
float dx = x - x0, dy = y - y0;
float val =
(1 - dx) * (1 - dy) * img.at<float>(y0, x0) +
(dx) * (1 - dy) * img.at<float>(y0, x1) +
(1 - dx) * (dy)*img.at<float>(y1, x0) +
(dx) * (dy)*img.at<float>(y1, x1);
sum += val;
}
}
sino.at<float>(ia, ir) = sum / nSamples;
}
}
normalize(sino, sino, 0, 1, NORM_MINMAX);
return sino;
}
/* -------------------------------------------------
* 3. 滤波反投影（FBP）
* 滤波器：ramp（|ω|）
* ------------------------------------------------- */
Mat filterRamp(const Mat& sino)
{
int nAngles = sino.rows;
int rBins = sino.cols;
Mat sino_f = sino.clone();
// 1D FFT 每一行
Mat planes[] = { Mat_<float>(sino), Mat::zeros(sino.size(), CV_32F) };
Mat complex;
merge(planes, 2, complex);
dft(complex, complex, DFT_ROWS);
// 构造 ramp 滤波器
vector<float> ramp(rBins, 0);
for (int i = 0; i <= rBins / 2; ++i)
ramp[i] = i; // |ω|
for (int i = rBins / 2 + 1; i < rBins; ++i)
ramp[i] = rBins - i;
// 频域逐行相乘
for (int ia = 0; ia < nAngles; ++ia)
{
for (int k = 0; k < rBins; ++k)
{
float re = complex.at<float>(ia, 2 * k);
float im = complex.at<float>(ia, 2 * k + 1);
complex.at<float>(ia, 2 * k) = re * ramp[k];
complex.at<float>(ia, 2 * k + 1) = im * ramp[k];
}
}
// 逆 FFT
idft(complex, complex, DFT_ROWS | DFT_SCALE);
split(complex, planes);
return planes[0];
}
Mat fbp(const Mat& sino, int w, int h)
{
int nAngles = sino.rows;
int rBins = sino.cols;
Mat sino_f = filterRamp(sino);
Mat recon = Mat::zeros(h, w, CV_32F);
Point2f center(w / 2.0f, h / 2.0f);
for (int ia = 0; ia < nAngles; ++ia)
{
float theta = ia * CV_PI / nAngles;
float ct = cos(theta), st = sin(theta);
for (int u = 0; u < w; ++u)
{
for (int v = 0; v < h; ++v)
{
float x = u - center.x;
float y = v - center.y;
float r = x * ct + y * st;
//卧槽, 这里加了旋转中心,rBins / 2.0f扩展的时候在这里改
//ir探测器像素小于0直接跳过不插了,q
int ir = cvRound(r + rBins / 2.0f);
//注意, 这里是单插值, 要把小数扔进去
if (ir >= 0 && ir < rBins)
recon.at<float>(v, u) += sino_f.at<float>(ia, ir);
}
}
}
// 归一化
recon *= CV_PI / nAngles;
normalize(recon, recon, 0, 1, NORM_MINMAX);
return recon;
}
/* -------------------------------------------------
* 4. 旋转图像
* OpenCV提供的旋转比较底层 转个图需要好几步 封装一下
* ------------------------------------------------- */
Mat rotate(Mat recon)
{
// 获取图片的尺寸
int width = recon.cols;
int height = recon.rows;
// 创建一个新的Mat对象来存储旋转后的图像，宽高互换后的大小（顺时针旋转90度）
cv::Mat dst(width, height, recon.type());
// 获取旋转矩阵（顺时针旋转90度）
cv::Point2f center((width - 1) / 2.0, (height - 1) / 2.0); // 计算中心点坐标，用于旋转中心对齐图像中心点
cv::Mat rot = cv::getRotationMatrix2D(center, -90, 1.0); // 顺时针旋转90度，缩放因子为1（无缩放）
// 应用仿射变换进行旋转（顺时针旋转90度）
cv::warpAffine(recon, dst, rot, dst.size()); // 使用仿射变换函数进行旋转，结果存储在dst中
return dst;
}
/* -------------------------------------------------
* 4. main：一条龙跑通
* ------------------------------------------------- */
int main()
{
// 1. 生成 phantom
Mat phantom = phantomSheppLogan(256);
imshow("Shepp-Logan Phantom", phantom);
//waitKey(0);
// 2. 前向投影()
Mat sino = radonParallel(phantom, 180);
imshow("Sinogram (180 views)", sino);
//waitKey(0);
//验证前向投影是否正确, 需要用CT重建算法重建出切片来对比原图
// 3. FBP 重建
Mat recon = fbp(sino, phantom.cols, phantom.rows);
// 4.旋转重接结果
Mat dst = rotate(recon);
imshow("FBP Reconstruction", dst);
waitKey(0);
// 可选：写盘
// imwrite("phantom.png", phantom*255);
// imwrite("sino.png", sino*255);
// imwrite("recon.png", recon*255);
return 0;
}
运行结果如下:
本文结束, 由于博主在CT上已经做出进展了, 重心要从算法转软件工程实现了, 后面的等角扇束FBP, 等距扇束FBP, 锥束FDK, 迭代算法ART等等简单讲讲, 主要是附上代码和公式推导(也有可能先搁置, 具体看情况), 其他的看黄力宇老师书吧. 在后面的迭代算法ART和压缩感知算法等还得再往后搁置, 如果我有空, 这些肯定是要做的.
浙公网安备 33010602011771号