一种X射线影像的自适应拼接方法及其拼接系统

摘要

本发明公开了一种X射线影像的自适应拼接方法，具体包括如下步骤：步骤1、获取连续两幅X射线影像的横向图像带；步骤2、分别计算两条横向图像带的方差并确定模板图像和移动图像；步骤3、配准模板图像和移动图像；步骤4、采用重叠区非线性融合法对图像进行融合，得到拼接后的新图像。本发明通过比较分别从两幅X射线影像图的重叠区部分选取的横向图像带的灰度值的方差，自适应选择比较清晰且所含信息较多的横向图像带作为模板和另一幅影像图进行配准，进而采用重叠区域非线性过度法将输入的影像图进行融合而实现全尺度拼接。

说明书

技术领域

本发明属于数字图像检索技术领域，具体涉及一种X射线影像的自适应拼接 方法及其拼接系统。本发明用于对X射线影像设备采集的不同区域的X射线影 像进行自动拼接。

背景技术

数字化时代的来临使得X射线影像的数字化成为趋势，X射线影像极大地 方便了医生对患者的病灶检查以及定位，对病情的预断和预后起到了极大的作 用。目前，数字图像处理算法如图像增强、边缘检测、轮廓提取、图像配准拼 接、图像编码等已广泛应用，给医生的工作带来了极大的方便的同时减少了误 诊率。现有的X射线影像设备（如CR、DR等）一次拍摄的图像尺寸有限（约 为17×17英寸），要获得全脊柱、全下肢或全身完整图像至少需要拍摄两幅有 一定重合区域的图像。拍摄图像时，为了不丢失信息影响医生工作，通常规定 两幅连续图像的重叠区的高度大于250个像素，为了更全面、更直观的了解患 者病灶及其周围部位情况，需要将这些图像进行无缝拼接得到超大视角的图像。

传统的数字图像拼接算法主要包括图像配准和图像融合两个步骤。首先， 图像配准算法主要分为基于特征的配准算法和基于灰度信息的配准算法，由于X 射线影像的特征点一般较少，基于特征的图像配准方法往往精度达不到医学图 像配准的精度要求，因此，目前医学图像配准的研究热点主要集中在基于灰度 信息的配准算法，而该算法中基于归一化互信息的配准算法因具有较好的鲁棒 性和准确性而得到了广泛的应用，但是，传统的归一化互信息法图像配准算法 在配准图像时固定地选取其中一幅图像重叠区域的部分作为模板图像与另一幅 图像进行配准，而两幅图像在重叠区域所含的图像信息是不一样的，这样导致 最终结果随模板选择的不同而配准结果也有较大的变化。其次，常见的图像融 合方法有平均值法和重叠区线性过度法，其中较为常用的方法是重叠区线性过 度法，它使用线性变换的权重来融合图像以使过渡区看起来较柔和，但是，该 方法无法充分消除两幅待配准图像边缘部分过曝、图像边缘部分过度变形等对 图像配准精度的影响。

发明内容

针对上述现有技术中存在的缺陷或不足，本发明的目的在于，提供一种应 用于X射线影像的自适应拼接方法及其拼接系统，本发明通过比较分别从两幅 X射线影像图的重叠区部分选取的横向图像带的灰度值的方差，自适应选择比 较清晰且所含信息较多的横向图像带作为模板和另一幅影像图进行配准，进而 采用重叠区域非线性过度法将输入的影像图进行融合而实现全尺度拼接。

为了达到上述目的，本发明采用如下技术方案予以解决：

一种X射线影像的自适应拼接方法，具体包括如下步骤：

步骤1、获取连续两幅X射线影像的横向图像带；

步骤2、分别计算两条横向图像带的方差并确定模板图像和移动图像；

步骤3、配准模板图像和移动图像；

步骤4、采用重叠区非线性融合法对图像进行融合，得到拼接后的新图像。

进一步的，所述步骤1的具体操作如下：

顺序读入两幅待拼接的X射线影像图像：第一幅影像图I₁和第二幅影像图 I₂；在第一幅影像图I₁上的与第二幅影像图I₂的重叠区域截取横向图像带I_u，在 第二幅影像图I₂上的与第一幅影像图I₁的重叠区域截取横向图像带I_d，保证横向 图像带I_u与横向图像带I_d的宽度相等且均等于待拼接X射线影像图像的宽度， 同时保证横向图像带I_u与横向图像带I_d的高度相等且均不大于250个像素。

进一步的，所述步骤2的具体操作如下：

1）按照公式6、公式7分别计算横向图像带I_u的方差和横向图像带I_d的 方差：

$\left(\begin{array}{c}{\mu }_u=\frac{1}{W_t\times H_t}{\Sigma }_{(x,y)\in I_u}{f}_u(x,y)\\ {\sigma }_u^2=\frac{1}{W_t\times H_t-1}{\Sigma }_{(x,y)\in I_u}{[f_u(x,y)-{\mu }_u]}^2\end{array}\right)---\left(6\right)$

$\left(\begin{array}{c}{\mu }_d=\frac{1}{W_t\times H_t}{\Sigma }_{(x,y)\in I_d}{f}_d(x,y)\\ {\sigma }_d^2=\frac{1}{W_t\times H_t-1}{\Sigma }_{(x,y)\in I_d}{[f_d(x,y)-{\mu }_d]}^2\end{array}\right)---\left(7\right)$

其中，f_u(x,y)为横向图像带I_u在坐标位置(x,y)处的灰度值；f_d(x,y)为横向图 像带I_d在坐标位置(x,y)处的灰度值；W_t为横向图像带I_u和横向图像带I_d的宽度； H_t为横向图像带I_u和横向图像带I_d的高度；

2）若大于选取影像图I₁为模板图像，影像图I₂为移动图像；否则， 选取横向图像带I_d为模板图像，影像图I₁为移动图像。

进一步的，所述步骤3的具体操作如下：

使模板图像在移动图像上逐像素平移，每平移一个像素，则在移动图像上 取与模板图像对应的图像部分，并计算移动图像与模板图像之间的归一化互信 息；直到移动图像上所有像素被遍历；比较得到归一化互信息最大时模板图像 所在的位置，记录在该位置时模板图像在垂直方向上距初始遍历位置的平移量 Δ，通过式8得到影像图I₁和影像图I₂的重叠区的高度D：

D＝Δ+H_t              （8）

式中，H_t为横向图像带I_u和横向图像带I_d的高度；

进一步的，所述步骤4的具体操作如下：

设影像图I₁的高度为H₁，影像图I₂的高度为H₂，创建一幅宽度为W_d、高度 为H₁+H₂-D的新图像，对于重叠区上方即高度从1到H₁-D的部分，将影像图I₁中的对应部分复制到新图像的对应位置；对于重叠区下方即高度从H₁+1到 H₁+H₂-D部分，将影像图I₂中的对应部分复制到新图像的对应位置；

对于重叠区内坐标位置为(x,y)的点，找到该点对应于影像图I₁的像素值 f₁(x,y)以及该点对应于影像图I₂的像素值f₂(x,y)，利用式9计算得到新图像的重 叠区内坐标位置为(x,y)的点的像素值f(x,y)，最终得到新图像。

$\left(\begin{array}{c}\omega =2^{-\frac{2}{3}}{(\frac{y}{D}-\frac{1}{2})}^{\frac{1}{3}}+\frac{1}{2}\\ f(x,y)=(1-\omega )\times f_1(x,y)+\omega \times f_2(x,y)\end{array}\right)---\left(9\right).$

应用上述X射线影像的自适应拼接方法的拼接系统，具体包括如下依次相 连接的模块：

图像输入及横向图像带获取模块：用于获取连续两幅X射线影像的横向图 像带；

模板图像和移动图像确定模块：用于分别计算两条横向图像带的方差并确 定模板图像和移动图像；

模板图像和移动图像配准模块：用于配准模板图像和移动图像；

图像融合模块：用于采用重叠区非线性融合法对图像进行融合，得到拼接 后的新图像。

进一步的，所述图像输入及横向图像带获取模块具体用于实现以下功能：

顺序读入两幅待拼接的X射线影像图像：第一幅影像图I₁和第二幅影像图 I₂；在第一幅影像图I₁上的与第二幅影像图I₂的重叠区域截取横向图像带I_u，在 第二幅影像图I₂上的与第一幅影像图I₁的重叠区域截取横向图像带I_d，保证横向 图像带I_u与横向图像带I_d的宽度相等且均等于待拼接X射线影像图像的宽度， 同时保证横向图像带I_u与横向图像带I_d的高度相等且均不大于250个像素。

进一步的，所述模板图像和移动图像确定模块具体用于实现以下功能：

1）按照公式6、公式7分别计算横向图像带I_u的方差和横向图像带I_d的 方差

$\left(\begin{array}{c}{\mu }_u=\frac{1}{W_t\times H_t}{\Sigma }_{(x,y)\in I_u}{f}_u(x,y)\\ {\sigma }_u^2=\frac{1}{W_t\times H_t-1}{\Sigma }_{(x,y)\in I_u}{[f_u(x,y)-{\mu }_u]}^2\end{array}\right)---\left(6\right)$

$\left(\begin{array}{c}{\mu }_d=\frac{1}{W_t\times H_t}{\Sigma }_{(x,y)\in I_d}{f}_d(x,y)\\ {\sigma }_d^2=\frac{1}{W_t\times H_t-1}{\Sigma }_{(x,y)\in I_d}{[f_d(x,y)-{\mu }_d]}^2\end{array}\right)---\left(7\right)$

其中，f_u(x,y)为横向图像带I_u在坐标位置(x,y)处的灰度值；f_d(x,y)为横向图 像带I_d在坐标位置(x,y)处的灰度值；W_t为横向图像带I_u和横向图像带I_d的宽度； H_t为横向图像带I_u和横向图像带I_d的高度；

2）若大于选取影像图I₁为模板图像，影像图I₂为移动图像；否则， 选取横向图像带I_d为模板图像，影像图I₁为移动图像。

进一步的，所述模板图像和移动图像配准模块具体用于实现以下功能：

使模板图像在移动图像上逐像素平移，每平移一个像素，则在移动图像上 取与模板图像对应的图像部分，并计算移动图像与模板图像之间的归一化互信 息；直到移动图像上所有像素被遍历；比较得到归一化互信息最大时模板图像 所在的位置，记录在该位置时模板图像在垂直方向上距初始遍历位置的平移量 Δ，通过式8得到影像图I₁和影像图I₂的重叠区的高度D：

D＝Δ+H_t               （8）

式中，H_t为横向图像带I_u和横向图像带I_d的高度；

进一步的，所述图像融合模块具体用于实现以下功能：

设影像图I₁的高度为H₁，影像图I₂的高度为H₂，创建一幅宽度为W_d、高度为 H₁+H₂-D的新图像，对于重叠区上方即高度从1到H₁-D的部分，将影像图I₁中 的对应部分复制到新图像的对应位置；对于重叠区下方即高度从H₁+1到 H₁+H₂-D部分，将影像图I₂中的对应部分复制到新图像的对应位置；

对于重叠区内坐标位置为(x,y)的点，找到该点对应于影像图I₁的像素值 f₁(x,y)以及该点对应于影像图I₂的像素值f₂(x,y)，利用式9计算得到新图像的重 叠区内坐标位置为(x,y)的点的像素值f(x,y)，最终得到新图像；

$\left(\begin{array}{c}\omega =2^{-\frac{2}{3}}{(\frac{y}{D}-\frac{1}{2})}^{\frac{1}{3}}+\frac{1}{2}\\ f(x,y)=(1-\omega )\times f_1(x,y)+\omega \times f_2(x,y)\end{array}\right)---\left(9\right).$

与传统的X射线拼接方法及系统相比较，本发明的优点如下：

第一，本发明分别计算两幅图像部分重叠区域所含信息，选择较清晰且所 含信息较多的部分图像重叠区作为模板与另一幅影像图像进行配准，能够自适 应地达到更精确的配准结果，避免了传统的归一化互信息法图像配准算法的最 终结果对模板选择的过度依赖性。

