CT-计算机断层扫描成像实验

1、第二章 CT-计算机断层扫描成像实验（系列实验二）射线成像实验室七月 13, 2006目录0引言21CT成像实验原理21.1概述21.2投影定理31.3卷积反投影重建算法41.4一种实际算法51.4.1推导与描述51.4.2框图72实验方案82.1概述82.2实验环境92.2.1硬件环境92.2.2软件环境102.3实验步骤102.3.1概述102.3.2具体步骤扫描数据处理122.4FAQ & Tips122.4.1工作目录是啥？122.4.2如何确定样品的起始位置和水平扫描的长度？122.4.3为什么扫描完成后要保存数据？132.4.4为什么图像多出一条

2、横贯全图的线？133附录：CTSystem软件使用说明书133.1概述133.2界面介绍133.2.1新建扫描项目133.2.2转台位置调整143.2.3调整能谱敏感区域143.2.4扫描属性153.2.5扫描163.2.6投影变换窗口173.3投影变换的输出184参考文献210 引言自七十年代初第一台电子计算机断层扫描装置问世以来，成像技术发展异常迅速，设备不断更新。以医学成像为例，已实现了三大飞跃，即脏器清晰图像的获得，把生化病理研究推向分子结构的水平和直接提供有关成像组织的化学成分的信息，步入了断层显像的新时代。计算机断层扫描和图像重建技术，是在不破坏物体情况下，将物体每一个断层面上的结

3、构和组份的分布情况显示出来的一种实验方法，都是利用计算机图像重建的方法来得到物体内部的信息。人们对射线成像的最早认识是从x光机开始的。医用x光机成像技术的发展和应用已有近百年的历史，它是利用x射线的物理性能和生物效应，来对人体器官组织进行检查。由于普通x光机只能把人体内部形态投影在二维平面上，因此会引起成像器官和骨骼等的前后重叠，造成影像模糊。为了克服这一缺点，英国ENI公司的工程师豪恩斯菲尔德(G.N.Hounsfield)运用了美国物理学家科马克(Cormack)于1963年发表的图像重建数学模型，推出了第一台x 射线计算机断层图像重建技术（X-CT）装置，并1977年9月在英国Ackin

4、son Morleg医院投入运行。1979年该技术的发明者Hounsfield和Cormack为此获得了诺贝尔医学奖。X-CT 的出现是X射线成像技术的一个重大突破。经过多代的发展，X-CT已获得广泛的应用。在医学上，目前已可用来诊断脊柱和头部损伤，颅内肿病，脑中血凝块，及肌体软组织损伤，胃肠疾病，腰部和骨盆恶性病变等等。目前X-CT除了广泛应用于临床诊断、生命科学和材料科学以外，还在工业和交通等方面也有重要的应用，例如，在线实时无损检测工业CT 等。1 CT成像实验原理1.1 概述数学上可以证明，通过对物体进行多次投影就可得到该物体的几何形状。CT的基本思想是：让一束射线投射在物体上，通过物

5、体对射线的吸收（多次投影）便可获得物体内部的物质分布信息。当强度为的一个窄束射线穿过吸收系数为 的物体时，其强度满足指数衰减关系式中为 射线所穿过物质层厚度。在实际情况中，所研究的物体往往不是由单一成分组成的，当物体由若干个不同成分组成时，物体内部各处的也将可能不同。在这样的物质中， 束穿过整个物件后的强度为式中为处的吸收率。CT系统通过改变一组射线路径，记录下对应出射强度的变化来分析物体内部的分布。在实际操作中，总是假定物体中的吸收系数是一个连续函数，通过射线测量方法和图像处理技术，将数学物理方程通过计算机解出函数。在计算机屏幕上，可用颜色或灰度来表示的大小，从而被扫描的物体的切面图像即可显

