Multiscale Texture Analysis and Color Coherence Vector Based Feature Descriptor for Multispectral Image Retrieval

Content Based Image Retrieval (CBIR) for remote sensing image data is a tedious process due to high resolution and complexity of image interpretation. Development of feature extraction technique is a major portion to represent the image content in an optimal way. In this paper, we propose a feature descriptor which combines the color coherent pixel information and GLCM texture features in multi scale domain. Curvelet transform is used to decompose the image into coarse and detail coe ﬃ cients. Then Gabor magnitude is computed for each coe ﬃ cient to improve the directional information. GLCM texture features are extracted from the Gabor magnitude response. The novel feature set by combining the CCV and GLCM using curvelet and Gabor ﬁlter is developed. Mahalanobis distance measure is used to ﬁnd the similarity between query and feature database. Average Normalized Modiﬁed Retrieval Rate (ANMRR) is computed to evaluate the performance with the state of art methods.


Introduction
Texture represent the pattern arrangement of the pixels in the image. It is an inherent property of any image. Especially in remote sensing applications, images represent the surface reflectance of the earth. Various sensors are used to capture the earth surface like Multi/Hyperspectral, RADAR and LiDAR [1]. Texture extraction in spatial domain using Local Binary Patterns (LBP) and Gray Level Co-Occurrence Matrix has been explained in [2,3,4]. Other local texture pattern features like Local Derivative Patterns (LDP), Local Tetra Patterns (LTrp) are reviewed by [5,6]. Recent advancements in texture features have been developed using multi-scale GLCM windows [7] and Local Contourlet Tetra Patterns [8]. Texture analysis in multi-scale and multi resolution analysis is very useful to extract more precise features. In proposed work, Curvelet transform has been used to analyze the input image in multi scale domain. Curvelet is directionally more efficient than other existing multi resolution transforms. It can preserve the tiny non linear edge information due to its anisotropic nature. Coarse and detail scale coefficients of curvelet are used to compute the Gabor magnitudes. Total four Gabor magnitude responses are be used to extract the GLCM properties. Applying the Gabor magnitudes on each scale of curvelet gives more directional features of the non linear edges. Texture information extracted from these coefficients can perform retrieval applications efficiently. To carry the spatial domain information, color coherent pixels are computed and merged to the texture feature descriptor. The proposed feature descriptor outperforms the state of art retrieval methods. It has been tested on the public remote sensing image retrieval dataset namely UCmerced. The paper has been organized as follows: Section 2 Related work, Section 3 Methodology,Section 4 Results & Evaluation and Section 5 Conclusion.

Importance of Texture Feature in Remote Sensing
In recent years, many texture-based feature descriptors have been developed for the purpose of CBIR in remote sensing. This section explaines comprehensive review of texture feature extraction in remote sensing image data. The importance of spatial information retrieval in theoretical perspective for model fit and model selection techniques has been explained by Mihai Datcu, et.al, in (1998) [9]. Seisuke Fukuda, et.al (1999) have developed a wavelet-based texture feature set for the polarimetric SAR image data, Where Downsampling of wavelet coefficients are omitted. The energy of the wavelet coefficients is considered as texture descriptor [10]. M. Chica-Olmo, et.al (2000) have proposed univariate and multivariate texture features by using spatial variability based variogram measures [11]. It will refine the spatial and spectral structure of the image. A multiband image will be given as an output to quantify the local spatial variability of radiometric data. Jean-Luc Starck, et.al (2002) have improved the curvilinear feature extraction in multiscale domain by using curvelet transform to overcome the wavelet drawbacks [12]. Second and higher order statistical measures have been applied to get the texture information in multiscale and multi resolution domains. Andrey V. Bogdanov, et.al (2005) have applied neural network based classification on GLCM features for the classification of the sea ice in multi-sensor data [13]. F. Kayitakire (2006), et.al, have derived Forest variables from the high spatial resolution data by using GLCM [14]. Interactive retrieval system is proposed by Marin Ferecatu, et.al (2007), using relevance feedback method to improve the accuracy [15].
Feature level fusion also improves classification accuracy. Vijaya V, et.al (2009) had developed a system for texture feature classification by fusing statistical and wavelet-based texture features [16]. Changren Zhu, et.al (2010) have proposed shape and texture based feature set for detecting the ship in spaceborne optical images. Compactness, Convexness, Rectangularity, and Eccentricity are used for the shape based feature set and Mean, Standard variance and Entropy are used for analyzing the statistical features. Local Multiple Patterns are used as local texture features [17]. Shufu Xie, et.al (2010) had developed a Local Gabor Xor Patterns (LGXP) for facial recognition [18]. Multi-resolution based feature extraction has succeeded in extracting robust features. Gholamreza Akbarizadeh, et.al (2012) have proposed Kurtosis wavelet energy feature for the texture recognition in SAR images [19]. 3D DWT is suitable for extracting structure and texture information from the Hyperspectral data. Hyperspectral image consists of continuous spectral bands that provide more accurate information. Yuntao Qian, et.al (2012) has proposed a 3D wavelet-based texture feature extraction and classified it using sparse linear regression [20]. Simple statistics, Gabor filter bank and color histogram has been used as feature descriptors for efficient image retrieval by Yi Yang, et.al (2012) [21]. Local spectral histogram features are generated by merging the local histograms. Jiangye Yuan, et.al, (2014) have computed local spectral histograms and Gabor filter based texture features to segment the image regions [22]. Morphological texture features based on structuring elements have also shown good results in CBIR techniques. Erosion, Dilation, Open, Close operations by structuring elements will improve the feature extraction efficiency. Erchan Aptoula (2014) have developed a morphological texture descriptor based on Circular Covariance Histogram (CCH) and Rotation Invariant Triplet (RIT) for Land Cover image retrieval [23]. Alaa Al-Hammami and Hisham Al-Rashdan (2010) have proposed an enhancement to the Color Coherence Vector (CCV) based on modification of the distance and angle [24]. Zhen-feng Shao, et.al have extracted Gabor filter based texture and color features for retrieval application. Unichrome feature of each channel of the image data using Gabor filter with orientation (u, v) and opponent feature with (u, v') are used as color features respectively [25]. GLCM is the best features to extract texture information of mono channel image data. For Multichannel or Hyperspectral image, Data clustering and sparse representation based GLCM has been proposed by Xin Huang, et. al (2014) [26]. Local Binary Patterns (LBP), Local Ternary Patterns (LTP), Local Terra Patterns (LTrP) and Local Derivative Patterns (LDP) are most popular texture extraction methods evaluated by [27,28,29]. Color and Texture descriptor based image retrieval system has been developed by using Block Difference Inverse Probability (BDIP) and Block Variation of Local Correlation Coefficients (BVLC) by Chandan Singh, et.al (2016) [30]. Spectral end members and texture features based on GLCM has been developed for hyperspectral image retrieval by Jing Zhang, et.al (2017) [31].

Efficiency of Curvelet Transform
The wavelet-based approach have failed to preserve nonlinear singularities due to the directional limitation. Anisotropic nature of curvelet transform helps in extracting curved singularities. It is a multiscale and multi-resolution transform which can process an image in several orientations at each scale. Curvelet-based SAR image classification has been performed by using the Histogram of Curvelet (HOC) [32]. The curvelet transform has been developed by Emmanuel J. Candes, David L. Donoho, et al. (2002) in order to overcome the drawbacks of wavelet [33]. Lindsay Semler et al. (2006) highlighted the curvelet-based texture descriptors over CT scan images. Energy, entropy, mean, standard deviation has been calculated by using the effective curvelet-based texture descriptors. The results have been compared with similar algorithms based on wavelet and ridgelet descriptors [34]. Tobias Gebacket al. (2009) have explained that when the number of coefficients that are required to represent the edge to a specified accuracy is regarded, curvelets dispense an almost ideal representation of C 2 edges. Decomposition of the image f into J feature levels have been furnished by the discrete curvelet transforms, with L j directions on each level, and K jl,1 × K jl,2 spatial shifts for each of these directions. [35]. Jianwei ma and Gerlindplonka (2010) have described that the compression based on a wavelet, removal of noise or structure extraction have become computationally inefficient for geometric features with line and surface singularities. The first-generation curvelets have been used for image denoising. The first-generation curvelets are also been used to intensify image contrast, render astronomical images and integrate satellite images.After introduction of the second-generation curvelets in 2004, the applications of curvelets escalated in several fields. Hybrid methods are more popular that are pertaining to discrete curvelet transform, which is formulated by incorporating curvelets with another technique for image processing [36] The experimental outcome shows that row-planted coffee field can be sectioned in a satellite image with the conjunction of structural and density features. The yield of coffee production can be determined from the data on physical area and density, which is provided by these features. [45].

Methodology
Remote sensing benchmark multispectral dataset UC Merced is considered for the experiment in the proposed work. The dataset is consisting of 21 classes of land use data for research purpose [46]. The images are subsampled from USGS National Map urban area image collection. Each image consists of 256 x 256 pixels. The resolution of each pixel is 1 foot. The input image bands will be separated to apply the 2D Fast Discrete Curvelet Transform in Unequally Spaced Fast Fourier Transform (USFFT). Curvelet transform will decompose the image into three coefficients which are coarse, detail, and fine scale coefficients. The detail scale coefficients again consist of multiple orientations at each resolution [47]. Figure. 1 shows the process of retrieval system.

Probabilistic Principle Component Analysis (PPCA)
PCA is a popular method for dimensionality reduction. It will increase the variance of the projected space. The high dimensional strongly correlated data {D p } p P will be converted into low dimensional uncorrelated data {D l } l L removing low variance data. This may lead to some data loss due to reducing the dimension to (P-L) Eigen subspace. PPCA is a method of finding principal axes by the probabilistic way [48]. Each variable of this model is calibrated by a specific probability density function (pdf). The uncorrelated input variable z with some additive noise e is represented by the equation form x=A.z+e. Where A is model parameter. If some probability densities are set to the generative input variable then it can be called as probabilistic model. PPCA is used to find the maximum likelihood parameter set Θ ≡ A, λ. Where A R (p×L) is www.astesj.com 272 Figure 1: Architecture of the proposed system the loading the matrix, z is pdf of latent variable and e is pdf of additive noise variable. The parameter Θ can be measured by Expectation and Maximization (EM) algorithm [49]. PPCA assumes that p(z) = ℵ(z : 0, I) and p(e) = ℵ(e : 0, λ.I). Gaussian pdf is similar to linear operations so that the linear transformation of the Gaussian will produce the new Gaussian distribution.

Curvelet Decomposition
The directionality and anisotropic natures are the main features of the curvelet transform which can help to extract more information from the image data. The curvelet transform decomposes the image into three components like coarse, detail and fine scale coefficients [50]. The input image I will be decomposed with the help of low pass and bandpass filters l 0 , ∆s respectively.
The default number of scales will be log 2(n) -3. Where n is length of the image. Curvelet coefficients are denoted as c {j} {l} {k} = f , ϕ j,l,k . Where j represents scale, l represents orientation and k represents the position (k 1 , k 2 ). The curvelet coefficients can be stored in a cell matrix C j,l . Here scale j is from coarse to finest and orientation l is starts from top left corner and rotates in clockwise [51]. The rotation angles θ j,l will be equidistant and sequence can be  15 and eight translations with different grids will be considered. The grid sizes will be 1/2 −j × 1/2 −j/2 . So the anisotropic scaling principle of the curvelet will be width = length2. Each subband is partitioned into squares of proper scales as, ∆ s I → (W Q ∆ s I) Q Q s . Where W Q is smooth window. Assume that window WQ (n1, n2) is supported within the sheared rectangle p j = (n 1 , n 2 ) : 0 ≤ n 1 − n 0 < L ( j), −L j /2 ≤ n 2 < L j /2 The localized dyadic squares Q = k 1 2 s , . By multiplying the function with corresponding window function will give the localized result near Q. Let (T Q I) denotes the operator which can transport and renormalizes I so the part of the input supported near Q will become the part of the output supported near to [0, 1] 2 . T Q I(n 1 , n 2 ) = 2 s I(2 s n 1 −k 1 , 2 s n 2 −k 2 ).The outcomes of dyadic squares are renormalized into single scale. g Q = (T Q ) ( − 1)(W Q ∆ s I), Q Q s . Each standardized dyadic square will be analyzed by ridgelet transform. It has basis element ρ λ that make an orthonormal basis for l 2 (R 2 ).α µ = (g Q , ρ λ )whereµ = (Q, λ).

Color Coherence Vector
Color histograms are the most popular features to extract color descriptors for CBIR. But the spatial information will not cover by the color histogram features. To consider the color and spatial information of the image Color Coherence Vector (CCV) [53,52] descriptor has been used. The first stage to compute the CCV is blurring the image with Gaussian filter. Then the discretization of the color space will be applied to divide into n distinct colors of the image. The input image has www.astesj.com 3 channels, discretization function will find unique colors in each channel to form V unique values of each channel. Table. ?? shows the pseudo code for computing CCV feature vector.

Gabor Magnitude
Magnitude response of image will be calculated by applying Gabor transform for an input image in multiple scales and different orientations [54]. Bandlimited filters can be used to represent optimal localization in both spatial and frequency domain. If the amplitude of the frequency spectrum goes to zero which exceeds the threshold frequency value. The Gabor transform can also be defined as a Short Time Fourier Transform (STFT) with a Gaussian kernel. The 2D Gabor filter can be mathematically defined as where u = ucosθ y + vsinθ y , v = −ucosθ y + vsinθ y , f x = f max 2 (x/2) , θ y = yπ/8. f x and θ y Parameters are the center frequency and orientation of the plane wave for the Gaussian kernel. The variables γ is the ratio of the center frequency and η is the sharpness constant. If the γ and η values are fixed then the center frequency f u will define the scale of the Gabor filter. Gabor filters are complex type, which is a combination of sine and cosine functions. Gabor filter is applied to the input gray image with wavelength 4 and 90-degree orientation. Wavelength denotes the sinusoidal carrier in pixels per cycle. Orientation defines the direction of the plane wave range in 0 to 360 degree. Convolution of the Gabor filter and input image will produce the resultant image in the time-frequency domain. It can be decomposed into real and imaginary parts [55]. The mathematical representations of the convolution and decomposition as follows where ) is the transformed image obtained after convoluting with Gabor filter ψ x,y (u, v). u, v are the spatial coordinates of the image. Where I is the input image. E x,y (u, v) and O x,y (u, v) are the even and odd functions for the sin and cosine functions and P is the phase response of the sinusoid carrier. From the decomposed sine and cosine functions of the Gabor filter, Magnitude and Phase responses will be calculated as follows The magnitude responses for the curvelet coefficients are considered for the extraction of Haralick texture features. The adjoint transform of the curvelet coefficients will be computed for each scale j and orientation l by applying 2D Inverse Fast Fourier Transform (IFFT).
g j,l [n 1 , n 2 ] = F(C j,l,k ) (8) The adjoint transform will be computed for coarse, detail and fine scale coefficients separately. Gabor transform will be applied on the resultant image at each scale to compute the magnitude response. The surface view of the magnitude response at each scale has shown in Figure. 2.

