CN104156940B - 图像处理方法及装置 - Google Patents

图像处理方法及装置 Download PDF

Info

Publication number
CN104156940B
CN104156940B CN201410162145.1A CN201410162145A CN104156940B CN 104156940 B CN104156940 B CN 104156940B CN 201410162145 A CN201410162145 A CN 201410162145A CN 104156940 B CN104156940 B CN 104156940B
Authority
CN
China
Prior art keywords
image
value
gradient
vector
coordinate
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
CN201410162145.1A
Other languages
English (en)
Other versions
CN104156940A (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.)
GE Medical Systems Global Technology Co LLC
Original Assignee
GE Medical Systems Global Technology Co LLC
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 GE Medical Systems Global Technology Co LLC filed Critical GE Medical Systems Global Technology Co LLC
Publication of CN104156940A publication Critical patent/CN104156940A/zh
Application granted granted Critical
Publication of CN104156940B publication Critical patent/CN104156940B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G06T3/14—
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/75—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
    • G06V10/751—Comparing pixel values or logical combinations thereof, or feature values having positional relevance, e.g. template matching
    • G06V10/7515—Shifting the patterns to accommodate for positional errors
    • 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/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • 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/10072—Tomographic images
    • G06T2207/10088—Magnetic resonance imaging [MRI]
    • 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

Abstract

本发明为图像处理方法及装置，多次反复进行：第1工序，按第1及第2图像T、R上的每个坐标p决定像素值的梯度涉及的特征向量n，若坐标p的像素值的梯度向量∇T（p）、∇R’（p）大小为阈值THNorm以上，则将以该梯度向量的大小归一化该梯度向量后的值作为坐标p的特征向量n，若梯度向量小于该阈值则将0向量作为坐标p的特征向量n；第2工序，按第1及第2图像T、R’上的互相对应的每个坐标p算出将坐标p的特征向量彼此的内积的绝对值N次方而得的相关值，求出包含这些相关值的累计值的评价值D；和第3工序，变更第2图像R’以增大评价值D。从而包含同一对象的多个图像之间的对位中，有效抑制噪声对梯度的影响。

Description

