Detecting grain boundaries in deformed rocks using a cellular
This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/copyright
Author's personal copy Computers & Geosciences 42 (2012) 136–142
Contents lists available at SciVerse ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Detecting grain boundaries in deformed rocks using a cellular automata approach Pece V. Gorsevski a,n, Charles M. Onasch a, John R. Farver b, Xinyue Ye c a
School of Earth, Environment and Society, Bowling Green State University, Bowling Green, OH 43403, USA Department of Geology, Bowling Green State University, Bowling Green, OH 43403, USA c Center for Regional Development and School of Earth, Environment and Society, Bowling Green State University, Bowling Green, OH 43403, USA b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 22 June 2011 Received in revised form 31 August 2011 Accepted 2 September 2011 Available online 4 October 2011
Cellular automata (CA) are widely used in geospatial dynamic modeling and image processing. Here, we explore the application of two-dimensional cellular automata to the problem of grain boundary detection and extraction in digital images of thin sections from deformed rocks. The automated extraction of boundaries, which contain rich sources of information such as shape, orientation, and spatial distribution of grains, involves a CA Moore’s neighborhood-based rules approach. The Moore’s neighborhood is a 3 3 matrix that is used for changing states by comparing differences between a central pixel and its neighbors. In this dynamic approach, the future state of a pixel depends upon its current state and that of its neighbors. The rules that are defined determine the future state of each cell (i.e., on or off) while the number of iterations to simulate boundaries detection are specified by the user. Each iteration outputs different detection scenarios of grain boundaries that can be evaluated and assessed for accuracy. For a deformed quartz arenite, an r2 of 0.724 was obtained by comparing manually digitized grains to model derived grains. The value of this proposed method is compared against a traditional manual digitization approach and a recent GIS-based method developed for this purpose by Li et al. (2007). & 2011 Elsevier Ltd. All rights reserved.
Keywords: Cellular automata Grain boundary Thin section Edge detection GIS
1. Introduction The examination of thin sections with a polarizing light microscope is a standard way of characterizing rocks. Analyzing the mineralogy, size, shape, and spatial distribution of grains within a rock provides clues about processes such as the crystallization, metamorphism, deformation, and weathering (Fueten and Goodchild, 2001; Choudhury et al., 2006; Li et al., 2007). Although textural properties of rock microstructures are often sufficient for a broad geologic description such as the nature and origin of a rock, understanding of processes that control textural evolution requires objective and accurate measurements of the size, shape, orientation, and position of the mineral grains in a rock followed by an appropriate quantitative analysis. The recognition of grain boundaries is regarded as a vital first step in textural analysis of thin sections of rocks. Typically, grain boundaries are extracted either by manual digitization (Fabbri, 1984; Higgins, 2000; Jerram et al., 2003; Boorman et al., 2004) or by the use of automated extraction using image processing or n
Corresponding author. Tel.: þ1 419 372 7201. E-mail addresses: [email protected] (P.V. Gorsevski), [email protected] (C.M. Onasch), [email protected] (J.R. Farver), [email protected] (X. Ye). 0098-3004/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2011.09.008
GIS-based techniques (Goodchild and Fueten, 1998; van den Berg et al., 2002; Perring et al., 2004; Barraud, 2006; Choudhury et al., 2006; Li et al., 2007; Obara, 2007a,b). Some of the major shortcomings of the manual digitization arise from the labor intensive nature of this technique and lack of reproducibility between different analysts. On the other hand, the automated extraction of grain boundaries is aimed at replacing current, labor intensive manual procedures with consistent, replicable, and robust procedures allowing grains to be identified objectively using computerassisted analysis of images. Some image-specific, computer-assisted techniques that use automated extraction of grain boundaries include microscopic image segmentation using watershed methods (Luumbreras and Serrat, 1996; Zhou et al., 2004; Barraud, 2006), constrained automated seeded region growing (CASRG) algorithm for automated grain edge detection (Choudhury et al., 2006), application of neural networks (Fueten and Mason, 2007), genetic programming (Ross et al., 2001), gradient filtering (Heilbronner, 2000), image classification techniques (Li et al., 2007), and edge detection using sequential image processing routines based on Gaussian smoothing, gradient calculation, thresholding, and image skeletonization (Goodchild and Fueten, 1998). Some of these techniques take advantage of a fully integrated geographic information system (GIS) environment (e.g., Li et al., 2007) that
Author's personal copy P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
makes it relatively simple to implement a range of analytical capabilities such as data management, digital image processing, statistical analysis, spatial modeling, and geoprocessing. The GIS framework allows different algorithms and methodologies to be applied depending on the question of interest and problem that needs to be solved. In particular, GIS integration of microscopederived datasets provides several advantages including access to the functionality of the GIS for display and file management, reduced development time, enhanced usability, and built-in dynamic modeling capabilities. This kind of toolset has great potential for rapid identification and quantification of mineral microstructures that can be compiled in a spatial grain boundary database and quantitatively evaluated for understanding crucial processes that control grain evolution. Edge detection is a fundamental part of image processing, particularly in the areas of detection and extraction of features such as grain boundaries. There are many edge detection (boundary detection) techniques and methods (Canny, 1986; Lindeberg, 1998; Goodchild and Fueten, 1998; Zhou et al., 2004; Li et al., 2007) in the form of diverse algorithms that are executed either in a spatial domain (physical dimensions of the image) or a frequency domain (frequency that occurs in the spatial domain image) (Gao, 2009). Unlike the frequency domain filters, such as Fourier transformation, that operate on a single band amplitude of low and high frequency curves (i.e., gray scale image), the spatial domain filters are specifically intended to efficiently detect multifrequency edges based on an altered value of an input pixel and values of its neighborhood pixels. This paper focuses on the spatial domain and presents a GISbased approach of edge detection (grain boundaries) in digital images of thin sections from deformed rocks using a cellular automata (CA) approach (von Neumann, 1966). CA approaches have been progressively used to model a great variety of dynamic systems in different application domains such as area zoning and land use change (He et al., 2005; Kocabas and Dragicevic, 2007), urban dynamics studies (Clarke and Gaydos, 1998; Batty et al., 1999; Xia et al., 2010), forest insect infestations and wildfire disturbances (Bone et al., 2006; Alexandridis et al., 2008), artificial intelligence (Fawcett, 2008), and various geological applications (Miyamoto and Sasaki, 1997; D’Ambrosio et al., 2006; Lalonde et al., 2010). However, very few image processing methods combine edge detection and CA based methods (Slatnia and ¨ and Gunay, ¨ Kazar, 2008; Bas- turk 2009). CA are mathematical idealizations of physical systems in which space and time are discrete. CA use a uniform and regular lattice characterized by local neighborhood interactions while physical quantities take on a finite set of possible values such as
137
‘‘on’’ or ‘‘off’’ for updating in discrete time steps. During an interaction, the local state of cell values is affected by its nearest neighborhood cell values from the previous time step, according to a function known as the local transitional rule. The collective effects of local interactions govern the global behavior of dynamic, complex patterns, which can model a range of nonlinear and stochastic spatial processes (Wolfram, 1984; Batty et al., 1999). In this study we perform edge detection in images of deformed rocks using a CA model. The intention of this study is to demonstrate the potential of this methodology, which is based on an automated edge detection time-step process. The proposed approach is illustrated using digital images of thin sections from sandstone quartz arenite and quartz wacke samples and its accuracy is quantitatively evaluated.
2. Modeling approach 2.1. Cellular automata CA are similar to partial differential equations, but are discrete systems with an ability to describe continuous dynamic systems. CA consider space, time, and automaton properties that have a finite and countable number of states. The CA space is an array of cells, which evolves in discrete time steps. At each time step, all cells simultaneously update their states conditioned on their current state and those of their immediate neighbors. For instance, the properties of the CA are based on a regular twodimensional lattice where each cell of the lattice has a discrete state while the dynamic behavior is described by neighboring rules that govern the state of the cells through time. Examples of two-dimensional CA’s include von Neumann’s neighborhood, Moore’s neighborhood, and extended Moore’s neighborhood (Fig. 1). The two-dimensional von Neumann’s neighborhood is plusshaped with four cells surrounding the central cell (x0, y0) and is defined as NvN ðx0 ,y0 Þ ¼ fðx,yÞ : 9xx0 9 þ9yy0 9 r rg,
ð1Þ
where r is the neighborhood radius, which is used to control the size of the neighborhood. For example, by setting the radius to one, the neighborhood comprises cells in the immediate proximity of the central cell and includes all the cells that have a radius of less than and equal to 1 where the straight-line distance is measured from the center of the cell. A radius equal to one comprises four immediate
Fig. 1. Examples of CA neighborhoods. The dark cell is the central cell (x0, y0) while the gray cells are the neighboring cells. The state of neighboring cells predicts the next state of the central cell (x0, y0) according to defined transition rules. In (a) r¼ 1 for the von Neumann neighborhood, (b) r¼ 1 for the Moore’s neighborhood, and (c) r ¼2 for the extended Moore’s neighborhood.
Author's personal copy 138
P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
neighbors plus the central cell. Most common values for the radii vary from zero to three. The two-dimensional Moore’s neighborhood is square-shaped with eight cells surrounding the central cell (x0, y0) and is defined by
represent the colors associated with the image pixels while the local transition function could be stated as ( 0, if 9ssi 9 o e, ð8Þi ¼ 1,n f ðs1 ,s2 ,. . .,sn Þ ¼ ð5Þ s, otherwise
NM ðx0 ,y0 Þ ¼ fðx,yÞ : 9xx0 9 rr,9yy0 9 r rg:
where e is a threshold value for discerning regions of contrasting brightnesses. Here the transition rule is based on comparison criteria of a central cell state with the state of its neighboring cells. If the difference is less than a given threshold value, e, then the state of the central cell at the future state of the test cell at time (tþ 1) will be zero, otherwise it remains unchanged. The proposed transition rule used in this study is discussed below.
ð2Þ
The Moore’s neighborhood is an enlargement of the von Neumann’s neighborhood. The Moore’s neighborhood with a radius equal to 1 uses the eight neighboring cells, including the diagonals, while the extended Moore’s neighborhood with a radius greater than one extends over the distance beyond the eight neighboring cells from the central cell (x0, y0). The definition of neighboring or transition rules is the most important phase in the CA model design because the rules identify the state of the tested cell. Thus, cellular automaton is represented by four elements: A ¼ ðX,S,N,f Þ:
ð3Þ
X is the m dimensional space of each cell x¼(x1, x2, x3, .., xm), S is a nonempty finite set where cells can take only one state at any time from a set of states s A S, N is the neighborhood (i.e., Moore’s neighborhood) where the state of any cell depends on the states of other cells in the neighborhood, and f is the state transition function rule. The future state of a cell at time (tþ1) from starting time (t) is driven by the transition rule as a function of (a) current state of the test cell, (b) current states of its neighborhood cells, and (c) transition rules that drive the process over time. The state transition rule f is represented by f EF St þ 1 ðcÞ ¼ f ðSt ðcÞ,St ðncÞ,transition rulesÞ, tþ1
ð4Þ t
(c) is the future state of the test cell at time (t þ1), S (c) where S is the current state of the test cell at time (t), and St(nc) is the state of the neighborhood. Thus, the transition rule takes the previous state of the test cell and the status of the neighborhood as an input and returns the future state of the test cell at time (t þ1). 2.2. Edge detection Edge detection is a fundamental tool in image processing in the areas of feature detection and feature extraction (Jensen, 2007; Lillesand et al., 2007). In the image, an edge can be regarded as the boundary between two dissimilar regions that may result from changes in intensity, color, or texture. The overall aim of the edge detection is to identify sharp brightness changes or discontinuities in the brightness level. For example, if a linear series of pixels recorded light intensities of 0, 0, 0, 1, 1, 1 an edge or discontinuity would be expected between pixels with recorded intensities of 0 and 1. However, edges are scale dependent and the significance of a particular physical change in an image depends on the nature of the image. Edge detection is usually a subjective task, which often requires an edge enhancement that reinforces the visual boundaries between regions of contrasting brightnesses (Jensen, 2007; Lillesand et al., 2007). The enhancement algorithms consist of a window that is systematically moved through the image comparing the central pixel value to the surrounding pixels to enhance observed thresholds. Detecting edges require the calculation of first- and second-order derivatives often through application of convolution filters. Convolution is central in detecting edges or object boundaries and common manipulation techniques include low-pass filters, which deemphasize local features, high-pass filters used for edges that emphasize local features, or Fourier analysis used for noise removal. However, in a model for edge detection of a digital image based on two-dimensional CA from Eq. (3), the states could
3. Grain boundary detection using cellular automata The grain boundary detection process was applied to two sandstones that have experienced a moderate amount of low temperature deformation: a quartz arenite (Lower Silurian Tuscarora Sandstone) and quartz wacke (Middle–Upper Ordovician Martinsburg Formation). Arenite was chosen because its monomineralic nature presents a challenge to grain boundary detection. The wacke has a polymineralic matrix consisting primarily of illite and quartz, which have a distinctly different grain size and texture than the quartz and rock fragment framework grains. The wacke also has undergone more deformation than the arenite (Onasch, 1983). The thin sections were photographed with a digital camera in cross-polarized light with the gypsum plate (1/4l) inserted. The gypsum plate was used to increase the difference in hue between grains. Viewed without the gypsum plate, the interference colors of the grains are various shades of gray, often with only subtle differences between adjacent grains. To improve the ability to differentiate between grains with similar optical orientations, and hence similar interference color, the sample was photographed at three different angles of stage rotation. This was accomplished by rotating the analyzer and polarizer together relative to a fixed stage. This avoids the problems associated with rotating the stage relative to fixed polarizer and analyzer, which requires that the images then be rotated to coincide with one another. Once the three images were collected, each was separated into red, green, and blue bands and the grain boundary extraction process was applied separately to each. The first step in the extraction process was to apply a 7 7pixel median filter to remove noise created by features present in many grains, such as inclusions, cracks, twins, or subgrains. After the removal of noise, the CA with a Moore’s neighborhood was implemented. The pseudo code for the edge detection algorithm is shown in Fig. 2. For each band, the gray-scale color of each cell is compared to that of its neighbors. Fig. 2 shows that edge detection implements transition rules, which are used when the difference between the cell color and each of the neighbors’ color is greater than a certain incremental threshold value, which suggests that the cell may be a part of an edge. The types of transition rules used in this study are so-called ‘‘totalistic’’ rules (Wolfram, 1983; Adachi et al., 2004), where the future state depends on the sum of the states of the neighborhood cells. For instance, cells with an assigned state of 0 are dead cells while cells with an assigned state of 1 are alive cells. CA transition rules here are implemented through the Moore’s neighborhood where the total number of neighborhood ranges is between 0 and 9 (i.e., maximum alive neighbor cells can be 9 and minimum 0). The cells stay alive as long as there are 4–5 living neighbors, whereas those with fewer than 4 die out from isolation and more than 5 from overpopulation. However, in our simulations we implemented several scenarios consisting of different rules and
Author's personal copy P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
139
Fig. 2. Edge detection algorithm used with Moore’s neighborhood.
Fig. 3. Grains identified using CA model where the threshold is: (a) e ¼ 3, (b) e ¼5, (c) e ¼ 12, and (d) e ¼ 17. The threshold represents the similarity tolerance for the image gray levels (the horizontal axis¼ 3 mm).
threshold values to explore different solution alternatives. Each scenario yielded a best solution of the grain boundaries quantified by comparison of manually digitized and modeled grains using r2 assessment criteria. Following extraction of the boundaries in the R, G, B bands for the images from each of the three analyzer–polarizer rotation positions, the nine images were recombined through the overlay function in the GIS to obtain the final boundaries. The thinning iterative algorithm, which preserves the topology and shape of the original objects, was used to reduce the number of cells, which represented raster linear features. These were then converted to vector features to yield the grain boundary polygons, which were then used for all subsequent analyses. It should be noted that although we used a total of nine images to yield the final grain boundaries, we found that even a single image yielded the majority of the grain boundaries. For our purposes, the last step involved the automated construction of the grain boundary database using a zonal statistics function in which parameters such as area, orientation, and centroids are calculated for each grain. The measures of zonal geometry are particularly useful for exploration and understanding of spatial patterns, such as preferred orientation, that emerge in deformed rocks. The edge detection and CA application was written and implemented using an ARC/INFO GRID module, while
the exploration of spatial patterns used a GeoDA (Anselin et al., 2006). The results from this approach are discussed next.
4. Results Fig. 3 shows the results obtained from the quartz arenite by this model for different threshold values. The red lines are the grain boundaries created by the edge detection process that define the grain polygons. Fig. 3 shows the effect of different threshold values on the size and shape of the grain polygons. For instance, larger polygons are generated when the threshold value is large and smaller polygons are generated when the threshold value is small. The quantitative evaluation of the CA modeled predictions involved a comparison of grain polygons with a validation dataset constructed by manual digitization. The areas of the modelderived and manually digitized polygons were compared at randomly selected locations within the images (Fig. 4). The comparison in Fig. 4 is based on the grain polygons generated from the quartz arenite sample (Fig. 3). If the areas of the modelderived and manually digitized polygons were the same at each location, the best-fit line would have a slope of one and an r2 of 1.0. Deviation from this indicates that the model is less than
Author's personal copy
Model-derived grain polygons
P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
Model-derived grain polygons
140
Polygon area (pixels in hundreds) R-squared = 0.7244 N = 180 ε = 12
100 60 20 0
140
20 40 60 20 80 100 120 Manually-digitized grain polygons
Model-derived grain polygons
Model-derived grain polygons
R-squared = 0.6413 N = 175
100 60 20
N = 145
60
40 60 80 100 120 20 Manually-digitized grain polygons 140
R-squared = 0.5871
100
R-squared = 0.5065
150
N = 138
100 50 0 40 60 80 100 120 20 Manually-digitized grain polygons
40 60 80 100 120 20 Manually-digitized grain polygons
Fig. 4. Relationships between the areas of model-derived and manually digitized grain polygons from four different simulations obtained by applying different transition functions.
ε = 20
1
ε = 19 ε = 18 ε = 17 ε = 16
0.8
ε = 15 ε = 14 ε = 13
0.61
ε = 12 ε = 11 ε = 10
0.42
ε=9 ε=8 ε=7 ε=6
0.23
ε=5 ε=4 ε=3
0.04 0
ε=2 ε=1 Validation
-0.16
Fig. 5. Correlation plot for manually digitized and model-derived grain boundaries. The plot shows that manually digitized grain boundaries have the highest correlation with thresholds e ¼ 12 and e ¼ 13 from the CA simulation.
perfect in identifying the grain boundaries. For the four scenarios tested that use different transition functions, the best relationship has an r2 of 0.724.
To evaluate the effect of changing the threshold value, a comparison was made between the areas of the model-derived and manually digitized polygons using randomly sampled points in the image.
Author's personal copy P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
Fig. 5 illustrates the resulting patterns from one simulation scenario (function for generating the graph was adapted from Zheng et al., 2006), where strong negative correlation is apparent for grain polygons obtained using lower thresholds and strong positive correlation for grain polygons obtained using higher thresholds. The colors in the cells represent a measure of correlation strengths where darker colors signify strong correlations and lighter colors signify weak correlations. The correlation within midrange thresholds is slightly weaker than for grain polygons obtained by higher thresholds, but midrange thresholds correlate strongly with the validation dataset. For instance, the darkest colors are associated with threshold values e ¼ 12 and e ¼13. Also, it appears that threshold values e ¼12 and e ¼ 13 are correlated. Fig. 3c shows the grain boundaries with threshold values e ¼12, which are used for exploring the spatial patterns discussed below.
141
A recent GIS-based method developed by Li et al. (2007) used the same dataset, but a different methodology to extract grain boundaries from the digital images of the thin section. The approach used a series of filtering techniques for enhancing boundaries and reducing noise, and combined unsupervised and supervised classification to support the detection of the grain boundaries. As a comparison, the final results from the CA approach showed a slightly better overall fit between the relationship of model-derived and manually digitized grain polygons than was obtained by Li et al. (2007). The highest relationships obtained by the CA approach yielded higher r-square (r2 ¼0.724 with N ¼180) than the GIS-based approach of Li et al. (2007), which yielded a slightly lower r-square (r2 ¼0.623 with N ¼214). Because the CA approach is flexible, exploring different rules may yield even better results through an optimization process where thresholds and transition rules were systematically investigated. Fig. 6 shows the results obtained from the quartz wacke using an optimum threshold value of e ¼7 determined by trial and error. Framework grains are clearly identified as are many small grains in the matrix. To investigate the effect of using fewer images, runs with all three images (3 images 3 bands/image) and one image (1 image 3 bands) were compared (Fig. 6b and c). The threeimage run, which took 234 min, yields many more matrix grains than the one-image run, which took 77 min. It also shows signs of segmenting some framework grains that have undulatory extinction.
5. Discussion
Fig. 6. (a) Original RGB image, (b) one-band solution, and (c) three-band solution (the horizontal axis¼ 0.75 mm).
Extracting grain boundaries from a digital image of a thin section is a necessary first step for many different types of analyses of rocks. The traditional approach using manual digitization significantly limits the number of samples that can be included, hence, the depth or breadth of the study. The various automated or semiautomated methods that have been developed in recent years offer considerable improvement over manual digitization, both in speed and in reproducibility. The CA method described here has several advantages over other computer-based methods. First, it is flexible in that it can be optimized for a particular type of rock by varying the transitional rules and threshold values. Once optimized, the extraction process can be applied to numerous samples with no operator intervention. As shown, a high degree of correlation between model-derived and manually digitized grain boundaries can be achieved. The deformed quartz arenite was expected to present a particularly difficult challenge to grain boundary extraction due to the large number of microfractures and subgrains; yet, the method yielded a high correlation between model-derived and manually derived grain polygons. The challenge presented by the quartz wacke is that it is a very heterogeneous rock, which is polymineralic and has distinctly different grain size populations. It has also experienced a greater amount of deformation. Visual comparison with the starting image shows that the method clearly identifies framework grains as well as many of the small grains in the matrix. The effect of the number of images used is seen in the detail of the grain boundaries detected. While using fewer images takes less run time, it yields less detail. Depending on the nature of the sample and the intended purpose, one image might yield sufficient detail. For example, in the quartz wacke, accurate determination of the framework grains can be achieved with only one image (Fig. 6b), so if only these grains were of interest, one image would suffice. However, if the matrix grains were of interest, more images would be needed to capture the necessary detail that emerges from different spectral responses captured by the sensor.
Author's personal copy 142
P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
6. Conclusions The methodology presented here is a simple method of edge detection based on cellular automata using a digital image of a thin section. The detection procedure has applicability to both monochromatic (i.e., gray level) and color images. Also, the initial work presents a 2-D image methodology, but this can be expanded to prestack and 3-D images. The selection of transition rules in this study was limited to a few exploratory rules, but the methodology offers a great deal of flexibility for formulating and testing different rules and neighborhoods that are practical for edge detection. Choosing optimal rules for a given edge detection task is important; however, this remains a topic for future work. After the grains are identified, the usefulness of the method is in its ability to construct an automated GIS database for features such as area, perimeter, thickness, and centroid geometry. For geometry type centroid, additional items that describe the ellipsoidal approximation such as major axis, minor axis, and orientation can also be computed. The spatial distribution of features is used to explore patterns that range from completely clustered to completely dispersed, which are useful for understanding processes such as deformation. The results obtained from the sample examined are promising when this method is compared against a traditional manual digitization approach and a recent GIS-based method developed for this purpose by Li et al. (2007). However, further comparison with other edge detection methods based on cellular automata is necessary. Possible future research direction could be extended to capture processes that have spatial and temporal dimensions. For instance, the spatial patterns (spatial dependence and spatial heterogeneity) of grains might evolve in a dynamic context, which is important to integrate the knowledge of form (space) and process (time). Nearby grains might exhibit greater degrees of similarity (or dissimilarity) over time because of some microlevel mechanisms in the rock, depending on the selected indicators such as orientation and shape. At the same time, spatial heterogeneity of the main processes with respect to location in the rock might also be unstable due to the expansion (or shrinking) of individual grains. As such, textural analysis of thin sections of rocks can be conducted in a space–time context. Acknowledgements The authors would like to thank the reviewers, especially Dr. Frank Fueten, who provided helpful suggestions and excellent additions to the article. References Adachi, S., Peper, F., Lee, J., 2004. Universality of hexagonal asynchronous totalistic cellular automata. Lecture Notes in Computer Science 3305, 91–100. Alexandridis, A., Vakalis, D., Siettos, C.I., Bafas, G.V., 2008. A cellular automata model for forest fire spread prediction: the case of the wildfire that swept through Spetses Island in 1990. Applied Mathematics and Computation 204, 191–201. Anselin, L., Syabri, I., Kho, Y., 2006. GeoDa: An introduction to spatial data analysis. Geographical Analysis 38, 5–22. Barraud, J., 2006. The use of watershed segmentation and GIS software for textural analysis of thin sections. Journal of Volcanology and Geothermal Research 154, 17–33. ¨ ¨ Bas-turk, A., Gunay, E., 2009. Efficient edge detection in digital images using a cellular neural network optimized by differential evolution algorithm. Expert Systems with Applications 36, 2645–2650. Batty, M., Xie, Y., Sun, Z., 1999. Modeling urban dynamics through GIS-based cellular automata. Computer, Environment and Urban Systems 23, 1–29. Bone, C., Dragicevic, S., Roberts, A., 2006. A fuzzy-constrained cellular automata model of forest insect infestations. Ecological Modelling192, 107–125. Boorman, S., Boudreau, A., Kruger, F.J., 2004. The lower zone-critical zone transition of the Bushveld complex: A quantitative textural study. Journal of Petrology 45, 1209–1235.
Canny, J., 1986. A Computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence 8, 679–698. Choudhury, K.R., Meere, P.A., Mulchrone, K.F., 2006. Automated grain boundary detection by CASRG. Journal of Structural Geology 28, 363–375. Clarke, K.C., Gaydos, L.J., 1998. Loose coupling a cellular automaton model and GIS: long-term growth prediction for San Francisco and Washington/Baltimore. International Journal of Geographical Information Science 12, 699–714. D’Ambrosio, D., Spataro, W., Iovine, G., 2006. Parallel genetic algorithms for optimising cellular automata models of natural complex phenomena: an application to debris flows. Computers & Geosciences 32, 861–875. Fabbri, A.G., 1984. Image Processing of Geological Data. Van Nostrand Reinhold, New York, 224 pp. Fawcett, T., 2008. Data mining with cellular automata. ACM SIGKDD Explorations Newsletter 10, 32–39. Fueten, F., Goodchild, J.S., 2001. Quartz c-axis orientation determination using the rotating polarizer microscope. Journal of Structural Geology 23, 895–902. Fueten, F., Mason, J., 2007. An artificial neural net assisted approach to editing edges in petrographic images collected with the rotating polarizer stage. Computers and Geosciences 33, 1176–1188. Gao, J., 2009. Digital Analysis of Remotely Sensed Imagery. McGraw-Hill, New York, 645 pp. Goodchild, J.S., Fueten, F., 1998. Edge detection in petrographic images using the rotating polarizer stage. Computers and Geosciences 24, 745–751. He, C., Zhang, Q., Li, Y., Li, X., Shi, P., 2005. Zoning grassland protection area using remote sensing and cellular automata modeling—a case study in Xilingol steppe grassland in northern. Journal of Arid Environments 63, 814–826. Heilbronner, R., 2000. Automatic grain boundary detection and grain size analysis using polarization micrographs or orientation images. Journal of Structural Geology 22, 969–998. Higgins, M.D., 2000. Measurement of crystal size distributions. American Mineralogist 85, 1105–1116. Jensen, J.R., 2007. Remote Sensing of the Environment: An Earth Resource Perspective, second ed. Prentice Hall, New Jersey. Jerram, D.A., Cheadle, M.J., Philpotts, A.R., 2003. Quantifying the building blocks of igneous rocks: are clustered crystal frameworks the foundation? Journal of Petrology 44, 2033–2051. Kocabas, V., Dragicevic, S., 2007. Enhancing a GIS cellular automata model of land use change: Bayesian networks, influence diagrams and causality. Transactions in GIS 11, 681–702. Lalonde, M., Tremblay, G., Je´brak, M., 2010. A Cellular Automata Breccia Simulator (CABS) and its application to rounding in hydrothermal breccias. Computers and Geosciences 36, 827–838. Li, Y., Onasch, C.M., Guo, Y., 2007. GIS-based detection of grain boundaries. Journal of Structural Geology 30, 431–443. Lillesand, T.M., Kiefer, R.W., Chipman, J.W., 2007. Remote Sensing and Image Interpretation, sixth ed. John Wiley and Sons, New York. Lindeberg, T., 1998. Edge detection and ridge detection with automatic scale selection. International Journal of Computer Vision 30, 117–154. Luumbreras, F., Serrat, J., 1996. Segmentation of petrographical images of marbles. Computers and Geosciences 22, 547–558. Miyamoto, H., Sasaki, S., 1997. Simulating lava flows by an improved cellular automata method. Computers and Geosciences 23, 283–292. Obara, B., 2007a. An image processing algorithm for the reversed transformation of rotated microscope images. Computers and Geosciences 33, 853–859. Obara, B., 2007b. A new algorithm using image colour system transformation for rock grain segmentation. Mineralogy and Petrology 91, 271–285. Onasch, C.M., 1983. Cleavage and microstructure in sandstones of the Martinsburg Formation, Maryland. American Journal of Science 283, 936–966. Perring, C.S., Barnes, S.J., Verrall, M., Hill, R.E.T., 2004. Using automated digital image analysis to provide quantitative petrographic data on olivine-phyric basalts. Computers and Geosciences 30, 183–195. Ross, B.J., Fueten, F., Yashkir, D.Y., 2001. Automatic mineral identification using genetic programming. Machine Vision and Applications 13, 61–69. Slatnia, S., Kazar, O., 2008. Generalised evolutionary cellular automata basedapproach for edge detection. International Review on Computers and Software 3, 424–428. van den Berg, E.H., Meesters, A.G.C.A., Kenter, J.A.M., Schlager, W., 2002. Automated separation of touching grains in digital images of thin sections. Computers and Geosciences 28, 179–190. von Neumann, J., 1966. Theory of Self-Reproducing Automata. University of Illinois Press, Champaign, IL. Wolfram, S., 1983. Statistical mechanics of cellular automata. Reviews of Modern Physics 55, 601–644. Wolfram, S., 1984. Universality and complexity in cellular automata. Physica D 10, 1–35. Xia, L., Xiaohu, Z., Yeh, A., Xiaoping, L., 2010. Parallel cellular automata for largescale urban simulation using load-balancing techniques. International Journal of Geographical Information Science 24, 803–820. Zheng, T., Salganik, M., Gelman, A., 2006. How many people do you know in prison?: using over dispersion in count data to estimate social structure in networks. Journal of the American Statistical Association 101, 409–423. Zhou, Y., Starkey, J., Mansinha, L., 2004. Segmentation of petrographic images by integrating edge detection and region growing. Computers and Geosciences 30, 817–831.
Author's personal copy Computers & Geosciences 42 (2012) 136–142
Contents lists available at SciVerse ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Detecting grain boundaries in deformed rocks using a cellular automata approach Pece V. Gorsevski a,n, Charles M. Onasch a, John R. Farver b, Xinyue Ye c a
School of Earth, Environment and Society, Bowling Green State University, Bowling Green, OH 43403, USA Department of Geology, Bowling Green State University, Bowling Green, OH 43403, USA c Center for Regional Development and School of Earth, Environment and Society, Bowling Green State University, Bowling Green, OH 43403, USA b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 22 June 2011 Received in revised form 31 August 2011 Accepted 2 September 2011 Available online 4 October 2011
Cellular automata (CA) are widely used in geospatial dynamic modeling and image processing. Here, we explore the application of two-dimensional cellular automata to the problem of grain boundary detection and extraction in digital images of thin sections from deformed rocks. The automated extraction of boundaries, which contain rich sources of information such as shape, orientation, and spatial distribution of grains, involves a CA Moore’s neighborhood-based rules approach. The Moore’s neighborhood is a 3 3 matrix that is used for changing states by comparing differences between a central pixel and its neighbors. In this dynamic approach, the future state of a pixel depends upon its current state and that of its neighbors. The rules that are defined determine the future state of each cell (i.e., on or off) while the number of iterations to simulate boundaries detection are specified by the user. Each iteration outputs different detection scenarios of grain boundaries that can be evaluated and assessed for accuracy. For a deformed quartz arenite, an r2 of 0.724 was obtained by comparing manually digitized grains to model derived grains. The value of this proposed method is compared against a traditional manual digitization approach and a recent GIS-based method developed for this purpose by Li et al. (2007). & 2011 Elsevier Ltd. All rights reserved.
Keywords: Cellular automata Grain boundary Thin section Edge detection GIS
1. Introduction The examination of thin sections with a polarizing light microscope is a standard way of characterizing rocks. Analyzing the mineralogy, size, shape, and spatial distribution of grains within a rock provides clues about processes such as the crystallization, metamorphism, deformation, and weathering (Fueten and Goodchild, 2001; Choudhury et al., 2006; Li et al., 2007). Although textural properties of rock microstructures are often sufficient for a broad geologic description such as the nature and origin of a rock, understanding of processes that control textural evolution requires objective and accurate measurements of the size, shape, orientation, and position of the mineral grains in a rock followed by an appropriate quantitative analysis. The recognition of grain boundaries is regarded as a vital first step in textural analysis of thin sections of rocks. Typically, grain boundaries are extracted either by manual digitization (Fabbri, 1984; Higgins, 2000; Jerram et al., 2003; Boorman et al., 2004) or by the use of automated extraction using image processing or n
Corresponding author. Tel.: þ1 419 372 7201. E-mail addresses: [email protected] (P.V. Gorsevski), [email protected] (C.M. Onasch), [email protected] (J.R. Farver), [email protected] (X. Ye). 0098-3004/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2011.09.008
GIS-based techniques (Goodchild and Fueten, 1998; van den Berg et al., 2002; Perring et al., 2004; Barraud, 2006; Choudhury et al., 2006; Li et al., 2007; Obara, 2007a,b). Some of the major shortcomings of the manual digitization arise from the labor intensive nature of this technique and lack of reproducibility between different analysts. On the other hand, the automated extraction of grain boundaries is aimed at replacing current, labor intensive manual procedures with consistent, replicable, and robust procedures allowing grains to be identified objectively using computerassisted analysis of images. Some image-specific, computer-assisted techniques that use automated extraction of grain boundaries include microscopic image segmentation using watershed methods (Luumbreras and Serrat, 1996; Zhou et al., 2004; Barraud, 2006), constrained automated seeded region growing (CASRG) algorithm for automated grain edge detection (Choudhury et al., 2006), application of neural networks (Fueten and Mason, 2007), genetic programming (Ross et al., 2001), gradient filtering (Heilbronner, 2000), image classification techniques (Li et al., 2007), and edge detection using sequential image processing routines based on Gaussian smoothing, gradient calculation, thresholding, and image skeletonization (Goodchild and Fueten, 1998). Some of these techniques take advantage of a fully integrated geographic information system (GIS) environment (e.g., Li et al., 2007) that
Author's personal copy P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
makes it relatively simple to implement a range of analytical capabilities such as data management, digital image processing, statistical analysis, spatial modeling, and geoprocessing. The GIS framework allows different algorithms and methodologies to be applied depending on the question of interest and problem that needs to be solved. In particular, GIS integration of microscopederived datasets provides several advantages including access to the functionality of the GIS for display and file management, reduced development time, enhanced usability, and built-in dynamic modeling capabilities. This kind of toolset has great potential for rapid identification and quantification of mineral microstructures that can be compiled in a spatial grain boundary database and quantitatively evaluated for understanding crucial processes that control grain evolution. Edge detection is a fundamental part of image processing, particularly in the areas of detection and extraction of features such as grain boundaries. There are many edge detection (boundary detection) techniques and methods (Canny, 1986; Lindeberg, 1998; Goodchild and Fueten, 1998; Zhou et al., 2004; Li et al., 2007) in the form of diverse algorithms that are executed either in a spatial domain (physical dimensions of the image) or a frequency domain (frequency that occurs in the spatial domain image) (Gao, 2009). Unlike the frequency domain filters, such as Fourier transformation, that operate on a single band amplitude of low and high frequency curves (i.e., gray scale image), the spatial domain filters are specifically intended to efficiently detect multifrequency edges based on an altered value of an input pixel and values of its neighborhood pixels. This paper focuses on the spatial domain and presents a GISbased approach of edge detection (grain boundaries) in digital images of thin sections from deformed rocks using a cellular automata (CA) approach (von Neumann, 1966). CA approaches have been progressively used to model a great variety of dynamic systems in different application domains such as area zoning and land use change (He et al., 2005; Kocabas and Dragicevic, 2007), urban dynamics studies (Clarke and Gaydos, 1998; Batty et al., 1999; Xia et al., 2010), forest insect infestations and wildfire disturbances (Bone et al., 2006; Alexandridis et al., 2008), artificial intelligence (Fawcett, 2008), and various geological applications (Miyamoto and Sasaki, 1997; D’Ambrosio et al., 2006; Lalonde et al., 2010). However, very few image processing methods combine edge detection and CA based methods (Slatnia and ¨ and Gunay, ¨ Kazar, 2008; Bas- turk 2009). CA are mathematical idealizations of physical systems in which space and time are discrete. CA use a uniform and regular lattice characterized by local neighborhood interactions while physical quantities take on a finite set of possible values such as
137
‘‘on’’ or ‘‘off’’ for updating in discrete time steps. During an interaction, the local state of cell values is affected by its nearest neighborhood cell values from the previous time step, according to a function known as the local transitional rule. The collective effects of local interactions govern the global behavior of dynamic, complex patterns, which can model a range of nonlinear and stochastic spatial processes (Wolfram, 1984; Batty et al., 1999). In this study we perform edge detection in images of deformed rocks using a CA model. The intention of this study is to demonstrate the potential of this methodology, which is based on an automated edge detection time-step process. The proposed approach is illustrated using digital images of thin sections from sandstone quartz arenite and quartz wacke samples and its accuracy is quantitatively evaluated.
2. Modeling approach 2.1. Cellular automata CA are similar to partial differential equations, but are discrete systems with an ability to describe continuous dynamic systems. CA consider space, time, and automaton properties that have a finite and countable number of states. The CA space is an array of cells, which evolves in discrete time steps. At each time step, all cells simultaneously update their states conditioned on their current state and those of their immediate neighbors. For instance, the properties of the CA are based on a regular twodimensional lattice where each cell of the lattice has a discrete state while the dynamic behavior is described by neighboring rules that govern the state of the cells through time. Examples of two-dimensional CA’s include von Neumann’s neighborhood, Moore’s neighborhood, and extended Moore’s neighborhood (Fig. 1). The two-dimensional von Neumann’s neighborhood is plusshaped with four cells surrounding the central cell (x0, y0) and is defined as NvN ðx0 ,y0 Þ ¼ fðx,yÞ : 9xx0 9 þ9yy0 9 r rg,
ð1Þ
where r is the neighborhood radius, which is used to control the size of the neighborhood. For example, by setting the radius to one, the neighborhood comprises cells in the immediate proximity of the central cell and includes all the cells that have a radius of less than and equal to 1 where the straight-line distance is measured from the center of the cell. A radius equal to one comprises four immediate
Fig. 1. Examples of CA neighborhoods. The dark cell is the central cell (x0, y0) while the gray cells are the neighboring cells. The state of neighboring cells predicts the next state of the central cell (x0, y0) according to defined transition rules. In (a) r¼ 1 for the von Neumann neighborhood, (b) r¼ 1 for the Moore’s neighborhood, and (c) r ¼2 for the extended Moore’s neighborhood.
Author's personal copy 138
P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
neighbors plus the central cell. Most common values for the radii vary from zero to three. The two-dimensional Moore’s neighborhood is square-shaped with eight cells surrounding the central cell (x0, y0) and is defined by
represent the colors associated with the image pixels while the local transition function could be stated as ( 0, if 9ssi 9 o e, ð8Þi ¼ 1,n f ðs1 ,s2 ,. . .,sn Þ ¼ ð5Þ s, otherwise
NM ðx0 ,y0 Þ ¼ fðx,yÞ : 9xx0 9 rr,9yy0 9 r rg:
where e is a threshold value for discerning regions of contrasting brightnesses. Here the transition rule is based on comparison criteria of a central cell state with the state of its neighboring cells. If the difference is less than a given threshold value, e, then the state of the central cell at the future state of the test cell at time (tþ 1) will be zero, otherwise it remains unchanged. The proposed transition rule used in this study is discussed below.
ð2Þ
The Moore’s neighborhood is an enlargement of the von Neumann’s neighborhood. The Moore’s neighborhood with a radius equal to 1 uses the eight neighboring cells, including the diagonals, while the extended Moore’s neighborhood with a radius greater than one extends over the distance beyond the eight neighboring cells from the central cell (x0, y0). The definition of neighboring or transition rules is the most important phase in the CA model design because the rules identify the state of the tested cell. Thus, cellular automaton is represented by four elements: A ¼ ðX,S,N,f Þ:
ð3Þ
X is the m dimensional space of each cell x¼(x1, x2, x3, .., xm), S is a nonempty finite set where cells can take only one state at any time from a set of states s A S, N is the neighborhood (i.e., Moore’s neighborhood) where the state of any cell depends on the states of other cells in the neighborhood, and f is the state transition function rule. The future state of a cell at time (tþ1) from starting time (t) is driven by the transition rule as a function of (a) current state of the test cell, (b) current states of its neighborhood cells, and (c) transition rules that drive the process over time. The state transition rule f is represented by f EF St þ 1 ðcÞ ¼ f ðSt ðcÞ,St ðncÞ,transition rulesÞ, tþ1
ð4Þ t
(c) is the future state of the test cell at time (t þ1), S (c) where S is the current state of the test cell at time (t), and St(nc) is the state of the neighborhood. Thus, the transition rule takes the previous state of the test cell and the status of the neighborhood as an input and returns the future state of the test cell at time (t þ1). 2.2. Edge detection Edge detection is a fundamental tool in image processing in the areas of feature detection and feature extraction (Jensen, 2007; Lillesand et al., 2007). In the image, an edge can be regarded as the boundary between two dissimilar regions that may result from changes in intensity, color, or texture. The overall aim of the edge detection is to identify sharp brightness changes or discontinuities in the brightness level. For example, if a linear series of pixels recorded light intensities of 0, 0, 0, 1, 1, 1 an edge or discontinuity would be expected between pixels with recorded intensities of 0 and 1. However, edges are scale dependent and the significance of a particular physical change in an image depends on the nature of the image. Edge detection is usually a subjective task, which often requires an edge enhancement that reinforces the visual boundaries between regions of contrasting brightnesses (Jensen, 2007; Lillesand et al., 2007). The enhancement algorithms consist of a window that is systematically moved through the image comparing the central pixel value to the surrounding pixels to enhance observed thresholds. Detecting edges require the calculation of first- and second-order derivatives often through application of convolution filters. Convolution is central in detecting edges or object boundaries and common manipulation techniques include low-pass filters, which deemphasize local features, high-pass filters used for edges that emphasize local features, or Fourier analysis used for noise removal. However, in a model for edge detection of a digital image based on two-dimensional CA from Eq. (3), the states could
3. Grain boundary detection using cellular automata The grain boundary detection process was applied to two sandstones that have experienced a moderate amount of low temperature deformation: a quartz arenite (Lower Silurian Tuscarora Sandstone) and quartz wacke (Middle–Upper Ordovician Martinsburg Formation). Arenite was chosen because its monomineralic nature presents a challenge to grain boundary detection. The wacke has a polymineralic matrix consisting primarily of illite and quartz, which have a distinctly different grain size and texture than the quartz and rock fragment framework grains. The wacke also has undergone more deformation than the arenite (Onasch, 1983). The thin sections were photographed with a digital camera in cross-polarized light with the gypsum plate (1/4l) inserted. The gypsum plate was used to increase the difference in hue between grains. Viewed without the gypsum plate, the interference colors of the grains are various shades of gray, often with only subtle differences between adjacent grains. To improve the ability to differentiate between grains with similar optical orientations, and hence similar interference color, the sample was photographed at three different angles of stage rotation. This was accomplished by rotating the analyzer and polarizer together relative to a fixed stage. This avoids the problems associated with rotating the stage relative to fixed polarizer and analyzer, which requires that the images then be rotated to coincide with one another. Once the three images were collected, each was separated into red, green, and blue bands and the grain boundary extraction process was applied separately to each. The first step in the extraction process was to apply a 7 7pixel median filter to remove noise created by features present in many grains, such as inclusions, cracks, twins, or subgrains. After the removal of noise, the CA with a Moore’s neighborhood was implemented. The pseudo code for the edge detection algorithm is shown in Fig. 2. For each band, the gray-scale color of each cell is compared to that of its neighbors. Fig. 2 shows that edge detection implements transition rules, which are used when the difference between the cell color and each of the neighbors’ color is greater than a certain incremental threshold value, which suggests that the cell may be a part of an edge. The types of transition rules used in this study are so-called ‘‘totalistic’’ rules (Wolfram, 1983; Adachi et al., 2004), where the future state depends on the sum of the states of the neighborhood cells. For instance, cells with an assigned state of 0 are dead cells while cells with an assigned state of 1 are alive cells. CA transition rules here are implemented through the Moore’s neighborhood where the total number of neighborhood ranges is between 0 and 9 (i.e., maximum alive neighbor cells can be 9 and minimum 0). The cells stay alive as long as there are 4–5 living neighbors, whereas those with fewer than 4 die out from isolation and more than 5 from overpopulation. However, in our simulations we implemented several scenarios consisting of different rules and
Author's personal copy P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
139
Fig. 2. Edge detection algorithm used with Moore’s neighborhood.
Fig. 3. Grains identified using CA model where the threshold is: (a) e ¼ 3, (b) e ¼5, (c) e ¼ 12, and (d) e ¼ 17. The threshold represents the similarity tolerance for the image gray levels (the horizontal axis¼ 3 mm).
threshold values to explore different solution alternatives. Each scenario yielded a best solution of the grain boundaries quantified by comparison of manually digitized and modeled grains using r2 assessment criteria. Following extraction of the boundaries in the R, G, B bands for the images from each of the three analyzer–polarizer rotation positions, the nine images were recombined through the overlay function in the GIS to obtain the final boundaries. The thinning iterative algorithm, which preserves the topology and shape of the original objects, was used to reduce the number of cells, which represented raster linear features. These were then converted to vector features to yield the grain boundary polygons, which were then used for all subsequent analyses. It should be noted that although we used a total of nine images to yield the final grain boundaries, we found that even a single image yielded the majority of the grain boundaries. For our purposes, the last step involved the automated construction of the grain boundary database using a zonal statistics function in which parameters such as area, orientation, and centroids are calculated for each grain. The measures of zonal geometry are particularly useful for exploration and understanding of spatial patterns, such as preferred orientation, that emerge in deformed rocks. The edge detection and CA application was written and implemented using an ARC/INFO GRID module, while
the exploration of spatial patterns used a GeoDA (Anselin et al., 2006). The results from this approach are discussed next.
4. Results Fig. 3 shows the results obtained from the quartz arenite by this model for different threshold values. The red lines are the grain boundaries created by the edge detection process that define the grain polygons. Fig. 3 shows the effect of different threshold values on the size and shape of the grain polygons. For instance, larger polygons are generated when the threshold value is large and smaller polygons are generated when the threshold value is small. The quantitative evaluation of the CA modeled predictions involved a comparison of grain polygons with a validation dataset constructed by manual digitization. The areas of the modelderived and manually digitized polygons were compared at randomly selected locations within the images (Fig. 4). The comparison in Fig. 4 is based on the grain polygons generated from the quartz arenite sample (Fig. 3). If the areas of the modelderived and manually digitized polygons were the same at each location, the best-fit line would have a slope of one and an r2 of 1.0. Deviation from this indicates that the model is less than
Author's personal copy
Model-derived grain polygons
P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
Model-derived grain polygons
140
Polygon area (pixels in hundreds) R-squared = 0.7244 N = 180 ε = 12
100 60 20 0
140
20 40 60 20 80 100 120 Manually-digitized grain polygons
Model-derived grain polygons
Model-derived grain polygons
R-squared = 0.6413 N = 175
100 60 20
N = 145
60
40 60 80 100 120 20 Manually-digitized grain polygons 140
R-squared = 0.5871
100
R-squared = 0.5065
150
N = 138
100 50 0 40 60 80 100 120 20 Manually-digitized grain polygons
40 60 80 100 120 20 Manually-digitized grain polygons
Fig. 4. Relationships between the areas of model-derived and manually digitized grain polygons from four different simulations obtained by applying different transition functions.
ε = 20
1
ε = 19 ε = 18 ε = 17 ε = 16
0.8
ε = 15 ε = 14 ε = 13
0.61
ε = 12 ε = 11 ε = 10
0.42
ε=9 ε=8 ε=7 ε=6
0.23
ε=5 ε=4 ε=3
0.04 0
ε=2 ε=1 Validation
-0.16
Fig. 5. Correlation plot for manually digitized and model-derived grain boundaries. The plot shows that manually digitized grain boundaries have the highest correlation with thresholds e ¼ 12 and e ¼ 13 from the CA simulation.
perfect in identifying the grain boundaries. For the four scenarios tested that use different transition functions, the best relationship has an r2 of 0.724.
To evaluate the effect of changing the threshold value, a comparison was made between the areas of the model-derived and manually digitized polygons using randomly sampled points in the image.
Author's personal copy P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
Fig. 5 illustrates the resulting patterns from one simulation scenario (function for generating the graph was adapted from Zheng et al., 2006), where strong negative correlation is apparent for grain polygons obtained using lower thresholds and strong positive correlation for grain polygons obtained using higher thresholds. The colors in the cells represent a measure of correlation strengths where darker colors signify strong correlations and lighter colors signify weak correlations. The correlation within midrange thresholds is slightly weaker than for grain polygons obtained by higher thresholds, but midrange thresholds correlate strongly with the validation dataset. For instance, the darkest colors are associated with threshold values e ¼ 12 and e ¼13. Also, it appears that threshold values e ¼12 and e ¼ 13 are correlated. Fig. 3c shows the grain boundaries with threshold values e ¼12, which are used for exploring the spatial patterns discussed below.
141
A recent GIS-based method developed by Li et al. (2007) used the same dataset, but a different methodology to extract grain boundaries from the digital images of the thin section. The approach used a series of filtering techniques for enhancing boundaries and reducing noise, and combined unsupervised and supervised classification to support the detection of the grain boundaries. As a comparison, the final results from the CA approach showed a slightly better overall fit between the relationship of model-derived and manually digitized grain polygons than was obtained by Li et al. (2007). The highest relationships obtained by the CA approach yielded higher r-square (r2 ¼0.724 with N ¼180) than the GIS-based approach of Li et al. (2007), which yielded a slightly lower r-square (r2 ¼0.623 with N ¼214). Because the CA approach is flexible, exploring different rules may yield even better results through an optimization process where thresholds and transition rules were systematically investigated. Fig. 6 shows the results obtained from the quartz wacke using an optimum threshold value of e ¼7 determined by trial and error. Framework grains are clearly identified as are many small grains in the matrix. To investigate the effect of using fewer images, runs with all three images (3 images 3 bands/image) and one image (1 image 3 bands) were compared (Fig. 6b and c). The threeimage run, which took 234 min, yields many more matrix grains than the one-image run, which took 77 min. It also shows signs of segmenting some framework grains that have undulatory extinction.
5. Discussion
Fig. 6. (a) Original RGB image, (b) one-band solution, and (c) three-band solution (the horizontal axis¼ 0.75 mm).
Extracting grain boundaries from a digital image of a thin section is a necessary first step for many different types of analyses of rocks. The traditional approach using manual digitization significantly limits the number of samples that can be included, hence, the depth or breadth of the study. The various automated or semiautomated methods that have been developed in recent years offer considerable improvement over manual digitization, both in speed and in reproducibility. The CA method described here has several advantages over other computer-based methods. First, it is flexible in that it can be optimized for a particular type of rock by varying the transitional rules and threshold values. Once optimized, the extraction process can be applied to numerous samples with no operator intervention. As shown, a high degree of correlation between model-derived and manually digitized grain boundaries can be achieved. The deformed quartz arenite was expected to present a particularly difficult challenge to grain boundary extraction due to the large number of microfractures and subgrains; yet, the method yielded a high correlation between model-derived and manually derived grain polygons. The challenge presented by the quartz wacke is that it is a very heterogeneous rock, which is polymineralic and has distinctly different grain size populations. It has also experienced a greater amount of deformation. Visual comparison with the starting image shows that the method clearly identifies framework grains as well as many of the small grains in the matrix. The effect of the number of images used is seen in the detail of the grain boundaries detected. While using fewer images takes less run time, it yields less detail. Depending on the nature of the sample and the intended purpose, one image might yield sufficient detail. For example, in the quartz wacke, accurate determination of the framework grains can be achieved with only one image (Fig. 6b), so if only these grains were of interest, one image would suffice. However, if the matrix grains were of interest, more images would be needed to capture the necessary detail that emerges from different spectral responses captured by the sensor.
Author's personal copy 142
P.V. Gorsevski et al. / Computers & Geosciences 42 (2012) 136–142
6. Conclusions The methodology presented here is a simple method of edge detection based on cellular automata using a digital image of a thin section. The detection procedure has applicability to both monochromatic (i.e., gray level) and color images. Also, the initial work presents a 2-D image methodology, but this can be expanded to prestack and 3-D images. The selection of transition rules in this study was limited to a few exploratory rules, but the methodology offers a great deal of flexibility for formulating and testing different rules and neighborhoods that are practical for edge detection. Choosing optimal rules for a given edge detection task is important; however, this remains a topic for future work. After the grains are identified, the usefulness of the method is in its ability to construct an automated GIS database for features such as area, perimeter, thickness, and centroid geometry. For geometry type centroid, additional items that describe the ellipsoidal approximation such as major axis, minor axis, and orientation can also be computed. The spatial distribution of features is used to explore patterns that range from completely clustered to completely dispersed, which are useful for understanding processes such as deformation. The results obtained from the sample examined are promising when this method is compared against a traditional manual digitization approach and a recent GIS-based method developed for this purpose by Li et al. (2007). However, further comparison with other edge detection methods based on cellular automata is necessary. Possible future research direction could be extended to capture processes that have spatial and temporal dimensions. For instance, the spatial patterns (spatial dependence and spatial heterogeneity) of grains might evolve in a dynamic context, which is important to integrate the knowledge of form (space) and process (time). Nearby grains might exhibit greater degrees of similarity (or dissimilarity) over time because of some microlevel mechanisms in the rock, depending on the selected indicators such as orientation and shape. At the same time, spatial heterogeneity of the main processes with respect to location in the rock might also be unstable due to the expansion (or shrinking) of individual grains. As such, textural analysis of thin sections of rocks can be conducted in a space–time context. Acknowledgements The authors would like to thank the reviewers, especially Dr. Frank Fueten, who provided helpful suggestions and excellent additions to the article. References Adachi, S., Peper, F., Lee, J., 2004. Universality of hexagonal asynchronous totalistic cellular automata. Lecture Notes in Computer Science 3305, 91–100. Alexandridis, A., Vakalis, D., Siettos, C.I., Bafas, G.V., 2008. A cellular automata model for forest fire spread prediction: the case of the wildfire that swept through Spetses Island in 1990. Applied Mathematics and Computation 204, 191–201. Anselin, L., Syabri, I., Kho, Y., 2006. GeoDa: An introduction to spatial data analysis. Geographical Analysis 38, 5–22. Barraud, J., 2006. The use of watershed segmentation and GIS software for textural analysis of thin sections. Journal of Volcanology and Geothermal Research 154, 17–33. ¨ ¨ Bas-turk, A., Gunay, E., 2009. Efficient edge detection in digital images using a cellular neural network optimized by differential evolution algorithm. Expert Systems with Applications 36, 2645–2650. Batty, M., Xie, Y., Sun, Z., 1999. Modeling urban dynamics through GIS-based cellular automata. Computer, Environment and Urban Systems 23, 1–29. Bone, C., Dragicevic, S., Roberts, A., 2006. A fuzzy-constrained cellular automata model of forest insect infestations. Ecological Modelling192, 107–125. Boorman, S., Boudreau, A., Kruger, F.J., 2004. The lower zone-critical zone transition of the Bushveld complex: A quantitative textural study. Journal of Petrology 45, 1209–1235.
Canny, J., 1986. A Computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence 8, 679–698. Choudhury, K.R., Meere, P.A., Mulchrone, K.F., 2006. Automated grain boundary detection by CASRG. Journal of Structural Geology 28, 363–375. Clarke, K.C., Gaydos, L.J., 1998. Loose coupling a cellular automaton model and GIS: long-term growth prediction for San Francisco and Washington/Baltimore. International Journal of Geographical Information Science 12, 699–714. D’Ambrosio, D., Spataro, W., Iovine, G., 2006. Parallel genetic algorithms for optimising cellular automata models of natural complex phenomena: an application to debris flows. Computers & Geosciences 32, 861–875. Fabbri, A.G., 1984. Image Processing of Geological Data. Van Nostrand Reinhold, New York, 224 pp. Fawcett, T., 2008. Data mining with cellular automata. ACM SIGKDD Explorations Newsletter 10, 32–39. Fueten, F., Goodchild, J.S., 2001. Quartz c-axis orientation determination using the rotating polarizer microscope. Journal of Structural Geology 23, 895–902. Fueten, F., Mason, J., 2007. An artificial neural net assisted approach to editing edges in petrographic images collected with the rotating polarizer stage. Computers and Geosciences 33, 1176–1188. Gao, J., 2009. Digital Analysis of Remotely Sensed Imagery. McGraw-Hill, New York, 645 pp. Goodchild, J.S., Fueten, F., 1998. Edge detection in petrographic images using the rotating polarizer stage. Computers and Geosciences 24, 745–751. He, C., Zhang, Q., Li, Y., Li, X., Shi, P., 2005. Zoning grassland protection area using remote sensing and cellular automata modeling—a case study in Xilingol steppe grassland in northern. Journal of Arid Environments 63, 814–826. Heilbronner, R., 2000. Automatic grain boundary detection and grain size analysis using polarization micrographs or orientation images. Journal of Structural Geology 22, 969–998. Higgins, M.D., 2000. Measurement of crystal size distributions. American Mineralogist 85, 1105–1116. Jensen, J.R., 2007. Remote Sensing of the Environment: An Earth Resource Perspective, second ed. Prentice Hall, New Jersey. Jerram, D.A., Cheadle, M.J., Philpotts, A.R., 2003. Quantifying the building blocks of igneous rocks: are clustered crystal frameworks the foundation? Journal of Petrology 44, 2033–2051. Kocabas, V., Dragicevic, S., 2007. Enhancing a GIS cellular automata model of land use change: Bayesian networks, influence diagrams and causality. Transactions in GIS 11, 681–702. Lalonde, M., Tremblay, G., Je´brak, M., 2010. A Cellular Automata Breccia Simulator (CABS) and its application to rounding in hydrothermal breccias. Computers and Geosciences 36, 827–838. Li, Y., Onasch, C.M., Guo, Y., 2007. GIS-based detection of grain boundaries. Journal of Structural Geology 30, 431–443. Lillesand, T.M., Kiefer, R.W., Chipman, J.W., 2007. Remote Sensing and Image Interpretation, sixth ed. John Wiley and Sons, New York. Lindeberg, T., 1998. Edge detection and ridge detection with automatic scale selection. International Journal of Computer Vision 30, 117–154. Luumbreras, F., Serrat, J., 1996. Segmentation of petrographical images of marbles. Computers and Geosciences 22, 547–558. Miyamoto, H., Sasaki, S., 1997. Simulating lava flows by an improved cellular automata method. Computers and Geosciences 23, 283–292. Obara, B., 2007a. An image processing algorithm for the reversed transformation of rotated microscope images. Computers and Geosciences 33, 853–859. Obara, B., 2007b. A new algorithm using image colour system transformation for rock grain segmentation. Mineralogy and Petrology 91, 271–285. Onasch, C.M., 1983. Cleavage and microstructure in sandstones of the Martinsburg Formation, Maryland. American Journal of Science 283, 936–966. Perring, C.S., Barnes, S.J., Verrall, M., Hill, R.E.T., 2004. Using automated digital image analysis to provide quantitative petrographic data on olivine-phyric basalts. Computers and Geosciences 30, 183–195. Ross, B.J., Fueten, F., Yashkir, D.Y., 2001. Automatic mineral identification using genetic programming. Machine Vision and Applications 13, 61–69. Slatnia, S., Kazar, O., 2008. Generalised evolutionary cellular automata basedapproach for edge detection. International Review on Computers and Software 3, 424–428. van den Berg, E.H., Meesters, A.G.C.A., Kenter, J.A.M., Schlager, W., 2002. Automated separation of touching grains in digital images of thin sections. Computers and Geosciences 28, 179–190. von Neumann, J., 1966. Theory of Self-Reproducing Automata. University of Illinois Press, Champaign, IL. Wolfram, S., 1983. Statistical mechanics of cellular automata. Reviews of Modern Physics 55, 601–644. Wolfram, S., 1984. Universality and complexity in cellular automata. Physica D 10, 1–35. Xia, L., Xiaohu, Z., Yeh, A., Xiaoping, L., 2010. Parallel cellular automata for largescale urban simulation using load-balancing techniques. International Journal of Geographical Information Science 24, 803–820. Zheng, T., Salganik, M., Gelman, A., 2006. How many people do you know in prison?: using over dispersion in count data to estimate social structure in networks. Journal of the American Statistical Association 101, 409–423. Zhou, Y., Starkey, J., Mansinha, L., 2004. Segmentation of petrographic images by integrating edge detection and region growing. Computers and Geosciences 30, 817–831.