第二，本发明在进行图像融合的过程中，选用重叠区非线性过渡融合法， 在将重叠区域进行融合的时权重的变化是非线性的，使得图像边缘部分在图像 融合时所占的权重更少，从而使得重叠区域更加平滑，能够更好地抑制图像边 缘区域过曝已经拍摄过程中抖动等影响，避免了重叠区线性过度法无法充分消 除图像边缘部分过曝和过度变形对配准精度的影响。

附图说明

图1是本发明的X射线影像的自适应拼接算法的流程图。

图2是本发明的X射线影像的自适应拼接系统的结构示意图。

图3是本发明的实施例1的效果对比图。其中，（a）中由上至下依次为待 拼接的三幅连续的X射线影像图像；（b）为采用本发明的方法进行拼接之后得 到的X射线影像图像。

图4是本发明的实施例1的步骤3计算得到的重叠高度D。图中重叠区域 均用水平直线标出。

图5是本发明的实施例1的步骤3中，模板图像和移动图像配准过程的示 意图。其中，（a）、（b）、（c）分别为将移动图像在模板图像上由下到上逐像素移 动并定位归一化互信息最大的位置的过程。（a）、（b）、（c）中左边的图中均包括 三条水平线，该三条水平线由上至下分别表示模板图像匹配过程的最终位置、 归一化互信息最大的位置、模板图像匹配的初始位置。

图6是直接进行图像拼接而未进行图像融合和按照本发明的步骤4进行了 重叠区非线性融合前后的效果对比。其中，（a）表示两幅连续图像直接拼接的 效果；（b）为按照本发明的步骤4进行了重叠区非线性融合前后的效果。

图7是本发明的实施例2的人体侧位X射线影像图像自适应拼接效果图。 其中，（a）中由上至下依次为待拼接的三幅连续的X射线影像图像；（b）为采 用本发明的方法进行拼接之后得到的X射线影像图像。

具体实施方式

首先，对本发明中涉及的相关技术术语介绍如下：

1、图像灰度的方差：即衡量图像所有像素点灰度值波动大小的一种量值， 方差越大，表明该图像的灰度值范围越宽，灰度值变化越大，即该图像所含有 的信息也就越多。在对X射线影像图进行配准时，可供选择的模板是加载的两 幅影像图像中上边图像的部分重叠区或下边图像的部分重叠区。选择哪个重叠 区取决于哪一幅影像图足够清晰，所含信息较多，即图像灰度值的方差越大。 当用较好的边缘进行配准时，能够达到更好的配准精度。

记f(x,y)为灰度图像I在坐标位置为(x,y)的灰度值，W为灰度图像的宽度， H为灰度图像的高度，μ为灰度图的均值，则灰度图像I的灰度的方差σ²的计算 公式如下：

$\left(\begin{array}{c}\mu =\frac{1}{W\times H}{\Sigma }_{(x,y)\in I}f(x,y)\\ {\sigma }^2=\frac{1}{W\times H-1}{\Sigma }_{(x,y)\in I}{[f(x,y)-\mu ]}^2\end{array}\right)---\left(1\right)$

2、归一化互信息（Normalized Mutual Information）：用于描述两个系统间的 统计相关性，或者是一个系统中所包含另一个系统中信息的多少，它可用熵 （Entropy）来描述。熵表达的是一个系统的复杂性或者是不确定性，当两幅图 像的空间位置达到完全一致时，其中一幅图像表达的关于另一幅图像的信息， 也就是对应像素灰度的互信息应为最大。对于离散的数字图像，两幅灰度图的 互信息可用归一化的联合直方图表示，对于两幅宽度均为W'，高度均为H'的同 样大小的图像A和B，记M_A(i)为图像A中灰度值为i的像素点的个数，M_B(j)为 图像B中灰度值为j的像素点的个数，M_AB(i,j)为图像A中灰度值为i同时在图像 B中灰度值为j的像素点个数，则图像A的灰度直方图h_A(i)和图像B的灰度直方 图h_B(j)以及它们的联合灰度直方图表示为：

