CT成像——等距扇束正反投影的简单介绍以及在matlab上的实现

Messiahx

已于 2022-07-10 21:42:10 修改

阅读量3.9k

点赞数 14

文章标签： matlab 图像处理健康医疗

于 2022-07-10 21:32:33 首次发布

本文链接：https://blog.csdn.net/Messiahx/article/details/125709097

版权

相信会点进来看这篇文章的朋友应该也是在学习断层成像方面的吧？

我接触CT成像和matlab不久，大概也就一个月，写出这些感觉还是充满了困难，在csdn上找，CT成像的相关且不用下载的代码也不多，想到学习路上的辛苦以及csdn论坛在深度学习等方向上给我的许多帮助，这里码出这一篇文章，希望能帮到需要的朋友。

为了预防别的方向的朋友点进来这里还是稍微介绍一下，CT机和X光这两个名词一定听说过吧。简单来说就是一个X光射线源发射X光，假定发射的X光是单能光(实际上并不是)，相同的物质拥有相同的衰减系数，不同的物质衰减系数往往不一样。这个衰减系数跟X光穿过一定厚度物体之后减少的光子数相关，在假定单能光的情况下会有:

$I_{out}=I_{in}\, \, *\, \, e^{-\int \mu_{s}ds }$ ，其中I是X光的强度，跟光子数有关，μ是衰减系数，一般分析时认为对于某一类物质，该物质的μ是不变的(实际也不是，μ跟光子能量有关，而一般射线源射出的是多能光，因此对于同一射线源射出的X光，同一物质的衰减系数是不同的)，s是穿过物体的距离。非常符合直觉的，从公式可以看出，X光的衰减跟穿过物体的距离(或者说厚度)以及物体的衰减系数有关。换句话说，X光的衰减里面包含了该条光方向物体厚度的信息，以及该物体材质的信息。一条光包含了一个方向的信息，那么朝物体的某一个平面发射出足够多角度的光将该平面从各个方向穿透，得到各方向厚度和材质的信息是不是就能够将该物体的平面重建出来了呢？

这就是所谓的断层成像，或者说CT成像。

前面不太明白，或者说想不太出来那个画面的可以去看看拉东变换和拉东逆变换，没接触过但是有兴趣的可以去看看傅里叶切片定理，该定理说明了根据足够多条不同角度的拉东变换可以很好地重建原图像。

回到现在的问题，等距扇束投影和反投影。

上图中的s是射线源，发射出扇形的射线穿过物体，被D2线上等距离s分布的检测器接收，其中O是物体中心，物体中心O到射线源s的距离为D，很显然各探测器与射线源相交的直线在过D'2线上的各点也是等距的，射线源与y轴的夹角为β(默认逆时针是正方向)，射线源过物体中心的线与扇形射线的夹角为γ，θ是射线垂线与x轴的夹角，t是射线垂线与原点的距离(求θ和t是为了将射线写成xcosθ+ysinθ=t的形式)。

想象射线源与探测器绕着物体中心逆时针旋转(β不断增大)得到各个角度射线的过程，这就是正投影的流程。

