A Fast Workﬂow for Automatically Extracting the Apparent Attitude of Fractures in 3-D Digital Core Images

: Fractures play a crucial role as ﬂuid conduits and reservoir spaces in reservoirs. The distribution and characteristics of fractures determine the presence of high-quality reservoirs. To accurately analyze and observe fracture parameters, three-dimensional (3-D) digital cores generated from computed tomography (CT) are utilized. However, the current process of extracting fracture properties from these digital cores is time-consuming and labor-intensive. This paper introduces a new, fast, and automatic workﬂow for extracting the apparent dip angle and direction of fractures from 3-D digital core images. The proposed workﬂow involves several steps. Firstly, two perpendicular cross-sections are obtained from the digital core and converted into binary images. Next, the coordinates of four fracture feature points within the core image are automatically extracted. The fracture plane is then ﬁtted using the least squares method based on the extracted coordinates. Finally, the apparent dip angle and direction of the fracture are calculated using the plane’s normal vector. By comparing and analyzing the proposed workﬂow with the original method, it becomes evident that the method proposed in this paper allows for quick, automated, and accurate extraction of the apparent dip angle and direction of fractures. The application of this workﬂow to extract fracture attitudes in 3-D micro-CT and full-hole digital core images signiﬁcantly enhances efﬁciency.


Introduction
Fractures are present in both conventional and unconventional reservoirs, occurring at various scales. Accurately identifying fractures is essential for enhancing reservoir fluid migration characteristics and increasing productivity. Therefore, studying the development of reservoir fractures is of great significance [1][2][3][4][5]. Numerous researchers worldwide have made considerable efforts in fracture research [6][7][8][9], but most studies have focused on qualitative analysis [10][11][12][13][14]. Although some scholars have used quantitative methods to study fractures, the number of such studies is limited. Therefore, it is crucial to quickly and effectively identify fracture development characteristics and quantitatively characterize reservoir fractures.
Previous studies primarily utilized a combination of core samples, thin sections, and the Formation Micro-Imager (FMI) log data to analyze fracture characteristics. These studies quantitatively described fracture parameters in the study area and predicted fractures using numerical simulation techniques [15][16][17][18][19][20]. However, thorough research of the literature has revealed several limitations associated with these methods. These limitations include weathering and denudation of the outcrop, susceptibility to structure and lithology influences, operator subjectivity, and time-consuming processes [21][22][23][24][25].

The Original Workflow
The original workflow includes data pre-processing and parameter extraction were referenced from a commercial software. As depicted in Figure 1, the data p cessing involves data trimming, filtering, and segmentation to address artifacts an present in the X-ray CT reconstructed data. Following the pre-processing stage, t ture space distribution is obtained, and, subsequently, the fracture apparent dip an direction are calculated. After applying X-rays to the samples, CT images of the samples can be obtain ure 2a). However, issues such as tilting and fuzzy edges may arise, necessitating th ping of the samples for easier subsequent processing (Figure 2b). The CT data ar to artifacts, noise, and other interference factors caused by the scanning X-ray. To m the influence of these factors, it is essential to filter the data before conducting th segmentation. Fractures typically have lower grayscale values compared to the m other components in the CT images. As illustrated in Figure 2b, the lower grays gions correspond to fractures. Hence, threshold segmentation can be employed to entiate fractures from other components. In some cases, additional methods, like hat algorithm or watershed segmentation, may be required to process regions wit mal grayscale changes. Through threshold segmentation, all pixels belonging to fr are categorized together, while the matrix or minerals are grouped separately. In 2c, the blue regions represent fractures, while the remaining areas represent the m

Data Pre-Processing in the Original Workflow
After applying X-rays to the samples, CT images of the samples can be obtained (Figure 2a). However, issues such as tilting and fuzzy edges may arise, necessitating the cropping of the samples for easier subsequent processing (Figure 2b). The CT data are prone to artifacts, noise, and other interference factors caused by the scanning X-ray. To mitigate the influence of these factors, it is essential to filter the data before conducting threshold segmentation. Fractures typically have lower grayscale values compared to the matrix or other components in the CT images. As illustrated in Figure 2b, the lower grayscale regions correspond to fractures. Hence, threshold segmentation can be employed to differentiate fractures from other components. In some cases, additional methods, like the top-hat algorithm or watershed segmentation, may be required to process regions with minimal grayscale changes. Through threshold segmentation, all pixels belonging to fractures are categorized together, while the matrix or minerals are grouped separately. In Figure 2c, the blue regions represent fractures, while the remaining areas represent the matrix.

