基于MATLAB的扇形束投影CT重建毕业作品.docx
《基于MATLAB的扇形束投影CT重建毕业作品.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的扇形束投影CT重建毕业作品.docx（45页珍藏版）》请在冰豆网上搜索。
基于MATLAB的扇形束投影CT重建毕业作品
任务书
I、题目：
基于matlab的扇形束投影CT重建
II、使用的原始资料（数据）及设计技术要求：
1、原始资料：
matlab2013软件1套；计算机
2、毕设要求：
（1）查找扇形束投影CT重建资料；
（2）掌握扇形束投影计算方法；
（3）掌握扇形束投影CT重建方法；
（4）数据处理。
III、工作内容及完成时间：
1、查阅扇形束投影CT重建的相关资料并撰写开题报告；03.09－04.10
2、研究仿真建模与扇形束投影计算方法；04.11－04.15
3、研究扇形束投影CT重建算法；04.16－04.30
4、研究等角型和等距型投影重建；05.04－05.15
5、对比分析；05.16－05.31
6、总结并撰写论文，答辩。
06.01－06.26
Ⅳ、参考文献：
[1]张朝宗.工业CT技术和原理[M].北京：
科学出版社,2009.
[2]庄天戈.CT原理与算法[M].上海:
上海交通大学出版社,1992.
[3]张聪哲,李晓苇,杨昆.扇束等角型CT与等距型CT的比较研究[J].CT理论与
应用研究,2013,22
（2）:
215-223.
[4]GaborT.Herman.Fundamentalsofcomputerizedtomography:
image
Reconstructionfromprojections[M].springer,2009.
[5]张全红,路宏年,杨民,傅健.用对称反投影及递归迭代实现扇束CT快速重建
[J].CT理论与应用研究,2004,13（4）:
16-19.
基于MATLAB的扇形束投影CT重建
摘要：
伴随着CT成像技术的逐步发展，且因为扫描速度慢，成像质量差等缺点，第一代成像系统中的平行束CT成像技术，已逐渐被图像重建易于实现和控制的扇形束CT成像技术取代，且扇形束CT成像技术又被分为等角型与等距型。
本论文主要通过MATLAB软件，先将扇形束投影数据重排为平行束投影数据，而后利用滤波反投影（FILTEREDBACKPROJECTION，FBP）重建算法，在改变检测对象、旋转增角、探测器间距、插值方式和滤波方式等条件的情况下重建图像。
之后再通过对比分析重建图像及其对应的峰值信噪比，来研究这些条件对等角型、等距型这两种检测方法下对成像质量的影响。
本论文对实际应用的参数选择具有积极的意义，可在此基础上得出在在实际应用中的建议参数。
在本研究中，采用控制变量法进行分析，得到以下的结论为：
重建图像的精度与以上参数均有较大联系，且当试验参数有变化时该定性关系也会产生较大变化；且总体来说，等距型的重建质量高于等角型的重建质量。
关键词：
扇形束图像重建CT等角型等距型
指导老师签名：
TheCTProjectionofFanbeamConstructionBasedonMatlab
Abstract:
AlongwiththegradualdevelopmentofCT,thefirstgenerationimagingsystem（parallelbeamCT）hasbeentakenplacedbythesecondgenerationimagingsystem（fanbeamCT）gradually,becausetheformerhasslowscanningspeedandpoorimagingqualityandthelatterismorelikelytobereachedandcontrolled.AndfanbeamCTimagingtechnologyincludeequiangularfanbeamCTandequidistantfanbeamCT,accordingtoitsFanSensorGeometry.
ThispapermainlyresearchimagereconstructionoffanbeamwithMATLAB.Firstly,itwillrearrangethefanbeamprojectiondatatogetprojectiondataofparallelbeam.ThenitwillreconstructtheimagewithFBPalgorithmsandresearchitbychangingitsparametersofobjectdetection,FanRotationIncrement,FanSensorSpacing,InterpolationandFilter.Thenitwillcontrastthereconstructionimagesandtheirpeaksignal-to-noiseratio（PSNR）tostudytheeffectsoftheseconditionsontheimagingqualityofequiangularfanbeamCTandequidistantfanbeamCT.Thisthesis’spositivesignificanceforpracticalapplication,istogetthesuggestionparametersinthisapplication.
Inthispaper,usingcontrolvariablemethodtocarryontheanalysis,wecangetthefollowingconclusions:
imagereconstructionaccuracyhasgreatercontactwiththeaboveparameters.Whentestparameterschange,thequalitativerelationswillhaveagreaterchange.Ingeneral,thereconstructionqualityofequiangularfanbeamCTisbetterthantheequidistantfanbeamCT’s.
Keywords:
fanbeamimagereconstructionCTequiangularequidistant
SignatureofSupervisor:
基于MATLAB的扇形束投影CT重建
1.引言
1.1课题的研究背景
图像可以被用来描述物理系统或物体内部某些特性的分布。
通常，图像由光线投射或反射过光学仪器后产生。
然而在实际产生图像时，有时则需要用不可见的辐射探测物间接测量而形成图像。
当X射线照射人体或工件后，它蕴含的能量会由于吸收与散射的作用而衰减掉一部分，射线被衰减的能量取决于它照射的物质的原子系数（或原子组成）、密度及X射线的能量频谱。
在医学检测时，医用X射线透射过人体各种密度不同的组织器官时，由于这种差异的存在，它被吸收的能量也不尽相同。
因而X射线探测器上接收到的发自不同角度、不同位置照射到的X射线的量是不相同的。
在普通的医学X射线照片中，我们往往可以从中发现骨质，就是因为骨质的物质属性与其它组织有差别，它的物质密度大，它对射线的衰减较为严重（吸收量大），因而在透射照片上感光较差而呈现浅色的图像。
医学上的CT通常利用多个方向的X射线投影值，从物体外部检测到的数据来重建物体横截面信息，进而获得人体内部组织的密度分布。
利用数学推导，CT系统就能够重新建立起探测器截面的断层截面的图像。
最后，CT系统得到的CT图像是用不同的灰度值来反映各部分（人体器官、工件的结构）对X射线的吸收能力，通常用黑影表示低密度的吸收区，而白色部分表示高密度的吸收区。
这样就能把不同的结构区分开来这是一种获取人体内部信息的极其有效手段，人类洞察物体内部结构的能力得到了极大的增强。
数学家Radon、物理学家A.M.Cormack、工程师GN.Hounsfield先后对CT理论与设备进行了大量研究，最终取得了实质性突破。
而后这项技术广泛应用于诊断医学与无损检测技术方面，它为诊断人体疑难疾病、工件内部缺陷提供了一种无损害的优秀方法。
1.2国内外研究状况
CT早期的理论研究可追溯至上个世纪初。
1917年，奥地利数学家J.Radon[1]率先开展了通过投影重建原始断层图像的研究，并推出了Radon变换方程与Radon反变换方程，由此奠定了用投影图像重建原始断层图像的理论基础。
1956年，物理学家A.M.Cormak与Bracewell进一步建立了投影图像的精确重建理论[2]，并将这一重建理论重建出太阳微波的发射图像[3]。
1985年，GT.Herman在他的专著《Imagereconstructionfromprojections》中更系统阐述了CT的理论基础[4]。
1967年至1970年间，英国EMI公司的工程师G.Housfield成功研制出第一台用于临床用的计算机断层成像扫描装置，并于1971年将其正式安装在伦敦的AikinsonMorley医院[5]。
此后，CT技术在短期获得极大发展。
CT成像早期的平行束的投影图像重建中，旋转的耗时较长、扫描的速度过慢、效率过低、重建出来的图像伪影较为严重。
在此基础上扇形束CT就展现出其独特的优势。
之后逐渐产生了第二代、第三代，而且还包括其后的第四、五代CT扫描装置的产生。
其中各代CT装置都有其明显的特点，如：
仅第一代使用的是非扇形束的扫描方式；第二代则是效率较低的旋转或平移方式；第三、四代CT则都是使用的连续旋转的运动方式来进行扫描；第五代CT则较前者的优势是可以实现快速CT重建。
以下为各代CT装置的原理示意图[4]：
图11一代CT原理示意图
图12二代CT原理示意图
图13三代CT原理示意图
图14四代CT原理示意图
图15五代CT原理示意图
扇形束投影重建算法大致分为两类：
一类是重排算法，即将视图中采集到的扇形束数据通过运算改变成平行束的扫描数据，进而通过平行束常用的滤波反投影方法（卷积反投影重建算法与Radon反变换两种）进行图像重建；另一种则是采用扇形束投影的直接重建算法，这种方法只是在扇形数据进行加权的基础上运用扇形束特有的重建算法进行反投影重建。
本论文中所采用的方法为第一种中的Radon反变换算法。
2.CT成像相关技术理论
2.1CT成像的物理原理
Radon变换这个起源于积分几何的变换公式可以建立物体投影的数据和物体的实际体素之间的联系。
这一变换得出了求解这些关系的正确的数学变换，最终使重建图像成为可能。
这部分内容包含的基本数学原理有Radon变换及其Radon反变换等。
在工程运应用时，物体内部断层图像可以由Radon空间的投影数据进行Radon反变换后可以转化得到。
在实际情况中，因为投影角度是不连续的，故无法对投影数据进行直接的Radon反变换运算。
由此，人们创造出很多种让离散数据进行Radon变换的方法来方便计算机执行重建运算。
。
对于进行透照的单色X射线，物质对它的衰减遵循Lambert吸收定律[1],[5],[6]，这个定律可以用如下的公式来表示：
（2.1）
对上式两边进行积分运算，可以知道另一种X射线衰减的形式，如下：
（2.2）
若X射线透复合材料，则上面等式应改写为如下更常用的形式：
（2.3）
在CT中，通过测试穿过物体的大量的辐射路径求传输比I0/I。
重建算法的输入的线积分为对传输比取对数和相反数求得：
（2.4）
投影数据p位于上面等式的右端。
图像重建就是从投影数据来求得物体内部断层上各点的物体衰减系数的过程。
2.2CT成像的数学原理基础----Radon变换
2.2.1投影过程和Radon变换
f（x,y）表示物体横断面的体素密度分布函数，物体的旋转投影的数据是f（x,y）通过Radon变换就得出的，是f沿射线所在直线的积分值，下面这个表达式[7]就是Radon变换的表达形式。
已知f（x,y）为广义函数，它的线积分表达式在xoy平面上的表示为：
（2.5）
上式中所有的s和ϕ均是已知的量，其中
（2.6）
则称pf（s,ϕ）为函数f（x,y）的Radon变换，记作Rf=pf（s,ϕ）
这里L是在xoy面上的直线，s是原点到直线L的距离，ϕ表示直线L的垂直线与x轴正方向所成的角度。
s和ϕ是直线L的位置参数，直线L由所给定的s和ϕ唯一确定（如图2-1）。
为方便运算，应通过绕原点逆时针旋转，使坐标系xoy旋转ϕ。
新的坐标系如图2-2所示，用以下式子就可以表示直线L的参数方程：
（2.7）
联系式（2.7），由对弧长的曲线积分的计算式可知：
f（x,y）的Radon变换可以表示如下。
（2.8）
图21f（x,y）沿L的积分
图22坐标系xoy和sot的关系
式（2.8）中pf（s,ϕ）表示断层的投影数据，ϕ为投影角度，s指探元所处位置，在Radon变换的过程中，物体内部投影的一维数据由二维的断层结构通过投影转化得到，这等价于把一个二维函数变换至Radon空间内。
2.2.2Radon反变换
CT是根据照射射线后的射线强度数据推出物体横断面上各体素射线衰减系数的分布。
Radon反变换就是通过一定部位内众多的投影路径的函数积分值来求解衰减系数的分布。
若不存在特殊的条件限制被积函数，在实际工程中不可以确定它，因为必须有无穷多个积分值，但这一点不能实现。
J.Radon在1917年创建Abel积分方程并解之，于是获得了求逆公式。
这一原理巧妙运用了平均值思想，被称之为Radon定理[8]。
Radon定理：
定义f（x,y）广义函数，
是绝对可积的，并且f是在以任意点（x,y）为中心,r为半径的圆周上积分的平均值，可得：
（2.9）
（2.10）
因此只是由其Radon变换来唯一决定函数f（x,y）的值，而且
（2.11）
其中
f（x,y;q）表示以（x,y）为中心，以q为半径的圆的切线p=xcosϕ+ysinϕ+q上函数f（x,y）的积分平均值，即：
（2.12）
通过以上Radon定理可以求出Radon反变换的解析式如下：
（2.13）
上式中
指代对
的偏导数,
是极坐标形式下的f（x,y）。
所以，在实际应用时，投影数据在进行Radon反变换后就变换出物体内部断层体素密度分布图像。
但因为投影角度是间断的，因而不能对投影数据直接进行Radon反变换。
为了解决这一矛盾，人们发明了很多种数据处理的方法将离散投影数据进行Radon变换，使计算机完成重建运算成为可能，其中主要是代数迭代、滤波反投影算法。
2.3扇形束投影的滤波反投影重建算法
2.3.1扇形束CT的分类
扇形束射线的形式有两种：
一类是等角型，另一类是等距型。
等角型是指在固定的X射线源的位置下，扇形束射线投影的数据是在等间角的位置获得的，这种布置方式下的探测器是在一段圆弧上均匀分布的。
等距型则是指在固定的X射线源的位置下，扇形束射线投影的数据是在等距离的位置获得的，这种布置方式下的探测器则是等距地分布在一条与源-旋转中心连线的直线上，而且在这种情况下的探测器的间角是不同的。
详细结构图请见图23、图24
表21平行束与扇形束反投影的程序描述
2.3.2等角型扇束的重建公式
为了不重排投影与损失空间分辨率，则需推导出合适的重建公式，使之能够直接采用扇形束的投影采样。
于是就推动了扇形束CT重建公式的发展[9]-[11]。
所有的扇形束都可以用两个参数γ和β唯一地确定表示，其中β为x轴y与中心射线所形成的夹角，γ是由任意射线与中心射线（经过旋转中心的射线）之间的夹角，如下图2-3所示。
β为表示在某次观测中被使用时，射线的投影角度。
γ为探测器角度，可以唯一确定扇形中任何一条射线的位置。
在关于平行投影重建的研究中，任何任意特定射线也可以且仅可以被两个参数（s和ϕ）确定（其中s是指射线到旋转中心的距离，ϕ表示投影的角度）。
倘若符合如下的要求，扇形束中的一次投影采样p（β,γ）就是平行束中的一次投影采样p（s,ϕ）的一部分：
ϕ=β+γ，
s=Dsinγ（2.14）
在上式中，D表示旋转中心距射线的距离。
与平行束产生的投影类似，扇形束的重建公式同样可以可通过等式（2.14）推导，下面是我们直接给出的重建公式。
图23等角射线束的参数关系
首先，对投影数据作预加权，其中使用的加权因子是cosγ；
然后对加权后的滤波，在这一步中使用的滤波器是
（2.15）
此处通过一个具有约束限制的ramp-filter[12]的滤波器来进行采样投影
（2.16）
可知这个滤波器中的采样值如下：
（2.17）
最后是进行加权反投影
（2.18）
在上式中，R表示射线源与旋转中心的距离。
（2.19）
（2.20）
2.3.3等距型扇束的重建公式
图2-4表示等距型探元扇形束投影产生的图示。
D1D2指探元所处的位置，S0B指代绝对位置被坐标原点确定的特定射线。
图24等距型示意图
在图2-4所标示的坐标系中，得到如下的重建公式：
（2.21）
在上式中，U表示为：
（2.22）
（2.23）
（2.24）
（2.25）
式（2.25）中S1表示为过待求点（r,θ）的一条射线。
由此可见，在使用等距型扇形束的FBP重建算法对断层图像重建时，需要先对投影数据进行加权处理，之后进行滤波，最后由重建点的具体位置计算出的投影地址（S1），作反投影重建。
2.3.4扇形束投影转化为平行束投影的原理
如图2-5所示,它是由射线源
发出的等角型扇束射线形成的投影。
物体
表示在图中小圆内部，射线源位于大圆上，D为它离圆心的距离，R为物体所处的圆的半径，
是X射线源，弧
为探测器所在弧线，
为中心射线。
扇形位置由该中心射线与y轴交角
所确定，同一扇形中的任何一条射线
由
确定。
由于
是该射线相对于
绕
的转角。
x-y为固定的坐标系，因此射线的绝对位置由
确定。
如果把这条射线看作是平行射线，自然也可被
所确定。
图25发散数据采集的集合结构
通过这种思想就：
将扇形束CT的投影数据转化为平行束投影数据，然后再通过平行束的重建算法得到断层重建的图像。
等距型的重排原理与之类似这里不做重复介绍。
2.4计算机模拟
2.4.1投影模型
本论文中模拟使用的模型为256×256大小的Shepp-Logan模型[13],[14]和自建模型。
Shepp-Logan脑图模型首先使用于医学CT，是专业人员共同认可的用来评价各类图像重建算法有效性大小的研究对象。
S-L脑图模型是由多个尺寸、朝向、位置、密度各不相同的椭圆构成，表征着一个大脑的断层图像。
为了检测结果的普适性以及进一步比较重建图像细节，自建了一个具有典型代表的的断层图像。
该图像同样由各种大小、方向、位置、密度（表示为灰度）各异的椭圆（包括圆）与矩形（包括正方形）组成，以便于更贴近于工业CT使用实际对象。
这两种模型如下图所示：
图26Shepp-Logan模型
图27自建模型
影响CT重建图像的质量的因素包含主要两类：
一是硬件的性能，即探测器的尺寸与透照布置的设计等；二是软件，对数据的修正以及重建算法的选择和优化。
本文在模拟仿真中对扇形束CT硬件参数进行适当的设置，使用的算法为滤波反投影重建算法（FBP）。
2.4.2投影采集几何参数
利用控制变量法进行仿真模拟的默认重建参数设置如下表2-2：
表22仿真模拟默认参数设置
参数名
默认值
探测器置方式
Arc
探测器间距（SS）
1
探测器放射线源距旋转中心的距离（D）
250
旋转增角（RI）
1
滤波方式
Shepp-Logan
插值方式
Linear
频率范围
0-1
输出尺寸
256*256
对重建图像质量的评价采用一个目前比较广泛接受且客观的评价标准---计算重建图像的峰值信噪比（PeakSignaltoNoiseRate,即PSNR，单位：
分贝）值。
PSNR是指在原图像与处理后的图像之间的均方误差比上（2*n-1）*2以后的进行对数运算处理后得到的值（n是每个采样值的比特数），其计算公式如下：
（2.26）
上式中，MSE为均方误差（MeanSquareError），即：
（2.27）
上式中，Framesize表示数据的总个数，In、Pn分别表示原断层图像的第n个点的灰度值与重建原断层图像的第n个点的灰度值。
由于重建得到的断层图像会一定程度上与原始断层图像存在差异，但为了得到经过处理后的重建图像的品质，经常通过这种量化的计算方式来作他的参考。
且该值越大，重建图像与原始图像越差异越小，反之亦然。
然而在实际使用中因为人的视觉感官的主观性的影响，该值的大小有时与人眼所观察到的视觉品质有一定差异，所以利用该值进行重建图像的好坏需适当结合实际重建图像进行分析。
具体程序的实现请参见附录。
3.参数改变时等角型与等距型的对比研究
在对比研究中，我们采用控制变量的研究方法，即在其他条件不变的情况下，研究某一参数对图像重建质量的影响。
此处使用的模型为Shepp-logan模型，使用的其他参数为默认参数（详见表22）。
3.1等角型下的对比研究
3.1.1旋转增角（RI）的改变
此处使用的模型为Shepp-logan模型，使用的其他参数为默认参数（详见表22）。
此处主要研究旋转增角变化对重建图像质量的影响。
图31重建图像随旋转增角的变化（SS=0.1）
图32重建图像随旋转增角的变化（SS=2）
图31、图32表示在上表中探测器间距为0.1和2时，重建图像随旋转增角的变化（从左至右，再从上至下）。
表31不同RI、SS时的PSNR值（单位：
dB，其他参数为默认值）
SS
RI
0.1
0.5
1
2
4
8
0.1
74.96739
72.75131
69.19522
65.71654
63.59127
63.08335
0.5
74.9523
72.75269
69.19527
65.71653
63.59128
63.08336
1
74.60759
72.73729
69.19518
65.71646
63.59094
63.08325
2
73.38194
72.44966
69.18553
65.71598
63.59045
63.08325
3
70.93889
70.94361
69.12055
65.71386
63.5926
63.0847
4
68.84553
69.17092
68.47797
65.70979
63.58935
63.08321
5
67.0
升级会员 