2017国赛A题省二等奖论文.pdf

1 CT 系统参数标定及成像 摘要 CT（Computed Tomography）可在不破坏样品的情况下利用样品对射线能量 的吸收特性对生物组织和工程材料样品进行断层成像， 获取样品内部的结构信息。 针对问题一， 首先将附件 1 模板的几何信息数据及附件 2 中模板的接收信息 利用 MATLAB 软件进行处理，通过 Radon 变换及基于傅里叶变换的中心切片定 理进行 X 射线 CT 系统的二维重建；其次根据重建的二维图像，找出旋转中心与 几何中心之间的距离关系，利用 MATLAB 编程求出旋转中心与几何中心的距离 10.8003m ，旋转中心的坐标 00 (,)xy 为( 9.2734,5.5363) ；然后将求探测器单元 之间的距离问题转化为求小圆的模板直径长度与小圆直径所占探测器的个数之 比，得出探测器单元之间的距离为0.2667d ,对该问题进行模型改进后，求得 1 0.2768d ，为最优解；最后通过投影最宽与最窄的位置之间的关系，解出每个 旋转方向的角度为0.9783 ， 180 个方向起始点为-56.74， 终止点为 119.35。 针对问题二及问题三，通过 Radon 反变换，利用 MATLAB 将附件 3 及附件 5 的不同未知介质的接收信息进行处理，分别得到两介质各自在正方形托盘的位 置图像。由图以及第一问，知未知介质的初始位置均偏离几何中心，对其校正使 其旋转中心与几何中心重合，利用 MATLAB 软件编程，得到某未知介质的吸收 率三维图像。通过取边界点求标准方程法确认问题二中的未知介质形状为椭圆； 问题三中的图像明显为不规则图形。 最后将附件 4 的十个位置坐标带入各自的程 序中，用 MATLAB 求出了问题二、问题三中给定的十个位置的吸收率k。 针对问题四，依据对问题一中模板的精度及稳定性分析，重建长方形模板， 使CT系统的旋转中心与托盘的几何中心重合， 观察放射角度与投影长度的关系， 最终整理出的图像数据与实际探测数据基本吻合， 说明其精确度较高，从图像显 示出的线性关系也证明了重建模板的良好稳定性。 关键词：Radon 变换 傅里叶变换 Radon 反变换 二维图像重建 吸收率 2 一 、问题重述 CT 可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组 织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。整个发射- 接收系统绕某固定的旋转中心逆时针旋转 180 次，对每一个 X 射线方向，在具 有 512 个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减 后的射线能量，并经过增益等处理后得到 180 组接收信息。 CT系统安装时往往存在误差， 因此需要对安装好的CT系统进行参数标定， 即借助于已知结构的样品（称为模板）标定 CT 系统的参数，并据此对未知结构 的样品进行成像。 建立相应的数学模型和算法，解决以下问题： （1）在正方形托盘上放置了两个由均匀固体介质组成的标定模板，根据附件 1 中该模板的数据、及附件 2 中的接收信息，确定 CT 系统旋转中心在正方形托盘 中的位置、探测器单元之间的距离以及该 CT 系统使用的 X 射线的 180 个方向。 （2）利用附件 3 中某未知介质的接收信息和（1）中得到的标定参数，确定该未 知介质在正方形托盘中的位置、 几何形状和吸收率等信息， 并利用附件 4 相应的 数据文件具体给出图 3 所给的 10 个位置处的吸收率。 （3）利用附件 5 中另一个未知介质的接收信息及（1）中得到的标定参数，给出 该未知介质的相关信息，并具体给出图 3 所给的 10 个位置处的吸收率。 （4）分析（1）中参数标定的精度和稳定性。在此基础上自行设计新模板、建立 对应的标定模型，以改进标定精度和稳定性，并说明理由。 二 、问题分析 2.1 问题一问题一 本问要求我们通过给定的均匀固体介质组成的标定模板的信息，求解关于 CT 系统的一系列标定参数。而题中附件二中给出的是某模板的接收信息，我们 需要将其扫描信息转化为模板的几何信息，并建立平行投影重建的系统模型， 即 可确定 CT 系统旋转中心在正方形托盘中的位置等参数。因此我们首先引进 3 Radon 变换及基于傅里叶变换的中心切片定理， 解决 X 射线 CT 系统的二维重建 问题，得到标定模板的二维重建图。其次，进行 CT 系统投影旋转中心位置的确 定。但由于 CT 系统在安装时往往存在误差等因素，使得旋转中心与几何中心不 重合，找出两者之间的距离关系并作图，利用 MATLAB 编程求出旋转中心与几 何中心的距离从而求出旋转中心的坐标。 然后选取图中所得到的小圆形结合附件 1 中相应的数据，利用比例求解法，求出探测器单元之间的距离，但这种人工手 动查询的方法存在测量误差， 得到的结果不准确， 因此在模型的改进中进一步提 出相应的改进模型，使得出的结果精度提高。 2.2 问题二及问题三：问题二及问题三： 题中给出了附件 3 利用上述 CT 系统得到的某未知介质的接收信息，要利用 这些数据， 则需要对这些数据进行处理变换， 本问通过Radon 反变换， 通过MATLAB 编程求解，即可以得到某介质在正方形托盘的位置。由于其旋转中心与几何中心 并不重合，则要对其进行校正，对校正后的图像进行三维转化，建立某未知介质 的吸收率三维图像，对图像进行坐标分析，可得出该图形的大体形状。并将题中 所给定的十个具体点的位置坐标带入相应代码中， 直接得出这十个位置的吸收率。 第二问与第三问所要解决的问题本质是一样的， 该问针对附件 5 的数据信息进行 分析求解。对比第二问与第三问，这两个不同的未知介质的形状及吸收率都存在 较大差异。 2.3 问题四问题四 经过问题一问题二的分析， 我们知道了产生误差的主要原因。题目中要分析 其精度和稳定性，并自己设计一套方案。因此我们将 CT 系统的旋转中心尽可能 的做到与托盘的几何中心重合， 它是重建集合坐标的原点位置， 其误差会影响到 CT 图像上的伪影。重新建立的模板是个长方形模板，有利于我们去观察在不同 角度下与投影长度的关系，分析数据的结果与 CT 系统照射方向相吻合进而证明 了其精度，从得到图像的线性关系可以看出该次检测没有波动并全部覆盖，因此 稳定性也得到了保证。 4 三 、模型假设 1、假设可将标定模板的质心看作质点； 2、X 射线在照射过程中忽略光的波动及光的衍射对实验的影响； 3、假设每一个 CT 系统使用的 X 射线均匀变动； 4、模板及模型重建过程中，假设所选模板的旋转中心及几何中心是重合的； 5、假设实验所收集的数据能客观反映实际情况。 四 、符号说明 序号序号 符号符号 符号说明符号说明 1 平面内直线到原点的距离 2 原点到直线的垂线与 x 轴的夹角 3 D 整个图像平面 4 n 探测器单元的个数 5 d 探测器单元之间的距离 6 S 像素坐标 7 Y 叠加平均后得到的像素值 8 m 旋转中心与几何中心的距离 五、模型的建立与求解 5.1 问题一问题一 5.1.1 平行投影二维重建（Radon 变换及傅里叶变换重建二维图像） Radon 变换是指将( , ) x y平面空间中的一条直线l :=xcossiny 做线积 分，映射成 Radon 空间的点( , ) 的变换。即只要知道了一个未知二维分布函数 的所有方向上的线积分，即可求出连续图像的 Radon 变换。 根据简单的狄拉克函数（Delta 函数） ： 1,0 ( ) 0,0 t t t 并结合直线方程可求出 Delta 函数的另一种表示： 5 1,cossin0 (cossin ) 0,cossin0 xy xy xy 即在线l上的点( , ) x y满足( )1x ，其他非l上的点则为 ( )0x 。 综上，可将 Radon 变换表示为： ( , )( , ) (cossin ) D Rf x yxydxdy D 为整个图像平面， ( , )f x y为图像上某一点( , )x y的像素点灰度值，特征函 数为狄拉克函数，为( , ) x y平面内直线到原点的距离，为原点到直线的垂线 与 x 轴的夹角。 Radon 变换可以理解为图像在变换域( , ) 空间的投影，( , ) 空间的每个点 对应图像空间的一条直线上的积分， 也可以理解为图像顺时针旋转角度之后在 水平轴上的投影， 特征函数使得积分沿着直线 =xcossiny 进行。( , ) 在 这里实际上就是直线的 Hough 变换1。若求 Radon 变换的逆变换,即可求出被测 物截面函数每一点的 ( , )f x y值。Radon 变换的正变换和逆变换恰好对应了 CT 技 术的射线投影过程和图像二维分布函数的重建过程（程序见附录 1） 。 图 1 对附件 1 的 0-90 度采 5 个方向进行 radon 变换 而在实际重建的过程中， 可依据中心切片定理，将所知的投影数据通过简单 6 的傅里叶反变换得到标定模板截面的评估， 即平行投影的一维傅里叶变换等同于 原始标定模板的二维傅里叶变换的一个切片2。首先定义代表截面函数的二维傅 里叶变换： 2 () ( , )( , ) jux vy F u vf x y edxdy 则某角度下一条投影 ( )P t 的傅里叶变换为： 2 ( )( ) jt Ft edt 若假设沿着0的直线方向，从平行投影的定义来看，可得到： 0( ) ( , )xf x y dy 将( , ) x y坐标系统旋转得到 ( ,)x y坐标系统，如上图 1，并将以上两式合并， 借助旋转坐标定义可得： 2( cos( )sin( ) ( )( , ) jxy Ff x y edxdy 综上所示，对 ( , )f x y在与 x 轴成 12 , k 角度的直线上的一系列投影逐个 进行傅里叶变换，可得傅里叶变换后 ( ,)F x y在与 x 轴成 12 , k 角度的直线 上的值，按理论来说，当投影线束无限多，即可求得 ( ,)F x y在频域上对应线束 的对应直线的值，从而可利用傅里叶逆变换得到原图 ( , )f x y，如图 2（程序见附 录 1） 。 图 2 几何信息的图像实现 7 5.1.2 CT 系统投影旋转中心位置的确定 1、 据中国专利技术提供的一种 CT 系统投影旋转中心定位方法3， 给出其详 细步骤： a、采集 360 度投影数据，得到投影正弦图； b、设定合适阈值，对投影正弦图进行图像分割，将空气背景部分的像素值 设为 0，其余不变； c、对于 B 步骤得到的分割后的正弦图，将各个投影角度的投影数据按照对 应的像素叠加平均，即投影正弦图的各行叠加求平均得到一维信号 Y = F (S)，其 中，S 为像素坐标，Y 为叠加平均后得到的像素值； d、利用二次多项式 Y= As2+Bs+C 采用最小二乘法拟合 Y = F (S)； e、计算抛物线 Y= As2+Bs+C 的对称中心横坐标值 S = |，此即为投影旋转中 心坐标-2a 值。 2、 按原理来说，探测器的中心即为正方形的几何中心，CT 系统的中心也 应该与正方形的几何中心相同，但由于 CT 系统在安装时往往存在误差、实验仪 器对测量的影响等也能产生误差影响，使得 CT 系统的中心与正方形中心存在位 置偏差。 用中心探测器即第 256 个探测器所旋转的路径交点作为旋转装置的旋转中 心，给出旋转中心与几何中心的距离关系。 假设取模板的几何中心为坐标原点， 分别以水平、竖直方向为横纵坐标建立 如图 3 的二维平面直角坐标系， 图 3 旋转中心与几何中心位置图像 8 利用 MATLAB 编程（程序见附录 1）求得：旋转中心与几何中心的距离 10.8003m ，旋转中心的坐标 00 (,)xy为( 9.2734,5.5363)。 5.1.3 探测器单元距离及 X 射线的方向 1、 分别将附件1及附件2中的数据通过MATLAB进行处理(程序见附录3)， 可分别得到题中所给的两个均匀固体介质组成的标定模板的几何信息示意图及 投影图， 将两个图进行对比，并用 X 射线从不同视角进行照射， 可得到两介质旋 转过程中任意方向上的边缘曲线图，图 4 中左右两图分别给出的是采用 1 度、 135 度方向照射所得图像。 图 4 分别采用 1 度及 135 度方向 X 射线照射投影图 利用经过部分照射的小圆的表现来求距离， 采用不同视角进行多次验证， 并 结合图 4，进行游标移动，可得小圆直径方向上占用的探测器个数为 30。由题可 知，小圆的半径为 4mm,则其直径为 8mm。由于探测器单元之间的距离可等价转 化为求小圆的模板直径长度与小圆直径所占探测器的个数之比， 得出探测器单元 之间的距离为 8 0.2667 30 d 。 2、要确定该 CT 系统使用的 X 射线的 180 个方向，首先假设每一个 CT 系 统使用的 X 射线均匀变动，投影最宽的位置与投影最窄的位置相差 90 度，且旋 转角度相差 92 个方向，则每个旋转方向的角度为 90 0.9783 92 。选择以 x 正 方向为 0（第 58 个方向是 0，每个方向转 0.9783）逆时针转动，180 个方向起 始点为-56.74，终止点为 119.35 9 5.2 问题二问题二 5.2.1 Radon 反变换原理 Radon 反变换4给出从投影重建的解，完成从图像空间到其他空间的变换， 该过程可借助傅里叶变换进行。 因为 ( , )f x y可用( , )F u v的 2-D 傅里叶反变换表示，写成极坐标的形式为： 0 ( , )|()exp( 2)f x ydq F qtjqp dp 上式中方括号内是| |()q F qt的 1-D 傅里叶反变换。 利用傅里叶变换的卷积定 理可得： 111 |()| ()Fq F qtFqFF qt 上式等号右边的第二项等于 Radon 变换 ( , ) f Rx y 1 ( , )|( , ) f f x ydFqRq t 将| |q的 1-D 傅里叶反变换表示为： 1111 sgn | sgn 2 2 q FqFqqFjqF j 利用微分性质，可将上式第 2 个等号右边的第一个反变换表示为： 1 2 ( )Fjp 利用柯西主值，可将上式第二个等号右边的第 2 个反变换表示为; 1 2 sgn11 () 22 q F jp 因此 1 2 11 ( )() 2 Fqp p 经整理得到 Radon 反变换： 22 0 111 ( , )( , )( )() 22 f x ydRp tp p 5.2.2 原理应用 利用附件二中的数据，通过 MATLAB 编程求解，并用竖坐标表示其吸收率， 建立三维图像， 可得出某未知介质的一系列的标定参数，如确定该未知介质在正 10 方形托盘中的位置、几何形状和吸收率等信息，如图 5（程序见附录 2） ： 图 5 附件二吸收率三维表示 该未知介质的标定参数可用来检验依据附件三所求的吸收率的结果。 由第一问的 180 个旋转方向可知， 其初始位置偏离几何中心，对附件三给出 的某未知介质的接收信息进行处理， 通过上文中的 Radon 反变换应用，得到介质 在正方形托盘的位置如图 6 所示，对该未知介质的旋转中心进行校正，使其旋转 中心与几何中心重合于一点，利用 MATLAB 软件编程，得到某未知介质的吸收率 三维图像，如图 7。 （相应程序附录 2） 图 6 介质在正方形托盘的位置 11 图 7 附件 3 吸收率三维表示 通过提取若干边界点（不少于 10 个） ，求出边界方程，与椭圆的标准方程进 行对比可知，该图形的几何形状为椭圆形，题中图 3 给出选取的十个位置，依据 附件 4 给出的这十个点的具体坐标， 带入相应的程序中，即可解出给定的十个位 置的吸收率k分别为：-0.0109，0.0075，-0.0095，0.4932，0.4910，0.4939， 0.4896，-0.0178，-0.0009，-0.0033。 5.3 问题三问题三 按照第二问中同样的方法求解， 对附件 5 的数据进行处理， 得到另一介质在 正方形托盘中的位置，如图 8。 图 8 另一介质在正方形托盘中的位置图像 12 从图 8 中可明显看出，该图形与第二问中的图形形状差别很大， 第二问中的 图形比较规则轮廓比较圆滑有规律，而附件 5 给出的图形形状是非常不规则的。 该问中另一未知介质的吸收率三维显示如图 9。 图 9 另一未知介质的吸收率三维显示 按照附件 4 给出的 10 个位置坐标，带入相应程序（见附录 2）求出该十个 点的吸收率k分别为： -0.0289， 0.1189， 0.6502， -0.0112， 1.2797， 0.9415， 1.2114， 0.8662，-0.0181，0.0310。 5.4 问题四问题四 5.4.1 精度及稳定性的分析 根据问题一中的求解过程及数据分析，造成误差的主要原因是 CT 的旋转中 心不与正方形托盘的几何中心重合5，从而造成探测器接收到的信息与物体的实 际信息存在误差， 还有些外界因素的干扰比如噪音也会使最后测量出数据信息存 在偏差。 在测量物体的的长度时， 需要我们自己用肉眼先去观察， 之后进行测量。 而且测量的数据值较小并且单一，得出的答案也就比较容易存在误差，那么其数 据的精度也就受到了影响。 正方形的托盘在扫描过程中保持不变， 动的是 CT 系统 X 射线发射器和探测 器并且运动时他们俩相对静止， 在旋转过程中总共旋转 180 次， 但每次旋转的角 度可能会各不相同， 这样它旋转的总角度就得不到保证，有可能会导致物体不能 被全部照射上去，因此其稳定性不高。 13 5.4.2 新模板的设计及分析 1、鉴于对原模板的分析，现做出一种改进方案，以改进精度和稳定性。首 先要尽可能地使 CT 系统的旋转中心与托盘的几何中心进行重合，并给系统设定 一种运行方式即转动 180 次， 每次转动的度数为一度。 使其对某种介质进行扫描， 不断地调试， 利用问题一和问题二求解过程中得到的结论， 使两中心尽可能的重 合。最后得到的模板如图 10 所示。 图 10 新定义的模板图像 对其进行不同角度求解投影长度以验证其精确度，其中我们选取 0 度，45 度， 90 度这三个方向， 对探测器所接收到信息进行处理， 得到放射角度与投影长 度的关系。最终整理出的图像数据与实际探测数据相吻合， 如图 11， 从而证明其 精确度是准确的。 从图像中也可以看到每条线都是线性的所以它的稳定性也得到 了保证。 图 11 发射角度与投影长度的关系 14 六、模型的评价 6.1 优点优点 1、在解决同一个问题时，给出不同的解决方案，比如在求探测器单元距离 时，对原始求解模型进行改进，求出不同情况下的距离，最终对不同的方案进行 对比，选取相对精确、更真实可靠、更有信服力的一种方案的结果作为最终的答 案。 2、模型求解过程中，进行了合理假设，使求解过程简化。 6.2 缺点缺点 1、本文所用数据进行了精度取舍，结果存在误差。 2、在论文求解过程中进行的假设不一定符合实际，也可能对实验结果造成 误差。 七、模型的检验与改进 7.1 模型的检验模型的检验 在上文对二维图像进行重建时，实际采用的是基于 Radon 变换及傅里叶变 换的中心切片定理， 在重建过程中直接求出了其二维重建图即为变换角度为一度 时的图像，缺少不同角度下的对比，可能存在一定的偶然性及可信度。下面我们 通过反 Radon 变换，取不同的变换精度，分别作出以 10、5、3、2、1 度进行旋 转的图像，如图 12（程序见附录 4） 。 15 图 12 反 radon 变换精度比较（间隔分别为 10，5，3，2，1， ） 由图 12 可明显看出，在变换精度为 1 度时的图像最接近题中所给的标定模 板，并与上文利用中心切片定理作出的二维重建图像基本一致，因此，上文所用 模型求解结果是正确的6。 7.2 模型的改进模型的改进 7.2.1 改进前言 针对问题一中对探测器单元距离的求解问题， 一方面， 我们利用规则的小圆 形图形的几何长度及其投影建立 MATLAB 投影坐标系，在一定程度上减小了由 于模板图形不规则所带来的不可避免的测量误差等。但另一方面，由于所测量的 16 小圆形的尺寸比较小， 在单位像素点的取值上也存在一定的误差，且在投影坐标 上进行游标移动时采用手动取点的方式， 综上使得所求小圆的直径距离与两选取 点的横坐标之差存在误差量。因此，关于求解探测器单元距离，在下文中给出另 外的解答。 7.2.2 改进详解 通过 MATLAB 编程，假设有数据信息时为逻辑 1，反之则为逻辑 0，并利用 累加的方法，求出投影数据最宽的位置为第 58 列，且有 289 个像素点；在投影 数据最窄的位置为第 150 列，且有 137 个像素点,如下图 13 所示。 （图 13 二维吸收率图像） 结合题中给出数据并假设：椭圆最宽的位置在第 58 个方向上，长度为 1 80Lmm,单元探测器个数为 1 289n ， 此时探测器单元距离为 1 d； 椭圆最窄的 位置在第 150 个方向上，长度为 2 30838Lmm （30mm 为椭圆的直径，8mm 为小圆形的直径） ，单元探测器个数为 2 137n ，此时探测器单元距离为 2 d。由 L d n 求解，得出 12 0.2768,0.2774dd ，利用最宽的位置球的数据更精确，存 在的误差更小一些，更接近于实际。 17 八、参考文献 1 游福成.数字图像处理M.北京:电子工业出版社,2011,6:99-128. 2 张铁.CT 成像的改进 Fourier 算法极其实验J.CT 理论与应用研究, 2000,9(4).10:7-11 3 一 种 Ct 系统 投影旋 转 中心 定位方 法 P.中国 .测量 装置的制造及其 应用 技 术.G01N23/04GK102539460SQ20121000367.2012 4 石治赫.CT 扫描中的数学拉东(Radon)变换J.首都师范大学学报(自然科学 版).2013.34(4):15-18 5 康晓月.CT 重建算法的比较研究D.中北大学.2011 6 王璟.离散 Radon 变换卷积分配性及在从投影直接提取 CT 图像边缘的应用J.核电子学 与探测技术,2004,24(7):359-362 九、附录 附录 1：问题一程序： A1 = load(2.dat); A1=A1.100000; A2 = sum(A10); im, A3 = max(A2); in, A4 = min(A2); imagesc(A1) hold on plot(A3,A4; A3 A4,380 275; 512 512,w); set(gca,xtick, A4-180:30:A4:30:180) set(gca,xticklabel, A4-180:30:A4:30:180-A4) theta0 = -A4 disp(distance) d = 80/sum(A1(:,A3)0) A5 = find(A1(:,A3)0); disp(ordinate) y0 = (256-(max(A5)+min(A5)/2 )*d 18 A6 = find(A1(:,A4)0); A6 = A6(A6100); disp(abscissa) x0 = -(256-(max(A6)+min(A6)/2 )*d xb = -50 50 50 -50 -50; yb = -50 -50 50 50 -50; t = linspace(0,2*pi,90); xi = 15*cos(t); yi = 40*sin(t); xi2 = 4*cos(t)+45; yi2 = 4*sin(t); figure fill(xb, yb,49/256,225/256,223/256) hold on fill(xi,yi,234/256,114/256,48/256, xi2, yi2,38/256,22/256,230/256) plot(x0,y0,bo); plot(x0,y0,b.); plot(0,0,bx); axis image text(x0-5,y0+3,sprintf(Center of rotation, x0, y0) text(-5,-3,sprintf(Geometric center, 0, 0) 附录 2：问题二以及问题三： textA4=load(4.dat); d = 0.2778; nhalf = 512/2; xc = -33.5*d; yc = 20*d; phantom = load(5.dat);%or 3.dat 19 phantom = zeros(100,180); phantom; zeros(100,180); imagesc(phantom) figure img = iradon(phantom,0:179+30); n = size(img,1); x, y = meshgrid(-n/2:n/2*d); imagesc(x(1,:), y(:,1), img) axis image hold on plot(-xc,yc,ko) plot(-xc,yc,k.) plot(0,0,kx) xp = -50 50 50 -50 -50; x0 = -50 50 50 -50 -50; yp = -50 -50 50 50 -50; y0 = -50 -50 50 50 -50; xp, yp = rotxy(xp, yp, xc, -yc, 0); plot(xp,yp,-w) hold on I=imcrop; figure imshow(I) m,n=size(I); mn=(m+n)/2; I = imresize(I,200/mn); K=zeros(10,2); textA4=textA4.*2; K(:,1)=1,2,3,4,5,6,7,8,9,10; K(1,2)=I(textA4(1,1),textA4(1,2); 20 K(2,2)=I(textA4(2,1),textA4(2,2) ; K(3,2)=I(textA4(3,1),textA4(3,2) ; K(4,2)=I(textA4(4,1),textA4(4,2); K(5,2)=I(textA4(5,1),textA4(5,2); K(6,2)=I(textA4(6,1),textA4(6,2); K(7,2)=I(textA4(7,1),textA4(7,2); K(8,2)=I(textA4(8,1),textA4(8,2); K(9,2)=I(textA4(9,1),textA4(9,2); K(10,2)=I(textA4(10,1),textA4(10,2); disp(absorptivity) K 函数 1 function x, y = rotxy(x0, y0, xc, yc, angle) angle = angle/180*pi; x = (x0-xc)*cos(angle) + (y0-yc)*sin(angle); y = (y0-yc)*cos(angle) - (x0-xc)*sin(angle); 函数 2 function Eabsorb = ctradon(phantom, nhalf, angle) imrot = imrotate(phantom,angle,crop); imrot = imrot(abs(-400.5:400.5)nhalf, :); Eabsorb = sum(imrot,2); 附录 3：问题四： I = zeros(256, 256); r, c = size(I); I(floor(1/5*r:4/5*r), floor(2/5*c:3/5*c) = 1; figure; subplot(3, 2, 1); imshow(I); title(original); 21 R, xt = radon(I, 0:1:180); subplot(3, 2, 2); plot(xt, R(:, 1); title(0radon change); subplot(3, 2, 3); plot(xt, R(:, 45); title(45radon change); subplot(