K Delibasis

K Kottari

I Maglogiannis

Automated Detection of Streaks in Dermoscopy Images

Keywords: Computer vision, Skin lesions, Image Segmentation, Streaks detection, Dermoscopy images, Melanoma classification

In this paper we present a novel algorithm for the detection of dark linear structures, which appear in digital dermoscopy images of skin lesions and they are called as streaks in relevant literature. The proposed algorithm is capable of detecting such linear structures using local image curvature information obtained by the Hessian matrix. A linear structure is characterized as streak, based on its geometric characteristics. The streak detection algorithm is applied to a number of dermoscopy images containing malignant and non-malignant skin lesions. In this work we propose also a new class of streak based image features and we investigate their value in image classification according to malignancy.

Introduction

Skin cancer is considered among the most frequent types of cancer and one of the most malignant tumors. According to a recent study, its incidence has increased faster than that of almost all other cancers among young adults [START_REF] Reed | Increasing Incidence of Melanoma Among Young Adults: An Epidemiological Study in Olmsted County, Minnesota[END_REF]. Several additional reports [START_REF] Stern | Prevalence of a history of skin cancer in 2007: results of an incidence-based model[END_REF], [START_REF] Rogers | Incidence estimate of nonmelanoma skin cancer in the United States, 2006[END_REF], [START_REF]Cancer Facts & Figures[END_REF] have proven the importance of the early diagnosis in skin cancer related diseases and especially the melanoma cases. Thus Dermoscopy or Dermatoscopy, which is the examination of skin lesions or disorders through the use of the dermatoscope, a medical instrument capable of magnifying the examined area of the skin, has been widely used. This technique can reveal many color and structural features of the skin lesion, which would be impossible to observe with a naked eye. In addition advanced computer based systems can be used to capture a digital image produced by the dermatoscope, then analyze it by checking its morphological and color characteristics and extract measurable features. Such features could feed an expert system for providing an evaluation about the malignancy of the lesion [START_REF] Maglogiannis | Overview of advanced computer vision systems for skin lesions characterization[END_REF], [START_REF] Korotkov | Computerized analysis of pigmented skin lesions: A review[END_REF].

Methodologies and rules applied in the classic dermoscopy, can also be utilized in computer based dermoscopy. More specifically, in the conventional procedure, the following diagnostic methods are mainly used: (i) ABCD rule, (ii) Menzies rule, (iii) 7-point checklist and (iv) Pattern and Texture Analysis. For instance the ABCD rule investigates the asymmetry (A), border (B), color (C), and differential structures (D) such as Pigment Network and Streaks of the lesion and defines the basis for a diagnosis by a dermatologist. According to the published overviews [START_REF] Maglogiannis | Overview of advanced computer vision systems for skin lesions characterization[END_REF], [START_REF] Korotkov | Computerized analysis of pigmented skin lesions: A review[END_REF].

It is evident that the emphasis has been on assessment of lesion size, shape, color, and texture of lesions. Differential structures, especially streaks, which are distinct dark linear structures, seem to have been neglected by the computer vision community, although they are considered essential for manual lesion categorization. Streaks are brownish-black linear structures of variable thickness, not clearly combined with pigment network lines, mostly at the periphery of the lesion. Irregular streaks are very likely to represent melanoma, especially when the streaks are distributed unevenly.

The detection of streaks in dermoscopy images has been recently reported in [START_REF] Sadeghi | Oriented Pattern Analysis for Streak Detection in Dermoscopy Images[END_REF], [START_REF] Sadeghi | Detection and Analysis of Irregular Streaks in Dermoscopic Images of Skin Lesions[END_REF].

Fig. 1. Examples of images with streaks

In this work, we concentrate on the detection/segmentation of streaks utilizing advanced image processing techniques. Our work involves a simpler approach than the one described in [START_REF] Sadeghi | Detection and Analysis of Irregular Streaks in Dermoscopic Images of Skin Lesions[END_REF]. Furthermore we attempt to evaluate the potential usefulness and diagnostic value of streaks for the differentiation of dermoscopy images into melanoma and non-malignant lesions. More specifically at a first stage we propose a novel algorithm for the segmentation of streaks in dermoscopy images, using the Hessian matrix to obtain local information about the image curvature. Dark linear structures are detected at different scales and their geometric characteristics are used to discriminate between not-streak structures. At a second stage we propose new features based on streaks existence, localization and distribution and investigate their value in computer based dermoscopy image analysis. The rest of the paper is structured as follows:

In Section 2 we describe the technical details of proposed methodology for streaks segmentation and features extraction. Section 3 presents the initial experimental results and Section 4 concludes the paper.

Methodology

The proposed algorithm involves a number of steps, which are illustrated in Fig. 2. Firstly, the dark linear structures are identified and segmented using hessian-base multiresolution approach and hysteresis thresholding. The segmented linear structures are labeled and their geometric characteristics (curvature, position and orientation) are then calculated. Those that satisfy the predefined criteria are labeled as streaks. A number of streak-based features are extracted and used for classification of each dermoscopy image according to malignancy. Each step of the proposed algorithm is presented in detail in the following subsections.

Segmentation of linear structures in images

In this subsection, we describe an algorithm for the segmentation of the dark linear structures in dermoscopy images. The Hessian matrix is used as a measure of linearity for each pixel (x,y) of image I. The use of the Hessian matrix for image feature extraction has been reported frequently in the literature [START_REF] Mikolajczyk | A Comparison of Affine Region Detectors[END_REF], [START_REF] Tuytelaars | Local Invariant Feature Detectors: A Survey[END_REF]. In [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], [START_REF] Aylward | Initialization, noise, singularities, and scale in height-ridge traversal for tubular object centerline extraction[END_REF], the eigenvalues of the Hessian matrix have been used to segment vessels. Let us define the Hessian matrix as following:

( ) ( ) ( ) ( ) ( )
, where the subscript denotes partial differentiation with respect to the corresponding image axis. The image partial derivatives in H are calculated using linear convolution (denoted by * ) with the partial derivatives of a 2D Gaussian function, with variance σ:

( )

2 2 2 2 2 1 , ; 2 x y G x y e σ σ πσ + - ⎛ ⎞ = ⎜ ⎟ ⎝ ⎠ (2)
I xx x, y;σ ( ) = σ γ I x, y ( ) * ∂ 2 G x, y;σ ( ) ∂x 2 (3)
Linear structure segmentation

Linear structure labeling

Identify streaks according to curvature, position and orientation

Streak-related features

Malignancy classification

Hessian-based dark line enhancment

Hysterisis thresholding

The σ parameter is used to select image features at different scales, based on matrix H, depending on the spatial resolution of the dermoscopy image. This is a wellreferenced technique [START_REF] Aylward | Initialization, noise, singularities, and scale in height-ridge traversal for tubular object centerline extraction[END_REF], [START_REF] Martinez-Perez | Segmentation of blood vessels from red-free and fluorescein retinal images[END_REF]. The selected values for σ are discussed in the result section. The image partial derivatives, as described in (3) are normalized with respect to σ by multiplying them with the term σ γ , as described in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], [START_REF] Lindeberg | Edge detection and ridge detection with automatic scale selection[END_REF], so that the resulting values for different values of σ are directly comparable. In this work, we set γ=1. The Hessian matrix has been frequently utilized to describe the local structure of an image. In [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] a number of types of image points are defined according to the eigenvalues of H, in 2 and in 3 dimensions. In [START_REF] Maglogiannis | Enhancing classification accuracy utilizing globules and dots features in digital dermoscopy[END_REF] the Hessian matrix is used to define the diffusivity function in order to enhance dark circular structures in dermoscopy images.

In this work we employ a technique described in [START_REF] Steger | Extracting lines using differential geometry and Gaussian smoothing[END_REF], in order to detect dark linear structures in dermoscopy images. Let λ 1 , λ 2 be the eigenvalues of H at any pixel (x,y) with a specified value of σ. Since H is symmetric, both λ 1 , λ 2 are real numbers. If the pixel (x,y) belongs to a dark linear structure, then one of the eigenvalues of its Hessian matrix is of a positive and high value, whereas the other eigenvalue approaches zero. The eigenvalues λ 1 , λ 2 can be calculated in closed form:

( ) ( ) 1 2 2 , 1 sgn 4 2 q D a q q b D T r D λ λ = = = - + - (4)
where a=1, Tr is the trace and D the determinant of H. Let λ max be the maximum eigenvalue. The values λ max for an exemplary dermoscopy image (Fig. 3a) for the pixels inside the lesion are shown in Fig. 3(b). Pixels with high value of λ max , are candidates for dark linear structures. The eigenvector, corresponding to the maximum eigenvalue, points to the direction of the maximum image curvature. If a pixel belongs to a linear structure, then the eigenvector with λ max is directed perpendicular to the axis of the linear structure. The eigenvector u with λ max is calculated as following:

( ) For any given pixel (x,y), we calculate the exact point coordinates, where the 1 st order directional derivative along u, vanishes and the 2 nd order directional derivative of the image, obtains high absolute value [16 Eqs.(23,24)]:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , , , , , , , ,
x y x y p x y p x y t x y u x y t x y u x y

= ( 6
)
where t is given by

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 , , , , , , , 2 , , , , , x x y
+ = - + + (7)
A pixel (x,y) is considered as belonging to a dark linear structure, if the vanishing of the derivative occurs inside the pixel, or equivalently:

( ) ( ) ( ) ( ) 2
, , , 0.5 , 0.5

x y p x y p x y t x y < ⇔ < , ( 8
)
where 2 stands for the distance metric norm 2. Vector u is illustrated for all pixels that satisfy λ max >0, in the case of a dermoscopy image (Fig. 4(a)), as well as in the case of a simple image containing dark vessels (Fig. 5, left column). We used an image of vessels in order to clearly illustrate the accuracy of the aforementioned algorithm. To demonstrate the use of this algorithm in the case of dermoscopy images, vectors u are superimposed on magnified portions of dermoscopic lesions image, Fig.

Multi-resolution approach

The approach described above is parameterized by the σ parameter of the Gaussian, as shown in Eqs. ( 2) and ( 3). The selection of σ affects the size of the detected linear structures, as described in the theory of scale-space [START_REF] Lindeberg | Edge detection and ridge detection with automatic scale selection[END_REF]. Thus, the above linear structure detection is applied to a given image at a different scales σ i , with i=1,2,…,n , producing different segmented images B i : The segmented images B i , produced by each σ i value, are further refined by applying hysteresis thresholding to the maximum eigenvalues of the image pixels. More specifically, high and low thresholds for λ max (TH, TL respectively) are selected and the following process is performed for each image B i . A non-zero pixel (x 0 ,y 0 ) of B i is preserved if the following condition is met:

( ) ( ) ( ) ( ) 2 , 1 , ; , , ; 0.5
λ max (x 0 ,y 0 ) ≥ TH OR ( λ max (x 0 ,y 0 ) ≥TL AND (x 0 ,y 0 ) is connected to pixel with λ max ≥ TH via pixels with λ max ≥TL)

Otherwise its value is set equal to 0. The final image with the linear segmented structures is obtained as the union of the segmented images B i for the different values of σ, that have undergone hysteresis thresholding:

1 2 1 ... n n i i B B B B B = = ∪ ∪ = U (10)

Following and labeling of linear structures

The previous algorithm creates a binary image, which contains pixels that belong to dark linear structures. The resulting image B is skeletonized using a morphological algorithm. An exemplary dermatoscopy image of a lesion with the segmented linear structures superimposed as black pixels is shown in Fig. 6, (a) before and (b) after skeletonization. The junction pixels and the beginning/ending pixels of each streak are identified. This is performed using the following algorithm, implemented in [START_REF] Kovesi | MATLAB and Octave functions for computer vision and image processing[END_REF].

Let us define the 3x3 neighborhood round each pixel as following:

v 2 v 1 v 8 v 3 v 0 v 7 v 4 v 5 v 6
A binary pixel with value v 0 >0 is a junction if the following condition is fulfilled

1 0 1 6 0 i i i v v ANDv + = - ≥ > ∑ (11)
A pixel with value v 0 >0 is a beginning/ending pixel if the following condition is satisfied [15]

8 1 0 1 2 0 i i i v v ANDv + = - ≥ > ∑ ( 12
)
The algorithm starts from a beginning/ending pixel and follows the segmented pixels until another beginning/ending pixel is found. The visited pixels are labeled as belonging to a linear structure. If a junction is found, the algorithm is re-initialized to track the new structure (as long as non-visited pixels are connected to the junction).

The linear structures outside the lesion of Fig. 3(a) are discarded (Fig. 6,c) and the remaining structures are labeled by the algorithm. The resulting labeling is shown in color in Fig. 6(d).

Selection of streaks

The geometric characteristics of streaks in dermoscopy images are summarized in [7, Introduction section], as following: dark linear structures, "co-radially oriented in the boundary", "darker than their neighborhood", "shorter than the 1/3 of the minor axis of the lesion and they should be longer than one percent of the major axis", and have low curvature. The segmented and labeled linear structures identified by the algorithm so far, are checked against these characteristics, as described below.

First the radial orientation of a linear structure is considered. In [START_REF] Sadeghi | Detection and Analysis of Irregular Streaks in Dermoscopic Images of Skin Lesions[END_REF] this characteristic is checked by assuming the shape of the lesion to be elliptic and by calculating the orientation of the streak with respect to the two foci of the ellipse. However, in this work we utilize a simpler approach: we calculate the apparent angle φ of the streak from the centroid of the skin lesion. Given the first (x 1 ,y 1 ) and last pixel (x N ,y N ) of a labeled linear structure (coordinates with respect to the centroid of the lesion), the apparent angle from the centroid is calculated as following:

1 1 1 tan tan n n y y x x ϕ - - = - ( 13
)
If φ is greater than a predefined threshold φ tol then the linear structure is discarded. The requirement for linear structures darker than the background is fulfilled by requiring positive eigenvalues of the H matrix. The upper and the low limit of the length of each linear structure require the calculation of the major and minor axis of the binary lesion, which is performed using the moments of inertia of the 2 nd order. An example of the calculation of the major axis is shown in Fig. 8. Thus a labeled linear structure is discarded if its length L is smaller than 0.01 of the length of the major axis of the lesion L 1 , or if its length is greater than 1/3 of the minor axis L 2 of the lesion:

( ) ( ) 2 2 2 2 1 1 1 0.01 , 3 N N L L ORL L L x x y y > < = - + - (15)
Since the curvature of a linear structure should be small to be considered as streak, its length L is calculated as the distance between beginning and ending point. L is compared to major axis length in order to fulfill the relevant rule [START_REF] Sadeghi | Oriented Pattern Analysis for Streak Detection in Dermoscopy Images[END_REF].

φ Fig. 8. An example of the calculation of the major axis of a binary lesion (continuous line). The centroid of the lesion is marked by a square.

In this study we seek to identify linear structures that resemble straight lines, thus structures with increased curvature need to be discarded. There are many approaches for curvature estimation of pixelized structures. We have utilized the simple but efficient approach implemented in [START_REF] Kovesi | MATLAB and Octave functions for computer vision and image processing[END_REF], described as following:

For each linear binary object containing N pixels {(x i ,y i )}, i=1,2,…,N, the equation of the line defined by its 1 st and last pixel (x 1 ,y 1 ), (x N ,y N ) is considered and the maximum distance d max of the distances d i of each pixel (x i ,y i ) from the line was calculated:

( )

1 1 max , max i i i N N i d x y y x C C y x y x d d δ δ = + - = - = (16)
If d max is greater than a predefined threshold tol then the linear structure is discarded.

tol d d > (17)
According to [START_REF] Sadeghi | Oriented Pattern Analysis for Streak Detection in Dermoscopy Images[END_REF], [START_REF] Sadeghi | Detection and Analysis of Irregular Streaks in Dermoscopic Images of Skin Lesions[END_REF], streaks in dermoscopy images appear close to the boundary of the lesion. In this work, this characteristic is implemented by calculating the distance transform (DT) of the binary lesion. A linear structure is discarded if the maximum distance of its two beginning/ending pixels from the boundary is greater than the (1/3) of the length of the minor axis L 2 of the lesion. Equivalently, (

Accepted streaks

The linear structures that remain after discarding structures according to criteria in Eqs. ( 14), ( 15), ( 17) and ( 18) are considered as streaks.

Extraction of streak-based features

In [START_REF] Kovesi | MATLAB and Octave functions for computer vision and image processing[END_REF] streaks are classified as regular or irregular according to their spatial distribution. In this work we define the angular and radial asymmetry of the spatial distribution of streaks with respect to the centroid of the lesion, as following. The centroid of each streak c i is located, defined as its median pixel. The azimuthial angle of the centroid θ i and its distance from the boundary r i =DT(c i ) (using the distance transform DT of the binary lesion) is calculated for each streak i=1,2,…,K. The histogram of θ i and r i is calculated (H θ and H r respectively) using 10 bins. The angular and radial asymmetry (A θ and A r respectively) are calculated as following:

( ) ( ) ( ) ( ) ( ) ( ) max min max max min max r r r r H H A H H H A H θ θ θ θ - = - = (19)
Both A θ and A r are normalized in the [0,1], with values approaching 1 indicating high asymmetry. The number of detected streaks, as well as the number of pixels in detected streaks is also calculated and used as a feature.

Experimental Results

The proposed algorithm was applied to 64 dermoscopy images with non-malignant lesions and to 35 images of malignant lesions, determined by histological analysis. All images have size of 632 × 387 pixels and the same magnification, resulting in spatial resolution of about 0.05 mm/pixel. They were acquired using the ELM Molemax II device at Hospital of Wien and at the department of Plastic Surgery and Dermatology in the General Hospital of Athens G. Gennimatas. The area of the lesion has been manually delineated for each one of the images. The same parameters of the algorithm were used for all images: d max =2 pixels, φ tol =2 deg, TH=0.1, TL=0.01 (with respect to maximum λ max value). The σ parameter, described in subsections 2.1 and 2.2, takes integer values of 1, 2 and 3. The detected streaks are shown in Fig. 10 for exemplar non-malignant dermoscopy images (left column) and on malignant lesions (right column). An initial effort to investigate the diagnostic value of streaks for the detection of malignancy in the skin lesions was undertaken. More specifically, the four streak-based features (described in Section 2.5) were used for classifying the dermoscopy images into 2 classes: "malignant" and "non-malignant". The ground truth of the specific lesion was determined by histological analysis as already mentioned. We have experimented with the following classifiers: the kNN (with k=5), the linear and quadratic discriminant, support vector machines (using Sequential Minimal Optimization), as well as feed forward neural networks with back error propagation, considering different number of neurons in the hidden layer. In addition to accuracy, sensitivity and specificity (true negative rate, complementary to false positive rate) are also presented. Sensitivity is of importance, since by its definition (true positive rate) is complementary to false negative rate, which is required to be low by any medical decision support system. The reported results prove the relevance of the streak-based features with skin lesion malignancy. The complexity of the proposed algorithm is analogous to the total number of pixels in the dermoscopy images, as well as analogous to max(σ i ) 2 , used in Eqs. ( 2), (3). Execution time was approximately equal to 2 seconds, for each dermoscopy image, with the parameters described above, using Matlab on an Intel(R) Core i5-2430 CPU @ 2.40 GHz Laptop with 4 GB Ram, under Windows 7 Home Premium.

Discussion and Conclusion

An algorithm for the detection of streaks in dermoscopy images has been presented. The proposed methodology segments dark linear structures using Hessian-based multi-resolution image features. Qualitative results show successful detection of streaks. The quantitative results, concerning the classification of lesion into malignant or nonmalignant, according to features related to the detected streaks, appear promising. Although the support vector machines achieved the best overall accuracy (75.11 %), the highest sensitivity (equivalently lowest false negative rate) was achieved by the quadratic discriminant classifier. It is in our future plans to continue this work by extending the dataset, as well as defining and utilizing more streak-based features, possibly in conjunction with other, already established features. Thus, we believe that the classification accuracy and sensitivity may increase further.

Fig. 2 .

 2 Fig. 2. Brief outline of the proposed algorithm

Fig. 3 .

 3 Fig. 3. The value of λ max for the pixels inside the lesion (b) of an exemplary dermoscopy image (a).

  u x y I x y u x y t x y I x y u x y I x y u x y u x y I x y u x y

4

 4

  (a), (b). Yellow dots in both examples indicate pixels that satisfy the criterion of Eq.(8).

Fig. 4 .

 4 Fig. 4. Vectors u is superimposed on the pixels that satisfy λ max >0, for magnified portions (a) and (b) of the dermoscopy image in Fig.3(a), using σ=3.

  Segmented images B 2 and B 3 (for σ=2 and σ=3 respectively) for the vessel image are shown in the right column of Fig.5.

Fig. 5 .

 5 Fig. 5. Illustration of calculation of eigenvector for maximum eigenvalue (left column) in the case of an image with dark vessels, for σ 2 and 3 (upper and lower line, respectively). Right column shows the resulting vessel segmentation.

Fig. 6 .

 6 Fig. 6. The segmented linear structures (black pixels) from a dermatoscopy image of a lesion, (a) before and (b) after skeletonization. Skeletonized structures inside the lesion (c). The labeled linear structures that lie inside the lesion are shown in color in (d).

Fig. 7 .

 7 Fig. 7. An exemplar dermoscopy image with the detected streaks. The apparent angle φ of a streak from the center of the lesion (marked by "*") is shown.

Fig. 9 .

 9 Fig. 9. The detected streaks from a dermoscopy image are shown, with the border of the lesion, as well as the distance transform of the border (encoded as shades of gray inside the lesion).

Fig. 10 .

 10 Fig. 10. Streaks detected by the proposed algorithm on non-malignant (a) and on malignant lesions (b). Images have been cropped, but are all in the same scale.

Table 1 .

 1 The classification results using the four streak-based features.The training set was 32 instances from non-malignant and 22 cases from malignant class. The achieved accuracy, sensitivity and specificity are shown in Table1, averaged over 50 executions.

	Classification Method	accuracy	sensitivity	specificity
	kNN, k=5		0.7022	0.600	0.7437
	kNN, k=7		0.7040	0.3953	0.8294
	Linear Discriminant	0.7249	0.7538	0.7131
	Quadratic Discriminant	0.7000	0.7877	0.6644
	Support Vector Machines	0.7511	0.6983	0.7731
	Neural Nets	Hidden 9x1	0.6051	0.5757	0.6847
	Neural Nets	Hidden 3x1	0.5829	0.5804	0.5894

Acknowledgments

This work has been supported by the 12CHN 204 Bilateral Greece-China Research Program of the Hellenic General Secretariat of Research and Technology and the Chinese Ministry of Research and Technology entitled "Personalization of melanoma therapeutic management through the fusion of systems biology and intelligent data mining methodologies-PROMISE," sponsored by the Program "Competitiveness and Entrepreneurship," Priority Health of the Peripheral Entrepreneurial Program of Attiki.


Ioannis Valavanis

Theodosios Goudas

Maria Michailidou

Ilias Maglogiannis

Heleni Loutrari

Aristotelis Chatziioannou

G P Livanos

M Simou

A novel image analysis methodology for the evaluation of angiogenesis in matrigel assays and screening of angiogenesis-modulating compounds

Keywords: angiogenesis, tubules, matrigel, RGB image, image analysis

Angiogenesis, with a prominent role in tumor growth, is a physiological process through which new blood vessels are generated from preexisting vasculature. The effect of angiogenesis inhibition test agents to cell differentiation can be evaluated using the matrigel differentiation assay, an environment simulating in vitro the natural angiogenic process under study. Using this assay, angiogenic activity can be evaluated by the quantification of tubules, junctions and other structures of interest in the matrigel images, i.e. microphotographs of the formed cellular networks under different fields of view. In this work, we present an image analysis workflow that uses separate channels in the RGB matrigel images, each contributing different information related to angiogenic activity (tubules, junctions) and other structures (cells, lumens). The workflow is successfully showcased on the basis of a series of matrigel images generated for the evaluation of angiogenesis inhibition compounds, and results are cross-validated with a manual expert scoring.

Introduction

Angiogenesis is a physiological complex process through which new blood vessels are generated from pre-existing vasculature. Endothelial cells, the building units of vasculature, form the internal lining of blood vessels. Endothelial cell invasion, migration, proliferation, differentiation and structural rearrangement are essential steps of the angiogenic process that are regulated by a balance between pro-and antiangiogenic factors. However under certain conditions, this balance is disrupted resulting in either poor vascularization, as in the case of ischemic disorders, or in abnormal vasculature, as in the case of tumors [START_REF] Potente | Basic and Therapeutic Aspects of Angiogenesis[END_REF].

The prominent role of angiogenesis in tumor growth and metastasis was first described by Folkman in 1971 [START_REF] Folkman | Isolation of a Tumor Factor Responsible for Angiogenesis[END_REF]. Since then, scientific literature has offered a massive information on the role of angiogenesis in the progression of several human diseases and new therapeutic strategies for their treatment have been proposed through targeting of the angiogenic process [START_REF] Bridges | The Angiogenic Process as a Therapeutic Target in Cancer[END_REF] [START_REF] Schmidt | Angiogenesis: a target in solid tumors, also in leukemia?[END_REF].

One of the challenges researchers initially faced was the selection of scientifically and technically efficient methods for assessing the endothelial response to candidate angiogenesis inhibitors [START_REF] El-Kenawi | Angiogenesis inhibitors in cancer therapy: mechanistic perspective on classification and treatment rationales[END_REF]. Ideally, an assay should be valid, methodologically straightforward, easily quantifiable, cost-effective and able to provide an environment simulating the natural angiogenic process under study. One of the assays that encases several of the afore-mentioned factors and has been extensively used is the capillarylike tube formation assay, otherwise known as matrigel differentiation assay. It is considered a suitable in vitro angiogenesis model as it examines the ability of endothelial cells to produce tubular structures when cultured on surfaces pre-coated with matrigel. It is a specialized commercially available polymer substance extracted from Engelbreth-Holm-Swarm murine sarcoma. Matrigel contains basement membrane components and mimics the conditions of the external tissue microenvironment of the endothelium [START_REF] Hughes | Matrigel: a complex protein mixture required for optimal growth of cell culture[END_REF]. This extracellular matrix activates differentiation pathways in endothelial cells enabling them to form networks that display tubular extensions of the cytoskeleton and dynamic rearrangements of the cells resembling vascular lumens. Microscopes are used to capture phase-contrast images of the angiogenesis environment in vitro [START_REF] Ko | In vitro Human Umbilical Vein Endothelial Cells (HUVEC) Tube-formation Assay[END_REF]. The cells can be treated with test agents and their effects on endothelial cell differentiation are evaluated using quantification of tube number and length, manually or by image processing software techniques applied to the images obtained by the microscope [START_REF] Doukas | Computer-supported angiogenesis quantification using image analysis and statistical averaging[END_REF].

Wimasis and Angiosys are two commercial software packages that are dedicated to angiogenesis quantification in matrigel images, each featuring certain specifications and capabilities [START_REF] Khoo | A comparison of methods for quantifying angiogenesis in the Matrigel assay in vitro[END_REF]. In order that matrigel images are used by Angiosys, they need to be in prior processed (skeletonized) by a dedicated image-processing tool (e.g. Adobe Photoshop combined with an Image Processing Kit plug-in as described in [START_REF] Khoo | A comparison of methods for quantifying angiogenesis in the Matrigel assay in vitro[END_REF]). Wimasis performs automatically the whole process of image processing and angiogenesis quantification. Most frequently extracted angiogenesis features include tubules length, number of tubules (separate paths in the angiogenic pathway), number of junctions (cross-points of tubules) and tubules area [START_REF] Khoo | A comparison of methods for quantifying angiogenesis in the Matrigel assay in vitro[END_REF]. In [START_REF] Gallo-Oller | Quantitative method for in vitro matrigel invasiveness measurement through image analysis software[END_REF], the red channel of matrigel image is used within the open source ImageJ software package to measure the area that cells cover in the image. A semi-automatic manual approach for measuring tubules and cells in matrigel images using the Scion Image, an image processing software freely available from the National Insitute of Health (NIH) website, is presented in [START_REF] Donovan | Comparison of three in vitro human 'angiogenesis' assays with capillaries formed in vivo[END_REF].

In this work, we present a fully automatic workflow, written in Java and utilizing the ImageJ 1.47 library, for matrigel image processing and analysis, in order to quantify the angiogenesis presented in matrigel images and extract other useful structures in these images (cells, lumens). Differing from other approaches, we neatly use separate channels in the RGB matrigel images, each contributing different information related angiogenic activity (tubules, junctions) and other structures (cells, lumens) in order to maximize the efficiency of capturing the information content within the matrigel images. We showcase successfully the use of the workflow for the analysis of a series of matrigel images, used for the evaluation of angiogenesis inhibition activity of several compounds, and cross-validate results with a manual expert scoring. The rest of the paper is structured as follows: Section 2 presents matrigel images acquisition, along with sample images, and the image-processing workflow, while Section 3 presents and discusses the results.

Materials and Methods

In this Section, matrigel images and their formulation using the matrigel capillarylike tube formation assay are firstly described. Next, three image processing subworkflows, each detecting different structures in the matrigel images, are described.

Matrigel capillary-like tube formation assay

Matrigel capillary-like tube formation assay was performed as previously described in [START_REF] Michailidou | Microvascular Endothelial Cell Responses in vitro and in vivo: Modulation by Zoledronic Acid and Paclitaxel?[END_REF]. Briefly, serum-starved human umbilical cord vein endothelial cells (HUVEC) were cultured in 24-well plates pre-coated with matrigel in presence of test compounds SC1117, SC1119, SC1123 and SC1124 (0.2 -2 μM). At the end point, cells were stained with Eosin-Haematoxylin and microphotographs, using a Leica microscope (20x magnification lens), of the newly formed cellular networks were obtained from 10 different fields of view [START_REF] Michailidou | Microvascular Endothelial Cell Responses in vitro and in vivo: Modulation by Zoledronic Acid and Paclitaxel?[END_REF]. Thus, a total of 10 matrigel images were extracted and analyzed per sample (each image was 2048x1532 pixels, 0.68 μm/pixel). Evaluation of the anti-angiogenic effects of the test compounds took place by image processing and analysis workflow described next, that measured various phenotypical characteristics of cellular differentiation such as number of cells, tubular length, cellular junctions where tubules may cross-point, and lumens (the inside of a tubular structure). Manual scoring of tubule number was also performed for crossvalidation of the image analysis data.

A sample matrigel image is presented in Fig. 1 (top). This image is the one chosen to showcase the results of the image processing analysis workflow in the next section.

Matrigel Image Processing and Features Extraction

Three image processing and analysis sub-workflows were implemented, each of which using a combination of popular image processing operators, processes a different channel in the RBG matrigel images and derives different features from the image. The image processing sub-workflows are described below.

-Blue Channel Processing for angiogenesis path and tubules' extraction: The workflow is presented in Fig. 2 and processes the blue channel containing the angiogenesis activity both in tubules and around/within junctions in order to segment the angiogenesis path. Image parts that are featured a closed contour of the angiogenesis path are characterized as lumens. Part of the angiogenic path beginning and closing with a junction or cell (see next workflows for their identification) is characterized as tubule. Furthermore, an automatic thresholding in the blue channel and the consequent calculation of the ratio of black to total image is used to derive an % angiogenesis quantification in the image. Blue channel can be processed differently in order to identify the cells in the matrigel image, as well. The corresponding sub-workflow is presented in Fig. 3.

-Red Channel Processing for junctions' extraction: Collagen found in the junction can be clearly visible in the red channel in matrigel images and the sub-workflow used for this scope is presented in Fig. 4 The results of each of the workflows described above are merged in a finally resulting image, where tubules, junctions and cells are shown (see Fig. 1, bottom for the resulting image obtained from the original matrigel image presented in Fig. 1, top). Numerical results provided include quantification of angiogenesis (%), number of lumens, number of cells, number of junctions, total length of angiogenesis path (number of pixels) and ratio of total area found in cells (%).

Results and Discussion

The matrigel images extracted to evaluate the inhibition activity of the test compounds SC1117, SC1119, SC1123 and SC1124 (0.2-2μM) and the control sample images (no inhibition compound applied) were analyzed. Each compound and the control sample were evaluated on the basis of the available images corresponding to the fields of the formed cellular network. The workflow could real-time process the images and identify successfully structures of interest: tubules, cells, junctions and lumens. Sample imaging results are presented in Fig. 1 (bottom), where coloring is used to annotated the identified structures. Other imaging results not shown here support the efficiency of the algorithm to capture the structures of interest in the matrigel images.

Numerical results for each sample, calculated as average values of the features extracted for all images corresponding to the 10 fields of views, are presented in Table 1 (see Section 2.2 for description and extraction of the features used). A careful observation of the entries in Table 1 shows that the application of the inhibition compound reduces for most cases the angiogenic activity as measured by the features extracted that actually measure this activity, i.e. quantification of angiogenesis, total length of angiogenesis path, number of junctions. This appears more intensively for the increased concentration for the inhibition compounds (0.2μΜ vs. 2μΜ). Number of cells and lumens complete the angiogenesis information and can be used additionally by the expert. Angiogenic activity reduction as measured by the quantification of angiogenesis (%) is presented graphically in Fig. 5, showing the effects of the various compounds and their concentration compared to the control sample. Manual scoring results for the evaluation of angiogenesis are presented in Fig. 6. Specifically, the tubule number in the matrigel images after the effects of the various compounds and their concentration, are presented as % of tubule number measured in the control sample. Manual scoring results show the same pattern of angiogenesis reduction as implied by results regarding the automatic quantification of angiogenesis as presented in Fig. 5. Specifically, both manual and automatic analysis show that a greater concentration for the inhibition compounds results in a greater reduction in angiogenesis and that SC117 compound inhibits mostly the angiogenesis among all compounds tested. Future work will include the evaluation and settings optimization of our image processing and analysis workflow using a larger dataset of images and the final evaluation of various natural or synthetic inhibition compounds on the basis of the matrigel capillary-like tube formation assay. A sophisticated analysis of the features extracted from a larger matrigel images dataset, will set the basis for a robust fully automated analysis workflow for the evaluation of endothelial response to candidate angiogenesis inhibitors. Quantification of angiogenesis (%)

Fig. 1 Fig. 2

 12 Fig. 1 Top: Original sample matrigel image. Bottom: Result of the image processing workflow (green corresponds to tubules, red color contours junction, and yellow contours cells)

Fig. 3

 3 Fig. 3 Processing blue channel in matrigel image for cells' extraction

Fig. 4

 4 Fig. 4 Processing red channel in matrigel image for junctions' extraction

Fig. 5

 5 Fig. 5 Quantification of angiogenesis (x 100 %) showing the effects of the various compounds and their concentration compared to the control sample.

Fig. 6

 6 Fig.[START_REF] Hughes | Matrigel: a complex protein mixture required for optimal growth of cell culture[END_REF] Manual scoring results for the evaluation of angiogenesis (tubule number in the matrigel images after the effects of the various compounds are presented as % of tubule number measured in the control sample)

Table 1 .

 1 Average values of features calculated for all 10 fields of views for each sample (control and five test compounds for concentrations 0.2μΜ and 2μΜ)

		Quantification of			Total length of
	Sample	angiogenesis (x 100 %)	Number of Junctions	Number of Cells	angiogenesis path (pixels)	Number of Lumens	Ratio of total area found in cells
	Control	0.096	21.769	71.692	34420	12.769	0.012
	SC117 0.2μM	0.062	14.909	79.455	31880	9.273	0.014
	SC117 2μM	0.056	9.600	92.200	30033	6.200	0.017
	SC119 0.2μM	0.091	25.000	100.600	34431	14.000	0.015
	SC119 2μM	0.062	12.500	124.500	32184	7.100	0.021
	SC123 0.2μM	0.072	20.100	87.900	35826	13.300	0.015
	SC123 2μM	0.070	18.818	91.000	34944	10.818	0.016
	SC124 0.2μM	0.071	21.800	59.500	30668	9.400	0.011
	SC124 2μM	0.063	16.700	64.400	33088	9.300	0.012

Acknowledgement

This work was supported by the Bilateral Greece-China Research Program of the Hellenic General Secretariat of Research and Technology and the Chinese Ministry of Research and Technology PROMISE (12CHN204, Personalization of melanoma therapeutic management through the fusion of systems biology and intelligent data mining methodologies).


Zohreh Alizadeh
email: z.m.alizadeh.b@gmail.com

Mohsen Afsharchi
email: afsharchim@znu.ac.ir

Armin Ghasem Azar
email: armin.ghasemazar@gmail.com

A Distributed Energy-Efficient Algorithm for Cellular Target Tracking in Wireless Sensor Networks

Keywords: Wireless sensor networks, target tracking, energy efficiency

Energy efficiency is one of the most challenging issues in wireless sensor networks, particularly in the target tracking. The main purpose of these networks is to preserve the distributed important sites from the targets who intend to destroy the sites in an environment. This paper proposes a distributed energy-efficient mechanism for sleep scheduling the sensors followed by a dynamic cellular clustering algorithm for tracking the targets. Probabilistic positions of the targets are predicted based on an improved particle filter procedure consecutively. A cell, including a subset of sensors, is constructed over time considering the predicted positions of targets. This cellular algorithm not only decreases the number of awake sensors, but also increases the tracking quality. Moreover, a concept named communication base has been proposed for alleviating the communication volume of the sensors. Utilizing communication bases, sensors are able to exchange information in maximum three hops. The experimental results demonstrate the capability of the provided algorithms and mechanism in optimizing the energy of sensors and increasing the tracking quality.

Introduction

A wireless sensor network is a network of spatially-distributed autonomous sensors for monitoring physical or environmental conditions to cooperatively disseminate their data through the network. Sensors are generally equipped with the capabilities of sensing, processing, and communicating. One of the most important issues in wireless sensor networks is energy consumption of sensors since these networks consist of a large number of battery-limited sensors which are typically used in the high-risk areas and thus the battery replacement is practically impossible [START_REF] Jiang | Probability-Based Prediction and Sleep Scheduling for Energy-Efficient Target Tracking in Sensor Networks[END_REF]. The importance of energy efficiency of sensors is further significant in target tracking, as one of the most important applications of wireless sensor networks [START_REF] Zhuang | Minimizing Energy Consumption With Probabilistic Distance Models in Wireless Sensor Networks[END_REF]. Target tacking in wireless sensors networks consists of accumulating information about target movement by sensors and delivering related information to other sensors or a base station to make an appropriate decision about the targets. Idle listening is a major source of energy waste for target detection and measurement [START_REF] Jiang | Probability-Based Prediction and Sleep Scheduling for Energy-Efficient Target Tracking in Sensor Networks[END_REF]. In particular, sleep scheduling is one of the most commonly-used mechanisms for reducing the energy consumption during idle listening in order to extend the network lifetime [START_REF] Wang | Cluster-Based Dynamic Energy Management for Collaborative Target Tracking in Wireless Sensor Networks[END_REF]. The idea of sleep scheduling is to put sensors in the sleep state for most of the time and only wake them up periodically. In certain cases, the sleep pattern of sensors may also be explicitly scheduled (i.e. forced to sleep or awaken on demand). Additionally, sensor clustering has been demonstrated to be an energy-efficient structure in sensor networks since information routing and relaying are only operated by special sensors namely cluster heads [START_REF] Fuemmeler | Smart Sleeping Policies for Energy-Efficient Tracking in Sensor Networks[END_REF]. In each time step, a number of the sensors create a cluster and a sensor with a specific qualification becomes cluster head. Cluster heads can aggregate and process the information sent by cluster members, thus reduce the network load and alleviate the bandwidth.

This paper gives an account for utilizing a multi-objective optimization technique for localizing the sensors in the environment as a precursor for the main proposed algorithm. The assumed rectangular-shape environment is composed of some entrance gates and several solid obstacles followed by a set of important sites. Minimizing the number of sensors while maximizing the network coverage in the environment are the assumed objectives of the used technique which are obviously in conflict with each other. Furthermore, a distributed energy-efficient algorithm for sleep scheduling the sensors has been proposed followed by a dynamic cellular clustering procedure for effectively tracking the targets. In naive clustering approaches, a cluster head is designated to each cluster which is responsible for gathering the information from sensors and reporting it to other cluster heads. However, considering the fact that cluster heads must be awake continuously, their limited energy will not last long. Hence, this paper has substituted the cluster head with the concept of cell while the role of cluster head has been modified with constant communication bases. Furthermore, an improved particle filter algorithm based on the radial basis function networks has been used for predicting the trajectory of the targets through the time [START_REF] Wang | An Improved Particle Filter for Target Tracking in Sensor Systems[END_REF]. In each time interval, a cell is created based on the information gathered by the sensors and probabilistic positions of the targets in the previous time interval.

The rest of this paper is organized as follows. Section 2 contains a review of the related work. In Section 3, the comprehensive description of the system model is presented. The proposed algorithm followed by its required procedures are members of Section 4. Section 5 demonstrates experimental results. Finally, conclusion and future work are drawn in Section 6.

Related Work

Recent developments in target tracking have heightened the need for effectively maximizing the target tracking quality [START_REF] Jiang | Probability-Based Prediction and Sleep Scheduling for Energy-Efficient Target Tracking in Sensor Networks[END_REF][START_REF] Zhuang | Minimizing Energy Consumption With Probabilistic Distance Models in Wireless Sensor Networks[END_REF]. Jiang et al. [START_REF] Jiang | Probability-Based Prediction and Sleep Scheduling for Energy-Efficient Target Tracking in Sensor Networks[END_REF] applied a probability-based prediction mechanism in order to increase the tracking accuracy and reduce the sensors' active time by predicting the position of each target more precisely. Zhuang et al. [START_REF] Zhuang | Minimizing Energy Consumption With Probabilistic Distance Models in Wireless Sensor Networks[END_REF] proposed an energy-efficient algorithm to accurately minimize the energy consumption of sensors in wireless sensor networks taking the probabilistic distance distributions into account. However, the number of hops needed to track the targets was high.

As one of the sleep scheduling approaches, proactive wake-up has been studied for awakening sensors pro-actively to prepare for the approaching target [START_REF] Wang | Cluster-Based Dynamic Energy Management for Collaborative Target Tracking in Wireless Sensor Networks[END_REF][START_REF] Fuemmeler | Smart Sleeping Policies for Energy-Efficient Tracking in Sensor Networks[END_REF]. Wang et al.

[3] represented a low-energy communication structure while the energy efficiency of wireless sensor networks is enhanced by cluster-based dynamic energy management. However, the role of cluster head was idle in the cluster which resulted in reducing the energy of cluster head and target loss. Furthermore, the sleeping policy introduced by Fuemmeler et al. [START_REF] Fuemmeler | Smart Sleeping Policies for Energy-Efficient Tracking in Sensor Networks[END_REF] performed well in energy conservation and tracking the location of the target within the accuracy of the radius of the sensor, but there was a trade-off between the energy savings and the tracking errors.

Particle filtering for target tracking in wireless sensor networks has already attracted some attention [START_REF] Djuric | Multiple Particle Filtering[END_REF][START_REF] Wälchli | Distributed Event Localization and Tracking With Wireless Sensors[END_REF]. Recently, Djuric et al. [START_REF] Djuric | Multiple Particle Filtering[END_REF] have proposed a class of interacting particle filters for multi-target tracking. This class of algorithms relied on splitting the statespace into lower dimensional sub-spaces in order to become computationally tractable, but did not guarantee that the particles were assigned proper weights. Additionally, the velocity of the targets has been assumed to be fixed. In real world, most of the times, the targets move in the environment with different velocities depending on the situations. DELTA was the distributed algorithm proposed by Wälchli et al. [START_REF] Wälchli | Distributed Event Localization and Tracking With Wireless Sensors[END_REF] which tracked the target solely at constant speed by dynamic clustering and selection of cluster head based on the light measurement. The advantage of DELTA was that the communication range of the sensors was higher than their sensing range. However, the main challenge of this approach was that it could only deal with targets' constant velocity.

As a result, according to the significant issues mentioned in the related work, followings are the main contributions of this paper:

• Proposing a new sensor cluster formation based on a new concept named cell; • Providing a cellular algorithm for optimizing the energy consumption of sensors;

• Combining a multi-objective evolutionary technique for optimally localizing the sensors with the proposed energy-efficient target tracking algorithm; • Employing a set of communication bases for effectively decreasing the communication volume of sensors by minimizing the hops in information exchange;

System Model

In this section, the proposed system model including its elements and running procedures will be described. It is assumed that there is an environment, e.g., airport, city, etc. The dimension of the environment is M × N , such that M, N ∈ R + . This environment has G ∈ N entrance gates and O ∈ N various-shaped planar obstacles in this environment which it is not possible to access into them. Moreover, there are T ∈ N important sites such as school, hospital, etc. In this paper, individuals or objects, including minatory purposes, are named as the targets.

Wireless Sensor Networks: The main purpose in this paper is to track and then, destroy the targets which are trying to reach the important sites. Sensors are scattered in the environment such that targets are tracked with a high accuracy and low loss rate. In this paper, sensors have been assumed as radar-based ones which invisibly sense the targets. One sensor begins to sense in its sensing radius and communicate with other sensors within its communication radius. However, each sensor has a cost and limited amount of energy. Therefore, it is desirable to use the minimum number of sensors.

Here, due to simplicity, the purchasing cost of sensors has not been considered. Equation (1) refers to the first objective which is reducing the number of sensors.

f (s) = min S. (1)
S ∈ N denotes the number of sensors. Each sensor s i , i ∈ (1, 2, ..., S), has a sensing radius R s i ∈ R + and a communication radius R c i ∈ R + . This paper assumes that communication radius of each sensor R c i ∈ R + is at least twice its sensing radius R s i ∈ R + . Next, (2) indicates the total overlapping area based on the pairwise overlap between each two independent sensors.

E = - 1 2 × L s i -L s j × 4 × (R s i ) 2 -L s i -L s j 2 -2 × (R s i ) 2 × tan -1     L s i -L s j 4 × (R s i ) 2 -L s i -L s j 2     + π × (R s i ) 2 .
(

) 2
In this paper, it has been assumed that all sensors have the same sensing and communicating radius. L s i is the location of sensor s i . As a result, (3) represents the second objective in the sensor placement problem, which implies maximizing the covered area. o i is the obstacle o i in the environment, where |o i | is its area.

g(s) = max S i=1 π × (R s i ) 2 + O i=1 |o i | -E . (3)
For expressing the function related to the second objective which is increasing the area covered by sensors, the main challenge is the distance between locations of each two sensors. Fig. 1 reflects the status of sensors with various possible assumptions. Fig. 1(a) demonstrates the distance between two sensors which is twice the sensing radius. It is obvious that these two sensors have no intersection with each other. Now, if another sensor is also placed next to them, a status similar to Fig. 1(b) is obtained. It is clear that the space between these three sensors is a blind space or non-covered space and ultimately increases the risk of losing the target. Furthermore, it is assumed that placing the senors inside/on obstacles and important sites is not allowed. As a result, based on the experiments done and according to Fig. 1(c), the distance between each two independent sensors should not be less that 1.5 times the sensing radius.

Since satisfying the mentioned objectives and constraints at the same time is impossible, multi-objective optimization techniques have been used [START_REF] Deb | Multi-Objective Optimization[END_REF]. Multi-objective optimization is an area of multiple criteria decision-making, where mathematical optimization problems involving more than one objective function should be optimized, simultaneously. Optimal decisions need to be taken in the presence of trade-offs between two or more conflicting objectives. Solving a multi-objective optimization problem necessitates computing all or a representative set of Pareto-optimal solutions. Here, each Pareto-optimal solution includes locations of a random number of sensors. Evolutionary algorithms are one of the most well-known meta-heuristic search mechanisms utilized for multi-objective optimization problems since their structure is free of search space and goal capacities [START_REF] Deb | Multi-Objective Optimization Using Evolutionary Algorithms[END_REF]. Evolutionary algorithms involve biological mechanisms such as reproduction, mutation, recombination, natural selection and survival of the fittest. The main advantage of evolutionary algorithms, when applied to solve multi-objective optimization problems, is the fact that they typically generate sets of solutions, allowing computation of an approximation of the entire Pareto-front. Each Pareto-front includes a set of Pareto-optimal solutions. Currently, most evolutionary multi-objective optimization algorithms apply Pareto-based ranking schemes such as the Non-dominated Sorting Genetic Algorithm-II (NSGA-II) [START_REF] Deb | A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF]. Applying this algorithm, S sensors with the same sensing radius are dispensed in the environment. Fig. 2 presents the arrangement of the sensors obtained from applying the NSGA-II including four obstacles and two important sites. The difference between Fig. 2(a) and Fig. 2(b) is the sensing radius and number of sensors. Increasing the sensing radius, the number of sensors decreases and consequently, the covered area also decreases. However, if the radius is too high or too low, the total covered area is close to the maximum. Due to the fact that the implementation and procedure of the NSGA-II is a pre-process for the main proposed algorithm, its description has been disregarded [START_REF] Deb | A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF]. Energy Consumption of a Sensor: Each sensor is able to do a sort of actions and consumes some energy for performing each of the actions. Table 1 includes the energy each sensor consumes for performing the actions. Each sensor is equipped with a battery which has a limited energy. The reason for using battery-powered sensors instead of battery-less sensors is their high prices. As the battery goes out of energy, the sensor will break down and result in network failure and information loss. The actions 1 to 5 outlined in Table 1 are the ordinary actions each sensor performs in the environment. Two additional actions have been taken into consideration in this paper, namely, waiting for receiving the messages and delay in doing the actions 1-5.

