=Paper=
{{Paper
|id=Vol-2391/paper22
|storemode=property
|title=Selection in a 3D microtomographic image the region with the highest quality
|pdfUrl=https://ceur-ws.org/Vol-2391/paper22.pdf
|volume=Vol-2391
|authors=Anton Kornilov,Ilia Safonov,Alena Goncharova,Ivan Yakimchuk
}}
==Selection in a 3D microtomographic image the region with the highest quality ==
https://ceur-ws.org/Vol-2391/paper22.pdf
Selection in a 3D microtomographic image the region with the
highest quality
A S Kornilov1,2, I V Safonov1, A V Goncharova2 and I V Yakimchuk1
1
Schlumberger Moscow Research, Pudovkina 13, Moscow, Russia, 119285
2
National Research Nuclear University MEPhI, Kashirskoe shosse 31, Moscow, Russia,
115409
e-mail: kranton94@mail.ru
Abstract. We present an algorithm for processing of X-ray microtomographic (micro-CT)
images that allows automatic selection of a sub-volume having the best visual quality for
further mathematical simulation, for example, flow simulation. Frequently, an investigated
sample occupies only a part of a volumetric image or the sample can be into a holder; a part of
the image can be cropped. For each 2D slice across the Z-axis of an image, the proposed
method locates a region corresponding to the sample. We explored applications of several
existing blind quality measures for an estimation of the visual quality of a micro-CT image
slice. Some of these metrics can be applied to ranking the image regions according to their
quality. Our method searches for a cubic area located inside regions belonging to the sample
and providing the maximal sum of the quality measures of slices crossing the cube across the
Z-axis. The proposed technique was tested on synthetic and real micro-CT images of rocks.
1. Introduction
In the last decade, X-ray microtomography (micro-CT) has been actively used for analyzing of rock,
soil, granular materials [1] and building their 3D models for a Digital Rock [2]. One of the
applications of these digital models is a flow modelling of a multiphase fluid [3]. A typical micro-CT
image can have size up to ~40003 voxels. Some systems allow getting even larger images. However,
the sample under study occupies only part of the image or can be in the holder. For mathematical
modelling, researchers often use only a sample fragment of a cubic form with sides from 400 to 1000
voxels. Thus, there is a task of cropping a fragment from a whole micro-CT image, which is suitable
and, in some sense, optimal for further analysis. Currently, an operator performs this procedure
manually. The quality of the micro-CT image varies in the image volume, individual parts of the
sample, for example, dense inclusions, can significantly affect the quality [4]. For a reasonable
fragment selection, the operator needs to view several hundred slices of a 3D image. The selection
process takes considerable time and depends on the subjective opinion of the operator. So, the
development of an automatic quick and formalized algorithm of the fragment selection for
mathematical modelling is an important problem.
The complexity of this problem is explained by the huge variety of images that can be on the input
of the algorithm. This diversity is caused not only by the variety of scanned samples but also by
V International Conference on "Information Technology and Nanotechnology" (ITNT-2019)
Image Processing and Earth Remote Sensing
A S Kornilov, I V Safonov, A V Goncharova and I V Yakimchuk
features of various microtomographic systems and parameters of image reconstruction algorithms. In
addition, certain tomography artifacts can be differently suppressed by various reconstruction
software. Frequently, only a part of an original reconstructed tomographic image relates to a sample.
And also the software of a micro-CT system reconstructs from shadow projections an image having a
bit depth of 8 or 16 bits per voxel.
Figure 1 demonstrates several examples of 2D slices perpendicular to the Z axis. These images
were obtained using different micro-CT systems. A granular material is placed in a cylindrical holder,
which is completely inside the field of view, non-reconstructed areas in the corners of the image are
filled by gray (figure 1a). A cylindrical sample of sandstone is completely in the field of view, the
areas in the corners are filled by gray (figure 1b). A cylindrical sandstone sample is completely in the
field of view but is off-center; unreconstructed areas are filled by black (figure 1c). Note, the holder is
not of interest for further analysis. In the following example, scanning is performed inside a sample,
unreconstructed corner areas are darker and blurred (figure 1d) than a center zone. A piece of an
arbitrary shape sample is completely in the field of view (figure 1e). A rectangular fragment is
cropped from an original reconstructed image (figure 1f). In addition, several dozen slices at the top
and bottom of a 3D image may not contain the sample. They can contain a holder or air.
a) b) c)
d) e) f)
Figure 1. Examples of slices of various micro-CT images.
One more challenge is a formalization of a criterion for a selection of a fragment for further
modelling. Which one of the many possible fragments to choose? It is reasonable to crop the fragment
optimal from the point of view of some criterion. What can the criterion be in this case? Ideally, the
criterion should depend on the objectives of the subsequent modeling. For example, to simulate the
fluid flow, the criterion may be based on porosity or coherence of the pore structure. However,
regardless of the methods and objectives of the subsequent analysis, we can say that we are interested
to determine the most qualitative fragment in terms of visual assessment.
An important requirement for the developed algorithm and software is a high speed of the micro-
CT image analysis. The processing time should be substantially less than the time required for the
selection of the suitable fragment manually. We propose a fast algorithm that automatically selects
from a micro-CT image a fragment of a given size that has the highest visual quality assessment.
2. Algorithm for selection of micro-CT image region
2.1. Main ideas
We propose to select a convex region related to the sample on each slice of the 3D image and evaluate
the image quality in this area. Then we determine whether it is possible to inscribe a cube of a given
V International Conference on "Information Technology and Nanotechnology" (ITNT-2019) 161
Image Processing and Earth Remote Sensing
A S Kornilov, I V Safonov, A V Goncharova and I V Yakimchuk
size into selected regions of consecutive slices and calculate the image quality estimation for the cube,
summing up the estimations of slices intersecting the cube. The result of the algorithm is the cube that
corresponds to a fragment of the 3D image with the highest image quality assessment. The proposed
approach allows organizing parallel processing of the slices to select the regions located inside the
sample and to assess the image quality in this area. Image handling in the slice-wise mode is one from
approaches to reduce the processing time of a tomographic image [5].
The method of selecting the region related to the sample assumes that there are areas of varying
intensity or filled with different textures in the image of the sample, between which there are visible
edges. We propose to apply the thresholding of the outcome of the variance filter [6] to detect the
edges between regions. At the same time, if on the image of the slice there are only the edges of the
holder, then they should be ignored. From the obtained convex regions for the stack of slices, we can
easily determine whether the cube of a given size can be inscribed in them. Moreover, we can find the
maximum possible size of the inscribed cube using an algorithm that determines the maximum
rectangle that fits inside the connected region [7].
Over the past decade, several no-reference (blind) image quality metrics have been proposed [8-
11]. Algorithms for calculating these quality metrics were trained and tuned on a set of natural 2D
photos. They can assess the quality of 2D grayscale images with a bit depth of 8 bits per pixel. It was
assumed that among the existing non-reference quality metrics we can find applicable for
characterizing the quality of the slices of micro-CT images.
To use mentioned no-reference image quality metrics, we should convert images from 16 bits per
voxel to 8 bits per voxel. Generally, intensities of images having 16 bits per voxel do not occupy the
entire dynamic range. Thus, we need to increase the contrast in converting the images to 8 bits per
voxel. Sometimes, it is advisable to increase the contrast in case of images with 8 bits per voxel too.
The intensities of all 3D image slices should be normalized in the same range. The parameters for
normalization are determined by means of analyzing the intensity histogram H of the central part of
the 3D image. To calculate the histogram H, we perform a preliminary pass through the central part of
the image, for which we indent 25% from each side of the image.
The lower bound l for intensities normalization is calculated using the statement [12]:
l= ( { } {
min min i H [i ] ≥ H 0 , min i ∑ k =0 H [ k ] ≥ C0 ,
i
}) (1)
where H0 is the threshold for the height of the histogram column for dark tones; С0 is the threshold for
the area of the histogram in dark tones, functionally С0 is equivalent to the percentile of H.
The upper limit u for intensities normalization is calculated with the formula:
u=
{ } {
max max i H [i ] ≥ H1 , max i ∑ k =i H [ k ] ≥ C1 ,
2n −1
}
(2)
where H1 is the threshold for the height of the histogram column for bright colors; С1 is the threshold
for the area of the histogram in bright colors; n is the bit depth, 8 or 16.
Normalization of the intensities of each slice I in the range [0, 255] is performed according to the
formula:
0 : I ( r, c ) ≤ l,
255 ( I ( r , c ) − l )
= I ( r, c ) : l < I ( r , c ) < u, (3)
( u − l )
255 : u ≤ I ( r , c )
where (r, c) is the slice pixel coordinates, r = 0 ... N – 1, c = 0 ... M – 1; N is the number of rows in the
image; M is the number of columns.
2.2. Region of sample selection in 2D slice
The algorithm for selecting the region belonging to the sample searches for a convex region, within
which there are distinguishable edges between regions of different intensity. It is important to find the
area inside the sample, but it is not necessary to accurately detect the boundary of the sample. Often,
V International Conference on "Information Technology and Nanotechnology" (ITNT-2019) 162
Image Processing and Earth Remote Sensing
A S Kornilov, I V Safonov, A V Goncharova and I V Yakimchuk
the noise level in micro-CT images is quite high. Therefore, many well-known edge-detection filters
are not applicable, as they react to noise. In addition, we are not interested in the edges themselves, but
areas where there are changes in intensity. The variance filter allows getting an image in which
relatively large values are in the pixels, around which there is a change in the intensity and relatively
small values in uniform areas.
The variance filter is defined as [6]:
2
1 W −1W −1 W W
= I var ( r , c ) 2 ∑ ∑
I r + i − , c + j − − IW ( r , c ) , (4)
W =i 0=j 0 2 2
where [ ] is the operator of taking the integer part; I(r, c) is the intensity of a pixel with coordinates
(r, c), r = [W / 2] ... N – [W / 2] – 1, c = [W / 2] ... M – [W / 2] – 1; the mean value IW within the
window is calculated as:
1 W −1W −1 W W
= IW ( r , c ) 2 ∑ ∑
I r + i − , c + j − . (5)
W =i 0=j 0 2 2
The window size depends on the slice size, and it is selected about 100 by 100 pixels. Figure 2b
shows the result of the variance filter for the slice in figure 1a. Thresholding of the outcome of the
variance filter allows obtaining a binary image, in which the pixels are set to value 1 when the
intensity varies in the local vicinity. Figure 2b shows the result of the thresholding. It should be noted
that after the thresholding, several connected regions may form. To merge them into a single
connected region, we use the morphological closing operation [13]. Several small areas can be located
near the boundaries of the image. For instance, they can occur on the outer border of the holder. To
exclude these undesirable regions from consideration we process only the connected region with the
largest area. For this region, we calculate the convex hull and fill it with values 1. The resulting region
protrudes approximately half to two-thirds of the size of the variance filter window from the sample
region. To ensure that the region is inside the sample region, a morphological erosion filter [13] is
applied with the size of the structural element equal to the size of the variance filter window.
Note that the image of the slice may be missing a fragment relating to the sample. The proposed
algorithm allows detecting such a situation. We make a conclusion about this situation based on the
area and shape of the largest detected region after the application of the closing operation. Also,
according to the results of the sample area selection on the slice, the algorithm produces
recommendations for the operator of the micro-CT system. For example, if the sample is significantly
displaced from the center of the slice, then it is probably possible to reconstruct an image of better
quality if you re-shoot the shadow projections by placing the sample in the center.
a) b) c)
Figure 2. Illustration of the sample region selection stages for the slice from figure 1a: (a) the result of
the variance filter; (b) the result of the thresholding of image from figure 2a; (c) a convex region
constructed along the outer boundary of region from figure 2b.
2.3. Analysis of the use of no-reference quality metrics for micro-CT images
At present, a sufficiently big number of no-reference image quality metrics has been proposed. Such
metrics make possible to evaluate the quality based on the image only without the reference. Part of
the metrics is designed to assess one factor affecting the quality, for example, the blurriness level [14].
V International Conference on "Information Technology and Nanotechnology" (ITNT-2019) 163
Image Processing and Earth Remote Sensing
A S Kornilov, I V Safonov, A V Goncharova and I V Yakimchuk
Others claim to universality. We evaluated the following algorithms for assessing of quality of slices
of micro-CT image: BIQI (Blind Image Quality Index) [8], BRISQUE (Blind Referenceless Image
Spatial Quality Evaluator) [9], NIQE (Natural Image Quality Evaluator) [10] and OG-IQA (Oriented
Gradients Image Quality Assessment) [11].
The considered algorithms were trained on natural photos, for instance, from the LIVE data set
[15]. This data set, in addition to undistorted photos, contains images that have been distorted. Several
types of distortion were used: additive white Gaussian noise, Gaussian blur, JPEG and JPEG2000
compression artifacts. However, this dataset is rather idealistic, because only one type of distortion
affects an image at the same time. After interviewing many observers, an estimate of quality for each
image from the set was calculated based on the pairwise comparison.
The BIQI algorithm [8] implements a two-step approach to assess the quality of the images. This
approach is based on a natural scene statistic (NSS) and the assumptions that natural images have
certain statistical properties, and distortions change these statistical characteristics in such a way that
the type and degree of distortion can be estimated with a high confidence level. The first step in BIQI
is an identification of distortion type. The second stage is an estimation of the strength of the
distortion. Features for learning are taken from the model of natural scenes in the wavelet domain [16].
The models for classification and regression are trained by support vector machine (SVM) on data
from the LIVE dataset. The BRISQUE method [9] uses the features of the NSS model calculated in
the spatial domain [17]. A single regression model for all distortion types is trained on images from
LIVE dataset using SVM method. The NIQE [10] technique constructs of a multivariate Gaussian
distribution (MVG) of statistical features based on the statistical model of natural scenes in the spatial
domain [17]. The quality of the estimated image is defined as the distance between the
multidimensional normal distribution of this image characteristics and the reference MVG obtained as
a result of training on the features of pristine undistorted images. The OG-IQA [11] algorithm
analyzes the correlation structure of the gradient orientation in the image. The idea is that the
orientation of the local gradients contains significant information about the naturalness of the image,
and when distortions appear, this information changes in a certain way. The AdaBoost algorithm is
used to train the regression model.
Initially, we assumed that these metrics can be used for an absolute quality assessment of micro-CT
images, i.e. regardless of the equipment used and the sample analyzed. However, from the
experiments, we saw that the metrics are not applicable for assessing the quality of images for samples
of different rocks, because frequently the results did not correspond to the expert’s opinion and the
physical parameters of the experiments. However, for a comparative quality assessment of slices from
one sample, the results were quite reasonable.
To study the applicability of these criteria for the comparative quality assessment of the slices from
one sample, we formed a set of 10 images obtained as follows. We obtained six images by scanning
the same sample with different exposure times. The noise level is directly dependent on the exposure
time. The highest noise level corresponds to the shortest exposure. For two of these six images, we
additionally reconstructed two images with a binning and a reduced number of projections. Binning is
the process of combining the charges of neighboring pixels on an image sensor into one “superpixel”.
In tomographic reconstruction, binning is lowering the spatial resolution due to summation the values
of neighboring pixels, which leads to an increase in the signal-to-noise ratio. Using fewer projections
for reconstruction causes noise-like artifacts.
For these 10 images, it is known how they are ranked according to their quality. For the central part
of these images, we calculated the quality assessment by several quality metrics. Table 1 contains the
number of errors in image rankings using considering no-reference quality metrics. The number of
errors is counted as the number of images that must be removed from the series to obtain a properly
sorted sequence by quality. The best result was shown by NIQE metrics, which made two mistakes
with images reconstructed with binning.
NIQE advantage was expected. This metrics is less associated with a fixed set of distortions of the
LIVE dataset and uses a simpler model which is less tended to overfit. Moreover, for NIQE we built
own MVG model of high-quality micro-CT images instead of photos from LIVE dataset. Further
adjustment of this model may lead to improved results.
V International Conference on "Information Technology and Nanotechnology" (ITNT-2019) 164
Image Processing and Earth Remote Sensing
A S Kornilov, I V Safonov, A V Goncharova and I V Yakimchuk
Table 1. The number of errors in the quality ranking of micro-CT images using non-reference
quality metrics.
Quality metrics Errors count
BIQI 4
BRISQUE 4
OG-IQA 6
NIQE 2
2.4. Algorithmic optimization
We achieve a high speed of handling of images due to parallel processing and the use of optimized
algorithms. An example of algorithmic optimization is the implementation of a variance filter via
integral images. In the integral image IntI, each pixel with coordinates (r, c) is equal to the sum of the
pixels of the image I located to the left and above (r, c) [18]:
IntI ( r , c ) = ∑ I ( r ′, c′ ). (6)
r ′≤ r
c′≤ c
IntI can be calculated in a single pass through the image I by the following statement:
IntI ( r ,=
c ) I ( r , c ) + I ( r , c − 1) + I ( r − 1, c ) − I ( r − 1, c − 1) . (7)
The sum of the pixels in the square window for the image I via the integral image IntI is calculated
as:
W − 1 W − 1 W W
S ( r , c ) = IntI r + ,c + + IntI r − , c −
2 2 2 2
(8)
W W − 1 W − 1 W
− IntI r − , c + − IntI r + 2 , c − 2 ,
2 2
where W is the window size; the coordinate (r, c) is in the center of the window.
With the help of integral images, we can quickly perform the operations of the summation and
averaging of rectangular image blocks. The area of integral images application is expanded by
calculating such images for squares and higher powers of pixel values [19]. If we additionally
calculate the integral image of IntI2 for I2 (via the formula (7)) and the sum of pixels in the square
window S2 (using formula (8)), then the variance filter can be calculated as:
1 S ( r, c )
2
= I var ( r , c ) S ( r , c ) − . (9)
W 2 W 2
2
Table 2 contains the processing time of a slice of 4000×4000 pixels by three various
implementations of variance filter. We estimated the time for a single-threaded C++ program on a
laptop with CPU Intel Core i7 2.6 GHz. The slowest approach is the implementation of the filter
directly via formulas (4) and (5) in several nested loops. It is much faster to calculate the variance
from the histogram, which during the window movement across the image is not recalculated
completely but is updated. Finally, the fastest way is to calculate the local variances via integral
images. The integral images IntI and IntI2 require a larger amount of memory in comparison with I,
but in our task, this is not a limitation.
Table 2. The variance filter processing time for a slice of 4000×4000 pixels.
Method Processing time, s
In several nested loops 4.3
Via histogram 1.1
Via integral images 0.5
3. Results and conclusions
We tested about one hundred 3D micro-CT images that had from 1000 to 2000 slices. The region
inside the sample was selected correctly for all cases. The results of the automatic selection of regions
for the slices from figure 1 are shown in figure 3.
V International Conference on "Information Technology and Nanotechnology" (ITNT-2019) 165
Image Processing and Earth Remote Sensing
A S Kornilov, I V Safonov, A V Goncharova and I V Yakimchuk
To test the algorithm for cropping of 3D fragment having the highest visual quality, we developed
a software tool for generating a synthetic image of the micro-CT slice (see figure 4). This tool allows
adjusting the sample size, the level of additive or multiplicative Gaussian white noise, the contrast, the
radius of blurring filter, and the intensity of the holder region. From these synthetic slices, we formed
a 3D image in which the noise level increased from the central slice to the top and bottom. For the
slices of this 3D image, we calculated inverted and normalized NIQE quality metrics Q. As expected,
the quality assessment Q is correlated with the noise level and is decreased from the central slice to top
and bottom. Our algorithm identifies the central region as having the highest quality because the sum
of quality assessments for the central slices is maximum (see figure 5). Background of the graph is a
cross-section of the generated 3D image.
a) b) c)
d) e) f)
Figure 3. The results of the algorithm for the selection of areas related to the sample for images from
figure 1.
a) b)
Figure 4. Synthetic image of the slice (a) and the result of the area related to the sample selection (b).
Outcomes for real micro-CT images of rocks and granular materials are in good agreement with
expert opinion, except regions damaged by X-ray tomography-specific local defects, for example, ring
artifacts and intensity distortions due to high-dense inclusions. The existing no-reference image
quality metrics can give a relative estimation of the noise level and sharpness for the slices from a
single image, but do not consider tomography-specific artifacts. Figure 6 demonstrates a plot with a
cross-section of normalized NIQE quality metrics for slices of real micro-CT image.
The processing time of a 4000×4000×2000 voxel 3D image takes less than 6 minutes on a laptop
with CPU Intel Core i7 2.6 GHz, 8 Gb RAM, SSD. The use of the proposed algorithm allows reducing
the operator’s spent time on the selection of a 3D micro-CT image fragment for further modelling and
to excludes subjective judgments from this stage.
V International Conference on "Information Technology and Nanotechnology" (ITNT-2019) 166
Image Processing and Earth Remote Sensing
A S Kornilov, I V Safonov, A V Goncharova and I V Yakimchuk
Figure 5. The graph of normalized NIQE quality estimation on the synthetic 3D image slices, where
the noise level increases from central slice to top and bottom.
Figure 6. The graph of normalized NIQE metrics for slices of the real 3D image.
4. References
[1] Safonov I, Yakimchuk I and Abashkin V 2018 Algorithms for 3D particles characterization
using X-ray microtomography in proppant crush test J. of Imaging 4 134
[2] Koroteev D A, Dinariev O, Evseev N, Klemin D V, Safonov S, Gurpinar O M, Berg S, van
Kruijsdijk C, Myers M, Hathon L A, de Jong H and Armstrong R 2013 Application of digital
rock technology for chemical EOR screening Proc. of the SPE Enhanced Oil Recovery Conf.
165258
[3] Dinariev O and Evseev N 2016 Multiphase flow modeling with density functional method
Computational Geosciences 20 835-856
[4] Buzug T M 2008 Computed Tomography: From Photon Statistics to Modern Cone-Beam CT
(New York: Springer Science & Business Media)
[5] Smelkina N A, Kolsanov A V, Chaplygin S S, Zelter P M and Khramov A G 2017 Pulmonary
emphysema recognition by CT scan Computer Optics 41(5) 726-731 DOI: 10.18287/2412-
6179-2017-41-5-726-731
[6] Fabijańska A 2011 Variance filter for edge detection and edge-based image segmentation
MEMSTECH
[7] van Droogenbroeck M and Piérard S 2011 Object descriptors based on a list of rectangles:
method and algorithm Int. Symp. on Mathematical Morphology and Its Applications to Signal
and Image Processing
[8] Moorthy A K and Bovik A C 2010 A two-step framework for constructing blind image quality
indices IEEE Signal processing letters 17 513-516
[9] Mittal A, Moorthy A K and Bovik A C 2012 No-reference image quality assessment in the
spatial domain IEEE Transactions on Image Processing 21 4695-4708
V International Conference on "Information Technology and Nanotechnology" (ITNT-2019) 167
Image Processing and Earth Remote Sensing
A S Kornilov, I V Safonov, A V Goncharova and I V Yakimchuk
[10] Mittal A, Soundararajan R and Bovik A C 2013 Making a “completely blind” image quality
analyzer IEEE Signal Processing Letters 20 209-212
[11] Liu L, Hua Y, Zhao Q, Huang H and Bovik A C 2016 Blind image quality assessment by
relative gradient statistics and adaboosting neural network Signal Processing: Image
Communication 40 1-15
[12] Safonov I V, Kurilin I V, Rychagov M N and Tolstaya E V 2018 Adaptive Image Processing
Algorithms for Printing (Singapore: Springer)
[13] Soille P 2004 Morphological Image Analysis: Principles and Applications (New York: Springer
Science & Business Media)
[14] Asatryan D G 2017 Image blur estimation using gradient field analysis Computer Optics 41(6)
957-962 DOI: 10.18287/2412-6179-2017-41-6-957-962
[15] Sheikh H R, Sabir M F and Bovik A C 2006 A statistical evaluation of recent full reference
image quality assessment algorithms IEEE Transactions on Image Processing 15 3440-3451
[16] Srivastava A, Lee A B, Simoncelli E P and Zhu S C 2003 On advances in statistical modeling of
natural images J. of mathematical imaging and vision 18 17-33
[17] Ruderman D L and Bialek W 1994 Statistics of natural images: Scaling in the woods Advances
in neural information processing systems
[18] Crow F C 1984 Summed-area tables for texture mapping ACM SIGGRAPH computer graphics
18 207-212
[19] Phan T, Sohoni S, Chandler D M and Larson E C 2012 Performance-analysis-based acceleration
of image quality assessment SSIAI 81-84
V International Conference on "Information Technology and Nanotechnology" (ITNT-2019) 168