CN1251147C

CN1251147C - 无失真的图象对比度增强的方法

Info

Publication number: CN1251147C
Application number: CNB018124224A
Authority: CN
Inventors: 布莱恩·G·詹姆斯
Original assignee: ASENTEC TECHNOLOGY Co Ltd
Current assignee: ASENTEC TECHNOLOGY Co Ltd
Priority date: 2000-07-07
Filing date: 2001-07-05
Publication date: 2006-04-12
Anticipated expiration: 2021-07-05
Also published as: US6633684B1; EP1299858B1; ATE500569T1; MXPA02012756A; WO2002005206A2; JP4175461B2; AU2001272251B2; AU7225101A; NZ523251A; EP1299858A2; KR20030012934A; KR100572444B1; JP2004503029A; WO2002005206A3; DE60144150D1; CA2352678A1; CN1440542A

Abstract

通过图象的点强度的对比度的局部增强可以增强如X光这样的图象。用第一低频上曲线拟合局部最大值，用第二独立的低频下曲线拟合局部最小值，并以位于其间的原始图象数据形成一个航路。为每点提取出航路局部最大强度和航路局部最小强度之间的局部范围。用航路的局部范围和图象的动态范围之间的比来换算每点，以便使在它和它相邻点之间的强度的变化最大。最好使用一个迭代移动平均技术来建立航路曲线。在一个最佳实施例中，将换算到航路之外的非正常值点以比动态范围高的精度临时存储。形成航路校正数据的直方图，具有大于动态范围的范围，并基本上包围所有非正常值点。在这个直方图校正中只剪裁最不正常的非正常值，并将所得的对于整个图象的范围界限换算到动态范围。

Description

无失真的图象对比度增强的方法

发明领域

本发明涉及用于提高数字X-光图象中的相邻数据之间的对比度的方法。更特别地，该方法确定相邻数字图象数据之间的实际对比度，并在不失真图象的情况下将每点的强度扩展到图象的动态范围。

发明背景

X-光图象典型地出自采用胶片照相术的标准X-光机器。在这些情况下，利用数字扫描技术将所得的X光图象转换成一个计算机文件。最近，出现了使用一组光敏传感器来直接捕获X光图象的数字版本的X光机器。X光图象被用作医学诊断工具。虽然有其他高级的相关成像处理，特别是CAT扫描和MRI-s，但因为X光图象相对便宜，所以它仍然被广泛使用，并且构成这种图象的大多数，并且非常有可能继续保持这种状态。本发明提高了X光图象对医生的用处。

部分地由于涉及X光成像的实际限制，难以提供既定义了类似肺的相邻软组织内的密度的变化又定义了类似骨的相邻致密组织内的变化的图象。变化由强度的改变来表明。每个数字图象与一个动态范围相关联。一旦被显影(对于胶片)或数字化地呈现为一个正象，致密物质的亮强度(通常描绘成白色)区域典型地占据动态范围的高端，低强度(黑色)占据动态范围的较低端。虽然成像方法确实可以捕获细微的变化，但这种变化之间的强度不能被人眼容易地检测到。当传统上在一个光盒中透照出的胶片图象被转换成数字图象并在一个数字显示器上显示时，这种情况则更糟。换句话说，黑区域内的变化和白区域内的变化不容易被区分。

因此，已知有很多用于提高数字X光图象中的对比度的方法，其中最著名的是对比度扩展。用于实现对比度扩展的手段的各种方法是几个专利申请的主题。

例如，在颁给Maack等的美国专利5,357,549中，提供了一种用于只在感兴趣的特定区域——例如胸部X光的肺部区域——中扩展图象强度的技术。Maack将这个称为动态范围压缩。Maack定位低频分量，确定均衡因子，并将其运用到图象来用于压缩低频分量；从而使得动态范围的剩余部分可用于图象强度的较高频率区域。这个方法不能增强已经选择的并且是诊断者立即感兴趣的多于一个的图象强度区域，损失了其他区域的数据。

颁给Fang的美国专利5,835,618在Maack基础上进行改进，使用一种用于同时增强在暗和亮强度区域中的图象的动态范围再映射的方法。这个再映射或校正技术实际上是平滑数据(例如通过低通滤波器)，确定数据中值，将已平滑的数据调整到中值，然后将该已平滑的已调整数据运用到原始数据。提供了已调整数据的两个曲线，每个曲线从中值偏移，一个向上偏移常数Δ1，一个向下偏移常数Δ2，用于建立上限和下限，并将数据分成两个不同的总体。常数Δ1、Δ2定义了一个范围。然后，将单独的加法或乘法算法首先运用到该范围内的原始数据上，其次运用到该范围外的数据上。例如，在加法选择中，以数据的原始和中值强度的差来调整该范围内的原始数据，这个差由用户定义的在0和1之间的控制参数来换算。然后，将已调整数据换算到图象的整个动态范围。采用一个不同的算法来调整该范围之外的数据。

不幸的是，Fang的方法对范围之内和之外的数据处理不同，产生了人为细节和其他失真结果，例如在范围边界产生了不相继。刚好在范围内的数据和刚好在范围外的数据可以导致非常不同的值，使较强信号失真。只要平滑图象从整个图象的中值偏离，调整算法就导致这些失真，偏离的幅度影响失真的幅度。Fang认识到，用户可以操纵动态压缩的程度。然而，为了使失真最小，用户必须操纵下面的每一个：校正的衰减、向上偏移和向下偏移。例如，所选择的范围越大，则失真被最小化得越多，并且细微细节损失得越多。较小范围可以增强较弱信号，但较强信号变得严重失真。也就是说，该过程在用户方面需要时间、努力和经验来管理三个参数来试图最小化图象失真同时最大化图象增强。

上述图象增强的每种方法导致原始数据的损耗或失真。损耗和失真呈现出假象，可以严重地损害一个放射学家或其他诊断者对图象的解释。理想地，如果强度变化在图象中的相邻数据之间存在，则在不引入假象的情况下它们之间的对比度应该被最大地增强以供诊断者检测，并且不管强度变化是在图象的亮区域还是暗区域中。

一种最优方法应该采用与它们性质无关的控制参数。相关变量需要用户作出折衷，牺牲一个结果的一部分以便实现另一个结果的一部分。不幸的是，在能够正确地实现之前，这种判断需要用户获得该技术背景的专业知识。此外，这样的控制参数需要是坚固的，使得一给定参数中的小改变导致可管理的改变，不会引起过分的结果，并且即使很差的选择也会给出“可忍受的”结果。

发明概述

本发明是一个在不使图象失真和不增加人为细节的情况下使图象的相邻点之间的细节最大化的方法。局部地，图象中的所有对比度都以相同方式被校正。

在优选实施例中，将该方法运用到相邻点之间的变化在数字上是显著的但对比度对于人眼来说太低而不能辨别的图象上。这种情况包括X光图象。在X光图象中有比一个人能猜到的更多的信息。本发明增强了这种信息，以便得到一个更有展示作用的处理后X光图象。本发明的成就是：

·提高了的对比度——用于帮助医生和其他诊断者的可视解释；

·处理后图象和输入图象的良好相关性——例如，一个胸部图象的处理后图象仍然类似原始胸部图象；

·使细节最大化——X光图象的每个部分均以最大可能可见性被最大化，在暗区域内和亮区域内的细微变化之间对比度提高。

·图象增强的自动化——这样它不要求终端用户具有在先的图象增强经验；以及

·避免失真和人为细节——这样避免了向图象解释这种已经有挑战性的任务增加错误信息。

