CN102646265B - 图像处理设备和方法 - Google Patents

图像处理设备和方法 Download PDF

Info

Publication number
CN102646265B
CN102646265B CN201110045243.3A CN201110045243A CN102646265B CN 102646265 B CN102646265 B CN 102646265B CN 201110045243 A CN201110045243 A CN 201110045243A CN 102646265 B CN102646265 B CN 102646265B
Authority
CN
China
Prior art keywords
image
filtering
image processing
iterative
processing equipment
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
CN201110045243.3A
Other languages
English (en)
Other versions
CN102646265A (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.)
Canon Medical Systems Corp
Original Assignee
Toshiba Corp
Toshiba Medical Systems Corp
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 Toshiba Corp, Toshiba Medical Systems Corp filed Critical Toshiba Corp
Priority to CN201110045243.3A priority Critical patent/CN102646265B/zh
Priority to JP2012034225A priority patent/JP5925515B2/ja
Priority to US13/402,435 priority patent/US9202261B2/en
Publication of CN102646265A publication Critical patent/CN102646265A/zh
Application granted granted Critical
Publication of CN102646265B publication Critical patent/CN102646265B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • 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

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)
  • Measuring And Recording Apparatus For Diagnosis (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本公开提供了图像处理设备和方法。该图像处理设备可以包括：滤波装置，用于对输入图像进行迭代滤波；以及迭代滤波停止装置，用于根据所述迭代滤波过程中每次滤波后得到的经滤波图像随迭代次数的变化速率来确定是否停止所述迭代滤波。

Description

图像处理设备和方法
技术领域
本公开涉及图像处理领域，具体而言，涉及用于对图像(如医学图像)降噪的图像处理设备和方法。
背景技术
图像降噪是图像处理的一个重要环节。例如，在医学图像处理中，由于设备的信号检测性能、成像算法不适当、扫描剂量低等诸如此类的原因，利用计算机断层扫描(Computed Tomography，CT)设备、核磁共振成像(Magnetic Resonance Imaging，MRI)设备、超声波(Ultrasound，UL)诊断装置等医疗设备获取的医学图像常常会包含噪声。这些图像中的噪声会严重影响医生的诊断，并且会影响后续图像处理(如图像分割、测量和体绘制(Volume Rending)等)的正确性。
发明内容
在下文中给出关于本公开的一些方面的简要概述，以便提供对于本公开的基本理解。该概述并非意图确定本公开的关键或重要部分，也不是意图限定本公开的范围。其目的仅仅是以简化的形式给出某些概念，以此作为稍后论述的更详细描述的前序。
根据本公开的一个方面，提供了一种图像处理设备。该图像处理设备可以包括：滤波装置，用于对输入图像进行迭代滤波；以及迭代滤波停止装置，用于根据所述迭代滤波过程中每次滤波后得到的经滤波图像随迭代次数的变化速率来确定是否停止所述迭代滤波。
根据本公开的另一方面，提供了一种图像处理方法。该图像处理方法可以包括：对输入图像进行迭代滤波；以及根据所述迭代滤波过程中每次滤波后得到的经滤波图像随迭代次数的变化速率来确定是否停止所述迭代滤波。
另外，本公开的实施例还提供了用于实现上述方法的计算机程序。
此外，本公开的实施例还提供了至少计算机可读介质形式的计算机程序产品，其上记录有用于实现上述方法的计算机程序代码。
附图说明
参照下面结合附图对本公开实施例的说明，会更加容易地理解本公开的以上和其它目的、特点和优点。附图中的部件只是为了示出本公开的原理。在附图中，相同的或类似的技术特征或部件将采用相同或类似的附图标记来表示。
图1是示出根据本公开的一个实施例的图像处理设备的结构的示意性框图；
图2是示出根据该实施例的图像处理方法的示意性流程图；
图3是示出根据本公开的一个具体实施例的图像处理设备的结构的示意性框图；
图4是示出根据该具体实施例的图像处理方法的示意性流程图；
图5是示出用于确定停止迭代滤波的时间的一个具体示例的示意性流程图；
图6(A)、6(B)和6(C)是分别示出对不同类型的图像进行迭代滤波的过程中经滤波图像与噪声图像的相关系数的变化曲线的示意图；
图7是示出根据本公开的实施例的应用场景的一个示例的示意图；
图8是示出根据本公开的另一具体实施例的图像处理设备的结构的示意性框图；
图9是示出根据本公开的另一具体实施例的图像处理设备的结构的示意性框图；
图10是示出图8所示的实施例的图像处理方法的示意性流程图；
图11是示出图9所示的实施例的图像处理方法的示意性流程图；
图12(A)、12(B)是分别示出在迭代滤波过程中经滤波图像和噪声图像的方差的变化曲线的示意图；
图13是示出根据本公开的另一具体实施例的图像处理设备的结构的示意性框图；
图14是示出根据该具体实施例的图像处理方法的示意性流程图；
图15是示出该具体实施例的一个具体应用示例的示意图；
图16是示出根据本公开的另一实施例的图像处理设备的结构的示意性框图；
图17是示出根据该实施例的图像处理方法的示意性流程图；
图18是示出了图像中的各个图像单元的梯度值分布的示例的示意图；及
图19是示出用于实现本公开的计算机的结构的示例性框图。
具体实施方式
下面参照附图来说明本公开的实施例。在本公开的一个附图或一种实施方式中描述的元素和特征可以与一个或更多个其它附图或实施方式中示出的元素和特征相结合。应当注意，为了清楚的目的，附图和说明中省略了与本公开无关的、本领域普通技术人员已知的部件和处理的表示和描述。
在进行图像处理时，通常可以对图像进行迭代滤波，以尽量滤除图像中的噪声。如果迭代的次数太少，则降噪效果可能不理想，而另一方面，如果迭代次数过多，则图像中的有用信息也会被平滑，从而导致图像模糊，影响后续的处理。本公开的一些实施例提供了能够在图像的迭代滤波过程中确定停止迭代滤波的适当时机的图像处理设备和方法。
应理解，根据本公开的实施例的设备和方法可用于对各种类型的图像进行降噪处理。例如，根据本公开的实施例的设备和方法可用于医学图像的降噪处理。这里所述的医学图像可以是根据利用医疗诊断成像装置获得的被检测者的数据而形成的图像。这里所述的医疗诊断装置包括但不限于：X射线成像诊断装置、超声波(UL)诊断成像装置、计算机断层扫描(CT)装置、磁共振成像(MRI)诊断成像装置、正电子发射断层扫描(Positron Emission Tomography，PET)装置等。图7是示出了利用医疗诊断成像装置获取医学图像的示意图。如图7所示，对于每种成像设备(例如上述医疗诊断装置)均有一些预定的或用户自定义的扫描协议。以CT为例，可以有诸如心电同步扫描(cardiac synchronization)、腹部4相增强扫描(abdominal 4 phases enhance)、动态扫描等协议，如图7中所示的协议1、协议2、......、协议N(N≥1)；利用某个协议，可以获得多个图像序列(如图7中所示的序列1、序列2、......、序列N等)。因此，待处理的图像可以是根据利用医疗诊断装置、基于相应的协议获取的患者数据所形成的医学图像。
图1是示出根据本公开的一个实施例的图像处理设备的结构的示意性框图，而图2是示出根据该实施例的图像处理方法的示意性流程图。
如图1所示，图像处理设备100可以包括滤波装置101和迭代滤波停止装置103。图像处理设备100可以采用如图2所示的方法进行图像滤波。下面参考图2所示的方法流程来描述图像处理设备100中的各个装置的功能。
滤波装置101用于对输入到该图像处理设备的输入图像进行迭代滤波(图2的步骤202)。根据实际应用场景，滤波装置101(即步骤202)可以采用任何适当的滤波方法，如非线性扩散(Non-Linear Diffusion，NLD)滤波方法、线性滤波方法、各向异性扩散(Anisotropic Diffusion，ATD)滤波方法等等，这里不作限定。
滤波装置101在每次滤波后将得到的经滤波图像输出到迭代滤波停止装置103，而迭代滤波停止装置103用于根据每次滤波后得到的经滤波图像来判断是否停止迭代滤波。具体地，迭代滤波停止装置103可以根据迭代滤波过程中经滤波图像随着迭代次数的增加而变化的速率来判断是否停止迭代滤波(图2的步骤206)。例如，如果迭代滤波停止装置103判断经滤波图像的变化速率下降为小于或等于预定的阈值，这就说明此时的滤波操作对图像的影响变得足够小，也就是说可被滤除的噪声变得足够少了，则可以指示滤波装置101停止迭代滤波，否则，则指示滤波装置101继续进行下一次迭代。
这里所述的经滤波图像的“变化”是指当前滤波得到的经滤波图像与之前滤波(“之前滤波”是指当前迭代之前进行的某一次迭代滤波(如上一次迭代等)，在此不作限定)得到的经滤波图像之间的差异。在对图像的迭代滤波过程中，随着迭代次数的增加，每次滤波得到的经滤波图像与之前滤波得到的经滤波图像相比的变化(即差异)会随着迭代次数而变小。也就是说，随着迭代次数的增加，经滤波图像的变化速率会变慢。图1和图2所示的设备和方法基于这一原理，根据经滤波图像的变化速率来确定停止迭代滤波的时间。利用这种方法或设备，可以有效地确定停止迭代滤波的适当时间，从而在滤除图像中的大部分噪声的同时保证图像中的有用信息不被模糊，另外，还可避免手工设置停止迭代滤波时间的麻烦，从而提高图像处理的自动化程度。
经滤波图像在迭代过程中的变化可以反映于其不同特征或特性，如每次迭代得到的经滤波图像与噪声图像的相关性在迭代过程中的变化、经滤波图像或噪声图像的某些图像特征(如方差或标准差)的值在迭代过程中的变化等等，但不限于此，例如还可以是上述各种变化的数学组合等。下文中将参考图3至图15来描述在迭代滤波过程中根据经滤波图像的变化速率来判断是否停止迭代滤波的图像处理设备和方法的一些更为具体的实施例。
图3是示出根据一个具体实施例的图像处理设备的结构的示意性框图，而图4是示出根据该具体实施例的图像处理方法的示意性流程图。
如上所述，在迭代滤波过程中，每次滤波得到的经滤波图像与之前滤波得到经滤波图像相比，其变化会变小，即变化速率会变慢。在迭代滤波的开始阶段，每次滤波滤除的噪声较多，而随着迭代次数的增加，图像中的噪声越来越少，而图像本身也逐渐平滑。滤波图像的变化速率与滤波图像的梯度成正比，随着图像本身逐渐平滑，图像的梯度随之变小，因此，经滤波图像的变化速率会变慢。相应地，随着经滤波图像中噪声含量的变化，其与所滤除的噪声图像之间的相关性也会发生变化。也就是说，所述相关性的变化也反映了经滤波图像的变化，相关性的变化速率也反映了经滤波图像的变化速率。图6(A)、6(B)和6(C)分别示出了这种相关性在迭代滤波过程中的示例性变化曲线。图6(A)、6(B)和6(C)中的曲线分别表示在对三个不同图像(在所图示的示例中，三个图像分别为采用UL、CT和MRI设备针对不同人体部位而获取的医学图像)进行迭代滤波过程中每次滤波得到的经滤波图像与相应的噪声图像的相关系数随时间(迭代次数)的变化趋势。如图6(A)和6(B)所示，在迭代滤波的开始阶段，经滤波图像与噪声图像的相关系数的值呈下降趋势，这表示越来越多的噪声被滤除。经过一段时间之后，相关系数的曲线到达谷底(图中用“×”表示的谷点T1、T2)。在谷点之后，经滤波图像与噪声的相关系数又呈现缓慢上升的趋势，这表示在谷点之后，随着迭代次数的增加，图像中的有用信息也被逐渐平滑掉，也就是说一些有用信息也被滤除，或者说进入了噪声图像，从而使得经滤波图像与滤除的噪声图像之间的相关性又有所增大。一种已有的方法即通过寻找上述谷点来确定停止迭代滤波的时间。在这种已有的方法中，将经滤波图像与噪声图像的相关系数最小的时间点确定为停止迭代滤波的时间点，相关描述可以参见Pavel Mrázek的论文“Selection of Optimal Stopping Time for Nonlinear Diffusion Filtering”(刊于“Scale-Space and Morphology in Computer Vision”，第三次国际会议，2001年)(下文中将该论文称为相关文献1)。在图6(A)所示的情况下，存在一个谷点T1，且该谷点的位置对应于正确的停止时间，因此，在这种情况下，利用上述已有方法可以确定正确的停止时间。但是，上述已有方法并不适用于图6(B)和6(C)所示的情况：例如，在图6(B)所示的曲线中，出现谷点T2的位置明显滞后于正确的停止时间，且在谷点T2之后，相关系数曲线保持平稳(即在谷点T2之后，相关系数值几乎没有变化)，因此，在这种情况下，如果利用上述已有方法，则所确定的时间点为谷点T2，而不是正确的停止时间；又如，图6(C)所示的曲线是单调下降的，没有谷点，因此，在这种情况下，如果利用上述已有方法，则无法找到相关系数最小的时间点，从而无法确定正确的停止时间。
本申请的发明人发现，如6(A)、6(B)和6(C)的三个曲线所示，不管是否有谷点，这三个曲线在迭代滤波的开始阶段都是单调下降的，且下降的速率越来越慢。图3和图4所示的设备和方法即利用经滤波图像与噪声图像之间的相关性在迭代滤波过程中下降的速率(而不是相关性最小的时间点)来确定停止迭代滤波的时间。
如图3所示，图像处理设备300可以包括滤波装置301和迭代滤波停止装置303，其中迭代滤波停止装置303可以包括噪声图像获取单元303-1、相关性获取单元303-2和停止判断单元303-3。图像处理设备300可以采用如图4所示的方法流程进行图像滤波。下面参考图4所示的方法流程来描述图像处理设备300中的各个装置和/或单元的功能。
滤波装置301与图1所示的滤波装置101相似，用于采用任何适当的方法对输入到该图像处理设备的输入图像进行迭代滤波(图4的步骤402)，这里不再重复。
滤波装置301在每次滤波后将得到的经滤波图像输出到噪声图像获取单元303-1。噪声图像获取单元303-1用于在每次滤波后根据所得到的经滤波图像和输入图像来估算噪声图像(图4中的步骤406-1)。
作为一个示例，可以计算每次滤波后得到的经滤波图像与输入图像之间的差值图像，作为对应的噪声图像。设Io表示输入图像，Ist表示当前滤波后得到的经滤波图像，则对应的噪声图像Int可以通过下式来估算：
Int＝Io-Ist                            (1)
应理解，每次滤波后，可以采用任何适当的其他方法(如取根据式(1)计算得到的差值图像的绝对值，即Int＝|Io-Ist|，这里不一一列举)来估算噪声图像，包括对上式(1)进行任何形式的数学变化，因此本公开并不局限于上述示例。
噪声图像获取单元303-1将得到的噪声图像输出到相关性获取单元303-2，而相关性获取单元303-2可以估算该噪声图像与当前迭代得到的经滤波图像的相关性(图4的步骤406-2)，并将估算结果输出到停止判断单元303-3。
作为一个示例，可以计算噪声图像与经滤波图像之间的相关系数，以表征两个图像之间的相关性。例如，可以采用下面的公式来计算两个图像之间的相关系数：
ρ t = E { [ I nt - E ( I nt ) ] [ I st - E ( I st ) ] } E [ I nt - E ( I nt ) ] 2 · E [ I st - E ( I st ) ] 2 - - - ( 2 )
其中，ρt表示当前迭代得到的经滤波图像Ist与对应的噪声图像Int之间的相关系数，E(·)表示均值函数。
应理解，还可以采用除相关系数之外的其他参数来表征两个图像之间的相关性，包括对上式(2)进行任何形式的数学变化，因此本公开并不局限于上述示例。
停止判断单元303-3用于根据每次滤波得到的经滤波图像与噪声图像之间的相关性在所述迭代滤波过程中下降的速率来确定是否停止迭代滤波(图4中的步骤406-3)。如果停止判断单元303-3判断相关性的下降速率已经小于或等于预定的值，则可以指示滤波装置301停止迭代滤波，否则，指示滤波装置301继续进行下一次迭代。
图5示出了根据表征相关性下降速率的特征值来确定迭代停止时间的一个示例。如图5所示，在步骤406-31中，停止判断单元303-3根据当前滤波得到的经滤波图像与噪声图像之间的相关性的值以及之前滤波后估算的相关性的值来计算用于表示相关性随迭代次数的增加而下降的速率的特征值。然后，在步骤406-32中，停止判断单元303-3判断该特征值与预定的阈值之间是否满足预定的条件，若是，则停止迭代滤波，否则，继续进行下一次迭代。
作为一个具体示例，所述特征值可以是经滤波图像与噪声图像的相关系性随迭代次数的变化曲线(如图6(A)、6(B)或6(C)所示的相关系数曲线)的倾角。例如，可以用下式来计算所述倾角θ：
θ = arctan ρ t - ρ t - Δt Δt · 180 π - - - ( 3 )
其中，ρt表示在时间t得到的经滤波图像Ist与对应的噪声图像Int之间的相关系数，ρt-Δt表示在时间t-Δt得到的经滤波图像Is(t-Δt)与对应的噪声图像In(t-Δt)之间的相关系数，Δt表示时间间隔。
作为另一具体示例，所述倾角还可以用弧度来表示：
θ rad = arctan ρ t - ρ t - Δt Δt - - - ( 4 )
当停止判断单元303-3确定计算得到的倾角θrad或θ小于或等于预定的阈值时，则指示滤波装置301停止迭代滤波，否则，指示滤波装置301继续进行下一次迭代。
作为另一具体示例，所述特征值可以是经滤波图像与噪声图像的相关系性随迭代次数的变化曲线(如图6(A)、6(B)或6(C)所示的相关系数曲线)的斜率。例如，可以用下式来计算所述斜率kslope：
k slope = ρ t - ρ t - Δt Δt - - - ( 5 )
其中，ρt表示在时间t得到的经滤波图像Ist与对应的噪声图像Int之间的相关系数，ρt-Δt表示在时间t-Δt得到的经滤波图像Is(t-Δt)与对应的噪声图像In(t-Δt)之间的相关系数，Δt表示时间间隔。
当停止判断单元303-3确定计算得到的斜率kslope小于或等于预定的阈值时，则指示滤波装置301停止迭代滤波，否则，指示滤波装置301继续进行下一次迭代。
应注意，用于表征经滤波图像与噪声图像之间的相关性在所述迭代滤波过程中下降的速率的特征值并不局限于上述相关系数曲线的倾角和斜率。可以采用能够表征相关性下降速率的任何其他特征值，作为所述特征值。例如，可以采用当前滤波后估算的相关系数ρt与之前滤波后估算的相关系数ρt-Δt之间的差(或该差的倒数)作为所述特征值；当判断上述差小于或等于预定的阈值时(或者当上述差的倒数大于或等于预定的阈值时)，停止迭代滤波，否则继续进行下一次迭代。此外，对上述各种特征值的计算方式还可以进行任何形式的数学变化。
另外，应理解，上述实施例/示例中所述的预定的阈值是根据实际应用场景来预先确定的，其可以是通过对不同类型(或特性)的图像进行分析而获得的实验值或经验值。针对不同的应用场景以及所采用的不同特征值参数，该阈值可以采取不同的值。例如，在所述输入图像为根据通过医疗诊断装置得到的数据而形成的医学图像的情况下，可以根据实验或经验针对每种扫描协议或成像设备的类型而预先确定相应的阈值。可选地，可以将针对多个扫描协议或成像设备的多个阈值(例如以配置文件等形式)保存在存储单元(该存储单元可以是内置于图像处理设备中的存储装置，也可以是位于图像处理设备外部并可由该图像处理设备访问的存储装置)中。这样，在迭代滤波过程中，可以根据输入图像的扫描协议或成像设备的类型来获取(如从存储单元中读取)相应的阈值，作为判断相关性下降速率的依据，从而进一步提高图像处理的自动化程度。
在参考图3-5描述的图像处理设备或方法中，通过利用经滤波图像与噪声图像的相关性下降的速率，能够有效地确定停止迭代滤波的时间，在滤除图像中的大部分噪声的同时，能够保证图像中的有用信息不被平滑掉。
图8和图9是分别示出根据另外两个具体实施例的图像处理设备的结构的示意性框图，图10和图11是分别根据这两个具体实施例的图像处理方法的示意性流程图。参考图3至图5描述的实施例利用经滤波图像与噪声图像的相关性下降的速率来表征迭代滤波过程中经滤波图像的变化速率。在图8至图11所示的实施例中，利用经滤波图像或噪声图像的预定图像特征值在迭代滤波过程中的变化速率来表征经滤波图像的变化速率。
这里所述的预定图像特征可以是能够反映在迭代滤波过程中经滤波图像的变化的任何图像特征，例如经滤波图像或噪声图像的方差或标准差、经滤波图像或噪声图像的像素的平均梯度、经滤波图像或噪声图像的灰度的动态范围(即最大灰度与最小灰度的差)等，这里不一一列举。
图12(A)、12(B)分别示出了在迭代滤波过程中经滤波图像和噪声图像的方差的示例性的变化曲线。如图12(A)所示，随着迭代次数的增加，经滤波图像的方差单调下降，且下降速率变慢。如图12(B)所示，随着迭代次数的增加，噪声图像的方差单调上升，且上升速率变慢。可以利用上文中描述的方法来估算噪声图像，这里不再重复。
下面，参考图8和图10来描述利用经滤波图像的预定图像特征值的变化速率来确定滤波停止时间的实施例。
如图8所示，图像处理设备800可以包括滤波装置801和迭代滤波停止装置803，其中，迭代滤波停止装置803可以包括图像特征计算单元803-4和停止判断单元803-3。图像处理设备800可以采用如图10所示的图像处理方法进行图像滤波。
与上文所述的滤波装置101和301相似，滤波装置801用于对输入到该图像处理设备的输入图像进行迭代滤波(图10的步骤1002)，这里不再重复。
滤波装置801在每次滤波后将得到的经滤波图像输出到图像特征计算单元803-4，由图像特征计算单元803-4来计算每次滤波得到的经滤波图像的预定图像特征值(如图10所示的步骤1006-1)，例如，计算经滤波图像的方差或标准差，并将计算结果输出到停止判断单元803-3。停止判断单元803-3根据在迭代滤波过程中所述图像特征(如方差或标准差)值随迭代次数的增加而变化的速率(如方差或标准差下降的速率)来确定是否停止迭代滤波。如果停止判断单元803-3判断所述图像特征值的变化速率小于或等于预定的阈值，则可以指示滤波装置801停止迭代滤波，否则，指示滤波装置801继续进行下一次迭代。
作为具体示例，停止判断单元803-3可以根据图像特征计算单元803-4的计算结果来计算用于表示预定图像特征值在迭代滤波过程中的变化速率的特征值。例如，该特征值可以是反映预定图像特征值随着迭代次数的增加而变化的曲线的倾角或斜率(例如，图12(A)所示的曲线的倾角或斜率，可以采用上述公式(3)-(5)来计算，这里不再重复)；或者，该特征值可以是当前滤波得到的经滤波图像的预定图像特征值与之前滤波得到的经滤波图像的预定图像特征值之间的差(或该差的倒数)，这里不一一列举。然后，停止判断单元803-3判断计算得到的特征值与预定的阈值之间是否满足预定的关系(如，在使用倾角、斜率或上述差作为特征值的情况下，判断特征值是否小于或等于预定的阈值；在使用上述差的倒数作为特征值的情况下，判断特征值是否大于或等于预定的阈值)，若是，则指示滤波装置801停止迭代滤波，否则，指示滤波装置801继续进行下一次迭代。与上述实施例/示例中相似，这里所述的阈值是根据实际应用场景来预先确定的。针对不同的应用场景以及所采用的不同特征参数，该阈值可以采取不同的值(例如，可以通过如图12(A)所示的图像特征值的变化曲线来预先确定阈值，这些曲线可以通过对不同类型(或特性)的图像进行分析而获得，这里不作详述)，这里不限定其具体数值。
下面，参考图9和图11来描述利用噪声图像的预定图像特征值的变化速率来确定滤波停止时间的实施例。
如图9所示，图像处理设备900可以包括滤波装置901和迭代滤波停止装置903。与图8所示的实施例相似，迭代滤波停止装置903也包括图像特征计算单元903-4和停止判断单元903-3，不同之处在于，迭代滤波停止装置903还包括噪声图像获取单元903-1。图像处理设备900可以采用如图11所示的图像处理方法进行图像滤波。
与上文所述的滤波装置101和301等相似，滤波装置901用于对输入到该图像处理设备的输入图像进行迭代滤波(图11的步骤1102)，这里不再重复。
滤波装置901在每次滤波后将得到的经滤波图像输出到噪声图像获取单元903-1。噪声图像获取单元903-1用于在每次滤波后根据所得到的经滤波图像和输入图像来估算噪声图像(图11中的步骤1106-1)。噪声图像获取单元903-1可以采用上文中参考图3和图4所描述的方法来估算噪声图像，这里不再重复。
噪声图像获取单元903-1将得到的噪声图像输出到图像特征计算单元903-4，由图像特征计算单元903-4来计算噪声图像的预定图像特征值(如图11所示的步骤1106-2)，例如，计算噪声图像的方差或标准差，并将计算结果输出到停止判断单元903-3。
停止判断单元903-3根据在迭代滤波过程中噪声图像的图像特征值(如方差或标准差)随迭代次数的增加而变化的速率(如方差或标准差上升的速率)来确定是否停止迭代滤波。如果停止判断单元903-3判断噪声图像的所述图像特征值的变化速率达到预定的程度，则可以指示滤波装置901停止迭代滤波，否则，指示滤波装置901继续进行下一次迭代。
作为具体示例，停止判断单元903-3可以根据图像特征计算单元903-4的计算结果来计算用于表示预定图像特征值在迭代滤波过程中的变化速率的特征值。例如，该特征值可以是反映预定图像特征值随着迭代次数的增加而变化的曲线的倾角或斜率(例如，图12(B)所示的曲线的倾角或斜率，可以采用上述公式(3)-(5)来计算，这里不再重复)；或者，该特征值可以是当前滤波得到的噪声图像的预定图像特征值与之前滤波得到的噪声图像的预定图像特征值之间的差(或该差的倒数)，这里不一一列举。然后，停止判断单元903-3判断计算得到的特征值与预定的阈值之间是否满足预定的关系(如，在使用倾角、斜率或上述差作为特征值的情况下，判断特征值是否小于或等于预定的阈值；在使用上述差的倒数作为特征值的情况下，判断特征值是否大于或等于预定的阈值)，若是，则指示滤波装置901停止迭代滤波，否则，指示滤波装置901继续进行下一次迭代。与上述实施例/示例中相似，这里所述的阈值是根据实际应用场景来预先确定的。针对不同的应用场景和所采用的不同特征参数，该阈值可以采取不同的值(例如，可以通过如图12(B)所示的图像特征值的变化曲线来预先确定阈值，这些曲线可以通过对不同类型(或特性)的图像进行分析而获得，这里不作详述)，这里不限定其具体数值。
在参考图8-11描述的图像处理设备或方法中，通过利用迭代滤波过程中经滤波图像或噪声图像的预定图像特征值的变化速率，能够有效地确定停止迭代滤波的时间，在滤除图像中的大部分噪声的同时，能够保证图像中的有用信息不被平滑掉。
作为具体示例，参考图8-11描述的图像处理设备或方法可以采用下式来计算图像的方差，作为所述预定图像特征值：
D = E [ I - E ( I ) ] 2 - - - ( 6 )
其中，I表示图像(如上述的经滤波图像或噪声图像)，D表示该图像的方差。
图13是示出根据另一具体实施例的图像处理设备的结构的示意性框图，图14是示出根据该具体实施例的图像处理方法的示意性流程图。上文描述的实施例分别利用噪声图像与经滤波图像之间的相关性或者噪声图像(或经滤波图像)的预定图像特征值在迭代滤波过程中的变化速率，而图13和图14所示的实施例则考虑将二者组合起来。如图6(A)、6(B)、6(C)所示，噪声图像与经滤波图像之间的相关性在迭代滤波的开始阶段单调下降，之后会缓慢上升(如图6(A)和6(B)所示)或者继续缓慢下降(如图6(C)所示)；而如图12(B)所示，噪声图像的方差在迭代过程中是单调上升的。如果将二者以适当的方式组合起来(下文中将参考图13-15描述将二者组合起来的具体实施例)，则所产生的组合特征值随着迭代次数的变化曲线也会反映迭代滤波过程中图像的变化。例如，在图6(B)所示的相关性曲线谷点延后的情况下，由于噪声图像的方差的单调上升，组合特征值的曲线会使谷点提前到适当的位置。又如，在图6(C)所示的相关性单调下降的情况下，由于噪声图像的方差的单调上升，组合特征值的曲线会出现谷点，而且该谷点对应正确的停止时间。图15(A)是示出了噪声图像与经滤波图像的相关系数、噪声图像的方差以及二者的组合特征值随迭代次数的变化的示意图。在图15(A)中，纵轴表示相关系数、方差及组合特征的取值，横轴表示迭代时间。如图15(A)所示，组合特征值的曲线在t＝0.4的时间点出现了谷点(极小值)。利用所述组合特征值随迭代次数的增加而出现的谷点(即极小值)，可以确定停止迭代滤波的时间。
为了方便，下文的描述以噪声图像的方差为所述预定图像特征、噪声图像与经滤波图像的相关系数表征所述相关性为例。应理解，在其他示例中，可以采用噪声图像的其他图像特征(如标准差等)或者采用经滤波图像的方差或标准差的倒数等，也可以采用反映经滤波图像与噪声图像之间的相关性的其他参数，这里不一一列举。
如图13所示，图像处理设备1300可以包括滤波装置1301和迭代滤波停止装置1303，其中，迭代滤波停止装置1303包括噪声图像获取单元1303-1、相关性获取单元1303-2、图像特征计算单元1303-4、组合单元1303-5和停止判断单元1303-3。图像处理设备1300可以采用如图14所示的图像处理方法进行图像滤波。
与上文所述的滤波装置101和301等相似，滤波装置1301用于对输入到该图像处理设备的输入图像进行迭代滤波(图14的步骤1402)，这里不再重复。
滤波装置1301在每次滤波后将得到的经滤波图像输出到噪声图像获取单元1303-1。噪声图像获取单元1303-1用于在每次滤波后根据所得到的经滤波图像和输入图像来估算噪声图像(图14中的步骤1406-1)。噪声图像获取单元1303-1可以采用上文中参考图3和图4所描述的方法来估算噪声图像，这里不再重复。
噪声图像获取单元1303-1将得到的噪声图像分别输出到相关性获取单元1303-2和图像特征计算单元1303-4。相关性获取单元1303-2估算该噪声图像与当前滤波得到的经滤波图像的相关性(图14的步骤1406-2)，并将估算结果输出到组合单元1303-5。例如，相关性获取单元1303-2可以计算噪声图像与经滤波图像之间的相关系数(如采用上文描述的方法，这里不再重复)，以表征两个图像之间的相关性。
图像特征计算单元1303-4计算噪声图像的预定图像特征值(如图14所示的步骤1406-4)，例如，计算噪声图像的方差(如采用上式(6))，并将计算结果输出到组合单元1303-5。
组合单元1303-5将相关性获取单元1303-2获得的相关性值与图像特征计算单元1303-4计算得到的方差值组合起来，以生成组合特征值(如图14所示的步骤1406-5)。例如，如果用C1表示单调上升的噪声图像的方差曲线，C2表示噪声图像与经滤波图像的相关性曲线，则组合特征值可表示为：
Ccombine＝f(C1)+g(C2)                        (7)
其中，Ccombine表示组合特征值，f(·)表示对曲线C1所做的任何数学变化，而g(·)表示对曲线C2所做的任何数学变化。
作为一个示例，组合单元1303-5可以计算相关性获取单元1303-2获得的相关性值与图像特征计算单元1303-4计算得到的方差值求和，并将该和作为组合特征值。但是，在一些情况下，经滤波图像与噪声图像的相关系数的取值范围往往不同于噪声图像的方差的取值范围，因此，相关性值与方差值的和曲线有可能没有上述谷点(极小值)，因此，需要进行调整，以使两个曲线相匹配。下文描述将相关性获取单元1303-2获得的相关性值与图像特征计算单元1303-4计算得到的方差值进行组合的一些示例。
作为一个示例，在进行滤波之前，图像处理设备(如滤波装置1301)可以首先将原始输入图像归一化，换言之，在图14所示的步骤1402之前，可以包括将输入图像归一化的步骤(图中未示出)。这样，在每次迭代后，噪声图像获取单元1303-1基于经过归一化的输入图像以及经滤波图像来估算噪声图像。例如，组合单元1303-5可以将利用归一化的输入图像而得到的所述相关性值(来自相关性获取单元1303-2)与所述方差值(来自图像特征计算单元1303-4)求和，并将该和作为组合特征值。
作为另一示例，图像特征计算单元1303-4在每次迭代后可以对噪声图像获取单元1303-1估计的噪声图像进行归一化处理(在该示例中，无需对原始输入图像进行归一化)，计算归一化的噪声图像的方差，并将归一化的噪声图像的方差输出到组合单元1303-5。可以通过将噪声图像的动态范围(例如，可以用[Nmin，Nmax]表示图像的动态范围)映射到[0，1]来进行归一化。针对不同类型的输入图像，每次迭代产生的噪声图像的动态范围不同，使得归一化不是基于一个相同的动态范围来进行，因此，可以指定一个标准动态范围，在每次迭代时都以标准动态范围映射到[0，1]之间。例如，针对某种类型的输入图像，可以设定其相应的标准动态范围为[Smin，Smax](应理解，针对不同的输入图像，可以采用不同的标准动态范围)。归一化运算可以用下式(8)来表示：
N ′ = N S max - S min + S min - - - ( 8 )
其中，N表示噪声图像的像素值，N’表示归一化后的像素值。
应理解，在实际应用中，可以根据实际需要来确定上述标准动态范围，而不应将本公开限于具体的示例。例如，在医学图像的情况下，标准动态范围可以根据获取输入图像的医疗诊断装置的不同类型来确定，或者可以根据输入图像的扫描协议等来确定，这里不一一列举。
在该示例中，相关性获取单元1303-2可以基于噪声图像获取单元1303-1估计的噪声图像(而不是经过归一化的噪声图像)来计算所述相关性。
作为另一示例，组合单元1303-5可以对计算得到的噪声图像的方差进行变换，然后将经过变换的方差与相关性获取单元1303-2得到的相关性值求和，并将该和作为组合特征值。例如，组合单元1303-5可以对图像特征计算单元1303-4计算得到的噪声图像的方差进行变换处理，然后将变换后的方差值与相关性获取单元1303-2得到的相关性值相加，将得到的和作为所述组合特征值。例如，组合单元1303-5可以通过下式对图像特征计算单元1303-4计算得到的噪声图像的方差进行变换：
V i ′ = a i ( V i - offset ) - - - ( 9 )
其中，Vi表示图像特征计算单元1303-4计算得到的噪声图像的方差；表示组合单元1303-5对方差值Vi进行变换后得到的变换值；i表示迭代时间，i＝1，2，...；offset表示平移系数，且V1表示第一次迭代后图像特征计算单元1303-4计算得到的噪声图像的方差；ai为随迭代次数增加而增大的系数，该系数大于1，例如，ai＝1.0+c×i；c表示大于零的预定常数。在式(9)中，通过求方差值Vi的平方根，使得变换后的方差值变大；通过平移系数offset，将变换后的方差曲线平移，使其从0开始；通过系数ai，增大了变换后的方差曲线的上升速率，且随着迭代次数的增加，该上升速率会增大。图15(B)是示出了经过上述变换的噪声图像的方差、噪声图像与经滤波图像的相关系数以及二者的和随迭代次数的变化的示意图。如图15(B)所示，与相关系数曲线相比，所得到的和曲线中的谷点被提前了。根据多次实验，如果将图15(B)所示的相关系数曲线的谷点作为停止迭代滤波的时间点，则在经过迭代滤波的输出图像中，有用信息会被模糊；而如果采用所示的和曲线中的谷点，则在经过迭代滤波的输出图像中，原输入图像中的有用信息会被保留，而大部分噪声则被滤除。
另外，上式(9)中的参数c是根据实际应用场景预先确定的值，例如可以是通过对不同类型(或特性)的图像进行分析而获得的实验值或经验值。根据上式(9)可知，c的值越小，则变换后的方差曲线越平缓，反之，c的值越大，则变换后的方差曲线的上升速率越快。c的值不可太小，否则，有可能导致组合单元1303-5得到的组合特征值随迭代次数的变化曲线单调下降，而不出现极小值。c的值也不可太大，否则，会导致组合单元1303-5得到的组合特征值随迭代次数的变化曲线上升速率过快，从而有可能使曲线中的谷点早于适当的迭代滤波停止时间点。因此，在对不同类型(或特性)的图像进行分析时，可以根据不同c值所导致的不同的变换方差曲线，来确定最适合的c值，使得变换方差曲线与相关性曲线的和曲线出现谷点并使谷点的位置对应适当的迭代停止时间点为宜。由于在不同的应用场景下c的值会不同，因此，这里不对c的值作具体的限定。
在实际应用中，可以针对每种类型(特性)的图像来选取预先设置的相应的c值。每种类型(如每种成像设备或扫描协议)的图像可以对应相同的c值，不同类型的图像可以对应不同的c值。作为具体示例，针对多个类型的图像，可以预先确定多个不同的c0值，并将这些值(例如以配置文件的形式)存储在存储单元中(该存储单元可以内置于图像处理设备中，也可以位于图像处理设备外部并可由图像处理设备访问)。这样，在实际的图像处理过程中，无需人工干预，即可对多种类型的图像进行处理，从而大大提高图像处理的自动化水平。
作为另一具体示例，组合单元1303-5可以采用下式对图像特征计算单元1303-4计算得到的方差进行变换：
V i ′ = α × V i - - - ( 9 a )
其中，Vi表示图像特征计算单元1303-4计算得到的噪声图像的方差；表示组合单元1303-5对方差值Vi进行变换后得到的变换值；α表示预定的变换参数。参数α是根据实际应用场景预先确定的值，例如，可以是通过对不同类型(或特性)的图像进行分析而获得的实验值或经验值，在对不同类型(或特性)的图像进行分析时，可以根据不同α值所导致的不同的变换方差曲线，来确定最适合的α值，使得变换方差曲线与相关性曲线的和曲线出现谷点并使谷点的位置对应适当的迭代停止时间点为宜，这里不作限定。在实际应用中，可以针对每种类型(特性)的图像来选取预先设置的相应的α值。每种类型的图像(如每种成像设备生成的图像或根据每种扫描协议生成的图像)可以对应相同的α值，不同类型的图像可以对应不同的α值。
组合单元1303-5将得到的组合特征值输出到停止判断单元1303-3。停止判断单元1303-3用于在所述组合特征值在迭代滤波过程中达到极小值(谷点)时，指示滤波装置1301停止迭代滤波(图14的步骤1406-3)。具体地，停止判断单元1303-3判断组合特征值是否达到极小值；若是，则指示滤波装置1301停止迭代滤波，否则，指示滤波装置1301继续进行下一次迭代。停止判断单元1303-3可以采用任何适当的方法来判断组合特征值是否到达谷点(极小值)，例如，可以将当前滤波得到的组合特征值与之前滤波得到的组合特征值相比较，如果前者大于或等于后者，则可以认为组合特征值呈现增大的趋势，从而可以停止迭代滤波。当然，停止判断单元1303-3还可以采用其他方法来确定组合特征值是否到达谷点，这里不一一列举。
在参考图13-15描述的图像处理设备或方法中，通过利用经滤波图像与噪声图像的相关系数以及噪声图像的预定图像特征值的组合特征值，能够有效地确定停止迭代滤波的时间，在滤除图像中的大部分噪声的同时，能够保证图像中的有用信息不被平滑掉。
本公开的一些实施例还提供了利用非线性扩散滤波方法进行迭代滤波的图像处理设备或方法。
P.Perona和J.Malik的论文“Scale-Space and Edge Detection UsingAnisotropic Diffusion”(刊于Proceedings of IEEE，Computer SocietyWorkshop on Computer Vision，1987年11月，第16-22页)(下文称为相关文献2)提出了非线性扩散滤波方法(应注意，该相关文献2将其提出的方法自称为各向异性扩散(ATD)，然而，该方法的本质是非线性扩散，而不是真正的各向异性扩散。尽管如此，本领域的一些论文和其他文献中仍然将相关文献2的方法称为各向异性扩散)，其中提出了下列等式(称为P-M等式)：
∂ I ∂ t = ▿ ( g * ▿ I ) - - - ( 10 )
在上式中，I表示图像，g表示扩散系数。扩散系数g为图像梯度|▽I|的函数：
g ( | ▿ I | ) = 1 1 + ( | ▿ I | / K ) 2 - - - ( 11 )
式(11)中的K称为梯度阈值。梯度阈值K是用于控制对图像中边缘的敏感性的参数。
图16-17所示的实施例提供了利用图像的梯度分布来确定在非线性扩散滤波中使用的梯度阈值K的图像处理设备和方法，其中，图16为示出根据该实施例的图像处理设备的结构的示意性框图，图17是示出根据该实施例的图像处理方法的示意性流程图。
如图16所示，图像处理设备1600可以包括滤波装置1601、梯度阈值确定装置1605和迭代滤波停止装置1603。图像处理设备1600可以采用如图17所示的图像处理方法进行图像滤波。
滤波装置1601用于基于非线性扩散滤波方法对输入图像进行迭代滤波(图17中的步骤1702)，并在每次迭代后将得到的经滤波图像输出到迭代滤波停止装置1603。滤波装置1601可以采用任何非线性扩散滤波算法进行滤波，这里不作限定。
迭代滤波停止装置1603用于根据经滤波图像在迭代滤波过程中随着迭代次数的增加而变化的速率来确定是否停止迭代滤波(图17中的步骤1706)。具体地，迭代滤波停止装置1603可以具有与上述实施例/示例中描述的迭代滤波停止装置103、303、803、903或1303相似的功能和结构(如采用上文参考图2-15描述的方法来确定迭代滤波停止时间)，这里不再重复。
梯度阈值确定装置1605用于在滤波装置1601进行每次滤波之前计算输入图像或上一次滤波得到的经滤波图像中的每一图像单元的梯度值(图17的步骤1710)，并根据图像中的所有图像单元的梯度分布来确定梯度阈值K(图17中的步骤1712)。在此，“图像单元的梯度值”是指该图像单元与其邻域相比的灰度变化值。
这里所述的图像单元可以是图像中的像素，如二维图像中的像素(Pixel)或三维图像中的像素(Voxel)等等。图18示出了图像中的各个图像单元的梯度值分布的一个示例。在图18中，横轴表示图像单元的梯度值，纵轴表示与梯度值对应的图像单元的个数。作为一个示例，梯度阈值确定装置1605可以获取与输入图像对应的比例值，并按照该比例值、根据梯度分布来计算梯度阈值K，使得图像中梯度值小于梯度阈值K的图像单元的个数与所有图像单元的个数的比接近所述比例值(当然，对于该比例可以做其他数学变化，例如让其反映图像中梯度值小于梯度阈值K的图像单元的个数与梯度值不小于梯度阈值K的图像单元的个数的比)。
作为一个具体示例，可以将图像单元的梯度值按照从小到大或者从大到小的顺序排序，并将排名与所述比例值对应的梯度值确定为梯度阈值，使得梯度值小于梯度阈值K的图像单元的个数与所有图像单元的个数的比接近所述比例值。作为另一具体示例，还可以计算图像中图像单元的梯度的直方图，并根据直方图来确定梯度阈值，使得梯度值小于梯度阈值K的图像单元的个数与所有图像单元的个数的比接近所述比例值。
这里所述的比例值可以是根据输入图像的类型或特性(如医学图像的成像设备或扫描协议等)而预先确定的值，可以用P来表示，且0＜P＜1。每种类型的图像可以对应相同的比例值，而不同的图像类型可以对应不同的比例值。应理解，在实际应用场景中，针对每种图像类型，可以通过分析多个该类型的图像来预先确定相应的比例值，从而在后续的非线性扩散滤波中使用，这里不作详述。例如，可以将预先确定的、与多个图像类型对应的多个比例值(例如以配置文件的形式)存储在存储单元中(该存储单元可以内置于图像处理设备中，也可以位于图像处理设备外部并可由图像处理设备访问)，在实际的非线性扩散滤波中，梯度阈值确定装置1605可以根据输入图像的类型或特性(如成像设备或扫描协议的类型等)来获取所存储的、与输入图像的类型或特性对应的比例值，而无需人工手动设置该参数，因而可以进一步提高图像处理的自动化水平。由于不同应用场景下该比例值的大小是不同的，因此，这里也不具体限定该比例值的数值。
梯度阈值确定装置1605可以采用任何适当的方法来计算图像中的图像单元的梯度值，并根据图像中的所有图像单元的梯度分布来确定梯度阈值。作为一个示例，假使图像为n维图像(n≥1)，可以基于下式来计算图像单元的梯度值：
| I grad ( D 1 , D 2 , . . . , D n ) | = I D 1 2 + I D 2 2 + . . . + I D n 2 - - - ( 12 )
其中，D1，D2，...，Dn表示图像的n个维度，|Igrad(D1，D2，...，Dn)|表示图像中的位于(D1，D2，...，Dn)的像素的梯度值，而(1≤k≤n)表示该像素在第Dk维上的导数。
在每次迭代中，梯度阈值确定装置1605将计算得到的梯度阈值输出到滤波装置1601，滤波装置1601根据梯度阈值确定装置1605确定的梯度阈值来进行非线性扩散滤波(图17中的步骤1702)。
根据本公开的实施例的图像处理方法和设备可以应用于各种图像的滤波处理。这些图像可以是一维或者多维的(如二维或三维等)。例如，根据本公开的实施例的图像处理方法和设备可应用于医学图像(如根据通过医疗诊断装置获得的数据而形成的图像)的降噪处理。作为一个示例，上述方法的各个步骤以及上述设备的各个组成模块和/或单元可以实施为医疗诊断装置(如X射线诊断装置、UL诊断装置、CT装置、MRI诊断装置或PET装置等)中的软件、固件、硬件或其组合，并作为该医疗诊断装置中的一部分。作为一个示例，可以在已有医疗诊断装置中实施根据本公开上述方法和/或设备，其中对已有医疗诊断装置的各组成部分作一定修改即可。作为另一示例，上述方法的各个步骤以及上述设备的各个组成模块和/或单元可以实施为独立于所述医疗诊断装置的装置。上述设备中各个组成模块、单元通过软件、固件、硬件或其组合的方式进行配置时可使用的具体手段或方式为本领域技术人员所熟知，在此不再赘述。
作为一个示例，上述方法的各个步骤以及上述设备的各个组成模块和/或单元可以实施为软件、固件、硬件或其组合。在通过软件或固件实现的情况下，可以从存储介质或网络向具有专用硬件结构的计算机(例如图19所示的通用计算机1900)安装构成用于实施上述方法的软件的程序，该计算机在安装有各种程序时，能够执行各种功能等。
在图19中，运算处理单元(即CPU)1901根据只读存储器(ROM)1902中存储的程序或从存储部分1908加载到随机存取存储器(RAM)1903的程序执行各种处理。在RAM 1903中，也根据需要存储当CPU 1901执行各种处理等等时所需的数据。CPU 1901、ROM 1902和RAM 1903经由总线1904彼此链路。输入/输出接口1905也链路到总线1904。
下述部件链路到输入/输出接口1905：输入部分1906(包括键盘、鼠标等等)、输出部分1907(包括显示器，比如阴极射线管(CRT)、液晶显示器(LCD)等，和扬声器等)、存储部分1908(包括硬盘等)、通信部分1909(包括网络接口卡比如LAN卡、调制解调器等)。通信部分1909经由网络比如因特网执行通信处理。根据需要，驱动器1910也可链路到输入/输出接口1905。可拆卸介质1911比如磁盘、光盘、磁光盘、半导体存储器等等根据需要被安装在驱动器1910上，使得从中读出的计算机程序根据需要被安装到存储部分1908中。
在通过软件实现上述系列处理的情况下，从网络比如因特网或存储介质比如可拆卸介质1911安装构成软件的程序。
本领域的技术人员应当理解，这种存储介质不局限于图19所示的其中存储有程序、与设备相分离地分发以向用户提供程序的可拆卸介质1911。可拆卸介质1911的例子包含磁盘(包含软盘(注册商标))、光盘(包含光盘只读存储器(CD-ROM)和数字通用盘(DVD))、磁光盘(包含迷你盘(MD)(注册商标))和半导体存储器。或者，存储介质可以是ROM 1902、存储部分1908中包含的硬盘等等，其中存有程序，并且与包含它们的设备一起被分发给用户。
本公开还提出一种存储有机器可读取的指令代码的程序产品。所述指令代码由机器读取并执行时，可执行上述根据本公开实施例的方法。
相应地，用于承载上述存储有机器可读取的指令代码的程序产品的存储介质也包括在本公开的公开中。所述存储介质包括但不限于软盘、光盘、磁光盘、存储卡、存储棒等等。
在上面对本公开具体实施例的描述中，针对一种实施方式描述和/或示出的特征可以用相同或类似的方式在一个或更多个其它实施方式中使用，与其它实施方式中的特征相组合，或替代其它实施方式中的特征。
应该强调，术语“包括/包含”在本文使用时指特征、要素、步骤或组件的存在，但并不排除一个或更多个其它特征、要素、步骤或组件的存在或附加。
在上述实施例和示例中，采用了数字组成的附图标记来表示各个步骤和/或单元。本领域的普通技术人员应理解，这些附图标记只是为了便于叙述和绘图，而并非表示其顺序或任何其他限定。
此外，本公开的方法不限于按照说明书中描述的时间顺序来执行，也可以按照其他的时间顺序地、并行地或独立地执行。因此，本说明书中描述的方法的执行顺序不对本公开的技术范围构成限制。
尽管上面已经通过对本公开的具体实施例的描述对本公开进行了披露，但是，应该理解，上述的所有实施例和示例均是示例性的，而非限制性的。本领域的技术人员可在所附权利要求的精神和范围内设计对本公开的各种修改、改进或者等同物。这些修改、改进或者等同物也应当被认为包括在本公开的保护范围内。

Claims (12)

1.一种图像处理设备，包括：
滤波装置，用于对输入图像进行迭代滤波；以及
迭代滤波停止装置，用于所述迭代滤波过程中每次滤波后根据得到的经滤波图像来确定是否停止所述迭代滤波，
所述迭代滤波停止装置还包括：
噪声图像获取单元，用于在所述滤波装置每次进行滤波后根据所得到的经滤波图像和所述输入图像来估算噪声图像；
相关性获取单元，用于估算所述噪声图像与所述经滤波图像的相关系数；
图像特征计算单元，用于计算所述噪声图像的方差或者标准差；
组合单元，用于组合所述相关系数与所述噪声图像的方差或者标准差，得到组合特征值；及
停止判断单元，被配置为：当所述组合特征值在所述迭代滤波过程中达到极小值时停止所述迭代滤波。
2.如权利要求1所述的图像处理设备，其中，所述停止判断单元根据所述相关系数在所述迭代滤波过程中下降的速率来确定是否停止所述迭代滤波。
3.如权利要求2所述的图像处理设备，其中，所述停止判断单元被配置为：
在停止所述迭代滤波的情况下，计算表示所述相关系数在所述迭代滤波过程中下降的速率的特征值，并判断所述特征值与预定阈值之间是否满足预定的关系。
4.如权利要求3所述的图像处理设备，其中，所述特征值为表示所述相关系数在所述迭代滤波过程中随迭代次数的增加而下降的曲线的倾斜度的值。
5.如权利要求3或4所述的图像处理设备，其中，所述预定阈值为根据输入图像的特性而预先确定的值。
6.如权利要求1所述的图像处理设备，其中，所述图像特征计算单元还被配置用于：对计算得到的所述噪声图像的方差或者标准差进行变换，得到经变换的方差或者标准差，并且
其中，所述组合单元被配置为：计算所述相关系数与经变换的方差或者标准差的和。
7.如权利要求1到4中任一项所述的图像处理设备，其中，所述滤波装置被配置为基于非线性扩散方法进行滤波。
8.如权利要求7所述的图像处理设备，还包括：
梯度阈值确定装置，用于在滤波装置进行每次滤波之前计算所述输入图像或上一次滤波得到的经滤波图像中的每一图像单元的梯度值，并根据所述图像中的所有图像单元的梯度分布来确定梯度阈值，
其中，所述滤波装置被配置为利用所述梯度阈值确定装置确定的梯度阈值来基于非线性扩散方法进行滤波。
9.如权利要求8所述的图像处理设备，其中，所述梯度阈值确定装置被配置为：
获取与所述输入图像对应酶比例值，并按照所述比例值、根据所述梯度分布来计算所述梯度阈值，其中，图像中梯度值小于所述梯度阈值的图像单元的个数与所有图像单元的个数的比接近所述比例值。
10.如权利要求9所述的图像处理设备，其中，所述比例值为根据所述输入图像的特性而预先确定的值。
11.如权利要求1到4中任一项所述的图像处理设备，其中，所述输入图像为根据通过医疗诊断装置获得的数据而形成的医学图像。
12.一种图像处理方法，包括：
对输入图像进行迭代滤波；以及
根据所述迭代滤波过程中每次滤波后得到的经滤波图像来确定是否停止所述迭代滤波，
在确定是否停止所述迭代滤波时，进而，
每次滤波处理后，根据得到的经滤波图像和所述输入图像来估算噪声图像；
计算所述噪声图像与所述经滤波图像的相关系数；
计算噪声图像的方差或者标准差；
组合该相关系数与该噪声图像的方差或者标准差，得到组合特征值，当所述组合特征值在所述迭代滤波过程中达到极小值时停止所述迭代滤波。
CN201110045243.3A 2011-02-22 2011-02-22 图像处理设备和方法 Active CN102646265B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201110045243.3A CN102646265B (zh) 2011-02-22 2011-02-22 图像处理设备和方法
JP2012034225A JP5925515B2 (ja) 2011-02-22 2012-02-20 画像処理装置、画像処理方法、及び画像処理プログラム
US13/402,435 US9202261B2 (en) 2011-02-22 2012-02-22 Image processing apparatus and method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201110045243.3A CN102646265B (zh) 2011-02-22 2011-02-22 图像处理设备和方法

Publications (2)

Publication Number Publication Date
CN102646265A CN102646265A (zh) 2012-08-22
CN102646265B true CN102646265B (zh) 2015-09-23

Family

ID=46652787

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201110045243.3A Active CN102646265B (zh) 2011-02-22 2011-02-22 图像处理设备和方法

Country Status (3)

Country Link
US (1) US9202261B2 (zh)
JP (1) JP5925515B2 (zh)
CN (1) CN102646265B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103728652B (zh) * 2013-12-18 2015-12-30 中国原子能科学研究院 一种核元素γ能谱小波降噪的方法、装置及核元素探测器
WO2016011489A1 (en) * 2014-07-23 2016-01-28 The University Of Sydney Thoracic imaging for cone beam computed tomography
WO2017132600A1 (en) * 2016-01-29 2017-08-03 Intuitive Surgical Operations, Inc. Light level adaptive filter and method
CN109171706B (zh) * 2018-09-30 2020-12-29 南京信息工程大学 基于分类匹配和分数阶扩散的心电信号去噪方法和系统
JP7246194B2 (ja) * 2019-01-28 2023-03-27 キヤノンメディカルシステムズ株式会社 医用画像処理装置、医用画像処理方法、およびプログラム
JP7391824B2 (ja) 2020-12-03 2023-12-05 株式会社東芝 情報処理装置、情報処理方法およびプログラム
CN116157823A (zh) * 2021-09-23 2023-05-23 京东方科技集团股份有限公司 一种显示设备系统及自适应增强画质的方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101916433A (zh) * 2010-08-10 2010-12-15 西安电子科技大学 基于偏微分方程的强噪声污染图像的去噪方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6687388B2 (en) * 2000-01-28 2004-02-03 Sony Corporation Picture processing apparatus
GC0000235A (en) 2000-08-09 2006-03-29 Shell Int Research Processing an image
US7004904B2 (en) * 2002-08-02 2006-02-28 Diagnostic Ultrasound Corporation Image enhancement and segmentation of structures in 3D ultrasound images for volume measurements
US20070031058A1 (en) * 2005-06-08 2007-02-08 Canamet Canadian National Medical Technologies Inc. Method and system for blind reconstruction of multi-frame image data
JP4653059B2 (ja) * 2006-11-10 2011-03-16 オリンパス株式会社 撮像システム、画像処理プログラム
JP5164645B2 (ja) * 2008-04-07 2013-03-21 アルパイン株式会社 カルマンフィルタ処理における繰り返し演算制御方法及び装置

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101916433A (zh) * 2010-08-10 2010-12-15 西安电子科技大学 基于偏微分方程的强噪声污染图像的去噪方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Selection of Optimal Stopping Time for Nonlinear Diffusion Filtering;Pavel Mrazek等;《International Journal of Computer Vision》;20030531;第52卷(第2-3期);第189-203页 *
基于噪声方差确定非线性扩散除噪声的最优停止时间;刘鹏等;《电子与信息学报》;20090930;第31卷(第9期);第2084-2087页 *
数字图像中椒盐噪声的滤波算法研究;郭红伟;《万方学位论文全文数据库》;20101222;第40-46页 *

Also Published As

Publication number Publication date
JP5925515B2 (ja) 2016-05-25
US9202261B2 (en) 2015-12-01
JP2012174275A (ja) 2012-09-10
CN102646265A (zh) 2012-08-22
US20120213451A1 (en) 2012-08-23

Similar Documents

Publication Publication Date Title
CN102646265B (zh) 图像处理设备和方法
Yang et al. Local statistics and non-local mean filter for speckle noise reduction in medical ultrasound image
CN102999884B (zh) 图像处理设备和方法
US9179881B2 (en) Physics based image processing and evaluation process of perfusion images from radiology imaging
CN101933045B (zh) 短轴后期增强心脏mri的自动3d分割
US8346011B2 (en) Reducing noise in an image
EP2380132B1 (en) Denoising medical images
Ukwatta et al. Three‐dimensional ultrasound of carotid atherosclerosis: semiautomated segmentation using a level set‐based method
Saha Tensor scale: A local morphometric parameter with applications to computer vision and image processing
US20130004085A1 (en) Image processing apparatus and method
US9299146B2 (en) Image processing device, image processing method and image processing program
US20130202177A1 (en) Non-linear resolution reduction for medical imagery
US9536318B2 (en) Image processing device and method for detecting line structures in an image data set
US20140355863A1 (en) Image processing apparatus, method and medical image device
Wen et al. A novel Bayesian-based nonlocal reconstruction method for freehand 3D ultrasound imaging
Lorza et al. Carotid artery lumen segmentation in 3D free-hand ultrasound images using surface graph cuts
Wang et al. Anisotropic diffusion filtering method with weighted directional structure tensor
US7711164B2 (en) System and method for automatic segmentation of vessels in breast MR sequences
US10417764B2 (en) System and methods for diagnostic image analysis and image quality assessment
US9196030B2 (en) System and method for determining a property of blur in a blurred image
US8157736B2 (en) System and method for feature detection in ultrasound images
Aviles et al. Robust cardiac motion estimation using ultrafast ultrasound data: a low-rank topology-preserving approach
JP5132559B2 (ja) デジタル画像のセグメント化方法およびコンピュータ読み取り可能なプログラム記憶装置
Neubert et al. Constrained reverse diffusion for thick slice interpolation of 3D volumetric MRI images
WO2013078451A1 (en) Methods of determining local spectrum at a pixel using a rotationally invariant s-transform (rist)

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
C41 Transfer of patent application or patent right or utility model
TR01 Transfer of patent right

Effective date of registration: 20160802

Address after: Japan Tochigi

Patentee after: Toshiba Medical System Co., Ltd.

Address before: Tokyo, Japan, Japan

Patentee before: Toshiba Corp

Patentee before: Toshiba Medical System Co., Ltd.