医学图像配准

6 3基于基准点的配准方法 三维医学图像的配准算法可以基于从图像提取的特征 或基于图像强度的相似性最优化 首先从三维医学图像提取一个表面 然后在该表面上提取基准线 可能的话再从这些基准线上选取基准点 脊线是表面上的轨迹 或轨迹的集合 该处的主曲率在相应的主方向达到局部极大 脊线的提取要通过计算图像强度的三阶导数完成 涉及微分几何的知识 6 3 1极值线与极值点 极值线 ExtremalLines 极值函数e是曲率k在t方向的方向导数 局部绝对值最大的曲率点是极值函数e的零交叉点 极值线定义为e1或e2的零交叉点轨迹 脊线是正最大曲率极大及负最大曲率极小的点轨迹 极值点 ExtremalPoints 极值点定义为三个隐含数表面的交点 这些表面是其中 f是等值面 这些极值点对刚体变换具有不变性 6 3 2极值点的自动提取方法 第1步 提取等值表面 在小立方体元的内部提取等值表面 该等值表面 f I 与小立方体元的三个边相交 图a 第2步 提取脊线 用小立方体元角点e1值线性内插计算得到等值表面三角形顶点的e1值 如果这些e1值符号相同 则在此小立方体元内没有的极值线 否则 在e1符号相反的小立方体元的两条边上寻找插值后e1 0的点 P1 P2 联接P1 P2的线段就是近似的极值线 图b 第3步 计算极值点的位置 由于P1 P2两点都在小立方体元的表面上 可以用表面内e2值做双线性内插得到P1 P2点上的e2值 如果这两个e2值符号相同 则在此小立方体元内没有极值点存在 否则 可以沿线段 P1 P2 对e2值做内插 插值后该线段上e2 0的点P就是所寻求的极值点 图c a b c 6 3 3基于极值点的刚体配准 对同一个病人用相同成像模式但不同位置的两幅图像进行配准 可以从这两幅图像上分别提取极值线及极值点 再设法寻找这两组极值点间的对应关系 最后 计算能使这些匹配线段重叠的最佳刚体变换 一种典型的技术是采用几何散列 Hashing 表与Hough变换相结合的方法 6 3 4仅依赖基准点位置的刚体配准 在许多情况下可以通过人为地在目标表面加记号 或者用交互式方法或自动识别图像中明显的解剖标志点 几何图形明显转折点或图像灰度阶跃变化点等 统称基准点 如果这类基准点数目足够多的话 可以用来引导图像配准 其特点是仅依靠基准点的空间坐标而无须计算那些复杂的微分几何特性 两幅图像的配准可以通过函数f对参数c R和t的最小化求解 令 的奇异值分解 SVD 为UDVT 这里U和V是归一正交矩阵 D是对角矩阵 则 6 4倒角匹配图像配准法 倒角匹配 ChamferMatching 是一种对分割特征进行匹配的经典图像配准方法 后来在从卫星图像中寻找湖泊的轮廓方面得到有效的应用 该方法的基本思想是将线条图与图像匹配 6 4 1代价函数与距离变换 平均代价函数 平方根代价函数 上述方程是倒角匹配算法的核心 关于F 并不直接使用图像的灰度值 而是先对输入图像分割 然后再对分割的结果做距离变换得到的 在距离变换中 像素值给出到分割特征的最短距离 所以代价函数是两幅分割的图像特征间平均距离的精确估算 对此函数的最小化可以实现两幅图像的配准 倒角匹配算法的关键部分是距离变换 DistanceTransformation 距离变换的类型决定了使用该变换后代价函数描述两幅图像特征间距离的准确程度 距离变换是将二值图转换为灰度图像的一种变换 每个象素的灰度值是该象素与最近背景象素间的距离 常用的方法为倒角 Chamfer 算法 倒角算法通过对图像两次扫描 实现距离变换 1 一维距离变换对图像的一列像素做距离变换 设二值图像像素0 背景 1 物体 第一次 左 右扫描 背景点保留0 物体点为到最近0的距离 第二次 右 左扫描 物体点取值 min 第一次结果 物体点到最近0的距离 得到的结果是用灰度表示的反映该像素与最近边界距离的图像 2 二维倒角变换将模板与图像卷积实现二值图像到灰度图像的转换 3x3正向模板3x3反向模板 第一次扫描 用正向模板 沿图象从左向右 自上向下遍历 在模板每一位置 模板元素与图象上对应象素值相加 得到一组5个临时和 取其最小值作为该象素之新灰度值 模板空的元素什么都不作 第二次扫描 模板从右向左 自底向上 运算方法相同 两次扫描后结果即为图象的距离变换灰度图 此3x3模板产生的距离变换灰度值为笛卡儿距离的3倍 经规范化 与真实笛卡儿距离最大偏差8 距离变换示例 二值的哑铃形状图像经过距离变换转化成灰度图像 数值的大小代表该像素到物体边界的距离 左图为原始图象 右图为使用3x3的模板得到的距离变换图象 距离变换示例 二值图像转换为灰度图像 左图给出带有噪声的二值图像的距离变换 原始图像是包含较多噪声干扰的局部重叠的两个圆盘状对象 图 b 和 c 分别是使用3x3模板和5x5模板得到的距离变换结果 由于原图噪声过多 所以直接进行距离变换的效果不理想 对原图做闭运算处理 起到了平滑和去噪的作用 得到可以接受的距离变换图象 d 6 4 2图像分割与代价函数的优化 采用何种分割方法取决于实际应用 倒角匹配算法只要求将图像F分割成二值图像作为距离变换的输入 还要从分割得到的图像G产生一组轮廓点 轮廓点可用简单的轮廓跟踪法或对二值图像扫描搜索非零元素 线条图中点的顺序无关紧要 倒角匹配算法对所用的分割算法并无特殊的约束 从效率因素考虑总是希望线条中的点数目相对少些为好 倒角匹配算法对分割质量具高度强健性 这意味可以使用 低质量 自动分割技术 即算法可以是全自动的 即使图像质量差些也无妨 典型优化算法 单纯形法 Simplex 特点是容易实现 能够给出良好特性 该方法在M维空间中移动具有M 1个点的单纯形 M是参数的个数 对单纯形做反射 扩展或收缩去捕获最小值 Powell法 6 5基于最大互信息的多模医学图像配准 互信息是信息理论中的一个基本概念 通常用于描述两个系统间的统计相关性 或者是在一个系统中包含的另一个系统的信息的多少 在多模医学图像配准问题中 虽然两幅图像来源于不同的成像设备 但是它们基于共同的人体解剖信息 所以当两幅图像的空间位置完全一致时 其中一幅图像表达的关于另一幅图像的信息 也就是对应像素灰度的互信息应为最大 互信息只是一个相似性测度 如何准确计算互信息和实现不同模式图像间互信息的最大化 尚存在许多重要的技术问题需要解决 实施最大互信息配准法的一些重要技术问题 包括不增加新数据点的格点采样子集不产生分数灰度值的PV插值技术出界点策略在搜索策略上采用无须计算梯度的Powell算法 医学图像配准方法可分为前瞻性和回顾性两种 如果在成像前 在病人的颅骨附近固定标志点或定位框架 成像后通过对齐标志点使两幅图像配准 称为前瞻性配准算法 前瞻性配准方法一般说来配准精度较高 但它是有创伤的 而且操作比较复杂 前瞻性和回顾性医学图像配准 如果成像时不采取任何措施 在成像后应用某种算法寻找两幅图像的空间对应关系 则称为回顾性配准 例如特征点 面 法 不变矩配准法 相关配准法等 特征点 面 法原理简单 应用面广 但是需要较多的人工介入 配准精度受特征点 面 提取精度的限制 前瞻性和回顾性医学图像配准 6 5 1配准原理对于在不同时间或 和不同条件下获取的两幅图像A X 和B X 配准 就是要定义一个相似性测度并寻找一个空间变换关系 使得经过该空间变换后两幅图像间的相似性达到最大 即使图像A上的每一个点在图像B上都有唯一的点与之相对应 并且这两点应对应同一解剖位置 式中S是相似性测度 T是空间变换 配准的过程可归结为寻求最佳空间变换 由于空间变换包含多个参数 这是一个多参数最优化问题 一般由迭代过程实现 增量 T可以用计算梯度的方法计算 6 5 2以互信息为相似性测度互信息是信息理论中的一个基本概念 通常用于描述两个系统间的统计相关性 或者是在一个系统中包含的另一个系统的信息的多少 可以用熵来描述 熵表达的是一个系统的复杂性或者是不确定性 系统A的熵定义为 两个系统的联合熵为 其中 如果 表示已知系统B时A的条件熵 那么 与 的差值 就代表了在系统B中包含的A的信息 即互信息 因此两个系统间的互信息可以用下式描述 用联合概率分布 和完全独立时的概率分布 间的广义距离来估计互信息 对于离散的数字图像 联合概率分布 可以用归一化的联合直方图表示 边缘概率分布 表示为 边缘概率分布 表示为 互信息配准的关键技术问题 1 采样子集 Sub sampling 关于待配准的两幅图像 可以是二维或三维的 的定义域是整幅图像大小 对于两幅图像的联合直方图及每幅图像的独立直方图的计算应使用图像的全部数据点灰度值 理论上 也可以使用任一幅图像的子集或超集 Super sampling 但是 对于大多数断层扫描医学图像来说 3D体积数据集所包含的数百万 乃至上千万的数据点极大地增加了计算的负担 就目前的计算机速度来说 会使3D体积图像的配准无法实时实现和临床应用 因此 实用中都是只采用部分数据点计算图像的直方图 称做采样子集 如何选取部分采样点也有两种方法 一个是随机采样法 另一个是从原有数据格点中有规律的每隔几点选取一点 采样子集技术固然能几十倍 甚至几百倍的加快计算速度 但是本来用直方图近似概率分布的误差会因过少的数据点 使误差进一步扩大 实验表明 对于CT图像在XY平面上的分辨率为512 512 在Z轴方向采样28 34层 MR图像在XY平面上的分辨率为256 256 在Z轴方向采样20 26层的图像 采样因子fx fy 4 fz 2会取得满意的配准结果 为进一步减少计算量 可将XY平面上的采样范围限制在人脑灰度重心 Cx Cy 附近一定的区域内 例如200 200mm 对于CT MR配准 这种策略还可以明显减少图像背景伪迹对互信息计算的影响 2 插值技术 一些学者首先通过再采样 Resampling 使两幅待配准图像具有相同的空间分辨率 然后进行配准 但是再采样过程中必然要应用插值技术 即引入新的数据点 常用的三线性插值法按新采样点到其各邻点的距离产生相应的权重 新采样点的灰度值则为各邻点的灰度值按相应权重插值结果 插值产生的非整数灰度值给直方图的计算带来很大的问题 一般数字化的医学图像的灰度值均为整数 例如0 255 将插值结果整数化的过程必然引入误差 配准过程中 浮动图f中的样本在空间变换下对应的参考图r的像素通常落在非整数坐标上 因此需要再次使用插值运算 显然 两次插值对于完全依赖灰度计算的互信息而言将产生较大的误差 因此 可以直接应用原始浮动图的X Y Z三个方向上的等间隔的采样子集 由于不产生新的数据点 避免了一次插值 在配准过程中 浮动图f的样本a在某种空间变换下对应的参考图r的像素为b 通常b的空间坐标与任意一个实际的参考图像素并不重叠 此时可采用不引入新采样点的三线性PV TrilinearPartialVolumeDistribution 插值算法 该算法不是通过邻居点确定b点的灰度 而是按照周围8个像素与b点的空间距离分配权重 使周围8个像素点灰度贡献于联合灰度分布统计 即 其中 r i 是8个邻居点的灰度 wi是权重 这样就又避免了一次插值 与普通三线性按像素灰度值逐次加1的统计方法不同 PV插值算法在各灰度统计结果上每次增加的是分数权值 这不仅使互信息的计算更为精确 而且对于小的空间变换增量 T 互信息的变化 S也更平滑 对于优化过程中的局部极值问题会有所缓解 三线性PV插值 Tri linearPartialVolumeDistribution 当浮动图f中的某样本点Pf经过一定的空间变换T后的对应点Pr落在参考图r之外时 我们称Pr点为出界点 对于整幅图像来说 经空间变换后的浮动图f 由两部分组成 f fo fn其中 fo r Tffn a b c d 这里 fo是f 与参考图相互重叠的部分 fn是出界点区 3 出界点 Outlier 策略 显然 互信息的计算必须考虑出界点 将出界点忽略 即在不同的迭代周期内使用不同个采样点数来计算互信息 或将这些出界点的灰度近似为零 实验结果表明对配准精度都有不良影响 令出界点的灰度等于距其最近的边界像素点的灰度 这样做相当于扩大了参考图的背景 同时保持优化过程中的样本数不变 计算的互信息值更为准确 此外 如果浮动图的个别扫描层面不在参考图的扫描范围内 可以不对这几层数据进行采样 这样不仅进一步减少了计算量 而且可减少出界点的数目 将采样点有限度地取在感性趣区 ROI 附近 例如右下图中的S 经空间变换后的S 不产生出界点 6 5 3多参数最优化算法基于最大互信息的配准过程实际上是一个多参数的最优化过程 因此 如何选择合适的优化策略直接关系到配准结果的精度和速度 图像空间变换的形式通常可分为四种 刚体变换 仿射变换 投影变换和曲线变换 人脑可以近似为刚体 通过6个空间变换参数 包括沿着3个坐标轴的平移量 和分别围绕着3个坐标轴的旋转角度 即可确定两幅图像间的空间位置关系 一经确定了刚体变换的形式 配准过程就转化为寻找6个空间变化参数的优化过程了 对于每两幅待配准的图像 首先定义一个统一的立体坐标系统 X轴沿着行扫描方向 对病人而言是从右到左 Y轴沿着列扫描方向 从前到后 Z轴沿着层扫描的方向 从颅底到颅顶 坐标原点定义在图像的灰度重心 即 其中g x y z 是位于坐标 x y z 的像素点的灰度 这样做相当于将两幅图像进行了粗配准 有利于减少后续的优化步骤 而且可以避免那些远离全局最优点的局部极值 选择其中一幅图像作为浮动图F 另一幅图像作为参考图R 从浮动图的空间坐标 到参考图的空间坐标 的刚体变换可以用下式描述 其中 和 为3 3的对角阵 分别代表图像F和R的像素大小 和 分别是两幅图像的中心 是3 3的旋转矩阵 t是平移向量 采用Powell多参数优化算法 该算法由于无需计算梯度 可以加快搜索最大互信息的速度 在每一维内使用Brent算法迭代地搜索和估计配准参数 从而使互信息不断增加 在Powell算法中 令6个参数的初值为零 初始搜索方向为共轭的单位向量 考虑到成像过程中 病人在XY平面的平移和旋转比其它方向的平移和旋转更明显 令优化过程的搜索顺序为 经验证 这种搜索顺序效率最高 精度最高 应当指出 Powell算法容易受到局部极值的干扰 多分辨的金字塔方式可以较好的解决这一问题 在配准过程的初始阶段试用较粗的图像分辨以提高程序运行的速度 而后提高采样率提高配准精度 该方法可以较好地避免局部极值 6 5 4配准结果的评估 在医学图像配准问题中 通常没有所谓的 金标准 但是通过前瞻性的 基于标志点的配准方法可以得到一个近似的标准结果 在Vanderbilt大学医疗中心进行神经外科手术的一些病人 颅骨上被固定定位标记并接受多模医学图像 CT MR PET 数据采集 成像后通过配准定位标记点得到用于回顾性算法评估的 金标准 配准算法的研究人员使用的是已经擦掉标志点的3D多模图像数据 完成配准工作后 将所得结果提交Vanderbilt大学接受评估 评估之前 一些感兴趣区 VolumeofInterest 一般为10个 由医学专家给出 这些区域通常就是神经外科手术中的敏感区 每个感兴趣区被定义在MR图像中 同时计算其中心c 而后应用前瞻性配准算法得到的 金标准 确定其在CT上的对应点c 再用待评估算法的配准结果确定MR中对应c 的点c 通过计算每一个原点c与对应点c 的距离 作为目标配准误差 TargetRegistrationError 简称TRE 并由此统计出相应配准算法的精度 回顾性图像配准算法评估项目是一种 双盲 性的研究过程 所谓 双盲 即评估人员不知道被评估的具体算法 而算法的研究人员也不知道 金标准 直到提交所有的配准结果 这样就使得对算法的评估更加真实 可靠 并且更符合临床实际 实验所采用的所有图像数据来源于美国田纳西州的Vanderbilt大学的 回顾性图像配准算法评估 项目 该项目受美国NIH支持 编号1R01NS33926 01 J MichaelFitzpatrick教授是该项目的主要负责人 Vanderbilt大学提供的七个病人的全套CT和MR图像的3D体积数据 每个病人有1套CT数据和6套MR数据 包括PD T1 T2和分别矫正过几何失真的PD rectified T1 rectified T2 rectified图像 代号Patient 006的病人没有T1 rectified数据 CT图像在XY平面上的分辨率为512 512 在Z轴方向采样28 34层 像素在X Y Z三个方向上的大小分别为0 653595mm 0 653595mm 4 0mm MR图像在XY平面上的分辨率为256 256 在Z轴方向采样20 26层 像素在X Y方向的大小为1 25 1 28mm Z方向的像素大小为4 0mm 以MR图像中的像素对角距作为一个像素大小 即 应用了来自七个病人的CT MR共41套数据 经Vanderbilt大学评估 所有配准结果与 金标准 相比的误差都小于4mm 全部达到了亚像素级的配准精度 6 5 5实验结果 下表是Vanderbilt大学用目标配准误差的均值 中位数和最大值给出的评估结果 表中的每一列代表一种模式对 例如 CT到T1像的配准用CT T1表示 在对一种算法进行评估时 首先对各个模式对 分别计算七个病人的十个感兴趣区的TRE 而后按照模式对分别统计相应的误差均值 中位数和最大值作为对该算法的评估结果 之所以在列出均值的同时给出中位数 是为了减少少数出界点的影响 表中rect 表示所用MR数据在配准前已经经过了几何失真的矫正 配准结果与 金标准 相比的误差 单位 mm 第1层第5层第10层第15层 CT图像MR图像经过配准变换后提取CT图像的边缘 的CT图像并加MR图像的结果 最大互信息配准法具有如下优点 人工干预少 只依赖于图像本身的信息 不需要任何假设或先验医学知识 也不需要对图像进行特征点提取 组织分类等预处理 是一种自动而有效的配准算法 精度高 可以达到亚像素级 可靠性高 对图像中的几何失真不敏感 不依赖于成像设备 可应用于多模医学图像配准 多模医学图像配准具有非常广泛的临床应用价值 将多模图像结合起来进行分析可以为临床诊断和治疗提供更加全面而互补的信息 6 6医学图像配准的评估 医学图像配准 特别是多模医学图像配准结果的评估一直是件很困难的事情 由于待配准的多幅图像基本上都是在不同时间或 和条件下获取的 所以没有绝对的配准问题 即不存在什么金标准 GoldStandard 只有相对的最优 某种准则下的 配准 在此意义上 最优配准与配准的目的有关 6 6 1准标配准误差 立体定向框架系统 StereotacticFrameSystems 包括立体定向参考框架 立体定向图像获取 探针或手术器械导向几部分 优点是定位准确 不易产生图像畸变 使用立体定向框架系统的体积图像数据可以用来评估其它配准方法的精度 使用人工记号作准标的方法很多 一种准标 FiducialMarks 是使用9根棍棒组成的3个方向的N字型结构 在CT测试时 棒内充以硫酸铜溶液 作PET测试则填充氟18 这样 在两组图像中都可见此N字型准标 从而可对图像准确空间定位 例如用在人脑表面嵌螺丝作标记 每人8个 的方法对多个病人做CT MR T1 T2及PD 和PET实测 得到多组数据 这些数据专门用于多模医学图像配准算法评估使用 准标配准误差 FiducialRegistrationError orFRE 定义为 6 6 2目标配准误差 目标配准误差 TargetRegistrationError TRE 的定义与FRE相似 式中 T为三维空间变换 p和q是是一组临床相关的解剖点在两幅图像中对应位置 TRE在两个方面优于FRE 例如 Vanderbilt大学的配准评估是选取10个神经外科手术敏感区的中心作为用于配准评估的点集 因此更适合临床应用 FRE的准标选择受物理条件限制 只能选在颅脑表面 与临床关注区域远些 再者 FRE经常会高估或低估配准误差 6 6 3配准评估数据集 1 回顾性图像配准评估RIRE数据库 TheRetrospectiveImageRegistrationEvaluationProject 美国田纳西州的Vanderbilt大学的 回顾性图像配准算法评估 项目 该项目受美国NIH支持 又称VanderbiltDatabase 基于标记的回顾式图像配准评估 提供的七个病人的41套CT和MR图像的三维体积数据 包括每个病人的1套CT体数据和6套MR体数据 PD T1 T2和分别矫正过几何失真的PD rectified T1 rectified T2 rectified图像 代号Patient 006的病人没有T1 rectified数据 用于PET MR配准的三维体数据共35套 新增DICOM兼容的数据格式 2 非刚体图像配准评估数据库 TheNon RigidImageRegistrationEvaluationProject NIREP NIREP配准数据集由美国爱荷华大学 UniversityofIowa 提供 用于个体内部或个体之间解剖变形的比较 包括16名正常被试的三维MR脑图像体数据 全部受试是右利手 其中男性8人 平均年龄32 1岁 标准差8 8岁 年龄范围22 49岁 女性8人 平均年龄32 6岁 标准差7 5岁 年龄范围23 47岁 每套数据分割出32个三维的灰质感兴趣区 ROI 包括FP SFG MFG IFG OrbG PreCG PostCG SPL IPL STG ITG TP等 分别位于额叶 顶叶 颞叶和枕叶 扣带回及脑岛 不包括小脑 下丘脑和脑干部分 NIREP使用多种不同测度评价配准的性能 相对覆盖率 Relative