在本发明的主要方面，通过使第一低频上部曲线或表面拟合于局部最大值，并使第二独立的低频下部曲线或表面拟合于局部最小值，其间的空间或体积形成一个航路(fairway)，通过相邻点强度之间的对比度的局部增强来增强图象。原始图象数据驻留在航路中。为每一点提取出航路局部最大强度和航路局部最小强度之间的局部范围。一个局部比例因数确定为该图象的局部范围和动态范围之间的比。然后将每点用其局部比例因数换算，以便使它和它相邻点之间的强度上的变化最大化。

更可取地，采用两维移动平均来确定低频曲线，来调节相邻列和相邻行之间的点强度的变化。

更可取地，甚至落在航路之外的非正常值也被增强，不是通过失真截断处理增强，而是通过根据非正常值的强度范围的确定来重新换算增强图象的直方图校正而增强。更具体地，局部换算的非正常值具有在图象的动态范围之外的强度，因此被以较高精度临时存储或保留。

接着，所有存储为高精度强度的数据，包括非正常值，都被放在一个具有大于图象的动态范围并且大到足够充分捕获几乎所有非正常值的预定或扩展范围的直方图中。进行直方图计数，并将预定剪裁率运用到直方图低端和高端。因为航路趋向在相邻点之间的强度中的大变化和小变化，非正常值可以出现在图象中航路之外几乎任何局部位置。其结果是，以这种方式选择的非正常值的最不正常值通常被广泛地分散，因此在剪裁时可最小化地影响图象增强。剪裁的图象点具有一个定义具有最小强度和新的最大强度的新范围的剪裁范围。对所有点进行第二次换算，这次是以确定为剪裁范围与图象的动态范围之比的比例因数来换算。

更可取地，采用两维移动平均来确定低频曲线，并且，更可取地，采用滤波框中值确定来确定低频曲线，这通过取在前的滤波框总和并仅仅减去落后点强度和加上超前点强度以获得一个新的滤波框总和，显著地减少了所需要的计算数量。仅仅通过用框中的点数除框总和来规格化每个框总和来完成移动平均。

附图简要说明

图1是一个胸部的X光，其中，许多细节是立即可见的，而其他细节在骨盆和脊柱的表面上均匀的亮区域中以及肺的均匀暗血管区域中看起来是丢失了；

图2是在应用了本发明的优选方法之后图1所示X光的数字增强图象；

图3显示了来自图1的X光的图象的一个样本列数据的数学选择的想象表示，以供使用本发明的方法处理。样本列被显示为从图1的图象提取出并放在一边，所提取的列以放大的宽度显示以便容易查看；

图4a和4b是相关图，显示了图3的样本列数据内的数据强度的变化。在图4a中，曲线图的从左到右对应于图象从下到上，图4b是对应的样本列的复制，还以放大的宽度显示以供查看；

图5a和5b是运用到图4a的数据上的简单的现有技术形式的对比度扩展的图。图5a显示了其上绘有一个平滑的低频中值的图4a的数据，图5b显示了根据自中值的局部数据变化重新换算的输出数据，该输出数据然后被扩展到全动态范围；

图6a-6c是相关图，是将本发明的方法运用到图4a的数据上的显示。图6a显示了航路的上下边界曲线。图6b显示了对应的输出数据。图6c是对应的图象列的复制，其以放大的宽度显示以供查看；

图7a和7b是其上绘有换算航路的图4的原始输入数据的图；图7a的航路是理想的手画草图；图7b的航路是采用本发明的优选移动平均处理计算出的；

图7c是其间带有原始数据的2D边界表面的想象和部分表示；

图8a和8b是分别从绘出低频上下边界曲线或平滑趋向函数产生的孤立图象；

图9a和9b是显示和比较原始输入数据和以比动态范围高的精度存储的第一航路校正换算后的数据的图；

图10a和10b是显示在确定上趋向函数中的初始步骤的图。图10a显示了叠加在原始输入数据上的低频中值，图10b显示了在中值上或大于中值的残留点；

图11a和11b显示了从图10b的数据的迭代收敛确定，其中，为在运用前面的中值之后剩余的残留数据迭代确定一个新的中值，以便将数据总体减少到一个新的更小的残留点数；

图12是显示在图象边缘反照的数据的图；

图13a和13f显示了分别用于将一个有效的移动平均技术运用到图象数据的伪代码和变量定义；

图13b-13e是用于提高求和的效率的移动平均技术的想象例示，图13b用于以一列递增的一个7点行的列求和；图13c使用一个7×7框进行逐列的行求和；以及，图13d和13e使用一个2D框，先是逐行的列求和，然后以一列递增；

图14a-14c显示了伪代码，图14d显示了用于示范运用图13a的移动平均数据并确定上下边界曲线的迭代方法的一个实施例的变量定义；

图15a显示了用于一个替换图14b的技术的改进的收敛技术的伪代码；

图15b是显示一个依据图15a的伪代码对图4a的数据的上边界曲线的更快收敛的迭代确定；

图16是显示采用图15a的更快收敛技术获得的航路的所得的上下边界曲线的图；

图17a和17b是相关图，更特别地，图17a显示了依据图16的航路的上下边界曲线，图17b显示了在第二直方图校正步骤之后的相应输出数据；

图18a和18b显示了伪代码，图18c显示了用于运用根据图13a、14a-14d确定的上下边界曲线和应用用于输出如图2的图象的校正图象数据的第一航路和第二直方图校正的变量定义；

图19、20a和20b显示了从噪声抑制的程度的变化产生的输出图象。具体地，

图19运用1.5％的最小边界间隔(surf_mdff)，

图20a和20b运用3％的边界曲线的最小间隔，并且是与图2相同的图象，是为了方便比较而再现的，以及

图21运用12％的最小边界曲线间隔；以及

图22-25显示了在运用到所示范的移动平均技术上时从滤波(imag_pcnt)程度的变化产生的输出图象。具体地，图22运用图象大小的5％的滤波器并且还是与图2相同的图象，是为了方便比较而再现的，图23运用7％的滤波器，图24运用10％的滤波器，以及，图25运用图象大小的20％的滤波器。

优选实施例的详细描述

图象特性

每个数字设备具有一个动态范围，该动态范围是对其记录相对能量波动的能力的度量。这里，在X光的情况下，荧光的强度从荧光屏直接发出并由一个电荷耦合(CCD)器件或扫描仪中的薄膜透视来捕获。这个动态范围通常由记录或扫描分辨率和数字存储装置(数据文件)来设置。一个典型的数字文件使用给出动态范围为0-255的8位存储器，。另外已知的还有使用给出动态范围为0-1024的10位存储器的文件、给出范围为4096的12位存储器的文件以及给出范围为0-65,535的16位存储器的文件。对于这个应用的目的，申请人使用一个8位存储器(0-255)描述了本发明的过程，该过程与实际的图象动态范围无关，所公开的范围不作为对本发明的限制。

本发明的一个目的是运用图象增强技术来使局部对比度最大化而同时避免与超过图象的动态范围相关联的失真。本发明认识到，使局部对比度最大化的一部分是将局部范围扩展到动态范围。从点到邻近点的相对对比度被保留在局部水平——所处理图象的每个有效尺寸的子部分必须利用所使用的技术的大多数动态范围。增强技术的另一部分是确定要扩展的局部范围。另一个部分是最小化或防止与在动态范围的极限值削波相关联的失真。削波发生在对比度扩展将一个强度值换算到最小值或最大值之外的一个理论值的地方。例如，对于动态范围0-255，在假定范围256-300内的换算后输出值只能被记录为255，因此这个区域内的细节丢失。

参考图1，提供了用于图象增强的图象10。运用本发明产生图2所示的输出图象30。

样本图象

图象10包括大约878×1007象素的数据阵列，每个象素以0-255强度级的样本动态分辨率存储。本公开文本的大部分是将本发明的应用例示为它只以一维形式运用到一个878数据列的样本上。如同这里也描述的，该方法实际上采用用于分析行和列的嵌套操作在两维情况下实施。

下面介绍一个一维中间处理列的例子。

参考图3，将从从左开始在图象大约2/3处选择出的一列数据11显示为在物理上是从图象10提取出的(只是图解目的)。

转到图4a，样本列11的原始数据绘成高频曲线14。图的左边对应于X光图象10的下面。横轴对应于从图象10的下面开始计数的点数的索引(0-1007)。纵轴是以动态范围的一个分数表示的每点的相对强度。对于图象10，值1.0表示尽可能地白(对于这个特定文件，是动态范围的上限255)，0表示尽可能地黑(动态范围的下限)。虽然在这个图象中使用黑到白的整个动态范围(0到1)，但这个特定列或数据子集只使用大约75％的动态范围。例如，在图象10的脊骨处的数据的相对强度是接近1.0。

动态范围的较差应用的例子

在图4b中，显示的样本列数据11对应于原始数据曲线14的方向。

在图5a和5b中，参考一个现有技术的结果，在现有技术中已知通过从图4a的原始数据11移去低频分量并将所得的范围扩展到全动态范围来保持局部对比度。在图5a中，将一个计算出的平滑低频曲线或趋向15叠加在真实数据14上，其中该曲线或趋向15是采用现有技术从真实数据14计算出的。

在图5b所示的这个现有技术中，输出曲线16是从原始数据14减去趋向15后，再将数据的最小值和最大值换算到动态范围的结果。这个简单的方法保留了局部细节，然而，图5b的残留数据图在动态范围是怎样局部使用方面是不够的。注意，在点330至350具有非常小的峰谷范围的在17处的数据、在580至620的在18处的数据以及从780横跨到850的在19的数据将被理想地校正到动态范围的相似使用。注意，在18处的数据的范围大约只使用动态范围的五分之一，而在17处的数据使用大约一半，在19处的数据使用几乎全部动态范围。这个方法对所有数据集合14的处理相同，因此没有使X光图象内的细节的可见性最大化。参考现有技术图5b，可以看到，由于在19处的大范围强度数据，使得在18处的数据的小强度范围只能被扩展到动态范围的大约25％。

改进应用

为了最优地扩展图象10的所有区域，甚至那些初始状态利用了动态范围的一小部分或利用了一大部分的区域，局部对比度必须被独立地处理。完成这一点的一种方式是识别并选择相邻数据的子部分，并充分独立于数据14的剩余部分来对其换算。

参考图7a，显示了三条曲线，图4a的输入或原始数据14，以及大略手画的上曲线20和下曲线21。曲线20、21代表一个提供原始数据14的理想边界的概念。这些边界曲线20、21具有低频率，并且更进一步地具有如下特性：它们相切于原始数据的局部极值；上边界相切于局部最大值(20_i，20_i+1…20_n)，下边界相切于局部最小值(21_i，21_i+1…21_n)。

两条曲线20、21的低频方面确保局部相对强度或局部对比度被保留，就象在图5a中的单个趋向曲线15一样。然而，与现有技术图5b及颁给Fang的美国专利5,835,618的现有技术系统不同，边界曲线20、21独立响应于实际数据14，而不是在中值15上。

上下边界曲线20、21一起形成一个边界航路22，该边界航路22可以类似于一个将最终被扩展到动态范围的限定数据范围。

现在可以看到，与现有技术图5b的仅仅25％不同，在580-620点之间的强度范围数据18和在780-850点之间的数据19都可以直接换算到动态范围的大约100％。

本发明的挑战在于确定这些上下边界曲线20、21要被放在哪里以及如何放、以及克服由计算机处理系统所带来的复杂局限性。说明书的剩余部分公开了用于确定这些边界曲线20、21的方法。

理想情况

如等式1所引入的数学校正(即图18a中的R)，运用边界曲线20、21调整每个原始输入数据点，每次调整一个。一个输入数据点是一个具有已知强度的象素或点，对于选自图3的列11，在每列11中从样本图象1007点或行选取。

……………等式1

其中：

·图象点是在一给定点的输入数据的值；

·下边界是在同一点的最小曲线的值；

·上边界是在同一点的最大曲线的值；

·动态范围界限是该文件或系统的动态范围界限(这里是256)；

以及

·输出数据是校正后的数据值。

这个第一校正的成功取决于确定上下边界曲线20、21的方法，该两曲线是由从原始输入图象数据14和两个预先确定的变量或参数确定的。

如果可以实现理想的手画情况，则没有数据点落在边界曲线20、21之外，等式1是一个全解，所有局部对比度被保留并与图象的其他区域实质无关地被换算。不幸的是，数字方法的局限性和收益递减的定律均使边界曲线将不是必然包围所有点。

参考图7b，计算出的曲线20、21与前面手眼估计的最佳上下边界曲线具有相似性。然而，不同于理想情况，有一些差别：即，在这个图中有一些落在计算出的航路22之外的非正常值点23；换句话说，虽然采用了这里所公开的最佳的另外的技术，使影响变得基本上不明显，但原始数据14中的点和图象10整体仍然稍微被削波。

为得到最佳结果，实施一个两步骤校正来处理这些计算出的航路22的实用性。采用包括使用直方图的另外的技术来处理非正常值，并避免其中包含的图象信息的丢失和失真。

例如，在两步骤校正之后，对于图3、4a的样本列，上曲线20的整个长度产生大约只占列数据14的总体的0.2％的削波数据的故障率。削波数据都成为最不重要的点，因为他们在整个校正图象30中分布，其结果是没有任何原始图象10的相邻点集被失真。因此，即使在第一步骤中边界曲线20、21本身可能有大约10％的故障率，但第二和最终校正将产生低得多的仅约0.2％的故障率。

第一航路校正

如果落在计算出的航路22之外的非正常值点23(例如，参见图9a，大约在点330和350处)被立即约束到动态范围，则图象数据完整性受到损失。这样的约束条件出现在使用计算机处理系统时，该系统利用声称与用于输入图象10具有相同限制精度的变量。例如，落在上边界之上的非正常值点假定在数学上换算到280。根据0-255的动态范围，这个非正常值点将被削波到255，相对强度数据将被丢失。因此，采用数字存储装置来避免在航路校正中的削波非正常值，可以使用一个高于动态范围的临时数据精度，并将数据保持在其中，直到运用一个另外的第二校正。

将等式1运用到图象10的原始数据，并采用具有高于图象动态范围界限的精度的变量类型来临时存储图象结果。例如，对于目前常用的用于所存储的图象数据的0-255(8位)、1023(10位)或4095(12位)的动态范围界限，一个标准的2字节(16位)的带符号整数变量就足够了。2字节的整数允许一个更大范围，到大约+/-32,000，并同样包括负数。因此，一个2字节的带符号整数可以保存非正常值而不产生任何失真。

参看图9b，将等式1运用到图3、4a的样本列数据，并将0-255的动态范围界限用作一个例子，所计算出的临时数据24的范围几乎为-50到350。除了被临时存储在更高精度的2字节整数变量中之外，所得的数据范围对于低于0和大于255的非正常值点将产生失真。

运用第二直方图校正，将图9b的这个航路校正的所得数据24适当地换算到实际动态范围。

第二直方图校正

图9b的临时高精度图象数据大部分被强制变回图象的动态范围。可以采用各种方式。

一种方法是采用一个强制故障率。简单地，计算临时数据24的直方图，并为了数据14的大多数而牺牲极少量级的非正常值点23。一个合理的第一途径是削波或剪裁大约1％；即，剪裁数据的底1％和顶1％，并将剩余点换算回图象的动态范围。

