CN102208109A - X射线图像和激光图像的异源图像配准方法 - Google Patents

Abstract

一种X射线图像和激光图像的异源图像配准方法，用X射线对物体照射，采集物体的X射线图像，该X射线图像是X射线透射物体得到的图像；用激光对同一物体进行照射，采集物体的激光图像，该激光图像是激光照射到物体表面反射得到的图像；所述配准方法的步骤包括：1)对采集到的X射线图像和激光图像进行预处理，来抑制图像噪声点的干扰；2)对激光图像边缘进行修正；3)分别提取X射线图像和激光图像的特征，并进行粗匹配；4)X射线图像和激光图像的精确配准。与现有技术相比，本发明先对规则形状的标定块的X射线与激光图像进行粗配准，然后进行精确配准，归一化互信息极值结果表明，X射线与激光图像配准精确有效。

Description

X射线图像和激光图像的异源图像配准方法

技术领域

本发明属于图像处理技术领域，具体是一种异源图像配准方法。所谓图像配准是指依据相似性度量准则得到图像间的坐标变换参数，使得从不同传感器、不同视角和不同时间获取的同一场景的两幅或多幅图像变换到同一坐标系下。

背景技术

近20年来，有各国研究者深入图像配准方法的研究，对于同源图像配准实用算法研究较为成熟和完善。Brown详细地总结了九十年代之前配准技术成果，Zitova详细总结了九十年代至今出现的相关技术和方法，这些图像配准研究的基本都是局限于同源图像配准的方法。相对于传统的同源图像而言，异源图像的最大特点在于其配准的是目标来自在不同传感器不同机理产生的特征图像，异源点之间的灰度值基本不具有共性信息，因此许多同源图像配准算法针对异源图像已无能为力，学者又针对各种异源图像配准开发的新算法，异源图像的来源主要是激光图像、可见光图像、红外图像、热红外图像、近红外图像，各种波段的光谱之间的配准，及医学图像上包括X射线、CT(Computed Tomograhy，计算机断层成像)、MRI(Magnetic Resonance Imaging，磁共振成像)、US(Ultrasonic Sound，超声成像)、CTA(Computed Tomographic Angiography计算机断层血管成像术)、MRA(Magnetic Resonance Angiography，磁共振血管成像术)、主要包括PET(positron Emission Tomography，正电子发射断层成像)、SPECT(Single Photon Emission-CT，单光子发射断层成像)等图像之间配准等。但基于激光与X射线的图像配准尚未有学者进行深入研究。

通过对各类图像配准方法的研究可以看出，通常不同的图像配准方法适用于不同类型的图像配准。尤其是异源图像配准，还没有任何一种通用的图像配准方法能适用于所有异源图像配准问题。

发明内容

为了解决现有技术中存在的上述问题，本发明提出一种X射线图像和激光图像的异源图像配准方法(本技术方案中：X射线图像是指物体经X射线照射得到的图像；激光图像是指物体经激光射线照射后，经CCD成像设备拍摄得到的图像)，具体技术方案如下：

一种X射线图像和激光图像的异源图像配准方法，用X射线对物体照射，采集物体的X射线图像，该X射线图像是X射线透射物体得到的图像；用激光对同一物体进行照射，采集物体的激光图像，该激光图像是激光照射到物体表面反射得到的图像；

所述配准方法的步骤包括：

1)对采集到的X射线图像和激光图像进行预处理，来抑制图像噪声点的干扰；

2)对激光图像边缘进行修正；

3)分别提取X射线图像和激光图像的特征，并进行粗匹配：

3-1)通过Hessian矩阵求取两类图像的特征点候选集合，提取合适的特征点；该特征点本质上就是图像的优选边、角等特征。具体求解方法如下表所示

3-2)对特征点所在局部图像的特征用向量形式描述，得到特征向量；

3-3)对两图像的特征向量和特征点分别进行两两比较，找出相互匹配的若干对特征点，建立两图像间的对应关系；

4)X射线图像和激光图像的精确配准：

用互信息表示的相似性测度对两图像进行精确配准，通过对完整的激光图像寻找一个仿射变换即旋转参数和平移参数，该仿射变换对应的X射线图和激光图像的互信息最大，此时X射线图像和完整的激光图像实现精确配准。