经过一些不太困难的数学分析我们可以发现：θ=β+γ，t=Dsinγ，因为是等距而不是等角(每一个相邻探测器夹角相等，也就是γ相等)扇形束，用s(这里的s是D'2上的s，并不是真正的探测器间距)表示θ和t。分析可以知道：

$\theta =\beta +arctan(\frac{s}{D})$

$t=\frac{Ds}{\sqrt{D^{2}+s^{2}}}$

到此重要的参数介绍完成，剩余的反投影推导过程有兴趣的可以自行推导或者自行搜索，大致的推导过程是根据平行束得到的反投影公式(具体搜索平行束反投影和傅里叶切片定理)，做极坐标转换，再将原式中对t和θ的积分根据上面给出的转换变成对s和β的积分，将所得到式子不断变形最后得到带 $h(t)=\int_{-\infty }^{+\infty}\left |w \right |e^{j2\pi wt }dw$ 的形式，也就是斜坡滤波在时域的表达式，最后就可以将积分写成带卷积的形式，得到：

其中：

因此反投影大致上可以分为三步：

首先将投影数据乘以权重 $\frac{D}{\sqrt{D^{2}+\sqrt{s^{2}}}}$ ，再对得到数据用斜坡滤波器h(s)做卷积，最后通过 $f(x,y)=\triangle \beta\sum \frac{1}{2U^{2}(x,y,\beta _{i}))}Q_{\beta _{1}}(s^{'})$ (连续域上的积分体现在离散域上就是累加)，得到最终重建结果，其中的累加的含义是射线源所有旋转角度穿过点(x,y)射线的累加，s'表示在β=βi下过点(x,y)和射线源的直线在D2’上的交点与O点连接的长度，根据这个长度插值得到投影值(因为仿真图片的像素是离散的，而显然并不是所有直线都能正好地穿过离散的像素点，不如说大部分都是穿过像素间)。

下面是代码部分：

clc;
clear;

tic
img=phantom(256);
beta=0:1:359;%射线源与探测器旋转角度，默认逆时针旋转
[M,N]=size(img);
dS=1;%探测器在物体中心线处的间隔（探测器实际上间隔会大于dS）
rayNum=ceil(sqrt(2)*M);%探测器数=射线数
distance=2*M;%射线源到物体中心的距离，SOD

R=Projection(img,beta,distance,rayNum,dS);%正投影，因为用了interp2比较慢
recon=FanBeamFBP(M,N,beta,R,rayNum,dS,distance);%滤波反投影

subplot(121);
imshow(R,[]);
subplot(122);
imshow(recon,[]);
toc

function u=U(D,x,y,beta)%公式中的U
u=(D+x*sin(beta)-y*cos(beta))/D;
end

function proj = Projection(img, beta,distance,rayNum,dS)%正投影过程
proj = zeros(rayNum, length(beta));
[M,~]=size(img);
M2=(M-1)/2;
[X, Y] = meshgrid(-M2 : 1 : M2);%图片坐标系重定位，以图像中心为(0,0)
ray_x=-(rayNum-1)/2 * dS:dS:(rayNum-1)/2 * dS;%各光线在物体中心坐标系中与射线源中心线垂直线上的各值，dS为各探测器在该线上的间隔
CY=-M2*sqrt(2):1:M2*sqrt(2);%垂直x轴线采样点的y,图像为正方形且采样间隔为1，因此设采样点数为根号2*M+1
CX=zeros([1,length(CY)]);%垂直x轴线采样点的x
for i=beta%360角全投影,默认逆时针旋转
proj_1=zeros(rayNum,1);%一个角度的投影
r=deg2rad(i);%角度转弧度
CX_1=CX.*cos(-r)+CY.*sin(-r);
CY_1=CY.*cos(-r)-CX.*sin(-r);%射线源旋转i度后采样点跟着旋转
for j=1:length(ray_x)%每个角度有ray_x根线
gamma=atan(ray_x(j)/distance);%gamma角，即射线与穿过源射线的夹角
theta=r+gamma;%cos(theta)*x+sin(theta)*y=t中的theta
t=(distance*ray_x(j))/sqrt(distance*distance+ray_x(j)*ray_x(j));%cos(theta)*x+sin(theta)*y=t中的t
bias_x=t*cos(theta);%扇形束采样点的x偏置
bias_y=t*sin(theta);%扇形束采样点的y偏置
CX_real=CX_1.*cos(-gamma)+CY_1.*sin(-gamma)+bias_x;
CY_real=CY_1.*cos(-gamma)-CX_1.*sin(-gamma)+bias_y;%采样点旋转gamma角并加上偏置得到射线采样点
inter=interp2(X, Y, img, CX_real, CY_real, 'linear',0)';%线性插值得到各采样点的值,如果可以最好不要使用interp函数，很慢
proj_1(j)=sum(inter,1);
end
proj(:,i+1)=proj_1;
end
end

function recon=FanBeamFBP(sizeM,sizeN,beta,R,rayNum,dS,distance)%反投影过程
M=sizeM;%重建图像的行数
N=sizeN;%重建图像的列数
ray_x=-(rayNum-1)/2 * dS:dS:(rayNum-1)/2 * dS;
ray_x=ray_x';
weight_d=ones(rayNum,1).*distance;
weight=weight_d./(sqrt(weight_d.*weight_d+ray_x.*ray_x));%得到权重
clear weight_d;
l=length(ray_x);
clear ray_x;
for i = 1:length(beta)
R(:,i)=R(:,i).*weight;%投影加权
end
clear weight;
t=linspace(-N/2,N/2-1,N);
filt=0.0085*(sinc(t)/2-sinc(t/2).^2/4);%得到滤波的卷积核
p=R;
R = zeros(l+2*N,360);
R(N+1:l+N,:) = p;
for n = 1:360
temp(:,n) = conv(R(:,n),filt);
end
R=temp(3*M/2:3*M/2+l,:);%投影滤波
clear p;
recon=zeros([M,N]);%重建图像
M2=(M-1)/2;
dbeta=deg2rad(1);
for i=beta%射线源旋转每一个角度
r=deg2rad(i);%角度转弧度
for x=1:M
for y=1:N%对于重建图像中的每一个点
sp=distance*((y-M2-1)*cos(r)+(x-M2-1)*sin(r))/(distance+(y-M2-1)*sin(r)-(x-M2-1)*cos(r));
%得到s'，这里的推导可以试着自己推一推，过(x,y)和射线源的点和D‘2(该线随β变化而变化)的交点与O连接线段的长度。再注意：这段循环中的x和y并不是正常的x和y，x表示行数而y表示列数，与我们习惯中的x和y的意义相反，因此推得公式的x，y要互换
if (sp>=-(l-1)/2)&&(sp<=(l-1)/2)%若s’没有超出探测器的范围
isp = ceil(sp); sp0 = isp-sp;
proj_sp=(1-sp0)*R(isp+(l-1)/2+1,i+1)+sp0*R(isp+(l-1)/2,i+1);%根据s‘的值用他前后的射线投影为他插值，这里没用interp1，能够很大程度上加快运行速度
recon(x,y)=recon(x,y)+dbeta*(1/2*power(U(distance,(y-M2-1),(x-M2-1),r),2))*proj_sp;%重建图像上的像素点
end
end
end
end
end