Georgios Siamantas

Konstantinos Gatsis

Antony Tzes
email: tzes@ece.upatras.gr

Mobile Robot-Assisted Cellular Environment Coverage

Keywords: Robotics, Control, Planning

The robotic coverage problem of a known rectangular cellular environment with obstacles is considered in this article. The robot can move only at directions parallel and perpendicular to the sides of the rectangle and can cover one cell at each time unit. A suboptimal minimum-time complete area coverage path-planning algorithm is presented. This algorithm minimizes a distance cost metric related to its current position and the quadtree-based decomposed blocks of the unexplored space. The efficiency of the suggested algorithm is shown through simulation studies and the implementation of a non-linear controller designed for a two-wheeled robot to achieve tracking of reference trajectories. For the implementation of the controller, the localization of the robotic vehicle was necessary and it was achieved via image processing.

Introduction

Robotic path-planning algorithms tackle the problem of finding a path in an environment that enables a robot to visit or sense with its sensors as much as or all if possible of this environment. Many applications of such algorithms appear in the literature [START_REF] Choset | Coverage for robotics-A survey of recent results[END_REF]. In the present work we consider a robot placed inside a cellular environment with obstacles. We seek to find a path that minimizes the time that takes the robot to visit each cell of this unexplored environment. The entry and exit points of this path can be different cells. The optimal solution for the suggested l 1 -distance cost metric rectilinear environment is NP-hard [START_REF] Deng | How to learn an unknown environment I: the rectilinear case[END_REF].

Similar works were presented by Zelinsky [START_REF] Zelinsky | Planning Paths of Complete Coverage of an Unstructured Environment by a Mobile Robot[END_REF] and Gabriely and Rimon [START_REF] Gabriely | Spanning-tree based coverage of continuous areas by a mobile robot[END_REF]. Zelinsky used the conventional wavefront algorithm to determine a coverage path (distance transform coverage-DTC). In DTC the starting and ending (goal) robot positions are specified in order to create an expanding field of increasing distance values from the end position extending to the whole grid. The robot then finds a cover path from start to goal performing a "pseudo-gradient ascent", on the distance field. Gabriely and Rimon describe a Spanning Tree Coverage (STC) algorithm. The STC works not on the original cellular environment but on a coarsen grid which is a tool-based approximation of the work-area and achieves a covering path in linear time O(N ), where N is the number of cells comprising the approximate area. The DTC method is very well suited for applications where we want to create a covering path in O(n) time (n the grid area) and finish in a specific cell in the grid. But it does not guarantee that this path will be as close to the o as possible. In this paper the proposed coverage method produces close to optimal paths most of the time, overcomes the need of specifying a goal position as required in the DTC algorithm and works directly in the original cellular environment grid in contrast to the STC algorithm.

The paper is organized as follows. The model of the environment, the robot motion and associated metrics and methods are described in Section 2. The proposed algorithm is provided in Section 3 followed by simulation results and robot experiments in Section 4, while conclusions are drawn in the last Section 5.

Problem statement

Preliminaries, Concepts

We consider a cellular environment to be a rectangular grid of n x n dimensions where n = 2 k , k > 0. The grid cells can be empty (free) and thus possible for the robot to visit them or occupied by an obstacle and thus inaccessible to the robot. The topology of the grid is dictated by the number and positions of the obstacles inside the grid. Each grid cell has coordinates of the form: (i, j), where 1≤ i ≤ n and 1≤ j ≤ n. The robot can move under the Manhattan or l 1 -metric motion model i.e. parallel or perpendicular to the sides of the rectangle and only one grid cell at each time unit, thus the set of actions (moves) that the robot can make is U = {(0, 1), (1, 0), (0, -1), (-1, 0)} where the robot's position is changed at the kth time instant to X k+1 (i, j) = X k (i, j) + U . Each empty cell has from 1 to 4 neighbor free cells depending on adjacent obstacles and on whether is located at a corner (1-2 neighbors), an edge (1-3 neighbors) or inside the grid (1-4 neighbors). Each empty cell is accessible from any other empty cell in the grid, i.e. obstacles do not form closed regions of empty space inside the grid. The following notation is used in the remaining of this article: C the set of grid cells, C f ree the set of free grid cells, C obs the set of cells occupied by obstacles, C f eas the set of cells that the robot can visit from the current position, C vis the set of already visited cells, C unvis the set of remaining unvisited cells, X i,j the grid position at coordinates (i, j), X cur the robot's current position coordinates, X next the robot's next move position coordinates, X I the robot's initial position coordinates.

