CN1488125A - 响应数字图象中的特征对准点阵的方法 - Google Patents

Abstract

一种用于生成考虑数字图象中的特征的点阵的系统和方法。该方法包括一初始化点阵的过程和一相对于图象特征优化该点阵的配置的过程。通过调整点阵点以使一个点阵点空间坐标的复合函数取极值来进行该优化过程。用原子解释这些点阵点。每个点阵点对一个原子势场贡献一个势函数。该图象代表另一个势场。该复合函数是在各点阵点处求值的原子势场和图象势场的加权和。

Description

响应数字图象中的 特征对准点阵的方法

技术领域

本发明涉及数字图象的分析，尤其涉及从这种图象建立计算网格。

背景技术

常常分析数字图象以得到用于进一步计算的网格。例如，分析地震图象以得到用来模拟地下储油层中的石油流动的地质网格。类似地，分析人脑中的医学图象以得到用来模拟动脉中的血液流动的网格。一种略微不同的例子是图象变形，其中从一幅图象导出的网格可用来有效地对相继图象中的差异编码。在所有这些应用中，分析图象以构建网格。

目前的图象和网格建立技术

具体地，这些应用通常包括下述顺序步骤：

(1)处理图象以增强感兴趣的特征。

(2)寻找图象中的对感兴趣的区域定界的曲线或曲面。

(3)用计算网格填充由这些区域定义的空间。

(4)在空间填充网格上模拟某过程。

对于医学成像应用，参见Cebral，J.R.和Lohner，R.的“从医学图象到CFD网格(From Medical Images to CFD Meches)”，Proceedings ofthe 8th International Meshing Roundtable，P.321-331，Sandia，1999。对于地震成像，参见Garrett，S.，Griesbach，S.，Johnson，D.，Jones，D.，Lo，M.Orr，W.和Sword，C.的“地球模型合成(Earth Model Synthesis)”，First Break，v.15，no.1，P13-20，1997.

出于各种原因，目前通常在很少考虑后继步骤的需要的情况下独立地处理该顺序中的每个步骤。例如，通常在步骤(2)中产生比可在步骤(4)所使用的空间填充网格中的表达更加详细地产生界定曲线或曲面。例行地在比实际上石油储藏层模拟中可使用的网格分辨率更高地绘制表示地质界面的地震反射。

分辨率上的不一致常常伴随数据结构的不一致。例如，步骤(2)中通过简单的直线段链表示的二维(2D)曲线可在步骤(3)中变成相对复杂的三角形网格。目前的这种不一致干扰该分析顺序，并且可能产生步骤(1)中要处理的图象和步骤(4)使用的网格之间的难以量化的不一致。

这些不一致是昂贵的。在地震图象分析的例子中，目前该顺序的一次迭代可能需要数月的工作。这种高成本使得难以在试图评估分析结果的不确切性时进行多次迭代。

分辨率和数据结构不一致的一个原因是实现分析的顺序。历史上，上面列出的这四个步骤可能是通过独立开发的由不同的专家使用的不同计算机程序完成的。现在，单个人或一个小组可以以该顺序完成所有四个步骤，并且存在使这些步骤更加一体化的要求。

不一致的另一个原因是方法学上的不足。模拟步骤(4)常常需要偏微分方程(PDE)的数值求解。过去，这些求解只对其中可按简单数组对网格元素编下标的结构化网格是可行的。把步骤(3)中产生的非结构化网格转换成结构化网格产生要处理的图象和模拟中使用的网格之间的进一步的不相容。但是，求解非结构化网格上的PDE的数值方法的最新进展认为可以避免该转换和所造成的不相容。

利用点阵的网格建立

在网格上进行的大多数计算的精度取决于网格元素的规则性。在带有接近等边三角形的高度规则三角形网格上进行的模拟要比在带有细长三角形的不规则网格上进行的模拟更准确。

这就是德洛内三角剖分流行的一个原因。给定一组代表2D网格的结点位置的点，当和所有其它可能的三角剖分比较时，德洛内三角剖分产生几乎接近等边的三角形。但是德洛内三角剖分独自并不保证规则网格。为此，必须细心地选择网格结点的位置。

(对于应用于网格建立问题的德洛内三角部分的综述，参见Bern，M.和Eppstein，D.的“网格生成和最优三角部分(Mesh Generationand Optimal Trangulation)”，Computing in Euclidean Geometry，Du，D.-Z.和Hwang，F.K.编，World Scientific，1995)。

在图象分析领域之外，已经很好地研究了选择网格的最优结点位置问题。该问题的一种解决办法基于观察到由简单晶体点阵定义的一组点的三角剖分产生高规则网格(参见Mello，U.T.和Cavalcanti.P.R.的“把晶体点阵作为样本的用于网格生成的点建立策略(A Point CreationStrategy for Mesh Generation Using Crystal Lattices as Templates)”。Proceedings of the 9th International Meshing Roundtable，P.253-261，Sandia，2000)。但是，均匀的点阵很少地准确地对准要建立网格的对象的边界。从而，一些解决办法从近似均匀的点阵开始，然后利用原子或气泡间物理力的数值模型自动地把网格结点(原子或气泡)移到更加优化的位置。(参见：Shimada，K.“基于物理的网格生成：通过气泡灌注的平面和体积的自动式三角剖分(Physicauy-Based Mesh Generation：Automated Triangulation of Surfaces and Volumes via BubblePacking)”，博士论文，麻省理工学院，1993；Shimada，K.和Gossard，DC.的“气泡网格：通过球填充的非簇几何的自动式三角形网格建立(Bubble Mesh：Automated Triangular Meshing of Non-ManifoldGeometry by Sphere Packing)”，in Proceedings of the ThirdSymposium on Solid Modeling and Applications，P409-419，1995。授予Itoh.T.，Inoue，K.，Yamada，A和Shimada，K.的美国6,124,857号专利“网格建立方法和设备(Meshing Method and Apparatus)”(2000年)把该工作扩展到四边形网格建立。另参见：Bossen，F.J.和Heckbert，P.S.的“用于不匀网格生成的可挠方法(A Pliant Method for AnisotropicMesh Generation)”，Proceedings of 5th International MeshingRoundtable，P63-74，Sandia，1996。)

在所有这些例子中，基于点阵的网格建立从准确定义要建立网格的对象的内、外边界的几何模型起步。从而，需要一种能够解决事先不知道这些边界的问题的网格建立系统和方法。

本说明书中引用的美国专利：

4,908,781，授予Levinthal，C.，和Fine，R.，Computing Device forCalculating Energy and Pairwise Central Forces of ParticleInteractions，1990.

5,596,511，授予Toyoda，S.，Ikeda，H，Hashimoto，E.，和Miyakawa，N.，Computing Method and Apparatus for a Many-Body Problem，1997.

6,124,857，授予Itoh，T.，Inoue，K.，Yamada，A.，和Shimada，K.，Meshing Method and Apparatus，2000.

本说明书中引用的其它参考资料

Bentley，J.L.，和Friedman，J.H.，Data Structures for RangeSearching，Computing Surveys，vol.11，no.4，1979.

Bern，M.，和Eppstein，D.，Mesh Generation and OptimalTriangulation，in Computing in Euclidean Geomaetry，Du，D.-Z.和Hwang.，F.K.eds.，World Scientific，1995.

Bossen，F.J.，和Heckbert，P.S.，A Pliant Method for Anisotropic MeshGeneration，Proceedings of the 5th International MeshingRoundtable，p63-74，Sandia，1996.

Byrd，R.H.，Nocedal，J.，和Schnabel，R.B.，Representations of Quasi-Newton Matrices and Their Use in Limited Memory Methods，TechnicalReport NAM-03，Northwestern University，Department of ElectricalEngineering and Computer Science，1996.

Cebral，J.R.，and Lhner，R.，From Medical Images to CFD Meshes，Proceedings of the 8th International Meshing Roundtable，p.321-331，Sandia，1999.

Garrett，S.，Griesbach，S.，Johnson，D.，Jones，D.，Lo，M.，Orr，W.，和Sword，C.，Earth Model Syntheeis，First Break，v.15，no.1，P.13-20，1997.

Mello，U.T.，和Cavalcanti，P.R.，A Point Creation Strategy for MeshGeneration Using Crystal Lattices as Templates，Proceedings of the 9thInternational Meshing Roundtable，p.253-261，Sandia，2000.

