Edge Detection from MRI and DTI Images with an Anisotropic Vector Field Flow Using a Divergence Map

The aim of this work is the extraction of edges from Magnetic Resonance Imaging (MRI) and Diffusion Tensor Imaging (DTI) images by a deformable contour procedure, using an external force field derived from an anisotropic flow. Moreover, we introduce a divergence map in order to check the convergence of the process. As we know from vector calculus, divergence is a measure of the magnitude of a vector field convergence at a given point. Thus by means level curves of the divergence map, we have automatically selected an initial contour for the deformation process. If the initial curve includes the areas from which the vector field diverges, it will be able to push the curve towards the edges. Furthermore the divergence map highlights the presence of curves pointing to the most significant geometric parts of boundaries corresponding to high curvature values. In this way, the skeleton of the extracted object will be rather well defined and may subsequently be employed in shape analysis and morphological studies.


MRI and DTI Images
Boundary extraction and image segmentation are widely used in the field of Medical and Biological Image Analysis, allowing researchers to make significant progresses in approaches that utilize images at all spatial scales, ranging from cellular, molecular or organ imaging.In this work we present an implementation of a method for feature extraction in MRI and DTI datasets [1].In recent years, the OPEN ACCESS DTI technique has been greatly developed, given that it can be used as a marker for white matter (WM) tracts' integrity through the detection of the spontaneous diffuse on motion of water molecules in the brain.WM water diffusion appears to be greater along fiber directions than in perpendicular ones, therefore, the diffusion process is strongly directional and it can be used to estimate the patterns of WM connectivity in human brain.WM diffusion is a three-dimensional process strongly anisotropic, because axonal membranes and myelin sheaths represent barriers to random water molecular motion.Along fiber directions the translational movement is about three to six times faster than perpendicularly.Consequently, it would be of great help to indicate the overall orientation of fiber bundles.As a result of this decreasing water mobility, the diffusion coefficient is smaller when measured perpendicularly to prevalent fiber directions.Anisotropic diffusion is more adequately characterized by a tensor D T rather than a single scalar D. The diffusion tensor D T is a symmetric array of nine elements that describes mobility along orthogonal axes In this way, the effects of anisotropy could be fully investigated providing a great amount of details about cerebral microstructures.As it is very difficult to display tensor data, for a synthetic representation, we can recur to an ellipsoid centered at each voxel, whose elongated shape describes locally the asymmetrical level of diffusion.
The Cartesian equation of the ellipsoid, with the center in the origin of the reference system, is given by: In Diffusion Tensor Imaging (DTI), the matrix D T could be determined at any voxel by means of a best-fitting procedure, using at least six noncollinear magnetic gradient directions in order to evaluate the six independent coefficients of the equation.The ellipsoid is a three-dimensional representation of mean diffusion distances covered in space by molecules at a given time interval.Its major principal axis is oriented in the direction of maximum diffusivity.In the case of isotropic diffusion, the ellipsoid will be a sphere with a radius proportional to the diffusion coefficient D. The diagonalization of the symmetric diffusion tensor D T yields three eigenvectors   representing major, medium, and minor principal axes, whereas the corresponding eigenvalues   estimate apparent diffusivity along their directions.Geometrically the level of anisotropic diffusion can be expressed by how much ellipsoidal shape differs from that of a sphere (Figure 1).Mathematically it may be evaluated by the degree to which the three tensor eigenvalues differ from one another.A common measure of diffusion motion is fractional anisotropy FA [2]: a scalar ranging from 0 (isotropic diffusion) to 1 (maximum anisotropic level), invariant for affine transformations.Grayscale images, or FA maps, may be generated encoding values from the unitary interval [0,1] to the gray color space [0,255].Dark regions correspond to isotropic diffusion with a spherical shape of the ellipsoid whereas bright regions are anisotropic zones with an ellipsoid of more or less elongated shape.In order to represent information about the directionality of maximum diffusivity, the three components of 1 v  are encoded in corresponding levels of red-green-blue color.The resulting maps, called Directionally Encoded Color Images (DEC) or directional maps, resume the degree of anisotropy and reveal the local structure of fiber directions.In this study, we also carry out edge detection using structural T1-weighted images.In this way, we could investigate possible variations in volume or topography of gray matter (GM) between individuals from populations in different physiological or pathological conditions, through voxel-based morphometry (VBM) analysis [3,4].The processed GM and FA images derived from a dataset composed by patients affected by Alzheimer disease and subjects of a control group [5].

