基于模板的颅脑图像基准点自动定位方法

王德兵

(安徽工贸职业技术学院计算机信息工程系，安徽 淮南 232007)

医学图像配准一直是研究者普遍关注的问题。现在，已开发出多种方法可用于医学图像配准。其中，基于标志点的图像配准在现实中最为常用，如在诸多商业图像导航系统中广泛使用此技术，而基于标志点的图像配准涉及不同图像(或物理空间)中对应点坐标的确定，以及使用这些对应点开展的几何变换估计，这些标记点可以是内在标志或外部标志。内在标志一般来源于自然存在的特征，如解剖标志点，而外部标志则为人工选择标志。基准点是基于标记的颅脑图像配准中最重要的参考指标之一，而颅脑图像中基准点位置的估算对于成功开展基于标记的图像配准十分重要。现有的商业计算机辅助手术软件系统要么不具备自动定位基准点的功能，要么声称具有此功能，但实际上并非如此。比如，Stealth Station(手术导航系统)就不具备此功能，史塞克(Stryker)在系统界面上设置“自动发现(Find automatically)”按键并声称具备此功能，然而实际显示所有的自动发现功能均不能正常工作，史塞克医疗器械有限公司承认确实不能使用这些功能，这就意味着基准标记仍需用户手动选择，用户必须观察所有图像切面并使用鼠标选择基准中心用于配准。在手动选择配准过程中，图像强度、对比度、切面间距，甚至是图像中不完整基准点均会使不同用户生成不同的鼠标定位结果。因此人工选择具有较大缺陷，并可能造成误差。自动基准点检测是图像配准的首要步骤，近年来研究者关于这一主题已开展了一些研究。文献[1]提出了一种基于强度的搜索和分类的自动标记定位技术，然而，这一方法需要关于标记大小和形状的先验信息，且仅适用于圆柱形标记。文献[2]开发了基于超声波模块的方法用于检测并定位皮下基准点，但该方法需使用附加硬件设备才能完成基准点的定准。文献[3]也提出了相应的方法，但要求基准点仅限于最佳尺寸大小。文献[4]提出了种子和条件膨胀方法，而该方法则需要高清晰图像。其他技术需要则需要增加硬件设备，如关节臂[5-6]、光学三角系统[7]、磁场数字转换器[8]等。本文充分利用图像空间的强度信息，提出了一种基于模板的基准点检测算法和决策系统，完成基准点的定准，克服了上述方法的缺陷。

1 检测系统结构

基准点检测过程如图1所示。

图1 基准点检测过程

通常情况下，使用扫描仪器获取的图像为轴位图像切面。切面间距(即两个连续切面间的距离)可根据模态和医生的要求固定，最常见的切面间距为2mm，一个图像序列可含有任意数目的图像切面。在大多数情况下，全脑图像序列含有75张切面，部分大脑图像序列含有25或52张切面，集中于肿瘤区域。一个图像序列可创建三维图像，即三维矩阵M(x，y，z)，M(x，y，z)的每个元素均是切面z中像素(x，y)的强度值。在大多数情况下，0≤x，y≤255， 0≤z≤25(或者0≤z≤50，0≤z≤75)。从M(x，y，z)中，提取了三种不同的图像视图：即轴位、冠状位和矢状位视图。轴位视图实际是输入图像序列，并对三个视图分别进行处理。在这一阶段使用图像处理技术是为了获取皮肤和基准点清晰的轮廓，如果一段轮廓的曲率具有基准点的特征，这一段轮廓被选为候选基准点。从三个视图中获取三个候选集，作为基于模板的决策系统的输入，再使用该系统测试轴位视图中各个候选基准点，以确定在相同的三维位置是否也有候选基准点。如果有，认为这个候选点是基准点，否则，这个候选点是边缘检测的误差。误差可能有诸多原因导致，例如，图像中不完整的基准点，或一段轮廓是从某个角度看起来像是基准点但却是皮肤的一部分。

