CN114565658A - 基于ct技术的孔隙尺寸计算方法及装置 - Google Patents
基于ct技术的孔隙尺寸计算方法及装置 Download PDFInfo
- Publication number
- CN114565658A CN114565658A CN202210047008.8A CN202210047008A CN114565658A CN 114565658 A CN114565658 A CN 114565658A CN 202210047008 A CN202210047008 A CN 202210047008A CN 114565658 A CN114565658 A CN 114565658A
- Authority
- CN
- China
- Prior art keywords
- pore
- voxel
- dimensional
- radius
- sample
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 239000011148 porous material Substances 0.000 title claims abstract description 254
- 238000005516 engineering process Methods 0.000 title claims abstract description 34
- 238000004364 calculation method Methods 0.000 title claims description 15
- 239000007787 solid Substances 0.000 claims abstract description 54
- 238000012545 processing Methods 0.000 claims abstract description 13
- 238000002591 computed tomography Methods 0.000 claims description 31
- 238000010586 diagram Methods 0.000 claims description 13
- 239000003086 colorant Substances 0.000 claims description 9
- 238000003384 imaging method Methods 0.000 claims description 8
- 239000002689 soil Substances 0.000 claims description 8
- 238000003491 array Methods 0.000 claims description 6
- 238000004891 communication Methods 0.000 claims description 6
- 238000000034 method Methods 0.000 abstract description 20
- 239000011800 void material Substances 0.000 description 4
- 238000012216 screening Methods 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 230000009191 jumping Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000009659 non-destructive testing Methods 0.000 description 1
- 239000000758 substrate Substances 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
- G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
- G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
- G01N23/046—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. computed tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2223/00—Investigating materials by wave or particle radiation
- G01N2223/03—Investigating materials by wave or particle radiation by transmission
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2223/00—Investigating materials by wave or particle radiation
- G01N2223/10—Different kinds of radiation or particles
- G01N2223/101—Different kinds of radiation or particles electromagnetic radiation
- G01N2223/1016—X-ray
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2223/00—Investigating materials by wave or particle radiation
- G01N2223/40—Imaging
- G01N2223/419—Imaging computed tomograph
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2223/00—Investigating materials by wave or particle radiation
- G01N2223/60—Specific applications or type of materials
- G01N2223/649—Specific applications or type of materials porosity
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10056—Microscopic image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Geometry (AREA)
- Life Sciences & Earth Sciences (AREA)
- Biochemistry (AREA)
- Pulmonology (AREA)
- Radiology & Medical Imaging (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- General Health & Medical Sciences (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
- Image Generation (AREA)
Abstract
本发明提供基于CT技术的孔隙尺寸计算方法及装置，能够形象直观的描述和模拟试样内部孔隙的形貌尺寸，并且精确地计算孔隙尺寸，从而更加准确地获取试样的真实孔隙信息。方法包括：步骤1.将试样原始图像进行图像处理和三维重构，得到由nx×ny×nz个体素构成的试样三维图像，将图像沿x方向按顺序以二维切片的图片保存，得到试样的一系列二维切片灰度图；步骤2.将切片灰度图进行二值化处理，得到仅包含试样固体体素和孔隙体素的二值化图片；步骤3.将所有二值化后的切片图进行叠加，得到试样的三维二值图像，该图像为(nx×ny×nz)的三维数组，记为数组A；步骤4.对孔隙空间中每一个孔隙体素，依次求解孔隙半径。
Description
技术领域
本发明属于测量技术领域，具体涉及基于CT技术的孔隙尺寸计算方法及装置。
技术背景
随着科学技术的发展和研究手段的不断丰富，X射线断层扫描成像(X-rayComputed Tomography，X-CT)等CT三维成像技术逐渐成为获取材料内部信息的有效手段。近十年来，利用CT技术获取并分析多孔介质的细观或微观结构特征成为一种趋势，其在三维成像、无损检测、科研等各领域均有广泛应用。
利用CT技术获取多孔介质的三维孔隙结构后，需要对孔隙结构进行分析，获取其基本参数，比如孔隙直径及其分布等。目前基于CT图像计算孔隙直径的方法主要有两类，一类是基于二维的图片，将孔隙区域转化为等面积的圆，孔隙区域的直径即为圆的直径。另一类是基于三维图像，将孔隙结构简化为由规则的球状孔和圆管状孔喉构成，其中大的孔按照一定的体积换算规则简化为球形，大孔的直径极为球的直径；而小孔隙也按照体积换算规则简化为圆管，小孔隙直径即为对应圆管直径。由于真实多孔介质中孔隙结构并非球形和圆管，因此有些方法还会在此基础上根据孔隙形状增加正方形、三角形等形状，再按照一定规则，比如赋值形状因子，将其它形状等效为球和圆管。然而，天然的多孔介质中孔隙结构非常复杂，形状极不规则，空间变异性大，以上这些方法根据当量面积和体积，将孔隙结构用规则形状进行等效，对孔隙结构做了简化处理，无法表征孔隙的真实形貌尺寸。
发明内容
本发明是为了解决上述问题而进行的，目的在于提供一种基于CT技术的孔隙尺寸计算方法及装置，能够形象直观的描述和模拟试样内部孔隙的形貌尺寸，并且精确地计算孔隙尺寸，从而更加准确地获取试样的真实孔隙信息。
本发明为了实现上述目的，采用了以下方案：
<方法>
本发明提供了基于CT技术的孔隙尺寸计算方法，其特征在于，包括以下步骤：
步骤1.使用三维成像显微***对待计算的多孔介质试样进行CT扫描，获取试样原始投影图像；将试样原始图像进行图像处理和三维重构，得到由nx×ny×nz个体素构成的试样三维图像，体素表示三维图像中的一个正方体结构(边长为一个像素)的三维像素点；将图像沿x方向按顺序以二维切片的图片保存，得到试样的一系列二维切片灰度图，每个切片图上包含ny×nz个体素；
步骤2.将切片灰度图进行二值化处理，得到仅包含试样固体体素和孔隙体素的二值化图片，将固体体素用数值a0表示，孔隙体素用数值a1表示，则每张切片图都可以用由a0和 a1构成的二维数组表示，每个二维数组包含ny×nz个数值元素，分别与ny×nz个体素一一对应；
步骤3.将所有二值化后的切片图进行叠加，得到试样的三维二值图像，该三维二值图像为(nx×ny×nz)的三维数组，记为数组A；数组A中的元素记为ip[i][j][k]，其中i，j，k分别代表元素所对应的体素在x，y，z方向的坐标，ip[i][j][k]＝a1时，对应的体素为孔隙； ip[i][j][k]＝a0时，对应的体素为固体；
步骤4.对孔隙空间中每一个ip[i][j][k]＝a1的孔隙体素，按从左至右、由下至上、由后至前的顺序，依次求解孔隙半径，具体包括以下子步骤：
步骤4-1.将当前待求解的孔隙体素的半径d初始值设定为1个体素尺寸，d＝1，以当前待求解的孔隙体素为中心，d为半径作球；
步骤4-2.对球体区域中的元素进行判断：
若球体区域中存在ip[i][j][k]＝a0的元素，表明所作球超出试样范围或含有固体体素，说明d是以当前待求解的孔隙体素为球心、与试样边界或固体边界内切球的最大半径，则球心体素对应的孔隙半径为d，令元素ip[i][j][k]＝d；然后，进入步骤4-3；
若球体区域中不存在ip[i][j][k]＝a0的元素，表明以d为半径所作球未超出试样范围并且不含有固体体素，进入步骤4-4；
步骤4-3.判断当前已求解的孔隙体素是否为最后一个待求解的孔隙体素，若是，则求解结束，若否，则将下一个孔隙体素设定为当前待求解的体素，返回步骤4-1继续求解；
步骤4-4.令d＝d+1，将球的半径扩大1个体素尺寸，然后返回步骤4-2；
通过以上步骤4-1～步骤4-2得到显示试样所有孔隙体素的孔隙半径信息的新的三维数组B，该数组B中ip[i][j][k]＝a0时表示该点为固体体素；ip[i][j][k]不等于a0时，其值代表该孔隙体素所对应的孔隙半径。
具体地，本发明提供的基于CT技术的孔隙尺寸计算方法，还可以具有以下特征：在步骤2中，a0＝0，a1＝1；即将固体体素用0表示，孔隙体素用1表示，则每张切片图都可以用一个元素由0和1构成的二维数组表示。
优选地，本发明提供的基于CT技术的孔隙尺寸计算方法，还可以具有以下特征：在步骤3中，统计孔隙体素的数量为np，统计固体体素数量为ns，按下列公式计算孔隙度
优选地，本发明提供的基于CT技术的孔隙尺寸计算方法，还可以包括：步骤5.为反映土壤孔隙性质的差异，依据步骤4得到的数组B中ip[i][j][k]将孔隙按半径大小进行分类区分、统计数量和确定分布情况。
优选地，本发明提供的基于CT技术的孔隙尺寸计算方法，还可以具有以下特征：步骤6.根据各体素与元素的对应关系和步骤4得到的数组B，生成反映真实孔隙信息的试样孔隙尺寸及分布情况三维模型图。
优选地，本发明提供的基于CT技术的孔隙尺寸计算方法，还可以具有以下特征：步骤 7.根据孔隙半径大小不同，在对应的试样孔隙尺寸及分布情况三维模型图上用不同颜色进行不同半径级别的区分，并能够整体或者分级分别显示不同半径级别的孔隙。
优选地，本发明提供的基于CT技术的孔隙尺寸计算方法，还可以具有以下特征：步骤 8.按照以下公式计算孔隙的平均半径：
式中，
为所有孔隙的平均半径，nx为试样的像素尺寸，
为半径为Dj的孔隙数目。
<装置>
进一步，本发明还提供了一种自动实现上述<方法>的基于CT技术的孔隙尺寸计算装置，其特征在于，包括：
扫描重构部，使用三维成像显微***对待计算的多孔介质试样进行CT扫描，获取试样原始投影图像；将试样原始图像进行图像处理和三维重构，得到由nx×ny×nz个体素构成的试样三维图像，体素表示三维图像中的一个三维像素点；将图像沿x方向按顺序以二维切片的图片保存，得到试样的一系列二维切片灰度图，每个切片图上包含ny×nz个体素；
二值化部，将切片灰度图进行二值化处理，得到仅包含试样固体体素和孔隙体素的二值化图片，将固体体素用数值a0表示，孔隙体素用数值a1表示，则每张切片图都可以用由a0和a1构成的二维数组表示，每个二维数组包含ny×nz个数值元素，分别与ny×nz个体素一一对应；
三维数组生成部，将所有二值化后的切片图进行叠加，得到试样的三维二值图像，该三维二值图像为(nx×ny×nz)的三维数组，记为数组A；数组A中的元素记为ip[i][j][k]，其中i，j，k分别代表元素所对应的体素在x，y，z方向的坐标，ip[i][j][k]＝a1时，对应的体素为孔隙；ip[i][j][k]＝a0时，对应的体素为固体；
空隙半径求解部，对孔隙空间中每一个ip[i][j][k]＝a1的孔隙体素，依次求解孔隙半径，具体包括以下子步骤：
步骤4-1.将当前待求解的孔隙体素的半径d初始值设定为1个体素尺寸，d＝1，以当前待求解的孔隙体素为中心，d为半径作球；
步骤4-2.对球体区域中的元素进行判断：
若球体区域中存在ip[i][j][k]＝a0的元素，表明所作球超出试样范围或含有固体体素，说明d是以当前待求解的孔隙体素为球心、与试样边界或固体边界内切球的最大半径，则球心体素对应的孔隙半径为d，令元素ip[i][j][k]＝d；然后，进入步骤4-3；
若球体区域中不存在ip[i][j][k]＝a0的元素，表明以d为半径所作球未超出试样范围并且不含有固体体素，进入步骤4-4；
步骤4-3.判断当前已求解的孔隙体素是否为最后一个待求解的孔隙体素，若是，则求解结束，若否，则将下一个孔隙体素设定为当前待求解的体素，返回步骤4-1继续求解；
步骤4-4.令d＝d+1，将球的半径扩大1个体素尺寸，然后返回步骤4-2；
通过以上步骤4-1～步骤4-2得到显示试样所有孔隙体素的孔隙半径信息的新的三维数组B，该数组B中ip[i][j][k]＝a0时表示该点为固体体素；ip[i][j][k]不等于a0时，其值代表该孔隙体素所对应的孔隙半径；
控制部，与扫描重构部、二值化部、三维数组生成部以及空隙半径求解部均通信相连，控制它们的运行。
优选地，本发明提供的基于CT技术的孔隙尺寸计算装置，还可以包括：孔隙差异确定部，与控制部通信相连，依据空隙半径求解部得到的数组B中ip[i][j][k]将孔隙按半径大小进行分类区分、统计数量和确定分布情况。
优选地，本发明提供的基于CT技术的孔隙尺寸计算装置，还可以包括：试样三维模型生成部，与控制部通信相连，根据各体素与元素的对应关系和数组B，生成反映真实孔隙信息的试样孔隙尺寸及分布情况三维模型图。
优选地，本发明提供的基于CT技术的孔隙尺寸计算装置，还可以包括：输入显示部，与扫描重构部、二值化部、三维数组生成部、空隙半径求解部以及控制部均通信相连，用于让用户输入操作指令，并进行相应显示。
优选地，本发明提供的基于CT技术的孔隙尺寸计算装置，还可以包括：输入显示部能够根据操作指令对扫描重构部获取的试样原始投影图像、重构的试样三维图像和二维切片灰度图进行相应显示，对二值化部的处理数据进行显示，对三维数组生成部生成的三维数组 A进行显示，对空隙半径求解部求解结果进行显示，并能够将数组A和数组B以网格图的形式进行显示和比较。
优选地，本发明提供的基于CT技术的孔隙尺寸计算装置，还可以具有这样的特征：输入显示部还能够根据孔隙半径大小不同，在对应的试样孔隙尺寸及分布情况三维模型上用不同颜色进行不同半径级别的区分，并能够根据操作指令以整体的形式在三维模型中显示以不同颜色表示的各半径级别的孔隙或者按不同半径分级分别独立地在三维模型中仅显示相应半径级别的孔隙。
发明的作用与效果
本发明通过对三维数组形成的孔隙空间进行搜索判断得到与固体或试样边界相切且不含固体体素的最大内切球，从而求解得到孔隙空间中各孔隙体素的孔隙半径，确定孔隙结构，具有操作简单可行成本低，适用面广的优点，能够在不损伤试样、不简化和改变试样结构形状的前提下，形象直观的表征试样中孔隙的真实形貌尺寸，而且能够精确地计算孔隙尺寸并分类，对于土壤，混凝土，水泥等材料的内部孔隙结构和力学性能研究具有重要意义。
附图说明
图1为本发明实施例中重建获得的试样三维图像；
图2为本发明实施例中扫描得到的试样某一处区域空隙的切片图(a)和对应的空隙半径求解结果图(b)；(b)中灰色像素代表固体体素，白色像素代表孔隙体素，数字代表该点孔隙半径；
图3为本发明实施例涉及的空隙半径求解的流程图；
图4为本发明实施例中得到的不同尺寸孔隙的数量占比图；
图5为本发明实施例中得到的含所有级别空隙的孔隙尺寸分布三维模型图；
图6为本发明实施例中得到的仅显示半径为1～5体素级别孔隙的孔隙尺寸分布三维模型图；
图7为本发明实施例中得到的仅显示半径为6～10体素级别孔隙的孔隙分布三维模型图；
图8为本发明实施例中得到的仅显示半径为11～90体素级别孔隙的孔隙分布三维模型图。
具体实施方式
以下结合附图对本发明涉及的基于CT技术的孔隙尺寸计算方法及装置进行详细地说明。
<实施例>
本实施例所提供的基于CT技术的孔隙尺寸计算方法包括以下步骤：
步骤1.选取一块多孔介质试样，并使用三维成像显微***对试样进行CT扫描，获取原始投影图像；将原始图像进行图像处理和三维重构，得到如图1所示由1000×1000×1000 个体素构成的试样三维图像。
步骤2.将试样三维图像沿x方向按顺序以二维切片的图片保存，得到试样的一系列二维切片灰度图，每个切片图上包含1000×1000个体素。将切片灰度图进行二值化处理，得到仅包含试样固体体素和孔隙体素的二值化图片，将固体体素用0表示，孔隙体素用1表示，则每张切片图都可以用一个元素由0和1构成的二维数组表示，每个二维数组包含1000× 1000个元素。
步骤3.将所有二值化后的切片图进行叠加，得到试样的三维二值图像，即为(1000× 1000×1000)的三维数组，记为A。将三维数组中每个体素对应的数字赋值给ip[i][j][k]，其中i，j，k分别代表该体素在x，y，z方向的坐标，ip[i][j][k]＝1时，该体素为孔隙，统计得到孔隙数量为80523805；ip[i][j][k]＝0时，该体素为固体，统计得到固体点数量为1000000000，根据上文公式1计算得到孔隙度为0.0805238。土壤孔隙度仅反映土壤孔隙“量”的间题，并不能说明土壤孔隙“质”的差别。即使两种土壤的孔隙度相同，如果大小孔隙的数量分配各异，土壤性质亦会有很大差异。
步骤4.按从左至右、由下至上、由后至前的顺序，对孔隙空间中每一个孔隙体素，即所有ip[i][j][k]＝1的体素，依次求解孔隙半径，具体包括以下子步骤：
(1)以孔隙体素(i，j，k)为中心，1个体素尺寸为半径作球(d0＝1)，设定体素(i，j，k)的半径为d0，并赋值数组A中对应的元素ip[i][j][k]＝d0。
(2)在体素(i，j，k)继续以d1＝d0+1个体素为半径作球，若所作球超出试样范围或含有固体体素(存在ip[i][j][k]＝0)，说明d0是以体素(i，j，k)为球心、与试样边界或固体边界内切球的最大半径，则球心体素(i，j，k)对应的孔隙半径为d0，改数组A中体素对应的元素ip[i][j][k]＝d0，并存储至数组B。
(3)若(2)中以d1为半径所作球未超出试样范围并且不含有固体体素，则设定d1为球心体素(i，j，k)的孔隙半径，并更新数组A中该点对应的元素ip[i][j][k]＝d1；然后以体素 (i，j，k)为球心、d2＝d1+1个体素尺寸为半径继续作球，若所作球超出试样范围或含有固体体素(存在ip[i][j][k]＝0)，说明d1是以体素(i，j，k)为球心、与试样边界或固体边界内切球的最大半径，则球心体素(i，j，k)对应的孔隙半径为d1，改数组A中体素对应的元素ip[i][j][k]＝d1，并存储至数组B。
(4)若(3)中以d2为半径所作球未超出试样范围并且不含有固体体素，则设定d2为球心体素(i，j，k)的孔隙半径，改数组A中该点对应的元素ip[i][j][k]＝d2，并存储至数组B；然后以体素(i，j，k)为球心、d3＝d2+1个体素尺寸为半径继续作球，若所作球超出试样范围或含有固体体素(存在ip[i][j][k]＝0)，说明d2是以体素(i，j，k)为球心、与试样边界或固体边界内切球的最大半径，则球心体素(i，j，k)对应的孔隙半径为d2，改数组A中体素对应的元素ip[i][j][k]＝d2，并存储至数组B。
(5)以此类推……直至求出体素(i，j，k)对应的孔隙半径；
(6)求出体素(i，j，k)对应的孔隙半径后，跳转至下一个孔隙体素重复以上步骤，直至求出整个试样所有孔隙体素对应的孔隙半径。
求解完整个试样后得到一个新的三维数组B，数组B中ip[i][j][k]＝0时表示该点为固体体素(与素组A中一样)；ip[i][j][k]不等于0时，其值代表该体素所对应的孔隙半径。
步骤5.为反映土壤孔隙性质的差异，依据步骤4得到的数组B中ip[i][j][k]将孔隙按半径大小进行分类区分、统计数量和确定分布情况。依据ip[i][j][k]将半径为1～5，6～10，11～90 个体素的孔隙进行分类统计得到孔隙尺寸的数量占比如下表1，根据上文公式2计算得到平均孔隙尺寸为17.6985体素。
表1.孔隙尺寸数量
| 孔隙半径(单位：体素) | 1～5 | 6～10 | 11～90 |
| 孔隙数量 | 11038170 | 2195767 | 67289868 |
根据求解结果得到孔隙尺寸占比如图4所示。
步骤6.根据各体素与元素的对应关系和求解结果生成孔隙尺寸分布图，生成反映真实孔隙信息的试样孔隙尺寸及分布情况三维模型图，如图5所示。进一步，将生成的三维模型图与试样三维图对比如图5和1所示：图中紫色(仅显示相应灰度)越深，代表孔隙半径越小，绿色(仅显示相应灰度)越深代表孔隙半径越大；可以看到计算得到的图5所示的孔隙尺寸分布与图1所示扫描得到的试样表面的孔隙分布一致。
步骤7.根据孔隙半径大小不同，在对应的试样孔隙尺寸及分布情况三维模型上用不同颜色进行不同半径级别的区分，并能够整体或者分级分别显示不同半径级别的孔隙。筛选出孔隙半径在1到5个像素点之间的孔隙并显示如图6所示；筛选出孔隙半径在6到10个像素点之间的孔隙并显示如图7所示；筛选出孔隙半径在11到90个像素点之间的孔隙并显示如图8所示；根据孔隙分布图可以看出，大孔径的孔隙主要分布在试样表面且大部分相互连通，小孔径的孔隙均匀密布整个试样。
进一步，本实施例还提供能够自动实现上述方法的基于CT技术的孔隙尺寸计算装置，该装置包括扫描重构部、二值化部、三维数组生成部、空隙半径求解部、孔隙差异确定部、试样三维模型生成部、输入显示部以及控制部。
扫描重构部使用三维成像显微***对待计算的多孔介质试样进行CT扫描，获取试样原始投影图像；将试样原始图像进行图像处理和三维重构，得到由nx×ny×nz个体素构成的试样三维图像，体素表示三维图像中的一个三维像素点；将图像沿x方向按顺序以二维切片的图片保存，得到试样的一系列二维切片灰度图，每个切片图上包含ny×nz个体素。
二值化部将切片灰度图进行二值化处理，得到仅包含试样固体体素和孔隙体素的二值化图片，将固体体素用数值a0表示，孔隙体素用数值a1表示，则每张切片图都可以用由a0和a1构成的二维数组表示，每个二维数组包含ny×nz个数值元素，分别与ny×nz个体素一一对应，元素与其对应的体素具有相同的坐标。
三维数组生成部将所有二值化后的切片图进行叠加，得到试样的三维二值图像，该三维二值图像为(nx×ny×nz)的三维数组，记为数组A；数组A中的元素记为ip[i][j][k]，其中i，j，k分别代表元素所对应的体素在x，y，z方向的坐标，ip[i][j][k]＝a1时，对应的体素为孔隙；ip[i][j][k]＝a0时，对应的体素为固体。
空隙半径求解部对孔隙空间中每一个ip[i][j][k]＝a1的孔隙体素，依次求解孔隙半径，具体包括以下子步骤：
步骤4-1.将当前待求解的孔隙体素的半径d初始值设定为1个体素尺寸，d＝1，以当前待求解的孔隙体素为中心，d为半径作球；
步骤4-2.对球体区域中的元素进行判断：
若球体区域中存在ip[i][j][k]＝a0的元素，表明所作球超出试样范围或含有固体体素，说明d是以当前待求解的孔隙体素为球心、与试样边界或固体边界内切球的最大半径，则球心体素对应的孔隙半径为d，令元素ip[i][j][k]＝d；然后，进入步骤4-3；
若球体区域中不存在ip[i][j][k]＝a0的元素，表明以d为半径所作球未超出试样范围并且不含有固体体素，进入步骤4-4；
步骤4-3.判断当前已求解的孔隙体素是否为最后一个待求解的孔隙体素，若是，则求解结束，若否，则将下一个孔隙体素设定为当前待求解的体素，返回步骤4-1继续求解；
步骤4-4.令d＝d+1，将球的半径扩大1个体素尺寸，然后返回步骤4-2；
通过以上步骤4-1～步骤4-2得到显示试样所有孔隙体素的孔隙半径信息的新的三维数组B，该数组B中ip[i][j][k]＝a0时表示该点为固体体素；ip[i][j][k]不等于a0时，其值代表该孔隙体素所对应的孔隙半径。
孔隙差异确定部依据空隙半径求解部得到的数组B中ip[i][j][k]将孔隙按半径大小进行分类区分、统计数量和确定分布情况，如表1和图4所示。
试样三维模型生成部根据各体素与元素的对应关系和数组B，生成反映真实孔隙信息的试样孔隙尺寸及分布情况三维模型图，如图5和6～8所示。
输入显示部用于让用户输入操作指令，并进行相应显示。例如，输入显示部能够根据操作指令对扫描重构部获取的试样原始投影图像、重构的试样三维图像(图1所示)和二维切片灰度图进行相应显示，对二值化部的处理数据进行显示，对三维数组生成部生成的三维数组A进行显示，对空隙半径求解部求解结果进行显示(图2(b)所示)，并能够将数组A和数组B以网格图的形式进行显示和比较(图2所示)。输入显示部还能够根据孔隙半径大小不同，在对应的试样孔隙尺寸及分布情况三维模型图上用不同颜色进行不同半径级别的区分(图5所示)，并能够根据操作指令以整体的形式在三维模型中显示以不同颜色表示的各半径级别的孔隙或者按不同半径分级分别独立地在三维模型中仅显示相应半径级别的孔隙 (图6～8所示)。
控制部与扫描重构部、二值化部、三维数组生成部、空隙半径求解部孔隙差异确定部、试样三维模型生成部、输入显示部均通信相连，控制它们的运行。
以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的基于CT技术的孔隙尺寸计算方法及装置并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。
Claims (10)
1.基于CT技术的孔隙尺寸计算方法，其特征在于，包括以下步骤：
步骤1.使用三维成像显微***对待计算的多孔介质试样进行CT扫描，获取试样原始投影图像；将试样原始图像进行图像处理和三维重构，得到由nx×ny×nz个体素构成的试样三维图像，体素表示三维图像中的一个三维像素点；将图像沿x方向按顺序以二维切片的图片保存，得到试样的一系列二维切片灰度图，每个切片图上包含ny×nz个体素；
步骤2.将切片灰度图进行二值化处理，得到仅包含试样固体体素和孔隙体素的二值化图片，将固体体素用数值a0表示，孔隙体素用数值a1表示，则每张切片图都可以用由a0和a1构成的二维数组表示，每个二维数组包含ny×nz个数值元素，分别与ny×nz个体素一一对应；
步骤3.将所有二值化后的切片图进行叠加，得到试样的三维二值图像，该三维二值图像为(nx×ny×nz)的三维数组，记为数组A；数组A中的元素记为ip[i][j][k]，其中i，j，k分别代表元素所对应的体素在x，y，z方向的坐标，ip[i][j][k]＝a1时，对应的体素为孔隙；ip[i][j][k]＝a0时，对应的体素为固体；
步骤4.对孔隙空间中每一个ip[i][j][k]＝a1的孔隙体素，依次求解孔隙半径，具体包括以下子步骤：
步骤4-1.将当前待求解的孔隙体素的半径d初始值设定为1个体素尺寸，d＝1，以当前待求解的孔隙体素为中心，d为半径作球；
步骤4-2.对球体区域中的元素进行判断：
若球体区域中存在ip[i][j][k]＝a0的元素，表明所作球超出试样范围或含有固体体素，说明d是以当前待求解的孔隙体素为球心、与试样边界或固体边界内切球的最大半径，则球心体素对应的孔隙半径为d，令元素ip[i][j][k]＝d；然后，进入步骤4-3；
若球体区域中不存在ip[i][j][k]＝a0的元素，表明以d为半径所作球未超出试样范围并且不含有固体体素，进入步骤4-4；
步骤4-3.判断当前已求解的孔隙体素是否为最后一个待求解的孔隙体素，若是，则求解结束，若否，则将下一个孔隙体素设定为当前待求解的体素，返回步骤4-1继续求解；
步骤4-4.令d＝d+1，将球的半径扩大1个体素尺寸，然后返回步骤4-2；
通过以上步骤4-1～步骤4-2得到显示试样所有孔隙体素的孔隙半径信息的新的三维数组B，该数组B中ip[i][j][k]＝a0时表示该点为固体体素；ip[i][j][k]不等于a0时，其值代表该孔隙体素所对应的孔隙半径。
2.根据权利要求1所述的基于CT技术的孔隙尺寸计算方法，其特征在于，还包括：
步骤5.为反映土壤孔隙性质的差异，依据步骤4得到的数组B中ip[i][j][k]将孔隙按半径大小进行分类区分、统计数量和确定分布情况。
3.根据权利要求1所述的基于CT技术的孔隙尺寸计算方法，其特征在于，还包括：
步骤6.根据各体素与元素的对应关系和步骤4得到的数组B，生成反映真实孔隙信息的试样孔隙尺寸及分布情况三维模型图。
4.根据权利要求1所述的基于CT技术的孔隙尺寸计算方法，其特征在于，还包括：
步骤7.根据孔隙半径大小不同，在对应的试样孔隙尺寸及分布情况三维模型上用不同颜色进行不同半径级别的区分，并能够整体或者分级分别显示不同半径级别的孔隙。
5.基于CT技术的孔隙尺寸计算装置，其特征在于，包括：
扫描重构部，使用三维成像显微***对待计算的多孔介质试样进行CT扫描，获取试样原始投影图像；将试样原始图像进行图像处理和三维重构，得到由nx×ny×nz个体素构成的试样三维图像，体素表示三维图像中的一个三维像素点；将图像沿x方向按顺序以二维切片的图片保存，得到试样的一系列二维切片灰度图，每个切片图上包含ny×nz个体素；
二值化部，将切片灰度图进行二值化处理，得到仅包含试样固体体素和孔隙体素的二值化图片，将固体体素用数值a0表示，孔隙体素用数值a1表示，则每张切片图都可以用由a0和a1构成的二维数组表示，每个二维数组包含ny×nz个数值元素，分别与ny×nz个体素一一对应；
三维数组生成部，将所有二值化后的切片图进行叠加，得到试样的三维二值图像，该三维二值图像为(nx×ny×nz)的三维数组，记为数组A；数组A中的元素记为ip[i][j][k]，其中i，j，k分别代表元素所对应的体素在x，y，z方向的坐标，ip[i][j][k]＝a1时，对应的体素为孔隙；ip[i][j][k]＝a0时，对应的体素为固体；
空隙半径求解部，对孔隙空间中每一个ip[i][j][k]＝a1的孔隙体素，依次求解孔隙半径，具体包括以下子步骤：
步骤4-1.将当前待求解的孔隙体素的半径d初始值设定为1个体素尺寸，d＝1，以当前待求解的孔隙体素为中心，d为半径作球；
步骤4-2.对球体区域中的元素进行判断：
若球体区域中存在ip[i][j][k]＝a0的元素，表明所作球超出试样范围或含有固体体素，说明d是以当前待求解的孔隙体素为球心、与试样边界或固体边界内切球的最大半径，则球心体素对应的孔隙半径为d，令元素ip[i][j][k]＝d；然后，进入步骤4-3；
若球体区域中不存在ip[i][j][k]＝a0的元素，表明以d为半径所作球未超出试样范围并且不含有固体体素，进入步骤4-4；
步骤4-3.判断当前已求解的孔隙体素是否为最后一个待求解的孔隙体素，若是，则求解结束，若否，则将下一个孔隙体素设定为当前待求解的体素，返回步骤4-1继续求解；
步骤4-4.令d＝d+1，将球的半径扩大1个体素尺寸，然后返回步骤4-2；
通过以上步骤4-1～步骤4-2得到显示试样所有孔隙体素的孔隙半径信息的新的三维数组B，该数组B中ip[i][j][k]＝a0时表示该点为固体体素；ip[i][j][k]不等于a0时，其值代表该孔隙体素所对应的孔隙半径；
控制部，与所述扫描重构部、所述二值化部、所述三维数组生成部以及所述空隙半径求解部均通信相连，控制它们的运行。
6.根据权利要求5所述的基于CT技术的孔隙尺寸计算装置，其特征在于，还包括：
孔隙差异确定部，与所述控制部通信相连，依据所述空隙半径求解部得到的数组B中ip[i][j][k]将孔隙按半径大小进行分类区分、统计数量和确定分布情况。
7.根据权利要求5所述的基于CT技术的孔隙尺寸计算装置，其特征在于，还包括：
试样三维模型生成部，与所述控制部通信相连，根据各体素与元素的对应关系和数组B，生成反映真实孔隙信息的试样孔隙尺寸及分布情况三维模型图。
8.根据权利要求5所述的基于CT技术的孔隙尺寸计算装置，其特征在于，还包括：
输入显示部，与所述扫描重构部、所述二值化部、所述三维数组生成部、所述空隙半径求解部以及所述控制部均通信相连，用于让用户输入操作指令，并进行相应显示。
9.根据权利要求8所述的基于CT技术的孔隙尺寸计算装置，其特征在于：
其中，所述输入显示部能够根据操作指令对所述扫描重构部获取的所述试样原始投影图像、重构的试样三维图像和二维切片灰度图进行相应显示，对所述二值化部的处理数据进行显示，对所述三维数组生成部生成的三维数组A进行显示，对所述空隙半径求解部求解结果进行显示，并能够将数组A和数组B以网格图的形式进行显示和比较。
10.根据权利要求8所述的基于CT技术的孔隙尺寸计算装置，其特征在于：
其中，所述输入显示部还能够根据孔隙半径大小不同，在对应的试样孔隙尺寸及分布情况三维模型图上用不同颜色进行不同半径级别的区分，并能够根据操作指令以整体的形式在三维模型中显示以不同颜色表示的各半径级别的孔隙或者按不同半径分级分别独立地在三维模型中仅显示相应半径级别的孔隙。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210047008.8A CN114565658B (zh) | 2022-01-14 | 基于ct技术的孔隙尺寸计算方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210047008.8A CN114565658B (zh) | 2022-01-14 | 基于ct技术的孔隙尺寸计算方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114565658A true CN114565658A (zh) | 2022-05-31 |
| CN114565658B CN114565658B (zh) | 2024-07-02 |
Family
ID=
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115254213A (zh) * | 2022-06-24 | 2022-11-01 | 中国计量大学 | 一种基于真实土壤孔隙网络的微流控芯片装置 |
| CN115455794A (zh) * | 2022-08-10 | 2022-12-09 | 武汉理工大学 | 基于连通孔隙划分计算区域的lbm并行优化方法、装置及存储介质 |
| CN116337708A (zh) * | 2022-12-06 | 2023-06-27 | 武汉理工大学 | 提取多孔材料孔隙特征和流体分布的方法与*** |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070070068A1 (en) * | 2005-09-27 | 2007-03-29 | Matthias John | Method and apparatus for processing of diagnostic image data |
| CN101556703A (zh) * | 2009-05-16 | 2009-10-14 | 中国石油大学(华东) | 基于连续切片图像的网络模型建立方法 |
| JP2016146132A (ja) * | 2015-02-09 | 2016-08-12 | 学校法人日本大学 | 形状特徴抽出方法、形状特徴抽出処理装置、形状記述方法及び形状分類方法 |
| US20180121579A1 (en) * | 2016-10-31 | 2018-05-03 | Bp Corporation North America Inc. | Direct numerical simulation of petrophysical properties of rocks with two or more immicible phases |
| CN110793898A (zh) * | 2019-10-22 | 2020-02-14 | 浙江大学 | 一种定量分析土柱中不同大小3d孔隙空间分布的方法 |
| CN111968089A (zh) * | 2020-08-15 | 2020-11-20 | 晋江市博感电子科技有限公司 | 一种基于最大内切球机制的l1中值骨架提取方法 |
| CN113129275A (zh) * | 2021-03-31 | 2021-07-16 | 中国矿业大学 | 一种基于岩土体材料数字图像三维结构表征方法 |
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070070068A1 (en) * | 2005-09-27 | 2007-03-29 | Matthias John | Method and apparatus for processing of diagnostic image data |
| CN101556703A (zh) * | 2009-05-16 | 2009-10-14 | 中国石油大学(华东) | 基于连续切片图像的网络模型建立方法 |
| JP2016146132A (ja) * | 2015-02-09 | 2016-08-12 | 学校法人日本大学 | 形状特徴抽出方法、形状特徴抽出処理装置、形状記述方法及び形状分類方法 |
| US20180121579A1 (en) * | 2016-10-31 | 2018-05-03 | Bp Corporation North America Inc. | Direct numerical simulation of petrophysical properties of rocks with two or more immicible phases |
| CN109891460A (zh) * | 2016-10-31 | 2019-06-14 | Bp北美公司 | 具有两个或更多个不混溶相的岩石的岩石物理性质的直接数值模拟 |
| CN110793898A (zh) * | 2019-10-22 | 2020-02-14 | 浙江大学 | 一种定量分析土柱中不同大小3d孔隙空间分布的方法 |
| CN111968089A (zh) * | 2020-08-15 | 2020-11-20 | 晋江市博感电子科技有限公司 | 一种基于最大内切球机制的l1中值骨架提取方法 |
| CN113129275A (zh) * | 2021-03-31 | 2021-07-16 | 中国矿业大学 | 一种基于岩土体材料数字图像三维结构表征方法 |
Non-Patent Citations (4)
| Title |
|---|
| SHOICHIRO HAMAMOTO等: "Pore network structure linked by X-ray CT to particle characteristics and transport parameters", 《SOILS AND FOUNDATIONS》 * |
| 杨文丽: "序列断层粗粒土图像三维重建及可视化技术研究", 《中国博士学位论文全文数据库 信息科技辑》 * |
| 林缅;江文滨;李勇;易智星;张召彬;: "页岩油(气)微尺度流动中的若干问题", 矿物岩石地球化学通报, no. 01, 10 January 2015 (2015-01-10) * |
| 胡沛裕;王树刚;宋双林;蒋爽;梁运涛;: "基于完全并行细化算法的孔隙网络中轴提取方法", 化工学报, no. 1 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115254213A (zh) * | 2022-06-24 | 2022-11-01 | 中国计量大学 | 一种基于真实土壤孔隙网络的微流控芯片装置 |
| CN115254213B (zh) * | 2022-06-24 | 2024-05-03 | 中国计量大学 | 一种基于真实土壤孔隙网络的微流控芯片装置 |
| CN115455794A (zh) * | 2022-08-10 | 2022-12-09 | 武汉理工大学 | 基于连通孔隙划分计算区域的lbm并行优化方法、装置及存储介质 |
| CN115455794B (zh) * | 2022-08-10 | 2024-03-29 | 武汉理工大学 | 基于连通孔隙划分计算区域的lbm并行优化方法、装置及存储介质 |
| CN116337708A (zh) * | 2022-12-06 | 2023-06-27 | 武汉理工大学 | 提取多孔材料孔隙特征和流体分布的方法与*** |
| CN116337708B (zh) * | 2022-12-06 | 2023-11-17 | 武汉理工大学 | 提取多孔材料孔隙特征和流体分布的方法与*** |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109697752B (zh) | 基于岩心ct图像孔隙信息提取定量表征岩心非均质性方法 | |
| CN107194872B (zh) | 基于内容感知深度学习网络的遥感图像超分辨率重建方法 | |
| CN110135311A (zh) | 一种基于三维岩心扫描图像的孔隙与孔喉识别***及方法 | |
| GB2398469A (en) | Identifying an object from images recorded at different positions and orientations | |
| WO2005112769A1 (en) | Nodule detection | |
| CN108205806A (zh) | 一种锥束ct图像三维颅面结构的自动解析方法 | |
| CN110728666B (zh) | 基于数字病理玻片进行慢性鼻窦炎的分型方法及其*** | |
| CN109345625A (zh) | 岩心图像自适应分区三维重建方法 | |
| CN111598803B (zh) | 一种基于变分辨率体素格网与稀疏卷积的点云滤波方法 | |
| CN113706603A (zh) | 页岩有机质孔隙连通性分类表征的方法 | |
| Landini | Applications of fractal geometry in pathology | |
| CN110264472B (zh) | 一种基于sd-oct的水果近皮下细胞无损成像方法 | |
| CN113160392A (zh) | 基于深度神经网络的光学建筑目标三维重建方法 | |
| CN112651955A (zh) | 一种肠道图像的识别方法及终端设备 | |
| CN103903255B (zh) | 一种超声图像分割方法和*** | |
| CN114511546A (zh) | 基于dbscan聚类和四象限的激光点云林木胸径获取方法 | |
| CN112071430B (zh) | 一种病理指标的智能预测*** | |
| CN117611542A (zh) | 一种基于胎儿宫内颅脑影像检测方法及*** | |
| Brabant et al. | EDART, a discrete algebraic reconstructing technique for experimental data obtained with high resolution computed tomography | |
| CN114565658B (zh) | 基于ct技术的孔隙尺寸计算方法及装置 | |
| CN114565658A (zh) | 基于ct技术的孔隙尺寸计算方法及装置 | |
| CN109829879B (zh) | 维管束的检测方法及装置 | |
| CN112002019B (zh) | 一种基于mr混合现实的模拟人物阴影的方法 | |
| CN110739051B (zh) | 利用鼻息肉病理图片建立嗜酸性粒细胞占比模型的方法 | |
| Bolshakov et al. | Investigation of the pore space structure by a scanning electron microscope using the computer program collector |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |