- Email: [email protected]
Automated grain boundary detection using the level set method
ARTICLE IN PRESS Computers & Geosciences 35 (2009) 267–275
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Automated grain boundary detection using the level set method$ Bibo Lu a,b, Min Cui c, Qiang liu a, Yangang Wang a,d, a
College of Mathematics, Jilin University, Changchun 130012, PR China School of Computer Science and Technology, Henan Polytechnic University, Jiaozuo 454000, P.R. China c Basin and Reservoir Research Centre, China University of Petroleum, Beijing 102249, PR China d Supercomputing Center, Computer Network Information Center, Chinese Academy of Sciences, Beijing 100190, P.R. China b
a r t i c l e in fo
abstract
Article history: Received 9 February 2007 Received in revised form 7 May 2008 Accepted 14 May 2008
In this paper, we detect the grain boundary by using a novel automated algorithm which is based on the level set method (the LSM). The LSM starts with only one initial closed curve which can be arbitrary in both shape and location in the microphotograph. The curve splits and merges repeatedly until it reaches the boundaries of the objects. The evolution depends on the structure of the image. This method requires simple and feasible constraints which ensure the exclusion of the objects unsuitable for reliable measurements and strain analysis. The LSM has the following three advantages: (1) with a single image as input, segmentation is automated without human intervention, (2) the grain boundary detected by the LSM is a closed curve, and (3) simple post-processing maintains the shape of the boundary. The segmentation and strain analysis by the LSM are demonstrated and the results are compared with those from hand-drawn method. & 2008 Elsevier Ltd. All rights reserved.
Keywords: Grain boundary detection The level set method Automated image segmentation Strain analysis
1. Introduction Strain analysis plays an important role in the study of structural geology, especially in the investigation of the tectonic history of a region. In theory, high sample intensity is necessary to provide enough raw data for strain analysis. However, in practice, the samples used for strain analysis seem to be insufficient, partly because collecting raw data for strain analysis is laborious and time consuming. Even though many methods have been developed, efficiency and accuracy are yet to be attained. In this paper, we focus on the problem of extracting data from a single digital microphotograph for strain analysis. We detect the grain boundary by using a novel automated algorithm which is based on the level set
$ Code available from server at http://www.iamg.org/CGEditor/index. htm. Corresponding author. Tel./fax: +86 43185167051. E-mail addresses: [email protected] (B. Lu), [email protected] (Y. Wang).
0098-3004/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2008.05.006
method (the LSM). Our analysis of the microphotograph of the mylonites of Daqing Mountain in Inner Mongolia, China, indicates that the LSM enables us to extract the data which are necessary for an accurate strain analysis. Traditional methods of measurement required sustained use of the polarizing microscope with skilled manipulation of the rotating stage and knowledge of the use of various graticules (Ramsay, 1967). An automatic method proposed by Goodchild and Fueten (1998) constructed a composite data set from a number of sequential images obtained from the polarizing microscope with automatic rotation stages. By using some digital images of the grain slice, geologists were able to detect the boundaries of clasts manually. They did this with the help of image processing software, such as Photoshop; and it is known that once the boundaries are obtained, other methods can be employed to make the measurements required for strain analysis. To improve the efficiency of the boundary detection, various semi-automated and automated segmentation methods have been proposed. Most of them are based
ARTICLE IN PRESS 268
B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
on image processing technology, particularly, image segmentation and boundary detection algorithms. These methods can be classified into two main categories: edgedriven methods and region-driven methods. The methods in the first category try to detect the grain boundary directly. Mlynarczuk (1999) proposed an algorithm based on the watershed method for segmenting grains from polarizing microscope images of dolomite and quartzite. Some traditional operators, such as the Robert operator and the Canny operator, were used to extract boundaries defined by those points with large gradient magnitude. Starkey and Samantaray (1993) combined the outcomes of a number of sequential filtered images to produce a final edge image. Lumbreras and Serrat (1996) proposed a two-step strategy to segment an image by performing an over-segmentation of the image and then using a region-merging procedure to get the final boundary. The key point of the region-merging procedure is to estimate two parameters from the information of a sequence of images of the same sample. Using a series of image processing routines, Goodchild and Fueten (1998) extracted the edges from a gradient image obtained by passing a simple gradient operator over images captured during the sampling process. Heilbronner (2000) developed an automatic method for boundary detection and grain size analysis by using polarizing micrographs or orientation images. Bartozzi et al. (2000) employed a detection-filtering algorithm to automatically rebuild a real boundary net, where boundaries are defined as high gradient features on an orientation contrast images. Recently, Obara and Kozˇusˇni´kova´ (2007) proposed a new method of detecting the calcite grains in marble. Based on a set of colored images filtered by a morphological alternating sequence filter, a composite image is constructed and the final grain boundary is determined by watershed segmentation technique. The initial boundary produced by the edge detection methods suffers from some imperfections. Such a boundary, which may be open or discontinuous, does not coincide with the real boundary. Consequently, a post-processing method, such as mathematical morphology, is required in order to obtain a closed boundary. However, the subjectivity of the postprocessing makes these methods unsuitable for strain analysis (Roy Choudhury et al., 2006). The region-driven methods are based mainly on a popular image segmentation technique, i.e., the region growing method, which identifies the region of the grain by absorbing all the points similar to the seed point. The major advantage of the region growing approach is that the boundary of the region is closed and this is preferred by both measurements and strain analysis. A critical step in region growing method is to select the proper seed, which often involves human interventions in order to help the seed grow. Sometimes, some holes occur in the identified clast. The size of the hole depends on the distribution of the noise and the area of the heterogeneity of color information within the clast. While the small holes can be filled by employing simple morphological operations, the large ones cannot, and remain almost the same as before. Post-processing affects the extracted boundaries and features, such as major and minor axes,
orientation, and centroid. A recent method proposed by Roy Choudhury et al. (2006), CASRG, identified a region which is a set of points of the grain. They selected some seeds manually and chose the optimal threshold separately for each grain. This method is very efficient and is accurate for low and high strain samples. However, this method requires that all the seeds be chosen manually, using the computer mouse. But where large or high stress samples are involved, selecting seeds by clicking the mouse is onerous. Zhou et al. (2004) proposed a segmentation method of petrographic images by integrating the edge detection method with that of the region growing. They employed a boundary detection method to get the initial edge whose information is used for the automated seed selection in the region growing method. The methods mentioned above are mainly based on the traditional edge detection and image segmentation technology in image processing. Human interventions are required during the segmentation, and the objectivity of the data is inevitably corrupted by the post-processing. In this paper, we propose a new automated grain boundary detection method: LSM. The LSM requires a little mathematical knowledge of partial differential equation (PDE). We show that the LSM produces impressive results in grain boundary detection as well as in strain analysis. The remaining part of the paper is organized as follows. Section 2 briefly introduces the LSM. In Section 3, we discuss the algorithm for evolving the level set and feasible post-processing. In Section 4, we demonstrate the performance of the LSM. Section 5 is discussion and conclusions. 2. The LSM The LSM aims at identifying grain boundaries and representing them as closed outlines. Unlike CASGR (Roy Choudhury et al., 2006) which manually identifies all the interesting objects one by one, the LSM identifies all the objects that stand out from average gray levels simultaneously. The LSM is a PDE-based segmentation model with the level set formulation. Recently, methods based on PDE have been widely used in image processing and analysis; and impressive results have been obtained by various PDE models. Using level set, Chan and Vese (2001) proposed a successful segmentation model (C–V model) that is based on the classical active contour model (Kass et al., 1988). The essential idea in using the C–V model is to evolve an initial contour governed by a PDE until it stops on the edges of the objects. To introduce the LSM, we first of all review the C–V model and then formulate the model in terms of the level set. The detailed algorithm will be presented in the next section. 2.1. The C–V model Let us begin our discussion with a simple case. Assume that the image u0 consists of one object and background with constant intensity respectively. The value of the
ARTICLE IN PRESS B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
object is u1 and the background is u2 . Denote the boundary of the object by C 0. Then we have u0 u1 inside the object and u0 u2 outside the object. Then let us consider the following fitting term: Z ju0 ðx; yÞ c1 j2 dx dy F 1 ðCÞ þ F 2 ðCÞ ¼ insideðCÞ Z þ ju0 ðx; yÞ c2 j2 dx dy, outsideðCÞ
where C is any other variable curve, and the constants c1 ; c2 are the averages of u0 inside and outside C, respectively. In this simple example, four possible cases are illustrated in Fig. 1. In the first case, the curve C is outside the object. Thus, we have F 1 ðCÞ40 and F 2 ðCÞ ¼ 0 (see Fig. 1(a)). In the second case, the curve C is inside the object. Thus, we have F 1 ðCÞ ¼ 0 but F 2 ðCÞ40 (see Fig. 1(b)). In the third case, since the curve C is both inside and outside the object, we have F 1 ðCÞ40 and F 2 ðCÞ40 (see Fig. 1(c)). In the fourth case, it is only when the curve C coincides with C 0 , the boundary of the object, that the above fitting term can reach the minimizer (see Fig. 1(d)); that is
269
generates two other problems: (1) how to represent the curve used to trace the boundary of the object and (2) how to evolve the curve until it stops on the boundary of the object. The LSM provides a good solution to the first problem. 2.2. Level set formulation of the model The LSM (Osher and Sethian, 1988) represents a curve as the zero level set of a high-dimensional function and allows for major simplifications. The LSM offers a natural representation of the contour of the object. Thus, it is capable of dealing with complex structures with many advantages, especially when the curve undergoes complicated topological changes. The LSM has the following advantages: (1) It is easy to implement numerically. (2) The outline of the object is a closed curve. (3) It can easily deal with the topological changes. such as merging and splitting. (4) It is able to compute geometric quantities easily.
inf C ðF 1 ðCÞ þ F 2 ðCÞÞ 0 F 1 ðC 0 Þ þ F 2 ðC 0 Þ. Therefore, the energy functional Fðc1 ; c2 ; CÞ is defined by Z ju0 ðx; yÞ c1 j2 dx dy Fðc1 ; c2 ; CÞ ¼ l1 insideðCÞ Z ju0 ðx; yÞ c2 j2 dx dy, (1) þ l2 outsideðCÞ
where l1 40 and l2 40 are fixed parameters. Now we convert the segmentation of the image into a minimization problem: inf Fðc1 ; c2 ; CÞ.
c1 ;c2 ;C
(2)
Solving the optimization problem is equivalent to segmenting the image; and the attempt to solve the problem
In the LSM, a closed curve is the zero level set of a function in a high-dimensional space. Let fðx; y; tÞ be a level set function defined on the whole image u0 . The level set function is determined by the position of the points. fðx; y; tÞo0 if the point is inside the curve and fðx; y; tÞ40 if the point is outside the curve. The interface between the region fðx; y; tÞ40 and the region fðx; y; tÞo0 is the curve C, which we call the zero level set because f ¼ 0 on this curve. The zero level set can be considered as the interface between a surface and a plane. The level set evolves in time under the control of the appointed model. Fig. 2 illustrates the change of the level set function and the corresponding curve. In our method, we use the level set to capture the contour of the object. With the level set, the minimization problem can be deduced into the following equation after some mathematical manipulations (see Chan and Vese, 2001 for more details):
qf ¼ j,fj½l1 ðu0 c1 Þ2 þ l2 ðu0 c2 Þ2 , qt
(3)
=qffiffi t is the partial derivative of f with respect to t where qfp , j,fj ¼ ðqf=qxÞ2 þ ðqf=qyÞ2 . Given an initial level set, the level set f evolves according to Eq. (3). The regularization term in the original model is dropped to maintain the primal grain boundary. 3. The LSM for grain boundary detection
Fig. 1. Four possible cases of position of curve.
The LSM algorithm of the automatic segmentation consists of two steps: initializing the curve and then evolving it by solving the level set equation. To introduce the algorithm, we consider the original microphotograph of Daqing Mountain as shown in Fig. 3(a). For the sake of simplicity, we, first of all, convert the original image into a gray image in Fig. 3(b). Detecting the grain boundary on the gray image gives satisfactory results which are shown in Figs. 4 and 5.
ARTICLE IN PRESS 270
B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
Y
=0
X
Y
=0
Fig. 3. (a) Original color microphotograph. (b) Gray microphotograph.
X
Fig. 2. (a) Initial level set and corresponding curve. (b) Evolution of level set and corresponding curve.
3.1. Initial curve selection
scheme which is a simple and efficient method for solving PDE. Denote h ¼ 1 the space step (one pixel), Dt the time step and let ðxi ; yj Þ ¼ ðih; jhÞ for 1pipM; 1pjpN, where n Nh Mh is the size of image. Let fi;j ¼ fðxi ; yj ; nDtÞ be an 0 approximation to fðx; y; tÞ for nX0 and f ¼ f0 . Then we nþ1 compute f as follows: n n n fnþ1 ¼ fi;j þ nt j,fi;j j½l1 ðu0;i;j c1 ðf ÞÞ2 i;j n
þ l2 ðu0;i;j c2 ðf ÞÞ2 , n
The first step involved in the LSM is to choose an initial curve. The initial curve can be anywhere in the image and of an arbitrary shape and size. In the LSM, we choose the initial curve to be a contour fðx; y; 0Þ which is a circle with a radius 18 length of the microphotograph. The circle is located in the center of the image. The initial level set is computed as a signed distance function f, the distance between every point in the image and the initial contour. The sign of a pixel is defined to be negative (positive) if it is inside (outside) the contour. The pixel on the initial curve is zero. The boundary of the region is just the curve C.
(4)
n
where c1 ðf Þ and c2 ðf Þ are the averages inside and outside the curve C in the nth iteration, respectively. The algorithm used to solve the PDE (see Appendix A) stops n nþ1 when f ¼ f or when it reaches the maximum number of iterations. The algorithm is implemented in Matlab. We set l1 ¼ l2 ¼ 1 as usual and nt ¼ 0:1. The maximum number of iterations is set to be 20. The algorithm converges in 90 s. The evolving contours are shown in Fig. 4. The final segmentation result is shown in Fig. 4(d). 3.3. Post-processing
3.2. Level set evolution After the initial curve has been selected, the next step is to solve the level set equation numerically. To discretize the equation in f, we use an explicit finite difference
The contour obtained above is very complicated, making it impracticable for direct strain analysis. Consequently, we need to select suitable grains for strain analysis. But making the selection is not an easy task because the noise disturbs the segmented results and
ARTICLE IN PRESS B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
271
Fig. 4. Segmentation of microphotograph, l1 ¼ 1, l2 ¼ 1, nt ¼ 0:1: (a) Initial contour, (b) five iterations, (c) 10 iterations, and (d) 20 iterations.
Fig. 5. (a) Grains identified by LSM in color microphotograph. (b) Contours of grains identified by LSM in corresponding binary image.
Fig. 6. (a) Hand-drawn boundaries of grains in color microphotograph. (b) Hand-drawn boundaries of grains in corresponding binary image.
there are too many other objects. Hence, we set constraints so that only the useful objects are selected. Intensity, area, and shape are dominant factors in grain identification. The LSM obtains the boundaries of the
objects on the basis of their intensity; so the criteria based on the area and those based on shape information should be adopted to identify the candidate grains and to abandon the unsuitable grains. Since the contour obtained
ARTICLE IN PRESS 272
B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
critical criterion for high strain mylonites. As for the flake minerals, we discard the objects whose ratio between the major axis and the minor axis is greater than the three times average. From what we have discussed in Section 3, we can summarize the principal steps involved in the LSM algorithm as follows: Step 1: Choose an initial curve. Step 2: Compute the signed distance function as the initial level set. Step 3: Solve the PDE and obtain the evolving curve.
by the LSM is a closed curve, the objects with an area ranging between two thresholds k1 and k2 are selected. In our implementation, we select the objects with an area between 1000 and 2000 pixels for the image in Fig. 3. The conglutination of several grains or other minerals such as mica may be segmented as an object with a large area. But quartz or something meaningless may result in contour with small area. The objects with either a large or a small area can be discarded after a simple calculation of the areas. As for the shape information, the ratio between the major axis and the minor axis is adopted as another
500
800
Centroid_Y (in pixels)
Centroid_X (in pixels)
400 Hand Drawn
Hand Drawn
600
400 r = 1.00
200
300 200 r = 0.99
100
0
0 0
200
400 LSM
600
800
0
90
1200
Orientation (in degrees)
100
200 300 LSM
400
500
Area (in pixels)
900
30
Hand Drawn
Hand Drawn
60
0 -30
r = 0.79
600 r = 0.095
300 -60 -90
0 -90 -60 -30
0 30 LSM
60
90
0
80
300
600 LSM
900
1200
30 Major Axis (in pixels)
Minor Axis (in pixels)
Hand Drawn
Hand Drawn
60
40
20
20
10
r = 0.96
r = 0.96
0
0 0
20
40 LSM
60
80
0
10
20 LSM
Fig. 7. Comparisons of features of LSM vs. hand-drawn regions.
30
ARTICLE IN PRESS B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
Step 4: Check whether the termination criterion term is met. If not, repeat Steps 2–4. Step 5: Remove the objects which are unsuitable for strain analysis.
LSM 40
4. Methodology for validation of the LSM
30
Frequency
The results produced by the LSM are shown in Fig. 5. The identified regions are similar to the boundaries of the grains. To further validate the LSM, we compare the results produced by the LSM and by the hand-drawn method. The hand-drawn boundary is shown in Fig. 6. A reliable strain calculation by the Rf =f method requires large scale examples. The photomicrograph in Fig. 3 is of such a high resolution that the effects of the LSM can be examined visually. This photomicrograph has more than 70 clasts, which are enough for reliable strain analysis. The statistical method is used to check correspondence between the various methods. Since the aim of boundary detection is to carry out strain analysis from the information extracted from the clasts, the features and the strain analysis results should be compared.
20
10
0 -80
-60
-40
-20 0 Orientation
20
40
60
40
60
80
Hand Drawn 40
4.1. Comparison of features
30
Frequency
Grains have many features, some of which are centroid, area, major and minor axes, and orientation. In order to validate the LSM, we compare the features mentioned above. We choose the diagonal line y ¼ x to represent a perfect match between two sets of measurements (see Fig. 7). Also, we use histograms to compare orientations (see Fig. 8). In Fig. 7, we can see that there are roughly equal numbers of points on either side of the line in all the six cases. This indicates that two methods do not systematically over- or underestimate the parameters with each other. Moreover, most points lie close to the line, indicating that the results produced by the LSM are accurate. The centroids exacted by the two methods are almost the same. The sampling is symmetrical in one sample. As for orientation, most of the points are distributed between 20 and 20 , and the correlation is relatively low. This is caused mainly by two factors. One is the artificial separation of 90 and 90 , which are in fact the same. If we adjust these points, we will make the correlation substantially high. The other is the deformation intensity. The sample, especially the one drawn by hand, is seriously deformed that the orientation is difficult to be identified. The major axis and the minor axis have the same correlation, but the latter deviates more from the diagonal line. This is caused by the minor axis, between 5 and 30 pixels (see Fig. 7(f)), which is shorter than the major axis, between 10 and 70 pixels (see Fig. 7(e)). In Fig. 8, we can get the range of the orientation. All of them are relatively concentrated and closed, indicating that they are similar.
273
20
10
0 -60
-40
-20
0
20
Orientation Fig. 8. Frequency histograms of orientation of LSM and hand-drawn regions.
Table 1 A comparison of finite strain estimates (Rs ; y; Ri ) from data acquired manually and from those derived from the LSM
Hand drawn The LSM
Rs
y
Ri
2.21 2.4
6.2 9.18
2.03 2.06
4.2. Comparison of strain analysis Results of strain calculation using the Rf =f method (Yu and Zheng, 1984; Ailleres and Champenois, 1994) are shown in Table 1 in which Rs and Ri estimated from both
the manual and the LSM data are virtually identical. The orientation angle y is not quite uniform as demonstrated in the previous subsection. The remarkable similarity among the strain analysis results can be attributed to the
ARTICLE IN PRESS 274
B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
fact that strain calculation is an averaging procedure and there is no relative bias between the two methods. Thus, combining this method with the averaging procedure can cancel relative errors. The comparisons of strain analysis using the two methods indicate that the LSM is feasible and accurate for strain analysis.
5. Discussion and conclusions In this paper, we have presented a new automated approach for segmenting microphotographs by the LSM. This technique is used to identify the objects whose intensity is different from the background. The LSM is a very elegant method implemented by using the finite difference method. We have demonstrated that the LSM can identify the clasts in a complex microphotograph. Simple post-processing is required to select the proper objects for strain analysis. Compared with other methods, the LSM has the following three advantages: automated segmentation, closure of the boundary of the grain, and simple post-processing that does not change the shape of the boundary. The major advantage of the LSM algorithm is that it does not require human intervention in segmentation. In the LSM, the initial curve can be arbitrary in both shape and location in the image without affecting the final results. Furthermore, the curve evolves according to the energy functional which only depends on the structure of the image. Other algorithms require human intervention, such as seeding (Roy Choudhury et al., 2006) Although the seeding procedure can be replaced by other algorithms, those algorithms are usually complicated (Zhou et al., 2004). Another remarkable advantage of the LSM is that the boundary of every grain, inherited from the property of the level set, is a closed curve. Some features of the data, such as area and centroid, can be easily extracted. In order to ensure the accuracy of the measurements as well as of the strain analysis, geologists prefer the closed contour to the open boundary as segmentation result. The post-processing procedure in the LSM is very natural for strain analysis. The constraints exclude the unsuitable grains, which are in contrast with the unbiased automated segmentation. The post-processing procedure can be easily implemented without changing the shape of the grain boundary. The constraints, mainly derived from the information such as area and shape, are determined by the prior knowledge of the mineral. The post-processing is the only human intervention in the whole process. However, the LSM has certain limitations. The first limitation is that the LSM can only segment the bright grains in the image. Grains of mid-gray scale level are all missed. Inaccurate segmentation may occur with the gray grains with bright inclusion. Still, pre-processing filter, such as closing by reconstruction, can be used to improve the performance of the LSM. Another limitation may be the selection of the segmented objects. Although the initial segmented results can show all the clasts, even the very small meaningless ones in the image, the selection is
a difficult task. The criterion we choose is but an attempt. In the future, the color information and more prior knowledge of the clasts may offer a better criterion. The third limitation is that when we have two grains of similar color placed side by side, the LSM may treat them as one large grain. Thus, the conglutination occurs. The conglutination is caused either by the fact of the regions having subtle difference or by their being connected by a path which may be perceptible with the LSM, so other constraints become necessary to avoid such mistakes. Some of the limitations mentioned above can be resolved by using the color information of the images taken with different polarization setups. But quartz with the c-axis perpendicular to the section is always dark. This is a serious challenge for the LSM. The CPU time needed to run the LSM on a workstation is about several minutes, because solving the level set equation takes time and the regions detected consist of all objects, instead of the selected ones, within the microphotograph. It is still important to improve the efficiency of the LSM by adopting an advanced numerical scheme on fast computers. The challenges and limitations motivate us to improve our method in future work.
Acknowledgments This paper was jointly supported by the National Science Foundation of China (10801061, 10531040 and 10671082) and the Major State Basic Research Development Program of China (973 Program No. G2005CB422107), doctor fund of Henan Polytechnic University. The authors would like to thank the anonymous reviewers and the editor for their constructive suggestions that helped to improve the manuscript. The authors also would like to express their deep thanks to Professor Jingxue Yin and Professor Yongkui Zou, under whose guidance this paper was completed.
Appendix A. The level set algorithm The level set algorithm is as follows: 0 0 Step 1. Calculate the initial level set f : fi;j is the distance between the pixel ði; jÞ and the initial level set. 0 The sign of fi;j is negative when the pixel is inside the curve and vice versa. If the initial curve is a circle with 0 center ðx0 ; y0 Þ and radius r, fi;j can be determined as follows: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f0i;j ¼ ði x0 Þ2 þ ðj y0 Þ2 r. Step 2. Calculate the two parameters c1 and c2 : The means of the pixels inside and that of those outside the curve, respectively. Step 3. Calculate the gradient norm j,fj: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! !2 u u qfn 2 qfni;j i;j n t j,fi;j j ¼ þ , qx qy
ARTICLE IN PRESS B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
where
qfni;j qx
¼
fniþ1;j
2
fni1;j
;
qfni;j qy
¼
fni;jþ1
2
fni;j1
.
nþ1
Step 4. Calculate fi;j : n n n fnþ1 ¼ fi;j þ nt j,fi;j j½l1 ðu0 ði; jÞ c1 ðf ÞÞ2 i;j n
þ l2 ðu0 ði; jÞ c2 ðf ÞÞ2 .
(5)
Step 5. Check the termination criterion: If the termination criterion is met, f ¼ 0 is the final curve; otherwise set n ¼ n þ 1 and repeat Steps 2–5. Appendix B. Supplementary data Supplementary data associated with this article can be found in the online version at 10.1016/j.cageo.2008. 05.006.
References Ailleres, L., Champenois, M., 1994. Refinements to the Fry method (1979) using image processing. Journal of Structural Geology 16 (9), 1327–1330. Bartozzi, M., Boyle, A.P., Prior, D.J., 2000. Automated grain boundary detection and classification in orientation contrast images. Journal of Structural Geology 22 (11), 1569–1579.
275
Chan, T., Vese, L., 2001. Active contours without edges. IEEE Transactions on Image Processing 10, 266–277. Goodchild, J., Fueten, F., 1998. Edge detection in petrographic images using the rotation polarizer stage. Computers & Geosciences 24, 745–751. Heilbronner, R., 2000. Automatic grain boundary detection and grain size analysis using polarization micrographs or orientation images. Journal of Structural Geology 22 (15), 969–981. Kass, M., Witkin, A., Terzopoulos, D., 1988. Snakes: active contour models. International Journal of Computer Vision 1, 321–331. Lumbreras, F., Serrat, J., 1996. Segmentation of petrographical images of marbles. Computers & Geosciences 22, 547–558. Mlynarczuk, M., 1999. Some remarks on the application of image analysis and image processing for the description of the geometrical structures of rock. Physicochemical Problems of Mineral Processing 33, 107–116. Obara, B., Kozˇusˇni´kova´, A., 2007. Utilisation of the image analysis method for the detection of the morphological anisotropy of calcite grains in marble. Computational Geosciences 11 (4), 275–281. Osher, S., Sethian, J.A., 1988. Fronts propagation with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics 79, 12–49. Ramsay, J.G., 1967. Folding and Fracturing of Rocks. McGraw-Hill, New York, 568pp. Roy Choudhury, K., Patrick, A.M., Mulchrone, K.F., 2006. Automated grain boundary detection by CASRG. Journal of Structural Geology 28 (3), 363–375. Starkey, J., Samantaray, A.K., 1993. Edge detection in petrographic images. Journal of Microscopy 172 (3), 263–266. Yu, H., Zheng, Y., 1984. A statistical analysis applied to the Rf =f method. Tectonophysics 110, 151–155. Zhou, Y., Starkey, J., Mansinha, L., 2004. Segmentation of petrographic images by integrating edge detection and region growing. Computers & Geosciences 30, 817–831.
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
Automated grain boundary detection using the level set method$ Bibo Lu a,b, Min Cui c, Qiang liu a, Yangang Wang a,d, a
College of Mathematics, Jilin University, Changchun 130012, PR China School of Computer Science and Technology, Henan Polytechnic University, Jiaozuo 454000, P.R. China c Basin and Reservoir Research Centre, China University of Petroleum, Beijing 102249, PR China d Supercomputing Center, Computer Network Information Center, Chinese Academy of Sciences, Beijing 100190, P.R. China b
a r t i c l e in fo
abstract
Article history: Received 9 February 2007 Received in revised form 7 May 2008 Accepted 14 May 2008
In this paper, we detect the grain boundary by using a novel automated algorithm which is based on the level set method (the LSM). The LSM starts with only one initial closed curve which can be arbitrary in both shape and location in the microphotograph. The curve splits and merges repeatedly until it reaches the boundaries of the objects. The evolution depends on the structure of the image. This method requires simple and feasible constraints which ensure the exclusion of the objects unsuitable for reliable measurements and strain analysis. The LSM has the following three advantages: (1) with a single image as input, segmentation is automated without human intervention, (2) the grain boundary detected by the LSM is a closed curve, and (3) simple post-processing maintains the shape of the boundary. The segmentation and strain analysis by the LSM are demonstrated and the results are compared with those from hand-drawn method. & 2008 Elsevier Ltd. All rights reserved.
Keywords: Grain boundary detection The level set method Automated image segmentation Strain analysis
1. Introduction Strain analysis plays an important role in the study of structural geology, especially in the investigation of the tectonic history of a region. In theory, high sample intensity is necessary to provide enough raw data for strain analysis. However, in practice, the samples used for strain analysis seem to be insufficient, partly because collecting raw data for strain analysis is laborious and time consuming. Even though many methods have been developed, efficiency and accuracy are yet to be attained. In this paper, we focus on the problem of extracting data from a single digital microphotograph for strain analysis. We detect the grain boundary by using a novel automated algorithm which is based on the level set
$ Code available from server at http://www.iamg.org/CGEditor/index. htm. Corresponding author. Tel./fax: +86 43185167051. E-mail addresses: [email protected] (B. Lu), [email protected] (Y. Wang).
0098-3004/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2008.05.006
method (the LSM). Our analysis of the microphotograph of the mylonites of Daqing Mountain in Inner Mongolia, China, indicates that the LSM enables us to extract the data which are necessary for an accurate strain analysis. Traditional methods of measurement required sustained use of the polarizing microscope with skilled manipulation of the rotating stage and knowledge of the use of various graticules (Ramsay, 1967). An automatic method proposed by Goodchild and Fueten (1998) constructed a composite data set from a number of sequential images obtained from the polarizing microscope with automatic rotation stages. By using some digital images of the grain slice, geologists were able to detect the boundaries of clasts manually. They did this with the help of image processing software, such as Photoshop; and it is known that once the boundaries are obtained, other methods can be employed to make the measurements required for strain analysis. To improve the efficiency of the boundary detection, various semi-automated and automated segmentation methods have been proposed. Most of them are based
ARTICLE IN PRESS 268
B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
on image processing technology, particularly, image segmentation and boundary detection algorithms. These methods can be classified into two main categories: edgedriven methods and region-driven methods. The methods in the first category try to detect the grain boundary directly. Mlynarczuk (1999) proposed an algorithm based on the watershed method for segmenting grains from polarizing microscope images of dolomite and quartzite. Some traditional operators, such as the Robert operator and the Canny operator, were used to extract boundaries defined by those points with large gradient magnitude. Starkey and Samantaray (1993) combined the outcomes of a number of sequential filtered images to produce a final edge image. Lumbreras and Serrat (1996) proposed a two-step strategy to segment an image by performing an over-segmentation of the image and then using a region-merging procedure to get the final boundary. The key point of the region-merging procedure is to estimate two parameters from the information of a sequence of images of the same sample. Using a series of image processing routines, Goodchild and Fueten (1998) extracted the edges from a gradient image obtained by passing a simple gradient operator over images captured during the sampling process. Heilbronner (2000) developed an automatic method for boundary detection and grain size analysis by using polarizing micrographs or orientation images. Bartozzi et al. (2000) employed a detection-filtering algorithm to automatically rebuild a real boundary net, where boundaries are defined as high gradient features on an orientation contrast images. Recently, Obara and Kozˇusˇni´kova´ (2007) proposed a new method of detecting the calcite grains in marble. Based on a set of colored images filtered by a morphological alternating sequence filter, a composite image is constructed and the final grain boundary is determined by watershed segmentation technique. The initial boundary produced by the edge detection methods suffers from some imperfections. Such a boundary, which may be open or discontinuous, does not coincide with the real boundary. Consequently, a post-processing method, such as mathematical morphology, is required in order to obtain a closed boundary. However, the subjectivity of the postprocessing makes these methods unsuitable for strain analysis (Roy Choudhury et al., 2006). The region-driven methods are based mainly on a popular image segmentation technique, i.e., the region growing method, which identifies the region of the grain by absorbing all the points similar to the seed point. The major advantage of the region growing approach is that the boundary of the region is closed and this is preferred by both measurements and strain analysis. A critical step in region growing method is to select the proper seed, which often involves human interventions in order to help the seed grow. Sometimes, some holes occur in the identified clast. The size of the hole depends on the distribution of the noise and the area of the heterogeneity of color information within the clast. While the small holes can be filled by employing simple morphological operations, the large ones cannot, and remain almost the same as before. Post-processing affects the extracted boundaries and features, such as major and minor axes,
orientation, and centroid. A recent method proposed by Roy Choudhury et al. (2006), CASRG, identified a region which is a set of points of the grain. They selected some seeds manually and chose the optimal threshold separately for each grain. This method is very efficient and is accurate for low and high strain samples. However, this method requires that all the seeds be chosen manually, using the computer mouse. But where large or high stress samples are involved, selecting seeds by clicking the mouse is onerous. Zhou et al. (2004) proposed a segmentation method of petrographic images by integrating the edge detection method with that of the region growing. They employed a boundary detection method to get the initial edge whose information is used for the automated seed selection in the region growing method. The methods mentioned above are mainly based on the traditional edge detection and image segmentation technology in image processing. Human interventions are required during the segmentation, and the objectivity of the data is inevitably corrupted by the post-processing. In this paper, we propose a new automated grain boundary detection method: LSM. The LSM requires a little mathematical knowledge of partial differential equation (PDE). We show that the LSM produces impressive results in grain boundary detection as well as in strain analysis. The remaining part of the paper is organized as follows. Section 2 briefly introduces the LSM. In Section 3, we discuss the algorithm for evolving the level set and feasible post-processing. In Section 4, we demonstrate the performance of the LSM. Section 5 is discussion and conclusions. 2. The LSM The LSM aims at identifying grain boundaries and representing them as closed outlines. Unlike CASGR (Roy Choudhury et al., 2006) which manually identifies all the interesting objects one by one, the LSM identifies all the objects that stand out from average gray levels simultaneously. The LSM is a PDE-based segmentation model with the level set formulation. Recently, methods based on PDE have been widely used in image processing and analysis; and impressive results have been obtained by various PDE models. Using level set, Chan and Vese (2001) proposed a successful segmentation model (C–V model) that is based on the classical active contour model (Kass et al., 1988). The essential idea in using the C–V model is to evolve an initial contour governed by a PDE until it stops on the edges of the objects. To introduce the LSM, we first of all review the C–V model and then formulate the model in terms of the level set. The detailed algorithm will be presented in the next section. 2.1. The C–V model Let us begin our discussion with a simple case. Assume that the image u0 consists of one object and background with constant intensity respectively. The value of the
ARTICLE IN PRESS B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
object is u1 and the background is u2 . Denote the boundary of the object by C 0. Then we have u0 u1 inside the object and u0 u2 outside the object. Then let us consider the following fitting term: Z ju0 ðx; yÞ c1 j2 dx dy F 1 ðCÞ þ F 2 ðCÞ ¼ insideðCÞ Z þ ju0 ðx; yÞ c2 j2 dx dy, outsideðCÞ
where C is any other variable curve, and the constants c1 ; c2 are the averages of u0 inside and outside C, respectively. In this simple example, four possible cases are illustrated in Fig. 1. In the first case, the curve C is outside the object. Thus, we have F 1 ðCÞ40 and F 2 ðCÞ ¼ 0 (see Fig. 1(a)). In the second case, the curve C is inside the object. Thus, we have F 1 ðCÞ ¼ 0 but F 2 ðCÞ40 (see Fig. 1(b)). In the third case, since the curve C is both inside and outside the object, we have F 1 ðCÞ40 and F 2 ðCÞ40 (see Fig. 1(c)). In the fourth case, it is only when the curve C coincides with C 0 , the boundary of the object, that the above fitting term can reach the minimizer (see Fig. 1(d)); that is
269
generates two other problems: (1) how to represent the curve used to trace the boundary of the object and (2) how to evolve the curve until it stops on the boundary of the object. The LSM provides a good solution to the first problem. 2.2. Level set formulation of the model The LSM (Osher and Sethian, 1988) represents a curve as the zero level set of a high-dimensional function and allows for major simplifications. The LSM offers a natural representation of the contour of the object. Thus, it is capable of dealing with complex structures with many advantages, especially when the curve undergoes complicated topological changes. The LSM has the following advantages: (1) It is easy to implement numerically. (2) The outline of the object is a closed curve. (3) It can easily deal with the topological changes. such as merging and splitting. (4) It is able to compute geometric quantities easily.
inf C ðF 1 ðCÞ þ F 2 ðCÞÞ 0 F 1 ðC 0 Þ þ F 2 ðC 0 Þ. Therefore, the energy functional Fðc1 ; c2 ; CÞ is defined by Z ju0 ðx; yÞ c1 j2 dx dy Fðc1 ; c2 ; CÞ ¼ l1 insideðCÞ Z ju0 ðx; yÞ c2 j2 dx dy, (1) þ l2 outsideðCÞ
where l1 40 and l2 40 are fixed parameters. Now we convert the segmentation of the image into a minimization problem: inf Fðc1 ; c2 ; CÞ.
c1 ;c2 ;C
(2)
Solving the optimization problem is equivalent to segmenting the image; and the attempt to solve the problem
In the LSM, a closed curve is the zero level set of a function in a high-dimensional space. Let fðx; y; tÞ be a level set function defined on the whole image u0 . The level set function is determined by the position of the points. fðx; y; tÞo0 if the point is inside the curve and fðx; y; tÞ40 if the point is outside the curve. The interface between the region fðx; y; tÞ40 and the region fðx; y; tÞo0 is the curve C, which we call the zero level set because f ¼ 0 on this curve. The zero level set can be considered as the interface between a surface and a plane. The level set evolves in time under the control of the appointed model. Fig. 2 illustrates the change of the level set function and the corresponding curve. In our method, we use the level set to capture the contour of the object. With the level set, the minimization problem can be deduced into the following equation after some mathematical manipulations (see Chan and Vese, 2001 for more details):
qf ¼ j,fj½l1 ðu0 c1 Þ2 þ l2 ðu0 c2 Þ2 , qt
(3)
=qffiffi t is the partial derivative of f with respect to t where qfp , j,fj ¼ ðqf=qxÞ2 þ ðqf=qyÞ2 . Given an initial level set, the level set f evolves according to Eq. (3). The regularization term in the original model is dropped to maintain the primal grain boundary. 3. The LSM for grain boundary detection
Fig. 1. Four possible cases of position of curve.
The LSM algorithm of the automatic segmentation consists of two steps: initializing the curve and then evolving it by solving the level set equation. To introduce the algorithm, we consider the original microphotograph of Daqing Mountain as shown in Fig. 3(a). For the sake of simplicity, we, first of all, convert the original image into a gray image in Fig. 3(b). Detecting the grain boundary on the gray image gives satisfactory results which are shown in Figs. 4 and 5.
ARTICLE IN PRESS 270
B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
Y
=0
X
Y
=0
Fig. 3. (a) Original color microphotograph. (b) Gray microphotograph.
X
Fig. 2. (a) Initial level set and corresponding curve. (b) Evolution of level set and corresponding curve.
3.1. Initial curve selection
scheme which is a simple and efficient method for solving PDE. Denote h ¼ 1 the space step (one pixel), Dt the time step and let ðxi ; yj Þ ¼ ðih; jhÞ for 1pipM; 1pjpN, where n Nh Mh is the size of image. Let fi;j ¼ fðxi ; yj ; nDtÞ be an 0 approximation to fðx; y; tÞ for nX0 and f ¼ f0 . Then we nþ1 compute f as follows: n n n fnþ1 ¼ fi;j þ nt j,fi;j j½l1 ðu0;i;j c1 ðf ÞÞ2 i;j n
þ l2 ðu0;i;j c2 ðf ÞÞ2 , n
The first step involved in the LSM is to choose an initial curve. The initial curve can be anywhere in the image and of an arbitrary shape and size. In the LSM, we choose the initial curve to be a contour fðx; y; 0Þ which is a circle with a radius 18 length of the microphotograph. The circle is located in the center of the image. The initial level set is computed as a signed distance function f, the distance between every point in the image and the initial contour. The sign of a pixel is defined to be negative (positive) if it is inside (outside) the contour. The pixel on the initial curve is zero. The boundary of the region is just the curve C.
(4)
n
where c1 ðf Þ and c2 ðf Þ are the averages inside and outside the curve C in the nth iteration, respectively. The algorithm used to solve the PDE (see Appendix A) stops n nþ1 when f ¼ f or when it reaches the maximum number of iterations. The algorithm is implemented in Matlab. We set l1 ¼ l2 ¼ 1 as usual and nt ¼ 0:1. The maximum number of iterations is set to be 20. The algorithm converges in 90 s. The evolving contours are shown in Fig. 4. The final segmentation result is shown in Fig. 4(d). 3.3. Post-processing
3.2. Level set evolution After the initial curve has been selected, the next step is to solve the level set equation numerically. To discretize the equation in f, we use an explicit finite difference
The contour obtained above is very complicated, making it impracticable for direct strain analysis. Consequently, we need to select suitable grains for strain analysis. But making the selection is not an easy task because the noise disturbs the segmented results and
ARTICLE IN PRESS B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
271
Fig. 4. Segmentation of microphotograph, l1 ¼ 1, l2 ¼ 1, nt ¼ 0:1: (a) Initial contour, (b) five iterations, (c) 10 iterations, and (d) 20 iterations.
Fig. 5. (a) Grains identified by LSM in color microphotograph. (b) Contours of grains identified by LSM in corresponding binary image.
Fig. 6. (a) Hand-drawn boundaries of grains in color microphotograph. (b) Hand-drawn boundaries of grains in corresponding binary image.
there are too many other objects. Hence, we set constraints so that only the useful objects are selected. Intensity, area, and shape are dominant factors in grain identification. The LSM obtains the boundaries of the
objects on the basis of their intensity; so the criteria based on the area and those based on shape information should be adopted to identify the candidate grains and to abandon the unsuitable grains. Since the contour obtained
ARTICLE IN PRESS 272
B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
critical criterion for high strain mylonites. As for the flake minerals, we discard the objects whose ratio between the major axis and the minor axis is greater than the three times average. From what we have discussed in Section 3, we can summarize the principal steps involved in the LSM algorithm as follows: Step 1: Choose an initial curve. Step 2: Compute the signed distance function as the initial level set. Step 3: Solve the PDE and obtain the evolving curve.
by the LSM is a closed curve, the objects with an area ranging between two thresholds k1 and k2 are selected. In our implementation, we select the objects with an area between 1000 and 2000 pixels for the image in Fig. 3. The conglutination of several grains or other minerals such as mica may be segmented as an object with a large area. But quartz or something meaningless may result in contour with small area. The objects with either a large or a small area can be discarded after a simple calculation of the areas. As for the shape information, the ratio between the major axis and the minor axis is adopted as another
500
800
Centroid_Y (in pixels)
Centroid_X (in pixels)
400 Hand Drawn
Hand Drawn
600
400 r = 1.00
200
300 200 r = 0.99
100
0
0 0
200
400 LSM
600
800
0
90
1200
Orientation (in degrees)
100
200 300 LSM
400
500
Area (in pixels)
900
30
Hand Drawn
Hand Drawn
60
0 -30
r = 0.79
600 r = 0.095
300 -60 -90
0 -90 -60 -30
0 30 LSM
60
90
0
80
300
600 LSM
900
1200
30 Major Axis (in pixels)
Minor Axis (in pixels)
Hand Drawn
Hand Drawn
60
40
20
20
10
r = 0.96
r = 0.96
0
0 0
20
40 LSM
60
80
0
10
20 LSM
Fig. 7. Comparisons of features of LSM vs. hand-drawn regions.
30
ARTICLE IN PRESS B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
Step 4: Check whether the termination criterion term is met. If not, repeat Steps 2–4. Step 5: Remove the objects which are unsuitable for strain analysis.
LSM 40
4. Methodology for validation of the LSM
30
Frequency
The results produced by the LSM are shown in Fig. 5. The identified regions are similar to the boundaries of the grains. To further validate the LSM, we compare the results produced by the LSM and by the hand-drawn method. The hand-drawn boundary is shown in Fig. 6. A reliable strain calculation by the Rf =f method requires large scale examples. The photomicrograph in Fig. 3 is of such a high resolution that the effects of the LSM can be examined visually. This photomicrograph has more than 70 clasts, which are enough for reliable strain analysis. The statistical method is used to check correspondence between the various methods. Since the aim of boundary detection is to carry out strain analysis from the information extracted from the clasts, the features and the strain analysis results should be compared.
20
10
0 -80
-60
-40
-20 0 Orientation
20
40
60
40
60
80
Hand Drawn 40
4.1. Comparison of features
30
Frequency
Grains have many features, some of which are centroid, area, major and minor axes, and orientation. In order to validate the LSM, we compare the features mentioned above. We choose the diagonal line y ¼ x to represent a perfect match between two sets of measurements (see Fig. 7). Also, we use histograms to compare orientations (see Fig. 8). In Fig. 7, we can see that there are roughly equal numbers of points on either side of the line in all the six cases. This indicates that two methods do not systematically over- or underestimate the parameters with each other. Moreover, most points lie close to the line, indicating that the results produced by the LSM are accurate. The centroids exacted by the two methods are almost the same. The sampling is symmetrical in one sample. As for orientation, most of the points are distributed between 20 and 20 , and the correlation is relatively low. This is caused mainly by two factors. One is the artificial separation of 90 and 90 , which are in fact the same. If we adjust these points, we will make the correlation substantially high. The other is the deformation intensity. The sample, especially the one drawn by hand, is seriously deformed that the orientation is difficult to be identified. The major axis and the minor axis have the same correlation, but the latter deviates more from the diagonal line. This is caused by the minor axis, between 5 and 30 pixels (see Fig. 7(f)), which is shorter than the major axis, between 10 and 70 pixels (see Fig. 7(e)). In Fig. 8, we can get the range of the orientation. All of them are relatively concentrated and closed, indicating that they are similar.
273
20
10
0 -60
-40
-20
0
20
Orientation Fig. 8. Frequency histograms of orientation of LSM and hand-drawn regions.
Table 1 A comparison of finite strain estimates (Rs ; y; Ri ) from data acquired manually and from those derived from the LSM
Hand drawn The LSM
Rs
y
Ri
2.21 2.4
6.2 9.18
2.03 2.06
4.2. Comparison of strain analysis Results of strain calculation using the Rf =f method (Yu and Zheng, 1984; Ailleres and Champenois, 1994) are shown in Table 1 in which Rs and Ri estimated from both
the manual and the LSM data are virtually identical. The orientation angle y is not quite uniform as demonstrated in the previous subsection. The remarkable similarity among the strain analysis results can be attributed to the
ARTICLE IN PRESS 274
B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
fact that strain calculation is an averaging procedure and there is no relative bias between the two methods. Thus, combining this method with the averaging procedure can cancel relative errors. The comparisons of strain analysis using the two methods indicate that the LSM is feasible and accurate for strain analysis.
5. Discussion and conclusions In this paper, we have presented a new automated approach for segmenting microphotographs by the LSM. This technique is used to identify the objects whose intensity is different from the background. The LSM is a very elegant method implemented by using the finite difference method. We have demonstrated that the LSM can identify the clasts in a complex microphotograph. Simple post-processing is required to select the proper objects for strain analysis. Compared with other methods, the LSM has the following three advantages: automated segmentation, closure of the boundary of the grain, and simple post-processing that does not change the shape of the boundary. The major advantage of the LSM algorithm is that it does not require human intervention in segmentation. In the LSM, the initial curve can be arbitrary in both shape and location in the image without affecting the final results. Furthermore, the curve evolves according to the energy functional which only depends on the structure of the image. Other algorithms require human intervention, such as seeding (Roy Choudhury et al., 2006) Although the seeding procedure can be replaced by other algorithms, those algorithms are usually complicated (Zhou et al., 2004). Another remarkable advantage of the LSM is that the boundary of every grain, inherited from the property of the level set, is a closed curve. Some features of the data, such as area and centroid, can be easily extracted. In order to ensure the accuracy of the measurements as well as of the strain analysis, geologists prefer the closed contour to the open boundary as segmentation result. The post-processing procedure in the LSM is very natural for strain analysis. The constraints exclude the unsuitable grains, which are in contrast with the unbiased automated segmentation. The post-processing procedure can be easily implemented without changing the shape of the grain boundary. The constraints, mainly derived from the information such as area and shape, are determined by the prior knowledge of the mineral. The post-processing is the only human intervention in the whole process. However, the LSM has certain limitations. The first limitation is that the LSM can only segment the bright grains in the image. Grains of mid-gray scale level are all missed. Inaccurate segmentation may occur with the gray grains with bright inclusion. Still, pre-processing filter, such as closing by reconstruction, can be used to improve the performance of the LSM. Another limitation may be the selection of the segmented objects. Although the initial segmented results can show all the clasts, even the very small meaningless ones in the image, the selection is
a difficult task. The criterion we choose is but an attempt. In the future, the color information and more prior knowledge of the clasts may offer a better criterion. The third limitation is that when we have two grains of similar color placed side by side, the LSM may treat them as one large grain. Thus, the conglutination occurs. The conglutination is caused either by the fact of the regions having subtle difference or by their being connected by a path which may be perceptible with the LSM, so other constraints become necessary to avoid such mistakes. Some of the limitations mentioned above can be resolved by using the color information of the images taken with different polarization setups. But quartz with the c-axis perpendicular to the section is always dark. This is a serious challenge for the LSM. The CPU time needed to run the LSM on a workstation is about several minutes, because solving the level set equation takes time and the regions detected consist of all objects, instead of the selected ones, within the microphotograph. It is still important to improve the efficiency of the LSM by adopting an advanced numerical scheme on fast computers. The challenges and limitations motivate us to improve our method in future work.
Acknowledgments This paper was jointly supported by the National Science Foundation of China (10801061, 10531040 and 10671082) and the Major State Basic Research Development Program of China (973 Program No. G2005CB422107), doctor fund of Henan Polytechnic University. The authors would like to thank the anonymous reviewers and the editor for their constructive suggestions that helped to improve the manuscript. The authors also would like to express their deep thanks to Professor Jingxue Yin and Professor Yongkui Zou, under whose guidance this paper was completed.
Appendix A. The level set algorithm The level set algorithm is as follows: 0 0 Step 1. Calculate the initial level set f : fi;j is the distance between the pixel ði; jÞ and the initial level set. 0 The sign of fi;j is negative when the pixel is inside the curve and vice versa. If the initial curve is a circle with 0 center ðx0 ; y0 Þ and radius r, fi;j can be determined as follows: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f0i;j ¼ ði x0 Þ2 þ ðj y0 Þ2 r. Step 2. Calculate the two parameters c1 and c2 : The means of the pixels inside and that of those outside the curve, respectively. Step 3. Calculate the gradient norm j,fj: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! !2 u u qfn 2 qfni;j i;j n t j,fi;j j ¼ þ , qx qy
ARTICLE IN PRESS B. Lu et al. / Computers & Geosciences 35 (2009) 267–275
where
qfni;j qx
¼
fniþ1;j
2
fni1;j
;
qfni;j qy
¼
fni;jþ1
2
fni;j1
.
nþ1
Step 4. Calculate fi;j : n n n fnþ1 ¼ fi;j þ nt j,fi;j j½l1 ðu0 ði; jÞ c1 ðf ÞÞ2 i;j n
þ l2 ðu0 ði; jÞ c2 ðf ÞÞ2 .
(5)
Step 5. Check the termination criterion: If the termination criterion is met, f ¼ 0 is the final curve; otherwise set n ¼ n þ 1 and repeat Steps 2–5. Appendix B. Supplementary data Supplementary data associated with this article can be found in the online version at 10.1016/j.cageo.2008. 05.006.
References Ailleres, L., Champenois, M., 1994. Refinements to the Fry method (1979) using image processing. Journal of Structural Geology 16 (9), 1327–1330. Bartozzi, M., Boyle, A.P., Prior, D.J., 2000. Automated grain boundary detection and classification in orientation contrast images. Journal of Structural Geology 22 (11), 1569–1579.
275
Chan, T., Vese, L., 2001. Active contours without edges. IEEE Transactions on Image Processing 10, 266–277. Goodchild, J., Fueten, F., 1998. Edge detection in petrographic images using the rotation polarizer stage. Computers & Geosciences 24, 745–751. Heilbronner, R., 2000. Automatic grain boundary detection and grain size analysis using polarization micrographs or orientation images. Journal of Structural Geology 22 (15), 969–981. Kass, M., Witkin, A., Terzopoulos, D., 1988. Snakes: active contour models. International Journal of Computer Vision 1, 321–331. Lumbreras, F., Serrat, J., 1996. Segmentation of petrographical images of marbles. Computers & Geosciences 22, 547–558. Mlynarczuk, M., 1999. Some remarks on the application of image analysis and image processing for the description of the geometrical structures of rock. Physicochemical Problems of Mineral Processing 33, 107–116. Obara, B., Kozˇusˇni´kova´, A., 2007. Utilisation of the image analysis method for the detection of the morphological anisotropy of calcite grains in marble. Computational Geosciences 11 (4), 275–281. Osher, S., Sethian, J.A., 1988. Fronts propagation with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. Journal of Computational Physics 79, 12–49. Ramsay, J.G., 1967. Folding and Fracturing of Rocks. McGraw-Hill, New York, 568pp. Roy Choudhury, K., Patrick, A.M., Mulchrone, K.F., 2006. Automated grain boundary detection by CASRG. Journal of Structural Geology 28 (3), 363–375. Starkey, J., Samantaray, A.K., 1993. Edge detection in petrographic images. Journal of Microscopy 172 (3), 263–266. Yu, H., Zheng, Y., 1984. A statistical analysis applied to the Rf =f method. Tectonophysics 110, 151–155. Zhou, Y., Starkey, J., Mansinha, L., 2004. Segmentation of petrographic images by integrating edge detection and region growing. Computers & Geosciences 30, 817–831.