E I ∈ R + (mJ/s
) denotes the energy the sensor consumes in sleep (Idle) state. Each sensor spends E S ∈ R + (mJ/s) when sensing the environment. E T ∈ R + (mJ/s) is the energy the sensors employ for changing their states from sleeping to sensing mode and vice versa. In addition, each sensor consumes an amount of energy to communicate with other sensors which corresponds to E C ∈ R + (mJ/s). Finally, E P ∈ R + (mJ/s) demonstrates the energy that sensors use for predicting the targets' location, velocity, and direction. Although doing each action types 1 to 5, mentioned in Table 1, has only a little delay, however, aggregation of these little delays may cause a tracking gap in the system. E D ∈ R + (mJ/s) implies the energy each sensor consumes when it has a delay in performing the mentioned actions. Furthermore, each sensor waits a limited time for receiving any message from its neighbors regarding the wireless communication and transmission delays. Sensors apply the E W ∈ R + (mJ/s) when they are waiting for receiving information from their neighbor sensors. The state diagram of energy variations for each sensor with respect to Table 1 can be seen in Fig. 3.

The loop in E D si states that the sensor may remain in delay state continuously. This state is also applicable to E T si and E S si . The edge between E S si and E P si also means that the sensor can predict the target after sensing the environment. As another instance, according to the edge between E T si and E C si , the sensor begins to communicate with other neighbor sensors after waking up. As mentioned before, most of the energy of sensors is consumed for sending messages and communicating between sensors. Since sensors have a limited sensing and communication radius, therefore, they need to send the information they have collected via other sensors in the environment. To decrease this communication volume, a new concept namely communication base is added to the environment as follows.

