# A Doubly Adaptive Algorithm for Edge Detection in 3D Images

Authors: Sagrario Lantarón; M. Dolores López; Javier Rodrigo
Journal: Mathematical Problems in Engineering
 (2016)

---

## 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.


---

## Body

## 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 R3 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 planex1x2 containing the projection of the image and a set of points in this region. Then, it considers straight lines perpendicular to the plane x1x2 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 x1x3 and x2x3 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. Section2 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

Letf be a general piecewise continuous function in R⊂Rd, a compact d-interval.We assume thatR=⋃i=1nCi, where Cii=1n is a series of connected sets with pairwise disjoint nonempty interiors and we call Γi be the boundary of Ci,  i=1,…,n.DefineΓ≡⋃i=1nΓi. We call ΓJ the set of points in Γ with jump discontinuities. Achieving a good approximation of ΓJ is our aim.Letv0,v1,…,vd be a set of d+1 affinely independent vectors in Rd and T thed-simplex generated by them. Consider a function f:T⊂Rd→R. We call LTf, the unique affine map such that(1)LTfvj=fvj,j=0,1,…,d.It is proved at [21, 22] that the behavior of the next average integral AIT(f) depends on the continuity or lack of continuity of f on T:(2)AITf≡∫Tfx-LTfxdxvT,where vT is the Lebesgue measure of T.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-d for functions defined on subsets of Rd. In this section, we give an outline of EDAS-d. Its validity for d=1, 2 has been established in [22]. The case d=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-d algorithm are summarized.Step 1.
Consider the initial partition(j=1), Pj≡Tii=1n of R=[a,b], where a,b∈Rd.Step 2.
Put into the setGj the good simplices of Pj and put into the set Bj the bad simplices of Pj.
T is a good simplex if AITf≤E1orT≤E2 and T≤E4, where T is the diameter of T, E1 is the maximum local error of the approximant L, E2 is the approximation error of the points in ΓJ, and E4 is a positive real parameter.Step 3.
Divide each bad simplex into two simplices, by splitting its largest edge. We have a new partitionj=j+1.Step 4.
Repeat Steps 2 and 3 untilBj=∅.Step 5.
G = G j is the searched partition. Choose the following subset of G:(3)AΓE3J≡T∈G:AITf>E1,ja≡maxvi∈VT⁡fvi-minvi∈VT⁡fvi>E3,where E3 is the minimum magnitude of jump reported and E2 is also a stopping criterion.In the images studied in Section3, a voxel is considered as an edge voxel if it contains some barycenter of the simplices in AΓE3J.In practice, we have implemented the above algorithm using a binary tree whose leaves correspond to good simplices. In the cased=1, the integrals ∫T‍f(x)-LTf(x)dx 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 inR3 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.Figure 1
Continuous representation (a) and discrete representation (b) of a 2D image.
(a) (b)Letn and m be positive integers with n<m; define(4)n,m≡n,n+1,…,mA 3D image is given by an n1×n2×n3 array I=Iijk with indices i,j,k∈1,n1×1,n2×1,n3. From a 3D image I, one can build a piecewise constant function I~, defined on [0.5,n1+0.5]×[0.5,n2+0.5]×[0.5,n3+0.5] in the following way:(5)I~x1,x2,x3≡I111,ifx1,x2,x3∈0.5,1.5×0.5,1.5×0.5,1.5,Ii11,ifx1,x2,x3∈i-0.5,i+0.5×0.5,1.5×0.5,1.5i>1,I1j1,ifx1,x2,x3∈0.5,1.5×j-0.5,j+0.5×0.5,1.5j>1,Iij1,ifx1,x2,x3∈i-0.5,i+0.5×j-0.5,j+0.5×0.5,1.5i,j>1,I11k,ifx1,x2,x3∈0.5,1.5×0.5,1.5×k-0.5,k+0.5k>1,Ii1k,ifx1,x2,x3∈i-0.5,i+0.5×0.5,1.5×k-0.5,k+0.5i,k>1,I1jk,ifx1,x2,x3∈0.5,1.5×j-0.5,j+0.5×k-0.5,k+0.5j,k>1,Iijk,ifx1,x2,x3∈i-0.5,i+0.5×j-0.5,j+0.5×k-0.5,k+0.5i,j,k>1.


