CT图像重建.docx
《CT图像重建.docx》由会员分享，可在线阅读，更多相关《CT图像重建.docx（13页珍藏版）》请在冰豆网上搜索。
CT图像重建
昆明理工大学信息工程与自动化学院学生实验报告
（2009—2010学年第一学期）
课程名称：
医学成像系统与放射治疗装置开课实验室：
32082008年12月24日
年级、专业、班
学号
姓名
成绩
实验项目名称
CT图像重建
指导教师
刘利军
教
师评语
教师签名：
年月日
、实验目的与意义
医学成像技术是生物医学工程专业的一门重要的专业课程，课程主要涉及X光仪器，CT仪器，MRI
仪器和核医学仪器的工作原理及成像方法。
其中CT算法的出现又为后来数字化医学成像技术的发展提供
了基础。
该门课程为生物医学工程专业的专业基础课。
CT技术是医学成像系统中的一种重要手段。
它通过特定的算法，利用计算机的高速运算功能，可以在短时间内快速呈现人体断层图像。
让学生练习CT图像的重建有助于学生理解CT算法的内容，熟悉数
字图像重建的过程。
同时也能培养学生的团队精神和解决实际问题的能力。
二、实验算法原理
1、MATLA处理数字图像的基本函数；
2、X-CT三维图像重建的基本算法。
CT图象重建有四种基本的算法：
矩阵法，迭代法，傅立叶算法，反投影算法.我们采用的方法为卷积反投影•卷积反投影有：
平行光束投影的卷积反投影算法，等角扇形光来投影的重建算法•
1）.平行光束投影的卷积反投影算法
从投影重建三维物体的图像，就是重建一个个横断面。
这样三堆图像的重建就归结为二维图象的重建。
二维图像的重建问题可以从数学上描述如下。
假定g（x,y）表示一个二维的未知函数，通过g（x,y）的直线称为光钱（见图2.1）。
沿光线g（x,y）的
积分称作光线积分。
沿相同方向的一组光线积分，就构成一个投影。
图2.1中垂直于直线CC'（与X轴
夹角为）的光线所形成。
图2.1g（x,y）在方向的投影P（t）
的投影P（t），称之为g（x,y）在方向的投影。
光线积分和投影在数学上可以定义如下：
在图2.1中直线AB的方程为：
XcosYsinti（2.1）
其中ti是AE到原点的距离，g（x,y）沿AB的积分为：
P（ti）ABg（x,y）dsg（x,y）（xcosysintjdxdy（22）
对于给定的，g（x,y）在方向的投影P（t）是t的函数。
如果g（x,y）在各个方向的投影已知，g（x,y）
就可以唯一确定。
下面就讨论卷积反投影重建算法。
假定投影方向，如图2.2，将坐标（x,y）旋转角（逆时针方向）形成坐标系（t,s）。
g（x,y）在（t,s）
坐标系中为g（t,s）。
图2.2傅立叶切片定理示意图
坐标系（t,s）与（x,y）之间的关系为：
tcossinx
ssincosy
显然
Ptg（t,s）ds
令S（w）为P（t）的傅立叶变换则
（2.3）
（2.4）
S（w）P（t）exp（j2wt）dt
g（t,s）exp（j2wt）dsdt
将上式变换到（x,y）坐标系中，注意到变换的可比行列式
（2.5）
II
cossin
sincos
（2.6）
从而得到：
S（w）
g（x,y）exo[j2（xcosysin）]|dxdy
g（x,y）exp[j2（uxvy）]dxdy
（2.7）
其中
cossin
（2.8）
若令g（x,y）的傅立叶变换为G（u,v），由（2.8）可知
S（）G（u,v）G（,）（2.9）
若g（x,y）的傅立叶变换为G（u,v）的极坐标表示。
这说明g（x,y）在方向的投影P（t）
傅立叶变换S（w）等于G（u,v）在与u轴成角的直线上的值。
这就是著名的傅立叶投影切片定理。
可见
在整个（u,v）平面G（u,v）可以利用各个方向的投影来得到，从而g（x,y）也可以通过求G（u,v）的傅立叶
反交换的办法求得:
II
变换到极坐标中
cossin
得到
经推导得
其中
若令
其中
当图像的频谱是有限带宽时，则上式变为
（2.17）
0
h（t）J|e）p（j2t）d
由于图象及其频谱都是离散采样的，假定图象采样间隔为，则根据采样定理01/2。
为了进行
数学处理，只需知道h（t）在有限带宽上的离散采样点的值•这样我们有
1/（42）
h（n）0（2.18）
222
1/n222
其中n为正负整数。
（2.18）的离散形式为
Q（n）P（m）h（nm）（2.19）
m
假定P（m）在m0,1,N1之外的值为0,则上式变为
其中n0,1,2,N1从而可见为确定P（t）的N个采样点上的Q（n）的值，需要使用h（n）的2N—1个点上的值，从n=—（N—1）到（N—1）。
为求得Q（n），利用傅立叶变换计算卷积是比较快的方法，为清除循环卷积的周期交叠效应，实际
上h（n）取2N个点,P（m）补0，使之有（2N—1）个元素，则P（t）在N个采样点上就避免了交叠，如果使用以2为基的FFT（快速傅立叶变换）算法，P（m）和h（n）都必需朴0至（2N一1）个元素，（2N—1）为大于等于2N-l的最小的2的整数幕。
计算Q（n）的过程可以写为
Q（n）[FFT（FFT[P（n）0]FFT（h（n）0]（2.22）
其中FFT和IFFT分别表示快速傅立叶变换和反变换，光滑窗是在滤波过程中加入的光滑因子，例如引用
汉明窗，有时可以改进重建效果。
对于各个方向的投影，得到Q（n）之后就可以由（2.22）来计算g（x,y）。
重建步骤可以归纳为:
第二步。
反投影，由（2.14）的近似形式
M
g（x,y）Qi（xcosiysinJ（2.23）
Mii
来计算g（x,y）的近似值0（x,y）。
M为投影个数j为投影方向角，他们均匀的分布在0~的范围内。
当计算Qi（xcosiysinJ时，txcosiysini，不一定在Q（n）的整离散点上，这就需要插
值求得，预先将Q（n）插值加密，即最靠近的点，可以提高计算速度。
2）.等角扇形光来投影的重建算法
几乎所有的快遗CT设备都是用的扇形光束。
这里叙述的是等角度光束投影，如图2.3，测量投影数据
的探测器等间距地分布在D1D2弧上，弧的半径为2D,D为光源到图像中心的距离。
在下文中，f（r,）
图象在极坐标中的表示。
R（）表示在方向角为的投影中位量角为的光线产生的投影数据。
通过中
心的光线其为0。
L表示从光源到像素（r,）的距离。
图2.3等扇形束投影重建算法中的变量
L（,,）,D2r22Drsin（）（2.24）
表示在方向角为的投影中通过像素（r,）的光线的位置角
1PE丄1cos（）
（,,）tantan（2.25）
SEDsin（）
图像f（r,）和扇形投影R（）有下述关系
其中r和是投影中光线的最大位置角，从而可以得到这种重建算法的执行步骤:
第一步：
投影的修改，假定投影的抽样间隔为，抽样数据R（n）通过下式进行修改，
R'（n）R（n[D^osn（2.27）
第二步：
卷积（滤波）将第一步修改了的投影与响应函数g（）进行卷积
'112
g（）-（）h（）（2.28）
2sin
Q（）R（n）g（nm）d（2.29）
其离散形式为：
M
Q（n）R（m）g（nm）（2.30）
m1
这里也和上节一样可以加进一光滑窗，改进重建效应。
第三步：
反投影，就是执行（2.30）的积分
f（r,）
21
（2.31）
2Q（）d
0L2（,,）'）
近似的有：
f（x,y）
2M1
（）
（2.32）
2Q（）
mi1L（,,）
其f（x,y）是图象f（x,y）的近似图像，xrcos，yrsin,M是投影数，假定投影均匀地分布在0~2内。
为求得滤波后的投影Q，对象素（x,y）的贡献，首先由（2.24）和（2.25）计算出L和’•所有Q对（x,y）的贡献加起来再乘以2/M就得f（x,y）。
四、实验程序代码
原图象代码：
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:
10x0=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（AA2*cos（theta-alpha）A2+BA2*sin（theta-alpha）A2-（R-xO*cos（theta）-yO*sin（theta））^2）;
NN=AA2*cos（theta-alpha）A2+BA2*sin（theta-alpha）A2;g（i,ax）=rho*2*A*B*MM/NN;
forj=1:
N/2
R=R+vstep;
forw=forw+1;
MM=sqrt（AA2*cos（theta-alpha）A2+BA2*sin（theta-alpha）A2-（R-x0*cos（theta）-y0*sin（theta））A2）;
NN=AA2*cos（theta-alpha）A2+BA2*sin（theta-alpha）A2;g（i,forw）=rho*2*A*B*MM/NN;
R=-R;
back=back-1;
MM=sqrt（AA2*cos（theta-alpha）A2+BA2*sin（theta-alpha）A2-（R-x0*cos（theta）-y0*sin（theta））A2）;
NN=AA2*cos（theta-alpha）A2+BA2*sin（theta-alpha）A2;g（i,back）=rho*2*A*B*MM/NN;
R=-R;
end
end
radon0（k,:
）=real（sum（g））;
end
%%%%%%%%%%%%generateR-Lfunction%%%%%%%%%%%%%%%%
radon1=[zeros（k,N）,radon0];
ax=N+1;
RL（ax）=1/（4*vstepA2）;
forw=ax+1;
back=ax-1;
fork=1:
N/2
n=2*k-1;
RL（forw）=-1/（n*pi*vstep）A2;
RL（back）=RL（forw）;
RL（forw+1）=0;
RL（back-1）=0;
forw=forw+2;
back=back-2;
end
fork=1:
（180/angle+1）
radon2（k,:
）=conv（radon1（k,:
）,RL）;
end
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/（p^2*vstepA2*（4*门人2-1））;
end
fork=1:
（180/angle+1）
radon2（k,:
）=conv（radon1（k,:
）,SL）;
end
radonf=radon2（:
2*N:
3*N）;
figure
（1）
subplot（321）
plot（1:
（2*N+1）,radon1（1,:
））
title（'投影函数（已补零）'）
subplot（323）
plot（1:
（2*N+1）,SL）
title（'S-L卷积函数'）
subplot（325）
plot（1:
（4*N+1）,radon2（1,:
））
title（'卷积结果'）Xradon1=fft（radon1（1,:
））;
subplot（322）
plot（1:
（2*N+1）,abs（Xradon1））
title（'频谱'）
XRL=fft（SL）;
subplot（324）
plot（1:
（2*N+1）,abs（XRL））
Xradon2=fft（radon2（1,:
））;
subplot（326）
plot（1:
（4*N+1）,abs（Xradon2））
%%%%%%%%%%%%%%iradon%%%%%%%%%%%%%%%%%%%%%%%%%%fork=1:
（180/angle+1）
theta=（k-1）*angle*pi/180;
C=N/2-（N-1）*（cos（theta）+sin（theta））/2;
fori=1:
N
R=（i-1）*cos（theta）+C;nO=floor（R）;
ifn0>0&&n0<（N+1）
dot=R-nO;
l（i,1,k）=（1-dot）*radonf（k,nO）+dot*radonf（k,nO+1）;
else
I（i,1,k）=nan;
end
forj=2:
N
R=R+sin（theta）;
n0=floor（R）;
ifn0>0&&n0<（N+1）
dot=R-n0;
I（i,j,k）=（1-dot）*radonf（k,n0）+dot*radonf（k,n0+1）;else
I（i,j,k）=nan;
end
end
end
end
Ifinal=sum（I,3）;
Gfinal=mat2gray（lfinal）;
Gfinal=imrotate（Gfinal,90）;
figure
（2）
imshow（Gfinal）
五、实验结果与分析
反投影一般步骤为：
原像重建后图像
取投影反投影重建
程序流程:
卷积反投影法结果:
重建图象原图象
六、心得体会
通过本次的实验，对CT图象重建的基本方法之一：
卷积反投影，有了进一步的认识，在实验的过程中，采用的图象是经典的Shep-Logen头模型，得到的结果与原图象相比，有一定的差异，但影响不大•有待进一步的改进算法•
七、参考文献
[1]《医用电子仪器及装备一一医学成像系统及放射治疗装置类》，唐庆玉主编，清华大学电机系生物医
学工程及仪器专业
[2]《医学成像系统》，高上凯主编，清华大学出版社
[3]G.T.赫尔曼著（严洪编译），由投影重建图象，科学出版社
[4]董雏申、吴世法、王天童，卷积反投影法从x光图像重建三维轴对称图象，全国图基科学会议论文集，
1989。
⑸G,N.HousfleldComputedMedicalImaging-Journal0fComputedAssistedTomography,4（5）,
655—674．1980
[6]赶荣椿等，《数字图像处理导论》，百北工业大学出版牡，西安，1995，p77—179
[7]庄天戈，《CT原理与算法），上海交通大学出版杜，1992,p31—33
[8]RadonJ，UberdieBestimmumgvonFunktionenmuchihreintegralwertelangsgewissser
MannnmgfaltigkeitenBerichteSachsischeAkademieder
Wissenschaften,Leipzig,Math．Phys．K1,1917.69：
262~267
（注：
专业文档是经验性极强的领域,无法思考和涵盖全面,素材和资料部分来自网络,供参考。
可复制、编制,期待你的好评与关注）
升级会员 