6、示出来。实际的扫描装置通常是由排列成一定角度的多组探测器构成的，这样在每一个位置就可以获得多组数据，从而节省了测量时间，提高了工作效率。共有三种信息收集方式：透射式CT（TCT）、放射式CT(ECT)、反射式CT。我们主要考虑前两种CT的成像原理。1.2 投影定理O样品坐标系扫描角时的临时坐标系Figure 1 坐标转换及射线位置示意图。图中虚线代表射线经过的路线；黑点代表我们感兴趣的一个点。我们需要求那里的吸收系数。投影定理或中心切片定理是图像重建算法的基础。设在角度，位置上的射线吸收大小为，即其中表示延Figure 1所示的虚线积分。再设物体的一个二维平面内吸收率分布为。那么投影定理在非衍

7、射源情况下，其内容为：按照的一维傅立叶变换是的二维傅立叶变换的一个（过原点的）切片。即该定理的具体证明可以参见参考文献1，或者参考沈激老师的核技术的笔记。1.3 卷积反投影重建算法那么，作为一个实际的CT系统而言，首要算法问题是如何从实验上只能测得的，算出我们所需的。CT的算法很多，常见的有：反投影重建算法（累加法）、滤波（卷积）反投影重建算法、直接傅立叶变换重建算法、迭代重建算法等。本节将介绍其中较为流行的卷积反投影重建算法。为了清楚起见，我们重新声明一下推导中将使用的变量。首先，在样品的一个断面上（如Figure 1所示），我们现在可以测得任意投影值（理想情况下为二维函数），如式所示。我们

8、把它们用视角和视线（Figure 1中虚线）距离原点的距离唯一标记。我们需要求的是该断面上，任意一点（或者用极坐标标记）上的吸收系数。的二维傅立叶变换为。是空间频率点在极坐标系中的坐标。由投影定理可知：其中表示二维傅立叶反变换，表示一维傅立叶变换。但是以上的傅立叶变换并不适合计算机运算。于是我们需要对它们进行进一步的化简：因为所以把和代入得其中函数被定义为第二个积分。它也是在该截面上定义的一个函数，且可以由算出。我们需要的只是它关于角度的积分。我们继续化简。忽略掉固定的系数，是两个函数乘积的傅立叶反变换，即可以表示为这两个函数分别的反变换的卷积。定义那么以上表述为其中与卷积，也是一次对的滤波。

9、滤波函数是。总之，为了求解，我们需要把与卷积，得到二维函数。然后把在反投影、延方向上累加即可。1.4 一种实际算法1.4.1 推导与描述以上算法仅限于位置和角度连续无限变化的理想情况，实际中显然无法直接使用。另一方面，滤波函数的响应函数不收敛（在0点趋向无穷）。因此，我们在离散化以上算法的基础上还要把有所改变。实际扫描过程中，得到的数据是在角度和位置方向上都等距分布的（由式）。为的位置取样间隔（每个角度上取样次，）。计算机读入，并把它储存在二维矩阵中。也要进行离散化，有限化处理（一种近似）。其中一种处理方法是R-L滤波函数1，即即，只去原滤波函数在采样频率范围中较低（）的部分。根据式可以算出，

10、对应的离散化空间频率响应函数：请注意，这里的标记可以同时取正值和负值。可以取，其中表示取整。根据式，将与离散卷积得到。当然，也是离散分布的。剩下最后一步就是求式表示的积分：首先，取一个矩阵，表示最后的图像。对于矩阵中的每取一个点，换算成极坐标下的。之后对于每一个角，算出（即所谓反投影），再从上插值得到。最后把按照不同累加起来，可得近似正比于的数值（式）。对于成像而言，这个正比的结果就足够了。另外有一个细节值得注意：的坐标是相对于原点即转轴取的。而我们实验中用到的CT装置无法自动给出转轴的位置（即上面的），因此需要在矩阵中分析得到。CT软件中固化的算法是这样计算的：首先求每个角度下的重心的位置（