2 方法

本研究关键是提出了从三维空间检测基于模板的自动检测基准点的算法，作为决策系统的核心，解决了基准点检测自动定准问题。通过大脑边缘图的构建，获取大脑轮廓；根据大脑轮廓中角点的不同曲率确定候选基准点；根据上述方法获取轴位、冠状位和矢状位三个视图的三个侯选集，作为决策系统的输入。本部分分别详细介绍了算法中使用到的模板，边缘检测方法，基于曲率的候选基准点检测方法和决策系统中基于模板的三维空间检测算法。

2.1 基准点模板

图2为研究中使用的基准标记模型。如将该模型放在头骨上进行扫描，可得到图3中所示的9个剖视图。基于剖视图，扫描平面会切割基准标记。例如，如果剖切面与z轴垂直，可得到剖视图3(g)。如果剖切面与x或y轴垂直，可分别得到剖视图3(h)和3(i)，其他6个剖视图是在剖切面与任一轴均不垂直的情况下获得的。

图2 基准标记模型

图3 剖切面在不同位置时的剖视图

将上述9个视图分为四类：视图(e)、 (f) 和 (g)属于第0类，视图(h)属于第1类，视图(i)属于第2类，所余视图属于第3类。本文方法仅使用第0、1、 2类视图。

对于图4中的三维图像，当对含有完整基准标记的连续图像序列开展边缘检测时，有以下两种可能(分别使用以下模板进行说明)。

图4 大脑顶部基准标记的多平面视图

模板1：至少有第0类的剖视图，如图3中的(a)和(b)。图4中第一行为基准标记的扫描图，第二行为所检测到的这些标记的边缘。所检测到的边缘为两个不被图像其他边缘完全包围的同心圆，这就保证了图像只能是头骨上设置的基准标记。

图5 基准标记轴位图及三幅连续切面中检测到的边缘

模板2：依次至少有一幅第0类剖视图，至少有一幅第2类剖视图和至少有一幅第1类剖视图组成的图像序列。在此图像序列前后，均无基准标记。在第一类剖视图中，边缘含有单峰，称之为“边缘类型1”；在第2类剖视图中，边缘含有双峰，称为“边缘类型2”；另外还有一种情况，边缘中无峰，称之为“边缘类型0”。这样，根据基准标记图像序列中边缘上的峰的数量，可分别从图6和图7中获得字符串“0111112221000” 和“01122222211100”。对于连续的相同数字，仅保留一个，将多余的都删除，以此方式可将图6和图7的字符串均收缩为“01210”。

图6 图3中基准标记的矢状位视图，13幅连续切面中检测到的边缘，以及边缘上峰的数目

图7 图3中基准标记的冠状位视图，14幅连续切面中检测到的边缘，以及边缘上峰的数目

以上模板1和2中，同心圆和峰值均在仅有一个基准点的区域内进行检测，研究中选择基准点最大可能尺寸为12mm。

2.2 边缘图构建

使用基于比值的阈值法[9]在原始图像中将前景对象从背景中分割出来，对于最佳阈值，两个区域像素数的比值会发生明显变化。从所定位的潜在阈值集中，可自动确定最终阈值的子集，第一个阈值将用于前景/背景分割。使用比值曲线进行分割的步骤如下：

生成直方图并使用高斯滤波器对其进行平滑处理，以抑制较小的强度变化。

计算比值曲线的二次倒数，获得归一变化率曲线，然后将曲线平滑处理，以消除较小的峰值和谷值，所得到的最大和最小值的位置决定了候选阈值的值，示例图像的分割结果如图8所示。

(a) 原始图

(b) 分割图图8 前景/背景分割

在分割后的前景中，使用递归形态侵蚀、膨胀和灰度重建技术融合较小间隙和消除多余噪声(常见于病例图像中)，然后使用Canny边缘检测技术以获取边缘[10](图9a)。