所述步骤1)的预处理方法是，制作平整的正方形标定块，通过中值滤波对标定块X射线和标定块激光图像进行噪声滤波预处理。

是所述步骤2)中，先对预处理后的激光图像进行二值化；再用数学形态学腐蚀运算，消除标定块激光图像边缘的毛刺；然后用膨胀运算，使图像边缘空洞的边界向外扩张，经过膨胀操作，把这些洞将被补上。

所述步骤3)中，先进行特征检测：

设图像I中，点)定义在尺度σ上的Hessian矩阵为：

其中， $g\left(\mathrm{\σ}\right)=\frac{1}{2{\mathrm{\π\σ}}^2}{e}^{-(x^2+y^2)/2{\mathrm{\σ}}^2};$

表示高斯滤波二阶导

与原始图像I＝(x，y)卷积；

表示高斯滤波二阶导

与原始图像I＝(x，y)卷积；

表示高斯滤波二阶导与原始图像I＝(x，y)卷积；

构建积分图像，点

处积分图像定义为由原点和构成的矩形区域内所有像素值的和：

所述二阶高斯滤波σ是通过在原始图像I上，扩大方框滤波模板的方框的大小来形成尺度不同的图像金字塔，得到不同的方框滤波模板值；

初始层的方框滤波模板值对应的尺度因子为σ；对应二阶高斯滤波σ、相应的尺度值为s＝σ；

继续求解Hessian矩阵的行列式值：

$\frac{{\left|L_{\mathrm{xy}}\left(1.2\right)\right|}_F{\left|L_{\mathrm{xx}}\left(0.9\right)\right|}_F}{{\left|L_{\mathrm{xx}}\left(1.2\right)\right|}_F{\left|L_{\mathrm{xy}}\left(0.9\right)\right|}_F}=0.912\≈0.9---\left(1.8\right)$

设

是矩阵的Frobenius范数，进一步求解得到Hessian矩阵的Δ表达式为

det(H)＝D_xxD_xy-(0.9D_xx)²    (1.9)；

再进行主方向确定：

设特征点所在的尺度值为s；

先以特征点为中心，计算以6s为半径的邻域内的所有点在x、y方向上的Haar小波响应，赋予与邻域内的点与特征点的距离大小成正比的高斯权重系数；

然后将每60°的范围内的响应叠加，来形成新的矢量，并遍历整个邻域，将最长矢量的方向定为该特征点的主方向；

最后进行特征点描述向量：

以特征点为中心将坐标轴旋转至主方向，按照主方向选取边长为20s的正方形区域，将该正方形区域划分成4×4的子区域；

在每一个子区域内，先计算5s×5s范围内的小波响应：

设d_x为相对于主方向的水平方向的Haar小波响应、d_y为相对于主方向的垂直方向的Haar小波响应，赋予这两个响应以权值；

再求每个子区域内的响应及响应绝对值之和：∑d_x，∑d_y，∑|d_x|，∑|d_y|，并求得每个子区域的四维分量的矢量V_sub＝(∑d_x，∑|d_x|，∑d_y，∑|d_y|)；对每一特征点，则形成的4×(4×4)＝64维的描述向量，然后进行向量的归一化处理；

经过以上计算，每一特征点都包含了三个信息：特征尺度值、特征点坐标和64维邻域向量，通过所记录的特征点信息，分别进行后续图像配准工作。

所述步骤4)中，互信息用熵来表示：

对于概率分布函数为p(a)的随机变量集A，其熵H(A)定义如下：

H(A)＝-∑p(a)logp(a)  a∈A    (1.10)

对两个图像，它们的离散的随机变量集分别是A和B，它们的边缘概率分布函数分别是p(a)和p(b)，联合概率分布函数是p_AB(a，b)，则随机变量A和B的联合熵定义如下：

H(A)＝-∑p(a)logp(a)log(b)  a∈Ab∈B    (1.11)

如果H(A/B)表示已知B是A的条件熵，那么H(A)与H(A/B)的差值，就代表了在B中包含的A的信息，即互信息；因此两个随机变量集间的互信息定义为：