Under the above environment formulation the minimum-time full area coverage path must intersect itself as less times as possible if not at all. One known way [START_REF] Choset | Coverage of known spaces: The boustrophedon cellular decomposition[END_REF] to achieve this (minimizing intersections) is to execute a serpentine like motion consisting of straight moves followed by 180 degree U-turns, the so called boustrophedon (ox turning-like) motion and try to avoid visiting already visited space. In order to achieve as much as possible boustrophedon motion in a known topology environment is to ensure that at each path step the sum of the distances from the next robot position to as many as possible landmarks in the environment is maximized. An equivalent to this maximization criterion is to minimize the sum of the distances from all the other possible move positions except from the current pick position, e.g. if the robot can move to four adjacent cells then pick one possible move at a time and measure the distance of the other three possible move cells to the landmarks and then move to the cell with the minimum sum of these distances. We used the minimization criterion instead of the maximization because in the l 1 metric cellular environment, the simulation studies indicated a reduced sensitivity when dealing with moves that are equidistant from the landmarks. In our environment formulation all the possible robot positions are countable, so it is possible to sum the distances from the next position to all free cell positions in the grid. Although this is possible it is time consuming for large grid sizes. An equivalent criterion is to decompose free cell space to blocks and then sum the distances from the next position to some cell inside these blocks. The distances should be weighted according to the size of each block to ensure that the mean block distance is approximately the same with that of summing the distances of all the cells inside the block. The concept is shown in Figure 1(a). If the current robot position is in cell C and we choose to move to cell N the two other possible moves are to cells R. We compute the sum of all distances q 1 , q 2 etc. from cells R to the decomposed free space blocks and repeat the same computation for each of the 3 possible moves. The criterion to move to the next cell is decided by the one that makes the sum of these distances minimum.

Quad-tree decomposition

An n x n grid matrix G is defined such as:

G(i, j) = 1 if X i,j ∈ C obs or X i,j ∈ C vis 0 if X i,j ∈ C unvis (1)
where the grid cells that are considered obstacles and the already visited cells (soft obstacles) take value 1 and all other unvisited cells take value 0. In this paper the Quad Tree Decomposition (QTD) algorithm is used for decomposing the environment grid matrix into rectangular blocks of similar value cells. The QTD is a well established method [START_REF] Samet | Region Representation: Quadtrees from boundary codes[END_REF] of decomposing an image into blocks that meet some criterion of homogeneity and has been used extensively in robotics [START_REF] Kambhampati | Multiresolution Path Planning for Mobile Robots[END_REF]. The motivation for decomposing the grid matrix into blocks is to reduce the number of cell distance calculations by some constant factor. The quad tree method decomposes the grid matrix into two parts. The first part is a tree structure that shows the location and size of each homogeneous region. The second part specifies the cell value at each leaf node of the tree, representing if it is an obstacle, visited space or unvisited space.

Quad tree decomposition is based on the successive subdivision of the grid matrix into four equal-sized quadrants. If the matrix does not consist entirely of 1s or entirely of 0s, then it is subdivided into quadrants, sub-quadrants, and so on, until blocks are obtained that consist entirely of 1s or entirely of 0s; i.e. each block entirely contains cells that belong to C obs and C vis or entirely contains cells that belong to C unvis . The grid quad tree can be characterized as a variable resolution data structure. QTD is an O(n) complexity process, where n is the grid area [START_REF] Kambhampati | Multiresolution Path Planning for Mobile Robots[END_REF]. The worst case for the QTD in terms of distance calculations occurs when the topology of the grid corresponds to a checkerboard pattern. The amount of distance calculations required is a function of the resolution (i.e., the number of levels in the quad tree), the number of obstacles and visited cells, and their relative positioning inside the grid. Thus, when calculating the sum of distances instead of dealing with all the free grid cells we have to deal only with the nodes of the QTD, a substantial reduction.

Shortest Path Computation