图像处理方法及装置
技术领域
本发明涉及包含同一对象的多个图像之间的对位技术。
背景技术
一直以来，作为摄像装置，已知例如磁共振装置（MR）、放射线断层摄影装置（CT）、超声波装置（US）等。在这些摄像装置中，其各摄像模式（modality）分别存在优点/缺点，仅某一特定的摄像模式形成的图像存在诊断中的精度不足的情况。因此，近年来，不仅使用特定的摄像模式形成的图像，而使用多个不同摄像模式形成的图像进行诊断，从而谋求诊断精度的提高的尝试越来越多。
在利用多个不同摄像模式形成的图像的诊断中，按各摄像模式而图像的坐标系有所不同。因此，校正这些坐标系的差异、起因于内脏器官的变动/变形的错位的技术、即图像间的对位（Registration）技术尤为重要。
作为摄像模式彼此不同的多个图像间的对位手法，最普通的是采用相互信息量（Mutual Information）的手法（例如，参照非专利文献1等）。该手法广义上是基于图像的亮度值的手法（intensity based method）。即，利用相互信息量来进行对位时，对象图像间亮度值存在关联性成为前提。
然而，在US图像中，产生声影（Acoustic shadow），高反射体的后方的亮度值比本来的值还下降。另外血管的亮度值也依据血管的走向而发生变化。因此，例如在MR图像与US图像的对位、或者CT图像与US图像的对位中，经常发生亮度值的关联性不足的状况。对于这样的图像，适合采用基于图像的特征的手法（feature based method），而不是基于如相互信息量这样的亮度值的手法。
作为基于图像的特征的手法的代表例，提出采用标准梯度场（归一化梯度场，Normalized Gradient Field：NGF）的手法（例如，参照非专利文献2等）。标准梯度场，在算出图像上的坐标上各方向x、y、z的1次偏微分、即梯度向量（Gradient Vector）之后，以该梯度向量的长度（Vector Norm）对该梯度向量进行了归一化（normalized）。即，标准梯度场不依赖于像素值或者亮度值的大小或梯度的大小，而是仅表示梯度的方向的特征量。假设在某两个图像中互相对应的位置上产生相同方向的标准梯度场，则能将这两个图像的位置视为对齐。因此，该手法中，通过对标准梯度场所示的方向的对齐程度进行最优化，能够进行对位。
标准梯度场n在数学上表述如下。
（1）式中，示图像I的坐标p上的梯度向量。另外，表示该梯度向量的长度（Norm）。坐标p由x、y、z的各分量构成。
从（1）式能理解到，标准梯度场n单独不能分别起因于构造物的梯度和起因于噪声的梯度。因此，为了抑制起因于噪声（noise）的梯度的影响，非专利文献2中提出如下手法。
在此，ε是按照图像的噪声量而任意选择的常数。现有的方法中，按照（4）式来设定噪声项ε，考虑噪声项而对梯度向量进行归一化。
成为对位的指标的评价函数D以下式（任何一种都可）来定义。
在（5）、（6）式中，“T ，R ”分别表示成为对位的对象的两个图像。（5）式算出两个标准梯度场n（R，p）、n（T，p）的内积的平方的累计值。（6）式算出两个标准梯度场n（R，p）、n（T，p）的外积的平方的累计值。
现有技术文献
专利文献
非专利文献1：IEEE Trans. on Med. Imaging，16：187－198、1997；
非专利文献2：Proceeding of SPIE Vol. 7261、72610G－1、2009。
发明内容
发明要解决的课题
上述的采用标准梯度场的图像的对位手法中，起因于噪声的向量在标准梯度场上取微小的值。然而，如（5）、（6）式这样的评价函数取所有坐标上的内积或外积的累计值，因此即便一个个都是微小的值，但是如果其绝对数较多就不能忽略对评价函数的影响。因此，存在不能充分地排除起因于噪声的梯度（gradient）的影响的情况。
因为这样的情况，期待有包含同一拍摄对象的2个图像间的对位中，能够更加有效地抑制起因于噪声的梯度的影响的对位手法。
用于解决课题的方案
第1观点的发明，提供一种图像处理方法，多次反复进行以下工序而进行第1及第2图像的对位：
第1工序，按照包含同一拍摄对象的所述第1及第2图像中的每个坐标，决定该坐标上的像素值的梯度所涉及的特征向量，其中当该坐标的像素值的梯度向量的大小为既定阈值以上时，将以该梯度向量的大小对所述梯度向量进行归一化后的值决定为该坐标的特征向量，当所述梯度向量的大小小于所述既定阈值时，将0（零）向量决定为该坐标的特征向量；
第2工序，按照所述第1及第2图像中互相对应的每个坐标，算出与计算该坐标的特征向量彼此的内积的绝对值的N（N为自然数）次方而得的值相当的相关值，求出包含按每个所述坐标算出的相关值的累计值的评价值；以及
第3工序，以使所述评价值变得更大的方式变更所述第2图像。
第2观点的发明提供一种图像处理装置，其中包括：
决定单元，按照包含同一拍摄对象的第1及第2图像中的每个坐标，当该坐标上的像素值的梯度向量的大小为既定阈值以上时，将以该梯度向量的大小对所述梯度向量进行归一化后的值决定为该坐标的特征向量，当所述梯度向量的大小小于所述既定阈值时，将0（零）向量决定为该坐标的特征向量；
算出单元，按照所述第1及第2图像中互相对应的每个坐标，算出与计算该坐标的特征向量彼此的内积的绝对值的N（N为自然数）次方而得的值相当的相关值，求出包含按每个所述坐标算出的相关值的累计值的评价值；以及
变更单元，以使所述评价值变得更大的方式变更所述第2图像，
多次反复进行所述决定单元、算出单元及变更单元所进行的处理，进行所述第1及第2图像的对位。
第3观点的发明提供所述梯度向量为图像中各坐标轴向的1次偏微分的、上述第2观点的图像处理装置。
第4观点的发明提供所述相关值为计算所述特征向量彼此的内积平方而得的值的、上述第2观点或第3观点的图像处理装置。
第5观点的发明提供所述评价值为将所述相关值的累计值以所述第1或第2图像的大小进行归一化而得的值的、上述第2观点至第4观点的任一个观点的图像处理装置。
第6观点的发明提供所述变更包含平行移动、旋转、及变形中的至少一种的、上述第2观点至第5观点的任一个观点的图像处理装置。
第7观点的发明提供所述决定单元基于对所述第1及第2图像分别进行构造物强调滤波处理和/或噪声降低处理而得到的图像，求出所述梯度向量的、上述第2观点至第6观点的任一个观点的图像处理装置。
第8观点的发明提供所述第1及第2图像为医用图像的、上述第2观点至第7观点的任一个观点的图像处理装置。
第9观点的发明提供所述第1及第2图像其摄像模式互相不同的、上述第8观点的图像处理装置。
第10观点的发明提供所述第1及第2图像的一个为超声波图像（US图像）的、上述第9观点的图像处理装置。
第11观点的发明提供所述决定单元反复进行算出单元及变更单元所进行的处理，直至进行该处理的次数达到既定数，或直至所述评价值实质上收敛的、上述第2观点至第10观点的任一个观点的图像处理装置。
第12观点的发明提供还包括推定所述第1及第2图像的至少一个的对象区域中的噪声量，并基于该噪声量决定所述阈值的单元的、上述第2观点至第11观点的任一个观点的图像处理装置。
第13观点的发明提供所述对象区域为所述第1图像和所述第2图像之间共同的解剖学特征点的周边区域的、上述第12观点的图像处理装置。
第14观点的发明提供所述对象区域为所述至少一个的图像的中央区域的、上述第12观点的图像处理装置。
第15观点的发明提供所述对象区域为病变部的周边区域的、上述第12观点的图像处理装置。
第16观点的发明提供所述对象区域为全部图像区域的、上述第12观点的图像处理装置。
第17观点的发明提供所述决定阈值的单元求出所述对象区域中各像素上的梯度强度，并基于该梯度强度的直方图来决定所述阈值的、上述第12至第16观点的任一个观点的图像处理装置。
第18观点的发明提供所述决定阈值的单元根据以下的（a）～（f）式的任一个，或多个的组合的平均值算出噪声量Nqt的、上述第12观点至第16观点的任一个观点的图像处理装置。
Nqt＝Mfmax＋σ×HWHML （a）
Nqt＝Mfmax＋σ×HWHMR （b）
Nqt＝Mfmax＋σ×(HWHML＋HWHMR)/2 （c）
Nqt＝σ×Mfmax （d）
Nqt＝σ×Mmom1 （e）
Nqt＝σ×Mmom2 （f）
其中，Mfmax是在所述直方图分布上出现的峰值之中、梯度强度最低的峰值给出众数值的梯度强度；
HWHML是从所述直方图分布上的Mfmax观看低值侧的半值半宽度；
HWHMR是从所述直方图分布上的Mfmax观看高值侧的半值半宽度；
Mmom1是与所述直方图分布上的梯度强度从0到HWHMR的范围的重心相当的梯度强度；
Mmom2是与所述直方图分布上的梯度强度从HWHML到HWHMR的范围的重心相当的梯度强度。
第19观点的发明提供一种程序（program），用于使计算机（computer）作为上述第2观点至第18观点的任一个观点的图像处理装置起作用。
第20观点的发明提供一种程序，用于使计算机作为上述第12观点的图像处理装置起作用。
发明效果
依据上述观点的发明，当梯度向量的大小小于既定阈值时，特征向量成为0向量。因此，特征向量会成为仅由主要的构造物形成的分量构成。图像之间的对位的评价值包含计算每个坐标的特征向量彼此的内积的绝对值的N次方而得的相关值的累计值，因此能有效地抑制起因于噪声的梯度向量的影响。
附图说明
图1是概略地示出本实施方式所涉及的图像处理装置的结构的图；
图2是对位运算的概念图；
图3是示出按现有的方法进行的标准梯度场的例子的图；
图4是示出按本提案的方法进行的带有阈值判定的标准梯度场的例子的图；
图5是本实施方式所涉及的图像处理装置的图像对位处理的流程图；
图6是示出按现有的方法和本提案的方法进行的图像对位的例子的图；
图7是基于目标图像T及对象图像R’的至少一个来推定图像的噪声量的处理的流程图；
图8是示出梯度强度图像中梯度强度的直方图的例子的图；
图9是示出梯度强度的直方图中解析的特征量的第1图；
图10是示出梯度强度的直方图中解析的特征量的第2图；
图11示出按现有的方法进行的标准梯度场的例子和按本提案的方法（由噪声量的自动推定结果决定阈值THNorm）进行的带有阈值判定的标准梯度场的例子；
图12示出按现有的方法和本提案的方法（由噪声量的自动推定结果决定阈值THNorm）进行的图像对位的其他例子。
具体实施方式
以下，对发明的实施方式进行说明。此外，发明并不由这些限定。
图1概略地示出本实施方式所涉及的图像处理装置1的结构。如该图所示，图像处理装置1包括图像取得部2、对位条件设定部3、对位运算部4、和图像输出部5。
图像取得部2取得摄像模式彼此不同的两个图像。通常，根据用户的操作输入并取得这两个图像被。图像取得部2将这两个图像之中的一个，在对位处理中设定为固定的目标图像T，将另一个，在对位处理中设定为变更的对象图像R。此外，这里以R表示取得的原来的对象图像，以R’表示对位处理中的对象图像，以R”表示对位处理后的对象图像。
对位条件设定部3根据用户（user）的操作，设定作为对位运算的条件的最优化参数（parameter）。关于最优化参数的设定，随后进行详细说明。
对位运算部4根据所设定的最优化参数进行对位运算。对位运算通过将对对象图像R’加入移动或变形等的微小变更的处理、和算出目标图像T与对象图像R’的对位的评价值的处理，沿其评价值改善的方向反复进行来完成。
在本实施方式中，如图1所示，对位运算部4包括评价部（决定单元、算出单元）41、最优化部（变更单元）42、变换部（变更单元）43、和插值部（变更单元）44。此外，这是将对位运算部图式化时经常采用的结构图，评价部41也被称为“Cost Function（价值函数）”，最优化部42也被称为“Optimizer”，变换部43也被称为“Transformer”，插值部44也被称为“Image Interpolator（图形插值器）”。
对位运算部4通过反复进行这些评价部41～插值部44的处理来进行对位运算。此外，评价部41是发明中的决定单元及算出单元的一个例子。另外，最优化部42、变换部43及插值部44是发明中的变更单元的一个例子。以下，一边参照图像一边详细说明这些评价部41～插值部44的各部分的处理。
图2示出对位运算的概念图。在图2中，作为成为对位的对象的目标图像T及对象图像R的例子，提示腹部（肝脏）的US图像及MR图像。
对于目标图像T及对象图像R’的每一个，评价部41按每个坐标p，算出x、y、z各方向的1次偏微分、即梯度向量。此时，在算出梯度向量之前，优选对于目标图像T及对象图像R’的每一个，进行强调血管等的构造物的滤波处理和降低噪声的滤波处理的至少一种处理。此外，作为这些滤波处理，除了一般的现有的构造物强调滤波处理、噪声降低滤波处理之外，特别是能够对US图像采用US图像特有的滤波处理。作为US图像特有的滤波处理，例如可以考虑非专利文献3即Real－Time Speckle Reduction and Coherence Enhancement inUltrasound Imaging via Nonlinear Anisotropic Diffusion，Khaled Z. Abd－Elmoniem，Student Member，IEEE，Abou－Bakr M. Youssef，and Yasser M. Kadah*，Member，IEEE（IEEE TRANSACTIONS ON BIO MEDICAL ENGINEERING，VOL. 49，NO. 9，SEPTEMBER 2002）中公开的滤波处理等。
接着，对于目标图像T及对象图像R’的每一个，按每个坐标p，基于先求出的梯度向量∇T（p）、∇R’（p），按照下式算出带有阈值判定的标准梯度场（特征向量）n（T，p）、n（R’，p）。
在此THNorm是针对梯度向量的长度（Norm）的阈值。对于全部的梯度向量算出向量的Norm，如果Norm小于阈值就视为起因于噪声的梯度向量，将带有阈值判定的标准梯度场设为0向量。由此带有阈值判定的标准梯度场中，仅提取起因于构造物的向量，换言之对位中提供有效的信息的向量。此外，阈值THNorm例如根据实验或经验来决定。另外，阈值THNorm既可为目标图像T与对象图像R’之间共同的值，也可为分别独立的值。
图3示出按现有的方法形成的标准梯度场的例子。另一方面，图4示出按本提案的方法（经验决定阈值THNorm）形成的带有阈值判定的标准梯度场的例子。本提案的方法中，已知除去了起因于噪声的梯度场，只有起因于能成为对位的目标的内脏器官内的主要构造物的梯度场起支配作用。
如果算出了特征向量，则按目标图像T及对象图像R’中互相对应的每个坐标p，算出表示梯度场的对齐程度的评价函数D（T，R’）的值（评价值）。评价函数D如下式定义。
即，评价函数D（T，R’）是将目标图像T及对象图像R’中的各坐标p上的“特征向量彼此的内积的平方”（相关值）的累计值、即加法值或者总和，以图像的大小、即图像的总像素数进行了归一化的函数。
假设某两个向量平行则内积值取最大值，假设垂直则内积值取最小值。因此，如果以使（9）式的评价函数D取最大值的方式将两个图像的位置关系最优化，则该两个图像的梯度场会接近平行，达成对位。
此外，评价部41进行的评价函数的值的算出是每次对象图像R’变更、更新时进行的。
最优化部42根据所算出的评价函数的值，重复进行对象图像R’的变换参数的决定和变更，将变换参数最优化。该变换参数是规定对象图像R’的平行移动、旋转、变形等的变更方向和变更量的参数。变换参数的最优化的规则例如如下。
首次是将变换参数移位（shift）到事先设定的方向。第2次以后，在评价函数D的值比前一次得到改善的情况下，将变换参数移位到与前一次相同的方向，在变差的情况下，将变换参数移位到相反方向。将变换参数移位时的步（step）宽从事先设定的“最大步宽”开始。而且，在将变换参数向相反方向移位时，将对该时点的步宽乘以事先设定的“驰豫系数”而得的步宽设为新的步宽。
此外，作为最优化部42进行的变换参数的最优化的方法，可以采用例如梯度下降（Gradient Descent）法、共轭梯度（Conjugate gradient）法、LBFGS法、高斯－牛顿（Gauss-Newton）法、列文伯格－麦夸特（Levenberg-Marquardt）法、Amoeba法等。
变换部43根据在最优化部42决定的变换参数，对对象图像R’进行变换处理。
此外，作为变换部43进行的对对象图像R’的变换处理的方法，可以采用例如EulerTransform（欧拉变换）、Affine Transform（仿射变换）、B-Spline Transform（B样条插值变换）、Versor Transform（规范化四元数变换）、Quaternion Transform（四元数变换）等。
插值部44对在变换部43中已进行变换处理的图像进行插值处理，算出各坐标上的像素的像素值。在对对象图像R’进行变换处理之后，除了1个像素宽度单位的平行移动等特殊的情况外，大多数情况下像素不会重叠到坐标的位置。因此，需要通过这样的插值处理来求出各坐标的像素值。若完成该插值处理，则对象图像R’会被更新。
此外，作为插值部44进行的插值处理的方法，可以采用例如Nearest NeighborInterpolation（最邻近插值）、Linear Interpolation（线性插值）、B-SplineInterpolation（B样条插值）等。
对位运算部4多次反复进行评价部41～插值部44的处理。而且，在评价函数的前一次的值与本次的值的差量成为事先设定的“容许误差（tolerance）”以下的情况下、步宽成为小于事先设定的“最小步宽”的情况下、或将变换参数移位的反复运算的次数达到事先设定的“反复运算上限次数”的情况下，结束该反复运算。
作为对位运算条件，对位条件设定部3设定最优化部42进行的将变换参数移位的反复运算中的最优化参数。作为最优化参数，包含上述“最大步宽”、“最小步宽”、“反复运算上限次数”、“容许误差”及“驰豫系数”等。
图像输出部5在对位运算结束后，输出目标图像T和对位的对象图像R”。
以下，对本实施方式所涉及的图像处理装置进行的图像对位处理的流程进行说明。图5是本实施方式所涉及的图像处理装置进行的图像对位处理的流程图。
步骤S1中，图像取得部2取得对位的两个图像。而且，将取得的这两个图像的一个设定为目标图像T，将另一个设定为对象图像R。
步骤S2中，对位条件设定部3设定最优化参数。作为最优化参数，包含最大步宽、最小步宽、反复运算的上限次数、对位的容许误差、驰豫系数等。
步骤S3中，最优化部42决定对象图像R’的变换参数（平行移动、旋转移动、变形等）。在评价函数D的值比前一次得到改善的情况下，将变换参数移位到与前一次相同的方向。相反变差的情况下，将变换参数移位到相反方向。首次是决定适当的变换参数。
步骤S4中，变换部43根据所决定的变换参数，对对象图像R’进行变换处理。
步骤S5中，插值部44对对象图像R’进行插值处理。
步骤S6中，评价部41算出评价函数D（T，R’）的值、即评价值。
步骤S7中，对位运算部4进行反复运算的结束判定。在反复次数达到既定数的情况下，或者评价函数的值实质上收敛的情况下（达到容许误差或达到小于变换参数的移位量最小步宽的情况下），结束反复运算，进入步骤S8。除此之外的情况下，返回步骤S3。
步骤S8中，图像输出部5将最终的变换参数适用于对象图像，输出对位后图像。
图6示出按现有的方法和本提案的方法（经验决定阈值THNorm）进行的图像对位的例子。在图6中，从左到右依次为轴向（Axial）、径向（Sagittal）、冠状方向（Coronal）的各截面，从上到下依次为目标图像（该例子为US图像）、对象图像（该例子为MR图像）、按现有的方法形成的对位图像、按本提案的方法形成的对位图像。图中的箭头示出现有的方法和本提案的方法中能显著地视觉辨认到对位精度的提高的部位。如图6所示，本提案的方法比现有的方法还提高对位的精度。另外，能容易判别对位的计算结果是否达到期望的精度。
（变形例）
上述阈值THNorm也可以根据基于目标图像T及对象图像R’的至少一个的噪声量的推定结果来决定。
图7是根据目标图像T及对象图像R’的至少一个推定图像的噪声量的处理的流程图。
步骤T1中，在目标图像T及对象图像R’的至少一个中，设定推定噪声量的对象区域。图像中不仅包含诊断对象的内脏器官，还包含其他的内脏器官，因此优选以排除这些的方式设定对象区域。另外，在内脏器官或者医用图像中，噪声量并非恒定而有随着图像位置而不同的情况，因此通过将推定噪声量的对象区域重点限定在想要诊断的范围（例如病变部的周边等），能够更加有效地推定噪声量。作为上述对象区域而切出的范围，可以设为例如位于图像中央的第1部分范围、位于病变部周边的第2部分范围、或者目标图像T与对象图像R’之间位于共同的解剖学特征点的周边的第3部分范围等。第1部分范围例如可以设为对应于FOV的范围之中面积相当于FOV整体的1/2～1/3的图像中央的范围等。第2部分范围例如可以设为肿瘤等的病变部周边体积成为50～100mm3的范围等。另外，第3部分范围例如可以设为骨头形状、血管形状或血管分支等共同的解剖学特征点的周边体积成为50～100mm3的范围等。对象区域的设定，既可以手动设定也可以自动设定。在自动设定的情况下，根据需要，执行检测病变部或解剖学特征点的处理。此外，即便不执行本步骤，也可以执行以后的处理，但是要推定正确的噪声量，优选执行本步骤。在不执行本步骤的情况下，上述对象区域设为图像整体的区域。
步骤T2中，对于在步骤T1中设定的对象区域进行平滑（smoothing）处理。平滑处理为高斯滤波器（Gaussian Filter）等已知的滤波器。
步骤T3中，对于在步骤T2中生成的平滑处理完成图像，算出各图像位置上的梯度强度（Gradient Magnitude），生成将对应于各图像位置的像素的像素值设为该图像位置上的梯度强度的梯度强度图像。如众所周知地，梯度强度能通过算出x、y、z各方向的1次偏微分、即梯度向量（Gradient Vector）的大小而获得。
步骤T4中，对于在步骤3中生成的梯度强度图像，算出各图像位置上的梯度强度的直方图，解析直方图的特征量。图8示出梯度强度的直方图的例子。图8的直方图中，横轴表示梯度强度，纵轴表示频度的相对次数。在起因于噪声的梯度的情况下，该梯度强度为相对低值，其次数在直方图上集中在低值。另一方面，在起因于构造物的梯度的情况下，该梯度强度取相对中值/高值，其次数大范围分布。因而，图8的直方图中，在梯度强度较低的区域出现的波峰，可以认为对应于起因于噪声的梯度。
在此解析的特征量如下（参照图9、图10）：
・在直方图分布上出现的几个峰值之中、梯度强度最低的峰值给出众数值的梯度强度（图9的Mfmax）；
・直方图分布上的从Mfmax观看低值侧的半值半宽度（图9的HWHML）；
・直方图分布上的从Mfmax观看高值侧的半值半宽度（图9的HWHMR）；
・直方图分布上的与梯度强度从0到HWHMR的范围的重心相当的梯度强度（图10的Mmom1）；
・直方图分布上的与梯度强度从HWHML到HWHMR的范围的重心相当的梯度强度（图10的Mmom2）。
步骤T5中，算出噪声量。噪声量Nqt可由以下（a）～（f）式的任一个或多个组合的平均值算出。
Nqt＝Mfmax＋σ×HWHML （a）
Nqt＝Mfmax＋σ×HWHMR （b）
Nqt＝Mfmax＋σ×(HWHML＋HWHMR)/2 （c）
Nqt＝σ×Mfmax （d）
Nqt＝σ×Mmom1 （e）
Nqt＝σ×Mmom2 （f）。
在此，σ为任意常数，优选为2.0～3.0左右的值。在进行后述的图像对位（ImageRegistration）的过程中，用户打算重点仅以较大的构造物进行对位的情况下，最好将σ设定为大一点的值，用户打算也包含较小的构造物而进行对位的情况下，优选将σ设定为小一点的值。
阈值THNorm根据该推定的噪声量Nqt决定，但在本例中，将噪声量Nqt按照原样用作为阈值THNorm。
此外，在目标图像T及对象图像R’共同使用阈值THNorm的情况下，通过（1）从目标图像T中的噪声量的推定结果决定阈值THNorm；（2）从目标图像R’中的噪声量的推定结果决定阈值THNorm；（3）综合使用目标图像T中的噪声量的推定结果和目标图像R’中的噪声量的推定结果决定阈值THNorm等的方法来决定阈值。另外，在以目标图像T和对象图像R’分别决定阈值THNorm的情况下，从目标图像T中的噪声量的推定结果决定目标图像T用的阈值THNorm，从对象图像R’中的噪声量的推定结果决定对象图像R’用的阈值THNorm。
图11示出按现有的方法形成的标准梯度场的例子和按本提案的方法（由噪声量的自动推定结果决定阈值THNorm）形成的带有阈值判定的标准梯度场的例子。该图中，示出算出左端的图像（腹部US图像）中的标准梯度场的例子。可由图11理解，在现有的方法中不仅包含源至构造物的梯度场而且包含源至噪声的梯度场，但是在本提案的方法中仅明确给出源至构造物的梯度场（参照图11的白线）。
图12示出按现有的方法和本提案的方法（由噪声量的自动推定结果决定阈值THNorm）进行的图像对位的其他例子。图12中的图像分别显示目标图像（该例子中腹部US图像）、对象图像（该例子中腹部MR图像）、按现有的方法形成的对位图像、目标图像与按现有的方法形成的对位图像的融合（Fusion）显示、按本提案的方法形成的对位图像、目标图像与按本提案的方法形成的对位图像的融合（Fusion）显示。US图像、MR图像都是3D扫描后的图像，但在图12中，仅显示代表性的1个截面。在现有的方法的情况下以白线表示的血管稍微难以视觉辨认，但是在本提案的方法中能清楚地视觉辨认。另外，如图12中的白虚线所示，在现有的方法的情况下，能观察到目标图像与对位图像之间尚未校正完的位置偏移（misregistration），而在本提案的方法中，能几乎无位置偏移地达到对位。
这样，依据本实施方式所涉及的图像处理装置，如以下所说明的那样，可以进行高精度的对位。
在现有的方法中，起因于噪声的梯度向量在标准梯度场上取微小的值。然而，如（5）、（6）式那样，评价函数取所有坐标上的内积或外积的累计值。因此，即便内积或外积的一个个都为微小的值，在其绝对数较多的情况下不能忽略给出评价函数的影响。另一方面，本提案的方法中，当梯度向量的大小小于阈值THNorm时，带有阈值判定的标准梯度场成为0向量。因此，带有阈值判定的标准梯度场仅由主要的构造物形成的分量构成，能有效地抑制起因于噪声的梯度向量的影响。因而，本提案的方法能期待对位精度的提高。
另外，现有的方法中，如（5）、（6）式那样，算出两个标准梯度场的内积或外积的平方的累计值，评价函数的值取决于图像的大小（总像素数）。因此，难以从评价函数的绝对值预测对位的精度。本提案的方法中，由于以图像的总像素数归一化评价函数，所以评价函数的绝对值和对位的精度具有相关性。假设对位的两个图像完全相同，则（8）式的评价函数D的值成为1.0。另外，假设对位的两个图像中，50％的像素相同，剩余50％的像素完全没有相似性，则（8）式的评价函数D的值成为0.5。因此，如果预先指定用户所希望的对位精度，就能简单地判别对位的计算结果是否满足期望的精度。
在用多模式的图像进行的诊断中两个图像不一定是同时期摄影的，还要考虑存在有意的时差的可能性。例如，可以设想到一个图像为手术前摄影，另一个图像为手术后摄影的状况。在该情况下，一个图像上存在的病变不一定在另一个图像中存在。因而仅在一个图像上存在起因于病变的梯度向量的情况下，优选将这样的梯度向量从评价函数的计算中排除。然而，现有的方法中，仅在一个图像上存在起因于病变的梯度向量的情况下，另一个图像上与该病变对应的坐标上，存在起因于噪声的梯度向量。因此，这些标准梯度场的内积值会具有一些值，会影响评价函数。相对于此本提案的方法中，另一个图像上的对应于该病变的坐标上，存在0向量，因此内积值成为0，对评价函数不产生影响。因而，本提案的方法能期待对位精度的提高。
另外，作为进行对位的两个图像的组合，除了MR图像和US图像的组合外，还能适用CT图像和US图像、或MR图像和CT图像的组合等，所有的摄像模式的图像。但是，按采用标准梯度场的现有的方法进行的对位手法，尤其在对位对象的图像中包含US图像的情况下为有效的手法，因此按本提案的方法进行的对位手法也对包含US图像的对位特别有效。
另外，如果如上述变形例那样想从噪声量的自动推定结果决定阈值THNorm时，能够按每个图像最优化阈值THNorm，能够期待更高精度的对位。
此外，在本实施方式中，评价函数包含特征向量彼此的内积的平方，但是也可以包含特征向量彼此的内积的绝对值的N次方（N为自然数）。
另外，评价函数优选以图像的大小对特征向量彼此的内积的平方或内积的绝对值的N次方进行归一化，但是也可以不是归一化的函数。在对位的图像的大小、即图像的总像素数几乎不变的情况下，这样也能保持一定程度的对位精度。
另外，本实施方式是将发明适用于摄像模式互相不同的图像彼此的对位的例子，但是摄像模式相同也能适用于摄像的时相互相不同的图像彼此的对位。作为这样的图像，例如可以考虑手术前后的图像或造影摄影中的早期相和后期相的图像等。另外，发明不仅能适用于医用图像，还能适用于工业用的图像。
另外，本实施方式为图像处理装置，但是令计算机起到这样的图像处理装置功能的程序也是发明的实施方式的另一个例子。
附图标记说明
1 图像处理装置；2 图像取得部；3 对位控制条件设定部；4 对位运算部；41 评价部；42 最优化部；43 变换部；44 插值部；5 图像输出部。