11、相对于第一个点）的坐标：之后对不同角度下的一组数，用拟合得到常数项系数的值就是可以证明：如果被测物总是处在CT的探测范围内（射线束足够密级，范围足够大），且被测物吸收度变化不剧烈，那么被测物的吸收度重心在方向的投影点将是的重心。因此，随着视角的变化，重心的投影将以转动中心为轴线，按照三角函数随变化。，即1.4.2 框图延方向变化的一维数组扫描得到出射强度，接合式式转轴相对位置式一维数组式根据式算出一维数组对于每个卷积，式产生一组待求吸收度的点的算出对应极坐标对于每个对于每个根据算出用一种插值方法求出。比如线性插值按照式，对于不同累加，得到对应于点的吸收率把所有按照一定规则转换为颜色，画在的直角

12、坐标上结束开始延和排列的二维矩阵（与可以变量互换。去掉范围以外的数）延和排列的二维矩阵Figure 2 1.4.1节算法的流程框图。其中，双线框表示二维数组；粗线框表示一维数组；普通方框表示一个标量或者一个步骤；虚线框表示循环执行。循环变量已在方框边标出。2 实验方案2.1 概述本次实验的目的是得到样品水平断面的吸收系数分布图。易见，实验的步骤是首先由硬件系统取得样品的一系列投影，而后按照一定的算法算出各点的吸收度后以图像方式显示。第二步则完全是软件问题。比如对于Figure 3(a)所示的样品，这套系统能给出如Figure 3(b)所示的截面吸收系数分布图。受源强、准直孔、探测器等多方面条件

13、制约，目前的空间分辨只有毫米量级，扫描时间却持续数小时。(a)(b)Figure 3 扫描结果图。(a)，样品的示意图：一个空心铝制圆柱内放有两个铝制小棒，扫描平面垂直于拄体；(b)扫描结果图。2.2 实验环境2.2.1 硬件环境2探测系统探测器三维运动载物平台用户软件步进电机驱动程序驱动接口卡传感器电机驱动器主放大器高压源多道驱动程序多道分析卡待测物三维运动载物平台系统源准直器Figure 4硬件系统框图。图中虚线箭头表示软件联系；实线箭头表示电路连接；粗空心箭头表示力学连接；粗实心箭头表示准直后的射线。为了制造所需的硬件系统，需要一套移动系统、一套探测系统和一个放射源。因为移动源和探测器相

14、对较为困难，因此本项目中采用了移动待测物体的方案。放射源采用241Am。射线准直后经过放在三维运动载物平台上的待测样品，射入碲锌镉半导体探测器。探测器输出的信号经过放大后进入微机多道分析器。放射源。电脑协调探测与移动样品。整个系统框图如Figure 4。样品所处的位置请见Figure 6。之所以使用高分辨率的能谱探测器和多道，而不是简单地使用计数器作为探测系统，是因为一方面，它能最大限度地去除本底和散射对计数的干扰 实验证明，正常工作下，噪声计数小于信号的0.1%。；另一方面，通过拟合全能峰来获得计数比较精确。但是，目前我们还没有探测器阵列来提高扫描效率，因此目前的扫描时间较长。2.2.2 软

15、件环境从Figure 4中也可以看出，软件部分分为两个驱动程序和用户软件两个部分。驱动程序部分主要负责将用户软件发出的命令翻译给硬件，并且将硬件给出的数据反馈给用户软件。用户软件负责控制扫描、数据缓存、生成图像和进行图像变换、插值放大等处理。驱动程序部分用Microsoft Visual C+.net 2003编写。而用户软件部分中，界面和核心算法分别在Microsoft Visual Basic 6和Matlab R14环境下实现。2.3 实验步骤2.3.1 概述实验上，我们首先通过等间隔地取和，记录透射强度，来推算（或的常数倍）。的步长大概是，需要扫描的范围；步长左右，扫描范围约几个厘米。

16、如Figure 5，扫描的顺序是先把样品转到一定角度，再依次从开始，逐一读取各个上的计数。之后再把样品转到下一个角度，重复上述步骤直到转过半圈。因为，射线束无法移动，以上步骤最终都用移动和转动样品动作等效完成的。探测器准直孔放射源准直孔水平扫描方向旋转方向样品初始位置一轮水平扫描结束时的位置Figure 5 样品与射线束间的几何关系示意图。点划线表示射线束；实线笑脸表示样品的初始位置；虚线笑脸表示完成一个角度的水平扫描后，样品的位置。2.3.2 具体步骤 扫描水平扫描正方向转轴放射源准直孔探测器准直孔可三个自由度运动的载物平台Figure 6 载物台实物图1) 启动电脑与电源、N

