CN103679650A - 岩心三维图像修复方法 - Google Patents

Abstract

本发明公开了一种岩心三维图像修复方法，包括以下步骤：（1）提取缺失信息的岩心三维图像即原始三维图像中待修复区域的初始边界，并设定模板大小；（2）对原始三维图像的已知区域进行区域划分；（3）计算原始三维图像的边界上所有体素的边界模板的优先级值；（4）搜索最大优先级值的边界模板，并保存相应的边界体素的坐标；（5）判断边界体素的区域类型，并设置边界模板的匹配块的全局搜索区域；（6）在相应的范围内搜索最优匹配块，并进行相应体素的灰度值填充；（7）更新或设置各参数信息；（8）重复步骤（4）～（7）直至填充完待修复区域。通过本发明修复岩心三维图像复原度极高，能够保持真实三维图像的结构和细节。

Description

岩心三维图像修复方法

技术领域

本发明涉及一种图像修复方法，尤其涉及一种岩心三维图像修复方法，属于三维图像修复技术领域。

背景技术

在石油地质分析过程中，通常需要获得三维岩心结构来定量地研究渗流的微观机理，分析储集层渗流性质和运移规律。随着CT（computed tomography，电子计算机X射线断层扫描技术）分辨率的提高，CT机已较多应用到岩石三维成像中，但CT的分辨率与样本尺寸是相矛盾的。为了得到准确的分析结果，需要获取高分辨率的孔隙结构图像，但扫描的样本尺寸受到很大限制，使岩样的代表性有所欠缺。因此很难得到可供分析的岩心高分辨率三维孔隙结构图像。针对这一问题，目前的主要解决方案是以二维薄片为基础进行三维重建，得到具有与薄片特性相近的三维图像。也有基于三维图像的三维重建，其原理与基于二维图像的三维重建相似。我们统一将此种基于局部特性重建出三维图像的方法称为三维数学建模。然而三维数学建模是基于局部特性，重建结果缺乏一定的实际代表性，同时存在随机性，且对于多相或非均质的岩心三维图像的重建存在较大的困难，重建效果与实际三维情况存在较大差距。

将大尺寸岩心样本切割为多个小样本，用CT机分别对小样本进行扫描，将获取的小样本图像拼接成较大的三维图像，能有效地解决这一问题。然而大尺寸样本在切割的过程中存在三维图像信息的缺失，且缺损信息较大，在小样本图像拼接成大样本图像的过程中需要恢复缺损的信息，实现大范围信息缺失的三维图像修复。因此本发明的目的就是提供一种可修复大范围信息缺失的岩心三维图像的方法，保持原有岩心三维图像的形态结构。

目前对于三维图像修复的研究还较少，大部分还是针对二维图像的修复。图像修复就是对图像的信息缺损区域进行信息填充的过程,其目的就是为了对有信息缺损的图像进行恢复,并且要使观察者无法察觉到图像曾经缺损或已被修复。目前的三维图像修复，主要是针对于医学MRI、CT图像小范围信息缺失的三维图像修复以及视频的修复。Satoru Morita（Three Dimensional ImageInpainting[C].In proceeding of:Advances in Natural Computation,SecondInternational Conference,Xi'an,China,2006,4222:752-761）将基于偏微分方程PDE的方法拓展为三维修复算法，并将该算法应用于医学图像和视频的修复，但是该方法只是针对小范围信息缺失的三维图像修复，且主要是对于结构的修复较好，而对于纹理部分却存在模糊现象。由于视频中帧内的相关度往往大于帧间的相关度，帧内信息较为重要，所以对于视频的修复主要是逐帧进行，较少考虑帧之间的关系。而对于岩心三维图像，其在各个方向上的相关性都一样，需同等考虑各个方向的信息。

发明内容

本发明的目的就在于为了解决上述问题而提供一种复原度极高的缺失信息的岩心三维图像修复方法。

本发明通过以下技术方案来实现上述目的：

1、本发明所述岩心三维图像修复方法包括以下步骤：

（1）提取缺失信息的岩心三维图像即原始三维图像中待修复区域Ω的初始边界δΩ⁰，并设定模板大小；

（2）对原始三维图像的已知区域Φ进行区域划分；

