Received: 11 January 2019 Revised: 11 July 2019 Accepted: 18 July 2019

DOI: 10.1002/nag.2993

RESEARCH ARTICLE

Digital measurement of 2D and 3D cracks in sandstones

through improved pseudo color image enhancement and 3D
reconstruction method

Zhi Zhao | Xiao‐Ping Zhou

School of Civil Engineering, Chongqing

University, Chongqing 400045, China
Summary
An integrated approach of improved pseudo color image enhancement (IPCE)
Correspondence method and three‐dimensional (3D) reconstruction technology is proposed.
Xiao‐Ping Zhou, School of Civil
Engineering, Chongqing University, The evolution characteristics of two‐dimensional (2D) and 3D spatial cracks
Chongqing 400045, China. in X‐ray computed tomography (CT) images of sandstones subjected to triaxial
Email: [email protected]
compression with confining pressure of 10, 20, and 30 MPa are investigated.
Funding information The equations of crack width, length, and dip angle are established based on
National Natural Science Foundation of the proposed digital damage ratio. Based on the pseudo color images, point
China, Grant/Award Numbers: 51679017
cloud maps with roughness and spatial images of cracks are investigated. The
and 51839009; Natural Basic Research
program 973 of China, Grant/Award numerical results show that this proposed method is effective to study and
Number: 2014CB046903 understand the fracturing properties of rocks.

KEYWORDS
3D reconstruction, digital damage ratio, improved pseudo color image enhancement algorithm,
triaxial compression test, X‐ray CT images

1 | INTRODUCTION

Cracks are formed under the natural complex tectonic movements or the artificial engineering activities, which can be
classified into microscopic cracks, mesoscopic cracks and macroscopic cracks.1,2 These cracks play a crucial role on the
structure integrity, the strength, damage mechanism, and failure behaviors of geomaterials. Therefore, it is of great impor-
tance to understand evolution behaviors of cracks. Recently, digital image analysis techniques (DIATs) based on X‐ray
computed tomography (CT) imaging are introduced to geotechnical engineering,3,4 which provide an effective way to bet-
ter understand the cracking behaviors and failure process of geomaterials under various stress loading conditions.5-8 For
example, Sun et al explored the patterns of porosity evolution during the bearing failure of cement paste backfill subjected
to uniaxial compression,9 Takano et al investigated the deformation and the failure of a soil under triaxial compression
using in‐situ X‐ray CT images and digital image correlation.10 In addition, DIATs based on X‐ray CT images can also be
applied to investigate the petrophysical parameters of porous geomaterials such as porosity and pore size. For example,
Mayo et al quantitatively analyzed the porosity characteristics of carbonates, shales, and coal using synchrotron CT scan-
ning technology.11 Houston et al investigated the pore size distribution of soils by using tomographic and synthetic 3D
images.12 Tian and Han investigated the evolution and spatial distribution of internal pore structures of concrete by using
CT scanning. Moreover, X‐ray CT techniques are applied to investigate the characteristics of cracks as well.13 Bai and Liu
quantitatively measured the growth of crack through the edge detection of CT image of coal rocks by wavelet transforma-
tion.14 Mostafavi et al simulated the initiation and damage process of cracks in the poly‐granular graphite using X‐ray CT
images.15 Based on damage mechanics theory and X‐ray CT tests, Ma and Chen proposed a coupled approach to
Int J Numer Anal Methods Geomech. 2019;1–20. wileyonlinelibrary.com/journal/nag © 2019 John Wiley & Sons, Ltd. 1
2 ZHAO AND ZHOU

quantitatively describe the crack characteristics of shale after hydration.16 Tung et al conducted a large‐scale simulation of
microcrack initiation and propagation using a continuum mechanics description in realistic, voxel‐based microstructures
of heterogeneous quasi‐brittle materials obtained from CT imaging techniques.17 Different from the DIATs, numerical
techniques are another way to study the cracking behaviors in geomaterials. For instances, Rabczuk et al proposed the
cracking‐particle method with meshfree techniques to investigate the cracking behaviors of materials, in which there do
not exist any “mesh” orientation bias, and their results showed excellent consistency with the experiment results.18,19
Huynh et al presented an effective computational approach based on polygonal extended finite element method (XFEM)
to analyze the crack propagation with complex geometry, in which the single and mixed‐mode stress intensity factors for
cracks are analyzed as well as the computational efficiency.20 Ren et al developed the dual‐horizon peridynamics (DH‐PD)
approach to study the crack propagation in heterogeneous materials. The DH‐PD enriches the peridynamics theory and
can be considered as a landmark, which leads to successful and accurate simulations with irregular material.21 Zhou
et al studied the fracture propagation in poroelastic medias based on the phase field modeling methods, and the numerical
results are in good agreement with the previous experimental and numerical results.22,23 However, although the previous
work on the characteristics, the cracking behaviors, and failure process of geomaterials are fruitful, they are still limited to
the two‐dimensional (2D) cracks and three‐dimensional (3D) surface cracks. Moreover, the numerical techniques are hard
to capture the isolated secondary cracks away from the primary cracks and to recognize the potential crack zones
(weak zones). In addition, the studies of the characteristics, cracking behaviors, and evolution process of 3D spatial cracks
are still in its infancy phase. 3D X‐ray CT imaging technology is rarely applied to investigate the features, cracking behav-
iors, and failure process of the 3D cracks in geomaterials.
Generally, with the development of imaging technologies, many imaging apparatuses (eg, X‐ray CT) can offer high‐
resolution image of geomaterials. However, some characteristics are inevitably hidden by the pixel or voxel with low con-
trast. Therefore, various algorithms are developed to enhance the image including Histogram equalization algorithm,24
Retinex method,25 Wavelet transform algorithm,26 Partial differential equations algorithm,27 and pseudo color image algo-
rithm.28 For examples, Liange et al proposed an adaptive contrast enhancement algorithm to enhance infrared images
based on double plateaus histogram equalization,29 Al‐Ameen and Sulong improved the CT image using tuned brightness
controlled single‐scale Retinex method,30 Sahu proposed a novel method to enhance infrared image, in which an additive
wavelet transform is used for the enhancement and filtration of homomorphic infrared images,31 Pu et al proposed a novel
conceptual formulation of the fractional‐order variational framework to enhance the multi‐scale nonlocal contrast of
images with texture preserved,32 Jiang et al proposed an adaptive pseudo color enhancement method for welding radio-
graphic images based on the hue, saturation, and intensity (HSI) color space and the self‐transformation of pixels.33 How-
ever, although these classical image enhancement algorithms can enhance the image quality by adjusting the contrast,
brightness, and frequency of images, they are limited to the enhancement of the local areas and relatively single image for-
mats (gray images or color images). Meanwhile, although the classical pseudo color enhancement method has the ability to
make the images vary from gray space to color space, it is limited to the relatively unitary color grades to distinguish similar
objects with adjacent gray levels in X‐ray CT images, which has not the ability to separate the complex micro‐
characteristics (eg, primary and secondary cracks) with similar gray levels in geomaterials. In other words, few of studies
on the characteristics of cracks are reported using the pseudo image enhancement algorithm based on X‐ray CT images.
In this paper, an improved pseudo color image enhancement (IPCE) algorithm is proposed to accentuate the features
of cracks hidden with low contrast in the rock X‐ray CT images. Moreover, based on the proposed IPCE algorithm, the
damage ratio, crack width, crack length, and crack dip angle from the X‐ray CT images are defined to quantitatively
describe the characteristics of 2D cracks. In addition, the characteristics and evolution process of 3D spatial cracks
are investigated using the 3D reconstruction technique.
This paper is organized as follows: The material and experiment methodology are described in Section 2. The integrated
approach of IPCE algorithm and 3D reconstruction technology are described in Section 3. The evolution laws of fracture
parameters are described in Section 4. The results and discussions are summarized in Section 5; brief summaries and con-
clusions are drawn in Section 6.

2 | MATERIAL A ND EXPERIMENT METHODOLOGY

2.1 | Specimen descriptions

The sandstone specimens are all obtained from the Mazongxiang region of Chongqing, Southwest China. The mine
region is consisted of Luoguanshan, Xishan, and Xindainzi anticlines with an average dip angle of 17°. X‐ray diffraction,
ZHAO AND ZHOU 3

microscope investigation, and automated Pemeameter‐Porosimeter analyses are conducted to obtain the mineral com-
ponents and common physical parameters. The experimental results (Figure 1) show that sandstone specimens are fea-
tured by quartz, rock debris, feldspar, and clay mineral with a relative content of 65.58%, 20.2%, 8.16%, and 6.06%,
respectively. The average porosity, permeability, unit weight, and specific surface area are, respectively, 27.256,
0.05 md, 26.27, and 0.421, as shown in Table 1.

2.2 | Triaxial compression test

The obtained sandstone specimens are coded as specimens 1 to 3 and are polished as the standard size (Φ50 × 100 mm).
A multi‐function triaxial rock testing system (Rock 600‐50 HT PLUS) is applied to investigate the failure process of sand-
stone specimens. The sandstone specimen is loaded by an initial small axial stress and confining pressure to fix the sam-
ple location. The sandstone specimens 1 to 3 are respectively tested with a constant average axial loading speed of
0.22 MPa/min under confining pressure of 10, 20, and 30 MPa, and the complete stress‐strain curves are obtained.
Figure 2 shows the sandstone specimens 1 to 3 before and after the triaxial compression tests.

2.3 | X‐ray CT measurements

Sandstone specimens 1 to 3 are scanned by the SIEMENS‐SOMATOM scope X‐ray CT scanner, which is a high‐
resolution imaging system with a spatial resolution of 0.35 × 0.35 mm, scanning time of 1500 ms for per slice, and a
scanning voltage of 130 kV and electricity of 81 μA. The sandstone specimens are scanned after the trixial compression
tests under confining pressure of 10, 20, and 30 MPa. Three types of X‐ray CT images of sandstones are obtained from
top to bottom with a slice interval thickness of 0.75 mm. In this paper, 158 × 3 X‐ray CT images are obtained with a
pixel size of 512 × 512 for each slice. 139 × 3 and 100 × 3 X‐ray CT images are selected and cropped as 141 × 141 pixels
to study the evolution process of 2D and 3D cracks, respectively.

2.4 | Experiment results

Triaxial compression and X‐ray CT scanning testing results of sandstone specimens 1 to 3 are depicted in this subsec-
tion. Specimens 1 to 3 are, respectively, failed under the axial stress of 157.585, 206.355, and 222.121 MPa, and the fail-
ure modes exhibit shear failure with a slight difference between the failure surfaces. In addition, to interpret the failure

FIGURE 1 Mineral components of sandstone: A, relative contents of mineral; B, microscope image of sandstone under polarized light
(200 × 50 μm) [Colour figure can be viewed at wileyonlinelibrary.com]

TABLE 1 Physical parameters of sandstone obtained from experiments

Size, mm γ, kN/m3 ϕ, % K, md SSA, m2/g

Φ50 × 100 26.27 27.256 0.05 0.421

Note. γ—unit weight; ϕ—porosity; K—average permeability; SSA—specific surface area.

4 ZHAO AND ZHOU

FIGURE 2 The intact and damaged

sandstone specimens 1 to 3 [Colour figure
can be viewed at wileyonlinelibrary.com]

process in detail, the stress‐strain curve of typical X‐ray CT images from sandstone specimen 1 is plotted in Figure 3.
Similar to the initiation, propagation, and coalescence of cracks in rock specimens with pre‐existing cracks,34 the failure
processes of the sandstone samples are described as follows: The specimen is firstly compacted with increasing axial
stress, as shown in Figure 3 (1 and 2). Next, microcracks initiate and propagate in a way of symmetric double shell
cracks, as shown in Figure 3 (3 and 4). Then, microcracks coalesce to form macrocracks, and macrocracks propagate
in a stable way as axial stress increases, as shown in Figure 3 (5 and 6). Fourthly, secondary cracks are initiated from
the tips of macrocracks, and macrocracks propagate in an unstable way when axial stress reaches the peak value
(157.585 MPa), as shown in Figure 3 (7 and 8). Finally, macrocracks propagate to the boundary of the specimen, leading
to the complete failure of the sandstone specimen, as shown in Figure 3 (9 and 10).

3 | T H EO R Y O F IP C E A N D 3 D R E C O N S T R U C T I O N A LG O R ITH M

3.1 | Image representation and color space

Any digital X‐ray CT image is represented as a 2D array with each single pixel (image element) considered as an ele-
ment of array. Namely, a rock X‐ray CT image with size of N × N can be written as
ZHAO AND ZHOU 5

FIGURE 3 Stress‐strain curve of sandstone specimen 1 with corresponding X‐ray computed tomography (CT) images at different stress
levels during the failure process subjected to triaxial compression [Colour figure can be viewed at wileyonlinelibrary.com]

f ð1; 1Þ f ð1; 2Þ ⋯ f ð1; N Þ




N − 1N − 1
f ð2; 1Þ f ð2; 2Þ ⋯ f ð2; N Þ




I Rock ¼ ∑ ∑ f ðx; yÞ ¼

; (1)

⋮ ⋮ ⋮

0 0






f ðN; 1Þ f ðN; 2Þ ⋯ f ðN; N Þ

where f (x,y) denotes the grayscale intensity value of image.

Therefore, the components in X‐ray CT images of rocks (minerals, pores, and cracks etc) are represented by the pixel
sub‐arrays of gray scales. Initially, the pixel arrays are stored in Digital Imaging and Communications in Medicine
(DICOM) images, which can be converted to the common image formats (BMP, JPG, and TIFF etc) for digital image
processing in computer. In this paper, X‐ray CT images are converted to BMP images.
Generally, the gray level shades of BMP image vary from the weakest intensity of 0 (black) to the strongest intensity
of 255 (white). However, the recognition range of gray value is limited to a few dozen by human vision,35 and one of the
major drawbacks associated with the gray transformation is the ignorance of tiny contents in gray intensity image. For-
tunately, human beings can recognize and distinguish thousands of colors. Thus, the recognition of different objects can
be improved by the utilization of colors. In other words, the components of rock X‐ray CT image can be well identified
by color image processing technique.
The colors in a color image are usually presented in the RGB color space, and the gray intensity value in a gray image
is presented in the gray level space. Figure 4A to C, respectively, show the RGB color space, the discriminability of
human vision to color differences, and the gray level space.
Usually, X‐ray CT images are presented in the gray level space with different components represented by different
gay intensity values. Although the classical image enhancement algorithms can enhance the image quality by adjusting
the frequency of image, they are limited to the enhancement of the local areas, which has not the ability to separate the
complex micro‐characteristics (eg, primary and secondary cracks) with similar gray levels in geomaterials. Thus, the
IPCE algorithm is proposed to well describe the crack characteristics in X‐ray CT images of rocks. The characteristics
of interesting regions can be enhanced and distinguished through the transformations between the gray intensity
images and color images by IPCE algorithm.
6 ZHAO AND ZHOU

FIGURE 4 The transformation between the RGB color and the gray level space, A, the RGB color space; B, the relationship between the color
differences and the disciminability of human vison to colors; C, the gray level space [Colour figure can be viewed at wileyonlinelibrary.com]

3.2 | Conceive of IPCE algorithm

In the IPCE approach, the crack segmentation threshold TC is firstly determined. Assuming the unique gray level of the
image is divided into two parts G0 → TC and GTC → 255 by TC. Thus, the crack segmentation threshold value is calculated
by

8
< G0→TC ¼ ff ðx; yÞjf ðx; yÞ ≤ TC g
>
GTc→255 ¼ ff ðx; yÞjf ðx; yÞ ≤ TC g ; (2)
>
:
TC ¼ ∐max ½jΓðG0→TC Þ − ΓðGTc→255 Þj

where Γ and ∐max are the variance function and the searching function of maximum, respectively.
Then, the gray levels of the image are divided into M gray level subsets ( f 1, f 2,⋯ f i⋯, f M) by the minimum, average
gray level, and crack segmentation threshold to construct the basic color transformation functions IR(x,y), IG(x,y), and
IB(x,y). Finally, a pseudo color image is reproduced by combining the outputs of three channel transformations
consisted of the red, green and blue components, as shown in Figure 5.
The construction of these three channel transformation functions is crucial in the IPCE approach. The transforma-
tion function is considered as a linear function due to the different gray values of components in a rock X‐ray CT image.
In this paper, the transfer functions are considered as piecewise linear function, which can be written as

FIGURE 5 Schematic maps of the result of improved pseudo color image enhancement (IPCE) approach, A, the gray image with gray
intensity; B, the pseudo color image with gray intensity [Colour figure can be viewed at wileyonlinelibrary.com]
ZHAO AND ZHOU 7

8
>
> 0 f ≤ TC
>
< f max − f min
IRðx; yÞ ¼ 255 − ×f TC < f ≤ f a ; (3)
>
> f a − 90
>
:
255 f > fa

8
>
> 0 f ≤ f min
>
>
>
> 90 − f f min < f ≤ TC
>
>
<
f − f min
IGðx; yÞ ¼ f þ max ×f TC < f ≤ f a ; (4)
>
> f − 90
>
>
a i
>
> −
>
> f f
: 255 − f þ
max min
×f f > fa
f a − 90 i

8
> 0 f ¼0
>
>
>
> 4f − 1 1 ≤ f ≤ f min
>
>
>
>
>
< 255 f min < f ≤ TC
IBðx; yÞ ¼ f max − f min ; (5)
>
> 255 − ×f TC < f ≤ f a
>
> f a − 90
>
> i
>
>
>
> f max − f min
: 4f − ×f f > fa
f a − 90 i

where f max, f min, and f a are, respectively, the maximum, minimum, and average value of gray intensity value except
the gray level subset of 0, f i is the average gray intensity in the ith level.
Figure 6 shows the flow chart of converting a gray image to a pseudo color image. The pseudo color images can be
reproduced by the combination of the red, green, and blue components, which are obtained by Equations (2) to (5).
Meanwhile, considering the region of interest (RIO) scale of 141 × 141 pixels, X‐ray CT gray images of sandstone spec-
imens 1 to 3 are converted to the pseudo color images to accentuate the induced cracks and sub‐cracks components, as
shown in Figure 7. Each element of the gray image array is colored by the IPCE algorithm, which is presented by dif-
ferent colors in the corresponding location. Figure 7A to C represents the original images of the sandstone specimens 1
to 3. Figure 7D to F shows the corresponding pseudo color images. Figure 7G to I shows the color images for Figure 7A
converted by the classical green color, the classical rainbow code and the classical hot metal scale code enhanced algo-
rithms, respectively. Obviously, the proposed IPCE algorithm shows better performance than the classical algorithms.
As shown in Figure 7D to F, the induced cracks by stress and underlying sub‐cracks are enhanced, in which the induced
cracks and sub‐cracks are colored by condensed black and blue, respectively. Meanwhile, the underlying propagation
directions of cracks are presented by scatter blue pixel. In addition, the evolution process of cracks can be defined
through the characteristics of pseudo color images.

FIGURE 6 The flow chart of converting gray image to pseudo color image
8 ZHAO AND ZHOU

FIGURE 7 The pseudo color images converted by the gray X‐ray computed tomography (CT) images of sandstone specimens 1 to 3 using
improved pseudo color image enhancement (IPCE), A to C, the original images; D to F, the pseudo color images converted by the proposed
IPCE; and G to I, the color images of Figure 7D from sandstone specimen 1 converted by the classical green color, the classical rainbow code,
and the classical hot metal scale code enhanced algorithms, respectively

3.3 | Estimation of damage ratio and porosity

Once the cracks and different components in X‐ray CT images are colored in different colors, the characteristics of the
cracks can be described. To obtain the quantitative characteristics of cracks, the digital damage ratio and porosity are
established using digital image processing principles.

3.3.1 | Definition of the damage ration and porosity

Damage ratio is one of the most important indexes to describe the damage characteristics in fracture mechanics. In gen-
eral, any of the void space ratio, porosity, and the damage variable can be considered as the percentage of the space of
corresponding structural components occupied in the whole specimen scale. In this paper, the digital porosity (P) and
digital damage ratio (D) are evaluated by the pixel numbers percentage of the pore space and the damage region (cracks)
in a N × N X‐ray CT image, which are, respectively, written as
ZHAO AND ZHOU 9

∑M ½f ðx; yÞ
Rp
P¼ ; (6)
N2

∑ M ½f ðx; yÞ
RC
D¼ ; (7)
N2

where M is the counters of the pixel in the damage region or the pore space (void space), Rc and Rp are pixel sets of the
damage region and the pore space, respectively.

3.3.2 | Determination of the damage region and pore space

Once the counters of the pixels in the damage region and pore space (void space) are determined, damage ratio and
porosity can be calculated by Equations (6) and (7), which is extremely affected by the estimated numbers of the pixels
in the void space. Therefore, the determination of the void space is crucial, which can be obtained using the image seg-
mentation method. Usually, the image segmentation method is used to partition an X‐ray CT image into two or multi-
ple segments.36 In fact, various image segmentation algorithms including thresholding,37 clustering,38 deformable
models,39 and atlas‐based algorithm40 are proposed to detach the connected objects. In this study, to obtained accurate
void space, the proposed IPCE algorithm is used to capture the damage region (Figure 8A), and the C‐V model is
applied to separate the pore space from the solid objects.41
When the damage region and pore space are, respectively, determined by the IPCE algorithm and C‐V model, the
damage ratio and porosity of sandstone specimens 1 to 3 can be obtained by Equations (6) and (7).

3.4 | Evaluation of crack parameters

Except for the damage ratio and the porosity, the crack width, crack length, and crack dip angle are defined as well.
Figure 8 schematically illustrates the determination process of the quantitative parameters to describe the crack char-
acteristics. The original gray intensity image is converted to pseudo color image by the proposed IPCE algorithm, as
shown in Figure 8A. An interesting area is selected to find the optimal ellipse to determine the crack parameters in
the damage region, as shown in Figure 8B,C. It is worth to notice that the optimal ellipse must contain pixels as much
as possible from the damage region, which must stratify the following relationship as much as possible

∑ M ½f ðx; yÞ
RE
Pd ¼ ≈ 1; (8)
∑ M ½f ðx; yÞ
RC

where RE is the set of the pixels in the ellipse, Pd is the pixel percentage of the ellipse occupied in the damage region.
During this process shown in Figure 8C, the centroid (Oc) of the selected damage region is firstly calculated to locate
the center of the ellipse, which is written as
8
>
> X c ¼ ∑f ðx; yÞ × x i =∑f ðx; yÞ
>
> Rc Rc
<
Y c ¼ ∑f ðx; yÞ × yj =∑f ðx; yÞ ; (9)
>
> Rc Rc
>
> :
:
Z c ¼ zk

where Xc, Yc, and Zc are the centroid coordinates, xi, yj, and zk are the coordinates of the pixel in the damage region.
Then, the major axis (AB) of the ellipse is determined by the maximal value of the distance between the centroid
pixel and the pixels in the damage region. Thus, the coordinate of the pixel along the semi major axis direction with
the maximal distance to the centroid pixel is written as
10 ZHAO AND ZHOU

FIGURE 8 An example of the crack determination [Colour figure can be viewed at wileyonlinelibrary.com]

8 " rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!#
>    2
>
> 2
< ½X a ; Y a  ¼ x i; yj j max ðX c −x i Þ þ Y c −yj
; (10)
>
> qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
>
: Lc ¼ 2 ðX −X Þ2 þ ðY −Y Þ2
c a c a

where Xa and Ya are the coordinates of the pixel along the semi major axis direction with the maximal distance to the
centroid pixel, Lc is the length of cracks.
Next, the minor axis (CD) is determined by the maximal value along with the vertical line of the semi‐major axis of
the ellipse in the damage region. Thus, the coordinate of the pixel along the semi minor axis direction with the maximal
distance to the centroid pixel can be written as
ZHAO AND ZHOU 11

8 " rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !#
>    2
>
> 2
>
> ½X b ; Y b  ¼ x i; yj j max ðX c −x i Þ þ Y c −yj
>
<
! ! ; (11)
>
> Oc A ·Oc B ¼ ðX c − X a Þ × ðX c − X b Þ þ ðY c − Y a Þ × ðY c − Y b Þ ¼ 0
>
> qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
>
>
: Wc ¼ 2 ðX −X Þ2 þ ðY −Y Þ2 Þ
c b c b

where Xb and Yb are the coordinates of the pixel along the semi‐minor axis direction with the maximal distance to the
! !
centroid pixel. Oc A and Oc B are, respectively, the vectors of the semi‐major axis and semi‐minor axis, Wc is the width of
cracks.
Finally, the angle between the semi‐major axis and the X‐axis can be written as
8  
< α ¼ arctan Y c − Y a
Xc − Xa ; (12)
:
Ac ¼ α

where α is the angle between the semi‐major axis and the X‐axis, Ac is the angle of cracks.

3.5 | 3D reconstruction of the spatial cracks

3D reconstruction techniques are effective tools to investigate the microstructures of rocks with the aid of X‐ray CT
images. In this paper, splatting reconstruction method,42 with the advantages of high speed for pixel information search
and reconstruction of 3D models, low memory demands, and progressive display of 3D reconstruction models, is mod-
ified to investigate the characteristics and evolution process of 3D spatial cracks based on X‐ray CT images.
When all the cracks are captured and accentuated by IPCE algorithm, the crack regions are extracted from each
pseudo X‐ray CT image. Let χ(V) denote the discrete voxel set of cracks, and Vi,Vj,Vk represent the weight coordinates
Lc3 ¼ −1:98222E11 × Dc3 6 þ 4:18249E10 × Dc3 5 − 3:16703E9 × Dc3 4
of digital voxel V in the þ 1:01701E8 × Dc3 3 − 1:24761E6 × Dc3 2 directions, respectively.
þ 4662:14232 × Dc3 − 0:66368
To conduct the reconstruction procedure, the discrete voxel set of cracks are firstly converted to vector voxel set by

ΦV ðx; y; zÞ ¼ ∑ ΩV x − V i ; y − V j ; z − V z χ ðV Þ; (13)
V ∈N

where ΦV(x,y,z) is vector voxel set of cracks, N is the subset of voxel in the non‐zero volume reconstruction kernel ΩV.
In this paper, the non‐zero volume reconstruction kernel is modified for the reconstruction of cracks, which is
expressed as
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
χ 1 ðV Þ−1 =2 þ χ θ ðV Þ−1 =2
ΩV ½χ 1 ðV Þ; χ θ ðV Þ ¼ ∑ ; (14)
θ χ 1 ðV Þ1=4 × χ θ ðV Þ1=4

where χθ(V) denotes the neighbor voxel towards the direction of θ = 0,45,90,135 in each XOY, XOZ, and YOZ plane
surrounded the converted voxel χ1(V).
Once all the voxel sets of cracks are crossed through by a ray and projected on the digital image plane, the projection
of this digital image can be written as

þ∞
N2

∏RV ¼ ∑ χ ðV n Þ ∫ ΩV x − V nx ; y − V ny ; z − V nz dR; (15)

n¼1 −∞

where ∏RV denotes the collected discrete data set of 3D digital voxel.
Thus, the 3D models of spatial cracks can be reconstructed using the modified reconstruction method and kernel
function, and the corresponding 3D point data set can be extracted by Equation (15).
12 ZHAO AND ZHOU

Figure 9A to C shows the real failure surface of sandstone specimens 1 to 3. Figure 9D to F shows the corresponding
reconstructed 3D spatial cracks. Apparently, the reconstructed 3D spatial cracks are all shear cracks. Figure 9D shows
the cross cracks in the 3D space under confining pressure of 10 MPa, which are presented as a shell cracks located at the
two sides of the center axis line. The distance between the cracks gradually decreases from top to bottom of the speci-
men. When the cracks propagate to the bottom of the specimen, the cracks coalesce with each other, which looks like
the upper half part of X. Figure 9E,F shows the shear cracks under confining pressure of 20 and 30 MPa, respectively. It
is easy to find that the cracks plotted in Figure 9E are similar to those in Figure 9F. However, the thickness of 3D crack
in Figure 9E is larger than that of 3D crack in Figure 9F, and the angle between the sheer failure surface and horizon
surface in Figure 9E is smaller than that in Figure 9F.

4 | EV OLUT I ON PR OCE S S ON THE B A SI S OF FR A C T UR E PA R A M E T E RS

In this study, the damage ratio, porosity, crack width, crack length, and crack angle are extracted from X‐ray CT images
of sandstone specimens 1 to 3 by Equations (6) to (12), as shown in Figure 10. The failure process of cracks and the char-
acteristics of fracture parameters are summarized in Section 4.1.

4.1 | Characterization of fracture parameters

Figure 10A shows the porosity of X‐ray CT images from the sandstone specimens 1 to 3 with the average value 26.548%,
27.348%, and 27.032%, respectively. Obviously, the porosity curves agree well with each other and are corresponding to
the experimental result (27.256%), proving the accuracy of Equation (6). Figure 10B to E respectively shows the damage
ratio, crack width, crack length, and crack angle distributions for sandstone specimens 1 to 3. Based on the damage ratio
distribution, the failure process of sandstone specimens are divided into five stages.

FIGURE 9 The three‐dimensional (3D) spatial cracks of sandstone specimens 1 to 3 [Colour figure can be viewed at wileyonlinelibrary.com]
ZHAO AND ZHOU 13

FIGURE 10 The distribution of fractural parameters, A, porosity; B, damage ratio; C, crack width; D, crack length; and E, crack angle

