CN103814384B - 基于图像的跟踪 - Google Patents

基于图像的跟踪 Download PDF

Info

Publication number
CN103814384B
CN103814384B CN201280028147.0A CN201280028147A CN103814384B CN 103814384 B CN103814384 B CN 103814384B CN 201280028147 A CN201280028147 A CN 201280028147A CN 103814384 B CN103814384 B CN 103814384B
Authority
CN
China
Prior art keywords
image
matrix
low
tissue
rank
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
CN201280028147.0A
Other languages
English (en)
Other versions
CN103814384A (zh
Inventor
余维川
梁天柱
周晓巍
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.)
Hong Kong University of Science and Technology HKUST
Original Assignee
Hong Kong University of Science and Technology HKUST
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 Hong Kong University of Science and Technology HKUST filed Critical Hong Kong University of Science and Technology HKUST
Publication of CN103814384A publication Critical patent/CN103814384A/zh
Application granted granted Critical
Publication of CN103814384B publication Critical patent/CN103814384B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/56—Details of data transmission or power supply
    • A61B8/565—Details of data transmission or power supply involving data transmission via a network
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/58—Testing, adjusting or calibrating the diagnostic device
    • A61B8/587—Calibration phantoms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/18—Image warping, e.g. rearranging pixels individually
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/246—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/246—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
    • G06T7/251—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments involving models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/277—Analysis of motion involving stochastic approaches, e.g. using Kalman filters
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/08—Detecting organic movements or changes, e.g. tumours, cysts, swellings
    • A61B8/0883—Detecting organic movements or changes, e.g. tumours, cysts, swellings for diagnosis of the heart
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/48—Diagnostic techniques
    • A61B8/485—Diagnostic techniques involving measuring strain or elastic properties
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10016—Video; Image sequence
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10132—Ultrasound image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30048—Heart; Cardiac

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Multimedia (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Medical Informatics (AREA)
  • Public Health (AREA)
  • Biomedical Technology (AREA)
  • Animal Behavior & Ethology (AREA)
  • Molecular Biology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Veterinary Medicine (AREA)
  • Surgery (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Data Mining & Analysis (AREA)
  • Quality & Reliability (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Artificial Intelligence (AREA)
  • Image Analysis (AREA)
  • Image Processing (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本文中描述了便于图像运动分析的系统和方法。根据第一图像运动分析技术，根据局部仿射模型来变形第一图像。比较变形后的第一图像和第二图像，寻找变形后的第一图像与第二图像之间匹配的图案。根据所述图案估计运动参数的值。根据第二图像运动分析技术，把图像序列转换成输入矩阵。输入矩阵的列对应于与图像序列相关的向量化的图像。利用低秩矩阵逼近输入矩阵，低秩矩阵的秩小于输入矩阵。检测低秩矩阵的一个或者多个离异值。

Description

基于图像的跟踪
相关申请的交叉引用
本申请要求于2011年6月9日提交的题为“An Affine Warping Method to Solvethe Feature Motion Decorrelation Problem in Ultrasound Image Based Tracking”的美国临时专利申请编号61/457,813的优先权。本申请还要求于2012年1月12日提交的题为“Method of Image Analysis Based on an Ultrasound Image Sequence forAutomatic Mitral Leaflet Tracking”的美国临时专利申请编号61/631,815的优先权。
技术领域
本公开主要涉及到图像分析，更具体地涉及图像运动分析的简易化。
背景技术
软组织的力学性质与各种病症相关。例如，癌组织往往比非癌组织更硬，肝硬化的肝组织比正常肝组织更硬，缺血性心脏往往表现出心肌组织的异常收缩/松弛。力学性质和病症之间的相关性，可以帮助病症诊断。
弹性成像是一种无创性的方法，用于测量软组织的力学性质，该方法有利于各种病症的诊断。在对软组织施加机械力的之前和之后的图像被采集。这些图像有相关性，组织的运动或变形可以通过运动跟踪算法进行推算。
一个可用于弹性成像的运动跟踪算法的例子是散斑跟踪。在成像模式为超声成像时，当发射的超声波被组织反射并相互干涉时形成散斑。散斑追踪是基于散斑图案在组织运动前后保持不变的假设。散斑图案保持不变的假设被用作运动跟踪算法的基础，但是，散斑图案在组织变形后是会发生变化的。由于组织变形后散斑图案的变化(“特征运动去相关”)，散斑跟踪算法不能揭示真正的组织运动。因此，散斑跟踪算法不能准确估计组织的运动以及不能正确地推断出力学性质。
运动跟踪与几何形状也可用于各种病症的诊断。例如，二尖瓣是位于心脏的左心房和左心室之间的薄瓣膜结构，用于控制血流方向。二尖瓣相关的疾病，如二尖瓣反流，是最常见的心脏瓣膜病。对心脏进行成像并且获取患者的具体二尖瓣几何特征，以及跟踪二尖瓣的运动，可以便于诊断心脏瓣膜病和/或协助瓣膜修复的手术治疗。
在各种成像模式中，实时三维超声心动图提供了一种无创性的方法来模拟二尖瓣的三维几何和捕捉它的快速运动。为了生成一个完整的二尖瓣模型，需要在整个超声心动图序列中跟踪二尖瓣瓣叶。为了跟踪诸如二尖瓣瓣叶之类的对象，在整个图像序列中该对象被定位并分割。二尖瓣瓣叶的定位和分割可以手动完成，但手动划分是劳动密集的，并且容易产生很大的方差。尤其对于三维图像方差很大，因为操作者每次仅可以显示和处理体数据的二维投影或截面。跟踪算法能够在自然图像分析中跟踪一定物体。然而，跟踪二尖瓣瓣叶是非常困难的，并且现在没有方法用于这样的跟踪。因为缺乏可靠的特征以及瓣膜快速且不规则的运动，很难跟踪二尖瓣。
以上背景描述仅仅是为了提供关于运动跟踪算法的背景信息的概述，并非是详尽的。其他的背景在回顾以下详细描述的非限制性实施例的一个或多个后有可能进一步被显现。
发明内容
下面提供了说明书的简化总结，以便提供对说明书的某些方面的基本了解。本概要并非说明书的广泛概述。它的目的既不是确定说明书的主要或关键要素，也不是划定说明书的具体实施例的任何范围，或者权利要求的任何范围。其唯一目的是以简化的形式提供说明书的一些构思，作为稍后提出的更详细说明的开头。
根据一个或多个实施例和对应的公开，与图像运动分析技术的简易化相关地描述了各种非限制性方面。非限制性实施例通过特征运动去相关来便于确定力学性质。另外的非限制性实施例通过不依赖于用户交互和训练数据的分割来便于物体跟踪。
根据一个非限制性的方面，描述了便于通过弹性成像确定组织力学性质的系统和方法。根据非限制性实施例，描述了包括图像形变组件和估计组件的系统。图像形变组件配置成根据局部仿射模型来变换图像从而产生变换后的第一图像。估计组件配置成根据预设的匹配测度来寻找变换后的第一图像与第二图像之间匹配的图案。估计组件还配置成基于该图案估计运动参数。所述系统包括存储各种组件的存储器和执行一个或多个组件的处理器。
在另一非限制性实施例中，描述了一种方法，该方法包括通过包括处理器的系统根据局部仿射模型而变换图像来产生变换后的第一图像。所述方法还包括通过系统根据预设的匹配测度来寻找变换后的第一图像与第二图像之间匹配的图案，以及通过系统根据该图案估计运动参数。
在另一个非限制性实施例中，描述了一种计算机可读存储介质，其上存储了计算机可执行指令，计算机可执行命令响应于执行而使得包括处理器的设备进行操作。该操作至少包括：接收第一图像和第二图像；初始化运动参数；根据局部仿射模型来变换第一图像；根据预设的匹配测度来寻找变换后的第一图像与第二图像之间匹配的图案；以及根据该图案估计运动参数。
根据另一非限制性方面，描述了不依赖于另外的用户交互和训练数据的便于自动物体跟踪的系统和方法。根据非限制性实施例，描述了一种系统，其包括转换组件，转换组件便于把图像序列转换成输入矩阵，输入矩阵的列对应于图像序列中的向量化的图像。系统还包含近似组件，其配置成利用低秩矩阵逼近输入矩阵，低秩矩阵的秩小于输入矩阵的秩。系统还包含检测组件，其配置成检测低秩矩阵的一个或者多个离群值。所述系统包含存储各种组件的存储器和便于执行一个或多个组件的处理器。
在另一非限制性实施例中，描述了一种方法，该方法包括通过包括处理器的系统把图像序列转换成输入矩阵。输入矩阵的列对应于图像序列中的向量化的图像。该方法还包括通过系统利用一个低秩矩阵逼近输入矩阵，低秩矩阵的秩小于输入矩阵的秩。该方法还包括通过系统检测低秩矩阵的一个或者多个离群值。
在另一非限制性实施例中，描述了一种计算机可读存储介质，其上存储了计算机可执行指令，计算机可执行指令响应于执行而使得包括处理器的设备进行操作。操作至少包括：把图像序列转换成输入矩阵，输入矩阵的列对应于图像序列中的向量化的图像；利用低秩矩阵逼近输入矩阵，低秩矩阵的秩小于输入矩阵的秩；以及检测低秩矩阵的一个或者多个离群值。
下面的描述和附图阐述了说明书的某些示例性方面。然而，这些方面表示可以采用说明书的各种实施例的各种方式的仅一些方式。说明书的其他方面在结合附图考虑时将从说明书的以下详细说明中变得显而易见。
附图说明
以下将参照附图对各种方面和实施例作进一步说明，其中相同的参考标记始终指代相同部分，附图中：
图1示出了根据本公开的实施例的通过特征运动去相关来便于图像分析的示例非限制性系统；
图2示出了根据本公开的实施例的便于补偿特征运动去相关的示例非限制性系统；
图3示出了根据本公开的实施例的优化运动参数以便于补偿特征运动去相关的示例非限制性系统；
图4是根据本公开的实施例的通过特征运动去相关来图像分析的方法的示例非限制性处理流程图；
图5是根据本公开的实施例的便于补偿特征运动去相关的方法的示例非限制性处理流程图；
图6是根据本公开的实施例的优化运动参数以便于补偿特征运动去相关的方法的示例非限制性处理流程图；
图7是根据本公开的实施例的三维超声仿真实验中的示例非限制性输入图像；
图8至图10是根据本公开的实施例的三维超声仿真实验中平均相关系数关于运动的示例非限制性变化曲线；
图11至图13是根据本公开的实施例的三维超声仿真实验中均方误差关于运动的示例非限制性曲线；
图14是根据本公开的实施例的二维仿组织胶状模型的示例非限制性输入图像；
图15是根据本公开的实施例的二维仿组织胶状模型在2％压缩下耦合滤波方法得到的相关系数；
图16是根据本公开的实施例的二维仿组织胶状模型在2％压缩下仿射变形方法得到的相关系数；
图17是根据本公开的实施例的二维仿组织胶状模型在2％压缩下压扩方法得到的相关系数；
图18是根据本公开的实施例的二维仿组织胶状模型在2％压缩下耦合滤波方法估计的轴向应变；
图19是根据本公开的实施例的二维仿组织胶状模型在2％压缩下仿射变形方法估计的轴向应变；
图20是根据本公开的实施例的二维仿组织胶状模型在2％压缩下压扩方法估计的轴向应变；
图21是根据本公开的实施例的二维仿组织胶状模型在5％压缩下耦合滤波方法得到的相关系数；
图22是根据本公开的实施例的二维仿组织胶状模型在5％压缩下仿射变形方法得到的相关系数；
图23是根据本公开的实施例的二维仿组织胶状模型在5％压缩下压扩方法得到的相关系数；
图24是根据本公开的实施例的二维仿组织胶状模型在5％压缩下耦合滤波方法估计的轴向应变；
图25是根据本公开的实施例的二维仿组织胶状模型在5％压缩下仿射变形方法估计的轴向应变；
图26是根据本公开的实施例的二维仿组织胶状模型在5％压缩下压扩方法估计的轴向应变；
图27至图29是根据本公开的实施例的三维超声仿真实验中平均相关系数关于运动的示例非限制性变化曲线；
图30是示出根据本公开的实施例的横向应变运动估计的均方误差的示例非限制性曲线；
图31是示出根据本公开的实施例的轴向-横向切应变运动估计的均方误差的示例非限制性曲线；
图32是示出根据本公开的实施例的横向旋转运动估计的均方误差的示例非限制性曲线；
图33是根据本公开的实施例的二维仿组织模型在2％压缩下仿射变形方法估计的相关系数；
图34是根据本公开的实施例的二维仿组织模型在2％压缩下仿射变形方法估计的轴向应变；
图35是根据本公开的实施例的二维仿组织模型在5％压缩下仿射变形方法估计的相关系数；
图36是根据本公开的实施例的二维仿组织模型在5％压缩下仿射变形方法估计的轴向应变；
图37示出了根据本公开的实施例的便于分割并跟踪图像序列内物体特征的示例非限制性系统；
图38示出了根据本公开的实施例的便于把图像序列转换成一个输入矩阵的示例非限制性系统；
图39示出了根据本公开的实施例的便于把输入矩阵建模成一个低秩矩阵和一组离群值的示例非限制系统；
图40是根据本公开的实施例的便于分析一系列图像的方法的示例非限制性处理流程图；
图41是根据本公开的实施例的便于检测对应于图像中某一特征的离群值的方法的示例非限制性处理流程图；
图42是根据本公开的实施例的通过迭代地使功函数最小化来便于分割并跟踪图像特征的方法的示例非限制性处理流程图；
图43是根据本公开的实施例的示出矩阵谱分析过程的可行性的示例非限制性示图；
图44是根据本公开的实施例的通过对二维超声心动图序列的矩阵谱分析来跟踪二尖瓣的示例非限制性结果；
图45是根据本公开的实施例的通过对三维超声心动图序列的矩阵谱分析来跟踪二尖瓣的示例非限制性结果；
图46示出了可以实现本文所述的各种实施例的示例计算环境；以及
图47示出了可以实现本文所述的各种实施例的计算机网络的示例。
具体实施方式
参考附图来描述本公开的各个方面或特征，附图中相同的参考编号始终用来指代相同元件。在本说明书中，为了提供对本公开的全面理解，阐述了各种具体细节。然而，应当理解的是，在没有这些具体细节的情况下，或者利用其他方法、组件、材料等，也可以实施本发明的某些方面。另外，以框图形式示出了公知结构和装置，以便于描述和说明各种实施例。
根据本公开中描述的一个或多个实施例，结合图像运动分析技术描述了各种非限制性方面。非限制性实施例便于通过特征运动去相关来确定力学性质。另外的非限制性实施例便于通过无需用户交互或训练数据的分割来物体跟踪。
特征运动去相关
根据本公开中描述的一个或多个实施例，与通过特征运动去相关来确定物体力学性质相关地描述了各种非限制性方面。特征运动去相关是通过仿射变形处理(affinewarping process)而解决的。根据局部仿射模型来变形物体的第一图像。将变形后的图像和物体的第二图像进行比较，从而得到在变形图像与第二图像之间匹配的图案。基于所述图案，估计运动参数的值。运动参数可以用来确定力学性质。
现参考附图，首先图1示出了根据本公开的实施例的通过特征运动去相关来便于图像分析的示例非限制性系统100。系统100可以便于通过弹性成像来确定物体的力学性质。弹性成像便于基于对物体施加机械力之前和之后获得的物体图像来确定力学性质。这些图像可以是超声图像、计算机断层扫描图像、磁共振成像的图像、或任何其他类型的图像或系列的图像。这些图像可以是二维图像或者三维图像。
系统100包括图像变形组件102和估计组件104。该系统还包括存储这些组件的存储器和便于执行一个或多个组件的处理器。
系统100通过图像变形组件102和估计组件104来解决特征运动去相关问题。在实施例中，特征运动去相关问题存在于利用超声的散斑跟踪来便于确定物体的运动(例如组织运动)的特性。为了解决特征运动去相关的问题，不同的模型可以用于模拟斑点变化的原因。
这样的变化可能是由于组织及其运动的复杂性而导致的。组织可能有不同的形状和力学性质(例如，通常癌组织比正常组织更硬，而来自不同器官的组织也具有不同力学性质)。因此，机械力下的各组织之间的相互作用是非常复杂的。超声图像特性也可能会引起图像变化。超声图像是当超声波反射后互相干涉时形成的。因此，当组织经历复杂的运动时，超声波会引起不同的干涉，产生了更多的变种散斑图案。
一个可以用来解决特征运动去相关问题的模型是利用压扩方法的模型。压扩方法试图将由于组织形变引起的图像变化建模成图像的压缩和平移。这种方法提供了组织形变的很好的近似，而且是描述并补偿由于组织形变引起图像变化的最好的方法之一，但是在组织形变很大时图像变化的补偿将失效。
另一个建模技术是耦合滤波方法。在耦合滤波方法中，复杂的组织运动和超声图像特性都被建模。局部仿射运动模型用来建模复杂组织运动，线性卷积模型用来描述超声波在成像过程中的干涉。
局部仿射模型用于模拟复杂的组织运动。例如，将组织的刚性/非刚性运动近似为局部区域中的仿射运动(例如，压缩或扩大)，而局部区域是在运动估计之前被预定义或者在运动估计期间被适应性调整。可以在时域或者频域内估计这样的局部运动。假设组织运动在局部区域中是仿射的，则复杂的组织运动问题可以通过以下局部仿射模型来排除：
xn＝MXn+T
这里Xn和xn分别是第n个组织散射点(根据实施例，是超声波反射组织单元)在组织运动前后的位置。M是用于仿射运动的矩阵，T是用于组织平移的矢量。在这个局部仿射模型中，可以建模各种组织运动，比如平移、压缩、扩大、剪切和旋转。
局部仿射模型中的局部区域可以被预设为足够小，使得复杂的组织运动在该区域内可以被近似成仿射运动，同时这个区域应该足够大，使得足够的散斑图案存在以识别匹配的图案。可替代地，可以在最优运动参数(在该情况下，最优的M和T)的估计期间适当调节局部区域。例如，局部区域起初可以是大的，使得产生最优M和T的粗略估计。然后，可以减小该局部区域，并产生精细估计。
基于局部仿射模型，组织的位移也为仿射的：
dn＝xn-Xn＝(M-I)Xn+T，
这里dn是第n个组织散射点的位移，I是单位矩阵，其大小与M相等。如果最优的运动参数M和T被准确估计,局部区域中的组织运动就可以被揭示。传统上，如果对dn关于Xn求导，应变值(例如，正应变和切应变)和旋转也可以被估计。可替代地，因为dn是Xn的仿射函数，应变值和旋转就是矩阵M的各项的组合，如无穷小应变理论中所述。以这样的方式，可以从最优运动参数中直接推理出应变和旋转。
线性卷积模型用于描述成像过程中的超声波的干涉。下面的线性卷积模型可以很好的描述散斑图案的形成：
这里Tn(X；Xn)是第n个组织散射点，H(X)表示超声系统的点扩散函数(PSF)，N是组织散射点总数，I(X；Xn)是射频超声图像。在原点处只有一个组织散射点存在时超声系统的点扩散函数可以被解释成线性系统响应函数。卷积模型假设每个组织散射点的贡献都是独立的，并且所有组织散射点加起来之后就产生最后的超声图像。
在如上所述的线性卷积模型中，每个组织散射点可以模拟成一个迪拉克(Dirac)函数：
Tn(X；Xn)＝anδ(X-Xn),
这里Xn是第n个组织散射点的位置，an是第n个组织散射点的反射系数，其表示有多少百分比的超声波被组织散射点反射。点扩散函数H(X)通常不具有解析表达式，其取决于超声探头(比如超声探头的形状和元件排列、发射的超声波的中心频率和带宽、组织的散射能力等等)。在一些文献中，点扩散函数被近似为Gabor函数。可替代地，通过基于不同成像设置求解声波方程来数值地计算出更准确的点扩散函数。
基于局部仿射运动模型和线性卷积模型，耦合滤波方法用如下公式描述组织运动读取前后的超声图像：
由于组织散射点位置的改变，超声图像中产生的干涉图案在组织运动前后可以不同。这将导致特征运动去相关问题。如果去相关可以得到补偿，则组织运动可以被可靠地估计。
通过耦合滤波方法而应用的补偿是基于组织运动前后获得的射频图像的如下等式关系：
这里I(MX+T；xn)表示组织运动后的超声图像的仿射变形的结果(即，应用仿射变形步骤之后的第二图像)，I(MX+T；xn)*H(X)表示组织运动后的超声图像的耦合滤波的结果(即，应用仿射变形步骤和耦合滤波步骤之后的第二图像)，I(X；Xn)*H(MX)表示组织运动前的超声图像的耦合滤波的结果(即，应用耦合滤波步骤之后的第一图像)。基于以上关系，可以得到满足以下等式的最优运动参数M和T：
然后从最优运动参数中可以直接或间接得到应变和旋转。
耦合滤波方法能够准确的估计组织运动，即使在组织形变很大的情况下。然而，为了补偿散斑图案的变化，在寻找到匹配图案之前，耦合滤波方法涉及仿射变形步骤和耦合滤波步骤两者，其很消耗时间。而且，耦合滤波方法只能用于在临床实践中通常很难得到的射频(RF)超声图像。另外，耦合滤波方法假设点扩散函数(PSF)在线性卷积模型中是已知的并且是不变量，在实际应用中一般不会这样。因此，需要一种能够应用于B超图像解决特征运动去相关问题的简单且鲁棒的方法。
系统100应用了耦合滤波的近似方法—仿射变形。系统100可以对两张二维或者三维、B超或者射频超声图像进行仿射变形，这两张图分别对应组织运动之前和之后的超声图像。应当注意的是，系统100也可以对超声图像序列进行仿射变形。例如，可以针对每对相邻帧估计小组织运动，并积累小组织运动以产生大变形的估计。可替代地，也可以直接使用第一帧作为基准帧来估计每个帧的组织运动。
系统100便于通过图像变形组件102和估计组件104来应用仿射变形。图像变形组件102可以根据局部仿射模型来变形图像，以产生变形后的第一图像。估计组件104配置成可以根据预设的匹配测度来寻找变形后的第一图像与第二图像之间匹配的图案。估计组件还配置成根据所述图案估计运动参数。
系统100使用的仿射变形不同于其他可以减弱特征运动去相关问题的方法。最相关的方法是耦合滤波方法。在耦合滤波方法中，局部仿射运动模型用来建模复杂组织运动，线性卷积模型用来描述当超声波被组织反射后干涉时超声图像的形成。相比之下，系统100使用的仿射变形方法只用了局部仿射模型用来描述组织运动。超声图像形成期间的超声波的干涉被忽略了。
系统100采用的仿射变形与耦合滤波方法的不同之处在于对散斑图案变化的补偿方式不同。在耦合滤波中，对两张超声图像中的一个进行了仿射变形，然后对两个超声图像应用两个不同的滤波器。相比之下，系统100使用仿射变形对一个超声图像仅进行仿射变形而没有进一步进行滤波步骤，这样减少了计算量。
系统100使用的仿射变形方法使用了对耦合滤波方法的如下近似：
I(MX+T；xn)＝I(X；Xn)
上述近似仅是近似关系。耦合滤波方法中的卷积操作被省略了。没有补偿掉的散斑图案变化可来源于点扩散函数项的变化。如果来自组织散射点项的图像变化远大于来自点扩散函数项的图案变化，以上近似能够提供很好的耦合滤波方法的近似。
仿射变形方法是基于I(MX+T；xn)＝I(X；Xn)。如在耦合滤波方法中，在寻找匹配图案之前向两个图像中的一个图像应用仿射变形步骤。没有应用滤波步骤。相应地，图像变形组件102根据局部仿射模型来变形图像，以产生变形后的第一图像，而不应用滤波方法，然后估计组件104根据预设的匹配测度来寻找变形后的第一图像与第二图像之间匹配的图案。
仿射变形方法相对于耦合滤波方法的另外一个好处在于其既不假设一致的点扩散函数也不需要关于点扩散函数的知识。该仿射变形方法不限制于射频图像，也能应用于B超图像。因为卷积运算被去掉，仿射变形方法的计算(比如通过处理器108计算)速度要快于耦合滤波方法的计算速度。为了进一步提高寻找最优运动参数的速度，也可以利用具有不同的处理器108(比如多核CPU、工作站、集群、GPU和GPU集群)的其他计算设备。
参考图2，示出了根据本公开的实施例的示例非限制性系统200，其便于通过仿射变形来补偿特征运动去相关。第一图像202和第二图像204被输入到系统200。在实施例中，第一图像202和第二图像204也可以是图像组。第一图像202和第二图像204可以对应于施加机械力前后得到的两张图像。需要说明的是，对于第一图像和第二图像在施加机械力前后的定义并不重要。比如第一图像也可以是施加机械力之后的图像，而第二图像也可以是施加机械力之前的图像。
图像变形组件102根据局部仿射模型来变形第一图像202，从而得到变形后的第一图像206(或者变形后的第一组图像)。变形后的第一图像206和第二图像204被输入到估计组件104。估计组件104(匹配组件208)根据预设的匹配测度来寻找变形后的第一图像206与第二图像204之间匹配的图案，并且(参数组件210)根据该图案估计运动参数。估计组件104可以估计组织运动，并且可以通过将反向运动场施加到该图像而恢复组织运动，来消除组织运动引起的组织图案变化。
匹配组件208可以使用一个或者多个匹配测度来寻找匹配图案。匹配测度可以包括变形后的第一图像206和第二图像204之间的相关系数、误差平方和、绝对误差和等等。
参数组件210可以便于估计运动参数。在实施例中，运动参数可以是下面等式中的M：
I(MX+T；xn)＝I(X；Xn)，
在另一实施例中，估计的运动参数也可以包含T。
参数组件210也可以包含便于根据一个或者多个限制条件来进行运动参数估计的限制组件。限制组件可以利用不同的数学模型来更好的估计组织运动，并且补偿由组织运动引起的散斑图案变化。
在实施例中，限制条件可以与被成像的物体的特性有关。限制条件可以是平滑度约束。当物体是生物组织时，平滑度约束确保相邻组织的运动是相似的。限制条件也可以是变形约束。当物体是生物组织时，变形约束假设可以从控制点的运动中插入其他位置处的组织运动。当物体是生物组织时，变形约束也可以与生物组织的不可压缩性相关。
现在参考图3，示出了根据本公开的实施例的示例非限制性系统300，其可以优化运动参数以便于补偿特征运动去相关。系统300接收第一图像或者第一组图像202以及第二图像或者第二组图像204的输入。第一图像202和第一图像204可以分别在施加机械力的前后获得。当物体是生物组织时，这些图像可以分别在组织运动前后采集得到。
系统300的输出是最优运动参数310或者一组最优运动参数(比如最优的M和T)。最优运动参数310可以对散斑图案变化进行最优补偿。正应变和切应变以及旋转角度可以根据无穷小应变理论而从最优运动参数M中估计出。在此情况下，由于应变和旋转是M的各项的简单组合，因此可以直接估计出组织运动。在射频图像中，因为有足够的高频成分作为独特特征，所以这样的方法可以得到很好的结果。但是在B超图像的运动估计中，更推荐传统的间接方法，即通过最优运动参数M和T来估计位移。应变和/或旋转可以通过对估计的位移求偏导得到。
初始化组件302对运动参数304a和304b进行初始化。系统300遍历一系列候选运动参数，其是通过由图像变形组件102变形第一图像202来得到变形后的第一图像206，并且把变形后的第一图像206和第二图像204输入到估计组件104。估计组件104利用每一个候选运动参数来补偿散斑图案变化，然后优化组件306利用一些标准决定补偿是不是最优的。如果当前的候选参数实现了最优补偿，该最优运动参数会作为最优运动参数310而被输出。
优化组件306可以使用定量测度来决定最优补偿是否出现。比如，散斑图案变化的补偿在各图案彼此最相似时得到最优。诸如绝对误差和(SAD)、误差平方和(SSD)、和/或相关系数(CC)之类的匹配测度可以用来寻找匹配图案以及用来衡量图案的相似度。
最优参数310的寻找可以通过由优化组件310在利用/不利用多尺度框架的情况下遍历附加限制约束(例如，组织不可压缩性约束、搜索空间约束、诸如时间限制或处理周期限制之类的时间阈值或计算阈值)下所有可能的候选运动参数来实现。而且，可以通过适合的启发式方法(例如，基于梯度的方法、贪婪算法、向一些候选子集施加优先级)的引导，在利用/不利用多尺度框架的情况下，在线确定候选运动参数集。
系统300可以通过对组织运动进行准确和鲁棒的估计来实现仿射变形方法的目标，即使组织运动很大。例如，最优运动参数310被准确的估计并且输出，以描述组织运动。然后，可以从最优运动参数310中直接或间接地产生对应变和旋转的估计。
可以看出，系统300使用的仿射变形方法能够得到跟耦合滤波方法相似的估计精度，即使组织运动大。因为系统300使用的仿射变形方法不需要如耦合滤波中对超声系统的点扩散函数(PSF)的卷积运算，所以需要的时间少于耦合滤波方法。而且，仿射变形方法不假设已知的、不变的点扩散函数(PSF)，并且可以应用于B超图像和射频图像两者，而耦合滤波方法只能应用于射频图像。
图4至图6示出了根据本公开的实施例的方法和/或流程图。为了简化说明，方法被描绘和描述为可以由包括处理器的系统执行的一系列动作。然而，根据本公开的动作可以按各种顺序和/或同时发生，并且可以具有本公开中没有呈现和描述的其他动作。此外，并非需要所有示出的动作来实现根据所公开主题的方法。另外，本领域技术人员将理解并领会到各方法可以通过状态图或事件可替代地被表示为一系列相互关联的状态。此外，还应认识到本说明书所公开的方法能够存储在制品上，以便于将所述方法传输和传递至硬件装置。
参考图4，示出了根据本公开的实施例的便于通过特征运动去相关来进行图像分析的方法400的流程图。方法400是仿射变形方法。在步骤402中，根据局部仿射模型来对图像进行变形，得到变形后的第一图像。在步骤404中，根据预设的匹配测度来寻找变形后的第一图像与第二图像之间匹配的图案，然后由系统基于该图案来估计运动参数。在步骤406中，根据所述图案估计运动参数。
参考图5，示出了根据本公开的实施例的便于补偿特征运动去相关的方法500的示例非限制性处理流程图。方法500是仿射变形方法。在步骤502中，接收第一图像或者第一组图像以及第二图像或者第二组图像。在步骤504中，根据局部仿射模型来变形第一图像，得到变形后的第一图像。在步骤506中，根据预设的匹配测度来寻找变形后的第一图像与第二图像之间匹配的图案。在步骤508中，根据所述图案估计运动参数。
参考图6，示出了根据本公开的实施例的优化运动参数以便于补偿特征运动去相关的方法600的示例非限制性处理流程图。方法600是仿射变形方法。在步骤602中，第一图像被接收。在步骤604中，第二图像被接收。在步骤606中，运动参数被初始化。在步骤608中，将第一图像变形得到变形后的第一图像。在步骤610中，基于变形后的第一图像和第二图像，根据预设的匹配测度来寻找变换后的第一图像与第二图像之间匹配的图案。在步骤612中，基于所述图案估计出运动参数。在步骤614中，判断该运动参数是否最优。在步骤616中，如果该运动参数是最优的，则输出该运动参数以便于确定被成像物体的运动和/或力学性质。在步骤618中，如果该运动参数不是最优的，则调节运动参数并且利用调节后的运动参数的输入来重复步骤608至步骤614继续寻找最优运动参数。
图7至图36示出了仿射变形方法的可行性，并用耦合滤波方法作为比较。仿射变形方法能够如耦合滤波方法那样实现可靠的应变估计，即使其不能实现严格相同的散斑图案。仿射变形方法能够处理及其困难的情况，并且其计算并不像耦合滤波方法那样密集。仿射变形方法能够应用于B超图像，而耦合滤波方法限制于射频图像。图7至图36的实验结果表明，即使当组织变形非常大时，仿射变形方法仍可以在补偿特征运动去相关方面提供跟耦合滤波方法相似的优良性能。
参考图7，示出了根据本公开的实施例的三维超声仿真实验中的示例非限制性输入图像700。三维仿真数据被用来演示仿射变形方法的结果，并且将方式变形方法的结果与耦合滤波方法的结果进行比较。
三维超声图像用Field II程序产生，该程序是广泛应用的超声图像仿真工具。超声探头频率设定为3MHz，探头的相对带宽被设置为0.5。在三维仿真中，使用了一个无切口(kerf)的16×16的多行线性阵列探头元件，同时使用动态聚焦来实现多深度处的高分辨率。图7示出了仿真三维超声图像的一个截面。
通常，在三维仿真中，2000个组织散射点均匀分布在一个10mm×10mm×10mm的目标区域内。分辨单元大小(例如，点扩散函数的半高全宽)约为0.6mm×1.8mm×1.8mm。在每个分辨单元中平均约有4个组织散射点。组织散射点的相关反射系数服从介于0和1之间的高斯分布。轴向、横向和纵向的采样率分别为20、10、10像素每毫米。三维超声图像的尺寸为201×101×01立体像素。
在构造第一超声图像体作为基准之后，根据如下方程模拟组织散射点运动后的位置：
xn＝MXn+T。
Field II程序被再次使用以产生第二超声图像体。三种三维仿射运动被模拟，包括横向压缩和扩张(即，主要沿横向的压缩/扩张)、沿轴向-横向平面的剪切(即，弹性轴垂直于电子束方向)、和横向旋转(即，旋转轴垂直于电子束方向)。选择这三种情况是因为散斑图案变化最大，而还有更多种局部仿射运动带来相似结果并且在此没有描述。在三种情况下，平移向量T被设为0，并且假设目标区域(ROI)的中心为坐标原点。用于以上三种2个方向上的三维仿射运动的M矩阵分别为：
横向压缩/扩张
沿轴向-横向平面的剪切
横向旋转
现在参考图8至图10，示出了根据本公开的实施例的在三维超声仿真实验中平均相关系数相对于运动的示例非限制性变化曲线。图8示出了针对横向压缩和扩张的情况的平均相关系数。图9示出了针对沿轴向横向平面的剪切的情况的平均相关系数。图10示出了针对横向旋转的情况的平均相关系数。在所有三种情况中，比较了仿射变形804、耦合滤波802和无补偿806.
包含25×37×37立体像素的窗口被用来寻找在其他图像体中的匹配图案。相关系数(CC)被用来作为匹配测度，对每一个窗口进行最高相关系数的均值的计算。一共使用了9×9×9＝729个窗口。
对于压缩/扩张情况(图8)，应用的正应变在-20％到20％之间变化，步长为1％。对于剪切情况(图9)，应用的切应变在-10％到10％之间变化，步长为1％。对于旋转情况(图10)，应用的旋转角度在-0.1到0.1弧度之间变化，步长为0.01弧度。三个情况都包括极大运动情形。
在横向压缩/扩张情况下，耦合滤波方法802得到的平均相关系数接近于1，说明耦合滤波方法可以补偿得到一致的散斑图案。在一些极端情况下(例如，针对轴向压缩/扩张和横向压缩/扩张的情况应用的正应变大于15％)，平均相关系数开始下降。而在沿轴向-横向平面的剪切和横向旋转的情况中，对于耦合滤波方法802，平均相关系数逐渐下降。这是由于超声系统的点扩散函数在不同空间位置处不保持相同。因此，由于点扩散函数的变化，可获得其他散斑图案变化。仿射变形方法804得到的平均相关系数逐渐下降，这是因为由于超声图像的干涉性质而形成的散斑图案变化没有得到补偿。没有补偿特征运动去相关806的平均相关系数急剧下降。
参考图11至图13，示出了根据本公开的实施例的在三维超声仿真实验中均方误差相对于运动的示例非限制性变化曲线。图11示出了针对横向压缩和扩张的情况的均方误差的梯度估计。图12示出了针对沿轴向-横向平面的剪切的情况的均方误差的梯度估计。图13示出了针对横向旋转的情况的均方误差的梯度估计。在所有三种情况中，比较了仿射变形804与耦合滤波802和无补偿806。
均方误差被计算为估计的运动参数和真实运动参数之间平方误差的均值。大小为25×37×37的窗口被用于寻找在其他图像体中的匹配图案，而运动估计的均方误差是针对所有窗口计算得到的。一共9×9×9＝729个窗口被计算。
对于压缩/扩张情况(图11)，应用的正应变在-20％到20％之间变化，步长为1％。对于剪切情况(图12)，应用的剪切应变在-10％到10％之间变化，步长为1％。对于旋转情况(图13)，应用的旋转角度在-0.1到0.1弧度之间变化，步长为0.01弧度。三个情况均包括极大运动情况。
在所有三种情况下，耦合滤波方法和仿射变形方法两者得到的均方误差都比较小，除了一些极端情况(比如，针对横向压缩/扩张情况的应用的正应变大于15％、针对沿轴向-横向平面剪切情况的应用的剪切应变大于8％、以及针对横向旋转情况的应用的旋转角度大于0.07弧度)，在这些情况下耦合滤波方法可能会失效。应注意，锯齿状(zig-zag)结构是由于离散最优化和离散输入图像体的性质。若特征运动去相关没有被补偿，则直接组织运动估计失效。
参考图14，示出了根据本公开的实施例的二维仿组织胶状模型的示例非限制性输入图像。二维模型数据实验被用来演示仿射变形方法的结果。仿组织胶状模型被受控制的机器压缩。模型中存在一个柱状区域，其比其他部分硬6倍，用于模拟乳腺病变。此外，注入液体的管道位于模型中，用于模拟血管，它会在压缩时产生最大的应变。
对模型执行2％压缩和5％压缩，同时压缩前后的两张超声图像被记录。仿射变形方法的性能与耦合滤波方法和压扩方法的性能进行比较。
图15至图17示出了耦合滤波方法(图15)、仿射变形方法(图16)和压扩方法(图17)在2％压缩下得到的示例性相关系数。包含113×31像素的窗口被用于寻找其他图像中的匹配图案，相关系数作为匹配测度。一共385×257＝98945个窗口被使用，每个窗口最大的相关系数被计算并显示出来。在2％压缩的情况下，组织变形和散斑图案变化都不大。
所有这三种方法都得到了高于0.9的平均相关系数，说明大部分的散斑图案变化已被补偿。耦合滤波方法、仿射变形方法和压扩方法得到的平均相关系数分别为0.96、0.95和0.92。耦合滤波方法和仿射变形方法相比于压扩方法能够减小更多的散斑图案变化。仿射变形方法与耦合滤波方法性能相似。
在管道部分，所有这三种方法得到的相关系数都显著下降。这是因为管道部分的组织运动不同于其他胶状区域，而且三种方法中使用的局部区域都不够小，所以管道部分中的组织运动的估计会被周围不同种类的组织运动干扰。更小的相关性窗口可以提高管道附近的估计，但是会牺牲胶状区域的估计精度，因为在胶状组织中将没有足够的散斑图案来进行鲁棒跟踪。
图18至图20示出了耦合滤波方法(图18)、仿射变形方法(图19)和压扩方法(图20)在2％压缩下得到的轴向应变估计。包含113×31像素的窗口被用于寻找其他图像中的匹配图案，并且显示了针对每个窗口估计的轴向应变。一共385×257＝98945个窗口被计算。在2％压缩的情况下，由于散斑图案变化不大，所以组织运动的鲁棒估计难度不大。
对于所有这三种方法，都在硬质夹杂物区域中检测到低轴向应变，表明所有这三种方法都可以区分具有不同硬度的组织区域。耦合滤波方法和仿射变形方法仅检测到部分管道区域，而压扩方法得到管道区域远厚于真实情况。这主要是因为三种方法中定义的局部区域都不够小，所以管道区域中的组织运动的估计会被周围其他区域中的不同组织运动干扰。耦合滤波方法和仿射变形方法的胶状区域中的轴向应变估计没有像压扩方法得到的结果那样平滑，这是主要因为耦合滤波方法和仿射变形方法使用了组织不可压缩性约束，以为了更准确的横向应变估计而限制了估计精度。
图21至图23示出了耦合滤波方法(图21)、仿射变形方法(图22)和压扩方法(图23)在5％压缩情况下得到的平均相关系数。包含113×31像素的窗口被用于寻找其他图像中的匹配图案，相关系数作为匹配测度。一共385×257＝98945个窗口被使用，每个窗口最大的相关系数被计算并显示出来。在5％压缩的情况下，组织变形和散斑图案变化都很大，其通常被认为是困难的。
耦合滤波方法和仿射变形方法得到的平均相关系数分别为0.87和0.83，而压扩方法得到的平均相关系数为0.65。相关系数的这种显著的下降说明压扩方法在大的形变情况下不足以补偿图像变化。
对于耦合滤波方法和仿射变形方法，在管道部分中相关系数都显著下降。硬质夹杂物区域周围的相关系数也下降。这是因为不同区域中的组织运动是不同的，而且该两种方法中定义的局部区域都不够小，所以一个区域中的组织运动的估计会被其他区域中的不同组织运动所干扰。
参照图24至图26，示出了耦合滤波方法(图24)、仿射变形方法(图25)和压扩方法(图26)在5％压缩情况下得到的示例性轴向应变估计。包含113×31像素的窗口被用于寻找其他图像中的匹配图案，并且显示了针对每个窗口估计的轴向应变。一共385×257＝98945个窗口被使用。在5％压缩的情况下，散斑图案变化很大，所以被认为难度很大。
压扩方法不能产生轴向应变的可靠估计，这是因为它不能有效的补偿散斑图案变化。对于耦合滤波方法和仿射变形方法，都在硬质夹杂物区域中检测到低轴向应变，说明这两种方法都可以区分具有不同硬度的组织区域，即使在极其困难的情况。耦合滤波方法和仿射变形方法可以检测到仅部分的管道区域。在硬圆柱区域周围也观察到轴向应变估计的小偏差。这主要是因为这两种方法中使用的局部区域都不够小，所以管道区域中的组织运动的估计会被其他区域中的不同组织运动所干扰。耦合滤波方法和仿射变形方法的胶状区域中的轴向应变估计不够平滑，这主要是因为该两种方法使用了组织不可压缩性约束，限制了轴向应变的估计精度，以为了更准确的横向应变估计。
此外，三维仿真和二维模型数据实验均被用来演示对B超图像应用仿射变形方法的可行性。仿真和模型数据实验中采用了相同的设置，并且比较了使用B超图像和射频超声图像的散斑跟踪方法的性能。
在三维仿真中，以下结果被比较：使用射频图像的仿射变形方法、使用B超图像的仿射变形方法、使用射频图像的直接相关方法、使用B超图像的直接相关方法应用于。使用射频图像的仿射变形方法的性能、使用B超图像的仿射变形方法的性能、使用射频图像的直接相关方法的性能、使用B超图像的直接相关方法的性能在以下三种三维组织运动情况下进行了比较：分别是横向压缩和扩张、沿轴向-横向平面的剪切、和横向旋转。所有这三种情况均包括极大运动情况。
参考图27至图29，示出了根据本公开的实施例的在三维超声仿真实验中平均相关系数相对于运动的示例非限制性变化曲线。图27对应于横向压缩和扩张；图28对应于沿轴向-横向平面的剪切；以及图29对应于横向旋转的情况。元素2702显示仿射变形方法应用于射频图像，元素2704显示仿射变形方法应用于B超图像，元素2706显示无补偿的相关方法应用于射频图像，元素2708显示无补偿的相关方法应用于B超图像。
两种方法在所有这三种情况下使用B超图像得到的平均相关系数都远高于射频图像上得到的相关系数，这部分是因为B超图像中的相位信息被去除了。当组织经历复杂运动时，射频图像中的相位信息改变很多，这导致在匹配时的相关系数的急剧下降。总之，在B超图像中的散斑图案变化远小于射频图像中的散斑图案变化。
参考图30至图32，示出了在上述三种情况下运动估计的均方误差相对于运动的变化曲线。图30对应于横向压缩和扩张；图31对应于沿轴向-横向平面的剪切；以及图32对应于横向旋转。元素2702显示仿射变形方法应用于射频图像，元素2704显示仿射变形方法应用于B超图像，元素2706显示无补偿的相关方法应用于射频图像，元素2708显示无补偿的相关方法应用于B超图像。
在所有这三种情况下，使用射频图像的仿射变形方法得到的均方误差接近于0，除了一些极端情况(例如，针对轴向-横向剪切的情况所应用的切应变大于8％，针对横向旋转的情况所应用的旋转角度大于0.07弧度)。然而，使用B超图像的仿射变形方法得到的均方误差则逐渐增大(例如，在横向压缩/扩张的情况)，这主要是因为B超图像中的散斑图案不足以用于鲁棒应变估计。如果特征运动去相关没有被补偿，则组织运动估计会失败，而使用B超图像的直接相关方法的性能略好于使用射频图像的直接相关方法，这是因为在B超图像中散斑图案变化要小得多。
在二维模型数据实验中，比较了仿射变形方法应用于B超图像和射频图像的结果。在仿组织胶状模型压缩的情况下比较了B超图像仿射变形和射频图像仿射变形的性能。模型分别被压缩2％和5％。在2％压缩的情况下，组织变形不大，所以组织运动的鲁棒估计难度不太大，而在5％压缩的情况下，组织变形很大，所以估计难度很大。
参考图33，示出了在二维仿组织模型被2％压缩下仿射变形方法应用于B超图像得到的示例非限制性相关系数估计。相比于图16中的在二维仿组织模型被压缩2％下仿射变形方法应用于射频图像得到的相关系数估计，图33中仿射变形方法应用于B超图像得到的平均相关系数在每一点都接近于1(比如即使在管道区域都接近于0.995)，而针对使用RF图像的仿射变形方法的管道区域中相关系数急剧下降(注意，针对使用B超图像的仿射变形方法，也观察到管道区域中的相关系数略低)。这部分是因为散斑图案变化在B超图像中远小于在射频图像中的变化。
参考图34，示出了在二维仿组织模型被2％压缩下仿射变形方法应用于B超图像得到的示例非限制性轴向应变估计。图16示出了在二维仿组织模型被压缩2％下仿射变形方法用于射频图像得到的轴向应变估计。图34和图16提供了仿射变形方法用于B超图像和射频图像的轴向应变估计的定性验证。
从使用射频图像的仿射变形方法的最优运动参数中直接推出轴向应变，而对于使用B超图像的仿射变形方法，散斑图案不足以用于鲁棒应变估计，并且从所估计的位移中间接推出轴向应变。使用B超图像和射频图像的仿射变形方法都可以正确地描述硬质夹杂物区域和管道部分，呈现了使用B超图像和RF图像的仿射变形方法之间的相似的估计效果。尽管在B超图像实验中，散斑图案不足以(即，没有足够精细纹理)用于鲁棒应变估计，但组织位移能够得到稳定估计。通过对估计的组织位移求导得到的轴向应变估计示出了相当好的性能。作为估计位移的梯度的胶状区域中的轴向应变估计不够光滑，部分因为求导运算会扩大估计中的噪声。
参考图35，示出了在二维仿组织模型被5％压缩下仿射变形方法应用于B超图像得到的示例非限制性相关系数。相比于图22中的在二维仿组织模型被5％压缩下仿射变形方法应用于射频图像得到的相关系数估计，图35中仿射变形方法应用于B超图像得到的平均相关系数在每一点都接近于1(比如即使在管道区域都接近于0.99)，而对于使用射频图像的仿射变形方法，管道区域中的相关系数急剧下降(注意，对于使用B超图像的仿射变形方法，管道区域中的相关系数也略小)。
参考图36，示出了在二维仿组织模型被5％压缩下仿射变形方法应用于B超图像得到的示例非限制性轴向应变估计。图25示出了在二维仿组织模型被压缩5％下仿射变形方法用于射频图像得到的轴向应变估计。图36和图25提供了仿射变形方法用于B超图像和射频图像的轴向应变估计的定性验证。
从使用射频图像的仿射变形方法的最优运动参数中直接推出轴向应变，而对于使用B超图像的仿射变形方法，从所估计的位移中推出轴向应变。仿射变形方法应用于B超图像和射频图像都可以正确地描述硬质夹杂物区域和管道部分，呈现了使用B超图像和RF图像的仿射变形方法之间的相似的应变估计效果，即使组织形变很大。在这样困难情况下，鲁棒地估计出组织位移，并且通过对估计的组织位移求导得到的轴向应变估计示出了相当好的性能。作为所估计的位移的梯度的胶状区域中的轴向应变估计不够光滑，部分因为求导运算会扩大估计噪声。
物体跟踪
根据本公开中描述的一个或者多个实施例，与不需要用户交互或者训练数据的图像分割的物体跟踪相关地描述了多个非限制性方面。物体可以通过图像序列来跟踪是通过以下步骤实现的：把图像序列转换成输入矩阵、用一个秩比输入矩阵较小的低秩矩阵逼近输入矩阵、以及检测低秩矩阵的离群值。离群值对应于被跟踪的物体。
现在参照附图，首先参考图37，示出了根据本公开的实施例的示例非限制性系统3700，它便于在图像序列中对物体的特征进行分割和跟踪。图像序列可以包括二维图像或三维体积。图像序列也可以是图像/子图像序列。图像序列可以是超声波图像序列、计算机断层扫描图像序列、磁共振图像序列或任何其他类型的图像序列。在实施例中，图像序列对应于超声心动图序列，被跟踪的物体特征是二尖瓣瓣叶。
通过将图像序列转换成输入矩阵来对图像序列进行处理。系统3700包含转换组件3702，其便于把图像序列转换成输入矩阵。输入矩阵的列序号对应于图像序列的帧序号。输入矩阵的列是与图像序列相关的向量化图像。在实施例中，输入矩阵的每一列是与图像序列相关的向量化图像。
输入矩阵通过低秩近似被一个低秩矩阵逼近。系统3700还包含近似组件3704，其进行低秩近似。近似组件3704用低秩矩阵来逼近输入矩阵，低秩矩阵的秩小于输入矩阵的秩。
输入矩阵可以被近似成低秩矩阵和一组离群值。低秩矩阵对应于图像序列中的背景部分。图像的背景部分对应于除了被跟踪物体特征之外的图像部分。低秩矩阵的离群值对应于被跟踪的物体的特征。系统3700还包含检测组件3706，其检测低秩矩阵中的一个或多个离群值，离群值对应于被跟踪的物体的特征。
系统3700还包括存储器3708，其存储组件3702-3706和处理器3708，处理器3708便于组件3702-3706中的一个或多个的执行。
参考图38，示出了根据本公开的实施例的示例非限制性系统(对应于转换组件3702)，它便于把图像序列3802转换成输入矩阵3804。转换组件3702便于将图像序列3802转换成输入矩阵3804。
图像序列3802具有帧序号{I1、I2、…In}。若输入矩阵3804对应于图像序列的帧序号{I1、I2、…In}，则图像序列3802被转换成具有列序号的输入矩阵3804。输入矩阵3804的每一列是来自图像序列3802的向量化后的图像。通过将图像像素的强度值按预设的顺序堆积成列向量，来对图像序列3802的各图像进行向量化。在本实施例中，预设的顺序是逐列扫描顺序。
参考图39，示出了根据本公开的实施例的示例非限制性系统(对应于近似组件3704)，它便于把输入矩阵3804建模成低秩矩阵3902和一组离群值3904。近似组件3704执行低秩近似，将输入矩阵3804近似为低秩矩阵3902和一组离群值3904，低秩矩阵3902的秩小于输入矩阵。
相比于物体的其他部分，更多的基向量被用来确定被跟踪的物体的特征。因为更多的基向量被用于确定物体的特征，所以低秩矩阵3902的离群值3904可以被认为是对应于被跟踪物体特征的位置。
图40至图42示出了根据本公开的实施例的方法和/或流程图。为了简化说明，方法被描绘和描述为可以由包括处理器的系统执行的一系列动作。然而，根据本公开的动作可以按各种顺序和/或同时发生，并且可以具有与本公开中没有呈现和描述的其他动作。此外，并非需要所有示出的动作来实现根据所公开主题的方法。另外，本领域技术人员将理解并领会到各方法可以通过状态图或事件可替代地被表示为一系列相互关联的状态。此外，还应当理解的是，本说明书中描述的方法能够被存储在制造的产品上，以便于将所述方法传送和转移到硬件装置。
参考图40，示出了根据本公开的实施例的方法4000的示例非限制性处理流程图，它便于对一系列图像进行分析。图像序列可以包括二维图像或三维体积。图像序列也可以是图像/子图像序列。图像序列可以是超声波图像序列、计算机断层扫描图像序列、磁共振图像序列或任何其他类型的图像序列。在实施例中，图像序列对应于超声心动图序列，方法4000便于对二尖瓣瓣叶进行跟踪。
在步骤4002中，图像序列被包括处理器的系统转换成输入矩阵。输入矩阵的列是与图像序列相关的向量化后的图像。在步骤4004中，通过系统，输入矩阵被低秩矩阵逼近。在步骤4006中，低秩矩阵的一个或多个离群值被系统检测。一个或多个离群值对应于图像序列中被跟踪的物体。
参考图41，示出了根据本公开的实施例的方法4100的示例非限制性流程图，它便于检测对应于图像特征的离群值。在步骤4102中，图像序列被处理器转换成输入矩阵。在步骤4104中，输入矩阵被分解成低秩矩阵和由离群值引起的强度变化之和。在步骤4106中，基于由离群值引起的强度变化，来检测低秩矩阵的离群值。
根据实施例，方法4100被用于在超声(超声心动图)图像序列中跟踪和分割二尖瓣瓣叶。从超声序列中跟踪二尖瓣瓣叶是很困难的，因为缺少可靠的特征，而且二尖瓣运动很快而且不规律。在超声图像中，二尖瓣瓣叶的强度和纹理值跟诸如心肌层的其他部分相似。而且，由于包括探针定位、声散斑、信号丢失、声阴影等的因素，诸如边缘的特征质量表现出大变化。
物体跟踪往往依赖于光流计算，以找到图像帧之间特征点的对应关系。但是，超声图像序列中的运动分析非常困难并且不可靠，因为特征运动去相关问题。在快速运动、大形变和低帧速的条件下，运动分析更加困难，对二尖瓣进行成像就是这种情况，尤其是对二尖瓣的三维成像。
方法4100提供可以不需要用户交互和训练数据的二尖瓣瓣叶自动跟踪方法。在步骤4102中，超声心动图像序列被处理器转换成输入矩阵。超声心动图序列的帧被向量化，向量化之后的图像对应输入矩阵的列。在步骤4104中，输入矩阵被分解成低秩矩阵和由离群值引起的强度变化之和。低秩矩阵对应于图像序列中的背景部分，比如心肌组织。由于二尖瓣瓣叶的运动相对于心肌运动需要更多的基向量来描述，所以低秩矩阵的离群值对应于二尖瓣瓣叶的像素。在步骤4106中，基于由离群值引起的强度变化来检测低秩矩阵的离群值，从而二尖瓣瓣叶的位置可以从心肌组织中分割并通过图像序列被跟踪。
超声心动图序列由n个图像组成。在4102中，n个图像的超声心动图序列被转换成输入矩阵，用D∈Rm×n表示，这里m和n分别是每个图像中的像素个数以及序列中的图像数量。
如果图像序列中不存在运动，则超声心动图序列中的所有图像都是相同的。因此，当不存在运动时，忽略图像中的噪声，那么D的秩等于1(rank(D)＝1)。
如果序列中存在运动，由于心脏运动存在因而该情况更现实，那么在心脏周期内图像和图像序列会变化，D的秩会大于1(rank(D)>1)。但是，由于帧之间存在相关性，矩阵D的矩阵谱(矩阵特征值的集合)中的能量分布会集中于有限多个主成分上。
D的能量谱的累积分布函数(cdf)定义如下：
这里r＝min(m,n)，σ1,…,σr是D的奇异值降序排列。超声心动图序列的cdf会在k>2的时候超过95％。
对应于不同区域的cdf是不同的。比如，包括二尖瓣瓣叶的瓣膜区域的相应cdf和与心肌层或血池对应的cdf会不一样。包括二尖瓣瓣叶的瓣膜区域的cdf比心肌层或血池的cdf更低。因此瓣膜区域的谱分布更分散。这是因为在一个心脏周期内与平滑运动的心肌层相比，二尖瓣瓣叶运动得很快而且组织形变大。
因此，描述由二尖瓣瓣叶的复杂运动所引起的图像变化需要更多的主成分。换句话而言，固定常数K，用一个秩为K的矩阵去拟合瓣膜区域所得到的残差会远大于拟合其他区域的残差。相应地，在4104中，输入矩阵被分解成低秩矩阵和由离群值引起的强度变化之和。低秩矩阵对应于图像序列中的背景部分，比如心肌组织或血池。由于更多的主成分被用来描述二尖瓣瓣叶的复杂运动所引起的图像变化，与心肌运动相比，需要更多的基向量来描述二尖瓣瓣叶运动，所以低秩近似中的离群值对应于二尖瓣瓣叶的像素。
换言之，在步骤4104中，输入矩阵被分解成低秩矩阵、由离群值引起的强度变化、以及随机噪声之和：
D＝B+E+ε，
这里D是输入矩阵,B是低秩矩阵,E是由一个多个离群值引起的强度变化，ε是随机噪声。
在实施例中，随机噪声可以被忽略，所以输入矩阵被分解成低秩矩阵和由离群值引起的强度变化之和：
D＝B+E。
在步骤4106中，可以基于由离群值引起的强度变化来检测低秩矩阵(B)的(具有强度变化E的)离群值，并且二尖瓣瓣叶从心肌中被分割并通过图像序列被跟踪。
可以通过比较输入矩阵的元素值的强度与低秩矩阵的相应值，来在步骤4106中检测离群值。当输入矩阵的元素值的强度和低秩矩阵的相应值不一样时，确定该元素是离群值。
根据本公开的实施例，图像特征的分割和物体的跟踪是通过能量函数的最小化和/或最优化来实现的。根据实施例，能量函数包含四个函数项：
第1项：描述输入矩阵和低秩矩阵之间的差异(比如，利用低秩矩阵拟合输入矩阵的残差平方和)；
第2项：衡量低秩矩阵的特性(比如，低秩矩阵的秩)；
第3项：测量离群值的绝对值特性(比如，离群值的个数)；
第4项：假设相邻的像素具有同样的类别(比如，离群值还是非离群值)。
应当理解的是，该能量函数可以具有可替代的或更多的函数项。此外，能量函数可以省略一个或多个项。例如，根据实施例，能量函数包括第1项、第2项和第3项，而不包括第4项。
最小化或者优化能量函数便于把输入矩阵近似成低秩矩阵并且检测低秩矩阵的离群值。最小化或者优化的实现可以由系统进行迭代，直至近似结果和检测结果之间到达收敛后完成。图42描述了一种迭代方法。
参考图42，示出了根据本公开的实施例的便于通过迭代地最小化能量函数来分割和跟踪图像的特征的方法4200的示例非限制性的处理流程图。在步骤4204中，图像序列被转换成输入矩阵。能量函数通过步骤4204至步骤4208中的迭代过程被最小化或最优化。在迭代过程中，在步骤4204，估计低秩矩阵以近似输入矩阵。根据实施例，这个估计可以通过矩阵补全来完成。在步骤4206中，低秩矩阵的离群值可以被检测。在实施例中，这个检测可以通过图切割方法完成。在步骤4208中，估计和检测被迭代直到达到收敛。在步骤4210中，可以从收敛后的离群值的位置中输出分割结果。
图43至图45示出了根据上述系统和方法在超声心动图序列中跟踪二尖瓣瓣叶的示例。图43示出了用矩阵谱分析过程来跟踪二尖瓣瓣叶的可行性。图44和图45示出了使用经胸超声心动图(TTE)数据的二尖瓣瓣叶跟踪的示例结果。
图43是根据本公开的实施例的示出了矩阵谱分析过程的可行性的示例非限制性示图4300。图43中的(A)示出了手动分割的瓣膜区域4302、心肌和血池区域(参考区域)4304以及整个图像(4306)。图43中的(B)示出了以上三个区域4302至4306的谱能量累积分布函数对k的曲线。
如图43中的(B)所示，与不同区域4302至4306对应的累积分布函数(cdf)也不一样。参考区域4304的cdf和整个图像的4306的cdf高于瓣膜区域4302的cdf。由于瓣膜区域4302的累积分布函数曲线低于心肌和血池区域的累积分布函数曲线，因此瓣膜区域的谱谱分布比心肌和血池区域更均匀。
瓣膜区域的谱分布与心肌和血池的频谱分布相比更均匀，这是因为二尖瓣瓣叶运动得很快而且组织形变大，相对而言心肌层在一个心脏周期内的运动平缓。因此，描述由二尖瓣瓣叶的复杂运动所引起的图像变化需要更多的主成分。这说明将图像表示为低秩矩阵和由离群值引起的强度变化之和适合于分离二尖瓣瓣叶。低秩矩阵对应于图像序列中的背景部分，比如心肌组织或血池。由于更多的主成分被用来描述由二尖瓣瓣叶的复杂运动引起的图像变化，因此与心肌运动相比，更多的基向量被用来描述二尖瓣瓣叶运动。
除了在图像中占据相对较小且相连的区域的二尖瓣瓣叶之外，超声心动图序列可以很好地被近似成低秩矩阵。因此，二尖瓣瓣叶分割问题可以描述成在低秩表达中检测离群值的问题。
描述中使用了以下符号。ij表示序列的第j帧中的第i个像素。S是描述离群值的位置的二值矩阵，即，如果ij是二尖瓣瓣叶像素则Sij＝1；反之则Sij＝0。PS(X)表示把矩阵X正投影到由S支持的矩阵的线性空间：
并且PS⊥(X)是互补投影，PS(X)+PS⊥(X)＝X。
描述中使用了三种矩阵X的范数。表示l1-范数。是弗罗贝尼乌斯(Frobenius)范数。||X||*表示核范数，也就是奇异值的求和。
任务是给出寻找低秩矩阵B的离群值的位置的对S的最优估计。由于B和S是未知的，需要同时估计它们。以下能量函数被最小化以估计B和S：
s.t. rank(B)≤K，
K是一个预设的常数。第一项降低非瓣膜区域的拟合残差，第二项调整二尖瓣瓣叶像素的稀疏性。
由于瓣叶像素在图像序列中空间上和时间上都会是相邻的，对S引入平滑正则化。考虑到图G＝(V,E)，其中V是表示序列中所有m×n像素的顶点的集合，E是在空间上和时间上连接相邻像素的边缘的集合。基于一阶马尔可夫随机场(MRFS)，采用以下能量函数来对S施加连续性：
这里A是G的加权节点-边关联矩阵，并且权重ωij,kl被定义为如下：
这里，DG是高斯平滑图像序列，μ是空间相邻的标准偏差，ωτ是时间连接的权重。引入空间上变化的ωij,kl利用边缘信息以改善分割。
为了使能量函数最小化可以被解决，利用核范数时对B的秩运算是松弛的，核范数被证明是秩运算的有效凸替代。加上平滑正则化并且以对偶形式写为后，能量函数的最终形式表示为如下：
以上定义的目标函数是非凸的并且同时包括连续变量和离散变量。对B和S的联合优化极其困难。因此，采用交替算法将对B和S的能量最小化分成两个步骤。
在B步骤中，S被固定为之前的估计值以下目标函数中的能量关于B被最小化：
这是一个矩阵补全问题，其可以例如使用SOFT-IMPUTE算法的多项式时间来解决。
在S步骤中，B被固定为之前的估计值以下目标函数中的能量关于S被最小化：
以上组合最优化问题可以使用图分割方法的多项式时间准确地解决。
图44和图45示出了矩阵谱分析过程的结果。图44是根据本公开的实施例的通过对二维超声心动图序列的矩阵谱分析方法来跟踪二尖瓣瓣膜的示例非限制性结果。图45是根据本公开的实施例的通过对三维超声心动图序列的矩阵谱分析方法来跟踪二尖瓣瓣膜的示例非限制性结果。
所有的图像序列都是在一个心脏周期内用飞利浦iE33超声成像系统和用于二维成像的S5-1探头(图44)和用于三维成像的X3-1探头(图45)采集得到。在具备3.4GHz英特尔i7CPU和8GB RAM的PC上运行的MATLAB上实现了该系统和方法，需要大约30秒来分析二维序列(图44)以及需要24.4分钟来分析三维序列(图45)。例如，可以通过C或其他编程语言的实现或者通过便于并行处理的GPU计算或其他计算系统的执行来改善效率。
图44示出了在两个二维超声心动图图像序列(A)和(B)上的两个示例结果。每一个序列包含记录了一个心脏周期的大约30帧，每一图像的大小是225×300像素。为了简化说明和描述，对于每一个序列(A)和(B)，显示了瓣膜张开过程4402和4404中的三帧。
二尖瓣瓣叶被检测为低秩矩阵的离群值(在帧4406和4408中被标注出来)。低秩矩阵捕获了心肌层的全局运动(如帧4410和4412所示)。
自动分割方法得到的瓣膜轮廓和瓣膜的手绘跟踪结果进行了定量比较。自动方法和手绘跟踪得到的瓣膜轮廓之间的平均绝对距离(MAD)被计算。第一个序列(A)上的平均绝对距离是0.89±0.71毫米，第二个序列(B)上的平均绝对距离是0.94±0.18毫米。(A)和(B)的平均绝对距离和成像系统的分辨率(～5毫米)在同一数量级上。
图45示出了从三维超声心动图序列中二尖瓣瓣叶跟踪的示例结果4500。三维的超声心动图序列包含17个体数据，每一个体数据具有80×80×80个立体像素。所选体数据的分割结果被显示。
整个二尖瓣瓣叶上的三个正交的二维截面(4502-4506)以体数据的形式被显示，二维截面上叠加有分割结果。通过表面渲染得到的二尖瓣瓣叶的三维模型4508也被显示。
计算设备示例
上述便于图像分析的系统和方法可以以软件、硬件或其组合来实现。图46和图47提供了用于上述系统和方法的硬件环境。图46示出了可以与上述系统和方法关联地使用的计算环境4600。图47示出了可以与上述系统和方法关联地使用的计算网络4700。应当理解的是，也使用人工智能来实现本文中描述的系统和方法。
参考图46，示出了可以实现一个或多个实施例的适合计算系统环境4600。计算系统环境4600仅是适合计算环境的示例，并不旨在限制任何实施例的使用或功能的范围。计算环境4600也不应被解释为具有与示例性操作环境4600中示出的各组件中的任何一个或组合相关的依赖性或要求。
参考图46，计算系统环境4600可以包括计算机4610，其可以是手持或非手持计算机。计算机4610可以仅仅能够与成像装置(例如，记录或转换图像的装置)连接。然而，计算系统环境4600可以是任何其他计算装置，其具备执行本文中描述的方法的处理器、以及显示器，例如，台式计算机、笔记本电脑、移动电话、移动上网装置、平板电脑等。计算机4610的组件可以包括但不限于处理单元4620、系统存储器4630、和系统总线4621，系统总线4621将包括系统存储器的各种系统组件耦接至处理单元4620。例如，本文中描述的方法可以存储在系统存储器4630中，并且由处理单元4620执行。
计算机4610还可以包括各种计算机可读介质，其可以是可以由计算机4610访问的任何可用介质。系统存储器4630可以包括诸如只读存储器(ROM)和/或随机存取存储器(RAM)之类的易失性和/或非易失性形式的计算机存储介质。通过示例和非限制性地，存储器4630还可以包括操作系统、应用程序、其他程序模块和程序数据。
计算装置通常包括各种介质，其可以包括计算机可读存储介质和/或通信介质，本文中使用的该两个术语如下所述彼此不同。
计算机可读存储介质可以是可以由计算机访问的任何可用存储介质，并且包括易失性和非易失性介质、可以拆卸和非可拆卸介质。通过示例和非限制性地，计算机可读存储介质可以按与用于诸如计算机可读指令、程序模块、结构数据或非结构数据之类的信息的存储的方法或技术相关联地实现。计算机可读存储介质可以包括但不限于RAM、ROM、EEPROM、闪存或其他存储技术；CD-ROM、数字多功能盘(DVD)或其他光盘存储；盒式磁带、磁带、磁盘存储或其他磁存储装置；或可以用来存储期望信息的其他有形和/或非临时性介质。一个或多个本地或远程计算装置可以经由访问请求、询问或其他数据检索协议来访问计算机可读存储介质，用于对于介质所存储的信息的各种操作。
通信介质通常将计算机可读指令、数据结构、程序模块或其他结构化或非结构化数据表现为诸如调制数据信号(例如，载波或其他传输机制)之类的数据信号，并且包括任何信息传递或传输介质。术语“调制数据信号”或信号指代这样的信号：其一个或多个特性被设置或改变，从而对一个或多个信号中的信息进行编码。通过示例和非限制性地，通信介质包括有线介质(诸如有线网络或直接有线连接)、以及无线介质(诸如声学、RF、红外和其他无线介质)。
用户可以通过输入装置4640将命令和信息输入计算机4610，例如选择图像的特征。监视器或其他类型显示装置(例如，触摸屏或虚拟显示器)也可以经由诸如输出接口4650之类的接口连接至系统总线4621。
计算机4610可以在使用对一个或多个其他远程计算机(诸如远程计算机4670)的逻辑连接在网络或分布环境下操作。远程计算机4670可以是个人计算机、服务器、路由器、网络PC、对端装置或其他公共网络节点、或者其他远程介质消耗或传输装置，并且可以包括与计算机4610相关的上述任何或所有组件。图46中示出的逻辑连接包括网络4671，诸如局域网(LAN)或广域网(WAN)，但也可以包括其他网络/总线。这样的网络环境在家、办公室、企业范围计算机网络、内联网和因特网中是常见的。
参考图47，示出了示例性网络或分布式计算环境4700的示意图。图46的计算机4610可以在图47的网络中操作。分布式计算环境包括计算对象4710、4712等、以及计算对象或装置4720、4722、4724、4726、4728等(包括应用程序4730、4732、4734、4736、4738所表示的程序、方法、数据存储、可编程逻辑等)。应当理解的是，对象4710、4712等、以及计算对象或装置4720、4722、4724、4726、4728等可以包括诸如远程控制器、PDA、音频/视频装置、移动电话、MP3播放器、笔记本电脑等的不同装置。
每个对象4710、4712等、以及计算对象或装置4720、4722、4724、4726、4728等可以通过通信网络4740直接或间接地与一个或多个其他对象4710、4712、以及计算对象或装置4720、4722、4724、4726、4728等进行通信。即使图47中示出为单个要素，网络4740可以包括将服务提供给图47的系统的其他计算对象和计算装置，并且/或者可以代表多个互连网络(未示出)。每个对象4710、4712等、或者4720、4722、4724、4726、4728等也可以包含诸如应用程序4730、4732、4734、4736、4738之类的应用程序，其可以使用API或其他对象、软件、固件和/或硬件，适用于与根据各种实施例提供的力学性质测量相关的各种组件进行通信。
存在多个支持分布式计算环境的系统、组件和网络构造。例如，计算系统可以通过有线或无线系统、通过本地网络或广泛分布网络连接在一起。通常，多个网络耦合至互联网，互联网针对广泛的分布式计算提供基础设施并且包含多个不同网络，但是任何网络基础设施可以用于对各种实施例中描述的技术的示例通信。
作为另外的非限制性示例，本文中描述的各种实施例适用于任何手持式、便携式和其他计算装置，并且各种计算对象可以被考虑为与本文中描述的各种实施例相关地使用，即，在装置可以请求基于指向的服务的任何地方。因此，图47中描述的通用远程计算机仅是一个示例，并且本主题公开的实施例可以利用任何具有网络/总线互用性和互动的客户端来实现。
尽管不是必须的，任何实施例可以部分地通过操作系统来实现，从而被装置或对象的服务开发人员使用，和/或包括在与可操作组件相关地操作的应用软件中。软件可以在诸如程序模块之类的由一个或多个计算机(例如，客户端工作站、服务器或其它装置)执行的计算机可执行指令的一般上下文中进行描述。本领域技术人员将理解网络互动可以利用各种计算机系统构造和协议来实现。
上述描述包括本主题公开的实施例的示例。当然，不可能描述每一个可想到的用于描述权利要求主题的组件或方法的组合，但应当理解的是，各种实施例的很多其他组合和排序是可能的。因此，权利要求主题旨在包含落入所附权利要求书的精神和范围内的所有这样的改变、变型和变化。此外，在每个处理中一些或全部处理块出现的顺序不应理解为限制性的。相反，应当理解的是，一些处理块可以按本公开中没有示出的各种顺序而执行。此外，包括摘要中的描述的本公开的示出实施例的上述描述不旨在穷尽的或者将所公开的实施例限制于所公开的精确形式。当为了示出目的在本公开中描述特定实施例和示例时，各种变型是可能的，其被考虑在相关领域技术人员可以理解的所述实施例和示例的范围内。
尤其，对于上述组件、模块、系统等所执行的各种功能，另有说明除外，用来描述这样的组件的术语旨在对应于执行上述组件具体功能(即，功能相等)的任何组件，即使结构上不等于执行本文中权力要求主题的示例方面功能的所公开的结构。上述系统、装置和电路相对于一些组件和/或块之间的互动而描述。应当理解的是，这样的系统、装置、电路和组件和/或块可以包括那些组件或指定子组件、一些指定组件或子组件、和/或附加组件，并且根据上述各种排列和组合。子组件也可以实现为可通信地耦合至其他组件的组件，而不是包括在主组件中的组件(分层的)。另外，应当注意，一个或多个组件可以结合成提供集合功能的一个组件，或者被划分成一些分离的子组件，并且可以提供任何一个或多个诸如管理层的中间层，以可通信地耦合至所述子组件，以便提供整体功能。本公开中描述的任何组件也可以与本领域技术人员已知的本公开中没有具体描述的一个或多个其他组件进行交互。
另外，本文中使用的词语“示例”或“示例性”意味着作为示例、例子或说明。本文中作为“示例性”描述的任何方面和涉及不必须被理解为对于其他方面或涉及是优选的或有益的。相反，单词“示例性”的使用是旨在呈现具体方式的构思。如在本申请中的使用，术语“或”旨在表示包括的“或”而不是排除的“或”。换言之，除非另有具体说明，或者从上下文中清楚得知，“X采用A或B”旨在表示任何自然包括性排列。换言之，若X采用A；X采用B；或X采用A和B，则“X采用A或B”满足上述任何情况。另外，本说明书和所附权利要求书中使用的冠词“一个”应通常理解为表示“一个或多个”，除非另有具体说明，或者从上下文清楚指示为单数形式。
另外，当相对于一些实施例中的仅一个实施例已公开一个方面时，这样的特征可以与其他实施例的一个或多个其他特征相结合，例如可以是针对任何给定的或特定的应用而期望的和有益的。此外，就详细说明书或权利要求书中使用的术语“包括”、“包括...的”、“具有”、“包含”及其变型、以及其他相似词语来说，这些术语旨在包括性的，在一定程度上与术语“包括”作为不排除任何附加或其他要素的开放性过渡词类似。

Claims (10)

1.一种图像运动分析的方法，包括：
通过包括处理器的系统把图像序列转换成输入矩阵，其中所述输入矩阵的列对应于与所述图像序列相关的向量化后的图像；
通过所述系统利用低秩矩阵来逼近所述输入矩阵，所述低秩矩阵的秩小于所述输入矩阵的秩；
通过所述系统检测所述低秩矩阵的一个或多个离群值；以及
通过所述系统把所述输入矩阵分解成所述低秩矩阵、由所述一个或多个离群值引起的强度变化、和随机噪声之和，其根据以下等式描述：
D＝B+E+ε，
其中D是所述输入矩阵,B所述低秩矩阵,E是由所述一个或多个离群值引起的强度变化，ε是所述随机噪声。
2.根据权利要求1所述的方法，还包括：所述系统忽略所述随机噪声。
3.根据权利要求1所述的方法，其中所述检测的步骤还包括：比较所述输入矩阵的元素值的强度和所述低秩矩阵的对应值。
4.根据权利要求3所述的方法，其中所述检测的步骤还包括：根据所述输入矩阵的元素值的强度和所述低秩矩阵的对应值之间的不同来识别所述一个或者多个离群值。
5.根据权利要求1所述的方法，还包括：通过所述系统使能量函数最小化以便于实现逼近的步骤和检测的步骤。
6.根据权利要求5所述的方法，其中所述最小化的步骤还包括：迭代地进行逼近步骤和检测步骤直到在逼近与检测之间达到收敛。
7.根据权利要求5所述的方法，其中所述最小化的步骤还包括：使具有第一项、第二项和第三项的能量函数最小化，所述第一项描述所述输入矩阵和所述低秩矩阵之间的差别，所述第二项衡量所述低秩矩阵的特性，所述第三项衡量所述一个或多个离群值的绝对值特性。
8.根据权利要求7所述的方法，其中所述最小化的步骤还包括：使具有第四项的能量函数最小化，所述第四项假设彼此相邻的像素具有同一类别。
9.根据权利要求8所述的方法，其中所述同一类别是指离群值或者非离群值中的一个。
10.根据权利要求1所述的方法，其中所述转换步骤还包括：将超声图像序列、计算机断层扫描图像序列、或磁共振成像序列转换成所述输入矩阵。
CN201280028147.0A 2011-06-09 2012-06-07 基于图像的跟踪 Active CN103814384B (zh)

Applications Claiming Priority (5)

Application Number Priority Date Filing Date Title
US201161457813P 2011-06-09 2011-06-09
US61/457,813 2011-06-09
US201261631815P 2012-01-12 2012-01-12
US61/631,815 2012-01-12
PCT/CN2012/000782 WO2012167616A1 (en) 2011-06-09 2012-06-07 Image based tracking

Publications (2)

Publication Number Publication Date
CN103814384A CN103814384A (zh) 2014-05-21
CN103814384B true CN103814384B (zh) 2017-08-18

Family

ID=47295429

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201280028147.0A Active CN103814384B (zh) 2011-06-09 2012-06-07 基于图像的跟踪

Country Status (3)

Country Link
US (1) US9390514B2 (zh)
CN (1) CN103814384B (zh)
WO (1) WO2012167616A1 (zh)

Families Citing this family (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9286691B2 (en) * 2009-04-17 2016-03-15 The Hong Kong University Of Science And Technology Motion estimation and compensation of feature-motion decorrelation
US9280711B2 (en) 2010-09-21 2016-03-08 Mobileye Vision Technologies Ltd. Barrier and guardrail detection using a single camera
US9959595B2 (en) * 2010-09-21 2018-05-01 Mobileye Vision Technologies Ltd. Dense structure from motion
KR20150015680A (ko) * 2013-08-01 2015-02-11 씨제이씨지브이 주식회사 특징점의 생성을 이용한 이미지 보정 방법 및 장치
WO2015075718A1 (en) * 2013-11-19 2015-05-28 Ramot At Tel-Aviv University Ltd. Method and system for tracking a region in a video image
CN103886615B (zh) * 2013-12-31 2016-04-20 华中科技大学 一种x射线血管造影图像中多运动参数的分离估计方法
US10310074B1 (en) * 2014-03-27 2019-06-04 Hrl Laboratories, Llc System and method for denoising synthetic aperture radar (SAR) images via sparse and low-rank (SLR) decomposition and using SAR images to image a complex scene
CN104268878B (zh) * 2014-09-28 2018-03-09 华侨大学 一种超声图像中运动探针检测和定位方法
KR102276339B1 (ko) 2014-12-09 2021-07-12 삼성전자주식회사 Ｃｎｎ의 근사화를 위한 학습 장치 및 방법
CN107205725B (zh) 2015-01-29 2020-06-30 皇家飞利浦有限公司 通过实时超声应变成像对心肌梗塞的评估
WO2016209922A1 (en) * 2015-06-22 2016-12-29 Board Of Trustees Of Michigan State University Shear viscosity imaging with acoustic radiation force
US10366298B2 (en) * 2015-11-10 2019-07-30 Shoou Jiah Yiu Method and system for identifying objects in images
CN105957117B (zh) * 2016-04-26 2018-09-14 东南大学 并行磁共振的图像重建方法、装置及并行磁共振成像系统
CN107437083B (zh) * 2017-08-16 2020-09-22 广西荷福智能科技有限公司 一种自适应池化的视频行为识别方法
KR101859392B1 (ko) * 2017-10-27 2018-05-18 알피니언메디칼시스템 주식회사 초음파 영상 기기 및 이를 이용한 클러터 필터링 방법
CN108416359A (zh) * 2018-03-09 2018-08-17 湖南女子学院 一种乐谱识别系统及识别方法
CN109584273B (zh) * 2018-11-21 2022-11-18 西安电子科技大学 一种基于自适应收敛参数的运动目标检测方法
IT201900007815A1 (it) * 2019-06-03 2020-12-03 The Edge Company S R L Metodo per il rilevamento di oggetti in movimento
CN111449684B (zh) * 2020-04-09 2023-05-05 济南康硕生物技术有限公司 心脏超声标准扫查切面快速获取方法及系统
CN112233104B (zh) * 2020-10-27 2022-12-20 广州大学 实时位移场和应变场检测方法、系统、装置和存储介质
US12106492B2 (en) * 2021-11-18 2024-10-01 Volkswagen Aktiengesellschaft Computer vision system for object tracking and time-to-collision

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1179223A (zh) * 1995-03-22 1998-04-15 德国Idt国际数字技术有限公司 基于多帧的数据流分割的方法和装置
CN101964050A (zh) * 2010-09-06 2011-02-02 中国海洋大学 一种基于模型定阶和信号消噪的模态参数识别方法
CN102024252A (zh) * 2010-12-10 2011-04-20 清华大学 基于矩阵秩最小化的恢复水下失真图像的重构方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6574353B1 (en) * 2000-02-08 2003-06-03 University Of Washington Video object tracking using a hierarchy of deformable templates
US6807536B2 (en) * 2000-11-16 2004-10-19 Microsoft Corporation Methods and systems for computing singular value decompositions of matrices and low rank approximations of matrices
US9400921B2 (en) 2001-05-09 2016-07-26 Intel Corporation Method and system using a data-driven model for monocular face tracking
USH2222H1 (en) * 2005-10-13 2008-08-05 The United States Of America As Represented By The Secretary Of The Air Force Normalized matched filter—a low rank approach
US7961952B2 (en) * 2007-09-27 2011-06-14 Mitsubishi Electric Research Laboratories, Inc. Method and system for detecting and tracking objects in images
US7957900B2 (en) * 2008-02-08 2011-06-07 Gaurav Chowdhary Tracking vehicle locations in a parking lot for definitive display on a GUI
CN101567051B (zh) 2009-06-03 2012-08-22 复旦大学 一种基于特征点的图像配准方法
US8774558B2 (en) * 2010-11-29 2014-07-08 Microsoft Corporation Rectification of characters and text as transform invariant low-rank textures
US8959128B1 (en) * 2011-11-16 2015-02-17 Google Inc. General and nested Wiberg minimization

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1179223A (zh) * 1995-03-22 1998-04-15 德国Idt国际数字技术有限公司 基于多帧的数据流分割的方法和装置
CN101964050A (zh) * 2010-09-06 2011-02-02 中国海洋大学 一种基于模型定阶和信号消噪的模态参数识别方法
CN102024252A (zh) * 2010-12-10 2011-04-20 清华大学 基于矩阵秩最小化的恢复水下失真图像的重构方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"运动目标检测与跟踪算法研究";施雪梅;《中国优秀硕士学位论文全文数据库 信息科技辑》;20101115;27、34-36页 *

Also Published As

Publication number Publication date
CN103814384A (zh) 2014-05-21
US9390514B2 (en) 2016-07-12
US20140112544A1 (en) 2014-04-24
WO2012167616A1 (en) 2012-12-13

Similar Documents

Publication Publication Date Title
CN103814384B (zh) 基于图像的跟踪
JP6993334B2 (ja) 自動化された心臓ボリュームセグメンテーション
Tavakoli et al. A survey of shaped-based registration and segmentation techniques for cardiac images
CN101790748B (zh) 对3d数字医学图像中的解剖实体进行分割的方法
US9245091B2 (en) Physically-constrained modeling of a heart in medical imaging
US8265363B2 (en) Method and apparatus for automatically identifying image views in a 3D dataset
RU2669680C2 (ru) Инициализация модели на основе классификации видов
Duchateau et al. Infarct localization from myocardial deformation: prediction and uncertainty quantification by regression from a low-dimensional space
Thamsen et al. Synthetic database of aortic morphometry and hemodynamics: overcoming medical imaging data availability
Giachetti On-line analysis of echocardiographic image sequences
Akgul et al. A coarse-to-fine deformable contour optimization framework
CN110517249A (zh) 超声弹性图像的成像方法、装置、设备及介质
Almeida et al. Left-atrial segmentation from 3-D ultrasound using B-spline explicit active surfaces with scale uncoupling
Laumer et al. Weakly supervised inference of personalized heart meshes based on echocardiography videos
CN101116104A (zh) 用于分割多维数据集中的表面的方法、系统和计算机程序
Barba-J et al. Segmentation and optical flow estimation in cardiac CT sequences based on a spatiotemporal PDM with a correction scheme and the Hermite transform
Mukherjee et al. Computing myocardial motion in 4-dimensional echocardiography
Wang et al. Deep learning for automated echocardiogram analysis
Bhan et al. Patient-specific cardiac computational modeling based on left ventricle segmentation from magnetic resonance images
JP2009273882A (ja) 定量的局所測定結果および形態構造の複合４次元プレゼンテーションのための方法および装置
Timoner Compact representations for fast nonrigid registration of medical images
CN109863534A (zh) 用于分割解剖结构的二维图像的方法和装置
Ge et al. Echoquan-net: direct quantification of echo sequence for left ventricle multidimensional indices via global-local learning, geometric adjustment and multi-target relation learning
KR102715428B1 (ko) 장기 기능에 대한 정보 제공 방법 및 이를 이용한 장기 기능에 대한 정보 제공용 디바이스
Pazokifard et al. Automatic 3D modelling of human diaphragm from lung MDCT images

Legal Events

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