CN111784724B - 改进型马尔科夫链蒙特卡洛二维岩石切片重构方法及系统 - Google Patents

Abstract

本发明公开了一种改进型马尔科夫链蒙特卡洛二维岩石切片重构方法及系统，对于岩石切片电镜扫描图像，计算其孔隙度；确定要用到的邻域系统及各邻域系统不同组合对应的条件概率；在偶数行从右向左迭代求值，奇数行从左向右迭代求值，每行的第一个和最后一个像素点由它的上一行的同列像素点通过竖向二邻域求得，得到重构路径；根据岩石切片图像表现出的各向异性的特征，考虑到条件概率也存在各向异性，故沿着重构路径，扫描二值化图像各方向的条件概率；利用蒙特卡洛计算赋值，得到重构图像。本发明使条件概率扫描方向与重构路径相匹配，并选取随机数重构图像，把各方向的影响作用均衡化，表征各向异性的同时提高重构图像与原图的变差函数的拟合度。

Description

改进型马尔科夫链蒙特卡洛二维岩石切片重构方法及系统

技术领域

本发明属于计算机图像处理技术领域，涉及一种改进型马尔科夫链蒙特卡洛二维岩石切片重构方法及系统。

背景技术

图像重构因高便利和低成本广泛应用于医学、地质等各领域，其目的是对未知区域进行特征预判和结构分析。为了使重构图像的特性更贴合实际，选择合适的成像方法是影响重构效果的直接因素。

现有技术中，Kejian Wu的文章“An Efficient Markov Chain Model for theSimulation of Heterogeneous Soil Structure”对MCMC方法进行了系统的阐述，利用马尔科夫链的思想，选取合适的邻域系统从某一固定方向遍历扫描岩石电镜切片图像，对其分布规律进行统计分析得条件概率，再利用蒙特卡洛算法，求取随机数解决计算问题，通过变差函数的拟合来判断重构结果的优劣。但是，路径会在重构过程中不断变换方向，而且为了让结果更贴合原图，会为条件概率设置权重，直接改变了统计信息，并且变差函数的拟合效果无法定量判断。因此，在岩石图像的各向异性普遍存在的情况下，本发明提到的改进方法更合理。

发明内容

针对现有技术中的缺陷和不足，本发明提供了一种改进型马尔科夫链蒙特卡洛二维岩石切片重构方法及系统，目的是针对岩石图像的各向异性这一广泛存在的特征，以及目前现有的技术未考虑由这一特性导致的条件概率也存在各向异性的问题。

为达到上述目的，本发明采取如下的技术方案：

一种改进型马尔科夫链蒙特卡洛二维岩石切片重构方法，包括以下步骤：

步骤一，对原图进行二值化，计算孔隙度：

针对大小为m*n的原图，通过最大类间方差求得图像的阈值，根据阈值将原图二值化，统计不同像素值像素点总数，计算孔隙度如下：

二值化图像中像素值为0的代表骨架，像素值为1的代表孔隙；上式中，num₁是像素值为1的像素点的总个数，num₀是像素值为0的像素点的总个数；

步骤二：确定邻域系统和条件概率计算公式：

在步骤一得到的二值化图像基础上确定要用到的邻域系统，列出邻域系统所有可能出现的组合，并求每个组合条件下当前所求像素点分别为0或1的概率；条件概率计算公式为：

该邻域条件下所求像素点为0的概率：

该邻域条件下所求像素点为1的概率：

上式中，N是领域系统中的某个组合，x_r是该组合中的像素点，x_(i,j)是当前所求像素点，num_r→0是在该组合条件下x_(i,j)为0的次数，num_r→1是x_(i,j)为1的次数；

步骤三：确定重构路径及条件概率扫描方向：

设置一张与原图大小相同的空图像，空图像中，在偶数行从右向左迭代求值，奇数行从左向右迭代求值，每行的第一个和最后一个像素点由它的上一行的同列像素点通过竖向二邻域求得，得到重构路径；根据重构路径，扫描步骤一得到的二值化图像各方向的条件概率；

步骤四：利用蒙特卡洛计算赋值：

把二值化图像中第一个像素点赋值空图像的第一个像素点，之后对第一行从左向右使用横向二邻域进行计算赋值，之后第二行最后一个像素值通过竖向二邻域，将第一行最后一个像素点作为条件，计算像素值，然后从第二行第n-1个像素点开始，用左向四邻域进行匹配计算，第三行从左向右迭代，之后的换行重复之前操作，直到最后一个像素点被确定；

步骤五：孔隙度满足条件停止循环，得到重构图像；

将孔隙度作为循环结束的条件，当重构图像的孔隙度P_KXD与原图孔隙度p_kxd满足下式，停止循环，得到重构图像；

p_kxd-0.05≤P_KXD≤p_kxd+0.05                (4)

否则，重新设置空图像重复步骤四和步骤五直至满足公式(4)，最终得到重构图像。

本发明还包括如下技术特征：

具体的，对步骤五得到的重构图像进行评价，方法如下：分别计算二值化图像和重构图像的变差值；

变差值计算公式如下：

其中，a是步长，即两点间距，在这里选取为1到50，N(a)是距离为a的两点的个数，x是像素值，i表示像素点坐标；

计算二值化图像和重构图像的变差函数曲线之间的距离平均值，计算公式如下：

式中，R(a)是重构图像的变差值，r(a)为二值化图像的变差值；通过数据表征二值化图像和重构图像的变差曲线的拟合度，取重构图像与二值化图像在同一步长下变差值的差值，将不同步长下的差值累加取平均值，越小越拟合。

本发明还提供一种基于马尔科夫链蒙特卡洛的岩石切片重构系统，包括：

孔隙度计算模块，用于通过最大类间方差求得原图的阈值，根据阈值将原图二值化，二值化后的图片中像素值为0的代表骨架，像素值为1的代表孔隙，统计像素值分别为0和1的像素点的总数，计算得到孔隙度；

领域系统及条件概率确定模块，用于在二值化图像基础上确定邻域系统，列出邻域系统所有可能出现的组合，并求每个组合条件下像素值分别为0或1的概率；

重构路径及各方向条件概率确定模块，用于在空图像中，在偶数行是从右向左迭代求值，奇数行是从左向右迭代求值，每行的第一个和最后一个像素点由它的上一行的同列像素点通过竖向二邻域求得，确定重构路径；根据重构路径，扫描步骤一得到的二值化图像各方向的条件概率；

计算赋值模块，用于利用蒙特卡洛依次迭代计算空图像中的每一个像素值；

获取重构图像模块，用于将孔隙度作为循环结束的条件，得到重构图像；

结果拟合度评价模块，用于通过数据表征二值化图像和重构图像的变差曲线的拟合度，取重构图像与二值化图像在同一步长下变差值的差值，将不同步长下的差值累加取平均值，越小越拟合。

本发明与现有技术相比，有益的技术效果是：

本发明提出的MCMC重构方法，依据岩石图像的各向异性，考虑条件概率也存在各向异性的特征，从不同方向扫描得该方向的条件概率，在迭代求值过程中，使条件概率扫描方向与重构路径相匹配，并选取随机数重构图像，把各方向的影响作用均衡化，在表征各向异性的同时，提高重构图像与原图像的变差函数的拟合度。

附图说明

图1是岩石电镜扫描原图像A与二值化示意图B；

图2是横向和竖向二邻域示意图；(灰色为邻域系统像素点)

图3是左向和右向四邻域示意图；(灰色为邻域系统像素点)

图4是蛇形重构路径示意图；

图5是现有技术方法和本发明成像结果示意图；

图6是变差函数拟合曲线示意图；

图7是本发明流程图。

以下对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。

具体实施方式

本发明针对岩石切片电镜扫描图像表现的各向异性的特征，利用马尔科夫链思想和蒙特卡洛采样方法，考虑到了因岩石图像的特征信息导致在不同方向条件概率表现出的各向异性，并利用迭代计算的方式进行成像显示。

实施例1：

如图1至图7所示，本实施例提供一种改进型马尔科夫链蒙特卡洛二维岩石切片重构方法，图7为本发明的流程图，本方法对条件概率的扫描方式进行各向异性的考量，使条件概率扫描方向与重构路径相匹配，并选取随机数重构图像，把各方向的影响作用均衡化，表征各向异性的同时提高重构图像与原图的变差函数的拟合度；具体包括以下步骤：

步骤一，对原图进行二值化，计算孔隙度：

图1左A为原图，针对大小为m*n的原图，通过最大类间方差求得图像的阈值，根据阈值将原图二值化，图1右B是二值化岩石图像，统计不同像素值像素点总数，计算孔隙度如下：

二值化图像中像素值为0的代表骨架，像素值为1的代表孔隙；上式中，num₁是像素值为1的像素点的总个数，num₀是像素值为0的像素点的总个数；具体的，在本实施例中，图1右B二值化岩石图像计算其孔隙率为17.68％；

步骤二：确定邻域系统和条件概率计算公式：

在步骤一得到的二值化图像基础上确定要用到的邻域系统，列出邻域系统所有可能出现的组合，并求每个组合条件下当前所求像素点分别为0或1的概率；条件概率计算公式为：

该邻域条件下所求像素点为0的概率：

该邻域条件下所求像素点为1的概率：