（3）计算原始三维图像的边界上所有体素q的边界模板的优先级值，所述体素q指三维数据场中的采样点；

（4）搜索最大优先级值的边界模板Ψ_p，并保存边界体素P的坐标；

（5）判断边界体素P的区域类型，并设置边界模板Ψ_p的匹配块的全局搜索区域；

（6）根据边界体素P的匹配中心因子MCF和最大距离MD在相应的范围内搜索最优匹配块Ψ_q，并进行相应体素的灰度值填充；

（7）更新或设置以下信息：置信度、边界面δΩ、新产生的边界体素P′的匹配中心因子MCF和最大距离MD，以及新产生的边界模板Ψ_p′的优先级值；

（8）重复步骤（4）～（7）直至填充完待修复区域Ω，即

上述方法的基本原理如下：

岩心CT图像具有灰度对比度较低、纹理单一、噪声较多、且结构不规则等特征。岩心的三维结构在三个维度上都保持一定的连续性，且三个方向的信息都同等重要。对于具体某一体素的填充，需同等考虑三个方向上的已知信息。通常岩心三维图像的结构复杂，但在一定的较小的范围内，其结构具有一定的相似性，且一般情况下岩心三维图像含有丰富的结构信息，因此可以根据三维图像已知区域的信息，利用三维图像修复的方法对缺失部分进行填补。在修复岩心三维图像的过程中，重要之处在于需要保持三维结构的连续性和完整性，但同时也得保持三维图像的纹理以及整体的结构特点，从而达到与原始图像的结构和细节相近的目的，且使得无法觉察出图像的缺损。同时，由于切割过程对岩心的三维图像造成了较大范围的信息缺失。因此可以利用基于样例的修复思想，根据三维图像的线性结构（在图像修复领域称为等照度面），对结构和纹理同时进行修复或填充，保持三维图像原有的整体结构和纹理特点。由于我们需要考虑图像各个方向上的信息，其样本块需要设立为三维体。三维空间中存在的往往是等照度面和边界面，而不再简单地以线的形式存在，若要较好地修复三维图像的结构，则修复过程中需通过保持并传播等照度面的趋势来修复三维图像的结构。因此，在我们的三维图像修复算法中，需要获得等照度面的信息以及其对于修复过程的影响方式。利用通常三维图像的结构和纹理都具有连续性这一特征，可以将传统的基于样例的修复方法中的全局匹配搜索限定为局部搜索，大大地减小搜索范围。由于通常的岩心CT序列图像包含了大量的随机噪声且灰度的对比度较低，为了较为准确地得到等照度面的信息，需首先对岩心CT序列图像进行简单的去噪和白平衡等预处理，使得岩心三维图像的结构更加清晰。若原始CT图像的结构较为清晰，则不需要预处理。

具体地，所述步骤（1）中，所述模板为一固定大小的立方体，边长为2winsize+1，winsize的大小根据原始三维图像的结构和纹理复杂度确定；

所述步骤（2）中，将已知区域Φ划分为各相（R₁，R₂，…R_n）和混合区域R₀，有n相（D₁，D₂，…D_n）的岩心，划分的区域种类为n+1类，划分规则为：若 $\begin{array}{cc}D\left(l\right)\∈D_m& \∀l\∈{\mathrm{\Ψ}}_q,\end{array}$ 则q∈R_m；若

D(k)≠D(l)，则q∈R₀；其中，q为某一体素，Ψ_q为样本块，R₀为混合区域，D_m为某一相，R_m为与D_m相对应的区域；

所述步骤（3）中，原始三维图像中的每一个体素q都存在一个灰度值和置信度C(q)，当体素没有被填充时灰度值为空，体素q的置信度C(q)表示q的灰度值的可信度，初始化为： $\begin{array}{cc}C\left(q\right)=0& \∀q\∈\mathrm{\Ω}\end{array}$ 和 $\begin{array}{cc}C\left(q\right)=1& \∀q\∈\mathrm{\Φ}\end{array};$

每一个边界体素P的模板即边界模板Ψ_p都有一个优先级值P(p)，其定义为：

P(p)＝C(p)D(p)

上式中，C(p)为边界模板Ψ_p的置信度，用于衡量体素P周围模块范围内可靠信息的数量，D(p)为数据项，表示等照度面在边界面δΩ处的强度，C(p)和D(p)的计算公式为：

$C\left(p\right)=\frac{{\mathrm{\Σ}}_{q\∈{\mathrm{\Ψ}}_p\∩\stackrel{\‾}{\mathrm{\Ω}}}C\left(q\right)}{\left|{\mathrm{\Ψ}}_p\right|}$

$D\left(p\right)=\frac{|\▿I_p^{\⊥}|\·\left|\sin \mathrm{\θ}\right|}{\mathrm{\α}}$

其中，

表示边界模板Ψ_p中已被填充的体素，|Ψ_p|表示边界模板Ψ_p的体积，α为归一化因子，

为边界体素P的梯度向量，该梯度包含方向和大小，

为等照度面，θ为边界体素P处的边界面δΩ与等照度面

之间的夹角，也等于梯度向量

与边界面δΩ的法线n_p之间的夹角；由于等照度面

的强度

等于梯度

的大小所以D(p)的计算公式为：

$D\left(p\right)=\frac{|\▿I_p|\·\left|\sin \mathrm{\θ}\right|}{\mathrm{\α}}$

$\sin \mathrm{\θ}=\frac{\▿I_p\·n_p}{|\▿I_p|*\left|n_p\right|}$

所述步骤（4）中，边界体素P是边界δΩ上的体素，边界δΩ为所有边界体素P的集合；

所述步骤（5）中，边界体素P的区域类型的判定方法为：根据模板范围内所有已知体素的相位判定该边界体素P属于哪一相或混合区域；根据边界体素P的区域类型设置边界模板Ψ_p的匹配块的全局搜索区域；

所述步骤（6）中，所述匹配中心因子MCF是边界体素P相应的匹配块局部搜索区域的中心体素；所述最大距离MD表示在边界体素P的匹配中心因子MCF不为空的情况下，局部搜索得到的最优匹配块Ψ_q与其边界模板Ψ_p的距离的最大范围；当局部搜索得到的匹配块与相应的边界模板Ψ_p的距离大于其最大距离MD时，该匹配块将不能作为其最优匹配块，需重新全局搜索其最优匹配块；某一边界体素P′的最大距离MD表示为：

MD_P′＝μ*d(Ψ_p,Ψ_q)

其中μ为比例参数，通常设定为1及其较为临近的值，边界体素P′为填充边界模板Ψ_p后新产生的边界体素，根据边界模板Ψ_p和其最优匹配块Ψ_q两个模板之间的距离，设置边界体素P′的最大距离MD；在设置新产生的边界体素P′的匹配中心因子MCF时，同步设置其最大距离MD；匹配准则为：

${\mathrm{\Ψ}}_q=\arg \underset{{\mathrm{\Ψ}}_{\hat{q}}\∈\mathrm{\Φ}}{\min }d({\mathrm{\Ψ}}_p,{\mathrm{\Ψ}}_{\hat{q}})$

其中两模板的距离

为边界模板Ψ_p和样本块

对应位置已知体素的灰度值之差的平方和；

得到最优匹配块Ψ_q后，用最优匹配块Ψ_q中已知的体素的灰度值填充边界模板Ψ_p相应位置上待修复的体素，此时已知区域Φ和待修复区域Ω也都相应地发生了改变；

所述步骤（7）中，边界模板Ψ_p中新填充的体素的置信度设置为边界模板Ψ_p的置信度，即是：

$\begin{array}{cc}C\left(q\right)=C\left(p\right)& \∀q\∈{\mathrm{\Ψ}}_p\∩\mathrm{\Ω}\end{array}$

若设置的模板是边长为2winsize+1大小的立方体，那么在搜索新的边界体素时，只需要以刚填充的边界体素P为中心，在边长为2(winsize+1)+1的立方体领域内搜索新的边界体素，并去除已填充的边界体素，而不需要重新搜索整个边界；对于填充边界模板Ψ_p后新产生的边界体素，需要设置其匹配中心因子MCF和最大距离MD；在更新边界模块Ψ_p的置信度和临时优先级值时，只需要对以边界体素P为中心、边长为4winsize+1大小的立方体邻域内的所有边界体素进行重新计算。

本发明的有益效果在于：

