CN104871208A - 用于对图像进行滤波的图像处理装置及方法 - Google Patents

Abstract

本发明涉及一种用于对图像进行滤波的图像处理装置。所述装置包括图像输入部(3)，所述图像输入部用于获得相同对象的第一图像和第二图像，所述第一图像和所述第二图像包括多个体素并且通过噪声协方差而互相关，每个体素具有包括信号值以及噪声值的体素值。联合双边滤波器(3)被提供为对所述第一图像进行滤波，其中，通过包括相对体素特异的权重的滤波函数来对所述第一图像的第一体素进行滤波，所述权重包括获得所述第一图像中的所述第一体素的体素值和第二体素的体素值以及所述第二图像中的第一体素和第二体素的体素值的可能性，假设所述第一图像的所述第一体素的信号值与所述第一图像的第二体素的信号值相同，并且所述第二图像的所述第一体素的信号值与所述第二图像的第二体素的信号值相同。所述经滤波的图像被提供在图像输出部(4)处。

Description

用于对图像进行滤波的图像处理装置及方法

技术领域

本发明涉及一种用于对图像进行滤波的图像处理装置以及相对应的方法。

背景技术

出现了不仅提供患者或对象的单幅图像而且提供所述患者或对象的在完全相同的几何结构中的两幅或多幅的新的成像模态。针对这的两个范例是：微分相位对比成像(DPCI；例如在US 2012/0099702 A1中描述的)，其中，复折射率的实部和虚部被测量；以及能谱CT(例如在US 2008/0253504A中描述的)，其中，由于光电效应、康普顿散射以及可能由于对比剂的X射线衰减被修复。

在WO 2011/145040 A1中，联合双边滤波已经被提出作为用于基于图像的去噪的公知的双边滤波的泛化。尤其是，图像处理装置被描述为包括：图像提供单元，所述图像提供单元用于提供对象的第一图像以及相同对象的第二图像；以及滤波单元，其基于在第一图像的图像值之间的第一相似度并且基于在第二图像的图像值之间的第二相似度，来对第一图像进行滤波。即使第一图像和第二图像中的一幅的图像值由于噪声而被干扰，这也允许基于图像值属于对象的相同部分的可能性(例如，属于对象的组织或骨材料的可能性，如果对象是人类或动物)来对第一图像进行滤波，从而例如改进了CT图像的边缘保持的基于图像的去噪。

然而，已知方法和设备并未认识噪声在各幅图像中能够是相关的。尤其是，对于DPCI而言，在相位对比图像内的噪声的经轴图像(其等价于具有在低频处的峰值的公知噪声功率谱)中有强的空间相关性，并且对于能谱CT而言，在光图像和康普顿图像中的相同空间位置处有噪声的强的反相关性。

还可能存在其中相同对象的两幅图像可利用但是处于不同空间分辨率处的其他情况，所述不同空间分辨率也由噪声中的不同空间相关性反映。一个范例可以是磁共振成像(MRI)，其中，在一些情况下，解剖图像(例如，T2加权的)以及参数图像(T1)被采集。

发明内容

本发明的目的是提供用于对图像进行滤波的改进的图像处理装置及方法，以及尤其用于联合基于图像的去噪的方法。

在本发明的第一方面，提出了一种用于对图像进行滤波的图像处理装置，所述装置包括：

-图像输入部，其用于获得相同对象的第一图像和第二图像，所述第一图像和第二图像包括多个体素并且通过噪声协方差而互相关，每个体素具有包括信号值以及噪声值的体素值，

-联合双边滤波器，其用于对所述第一图像进行滤波，其中，通过包括相对体素特异的权重的滤波函数对所述第一图像的第一体素进行滤波，所述权重包括获得所述第一图像中的所述第一体素的体素值和第二体素的体素值以及所述第二图像中的第一体素的体素值和第二体素的体素值的可能性，假设所述第一图像的所述第一体素的信号值与所述第一图像的第二体素的信号值相同，并且所述第二图像的所述第一体素的信号值与所述第二图像的第二体素的信号值相同，

-图像输出部，其用于提供经滤波的图像。

在本发明的又一方面中，提供了一种对应的图像处理方法。

在本发明的又一方面中，提供了一种包括程序代码模块的计算机程序，当所述计算机程序在计算机上运行时，所述程序代码模块用于令所述计算机执行处理方法的步骤，并且提供了一种存储有计算机程序产品的非暂态计算机可读记录介质，所述计算机程序产品在由处理器运行时，令本文公开的方法被执行。

从属权利要求中定义了本发明的优选实施例。应理解，要求保护的方法、计算机程序以及介质具有与要求保护并且在从属权利要求中定义的装置相似和/或相同的优选实施例。

根据本发明，提出了一种用于改进的联合双边滤波的方法。所提出的方法基本上包含两项：距离加权和范围加权。该第一项是纯粹几何促动的，并且不被改变。范围滤波引入基于图像值的加权，所述加权被提出为根据本发明而被改变。尤其是在图像之间以及在不同图像位置之间的噪声相关性根据本发明被考虑，以提供改进的基于图像的去噪。

在该背景下，“互相关”应当被理解为使得相关性能够被限制到在相同图像的体素之间的相关性，或者其可以被限制到在不同图像中的相同体素位置之间的相关性。

在优选的实施例中，所述滤波器被配置为使用可能性，所述可能性包括描述所述噪声协方差的噪声协方差矩阵。通过使用所述协方差，能够解决例如其中噪声在一个方向(例如，在针对DPCI的面内)中高度相关但是在z方向中高度不相关的情况。在这样的情况下，常规(联合)双边滤波表现不佳，因为范围滤波常数能够被调谐到仅面内或面外噪声特性。

优选地，所述滤波器被配置为使用以下形式的协方差矩阵：

$C=\left(\begin{array}{cccc}{c}_{11}& c_{12}& c_{13}& c_{14}\\ c_{12}& c_{22}& c_{23}& c_{24}\\ c_{13}& c_{23}& c_{33}& c_{34}\\ c_{14}& c_{24}& c_{34}& c_{44}\end{array}\right)$

或者

$C=\left(\begin{array}{cccc}{c}_{11}& c_{12}& c_{13}& c_{14}\\ c_{12}& c_{11}& c_{14}& c_{13}\\ c_{13}& c_{14}& c_{33}& c_{34}\\ c_{14}& c_{13}& c_{34}& c_{33}\end{array}\right)$