I(A，B)＝H(A)+H(B)-H(A，B)

＝H(A)-H(A/B)

＝H(B)-H(B/A)    (1.12)

在图像配准中，当两幅图像的空间位置达到完全一致时，其中一幅图像A表达的关于另一幅图像B的信息，也就是对应像素灰度的互信息应为最大；用联合概率分布和完全独立时的概率分布间的广义距离来估计互信息：

$I(A,B)=\mathrm{\Σp}(a,b)\log \frac{p(a,b)}{p\left(a\right)p\left(b\right)}---\left(1.13\right)$

对于待配准的两幅图像，可以认为它们是关于图像灰度的两个随机变量集参考图像和浮动图像，a和b是两幅图像中相关的像素灰度值，a和b通过图像之间的坐标变换相联系，对于离散的数字图像，联合概率分布p_AB(a，b)可以用归一化的联合直方图表示：

$p_{\mathrm{AB}}(i,j)=\frac{h(i,j)}{{\mathrm{\Σ}}_{i,j}h(i,j)}---\left(1.14\right)$

边缘概率分布p_A(b)表示为：

$p_A\left(i\right)=\underset{j}{\mathrm{\Σ}}{p}_{\mathrm{AB}}(i,j)---\left(1.15\right)$

边缘概率分布p_B(b)表示为：

$p_B\left(j\right)=\underset{i}{\mathrm{\Σ}}{p}_{\mathrm{AB}}(i,j)---\left(1.16\right)$

则有：

$I(A,B)=\underset{i,j}{\mathrm{\Σ}}{p}_{\mathrm{AB}}(i,j)\log \frac{p_{\mathrm{AB}}(i,j)}{p_A\left(i\right)\·p_B\left(j\right)}---\left(1.17\right);$

接下来寻找一个合适的仿射变换，使得一幅X射线图像经过此仿射变换后和激光图像的互信息最大。

所述步骤3)的先进行特征检测中，关于方框滤波的说明如下：

采用方框滤波近似代替二阶高斯滤波，采用积分图像加速卷积以提高计算速度。在原始图像上，通过扩大方框的大小来形成尺度不同的图像金字塔，9×9的方框滤波模板值如图3所示。

图6中色底数字表示盒滤波器模板大小，上述9×9滤波模板即为尺度空间的初始层，后续层通过依次放大滤波器的大小来计算。

滤波器大小扩展为15×15，21×21，27×27等。根据3D邻域内求极值点的运算可知，每组中初始层和最末层不包含极值点，因为它们只作比较用。其它组的滤波器大小类似构造，对每个新组，滤波器大小的步长成倍增加，依次为6、12、24、48，N×N大小的滤波器对应尺度s＝σ＝1.2×N/9。将尺度空间的每个点与本尺度的相邻位置以及相邻尺度的对应位置共26个邻域进行比较，得到的局部极值为侯选特征点，然后在尺度空间和图像空间中进行插值计算，得到稳定的特征点位置及尺度。

与现有技术相比，本发明先对规则形状的标定块的X射线与激光图像进行粗配准，然后进行精确配准，归一化互信息极值结果表明，X射线与激光图像配准精确有效。

附图说明

图1是X射线图像与激光图像由粗到精的配准过程示意图；

图2是粗配准算法特征匹配流程示意图；

图3是9×9方框滤波模板(其中灰色部分模板值为0)；

图4是实施例步骤3-1)中尺度空间方框滤波大小示意图；

图5是本方法应用场景的激光器和X射线机的几何位置示意图，其中1是激光器、2是CCD相机、3是X光发射器、4是被测物。

图6是尺度空间盒滤波器大小示意图。

具体实施方式

下面结合本方法在肉图像异物提取前需图像配准为例，对本方法作进一步说明：

基于X射线和激光图像特征及互信息的由粗到精的图像配准

1、X射线肉图像和激光肉图像特点

肉图像异物提取前需图像配准，特征图像源来自X射线相机与激光相机，两种图像虽然来自为常见的传感器，但由于两者的成像机制不同，适合肉的图像配准方法的比较难选择，肉图像配准的难点是如何找出肉X射线图像与激光图像所必须的匹配特征对，以确定图像变换模型。

X射线图像及激光图像中肉图像由于成像机制差异灰度对比度并不显著，此外肉X射线图像与激光图像的肉灰度与背景之间的对比度也有差异。在X射线图像中，肉与背景对比度相对较高，肉灰度值由肉的厚度和密度决定，灰度对比度可以在相对大的范围内变化；在激光图像中，肉灰度对比度完全由肉的厚度决定，灰度对比度相对X射线图像在较小的范围内变化；再者肉表面较为光滑，肉X图像与激光图像灰度梯度变化平滑，肉图像中没有明显的区域特征点，基于特征点的图像配准方法受到很大的局限；另外激光图像在运动过程中采集，受到系统的热噪声、震动噪声等的干扰，激光采集的肉图像相对X射线比真实肉在成像形态上误差大，且图像中肉边缘与背景的灰度差异不大，有时很难得到强的闭合完整边缘，肉边缘形状非规则，对肉边缘提取需要大量的运算时间，从射线图像、激光图像中提取肉一致特征困难更大，误匹配可能会出现。

本方法采用各种算法对肉X射线图像及激光图像进行配准大量试验，大多数现有图像配准算法用于肉X射线/激光图像上的效果并不好。例如，基于角点特征的hariss配准方法，仅适用于边界分明或灰度级别特征变换较大图像，该算法不能直接处理肉X射线图像及激光图像。

对于本方法采集系统中固定视角的X射线和激光图像的配准来说，两个视场的相对位置关系固定，只需做一次配准得到固定的平面变换关系后，X射线和激光图像可以自动实现配准，因此本方法中X射线和激光图像配准的对象将灰度特性上不显著的肉对象改为形状规则的标定物来进行校正，算法需求主要体现在配准精度和可靠性等方面，配准速度及实时性并不重要，从而减少配准算法的强度要求，达到X射线和激光图像配准的目的。

2、基于特征点与互信息相结合的由粗到精的配准方法

综合考虑上述因素，本方法提出采用一种基于X射线图像及激光图像特征点与互信息相结合图像配准方法。

首先选用已知规则且表面平整的标定块，通过对基于X射线图像及激光图像标定块特征点的进行配准，图像配准误差较大仅是粗配准；再选用有曲面变化标定块，将经过粗配准的X射线图像及激光图像利用互信息配准的方法进行图像的精确配准。

假设表面平整的标定块A，其激光图像为参考图像R_A的主轴方向即方向角为d_R，浮动图像F_A的主轴方向角为d_F，则两图像间的相对旋转角为：

θ_A＝d_F-d_R    (1.4)

设参考图像R_A和浮动图像F_A的质心坐标分别为C_R(x_r，y_r)和C_F(x_f，y_f)，则两图像在x，y方向的平移参数tx_A，ty_B、用下式计算：

tx_A＝x_f-x_r

(1.5)

ty_B＝y_f-y_r

可得到上述粗配准的变换参数T_A＝[θ_A，tx_A，ty_A]，同样有曲面变化标定块B，其激光图像为参考图像R_B，浮动图像为F_B，将变换参数T_A应用到浮动图像F_B，得到粗配准后的浮动图像T_A(F_B)，将参考图像R_B和经过粗配准的浮动图像T_A(F_B)按照基于互信息的图像配准流程进行下一步配准，得到配准参数T_B＝[θ_B，tx_B，ty_B]，最终实际配准结果为T_A+T_B。基本流程如图1所示。

粗配准中特征点集中了图像上很多重要的局部特征，具有旋转不变性优点，但激光图像及X射线图像存在干扰，容易丢失图像数据信息，配准离真正的配准位置还有偏差，互信息区域精确配准，是在变换参数的寻优过程中，每进行一次搜索，都重新计算一次互信息，以得到精确的配准图像。

3基于图像特征点的粗配准算法

过程包括：(1)通过Hessian矩阵求取图像的特征点候选集合，提取合适的特征点；(2)对特征点局部特征用向量形式描述；(3)通过两待匹配图像附带特征向量、特征点两两比较找出相互匹配的若干对特征点，也就建立了图像间的对应关系。

3-1)特征检测：

本特征检测是基于尺度空间的理论。图像I中)处的点定义在尺度σ上的Hessian矩阵为：

