Abstract
This paper proposes a new algorithm (DA3DED) for edge detection in 3D images. DA3DED is doubly adaptive because it is based on the adaptive algorithm EDAS-1 for detecting edges in functions of one variable and a second adaptive procedure based on the concept of projective complexity of a 3D image. DA3DED has been tested on 3D images that modelize real problems (composites and fractures). It has been much faster than the 1D edge detection algorithm for 3D images derived from EDAS-1.
1. Introduction
Three-dimensional imaging is applied in many topics such as medicine and materials science among others. 3D medical imaging models structures of the human body. This is important in many fields as image guided surgery, assessment of the quality of bones, and so forth; see, for example, [1–3]. The complex geometric structure of some materials (composites, foams, etc.) can be treated with 3D images: 3D images of materials provide data such as the 3D connectivity of a structure, distribution of particles, etc. These data can be used in simulations to compute macroscopic material properties, what can be considered the basis of virtual material design [4].
The detection of edges is an essential objective in image processing. Besides the applications aforementioned, it is useful in the analysis and study of 3D Satellite images, among others. There are a lot of procedures for determining 3D edges. These procedures can be classified in Direct Methods and Indirect Methods [5–7]. Among the direct methods we can find 1D edge detection methods (1D3DED) [8, 9], spatial difference filters [10–12], polynomial fitting method [13], and deformable models [14–20]. We can add to the direct methods the adaptive splitting methods, which are based on the adaptive splitting of the domain of the function guided by the value of an average integral. In [21], an algorithm (EDAS-3) to approximate the jump discontinuity set of functions defined on subsets of based in these methods was proposed. This algorithm overcomes most inconveniences presented by deformable models. It exhibits effective edge detection in 3D models with different interface topologies. Also, EDAS-3 can obtain the edges with a prefixed accuracy.
The present paper shows an accurate algorithm (DA3DED) that can obtain a point based representation of the jump discontinuity set of a 3D image. The algorithm is a type of 1D edge detection method (1D3DED) and adaptive splitting method. It begins by considering a region in the plane containing the projection of the image and a set of points in this region. Then, it considers straight lines perpendicular to the plane and computes the edge points lying on these straight lines. This task is accomplished using the 1D edge detector EDAS-1 [22]. If the number of edges detected in the region is high or the distance between them is small, the algorithm concludes that the region under consideration is complex and performs a subdivision process; otherwise, the region is not divided and the algorithm studies a contiguous zone. The process is repeated for the and planes. The edge points obtained are stored in a file. In this way, complex regions are crossed by many lines and simple regions are crossed by few lines. The resulting algorithm is doubly adaptive because the “straight line step” is performed in an adaptive way by EDAS-1 and the second step performs an adaptive process based on the complexity of the image. The algorithm uses less CPU time than previous algorithm [21]. We will give a C-like pseudocode description of the constituent algorithms of DA3DED.
Experimental results show a good behavior of DA3DED on practical 3D problems that model different types of composite materials and fractures.
The outline of the paper is as follows. Section 2 states some mathematical preliminary results and details the conceptual and algorithmic constituents of DA3DED: EDAS-1 algorithm, discrete and continuous representation of 3D images, projective complexity, subdivision procedures of discrete rectangles, and split criterion. It also describes the algorithm DA3DED. Section 3 shows computational experiments on models of complex materials. Section 4 provides some concluding remarks.
2. Materials and Methods
2.1. Mathematical Preliminaries
Let be a general piecewise continuous function in , a compact -interval.
We assume that , where is a series of connected sets with pairwise disjoint nonempty interiors and we call be the boundary of , .
Define . We call the set of points in with jump discontinuities. Achieving a good approximation of is our aim.
Let be a set of affinely independent vectors in and the d-simplex generated by them. Consider a function . We call , the unique affine map such that
It is proved at [21, 22] that the behavior of the next average integral depends on the continuity or lack of continuity of on :where is the Lebesgue measure of .
This result is the basis of the algorithm described in the next section.
For more details about the mathematical preliminaries of the problem, see [21, 22].
2.2. Constituents of the Algorithm DA3DED
DA3DED is based on several concepts and algorithms that are described in this subsection.
2.2.1. Algorithm EDAS-1
EDAS-1 is an algorithm to approximate the jump discontinuity set of functions of one variable. It is a particular case of the algorithm EDAS- for functions defined on subsets of . In this section, we give an outline of EDAS-. Its validity for 1, 2 has been established in [22]. The case has been studied in [21].
More details about the implementation and performance of EDAS-1 can be found in [22].
Next, the steps of the EDAS- algorithm are summarized.
Step 1. Consider the initial partition , of , where .
Step 2. Put into the set the good simplices of and put into the set the bad simplices of .
is a good simplex if and , where is the diameter of , is the maximum local error of the approximant , is the approximation error of the points in , and is a positive real parameter.
Step 3. Divide each bad simplex into two simplices, by splitting its largest edge. We have a new partition .
Step 4. Repeat Steps 2 and 3 until .
Step 5. is the searched partition. Choose the following subset of :where is the minimum magnitude of jump reported and is also a stopping criterion.
In the images studied in Section 3, a voxel is considered as an edge voxel if it contains some barycenter of the simplices in .
In practice, we have implemented the above algorithm using a binary tree whose leaves correspond to good simplices. In the case , the integrals have been computed using the Gauss-Legendre integration formula [23]. More details about the implementation and performance of EDAS-1 can be found in [22].
2.2.2. Continuous and Discrete Representations of a 3D Image
The algorithm DA3DED uses two different mathematical representations of a 3D image. It considers that the image is a function defined on a box in when it computes edge points using EDAS-1. This continuous representation is also useful to define the concept of “projective complexity”; we will see this in the next subsection. In the remaining steps, DA3DED uses the usual discrete representation. Figure 1 shows the difference between both representations in the case of a 2D image. Below we define these concepts.
(a)
(b)
Let and be positive integers with ; defineA 3D image is given by an array with indices . From a 3D image , one can build a piecewise constant function , defined on in the following way:
2.2.3. Projective Complexity of a 3D Image
Several complexity measures for 2D images have been proposed [24, 25]. Some definitions are based on statistical gray level features (global versus local histograms). Other definitions are based on spatial distributions of gray levels (edges, symmetry, and texture). Edges highlight image zones in which there are abrupt changes in luminosity level usually associated with surface discontinuities. The rationale is straightforward: images with more borders contain more objects (or objects with more facets), thus resulting in a greater perceived complexity. This has led to the following edge based measure of complexity. Complexity is measured by the number of edges per unit area in the image [26]. In the case of 3D images, we can extend this definition; for example, complexity can be measured by the number of edge voxels per unit volume. This definition is not suitable for adaptive algorithms based on 1D edge detection methods. Therefore, we propose a new definition consistent with the above one.
Consider a continuously defined function . Let and . Let be a fixed point in . Call the number of edge points of the function of one variable when . The projective complexity of on the set is defined as
Let be a positive integer. To approximate the above quantity, we subdivide into identical squares. Denote the center of the squares by , . Then, (6) can be approximated by
Observe that does not depend on the size of the square . We can consider as a vector with components and write (7) asIn each experiment, the value of has been almost constant. Then, using the norm equivalence, we can define the following measure of complexity:Although we have considered sets in the plane , similar definitions can be stated for sets in the planes or . The above definitions are applicable to 3D images because they can be extended to continuously defined functions using the procedure detailed in Section 2.2.2.
2.2.4. Subdivision and Selection Procedures
Consider a 3D image given by an array with indices .
The doubly adaptive algorithm performs split and point selection operations. These procedures are defined on the set of indices; that is, they consider discrete rectangles. First, we describe how a discrete rectangle is divided into four almost equal discrete subrectangles.
Consider the discrete rectangle with and . The subdivision is performed computing and ( denotes the floor function that returns the greatest integer less than or equal to ). The resulting discrete rectangles areSuppose now that we want to select points from , distributed regularly. The selected points are obtained as the setwhere the vectors and are given by the following algorithm: ; ; ; ; ; ; ; ; Observe that returns the integer nearest to . The above algorithm obtains a subset of whose points are distanced approximately by . If or , the subdivision is not possible.
2.2.5. Split Criterion
The set of indices of a 3D image is given byThe projections of this set on the coordinate planes are
Consider a set , . Suppose that is a subset of . The split criterion for is defined by the following procedure.
Split Criterion
Step 1. If return 0.
Step 2. Select points from , distributed regularly (Section 2.2.4). The result is the set .
Step 3. Consider the straight lines perpendicular to the plane containing , passing through each point in .
Step 4. Apply EDAS-1 to each one of these lines to obtain the number of edge points corresponding to each point ; see Figure 2. Compute the minimum distance between these edge points . In this step, all the edge points found are stored in a file.
Step 5. Compute
Step 6. If we call the limit of the maximum number of edge points and the limit of the minimum distance between edge points, then
if ) return 1; else return 0.
A discrete rectangle is subdivided when the measure of projective complexity corresponding to the discrete rectangle is greater than a fixed beforehand limit or when the minimum distance between edge points is less than a fixed threshold. In this last case, it is possible that the discontinuity surface is nonsmooth and this makes it necessary to consider small discrete rectangles to obtain a detailed description of the jump discontinuity set. Given a discrete rectangle , the split criterion is implemented with the function . takes the value 1 if the split criterion is satisfied and 0 in the contrary case.
2.3. Doubly Adaptive Algorithm
Consider the discrete rectangle . Let be a discrete rectangle contained in . If , we call a good discrete rectangle; otherwise, is called bad.
Doubly Adaptive Algorithm for 3D Edge Detection (DA3DED)
Step 1. Consider . Divide into four discrete rectangles using the procedure described in Section 2.2.4. Apply the split criterion to the resulting rectangles. The good rectangles are put into the set and the bad rectangles are put into the set .
Step 2. At each step we have a set of good discrete rectangles and a set of bad discrete rectangles . Divide each bad discrete rectangle into four discrete rectangles using the procedure described in Section 2.2.4. Test whether these children are good or bad to obtain the sets and .
Step 3. The algorithm stops if . If the stopping criterion is not satisfied, go to Step 2.
Step 4. Repeat the above steps for and .
DA3DED has been implemented using a quadtree. Below we detail this algorithm.
2.3.1. Algorithm to Generate the Quadtree
In this subsection, we provide a pseudocode of the algorithm used by DA3DED to generate the quadtree. We start by defining a structure of type node (each node is associated with a discrete rectangle ): struct node int ; int ; int ; int ; int ; int ; int ; int ; ,where, , , , and are the coordinates of the corners of the discrete rectangle: : indicator variable that takes the value 1 if the algorithm has applied the split criterion to the rectangle and 0 otherwise. : number of the node parent of the current one. : th child of the current node (−1 in case of leaf nodes). is the left child and is the right child. : determines the subset of used by the split criterion.
The algorithm uses alternately two distinct values of : and , in order to avoid the repetition of computations. In practice, these values are almost equal. We denote by the discrete rectangle associated with the current node .
Algorithm to Generate the Quadtree. See Algorithm 1.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
3. Results and Discussion
In this section, we have tested DA3DED with several synthetic material models. The results have been compared with those obtained with other algorithms: 3D Sobel, 3D Prewitt, and 3D EDAS-1. The test models have been the following:(i)Internal inclusions within a material(ii)Short fiber composites(iii)Fracture models
Materials with internal inclusions are also called particle composites (filled materials). For example, internal pores in aluminium alloys can be obtained by introducing hydrogen in the melt metal. Since the solubility of the hydrogen decreases with the temperature, when the metal solidifies, the gas is rejected and forms gas pores which are easily visualized using X-ray tomography [27]. Icosahedral inclusions have been described in [28]. Short fiber composites are a kind of composite materials with short fibers of materials included in a matrix material. For example, glass or carbon fibers are used to reinforce polymers. The habitual fibers present a circular cross section. Nowadays composite research is being conducted with more exotic fibers: hollow glass fibers and triangular glass fibers. In the case of fibers with triangular cross sections it is difficult to compute the fiber orientations from 2D images. This makes necessary 3D imaging [29]. Finally, we have considered a model of fractured material. The interface of a fracture model is rather difficult to approximate because it presents acute dihedral angles. In fact, many procedures to detect 3D edges assume that the jump discontinuity set is smooth.
To test the algorithms, we have designed examples having different complexity in different zones. The box containing the composite has been divided into eight equal cubes. The inclusions and fibers have been randomly placed in each cube with the following distribution:(i)Spherical inclusions: 60 spheres (with radius 8) in one of the cubes, each one of the remaining parts contains 4 spheres(ii)Icosahedral inclusions: 60 icosahedra in one of the cubes, each one of the remaining parts contains 4 icosahedra(iii)Cylindrical fibers: 30 cylinders in one of the cubes, each one of the remaining parts contains 2 cylinders(iv)Prismatic fibers: 30 triangular prisms in one of the cubes, each one of the remaining parts contains 2 triangular prismsThe 3D images used in the experiments have been generated in two steps. First, we have considered a continuous function with value 1 inside and 0 outside the bodies. In the case of icosahedral inclusions, prismatic fibers, and fracture model, we have used the algorithm proposed in [30] to find the distance from a point to a polyhedron. Then, the above continuous function has been discretized using a mesh and taking its value at the center of each cell.
The experimental environment has been the following:(i)CPU Intel(R) Core(TM) i7-4500U CPU 1.8 GHz.(ii)Running Software Microsoft Visual C++ 2012.
The test 3D images exhibit a complex structure containing several three-dimensional objects. In complex problems, it may be difficult to use deformable model algorithms. Therefore, we have compared DA3DED with 3D difference filters (Sobel, Prewitt) which are suitable for arbitrary (unstructured) images. In the experiments with the 3D Sobel method, the following mask to obtain the partial derivative of the image intensity with respect to the variable was adopted:
The corresponding mask used for the 3D Prewitt method is
We have used the MATLAB notation for matrices. The masks for the derivatives with respect to and can be obtained from the above ones [10, 31].
The results obtained with the Prewitt method are similar to those obtained with the Sobel method, we have not included them to save space.
3D EDAS-1 consists of applying EDAS-1 to the straight lines, perpendicular to the plane , passing through each one of the points of . This process is performed for each coordinate plane.
The numerical results for 3D EDAS-1 and DA3DED are reported in Tables 1 and 2. We call the total number of intervals and the total number of intervals containing jumps, generated in the applications of EDAS-1. is a parameter of adaptive cubature; see [22]. The results in Table 2 have been obtained with the same values of , , than those in Table 1.
The synthetic models are shown in Figures 3, 6, 9, 12, and 15. A graphic comparison of the results obtained with 3D Sobel, 3D EDAS-1, and DA3DED is shown in Figures 4, 5, 7, 8, 10, 11, 13, 14, 16, and 17.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
4. Conclusions
Fast and accurate edge detection in 3D images is necessary in many applications. Some typical examples include image guided surgery, 3D assisted face recognition, safety and security applications [32], seismic imaging, and satellite images with natural color [33].
The problem is challenging because 3D images involve an enormous amount of data (voxels) that must be processed. Algorithms that provide a simple point based representation of the jump discontinuity set (3D filters) process all the voxels in an image.
However, in most images the complexity varies from a region to another. The proposed algorithm (DA3DED) is based on this fact. The depth of the treatment is proportional to the complexity of the region into consideration.
DA3DED has been tested on 3D images that modelize real problems (composites, fractures). It has been much faster than the 1D edge detection algorithm for 3D images derived from EDAS-1.
The doubly adaptive algorithm provides a thorough description of the edges in complex zones. Regions with low complexity are described in less detail, but in this case the discontinuity surface can be effectively reconstructed using techniques such as interpolation.
Competing Interests
The authors declare that they have no competing interests.