#### 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 functionf:a,b3→R. Let α,β⊂a,b and S=α,β2. Let x1,x2∈S be a fixed point in S. Call ex1,x2 the number of edge points of the function of one variable fx1,x2,x when x∈a,b. The projective complexity of f on the set K=S×a,b is defined as(6)PCf,K≡∫Sex1,x2dx1dx2vS.Letd be a positive integer. To approximate the above quantity, we subdivide S into d2 identical squares. Denote the center of the squares by (x1i,x2j),  1≤i,j≤d. Then, (6) can be approximated by(7)PCaf,K=∑i,j=1dex1i,x2jβ-α/d2β-α2=∑i,j=1dex1i,x2jd2.Observe thatPCaf,K does not depend on the size of the square S. We can consider e as a vector with d2 components and write (7) as(8)PCaf,K=e1d2.In each experiment, the value of d has been almost constant. Then, using the norm equivalence, we can define the following measure of complexity:(9)PC∗f,K=e∞.Although we have considered sets in the plane x1x2, similar definitions can be stated for sets in the planes x1x3 or x2x3. 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 ann1×n2×n3 array I=Iijk with indices (i,j,k)∈1,n1×1,n2×1,n3.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 rectangleR=n,m×p,q with n<m and p<q. The subdivision is performed computing r=n+m-n/2 and s=p+q-p/2 (x denotes the floor function that returns the greatest integer less than or equal to x). The resulting discrete rectangles are(10)R1=n,r×p,s,R2=n,r×s+1,q,R3=r+1,m×s+1,q,R4=r+1,m×p,s.Suppose now that we want to select d+12 points from R, distributed regularly. The selected points are obtained as the set(11)N=h0,h1,…,hd×v0,v1,…,vd,where the vectors h and v are given by the following algorithm:left=n-0.5;right=m+0.5;up=q+0.5;down=p-0.5;h[0]=n;h[d]=m;v[0]=p;v[d]=q;for(i=1;i<d;i++)ht=left+i(right-left)/d;vt=down+i(up-down)/d;hi=ht+0.5;  if  ht-ht==0.5h[i]=h[i]-1;v[i]=⌊vt+0.5⌋;  if  vt-vt==0.5v[i]=v[i]-1;Observe that x+0.5 returns the integer nearest to x. The above algorithm obtains a subset of [[n,m]]([[p,q]]) whose points are distanced approximately by n-m/dq-p/d. If d>n-m or d>q-p, the subdivision is not possible.


#### 2.2.5. Split Criterion

The set of indices of a 3D image is given by(12)M=1,n1×1,n2×1,n3.The projections of this set on the coordinate planes are(13)Ra=1,n1×1,n2,Rb=1,n1×1,n3,Rc=1,n2×1,n3.Consider a setRα,  α=a,b,c. Suppose that R=n,m×p,q is a subset of Rα. The split criterion for R is defined by the following procedure.Split CriterionStep 1.
Ifd>n-m∥d>q-p return 0.Step 2.
Selectd+12 points from R, distributed regularly (Section 2.2.4). The result is the set N.Step 3.
Consider the straight lines perpendicular to the plane containingRα, 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 pointk∈NNEPk; see Figure 2. Compute the minimum distance between these edge points MINDk. In this step, all the edge points found are stored in a file.Figure 2
EDAS-1 obtains edge points (red) lying on the straight lines perpendicular to the considered plane.Step 5.
Compute(14)MAXNEP=maxk∈N⁡NEPk,MIND=mink∈N⁡MINDk.Step 6.
If we callTNEP the limit of the maximum number of edge points and TMIND the limit of the minimum distance between edge points, thenif(MAXNEP≥TNEP∥MIND≤TMIND) return 1; else return 0.A discrete rectangle is subdivided when the measure of projective complexity corresponding to the discrete rectanglePC∗ 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. SPCR takes the value 1 if the split criterion is satisfied and 0 in the contrary case.


### 2.3. Doubly Adaptive Algorithm

