CN110223376B - 一种基于单幅堆积颗粒材料图像的三维颗粒重建方法 - Google Patents

Abstract

本发明公开了一种基于单幅堆积颗粒材料图像的三维颗粒重建方法，依次通过获取堆积颗粒材料的颗粒粒径列表；基于单幅堆积颗粒材料的彩色图像提取颗粒的特征轮廓；利用已知颗粒特征轮廓，生成具有不同尺寸的颗粒表面散点模型；以及基于颗粒的表面散点模型，按照颗粒粒径列表的要求对三维堆积状态的颗粒进行重建；该基于单幅堆积颗粒材料图像的三维颗粒重建方法通过单幅图片中颗粒的特征轮廓实现重建数字颗粒模型，并保证生成的颗粒在整体上具有真实颗粒的特征，且在操作上具有极大的便捷性，具备批量生成三维颗粒模型的优点。

Description

一种基于单幅堆积颗粒材料图像的三维颗粒重建方法

技术领域

本发明涉及模型构建技术领域，特别涉及一种基于单幅堆积颗粒材料图像的三维颗粒重建方法。

背景技术

胶结颗粒材料中骨料的数字重建是正确建立胶结颗粒材料细观数字模型的前提。现行的骨料重建方法有很多，其中最具代表性的包括：

1)使用球体、椭球体和凸多面体等通过数学和几何学规则构建出的骨料颗粒模型；该方法的优点是构造方法简单、快捷，构造出的颗粒模型的参数易于控制；缺点是只能够构造具有一定规则形状的颗粒，与实际的颗粒形状有较大的差距；

2)通过微焦点CT设备，对材料进行全方位的断层扫描获取相关的CT图像序列，再利用数字图像处理技术处理CT图像序列，最后得到三维的骨料模型。此方法制备的骨料模型在外形上比较接近真实的骨料颗粒，但是CT设备昂贵且操作繁琐，并且在数字图像处理的过程中会丢失部分数据。

发明内容

本发明的目的是提供一种基于三维堆积颗粒的单幅二维图像实现三维颗粒准确且便捷地对颗粒进行重建的基于单幅堆积颗粒材料图像的三维颗粒重建方法。

为此，本发明技术方案如下：

一种基于单幅堆积颗粒材料图像的三维颗粒重建方法，步骤如下：

S1、获取堆积颗粒材料的颗粒粒径列表；

S2、基于单幅堆积颗粒材料的彩色图像提取其上颗粒的特征轮廓；

S3、利用步骤S2得到的颗粒特征轮廓，生成具有不同尺寸的颗粒表面散点模型；其中，每个颗粒表面散点模型的生成方法为：

1)根据待重建颗粒的具体轮廓要求，在颗粒轮廓数据库中随机选取两个符合要求的颗粒轮廓，并分别指定为导向轮廓和移动轮廓；

2)将导向轮廓水平放置于三维虚拟空间中，并将轮廓A旋转至其最长轴与x轴平行处；将最长轴等距划分为n个节点P_i，获取导向轮廓上与节点P_i的在x轴和z轴上坐标值相同的两个点：点M_i和点N_i之间的距离H；

3)将移动轮廓竖直放置于三维虚拟空间中，并以其重心坐标为节点，获取过其重心且与y轴平行的轴线与移动轮廓相交的两点：点O_j和点P_j之间的距离D，以计算缩放系数

并利用缩放矩阵按照比例缩放移动轮廓；接着，计算点M_i或点N_i与点O_j或点P_j之间的距离差值，通过以该距离差值为基础的平移矩阵将移动轮廓移动到导向轮廓的每个节点上，即形成颗粒表面散点模型；

S4、基于步骤三得到的各种尺寸的颗粒的表面散点模型，按照符合步骤一生成的三维堆积状态的颗粒的颗粒粒径列表的要求对三维堆积状态的颗粒进行重建。

2、根据权利要求1所述的基于单幅堆积颗粒材料图像的三维颗粒重建方法，其特征在于，步骤S1的具体实施步骤为：

S101、对堆积颗粒材料或部分堆积颗粒材料进行颗粒分析试验，得到原始颗粒级配曲线；

S102、设定待生成颗粒的总个数N，以及其极限粒径mind和maxd，将原始级配曲线分为若干个粒组，

S103、根据颗粒分析试验后得到的每个粒组所占体积百分比，计算出每个粒组中待生成颗粒的总个数；

S104、以每个粒组中颗粒直径的上、下界的粒径为界限，根据步骤S103计算出的该粒组中待生成颗粒的个数，线性插值出该粒组内粒径从小到大且满足步骤S101得到的颗粒级配曲线粒径变化趋势的每个颗粒的粒径；