Communication Bases: Sensors are not able to communicate with all other sensors due to their limited communication radius. In other words, if a sensor needs to communicate and send messages to other sensors which are far from it, it has to send the message through other neighbor sensors that probably results in wasting the energy and loosing the packets across the path. For this purpose, B ∈ N communication bases is utilized. The bases are connected in a complete graph. As Fig. 4 illustrates, communication bases are scattered in the environment as follows: the environment is divided into K ∈ N equal parts and one communication base is randomly allocated to each part. Then, it is checked if all the sensors in each part are in communication radius of their related communication base or not. If not, another communication base is randomly assigned to that part and this process is repeated until all sensors are connected to the communication bases. The next step is to merge the communication bases to have the minimal number. For each two communication bases, their average Euclidean distance is calculated and a new communication base is placed on that point. It is checked whether all sensors, related to two previous communication bases, can be connected to the new one or not. If so, previous bases are replaced by the new one. This process continues for each two communication bases in the same part. Afterwards, the merging process will be applied on the all survived and newly obtained bases in different parts until the minimal number of bases is obtained. As a result, in each area, depending on the type and number of obstacles as well as the communication radius of the sensors, a number of sensors are connected to each communication base. The role of communication bases in energy conservation is critical. Utilizing the communication bases, messages including the information of targets can be exchanged among sensors in at most three hops. Unmanned Aerial Vehicles: As mentioned before, there are some entrance gates in the environment. Targets enter the environment from these gates and attempt to destroy one or more of the important sites as well in unknown time intervals. The task of wireless sensors is to predict the location, direction, and velocity of the targets with a probability. As previously mentioned, the main goal is to have all the sites at full health. Therefore, the Unmanned Aerial Vehicle (UAV) technology is used to protect them and destroy the targets [START_REF] Chandler | UAV Cooperative Control[END_REF]. For this purpose, a constant and independent UAV station is considered which is connected to all communication bases. Suppose the information including the presence of a target is delivered to the UAV station in time interval t. Depending on the information sent, this station orders to one of UAVs inside to fly over the environment. Since all the information received by the UAV station is approximate and probabilistic, the UAV flying will be also probabilistic. The UAV flies with the maximum speed until it reaches near the target. Reaching near the target, it decreases its velocity and integrates its spatial information with the information sent to the UAV station before and the new information obtained during the flight. Since the target's predicted position is probabilistic and also, there are many sensors around the target, the destruction should be done with a high accuracy. It is assumed that the probability of destroying the target after striking the missile is one. This assumption preserves the sensors from the damage.