Consider the discrete rectangleRα=1,ni×1,nj. Let R be a discrete rectangle contained in Rα. If SPCR=0, we call R a good discrete rectangle; otherwise, R is called bad.Doubly Adaptive Algorithm for 3D Edge Detection (DA3DED)Step 1.
ConsiderRα=1,n1×1,n2. 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 G1 and the bad rectangles are put into the set B1.Step 2.
At each stepj we have a set of good discrete rectangles Gj and a set of bad discrete rectangles Bj. 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 Gj+1 and Bj+1.Step 3.
The algorithm stops ifBj=∅. If the stopping criterion is not satisfied, go to Step 2.Step 4.
Repeat the above steps forRb=1,n1×1,n3 and Rc=1,n2×1,n3.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 rectangleR):struct nodeintn; int m; int p; int q;intspc;intpt;intch[4];intd;Q,where, n, m, p, and q are the coordinates of the corners of the discrete rectangle:spc: indicator variable that takes the value 1 if the algorithm has applied thesplit criterion to the rectangle and 0 otherwise.pt: number of the node parent of the current one.ch[j]: jth child of the current node (−1 in case of leaf nodes). ch[0] is the left child and ch[3] is the right child.d: determines the subset of R used by the split criterion.The algorithm uses alternately two distinct values ofd: d1 and d2, in order to avoid the repetition of computations. In practice, these values are almost equal. We denote by Rc the discrete rectangle associated with the current node c.Algorithm to Generate the Quadtree. See Algorithm 1.Algorithm 1
(i) Input of data:nk, nl, 1≤k<l≤3, E1, E2, TNEP, TDMIN, d1, d2, po=0. (ii) Initialize the node 0: Q0.n=1, Q0.m=nk;  Q0.p=1; Q0.q=nl; Q0.spc=1; Q0.pt=-1; (iii) Split the node 0 into four discrete rectangles and initialize the nodes 1, 2, 3 and 4. Qj.n, Qj.m, Qj.p, Qj.q, j=1,…,4, are obtained according to the method described in Section 2.2.4. for(j=1;j≤4;j++)  Qj.spc=0;Qj.pt=0;Qj.d=d1; (iv) Complete the members of the node 0, for(j=0;j≤3;j++)Q0.ch[j]=j+1; (v) Initialize the current node and the counter of nodesc=1; nn=4; (vi) while (po<4)  rch=Qc.ch[3]; (right child of the current node) if(Q[c].spc==0) (the split criterion has not been applied to the current node)  ifSPCRc==1 (the current node satisfies the split criterion)  SplitRc into four discrete rectangles, computeQ[nn+j].n, Q[nn+j].m, Q[nn+j].p, Q[nn+j].q, j=1,…,4, by the method described in Section 2.2.4. if(Q[c].d==d1)DV=d2; elseDV=d1; for (j=1;j≤4;j++) {Q[nn+j].spc=0;Q[nn+j].pt=c; Qnn+j.d=DV;Q[c].ch[j-1]=nn+j;} Q[c].spc=1; c=nn+1; (new current node) nn=nn+4;  else (the current node does not satisfy the split criterion)  Qc.spc=1; for(j=0;j≤3;j++)Q[c].ch[j]=-1; c=Q[c].pt; if(c==0)po=po+1;   else (the split criterion has been applied to the current node)  if(Q[rch].spc==1) (all the childs have been visited. Go to the parent)  c=Q[c].pt; if(c==0)po=po+1;  else (there exist childs not visited)  for(j=1;j≤3;j++) if (Q[Q[c].ch[j]].spc==0)break; c=Q[c].ch[j]; if(c==0)po=po+1;


## 2.1. Mathematical Preliminaries

