一、CT图像重建原理
在CT成像中，物体对X射线的吸收起主要作用，在均匀物体中，X射线的衰减服从指数规律。在X射线穿透人体器官或组织时，由于人体器官或组织是由多种物质成分和不同的密度构成的，所以各点对X射线的吸收系数是不同的。将沿着X射线束通过的物体分割成许多小单元体（体素），令每个体素的厚度均为L。设L足够小，使得每个体素均匀，每个体素的吸收系数为常值，如果X射线的入射强度I0、透射强度I和物体体素的厚度L均为已知，沿着X射线通过路径上的吸收系数之和μ1＋μ2＋……＋μn就可计算出来。
为了建立CT图像，必须先求出每个体素的吸收系数μ1、μ2、μ3……μn。为求出n个吸收系数，需要建立如上式那样n个或n个以上的独立方程。CT成像装置从不同方向上进行多次扫描，来获取足够的数据建立求解吸收系数的方程。
吸收系数是一个物理量，它是CT影像中每个像素所对应的物质对X射线线性平均衰减量大小的表示。再将图像面上各像素的CT值转换为灰度，就得到图像面上的灰度分布，就是CT影像。
CT重建过程可以采用直接反投影和卷积反投影来实现。卷积反投影重建图像时，先把由检测器上获得的原始数据与一个滤波函数进行了卷积运算，得到各方向卷积的投影函数；然后再把它们从各方向进行反投影，即按其原路径平均分配到每一矩阵元上，进行叠加后得到每一矩阵元的CT值；再经过适当处理后就可以得到被扫描物体的断层图像，卷积反投影可消除单纯的反投影产生的边缘失锐效应，补偿投影中的高频成分和降低投影中心密度，并保证重建图像边缘清晰和内部分布均匀。
二、利用MATLAB2016a进行CT图像重建
1、仿真头部模型数据的获取
利用MATLAB中P = phantom(def, n)函数可以方便地直接获取Shepp-Logan头模型，其中字符型参数def在X射线断层扫描中研究者主要采用'Shepp-Logan'以及默认情况为'Modified Shepp-Logan' (default)，标量参数n指明矩阵P的行列数目，默认情况下为256。这样我们很容易获取模型的数据从而避免了利用C++等其他程序语言进行数据获取方式的不便。C++获取SheepLogan函数的思路为遍历图像，检测图像的每一个像素点是否依次在某一椭圆中，若在椭圆外，则置为背景色，在其他任何椭圆中，则置为该椭圆对应的灰度值。
P =phantom('Modified Shepp-Logan',200);
2、利用反投影法进行图像重建
图像重建的算法采取的是滤波反投影法，先通过X射线束在体层的各个方向上对各体素进行扫描，测得一系列的投影值，然后，把各个方向的投影值沿原路径反方向投影回与原体素空间位置一样的体素上，得到该体素在各方向上反投影值的总和，通过计算机运算，求出各体素值，实现图像的重建。
利用radon变换函数进行各个角度所获取的Shepp-Logan头像投影，也就是将图像在某一个方向上做线性积分（或累计求和）。如果将图像看成是二维函数f(x, y)，那么其投影就是在特定方向的线性积分，比如f(x, y)在垂直方向的线性积分就是其在x轴上的投影。函数R = radon(I, theta)中参数I为上节获取的灰度图像，theta为指定的投影角度。
theta0=0:36:144;
[R0,xp] = radon(P,theta0);%5个投影角度;
theta1=0:18:162;
[R1,xp] = radon(P,theta1);%10个投影角度;
theta2 = 0:10:170;
[R2,xp] = radon(P,theta2);%18个投影角度;
theta3 = 0:5:175;
[R3,xp] = radon(P,theta3);%36个投影角度;
theta4 = 0:2:178;
[R4,xp] = radon(P,theta4);%90个投影角度;
I0=iradon(R0,36);
I1=iradon(R1,18);
I2=iradon(R2,10);
I3=iradon(R3,5);
I4=iradon(R4,2);
3、滤波反投影图像重建过程
利用MATLAB中的iradon函数对上步获得的投影数据进行滤波反投影重建，获得Shepp-Logan模型的重建图像。
figure;
subplot(3,2,1),imshow(P)
title('(a)Shepp-Logan head phantom image')
subplot(3,2,2),imshow(I0)
title('(b)Backprojection from 5 projections');
subplot(3,2,3),imshow(I1)
title('(c)Backprojection from 10 projections');
subplot(3,2,4),imshow(I2)
title('(d)Backprojection from 18 projections');
subplot(3,2,5),imshow(I3)
title('(e)Backprojection from 36 projections');
subplot(3,2,6),imshow(I4)
title('(f)Backprojection from 90 projections');
三、X-CT图像重建MATLAB仿真结果
