Coherence Enhancement Based on Recursive Anisotropic Scale-Space with Adaptive Kernels

: The reduction of unnecessary details is important in a variety of imaging tasks. Image denoising can be generally formulated as a diffusion process that iteratively suppresses undesirable image features with high variance. We propose a recursive diffusion process that simultaneously computes the local geometrical property of the image features and determines the size and shape of the diffusion kernel, leading to an anisotropic scale-space. In the construction of the proposed anisotropic scale-space, image features due to undesirable noise are suppressed while signiﬁcant geometrical image features such as edges and corners are preserved across the scale-space. The diffusion kernels are iteratively determined based on the local geometrical properties of the image features. We demonstrate the effectiveness and robustness of the proposed algorithm in the detection of curvilinear features using simple yet effective synthetic and real images. The algorithm is quantitatively evaluated based on the identiﬁcation of ﬁssures in lung CT imagery. The experimental results indicate that the proposed algorithm can be used for the detection of linear or curvilinear structures in a variety of images ranging from satellite to medical images.


Introduction
The extraction of significant features from images plays an important role in a variety of applications ranging from low-level image processing to high-level computer vision techniques. The reconstruction of objects of interest is often required to interpret semantic structure from visual information. Image formation is, in general, a complex process involving the interplay between the physical environment and observation conditions. However, it is often necessary to discard nuisance factors such as undesirable intensity perturbations and geometrical distortions. This naturally leads to the analysis of imagery at multiple scales based on scale-space theory [1][2][3][4]. Scale-space theory provides an analytic framework in which a signal is represented at multiple scales constructed by a family of scale parameters. A typical scale-space representation is the Gaussian scale-space developed based on the solution of the heat equation [5][6][7][8]. This linear scale-space has numerous attractive properties, namely, scale-space axioms that include linearity, semi-group structure, existence of an infinitesimal generator, non-creation of local extrema, rotational symmetry, and scale invariance [9][10][11][12]. Several scale-space algorithms that employ wavelets [13][14][15] and a variety of kernels such as local phases forming monogenic signals [16] and Beltrami flows [17] were proposed. The analysis of visual information via scale-space is often beneficial because image features at a variety of scales coexist in a single image. It is desirable for features at a specific scale to be characterized at the corresponding scale. However, the feature scales are generally not known in advance for correspondences to be established between the scales and features at the characterization stage. This necessitates an alternative approach in which the analysis and feature characterization are performed in an alternative manner. One of the major drawbacks of linear scale-space approaches is the dislocation of significant image features such as edges and corners. This occurs because insignificant perturbations in appearance are suppressed while the scale increases. There is a trade-off between the degree of geometrical feature distortions and robustness against insignificant perturbations in scale-space construction [18]. It is often desirable to suppress nuisance factors such as noise while preserving geometrical properties that serve as characteristic indicators for recognizing objects of interest [19]. Thus, there have been several of proposals for overcoming the limitations arising from the isotropy of linear operations during the construction of scale-space by using anisotropic operations [20][21][22][23][24][25] in which the diffusion process is designed to consider local geometrical properties such as the direction and contrast of linear features. In addition to the linear features, curvilinear feature characterization based on structure tensors was developed; in this case, the orientations of the local features are used in an isotropic manner in the diffusion computation [26][27][28][29]. Curvilinear feature characterization was applied to the detection of curvilinear structures such as vessels [30][31][32][33], and for tubular structure segmentation [34] and fingerprint recognition [35][36][37]. Several approaches for anisotropic diffusion were developed based on structure-adaptive anisotropic filters [38,39] and anisotropic transformations [40]. However, static kernels have still been employed in conventional anisotropic diffusion to compute structure tensors for measuring local feature orientations. This motivates our proposed anisotropic diffusion process in scale-space based on varying the kernel shapes. In this paper, we propose an anisotropic scale-space based on data-driven kernels that adjust their sizes and shapes adaptively. The diffusion process is performed recursively based on the modified kernels, which are designed to consider the orientations of local geometrical features. The resultant desirable anisotropic properties improve the characterization of local features while suppressing insignificant details. In summary, our proposed algorithm is developed based on the scale-space theory; however our algorithm proceeds to apply anisotropic diffusion process based on the adaptive kernels that are alternatively designed to take into consideration of local geometrical features whereas the conventional anisotropic diffusion approaches are developed using static smoothing kernel for computing structure tensor in the analysis eigen-system.
We present preliminaries in Section 2 followed by the proposed algorithm in Section 3. Section 4 describes detection of curvilinear features based on the recursive anisotropic scale-space. The robustness and effectiveness of the proposed method are numerically demonstrated in Section 5. The conclusion follows in Section 6.

