一种图像重构后的网格化方法

1.本发明属于图像处理技术领域，具体涉及一种拍摄、ct扫描、核磁、电子显微等方式获取的物体数值图像重构后的网格化方法，包括图像处理、重构、网格化三大部分内容。


背景技术：

2.在物理实验中，存在实验室样品制备困难、设备制造或实验成本高、微观机理分析困难等问题，随着计算机技术的飞速发展，数值模拟技术逐渐成为与理论分析、实验室实验、现场实验同等重要的实验分析方法。而几何网格的构建是数值模拟的基础，网格的类型、连接方式、精度等决定了数值模拟的本构方程的迭代传递、节点间计算数值插值、计算精度等。
3.随着基础科学、工程科学等的逐渐深入，逐渐要求对实物实际结构下的固体力学、流体流动、传热、传质、化学反应、电磁特性等进行研究，以贴合实际工程复杂性需要，数值模拟技术在微观机理揭示方面的优势得以发挥。过去采用各种数学函数的方式构建复杂的微观几何模型，但仍然与真实的物体相去甚远，因此近年来兴起的技术利用拍摄、ct扫描、核磁、电子显微等技术手段获取真实物体断层图像，再通过重构技术实现真实物体的数值化，能够解决数值模拟对真实物体结构的需求问题。但要用于数值模拟实验研究尚且需要对构建的数值模型进行必要的网格化处理，以便于后续数值模拟软件几何模型的构建与迭代计算。
4.中国专利文献cn 114494637 a于2022年05月13日公开了一种基于结构体矩阵的砂岩三维真实模型重构方法，它包括运用x-射线成像设备获取砂岩标准试件几何拓扑结构图像数据；将图像数据转换成离散三维数字化矩阵；运用chan-vese模型对砂岩数字矩阵进行两相分割；建立结构体矩阵模式库：基于砂岩数据的结构体矩阵确立对应的中心体素单元模式和坐标信息；通过结构体矩阵反演得到几何拓扑结构数据和线性插值，得到砂岩两相微观结构的点坐标和方向矩阵，从而实现砂岩三维模型重构过程和结果可视化。
5.该专利仅停留在砂岩三维真实模型重构以完成数字岩心的可视化，没有对数字岩心进行网格化处理，不能进行后续科研分析。现有技术中，尽管avizo图像处理软件在提取孔隙网络模型(一种球-棍模型)过程中提供了网格化处理，但该软件生成的网格由于格式、网格精度等问题很难适用于其他商业化软件。而在电学特性、模量计算等领域通常利用开源软件对数字岩心直接进行相关计算，不仅需要研究者对网格进行节点连接编程处理，还需要对数值模拟部分进行编程，这不利于固体力学、流体流动、传热等多场耦合数值模拟研究。因此，急需提供一种基于图像的网格化处理方法，提高数值模型的适用性，用于微观尺度的多场耦合数值模拟研究，


技术实现要素：

6.针对现有技术存在的问题，本发明所要解决的技术问题就是提供一种图像重构后的网格化方法，它能解决数值模型网格化的问题，建立实物样品的图像处理、重构、网格化
一套完整流程。
7.本发明所要解决的技术问题是通过这样的技术方案实现的，它包括
8.步骤1、对物体样品进行ct扫描、核磁共振、电子显微实验，按数值模拟需要获取物体表面或者内部断层的几何结构图像，并通过三维重构将扫描数据处理为不同切面的物体图像；
9.步骤2、通过图像预处理对物体断层图像进行图像切割、对比度与亮度调整、滤波降噪处理，将图像中物体以外边角、扫描噪声扰动这些不利于网格化因素预先去除；
10.步骤3、对预处理后的图像进行灰度分析，通过双阈值分割将物体分为多物相，物相包括骨架颗粒和孔隙空间，并将分割好的图像切片进行三维堆叠重构，再次进行滤波降噪处理，并保存为原始xxx.raw格式三维图像；或者进一步对三维堆叠重构后的岩心生成填充物；
11.将堆叠重构或者生成填充物的数字岩心再次切割成lengthx
×
lengthy
×
lengthz个立体单元的小岩心，按x，y，z依次递增的排列顺序将岩心数据储存在相应的岩心文件中；
12.步骤4、matlab网格子程序中导入步骤3生成的岩心文件，岩心文件为按x、y、z依次递增排列的一维数组，设置网格像素尺寸(lengthx，lengthy，lengthz)，将一维数组重新组合成三维计算矩阵，并通过三维矩阵切片图像可视化处理判断重构后的数据正确性；
13.步骤5、对计算矩阵进行网格单元编号和生成节点坐标
14.网格子程序从三维计算矩阵(1，1，1)位置到(lengthx，lengthy，lengthz)依次进行编号，每个网格单元生成8个节点，最终生成的节点矩阵大小为(lengthx+1，lengthy+1，lengthz+1)，并计算节点坐标＝节点矩阵坐标
×
分辨率，将计算好的节点坐标按单元编号对应排序；节点坐标被储存在了na(lengthx+1，lengthy+1，lengthz+1)的矩阵中，网格单元编号则按不同的物相储存在cell1，cell2，cell3中；
15.步骤6、matlab网格子程序对网格单元按骨架颗粒、孔隙空间或填充物质的物相进行面网格计算，网格单元的六个面网格的连接格式为(n1 n2 n3 n4 c1 c2)，ni为构成面网格的节点及连接顺序，ci为面网格临近网格单元；二维网格单元节点连接格式(n1 n2 c1 c2)；按不同的物相将面网格储存在facem1，facem2或facem3中；
16.步骤7、节点坐标和面网格计算完成进行网格文件编写，编写格式为标准fluent网格格式，将生成的网格文件保存为xxx.msh文件；在hypermesh中将xxx.msh文件转化为应用软件接受的网格格式。
17.本发明的技术效果是：
18.本发明能实现各类图像的重构及其网格化，实现利用图像的真实物体微观结构的固体力学、流体力学、热力学、化学、电磁等的数值模拟研究工作。
附图说明
19.本发明的附图说明如下：
20.图1为图像重构与网格化处理过程示意图；
21.(a)、ct扫描岩心的断层切片图；
22.(b)、图像预处理切片；
23.(c)、三维重构图；
24.(d)、hypermesh网格可视化；
25.(e)、图像网格化后数值模拟应用图；
26.图2为孔隙填充水合物的流程图；
27.图3为图2中的循环搜索及生成水合物的子程序流程图；
28.图4为50
×
50
×
50三维三相网格化程序运行的界面图。
具体实施方式
29.下面结合附图和实施例对本发明作进一步说明：
30.实施例，包括以下步骤：
31.步骤1、对物体样品进行ct扫描、核磁共振、电子显微实验，按数值模拟需要获取物体表面或者内部断层的几何结构图像，并通过三维重构将扫描数据处理为不同断层切片图像。
32.如图1(a)所示，用479张1000
×
1000像素的ct几何结构图像，将这些几何结构图像重构所得三维图经切割的ct断层切片图。参见文献“基于x-ray ct和有限元方法的沥青混合料三维重构与数值试验研究”，万成，华南理工大学，第20～44页，2010年4月记载了的一种ct扫描重构图像的方法。
33.步骤2、通过图像预处理对断层切片图像进行图像切割以去除图像外侧非岩心部分、对比度与亮度调整、滤波降噪处理，将扫描噪声扰动的不利于网格化因素预先去除。
34.如图1(b)所示，本实施例通过imagej、avizo等图像处理软件将断层切片图像切割成长
×
宽为331
×
331的切片序列，将样品外侧非岩心部分去除；然后对图片进行亮度、对比度等调整工作，提高岩心骨架颗粒、孔隙、填充物间的分辩度，同时对图像进行滤波降噪工作，以便于对岩心进行分割。
35.步骤3、对预处理后的图片进行灰度分析，通过双阈值分割将物体分为多物相，物相包括骨架颗粒和孔隙空间，并将分割好的图像切片按循序叠放在一起进行三维重构为三维数字岩心，再次进行滤波降噪工作，除去因阈值分割引起的不必要孔隙结构，并保存为原始xxx.raw格式三维图像。
36.参见文献“pore fractal characteristics of hydrate-bearing sands and implications to the saturated water permeability”，zhun zhang et al，jgr(journal of geophysical research)：solid earth，125，p5～6，mar，2020，(“水合物储层砂岩孔隙分形特性及其对饱和水渗透率影响”，张准，地球物理研究杂志，第125期，5～6页，2020年3月)记载的一种双阈值分割方法。所述图像分割方法，采用双阈值分割的方法将骨架颗粒和孔隙空间进行分割，分割后的图像二值化处理以便后续操作，即将骨架颗粒记为数值“0”，孔隙空间记为数值“1”。
37.如图1(c)所示，在实施例中，采用双阈值分割，通过对步骤2预处理图片的灰度值进行阈值分割，分割后将岩心分割为骨架颗粒和孔隙两个部分。分割完成后，将ct断层切片图像序列按切片排列顺序堆叠形成三维的数字岩心，将岩心数据保存为三维a1particle.raw的原始图像数据。
38.再进一步，根据孔隙填充物的生成饱和度，通过不考虑生成条件下随机地进行像素点的替换或考虑生成条件下人为控制像素点的替换(条件由具体的孔隙填充物生成机理
决定)等方式，对三维重构后生成填充物。
39.通过像素点替换计算用于生成不同饱和度的填充样品数值生成。填充物相饱和度＝填充物相像素数
÷
(填充物相像素数+剩余孔隙像素数)，该填充物生成步骤用于克服实验室生成多相多孔介质样品饱和度有限的问题，可以实现目标范围饱和度的精准控制。
40.在实施例中，通过分水岭方法分割孔隙空间，将孔隙按喉道、孔隙分割成一定大小的孔隙区块，计算孔隙区块的几何中心坐标和平均孔隙半径并储存；另外对骨架空间进行分水岭分割，获取骨架中心坐标和骨架颗粒半径。参见文献“versatile and efficient pore network extraction method using marker-based watershed segmentation”，jeff t.gostick，physical review e 2017(“基于标记分水岭分割的通用高效孔隙网络提取方法”，jeff t.gostick，物理评论e，第96期，第1～13页，2017年8月)记载了一种分水岭方法分割孔隙空间，通过计算欧氏距离、“峰”拾取、“峰”移除、邻近“峰”合并等基本操作，最终将孔隙按喉道、孔隙分割成一定大小的孔隙区块，并计算孔隙区块的几何中心坐标、平均孔隙半径并储存；另外对骨架空间进行分水岭分割，获取骨架中心坐标、骨架颗粒半径。
41.基于matlab编写孔隙填充物程序，生成填充物包括以下步骤：
42.1)、导入xxx.raw文件、孔隙信息文件、骨架信息文件，将图像数据读取为三维矩阵，形成计算矩阵a(lenx,leny,lenz)，为生成填充物提供基础的骨架和孔隙数据；
43.2)、计算孔隙搜索范围，该搜索范围为自孔隙中心向外半径为ratio*plr的“球”形范围；并以孔隙中心循环搜索孔隙像素；
44.3)、在孔隙搜索范围内的孔隙中生成填充物：孔隙填充物由孔隙中心向骨架颗粒延伸，直至填充物达到与骨架颗粒最小间隙不小于一个体素；
45.4)、一个孔隙的填充物生成完成后进入下一个孔隙搜索，所有孔隙搜索完成后计算填充物饱和度；当填充物饱和度未达到预设饱和度时，调整孔隙搜索范围和预设孔隙半径阈值plr0直至满足预设填充物饱和度；当填充物饱和度到达预设饱和度后，填充物生成完成。
46.如图2所示，填充水合物生成程序按以下步骤：
47.在s101，导入a1particle.raw文件的图像数据；
48.在s102，导入孔隙信息文件和骨架信息文件；
49.所述孔隙信息包含孔隙区块的几何中心坐标和平均孔隙半径；所述骨架信息包含骨架中心坐标和骨架颗粒半径。
50.在s103，设置三维图像的尺寸lenx,leny,lenz；
51.在s104，将s101的图像数据读取为三维矩阵，形成计算矩阵a(lenx,leny,lenz)；
52.在s105，设置预期水合物饱和度sh0和孔隙半径阈值plr0，以控制水合物生成饱和度和确定不同大小孔隙中能否生成孔隙填充物；
53.在s106，利用孔隙信息，并在以孔隙中心向外半径为ratio*plr的“球”形范围内循环搜索孔隙像素；
54.在s107，在孔隙中生成水合物，将孔隙像素值替换为“7”；
55.在s108，计算水合物饱和度sh＝水合物相体积/(水合物相体积+剩余孔隙体积)；
56.在s109，判断水合物饱和度sh是否达到预设饱和度sh0，若是，则执行s111；否则，执行s110；
57.在s110，调整孔隙搜索范围ratio*plr和预设孔隙半径阈值plr0，返回s105；
58.调整水合物极限饱和度在40％以下，满足实际孔隙填充型天然气水合物饱和度极限要求(20％～40％)。
59.在s111，孔隙填充水合物完成。
60.上述s106与s107中循环搜索及生成水合物的子程序，如图3所示，该子程序，按以下步骤：
61.在s201，按孔隙中心坐标的不同，搜索各孔隙像素；
62.在s202，判断孔隙半径是否大于阈值plr0，若是，则执行s203；否则，返回s201；
63.在s203，计算孔隙的搜索半径ratio*plr；
64.在s204，以孔隙中心坐标为“球心”，搜索半径ratio*plr的“球”形范围内的孔隙像素。本步骤可以极大的降低水合物生成循环和判断的计算量。
65.在s205，判断搜索范围内孔隙像素是否与骨架材料相接触，如是，则执行s207；否则执行s206；
66.水合物和骨架颗粒间至少保留一个像素的间隔。
67.在s206，生成填充水合物物，即将孔隙像素值替换为“7”，进入s208；
68.在s207，继续循环下一个孔隙像素；
69.在s208，判断半径ratio*plr的“球”范围内的孔隙像素是否全部搜索完成，如是，则执行s209，否则返回s204；
70.在s209，半径ratio*plr“球”范围内的孔隙像素搜索循环结束，进入下一个孔隙循环106。
71.利用图2所示的水合物生成程序，在不考虑生成物理化学等条件的前提下，通过在孔隙中间向骨架颗粒搜索孔隙像素，当孔隙空间不与骨架颗粒接触时将孔隙空间像素替换为填充物，直到填充物饱和度满足预设饱和度则停止，最终生成12.52％、18.66％、27.07％、31.39％等水合物饱和度的数值岩心。
72.为了使生成的网格满足后续comsol多场耦合数值模拟计算的内存与cpu需要(本实例多场耦合数值模拟计算机为24核cpu、128g内存)，将本步骤制得的数字岩心再次切割成50
×
50
×
50的12.5万个立体单元的小岩心(受制于comsol计算规模，难以切割成更小岩心)，按x，y，z依次递增的排列顺序将岩心数据储存在相应的岩心文件(txt文本文件)中。
73.步骤4、matlab网格子程序中导入步骤3生成的岩心文件，岩心文件是按x、y、z依次递增排列的一维数组(一维数组是指在内存中连续储存的具有相同类型的一组数据的集合)，设置网格像素尺寸(lengthx，lengthy，lengthz)，网格子程序将一维数组重新组合成三维计算矩阵，并通过ct扫描切片图像可视化处理判断步骤3重构后的数据正确性(阈值分割的正确性)。
74.说明：此处的计算矩阵与步骤s104中“三维计算矩阵”的形式上相同，但包含信息不同，步骤s104的计算矩阵主要包含孔隙和骨架信息，用于孔隙填充物的填充，本处的计算矩阵包含孔隙/骨架和生成的填充物信息，用于计算网格。
75.作为实施例，将岩心文件数据导入matlab网格子程序中，设置岩心试样的网格像素尺寸[x y z]为：[50 50 50]；三相物质像素值[m1 m2 m3]：[1 7 3]，即标记每种物相的像素值，本案例骨架颗粒被标记为“1”，填充物被标记为“7”，孔隙空间为“3”；然后输入像素
比例尺，即一个像素所代表的实际物理尺寸，本实施例中，网格计算时采用像素比例[s1 s2 s3]：[1 1 1]。网格化程序的初始化数据输入如下：
[0076]
clear all；
[0077]
size1＝input(“请输入网格像素尺寸，格式[x,y,z]:”)；％试样像素尺寸
[0078]
lenx＝size1(1)；leny＝size1(2)；lenz＝size1(3)；
[0079]
im＝load(“hydrate.txt”)；％读取ct岩心数据
[0080]
mass＝input(“请输入三相物质像素值，格式为[m1,m2,m3]，单项两项充0：”)；％不同相的像素标定
[0081]
m1＝mass(1)；m2＝mass(7)；m3＝mass(3)；
[0082]
scal＝input(“请输入比例尺(一像素代表多少物理尺寸)[s1,s2,s3]，默认值1：”)；％物理尺寸比例
[0083]
s1＝scal(1)；s2＝scal(1)；s3＝scal(1)；
[0084]
sum1＝sum(im(:)＝＝m1)；
[0085]
sum2＝sum(im(:)＝＝m2)；
[0086]
sum3＝sum(im(:)＝＝m3)；
[0087]
a＝reshape(im,lenx,leny,lenz)；
[0088]
subplot(1,2,1)；
[0089]
imshow(a(:,:,2),[])；
[0090]
运行网格化程序，运行界面如图4所示：matlab网格子程序首先会将读取的一维(125000)岩心文件重新组合成三维a(50,50,50)的岩心计算矩阵。并显示了计算矩阵a的第2张切片，该切片显示结果与原始切片孔隙和骨架颗粒符合良好，验证了重构后的数据正确性。
[0091]
步骤5、对计算矩阵进行网格单元编号和生成节点坐标
[0092]
网格单元是指组成网格的最小的结构，由步骤3按50
×
50
×
50切割小岩心，再由步骤4设置网格像素尺寸为[50 50 50]，则1个小岩心对应1个网格单元。
[0093]
网格子程序从三维计算矩阵(1，1，1)位置到(50，50，50)依次进行编号，编号按x,y,z递增对正方体网格单元的八个顶点从1到132651＝51
×
51
×
51依次进行编号，每个网格单元生成8个节点(8节点围成结构化网格)，最终生成的节点矩阵大小为(51，51，51)，并计算节点坐标＝节点矩阵坐标
×
分辨率，即节点坐标(x，y，z)＝index(x，y，z)
×
scale(x0，y0，z0)，将计算好的节点坐标按网格单元编号对应排序。
[0094]
matlab网格子程序按x,y,z循环进行计算网格节点和单元编号如下：
[0095][0096]
[0097]
节点坐标被储存在了na(51，51，51)的矩阵中，网格单元编号则按不同的物相(骨架、孔隙、填充物)储存在cell1，cell2，cell3(分别对应物质像素值[m1 m2 m3])中，其大小随实际不同物相的单元数量决定。
[0098]
步骤6、matlab网格子程序对网格单元按骨架颗粒、孔隙空间、填充物质三物相进行面网格计算，所述面网格是指构成三维网格单元的六个面和二维网格单元的四个边；使不同物相网格相互独立便于后续数值模拟处理。
[0099]
面网格的节点连接格式为(n1，n2，n3，n4，c1，c2)的结构化网格，ni为构成面网格的节点及连接顺序，ci为面网格临近网格单元。例如网格单元1(8节点：1,2,4,5,7,8,10,11(按坐标递增排序))和网格单元2(8节点：2,3,5,6,8,9,11,12)是相邻的两个单元，则面(2,5,8,11)与临近单元1和2，其四个节点按正方形四顶点依次连接顺序为(2,5,11,8)，该网格记作(2,5,11,8,1,2)。对网格单元2来说，该面网格的节点连接格式为(2,5,11,8,2,1)。
[0100]
在数字岩心表面的二维网格单元提取四个边计算，二维网格单元节点连接格式是去掉n3，n4。
[0101]
按不同的物相将面网格储存在facem1，facem2，facem3中，直至网格单元所有的面网格计算完成。
[0102]
步骤7、在节点坐标和面网格计算完成后，进行网格文件编写，编写格式为标准fluent网格格式，包括网格维度信息、节点信息、面网格单元信息、边界面网格信息和网格名称，将生成的网格文件保存为xxx.msh文件。在hypermesh中将xxx.msh文件转化为comsol、abaqus、ansys等数值模拟网格即可进行后续数值模拟工作。
[0103]
网格维度信息有二维和三维，节点信息包括节点坐标和节点数量，面网格单元信息包括面网格连接，分类和数量。
[0104]
将节点坐标、面网格信息按fluent网格格式(该格式已公开)编写到.msh文件中，网格文件编写如下：
[0105]
fid＝fopen(strcat('outname.msh'),'w')；％网格名称声明
[0106]
fprintf(fid,'％s％d％s\n','(',0,'"created by wang form ct raw date")')；％网格作者信息声明
[0107]
fprintf(fid,'％s％d％d％s\n','(',2,3,')')；％网格几何信息
[0108]
fprintf(fid,'％s％d％s\n','(',0,'"note section:")')；％％节点信息
[0109]
fprintf(fid,'％s％d％s％d％2x％2x％d％d％s\n','(',10,'(',0,1,(lenx+1)*(leny+1)*(le nz+1),0,3,'))')；
[0110]
fprintf(fid,'％s％d％s％d％2x％2x％d％d％s\n','(',10,'(',1,1,(lenx+1)*(leny+1)*(le nz+1),1,3,')')；
[0111]
fprintf(fid,'％s\n','(')；
[0112]
fprintf(fid,'％d％d％d\n',(i-1)*s1,(j-1)*s2,(k-1)*s3)；
[0113]
fprintf(fid,'％s\n','))')；％％
[0114]
fprintf(fid,'％s％d％s\n','(',0,'"faces:")')；％％面网格声明
[0115]
fprintf(fid,'％s％d％s％2x％2x％2x％2x％2x％s\n','(',13,'(',0,1,lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+lenf7,0,0,'))')；％face
[0116]
％bar2＝waitbar(0,'interior1编写中')；
[0117]
fprintf(fid,'％s％d％s\n','(',0,'"interior faces of zone void")')；
[0118]
if lenf1＝＝0
[0119]
else
[0120]
fprintf(fid,'％s％d％s％2x％2x％2x％2x％2x％s\n','(',13,'(',7,1,lenf1,2,4,')')；％inte rior face
[0121]
fprintf(fid,'％s\n','(')；
[0122]
％for i＝1:lenf1
[0123]
fprintf(fid,'％x％x％x％x％x％x\n',interior1)；
[0124]
str1＝['网格计算中
[0125]
',num2str(((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenf1)*100/((lenx+1)*(le ny+1)*(lenz+1)+lenx*leny*lenz+lenx*leny*lenz*3+lenx*leny+lenx*lenz+leny*le nz)),'％']；
[0126]
waitbar(((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenf1)/((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenx*leny*lenz*3+lenx*leny+lenx*lenz+leny*lenz),bar1,str1)；
[0127]
％str2＝['interior1编写中...',num2str((i)/lenf1),'％']；
[0128]
％waitbar((i)/lenf1,bar2,str2)；
[0129]
％end
[0130]
fprintf(fid,'％s％s\n',')',')')；
[0131]
end％％
[0132]
……
[0133]
％bar2＝waitbar(0,'wall编写中')；％％边界网格声明
[0134]
if lenf7＝＝0
[0135]
else
[0136]
fprintf(fid,'％s％d％s\n','(',0,'"wall")')；
[0137]
fprintf(fid,'％s％d％s％2x％2x％2x％2x％2x％s\n','(',13,'(',13,lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+1,lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+lenf7,3,4,')')；％wall
[0138]
fprintf(fid,'％s\n','(')；
[0139]
％for i＝1:lenf7
[0140]
fprintf(fid,'％2x％2x％2x％2x％2x％2x\n',wall)；
[0141]
str1＝['网格计算中
[0142]
',num2str(((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+lenf7)*100/((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenx*leny*lenz*3+lenx*leny+lenx*lenz+leny*lenz)),'％']；
[0143]
waitbar(((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenf1+lenf2+lenf3+lenf4+lenf5+lenf6+lenf7)/((lenx+1)*(leny+1)*(lenz+1)+lenx*leny*lenz+lenx*leny*lenz*3+lenx*leny+lenx*lenz+leny*lenz),bar1,str1)；
[0144]
％ str2＝['walls编写中...',num2str((i)/lenf7),'％']；
[0145]
％ waitbar((i)/lenf7,bar2,str2)；
[0146]
％end
[0147]
fprintf(fid,'％s％s\n',')',')')；
[0148]
end％％
[0149]
最终生成的网格文件包含网格声明、网格维度信息、节点坐标、面网格信息、网格结构信息等，网格文件如下：
[0150]
(0"created by wang form ct raw date")％作者声明
[0151]
(2 3)％网格维度信息
[0152]
(0"note section:")％节点声明
[0153]
(10(0 1 2062b 0 3))％％节点信息
[0154]
(10(1 1 2062b 1 3)
[0155]
(
[0156]
0 0 0
[0157]
1 0 0
[0158]
2 0 0
[0159]
3 0 0
[0160]
4 0 0
[0161]
5 0 0
[0162]
6 0 0
[0163]
7 0 0
[0164]
8 0 0
[0165]
9 0 0
[0166]
……
[0167]
46 50 50
[0168]
47 50 50
[0169]
48 50 50
[0170]
49 50 50
[0171]
50 50 50
[0172]
))％％
[0173]
0"faces:")％％面网格信息
[0174]
(13(0 1 5d624 0 0))
[0175]
(0"interior faces of zone void")(13(7 1 1802d 2 4)(
[0176]
a60 a93 14bc 1489 2 1
[0177]
1488 1489 14bc 14bb 359 1a61 a94 14bd 148a 3 2a93 a94 14bd 14bc 18 2
[0178]
1489 148a 14bd 14bc 35a 2a62 a95 14be 148b 4 3a94 a95 14be 14bd 19 3
[0179]
148a 148b 14be 14bd 35b 3
……
[0180]
))
[0181]
……
[0182]
(0"wall")％％边界网格(13(d 59b8d 5d624 3 4)(
[0183]
1 34a5d a2a 0af30
[0184]
1 2a2b a2a 0af30
[0185]
1 2 35 34 0af30
[0186]
2 3a2c a2b 0af31
[0187]
2 3 36 35 0af31
[0188]
3 4a2d a2c 0af32
[0189]
3 4 37 36 0af32
[0190]
……
[0191]
))
[0192]
(0"cells:")％％网格单元声明(12(0 1 1e848 0 0))
[0193]
(12(4 1 99ea 1 4))
[0194]
(12(5 99eb af2f 1 4))
[0195]
(12(6af30 1e848 1 4))
[0196]
(0"zone sections:")％％网格各相信息声明
[0197]
(45(4fluid void)())
[0198]
(45(5solid hydrate)())
[0199]
(45(6solid solid)())
[0200]
(0"faces sections:")％％面网格信息声明
[0201]
(45(7interior vf)())
[0202]
(45(8interior hf)())
[0203]
(45(9interior sf)())
[0204]
(45(a interface vs)())
[0205]
(45(b interface vh)())
[0206]
(45(d wall walls)())％％边界网格信息声明
[0207]
图1(d)为hypermesh网格处理软件的网格显示图，本实施例在hypermesh中对网格进行了尺寸调整和网格格式转换工作，并将转换后的.nas格式网格文件导入进comasol中开展后续的数值模拟工作。
[0208]
本发明的核心点是图像重构后的网格化，这种网格化的具体应用如图1(e)所示，为comsol多场耦合数值模拟应用的几何显示图，通过comsol导入网格文件并进行了流-固耦合数值模拟研究。

技术特征：
1.一种图像重构后的网格化方法，其特征是，包括以下步骤：步骤1、对物体样品进行ct扫描、核磁共振、电子显微实验，按数值模拟需要获取物体表面或者内部断层的几何结构图像，通过三维重构将扫描数据处理为不同切面的物体图像；步骤2、通过图像预处理对物体断层图像进行图像切割、对比度与亮度调整、滤波降噪处理，将图像中物体以外边角、扫描噪声扰动这些不利于网格化因素预先去除；步骤3、对预处理后的图像进行灰度分析，通过双阈值分割将物体分为多物相，物相包括骨架颗粒和孔隙空间，并将分割好的图像切片进行三维堆叠重构，再次进行滤波降噪处理，并保存为原始xxx.raw格式三维图像；或者进一步对三维堆叠重构后的岩心生成填充物；将堆叠重构或者生成填充物的数字岩心再次切割成lengthx
×
lengthy
×
lengthz个立体单元的小岩心，按x，y，z依次递增的排列顺序将岩心数据储存在相应的岩心文件中；步骤4、matlab网格子程序中导入步骤3生成的岩心文件，岩心文件为按x、y、z依次递增排列的一维数组，设置网格像素尺寸(lengthx，lengthy，lengthz)，将一维数组重新组合成三维计算矩阵，并通过三维矩阵切片图像可视化处理判断重构后的数据正确性；步骤5、对计算矩阵进行网格单元编号和生成节点坐标网格子程序从三维计算矩阵(1，1，1)位置到(lengthx，lengthy，lengthz)依次进行编号，每个网格单元生成8个节点，最终生成的节点矩阵大小为(lengthx+1，lengthy+1，lengthz+1)，并计算节点坐标＝节点矩阵坐标
×
分辨率，将计算好的节点坐标按单元编号对应排序；节点坐标被储存在了na(lengthx+1，lengthy+1，lengthz+1)的矩阵中，网格单元编号则按不同的物相储存在cell1，cell2，cell3中；步骤6、matlab网格子程序对网格单元按骨架颗粒、孔隙空间或填充物质的物相进行面网格计算，网格单元的六个面网格的连接格式为(n1 n2 n3 n4 c1 c2)，ni为构成面网格的节点及连接顺序，ci为面网格临近网格单元；二维网格单元节点连接格式(n1 n2 c1 c2)；按不同的物相将面网格储存在facem1，facem2或facem3中；步骤7、节点坐标和面网格计算完成进行网格文件编写，编写格式为标准fluent网格格式，将生成的网格文件保存为xxx.msh文件；在hypermesh中将xxx.msh文件转化为应用软件接受的网格格式。2.根据权利要求1所述的图像重构后的网格化方法，其特征是，在步骤3中，岩心生成填充物按以下步骤：1)、导入xxx.raw文件、孔隙信息文件、骨架信息文件，将图像数据读取为三维矩阵，形成计算矩阵a(lenx,leny,lenz)，为生成填充物提供基础的骨架和孔隙数据；2)、计算孔隙搜索范围，该搜索范围为自孔隙中心向外半径为ratio*plr的“球”形范围；并以孔隙中心循环搜索孔隙像素；3)、在孔隙搜索范围内的孔隙中生成填充物：孔隙填充物由孔隙中心向骨架颗粒延伸，直至填充物达到与骨架颗粒最小间隙不小于一个体素；4)、一个孔隙的填充物生成完成后进入下一个孔隙搜索，所有孔隙搜索完成后计算填充物饱和度；当填充物饱和度未达到预设饱和度时，调整孔隙搜索范围和预设孔隙半径阈值plr0直至满足预设填充物饱和度；当填充物饱和度到达预设饱和度后，填充物生成完成。3.根据权利要求1或2所述的图像重构后的网格化方法，其特征是：所述lengthx、
lengthy、lengthz均为50。

技术总结
本发明公开了一种图像重构后的网格化方法，包括以下步骤，1、对物体样品几何结构图像进行三维重构，将扫描数据处理为不同切面的物体图像；2、图像预处理；3、将物体分割为多物相，将分割好的图像切片进行三维堆叠重构，或者进一步对三维堆叠重构后的岩心生成填充物；数字岩心再次切割成lengthX


技术研发人员：罗永江 王兴 张振宇 刘子豪 郭坤勇 胡千庭
受保护的技术使用者：重庆大学
技术研发日：2023.06.26
技术公布日：2023/9/13