上式中，N是领域系统中的某个组合，x_r是该组合中的像素点，x_(i,j)是当前所求像素点，num_r→0是在该组合条件下x_(i,j)为0的次数，num_r→1是x_(i,j)为1的次数；

本实施例中，图2分别为横向和竖向二邻域，它的条件概率计算公式为：

其中，N₂是二邻域的某个组合，x_r是该组合中的像素点，x_(i,j)是当前所求像素点，num_r→0是在该组合条件下x_(i,j)为0的次数，num_r→1是x_(i,j)为1的次数；左边表达式计算该邻域条件下所求像素点为0的概率，右边表达式计算该邻域条件下所求像素点为1的概率；

图3为左向和右向四邻域，它的条件概率计算公式为：

其中，N₄是二邻域的某个组合，x_r是该组合中的像素点，x_(i,j)是当前所求像素点，num_r→0是在该组合条件下x_(i,j)为0的次数，num_r→1是x_(i,j)为1的次数；左边表达式计算该邻域条件下所求像素点为0的概率，右边表达式计算该邻域条件下所求像素点为1的概率；相关解释与二邻域类似；

步骤三：确定重构路径及条件概率扫描方向：

图4为蛇形路径的原理图，设置一张与原图大小相同的空图像，空图像中，在偶数行是从右向左迭代求值，奇数行是从左向右迭代求值，每行的第一个和最后一个像素点由它的上一行的同列像素点通过竖向二邻域求得，得到重构路径；根据重构路径，扫描步骤一得到的二值化图像各方向的条件概率；

步骤四：利用蒙特卡洛计算赋值：

把二值化图像中第一个像素点赋值空图像的第一个像素点，之后对第一行从左向右使用横向二邻域进行计算赋值，之后第二行最后一个像素值通过竖向二邻域，将第一行最后一个像素点作为条件，计算像素值，然后从第二行第n-1个像素点开始，用左向四邻域进行匹配计算，第三行从左向右迭代，之后的换行重复之前操作，直到最后一个像素点被确定；

步骤五：孔隙度满足条件停止循环，得到重构图像；

将孔隙度作为循环结束的条件，当重构图像的孔隙度P_KXD与原图孔隙度p_kxd满足下式，停止循环，否则，重新设置空图像重复步骤四和步骤五；

p_kxd-0.05≤P_KXD≤p_kxd+0.05              (4)

具体的，图5A是现有技术中使用对二值化图像扫描一次得到的条件概率的重构图像；图5B是本实施例中使用与重构方向相同的条件概率的重构图像；这两张图像的孔隙度分别为14.21％和19.36％，都达到了循环结束的条件。

步骤六：分别计算二值化图像和重构图像的变差值；

变差值计算公式如下：

其中，a是步长，即两点间距，在这里选取为1到50，N(a)是距离为a的两点的个数，x是像素值，i表示像素点坐标；

计算二值化图像和重构图像的变差函数曲线之间的距离平均值，计算公式如下：

式中，R(a)是重构图像的变差值，r(a)为二值化图像的变差值；通过数据表征二值化图像和重构图像的变差曲线的拟合度，取重构图像与二值化图像在同一步长下变差值的差值，将不同步长下的差值累加取平均值，越小越拟合。

图6是三张图像的变差函数，标点是原图变差值，红线是现有技术图5A，蓝线是本实施例图5B，可以看到蓝线更加拟合标点，它们与标点之间的平均距离分别为0.0208(图5A)和0.0073(图5B)，在重构过程中使用与重构路径同方向的条件概率，即使用本发明方法的成像效果更好。

实施例2：

本实施例提供一种基于马尔科夫链蒙特卡洛的岩石切片重构系统，包括：

孔隙度计算模块，用于通过最大类间方差求得原图的阈值，根据阈值将原图二值化，二值化后的图片中像素值为0的代表骨架，像素值为1的代表孔隙，统计像素值分别为0和1的像素点的总数，计算得到孔隙度；

领域系统及条件概率确定模块，用于在二值化图像基础上确定邻域系统，列出邻域系统所有可能出现的组合，并求每个组合条件下像素值分别为0或1的概率；

重构路径及各方向条件概率确定模块，用于在空图像中，在偶数行是从右向左迭代求值，奇数行是从左向右迭代求值，每行的第一个和最后一个像素点由它的上一行的同列像素点通过竖向二邻域求得，确定重构路径；根据重构路径，扫描步骤一得到的二值化图像各方向的条件概率；

计算赋值模块，用于利用蒙特卡洛依次迭代计算确定空图像的每个点的像素值；

获取重构图像模块，用于将孔隙度作为循环结束的条件，得到重构图像；