$\left(\begin{array}{c}{h}_A\left(i\right)=\frac{M_A\left(i\right)}{W^{\prime }\times H^{\prime }}\\ h_B\left(j\right)=\frac{M_B\left(j\right)}{W^{\prime }\times H^{\prime }}\\ h_{\mathrm{AB}}(i,j)=\frac{M_{\mathrm{AB}}(i,j)}{W^{\prime }\times H^{\prime }}\end{array}\right)---\left(2\right)$

利用式3计算图像A和图像B的边界熵H(A)和H(B)以及两幅图像的联合熵 H(A,B)：

$\left(\begin{array}{c}H\left(A\right)=-\underset{i}{\Sigma }{h}_A\left(i\right)\log \left(h_A\left(i\right)\right)\\ H\left(B\right)=-\underset{j}{\Sigma }{h}_B\left(j\right)\log \left(h_B\left(j\right)\right)\\ H(A,B)=-\underset{i}{\Sigma }\underset{j}{\Sigma }{h}_{\mathrm{AB}}(i,j)\log \left(h_{\mathrm{AB}}(i,j)\right)\end{array}\right)---\left(3\right)$

利用式4计算图像A和图像B之间的归一化互信息：

$I(A,B)=\frac{H\left(A\right)+H\left(B\right)}{H(A,B)}---\left(4\right)$

3、重叠区非线性融合：为了使融合之后两幅图像重合区更加光滑而采用的 区别于线性权重的图像融合方法。在对连续的两幅影像图进行配准之后需要将 它们之间的重叠区进行融合，传统的图像融合在计算融合之后的图像的灰度值 时，认为距离被拼接图像的中心越远，相应的权重越小，但是权重是线性过渡 的，如果采用非线性的过渡权重，能够更好地利用原始图像中远离图像边缘部 分的灰度信息。在对两幅连续影像图I₁和I₂利用归一化互信息方法进行配准之 后，可得到它们之间重叠区域的高度D，对于重叠区内坐标位置为(x,y)的点， 找到该点对应于影像图I₁的像素值f₁(x,y)以及该点对应于影像图I₂的像素值 f₂(x,y)，利用式5计算非线性权重ω，得到融合图像重叠区内坐标位置为(x,y)的 点的像素值f(x,y)：

$\left(\begin{array}{c}\omega =2^{-\frac{2}{3}}{(\frac{y}{D}-\frac{1}{2})}^{\frac{1}{3}}+\frac{1}{2}\\ f(x,y)=(1-\omega )\times f_1(x,y)+\omega \times f_2(x,y)\end{array}\right)---\left(5\right)$

如图1所示，本发明的X射线影像的自适应拼接方法，具体包括如下步骤：

步骤1、获取连续两幅X射线影像的横向图像带。具体操作如下：

顺序读入两幅待拼接的X射线影像图像：第一幅影像图I₁和第二幅影像图I₂（本发明中认为待拼接的连续的X射线影像图像的宽度相等）。在第一幅影像图 I₁上的与第二幅影像图I₂的重叠区域截取横向图像带I_u，在第二幅影像图I₂上的 与第一幅影像图I₁的重叠区域截取横向图像带I_d，保证横向图像带I_u与横向图像 带I_d的宽度相等且均等于待拼接X射线影像图像的宽度，同时保证横向图像带I_u与横向图像带I_d的高度相等。

步骤2、分别计算两条横向图像带的方差并确定模板图像。具体操作如下：

1）按照公式6、公式7分别计算横向图像带I_u的方差和横向图像带I_d的 方差

$\left(\begin{array}{c}{\mu }_u=\frac{1}{W_t\times H_t}{\Sigma }_{(x,y)\in I_u}{f}_u(x,y)\\ {\sigma }_u^2=\frac{1}{W_t\times H_t-1}{\Sigma }_{(x,y)\in I_u}{[f_u(x,y)-{\mu }_u]}^2\end{array}\right)---\left(6\right)$

$\left(\begin{array}{c}{\mu }_d=\frac{1}{W_t\times H_t}{\Sigma }_{(x,y)\in I_d}{f}_d(x,y)\\ {\sigma }_d^2=\frac{1}{W_t\times H_t-1}{\Sigma }_{(x,y)\in I_d}{[f_d(x,y)-{\mu }_d]}^2\end{array}\right)---\left(7\right)$