The Distributed Cellular Algorithm

Initializing all the described procedures/elements is known as a pre-process for the main tracking algorithm. As the main parts of the proposed distributed cellular algorithm, two novel concepts named boundary sensors and sensors cell have been announced.

Boundary Sensors

As pictured in Fig. 2(a), a number of sensors are close to the entrance gates. Therefore, there is a unique set of boundary sensors in front of each entrance gate. The boundary sensors on the one hand are those which include the relevant gate in their sensing radius. Apparently, to be able to predict the direction of targets, these boundary sensors should always be awake and sense the environment, as there is no information of when and which gate the targets will enter the environment from. On the other hand, keeping all these boundary sensors awake is time and energy consuming. In this paper, it is assumed that the boundary sensors wake up periodically in specific time intervals and start sensing as long as a sensor detects a target. This approach will be a great help to reduce the energy consumption and eventually increase the reliability level of protecting the important bases. Fig. 5(a) is an example of boundary sensors in the environment. To be able to firstly, adapt this periodically sleeping/awaking strategy for boundary sensors and then, track the targets with maximum efficiency and accuracy, a concept called sensors' cell is proposed.

Sensors' Cell

In this section, the concept of the cell and its applicability in heightening the target tracking accuracy is introduced. As an important note, the targets are unaware of the Each target moves into the environment from one of the gates. Once it appears in the sensing radius of at least one boundary sensor, this sensor will sense and predict the position, velocity, and direction of the target. Obviously, the target will continue to move and leave the sensor's sensing radius. Therefore, firstly, to predict the position of the target with higher probability, and secondly, not to lose the target, a preliminary cell among corresponding boundary sensors is created. Fig. 5 indicates how the sensors' cell is constructed.