S105、保存所有粒组内所有颗粒的粒径，生成颗粒粒径列表。

进一步地，步骤S2的具体实施步骤为：

S201、将单幅堆积颗粒材料的彩色图像转变为灰度图像；

S202、采用Sobel算子分别对灰度图像从水平和竖直两个方向进行滤波，并求取其梯度幅值矩阵；

S203、在经过步骤S201处理得到的灰度图像中分别对背景对象和目标对象进行标记；其中，目标对象指图像中的颗粒所在区域，背景对象指图像中颗粒所在区域与非颗粒所在区域之间的分界线；

S204、对经过步骤S202处理得到的梯度幅值矩阵进行修正，即将目标对象和背景对象对应的点在梯度幅值矩阵中修正为一个最小值Inf，完成对梯度幅值矩阵的修正；

S205、从经过步骤S204得到的图像中提取颗粒的轮廓，并利用MPP算法对颗粒的轮廓进行优化，并计算每一个颗粒轮廓的圆度、椭圆度和综合形状指数；

S206、根据步骤S205计算得到的每个轮廓的圆度、椭圆度和综合形状指数，筛选出符合预期的颗粒轮廓，并生成颗粒轮廓数据库。

进一步地，在步骤S203中，对目标对象进行标记的具体实施步骤为：

1)以经过步骤S201处理得到的灰度图像作为掩模图像，以对掩模图像进行腐蚀运算后的图像作为样本图像，进行第一次重建运算；

2)对经过第一次重建运算后得到的图像进行求补运算，并将经过求补运算后得到的图像作为掩模图像；对经过第一次重建运算后得到的图像依次进行膨胀和求补运算，并将结果图像作为样本图像，进行第二次重建运算；然后对经过第二次重建运算后得到的图像进行求补运算；

其中，重建运算为利用样本图像f重建掩模图像g为Rg(f)，具体为：获取样本图像f；创建结构元素：B＝ones(3)；其中，ones是MATLAB中的系统函数，(3)表示选取3×3的矩阵；重复f_k+1＝(f_k+B)∩g，直至f_k+1＝f_k；其中，样本图像f为掩模图像g的一个子集；求补运算过程为对图像的灰度值进行取反运算，如果像素灰度为0，求补运算就是255，如果像素灰度为10，求补运算为245，也就是原始像素灰度值与其补集灰度值之和为255。

3)利用区域连通性经过步骤1)和步骤2)得到的图像进行处理：确定图像内由灰度值为0的像素包围出的多个局部区域，并将每个局部区域内出现的极大值的像素标记为1，其余像素标记为0，得到完成对目标对象进行标记的二值图像；

4)对经过上述步骤3)处理得到图像进行先开运算再闭运算的处理，完成对目标图像的标记。

进一步地，在步骤S203中，对背景对象进行标记的具体实施步骤为：

1)通过自动阈值分割算法对经过步骤S201处理得到的灰度图像进行处理，得到该灰度图像的最佳分割阈值；

2)根据经过步骤1)得到的最佳分割阈值T^*对灰度图像进行二值化处理：将灰度图像中灰度值为大于T^*的像素位置处的数值定义为1，即标识为目标对象；将灰度图像中灰度值为小于T^*的像素位置处的数值定义为0，即标识为非目标对象区域；

3)对经过步骤2)处理得到的二值图像进行距离变换，得到距离矩阵；

4)将距离矩阵代入watershed函数中进行分水岭变换后，所得分水岭脊线即为背景图像。

进一步地，在步骤S3中，每个颗粒表面散点模型生成前，对导向轮廓和移动轮廓进行如下处理：将导向轮廓和移动轮廓从笛卡尔直角坐标系下的坐标转换为极坐标系下的坐标，然后对极坐标列表依次进行等距线性插值和采用高斯加权移动平均法去除波动点的处理后再转换为笛卡尔直角坐标系下的坐标。

进一步地，在步骤S3的颗粒表面散点模型生成过程中，导向轮廓上与节点P_i的在x轴上坐标值相同的两个点之间的距离H的具体计算方法为：将轮廓A水平放置于三维虚拟空间Ω中，并去除其中重复的坐标点；计算轮廓A最长轴长度L以及其与x轴的夹角ω；根据夹角ω将轮廓A旋转至其最长轴与x轴平行处，并将最长轴等距划分为n个节点P_i，每个节点给出其对应坐标；将轮廓A以最长轴的端点为界限分为A₁和A₂两个部分，然后根据节点P_i的横坐标值，找出在A₁中距离其最近点A、B，在A₂中距离其最近点A’、B’，由A、B和A’、B’线性插值出导向轮廓上与节点P_i在x轴与z轴上坐标值相同的点M_i和N_i的坐标，进而根据点M_i和点N_i的y坐标值，计算点M_i和点N_i之间的距离记为H。

