基于局部特征射影变换的组织切片图像非刚性配准

赖明珠， 段志鸣*

(1. 海南师范大学数学与统计学院， 海口 571158； 2. 海南师范大学数据科学与智慧教育教育部重点实验室， 海口 571158)

图像配准的主要任务是寻找同一场景的不同图像之间的变换关系。这种不同可以是视角或时间的不同，也可以是由不同拍摄工具导致的模式的不同，最终得到图像中各个点的准确对应，实现图像间的信息融合[1]。

随着医学成像设备与方法不断地涌现，为了结合各种成像方法的优势，提高诊断的准确率，需要综合各个影像之间的信息，这使得医学图像配准技术在临床上具有重要的应用价值，也成为医学图像处理领域的研究热点[2]。

在病理分析过程中，医生往往需要分析患者不同染色下的组织病理切片图像来获得综合的判断，因此图像配准技术作为一种图像信息结合的技术，对于辅助医生诊断甚至恢复病理组织的三维结构起到了重要的作用。但是，因为组织切片图像中存在不均匀对比度、非线性的强度差、大量无纹理区域，以及图像之间具有局部非刚性的变形，使得组织切片图像配准依旧面临着很大的问题。

目前，主要的配准方法可以分为两类，区域配准法(全局方法)与特征配准法(局部方法)[3]。区域配准法以像素强度和最小均方误差、互信息、相位相关等优化目标为基础，广泛地应用于图像的时间配准以及多模态配准[4]。但是，对于组织切片图像中不均匀不一致的对比度以及大的无纹理区域，区域配准法并不能达到理想的效果。基于特征的方法假设在多幅待配准图像中都存在一系列可获得的对应的点，配准过程就是计算对应的点最大化相似度量。与基于区域的配准方法相比，特征配准法对于由于解剖差异以及光照造成的局部像素强度变化具有更好的鲁棒性，因此也更多地应用于医学图像的配准过程[5]。

在组织切片图像配准领域，Wang等[6-7]提出了一种综合了区域配准法与特征配准法优点的全自动的图像配准方法，首先应用特征配准法进行粗配准，之后应用区域配准法进行精细配准，该方法可以很好地解决因为形态失真、染色变化、染色伪影和组织损失导致的配准失效问题，但是对于不同的染色所获得的图像，该方法效果较差且不稳定；Obando等[8]使用K近邻方法根据灰度信息对图像公共信息进行提取，并通过B样条插值进行非刚性配准，该方法对于提升处理大数据量的图像具有很好的优势，且可以由粗到细进行配准，但是无法有针对性地对组织切片的局部变形进行配准；Trahearn等[9]提出了基于特征的配准方法，该方法在局部细化配准过程中以细胞核、脂肪粒等特征进行配准，因此不具有普遍性，适用范围有限；Liu等[10]在特征配准法的基础上，将特征点分为强匹配点与弱匹配点，并用强匹配点修正弱匹配点完成配准；Zhang等[11]提出了基于全局刚性配准的基础上，进一步进行全局的非刚性配准完成优化的方法，此外，也有学者应用深度学习的算法就特定的组织切片提出了专用的配准方法[12-13]。目前的研究方法主要聚焦于对粗配准结果的进一步细化上，虽然现有的基于图像灰度信息或特征信息的处理方法可以一定程度上消除图像局部变形的影响，但是未考虑到由于组织切片图像的特殊性，局部变形难以描述且各局部区域间的变换与图像整体的变换不一致等问题。

基于以上研究现状，考虑到图像各个区域之间存在局部的不确定性变形，现研究基于局部特征的图像全局非刚性配准。首先，对图像进行预处理，对预处理后的组织图像计算其尺度不变特征并进行配对；其次，根据配对的对应点坐标，应用随机抽样一致(random sample consensus，RANSAC)方法去除离群的配对点并计算图像对的全局射影变换参数；最后，在全局射影变换的基础上，采用K近邻算法将图像分成若干区域并对每个对应区域进行局部射影配准进而对全局配准进行修正，根据已有数据库中的手动标记的特征点评价配准方法的精度。

1 病理组织切片图像的预处理

显微镜下拍摄得到的组织切片图片中存在很大的问题无法直接应用，主要体现在对比度的不均匀、不一致以及大量特征稀少的无纹理区域这两个方面。尤其对于基于特征的配准方法，其主要依赖于组织的轮廓的特征，这使得突出组织轮廓结构显得尤为重要。所以要对图像进行预处理，在获得对比度较好且低噪声的组织切片图片之后对病理组织轮廓进行提取，为下一步特征提取做准备。

1.1 切片图像对比度均衡化