17、IM机箱电源、高压源和3个步进电机电源。缓慢加压（10V/s）到400V。2) 把样品放在载物台上，尽量使物体的几何中心与转轴（载物台中心，如Figure 6）重合。3) 启动CT软件CTSystem，进入主界面，见Figure 7。其中，扫描类型选择“截面扫描”。最好选择统一的工作目录，比如“E:CT实验work”，再取一个工作名（最好不要使用特殊字符）。最后，程序会自动建立工作目录：“E:CT实验work你输入的工作名”，并在那里记录所有扫描数据。4) 调整转台位置（Figure 8）。让样品水平方向上处在，探测器准直孔和放射源准直孔连线之右（从放射源的角度看，因为水平扫描是从右向左进行的

18、）；垂直方向上，该连线应该处于我们所需的样品断层上。5) 启动放射源。6) 调整能谱敏感区域（Figure 9）。用“W”、“S”、“A”、“D”键移动两个游标。游标所夹的区域就是敏感能谱区域。工作时，只有能量在这个区域内的光子才被当作有效计数。能窗最好包括主光电峰，和两个X射线逃逸峰。7) 设置扫描参数（Figure 10）。角度的间隔（步长）大概是；X方向间隔左右。8) 开始扫描。9) （等待，找点娱乐活动，但最好不要乱动扫描中的电脑。）10) 扫描结束后保存一份报告，用于下一节的数据分析。11) 缓慢降低高压（10V/s）至0V，关闭NIM机箱电源、3个步进电机电源和放射源。12) 投影

19、变换（Figure 12）。（这台工控机比较老，喝杯咖啡，稍等一会）结束后会自动生成一份结果图到工作目录。当然也可以手工保存。13) 如有必要可以再次修改扫描参数，重新投影变换。 数据处理本实验要求学生自主完成一套CT成像算法，并要求将该算法实际编写出来。最后，利用你编的程序，把上节步骤12中得到的数据（透射计数值组成的数组和其他必要参数），处理成CT图像，并与步骤12的结果相比较。你可以使用1.4.1节的算法和1.4.2节的流程。当然如果主动调研或开发新的算法更好（CT算法还有很多种。即使对于滤波反投影算法，滤波函数也有很多种取法1）。这个过程中你可能会用到：1) 参考文献1中

20、介绍了数种算法。它可以在科大超新数字图书馆下载到。2) 对于实际编程，C、FORTRAN等“裸”语言效率较高，但你可能会发现十分繁琐。不过，能找到适于它们的计算函数库。这可以简化你的工作。3) 如果你使用上述语言编程，可以把结果以数组形式输出，再用EXCEL、ORIGIN等软件作图。4) 一种十分高效的开发手段是在MATLAB，MATHMATICA，LABVIEW等环境实践你的算法。宽松的运算规则和极其丰富的库函数会大大简化你的工作。实验中使用的CTSystem就是调用MATLAB引擎实现成像算法的。5) 知之为知之，不知google知。另外以上软件都可以在科大天狼搜索找到高速下载。2.4 F

21、AQ & Tips2.4.1 工作目录是啥？(For Win9X Ed.)新建一个扫描前，软件会提示简历一个工作目录，这个目录将被用来：1) 记录所有数据点的能谱，文件名是“你输入的工作名的前四个字母坐标.dat”2) 保存其他文件的默认目录2.4.2 如何确定样品的起始位置和水平扫描的长度？为了取得完整的信息，在每个角度上，都要让整个样品都得到扫描。即，在每个角度扫描开始前，样品应该处于射线束的右边（Figure 5）；而完成时，应处于其左边。因此，扫描前最好：1) 尽量让样品的几何中心与转轴重合。2) 让样品在每次水平扫描起始和终止位置上至少距离射线束2mm。3) 利用节步骤