在检测到的所有边缘中，使用射线搜索程序选择大脑的轮廓：1)对每个边缘目标的像素进行编号并选择上半部分作为候选2)从图像四个角向中心设置四条射线3)与任一射线相交的第一条边缘即大脑的轮廓。在开展此程序前，需使用高斯平滑滤波器消除较小的物体和噪声(图9b)。

(a) 初始边缘图

(b) 轮廓选择图9 轮廓选择

2.3 基于曲率的候选基准点检测

大脑轮廓识别后，检测轮廓内的角点。单尺度特征检测不能同时检测细微和粗略特征，而多尺度特征检测[11]可解决这一问题。研究中使用自适应曲率阈值，而非单一的全局阈值。文中，曲率定义为：

在固定低尺度下计算曲率，以保留所有真实角点。曲率的所有局部最大值都被视为候选角点，包括错误角点，使用自适应局部阈值和角点角度移除初始候选角点中的错误角点和边界噪声。角点检测结果如图10所示。

图10 角点检测

根据角点位置，使用K-均值聚类将所选择的角点聚类至若干个类别中，再进行多项式拟合，搜索可将数据拟合至各个类别的三级多项式P(x)的系数，绘制各组的拟合曲线，如图11所示。选择各条曲线的中心点，投射至大脑轮廓上，所得到的点即为基准标记的中心(见图12)。

图11 曲线拟合

图12 最终结果

2.4 基于模板的三维空间基准点检测

从三个视图获取三个候选基准点集合后，使用决策系统确定基准点的位置。本节给出了决策系统及基准点检测算法的规则。对于从DICOM图像序列重建的三维图像“IMG”(包含N切面)，令I(a，i)，I(s，j)，和I(c，k)表示轴位切面图像i，矢状位切面图像j，和冠状位切面图像k。

三维空间基准点检测算法描述如下：

1. 定义Ma，Ms和Mc为与“IMG”具有相同尺寸的矩阵，且所有元素值设为0。

2. 令i=1。

3. 对I(a，i)开展边缘检测。

4.开展基于曲率的候选基准点检测。如果“模板1”存在，令Ma(x，y)=1，其中 (x，y)是同心圆中心坐标。

5.遍历整个大脑边界。对于所遇到的每个“边缘类型2”，对I(a，i-1)，I(a，i-2)，…， 和I(a，i+1)，I(a，i+2)等开展边缘检测，在此图像序列中相同大脑边界位置寻找峰数字符串，如果有“模板2”，令Ma(x，y)=2，其中 (x，y)是两个峰之间谷值坐标。

6.如果既没有“模板1”也没有“模板2”，并且不存在i=N，则令i=i+1，返回第一步。

7. 对Ms中所有值为1的元素，即Ms(x，y，z)=1，进行以下操作：

a) 令j=x，使用I(s，j)替换I(a，i)且使用Ms替换Ma，重复第3～5步操作；

b)令k=y，使用I(c，j)替换I(a，i)且使用Mc替换Ma，重复第3～5步操作。

8. 如操作7(a) 和 7(b)中均有“模板2”，则在坐标(x，y，z)处有基准标记。

9. 对Ms中所有值为2的元素，即Ms(x，y，z)，=2，开展以下操作：

a) 令j=x，使用I(s，j)替换I(a，i)且使用Ms替换Ma，重复第3～5步操作；

b)令k=y，使用I(c，j)替换I(a，i)且使用Mc替换Ma，重复第3～5步操作。

10.如操作9(a) 和 9(b)中有“模板1”或“模板2”，则在坐标(x，y，z)处有基准标记。

发现检测边缘以找到“模板2”的峰值并非易事，通过多种边缘检测滤波器如“Sobel”，“Laplacian”, “Prewiit”, “Zero-cross”, “Canny”等检测发现，其中使用“Laplacian”滤波器获取结果最优，所有结果均是基于此边缘检测算法。同时，尽管“Laplacian”在边缘检测中仍然不够完善，但通过考虑连续图像序列而非单幅图像，仍可以成功发现“模板2”并最终找到基准标记。