这些协方差矩阵已经被示出为是有利的。第二协方差矩阵更容易操作。

以下特殊情况是尤其有利的：对于微分相位对比成像，在经重建的折射率的虚部(第一图像)与实部(第二图像)之间几乎没有任何噪声相关性，造成(命名涉及C的一般版本，即，上面示出的两个版本的较上的协方差矩阵)c₁₃≈c₂₃≈c₁₄≈c₂₄≈0，然而在折射率的实部中在经轴图像内存在强的噪声相关性，即，c₃₄＞0。协方差矩阵当然取决于在体素位置上，即如果第二图像的第二体素在不同经轴图像中，那么噪声是不相关的，即c₃₄≈0。

对于能谱成像，在每幅图像(即，康普顿图像以及光电图像)内仅存在较小噪声空间相关性，实现c₁₂≈c₃₄≈0，并且在不同图像中的不同位置之间也没有相关性，但是在相同位置处在两幅图像中的图像值之间有强的负相关性，并且c₁₃＜0，c₂₄＞0。

此外，在实施例中，所述滤波器被配置为使用以下形式的可能性：

$p({\mathrm{\Δ}}_a,{\mathrm{\Δ}}_b|a_1=a_2,b_1=b_2)\∝F\exp \left(\frac{1}{4}\right({\mathrm{\Δ}}_a^2(s_{12}-s_{11})+2{\mathrm{\Δ}}_a{\mathrm{\Δ}}_b(s_{14}-s_{13})+{\mathrm{\Δ}}_b^2(s_{34}-s_{33})\left)\right)$

其中，

$F=\frac{\sqrt{\mathrm{\π}(s_{11}+s_{12})}}{\sqrt{(s_{11}+s_{12})(s_{33}+s_{34})}-{(s_{13}+s_{14})}^2}$

$S=\left(\begin{array}{cccc}{s}_{11}& s_{12}& s_{13}& s_{14}\\ s_{12}& s_{11}& s_{14}& s_{13}\\ s_{13}& s_{14}& s_{33}& s_{34}\\ s_{14}& s_{13}& s_{34}& s_{33}\end{array}\right)$

S是所述协方差矩阵的逆矩阵，

a₁和a₂是所述第一图像的所述第一体素值和第二体素值，

b₁和b₂是所述第二图像的所述第一体素值和第二体素值，并且

Δ_a＝a₁-a₂，Δ_b＝b₁-b₂。

更进一步，所述第一图像的所述经滤波的第一体素值被获得为：

其中，

$N_i=\underset{j\∈N_i}{\Σ}{w}_d(i,j)p(a_i-a_j,b_i-b_j|a_i=a_j,b_i=b_j)$

并且，

w_d是标准(联合)双边滤波中的(常用)距离权重。可以使用随着增加的距离|i-j|而单调地减小的任何函数，例如：

$w_d(i,j)=e^{-{(i-j)}^2/2{\mathrm{\σ}}^2}$

或者

$w_d(i,j)=\frac{1}{{(i-j)}^2+\mathrm{\ϵ}}$

在实际情况下，所述第一图像和第二图像中的至少一幅包括非白噪声，在所述第一图像与第二图像之间存在噪声相关性，和/或在所述第一图像和/或所述第二图像的所述第一体素和第二体素处的噪声噪声协方差相同。

优选地，相同滤波器被应用于所述第二图像。此外，提出的想法能够一直扩展到具有三幅(或更多幅)图像的情况。

根据另一个实施例，所述图像输入部被配置为提供所述第一图像以及所述第二图像，使得它们示出所述对象的不同性质。优选地，所述图像输入部被配置为提供能谱计算机断层摄影图像，其中，所述第一图像指示所述对象的第一物理性质并且所述第二图像指示所述对象的第二物理性质，所述第一物理性质和所述第二物理性质是彼此不同的。

例如，所述第一图像和所述第二图像中的一幅是康普顿图像，并且所述第一图像和所述第二图像中的另一幅是光电图像。在另一个实施例中，所述第一图像和所述第二图像中的一幅是康普顿图像或光电图像，并且所述第一图像和所述第二图像中的另一幅是康普顿图像与光电图像的组合。

更进一步，所述图像输入部被配置为提供微分相位对比计算机断层摄影图像，其中，所述第一图像和所述第二图像中的一幅指示所述对象的吸收，并且所述第一图像和所述第二图像中的另一幅指示由所述对象引起的辐射的相位移位，所述辐射穿过所述对象。

附图说明

本发明的这些和其他方面将根据下文描述的(一个或多个)实施例变得显而易见并且将参考下文描述的(一个或多个)实施例得到阐述。在以下附图中：

图1示意性并且示范性地示出了根据本发明的图像处理装置的实施例，

图2示意性并且示范性地示出了图像输入单元的实施例，

图3示意性并且示范性地示出了多色辐射源的发射谱，

图4示意性并且示范性地示出了在显示的能量范围内的光电效应的横截面、康普顿效应以及具有k边的材料的总衰减的能量依赖性，

图5示意性并且示范性地示出了微分相位对比计算机断层摄影系统的实施例的元件，并且

图6示意性并且示范性地示出了图示用于处理图像的图像处理方法的实施例。

具体实施方式

图1示意性地示出了根据本发明的图像处理装置1的实施例。其包括图像输入部2，所述图像输入部用于获得相同对象的第一图像和第二图像，第一图像和第二图像包括多个体素并且通过噪声协方差互相关，每个体素具有包括信号值和噪声值的体素值。装置1还包括用于对第一图像进行滤波的双边滤波器3，其中，第一图像的第一体素由包括相对体素特异的权重的滤波函数进行滤波。所述权重包括获得在第二图像中的所述第一体素的体素值和第二体素的体素值的可能性以及在第二图像中的第一体素的体素值和第二体素的可能性，假设第一图像的所述第一体素的信号值与第二图像的第一体素的信号值相同，并且第一图像的第二体素的信号值与第二图像的第二体素的信号值相同。装置1还包括用于提供所述经滤波的图像的图像输出部4，所述图像输出部例如可以是用于显示经滤波的图像的显示器。

