CN108701360B - 图像处理系统和方法 - Google Patents

图像处理系统和方法 Download PDF

Info

Publication number
CN108701360B
CN108701360B CN201680081685.4A CN201680081685A CN108701360B CN 108701360 B CN108701360 B CN 108701360B CN 201680081685 A CN201680081685 A CN 201680081685A CN 108701360 B CN108701360 B CN 108701360B
Authority
CN
China
Prior art keywords
image
registration
function
volume
iteration
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201680081685.4A
Other languages
English (en)
Other versions
CN108701360A (zh
Inventor
S·卡布斯
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Koninklijke Philips NV
Original Assignee
Koninklijke Philips NV
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Koninklijke Philips NV filed Critical Koninklijke Philips NV
Publication of CN108701360A publication Critical patent/CN108701360A/zh
Application granted granted Critical
Publication of CN108701360B publication Critical patent/CN108701360B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G06T7/337Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving reference images or patches
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/149Segmentation; Edge detection involving deformable models, e.g. active contour models
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30061Lung

Abstract

用于迭代地计算图像配准或图像分割的系统和方法。配准和分割计算由包括相似性量度分量的优化函数驱动,基于在迭代期间监测图像位置处的体积元素的体积变化来相对减轻所述相似性量度分量对所述迭代计算的影响。还有用于量化配准误差的系统和相关方法5。这包括将一系列边缘检测器应用于输入图像并将相关的滤波器响应组合成组合响应。该一系列滤波器使用滤波器参数进行参数化。然后找到组合响应的极值,并且然后返回与所述极值相关联的滤波器参数作为输出。该滤波器参数涉及给定图像位置10处的配准误差。

Description

图像处理系统和方法
技术领域
本发明涉及图像处理系统,图像处理方法,计算机程序单元和计算机可读介质。
背景技术
在各种临床应用中需要图像配准(也称为图像对齐或图像融合),例如,针对来自不同协议、模态、患者定位、采集时间或对比度状态的图像数据的空间对齐。
在临床工作站中,在察看器功能中采用图像配准。这些功能允许并排查看不同的图像数据集。用户通过鼠标点击或以其他方式选择第一图像数据集中的位置,并且然后第二图像数据集自动滚动到解剖学上对应的位置。
图像配准的其他应用可以在灌注分析系统或时间序列分析系统中找到。在荧光透视成像系列中采集的帧相对于参考帧对齐,从而然后执行灌注参数估计。
再其他的图像配准应用是多模态肝脏分析、心脏分析、前列腺分析和一般的随访分析。这些应用中的大多数需要非刚性或非仿射配准,即,搜索位移的密集矢量场而不是针对刚性或仿射变换的仅六个或十二个参数。
另一个重要的图像处理任务是自动图像分割,其中人们希望找到感兴趣的图像结构,例如器官描绘。
遗憾的是,已经发现输入图像中的特定缺陷会导致配准或分割算法的严重表现不佳,所述分割算法算法在其他情况下己知会产生准确的结果。
而且,关于图像配准,期望一种检测次优精度的机制,以确保临床实践中针对应用的限定水平的准确度。这是因为已经发现准确的图像对齐对于大量后续图像处理步骤或可视化意义重大。
可以通过几种方式来估计图像配准的准确性。定性测量(例如,视觉排序)倾向于比定量测量(例如,标注的标记,功能验证)耗时更少(但不太准确)。
发明内容
因此,可能需要用于图像配准或图像分割的替代系统或方法,或者用于支持配准任务的系统或方法。
本发明的目标通过独立权利要求的主题得以解决,其中,在从属权利要求中并入了另外的实施例。应当注意,以下描述的本发明的方面同样地应用于图像处理方法计算机程序单元以及计算机可读介质。
根据本发明的第一方面,提供了一种图像处理系统,包括:
输入端口,其用于接收两幅或更多幅输入图像以进行配准,包括源图像和目标图像;
迭代求解器,其被配置为通过迭代通过至少一个中间配准变换T(N-j)来求解最终配准变换T(N),其中,所述迭代求解器由包括至少一个函数分量的优化函数驱动,所述函数分量基于测量目标图像与通过将所述至少一个中间配准变换应用于所述源图像能够获得的变换图像之间的相似性的相似性量度,
体积监测器,其被配置为在迭代周期中监测当给定的迭代周期中应用各自的中间配准变换时所述源图像中的图像位置周围的预定义的邻域的体积变化;
体积评估器,其被配置为针对接受条件来评估所述体积变化;以及
校正器,其被配置为当由所述图像评估器发现所述变化违反所述接受条件时,在后续的迭代周期中减轻或至少相对减轻所述第一函数分量对所述图像位置的影响。所述接受条件基于所述图像位置周围的邻域的体积变化。取决于输入图像的维度,本文中使用的“体积”包括3D(三维)中的体积或2D(二维)中的“体积”(通常称为“区域”)。
换句话说,所提出的配准系统包括自校正特征,其中响应于迭代期间个体图像位置处的体积变化,基于目标(或代价)函数自动减轻更新函数的至少一部分。减轻包括优选地完全抑制,但是对于一些实施例,部分抑制可能是足够的。已经发现该自校正特征使得所提出的配准算法对输入图像中编码的伪影或异常更加鲁棒。
根据一个实施例,在迭代的剩余部分期间,所述第一函数分量对于所述图像位置保持减轻。
代价或目标函数通常(但必然)包括多于一个的函数分量。特别地,所述代价函数包括,除了针对相似性的所述函数分量之外,还有用于增强特定期望的图像属性(例如平滑度等)的正则化器分量。但是在一个实施例中,所述校正器被配置为仅减轻针对所述相似性量度的所述函数分量而不减轻所述正则化器。所述正则化器在整个迭代期间保持有效。通过分别乘以<1的适当的因子或因子“0”,减轻可以是仅向下加权或者完全消除。
代价函数的函数分量通常在可从代价函数导出的更新函数中具有对应的功能对应物。更新函数用于迭代,并且在一个实施例中,它是更新函数中的对应函数分量,其由校正器减轻,以便控制迭代中的相似性量度的影响。
与上述“直接”减轻相反,可以通过对正则化器或其他项进行向上加权或放大,同时相似性量度分量保持不变,来实现相对减轻。简而言之,相似性量度分量是向下加权而正则化器保持不受影响,或者正则化器(或至少一个其他项)被放大或向上加权而相似性量度分量保持不受影响。
在迭代期间考虑的变换是非仿射类型的并且尤其包括变形和/或剪切。这些类型的非仿射变换可以由矢量场表示,根据所述矢量场可以推导出局部体积元素变化。体积元素例如是以各图像位置为中心的图像元素的邻域(例如,球体或圆形或任何其他形状)。
根据第二方面,提供了一种图像处理方法,包括:
接收两幅或更多幅输入图像以进行配准,包括源图像和目标图像;
通过迭代通过至少一个中间配准变换T(N-j)来迭代地求解最终配准变换T(N),其中,所述迭代求解器由包括至少一个函数分量的优化函数驱动,所述函数分量基于测量目标图像与通过将所述至少一个中间配准变换应用于所述源图像能够获得的变换图像之间的相似性的相似性量度,
在迭代周期中监测当在给定的迭代周期中应用各自的中间配准变换时所述源图像中的图像位置周围的预定义的邻域的体积变化;
针对接受条件来评估所述体积变化;并且
当图像评估器(IVE)发现所述变化违反所述接受条件时在后续的迭代周期中相对减轻第一函数分量对所述图像位置的影响。
可以在基于模型的分割算法中以类似的方式使用所述自校正特征。换句话说,根据第三方面,提供了一种图像处理系统,包括:
输入端口,其用于接收输入图像以用于针对基于模型的分割;
迭代求解器,其被配置为通过迭代通过初始模型的至少一个中间几何状态来求解最终分割,其中,所述迭代求解器由包括至少一个函数分量的优化函数驱动,所述至少一个函数分量包括基于图像特征量度的第一函数分量,所述图像特征量度在所述迭代期间对所述输入图像中的至少一个图像位置中的至少一个图像特征进行量度,
模型监测器,其被配置为在迭代期间监测模型的几何变化;
几何变化评估器,其被配置为针对接受条件来评估所述几何变化;以及
校正器,其被配置为当由所述几何变化评估器发现所述几何变化违反所述接受条件时,在后续的迭代中相对减轻所述第一函数分量对所述图像位置的影响。
根据一个实施例,所述接受条件基于所述模型的变形。
同样,类似于配准实施例,所述代价函数优选地包括额外的函数分量,所述额外的函数分量对模型的属性进行量度并评估在迭代期间模型如何被改变(例如变形)。同样,在一个实施例中并且为了实现相对减轻,所述校正器将仅减轻图像特征分量并且不减轻对所述模型的属性进行量度的函数分量。同样,减轻可以是仅向下加权或完全消除。相反,图像特征分量保持不变,但现在针对模型属性的函数分量被放大或向上加权。可以通过在迭代的更新函数中减轻或放大相应的函数分量来再次实现减轻动作。
根据第三方面,提供了一种图像处理方法,包括:
接收输入图像以用于基于模型的分割;
通过迭代通过初始模型的至少一个中间几何状态来迭代地求解最终分割,其中,所述迭代求解器由包括至少一个函数分量的优化函数驱动,所述至少一个函数分量基于在迭代期间对所述输入图像中的至少一个图像位置中的至少一个图像特征进行量度的图像特征量度,
在迭代期间监测所述模型的几何变化;
针对接受条件来评估所述几何变化;并且
当由所述几何变化评估器发现所述几何变化违反所述接受条件时在后续的迭代中相对减轻所述第一函数分量对所述图像位置的影响。
在上述实施例中,根据较早的迭代步骤,相对于体积或几何形状来取得体积或几何形状的变化。具体地,针对各图像位置定义初始邻域或几何结构,并且在一个或多个后续的迭代步骤中监测变化。
根据第五方面,提供了一种用于评估配准结果的图像处理系统,包括:
输入端口,其用于接收两幅或更多幅输入图像;
测试模块,其被配置为将一系列参数化边缘检测滤波器应用于所述两幅或更多幅图像,以根据输入图像产生一系列不同的滤波器响应,其中,该系列中的检测滤波器利用滤波器参数来参数化,所述滤波器参数与各个边缘检测滤波器的特性有关;
组合器,其被配置为将来自两幅或更多幅图像的各个响应组合成组合响应,所述组合响应作为所述滤波器参数的函数非单调地变化;
极值查找器,其被配置为查找作为所述滤波器参数的函数所述组合响应的极值;以及
输出端口,其被配置为输出来自滤波器参数中的与所述极值相关联的滤波器参数,如此输出的滤波器参数与图像位置处的配准结果的质量度量有关。
根据一个实施例,所述组合器操作以逐点地将各个响应相乘,以沿着所述图像的至少一个空间维度来形成每图像位置的权重。更具体地并且根据一个实施例,所述组合器操作以对所述权重求和以获得组合响应。
根据一个实施例,所述边缘检测滤波器中的至少一个基于二阶或更高阶图像导数。
更具体地,根据一个实施例,所述滤波器参数涉及二阶或更高阶导数的导数核。更加具体地并且根据一个实施例,所述导数核是基于高斯分布的,并且其中,所述滤波器参数指定所述高斯分布的标准偏差。
已经观察到以此方式构造组合响应(表面或阵列)以产生针对噪声鲁棒的配准准确度,并且能够应对两幅输入图像源于不同成像模态的情况。可以最小化或完全避免假阳性响应。
根据一个实施例,所述系统包括可视化器,所述可视化器被配置为在显示设备上可视化每图像位置的输出质量度量。换句话说,可以生成并输出精度图,使得用户可以确定配准确实产生良好对应性的位置和不能产生良好对应性的位置。通过目视检查可以找到具有相对较高配准的位置。换句话说,本文中提出的系统允许在空间上求解配准错误的量。
根据第六方面,提供了一种
用于评估配准结果的图像处理方法,包括:
接收两幅或更多幅输入图像;
将一系列参数化的边缘检测滤波器应用于所述两幅或更多幅图像,以根据输入图像产生一系列不同的滤波器响应,其中,所述系列中的所述检测滤波器利用滤波器参数来参数化,所述滤波器参数与各个边缘检测滤波器的特性有关;
将来自两幅或更多幅图像的各个响应组合成组合响应,所述组合响应作为所述滤波器参数的函数非单调地变化;
根据所述滤波器参数找到所述组合响应的极值;并且
输出来自所述滤波器参数中的与所述极值相关联的滤波器参数,如此输出的滤波器参数与图像位置处的配准结果的质量度量有关。
所提出的配准精量度化可以应用于真实的患者数据,而不是像过去那样仅仅用于测试数据。而且,所提出的方法产生准确度图,从而提供每图像位置估计。因此,准确度评估可以潜在地基于整个数据集而不是仅基于相对小的一组图像位置或观察平面。简而言之,所提出的方法允许局部的定量准确度估计。准确度估计可以在合适的长度尺寸上局部地表示,例如以mm或cm为单位等。
所提出的用于配准或分割和配准精量度化的方法和系统可用于任何成像模态或成像模态的组合。特别地,所提出的方法和系统可以用于非周期性图像数据集,即,在采集时间上不是周期性的图像数据集。
附图说明
现在将参考附图描述本发明的示范性实施例,其中,
图1示出了图像处理系统的框图;
图2示出了图像处理方法的流程图;
图3更详细地示出了根据图1的图像处理系统的一部分;
图4示出了第二处理方法;
图5示出了可以应用图4的所提出的方法的示例性测试数据;
图6示出了由根据流程图4的方法当应用于根据图5的测试数据时生成的数值数据的图;
图7示出了通过将流程图图4的方法应用于根据图5的图像数据生成的另外的数值数据;并且
图8和9示出了根据图4所提出的方法可获得的输出数据的示例性可视化。
具体实施方式
参考图1,示出了成像处理装置100的示意性框图。
从广义上讲,图像处理装置(或“超级”系统)IPS包括在图1的左侧示出的图像配准或分割系统以及在右侧示出的配准评估系统。这样,两个系统本身就是独立的系统,但是当然可以结合使用,如图1中所示。更具体地,这里应理解,每个系统可以与另一个系统分开操作,但是如果需要,可以被布置为在超级系统中协同操作,如图1所示。
图像处理系统IPS布置对由成像装置IM采集的示意性示出为A、B的数字图像进行操作。本文设想了任何类型的医学或非医学成像设备,并且特别包括以下非详尽的成像模态列表:诸如CT或C型臂的X射线成像装置,MRI或(SPECT)或PET成像装置,以及超声装置以及其他。尽管可以通过相同的成像模态采集两幅图像A、B,但是本文还设想了多模型场景,其中图像已经通过不同的成像模态采集。例如,图像A、B可以不必都是X射线图像。在针对模态-模态成像的示例中,一幅图像A可以是X射线图像,而另一幅图像B是超声图像。本文中还设想图像模态的任何其他组合。尽管在上文和下文中,参考了“两幅图像”A、B,但这些应该在概念上解释为仅仅是方便的简写,因为当然,系统能够在两幅以上的图像上操作。
针对图像数据的图像A、B包括图像元素。取决于图像数据是分别是2D还是3D,图像元素也称为像素或体素。任何图像元素由2D或3D中的位置定义,其可以通过2D中的坐标(i,j)或3D中的(i,j,k)方便地描述。除了该位置信息之外,还有图像值被分配给这些位置中的任何一个。图像值在数字上表示物理特性,例如X射线图像中的X射线衰减。如果需要可视化,则图像值可转换为灰度值或颜色值。单幅图像值-有时称为“强度”-表示在图像采集时考虑的特定属性的质量或数量。例如,在X射线成像中,图像值或强度(通常映射到从白色到黑色的范围)指示不同程度的X射线衰减。衰减是由X辐射的射线在穿过图像样本(例如患者)之后经过所经历的衰减,并且由探测器元件沿着所述射线的方向“看到”或探测到。图像值本身通过A/D转换经由合适的电路而从X射线探测器的所述探测器元件处的信号导出。尽管在其他成像模态中,这种“信号链”和所使用的检测器技术可能不同,但一般原则仍然是每个图像值表示要被成像的对象与探询信号的相互作用的程度或水平。探询信号例如是X射线辐射,但也可以是超声或核辐射,这取决于所使用的成像模态。
成像值的总和,特别是它们在检测到的图像数据上的变化,产生对图像结构进行编码的对比度。这些图像结构表示成像对象的各个方面。这些结构包括相对均匀的成像值的平台,更有趣的是,这些平台之间的边缘。边缘携带重要的图像信息,因为这些信息例如向查看者提供解剖学细节,例如器官的表面或边界,例如肺或肝脏或动物或无生命对象中的任何其他感兴趣区域。由采集成像装置IM使用的探询信号可以是穿透性的,例如X射线,但是不一定在所有实施例中都如此,因为本文还设想了利用可见光或不可见光的光学成像。
在一个实施例中,所提出的成像系统有助于图像处理中经常面临的重要任务,即配准的图像处理。在(图像)配准中,人们希望建立两幅(或更多幅)输入图像之间的对应关系,比如图像A、B。再次,为了简单起见,我们将以下内容限制为考虑两幅图像,并且理解它可以应用于多于两幅图像。更具体地,在配准中,人们希望在两幅图像中找到对应的图像结构。例如,可以在两幅图像中不同地表示同一器官的边界。可能会是这种情况,因为器官本身在形状或位置上改变,而两幅图像是在不同的采集时间获得的,或者因为成像模态不同。因此,图像配准本质上是一种工具,利用其,用户可以更容易地“重新找到”两幅图像中的相同结构。
现在更详细地讨论成像处理系统的配准功能,可以将配准问题定义为找到从“源图像”的变换T,即图像A到“目标图像”,比如图像B。标示“目标”或“源”图像完全是任意的,并且可以由用户指定。一旦找到变换,就允许将A变换为B,即将T应用于A产生B,简称T(A)=B。然而,确切的对应关系通常不会在现实中实现,并且希望至少找到变换T,其以最佳方式近似目标图像,可通过适当的量度来定义。作为说明,如果使用相同的模态,则对应关系可以关于差异被定义为|T(A)-B|。然而,应当理解,当使用不同的模态时,需要更一般的方法(例如互信息等),因为可能不可能简单地采集差异。
换句话说,图像配准可以被表达为优化问题,其中人们希望通过优化有时被称为代价函数的目标函数来计算配准变换T作为来自给定图像数据A、B(其中,x表示特定图像位置)的配准结果。然后,最终结果T是使目标函数最小化或最大化的变换。在下文中,优化问题被公式化为最小化问题,其中人们希望使代价最小化,但是这被理解为不对下文进行限制,因为本文还设想了最大化问题,其中任务是最大化目标函数,其在该上下文中不再称为代价函数,而是作为“效用”函数。
优化问题通常可以转换为迭代算法。这可以通过应用诸如共轭梯度或牛顿-拉夫森或其他的数值技术来实现迭代过程,给定针对配准变换的初始估计T0,迭代通过一个或多个中间变换TN-j,直到人们到达最终结果TN。迭代次数j(可以容易地达到数百或数千),将通常取决于希望实现的准确度。在其他情况下,人们可以在已经执行预定数量的迭代之后简单地中止迭代。
更具体地,图像配准可以被制定为迭代地最小化所谓的“内力”和“外力”的联合加和。从概念上讲,外力来自所选择的相似性量度(例如,平方差的和,(局部)互相关,互信息,归一化梯度场),并且因此基于基于图像的特征(强度,边缘等)。
另一方面,内力被链接到正则化项(或“变形模型”),所述正则化项在更简单的设置中惩罚位移的空间梯度或者并入更复杂的项,例如来自物理驱动的弹性模型。从初始变形(例如,零位移)开始,在每个迭代步骤中,内力和外力的联合加和减小。如果内力和外力处于平衡状态,即达到平衡点时,则达到收敛。由于代价函数由外部和内部两个分量组成,因此优化任务是在两个(或更多)准则之间进行协商的任务之一:i)解(Tj)相对于根据图像A、B的“输入数据”的相似度(外力),以及ii)解(Tj)在多大程度上符合所需的性质,例如由代价函数的(一个或多个)内力分量定义的平滑度等。
然而,在上述公式的配准任务中,图像伪像或病理变化可能扰乱平衡点,导致在迭代期间提出非现实的大的局部变形。例如,具有高图像强度和锐利边缘的条纹伪影(由金属或造影剂引起)可导致大的外力。类似地,缺失的图像结构(例如,由于信号消失或病理变化)可导致显示局部压缩或甚至折叠的变形。
更详细地并且参考图像处理系统IPS中的配准功能,建议自动检测和抑制在迭代期间引起的不希望的局部变形,从而实现针对共同扰乱迭代配准算法的收敛特性的因素的更好的鲁棒性。具体地,图像处理系统包括输入端口IN,其中,接收要配准的两幅或更多幅输入图像A、B。图像可以直接从成像模态IM或从多个不同的图像模态接收,但也可以从针对图像的数据库或其他存储设施中检索,例如,如在医院信息系统(HIS)中使用的PACS(图片存档及通信系统)。然后将输入图像A、B传递给迭代求解器SOLV-REG,迭代求解器基于代价函数实现迭代配准变换。求解器从初始变换T0迭代通过一个或多个中间配准TN-j,一旦达到足够的收敛水平就到最终变换TN。迭代求解器由优化函数f驱动。优化函数f是至少两个函数分量的加和。一个分量是正则化项,而另一个分量D是相似性量度。
更正式地,在本文中在一个实施例中我们假设一种配准方法,其中(至少)两个函数的总和将被最小化,其中一个函数D对应于相似性量度而另一个函数S对应于正则化项,
D[u]+S[u]→min (1)
其中,u的范围在配准变换T的空间上。u是变形矢量场,可以关于其来描述变换T。这是因为通常不能获得针对T的以“公式”形式的封闭解析表示。项T和u在本文中可互换地使用,因为一个定义另一个。然后,优化这两个函数D、S的联合和作为代价函数可以然后被示出为等效于迭代求解偏微分方程组(x表示图像位置):
Lu(i+1)(x)=u(i)(x)+f(u(i)(x)) (2)
利用适当的边界条件,并且i是迭代步骤。更新或优化函数(2)可以通过应用微积分方法从优化问题(1)导出,并且PDF问题(2)可以通过有限差分方法、通过有限元方法、通过有限体积方法或通过这些方法的混合或通过其他方法来求解,其中,(2)变为矩阵方程组。
更新函数(2)告诉我们如何迭代通过可能的变换的空间。L定义了源自正则化器S的微分算子,u(i)定义了在迭代步骤i中计算的变形矢量场(DVF),并且f定义了源自相似性量度D的函数。正则化项(由L表示)被链接到内力,相似性量度(由力项f表示)被链接到外力。可以看出,本文假设代价函数(1)的结构使得它在功能上可分离成针对外力和内力的不同函数分量D,S,并且本文中还假设从(1)导出的更新函数(2)“继承”该可分离性,以类似地产生分别与D,S相关联的两个有关的函数分量f,L。作为解释性说明,在此将理解,代价函数以及更新函数当然可以包括比针对内力和外力的函数分量更多的函数分量,以便考虑其他优化准则。
相似性量度对目标图像与迭代期间产生的中间结果有多不同进行量度。更具体地,在迭代期间,将中间结果TN-j重复地应用于源图像A,并且然后通过合适的量度来评估该图像TN-j(A)和目标图像B的相似度。合适的相似性量度是A与TN-j(A)之间的最小平方差或其他,特别是不基于差异运算的相似性量度。
在每次或(任何第m次(m>2))次迭代之后(或随机地),(图像)体积元素监测器IVM操作以监测当在迭代周期中应用中间变换TN-j时给定位置x处的图像元素的体积如何变化。优选地,体积元素监测器IVM操作以监测每一个图像位置x,但是设想监测限于图像位置的真子集实施例。这些图像位置可以是用户可定义的。
(图像)体积元素评估器IVE从图像体积元素监测器IVM接收监测结果,并通过将每个图像位置x的监测的变化与预设条件进行比较来评估或评价体积元素的体积变化。特别地,在图像配准过程期间,存在对x周围的各个邻域的体积变化的逐体素或逐像素方式自动监测。每个图像元素由这样的邻域定义。例如,如果对于特定体素,体积变化落在信任区间或其他接受条件之外,则针对该体素位置x标记该体积变化。例如,可以通过要求对于所有体积变化存在可允许的压缩因子<2且扩展小于100%来定义信任区间。
当图像评估器IVE发现,对于给定位置x,如此监测的改变违反所述条件,则将标志信号传递到校正器模块C。在接收到该信号时,校正器模块操作以针对给定图像元素位置x相对减轻(即,完全或至少部分地禁用或抑制),更新函数中的在功能上与相似性量度D相关的函数分量。更新函数中的该相似性量度分量被指定为fin(2)。特别地,对于剩余的迭代,相似性量度分量对于该图像位置x保持无效。换句话说,修改该像素的迭代公式,以便从代价函数中排除与相似性量度有关的函数分量。换句话说,对于该像素的优化问题现在已经被局部地改变为,对于该像素x,不考虑相似性量度。
通过全部或部分禁用或抑制的自动(相对)的减轻可以通过将(2)中的项f(表示所述外力)逐体素地与以掩模矩阵或掩模图像M以及M中与图像B中的图像位置x相对应的行和列相乘来实现。掩模图像M在全部位置以“1s”初始化,并且针对任何图像位置x,1s被设置为0,在自动监测步骤期间设置针对其的标志。换句话说,伪影引起的“大”外力得以相对减轻(或者在一些实施例中完全关闭),并且在剩余的图像配准过程中,找到新的平衡点。新的平衡点表示不存在伪影或病理变化的场景。更正式地说,在此提出将(2)修改为:
Lu(i+1)(x)=u(i)(x)+M(x)f(u(i)(x)) (3)
在M是掩模矩阵M的情况下,通过抑制或禁用更新函数分量f的函数动作来选择性地相对减轻。优选地,函数分量f的动作被完全抑制,但是在一些应用场景中,仅通过乘以小的因子(<<1)来减轻就足够了。应当理解,(2)、(3)中的更新函数的公式仅仅是示例性的,因为公式对应于“首先优化。然后离散”方法。在该方法中,人们首先通过在更新函数方面通过任何方法获得算法指令来优化代价函数,并且然后将其离散化以最终达到针对更新函数的计算机可实现的形式。但是在本文中还在替代实施例中设想了反过来的方法“首先离散化然后优化”,并且对于相同的代价函数,该方法可以得到针对(2)和(3)的完全不同的形式。
在迭代结束时,在输出端口OUT输出最终配准结果,即最终变换TN(=uN)。这里考虑的变换是非刚性、非仿射类型的,以允许考虑图像元素位置x处的剪切、压缩或扩张变化。变换T本身通常关于变形矢量场(DVF)来描述,变形矢量场继而可以以合适的函数表达来描述,例如矢量集合或者作为一个或多个矩阵或任何其他数据结构。针对最终配准变换结果TN的该函数表达式可以应用于源图像A以便将A“配准”到“B”上。换句话说,TN(A)是对B的估计。取决于配准的质量,预计TN(A)在视觉上与目标图像B大致对应。
现在参考图2左侧的流程图,其中更详细地解释了迭代配准算法的操作。这里将理解,对配准方法的不同步骤的以下解释不一定绑定到根据图2的体系结构,因此,本领域技术人员可以理解以下步骤本身构成教导。
在步骤S210,接收用于配准的两幅或更多幅图像A、B。
步骤S220包括通过朝向最终结果TN迭代通过一个或多个中间结果来迭代地求解最终配准变换。迭代由(1)驱动并且实现为(3)。迭代在迭代公式寻求改进代价函数的意义上被“驱动”。换句话说,下一个配准的结果Tl不比上一个配准结果Tk(l>k)引起更多的代价(并且实际上通常将降低代价)。根据一个实施例,针对配准结果的迭代求解可以如上根据(3)那样地实现如下。
在步骤S230中,监测图像体积元素在给定图像位置x处的变化。这通常包括针对每个被监测的图像位置,具有合适尺寸的邻域(球体或立方体等)的定义,所述尺寸可以是用户可调节的。当在迭代周期期间将中间变换配准结果应用于那些图像位置x时,在所述位置x处将发生相应体积元素的相应图像体积变化,并且可以通过向量场按x来描述这些局部体积变化。这些矢量场可以描述相应(按x)邻域的局部体积变化,例如压缩或扩张,或剪切,并且这些变化可以通过适当的量度来量化,例如比例因子或叉积值,以对剪切等进行量化。量换句话说,针对它们的定义使用的量度和几何工具最终将取决于监测的变化的性质。
在步骤S240,通过在数值上比较量化的变化与预定义的(并且可能是用户可调节的)阈值来建立步骤S230处的按x的评估。更具体地,在小的或中等的变形的情况下,在收敛之后,我们具有M≡1,因为不需要抑制。然而,在导致临时大的局部体积变化的伪影或病理变化的情况下,掩模M被局部地设置为0,这意味着要优化的整体函数被修改。对于当前迭代步骤i和对于所有后续的迭代步骤j>i,这继而将当前位置x处的力项f归零。在数学上说,在下一个迭代步骤中,以当前变形矢量场u(i)作为起始点来开始新的优化过程。由于M(x)=0,对于具有不可靠的高体积变化的所有先前监测的位置x,因而取消了针对那些图像元素位置x的外力。这意味着内力是所述位置x处在迭代优化中唯一的“活动参与者”。在优化过程中仅考虑内力(和可能的其他函数分量),同时不再考虑外力。基于内力(但不是外力)在给定位置x处进行优化将导致围绕该图像位置x的平滑内插变形。
在步骤S250,相对减轻了对于在所考虑的图像位置处的后续的迭代的相似性量度的影响。在步骤S240的评估已经确定已经违反了预设接受条件的条件下,相对减轻了该效果。通过禁用迭代算法中的相应函数分量或从代价函数导出的更新函数,可以相对减轻相似性量度的影响。通过减轻相似性分量可以实现相对减轻,同时正则化项不受影响。替代地并且与此相反,可以放大正则化器,同时相似性量度分量保持不受影响。作为进一步的变型,还设想减轻相似性量度分量并放大正则化器分量以使相对减轻更加明显。
首先图示通过禁用或抑制或减小加权操作来减轻相似性量度分量的实施例,这可以通过定义最初填充有单位元素的掩码矩阵来实现。每个条目或“标志”对应于要配准的图像中的图像位置,因此掩模矩阵M是2D或3D,取决于图像的维度。如果评估显示违反了接受条件,则掩模矩阵中的相应位置变为零。然后将迭代的更新函数与该矩阵相乘,使得针对给定位置的相似性项由于在掩模矩阵M中的相应位置处的零项而被取消。将理解,要减轻的外力分量不一定被集成到单个函数单元中(例如在(2)中示意性地表示为f),而是外力分量可以在更新函数中在多于一个函数子分量中展开(f=f1,f2,f3,..)。然后在这种情况下,掩模矩阵包括多个矩阵M=M1,M2,M3,...,每个矩阵作用于相应的分量f1,f2,f3上,以对至少一些执行减轻动作,优选地,在所有子分量上。如上所述,代价函数以及因此更新函数可以包括除外力和内力之外的分量以便考虑其他优化准则,并且这些其他准则不受所提出的减轻动作的影响。也就是说,本文中设想减轻动作以仅作用于源于代价函数中的相似性量度的(更新函数的)外力分量。
上面的公式(3)可以由偏微分方程组(偏微分方程)(PDE)表示。在下文中,呈现伪代码,其形成用于上述图像处理方法的一个实施例。更具体地,根据一个实施例,所述方法可以实现为:
针对所有体素位置M(x)=1,设置掩模M(x)=1
对于所有的i=1,2,...直到收敛,循环进行以下操作:
通过求解决偏微分方程组来计算u(i)
计算局部体积变化V(x):=det(Jac(x+u(i)(x)))
针对所有体素位置x
如果V(x)在信任区间之外,则设置M(x)=0
结束
其中,在第一个迭代步骤中,u(o)≡0(或任何给定的,先验已知的,初始位移)。在上文中,算子Jac(.)表示雅可比行列式(矩阵),而det(.)是表示取雅可比矩阵的行列式的算子。
这里应理解,上述伪代码仅仅是实现所提出的方法的许多实施例中的一个,并且不应被解释为限制本文中所描述的内容。
在替代实施例中,可以以与上述实施例相反的方式来实现相似性量度分量的影响的相对减轻。换句话说,相似性量度部件保持不受影响,而是相对于相似性量度部件放大或加权的正则化器。在这种情况下,掩模M在(3)中重新排列以对项L进行操作而不对项f进行操作。此外,伪代码将需要被修改,因为M中的条目现在默认处处为“1s”,但现在在违反接受条件的图像位置处设置为因子>1。
如上所述,人们也可以分别通过在f和L处使用相应的掩模矩阵来组合这些实施例,减轻f和放大L以增加f的相对减轻的有效性。
所提出的图像处理方法和用于配准的图像处理系统的一种可能的应用场景可以描述如下:所提出的方法的自校正特征(即,在迭代期间响应于个体图像位置处的变化自动抑制更新函数的部分)可以在一个实施例中适用于剪切敏感的呼吸运动估计。这是因为已知肺在裂隙和胸膜的呼吸期间经历滑动运动。最新的配准算法返回变形矢量场,其中跨胸膜或裂隙的任何潜在剪切被平滑。自动剪切检测可以用于局部地减小内力的权重,从而允许局部地滑动,同时仍然在全局上满足正则化项所规定的物理模型。可以通过形成局部变形矢量场的叉积来实现自动剪切监测。
如前面简要指出的,图像处理系统IPS可以被布置为执行分割而不是先前描述的配准。虽然分割和配准是替代操作,但是本文可以设想这样的实施例,其中功能,配准和分割都集成在同一图像处理系统中,并且用户可以选择要执行哪个操作-分割或配准。
简要地返回参考图1,在分割的场景中,图像处理系统IPS包括代替配准求解器SOLV-REG的迭代分段求解器SOLV-SEG。
与配准应用的情况一样,分割任务同样可以关于具有两个函数分量的代价函数的优化问题来制定。因此,基于监测每个图像位置的个体图像值或图像元素,在迭代期间禁用某些函数分量的关于配准的上述原理也可以应用于图像分割。
更具体地,在基于模型的分割中,希望使用几何模型,例如网格模型,其表示希望分割的对象的粗略形状。然后将该模型几何地改变为不同的状态,例如压缩,平移等,以适配到要被分割的图像中的图像结构。然后,任务是找到最佳图像结构,即最适合模型同时包括正确的图像特征的图像结构,例如在正确的位置具有边缘。这可以再次表述为联合优化函数。也就是说,在基于模型的分割中,通过联合优化外力(基于图像特征)和内力(网格正则化),基于网格的形状表示被迭代地适应于例如器官表面。网格的大的局部压缩、扩张、弯曲或折叠可以被检测,并且被用于局部抑制外力。
然后,针对SOLV-SEG的迭代分割在初始几何状态处针对分割模型进行迭代,通过所述模型的中间状态朝向模型的最终状态。在该场景中模型的最终状态意味着图像中的位置和变形,其正确地概述了要被分割的图像结构。例如,在心脏图像中,人们可能希望分割表示植入心脏瓣膜的图像结构。因为心脏瓣膜的类型是已知的,因而其预期的形状和尺寸等也是已知的。然后这构成了所述模型,其正在被拟合(拉伸,压缩或者“可容许的”扭曲)到心脏图像,从而找到表示瓣膜的图像结构。
在迭代结束之后,然后可以将变形的(并且可能移位的模型)叠加在输入图像上以在视觉上勾勒出瓣膜形状。
对应于之前在图1和图2(的左侧部分)中讨论的图像配准实施例,在图像处理系统IPS的分割实施例中,现在包括监测几何变化的模型监测器MM(而不是IVM),所述几何变化是所述初始模型在迭代期间所经受的。代替图像评估器IVE,然后存在模型几何变化评估器MCE。这被配置为评估几何变化并将这些变化与接受条件进行比较。接受条件定义了在当前分割环境中接受的允许的几何变化的类型。如果在监测期间出现在迭代展开时提出不允许的模型改变,则向校正器通知此情况,并且对于相应的图像位置,相对且选择性地减轻更新函数中的图像特征分量,如上在配准场景中所述。更具体地,禁用更新功能中的分量,所述分量对应于代价函数中的已从其导出更新函数的图像特征分量。换句话说,在迭代期间相对于图像特征分量优先考虑模型分量。
现在返回参考图2,但现在参考其右侧,这示出了用于基于模型的分割的图像处理方法的流程图。
在步骤S210',接收要分割的输入图像A。
在步骤S220',执行迭代分割算法以求解最终配准分割结果。在此过程中,迭代通过初始模型的一个或多个中间几何状态。如前所述,在配准的上下文中,迭代求解器由用于分割任务的优化函数驱动。优化函数包括至少一个(优选地两个)函数分量,其中一个函数分量在迭代期间测量图像中的图像特征。这些特征是在相对于模型的当前几何状态(位置、压缩、扩张等)的不同位置处测量的。以此方式,试图使模型与周围的图像信息协调。更具体地说,一方面,模型的边界应该与图像中的边缘结构相符,并且这通过使用目标函数中的图像特征分量来建立。另一方面,要应用于模型的几何变化不得违反特定正则化条件。例如,禁止不切实际的几何变化,例如将模型分解为两个。换句话说,几何要求也需要被观察,并且这通过形成第二函数分量的正则化项表示。如果存在冲突,则本文提出的方法优选观察正则化要求并放弃图像特征贡献。
更具体地,在步骤S230',监测模型经历的几何变化。几何状态变化包括体积变化和/或拓扑考虑,例如连通性和/或其他量度。
在步骤S240',将那些变化与针对那些几何状态变化的接受条件进行比较。例如,这些可包括最大变形量(剪切等),允许的最大体积变化量(如上面讨论的配准实施例中)和/或是否保持连通性,或者是能够定义几何状态的任何其他合适的、可量化的几何条件。
在步骤S250',如果在先前的步骤S240'发现那些条件已被违反,则相对减轻图像特征分量或其在迭代更新函数中的影响。在足够次数的迭代之后,获得模型的当前几何状态的形式的最终结果。这包括其相对于图像坐标系的形状和位置。迭代结束后成像模型的形状和位置定义了分割。换句话说,在其当前几何状态(在迭代结束时)由模型限定的图像部分内的所有图像元素然后被定义为被认为构成感兴趣的图像结构。相对减轻的概念类似于先前在配准实施例中所解释的。
在根据图2的针对配准或分割的实施例中,设想了替代实施例,其中,代价函数仅包括一个(单个)函数分量,即,相似项和图像特征项。特别是在这种情况下,减轻仅被应用于该单个项。特别地,尽管在本文中是优选的,但是正则化项的存在并不针对所有实施方案都是强制性的。结合图1简要提及的配准精度评估系统REG-EVAL在图3中更详细地给出。
本文中提出的配准评估子系统REG-EVAL允许自动估计配准精度,并且根据图像元素评估配准精度。换句话说,配准精度在各个像素或体素位置处量化,并且因此可以针对图像的所有图像元素示出,但是在本文中在一些实施例中也设想由所提出的评估器仅处理图像的一部分或其子集。评估器的输出是配准精度图,即具有与配准到图像T(A)上相同尺寸的图像,所述图中的每幅图像像素或体素指示适当尺寸(例如毫米)的精度、厘米或微米。
所提出的图像评估器的构造(将在下面更详细地解释)使得其响应在人类观察者可能选择以来估计配准的准确性的图像区域敏感。更具体地,边缘状结构(例如血管、器官或肿瘤描绘)是这样的图像位置,即其中通过构造,所提出的评估器将返回非零响应,而在具有均匀图像强度的区域中(在其处可以预计甚至人类观察者在判断准确性上也有困难)将返回零响应。更具体地,如果估计了边缘之间的不匹配或非最佳匹配,则返回非零响应。当(i)针对该图像元素的配准是精确的或当(2)这是均匀强度的区域时,预期零值作为响应。如本文所使用的,本文使用的“零响应”的概念表示在给定位置未发现错配或配准错误的信息的响应。
然而,更具体地说,设B和T(A)是比较的两幅图像-其中B是目标图像,例如随诊图像,而图像T(A)是在配准到源图像上后经变形的基线图像A。这些图像在输入端口IN处被接收。然后,任务是将每图像元素量化到T(A)与B类似的程度。通常,由于不准确性,并且确切地是所提出的评估器对所响应的这些不准确性,B和T(A)中的各图像边缘可能在两幅图像B、T(A)中不在相同位置处延伸。这是通过以下方式来实现的:所提出的评估器使用通过一些参数(即σ2)参数化的一系列边缘滤波器EF。诸如边缘检测滤波器、Sobel滤波器或二阶图像导数运算符等的图像滤波器通常将在边缘图像结构处返回响应。该响应可以图形化地显示为通常在感兴趣的图像结构处具有峰值特征的响应面。族中的个体边缘滤波器EF在它们的响应特性方面不同,其由参数σ2量度。例如,在一个实施例中,各个峰值的响应峰值宽度不同。更具体地,并且在一个实施例中,边缘滤波器是具有高斯内核的二阶图像导数和对应于内核中使用的高斯分布的方差σ2。然后,这相当于在所述参数σ2变化时返回的峰的不同的峰宽度。测试模块TE操作以将系列EFσ2中的每个边缘滤波器EF分别应用于两幅输入图像B和T(A)。因此,对于每个边缘滤波器EF,它产生多个响应面(或(一个或多个)响应阵列)对,其中存在针对输入图像B和T(A)的相应响应面。在一些实施例中,代替二阶导数或者除了二阶导数之外,可以使用高于二阶的图像导数。
非常宽泛地,然后通过组合器Σ将多个响应面对组合成组合响应面。这种组合响应面的配置使其随着滤波器参数σ2非单调地变化。然后,极值查找器SE操作以在组合响应面中找到作为滤波器参数的函数的极值。应当理解,作为实际情况,仅使用有限数量的滤波器元件(并且因此滤波器参数),并且因此仅在离散采样点的数量上精确地知道组合响应面。因此,进行插值以构建针对这些控制点的连续表面,然后检查该连续表面的极值点。换句话说,如此定义的滤波器参数值因此不一定对应于用于获得滤波器响应面的任何初始多个边缘滤波器。
输出端口OUT输出在其处假定组合响应曲面的极值的(可能是插值的)滤波器参数。针对每个图像像素或更一般地图像元素重复该操作,以产生不同的滤波器参数。这样返回的每个滤波器参数然后是相应图像位置处的配准精度的量度。
在输出端口OUT处产生的个个参数值σ2,其可以被称为准确度度分数,然后可以在显示设备MT上通过查看器VIS的操作来可视化。
更具体地,现在将以广义的术语描述所提出的评估器的实施例,并且稍后将参考图4的流程图更详细地描述该实施例。宽泛地说,通过计算T(A)和B中的二阶导数来分析图像边缘的未对准。图像边缘由二阶导数图像中的正峰值和负峰值形式的响应面表示。更具体地,图像边缘由梯度矢量表示,其中,其分量中的至少一个表现出具有所述正峰值和负峰值的响应面。正峰值与负峰值之间的空间距离取决于导数核–例如,选择的高斯分布的标准偏差越大,两个峰值之间的空间距离就越大。假设T(A)和B中的图像边缘未对齐,对于一定大小的导数核,来自图像T(A)的二阶导数的负峰值与来自图像B的二阶导数的正峰值匹配。如果未对齐较小(较大),则导数核需要较小(较大)以匹配两个峰值。可以通过计算A和B的二阶导数的标量积,然后进行阈值处理和积分来测量峰匹配的良好性。如果峰值处于相同的空间位置,则得到的值最小。计算针对各种导数核大小的匹配良好性允许确定最小化标量乘积的内核大小,并且因此允许自动确定图像T(A)与图像B中的未对齐的图像边缘之间的空间距离。
计算的空间距离(或配准误差或精度估计)可以被颜色编码并覆盖到图像B上和/或(在变换之后)覆盖到图像A上。
应当理解,在刚刚描述的示例中,组合响应面通过应用标量积并且然后在图像位置x上对标量乘积的值进行积分并且相对于零或任何其他参考基线值进行阈值处理而形成。仅将负的标量乘积值相加。然后,在这种情况下,组合响应面的物理意义可以被认为是各个峰之间的空间位置,其被表示为峰之间匹配的良好性。
现在参考图4的流程图,其中解释了用于配准精量度化的成像处理方法。同样,如上面在图2的流程图中,所解释的方法步骤和实施例不一定与图3中所示的体系结构相关联,并且以下方法步骤本身构成教导。
在步骤S410,接收两幅(或更多幅)输入图像T(A)和B。图像直接从成像装置接收,或者根据用户或协议请求从存储设施(例如PACS或任何其他存储器或数据库系统)检索。
在步骤S420,分别将一系列参数化边缘检测滤波器应用于所述两幅输入图像。以此方式,产生每幅输入图像和滤波器参数的参数化的一系列不同滤波器响应的。边缘检测滤波器通过滤波器参数来参数化。该参数涉及相应边缘滤波器的特性,例如峰值响应的宽度等。
在步骤S430,根据两幅输入图像的响应面对中的信息在组合操作中被合并到组合响应面中。该组合响应或组合面被构造成作为所述滤波器参数的函数非单调地变化。
在步骤S440,建立作为滤波器参数的函数的组合响应面的极值。然后返回与组合响应面中的该极值相关联的滤波器参数。可以通过从来自边缘滤波器的集合的现有参数进行插值来获得该参数。
然后,用于量化配准精度的量度基于这样返回的滤波器参数,所述参数对于不同的图像位置可以是不同的。换句话说,将对图像中的每个图像元素或其子集重复先前的步骤。滤波器参数可以直接在屏幕上显示,或者可以转换为长度尺寸,例如毫米或微米。
在下文中,现在将以测试图像数据的示例描述用于图像配准量化的上述方法的实施例。在图5的窗格A)和B)中示出了这样的测试图像数据的示例。为了便于表示,测试数据被认为是二维的,其中,图5中的水平x和垂直y轴分别表示列i和行j。在以下意义上测试数据是人为的,即只有一条边缘在图像中垂直延伸,但是在不同的位置,因此表示配准不准确。
在执行所提出的方法的过程中获得的数值数据显示在图6的图的集合中。每行(在图6中了示出四行)对应于针对高斯情况的特定的滤波器参数σ2,具体地,σ2=4、6、8和10,分别从顶部到底部。为了进一步简化方法步骤的图形说明,在图中仅示出了沿x轴的一维呈现。因此,出于图6、7中的说明目的,响应“面”和组合响应“面”是1D中的曲线。
图6第一列中的图表是相同的,并且仅示出两个测试输入图像的概况,其中,边缘出现在不同的图像位置。第二和第三列(从左到右)示出了二阶导数的相应构造。为了便于图形显示,第二列仅示出一阶导数的X分量,并且第三列仅示出二阶导数算子的纯偏分的X分量。最后一列示出了标量乘积的逐点取得,并且该操作是形成图7所示的组合响应函数(也称为“指标函数”)的一部分。具体地,图7中的组合响应面是通过在相对于零进行阈值处理之后进行积分(沿x轴)从标量积获得的。标量积可以被理解为仅用于逐点操作的一个实施例,以通过所考虑的图像位置沿着(每个)空间维度将空间信息折叠成单个权重。然后,优选地沿着每个空间维度,对权重进行加和。在一个实施例中,如在标量乘积实施例中那样,权重是有正负号的利用阈值处理来实现求和。然后对每个图像位置重复这些步骤。
被视为在滤波器参数σ2的函数的极值被看作是在10的量级,并且在该示例正是是该值被返回作为在给定的图像位置处的配准确度量化。换句话说,对于约为10的针对σ2的值,针对图5中的图像A)负峰值(以实线显示在图6中,第3列,第3行)和在图5中的图像B)中的正峰值(以虚线示出)出现在大致相同的体素位置处。
可以看出,对于均匀区域,将返回“零”响应,因为如果两个峰没有“交叠”,则标量乘积为零(针对两幅输入图像分别显示为实线和虚线)。因此,通过构造,所提出的方法将仅在边缘足够接近地靠近(即,通过改变响应面中的峰的宽度可以实现至少部分交叠)的图像位置处产生非平凡响应。因此,如果我们检查图6中的第三列,可以看出,仅由于峰之间的相对紧密的空间距离,人们可以通过改变滤波器参数来在峰之间建立交叠。
继续参考图6,基于二阶边缘检测来形成组合响应曲线的步骤可以在一个实施例中实现如下:
为了更好地理解,我们首先假设具有单个边缘结构的图像的情况。然后,对于两个给定的输入图像A和B(为了符号方便我们删除较早的符号T(A)并且现在仅将该图像称为“图像A”),对于二阶导数图像牌子EF的情况,我们建议使用以下子例程来计算边缘响应曲线,通过σ2进行参数化。
首先,针对两幅输入图像A、B计算关于x、y和z的一阶常规导数。这得到中间图像Ax,Ay,Az,Bx,By和Bz,即,在从体素位置(i,j,k)到实数集合的映射中。优选地,对一阶常规导数执行归一化,得到图像Ax n,Ay n,Az n,Bx n,By n和Bz n。图6的第2列,是Ax n和Bx n的示例(在图6中仅示出了x分量)。
接下来,根据Ax n和Bx n计算两幅图像相对于xx,xy和xz的二阶常规导数,得到Axx,Axy,Axz,Bxx,Bxy和Bxz。再次,我们将这些归一化并存储Axx n和Bxx n。图6的第3列是仅x分量的可视化。
接下来,针对两幅图像来根据Ay n和By n计算关于yx,yy和yz的二阶常规导数,得到Ayx,Ayy,Ayz,Byx,Byy和Byz。我们再次归一化和存储Ayy n和Byy n
接下来根据Az n和Bz n我们计算针对两幅图像的关于zx,zy和zz的二阶常规导数,得到Azx,Azy,Azz,Bzx,Bzy和Bzz。归一化产然后得到Azz n和Bzz n并且它们然后被存储。
这些操作然后结束计算针对给定的σ2的响应曲线,并且针对一范围的参数σ2重复以上计算。
然后,组合响应曲线的计算包括针对每个体素位置(i,j,k)计算S(i,j,k):=Axx n*Bxx n+Ayy n*Byy n+Azz n*Bzz n。针对此的示例性值示于图6的第4列中示出(再次仅针对x分量)。
组合响应曲线的计算还包括:对S的所有小于0的值进行求和,并且该结果被存储为与σ2相关联矢量指标函数的求和的值(σ2),以便形成组合响应函数。人们然后确定指标函数最小的σ2(参见图7中的示例)。然后存储最优σ,任选地与其函数值一起存储在当前体素位置的输出图像中。
归一化步骤允许使滤波器响应即使对于来自不同模态的图像也是可比较的。在图像导数计算的每个阶段可以应用或不应用归一化,但是如上所解释,在每个阶段应用归一化是优选的。根据一个实施例,归一化可以如E.Haber等人在“Bildverarbeitung für dieMedizin 2005”,Springer(2005),第350-54页中的“Beyond Mutual Information:ASimple and Robust Alternative”中所进行的。
如上所述,上述子程序适用于具有单个(突出)边缘结构的图像,如图5,A),B)中的测试图像中。对于具有更复杂结构的真实世界图像,上述子程序应用于输入图像的子集,然后优选地针对所有体素位置重复。
因此,本发明提出,针对两幅给定图像A和B,针对每个体素位置定义以相应体素位置为中心的足够小的邻域。在一个实施例中,选择21×21×21的尺寸,但这仅仅是示例性的。然后,执行上述子例程,得到(i)表示为“σ2-opt”的最佳σ2,以及(ii)表示为“指标函数(σopt)”的针对σ2-opt的组合响应功能的函数的值。
值σ-opt然后可以在当前体素位置被存储在一个第一输出图像Σ中,并且任选地,人们可以将指标函数(σ2-opt)存储在第二输出图像F的当前体素位置处。在不同的位置i,j,k的σ2-opt的集合Σ可以然后被转换成空间尺寸,如毫米,空间尺寸然后可以以彩色编码来可视化,作为针对估计的估计配准误差的估计的空间分辨图。图8(A-C)和图9是针对颜色或灰度值编码可视化的示例。任选地,人们可以执行Σ的平滑或聚类。
任选地,误差图可以叠加在原始图像数据集A或B上。
在另外的可视化实施例中,在输出图像Σ中,将所有条目设置为零,其中在第二输出图像F中,各个值高于特定阈值。更精确地说,对于每个体素位置(i,j,k),如果F(i,j,k)>阈值,则Σ(i,j,k)=0。该阈值处理允许实现更好的噪声鲁棒性。
在一个实施例中,可视化包括引导功能,以引导观看者通过图像数据集到具有最大估计配准误差的图像区域,并且首先显示这些区域。
在一个实施例中,根据估计的配准错误的大小,用不同的参数设置重新运行配准,以便自动改进配准结果,并且这可以根据用户请求任选地重复。
尽管仅针对一个图像尺寸(即x维度)解释了图4中的方法的操作,但是操作被理解为通常针对每个空间方向x,y,z被应用。在上文中,标量乘积操作被广义地解释为将响应曲线的任何逐点组合包括到权重以及沿一个、多个或所有空间维度的那些权重的加和。
如果在相同位置处发生强度变化,则上述方法中的准确度估计基于两幅图像彼此相似的假设。对于以不同模态采集的图像并因此被编码的不同类型的解剖结构,对于两幅图像,具有改变的强度的图像位置的集合可能是不同的。例如,可能的情况是,诸如表示肝脏病变的特定图像结构仅在由第一模态采集的图像中被编码而不在由第二模态采集的图像中被编码。对于我们的假设,这意味着一幅图像中缺失的强度变化(“边缘”)并不一定意味着图像未对齐。因此,可以说,所提出的配准精量度化器遵循以下规则“如果强度改变–如果有的话–发生在相同的位置,则两幅图像彼此相似”。换句话说,可以避免“假阳性”响应。仅在一幅图像中的位置处存在边缘结构而不存在另一幅图像中的边缘结构将不会产生误对齐响应。这可以在图6中很容易看出,其中,组合响应曲线基于逐点标量积计算:如果边缘不在同一位置或在所述位置周围的特定邻域之外,这将产生零。
现在参考图8和9,其示出了通过将所提出的方法应用于更复杂的图像数据而获得的示例性图像。
图8示出了与图5中的简单图像相比具有线性增加的移位的合成测试数据A),B),其中,边缘仅具有恒定的移位。图8A)和B)是不完美配准输出的模拟,其中边缘仅在左图像边界处完美匹配。朝向正确的图像边界,配准不匹配呈线性地增加。已经根据上述方法估计了配准不匹配。估计的配准误差在图8C)中示出。
图9图示了将图4中的方法应用于实际测试数据。图9中的误差图像基于一对冠状CT肺部图像数据集(在图8,D)呼气状态下和图8,E)吸气状态下采集),并且在迭代配准方案中成功配准,导致变形的吸气图像图8,F),其几乎等于呼气图像。
出于概念验证的目的,重复根据图8D)-F)的配准,但是这次分别排除最后30个,最后40个和最后50个迭代步骤,因此产生三个人为变差的结果。然后针对这三个集合中的每个集合来量化配准精度,其中图9的右上方示出了最后30个迭代步骤被丢弃时的配准误差,图9的左下方示出了最后40个迭代步骤被丢弃时的配准误差,并且右下角示出了最后50个迭代步骤骤被丢弃时的配准误差。未对齐的量是灰度值编码的,并且可以看出,与由于迭代步骤的数量减少的较差的配准精度直接相关。图9中的左上方示出,作为参考图像,实际上没有配准误差,对于允许配准算法运行直到收敛的情况。
图1和图3中的各种部件和模块可以如图所示布置为单独的功能单元或模块,但也可以合并为更少的,特别是单个单元。例如,在一个实施例中,部件被形成为单个软件套件中的不同例程,其可以由通用计算机执行,例如与特定成像器相关联的工作站或者与具有一组成像器的网络相关联的工作站。还设想了除软件之外的实现,例如通过编程一个或多个现场可编程门阵列FPGA来执行各种功能。这些部件也可以实现为硬连线专用IC(集成电路)。
在本发明的另一示范性实施例中,提供了一种计算机程序或计算机程序单元,其特征在于,其适于在合适的系统上执行根据前述实施例中的一个的方法的方法步骤。
计算机程序单元因此可以存储在计算单元上,其也可以是本发明的实施例的部分。该计算单元可以适于执行上述方法的步骤或引起上述方法的步骤的执行。此外,其也可以适于操作以上描述的装置的部件。所述计算单元可以适于自动地操作和/或执行用户的命令。计算机程序可被加载到数据处理器的工作存储器中。数据处理器因此可以被装备为实施本发明的方法。
本发明的该示例性实施例覆盖了从最开始使用本发明的计算机程序和借助于更新将现有程序转变为使用本发明的程序的计算机程序。
另外,计算机程序单元可以能够提供所有必要的步骤来完成如以上所描述的方法的示范性实施例的流程。
根据本发明的另一范例性实施例,提出了一种计算机可读介质,诸如CD-ROM,其中,所述计算机可读介质具有存储在其上的计算机程序单元,所述计算机程序单元由前部分所描述。
计算机程序可以存储和/或分布在与其它硬件一起或作为其它硬件的部分来提供的合适的介质(特别是,但不一定是非瞬态介质)中,例如光存储介质或固态介质,但也可以用其它形式来发布,例如经由互联网或者其它有线或无线电信系统。
但是,也可通过类似万维网的网络提供计算机程序,并且能够从这样的网络将计算机程序下载到数据处理器的工作存储器中。根据本发明的另外的示范性实施例,提供了一种用于使得计算机程序单元可供下载的介质,所述计算机程序单元被布置为执行本发明的先前描述的实施例中的一个。
必须指出,本发明的实施例参考不同主题进行描述。具体而言,一些实施例是参考方法型权利要求描述的,而其他实施例是参照设备型权利要求描述的。然而,本领域技术人员以上和以下描述可以得出,除非另行指出,除了属于同一类型的主题的特任的任何组合之外,涉及不同主题的特征之间的任何组合也被认为由本申请公开。然而,所有特征能够被组合,提供超过所述特征的简单加和的协同效应。
尽管已经在附图和前面的描述中详细图示和描述了本发明,但是这样的图示和描述应当被认为是说明性或示范性的,而非限制性的。本发明不限于所公开的实施例。本领域技术人员通过研究附图、公开内容以及从属权利要求,在实践请求保护的本发明时能够理解并且实现对所公开的实施例的其他变型。
在权利要求书中,词语“包括”不排除其他元件或步骤,并且词语“一”或“一个”不排除多个。单个处理器或其他单元可以实现在权利要求中记载的若干项目的功能。尽管在互相不同的从属权利要求中列举了特定措施,但是这并不指示不能有利地使用这些措施的组合。权利要求书中的任何附图标记不应被解释为对范围的限制。

