专利名称：医学图像中的3d局部结构的鲁棒比例空间分析方法
相关美国申请的交叉引用本申请要求来自2003年7月18日收录的Okada等人的美国临时申请No.60/488,603的“Robust Scale-Space Analysis of 3D Local Structures in medicalImages”的优先权，其内容在此被引入以供参考。
背景技术：
在测定体积的医学图像分析中存在的问题是跨越各种比例来表征肿瘤的3D局部结构，因为肿瘤的大小和形状实际上非常多样化。肿瘤的这种基础比例同样提供了有用的临床信息，与恶性的概率高度关联。有许多以前提出的解决这个问题的方法。然而，当目标形状大大不同于孤立高斯形状时，这些现有的技术方法容易对信号噪音敏感且它们的精确度降低。在医学领域中，由于很多肿瘤呈现为不规则形状这些约束在噪音背景信号中很强。
发明内容
本发明涉及协方差矩阵的鲁棒估计方法，其描述了感兴趣的结构的分布和三维(3D)取向。给定输入信号，利用每个分析带宽的所有可用数据点，为扩展的均值移动向量进行基于均值移动(mean shift-based)的梯度上升。收敛到相同点的数据点被一起分组，形成一组局部结构候选者。这些收敛候选者可以被解释为信号的空间极值。接着，对于每个候选者，基础比例通过约束的最小二乘法估计协方差矩阵而确定。最后，对于每个候选者，跨越分析比例来进行稳定性测试，产生每个局部目标的最佳比例估计。结果，可以找出能够在空间上发生变化的信号的局部比例。
这些方法利用算法，这里被称为基于均值移动的带宽选择算法，用于分析一般离散信号和估计完全参数化的协方差矩阵。用本发明的方法，可以解决表示图像中的局部结构的问题。基于均值移动的带宽选择算法的鲁棒模式和比例选择与完全参数化的协方差矩阵的考虑同样能够以更灵活和鲁棒的方式估计肿瘤的比例，减轻了先前方法的上述缺点。由于医学领域中的许多目标对象处理复杂的3D结构，本发明的方法可以被用于许多不同的应用场合。
附图1是描述本发明的优选实施例的方法的流程图。
附图2是描述本发明的另一个优选实施例的方法的流程图。
附图3-4是描述利用本发明的优选方法来进行分析的结果的图像。
附图5描述了进行本发明的优选方法的计算机系统的优选实施例。
优选实施例的详细说明1、连续信号的固定带宽均值移动向量医学图像可以由在n个d维点xi计算的D维连续信号fRD→R表示，且与每个点xi相关的不确定性可以由D×D矩阵Hi表示，其中i∈1，...n。该矩阵Hi被称为带宽矩阵。该信号可以具有一个或多个极值。信号的极值可以与肿瘤或其它目标对象的位置相关。现在参考附图1，在步骤11，在由f表示的医学图像的分析中的第一步是确定信号的空间极值μ。为了找到该极值，可以首先定义函数，m(x；H)，其中x是对应信号测量的空间位置，H是相应的带宽，这里被称为扩展的均值移动向量，m(x;H)≡∫μΦ(x-μ;H)f(μ)dμ∫Φ(x-μ;H)f(μ)dμ-x·---(3)]]>其中高斯内核Φ(x-μ；H)可以被定义为exp[-12D2(x,μ;H)],]]>其D2(x，μ；H)＝(x-μ)TH-1(x-μ)，且H-1是带宽矩阵的加权调和平均值，H-1(x)=Σi=1nwi(x)Hi.]]>权可以被定义为wi(x)=1|Hi|1/2exp[-12D2(x,μ;Hi)]Σi=1n1|Hi|1/2exp[-12D2(x,μ;Hi)],]]>且可以进行归一化。
方程式(3)可以被如下地对给定固定分析带宽H的f的空间极值μ进行定位。首先，进行极值的估计μ1，然后接着从方程式(3)对这个极值计算m1(x；H)。如果y1被用于表示对于μ1的初始估计的方程式(3)的第一项的结果，那么，对于下一个方程式(3)的迭代，x被替换为y1，μ1被替换为m1(x ；H)，被表示为μ2。这个处理过程可以被重复，每次以前一个迭代的第一项的结果替换方程式(3)的第二项，且根据m(x；H)的前一个估算值估算第一项。对于每个方程式(3)的迭代k，结果误差会收敛到零。其扩展的均值移动向量m(x；H)充分接近零的μk的值可以被视为信号f的极值。信号的数据空间可以通过分组收敛到同样极值的数据点来划分。
2、协方差矩阵的约束的最小二乘解附图1的步骤12，下一步是估计肿瘤的D维分布和取向，其在步骤11中作为空间极值的中心μ被找到。D维局部表面的几何信息可以由在极值估计的协方差矩阵∑表征。
当f可以通过高斯进行近似计算时，协方差∑可以由m(x)≈-H(∑+H)-1(X-μ)。这可以按下面简单的形式改写，∑H-1mi≈bi≡μ-xi-mi(4)考虑收敛到极值μ的所有轨道点{xii＝1，...，tu}，过完备线性方程组可以被表示为A∑≈B， (5)其中A=(m1;...;mtu)TH-T,]]>B=(bi;...;btu)T,]]>且其中在RD×D中∑是对称的正定矩阵。协方差可以由方程式(5)的约束的最小二乘解来估计。这个解产生以下的闭合形式，Σ=UPΣP-1UQ~ΣQ~UQ~TΣP-1UPT;---(6),]]>其中该解涉及以下对称Schur分解ATA=UPΣP2UPT,]]>且BTB≡Q，其中Q~=ΣPUPTQUPΣP=UQ~ΣQ2UQ~T.]]>这个闭合形式可以通过确定区域的唯一最小者而被发现，g(Y)≡||AY-BY-T||F2,]]>其中∑＝YYT。
3、比例选择以上两步可为每个分析带宽H产生一对中心位置和协方差估计{μh；∑h}。下一步的步骤13，涉及寻找跨越带宽范围分析的目标结构的最佳估计值。这个最佳估计值可以通过利用Jensen-Shannon偏差的形式而被发现，JS(h)=12log|12a+1Σj=h-ah+aΣj|Πj=h-ah+a2a+1|Σj|+12Σj=h-ah+a(μj-μ)T[Σj=h-ah+aΣj]-1(μj-μ)]]>
(7)其中μ=12a+1Σh-ah+aμj.]]>给出邻近参数a，可以为每个分析带宽h计算偏差。跨越带宽的偏差JS(h)的极值可以提供最后的比例估计，其在整个比例范围内最稳定。
所述的稳定性测试要求一组分析带宽先验(a prioir)。在本发明的一个实施例中，H＝hI，且h是随恒定步骤变化的。为了获得比例选择的更高性能，优选具有更密集的分布的分析带宽。然而，这样密集的采样会禁止地扩大检索空间，特别是当考虑到完全参数化的H时。
4、局部多比例分析的算法在某些应用场合，所要解决的任务是表示由其它装置提供的粗略位置的局部结构的比例。一个例子是肿瘤的结构分析，其在测定体积图像中的位置由放射学家人工地提供。在这种情况中最简单的策略是仅仅从给定的标记符的点进行均值移动迭代。收敛点充当肿瘤中心估计，且从标记符的所有的轨道点被用于估计比例。当提供的位置被不确定性所污染时和当迭代收敛太快时，这个自然的策略会失效，这使方程式(4)不充分完成且秩不完整。这些问题可以通过以下步骤解决，在附图2中所述。首先，在步骤21，考虑从标记符的近邻区域采样的一组开始点。在步骤22，在进行均值移动迭代后，大多数开始点收敛到的点充当位置估计μ。接着，在步骤23，在估计μ周围的规则的采样可以被进行。在步骤24，比例估计值∑可以通过利用收敛到μ的所有轨道点求解方程式(4)而产生。在步骤25，方程式(7)的相同的稳定性测试可以被用于最后的估计。
6、实例在第5部分中所述的局部多比例分析算法的3D域实现是利用由14位患者显示肺部肿瘤的高分辨度计算机X射线断层造影术(HRCT)图像鉴定的。使用了具有0.25间隔h＝(0.252；...；112)且a＝1的44种分析比例的总和。提供了肿瘤的粗略的位置。作为前处理，利用标记符抽取了32×32×32大小的感兴趣的量。附图3和4示出了得到的中心和部分实心的根瘤(nodule)的例子，其根瘤的几何形状大大偏离简单的高斯结构。对于复杂情况的肿瘤位置、分布和3D取向的正确估计证明了本发明的方法的有效性。
可以理解的是，本发明能够以硬件、软件、固件、特定用途处理器或其结合的各种形式来实现。在一个实施例中，本发明能够以软件的方式被实现为在计算机可读程序存储设备上有形体现的应用程序。该应用程序可以被加载到包括任何适当体系结构的机器，以及通过包括任何适当体系结构的机器进行。
现在参考附图5，根据本发明的实施例，用于实现本发明的计算机系统101可以包括，特别是，中央处理单元(CPU)102、存储器103和输入/输出(I/O)接口104。计算机系统101一般通过I/O接口104被连接到显示器105和各种输入设备106，例如鼠标和键盘。支持电路可以包括电路，例如高速缓冲存储器、电源、时钟电路和通信总线。存储器103可以包括随机存取存储器(RAM)、只读存储器(ROM)、磁盘驱动器、磁带机等，或者它们的结合。本发明可以被实现为例行程序107，其被存储在存储器103中并由CPU102进行来处理来自信号源108的信号。照此，计算机101是通用计算机系统，其在当进行本发明的例程107时成为特定用途计算机系统。
计算机系统101还包括操作系统和微指令代码。这里描述的各种处理和功能可以是微指令代码的一部分或应用程序的一部分(或它们的结合)，其通过操作系统而进行。另外，各种其它的外围设备可以被连接到计算机平台，例如附加数据存储设备和打印设备。
可以进一步理解的是，因为在相应附图中描述的一些组成系统组件和方法步骤能够以软件实现，在系统组件(或处理步骤)之间的实际连接可以根据本发明被编程的方式而不同。这里提供的本发明的给出的教导，相关领域的一个普通技术人员将可以考虑本发明的这些和类似的实现或配置。
尽管已经参考优选实施例详细描述了本发明，本领域普通技术人员可以理解，可以进行各种修改和替代而并不脱离如在所附权利要求中所述的本发明的构思和范围。
权利要求
1.一种分析医学图像中的三维结构的方法，用于确定肿瘤的位置、形状和取向，所述方法包括步骤为一组带宽H寻找D维空间信号f的至少一个空间极值μ；通过为每个极值μ估计信号f的协方差∑来估计每个极值μ附近的信号f的D维分布和取向；以及从至少一个极值μ和协方差∑寻找μ和∑的最稳定的估计值，来从信号f寻找目标结构的最佳估计值。
2.根据权利要求1的方法，其中信号f的空间极值μ通过求解具有第一项和第二项的方程式而找到，上述第一项和第二项是位于该方程式右侧并从左到右计的项，所述方程式是m(x;H)≡∫μΦ(x-μ;H)f(μ)dμ∫Φ(x-μ;H)f(μ)dμ-x]]>其中m(x；H)是扩展的固定带宽均值移动向量，x是对应信号测量的空间位置，且Φ(x-μ；H)是带宽H的高斯内核。
3.根据权利要求2的方法，其中Φ(x-μ；H)可以被定义为exp[-12D2(x,μ;H)],]]>其中D2(x，μ；H)＝(x-μ)TH-1(x-μ)，且H-1是带宽矩阵的加权调和平均值，H-1(x)=Σi=1nwi(x)Hi,]]>且其中权可以被定义为wi(x)=1|Hi|12exp[-12D2(x,μ;Hi)]Σi=1n1|Hi|12exp[-12D2(x,μ;Hi)],]]>且可以进行归一化。
4.根据权利要求2的方法，其中寻找信号f的空间极值μ包括步骤进行极值的估计，并为这个极值估算扩展的固定带宽均值移动向量，其中y1被用来为μ1的初始估计值表示m(x；H)的第一项；估算扩展的固定带宽均值移动向量m(x；H)，以y1替代第二项，且以m(x；H)的前一个计算替代μ1的初始估计值；重复所述扩展的固定带宽均值移动向量m(x；H)的估算，每次以前一个迭代的第一项的替换第二项，且根据m(x；H)的前一个估算值估算第一项，直到μk的值被找到，其扩展的均值移动向量m(x；H)充分接近零；通过对收敛到相同极值的所述数据点进行分组，划分信号的数据点。
5.根据权利要求1的方法，其中估计协方差的步骤包括定义协方差的步骤，通过A∑＝B，其中A=(m1;...;mtu)TH-T]]>是D×D维的正定矩阵，且B=(bi;...;btu)T,]]>其中∑H-1mi≈bi≡μ-xi-mi；并且计算方程式Σ*=UPΣP-1UQ~ΣQ~UQ~TΣP-1UPT;]]>其中ATA=UPΣP2UPT,]]>且BTB≡Q，其Q~=ΣPUPTQUPΣP=UQ~ΣQ~2UQ~T]]>是对称Schur分解。
6.根据权利要求1的方法，其中μ和∑的最稳定的估计值通过为每个带宽值h附近的近邻区域a计算Jenden-Shannon偏差并寻找所述偏差的极值而找到，其通过JS(h)=12log|12a+1Σj=h-ah+aΣj|2a+1Πj=h-ah+a|Σj|+12Σj=h-ah+a(μj-μ)T[Σj=h-ah+aΣj]-1(μj-μ)]]>定义，其中μ=12a+1Σh-ah+aμj·]]>
7.一种分析医学图像中的三维结构的方法，用于确定其位置已经被标记的肿瘤的形状和取向，所述方法包括步骤从标记符的近邻区域采样一组开始点；通过进行平均移动基础梯度上升迭代，为一组带宽H估计极值μ；在极值估计μ的近邻区域中，进行点的规则的采样；通过估计极值μ的点的采样的协方差∑来估计极值μ附近的分布和取向；以及通过计算Jensen-Shannon偏差来估计μ和∑的稳定性。
8.一种确定在医学图像中的肿瘤的位置、形状和取向的方法，所述方法包括步骤通过具有第一项和第二项的方程式，定义扩展的固定带宽均值移动向量m(x；H)，上述第一项和第二项是位于该方程式右侧并从左到右计的项，所述方程式是m(x;H)≡∫μΦ(x-μ;H)f(μ)dμ∫Φ(x-μ;H)f(μ)dμ-x,]]>其中x是对应信号测量的空间位置，且其中Φ(x-μ；H)可以被定义为exp[-12D2(x,μ;H)],]]>其D2(x，μ；H)＝(x-μ)TH-1(x-μ)，且H-1是带宽矩阵的加权调和平均值，H-1(x)=Σi=1nwi(x)Hi,]]>且其中权可以被定义为wi(x)=1|Hi|12exp[-12D2(x,μ;Hi)]Σi=1n1|Hi|12exp[-12D2(x,μ;Hi)],]]>且可以进行归一化Φ(x-μ；H)；进行极值的估计，并为这个极值估算扩展的固定带宽均值移动向量，其中y1被用来为μ1的初始估计值表示m(x；H)的第一项；估算扩展的固定带宽均值移动向量m(x；H)，以y1替代第二项，且以m(x；H)的前一个估计值替代μ1的初始估计值；重复所述扩展的固定带宽均值移动向量m(x；H)的计算，每次以前一个迭代的第一项的替换第二项，且根据m(x；H)的前一个估算值估算第一项，直到μk的值被找到，其扩展的均值移动向量m(x；H)充分接近零；通过对收敛到相同极值的所述数据点进行分组，划分信号的数据点；通过定义协方差为A∑＝B，估计信号的D维分布和取向，其中A=(m1;...;mtu)TH-T]]>是D×D维的正定矩阵，且B=(bi;...;btu)T,]]>其中∑H-1mi≈bi≡μ-xi-mi；并且估算方程式Σ*=UPΣP-1UQ~ΣQ~UQ~TΣP-1UPT;]]>其中ATA=UPΣP2UPT,]]>且BTB≡Q，其Q~=ΣPUPTQUPΣP=UQ~ΣQ~2UQ~T]]>是对称Schur分解；以及从至少一个极值μ和协方差∑中查找μ和∑的最稳定的估计值，通过为每个带宽值h附近的近邻区域a计算Jenden-Shannon偏差并寻找所述偏差的极值，其通过JS(h)=12log|12a+1Σj=h-ah+aΣj|Πj=h-ah+a2a+1|Σj|+12Σj=h-ah+a(μj-μ)T[Σj=h-ah+aΣj]-1(μj-μ)]]>定义，其中μ=12a+1Σh-ah+aμj,]]>来从信号f寻找目标结构的最佳估计值。
9.一种通过计算机可读的程序存储设备，其有形地体现了通过计算机可进行的指令的程序，来进行分析医学图像中的三维结构的方法的步骤，为了确定肿瘤的位置、形状和取向，所述方法包括步骤为一组带宽H寻找D维空间信号f的至少一个空间极值μ；通过为每个极值μ估计信号f的协方差∑来估计每个极值μ附近的信号的D维分布和取向；以及从至少一个极值μ和协方差∑来寻找μ和∑的最稳定的估计值，以从信号f寻找目标结构的最佳估计值。
10.根据权利要求9的计算机可读程序存储设备，该方法步骤进一步包括通过求解具有第一项和第二项的方程式寻找信号f的空间极值μ，上述第一项和第二项是位于该方程式右侧并从左到右计的项，所述方程式是m(x;H)≡∫μΦ(x-μ;H)f(μ)dμ∫Φ(x-μ;H)f(μ)dμ-x,]]>其中m(x ；H)是扩展的固定带宽均值移动向量，x是对应信号测量的空间位置，且Φ(x-μ；H)是带宽H的高斯内核。
11.根据权利要求10的计算机可读程序存储设备，其中Φ(x-μ；H)可以被定义为exp[-12D2(x,μ;H)],]]>其中D2(x，μ；H)＝(x-μ)TH-1(x-μ)，且H-1是带宽矩阵的加权调和平均值，H-1(x)=Σi=1nwi(x)Hi,]]>且其中权可以被定义为wi(x)=1|Hi|12exp[-12D2(x,μ;Hi)]Σi=1n1|Hi|12exp[-12D2(x,μ;Hi)],]]>并且可以进行归一化。
12.根据权利要求10的计算机可读程序存储设备，寻找信号f的空间极值μ的方法步骤包括进行极值的估计，并为这个极值估算扩展的固定带宽均值移动向量，其中y1被用来为μ1的初始估计值表示m(x；H)的第一项；估算扩展的固定带宽均值移动向量m(x；H)，以y1替代第二项，且以m(x；H)的前一个计算替代μ1的初始估计值；重复所述扩展的固定带宽均值移动向量m(x；H)的估算，每次以前一个迭代的第一项替换第二项，且根据m(x；H)的前一个估算值估算第一项，直到μk的值被找到，其扩展的均值移动向量m(x；H)充分接近零；和通过对收敛到相同极值的所述数据点进行分组，划分信号的数据点。
13.根据权利要求9的计算机可读程序存储设备，其中估计协方差的方法步骤包括定义协方差，通过A∑＝B，其中A=(m1;...;mtu)TH-T]]>是D×D维的正定矩阵，且B=(bi;...;btu)T,]]>其中∑H-1mi≈bi≡μ-xi-mi；并且计算方程式Σ*=UPΣP-1UQ~ΣQ~UQ~TΣP-1UPT;]]>其中ATA=UPΣP2UPT,]]>且BTB≡Q，其Q~=ΣPUPTQUPΣP=UQ~ΣQ~2UQ~T]]>是对称Schur分解。
14.根据权利要求9的计算机可读程序存储设备，其中寻找μ和∑的最稳定的估计值的方法步骤包括为每个带宽值h附近的近邻区域a计算Jenden-Shannon偏差和寻找所述偏差的极值，其通过JS(h)=12log|12a+1Σj=h-ah+aΣj|Πj=h-ah+a|Σj|2a+1+12Σj=h=ah+a(μj-μ)T[Σj=h-ah+aΣj]-1(μj-μ)]]>定义，其中μ=12a+1Σh-ah+aμj·]]>
全文摘要
本发明涉及确定在医学图像中的肿瘤的位置、形状和取向的方法，包括通过为一组带宽H进行(11)基于均值移动的梯度上升迭代，为一组带宽H寻找D维空间信号f的多个空间极值μ，以及接着通过为每个极值μ估计信号f的协方差∑，确定(12)每个极值μ附近的信号的D维分布和取向。μ和∑的最佳估计值通过在全组的μ和∑上进行(13)Jensen－Shannon偏差确定。
文档编号G06T7/00GK1826621SQ200480020708
公开日2006年8月30日 申请日期2004年7月19日 优先权日2003年7月18日
发明者K·奥卡达, D·科马尼丘, L·博戈尼 申请人:西门子共同研究公司, 美国西门子医疗解决公司