因此，如下面以及图18a中更详细公开的，第二校正使用直方图来集合图9b的航路校正临时数据24，并剪裁数据的最低和最高1％。然后将较高精度的临时变量中的点值的剩余的剪裁后范围再次换算回较低精度的动态范围0到255。所得的输出数据34显示在图6b中，所有点的整个输出数据被表示在图2的输出图象30中。在样本列11的输出数据34上的影响被显示在图6b中，其中，对比度显著增大，并且，对于相同的动态范围，不考虑这些点在图象10的实质低强度(暗)区域或实质高强度(亮)区域中是否是初始相邻的。比较图6b的代表增强后数据34的列31与图4b的原始列11和数据14。

注意，在图6a中，航路的上下边界曲线20、21定义了都近似相等的局部动态范围，即使它们包围了在点330到350出现的强度数据17的中间范围、在580到620的数据18处的小范围和在780到850的数据19处的最大局部范围之一中的每一个。在应用中，发现剪裁率远小于1％。

这个只强制剪裁存储在临时高精度数据中的少量最高和最低值的方法意味着在数据点的列中只有少量稀疏间隔的点被影响。因此，在图象质量方面没有太多牺牲。简单地，少量在理论上甚至应该更亮的孤立的极亮(高强度)点被剪裁。然而，保存这些少量的点是不经济的，其成本是对于整个图象的整体更低的对比度。

确定上下边界曲线

对于第一航路校正，运用一些准则来确定极端的上下边界曲线20、21。这些准则定义在某种意义上试图最佳复制图7a的手画曲线的规则。准则包括：

·形成上下边界曲线20、21的最大和最小曲线由低频分量组成；

·每个曲线20、21的频率足够低(足够平滑)到保留局部细节并且仍然足够高到分离感兴趣的子区域(例如，以保留图3的肋骨内的对比度、肋骨之间的肺部内的对比度以及肋骨和肺之间的对比度)；

·最大和最小曲线20、21相切于局部极值数据点20₁-20_n，21_i-21_n；以及

·计算出的航路22的宽度类似于图7a中的宽度，但要理解，一些点23将在航路22之外。

以最简单的术语，并使用一个数字方法或其他方法，使两个边界曲线20、21拟合于原始输入数据。一个具有低阶或低频率的上曲线20尽可能地与数据14的最大强度拟合。同样具有低阶或低频率的第二曲线21尽可能地与数据14的最小强度拟合。忽略任何非正常值23，在两个边界曲线20、21之间的每个点可以被换算到图象的动态范围。因为采用用于确定上下边界曲线20、21的任何数字方法具有非正常值23都是不可避免的，因此提供上述第二直方图技术来使其对图象增强的影响最小化。

确定上下边界曲线20、21的一种方法是运用移动平均。这仅仅意味着产生一个趋向15，该趋向15由一系列点组成，每点是其周围或相邻数据的平均。一维(1D)分析仅仅是将包含在一预定间隔内的强度数据14求和，然后在样本列11上下对每点重复该过程。2D分析对所分析的每点周围的一个区域内的所有强度数据进行平均。出于这个公开的目的，实现作为一个长方形框的区域的概念，在其中心运用平均或中值。这个框或滤波框的大小为在该点的任一侧一指定数目的行(与在单独的1D分析中的情况一样)以及在该点的任一侧的数个列。

在一个实施例中，通过一个接一个地选择图象的单个列11中的每一个来将移动平均运用到图1的图象是很方便的。可以采用其他方法，例如逐行分析，只影响计算的顺序。

为所选择的列11确定上下边界曲线20、21。2D分析的使用确保不仅在列11中的相邻点之间保持相继性，还保持相邻列11、11、11…之间的相继性。事实上，采用2D分析，例示的逐列分析提供了与逐行分析相同的效果。

每列11由多个相邻点组成，每个与列相交的行中一个点，每点与(动态范围内的)一个强度值相关联。实施第一移动平均，以便使低频或平滑曲线15通过原始强度数据的中值。这个中值曲线将原始数据14分成点的残留上部总体40，那些具有比第一中值曲线15大的强度的点。它还将原始数据14分成点的残留下部总体41，那些具有比第一中值曲线15小的强度的点。

然后，再次采用移动平均技术执行一个迭代过程，以便为点的残留上部和残留下部总体40、41中的每一个计算相继的新中值曲线。随着总体40、41在每次迭代时变得越来越小，每个相继的中值曲线分别越来越向原始数据14的最大和最小强度20_i-20_n和21_i-21_n移动。

应该理解，可以使用其他曲线拟合技术，并且例示的移动平均技术可以进一步被优化。下面公开了一些其他选择。

更详细地，如同运用到仅仅上边界的确定中一样，参看图10a-12b，将一个移动平均技术运用到图4a的原始数据上。参看图10a，计算第一中值低频曲线15并叠加到图4a的原始数据14上。其结果是一部分低频曲线24_i；部分是因为只绘出了大于低频曲线24_i的点的残留总体40。

在这个实施例中，以相同方式实施对上边界曲线20和下边界曲线21的确定，迭代地确定相继中值曲线并减少各自的数据总体40、41。对上边界曲线20的描述同样适用于下边界曲线21。

为了确定上边界曲线20，第一个步骤是替换低于计算出的中值低频分量数据的输入数据。参看图10b，绘出了所得的数据24_i，并且，为了清楚起见，从图中省略显示了较低点(用于确定下边界曲线21)。换句话说，第一个步骤是识别残留上部总体40_i，并生成一个在逐点基础上由输入数据和先前计算出的低频分量数据中的较大者组成的人为图象，在这种情况下是第一中值曲线15。所得的人为数据的总体40_i在局部峰值具有细节，在局部波谷处只具有低频分量。

从人为数据的这个残留总体40_i，计算一个相继的新中值曲线24_i+1。计算由真实数据局部峰值和低频中值曲线24_i+1的波谷中的较大值在逐点基础上组成的这个残留总体的新低频分量，并再次生成由残留总体的最大值和其低频分量组成的第二残留总体40_i+1。

如图11a所示，在一个从i到n次迭代的迭代循环中重复这个过程，随着计算出每个相继的中值曲线24_i+1，每个相继的中值曲线退回成为一个先前中值曲线24_i。

可以使用相同的迭代方法来找到下边界曲线21。唯一的差别是第一和后续残留下部总体40_i到40_n(人为图象)被计算为输入数据和相继曲线26_i到26_n的相继低频分量中的较小值。以这种方式找到上下边界曲线产生一个如图7b所示的适当的计算出的航路22。

在一维中，对于样本列11，如图7b所示，上下边界曲线20、21将大量的原始数据14夹在中间。参看图7c，对于计算实际发生的2D图象，可以为原始数据14的每个相继分析的集合将列数据14的航路22在第二维(相邻列)中扩展，并类似于由每个相邻上边界曲线20，20，20…(20+)和每个相邻下边界曲线21，21，21…(21+)组成的两个平滑表面之间的空间。上下边界表面20+、21+的图象例子分别如图8a和图8b所示。

术语边界曲线20、21和表面20+、21+在这里可交换地使用，并且通常反映分析过程的前后关系，例如图象是在逐行1D曲线基础上还是出于效率的目的在同时定义表面的2D中分析的。

在这个2D环境中，可以以稍微修改的形式重新陈述准则，即：最大和最小边界表面20+、21+由低频分量组成；最大和最小表面相切于局部极值数据点(20_i-20_n或21_i-21_n)；以及，分开表面的航路22接近数据的极值范围，只有少量的非正常值23落在表面之间的空间之外。

参看图11b，最后迭代40_n以一个满足上述每个准则的最大边界曲线20结束。

确定低频曲线的阶数

如上所述，一种确定数据集合的低频分量的常规方法是使用移动平均技术或滤波器，简单陈述为：考虑数据(图象)的一个子集；计算这个子集的平均值(平均强度)；以及，将第二图象数据阵列中的这个值放在子集的中心(中值强度)。

