生物医学工程博士论文答辩.ppt

浙江大学博士学位论文答辩 答辩人 徐茂林导师 汪元美教授浙江大学生物医学工程与仪器科学学院2004 3 11下午2 00 论文题目 透射及散射断层成像中若干反演算法的研究ANewApproachtoReconstructionAlgorithmsforTransmissionandScatteringTomography 说明 本课题承蒙国家自然科学基金 No 69972047 资助 文中所用的实测数据由中国计量科学研究院提供 下图是中国计量科学研究院研制的工业CT机 工业CT机 本文立足于二维由投影重建图像理论 着重讨论了三大块内容 透射断层成像中变换法和级数展开法这两大基本算法的改进和创新 包括滤波反投影算法 FBP 的关键 新滤波函数的设计 图像局部重建的新方法以及FBP算法的误差估计等 对级数展开法中迭代算法做出改进 提出新的对称块迭代算法 在散射断层成像中 讨论了类似于CT中的滤波反投影算法 对S Norton提出的背散射成像的解析重建 论文主要内容 公式 基于散射的物理过程重新作了推导 使探测器响应函数有了明确的物理意义 针对散射问题由于射线 衰减 所导致的问题的 非线性性 及 不适定性 给出基于散射能量谱的行扫描Compton背散射成像模型 一个大型非线性方程组 用基于多目标优化的神经网络算法求得重建问题的解 此外 用滤波反投影方法给出了Laplace方程一般边值问题解的解析表达式 提供了势函数的一种快速求解算法 论文主要内容框架图 CT的数学基础 记被测物体对X射线的线性衰减系数为 根据Beer定理有 CT图像重建就是根据以上各个方向的线积分值 即射线投影 来确定物体内部的组织结构 即反演分布 A 反演算法 变换法用给定的一系列解析式 算子或函数变换 反演公式 A 来重建图像 常用Randon求逆公式的 离散化 形式 级数展开法直接对公式 A 离散化 得到一大型稀疏线性方程组 其中 为投影数据 为待求图像向量 级数展开法就是由给定的投影数据Y 估计出图像向量X 要求所估X与一组 个基图象 b1 b2 bJ 的组合能足够逼近所希望重建的图象 射线投影可由Radon变换 来描述 其中 各参数如下图所示 Radon运用平均值的思想 建立Abel积分方程获得求逆公式 Radon变换示意图 Radon逆变换公式中的广义积分是发散的 不能直接用来进行图象重建 需要正则化处理 由Fourier中心切片定理 Radon逆变换公式亦可写为 对滤波器作正则化处理 可获得CT图像重建的变换法 1 滤波反投影算法 filterbackprojection简写成FBP 其中 称为滤波函数 W R 称为窗函数 2 卷积反投影算法 convolutionbackprojection简写成CBP 滤波反投影算法 采用先修正 后反投影的方法 对每个投影函数作滤波处理 得到修正的投影函数 然后在对修正的投影函数作反投影运算 关键问题是投影函数的修正即 乘以权重因子 R W R 卷积反投影算法 用卷积方法修正投影函数 然后再作直接反投影重建图像的方法 优点 能快速实现 精度高 一般要求 完全投影数据 均匀采样 一 滤波反投影算法的改进及其它应用 新滤波函数及其时频特性分析 1 新滤波函数的表达式及时频特性曲线图 新滤波函数的频域曲线图 与几个常用滤波函数的频域曲线比较 A 4 A 4 A 2 A 2 黑色曲线 矩形窗红色曲线 新滤波函数蓝色曲线 S L滤波函数绿色曲线 Hamming滤波函数 时域特性曲线图 新滤波函数 矩形滤波函数 S L滤波函数 Hamming滤波函数 时域特性曲线比较 2 时域特性曲线图的数值分析 以下将精细地讨论新滤波函数的特性 给出其主瓣 旁瓣的幅值和宽度 并与常用滤波函数作比较 求解上式 可依次得出主 旁瓣的宽度和幅值 列表如下 各滤波函数主瓣 旁瓣的幅值和宽度 注 mk为滤波函数的k阶矩 定义为 讨论 滤波函数对重建误差的影响主要有以下三个因素 主瓣 临近旁瓣及远处旁瓣 主瓣越高越窄 空间分辨率越好 旁瓣越小 数值精度越高 通常密度分辨率就越好 从时域特性图中及上表中可以看出 新滤波函数的主瓣高而窄 尤其是远处旁瓣几乎为零 能够在保证空间分辨率的同时 相对地提高密度分辨率 3 数字实验结果 模拟数据 模拟数据重建图像评价测量值比较 实测数据 图3 1 9 采用R L滤波函数 图3 1 10 采用S L滤波函数 图3 1 11 采用Hamming滤波函数 图3 1 12 采用新滤波函数 局部重建 1 新滤波函数用于局部重建的特性 由于滤波函数是一偶函数 作分部积分后 有 根据Fourier分析中Riemann定理 有 其中 因此 当S较大时 成立如下的渐近近似式 对新滤波函数以及其对应的窗函数有以下结论 并且 于是得时域中新滤波函数的渐进展开为 根据 式 对于任何滤波函数均有 再由 式不难看出 其他常用滤波函数的衰减都不会超过2阶 综上 在s较大时 新滤波函数是2阶衰减的 而其他滤波函数由于H A 不为零则是1阶衰减的 因此 新滤波函数远处旁瓣的衰减最快 时域曲线局部区域放大比较图 新滤波函数的时域特性曲线图局部放大 Hamming滤波函数的时域特性曲线图局部放大 2 局部重建原理 选择适当的滤波函数 使得滤波函数的旁瓣迅速衰减 譬如在某一有限区间外很小 以至于几乎为零 则可以通过如下的公式作重建 是示性函数 与滤波反投影公式作比较 得出结论 采用满足如上条件的滤波函数 重建区域内任意一点的图像函数值 只需用到该点某一领域内的Radon投影数据 而不必用到全部投影数据 局部重建采用数据范围示意图 本文给出的新滤波函数即具有上述旁瓣迅速衰减的特性 可以根据上图所示的投影数据采集方式作局部重建 3 局部重建实验结果 模拟数据 Shepp Logan仿真图像 保持原图图像分辨率 放大256 256图象分辨率 注 数据量为301 200 局部重建区域半径为全局区域半径的0 5倍 实测数据 全局重建图像 保持原图图像分辨率 放大256 256图象分辨率 注 这是用陶磁叶片的实测数据进行的图象重建实验 数据量为241 576 重建区域半径为全局区域半径的0 36倍 时间约为2 3s 误差估计 1 对相当广泛的目标函数 算法有如下误差估计 分析上式可以得出 并计算得 因此 采用新滤波函数重建时 算法本身产生的误差是高阶小量 2 局部重建的二次误差估计 引进径向对称函数 故采用上述方法作局部重建产生的二次误差为 结论 CBP算法是CT中标准的方法 关于它的理论和算法已有很多研究 同时 人们认为CBP没有局部性 似乎不可能直接用来作局部重建 本文证明 无论从理论或数据实验来看 直接用CBP实现图像局部重建是可能的 所谓直接用CBP是指 就用CBP的计算公式 算法 以至于计算程序 但是 要恰当选取其中的滤波函数 本文选取一种新的滤波函数 它的旁瓣衰减很快 经证明在可能选取的滤波函数中 它的衰减是最快的 用新滤波函数 我们忽略远离计算点的数据 而实现图像局部重建 并给出了误差估计 本文建议的图像局部重建方法简单 快速 与全部数据用CBP重建的图像比较仍有较高的空间和密度分辨率 当ROI取的足够大时 重建的正是整个 上的图像 滤波反投影算法的其它应用 考虑Laplace方程的边值问题 其中D是R2中有界单连通区域 边界适当光滑 由于势函数u生成的梯度场是一势场 其积分与路径无关 因此 由边界条件则可以利用Radon求逆公式得到其共轭函数 利用Cauchy Riemann条件 选择适当的积分路线 得到势函数 M0为D中任一固定点 M为D中的流动点 作数字模拟实验 取如下模拟函数 对应的共轭函数的精确解为 在给定这些模拟函数满足的Laplace方程边值条件后 利用滤波反投影方法求得其共轭函数 记解为V 共轭函数精确解v 滤波反投影算法解v v v 误差灰度图 利用Cauchy Riemann条件求得相应的势函数u 用三维立体图分别描述出势函数及其共轭函数解 共轭函数解v 空间坐标图 势函数解u 空间坐标图 注 用滤波反投影法求得共轭函数的解较为精确 误差r为0 001829 计算所用时间也较短 求得空间分辨率为256 256的解的BMP图大约需1秒 二 级数展开法的若干研究 经典ART和SIRT方法的迭代格式 数字化图像的Radon变换可表示为 则图像重建问题描述成向量的形式 为投影数据 为待求图像向量 目的 寻找最优的图象向量 优点 适用于更普通的数据采集方法 如 沿曲线或沿区域 带状区域 的积分 容易利用先验条件 适于处理不完全数据问题 缺点 数据存储量大 重建速度慢 精度不高 1 加法ART算法迭代格式 2 SIRT法采用逐点校正 其迭代格式为 对称块迭代算法 1 原理简介 利用几何对称性简化运算并改变迭代方式 这种对称性描述为 即在平面XOY上任给一条投影射线 总可以找到另外七条投影线而形成一几何对称图形 这一对称射线组所穿过的八个对称象素所形成的图形成八卦状 见上图 结构十分简单 这一对称性可以用下面的群来描述 记分别是的关于轴和轴的翻转变换 分别是关于直线和的翻转变换 是恒等变换 构造如下集合 将中的乘法定义为变换的合成运算 那么由翻转变换的性质可知 是一个Abel群 而且是对称射线组的变换群 变换群不仅是保距的 而且是保像素格不变的 不计像素值 或者说像素格在变换群下仍然是像素格 只是位置不同 一般的变换群例如旋转变换群是保距的 但并不保持像素格不变 若将在变换群下的不变性看作平面上射线之间的一种关系 那么它是一等价关系 它将平面上所有的射线分成了等价类 每一类就是一对称射线组 于是 在平面上 每一射线都属于某一对称射线组 即所有对称射线组覆盖了全平面的射线 而且两个对称射线组 除对称轴外 不含同一射线 否则这两个对称射线组就是同一对称射线组 这样 就赋与了平面上的所有射线以简明的结构 便与计算及分析 这七条射线构造如下 若射线L1穿过的像素点为j1 坐标位置为 N是图象向量X的行 列 象素个数 则其余对称象素的位置分别为 利用这一对称性 同时可以得到八个Radon数据及伪投影数据 则迭代式中所需的计算元素均已知 而计算量和时间均节省了约七倍 大大加快了代数重建的速度 在这一基础上 提出对称块迭代算法 2 对称块迭代算法的迭代格式 3 对称ART算法 SART 方程的循环方式变为 即迭代是按对称射线组的顺序进行的 其中 是第k个对称射线组的指标集 是对称射线组的个数 M是总射线数 4 对称块ART算法 SBIRT 这里 而是关于对称块指标集的控制序列 即控制对称块的迭代顺序 表示第次迭代的权系数 3 数字实验 Shepp Logan仿真头像模型 对称块迭代算法重建结果 经典ART算法重建结果 SIRT算法重建结果 SART算法重建结果 重建误差及时间比较 注 重建图像大小为256 256 投影数据量为301 200 即在200个角度下取301条投影 4 实测数据 a SBIRT b SART c ART 这是用陶磁叶片的实测数据进行的图象重建实验 所给的平行投影数据量是576 281 图 a b c 是分别由算法SBIRT SART及ART重建的256 256T图象 迭代次数均为6次 迭代时间分别为95 57及301秒 实验环境与模拟数据重建实验相同 结果同样表明 对称迭代算法重建速度快 而且重建图象具有较高的分辨率 结论 基于平面上射线和象素的几何对称性 提出了对称ART算法及对称块迭代算法 SBIRT 模拟和实测数据实验表明 基于这种对称性结构的迭代算法 SART及SBIRT 使重建图象的各项性能指标都优于经典迭代算法 特别是速度快 值得注意的是这一对称性只是平面上射线和象素的几何性质 与重建图象函数的分布无关 因此 基于这一对称性的重建算法 对于所有重建问题都适用 传统ART算法迭代过程 对称SART算法迭代过程 对称块SBIRT算法迭代过程 三 背散射成像的多准则神经网络优化重建 Compton散射成像原理 Compton散射成像原理示意图 上图是Compton散射成像示意图 设有一单能光辐射源 其辐射强度为 发射的光子能量为E0 由该高能源辐射的光子经准直器后成为一束细光束照射到一样品上 光子与样品物质核外电子发生康普顿散射相互作用 散射后的光子经准直器 由探测器进行计数 探测器的光子计数率为 这里是入射光线方向上的小线元 dS是入射光束的横截面积 d 为探测器对散射点所张开的立体角 是Klein Nishina微分散射截面 为测量点的电子密度 其中的大小由放射源准直器的孔径决定 及的大小由探测器准直器决定 称为衰减因子 Ls和Ld分别为入射光线及散射光线的路径 此时 探测器的响应 即光子计数率 不仅与散射点的电子密度有关 而且还与入射光线和散射光线经过各点的电子密度有关 这样就使Compton散射成像问题成为一个 非线性问题 基于能量谱的行扫描背散射成像模型 行扫描示意图 设光源为单色高能 光 并选择广角理想点探测器 光源从左到右扫描 若不考虑多重散射效应及背景辐射 则入射光子经第i行物质电子散射后的能谱可写成如下形式 其中 是第i扫描行第k个能级的光子计数率 是探测器第k个能级的能量 表示第 i j 个象素的电子密度 N是被测截面图象在一维方向的象素个数 设分成N N的网格 是探测器在第i扫描行所能分辩的能级数 G是一个常数 它是由系统的几何形状和探测器的效率决定的 是Kronecker符号 表示第 i j 个象素对第k个能级有贡献时为1 否则为0 是第 i j 个象素的有效散射体积 是探测器对第 i j 个象素所张的立体角 若以表示探测器表面在散射方向的投影值 表示探测器到第 i j 个象素的距离 则 现将第i扫描行的电子密度及能谱分别记为向量及 则可写为如下矩阵形式 若记 即是由 构成的准对角矩阵 则成像的数学模型可表为 记J N N是重建图象的象素总数 是能谱的能级总数 可见是J维向量 是维向量 而是K J维矩阵 其元素记为 这是一个大型非线性方程组 重建问题即是要从所测能谱数据由上式反演出电子密度来来 多准则神经网络优化重建 如果考虑到多重散射效应及背景辐射等多种随机误差 成像模型应修正为 这里表示随机误差向量 由于指数衰减因子的存在 上面方程式投影矩阵的条件数一般很大 因此是不适定的 任何测量误差或模型误差在重建过程中将会被放大 为了获得稳定的解 必须用一些适当的先验条件对解空间加以限制 也就是说 我们只能求出该方程在某些准则下的优化解来 上述图像重建问题的正则化求解可表为在如下定义的可行域中 求使得下式 成立的图像向量 于是图像重建问题转化为一个多目标优化问题 该多目标优化问题可用线性加权法转化为如下单目标优化问题 即 由Kuhn Tucker条件可知 此单目标优化问题的最优解是上述多目标优化问题的非劣解 对图像重建问题来说 可选用以下三个目标函数 3 能量函数 它可使重建图象具有极小峰值 2 测量数据与再投影数据的加权平方误差函数 1 图像的负熵函数 这里是正则化常数 使最小等价于使 最大 它可保持重建图像的全局平滑性 是正则化常数 是正则化常数 是数据中噪声的方差 下面用人工神经网络 ArtificialNeuralNetworks简记为ANN 方法求解这一优化问题 人工神经网络是对人类大脑系统一阶特性的一种描述 它由大量的处理单元即人工神经元 通过适当的方式互连构成 是一个大规模的非线性自适应系统 简单地讲 它是一个数学模型 可以用电子线路来实现 也可以用计算机程序来模拟 是人工智能研究的一种方法 从数学的角度来看 这对应于一个微分动力学系统 即对应于一个常微分方程组 如果神经网络系统的能量函数还满足Lyapunov函数条件 则该微分动力系统的平衡稳定点即为能量函数的极小值点 现在将优化问题映射到一神经网络 相应的归一化图象变元映射为神经元输出变元 这一神经网络的能量函数由下式给出 这里J是神经元的总数 是增益函数 第一项是神经元之间的交互能量 第二项与约束惩罚有关 神经元受抑制的非线性输入输出关系可用n元函数表示 而 这里 时 否则 是一个较大的惩罚参数 而 第三项是调节项 网络演化的动力学特性由非线性常微分方程组来描述 可以论证该网络系统是渐近稳定的 即是上述动力学方程组的Liapunov函数 亦即 且当且仅当时 于是该网络系统是渐近稳定的 即网络最终收敛到它的稳定态 它对应于的全局最小值 这里与前述单目标优化问题的罚函数具有相同的形式 忽略调节项 因此当参数较大时 网络的稳态解即是优化问题的非劣解 实验结果 注 假设该头像模型所代表的截面大小为50mm 50mm 重建中被划分成64 64的网格 b 是多准则