According to Fig. 5(b), sensors, which have a gate in their sensing radius, are connected to each other via their communication radius and they create a complete graph. This graph is represented by G = (V, E). V includes the set of boundary sensors and E is the set of edges which exhibits the communications between these sensors. To reduce the communication volume and have a planar graph, the convex hull of the vertices is created. The convex hull of a set of points in the plane is the smallest convex polygon which includes all points within or on the perimeter [START_REF] De Berg | Computational Geometry[END_REF]. Fig. 5(c) depicts the convex hull of each set of boundary sensors. The target will certainly enter inside one of the preliminary created cells. Additionally, this target will definitely get closer to one of the sensors than the others. In other words, as mentioned before, in the first step, all boundary sensors wake up and sleep periodically until a target is recognized. As soon as a target is detected, the sensor, which is closest to the target, continues its work (sensing) and the others wait and remain in awake mode. In later steps, the target definitely continues its way and goes out of the sensor's sensing radius. As a result, the sensor sends the collected information to other sensors of the cell, and then, next sensor begins to sense based on the received information, and other sensors stay in awake mode. The information tables are updated in each time interval and target path is obtained. Sensors inside the cell utilize the probabilistic direction and inform the other sensors of the neighborhood cell. These neighbor sensors create the new cell to continue tracking in that cell. Next, the procedure of creating the new cells based on the exchanged information will be described in detail.

Creating the New Cells: Considering the fact that the targets move continuously and leave the cells consecutively, the preliminary-created boundary cells should be updated in order that the targets are not lost and the tracking procedure continues. In each time interval, the active sensor in the current cell broadcasts an alert message to the sensors in the same cell and also its neighbor sensors which are not in the same cell. Fig. 6 indicates how to predict and create a new cell. According to target's estimated position and direction, sensors update their table and compare the prior position with the new estimated position and decide on which sensor the target is going to get close to. Utilizing the target's estimated direction, the new cell is created and this process is repeated until the target is tracked and destroyed.

Tracking

Efficiently tracking the targets is an indispensable task in the target tracking problem. To track a target appropriately, it is essential to estimate its position and direction to create the next cell before the target leaves the current cell. Considering the fact that the target's movement is non-linear, it is difficult to estimate its exact position. Hence, an efficient method, which is able to cope with target's non-deterministic movement, should be utilized. Here, the concept of a well-known probabilistic tracking method, i.e. particle filter (PF), is employed [START_REF] Arulampalam | A Tutorial on Particle Filters for Online Nonlinear/non-Gaussian Bayesian Tracking[END_REF]. PF is a sequential Monte Carlo method based on representing the probability densities. This method uses an online posterior density estimation algorithm that predicts the posterior density of the state-space by directly implementing the Bayesian recursion concept. Nonetheless, the pure PF method cannot overcome the degeneracy of particles and accumulation of estimation errors. Therefore, an improved version of this algorithm using Radial-Basis Function Network (RBFN) has been engaged with the proposed distributed cellular algorithm [START_REF] Wang | An Improved Particle Filter for Target Tracking in Sensor Systems[END_REF].

This improved version of PF employs the RBFN in the sampling step for dynamically constructing the process model from observing and updating the value of each par-ticle. With the sampling step in RBFN, the PF can give an accurate proposal distribution and maintain the convergence of a sensor system. RBFN is a three-layer feed-forward neural network which is embedded with several radial-basis functions. This network is characterized by an input layer, a single layer of non-linear processing neurons, and an output layer. The reason for applying RBFN is to dynamically construct the process model from observations and update the value of each particle. With the RBFN sampling step, this approach provides an accurate proposal distribution of a sensor. Fig. 7 demonstrates the changes in target tracking accuracy compared to the changes in the number of particles used in the PF. At each step of the algorithm, in which a cell is created and sensing and tracking operations are done inside that cell, all the information is circulated among communication bases. The reason is that since only communication bases are in contact with the UAV station, therefore, it is essential to send all information regarding the direction and velocity of all the targets to the UAV station. For instance, suppose tracking is being done in a cell and simultaneously a UAV flying over the environment is tracking the target. If another target enters the environment at that moment, a new instance of the tracking algorithm will be created. Hence, the communication bases should send the new information to UAVs station in order that new UAVs fly or the flying UAVs change their directions. In this paper, for simplicity, no constraint has been considered on the number of UAVs. Moreover, considering the communication protocols between all elements of the assumed problem is out of scope of this paper. As a final assumption for the tracking procedure, the targets are unable to destroy the UAVs.

Experimental Setup and Analysis

In this section, first, the experimental setup will be described, and subsequently, the experimental results will be clarified precisely.

Experimental Setup

The system model and proposed algorithms have been implemented with Matlab R2014b on a personal computer with an Intel Core i7 2.0 GHz CPU and 6 GB of memory. The dimension for the assumed environment has been done in the context of M = 600 m and N = 400 m. It has been assumed that there are G = 4 entrance gates (one in each side). Moreover, there are O = 8 two-dimensional planar and solid obstacles, and T = 10 important sites in the environment. The sensing and communication radius of all the sensors have been considered as R s = 35 meters and R c = 70 meters, respectively. Considering the communication radius twice the sensing radius represents the capability of the direct communication for the sensors which have a common sensing area. As the setting of evolutionary algorithm applied for sensors placement, the tournament selection operator (tournament size=2), linear crossover (probability=80%) and exchange mutation (probability=20%) have been set. Moreover, the generation and population sizes are equal to 400 and 100, respectively. The selected numbers here are the average of various tests. Therefore, S = 83 sensors are dispensed in the environment followed by B = 11 communication bases with respect to the assumed setting.

According to the intended purposes of this paper, the concept and characteristics of a sensor, namely Intel Imote2 [START_REF] Nachman | Imote2: Serious Computation at the Edge[END_REF] has been utilized. The reason for using this sensor is its capabilities in facilitating the tracking accuracy and quality. This sensor includes a CPU, a number of internal memories, and a communication device which uses the Zigbee radio communication to communicate between each two sensors. In addition, the sensor has a tiny OS that allows us to mount various algorithms on it. Followings are the corresponding amounts of energy consumed by each sensor for performing actions mentioned in Table 1.