本发明同等考虑了三维图像各个方向的信息，能够修复由于切割造成的较大范围信息缺失的岩心三维图像，利用三维空间中的等照度面信息能够较好地修复三维图像的结构和细节，且不存在边界模糊现象，同时结合了局部搜索和全局搜索的优点，大大地减少了匹配过程中的计算量，而且保证了匹配的正确性，修复后能够得到结构和纹理连续、完整的岩心三维图像，保持了真实三维图像的结构和细节，使得获得的岩心三维图像与实际三维结构极为相似，解决了以往基于局部特性进行三维重建导致获得的三维图像缺乏一定实际代表性以及重建效果与实际三维结构存在较大差距的问题。

附图说明

图1-1是本发明实施例中无信息缺失的真实岩心CT序列图像之一；

图1-2是本发明实施例中无信息缺失的真实岩心CT序列图像之二；

图1-3是本发明实施例中无信息缺失的真实岩心CT序列图像之三；

图2-1是本发明实施例中存在信息缺失的原始岩心CT序列图像之一；

图2-2是本发明实施例中存在信息缺失的原始岩心CT序列图像之二；

图2-3是本发明实施例中存在信息缺失的原始岩心CT序列图像之三；

图3-1为图1-1中的真实岩心CT图像经过去噪和白平衡预处理后的图像；

图3-2为图1-2中的真实岩心CT图像经过去噪和白平衡预处理后的图像；

图3-3为图1-3中的真实岩心CT图像经过去噪和白平衡预处理后的图像；

图4-1为图2-1中的原始岩心CT图像经过去噪和白平衡预处理后的图像；

图4-2为图2-2中的原始岩心CT图像经过去噪和白平衡预处理后的图像；

图4-3为图2-3中的原始岩心CT图像经过去噪和白平衡预处理后的图像；

图5-1为图4-1中预处理后的原始岩心CT图像经过修复后的图像；

图5-2为图4-2中预处理后的原始岩心CT图像经过修复后的图像；

图5-3为图4-3中预处理后的原始岩心CT图像经过修复后的图像；

图6-1为与图3-1、图3-2、图3-3对应的经过去噪和白平衡预处理后的真实三维图像；

图6-2为与图4-1、图4-2、图4-3对应的经过去噪和白平衡预处理后的原始三维图像；

图6-3为与图5-1、图5-2、图5-3对应的经过修复后的岩心三维图像。

具体实施方式

下面结合具体实施例和附图对本发明作进一步说明：

实施例：

为了使本发明所述修复方法更加便于理解和接近于真实应用，下面从原始岩心样本的切割开始到图像修复完成为止进行整个流程的整体说明，其中包括本发明的核心修复方法：

（1）记录下原始岩心样本的尺寸，对其进行切割。

（2）使用CT机依次扫描切割后的小岩心样本，得到128张无信息缺失的真实岩心CT序列图像并存储在计算机中，所述128张真实岩心CT序列图像以图1-1、图1-2、图1-3示出的三张连续图像为代表，如图1-1、图1-2和图1-3所示，图片大小为128*128，同时记录下各样本相对于原始岩心样本所处的位置。

（3）使用岩心三维图像修复软件读取步骤（2）中扫描得到的真实岩心CT序列图像，并设定序列图像中间15*128大小的部分为由于切割造成的缺失部分，并将其灰度值设为100，得到128张存在信息缺失的原始岩心CT序列图像并存储在计算机中，所述128张原始岩心CT序列图像以图2-1、图2-2、图2-3示出的三张连续图像为代表，如图2-1、图2-2和图2-3所示，图片大小为128*128，中间15*128大小的部分为缺失部分。

（4）使用岩心三维图像修复软件读取步骤（2）中得到的真实岩心CT序列图像，并使用岩心三维图像修复软件的批处理功能对其进行去噪和白平衡预处理，增加图像的对比度，使其更加清晰，得到128张对比度更高、更清晰的无信息缺失的真实岩心CT序列图像并存储在计算机中，以图3-1、图3-2、图3-3示出的三张连续图像为代表；对应地，如图6-1所示，将经过去噪和白平衡预处理后的真实岩心CT序列图像组成的真实三维图像也存储在计算机中。