Claims (4)

1.一种图像处理系统,包括:
输入端口(IN),其用于接收两幅或更多幅输入图像以进行配准,所述两幅或更多幅输入图像包括源图像和目标图像;
迭代求解器(SOLV-REG),其被配置为通过迭代通过至少一个中间配准变换T(N-j)来求解最终配准变换T(N),其中,所述迭代求解器由包括至少一个第一函数分量的优化函数驱动,所述至少一个第一函数分量基于对所述目标图像与通过将所述至少一个中间配准变换应用于所述源图像而能够获得的变换图像之间的相似性进行量度的相似性量度,
体积监测器(IVM),其被配置为在迭代周期中监测当在给定的迭代周期中应用各自的中间配准变换时所述源图像中的图像位置周围的邻域的体积变化,所述邻域是针对所述源图像中的所述图像位置预定义的并且在迭代周期中被监测;
体积评估器(IVE),其被配置为针对接受条件来评估所述体积变化;以及
校正器(C),其被配置为当由所述体积评估器(IVE)发现所述变化违反所述接受条件时,在后续的迭代周期中减轻或至少相对减轻所述第一函数分量对所述图像位置的影响。
2.根据权利要求1所述的图像处理系统,其中,在所述迭代的剩余部分期间,所述第一函数分量对所述图像位置保持减轻。
3.一种图像处理方法,包括:
接收两幅或更多幅输入图像以进行配准,所述两幅或更多幅输入图像包括源图像和目标图像;
通过迭代通过至少一个中间配准变换T(N-j)来迭代地求解最终配准变换T(N),其中,所述迭代地求解由包括至少一个第一函数分量的优化函数驱动,所述至少一个第一函数分量基于对所述目标图像与通过将所述至少一个中间配准变换应用于所述源图像而能够获得的变换图像之间的相似性进行量度的相似性量度,
在迭代周期中监测当在给定的迭代周期中应用各自的中间配准变换时所述源图像中的图像位置周围的邻域的体积变化,所述邻域是针对所述源图像中的所述图像位置预定义的并且在迭代周期中被监测;
针对接受条件来评估所述体积变化;并且
当发现所述变化违反所述接受条件时,在后续的迭代周期中减轻或者至少相对减轻所述第一函数分量对所述图像位置的影响。
4.一种计算机可读介质,具有存储于其上的用于控制根据权利要求1-2中的任一项所述的系统的计算机程序单元,所述计算机程序单元当由处理单元(PU)运行时适于执行根据权利要求3所述的方法的步骤。
CN201680081685.4A 2015-12-15 2016-12-07 图像处理系统和方法 Active CN108701360B (zh)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
EP15200216.8 2015-12-15
EP15200216 2015-12-15
PCT/EP2016/079969 WO2017102468A1 (en) 2015-12-15 2016-12-07 Image processing systems and methods

