A Doubly Adaptive Algorithm for Edge Detection in 3 D Images

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.


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][2][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][6][7].Among the direct methods we can find 1D edge detection methods (1D3DED) [8,9], spatial difference filters [10][11][12], polynomial fitting method [13], and deformable models [14][15][16][17][18][19][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 R 3 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  1  2 containing the projection of the image and a set of points in this region.Then, it considers straight lines perpendicular to the plane  1  2 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  1  3 and  2  3 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.

Mathematical Preliminaries.
Let  be a general piecewise continuous function in  ⊂ R  , a compact -interval.
We assume that  = ⋃  =1   , where {  }  =1 is a series of connected sets with pairwise disjoint nonempty interiors and we call Γ  be the boundary of   ,  = 1, . . ., .
Define Γ ≡ ⋃  =1 Γ  .We call Γ  the set of points in Γ with jump discontinuities.Achieving a good approximation of Γ  is our aim.
Let {k 0 , k 1 , . . ., k  } be a set of  + 1 affinely independent vectors in R  and  the d-simplex generated by them.Consider a function  :  ⊂ R  → R. We call   , the unique affine map such that It is proved at [21,22] that the behavior of the next average integral AI  () depends on the continuity or lack of continuity of  on : where V() 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].

Constituents of the Algorithm DA3DED
. DA3DED is based on several concepts and algorithms that are described in this subsection.

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 R  .In this section, we give an outline of EDAS-.Its validity for  =1, 2 has been established in [22].The case  = 3 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 Step 2. Put into the set G  the good simplices of   and put into the set B  the bad simplices of   .
is a good simplex if (AI  () ≤  1 or || ≤  2 ) and || ≤  4 , where || is the diameter of ,  1 is the maximum local error of the approximant ,  2 is the approximation error of the points in Γ  , and  4 is a positive real parameter.
Step 3. Divide each bad simplex into two simplices, by splitting its largest edge.We have a new partition  =  + 1.
Step 5.  = G  is the searched partition.Choose the following subset of : where  3 is the minimum magnitude of jump reported and  2 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 AΓ   3 .In practice, we have implemented the above algorithm using a binary tree whose leaves correspond to good simplices.In the case  = 1, 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].

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 R 3 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.
Let  and  be positive integers with  < ; define A 3D image is given by an  ( Let  be a positive integer.To approximate the above quantity, we subdivide  into  2 identical squares.Denote the Mathematical Problems in Engineering center of the squares by ( 1 ,  2 ), 1 ≤ ,  ≤ .Then, (6) can be approximated by Observe that PC  (, ) does not depend on the size of the square .We can consider  as a vector with  2 components and write (7) as In 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  1  2 , similar definitions can be stated for sets in the planes  1  3 or  2  3 .
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.

Subdivision and Selection
Procedures.Consider a 3D image given by an 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 R = [[, ]] × [[, ]] with  <  and  < .The subdivision is performed computing  =  + ⌊( − )/2⌋ and  =  + ⌊( − )/2⌋ (⌊⌋ denotes the floor function that returns the greatest integer less than or equal to ).The resulting discrete rectangles are Suppose now that we want to select ( + 1) 2 points from R, distributed regularly.The selected points are obtained as the set where the vectors ℎ and V are given by the following algorithm: If  >  −  or  >  − , the subdivision is not possible.

Split Criterion. The set of indices of a 3D image is given by
The projections of this set on the coordinate planes are Consider a set R  , ( = , , ).Suppose that R = [[, ]] × [[, ]] is a subset of R  .The split criterion for R is defined by the following procedure.
Step 3. Consider the straight lines perpendicular to the plane containing R  , passing through each point in N.
Step 4. Apply EDAS-1 to each one of these lines to obtain the number of edge points corresponding to each point k ∈ N(NEP(k)); see Figure 2. Compute the minimum distance between these edge points (MIND(k)).In this step, all the edge points found are stored in a file.Step 6.If we call TNEP the limit of the maximum number of edge points and TMIND the limit of the minimum distance between edge points, then if ((MAXNEP ≥ TNEP) ‖ (MIND ≤ TMIND)) return 1; else return 0.
A discrete rectangle is subdivided when the measure of projective complexity corresponding to the discrete rectangle (PC * ) 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 R, the split criterion is implemented with the function SPC.SPC(R) takes the value 1 if the split criterion is satisfied and 0 in the contrary case.

Doubly Adaptive Algorithm for 3D Edge Detection (DA3DED)
Step Divide R  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 G 1 and the bad rectangles are put into the set B 1 .
Step 2. At each step  we have a set of good discrete rectangles G  and a set of bad discrete rectangles B  .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 G +1 and B +1 .
Step 3. The algorithm stops if B  = 0.If the stopping criterion is not satisfied, go to Step 2.
Step 4. Repeat the above steps for DA3DED has been implemented using a quadtree.Below we detail this algorithm.

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 R): where, , , , and  are the coordinates of the corners of the discrete rectangle: spc: indicator variable that takes the value 1 if the algorithm has applied the split criterion to the rectangle and 0 otherwise.The algorithm uses alternately two distinct values of :  1 and  2 , 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.

Results and Discussion
In this section, we have tested DA3DED with several synthetic material models.The results have been compared with Mathematical Problems in Engineering 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 prisms The 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 200 × 200 × 200 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.
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  1 was adopted: We have used the MATLAB notation for matrices.The masks for the derivatives with respect to  2 and  3 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.The numerical results for 3D EDAS-1 and DA3DED are reported in Tables 1 and 2. We call TNI the total number of intervals and TNJI the total number of intervals containing jumps, generated in the applications of EDAS-1. 5 is a parameter of adaptive cubature; see [22].The results in Table 2 have been obtained with the same values of   ,  = 1, . . ., 5, than those in Table 1.
The synthetic models are shown in Figures 3, 6

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.

Figure 1 :
Figure 1: Continuous representation (a) and discrete representation (b) of a 2D image.

Figure 2 :
Figure 2: EDAS-1 obtains edge points (red) lying on the straight lines perpendicular to the considered plane.
Let [, ] ⊂ [, ] and  = [, ] 2 .Let ( 1 ,  2 ) ∈  be a fixed point in .Call ( 1 ,  2 ) the number of edge points of the function of one variable ( 1 ,  2 , ) when  ∈ [, ].The projective complexity of  on the set  =  × [, ] is defined as [26] 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.

Table 1 :
Performance of 3D EDAS-1 on the test models.

Table 2 :
Performance of DA3DED on the test models.3D EDAS-1 consists of applying EDAS-1 to the straight lines, perpendicular to the plane     , passing through each one of the points of [[1,   ]] × [[1,   ]].This process is performed for each coordinate plane.