（5）使用岩心三维图像修复软件读取步骤（3）中得到的原始岩心CT序列图像，并使用岩心三维图像修复软件的批处理功能对其进行去噪和白平衡预处理，增加图像的对比度，使其更加清晰，得到128张对比度更高、更清晰的存在信息缺失的原始岩心CT序列图像并存储在计算机中，以图4-1、图4-2、图4-3示出的三张连续图像为代表；对应地，如图6-2所示，将经过去噪和白平衡预处理后的原始岩心CT序列图像组成的原始三维图像也存储在计算机中。

（6）提取步骤（5）得到的原始三维图像中待修复区域Ω的初始边界δΩ⁰，并设定模板大小，该模板为一固定大小的立方体，边长为2winsize+1，winsize的大小根据原始三维图像的结构和纹理复杂度确定，本例设为9*9*9大小的立方体；

（7）对原始三维图像的已知区域Φ进行区域划分：根据岩心纹理单一、相位较少的特征，依据三维图像中的相以及模板大小的邻域内的体素值对原始三维图像的已知区域Φ进行区域的划分，将已知区域Φ划分为各相（R₁，R₂，…R_n）和混合区域R₀，有n相（D₁，D₂，…D_n）的岩心，划分的区域种类为n+1类；

划分规则为：

若 $\begin{array}{cc}D\left(l\right)\∈D_m& \∀l\∈{\mathrm{\Ψ}}_q\end{array},$ 则q∈R_m；若

D(k)≠D(l)，则q∈R₀，其中，q为某一体素，Ψ_q为样本块，R₀为混合区域，D_m为某一相，R_m为与D_m相对应的区域；

具体就是对于某一体素q，以其为中心点，模板大小邻域内即样本块Ψ_q的所有体素若全为某一相D_m（岩石相、孔隙相等），则将体素q划分为这一相对应的区域R_m；若不止一相，则该体素属于混合区域R₀；本实施例中，原始三维图像的已知区域Φ含有岩石相和孔隙相，被划分为3个区域。

（8）计算原始三维图像的边界上所有体素q的边界模板的优先级值，所述体素q指三维数据场中的采样点；原始三维图像中的每一个体素q都存在一个灰度值Gray-level（当体素没有被填充时为空）和置信度C(q)。体素q的置信度C(q)表示q的灰度值的可信度，初始化为： $\begin{array}{cc}C\left(q\right)=0& \∀q\∈\mathrm{\Ω}\end{array}$ 和 $\begin{array}{cc}C\left(q\right)=1& \∀q\∈\mathrm{\Φ}\end{array}.$

边界体素P是边界δΩ上的体素，边界δΩ为所有边界体素P的集合，每一个边界体素P的边界模板Ψ_p都有一个优先级值P(p)，其定义为：

P(p)＝C(p)D(p)

其中C(p)为边界模板的置信度，用于衡量体素P周围模块范围内可靠信息的数量，D(p)为数据项，表示等照度面在边界δΩ处的强度，C(p)和D(p)的计算公式为：

$C\left(p\right)=\frac{{\mathrm{\Σ}}_{q\∈{\mathrm{\Ψ}}_p\∩\stackrel{\‾}{\mathrm{\Ω}}}C\left(q\right)}{\left|{\mathrm{\Ψ}}_p\right|}$

$D\left(p\right)=\frac{|\▿I_p^{\⊥}|\·\left|\sin \mathrm{\θ}\right|}{\mathrm{\α}}$

其中，

表示边界模板Ψ_p中已被填充的体素，|Ψ_p|表示边界模板Ψ_p的体积，α为归一化因子，为边界体素P的梯度向量，该梯度包含方向和大小，

为等照度面，θ为边界体素P处的边界面δΩ与等照度面

之间的夹角，也等于梯度向量

与边界面δΩ的法线n_p之间的夹角。说明：当等照度面的影响相同时，算法会优先选择那些含有可靠信息较多的边界体素P进行填充，对于那些更加靠近目标区域中心的体素将后被填充。

由于等照度面

的强度

等于梯度

的大小

D(p)的计算公式为：

$D\left(p\right)=\frac{|\▿I_p|\·\left|\sin \mathrm{\θ}\right|}{\mathrm{\α}}$