Shimada，K.，Physically-Based Mesh Generation：AutomatedTriangulation of Surfaces and Volumes via Bubble Packing，Ph.D.thesis，Massachusetts Institute of Technology，1993.

Shimada，K.，和Gossard，D.C.，Bubble Mesh：Automated TriangularMeshing of Non-Manifold Geometry by Sphere Packing：Proceedings ofthe Third Symposium on Solid Modeling and Applications，ACMPress，p.409-419，1995.

发明内容

依据本发明的各实施例，一种用于生成考虑数字图象中的特征的点阵的方法包括初始化点阵的过程以及相对于图象特征优化点阵的对准的过程。通过调整点阵点使该点阵点的空间坐标的一个复合函数取极值进行该优化过程。该复合函数可以是点之间的成对(pair-wise)距离的第一函数和从图象导出的第二函数(例如，从点阵点附近图象的采样值算出)的加权组合。

可以从优化的点阵生成网格。该网格可应用于各种各样的应用中的任一个中。例如，该网格可用于模拟和该图象相关的一个过程。作为另一个例子，可以利用该网格压缩图象序列(例如视频流)中的该图象或后续图象。

在一些实施例中，可以利用优化的点阵在生成或不生成中间网格的情况下编码和该图象相关的信息。可设想任何各种类型的相关信息。

至少本发明的一些实施例的一个目的是一种能用下述更加有效的顺序替代前面描述(在对现有技术的说明中)的图象分析顺序的方法：

(1)处理图象以增加感兴趣的特征。

(2)用和图象特征对准的计算网格填充空间。

(3)在该空间填充网格上模拟某过程。

代替寻找图象内的区域边界然后对这些区域建立网格，简单地建立一个对准这些边界的网格。

本发明的一些实施例的其它目的是利用以下方法建立的：

.高规则网格，适用于进一步计算，例如求解偏微分方程；以及

.渐变网格，其中网格结点的密度随图象复杂性或随用户定义的函数变化；以及

.3维网格，在可得到3维图象的情况下。

从对下面的说明和各附图的研究，本发明的各种实施例的这些和其它目的和优点会变得清楚。

附图说明

图1示出本发明方法的主要组成部分之间的数据流动。

图2示出任何二个原子之间归一化距离的成对的力函数(图2a)和势函数(图2b)。

图3示出一个原子的标称距离4(图3a)和8(图3b)对一个采样原子势场的贡献。

图4示出三角形化后的原子的六边形二维点阵。

图5是原子的面心立方系三维点阵。

图6是地震图象，已经为了增强地下地质的断层、不连续做了处理。

图7是标称距离函数，其代表点阵中在原子间所需的可变间隔。

图8是在点阵中的原子的初始伪规则分布；每个原子用一个其直径和该原子位置处算出的标称距离函数成比例的圆来标记。

图9是和地震图象的伪规则初始点阵对应的网格。

图10是和地震图象的伪随机初始点阵对应的网格。

图11是和地震图象的优化点阵对应的网格。

图12是该网格中的和该地震图象中的特征最对准的各边缘。

图13是和优化点阵对应的Voronoi网格，其中用地震图象中的特征推斥原子。

图14是人的头部的磁共振图象(MRI)。

图15为增强边缘特征处理后的人的头部图象。

图16是和人头部图象中的特征对齐的网格。

图17是该网格中各条最对准人头部图象中的特征的边缘。

图18是人躯干的磁共振图象(MRI)。

图19为增强边缘特征处理后的人躯干图象。

图20是和人躯干图象中的特征对齐的网格。

图21是该网格中最对准人躯干图象中的特征的各边缘。

图22是用来实现本发明的方法的计算设备的组成部分。

图23是本发明方法的一实施例。

尽管本发明易于有各种修改和替代形式，但仍在附图中以例子的方式示出各特定实施例并在文中予以详细说明。然而应该理解，各附图以及对它们的详细说明不意味着把本发明限制在所公开的具体形式上，恰恰相反，本发明覆盖由附后权利要求书定义的本发明的精神和范围内的所有修改、等同和替代。

具体实施方式

本文说明一种在考虑数字图象中的特征情况下生成点阵的方法。该方法包括初始化点阵并然后优化点阵。在一组实施例中，优化点阵的点趋于对准(吻合)图象中的各特征。在第二组实施例中，优化点阵的点趋于沿边地对准，而不是对准图象中的特征。在另一些实施例中，优化点阵的点可趋向对准某些特征但并不对准同一图象中的其它特征。

以下，把点阵的点称为原子。利用晶体的原子结构说明点阵的初始化过程。类似地，利用使原子点阵的势能为最小说明点阵优化。这些物理模型使得更容易理解本文中的说明。但是本发明的语境下的原子仅仅是一个点。

图1更详细地说明本发明的方法的主要组成部分之间的数据源。数据输入部分100产生图象，该图象已经另外得到处理以增强感兴趣的特征。这里，该图象或者是作为该处理的结果直接得到的或者是从计算机存储器装入的。点阵初始化器200建立一个覆盖该图象的初始原子点阵。点阵优化器300移动原子以使点阵对准该图象内的特征。数据输出部分400使用该优化的原子点阵，或者把它传到某个其他处理或者在计算机存储器存储它。

点阵优化器300包括一个类函数最小化器310和一个势能计算机320。该类最小化器310迭代地搜索多变量函数的最小值。在该搜索期间，该最小化器需要反复地计算该函数以及它相对于每个变量的偏导数。这里，这些变量是点阵中原子的空间坐标。给定原子坐标以及图象，势能计算机320计算点阵势能和势能及其相对于每个原子坐标的偏导数。

计算势能

二维(2D)空间中的原子具有x和y坐标，而三维(3D)空间中的原子具有x、y和z坐标。向量x代表2D(或3D)空间中点的x和y(或x、y和z)空间坐标。给定二个位置为x_i和x_j的原子，|x_i-x_j|代表它们之间的距离。在本优选实施例中，用来计算该原子间距离的向量范数|·|是欧几里得范数。但是，可以替代地采用各种其它任何范数。

成对势函数。

出于计算效率的考虑，利用一个简单的成对力函数作为原子间相互作用的模型，从而一个原子的邻居施加到该原子上的总力简单地为它们中的每一个施加的力的和。即使在这样的简化下，对该成对力函数仍存在许多合理选择。

为了避免二个或更多的原子具有相同或几乎相同的坐标，如果它们互相过分接近则它们之间的力可以是推斥的(正的)。类似地，为了防止出现大的无原子的空间，如果二个原子彼此过分远离则它们之间的力可以是吸引的(负的)。为了便于数值计算，力可以是有界的。为了防止点阵中的每个原子向每个其它原子施加力，超出截止距离后力可能为零。另外，力函数作为原子间距离的函数可以是连续的。Shimada(1993)提出的力函数具有这些特性。从而，在一组实施例中，我们采用Shimada的力函数，如后面说明那样。然而，可以采用各种和这些特性相容的任何力函数。

令d为二个原子之间的标称距离，在这种距离下力从推斥变为吸引。从而，位于x_i和x_j处的二个原子之间的力可由下列的三次多项式给出：

其中u是二个原子之间按如下定义的归一化距离

$u\≡\frac{|x_i-x_j|}{d}$

该多项式函数中的各系数保证该力是有界和连续的，它在u-1和 $u\&GreaterEqual;\frac{3}{2}$ 时等于零，并且在O≤u＜1之间为正，和在 $1<u<\frac{3}{2}$ 之间为负，图2a示出该力函数。

原子上的力通常是一个向量。这里，该向量的方向由f(u)的符号以及二个原子的位置x_i和x_j暗指。

利用标量的势要比利用向量力的多个分量要方便得多。从而，按照周知的惯例，我们把该力定义为一个标量势的负梯度：

积分常数 选择成使当 $u=\frac{3}{2}$ 时φ(u)是连续的。图2b示出该势函数。如所期望那样，势函数φ(u)在归一化距离u＝1处取最小值，在此处力函数f(u)为零。

势能和势场

给定归一化距离u的势函数φ(u)，我们把原子势能定义为成对势的下述和：

$A=A(x_1,x_2,...,x_n)\&equiv;\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&phi;}\left[\frac{|x_i-x_j|}{d\left(x_j\right)}\right],----\left(1\right)$

