1、一实验目的与意义医学成像技术是生物医学工程专业的一门重要的专业课程,课程主要涉及 X光仪器，CT仪器，MRI 仪器和核医学仪器的工作原理及成像方法。其中 CT算法的出现又为后来数字化医学成像技术的发展提供了基础。该门课程为生物医学工程专业的专业基础课。CT技术是医学成像系统中的一种重要手段。它通过特定的算法，利用计算机的高速运算功能，可以在短时间内快速呈现人体断层图像。让学生练习 CT图像的重建有助于学生理解 CT算法的内容，熟悉数字图像重建的过程。同时也能培养学生的团队精神和解决实际问题的能力。二实验原理利用滤波反投影图像重建算法，实现对 Shep-Logen头模型的重建，主要用到的原理有傅
2、立叶切片定理，快速傅立叶变换 FFT，以及滤波函数的设计。1.MATLAB处理数字图像的基本函数；2.X-CT三维图像重建的基本算法。三实验主要内容1.深入理解 X-CT成像的基本原理与图像重建步骤2.分析 X-CT图像的格式，并获取特定部的 X-CT图像数据3.研究 X-CT三维图像重建的基本算法（矩阵法、迭代法、傅里叶变换法、反投影法） ，其中重点研究反投影法（滤波反投影法和卷积反投影法）4.学习 MATLAB处理数字图像的基本函数；5.结合 X-CT图像数据与三维图像重建算法重建 XCT 图像四实验准备1.X-CT成像的基本原理、图像重建步骤及三维图像重建的基本算法相关知识X-CT（X-
3、ray computed tomography, X-CT）是运用扫描并采集投影的物理技术，以测定 X 射线在人体内的衰减系数为基础，采用一定算法，经计算机运算处理，求解出人体组织的衰减系数值在某剖面上的二维分布矩阵，再将其转为图像上的灰度分布，从而实现建立断层解剖图像的现代医学成像技术，X-CT 成像的本质是衰减系数成像。X-CT图像重建的原理源于 X射线通过介质时衰减的物理规律。根据扫描所获的投影值来求解成像剖面（实为断层）上衰减系数的分布，是选择数学方法的基本思想。理想单能窄束 X射线透射各向同性均匀连续介质时，强度衰减的物理规律符合朗伯（Lambert）定律 0xIe对公式两边同取对数
4、并整理可得1lnIx重建 CT像的重要环节就是从这一基本关系出发，通过对受检体的扫描，测出足够的投影值，再运用一定算法对投影值进行处理，确定各体素衰减系数 值的二维分布矩阵。X-CT三维图像重建的基本算法包括：矩阵法、迭代法、傅里叶变换法、反投影法，其中反投影法包括滤波反投影法和卷积反投影法。在实用 X-CT机中，图像重建的算法采取的是滤波反投影法。反投影法又称总和法，先通过 X射线束在体层的各个方向上对各体素进行扫描，测得一系列的投影值，然后，把各个方向的投影值沿原路径反方向投影回与原体素空间位置一样的体素上，得到该体素在各方向上反投影值总和，通过计算机运算，求出各体素 值，实现图像的重建。
5、2. X-CT三维图像重建的基本算法X-CT三维图像重建的基本算法包括：矩阵法、迭代法、傅里叶变换法、反投影法，其中反投影法包括滤波反投影法和卷积反投影法。矩阵法在实际中不常用，下面主要介绍反投影法、迭代法和傅里叶变换法。反投影法在实用 X-CT机中，图像重建的算法采取的是滤波反投影法。反投影法又称总和法，先通过 X射线束在体层的各个方向上对各体素进行扫描，测得一系列的投影值，然后，把各个方向的投影值沿原路径反方向投影回与原体素空间位置一样的体素上，得到该体素在各方向上反投影值的总和，通过计算机运算，求出各体素 值，实现图像的重建。下面用四体素（设 1234,）矩阵的重建对反投影法作定性说明。
6、对四体素矩阵作 0、45、90、135投影（即扫描） ，再将投影值反投回原矩阵的对应位置（即扫描通过的各个体素点）上，即可将原矩阵中的四体素的特征参数 值解出。其过程如下图运算中的基数（cardinal number）等于所有体素的特征参数的总和，这个总和也等于任一方向上投影值的总和。此算法由计算机执行。反投影重建的缺点是会出现图像的边缘失锐（即一种伪像）现象。下图定性地说明了边缘失锐的现象和产生此现象的原因。为了消除反投影法产生的图像的边缘失锐，在实际中采用的算法是滤波反投影法（filtered back projection） 。此方法是把获得的投影函数作卷积处理，即人为设计一种滤波函数，
7、用它对投影函数进行改造（卷积） ，之后把这些改造过的投影函数进行反投等处理，就可以达到消除星状伪影的目的。用一个滤波函数 h(x)对投影函数 (,)pR进行卷积计算的例子就是一种滤波方法。滤波效果的好坏取决于滤波函数形式的选择。滤波反投影法的另一优点是每一次投照结束，就可以通过计算机对投影函数作数学处理，待扫描结束后，数据的处理求解随之很快完成，所以图像重建的速度很快。迭代法迭代算法的一般原理迭代法是数值计算中的一类重要方法，它是求解方程或方程组的一类基本方法，在科学与工程计算中经常使用，尤其是由于计算机的普遍使用，使迭代算法的应用更为广泛。迭代法是一种逐次逼近的方法，它的基本思想是使用某个固
8、定公式反复校正值的近似值，从而得到一个近似根序列xk，使得该序列的极限就是函数 f(x)=0的解，求出的值是数值解而不是解析解。对于非线性方程 f(x)=0，用迭代法求根的具体做法是：先将方程 f(x)=0改写成便于迭代的等价形式 x=g(x)，然后取初值 x0没计算，xk+1=g(xk)（k=0,1,2） ，最后得到一个近似值 。x3.分析 X-CT图像的格式，并获取特定部的 X-CT图像数据X-CT图像像是灰度图像，一个 CT值应对应图像平面上某一级灰度。4.学习 MATLAB的基本使用方法及处理数字图像的基本函数MATLAB支持五种图像类型，即索引图像、灰度图像、二值图像、RGB 图像和
9、多帧图像阵列；支持 BMP、GIF、HDF、JPEG、PCX、PNG、TIFF、XWD、CUR、ICO 等图像文件格式的读，写和显示。MATLAB 对图像的处理功能主要集中在它的图像处理工具箱（Image Processing Toolbox）中。图像处理工具箱是由一系列支持图像处理操作的函数组成，可以进行诸如几何操作、线性滤波和滤波器设计、图像变换、图像分析与图像增强、二值图像操作以及形态学处理等图像处理操作。五实验流程图提取图象数据初始化产生R-L函数产生S-L函数进行radon 反变换显示重建图象结果结束六实验源程序及结果第一步 先利用 MATLAB里的 phantom变换生成头模型：1
10、.MATLAB程序p=phantom(256);imshow(p)title(头颅原图像)2.结果生成的头模型如下：第二步 先利用 radon函数生成 0到 180度 19个方向上的投影数据：1.程序%radon函数生成 0到 180度 19个方向上的投影数据theta1=0:10:180;R1,xp = radon(P,theta1); %radon变换 theta2 = 0:10:170; %总共有 10个投影角度;R2,xp = radon(P,theta2); %总共有 18个投影角度;theta3 = 0:5:175; R3,xp = radon(P,theta3); %总共有 36
11、个投影角度;theta4 = 0:2:178; R4,xp = radon(P,theta4); %总共有 90个投影角度;figure;subplot(2,2,1)imagesc(theta1,xp,R1) %radon变换图像colormap(bone); xlabel(Angle - theta (degrees);ylabel(Position - xprime (pixels);title(总共有 10个投影角度时);subplot(2,2,2)imagesc(theta2,xp,R2) %radon变换图像colormap(bone); xlabel(Angle - theta (
12、degrees);ylabel(Position - xprime (pixels);title(总共有 18个投影角度时);subplot(2,2,3)imagesc(theta3,xp,R3) %radon变换图像colormap(bone); xlabel(Angle - theta (degrees);ylabel(Position - xprime (pixels);title(总共有 36个投影角度时);subplot(2,2,4)imagesc(theta4,xp,R4) %radon变换图像colormap(bone); xlabel(Angle - theta (degree
13、s);ylabel(Position - xprime (pixels);title(总共有 90个投影角度时);2.运行结果3.关于 radon函数radon函数是用来处理计算制定方向上图像矩阵的投影。二元函数 f(x,y)的投影是在某一方向上的线积分，例如 f(x,y)在垂直方向上的线积分是 f(x,y)是在 x方向上的投影，在水平方向上的线积分是在 y方向上的投影。第三步 编写反 radon变换的程序：1.程序I1 = iradon(R1,20);I2 = iradon(R2,10);I3 = iradon(R3,5);I4 = iradon(R4,2);figure;subplot(2
14、,2,1),imshow(I1)title(10个投影重建);subplot(2,2,2),imshow(I2)title(18个投影重建);subplot(2,2,3),imshow(I3)title(36个投影重建);subplot(2,2,4),imshow(I4)title(90个投影重建);2.运行结果3.关于 iradon函数：iradon（R，theta）在二值数组 R中从映射数据重建图像。R 的每一行数据是平行束投影数据。irandon 假设旋转中心是映射的中心点，定义为 ceil(size(R,1)/2)。theta 描述的是映射中需要的角度。4.实验结果分析以上各图用不同部
15、分的逆 radon变换来重构图像。可以发现左上角第一幅图重构图像效果最差，因为它用来重构图像的投影最少，第二幅图重构图像效果比第一幅图好，第三幅图比前两幅好，第四幅图在四幅图中最好，因为随着旋转角度变小，投影个数增大，由于第一幅图重构图像投影太少，重构图像有很多的虚假点，为了避免这种情况，可以增加重构图像的投影角度的数目。第四步 与 MATLAN中自带函数 iradon做对比：1.程序：p=phantom(256);rethal=0:10:180;R1,xp=radon(p,rethal);ir=iradon(R1,rethal,hamming);figureimshow(ir)title(M
16、ATLAN中自带的 19个投影重建图象);2.运行结果：3. 实验结果分析由实验结果可以看出反 radon变换生成的图像内部细节不清楚。经分析不难发现，这并不是理论上有问题，而是在反 radon变换生成的图像里，头模型内部细节变化很微小，人眼无法辨别，故看不清内部细节。七实验结论参照实验指导书，选用反 radon变换，即在步骤 3生成图像的时候无法看清内部细节。但是这样的现象并不是理论上的问题，经过分析不难发现，在生成的图像 ir里边，头模型内部的细节变化很微小，一般在 0.01-0.03之间。变化本身就很微小的，当它反应到图像上的灰度级变化的时候，人眼根本无法分辨很小灰度级的变化，因此，看不
17、到步骤 3生成图像细节内部。为了能看到生成图像的细节内部，我们必须采用基于内插变化的滤波投影算法。而MATLAB中的 iradon函数，正是基于这个算法，才得出了具有明显细节变化的灰度图象。即步骤 4生成的图像。八、总结及心得体会通过对树种算法的实现，理解到在理论和计算机实现时，并不是那么简单，涉及到的内容不仅仅是滤波投影算法的理论，而且涉及到图像显示时，人眼结构对灰度级变化的敏感性。在这次的实验中，我收获很大也遇到了很多的问题，最终通过和其他同学讨论也解决了这些问题。通过这次的实验我了解了用 matlab实现三维 CT图像重建的强大功能和一些基本的重建算法，并基本理解了这些算法的原理及步骤，
18、为以后使用 MATLAB进行图像重建打下一定基础，对以后的工作也有积极的作用。通过合作、讨论也培养了团队精神和钻研问题、解决问题的能力，同时也增强了我面对困难的信心，为即将面临就业压力的我们起了积极作用。八、参考资料1.教科书 医用电子仪器及装备医学成像系统及放射治疗装置类 唐庆玉 主编 清华大学电机系生物医学工程及仪器专业2.医学成像系统 高上凯 主编 清华大学出版社3.医学图像处理与分析 罗述谦 周果宏 编著 科学出版社4.计算机图象处理的数学和算法基础刘丹 编著 国防工业出版社5. GT。赫尔曼著(严洪编译)，由投影重建图象，科学出版社 6.董雏申、吴世法、王天童，卷积反投影法从x光图像重建三维轴对称图象，全国图基科学会议论文集，l989。7. 赶荣椿等，数字图像处理导论，百北工业大学出版牡，西安，1995，p771798. 庄天戈，CT原理与算法)，上海交通大学出版杜，1992，p31 33