第一图像和第二图像能够是示出相同对象的二维、三维或四维图像，所述相同对象能够是整个对象或仅是对象的部分，例如仅是对象的切片。如果对象是诸如人的心脏的运动对象，则第一图像和第二图像示出优选地在相同运动状态的对象。该实施例，第一图像和第二图像是基于投影数据，尤其是基于计算机断层摄影投影数据来重建的，其中，第一图像和第二图像根据相同的投影数据重建。

图像输入单元2适于提供第一图像和第二图像，使得所述第一图像和第二图像示出对象的不同性质，尤其是不同物理性质。在实施例中，图像处理单元2适于提供能谱计算机断层摄影图像，其中，第一图像指示对象的第一物理性质并且第二图像指示对象的第二物理性质，其中，第一物理性质和第二物理性质彼此不同。例如，第一图像和第二图像中的一幅能够是康普顿图像并且第一图像和第二图像中的另一幅能够是光电图像，或者第一图像和第二图像中的一幅能够是康普顿图像或光电图像并且第一图像和第二图像中的另一幅能够是康普顿图像和光电图像的组合。

图像输入单元2能够是其中第一图像和第二图像被存储的存储单元，或者图像输入单元2还能够是允许经由数据连接来接收第一图像和/或第二图像的图像接收单元。图像输入单元2还能够是生成第一图像和第二图像的成像模态，如能谱计算机断层摄影系统。在下文中将参考图2示意性并且示范性地描述这样的能谱计算机断层摄影系统。

能谱计算机断层摄影系统包括机架14，所述机架围能够围绕旋转的轴R进行旋转，所述旋转的轴R平行于z方向延伸。辐射源15，例如X射线管，被安装在机架14上。在该实施例中，辐射源15发射多色辐射。辐射源15被提供具有准直器设备16，所述准直器设备从由辐射源15发射的辐射形成圆锥形辐射射束17。在其他实施例中，准直器设备16能够适于形成具有另一形状，例如具有扇形形状的辐射射束。

辐射穿过在圆柱形检查区5中的感兴趣区域中的对象(未在示出图2中)，例如患者或技术性对象。在已经穿过了感兴趣区域之后，辐射射束17入射在探测设备6上，所述探测设备在该实施例中具有二维探测表面，其中，探测设备6被安装在机架14上。在另一个实施例中，探测设备6能够包括一维探测表面。

探测设备6优选地包括能量解析探测器元件。所述能量解析元件优选地是光子计数探测器元件，其根据对入射光子进行计数的原则工作并且输出示出在不同能量窗中的光子的数量的信号。这样的能量解析探测设备例如描述在Llopart X等有的“First test measurements of a 64k pixel readout chipworking in a single photon counting mode”，Nucl.Inst.and Meth.A，509(1-3)、第157至163页(2003年)中并且在Llopart,X等的“Medipix2:A64-k pixel readout chip with 55μm square elements working in a single photoncounting mode”，IEEE Trans.Nucl.Sci，49(5)、第2279至2283页(2002年)中。优选地，能量解析探测器元件被适配为使得每个探测器元件提供针对至少两个不同能量窗的至少两个能量解析的探测信号。然而，具有甚至更加高的能量分辨率是更有利的，以便增强成像系统的敏感度和噪声鲁棒性。

机架14由电动机7以优选地固定的但是可调节的角速度被驱动。另外的电动机8被提供为平行于旋转轴或z轴的方向移动对象，例如被布置在检查区5中的患者台上的患者。这些电动机7、8例如由控制单元9控制，使得辐射源15和检查区5，尤其是感兴趣区域，沿螺旋形轨迹相对于彼此移动。还能够是对象或检查区5，尤其是感兴趣区域，并不移动，并且辐射源15被旋转，即辐射源15相对于感兴趣区域沿圆形轨迹行进。由探测设备6采集的数据在该实施例中是能量依赖的探测数据，所述能量依赖的探测数据被提供到重建单元10，以根据所提供的能量依赖的探测数据来重建感兴趣区域的至少第一图像和第二图像。重建单元10也优选地由控制单元9所控制。经重建的图像被提供到滤波单元3，所述滤波单元用于对至少第一图像进行滤波。

重建单元10包括计算单元12，所述计算单元用于根据探测数据来确定至少两个衰减分量。探测数据的衰减分量由仅一个或由若干个物理效应所引起，例如，一个衰减分量能够是由康普顿效应引起的分量，并且另一个衰减分量能够是由光电效应引起的分量。另外的衰减分量能够是由存在于感兴趣区域中的k边引起的分量。备选地或者额外地，衰减分量能够是例如由在感兴趣区域内的特定材料的吸收所引起的分量。例如，一个衰减分量例如能够是由特定材料的吸收引起的分量，并且另一个衰减分量能够是由另一种材料的吸收引起的分量。

重建单元10还包括反向投影单元13，所述反向投影单元用于对探测数据的计算的衰减分量进行反向投影，使得对应的衰减分量图像被重建。例如，康普顿图像能够被重建为第一图像和第二图像中的一幅，并且光电图像能够被重建为第一图像和第二图像中的另一幅。

对计算单元12的输入是多个，在该实施例中至少三个，能量箱的能量解析的探测数据d_i。这些探测数据d_i示出了第i个能量箱e_i的谱敏感度D_i(E)。此外，多色辐射源15的发射谱T(E)通常是已知的或能够被测量。多色辐射源的这样的发射谱T(E)的范例示意性并且示范性地示出在图3中。在计算单元12中，对探测数据d_i的生成被建模为具有能谱P(E)的光电效应、具有能谱B(E)的康普顿效应及具有能谱K(E)的k边的线性组合。

能谱P(E)、B(E)和K(E)示范性并且示意性地示出在图4中。

对探测数据的生成能够由以下方程式组所建模：

d_i＝∫dET(E)D_i(E)exp[-(ρ_光电P(E)+ρ_康普顿B(E)+ρ_k边K(E))]，    (1)

其中，ρ_光电、ρ_康普顿和ρ_k边分别是光电分量、康普顿分量和k边分量的密度长度乘积。