Publications (2)

Publication Number Publication Date
CN108701360A CN108701360A (zh) 2018-10-23
CN108701360B true CN108701360B (zh) 2023-05-26

Family

ID=55024807

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201680081685.4A Active CN108701360B (zh) 2015-12-15 2016-12-07 图像处理系统和方法

Country Status (4)

Country Link
US (1) US11538176B2 (zh)
EP (1) EP3391337A1 (zh)
CN (1) CN108701360B (zh)
WO (1) WO2017102468A1 (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109964249A (zh) * 2016-09-21 2019-07-02 皇家飞利浦有限公司 用于对身体部分的自适应轮廓勾画的装置
US11080858B2 (en) * 2016-12-12 2021-08-03 Koninklijke Philips N.V. Method and apparatus for segmenting two dimensional images of an anatomical structure
US20190045213A1 (en) * 2017-08-03 2019-02-07 Intel Corporation Reference frame reprojection for improved video coding
CN111815684B (zh) * 2020-06-12 2022-08-02 武汉中海庭数据技术有限公司 基于统一残差模型的空间多元特征配准优化方法及装置
CN113327274B (zh) * 2021-04-15 2024-01-30 中国科学院苏州生物医学工程技术研究所 集成分割功能的肺ct图像配准方法及系统

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101038669A (zh) * 2007-04-12 2007-09-19 上海交通大学 全局异常信号环境下基于联合显著图的鲁棒图像配准方法
CN102135606A (zh) * 2010-12-13 2011-07-27 电子科技大学 一种基于knn分类算法的mr图像灰度不均匀性校正分割方法
EP2369551A1 (en) * 2010-03-25 2011-09-28 Emory University Imaging system and method
WO2014001959A1 (en) * 2012-06-27 2014-01-03 Koninklijke Philips N.V. Image quality driven non-rigid image registration
CN103871056A (zh) * 2014-03-11 2014-06-18 南京信息工程大学 基于各向异性光流场与去偏移场的脑部mr图像配准方法
CN103871063A (zh) * 2014-03-19 2014-06-18 中国科学院自动化研究所 一种基于点集匹配的图像配准方法
CN104091337A (zh) * 2014-07-11 2014-10-08 北京工业大学 一种基于PCA及微分同胚Demons的变形医学图像配准方法
CN104732529A (zh) * 2015-03-05 2015-06-24 北京空间机电研究所 一种遥感图像形状特征配准方法
CN104835112A (zh) * 2015-05-07 2015-08-12 厦门大学 一种肝脏多相期ct图像融合方法

Family Cites Families (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB0227887D0 (en) * 2002-11-29 2003-01-08 Mirada Solutions Ltd Improvements in or relating to image registration
US7155047B2 (en) 2002-12-20 2006-12-26 General Electric Company Methods and apparatus for assessing image quality
JP4705104B2 (ja) 2004-08-09 2011-06-22 ブラッコ・シュイス・ソシエテ・アノニム 複数のマスクに基づく医療画像処理のためのイメージ登録方法および装置
WO2006051488A1 (en) 2004-11-10 2006-05-18 Koninklijke Philips Electronics N.V. System and method for registration of medical images
US8126230B2 (en) * 2006-04-20 2012-02-28 Koninklijke Philips Electronics N.V. Method of motion correction for dynamic volume alignment without timing restrictions
US8358818B2 (en) * 2006-11-16 2013-01-22 Vanderbilt University Apparatus and methods of compensating for organ deformation, registration of internal structures to images, and applications of same
US8098911B2 (en) * 2006-12-05 2012-01-17 Siemens Aktiengesellschaft Method and system for registration of contrast-enhanced images with volume-preserving constraint
US20090074276A1 (en) 2007-09-19 2009-03-19 The University Of Chicago Voxel Matching Technique for Removal of Artifacts in Medical Subtraction Images
US9008462B2 (en) * 2011-02-17 2015-04-14 The Johns Hopkins University Methods and systems for registration of radiological images
KR101141312B1 (ko) 2011-05-16 2012-05-04 동국대학교 산학협력단 영상의 융합 기법을 이용한 의료용 혈관영상 처리방법
US8655040B2 (en) * 2012-03-01 2014-02-18 Empire Technology Development Llc Integrated image registration and motion estimation for medical imaging applications
CN102722890B (zh) * 2012-06-07 2014-09-10 内蒙古科技大学 基于光流场模型的非刚性心脏图像分级配准方法
CN104395933B (zh) * 2012-06-27 2017-05-17 皇家飞利浦有限公司 运动参数估计
JP5889265B2 (ja) 2013-04-22 2016-03-22 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 画像処理方法および装置並びにプログラム
JP6145178B2 (ja) * 2013-10-18 2017-06-07 コーニンクレッカ フィリップス エヌ ヴェKoninklijke Philips N.V. 医療画像の位置合わせ
CN103632338B (zh) 2013-12-05 2016-08-31 鲁东大学 一种基于匹配曲线特征的图像配准评估法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101038669A (zh) * 2007-04-12 2007-09-19 上海交通大学 全局异常信号环境下基于联合显著图的鲁棒图像配准方法
EP2369551A1 (en) * 2010-03-25 2011-09-28 Emory University Imaging system and method
CN102135606A (zh) * 2010-12-13 2011-07-27 电子科技大学 一种基于knn分类算法的mr图像灰度不均匀性校正分割方法
WO2014001959A1 (en) * 2012-06-27 2014-01-03 Koninklijke Philips N.V. Image quality driven non-rigid image registration
CN103871056A (zh) * 2014-03-11 2014-06-18 南京信息工程大学 基于各向异性光流场与去偏移场的脑部mr图像配准方法
CN103871063A (zh) * 2014-03-19 2014-06-18 中国科学院自动化研究所 一种基于点集匹配的图像配准方法
CN104091337A (zh) * 2014-07-11 2014-10-08 北京工业大学 一种基于PCA及微分同胚Demons的变形医学图像配准方法
CN104732529A (zh) * 2015-03-05 2015-06-24 北京空间机电研究所 一种遥感图像形状特征配准方法
CN104835112A (zh) * 2015-05-07 2015-08-12 厦门大学 一种肝脏多相期ct图像融合方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于自由形变的3D非线性医学图像配准;李文龙等;《中国医学影像技术》;20111220(第12期);第2536-2540页 *

Also Published As

Publication number Publication date
US11538176B2 (en) 2022-12-27
CN108701360A (zh) 2018-10-23
US20180336689A1 (en) 2018-11-22
EP3391337A1 (en) 2018-10-24
WO2017102468A1 (en) 2017-06-22

Similar Documents

Publication Publication Date Title
CN108701360B (zh) 图像处理系统和方法
Vishnevskiy et al. Isotropic total variation regularization of displacements in parametric image registration
US8326086B2 (en) Elastic image registration
US20140323845A1 (en) Automated 3-d orthopedic assessments
US10867423B2 (en) Deformation field calculation apparatus, method, and computer readable storage medium
Bağcı et al. The role of intensity standardization in medical image registration
US20170206670A1 (en) Image processing apparatus, image processing method, and non-transitory computer-readable storage medium
US9224204B2 (en) Method and apparatus for registration of multimodal imaging data using constraints
JP2005197792A (ja) 画像処理方法、画像処理装置、プログラム、記憶媒体及び画像処理システム
US7616783B2 (en) System and method for quantifying motion artifacts in perfusion image sequences
US20190035094A1 (en) Image processing apparatus, image processing method, image processing system, and program
Tenbrinck et al. Histogram-based optical flow for motion estimation in ultrasound imaging
KR102036834B1 (ko) 이미지 처리 방법
CN115861656A (zh) 用于自动处理医学图像以输出警报的方法、设备和系统
US10810717B2 (en) Image processing apparatus, image processing method, and image processing system
JP2016002251A (ja) 画像処理装置、画像処理方法、およびプログラム
US8867809B2 (en) Image processing method
CN113554647B (zh) 医学图像的配准方法和装置
US11138736B2 (en) Information processing apparatus and information processing method
JP2023510246A (ja) 医用画像における腫瘍体積の変化の測定
RU2478337C2 (ru) Способ определения контура сердца на флюорографических снимках
JP5706933B2 (ja) 処理装置および処理方法、プログラム
US11704795B2 (en) Quality-driven image processing
Abidi Deformable image registration methodologies for thoracic images seqouences
CN116612163A (zh) 用于将解剖图像配准到功能图像的深度学习

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant