Petros Fotios
email: petros@channelsight.com

Andrian Moroi
email: andrian@channelsight.com

George Bagropoulos
email: george@channelsight.com

Kieran Dun- Don
email: kieran@channelsight.com

Content Based Image Retrieval and its Application to Product Recognition

Keywords: EMD, CBIR, Product Recognition, Image Descriptor, 1

Product Recognition is a challenging problem in many practical applications. This paper presents a new approach for product recognition. By utilizing a set of crawlers our task is to extract informative content from web pages and automatically recognize products found on web pages. A set of images is extracted from each web page and then a new "content-based" image retrieval technique is performed to rank the images from our product catalog. The proposed content-based image retrieval technique utilizes the Empirical Mode Decomposition and processes the first extracted component of the source image. This component maintains the highest local spatial variations of the source image. An adaptive local-threshold technique is applied for the extraction of edges. A quantized and normalized histogram is created for the representation of images. Simulation results reveal that the proposed method is a promising tool for the challenge task of product recognition.

Introduction

Product Recognition (PR) is a challenging task. To date, an increasing number of internet users buy products from well-known e-commerce web sites such as Amazon, Newegg and Walmart. Moreover, the interest of manufacturers for product coverage in web sites has grown in recent years. Many companies invest in information technology in order to increase their sales by means of implementing efficient product search services. One of the most popular approaches is to offer to clients straightforward web-services. These services are dedicated to find reviews, specifications availability and prices of products from different retailers with the ultimate goal to purchase them.

Today the tremendously evolution in portable electronic devices like smartphones and tablets, provides a perfect environment where more and more users tend to shop online. Using these devices several web-services are implemented, having as goal to provide an attractive and alternative environment to proceed to shopping actions. Current product-image search systems on the market include Google Goggles [START_REF]Google goggles[END_REF], Vuforia [START_REF]Vuforia[END_REF] and Flow [START_REF]Flow[END_REF].

Product Image Retrieval (PIR) is an emerging field of content-based image retrieval (CBIR). Several algorithms and systems have been introduced in literature for this new research topic [START_REF] Jing | Pagerank for product image search[END_REF][START_REF] Tsai | Mobile product recognition[END_REF][START_REF] Zhang | Mobile product recognition with efficient Bag-of-Phrase visual search[END_REF][START_REF] Chai | Re-ranking using compression-based distance measure for content-based commercial product image retrieval[END_REF][START_REF] He | Mobile product search with bag of hash bits[END_REF][START_REF] Wang | When Amazon meets Google: Product visualization by exploring multiple Web sources[END_REF][START_REF] Eason | On certain integrals of Lipschitz-Hankel type involving products of Bessel functions[END_REF]. From these studies it can be seen that the product recognition problem falls into two categories: a) systems where users get online information regarding a product by snapping a picture of a product with a camera-device (smartphone or tablet) and b) systems where a vertical search is performed based on domain-specific search such as "Philips TVs" domain. In the latter case, focused crawlers/parsers extract the informative content from web product-pages and try to recognize the product by taking into consideration text content and corresponding product images that appear in the web product page [START_REF] Eason | On certain integrals of Lipschitz-Hankel type involving products of Bessel functions[END_REF]. Almost all the existing PR approaches involve CBIR techniques. A wide variety of CBIR have already been proposed in the literature (e.g. [START_REF] Liu | Content-based image retrieval using color difference histogram[END_REF][START_REF] Chen | Content-Based Image Retrieval with LIRe and SURF on a Smartphone-Based Product Image Database[END_REF][START_REF] Bay | Speeded-up robust features (SURF)[END_REF][START_REF] Van Zijl | Content-Based Image Retrieval with Cellular Automata[END_REF][START_REF] Banda | When too similar is bad: A practical example of the solar dynamics observatory content-based image-retrieval system[END_REF][START_REF] Chatzichristofis | Combining color and spatial color distribution information in a fuzzy rule based compact composite descriptor[END_REF]).

In this work a new CBIR method is proposed for the PR challenging task. Considering the scenario that customer seeks for a product, our goal is to create a search engine that recognizes products from websites with high commercial transactions. Images and description of products constitute the informative content of webpages. The downloaded images from retailers' web pages are represented by a vector of visual features. The proposed CBIR algorithm retrieves similar images from a database which includes images for several products from a large set of manufacturers.

The proposed CBIR method utilizes the Empirical Mode Decomposition (EMD) in order to extract the edges of images. The highest local spatial variations are localized with the representation of images in the Hue-Saturation-Value (HSV) color space. EMD is an adaptive filtering method introduced by Huang and et al. [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] in 1998. It does not require any a priori assumption regarding the basis and it is based on local characteristics of the signal (extrema). Some modifications of the EMD applied for Image Processing can be found in [START_REF] Liang | Edge detection via window empirical mode decomposition[END_REF][START_REF] Guanlei | Improved bi-dimensional EMD and Hilbert spectrum for the analysis of textures[END_REF][START_REF] Bhuiyan | Fast and adaptive bidimensional empirical mode decomposition using order-statistics filter based envelope estimation[END_REF]. After the extraction of the first component a local-threshold adaptive approach [START_REF] Gatos | Adaptive degraded document image binarization[END_REF] is performed in order to extract qualitative edges and interested points of images. Then a quantized color histogram is created. Twenty five [START_REF] Sauvola | Adaptive document image binarization[END_REF] perceptual colors are selected based on the color image quantization method proposed in [START_REF] Wang | Robust image retrieval based on color histogram of local feature regions[END_REF].

The rest of this paper is organized as follows. In Section 2, an overview of the system architecture is introduced. The proposed CBIR is presented in Section 3. In Section 4, performance results are shown. Section 5 concludes the paper.

2

System Architecture

Overview of the System

The World Wide Web (www) has become an important data resource. The rapid growth of its size reveals the need for effective web search engines. However, the task of extracting information from www is non-trivial. Conventional techniques for extracting information include web crawlers.

Crawler is a program that downloads web pages and extracts raw information regarding the downloaded web pages [START_REF] Eason | On certain integrals of Lipschitz-Hankel type involving products of Bessel functions[END_REF]. In general web crawler is a keystone in many web services (e.g. shopping search engines, Internet Archives). In our system a PR search engine is implemented taking into consideration the vertical search approach. In the vertical search scenario focused web crawlers are utilized optimized for specific domain. Focused web crawlers which seek and maintain pages related to a specific topic, provide a narrow segment of the www, and thus they constitute a robust and effective component for the design of a potential high-performance PR system.

Since our goal is the PR challenging task, it is clear that an efficient approach requires the integration of both textual and visual content. A brief description of the overall workflow of our system is shown in Fig. 1. From Fig. 1 it can be seen that the PR system incorporates four major components: a) the Scheduler, b) the Focused Crawler, c) the HTML-Parser, d) the Search UI component and e) the PR Unit. The Scheduler is responsible to create crawl sessions devoted to particular websites. The Scheduler triggers Focused Crawlers to navigate to websites and extract the product pages from online retailers. The HTML-Parser processes the source code from a product page and selects the informative content including product title, price, description, availability, metadata and product images. The Validation Unit verifies the matching decision produced from the CBIR algorithm between a product page and a product that appears in the Product Catalog database through a text analysis process. The Search UI component is the graphical user interface. It includes features such as product querying, product purchasing, product description, and list of prices regarding an item found in different retailers.

The PR unit is the most intelligent part of the system. This component incorporates the proposed CBIR technique and analyzes images in order to represent them in an appropriate feature space where the similarity between images is more effective. In general CBIR process implies several challenging issues, including analysis of lowlevel features (e.g. color, shape, texture, spatial layout, etc.) and creation of feature vectors. The size of the feature vector affects the quality and the speed of a CBIR system. Basically, the key to attain efficient CBIR is to compromise these factors. Thus, a real-time and accurate CBIR system, for the representation of images, requires a small size of feature vector, while at the same time it has to maintain its accuracy through a large set of features.

Proposed Method for PIR

Let C be the color product image with width w and height h. Considering the HSV color space the Value channel V of the color image C is selected as the source image for the edge detection method. The Value component of the HSV color space represents the brightness of an image and it's vital to the definition of edge points. The reason for the selection of the HSV color space for the edge detection process is the fact that HSV corresponds similar to the way humans perceive color [START_REF] Herodotou | A color segmentation scheme for object-based video coding[END_REF] and it is well-known, that in comparison to RGB, it is more accurate in shadow distinction. However, it should be mentioned that the selection of the HSV color is not panacea and different color spaces may also be suggested in this step. For the edge detection process a preprocessing analysis of the grayscale image V is performed by utilizing the EMD to filter the image and extract a more informative component. The application of the EMD for the definition of edge points is described below.

Bidimensional Empirical Mode Decomposition BEMD

EMD is a proper signal processing technique for analyzing nonlinear and nonstationary data. It is an adaptive technique and the extraction of the basis is attained based on local characteristics (extrema) of the analyzed signal. Bearing in mind the Hilbert transform and the concept of instantaneous frequency the cooperation of these definitions with the EMD forms a powerful tool for time-frequency analysis, named Hilbert-Huang-Transform (HHT) [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. The EMD decomposes a signal into its intrinsic mode functions (IMFs). These components exhibit well behavior under the Hilbert transform.

In recent years, the EMD was applied to analyze also two-dimensional data such as images. Proposed versions of the BEMD can be found in the literature [START_REF] Liang | Edge detection via window empirical mode decomposition[END_REF][START_REF] Guanlei | Improved bi-dimensional EMD and Hilbert spectrum for the analysis of textures[END_REF][START_REF] Bhuiyan | Fast and adaptive bidimensional empirical mode decomposition using order-statistics filter based envelope estimation[END_REF]. In the proposed CBIR algorithm a modification of Window Empirical Mode Decomposition (WEMD) [START_REF] Liang | Edge detection via window empirical mode decomposition[END_REF] is used in combination with the improved bidimensional EMD (IEMD) [START_REF] Guanlei | Improved bi-dimensional EMD and Hilbert spectrum for the analysis of textures[END_REF]. The proposed approach of the BEMD is given below:

1. Initialization. Set r 0 =V, i=1 where V is the grayscale image from the Value channel of the color product image. 2. Extract all extrema points (maxima-minima) in r i-1 3. Calculate the mean signal m i-1 described as follows:

(a) Initialization: Set the maximal size N of the window and set its original size M to 3. (b) If the number of maxima is equal with the number of minima in the window M, then compute the mean value from all pixels appeared in the window of size M, and turn to step (d); 2)-( 5) until i=n where n equals to the number of IMFs.

m i-1 ; 4. Compute imf i-1 =r i-1 -m i-1 , r i = m i-1 ,i=i+1; 5. Repeat steps (
The summation of all the IMFs and the residue image R can reconstruct the original grayscale image V.

(1)

The maximum size N of the window can be determined based on the maxima and minima maps created from the source image. The Euclidean distance to the nearest other local maximum (minimum) point is estimated and denoted as D i (D j ) and the maximum size N of the window is computed by the following equation:

, . (2)
where min/max{} denotes the minimum/maximum value of the elements in the ar-ray{}.

The extrema points are defined based on the method proposed in [START_REF] Guanlei | Improved bi-dimensional EMD and Hilbert spectrum for the analysis of textures[END_REF]. According to this method five different types of extrema are proposed. Assuming a square mask window S of odd size, the center pixel S i,j is an extrema if one of the following conditions are fulfilled:

 First type of extrema:

If S i,j =S i-1 , j =S i+1,j and

(3)

S i,j is a local maximum. If S i,j =S i-1,j =S i+1,j and (4)
S i,j is a local minimum.

 Second type of extrema:

If S i,j =S i,j-1 =S i,j+1 and (5)
S i,j is a local maximum. If S i,j =S i,j-1 =S i,j+1 and (6)
S i,j is a local minimum.

𝑆 𝑙,𝑗 > 𝑆 𝑙,𝑗 +𝑞 𝑞=±1 𝑖+1 𝑙=𝑖-1 𝑆 𝑙,𝑗 < 𝑆 𝑙,𝑗 +𝑞 𝑞=±1 𝑖+1 𝑙=𝑖-1 𝑆 𝑖,𝑙 > 𝑆 𝑖+𝑞,𝑙 𝑞=±1 𝑗 +1 𝑙=𝑗 -1 𝑆 𝑖,𝑙 < 𝑆 𝑖+𝑞,𝑙 𝑞=±1 𝑗 +1 𝑙=𝑗 -1
 Third type of extrema:

If S i,j =S i-1,j+1 =S i+1,j-1 and (7) S i,j is a local maximum. If S i,j =S i-1,j+1 =S i+1,j-1 and (8)
Si,j is a local minimum.

 Fourth type of extrema:

If S i,j =S i-1,j+1 =S i+1,j-1 and (9) S i,j is a local maximum. If S i,j =S i-1 , j-1 =S i+1 , j-1 and (10)
S i,j is a local minimum.

 Fifth type of extrema:

If (11) S i,j is a local maximum. If (12)
S i,j is a local minimum.

Edge detection points

In this section the edges of the product image are defined. Once the first IMF component (IMF1) is extracted, through the application of BEMD to the Value channel, a local edge detection technique is used. The applied edge detection algorithm is an adaptive technique introduced in [START_REF] Gatos | Adaptive degraded document image binarization[END_REF]. This method deals with degradations appeared

𝑆 𝑖+𝑑,𝑗 -𝑑 > 𝑆 𝑖+𝑑+𝑞,𝑗 -𝑑 +𝑙 +1 𝑙=-1 +1 𝑞=-1 𝑙≠-𝑞 +1 𝑑 =-1 𝑆 𝑖+𝑑,𝑗 -𝑑 < 𝑆 𝑖+𝑑+𝑞,𝑗 -𝑑 +𝑙 +1 𝑙=-1 +1 𝑞=-1 𝑙≠-𝑞 +1 𝑑 =-1 𝑆 𝑖+𝑑 ,𝑗 +𝑑 > 𝑆 𝑖+𝑑 +𝑞,𝑗 +𝑑+𝑙 +1 𝑙=-1 +1 𝑞=-1 𝑙≠𝑞 +1 𝑑 =-1 𝑆 𝑖+𝑑 ,𝑗 +𝑑 < 𝑆 𝑖+𝑑 +𝑞,𝑗 +𝑑+𝑙 +1 𝑙=-1 +1 𝑞=-1 𝑙≠𝑞 +1 𝑑 =-1 𝑆 𝑙,𝑗 > 𝑆 𝑖+𝑙,𝑗 +𝑞 +1 𝑞=-1 +1 𝑙=-1 𝑆 𝑙,𝑗 < 𝑆 𝑖+𝑙,𝑗 +𝑞 +1 𝑞=-1 +1 𝑙=-1
due to shadows, non-uniform illumination, low contrast and high signal-dependent noise. Moreover, it runs automatically without requiring the adjustment of a set of design parameters in each execution. It is worthy to mention that this method belongs to local-threshold techniques. In comparison to global-threshold methods, these techniques provide better performance since they utilize the local area information in order to estimate a threshold value per pixel. Each pixel is classified based on its threshold value into foreground or background class. A flow chart with the major steps of the edge detection algorithm is presented in Fig. 2.

BEMD Sauvola's thresholding Interpolation

Classification of pixels based on contrast criterion BS(x,y)-imf(x,y)>d(BS(x,y))

Post-processing Noise removal

(imf) (B) (BS) (FT) Final binary image

Fig. 2. Flow chart of the edge detection process.

In the first step of the original method a preprocessing filtering stage of the grayscale is required by utilizing a low-pass Wiener filter [START_REF] Jain | Fundamentals of Digital Image Processing[END_REF]. This filter is usually used for image restoration. However, in this paper this step is skipped since a preprocessing filtering is already achieved through the extraction of the first IMF component. The second step of the algorithm implements a rough classification of pixels into foreground and background clusters. This is attained via the Sauvola's approach [START_REF] Sauvola | Adaptive document image binarization[END_REF]. The initial local threshold estimation is given by: [START_REF] Bay | Speeded-up robust features (SURF)[END_REF] where m is the mean value of the rectangular window, T is the threshold value for the center pixel of the window and s is the variance value. R is a constant fixed to 128 and k is set to 0.2 according to [START_REF] Gatos | Adaptive degraded document image binarization[END_REF]. Assuming the first IMF image imf the binary image B is extracted where 1s' represent the foreground pixels.

A background surface BS is then computed by following two rules: For zero pixels in image B the corresponding value at BS image equals to imf corresponding values. For the foreground pixels the corresponding values of the BS image are computed via interpolation described in the following equations: [START_REF] Van Zijl | Content-Based Image Retrieval with Cellular Automata[END_REF] otherwise [START_REF] Banda | When too similar is bad: A practical example of the solar dynamics observatory content-based image-retrieval system[END_REF] The window size of dx x dy is set to 5x5.

The next step of the process involves the final thresholding phase, where the binary image FT is estimated taking into account the imf image and the background surface BS. The foreground pixels are determined based on the distance value between the imf image and the background surface BS. According to this method, a threshold d is

𝑇 = 𝑚 + (1 -𝑘 1 - 𝑠 𝑅 ) 𝐵𝑆 𝑖, 𝑗 = 𝑖𝑚𝑓 𝑖, 𝑗 𝑖𝑓 𝐵 𝑖, 𝑗 = 0 𝐵𝑆 𝑖, 𝑗 = (𝑖𝑚𝑓 𝑖𝑥, 𝑖𝑦 1 -𝐵 𝑖𝑥, 𝑖𝑦 ) 𝑗 +𝑑𝑦 𝑖𝑦 =𝑗 -𝑑𝑦 𝑖+𝑑𝑥 𝑖𝑥 =𝑖-𝑑𝑥 1 -𝐵 𝑖𝑥, 𝑖𝑦 𝑗 +𝑑𝑦 𝑖𝑦 =𝑗 -𝑑𝑦 𝑖+𝑑𝑥 𝑖𝑥 =𝑖-𝑑𝑥
proposed and a foreground pixel is recognized if the calculated distance value exceeds the threshold d value. The final thresholding binary image at each point (x,y) is given by the following equations: [START_REF] Chatzichristofis | Combining color and spatial color distribution information in a fuzzy rule based compact composite descriptor[END_REF] The threshold d for each pixel simulates a logistic sigmoid function and it is given by the following equation: [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] where the variable b represents the average background value, the parameter q is equal to 0.6, the p 1 is set to 0.5 and p 2 to 0.8 according to [START_REF] Gatos | Adaptive degraded document image binarization[END_REF] after extensive experimental work. The variable δ is the average distance between foreground and background values of the pixels and it is expressed by the following formula: [START_REF] Liang | Edge detection via window empirical mode decomposition[END_REF] The final step of the edge detection method includes the post-processing enhancement of the binary image by applying shrink and swells filters. These filters aim to remove isolated pixels and fill potential breaks and gaps. A detailed description of the algorithm can be found in [START_REF] Schilling | Fundamentals of Robotics Analysis and Control[END_REF].

Interested Points Detection

In the proposed method neurobiological findings are taken into consideration for the design of the PIR method. It is well-known that the ganglion cells (GC) constitute the unique output of the retina to the primary visual cortex. In this early neuronal vision component (retina) a preprocessing step is responsible to filter the initial visual information. This phase is required to minimize the amount of the transmitted information due to the limited bandwidth of the channel between retina and primary visual cortex. In retina the GCs and the connected receptive fields (RF) are the major components. The RFs of GCs have a center-surround mechanism. The RF consists of two concentric areas, the center and the surround that have antagonistic effects on the response of a GC. Consequently, this mechanism sends information regarding two different regions and the brain is responsible to run the "filling-in" mechanism in order to create the scene. Inspired by the biological visual system of mammals, only the interested points of images are taken into consideration for the next processing step. Thus, motivated by the biological vision systems, the considered interested points are the edges of images and the surrounded pixels of them. Below follows a detailed description of the selection of the interested points: Let L={P 1 (x 1 ,y 1 ) , P 2 (x 2 ,y 2 ), …, P n (x n ,y n )} be a set of n edge points extracted from the edge detection algorithm. Suppose that x i , y i are the x-coordinate and y-coordinate respectively of the ith point. Assume a window of constant size k x k. We slide the window to each edge point and seek for the interested points in the (k 2 -1) neighbor points. For each edge point a unique window k x k is associated containing k 2

𝐵(𝑥, 𝑦)

𝑦 𝑥 tial interested points. It is reasonable to assume that the edge points belong to the set of the interested points. Thus, our investigation is located only to the rest (k 2 -n) neighbor points. n equals to the number of edge points in the window. In our analysis the constant size k of the window is set to 3. The entropy is utilized as criterion for the selection of the rest interested points. Let f 1 , f 2 , …, f m be the observed gray-level frequencies in the sliding window of the imf image. The percentage of occurrence of a specific gray level p i equals to: [START_REF] Guanlei | Improved bi-dimensional EMD and Hilbert spectrum for the analysis of textures[END_REF] The entropy H of a set of pixels appeared in the window is defined as: [START_REF] Bhuiyan | Fast and adaptive bidimensional empirical mode decomposition using order-statistics filter based envelope estimation[END_REF] where s represents the number of pixels in the set. For each gray level appeared at the window the algorithm removes the corresponding pixels and estimates the entropy of the subset that includes the remaining pixels. Our target is to find the gray level that produces the maximum entropy of the subset of pixels in the window and therefore to select only the subset of pixels that it is more heterogeneous.

Normalized Quantized Histogram

Color is one of the essential image indexing features in CBIR. The description of color information in images can be expressed by the means of color histograms [START_REF] Wang | Robust image retrieval based on color histogram of local feature regions[END_REF], color moments [START_REF] Stricker | Color indexing with weak spatial constraints[END_REF] and color correlograms [START_REF] Huang | Image indexing using color correlograms[END_REF]. In our approach a normalized quantized histogram is proposed to create the feature vector of an image. Bearing in mind the performance of our model, the computational time of the retrieval process is one of the keystone design parameters. Thus, it is appropriate to reduce the color space in order to provide fast indexing even in the cases of use of big databases. This leads to incorporation of a small set of colors. In the present study 25 perceptual colors are selected from the RGB color space. These 25 colors divide uniformly the RGB space and are proposed in [START_REF] Wang | Robust image retrieval based on color histogram of local feature regions[END_REF]. Table 1. presents the quantized RGB color palette. It is worthy to mention that additional features may be included in the proposed approach. However, our goal is to investigate the performance of the model by taking into account only the color of the interested points avoiding the potential noise from the rest pixels (background) and maintaining the small size of the feature vector. After the recognition of the interested points of the image, the quantized image Q is created and the normalized quantized histogram is calculated. Given a color space containing L (25) color bins, the color histogram of image Q is expressed by: [START_REF] Gatos | Adaptive degraded document image binarization[END_REF] where n x is the total number of interested points in the xth color bin and N the total number of interested points of the examined image.

Experimental results

In this paper a robust PIR is proposed based on color histogram of the interested points of product images. A brief description of the process can be summarized as follows:

1. The RGB query image is transformed to HSV color space. The Value channel of the query image is selected for the preprocessing filtering step. 2. The BEMD is applied to the Value channel and the first IMF component is extracted. 3. An adaptive local threshold method is utilized for the recognition of edge points. 4. The entropy is selected as the criterion for the definition of interested points. 5. The query image is quantized by utilizing the color lookup table. 6. The normalized color histogram of the interested points is created. 7. The similarity between color images is computed by using the color histogram of interested points. The Tanimoto coefficient [START_REF] Duda | Pattern Classification[END_REF] is used as a metric of the distance.

The proposed low level feature vector has been incorporated in the product recognition platform named "Agora" which is the outcome from the cooperation of R&D and Software department of ChannelSight. Our tests derived from a subset of our database [START_REF]ChannelSight company information[END_REF] containing currently 48.071 Philips product images for the region of Great Britain. The database contains images from different product categories (e.g. ironing, cooking, headphones, etc.) and regions. The PR system can be applied to several types of products from different brands. In our PR system a text analysis tool is also integrated to validate and recognize products from retrieved images provided by the proposed PIR algorithm. However, the scope of the proposed paper is to examine the effectiveness of the PIR technique. Only the first 10 images from the retrieval process are taken into consideration for the recognition of query products and the product is considered as recognized only if an image related to the query product is included in the first ten retrieved images.

𝐻 𝑥 = 𝑛 𝑥 𝑁 (𝑥 = 1,2, . . . , 𝐿)

Eight categories from Philips's product catalog are selected and a set of query images for each category is created from products appeared in Amazon retailer. Table 2 presents the selected categories, the number of existing images in Philips's product catalog and the percentage of retrieved products. From the experimental results it can be seen that the proposed method is resilient and can retrieve products that have some variations in color like the first pair in Fig. 3 and the last pair of images in Fig. 5. The color histogram, which is utilized as a feature for indexing images, is a statistical feature with low computational cost. It is also invariant to rotation and scale of image scene and to any displacement of objects in the image. Experimental results reveal the effectiveness of the proposed method in those cases. Moreover, it is worthy to mention the strength of the proposed technique in cases like the first and the penultimate example in Fig. 5. From this figure it can be seen that the stamp "Promix" in the first example and the shadow of product in the third pair of images does not drastically affect the color histogram, since the only pixels that are taken into account for the creation of the normalized and quantized histogram are the interested points of an object and not the entire object. Fig. 6 presents the edges, the quantized images and the extracted interested points of the penultimate pair of images in Fig. 5. Fig. 7 presents examples in different cases such as rotation, viewpoint changes and illumination variations. Fig. 8 shows the normalizedquantized color histograms. The first histogram refers to the Amazon's image while the second histogram refers to the image appeared in the Philips product catalog.

From the second histogram it can be seen that the quantity of black color is increased, however, the black color is still limited to low values due to the introduction of notion of interested points. The orange colored region reveals the quantity of white in the tested images. The total retrieval time is not constant and depends on the number of products included in the database for the specified category. It should be mentioned that the images that appeared in the product catalog are already represented as vectors and the only time load is due to the estimation of Tanimoto distances between query image and product images included in specified product category.

To validate the effectiveness of the proposed method, a comparison study between the proposed method and a number of existing CBIR in the literature is given in Fig. 9. The Scalable Color Descriptor (SCD), the Edge Histogram (EH) and the Fuzzy Color and Texture Histogram (FCTH) [START_REF] Chatzichristofis | Combining color and spatial color distribution information in a fuzzy rule based compact composite descriptor[END_REF] are used. SCD and EH are part of MPEG-7 specification for multimedia content description [START_REF] Manjunath | Color and texture descriptors[END_REF]. In our tests, in the case of SCD 256 coefficients are used and the number of discarded bitplanes was set to 0. The size of feature vector for SCD and EH is 256 and 80 respectively. The FCTH method provides a set of 192 features. Fig. 9 reveals the performance of each CBIR in the examined Philips product categories. From this figure it can be seen that the product recognition performance of the CBIR methods is similar except the case of Headphones where the database is larger. In this case the proposed method and the SCD provide better results. However, it should be mentioned that the performance of the proposed method is attained with the smallest size of features.

Conclusion

In this paper a new PIR technique is proposed for the recognition of products in the ecommercial web sites. The proposed method filters the Value channel of HSV color space and extracts the first IMF through the BEMD. The first IMF provides the highest frequency content, thus it is the most informative component. An adaptive algorithm with local thresholds is applied to the first IMF for the detection of qualitative edges. The interested points are defined following the condition of maximum entropy and the color histogram is created based on a subset of perceptual twenty five colors. The experimental results reveal the effectiveness of the proposed method. The PIR returns the associated product providing an average correct recognition over 90%. The execution time of query images with different sizes reveal that the proposed method is a real-time technique while the short feature vector is in agreement with requirements for design a fast image retrieval system. Further research involves the comparison of the proposed technique with conventional image retrieval methods in the application of product image retrieval. Moreover, future work will focus on the challenge of creating an effective and robust distributed system for product recognition, integrating also adaptive techniques that combine text analysis algorithms and CBIR methods.

Fig. 1 .

 1 Fig. 1. Workflow of the PR system.

  𝑚𝑖𝑛 𝐷 𝑖 min 𝐷 𝑗 (c) Set M=M+2. If M<N, turn to step (b), otherwise compute the mean value from all pixels appeared in the window of size N; (d) Turn to the next pixel, set M=3, turn to step (b), repeat until all pixels of the image are processed, then compute the mean surface

Fig. 3 .

 3 Fig. 3. Results from Ironing and Face-Shavers category. The left image of each pair is from Philips's catalog.

Fig. 4 .

 4 Fig. 4. Results from Cooking and Headphones category.

Fig. 5 .

 5 Fig. 5. Results from Food-Preparation category.

Fig. 6 .

 6 Fig. 6. The quantized images, the edges and the interested points.

Fig. 7 .

 7 Fig. 7. Results from Audio, Coffee and Hair removal category.

Fig. 8 .

 8 Fig. 8. Normalized and quantized color histograms for the query and retrieved image in the penultimate pair of images in Fig. 7.

Fig. 9 .

 9 Fig. 9. Product Recognition results of the examined methods in five product categories of Philips.

Table 1 .

 1 Color-lookup table.

	Color	R	G	B	Color	R	G	B
	Black	0	0	0	Plum	146 109	0
	Sea green	0	182	0	Teal	146 182 170
	Light green	0	255 170	Dark red	182	0	0
	Olive green	36	73	0	Magenta	182	73	170
	Aqua	36	146 170 Yellow green 182 182	0
	Bright green	36	255	0	Flouro green 182 255 170
	Blue	73	36	170	Red	219	73	0

Table 2 .

 2 Details of the results.

	Categories	Amazon	Retrieved	Percentage of	Images in
		Products	Products	correct retrieval	Product Catalog
	Ironing	37	34	91.89 %	327
	Shavers	23	21	91.30 %	266
	Cooking	11	11	100 %	63
	Food Preparation	26	25	96.15 %	149
	Coffee	21	21	100 %	480
	Headphones	91	89	97.80 %	533
	Audio	27	26	96.29 %	351
	Hair Removal	23	22	95.65 %	137

Table 3 .

 3 Execution time for image representation.

	Size of image	64 x 64	128 x 128	256 x256	512 x 512
	Execution time (s)	0.031214	0.056614	0.186249	0.678315


Samaneh Mazaheri
email: mazaheri.samaneh@yahoo.com

Suhaiza Puteri

Rahmita Sulaiman

Mohd Zamrin Wirza

Fatimah Dimon

Rohollah Moosavi Khalid

Tayebi

Puteri Suhaiza Sulaiman

Rahmita Wirza

Mohd Zamrin Dimon

Fatimah Khalid

Rohollah Moosavi Tayebi

Uncertainty Estimation for Improving Accuracy of Non-rigid Registration in Cardiac Images

Keywords: Uncertainty, Multimodality Image Registration, Hybrid, Featurebased, Intensity-based, Echocardiography, Computed Tomography (CT)

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Echocardiography imaging is an attractive alternative to MRI and CT during surgery due to ease of use, its safety, minimal disruption of the procedure, comparatively low cost, and lack of compatibility problems between echocardiography imaging and standard operating theater equipment. However, due to inferior overall image quality of echocardiography images than MRI or CT images, clinicians often have problem in assessing target organs in the clinical applications. To overcome this limit of echocardiography imaging, there have been several efforts to display or register an echocardiography image with its corresponding high-quality CT (or MRI) image to provide clear information related to anatomy and particular target [1, 2, 3, and 4]. Some structure target are often nearly invisible in an echocardiography image, the corresponding CT image is useful for providing information to clinicians for diagnosis and planning of intervention. To obtain a CT image correctly aligned to an echocardiography image for the heart, non-rigid registration between CT and echocardiography images is essential [START_REF] Lee | Non-rigid registration between 3D ultrasound and CT images of the liver based on intensity and gradient information[END_REF].

Registering pre-operative CT scan to echocardiography image of the patient heart can assist the guidance of surgical process and instruments, during cardiac surgery [START_REF] Huang | Dynamic 2D Ultrasound and 3D CT Image Registration of the Beating Heart[END_REF]. Also it can provide image-guidance for procedures that involves both anatomical and functional imaging.

The present paper briefly reviews the two-level non-rigid registration algorithm between CT images and echocardiography of the heart, and then by using uncertainty function tries to improve the registration accuracy through considering local deformations, which cannot be compensated using a global registration algorithm. In the registration technique, segmentation result information and intensity value information of the ventricles, are used for registration procedure.

Uncertainty is inherently an important issue to all image analysis tasks. However, few common algorithms estimate a measure of uncertainty for the calculated results. Automatic error detection would be very valuable in the practical use of medical image analysis tools, as it could provide the clinicians with more confidence on whether or not to use a computer-generated analysis. Registration uncertainty is important information to convey to a surgeon when surgical decisions are taken based on registered image data. However, conventional non-rigid registration methods only provide the most likely deformation [START_REF] Risholm | Summarizing and Visualizing Uncertainty in Non-rigid Registration[END_REF].

There is a level of uncertainty associated with registering images that have any differences. A confident registration with a measure of uncertainty is critical for many change detection applications such as medical diagnostics. Several notable papers have attempted to quantify uncertainty in image registration in order to compare results [START_REF] Simonson | A Statistics Based Approach to Binary Image Registration with Uncertainty Analysis[END_REF][START_REF] Domokos | Parametric estimation of affine deformations of binary images[END_REF].

In this paper, we describe a method for non-rigid registration of cardiac images and quantify the uncertainty of the estimated transformation. The registration method was validated on a clinical dataset with results comparable to previously published methods, but with the added benefit of also providing uncertainty estimates which may be important to take into account during cardiac surgery procedures. Improvements are demonstrated in accuracy and robustness over non-rigid registration between CT and echocardiography images and it is shown how to summarize this uncertainty in terms of uncertainty quantification and user interaction. We also emphasis that non-rigid registration results come with a level of uncertainty and conveying the uncertainty, and not only the most likely estimation, is clinically important.

Previous works in our lab towards the goal of Computer Assisted Medical Research has included the develop an automatic coronary arterial tree extraction in angiograms [START_REF] Moosavi Tayebi | A fast and accurate method for automatic coronary arterial tree extraction in angiograms[END_REF], cardiac ultrasound fusion system development [START_REF] Mazaheri | Hybrid Pixel-based Method for Cardiac Ultrasound Fusion Based on Integration of PCA and DWT[END_REF], Wavelet enhancement for xray angiogram [START_REF] Moosavi Tayebi | Using Wavelet for X-ray Angiography Enhancement[END_REF], CT angiography components categorization and coronary artery enhancement [START_REF] Moosavi Tayebi | Cardiac Components Categorization and Coronary Artery Enhancement in CT Angiography[END_REF], a hybrid non-rigid registration method based on pre-and intra-operative cardiac images [START_REF] Simonson | A Statistics Based Approach to Binary Image Registration with Uncertainty Analysis[END_REF], and some surveys including review on segmentation approaches [START_REF] Mazaheri S | Echocardiography Image Segmentation: A Survey[END_REF][START_REF] Mazaheri S | A Review of Ultrasound and Computed Tomography Registration Approaches[END_REF], and review on registration of cardiac images [START_REF] Moosavi Tayebi | Coronary artery segmentation in angiograms with pattern recognition techniques -A survey[END_REF].

In this paper, we are trying to quantify the uncertainty measurement in the hybrid registration method which is presented in [START_REF] Mazaheri S | Hybrid Non-rigid Registration Framework for CT and Echocardiography Images[END_REF] in details. This paper is organized as following; Section 2, refers to describe the hybrid registration procedure briefly; in Section 3, uncertainty definition and evaluation is explained; Section 4, registration uncertainty quantification is explored and finally in Section 5, some future directions are presented and the paper is concluded.

Hybrid Non-Rigid Registration

In this work, interested region in two different modality images are segmented by a level set model which is proposed in [START_REF] Mazaheri S | A Robust and Accurate Segmentation of Echocardiographic Images through a Geometrically Constrained Level Set[END_REF], in first modality -echocardiography-, and in the second modality -computed tomography-by k-means clustering technique [START_REF] Ma | Segmentation of Small Animal Computed Tomography Images Using Original CT Values and Converted Grayscale Values[END_REF], and then images are registered together with the two-level registration method presented in [START_REF] Mazaheri S | Hybrid Non-rigid Registration Framework for CT and Echocardiography Images[END_REF].

The extracted features or segmented contours are represented by two different function sets of images. An objective function which is described based on segmentation result defines the mapping between the segmented contours. The presented registration method uses two relatively separated discontinuity sets to represent the contour sets of the associated images to tackle the registration problem between these contour sets.

For the ambiguity problem of the correspondence, we apply the idea of consistent registration [START_REF] Christensen | Consistent Image Registration[END_REF][START_REF] Johnson | Consistent Landmark and Intensity-Based Image Registration[END_REF] to simultaneously estimate the forward and reverse transformations and to constrain one transformation to be the inverse of the other one. In this case, the contour sets of images have equal influence on the registration. Thus, the registration technique is one-to-one in the way that it allows to determine one-to-one correspondences between the contour sets of the images.

Symmetric one-to-one contour matching is not only more sound in mathematical sense, but also very important in determining the one-to-one correspondence of the same anatomical structures in two different images. The ability to correlate echocardiography images with pre-operative dynamic CT images would provide the surgeon with complementary information captured in two different modalities, and greatly facilitate the interpretation of the images.

Registration and segmentation of the ventricular myocardium in the images are important steps to extract useful qualitative or quantitative information about myocardial intervention. The hybrid non-rigid registration algorithm uses intensity value information and segmentation result information through a level set based segmentation to register echocardiography and CT images.

Different similarity measures are considered based on the correlation of the intensity value information and segmentation result information between echocardiography and CT images, according to the image characteristics of the features, including ventricle endocardial and ventricle surfaces in two images.

To represent the similarity measure between echocardiography and CT images, two objective functions based on the segmentation result information and intensity value information are adopted. The registration technique conducts the segmentation-based and intensity-based registration procedures according to corresponding similarity measure. The images have been preprocessed in order to increase the relationship between CT and echocardiography pixels. Image pre-processing is necessary to make the modalities as close as possible (interface enhancing in CT scans, filtering, and artefact and acoustic shadow removal and filtering in echocardiography images). The images have been preprocessed in order to increase the relationship between CT and echocardiography pixels.

The discussed non-rigid registration is less restrictive regarding the structure which is due to the faster implementation and inherently smooth nature of the basic functions. No search space or point to point mapping is needed either. It is possible to achieve smooth results and handle multiple point set models at the same time.

Uncertainty Definition and Evaluation

An important, but somehow neglected aspect of non-rigid registration is: how is it possible to quantify and visualize the registration uncertainty. The significance can be simply realized in the range of surgeries. Critical functional regions of the target are usually determined in the pre-operative domain. Current image guided navigation systems employ rigid registration to build the mapping between the intra-and pre-operative domains, but there is a move towards non-rigid registration recently, as we described in details in [START_REF] Mazaheri S | Hybrid Non-rigid Registration Framework for CT and Echocardiography Images[END_REF].

Regular registration techniques gives the most likely deformation and hence the most likely position of a functional region. However, non-rigid registration can behave randomly in the way that small changes of data parameters can produce completely different result.

In addition, metrics like resection of tissue and degraded intra-operative image quality can all contribute to increasing in the registration uncertainty. Thus, for a surgeon, the uncertainty in the estimated position of functional regions can be exactly as significant information as the most likely estimate of the position.

Registration uncertainty can be valuable information to portray to a surgeon/interventionist, especially when important clinical decisions are based on the registered data. In surgeries, target points are defined in the pre-procedural image and transformed into the intra-procedural image space to guide the surgeon during the procedures. It is important to understand that there is an uncertainty associated with the estimated transformation and consequently the transformed points.

In this paper, we define the registration uncertainty and explore its use for diagnostic purposes. We evaluated the registration accuracy via uncertainty estimation on images from 16 subjects in two cardiac phases, systolic and diastolic, acquired from UiTM hospital in Malaysia under supervision of cardiac surgeon.

Since, results indicate that the described registration system is capable of producing segmentation results with high robustness and high accuracy with minimal user interaction across all subject groups, high registration uncertainty will indicate functional abnormality of the subject. We show that registration uncertainties are good predictors for the functional abnormality of the subject.

Several challenges are common in inter-subject registration, which can cause the occurrence of local uncertainties of the transformations [START_REF] Heinrich | Uncertainty Estimates for Improved Accuracy of Registration-Based Segmentation Propagation using Discrete Optimization[END_REF]:

 Missing one-to-one correspondences of anatomies across subjects  Mismatch of structures due to local minima of the registration  Untrustworthy image information due to acquisition noise and artefacts No matter how robust a registration technique is, it is important to have the ability to alert the user if the uncertainty of the registration quality is high. High uncertainty can either be a sign of an unreliable registration result or of an abnormal cardiac anatomy. The integration of feature-and intensity-based registration has shown the ability to achieve both good robustness and accuracy [START_REF] Mazaheri S | Hybrid Non-rigid Registration Framework for CT and Echocardiography Images[END_REF]. Back to the important, but not fully explored topic in the domain of non-rigid registration; the question which mentioned at the beginning of the section; how is it possible to quantify and visualize the registration uncertainty? Since the hybrid registration combines information from two sources which they are intensity value and segmentation result information, this question can be further divided into uncertainty arising from the segmentation result [START_REF] Risholm | Summarizing and Visualizing Uncertainty in Non-rigid Registration[END_REF][START_REF] Kybic | Bootstrap resampling for image registration uncertainty estimation without ground truth[END_REF] and uncertainty arising from abnormal intensity, which is related to functional abnormality of the subject. The registration uncertainty which comes from intensity relates to abnormal intensity like ischemic cardiomyopathy, while registration uncertainty which comes from segmented region (geometry distribution) corresponds well to abnormal geometry like dilated cardiomyopathy.

Reinitialize the Atlas Position

Fig. 2. Workflow of the uncertainty estimation framework

We also explore the potential of registration uncertainty in improving the robustness: If and only if the uncertainty of the registration is high, the system will ask the user to input additional landmarks to help better initialize the estimated transformation for registration. The landmarks include apex, center of mitral valve, center of left ventricle, center of right ventricle and two right ventricle insertion points (anterior and inferior insertion points of right ventricle) and center of basal plane. [START_REF] Huang | Dynamic 3D ultrasound and MR image registration of the beating heart[END_REF] Registration uncertainty

Uncertainty Definition

Given two images, and , we can estimate a transformation which maps image to , so that a pixel of ( ) correspond to ( ( )) and their intensity values should be similar. Using a probabilistic formulation for the image registration problem [START_REF] Risholm | Summarizing and Visualizing Uncertainty in Non-rigid Registration[END_REF], the uncertainty of a transformation at point can be modeled by the following equation:

( ( )|( ( ), )) = 1 - ( ( ), ) ( ) ( ( )) (( ( ), )) (1)
We model the likelihood term (( ( ), )| ( )) as a normal distribution of the intensity difference between transformed and estimated. Comparably, the prior of the transformation, ( ( )), is modeled as a Rician distribution of the Jacobian determinant of the transformation [START_REF] Rohlfing | Volume-preserving nonrigid registration of MR breast images using free-form deformation with an incompressibility constraint[END_REF]. The distribution is estimated based on the inversion technique proposed in [START_REF] Koay | Analytically exact correction scheme for signal extraction from noisy magnitude MR signals[END_REF]. The Rician distribution is a non-negative and asymmetric distribution which approximates the distribution of the Jacobian determinant well for a given transformation. Finally, ( ( ),

) can be modelled as a constant term.

Uncertainty Quantification and User Interaction

For each pixel in all images, its registration uncertainty can be evaluated and visualized using equation [START_REF] Wein | Automatic non-linear mapping of pre-procedure CT volumes to 3D ultrasound[END_REF]. We can further define the registration uncertainty of a given region by averaging over the region. The quantification of uncertainty can be used to inform the user about how reliable the registration results are. Based on results from the uncertainty analysis, we can design a system that detects functional abnormality with high uncertainty. High uncertainty can either be a sign of an unreliable registration result (registration failure) as well as abnormal cardiac anatomy.

Fig. 3. LV and RV differences; difference between transformed LV and RV, and the intra-procedural CT is shown in green color, while difference between transformed LV and RV, and the pre-procedural echocardiography is shown in pink.

In the cases, we detect registration failure, the user is asked to define 6 landmarks (apex, center of left ventricle, anterior and inferior insertion points of right ventricle, center of right ventricle and center of basal plane). These landmarks are also defined in the atlas.

By introducing knowledge about these additional 6 landmarks, the atlas-to-image registration can be initialized more accurately and all registrations performed correctly. When the registration is considered successful, the uncertainty relates to abnormality of the patient's cardiac anatomy. This is possibly due to the fact that our segmentation algorithm using GCL is designed to segment pathological images well [START_REF] Mazaheri S | A Robust and Accurate Segmentation of Echocardiographic Images through a Geometrically Constrained Level Set[END_REF] and kmeans clustering [START_REF] Ma | Segmentation of Small Animal Computed Tomography Images Using Original CT Values and Converted Grayscale Values[END_REF]. To examine if the uncertainty correlates to the abnormality of the patients, as mentioned in previous section, registration uncertainty which comes from intensity relates to abnormal intensity like ischemic cardiomyopathy; while registration uncertainty which comes from segmented region or geometry distribution corresponds to abnormal geometry like dilated cardiomyopathy. Figure 4 shows that registration uncertainty is a very good predictor for separating ischemic cardiomyopathy, dilated cardiomyopathy, and normal subjects.

The accuracy is significantly improved compared to previous registration methods by utilizing both pre-and intra-image information. We define a system that detects either registration failures or abnormality of the patient's cardiac anatomy using registration uncertainties.

Conclusion

Registration uncertainty is important information to convey to a user when automatic image analysis is performed. Uncertainty information may be used to provide additional diagnostic information to traditional analysis of cardiac function. Cardiac pathology is not always easily detectable in images, but likely to be detected by registration uncertainty. Since, we can define uncertainty for every part of the cardiac anatomy, it is desirable to investigate if the relationship between uncertainty and abnormality could help to detect these pathologies automatically.

The registration uncertainty has the potential to convey information on our confidence in the resulting transformation. To form an uncertainty map describing the registration result, we calculate uncertainty of the transformation at every pixel.

While the particular registration uncertainty proposed in this paper as a preliminary outcome, it is believed, firstly, that in the future, with increasing computational power, more highly developed algorithms and better modeling, it will be practical to get dense information related to uncertainty of the estimated transformation, and secondly that it will be important to share that information with the surgeon in a purposeful and valid way.

Fig. 1 .

 1 Fig. 1. Workflow for hybrid echocardiography-CT registration

Fig. 4 .

 4 Fig. 4. This figure shows registration uncertainty

Acknowledgment

This research was supported by Ministry of Science, Technology and Innovation, Malaysia, under science Fund no. UPM0007353.


Jerzy Martyna

PSO-Least Squares SVM for Clustering in Cognitive Radio Sensor Networks

Keywords: least squares support vector machines, particle swarm optimisation, cognitive radio technology, wireless sensor networks, clustering problem

In this paper, a solution for a cluster formation in cognitive radio networks is presented. The solution features a network-wide energy consumption model for these networks. The particle swarm optimisation (PSO) and least squares support vector machines (LS-SVMs) have been transformed into our clustering problem. The obtained results show that the given hybrid AI system provides a good estimate of a cluster formation. Through extensive simulations, we observed that the PSO-LS-SVM method can be effectively used under various spectrum characteristics. Moreover, the formed clusters are reliable and stable in a dynamic frequency environment.

Introduction

Cognitive radio (CR) technology belongs to a key technology that enables the improvement of the spectrum use in a dynamic manner. The Federal Communications Commission (FCC) reported that a large portion of the assigned spectrum is used sporadically, leading to a significant amount of the unused spectrum [START_REF]Frequency Allocation[END_REF].

CR devices allow one to sense the environment (cognitive capability), analyse and study the detected information, and adapt to the environment (reconfigurable capabilities). Moreover, CR technology can be integrated into wireless sensor networks [START_REF] Akyildiz | Wireless Sensor Networks: A Survey[END_REF]. In this integration, nodes can communicate by tuning their radios to vacant radio channels through a collision-free band instead of through the heavily crowded unlicensed bands. The cognitive radio sensor networks [START_REF] Akan | Cognitive Radio Sensor Networks[END_REF][START_REF] Yau | Cognitive Radio-based Wireless Sensor Networks: Conceptual Design and Open Issues[END_REF] (CRSN) allows us to implement multimedia applications with high bandwidth requirements. Moreover, the CRSN gathers all of the sensed data of observed phenomena.

A two-tiered CRSN is composed of a two-level hierarchy. The first level consists of ordinary sensors according to the application-specific objectives. All ordinary sensors are grouped into clusters. All cluster head sensor nodes that are responsible for data aggregation and sending it to the base station (BS) belong to the second level. It is obvious that each clustering technique used in CRSN must reduce the number of exchanged messages between the ordinary sensors and the cluster head nodes, and minimise the distance between the cluster head nodes to reduce transmission power.

Many solutions for the realisation of clustering in the CRSN have been proposed in the recent literature. Among others, in the paper by Younis et al. [START_REF] Younis | Heed: a Hybrid, Energy-efficient, Distributed Clustering Approach for Ad Hoc Sensor Networks[END_REF], an election clustering method based on nodes residual energy as a primary clustering parameter and node proximity toits neighbours as secondary parameter was proposed. Unfortunately, it was designed for fixed channel settings and can not be used in the CRSN. A deterministic clustering scheme was proposed in the paper by Bradonji and Lazos [START_REF] Bradonji | Graph-based Criteria for Spectrum-aware Clustering in Cognitive Radio Networks[END_REF]. In this approach, clusters are realised by the maximising the summation of idle channels with the number of nodes within the cluster. According to the authors' statement, this solution is stable, but is a NP-hard and violates design constraints for energy and computation. A specific aware clustering method has been presented in the paper by Zhang et al. [START_REF] Zhang | Distributed Spectrum-aware Clustering in Cognitive Radio Sensor Networks[END_REF]. The proposed method demands introducing the k-mean clustering algorithm. In this method each node is treated as cluster head node and then merges cluster heads within each iteration until the number of cluster heads reaches an optimal number that was obtained a priori by the theoretical analysis. Applying artificial intelligence (AI) can be an alternative approach to solving the clustering problem in the CRSN. Cognitive radio sensor networks, that deal with the intelligent assignment and use of the radio spectrum, are the natural study research area for the use of AI techniques. Support vector machines (SVMs) have been developed as an alternative that avoids artificial neural network (ANN) limitations [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF], [START_REF] Vapnik | Statistical Learning Theory[END_REF]. SVM compute globally optimal solutions, unlike tend to fall into local minima [START_REF] Cherkassky | Practical Solution of SVM Parameters and Noise Estimation for SVM Regression[END_REF]. Least squares support vector machines (LS-SVMs), introduced by Suykens et al. [START_REF] Suykens | Least Squares Support Vector Machine Classifier[END_REF], [START_REF] Suykens | Sparse Approximation Using Least Squares Support Vector Machines[END_REF], simplify the training process of standard SVM in a great extent by replacing the inequality constraints with equality ones. A particle swarm optimisation (PSO) [START_REF] Kennedy | Particle Swarm Optimization[END_REF] belongs to stochastic optimisation technique influenced by the social behaviour of bird flocking or fish schooling. The use of a PSO makes it possible to optimise the performance of the classifier LS-SVM method for clustering in the CRSN networks. Use of the PSO-LS-SVM methods permit to build a novel bybrid artificial intelligence (AI) system.

The main goal of this paper is to introduce a new hybrid AI system based on PSO-LSSVM techniques for clustering in CRSNs. In the proposed method, the structure of the CRSN is mapped into the network energy consumption model. Thus, a spectrum-aware clustering process is efficiently realised. We show that suggested method achieves good performance efficiency: the clusters are organized in such a way that the total communication power is minimised. This paper is organised as follows. In Section 2, we present an energy-efficient spectrum aware consumption model for CRSNs. Section 3 introduces the hybrid AI system based on the PSO and LS-SVM methods. In Section 4, we evaluate the effectiveness of the proposed method through intensive simulations. Section 5 concludes the paper.

Model of the network energy consumption

In this section, we present the network energy consumption model. We take into account data transmission, data signalling and data sensing in the CRSN.

The data transmission involves intracluster and intercluster communications. In the first case, all CRSN nodes belonging to the cluster send their data to the cluster heads [START_REF] Zhang | Distributed Spectrum-aware Clustering in Cognitive Radio Sensor Networks[END_REF]. Thus, the sum of transmission power of all cluster members is given by:

P intra = K k=1 N k i=1,i =j P i (n k i ) (1)
where N k is the number of CRSN nodes in the k-th cluster, 1 ≤ k ≤ K. P i is the power needed to reading data transmission from the CRSN node to the cluster head node. The distance between cluster members and their cluster head are short within the cluster and, assuming a free-space channel model as an appropriate channel model, we can rewrite the sum transmission power as:

P intra = C 0 P r K k=1 N k i=1,i =j d 2 (n k i , n k j ) (2)
where P r is the minimal received power required for the cluster head to correctly decode the transmitted information, C 0 is constant loss factor, and d(n k i , n k j ) is the Euclidean distance between the i-th and j-th nodes in the k-th cluster.

As the cluster head node, we can select the number of nodes with maximal residual energy. As proven by Zhang [START_REF] Zhang | Energy-efficient Spectrum-aware Clustering for Cognitive Radio Sensor Networks[END_REF], intracluster energy achieves the upper bound if each node will have an equal probability of becoming cluster head. Hence, the sum for intracluster power consumption is given by:

P intra = C 0 P r K k=1 N k j=1 1 N k N k i=1,i =j d 2 (n k i , n k j ) = 2C 0 P r K k=1 N k i=1 d 2 (n k i , center(k)) (3)
where center(k) = ( 1

N k N k i=1 x k i , 1 N k N k i=1 y k i )
is the center of the k-th cluster. The total network wide energy consumption is formulated as follows:

E network = P sen t sen + P sig t sig + (KC 0 P r d 2 max + 2C 0 P r K k=1 N k j=1 d 2 (n k i , center(k)))t data (4)
where P sen , P sig is the energy devoted to sensing and signalling, t sen , t sig , t data are the sensing, signalling, and data transmission times, respectively, d max is the transmission range of the CRSN node.

3 Hybrid AI System Based on the PSO and LS-SVM Methods for Clustering in CRSN System

Particle Swarm Optimization

The particle swarm optimization (PSO) is a swarm intelligence method that models social behaviour of organisms such as bird flocking and fish schooling [START_REF] Kennedy | Particle Swarm Optimization[END_REF] to guide swarms of particles toward the most promising regions of the search space. PSO is conceptually similar to the crossover operation used by genetic algorithms. The main difference of particle swarm optimization method from the evolutionary computing is that flying potential through hyperspace are accelerating toward "better" solutions, while in evolutionary computation schemes operate on potential solutions which are represented as locations in hyperspace [START_REF] Kennedy | Matching Algorithms to Problems, An Experimental Test of the Particle Swarm and Some Genetic Algorithms on the Multimodal Problem Generator[END_REF].

Each particle represents a candidate position (i.e., solution). A particle is considered as a point in a D-dimension space, and its status is characterized according to its position and velocity. Thus, the D-dimensional for the particle i at the t-th iteration can be represented as x t i = {x t i1 , x t i2 , . . . , x t iD }. The velocity (i.e., distance change) for particle i at iteration t can be defined by v t i = {v t i1 , v t i2 , . . . , v t iD }. Let p t i = {p t i1 , p t i2 , . . . , p t iD } be the best solution that particle i has obtained until iteration t and p t g = {p t g1 , p t g2 , . . . , p t gD } be the best solution obtained from p t i in the population at iteration t. To search for the optimal solution, each particle changes its velocity according to the cognition and social parts as follows:

V t id = V t-1 id + c 1 r 1 (P t id -x t id ) + c 2 r 2 (P t gd -x t id ), d = 1, 2, . . . D (5)
where c i indicates the cognition learning factor; c 2 indicates the social learning factor, and r 1 , r 2 are random numbers uniformly distributed in U (0, 1). Each particcle moves to a new potential solution based on the equation

X t+1 id = X t id + V t id , d = 1, 2, . . . , D (6)
The iteration is terminated if the number of iteration achieves the given of number of maximum iteration.

LS-SVM Classifiers

Least squares support vector machines (LS-SVM) [START_REF] Suykens | Least Squares Support Vector Machine Classifier[END_REF], [START_REF] Suykens | Sparse Approximation Using Least Squares Support Vector Machines[END_REF] are least squares versions of support vector machines (SVM) [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF], [START_REF] Vapnik | Statistical Learning Theory[END_REF], which are set of related supervised learning methods that analyse data and recognize patterns, and which are used for classification and regression analysis.

Let {(x i , y i )} N i=1 be a given training set, with the input data x i ∈ R n and output data y i ∈ R with class labels y i ∈ {-1, +1} and the linear classifier

y(x) = sign[w T x + b] (7)
When the data of the two classes are separable we have the original SVM classifier [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF], [START_REF] Vapnik | Statistical Learning Theory[END_REF] that satisfies the following conditions:

w T φ(x i ) + b ≥ +1 if y i = 1 w T φ(x i ) + b ≤ -1 if y i = -1 (8)
These two sets of inequalities can be combined into one single set as follows

y i [w T φ(x i ) + b] -1 ≥ 0, i = 1, 2. . . . , N (9)
where φ: R n → R m is the feature mapping the input space to a usually high dimensional feature space. The data points are linearly separable by a hyperplane defined by the pair (w ∈ R m , b ∈ R). Thus, the classification function is given by

f (x) = sign{w T φ(x) + b} (10)
Instead of estimating with the help of the feature map we work with a kernel function in the original space given by

K(x i , x j ) = φ(x i ) T • φ(x j ) (11)
In order to allow for the violation of Eq.( 9), we introduce slack variables ξ i such that

y i [w T φ(x i ) + b] ≥ 1 -ξ i , ξ i > 0, i = 1, 2, . . . , N (12)
The following minimization problem is accounted for as follows:

min w,b,ξ J(w, b, ξ) = 1 2 || w || 2 +C N i=1 ξ i subject to y i [w T φ(x i ) + b] ≥ 1 -ξ i , ξ i > 0, i = 1, 2, . . . , N, C > 0 ( 13
)
where C is a positive constant parameter used to control the tradeoff between the training error and the margin.

The dual problem of the system (13), obtained as a result of Karush-Kuhn-Tucker (KKT) condition [START_REF] Kuhn | Nonlinear Programming[END_REF], leads to a well-known convex quadratic programming (QP). The solution of the QP problem is slow for large vectors and it is difficult to implement in the on-line adaptive form. Therefore, a modified version of the SVM called the Least Squares SVM (LS-SVM) was proposed by Suykens et al [START_REF] Suykens | Least Squares Support Vector Machine Classifier[END_REF], [START_REF] Suykens | Least Squares Support Vector Machines[END_REF].

In the LS-SVM method, the following minimization problem is formulated min w,b,e J(w, b, e) = 1 2

w T w + γ 1 2 W k=1 e 2 k subject to y k [w T φ(x k ) + b] = 1 -e k , k = 1, 2, . . . , N (14)
The corresponding Lagrangian for Eq. ( 13) is given by

L(w, b, e; α) = J(w, b, e) - N k=1 α k {y k [w T φ(x k ) + b] -1 + e k } ( 15
)
where the α k are the Langrange multipliers. The optimality condition leads to the following (N + 1)

× (N + 1) linear system 0 Y T Y Ω * + γ -1 I b α = 0 1 ( 16
)
where

Z = [φ(x 1 ) T y 1 , . . . , φ(x N ) T y N ] Y = [y 1 , . . . , y N ] 1 = [1, . . . , 1] α = [α 1 , . . . , α N ]
and Ω * = ZZ T . Due to the application of Mercer's condition [START_REF] Mercer | Functions of positive and negative type and their connection with the theory of integral equations[END_REF] there exists a mapping and an expansion

Ω * kl = y k y l φ(x k ) T φ(x l ) = y k y l K(x k , x l ) (17)
Thus, the LS-SVM model for the function estimation is given

y(x) = N k=1 α k y k • K(x, x k ) + b (18)
where parameters α k and b are based on the solution to Eqs. ( 15) and ( 16).

Mixtures of kernels

Each kernel function is characterized by its advantages and disadvantages. For instance, a radial basis function (RBF) kernel K(x;

x i ) = exp{-| x -x i | 2 /σ 2 },
where σ is the width of the radial basis function, is a typical local kernel in which only the data that are close have an influence on the kernel values. The polynomial kernel (see [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]) such as K(x; x i ) = [x • x i + 1] q , where q is the kernel parameter which defines the degree of the polynomial to be used, guarantees the influence of all the data points that are far away from each other. Therefore, the mixture of these kernels gives a better performance than any single kernel.

As was defined by Smits et al. [START_REF] Smits | Improved SVM Regression Using Mixtures of Kernels[END_REF], an exemplary mixture of the RBF and polynomial kernels is given by

K mix = ρ • K poly + (1 -ρ)K RBF ( 19
)
where ρ is the mixing coefficient treated as a constant scalar.

LS-SVM transformed into a clustering problem in CSRN

The input data are described by a coordinate (r, z) and the output data is the energy value. The LS-SVM method can be transformed into a clustering problem by the use of the same set of input vectors but different sets of labels. We can transform the Eq. ( 16) into a system

0 1 T 1 Ω b α = 0 Y (20)
where

Y = [y 1 ; . . . , ; y N ], Ω = K(x i , y i ) + γ -1 I, 1 T = [1 1 ; . . . ; 1 N ], α = [α 1 ; . . . ; α N ].
Thus, the Ω is given by

Ω = K(x i , x j ) + γ -1 I (21)
The solution of Eq. ( 21) gives the values

   b = 1 T Ω -1 Y 1 T Ω -1 1 α = Ω -1 (Y -b1) (22)
By setting

A = Ω -1 and B = 1 T Ω -1 1 T Ω -1 1 we obtain b = BY α = A(Y -b1) (23)
where A and B are precalculated matrixes that depend only on the input vector (x k ) but not on the vector y k .

The sensors are usually correlated due to the high probability that the adjacent pixels will contain the sensors in the cluster. We assume that the sensor field is two-dimensional and the image energy distribution of the sensor field on the surface, known as the point spread funcion (PSF), can be approximated by the Gaussian PSF. On the other hand, the center point of the PSF corresponds to the measured sensor position.

In our approach based on the mixtures of kernels, we take into consideration two sets of indexes of the neighbourhood of the sensor image that satisfies the condition of linearly separable patterns. We recall that the linear separability

Experimental Results

An experiment to verify clustering efficiency of PSO-LS-SVM approach was performed for randomly deployed Pus and CRSN nodes in a 200 m x 200 m square area. We assumed that each PU randomly operates in four channels in which the neighbouring Pus are inactive. In the simulation, we supposed that the distance between a transmitting and receiving radio pair is normally distributed with a mean of 30 and variance of 15.

We used a link gain of h = h i /d, where d is the distance from transmitter to receiver. If the distance is was less than 10, we set link gain to one.

The proposed PSO-LS-SVM approach for clustering in CRSN is composed of two main phases. PSO algorithm was applied as a technique to improve the setting of the parameter values of LS-SVM. Clustering using LS-SVM method consists of two phases: training phase and calculate fitness phase. A block diagram of the algorithm is depicted in Fig. 1.

The PSO-LS-SVM algorithm was developed based on Matlab R2014b and LS-SVMlab 1.8 for the optimization of model. Through training, the parameter values of the PSO and LS-SVM approach were set as follows. The cognition learning factor c 1 and the social learning factor c 2 were set to 1.8. The number of particles and generations were found to be 8 and 150. Since the PSO-LS-SVM approach has a larger solution space, the number of solution evaluated was raised to 1000. The optimized parameters for LS-SVMs inlude kernel parameters and regularization parameter. If radial basis function (RBF) or polynomial (Pol) function is chosen as a kernel function, we used the vector ν = (γ, σ 1 , σ 2 , . . . , σ Ninput ) and ν = (γ, σ, δ), respectively, where γ is the regularization parameter, δ and σ i (i = 1, 2, . . . , N input ) are kernel parameters, N input is the dimension number of input patterns. To perform the clustering evaluation, six PUs are randomly deployed in our sensor field. Furthermore, the sum of the squared errors (SSE) of the clustering results for various numbers of clusters was taken into account for each scenario. In Fig. 2, we compare the average SSE with respect to CRSNs of different sizes for all of the studied number of clusters. The experiment shows that the method for a given CRSN has clustering effectiveness with the minimal average SSE, when the cluster number is equal to five or six. Fig. 3 plots the average cluster size versus the number of nodes in the CRSN and given number of active PUs. As we can see, all PUs behave similarly in terms of the number of clusters and the cluster size for different numbers of PUs for various CRSN sizes. However, the small number of PUs gives a greater average cluster size.

To compare the efficiency of LS-SVM and PSO with various kernels, we used the Rand index. Rand index measures the percentage of decisions that are correct, namely

RI = T P + T N T P + F P + F N + T N × 100% (28)
where T P is the number of true positives, T N is the number of true negatives, F P is the number of false positives, and F N is the number of false negatives.

Table 1 presents the values of Rand index for the PSO-LS-SVM approach under different circumstances. As can be seen from Table 1, the efficiency of PSO-LS-SVM with the mixture of kernels is better than other models.

Conclusion

In this paper, a new clustering method for cognitive radio sensor networks is introduced. The applied methodology used in clustering has been found in the network energy consumption model for CRSNs. Three parts of this model were used: multi-task sensing, signalling and data transmission. Using the energy distribution of the sensor field surface and applying the PSO-LS-SVM system, we have been able to solve the clustering problem in the CRSN. The presented hybrid artificial intelligence system is shown to satisfy CRSN energy constraints through extensive simulation experiments. Our PSO-LS-SVM method showed that the proposed solution offers several advantages with respect to classical clustering schemes, such as k-means clustering algorithm, etc. Moreover, this method incorporates the spectrum-sensing capability and spectrumaware sensing constraint in the nodes and their localisation in the network. The presented method achieves comparable effectiveness with the other clustering method in cognitive radio sensor networks. Future research directions include

  For any data block, the input points are in the form {(r o + dr, z o + dz) : | dr |< m, | dz |< n}. All such sets of points can be translated to the same set {(dr, dz) :| dr |< m, | dz |< n} by subtracting (r o , z o ) from all the vectors, m and n are the half numbers of the horizontal and vertical pixels of the blocks.

Fig. 1 .

 1 Fig. 1. Block diagram of PSO-LS-SVM-based algorithm.

Fig. 2 .

 2 Fig. 2. Average sum of the squared errors (SSE) versus size of the CRSN.

Fig. 3 .

 3 Fig. 3. Average cluster size versus size of the CRSN.

Table 1 .

 1 Efficiency comparison between various models obtained by PSO-LS-SVM.

	Model	RI (%)
	LS-SVM P oly	61.8
	LS-SVM RBF	65.4
	LS-SVM mix	67.5
	PSO-LS-SVM P oly	83.6
	PSO-LS-SVM RBF	88.2
	PSO-LS-SVM mix	89.3

requires that in order to be classified, the patterns must be separated from each other to ensure that decision surfaces should consist of hyperplans.

The LS-SVM with the RBF and polynomial kernels transformed into a sensor clustering problem has the fitted energy intensity surface function over the constant vector space as follows

where (r, z) are the coordinates of the pixels. The function g(r, z) gives the corresponding energy intensity value, b and α are obtained as a solution of Eq. (23).

A multi-class formulation of the LS-SVM transformed into a sensor clustering problem

In comparison with the standard SVM method, the LS-SVM has a lower computational complexity and memory requirements. Nevertheless, in certain situations, such as the classification of several characters, clusters, etc., a multi-class classifiction is very suitable.

In the multi-class formalism we now use multiple output values y i with i = 1, . . . , n y , where n y defines the number of output values [START_REF] Vapnik | Statistical Learning Theory[END_REF]. Thus, in the primal weight space the multi-class classification system possesses the following binary classifiers

with mappings on a high dimensional feature space φ (i) (.) : R n → R n h , i = 1, 2, . . . , n y , with dimensions n h1 , n h2 , . . . , n hn y .

By the extension of Eq. ( 23) to a multi-class problem, we obtain

where matrix Y is given by

2 , . . . , y

and vector b M = [b (1) , . . . , b (ny) ], where N is the number of pixels in one processing block. the development of clustering methods for CRSN with energy constraints by allowing each node to select its role without the need for cooperation between nodes.