最直接的二维(2D)实现是使用一个在行和列尺寸的每一个中具有奇数点的滤波框，使得中值点在各个数据阵列中容易确定和索引。为了方便起见，公开了一个长方形滤波框；计算框内的点的平均值，并将结果存储在位于框的中心的点。将移动平均滤波框的中心移动一点，然后为下一个输出点计算一个新的平均值。有效算法确保所重复的计算量被最小化，在这个意义上，可能永远不会是框的物理表现。为图象中的所有点执行这个过程。可以使用其他形状，例如圆圈，并且同样可以使用加权的平均方案。

框的尺寸确定了阶数或频率是低还是高。例如，对于原始输入数据的例示列的1007个数据点，一个503点的滤波框假设产生两个中值点或一个零频率的线性曲线，混合图象10的所有数据并弄暗任何对比度变化。如果该框仅仅是一个或三个点的宽度，中值将与原始数据实质相同，从而产生一个根本不保留局部数据的高频曲线。滤波器的尺寸的选择将在下面更详细地讨论。

移动平均滤波器分析具有几个实际的弱点，包括：一个是没有经过边缘的数据，因此需要人为技术来计算边缘附近的低频分量；以及，另一个是由于许多冗余加法而产生的高计算无效性。

参看图12和14c，在本发明中处理图象边缘的方式是在端点或图象边界反照数据。提供一个比图象数据的大小大的临时数据阵列(图象)。将图象数据存储在阵列中，使得在图象数据的开始和结束都有足够空间用于反照数据。反照的点的数量取决于移动平均滤波框的尺寸。最少必须反照滤波框尺寸(col_pad和row_pad的大小)的1/2的点，以将滤波框放在图象边缘，并且仍然计算一个中值，将其存储在图象边缘。

更特别地，参看图12，除了有更大范围的点数之外，初始输入数据与前面的图相同；对于图1的878×1007点图象，进入低于点1的负数和图象边缘1007之外。反照数据44显示在点1之前和点1007之后，由实的垂直线51指示。额外的人为数据正好是真实数据14的逐点反照数据，在这个图上真实数据14的范围是从点1到1007。

这种反照是一种处理图象边缘的实用方式，因为所得的人为反照数据44与真实数据14具有类似的特性，包括：真实数据和人为数据或人为数据和人为数据之间的点对点对比度与在实际输入的原始数据中的对比度一样；以及，输入数据和人为数据的基本统计质量类似。

移动平均滤波器

如图13a和13b的样本Visual Basic(微软公司)代码中所阐述的，通过重复使用一次计算出来并多次运用的值，可以使通常低效率的移动平均变得更有效。当用于计算上下边界曲线20、21时，这个确定在图14a的伪代码中也被称为子例程effc_2dma。以满足所预期的滤波框大小的至少1/2的足够数目的行和列在图象10的边缘的全部周围产生反照数据44。

一旦将图象边缘反照，则可以运用移动平均技术来找到低频边界曲线20、21。如果机械地运用，一个移动平均技术在计算上是低效率的，因为有许多重复的计算(加法)。当移动平均技术要用在一个迭代过程中时这是所特别关心的。

参看图13a的伪代码，在A，将整个图象尺寸初始化为rmin和rmax，包括尺寸的反照数据rows-pad。在B，将滤波框列尺寸初始化为c2bg和c3en。运用一个带有均匀加权(1.0)的长方形移动平均滤波框。

注意，在移动滤波器的最低效应用中，对一行的子集(1D)或一框的点(2D)、滤波框的强度值求和。在这样一个分析中，点的周围滤波框以一点递增，对每点求和，并计算一个中值。操作的数目等于点加法的数目和获得中值的除法。

通过承认大量滤波框平均计算是重复的来提高效率。

参看图13b，在已知的和基本的概念中，在1D中，对于一行图象点IIIIII…，对具有示范尺寸的1行和七列点BBBBBBB的滤波框求和(7个加法)，并除以7以确定中值，如在BBBBBBB的中心的高亮方框中所示。当滤波器递增一列时，从总和中减去最左列中的点的值，并加上最右列中的点的值，然后全部再除以7。结果是只有2个算术操作和一个除法。对于更大的滤波框，这导致节省很大的计算量。

申请人已经改进了这个已知的方法来运用到2D中的工作中，并且可以将其运用到整个图象10，同时避免这些非常大的点阵列中的重复计算。

可以证明，一个行和列方法或一个列和行方法可以产生相同的分析。在图13c中，对相邻行中的多个滤波器子集点求和(7个高，3个在中心之下，3个在中心之上)，然后将所得的行总和在多个滤波器列(7个宽)上求和。在图13d中，对相邻列中的多个滤波器子集点求和(7个宽，3个在中心左边，3个在中心右边)，然后将所得的列总和在多个滤波器行(7个高)上求和，这产生相同的计算出的中值。

图13e显示了如图13b所述的用于一维的有效计算方法的逐行遍历(traverse)的2D应用。

图13d、13e所示的方法只是可视地显示在2D中运用的一个移动平均技术。实际上，如同可以在图13a中C和D处的代码看到的，将分析实际分级成两个迭代步骤和一个规格化步骤。

在第一迭代步骤中，在C，为滤波框尺寸内的多个点计算一维滤波框列总和，并存储在阵列(stot)中。对于整个图象，采用有效的1D方法计算每个相继的滤波框列总和。对于滤波框列总和，用中心点索引阵列(stot)。

在第二迭代步骤中，在D，对于整个图象，执行逐点移动平均分析，有效地同时确定图象中的每一个移动平均曲线。对于每个和所有点，如图13e所示，对于一点的急射(snap shot)，对该点之下、在该点以及该点之上的滤波框列总和(stot)进行求和一一将列总和从stot阵列中拉出来。用滤波框点数或规格化步骤进行的除法还未完成。对于每个相继点，通过减去每个先前的滤波框列总和stot并加上每个相继总和作为计算循环遍历，无论是在逐列还是逐行基础上，运用同一高效概念。

最后，图象中的所有点都用滤波框中的点数(npts)来规格化。

更特别地，参看图13a中的代码，在C，在第一嵌套循环，对于irow和ico2，为第一样本列11的滤波框对相邻列的图象数据求和(滤波框列总和或部分总和)。第一嵌套循环(ico2)为样本列方向(滤波框宽度)的任一侧的列中并在输入数据的第一列的一特定行位置的点计算部分总和。第二嵌套循环(irow)然后逐行并沿图象界限之间的列来递增部分总和例程。

接着，在C，在第二嵌套循环，对于ico1和irow，为每个相继列计算部分总和。在这里，效率如上所述被提高；当从第一列移动到最后一列、并重新计算部分总和时，新总和仅仅是最后一个总和，减去一个结尾点并加上一个在前点。更具体地，在这个第二嵌套循环的最里面的循环中，由于中心点只移动一个位置，下一个递增滤波框的部分求和与舍弃先前的第一单元(ico5)和增加当前最后一个单元(ico4)相同。以这种方式，每个部分求和只实现一次，而不是每次移动框时都实现。对于部分求和的操作总数非常接近于对于图象中每点的一次加法和一次减法。

在D，执行滤波框求和。注意，不需要一个耗费计算的单独步骤来获得滤波框行尺寸的部分和，因为它们在C完成并在一个阵列(stot)中索引。在D的第一嵌套循环通过合并来自C的部分总和与跨在第一行上的行的滤波框中相邻点的总和来计算滤波框。

在D的第二嵌套循环然后为剩下的列和剩下的行相继计算全部滤波框总和。这个第二步骤再次被尽可能高效地完成。由于滤波框的中心点只移动一个位置，新总和与舍弃先前的第一部分总和(ico5)和增加当前最后一个部分总和(ico4)相同。