-Action 1: 1 mJ/s -Action 2: 13.4 mJ/s -Action 3: 2.6 mJ/s -Action 4: 11.5 mJ/s -Action 5: 16.7 mJ/s -Action 6: 9.6 mJ/s -Action 7: 0.5 mJ/s Obviously, the energy consumption of sensors in off mode is zero. The default amount of energy for each sensor in the environment is equal to 1800 mJ. The amounts mentioned for Types 6 and 7 are based on the same simulation running for the proposed algorithm. In this paper, to simplify the model, the energy of communication bases has been considered infinite. Trying to dispense the communication bases in the environment, the environment is divided into K = 4 parts. Furthermore, it has been assumed that the initial places of UAVs are as a point at one corner of the environment and their energy is infinite. Another point about the UAVs is their velocity. In this study, the velocity for each UAV has been selected from the range [START_REF] Zhuang | Minimizing Energy Consumption With Probabilistic Distance Models in Wireless Sensor Networks[END_REF][START_REF] De Berg | Computational Geometry[END_REF] m/s. Now, selecting the proper velocity depends upon the received information including the target's direction and velocity. However, in general, UAVs upon receiving the information, make their attempt to fly as quickly as possible to reach near the target, and then, reduce their velocity gradually in order to track the target more accurately.

The sensors and UAVs are not aware of when and where the targets enter the environment from. In other words, the targets know how many they are, in which time interval they enter the environment, which gate they enter the environment from, and which important base they will go through. Furthermore, it is assumed that the total number of targets is 25 in the simulation framework. However, the moving direction is randomly determined in each time interval. Finally, the targets select their desired velocity from the range [0.5, 10] m/s.

Experimental Results

There are some studies in the field of centralized target tracking in which all sensors are usually awake and sense the environment [START_REF] Jiang | Probability-Based Prediction and Sleep Scheduling for Energy-Efficient Target Tracking in Sensor Networks[END_REF][START_REF] Zhang | DCTC: Dynamic Convoy Tree-Based Collaboration for Target Tracking in Sensor Networks[END_REF][START_REF] Cao | On Collaborative Tracking of a Target Group Using Binary Proximity Sensors[END_REF]. The aim of those studies is to firstly update the position of targets continuously and then transfer it among all the sensors. The updating cost in centralized mode depends on the network size involving sending/receiving many messages within the network. In contrast to this approach, the proposed distributed algorithm needs to send and receive less messages. There is no need to exchange the information among all the sensors and even the neighbors since the targets choose a dynamic path for moving. This issue is also helpful for energy efficiency of sensors. In other words, since the prediction and tracking in each step is performed inside a cell, other sensors and even its neighbors are able to go to sleep mode which will also reduce the energy consumption. Thus, each sensor only needs to communicate with its neighbors inside the cell or receive information from them.

Suppose a sensor goes out of work for any reason, the algorithm has been designed in a way that it does not have much effect on tracking quality and associated operations. In other words, deleting a sensor from the network as well as inside the cell, one of its neighboring sensors will be replaced and the cell still lasts. Merging two cells, when a sensor or sensors are lost, imposes a cost consisting O (n V ), where n V is the number of sensors used in the merge operation. In addition, the network will inflict an extra cost for updating the tracking information when the target moves (from one cell to another cell). As a result, the cost equals to the number of cells the targets pass. It is beneficial to note that the 1000 particles have been used in the simulations.

To assess the efficiency of the proposed algorithm, it has been compared with three well-known algorithms, VigilNet [START_REF] He | VigilNet: An Integrated Sensor Network System for Energy-Efficient Surveillance[END_REF], Forms [START_REF] Sarkar | Differential Forms for Target Tracking and Aggregate Queries in Distributed Networks[END_REF], and iTOA (improved Time Of Arrival measurement model) [START_REF] Xu | Target Tracking and Mobile Sensor Navigation in Wireless Sensor Networks[END_REF]. VigilNet algorithm has been chosen for the similar concepts used in this paper, its high speed in tracking operations, and specifications of the sensors. Forms has been selected due to using an idea close to the one proposed in this paper. Finally, the reason for choosing iTOA is its similarity of prediction and tracking parameters used followed by the concept of communication base. The results obtained are the average of ten experiments done on the whole system. Three parameters including total target tracking and destroying time (I T ), energy efficiency of the sensors (EE), and tracking quality (Q T ) have been utilized in order to evaluate the performance of the proposed algorithm. Regarding the tracking quality, it can be noted that the purpose is to investigate how successful the system is against all difficulties including target tracking error, low detection probability, sensor destruction, and eventually the obstacles, which exist in the environment.

Hereinafter, the term DEECTA (Distributed Energy-Efficient Cellular Target tracking Algorithm) will refer to the proposed algorithm. In the paper, Fig. 8(a) represents the total target tracking and destroying time compared with other algorithms in terms of variations in the target tracking velocity. The calculation of IT will continue since the simulation has started until all targets are destroyed. The numbers obtained for IT in Fig. 8(a) are the average times consumed for tracking all targets. As a final note for this evaluation, in spite of the fact that assigning a greater velocity to the UAV will decrease the efficiency and even the destruction of target will be harder, however, the proposed algorithm works better than other algorithms even with a large number of target miss- In this assessment, the performance of algorithm is evaluated with changes in UAVs velocity. It can be seen that when the UAVs' velocity fluctuates between 6 to 12 m/s, the amount of energy consumed by sensors is reduced. Furthermore, Fig. 8(b) indicates that the algorithm Forms utilizes more energy than other algorithms. The reason is the large number of sensors used in its tracking operation. The reason for superiority of the DEECTA algorithm is that at each step of the algorithm, a few number of sensors, that is, only one sensor in each cell, is in communication with others. In Fig. 9(a) , the impact of the sensor density on the performance of the algorithm in relation to energy consumption has been discussed with different sensing radius. If the sensing radius is large, then the density will be low. This in turn could increase the error rate of localization and target tracking since many small number of sensors are involved in tracking. Now, comparing the impact of density with the amount of energy consumed, the authors come to this conclusion that when the sensors density increases, energy consumption declines. The reason is the limited participation of sensors due to the idea of cell for tracking the targets. However, the proposed algorithm saves large quantity of sensor energies because of the reduction in the amount of communication. involved. Finally, Fig. 9(b) represents the tracking quality compared to the different UAVs' speeds. In this case, the number of sensors involved in tracking will be important. Note that when the number of sensors is too high or low, tracking will be done slower than if average number of sensors are involved. In addition, changing the sensing radius of the sensors, as depicted in Fig. 9(a), influences the number of utilized sensors and consecutively the number of communication bases.

As a final point, the effectiveness of using communication bases on energy efficiency and tracking accuracy will be noted. Using these bases, energy efficiency of sensors increases by 48% compared to the case in which communication bases are not used. Also, the communication volume decreases by 79% using these bases. If communication bases did not exist in the environment, sensors had to communicate more in order to transmit their information to the UAV base which would also result in delay in tracking, rapid reduction of sensors' energy and increase in information loss. As a

Conclusion and Future Work

In this paper, the problem of multi-sensor multi-target tracking using dynamic clustering and sleep scheduling in wireless sensor networks is investigated. As the main contribution, a distributed cellular cluster formation algorithm has been proposed to tackle with the problem of tracking and destroying the targets to protect the important sites in the environment. As a pre-process for the proposed algorithm, a multi-objective evolutionary algorithm has been used for the sensor placement with two conflicting objectives including minimizing the number of sensors and maximizing the area coverage. Next, a new concept named cell is proposed to minimize the number of awake sensors participated in tracking targets. Applying this method, only the sensors in the cell are turned on periodically and there is no need to wake up all sensors. Moreover, some communication bases have been used in order to decrease the sensors' amount of communications and also number of messaging hops. Maximizing the tracking quality motivates the authors to use an improved version of the particle filter algorithm based on radial basis function networks. Finally, in order to be able to destroy the targets moving in the environment, a UAV station including a set of UAVs, has been lodged in the environment. The results include reduction in communication volume, energy consumption, and also improvement of reliability and tracking quality.

As multiple works which can be done in the future, the possibility of applying the learning approaches for UAVs and sensors, considering a finite energy source for communication bases in order to propose a better and robust cluster formation and tracking mechanisms, improved target tracking procedure, and sleep scheduling algorithm for sensors can be mentioned.

Fig. 1 :

 1 Fig. 1: The status of sensors with different distances from each other

Fig. 2 :

 2 Fig. 2: An example of localizing the sensors through the NSGA-II

Fig. 3 : 6

 36 Fig. 3: The state diagram of energy variations for each sensor

Fig. 4 :

 4 Fig. 4: An example of placing the communication bases in the environment

Fig. 5 :

 5 Fig. 5: Forming the preliminary cell among boundary sensors

Fig. 6 :

 6 Fig. 6: Constructing a new cell based on the predicted trajectory of the target

Fig. 7 :

 7 Fig. 7: Prediction accuracy versus number of particles

Fig. 8 :

 8 Fig. 8: Evaluations of the algorithm considering the changes in the tracking speed of the UAVs

Fig. 9 :

 9 Fig. 9: Energy efficiency of the sensors in a basis of increasing the number of sensors combined with the tracking quality evaluation

Table 1 :

 1 Actions and corresponding energy consumption of the sensors