由于至少三个探测信号d₁、d₂、d₃可用于所述至少三个能量箱e₁、e₂、e₃，所以至少三个方程式的组被形成为具有三个未知量，所述三个未知量是三个密度长度乘积，因此能够在计算单元12中利用已知数值方法来对所述三个未知量进行求解。如果多于三个能量箱可利用，则优选地使用将测量值的噪声统计考虑在内的最大似然方法。通常，三个能量箱是足够的。然而，为了增加敏感度和噪声鲁棒性，优选地具有针对更多能量箱的更多探测数据。针对每个衰减分量，即针对每个密度长度乘积，能够通过对各自的密度长度乘积进行反向投影来重建对象的衰减分量图像。备选地或额外地，不同密度长度乘积的图像，或者在反向投影之前，不同密度长度乘积能够被组合以生成组合图像，所述组合图像是由组合的衰减分量生成的对象的图像。

在另一个实施例中，图像输入单元2能够适于提供微分相位对比计算机断层摄影图像，其中，第一图像和第二图像中的一幅指示对象的吸收，并且第一图像和第二图像中的另一幅指示由对象引起的穿过对象的辐射的相位移位。尤其是，第一图像优选地示出对象的总衰减系数，并且第二图像优选地示出对象的折射率的实部。同样，这些第一图像和第二图像能够被存储在存储单元中，所述存储单元可以是图像输入单元2，或者能够经由图像接收单元来接收这些第一图像和第二图像，所述图像接收单元也可以是图像输入单元2。然而，图像输入单元2也能够是微分相位对比计算机断层摄影系统，所述微分相位对比计算机断层摄影系统生成第一图像和第二图像，并且将在下文中参考图5示范性并且示意性地描述所述微分相位对比计算机断层摄影系统。

图5示出了微分相位对比计算机断层摄影系统的示范性实施例的三维表示。相当大的X射线源32被布置为邻近于源光栅34。因为X射线源32关于发射的辐射的波长由于其尺寸而可以被认为是不相干的，因而源光栅34被采用以提供多个单个相干X射线源。

X辐射35从X射线源32在光学轴37的方向中发出，可能构成X射线的扇形射束或锥形射束。X射线射束的各自的形状未描绘在图5中。X辐射35到达对象36，穿过对象36并且随后到达射束分离器光栅38。射束分离器光栅38的沟或间隙变更通过的电磁辐射关于射束分离器光栅的固体区(阻挡区域)的相位。相应地，执行经过φ的，尤其是经过π的相位移位。

分析器光栅被布置在射束分离器光栅38与X射线探测器42之间。在X射线探测器的方向中的来自射束分离器光栅38的多个波到达分析器光栅40，并且随后在X射线探测器42的表面上产生强度调制图样。

通过将分析器光栅40对射束分离器光栅38进行移位，因此相对于彼此移动光栅，尤其是利用射束分离器光栅38或分析器光栅40的光栅周期的分数，多个发送图像可由图像探测器42所获得。因此，在对发送图像的采集期间，所谓的相位步进被执行，其中，光栅38、40的相对位置例如在一个周期上以四至十步地变化，并且针对每个相对位置测量发送图像。根据该发送图像系列，通过针对每个探测器像素对所测量的强度进行简单的平均来获得常规吸收图像。如在常规X射线发送成像中，该测量涉及按比尔定律的沿X射线路径的对总衰减系数的线积分。通过分析在相位步进周期上测量的强度的动态，来获得与折射率的实部有关的信息。强度优选地利用分析器光栅40的光栅周期周期性地变化。通常，能够假设正弦变化。弯曲的相位涉及撞击射束分离器光栅38的相位波前的梯度。物理上，相位波前的该梯度与对象的折射率的实部的投影的一阶导数线性相关。因此能够基于弯曲的相位来确定对象的折射率的实部的投影。通过将对象的折射率的实部的所确定的投影进行反向投影，示出对象的折射率的实部的图像能够被重建。微分相位对比计算机断层摄影系统因此例如能够重建示出对象的总衰减系数的第一图像以及示出对象的折射率的实部的第二图像。对已知微分相位对比计算机断层摄影技术的更加详细的描述例如被公开在Timm Weitkamp等人的文章“X-ray phase imaging with a gratinginterferometer”，OPTICS EXPRESS，第13卷、第16号、第6296至6304页(2005年8月)中，通过引用将所述文章整体并入本文。

在WO 2011/145040 A1中，联合双边滤波已经被提出作为用于基于图像的去噪的公知的双边滤波的泛化。在该方法中，滤波单元优选地适于执行滤波，使得通过以下滤波方程来对第一图像去噪：

${\hat{u}}^{\left(1\right)}\left({\stackrel{\→}{x}}_i\right)=\frac{\underset{{\stackrel{\‾}{x}}_j\∈\mathrm{\Δ}i}{\Σ}{L}_1({\stackrel{\→}{u}}^{\left(1\right)}\left({\stackrel{\→}{x}}_i\right)-{\stackrel{\→}{u}}^{\left(1\right)}\left({\stackrel{\→}{x}}_j\right))L_2({\stackrel{\→}{u}}^{\left(2\right)}\left({\stackrel{\→}{x}}_i\right)-{\stackrel{\→}{u}}^{\left(2\right)}\left({\stackrel{\→}{x}}_j\right))w_d({\stackrel{\→}{x}}_i-{\stackrel{\→}{x}}_j)u^{\left(1\right)}\left({\stackrel{\→}{x}}_j\right)}{\underset{{\stackrel{\‾}{x}}_j\∈\mathrm{\Δ}i}{\Σ}{L}_1({\stackrel{\‾}{u}}^{\left(1\right)}\left({\stackrel{\‾}{x}}_i\right)-{\stackrel{\‾}{u}}^{\left(1\right)}\left({\stackrel{\→}{x}}_j\right))L_2({\stackrel{\→}{u}}^{\left(2\right)}\left({\stackrel{\→}{x}}_i\right)-{\stackrel{\→}{u}}^{\left(2\right)}\left({\stackrel{\→}{x}}_j\right))w_d({\stackrel{\→}{x}}_i-{\stackrel{\→}{x}}_j)},---\left(2\right)$

其中，指代针对在第一图像的位置处的第一图像元素的平均图像值，Δ_i指代包括在第一位置处的第一图像元素的第一区域，指代在第一区域Δ_i内的第一区域图像元素的位置，指代在a)第一区域图像元素的图像值以及任选地围绕在位置处的各自的第一区域元素的图像元素的图像值与b)在位置处的第一图像元素的图像值以及任选地围绕第一图像的第一图像元素的图像元素的图像值之间的第一相似度。