进一步地，在步骤S3的颗粒表面散点模型生成过程中，将移动轮廓竖直放置于三维虚拟空间Ω中，并以其重心G坐标为节点；经过重心G作一条平行于y轴的轴线，该轴线与移动轮廓有两个交点：点O_j和点P_j；找到y坐标值较大的交点P_j，并以该交点找到在轮廓B中距离其最近的两个点：点C和点D，由点C和点D插值出位于移动轮廓上且与重心G具有相同x坐标值和z坐标值的点在y轴上坐标值，同理，采用相同方法获得点O_j的坐标，根据点O_j和点P_j在y轴上的坐标值的差值，即得到移动轮廓上过重心G作一条平行于y轴的轴线的长度D；计算缩放系数

利用缩放矩阵按照比例缩放移动轮廓；接着，计算点M_i或点N_i与点O_j或点P_j之间的距离差值，通过以该距离差值为基础的平移矩阵将移动轮廓移动到导向轮廓的每个节点上，颗粒表面散点模型形成。

该基于单幅堆积颗粒材料图像的三维颗粒重建方法依次通过根据堆积颗粒材料的级配曲线生成待制备颗粒的粒径列表，再由一张堆积颗粒材料的二维图片提取其中颗粒的特征轮廓，最后随机选取两个特征轮廓分别作为移动轮廓和导向轮廓按照生成的颗粒粒径要求重建三维颗粒三个步骤实现。

与现有技术相比，该基于单幅堆积颗粒材料图像的三维颗粒重建方法通过单幅图片中颗粒的特征轮廓实现重建数字颗粒模型，并保证生成的颗粒在整体上具有真实颗粒的特征，且在操作上具有极大的便捷性，具备批量生成三维颗粒模型的优点。

附图说明

图1为本发明的基于单幅堆积颗粒材料图像的三维颗粒重建方法的流程图；

图2为本发明的实施例中经过步骤S1的颗粒分析试验得到的原始级配曲线和插值后得到级配曲线的图像；

图3为本发明的基于单幅堆积颗粒材料图像的三维颗粒重建方法中提取图像上黄豆颗粒的特征轮廓的流程图；

图4为本发明的实施例中经过步骤S201处理得到的原始灰度图像；

图5为本发明的实施例中经过步骤S201处理得到的梯度幅值矩阵的对应图像；

图6为本发明的实施例中将经过步骤S2031～S2033得到的二值图像与原始图像叠加后得到的图像；

图7为本发明的实施例中经过步骤S2034～S2037得到的标记有背景图像的二值图像；

图8为本发明的实施例中经过步骤S204将梯度幅值矩阵经分水岭变换后得到的图像；

图9为本发明的实施例中在步骤S204得到的图像中粗提取出某一黄豆颗粒轮廓的图像；

图10为本发明的实施例中将图9经过MPP算法优化后得到的黄豆颗粒轮廓的图像；

图11(a)为本发明的实施例中的步骤S302中轮廓B自笛卡尔直角坐标系下的坐标转换为在极坐标系下的坐标的示意图；

图11(b)为本发明的实施例中的步骤S302中轮廓B的极坐标经过插值后得到的极坐标的示意图；

图11(c)为本发明的实施例中的步骤S302中轮廓B经过插值后得到的极坐标经过去除波动点处理后得到极坐标的示意图；

图12为本发明的实施例中的步骤S302中将轮廓B在极坐标系下的坐标转换为在笛卡尔直角坐标系下的坐标的示意图；

图13为本发明的实施例中的步骤S3031中将轮廓A水平放置于三维虚拟空间Ω中的示意图；

图14为本发明的实施例中的步骤S3032中将轮廓B竖直放置于三维虚拟空间Ω中的示意图；

图15为本发明的实施例中经过步骤S3032得到的一颗黄豆颗粒的表面散点模型的示意图；

图16为本发明的实施例中将黄豆颗粒的表面散点模型离散化为若干个三角形网片结构的示意图。

具体实施方式

下面结合附图及具体实施例对本发明做进一步的说明，但下述实施例绝非对本发明有任何限制。

以下以三维堆积状态的黄豆颗粒为例，对该基于单幅堆积颗粒材料图像的三维颗粒重建方法的具体过程进行详细说明。

如图1所示，该基于单幅堆积颗粒材料图像的三维颗粒重建方法依次通过生成颗粒粒径列表、基于单幅堆积颗粒材料图像提取颗粒特征轮廓、重建三维骨料颗粒这三个步骤实现。