Generalized Gradient Vector Flow and Divergence Map
An active contour or snake is a curve defined within a given 2D image I(x, y) and subjected to modifications under the action of forces, until the evolving curve fits well into a final contour [6][7][8].In traditional parametric models, a snake is expressed explicitly by equations . The final shape of the contour to be extracted will be such as to minimize an energy functional that is the sum of an internal energy associated with it and an external energy related to a potential function, thus we have: the first term represents the internal energy that is in relation to the degree of flexibility of the active contour, so expressed: where α(s) is a function that controls the contour tension, while β(s) regularizes its rigidity.The external energy E x s is the energy associated with an external conservative force field, given by the gradient of a potential energy function   , deriving from the intensity image I(x, y), for example as follows: where  is the gradient operator.By using a variational approach [9], a contour that minimizes the total energy must satisfy the Euler-Lagrange equation: where are the external forces.From minimization of energy we derive a resolution of a static problem.By introducing a time-variable t, we may realize a deformable model able to create a geometrical shape that evolves over time.In this way, neglecting the inertial term and thus the second order derivatives, the active contour )) , ( where elasticity and rigidity are considered as constant functions and   s x 0  defines an initial curve. When the solution ) , ( t s x  is stabilizing, time derivatives will become null and we carry out the solution of Equation (2).In order to reduce the considerable sensitivity of this model to initial contours, edge detection may be performed using a different class of external forces, the GGVF force field or Generalized Gradient Vector Flow obtained by solving a diffusion problem [10][11][12].
In the GGVF framework the external force field  can be found as a solution of the following diffusion equation [11]: where 2   is the Laplacian operator, ) , ( y x f  is the gradient of the edge map f(x, y), derived from the gradient of the brightness function I(x, y), for example as or with any other edge detector.The field vectors point to the closest object boundaries with norms significantly different from zero in proximity of them.The functions   f g  and   f h  are non-negative and generally not uniform,   f g  is monotonically decreasing, since the vector field   will be weakly variable far from the edges where image intensities are expected to be rather uniform.On the other hand,   f h  should be monotonically increasing thus, close to boundaries, the vector field   should have a trend nearly equal to f  .In the generation of deformable contours, the main drawbacks are a weak convergence of models towards edges, specially in regions with highly variable concavities, the initialization problem, i.e., the excessive influence of shape and initial position of the active contour and the capture range, i.e., the size of area inside which an active contour can be initialized [13][14][15].We would like to introduce a method that tries to reduce the initialization problem because, as we can see for the test image of Figure 2, the convergence of arbitrary initial contours for the GGVF leads only partially to the expected results, giving seemingly incomprehensible outcomes.In this work, we suggest a more general analysis for the diffusion process of the external force field v  .To this end, we could see the GGVF Equation (4) as a special case of the following generalized parabolic equation [12,16]: where div is the divergence operator,   is a term generating an external force for the diffusion process of the field ) , , ( t y x v  , g(.) is the conduction coefficient that must be null or tend to zero at boundaries.If the diffusivity function g(.) is monotonically decreasing to zero, diffusion will be stopped across the edges to be extracted, so the vector flow can take place inside or outside the region.
We obtain the GGVF Equation ( 4) with initial condition   f , considering: where k is a constant positive value.If we have an anisotropic flow of the vector field v  with a null external force, i.e.,   0 , the Equation (5) will become: In this case we achieve results quite similar to those of the GGVF field, given that the term used in Equation ( 4), tends to zero either near edges, for the chosen initial conditions, or far from them, since the function   is not significant.Therefore, it is very reasonable to expect that the results of Equation ( 4) are very similar to those of Equation ( 5) with   0 . However, now considering the problem through the parabolic Equation ( 7), we could interpret the process as a field flow from boundaries toward the inside or outside without crossing edges and pointing to them because they are initially equal to the edge map gradient.As a consequence, the vector field can capture object boundaries from either side, as we could note in Figure 3.As we know from vector calculus, the divergence of a vector field is a measure of the field convergence at a given point by means signed scalar values; in other words, the divergence is the amount of field flux entering or leaving a point.Then, if we compute the divergence of ) , , ( t y x v  , we would obtain negative values in correspondence to object boundaries towards which the vector field converges (sinks), whereas positive values would outline regions from which the vector field spread out (sources), see Figure 3. Therefore, by evaluating the divergence of the force field v  , we could analyze the main features of its convergence [16].
From now on we name divergence map the grayscale image (Figure 4) reproducing divergence values of v  at different times t.Through careful analysis, we may point out that it is characterized by a gray background with divergence values near zero, dark curves with negative divergence in correspondence to edges towards which the vector field converges, and a system of light curves with positive values, defining regions from which the vector field comes out.Moreover, the sides of areas that delimit parts of image inside which the field is null (see enlarged image in the red box of Figure 3), collapse each other, especially in conjunction with long and deep concavities.Therefore, these curves, with positive divergence values, demarcate edge sections with high curvatures (Figure 4) and a great geometrical significance which concur in forming the skeleton of our figure.Upon variation of the vector field v  in time, there is a consequent variation of the related divergence map.As can be seen in Figure 5, where a circular initial contour is superimposed on the divergence map, the parts of the deformable curve that are positioned in areas from which the vector field diverges in the direction of boundaries can be pushed towards them.On the contrary, those traits that are inside the regions where the vector field is null remain trapped into their interior.Consequently, the divergence analysis of the external force field v  can be very useful in feature extraction, since it allows to delimit regions from which the field flow originated.Furthermore, we could note that the divergence map put in evidence the presence of curves pointing to the most significant geometric parts of boundaries with high curvature values.In this way, the geometrical shape of the extracted objects with their most significant characteristics will result well defined.After these considerations, in this paper, we suggest to use the map contour of divergence to select automatically an initial curve for the deformation process.To this end, given an image I(x, y), at first it has been evaluated v  and its divergence map in the gray color space C = [0,255], as follows: then we have generated a contour map of the field divergence I D , in order to visualize its level sets, i.e., the set of points in the image domain D where divergence is constant (Figure 6a).Subsequently, a level curve of high intensity, corresponding to high divergence values has been automatically selected as initial contour (Figure 6b).In this way, the chosen curve certainly encloses areas from which the field diverges, thereby pushing the deformable contour towards object edges, where sinks or zones of maximum field convergence are localized.As can be seen in Figure 6b, edges have been detected correctly.

Results and Discussion
If now we consider in Equation ( 5) an external force   to the steady state will result faster than in GGVF solution.To this end we have used the following anisotropic diffusion equation: with an external force for example as . In this way we will be able to speed up the convergence process, because   is approximately null near edges but increases as moving away from them.An efficient numerical scheme to approximate the equilibrium solution of the diffusion Equation ( 8) may be the one proposed by Perona-Malik [17] that recurs to the 4-neighbors discretization of the Laplacian operator: with the conduction function g(.) numerically evaluated as follows: through a gradient approximation of f  along horizontal and vertical directions, so we have The vector field so determined, from now on called AVF field, will be used for image processing in the subsequent applications.In Figure 7a,b, we could compare the AVF field derived from Equation ( 8) with GGVF field obtained through Equation ( 7), for the same number of iterations.This model could provide us an extension of the capture range, specifically an enlargement of areas inside which any initial contour may be placed.These regions result as well defined in the divergence map because it is delimited by curves corresponding to opposite values of divergence.As a consequence, we have a complete characterization of the capture range for a vector field that will be of great help in snake initialization procedures.These areas are varying according to the number of iterations after which the flow process has been interrupted in the numerical resolution of both Equations ( 7) and (8).Consequently, the evaluation of the AVF vector field will be less time consuming compared to the GGVF because we could choose a reduced number of iterations in order to have a comparable capture range size and to achieve similar results.Moreover, by increasing the number of iterations, the sides of areas that delimit parts of image inside which the field is null, collapse each other, especially in conjunction with deep and narrow concavities, as can be seen in Figure 8b.Therefore, the curves of high divergence allow us to clearly identify those parts of detected edges with high curvatures that result to be the most significant from the geometrical and morphological point of view.The first processed GM image (Figure 11) refers to a subject of the control group (see the first image of Figure 12).The force field has been evaluated from Equation ( 9) after eight iterations, the overlapping green curves represent initial contours that are automatically selected from the contour map of divergence.In Figures 13 and 14 we have treated GM images derived from two patients affected by Alzheimer's disease and they refer respectively to the second and third images which are shown in Figure 12.Once edges are detected, we produce a boundary representation of gray matter that can be used for an automatic analysis of shapes from the geometrical, metrical or morphological point of view.Moreover, we test the AVF method with T1-weighted MRI images of our dataset.The initial positions of the active contours are shown in green overlaid on the real image, whereas the outlines of the final contours are drawn in red.Many details of borders are captured both for gray matter and lateral ventricle profiles, whereas the interhemispheric fissure contours results are not completely retrieved, as we can see in Figure 15.For what concerns DTI images, this approach is addressed to demarcate the areas with higher diffusivity of water molecules.We could note the high level of noise present in them; in order to perform a noise reduction, images are initially pre-processed applying a median digital filter (Figure 16).As previously mentioned, the most significant information of a grayscale FA map are taken from zones with high intensity levels corresponding to preferential directions of water diffusion, presumably along axonal bundles.Thus after analyzing the intensity distribution of image pixels (Figures 17 and 18

Conclusions
In this paper we have proposed a novel approach for detecting edges from MRI and DTI images by a deformable contour procedure, using an external force field derived from an anisotropic flow.We have suggested a method to overcome some limitations of traditional parametric snakes, especially in regards to the initialization problem.We have derived an external force field for the deformation process of the active contour as a resolution of an anisotropic diffusion problem.The force field should be considered in relation to its divergence map because edge extraction will turn out well only if areas from which the vector flow comes out are completely included inside the chosen initial contour.As an automatic initialization procedure, we have proposed the selection of level curves by the divergence contour map.This technique does not require the user to input manually any initial curve close to the edges, so dealing with a large number of images may become less difficult and fastidious to do.In such a way, it is possible to obtain promising results able to explore the spatial distribution of white matter bundles for DTI images, whereas edge extraction from conventional MRI images allow us to thoroughly investigate the morphology of the human brain both for pathological and non-pathological subjects.
left, anterior-posterior, inferior-superior directions, respectively, as well as in coupled different space directions anterior, right-inferior, anterior-inferior, so we have the symmetric positive-definite matrix:

Figure 2 .
Figure 2. Edge extraction for a test image with arbitrary contours.

Figure 4 .
Figure 4. Test image and its divergence map.

Figure 5 .
Figure 5.An arbitrary contour and divergence map.

Figure 8 .
Figure 8.(a) Divergence map of AVF field; (b) Divergence map of AVF field after 60 iterations after 180 iterations.

Figure 9 .
Figure 9. FA and directional Maps derived by the Software Tools FSL.

Figure 10 .
Figure 10.Segmentation of a gray matter image (GM), derived by the Software Tools FSL.

Figure 11 .
Figure 11.Edge extraction with AVF force field after 8 iterations.

Figure 12 .
Figure 12.GM images for a subject of the control group and two patients affected by Alzheimer's disease.

Figure 13 .
Figure 13.Edge extraction with AVF force field after 8 iterations.

Figure 14 .
Figure 14.Edge extraction with AVF force field after 8 iterations.
), we have obtained an image segmentation by using thresholds ranging from 180 to 210.Subsequently, edge extraction has been performed through an automatic selection of an initial contour from the level curves of the divergence map (Figures18-20).

Figure 17 .
Figure 17.Intensity distribution of image pixels for a FA map.

Figure 18 .
Figure 18.DTI image extraction for a non-pathological subject.

Figure 19 .
Figure 19.DTI image extraction for a non-pathological subject.

Figure 20 .
Figure 20.DTI image extraction for a pathological subject.