函数指代在以下a)与b)之间的第二相似度：a)第二图像的第二区域图像元素在第二区域内的位置处的图像值，所述第二区域对应于第一区域；以及任选地，第二区域元素的围绕在位置处的图像元素的图像值；b)第二图像的第二图像元素在位置处的的图像值，所述第二图像元素对应于在第一图像中的相同位置处的第一图像元素；以及任选地，围绕第二图像元素的图像元素的图像值。函数是基于在位置处的图像元素与位置处的图像元素之间的距离的距离函数，并且指代在第一图像中的位置处的未经滤波的图像值。

第一相似度L₁、第二相似度L₂以及距离函数w_d的组合除以示出在方程式(2)中的分母能够被认为是用于执行根据方程式(2)的加权平均的权重。因此，根据方程式(2)，滤波单元3优选地适于对第一图像的图像元素的图像值进行加权平均，其中，为了确定针对第一图像的第一图像元素的平均图像值，在包括第一图像元素的第一区域内的第一区域图像元素的图像值被加权并且被平均，其中，利用取决于第一相似度并且取决于第二相似度的权重对第一区域图像元素的图像值进行加权，第一相似度是以下a)与b)之间的相似度：a)第一区域图像元素的图像值以及任选地围绕各自的第一图像元素的图像元素的图像值与b)第一图像元素的图像值以及围绕第一图像的第一图像元素的图像元素的图像值，第二相似度是在以下a)与b)之间的相似度：a)第二图像的第二区域内的第二区域图像元素的图像值，所述第二区域对应于所述第一区域；以及任选地，围绕各自的第二区域图像元素的图像元素的图像值；b)第二图像的第二图像元素的图像值，所述第二图像元素对应于第一图像元素；以及任选地，围绕第二图像元素的图像元素的图像值。

向量能够被认为是以为中心围绕的图像片块。例如，如果围绕各自的第一图像元素的图像元素的图像元素的任选的图像值以及围绕各自的第二图像元素的图像元素的任选的图像值被考虑，则在二维情况下，图像片块可以是以在位置处的像素为中心围绕在位置处的像素的9x9子图像。片块也能够具有另一尺寸。例如，片块能够具有5x5或7x7的尺寸。此外，三维片块也能够被使用。如果所述任选的图像值不被考虑，即如果仅在位置处与位置处的单个图像值之间的差异被考虑，则各自的向量分别仅包括在位置处或处的各自的单个图像值。片块尺寸则是1x1。

如果基于多于两幅图像来执行滤波，则滤波单元3优选地适于根据以下滤波方程来对第一图像进行滤波：

${\hat{u}}^{\left(1\right)}\left({\stackrel{\→}{x}}_i\right)=\frac{\underset{{\stackrel{\‾}{x}}_j\∈{\mathrm{\Δ}}_i}{\Σ}\underset{k=1...N}{\Π}{L}_k({\stackrel{\→}{u}}^{\left(k\right)}\left({\stackrel{\→}{x}}_i\right)-{\stackrel{\→}{u}}^{\left(k\right)}\left({\stackrel{\→}{x}}_j\right))K({\stackrel{\→}{x}}_i-{\stackrel{\→}{x}}_j)u^{\left(1\right)}\left({\stackrel{\→}{x}}_j\right)}{\underset{{\stackrel{\‾}{x}}_j\∈{\mathrm{\Δ}}_i}{\Σ}\underset{k=1...N}{\Π}{L}_k({\stackrel{\→}{u}}^{\left(k\right)}\left({\stackrel{\→}{x}}_i\right)-{\stackrel{\→}{u}}^{\left(k\right)}\left({\stackrel{\→}{x}}_j\right))K({\stackrel{\→}{x}}_i-{\stackrel{\→}{x}}_j)},---\left(3\right)$

其中，N指代对第一图像进行滤波所考虑的图像的数量，并且指数k指示各自的元素属于第k图像。

第k相似度能够由以下方程式定义：

$L_k({\stackrel{\→}{u}}^{\left(k\right)}\left({\stackrel{\→}{x}}_i\right)-{\stackrel{\→}{u}}^{\left(k\right)}\left({\stackrel{\→}{x}}_j\right))=e^{-\left|\right|{\stackrel{\‾}{u}}^{\left(k\right)}\left({\stackrel{\‾}{x}}_i\right)-{\stackrel{\‾}{u}}^{\left(k\right)}\left({\stackrel{\‾}{x}}_j\right)|{|}^2/2{\mathrm{\σ}}_L^{\left(k\right)}\left({\stackrel{\‾}{x}}_i\right){\mathrm{\σ}}_L^{\left(k\right)}\left({\stackrel{\‾}{x}}_j\right)},---\left(4\right)$

其中，和分别指代在位置处和处的局部噪声估计。

也能够在被认为实际上是同质区域的区域内通过对标准差的确定来人工地获得噪声估计。也能够通过开始于原始数据中的噪声的经典噪声传播来获得空间解析的噪声估计，所述经典噪声传播已知至少为良好的近似。已知的经典噪声传播例如被公开在Adam Wunderlich等人的文章“Imagecovariance and lesion detectability in direct fan-beam x-ray computedtomography”，Institute of Physics and Engineering in Medicine，第2471至2493页(2008年)中，通过引用将所述文章整体并入本文。相同的原理能够应用于微分相位对比计算机断层摄影以及能谱计算机断层摄影。

代替于局部噪声估计，全局噪声估计也能够被使用，尤其是如果噪声不在图像上大地变化。

上面提到的方程式确保在以下情况下第一区域图像元素的图像值的权重减小：如果a)第一区域图像元素与第一图像元素之间的第一相似度和b)第二区域图像元素与第二图像元素之间的第二相似度中的一个减小，并且如果a)第一区域图像元素与第一图像元素之间的第一相似度和b)第二区域图像元素与第二图像元素之间的第二相似度中的另一个保持不变或者减小。

相似度提供两个图像值或者两个图像片块多么可能相同的估计，尤其是如果各自的图像将被去噪的话。因此定义在方程式(4)中的相似度也能够被认为是图像值相似的概率的种类。对于比噪声小得多的差异，可能性或概率接近于一，并且对于大于噪声水平的差异的相似度落到零。尤其是，如果具有大于1x1的尺寸的图像片块被考虑，因为相似度的高斯原型，如果围绕和的图像片块在N图像中的任一个中在统计上不同，则在位置处的图像值的得到的权重将是小的。因此能够确保，在第一图像中的统计学显著的边上没有平滑被执行。滤波因此是边缘保持的。

