A three-dimensional modeling method of reef limestone microstructure
By normalizing gray values, calculating global thresholds, and segmenting using Gaussian mixture models, combined with morphological operations and graph cut algorithms, the problem of insufficient accuracy and stability in reef limestone processing using traditional digital core technology has been solved, achieving efficient and robust 3D digital core modeling.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- SANYA SCI & EDUCATION INNOVATION PARK WUHAN UNIV OF TECH
- Filing Date
- 2026-03-26
- Publication Date
- 2026-07-03
Smart Images
Figure CN121904291B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of three-dimensional modeling technology for reef limestone, specifically a three-dimensional modeling method for the microstructure of reef limestone. Background Technology
[0002] With the development of marine underground spaces, complex reef limestone strata are frequently encountered in the design and construction of many geotechnical engineering projects. Reef limestone is a special type of island and reef rock mass. The unique composition and development environment result in numerous pores, varied structures, and significant differences in strength, making it quite different from terrestrial rocks. To ensure efficient and safe development of marine underground spaces, research on the physical and mechanical properties of reef limestone is necessary to provide effective guidance and reference for engineering construction.
[0003] Digital core technology (DCore) is based on 3D grayscale volumes acquired by industrial / micron-level CT scans. Through image processing, volumetric data modeling, and parameterization, it transforms the pore-matrix structure inside the actual core into a computable digital model, which is then used for seepage and multiphysics simulations. Reef limestone, as an important type of carbonate rock, exhibits strong heterogeneity and dual-porosity characteristics. It has a large pore-throat scale span, significant grayscale overlap, prominent ring artifacts and partial voxel effects, and easily develops boundary regions connected to the outside air during sample preparation and scanning. Due to these characteristics, applying traditional DCore techniques to reef limestone presents a dual challenge of accuracy and stability: on the one hand, it is necessary to suppress missegmentation while maintaining pore-throat connectivity and geometric details; on the other hand, traditional DCore techniques rely primarily on fixed or empirical thresholds, lacking adaptability and robustness to statistical fluctuations across slices. Summary of the Invention
[0004] The purpose of this invention is to provide a three-dimensional modeling method for the microstructure of reef limestone in order to address the shortcomings of existing technologies. While maintaining the topological connectivity of pore throats and detail resolution, the method outputs digital cores and associated parameters that can be directly used for seepage and multiphysics simulations through statistically robust strength calibration, sample shell mask construction, reliable seed generation, and constrained graph model solving, thereby solving the problems existing in the above-mentioned existing technologies.
[0005] To achieve the above objectives, the present invention adopts the following technical solution:
[0006] A three-dimensional modeling method for the microstructure of reef limestone includes the following steps: S1, scanning and acquiring the original grayscale image of the core sample, and preprocessing the original grayscale image; S2, extracting the sample shell mask based on the preprocessed grayscale image; S3, establishing a Gaussian mixture model of the sample grayscale distribution and determining the segmentation threshold; S4, classifying and segmenting the volume data according to the segmentation threshold to obtain pore candidate regions, matrix candidate regions, and unresolved regions; S5, performing morphological closing operations on the pore candidate regions and matrix candidate regions to obtain the set of pixels belonging to pores and the set of pixels belonging to pores. S6. For unresolved regions, a graph-cut-based energy minimization algorithm is used for classification and segmentation. In S2, the process of extracting the sample shell mask is as follows: S21. Calculate and select a global threshold to identify air candidate regions; S22. Perform morphological opening operation on the air candidate regions to obtain the mask image; S23. Perform morphological reconstruction using the image boundary as a seed to obtain the initial shell mask; S24. Perform morphological closing operation on the initial shell mask; S25. Retain the connected region with the largest number of voxels as the final sample shell mask.
[0007] Furthermore, in S1, the preprocessing includes grayscale value normalization. The grayscale value normalization process is as follows: The grayscale value of the low quantile is set to... The gray value of the high quantile is grayscale value for each voxel The linear stretching normalization is calculated using the following formula:
[0008]
[0009] in, For normalized volume data, Represent the position of each voxel; map the normalized volume data to... The grayscale value outside the specified range is truncated.
[0010] Furthermore, in S21, the Otsu method is applied to automatically calculate the global threshold, the process of which is as follows: For any CT slice, the grayscale value is... Quantified as The grayscale values were analyzed at several levels, and the frequency of each grayscale level was counted to obtain a normalized histogram. It satisfies the following formula:
[0011]
[0012] For a certain candidate threshold Based on this gray level, pixels are divided into two categories: low gray level category 1 Grayscale High grayscale category Grayscale The intra-class weights, intra-class average gray levels, and global average gray levels of low-gray-level pixels and high-gray-level pixels are calculated separately, as follows:
[0013] Low grayscale Intra-class weight The calculation formula is as follows:
[0014]
[0015] High grayscale Intra-class weight The calculation formula is as follows:
[0016]
[0017] Low grayscale Intra-class average gray level The calculation formula is as follows:
[0018]
[0019] High grayscale Intra-class average gray level The calculation formula is as follows:
[0020]
[0021] Global average gray level The calculation formula is as follows:
[0022]
[0023] Calculate the candidate threshold Inter-class variance The calculation formula is as follows:
[0024]
[0025] For possible candidate thresholds Calculate the inter-class variance one by one Selecting classes with different variances Maximum candidate threshold As a global threshold The calculation formula is as follows:
[0026]
[0027] In step S22, the morphological opening operation is performed as follows: Let the air candidate region image be... The neighborhood range is The formula for calculating the opening operation is as follows:
[0028]
[0029] in, For mask images;
[0030] In step S23, the morphological reconstruction process is as follows: Let the image boundary seed be the labeled image. Select structural elements Then there is an expansion operator. The calculation formula is as follows:
[0031]
[0032] Let the initial state Then there is an iteration to the th When walking, When the iteration reaches a certain point that satisfies... When the iteration stops, the final result of the iteration is the external air mask. External air mask Take the complement to obtain the initial shell mask. ;
[0033] In S24, the morphological closing operation is performed as follows: Selecting a structuring element. Then the formula for calculating the closing operation is as follows:
[0034]
[0035] in, This is a binary mask of the shell after morphological closing operations.
[0036] In S25, the process for determining the connected region with the largest number of voxels is as follows: The shell binary mask after morphological closing operation processing... Perform 3D connected component labeling, assigning a connected component label to each foreground voxel as follows:
[0037]
[0038] in, The total number of connected regions; for each label The number of voxels in the corresponding connected regions is counted using the following formula:
[0039]
[0040] in, For the first The size of the connected region corresponding to each label; making The largest tag The corresponding region is the connected region with the largest number of voxels.
[0041] Further, in S3, the process of establishing a Gaussian mixture model of sample grayscale distribution and determining the segmentation threshold is as follows: S31, extract slices at equal intervals from the volume data; S32, randomly sample the grayscale values of the shell mask region in each slice; S33, summarize all sampling results to form an overall grayscale training sample set; S34, establish a Gaussian mixture model; S35, perform maximum likelihood fitting on the grayscale training sample set to obtain model parameters and plot the posterior probability curve; S36, determine the segmentation threshold based on the posterior probability curve.
[0042] Furthermore, in S31, the process of extracting slices is as follows: Let the total number of slices be... The number of representative slices required for training and statistics is Calculate the equal interval step size The calculation formula is as follows:
[0043]
[0044] The slice number is automatically determined based on equal intervals, calculated using the following formula:
[0045]
[0046] in, For the first The serial number of the representative slice, ;
[0047] In step S32, the random sampling process is as follows: setting the sampling ratio. Calculate the first Number of sampling points per slice The calculation formula is as follows:
[0048]
[0049] in, The total number of pixels within the shell mask area; Pixels Number it, and then Randomly generated Let there be a set of unique integer indices, and denote the grayscale values corresponding to these indices as... , as the grayscale sample of the slice;
[0050] In S33, the unified calculation formula is as follows:
[0051]
[0052] in, This is the overall grayscale training sample set; For the merged first Individual sample grayscale values; grayscale training sample set Total number of samples The sum of the number of sampling points for each slice is calculated using the following formula:
[0053]
[0054] In S34, the Gaussian mixture model is established as follows:
[0055]
[0056] in, These correspond to three components: pores, grayscale overlap area, and matrix, respectively. For the first The mixing weight of each Gaussian component; For the first The mean of the Gaussian components; For the first The variance of each Gaussian component;
[0057] In S35, the maximum likelihood fitting process is as follows:
[0058] The EM algorithm is used to process the grayscale training sample set. Maximum likelihood fitting, i.e., maximizing the log-likelihood function, is performed as follows:
[0059]
[0060] The iterative steps of the EM algorithm are as follows: Calculate the posterior probability: given the current parameters... In the case of calculating each sample Belongs to the The posterior probability of each Gaussian component is calculated using the following formula:
[0061]
[0062] Update model parameters: Re-estimate the effective sample size and parameters of each Gaussian component using the above posterior probabilities, as shown in the following formula:
[0063]
[0064] After completing one iteration, a new set of parameters is obtained, and then the next iteration continues. After each iteration, the overall log-likelihood value is calculated based on the current parameters. The iteration stops when the goodness of fit changes very little between two adjacent iterations. If the number of iterations reaches the set upper limit before the algorithm converges, the iteration is terminated early. The three sets of parameters obtained at the end of the iteration are... These are the final model parameters;
[0065] In S36, the process of determining the segmentation threshold is as follows: Select the posterior probability threshold according to the posterior probability curve. When a gray value falls into the lowest gray component and its posterior probability is not less than the posterior probability threshold, the gray value is defined as the porosity threshold. When a gray value falls into the highest gray component and its posterior probability is not less than the posterior probability threshold, the gray value is defined as the matrix threshold.
[0066] Furthermore, in S4, the classification and segmentation process is as follows:
[0067] For each slice, its pixel grayscale values are iterated; when the grayscale value is less than or equal to the pore threshold and the pixel is located within the shell mask, it is marked as a pore voxel; when the grayscale value is greater than or equal to the matrix threshold and the pixel is located within the shell mask, it is marked as a matrix voxel; when the pore threshold is less than or equal to the grayscale value less than or equal to the matrix threshold and the pixel is located within the shell mask, it is marked as a pending voxel; thus, the slice is classified into pore candidate regions, matrix candidate regions, and pending regions.
[0068] Further, in S6, the segmentation and classification process is as follows: S61, define the region of interest, and construct a corresponding node for each pixel within the region of interest; S62, add neighborhood edges for adjacent pixels and assign a smoothing weight to each neighborhood edge; S63, for each pixel within the region of interest, calculate its negative log-likelihood cost belonging to the pores and matrix, and use it as the connection weight for connecting the pixel to the background terminal and the foreground terminal, respectively; S64, establish an energy function based on the smoothing weight and connection weight, apply the maximum flow algorithm to extract the minimum cut, and perform optimal segmentation of the voxels in the unresolved region.
[0069] Furthermore, in S61, the process of establishing the region of interest is as follows:
[0070] Calculate the minimum bounding rectangle of the shell mask in the current slice, and extend the safety boundary to a predetermined width in all directions. The local rectangular region determined in this way is denoted as the region of interest (ROI).
[0071] In S62, the smoothing weight The calculation formula is as follows:
[0072]
[0073]
[0074] in, The difference in grayscale between adjacent pixels; Pixels and grayscale value; For smoothing coefficients; The standard deviation of grayscale values within the region of interest (ROI);
[0075] In S63, the formula for calculating the negative log-likelihood cost of pores is as follows:
[0076]
[0077] in, The cost of this pixel belonging to the aperture. This represents the likelihood probability that the pixel belongs to the corresponding category. The parameters are Gaussian distribution parameters for pores;
[0078] The formula for calculating the negative log-likelihood cost of the matrix is as follows:
[0079]
[0080] in, The cost of this pixel belonging to the matrix. These are the parameters of a matrix-type Gaussian distribution.
[0081] Furthermore, in S64, before performing optimal segmentation on the voxels in the unresolved region, the set of pore pixels is connected to the foreground terminal with an edge whose weight tends to infinity, forcing it to be connected to the background in the optimal segmentation so that it is always classified as pore; the set of matrix pixels is connected to the background terminal with an edge whose weight tends to infinity, forcing it to be connected to the foreground in the optimal segmentation so that it is always classified as matrix.
[0082] Compared with the prior art, the beneficial effects of the present invention are:
[0083] The proposed 3D modeling method for ultra-large CT core images can effectively process ultra-large-scale CT core image data. It employs GMM dual-threshold segmentation and shell mask construction technology to accurately extract core regions, avoiding the missegmentation of air at core boundaries as pores, thus completely preserving the true boundaries of the core. Furthermore, it generates reliable seed regions through closure operations and combines region of interest extraction and graph cut algorithms for precise pore identification. The automatic threshold segmentation results are stable and reliable, significantly reducing manual intervention. Compared with traditional core modeling techniques, this method has higher robustness and accuracy in identifying complex core pore structures. In addition, this invention optimizes resource usage for processing massive CT image data, employing parallel computing techniques such as memory mapping and block mapping to greatly reduce memory consumption. Even with ultra-large data, it can smoothly complete 3D reconstruction without being limited by computer memory and performance, significantly improving processing efficiency. Simultaneously, the 3D digital core model obtained by this invention is realistic, reliable, and structurally stable, and can be exported to a standardized 3D data format, facilitating integration with other software platforms. It possesses good format compatibility and engineering applicability, meeting the modeling needs of research involving coupling different scales and physical fields. Attached Figure Description
[0084] Figure 1 This is a schematic diagram of a CT slice after grayscale normalization in this invention; Figure 2 This is a schematic diagram of the Gaussian distribution fitting curve in this invention; Figure 3 This is a schematic diagram of a CT slice after dual-threshold segmentation in this invention; Figure 4 This is a schematic diagram of a CT slice after performing morphological closing operations in this invention; Figure 5 This is a schematic diagram of a CT slice after detailed segmentation in this invention. Detailed Implementation
[0085] This embodiment provides an efficient 3D modeling method for the microstructure of reef limestone using CT image data with extremely large voxels. Taking core CT scan grayscale data of approximately 2050×2050×1000 voxels as an example, memory mapping technology is used to map the original CT data of billions of voxels to the memory address space for parallel reading and processing. This enables the entire process of image segmentation and 3D model reconstruction to be completed in a normal computer environment. The main steps of this method are as follows:
[0086] S1. Continuous cross-sectional grayscale images of core samples are obtained using industrial CT scanning and stacked to form a three-dimensional raw data file. The obtained raw format data is imported into a computing platform for preprocessing to make the image brightness uniform, reduce noise, and enhance the contrast and resolution between pores and matrix.
[0087] Specifically, preprocessing includes grayscale value normalization, and the grayscale value normalization process is as follows:
[0088] Set the gray value of the low quantile to The gray value of the high quantile is For each voxel gray value Linear stretching normalization is performed according to the following formula:
[0089]
[0090] in, The data is the normalized volume data; This represents the position of each voxel; preferably, the low quantile is 1% and the high quantile is 99%, which effectively reduces the influence of extreme noise and artifacts, while better preserving the gray-level contrast between pores and the matrix, which is beneficial for subsequent threshold estimation and classification. Robust intensity normalization is achieved by calculating the low and high quantiles of gray values.
[0091] All voxel grayscale values are mapped to the [0,1] interval, and values outside the range are truncated. Specifically, to ensure the stability of the results, if the calculated value is less than 0, it is directly set to 0, and if it is greater than 1, it is directly set to 1, which is equivalent to "sandwiching" all values in the [0,1] interval.
[0092] After the above grayscale value normalization process, the normalized volume data is obtained. Its numerical range is stable within [0,1]. Subsequent steps such as shell mask extraction, threshold estimation, and classification are all uniformly based on this. As input. The CT slices after grayscale normalization are as follows. Figure 1 As shown, the grayscale contrast is improved, laying the foundation for subsequent segmentation. A quantile-driven histogram estimation strategy is adopted, using a pair of robust quantiles (e.g., the 1st and 99th percentiles) of the entire data as low / high reference values. The volume data is linearly compressed into the computational domain, effectively reducing the impact of extreme noise and ring artifacts on subsequent threshold estimation. Quantile estimation is achieved by accumulating histograms piece by piece (e.g., 4096 grayscale bins), which is a single-pass, constant-memory streaming computation, suitable for long Z-axis and high bit depth data.
[0093] Furthermore, after grayscale normalization, noise reduction filtering can be performed as needed to provide a more stable input for subsequent shell mask extraction and threshold estimation. Specifically, if the acquired CT images in actual applications have obvious noise or artifacts, any mature noise reduction algorithm (such as median filtering, bilateral filtering, mean filtering, etc.) can be selected for noise reduction filtering; if the image quality is already good, no additional noise reduction processing is required, and the process can proceed directly to the next step.
[0094] S2. In order to remove background air interference from outside the sample, a threshold-based and morphology reconstruction method is used to extract the shell mask of the sample.
[0095] S21. Apply the Otsu method to a representative CT slice after normalization to automatically calculate the global threshold and preliminarily identify obviously dark voxels as air candidate regions.
[0096] Specifically, the calculation process for the global threshold is as follows:
[0097] For any CT slice, the grayscale value will be... Quantified as The grayscale values were analyzed at several levels, and the frequency of each grayscale level was counted to obtain a normalized histogram. It satisfies the following formula:
[0098]
[0099] For a certain candidate threshold Using this gray level as the boundary, the pixels are divided into the following two categories:
[0100] Low grayscale Grayscale ;
[0101] High grayscale Grayscale .
[0102] For the two categories mentioned above, calculate the intra-class weight (proportion) and intra-class average gray level, and then calculate the global average gray level. The calculation process is as follows:
[0103] Low grayscale Intra-class weight The calculation formula is as follows:
[0104]
[0105] High grayscale Intra-class weight The calculation formula is as follows:
[0106]
[0107] Low grayscale Intra-class average gray level The calculation formula is as follows:
[0108]
[0109] High grayscale Intra-class average gray level The calculation formula is as follows:
[0110]
[0111] Global average gray level The calculation formula is as follows:
[0112]
[0113] Then calculate the candidate threshold. Inter-class variance The calculation formula is as follows:
[0114]
[0115] For all possible candidate thresholds Calculate the inter-class variance one by one Selecting classes with different variances Maximum candidate threshold As a global threshold The calculation formula is as follows:
[0116]
[0117] global threshold Often falling precisely near the "valley" between the two main peaks of the histogram, this allows for the automatic segmentation of images into two main categories of regions with distinct gray levels and appropriate proportions, without relying on human experience. By employing Otsu's method as the automatic global threshold selection method, this technique effectively distinguishes different gray-level regions and improves the stability of subsequent segmentation. For each CT slice, the global threshold is calculated once using Otsu's method to identify air candidate regions, eliminating the need for manual setting.
[0118] S22. Perform appropriate morphological opening operations on the candidate air regions to remove contiguous noise, making the air distribution purer, and obtain the mask image.
[0119] Specifically, the calculation process for the opening operation is as follows:
[0120] The "air candidate regions" obtained by threshold segmentation are represented as a binary image. Pixels within the air candidate region are recorded as 1, and those at other locations are recorded as 0. For ease of explanation, this embodiment uses a small structural element with a radius of approximately 1 pixel. (For example, a small circular or square template of 3×3) is used to describe the "neighborhood range".
[0121] Corrosion is defined as follows:
[0122]
[0123] The expansion is defined as follows:
[0124]
[0125] Based on this, the formula for calculating the opening operation is as follows:
[0126]
[0127] By using opening operations, very small, isolated noise and burrs can be removed without significantly altering the overall shape of the main area.
[0128] In this embodiment, the opening operation is performed only once for the air candidate regions of each slice, and the use of structuring elements with small radii is limited. This is known as "moderate opening operation". Moderate opening operation can effectively remove isolated noise and small artifacts in small areas, making the external air area cleaner, while not having a significant "over-smoothing" effect on the real sample contour, thus balancing noise reduction and preservation of geometric details.
[0129] S23. Perform morphological reconstruction (i.e. flooding filling starting from the edge) using the image boundary as a seed, retain the air regions that are only connected to the slice boundary, and then subtract these external air regions from the whole image to obtain the initial shell mask.
[0130] Specifically, the morphological reconstruction process is as follows:
[0131] Let the binary image of the air candidate after morphological opening be the mask image. In the mask image In the image, pixels whose grayscale features match air and are preserved after the opening operation are assigned a value of 1, while other pixels are assigned a value of 0; the image boundary seed is set to a marker image. In the labeled image In the image, pixels located on the four edges of the image and whose grayscale value can be determined to be air are assigned a value of 1, while all other pixels are assigned a value of 0.
[0132] Terrain Bloating (Single-Step Expansion): Select a 3×3 structuring element (e.g., a square template with eight neighboring regions) for the current labeled image Perform a standard expansion, that is This causes it to expand outwards in space, and then overlaps with the mask image. By finding the intersection, we obtain the dilation operator. dilation operator The expression is as follows:
[0133]
[0134] Terrain expansion ensures that the expansion process only occurs within... It is conducted within permitted areas and will not enter areas that cannot belong to air, thus ensuring the expansion of connectivity.
[0135] Iterative reconstruction: Let the initial state Then there is the first When walking, That is, in each iteration, the mask image is... Under the constraint, the currently known air region is expanded outwards by one layer. When the iteration reaches a certain point that satisfies... When the iteration stops, no more connected components can be obtained. Record the final result. ,Right now , "Being located in the air candidate region" Inside, and at the outer boundary, seeds The set of all connected pixels, i.e., the outer air mask. Within the entire image range, for... By taking the complement set, the initial two-dimensional shell mask of the sample can be obtained. That is, all connected regions that "do not belong to the outside air", among which, This indicates the area where the sample is located (including the rock matrix and its internal pores). This is external air. Because this step excludes air that is connected to the boundary, it ensures that internal pores are not mistakenly removed as background.
[0136] S24, In order to make the shell mask To make it more coherent, this embodiment uses the initial shell mask. Perform a morphological closing operation (dilation followed by erosion) to fill small holes and smooth the edges.
[0137] Specifically, the morphological closing operation is performed as follows:
[0138] In order to appropriately fill the small holes inside the shell and smooth the boundaries without significantly changing the overall shape, this embodiment selects a circular (or approximately circular) structural element with a radius of about 2 pixels. The size of the structuring element can be adjusted appropriately according to the spatial resolution and sample scale, but it should generally not be too large to avoid excessive smoothing of details.
[0139] The formula for calculating morphological closing is as follows:
[0140]
[0141] in, This is a binary mask of the shell after morphological closing operations.
[0142] First, perform morphological closing operations on the set. The process involves expansion followed by corrosion to fill small, narrow cracks and holes in the shell mask, eliminate minor depressions at the boundaries, and make the shell profile more coherent and smooth.
[0143] S25. After the closing operation, combined with the connected component labeling, only the connected region with the largest number of voxels is retained as the final shell body to eliminate isolated noise blocks and prevent mislabeling of scattered small blocks.
[0144] Specifically, let the shell binary mask after morphological closing operation be... Perform 3D connected component labeling. Preferably, use 3D 26-connectivity (i.e., in 3D voxel space, all voxels that are adjacent vertically, horizontally, and diagonally are considered connected). Alternatively, 8-connectivity can be used as a simplified implementation on a 2D plane slice by slice. After labeling, for each foreground voxel ( Assign a connected component label as follows:
[0145]
[0146] in, This represents the total number of connected regions detected.
[0147] For each tag The number of voxels corresponding to these voxels is as follows:
[0148]
[0149] in, For the first The size of the connected component corresponding to each label.
[0150] Connected regions with the largest number of voxels That is, to choose to make The largest tag The corresponding connected regions are used as the main shell of the sample to automatically remove some scattered small pieces, noise, or isolated areas of non-target samples.
[0151] Based on the selected labels Construct the final shell mask :when hour, ,otherwise That is, only the voxel corresponding to the largest connected component is retained as the sample shell, and all other small blocks are regarded as background or invalid regions.
[0152] Furthermore, depending on the actual data, a second light etching can be selectively performed to remove local edge burrs, further improving the regularity of the boundary, thereby obtaining a complete and clean sample shell mask.
[0153] Since performing region growing or map cutting directly across the entire image area can easily disrupt internal connectivity due to the intrusion of external air, this embodiment specifically constructs a "sample shell mask" before the three-class classification. By applying an adaptive threshold to each slice to obtain candidate air regions, and combining morphological reconstruction and connectivity analysis, only external air voxels connected to the image boundary are removed, while retaining the true pores inside the sample. The shell mask thus strictly limits the spatial domain of subsequent decisions, eliminating the systematic interference caused by the connection between the reef limestone sample and the outside world during the preparation / scanning stage. This pre-constraint is the key difference of this embodiment for the reef limestone scene.
[0154] S3. After obtaining the final sample region mask, the sample grayscale distribution is modeled using a Gaussian Mixture Model (GMM) to automatically determine the segmentation threshold.
[0155] S31. Extract several slices from the three-dimensional volume data at uniform intervals. The specific extraction process is as follows:
[0156] Define the total number of slices and the number of targets extracted: Assume the reconstructed 3D CT volume data is in There are a total of directions Zhang slice ( In practical applications, the number of representative slices needed for training and statistics should be pre-defined. (For example , The size can be adjusted within a reasonable range based on the sample size and the amount of calculation.
[0157] Calculate the equal-interval step size: To ensure uniform sampling across the entire sample height, an approximate equal-interval step size is calculated as follows:
[0158]
[0159] Soon The slices are roughly divided into: There are several intervals, each with a thickness of approximately [missing information]. Zhang slice.
[0160] Automatically determine slice numbers at equal intervals: [The sentence is incomplete and likely refers to a specific slice number.] The serial number of the representative slide ( The definition is as follows:
[0161]
[0162] That is, a slice is selected at the "middle position" of each interval as a representative to ensure that the top, middle and bottom areas are covered, and the slices are relatively evenly distributed and not concentrated in a certain area.
[0163] By using the above-mentioned equidistant extraction method, a set of representative slices with reasonable selection and uniform spatial distribution can be automatically obtained without increasing human intervention. This avoids the instability caused by completely random extraction, as well as the subjectivity and workload of manual selection of each slice.
[0164] S32. In each slice, only the internal region of the sample (i.e., the shell mask) is considered. The grayscale values of a portion of the sample are randomly sampled. The specific process of random sampling is as follows:
[0165] Limit the sampling range and only take points within the sample: for the first A representative slice, only those already masked by the shell. Sampling is performed within the marked "sample interior region," meaning only samples meeting the criteria are sampled. The pixels in this internal region are included in the statistics; the total number of pixels in this internal region is denoted as . (That is, the number of pixels belonging to the sample in this slice is) This avoids mixing external air or invalid background into the training data.
[0166] Determine the sample size by sampling at a fixed ratio: Set a sampling ratio. (For example That is, about 5% of the pixels inside the sample of the slice are taken as grayscale samples; The ratio can be adjusted within a certain range depending on the sample size and computational cost, but generally a small proportion is sufficient to obtain a representative statistical sample. Number of sampling points per slice for:
[0167]
[0168] Sampling method: First, the internal region of the slice Pixels Number it, and then Randomly generated A set of unique integer indices Record the gray values corresponding to these indices as , which serves as the grayscale sample for this slice.
[0169] Furthermore, equivalent methods such as regular grid sampling can also be used, for example, "every..." The above uniform random sampling method can be replaced by the method of "taking one sample point per pixel" or "sampling on fixed grid nodes", as long as it can cover the internal area of the sample and ensure that the sample is representative.
[0170] S33. For each representative slice The sampling results of all representative slices within the sample's internal region are as follows. The data is then aggregated to form a unified grayscale training sample set. The calculation formula is as follows:
[0171]
[0172] in, For the merged first Each sample gray value, in other words, each The normalized gray value of a voxel reflects the distribution of the three microstructures—pores, transitional grays, and matrix—in the gray space; the total number of samples in the training sample set. The sum of the number of sampling points for each slice is calculated using the following formula:
[0173]
[0174] That is, after "stitching together" the grayscale samples of all representative slices, a total of 1 sample point.
[0175] S34. To characterize the above grayscale distribution, this embodiment assumes that it is composed of a weighted superposition of three Gaussian distributions. The Gaussian mixture model of the sample grayscale distribution is then established as follows:
[0176]
[0177] in, These correspond to three components (porosity, grayscale overlap / transition zone, and matrix). For the first The mixing weights of each Gaussian component satisfy the following conditions: ; For the first The mean of the Gaussian components; For the first The variance of each Gaussian component; the one-dimensional Gaussian distribution is represented as follows:
[0178]
[0179] S35. To obtain the above model parameters The EM algorithm is used to train the samples. Maximum likelihood fitting, i.e., maximizing the log-likelihood function, is performed as follows:
[0180]
[0181] The iterative steps of the EM algorithm are as follows:
[0182] Calculate the posterior probability (E-step): given the current parameters In the case of calculating each sample Belongs to the The posterior probability of each Gaussian component is calculated using the following formula:
[0183]
[0184] Update model parameters (M-step): Re-estimate the effective sample size and parameters of each component using the above posterior probabilities (using the current soft partitioning results, re-estimate the position, width, and proportion of each Gaussian curve to better match the actual grayscale distribution). The calculation formula is as follows:
[0185]
[0186] After completing one round of E-step and M-step, a new set of parameters is obtained, and then the next round of iteration begins. At the end of each iteration, the overall log-likelihood value is calculated based on the current parameters. If the change in goodness of fit between two consecutive iterations is very small (e.g., less than a preset threshold), it indicates that the improvement from continued iteration is negligible, and the algorithm is considered to have converged, and iteration can be stopped. If the number of iterations reaches the set upper limit (e.g., 100) before the algorithm converges, it is terminated early, and the current parameters are used as the result. The three sets of parameters obtained at the end of the iteration are... , which are the final parameters of the GMM used for subsequent double threshold estimation and classification, and the posterior probability curve is plotted.
[0187] The EM algorithm was used to perform maximum likelihood fitting on the three-component Gaussian mixture model to obtain the mixture weights, mean, and variance of each component. The three components can be respectively identified as the low-gray-level porosity component, the intermediate-gray-level overlapping or transitional component, and the high-gray-level matrix component, thus providing a statistical basis for subsequent threshold selection and segmentation.
[0188] Furthermore, other maximum likelihood optimization methods or Bayesian inference-based methods can be used to fit the GMM.
[0189] Furthermore, depending on the actual lithology or grayscale distribution characteristics, a two-component GMM (e.g., distinguishing only between pores and matrix) or a GMM with more than three components (further subdividing different types of matrix or pores) can also be used.
[0190] S36. Based on the posterior probability curves of each component of the GMM, this embodiment selects a posterior probability threshold of 0.90 (90%) to determine the dual thresholds for segmentation: when a gray value falls into the lowest gray component (corresponding to the pore class) and its posterior probability... At that time, this gray level is defined as the upper limit of pore gray level. When the gray value falls into the highest gray component (corresponding to the matrix class) and its posterior probability... At that time, this gray level is defined as the lower limit of the matrix gray level. .
[0191] From this obtained and Two thresholds divide the grayscale range into three stages: below The pore size is higher than that. The range of grayscale values is the matrix range, while the grayscale values between the two thresholds belong to the "uncertain" range to be classified.
[0192] Figure 2 The GMM fitting curve is illustrated, with the vertical dashed line indicating the calculated porosity threshold. and matrix threshold By using the posterior probability curves of the two components of the model as a basis, a double threshold is automatically given at the intersection of the posterior threshold p. and This allows for adaptive and reproducible threshold pairs across different slices and samples, avoiding cross-slice drift caused by manual parameter tuning.
[0193] S4. Using the above-mentioned dual thresholds, perform preliminary three-class segmentation on the CT body data.
[0194] Specifically, for each slice, iterate through its pixel grayscale values: when the grayscale value is... And this pixel is located in the sample mask ( When the value is within a certain range, it is labeled as a pore voxel; when the grayscale value is within a certain range, it is labeled as a pore voxel. When located within the sample, it is labeled as a matrix voxel; while for voxels within the sample whose grayscale values are between the two thresholds, they are temporarily labeled as undetermined voxels.
[0195] After dual-threshold segmentation, the core slices were divided into three regions: such as Figure 3 As shown, the candidate pore region (black area), the candidate matrix region (white area), and the undecided gray area are represented. The dual-threshold initial classification method is essentially an improved threshold segmentation method. By setting upper and lower thresholds, it effectively identifies the most clearly defined pore and matrix regions, while leaving the difficult-to-distinguish transitional regions for subsequent processing, thus significantly reducing the workload of manually adjusting the thresholds.
[0196] S5. To improve the robustness of the unknown region determination, morphological closing operation is performed on the binary results obtained from the initial classification to eliminate isolated small holes and fragmented noise, obtain the pixel set belonging to the pores and the pixel set belonging to the matrix, and generate a geometrically coherent "reliable seed" with higher confidence.
[0197] Specifically, such as Figure 4 As shown, closed-loop filling processes are performed on the candidate pore and matrix regions, respectively. More specifically, the regions marked as pores (black) are expanded by an expansion operation to enlarge their boundaries and fill any small cavities or breaks that may exist inside the pores. Subsequently, an erosion operation with an equivalent structural element is applied to eliminate the slight over-expansion caused by the expansion, without breaking the connected pore domains, thus preserving the overall shape after expansion. After the closing operation, the originally scattered or noise-induced small pore fragments are connected into larger connected pieces, resulting in reliable pore seed regions. Similarly, a closing operation is performed on the initial classification region of the matrix (white) to obtain a reliable matrix seed region. Reliable seeds are forced to mutually exclude each other, serving as a hard constraint in the subsequent energy model to suppress the mutual invasion and propagation from pores to the matrix or from the matrix to pores at the source.
[0198] The closing operation cleans up the initial results: it fills in small holes inside the pores, connects narrow gaps, removes isolated noise points, and makes the seed region of reliable pores clearer. and reliable substrate seed area These form continuous block-shaped regions, which are more geometrically connected and regular, making them easier to use for subsequent fine segmentation steps such as graph cut or region growing.
[0199] Reliable pore seeds Reliable matrix seeding will be achieved by confining the sample within the mask area and excluding any voxels already labeled as matrix. Similarly, the sample range is limited, and marked pores and those already belonging to the sample are excluded. The voxels generated from this. and These represent the sets of pixels that can be confidently identified as belonging to the pores and the matrix, respectively, and will serve as hard constraint seed points in the subsequent graph cut subdivision algorithm.
[0200] S6. For uncertain regions that are still unclassified, a graph cut-based energy minimization algorithm is used for fine segmentation. Specifically, to ensure efficiency, a graph cut model is constructed only within the region of interest of each slice. The goal of the graph cut model is to divide the graph into two disjoint subsets (foreground and background) using a minimum cut algorithm, thereby minimizing the total energy (the sum of the cut weights).
[0201] S61. Calculate the sample shell mask in the current slice. The smallest bounding rectangle of part 1 is defined, and a safety boundary of a predetermined width (halo, e.g., 16 pixels) is extended in all directions. The local rectangular region defined by this is denoted as the Region of Interest (ROI). The ROI covers the transition zone near the real boundary, ensuring the integrity of the topology at the boundary, and significantly reduces the number of nodes and edges that need to participate in graph construction. This causes the peak memory and time for graph cut solutions to decrease linearly with the ROI area. By constructing the graph model and performing segmentation only for the pixels within the ROI, the computational load can be significantly reduced. At the same time, the margin of the safety boundary ensures that the segmentation is not affected by boundary clipping. Within the ROI, a corresponding node is constructed for each pixel, and undirected neighborhood edges are added according to its connection relationship with neighboring pixels (in this embodiment, four-neighborhood is used, that is, each pixel node is connected to its top, bottom, left, and right neighboring pixels).
[0202] S62. In this planar mesh diagram, each neighborhood edge is assigned a smooth weight. Smoothing weights are used to characterize the "smoothing cost" when adjacent pixels are segmented into different categories. It reflects the similarity of adjacent pixels belonging to the same region.
[0203] Smoothing weights The setting and the grayscale difference between adjacent pixels Regarding this, let the standard deviation of the grayscale values within the ROI be... Then it has the following definition:
[0204]
[0205] in, Pixels and grayscale value; The smoothing coefficient (in this embodiment, it is taken as...) ).
[0206] As can be seen from the above formula, when the gray level difference between adjacent pixels... When I was very young, Approaching the maximum value means that in minimizing energy, pixels with similar gray levels tend to be grouped into the same category; conversely, if adjacent pixels have different gray levels... Very large, then Approaching 0 indicates that the cost of classifying pixels into different categories is very low, and graph cut allows for free cutting of boundaries at this high gradient.
[0207] By using this exponentially decaying weighting, the smoothing term of the graph cut can adapt to the gray-level unevenness of CT images: at the interface between the matrix and pores, where the gray-level contrast is large, It is very small, which helps to preserve actual edge details, while the gray-scale variation is small in homogeneous regions inside the rock. The larger the value, the smoother and more consistent the segmentation results in these regions.
[0208] S63. The graph cut model also incorporates grayscale probability information from the GMM model as data items. Specifically, it uses Gaussian distribution parameters for pore and matrix classes obtained from previous GMM fitting. and For each pixel within the ROI, calculate its negative log-likelihood cost for belonging to the pores and matrix, using the following formula:
[0209]
[0210] in, The cost of this pixel belonging to the aperture. The cost of this pixel belonging to the matrix. This represents the likelihood probability that the pixel belongs to the corresponding category.
[0211] The negative log-likelihood costs of the pores and matrix are appropriately scaled and used as connection weights for the pixel to connect to the background terminal (representing the pore class) and the foreground terminal (representing the matrix class), respectively. The connection weights reflect the probability that the pixel belongs to the foreground or the background. As a result, pixels with grayscale values closer to the pore distribution will incur a higher energy penalty when cut into the matrix, while pixels with grayscale values closer to the matrix distribution will have a higher cost when cut into pores.
[0212] In this embodiment, a graph model consisting of two types of energy terms is established within the ROI of each slice. The data term (corresponding to the terminal edge) adopts the negative log-likelihood of the two components of the Gaussian mixture model to measure the cost of a pixel belonging to "pore" or "matrix" respectively. The smoothing term (corresponding to the neighborhood edge) adopts the Potts model of four neighborhoods. This form provides strong consistency constraints in homogeneous regions and automatically weakens at strong gradients to preserve edges.
[0213] S64. The connection weights between the source (representing "foreground") / sink (representing "background") terminals and pixel nodes in the graph cut model, together with the smoothing weights of the adjacent edges mentioned above, constitute the energy function of graph cut segmentation. The graph cut algorithm achieves optimal binary labeling with global energy minimization through minimum cut-maximum flow solution, and optimally segments voxels in the unresolved region.
[0214] To ensure that the segmentation results strictly adhere to the constraints of the reliable seed region, the seed points are assigned infinitely large weights in the graphical model: this is used to define the pore seed region. All pixels within the matrix are connected to the foreground endpoint (matrix) by an edge with a weight approaching infinity. This means that if these nodes are incorrectly segmented into the matrix, an unacceptable cost will be incurred, thus forcing them to be connected to the background (pores) in the optimal solution, i.e., always segmented as pores. Similarly, the matrix seed region... The inner pixel is connected to the background terminal with infinite weights to force it to remain as the matrix.
[0215] Furthermore, for the ROI edge area outside the shell mask (i.e., the original...) The portion that is actually the external air is also given infinite weight in this embodiment, thereby ensuring that all pixels outside the sample are determined to be background pores.
[0216] With the above settings, the image cut optimization process will strictly follow the guidance of known reliable regions, and will only automatically classify pixels in "uncertain" regions as pores or matrix. Figure 5The results of the slice segmentation after graph cut subdivision are shown. The gray blurred areas have been automatically assigned as pores (black) or matrix (white) based on the principle of minimizing neighborhood grayscale and overall energy. The pore / matrix interface is smoother and more accurate, without over-segmentation or under-segmentation. By performing the above graph cut segmentation process on each slice, the three-dimensional segmentation of the entire volume data can be completed.
[0217] Furthermore, this embodiment employs multi-process parallel processing to accelerate the processing of each slice: for example, four parallel sub-processes are started to process different slice regions respectively, while the main process collects the results, thereby fully utilizing the multi-core CPU to improve computational efficiency. After all slices are processed, the segmentation results are reconstructed into a three-dimensional binary core model. In this three-dimensional model, voxel values "1" represent pore space, "2" represent rock matrix, and "0" represent the external background of the sample, achieving a clear separation of the internal pores and matrix of the rock. Based on the final three-dimensional segmentation data, parameters such as the porosity of the sample can also be automatically calculated; for example, the overall porosity can be obtained by calculating the volume fraction of pore voxels in all slices. In addition, this embodiment saves the progressively generated process output files during the segmentation process, facilitating result verification and further utilization. For example, the program uses the above sample mask... The preliminary three-class classification results, seed regions, and final segmentation results are exported and saved as 3D data files. The time consumed for each slice in each step (normalization, masking, threshold classification, seed generation, graph cut subdivision, etc.) is recorded (output time log file), along with GMM model parameters and posterior curve data (output JSON and NPZ files). By retaining these procedural model files, users can check the intermediate results at each step to verify the effectiveness of the method or use them for subsequent analysis.
[0218] Finally, the obtained binary digital core model is exported to a commonly used 3D simulation mesh format, such as STL or VTK. A typical approach is to perform isosurface extraction algorithms (e.g., Marching Cubes) on the binary volume data to extract the interface between the pores and matrix as a set of triangular facets and save it in STL format; or directly convert the voxel data to a VTK structured mesh, retaining the classification identifier of each voxel. The exported STL model can be used to construct surface meshes for pore space (e.g., for defining boundary conditions in hydrodynamic simulations), while the VTK format volume data can be directly imported into numerical simulations such as finite difference and lattice Boltzmann. Regardless of the format used, the 3D digital core model obtained after the above segmentation and processing realistically reproduces the internal pore structure in the original CT data and is more regular and cleaner due to the removal of noise and small artifacts. By exporting this digital model in a general simulation format, a connection between CT images and numerical simulations is achieved, which can be used in subsequent fields such as seepage coupling simulation and rock mechanics analysis, providing strong support for the study of the internal microstructure and properties of rocks.
[0219] This embodiment employs a collaborative path of "statistically robust dual threshold determination—sample shell mask—reliable seed—constrained graph model—ROI and halo control scaling—memory mapping and parallel pipeline" to address the external air connectivity interference, gray-level overlap, and topological connectivity requirements of reef limestone CT. It provides a systematic solution that can operate stably in ordinary computing environments without reducing resolution. Compared to the naive approach of directly transferring general segmentation algorithms to different objects, the various stages of this embodiment are coupled together: the shell mask confines the decision domain to a physically reasonable space; the closed-loop seed provides high-confidence constraints; the data items derived from GMM and the adaptive smoothing items form an energy model within the ROI that balances noise suppression and edge preservation; and finally, a closed loop with procedural data and parameterized output meets the engineering requirements of digital core analysis from image to computable model. This solution balances segmentation accuracy, topology preservation, computational constraints, and parameterization availability, and is suitable for data processing scenarios involving reef limestone and carbonate rock masses with similar external connectivity and gray-level overlap characteristics.
Claims
1. A method of three-dimensional modeling of reef limestone microstructure, characterized in that, Includes the following steps: S1. Scan and acquire the original grayscale image of the core sample, and preprocess the original grayscale image; S2. Based on the preprocessed grayscale image, extract the sample shell mask; S3. Establish a Gaussian mixture model of sample grayscale distribution and determine the segmentation threshold; S4. Classify and segment the volume data using a segmentation threshold to obtain candidate regions for pores, candidate regions for matrix, and unresolved regions. S5. Perform morphological closing operations on the pore candidate region and the matrix candidate region to obtain the pixel set belonging to the pore and the pixel set belonging to the matrix. S6. For unresolved areas, a graph-cut-based energy minimization algorithm is used for classification and segmentation; In step S2, the process of extracting the sample shell mask is as follows: S21. Calculate and select a global threshold to identify air candidate regions; S22. Perform morphological opening operation on the air candidate region to obtain the mask image; S23. Perform morphological reconstruction using the image boundary as a seed to obtain the initial shell mask; S24. Perform morphological closing operation on the initial shell mask; S25. The connected region with the largest number of voxels is retained as the final sample shell mask; In step S21, the global threshold is automatically calculated using the Otsu method, as follows: For any CT slice, the grayscale value will be... Quantified as The grayscale values were analyzed at several levels, and the frequency of each grayscale level was counted to obtain a normalized histogram. It satisfies the following formula: , For a certain candidate threshold Using this gray level as the boundary, the pixels are divided into the following two categories: Low grayscale Grayscale ; High grayscale Grayscale ; The intra-class weights, intra-class average gray levels, and global average gray levels of low-gray-level pixels and high-gray-level pixels are calculated separately. The calculation process is as follows: Low grayscale Intra-class weight The calculation formula is as follows: , High grayscale Intra-class weight The calculation formula is as follows: ， Low grayscale Intra-class average gray level The calculation formula is as follows: ， High grayscale Intra-class average gray level The calculation formula is as follows: , Global average gray level The calculation formula is as follows: ， Calculate the candidate threshold Inter-class variance The calculation formula is as follows: , For possible candidate thresholds Calculate the between-class variances one by one Selecting classes with different variances Maximum candidate threshold As a global threshold The calculation formula is as follows: ， In S22, the morphological opening operation is performed as follows: Let the air candidate region image be The neighborhood range is The formula for calculating the opening operation is as follows: , in, For mask images; In S23, the morphological reconstruction process is as follows: Let the image boundary seed be the labeled image. Select structural elements Then there is an expansion operator. The calculation formula is as follows: , Let the initial state Then there is an iteration to the th When walking, ; When the iteration reaches a certain point that satisfies... When the iteration stops, the final result of the iteration is the external air mask. External air mask Take the complement to obtain the initial shell mask. ; In S24, the morphological closing operation is performed as follows: Selecting structural elements Then the formula for calculating the closing operation is as follows: , in, This is a binary mask of the shell after morphological closing operations. In S25, the process for determining the connected region with the largest number of voxels is as follows: Binary mask of the shell after morphological closing operation Perform 3D connected component labeling, assigning a connected component label to each foreground voxel as follows: ， in, This represents the total number of connected regions. For each tag The number of voxels in the corresponding connected regions is counted using the following formula: ， in, For the first The size of the connected region corresponding to each label; make The largest tag The corresponding region is the connected region with the largest number of voxels.
2. The three-dimensional modeling method for the microstructure of reef limestone according to claim 1, characterized in that, In step S1, the preprocessing includes grayscale value normalization, and the grayscale value normalization process is as follows: Let the gray value of the low quantile be... The gray value of the high quantile is grayscale value for each voxel The linear stretching normalization is calculated using the following formula: ， in, For normalized volume data, Indicates the position of each voxel; Map the normalized volume data to The grayscale value outside the specified range is truncated.
3. The three-dimensional modeling method for the microstructure of reef limestone according to claim 1, characterized in that, In step S3, the process of establishing a Gaussian mixture model of the sample grayscale distribution and determining the separation threshold is as follows: S31. Extract slices from the volume data at equal intervals; S32. Randomly sample the grayscale values of the shell mask area in each slice; S33. All sampling results are summarized in a unified manner to form an overall grayscale training sample set; S34. Establish a Gaussian mixture model; S35. Perform maximum likelihood fitting on the gray training sample set, obtain the model parameters, and plot the posterior probability curve. S36. Determine the segmentation threshold based on the posterior probability curve.
4. The three-dimensional modeling method for the microstructure of reef limestone according to claim 3, characterized in that, In step S31, the process of extracting slices is as follows: Let the total number of slices be The number of representative slices required for training and statistics is Calculate the equal interval step size The calculation formula is as follows: ， The slice number is automatically determined based on equal intervals, calculated using the following formula: ， in, For the first The serial number of the representative slice, ; In step S32, the random sampling process is as follows: Set sampling ratio Calculate the first Number of sampling points per slice The calculation formula is as follows: , in, The total number of pixels within the shell mask area; Will Pixels Number it, and then Randomly generated Let there be a set of unique integer indices, and denote the grayscale values corresponding to these indices as... , as the grayscale sample of the slice; In S33, the unified calculation formula is as follows: ， in, This is the overall grayscale training sample set; For the merged first Individual sample grayscale values; Gray training sample set Total number of samples The sum of the number of sampling points for each slice is calculated using the following formula: , In S34, the Gaussian mixture model is established as follows: , in, These correspond to three components: pores, grayscale overlap area, and matrix, respectively. For the first The mixing weight of each Gaussian component; For the first The mean of the Gaussian components; For the first The variance of each Gaussian component; In S35, the maximum likelihood fitting process is as follows: The EM algorithm is used to process the grayscale training sample set. Maximum likelihood fitting, i.e., maximizing the log-likelihood function, is performed as follows: ， The iterative steps of the EM algorithm are as follows: Calculate the posterior probability: given the current parameters In the case of calculating each sample Belongs to the The posterior probability of each Gaussian component is calculated using the following formula: ， Update model parameters: Re-estimate the effective sample size and parameters of each Gaussian component using the above posterior probabilities, as shown in the following formula: ， ， After completing one round of iteration, a new set of parameters is obtained, and then the next round of iteration continues; After each iteration, the overall log-likelihood value is calculated based on the current parameters. The iteration stops when the goodness of fit changes very little between two consecutive iterations. If the number of iterations reaches the set upper limit before the algorithm converges, the iteration is terminated early. The three sets of parameters obtained at the end of the iteration are then used to determine the final result. These are the final model parameters; In step S36, the process of determining the segmentation threshold is as follows: The posterior probability threshold is selected based on the posterior probability curve. When a gray value falls into the lowest gray component and its posterior probability is not less than the posterior probability threshold, the gray value is defined as the porosity threshold. When a gray value falls into the highest gray component and its posterior probability is not less than the posterior probability threshold, the gray value is defined as the matrix threshold.
5. The three-dimensional modeling method for the microstructure of reef limestone according to claim 1, characterized in that, In step S4, the classification and segmentation process is as follows: For each slice, its pixel grayscale values are iterated; when the grayscale value is less than or equal to the pore threshold and the pixel is located within the shell mask, it is marked as a pore voxel; when the grayscale value is greater than or equal to the matrix threshold and the pixel is located within the shell mask, it is marked as a matrix voxel; when the pore threshold is less than or equal to the grayscale value less than or equal to the matrix threshold and the pixel is located within the shell mask, it is marked as a pending voxel; thus, the slice is classified into pore candidate regions, matrix candidate regions, and pending regions.
6. The three-dimensional modeling method for the microstructure of reef limestone according to claim 1, characterized in that, In step S6, the segmentation and classification process is as follows: S61. Define the region of interest, and construct a corresponding node for each pixel within the region of interest; S62. Add neighborhood edges to adjacent pixels and assign a smooth weight to each neighborhood edge; S63. For each pixel in the region of interest, calculate its negative log-likelihood cost belonging to the pores and matrix, and use it as the connection weight for connecting the background terminal and the foreground terminal, respectively. S64. Based on the smoothing weight and connection weight, establish an energy function, apply the maximum flow algorithm to extract the minimum cut, and perform optimal segmentation of voxels in the unresolved region.
7. The three-dimensional modeling method for the microstructure of reef limestone according to claim 6, characterized in that, In step S61, the process of establishing the region of interest is as follows: Calculate the minimum bounding rectangle of the shell mask in the current slice, and extend the safety boundary to a predetermined width in all directions. The local rectangular region determined in this way is denoted as the region of interest (ROI). In S62, smoothing weights The calculation formula is as follows: , in, The difference in grayscale between adjacent pixels; Pixels and grayscale value; For smoothing coefficients; The standard deviation of grayscale values within the region of interest (ROI); In S63, the formula for calculating the negative log-likelihood cost of pores is as follows: ， in, The cost of this pixel belonging to the aperture. This represents the likelihood probability that the pixel belongs to the corresponding category. The parameters are Gaussian distribution parameters for pores; The formula for calculating the negative log-likelihood cost of the matrix is as follows: , in, The cost of this pixel belonging to the matrix. These are the parameters of a matrix-type Gaussian distribution.
8. The three-dimensional modeling method for the microstructure of reef limestone according to claim 7, characterized in that, In step S64, before performing optimal segmentation on the voxels in the unresolved region, the set of pore pixels is connected to the foreground terminal by an edge with a weight approaching infinity, forcing it to be connected to the background in the optimal segmentation so that it is always classified as pores; the set of matrix pixels is connected to the background terminal by an edge with a weight approaching infinity, forcing it to be connected to the foreground in the optimal segmentation so that it is always classified as matrix.