Letf be a general piecewise continuous function in R⊂Rd, a compact d-interval.We assume thatR=⋃i=1nCi, where Cii=1n is a series of connected sets with pairwise disjoint nonempty interiors and we call Γi be the boundary of Ci,  i=1,…,n.DefineΓ≡⋃i=1nΓi. We call ΓJ the set of points in Γ with jump discontinuities. Achieving a good approximation of ΓJ is our aim.Letv0,v1,…,vd be a set of d+1 affinely independent vectors in Rd and T thed-simplex generated by them. Consider a function f:T⊂Rd→R. We call LTf, the unique affine map such that(1)LTfvj=fvj,j=0,1,…,d.It is proved at [21, 22] that the behavior of the next average integral AIT(f) depends on the continuity or lack of continuity of f on T:(2)AITf≡∫Tfx-LTfxdxvT,where vT is the Lebesgue measure of T.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-d for functions defined on subsets of Rd. In this section, we give an outline of EDAS-d. Its validity for d=1, 2 has been established in [22]. The case d=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-d algorithm are summarized.Step 1.
Consider the initial partition(j=1), Pj≡Tii=1n of R=[a,b], where a,b∈Rd.Step 2.
Put into the setGj the good simplices of Pj and put into the set Bj the bad simplices of Pj.
T is a good simplex if AITf≤E1orT≤E2 and T≤E4, where T is the diameter of T, E1 is the maximum local error of the approximant L, E2 is the approximation error of the points in ΓJ, and E4 is a positive real parameter.Step 3.
Divide each bad simplex into two simplices, by splitting its largest edge. We have a new partitionj=j+1.Step 4.
Repeat Steps 2 and 3 untilBj=∅.Step 5.
G = G j is the searched partition. Choose the following subset of G:(3)AΓE3J≡T∈G:AITf>E1,ja≡maxvi∈VT⁡fvi-minvi∈VT⁡fvi>E3,where E3 is the minimum magnitude of jump reported and E2 is also a stopping criterion.In the images studied in Section3, a voxel is considered as an edge voxel if it contains some barycenter of the simplices in AΓE3J.In practice, we have implemented the above algorithm using a binary tree whose leaves correspond to good simplices. In the cased=1, the integrals ∫T‍f(x)-LTf(x)dx 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 inR3 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.Figure 1
Continuous representation (a) and discrete representation (b) of a 2D image.
(a) (b)Letn and m be positive integers with n<m; define(4)n,m≡n,n+1,…,mA 3D image is given by an n1×n2×n3 array I=Iijk with indices i,j,k∈1,n1×1,n2×1,n3. From a 3D image I, one can build a piecewise constant function I~, defined on [0.5,n1+0.5]×[0.5,n2+0.5]×[0.5,n3+0.5] in the following way:(5)I~x1,x2,x3≡I111,ifx1,x2,x3∈0.5,1.5×0.5,1.5×0.5,1.5,Ii11,ifx1,x2,x3∈i-0.5,i+0.5×0.5,1.5×0.5,1.5i>1,I1j1,ifx1,x2,x3∈0.5,1.5×j-0.5,j+0.5×0.5,1.5j>1,Iij1,ifx1,x2,x3∈i-0.5,i+0.5×j-0.5,j+0.5×0.5,1.5i,j>1,I11k,ifx1,x2,x3∈0.5,1.5×0.5,1.5×k-0.5,k+0.5k>1,Ii1k,ifx1,x2,x3∈i-0.5,i+0.5×0.5,1.5×k-0.5,k+0.5i,k>1,I1jk,ifx1,x2,x3∈0.5,1.5×j-0.5,j+0.5×k-0.5,k+0.5j,k>1,Iijk,ifx1,x2,x3∈i-0.5,i+0.5×j-0.5,j+0.5×k-0.5,k+0.5i,j,k>1.


### 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 functionf:a,b3→R. Let α,β⊂a,b and S=α,β2. Let x1,x2∈S be a fixed point in S. Call ex1,x2 the number of edge points of the function of one variable fx1,x2,x when x∈a,b. The projective complexity of f on the set K=S×a,b is defined as(6)PCf,K≡∫Sex1,x2dx1dx2vS.Letd be a positive integer. To approximate the above quantity, we subdivide S into d2 identical squares. Denote the center of the squares by (x1i,x2j),  1≤i,j≤d. Then, (6) can be approximated by(7)PCaf,K=∑i,j=1dex1i,x2jβ-α/d2β-α2=∑i,j=1dex1i,x2jd2.Observe thatPCaf,K does not depend on the size of the square S. We can consider e as a vector with d2 components and write (7) as(8)PCaf,K=e1d2.In each experiment, the value of d has been almost constant. Then, using the norm equivalence, we can define the following measure of complexity:(9)PC∗f,K=e∞.Although we have considered sets in the plane x1x2, similar definitions can be stated for sets in the planes x1x3 or x2x3. 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 ann1×n2×n3 array I=Iijk with indices (i,j,k)∈1,n1×1,n2×1,n3.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 rectangleR=n,m×p,q with n<m and p<q. The subdivision is performed computing r=n+m-n/2 and s=p+q-p/2 (x denotes the floor function that returns the greatest integer less than or equal to x). The resulting discrete rectangles are(10)R1=n,r×p,s,R2=n,r×s+1,q,R3=r+1,m×s+1,q,R4=r+1,m×p,s.Suppose now that we want to select d+12 points from R, distributed regularly. The selected points are obtained as the set(11)N=h0,h1,…,hd×v0,v1,…,vd,where the vectors h and v are given by the following algorithm:left=n-0.5;right=m+0.5;up=q+0.5;down=p-0.5;h[0]=n;h[d]=m;v[0]=p;v[d]=q;for(i=1;i<d;i++)ht=left+i(right-left)/d;vt=down+i(up-down)/d;hi=ht+0.5;  if  ht-ht==0.5h[i]=h[i]-1;v[i]=⌊vt+0.5⌋;  if  vt-vt==0.5v[i]=v[i]-1;Observe that x+0.5 returns the integer nearest to x. The above algorithm obtains a subset of [[n,m]]([[p,q]]) whose points are distanced approximately by n-m/dq-p/d. If d>n-m or d>q-p, the subdivision is not possible.


### 2.2.5. Split Criterion