距离函数，其也能够被认为是空间加权函数，能够由以下方程式定义：

$w_d({\stackrel{\→}{x}}_i-{\stackrel{\→}{x}}_j)=e^{-|{\stackrel{\→}{x}}_i-{\stackrel{\→}{x}}_j|/2{\mathrm{\σ}}_{w_d}^2},---\left(5\right)$

或者

其中，参数控制在平均期间考虑的近邻的有效尺寸。

在下文中，用于如本文提出的联合双边滤波的另外的改进的方法被概述。双边滤波基本上包含两项：距离加权(在上面的方程式(5)和(6)中被称为“w_d”)和范围加权(在上面的方程式(4)中被称为“L”)。距离加权是纯粹几何促动的，并且相比于上面描述的已知方法并不被改变。范围滤波引入基于图像值的加权，所述加权被提出为相比于上面描述的已知方法根据本发明而被改变。

现在应当针对两幅图像说明范例。双边滤波的基本任务是估计两个不同像素属于相同问题类的概率，即底层真实像素值相同的可能性。该估计是基于对四个(或更多个)图像样本的评估的，所述样本是在两幅图像中的相同两个不同位置处的样本。这些样本被指代为a₁、a₂、b₁、b₂(a_i和b_i对应于上面提及的和并且这些用在下文中以便实现更紧凑的标记)，并且它们形成向量(a₁、a₂、b₁、b₂)。在这些样本上的二阶统计通过(新引入的)协方差矩阵C来描述这些测量结果中的噪声。

$C=\left(\begin{array}{cccc}{c}_{11}& c_{12}& c_{13}& c_{14}\\ c_{12}& c_{22}& c_{23}& c_{24}\\ c_{13}& c_{23}& c_{33}& c_{34}\\ c_{14}& c_{24}& c_{34}& c_{44}\end{array}\right)---\left(7\right)$

其中，对角元素描述样本的噪声方差，并且非对角元素描述在图像内的噪声的相关性(c₁₂、c₃₄)或在图像之间的噪声的相关性(c₁₃、c₁₄、c₂₃、c₂₄)。矩阵C是对称并且正定的。利用该协方差矩阵，能够写出在测量的向量(a₁、a₂、b₁、b₂)中获得噪声实现(n₁、n₂、n₃、n₄)的可能性：

$p(n_1,n_2,n_3,n_4)\∝\exp \left(\frac{1}{2}\right(n_1,n_2,n_3,n_4)C^{-1}{(n_1,n_2,n_3,n_4)}^T)---\left(8\right)$

不能够在测量值中区分噪声成分与信号成分。因此，检查在a₁＝a₂和b₁＝b₂的假设下多么可能观测实际测量值(a₁、a₂、b₁、b₂)。针对给定的统计这能够容易地被回答，因为在该假设下，能够通过计算差异Δ_a＝a₁–a₂，Δ_b＝b₁–b₂并且在产生这些差异的(n₁、n₂、n₃、n₄)的所有组合上进行积分消除测量值中的信号成分。为了标记的紧凑，设置S＝C^-1。能够根据以下计算所期望的可能性：

$p({\mathrm{\Δ}}_a,{\mathrm{\Δ}}_b|a_1=a_2,b_1=b_2)\∝\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\exp \left(\frac{1}{2}\left(n_1,{\mathrm{\Δ}}_a-n_1,n_3,{\mathrm{\Δ}}_b-n_3\right)S{\left(n_1,{\mathrm{\Δ}}_a-n_1,n_3,{\mathrm{\Δ}}_b-n_3\right)}^T\right){\mathrm{dn}}_1{\mathrm{dn}}_3---\left(9\right)$

该边缘可能性需要被插入为联合双边滤波方程式的联合可能性项，实现：

针对最一般的协方差矩阵，对积分的明确计算变得相当复杂。然而，一旦协方差矩阵C数值上已知，计算转置矩阵S以及针对实际测量值的可能性就成为简单任务。

两个简单情况应当被考虑：如果四个样本的噪声不相关，那么协方差矩阵及其转置二者都是对角的，并且期望的可能性仅是高斯函数的乘积。另一感兴趣情况是以下形式的协方差矩阵：

$C=\left(\begin{array}{cccc}{c}_{11}& c_{12}& c_{13}& c_{14}\\ c_{12}& c_{11}& c_{14}& c_{13}\\ c_{13}& c_{14}& c_{33}& c_{34}\\ c_{14}& c_{13}& c_{34}& c_{33}\end{array}\right)---\left(11\right)$

几何地，该协方差对其中在两个图像位置处的图像方差相同(即，c₁₁＝c₂₂，c₃₃＝c₄₄)并且在图像之间的协方差关于两个图像位置对称(c₁₃＝c₂₄，c₁₄＝c₂₃)的情况进行建模。在典型的情况下，这些关系近似地被满足，因为在重建的CT图像中的方差在空间中仅缓慢地变化。

该增加的对称也被反映在转置矩阵S中，其读出为：

$S=\frac{1}{N}\left(\begin{array}{cccc}{c}_{34}{H}_6+c_{11}{H}_2-c_{33}{H}_3& c_{33}{H}_6-c_{12}{H}_2-c_{34}{H}_3& c_{13}(c_{13}^2-H_5)+c_{14}{H}_1& c_{14}(c_{14}^2-H_5)+c_{13}{H}_1\\ c_{33}{H}_6-c_{12}{H}_2-c_{34}{H}_3& c_{34}{H}_6+c_{11}{H}_2-c_{33}{H}_3& c_{14}(c_{14}^2-H_5)+c_{13}{H}_1& c_{13}(c_{13}^2-H_5)+c_{14}{H}_1\\ c_{13}(c_{13}^2-H_5)+c_{14}{H}_1& c_{14}(c_{14}^2-H_5)+c_{13}{H}_1& c_{12}{H}_6+c_{33}{H}_4-c_{11}{H}_3& c_{11}{H}_6-c_{34}{H}_4-c_{12}{H}_3\\ c_{14}(c_{14}^2-H_5)+c_{13}{H}_1& c_{13}(c_{13}^2-H_5)+c_{14}{H}_1& c_{11}{H}_6-c_{34}{H}_4-c_{12}{H}_3& c_{12}{H}_6+c_{33}{H}_4-c_{11}{H}_3\end{array}\right)$

其中，

N＝det(C)，

H₁＝c₁₂c₃₃+c₁₁c₃₄-c₁₄c₁₃，

$H_2=c_{33}^2-c_{34}^2,$

$H_3=c_{14}^2+c_{13}^2,$

$H_4=c_{11}^2-c_{12}^2,$

H₅＝c₁₁c₃₃+c₁₂c₃₄，

H₆＝2c₁₃c₁₄。

显而易见，协方差矩阵的对称性也存在于其转置中，所述转置因此能够被写为：

$S=\left(\begin{array}{cccc}{s}_{11}& s_{12}& s_{13}& s_{14}\\ s_{12}& s_{11}& s_{14}& s_{13}\\ s_{13}& s_{14}& s_{33}& s_{34}\\ s_{14}& s_{13}& s_{34}& s_{33}\end{array}\right)---\left(13\right)$

在该假设下，用于对向量(a₁、a₂、b₁、b₂)进行测量的期望的可能性是：

$p({\mathrm{\Δ}}_a,{\mathrm{\Δ}}_b|a_1=a_2,b_1=b_2)\∝F\exp \left(\frac{1}{4}\right({\mathrm{\Δ}}_a^2(s_{12}-s_{11})+2{\mathrm{\Δ}}_a{\mathrm{\Δ}}_b(s_{14}-s_{13})+{\mathrm{\Δ}}_b^2(s_{34}-s_{33})\left)\right)---\left(14\right)$

其中，

$F=\frac{\sqrt{\mathrm{\π}(s_{11}+s_{12})}}{\sqrt{(s_{11}+s_{12})(s_{33}+s_{34})}-{(s_{13}+s_{14})}^2}---\left(15\right)$

应注意，在Eq.(14)中给出为被插入在滤波器Eq.(16)中的联合可能性不具有如现有技术中实现的在两幅图像上导出的独立相似度的乘积的形式，但是有额外的项，即差异的乘积Δ_aΔ_b。

总而言之，提出的在两幅图像a和b上的滤波能够被写出为：

其中，

$N_i=\underset{j\∈N_i}{\Σ}{w}_d(i,j)p(a_i-a_j,b_i-b_j|a_i=a_j,b_i=b_j)---\left(17\right)$

作为体素特异的归一化因子，p是方程式(14)的可能性，w_d(i,j)是双边滤波中的常用距离加权(在方程式(5)和(6)中被称为两个求和在体素i的近邻中的全部体素上进行。第二图像b相应地被滤波。

总而言之，对图像的已知双边滤波将体素特异的权重用于局部平均：

$N_i=\underset{j\∈N_i}{\Σ}{w}_d(i,j)e^{-{(a_i-a_j)}^2/2{\mathrm{\σ}}_n^2}---\left(18\right)$

其中，w_d是距离依赖性权重(例如，高斯)，σ是图像a中的(局部)噪声水平。

对图像的已知联合双边滤波将体素特异的权重用于局部平均：

$N_i=\underset{j\∈N_i}{\Σ}w(i,j)e^{-{(a_i-a_j)}^2/2{\mathrm{\σ}}_a^2}{e}^{-{(b_i-b_j)}^2/2{\mathrm{\σ}}_b^2}$

其中，特意地将具有噪声水平σ_a和σ_b的简单高斯相似度测量用作如已经在Eq.(4)中提及的相似度测量L₁和L₂，以清楚示出针对本发明的不同。

相比之下，提出的对图像的双边滤波使用与来自所有图像的边有关的信息以及完全二阶噪声信息来获得：

$N_i=\underset{j\∈N_i}{\Σ}{w}_d(i,j)p(a_i-a_j,b_i-b_j|a_i=a_j,b_i=b_j)$

$p(a_i-a_j,b_i-b_j|a_i=a_j,b_i=b_j)\∝\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\exp \left(\frac{1}{2}\left(n_1,{\mathrm{\Δ}}_a-n_1,n_3,{\mathrm{\Δ}}_b-n_3\right)S_{ij}{\left(n_1,{\mathrm{\Δ}}_a-n_1,n_3,{\mathrm{\Δ}}_b-n_3\right)}^T\right){\mathrm{dn}}_1{\mathrm{dn}}_3$

其中，Δ_a和Δ_b分别是a_i-a_j和b_i-b_j的简易形式，并且S_ij是针对图像样本(a_i、a_j、b_i、b_j)的噪声协方差矩阵的转置。

尽管已经在附图和前文描述中详细图示并描述了本发明，但是这样的图示和描述应被视为说明性或示范性的，而非限制性的；本发明不限于所公开的实施例。本领域技术人员通过研究附图、公开内容以及权利要求书，在实践要求保护的本发明时，能够理解并实现对所公开的实施例的其他变型。

在权利要求书中，词语“包括”不排除其他元件或步骤，并且词语“一”或“一个”不排除多个。单个元件或其他单元可以完成权利要求书中所记载的若干个项目的功能。尽管在互不相同的从属权利要求中记载特定措施，但是这并不指示不能有利地使用这些措施的组合。

计算机程序可以存储/分布在适合的非暂态介质上，例如与其他硬件一起被提供或作为其他硬件的部分被提供的光学存储介质或固态介质，但是计算机程序也可以以其他形式分布，例如经由因特网或其他的有线或无线的电信系统分布。

权利要求书中的任何附图标记都不应被解释为对范围的限制。

Claims (13)

1.一种用于对图像进行滤波的图像处理装置，所述装置包括：

-图像输入部(3)，其用于获得相同对象的第一图像和第二图像，所述第一图像和所述第二图像包括多个体素并且通过噪声协方差而互相关，每个体素具有包括信号值以及噪声值的体素值，

-联合双边滤波器(3)，其用于对所述第一图像进行滤波，其中，通过包括相对体素特异的权重的滤波函数对所述第一图像的第一体素进行滤波，所述权重包括获得所述第一图像中的所述第一体素的体素值和第二体素的体素值以及所述第二图像中的第一体素的体素值和第二体素的体素值的可能性，假设所述第一图像的所述第一体素的信号值与所述第一图像的第二体素的信号值相同并且所述第二图像的所述第一体素的信号值与所述第二图像的第二体素的信号值相同，

-图像输出部(4)，其用于提供经滤波的图像。

2.如权利要求1所述的图像处理装置，

其中，所述滤波器(3)被配置为使用可能性，所述可能性包括描述所述噪声协方差的噪声协方差矩阵。

3.如权利要求2所述的图像处理装置，

其中，所述滤波器(3)被配置为使用以下形式的协方差矩阵：

$C=\left(\begin{array}{cccc}{c}_{11}& c_{12}& c_{13}& c_{14}\\ c_{12}& c_{22}& c_{23}& c_{24}\\ c_{13}& c_{23}& c_{33}& c_{34}\\ c_{14}& c_{24}& c_{34}& c_{44}\end{array}\right)$

或者

$C=\left(\begin{array}{cccc}{c}_{11}& c_{12}& c_{13}& c_{14}\\ c_{12}& c_{11}& c_{14}& c_{13}\\ c_{13}& c_{14}& c_{33}& c_{34}\\ c_{14}& c_{13}& c_{34}& c_{33}\end{array}\right).$

4.如权利要求2所述的图像处理装置，

其中，所述滤波器(3)被配置为使用以下形式的可能性：

$p(a_i-a_j,b_i-b_j|a_i=a_j,b_i=b_j)\∝\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\underset{-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\∫}}\exp \left(\frac{1}{2}(n_1,{\mathrm{\Δ}}_a-n_1,n_3,{\mathrm{\Δ}}_b-n_3)S_{ij}{(n_1,{\mathrm{\Δ}}_a-n_1,n_3,{\mathrm{\Δ}}_b-n_3)}^T\right){\mathrm{dn}}_1{\mathrm{dn}}_3$