其中x₁，x₂…，x_n是点阵中n个原子的坐标，而d(x)是位置x的标称原子间距离函数。标称距离函数d(x)不必是常数；但是，为了确保平滑的渐变点阵，我们要求它是平滑的。具体地说，我们要求|d|＜＜1，从而对于小于该势函数φ(u)的截止距离 的|x_i-x_j|/d，d(x_i)≈d(x_j)。此外，因子1/2补偿在总势能A的定义中φ[|x_i-x_j|/d(x_j)]≈φ[|x_j-x_i|/d(x_i)]出现二次。

还利用原子势场定义原子势能A

$a\left(x\right)\&equiv;\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\mathrm{\&phi;}\frac{|x-x_j|}{d\left(x_j\right)}],----\left(2\right)$

从而

$A=A(x_1,x_2,...x_n)\&equiv;\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}a\left(x_i\right).$

换言之，把原子势能定义为在原子坐标处计算原子势场得出的值之和的一半。

类似地，定义图象势能

$B=B(x_1,x_2,...,x_n)\&equiv;\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}b\left(x_i\right),----\left(3\right)$

其中b(x)是基于输入图象的图象势场。

在许多语境下，图象势场简单地是一幅图象(或图象的平滑版本)，典型地由存储在计算机存储器内的2D(或3D)数值数组表示。这里，我们用术语“势场”强调(并且后面利用)原子和图象势场之间的相似性。

在一组实施例中，假定图象已被处理过，从而图象势场得到一个在感兴趣的特征内的最小值(例如b(x)≈-1)和一个在不感兴趣的区域内的最大值(例如b(x)≈0)。在这种情况下，使图象势能B为最小等同于把原子移动到和该图象中感兴趣的特征对应的最小值内。

在第二组实施例中，假定图象已被处理过，从而图象势场得到一个在感兴趣的特征内的最大值(例如b(x)≈0)和一个在不感兴趣的区域内的最小值(例如b(x)≈0)。在这种情况下，使图象势能B为最小等同于从和图象中感兴趣的区域对应的最大值中移出原子。

在第三组实施例中，图象势场可沿图象特征的第一子集得到一个最大值(例如b(x)≈1)，沿图象特征的第二子集得到一个最小值(例如b(x)≈-1)，并且在这些图象特征之外的中间值(例如b(x)≈0)。

更一般地，我们移动原子以使原子和图象势能的下述加权和为最小：

      P＝P(x₁，x₂，…x_n)≡(1-β)A+βB，    (4)把它称为总势能。标量因子β确定A和B对总势能P的相对贡献。当β＝0时，原子趋向一个不和该图象对齐的完美规则点阵。当β＝1时，原子只对图象采样值敏感；即，原子只根据它们向图象中的各最小值吸引和/或从图象中的各最大值推斥移动。这样，原子趋于在图象中会聚到各最小值并腾空各最大值，由此产生高度不规则点阵。典型地，选择

以得到一个接近规则并且考虑图象各特征的点阵。

根据势场a(x)和b(x)，总势能为

$P=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{1}{2}(1-\mathrm{\&beta;})a\left(x_i\right)+\mathrm{\&beta;b}\left(x_i\right).$

根据定义为

      p(x)≡(1-β)a(x)+βb(x)，          (5)的总势场，总势能为

$P=\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}p\left(x_i\right)+\mathrm{\&beta;b}\left(x_i\right).----\left(6\right)$

类似于标称距离函数d(x)，式(5)和(6)中的标量因子β可以是位置x的平滑变化函数。(至于d(x)，平滑性意指β(x)的导数是可忽略不计的。)这种生成使得能在空间上改变点阵规则性和对图象特征的灵敏性(例如吸引或排斥)之间的平衡。点阵对图象特征的灵敏性在图象的一个部分中可能要比某个其它部分中更为重要。为简单起见，令β代表一个常数标量因子。

总势能P是一个带有多个局部最小值的原子坐标x₁、x₂、…，x_n的非二次函数。(例如，请注意任何二个原子的坐标互换不改变P。)从而对最小值(通常接近原始点阵坐标)的搜索必须是迭代的。在有效迭代搜索中，必须相对于原子坐标计算P的偏导数。例如，考虑P相对于第i个原子的x坐标的变化：

$\frac{\&PartialD;P}{\&PartialD;x_i}=(1-\mathrm{\&beta;})\frac{\&PartialD;A}{\&PartialD;x_i}+\mathrm{\&beta;}\frac{\&PartialD;B}{\&PartialD;x_i}.----\left(7\right)$

在计算原子势能的偏导数A/xi中，回想式(1)的双重求和中项φ[|x_i-x_j|/d(x_j)]≈φ[|x_i-x_j|/d(x_i)]出现二次。从而，

$\frac{\&PartialD;A}{\&PartialD;x_i}=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&phi;}}^{\&prime;}\left[\frac{|x_i-x_j|}{d\left(x_j\right)}\right]\frac{1}{d\left(x_j\right)}\frac{x_i-x_j}{|x_i-x_j|}----\left(8\right)$

$=\frac{\&PartialD;a}{\&PartialD;x}\left(x_i\right)$

$\frac{\&PartialD;B}{\&PartialD;x_i}=\frac{\&PartialD;b}{\&PartialD;x}\left(x_i\right)----\left(9\right)$

$\frac{\&PartialD;P}{\&PartialD;x_i}=\frac{\&PartialD;p}{\&PartialD;x}\left(x_i\right)----\left(10\right)$

可以容易地得到相对于每个原子的y(以及三维中的z)坐标的偏导数的类似结果。

计算

为了计算总势能P，需要计算它的组成部分A和B。为了根据式(3)计算图象势能B，必须计算每个原子位置x＝x_i处的图象势场b(x)。图象典型地是均匀采样的，对b(x_i)的最简单和最有效的逼近是离点x_i最近的图象采样处的图象势场b(x)的值。更准确的逼近(内插)是可能的，但在本文中示出的所有例子中均采用该简单和快速的最近邻近内插法。式(3)暗指计算B的成本(计算复杂性)是O(n)，其中n是原子的数量。

相反，式(1)中的双重求和暗指计算A的最直接方法的成本是O(n²)。在实际应用中，n是足够大的，从而计算A的O(n²)成本会远大于计算B的O(n)成本。为了降低计算A的成本，利用对势函数φ(u)的设计，对于大于截止距离

的归一化距离u其为零。只有那些离位置x＝xi处的原子最近的原子才对该位置上的原子势场a(x_i)有贡献。

该观察引出确定哪些邻原子在各个位于x＝xi处的原子的 $\frac{3}{2}d\left(x_i\right)$ 距离之内的问题。该问题的求解是非直接的，因为在优化点阵的期间原子反复地移动。

例如，如果建立相邻原子列表，每个原子一个表，一旦原子移动必须更新这些表(或者至少检查它们是否需要更新)。对于密度接近常量的点阵，利用简单的数据结构建立并更新这些表的成本是O(mn)，其中m是截止距离内的邻原子的平均数量。对于密度变化的点阵，需要更复杂的数据结构并且成本变成O(mnlogn)。(参见，例如，Bentley，J.L.，和Friedman，J.H.的“用于范围搜索的数据结构(Data Stricture forRange Searching)”，Computing Surveys，vol.11，no.4，1979.另参见授予Levinthal，C.和Fine，R.的美国4,908,781号专利“用于计算粒子交互作用的能量和成对中心力的计算设备(Computing Device for CalculatingEnergy and Pairwise Central FOrces of Particle Interactions)”，1990，以及授予Toyoda，S.，Ikeda，H，Hashimoto，E.和Miyakawa的美国5,596,511号专利“用于多体问题的计算方法和设备”(ComputingMethod and Apparatus for a Many-Body Problem)，1997。)

利用原子势场a(x)对原子势能A的表达式提出一种较简单的解决办法。把式(2)解释为一种计算象图象势场b(x)那样采样的原子势场a(x)的方法。具体地，用一个维数和用来表示图象b(x)的数组的维数相同的二维或三维数组表示a(x)。首先对所有采样的x把a(x)初始化为零。接着，对于每个位于位置x＝x_j的原子，累计一个采样势函数φ[|x_i-x_j|/d(x_j)]。在空间上把该累计限制在以位置x_j为中心以 (x_j)为半径的圆(或球)内的采样上，在此采样的势函数的贡献是非零的。

图3a和3b示出二个这样的标称距离分别为d＝4和d＝8的势函数。黑和白之间的灰度分别对应-0.05和0.05之间的采样函数值。原子势场a(x)简单地是许多这样的函数的累计。出于计算效率考虑，这些采样势函数可以是预先计算的并且可对不同的标称距离d制表。接着，给定任何位置x的d(x)，可以高效地通过查表(或者查表并插值)确定势函数的适当值。