The set of indices of a 3D image is given by(12)M=1,n1×1,n2×1,n3.The projections of this set on the coordinate planes are(13)Ra=1,n1×1,n2,Rb=1,n1×1,n3,Rc=1,n2×1,n3.Consider a setRα,  α=a,b,c. Suppose that R=n,m×p,q is a subset of Rα. The split criterion for R is defined by the following procedure.Split CriterionStep 1.
Ifd>n-m∥d>q-p return 0.Step 2.
Selectd+12 points from R, distributed regularly (Section 2.2.4). The result is the set N.Step 3.
Consider the straight lines perpendicular to the plane containingRα, 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 pointk∈NNEPk; see Figure 2. Compute the minimum distance between these edge points MINDk. In this step, all the edge points found are stored in a file.Figure 2
EDAS-1 obtains edge points (red) lying on the straight lines perpendicular to the considered plane.Step 5.
Compute(14)MAXNEP=maxk∈N⁡NEPk,MIND=mink∈N⁡MINDk.Step 6.
If we callTNEP the limit of the maximum number of edge points and TMIND the limit of the minimum distance between edge points, thenif(MAXNEP≥TNEP∥MIND≤TMIND) return 1; else return 0.A discrete rectangle is subdivided when the measure of projective complexity corresponding to the discrete rectanglePC∗ 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. SPCR takes the value 1 if the split criterion is satisfied and 0 in the contrary case.


## 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-d for functions defined on subsets of Rd. In this section, we give an outline of EDAS-d. Its validity for d=1, 2 has been established in [22]. The case d=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-d algorithm are summarized.Step 1.
Consider the initial partition(j=1), Pj≡Tii=1n of R=[a,b], where a,b∈Rd.Step 2.
Put into the setGj the good simplices of Pj and put into the set Bj the bad simplices of Pj.
T is a good simplex if AITf≤E1orT≤E2 and T≤E4, where T is the diameter of T, E1 is the maximum local error of the approximant L, E2 is the approximation error of the points in ΓJ, and E4 is a positive real parameter.Step 3.
Divide each bad simplex into two simplices, by splitting its largest edge. We have a new partitionj=j+1.Step 4.
Repeat Steps 2 and 3 untilBj=∅.Step 5.
G = G j is the searched partition. Choose the following subset of G:(3)AΓE3J≡T∈G:AITf>E1,ja≡maxvi∈VT⁡fvi-minvi∈VT⁡fvi>E3,where E3 is the minimum magnitude of jump reported and E2 is also a stopping criterion.In the images studied in Section3, a voxel is considered as an edge voxel if it contains some barycenter of the simplices in AΓE3J.In practice, we have implemented the above algorithm using a binary tree whose leaves correspond to good simplices. In the cased=1, the integrals ∫T‍f(x)-LTf(x)dx 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 inR3 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.Figure 1
Continuous representation (a) and discrete representation (b) of a 2D image.
(a) (b)Letn and m be positive integers with n<m; define(4)n,m≡n,n+1,…,mA 3D image is given by an n1×n2×n3 array I=Iijk with indices i,j,k∈1,n1×1,n2×1,n3. From a 3D image I, one can build a piecewise constant function I~, defined on [0.5,n1+0.5]×[0.5,n2+0.5]×[0.5,n3+0.5] in the following way:(5)I~x1,x2,x3≡I111,ifx1,x2,x3∈0.5,1.5×0.5,1.5×0.5,1.5,Ii11,ifx1,x2,x3∈i-0.5,i+0.5×0.5,1.5×0.5,1.5i>1,I1j1,ifx1,x2,x3∈0.5,1.5×j-0.5,j+0.5×0.5,1.5j>1,Iij1,ifx1,x2,x3∈i-0.5,i+0.5×j-0.5,j+0.5×0.5,1.5i,j>1,I11k,ifx1,x2,x3∈0.5,1.5×0.5,1.5×k-0.5,k+0.5k>1,Ii1k,ifx1,x2,x3∈i-0.5,i+0.5×0.5,1.5×k-0.5,k+0.5i,k>1,I1jk,ifx1,x2,x3∈0.5,1.5×j-0.5,j+0.5×k-0.5,k+0.5j,k>1,Iijk,ifx1,x2,x3∈i-0.5,i+0.5×j-0.5,j+0.5×k-0.5,k+0.5i,j,k>1.


## 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 functionf:a,b3→R. Let α,β⊂a,b and S=α,β2. Let x1,x2∈S be a fixed point in S. Call ex1,x2 the number of edge points of the function of one variable fx1,x2,x when x∈a,b. The projective complexity of f on the set K=S×a,b is defined as(6)PCf,K≡∫Sex1,x2dx1dx2vS.Letd be a positive integer. To approximate the above quantity, we subdivide S into d2 identical squares. Denote the center of the squares by (x1i,x2j),  1≤i,j≤d. Then, (6) can be approximated by(7)PCaf,K=∑i,j=1dex1i,x2jβ-α/d2β-α2=∑i,j=1dex1i,x2jd2.Observe thatPCaf,K does not depend on the size of the square S. We can consider e as a vector with d2 components and write (7) as(8)PCaf,K=e1d2.In each experiment, the value of d has been almost constant. Then, using the norm equivalence, we can define the following measure of complexity:(9)PC∗f,K=e∞.Although we have considered sets in the plane x1x2, similar definitions can be stated for sets in the planes x1x3 or x2x3. 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 ann1×n2×n3 array I=Iijk with indices (i,j,k)∈1,n1×1,n2×1,n3.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 rectangleR=n,m×p,q with n<m and p<q. The subdivision is performed computing r=n+m-n/2 and s=p+q-p/2 (x denotes the floor function that returns the greatest integer less than or equal to x). The resulting discrete rectangles are(10)R1=n,r×p,s,R2=n,r×s+1,q,R3=r+1,m×s+1,q,R4=r+1,m×p,s.Suppose now that we want to select d+12 points from R, distributed regularly. The selected points are obtained as the set(11)N=h0,h1,…,hd×v0,v1,…,vd,where the vectors h and v are given by the following algorithm:left=n-0.5;right=m+0.5;up=q+0.5;down=p-0.5;h[0]=n;h[d]=m;v[0]=p;v[d]=q;for(i=1;i<d;i++)ht=left+i(right-left)/d;vt=down+i(up-down)/d;hi=ht+0.5;  if  ht-ht==0.5h[i]=h[i]-1;v[i]=⌊vt+0.5⌋;  if  vt-vt==0.5v[i]=v[i]-1;Observe that x+0.5 returns the integer nearest to x. The above algorithm obtains a subset of [[n,m]]([[p,q]]) whose points are distanced approximately by n-m/dq-p/d. If d>n-m or d>q-p, the subdivision is not possible.


## 2.2.5. Split Criterion

The set of indices of a 3D image is given by(12)M=1,n1×1,n2×1,n3.The projections of this set on the coordinate planes are(13)Ra=1,n1×1,n2,Rb=1,n1×1,n3,Rc=1,n2×1,n3.Consider a setRα,  α=a,b,c. Suppose that R=n,m×p,q is a subset of Rα. The split criterion for R is defined by the following procedure.Split CriterionStep 1.
Ifd>n-m∥d>q-p return 0.Step 2.
Selectd+12 points from R, distributed regularly (Section 2.2.4). The result is the set N.Step 3.
Consider the straight lines perpendicular to the plane containingRα, 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 pointk∈NNEPk; see Figure 2. Compute the minimum distance between these edge points MINDk. In this step, all the edge points found are stored in a file.Figure 2
EDAS-1 obtains edge points (red) lying on the straight lines perpendicular to the considered plane.Step 5.
Compute(14)MAXNEP=maxk∈N⁡NEPk,MIND=mink∈N⁡MINDk.Step 6.
If we callTNEP the limit of the maximum number of edge points and TMIND the limit of the minimum distance between edge points, thenif(MAXNEP≥TNEP∥MIND≤TMIND) return 1; else return 0.A discrete rectangle is subdivided when the measure of projective complexity corresponding to the discrete rectanglePC∗ 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. SPCR takes the value 1 if the split criterion is satisfied and 0 in the contrary case.


## 2.3. Doubly Adaptive Algorithm