$\sin \mathrm{\θ}=\frac{\▿I_p\·n_p}{|\▿I_p|*\left|n_p\right|}$

当等照度面

的强度越强，且与边界面δΩ越接近垂直状态，数据项D(p)越大，在置信度相同的情况下，算法会优先选择这样的边界体素P进行修复。

（9）搜索最大优先级值的边界模板Ψ_p，并保存边界体素P的坐标。

（10）判断边界体素P的区域类型，并设置边界模板Ψ_p的匹配块的全局搜索区域；

边界体素P的区域类型的判定方法为：根据模板范围内所有已知体素的相位判定该边界体素P属于哪一相或混合区域，该方法与步骤（7）中的区域划分方法是一样的。根据边界体素P的区域类型设置边界模板Ψ_p的匹配块的全局搜索区域GSA（Global Search Area），若该边界体素的区域类型为R_m，则边界模板Ψ_p的GSA为相应的区域R_m，若其区域类型为R₀，则其GSA为整个三维图像。

（11）根据边界体素P的匹配中心因子MCF和最大距离MD在相应的范围内搜索最优匹配块Ψ_q，并进行相应体素的灰度值填充；

基于通常三维图像的结构和纹理都具有连续性，因此采用匹配中心因子MCF来限定匹配块搜索范围。匹配中心因子MCF是边界体素P相应的匹配块局部搜索区域的中心体素，将匹配块搜索限定为局部搜索。当MCF存在时，匹配块局部搜索区域LSA（Local Search Area）为以MCF为中心，边长为两倍模板边长大小（即边长为4winsize+1）的立方体区域，因此，匹配块的搜索区域为GSA∩LSA。初始化时，各边界体素的MCF设为空（NULL）。当边界体素的MCF为NULL时，该边界体素的匹配块搜索范围为该边界体素所对应的GSA。若样本块Ψ_q为边界模板Ψ_p的最优匹配块，体素P′为Ψ_p填充后新产生的边界体素，那么将Ψ_q的中心体素Q设置为边界体素P′的匹配中心因子。依次类推，P′的匹配中心因子为Q，且其最优匹配块Ψ_p′的中心体素为Q′，则修复P′后新产生的边界体素P″的匹配中心因子为Q′。随着匹配中心因子的不断更新，LSA也相应地发生移动，保持了等照度面的延伸趋势。

为了防止修复过程中匹配搜索陷入某一个局部中，以及在突变情况下（如尖锐的角，结构的消失和出现）局部搜索造成的错误，算法中设置了最大距离作为跳出局部搜索的条件限制。最大距离MD（Maximum Distance）表示在边界体素P的匹配中心因子不为空的情况下，局部搜索得到的最优匹配块与其边界模板的距离的最大范围。即是当局部搜索得到的最优匹配块与相应的边界模板的距离大于其最大距离MD时，表明此处的图像结构极可能发生了突变，该匹配块将不能作为其最优匹配块，需在其对应的GSA内全局搜索其最优匹配块。边界体素P′的最大距离表示为：

MD_P′＝μ*d(Ψ_p,Ψ_q)

其中μ为比例参数，通常设定为1及其较为临近的值，边界体素P′为填充边界模板Ψ_p后新产生的边界体素，根据边界模板Ψ_p和其最优匹配块Ψ_q两个模板之间的距离，设置边界体素P′的最大距离MD；在设置新产生的边界体素P′的匹配中心因子MCF时，同步设置其最大距离MD；匹配准则为：

${\mathrm{\Ψ}}_q=\arg \underset{{\mathrm{\Ψ}}_{\hat{q}}\∈\mathrm{\Φ}}{\min }d({\mathrm{\Ψ}}_p,{\mathrm{\Ψ}}_{\hat{q}})$

其中两模板的距离

为边界模板Ψ_p和样本块

对应位置已知体素的灰度值之差的平方和；

得到最优匹配块Ψ_q后，用最优匹配块Ψ_q中已知的体素的灰度值填充边界模板Ψ_p相应位置上待修复的体素，此时已知区域Φ和待修复区域Ω也都相应地发生了改变；