Claims (18)

1.一种图像处理方法，多次反复进行以下工序以进行第1及第2图像的对位：
第1工序，按照包含同一拍摄对象的所述第1及第2图像中的每个坐标，决定该坐标上的像素值的梯度所涉及的特征向量，其中当该坐标上的像素值的梯度向量的大小为既定阈值以上时，将以该梯度向量的大小对所述梯度向量进行归一化后的值决定为该坐标上的特征向量，当所述梯度向量的大小小于所述既定阈值时，将0（零）向量决定为该坐标上的特征向量；
第2工序，按照所述第1及第2图像中互相对应的每个坐标，算出与计算该坐标上的特征向量彼此的内积的绝对值的N次方而得的值相当的相关值，求出包含按每个所述坐标算出的相关值的累计值的评价值，其中N为自然数；以及
第3工序，以使所述评价值变得更大的方式变更所述第2图像。
2.一种图像处理装置，其中包括：
决定单元，按照包含同一拍摄对象的第1及第2图像中的每个坐标，当该坐标上的像素值的梯度向量的大小为既定阈值以上时，将以该梯度向量的大小对所述梯度向量进行归一化后的值决定为该坐标上的特征向量，当所述梯度向量的大小小于所述既定阈值时，将0（零）向量决定为该坐标上的特征向量；
算出单元，按照所述第1及第2图像中互相对应的每个坐标，算出与计算该坐标上的特征向量彼此的内积的绝对值的N次方而得的值相当的相关值，求出包含按每个所述坐标算出的相关值的累计值的评价值，其中N为自然数；以及
变更单元，以使所述评价值变得更大的方式变更所述第2图像，
多次反复进行所述决定单元、算出单元及变更单元所进行的处理，进行所述第1及第2图像的对位。
3.如权利要求2所述的图像处理装置，其中所述梯度向量为图像中各坐标轴向的1次偏微分。
4.如权利要求2所述的图像处理装置，其中所述相关值是计算所述特征向量彼此的内积的平方而得的值。
5.如权利要求2所述的图像处理装置，其中所述评价值是以所述第1或第2图像的大小对所述相关值的累计值进行归一化而得的值。
6.如权利要求2所述的图像处理装置，其中所述变更包含平行移动、旋转及变形之中的至少一种。
7.如权利要求2所述的图像处理装置，其中所述决定单元基于对所述第1及第2图像分别进行构造物强调滤波处理和/或噪声降低处理而得到的图像，求出所述梯度向量。
8.如权利要求2所述的图像处理装置，其中所述第1及第2图像为医用图像。
9.如权利要求8所述的图像处理装置，其中所述第1及第2图像中摄像模式互相不同。
10.如权利要求9所述的图像处理装置，其中所述第1及第2图像的一个为超声波图像、即US图像。
11.如权利要求2所述的图像处理装置，其中反复进行所述决定单元、算出单元及变更单元所进行的处理，直至进行该处理的次数达到既定数，或直至所述评价值实质上收敛。
12.如权利要求2所述的图像处理装置，其中还包括推定所述第1及第2图像的至少一个的对象区域中的噪声量，并基于该噪声量决定所述阈值的单元。
13.如权利要求12所述的图像处理装置，其中所述对象区域为所述第1图像与所述第2图像之间共同的解剖学特征点的周边区域。
14.如权利要求12所述的图像处理装置，其中所述对象区域为所述至少一个图像的中央区域。
15.如权利要求12所述的图像处理装置，其中所述对象区域为病变部的周边区域。
16.如权利要求12所述的图像处理装置，其中所述对象区域为全部图像区域。
17.如权利要求12所述的图像处理装置，其中所述决定阈值的单元求出所述对象区域中的各像素上的梯度强度，并基于该梯度强度的直方图来决定所述阈值。
18.如权利要求17所述的图像处理装置，其中所述决定阈值的单元根据以下的（a）～（f）式的任一个或多个的组合的平均值来算出噪声量Nqt：
Nqt＝Mfmax＋σ×HWHML （a）
Nqt＝Mfmax＋σ×HWHMR （b）
Nqt＝Mfmax＋σ×(HWHML＋HWHMR)/2 （c）
Nqt＝σ×Mfmax （d）
Nqt＝σ×Mmom1 （e）
Nqt＝σ×Mmom2 （f）
其中，Mfmax是在所述直方图分布上出现的峰值之中、梯度强度最低的峰值给出众数值的梯度强度；
σ为任意常数，
HWHML是从所述直方图分布上的Mfmax观看低值侧的半值半宽度；
HWHMR是从所述直方图分布上的Mfmax观看高值侧的半值半宽度；
Mmom1是与所述直方图分布上的梯度强度从0到HWHMR的范围的重心相当的梯度强度；
Mmom2是与所述直方图分布上的梯度强度从HWHML到HWHMR的范围的重心相当的梯度强度。
CN201410162145.1A 2013-04-22 2014-04-22 图像处理方法及装置 Active CN104156940B (zh)

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
JP2013-088963 2013-04-22
JP2013088963 2013-04-22
JP2013-230466 2013-11-06
JP2013230466A JP5889265B2 (ja) 2013-04-22 2013-11-06 画像処理方法および装置並びにプログラム