S是所述协方差矩阵的逆矩阵，

a₁和a₂是所述第一图像的第一体素值和第二体素值，

b₁和b₂是所述第二图像的第一体素值和第二体素值，并且

Δ_a＝a₁-a₂，Δ_b＝b₁-b₂。

5.如权利要求2所述的图像处理装置，

其中，所述滤波器(3)被配置为使用以下形式的可能性：

$p({\mathrm{\Δ}}_a,{\mathrm{\Δ}}_b|a_1=a_2,b_1=b_2)\∝F\exp \left(\frac{1}{4}\right({\mathrm{\Δ}}_a^2(s_{12}-s_{11})+2{\mathrm{\Δ}}_a{\mathrm{\Δ}}_b(s_{14}-s_{13})+{\mathrm{\Δ}}_b^2(s_{34}-s_{33})\left)\right)$

其中，

$F=\frac{\sqrt{\mathrm{\π}(s_{11}+s_{12})}}{\sqrt{(s_{11}+s_{12})(s_{33}+s_{34})}-{(s_{13}+s_{14})}^2},$

$S=\left(\begin{array}{cccc}{s}_{11}& s_{12}& s_{13}& s_{14}\\ s_{12}& s_{11}& s_{14}& s_{13}\\ s_{13}& s_{14}& s_{33}& s_{34}\\ s_{14}& s_{13}& s_{34}& s_{33}\end{array}\right),$

S是所述协方差矩阵的逆矩阵，

a₁和a₂是所述第一图像的第一体素值和第二体素值，

b₁和b₂是所述第二图像的第一体素值和第二体素值，并且

Δ_a＝a₁-a₂，Δ_b＝b₁-b₂。

6.如权利要求2所述的图像处理装置，

其中，所述第一图像的经滤波的第一体素值被获得为：

$\widetilde{a_i}=\frac{1}{N_i}\underset{j\∈N_i}{\Σ}{a}_j{w}_d(i,j)p(a_i-a_j,b_i-b_j|a_i=a_j,b_i=b_j)$

其中，

$N_i=\underset{j\∈N_i}{\Σ}{w}_d(i,j)p(a_i-a_j,b_i-b_j|a_i=a_j,b_i=b_j)$

并且，

w_d是双边滤波中的距离权重。

7.如权利要求1所述的图像处理装置，

其中，所述图像输入部(2)被配置为提供所述第一图像和所述第二图像，使得它们示出所述对象的不同性质。

8.如权利要求7所述的图像处理装置，

其中，所述图像输入部(2)被配置为提供能谱计算机断层摄影图像，其中，所述第一图像指示所述对象的第一物理性质并且所述第二图像指示所述对象的第二物理性质，所述第一物理性质和所述第二物理性质是彼此不同的。

9.如权利要求7所述的图像处理装置，

其中，所述第一图像和所述第二图像中的一幅是康普顿图像，并且所述第一图像和所述第二图像中的另一幅是光电图像。

10.如权利要求7所述的图像处理装置，

其中，所述第一图像和所述第二图像中的一幅是康普顿图像或光电图像，并且所述第一图像和所述第二图像中的另一幅是康普顿图像与光电图像的组合。

11.如权利要求1所述的图像处理装置，其中，所述图像输入部(2)被配置为提供微分相位对比计算机断层摄影图像，其中，所述第一图像和所述第二图像中的一幅指示所述对象的吸收，并且所述第一图像和所述第二图像中的另一幅指示由所述对象引起的辐射的相位移位，所述辐射穿过所述对象。

12.一种对图像进行滤波的图像处理方法，所述方法包括：

-获得相同对象的第一图像和第二图像，所述第一图像和所述第二图像包括多个体素并且通过噪声协方差而互相关，每个体素具有包括信号值以及噪声值的体素值，

-对所述第一图像进行联合双边滤波，其中，通过包括相对体素特异的权重的滤波函数来对所述第一图像的第一体素进行滤波，所述权重包括获得所述第一图像中的所述第一体素的体素值和第二体素的体素值以及所述第二图像中的第一体素的体素值和第二体素的体素值的可能性，假设所述第一图像的所述第一体素的信号值与所述第一图像的第二体素的信号值相同并且所述第二图像的所述第一体素的信号值与所述第二图像的第二体素的信号值相同，

-提供经滤波的图像。

13.一种包括程序代码模块的计算机程序，所述程序代码模块用于，当所述计算机程序在计算机上运行时，令所述计算机执行如权利要求12所述的方法的步骤。