（12）更新或设置以下信息：置信度、边界面δΩ、新产生的边界体素P′的匹配中心因子MCF和最大距离MD，以及新产生的边界模板Ψ_p′的优先级值；边界模板Ψ_p中新填充的体素q的置信度设置为边界模板Ψ_p的置信度，即是：

$\begin{array}{cc}C\left(q\right)=C\left(p\right)& \∀q\∈{\mathrm{\Ψ}}_p\∩\mathrm{\Ω}\end{array}$

若设置的模板是边长为2winsize+1大小的立方体，那么在搜索新的边界体素时，只需要以刚填充的边界体素P为中心，在边长为2(winsize+1)+1的立方体领域内搜索新的边界体素，并去除已填充的边界体素，而不需要重新搜索整个边界；对于填充边界模板Ψ_p后新产生的边界体素P′，需要设置其匹配中心因子MCF和最大距离MD；在更新边界模块Ψ_p的置信度和临时优先级值时，只需要对以边界体素P为中心、边长为4winsize+1大小的立方体邻域内的所有边界体素进行重新计算。

（13）重复步骤（9）～（12）直至填充完待修复区域Ω，即

修复后的岩心三维图像如图6-3所示。

如图5-1、图5-2、图5-3和图6-3所示，可以直观清楚的看到修复结果在各个方向上的结构和纹理都具有完整性和连续性，保持了原始三维图像中已知部分的结构和纹理特征。将图6-3和图6-1对比，并将图5-1、图5-2、图5-3与3-1、图3-2、图3-3对比，可以看出修复后的三维图像与预处理后的真实三维图像在结构和纹理上都极为相似，即说明本发明的岩心三维图像修复效果是较好的。

上述步骤中，步骤（6）～（13）为本发明所述修复方法的主要步骤。

本实施例中，我们将真实的岩心CT序列图像人为地去除中间一部分，模拟切割对三维图像造成的信息缺失，作为存在信息缺失的原始三维图像进行修复，并将其修复结果与真实岩心CT序列图像进行直观地对比观察，从而验证了本发明的可靠性。

上述实施例只是本发明的较佳实施例，并不是对本发明技术方案的限制，只要是不经过创造性劳动即可在上述实施例的基础上实现的技术方案，均应视为落入本发明专利的权利保护范围内。

Claims (2)

1.岩心三维图像修复方法，其特征在于：包括以下步骤：

（1）提取缺失信息的岩心三维图像即原始三维图像中待修复区域Ω的初始边界δΩ⁰，并设定模板大小；

（2）对原始三维图像的已知区域Φ进行区域划分；

（3）计算原始三维图像的边界上所有体素q的边界模板的优先级值，所述体素q指三维数据场中的采样点；

（4）搜索最大优先级值的边界模板Ψ_p，并保存边界体素P的坐标；

（5）判断边界体素P的区域类型，并设置边界模板Ψ_p的匹配块的全局搜索区域；

（6）根据边界体素P的匹配中心因子MCF和最大距离MD在相应的范围内搜索最优匹配块Ψ_q，并进行相应体素的灰度值填充；

（7）更新或设置以下信息：置信度、边界面δΩ、新产生的边界体素P′的匹配中心因子MCF和最大距离MD，以及新产生的边界模板Ψ_p′的优先级值；

（8）重复步骤（4）～（7）直至填充完待修复区域Ω，即

2.根据权利要求1所述的岩心三维图像修复方法，其特征在于：

所述步骤（2）中，将已知区域Φ划分为各相（R₁，R₂，…R_n）和混合区域R₀，有n相（D₁，D₂，…D_n）的岩心，划分的区域种类为n+1类，划分规则为：若

则q∈R_m；若

D(k)≠D(l)，则q∈R₀；其中，q为某一体素，Ψ_q为样本块，R₀为混合区域，D_m为某一相，R_m为与D_m相对应的区域；

所述步骤（3）中，原始三维图像中的每一个体素q都存在一个灰度值和置信度C(q)，当体素没有被填充时灰度值为空，体素q的置信度C(q)表示q的灰度值的可信度，初始化为： $\begin{array}{cc}C\left(q\right)=0& \∀q\∈\mathrm{\Ω}\end{array}$ 和 $\begin{array}{cc}C\left(q\right)=1& \∀q\∈\mathrm{\Φ}\end{array};$

每一个边界体素P的模板即边界模板Ψ_p都有一个优先级值P(p)，其定义为：

P(p)＝C(p)D(p)

上式中，C(p)为边界模板Ψ_p的置信度，用于衡量体素P周围模块范围内可靠信息的数量，D(p)为数据项，表示等照度面在边界面δΩ处的强度，C(p)和D(p)的计算公式为：

$C\left(p\right)=\frac{{\mathrm{\Σ}}_{q\∈{\mathrm{\Ψ}}_p\∩\stackrel{\‾}{\mathrm{\Ω}}}C\left(q\right)}{\left|{\mathrm{\Ψ}}_p\right|}$

$D\left(p\right)=\frac{|\▿I_p|\·\left|\sin \mathrm{\θ}\right|}{\mathrm{\α}}$

$\sin \mathrm{\θ}=\frac{\▿I_p\·n_p}{|\▿I_p|*\left|n_p\right|}$

其中，

表示边界模板Ψ_p中已被填充的体素，|Ψ_p|表示边界模板Ψ_p的体积，α为归一化因子，为边界体素P的梯度向量，该梯度包含方向和大小，大小等于等照度面强度

θ为边界体素P处的边界面δΩ与等照度面之间的夹角，也等于梯度向量

与边界面δΩ的法线n_p之间的夹角；

所述步骤（4）中，边界体素P是边界δΩ上的体素，边界δΩ为所有边界体素P的集合；

所述步骤（5）中，边界体素P的区域类型的判定方法为：根据模板范围内所有已知体素的相位判定该边界体素P属于哪一相或混合区域；根据边界体素P的区域类型设置边界模板Ψ_p的匹配块的全局搜索区域；

所述步骤（6）中，所述匹配中心因子MCF是边界体素P相应的匹配块局部搜索区域的中心体素；所述最大距离MD表示在边界体素P的匹配中心因子MCF不为空的情况下，局部搜索得到的最优匹配块Ψ_q与其边界模板Ψ_p的距离的最大范围；当局部搜索得到的匹配块与相应的边界模板Ψ_p的距离大于其最大距离MD时，该匹配块将不能作为其最优匹配块，需重新全局搜索其最优匹配块；某一边界体素P′的最大距离MD表示为：

MD_P′＝μ*d(Ψ_p,Ψ_q)

其中μ为比例参数，通常设定为1及其较为临近的值，边界体素P′为填充边界模板Ψ_p后新产生的边界体素，根据边界模板Ψ_p和其最优匹配块Ψ_q两个模板之间的距离，设置边界体素P′的最大距离MD；在设置新产生的边界体素P′的匹配中心因子MCF时，同步设置其最大距离MD；匹配准则为：

${\mathrm{\Ψ}}_q=\arg \underset{{\mathrm{\Ψ}}_{\hat{q}}\∈\mathrm{\Φ}}{\min }d({\mathrm{\Ψ}}_p,{\mathrm{\Ψ}}_{\hat{q}})$

其中两模板的距离

为边界模板Ψ_p和样本块对应位置已知体素的灰度值之差的平方和；

得到最优匹配块Ψ_q后，用最优匹配块Ψ_q中已知的体素的灰度值填充边界模板Ψ_p相应位置上待修复的体素，此时已知区域Φ和待修复区域Ω也都相应地发生了改变；

所述步骤（7）中，边界模板Ψ_p中新填充的体素的置信度设置为边界模板Ψ_p的置信度，即是：

$\begin{array}{cc}C\left(q\right)=C\left(p\right)& \∀q\∈{\mathrm{\Ψ}}_p\∩\mathrm{\Ω}\end{array}$

若设置的模板是边长为2winsize+1大小的立方体，那么在搜索新的边界体素时，只需要以刚填充的边界体素P为中心，在边长为2(winsize+1)+1的立方体领域内搜索新的边界体素，并去除已填充的边界体素，而不需要重新搜索整个边界；对于填充边界模板Ψ_p后新产生的边界体素，需要设置其匹配中心因子MCF和最大距离MD；在更新边界模块Ψ_p的置信度和临时优先级值时，只需要对以边界体素P为中心、边长为4winsize+1大小的立方体邻域内的所有边界体素进行重新计算。

Images