其中， $g\left(\mathrm{\σ}\right)=\frac{1}{2{\mathrm{\π\σ}}^2}{e}^{-(x^2+y^2)/2{\mathrm{\σ}}^2}$

表示高斯滤波二阶导

与原始图像I＝(x，y)卷积；

表示高斯滤波二阶导

与原始图像I＝(x，y)卷积；

表示高斯滤波二阶导

与原始图像I＝(x，y)卷积；

构建积分图像，点

处积分图像

定义为由原点和

构成的矩形区域内所有像素值的和：

在原始图像上，通过扩大方框的大小来形成尺度不同的图像金字塔，9×9的方框滤波模板值如图3所示。

对应二阶高斯滤波σ＝1.2、相应的尺度值为s＝σ＝1.2，D_xx、D_xy、D_yy分别代表方框滤波模板与图像卷积后所得的值；继续求解Hessian矩阵的行列式值：

$\frac{{\left|L_{\mathrm{xy}}\left(1.2\right)\right|}_F{\left|L_{\mathrm{xx}}\left(0.9\right)\right|}_F}{{\left|L_{\mathrm{xx}}\left(1.2\right)\right|}_F{\left|L_{\mathrm{xy}}\left(0.9\right)\right|}_F}=0.912\≈0.9---\left(1.8\right)$

其中

是矩阵的Frobenius范数；进一步求解得到Hessian矩阵的Δ表达式为

det(H)＝D_xxD_xy-(0.9D_xx)²

(1.9)

构建尺度金子塔：每阶选4层尺度图像，3阶的构建参数如图4所示。图4中，方框内的数字代表滤波模板的大小，如果原始图像远远大于该模板值大小的时候，可继续增加模板阶次(Octave)。上述9×9的滤波为初始层，对应的尺度因子为σ＝1.2，后续层通过依次放大滤波方框的大小来计算，滤波大小依次为9×9、15×15、21×21、27×27等，当尺度较大时，滤波大小的步长也该相应增加。例如：当滤波模板的大小为N×N时，对应的尺度s＝1.2×N9；利用Hessian矩阵得到极值后，接下来选取3×3×3立体邻域内进行非极大值拟制，只有当该点比上一尺度、本尺度周围、下一尺度总共26个邻域像素值都小或者都大的时候，才将该点确定为候选特征点，然后在尺度空间和图像空间中进行插值运算，得到稳定的特征点位置及所在的尺度值。

3-2)主方向确定

为确保特征点的旋转不变性，先对以特征点为中心，计算以6s(s为该特征点所在的尺度值)为半径的邻域内的所有点在x、y方向上的Haar小波(小波边长设为4s)响应，赋予与这些邻域点离特征点的距离大小成正比的高斯权重系数；然后将每60°的范围内的响应叠加以形成新的矢量，并遍历整个圆形区域，将最长矢量的方向定为该特征点的主方向，从而获得每个特征点的主方向。

3-3)特征点描述向量

确定了特征点的主方向后，以特征点为中心将坐标轴旋转至主方向，按照主方向选取边长为20s的正方形区域，将该区域划分成4×4的子区域，在每一个子区域内，计算5s×5s范围内的小波响应，设定：

d_x--相对于主方向的水平方向的Haar小波响应

d_y--相对于主方向的垂直方向的Haar小波响应

赋予上述两个响应以权值，增加其鲁棒性、减小误差。再求取每个子区域内的响应及响应绝对值之和：∑d_x，∑d_y，∑|d_x|，∑|d_y|，并求得每个子区域的四维分量的矢量V_sub＝(∑d_x，∑|d_x|，∑d_y，∑|d_y|)。因此，对每一特征点，则形成4×(4×4)＝64维的描述向量，再进行向量的归一化处理。

经过以上计算，每一特征点都包含了三个信息：特征尺度值、特征点坐标和64维邻域向量。通过所记录的特征点信息，分别进行后续图像配准工作。

4、基于最大互信息法的精确配准方法

本例中的互信息用熵来表示，表达的是一个系统的复杂性或不确定性，从而避免了因分割带来的误差，因具有精度较高、稳健性强、不需要预处理而能实现自动配准的特点。