所获得的组织切片图片为红绿蓝(red-green-blue，RGB)彩色图片，如图1为将彩色图片直接处理为灰度图以及将三色通道分别处理为灰度图的结果，可以看出，绿色通道对比度更好，组织四周轮廓特征更加明显，所以选择绿色通道的图像进行后续分析。

图1 各个颜色通道切片图像

目前，常见的对比度增强方法包括空域增强与频域增强两种。空域增强属于一种直接增强的方式，对图像的色彩空间直接进行处理，而频域增强属于一种间接增强的方式，通过各种滤波手段对图像进行增强[14]。故采用空域增强与频域增强结合的方法对图像进行预处理。

图2中实线为原始绿色通道图片的直方图，可以发现，图像灰度十分集中，若采用全局直方图均衡化的方法，会使得特征区域图像灰度变化过大而失真，因此采用局部直方图均衡化的方法[14]。应用10×10的模板对图像对应区域进行直方图均衡化。经过局部直方图优化后结果如图2中虚线所示，图片的灰度在高灰度区域分布更加均匀，而低灰度区域改变较小，使得图像的组织轮廓特征与周围环境的对比度得到了提高的同时减少信息损失。

图2 直方图增强前后图像直方图

如图3(a)所示，经过直方图优化后，组织图像的对比度虽然得到了增强，但是图像背景也得到了增强，对图像特征的提取造成了干扰。为了进一步降低高频噪声的干扰，在直方图增强的基础上对图像进行滤波，因为组织切片最重要的特征为轮廓特征，而一般的高斯滤波会导致边缘信息的丢失，在降噪的同时破坏组织轮廓信息，而双边滤波可以在高斯滤波的同时，考虑图像像素值的影响，可以更好地保留边缘信息[15]。所以，采用双边滤波对图像进行频域增强。

图3 增强后组织切片图像

如图3(b)所示，经过双边滤波后，图像的边界特征被保留的同时，背景的噪声也被一定程度上消除，为下一步组织轮廓特征的提取提供了良好的前提条件。

1.2 切片图像轮廓特征的提取

经过预处理的图片，其组织轮廓特征与背景非纹理区域对比明显，但是背景中的非纹理区域对于下一步配准依旧造成了阻碍，所以需要针对组织的边界特征进行单独的增强。

在血管图像分割领域，Chaudhuri等[16]提出了匹配滤波对血管特征进行分割，故沿用这种方法对组织轮廓进行分割。

首先假设沿着组织边界截面的灰度分布符合高斯分布。由此可以构造沿着各个方向的高斯匹配滤波函数，对图像进行滤波，进而模拟组织轮廓的走向。进一步地，采用与组织轮廓走向角度一致的高斯滤波器所构造的模板与原图像进行卷积时，可以获得较大的卷积值。按照上述思路，可以构造一系列高斯匹配滤波模板对组织切片图像进行卷积运算，通过提取相应最大的点，完成对组织切片轮廓特征的提取。

高斯滤波函数为

(1)

式(1)中：L为假定的该方向下组织轮廓段的长度；θ为模板的旋转角度；σ为组织轮廓横截面的延展范围。因为假定组织轮廓方向沿着模板的y轴，所以模板要相应地旋转。

为了不改变原本组织切片图像的灰度特性，需要保证卷积模板各个系数的平均值为零，该模板均值为

(2)

所以，最终用于匹配滤波的模板为

K′θ=K(x,y)-m

(3)

针对所处理的图像，选择L为10～20，σ为8～20，并将圆周360°方向均分为12份，沿着12个方向对图像进行滤波。最终处理后的图像如图4所示，其背景区域基本去除，组织轮廓与主要纹理突出，可以进行下一步的特征提取工作。

图4 组织轮廓提取后图像

2 图像对的特征提取与配对

2.1 图像的特征提取

目前，常见的图像特征提取方法包括尺度不变特征转换(scale-invariant feature transform，SIFT)特征[17]、加速稳健特征(speeded up robust feature，SURF)[18]、Harris角点特征[19]等。采用SIFT特征作为配准所用的特征。

SIFT特征是一种广泛应用的特征描述子，它对图像尺度缩放与旋转具有不变性，对光照变化以及相机视点的变化也具有部分不变性，因此对于存在旋转变换以及光照影响的组织切片图片，提取SIFT特征作为图像配准的依据可以提供很高的准确性与稳定性。

提取SIFT特征的第一步为构建尺度空间，SIFT算法构建尺度空间的方式包括用不同尺度高斯核函数对原图像进行滤波以及图像的降采样两步，从而使得原图像细节特征逐渐减少来模拟大尺度特征的观察方式。

首先，一个图像的尺度空间定义为

L(x,y,σ)=G(x,y,σ)*I(x,y)

(4)