1. The compaction of natural pre‐existing fissures (I). The natural pre‐existing fissures in sandstone specimens are
compressed and closed under axial stress. No apparent damage characteristics are observed. The damage ratio,
crack width, crack length, and crack angle are all 0.
2. The initiation and propagation of microcracks (II). Microcracks in sandstone specimens initiate, propagate with
increasing axial stress. The damage ratio, crack width, crack length, and crack angle increase gradually with the
nearly same tendency and increase rate for each specimen. However, the tendencies of crack width and length from
specimen 1 have large difference from specimens 2 and 3, which may be caused by the difference of crack patterns,
such as double shell cracks for specimen 1, oblique cracks with different angle for specimens 2 and 3, as shown in
Figures 7 and 9.
14 ZHAO AND ZHOU

3. Microcracks coalescence, the formation, and stable propagation of macrocracks (III). Microcracks coalesce and
macrocracks form and propagate in a stable way with successively increasing axial stress. The damage ratios and
crack angle for specimens 1 to 3 keep almost constant, which implies the stable propagation of macrocracks. How-
ever, the tendencies of crack width and length increase first, then fluctuate slightly, which may be caused by the
transition of microcracks to macrocracks.
4. The unstable propagation of macrocracks (IV). Macroscopic cracks gradually propagate in an unstable way with
increasing axial stress. The damage ratio and crack angle fluctuate slightly, and the crack width and crack length
decrease gradually, which may be caused by unstable propagation of macrocracks.
5. The failure of sandstone specimens (V). In the final stage of the evolution process of cracks, the damage ratio, crack
width, crack length, and crack angle fluctuate sharply due to unstable propagation of macrocracks, leading to the
complete failure of rocks.

5 | R E S U L T S AN D D I S C U S S I O N S

5.1 | Color spectrum analysis of cracks

X‐ray CT gray images of sandstone specimens 1 to 3 are converted to the pseudo color images using IPCE algorithm to
accentuate the induced cracks and sub‐cracks components. Figure 11 shows the color analysis of the induced cracks and
sub‐cracks (eg, Figure 7) represented by black and blue, respectively. Obviously, different components in the X‐ray CT
images of sandstone specimens are distinguished by color grades, which are generated by the IPCE algorithm. Figure 11

FIGURE 11 The components (color spectrum) analysis of pseudo color image converted by improved pseudo color image enhancement
(IPCE) of sandstone specimens 1 to 3, A, the distribution of total image components; B, the induced crack components (black in pseudo
color images); C, the sub‐crack and penetration components (blue in in pseudo color images)
ZHAO AND ZHOU 15

FIGURE 12 Pseudo color images of sandstone specimens 1 to 3 enhanced by the improved pseudo color image enhancement (IPCE)
algorithm in each stage of the failure process [Colour figure can be viewed at wileyonlinelibrary.com]
16 ZHAO AND ZHOU

A shows the total statistical components of the X‐ray CT images from sandstone specimens 1 to 3, which is considered
as a complete color spectrum for sandstone specimens. Figure 11B shows the statistical induced crack components col-
ored by black. Figure 11C shows the statistical perforation components (sub‐cracks) colored by blue. Thus, based on the
variation of color components, the evolution process of 2D cracks can be concluded, which is apparently corresponding
to the description in Section 4.1. In addition, the evolution process considering sub‐cracks can be interpreted more
detailed by accentuating the sub‐cracks using IPCE algorithm, which is remarkable evidence for the accuracy and reli-
ability of the proposed algorithm.

5.2 | Characterization of 2D cracks based on IPCE

Figure 12 shows the implementation result obtained by IPCE algorithm on X‐ray CT images from each stage of the evo-
lution process of sandstone specimens 1 to 3 with confining pressure of 10, 20, and 30 MPa, respectively. Obviously, the
characteristics of 2D induced cracks by stress and sub‐cracks are enhanced.
Initially, the pre‐existing defects and pores in sandstone specimens 1 to 3 are compacted and closed under axial
stress, as scattered blue pixels shown in Figure 12A to C. No apparent damage occurs.
Then, with continuous increase of axial stress, induced microcracks, such as double shell cracks for specimen 1,
oblique cracks with different inclinations for specimens 2 and 3, initiate and propagate, as shown in Figure 12D to F.
Some concentrated black pixels form the cracks, and the other blue pixels show the potential growth direction
of cracks.
Microcracks coalesce to form macrocracks, and macrocracks propagate stably with increasing axial stress, as shown
in Figure 12G to I. The macrocracks are featured by the concentrated black pixels and are surrounded by the condensed
blue pixels in the potential damage region.Next, with successive increase of axial stress, macrocracks propagate gradu-
ally in an unstable way. Moreover, sub‐cracks initiated from different locations of macrocracks in sandstone specimens
1 to 3, which are featured by some concentrated blue pixels and the other blue pixels, become gradually dispersive, as
shown in Figure 12J to L.
Finally, macrocracks propagate to the boundary of sandstone specimens, leading to the failure of sandstone, and the
blue pixels get scattered again, as shown in Figure 12M to O.

5.3 | Spatial characterization of 3D crack

Once the pseudo color images of sandstone specimens 1 to 3 are obtained, the corresponding datasets for 3D cracks can
be extracted by Equation (15). To interpret the evolution process of 3D cracks clearly and intuitively, the 3D cracks
datasets from each stage of the failure process of sandstone specimens 1 to 3 are extracted. The point cloud maps of
3D cracks with roughness are constructed to investigate the evolution process of 3D cracks.At first, the pre‐existing
defects are closed with uniformly distributed point clouds and similar roughness, which implies no apparent damage
occurrence in sandstone specimens 1 to 3, as shown in Figure 13A to C.
Then, primary induced microcracks, such as the symmetrical double shell cracks in specimen 1, oblique
cracks with the small inclination in specimen 2, and oblique cracks with the large inclination in specimen 3, are
initiated from the edge of specimen and propagate towards the other side of specimens. The roughness is nonuni-
form, it is the lowest for primary microcracks zone, followed by the adjacent potential damage zone, and it is the
largest for the adjacent undamage zone, which are respectively marked by blue and white, as shown in Figure 13
D to F.
These microcracks coalesce to form macrocracks, and macrocracks propagate in a stable way as axial stress increases.
Point cloud maps are featured by the lowest roughness in macrocrack zones and by uniform roughness outside the
factured zones, as shown in Figure 13G to I.
Next, macrocracks propagate towards the boundary of samples in an unstable way as axial stress continues to
increase. In addition, sub‐cracks initiate and propagate. In specimen 1, sub‐cracks propagate and coalesce with the
symmetrical double shell cracks. In specimen 2, sub‐cracks with a small inclination initiate from the upper edge
of macrocracks and propagate along the direction of macrocracks. In specimen 3, sub‐cracks initiate from the middle
edge of macrocracks and propagate towards the boundary of the specimen, as shown in Figure 13J to L. The rough-
ness in these figures is similar to that in Figure 13G to I, but the roughness in the potential damage zones is more
obvious.
ZHAO AND ZHOU 17

FIGURE 13 Three‐dimensional (3D) point cloud maps of cracks with roughness in sandstone specimens 1 to 3 at each stage of the failure
process [Colour figure can be viewed at wileyonlinelibrary.com]
18 ZHAO AND ZHOU

FIGURE 14 Progressive display of the spatial three‐dimensional (3D) cracks [Colour figure can be viewed at wileyonlinelibrary.com]

Finally, macrocracks propagate to the boundary of samples, leading to the complete failure of the samples. Point
cloud maps in the factured zones are incomplete with the lowest roughness, and the roughness outside the factured
zones become uniform again, as shown in Figure 13M to O.
Figure 14 shows the progressive characteristics of 3D spatial cracks in each stage. The cracks in sandstone specimen 1
are the double shell cracks, which look like the upper half part of X, as shown in Figure 14A to E. The cracks in sand-
stone specimens 2 and 3 are both oblique cracks. The cracks in sandstone specimen 2 are the single induced crack with
a sub‐crack presented in the bottom of the specimen, as shown in Figure 14F to J. The cracks in sandstone specimen 3
are the induced crack with a sub‐crack in the lower surface of crack, as shown in Figure 14K to O. Moreover, the char-
acteristics and evolution process of 2D cracks (Figure 12) and reconstructed 3D spatial cracks (Figures 13 and 14) can be
verified reciprocally. In total, although the surface characteristics of the cracks have slight difference, all the cracks are
generated by the shear failure of the specimens.

6 | SUMMARY AND CONCLUSIONS

In this paper, an improved pseudo color image algorithm is proposed to improve the characteristics of cracks hidden
with low contrast. Moreover, the digital damage ratio is defined to describe the failure process of the cracks. Based
on the damage region (cracks), the crack width, crack length, and crack angle are defined to quantitatively describe
the characteristics of cracks. On the basis of the integrated method of the IPCE algorithm and 3D reconstruction tech-
niques, the detailed evolution process of 2D and spatial 3D cracks in the sandstone specimens are investigated. The
main conclusions can be drawn as follows:
ZHAO AND ZHOU 19

1. The proposed IPCE algorithm is excellent and reliable to accentuate the crack features, and the statistical color
spectrum information can be a fine factor to describe the evolution process of cracks.
2. The pore scale and crack parameters can be estimated using IPCE algorithm, and the defined crack width, crack
length, and crack angle can be fine factors to describe the crack characteristics. The evolution process of the cracks
is divided into five stages based on the damage ratio, which is in good agreement with the various stages of
experiments.
3. Based on the pseudo color images, point cloud maps with roughness and spatial images of cracks, the characteris-
tics, and evolution process of 2D and 3D cracks can be investigated.

ACK NO WLE DGE MEN TS

This work was supported by the National Natural Science Foundation of China (grant nos. 51839009 and 51679017) and
the Natural Basic Research program 973 of China (grant no. 2014CB046903).

ORCID
Xiao‐Ping Zhou https://orcid.org/0000-0003-1551-6504

R EF E RE N C E S
1. Yi SZ. Fracture and Damage Theories and Their Applications. Beijing (in Chinese): Tsinghua University Press; 1992.
2. Yang W. Macro and Micro Fracture Mechanics. Beijing (in Chinese): National Defence Industry Press; 1995.
3. Trawiński W, Bobiński J. Tejchman J two‐dimensional simulations of concrete fracture at aggregate level with cohesive elements based
on X‐ray μct images. Eng Fract Mech. 2016;168:204‐226.
4. Lublóy É, Balázs GL, Kapitány K, Barsi Á. Ct analysis of core samples from fire‐damaged concrete structures. Mag Concrete Res.
2017;69(15):1‐9.
5. Carlson WD. Three‐dimensional imaging of earth and planetary materials. Earth Planet Sc Lett. 2006;249(3):133‐147.
6. Kerckhofs G, Schrooten J, Van Cleynenbreugel T, Lomov SV, Wevers M. Validation of x‐ray microfocus computed tomography as an
imaging tool for porous structures. Rev Sci Instrum. 2008;79(1):489.
7. Suzuki T, Shiotani T, Ohtsu M. Evaluation of cracking damage in freeze‐thawed concrete using acoustic emission and x‐ray ct image.
Construct Build Mater. 2017;136:619‐626.
8. Wang Y, Li CH, Hao J, Zhou RQ. X‐ray micro‐tomography for investigation of meso‐structural changes and crack evolution in longmaxi
formation shale during compressive deformation. J Petrol Sci Eng. 2018;164:278‐288.
9. Sun W, Hou K, Yang Z, Wen Y. X‐ray ct three‐dimensional reconstruction and discrete element analysis of the cement paste backfill pore
structure under uniaxial compression. Construct Build Mater. 2017;138:69‐78.
10. Takano D, Lenoir N, Otani J, Hall SA. Localised deformation in a wide‐grained sand under triaxial compression revealed by x‐ray tomog-
raphy and digital image correlation. Soils Found. 2015;55(4):906‐915.
11. Mayo S, Josh M, Nesterets Y, et al. Quantitative micro‐porosity characterization using synchrotron micro‐ct and xenon k‐edge subtraction
in sandstones, carbonates, shales and coal. Fuel. 2015;154:167‐173.
12. Houston AN, Otten W, Falconer R, Monga O, Baveye PC, Hapca SM. Quantification of the pore size distribution of soils: assessment of
existing software using tomographic and synthetic 3d images. Geoderma. 2017;299:73‐82.
13. Tian W, Han N. Pore characteristics (>0.1 mm) of non‐air entrained concrete destroyed by freeze‐thaw cycles based on CT scanning and
3D printing. Cold Reg Sci Technol. 2018;151:314‐322.
14. Bai C, Liu X. The edge detection technology of CT image for study the growth of rock crack. Isecs International Colloquium on Computing,
Communication, Control, and Management Proceedings. 2009;1(4):286‐288.
15. Mostafavi M, Baimpas N, Tarleton E, et al. Three‐dimensional crack observation, quantification and simulation in a quasi‐brittle material.
Acta Mater. 2013;61(16):6276‐6289.
16. Ma T, Chen P. Study of meso‐damage characteristics of shale hydration based on CT scanning technology. Petrol Explor Dev.
2014;41(2):249‐256.
17. Tung NT, Yvonnet J, Bornert M, Chateau C, Bilteryst F, Steib E. Large‐scale simulations of quasi‐brittle microcracking in realistic highly
heterogeneous microstructures obtained from micro ct imaging. Extreme Mechanics Letters. 2017;17:50‐55.
18. Rabczuk T, Zi G, Bordas S, Nguyen‐Xuan H. A simple and robust three‐dimensional cracking‐particle method without enrichment.
Comput Methods Appl Mech Eng. 2010;199(37–40):2437‐2455.
20 ZHAO AND ZHOU

19. Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int J Numer Methods Eng.
2004;61(13):2316‐2343.
20. Huynh HD, Nguyen MN, Cusatis G, Tanaka S, Bui TQ. A polygonal XFEM with new numerical integration for linear elastic fracture
mechanics. Engineering Fracture Mechanics. 2019;213:241‐263.
21. Ren H, Zhuang X, Rabczuk T. Dual‐horizon peridynamics: a stable solution to varying horizons. Comput Methods Appl Mech Eng.
2017;318:762‐782.
22. Zhou S, Rabczuk T, Zhuang X. Phase field modeling of quasi‐static and dynamic crack propagation: COMSOL implementation and case
studies. Advances in Engineering Software. 2018;122:31‐49.
23. Zhou S, Zhuang X, Rabczuk T. A phase‐field modeling approach of fracture propagation in poroelastic media. Engineering Geology.
2018;240:189‐203.
24. Kim YT. Contrast enhancement using brightness preserving bi‐histogram equalization. IEEE T Consum Electr. 1997;43(1):1‐8.
25. Jobson DJ, Rahman Z, Woodell GA. A multiscale retinex for bridging the gap between color images and the human observation of scenes.
IEEE T Image Process. 2002;6(7):965‐976.
26. Hsieh CT, Lai E, Wang YC. An effective algorithm for fingerprint image enhancement based on wavelet transform. Pattern Recogn.
2003;36(2):303‐312.
27. Yi L. Fourth‐order partial differential equations for image enhancement. Appl Math Comput. 2006;175(1):430‐440.
28. Cho MK. Study of gray image pseudo‐color processing algorithms. Proceedings of SPIE—The International Society for Optical Engineering.
2012;19(11):19.
29. Liang K, Ma Y, Xie Y, Zhou B, Wang R. A new adaptive contrast enhancement algorithm for infrared images based on double plateaus
histogram equalization. Infrared Phys Techn. 2012;55(4):309‐315.
30. Al‐Ameen Z, Sulong G. A new algorithm for improving the low contrast of computed tomography images using tuned brightness con-
trolled single‐scale retinex. Scanning. 2015;37(2):116‐125.
31. Sahu A. Infrared image enhancement using wavelet transform. Int J Eng Res Appl. 2012;2(2):40‐47.
32. Pu YF, Siarry P, Chatterjee A, et al. A fractional‐order variational framework for retinex: fractional‐order partial differential equation
based formulation for multi‐scale nonlocal contrast enhancement with texture preserving. IEEE T Image Process. 2017;27(3):1214‐1229.
33. Jiang H, Zhao Y, Gao J, Gao Z. Adaptive pseudo‐color enhancement method of weld radiographic images based on hsi color space and
self‐transformation of pixels. Rev Sci Instrum. 2017;88(6):065106.
34. Zhou S, Zhuang X, Zhu H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theoretical and
Applied Fracture Mechanics. 2018;96:174‐192.
35. Abidi BR, Zheng Y, Gribok AV, Abidi MA. Improving weapon detection in single energy x‐ray images through pseudocoloring. IEEE t
SYST Man Cy c. 2006;36(6):784‐796.
36. Hawkins DM, Krooden JAT. A review of several methods of segmentation. Geomathematical and Petrophysical Studies in Sedimentology.
1979;117‐126.
37. Gentil F, Parente M, Martins P, Garbe C. Paã§O J, Ferreira AJM, Tavaresc JMRS, Jorge RN. The influence of muscles activation on the
dynamical behaviour of the tympano‐ossicular system of the middle ear. Comput Method Biomec. 2013;16(4):392‐402.
38. Shi L, Wang D, Chu WC, et al. Automatic mri segmentation and morphoanatomy analysis of the vestibular system in adolescent idio-
pathic scoliosis. Neuroimage. 2011;54(1):S180‐S188.
39. Bradshaw AP, Curthoys IS, Todd MJ, et al. A mathematical model of human semicircular canal geometry: a new basis for interpreting
vestibular physiology. Jaro‐J Assoc Res Oto. 2011;11(2):145‐159.
40. Noble JH, Rutherford RB, Labadie RF, Majdani O, Dawant BM. Modeling and segmentation of intra‐cochlear anatomy in conventional
CT. SPIE Medical Imaging International Society for Optics and Photonics. 2010;7623:762302.
41. Chan TF, Vese LA. Active contours without edges. IEEE Press. 2001;10:266‐277.
42. Westover L. Footprint evaluation for volume rendering. ACM SIGGRAPH Computer Graphics. 1990;24(4):367‐376.

How to cite this article: Zhao Z, Zhou X‐P. Digital measurement of 2D and 3D cracks in sandstones through
improved pseudo color image enhancement and 3D reconstruction method. Int J Numer Anal Methods Geomech.
2019;1–20. https://doi.org/10.1002/nag.2993