专利名称：用于医学图像配准的通用变换滤波合成的制作方法
技术领域：
本发明一般涉及图像配准。更特别地，本发明提出一种有效的配准技术，用于高精确性、计算效率性和可靠性地匹配数字图像，特别是医学图像。
背景技术：
图像配准是对从不同的视角和/或通过不同的传感器在不同时刻拍摄的相同景象的两个或多个图像进行叠加的过程。它几何学上使至少两个图像的内容对准，该两个图像分别称为参考图像和传感图像。常规配准公式可表述为给定两个形状、一个输入D(在传感图像中)和一个目标S(在参考图像中)，以及一个相异标准，在S中找到和D的任何点关联的最佳空间转换T，并且最小化转换的形状△＝T(D)和目标S之间的相异标准。可以按照形状△和S(基于特征的配准)的轮廓定义相异标准或者在由这些轮廓(基于区域的配准)确定的整个区域中定义相异标准。
图像配准是所有图像分析任务中的关键步骤，其中从各种数据源的组合中获得最终信息。典型地，远程传感、环境监视、天气预报等领域中需要配准。在医学领域，它可用于例如组合计算机X线断层摄影术(CT)和核磁共振(NMR)数据以获得关于病人的更完整的信息，诸如监视肿瘤生长、检验治疗效率、对病人的数据和解剖图进行比较等等。
在《图像和视觉计算》(Image and Vision Computing)21(2003)第977-1000页，由Barbara Zitova和Jan Flusser撰写的“Image registration methodsa survey”中，回顾了大量的配准方法。大多数方法包括以下阶段-特征检测手工或优选地自动检测突出的和区别性的对象(封闭边界区域、边、轮廓、交汇点、角等)。此阶段非常重要，因为它能显著地影响所选的配准算法的性能；-特征匹配在传感图像中检测到的特定特征与在参考图像中检测到的用于表示该相同特征的象素之间建立对应关系；-转换模型评估评估将传感图像匹配到参考图像的转换，其通常被称为映射函数。该映射函数基于图像域Ω来定义，Ω对于两个图像是通用的。使用预先建立的特征对应关系来计算映射函数的参数。此阶段可能需要大量的计算时间，特别是当将配准的图像中的内容具有熟练的非刚性运动时。这例如可以联系到图像中的对象的变形(与简单转换或旋转相对的)，和/或联系到不同的视角。可使用相异函数的最优化来评估映射函数，相异函数还称为目标函数；-图像再取样和转换通过使用适当的内插法技术的映射函数转换传感图像。
在一些方法中，第二和第三阶段并入第一步骤中。
在《计算机视觉和图像理解》(Computer Vision and Image Understanding)第142-165页，2003，由N.Paragios等撰写的“Non rigid registration using distancefunctions(使用距离功能的非刚性配准)”中，揭示了为了改进配准，在并入的特征匹配和转换模型评估阶段之前将距离转换应用到所检测的特征。该距离转换，也称为距离映射DM(x，y)，定义为Φc(x,y)=0(x,y)∈C+DM((x,y),C)(x,y)∈Rc-DM((x,y),C)(x,y)∈Ω-Rc]]>其中C是图像域Ω中的给定特征，Rc是C的凸包，而DM((x，y)，C))是在点阵位置(x，y)和特征C之间的最小欧几里德(Euclidean)距离。
当使用梯度下降方法来最小化相异标准时，距离转换提供了一适宜的特征空间，其允许大的捕捉范围(其中可以比较相似特征的距离)和定位精确性。
不仅如此，尽管很有前景，但是使用距离转换也具有一些缺点。特别地，它需要从图像中提取特征(第一步骤)，这不是总能以一种可靠的方式进行的。而且，距离映射(功能ED)的提取是非线性图像处理操作，众所周知其对于噪声是敏感的。
发明内容
本发明的一个目的是提供一种方法，在保持大捕捉范围和定位精确性的优点的同时避免距离转换的缺点。
相应的，本发明提供一种根据权利要求1的装置、根据权利要求11的方法和根据权利要求17的计算机程序产品。
本发明利用了特定滤波核心(filter kemel)类型，其在保持配准方法的大捕捉范围的同时保证了定位精确性。这种滤波核心表现为在其中心周围有尖锐的峰值，并且实质上具有随着与核心的原点之间距离的增加有一个指数式衰减或幂反比定律(inverse power law)的性质。使用这种滤波核心过滤传感图象和参考图象使得在保持被配准在一起的特征的细节和使它们充分地变模糊以允许大捕捉范围之间提供了一个良好的折衷。
本发明的其它特征和优点将在此后的描述中结合附图进一步展现，其中图1是根据本发明的配准方法的流程图；图2是表示不同滤波核心的图表；图3是示出传感图象特征的卷积的示意图，该传感图象具有本发明中使用的滤波核心；以及图4是示出根据本发明的教导编程的通用计算机的示意图。
具体实施例方式
本发明用于解决两个或多个图象的配准。尽管本发明是以软件实现的方式阐述，它还可通过硬件组件来实施，例如在计算机系统的图象卡中。
在下面的描述中，仅作出将两个图象配准在一起的参考。图像处理领域的技术人员可将本发明轻易地扩展到大量初始化图象的应用。
现在参考附图，特别是参考其中的图1，示出了根据本发明的配准方法的示意性框图。整个方案包括由医学2D或3D图像D或者传感图象的获取组成的初始化步骤10，上述图象将用于与在步骤12中检索的参考图象S进行配准。输入图象S或者是之前获取的或者是例如从数据库中取出的。如果需要的话，初始化步骤10可包括对该图象中的一个或两个进行数字化。
在可选的第二步骤20中，在图像D和参考图象S中检测特征。这就提供了特征增强图象EDI(D)和EDI(S)。所检测到的特征可以是图象中描述的对象的边界(然后，在下面的描述中称EDI(D)和EDI(S)为边界检测图象)。它们还可包括棱线或者诸如例如血管的管状对象的中心线。
使用诸如局部方差方法的本领域已知技术创建特征增强或边界检测图象。源图象须经边界检测以便检测其包含的对象的轮廓。因此，边界检测图像中的像素值说明源图像中感兴趣区(ROI)中的特征。它们表示特征凸出量，其可以是像素强度值、像素强度中的局部梯度或者是源图像中关于特征强度的任何适当的数据。
根据本发明在配准方法中使用的滤波核心，第二步骤保持可选的状态，并且如下所述，第二步骤精确到足以能够避免对提取特征的需要。
在第三步骤30中，将各向同性的低通滤波L施加到边界检测图像EDI(D)和EDI(S)、或传感图像D和参考图像S中。已知配准方法的距离映射由线性卷积代替，已知配准方法的距离映射由非线性操作产生。
该操作需要各向同性滤波核心32。这种核心的一般形状在图2中通过曲线B示出。这种滤波核心必须显示尖锐的中心峰值，并带有向远处相对缓慢的衰减，以便结合定位精确性和清晰度，还带有一大捕捉范围。离开原点，滤波核心可表现为r的指数衰变或幂反比定律，r为到核心中心的距离。
图2中的曲线A示出各向同性的高斯滤波核心，其常用于图像分析。它没有原点处的核心B那样尖锐，而且它还在大的距离衰减得更快。在图2中，曲线A和B描述了各向同性的滤波核心，其峰值均具有相同的有效宽度W=Σrr2.L(r)ΣrL(r).]]>中心锐度比使用的各向同性的滤波核心L(r)缓慢的衰减特征是本发明的一些实施例，并且可如下量化半宽度(即，在r＝±W2)的L(r)的斜率大于具有相同有效宽度W的各向同性高斯核心A的斜率的至少三倍，而两倍宽度(即，在r＝2.W)的L(r)的斜率小于具有相同有效宽度W的各向同性高斯核心A的斜率的至少三倍。
组合了图3(如后所述)的过滤实施的这种尖锐的的滤波核心(图2示出)允许当在大距离上模糊所检测的特征以扩大捕捉范围时聚焦在这些特征上，从而减少了在特征提取中对噪声或错误的敏感度。这种滤波核心引入了与距离映射相对的特征的平滑阈值。
设计了改进的具有将与核心的原点的非零距离(半径r为到核心中心的距离)表示为类似exp(-kr)/rn的核心的各向同性的滤波核心，而不是高斯滤波的传统的exp(-r2/2σ2)表示，n为≥0的整数。这种核心对于与特征的定位刻度s相比的小距离是很明显的，并且根据对于范围从该刻度s到ηs的距离的上述定律应该更平缓，其中η是适应于图像尺寸的参数，典型地η≈10。系数k的值也适应于所需要的定位刻度s。
这种各向同性的滤波核心L(r)可以使用一组带有不同离散核心尺寸σ的高斯曲线，作为高斯滤波(对于d维图像，d是大于1的整数)的连续分布的近似来推导，每个核心被给出一个权重g(σ)。得出的滤波具有一个等于高斯核心的权重和的核心L(r)=Σσg(σ)·e-r2/σ2σd---(1)]]>为了计算效率，使用了对每个分辨率级别都具有一个或多个单一σ高斯曲线(递归无限脉冲响应或IIR)的多分辨率锥体。
实际上，使用上述定义的核心过滤任何图像可通过如下步骤执行首先使用多种不同尺寸σ的标准高斯核心过滤它，然后通过向每个使用尺寸σ的核心过滤的单独的图像给出权重g(σ)来线性地组合得到的多个过滤后的图像。
当应用于边界检测图像EDI(D)和EDI(S)、或传感图像D和参考图像S中时，变量σ的每个高斯滤波首先施加到这些图像中的一个以生成单独的过滤图像，从初始图像中生成多个单独的过滤图像。得到的过滤图像(带有上述定义的核心)从单独的过滤图像和使用权重g(σ)的高斯滤波的加权组合中获得。
其它方法(计算上需要更多代价的)可用于这种滤波合成(例如，傅里叶域，基于解决适当的部分差分等式，等等)。
图3中示出了这种滤波核心如何施加到边界检测图像。已经确定被配准的对象281(从传感图像D或参考图像S中)的轮廓280。为了从边界检测图像L(EDI(D(p)))和L(EDI(S(p)))中计算过滤图像F(p)，在像素p周围定义了窗口win(p)(此处窗口是圆形的，而p是它的中心)，也为win(p)内部的所有像素j定义了各向同性空间分布Wj(p)。在p点该空间分布是最大值，并且对于属于以p为中心的圆周的所有j像素都是相同的值。在win(p)以外，空间分布是空值。关于win(p)的空间分布与等式(1)的滤波核心相应。win(p)的窗口大小是为滤波核心选择的参数的函数。因此我们有F(p)Σj∈win(p)ED(j).Wj(p)---(2)]]>其中根据(1)Wj(p)=L(r),]]>r是在像素j和p之间的欧几里德距离，ED(p)是在边界检测图像EDI(D)或EDI(S)中像素p的像素值。一般来说，对于不属于边界的所有像素j，ED(j)＝0。
当直接将滤波核心L施加到初始化参考图像和传感图像时，相同的阐述保持有效，即，略去在前的特征提取。
在配准方法的第四步骤40，确定映射函数T。此步骤中可使用确定映射函数的本领域已知技术。例如，可使用前面提到的由Paragios等编写的文章中描述的映射函数。
选择有效的特征空间之后，即，使用距离映射(本发明中其由使用等式(1)的特定滤波核心的卷积替代)之后，该映射函数确定包括全局线性配准模型(刚性、仿射，等等)和定位变形的集合。
如上所述，配准实际上是寻找传感图像D和参考图像S之间的空间变换，其对于图像域Ω最小化给定的相异准则。在Paragios等的文章中，对于Ω的被最小化的相异准则是Φ=Σp∈Ω[D(p)-S(p^)]2---(3)]]>其中p^=A(p)+U(p)=T(p);]]>A是说明D和S之间刚性变形的仿射变换，并包括旋转、平移和比例系数；U是说明定位非刚性变形的变换并且具有仅用于被配准的形状的非刚性部分的值。
T是映射函数。
诸如梯度下降方法的迭代处理用于恢复最优配准参数。
这种方法的改进包括引入缓慢变换的权重α、β和γ，其中α2+β2＝1，以便最小化下述改进后的相异准则的归一化的卷积的值Φ′=Σp∈Ω[α(p).D(p)-β(p^).S(p^)+γ]2---(4)]]>其中α(p)和β(p)是缓慢变换的权重，说明图像D和S之间的不同比例系数。这些权重是归一化的，即α2+β2＝1。它们的确定使用如此后描述的具有相对大的核心的低通滤波LP；γ是独立于P的偏移值。
可以以传统的方式通过D(p)和S(T(p))的阶(order)1和2的统计矩来进行权重参数α和β的最优选择。为了允许这些权重参数的一定量的空间变化，建议通过引进以每个像素P为中心的窗口winLP(p)来局部地计算D(p)和S(T(p))的阶(order)1和2的统计矩，在该窗口winLP(p)上定义了相应于上述低通滤波核心LP的空间分布WLP(u)⟨M1(p)⟩=Σu∈winLP(p)M1(p+u).WLP(u),---(5)]]>⟨M12(p)⟩=Σu∈winLP(p)M12(p+u).WLP(u),---(6)]]>⟨M1.M2(p)⟩=Σu∈winLP(p)M1(p+u).M2(p+u)WLP(u)---(7)]]>其中u是窗口winLP(p)中得到的像素；而M1(p)是D(p)或者S(T(p))，而M2是此两者中另一个。
核心LP的大小足够说明图像描述的场景的局部失真的有限范围。它例如可以是高斯滤波，基于允许的失真选择其变化。
诸如梯度下降方法的迭代处理用于恢复最优配准参数。它包含Φ′的最小值的每一迭代步骤的计算。在初始化步骤中，α＝β＝1/而γ＝0。在每次循环中更新α、β和γ的值。
使用了图3中示出的滤波核心，结合低通滤波核心LP的使用，根据本发明的分割方法实际上基于灰度级别变形(边界强度数据的使用)和几何变形(缓慢变化的权重α、β和γ)的连续使用。
在图1的第五步骤50中，使用本领域已知技术形成了合成数据图像。它包括参考图像和传感图像，通过使用了适当的内插技术的映射函数T的方法对图像进行变换。
本发明还提供了一种装置，用于配准包括至少两维的像素数据组的图像，还包括用于接收输入图像D的获取装置和用于存储参考图像S的存储装置，无论是提前获取还是从数据库中获取，最优特征检测可以从输入图像和参考图像中提供例如边界检测图像。本发明的装置进一步包括可以实施上述方法的处理。
本发明可使用根据本发明的教导进行了编程的传统的通用数字计算机或微处理器来方便地实施。
图4是根据本发明的计算机系统300的框图。计算机系统300可以包括CPU(中央处理器)310、存储器320、输入设备330、输入/输出传输通道340和显示设备350。其他设备，如附加的磁盘驱动器、附加的存储器、网络连接…可以包括在该系统中，但是此处没有显示。
存储器320包括包含了被配准的传感图像和参考图像的数据文件。存储器320可进一步包括由CPU310执行的计算机程序产品。该程序包括根据本发明执行上述方法来配准图像的指令。输入设备用于从用户那里接收指令，例如无论是否提供边界检测图像。输入/输出通道可以用于从外部传感器装置接收将要存储于存储器320中的传感图像D，以及向其他装置发送配准图像(输出图像)。显示设备用于显示包括从传感图像和参考图像生成的配准图像的输出图像。
权利要求
1.一种用于配准包括至少是二维的数据组的图像的装置，该装置包括至少获得第一图像(10)和第二图像(12)的装置；用于使用低通滤波核心(32)对所述第一和所述第二图像进行滤波的装置，该低通滤波核心(32)具有尖锐的中心峰值和远离所述中心峰值的缓慢衰减；和用于确定将所述滤波后的第一图像变换到所述滤波后的第二图像的映射函数(40)的装置。
2.如权利要求1的装置，其中所述低通滤波核心是各向同性的并且具有比各向同性高斯滤波核心基本更尖锐的峰值，该高斯滤波核心具有与所述低通滤波核心相同的有效宽度，而且所述低通滤波核心和高斯滤波核心相比在有效宽度之外减少得基本更缓慢。
3.如权利要求2的装置，其中在离开原点W/2的距离，所述低通滤波核心具有的斜率大于所述高斯核心的斜率的至少三倍，和在离开原点2.W的距离，所述低通滤波核心具有的斜率小于所述高斯核心的斜率的至少三倍，其中W表示所述核心的有效宽度。
4.如上述任一权利要求的装置，其中所述低通滤波核心基本上具有离开中心峰值的指数衰减或幂反比定律的性质。
5.如上述任一权利要求的装置，其中所述低通滤波核心是各向同性滤波核心，该各向同性滤波核心被定义为不同变量σ的高斯滤波的和L(r)=Σσg(σ)·e-r2/σ2σd,]]>d是所述第一和第二图像的维数，r是离开所述滤波核心中心的距离，并且每个变量σ的高斯滤波被分配一权重g(σ)。
6.如上述任一权利要求的装置，其中所述每个变量σ的高斯滤波首先分别施加到所述第一和第二图像以生成第一和第二批单独的过滤图像，并且所述第一和第二过滤图像分别从具有使用所述权重g(σ)的所述高斯滤波的所述第一和第二批单独的过滤图像的加权组合中获得。
7.如上述任一权利要求的装置，进一步包括用于将确定的映射函数施加到所述第一图像以产生与所述第二图像进行配准的图像的装置。
8.如上述任一权利要求的装置，其中用于获取第一和第二图像的装置包括用于处理两个输入图像以提供作为基于该两个输入图像的特征增强图像的所述第一和第二图像的装置。
9.如上述任一权利要求的装置，其中确定映射函数的装置包括用于最小化相异准则的装置，该相异准则通过具有缓慢变化的权重的变换的过滤第一图像和过滤第二图像之间的加权差的平方来定义。
10.一种医学成像系统，其包括用于获取至少两个表示身体器官的输入图像(10，12)的装置，和如上述任一权利要求所述的用于配准所述图像的装置。
11.一种配准包含至少是二维的数据组的图像的方法，该方法包括以下步骤至少获得第一图像(10)和第二图像(12)；使用低通滤波核心(32)对所述第一和所述第二图像进行滤波，该低通滤波核心(32)具有尖锐的中心峰值和远离所述中心峰值的缓慢衰减；和确定将所述滤波后的第一图像变换到所述滤波后的第二图像的映射函数(40)。
12.如权利要求10的方法，其中所述低通滤波核心是各向同性的并且具有比各向同性高斯滤波核心基本上更尖锐的峰值，该高斯滤波核心具有与所述低通滤波核心相同的有效宽度，而且所述低通滤波核心和高斯滤波核心相比在有效宽度之外减少得基本更缓慢。
13.如权利要求11的方法，其中在离开原点W/2的距离，所述低通滤波核心具有的斜率大于所述高斯核心的斜率的至少三倍，和在离开原点2.W的距离，所述低通滤波核心具有的斜率小于所述高斯核心的斜率的至少三倍，其中W表示所述核心的有效宽度。
14.如权利要求11-13之一的方法，其中所述低通滤波核心基本上具有离开中心峰值的指数衰减或幂反比定律的性质。
15.如权利要求11-14之一的方法，其中所述低通滤波核心是各向同性滤波核心，该各向同性滤波核心被定义为不同变量σ的高斯滤波的和L(r)=Σσg(σ)·e-r2/σ2σd,]]>d是所述第一和第二图像的维数，r是离开所述滤波核心中心的距离，并且每个变量σ的高斯滤波被分配一权重g(σ)。
16.如权利要求11-15之一的方法，其中确定映射函数的步骤包括最小化相异准则，该相异准则通过具有缓慢变化的权重的变换的过滤第一图像和过滤第二图像之间的加权差的平方来定义。
17.一种计算机程序产品，用于在图像处理装置的处理单元中执行，该计算机程序产品包括当在该处理单元中运行该程序产品时执行根据权利要求11-16之一的分割方法的指令。
全文摘要
使用低通滤波核心过滤将被配准的第一和第二图像，该低通滤波核心具有尖锐的中心峰值和远离该中心峰值的缓慢衰减。该装置确定将过滤的第一图像转变成过滤的第二图像的映射函数。
文档编号G06T7/00GK101052993SQ200580037680
公开日2007年10月10日 申请日期2005年10月25日 优先权日2004年11月5日
发明者S·马克拉姆－埃贝德 申请人:皇家飞利浦电子股份有限公司