结果拟合度评价模块，用于通过数据表征二值化图像和重构图像的变差曲线的拟合度，取重构图像与二值化图像在同一步长下变差值的差值，将不同步长下的差值累加取平均值，越小越拟合。

Claims (3)

1.一种改进型马尔科夫链蒙特卡洛二维岩石切片重构方法，其特征在于，具体包括以下步骤：

步骤一，对原图进行二值化，计算孔隙度：

针对大小为m*n的原图，通过最大类间方差求得图像的阈值，根据阈值将原图二值化，统计不同像素值像素点总数，计算孔隙度如下：

二值化图像中像素值为0的代表骨架，像素值为1的代表孔隙；上式中，num₁是像素值为1的像素点的总个数，num₀是像素值为0的像素点的总个数；

步骤二：确定邻域系统和条件概率计算公式：

在步骤一得到的二值化图像基础上确定要用到的邻域系统，列出邻域系统所有可能出现的组合，并求每个组合条件下当前所求像素点分别为0或1的概率；条件概率计算公式为：

该邻域条件下所求像素点为0的概率：

该邻域条件下所求像素点为1的概率：

上式中，N是领域系统中的某个组合，x_r是该组合中的像素点，x_(i,j)是当前所求像素点，num_r→0是在该组合条件下x_(i,j)为0的次数，num_r→1是x_(i,j)为1的次数；

步骤三：确定重构路径及条件概率扫描方向：

设置一张与原图大小相同的空图像，空图像中，在偶数行从右向左迭代求值，奇数行从左向右迭代求值，每行的第一个和最后一个像素点由它的上一行的同列像素点通过竖向二邻域求得，得到重构路径；根据重构路径，扫描步骤一得到的二值化图像各方向的条件概率；

步骤四：利用蒙特卡洛计算赋值：

把二值化图像中第一个像素点赋值空图像的第一个像素点，之后对第一行从左向右使用横向二邻域进行计算赋值，之后第二行最后一个像素值通过竖向二邻域，将第一行最后一个像素点作为条件，计算像素值，然后从第二行第n-1个像素点开始，用左向四邻域进行匹配计算，第三行从左向右迭代，之后的换行重复之前操作，直到最后一个像素点被确定；

步骤五：孔隙度满足条件停止循环，得到重构图像；

将孔隙度作为循环结束的条件，当重构图像的孔隙度P_KXD与原图孔隙度p_kxd满足下式，停止循环，得到重构图像；

p_kxd-0.05≤P_KXD≤p_kxd+0.05                (4)

否则，重新设置空图像重复步骤四和步骤五直至满足公式(4)，最终得到重构图像。

2.如权利要求1所述的改进型马尔科夫链蒙特卡洛二维岩石切片重构方法，其特征在于，对所述步骤五得到的重构图像进行评价，方法如下：分别计算二值化图像和重构图像的变差值；

变差值计算公式如下：

其中，a是步长，即两点间距，在这里选取为1到50，N(a)是距离为a的两点的个数，x是像素值，i表示像素点坐标；

计算二值化图像和重构图像的变差函数曲线之间的距离平均值，计算公式如下：

式中，R(a)是重构图像的变差值，r(a)为二值化图像的变差值；通过数据表征二值化图像和重构图像的变差曲线的拟合度，取重构图像与二值化图像在同一步长下变差值的差值，将不同步长下的差值累加取平均值，越小越拟合。

3.一种改进型马尔科夫链蒙特卡洛的岩石切片重构系统，其特征在于，包括：

孔隙度计算模块，用于通过最大类间方差求得原图的阈值，根据阈值将原图二值化，二值化后的图片中像素值为0的代表骨架，像素值为1的代表孔隙，统计像素值分别为0和1的像素点的总数，计算得到孔隙度；

领域系统及条件概率确定模块，用于在二值化图像基础上确定邻域系统，列出邻域系统所有可能出现的组合，并求每个组合条件下像素值分别为0或1的概率；

重构路径及各方向条件概率确定模块，用于在空图像中，在偶数行是从右向左迭代求值，奇数行是从左向右迭代求值，每行的第一个和最后一个像素点由它的上一行的同列像素点通过竖向二邻域求得，确定重构路径；根据重构路径，扫描步骤一得到的二值化图像各方向的条件概率；

计算赋值模块，用于利用蒙特卡洛依次迭代计算空图像中的每一个像素值；

获取重构图像模块，用于将孔隙度作为循环结束的条件，得到重构图像；

结果拟合度评价模块，用于通过数据表征二值化图像和重构图像的变差曲线的拟合度，取重构图像与二值化图像在同一步长下变差值的差值，将不同步长下的差值累加取平均值，越小越拟合。

Images