Extraction Method in the Original Workflow: Moments of Inertia Method
After the pre-processing, the fracture dip angle and direction can be obtained usi the moments of inertia method (MIM). The MIM method is commonly used to calcul the global characteristics of a given shape. The main step in the MIM for extracting t apparent dip angle is to compute the first-order moments and second-order moments the fracture surface. The first-order moments determine the centroid or center of mass. the discrete case, the first-order moments are defined as follows: where is the area and ( , , ) is a point in the object. The second-order moments are defined in the discrete case as: The fracture's normal vector points to the direction of its major inertia axis. It is d termined by finding the eigenvector corresponding to the largest eigenvalue of the iner matrix:

Extraction Method in the Original Workflow: Moments of Inertia Method
After the pre-processing, the fracture dip angle and direction can be obtained using the moments of inertia method (MIM). The MIM method is commonly used to calculate the global characteristics of a given shape. The main step in the MIM for extracting the apparent dip angle is to compute the first-order moments and second-order moments of the fracture surface. The first-order moments determine the centroid or center of mass. In the discrete case, the first-order moments are defined as follows: where A(X) is the area and (x i , y j , z k ) is a point in the object. The second-order moments are defined in the discrete case as: Processes 2023, 11, 2517

of 16
The fracture's normal vector points to the direction of its major inertia axis. It is determined by finding the eigenvector corresponding to the largest eigenvalue of the inertia matrix: The definition of the apparent dip angle of the fracture is the angle, ϕ, between the fracture's normal vector, v, and the coordinate Z-axis, and the apparent dip direction is the angle, θ, between the projection of v (vector v') on plane XY and the X-axis (Figure 3). The dip angle, ϕ, ranges from 0 • to 90 • , and the dip direction, θ, ranges from 0 • to 360 • . The definition of the apparent dip angle of the fracture is the angle, , between the fracture's normal vector, v, and the coordinate Z-axis, and the apparent dip direction is the angle, , between the projection of v (vector v') on plane XY and the X-axis (Figure 3). The dip angle, , ranges from 0° to 90°, and the dip direction, , ranges from 0° to 360°.

The New Workflow
A new comprehensive workflow for extracting the dip angle and direction of fractures has been established in this study. Figure 4 illustrates the flowchart consisting of two main parts: data pre-processing and parameter extraction. When comparing our workflow to the original one, the pre-processing section of this workflow is simpler and more automated. To extract the fracture parameters using the method described in this paper, the first step is to obtain two orthogonal 2-D images of the core. Next, the coordinates of all the fractures in the 2-D images are automatically extracted. Finally, the least squares method is utilized to obtain the fitted fracture plane and extract the fracture's apparent dip angle and direction. In the new workflow, in order to obtain a more accurate fitting plane for the fractures, several steps are undertaken. Firstly, two non-parallel planes are selected to extract the

The New Workflow
A new comprehensive workflow for extracting the dip angle and direction of fractures has been established in this study. Figure 4 illustrates the flowchart consisting of two main parts: data pre-processing and parameter extraction. When comparing our workflow to the original one, the pre-processing section of this workflow is simpler and more automated. To extract the fracture parameters using the method described in this paper, the first step is to obtain two orthogonal 2-D images of the core. Next, the coordinates of all the fractures in the 2-D images are automatically extracted. Finally, the least squares method is utilized to obtain the fitted fracture plane and extract the fracture's apparent dip angle and direction. The definition of the apparent dip angle of the fracture is the angle, , between the fracture's normal vector, v, and the coordinate Z-axis, and the apparent dip direction is the angle, , between the projection of v (vector v') on plane XY and the X-axis ( Figure 3). The dip angle, , ranges from 0° to 90°, and the dip direction, , ranges from 0° to 360°.

The New Workflow
A new comprehensive workflow for extracting the dip angle and direction of fractures has been established in this study. Figure 4 illustrates the flowchart consisting of two main parts: data pre-processing and parameter extraction. When comparing our workflow to the original one, the pre-processing section of this workflow is simpler and more automated. To extract the fracture parameters using the method described in this paper, the first step is to obtain two orthogonal 2-D images of the core. Next, the coordinates of all the fractures in the 2-D images are automatically extracted. Finally, the least squares method is utilized to obtain the fitted fracture plane and extract the fracture's apparent dip angle and direction. In the new workflow, in order to obtain a more accurate fitting plane for the fractures, several steps are undertaken. Firstly, two non-parallel planes are selected to extract the  In the new workflow, in order to obtain a more accurate fitting plane for the fractures, several steps are undertaken. Firstly, two non-parallel planes are selected to extract the coordinates of the fracture feature points. Due to the lack of actual data, the orthogonal planes XZ and YZ in the 3-D data are chosen. Next, the orthogonal 2-D images are subjected to separate image binarization processing. Additionally, the automatic picking of fracture coordinates is crucial in our method. This is accomplished by traversing each connected domain, and all the fracture coordinates are automatically extracted using the minimum external rectangle algorithm. The algorithm flow is as follows:

1.
Firstly, the binary image is analyzed for connected domains, specifically identifying the fracture region; 2.
Then, the minimum external rectangle of each connected domain is extracted. As shown in Figure 5, the coordinates of the upper left corner of the external rectangle, as well as the length and width of the rectangle, are recorded. The points A and E represent the upper left corners of the minimum external rectangles, with coordinates ( x A , y A ) and (x E , y E ), respectively. The variables a and c represent the width of each rectangle, while b and d represent the length; 3.
Finally, the tilt direction of the connected domain is determined. If it is sloping downward, the coordinates for point A are (x A , y A ), and the coordinates for point B are (x A + b, y A − a) (as shown in Figure 5a). If it is sloping upward, the coordinates for point C are (x E , y E − c), and the coordinates for point D are (x E + d, y E ) (as shown in Figure 5b).
rocesses 2023, 11, x FOR PEER REVIEW 6 coordinates of the fracture feature points. Due to the lack of actual data, the ortho planes XZ and YZ in the 3-D data are chosen. Next, the orthogonal 2-D images are jected to separate image binarization processing. Additionally, the automatic picki fracture coordinates is crucial in our method. This is accomplished by traversing each nected domain, and all the fracture coordinates are automatically extracted using the imum external rectangle algorithm. The algorithm flow is as follows: 1. Firstly, the binary image is analyzed for connected domains, specifically identi the fracture region; 2. Then, the minimum external rectangle of each connected domain is extracte shown in Figure

Extraction Method in the New Workflow: Least Squares Method
After extracting the coordinates of the fracture feature points, the fracture pla fitted using the least squares method. The fracture coordinates are used as inputs fo fitting. The angle between the normal vector of the plane and the positive direction Z-axis represents the apparent dip angle of the fracture, and the angle between the jection of the normal vector on plane XY and the X-axis represents the apparent dip d tion of the fracture. In Figure 6, the blue rectangles represent orthogonal 2-D planes inclined rectangle P represents the actual fracture plane. The four black dots represen feature points of the fracture in the orthogonal plane. The letter v denotes the plane' mal vector, and represents the apparent fracture dip angle. The letter v' denote projection vector of the plane's normal vector, and represents the apparent fractur direction.

Extraction Method in the New Workflow: Least Squares Method
After extracting the coordinates of the fracture feature points, the fracture plane is fitted using the least squares method. The fracture coordinates are used as inputs for this fitting. The angle between the normal vector of the plane and the positive direction of the Z-axis represents the apparent dip angle of the fracture, and the angle between the projection of the normal vector on plane XY and the X-axis represents the apparent dip direction of the fracture. In Figure 6, the blue rectangles represent orthogonal 2-D planes. The inclined rectangle P represents the actual fracture plane. The four black dots represent the feature points of the fracture in the orthogonal plane. The letter v denotes the plane's normal vector, and ϕ represents the apparent fracture dip angle. The letter v' denotes the projection vector of the plane's normal vector, and θ represents the apparent fracture dip direction. Let the equation of the plane be: , , For a series of n points (  3): , , , 0, 1, … , 1. The n points are used to fit the plane equation, and ∑ a x a y a z must be minimized. To minimize , 0, 0,1,2.
So, it is equal to this: Rewrite to matrix form, then: The parameters , , can be obtained by solving Equation (14), and the plane equation can be obtained by substituting them into Equation (12).
According to the plane equation, the plane normal vector, 1 , then the equation of the apparent dip angle, , is: And the equation of the apparent dip direction, , is: To verify the effectiveness of this method, two plate fracture models, model 1 and model 2, were established (Figures 7a and 8a). These models have dimensions of 200 × 200 × 200 and have dip angles of 30° and 55°, respectively. And the dip directions of both Let the equation of the plane be: For a series of n points (n ≥ 3): (x i , y i , z i ), i = 0, 1, . . . , n − 1.
The n points are used to fit the plane equation, and S = ∑ n−1 i=0 (a 0 x + a 1 y + a 2 − z) 2 must be minimized. To minimize S, ∂ S ∂ a k = 0, k = 0, 1, 2.
So, it is equal to this: Rewrite to matrix form, then: The parameters a 0 , a 1 , a 2 can be obtained by solving Equation (14), and the plane equation can be obtained by substituting them into Equation (12).
According to the plane equation, the plane normal vector, v = [−a 0 −a 1 1], then the equation of the apparent dip angle, dip, is: And the equation of the apparent dip direction, dir, is: To verify the effectiveness of this method, two plate fracture models, model 1 and model 2, were established (Figures 7a and 8a). These models have dimensions of 200 × 200 × 200 and have dip angles of 30 • and 55 • , respectively. And the dip directions of both models are 0 • . For both model 1 and model 2, slices at plane XZ and YZ (specifically, slice 100) were selected, and the coordinates at the edges of the fractures were extracted (Figures 7b and 8b). In these figures, the blue region represents the fracture, and the red rectangles represent the XZ and YZ planes. The red circles indicate the coordinates of the feature points of fractures. Figure 7c shows the 3-D shape diagram of model 1's fracture, fitted using the proposed method. The blue plane represents the fitted fracture, while the orange points represent the locations of the extracted fracture feature points. The apparent dip angle of model 1's fracture was determined to be 29.81 • , and the apparent dip direction is extremely close to 0 • . Figure 8b depicts the specific extraction positions of the fracture coordinate points in the 100th slice of the XZ and YZ planes of model 2. Figure 8c illustrates the 3-D shape diagram of model 2's fracture, fitted using the least squares method. The apparent dip angle of model 2's fracture was found to be 55.07 • , and the apparent dip direction is extremely close to 0 • . The absolute errors for these two models' dip angles were calculated to be 0.19 and 0.07, while the relative errors were determined to be 0.64% and 0.12%, respectively. This result shows that when the fractures are relatively smooth, the error in extracting the fracture dip angle and direction using this method is negligible. The application of this method to both model 1 and model 2 demonstrates its efficacy in smooth fractures. In the subsequent section, this workflow will be applied to practical field data.  ures 7b and 8b). In these figures, the blue region represents the fracture, and the red rectangles represent the XZ and YZ planes. The red circles indicate the coordinates of the feature points of fractures. Figure 7c shows the 3-D shape diagram of model 1's fracture, fitted using the proposed method. The blue plane represents the fitted fracture, while the orange points represent the locations of the extracted fracture feature points. The apparent dip angle of model 1's fracture was determined to be 29.81°, and the apparent dip direction is extremely close to 0°. Figure 8b depicts the specific extraction positions of the fracture coordinate points in the 100th slice of the XZ and YZ planes of model 2. Figure 8c illustrates the 3-D shape diagram of model 2's fracture, fitted using the least squares method. The apparent dip angle of model 2's fracture was found to be 55.07°, and the apparent dip direction is extremely close to 0°. The absolute errors for these two models' dip angles were calculated to be 0.19 and 0.07, while the relative errors were determined to be 0.64% and 0.12%, respectively. This result shows that when the fractures are relatively smooth, the error in extracting the fracture dip angle and direction using this method is negligible. The application of this method to both model 1 and model 2 demonstrates its efficacy in smooth fractures. In the subsequent section, this workflow will be applied to practical field data.  models are 0°. For both model 1 and model 2, slices at plane XZ and YZ (specifically, slice 100) were selected, and the coordinates at the edges of the fractures were extracted (Figures 7b and 8b). In these figures, the blue region represents the fracture, and the red rectangles represent the XZ and YZ planes. The red circles indicate the coordinates of the feature points of fractures. Figure 7c shows the 3-D shape diagram of model 1's fracture, fitted using the proposed method. The blue plane represents the fitted fracture, while the orange points represent the locations of the extracted fracture feature points. The apparent dip angle of model 1's fracture was determined to be 29.81°, and the apparent dip direction is extremely close to 0°. Figure 8b depicts the specific extraction positions of the fracture coordinate points in the 100th slice of the XZ and YZ planes of model 2. Figure 8c illustrates the 3-D shape diagram of model 2's fracture, fitted using the least squares method. The apparent dip angle of model 2's fracture was found to be 55.07°, and the apparent dip direction is extremely close to 0°. The absolute errors for these two models' dip angles were calculated to be 0.19 and 0.07, while the relative errors were determined to be 0.64% and 0.12%, respectively. This result shows that when the fractures are relatively smooth, the error in extracting the fracture dip angle and direction using this method is negligible. The application of this method to both model 1 and model 2 demonstrates its efficacy in smooth fractures. In the subsequent section, this workflow will be applied to practical field data.

Application in Field Data and Results
To verify the reliability and practicality of this workflow, we conducted comparison experiments using micro-CT and full-hole digital core samples. In these experiments, we utilized the MIM method and the least squares method, respectively.

Practical Data
In this paper, we selected one micro-CT sample and six full-hole digital core samples for fracture analysis. Sample 1 consists of micro-CT data with irregular fractures, while samples 2-7 comprise CT data of the 3-D full-hole digital core. Please refer to Table 1 for more detailed information about all the samples. Additionally, Figure 9 displays the micro-CT digital core image of sample 1, while Figure 10a-f depict the original 3-D full-hole digital core images of samples 2-7.

Application in Field Data and Results
To verify the reliability and practicality of this workflow, we conducted comparison experiments using micro-CT and full-hole digital core samples. In these experiments, we utilized the MIM method and the least squares method, respectively.

Practical Data
In this paper, we selected one micro-CT sample and six full-hole digital core samples for fracture analysis. Sample 1 consists of micro-CT data with irregular fractures, while samples 2-7 comprise CT data of the 3-D full-hole digital core. Please refer to Table 1 for more detailed information about all the samples. Additionally, Figure 9 displays the micro-CT digital core image of sample 1, while Figure 10a-f depict the original 3-D full-hole digital core images of samples 2-7.

Application in Field Data and Results
To verify the reliability and practicality of this workflow, we conducted comparison experiments using micro-CT and full-hole digital core samples. In these experiments, we utilized the MIM method and the least squares method, respectively.

Practical Data
In this paper, we selected one micro-CT sample and six full-hole digital core samples for fracture analysis. Sample 1 consists of micro-CT data with irregular fractures, while samples 2-7 comprise CT data of the 3-D full-hole digital core. Please refer to Table 1 for more detailed information about all the samples. Additionally, Figure 9 displays the micro-CT digital core image of sample 1, while Figure 10a-f depict the original 3-D full-hole digital core images of samples 2-7.

Micro-Fracture in Micro-CT Sample
Firstly, the fracture morphology of sample 1 was extracted using the original method ( Figure 11a). The apparent dip angle of the fracture was calculated as 65.28° using MIM, and the apparent dip direction of the fracture was calculated as 37.11°. Subsequently, frame 400 on the XZ and YZ planes of the sample was selected to extract the feature point coordinates of fracture (Figure 11b). When dealing with irregular fractures, two ideas were considered for selecting the feature points: (1) selecting a feature point closer to the edge based on the connectivity domain, or (2) using the medial-axis transformation or refinement algorithm to extract the skeleton of the fractures and then extracting the feature points of the skeleton. With the new workflow, which used the least squares method, a fitted 3-D fracture morphology was obtained with the dip angle of 64.87° and dip direction of 35.73° for sample 1 (Figure 11c).

Macro-Fractures in Full-Hole Digital Cores
We conducted a comparative analysis of the accuracy between the new method and the original method for extracting the apparent fracture dip angle and direction. In accordance with the methodology, we utilized the original workflow to extract the fracture

Micro-Fracture in Micro-CT Sample
Firstly, the fracture morphology of sample 1 was extracted using the original method (Figure 11a). The apparent dip angle of the fracture was calculated as 65.28 • using MIM, and the apparent dip direction of the fracture was calculated as 37.11 • . Subsequently, frame 400 on the XZ and YZ planes of the sample was selected to extract the feature point coordinates of fracture (Figure 11b). When dealing with irregular fractures, two ideas were considered for selecting the feature points: (1) selecting a feature point closer to the edge based on the connectivity domain, or (2) using the medial-axis transformation or refinement algorithm to extract the skeleton of the fractures and then extracting the feature points of the skeleton. With the new workflow, which used the least squares method, a fitted 3-D fracture morphology was obtained with the dip angle of 64.87 • and dip direction of 35.73 • for sample 1 (Figure 11c).

Micro-Fracture in Micro-CT Sample
Firstly, the fracture morphology of sample 1 was extracted using the original method (Figure 11a). The apparent dip angle of the fracture was calculated as 65.28° using MIM, and the apparent dip direction of the fracture was calculated as 37.11°. Subsequently, frame 400 on the XZ and YZ planes of the sample was selected to extract the feature point coordinates of fracture (Figure 11b). When dealing with irregular fractures, two ideas were considered for selecting the feature points: (1) selecting a feature point closer to the edge based on the connectivity domain, or (2) using the medial-axis transformation or refinement algorithm to extract the skeleton of the fractures and then extracting the feature points of the skeleton. With the new workflow, which used the least squares method, a fitted 3-D fracture morphology was obtained with the dip angle of 64.87° and dip direction of 35.73° for sample 1 (Figure 11c).

Macro-Fractures in Full-Hole Digital Cores
We conducted a comparative analysis of the accuracy between the new method and the original method for extracting the apparent fracture dip angle and direction. In accordance with the methodology, we utilized the original workflow to extract the fracture

Macro-Fractures in Full-Hole Digital Cores
We conducted a comparative analysis of the accuracy between the new method and the original method for extracting the apparent fracture dip angle and direction. In accordance with the methodology, we utilized the original workflow to extract the fracture dip angle and direction. However, due to some of the core being reconstructed after imaging scanning, the resulting 3-D digital data became distorted (see Figure 10b), requiring adjustments for further analysis. The extraction of fracture space distribution characteristics was then necessary. In this study, we consistently employed the threshold segmentation method, which is commonly used. This method is based on grayscale values; however, the grayscale values of fractures and edges are quite similar. Therefore, the initial step involved cropping out redundant parts of the data (see Figure 12a). Subsequently, we obtained the fracture space distribution through threshold segmentation, as shown in Figure 12b, where the blue region indicates fractures. Finally, we calculated the dip angle and direction using the MIM (see Table 2).
Processes 2023, 11, x FOR PEER REVIEW 11 of 17 dip angle and direction. However, due to some of the core being reconstructed after imaging scanning, the resulting 3-D digital data became distorted (see Figure 10b), requiring adjustments for further analysis. The extraction of fracture space distribution characteristics was then necessary. In this study, we consistently employed the threshold segmentation method, which is commonly used. This method is based on grayscale values; however, the grayscale values of fractures and edges are quite similar. Therefore, the initial step involved cropping out redundant parts of the data (see Figure 12a). Subsequently, we obtained the fracture space distribution through threshold segmentation, as shown in Figure 12b, where the blue region indicates fractures. Finally, we calculated the dip angle and direction using the MIM (see Table 2).
(a) (b)    To minimize errors and ensure accuracy, the new workflow requires the selection of XZ and YZ planes of cropped core data for extracting fracture feature point coordinates. Firstly, the CT image undergoes pre-processing, resulting in effective orthogonal section 2-D images, as shown in Figure 13a. Each fracture, represented as fractures 1-4, can be observed from top to bottom. Next, the 2-D orthogonal core image undergoes image binarization, as depicted in Figure 13c,d. Following this, the fracture feature point coordinates are automatically extracted using the minimum external algorithm. Lastly, the proposed method is employed to fit the 3-D shape of the fracture and extract the apparent dip angle and direction. Figure 13e-h present the fitted planes of fractures 1-4, respectively. To minimize errors and ensure accuracy, the new workflow requires the selection of XZ and YZ planes of cropped core data for extracting fracture feature point coordinates. Firstly, the CT image undergoes pre-processing, resulting in effective orthogonal section 2-D images, as shown in Figure 13a. Each fracture, represented as fractures 1-4, can be observed from top to bottom. Next, the 2-D orthogonal core image undergoes image binarization, as depicted in Figure 13c,d. Following this, the fracture feature point coordinates are automatically extracted using the minimum external algorithm. Lastly, the proposed method is employed to fit the 3-D shape of the fracture and extract the apparent dip angle and direction. Figure 13e-h present the fitted planes of fractures 1-4, respectively. We processed seven sets of data using both the original workflow and the new workflow. The results are presented in Table 2, where, for both dip angle and dip direction, 1 and 2 refer to the results calculated using the original workflow and the new workflow, respectively. We processed seven sets of data using both the original workflow and the new workflow. The results are presented in Table 2, where, for both dip angle and dip direction, 1 and 2 refer to the results calculated using the original workflow and the new workflow, respectively.

Discussion
For the study of practical fractures, we selected one micro-CT digital core and six 3-D full-hole digital cores. We applied the MIM and the least squares method to extract the fracture dip angle and direction, allowing for a comparative analysis. Figure 14 shows the comparison between the results of the original method and the new method. In Figure 14, the color of the dots reflects the dip angles, and the horizontal and the vertical axis represent the results calculated using the MIM and the least squares method, respectively. The comparison of dip angles is presented in Figure 14a. Throughout all of the specimens, the calculated dip angles fall within a range of 3 • of error. Furthermore, the largest absolute error value among these fractures, as indicated in Table 2, is 2.15 • . The average absolute error and relative error are 0.69 and 17.79%, respectively. In Figure 14b, the results analysis of the dip directions is presented. Throughout all of the specimens, when the fracture dip angle is less than 5 • , the error in the calculated value of the dip direction will be larger, such as fracture 1 of sample 2, for which the absolute error is 145.88 • while the dip angle 1 is 1.16 • , and fracture 4 of sample 6, for which the absolute error is −77.06 • while the dip angle 1 is 0.89 • . This is because when the dip angle is small, which means that the fracture is almost horizontal, the dip direction can change easily when the dip angle changes a little. But for the near-horizontal fractures, the dip direction is not important. They are always considered as horizontal fractures. These comparison results demonstrate the accuracy of the method proposed in this paper for extracting the apparent attitude of fractures in practical cases. The original workflow, using the MIM, is capable of accurately determining the dip angle and direction for fractures in 3-D core data. However, it necessitates processing the entire 3-D digital core data. This process involves several challenges. Firstly, the CT images often contain artifacts, noise, and other sources of interference, which necessitates cropping the data prior to subsequent processing. Secondly, achieving accurate threshold segmentation for the entire 3-D data is difficult, and inaccuracies in this step can impact the precision of the parameter extraction results. Lastly, the MIM method for extracting fracture parameters requires skilled researchers to manually carry out each step of the procedure [44,45]. Consequently, the original workflow encounters issues related to high memory consumption, complex processing procedures, and slow processing speeds. In contrast, the new workflow only requires the extraction of fracture location information from 2-D image slices, eliminating the need for cropping and threshold segmentation of the entire 3-D data. Consequently, the new workflow significantly improves efficiency while maintaining an acceptable accuracy of parameter extraction.
Although we have developed a new, efficient workflow for extracting fracture dip The original workflow, using the MIM, is capable of accurately determining the dip angle and direction for fractures in 3-D core data. However, it necessitates processing the entire 3-D digital core data. This process involves several challenges. Firstly, the CT images often contain artifacts, noise, and other sources of interference, which necessitates cropping the data prior to subsequent processing. Secondly, achieving accurate threshold segmentation for the entire 3-D data is difficult, and inaccuracies in this step can impact the precision of the parameter extraction results. Lastly, the MIM method for extracting fracture parameters requires skilled researchers to manually carry out each step of the procedure [44,45]. Consequently, the original workflow encounters issues related to high memory consumption, complex processing procedures, and slow processing speeds. In contrast, the new workflow only requires the extraction of fracture location information from 2-D image slices, eliminating the need for cropping and threshold segmentation of the entire 3-D data. Consequently, the new workflow significantly improves efficiency while maintaining an acceptable accuracy of parameter extraction.
Although we have developed a new, efficient workflow for extracting fracture dip angles and directions from full-hole digital cores, there are limitations associated with solely relying on core data for information extraction. Once cores are removed from the wellbore, the true location and orientation information is lost. Consequently, the dip angles and directions we obtain are only apparent angles and directions and do not reflect the original attitudes of the fractures. To address this issue, the most effective approach is to combine the data from core analysis with micro-resistivity imaging logging data, which are commonly used to calculate fracture attitudes (dip angle and dip direction). By utilizing imaging logging tools, we can obtain a resistivity map of the wellbore wall. Plate fractures in this map manifest as sinusoidal curves in the image logs. With fracture extraction and borehole deviation correction, we can derive the true orientations of the fractures. By carefully comparing the fractures observed in the cores and the imaging logging maps, we can spatially reconstruct the cores [46,47]. Micro-resistivity imaging logging enables us to determine the true attitudes and morphology of the fractures on the wellbore wall. However, it does not allow for direct visualization of the fracture faces, and sometimes it can be challenging to identify fractures due to the influence of bedding planes. Additionally, distinguishing an open fracture from a shale-filled fracture based solely on the imaging logs can be difficult. In contrast, digital core techniques enable us to directly identify open fractures, examine the morphology of the fracture faces, and obtain spatial distributions of the fractures. By combining these two methods, we can achieve a more accurate characterization of fractures.
Moreover, further research should be conducted to address the challenges associated with irregular fractures. This research would enhance the accuracy of calculating fracture dip angles and directions for fractures with diverse morphologies and increase the applicability of the method.

Conclusions
In this paper, we propose a novel and efficient automatic workflow for extracting the apparent dip angle and direction of fractures in 3-D digital core images using the least squares method. The workflow consists of several steps. Firstly, we extract the coordinates of fractures feature points from the edges of two orthogonal planes, namely, the XZ and YZ surfaces. Next, we automatically extract all the fracture coordinates within the core image. Then, we fit a plane to the fractures using the least squares method. Finally, we determine the apparent dip angle and direction of the fracture by calculating the angle between the normal vector of the fitting plane and the positive direction of the Z-axis and obtain the fractures' apparent dip direction by calculating the angle between the projection of the normal vector on plane XY and the X-axis.
The method presented in this paper enables accurate determination of the fracture apparent dip angle and direction. The error in the apparent dip angle and direction is less than 1 • when extracting the plate fracture model with known angles. Moreover, when applied to micro-CT and 3-D full-hole digital core images, the difference in calculating the dip angle between the original workflow and the new workflow is less than 3 • , and the dip direction extraction results between the original workflow and the new workflow are more accurate when the fracture has a higher dip angle. Notably, this new workflow not only accurately extracts the apparent dip angle and direction of fractures but also facilitates automatic data processing, significantly enhancing efficiency.
The dip angles and directions of the fractures obtained in this study are considered apparent angles and directions rather than true dip angles and directions because the spatial locations of the cores are not accurately determined. However, by combining fullhole digital cores with imaging logging data, it is possible to restore the cores to their true locations and accurately characterize the true orientations of the fractures. Therefore, it is