22、4)的转台调整界面，实验扫描位置，并确定水平扫描长度。4) 如果在扫描过程中，发现样品超出了扫描范围，最好能重新扫描。2.4.3 为什么扫描完成后要保存数据？保存数据。数据文件（投影值列表）可以被本软件再次读入，以备以后使用。比如，需要实验新的成像参数。2.4.4 为什么图像多出一条横贯全图的线？测试过程中发现，计数值有可能随机性地，非正常地偏低。这种现象出现的概率大概是0.1%的量级。相应地，该点的会偏高。对应在最后结果图上显示出一条横贯全图的直线。解决方法是使用3.2.6节介绍的修复数据功能：首先利用MATLAB SERVER输出的原始数据曲面寻找到这个坏点（计数明显偏低）。再把它的X,Y

23、坐标输入投影变换窗口，再修复数据即可。3 附录：CTSystem软件使用说明书3.1 概述CTSystem是专门为CT教学实验开发的软件系统。目的是控制硬件系统采集数据后分析，并最终绘制。它包括Turbo C下开发的MCA驱动、Visual C+下编写的步进电机驱动、用Visual Basic语言编写的用户界面和Matlab环境下运行的成像核心算法几部分。当然用户能看到的，只有界面而已。而在这个方面，我们尽量注重各个细节，力求方便使用。本软件固化了两种扫描方式：投影扫描和截面扫描。投影扫描类似于普通X光片，测量并显示，物体的吸收系数延某个方向投影（积分）到二维平面上得到的投影值分布。而截面扫描

24、指探测物体截面上的吸收系数分布，类似于医院用的CT机。3.2 界面介绍3.2.1 新建扫描项目启动CTSystem后，软件会首先要求用户输入成像的类型（投影扫描或截面扫描）、工作目录与工作的名称。如果在工作名称的文本框内输入了不能接受的特殊字符，软件会自动把它替换为下划线。之后点击确定后会在你指定的目录下建立工作目录，并准备响应的工具栏；如果点击取消&调整转台则会跳出Figure 8所示的界面。这个界面如下：Figure 7 新建扫描项目3.2.2 转台位置调整Figure 8 转台位置调整界面通过这个界面，你可以直接控制样品台的3自由度运动。你有两种方式控制转台：其一，用六个自锁按钮直接控制

25、。这种方法适合两个人合作时使用。其二，在“移动”框中输入所需的距离。这种方式比较精确，且可以一个人操控。三个轴的正方向分别是向右（从探测器看）、向上和从上往下看的顺时针方向转动。“已移动”框会显示相对于起始位置的位移。如果样品台已运转到极限位置，软件会弹出对话框告知。另外目前加载的步进电机驱动程序只能运行在WIN98下的。如果直接在NT系统下使用会报错（需要再增加附加的驱动WINIO.SYS）。3.2.3 调整能谱敏感区域我们不能直接把能谱的累计计数值作为某个点的透射强度。因为241Am放射源出射的射线包含多种能量。物质对不同能量的射线吸收度也不同。因此，把所有能量一概而论会造成最后图像的模糊

26、。Figure 9所示的241Am能谱中，最右边的主峰是59.5keV光电峰，之后的两个小峰是59.5keV射线的X射线逃逸峰。它们在一起代表了所有射入探测器的59.5keV光子数。其他的几个峰是放射源的低能X射线或者射线，一般不予以考虑。我们需要给出一个敏感区域以选择可用于计数的光子能量。这份工作将在本界面下完成。进入界面后，首先设置测量时间，再点击启动，指定时间后会绘出一张能谱（最好这个过程中不要放入样品）。之后，可以用A、D 可以控制红色的游标1；W、S控制绿色的游标2。两个游标所夹的区域就是计数的敏感区。保存敏感区域可以记录它并应用于之后的扫描。手工测试(DOS)指进入专用的DOS版界