其中，f_u(x,y)为横向图像带I_u在坐标位置(x,y)处的灰度值；f_d(x,y)为横向图 像带I_d在坐标位置(x,y)处的灰度值；W_t为横向图像带I_u和横向图像带I_d的宽度； H_t为横向图像带I_u和横向图像带I_d的高度；

2）若大于选取影像图I₁为模板图像，影像图I₂为移动图像；否则， 选取横向图像带I_d为模板图像，影像图I₁为移动图像；

步骤3、配准模板图像和移动图像。具体操作如下：

使模板图像在移动图像上逐像素平移，每平移一个像素，则在移动图像上 取与模板图像对应的图像部分，并计算移动图像与模板图像之间的归一化互信 息；直到移动图像上所有像素被遍历；比较得到归一化互信息最大时模板图像 所在的位置，记录在该位置时模板图像在垂直方向上距初始遍历位置的平移量 Δ，通过式8得到影像图I₁和影像图I₂的重叠区的高度D：

D＝Δ+H_t             （8）

式中，H_t为横向图像带I_u和横向图像带I_d的高度；

步骤4、采用重叠区非线性融合法对图像进行融合。具体操作如下：

设影像图I₁的高度为H₁，影像图I₂的高度为H₂，创建一幅宽度为W_d、高度为 H₁+H₂-D的新图像，对于重叠区上方即高度从1到H₁-D的部分，将影像图I₁中 的对应部分复制到新图像的对应位置；对于重叠区下方即高度从H₁+1到 H₁+H₂-D部分，将影像图I₂中的对应部分复制到新图像的对应位置；

对于重叠区内坐标位置为(x,y)的点，找到该点对应于影像图I₁的像素值 f₁(x,y)以及该点对应于影像图I₂的像素值f₂(x,y)，利用式9计算得到新图像的重 叠区内坐标位置为(x,y)的点的像素值f(x,y)，最终得到新图像。

$\left(\begin{array}{c}\omega =2^{-\frac{2}{3}}{(\frac{y}{D}-\frac{1}{2})}^{\frac{1}{3}}+\frac{1}{2}\\ f(x,y)=(1-\omega )\times f_1(x,y)+\omega \times f_2(x,y)\end{array}\right)---\left(9\right)$

为了说明本发明的方法的有效性，发明人进行了X射线影像图像拼接测试。 采用的操作系统为Win764、CPU为Intel(R)Core(TM)i5-2300CPU2.80GHz。

实施例1：

本实施例对连续的三幅人体正面的X射线图像进行配准融合，如图3（a） 所示。分别选取第一副X射线图像的下边缘重叠区的一条横向图像带和第二幅 X图像的上边缘重叠区的一条横向图像带，两个横向图像带的高度相等且均为 250个像素；然后比较该两条图像带灰度值的方差，选取方差较大的图像带作为 模板图像和另一幅图像进行比较。在该实施例中，选取第二副X射线图像的上 边缘重叠区作为模板图像，选取第一幅图像为移动图像，将模板图像与第一幅 图像利用步骤3的归一化互信息比较进行配准，最后采用重叠区域非线性融合 法对配准图像进行融合。第二幅与第三幅图像配准过程同上，配准结果如图3 （b）所示。由图3可以看到，使用步骤4进行了重叠区域融合之后，待拼接图 像边缘处的过度曝光的情况被改善，得到图3（b）所示的过渡柔和的拼接后的 图像，可知本发明的方法有效地抑制了边缘部分过度曝光或者拍摄抖动等影响。

实施例2：

本实施例与实施例相同，区别仅在于本实施例对人体侧位的X射线影像图 像自适应拼接，拼接前图像与拼接后的图像的对比效果见图7。

发明人通过对100组连续的X射线影像图像进行拼接测试，得到的结果稳 定可靠，平均处理时间为279.6ms，可以看出，本发明的方法不仅能够准确柔和 地拼接X射线影像，且处理速度快，运行效率高。