Consider the discrete rectangleRα=1,ni×1,nj. Let R be a discrete rectangle contained in Rα. If SPCR=0, we call R a good discrete rectangle; otherwise, R is called bad.Doubly Adaptive Algorithm for 3D Edge Detection (DA3DED)Step 1.
ConsiderRα=1,n1×1,n2. 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 G1 and the bad rectangles are put into the set B1.Step 2.
At each stepj we have a set of good discrete rectangles Gj and a set of bad discrete rectangles Bj. 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 Gj+1 and Bj+1.Step 3.
The algorithm stops ifBj=∅. If the stopping criterion is not satisfied, go to Step 2.Step 4.
Repeat the above steps forRb=1,n1×1,n3 and Rc=1,n2×1,n3.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 rectangleR):struct nodeintn; int m; int p; int q;intspc;intpt;intch[4];intd;Q,where, n, m, p, and q are the coordinates of the corners of the discrete rectangle:spc: indicator variable that takes the value 1 if the algorithm has applied thesplit criterion to the rectangle and 0 otherwise.pt: number of the node parent of the current one.ch[j]: jth child of the current node (−1 in case of leaf nodes). ch[0] is the left child and ch[3] is the right child.d: determines the subset of R used by the split criterion.The algorithm uses alternately two distinct values ofd: d1 and d2, in order to avoid the repetition of computations. In practice, these values are almost equal. We denote by Rc the discrete rectangle associated with the current node c.Algorithm to Generate the Quadtree. See Algorithm 1.Algorithm 1
(i) Input of data:nk, nl, 1≤k<l≤3, E1, E2, TNEP, TDMIN, d1, d2, po=0. (ii) Initialize the node 0: Q0.n=1, Q0.m=nk;  Q0.p=1; Q0.q=nl; Q0.spc=1; Q0.pt=-1; (iii) Split the node 0 into four discrete rectangles and initialize the nodes 1, 2, 3 and 4. Qj.n, Qj.m, Qj.p, Qj.q, j=1,…,4, are obtained according to the method described in Section 2.2.4. for(j=1;j≤4;j++)  Qj.spc=0;Qj.pt=0;Qj.d=d1; (iv) Complete the members of the node 0, for(j=0;j≤3;j++)Q0.ch[j]=j+1; (v) Initialize the current node and the counter of nodesc=1; nn=4; (vi) while (po<4)  rch=Qc.ch[3]; (right child of the current node) if(Q[c].spc==0) (the split criterion has not been applied to the current node)  ifSPCRc==1 (the current node satisfies the split criterion)  SplitRc into four discrete rectangles, computeQ[nn+j].n, Q[nn+j].m, Q[nn+j].p, Q[nn+j].q, j=1,…,4, by the method described in Section 2.2.4. if(Q[c].d==d1)DV=d2; elseDV=d1; for (j=1;j≤4;j++) {Q[nn+j].spc=0;Q[nn+j].pt=c; Qnn+j.d=DV;Q[c].ch[j-1]=nn+j;} Q[c].spc=1; c=nn+1; (new current node) nn=nn+4;  else (the current node does not satisfy the split criterion)  Qc.spc=1; for(j=0;j≤3;j++)Q[c].ch[j]=-1; c=Q[c].pt; if(c==0)po=po+1;   else (the split criterion has been applied to the current node)  if(Q[rch].spc==1) (all the childs have been visited. Go to the parent)  c=Q[c].pt; if(c==0)po=po+1;  else (there exist childs not visited)  for(j=1;j≤3;j++) if (Q[Q[c].ch[j]].spc==0)break; c=Q[c].ch[j]; if(c==0)po=po+1;


