Adaptive-window indicator kriging A thresholding method for ...
Computers & Geosciences 54 (2013) 239–248
Contents lists available at SciVerse ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Adaptive-window indicator kriging: A thresholding method for computed tomography images of porous media A.N. Houston a,n, W. Otten a, P.C. Baveye a,b, S. Hapca a a b
SIMBIOS Centre, Kydd Building, Abertay University, 40 Bell Street, Dundee DD1 1HG, UK Laboratory of Soil and Water Engineering, Department of Civil and Environmental Engineering, Rensselaer Polytechnic Institute, 110 8th street, Troy, NY 12180, USA
a r t i c l e i n f o
abstract
Article history: Received 25 June 2012 Received in revised form 19 November 2012 Accepted 20 November 2012 Available online 29 November 2012
Image segmentation is a crucial step in understanding the structure of porous materials, subsequent analyses being profoundly dependent upon segmentation accuracy. Computed tomography images of naturally occurring heterogeneous materials such as soils are particularly challenging to segment reliably, due to the prevalence of partial volume effects, noise, and other artefacts induced during the image acquisition process. As a result, boundaries between classes of objects, typically pore versus solid in the case of soil, are difficult to identify. Indicator kriging can address these problems in a robust fashion but the computational cost can be very great under some circumstances; this is particularly true under conditions that occur regularly in images of soil. The kriging window size parameter is decisive in obtaining a good quality result at reasonable cost, but is difficult to estimate for an image exhibiting significant heterogeneity. This work demonstrates that, allowing the kriging window size to adapt locally throughout the image provides a very efficient solution. Moreover, it is shown that this can be achieved using a conceptually simple mechanism that involves negligible extra processing cost. The adaptive-window indicator kriging method described in this study achieves easily an order of magnitude improvement in computational performance over a fixed size window implementation without sacrificing quality. In addition, it is shown that, by improving the locality of estimation, the new method is robust when applied to soil images. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Local segmentation X-ray Microtomography Indicator kriging Soil
1. Introduction Over the last decades, numerous segmentation methods have been developed to transform greyscale images into binary formats, in order to allow the identification of meaningful patterns. Among the most popular approaches are the k-means clustering algorithm (Ridler and Calvard, 1978) and Otsu’s (1979) method. These are global methods, i.e., applied to whole images at once, without any particular attention to local heterogeneity. A common assumption in such global algorithms is that image elements (i.e., voxels in 3D or pixels in 2D) belong to two distinct classes determined by greyscale value. Yet, in a number of situations, greyscale values are not distinct enough across the whole image to be treated as two different classes, in which case global approaches are not able to produce sensible results. Such a situation occurs, for example when images of two phase systems are obtained at resolutions commensurate with the feature size of these systems. This causes significant ‘‘partial volume’’ effects, where image elements represent
n
Corresponding author. Tel.: þ44 7941 930062. E-mail addresses: [email protected], [email protected] (A.N. Houston). 0098-3004/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cageo.2012.11.016
some mixture of two phases resulting in greyscale values intermediate between those of the two phases. In addition, systems that are spatially heterogeneous at scales larger than the image resolution, or that are difficult to uniformly illuminate during the acquisition process, cannot always accommodate a global threshold. Some of the most severe problems along these lines arise in connection with three-dimensional (3D) images obtained by application of X-ray computed tomography (CT) to natural porous media, in particular soils. These systems are notoriously heterogeneous within an extremely wide range of scales (Ketcham and Carlson, 2001; Kaestner et al., 2008; Kravchenko et al., 2010; Pajor et al., 2010; Hapca et al., 2011). Over the last decade, major technological advances in imaging techniques, among which X-ray CT and micro-fluorescence spectroscopy, have afforded the possibility to investigate and reconstruct the 3D architecture of natural porous media at micro-metric resolutions, yet this is still not adequate to capture the size of many of the solid constituents (e.g., clay minerals, silt particles) of these systems. This leads to voxels with greyscale values intermediate between those of the solid constituents of the soils, or of pore space. In addition, the presence of organic matter, often fine-textured and distributed non-uniformly in the pores, as well as the presence of
240
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
high density minerals can cause further spatial heterogeneity in the CT scans. Not unexpectedly, under these conditions, global algorithms have been observed to lead to unsuitable thresholds ¨ (Iassonov et al., 2009; Baveye et al., 2010; Schluter et al., 2010; Wang et al., 2011). The difficulty of classifying an image in which many elements consist of a mixture of phases can be addressed by giving consideration to the local neighbourhood of each problematic image element, the use of this additional information within the local neighbourhood allowing more reliable classification decisions to be made. Among the existing local segmentation techniques available, the indicator kriging method, first proposed by Oh and Lindquist (1999), has been reported as a robust procedure for bi-level segmentation of 3D images of porous media (Sezgin and Sankur, 2004; Prodanovic´ et al., 2006; Iassonov et al., 2009; Wang et al., 2011). In a review of volume image segmentation methods applicable to geo-materials (Iassonov et al., 2009) it was ranked second of 14 methods in terms of quality, while being significantly less computationally demanding and also requiring less operator intervention than the highest quality method. Similarly, in a recent study by Wang et al. (2011), indicator kriging was found to have a better performance than four other methods, when tested on 3D simulated images of soils. Local segmentation algorithms have been devised based on all manner of mathematical treatments of the local neighbourhood, but indicator kriging differs from most in making use of a geostatistical model that is specifically deduced for each image. The method relies in the first instance on a partial classification, followed by a geostatistical method applied to the local neighbourhood of each unclassified image element to determine class membership. The original method developed by Oh and Lindquist (1999) uses a moving window of constant size to apply ordinary indicator kriging for each unclassified image element in the image. As indicated by the authors (Oh and Lindquist 1999), the assumption of a kriging window of constant and relatively small size is computationally acceptable for homogeneous images. However, for heterogeneous systems, a small window may not adequately represent the local spatial correlation, leading to classification failure. To overcome this problem in practice, the algorithm of Oh and Lindquist (1999) needs to be implemented with a large kriging window, which is exceedingly inefficient computationally when high definition images are segmented. This issue is particularly relevant for 3D images of soils. On-going technological advances make it likely that images to be segmented in the near future will involve hundreds of billions of voxels. In this case the computational efficiency of a segmentation method is a key determinant of its applicability. In this general context, the objective of this article is to extend the indicator kriging algorithm of Oh and Lindquist (1999) such that the size of the kriging window adapts locally to spatial conditions throughout an image. The proposed adaptive-window indicator kriging method is intended to improve the robustness and computational performance of the original method when applied to images of highly heterogeneous structures. This extension is described in Section 2 immediately following this introduction. Afterwards, a section detailing the application of the method to soil images is presented, followed by a quantitative comparison of the original (as developed by Oh and Lindquist, 1999) and the extended method, from both computational and segmentation quality perspectives.
2. Adaptive-window indicator kriging: tTheoretical approach Image thresholding by indicator kriging as described by Oh and Lindquist (1999) makes use of a geostatistical procedure to estimate
class membership of image elements by minimizing the spatial variance of indicator variable. The method relies on an initial partial classification of the image via a two level or hysteresis threshold to create a set of indicator data. A model semivariogram can then be determined for this empirical indicator data and the model is globally applied to resolve unclassified image elements. Specifically, for each unclassified image element, a system of constrained linear equations (ordinary kriging system) is formed to describe the local spatial pore-solid correlation structure, based on the semivariogram model. Each solved system of equations provides a set of interpolation weights, yielding a minimum variance estimate of class membership probability; membership is then assigned to the most probable class. In their 1999 publication, Oh and Lindquist propose a spherical kriging window with a radius comparable to the correlation length in the indicator data. According to this approach a kriging window of constant radius is applied throughout an image and all pre-classified image elements (indicator data points) lying within the window are used to form a kriging system. Henceforth this strategy is referred to as ‘‘fixed window’’ kriging, whose sole parameter is the radius of the window. Being relatively easy to
Fig. 1. (a) Illustration of a 2D circular window pattern of radius 3. Elements are labelled with their Euclidean distance from the central estimation point; border colour is assigned according to this distance. (b) Example of a partly classified image region with elements assigned colured borders as in Fig. 1(a). Pre-classified image elements (i.e., data points for kriging) are shown in black for pore and white for solid. Unclassified elements are shown in grey and the estimation point (the centre of the window) is marked with a cross. See Table 1 for a summary explanation of adaptive versus fixed window processing. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
Table 1 Adaptive window mechanism applied to the 2D problem illustrated in Fig. 1 showing, left to right, the progressive growth in number of data points with increasing window size. A fixed window of radius 3 would use the cumulative information in the right-most column irrespective of local conditions. An adaptive window advances column by column checking the cumulative total per class against targets. If the lower and upper targets per class were 4 and 8 respectively then evidence collection would end at radius 2 because the upper target is met for the pore class. If the upper target were raised to 16 then evidence collection would extend to radius 2O2 where both lower targets are met. Class
Pore Solid Unclassified Total
Tally
Additional Cumulative Additional Cumulative Additional Cumulative Additional Cumulative
Radius (image elements) 1
O2
2
O5
2O2
3
þ3 3 þ0 0 þ1 1 þ4 4
þ3 6 þ0 0 þ1 2 þ4 8
þ2 8 þ1 1 þ1 3 þ4 12
þ2 10 þ1 2 þ5 8 þ8 20
þ2 12 þ2 4 þ0 8 þ4 24
þ0 12 þ2 6 þ2 10 þ4 28
implement in a programming language, this is a common strategy for kriging digital images. However, as noted in the introduction, it fails to account for geo-materials (in particular soils) that are not well characterised by a single correlation length; addressing this problem requires local adaptation of the kriging window. There are potentially many ways in which an adaptive window might be implemented and choosing a good strategy depends largely upon the application. The method described following is specific to indicator kriging of digital images. The essential principles are, first: exploit the coherent lattice structure of the image to obtain data points organised into groups, all members of a group being equidistant from the estimation point (see Fig. 1a for a 2D example). Second, neither the window radius nor the number of data points is constrained; instead lower and upper target numbers of points per class of indicator datum are used as a guide. When suitable combined targets are achieved then classification by kriging proceeds. Third, recognising the possibility of classification failure, multiple attempts per estimation point are supported. Each additional classification attempt includes further data points in the kriging system i.e., the adaptive window procedure is resumed as necessary. The criteria for satisfactory evidence used in this study are that either both indicator classes must meet the lower target number of points, or the upper target must be met for at least one class. For example, setting lower and upper targets to 4 and 16 respectively means that search cannot end until at least 8 data points (4 or more per class) have been found. If only 3 points of one class are found then search continues until at least 16 points of the other class are found giving a total of at least 19 points. Bearing in mind that data points are found in groups by lag distance, the number of points per class may exceed the target, hence the size of kriging system is theoretically unlimited. In practice, targets are significantly exceeded only when an estimation point has few data points nearby; greater lag distances tend to include many data points. Where data points of all indicator classes are abundant, it is clear that low targets should be met quickly, resulting in a kriging system of modest size. Contrastingly, where a single class of indicator datum dominates, the upper target is more relevant hence kriging systems are typically larger. An example of the adaptive window mechanism as applied to a 2D problem is shown in Fig. 1, and Table 1.
3. Materials Five soil samples were used for method development and testing, for which images were collected using an X-Tek HMX CT
241
Table 2 Settings used for CT scanning and reconstruction of images used in the study. Sample
Image Image Image Image Image
1 2 3 4 5
Energy-flux set
Filter
Proj
Res.
kV
mA
Element
mm
#
mm
110 110 150 150 150
93 93 50 50 47
Al Al Al Al Al
0.25 0.25 2 2 2
3010 3010 1200 1200 1200
0.032 0.032 0.035 0.035 0.035
scanner and then tomographic reconstruction was achieved using CTPro version 1.2 (Nikon Metrology X-Tek Systems Ltd, Tring, UK). The image acquisition settings that were used in this study are listed in Table 2. Noise reduction was achieved as part of the filtered back-projection reconstruction algorithm. This was judged to suppress impulse noise to an acceptable degree without excessive distortion of the image intensity distribution and while preserving edge information. Non linear filtering is applied as part of the segmentation method, but follows pre-classification and does not affect greyscale images. A single region of interest, 2563 voxels in size, was selected from each reconstructed image. The voxels in this region of interest were linearly mapped from the 32-bit floating-point representation of the original tomographic reconstruction to an unsigned 8-bit representation to allow information (both images and histograms) to be displayed directly using simple computer graphics techniques. As the intensity ranges in the original images were very wide, a statistical method for outlier removal based on quartiles (Zar, 1999) was applied to determine a suitably smaller interval for mapping. Values within this interval were linearly mapped to the output range while values outside the interval were clamped to the limits of the output range. This allows the limited precision of 8-bit image elements to encode in sufficient detail the part of the image intensity histogram that is of interest.
4. Implementation and assessment of adaptive-window indicator kriging The overall segmentation procedure can be split into three successive steps: pre-classification, followed by spatial variance modelling and finally classification by kriging. The method was implemented using Microsoft Visual Cþþ 2008 and a software package for Windows is available on request from the first author. Computational performance was quantified by measuring the elapsed time for the kriging stage of processing only; the assessment was carried out under Windows7 using a single Intel Xeon E3-1240 (3.3 GHz) processor with 16 GB of 1600 MHz DDR3 memory. 4.1. Pre-classification of pore and solid sub-populations The objective of pre-classification is to obtain a representative sample of the spatial distribution of the two classes of image elements; this information therefore influences the overall segmentation result at a fundamental level. The pre-classification method proposed by Oh and Lindquist (1999) uses an interval defined by an upper and lower threshold of the image intensity. Two suggestions for automatically determining the interval from the intensity histogram of an image were given: (1) fitting of a bi-normal model using the expectation maximization algorithm; and (2) using an entropy function. It was found that neither of these methods works satisfactorily in the general case of soil CT images (Prodanovic´ et al., 2006; Iassonov et al., 2009). The image intensity histogram frequently
242
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
exhibits a skewed unimodal character (Fig. 2); hence, the automatic identification of two distinct modes is typically unstable and sometimes impossible. An attempt was made to resolve this problem by removing boundary image elements (identified by gradient magnitude) from the image histogram before pattern analysis. Unfortunately the histogram shape was not significantly improved and pattern analysis remained unreliable. Despite these problems, model fitting can be used to guide decisions made by a skilled operator. Given a Gaussian mixture model for example, thresholds may be deduced from the higher intensity inflection point of the modelled pore distribution versus the lower intensity inflection point of the modelled solid distribution. By weighing all the available evidence (perhaps including conflicting models) an operator can reach a more reasonable decision than might be obtained automatically.
4.2. Modelling spatial variance Spatial variance models (and hence kriging systems) based on the indicator semivariogram (as described in Goovaerts, 1997) were used in this study. The ‘‘classical’’ empirical semivariogram estimator proposed by Matheron (1963) was first calculated for the entire image. To constrain computational cost, the lag distance was limited to 12 image elements; this was found to be sufficient for all samples in this study. A theoretical indicator semivariogram model was then derived to fit the empirical data. When kriging, the semivariogram model is subject to mathematical constraints, such as the positive definiteness of the kriging matrix within the kriging system, hence only authorized models (Goovaerts, 1997) are fitted to the empirical semivariogram. Many authorized models exist and fitting these to empirical data is typically formulated as a residual error minimisation problem. In order to provide a good fit at short distances, a weighted least squares approach was applied using weights proportional to nðhÞg^ ðhÞ=g2 ðhÞ, where nðhÞis the number of pairs of image elements at lag distance h, g^ ðhÞis the empirical variogram and gðhÞ is the modelled variogram (Webster and Oliver, 2007). Furthermore, simulated annealing, as described in Lark and Papritz (2003) was used to search the parameter space of each candidate model. This approach was selected over other numerical methods because it is independent of algebraic form (i.e., no knowledge of model properties or ‘‘patterns’’ is required) and facilitates constraining of parameter values (such as positiveness of some model parameters), making it widely applicable. The types of indicator semivariogram model supported at the moment by the
newly developed software are limited to the simple isotropic forms of spherical, exponential and Gaussian models; the ‘‘nesting’’ of models is not supported, nor is any form of anisotropy.
4.3. Classification by indicator kriging Given an indicator semivariogram model and a partly classified image, classification can be completed by indicator kriging. The solution of a kriging system is a vector of weights that are summed by indicator class and standardised to yield the membership probability for the unclassified element. When both classes are equally probable (or very nearly so) the reliability of classification is questionable. In order to avoid classification under such circumstances, a ‘‘dead-zone’’ can be specified (e.g., the probability interval 0.499–0.501), which causes kriged probabilities falling within this range to be rejected. Another problem that might occur is the presence of negative weights suggesting a certain clustering pattern within the indicator data points used to formulate the kriging system (Stein, 2002). While negative weights are not disastrous, they do not have a valid physical meaning and should not be used in obtaining a classification. To overcome this problem Oh and Lindquist (1999) recommend the correction scheme described by Deutsch (1996), and this approach was adopted for the present study. In rare circumstances, however, this can lead to the rejection of all weights, in which case a valid solution cannot be obtained. For a fixed-window implementation there is no way to resolve a classification failure within a single pass through the image. In the original manuscript, Oh and Lindquist (1999) do not give any indication of how this problem should be addressed. In this study, complete classification using a fixed-window was achieved by applying indicator kriging iteratively in a number of passes over the image. The semivariogram model, as well as all previous classifications remain unchanged during each subsequent pass. This is therefore an approximation that feeds back the results of previous kriging systems into new systems and is consequently vulnerable to numerical error. For an adaptive window, it is relatively easy to collect further indicator data points by resuming the search. This allows another (larger) kriging system to be formulated, perhaps yielding a valid solution. Applied iteratively, this progressive solution strategy leads to an incremental ‘‘search and solve’’ algorithm, extending the adaptive window strategy described in Section 2. The consequence is that progressively larger window sizes are employed until either a valid solution is
Fig. 2. (a) Grey scale (unclassified) image slice mapped from tomographic reconstruction of X-ray attenuation coefficients and (b) the intensity histogram of the entire volume image.
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
obtained or some practical limitation (typically storage and/or processing requirements) is encountered. 4.4. Majority filtering Oh and Lindquist (1999) recommend additional processing in the form of a majority filter applied following pre-classification and again following kriging. This filter is applied only to pre-classified image elements (unclassified elements remain unchanged) and can change pre-classified image elements to the unclassified state. Thus, poorly correlated elements (whether pore or solid) resulting from pre-classification may subsequently be assigned class membership by kriging. The filter parameters used in this study were based on a 60% majority threshold within a neighbourhood of radius one (i.e., the nearest 26 voxels in 3D), as suggested by Oh and Lindquist (1999).
243
indicates topological structure, and increases in numerical value with the number of distinct objects present in the image, but decreases with the number of ‘‘holes’’ penetrating fully through objects. In a very general sense, the IMC may be considered to provide an estimate of the ‘‘ruggedness’’ of object surfaces, while the ITC measures ‘‘fragmentation’’ of objects. Another measure that is used to compare segmentation results is the pore space connectivity, defined here as the axial spanning volume fraction of the segmented pore phase. This gives the fraction of all pore voxels that are members of a cluster connecting opposite faces of a cubic region of interest. The local connectivity of each pore element is determined based on the six-connected neighbourhood criterion, i.e., only face–face voxel connections are considered. Finally, a measure of classification error is also estimated based on the intra-class variance calculated on the original greyscale image (Otsu, 1979, Zhang et al., 2008, Hapca et al., 2012).
4.5. Measures for evaluation of segmentation performance In order to evaluate the segmented images, functional measures including the 3D Minkowski functionals are used in addition to global statistics. The 3D Minkowski functionals have been successfully used in the past for the physical characterisation of soil structure (Vogel et al., 2010, Falconer et al., 2012). They are computed for each segmented image according to the prescription ¨ in Ohser and Mucklich (2000), the first two measures being the pore volume fraction (porosity) and surface area, while the other two measures being named integral mean curvature (IMC) and integral total curvature (ITC). The IMC gives an indication of the extent to which 3D spatial structure departs from a purely planar configuration (which has zero mean curvature) and is a signed quantity such that convexity is positive and concavity negative. The ITC
5. Results and discussion 5.1. Accuracy of semivariogram models For each image an empirical indicator semivariogram was calculated and modelled using the approach described previously. As illustrated in Fig. 3, all of the images in this study are in general well represented by the Gaussian model. Note that the sill of an indicator variogram should be no larger than 0.25 as this is the upper bound for the variance of a Bernoulli random variable (Papritz, 2009). In addition, when a Gaussian model is fitted to the indicator variogram, a positive nugget constraint is implemented, for the model to be consistent with the triangular inequality
Fig. 3. Semivariogram models are shown fitted against empirical data for each sample image. In each case the Gaussian authorised model has been selected as providing the best fit. Near the origin, the fit is limited by the requirement for a positive ‘‘nugget’’ (y-axis intercept) value. The model parameter values are also shown for each sample image.
244
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
estimation points) increase very slowly versus increased adaptive window parameter, while the standard deviation decreases slightly. This indicates that the quantity of evidence collected by using an adaptive window is on average relatively stable across the image as a whole. This is further supported by the statistics on kriging matrix rank which exhibit very similar patterns. An assessment of kriging system solutions, in terms of the occurrence of negative weights, is presented in Table 4. The mean fraction of negative weights measured over all kriging systems is seen to increase steadily with matrix rank. Although the fractions for adaptive versus fixed window processing overlap, the larger matrix rank of fixed windows means that a greater absolute number of negative weights occur in fixed window processing. This suggests that an adaptive window avoids processing redundant information, and so is more efficient.
requirement as indicated by Matheron (1989). Although the Gaussian model proved to provide the best fit, some limitations in the accuracy of fit at small lag distances are apparent, which might be due to extremely small values of the empirical indicator semivariogram at these lag distances. It therefore seems that a nested model formulation (which is not covered in the present implementation) might be suited to accurately model both short and long range behaviour. 5.2. Robustness of kriging systems For the adaptive-window approach the parameters controlling the size of kriging systems are the lower and upper target number of indicator data points (i.e., pre-classified image elements) per class of indicator datum. In this study, lower targets from 1 to 10 inclusive were evaluated with upper targets remaining constant at 16. For each sample image, the 10 segmentations show the expected progression of kriging matrix rank and other statistics as seen in Table 3. For comparison, fixed window results are also shown in Table 3, but due to limited processing capacity, only two segmentations per sample image were evaluated; these used constant window radii of 3 and 4, respectively. In all the images tested, both mean and maximum radii (measured over all
5.3. Effect of adaptive sampling on functional measures Changing adaptive sampling parameters so as to increase the effective kriging window size (henceforth EKWS) has the effect, on average, of increasing both matrix rank and spatial extent of kriging systems. As illustrated in Table 5, this increase causes both surface area and mean curvature to decrease monotonically, in other words
Table 3 Statistical properties of kriging systems including computational performance are shown for both fixed and adaptive window parameter ranges. The effective window size is determined per kriging system as the farthest lag distance incorporated into the system. All of the values presented in this table are calculated over all kriging systems formed within a particular sample image and parameter combination. Note that for a fixed window, the small variability seen in lag distance is caused by a local absence of data points in some cases. Image Sample
Image 1
Image 2
Image 3
Image 4
Image 5
Window parameter
Effective window size
Kriging matrix rank
Adaptive low target
Fixed radius
Mean
StDev
1 2 3 ^ 10 – –
– – –
2.01 2.21 2.32
0.58 0.56 0.54
5.92 6.16 6.16
– 3 4
2.62 2.99 3.92
0.50 0.11 0.20
1 2 3 ^ 10 – –
– – –
1.72 1.94 2.06
– 3 4
1 2 3 ^ 10 – –
Max
Mean
Elapsed time
Mean rate
All kriging attempts (seconds)
(classifica-tions per second)
StDev
Max
10.31 14.63 17.17
6.42 6.64 6.29
36 40 40
1.98 2.54 2.84
2.5Eþ 06 1.9Eþ 06 1.7Eþ 06
6.40 3.00 4.00
26.31 47.08 106.76
4.86 20.14 41.09
47 123 240
5.30 32.00 352.24
9.3Eþ 05 1.5Eþ 05 1.4Eþ 04
0.50 0.47 0.46
5.20 5.39 5.39
9.48 14.25 16.98
6.39 6.79 6.51
38 40 41
3.09 4.45 5.79
1.7Eþ 06 1.2Eþ 06 9.3Eþ 05
2.39 2.97 3.91
0.37 0.17 0.26
5.83 3.00 4.00
27.14 56.84 129.09
5.40 20.00 38.41
49 123 257
12.71 116.69 1269.64
4.3Eþ 05 4.6Eþ 04 4.3Eþ 03
– – –
1.96 2.16 2.27
0.57 0.55 0.53
6.00 6.08 6.08
10.03 14.35 16.90
6.31 6.52 6.21
37 38 39
3.12 4.51 5.59
2.0E þ 06 1.4Eþ 06 1.1Eþ 06
– 3 4
2.59 2.99 3.93
0.47 0.10 0.18
6.40 3.00 4.00
26.50 47.78 107.77
5.02 20.60 40.84
48 123 257
11.69 75.29 815.54
5.3Eþ 05 8.2Eþ 04 7.5Eþ 03
1 2 3 ^ 10 – –
– – –
2.56 2.75 2.85
0.96 0.95 0.95
11.87 12.04 12.04
10.81 14.79 17.15
6.34 6.38 5.83
35 37 39
3.07 4.04 4.56
2.2Eþ 06 1.7Eþ 06 1.5Eþ 06
– 3 4
3.12 3.00 3.93
0.98 0.06 0.15
12.12 3.00 4.00
25.17 40.33 88.03
4.28 18.91 40.03
48 123 257
7.52 35.90 351.53
9.2Eþ 05 1.9Eþ 05 2.0E þ 04
1 2 3 ^ 10 – –
– – –
1.78 2.01 2.13
0.56 0.52 0.50
5.10 5.48 5.48
9.64 14.34 17.01
6.39 6.71 6.37
38 39 40
2.55 3.69 4.71
2.1Eþ 06 1.4Eþ 06 1.1Eþ 06
– 3 4
2.46 2.97 3.90
0.42 0.18 0.29
5.92 3.00 4.00
26.78 53.45 120.88
5.35 21.33 41.47
49 123 257
9.95 83.30 898.45
5.3Eþ 05 6.3Eþ 04 5.9Eþ 03
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
245
Table 4 The number of classification failures following the first kriging pass are shown for both fixed and adaptive window parameter ranges. Window parameter
Image sample Image 1
Adaptive low target 1 2 3 ^ 10 – –
Image 2
Image 3
Image 4
Image 5
Single Fixed pass radius failures – 0 – 0 – 0
Single Negative weights pass mean proportion failures 29.0% 0 30.6% 0 31.7% 0
Single Negative weights pass mean proportion failures 33.1% 0 35.2% 0 37.9% 0
Single Negative weights pass mean proportion failures 28.0% 0 29.0% 0 30.0% 0
Single Negative weights pass mean proportion failures 39.2% 0 42.3% 0 44.2% 0
Negative weights mean proportion 44.0% 47.1% 48.3%
– 3 4
38.4% 37.7% 45.8%
45.5% 44.9% 56.7%
36.6% 37.1% 45.5%
46.9% 42.6% 51.1%
49.5% 49.1% 50.0%
0 8595 436
0 376 7
0 8365 251
0 373941 154163
0 1680 16
Table 5 Functional measures of segmented images are shown for both fixed and adaptive window parameter ranges. Image sample
Image 1
Image 2
Image 3
Image 4
Image 5
Window parameter
Minkowski functionals Connectivity
Intra-class variance
3205.19 2635.14 1835.17
87.97% 87.35% 87.52%
1142.8 1146.47 1149.14
168.72 175.08 150.45
2791.43 2801.49 2137.67
87.72% 88.68% 82.29%
1158.62 1164.16 1185.83
27.40 27.01 26.66
422.38 387.71 361.98
17578.70 3472.20 4790.83
97.92% 97.96% 98.06%
1194.13 1197.26 1201.02
0.2877 0.2852 0.2823
24.91 23.49 22.41
304.09 278.08 244.75
6423.09 11685.40 12427.10
98.21% 98.44% 98.55%
1218.80 1239.81 1248.32
– – –
0.2771 0.2763 0.2758
24.39 23.95 23.55
357.11 343.51 327.67
13097.30 12936.20 14213.00
97.14% 97.25% 97.22%
1344.56 1348.99 1352.75
– 3 4
0.2749 0.2750 0.2711
22.63 22.18 20.53
302.68 305.04 265.20
10820.30 11754.30 9318.88
97.12% 97.49% 97.44%
1367.09 1382.18 1412.70
– – –
0.2732 0.2731 0.2732
13.85 13.85 13.76
181.95 181.41 177.95
1314.92 1215.32 1270.42
96.29% 96.31% 96.35%
883.96 886.40 888.62
– 3 4
0.2724 0.2731 0.2701
13.17 13.26 12.41
166.85 178.24 161.35
1078.11 529.78 638.92
96.46% 96.45% 96.67%
897.15 895.48 917.08
– – –
0.2516 0.2483 0.2475
20.89 20.30 19.96
380.71 334.97 308.12
37767.70 21451.40 12862.10
95.23% 95.53% 95.65%
1839.43 1846.57 1856.52
– 3 4
0.2454 0.2452 0.2413
18.74 18.26 17.06
254.88 251.17 193.76
14673.90 12057.90 5979.13
95.56% 95.79% 95.91%
1885.74 1899.21 1918.66
Adaptive low target
Fixed radius
Volume
Surface
IMC
1 2 3 ^ 10 – –
– – –
0.0981 0.0977 0.0974
10.02 9.83 9.68
202.70 192.40 183.01
– 3 4
0.0969 0.0972 0.0949
9.30 9.22 8.54
1 2 3 ^ 10 – –
– – –
0.2897 0.2891 0.2892
– 3 4
1 2 3 ^ 10 – – 1 2 3 ^ 10 – – 1 2 3 ^ 10 – –
the structure of the pore–solid interface becomes somewhat smooth (at the finest length scales) as the average size of the kriging system grows. This is consistent with the interpretation of kriging as an interpolation method. This interpretation is also supported by the fact that, for the soil samples used in this study, the porosity generally tends to decrease with larger values of EKWS as the dominant solid phase begins to mask the relatively sparse evidence of pore phase (Figs. 4 and 5). Given that pre-classified pore elements are relatively rare and clustered together, it can be argued that
ITC
extending the radius tends to favour solid elements. Over short ranges of EKWS there is no simple pattern of porosity change, which is a consequence of the different correlation structures of the images. The total curvature shows obvious change with EKWS in every case indicating that the number and topology of individual pore clusters is changing. In two cases the total curvature changes sign, a positive to negative transition indicating that the connection of previously isolated pore clusters contains ‘‘bridges’’ of solid material penetrating through the pore space. A negative to positive
246
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
Fig. 4. (a) Pre-classified image slice obtained by applying a hysteresis threshold to the grey scale image slice shown in Fig. 2. Pore elements are shown in black, solid elements in white and unclassified elements in purple. (b) Complete classification by adaptive-window indicator kriging of the image slice shown in (a).
Fig. 5. Volume rendering of 3D sample image 1 before and after segmentation. On the left is the grey-scale image, on the right is the binary image segmented using adaptive window indicator kriging (parameters used were lower and upper targets of 4 and 16, respectively).
transition indicates that previously existing solid ‘‘bridges’’ have been removed leaving a more simply connected pore space. The connectivity measure appears to be stable with respect to changes in the EKWS, with less than 1% fluctuations. This suggests that, despite the expected dominance of solid phase over longer ranges, there are perhaps cases where the evidence of pores at these ranges is sufficient to tip the balance in favour of pore classification. Finally, in all the tested images, an increase in the EKWS produced an increase in the intra-class variance measure (Table 5), which is explained by the increased degree of spatial averaging that tends to occur in larger kriging systems. If ICV is interpreted as a measure of classification error, the smallest EKWS therefore produces better segmentation results. 5.4. Fixed versus adaptive kriging window: eEvaluation of performance The adaptive window produces a substantial reduction in average matrix rank compared to a fixed window, greatly improving computational performance. The processing time for indicator kriging using an adaptive window ranged between 2 and 13 s
(Table 3), while for a fixed window the time lies between 32 s and 22 min. In the fixed window case at least two kriging passes were required to completely classify each image. Table 4 shows the number of classification failures at the end of the first kriging pass, for sample image 4 this amounts to more than 2% of the image for a fixed window of radius 3, falling to less than 1% at radius 4. In order to avoid such problems, Iassonov et al. (2009), implemented indicator kriging with a fixed window of radius 5. The samples used in this study however require a window radius exceeding 5, in particular sample 4 requires a window of size 12 to achieve complete classification in one pass. Such a large window size implies a processing cost approaching 27 times (123/43) that of a radius 4 fixed window and this was judged to be impractical due to its very inefficient use of computing resources. In general, the two methods produced very similar results in terms of functional measures with a difference in the estimated pore volume of less than 0.1% and a difference in connectivity measure of less than 1%. Compared to the adaptive window method, the fixed window approach has the tendency to reduce the pore volume, pore surface and the IMC measures, while the
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
ITC and connectivity measures have larger values (Table 5), a trend that is consistent throughout the five samples. With the increase in the EKWS, the adaptive window results seem to converge toward fixed window results, this trend confirming the robustness of the adaptive window method. Additionally, for all samples the calculated intra-class variance is consistently smaller under the adaptive window approach; again this favours the newly developed method. The adaptive window approach developed in this study, where window size is controlled by the lower and upper target number of data points per class, appears well suited to indicator kriging of digital images. In wider applications of kriging it is possible that data points may number only a few thousand-. This permits a ‘‘top-down’’ approach in which the entire data is searched per estimation point, first finding all points within some cutoff radius. According to Webster and Oliver (2007), between 4 and 20 local data points are required to obtain a good quality estimate by kriging and so, pragmatically, up to 20 data points lying near the estimation point might be selected. If the cutoff radius does not admit sufficient data points this mechanism fails, hence the radius may be set generously, increasing the average workload. This evidence selection strategy is weak in the case of applying indicator kriging to 3D images of soil. First, the total of data plus estimation points numbers in the millions to billions of which estimation points typically represent 10–30%. This scale of problem means that ‘‘bottom-up’’ search (i.e., gradual outward expansion from each estimation point) is essential for acceptable performance. Second, a fixed cutoff radius is simply not efficient given the heterogeneity typical of soil images. Considering the number of estimation points, the work involved in acquiring data points within a fixed radius tends to become a limiting factor. Third, selecting a specific maximum number of data points assumes the existence of an objective criterion e.g., that all points can be uniquely ordered by distance from the estimation point. For the exact gridding of digital images this is untrue as many equidistant data points are typically found per lag distance. Finally, the loss of subtle data variations caused by transforming continuous-valued samples into indicators, in conjunction with significant clustering, means that extra care is due in the selection of evidence. Constraining data points to some specified number is questionable when roughly equal proportions of solid versus pore classes are found locally and when adding a few extra data points might alter the balance enough to yield a different classification. It might be considered that these difficulties could be addressed by re-sampling the densely gridded image data into a sparse form through the use of some de-clustering method. While this strategy might be useful for some applications, it would defeat objectives of subsequent analyses applied to the segmented images. The high density of information in a CT image of soil is essential to accurately locating the solid-pore interface so that valid estimates of soil characteristics such as surface area, curvature and connectivity can be obtained.
6. Conclusions Motivated by the computational challenges imposed by the 3D images of soils, the present study proposes an extension of the commonly used indicator-kriging algorithm of Oh and Lindquist (1999), which allows for the size of the kriging window to adapt according to the local spatial conditions in an image. It was found that the proposed extension produces segmentation results very close to the fixed window method while greatly reducing computational cost. Through improved locality of estimation, the adaptive-window approach is clearly advantageous for images containing imbalanced proportions of material phases and/or a
247
high degree of spatial heterogeneity. In addition, the adaptive window capability of greatly extending the spatial radius, as required by the local image conditions, means that a valid solution can be achieved via kriging. A fixed window may be appropriate for relatively homogeneous spatial distributions of material but performs quite poorly for soil. Large fixed window sizes are prohibitively costly whereas the likelihood of classification failure is high for a small fixed window size. If the latter problem is remedied by further processing (such as using previously kriged classifications as input for new kriging systems) then the consistent quality of results might be difficult to guarantee. A number of potential improvements of the method in the future include increasing the efficiency of empirical semivariogram measurement and the accuracy of semivariogram modelling. In the current implementation, the collection of the empirical semivariogram is a limiting factor on overall computational performance and using a Monte-Carlo sub-sampling approach could help in addressing this problem. Improving the accuracy of semivariogram modelling is more challenging, as it appears that a nested model formulation may be required, implying that an automated fitting method is difficult to achieve. Another possible goal is support for anisotropic models; this however implies a local analysis of spatial variance, which introduces either localised model-boundary effects or potentially difficult interpolation problems. Further ways to improve the adaptive window indicator kriging method will be investigated in future studies.
References Baveye, P.C., Laba, M., Otten, W., Bouckaert, L., Dello Sterpaio, P., Goswami, R.R., Grinev, D., Houston, A.N., Hu, Y., Liu, J., Mooney, S., Pajor, R., Sleutel, S., Tarquis, A., Wang, W., Wei, Q., Sezgin, M., 2010. Observer-dependent variability of the thresholding step in the quantitative analysis of soil images and X-ray microtomography data. Geoderma 157, 51–63. Deutsch, C.V., 1996. Corrections for negative weights in ordinary kriging. Computers & Geosciences 22 (7), 765–773. Falconer, R., Houston, A.N., Otten, W., Baveye, P., 2012. Emergent behaviour of fungal dynamics: influence of soil architecture and moisture distribution. Soil Science 177 (2), 111–119. Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Oxford University Press, New York 496pp. Hapca, S.M., Wang, Z.X., Otten, W., Wilson, C., Baveye, P.C., 2011. Automated statistical method to align 2D chemical maps with 3D X-ray computed microtomographic images of soils. Geoderma 164, 146–154. Hapca, S.M., Houston, A.N., Otten, W., Baveye, P.C., 2012. New objective segmentation method based on minimizing locally the intra-class variance of greyscale images. Vadose Zone Journal, submitted. Iassonov, P., Gebrenegus, T., Tuller, M., 2009. Segmentation of X-ray computed tomography images of porous materials: a crucial step for characterization and quantitative analysis of pore structures. Water Resources Research 45, 9, art. no. W09415. Kaestner, A., Lehmann, E., Stampanoni, M., 2008. Imaging and image processing in porous media research. Advances in Water Resources 31, 1174–1187. Ketcham, R.A., Carlson, W.D., 2001. Acquisition, optimization and interpretation of X-ray computed tomographic imagery: applications to the geosciences. Computers & Geosciences 27, 381–400. Kravchenko, A., Falconer, R., Grinev, D., Otten, W., 2010. Fungal colonization in soils of contrasting managements: modelling fungal growth in 3D pore volumes of undisturbed soil samples. Ecological Applications 21, 1202–1210. Lark, M.R., Papritz, A., 2003. Fitting a linear model of coregionalization for soil properties using simulating annealing. Geoderma 115, 245–260. Matheron, G., 1963. Principles of geostatistics. Economic Geology 58, 1246–1266. Oh, W., Lindquist, W., 1999. Image thresholding by indicator kriging. IEEE Transactions on Pattern Analysis and Machine Intelligence 21, 590–602. ¨ Ohser, J., Mucklich, F., 2000. Statistical Analysis of Microstructures in Materials Science. Wiley, Chichester 375 pp.. Otsu, N., 1979. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man and Cybernetics 9, 62–66. Pajor, R., Falconer, R., Hapca, S., Otten, W., 2010. Modelling and quantifying the effect of heterogeneity in soil physical conditions on fungal growth. Biogeosciences 7, 3731–3740. Prodanovic´, M., Lindquist, W.B., Seright, R.S., 2006. Porous structure and fluid partitioning in polyethylene cores from 3D X-ray microtomographic imaging. Journal of Colloid Interface Science 298, 282–297.
248
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
Ridler, T.W., Calvard, S., 1978. Picture thresholding using an iterative selection method. IEEE Transactions on Systems, Man and Cybernetics SMC-8, 630–632. ¨ Schluter, S., Weller, U., Vogel, H.-J., 2010. Segmentation of X-ray microtomography images of soil using gradient masks. Computers & Geosciences 36, 1246–1251. Sezgin, M., Sankur, B., 2004. Survey over image thresholding techniques and quantitative performance evaluation. Journal of Electronic Imaging 13, 146–168. Stein, M.L., 2002. The screening effect in kriging. The Annals of Statistics 30 (1), 298–323. ¨ Vogel, H.J., Weller, U., Schluter, S., 2010. Quantification of soil structure based on Minkowski functions. Computers & Geosciences 36, 1236–1245.
Wang, W., Kravchenko, A.N., Smucker, A.J.M., Rivers, M.L., 2011. Comparison of image segmentation methods in simulated 2D and 3D microtomographic images of soil aggregates. Geoderma 162, 231–241. Webster, R., Oliver, M., 2007. Geostatistics for Environmental Scientists. Statistics in Practice Series, second ed Wiley 315 pp. Zar, J.H., 1999. Biostatistical Analysis. Prentice Hall 663 pp. Zhang, H., Fritts, J.E., Goldman, S.A., 2008. Image segmentation evaluation: a survey of unsupervised methods. Computer Vision and Image Understanding 110, 260–280.
Contents lists available at SciVerse ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Adaptive-window indicator kriging: A thresholding method for computed tomography images of porous media A.N. Houston a,n, W. Otten a, P.C. Baveye a,b, S. Hapca a a b
SIMBIOS Centre, Kydd Building, Abertay University, 40 Bell Street, Dundee DD1 1HG, UK Laboratory of Soil and Water Engineering, Department of Civil and Environmental Engineering, Rensselaer Polytechnic Institute, 110 8th street, Troy, NY 12180, USA
a r t i c l e i n f o
abstract
Article history: Received 25 June 2012 Received in revised form 19 November 2012 Accepted 20 November 2012 Available online 29 November 2012
Image segmentation is a crucial step in understanding the structure of porous materials, subsequent analyses being profoundly dependent upon segmentation accuracy. Computed tomography images of naturally occurring heterogeneous materials such as soils are particularly challenging to segment reliably, due to the prevalence of partial volume effects, noise, and other artefacts induced during the image acquisition process. As a result, boundaries between classes of objects, typically pore versus solid in the case of soil, are difficult to identify. Indicator kriging can address these problems in a robust fashion but the computational cost can be very great under some circumstances; this is particularly true under conditions that occur regularly in images of soil. The kriging window size parameter is decisive in obtaining a good quality result at reasonable cost, but is difficult to estimate for an image exhibiting significant heterogeneity. This work demonstrates that, allowing the kriging window size to adapt locally throughout the image provides a very efficient solution. Moreover, it is shown that this can be achieved using a conceptually simple mechanism that involves negligible extra processing cost. The adaptive-window indicator kriging method described in this study achieves easily an order of magnitude improvement in computational performance over a fixed size window implementation without sacrificing quality. In addition, it is shown that, by improving the locality of estimation, the new method is robust when applied to soil images. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Local segmentation X-ray Microtomography Indicator kriging Soil
1. Introduction Over the last decades, numerous segmentation methods have been developed to transform greyscale images into binary formats, in order to allow the identification of meaningful patterns. Among the most popular approaches are the k-means clustering algorithm (Ridler and Calvard, 1978) and Otsu’s (1979) method. These are global methods, i.e., applied to whole images at once, without any particular attention to local heterogeneity. A common assumption in such global algorithms is that image elements (i.e., voxels in 3D or pixels in 2D) belong to two distinct classes determined by greyscale value. Yet, in a number of situations, greyscale values are not distinct enough across the whole image to be treated as two different classes, in which case global approaches are not able to produce sensible results. Such a situation occurs, for example when images of two phase systems are obtained at resolutions commensurate with the feature size of these systems. This causes significant ‘‘partial volume’’ effects, where image elements represent
n
Corresponding author. Tel.: þ44 7941 930062. E-mail addresses: [email protected], [email protected] (A.N. Houston). 0098-3004/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cageo.2012.11.016
some mixture of two phases resulting in greyscale values intermediate between those of the two phases. In addition, systems that are spatially heterogeneous at scales larger than the image resolution, or that are difficult to uniformly illuminate during the acquisition process, cannot always accommodate a global threshold. Some of the most severe problems along these lines arise in connection with three-dimensional (3D) images obtained by application of X-ray computed tomography (CT) to natural porous media, in particular soils. These systems are notoriously heterogeneous within an extremely wide range of scales (Ketcham and Carlson, 2001; Kaestner et al., 2008; Kravchenko et al., 2010; Pajor et al., 2010; Hapca et al., 2011). Over the last decade, major technological advances in imaging techniques, among which X-ray CT and micro-fluorescence spectroscopy, have afforded the possibility to investigate and reconstruct the 3D architecture of natural porous media at micro-metric resolutions, yet this is still not adequate to capture the size of many of the solid constituents (e.g., clay minerals, silt particles) of these systems. This leads to voxels with greyscale values intermediate between those of the solid constituents of the soils, or of pore space. In addition, the presence of organic matter, often fine-textured and distributed non-uniformly in the pores, as well as the presence of
240
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
high density minerals can cause further spatial heterogeneity in the CT scans. Not unexpectedly, under these conditions, global algorithms have been observed to lead to unsuitable thresholds ¨ (Iassonov et al., 2009; Baveye et al., 2010; Schluter et al., 2010; Wang et al., 2011). The difficulty of classifying an image in which many elements consist of a mixture of phases can be addressed by giving consideration to the local neighbourhood of each problematic image element, the use of this additional information within the local neighbourhood allowing more reliable classification decisions to be made. Among the existing local segmentation techniques available, the indicator kriging method, first proposed by Oh and Lindquist (1999), has been reported as a robust procedure for bi-level segmentation of 3D images of porous media (Sezgin and Sankur, 2004; Prodanovic´ et al., 2006; Iassonov et al., 2009; Wang et al., 2011). In a review of volume image segmentation methods applicable to geo-materials (Iassonov et al., 2009) it was ranked second of 14 methods in terms of quality, while being significantly less computationally demanding and also requiring less operator intervention than the highest quality method. Similarly, in a recent study by Wang et al. (2011), indicator kriging was found to have a better performance than four other methods, when tested on 3D simulated images of soils. Local segmentation algorithms have been devised based on all manner of mathematical treatments of the local neighbourhood, but indicator kriging differs from most in making use of a geostatistical model that is specifically deduced for each image. The method relies in the first instance on a partial classification, followed by a geostatistical method applied to the local neighbourhood of each unclassified image element to determine class membership. The original method developed by Oh and Lindquist (1999) uses a moving window of constant size to apply ordinary indicator kriging for each unclassified image element in the image. As indicated by the authors (Oh and Lindquist 1999), the assumption of a kriging window of constant and relatively small size is computationally acceptable for homogeneous images. However, for heterogeneous systems, a small window may not adequately represent the local spatial correlation, leading to classification failure. To overcome this problem in practice, the algorithm of Oh and Lindquist (1999) needs to be implemented with a large kriging window, which is exceedingly inefficient computationally when high definition images are segmented. This issue is particularly relevant for 3D images of soils. On-going technological advances make it likely that images to be segmented in the near future will involve hundreds of billions of voxels. In this case the computational efficiency of a segmentation method is a key determinant of its applicability. In this general context, the objective of this article is to extend the indicator kriging algorithm of Oh and Lindquist (1999) such that the size of the kriging window adapts locally to spatial conditions throughout an image. The proposed adaptive-window indicator kriging method is intended to improve the robustness and computational performance of the original method when applied to images of highly heterogeneous structures. This extension is described in Section 2 immediately following this introduction. Afterwards, a section detailing the application of the method to soil images is presented, followed by a quantitative comparison of the original (as developed by Oh and Lindquist, 1999) and the extended method, from both computational and segmentation quality perspectives.
2. Adaptive-window indicator kriging: tTheoretical approach Image thresholding by indicator kriging as described by Oh and Lindquist (1999) makes use of a geostatistical procedure to estimate
class membership of image elements by minimizing the spatial variance of indicator variable. The method relies on an initial partial classification of the image via a two level or hysteresis threshold to create a set of indicator data. A model semivariogram can then be determined for this empirical indicator data and the model is globally applied to resolve unclassified image elements. Specifically, for each unclassified image element, a system of constrained linear equations (ordinary kriging system) is formed to describe the local spatial pore-solid correlation structure, based on the semivariogram model. Each solved system of equations provides a set of interpolation weights, yielding a minimum variance estimate of class membership probability; membership is then assigned to the most probable class. In their 1999 publication, Oh and Lindquist propose a spherical kriging window with a radius comparable to the correlation length in the indicator data. According to this approach a kriging window of constant radius is applied throughout an image and all pre-classified image elements (indicator data points) lying within the window are used to form a kriging system. Henceforth this strategy is referred to as ‘‘fixed window’’ kriging, whose sole parameter is the radius of the window. Being relatively easy to
Fig. 1. (a) Illustration of a 2D circular window pattern of radius 3. Elements are labelled with their Euclidean distance from the central estimation point; border colour is assigned according to this distance. (b) Example of a partly classified image region with elements assigned colured borders as in Fig. 1(a). Pre-classified image elements (i.e., data points for kriging) are shown in black for pore and white for solid. Unclassified elements are shown in grey and the estimation point (the centre of the window) is marked with a cross. See Table 1 for a summary explanation of adaptive versus fixed window processing. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
Table 1 Adaptive window mechanism applied to the 2D problem illustrated in Fig. 1 showing, left to right, the progressive growth in number of data points with increasing window size. A fixed window of radius 3 would use the cumulative information in the right-most column irrespective of local conditions. An adaptive window advances column by column checking the cumulative total per class against targets. If the lower and upper targets per class were 4 and 8 respectively then evidence collection would end at radius 2 because the upper target is met for the pore class. If the upper target were raised to 16 then evidence collection would extend to radius 2O2 where both lower targets are met. Class
Pore Solid Unclassified Total
Tally
Additional Cumulative Additional Cumulative Additional Cumulative Additional Cumulative
Radius (image elements) 1
O2
2
O5
2O2
3
þ3 3 þ0 0 þ1 1 þ4 4
þ3 6 þ0 0 þ1 2 þ4 8
þ2 8 þ1 1 þ1 3 þ4 12
þ2 10 þ1 2 þ5 8 þ8 20
þ2 12 þ2 4 þ0 8 þ4 24
þ0 12 þ2 6 þ2 10 þ4 28
implement in a programming language, this is a common strategy for kriging digital images. However, as noted in the introduction, it fails to account for geo-materials (in particular soils) that are not well characterised by a single correlation length; addressing this problem requires local adaptation of the kriging window. There are potentially many ways in which an adaptive window might be implemented and choosing a good strategy depends largely upon the application. The method described following is specific to indicator kriging of digital images. The essential principles are, first: exploit the coherent lattice structure of the image to obtain data points organised into groups, all members of a group being equidistant from the estimation point (see Fig. 1a for a 2D example). Second, neither the window radius nor the number of data points is constrained; instead lower and upper target numbers of points per class of indicator datum are used as a guide. When suitable combined targets are achieved then classification by kriging proceeds. Third, recognising the possibility of classification failure, multiple attempts per estimation point are supported. Each additional classification attempt includes further data points in the kriging system i.e., the adaptive window procedure is resumed as necessary. The criteria for satisfactory evidence used in this study are that either both indicator classes must meet the lower target number of points, or the upper target must be met for at least one class. For example, setting lower and upper targets to 4 and 16 respectively means that search cannot end until at least 8 data points (4 or more per class) have been found. If only 3 points of one class are found then search continues until at least 16 points of the other class are found giving a total of at least 19 points. Bearing in mind that data points are found in groups by lag distance, the number of points per class may exceed the target, hence the size of kriging system is theoretically unlimited. In practice, targets are significantly exceeded only when an estimation point has few data points nearby; greater lag distances tend to include many data points. Where data points of all indicator classes are abundant, it is clear that low targets should be met quickly, resulting in a kriging system of modest size. Contrastingly, where a single class of indicator datum dominates, the upper target is more relevant hence kriging systems are typically larger. An example of the adaptive window mechanism as applied to a 2D problem is shown in Fig. 1, and Table 1.
3. Materials Five soil samples were used for method development and testing, for which images were collected using an X-Tek HMX CT
241
Table 2 Settings used for CT scanning and reconstruction of images used in the study. Sample
Image Image Image Image Image
1 2 3 4 5
Energy-flux set
Filter
Proj
Res.
kV
mA
Element
mm
#
mm
110 110 150 150 150
93 93 50 50 47
Al Al Al Al Al
0.25 0.25 2 2 2
3010 3010 1200 1200 1200
0.032 0.032 0.035 0.035 0.035
scanner and then tomographic reconstruction was achieved using CTPro version 1.2 (Nikon Metrology X-Tek Systems Ltd, Tring, UK). The image acquisition settings that were used in this study are listed in Table 2. Noise reduction was achieved as part of the filtered back-projection reconstruction algorithm. This was judged to suppress impulse noise to an acceptable degree without excessive distortion of the image intensity distribution and while preserving edge information. Non linear filtering is applied as part of the segmentation method, but follows pre-classification and does not affect greyscale images. A single region of interest, 2563 voxels in size, was selected from each reconstructed image. The voxels in this region of interest were linearly mapped from the 32-bit floating-point representation of the original tomographic reconstruction to an unsigned 8-bit representation to allow information (both images and histograms) to be displayed directly using simple computer graphics techniques. As the intensity ranges in the original images were very wide, a statistical method for outlier removal based on quartiles (Zar, 1999) was applied to determine a suitably smaller interval for mapping. Values within this interval were linearly mapped to the output range while values outside the interval were clamped to the limits of the output range. This allows the limited precision of 8-bit image elements to encode in sufficient detail the part of the image intensity histogram that is of interest.
4. Implementation and assessment of adaptive-window indicator kriging The overall segmentation procedure can be split into three successive steps: pre-classification, followed by spatial variance modelling and finally classification by kriging. The method was implemented using Microsoft Visual Cþþ 2008 and a software package for Windows is available on request from the first author. Computational performance was quantified by measuring the elapsed time for the kriging stage of processing only; the assessment was carried out under Windows7 using a single Intel Xeon E3-1240 (3.3 GHz) processor with 16 GB of 1600 MHz DDR3 memory. 4.1. Pre-classification of pore and solid sub-populations The objective of pre-classification is to obtain a representative sample of the spatial distribution of the two classes of image elements; this information therefore influences the overall segmentation result at a fundamental level. The pre-classification method proposed by Oh and Lindquist (1999) uses an interval defined by an upper and lower threshold of the image intensity. Two suggestions for automatically determining the interval from the intensity histogram of an image were given: (1) fitting of a bi-normal model using the expectation maximization algorithm; and (2) using an entropy function. It was found that neither of these methods works satisfactorily in the general case of soil CT images (Prodanovic´ et al., 2006; Iassonov et al., 2009). The image intensity histogram frequently
242
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
exhibits a skewed unimodal character (Fig. 2); hence, the automatic identification of two distinct modes is typically unstable and sometimes impossible. An attempt was made to resolve this problem by removing boundary image elements (identified by gradient magnitude) from the image histogram before pattern analysis. Unfortunately the histogram shape was not significantly improved and pattern analysis remained unreliable. Despite these problems, model fitting can be used to guide decisions made by a skilled operator. Given a Gaussian mixture model for example, thresholds may be deduced from the higher intensity inflection point of the modelled pore distribution versus the lower intensity inflection point of the modelled solid distribution. By weighing all the available evidence (perhaps including conflicting models) an operator can reach a more reasonable decision than might be obtained automatically.
4.2. Modelling spatial variance Spatial variance models (and hence kriging systems) based on the indicator semivariogram (as described in Goovaerts, 1997) were used in this study. The ‘‘classical’’ empirical semivariogram estimator proposed by Matheron (1963) was first calculated for the entire image. To constrain computational cost, the lag distance was limited to 12 image elements; this was found to be sufficient for all samples in this study. A theoretical indicator semivariogram model was then derived to fit the empirical data. When kriging, the semivariogram model is subject to mathematical constraints, such as the positive definiteness of the kriging matrix within the kriging system, hence only authorized models (Goovaerts, 1997) are fitted to the empirical semivariogram. Many authorized models exist and fitting these to empirical data is typically formulated as a residual error minimisation problem. In order to provide a good fit at short distances, a weighted least squares approach was applied using weights proportional to nðhÞg^ ðhÞ=g2 ðhÞ, where nðhÞis the number of pairs of image elements at lag distance h, g^ ðhÞis the empirical variogram and gðhÞ is the modelled variogram (Webster and Oliver, 2007). Furthermore, simulated annealing, as described in Lark and Papritz (2003) was used to search the parameter space of each candidate model. This approach was selected over other numerical methods because it is independent of algebraic form (i.e., no knowledge of model properties or ‘‘patterns’’ is required) and facilitates constraining of parameter values (such as positiveness of some model parameters), making it widely applicable. The types of indicator semivariogram model supported at the moment by the
newly developed software are limited to the simple isotropic forms of spherical, exponential and Gaussian models; the ‘‘nesting’’ of models is not supported, nor is any form of anisotropy.
4.3. Classification by indicator kriging Given an indicator semivariogram model and a partly classified image, classification can be completed by indicator kriging. The solution of a kriging system is a vector of weights that are summed by indicator class and standardised to yield the membership probability for the unclassified element. When both classes are equally probable (or very nearly so) the reliability of classification is questionable. In order to avoid classification under such circumstances, a ‘‘dead-zone’’ can be specified (e.g., the probability interval 0.499–0.501), which causes kriged probabilities falling within this range to be rejected. Another problem that might occur is the presence of negative weights suggesting a certain clustering pattern within the indicator data points used to formulate the kriging system (Stein, 2002). While negative weights are not disastrous, they do not have a valid physical meaning and should not be used in obtaining a classification. To overcome this problem Oh and Lindquist (1999) recommend the correction scheme described by Deutsch (1996), and this approach was adopted for the present study. In rare circumstances, however, this can lead to the rejection of all weights, in which case a valid solution cannot be obtained. For a fixed-window implementation there is no way to resolve a classification failure within a single pass through the image. In the original manuscript, Oh and Lindquist (1999) do not give any indication of how this problem should be addressed. In this study, complete classification using a fixed-window was achieved by applying indicator kriging iteratively in a number of passes over the image. The semivariogram model, as well as all previous classifications remain unchanged during each subsequent pass. This is therefore an approximation that feeds back the results of previous kriging systems into new systems and is consequently vulnerable to numerical error. For an adaptive window, it is relatively easy to collect further indicator data points by resuming the search. This allows another (larger) kriging system to be formulated, perhaps yielding a valid solution. Applied iteratively, this progressive solution strategy leads to an incremental ‘‘search and solve’’ algorithm, extending the adaptive window strategy described in Section 2. The consequence is that progressively larger window sizes are employed until either a valid solution is
Fig. 2. (a) Grey scale (unclassified) image slice mapped from tomographic reconstruction of X-ray attenuation coefficients and (b) the intensity histogram of the entire volume image.
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
obtained or some practical limitation (typically storage and/or processing requirements) is encountered. 4.4. Majority filtering Oh and Lindquist (1999) recommend additional processing in the form of a majority filter applied following pre-classification and again following kriging. This filter is applied only to pre-classified image elements (unclassified elements remain unchanged) and can change pre-classified image elements to the unclassified state. Thus, poorly correlated elements (whether pore or solid) resulting from pre-classification may subsequently be assigned class membership by kriging. The filter parameters used in this study were based on a 60% majority threshold within a neighbourhood of radius one (i.e., the nearest 26 voxels in 3D), as suggested by Oh and Lindquist (1999).
243
indicates topological structure, and increases in numerical value with the number of distinct objects present in the image, but decreases with the number of ‘‘holes’’ penetrating fully through objects. In a very general sense, the IMC may be considered to provide an estimate of the ‘‘ruggedness’’ of object surfaces, while the ITC measures ‘‘fragmentation’’ of objects. Another measure that is used to compare segmentation results is the pore space connectivity, defined here as the axial spanning volume fraction of the segmented pore phase. This gives the fraction of all pore voxels that are members of a cluster connecting opposite faces of a cubic region of interest. The local connectivity of each pore element is determined based on the six-connected neighbourhood criterion, i.e., only face–face voxel connections are considered. Finally, a measure of classification error is also estimated based on the intra-class variance calculated on the original greyscale image (Otsu, 1979, Zhang et al., 2008, Hapca et al., 2012).
4.5. Measures for evaluation of segmentation performance In order to evaluate the segmented images, functional measures including the 3D Minkowski functionals are used in addition to global statistics. The 3D Minkowski functionals have been successfully used in the past for the physical characterisation of soil structure (Vogel et al., 2010, Falconer et al., 2012). They are computed for each segmented image according to the prescription ¨ in Ohser and Mucklich (2000), the first two measures being the pore volume fraction (porosity) and surface area, while the other two measures being named integral mean curvature (IMC) and integral total curvature (ITC). The IMC gives an indication of the extent to which 3D spatial structure departs from a purely planar configuration (which has zero mean curvature) and is a signed quantity such that convexity is positive and concavity negative. The ITC
5. Results and discussion 5.1. Accuracy of semivariogram models For each image an empirical indicator semivariogram was calculated and modelled using the approach described previously. As illustrated in Fig. 3, all of the images in this study are in general well represented by the Gaussian model. Note that the sill of an indicator variogram should be no larger than 0.25 as this is the upper bound for the variance of a Bernoulli random variable (Papritz, 2009). In addition, when a Gaussian model is fitted to the indicator variogram, a positive nugget constraint is implemented, for the model to be consistent with the triangular inequality
Fig. 3. Semivariogram models are shown fitted against empirical data for each sample image. In each case the Gaussian authorised model has been selected as providing the best fit. Near the origin, the fit is limited by the requirement for a positive ‘‘nugget’’ (y-axis intercept) value. The model parameter values are also shown for each sample image.
244
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
estimation points) increase very slowly versus increased adaptive window parameter, while the standard deviation decreases slightly. This indicates that the quantity of evidence collected by using an adaptive window is on average relatively stable across the image as a whole. This is further supported by the statistics on kriging matrix rank which exhibit very similar patterns. An assessment of kriging system solutions, in terms of the occurrence of negative weights, is presented in Table 4. The mean fraction of negative weights measured over all kriging systems is seen to increase steadily with matrix rank. Although the fractions for adaptive versus fixed window processing overlap, the larger matrix rank of fixed windows means that a greater absolute number of negative weights occur in fixed window processing. This suggests that an adaptive window avoids processing redundant information, and so is more efficient.
requirement as indicated by Matheron (1989). Although the Gaussian model proved to provide the best fit, some limitations in the accuracy of fit at small lag distances are apparent, which might be due to extremely small values of the empirical indicator semivariogram at these lag distances. It therefore seems that a nested model formulation (which is not covered in the present implementation) might be suited to accurately model both short and long range behaviour. 5.2. Robustness of kriging systems For the adaptive-window approach the parameters controlling the size of kriging systems are the lower and upper target number of indicator data points (i.e., pre-classified image elements) per class of indicator datum. In this study, lower targets from 1 to 10 inclusive were evaluated with upper targets remaining constant at 16. For each sample image, the 10 segmentations show the expected progression of kriging matrix rank and other statistics as seen in Table 3. For comparison, fixed window results are also shown in Table 3, but due to limited processing capacity, only two segmentations per sample image were evaluated; these used constant window radii of 3 and 4, respectively. In all the images tested, both mean and maximum radii (measured over all
5.3. Effect of adaptive sampling on functional measures Changing adaptive sampling parameters so as to increase the effective kriging window size (henceforth EKWS) has the effect, on average, of increasing both matrix rank and spatial extent of kriging systems. As illustrated in Table 5, this increase causes both surface area and mean curvature to decrease monotonically, in other words
Table 3 Statistical properties of kriging systems including computational performance are shown for both fixed and adaptive window parameter ranges. The effective window size is determined per kriging system as the farthest lag distance incorporated into the system. All of the values presented in this table are calculated over all kriging systems formed within a particular sample image and parameter combination. Note that for a fixed window, the small variability seen in lag distance is caused by a local absence of data points in some cases. Image Sample
Image 1
Image 2
Image 3
Image 4
Image 5
Window parameter
Effective window size
Kriging matrix rank
Adaptive low target
Fixed radius
Mean
StDev
1 2 3 ^ 10 – –
– – –
2.01 2.21 2.32
0.58 0.56 0.54
5.92 6.16 6.16
– 3 4
2.62 2.99 3.92
0.50 0.11 0.20
1 2 3 ^ 10 – –
– – –
1.72 1.94 2.06
– 3 4
1 2 3 ^ 10 – –
Max
Mean
Elapsed time
Mean rate
All kriging attempts (seconds)
(classifica-tions per second)
StDev
Max
10.31 14.63 17.17
6.42 6.64 6.29
36 40 40
1.98 2.54 2.84
2.5Eþ 06 1.9Eþ 06 1.7Eþ 06
6.40 3.00 4.00
26.31 47.08 106.76
4.86 20.14 41.09
47 123 240
5.30 32.00 352.24
9.3Eþ 05 1.5Eþ 05 1.4Eþ 04
0.50 0.47 0.46
5.20 5.39 5.39
9.48 14.25 16.98
6.39 6.79 6.51
38 40 41
3.09 4.45 5.79
1.7Eþ 06 1.2Eþ 06 9.3Eþ 05
2.39 2.97 3.91
0.37 0.17 0.26
5.83 3.00 4.00
27.14 56.84 129.09
5.40 20.00 38.41
49 123 257
12.71 116.69 1269.64
4.3Eþ 05 4.6Eþ 04 4.3Eþ 03
– – –
1.96 2.16 2.27
0.57 0.55 0.53
6.00 6.08 6.08
10.03 14.35 16.90
6.31 6.52 6.21
37 38 39
3.12 4.51 5.59
2.0E þ 06 1.4Eþ 06 1.1Eþ 06
– 3 4
2.59 2.99 3.93
0.47 0.10 0.18
6.40 3.00 4.00
26.50 47.78 107.77
5.02 20.60 40.84
48 123 257
11.69 75.29 815.54
5.3Eþ 05 8.2Eþ 04 7.5Eþ 03
1 2 3 ^ 10 – –
– – –
2.56 2.75 2.85
0.96 0.95 0.95
11.87 12.04 12.04
10.81 14.79 17.15
6.34 6.38 5.83
35 37 39
3.07 4.04 4.56
2.2Eþ 06 1.7Eþ 06 1.5Eþ 06
– 3 4
3.12 3.00 3.93
0.98 0.06 0.15
12.12 3.00 4.00
25.17 40.33 88.03
4.28 18.91 40.03
48 123 257
7.52 35.90 351.53
9.2Eþ 05 1.9Eþ 05 2.0E þ 04
1 2 3 ^ 10 – –
– – –
1.78 2.01 2.13
0.56 0.52 0.50
5.10 5.48 5.48
9.64 14.34 17.01
6.39 6.71 6.37
38 39 40
2.55 3.69 4.71
2.1Eþ 06 1.4Eþ 06 1.1Eþ 06
– 3 4
2.46 2.97 3.90
0.42 0.18 0.29
5.92 3.00 4.00
26.78 53.45 120.88
5.35 21.33 41.47
49 123 257
9.95 83.30 898.45
5.3Eþ 05 6.3Eþ 04 5.9Eþ 03
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
245
Table 4 The number of classification failures following the first kriging pass are shown for both fixed and adaptive window parameter ranges. Window parameter
Image sample Image 1
Adaptive low target 1 2 3 ^ 10 – –
Image 2
Image 3
Image 4
Image 5
Single Fixed pass radius failures – 0 – 0 – 0
Single Negative weights pass mean proportion failures 29.0% 0 30.6% 0 31.7% 0
Single Negative weights pass mean proportion failures 33.1% 0 35.2% 0 37.9% 0
Single Negative weights pass mean proportion failures 28.0% 0 29.0% 0 30.0% 0
Single Negative weights pass mean proportion failures 39.2% 0 42.3% 0 44.2% 0
Negative weights mean proportion 44.0% 47.1% 48.3%
– 3 4
38.4% 37.7% 45.8%
45.5% 44.9% 56.7%
36.6% 37.1% 45.5%
46.9% 42.6% 51.1%
49.5% 49.1% 50.0%
0 8595 436
0 376 7
0 8365 251
0 373941 154163
0 1680 16
Table 5 Functional measures of segmented images are shown for both fixed and adaptive window parameter ranges. Image sample
Image 1
Image 2
Image 3
Image 4
Image 5
Window parameter
Minkowski functionals Connectivity
Intra-class variance
3205.19 2635.14 1835.17
87.97% 87.35% 87.52%
1142.8 1146.47 1149.14
168.72 175.08 150.45
2791.43 2801.49 2137.67
87.72% 88.68% 82.29%
1158.62 1164.16 1185.83
27.40 27.01 26.66
422.38 387.71 361.98
17578.70 3472.20 4790.83
97.92% 97.96% 98.06%
1194.13 1197.26 1201.02
0.2877 0.2852 0.2823
24.91 23.49 22.41
304.09 278.08 244.75
6423.09 11685.40 12427.10
98.21% 98.44% 98.55%
1218.80 1239.81 1248.32
– – –
0.2771 0.2763 0.2758
24.39 23.95 23.55
357.11 343.51 327.67
13097.30 12936.20 14213.00
97.14% 97.25% 97.22%
1344.56 1348.99 1352.75
– 3 4
0.2749 0.2750 0.2711
22.63 22.18 20.53
302.68 305.04 265.20
10820.30 11754.30 9318.88
97.12% 97.49% 97.44%
1367.09 1382.18 1412.70
– – –
0.2732 0.2731 0.2732
13.85 13.85 13.76
181.95 181.41 177.95
1314.92 1215.32 1270.42
96.29% 96.31% 96.35%
883.96 886.40 888.62
– 3 4
0.2724 0.2731 0.2701
13.17 13.26 12.41
166.85 178.24 161.35
1078.11 529.78 638.92
96.46% 96.45% 96.67%
897.15 895.48 917.08
– – –
0.2516 0.2483 0.2475
20.89 20.30 19.96
380.71 334.97 308.12
37767.70 21451.40 12862.10
95.23% 95.53% 95.65%
1839.43 1846.57 1856.52
– 3 4
0.2454 0.2452 0.2413
18.74 18.26 17.06
254.88 251.17 193.76
14673.90 12057.90 5979.13
95.56% 95.79% 95.91%
1885.74 1899.21 1918.66
Adaptive low target
Fixed radius
Volume
Surface
IMC
1 2 3 ^ 10 – –
– – –
0.0981 0.0977 0.0974
10.02 9.83 9.68
202.70 192.40 183.01
– 3 4
0.0969 0.0972 0.0949
9.30 9.22 8.54
1 2 3 ^ 10 – –
– – –
0.2897 0.2891 0.2892
– 3 4
1 2 3 ^ 10 – – 1 2 3 ^ 10 – – 1 2 3 ^ 10 – –
the structure of the pore–solid interface becomes somewhat smooth (at the finest length scales) as the average size of the kriging system grows. This is consistent with the interpretation of kriging as an interpolation method. This interpretation is also supported by the fact that, for the soil samples used in this study, the porosity generally tends to decrease with larger values of EKWS as the dominant solid phase begins to mask the relatively sparse evidence of pore phase (Figs. 4 and 5). Given that pre-classified pore elements are relatively rare and clustered together, it can be argued that
ITC
extending the radius tends to favour solid elements. Over short ranges of EKWS there is no simple pattern of porosity change, which is a consequence of the different correlation structures of the images. The total curvature shows obvious change with EKWS in every case indicating that the number and topology of individual pore clusters is changing. In two cases the total curvature changes sign, a positive to negative transition indicating that the connection of previously isolated pore clusters contains ‘‘bridges’’ of solid material penetrating through the pore space. A negative to positive
246
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
Fig. 4. (a) Pre-classified image slice obtained by applying a hysteresis threshold to the grey scale image slice shown in Fig. 2. Pore elements are shown in black, solid elements in white and unclassified elements in purple. (b) Complete classification by adaptive-window indicator kriging of the image slice shown in (a).
Fig. 5. Volume rendering of 3D sample image 1 before and after segmentation. On the left is the grey-scale image, on the right is the binary image segmented using adaptive window indicator kriging (parameters used were lower and upper targets of 4 and 16, respectively).
transition indicates that previously existing solid ‘‘bridges’’ have been removed leaving a more simply connected pore space. The connectivity measure appears to be stable with respect to changes in the EKWS, with less than 1% fluctuations. This suggests that, despite the expected dominance of solid phase over longer ranges, there are perhaps cases where the evidence of pores at these ranges is sufficient to tip the balance in favour of pore classification. Finally, in all the tested images, an increase in the EKWS produced an increase in the intra-class variance measure (Table 5), which is explained by the increased degree of spatial averaging that tends to occur in larger kriging systems. If ICV is interpreted as a measure of classification error, the smallest EKWS therefore produces better segmentation results. 5.4. Fixed versus adaptive kriging window: eEvaluation of performance The adaptive window produces a substantial reduction in average matrix rank compared to a fixed window, greatly improving computational performance. The processing time for indicator kriging using an adaptive window ranged between 2 and 13 s
(Table 3), while for a fixed window the time lies between 32 s and 22 min. In the fixed window case at least two kriging passes were required to completely classify each image. Table 4 shows the number of classification failures at the end of the first kriging pass, for sample image 4 this amounts to more than 2% of the image for a fixed window of radius 3, falling to less than 1% at radius 4. In order to avoid such problems, Iassonov et al. (2009), implemented indicator kriging with a fixed window of radius 5. The samples used in this study however require a window radius exceeding 5, in particular sample 4 requires a window of size 12 to achieve complete classification in one pass. Such a large window size implies a processing cost approaching 27 times (123/43) that of a radius 4 fixed window and this was judged to be impractical due to its very inefficient use of computing resources. In general, the two methods produced very similar results in terms of functional measures with a difference in the estimated pore volume of less than 0.1% and a difference in connectivity measure of less than 1%. Compared to the adaptive window method, the fixed window approach has the tendency to reduce the pore volume, pore surface and the IMC measures, while the
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
ITC and connectivity measures have larger values (Table 5), a trend that is consistent throughout the five samples. With the increase in the EKWS, the adaptive window results seem to converge toward fixed window results, this trend confirming the robustness of the adaptive window method. Additionally, for all samples the calculated intra-class variance is consistently smaller under the adaptive window approach; again this favours the newly developed method. The adaptive window approach developed in this study, where window size is controlled by the lower and upper target number of data points per class, appears well suited to indicator kriging of digital images. In wider applications of kriging it is possible that data points may number only a few thousand-. This permits a ‘‘top-down’’ approach in which the entire data is searched per estimation point, first finding all points within some cutoff radius. According to Webster and Oliver (2007), between 4 and 20 local data points are required to obtain a good quality estimate by kriging and so, pragmatically, up to 20 data points lying near the estimation point might be selected. If the cutoff radius does not admit sufficient data points this mechanism fails, hence the radius may be set generously, increasing the average workload. This evidence selection strategy is weak in the case of applying indicator kriging to 3D images of soil. First, the total of data plus estimation points numbers in the millions to billions of which estimation points typically represent 10–30%. This scale of problem means that ‘‘bottom-up’’ search (i.e., gradual outward expansion from each estimation point) is essential for acceptable performance. Second, a fixed cutoff radius is simply not efficient given the heterogeneity typical of soil images. Considering the number of estimation points, the work involved in acquiring data points within a fixed radius tends to become a limiting factor. Third, selecting a specific maximum number of data points assumes the existence of an objective criterion e.g., that all points can be uniquely ordered by distance from the estimation point. For the exact gridding of digital images this is untrue as many equidistant data points are typically found per lag distance. Finally, the loss of subtle data variations caused by transforming continuous-valued samples into indicators, in conjunction with significant clustering, means that extra care is due in the selection of evidence. Constraining data points to some specified number is questionable when roughly equal proportions of solid versus pore classes are found locally and when adding a few extra data points might alter the balance enough to yield a different classification. It might be considered that these difficulties could be addressed by re-sampling the densely gridded image data into a sparse form through the use of some de-clustering method. While this strategy might be useful for some applications, it would defeat objectives of subsequent analyses applied to the segmented images. The high density of information in a CT image of soil is essential to accurately locating the solid-pore interface so that valid estimates of soil characteristics such as surface area, curvature and connectivity can be obtained.
6. Conclusions Motivated by the computational challenges imposed by the 3D images of soils, the present study proposes an extension of the commonly used indicator-kriging algorithm of Oh and Lindquist (1999), which allows for the size of the kriging window to adapt according to the local spatial conditions in an image. It was found that the proposed extension produces segmentation results very close to the fixed window method while greatly reducing computational cost. Through improved locality of estimation, the adaptive-window approach is clearly advantageous for images containing imbalanced proportions of material phases and/or a
247
high degree of spatial heterogeneity. In addition, the adaptive window capability of greatly extending the spatial radius, as required by the local image conditions, means that a valid solution can be achieved via kriging. A fixed window may be appropriate for relatively homogeneous spatial distributions of material but performs quite poorly for soil. Large fixed window sizes are prohibitively costly whereas the likelihood of classification failure is high for a small fixed window size. If the latter problem is remedied by further processing (such as using previously kriged classifications as input for new kriging systems) then the consistent quality of results might be difficult to guarantee. A number of potential improvements of the method in the future include increasing the efficiency of empirical semivariogram measurement and the accuracy of semivariogram modelling. In the current implementation, the collection of the empirical semivariogram is a limiting factor on overall computational performance and using a Monte-Carlo sub-sampling approach could help in addressing this problem. Improving the accuracy of semivariogram modelling is more challenging, as it appears that a nested model formulation may be required, implying that an automated fitting method is difficult to achieve. Another possible goal is support for anisotropic models; this however implies a local analysis of spatial variance, which introduces either localised model-boundary effects or potentially difficult interpolation problems. Further ways to improve the adaptive window indicator kriging method will be investigated in future studies.
References Baveye, P.C., Laba, M., Otten, W., Bouckaert, L., Dello Sterpaio, P., Goswami, R.R., Grinev, D., Houston, A.N., Hu, Y., Liu, J., Mooney, S., Pajor, R., Sleutel, S., Tarquis, A., Wang, W., Wei, Q., Sezgin, M., 2010. Observer-dependent variability of the thresholding step in the quantitative analysis of soil images and X-ray microtomography data. Geoderma 157, 51–63. Deutsch, C.V., 1996. Corrections for negative weights in ordinary kriging. Computers & Geosciences 22 (7), 765–773. Falconer, R., Houston, A.N., Otten, W., Baveye, P., 2012. Emergent behaviour of fungal dynamics: influence of soil architecture and moisture distribution. Soil Science 177 (2), 111–119. Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Oxford University Press, New York 496pp. Hapca, S.M., Wang, Z.X., Otten, W., Wilson, C., Baveye, P.C., 2011. Automated statistical method to align 2D chemical maps with 3D X-ray computed microtomographic images of soils. Geoderma 164, 146–154. Hapca, S.M., Houston, A.N., Otten, W., Baveye, P.C., 2012. New objective segmentation method based on minimizing locally the intra-class variance of greyscale images. Vadose Zone Journal, submitted. Iassonov, P., Gebrenegus, T., Tuller, M., 2009. Segmentation of X-ray computed tomography images of porous materials: a crucial step for characterization and quantitative analysis of pore structures. Water Resources Research 45, 9, art. no. W09415. Kaestner, A., Lehmann, E., Stampanoni, M., 2008. Imaging and image processing in porous media research. Advances in Water Resources 31, 1174–1187. Ketcham, R.A., Carlson, W.D., 2001. Acquisition, optimization and interpretation of X-ray computed tomographic imagery: applications to the geosciences. Computers & Geosciences 27, 381–400. Kravchenko, A., Falconer, R., Grinev, D., Otten, W., 2010. Fungal colonization in soils of contrasting managements: modelling fungal growth in 3D pore volumes of undisturbed soil samples. Ecological Applications 21, 1202–1210. Lark, M.R., Papritz, A., 2003. Fitting a linear model of coregionalization for soil properties using simulating annealing. Geoderma 115, 245–260. Matheron, G., 1963. Principles of geostatistics. Economic Geology 58, 1246–1266. Oh, W., Lindquist, W., 1999. Image thresholding by indicator kriging. IEEE Transactions on Pattern Analysis and Machine Intelligence 21, 590–602. ¨ Ohser, J., Mucklich, F., 2000. Statistical Analysis of Microstructures in Materials Science. Wiley, Chichester 375 pp.. Otsu, N., 1979. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man and Cybernetics 9, 62–66. Pajor, R., Falconer, R., Hapca, S., Otten, W., 2010. Modelling and quantifying the effect of heterogeneity in soil physical conditions on fungal growth. Biogeosciences 7, 3731–3740. Prodanovic´, M., Lindquist, W.B., Seright, R.S., 2006. Porous structure and fluid partitioning in polyethylene cores from 3D X-ray microtomographic imaging. Journal of Colloid Interface Science 298, 282–297.
248
A.N. Houston et al. / Computers & Geosciences 54 (2013) 239–248
Ridler, T.W., Calvard, S., 1978. Picture thresholding using an iterative selection method. IEEE Transactions on Systems, Man and Cybernetics SMC-8, 630–632. ¨ Schluter, S., Weller, U., Vogel, H.-J., 2010. Segmentation of X-ray microtomography images of soil using gradient masks. Computers & Geosciences 36, 1246–1251. Sezgin, M., Sankur, B., 2004. Survey over image thresholding techniques and quantitative performance evaluation. Journal of Electronic Imaging 13, 146–168. Stein, M.L., 2002. The screening effect in kriging. The Annals of Statistics 30 (1), 298–323. ¨ Vogel, H.J., Weller, U., Schluter, S., 2010. Quantification of soil structure based on Minkowski functions. Computers & Geosciences 36, 1236–1245.
Wang, W., Kravchenko, A.N., Smucker, A.J.M., Rivers, M.L., 2011. Comparison of image segmentation methods in simulated 2D and 3D microtomographic images of soil aggregates. Geoderma 162, 231–241. Webster, R., Oliver, M., 2007. Geostatistics for Environmental Scientists. Statistics in Practice Series, second ed Wiley 315 pp. Zar, J.H., 1999. Biostatistical Analysis. Prentice Hall 663 pp. Zhang, H., Fritts, J.E., Goldman, S.A., 2008. Image segmentation evaluation: a survey of unsupervised methods. Computer Vision and Image Understanding 110, 260–280.