步骤一、生成三维堆积状态的黄豆颗粒的颗粒粒径列表；

如图1所示，该步骤的具体实施方式如下：

S101、通过颗粒分析试验，得到如图2中的曲线a所示的三维堆积状态下的黄豆颗粒的颗粒级配曲线，即原始颗粒级配曲线；

S102、设定待生成颗粒的总个数N，以及其极限粒径mind和maxd，将原始级配曲线分为若干个粒组，

如表1所示的将原始级配曲线分为九个粒组为例进行说明；

表1：

由于颗粒分析试验中，采用八种尺寸筛对黄豆颗粒进行过筛，因此将全部黄豆颗粒分为九个粒组，具体分别为[mind,5],[5,10],[10,20],[20,40],[40,60],[60,100],[100,200],[200,300],[300,maxd]；

S103、根据颗粒分析试验后得到的每个粒组所占体积百分比，计算出每个粒组中待生成颗粒的总个数；

S104、以每个粒组中颗粒直径的上、下界的粒径为界限，根据步骤S103计算出的该粒组中待生成颗粒的个数，线性插值出该粒组内每个颗粒的粒径，即每个粒组内每个颗粒的粒径均不同，且从小到大呈现与步骤S101得到的颗粒级配曲线粒径变化趋势一致，即如图2中的曲线b所示的插值后的级配曲线；

S105、保存所有粒组内所有颗粒的粒径，生成颗粒粒径列表，用于后续建立推及材料模型用。

步骤二、如图3所示，以单幅三维堆积状态的黄豆颗粒的彩色图像为基础，利用matlab软件提取图像上黄豆颗粒的特征轮廓；具体实施步骤如下：

S201、获取一幅如图4所示的三维堆积状态的黄豆颗粒的彩色图像，并将盖彩色图像转换为如图4所示的灰度图像；

彩色图像转换为灰度图像具体转换过程如下：

由于在RGB色彩空间中，一张彩色图像上的每个像素均可以分别提取出三个颜色分量，分别是R分量、G分量和B分量，它们的组合构成了我们平时看到图片的颜色；而灰度图像是通过图片中各像素亮度值的不同，展示图片的具体信息；因此，首先采用公式(1)的加权平均法，将彩色图像上每个像素的R分量、G分量和B分量带入公式(1)中，得到灰度值，进而实现将彩色图像转化为灰度图像；

Y＝0.2989*R+0.5870*G+0.1140*B 式(1)，

其中，Y是像素亮度值，R、G、B分别为彩色图像中每个像素的红色分量、绿色分量和蓝色分量。

S202、为了有效提取灰度图像中颗粒的轮廓信息，采用Sobel算子公式(2)分别对灰度图像从水平和竖直两个方向进行滤波，然后根据公式(3)和公式(4)求取其梯度幅值矩阵；

其中，Sobel算子公式(2)包括公式(2-1)和公式(2-2)，具体为：

其中，S_x为X方向的Sobel算子模板，S_y为Y方向的Sobel算子模板；

其中，A为表示原始图像中局部矩阵，G_x为经X方向的Sobel算子检测后的图片，G_y为经Y方向的Sobel算子检测后的图片；

其中，G为经Sobel算子检测后图片的梯度幅值矩阵；

经过步骤S202后，将灰度图像处理为如图5所示的图像。

S203、在经过步骤S201处理得到的灰度图像中标记背景对象和目标对象；其中，目标对象指图像中的黄豆区域，背景对象指图像中黄豆区域与非黄豆区域之间的分界线；

目标对象的具体标记方法为：

S2031、以经过步骤S201处理得到的灰度图像作为掩模图像，以对掩模图像进行腐蚀运算后的图像作为样本图像，，将原灰度图像变暗一些，使原灰度图像内一些较亮的细节被消除，从而可以避免灰度值较大的噪点对于轮廓提取的影响；然后对经过腐蚀运算后得到的图像进行第一次重建运算；

对经过第一次重建运算后得到的图像进行求补运算，并将经过求补运算后得到的图像作为掩模图像；对经过第一次重建运算后得到的图像依次进行膨胀和求补运算，并将结果图像作为样本图像，进行第二次重建运算；然后对经过第二次重建运算后得到的图像进行求补运算；

其中，重建运算指利用样本图像f重建掩模图像g为Rg(f)，具体为：

1)获取样本图像f；

2)利用MATLAB的系统函数创建结构元素：B＝ones(3)； 式(5)；

3)重复f_k+1＝(f_k+B)∩g，直至f_k+1＝f_k；其中，样本图像f为掩模图像g的一个子集；

求补运算过程为对图像的灰度值进行取反运算；即在本实施例中，如果像素的灰度值为0，其补集就是255；如果像素灰度值为10，其补集为245；也就是说，根据原始像素灰度值，以及原始像素灰度值与其补集灰度值之和为255这一条件，将像素的灰度值对应修改为其补集；

S2032、利用区域连通性经过步骤1)和步骤2)得到的图像进行处理：确定图像内由灰度值为0的像素包围出的多个局部区域，并将每个局部区域内出现的极大值的像素标记为1，其余像素标记为0，得到完成对目标对象进行标记的二值图像；

S2033、对经过上述步骤S2032处理后的图像依次通过公式(6)和公式(7)进行先开运算再闭运算的处理；

其中，F为步骤S2032处理得到的灰度图像的图像函数，S为对黄豆区域进行标记前预设的结构元素；通过该步骤不仅能够消除图片中过于明亮且比结构元素小的细节对于图片分析的影响，同时也能够保持图片灰度值不发生全局性的变化；如图6所示为经过该步骤后得到的图像与原始图像重叠后的图像，由图7可以看出，经过上述步骤的对目标对象进行标记，实现了对图像中黄豆区域的准确标记。

背景对象的具体标记方法为：

S2034、对经过步骤S201处理得到的灰度图像按照公式(8)～公式(15)进行自动阈值分割算法，得到该灰度图像的最佳分割阈值；具体地，

设灰度图像中灰度最大值为L，图像总像素个数设为N，则图像中灰度值为i的点出现的概率为：

P_i＝n_i/N 式(8)，

其中，n_i表示灰度值为i的像素的个数，

图像中灰度值的平均值为：

C₁出现的概率为：

C₂出现的概率为：

C₁类某像素出现的概率为：

C₂类某像素出现的概率为：

进而以求得类间方差

为：

因此，最佳阈值T^*则满足：

即得到最佳分割阈值T^*；

S2035、根据经过步骤S2034得到的最佳分割阈值T^*对灰度图像进行二值化处理：将灰度图像中灰度值为大于T^*的像素位置处的数值定义为1，即标识为目标对象；将灰度图像中灰度值为小于T^*的像素位置处的数值定义为0，即标识为非目标对象区域；

S2036、对经过步骤S2035处理得到的二值图像进行距离变换，即按照公式(16)求出每一个像素距离其最近的非零像素的欧氏距离，并将距离值替换为原二值图像内该像素内的数值，得到距离矩阵；

其中，(x_i,y_i)表示目标像素坐标，(x_j,y_j)表示距离目标像素最近的非0像素的坐标。

S2037、再对上述步骤S3得到的距离矩阵代入MATLAB中的系统函数：watershed函数中进行分水岭变换后，得到分水岭脊线，即背景图像，具体如图7所示的图像中的白色连线构成的背景对象；

S204、以经过步骤S2031～步骤S2033处理得到针对黄豆颗粒区域作为目标对象的处理后图像，以及经过步骤S2034～步骤S2037处理得到针对非黄豆颗粒区域作为背景对象的处理后图像为依据，对经过步骤S202处理得到的梯度幅值矩阵进行修正，即将目标对象和分水岭脊线对应的点在梯度幅值矩阵中修正为一个最小值Inf，如图8所示为梯度幅值矩阵经过上述分水岭变换步骤即得到的图像；

S205、在经过步骤S204得到的图像中提取黄豆颗粒轮廓，接着利用MPP算法将如图9所示的黄豆颗粒粗提取轮廓优化为如图10所示的黄豆颗粒轮廓，并计算每个黄豆颗粒轮廓的圆度、椭圆度和综合形状指数；具体地，

利用MPP算法提取经过步骤S204处理得到的图像中的每一个黄豆颗粒的轮廓，并根据每个黄豆颗粒轮廓计算出对应黄豆颗粒的圆度、椭圆度和综合形状指数，具体计算方式如下：

(1)黄豆颗粒圆度按照公式(17)计算得到：

其中，S为颗粒轮廓的包围面积，L为轮廓的周长；

(2)黄豆颗粒椭圆度按照公式(18)计算得到：

其中，L_a为颗粒轮廓的长度，L_b为颗粒轮廓的宽度；

(3)黄豆颗粒综合形状指数按照公式(19)计算得到：

其中，L_a为颗粒轮廓的长度，L_b为颗粒轮廓的宽度。

S206、根据步骤S205计算得到的每个轮廓的圆度、椭圆度和综合形状指数，筛选出符合预期的颗粒轮廓，并生成颗粒轮廓数据库便于后续调用。

步骤三、生成颗粒表面散点模型；