因此，现在在C和D之间，全部2D求和的操作总数非常接近于对于图象10中每点的两次加法和两次减法。

最后，在E，当我们想知道求和数据的平均值而不仅仅是求和本身时，我们必须用作为框面积的所求和的点数(npts)除以强度总和。

采用上述技术，不需要在物理上保持一个强度的中间滤波框阵列，也不需要计算或操作来产生它们。移动平均所需的全部操作近似等于图象中每点的仅仅4次加法/减法和1次除法——不管要平均多少点。因此，在81×101点(图象尺寸的大约10％)的移动平均长方形框的情况下，这对于这种应用不是不合理的，普通低效移动平均技术将需要图象中每点8,181次加法和1次除法——而与之相对的是4次加法和1次除法现在是可能的。

因此，通过引入一个额外的保存图象数据阵列(stot)来用于部分和，解决了2D移动平均滤波器的计算低效问题。

为图象中的每一个和所有其他列重复这个技术。

迭代技术

如图14a-14c的样本代码和图14d中的定义所阐述的，设置移动平均滤波框尺寸，计算移动平均中值曲线24_i-24_n和26_i-26_n，并向上和向下迭代以形成极值上下边界曲线20、21以及极值曲线或极值表面20+、21+。

滤波框尺寸可以使用预定缺省值。如后面所示，图象10的增强对于框大小的非常大的变化是相当不敏感的。对于一特定类型的图象，缺省值产生可接受的结果。

参看图14a，在F，参考滤波框大小的一种方便方式是将其设置为图象10的尺寸的一个分数。从图象的行数和列数计算框参数rows_pad和cols_pad；简单地，计算为用2除尺寸的百分比，以在框中心的两侧都提供1/2框尺寸。注意，这也是在端点需要被反照的数据点44的数目。

在G，定义阵列界限以保存图象数据的各种临时版本。为了方便，选择变量cmin、cmax、rmin和rmax，以使得参考原始大小的图象10或反照大小的数据阵列中的同一数据点的索引值相同。应该理解，原始图象的存储、读和写以及如何从数据文件或图象阵列提取出列、行和强度数据是系统和程序设计语言所特定的，其方式是本领域技术人员公知的，这里不详细描述。

如图14a所定义的，image_wrk1包含原始数据，image_wrk2包含当前上边界数据，image_wrk3包含当前下边界数据，image_wrk4保存用于高效处理的滤波框总和。这些阵列被初始化到反照的尺寸。

在H，对image_wrk2被初始化，以保存反照原始数据，然后，例程data_mirr和effc_2dma计算反照数据的移动平均，并将其放在image_wrk2中，然后复制到image_wrk3。换句话说，在这个第一阶段，上下极值表面20+、21+的当前近似在这个较早点是相同的，简单地等于移动平均中值曲线15。为了高效的缘故，这个子例程只调用一次。

在I，将一个用于非正常值23的可容忍的故障率frat设置在一个任意值5％(高于或低于)。直方图校正的本质是只要故障率在迭代期间不是很大，就不会有任何问题。还有处理噪声的上表面最小和下表面最大参数surf_mdff。这个参数在后面更详细地描述。

在J，开始迭代。

在J1，迭代的第一步骤，以数据的残留总体40和上表面20+的最后一次迭代50_i中的较大者装载上表面迭代阵列imag_wrk2。以数据的残留总体41和下表面21+的先前迭代51_i中的较小者装载下表面迭代阵列imag_wrk3。就是这个第一步骤强迫两个阵列中的数据成为上下表面20+、21+。

在J2，迭代的第二步骤，将两个阵列中的数据进行反照。

在J3，迭代的第三步骤，找到人为数据集合的残留上部和下部总体40、41的低频分量50、51。从子例程effc_2dma输出极值上下曲线20的两个最新近似，并存储在阵列imag_wrk2和imag_wrk3中。

在J4，第四步骤，需要最小和最大表面20+、21+至少具有一个指定的最小间隔surf_mdff。一旦表面之间的间隔在逐点基础上小于最小值，则重新调整表面，使得它们以最小值间隔开；否则无意义的变化、甚至可能是噪声最后会被换算到整个动态范围。这个间隔是对称实现的，到每个表面一半。在嵌套循环结束处的四个条件“If”语句在先前调整将一个表面推出到动态范围之外这一未必发生、但有可能的事件中重新调整表面。在这个步骤之后，这两个表面将在处处都至少间隔最小值。

这个第四间隔步骤J4也帮助迭代的收敛，因为局限于这个校正的数据不太可能会在表面20+、21+之间的调整后空间之外。

在J5，第五步骤，两个表面20+、21+之间的空间之外的非正常值23的数目建立一个故障率。因此，该过程知道在极值表面的当前尝试不能包含数据的程度。

在J6，迭代的第六和最后一个步骤，确定迭代是否继续或当前的上下边界曲线50_n、51_n是否足够。简单地，如果迭代标志(itrf)被设置为no，则迭代结束，否则，循环重复J1的do loop。如果满足两个条件中的一个：如果正负两侧的故障率都小于可接受的故障率；或者，如果迭代数太高，则itrf标志将只被设置到no，发信号通知收敛。在这种情况下，使用任意数10。需要再次指出的是，由于第二校正，故障率在这里不是那么重要。在不是一直迭代的情况下仍然在曲线拟合的质量和所花费的努力量之间找到一个平衡。

一旦完成了循环，极值上下边界表面50_n、20+和51_n、21+是已知的，可以将其本身作为图象文件输出(如图8a，8b所示)，或为了质量控制目的作为其他计算机文件输出。用于输出这种图象以供显示或存储的过程在这里没有公开，因为这是本领域任何一个技术人员所公知的，它取决于诸如所用的特定图象格式和显示硬件等变量。既然建立了上下边界曲线，可以执行第一航路和第二直方图校正，以增强图1的图象并产生图2的图象。

改进的收敛

如同可以从图11a的多个曲线看出的，迭代解决方案收敛得稍微有些慢。存在一个计算上的动机以获得更快的收敛——每一遍迭代移动得更多。

在一个修改的方法中，如图15a所详细显示的，为下一个迭代中值曲线60分析的数据不仅仅象图11b中一样是残留数据总体40、41，而是包括数据的放大。将残留数据40、41和当前迭代的曲线50、51的极值点之间的差乘以迭代遍数。这有助于更快的收敛，因为每遍对于重复的故障点给出越来越多的校正。有检查来确保没有额外的校正在动态范围的20％之外并且没有值被动态范围之外的值所替换。

参看图15a-17，当运用到图4a的原始数据时，与前面的10次迭代24_i-24_n相比，这个修改的方法只用4次迭代60₁-60₄来收敛(图15b)，并且前面的10次迭代尽管迭代次数更高(10)，仍然产生大于5％的故障率。如图16所示，从该增强的收敛技术产生的航路22不象前面的10次迭代的情况中那样光滑，因为这是较少迭代的自然结果，但这是可接受的。如图17所示，如图6a和6b相比，后续第二直方图校正结果34与前面非常类似。哪一收敛给出更好的结果，在视觉上是不明显的，但显然两者都给出了一个非常好比较的效果，其中动态范围的大部分是局部使用的。

作为替换，非迭代统计技术再次开始地计算以用于找到极值表面、原始输入数据的移动平均值。然后，计算原始数据14和移动平均中值曲线15的差。最后，采用原始数据的相同大小的滤波框作为移动平均来计算图5b的残留数据的标准偏差。采用与图13a的有效2D移动平均滤波器相同的必要逻辑，但修改到适合一个移动标准偏差，可以有效地确定每一点的标准偏差的局部测量值。然后计算极值表面，作为每点的移动平均值+/-(常数*移动标准偏差)。将这个常数设置到一个将典型地产生类似1％故障率的数。根据经验，一个合适的常数看来似乎具有接近2.0的值。注意，在真正的高斯(Gaussian)数据或标准的钟形曲线中，大约98％的数据落在中值的3个标准偏差之内。

