图像噪声的降低的制作方法
专利名称：图像噪声的降低的制作方法
图像噪声的降低相关申请本申请要求根据美国专利法35USC 119(e)的申请日为2007年12月25日的美国 临时专利申请61/016，578的权益。上述专利申请文件的所有内容都以其整体引入在此作为参考。发明的技术领域与
背景技术：
本发明，在它的一些实施例中，涉及一种降低图像噪声的方法，尤其是，但不完全 是，涉及一种采用非线性过滤器降低医学图像噪声的方法。下列出版物和发明涉及图像噪声降低方法、图像采集和/或计算机视觉。US 2007053477——锥形束CT成像与扇束CT成像的全球去噪方法与装置KR20050031210——图像去噪的方法与装置JP2000050109——用于去除噪声的非线性图像过滤器US6459755——用于管理低剂量CT扫描的方法与装置
US2003099405——降低CT剂量过滤器，带有计算有效实施EP1774837——降低活性剂量的装置与方法JP200139874——用于MRI的磁场发生器W02007047599——用于高增益磁共振的方法与装置G01V3/00, G01R33/34Optimal Mass Transport for Registration and Warping, International Journal of Computer Vision, Volume 60, Issue 3(December 2004), Pages 225-240, Steven Haker, Lei Zhu, Allen Tannenbaum, Sigurd Angenent A Metric for Distributions with Applications to Image Databases, ICIP 1998, Pages 59-66, Rubner Yossi, Tomasi Carlo, Guibas, J. Leonidas. Shape Matching and Object Recognition Using Shape Contexts, IEEE T-PAMI, Volume 24, No.4, (April 2002), Belongie Serge, Jitendra Malik, Puzicha Jan. Matching 3D Models with Shape Distributions, Proceedings of the International Conference on Shape Modeling & Applications 2001, Pages 154—166, Robert Osada, Thomas Funkhouser, Bernard Chazelle, and David Dobkin P. J. Burt, E. H. Adelson, " The Laplacian Pyramid as a Compact Image Code, " IEEE Trans, on Communications, pp.532__540, April 1983Iddo Drori,Daniel Cohen-Or,Hezy Yeshurun, ACM Transactions on Graphics 22(3)，(Proc. of SIGGRAPH 2003)，303-312.John Goutsias and Henk J. A. M. Hei jmans, "Nonlinear Multiresolution Signal Decomposition Schemes-Part I :Morphological Pyramids，，，IEEE Trans, on Image Processing, Vol. 9, No. 11,November 2000.John Goutsias and Henk J.A. M. Heijmans, "Nonlinear Multiresolution Signal Decomposition Schemes-Part II :Morphological Wavelets", IEEE Trans, onlmage Processing, Vol. 9, No. 11, November 2000.
Jean Serra, Image Analysis and Mathematical Morphology, 1982.
发明内容
本发明的一些实施例在一个方面涉及一种降低图像中的像素噪声的方法，通过估 计该像素的没有噪声的真实灰度值，基于其他类似像素的灰度值，识别为它们相邻像素的 类似特征。这些特征包括取决于在邻域的像素的灰度的分布函数的特征，或者相对于彼此 旋转或缩放的比较邻域的特征，或者参与应用变换或过滤到邻域的特征，以致优先选择中 间尺度范围的结构，或者仅是某些方向的特征。根据本发明的一个示例性实施例，提供了一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻 域的对应特征的值；c)基于图像的检索像素的原始灰度值或者变换灰度值，计算对于正被检验的每个 像素的降低的噪声灰度值，相对于检索像素的一个或多个特征值，正被检验的像素具有更 大敏感性的类似特征值；其中，所述计算至少一个特征值的步骤包括计算在邻域的像素的原始灰度值或 变换灰度值的特征性分配，除了计算在该邻域的所有像素的平均灰度值之外。在本发明的一个示例性实施例中，计算特征值的步骤包括计算第二时刻分布或 更高时刻分布。可选地，在邻域的像素的灰度值是对于那些像素的图像的原始灰度值。在本发明的一个示例性实施例中，在邻域的像素的灰度值是由过滤器变换的图像 的灰度值。可选地，所述过滤器是高斯过滤器。在本发明的一个示例性实施例中，计算特征值的步骤还包括计算所述像素在第 二邻域的灰度值的至少一个特征性分布，该第二邻域是所述邻域的适当的子集。根据本发明的一个示例性实施例，提供了一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻 域的对应特征的值；以及c)计算对于正被检验的每个像素的降低的噪声灰度值，基于它的检索像素的原始 灰度值或者变换灰度值，检索像素的一个或多个特征值比正被检验的像素的类似特征值具 有更大的敏感性；其中，所述计算至少一个特征的值的步骤包括评估在邻域的图像的变换，或者评 估在邻域的图像上的过滤器的效果，以致优先选择了在邻域的尺寸和一些像素尺寸之间的 中间尺度范围的结构，或者一些方向而不是其他方向的结构，或者兼有这两种结构。可选 地，计算特征值的步骤包括评估线形过滤器对于该图像的响应。可选地，所述线性过滤器 是一个小波滤波器。可选地，所述线性过滤器是一个方向过滤器。可选地，计算特征值的步 骤包括评估在超过一个方向上的过滤器的效果。在本发明的一个示例性实施例中，计算特征值的步骤包括评估在所述邻域的图 像的非线性形态变换。可选地，所述的形态变换是多尺度变换。可选地，所述的形态变换是 形态算子。
在本发明的一个示例性实施例中，在所述邻域的方向的角上所述特征的值的相关 性具有对称性，以致如果所述邻域是旋转至少一个角度，对于在该邻域的图像像素的任意 灰度值，所述特征值是不变的。在本发明的一个示例性实施例中，对于至少一个特征，正被检验的像素的特征值 与检索像素的对应特征值是以它们在不同相对方向或尺度或两者兼有的邻域来评估的。在本发明的一个示例性实施例中，所述减少的噪声灰度值是基于所述检索像素的 加权平均灰度值。在本发明的一个示例性实施例中，所述加权平均值采用的权重取决于在检索像素 的特征值与正被检验的像素的对应的特征值之间的差的测算。在本发明的一个示例性实施例中，所述的方法还包括a)识别残留图像的相对更多结构部分，表示从原始图像到减噪图像的变化；以及b)存储所述结构部分到减噪图像，以产生一个改善的减噪图像。在本发明的一个示例性实施例中，所述的方法还包括估计一个或多个检索像素 的噪声水平，其中，所降低的噪声灰度值是低相关性或不相关于具有更高的估计噪声水平 的检索像素。在本发明的一个示例性实施例中，所述的方法还包括a)在图像中分配像素的一个补丁到在这些补丁的向量空间内的复数个簇之一，该 补丁表示为灰度的向量；b)近似所述补丁作为所降低的一组非相关参数的函数，参数量少于在所述补丁内 的像素量，与所述簇关联；以及c)根据在所述补丁的实际灰度与作为所降低的一组非相关参数的函数的所述补 丁的近似之间的差，校正所述补丁的像素的灰度。可选地，所述函数是所降低的一组基向量 的线性组合，非相关参数是在该线性组合内的基向量的系数。根据本发明的一个示例性的实施例，提供了一种获得带有改变的图像采集参数的 医学图像的方法，包括a)确定一组一个或多个图像采集参数，这些参数会在图像降噪后达到想要的噪声 水平；b)获得带有所述采集参数的图像；以及c)根据前述任意权利要求所述的方法降低该图像的噪声。可选地，所述获得图像的步骤包括获得CT图像；以及所述确定一组一个或多个 图像采集参数的步骤包括确定一组与该图像未降低噪声时的情况相比具有降低的X-射 线剂量的参数。可选地，所述获得图像的步骤包括获得MRI图像；以及所述确定一组一个 或多个图像采集参数的步骤包括确定一组与该图像未降低噪声时相比具有降低的静磁场 或减少的采集时间的参数。除非另作定义，这里所采用的素有技术术语和/或科学术语都具有本发明所涉及 领域的技术人员所熟知的普通含义。与这里所描述的方法和材料的类似或等同的方法和材 料都可被用于本发明的实施例的实施或测试，下面给出了示例性的方法和/或材料。在发 生冲突时，本专利说明书包括这些定义，都会控制。另外，所有材料、方法和例子都仅是示例 性的，而不试图作为对本发明的限制。
6
本发明的实施方式的方法和/或系统的执行，可涉及手动、自动或两者结合地执 行或完成所选择的任务。而且，根据本发明所述的方法和/或系统的实际设备和装置的实 施例，可通过硬件、软件或固件或它们的结合，采用操作系统来执行几种选择的任务。例如，用于执行根据本发明所述的实施例的选择任务的硬件可以是以芯片或电路 的方式实施。对于软件，本发明所述的实施例的选择任务可采用任意合适的操作系统通过 计算机来执行多个软件指令来实施。在本发明的一个示例性实施例中，根据这里所描述的 方法和/或系统的示例的一个或多个任务是通过数据处理器来执行，例如，用于执行多个 指令的计算平台。可选地，该数据处理器包括用于存储指令和/或数据的易失性存储器，和 /或用于存储指令和/或数据的非易失性存储器，例如，磁性硬盘和/或可移动介质。可选 地，还提供了网络连接。显示器和/或用户输入设备例如键盘或鼠标也可选地提供。附图
简要说明在这里仅通过实施例并附图，来描述本发明的一些实施方式。现在特别详细参考 有关附图，需要强调的是，通过实施例的方式来显示特别内容，用于说明本发明的实施例的 示例性讨论的目的。在这方面，本说明书采用附图使本领域技术人员能清楚了解本发明的 实施方式。在附图中图IA示意性地显示了现有技术的没有噪声的两维图像，而图IB示意性地显示了 现有技术的带有噪声以及带有选择像素和邻域的相同图像；图2A示意性地显示了图IB的图像以及选择的像素和邻域，以及显示了现有技术 中的其他类似于所选择的像素的其他像素；图2B示意性地显示了图IB的图像以及选择的像素和邻域，以及显示了根据本发 明的一个示例性实施例的类似于所选择的像素的其他像素；图3是根据本发明的一个示例性实施例所述的一种降低图像噪声的方法的流程 图；图4A是采用相对低的X-射线剂量获得的噪声CT图像；图4B显示了采用图3的方法降低噪声的图4A的图像；图4C是类似于图4A中图像的一个低噪声CT图像，但该图是采用相对高的X-射 线剂量获得的；图5是根据本发明的一个示例性实施例所述的一种方法的流程图，该方法采用遗 传算法找到更有效的特征组加权到图3所采用的方法；图6是根据本发明的一个示例性实施例所述的一种获得图像的方法的流程图，该 方法采用改变的采集参数以补偿噪声降低，因为该图像可被变换为采用图3所示的降噪方 法的具有正常噪声水平的图像；图7是根据本发明的一个示例性实施例所述的在任意结构降噪后检查残留图像 的方法的流程图，该结构是从原始图像不正确地移除的，并恢复高结构化的部分以降低图
像噪声；图8是根据本发明的一个示例性实施例所述的一个方法的流程图，该方法找到一 个相关的矩阵，采用该矩阵找到在图3所示方法中的距离测度；以及图9是根据本发明的一个示例性实施例所述的一个方法的流程图，该方法估计一个图像的特定体素的噪声。优选实施方式详述本发明在一些实施例中涉及一种降低图像噪声的方法，尤其是，但不完全是涉及 一种采用非线性过滤器的降低图像噪声的方法。本发明的一些实施例的一个方面涉及一种降低图像噪声的方法，在该方法中，比 较不同像素的邻域，以便找到具有类似邻域的像素。采用一个或多个特征来决定这个类 似性，包括参与评估在邻域中原始灰度或变换灰度的特征性分布函数的特征的至少一个特 征。可选地，该特征性分布函数是矩分布函数，高于第一矩(平均灰度值)，例如第二矩(标 准偏差)、第三矩(偏斜度)或者更高矩的分布函数。可选地，所述分布函数是该图像的原 始灰度值的分布函数。此外，所述分布函数也可是变换后图像(例如，由高斯过滤器过滤后 的图像)的灰度的分布函数。可选地，可找到该分布函数的特征，既用于在邻域的所有像 素，也仅用于位于第二邻域的像素，该第二邻域是所述邻域的一部分。本发明的一些实施例的一个方面涉及一种降低图像噪声的方法，在该方法中，比 较不同像素，采用以下特征之一参与对该图像应用一个变换的特征，或者参与评估在图像 上过滤的效果的特征，以致在邻域尺寸和一些像素尺寸之间的中间尺度范围的结构，或者 优先选择一些方向而不是其他方向的结构。可选地，该过滤器是线性过滤器，例如小波滤波 器。可选地，该过滤器是方向性的。可选地，该变换是非线性形态变换，例如，多尺度变换， 或者形态算子。在本发明的一些实施例中，评估参与评估方向性过滤器在超过一个方向上 的反应的特征。可选地，该特征具有与角度相关的对称度，以致如果所述邻域是旋转至少一 个角度，例如180度或90度，不论该图像看起来像什么，所述特征值是不变的。本发明的一些实施例的一个方面涉及一种降低图像噪声的方法，在该方法中，比 较不同像素的邻域，采用至少一个特征值，该特征值是对于这些像素之一计算的，该像素的 邻域相对于其他像素的邻域旋转或缩放。在本发明的一些实施例中，邻域是定义用于像素群，而不是用于个别像素，而噪声 降低是作用在像素群，而不是作用在个别像素。此外，邻域可被定义在位于像素之间的体 素，而噪声降低是作用在该体素，这些体素的灰度值是由接近它们的像素的插值来定义的。本发明的一些实施例的一个方面涉及一种降低图像噪声的方法，在该方法中，在 图像的噪声降低后，检查残留的图像(在降噪前后的图像的差别)，看它是否具有结构元 件，如果有的话，结构元件被存储在该图像中。在本发明的一些实施例中，对于不同像素估计噪声水平，以及当像素无噪声的真 实灰度值是采用类似像素的灰度值来估计时，越多噪声的像素对所估计的灰度值的影响越 少。在本发明的一些实施例的一个方面涉及一种噪声降低系统，在该系统中，像素的 噪声分量是通过对一组测试图像的检查(第一检查)来估计的，将在这些图像上找到的像 素的补丁分类为在这些补丁的理论上的向量空间的簇，将这些补丁作为向量考虑，它们的 分量是它们的像素的灰度值。对于每个簇，发现降噪的一组独立参数(数量少于在该补丁 中的像素)，例如，基向量的一组降噪的线性系数，这提供了大多数补丁在该簇内的好的近 似，对于测试图像。当对一个噪声图像进行降噪处理时，围绕每个像素的补丁被分类为属于 其中一个簇，而它的真实值，无噪声，是采用对于该簇的降噪参数组来估计的。
在本发明的一些实施例的一个方面涉及选择用于采集图像的改变的参数，例如， 较低X-射线剂量的CT图像，或者较低静磁场或较短采集时间的MRI图像，利用降噪方法 来最终获得足够质量的图像以满足诊断需求，例如，采用更高的X-射线剂量或更高的静磁 场，没有噪声降低。举例来说，这里所描述的方法可被应用在图像采集设备或它的工作站(例如，CT 机、MRI机)上，应用在图像处理站和/或通过网络连接到远程位置或远程服务器。为了更好地理解本发明的一些实施例，正如在图2B至图9所示，参考传统(也就 是，现有技术)的图像降低噪声方法如在图IA至图2A中的操作。图IA显示了二维图像100，它包括一个像素数组，每组具有一个数值，该数值是在 黑与白之间的灰度值。在CT图像中的灰度值表示成像物体的实际密度的方便映射，通常表 示为豪森菲尔德单位(HU)。例如，在大脑的CT图像中，该图像通常显现0HU，这表示水的密 度，被映射为黑色，显现70HU，被映射为白色。总的来说，在图像处理文献中，术语“像素(pixel) ”是用于表示二维图像的一个元 素，而“体素(voxel)”是用于表示二维图像的一个元素。因为这里所描述的方法总体上既 可用于二维图像，也可用于三位图像，这里所采用的术语“像素”和“体素”不应理解为将本 说明书限制为二维图像或三维图像的例子。除非另有特别说明，这里所采用的术语“像素” 和“体素”应被理解为应用于任意例子的通用术语，它们通常可交替使用。这里所使用的术语“灰度值(grey value) ”不仅是指黑白图像的亮度，还是指彩 色图像中任意色彩变量的程度，例如，在彩色图像中的红、绿或蓝，或者彩色图像的亮度或 饱和度。在诸如CT或MRI图像的医学图像中，通常只有单一密度变量，例如1\或1~2加权密 度，它被映射到灰度图像的亮度，而在本例中，“灰度值”是特别倾向于，但这里所述的方法 并不限制于它们对医学图像或黑白图像的适用性。这里所描述的降低噪声方法法可被特别 用于医学图像，因为医学图像通常具有相对高的噪声水平，由于在噪声水平和图像采集参 数(例如，X-射线剂量或MRI采集时间)之间通常有折中，施加经济或安全惩罚以用于降 低噪声。而且，因为医学图像总体上没有在“光照(lighting)”差别，一个像素的邻域的特 征通常是它的真实灰度值的号的指示。这是特别真实的，因为医学图像倾向于具有类似结 构，在该图像的不同部分重复出现，有时带有改变的尺度或方向。图像100包括亮区102和暗区104，两区之间具有相当尖锐的边界。在图IB中， 图像108是已加入噪声的图像100。在现有技术中，噪声有时是由像素的灰度值与相邻像素 的灰度值的平均来降低的，处于最接近位置的像素给该像素最大的权重。这在均勻区域做 的很好，没有细节，例如在图像108中的亮区102和暗区104，但会导致在这两区之间的边 界的污点。另一个现有技术的降低噪声方法，双边过滤器，它是非线性过滤器，试图避免由 于像素i的灰度值Ii与其他像素j的灰度值Ij的平均导致的问题。例如，当操作位于(Xi， Yi)的特定像素i时，对位于(Xj，Yj)的另一个体素j的灰度值给予权重Wj,由下式给出 这里，dp是欧几里得距离，在空间中两个像素之间，而Ili-IjI可被认为是在这两 个像素之间的理论上“距离”，它们互相类似的程度的测度。对于像素i的新灰度值是由下式定义的
这里N是围绕像素i的搜索窗口，而总和是对于在该搜索窗口
得所有像素j。用于降低噪声的另一类型的非线性过滤器是在以下文献中定义的L. Rudin, S. Osher, and Ε. Fatemi, " Nonlinear total variation based noise removal algorithms, “ Physica D60，259-268 (1992)。在非本地过滤器装置中，现有技术的进一步发展，两个像素的类同之处取决于这 两个像素的邻域的像素与像素的比较。例如，为降低像素i (在图IB中标记为110)的噪声 水平，邻域Mi (在图IB中标记为112)，是围绕像素110来限定的。接着搜索其他像素j，其 具有围绕每个检索像素j的相同尺寸和形状的邻域Mj,在邻域112的像素与每个检索像素 j的邻域的相应像素之间找到均方误差MSE(Mi,Mj)。对于检索像素的均方误差是小的，在它 们的邻域与像素110的邻域之间，给出最大权重Wj，当平均检索像素的灰度值时，获得对于 像素110的降噪灰度值。该权重Wj由下式给出 对于像素i的新值是由下式确定 二卞J 7
ZjjeN J图2A显示了图像200，类似于在图IB中的图像108，具有一组像素202，它们具有 类似于邻域112或像素110的邻域。每个像素202具有类似的邻域，因为像素202都是离 开在亮区102和暗区104之间的边界同样的距离，方向都是在几乎相同的方向。在其他现有技术的降低噪声方法中，采用非线性过滤器，两个邻域的类似度是 基于在该邻域的所有像素的平均灰度值，或者在该邻域的像素的灰度值梯度的方向上， 如下文所述Mahmoudi, M. and Sapiro, G. , “ Fast image and video denoising via nonlocal means of similar neighborhoods, " IEEE, Signal Proc. , Vol. 12, no.12, pp. 839-842, Dec. 2005。在由 A. Heidarzadeh and A. N. Avanaki, "An Enhanced Nonlocal Means Algorithm for Image Denoising, "9th ISSPA，Feb. 2007所描述的不同方法中，两个 邻域的相似度取决于这两个邻域的二进制边缘图的均方误差，采用Carmy边缘检测器来确 定，还取决于在这两个邻域的原始图像的均方误差。在详细解释本发明的至少一个实施例之前，需要明确的是，本发明并不限制于以 下说明所提出的应用细节。本发明能采用其他实施例或者以不同的方式来实施或执行。图2B显示了图像204，类似于图像108。根据本发明的一个实施例所述，对于在两 个邻域之间的相似度，采用不同标准计算权重Wj，找到较好的一组检索像素206，它们具有 足够接近的类似于像素110的邻域112的邻域。在图2B所示的特殊例子中，标准(将在下 面详细描述)不依赖于这些邻域的相对方向，以致所有与暗区104有相同距离的像素(如 像素110)都由接近的类似邻域112的邻域，根据这些标准。与采用非本地装置方法带有高 权重的检索像素202相比，以高权重放大的检索像素206组可允许进一步降低噪声，因为有 更多像素来平均灰度值。在本发明的一些实施例中，对于在两个邻域之间的相似度的标准 可取决于这两个邻域的相对方向，而在本发明的这些或其他实施例中，具有高权重的检索像素的质量可以不必大于现有技术方法，但这些检索像素的质量可以更好，在这个意义上， 它们提供了比像素110的真实灰度值更好的估计。图3显示了根据本发明一个实施例所述的一种降低图像噪声的方法的流程图 300。流程图300所示方法是图2B所示方法的概括，具有用于在邻域之间的相似性的不同 标准。在步骤302，获得有噪声的图像。噪声降低算法在一次检查一个像素，初始在步骤 304中设定像素i等于1。在步骤306中，考量像素i，在步骤308中找到像素i的特征向量 F1。特征向量是一个或多个特征的排序的一组值，每个值取决于所考量的像素的灰度值，和 /或在周围邻域内的其他像素的灰度值。该邻域不需要是邻接的。被考量的像素的坐标， 例如Xi和Yi (在二维图像)，或Xi、yi和Zi (在三维图像)，也可被处理为特征。本领域已知 的特征的例子包括像素i的灰度值，采用上述的双边过滤器，以及在围绕像素i的特定尺 寸的邻域内的每个像素的灰度值，采用非本地过滤器装置。如上所述，本领域已知的其他特 征包括在像素i的邻域内的所有像素的平均灰度值，在像素i的邻域内的灰度值梯度的方 向，以及在像素i的邻域的二进制边缘图的每个像素的灰度值，正如采用Carmy边缘检测器 确定的。正如将在下面详细描述的那样，根据本发明的实施例，可定义宽的变化范围的其他 特征。从步骤310起，检查一组检索像素，标记为检索像素j，以便找到具有与像素i类 似特征值的像素。检索像素j的灰度值非常近似于像素I，这将最大程度地有助于估计像 素i的没有噪声的真实灰度值。在步骤310，指数j被初始设为等于1。在步骤312，考量 检索像素j。检索像素可选地包括在该图像中的所有像素，或除了像素i以外的所有像素。 另外，检索像素仅包括在该图像内的像素的子集，例如，仅包括在检索窗口内围绕像素i的 像素，或仅包括随机选择的一些像素，或者在该检索窗口内规则间距的像素，和/或具有足 够接近于像素i的灰度值的像素。可选地，例如在医学图像中，该图像是分割为不同类型的 组织，采用任意已知的分割技术，仅或优选地，从如像素i的同类组织的像素中选择检索像
ο此外，检索像素可从其他图像的像素字典中选择，其他图像被期望为类似于本图 像。例如，如果本图像是一个医学图像，字典包括来自从同一病人身体的相同部分制备的早 期图像的像素，或者来自其他病人身体的相同部分的像素。在步骤314，评估对于检索像素j的特征向量F2。特征向量F2是一个或多个特征 的排序的一组值，每个值对应于特征向量F1的一个特征值。可选地，以同样方式定义在F1 和F2中的对应特征，采用在像素i和像素j的邻域的对应像素的灰度值。在本发明的一些 实施例中，在F1和F2中的对应特征值是不同定义的，例如，围绕其中一个像素的邻域可被指 向不同的角度，或以不同尺寸缩放，相对于围绕其他像素的邻域，如果需要，在计算特征值 时插入灰度值。在任何情况下，在中的对应特征都是可选地以类似的足够方式来定 义，以致能够比较它们，也能采用它们的值的差来计算在像素i和像素j之间的理论距离测 度，可知它们之间的相似度如何，用于降低噪声的目的。如果检索像素j是从之前存储的检索像素字典中取得，而不是从正被检查的图像 中取得，则对于像素j的特征向量F2,或者它的一些分量，也被存储在该字典中，在使用时不 需要每次计算。类似地，如果检索像素j之前是作为对于另一像素i的检索像素来使用的， 则它的特征向量&可选地是存储在内存中，不需要再次计算。可选地，对于在图像中的所有
11像素，特征向量F2是提前评估的，并存储在内存中，因此F2不需要在循环中对于检索像素j 和像素i而被评估。以至于F2的特征值是以相同方式定义的，作为F1的对应特征值，特征 向量F2或它的一些分量，也可从内存中取回，而不是再次计算，如果检索像素j之前采用为 在步骤306中检查的像素i的话。在步骤316中，距离测度(!(F1, F2)是可选地计算的，它是反映像素j与像素i的相 似度的理论距离，由它们的灰度值和它们邻域的灰度值来定义，也可能由它们的位置来定 义。距离测度d取决于构成特征向量F1和F2的每个对应特征值的差异。如果特征向量F1 和F2的每个特征值具有由F1 = ( ,…，力)与巧=( ，...，//)给出的k分量(特征值)， 则该距离测度可由下式定义J(F15F1) = [a,I/,1 -f^+a2\f^-f22\"+... + ak\fk - fk f这里，(Q1, α2, ... Qk)是加权向量，给予在计算距离测度中用于不同特征的权 重。参数β通常是同一级的正数，通常设为等于2，这使(!(F1, F2),正交分量的欧几里得距 离，每个等于在两个像素i和j的特征值之间的加权绝对差。正如下面所描述的，加权向量 (a” a2，... ak)是可选地采用遗传算法找到的，这试图找到一个最佳的加权向量以使降 低噪声方法的效率最大化。对于(!(F1, F2)的另一种表达，考虑在不同特征值之间的相关性，例如在一个邻域 中的不同像素的灰度值之间，在下面图8的说明中进行讨论。对于(!(F1, F2)的表达可包括 交叉项，例如(Zr11 - /rI2KZr21 - &2)，可提供对于在不同邻域之间的相似度的更有用的测度， 在不同特征值是相关的情形下。在步骤318中，对于像素j的权重Wj是可选地从d (F1，F2)计算的，并存储在内存。 当像素i和j的邻域彼此最大相似时，也就是，当d是小的时，权重Wj是最大的，而当d是 大的时，Wj是小的。例如，Wj = eXp(-d2/0N)。如果特征值仅取决于像素和它的邻域的灰 度值，而不取决于该像素的位置，则Wj是是由Wj = exp (-d2/ σ N-dp2/ σ ρ)定义的，这里dp是 在像素i和j之间物理距离的测度，例如，欧几里得距离。这里，。N和％是决定伴随在像 素i和j之间理论距离d和空间距离dp的增加而Wj减弱的程度的参数。此外，Wj具有与 d与dp的不同的相关性，但仍随着d与dp的值的增大而减弱。可选地，为节省计算时间，或 者为增强性能，权重Wj是被设为零，当它小于某些极限时，或当d和/或dp大于某些极限 时。在步骤320，检索像素j被增加1，以检查下一个检索像素。在步骤322，确定所有 检索像素是否已经被检查完。如果否，在步骤312考量下一个检索像素。如果所有检索像 素已经被检查完，检索像素j的加权平均灰度值(由Wj加权)将被计算。在步骤326，估计像素i的没有噪声的真实灰度值，基于检索像素的灰度值，可选 地也基于像素i的原始灰度值，检索像素具有对于估计的真实灰度值的更大影响，如果它 们被视为更近似于像素i，基于所具有的类似特征值。例如，特征值的相似度被用于计算理 论的距离测度(KF1, F2)，如前所述，每个检索像素j被分配一个权重Wj，该权重基于它的从 像素i的距离测度，并从检索像素j的加权平均灰度值(由Wj加权)中找到像素i的估计 的真实灰度值。该平均可以是平均数、中值、模式、去除轮廓的平均数，或任意其他类型的平 均。
12
另外，像素i的真实灰度值的估计是以不同方式计算的从检索像素的灰度值、检 索像素的特征向量F2，以及像素i的特征向量Fp例如，检索像素被分为不同级别，表示不 同组织类型，基于它们的特征向量F2的簇，只有在同类的检索像素(如像素i)被用于估计 像素i的真实灰度值，或者具有比像素i的估计的真实灰度值更大的效果。另外，只有顶部 的几个检索像素j具有接近于F1的特征向量F2 (通过一些测量)，这几个检索像素被用于 估计像素i的真实灰度值。可选地，代替采用检索像素的平均灰度值，从基于几个检索像素 的灰度的检查表找到像素i的估计的真实灰度值。可选地，校正的灰度值是像素i的原始灰度值与检索像素的加权平均的线性结 合。可选地，不明确考虑像素i的原始灰度值，但像素i自身被处理类似另一个检索像素，并 包括在加权平均。在本例中，如果F2的特征值是以与F1的对应特征值相同的方式来定义， 对于像素i自身的权重Wj将是1，如果F2的特征值是以不同的方式(例如，邻域旋转或缩 放)来定义，则对于像素i自身的权重Wj将会小于1。需要明确的是，这里所指的像素的灰度值不必然是该图像的原始灰度值，还可以 是变换的图像或过滤的图像的灰度值，例如，高斯过滤后的图像，其O等于或小于仅几个 像素的宽度。在步骤328中，像素i增加1，在步骤330中，将确定是否仍留有任何像素需要考 量。如果有，将回到步骤306考量下一个像素i。如果没有，该流程在步骤332结束，带有降 低噪声的图像，采用在步骤326中找到的校正的灰度值，作为输出。图4A显示了噪声图像400，头部切片的CT图像，以说明在图3中略述的方法。该 图像比普通图像更多噪声，因为它是采用减少的X-射线剂量获得的。图4B显示了降低噪 声的图像402，是采用图3所示的方法从图像400获得的，带有一组特征和加权向量，将在下 面描述。为了比较，图4C显示了一个低噪声的图像404，采用正常的对这类图像的X-射线 剂量获得的。减少噪声的图像402具有比原始图像400相当低的噪声，可看到更多细节，尤 其是在脑中，在不同组织之间有相对低的对比。图像402在质量上显得更接近与低噪声图 像404，与图像400相比。特征的示范类型特征的集中类型可被用在特征向量F1和&。在本发明的一些实施例中，计算一个或多个特征值的步骤包括找到在邻域内像
素的灰度值的分布特征。可选地，该特征值是灰度值的分布的矩，或者一个或多个矩的函
数，其中，分布的第一矩是平均数，第二矩是标准偏差，第三矩是偏斜度，等等。分布的第k
(γ/fc
矩，k> 1，可被定义为Mi=女其中In是在该邻域内第η像素的灰度值，
V η乂，
总和是在该邻域内N个像素的和，而M1是第一矩，也就是灰度值的平均数。另外，特征值 是，或者取决于，分布的次序统计，其灰度值对应于该分布的给定百分数。例如，特征值是 中间灰度值，它是在50%百分数处的灰度值。另外，采用不同百分数的灰度值，例如25%、 37. 5%、62. 5%或75%。可选地，采用中间百分数，例如在25%至75%中间，它具有潜在优 势特征值将是作为一个整体的邻域的特征，而不仅是在该邻域的几个轮廓像素。可选地， 如果从字典中选择检索像素，该字典包括其他图像的检索像素，带有不同的标准化的灰度 值，则这两个图像的灰度值是标准化的，因此它们可做有意义的比较，例如，在基于次序统计的特征的比较。仅取决于在邻域内像素的灰度值的分布特征的一个特征，尤其是如果该邻域是正 方形或者是相当各向同性的形状，具有潜在的优势该特征值对于在图像中的结构的方向 是相对不敏感的。例如，采用在图2B中图像204这样的一个特征，很可能产生一组类似像 素206的像素，这些像素具有接近于像素110的特征值，因为该特征值将主要取决于该像素 与暗区104的距离，而不取决于在暗区与亮区之间的边缘的局部方向。在另一方面，如果已 知特定的一部分图像具有面向某一特定方向边缘或纹理，例如从身体组织分割图，则采用 对结构方向敏感的特征是有优势的。可选地，从在邻域内像素的原始灰度值的分布中未找到该特征值，但从该图像已 被光滑化或其他方式处理之后的灰度值分布中可找到该特征值。在评估该特征值之前光滑 处理该图像，具有潜在的优势该特征值可更多地取决于在邻域的图像的结构特征，并对在 邻域的噪声更低敏感。可选地，对于这里所描述的任意类型的特征，所述光滑或其他图像处 理时在评估该特征值之前完成的，而不仅用于取决于灰度值分布的那些特征。例如，可以通 过高斯过滤器、二进制过滤器或者全变差过滤器(如前面引用的Rudin等人的文献)来完 成平滑处理。可选地，该平滑以某种方式完成，在最大尺寸规模的用于该特征的邻域内不会 平滑掉大多数结构，或者甚至用在最小尺寸规模的邻域内。例如，如果以宽度参数ο使用 高斯过滤器，则σ可以是比邻域的最大尺寸或者邻域的最小尺寸更小，或者至少不太大。 另外，以某种方式完成平滑处理，该方式有效地平滑去除邻域的所有空间结构，而特征值是 不在该邻域内的结构的测量，但特征值是围绕该邻域的更大尺度的结构的测量，或者是围 绕该邻域的平均梯度的测量。可选地，在找到灰度值的分布之前，该图像是以不同的方式来处理的。例如，应用 导数算子到该图像，代替每个像素的灰度值，通过在特定方向与该图像的导数成正例的值， 或者通过与该图像的梯度幅度成比例的值。如果完成，则例如，在邻域内像素的平均分布值 将是在该邻域内平均梯度的测量。可选地，该图像是在找到该梯度之前被平滑的，足够平滑 以致在该邻域的大多数像素具有几乎相同的梯度，使该特征值对噪声低敏感。在本发明的一些实施例中，计算特征值的步骤包括将变换或过滤应用到该图像， 至少在一个邻域内，优先选择在该邻域的最大尺寸与几个像素之间的中间尺度范围内的结 构。另外，该变换或过滤优先选择一些方向而不是其他方向的结构。以这种方式定义的特 征可以是在该图像中有利于挑选期望具有在特定范围(例如，血管)内的尺寸和/或方向 的结构，而忽视由于噪声导致的精细的密度变化。本发明的这些实施例可采用任意大变化范围的过滤器或变换器(线性或非线性) 所应用的特征，这些已经被应用于诸如计算机手写识别或图像的自动分类对象，因此它们 能不依赖于文字描述而被检索到，但这些技术尚未用于降低图像噪声。这些特征会依赖于在邻域内图像对于小波滤波器的响应，例如Meyer或Gabor滤 波器、Laplacian与高斯金字塔过滤器，或者现有技术中的任意其他线性过滤器。这些过滤 器可以对在邻域内的具有特定方向和/或特定尺度的结构非常敏感。可选地，该过滤器仅 应用到该邻域。另外，过滤器可应用于比该邻域更大的区域，例如，应用在该邻域内一个或 多个像素的灰度值上，在该图像已被过滤后。这些选项和另外的方式也可应用到这里所描 述的任意其他类型的参与对图像像素进行过滤或变换的特征。
另外，该特征值取决于图像对高斯过滤器或其他平滑过滤器在不同尺寸参数ο工 与02之间的响应的差别的。在两个这样的过滤器之间的差倾向于选择在O1与O2之间 的中间尺度内结构，但不取决于这些结构的方向，如果该过滤器是各向同性的。以此方式定 义的特征是特别有用的，如果该图像具有在许多不同方向上类似结构的话。在本发明的一些实施例中，特征可取决于在邻域内图像对非线性变换(例如形态 多尺度变换或形态算子)的响应。例如，该特征值取决于正被检查的像素的灰度值，或者在 邻域的特定像素，在采用特定尺度参数将非线性多尺度变换应用到图像之后。可选地，该特 征值取决于像素的灰度值，采用两个或更多不同尺度参数，例如，对于两个不同尺度参数的 像素的灰度值的差。形态多尺度变换的例子包括形态小波和形态金字塔，例如描述在以下 文献:E. H. Adelson, C. H. Anderson, J. R. Bergen, P. J. Burt, and J. M. Ogden, " Pyramid Methods in Image Processing, " RCA Engineer29, no. 6, Nov. -Dec. 1984, pp. 33-41, 或者"Nonlinear Multiresolution Signal Decomposition Schemes-Part I Morphological Pyramids, " John Goutsias, and Henk J. A. M. Heijmans. IEEE TRANSACTIONS ON IMAGE PROCESSIN G，VOL. 9，NO. 11，NOVEMBER 2000。该特征也可取决于在应用形态算子后的像素的灰度值。将形态算子应用到图像以 增强或抽取特定结构。形态算子的一个例子是顶帽变换，在输入图像与它的由结构元素打 开的形态之间的差。这样一个算子将揭示在暗背景上的光亮细节，带有控制所检测特征的 尺寸的结构原件的尺寸。可定义一个类似的算子，它从白色背景中抽取暗结构。关于形状匹配与图像变形的文献包括宽范围的技术，可被用于使图像的形状特征 化，以及任意这些方法可被用于定义邻域的特征，在应用如前所述的形态变换或形态算子 Γ$ ^。 ^iJiySiSiEiirth Mover' sYossi Rubner, Carlo Tomasi, and Leonidas J. Guibas, " A Metric for Distributions with Applications to Image Databases, " Proceedings of the 1998International Conference on Computer Vision, Bombay, India ；用于图像变形的 Kantorovich—Wasserstein 度量，例如描述于 Steven Haker, Lei Zhu, Allen Tannenbaum and Sigurd Angenent, " Optimal Mass Transport for Registration and Warping, " International Journal of Computer Vision 60 (3)，225-240 (2004)；由以下文献定义的形状签名Robert Osada, Thomas Funkhouser,Bernard Chazelle and David Dobkin, " Matching 3D Models with Shape Distributions, “ ACM Transactions on Graphics 21，807-832 (2002)；以及由以下文献 定义的形状匹配的度量Serge Belongie, Jitendra Malik, and Jan Puzicha, “ Shape Matching and Object Recognition Using Shape Contexts, " IEEE Transactions on Pattern Analysis and Machine Intelligence 24，509-522 (2002)；所有上述文件都引入 作为参考。在本发明的一些实施例中，对于正被检查的像素i与检索像素j的相应特征值是 以由几何变换所改变的邻域来计算的，从一个邻域到另一个邻域。例如，两个邻域是处于不 同的相对方向和/或尺度。另外，其中一个邻域可以是相对于另一个邻域的镜像。例如，如 果找到像素i的特征值的算法使用了在该邻域的像素的灰度值，该像素与像素i在+X方 向有一特定距离，则像素j的特征值是采用替代的像素来计算的，该像素与像素j在+y方 向(90度旋转)有相同距离，或者在-χ方向(反射)有相同距离，或者在+χ方向(尺度改变)有两倍距离，或者在+y方向(旋转加尺度改变)有两倍距离，等等。可选地，旋转的和 /或缩放的邻域的灰度值是在计算该特征值之前插入的，尤其是如果旋转角度不是90度的 整数倍(在像素的笛卡尔网格的情形下)或者缩放因子不是整数。可选地，这些像素是安 排在一个三角形或者六角形网格内，或者更复杂的模格内。采用以这种方式定义的特征将是尤为有用的如果无需旋转和/或缩放和/或反 射而定义相同的特征，以及用于各种不同的旋转角度和/或缩放因子，相同的权重给出作 为结果的特征。这可导致一个距离测度，它取决于在图像中的结构的方向和/或尺度，至少 对于一些方向和/或尺度的范围。这将是有优势的如果该图像包括带有不同方向或尺度 的类似结构，或者图像的互相的镜像。优化特征的加权向量的遗传算法以上面所述的宽变化范围的可能特征，是很难系统地找到最佳的特征组以及相应 的权重α，该权重给出降低噪声的最好结果，对于特定的图像形态、特定的身体部分或疾病 状态、特定的对比材料的类型和浓度，或者其他图像特征。当不同特征能协同相互作用时， 是特别真实的，特定的一组特征是怎样有效工作的。图5显示了流程图500，显示了一种 遗传算法，它能被用于找到最佳的或至少较高的对于一组特征的加权向量。这里，“较高的 (superior)，，是指在特定的图像级别中对于降低噪声更有效的。需要明确的是，本领域已知的任何优化方法都可被用来找到加权向量。采用遗传 算法是有优势的，因为该组可能的特征具有大量的模块部件、个体特征，它们协同相互作 用，有点类似在有机体中的基因。在步骤502中，选择并排序一组特征。可选地，该组特征包含非常大量的特征，期 望得到遗传算法的结果，可被用于消除许多特征。以至于一些权重的找到率非常低，这将表 明相应的特征对于降低噪声并不是非常有用，而它们将从列表中清除，节省了运行降噪软 件的计算机时间。在步骤504中，选择噪声测试图像，以及等效的低噪声测试图像，以及更低噪声或 没有噪声的图像。将测试特定顺序的特征组和对应的加权向量的效力，通过将降噪算法应 用到噪声测试图像，计算输出的图像是如何接近于低噪声测试图像。如果该测试图像是从 与降噪软件能工作的图像相同的图像形态中获得时，遗传算法可运行得最好。图像的性质 和噪声的性质都可能影响特定的特征组和相应的加权向量的效力。可选地，该噪声测试图 像是在与低噪声测试图像在相同时间内获得的，但分别以不同的采集参数获得，例如，以低 X-射线剂量获得的CT图像，或者以较短的采集时间获得MRI图像。另外，噪声图像是通过 加入模拟噪声到低噪声图像而获得的。在步骤506中，创建加权向量的起始群域，例如采用随机数字产生器。每个加权向 量具有对于在顺序组中的每个特征的一个分量。可选地，该群域是足够大的，足以在可能的 加权向量的相空间内扩散，这是不可能的遗传算法将陷入局部最大效力，这是比全局最大 效力更坏。在另一方面，该群域是足够大的以致用于运行遗传算法的计算时间实际上不长。 以至于特别当测试大量特征时，很难确定需要使用的群域有多大，以及很难满意上述两种 情形。以至于可通过限制特征的数量来避免采用大量特征的计算困难，在群域的任意一个 成员中非零加权，假设合并多个不同特征的协同优势能获得，无需使用比在组中所有特征 的小片段更多的特征。验证这样一个假设，并获得一个想法什么构成足够的群域尺寸，可通过测试不同起始群域来完成，看最终结果是否非常不同。在步骤508中，评估每个加权向量的降噪效力，通过在噪声测试图像上运行降噪 算法，并将该结果与低噪声测试图像比较。加权向量的效力测量是，例如，在从降噪算法获 得图像与低噪声测试图像之间的均方误差。可使用均方误差的负数或倒数，如果想要定义 该效力测量，以致数字平均数越高，加权向量越有效。可选地，采用误差的视觉或生理学 测量，以测量该加权向量的效力；参见James L. Mannos and David J. Sakrison, “ The effect of a visual fidelity criterion on the encoding of image, " IEEE Trans, on Information Theory 20, no. 4, p.525 (July 1974)。在步骤520中，通过在群域中的所有加权向量中找到具有指示最大效力的效力测 量的加权向量，找到最大效力的加权向量。在步骤512中，确定遗传算法是否已经集中，或者是否已经找到效力足够的加权 向量。例如，通过比较最大效力的加权向量的效力测量与前代的数量，或者在前代的几个 量，看在每代中是否改变很大。另外，大多数有效加权向量的效力测量是与效力的一些绝对 标准进行比较。如果该遗传算法还未集中，则在步骤514创建新一代的加权向量。例如，通 过有更高效加权向量产生比低效加权向量更多的子孙，从而创建新一代，最低效力的加权 向量根本不会产生。可选地，另外，加权向量经历一些突变，不同加权向量的相同分量可不 同程度地交叉。采用有效的遗传算法可完成，也可采用本领域已知的其他遗传算法。在步骤508，评估新一代的每个加权向量，该程序继续用于后续的代，直到效力测 量开始聚集。一旦满意聚集标准，例如基于在对于过去几代的效力测量上的相对改变，流程 终止于步骤516，产生最大效力的加权向量作为输出。示例性特征组和加权向量表1显示了一组排序的特征，以及采用如图5所示的遗传算法找到的加权向量的 相应分量。特征列表不是一个全面的覆盖这里所述的不同的可能特征的列表，但包括了特 定类型的特征，当测试其他类型作为替代特征或者附加特征时，可导致更有效的降低噪声。 需要注意的是，在表1中大范围的顺序的权重幅度反映了特征值的相对幅度，它们是非标 准化的，不同特征的相对有效性也是如此，以致带有最高权重的特征是不需要作为最重要 的特征包括在内的。显示在图4A、图4B和图4C中的图像是全三维图像在χ-y平面的截面 图，采用表1所示的特征组合加权向量来降低它的噪声。全3D图像的体素在ζ-方向是比 在χ和y方向更长的。显示在图4A、图4B和图4C中的截面图在χ和y方向是约450乘500 像素。在表1中，“图像自身”是指像素的灰度值，以豪森菲尔德单位表示。不同顺序统计 特征的平均灰度值，以豪森菲尔德单位表示，特别是百分数，特别对于在邻域的灰度分布， 在邻域的中心带有检索像素。“均方误差”是在邻域内灰度分布的标准导数，以豪森菲尔德 单位表示。高斯过滤器特征涉及灰度响应的差分，以豪森菲尔德单位表示，带有不同参数的 两个高斯过滤器的中心像素，在毫米级测量。在涉及梯度的特征中，梯度单位是豪森菲尔德 单位每毫米。Ix、Iy和Iz分别是指梯度的χ分量、y分量和ζ分量，而梯度额是指梯度的幅 度。表1 特征与加权向量 诜择降低噪声时考虎的图像参数在大多数图像形态中，获得低噪声的图像是可能的，但不够经济或不够安全。例 如，在CT成像中，可通过采用较高X-射线剂量来降低噪声，但可能会增加成本，而高剂量的 X-射线会致癌。可以选择图像的噪声水平，例如，权衡以下方面要么获得用于诊断医疗问 题的较好图像，因而增加病人能被有效治疗的机会，要么增加病人患上癌症的机会，如果采 用高剂量的X-射线的话。对于MRI，可通过较高静磁场来获得较低噪声。已知高的静磁场对身体健康的不良 影响，高静磁场的MRI设备的造价也更高，只有少数医院能承受得起，也只有少数病人能获 得这样的待遇。全身治疗型MRI设备通常所用的最高磁场为约3特斯拉的孔磁场，而具有 更高磁场的设备将很难或不可能以已知技术的任何代价来建造，虽然小型设备有时可采用
19更高的磁场。带有1.5或2特斯拉的孔磁场的MRI设备是明显低廉的。在MRI图像中的噪 声也可通过采取更长的图像采集时间来降低，但该采集时间受到大多数病人所能承受的时 间的限制，而且增加采集时间也会降低病人通过昂贵MRI设备的机会，最终限制能接受MRI 检查的病人的数量。最后，包括CT和MRI，都可通过较低来获得较低噪声，但通常分辨率越低，噪声越 高，而高噪声会影响该图像用于医疗状况的诊断。图6显示了流程图600，描述了在CT、MRI或任意其他图像形态中选择图像参数的 方法，该方法补偿了采用如图3所述的方法来降低噪声的降噪软件的能力。流程图600所 示的方法并不限于采用如图3所述的降噪方法，而可采用本发明所述的其他方法。也可采 用现有的降噪方法，例如在图IA和图IB中描述的方法，或者目前已知或将来可知的其他降 噪方法。在步骤602，获得标准采集参数，例如CT的X-射线剂量和分辨率，或MRI的磁场、 采集时间和分辨率，都将产生想要的噪声水平。这些参数可以是特定成像设备的标准操作 指令的一部分，但它们也可取决于是否存在造影剂，取决于该图像是否被制作用于健康筛 选目的，较低质量的图像也可被证明，或用于检查已知的医疗问题，以及取决于诊断特定症 状所需的分辨率水平和噪声水平。在步骤604中，找到新的采集参数，例如采用检查表，在采用降噪软件后将给予相 同的噪声水平，至少在该图像的那些部分，噪音水平是诊断的关键。可创建这样一个检查 表，例如，通过模拟或通过以真图像进行测试，对于特定的图像形态，采用对降低噪声最有 用的特征组，作为标准采集参数的函数，或等同于作为想要的最终噪声水平的函数。该测试 或模拟也可以对于一个范围的分辨率来完成，或者它能简单地假设相对噪声水平与体素 体积或像素面积的反平方根成比例，这将是真的，如果在不同像素的噪声是无关联的。在步骤606，采用新参数获得一个图像，该图像具有比想要的噪声水平更多的噪 声，但将在CT图像或其他χ-射线图像中采用较低X-射线剂量，或者MRI图像中采用较低 磁场和/或较短的采集时间。在步骤608中，采用降噪软件降低在图像中的噪声，导致获得 具有相同的噪声水平的图像，或者用于医疗诊断的相同有用性的图像，与由无降噪的标准 采集参数来产生图像一样。采用残余图像图7显示了根据本发明的示意性实施例所述的流程图700，用于采用一个残余图 像，该残余图像介于原始图像与降噪图像之间，恢复了原始图像的一些非噪声部分，这些部 分本已由降噪软件去除。流程图700对于评估是有用的，通过反馈，改善了降低噪声的方 法。在步骤702中，获得原始噪声图像。在步骤704中，采用如图3所述的降噪方法从该图 像中去除噪声，或者采用这里所述的或现有技术中的任意其他降噪方法。在步骤706中，通 过从原始图像中减去降噪图像，产生残余图像。理论上，该残余图像将是纯噪声的。然而实 际上，在该残余图像中通常有一些低的对照结构，显示降低噪声不会仅去除噪声，还会去 除在原始图像中的一些结构，作为降低噪声的结果。需要注意，取决于图像形态，即使真正噪声也不同于在图像中点到点的振幅或其 他特征，导致在残余图像中可能大尺度的结构，即使仅从原始图像中去除噪声。为避免大尺 度结构从原始图像中不正确地去除，图7所示方法一次仅作用于图像的一部分，在该部分内的噪声是均一的，或者噪声水平是相对低的。在步骤708中，采用非线性保留边缘过滤器来平滑残余图像和定位边缘。该过 滤器降低在残余图像中的噪声，使边缘平滑为平行于该边缘的表面，从原始图像中相应 增强该残余图像所包含的任意结构。合适的过滤器包括非线性的各向异性流过滤器，例 如Beltrami流和连贯增强过滤器，例如，描述在以下文献J.Weickert，‘‘ Anisotropic diffusion processing, “ Teubner 1998, and J. ffeickert, " Coherence enhancing diffusion filtering, " International Journal of Computer Vision 31, no. 2-3, pp. 111-128(1999)。在步骤710中，“结构(structureness) ”参数是评价为在已过滤的残余图像中的 位置的函数，它具有较高的值，如果该图像在该点具有一致的结构。有许多参数可被用于这 个目的。例如，结构张量的特征值，或Hessian矩阵的特征值，都可被采用。在步骤712，已 过滤的残余图像的具有较高“结构”值的部分将恢复到降噪图像，以产生改良的降噪图像。在距离测度中采用相关矩阵图8显示在图像中用于降低噪声方法的流程图800，对于距离测度(!(F1, F2)采用 比上面所给出的不同的表达。在某些方面，新的表达式是上面给出的表达式的总结，它考虑 了在不同特征值之间的可能的关联性。我们认为，仅当β =2，上述表达式变成
_1] Ci(FxiF1) = [ ,(//-Z12)2 -Jff +... + ,(/；-fl)f对于(!(F1, F2)的新表达式还包括在不同特征之间的交叉项，考虑在不同特征值之 间的相关性。对于(KF1, F2)的新表达式的潜在优势是它保持不变，如果该组特征值是由 一组相同数量的新特征值所代替，每个新特征值被定义为旧特征值的独立的线性结合。因 为新的特征值组会包含在旧的特征值组中的所有信息，对于(KF1, F2)的新表达式在某些意 义上仅取决于在该特征值中的信息，而不取决于定义该特征的特别方法。对于(KF1, F2)的 新表达式的另一个潜在优势是加权向量因子α，以及交叉项因子，都取决于在图像中的 位置，但不需要是相同位置，反应了在该图像的不同位置上的噪声的不同水平，以及不同的 噪声关联性。例如，在一个CT图像中，噪声水平总体上不同于在图像的不同部分，由于通过 I-D或2-D映射来形成这种图像。在MRI图像中，噪声水平可以是不均勻的，取决于RF域 的局部密度，也取决于在静磁场中可能的多相形。从图像的较大噪声部分给予像素较小的 权重，采用类似像素的灰度值来估计真正的灰度值将是什么，对于特定的像素，将会没有噪 声。在图8中找到的距离测度取决于在不同特征值之间的相关矩阵。对于本例，它将 被显示了 特征值是个别像素在该像素的邻域内的灰度值，特征正被找到。在本例中，对于 相等的权重α，上述的CKF1, F2)的上述表达是在非局部平均数方法中采用的均方误差。在邻域内第m个像素和在邻域内第η个像素之间矩阵系数Σ 由下式表达Σ mn = E[(Im-y) (Ιη-μ)]这里Im和In都是在邻域内第m个像素和第η个像素的灰度值，μ是在该邻域的 所有像素的平均灰度值，而E表示期望值。可选地，通过测量在均一图像内的灰度水平，以 致在灰度水平与平均灰度水平之间的局部系数都只归因于噪声。将邻域的灰度值表达为向 量N，相关矩阵Σ ■可在向量表达法中表达为
21
Σ= Ε[(Ν-μ)Τ(Ν-μ)]这里，μ表示一个向量，该向量的元素都等于μ。在步骤802中，对于图像的不同区域，找到平均灰度值μ。在步骤804中，对于该 图像的不同区域，找到相关矩阵的系数，例如前面定义的Σ mn。在本例中，这些区域非常小 于整个图像，但非常大于一个像素。可选地，系数Σ 是假设为以位置平滑变动，并采用插 入来再不同区域之间评估它们。在像素i和检索像素j之间的距离测度是由下式给出的dE (Ni, Nj) = [(Ni-Nj)1 · Σ · (Ni-Nj)1/2这里，Ni和Nj都是向量，分别表示在围绕像素i和像素j的邻域内像素的灰度 值。注意，i和j并不是标记Ni和Nj元素的指数，但它们是标记两个像素的指数，理论上 这两个像素互相远离，基于它们的邻域的相似性而被定义，例如，像素i是在降噪算法中被 检查，而检索像素j。Σ 1是相关矩阵Σ的倒数。在步骤806中，计算Σ的倒数。可选地，当该邻域具有大量像素时，为减少计算时 间，该矩阵因而是非常大的，其倒数是由下式近似的Σ 1^ZnVnVn7
η这里，λη是第η个特征值的倒数，而Vn是Σ的相应的特征向量。如果仅采用一些 最大的特征值和它们对应的特征向量，Σ的倒数可被快速找到。接着，距离测量是由下式接 近dz(Ni, Nj) (ΣΛ [(Ni - Nj).Vn]2)1/2
η这里，总和是超过特征向量。这个距离测度类似于之前定义在像素i和j之间的 理论上距离测度，对于特征值Vn之一，带有每个定义为Ni-Vn的特征，也就是在邻域内像素 的灰度值的特殊线性结合，由相应的特征值入 给出权重α。估计像素的噪声水平如前所述，因为噪声水平通常在图像中的位置间变动，在噪声区的检索像素j贡 献了比像素i相对少的估计的真正灰度值，例如，通过给它们比在较低噪声区的检索像素 更低的权重Wj。可选地，给出区域的噪声水平是通过在上述同一仿真的图像内在接近的像 素之间的关联的测量在来估计的。另一个方法，不仅提供了给定像素的噪声幅度的估计，而 且提供了该噪声的记号的估计，将在这里描述。对于较高的噪声检索像素，每个像素的噪声 的估计并不仅用于降低权重Wj，采用图3所示的降低噪声的方法，但也可直接用于降低噪 声，通过从每个像素减去估计的噪声。最终得到的图像时由它自身处理为降噪图像，或者它 是结合的，例如，通过采取加权平均，由图3所示的方法找到已降噪的图像，或通过任意其 他降低噪声的方法。图9所示的流程图900描述估计体素的噪声的方法。该方法是基于这样的想法 在无噪声图像中，几乎所有像素的补丁(合适尺寸)可能仅占有该尺寸的补丁的全部向量 空间的非常小的片段，将每个补丁当作灰度值的向量。一旦限定了向量空间所占有的部分， 任何从该部分的偏差都会假设是由噪声引起的。在步骤902中，检查一组一个或多个低噪声测试图像，可选地类似于已被使用的
22图像，例如，采用相同成像形态，可选地带有类似的图像采集参数，和/或看身体的相同部 分，如果它们是医学图像的话。在图9所示的例子中，图像是3-D的，因此像素是指体素，但 仍可采用2-D图像的方法。在每个测试图像中，检查一组尺寸为N1XN2XN3的体素的补丁， 可选地包括在该测试图像中的所有组的带有这些尺寸的补丁。每个补丁可表示为在一个理 论上N1XN2XN3维度向量空间内灰度值的一个向量，它的尺寸表示一个N1XN2XN3补丁的 体素之一的可能的灰度值。这些图像的补丁被分类为在这个向量空间的簇，采用现有技术 中已知的任意聚类算法。簇的数量是，例如，20至100，或者100至500，或者500至2000， 或者低于20，或者高于2000。采用大量的簇需要检查更多的补丁以获得关于向量空间是如 何分为簇的好的统计表，因而需要更多计算时间来对补丁进行分类，但采用过少的簇会导 致这些簇覆盖太多的向量空间，对于降低噪声是无用的。本发明人已经发现，采用数百个簇 可做得很好，当补丁尺寸在5至10体素，或10至15体素，或低于5体素，或高于15体素。 特别地，11x11x11体素的补丁表现为非常好。虽然在本例中给出的上述补丁是在矩形阵列的形状内，但这些补丁并不需要这个 形状，也可以有任意定义的关于体素的其他形状。向量空间的维度是在补丁内体素的数量。在步骤904中，检查每个簇以找到有限组的独立参数，明显少于N1XN2XN3需要覆 盖全向量空间的基向量的线性系数，这相当好地近似于覆盖在该簇内的补丁。例如，在该簇 内的每个补丁，或者最大量的补丁，是非常近似为灰度值的一个向量，该向量是固定用于在 该簇内的所有补丁，加上有限组的基向量的线性结合，带有不同于在该簇内的不同补丁的 系数。另外，在该簇内的大多数补丁是有灰度值的固定向量加上有限数量的独立参数的一 个函数(并不必需线性函数)来近似的。基向量系数的数量或者其他独立参数，对于所有簇 不需要是相同的。基向量或者独立参数的数量(例如，用于给出的簇)通常小于&乂队乂队 的三倍或小于N1XN2XN3的一半或小于N1XN2XN3的四分之一或小于N1XN2XN3的八分之 一，或者大于N1XN2XN3W三倍。一些簇可具有比典型数量更多基矢量或更多独立参数。目 标是识别相位空间的小片段，该片段提供了对在该簇内大多数补丁的很好的近似。表示补 丁向量的非线性方法包括但不限于核PCA、Laplacian eigenmaps和扩散图。在步骤906中，获得噪声降低的图像。一旦补丁的簇和基函数或者参数组已经在 步骤902和904中找到，不得不再次找到它们以用于处理附加的图像，但该程序可在步骤 906中开始。在步骤908中，图像体素i开始一个循环，初始体素i的指数设等于1。在步骤910 中检查围绕第一体素的补丁(相对于它们的体素质疑来定义这些补丁)，并决定该补丁归 属于哪个簇。另外，也可采用图3所示的方法，一起检查群体的体素，取代每次检查一个体 素的方式，而这些补丁是相对于它们的体素的群来定义的，或者相对于插入在体素之间的 虚拟体素来定义。如果这些补丁挑选得合适，则测试图像是足以代表正在处理的图像，这些 测试图像是没有许多噪声的，但充满了定义好的结构，而正被处理的图像也不由噪声支配， 则正被处理的图像的大多数补丁可干净地落入一个簇或另一个镞。如果一个或多个这些情 形不够满意，该方法仍是有用的，只是在降低噪声方面的效力低。在步骤912中，找到该簇的对于基函数的系数或其他独立参数，提供了最好的近 似，或至少理论上最好的近似，到体素i的补丁。对于大多数补丁的近似越好，降低噪声的 效果越好，因为本方法假设任何在观察的灰度值与它的近似之间的差是仅由于噪声引起的。在步骤914中，估计体素i的灰度值的噪声分量，通过考虑在体素i的实际灰度值与体 素i将有的灰度值之间的差，采用根据簇的固定向量加上有限组的基函数(或有限组的独 立参数的函数)的线性结合的近似。这个程序的基本原理是在缺乏噪声的情况下，补丁可 由近似来很好地描述，在测试图像中找到的在该簇内的所有或大多数补丁都运作良好，因 此任意差是很可能主要由于噪声引起的。在步骤916中，体素i被增加1。在步骤918中，确定在该图像中的所有体素是否 都已经被检查过。如否，在步骤910开始下一个体素的检查。如果所有体素都已经被检查 过，则在步骤920结束该流程，输出估计的每个体素的噪声。可选地，估计的噪声被用于降 低该图像的噪声，要么直接从该图像中减去所估计的噪声，要么间接地例如通过估计的噪 声水平影响权重Wj，当采用图3所示的降噪方法，或通过采用两种方法的结合。这里所用的术语“约”是指士 10%。术语“包括”、“包含”、“具有”以及它们的连接词是指“包括但不限于”。这个词组
包含术语“包括有”和“主要含有”。词“主要含有”意指其成分或方法可包括附加成分和/或步骤，但这些附加成分和 /或步骤不会从本质上改变所要求的成分或方法的基础特征和创新特征。这里所用的“一个”和“所述的”包括复数，除非另有明确说明。例如“一个化合物” 或“至少一个化合物”可包括复数个化合物，包括它们的混合物。这里所用的词“示范的(exemplary)”是指“用作一个例子或例证”。任何描述为 “示范的”实施例并不是必须作为比其他实施例更好或更有优势的例子，和/或排除从其他 实施例中引入的特征。这里所用的词“可选地(optionally) ”是指“以一些实施例来提供，不以其他实施 例来提供”。本发明的任何特定实施例可包括多个“可选的”特征，除非这些特征发生冲突。通过本申请，本发明的不同实施例可在一个范围内展示。需要明确的是，在该范围 内的描述仅是为了方便和简短，而不应视为对本发明的范围固定的限制。相应地，描述的范 围应当被视为具有特别揭示了所有可能的子范围以及在该范围内的个别数值。例如，例如 从1至6的范围的描述应当被视为也揭示了例如1至3、1至4、1至5、2至4、2至6、3至6 等子范围，以及在该范围内的单个数字，例如，1、2、3、4、5和6。不论范围有多宽，都可这样 来应用。无论一个数字范围在这里如何指示，它是指包括在该指定范围内任意引用的数字 (部分或整体)。这里所用的词组“归类”或在第一指示数字与第二指示数值“之间归类”以 及将第一指示数字“归类”到，或从第一指示数字“归类”到第二指示数字是可交换使用的， 它们的意思是指包括第一和第二指示数字以及在两个数字之间的部分或全部数字。需要明确的是，为表述清晰，本发明的特定特征是描述在分开的实施例中的，但也 可在单个实施例中结合提供。相反，为表述清晰，本发明的不同特征是描述在单个实施例中 的，但也可分开提供或者在任意合适的子组合中提供，或者在本发明的任意其他描述的实 施例中以合适的方式提供。描述在不同实施例中的特定特征不应视为那些实施例的必要特 征，除非该实施例缺乏那些特征元素就无法操作。虽然本发明已经结合特定实施例进行了描述，对于本领域技术人员显然还有许多 变动、修饰和改变。相应地，所有这些变动、修饰和改变都将落入本发明的精神与所附的权利要求的范围之内。 在本说明书中提及的所有出版物、专利和专利申请都在这里整体引入作为参考， 即使每个个别的出版物、专利或专利申请时特殊的和单独指示为引入作为参考。另外，在本 申请中的任意参考文献的引用或鉴定都不应认为该参考文献是本发明的现有技术。以至于 所用的章节标题，它们也不应视为必要的限制。
权利要求
一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻域的对应特征的值；以及c)基于图像的检索像素的原始灰度值或者变换灰度值，计算对于正被检验的每个像素的降低的噪声灰度值，相对于检索像素的一个或多个特征值，正被检验的像素具有更大敏感性的类似特征值；其中，所述计算至少一个特征值的步骤包括计算在邻域的像素的原始灰度值或变换灰度值的特征性分配，除了计算在该邻域的所有像素的平均灰度值之外。
2.根据权利要求1所述的方法，其特征在于，所述计算特征值的步骤包括计算第二矩 分布或更高矩分布。
3.根据权利要求1所述的方法，其特征在于所述在邻域的像素的灰度值是对于那些 像素的图像的原始灰度值。
4.根据权利要求1所述的方法，其特征在于所述在邻域的像素的灰度值是由过滤器 变换的图像的灰度值
5.根据权利要求4所述的方法，其特征在于所述过滤器是高斯过滤器。
6.根据权利要求1所述的方法，其特征在于，所述计算特征值的步骤还包括计算所述 像素在第二邻域的灰度值的至少一个特征性分布，该第二邻域是所述邻域的适当的子集。
7.一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻域的 对应特征的值；c)计算对于正被检验的每个像素的降低的噪声灰度值，基于它的检索像素的原始灰度 值或者变换灰度值，检索像素的一个或多个特征值比正被检验的像素的类似特征值具有更 大的敏感性；其中，所述计算至少一个特征的值的步骤包括评估在邻域的图像的变换，或者评估在 邻域的图像上的过滤器的效果，以致优先选择在邻域的尺寸和一些像素尺寸之间的中间尺 度范围的结构，或者一些方向而不是其他方向的结构，或者兼有这两种结构。
8.根据权利要求7所述的方法，其特征在于，所述计算特征值的步骤包括评估线形过 滤器对于所述图像的响应。
9.根据权利要求8所述的方法，其特征在于所述线性过滤器是一个小波滤波器。
10.根据权利要求8所述的方法，其特征在于所述线性过滤器是一个方向过滤器。
11.根据权利要求10所述的方法，其特征在于，所述计算特征值的步骤包括评估在超 过一个方向上的过滤器的效果。
12.根据权利要求7所述的方法，其特征在于，所述计算特征值的步骤包括评估在所 述邻域的图像的非线性形态变换。
13.根据权利要求12所述的方法，其特征在于所述的形态变换是多尺度变换。
14.根据权利要求12所述的方法，其特征在于所述的形态变换是形态算子。
15.根据权利要求7所述的方法，其特征在于在所述邻域的方向的角上所述特征的值的相关性具有对称性，以致如果所述邻域是旋转至少一个角度，对于在该邻域的图像像素 的任意灰度值，所述特征值是不变的。
16.根据权利要求7所述的方法，其特征在于对于至少一个特征，正被检验的像素的 特征值与检索像素的对应特征值是以它们在不同相对方向或尺度或两者兼有的邻域来评 估的。
17.根据权利要求7所述的方法，其特征在于所述减少的噪声灰度值是基于所述检索 像素的加权平均灰度值。
18.根据权利要求1所述的方法，其特征在于所述加权平均值采用的权重取决于在检 索像素的特征值与正被检验的像素的对应的特征值之间的差的测算。
19.根据权利要求1所述的方法，还包括a)识别残留图像的相对更多结构部分，表示从原始图像到减噪图像的变化；以及b)存储所述结构部分到减噪图像，以产生一个改善的减噪图像。
20.根据权利要求1所述的方法，还包括估计一个或多个检索像素的噪声水平，其中， 所降低的噪声灰度值是低相关性或不相关于具有更高的估计噪声水平的检索像素。
21.根据权利要求1所述的方法，还包括a)在图像中分配像素的一个补丁到在这些补丁的向量空间内的复数个簇之一，该补丁 表示为灰度的向量；b)近似所述补丁作为所降低的一组非相关参数的函数，参数量少于在所述补丁内的像 素量，与所述簇关联；以及c)根据在所述补丁的实际灰度与作为所降低的一组非相关参数的函数的所述补丁的 近似之间的差，校正所述补丁的像素的灰度。
22.根据权利要求21所述的方法，其特征在于所述函数是所降低的一组基向量的线 性组合，而非相关参数是在该线性组合内的基向量的系数。
23.一种获得带有改变的图像采集参数的医学图像的方法，包括a)确定一组一个或多个图像采集参数，这些参数会在图像降噪后达到想要的噪声水平；b)获得带有所述采集参数的图像；以及c)根据前述任意权利要求所述的方法降低该图像的噪声。
24.根据权利要求23所述的方法，其特征在于所述获得图像的步骤包括获得CT图像；以及所述确定一组一个或多个图像采集参数的步骤包括确定一组与该图像未降低噪声时 的情况相比具有降低的χ-射线剂量的参数。
25.根据权利要求23所述的方法，其特征在于所述获得图像的步骤包括获得MRI图像；以及所述确定一组一个或多个图像采集参数的步骤包括确定一组与该图像未降低噪声时 的情况相比具有降低的静磁场或减少的采集时间的参数。
全文摘要
一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻域的对应特征的值；以及c)基于图像的检索像素的原始灰度值或者变换灰度值，计算对于正被检验的每个像素的降低的噪声灰度值，相对于检索像素的一个或多个特征值，正被检验的像素具有更大敏感性的类似特征值；其中，所述计算至少一个特征值的步骤包括计算在邻域的像素的原始灰度值或变换灰度值的特征性分配，除了计算在该邻域的所有像素的平均灰度值之外。
文档编号H04N1/38GK101919230SQ200880123648
公开日2010年12月15日 申请日期2008年12月25日 优先权日2007年12月25日
发明者伊利兰·达昂, 兹维·德维尔, 埃泽尔·巴拉维夫, 塔尔·克尼希, 盖伊·罗斯曼 申请人:梅迪奇视觉-脑科技有限公司
图像噪声的降低相关申请本申请要求根据美国专利法35USC 119(e)的申请日为2007年12月25日的美国 临时专利申请61/016，578的权益。上述专利申请文件的所有内容都以其整体引入在此作为参考。发明的技术领域与
背景技术：
本发明，在它的一些实施例中，涉及一种降低图像噪声的方法，尤其是，但不完全 是，涉及一种采用非线性过滤器降低医学图像噪声的方法。下列出版物和发明涉及图像噪声降低方法、图像采集和/或计算机视觉。US 2007053477——锥形束CT成像与扇束CT成像的全球去噪方法与装置KR20050031210——图像去噪的方法与装置JP2000050109——用于去除噪声的非线性图像过滤器US6459755——用于管理低剂量CT扫描的方法与装置
US2003099405——降低CT剂量过滤器，带有计算有效实施EP1774837——降低活性剂量的装置与方法JP200139874——用于MRI的磁场发生器W02007047599——用于高增益磁共振的方法与装置G01V3/00, G01R33/34Optimal Mass Transport for Registration and Warping, International Journal of Computer Vision, Volume 60, Issue 3(December 2004), Pages 225-240, Steven Haker, Lei Zhu, Allen Tannenbaum, Sigurd Angenent A Metric for Distributions with Applications to Image Databases, ICIP 1998, Pages 59-66, Rubner Yossi, Tomasi Carlo, Guibas, J. Leonidas. Shape Matching and Object Recognition Using Shape Contexts, IEEE T-PAMI, Volume 24, No.4, (April 2002), Belongie Serge, Jitendra Malik, Puzicha Jan. Matching 3D Models with Shape Distributions, Proceedings of the International Conference on Shape Modeling & Applications 2001, Pages 154—166, Robert Osada, Thomas Funkhouser, Bernard Chazelle, and David Dobkin P. J. Burt, E. H. Adelson, " The Laplacian Pyramid as a Compact Image Code, " IEEE Trans, on Communications, pp.532__540, April 1983Iddo Drori,Daniel Cohen-Or,Hezy Yeshurun, ACM Transactions on Graphics 22(3)，(Proc. of SIGGRAPH 2003)，303-312.John Goutsias and Henk J. A. M. Hei jmans, "Nonlinear Multiresolution Signal Decomposition Schemes-Part I :Morphological Pyramids，，，IEEE Trans, on Image Processing, Vol. 9, No. 11,November 2000.John Goutsias and Henk J.A. M. Heijmans, "Nonlinear Multiresolution Signal Decomposition Schemes-Part II :Morphological Wavelets", IEEE Trans, onlmage Processing, Vol. 9, No. 11, November 2000.
Jean Serra, Image Analysis and Mathematical Morphology, 1982.
发明内容
本发明的一些实施例在一个方面涉及一种降低图像中的像素噪声的方法，通过估 计该像素的没有噪声的真实灰度值，基于其他类似像素的灰度值，识别为它们相邻像素的 类似特征。这些特征包括取决于在邻域的像素的灰度的分布函数的特征，或者相对于彼此 旋转或缩放的比较邻域的特征，或者参与应用变换或过滤到邻域的特征，以致优先选择中 间尺度范围的结构，或者仅是某些方向的特征。根据本发明的一个示例性实施例，提供了一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻 域的对应特征的值；c)基于图像的检索像素的原始灰度值或者变换灰度值，计算对于正被检验的每个 像素的降低的噪声灰度值，相对于检索像素的一个或多个特征值，正被检验的像素具有更 大敏感性的类似特征值；其中，所述计算至少一个特征值的步骤包括计算在邻域的像素的原始灰度值或 变换灰度值的特征性分配，除了计算在该邻域的所有像素的平均灰度值之外。在本发明的一个示例性实施例中，计算特征值的步骤包括计算第二时刻分布或 更高时刻分布。可选地，在邻域的像素的灰度值是对于那些像素的图像的原始灰度值。在本发明的一个示例性实施例中，在邻域的像素的灰度值是由过滤器变换的图像 的灰度值。可选地，所述过滤器是高斯过滤器。在本发明的一个示例性实施例中，计算特征值的步骤还包括计算所述像素在第 二邻域的灰度值的至少一个特征性分布，该第二邻域是所述邻域的适当的子集。根据本发明的一个示例性实施例，提供了一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻 域的对应特征的值；以及c)计算对于正被检验的每个像素的降低的噪声灰度值，基于它的检索像素的原始 灰度值或者变换灰度值，检索像素的一个或多个特征值比正被检验的像素的类似特征值具 有更大的敏感性；其中，所述计算至少一个特征的值的步骤包括评估在邻域的图像的变换，或者评 估在邻域的图像上的过滤器的效果，以致优先选择了在邻域的尺寸和一些像素尺寸之间的 中间尺度范围的结构，或者一些方向而不是其他方向的结构，或者兼有这两种结构。可选 地，计算特征值的步骤包括评估线形过滤器对于该图像的响应。可选地，所述线性过滤器 是一个小波滤波器。可选地，所述线性过滤器是一个方向过滤器。可选地，计算特征值的步 骤包括评估在超过一个方向上的过滤器的效果。在本发明的一个示例性实施例中，计算特征值的步骤包括评估在所述邻域的图 像的非线性形态变换。可选地，所述的形态变换是多尺度变换。可选地，所述的形态变换是 形态算子。
在本发明的一个示例性实施例中，在所述邻域的方向的角上所述特征的值的相关 性具有对称性，以致如果所述邻域是旋转至少一个角度，对于在该邻域的图像像素的任意 灰度值，所述特征值是不变的。在本发明的一个示例性实施例中，对于至少一个特征，正被检验的像素的特征值 与检索像素的对应特征值是以它们在不同相对方向或尺度或两者兼有的邻域来评估的。在本发明的一个示例性实施例中，所述减少的噪声灰度值是基于所述检索像素的 加权平均灰度值。在本发明的一个示例性实施例中，所述加权平均值采用的权重取决于在检索像素 的特征值与正被检验的像素的对应的特征值之间的差的测算。在本发明的一个示例性实施例中，所述的方法还包括a)识别残留图像的相对更多结构部分，表示从原始图像到减噪图像的变化；以及b)存储所述结构部分到减噪图像，以产生一个改善的减噪图像。在本发明的一个示例性实施例中，所述的方法还包括估计一个或多个检索像素 的噪声水平，其中，所降低的噪声灰度值是低相关性或不相关于具有更高的估计噪声水平 的检索像素。在本发明的一个示例性实施例中，所述的方法还包括a)在图像中分配像素的一个补丁到在这些补丁的向量空间内的复数个簇之一，该 补丁表示为灰度的向量；b)近似所述补丁作为所降低的一组非相关参数的函数，参数量少于在所述补丁内 的像素量，与所述簇关联；以及c)根据在所述补丁的实际灰度与作为所降低的一组非相关参数的函数的所述补 丁的近似之间的差，校正所述补丁的像素的灰度。可选地，所述函数是所降低的一组基向量 的线性组合，非相关参数是在该线性组合内的基向量的系数。根据本发明的一个示例性的实施例，提供了一种获得带有改变的图像采集参数的 医学图像的方法，包括a)确定一组一个或多个图像采集参数，这些参数会在图像降噪后达到想要的噪声 水平；b)获得带有所述采集参数的图像；以及c)根据前述任意权利要求所述的方法降低该图像的噪声。可选地，所述获得图像的步骤包括获得CT图像；以及所述确定一组一个或多个 图像采集参数的步骤包括确定一组与该图像未降低噪声时的情况相比具有降低的X-射 线剂量的参数。可选地，所述获得图像的步骤包括获得MRI图像；以及所述确定一组一个 或多个图像采集参数的步骤包括确定一组与该图像未降低噪声时相比具有降低的静磁场 或减少的采集时间的参数。除非另作定义，这里所采用的素有技术术语和/或科学术语都具有本发明所涉及 领域的技术人员所熟知的普通含义。与这里所描述的方法和材料的类似或等同的方法和材 料都可被用于本发明的实施例的实施或测试，下面给出了示例性的方法和/或材料。在发 生冲突时，本专利说明书包括这些定义，都会控制。另外，所有材料、方法和例子都仅是示例 性的，而不试图作为对本发明的限制。
6
本发明的实施方式的方法和/或系统的执行，可涉及手动、自动或两者结合地执 行或完成所选择的任务。而且，根据本发明所述的方法和/或系统的实际设备和装置的实 施例，可通过硬件、软件或固件或它们的结合，采用操作系统来执行几种选择的任务。例如，用于执行根据本发明所述的实施例的选择任务的硬件可以是以芯片或电路 的方式实施。对于软件，本发明所述的实施例的选择任务可采用任意合适的操作系统通过 计算机来执行多个软件指令来实施。在本发明的一个示例性实施例中，根据这里所描述的 方法和/或系统的示例的一个或多个任务是通过数据处理器来执行，例如，用于执行多个 指令的计算平台。可选地，该数据处理器包括用于存储指令和/或数据的易失性存储器，和 /或用于存储指令和/或数据的非易失性存储器，例如，磁性硬盘和/或可移动介质。可选 地，还提供了网络连接。显示器和/或用户输入设备例如键盘或鼠标也可选地提供。附图
简要说明在这里仅通过实施例并附图，来描述本发明的一些实施方式。现在特别详细参考 有关附图，需要强调的是，通过实施例的方式来显示特别内容，用于说明本发明的实施例的 示例性讨论的目的。在这方面，本说明书采用附图使本领域技术人员能清楚了解本发明的 实施方式。在附图中图IA示意性地显示了现有技术的没有噪声的两维图像，而图IB示意性地显示了 现有技术的带有噪声以及带有选择像素和邻域的相同图像；图2A示意性地显示了图IB的图像以及选择的像素和邻域，以及显示了现有技术 中的其他类似于所选择的像素的其他像素；图2B示意性地显示了图IB的图像以及选择的像素和邻域，以及显示了根据本发 明的一个示例性实施例的类似于所选择的像素的其他像素；图3是根据本发明的一个示例性实施例所述的一种降低图像噪声的方法的流程 图；图4A是采用相对低的X-射线剂量获得的噪声CT图像；图4B显示了采用图3的方法降低噪声的图4A的图像；图4C是类似于图4A中图像的一个低噪声CT图像，但该图是采用相对高的X-射 线剂量获得的；图5是根据本发明的一个示例性实施例所述的一种方法的流程图，该方法采用遗 传算法找到更有效的特征组加权到图3所采用的方法；图6是根据本发明的一个示例性实施例所述的一种获得图像的方法的流程图，该 方法采用改变的采集参数以补偿噪声降低，因为该图像可被变换为采用图3所示的降噪方 法的具有正常噪声水平的图像；图7是根据本发明的一个示例性实施例所述的在任意结构降噪后检查残留图像 的方法的流程图，该结构是从原始图像不正确地移除的，并恢复高结构化的部分以降低图
像噪声；图8是根据本发明的一个示例性实施例所述的一个方法的流程图，该方法找到一 个相关的矩阵，采用该矩阵找到在图3所示方法中的距离测度；以及图9是根据本发明的一个示例性实施例所述的一个方法的流程图，该方法估计一个图像的特定体素的噪声。优选实施方式详述本发明在一些实施例中涉及一种降低图像噪声的方法，尤其是，但不完全是涉及 一种采用非线性过滤器的降低图像噪声的方法。本发明的一些实施例的一个方面涉及一种降低图像噪声的方法，在该方法中，比 较不同像素的邻域，以便找到具有类似邻域的像素。采用一个或多个特征来决定这个类 似性，包括参与评估在邻域中原始灰度或变换灰度的特征性分布函数的特征的至少一个特 征。可选地，该特征性分布函数是矩分布函数，高于第一矩(平均灰度值)，例如第二矩(标 准偏差)、第三矩(偏斜度)或者更高矩的分布函数。可选地，所述分布函数是该图像的原 始灰度值的分布函数。此外，所述分布函数也可是变换后图像(例如，由高斯过滤器过滤后 的图像)的灰度的分布函数。可选地，可找到该分布函数的特征，既用于在邻域的所有像 素，也仅用于位于第二邻域的像素，该第二邻域是所述邻域的一部分。本发明的一些实施例的一个方面涉及一种降低图像噪声的方法，在该方法中，比 较不同像素，采用以下特征之一参与对该图像应用一个变换的特征，或者参与评估在图像 上过滤的效果的特征，以致在邻域尺寸和一些像素尺寸之间的中间尺度范围的结构，或者 优先选择一些方向而不是其他方向的结构。可选地，该过滤器是线性过滤器，例如小波滤波 器。可选地，该过滤器是方向性的。可选地，该变换是非线性形态变换，例如，多尺度变换， 或者形态算子。在本发明的一些实施例中，评估参与评估方向性过滤器在超过一个方向上 的反应的特征。可选地，该特征具有与角度相关的对称度，以致如果所述邻域是旋转至少一 个角度，例如180度或90度，不论该图像看起来像什么，所述特征值是不变的。本发明的一些实施例的一个方面涉及一种降低图像噪声的方法，在该方法中，比 较不同像素的邻域，采用至少一个特征值，该特征值是对于这些像素之一计算的，该像素的 邻域相对于其他像素的邻域旋转或缩放。在本发明的一些实施例中，邻域是定义用于像素群，而不是用于个别像素，而噪声 降低是作用在像素群，而不是作用在个别像素。此外，邻域可被定义在位于像素之间的体 素，而噪声降低是作用在该体素，这些体素的灰度值是由接近它们的像素的插值来定义的。本发明的一些实施例的一个方面涉及一种降低图像噪声的方法，在该方法中，在 图像的噪声降低后，检查残留的图像(在降噪前后的图像的差别)，看它是否具有结构元 件，如果有的话，结构元件被存储在该图像中。在本发明的一些实施例中，对于不同像素估计噪声水平，以及当像素无噪声的真 实灰度值是采用类似像素的灰度值来估计时，越多噪声的像素对所估计的灰度值的影响越 少。在本发明的一些实施例的一个方面涉及一种噪声降低系统，在该系统中，像素的 噪声分量是通过对一组测试图像的检查(第一检查)来估计的，将在这些图像上找到的像 素的补丁分类为在这些补丁的理论上的向量空间的簇，将这些补丁作为向量考虑，它们的 分量是它们的像素的灰度值。对于每个簇，发现降噪的一组独立参数(数量少于在该补丁 中的像素)，例如，基向量的一组降噪的线性系数，这提供了大多数补丁在该簇内的好的近 似，对于测试图像。当对一个噪声图像进行降噪处理时，围绕每个像素的补丁被分类为属于 其中一个簇，而它的真实值，无噪声，是采用对于该簇的降噪参数组来估计的。
在本发明的一些实施例的一个方面涉及选择用于采集图像的改变的参数，例如， 较低X-射线剂量的CT图像，或者较低静磁场或较短采集时间的MRI图像，利用降噪方法 来最终获得足够质量的图像以满足诊断需求，例如，采用更高的X-射线剂量或更高的静磁 场，没有噪声降低。举例来说，这里所描述的方法可被应用在图像采集设备或它的工作站(例如，CT 机、MRI机)上，应用在图像处理站和/或通过网络连接到远程位置或远程服务器。为了更好地理解本发明的一些实施例，正如在图2B至图9所示，参考传统(也就 是，现有技术)的图像降低噪声方法如在图IA至图2A中的操作。图IA显示了二维图像100，它包括一个像素数组，每组具有一个数值，该数值是在 黑与白之间的灰度值。在CT图像中的灰度值表示成像物体的实际密度的方便映射，通常表 示为豪森菲尔德单位(HU)。例如，在大脑的CT图像中，该图像通常显现0HU，这表示水的密 度，被映射为黑色，显现70HU，被映射为白色。总的来说，在图像处理文献中，术语“像素(pixel) ”是用于表示二维图像的一个元 素，而“体素(voxel)”是用于表示二维图像的一个元素。因为这里所描述的方法总体上既 可用于二维图像，也可用于三位图像，这里所采用的术语“像素”和“体素”不应理解为将本 说明书限制为二维图像或三维图像的例子。除非另有特别说明，这里所采用的术语“像素” 和“体素”应被理解为应用于任意例子的通用术语，它们通常可交替使用。这里所使用的术语“灰度值(grey value) ”不仅是指黑白图像的亮度，还是指彩 色图像中任意色彩变量的程度，例如，在彩色图像中的红、绿或蓝，或者彩色图像的亮度或 饱和度。在诸如CT或MRI图像的医学图像中，通常只有单一密度变量，例如1\或1~2加权密 度，它被映射到灰度图像的亮度，而在本例中，“灰度值”是特别倾向于，但这里所述的方法 并不限制于它们对医学图像或黑白图像的适用性。这里所描述的降低噪声方法法可被特别 用于医学图像，因为医学图像通常具有相对高的噪声水平，由于在噪声水平和图像采集参 数(例如，X-射线剂量或MRI采集时间)之间通常有折中，施加经济或安全惩罚以用于降 低噪声。而且，因为医学图像总体上没有在“光照(lighting)”差别，一个像素的邻域的特 征通常是它的真实灰度值的号的指示。这是特别真实的，因为医学图像倾向于具有类似结 构，在该图像的不同部分重复出现，有时带有改变的尺度或方向。图像100包括亮区102和暗区104，两区之间具有相当尖锐的边界。在图IB中， 图像108是已加入噪声的图像100。在现有技术中，噪声有时是由像素的灰度值与相邻像素 的灰度值的平均来降低的，处于最接近位置的像素给该像素最大的权重。这在均勻区域做 的很好，没有细节，例如在图像108中的亮区102和暗区104，但会导致在这两区之间的边 界的污点。另一个现有技术的降低噪声方法，双边过滤器，它是非线性过滤器，试图避免由 于像素i的灰度值Ii与其他像素j的灰度值Ij的平均导致的问题。例如，当操作位于(Xi， Yi)的特定像素i时，对位于(Xj，Yj)的另一个体素j的灰度值给予权重Wj,由下式给出 这里，dp是欧几里得距离，在空间中两个像素之间，而Ili-IjI可被认为是在这两 个像素之间的理论上“距离”，它们互相类似的程度的测度。对于像素i的新灰度值是由下式定义的
这里N是围绕像素i的搜索窗口，而总和是对于在该搜索窗口
得所有像素j。用于降低噪声的另一类型的非线性过滤器是在以下文献中定义的L. Rudin, S. Osher, and Ε. Fatemi, " Nonlinear total variation based noise removal algorithms, “ Physica D60，259-268 (1992)。在非本地过滤器装置中，现有技术的进一步发展，两个像素的类同之处取决于这 两个像素的邻域的像素与像素的比较。例如，为降低像素i (在图IB中标记为110)的噪声 水平，邻域Mi (在图IB中标记为112)，是围绕像素110来限定的。接着搜索其他像素j，其 具有围绕每个检索像素j的相同尺寸和形状的邻域Mj,在邻域112的像素与每个检索像素 j的邻域的相应像素之间找到均方误差MSE(Mi,Mj)。对于检索像素的均方误差是小的，在它 们的邻域与像素110的邻域之间，给出最大权重Wj，当平均检索像素的灰度值时，获得对于 像素110的降噪灰度值。该权重Wj由下式给出 对于像素i的新值是由下式确定 二卞J 7
ZjjeN J图2A显示了图像200，类似于在图IB中的图像108，具有一组像素202，它们具有 类似于邻域112或像素110的邻域。每个像素202具有类似的邻域，因为像素202都是离 开在亮区102和暗区104之间的边界同样的距离，方向都是在几乎相同的方向。在其他现有技术的降低噪声方法中，采用非线性过滤器，两个邻域的类似度是 基于在该邻域的所有像素的平均灰度值，或者在该邻域的像素的灰度值梯度的方向上， 如下文所述Mahmoudi, M. and Sapiro, G. , “ Fast image and video denoising via nonlocal means of similar neighborhoods, " IEEE, Signal Proc. , Vol. 12, no.12, pp. 839-842, Dec. 2005。在由 A. Heidarzadeh and A. N. Avanaki, "An Enhanced Nonlocal Means Algorithm for Image Denoising, "9th ISSPA，Feb. 2007所描述的不同方法中，两个 邻域的相似度取决于这两个邻域的二进制边缘图的均方误差，采用Carmy边缘检测器来确 定，还取决于在这两个邻域的原始图像的均方误差。在详细解释本发明的至少一个实施例之前，需要明确的是，本发明并不限制于以 下说明所提出的应用细节。本发明能采用其他实施例或者以不同的方式来实施或执行。图2B显示了图像204，类似于图像108。根据本发明的一个实施例所述，对于在两 个邻域之间的相似度，采用不同标准计算权重Wj，找到较好的一组检索像素206，它们具有 足够接近的类似于像素110的邻域112的邻域。在图2B所示的特殊例子中，标准(将在下 面详细描述)不依赖于这些邻域的相对方向，以致所有与暗区104有相同距离的像素(如 像素110)都由接近的类似邻域112的邻域，根据这些标准。与采用非本地装置方法带有高 权重的检索像素202相比，以高权重放大的检索像素206组可允许进一步降低噪声，因为有 更多像素来平均灰度值。在本发明的一些实施例中，对于在两个邻域之间的相似度的标准 可取决于这两个邻域的相对方向，而在本发明的这些或其他实施例中，具有高权重的检索像素的质量可以不必大于现有技术方法，但这些检索像素的质量可以更好，在这个意义上， 它们提供了比像素110的真实灰度值更好的估计。图3显示了根据本发明一个实施例所述的一种降低图像噪声的方法的流程图 300。流程图300所示方法是图2B所示方法的概括，具有用于在邻域之间的相似性的不同 标准。在步骤302，获得有噪声的图像。噪声降低算法在一次检查一个像素，初始在步骤 304中设定像素i等于1。在步骤306中，考量像素i，在步骤308中找到像素i的特征向量 F1。特征向量是一个或多个特征的排序的一组值，每个值取决于所考量的像素的灰度值，和 /或在周围邻域内的其他像素的灰度值。该邻域不需要是邻接的。被考量的像素的坐标， 例如Xi和Yi (在二维图像)，或Xi、yi和Zi (在三维图像)，也可被处理为特征。本领域已知 的特征的例子包括像素i的灰度值，采用上述的双边过滤器，以及在围绕像素i的特定尺 寸的邻域内的每个像素的灰度值，采用非本地过滤器装置。如上所述，本领域已知的其他特 征包括在像素i的邻域内的所有像素的平均灰度值，在像素i的邻域内的灰度值梯度的方 向，以及在像素i的邻域的二进制边缘图的每个像素的灰度值，正如采用Carmy边缘检测器 确定的。正如将在下面详细描述的那样，根据本发明的实施例，可定义宽的变化范围的其他 特征。从步骤310起，检查一组检索像素，标记为检索像素j，以便找到具有与像素i类 似特征值的像素。检索像素j的灰度值非常近似于像素I，这将最大程度地有助于估计像 素i的没有噪声的真实灰度值。在步骤310，指数j被初始设为等于1。在步骤312，考量 检索像素j。检索像素可选地包括在该图像中的所有像素，或除了像素i以外的所有像素。 另外，检索像素仅包括在该图像内的像素的子集，例如，仅包括在检索窗口内围绕像素i的 像素，或仅包括随机选择的一些像素，或者在该检索窗口内规则间距的像素，和/或具有足 够接近于像素i的灰度值的像素。可选地，例如在医学图像中，该图像是分割为不同类型的 组织，采用任意已知的分割技术，仅或优选地，从如像素i的同类组织的像素中选择检索像
ο此外，检索像素可从其他图像的像素字典中选择，其他图像被期望为类似于本图 像。例如，如果本图像是一个医学图像，字典包括来自从同一病人身体的相同部分制备的早 期图像的像素，或者来自其他病人身体的相同部分的像素。在步骤314，评估对于检索像素j的特征向量F2。特征向量F2是一个或多个特征 的排序的一组值，每个值对应于特征向量F1的一个特征值。可选地，以同样方式定义在F1 和F2中的对应特征，采用在像素i和像素j的邻域的对应像素的灰度值。在本发明的一些 实施例中，在F1和F2中的对应特征值是不同定义的，例如，围绕其中一个像素的邻域可被指 向不同的角度，或以不同尺寸缩放，相对于围绕其他像素的邻域，如果需要，在计算特征值 时插入灰度值。在任何情况下，在中的对应特征都是可选地以类似的足够方式来定 义，以致能够比较它们，也能采用它们的值的差来计算在像素i和像素j之间的理论距离测 度，可知它们之间的相似度如何，用于降低噪声的目的。如果检索像素j是从之前存储的检索像素字典中取得，而不是从正被检查的图像 中取得，则对于像素j的特征向量F2,或者它的一些分量，也被存储在该字典中，在使用时不 需要每次计算。类似地，如果检索像素j之前是作为对于另一像素i的检索像素来使用的， 则它的特征向量&可选地是存储在内存中，不需要再次计算。可选地，对于在图像中的所有
11像素，特征向量F2是提前评估的，并存储在内存中，因此F2不需要在循环中对于检索像素j 和像素i而被评估。以至于F2的特征值是以相同方式定义的，作为F1的对应特征值，特征 向量F2或它的一些分量，也可从内存中取回，而不是再次计算，如果检索像素j之前采用为 在步骤306中检查的像素i的话。在步骤316中，距离测度(!(F1, F2)是可选地计算的，它是反映像素j与像素i的相 似度的理论距离，由它们的灰度值和它们邻域的灰度值来定义，也可能由它们的位置来定 义。距离测度d取决于构成特征向量F1和F2的每个对应特征值的差异。如果特征向量F1 和F2的每个特征值具有由F1 = ( ,…，力)与巧=( ，...，//)给出的k分量(特征值)， 则该距离测度可由下式定义J(F15F1) = [a,I/,1 -f^+a2\f^-f22\"+... + ak\fk - fk f这里，(Q1, α2, ... Qk)是加权向量，给予在计算距离测度中用于不同特征的权 重。参数β通常是同一级的正数，通常设为等于2，这使(!(F1, F2),正交分量的欧几里得距 离，每个等于在两个像素i和j的特征值之间的加权绝对差。正如下面所描述的，加权向量 (a” a2，... ak)是可选地采用遗传算法找到的，这试图找到一个最佳的加权向量以使降 低噪声方法的效率最大化。对于(!(F1, F2)的另一种表达，考虑在不同特征值之间的相关性，例如在一个邻域 中的不同像素的灰度值之间，在下面图8的说明中进行讨论。对于(!(F1, F2)的表达可包括 交叉项，例如(Zr11 - /rI2KZr21 - &2)，可提供对于在不同邻域之间的相似度的更有用的测度， 在不同特征值是相关的情形下。在步骤318中，对于像素j的权重Wj是可选地从d (F1，F2)计算的，并存储在内存。 当像素i和j的邻域彼此最大相似时，也就是，当d是小的时，权重Wj是最大的，而当d是 大的时，Wj是小的。例如，Wj = eXp(-d2/0N)。如果特征值仅取决于像素和它的邻域的灰 度值，而不取决于该像素的位置，则Wj是是由Wj = exp (-d2/ σ N-dp2/ σ ρ)定义的，这里dp是 在像素i和j之间物理距离的测度，例如，欧几里得距离。这里，。N和％是决定伴随在像 素i和j之间理论距离d和空间距离dp的增加而Wj减弱的程度的参数。此外，Wj具有与 d与dp的不同的相关性，但仍随着d与dp的值的增大而减弱。可选地，为节省计算时间，或 者为增强性能，权重Wj是被设为零，当它小于某些极限时，或当d和/或dp大于某些极限 时。在步骤320，检索像素j被增加1，以检查下一个检索像素。在步骤322，确定所有 检索像素是否已经被检查完。如果否，在步骤312考量下一个检索像素。如果所有检索像 素已经被检查完，检索像素j的加权平均灰度值(由Wj加权)将被计算。在步骤326，估计像素i的没有噪声的真实灰度值，基于检索像素的灰度值，可选 地也基于像素i的原始灰度值，检索像素具有对于估计的真实灰度值的更大影响，如果它 们被视为更近似于像素i，基于所具有的类似特征值。例如，特征值的相似度被用于计算理 论的距离测度(KF1, F2)，如前所述，每个检索像素j被分配一个权重Wj，该权重基于它的从 像素i的距离测度，并从检索像素j的加权平均灰度值(由Wj加权)中找到像素i的估计 的真实灰度值。该平均可以是平均数、中值、模式、去除轮廓的平均数，或任意其他类型的平 均。
12
另外，像素i的真实灰度值的估计是以不同方式计算的从检索像素的灰度值、检 索像素的特征向量F2，以及像素i的特征向量Fp例如，检索像素被分为不同级别，表示不 同组织类型，基于它们的特征向量F2的簇，只有在同类的检索像素(如像素i)被用于估计 像素i的真实灰度值，或者具有比像素i的估计的真实灰度值更大的效果。另外，只有顶部 的几个检索像素j具有接近于F1的特征向量F2 (通过一些测量)，这几个检索像素被用于 估计像素i的真实灰度值。可选地，代替采用检索像素的平均灰度值，从基于几个检索像素 的灰度的检查表找到像素i的估计的真实灰度值。可选地，校正的灰度值是像素i的原始灰度值与检索像素的加权平均的线性结 合。可选地，不明确考虑像素i的原始灰度值，但像素i自身被处理类似另一个检索像素，并 包括在加权平均。在本例中，如果F2的特征值是以与F1的对应特征值相同的方式来定义， 对于像素i自身的权重Wj将是1，如果F2的特征值是以不同的方式(例如，邻域旋转或缩 放)来定义，则对于像素i自身的权重Wj将会小于1。需要明确的是，这里所指的像素的灰度值不必然是该图像的原始灰度值，还可以 是变换的图像或过滤的图像的灰度值，例如，高斯过滤后的图像，其O等于或小于仅几个 像素的宽度。在步骤328中，像素i增加1，在步骤330中，将确定是否仍留有任何像素需要考 量。如果有，将回到步骤306考量下一个像素i。如果没有，该流程在步骤332结束，带有降 低噪声的图像，采用在步骤326中找到的校正的灰度值，作为输出。图4A显示了噪声图像400，头部切片的CT图像，以说明在图3中略述的方法。该 图像比普通图像更多噪声，因为它是采用减少的X-射线剂量获得的。图4B显示了降低噪 声的图像402，是采用图3所示的方法从图像400获得的，带有一组特征和加权向量，将在下 面描述。为了比较，图4C显示了一个低噪声的图像404，采用正常的对这类图像的X-射线 剂量获得的。减少噪声的图像402具有比原始图像400相当低的噪声，可看到更多细节，尤 其是在脑中，在不同组织之间有相对低的对比。图像402在质量上显得更接近与低噪声图 像404，与图像400相比。特征的示范类型特征的集中类型可被用在特征向量F1和&。在本发明的一些实施例中，计算一个或多个特征值的步骤包括找到在邻域内像
素的灰度值的分布特征。可选地，该特征值是灰度值的分布的矩，或者一个或多个矩的函
数，其中，分布的第一矩是平均数，第二矩是标准偏差，第三矩是偏斜度，等等。分布的第k
(γ/fc
矩，k> 1，可被定义为Mi=女其中In是在该邻域内第η像素的灰度值，
V η乂，
总和是在该邻域内N个像素的和，而M1是第一矩，也就是灰度值的平均数。另外，特征值 是，或者取决于，分布的次序统计，其灰度值对应于该分布的给定百分数。例如，特征值是 中间灰度值，它是在50%百分数处的灰度值。另外，采用不同百分数的灰度值，例如25%、 37. 5%、62. 5%或75%。可选地，采用中间百分数，例如在25%至75%中间，它具有潜在优 势特征值将是作为一个整体的邻域的特征，而不仅是在该邻域的几个轮廓像素。可选地， 如果从字典中选择检索像素，该字典包括其他图像的检索像素，带有不同的标准化的灰度 值，则这两个图像的灰度值是标准化的，因此它们可做有意义的比较，例如，在基于次序统计的特征的比较。仅取决于在邻域内像素的灰度值的分布特征的一个特征，尤其是如果该邻域是正 方形或者是相当各向同性的形状，具有潜在的优势该特征值对于在图像中的结构的方向 是相对不敏感的。例如，采用在图2B中图像204这样的一个特征，很可能产生一组类似像 素206的像素，这些像素具有接近于像素110的特征值，因为该特征值将主要取决于该像素 与暗区104的距离，而不取决于在暗区与亮区之间的边缘的局部方向。在另一方面，如果已 知特定的一部分图像具有面向某一特定方向边缘或纹理，例如从身体组织分割图，则采用 对结构方向敏感的特征是有优势的。可选地，从在邻域内像素的原始灰度值的分布中未找到该特征值，但从该图像已 被光滑化或其他方式处理之后的灰度值分布中可找到该特征值。在评估该特征值之前光滑 处理该图像，具有潜在的优势该特征值可更多地取决于在邻域的图像的结构特征，并对在 邻域的噪声更低敏感。可选地，对于这里所描述的任意类型的特征，所述光滑或其他图像处 理时在评估该特征值之前完成的，而不仅用于取决于灰度值分布的那些特征。例如，可以通 过高斯过滤器、二进制过滤器或者全变差过滤器(如前面引用的Rudin等人的文献)来完 成平滑处理。可选地，该平滑以某种方式完成，在最大尺寸规模的用于该特征的邻域内不会 平滑掉大多数结构，或者甚至用在最小尺寸规模的邻域内。例如，如果以宽度参数ο使用 高斯过滤器，则σ可以是比邻域的最大尺寸或者邻域的最小尺寸更小，或者至少不太大。 另外，以某种方式完成平滑处理，该方式有效地平滑去除邻域的所有空间结构，而特征值是 不在该邻域内的结构的测量，但特征值是围绕该邻域的更大尺度的结构的测量，或者是围 绕该邻域的平均梯度的测量。可选地，在找到灰度值的分布之前，该图像是以不同的方式来处理的。例如，应用 导数算子到该图像，代替每个像素的灰度值，通过在特定方向与该图像的导数成正例的值， 或者通过与该图像的梯度幅度成比例的值。如果完成，则例如，在邻域内像素的平均分布值 将是在该邻域内平均梯度的测量。可选地，该图像是在找到该梯度之前被平滑的，足够平滑 以致在该邻域的大多数像素具有几乎相同的梯度，使该特征值对噪声低敏感。在本发明的一些实施例中，计算特征值的步骤包括将变换或过滤应用到该图像， 至少在一个邻域内，优先选择在该邻域的最大尺寸与几个像素之间的中间尺度范围内的结 构。另外，该变换或过滤优先选择一些方向而不是其他方向的结构。以这种方式定义的特 征可以是在该图像中有利于挑选期望具有在特定范围(例如，血管)内的尺寸和/或方向 的结构，而忽视由于噪声导致的精细的密度变化。本发明的这些实施例可采用任意大变化范围的过滤器或变换器(线性或非线性) 所应用的特征，这些已经被应用于诸如计算机手写识别或图像的自动分类对象，因此它们 能不依赖于文字描述而被检索到，但这些技术尚未用于降低图像噪声。这些特征会依赖于在邻域内图像对于小波滤波器的响应，例如Meyer或Gabor滤 波器、Laplacian与高斯金字塔过滤器，或者现有技术中的任意其他线性过滤器。这些过滤 器可以对在邻域内的具有特定方向和/或特定尺度的结构非常敏感。可选地，该过滤器仅 应用到该邻域。另外，过滤器可应用于比该邻域更大的区域，例如，应用在该邻域内一个或 多个像素的灰度值上，在该图像已被过滤后。这些选项和另外的方式也可应用到这里所描 述的任意其他类型的参与对图像像素进行过滤或变换的特征。
另外，该特征值取决于图像对高斯过滤器或其他平滑过滤器在不同尺寸参数ο工 与02之间的响应的差别的。在两个这样的过滤器之间的差倾向于选择在O1与O2之间 的中间尺度内结构，但不取决于这些结构的方向，如果该过滤器是各向同性的。以此方式定 义的特征是特别有用的，如果该图像具有在许多不同方向上类似结构的话。在本发明的一些实施例中，特征可取决于在邻域内图像对非线性变换(例如形态 多尺度变换或形态算子)的响应。例如，该特征值取决于正被检查的像素的灰度值，或者在 邻域的特定像素，在采用特定尺度参数将非线性多尺度变换应用到图像之后。可选地，该特 征值取决于像素的灰度值，采用两个或更多不同尺度参数，例如，对于两个不同尺度参数的 像素的灰度值的差。形态多尺度变换的例子包括形态小波和形态金字塔，例如描述在以下 文献:E. H. Adelson, C. H. Anderson, J. R. Bergen, P. J. Burt, and J. M. Ogden, " Pyramid Methods in Image Processing, " RCA Engineer29, no. 6, Nov. -Dec. 1984, pp. 33-41, 或者"Nonlinear Multiresolution Signal Decomposition Schemes-Part I Morphological Pyramids, " John Goutsias, and Henk J. A. M. Heijmans. IEEE TRANSACTIONS ON IMAGE PROCESSIN G，VOL. 9，NO. 11，NOVEMBER 2000。该特征也可取决于在应用形态算子后的像素的灰度值。将形态算子应用到图像以 增强或抽取特定结构。形态算子的一个例子是顶帽变换，在输入图像与它的由结构元素打 开的形态之间的差。这样一个算子将揭示在暗背景上的光亮细节，带有控制所检测特征的 尺寸的结构原件的尺寸。可定义一个类似的算子，它从白色背景中抽取暗结构。关于形状匹配与图像变形的文献包括宽范围的技术，可被用于使图像的形状特征 化，以及任意这些方法可被用于定义邻域的特征，在应用如前所述的形态变换或形态算子 Γ$ ^。 ^iJiySiSiEiirth Mover' sYossi Rubner, Carlo Tomasi, and Leonidas J. Guibas, " A Metric for Distributions with Applications to Image Databases, " Proceedings of the 1998International Conference on Computer Vision, Bombay, India ；用于图像变形的 Kantorovich—Wasserstein 度量，例如描述于 Steven Haker, Lei Zhu, Allen Tannenbaum and Sigurd Angenent, " Optimal Mass Transport for Registration and Warping, " International Journal of Computer Vision 60 (3)，225-240 (2004)；由以下文献定义的形状签名Robert Osada, Thomas Funkhouser,Bernard Chazelle and David Dobkin, " Matching 3D Models with Shape Distributions, “ ACM Transactions on Graphics 21，807-832 (2002)；以及由以下文献 定义的形状匹配的度量Serge Belongie, Jitendra Malik, and Jan Puzicha, “ Shape Matching and Object Recognition Using Shape Contexts, " IEEE Transactions on Pattern Analysis and Machine Intelligence 24，509-522 (2002)；所有上述文件都引入 作为参考。在本发明的一些实施例中，对于正被检查的像素i与检索像素j的相应特征值是 以由几何变换所改变的邻域来计算的，从一个邻域到另一个邻域。例如，两个邻域是处于不 同的相对方向和/或尺度。另外，其中一个邻域可以是相对于另一个邻域的镜像。例如，如 果找到像素i的特征值的算法使用了在该邻域的像素的灰度值，该像素与像素i在+X方 向有一特定距离，则像素j的特征值是采用替代的像素来计算的，该像素与像素j在+y方 向(90度旋转)有相同距离，或者在-χ方向(反射)有相同距离，或者在+χ方向(尺度改变)有两倍距离，或者在+y方向(旋转加尺度改变)有两倍距离，等等。可选地，旋转的和 /或缩放的邻域的灰度值是在计算该特征值之前插入的，尤其是如果旋转角度不是90度的 整数倍(在像素的笛卡尔网格的情形下)或者缩放因子不是整数。可选地，这些像素是安 排在一个三角形或者六角形网格内，或者更复杂的模格内。采用以这种方式定义的特征将是尤为有用的如果无需旋转和/或缩放和/或反 射而定义相同的特征，以及用于各种不同的旋转角度和/或缩放因子，相同的权重给出作 为结果的特征。这可导致一个距离测度，它取决于在图像中的结构的方向和/或尺度，至少 对于一些方向和/或尺度的范围。这将是有优势的如果该图像包括带有不同方向或尺度 的类似结构，或者图像的互相的镜像。优化特征的加权向量的遗传算法以上面所述的宽变化范围的可能特征，是很难系统地找到最佳的特征组以及相应 的权重α，该权重给出降低噪声的最好结果，对于特定的图像形态、特定的身体部分或疾病 状态、特定的对比材料的类型和浓度，或者其他图像特征。当不同特征能协同相互作用时， 是特别真实的，特定的一组特征是怎样有效工作的。图5显示了流程图500，显示了一种 遗传算法，它能被用于找到最佳的或至少较高的对于一组特征的加权向量。这里，“较高的 (superior)，，是指在特定的图像级别中对于降低噪声更有效的。需要明确的是，本领域已知的任何优化方法都可被用来找到加权向量。采用遗传 算法是有优势的，因为该组可能的特征具有大量的模块部件、个体特征，它们协同相互作 用，有点类似在有机体中的基因。在步骤502中，选择并排序一组特征。可选地，该组特征包含非常大量的特征，期 望得到遗传算法的结果，可被用于消除许多特征。以至于一些权重的找到率非常低，这将表 明相应的特征对于降低噪声并不是非常有用，而它们将从列表中清除，节省了运行降噪软 件的计算机时间。在步骤504中，选择噪声测试图像，以及等效的低噪声测试图像，以及更低噪声或 没有噪声的图像。将测试特定顺序的特征组和对应的加权向量的效力，通过将降噪算法应 用到噪声测试图像，计算输出的图像是如何接近于低噪声测试图像。如果该测试图像是从 与降噪软件能工作的图像相同的图像形态中获得时，遗传算法可运行得最好。图像的性质 和噪声的性质都可能影响特定的特征组和相应的加权向量的效力。可选地，该噪声测试图 像是在与低噪声测试图像在相同时间内获得的，但分别以不同的采集参数获得，例如，以低 X-射线剂量获得的CT图像，或者以较短的采集时间获得MRI图像。另外，噪声图像是通过 加入模拟噪声到低噪声图像而获得的。在步骤506中，创建加权向量的起始群域，例如采用随机数字产生器。每个加权向 量具有对于在顺序组中的每个特征的一个分量。可选地，该群域是足够大的，足以在可能的 加权向量的相空间内扩散，这是不可能的遗传算法将陷入局部最大效力，这是比全局最大 效力更坏。在另一方面，该群域是足够大的以致用于运行遗传算法的计算时间实际上不长。 以至于特别当测试大量特征时，很难确定需要使用的群域有多大，以及很难满意上述两种 情形。以至于可通过限制特征的数量来避免采用大量特征的计算困难，在群域的任意一个 成员中非零加权，假设合并多个不同特征的协同优势能获得，无需使用比在组中所有特征 的小片段更多的特征。验证这样一个假设，并获得一个想法什么构成足够的群域尺寸，可通过测试不同起始群域来完成，看最终结果是否非常不同。在步骤508中，评估每个加权向量的降噪效力，通过在噪声测试图像上运行降噪 算法，并将该结果与低噪声测试图像比较。加权向量的效力测量是，例如，在从降噪算法获 得图像与低噪声测试图像之间的均方误差。可使用均方误差的负数或倒数，如果想要定义 该效力测量，以致数字平均数越高，加权向量越有效。可选地，采用误差的视觉或生理学 测量，以测量该加权向量的效力；参见James L. Mannos and David J. Sakrison, “ The effect of a visual fidelity criterion on the encoding of image, " IEEE Trans, on Information Theory 20, no. 4, p.525 (July 1974)。在步骤520中，通过在群域中的所有加权向量中找到具有指示最大效力的效力测 量的加权向量，找到最大效力的加权向量。在步骤512中，确定遗传算法是否已经集中，或者是否已经找到效力足够的加权 向量。例如，通过比较最大效力的加权向量的效力测量与前代的数量，或者在前代的几个 量，看在每代中是否改变很大。另外，大多数有效加权向量的效力测量是与效力的一些绝对 标准进行比较。如果该遗传算法还未集中，则在步骤514创建新一代的加权向量。例如，通 过有更高效加权向量产生比低效加权向量更多的子孙，从而创建新一代，最低效力的加权 向量根本不会产生。可选地，另外，加权向量经历一些突变，不同加权向量的相同分量可不 同程度地交叉。采用有效的遗传算法可完成，也可采用本领域已知的其他遗传算法。在步骤508，评估新一代的每个加权向量，该程序继续用于后续的代，直到效力测 量开始聚集。一旦满意聚集标准，例如基于在对于过去几代的效力测量上的相对改变，流程 终止于步骤516，产生最大效力的加权向量作为输出。示例性特征组和加权向量表1显示了一组排序的特征，以及采用如图5所示的遗传算法找到的加权向量的 相应分量。特征列表不是一个全面的覆盖这里所述的不同的可能特征的列表，但包括了特 定类型的特征，当测试其他类型作为替代特征或者附加特征时，可导致更有效的降低噪声。 需要注意的是，在表1中大范围的顺序的权重幅度反映了特征值的相对幅度，它们是非标 准化的，不同特征的相对有效性也是如此，以致带有最高权重的特征是不需要作为最重要 的特征包括在内的。显示在图4A、图4B和图4C中的图像是全三维图像在χ-y平面的截面 图，采用表1所示的特征组合加权向量来降低它的噪声。全3D图像的体素在ζ-方向是比 在χ和y方向更长的。显示在图4A、图4B和图4C中的截面图在χ和y方向是约450乘500 像素。在表1中，“图像自身”是指像素的灰度值，以豪森菲尔德单位表示。不同顺序统计 特征的平均灰度值，以豪森菲尔德单位表示，特别是百分数，特别对于在邻域的灰度分布， 在邻域的中心带有检索像素。“均方误差”是在邻域内灰度分布的标准导数，以豪森菲尔德 单位表示。高斯过滤器特征涉及灰度响应的差分，以豪森菲尔德单位表示，带有不同参数的 两个高斯过滤器的中心像素，在毫米级测量。在涉及梯度的特征中，梯度单位是豪森菲尔德 单位每毫米。Ix、Iy和Iz分别是指梯度的χ分量、y分量和ζ分量，而梯度额是指梯度的幅 度。表1 特征与加权向量 诜择降低噪声时考虎的图像参数在大多数图像形态中，获得低噪声的图像是可能的，但不够经济或不够安全。例 如，在CT成像中，可通过采用较高X-射线剂量来降低噪声，但可能会增加成本，而高剂量的 X-射线会致癌。可以选择图像的噪声水平，例如，权衡以下方面要么获得用于诊断医疗问 题的较好图像，因而增加病人能被有效治疗的机会，要么增加病人患上癌症的机会，如果采 用高剂量的X-射线的话。对于MRI，可通过较高静磁场来获得较低噪声。已知高的静磁场对身体健康的不良 影响，高静磁场的MRI设备的造价也更高，只有少数医院能承受得起，也只有少数病人能获 得这样的待遇。全身治疗型MRI设备通常所用的最高磁场为约3特斯拉的孔磁场，而具有 更高磁场的设备将很难或不可能以已知技术的任何代价来建造，虽然小型设备有时可采用
19更高的磁场。带有1.5或2特斯拉的孔磁场的MRI设备是明显低廉的。在MRI图像中的噪 声也可通过采取更长的图像采集时间来降低，但该采集时间受到大多数病人所能承受的时 间的限制，而且增加采集时间也会降低病人通过昂贵MRI设备的机会，最终限制能接受MRI 检查的病人的数量。最后，包括CT和MRI，都可通过较低来获得较低噪声，但通常分辨率越低，噪声越 高，而高噪声会影响该图像用于医疗状况的诊断。图6显示了流程图600，描述了在CT、MRI或任意其他图像形态中选择图像参数的 方法，该方法补偿了采用如图3所述的方法来降低噪声的降噪软件的能力。流程图600所 示的方法并不限于采用如图3所述的降噪方法，而可采用本发明所述的其他方法。也可采 用现有的降噪方法，例如在图IA和图IB中描述的方法，或者目前已知或将来可知的其他降 噪方法。在步骤602，获得标准采集参数，例如CT的X-射线剂量和分辨率，或MRI的磁场、 采集时间和分辨率，都将产生想要的噪声水平。这些参数可以是特定成像设备的标准操作 指令的一部分，但它们也可取决于是否存在造影剂，取决于该图像是否被制作用于健康筛 选目的，较低质量的图像也可被证明，或用于检查已知的医疗问题，以及取决于诊断特定症 状所需的分辨率水平和噪声水平。在步骤604中，找到新的采集参数，例如采用检查表，在采用降噪软件后将给予相 同的噪声水平，至少在该图像的那些部分，噪音水平是诊断的关键。可创建这样一个检查 表，例如，通过模拟或通过以真图像进行测试，对于特定的图像形态，采用对降低噪声最有 用的特征组，作为标准采集参数的函数，或等同于作为想要的最终噪声水平的函数。该测试 或模拟也可以对于一个范围的分辨率来完成，或者它能简单地假设相对噪声水平与体素 体积或像素面积的反平方根成比例，这将是真的，如果在不同像素的噪声是无关联的。在步骤606，采用新参数获得一个图像，该图像具有比想要的噪声水平更多的噪 声，但将在CT图像或其他χ-射线图像中采用较低X-射线剂量，或者MRI图像中采用较低 磁场和/或较短的采集时间。在步骤608中，采用降噪软件降低在图像中的噪声，导致获得 具有相同的噪声水平的图像，或者用于医疗诊断的相同有用性的图像，与由无降噪的标准 采集参数来产生图像一样。采用残余图像图7显示了根据本发明的示意性实施例所述的流程图700，用于采用一个残余图 像，该残余图像介于原始图像与降噪图像之间，恢复了原始图像的一些非噪声部分，这些部 分本已由降噪软件去除。流程图700对于评估是有用的，通过反馈，改善了降低噪声的方 法。在步骤702中，获得原始噪声图像。在步骤704中，采用如图3所述的降噪方法从该图 像中去除噪声，或者采用这里所述的或现有技术中的任意其他降噪方法。在步骤706中，通 过从原始图像中减去降噪图像，产生残余图像。理论上，该残余图像将是纯噪声的。然而实 际上，在该残余图像中通常有一些低的对照结构，显示降低噪声不会仅去除噪声，还会去 除在原始图像中的一些结构，作为降低噪声的结果。需要注意，取决于图像形态，即使真正噪声也不同于在图像中点到点的振幅或其 他特征，导致在残余图像中可能大尺度的结构，即使仅从原始图像中去除噪声。为避免大尺 度结构从原始图像中不正确地去除，图7所示方法一次仅作用于图像的一部分，在该部分内的噪声是均一的，或者噪声水平是相对低的。在步骤708中，采用非线性保留边缘过滤器来平滑残余图像和定位边缘。该过 滤器降低在残余图像中的噪声，使边缘平滑为平行于该边缘的表面，从原始图像中相应 增强该残余图像所包含的任意结构。合适的过滤器包括非线性的各向异性流过滤器，例 如Beltrami流和连贯增强过滤器，例如，描述在以下文献J.Weickert，‘‘ Anisotropic diffusion processing, “ Teubner 1998, and J. ffeickert, " Coherence enhancing diffusion filtering, " International Journal of Computer Vision 31, no. 2-3, pp. 111-128(1999)。在步骤710中，“结构(structureness) ”参数是评价为在已过滤的残余图像中的 位置的函数，它具有较高的值，如果该图像在该点具有一致的结构。有许多参数可被用于这 个目的。例如，结构张量的特征值，或Hessian矩阵的特征值，都可被采用。在步骤712，已 过滤的残余图像的具有较高“结构”值的部分将恢复到降噪图像，以产生改良的降噪图像。在距离测度中采用相关矩阵图8显示在图像中用于降低噪声方法的流程图800，对于距离测度(!(F1, F2)采用 比上面所给出的不同的表达。在某些方面，新的表达式是上面给出的表达式的总结，它考虑 了在不同特征值之间的可能的关联性。我们认为，仅当β =2，上述表达式变成
_1] Ci(FxiF1) = [ ,(//-Z12)2 -Jff +... + ,(/；-fl)f对于(!(F1, F2)的新表达式还包括在不同特征之间的交叉项，考虑在不同特征值之 间的相关性。对于(KF1, F2)的新表达式的潜在优势是它保持不变，如果该组特征值是由 一组相同数量的新特征值所代替，每个新特征值被定义为旧特征值的独立的线性结合。因 为新的特征值组会包含在旧的特征值组中的所有信息，对于(KF1, F2)的新表达式在某些意 义上仅取决于在该特征值中的信息，而不取决于定义该特征的特别方法。对于(KF1, F2)的 新表达式的另一个潜在优势是加权向量因子α，以及交叉项因子，都取决于在图像中的 位置，但不需要是相同位置，反应了在该图像的不同位置上的噪声的不同水平，以及不同的 噪声关联性。例如，在一个CT图像中，噪声水平总体上不同于在图像的不同部分，由于通过 I-D或2-D映射来形成这种图像。在MRI图像中，噪声水平可以是不均勻的，取决于RF域 的局部密度，也取决于在静磁场中可能的多相形。从图像的较大噪声部分给予像素较小的 权重，采用类似像素的灰度值来估计真正的灰度值将是什么，对于特定的像素，将会没有噪 声。在图8中找到的距离测度取决于在不同特征值之间的相关矩阵。对于本例，它将 被显示了 特征值是个别像素在该像素的邻域内的灰度值，特征正被找到。在本例中，对于 相等的权重α，上述的CKF1, F2)的上述表达是在非局部平均数方法中采用的均方误差。在邻域内第m个像素和在邻域内第η个像素之间矩阵系数Σ 由下式表达Σ mn = E[(Im-y) (Ιη-μ)]这里Im和In都是在邻域内第m个像素和第η个像素的灰度值，μ是在该邻域的 所有像素的平均灰度值，而E表示期望值。可选地，通过测量在均一图像内的灰度水平，以 致在灰度水平与平均灰度水平之间的局部系数都只归因于噪声。将邻域的灰度值表达为向 量N，相关矩阵Σ ■可在向量表达法中表达为
21
Σ= Ε[(Ν-μ)Τ(Ν-μ)]这里，μ表示一个向量，该向量的元素都等于μ。在步骤802中，对于图像的不同区域，找到平均灰度值μ。在步骤804中，对于该 图像的不同区域，找到相关矩阵的系数，例如前面定义的Σ mn。在本例中，这些区域非常小 于整个图像，但非常大于一个像素。可选地，系数Σ 是假设为以位置平滑变动，并采用插 入来再不同区域之间评估它们。在像素i和检索像素j之间的距离测度是由下式给出的dE (Ni, Nj) = [(Ni-Nj)1 · Σ · (Ni-Nj)1/2这里，Ni和Nj都是向量，分别表示在围绕像素i和像素j的邻域内像素的灰度 值。注意，i和j并不是标记Ni和Nj元素的指数，但它们是标记两个像素的指数，理论上 这两个像素互相远离，基于它们的邻域的相似性而被定义，例如，像素i是在降噪算法中被 检查，而检索像素j。Σ 1是相关矩阵Σ的倒数。在步骤806中，计算Σ的倒数。可选地，当该邻域具有大量像素时，为减少计算时 间，该矩阵因而是非常大的，其倒数是由下式近似的Σ 1^ZnVnVn7
η这里，λη是第η个特征值的倒数，而Vn是Σ的相应的特征向量。如果仅采用一些 最大的特征值和它们对应的特征向量，Σ的倒数可被快速找到。接着，距离测量是由下式接 近dz(Ni, Nj) (ΣΛ [(Ni - Nj).Vn]2)1/2
η这里，总和是超过特征向量。这个距离测度类似于之前定义在像素i和j之间的 理论上距离测度，对于特征值Vn之一，带有每个定义为Ni-Vn的特征，也就是在邻域内像素 的灰度值的特殊线性结合，由相应的特征值入 给出权重α。估计像素的噪声水平如前所述，因为噪声水平通常在图像中的位置间变动，在噪声区的检索像素j贡 献了比像素i相对少的估计的真正灰度值，例如，通过给它们比在较低噪声区的检索像素 更低的权重Wj。可选地，给出区域的噪声水平是通过在上述同一仿真的图像内在接近的像 素之间的关联的测量在来估计的。另一个方法，不仅提供了给定像素的噪声幅度的估计，而 且提供了该噪声的记号的估计，将在这里描述。对于较高的噪声检索像素，每个像素的噪声 的估计并不仅用于降低权重Wj，采用图3所示的降低噪声的方法，但也可直接用于降低噪 声，通过从每个像素减去估计的噪声。最终得到的图像时由它自身处理为降噪图像，或者它 是结合的，例如，通过采取加权平均，由图3所示的方法找到已降噪的图像，或通过任意其 他降低噪声的方法。图9所示的流程图900描述估计体素的噪声的方法。该方法是基于这样的想法 在无噪声图像中，几乎所有像素的补丁(合适尺寸)可能仅占有该尺寸的补丁的全部向量 空间的非常小的片段，将每个补丁当作灰度值的向量。一旦限定了向量空间所占有的部分， 任何从该部分的偏差都会假设是由噪声引起的。在步骤902中，检查一组一个或多个低噪声测试图像，可选地类似于已被使用的
22图像，例如，采用相同成像形态，可选地带有类似的图像采集参数，和/或看身体的相同部 分，如果它们是医学图像的话。在图9所示的例子中，图像是3-D的，因此像素是指体素，但 仍可采用2-D图像的方法。在每个测试图像中，检查一组尺寸为N1XN2XN3的体素的补丁， 可选地包括在该测试图像中的所有组的带有这些尺寸的补丁。每个补丁可表示为在一个理 论上N1XN2XN3维度向量空间内灰度值的一个向量，它的尺寸表示一个N1XN2XN3补丁的 体素之一的可能的灰度值。这些图像的补丁被分类为在这个向量空间的簇，采用现有技术 中已知的任意聚类算法。簇的数量是，例如，20至100，或者100至500，或者500至2000， 或者低于20，或者高于2000。采用大量的簇需要检查更多的补丁以获得关于向量空间是如 何分为簇的好的统计表，因而需要更多计算时间来对补丁进行分类，但采用过少的簇会导 致这些簇覆盖太多的向量空间，对于降低噪声是无用的。本发明人已经发现，采用数百个簇 可做得很好，当补丁尺寸在5至10体素，或10至15体素，或低于5体素，或高于15体素。 特别地，11x11x11体素的补丁表现为非常好。虽然在本例中给出的上述补丁是在矩形阵列的形状内，但这些补丁并不需要这个 形状，也可以有任意定义的关于体素的其他形状。向量空间的维度是在补丁内体素的数量。在步骤904中，检查每个簇以找到有限组的独立参数，明显少于N1XN2XN3需要覆 盖全向量空间的基向量的线性系数，这相当好地近似于覆盖在该簇内的补丁。例如，在该簇 内的每个补丁，或者最大量的补丁，是非常近似为灰度值的一个向量，该向量是固定用于在 该簇内的所有补丁，加上有限组的基向量的线性结合，带有不同于在该簇内的不同补丁的 系数。另外，在该簇内的大多数补丁是有灰度值的固定向量加上有限数量的独立参数的一 个函数(并不必需线性函数)来近似的。基向量系数的数量或者其他独立参数，对于所有簇 不需要是相同的。基向量或者独立参数的数量(例如，用于给出的簇)通常小于&乂队乂队 的三倍或小于N1XN2XN3的一半或小于N1XN2XN3的四分之一或小于N1XN2XN3的八分之 一，或者大于N1XN2XN3W三倍。一些簇可具有比典型数量更多基矢量或更多独立参数。目 标是识别相位空间的小片段，该片段提供了对在该簇内大多数补丁的很好的近似。表示补 丁向量的非线性方法包括但不限于核PCA、Laplacian eigenmaps和扩散图。在步骤906中，获得噪声降低的图像。一旦补丁的簇和基函数或者参数组已经在 步骤902和904中找到，不得不再次找到它们以用于处理附加的图像，但该程序可在步骤 906中开始。在步骤908中，图像体素i开始一个循环，初始体素i的指数设等于1。在步骤910 中检查围绕第一体素的补丁(相对于它们的体素质疑来定义这些补丁)，并决定该补丁归 属于哪个簇。另外，也可采用图3所示的方法，一起检查群体的体素，取代每次检查一个体 素的方式，而这些补丁是相对于它们的体素的群来定义的，或者相对于插入在体素之间的 虚拟体素来定义。如果这些补丁挑选得合适，则测试图像是足以代表正在处理的图像，这些 测试图像是没有许多噪声的，但充满了定义好的结构，而正被处理的图像也不由噪声支配， 则正被处理的图像的大多数补丁可干净地落入一个簇或另一个镞。如果一个或多个这些情 形不够满意，该方法仍是有用的，只是在降低噪声方面的效力低。在步骤912中，找到该簇的对于基函数的系数或其他独立参数，提供了最好的近 似，或至少理论上最好的近似，到体素i的补丁。对于大多数补丁的近似越好，降低噪声的 效果越好，因为本方法假设任何在观察的灰度值与它的近似之间的差是仅由于噪声引起的。在步骤914中，估计体素i的灰度值的噪声分量，通过考虑在体素i的实际灰度值与体 素i将有的灰度值之间的差，采用根据簇的固定向量加上有限组的基函数(或有限组的独 立参数的函数)的线性结合的近似。这个程序的基本原理是在缺乏噪声的情况下，补丁可 由近似来很好地描述，在测试图像中找到的在该簇内的所有或大多数补丁都运作良好，因 此任意差是很可能主要由于噪声引起的。在步骤916中，体素i被增加1。在步骤918中，确定在该图像中的所有体素是否 都已经被检查过。如否，在步骤910开始下一个体素的检查。如果所有体素都已经被检查 过，则在步骤920结束该流程，输出估计的每个体素的噪声。可选地，估计的噪声被用于降 低该图像的噪声，要么直接从该图像中减去所估计的噪声，要么间接地例如通过估计的噪 声水平影响权重Wj，当采用图3所示的降噪方法，或通过采用两种方法的结合。这里所用的术语“约”是指士 10%。术语“包括”、“包含”、“具有”以及它们的连接词是指“包括但不限于”。这个词组
包含术语“包括有”和“主要含有”。词“主要含有”意指其成分或方法可包括附加成分和/或步骤，但这些附加成分和 /或步骤不会从本质上改变所要求的成分或方法的基础特征和创新特征。这里所用的“一个”和“所述的”包括复数，除非另有明确说明。例如“一个化合物” 或“至少一个化合物”可包括复数个化合物，包括它们的混合物。这里所用的词“示范的(exemplary)”是指“用作一个例子或例证”。任何描述为 “示范的”实施例并不是必须作为比其他实施例更好或更有优势的例子，和/或排除从其他 实施例中引入的特征。这里所用的词“可选地(optionally) ”是指“以一些实施例来提供，不以其他实施 例来提供”。本发明的任何特定实施例可包括多个“可选的”特征，除非这些特征发生冲突。通过本申请，本发明的不同实施例可在一个范围内展示。需要明确的是，在该范围 内的描述仅是为了方便和简短，而不应视为对本发明的范围固定的限制。相应地，描述的范 围应当被视为具有特别揭示了所有可能的子范围以及在该范围内的个别数值。例如，例如 从1至6的范围的描述应当被视为也揭示了例如1至3、1至4、1至5、2至4、2至6、3至6 等子范围，以及在该范围内的单个数字，例如，1、2、3、4、5和6。不论范围有多宽，都可这样 来应用。无论一个数字范围在这里如何指示，它是指包括在该指定范围内任意引用的数字 (部分或整体)。这里所用的词组“归类”或在第一指示数字与第二指示数值“之间归类”以 及将第一指示数字“归类”到，或从第一指示数字“归类”到第二指示数字是可交换使用的， 它们的意思是指包括第一和第二指示数字以及在两个数字之间的部分或全部数字。需要明确的是，为表述清晰，本发明的特定特征是描述在分开的实施例中的，但也 可在单个实施例中结合提供。相反，为表述清晰，本发明的不同特征是描述在单个实施例中 的，但也可分开提供或者在任意合适的子组合中提供，或者在本发明的任意其他描述的实 施例中以合适的方式提供。描述在不同实施例中的特定特征不应视为那些实施例的必要特 征，除非该实施例缺乏那些特征元素就无法操作。虽然本发明已经结合特定实施例进行了描述，对于本领域技术人员显然还有许多 变动、修饰和改变。相应地，所有这些变动、修饰和改变都将落入本发明的精神与所附的权利要求的范围之内。 在本说明书中提及的所有出版物、专利和专利申请都在这里整体引入作为参考， 即使每个个别的出版物、专利或专利申请时特殊的和单独指示为引入作为参考。另外，在本 申请中的任意参考文献的引用或鉴定都不应认为该参考文献是本发明的现有技术。以至于 所用的章节标题，它们也不应视为必要的限制。
权利要求
一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻域的对应特征的值；以及c)基于图像的检索像素的原始灰度值或者变换灰度值，计算对于正被检验的每个像素的降低的噪声灰度值，相对于检索像素的一个或多个特征值，正被检验的像素具有更大敏感性的类似特征值；其中，所述计算至少一个特征值的步骤包括计算在邻域的像素的原始灰度值或变换灰度值的特征性分配，除了计算在该邻域的所有像素的平均灰度值之外。
2.根据权利要求1所述的方法，其特征在于，所述计算特征值的步骤包括计算第二矩 分布或更高矩分布。
3.根据权利要求1所述的方法，其特征在于所述在邻域的像素的灰度值是对于那些 像素的图像的原始灰度值。
4.根据权利要求1所述的方法，其特征在于所述在邻域的像素的灰度值是由过滤器 变换的图像的灰度值
5.根据权利要求4所述的方法，其特征在于所述过滤器是高斯过滤器。
6.根据权利要求1所述的方法，其特征在于，所述计算特征值的步骤还包括计算所述 像素在第二邻域的灰度值的至少一个特征性分布，该第二邻域是所述邻域的适当的子集。
7.一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻域的 对应特征的值；c)计算对于正被检验的每个像素的降低的噪声灰度值，基于它的检索像素的原始灰度 值或者变换灰度值，检索像素的一个或多个特征值比正被检验的像素的类似特征值具有更 大的敏感性；其中，所述计算至少一个特征的值的步骤包括评估在邻域的图像的变换，或者评估在 邻域的图像上的过滤器的效果，以致优先选择在邻域的尺寸和一些像素尺寸之间的中间尺 度范围的结构，或者一些方向而不是其他方向的结构，或者兼有这两种结构。
8.根据权利要求7所述的方法，其特征在于，所述计算特征值的步骤包括评估线形过 滤器对于所述图像的响应。
9.根据权利要求8所述的方法，其特征在于所述线性过滤器是一个小波滤波器。
10.根据权利要求8所述的方法，其特征在于所述线性过滤器是一个方向过滤器。
11.根据权利要求10所述的方法，其特征在于，所述计算特征值的步骤包括评估在超 过一个方向上的过滤器的效果。
12.根据权利要求7所述的方法，其特征在于，所述计算特征值的步骤包括评估在所 述邻域的图像的非线性形态变换。
13.根据权利要求12所述的方法，其特征在于所述的形态变换是多尺度变换。
14.根据权利要求12所述的方法，其特征在于所述的形态变换是形态算子。
15.根据权利要求7所述的方法，其特征在于在所述邻域的方向的角上所述特征的值的相关性具有对称性，以致如果所述邻域是旋转至少一个角度，对于在该邻域的图像像素 的任意灰度值，所述特征值是不变的。
16.根据权利要求7所述的方法，其特征在于对于至少一个特征，正被检验的像素的 特征值与检索像素的对应特征值是以它们在不同相对方向或尺度或两者兼有的邻域来评 估的。
17.根据权利要求7所述的方法，其特征在于所述减少的噪声灰度值是基于所述检索 像素的加权平均灰度值。
18.根据权利要求1所述的方法，其特征在于所述加权平均值采用的权重取决于在检 索像素的特征值与正被检验的像素的对应的特征值之间的差的测算。
19.根据权利要求1所述的方法，还包括a)识别残留图像的相对更多结构部分，表示从原始图像到减噪图像的变化；以及b)存储所述结构部分到减噪图像，以产生一个改善的减噪图像。
20.根据权利要求1所述的方法，还包括估计一个或多个检索像素的噪声水平，其中， 所降低的噪声灰度值是低相关性或不相关于具有更高的估计噪声水平的检索像素。
21.根据权利要求1所述的方法，还包括a)在图像中分配像素的一个补丁到在这些补丁的向量空间内的复数个簇之一，该补丁 表示为灰度的向量；b)近似所述补丁作为所降低的一组非相关参数的函数，参数量少于在所述补丁内的像 素量，与所述簇关联；以及c)根据在所述补丁的实际灰度与作为所降低的一组非相关参数的函数的所述补丁的 近似之间的差，校正所述补丁的像素的灰度。
22.根据权利要求21所述的方法，其特征在于所述函数是所降低的一组基向量的线 性组合，而非相关参数是在该线性组合内的基向量的系数。
23.一种获得带有改变的图像采集参数的医学图像的方法，包括a)确定一组一个或多个图像采集参数，这些参数会在图像降噪后达到想要的噪声水平；b)获得带有所述采集参数的图像；以及c)根据前述任意权利要求所述的方法降低该图像的噪声。
24.根据权利要求23所述的方法，其特征在于所述获得图像的步骤包括获得CT图像；以及所述确定一组一个或多个图像采集参数的步骤包括确定一组与该图像未降低噪声时 的情况相比具有降低的χ-射线剂量的参数。
25.根据权利要求23所述的方法，其特征在于所述获得图像的步骤包括获得MRI图像；以及所述确定一组一个或多个图像采集参数的步骤包括确定一组与该图像未降低噪声时 的情况相比具有降低的静磁场或减少的采集时间的参数。
全文摘要
一种降低图像噪声的方法，包括a)对于在图像中正被检验的每个像素，选择一组检索像素；b)计算每个检索像素的邻域的一个或多个特征的值，并计算正被检验的像素的邻域的对应特征的值；以及c)基于图像的检索像素的原始灰度值或者变换灰度值，计算对于正被检验的每个像素的降低的噪声灰度值，相对于检索像素的一个或多个特征值，正被检验的像素具有更大敏感性的类似特征值；其中，所述计算至少一个特征值的步骤包括计算在邻域的像素的原始灰度值或变换灰度值的特征性分配，除了计算在该邻域的所有像素的平均灰度值之外。
文档编号H04N1/38GK101919230SQ200880123648
公开日2010年12月15日 申请日期2008年12月25日 优先权日2007年12月25日
发明者伊利兰·达昂, 兹维·德维尔, 埃泽尔·巴拉维夫, 塔尔·克尼希, 盖伊·罗斯曼 申请人:梅迪奇视觉-脑科技有限公司
最新回复(0)