CN101224114A - 一种基于尺度空间分解的x射线图像超动态范围重建方法 - Google Patents

Abstract

本发明公开了一种基于尺度空间分解与重构的X射线图像超动态范围重建方法，该方法的步骤包括：拍摄同一成像对象在不同成像条件下的两幅常规X线图像，这两幅图像分别显示成像对象不同感兴趣区域的信息。首先对这两幅图像进行基于尺度空间的分解，然后通过一定的重构规则对两幅图像的分解分量进行重构，再此基础上重建一幅动态范围扩展的图像。该方法能有效扩展小动态范围X线成像设备的成像动态范围，具有很大的实际应用价值。

Description

一种基于尺度空间分解的X射线图像超动态范围重建方法

技术领域

本发明属于图像处理领域，涉及X射线图像的处理，尤其涉及X射线图像动态范围的改善，具体涉及一种基于尺度空间分解与重构的X射线图像处理方法。

背景技术

X射线成像技术在医疗诊断和工业领域有着广泛应用。数字化X线影像技术是计算机数字图像处理技术与X射线放射技术相结合而形成的一种先进的X线成像技术。对医学诊断中常用的胸部X线成像而言，由于胸部组织结构的密度和厚度存在着较大的差异，传统X线平片往往难以同时显示不同部位组织结构的病变，而数字化X线摄影的最大优点就在于它只需一次曝光，通过调谐和空间频率的处理，可以获得多幅图像，分别显示不同的组织结构。从20世纪80年代起，数字化X射线成像技术CR及DR逐渐取代了传统的胶片成像技术。CR与DR的共同点都是将X线影像信息转化为数字影像信息。现在普遍应用的DR主要是采用平板探测器(FPD)对X线产生的图像信号进行扫描和直接读出，成像原理是先将X线信号转变为可见光，通过光电2极管组成的藻膜层(TFT)进行聚集，再由专门的读出电路直接读出送入计算机系统进行处理。目前FPD每个成像单元的面积可达127μm×127μm，可使成像具有很宽的动态范围，但这种DR系统价格昂贵，极大地限制了其普及和推广。由于碘化铯钠(CsI:Na)具有将X光转化为可见光的能力，用碘化铯钠加光学镜头耦合CCD从理论上讲可以实现DR探测系统的功能，这种DR系统被称为CCD-DR。CCD每一个光敏元件称为一个像素，它是成像的基本单元，其面积可达到15μm×15μm。由于CCD成像单元面积大小的限制，CCD-DR设备的成像动态范围很小。因此，如何使CCD-DR系统获得较宽成像动态范围的X射线图像是非常有意义的问题。

目前，国内外还未针对X射线图像提出扩展其动态范围的方法。

本发明提出一种基于尺度空间分解的X射线图像超动态范围重建方法。该方法针对成像动态范围小的X线设备，通过在两个不同球管电压下拍摄的同一成像对象的两幅X射线图像(这里称为原始图像)重建该成像对象具有较大动态范围的图像。首先利用尺度空间分解的方法将两幅原始图像进行分解，再根据一定的重构规则将两幅原始图像分解后的分解分量进行重构，在此基础上重建出一幅具有宽动态范围的图像。

该方法的优点是针对同一成像对象在不同成像条件下获得的两幅图像，用基于尺度空间分解的方法，通过系数重构获得较宽动态范围的图像，算法简单易行。所得的重建图像有较宽的动态范围，较高的空间分辨率和密度分辨率。利用这种图像超动态范围重建方法，可以扩展CCD-DR系统的成像动态范围，从而使这种廉价的DR设备得以广泛应用。

发明内容

本发明的目的在于，提出一种基于尺度空间分解的X射线图像超动态范围重建方法，该方法从图像多尺度分解的理论出发，针对动态范围小的X线成像设备构造。

为了实现上述目的，本发明采取如下技术解决方案：

一种基于尺度空间分解的X射线图像超动态范围重建方法，其特征在于，该方法具体的步骤是：

A、对同一成像对象放置在两个不同球管电压下，该两个球管电压分别选择成像对象的两个不同感兴趣区域，在两个不同球管电压下，以合适的曝光条件分别曝光，得到两幅图像，该两幅图像分别显示成像对象的不同区域信息；