S301、根据待重建的黄豆颗粒的具体轮廓要求，在颗粒轮廓数据库中随机选取两个符合要求的黄豆颗粒轮廓，并分别指定为轮廓A和轮廓B；其中，轮廓A作为导向轮廓，在颗粒重建时负责引导移动轮廓移动且制约着颗粒在xoy面上的投影的形状；轮廓B作为移动轮廓，在颗粒重建时负责构建颗粒的表面；

S302、将轮廓A和轮廓B在笛卡尔直角坐标系下的坐标转换为极坐标系下的坐标；然后对极坐标列表进行等距线性插值，获得插值后的坐标列表；接着对该坐标列表采用高斯加权移动平均法获得去除波动点以后的坐标列表；最后再将轮廓A和轮廓B从极坐标系下的坐标转换为笛卡尔直角坐标；

该过程在对直角坐标系与极坐标系进行转化的具体方法为：先将轮廓平移至直角坐标原点处，以坐标原点为极坐标的极点；然后以坐标与极点之间的距离作为极径，就可以将直角坐标系转换为极坐标系；反之，将极坐标列表以极径为y轴，角度为x轴展示在笛卡尔直角坐标系中，即实现极坐标转化为直角坐标；

以轮廓B为例，如图11(a)所示为轮廓B自笛卡尔直角坐标系下的坐标转换为在极坐标系下的坐标的示意图；如图11(b)所示为轮廓B以极坐标系下的坐标为基础，经过插值后得到的极坐标的示意图；如图11(c)所示为轮廓B的经过插值后得到的极坐标经过去除波动点处理后得到极坐标的示意图；如图12所示为将轮廓B在极坐标系下的坐标转换为在笛卡尔直角坐标系下的坐标的示意图；其中，图12中的多个散点为轮廓B初始状态下在笛卡尔直角坐标系下的坐标散点，图12中的平滑曲线为轮廓B经过转换为极坐标后进行的一系列处理后再次转换为笛卡尔直角坐标系下坐标后绘制得到的曲线；

S303、以轮廓A为导向轮廓，轮廓B为移动轮廓，构建颗粒表面散点模型；具体实施步骤如下：

S3031、如图13所示，将轮廓A水平放置于三维虚拟空间Ω中，并去除其中重复的坐标点；计算轮廓A最长轴长度L以及其与x轴的夹角ω；根据夹角ω将轮廓A旋转至其最长轴与x轴平行处，并将最长轴等距划分为n个节点P_i，每个节点给出其对应坐标；将轮廓A以最长轴的端点为界限分为A₁和A₂两个部分，然后根据节点P_i的横坐标值，找出在A₁中距离其最近点A、B，在A₂中距离其最近点A’、B’，由A、B和A’、B’线性插值出导向轮廓上与节点P_i在x轴与z轴上坐标值相同的点M_i和N_i的坐标，进而根据点M_i和点N_i的y坐标值，计算点M_i和点N_i之间的距离记为H。

S3032、如图14所示，将轮廓B竖直放置于三维虚拟空间Ω中，并以其重心G坐标为节点；经过重心G作一条平行于y轴的轴线，该轴线与轮廓B有两个交点：点O_j和点P_j；找到y坐标值较大的交点P_j，并以该交点找到在轮廓B中距离其最近的两个点：点C和点D，由点C和点D插值出位于轮廓B上且与重心G具有相同x坐标值和z坐标值的点在y轴上坐标值，同理，采用相同方法获得点O_j的坐标，根据点O_j和点P_j在y轴上的坐标值的差值，即得到轮廓B上过重心G作一条平行于y轴的轴线的长度D；计算缩放系数

利用缩放矩阵按照比例缩放轮廓B；接着，计算点M_i或点N_i与点O_j或点P_j之间的距离差值，通过以该距离差值为基础的平移矩阵将移动轮廓移动到导向轮廓的每个节点上，这就形成了如图15所示的颗粒表面散点模型；

S3033、如图16所示，将经过步骤S3032处理得到的黄豆颗粒的颗粒表面散点模型离散化为若干个三角形网片，并以二进制存储和标准ASII码对点与面进行存储，形成STL文件；进而通过后续读取相关STL文件即可再现或存储三维模型。其中，如图16所示为经过步骤3033处理后得到的四种黄豆颗粒的三维模型。

步骤四、基于步骤三得到的各种尺寸的黄豆颗粒的表面散点模型，按照符合步骤一生成的三维堆积状态的黄豆颗粒的颗粒粒径列表的要求对三维堆积状态的黄豆颗粒进行重建。

