CT图像重建Word格式.docx
《CT图像重建Word格式.docx》由会员分享，可在线阅读，更多相关《CT图像重建Word格式.docx（15页珍藏版）》请在冰豆网上搜索。
让学生练习CT图像的重建有助于学生理解CT算法的内容，熟悉数字图像重建的过程。
同时也能培养学生的团队精神和解决实际问题的能力。
二、实验算法原理
1、MATLAB处理数字图像的基本函数；
2、X-CT三维图像重建的基本算法。
CT图象重建有四种基本的算法:
矩阵法,迭代法,傅立叶算法,反投影算法.我们采用的方法为卷积反投影.卷积反投影有:
平行光束投影的卷积反投影算法,等角扇形光来投影的重建算法.
1）.平行光束投影的卷积反投影算法
从投影重建三维物体的图像，就是重建一个个横断面。
这样三堆图像的重建就归结为二维图象的重建。
二维图像的重建问题可以从数学上描述如下。
假定表示一个二维的未知函数，通过的直线称为光钱（见图2．1）。
沿光线的积分称作光线积分。
沿相同方向的一组光线积分，就构成一个投影。
图2．1中垂直于直线（与X轴夹角为）的光线所形成。
图2.1在方向的投影
的投影，称之为在方向的投影。
光线积分和投影在数学上可以定义如下：
在图2.1中直线AB的方程为：
（2.1）
其中是AB到原点的距离，沿AB的积分为：
（2.2）
对于给定的，在方向的投影是t的函数。
如果在各个方向的投影已知，就可以唯一确定。
下面就讨论卷积反投影重建算法。
假定投影方向，如图2.2，将坐标旋转角（逆时针方向）形成坐标系。
在坐标系中为。
图2.2傅立叶切片定理示意图
坐标系与之间的关系为：
（2.3）
显然
（2.4）
令为的傅立叶变换则
（2.5）
将上式变换到坐标系中，注意到变换的可比行列式
（2.6）
从而得到：
（2.7）
其中
（2.8）
若令的傅立叶变换为，由（2.8）可知
（2.9）
若的傅立叶变换为的极坐标表示。
这说明在方向的投影
傅立叶变换等于在与u轴成角的直线上的值。
这就是著名的傅立叶投影切片定理。
可见在整个平面可以利用各个方向的投影来得到，从而也可以通过求的傅立叶反交换的办法求得：
（2.10）
变换到极坐标中
，
得到
（2.11）
经推导得
（2.12）
若令
（2.13）
则
（2.14）
（2.13）式右端是两频谱函数和的乘积的傅立叶反变换。
是投影
傅立叶变换。
若的傅立叶反变换为，则根据卷积定理有：
（2.15）
或
（2.16）
当图像的频谱是有限带宽时，则上式变为
（2.17）
由于图象及其频谱都是离散采样的，假定图象采样间隔为，则根据采样定理。
为了进行数学处理，只需知道h（t）在有限带宽上的离散采样点的值．这样我们有
（2.18）
其中n为正负整数。
（2.18）的离散形式为
（2.19）
假定在之外的值为0，则上式变为
（2.20）
（2.21）
其中从而可见为确定的N个采样点上的的值，需要使用的2N—1个点上的值，从n=一（N—1）到（N—1）。
为求得，利用傅立叶变换计算卷积是比较快的方法，为清除循环卷积的周期交叠效应，实际上取2N个点，补0，使之有（2N—1）个元素，则在N个采样点上就避免了交叠，如果使用以2为基的FFT（快速傅立叶变换）算法，和都必需朴0至（2N一1）个元素，（2N一1）为大于等于2N—l的最小的2的整数幂。
计算的过程可以写为
（2.22）
其中FFT和IFFT分别表示快速傅立叶变换和反变换，光滑窗是在滤波过程中加入的光滑因子，例如引用汉明窗，有时可以改进重建效果。
对于各个方向的投影，得到之后就可以由（2.22）来计算。
重建步骤可以归纳为：
第一步：
卷积，也称滤波，由（2．22）对每个方向计算。
第二步。
反投影，由（2．14）的近似形式
（2.23）
来计算的近似值。
M为投影个数为投影方向角，他们均匀的分布在0~的范围内。
当计算时，，不一定在的整离散点上，这就需要插值求得，预先将插值加密，即最靠近的点，可以提高计算速度。
2）.等角扇形光来投影的重建算法
几乎所有的快遗CT设备都是用的扇形光束。
这里叙述的是等角度光束投影，如图2.3，测量投影数据的探测器等间距地分布在弧上，弧的半径为2D，D为光源到图像中心的距离。
在下文中，图象在极坐标中的表示。
表示在方向角为的投影中位量角为的光线产生的投影数据。
通过中心的光线其为0。
L表示从光源到像素的距离。
图2.3等扇形束投影重建算法中的变量
（2.24）
表示在方向角为的投影中通过像素的光线的位置角
（2.25）
图像和扇形投影有下述关系
（2.26）
其中和是投影中光线的最大位置角，从而可以得到这种重建算法的执行步骤：
投影的修改，假定投影的抽样间隔为，抽样数据通过下式进行修改，
（2.27）
第二步：
卷积（滤波）将第一步修改了的投影与响应函数进行卷积
（2.28）
（2.29）
其离散形式为：
（2.30）
这里也和上节一样可以加进一光滑窗，改进重建效应。
第三步：
反投影，就是执行（2.30）的积分
（2.31）
近似的有：
（2.32）
其是图象的近似图像，，，M是投影数，假定投影均匀地分布在0~内。
为求得滤波后的投影，对象素的贡献，首先由（2.24）和（2.25）计算出和．所有对的贡献加起来再乘以就得。
四、实验程序代码
原图象代码:
p=phantom（256）;
imshow（p）
卷积反投影法代码：
H=[000.920.6990*pi/1801;
...
0-0.01840.8740.662490*pi/180-0.98;
0.2200.310.1172*pi/180-0.2;
-0.2200.410.16108*pi/180-0.2;
00.350.250.2190*pi/1800.1;
00.10.0460.04600.2;
0-0.10.0460.04600.2;
-0.08-0.6050.0460.02300.1;
0-0.6050.0230.02300.1;
0.06-0.6050.0460.02390*pi/1800.1];
angle=1;
N=100;
vstep=2/N;
ax=N/2+1;
fork=1:
（180/angle+1）
theta=（k-1）*angle*pi/180;
fori=1:
10
x0=H（i,1）;
y0=H（i,2）;
A=H（i,3）;
B=H（i,4）;
alpha=H（i,5）;
rho=H（i,6）;
R=0;
forw=ax;
back=ax;
MM=sqrt（A^2*cos（theta-alpha）^2+B^2*sin（theta-alpha）^2-（R-x0*cos（theta）-y0*sin（theta））^2）;
NN=A^2*cos（theta-alpha）^2+B^2*sin（theta-alpha）^2;
g（i,ax）=rho*2*A*B*MM/NN;
forj=1:
N/2
R=R+vstep;
forw=forw+1;
g（i,forw）=rho*2*A*B*MM/NN;
R=-R;
back=back-1;
g（i,back）=rho*2*A*B*MM/NN;
end
radon0（k,:
）=real（sum（g））;
end
%%%%%%%%%%%%generateR-Lfunction%%%%%%%%%%%%%%%%
radon1=[zeros（k,N）,radon0];
ax=N+1;
RL（ax）=1/（4*vstep^2）;
forw=ax+1;
back=ax-1;
n=2*k-1;
RL（forw）=-1/（n*pi*vstep）^2;
RL（back）=RL（forw）;
RL（forw+1）=0;
RL（back-1）=0;
forw=forw+2;
back=back-2;
radon2（k,:
）=conv（radon1（k,:
）,RL）;
radonf=radon2（:
2*N:
3*N）;
%%%%%%%%%%%%generateS-Lfunction%%%%%%%%%%%%%%%%%%
radon1=[zeros（k,N）,radon0];
forv=1:
（2*N+1）
n=v-N-1;
SL（v）=-2/（pi^2*vstep^2*（4*n^2-1））;
fork=1:
）,SL）;
figure
（1）
subplot（321）
plot（1:
（2*N+1）,radon1（1,:
））
title（'
投影函数（已补零）'
）
subplot（323）
（2*N+1）,SL）
S-L卷积函数'
subplot（325）
（4*N+1）,radon2（1,:
卷积结果'
Xradon1=fft（radon1（1,:
））;
subplot（322）
（2*N+1）,abs（Xradon1））
频谱'
XRL=fft（SL）;
subplot（324）
（2*N+1）,abs（XRL））
Xradon2=fft（radon2（1,:
subplot（3
升级会员 