27、面读取能谱。敏感区信息会生成一份相关的报告。Figure 9 调整能谱敏感区域3.2.4 扫描属性有两个窗口用于输入，扫描和成像中的必要参数。如Figure 10用于配置截面扫描（投影扫描窗口也大同小异）。其中需要说明的是：1) 角度方向插值数：这个参数描述角度方向滤波的大小。插值数越大，角度方向上，图像就越模糊，随之噪声也会变小。这种插值是针对与转轴同心的圆形的样品；对于其他形状，有可能会造成图像失真。2) 线度插值数：表示最后图像的宽度（象素），与X步数间的比值。相当于线度上的插值个数。3) 滤波函数：如1.4.1节所述，成像过程中需要一个实际可用的滤波函数。可供选择的有：HRL（见式）、

28、HSL（见参考文献1。它压低了图像的高频率分量，噪声更小但是会忽略一部分样品的细节结构）和Accumulation（实际上采用这种滤波函数，相当于使用早期CT累加法成像。关于累加法，具体请见参考文献1）。4) 插值函数：用于向插值的方法。具体内容请参考文献3。Figure 10 设置扫描参数3.2.5 扫描投影扫描和截面扫描使用了相同的扫描窗口，如Figure 11。扫描的同时，软件会把目前得到的数值点图形化，即以与数值相关的颜色标记在图像窗口内。在截面扫描下，横向代表水平方向；纵向代表不同角度。数值与颜色的关系由一条可以调整的数值-颜色映射曲线决定。这是一条三段线性函数。其他功能包括：1)

29、你可以使用开始、暂停、继续等按钮调节扫描进程。但是每条命令只有在上一个数据点处理完毕后才会生效。2) 详细信息：给出一份实时的报告，包括扫描的参数和所有的扫描结果。3) 保存（原始）图像：保存图像窗口的图片4) 电机降速：正常电机转速是每秒1000步；按下此按钮后变为每秒100步。如果机械部分发出过大的振动，可以考虑用这种方法改善。5) 低计数停机：探测器可能会因为过大的机械振动或电压振动而自动保护停机。这时会产生一系列的计数为0的错误数据点。为了防止这种情况的发生，可以启用低计数停机。这样在一个点的计数值小于指定值（通过滚动滑块）时，软件会自动停下来，等待处理。6) 后台文件的保存：正如前面

30、所述，每个数据点的能谱都会被储存在工作目录中。文件名是“工作名的前4个字符+2位水平坐标+2位角度或高度坐标.dat”；格式是纯文本文件。进度统计区域拖动以调整映射曲线点击这里可以调换颜色扫描数据的矩阵每个数据点的值以颜色的方式显示显示何种颜色由右边的颜色映射决定Figure 11 扫描窗口及注释3.2.6 投影变换窗口Figure 12 投影变换窗口如Figure 12所示，这个截面很简单，主要是连接后台的MATLAB环境（详见下节）。它的主要功能是：1) 其他选项的功能请见3.2.4节。2) 启动：启动MATLAB SEVER、注入扫描数据并计算。这个过程比较耗时，如果计算中点击窗口，回提

31、示MATLAB忙。3) 重新计算：修改其他选项后，用这个按钮重新计算。4) 修复数据：把指定的点的值改为周围的点插值。有些点的计数非正常地偏低。我们需要用这条命令手工修复它们。5) MATLAB命令行：启动后，会显示MATLAB COMMAND WINDOW。这里可以执行所有MATLAB语句。也可以用来修改算法。3.3 投影变换的输出投影变换窗口会输出近十张函数图和其他附加信息。本节将告诉你如何解读它们，并且如何利用它们改进自己的设置。还是以Figure 3(a)所示的实物为例。运行反投影变换中，程序将自动输出Figure 13到Figure 17，还有Figure 3(b)。Figure 13 原始数据（透射强度），即的二维分布。如果出现2.4.4节所述的异常，这里将会看到一个明显低于周围值的坏点，并且可以查到它的X、Y坐标（用于修复功能）。Figure 14 在扫描不同角度时，与随角度的变化（相对值）。如果，变化显著（如本图）则说明扫描过程出现了异常。比如本图中持续走低是由探测器工作状态漂移导致。Figure 15 （见式）与拟合曲线。同时信息窗口也会给出具体的拟合结果，如“the