利用三维扫描仪对三维堆积状态下的黄豆颗粒进行扫描，并将通过其得到的STL格式文件与采用本实施例的重建方法获得的黄豆颗粒的STL格式文件进行比较，两种方法得到的STL格式文件反映的黄豆颗粒的圆度、椭圆度和综合形状指数三种颗粒信息具有统计上的一致性，因此，本申请公开的基于单幅堆积颗粒材料图像的三维颗粒重建方法能够满足替代目前利用三维扫描仪实现对颗粒堆积材料进行三维重建的要求。

Claims (8)

1.一种基于单幅堆积颗粒材料图像的三维颗粒重建方法，其特征在于，步骤如下：

S1、获取堆积颗粒材料的颗粒粒径列表；

S2、基于单幅堆积颗粒材料的彩色图像提取其上颗粒的特征轮廓；

S3、利用步骤S2得到的颗粒特征轮廓，生成具有不同尺寸的颗粒表面散点模型；其中，每个颗粒表面散点模型的生成方法为：

1)根据待重建颗粒的具体轮廓要求，在颗粒轮廓数据库中随机选取两个符合要求的颗粒轮廓，并分别指定为导向轮廓和移动轮廓；

2)将导向轮廓水平放置于三维虚拟空间中，并将轮廓A旋转至其最长轴与x轴平行处；将最长轴等距划分为n个节点P_i，获取导向轮廓上与节点P_i的在x轴和z轴上坐标值相同的两个点：点M_i和点N_i之间的距离H；

3)将移动轮廓竖直放置于三维虚拟空间中，并以其重心坐标为节点，获取过其重心且与y轴平行的轴线与移动轮廓相交的两点：点O_j和点P_j之间的距离D，以计算缩放系数

并利用缩放矩阵按照比例缩放移动轮廓；接着，计算点M_i或点N_i与点O_j或点P_j之间的距离差值，通过以该距离差值为基础的平移矩阵将移动轮廓移动到导向轮廓的每个节点上，即形成颗粒表面散点模型；

S4、基于步骤三得到的各种尺寸的颗粒的表面散点模型，按照符合步骤一生成的三维堆积状态颗粒的颗粒粒径列表的要求对三维堆积状态的颗粒进行重建。

2.根据权利要求1所述的基于单幅堆积颗粒材料图像的三维颗粒重建方法，其特征在于，步骤S1的具体实施步骤为：

S101、对堆积颗粒材料或部分堆积颗粒材料进行颗粒分析试验，得到原始颗粒级配曲线；

S102、设定待生成颗粒的总个数N，以及其极限粒径mind和maxd，将原始级配曲线分为若干个粒组，

S103、根据颗粒分析试验后得到的每个粒组所占体积百分比，计算出每个粒组中待生成颗粒的总个数；

S104、以每个粒组中颗粒直径的上、下界的粒径为界限，根据步骤S103计算出的该粒组中待生成颗粒的个数，线性插值出该粒组内粒径从小到大且满足步骤S101得到的颗粒级配曲线粒径变化趋势的每个颗粒的粒径；

S105、保存所有粒组内所有颗粒的粒径，生成颗粒粒径列表。

3.根据权利要求1所述的基于单幅堆积颗粒材料图像的三维颗粒重建方法，其特征在于，步骤S2的具体实施步骤为：

S201、将单幅堆积颗粒材料的彩色图像转变为灰度图像；

S202、采用Sobel算子分别对灰度图像从水平和竖直两个方向进行滤波，并求取其梯度幅值矩阵；

S203、在经过步骤S201处理得到的灰度图像中分别对背景对象和目标对象进行标记；其中，目标对象指图像中的颗粒所在区域，背景对象指图像中颗粒所在区域与非颗粒所在区域之间的分界线；

S204、对经过步骤S202处理得到的梯度幅值矩阵进行修正，即将目标对象和背景对象对应的点在梯度幅值矩阵中修正为一个最小值Inf，完成对梯度幅值矩阵的修正；

S205、从经过步骤S204得到的图像中提取颗粒的轮廓，并利用MPP算法对颗粒的轮廓进行优化，并计算每一个颗粒轮廓的圆度、椭圆度和综合形状指数；

S206、根据步骤S205计算得到的每个轮廓的圆度、椭圆度和综合形状指数，筛选出符合预期的颗粒轮廓，并生成颗粒轮廓数据库。

4.根据权利要求3所述的基于单幅堆积颗粒材料图像的三维颗粒重建方法，其特征在于，在步骤S203中，对目标对象进行标记的具体实施步骤为：

1)以经过步骤S201处理得到的灰度图像作为掩模图像，以对掩模图像进行腐蚀运算后的图像作为样本图像，进行第一次重建运算；