第二直方图校正

参看图18a中的编码和图18b中的定义，如上面所讨论的，将数据校正分成2个步骤考虑到了极值表面的估计中的灵活性。表面20+、21+的近似曲线20、21是可接受的。甚至在第一步骤中在极值表面之外的高百分比的数据23也不会产生严重的数据削波。

这是与现有技术的很重要的不同，现有技术运用条件语句以便在图象数据高于一个阈值或低于阈值时不同地处理图象数据；在某种意义上是创建两个单独的数据总体和两个不同的校正。

在本发明中，在第一航路校正时，作出换算范围的初步评价，其识别并确保在直方图及其扩展的换算范围的使用期间甚至非正常值23也不被丢失而是继续被保留并保存在主数据总体14中。最后，一旦尽可能多的数据点已经被包络，则只有最不正常的非正常值最终被剪裁。

在R，航路校正以伪代码的形式显示，其中，小于下曲线imag_wrk3的原始数据imag_wrk1用航路的范围(imag_wrk2-imag_wrk3)换算到动态范围(drim)，将校正的数据存储在imag_wrk4中。

第二直方图校正允许对确切地有多少数据要被削波或剪裁进行控制。并且，由于来自第一航路校正的数据实际上已经是残留的，只有少量要被剪裁的数据也在整个图象中稀疏地分布，不会导致图象中可注意到的失真。应用等式1的航路校正(在R)通过临时存储在比动态范围所需的更高精度的变量中来保留非正常值信息。这个航路校正也考虑了负数。

将下极值曲线21之下的数据作为负数计算，将上极值曲线20之上的值作为输入图象文件/数据的动态范围之上的数计算。然后在通过比较9b和6b可以看出具有与输出数据非常类似的统计质量的数据上执行直方图校正。采用一个比可能发生的更宽范围的事件来作出第一航路校正后的数据的直方图。根据经验，已经发现，由航路校正计算出的值的范围不太可能超过动态范围的三倍，将这些限制任意地用于设置直方图尺寸。

然后通过循环所有第一航路校正后的数据并对每个强度值的出现次数计数来进行直方图计数。

再参看图18a，在S，对直方图正常计数，其中，航路校正数据位于所期望的范围-drlm到drl2内(输出动态范围的三倍)。虽然不太可能有航路校正数据值在这个范围之外，但仍然处理这些可能性。仅仅将范围扩展到3*drlm之外是没有保证的。

因此，在T，在数据点位于所期望的范围之外这样的不太可能的事件中，在直方图的极值处递增直方图计数。预期在极值处的直方图计数在实际上的所有情况中将为零，在甚至差不多例外的情况中非常接近零。Select Case编码在第一逻辑校正Case语句跳出计数例程。进行直方图的累计，直到该值到达在步骤U设置到图象数据的1％的可接受的剪裁(故障)率。

通过使用累计(running total)来找到剪裁点。如在V所示，对于下剪裁点，循环在直方图的低端开始，并向上逐步发展到累计高于可接受的剪裁率。然后用一个值向下重新调整剪裁点的位置，以便使累计刚好低于可接受的剪裁率。将这个值约束到(-drlm)或更大，而不管这样一种情况出现的可能性很遥远。

在W，除了循环在直方图的顶部开始并且向下发展之外，将同一方法用于定位最大或上剪裁点。

一旦已知上下剪裁点的值，则作出从可接受的剪裁界限到图象动态范围界限的直线调整的简单直方图校正。

在一个实施例X1中所示的最简单的情况中，将第一阶校正数据从剪裁界限换算到动态范围。最后，将动态范围之外的任何数据点剪裁到动态范围界限的全部范围。这是本发明的全部校正。

每个采样列的新换算的强度数据中的每一个可以被重新集合，输出图象实际上不能与图2中的图象区别开来。

虽然上面在X1本身是一个非常合适的解决方案，但根据经验发现，先前确定的剪裁点最好被设置到一个不同于刚好设置到动态范围界限的范围，因为人的视觉在较暗区域不太敏感。下面的修改将剪裁点设置到动态范围的20％和95％，这导致在图2的图象中找到有利的结果。这还导致增大数据完整性。图2是只有大约0.2％或大约880000点中的大约1800点最终被剪裁的结果。通常，用于对剪裁数据计数的代码不会被包括在用于提高计算效率中。

在这一点上，全部第一和第二校正图象位于阵列imag_wrk4中，可以作为一个标准图象文件或任何其他类型的计算机文件被输出到介质、例如磁盘上。并且，如果被包含在一个带有GUI界面的软件包中，所处理的数据可以从电子(RAM)存储器显示到监视器或其他显示设备上。

图象校正参数

本发明利用输入图象和两个数字参数。这些参数先前已经在样本代码中提到，它们是：

·surf_mdff，保存表面之间的最小可允许差，作为信噪比S/N管理问题(图14b，J4)；以及

·imag_pcnt，保存要用于移动平均处理(图14b，F)的(两维中的)图象大小的百分比。

如果参数surf_mdff被设置为值8，这意味着两个极值表面20+、21+将被约束为使得它们永远不会小于8个单位的强度差。如果surf_mdff用动态范围的百分比表达，则0-255的3％变为8。

如果参数imag_pcnt被设置为10％，图象的大小是2000×1500点，则移动平均滤波框尺寸将被设置为200×150点，将所有这些点进行平均，以确定一单个低频点的值。注意，每点是30,000个加法，因此需要一个高效的算法。更确切地，滤波框将被设置到151×201点。

虽然这两个参数本身在概念上是足够简明的了，但申请人重想起用户不必应付这种任务是一个所希望的准则。

这两个参数中的每一个都是与另一个无关的。参数surf_mdff管理噪声问题。在这个意义上，噪声被考虑为大约是在每点的数字强度值的随机波动。事实上，surf_mdff是一个用动态范围的百分比表达的关于任何给定点的容限。另一方面，参数imag_pcnt代表要在一个瞬时考虑的点的区域(滤波框)。以笛卡尔坐标表示，imag_pcnt是对x和y轴(域)的控制，surf_mdff是对z轴(范围)的控制。

可以为具有已知和历史噪声分量的图象提供缺省值，并因此完全从用户手中取得。否则，一般指南可以帮助非技术用户从较高级别选择来指定参数，例如图象中的噪声是干净、显著还是严重的，对于这些，surf_mdff的值可以为用户指定在动态范围限制的3％、6％(典型缺省)或12％。

参数数字值由一个3个单元的小的查找表和一个或两个乘法在程序中内部分配。这是特别有用的方法，其中将本发明编码成一个GUI，并交互地完成该过程。用户可以马上估计噪声电平，然后从一个带有上述选择的下拉菜单中进行选择。该参数的坚固性使得三个选择的这个有限的选择应该是非常合适的。注意，这个方案的细节不是限制性的，而是指出实际实现的方式。

关于框大小参数imag_pcnt，其选择和幅度主要是图象内的细节的问题。图象中的细节越多，移动平均滤波框必须越小，以便可以在更小的区域上使用全部动态范围。由于该参数也是坚固的，具有一个低增益，可以对其使用一个类似的方案。在X光图象的环境中，下面的值是根据经验确定的；对于胸部X光——5％，对于手——10％(典型缺省)，或者，对于腿骨，细节级别在图象大小的20％。

用户可以通过仅仅知道研究对象是什么来估计该参数。该参数的坚固性使得三个选择的这个有限的选择应该是非常合适的。

还有一种情况是，如果该过程不在GUI中执行，因为需要这么少的专家经验，并且输入的“批处理”风格也可以被做得很容易，参数估计的这个直接的形式将是特别有用的。

在一个极端，可以将更细级别的选择实现为一个用于技术上内行或严谨的选项，而较粗级别的选项可用于每一个用户。用户输入参数的缺省评价也被找到来产生可接受的结果而不用实际咨询用户。

在任何情况下，以人类可理解的方式实现变量是容易的，使得用户不需要增大他们对该过程的物理和数学知识。

参看图19-21，展示了作为以对surf_mdff的值的不同选择处理同一输入图象(图1)的结果的图象。图19使用1.5％的surf_mdff值，图20a和20b相同，都在3％，图21使用12％。注意，图2、20a和20b是相同的，只是为了比较的目的而复制的。

注意到该参数的双倍才产生了一个显著的效果，这清楚地证明了该参数的坚固本质。分别在1.5和3％的图19、20a的图象实际上是不能区分的。同样，图20b和21实际上是相互不能区分的。申请人知道关于参数可以变化得多大是有限制的。

增大的surf_mdff值更好地控制随机噪声放大，但这也限制更多细微效果的放大。在当前例示图象的情况下，主观分析发现图2(与图20a、20b相同)是较好的校正——随机噪声几乎不引人注意，而细微方面仍然非常清楚。当surf_mdff值较小时，噪声放大变得较不合乎需要，而当surf_mdff值较高时，细微方面的放大减小，与在图21中的情况一样。

然而还可以看到，所有三个图象都代表对图1的一个惊人的改进，因此，surf_mdff值的任一个合理的选择都将产生一个良好的结果。

参看图22-25，展示了作为以对imag_pcnt的不同参数选择处理同一输入图象(图1)的结果的图象。参数surf_mdff值保持为常数，但imag_pcnt值变化。再次使用图2作为参考，并具有图象大小的5％的值。出于比较的目的，将图2复制为图22。图23使用在7％的滤波框，图24在10％，图25在图象尺寸的20％。

通过主观观察到图22(图2)、23或图24中的任一个都是可接受的结果，再次清楚地证明第二个框大小参数的坚固本质。最后，在图25的20％，结果与其他相比比较暗淡。然而，图25比图1仍然有很大提高，必须记住，使用这个X光的20％的值将代表用户可能选择的最坏的情况，如果允许这样的话。可以预期，诊断者不会混淆胸部和腿部图象之间的细节中的差别。

申请人选择的框大小参数设置是5％(图2，22)，建议的缺省设置是10％(图24)，这二者与中间值7％(图23)一样都具有一个非常好的结果。

Claims

1.一种用于使选自一个数字图象的数据总体的相邻点之间的对比度最大化的方法，每个图象点由一个位置和一个强度来定义，图象具有一个强度动态范围，所述方法包括如下步骤：

确定第一低频趋向函数，该函数是对数据总体的最大强度的曲线拟合；

确定一个独立于第一趋向函数的第二低频趋向函数，该函数是对数据总体的最小强度的曲线拟合；

为由第一和第二趋向函数界定的数据总体建立一个最大和最小航路；以及

对于数据总体中的每点，

从航路提取局部最小和最大强度的一个范围；

将一个局部比例因数定义为图象的动态范围和提取出的局部范围之间的比；以及

用局部比例因数换算所述点。

2.如权利要求1所述的方法，其中，每一图象点相对于局部航路的值被换算到一动态范围，该动态范围为

3.如权利要求1所述的方法，其中，对于数据总体存在一个或更多非正常值点，所述非正常值点具有在航路之外的强度，所述方法进一步包括如下步骤：

将一个或多个非正常值点确认为在最大和最小航路之外的点；

用局部比例因数换算一个或多个非正常值点中的每一个点，其结果是换算出的强度在图象的动态范围之外；

以一个大于动态范围的精度存储非正常值点的换算后强度，以防止丢失强度信息；

形成整个数据总体的一个强度直方图，所述直方图建立一个预定义的较低强度和一个预定义的较高强度，用于形成一个大于动态范围并包围基本上所有非正常值的强度的范围；

剪裁具有低于预定义的较低强度和高于预定义的较高强度的强度的点的强度直方图，用于建立一个具有在最小和最大剪裁强度之间的剪裁范围的剪裁数据总体；以及

用剪裁范围与图象的动态范围的比来换算剪裁数据总体的每点。

4.如权利要求3所述的方法，其中，将数据总体调整到适合于人的视觉的更大强度，进一步包括：

将剪裁的数据总体的每点换算到一个小于图象的动态范围的输出范围；以及

将每个已换算的点偏移一个小于动态范围和输出范围之间的差的增量强度值，以便使换算后的点位于一个仍然在动态范围内的较高强度范围内。

5.如权利要求1所述的方法，其中，第一和第二趋向函数中的每一个的确定包括如下步骤：

确定代表整个数据总体的第一函数；

迭代地为具有大于第一函数或先前相继上部函数的较大者的强度的数据总体的点的残留子集确定一个相继上部函数，向上收敛，直到少于预定义数目的点大于相继上部函数，收敛的相继上部函数形成第一趋向函数；以及

迭代地为具有小于第一函数或先前相继下部函数的较小者的强度的数据总体的点的残留子集确定一个相继下部函数，向下收敛，直到少于预定义残留数目的点小于相继下部函数，收敛的相继下部函数形成第二趋向函数。

6.如权利要求5所述的方法，其中，上部函数和下部函数中的一个或两个的收敛由下列步骤提高：

确定点的残留子集和迭代中每个相继函数之间的差；以及

放大所述差，并将其加到点的一个子集，以形成点的一个放大子集；以及

将相继函数的迭代确定运用到所述点的放大子集。

7.如权利要求6所述的方法，其中，通过每次迭代时递增一个计数值并用计数值乘以所述差来放大所述差。

8.如权利要求5所述的方法，其中，用来确定第一趋向函数和第二趋向函数中的每一个的一个或多个函数是移动平均。

9.如权利要求8所述的方法，其中，用于第一中值曲线、相继上部曲线和相继下部曲线中的每一个的函数是移动平均。

10.如权利要求8所述的方法，其中，移动平均具有一个有预定大小的列尺寸和行尺寸的滤波框，并通过下列步骤为图象在两维中优化：

在图象的边缘以与滤波框列和行尺寸互补的许多点扩展图象；

对于滤波框的大小，及关于沿第一列或行尺寸的每个图象点，确定点的一个一维子集的强度的第一总和；

存储所述子集总和，并指向图象中的每点；

对于滤波框的大小，及关于沿一个行或列尺寸的每个图象点，确定子集总和的第二总和；以及

通过除以滤波框中的点数来规格化第二总和。

11.如权利要求10所述的方法，其中，图象具有指定尺寸，滤波框尺寸被设置到图象尺寸的一个百分比，用以调整第一和第二趋向函数中的一个或两个的频率。

12.如权利要求11所述的方法，其中，滤波框尺寸在图象尺寸的5-20％之间。

13.如权利要求12所述的方法，其中，图象是胸部的X光，滤波框尺寸为5％。

14.如权利要求10所述的方法，其中，航路图象的局部范围被约束到不窄于图象动态范围的一个百分比以抑制噪声。

15.如权利要求14所述的方法，其中，航路图象的最小范围被约束到3-12％之间。

16.如权利要求15所述的方法，其中，航路图象的最小范围被约束到6％。

17.如权利要求10所述的方法，其中，在图象边缘通过在边缘反照图象数据来扩展图象。

18.如权利要求10所述的方法，其中，

图象具有指定尺寸，滤波框尺寸被设置到图象尺寸的一个百分比，用于调整第一和第二趋向函数中的一个或两个的频率；

航路图象的局部范围被约束到不窄于图象动态范围的一个百分比以抑制噪声；以及

频率和噪声抑制百分比可由诊断者调整。