优选实施例

上述分析建议了二种很不相同的计算总势能和它的偏导数的算法。本发明的优选实施例利用式(2)和(5)计算总势场p(x)，并接着利

用式(6)计算总势能和利用式(10)计算它的偏导数。下面的伪代码列表详细说明该算法。

算法1：计算P， $\frac{\&PartialD;p}{\&PartialD;x_i},\frac{\&PartialD;p}{\&PartialD;y_i},$ 和

101：初始化总势场p(x)＝βb(x)

102：对所有的原子位置x_j＝x₁，x₂，…x_n{

103：  对所有 $|x-x_j|<\frac{3}{d}d\left(x_j\right)$ 的x{

104：    累计p(x)＝p(x)+(1-β)φ[|x-x_j|/d(x_j)]

105：     }

106：  }

107：初始化总势能P＝0

108：  对所有原子位置x＝x₁，x₂，…x_n{

109：    累计 $P=P+\frac{1}{2}[p\left(x_i\right)+\mathrm{\&beta;b}\left(x_i\right)]$

110：    计算 $\frac{\&PartialD;p}{\&PartialD;x_i}=\frac{1}{2}[p(x_i+1,y_i,z_i)-p(x_i-1,y_i,z_i)]$

111：    计算 $\frac{\&PartialD;p}{\&PartialD;y_i}=\frac{1}{2}[p(x_i,y_i+1,z_i)-p(x_i,y_i-1,z_i)]$

112：    计算 $\frac{\&PartialD;p}{\&PartialD;z_i}=\frac{1}{2}[p(x_i,y_i,z_i+1)-p(x_i,y_i,z_i-1)]$

113：  }

行101到106计算类似于图象势场b(x)采样的总势场p(x)。给定该场，行107至113计算总势能P和它的偏导数 $\frac{\&PartialD;p}{\&PartialD;x_i},\frac{\&PartialD;p}{\&PartialD;y_i},$ 和

在行109，可以如前面所述通过选择最近图象采样位置处的对应场值或者通过任何所需的插值方法，近似原子位置x_i处的总势场值和图象势场值。在行110-112利用简单的中心有限差分近似从总势场计算各偏导数；当然，可采用替代的(例如高阶)对各导数的数值近似。为明确性起见，该列表取三维坐标空间。对于二维空间，简单地不计z坐标和偏导数

即可。

假定原子位置和标称距离函数d(x)一致，算法1的计算成本为O(N)，其中N为图象中的样本数量。回想象对图象那样对总势场取样，并且每个原子对总势场中的那些最靠近该原子的样本贡献空间受限的势函数(类似图3)。从而，累计来自所有原子的作用和该场中的样本N的数量成比例。

该算法的成本和常规图象处理的成本不相上下，并且它所要求的数据结构不比表示图象的简单数组更复杂。另外，非常量标称距离函数d(x)和常量标称距离d的成本是一样的。

替代实施例

本发明的一种替代实施例不计算总势场。相反，它利用式(1)、(3)和(4)计算总势能，并且利用式(7)、(8)和(9)计算它的各偏导数。下面的伪代码列表详细地说明该算法。

算法2：计算P， $\frac{\&PartialD;p}{\&PartialD;x_i},\frac{\&PartialD;p}{\&PartialD;y_i},$ 和

201：初始化总势能P＝0

202：对所有的原子位置x_j＝x₁，x₂，…x_n{

203：  累计P＝P+βb(x_i)

204：  初始化 $\frac{\&PartialD;p}{\&PartialD;x_i}=\frac{1}{2}\mathrm{\&beta;}[b(x_i+1,y_i,z_i)-b(x_i-1,y_i,z_i)]$

205：  初始化 $\frac{\&PartialD;p}{\&PartialD;y_i}=\frac{1}{2}\mathrm{\&beta;}[b(x_i,y_i+1,z_i)-b(x_i,y_i-1,z_i)]$

206：  初始化 $\frac{\&PartialD;p}{\&PartialD;z_i}=\frac{1}{2}\mathrm{\&beta;}[b(x_i,y_i,z_i+1)-p(x_i,y_i,z_i-1)]$

207：对所有 $|x_i-x_j|<\frac{3}{2}d\left(x_i\right)$ 的原子位置x_j{

208：  累计 $P=P+\frac{1}{2}(1-\mathrm{\&beta;})\mathrm{\&phi;}\left[\right|x_i-x_j|/d\left(x_i\right)]$

209：  计算Δ＝(1-β)φ’[|x_i-x_j|/d(x_i)]/[d(x_i)|x_i-x_j|]

210：  累计 $\frac{\&PartialD;p}{{\&PartialD;x}_i}=\frac{\&PartialD;p}{\&PartialD;x_i}+\mathrm{\&Delta;}(x_i-x_j)$

211：  累计 $\frac{\&PartialD;p}{\&PartialD;y_i}=\frac{\&PartialD;p}{\&PartialD;y_i}+\mathrm{\&Delta;}(y_i-y_j)$

212：  累计 $\frac{\&PartialD;p}{\&PartialD;z_i}=\frac{\&PartialD;p}{\&PartialD;z_i}+\mathrm{\&Delta;}(z_i-z_j)$

213：  }

214：}

对于每个原子，行203和行208分别将每个邻原子的图象势能B和原子势能A累计到总势能P中。类似地，行204至206和210至212分别把图象和原子势能的各偏导数的贡献累计到总势能的偏导数中。

行207意味着使用能够快速确定最靠近位于x＝x_i的原子的那些邻原子的辅助数据结构。一旦原子移动，必须重构或者以某种方式更新该数据结构。在任何迭代地搜索使总势能为最小的原子坐标中，原子反复移动。从而，保持该数据结构的成本是相当大的。对于最有效率的数据结构，用于非常量标称距离函数d(x)的成本要高于用于常量标称距离d的成本。

点阵初始化

如前面所述，总势能P是一个原子坐标的带有许多局部最小值的非二次函数。在优化点阵期间，在搜索最小值中迭代地移动原子。实践中，既不寻求也不找出全局最小值。相反，找出一个接近原始点阵的优化原子点阵。从而，希望初始点阵

.使原子势能(局部)最小，

.高度规则，以及

.和标称距离函数d(x)相容。

常量标称距离

对于常量标称距离，我们可以容易地建立一个具有这些特性的初始点阵。图4示出一个这样的用于二维空间的点阵。在该理想点阵中，可以把这些原子连接为形成等边三角形。任何一个原子和它的六个最近邻原子之间的距离简单地是该常量标称距离d。在此距离下，任何一个原子对另一个原子施加的力准确地为零，并且原子势能局部为最小。

图5示出一个用于三维空间的规则点阵。这是一个面心立方系(FCC)点阵，其中原子排列在看起来象图4的二维点阵的各水平层中，但是每个层略微偏离以填充到上下各层的洞中。为了使图5清晰，用灰度级不同的球描绘不同层中的原子。任何原子和它的十二个最近邻原子的距离等于该常量标称距离d。

和能够完全填充二维空间的等边三角形(如图4中所示)不同，等边四面体不能填充三维空间。尽管如此，FCC点阵中的原子可以三角形化以得到高规则四面体。

可变标称距离

对于非常量标称距离函数d(x)，初始排列原子更困难。

