2017国赛A题省二等奖论文

1 CT系统参数标定及成像 摘要 CT Computed Tomography 可在不破坏样品的情况下利用样品对射线能量 的吸收特性对生物组织和工程材料样品进行断层成像 获取样品内部的结构信息 针对问题一 首先将附件1模板的几何信息数据及附件2中模板的接收信息 利用MATLAB软件进行处理 通过Radon变换及基于傅里叶变换的中心切片定 理进行X射线CT系统的二维重建 其次根据重建的二维图像 找出旋转中心与 几何中心之间的距离关系 利用MATLAB编程求出旋转中心与几何中心的距离 10 8003m 旋转中心的坐标 0 0 x y 为 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 xcos siny 做线积 分 映射成Radon空间的点 的变换 即只要知道了一个未知二维分布函数 的所有方向上的线积分 即可求出连续图像的Radon变换 根据简单的狄拉克函数 Delta函数 1 0 0 0 tt t 并结合直线方程可求出Delta函数的另一种表示 5 1 cos sin 0 cos sin 0 cos sin 0 x yx y x y 即在线l上的点 x y 满足 1x 其他非l上的点则为 0 x 综上 可将Radon变换表示为 cos sin D R f x y x y dxdy D为整个图像平面 f x y 为图像上某一点 x y 的像素点灰度值 特征函 数 为狄拉克函数 为 x y 平面内直线到原点的距离 为原点到直线的垂线 与x轴的夹角 Radon变换可以理解为图像在变换域 空间的投影 空间的每个点 对应图像空间的一条直线上的积分 也可以理解为图像顺时针旋转 角度之后在 水平轴上的投影 特征函数 使得积分沿着直线 xcos siny 进行 在 这里实际上就是直线的Hough变换 1 若求Radon变换的逆变换 即可求出被测 物截面函数每一点的 f x y 值 Radon变换的正变换和逆变换恰好对应了CT技 术的射线投影过程和图像二维分布函数的重建过程 程序见附录1 图1 对附件1 的0 90度采5个方向进行radon变换 而在实际重建的过程中 可依据中心切片定理 将所知的投影数据通过简单 6 的傅里叶反变换得到标定模板截面的评估 即平行投影的一维傅里叶变换等同于 原始标定模板的二维傅里叶变换的一个切片 2 首先定义代表截面函数的二维傅 里叶变换 2 j ux vyF u v f x y e dxdy 则某角度下一条投影 P t 的傅里叶变换为 2 j tF t e dt 若假设沿着 0 的直线方向 从平行投影的定义来看 可得到 0 x f x y dy 将 x y 坐标系统旋转得到 x y 坐标系统 如上图1 并将以上两式合并 借助旋转坐标定义可得 2 cos sin j x yF f x y e dxdy 综上所示 对 f x y 在与x轴成 1 2 k 角度的直线上的一系列投影逐个 进行傅里叶变换 可得傅里叶变换后 F x y 在与x 轴成 1 2 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 旋转中心的坐标 0 0 x y 为 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 266730d 2 要确定该CT系统使用的X射线的180个方向 首先假设每一个CT系 统使用的X射线均匀变动 投影最宽的位置与投影最窄的位置相差90度 且旋 转角度相差92个方向 则每个旋转方向的角度为 90 0 978392 选择以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 y d q F qt j qp dp 上式中方括号内是 q F qt 的1 D傅里叶反变换 利用傅里叶变换的卷积定 理可得 1 1 1 F q F qt F q F F qt 上式等号右边的第二项等于Radon变换 fR x y 1 ff x y d F q R qt 将 q 的1 D傅里叶反变换表示为 1 1 1 1 sgn sgn 2 2 qF q F q q F j q F j 利用微分性质 可将上式第2个等号右边的第一个反变换表示为 1 2 F j p 利用柯西主值 可将上式第二个等号右边的第2个反变换表示为 1 2 sgn 1 1 2 2 qF j p 因此 1 21 1 2F q p p 经整理得到Radon反变换 2 2 0 1 1 1 2 2fx y d R pt p 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 80L mm 单元探测器个数为 1 289n 此时探测器单元距离为 1d 椭圆最窄的 位置在第150个方向上 长度为 2 30 8 38L mm 30mm为椭圆的直径 8mm 为小圆形的直径 单元探测器个数为 2 137n 此时探测器单元距离为 2d 由 Ld n 求解 得出 1 20 2768 0 2774d d 利用最宽的位置球的数据更精确 存 在的误差更小一些 更接近于实际 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 A1 0 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 A6 100 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 0 radon change subplot 3 2 3 plot xt R 45 title 45 radon change subplot 3 2 4 plot xt R 90 title 90 radon change subplot 3 2 5 plot xt R 135 title 135 radon change subplot 3 2 6 plot xt R 180 title 180 radon change 附录4 补充 plot sum textA1 subplot 3 2 1 imshow textA1 xlabel Original image subplot 3 2 6 textA11 imrotate textA1 90 plot sum textA11 xlabel theta 90 subplot 3 2 2 plot sum textA1 xlabel theta 90 subplot 3 2 3 textA12 imrotate textA1 30 plot sum textA12 xlabel theta 30 subplot 3 2 4 textA13 imrotate textA1 45 plot sum textA13 xlabel theta 45 subplot 3 2 5 textA14 imrotate textA1 60 plot sum textA14 xlabel 22 theta 60 subplot 3 2 1 imshow textA1 xlabel Original image subplot 3 2 2 textA11 iradon textA2 10 imshow textA11 xlabel spac ing 10 subplot 3 2 3 textA12 iradon textA2 5 imshow textA12 xlabel spaci ng 5 subplot 3 2 4 textA13 iradon textA2 3 imshow textA13 xlabel spaci ng 3 subplot 3 2 5 textA14 iradon textA2 2 imshow textA14 xlabel spaci ng 2 subplot 3 2 6 textA15 iradon textA2 1 imshow textA15 xlabel spaci ng 1 r 1 512 subplot 1 2 1 x 62 c textA2 x plot r c subplot 1 2 2 x 150 c textA2 x plot r c xlabel direction 135 plot r c n 30 d 8 n Z Eabsorb h pcolor Z 100000 set h LineStyle none c h contourf Z set h LineStyle none