In order to measure distances in the l 1 -metric grid we used an algorithm proposed by [START_REF] Lee | An Algorithm for path connections and its applications[END_REF] that it is widely used in routing two terminal nets in circuit board design problems. This algorithm has O(n) time and space complexity, where n is the area of the grid matrix and guarantees of finding an optimal solution if one exist. The used algorithm consists of two phases [START_REF] Sherwani | Algorithms for VLSI Physical Design Automation[END_REF]. The propagation phase is based on breadth-first search. The search resembles a wave propagating from the starting cell. The start cell is labeled '0' and the wave front propagates to all the unblocked vertices adjacent to it. Every unblocked cell adjacent to the start cell is marked with a label '1' corresponding to its l 1 distance from start. Then, every unblocked cell adjacent to the cells with distance 1 is marked with a label '2' because its distance from the start cell is 2, and so on. This process continues until the target cell is reached or no further expansion of the wave can be carried out. An example of the algorithm is shown in Figure 1(b). Due to the breadth-first nature of the search, this algorithm is guaranteed to find a path between a source and a target cell, if one exists. In addition, it is guaranteed to be the shortest path between these cells [START_REF] Sherwani | Algorithms for VLSI Physical Design Automation[END_REF]. The second phase is the retrace phase where from the target cell we follow back the cell labels in descending distance order till we reach the starting cell. In our formulation for the sum of distances computation the propagation phase is computed once in each robot's move. The retrace phase is not needed because we only need to find the minimum distance path between the starting and ending cells if such a path exists.

Environment Coverage

Let B the set of all QTD blocks, B f ree the set of QTD blocks of cells with value 0, B occ the set of QTD blocks of cells with value 1, B i,j the position of a block as the upper left cell position X i,j of this block, A B i,j the area of a block B i,j given from the QTD. d m,n k,l is the l 1 -metric distance from X k,l to X m,n and is returned by the aforementioned in Subsection 2.3 algorithm. If there is no path connecting X k,l to X m,n because of the presence of obstacles and/or cells that are already visited then the implemented algorithm returns a very large value to bias the next move properly.

In order to choose which move to select the following quantity is defined:

D m,n i,j = ∀X k,l ∈{C f eas \{Xi,j }} d m,n k,l . (2)
This quantity reflects the definition of the sum of the distances from some cell X m,n to all the other possible next moves except the one for which we compute the sum.

Similarly, let:

∆ i,j = ∀Bm,n∈B f ree A B m,n D m,n i,j , (3)
where ∆ i,j is the sum of the previously defined D distances to all the free block positions scaled by each block's area as explained in section II.

If the set M is defined as:

{ ∆ i,j | X i,j ∈ C f eas }, (4)
then the position coordinates of the robot's next move in order to fulfill the minimization criterion of the proposed algorithm can be set as:

X next = arg inf Xi,j M. (5)
In the algorithm Grid Cover shown below function QuadTreeDecomp() performs quad tree decomposition on the given grid. Function FindUnvisitedCell-Blocks() finds in the quad tree structure the QTD blocks of cells with value 0 i.e. that correspond to unvisited cells. Function FindAdjacentUnvisitedCells() finds the adjacent to X cur unvisited cells that are not obstacles. Function GridIs-Covered() checks if the grid is covered by the robot. Function FindNearestUn-visitedCellTo() finds the nearest to X cur unvisited cell (even if it is necessary to pass over already visited cells) thus implementing a backtracking mechanism that guarantees the grid's full coverage by the robot. Function FindMinDis-tancePathTo() finds a minimum distance path to that nearest cell using the suggested shortest path algorithm.

Algorithm Grid Cover (G, X I ) Input: A grid matrix G; the starting cell position X I . Output: An array Path with the covering path cell positions. 1.

X cur ← X I 2. S ← QuadTreeDecomp(G) 3. B f ree ← FindUnvisitedCellBlocks(S) 4. C f eas ← FindAdjacentUnvisitedCells(X cur ) 5. if sizeof(C f eas ) = 1 then X N ext ← C f eas ; goto 8. 6. if sizeof(C f eas ) > 1 then
With C f eas and B f ree compute equations ( 2)-( 4) Set X N ext according to equation ( 5); goto 8.

7. if sizeof (C f eas ) = 0 then if GridIsCovered() = True then return Path else X near ← FindNearestUnvisitedCellTo(X cur ) Path ← Path + FindMinDistancePathTo(X near ) X N ext ← X near 8. Update Path, G, X cur with X N ext ; goto 2.

Simulations and Experiments

In order to evaluate the performance of the algorithm we define as M N P the minimum number of path steps, as AN P the algorithm number of path steps and as a performance index P I opt = M N P/AN P . The numerator refers to the number of robot moves needed to cover optimally the whole grid in minimum time and the denominator refers to the number of moves proposed by the algorithm for full area coverage. So the grid will be covered optimally if P I opt =1. Several tests were performed with different grid areas, obstacle topologies and robot start positions and the algorithm performed acceptably with P I opt > 0.87 most of the time. The proposed algorithm shows O(n 2 ) time complexity where n corresponds to the grid area. This is because wavefront propagation algorithms have an O(n) time complexity and that, without taking QTDs into account, it takes at most (n 2 + n)/2 distance propagation related operations to cover the whole area. DTC computed for all possible end positions in the grid (full DTC) has also O(n 2 ) time complexity.

In Figure 2 we see a comparison of the performance of the proposed algorithm with full DTC. In an 8x8 grid 100 simulation experiments were performed in MATLAB with 15 uniformly randomly distributed obstacles each time. In Figure 2(a) we see that the proposed algorithm (solid line) creates cover paths that have approximately the same number of steps as the minimum-step paths and considerably faster paths than the maximum-step ones created by the full DTC algorithm. Also in Figure 2(b) we see the execution times in the same experiment. The proposed algorithm (dashed line) was almost 7 times faster than full DTC. In Figure 3 we see the Amigo differential drive two-wheeled robot that was used to test the proposed algorithm in real time. In order to do this a nonlinear controller was implemented, based on [START_REF] Klancar | Mobile Robot Control on a Reference Path[END_REF], in order to achieve tracking of the required by the algorithm reference trajectories. For the implementation of this controller, the localization of the robotic vehicle was necessary and it was achieved via image processing during the experiments. The robot was recognized in images taken by a surveillance camera Figure 3(b), and after appropriate geometric analysis, it was possible to estimate the position and the orientation of the robot. In order to make it easier to recognize the robot's position and orientation inside the image, the robot's external surface was covered by white material and three small black spheres were placed on it, as shown in Figure 3(a). As shown in Figure 3(c) the robot followed the reference path computed by the algorithm with acceptable accuracy.

Conclusions

In this work a suboptimum complete area coverage algorithm is presented for rectangular l 1 -metric grid environments. The start and finish positions may not be the same and the grid may have an arbitrary number of obstacles in arbitrary positions. The algorithm showed good performance in various grid topologies with acceptable time complexity taking into account its simplicity and hardness of the problem. The algorithm works on the original cellular environment and not on a coarsen approximation grid like the STC algorithm. It also guarantees the grid's full coverage by the robot because of its implemented backtracking mechanism. For comparison we performed various experiments with random topologies and saw that the proposed algorithm produced equally fast cover paths with the minimum-step paths produced by the full DTC algorithm but with considerably less computation time. The full DTC algorithm is an extension of the original DTC algorithm that we derived in order to compare it with the proposed algorithm. In the original DTC the goal position must be specified and produces cover paths with varying step sizes depended on this position. Finally the performance of the algorithm was validated in real time experiments with a differential drive two-wheeled robot moving according to the algorithm output in a laboratory test environment.

Fig. 1 .

 1 Fig. 1. (a) Sum of distances computation (b) The propagation and retrace phases.

Fig. 2 .

 2 Fig. 2. Simulation results. The gray region represents the minimum and maximum path steps generated by full DTC.

Fig. 3 .

 3 Fig. 3. (a) Robot in test environment (b) Camera (c) Reference and real robot paths.


Stamos Katsigiannis
email: stamos@di.uoa.gr

Eystratios G Keramidas
email: e.keramidas@di.uoa.gr

Dimitris Maroulis

Contourlet Transform for Texture Representation of Ultrasound Thyroid Images

Keywords: contourlet transform, ultrasound images, feature extraction, thyroid, feature selection

Texture representation of ultrasound (US) images is currently considered a major issue in medical image analysis. This paper investigates the texture representation of thyroid tissue via features based on the Contourlet Transform (CT) using different types of filter banks. A variety of statistical texture features based on CT coefficients, have been considered through a selection schema. The Sequential Float Feature Selection (SFFS) algorithm with a k-NN classifier has been applied in order to investigate the most representative set of CT features. For the experimental evaluation a set of normal and nodular ultrasound thyroid textures have been utilized. The maximum classification accuracy was 93%, showing that CT based texture features can be successfully applied for the representation of different types of texture in US thyroid images.

Introduction

Modern medical ultrasonography (US) presents a unique set of advantages including real-time data acquisition, low cost, absence of any side effects and high resolution imaging. Thus, US has become the most common imaging modality for certain types of medical examinations, including that of the thyroid gland. Currently a challenge for the image analysis community is the application of automated or semi-automated computational methods on US images, for the computer aided diagnosis (CAD) of different types of diseases. A key issue for any CAD approach remains the texture representation method utilized.

Texture representation is a fundamental issue in image analysis and computer vision. It has been extensively investigated in the literature for more than three decades. Numerous approaches have been proposed dealing with textural features extraction which can be divided into four main categories [START_REF] Tuceryan | Texture Analysis, The Handbook of Pattern Recognition and Computer Vision[END_REF] statistical, signal processing, model-based and geometrical. Signal processing approaches have drawn much attention, resulting in the proposal of a variety of texture representation methods, including the power spectral methods using the Fourier spectrum, the Discrete Cosine Transform (DCT), the Discrete Hartley Transform (DHT) [START_REF] Manian | Texture Classification Using Logical Operators[END_REF], and more recently, the Gabor filters [START_REF] Pichler | A comparison of texture feature extraction using adaptive Gabor filtering, pyramidal and tree structured wavelet transforms[END_REF], the Haar [START_REF] Lonnestad | A new set of texture features based on the Haar transform[END_REF] basis functions, the Discrete Wavelet Transform (DWT) and the Contourlet Transform (CT) [START_REF] Do | Contourlets: A Directional Multiresolution Image Representation[END_REF] [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF].

Although the Discrete Wavelet Transform has been successfully applied for a wide range of image analysis problems, for two dimensions it tends to ignore the smoothness along contours [START_REF] Do | Contourlets: A Directional Multiresolution Image Representation[END_REF]. In addition, the DWT provides only limited directional information which is an important aspect of multidimensional signals [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF]. These limitations have been partially addressed by the CT which can efficiently approximate a smooth contour at multiple resolutions. Additionally in the frequency domain, the CT offers a multiscale and directional decomposition, providing anisotropy and directionality, features missing from the DWT [START_REF] Do | Contourlets: A Directional Multiresolution Image Representation[END_REF][6] [START_REF] Po | Directional multiscale modelling of images using the contourlet transform[END_REF]. The CT has been successfully used in a variety of texture analysis applications, including SAR and natural image classification [START_REF] Liu | Minimum Distance Texture Classification of SAR Images in Contourlet Domain[END_REF], content-based image retrieval [START_REF] Srinivasa Rao Ch | Content Based Image Retrieval using Contourlet Transform[END_REF], image denoising [START_REF] Tsakanikas | Improving 2-DE gel image denoising using contourlets[END_REF], despeckling of images [START_REF] Varshney | Despeckling Synthetic Aperture Radar Imagery using the Contourlet Transform, Application of signal processing[END_REF], image compression [START_REF] Karras | Image Compression Using the Wavelet Transform on Textural Regions of Interest[END_REF], etc.

The aim of this study is to investigate the performance of the Contourlet Transform (CT) for the representation of medical ultrasound (US) textures of the thyroid gland. A set of statistical features calculated from the CT coefficients are evaluated through a supervised classification schema on real thyroid US images. Additionally, a feature selection phase has been applied through the Sequential Float Feature Selection (SFFS) algorithm, for the extraction of the most representative set of CT features.

The rest of this paper is organized in three sections. Section 2 presents the proposed methodology including the CT and the above mentioned feature extraction schema. This methodology is evaluated through an experimental study on real thyroid ultrasound images presented in section 3, whereas conclusions are presented in section 4.

Methodology

The Contourlet Transform

The CT is a directional multiresolution image representation scheme proposed by Do and Vetterli [START_REF] Do | Contourlets: A Directional Multiresolution Image Representation[END_REF] which is effective in representing smooth contours in different directions of an image [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF], thus providing directionality and anisotropy [START_REF] Do | Contourlets: A Directional Multiresolution Image Representation[END_REF] [START_REF] Srinivasa Rao Ch | Content Based Image Retrieval using Contourlet Transform[END_REF]. The method utilizes a double filter bank, in order to obtain a sparse expansion of typical images containing smooth contours. In this filter bank, first the Laplacian Pyramid (LP) is used to detect the point discontinuities of the image and then a Directional Filter Bank (DFB) to link point discontinuities into linear structures. This scheme results in an image expansion that uses basic elements like contour segments and thus it is named CT. The separation of directional and multiscale decomposition stages provides a fast and flexible transform with computational complexity O(N) for Npixel images when using Finite Impulse Response (FIR) filters [6][8].

The LP, introduced in [START_REF] Burt | The Laplacian Pyramid as a Compact Image Code[END_REF], provides a way to obtain multiscale decomposition. In each decomposition level, it creates a downsampled lowpass version of the original image and a bandpass image. A coarse image with the lower frequencies and a more detailed image with the supplementary high frequencies containing the point discontinuities are obtained. This scheme can be iterated continuously in the lowpass image, restricted only from the size of the original image.

The next step consists of the DFB which is a 2-D directional filter bank proposed by Bamberger and Smith, that can be maximally decimated while achieving perfect reconstruction [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF] [START_REF] Bamberger | A filter bank for the directional decomposition of images: Theory and design[END_REF]. The original DFB is efficiently implemented via an l-level binary tree leading to 2 l subbands with wedge-shaped frequency partitioning. However, in the CT, a new method is applied [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF] which avoids modulating the input image and uses a simpler rule for expanding the decomposition tree [START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coefficients[END_REF]. The simplified DFB consists of two stages. First, a two-channel quincunx filter bank [START_REF] Vetterli | Multidimensional subband coding: Some theory and algorithms[END_REF] that divides the 2-D spectrum into vertical and horizontal directions, and then a shearing operator that just reorders the samples.

Bandpass images from the LP decomposition are fed into a DFB in order to obtain the directional information. This scheme can be iterated on the coarse image and the combined result is a double iterated filter bank, named contourlet filter bank [START_REF] Liu | Minimum Distance Texture Classification of SAR Images in Contourlet Domain[END_REF], which decomposes images into directional subbands at multiple scales.

Feature extraction

The set of statistical texture features in the contourlet domain evaluated in this study, have been proposed by Liu [8]. This set consists of the first order statistical measures that follow: i) Energy. The energy E jk of a subband image I jk from the CT decomposition is defined as

∑∑ = = ⋅ = jk jk N n M m jk jk jk jk n m I N M E 1 1 2 )] , ( [ 1 .
(

) 1
Where I jk is the subband image of the kth direction in the jth level. M jk is the row size and N jk the column size of the subband image I jk .

ii) Standard deviation. Provides the means to capture the scale of the diversity of the image. The standard deviation S jk of the subband image I jk is defined as

∑∑ = = - ⋅ = jk jk N n M m jk jk jk jk jk I n m I N M S 1 1 2 _ ) ) , ( ( 1 .
(

) 2
∑∑ = = ⋅ = jk jk N n M m jk jk jk jk n m I N M I 1 1 _ ) , ( 1
.
(

) 3
iii) Information entropy. For texture images, it represents the complexity of the texture information. The information entropy H jk of a subband image I jk is defined as

∑∑ = = ⋅ - = jk jk N n M m jk jk jk n m p n m p H 1 1 ) , ( log ) , ( . (4)
Using the above features, the feature vector for the subband image of the kth direction in the jth level is defined as f jk =(E jk , S jk , H jk ). If the CT is defined as J level and for the jth (j=1,…, J) level a K j bands DFB is applied, then a total number o K subband images is obtained, where

∑ = = J j j K K 1 . (6)
By combining and rearranging all the feature vectors of the K subband images, the feature vector F={ f i }, i=1,2,…,K of the input image is obtained and normalized as proposed in [START_REF] Srinivasa Rao Ch | Content Based Image Retrieval using Contourlet Transform[END_REF].

Feature Selection

Feature selection offers more than one significant advantages, including reduction of computational complexity, improved generalization ability and robustness against outliers. A widely adopted algorithm for feature selection is the sequential floating forward selection (SFFS) algorithm [START_REF] Pudil | Statistical approach to pattern recognition: Theory and practical solution by means of PREDITAS system[END_REF] which has been used in a broad range of applications [START_REF] Ververidis | Fast Sequential Floating Forward Selection applied to emotional speech features estimated on DES and SUSAS data collections[END_REF]. The idea behind the SFFS algorithm consists of consecutive forward selection steps, followed by a number of backward steps as long as the resulting subsets are better than the previously evaluated ones at the same level.

Starting from an initially empty set of features, at each forward step an additional feature for which the classification accuracy is maximized is selected. Respectively, at each backward step the maximum subset that results in an improved classification accuracy is being selected.

Experimental evaluation

Supervised classification experiments have been carried out to evaluate the performance of the proposed texture representation approach. The image set consists of real medical ultrasound images of the thyroid gland. A total of 72 thyroid ultrasound images were obtained from examinations performed on 43 patients, using a Philips HDI 5000 sonographic imaging system, with a spatial resolution of 470×470 pixels and amplitude resolution of 8 bits. During these examinations the parameters of the sonograph were kept the same. This set of ultrasound images includes hypogenic thyroid nodules, as it has been diagnosed by expert physicians. From each ultrasound image, an equal number of healthy and nodular sample blocks have been selected. The total number of non overlapping 32x32 pixel blocks (Fig. 2) resulting from this process was 200. The classification task was implemented by means of the non parametric and generally effective k-NN classification approach [24]. The distance measure used by the k-NN was the Euclidean and parameter k was {3, 5}. For all experiments conducted, the classification accuracy for all experiments was estimated by 10-fold cross validation [START_REF] Theodoridis | Pattern Recognition[END_REF].

Due to the small size of thyroid nodules, the sample blocks selected are also small (32x32 pixels) and they did not allow decomposition of more than three levels with the LP and six levels with the DFB. The first LP decomposition level supports up to six levels of DFB decomposition, decreased by one for every extra LP decomposition level. The filters applied for the LP were the filters Burt, 5-3 and 9-7 and for the DFB the 5-3 and 9-7 filters. All possible combinations of filters and decomposition levels were tested.

Detailed results about the maximum classification accuracy obtained using the k-NN classifier and different filter combinations are shown in Fig. 3. From this figure it can be noticed that the maximum classification accuracy is 77%. This resulted from the filter 9-7 for the LP and the DFB, with two levels of LP decomposition, decomposed into sixteen and two directional subbands respectively, from finer to coarser scale.

The application of the SFFS selection algorithm led to improved classification results for every experimental setup as shown in Fig. 4. In this case the maximum classification accuracy was 93% using only 36 out of 205 features and the confusion matrix is presented in Table 1. This accuracy has been obtained via 5-3 filters for the LP and the 9-7 filters for the DFB, with two levels of LP decomposition, and then decomposed into sixty four and four directional subbands respectively from finer to coarser scale.

This result shows that a significantly smaller subset of CT features is necessary and sufficient to describe effectively the thyroid US texture.

Conclusions

In this study a methodology for the texture representation of thyroid tissue in US images has been investigated, utilizing features based on the Contourlet Transform (CT) and various types of filter banks. The experimental evaluation through supervised classification on real US thyroid images led to promising results. Furthermore, through a feature selection phase the maximum classification accuracy reached 93% for a significantly smaller set of features. These results are considered to provide evidence for the effectiveness of CT texture representation of US thyroid images. As overall conclusion, it can be argued that the combinations of 5-3 and 9-7 filters resulted in better classification accuracy for the second level of LP decomposition. Future work could include US images of higher resolution and the evaluation of different types of statistical features.

Fig. 1 .

 1 Fig. 1. The contourlet filter bank. Images l i and h i are the lowpass and detailed image, respectively, obtained from the ith level of decomposition with the LP. In every level of LP decomposition, h i is fed into the DFB.

Fig. 2 .

 2 Fig. 2. Sample images from experimental dataset. (a) Normal and (b) nodular thyroid tissue.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Maximum classification accuracy of the k-NN classifier for each filter combination at each level of LP decomposition (first, second, third). The first filter of the filter combinations refers to the one used for the LP and the second refers to the one used for the DFB.

Table 1 .

 1 Confusion matrix for the setup which gave the best results (93% accuracy).

	Correct
		Nodular	Healthy
	Prediction
	Nodular	45,6 %	4.4 %
	Healthy	2.1 %	47.9 %

Acknowledgement

This work was partially supported by the National and Kapodistrian University of Athens Special Account of Research Grants. We would also like to thank EUROMEDICA S.A. Greece for the provision of part of the medical images.