第一个复杂性是，除非另外规定必须计算函数d(x)。作为如何计算该函数的一个例子，研究图6中显示的图象(即图象势场b(x)。该图象是为增强地下地质中的断层、不连续性做过处理的三维地质图象(未示出)的水平切片。该图象中的暗线性特征代表和该水平切片相交的断层踪迹。这些断层是近似垂直的，几乎和该水平切片正交。

图7示出平滑该地震图象以得到对标称距离函数d(x)的估计的结果。图中最暗的区域对应于最小距离d＝6个采样，而最淡区域区域对应于最大距离d＝18个采样。根据在该地震图象中观察到的细节等级，明确地规定这些最小和最大距离。在该图象的中左部分中距离较小，而在右下部中距离较大。平均距离约为9.7个采样。

在一些应用中，可以明确地规定或者利用一个用于图象绘制的计算机程序交互地建立标称距离函数d(x)。对于本发明，计算函数d(x)的实际方式不重要。如前面所述，我们只要求该函数是平滑的，即，|d|＜＜1。

第二个复杂性是在点阵中和标称距离函数d(x)一致地排列原子。这里，我们说明二种用于该排列的算法。

优选实施例

本发明的该优选实施例使用一种生成伪规则点阵的算法。下面的伪代码列表详细说明该算法。

算法3：初始化伪规则点阵

301：初始化布尔标志图象状数组ω(x)＝假

302：建立一个空的原子列表list

303：建立一个空的原子地点队列queue

304：把该图象中心的位置x_i添加到queue中

305：当queue不空{

306：  从queue得到并取出第一地点x_i

307：  若x_i位于该图象的坐标边界内{

308：    设置Sphere＝以x_i为中心γd(x_i)为直径的球区

309：    若，对于Sphere内的所有采样，ω(x)＝假{

310：       对Sphere内的所有采样，设置ω(x)＝真

311：       把坐标为x_i的原子加到list上

312：       在queue的尾部添加邻原子的理想位置

313：       }

314：   }

315：}

行301中初始化的数组ω(x)是一个临时工作数组，其维数和图象的维数相同。它的唯一目的是当由该算法生成原子时在该点阵中标记原子的位置。(行309中的检查确保这样标记的位置不会被再次标记。)在本发明的该优选实施例中，该数组可以和算法1中用来存储总势场p(x)的数组为同一个数组，从而不需要附加的存储器。图8示出数组ω(x)的一个例子，它是从图7中所示的标称距离函数计算出的。

图8中的圆形区域以某个原子位置为中心，并且其直径和该原子位置处的标称距离函数d(x)的值成比例。比例常数是该算法的行308中的因子γ。通过把该因子选择成小于1，允许该初始点阵中的一些原子要比该标称距离函数d(x)所蕴含的更加彼此接近并知道其它原子会彼此更加远离。以试验方法确定因子γ＝0.8产生和平滑距离函数d(x)一致的伪规则点阵。

行312中的理想地点是图4和5中示出的规则点阵的邻原子的位置。(这些邻原子的距离不按γ标定。)从而，对于常量d，算法3产生一个和它们中的一个类似的规则点阵。对于非常量的d(x)，它产生一个伪规则点阵。

在任何情况下，把理想地点放到queue中的处理造成该点阵从放到queue中的第一地点向外生长。从而，该第一地点充当一个点阵从其生长的种子。算法3的行304把该第一地点选为该图象的中心。可以采用替代的种子位置。例如，如果需要一个用于模拟液体流动的网格，则可把一个或多个流体源的位置作为点阵生长的种子。作为另一个例子，可以把图象中的特征的质心选为种子位置。

图9中示出的网格是利用算法3以及图7中所示的标称距离函数建立的伪规则初始点阵的德洛内三角剖分。尽管该初始点阵的网格中的大部分三角形是高度规则的，但所示出的整个德洛内三角剖分在点阵点的凸壳附近产生一些细长三角形。在相继的计算中可以简单地不计这些网格边界附近的不规则三角形。

替代实施例

本发明的一替代实施例采用一种完全不同的生成伪随机初始点阵的算法。下面的伪代码列表详细说明该算法。

算法4：初始化伪随机点阵

401：建立一个空的原子列表list

402：对通过该图采样的所有x{

403：  设置d＝d(x)

404：  若二维，设置 $\mathrm{\&rho;}=2/\sqrt{3}{d}^2$

405：  若三维，设置 $\mathrm{\&rho;}=2/\sqrt{2}{d}^3$

406：  产生一个在[0:1]中均匀分布的伪随机数r

407：  若r＜ρ，把x处的原子加到list中

408：}

该算法行404(用于二维图象)或行405(用于三维图象)计算和标称距离函数d(x)的值对应的标称点阵密度ρ。在二维空间中，点阵密度和原子之间的距离的平方成反比；在三维空间中，它和距离的立方成反比。比例常数

(行404)和

(行405)分别与图4和5中示出的理想点阵对应。

点阵密度是在任何对该图象采样的位置上存在一个原子的概率。行406和407利用一个伪随机数发生器按该概率对该点阵添加一个原子。该算法使用平常的用于在间隔[0:1]上均匀分布地生成一个伪随机数的计算机软件。

该伪随机算法要比伪规则算法3更简单和更快。它还产生完全各向同性的点阵，不具有任何由伪规则点阵展示的对称性上的优选方向或平面。取决于应用，这可以是或不是优点。

遗憾的是，伪随机初始点阵是很不规则的。这种点阵里的原子不展示几何模式，这使得术语“伪随机点阵”换句话是矛盾修饰。

图10示出和利用算法4为该地震图象产生的伪随机初始点阵相对应的网格。尽管该网格统计上和图7中示出的标称距离函数一致，但它是很不规则的，带有许多离等边差得很远的三角形。点阵优化会使该伪随机初始点阵列规则，但要比优化图9中示出伪规则初始点阵需要更多的工作(更多迭代)。

点阵优化

点阵优化器移动初始点阵里的原子以得到优化的点阵。在一组实施例中，优化后的点阵是规则的并且要比初始点阵更加对准图象特征。换言之，在优化的点阵里，远离图象特征的原子趋于朝向考虑标称距离函数的伪规则结构，而靠近图象特征的原子趋于和这些特征重合。在另一组实施例中，优化后的点阵是规则的并且要比初始阵列更纯地沿着图象特征。该点阵优化器使用平常的用来使带有许多变量的任意函数为最小的计算机软件。该点阵优化器应用该软件以使作为原子坐标的函数的点阵总势能为最小。

可对类函数求最小器(minimizer)作出各种各样的选择。在本优选实施例中，选择存储器有限的Broyden-Fletcher-Goldfarb-Shanno(L-BFGS)求最小器。(例如参见，Byra，R.H.，Nocedal，J.和Schnabel，R.B.的“伪牛顿矩阵表达式和它们在存储器有限方法中的使用(Representations of Quasi-Newton Matrices and Their Use in LimitedMemory Methods)”，技术报告NAM-03，西北大学，电子工程和计算机科学系，1996。尤其请注意第17页上的简单双循环递归。)和其它求最小器一样，L-BFGS法在搜索最小值中迭代地评估某函数和它的各个偏导数。

和另一种周知的共轭梯度法相比，L-BFGS求最小器需要较多的计算机存储器但函数求值较少。但是，当和表示图象所需的存储器相比时，要增加的存储量是不多的。另外，每次函数求值(计算点阵总势能)的成本明显地比该求最小器进行的其它计算更为节约。从而，L-BFGS求最小器很适用于本发明。

优选实施例

本发明的该优选实施例使用下述点阵优化算法。

算法5：优化点阵

501：得到初始点阵原子坐标x₁、x₂、…，x_n

502：建立势能计算机

503：建立求最小器

504：计算初始点阵总势能P

505：循环{

506：  设置P_o＝P

507：  随机扰动x₁、x₂、…，x_n

508：  循环{

509：    设置Pi＝P

510：    通过调整x₁、x₂、…，x_n使求最小器减小P

511：   }当P_i-P＞ε|P_i|

512：   }当P_o-P＞|P_o|

该算法在行501从一个初始原子点阵，例如通过算法3产生的伪规则点阵或者通过算法4生成的伪随机点阵开始。它接着(行502)建立一个负责计算总势能P以及它的各偏导数的势能计算。它接着(行503)建立一个利用该势能计算机使P最小的求最小器。(建立势能计算机和求最小器包括分配存储器以及对一些变量和表格的初始化。)它接着(行504)计算该初始点阵的总势能P。

该算法的剩余部分由二个嵌套循环构成。开始于行508的内循环使该求最小器调整原子坐标x₁、x₂、…，x_n以减小总势能P。在由行511中的小阈值ε判定P的减小变成不显著之前，继续该循环。典型的阀值为ε＝0.001。

回想该总势能是一个具有许多局部最小值的函数。该内循环从目前的原子坐标开始，并且趋于朝向和这些坐标最近的最小值。观察到该局部最小值可能具有比附近的另一个最小值大的总势能。

从行505开始的外循环使该算法从一个局部最小值移动到另一个最小值，直至总势能P上的减小是不显著的。行507中原子各坐标的随机扰动是小的，典型地对每个原子位置x_i不大于标称距离d(x_i)的10％。利用一个平常的伪随机数发生器计算这些扰动。内部最小化循环的相继迭代典型地产生总势能P的显著减小。

算法5中的内、外循环采用相同的收敛检查。当总势能P的减小变成不显著时二个循环都停止。替代的收敛准则在数值优化中是寻常的，而在此的选择不是重要的。例如，可以当原子各坐标的最大改变小于某阈值时结束这些循环。

图11示出优化用于该地震图象的初始伪规则点阵的结果，其中把该图象处理成沿着图象特征具有接近-1的值和离开图象特征为0的值。回想式(4)，其中指出总势能P是原子势能和图象势能的加权和。在该例子中，采用图象加权β＝0.3。所产生的优化点阵既是高度规则的(同样，忽略点阵点包壳附近的三角形)又是良好对准图象特征的。

对于具有小的或窄的特征的图象，利用该图象的略加平滑的版本进行算法5的前几次迭代是有用的。最初的这几次迭代能使原子向那些由于它们离初始点阵位置过远在别的情况下可能被错过的特征移动。原子最初被吸向被平滑的特征，然后在相继的迭代中被吸向在原始未平滑图象中的高分辨率特征。

曲线和曲面提取

在空间填充网格中，隐含地出现曲线(二维中)或曲面(三维中)。例如，在图11中示出的二维网格中，任何相连的三角形直边序列代表一条曲线。在三维四面体网格里，任何相连的四面体的三角形面的组合代表一个曲面。当然，对大部分的这种曲线或曲面是不感兴趣的，因为它们不对准图象特征。

图12突出那些看起来和该地震图象中的断层最对准的三角形边。简单地选择图中示出的这些边，因为它们下面的所有图象采样值低于规定的阈值(-0.2)。未尝试促成各条边之间的连接。尽管如此，许多边是连接的并且的确形成连续曲线。

空间填充网格要比该例子中所采用的简单定阈值更有利于精巧的曲线或曲面提取算法。作为本发明的一个后续，期待开发新的提取算法。和目前通常使用的提取(或分割)算法不同，这些新的算法会产生保证和空间填充网格一致的曲线或曲面。

替代的点阵和网格

在前面的例子中，图象势场定义成向图象特征吸引原子。在一些情况中，希望建立从图象特征推斥原子的伪规则点阵。

可以利用优化算法5得到这样的点阵，该算法带有相同的原子势函数，但另外带有一个沿着图象特征达到最大值(例如1)并且离开图象特征达到最小值(例如0)的图象势场。图13示出这样的对该地震图象优化的结果。这种优化点阵中的原子趋于排在图象特征的旁边(而不是在它们的上面)。

图13还示出从该优化点阵产生的Voronoi网格。Voronoi网格是德洛内三角剖分的对偶。这样，各原子位于网格元素内而不是位于网格元素的顶点处。由于该优化点阵中的原子趋于排在图象特征的旁边，Voronoi网格元素的边界趋于排在这些特征上。

Voronoi网格具有各种各样的应用包括流体流动模拟。Voronoi网格通常导致对偏微分方程的有限差分求解，而单形的三角形/四面体网格导致有限元求解。目前广泛应用这二种求解。在这二种类型中，希望使网格元素边界对准图象特征，因为这些特征常常与物理特性上的不连续性对应。

更一般地，可能希望建立一个伪规则点阵，其中原子被图象中一些特征吸引并且被其它特征推斥。例如，按局部密度和人口密度对应的正常模式配置通信收发机是有好处的。还可能希望把这些收发机放在高海拔的区域(例如山脊的顶上)并且避免放在低海拔的区域(例如峡谷或河谷)。这样，图象势场可包括各最大值和各最小值以及中间平稳段。

其它例子

本发明的方法还可应用于医学图象。图14示出人头部的磁共振图象(MRI)。可以从国家医学实验室免费地得到作为可视人体工程的一部分的该数字图象。

图15示出为增强原始MRI中边缘、不连续性而处理后的图象。采用了简单和周知的Prewitt边缘增强算法。可以把该边缘增强后的图象用作为本发明方法中的图象势场b(x)。

图16示出为对准图15的图象而优化的网格。对该地震图象，通过平滑图象计算非常量的标称距离函数d(x)。标称距离值的范围从最小d＝4到最大d＝12。其图象标量因子为β＝0.4。

图17突出该网格中和图象中的各边对得最准的线段。利用与用于该地震图象的一样的简单定阈值算法选择示出的这些线段。

图18、19、20和21提供另一个医学成像例子。和头部图象一样，可以从国家医学实验室免费得到人的躯干的图象。除了标称距离值的范围从最小d＝5到最大d＝15之外，对该图象的图象处理以及点阵优化是和对头部图象的处理和优化一样的。

计算机系统、存储器媒体以及一种方法实施例

注意，本发明的用于生成考虑数字图象中特征的伪规则点阵的方法可以在各种各样的计算机系统中的任一种上实现，例如在台式计算机、小型计算机、工作站、多处理机系统、各种类型的并行处理机、分布式计算网络、等等中实现。本发明的方法可以在一个或多个存储在诸如CD-ROM、磁盘、磁泡存储器、半导体存储器(例如各种类型的RAM或ROM中的任一种)的各种存储器的任一种上的软件程序或模块中实现。另外，可以在各种各样的诸如光纤、金属线、自由空间和/或通过任何网络如因特网和/或PSTN(公共交换电话网)中的任一种载体媒体上发送上述一个或多个软件程序和/或它们产生的结果。

图22示出一种可通过操作来实现本发明的点阵生成方法的计算机系统的实施例500。计算机系统500可包括一个处理器502，存储器(例如随机存取存储器506或/或非易失性存储器部件504)，一个或多个输入部件508，一个或多个显示部件510，以及一个或多个接口部件512。这些组成部分子系统可以根据各种配置中的任一种互连。非易失性存储器部件504可包括诸如磁带机、磁盘机、半导体ROM或EEPROM等部件。输入部件508可包括诸如键盘、鼠标、数字化垫、跟踪球、触敏垫和/或光笔等部件。显示部件510可包括诸如监视器、投影仪、头戴式显示器等部件。接口部件512可配置成从一个或多个采集部件和/或通过网络从一个或多个远程计算机或存储部件获取数字图象。

取决将成象的目标或处理的类型可使用各种采集部件中的任何部件。该一个或多个采集部件可以感测各种形式的机械能(例如声能，位移和/或应力/应变)和/或电磁能(例如光能，无线电波能，电流和/或电压)中的任何形式。

处理器502可配置成从RAM 506和/或非易失性存储器部件504读程序指令和/或数据，并且把计算结果存入RAM 506和/或非易失性存储器部件504。程序指令指导处理器502根据本文中所说明的各方法实施例的任何组合来操作输入的图象。可以通过各种手段的任一种向计算机系统500提供输入图象。例如，可以利用一个或多个接口部件512把输入图象采集到非易失性存储器504和/或RAM 506中。作为另一个例子，可以通过装在非易失性存储器部件504上/中的盘或带等存储媒介向计算机系统500提供输入图象。在这种情况下，可以通过计算机系统500或通过其他计算机系统预先把输入图象记录到存储器媒体上。

注意，输入图象不必必须是通过采集部件得到的原始传感器数据。例如，输入图象可以是对一组原始传感器数据的一次或多次预处理操作的结果。该一次或多次预处理操作可以由计算机系统500和/或一个或多个其它计算机完成。此外，输入图象可以完全独立于传感器数据，例如在通过利用CAD软件包由设计者产生图象的情况中。

图23示出一种用来产生伪规则的并且考虑数字图象中的特征的点阵的方法的一个实施例600。在步骤605中，可把图象(例如目标和/或过程的图象)采集到可本地访问的存储器，例如计算机系统500的非易失性存储器504和/或RAM 506中。在步骤610中，可预处理该图象以便增强或者揭示该图象中感兴趣的特征。如果在步骤605中得到的图象里这些感兴趣的特征是足够清楚的，则可以省掉步骤610。(如前面所述，步骤605中得到的图象可能已经受过预处理操作。)

本发明预期在步骤610的预处理中使用任何所需的处理技术或者组合处理技术。例如，图象可能受到边缘检测、卷积滤波、傅里叶变换、非线性滤波、阈值处理、等等。该预处理可产生象该数字图象那样采样的图象势场。另外，可对该图象进行预处理操作以产生标称距离函数d(x)。预处理可以是自动的或者是响应用户输入进行的。例如，用户可以增强图象中感兴趣的特征和/或区域。

在步骤615中，可以在该图象覆盖的空间中初始化与该标称距离函数d(x)相一致的点阵。可以通过算法3或算法4的方法初始化该点阵。在一替代实施例中，可以通过执行若干次标定参数β为0的点阵优化迭代从诸如矩形点阵或等边点阵的廉价点阵来产生和该标称距离函数相一致的初始点阵。但是，为了避免撕裂该点阵，可以把该低廉价点阵中的原子间距离设置成等于该标称距离函数的最小值。在另一替代实施例中，可以通过在标定参数β的值不为零并且用距离函数d(x)代替图象势函数b(x)的情况下执行若干次点阵优化迭代来从廉价点阵产生初始点阵。

在步骤620，通过算法5优化初始点阵。优化后的点阵可用于各种应用。例如，可三角形化该点阵以产生一个空间填充网格，并且该网格可用于执行基于网格的模拟(例如储层模拟)、图象编码算法、信号分析方法、等等。在一些实施例中，优化后的点阵可不必介入网格生成步骤由应用来使用。

结论、衍生和范围

上面提供的各例展示本发明方法在产生伪规则的并且响应于数字图象中的特征的点阵中的应用性。优化后的点阵可以具有和图象中的空间变化细节程度相符的可变密度。

这些例子还示出这种点阵的三角形化产生良好对准图象特征的空间填充网格。该网格中的高度规则元素使它成为后继的计算，例如对流体流动的模拟的适用框架。

点阵框架和图象的组合可以充当改进型特征提取算法的基础。

尽管上面示出的例子展示二维图象，本发明方法同样良好地应用于三维图象。在说明中当必要时突出了二维和三维方程及算法之间的差别。本发明方法自然地推广到N维空间，其中N是大于零的任何整数。

各例中的图象是从矩形二维网格上均匀采样的。但是，预期可有各种各样的以各种替代图象拓扑运行本方法的实施列。例如，可以在球的表面上采样二维图象。从而，在一实施例中，本发明方法可以产生在这种球面上的、和该球面上定义的距离函数一致的、并且还和该球面上展示的图象特征对准的伪规则点阵。可以预期不同维数下的各种拓扑。

本发明的方法是按最小化问题提出的。一个明显的数字事实是，函数f的最小化等同于它的负数-f的最大化。从而，预期存在替代实施例，其中通过使由点阵中的点(即原子)之间的原子间距离第一函数和点阵点的位置第二函数所组合成的复合函数为最大进行点阵优化操作。

尽管上面的说明带有许多特定细节，但这些细节不应构成对本发明的限制。相反，它们只是对本发明的当前的各优选实施例提供示意说明，例如，简单的多项式成对势函数可以由另一个具有类似性质的函数代替。类似地，除了上面给出的二种算法外，可以利用各种算法初始化原子点阵。

本发明的范围应由附后权利要求书和其法律等同物确定，而不是由前面提供的实施例确定。

Claims (58)

1.一种计算机实现的用于产生考虑N维数字图象中的特征的点的N维点阵的方法，其中N是大于零的整数，该方法包括：

在通过所述图象采样的空间子集中初始化所述点的所述点阵；

通过移动所述点以使所述点的空间坐标的一个复合函数取极值来优化所述点阵；

其中所述复合函数是所述点之间的成对距离的第一函数和所述点附近的所述图象采样值的第二函数的加权组合；以及

其中所述优化的点阵可用于产生一个用来模拟与所述图象相关的某过程的网格。

2.如权利要求1的方法，其中所述加权组合的加权被选择成平衡所述点阵的所需的规则性程度和所述点阵对所述图象中的所述特征的所需的响应性程度。

3.如权利要求1的方法，还包括从该点阵产生一个空间填充网格。

4.如权利要求3的方法，其中该网格包括N维网格元素，该方法还包括识别所述网格元素子集的边界组合，并且其中所述边界组合对应于所述图象中感兴趣区域之间的边界。

5.如权利要求1的方法，其中所述初始化使用一个距离函数，该距离函数用于为所述点阵中的每个点规定从所述点阵中的所述点到所述点阵中靠近所述点的其它点的优选距离。

6.如权利要求5的方法，其中所述初始化点阵包括：

从一队列中读出第一试用点阵点；

检查所述第一试用点阵点周围的第一邻域以判定在该第一邻域内是否不包含任何处于该点阵当前状态下的点阵点；以及

如果该第一邻域不包含任何处于该点阵当前状态下的点阵点，则执行包含如下步骤的更新操作：

把该第一试用点阵点添加到该点阵上，以及

把新的试用点阵点增补到该队列中，其中对所述新的试用点阵点分配坐标位置，该坐标位置与相对于该第一试用点阵点以及该距离函数为完全规则的点阵对应。

7.如权利要求6的方法，还包括：

在该队列中存储一个或多个种子点，以及

反复地进行所述读出、检查和对该更新操作的条件执行，直至填充所述空间的该子集。

8.如权利要求5的方法，其中所述初始化该点阵包括随机地在通过所述图象采样的位置上对所述点阵添加点，其中在所述采样位置之一上对所述点阵添加一个点的概率由在该一个采样位置上计算的距离函数决定。

9.如权利要求1的方法，其中所述优化所述点的点阵包括计算该第一函数的值，其中所述计算该第一函数的值包括：

在一点阵点的邻域中的各图象采样位置上计算函数范式以得到本地函数值；

在一个表示原子势场的数组中向对应的采样添加本地函数值；

为该点阵的每个点重复进行所述计算和所述添加；

根据表示原子势场的数组中的一个或多个采样，计算和每个点阵点对应的原子势场值；以及

为每个点阵点添加原子势场值以确定该第一函数的值。

10.如权利要求9的方法，其中该函数范式包括与该点阵点的距离的函数除以在该点阵点的邻域中计算的标称距离函数。

11.如权利要求1的方法，其中所述优化所述点的点阵还包括计算所述复合函数相对于所述点的每个所述空间坐标的导数，以及根据所述导数确定对所述点的所述空间坐标的一组更新值。

12.如权利要求1的方法，其中该图象沿着所述图象中的所述特征的至少一个子集得到最小值。

13.如权利要求1的方法，其中该图象沿着所述图象中的所述特征的至少一个子集得到最大值。

14.一种计算机实现的用于生成对象的表达以供模拟的方法，该方法包括：

获得该对象的图象；

在通过该图象采样的空间中初始化点的点阵；以及

一次或多次调整这些点在该空间中的位置以使这些点的空间坐标的一个复合函数取极值；

其中该复合函数是点之间的成对距离的第一函数和点附近的所述图象的采样值的第二函数的加权组合；以及

其中所述一次或多次调整位置产生一个最终点阵，该点阵用于生成一个用来模拟与所述对象相关的某过程的网格。

15.如权利要求14的方法，其中所述对象包括一种或多种流体的地下储层，并且其中模拟包括基于所述网格对所述流体的流动的模拟。

16.如权利要求14的方法，其中所述成对距离之一对应于一对所述点阵点，并且其中成对距离包括(a)该对点阵中各点之间的距离与(b)在该对点阵点的邻域中算出的标称距离函数值的比。

17.一种计算机实现的用于从图象生成点阵的方法，该方法包括：

在通过该图象采样的空间的至少一部分中初始化一个点阵；以及

在该空间中一次或多次调整该点阵中的点的位置以向总势能的极值收敛，其中所述总势能包括该点阵的原子势能和该点阵的图象势能的组合；

其中该图象势能是从该图象和点阵点的位置计算的；以及

其中所述一次或多次调整位置产生一个用于编码与该图象相关的信息的最终点阵。

18.如权利要求17的方法，其中该图象势能包括图象势场值的和，并且其中每个所述图象势场值对应于在一个点阵点处对图象势场求得的值。

19.如权利要求18的方法，其中该图象包括感兴趣的特征，该方法还包括通过处理该图象产生该图象势场，以便该图象势场沿着一个或多个所述感兴趣的特征得到第一势值，以及离开所述感兴趣的特征得到第二势值。

20.如权利要求19的方法，其中第一势值是最大势值。

21.如权利要求19的方法，其中第一势值是最小势值。

22.如权利要求17的方法，其中该原子势能包括各对所述点阵点的原子势的和。

23.如权利要求22的方法，其中所述各对所述点阵点中的一对的原子势是该对点之间的归一化距离的函数，并且其中该归一化距离等于该对点之间的未归一化距离与一标称距离函数的本地值的比。

24.如权利要求23的方法，还包括：

平滑该图象以产生一个平滑图象；以及

基于该平滑图象的值对该标称距离函数分配值。

25.如权利要求17的方法，还包括从该最终点阵生成一个网格，其中该网格能用于编码与该图象相关的信息。

26.如权利要求25的方法，还包括利用该网格对某物理过程进行模拟。

27.如权利要求17的方法，其中所述组合包括线性组合αA+βB，其中A是原子势能，B是图象势能，并且其中系数α和β控制所述点阵相对于图象中的特征的位置灵敏性程度以及所述点阵的规则性程度。

28.如权利要求17的方法，还包括：

获得一系列图象；以及

在所述系列图象的每幅所述图象上进行所述初始化和调整以产生一系列对应的最终点阵。

29.如权利要求17的方法，还包括：

获得相继的图象；

利用该最终点阵对相继图象中的差异编码以生成图象编码；以及

在传输媒体上发送该图象编码。

30.一种计算机实现的用于确定适用于模拟的点阵的方法，该方法包括：

在图象上运算以产生图象势场，其中所述图象势场沿着所述图象中的感兴趣特征得到极值势值；

在和所述图象对应的空间的一个区域内初始化点阵；以及

在所述区域内反复调整所述点阵的位置以使总势能取极值；

其中所述总势能包括原子势能和图象势能的组合；

其中该原子势能包括各对所述点阵点的原子势的和；

其中该图象势能包括图象势的和，其中每个所述图象势对应于在一个点阵点处对图象势场求得的值。

31.如权利要求30的方法，其中所述在第一图象上运算以产生图象势场包括检测所述第一图象中的边，并且其中所述感兴趣的特征包括所述第一图象中的边。

32.如权利要求30的方法，其中所述在图象上运算以产生图象势场包括：

识别图象中感兴趣的特征；

向和所述感兴趣的特征的至少一个子集对应的第一组空间位置处的图象势场分配第一常数值；以及

向不和所述感兴趣的特征对应的第二组空间位置处的图象势场分配第二常数值。

33.如权利要求32的方法，其中所述第一常数值是所述图象势场的最大值。

34.如权利要求32的方法，其中所述第一常数值是所述图象势场的最小值。

35.如权利要求32的方法，其中所述在图象上运算以产生图象势场还包括对通过所述分配第一常数值和第二常数值产生的图象势场的第一版本进行平滑以得到该图象势场的第二版本。

36.如权利要求35的方法，其中该图象势场的所述第二版本用于为所述反复调整位置的第一组重复计算图象势，并且其中该图象势函数的所述第一版本用于为所述反复调整位置的第二组重复计算原子图象势。

37.如权利要求30的方法，还包括：

在所述空间中随机地移动所述点阵点的位置；执行所述随机移动和所述反复调整的多次迭代以得到一系列总势能的局部极值。

38.如权利要求30的方法，其中所述初始化使用一个用于规定从所述点阵中的每个点到所述点附近的其它点的优选距离的距离函数，并且其中所述优选距离是所述区域内的位置的函数。

39.如权利要求38的方法，其中所述初始化在空间的该区域里的点阵包括：

从一队列中读出第一试用点阵点；

检查所述第一试用点阵点周围的第一邻域以判定该第一邻域是否不含有所述点阵的点；以及

若该第一邻域不含有所述点阵的点，则执行以下更新操作：

把该第一试用点阵点添加到该点阵上，以及

对该队列添加新的试用点阵点，其中对所述新的试用点阵点分配坐标位置，该坐标位置与相对于该第一试用点阵点和该距离函数为完全规则的点阵对应。

40.如权利要求39的方法，还包括：

在该队列中存储一个或多个种子点；和

反复地执行所述读出、检查并且条件地执行更新操作，直至填充所述空间的子集。

41.如权利要求39的方法，其中围绕所述第一试用点阵点的第一邻域是一个N维的球，该球的直径和在该第一试用点阵点处计算的距离函数成比例，其中N是和所述图象对应的所述空间的维数。

42.如权利要求30的方法，其中所述初始化在空间的该区域中的点阵包括：

根据一距离函数计算该区域内一候选位置处的本地概率值；

生成一随机数值；

判定该随机数值是否满足相对于该本地概率值的一个不等式条件；以及

若该随机数值满足该相对于本地概率值的不等式条件，则条件地把该候选位置添加到该点阵上。

43.如权利要求42的方法，其中所述初始化点阵还包括：

移动该候选位置经过一组通过所述第一图象采样的图象位置；以及

对每个所述图象位置进行所述计算、生成、判定和条件地添加。

44.如权利要求30的方法，其中该组合的形式为(1-β)A+βB，其中A为原子势能，B为图象势能，并且其中参数β取自零至1的范围内的任何值。

45.如权利要求30的方法，还包括三角形化该点阵以从该点阵生成一个网格，其中所述网格包括网格元素的组合，并且其中网格元素的顶点是各点阵点。

46.如权利要求45的方法，其中所述三角形化点阵包括对该点阵进行德洛内三角剖分。

47.权利要求45的方法，还包括识别所述网格元素的靠近一个或多个所述感兴趣的特征的边界组合。

48.如权利要求45的方法，还包括在该网格上进行模拟，其中该模拟产生表示物理系统特性的输出。

49.如权利要求48的方法，其中该物理系统包括地下储层，该地下储层包含一种或多种流体。

50.一种被配置成确定适用于模拟的点阵的计算机系统，该计算机系统包括：

一个处理器；

一个和所述处理器耦接并且被配置成存储程序指令的存储器；以及

一个用于向所述存储器提供图象的输入端口；

其中所述处理器被配置成从所述存储器读出并执行程序指令，其中，响应对程序指令的所述执行，该处理器可通过操作完成：

在该图象上运算以生成图象势场，其中所述图象势场沿着所述图象中感兴趣的特征得到一个或多个极值势值；

在和所述图象对应的空间区域中初始化点阵；以及

在所述区域内反复调整所述点阵的位置以使总势能取极值；

其中所述总势能包括原子势能和图象势能的组合；

其中该原子势能包括各对所述点阵点的原子势的和；

其中该图象势能包括图象势的和，以及

其中每个所述图象势和在一个点阵点处对图象势场求得的值对应。

51.一种配置成存储程序指令的计算机可读存储媒体，其中可执行程序指令以实现：

初始化在通过图象采样的空间子集中的点的点阵；以及

通过移动所述点以使所述点的空间坐标的一个复合函数取极值来优化所述点阵；

其中所述复合函数是所述点之间的成对距离的第一函数和所述点附近的所述图象采样值的第二函数的加权组合；以及

其中所述优化的点阵用于生成一个用来模拟和所述图象相关的某过程的网格。

52.如权利要求51的存储媒体，其中还可执行该程序指令以从该点阵生成一个空间填充网格。

53.如权利要求52的存储媒体，其中该网格包括N维网格元素，其中还可执行该程序指令以识别所述网格元素一子集的边界组合，并且其中所述边界组合对应于所述图象中感兴趣区域之间的边界。

54.如权利要求51的存储媒体，其中所述初始化使用一个用于规定从所述点阵的每个点到所述点附近的其它点的优选距离的距离函数；并且其中所述优选距离是所述空间中的位置的函数。

55.如权利要求54的存储媒体，其中所述初始化点阵包括：

从一队列中读出第一试用点阵点；

检查所述第一试用点阵点周围的第一邻域以判定该第一邻域是否不含有所述点阵的点；以及

若该第一邻域不含有所述点阵的点，则执行以下更新操作：

把该第一试用点阵点添加到该点阵上，以及

对该队列添加新的试用点阵点，其中向所述新的试用点阵点分配坐标位置，该坐标位置与相对于该第一试用点阵点和该距离函数为完全规则的点阵对应。

56.如权利要求54的存储媒体，其中所述初始化该点阵包括在通过所述图象采样的位置上随机地对所述点阵添加点，其中在所述采样位置之一处对所述点阵添加一个点的概率由在该一个所述采样位置上计算的距离函数确定。

57.如权利要求51的存储媒体，其中所述优化所述点的点阵包括反复计算该第一函数的值，其中所述计算该第一函数的值包括：

在一个点阵点的邻域中的各图象采样位置上计算函数范式以得到各本地函数值；

在一个表示原子势场的数组中向对应的采样添加本地函数值；

对该点阵的每个点重复进行所述计算和所述添加；

根据表示原子势场的该数组中的一个或多个采样，计算和每个点阵点对应的原子势场值；以及

为每个点阵点添加原子势场值以确定该第一函数的值。

58.如权利要求57的存储媒体，其中该函数范式包括对该点阵点的距离的函数除以在该点阵点的邻域计算的标称距离函数。

Images