=Paper=
{{Paper
|id=Vol-2210/paper37
|storemode=property
|title=Unsupervised segmentation of ceramic proppant particles in 3D microCT images
|pdfUrl=https://ceur-ws.org/Vol-2210/paper37.pdf
|volume=Vol-2210
|authors=Ekaterina Serkova,Ilia Safonov,Ivan Yakimchuk,Victoria Evstefeeva
}}
==Unsupervised segmentation of ceramic proppant particles in 3D microCT images==
https://ceur-ws.org/Vol-2210/paper37.pdf
Unsupervised segmentation of ceramic proppant particles in
3D microCT images
E P Serkova1,2, I V Safonov1, I V Yakimchuk1 and V Yu Evstefeeva1,2
1
Schlumberger Moscow Research, Pudovkina 13, Moscow, Russia, 119285
2
Moscow State University, Leninskie Gory 1, Moscow, Russia, 119991
Abstract. Oil and gas industry uses different types of ceramic proppants in millions of
kilograms per year. X-ray microtomography (microCT) imaging can be applied for
investigation of quality of the material. For analysis, it is necessary to segment spherical
contacting particles of proppant. We apply a marker-controlled watershed for segmentation.
The method of markers detection has several parameters which have crucial influence on
segmentation outcome. To optimize segmentation quality, we propose unsupervised (non-
reference) measure based on a compactness of 3D connected regions, where compactness is
calculated via central geometric moments of second order. In addition, we demonstrate
advantages of our technique for compactness estimation over the method based on ratio of
surface area and volume of a region.
1. Introduction
Oil and gas industry uses various types of proppants (from phrase “propping agent”) in hydraulic
fracturing technology. Annually in the world hundreds of thousands of tons of ceramic proppant are
produced, the particles of which are granules of spherical shape with a size of about 1 mm.
Mechanical strength and conductivity are the most important attributes of proppant pack for optimal
fracturing job design. Crush test is one of conventional procedures for measurement of proppant
characteristics. Proppant grains crush test should be performed at axial stress up to 100 MPa. Some
particles are crushed during the test. The fraction of the crushed particles depends on stress and quality
of the proppant. For some proppants the crush-rate is only a few percent. X-ray micro-tomography
(microCT) makes capable measurement of morphometric characteristics of each particle as well as
their fragments in initial state and after stress.
Figure 1 shows example of real reconstructed 3D microCT image of proppant pack before stress.
The grayscale image has size 4000х4000х2000 voxels and 8 bit-depth. The image was obtained with
SkyScan 1172 microCT system (Bruker MicroCT, Belgium). One can see proppant particles in the
image. Also fragments of crushed particles can be found in images scanned after loading. The usual
image of proppant pack contains several hundreds of particles. Intensities of particle regions are much
lighter than dark background. However, there is many contacting regions of granules and their
fragments. Separation of the contacting regions is a challenging task.
Marker-controlled watershed for distance map is a traditional method for segmentation of touching
regions. An algorithm for detection of markers has several parameters, which influence to
segmentation quality considerably. Since the fraction of broken particles can be equal to several
percent, even single segmentation errors leads to bias of an assessment of quality of the material.
Selecting the segmentation parameters manually by the operator is a long and unobvious process.
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018)
Image Processing and Earth Remote Sensing
E P Serkova, I V Safonov, I V Yakimchuk and V Yu Evstefeeva
Figure 1. Example of reconstructed 3D microCT image of proppant particles.
In addition, it is quite difficult to detect segmentation errors in a 3D image visually. Our
preliminary experiments have shown the optimal segmentation parameters vary from image to image;
it is impossible to set parameters once based on a previously processed sample. Therefore, it is
important to develop an unsupervised (non-reference) quality metrics that allow selection of the
parameters automatically with aim to minimize the number of errors of segmentation of proppant
particles.
The main contribution of the paper is proposed unsupervised metrics of segmentation of 3D
spherical particles by means of maximization of average compactness of segmented regions, where
compactness is calculated via second order central geometrical moments.
2. Segmentation of particles
Figure 2 demonstrates flow-chart of segmentation algorithm. Reconstructed grayscale 3D image is
downsampled to image G having size 1000x1000x500 voxels and 8-bits depth. It is necessary to
reduce requirements to memory space and decrease processing time. Lighter voxels of ceramic
particles differ from dark voxels of background and holder. Histogram of intensities of the image is bi-
modal. For images with such histograms thresholding by Otsu algorithm [1] is a good solution to
distinguish voxels of particles from background. After thresholding we obtain binary image T, where
voxels of solid are designated by 1 and voxels of voids are designated by 0. It is required to split
regions formed by touching particles in T image.
Conventional way for separation of overlapping or contacting convex regions without holes is an
application of watershed algorithm to the inverted outcome of the geodesic distance transform [2]. The
general idea of watershed algorithm is the following: image is considered as a geological relief; a
water source is placed in each regional minimum in the relief, to flood the entire relief from sources,
and build barriers in the place where different water sources meet; the resulting set of barriers
constitutes a watershed by flooding [3]. For volumetric images watershed algorithm operates
identically to 2D one.
Before application of distance transform we need to fill holes, which are pores in particles. Several
pores are open, they connect with background voids. That is why filling of holes in 3D keeps these
pores unchanged. It is required to perform filling of holes for 2D slices. Theoretically, the filling of
holes should be done for slices in all three mutually perpendicular directions. Such approach ensures
that the open pores, which are penetrate a particle and are parallel to image axes, are filled. However,
in practice open pores are tortuous, so, it is enough to fill the holes for 2D slices in only one direction.
Binary image Tf is outcome of slice-by-slice filling of holes for T image. It is worth to note, filling of
holes in 2D can lead to filling of space between touching particles. Fortunately, it happens seldom and
it has no negative impact on the next processing steps.
On the next stage, geodesic distance transform builds the distance map D by calculation for each
voxel of Tf image Euclidean distance to the nearest voxel equal to zero. Inverted D plays the role of
relief for watershed algorithm. Local minima on (-D) are basins origins. Frequently there are several
local minima inside a connected region because even survived particles have non-ideal convex shape,
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 283
Image Processing and Earth Remote Sensing
E P Serkova, I V Safonov, I V Yakimchuk and V Yu Evstefeeva
a lot of fragments of crushed particles are concave. It leads to over-segmentation. To avoid over-
segmentation due to the huge number of local minima the marker-controlled watershed is used, where
markers play role of water sources [4].
G
Thresholding by Otsu
T
Filling of holes for 2D slices
Tf
Geodesic distance transform
D
Grayscale morphological reconstruction by dilation
R
Maximal filter
Rf
Markers detection: R = = Rf
M
Marker-controlled watershed for inverted D
Lw
Masking of labels by T
L
Figure 2. Flow-chart of our segmentation technique.
There is a plenty of approaches for markers generation. We tried several of them and found out the
following steps providing the best segmentation outcome. The first step is grayscale morphological
reconstruction by dilation [5]. Morphological reconstruction by dilation uses two images: “seed”
image (D – delta_h), which specifies the values that spread, and “mask” image D, which gives the
maximum allowed value at each voxel. The mask image limits the spread of high-intensity values. The
resulting reconstructed image R looks exactly like seed image, but with the peaks cut off, where
delta_h determines height of peaks. Difference (D – R) is referenced as ‘h-dome’ operator, ((D – R) >
0) is referenced as ‘h-maxima’ operator.
Figure 3 demonstrates example of markers generation for simple image (see figure 3a) containing
two contacted regions. Figure 3b shows distance map. Figure 3c shows plot of profile for the distance
map along axis K. One can see four local maxima in the plot, where three local maxima relate to the
right region. We need to suppress redundant local maxima. Adjustment of delta_h mitigates issue of
several local maxima for some particles, for example two local maxima of central peak are combined
to single region. However, increasing of delta_h cannot solve the issue completely, e.g. the rightmost
peak has local maximum that should be ignored. Application of maximal filter to image R and
comparison of filtration outcome Rf with R employ to suppress the most of redundant local maxima.
Binary image M containing markers is: M = (R == R f) (see figure 3c). The maximal filter uses cubical
structural element with max_filter_size size. These three parts of figure 3 illustrate meaning of
parameters delta_h and max_filter_size.
The next stage is marker-controlled watershed for inverted D and markers from M, where each
connected region has unique label. Image Lw is outcome of segmentation by watershed. Finally, to
obtain labelled particles in T image it is necessary to multiply Lw by T.
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 284
Image Processing and Earth Remote Sensing
E P Serkova, I V Safonov, I V Yakimchuk and V Yu Evstefeeva
K
(a) (b)
max_filter_size
delta_h
delta_h
Distance
max_filter_size
delta_h
K
(c)
Figure 3. Illustration of parameters for generation of markers.
3. Unsupervised segmentation quality metrics
Segmentation outcome depends on parameters delta_h and max_filter_size. Decrease of both leads to
over-segmentation, because regions of granules split in fragments. Increase of both leads to under-
segmentation, because neighbouring particles and fragments merge. How to set the parameters
properly? Subjective choice based on visual analysis of segmentation result is a troublesome due to 3D
nature of data. It is hard to find optimal parameters visually. Paper [6] analyses supervised quality
measures for segmentation such as Global Consistency Error and Rand Index. However, we have no
ground truth to be able to apply supervised metrics. How to automatically choose the optimal
parameter values, based only on the analysed image? Survey [7] describes a few dozen of
unsupervised metrics for estimation of segmentation quality of 2D images. All of them are not
universal, but depend on the task being solved. When segmenting regions have more or less the same
shape, the shape factor can be used as metrics.
We know a-priori that studied particles have a rounded shape. We use the fact for formulation of
criterion of unsupervised segmentation. A sphere is the most compact body in three-dimensional
space. In the next section, we describe a compactness for characterizing the closeness of a region
shape to a sphere. Maximal average compactness of segmented regions corresponds to the best
segmentation:
𝑁
1
𝑄 = ∑ 𝐶𝑖 ,
𝑁
𝑖
where N is the number of segmented regions having volume greater then Vt; Ci is compactness of i-th
region. The threshold Vt is introduced to exclude from consideration too small fragments. In addition,
it is hard issue: what is compactness for regions consisting of only a few voxels?
Let’s consider how Q changes in the case of improper segmentation. Over-segmentation is
splitting of spherical regions. It leads to decreasing of average compactness Q. Under-segmentation is
merging of spherical regions. As rule, it leads to decreasing of average compactness too. It is worth to
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 285
Image Processing and Earth Remote Sensing
E P Serkova, I V Safonov, I V Yakimchuk and V Yu Evstefeeva
note, proposed segmentation criterion is correct in assumption of existing of at least several tens of
survived spherical particles in image after loading. If almost all particles are crushed then the metrics
is ineffective.
4. Compactness
Many 3D shape factors are a natural extension of corresponding measures for 2D images.
Compactness in 2D is sometimes called the circularity or the roundness. Shape compactness in 2D is
generally understood as the degree to which a given shape differs from a region bounded by a circle.
There are plenty publications, which describe compactness as ratio of area of region to squared
perimeter. Analogous definition of compactness for 3D regions is ratio of squared volume to area of
surface of a region in cubic power [9]:
36𝜋𝑉 2
𝐶𝑎 = ,
𝑆3
where V is volume of 3D region, S is area of surface.
Paper [8] performs detailed analysis of drawbacks of approaches, which use perimeter and surface
area for calculation of compactness for 2D and 3D images correspondingly. First, a perimeter and a
surface area are difficult to calculate invariantly to rotation and quite precisely due to the digital nature
of images. Second, the surface of the regions is not ideal, there are various types of noises that arise
during the registration and segmentation of regions. Third, the evaluation of compactness is seriously
affected by holes in regions.
Bribiesca [11] describes compactness for 3D regions that differs from 𝐶𝑎 , but it is based on volume
and area of surface also. The paper [11] demonstrates Bribiesca’s is less sensitive to small distortions
of surface in comparison with 𝐶𝑎 .
There are several approaches for calculation of surface area. For example, the paper [9] considers
marching cubes algorithm [10] for estimation of outer surface area. A surface area can be estimated as
the number of voxels of the outer shell, where the shell is the difference in the region and the result of
its erosion with a 3x3x3 structural element in the form of a cube (26-connectivity), or a ball (18-
connectivity), or a cross (6-connectivity). Surely, all pores should be filled in advance. Which method
is preferable?
Zunic et al. [12] proves properties of 3D shapes compactness based on second order central
geometric moments and demonstrates its advantages over other algorithms for compactness
calculation:
35/3 𝜇000 5/3
𝐶𝑚 = ,
5(4𝜋)2/3 𝜇200 + 𝜇020 + 𝜇002
where central geometric moments are:
𝜇𝑝𝑞𝑟 = ∑ ∑ ∑ 𝐼(𝑥, 𝑦, 𝑧)(𝑥 − 𝑚100 )𝑝 (𝑦 − 𝑚010 )𝑞 (𝑧 − 𝑚001 )𝑟 ,
𝑥 𝑦 𝑧
where geometric moments are:
𝑚𝑝𝑞𝑟 = ∑ ∑ ∑ 𝐼(𝑥, 𝑦, 𝑧)𝑥 𝑝 𝑦 𝑞 𝑧 𝑟 ,
𝑥 𝑦 𝑧
where 𝐼(𝑥, 𝑦, 𝑧) is indicator function of particle region, I equals one for voxels belonging the region,
and I equals zero otherwise. Zero order geometric moments equal to volume of particle: 𝑚000 =
𝜇000 = 𝑉. First order geometric moments are centroids. Compactness lies in the range from 0 to 1.
Compactness of ideal sphere equals 1.
Paper [12] is based on outcomes from Mamistvalov‘s theory. Mamistvalov published in English
mathematical theory for recognition of n-dimensional solids via geometric moment invariants in 1998
[13]. In Russian, it was done in 1974. Unfortunately, those papers are not widely known. Ratio
𝜇200 +𝜇020 +𝜇002
𝜇 5/3 is so-called first 3D moment invariant based on second order central geometric
000
moments. The invariant is constant or slightly changes for origin translation, uniform scaling and
rotation. 𝐶𝑚 is the inverted first 3D moment invariant multiplied by the constant, so, it is an invariant
for translation, rotation and scaling as well. In general, geometric moments have been widely used in
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 286
Image Processing and Earth Remote Sensing
E P Serkova, I V Safonov, I V Yakimchuk and V Yu Evstefeeva
statistics for description of the form of a probability density function and in classic rigid-body
mechanics to measure the mass distribution of a body. Usage of moment invariants is a prospective
direction in an image processing for creation of shape factors as well as for development of algorithms
for shapes analysis and recognition. According to monography [14], last decade an interest to shape
analysis and classification via moment invariants grows. Not only geometric invariants are applied,
but orthogonal moments such as Gaussian-Hermite, Zernike, Chebyshev, Legendre and Fourier-
Mellin.
Let’s consider which method of computing the compactness is better suited for our problem. We
calculate the compactness coefficients for spheres of different radii with different random noises such
as protuberances and cavities on a surface. Figure 4 shows a slice of the sphere with typical noise on
the surface. In the next experiment, we determine which method for compactness calculation allows us
to distinguish sphere from regions formed by cutting off segments of different sizes from the sphere
better.
Figure 4. Slice of sphere with cavities and outliers on a surface.
Plots in figure 5 show various compactnesses depending on radius of sphere: 𝐶𝑚 ; 𝐶𝑎 , where
surface area is calculated by marching cubes algorithm; 𝐶𝑎 , where surface area is calculated by means
of erosion, Bribiesca’s compactness; theoretical ideal case, where compactness equals 1. 𝐶𝑎 for both
considered approaches are very different from 1 and have significant fluctuations. 𝐶𝑚 is close to 1
except for spheres having radius less than 10. Bribiesca’s compactness is almost 1 for all range of
considered radiuses.
Figure 5. Compactness of spheres of different radiuses with noisy surface.
Plots in figure 6 show compactnesses for bodies from hemisphere (L=60) to sphere (L=0). Figure 7
illustrates meaning of L parameter. 𝐶𝑚 and 𝐶𝑎 allow to distinguish sphere and sphere with clipped
segment. Bribiesca’s compactness is too similar even for sphere and hemisphere.
We conclude, application of compactness 𝐶𝑚 calculated via second order central geometrical
moments has obvious advantages in comparison with well-known approaches based on surface area
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 287
Image Processing and Earth Remote Sensing
E P Serkova, I V Safonov, I V Yakimchuk and V Yu Evstefeeva
and volume of 3D region. In addition, in contrast to 𝐶а , 𝐶𝑚 has robustness to presence of the modest
number of open and close pores.
Figure 6. Compactness when the shape of the body changes from the hemisphere to the sphere.
Figure 7. Illustration of meaning of L parameter for figure 6.
5. Results and discussion
For estimation of benefits of proposed unsupervised segmentation technique, we processed five 3D
images with fixed parameters delta_h and max_filter_size as well as with parameters obtained by
maximization of 𝑄 metrics. It is worth to note, those fixed parameters were optimal according to 𝑄
criterion for two images segmented previously. We apply gradient descent algorithm for looking for
maximal 𝑄. Figure 8 shows 𝑄 metrics depending on parameters delta_h and max_filter_size.
Figure 8. Segmentation quality metrics Q depending on parameters delta_h and max_filter_size.
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 288
Image Processing and Earth Remote Sensing
E P Serkova, I V Safonov, I V Yakimchuk and V Yu Evstefeeva
Table 1. The number of erroneously segmented particles.
Image Fixed parameters Parameters according to
maximal Q
1 37 0
2 37 0
3 54 0
4 44 2
5 38 1
Table 1 contains for each tested image the number of erroneously segmented particles by
segmentation with fixed parameters and by proposed unsupervised segmentation. Total number of
particles in each image is 615. Segmentation with fixed parameters leads to 30-50 errors. Figure 9
shows slice of segmented and labeled image by segmentation with fixed parameters. One can see, one
ellipsoidal particle was segmented as two regions, ten pairs of neighboring particles were pairwise
combined into one region. Unsupervised segmentation had no errors, or just one or two errors in the
worst case. Figure 10 demonstrates the slice of 3D image, where unsupervised segmentation was
applied. Segmentation outcome is fully correct.
So, proposed unsupervised criterion for the segmentation of 3D images of ceramic proppant allows
providing high quality segmentation in automatic mode.
Figure 9. Slice of segmented 3D image with Figure 10. Slice of segmented 3D image with
several separated and merged particles; fixed proper segmentation; parameters according to
parameters were used for segmentation. maximal Q were selected.
6. References
[1] Otsu N A 1979 Threshold selection method from gray-level histograms IEEE Trans. Syst. Man
Cybern. 9(1) 62-66
[2] Vincent L and Soille P 1991 Watersheds in digital spaces: an efficient algorithm based on
immersion simulations IEEE Trans. Pattern Anal. Mach. Intell. 6 583-598
[3] Gonzalez R C and Woods R E 2012 Digital image processing
[4] Safonov I V, Mavrin G N and Kryzhanovsky K A 2006 Segmentation of convex cells with
partially undefined boundaries Pattern Recogn. Im. Analysis 16(1) 46-49
[5] Vincent L 1993 Morphological grayscale reconstruction in image analysis: Applications and
efficient algorithms IEEE Trans. Image Process 2(2) 176-201
[6] Myasnikov E V 2017 Hyperspectral image segmentation using dimensionality reduction and
classical segmentation approaches Computer Optics 41(4) 564-72 DOI: 10.18287/2412-6179-
2017-41-4-564-572
[7] Zhang H, Fritts J E and Goldman S A 2008 Image segmentation evaluation: A survey of
unsupervised methods Comput. Vis. Image Underst. 110(2) 260-280
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 289
Image Processing and Earth Remote Sensing
E P Serkova, I V Safonov, I V Yakimchuk and V Yu Evstefeeva
[8] Montero R S and Bribiesca E 2009 State of the art of compactness and circularity measure
Internat. Math. Forum 4(27) 1305-1335
[9] Zhao B and Wang J 2016 3D quantitative shape analysis on form, roundness, and compactness
with μCT Powder Technol. 291 262-275
[10] Lorensen W E and Cline H E 1987 Marching cubes: A high resolution 3D surface construction
algorithm ACM Siggraph Comp. Graphics 21(4) 163-169
[11] Bribiesca E 2008 An easy measure of compactness for 2D and 3D shapes Pattern Recognit.
41(2) 543-554
[12] Žunić J, Hirota K and Martinez-Ortiz C 2012 Compactness measure for 3d shapes ICIEV 1180-
1184
[13] Mamistvalov A G 1998 N-dimensional moment invariants and conceptual mathematical theory
of recognition n-dimensional solids IEEE Trans. Pattern. Anal. Mach. Intell. 20(8) 819-831
[14] Flusser J, Suk T and Zitová B 2016 2D and 3D Image Analysis by Moments (John Wiley &
Sons)
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 290