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

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

我接触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

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

【转】一篇文章完整了解CT成像技术（完整版）
1．CT的发明与发展 1.1 CT的发明 CT是计算机断层摄影术(Computed Tomography,CT)的简称,是继1895年伦琴发现X线以来,医学影像学发展史上的一次革命. CT的发明可以追 ...
python框架下的CT成像重建后处理（CIL库）
CT重建 1 python重建库 1.1 scikit-image库(平行束) 1.2 ODL库(数据格式奇怪,失败.) 2 MATLAB重建库 2.1 iradon重建平行束 2.2 ifanbea ...
口腔ct重建服务器原理,口腔CT成像的迭代重建算法研究
摘要: 自CT技术诞生以来,它就被认为是二十世纪影响人类发展的十大技术之一.其应用更是涉及了多个领域,包括医学,生物,工业,安全等等.特别是在医学方面,CT更是成为了医学检测的主要手段.而口腔CT(D ...
(Python)卫星RPC有理多项式模型读取与正反投影坐标计算原理与实现
(Python)卫星RPC有理多项式模型读取与正反投影坐标计算原理与实现文章目录 (Python)卫星RPC有理多项式模型读取与正反投影坐标计算原理与实现摘要 RPC几何定位模型介绍 RPC模型库 ...
【医学图像处理】CT成像技术之CT剂量
[医学图像处理]CT成像技术之CT剂量剂量单位我们可以将不同的剂量测量视为逐步进行,每次将额外的变量添加到方程式中. CT剂量指数(CTDI) 首先,我们通过单机架旋转来测量检测器的剂量,以得出C ...
扇束重建的过程（弧形探测器）
本部分与<扇束重建的过程(平板探测器)>基本相同,只需要替换其中的某些变量即可. %% 清除内存 clc; clear; close all;%% 生成头模型以及设置头模型的参数 N=25 ...
计算机数学方法记录图像原理,【从业必读】CT成像基本原理
原标题:[从业必读]CT成像基本原理 CT成像基本原理是用X线束对人体检查部位一定厚度的层面进行扫描,由探测器接收透过该层面的X线,转变为可见光后,由光电转换器转变为电信号,再经模拟/数字转换器(an ...
【医学图像处理】CT成像技术之CT伪影
[医学图像处理]CT成像技术之CT伪影导致CT图像伪影的因素可以分类为: 物理因素病人属性扫描器因素螺旋和多层切片因素物理因素光束硬化 X射线束具有不同能量的光子,这些能量在平均" ...
【医学图像处理】CT成像技术之CT设备
[医学图像处理]CT成像技术之CT设备本系列重点介绍CT成像技术,并将介绍用于获取图像的设备,如何形成和显示图像,影响图像质量的因素以及如何测量剂量. 关注公众号,同样精彩内容本章节重点介绍CT成 ...

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

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

最新文章

热门文章