Texture Analysis
To identify local or global pixel patterns, texture features are very important. GLCM plays a vital role in texture analysis. Haralick, et.al, (1973) derived statistical texture features based on GLCM. GLCM matrix consists value i with a sum of the frequency of the pixel value at a pixel value u with a pixel value y in a specific spatial orientation at v. Mathematically it can be represented as follows The texture is an inherent property for all kind of surfaces. It can give important information about the pixel arrangement structure of each surface and spatial relationship. The smallest differences in texture for large image databases is very difficult. GLCM properties will be helpful to detect the texture features in local and global areas. Some of the GLCM properties derived by Haralick [4] are considered for the experiment.

Similarity Measure
In this work, the Mahalanobis distance measure is used as a distance metric. Suppose if there are two groups G 1 and G 2 with some characteristics. these characteristics are considered as the feature vectors representing the objects of each group www.astesj.com 274  [56]. the feature vector x with dimension p have the same variation about its mean within either group. The mathematical representation of Mahalanobis distance is as follows where µ 1 andµ 2 are the mean of two vectors S is covariance matrix.

Results and Evaluation
The features have been extracted using the process explained in the section2. The experiment has been conducted on UC Merced Land use dataset which contains 21 classes with 50 samples of each class. Each image has 256 × 256 pixels. The spatial resolution of the data is 0.3m. The feature descriptor size for each image is 1 × 24. color coherence vector feature length is 1 × 8 because in this experiment only coherent pixel information is considered as a feature. The color coherence vector feature is computed for the multispectral image directly. Further feature extraction process the input image will be reduced to a single band by applying PPCA. The resulting of PPCA image will be decomposed into coarse and detail coefficients using curvelet transform. The coarse coefficient in one scale and detail coefficients at three different scales are used to extract the texture information. Obtained Gabor magnitude response of each coefficient to compute the GLCM properties. Contrast, Correlation, Energy and Homogeneity features have been extracted from each coefficient. The feature length of texture features is 1 × 16. The color and texture features will be combined together to form a single feature descriptor of an image. Feature descriptor of each image in the database will be computed and stored as a feature base. While the query image submitted to the retrieval system the same features will be extracted to describe the content of the image and matched with the feature base using Mahalanobis similarity metric. The feature obtained lower distance is considered as the most relevant image to the query. In this experiment, top 10 relevant images are shown as a result. Figure.  sures will not considered the position of the relevant image appears among the top retrieved results. The retrieval system performance is evaluated by using Average Normalized Modified Retrieval Rate (ANMRR) metric [57,58]. ANMRR The performance metric ANMRR as shown below where R(I) is position of relevant image I retrieved, I(q) is retrieved images for query q.
, ∀q]} Average Rank will be calculated as follows where TG(q) is the total ground truth images in the database to the query q. TG(q) is influencing the average rank in the above equation. to reduce this factor Modified Retrieval rate (MRR) is computed as Aerage R(q) − 0.5.[1 + T G(q)]. Normalization will be applied to make the value exist between 0 and 1.
Finally the ANMRR is derived for the Total number of queries TQ is follows Table 2 shows the ANMRR values for the sample retrieval results with ten queries from each class. we have tested for the four classes in the database. Retrieval performance was compared with the five existing popular texture feature extraction methods. The proposed method is developed by combining the CCV and texture features using GLCM from the Gabor magnitude of the curvelet sub-bands in different scales. Extracting Gabor magnitude from the curvelet sub-bands in multiple scales will improve the directional features of the input image. Curvelet is capable of extracting anisotropic edge information in various orientations at different scales. The resultant of the curvelet decomposed image will have one coarse and three detail coefficients as explained in the section 3. The comparison of the proposed method with existing literature has been shown in Figure 4.

Conclusion
This paper develops a feature extraction method which is a combination of color coherent pixel information and multi   scale texture features. The proposed feature descriptor is extracted using Haralick texture features from curvelet coefficients. To extract the robust directional infornamtion, Gabor magnitude responses are computed from the curvelet coefficients. Haralick texture features extracted from the resultant gabor magnitude responses are shown best retrieval rate than other state of art feature extraction techniques. It improves the retrieval accuracy of the CBIR search engines. Retrieval system was tested on various classes of remote sensing image data. ANMRR performance metric gives the best average retrieval rate for a set of queries from various classes. Four classes have been tested in the dataset. This work will be enhance by adopting machine learning techniques for the target classification for high resolution imagery data.