式(4)中：*表示卷积操作；G(x,y,σ)为一个可变尺度的高斯核函数，σ为高斯核；I(x,y)为输入图像。首先通过一系列不同尺度的高斯核函数对原图像进行卷积运算得到一系列尺度与原图像相同的图像作为高斯金字塔的第一层，之后对该层进行降采样即得到下一层，以此类推，构建高斯金字塔。

获得高斯金字塔后，为了在尺度空间内检测关键点，需要进一步构建高斯差分金字塔，SIFT算法所用差分方式为令同一层内相邻尺度空间的图像相减，对于某一层内第k个图片，可通过式(5)进行计算，即

D(x,y,σ)=L(x,y,kσ)-L[x,y,(k-1)σ]

(5)

尺度空间内的关键点为高斯差分金字塔内的局部极值点，为了寻找该关键点，每一尺度空间不仅要与其自身周围8个像素点进行比较，还要与相邻两个尺度空间内对应位置的9个像素点进行比较，从而获得尺度空间和图像空间的极值点。

为了使得描述符具有尺度不变性，需要对关键点分配相应的基准方向。对于在高斯差分金字塔中检测到的关键点，计算以其为中心3σ邻域内像素的梯度和方向，具体计算公式如下。

m(x,y)={[L(x+1,y)-L(x-1,y)]2+

(6)

(7)

之后采用梯度直方图计算关键点的主方向，梯度直方图将0～360o的方向分为36份，直方图的峰值方向即为关键点的主方向。为了使得该特征对光照以及视角变化具有部分不变性，需要进一步定义描述关键点的向量。首先，将关键点附近的邻域分为4×4的子区域，每个区域内将0～360°的方向分为8份，计算8个方向的梯度值，进而可以得到一个4×4×8的128维向量描述符用于下一步的匹配。

2.2 特征点的匹配

在获得两幅图片的特征点集之后，需要对两幅图中相应的特征点进行匹配。对于高维空间点的匹配，一般可以采用最优结点优先(best-bin-first,BBF)算法[20]进行预匹配。首先，假设图像I1中的所有SIFT特征描述符的集合为DES1，图像I2中的所有SIFT特征描述符的集合为DES2，则对于DES1中的任意一个向量des，都可通过式(8)获得一系列距离值，即

Dis={des·dess|dess∈DES2}

(8)

式(8)中：·表示向量的点积，该集合包含了des与图像I2中所有描述符的距离。其中，Dis内最小的元素所对应的dess为des最邻近的点。之后可以比较最临近点的距离与第二临近点的差距，如果最临近点比第二临近点的距离小得多，则认为完成了匹配，否则可以认为该点不匹配。

通过以上方法，最后可以得到图5所示若干匹配点，可以发现，这些匹配点中存在大量离群值需要去除，因此无法直接应用。

图5 BBF匹配的特征点

2.3 组织切片图像的全局射影变换估计

因为组织切片可以近似地看成是一个平面，所以应用平面射影变换模型是可行的。

设所获得的图像对中对应点的齐次坐标分别为p=[x,y,1]T，q=[x′,y′,1]T，则射影变换矩阵H满足：

zq=Hp

(9)

(10)

式中：z为p坐标的比例系数。

对于之前所获得的一系列匹配点，采用RANSAC算法估算变换矩阵可以剔除图5中的离群匹配点[21]。已知，对于上述方程，最少需要4个匹配点的坐标，所以首先随机选取4个点估计射影变换矩阵H，之后通过所得的H计算图像I1中所有的特征点坐标在I2中的对应点坐标，并计算这些点与之前已有所得I2中特征点的距离，统计在预先规定好的距离范围内的点，通过不断地迭代，找到匹配效果最好的射影变换矩阵H，同时剔除离群配对点，得到图6所示的匹配特征点。

图6 RANSAC提取后的特征点

2.4 组织切片图像对的局部射影变换估计

实际采集到的一系列组织切片图片，在采集以及染色过程中已经发生了局部的非线性变形，因此全局的射影变换配准方法无法达到对这些变形的精确计算，为此，需要一种基于全局射影变换的局部射影配准方法，用以降低因为局部变形导致的全局配准误差[22]。

在进行局部射影变换配准之前，首先对之前全局射影变换配准过程中通过RANSAC计算剔除离群值的特征点按照其所在的图像区域进行分类，进而将组织图片分成若干个区域，用以计算任意一个特征点在图片中所属区域。选用K近邻算法对图像中的特征点进行分类，以一幅图像中通过RANSAC计算剔除离群值的特征点作为训练集，通过计算各个特征点像素坐标之间的欧式距离将图片分成K个区域。如首先假设将切片图片分为3个区域进行训练，最终训练结果如图7所示，其中3种颜色分别表示特征点所属的3类区域，接下来图像中的任意一点可以根据其与图7中3类点的距离关系进行分类，由此完成图像区域的划分，目前尚无对图片分类数量确定的具体方法，预先假设所有图片的分类数K=3。

图7 对特征点进行分类(K=3)

图像完成区域划分之后，使用前文最优结点优先算法得到的特征点对，根据其所在的区域再次利用RANSAC算法计算其所属区域的射影变换矩阵，由此可以得到K个射影变换矩阵，即图像的变换矩阵为

H={H1,H2,…,HK}

(11)

对于图像I1中处于第i个区域内的点pi，其对应在图像I2中对应区域的点的坐标qi为

ziqi=Hipi，i=1,2,…,K

(12)

式(12)中：zi为点pi坐标的比例系数。

由于通过K近邻算法得到的区域内特征点的个数是未知的，若区域内特征点的个数过少，则所计算的局部变换矩阵会存在较大误差。而且，局部变换矩阵是根据全局变换矩阵进行的调整，所以全局变换矩阵H与局部变换矩阵Hi的差异一般不会很大，为此通过计算矩阵H与Hi之间的差的行列式值来评价全局变换矩阵与局部变换矩阵之间的差异，当差异过大时，将局部变换矩阵Hi替换成全局变换矩阵H，从而剔除离群矩阵。

3 方法评估验证

3.1 数据集

所用数据集来自2019年IEEE国际生物医学成像研讨会所提供的ANHIR数据集[23]，该数据集提供了不同类型组织(病变、肺叶、乳腺)的高分辨率连续组织切片图像，且一组中每个切片都选择不同的染色方式进行染色，并由领域内专家手动标记好有参考意义的特征点用以评价算法的配准精度。

选择小鼠的肺叶组织切片图像与乳腺组织切片图像用作本方法的评估，两类图像均使用蔡司Axio Imager M1显微镜获得。肺叶组织切片图像像素大小1.274 μm/像素，两幅待配准组织图像的染色方式分别为苏木精-伊红染色与cd31免疫荧光染色。乳腺组织切片图像像素大小2.294 μm/像素，两幅待配准组织图像的染色方式分别为苏木精-伊红染色与抗雌激素受体染色。每对图像中共有超过80个由领域内专家手动标记的对应点，用于评估算法的准确性。

3.2 实验验证

为了验证两种配准模型的配准效果，分别将数据集中第一张图(Source图)上的人工标记点的坐标代入已经得出的全局射影变换矩阵与局部射影变换矩阵中，并计算所得的对应特征点坐标与第二张图(Target图)中专家标记的真实对应点坐标之间的距离。

利用小鼠肺叶组织切片与乳腺组织切片序列用以评价配准精度，计算结果如图8和图9所示，其中，绿色连线为数据库中提供的由专家手动标记的两幅图中特征点之间的对应关系，红色连线为全局与局部算法计算出的Source图中的点与专家标记出的在Target图中的特征点的对应关系。之后，应用专家所标记出的在Target图中的特征点坐标与算法计算出的对应点坐标之间的距离值用于对配准精度进行评价。

通过对图8和图9的直接观察可以发现，由于局部区域之间存在变形，全局方法存在误差分布不均的情况，对于部分局部区域存在较大的误差，而采用局部配准的方法可以对这些局部区域进行修正，从而提高整体的配准精度。

图8 对肺叶组织人工标记点进行匹配

图9 对乳腺组织人工标记点进行匹配

进一步地，可以规定一个阈值，若误差距离小于该阈值，则认为配准成功。为了避免阈值选择的任意性，可以计算不同阈值下配准的成功率，构造出一连续的单调曲线，也为以后根据精度选择配准方法提出依据果。结如图10所示，其中图10(a)～图10(d)为4个不同的小鼠肺叶组织切片图像对的配准精度，图10(e)、图10(f)为两个不同的小鼠乳腺组织切片图像对的配准精度。

图10 不同组织切片图像配准成功率

3.3 实验分析

由实验结果可知，采用局部配准方法对于配准过程中的由于局部变形引起的误差起到了一定的修正作用，由于图像变形的随机性，所以对于不同图像的优化结果并不相同，对于局部变形较小的组织图像对，两条曲线几乎重合，说明各个区域之间的局部变换矩阵与全局变换矩阵差异不大，且由于选取的特征为尺度不变特征，该方法更具有普遍性。

4 结论

为提高病理组织切片图像的配准精度，在传统全局配准的基础上，通过应用已知特征点对组织进行划分并对不同区域进行局部配准，进而完成全局的非刚性配准。通过比较全局方法与局部方法在不同误差阈值下的准确率，证明了该方法对于局部存在变形的组织图像配准起到了误差修正作用。但是该方法依旧具有局限性：对于基于特征的方法，主要依靠足够的正确特征点进行配准矩阵的计算，所以未来可以通过进一步优化特征点的提取与匹配方法来提高精度。