3 实验结果

尽管模板匹配方法在发现和定位标记中非常有用，但使用该方法时需将标记的基准点定义为其形心。该方法首先将标记从背景中分割出来，然后计算其强度加权形心。由于所使用的标记在CT和MR图像中较亮，在空气背景上无法成像。因此，分割时将会遇到以下问题：使用阈值法将图像分割为标记和背景体元，和与区域生长相关的标记体元。

因此本研究使用自适应阈值法RaBAM[12-13]。Sahoo 和 Glasbey对现有的阈值算法进行评阅和对比发现这些技术通常假设前景和背景体元来自两个主要模式。因为标记的尺寸与图像切面的厚度在同一数量级，许多含有标记的体元会遭受部分容积效应，即这些体元包含标记和空气的混合物。因此，图像在标记周围区域的强度直方图并没有两个主要模式。通常，强度直方图有一个对应背景噪声的主要模式和对应前景标记的广泛分布区域。如果含有组织或者CT图中的标志杆，背景模式也有可能是分散的。对于算法的第一部分，我们发现RaBAM算法可以胜任。

本研究使用六位病人的42套医学图像(每人一套CT和六套MR图像，正向和反向读出梯度T1，PD，T2)开发并测试了所提出的自动定位技术。每位病人使用4个标记，共计168幅标记图像，每套CT图像和MR图像平均处理时间分别为87s(时间范围82～91s)和18秒(9～48s)。

从24幅CT图像中识别出23个候选基准点，即算法对CT图像的标记误判率为4.66%。在144幅MR图像中选择具有最高平均强度的四个标记后，从144个真实标记中识别出131个，算法对MR图像的标记误判率为9%。

在含有45种MR模态和11种CT模态的18种案例下测试了本文的方法。测试中仅使用常规MR图像(切面尺寸：256*256或 512*512像素)和CT图像(切面尺寸：512*512像素)，所有的图像数据表示为DICOM图像序列，每个文件均为轴向扫描切面的图像文件。经测试，得出了理想结果。因文章篇幅限制，表1给出了两种案例的具体结果。

表1 案例图像信息和基准点测试结果

使用CT或MR扫描仪获取的图像质量的差异，并非所有模态下的基准标记均是可见的。例如，尽管MR的T1和T2模态下全部的12个基准标记均可见，但在案例1中的CT模态下，仅能通过视觉识别出9个基准标记。

并没有对CT模态下漏掉的3个基准标记进行重建，原因如下：1)认为所漏掉的标记在CT模态下并不存在，因为关于它们并没有可用的信息，如像素强度等；2)因为基准标记检测的目的在于实现自动配置而非相反目的，不能使用其他模态下的基准标记位置信息对它们进行CT模态重建。

表1中最后一行为实验结果。在两个案例中，本文方法均检测出了所有可见基准标记，没有基准标记漏检，检测到的标记位置均无误。

在全部55个模态下(MR和CT)，共计430个基准点(其中107个在CT图像中)，仅在一个MR模态下的一个基准点未能检测出来，未检测到的基准点的轴位视图如图13所示，该基准标记并没有位于图像的头骨部分，因此对此图像序列的边缘检测没有生成任何可识别模板。

图13 本文方法未能检测到的基准标记

4 结论

本文提出了一种全新的全自动基准点检测算法，可用于不同的图像模态中，经18个案例测试，实验结果准确，误判率很低，这个方法有助于实现基于基准点的全自动图像配准。该方法和人眼识别具有相同的性能，不受像素强度值，基准标记位置和方位的影响。后续关于基准点检测的研究将关注以下方面：将基准点检测与图像配准向结合； 测试基准点检测方法对iMRI图像的应用。