Preliminary
Let I : Ω → R be an image whose spatial domain is denoted by Ω. The conventional Gaussian scale-space is constructed by solving the heat equation defined by where t ≥ 0 is a scale parameter, u : Ω × R + → R represents the solution of the heat equation, and R + denotes the non-negative temporal domain. The initial condition u(x, 0) of the heat equation for constructing the Gaussian scale-space is defined by the given image I. Such a construction of the Gaussian scale-space based on the heat equation is often useful for characterizing the global properties of objects present in an image. However, the isotropic diffusion process in the construction of scale-space leads to the dislocation of critical points and undesirable geometrical distortion of the objects. Thus, algorithms that consider the local geometrical properties of image features in the determination of the degree and direction of the diffusion process were developed. A nonlinear diffusion algorithm was proposed to avoid the blurring and localization limitations that stem from the linear filtering [25]. An inhomogeneous process is applied to reduce the diffusivity at the locations with larger contrast. A typical measure of contrast is the norm of the image gradient ∇u 2 2 . The diffusion process is modified to consider the local image structure characterized by the image gradient in the following diffusion equation: where the diffusion is controlled by the function g : R + → R + defined by λ > 0 is a control parameter that determines the degree of diffusion. The diffusivity function g is designed to preserve significant geometrical image features such as edges and corners. It leads to an anisotropic behavior in which blurring occurs preferentially along linear features such as edges while blurring across edges is suppressed. The diffusivity function g determines the degree of diffusion based on the local image features using a scalar value, but it is desirable to also consider the orientation of the diffusion in addition to its degree. A more sophisticated diffusion process that considers a vector-valued structure tensor D to determine the direction of the diffusion process based on Gaussian smoothing is defined by Here, R is the rotation matrix defining the local coordinate system (λ 1 , λ 2 ) aligned with the gradient vector defined by x 2 ) denotes the coordinates in the image domain x ∈ Ω and I σ denotes the blurred image with a Gaussian kernel K σ defined by Here, K σ is a Gaussian kernel with the standard deviation σ, * denotes the convolution operator, and N is defined by The structure tensor matrix S 0 is defined by where ⊗ denotes the outer product operator and T denotes the transpose operator. The orthonormal basis of eigenvectors λ 1 and λ 2 has the geometrical properties λ 1 ∇I σ and λ 2 ⊥ ∇I σ . The structure tensor S ρ at scale ρ ≥ 0 can be obtained by smoothing S 0 (∇I σ ) with a Gaussian kernel K ρ as follows: where the Gaussian kernel K ρ in the spatial domain x ∈ Ω is defined by ρ denotes the standard deviation of the kernel K ρ . Please note that the diffusion process in the computation of the structure tensor S ρ (∇I σ ) at scale ρ is characterized by the isotropic Gaussian kernel K ρ , where the eigenvectors are computed from the blurred structure tensor based on the Gaussian kernel.

Recursive Anisotropic Diffusion
In this section, we propose a recursive anisotropic diffusion process based on an anisotropic Gaussian kernel that is determined using the orientation of the local features. Let G(x; σ 1 , σ 2 ) be a Gaussian kernel in the spatial domain x = (x 1 , x 2 ) ∈ Ω with the standard deviations σ 1 and σ 2 as follows: ∈ Ω, and σ 1 and σ 2 are the standard deviations in the x 1 and x 2 -directions, respectively. Let J (σ 1 ,σ 2 ) be an image blurred by an anisotropic Gaussian kernel G(x; σ 1 , σ 2 ) defined by We propose an anisotropic diffusion process based on the structure tensor computed from the blurred image J(x; σ 1 , σ 2 ) with an anisotropic Gaussian kernel G(x; σ 1 , σ 2 ) as follows: where t is a scale parameter based on the anisotropic Gaussian kernel G(σ 1 , σ 2 ). The shape of the Gaussian kernel G(x; σ 1 , σ 2 ) is recursively determined by the structure tensor computed from the blurred image J(x; σ 1 , σ 2 ) in such a way that the local geometrical features are recursively characterized. In the analysis of local image features, the structure tensor S t (x) at scale t is defined by where J t is the image diffused by the anisotropic Gaussian kernel G(σ 1 , σ 2 ). Eigen-analysis is performed for extracting local image features as follows: where v 1 and v 2 are the eigenvectors, and λ 1 and λ 2 are their corresponding eigenvalues. Then, the anisotropic smoothing kernel G(σ 1 , σ 2 ) is determined based on the eigenvalues λ 1 and λ 2 as follows: A local geometric feature obtained from the image is characterized by the eigenvalues λ 1 and λ 2 computed from the structure tensor S t at scale t as proposed in [29]. Small values of both λ 1 and λ 2 indicate a flat region and a large difference between λ 1 and λ 2 indicates the existence of a strong linear feature such as an edge. The diffusion process is performed along the direction of a larger eigenvector in proportion to the corresponding eigenvalue as follows: The structure tensor S t at scale t is computed based on the blurred image J t with an anisotropic kernel G(σ 1 , σ 2 ). The shape of the kernel is recursively determined based on the local features characterized by the eigenvalues of the structure tensor. Thus, the proposed algorithm is designed to improve the characterization power with regard to analyzing the local geometries of image features while being robust against insignificant intensity perturbations. This recursive anisotropic scale-space with improved geometrical features can be useful in further higher-level recognition tasks such as segmentation and recognition. The algorithm of the proposed recursive anisotropic diffusion process is presented as Algorithm 1.

Detection of curvilinear structures
In this section, we present an algorithm to detect curvilinear structures based on the representation obtained by the proposed recursive anisotropic diffusion process. Let u : Ω → R be a representation obtained from the proposed algorithm described in Algorithm 1, where the anisotropic diffusion process is recursively applied. The detection of curvilinear structures is often desired in a variety of imagery domains such as those involving medical images, satellite images, electronic designs, and maps. Features of interest forming curvilinear structures are better characterized when insignificant details are suppressed so that these curvilinear structures can be more easily detected. An algorithm for detecting curvilinear structures is developed based on the analysis of local image features [41][42][43], where the directional derivatives are used to determine the orientations of the curvilinear features as proposed in [4]. Let γ : C 1 → Ω be a parameterized curve from t ∈ [0, 1] ⊂ C 1 . The magnitude of the second derivative γ (t) at the point where γ (t) = 0 can be used as a saliency measure to characterize linear features. Positive contrast due to bright lines on a dark background is characterized by γ (t) 0, and negative contrast due to dark lines on a bright background is characterized by γ (t) 0. The characteristic kernel f (x) for linear features can be defined by where h > 0 is a function value and the kernel has the same shape as a Heaviside function with a threshold w > 0. It is often useful to relax the shape of the kernel for computational convenience and numerical stability by using the characteristic kernel relaxed with respect to linear features defined by which is a round bell-shaped function with a threshold w > 0. The convolution of the kernel function f with a Gaussian kernel is then used to compute the responses in the detection of linear features. Denoting the first and second derivatives in the x 1 and x 2 directions as d x 1 , d x 2 , d x 1 x 1 , d x 2 x 2 , respectively, we have where G x 1 , G x 2 and G x 1 x 1 , G x 2 x 2 denote the first and second derivatives of the Gaussian kernel G with respect to the x 1 and x 2 directions, respectively. The direction of the linear features present on the curve γ(t) can be obtained by the eigen-analysis of the following Hessian matrix: The normal direction to the curve γ(t) can be obtained from the eigenvectors of the Hessian matrix H.
The absolute values of the eigenvalues are used to detect linear features, which tend to have local eigenvalue maxima along the normal direction of the curve. The linear features can be detected by the saliency measure θ proposed in [41]. It is defined as It is likely that a stronger linear feature will have a larger saliency value. This can be used as an indicator in the detection of curvilinear features.

Experimental Results
In this section, we provide quantitative and qualitative results for curvilinear structure detection using satellite images and lung CT images. Please note that the number of data used for the evaluation of the proposed algorithm is limited due to the short of availability in manual annotations given by medical experts such as radiologists. Our proposed recursive anisotropic diffusion algorithm constructs representations that improve the linear features and extract the curvilinear structures. We qualitatively and quantitively compared the algorithm against isotropic diffusion using a Gaussian kernel, anisotropic diffusion based on the scalar diffusivity function [25], and anisotropic diffusion based on the vector-valued diffusivity function with enhanced coherence features [44]. Some of the example data are presented in Figure 1, which shows (a) an aerial image as well as (b) sagittal plane, (c) axial plane, and (d) coronal plane lung volume CT images from left to right. To demonstrate the scale-space developed by the proposed recursive anisotropic diffusion process, a series of diffused images at varying scales are presented in Figure 2, where (a) the original image, and the diffused images at scales (b) t = 10, (c) t = 30, (d) t = 50, and (e) t = 100 are shown from left to right. The qualitative visual comparisons of the diffusion process for the aerial image and for the sagittal, axial, and coronal plane lung CT images are presented from top to bottom in Figure 3. In this figure, (a) shows the original image, and (b)-(e) display the diffused images obtained via the isotropic scheme using a Gaussian kernel, the anisotropic scheme proposed by Perona-Malik [25], the coherence-enhancing scheme proposed by Weickert [44], and our proposed algorithm using recursive anisotropic diffusion, respectively, from left to right. As shown in the figures, a variety of linear features corresponding to roads and anatomical partitions are effectively characterized in the diffusion images obtained from our proposed algorithm.
We apply a curvilinear detection algorithm based on the saliency measure on the characteristic representations with respect to linear features obtained via isotropic diffusion using a Gaussian kernel, anisotropic diffusion using a scalar diffusivity function, anisotropic diffusion using a vector-valued diffusivity function, and our recursive anisotropic diffusion algorithm. The detection of curvilinear features is visually presented in Figure 4, where the detection results are shown on the bottom row, the ground truth on the middle row, and the obtained diffused images on the top row. Figure 4 shows the diffusion images that were obtained via (a) isotropic diffusion, (b) the anisotropic diffusion method developed by Perona-Malik, (c) the coherence-enhancing method developed by Weickert, and (d) our proposed diffusion algorithm, from left to right. The ground truth for the detection of fissures in lung CT images was obtained by the radiologist with manual annotations. For the detection of roads in the satellite images, the ground truth was obtained by human annotations based on the assumption that the central lines between relatively large blocks are likely to be roads. As shown in the detection results, our proposed algorithm has better accuracy and fewer false positives.   The results of fissure detection in the sagittal plane lung CT image are presented in Figure 5, where the diffused images are shown in the top row, the ground truth for the fissure to be detected is shown in the middle row, and the detection results for the fissure are shown in the bottom row. The diffused images shown in Figure 5 were obtained via (a) isotropic diffusion, (b) anisotropic diffusion method proposed by Perona-Malik, (c) coherence-enhancing diffusion method proposed by Weickert, and (d) our proposed algorithm, from left to right. As shown in the detection results, our algorithm outperforms the other algorithms in terms of both accuracy and sensitivity.
The result of the detection of fissures in the coronal and axial plane lung CT images is visually presented in Figure 6, where the ground truths and detection results are shown together.
The quantitative ROC analysis results for the detection of (a) roads in the aerial image and fissures in the (b) coronal, (c) axial, and (d) sagittal plane lung CT images in Figure 7 show that our algorithm outperforms the others. In the analysis of ROC curves, false positive indicates the detection of pixels that should not be categorized into any curvilinear feature and false negative indicates the miss of pixels that should be detected as curvilinear features. The precision and recall are computed for the detection of curvilinear features corresponding to roads in aerial images and fissures in coronal, axial, and sagittal plane CT images and presented in Table 1. Our algorithm consistently exhibits better results for both the precision and recall across different detection thresholds θ. The F-1 score is also computed for the detection of curvilinear features in Table 2 where the higher value indicate better performance. the accuracy was obtained based on the detection scheme that is determined by the overlap within a circular neighbourhood with respect to each location under determination. For each point x ∈ Ω, its associated circular neighbourhood is defined by B x = {y| x − y ≤ 5}. Thus, the detection is determined as true when it has an overlap with respect to the circular neighbourhood of the reference.

Conclusions
We presented an algorithm for detecting curvilinear features present in images. Linear features are detected in the image representations obtained by a recursive anisotropic diffusion process that was developed to improve the characterization power of linear features while suppressing insignificant details. The results demonstrate that our algorithm is a useful characterization scheme, especially for the detection of curvilinear features. The local geometrical image features are analyzed using the eigenvectors and corresponding eigenvalues of the structure tensor obtained in an anisotropic scale-space. The image features due to undesirable noise are suppressed while significant geometrical image features such as edges and corners are preserved across the scale-space in the proposed anisotropic scale-space. The structure tensor, which takes local geometrical features into account, is further enhanced for the analysis of local image structures by using anisotropic diffusion kernels that are recursively determined based on the eigenvalues of the structure tensor. The effectiveness and robustness of the recursive anisotropic diffusion process in the detection of curvilinear features was demonstrated. The quantitative and qualitative results obtained using satellite and medical images indicate that the proposed algorithm outperforms conventional anisotropic diffusion algorithms and provides better accuracy.