2)对经过第一次重建运算后得到的图像进行求补运算，并将经过求补运算后得到的图像作为掩模图像；对经过第一次重建运算后得到的图像依次进行膨胀和求补运算，并将结果图像作为样本图像，进行第二次重建运算；然后对经过第二次重建运算后得到的图像进行求补运算；

其中，重建运算为利用样本图像f重建掩模图像g为Rg(f)，具体为：获取样本图像f；创建结构元素：B＝ones(3)；重复f_k+1＝(f_k+B)∩g，直至f_k+1＝f_k；其中，样本图像f为掩模图像g的一个子集；求补运算过程为对图像的灰度值进行取反运算，如果像素灰度为0，求补运算后就是255，如果像素灰度为10，求补后为245，也就是原始像素灰度值与其补集灰度值之和为255；

3)利用区域连通性经过步骤1)和步骤2)得到的图像进行处理：确定图像内由灰度值为0的像素包围出的多个局部区域，并将每个局部区域内出现的极大值的像素标记为1，其余像素标记为0，得到完成对目标对象进行标记的二值图像；

4)对经过上述步骤3)处理得到图像进行先开运算再闭运算的处理，完成对目标图像的标记。

5.根据权利要求3所述的基于单幅堆积颗粒材料图像的三维颗粒重建方法，其特征在于，在步骤S203中，对背景对象进行标记的具体实施步骤为：

1)通过自动阈值分割算法对经过步骤S201处理得到的灰度图像进行处理，得到该灰度图像的最佳分割阈值；

2)根据经过步骤1)得到的最佳分割阈值T^*对灰度图像进行二值化处理：将灰度图像中灰度值为大于T^*的像素位置处的数值定义为1，即标识为目标对象；将灰度图像中灰度值为小于T^*的像素位置处的数值定义为0，即标识为非目标对象区域；

3)对经过步骤2)处理得到的二值图像进行距离变换，得到距离矩阵；

4)将距离矩阵代入watershed函数中进行分水岭变换后，所得分水岭脊线即为背景图像。

6.根据权利要求1所述的基于单幅堆积颗粒材料图像的三维颗粒重建方法，其特征在于，在步骤S3中，每个颗粒表面散点模型生成前，对导向轮廓和移动轮廓进行如下处理：将导向轮廓和移动轮廓从笛卡尔直角坐标系下的坐标转换为极坐标系下的坐标，然后对极坐标列表依次进行等距线性插值和采用高斯加权移动平均法去除波动点的处理后再转换为笛卡尔直角坐标系下的坐标。

7.根据权利要求1所述的基于单幅堆积颗粒材料图像的三维颗粒重建方法，其特征在于，在步骤S3的颗粒表面散点模型生成过程中，导向轮廓上与节点P_i的在x轴上坐标值相同的两个点之间的距离H的具体计算方法为：将轮廓A水平放置于三维虚拟空间Ω中，并去除其中重复的坐标点；计算轮廓A最长轴长度L以及其与x轴的夹角ω；根据夹角ω将轮廓A旋转至其最长轴与x轴平行处，并将最长轴等距划分为n个节点P_i，每个节点给出其对应坐标；将轮廓A以最长轴的端点为界限分为A₁和A₂两个部分，然后根据节点P_i的横坐标值，找出在A₁中距离其最近点A、B，在A₂中距离其最近点A’、B’，由A、B和A’、B’线性插值出导向轮廓上与节点P_i在x轴与z轴上坐标值相同的点M_i和N_i的坐标，进而根据点M_i和点N_i的y坐标值，计算点M_i和点N_i之间的距离记为H。

8.根据权利要求1所述的基于单幅堆积颗粒材料图像的三维颗粒重建方法，其特征在于，在步骤S3的颗粒表面散点模型生成过程中，将移动轮廓竖直放置于三维虚拟空间Ω中，并以其重心G坐标为节点；经过重心G作一条平行于y轴的轴线，该轴线与移动轮廓有两个交点：点O_j和点P_j；找到y坐标值较大的交点P_j，并以该交点找到在轮廓B中距离其最近的两个点：点C和点D，由点C和点D插值出位于移动轮廓上且与重心G具有相同x坐标值和z坐标值的点在y轴上坐标值，同理，采用相同方法获得点O_j的坐标，根据点O_j和点P_j在y轴上的坐标值的差值，即得到移动轮廓上过重心G作一条平行于y轴的轴线的长度D；计算缩放系数

利用缩放矩阵按照比例缩放移动轮廓；接着，计算点M_i或点N_i与点O_j或点P_j之间的距离差值，通过以该距离差值为基础的平移矩阵将移动轮廓移动到导向轮廓的每个节点上，颗粒表面散点模型形成。

Images