## 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 rectangleR):struct nodeintn; int m; int p; int q;intspc;intpt;intch[4];intd;Q,where, n, m, p, and q are the coordinates of the corners of the discrete rectangle:spc: indicator variable that takes the value 1 if the algorithm has applied thesplit criterion to the rectangle and 0 otherwise.pt: number of the node parent of the current one.ch[j]: jth child of the current node (−1 in case of leaf nodes). ch[0] is the left child and ch[3] is the right child.d: determines the subset of R used by the split criterion.The algorithm uses alternately two distinct values ofd: d1 and d2, in order to avoid the repetition of computations. In practice, these values are almost equal. We denote by Rc the discrete rectangle associated with the current node c.Algorithm to Generate the Quadtree. See Algorithm 1.Algorithm 1
(i) Input of data:nk, nl, 1≤k<l≤3, E1, E2, TNEP, TDMIN, d1, d2, po=0. (ii) Initialize the node 0: Q0.n=1, Q0.m=nk;  Q0.p=1; Q0.q=nl; Q0.spc=1; Q0.pt=-1; (iii) Split the node 0 into four discrete rectangles and initialize the nodes 1, 2, 3 and 4. Qj.n, Qj.m, Qj.p, Qj.q, j=1,…,4, are obtained according to the method described in Section 2.2.4. for(j=1;j≤4;j++)  Qj.spc=0;Qj.pt=0;Qj.d=d1; (iv) Complete the members of the node 0, for(j=0;j≤3;j++)Q0.ch[j]=j+1; (v) Initialize the current node and the counter of nodesc=1; nn=4; (vi) while (po<4)  rch=Qc.ch[3]; (right child of the current node) if(Q[c].spc==0) (the split criterion has not been applied to the current node)  ifSPCRc==1 (the current node satisfies the split criterion)  SplitRc into four discrete rectangles, computeQ[nn+j].n, Q[nn+j].m, Q[nn+j].p, Q[nn+j].q, j=1,…,4, by the method described in Section 2.2.4. if(Q[c].d==d1)DV=d2; elseDV=d1; for (j=1;j≤4;j++) {Q[nn+j].spc=0;Q[nn+j].pt=c; Qnn+j.d=DV;Q[c].ch[j-1]=nn+j;} Q[c].spc=1; c=nn+1; (new current node) nn=nn+4;  else (the current node does not satisfy the split criterion)  Qc.spc=1; for(j=0;j≤3;j++)Q[c].ch[j]=-1; c=Q[c].pt; if(c==0)po=po+1;   else (the split criterion has been applied to the current node)  if(Q[rch].spc==1) (all the childs have been visited. Go to the parent)  c=Q[c].pt; if(c==0)po=po+1;  else (there exist childs not visited)  for(j=1;j≤3;j++) if (Q[Q[c].ch[j]].spc==0)break; c=Q[c].ch[j]; if(c==0)po=po+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 modelsMaterials 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 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.(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 variablex1 was adopted:(15)MX:,:,1=-202;-303;-202,MX:,:,2=-303;-606;-303,MX:,:,3=-202;-303;-202.The corresponding mask used for the 3D Prewitt method is(16)MX:,:,1=-101;-101;-101,MX:,:,2=-101;-101;-101,MX:,:,3=-101;-101;-101.We have used the MATLAB notation for matrices. The masks for the derivatives with respect tox2 and x3 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 planexixj, passing through each one of the points of 1,ni×1,nj. This process is performed for each coordinate plane.The numerical results for 3D EDAS-1 and DA3DED are reported in Tables1 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. E5 is a parameter of adaptive cubature; see [22]. The results in Table 2 have been obtained with the same values of Ei,  i=1,…,5, than those in Table 1.Table 1
Performance of 3D EDAS-1 on the test models.
3D image E 1 E 2 E 3 E 4 E 5 TNI TNJI CPU (s) Spherical inclusions 1 0 - 1 1 0 - 2 1 0 - 2 10 10 4770465 93447 109.8 Icosahedral inclusions 1 0 - 1 1 0 - 2 1 0 - 2 10 10 6123908 229703 156.7 Cylindrical fibers 1 0 - 1 1 0 - 2 1 0 - 2 10 10 4854444 102482 111.8 Prismatic fibers 1 0 - 1 1 0 - 2 1 0 - 2 10 10 4940281 112624 115.1 Fracture model 1 0 - 1 1 0 - 2 1 0 - 2 10 10 8297400 452000 234.4Table 2
Performance of DA3DED on the test models.
3D image TMIND TNEP d 1 d 2 TNI TNJI CPU (s) Spherical inclusions 3 25 10 9 2137515 71586 51.3 Icosahedral inclusions 3 25 10 9 4339122 192821 110.4 Cylindrical fibers 3 25 10 8 2774380 89438 65.8 Prismatic fibers 4 25 6 5 3369707 108608 79.5 Fracture model 4 17 8 7 3983960 338064 124.9The synthetic models are shown in Figures3, 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.Figure 3
Spherical inclusions.Figure 4
Spherical inclusions (intersection with the planex1=50). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(a) (b) (c) (d)Figure 5
Spherical inclusions (intersection with the planex1=150). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(a) (b) (c) (d)Figure 6
Icosahedral inclusions.Figure 7
Icosahedral inclusions (intersection with the planex1=75). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(a) (b) (c) (d)Figure 8
Icosahedral inclusions (intersection with the planex1=150). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(a) (b) (c) (d)Figure 9
Cylindrical fibers (circular cross section).Figure 10
Cylindrical fibers (intersection with the planex1=50). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(a) (b) (c) (d)Figure 11
Cylindrical fibers (intersection with the planex1=110). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(a) (b) (c) (d)Figure 12
Prismatic fibers (triangular cross section).Figure 13
Prismatic fibers (intersection with the planex1=50). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(a) (b) (c) (d)Figure 14
Prismatic fibers (intersection with the planex1=150). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(a) (b) (c) (d)Figure 15
Fracture model.Figure 16
Fracture model (section perpendicular to thex1-axis). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(a) (b) (c) (d)Figure 17
Fracture model (intersection with the planex3 = 175). (a) Exact intersection, (b) 3D Sobel, (c) 3D EDAS-1, and (d) DA3DED.
(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.



---

*Source: 4904279-2016-12-01.xml*