Publications (2)

Publication Number Publication Date
CN104156940A CN104156940A (zh) 2014-11-19
CN104156940B true CN104156940B (zh) 2017-05-17

Family

ID=51729035

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201410162145.1A Active CN104156940B (zh) 2013-04-22 2014-04-22 图像处理方法及装置

Country Status (3)

Country Link
US (1) US9256916B2 (zh)
JP (1) JP5889265B2 (zh)
CN (1) CN104156940B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107004259A (zh) * 2014-12-10 2017-08-01 皇家飞利浦有限公司 多对比度成像中的统计加权正则化
CN107211090B (zh) * 2015-02-02 2020-02-28 富士胶片株式会社 操作装置、跟踪系统、操作方法及介质
JP6775294B2 (ja) * 2015-11-12 2020-10-28 ジーイー・メディカル・システムズ・グローバル・テクノロジー・カンパニー・エルエルシー 画像処理装置および方法並びにプログラム
WO2017102468A1 (en) 2015-12-15 2017-06-22 Koninklijke Philips N.V. Image processing systems and methods
US11350995B2 (en) 2016-10-05 2022-06-07 Nuvasive, Inc. Surgical navigation systems and methods
CN108280850A (zh) * 2018-03-05 2018-07-13 北京中科嘉宁科技有限公司 一种基于B样条和Levenberg-Marquardt优化的快速图像配准方法
US11430140B2 (en) * 2018-09-18 2022-08-30 Caide Systems, Inc. Medical image generation, localizaton, registration system
JP7290274B2 (ja) 2019-07-04 2023-06-13 東芝エネルギーシステムズ株式会社 荷電粒子の出射制御装置、方法及びプログラム
CN113627446B (zh) * 2021-08-18 2023-10-31 成都工业学院 基于梯度向量的特征点描述算子的图像匹配方法及系统
WO2023054089A1 (ja) * 2021-10-01 2023-04-06 ソニーグループ株式会社 映像処理システム、並びに、医療情報処理システムおよび動作方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102231191A (zh) * 2011-07-17 2011-11-02 西安电子科技大学 基于asift的多模态图像特征提取与匹配方法

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004227519A (ja) * 2003-01-27 2004-08-12 Matsushita Electric Ind Co Ltd 画像処理方法
JP4617883B2 (ja) * 2004-01-06 2011-01-26 ソニー株式会社 画像処理装置および方法、プログラム並びに記録媒体
US8090429B2 (en) 2004-06-30 2012-01-03 Siemens Medical Solutions Usa, Inc. Systems and methods for localized image registration and fusion
JP4559844B2 (ja) * 2004-12-27 2010-10-13 株式会社東芝 画像処理装置及び画像処理方法
US20070167784A1 (en) 2005-12-13 2007-07-19 Raj Shekhar Real-time Elastic Registration to Determine Temporal Evolution of Internal Tissues for Image-Guided Interventions
KR20070110965A (ko) 2006-05-16 2007-11-21 주식회사 메디슨 초음파 영상과 외부 의료영상의 합성 영상을디스플레이하기 위한 초음파 시스템
EP2081494B1 (en) 2006-11-16 2018-07-11 Vanderbilt University System and method of compensating for organ deformation
US8098911B2 (en) 2006-12-05 2012-01-17 Siemens Aktiengesellschaft Method and system for registration of contrast-enhanced images with volume-preserving constraint
IL191615A (en) * 2007-10-23 2015-05-31 Israel Aerospace Ind Ltd A method and system for producing tie points for use in stereo adjustment of stereoscopic images and a method for identifying differences in the landscape taken between two time points
EP2131212A3 (en) 2008-06-05 2011-10-05 Medison Co., Ltd. Non-Rigid Registration Between CT Images and Ultrasound Images
US8144947B2 (en) * 2008-06-27 2012-03-27 Palo Alto Research Center Incorporated System and method for finding a picture image in an image collection using localized two-dimensional visual fingerprints
GB0913930D0 (en) * 2009-08-07 2009-09-16 Ucl Business Plc Apparatus and method for registering two medical images
JP2013020294A (ja) * 2011-07-07 2013-01-31 Sony Corp 画像処理装置と画像処理方法およびプログラムと記録媒体
JP5591309B2 (ja) * 2012-11-29 2014-09-17 キヤノン株式会社 画像処理装置、画像処理方法、プログラム、及び記憶媒体

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102231191A (zh) * 2011-07-17 2011-11-02 西安电子科技大学 基于asift的多模态图像特征提取与匹配方法

Also Published As

Publication number Publication date
JP5889265B2 (ja) 2016-03-22
US20140314295A1 (en) 2014-10-23
CN104156940A (zh) 2014-11-19
US9256916B2 (en) 2016-02-09
JP2014225218A (ja) 2014-12-04

Similar Documents

Publication Publication Date Title
CN104156940B (zh) 图像处理方法及装置
JP4774501B2 (ja) 後処理を通じた２次元超音波映像の画質改善方法及び装置
CN104091337B (zh) 一种基于PCA及微分同胚Demons的变形医学图像配准方法
CN109949349B (zh) 一种多模态三维图像的配准及融合显示方法
Luo et al. Vision-based surgical field defogging
Xing et al. PDE-based spatial smoothing: a practical demonstration of impacts on MRI brain extraction, tissue segmentation and registration
JP4274400B2 (ja) 画像の位置合わせ方法および装置
Sindhwani et al. Semi‐automatic outlining of levator hiatus
JP4708740B2 (ja) 画像処理装置及び画像処理方法
CN103295195A (zh) 低对比度图像的脉管增强的方法及其系统
Fourcade et al. Deformable image registration with deep network priors: a study on longitudinal PET images
Hsia et al. A 3D endoscopic imaging system with content-adaptive filtering and hierarchical similarity analysis
JP2004302581A (ja) 画像処理方法並びに画像処理装置
Qian et al. A cascade-network framework for integrated registration of liver DCE-MR images
US8139830B2 (en) System and method for automated alignment of leg volumes in whole-body magnetic resonance scans
Amir-Khalili et al. Real-time extraction of local phase features from volumetric medical image data
CN106570837A (zh) 基于高斯多尺度空间的脑组织mri图像偏移场校正方法
Li et al. A medical image registration method based on weighted mutual information
Mehr et al. Deep Learning-Based Ultrasound Image Despeckling by Noise Model Estimation.
Chen et al. Unsupervised image-to-image translation in multi-parametric MRI of bladder cancer
KR101341576B1 (ko) 등고선 기반 관심영역 결정방법 및 장치
Chen et al. Distance transforms in multichannel MR image registration
Qiu et al. Stabilization algorithm based on improved motion model for jittery video in minimally invasive surgery
Ghosh¹ et al. Check for updates Enhancing Sustainable Development in Medical Image Fusion Using Proposed (EBCE) with DWT
Kumar et al. Motion correction strategies for interventional angiography images: a comparative approach

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