对于概率分布函数为p(a)的随机变量集A，其熵H(A)定义如下：

H(A)＝-∑p(a)logp(a)  a∈A    (1.10)

对两个离散的随机变量A和B，它们的边缘概率分布函数分别是p(a)和p(b)，联合概率分布函数是p_AB(a，b)，则随机变量A和B的联合熵定义如下：

H(A)＝-∑p(a)logp(a)log(b)  a∈A    b∈B    (1.11)

如果H(A/B)表示已知系统B是A的条件熵，那么H(A)与H(A/B)的差值，就代表了在系统B中包含的A的信息，即互信息。因此两个系统间的互信息定义为：

I(A，B)＝H(A)+H(B)-H(A，B)

＝H(A)-H(A/B)

＝H(B)-H(B/A)    (1.12)

在图像配准中，虽然两幅图像来源于不同的成像设备，但是它们基于共同的图像原信息，所以当两幅图像的空间位置达到完全一致时，其中一幅图像A表达的关于另一幅图像B的信息，也就是对应像素灰度的互信息应为最大。用联合概率分布和完全独立时的概率分布间的广义距离来估计互信息：

$I(A,B)=\mathrm{\Σp}(a,b)\log \frac{p(a,b)}{p\left(a\right)p\left(b\right)}---\left(1.13\right)$

对于待配准的两幅图像，可以认为它们是关于图像灰度的两个随机变量集参考图像和浮动图像，a和b是两幅图像中相关的像素灰度值，a和b通过图像之间的坐标变换相联系，对于离散的数字图像，联合概率分布p_AB(a，b)可以用归一化的联合直方图表示：

$p_{\mathrm{AB}}(i,j)=\frac{h(i,j)}{{\mathrm{\Σ}}_{i,j}h(i,j)}---\left(1.14\right)$

边缘概率分布p_A(b)表示为：

$p_A\left(i\right)=\underset{j}{\mathrm{\Σ}}{p}_{\mathrm{AB}}(i,j)---\left(1.15\right)$

边缘概率分布p_B(b)表示为：

$p_B\left(j\right)=\underset{i}{\mathrm{\Σ}}{p}_{\mathrm{AB}}(i,j)---\left(1.16\right)$

则有：

$I(A,B)=\underset{i,j}{\mathrm{\Σ}}{p}_{\mathrm{AB}}(i,j)\log \frac{p_{\mathrm{AB}}(i,j)}{p_A\left(i\right)\·p_B\left(j\right)}---\left(1.17\right)$

接下来寻找一个变换使得一幅X射线图像经过此变换后和激光图像的互信息最大，可以采用仿射变换，即在二维空间中寻找二个方向上的平移值及旋转角度。

所述仿射变换的原理是：若变换S：Rⁿ→Rⁿ，S(x)＝T(x)+a，T是非奇异线性变换，a∈Rⁿ，则变换S称为仿射变换。二维欧氏空间上的仿射变换为：

S(X)＝T(X)+A

(2.2)

式中：

$X=\left[\begin{array}{c}x\\ y\end{array}\right]$ $T=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$ $A=\left[\begin{array}{c}e\\ f\end{array}\right]$ a、b、c、d、e、f皆为实数。

仿射变换具有平行线转换成平行线和有限点映射到有限点的一般特性。平移、旋转等是二维仿射变换的特例。

平移：

设坐标点(x，y)经平移Δx和Δy后的坐标为(x′，y′)，仿射变换形式为：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]$

(2.3)

式中 $t=\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]$ 是平移矩阵；

旋转

坐标点(x，y)经旋转θ角后的坐标为(x′，y′)，反射变换形式为：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]$

(2.4)

式中 $R=\left[\begin{array}{cc}\cos \mathrm{\θ}& \sin \mathrm{\θ}\\ -\sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]$ 是旋转矩阵。

本方法中，图像变换之间只有旋转和平移，激光图像向X射线图像的仿射变换可用式(2.5)表示：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=\left[\begin{array}{cc}{a}_{11}& 0\\ 0& a_{22}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]+\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]$

(2.5)

其中(x，y)和(x′，y′)分别代表激光图像和X射线图像相匹配的像素点。只要给定两组(x，y)和(x′，y′)就可以求得仿射变换的四个参数a₁₁，a₂₂，Δx，Δy，那么就可以实现激光图像和X射线图像的配准。

5、判断X射线图像与激光图像配准

为判断X射线图像与激光图像配准，本方法以配准的度量函数即相似性测度，来判断两幅图像是否完全对齐，衡量经过空间变换后激光图像与X射线图像在是否实现了空间一致性。

由于X射线图像与激光图像属于异源图像配准，两图像其灰度属性存在非常大的差异，为此，本方法采用统计型相似性测度来检测图像配准的准确性，本方法利用图像的灰度特性，对两幅图像重叠区域内根据像素灰度值直接计算相似性测度函数，而不用事先选取控制点或提取图像特征，使得肉图像配准相似度测度可判。

Claims (5)

1.一种X射线图像和激光图像的异源图像配准方法，其特征是用X射线对物体照射，采集物体的X射线图像，该X射线图像是X射线透射物体得到的图像；用激光对同一物体进行照射，采集物体的激光图像，该激光图像是激光照射到物体表面反射得到的图像；

所述配准方法的步骤包括：

1)对采集到的X射线图像和激光图像进行预处理，来抑制图像噪声点的干扰；

2)对激光图像边缘进行修正；

3)分别提取X射线图像和激光图像的特征，并进行粗匹配：

3-1)通过Hessian矩阵求取两类图像的特征点候选集合，提取合适的特征点；

3-2)对特征点所在局部图像的特征用向量形式描述，得到特征向量；

3-3)对两图像的特征向量和特征点分别进行两两比较，找出相互匹配的若干对特征点，建立两图像间的对应关系；

4)X射线图像和激光图像的精确配准：

用互信息表示的相似性测度对两图像进行精确配准，通过对完整的激光图像寻找一个仿射变换即旋转参数和平移参数，该仿射变换对应的X射线图和激光图像的互信息最大，此时X射线图像和完整的激光图像实现精确配准。

2.根据权利要求1所述的X射线图像和激光图像的异源图像配准方法，其特征是所述步骤1)的预处理方法是，制作平整的正方形标定块，通过中值滤波对标定块X射线和标定块激光图像进行噪声滤波预处理。

3.根据权利要求2所述的X射线图像和激光图像的异源图像配准方法，其特征是所述步骤2)中，先对预处理后的激光图像进行二值化；再用数学形态学腐蚀运算，消除标定块激光图像边缘的毛刺；然后用膨胀运算，使图像边缘空洞的边界向外扩张，经过膨胀操作，把这些洞将被补上。

4.根据权利要求3所述的X射线图像和激光图像的异源图像配准方法，其特征是所述步骤3)中，

先进行特征检测：

设图像I中，点

)定义在尺度σ上的Hessian矩阵为：

其中， $g\left(\mathrm{\σ}\right)=\frac{1}{2{\mathrm{\π\σ}}^2}{e}^{-(x^2+y^2)/2{\mathrm{\σ}}^2};$

表示高斯滤波二阶导

与原始图像I＝(x，y)卷积；

表示高斯滤波二阶导与原始图像I＝(x，y)卷积；

表示高斯滤波二阶导

与原始图像I＝(x，y)卷积；

构建积分图像，点

处积分图像定义为由原点和

构成的矩形区域内所有像素值的和：

所述二阶高斯滤波σ是通过在原始图像I上，扩大方框滤波模板的方框的大小来形成尺度不同的图像金字塔，得到不同的方框滤波模板值；

初始层的方框滤波模板值对应的尺度因子为σ；对应二阶高斯滤波σ、相应的尺度值为s＝σ；

继续求解Hessian矩阵的行列式值：

$\frac{{\left|L_{\mathrm{xy}}\left(1.2\right)\right|}_F{\left|L_{\mathrm{xx}}\left(0.9\right)\right|}_F}{{\left|L_{\mathrm{xx}}\left(1.2\right)\right|}_F{\left|L_{\mathrm{xy}}\left(0.9\right)\right|}_F}=0.912\≈0.9---\left(1.8\right)$

设是矩阵的Frobenius范数，进一步求解得到Hessian矩阵的Δ表达式为

det(H)＝D_xxD_y-(0.9D_xx)²    (1.9)；

再进行主方向确定：

设特征点所在的尺度值为s；

先以特征点为中心，计算以6s为半径的邻域内的所有点在x、y方向上的Haar小波响应，赋予与邻域内的点与特征点的距离大小成正比的高斯权重系数；

然后将每60°的范围内的响应叠加，来形成新的矢量，并遍历整个邻域，将最长矢量的方向定为该特征点的主方向；

最后进行特征点描述向量：

以特征点为中心将坐标轴旋转至主方向，按照主方向选取边长为20s的正方形区域，将该正方形区域划分成4×4的子区域；

在每一个子区域内，先计算5s×5s范围内的小波响应：

设d_x为相对于主方向的水平方向的Haar小波响应、d_y为相对于主方向的垂直方向的Haar小波响应，赋予这两个响应以权值；

再求每个子区域内的响应及响应绝对值之和：∑d_x，∑d_y，∑|d_x|，∑|d_y|，并求得每个子区域的四维分量的矢量V_sub＝(∑d_x，∑|d_x|，∑d_y，∑|d_y|)；对每一特征点，则形成的4×(4×4)＝64维的描述向量，然后进行向量的归一化处理；

经过以上计算，每一特征点都包含了三个信息：特征尺度值、特征点坐标和64维邻域向量，通过所记录的特征点信息，分别进行后续图像配准工作。

5.根据权利要求1所述的X射线图像和激光图像的异源图像配准方法，其特征是所述步骤4)中，互信息用熵来表示：

对于概率分布函数为p(a)的随机变量集A，其熵H(A)定义如下：

H(A)＝-∑p(a)logp(a)  a∈A    (1.10)

对两个图像，它们的离散的随机变量集分别是A和B，它们的边缘概率分布函数分别是p(a)和p(b)，联合概率分布函数是p_AB(a，b)，则随机变量A和B的联合熵定义如下：

H(A)＝-∑p(a)logp(a)log(b)  a∈Ab∈B    (1.11)

如果H(A/B)表示已知B是A的条件熵，那么H(A)与H(A/B)的差值，就代表了在B中包含的A的信息，即互信息；因此两个随机变量集间的互信息定义为：

I(A，B)＝H(A)+H(B)-H(A，B)

＝H(A)-H(A/B)

＝H(B)-H(B/A)    (1.12)

在图像配准中，当两幅图像的空间位置达到完全一致时，其中一幅图像A表达的关于另一幅图像B的信息，也就是对应像素灰度的互信息应为最大；用联合概率分布和完全独立时的概率分布间的广义距离来估计互信息：

$I(A,B)=\mathrm{\Σp}(a,b)\log \frac{p(a,b)}{p\left(a\right)p\left(b\right)}---\left(1.13\right)$

对于待配准的两幅图像，可以认为它们是关于图像灰度的两个随机变量集参考图像和浮动图像，a和b是两幅图像中相关的像素灰度值，a和b通过图像之间的坐标变换相联系，对于离散的数字图像，联合概率分布p_AB(a，b)可以用归一化的联合直方图表示：

$p_{\mathrm{AB}}(i,j)=\frac{h(i,j)}{{\mathrm{\Σ}}_{i,j}h(i,j)}---\left(1.14\right)$

边缘概率分布p_A(b)表示为：

$p_A\left(i\right)=\underset{j}{\mathrm{\Σ}}{p}_{\mathrm{AB}}(i,j)---\left(1.15\right)$

边缘概率分布p_B(b)表示为：

$p_B\left(j\right)=\underset{i}{\mathrm{\Σ}}{p}_{\mathrm{AB}}(i,j)---\left(1.16\right)$

则有：

$I(A,B)=\underset{i,j}{\mathrm{\Σ}}{p}_{\mathrm{AB}}(i,j)\log \frac{p_{\mathrm{AB}}(i,j)}{p_A\left(i\right)\·p_B\left(j\right)}---\left(1.17\right);$

接下来寻找一个合适的仿射变换，使得一幅X射线图像经过此仿射变换后和激光图像的互信息最大。

Images