B、将拍摄到的两幅图像分别进行基于尺度空间的图像分解，将两幅图像分解为控制图像整体灰度范围的大区域信号和反应图像细节的小尺度信号；其中大区域信号分解后得到大尺度分量即低频分量，主要包含图像的背景信息，小尺度信号分解后得到小尺度分量即高频分量，反映图像的细节信息；

C、在低频分量和高频分量所对应的频率上进行重组，由重组后的分解分量重建一幅动态范围较宽的重建图像。

图像多尺度分解的方法有很多，可以参考如下文献和专利：

1、中国专利申请(申请号：200610041726.5)基于尺度空间分解与重构的X射线图像均衡显示处理方法；

2、Mary Couwenhoven，Robert Senn，David Foos，et al.“Enhancement Methodthat Provides Direct and Independent Control of Fundamental Attributes of ImageQuality for Radiographic Imagery”，Proc.Of SPIE Medical Imaging 2004，Vol.5367。

3、Pietro Perona and Jitendra Malik，“Scale-Space and Edge Detection UsingAnisoptropic Diffusion”，IEEE Transactions on Pattern Analysis and MachineIntelligence，Vol.12，NO.7，July 1990。

4、Ana Petrovic，Oscar Divorra Escoda and Pierre Vandergheynst，“Multiresolution。Segmentation of Natural Images：From Linear to NonlinearScale-Space Representations”，IEEE Transactions on Image Processing，Vol.13，NO.8。

5、Min Zhang.Xuanqin Mou.A Novel Contrast Equalization Method for ChestRadiograph.Proc.of SPIE.Vol.6144 61446R-1。

6、中国专利申请(申请号：200610041726.5)基于尺度空间分解与重构的X射线图像均衡显示处理方法。

7、Mary Couwenhoven，Robert Senn，David Foos，et al.“Enhancement Methodthat Provides Direct and Independent Control of Fundamental Attributes of ImageQuality for Radiographic Imagery”，Proc.Of SPIE Medical Imaging 2004，Vol.5367。

8、Pietro Perona and Jitendra Malik，“Scale-Space and Edge Detection UsingAnisoptropic Diffusion”.IEEE Transactions on Pattern Analysis and MachineIntelligence，Vol.12，NO.7，July 1990。

9、Ana Petrovic，Oscar Divorra Escoda and Pierre Vandergheynst，“Multiresolution。Segmentation of Natural Images：From Linear to NonlinearScale-Space Representations”，IEEE Transactions on Image Processing，Vol.13，NO.8。

10、Min Zhang.Xuanqin Mou.A Novel Contrast Equalization Method for ChestRadiograph.Proc.of SPIE.Vol.6144 61446R-1。

本申请参考文献5所述的图像多尺度分解与重建方法(但多尺度分解方法并非仅限于这一种)，该方法具体的步骤是：

a.基于尺度空间的图像分解

基于尺度空间的图像分解，可将图像分解为控制图像整体灰度范围的大区域信号和反应图像细节的小尺度信号。该分解通过模板卷积实现，分解过程具有可逆性，因而可以恢复得到原来的图像。

设分解层数为N，两幅原始图像为A、B，分解尺度依次为：σ₀，σ₁，σ₂，…，σ_N-1。

第1步：用尺度为σ₀的尺度函数G₀与原图像F₀(这里F₀代表A或B)进行模板卷积，得到F₁，F₀与F₁之差为D₀。

F₁(x，y)＝F₀(x，y)*G₀(x，y)；

D₀(x，y)＝F₀(x，y)-F₁(x，y)；    (1)

第i步(2≤i≤N)：用尺度为σ_i-1的尺度函数G_i-1对图像D_i-2进行平滑，得到F_i，D_i-2与F_i之差为D_i-1。

F_i(x，y)＝D_i-2(x，y)*G_i-1(x，y)；

D_i-1(x，y)＝D_i-2(x，y)-F_i(x，y).    (2)

b.分解分量重构

对原图像F₀，每一次分解过程都可以看作将待分解图像D_i-1(1≤i≤N-1)分解为D_i-1的高频部分D_i与低频部分F_i+1的过程。保留每次分解获得的F_i+1，以及D_N-1。设图像A和B经多尺度分解后需保留的分解分量分别为F₁ ^A，F₂ ^A，…，F_N ^A，D_N-1 ^A和F₁ ^B，F₂ ^B，…，F_N ^B，D_N-1 ^B。

由于低频分解分量F₁ ^A，F₂ ^A和F₁ ^B，F₂ ^B包含的主要是图像的背景信息，所以用加权平均法对A、B的低频分解分量F₁ ^A，F₂ ^A和F₁ ^B，F₂ ^B分别进行重构。

高频分量(这里指F₃ ^A，…，F_N ^A，D_N-1 ^A和F₃ ^B，…，F_N ^B，D_N-1 ^B)反映图像的细节信息。高频分量中数值较大的点代表该点为原图上变化较剧烈的点。对于A、B两幅图像分解后高频分量的重构，可分别采用相关性原则。

c.由重构分解分量重建图像

重构的分解分量为F₁，F₂，F₃，…，F_N，D_N-1，由这组分解分量可重建目标图像——具有大动态范围的X射线图像。为获得对重建图像的均衡显示效果，选择合适的均衡系数α_i(0＜α_i≤1)，分别对各分解分量进行不同程度的抑制或增强，加上均衡效果后，重建过程可表示为：

$F_0^{\′}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_i(x,y)\×{\mathrm{\α}}_i+D_{N-1}(x,y)\×{\mathrm{\α}}_{N+1}$

本发明的基于尺度空间分解的X射线图像超动态范围重建方法，根据此方法，对同一成像对象在不同成像条件下获得的两幅动态范围较窄的图像进行基于尺度空间的分解。将两幅图像的分解分量，在对应的频率上进行重组，由重组后的分解分量重建一幅动态范围较宽的重建图像。该方法的优点是将简单的基于尺度空间的分解方法应用到其中，在拓宽X射线图像动态范围、不降低图像信噪比的同时，能够有效地保持图像的灰度、细节和边缘信息。

通过实验证明，本发明具有良好的扩展X射线图像动态范围的效果。这一方法可解决采用碘化铯钠加光学镜头耦合CCD获取图像的DR设备成像动态范围小的问题。

附图说明

图1为X射线图像超动态范围重建方法的示意图

图2为基于尺度空间的图像分解示意图

图3为图像重建示意图

图4为X射线图像超动态范围重建方法的实验1的原图，实验原图的空间分辨率为1024*1024，灰度动态范围为14比特。

图5为使用本发明所述方法对实验1的原图处理的结果

图6为X射线图像超动态范围重建方法的实验2的原图，实验原图的空间分辨率为1024*1024，灰度动态范围为14比特。

图7为使用本发明所述方法对实验2的原图处理的结果

为了更清楚的理解本发明，下面结合附图的实验仿真结果作进一步的详细描述。

具体实施方式

本发明公开的X射线图像超动态范围重建方法，该方法的步骤包括：拍摄同一成像对象在不同成像条件下获得的两幅图像，这两幅图像成像动态范围小，分别包含成像对象不同区域的信息。首先对这两幅图像进行基于尺度空间的分解，然后通过一定的重构规则对两幅图像相应的分解分量进行重构，再对重构分量选择合适的均衡系数重建一幅具有宽动态范围的图像。其算法流程图如图1所示。分解过程可选择高斯滤波算子，其算法流程参见图2。对重构的分解分量选择合适的均衡系数，可压缩大尺度信号，增强小尺度信号，达到对重建图像的均衡显示。重建过程参见图3。

申请人分别用特制的对比度模体和人形模体作为成像对象拍摄的X线图像验证了本

发明的方法。

该方法具体的流程如下：

1.基于尺度空间的图像分解

设同一成像对象在不同成像条件下获得的两幅动态范围较窄的图像分别为A和B。这两幅图像是用同一X线成像设备在两个不同成像电压下获得的。

设分解层数为N，两幅原始图像为A、B，分解尺度依次为：σ₀，σ₁，σ₂，…，σ_N-1。由于唯一满足尺度空间中尺度定理的滤波器是高斯函数，所以选择高斯函数为尺度核函数。对于分解尺度的选择，考虑到小尺度滤波边缘定位准确对噪声较敏感，大尺度滤波边缘抗噪声能力强的特点，按照从大尺度到小尺度的顺序，及，σ_i+1＜σ_i，来分解不同尺度的图像内容。

第1步：用尺度为σ₀的尺度函数G₀与原图像F₀(这里F₀代表A或B)进行模板卷积，得到F₁，F₀与F₁之差为D₀。分解过程如图2(a)所示。

F₁(x，y)＝F₀(x，y)*G₀(x，y)；    (1)

D₀(x，y)＝F₀(x，y)-F₁(x，y)；

第i步(2≤i≤N)：用尺度为σ_i-1的尺度函数G_i-1对图像D_i-2进行平滑，得到F_i，D_i-2与F_i之差为D_i-1。分解过程如图2(b)所示。

F_i(x，y)＝D_i-2(x，y)*G_i-1(x，y)；    (2)

D_i-1(x，y)＝D_i-2(x，y)-F_i(x，y)；

2、分解分量的重构

对原图像F₀，每一次分解过程都可以看作将待分解图像D_i-1(1≤i≤N-1)分解为D_i-1的高频部分D_i与低频部分F_i+1的过程。保留每次分解获得的F_i+1，以及D_N-1。设图像A和B经多尺度分解后需保留的分解分量分别为F₁ ^A，F₂ ^A，…，F_N ^A，D_N-1 ^A和F₁ ^B，F₂ ^B，…，F_N ^B，D_N-1 ^B。

1)低频分量的重构

由于低频分解分量F₁ ^A，F₂ ^A和F₁ ^B，F₂ ^B包含的主要是图像的背景信息，可以采用但并不仅限于用加权平均法对A、B的低频分解分量F₁ ^A，F₂ ^A和F₁ ^B，F₂ ^B分别进行重构。加权平均法的具体计算公式为：

$F_k^C(i,j)=\mathrm{\α}\·F_k^A(i,j)+(1-\mathrm{\α})\·F_k^B(i,j),k=\mathrm{1,2}---\left(3\right)$

F_k ^C(i，j)为重建图像C相应的分解分量，α为权值，取值范围为[0，1]。

2)高频分量的重构

高频分量(这里指F₃ ^A，…，F_N ^A，D_N-1 ^A和F₃ ^B，…，F_N ^B，D_N-1 ^B)反映图像的细节信息。高频分量中数值较大的点代表该点为原图上变化较剧烈的点。对于A、B两幅图像分解后高频分量的重构，可分别采用但并不仅限于相关性原则。以A、B两幅图第M(3≤M≤N)层分解的高频分量F_M ^A、F_M ^B为例说明高频分解分量重构的相关性方法。

首先分别计算高频分量F_M ^A、F_M ^B中任一象素点的区域能量。设(m，n)为高频分量中的任一象素点，f(m，n)为其象素值，以该点为中心选取一个小区域(本文中选择3×3的区域)计算区域能量E(m，n)。具体计算公式为：

$E(m,n)=\underset{i=-1}{\overset{1}{\mathrm{\Σ}}}\underset{j=-1}{\overset{1}{\mathrm{\Σ}}}{f}^2(m+i,n+j)---\left(4\right)$

用C_AB(m，n)表示F_M ^A和F_M ^B在以象素点(m，n)为中心的3×3区域内的相关性，C_AB(m，n)的计算公式如(5)所示：

$C_{\mathrm{AB}}(m,n)=\frac{2\underset{i=-1}{\overset{1}{\mathrm{\Σ}}}\underset{j=-1}{\overset{1}{\mathrm{\Σ}}}{F}_M^A(m+i,n+j)*F_M^B(m+i,n+j)}{E^A(m,n)+E^B(m,n)}---\left(5\right)$

选择阈值λ，当|C_AB(m，n)|≤λ时，直接选择局域能量较大的象素值作为重建图像对应分解分量的象素值。当|C_AB(m，n)|＞λ时，重建图像对应分解分量的象素值可通过对F_M ^A、F_M ^B象素值的加权运算获得。设F_M ^A、F_M ^B的加权因子分别为W^A、W^B，由式(6)确定：

$W^A(m,n)=\left\{\begin{array}{cc}0.5*(1-\frac{1-C_{\mathrm{AB}}(m,n)}{1-\mathrm{\&lambda;}})& E^A(m,n)<E^B(m,n)\\ 0.5*(1+\frac{1-C_{\mathrm{AB}}(m,n)}{1-\mathrm{\&lambda;}})& E^A(m,n)\&GreaterEqual;E^B(m,n)\end{array}\right.---\left(6\right)$

W^B(m，n)＝1-W^A(m，n)

3、由重构分量重建图像

对原始图像F₀的每一次分解过程都可以看作将待分解图像D_i-1(1≤i≤N-1)分解为D_i-1的高频部分D_i与低频部分F_i+1的过程。保留每次分解获得的F_i+1，以及D_N-1，即可恢复原始图像F₀，恢复得到的图像用F₀′表示。重建过程如图3所。

F₀(x，y)＝D₀(x，y)+F₁(x，y)

D₀(x，y)＝D₁(x，y)+F₂(x，y)

D₁(x，y)＝D₂(x，y)+F₃(x，y)



D_N-2(x，y)＝D_N-1(x，y)+F_N(x，y)

$F_0^{\&prime;}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_i(x,y)+D_{N-1}(x,y)---\left(7\right)$

本发明中将重构的分解分量记做F₁，F₂，F₃，…，F_n，D_N-1，由这组分解分量根据式(7)可重建目标图像--动态范围扩展的X射线图像。为获得对重建图像的均衡显示效果，还可选择合适的加权系数α_i(0＜α_i-1≤α_i≤1，i＝1∶N+1)，分别对各分解分量进行不同程度的抑制或增强，此时，重建过程可表示为：

$F_0^{\&prime;}(x,y)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_i(x,y)\&times;{\mathrm{\&alpha;}}_i+D_{N-1}(x,y)\&times;{\mathrm{\&alpha;}}_{N+1}---\left(8\right)$

实验中，首先用特制的对比度测试模体和模拟人体胸部不同部位的阶梯型树脂模体作为成像对象，采用基于影像增强器成像链的X光机分别拍摄两幅动态范围不同的图像作为原始图像(参见图4)验证本发明的重建方法。实验结果如图5所示。从图中可以看出，两幅原始图像均无法同时显示成像对象三个部分的细节，而重建图像可以显示成像对象三个部分中的细节，重建图像的动态范围比原始图像的动态范围大，但灰度范围与原始图像不同。

在随后实验中，用人形模体在同一X光机上做了相同实验。由于实验所用人形模体的肺部没有明显的细节，将对比度模体放置在人形模体肺野下方，作为肺部的细节。将人形模体放置在检测床上，在50kVp，100mA，63ms和75kVp，100mA，16ms的条件下分别曝光，获得一幅低电压图像和一幅高电压图像。低电压图像只能显示心脏和肺野部分，高电压图像只能显示心脏和纵膈部分，如图6所示。用本发明所述方法得到重建图像。实验结果如图7所示。重建图像能同时显示心脏、肺野和纵膈三部分信息。

综上所述，本发明的X射线图像超动态范围重建方法针对X射线图像自身特点，通过实际应用表明，能够兼顾图像整体对比度调整及动态范围增大，从而有利于医用诊断及工业使用，具有实际应用价值。

Claims (3)

1.一种基于尺度空间分解的X射线图像超动态范围重建方法，其特征在于，该方法具体的步骤是：

A、对同一成像对象放置在两个不同球管电压下，该两个球管电压分别选择成像对象的两个不同感兴趣区域，在两个不同球管电压下，以合适的曝光条件分别曝光，得到两幅图像，该两幅图像分别显示成像对象的不同区域信息；

B、将拍摄到的两幅图像分别进行基于尺度空间的图像分解，将两幅图像A、B分解为控制图像整体灰度范围的大区域信号和反应图像细节的小尺度信号；其中大区域信号分解后得到大尺度分量即低频分量，主要包含图像的背景信息，小尺度信号分解后得到小尺度分量即高频分量，反映图像的细节信息；

C、在低频分量和高频分量所对应的频率上进行重组，由重组后的分解分量重建一幅动态范围较宽的重建图像。

2.如权利要求1所述的方法，其特征在于，所述的两幅图像的基于尺度空间的图像分解多尺度分解的方法选择高斯函数为尺度核函数，具体包括下列步骤：

设分解层数为N，两幅原始图像为A、B，分解尺度依次为：σ₀，σ₁，σ₂，…，σ_N-1；

第1步：用尺度为σ₀的尺度函数G₀与原图像F₀进行模板卷积得到F₁，F₀代表两幅原始图像A或B，F₀与F₁之差为D₀；如下式所示：

F₁(x，y)＝F₀(x，y)*G₀(x，y)；

D₀(x，y)＝F₀(x，y)-F₁(x，y)；    (1)

第i步：用尺度为σ_i-1的尺度函数G_i-1对图像D_i-2进行平滑，得到F_i，D_i-2与F_i之差为D_i-1；如下式所示：

F_i(x，y)＝D_i-2(x，y)*G_i-1(x，y)；

D_i-1(x，y)＝D_i-2(x，y)-F_i(x，y).    (2)

3.如权利要求1所述的方法，其特征在于，所述的低频分量和高频分量用加权平均法对两幅图像A、B的低频分解分量进行重构，具体计算公式为：

$F_k^C(i,j)=\mathrm{\&alpha;}\&CenterDot;F_k^A(i,j)+(1-\mathrm{\&alpha;})\&CenterDot;F_k^B(i,j),k=\mathrm{1,2}$

式中，F_k ^C(i，j)为重建图像C相应的分解分量，α为权值，取值范围为[0，1]；

采用相关性原则对两幅图像A、B分解的高频分量进行重构，以两幅图像A、B的第M层分解的高频分量F_M ^A、F_M ^B为例说明高频分解分量重构的具体方法：

首先分别计算高频分量F_M ^A、F_M ^B中任一象素点的区域能量。设(m，n)为高频分量中的任一象素点，f(m，n)为其象素值，以该点为中心选取一个小区域(本文中选择3×3的区域)计算区域能量E(m，n)。具体计算公式为：

$E(m,n)=\underset{i=-1}{\overset{1}{\mathrm{\&Sigma;}}}\underset{j=-1}{\overset{1}{\mathrm{\&Sigma;}}}{f}^2(m+i,n+j)$

用C_AB(m，n)表示F_M ^A和F_M ^B在以象素点(m，n)为中心的3×3区域内的相关性，C_A ^B(m，n)的计算公式为：

$C_{\mathrm{AB}}(m,n)=\frac{2\underset{i=-1}{\overset{1}{\mathrm{\&Sigma;}}}\underset{j=-1}{\overset{1}{\mathrm{\&Sigma;}}}{F}_M^A(m+i,n+j)*F_M^B(m+i,n+j)}{E^A(m,n)+E^B(m,n)}$

选择阈值λ，当|C_AB(m，n)|≤λ时，直接选择局域能量较大的象素值作为重建图像对应分解分量的象素值。当|C_AB(m，n)|＞λ时，重建图像对应分解分量的象素值可通过对F_M ^A、F_M ^B象素值的加权运算获得。设F_M ^A、F_M ^B的加权因子分别为W^A、W^B，由下式确定：

$W^A(m,n)=\left\{\begin{array}{cc}0.5*(1-\frac{1-C_{\mathrm{AB}}(m,n)}{1-\mathrm{\&lambda;}})& E^A(m,n)<E^B(m,n)\\ 0.5*(1+\frac{1-C_{\mathrm{AB}}(m,n)}{1-\mathrm{\&lambda;}})& E^A(m,n)\&GreaterEqual;E^B(m,n)\end{array}\right.$

W^B(m，n)＝1-W^A(m，n)。

Images