# Characterization of Complex Image Spatial Structures Based on Symmetrical Weibull Distribution Model for Texture Pattern Classification

Authors: Jinping Liu; Jiezhou He; Zhaohui Tang; Pengfei Xu; Wuxia Zhang; Weihua Gui
Journal: Complexity
 (2018)

---

## Abstract

Texture pattern classification has long been an essential issue in computer vision (CV). However, texture is a kind of perceptual concept of human beings in scene observation or content understanding, which cannot be defined or described clearly in CV. Visually, the visual appearance of the complex spatial structure (CSS) of texture pattern (TP) generally depends on the random organization (or layout) of local homogeneous fragments (LHFs) in the imaged surface. Hence, it is essential to investigate the latent statistical distribution (LSD) behavior of LHFs for distinctive CSS feature characterization to achieve good classification performance. This work presents an image statistical modeling-based TP identification (ISM-TPI) method. It firstly makes a theoretical explanation of the Weibull distribution (WD) behavior of the LHFs of the imaged surface in the imaging process based on the sequential fragmentation theory (SFT), which consequently derives a symmetrical WD model (SWDM) to characterize the LSD of the TP’s SS. Multidirectional and multiscale TP features are then characterized by the SWDM parameters based on the oriented differential operators; in other words, texture images are convolved with multiscale and multidirectional Gaussian derivative filters (GDFs), including the steerable isotropic GDFs (SIGDFs) and the oriented anisotropic GDFs (OAGDFs), for the omnidirectional and multiscale SS detail exhibition with low computational complexity. Finally, SWDM-based TP feature parameters, demonstrated to be directly related to the human vision perception system with significant physical perception meaning, are extracted and used to TP classification with a partial least squares-discriminant analysis- (PLS-DA-) based classifier. The effectiveness of the proposed ISM-TPI method is verified by extensive experiments on three texture image databases. The classification results demonstrate the superiority of the proposed methods over several state-of-the-art TP classification methods.


---

## Body

## 1. Introduction

Texture is an important visual perception cue ubiquitously existing in almost all-natural images. TP classification (TPC) plays an essential role in a variety of image analysis and computer vision (CV) tasks, such as image segmentation [1], scene classification, object recognition and image retrieval [2], and many others more which cannot be enumerated exhaustively [3]. Developing a TPC or identification system mainly involves two steps, designing effective algorithm for TP feature extraction [2] and training a classifier for TPC [4, 5].A great many texture feature extraction methods can be found in the literature reviews, e.g., the early comparative study work [6], a more recent review of the texture descriptors [7], and a specific survey of the texture feature extraction methods for medical image analysis [8]. However, it is worth noting thattexture is a kind of subjective and abstract visual concept of human beings, including multiple perceptual properties of the observed object surface [9], such as coarseness, roughness, regularity, directionality, contrast, etc. but cannot be defined or described effectively in the CV. Hence, the TPC is still a challenging task due to the unclear definition of the concept TP, the wide ranges of natural texture types, and the large intraclass variations in TPs [10]. To facilitate CV application or digital image processing and analysis, researchers have proposed many approximation methods to describe these perceptions depending upon their specific definition of TP.There has much interest in use of local (invariant) descriptors for texture analysis in CV tasks, such as region segmentation, texture pattern classification, and even the real applications, e.g., local descriptors-based face detection or recognition [12]. These local structure descriptors can be mainly divided into two types, the sparse descriptor based on the detectable interest points, e.g., the scale invariant feature transform (SIFT) [13], descriptor and its variants [13, 14], and the dense descriptor based on the whole image pixel by pixel, e.g., the local binary pattern (LBP) descriptor [15], the pairwise rotation invariant cooccurrence local binary pattern (PRI-CoLBP) [16], and Weber Local Descriptor (WLD) [12]. In addition, some researchers combined both the ideals of both the sparse and the dense descriptors in a crossway. For example, Li et al. used SIFT in a dense sampling way [17] and Heikkilä et al. extracted the LBP base on the interest regions in a sparse way [18]. All these methods attempt to achieve discriminant features from local image patches by choosing a limited subset of the texture features [3] in the original pixel space or in the filtering response space. Even though there has been ample empirical evidence that these methods can achieve good classification performance on some specific texture databases, it concluded in a comparative study [19] that “no single approach did perform best or very close to the best for all images, thus no single approach may be selected as the clear “winner” ”for TPC.Visually, TPs, especially the fully microtextons-occupied complex TPs, are highly structured with the spatial configurations of LHFs. In other words, real-world TPs are distinct but share some random-distributed elements [20]. Randomly configured or distributed LHFs or the frequency of occurrence of LHFs connects the complexity of TPs [21, 22], which brings about the visual varieties of the SS of TPs. Hence, how to characterize the distinctive organizations of the TP’s SS is greatly important for TPC. This paper mainly concerns the characterization of the essential structural feature of the image surface appearance for TPC.Since the visual appearance of a region in the imaged TP is determined by the random layouts of LHFs, the distinctive TP characteristics are inevitably related to the statistical methods [23]. The early widespread methods are the local or global statistics of the images pixel intensities or the image filter responses, such as the first-order and second-order statistics, and cooccurrence matrixes-based statistics, e.g., gray level cooccurrence or difference histograms [24], gray level cooccurrence matrix (GLCM), gray level run length matrix (GLRM) [25], and multivariate image analysis (MIA) [26], as well as their variants [27]. It demonstrated that the statistical measures-based texture research method was a promising tool in computer texture analysis. As stated by Varma and Zisserman in [28], the success of the canonical texture research areas, such as synthesis, classification, segmentation, and compression, was largely due to learning a statistical representation of the filter bank responses or the image pixel intensity values.The statistics-based texture analysis methods directly compute some local or global statistics (e.g., low-order moments) as the texture features in the filter bank response domain or in the original pixel intensity domain, whereas they do not give any prior assumption of the latent statistical distribution (LSD) of the SS of the TP. There is evidence that statistics-based texture features cannot achieve distinctive description for some complex and confusing TP analysis, which may even achieve the same measures from distinct TPs. For instance, it could obtain the same mean and variance in high probability from two types of texture images sampling from different distribution models [21], seen from Figure 1, where it displays two artificial images with the same statistics (zero means and unit variance), but exhibits distinctly different visual appearance. Hence, these statistical measures may result in low classification accuracies, since they cannot offer meaningful description of the LSD of the SS characteristics of texture images.Figure 1
Four kinds of different artificial texture images sampling from different LSD model but both with zero means and unit variance. (a) An artificial black and white stripe texture image sampling from a uniform distribution model. (b) A Gaussian distributed 2D random texture image. (c) A half-normal distributed noise image. (d) A fractal texture image with the fractal dimension of 0.8.
(a) (b) (c) (d)Some researchers attempt to discover the relation between the visual appearance and the LSD of the texture image’s SS in biovisual perception experiments or computational textural modeling methods [23, 29]. Considerable efforts have devoted to the statistical distribution modeling- (SDM-) based methods to the LSD characterization to learn a fuller representation of the distribution models based on the random field constraint, especially integrated with the prevalent multichannel and multiscale image analysis, such as dyadic Wavelet transform (DWT) and Gabor Wavelet transform (GWT) [30]. Many useful statistical models are introduced to approximate the LSDs of TPs in the filter bank response domains [31, 32].In particular, researchers have recommended some sharp-peaked with fat and long-tailed or skew-distributed models, e.g., Gaussian scale mixture model (GSMM) [33], Mixture Gaussian model (MGM) [34], Gamma distribution model (GDM) [30], Student’s t distribution model (STDM) [35], etc., to characterize the marginal statistical distribution (MSD) of texture images in the multiscale image transform domain, owing to the sparse or skew behavior of the image filter bank responses (IFBRs). In other words, the MSDs of the IFBRs are highly kurtotic in the center of the distribution but long-tailed on both sides of the distribution or exhibit seriously left-skewed or right-skewed distribution behavior.MSD models can be extended to multivariate statistical distribution (MvSD) models, e.g., the joint distribution, representing the statistical dependence of IFBRs or image pixel intensity values in the specific distance and fixed orientation, which were subsequently investigated by augmenting the aforementioned MSD model arguments. Besides the applications of image statistical modeling in TPC, MSD, and MvSD models can be also used as a prior to improve the capabilities of the image processing and analysis technologies, such as image denoising [36], image retrieval [37], object recognition [38], image segmentation [39], etc.SDM-based TP feature description methods assume the image pixels as samples from an empirical distribution model (EDM) with some unknown model parameters, which can be used to describe the relationship between a certain image pixel and its adjacent pixels by using a certain probability model, e.g., Markov random field (MRF) model and simultaneous autoregression (SAR) model [40]. Although the SDM-based method provides a meaningful approach for the LSD characterization, current methods have taken few analyses of the exact statistical distribution of the SS of TP in theory; namely, the MSD and MvSD models selection for TP feature description depends on the researchers’ experience and greatly related to the training samples. It would achieve poor classification performance provided the unreasonable empirical model is introduced. In addition, the commonly used MSD and MvSD models do not take good use of human visual perception (HVP) characteristics, which lacks of physical perception meaning. In other words, these methods are not based on the mechanism of the human visual system, and therefore they can hardly represent multiple textural properties salient to the visual perception [9], and they need a lot of image samples for the reasonable statistical model learning. Hence, common SDM-based texture analysis methods are basically TP-specific [41] and they would have poor performance on the classification of the TPs with large intraclass variability and great interclass similarity.To obtain distinctive SS feature of TPs and achieve high classification accuracy, this paper investigates a novel image statistical modeling-based texture pattern identification (ISM-TPI) method for TPC. Firstly, it made a theoretical analysis of the LSD behavior of the SS of TP. It reveals that the spatial structure of LHFs in the image surface obeys a typical Weibull distribution (WD) by introducing the sequential fragmentation theory (SFT). A symmetrical Weibull distribution (SWD) model is derived based on the SFT and the SWD model parameters of the multidirectional and multiscale SS of TP are extracted by the local differential operators, including a kind of fast steerable isotropic Gaussian derivative filters and a kind of oriented anisotropic Gaussian derivative filters. And then, TPs can be categorized with high accuracy based on the SWDM-based features with a partial least square-discriminant analysis classifier. The effectiveness of the proposed ISM-TPI method is verified by extensive experiments on three different texture image databases, which achieve relatively superior performance compared with the commonly used determination methods.


## 2. Statistical Modeling of CSSs in Texture Patterns

### 2.1. Theoretical Analysis of the WD Behavior of the LSD of Randomly Arranged Fragments

During visual imaging or scene observation process, local homogenous regions or objects in the imaging scene will fragment the field of view (FOV) into separated parts, e.g., the inside-object regions and outside-object regions [42]. If we enhance the imaging resolution, these inside-object regions and outside-object regions will be fragmented again by the presence of more LHFs or more complicated objects and finer details, and so on, until in any direction, the limitation of the FOV is reached. When the FOV is sufficiently fragmented by LHFs, the imaging scene becomes stable and the addition of any new details will not significantly change its appearance and consequently it leads to a specific visual appearance of TP. On the contrary, when reducing the resolution of the visual sensor, the adjacent local homogeneous structures will recombine together and constitute relatively coarser fragments. Geusebroek [42, 43] has demonstrated that this imaging process that fragmenting the visual surface details into LHFs can be simulated or analogized by the sequential fragmentation (SF) process well known in ore milling process, which is used to describe the distribution of the ore-particle size or mass in the ball grinding mill [44]. Hence, the distribution of the local fragments or particles in images can be predicted by the SF theory (SFT).In terms of SFT [45], the statistical distribution of LHFs based on the fragmentation process obeys the power-law distribution, which can be expressed by the following formula:(1)Fx′→x=xβλ-1where x′→x represents a single fragmentation process from a coarse fragment structure x′ to a fine fragment structure x and parameter β in the ore grinding process represents the average mass (or size given the ore density) of the crushed ore particles, whereas in the texture image analysis it denotes the average illumination intensity of the LHF x; λ is a shape parameter, satisfying 1≥λ≥0.With the increasing of the resolution power of visual sensor, more local fragment details are introduced into the imaging scene, which will divide the image into smaller LHFs (regions or particles). Hence, the histogram distribution of the illumination contrast is the integrating over the power-law-distributed LHFs [43], given by(2)nx=c∫x∞nx′Fx′→xdx′Substitute (1) into (2), and let c=1/β; we can obtain nx by solving the following equation:(3)nx=xβλ-1∫x∞nx′dx′βBy solving (3), it can be found that nx obeys the typical WD law; namely,(4)nx=NTxβλ-1e-1/λx/βλwhere NT=1/∫-∞∞nxdx is a normalized constant to make sure that nx is a reasonable distribution function.Since the resolution power of the visual sensor is limited, the fragmentation process will inevitably terminate and exhibit a stable TP; the statistical distribution of the SS of LHFs corresponds to the accumulated distribution of the fragments whose local contrast larger thanx, which is demonstrated obeying a symmetrical Weibull distribution model (SWDM) given by the following formula:(5)N>x=∫x∞nx′dx′=Ce-1/λx/βλTo obtain the generalization characterization of TPs, the SWDM can be modified a little by introducing a location parameterμ, indicating the distribution location (center), namely, the probability density of the SS of TPs can be described by a triparameter SWDM [25, 35, 46],(6)fx∣μ,β,λ=Ce-1/λx-μ/βλwhere C=1/∫-∞∞e-1/λx-μ/βλdx=λ/2λ1/λβΓ1/λ is a normalized constant and Γx=∫0∞tx-1e-tdt is the Gamma function.


### 2.2. Relations to the HVP System

It has demonstrated that SWDM can completely characterize the spatial layout of the stochastic textures [43]. In the SWDM depicted in (6), μ is a location parameter, representing the global reflectance of texture image, and it can be derived from the shape-from-shading method [47]; β is a scale parameter, directly related to the illumination contrast of the LHFs; λ is the shape parameter, reflecting the size of LHFs in texture image.The visual appearance of TP varies with the statistical distribution of LHFs. The psychophysical evidence indicates that two TPs are often difficult to discriminate when they produce similar distribution of responses of a bank of orientation and spatial-frequency selective linear filters [9]. From regular texture to pure random TP, the distribution of the SS of TP gradually changes from a power-law distribution to an exponential distribution and then to a typical Gaussian distribution. SWDM is a generalized model for the effective expression of LSDs of various forms by adjusting the model parameters [25, 35]. For instance, when λ=1, a SWDM will turns into an exponential distribution with the mean of β, when λ≃3.6, it is basically close to the Gaussian distribution. In addition, some studies have demonstrated that the shape parameter λ is directly related to the fractal dimension Df of a texture image. The fractal dimension can be approximated by Df=-3λ [46].Hence, it reveals that SWDM generate a complete orthogonal basis for stochastic textures description and SWDM parameters are directly related to the HVP properties [43, 46]. The commonly used HVP properties for texture perception includecoarseness,regularity,contrast,roughness, anddirectionality [9]. SWDM parameters can be used to make a psychophysical explanation of these perceptual properties [43], namely,coarseness,regularity,contrast,roughness, anddirectionality.Coarseness is a fragment size-related perceptual property. Generally, the larger the size of basic fragments in texture pattern, the coarser it is felt of the HVP system. In SWDM,coarseness can be reflected by the shape parameter λ and the smaller λ means the coarser TP. However,Coarseness is an observation scale-related parameter [51]. As a psychophysical evidence, if we magnify a texture image but without resolution enhancement of visual sensor, the HVP system will adjust the observation perception and keep the scale invariance with no more extra details addition in the FOV. SWDM is demonstrated to be well fits this scale invariance. In other words, the shape parameter λ remains unchanged under this circumstance. Alternatively, if the texture image is magnified with an increase of the resolving power, the shape parameter λ becomes smaller for the perception ofcoarser TP.Regularity is a layout variation-related HVP characteristic. In general, a finer TP is felt to be more regular. As analyzed above,coarseness (orfineness) can be reflected by the shape parameter λ; henceregularitycan be also expressed by the shape parameter λ. In the extreme case, a few particles or even just one fragment exists in the FOV, the LSD of the SS of TP rejects the SWDM because the real distribution obeys either a power-law distribution or a pure regular fragment (just one or two fragments fully occupy the entire FOV) and the exhibited statistical distribution is often multimodal. This phenomenon can be reflected with the over-fitting of the maximum likelihood estimation (MLE) of the shape parameter, e.g., λ≫2 [43].Contrast is illumination variation-related HVP characteristics [51]. It results from the illumination intensity and the surface height variations of the observation objects [47]. The global illumination variation is reflected by the location parameter μ, which indicates the inhomogeneous illumination surface, whereas the surface height variations can be reflected by the scale parameter β, thewidth of the SWDM according to the SFT. Hence, the perceptualcontrast can be indicated by two SWDM parameters, μ and β [43].Roughness is originally a tactile property in surface perception. Though it seems like a three-dimensional perceptual property, human beings can perceive it based on a two-dimensional image.Roughness is a HVP property results from the height variations of a certain homogeneous fragment. The greater the height variations of the fragments in the FOV, therougher perception it is felt. As analyzed above, β is an indicator of the height variations of TP, and the shape parameter λ is the reflectance ofcoarseness of the fragments in texture image. Hence,roughness can be reflected by the shape parameter λ and scale parameter β, under a special illumination condition, indicated by the location parameter μ [43].Directionality is a global sense of the whole FOV and caused by the shapes of the basis texture elements as well as their placement rules. Though the SWDM parameters do not contain the direct shape information, it reveals that the anisotropy of grains can be described with the dominant direction of the shape parameter λ[43–45]. Anisotropy in texture shadows (or contrast) can be reflected by the dominant orientation on the scale parameter β. If we fully consider the overall SS information of TP, the SWDM parameters can effectively describe the HVP characteristicdirectionality with the corresponding dominant orientation information.


### 2.3. Parameter Estimation of SWDM

The commonly used model parameter estimation method is maximum likelihood estimation (MLE). Given a sample setX=x1,x2,…,xN samples from a SWDM, the corresponding logarithmic likelihood function (LLF) log⁡LX∣μ,β,λ indicates how well the statistical model fit the training samples X. We can achieve the model parameter by maximizing LLF; namely,(7)μ^,β^,λ^=arg⁡maxμ,β,λ⁡log⁡LX∣μ,β,λ=arg⁡maxμ,β,λ⁡log⁡∏i=1Nλ2λ1/λβΓ1/λexp⁡-1λxi-μβλAfter expanding and simplifying the LLF, we obtain(8)log⁡LX∣μ,β,λ=Nlog⁡λ-1λlog⁡2λ-log⁡β-log⁡Γ1λ-∑i=1N1λxi-μβλThe SWDM parameters can be estimated by computing the partial derivative oflog⁡LX∣μ,β,λ with respect to μ,β,λ and setting each of the derivative to be equal to zero, respectively. Hence, we can obtain the following equation set:(9)∂∂βlog⁡LX∣μ,β,λ=-Nβ+1β∑i=1Nxi-μβλ=0(10)∂∂λlog⁡LX∣μ,β,λ=1λ2Nλ+Nlog⁡λ+NΦ1λ-1λ∑i=1Nxi-μβλlog⁡xi-μβ=0(11)∂∂μlog⁡LX∣μ,β,λ=∑i=1xi≥μNxi-μβλ-∑i=1xi<μNxi-μβλ=0where Φx is the digamma function,(12)Φx=dlog⁡Γxdx=Γ′xΓxHowever, we cannot achieve the close-form solution of the SWDM parameters by solving the above equation set (9)-(11). In this work, we use the gradient-based optimization method, Newton-Raphson iterative optimization scheme, to estimate the SWDM parameter.Firstly, we rewrite the three parameters of SWDM in the column vector form, namely,θ=μ,λ,βT. Then we update the parameter estimation results step by step based on the following updating rule:(13)θn+1=θn-H-1∇log⁡L~X∣μ,β,λwhere θn denotes the parameter estimation result of the nth step, ∇log⁡L~X∣θ is the first-order derivatives of log⁡L~X∣θ, and log⁡L~X∣θ represents the negative of the LLF log⁡LX∣μ,λ,β,(14)∇log⁡L~X∣θ=1βλ∑i=1xi<μNμ-xiλ-∑i=1xi≥μNxi-μλ1λ∑i=1Nxi-μβλlog⁡xi-μβ-N-Nlog⁡λλ-NλΦ1λNβ-1β∑i=1Nxi-μβλand H is the Hessian matrix and H=∇∇log⁡L~X∣θ.Given initial parameter estimationθ0, we repeat the parameter estimation by (13) until the estimation results converge, then we achieve the best estimation result of the SWDM parameters θ^, where the derivative vector and the Hessian matrix are as follows:(15)H=λβλ∑i=1xi<μNμ-xiλ-1+∑i=1xi≥μNxi-μλ-11βλ∑i=1xi<μNμ-xiλlog⁡μ-xiβ-∑i=1xi≥μNxi-μλlog⁡xi-μβ-λβλ+1∑i=1xi<μNμ-xiλ-∑i=1xi≥μNxi-μλ1βλ∑i=1xi<μNμ-xiλlog⁡μ-xiβ-∑i=1xi≥μNxi-μλlog⁡xi-μβ1λ∑i=1Nxi-μβλlog2⁡xi-μβ-log⁡xi-μ/βλ+Nλ2log⁡λ-1+Φ1λ-1λψ1λ-λβλ+1∑i=1xi<μNμ-xiλ-∑i=1xi≥μNxi-μλ-λβλ+1∑i=1xi<μNμ-xiλ-∑i=1xi≥μNxi-μλ-1β∑i=1Nxi-μβλlog⁡xi-μβ-Nβ2+λ+1βλ+2∑i=1Nxi-μλwhere ψz is the trigamma function and ψz=d2Γz/dz2.The convergence rate of the iterative estimation process is related to the initial model parameterθ0. To fasten the estimation process, the initial model parameter θ0 is assigned as an empirical value which is close to the extremum point; namely, θ0=μ0,λ0,β0T, where(16)μ0=medianXλ0=221/N∑i=1Nxi-x-4/1/N∑i=1Nxi-x-2-31/N∑i=1Nxi-x-4/1/N∑i=1Nxi-x-2-3+1β0=∑i=1Nxi-μ0λ0Nλ0It reveals that computing the inverse of the Hessian matrix is a time-consuming task. Actually, the location parameterμ is nearly the mean value of the sampling points and we can find the relation of the shape parameter λ and scale parameter β by the partial derivative equation (9); namely,(17)β=∑i=1Nxi-μλN1/λHence, we subtract a mean value from each sampling points and get a new corrected sample setX~, which still obeys the SWDM with the same scale parameter and shape parameter but symmetric around the y-axis. So, if we make a rough estimation of the location parameter with the mean value of the samples and take into account the relation equation of the scale parameter and shape parameter, namely, substituting (17) into (10), we can achieve a univariate equation as follows:(18)1+log⁡λλ+1λΦ1λ-1λ∑i=1Nx~iλ∑i=1Nx~iλlog⁡Nx~i∑i=1Nx~iλ=0wherex~i represents the samples in X~. We can solve (18) by the Newton-Raphson scheme and then based on the relation formula (17), we can obtain the scale parameter conveniently.The fitting goodness of the SWD model is evaluated by the statisticsχ2 and the Kullback-Leibler divergence (KLD), which is as follows:(19)χ2=∑i=1nhxi-fxi∣θ^fxi∣θ^(20)KLD=∑i=1nfxi∣θ^log⁡fxi∣θ^hxiwhere fxi∣θ^ and hxi are the expected SWDM probability and empirical probability over xi, respectively. The lower value of χ2 and KLD statistics, the more precise statistical model we achieve.


## 2.1. Theoretical Analysis of the WD Behavior of the LSD of Randomly Arranged Fragments

During visual imaging or scene observation process, local homogenous regions or objects in the imaging scene will fragment the field of view (FOV) into separated parts, e.g., the inside-object regions and outside-object regions [42]. If we enhance the imaging resolution, these inside-object regions and outside-object regions will be fragmented again by the presence of more LHFs or more complicated objects and finer details, and so on, until in any direction, the limitation of the FOV is reached. When the FOV is sufficiently fragmented by LHFs, the imaging scene becomes stable and the addition of any new details will not significantly change its appearance and consequently it leads to a specific visual appearance of TP. On the contrary, when reducing the resolution of the visual sensor, the adjacent local homogeneous structures will recombine together and constitute relatively coarser fragments. Geusebroek [42, 43] has demonstrated that this imaging process that fragmenting the visual surface details into LHFs can be simulated or analogized by the sequential fragmentation (SF) process well known in ore milling process, which is used to describe the distribution of the ore-particle size or mass in the ball grinding mill [44]. Hence, the distribution of the local fragments or particles in images can be predicted by the SF theory (SFT).In terms of SFT [45], the statistical distribution of LHFs based on the fragmentation process obeys the power-law distribution, which can be expressed by the following formula:(1)Fx′→x=xβλ-1where x′→x represents a single fragmentation process from a coarse fragment structure x′ to a fine fragment structure x and parameter β in the ore grinding process represents the average mass (or size given the ore density) of the crushed ore particles, whereas in the texture image analysis it denotes the average illumination intensity of the LHF x; λ is a shape parameter, satisfying 1≥λ≥0.With the increasing of the resolution power of visual sensor, more local fragment details are introduced into the imaging scene, which will divide the image into smaller LHFs (regions or particles). Hence, the histogram distribution of the illumination contrast is the integrating over the power-law-distributed LHFs [43], given by(2)nx=c∫x∞nx′Fx′→xdx′Substitute (1) into (2), and let c=1/β; we can obtain nx by solving the following equation:(3)nx=xβλ-1∫x∞nx′dx′βBy solving (3), it can be found that nx obeys the typical WD law; namely,(4)nx=NTxβλ-1e-1/λx/βλwhere NT=1/∫-∞∞nxdx is a normalized constant to make sure that nx is a reasonable distribution function.Since the resolution power of the visual sensor is limited, the fragmentation process will inevitably terminate and exhibit a stable TP; the statistical distribution of the SS of LHFs corresponds to the accumulated distribution of the fragments whose local contrast larger thanx, which is demonstrated obeying a symmetrical Weibull distribution model (SWDM) given by the following formula:(5)N>x=∫x∞nx′dx′=Ce-1/λx/βλTo obtain the generalization characterization of TPs, the SWDM can be modified a little by introducing a location parameterμ, indicating the distribution location (center), namely, the probability density of the SS of TPs can be described by a triparameter SWDM [25, 35, 46],(6)fx∣μ,β,λ=Ce-1/λx-μ/βλwhere C=1/∫-∞∞e-1/λx-μ/βλdx=λ/2λ1/λβΓ1/λ is a normalized constant and Γx=∫0∞tx-1e-tdt is the Gamma function.


## 2.2. Relations to the HVP System

It has demonstrated that SWDM can completely characterize the spatial layout of the stochastic textures [43]. In the SWDM depicted in (6), μ is a location parameter, representing the global reflectance of texture image, and it can be derived from the shape-from-shading method [47]; β is a scale parameter, directly related to the illumination contrast of the LHFs; λ is the shape parameter, reflecting the size of LHFs in texture image.The visual appearance of TP varies with the statistical distribution of LHFs. The psychophysical evidence indicates that two TPs are often difficult to discriminate when they produce similar distribution of responses of a bank of orientation and spatial-frequency selective linear filters [9]. From regular texture to pure random TP, the distribution of the SS of TP gradually changes from a power-law distribution to an exponential distribution and then to a typical Gaussian distribution. SWDM is a generalized model for the effective expression of LSDs of various forms by adjusting the model parameters [25, 35]. For instance, when λ=1, a SWDM will turns into an exponential distribution with the mean of β, when λ≃3.6, it is basically close to the Gaussian distribution. In addition, some studies have demonstrated that the shape parameter λ is directly related to the fractal dimension Df of a texture image. The fractal dimension can be approximated by Df=-3λ [46].Hence, it reveals that SWDM generate a complete orthogonal basis for stochastic textures description and SWDM parameters are directly related to the HVP properties [43, 46]. The commonly used HVP properties for texture perception includecoarseness,regularity,contrast,roughness, anddirectionality [9]. SWDM parameters can be used to make a psychophysical explanation of these perceptual properties [43], namely,coarseness,regularity,contrast,roughness, anddirectionality.Coarseness is a fragment size-related perceptual property. Generally, the larger the size of basic fragments in texture pattern, the coarser it is felt of the HVP system. In SWDM,coarseness can be reflected by the shape parameter λ and the smaller λ means the coarser TP. However,Coarseness is an observation scale-related parameter [51]. As a psychophysical evidence, if we magnify a texture image but without resolution enhancement of visual sensor, the HVP system will adjust the observation perception and keep the scale invariance with no more extra details addition in the FOV. SWDM is demonstrated to be well fits this scale invariance. In other words, the shape parameter λ remains unchanged under this circumstance. Alternatively, if the texture image is magnified with an increase of the resolving power, the shape parameter λ becomes smaller for the perception ofcoarser TP.Regularity is a layout variation-related HVP characteristic. In general, a finer TP is felt to be more regular. As analyzed above,coarseness (orfineness) can be reflected by the shape parameter λ; henceregularitycan be also expressed by the shape parameter λ. In the extreme case, a few particles or even just one fragment exists in the FOV, the LSD of the SS of TP rejects the SWDM because the real distribution obeys either a power-law distribution or a pure regular fragment (just one or two fragments fully occupy the entire FOV) and the exhibited statistical distribution is often multimodal. This phenomenon can be reflected with the over-fitting of the maximum likelihood estimation (MLE) of the shape parameter, e.g., λ≫2 [43].Contrast is illumination variation-related HVP characteristics [51]. It results from the illumination intensity and the surface height variations of the observation objects [47]. The global illumination variation is reflected by the location parameter μ, which indicates the inhomogeneous illumination surface, whereas the surface height variations can be reflected by the scale parameter β, thewidth of the SWDM according to the SFT. Hence, the perceptualcontrast can be indicated by two SWDM parameters, μ and β [43].Roughness is originally a tactile property in surface perception. Though it seems like a three-dimensional perceptual property, human beings can perceive it based on a two-dimensional image.Roughness is a HVP property results from the height variations of a certain homogeneous fragment. The greater the height variations of the fragments in the FOV, therougher perception it is felt. As analyzed above, β is an indicator of the height variations of TP, and the shape parameter λ is the reflectance ofcoarseness of the fragments in texture image. Hence,roughness can be reflected by the shape parameter λ and scale parameter β, under a special illumination condition, indicated by the location parameter μ [43].Directionality is a global sense of the whole FOV and caused by the shapes of the basis texture elements as well as their placement rules. Though the SWDM parameters do not contain the direct shape information, it reveals that the anisotropy of grains can be described with the dominant direction of the shape parameter λ[43–45]. Anisotropy in texture shadows (or contrast) can be reflected by the dominant orientation on the scale parameter β. If we fully consider the overall SS information of TP, the SWDM parameters can effectively describe the HVP characteristicdirectionality with the corresponding dominant orientation information.


## 2.3. Parameter Estimation of SWDM

The commonly used model parameter estimation method is maximum likelihood estimation (MLE). Given a sample setX=x1,x2,…,xN samples from a SWDM, the corresponding logarithmic likelihood function (LLF) log⁡LX∣μ,β,λ indicates how well the statistical model fit the training samples X. We can achieve the model parameter by maximizing LLF; namely,(7)μ^,β^,λ^=arg⁡maxμ,β,λ⁡log⁡LX∣μ,β,λ=arg⁡maxμ,β,λ⁡log⁡∏i=1Nλ2λ1/λβΓ1/λexp⁡-1λxi-μβλAfter expanding and simplifying the LLF, we obtain(8)log⁡LX∣μ,β,λ=Nlog⁡λ-1λlog⁡2λ-log⁡β-log⁡Γ1λ-∑i=1N1λxi-μβλThe SWDM parameters can be estimated by computing the partial derivative oflog⁡LX∣μ,β,λ with respect to μ,β,λ and setting each of the derivative to be equal to zero, respectively. Hence, we can obtain the following equation set:(9)∂∂βlog⁡LX∣μ,β,λ=-Nβ+1β∑i=1Nxi-μβλ=0(10)∂∂λlog⁡LX∣μ,β,λ=1λ2Nλ+Nlog⁡λ+NΦ1λ-1λ∑i=1Nxi-μβλlog⁡xi-μβ=0(11)∂∂μlog⁡LX∣μ,β,λ=∑i=1xi≥μNxi-μβλ-∑i=1xi<μNxi-μβλ=0where Φx is the digamma function,(12)Φx=dlog⁡Γxdx=Γ′xΓxHowever, we cannot achieve the close-form solution of the SWDM parameters by solving the above equation set (9)-(11). In this work, we use the gradient-based optimization method, Newton-Raphson iterative optimization scheme, to estimate the SWDM parameter.Firstly, we rewrite the three parameters of SWDM in the column vector form, namely,θ=μ,λ,βT. Then we update the parameter estimation results step by step based on the following updating rule:(13)θn+1=θn-H-1∇log⁡L~X∣μ,β,λwhere θn denotes the parameter estimation result of the nth step, ∇log⁡L~X∣θ is the first-order derivatives of log⁡L~X∣θ, and log⁡L~X∣θ represents the negative of the LLF log⁡LX∣μ,λ,β,(14)∇log⁡L~X∣θ=1βλ∑i=1xi<μNμ-xiλ-∑i=1xi≥μNxi-μλ1λ∑i=1Nxi-μβλlog⁡xi-μβ-N-Nlog⁡λλ-NλΦ1λNβ-1β∑i=1Nxi-μβλand H is the Hessian matrix and H=∇∇log⁡L~X∣θ.Given initial parameter estimationθ0, we repeat the parameter estimation by (13) until the estimation results converge, then we achieve the best estimation result of the SWDM parameters θ^, where the derivative vector and the Hessian matrix are as follows:(15)H=λβλ∑i=1xi<μNμ-xiλ-1+∑i=1xi≥μNxi-μλ-11βλ∑i=1xi<μNμ-xiλlog⁡μ-xiβ-∑i=1xi≥μNxi-μλlog⁡xi-μβ-λβλ+1∑i=1xi<μNμ-xiλ-∑i=1xi≥μNxi-μλ1βλ∑i=1xi<μNμ-xiλlog⁡μ-xiβ-∑i=1xi≥μNxi-μλlog⁡xi-μβ1λ∑i=1Nxi-μβλlog2⁡xi-μβ-log⁡xi-μ/βλ+Nλ2log⁡λ-1+Φ1λ-1λψ1λ-λβλ+1∑i=1xi<μNμ-xiλ-∑i=1xi≥μNxi-μλ-λβλ+1∑i=1xi<μNμ-xiλ-∑i=1xi≥μNxi-μλ-1β∑i=1Nxi-μβλlog⁡xi-μβ-Nβ2+λ+1βλ+2∑i=1Nxi-μλwhere ψz is the trigamma function and ψz=d2Γz/dz2.The convergence rate of the iterative estimation process is related to the initial model parameterθ0. To fasten the estimation process, the initial model parameter θ0 is assigned as an empirical value which is close to the extremum point; namely, θ0=μ0,λ0,β0T, where(16)μ0=medianXλ0=221/N∑i=1Nxi-x-4/1/N∑i=1Nxi-x-2-31/N∑i=1Nxi-x-4/1/N∑i=1Nxi-x-2-3+1β0=∑i=1Nxi-μ0λ0Nλ0It reveals that computing the inverse of the Hessian matrix is a time-consuming task. Actually, the location parameterμ is nearly the mean value of the sampling points and we can find the relation of the shape parameter λ and scale parameter β by the partial derivative equation (9); namely,(17)β=∑i=1Nxi-μλN1/λHence, we subtract a mean value from each sampling points and get a new corrected sample setX~, which still obeys the SWDM with the same scale parameter and shape parameter but symmetric around the y-axis. So, if we make a rough estimation of the location parameter with the mean value of the samples and take into account the relation equation of the scale parameter and shape parameter, namely, substituting (17) into (10), we can achieve a univariate equation as follows:(18)1+log⁡λλ+1λΦ1λ-1λ∑i=1Nx~iλ∑i=1Nx~iλlog⁡Nx~i∑i=1Nx~iλ=0wherex~i represents the samples in X~. We can solve (18) by the Newton-Raphson scheme and then based on the relation formula (17), we can obtain the scale parameter conveniently.The fitting goodness of the SWD model is evaluated by the statisticsχ2 and the Kullback-Leibler divergence (KLD), which is as follows:(19)χ2=∑i=1nhxi-fxi∣θ^fxi∣θ^(20)KLD=∑i=1nfxi∣θ^log⁡fxi∣θ^hxiwhere fxi∣θ^ and hxi are the expected SWDM probability and empirical probability over xi, respectively. The lower value of χ2 and KLD statistics, the more precise statistical model we achieve.


## 3. CSS Feature Characterization

### 3.1. Local Differential-Based SS Characterization

Digital image is a collect of discrete sampling points from a continuous two-dimensional function. Hence, the local SS at any pixel pointx,y around a given original point 0,0 in a visual image I can be determined by the Taylor expansion formula [42]. The k-order Taylor approximation I^x,y can be computed as follows [52]:(21)I^x,y=I^0,0+∑n=1k1n!∑i=0nniIxiyn-ixiyn-iwhere ni=n!/i!n-i! is the combinatorial number and Ixiyn-i is the n-order image derivatives of function I, namely, Ixiyn-i=∂n/∂xi∂yn-iI0,0, which involves the important SS information of I.The above Taylor expansion formula indicates that the image pixel observation value can be approximated by the weighting summation of the differential structures, representing the SS of a texture image, under a certain observation scale. Since the observed image is smoothed by a Gaussian kernel, the differential structureIxiyn-i is equivalent to the convolution operation between image I and a Gaussian derivative filter (GDF) [52]; namely,(22)Ixmynx,y=∂n+m∂xn∂ymIx,y∗Gx,y,σ=Ix,y∗Gxmynx,y,σwhere the symbol ∗ denotes the 2D convolution operation, Gx,y,σ is the Gaussian kernel function (GKF) with the scale σ, and Gxmyn represents the (m+n) order GDF with m-order partial derivative and n-order partial derivative over x and y direction, respectively, m≥0, n≥0.


### 3.2. Steerable Isotropic GDF (SIGDF)

The isotropic GKF-based GDF,Gxmyn, can achieve the spatial structure characteristics of a texture image in its corresponding x or y direction depending on the derivative orders of the isotropic GKF. Since the GKF is circularly symmetric, we call the isotropic GKF-based GDF as isotropic GDF (IGDF). Intuitively, many TPs exhibit obvious directional behavior. Hence, the directional information is necessary to be incorporated into the IGDF so that the SS characteristics of TP in any direction can be achieved, which will benefit TPC.We adopt a kind of steerable IGDF [53] (SIGDF) for SS characterization, which can achieve multidirectional SS feature of TP efficiently. To facilitate description, we denote Gκ,σ as a κ-order IGDF with the scale parameter σ. Suppose fθ represents the directional operation of the function f at the orientation θ and Gθκ,σ represents the result of Gκ,σ after rotating angle θ. According to the research results of Freeman [53], θ-oriented GDF Gθκ,σ can be computed by the following formula:(23)Gθκ,σx,y=∑i=1MkiθGθiκ,σx,yAs can be seen from (23), the θ-oriented IGDF can be constructed through the weighting summation of some basic IGDF on the specific orientations. The computational schematic can be seen in Figure 2. The key issue in this expression is to determine the number of basic IGDFs and the corresponding weighting value for each of the basis IGDF.Figure 2
Computational schematic of the filtering response of SIGDF of a specific directionθ.We transform theGκ,σ from Cartesian coordinates into the polar coordinates to obtain the number of IGDF bases, which is supposed to be M. Suppose the corresponding polar coordinates expression formula of Gκ,σx,y is Gκ,σr,ϕ, where r=x2+y2, ϕ=anglex,y, which can be also expressed with the Fourier series; namely,(24)FGκ,σr,ϕ=∑n=0Nanreinϕwhere F denotes the Fourier decomposition. Comparing (24) with (23), it indicates that the number M is equivalent to the number of nonzero harmonic components in the Fourier series of Gκ,σr,ϕ, in other words, Mis equal to the sum of nonzero number in anr. After the determination of the M-value, it is easy to achieve the weighting coefficient kiθ, which is in fact the solution of the following equation by solving the Fourier transformation under the polar coordinates:(25)1eiθ⋯eiNθ=11⋯1eiθ1eiθ2⋯eiθM⋮⋮⋱⋮eiNθ1eiNθ2⋯eiNθMk1θk2θ⋮kjθNamely,kiθ can be obtained by solving (25) after choosing the appropriate GDF bases. Then we can obtain a steerable IGDF on any direction. For instance, the steerable first-order IGDF and second-order IGDF can be computed as follows:(26)Gθ1,σ=cos⁡θG0∘1,σ+sin⁡θG90∘1,σ(27)Gθ2,σ=cos2⁡θG0∘2,σ+sin⁡θcos⁡θG45∘2,σ+sin2⁡θG90∘2,σFigure 3 displays the oriented filters of first-order, second-order, and fourth-order IGDF, Gθ1,σ, Gθ2,σ, and Gθ4,σ where θ ∈ [0~180∘].Figure 3
SIGDF withσ2=4, θ=kπ/10, and k=0,1,2,…,9. (a) Gθ1,σ, (b) Gθ2,σ, and (c) Gθ4,σ.
(a) (b) (c)


### 3.3. Oriented Anisotropic GDF (OAGDF)

Many TPs have obvious geometric anisotropy visually, whose local linear structures, edge structures, and some specific spatial layout cannot be characterized effectively by the IGDF. Hence, it should adopt the anisotropic GDF (AGDF) to characterize various complicated anisotropic SS of TPs.The commonly used two-variable elongated Gaussian kernel function (eGKF) in the two-dimensional Euclidean space is given by [54, 55](28)EGσ,ρx=12πσ2exp⁡-12σ2xTρ200ρ-2xwhere x=x,yT is the 2D coordinate, σ is referred as the scale parameter, and ρ is the anisotropic factor, satisfying ρ≥1. It is analogous to the SIGDF, the directional information should be introduced; hence the oriented AGKFs are generated by rotating the eGKF, given by(29)EGσ,ρ,θx=12πσ2exp⁡-12σ2RθxTρ200ρ-2Rθx=12πσ2exp⁡-12σ2xTR-θρ200ρ-2Rθxwhere Rθ=cos⁡θsin⁡θ-sin⁡θcos⁡θ is the rotation matrix and θ is the orientation of AGKF. Geusebroek [56] presented an alternative expression of AGKF by setting different standard deviations in two orthonormal directions (e.g., in (28) and (29), they used σx and σy with respect to x and y, rather than a universal scale σ); we adopt the abovementioned definition of AGKF with same scale parameter σ inspired by Shui’s work [55], because they have derived an explicit noise suppression formula and it is much easier to compute.AGDF is the first-order partial derivative of eGKF with respect to the first variablex, EGσ,ρ′x=∂EGσ,ρx/∂x=-ρ2x/σ2EGσ,ρx, and the oriented AGDF (OAGDF) is defined as the rotating expression of the AGDF, which reflects the SS variation of TP along the direction θ,(30)∇EGσ,ρ,θx=∂EGσ,ρ,θx∂θ=∂EGσ,ρRθx∂x=-ρ2cos⁡θ,sin⁡θxσ2EGσ,ρ,θxIt is worth noting that the OAGDF defined in (30) is not steerable (because ρ>1), so that the computation of the OAGDF in an arbitrary direction θ cannot be obtained by the weighting summation of some fixed GKF bases of fixed orthonormal directions [54] in analogous to the generation of the SIGDF as denoted in Figure 3.To facilitate application, it reveals that the AGKF can be decomposed into two Gaussian line filters of nonorthogonal directions in literature [56]. The main thought came from the property of Fourier transform or the convolution theorem; namely, the Fourier transform of a convolution is the point-wise product of the Fourier transforms. Hence, an eGKF can be separated as the convolution of two nonorthogonal Gaussian line filters in accordance with the product components of the eGKF in the Fourier domain.(31)EGσ,ρ,θx=12πσuσvexp⁡-12xcos⁡θ+ysin⁡θ2σu2+-xsin⁡θ+ycos⁡θ2σv2where σu=ρ-2σ2 and σv=σ2ρ2. The Fourier transform of Gσ,ρ,θx is(32)FEGσ,ρ,θxωx=exp⁡-12ωu2σu2+ωv2σv2where ωx=ωx,ωyT, ωu=ωxcos⁡θ+ysin⁡θ, and ωv=-ωxsin⁡θ+ycos⁡θ.The quadratic form of the exponential coefficients inFEGσ,ρ,θxωx can be rewritten as the following form: (33)ωu2σu2+ωv2σv2=a11ωx2+2a12ωxωy+a22ωy2To keep the response remaining the same, we get the corresponding weighted coefficients in (33) as(34)a11=σu2cos2⁡θ+σv2sin2⁡θa12=σu2-σv2cos⁡θsin⁡θa22=σv2cos2⁡θ+σv2sin⁡θHence the Fourier transform ofEGσ,ρ,θx can be rewritten as (35)FEGσ,ρ,θxωx=exp⁡-12a11-a122a22ωx2+a22ωy+a12a22ωx2According to the property of the exponential function,FEGσ,ρ,θxωxis a product of two Gaussian function; hence we can get the convolution form of EGσ,ρ,θx by the inverse transform; namely,(36)EGσ,ρ,θx=Gx,σx∗Gx2,σφwhere Gx,σx is the standard Gaussian kernel function with zero mean and σxSTD, e.g., Gx,σx=1/2πσxexp⁡-1/2x2/σx2, and x2 is the 2-norm of vector x, σx=a11-a122/a22=σuσv/σv2cos2⁡θ+σu2sin2⁡θ, σφ=a22=σv2cos2⁡θ+σu2sin2⁡θ, tan⁡φ=σv2cos2⁡θ+σu2sin2⁡θ/σu2-σv2cos⁡θsin⁡θ, and r=x2+xtan⁡φ2.It indicates that the AGKF can be decomposed into a 1D Gaussian filter with the standard deviation (STD)σx and another 1D Gaussian filter with STD σφ at the orientation of φ. It is worth noting that the orientation φ is not orthogonal to the θ orientation. Hence the OAGDF can be also computed as(37)∇EGσ,ρ,θx=∂EGσ,ρ,θx∂θ=∂Gx,σx*Gx2,σφ∂θ=∂Gx,σx∂θ*Gx2,σφWe display a set of AGKFs with their corresponding OAGDFs, which exhibit obvious anisotropic characteristic in Figure4. It can be used to extract the imperceptible changes in texture images along the orientation π/2+θ. It reveals that OAGDF with large-scale and large anisotropic factor can effectively extract the grayscale variation characteristics along the direction θ at a high signal-to-noise ratio (SNR); namely, the OAGDF can effectively characterize the SS of texture images, especially the edges along the direction π/2+θ, which involves the essential layout property of the SS of TP. The directional local SS of TP I along the orientation π/2+θ is computed by the convolution operator.Figure 4
Anisotropic GKFs with their corresponding OAGDFsσ2=9, ρ2=4, θ=kπ/10, and k=0,1,2,…,9. (a) AGKFs and (b) OAGDFs.
(a) (b)


## 3.1. Local Differential-Based SS Characterization

Digital image is a collect of discrete sampling points from a continuous two-dimensional function. Hence, the local SS at any pixel pointx,y around a given original point 0,0 in a visual image I can be determined by the Taylor expansion formula [42]. The k-order Taylor approximation I^x,y can be computed as follows [52]:(21)I^x,y=I^0,0+∑n=1k1n!∑i=0nniIxiyn-ixiyn-iwhere ni=n!/i!n-i! is the combinatorial number and Ixiyn-i is the n-order image derivatives of function I, namely, Ixiyn-i=∂n/∂xi∂yn-iI0,0, which involves the important SS information of I.The above Taylor expansion formula indicates that the image pixel observation value can be approximated by the weighting summation of the differential structures, representing the SS of a texture image, under a certain observation scale. Since the observed image is smoothed by a Gaussian kernel, the differential structureIxiyn-i is equivalent to the convolution operation between image I and a Gaussian derivative filter (GDF) [52]; namely,(22)Ixmynx,y=∂n+m∂xn∂ymIx,y∗Gx,y,σ=Ix,y∗Gxmynx,y,σwhere the symbol ∗ denotes the 2D convolution operation, Gx,y,σ is the Gaussian kernel function (GKF) with the scale σ, and Gxmyn represents the (m+n) order GDF with m-order partial derivative and n-order partial derivative over x and y direction, respectively, m≥0, n≥0.


## 3.2. Steerable Isotropic GDF (SIGDF)

The isotropic GKF-based GDF,Gxmyn, can achieve the spatial structure characteristics of a texture image in its corresponding x or y direction depending on the derivative orders of the isotropic GKF. Since the GKF is circularly symmetric, we call the isotropic GKF-based GDF as isotropic GDF (IGDF). Intuitively, many TPs exhibit obvious directional behavior. Hence, the directional information is necessary to be incorporated into the IGDF so that the SS characteristics of TP in any direction can be achieved, which will benefit TPC.We adopt a kind of steerable IGDF [53] (SIGDF) for SS characterization, which can achieve multidirectional SS feature of TP efficiently. To facilitate description, we denote Gκ,σ as a κ-order IGDF with the scale parameter σ. Suppose fθ represents the directional operation of the function f at the orientation θ and Gθκ,σ represents the result of Gκ,σ after rotating angle θ. According to the research results of Freeman [53], θ-oriented GDF Gθκ,σ can be computed by the following formula:(23)Gθκ,σx,y=∑i=1MkiθGθiκ,σx,yAs can be seen from (23), the θ-oriented IGDF can be constructed through the weighting summation of some basic IGDF on the specific orientations. The computational schematic can be seen in Figure 2. The key issue in this expression is to determine the number of basic IGDFs and the corresponding weighting value for each of the basis IGDF.Figure 2
Computational schematic of the filtering response of SIGDF of a specific directionθ.We transform theGκ,σ from Cartesian coordinates into the polar coordinates to obtain the number of IGDF bases, which is supposed to be M. Suppose the corresponding polar coordinates expression formula of Gκ,σx,y is Gκ,σr,ϕ, where r=x2+y2, ϕ=anglex,y, which can be also expressed with the Fourier series; namely,(24)FGκ,σr,ϕ=∑n=0Nanreinϕwhere F denotes the Fourier decomposition. Comparing (24) with (23), it indicates that the number M is equivalent to the number of nonzero harmonic components in the Fourier series of Gκ,σr,ϕ, in other words, Mis equal to the sum of nonzero number in anr. After the determination of the M-value, it is easy to achieve the weighting coefficient kiθ, which is in fact the solution of the following equation by solving the Fourier transformation under the polar coordinates:(25)1eiθ⋯eiNθ=11⋯1eiθ1eiθ2⋯eiθM⋮⋮⋱⋮eiNθ1eiNθ2⋯eiNθMk1θk2θ⋮kjθNamely,kiθ can be obtained by solving (25) after choosing the appropriate GDF bases. Then we can obtain a steerable IGDF on any direction. For instance, the steerable first-order IGDF and second-order IGDF can be computed as follows:(26)Gθ1,σ=cos⁡θG0∘1,σ+sin⁡θG90∘1,σ(27)Gθ2,σ=cos2⁡θG0∘2,σ+sin⁡θcos⁡θG45∘2,σ+sin2⁡θG90∘2,σFigure 3 displays the oriented filters of first-order, second-order, and fourth-order IGDF, Gθ1,σ, Gθ2,σ, and Gθ4,σ where θ ∈ [0~180∘].Figure 3
SIGDF withσ2=4, θ=kπ/10, and k=0,1,2,…,9. (a) Gθ1,σ, (b) Gθ2,σ, and (c) Gθ4,σ.
(a) (b) (c)


## 3.3. Oriented Anisotropic GDF (OAGDF)

Many TPs have obvious geometric anisotropy visually, whose local linear structures, edge structures, and some specific spatial layout cannot be characterized effectively by the IGDF. Hence, it should adopt the anisotropic GDF (AGDF) to characterize various complicated anisotropic SS of TPs.The commonly used two-variable elongated Gaussian kernel function (eGKF) in the two-dimensional Euclidean space is given by [54, 55](28)EGσ,ρx=12πσ2exp⁡-12σ2xTρ200ρ-2xwhere x=x,yT is the 2D coordinate, σ is referred as the scale parameter, and ρ is the anisotropic factor, satisfying ρ≥1. It is analogous to the SIGDF, the directional information should be introduced; hence the oriented AGKFs are generated by rotating the eGKF, given by(29)EGσ,ρ,θx=12πσ2exp⁡-12σ2RθxTρ200ρ-2Rθx=12πσ2exp⁡-12σ2xTR-θρ200ρ-2Rθxwhere Rθ=cos⁡θsin⁡θ-sin⁡θcos⁡θ is the rotation matrix and θ is the orientation of AGKF. Geusebroek [56] presented an alternative expression of AGKF by setting different standard deviations in two orthonormal directions (e.g., in (28) and (29), they used σx and σy with respect to x and y, rather than a universal scale σ); we adopt the abovementioned definition of AGKF with same scale parameter σ inspired by Shui’s work [55], because they have derived an explicit noise suppression formula and it is much easier to compute.AGDF is the first-order partial derivative of eGKF with respect to the first variablex, EGσ,ρ′x=∂EGσ,ρx/∂x=-ρ2x/σ2EGσ,ρx, and the oriented AGDF (OAGDF) is defined as the rotating expression of the AGDF, which reflects the SS variation of TP along the direction θ,(30)∇EGσ,ρ,θx=∂EGσ,ρ,θx∂θ=∂EGσ,ρRθx∂x=-ρ2cos⁡θ,sin⁡θxσ2EGσ,ρ,θxIt is worth noting that the OAGDF defined in (30) is not steerable (because ρ>1), so that the computation of the OAGDF in an arbitrary direction θ cannot be obtained by the weighting summation of some fixed GKF bases of fixed orthonormal directions [54] in analogous to the generation of the SIGDF as denoted in Figure 3.To facilitate application, it reveals that the AGKF can be decomposed into two Gaussian line filters of nonorthogonal directions in literature [56]. The main thought came from the property of Fourier transform or the convolution theorem; namely, the Fourier transform of a convolution is the point-wise product of the Fourier transforms. Hence, an eGKF can be separated as the convolution of two nonorthogonal Gaussian line filters in accordance with the product components of the eGKF in the Fourier domain.(31)EGσ,ρ,θx=12πσuσvexp⁡-12xcos⁡θ+ysin⁡θ2σu2+-xsin⁡θ+ycos⁡θ2σv2where σu=ρ-2σ2 and σv=σ2ρ2. The Fourier transform of Gσ,ρ,θx is(32)FEGσ,ρ,θxωx=exp⁡-12ωu2σu2+ωv2σv2where ωx=ωx,ωyT, ωu=ωxcos⁡θ+ysin⁡θ, and ωv=-ωxsin⁡θ+ycos⁡θ.The quadratic form of the exponential coefficients inFEGσ,ρ,θxωx can be rewritten as the following form: (33)ωu2σu2+ωv2σv2=a11ωx2+2a12ωxωy+a22ωy2To keep the response remaining the same, we get the corresponding weighted coefficients in (33) as(34)a11=σu2cos2⁡θ+σv2sin2⁡θa12=σu2-σv2cos⁡θsin⁡θa22=σv2cos2⁡θ+σv2sin⁡θHence the Fourier transform ofEGσ,ρ,θx can be rewritten as (35)FEGσ,ρ,θxωx=exp⁡-12a11-a122a22ωx2+a22ωy+a12a22ωx2According to the property of the exponential function,FEGσ,ρ,θxωxis a product of two Gaussian function; hence we can get the convolution form of EGσ,ρ,θx by the inverse transform; namely,(36)EGσ,ρ,θx=Gx,σx∗Gx2,σφwhere Gx,σx is the standard Gaussian kernel function with zero mean and σxSTD, e.g., Gx,σx=1/2πσxexp⁡-1/2x2/σx2, and x2 is the 2-norm of vector x, σx=a11-a122/a22=σuσv/σv2cos2⁡θ+σu2sin2⁡θ, σφ=a22=σv2cos2⁡θ+σu2sin2⁡θ, tan⁡φ=σv2cos2⁡θ+σu2sin2⁡θ/σu2-σv2cos⁡θsin⁡θ, and r=x2+xtan⁡φ2.It indicates that the AGKF can be decomposed into a 1D Gaussian filter with the standard deviation (STD)σx and another 1D Gaussian filter with STD σφ at the orientation of φ. It is worth noting that the orientation φ is not orthogonal to the θ orientation. Hence the OAGDF can be also computed as(37)∇EGσ,ρ,θx=∂EGσ,ρ,θx∂θ=∂Gx,σx*Gx2,σφ∂θ=∂Gx,σx∂θ*Gx2,σφWe display a set of AGKFs with their corresponding OAGDFs, which exhibit obvious anisotropic characteristic in Figure4. It can be used to extract the imperceptible changes in texture images along the orientation π/2+θ. It reveals that OAGDF with large-scale and large anisotropic factor can effectively extract the grayscale variation characteristics along the direction θ at a high signal-to-noise ratio (SNR); namely, the OAGDF can effectively characterize the SS of texture images, especially the edges along the direction π/2+θ, which involves the essential layout property of the SS of TP. The directional local SS of TP I along the orientation π/2+θ is computed by the convolution operator.Figure 4
Anisotropic GKFs with their corresponding OAGDFsσ2=9, ρ2=4, θ=kπ/10, and k=0,1,2,…,9. (a) AGKFs and (b) OAGDFs.
(a) (b)


## 4. Texture Pattern Classification

The proposed ISM-TPI mainly includes the texture image feature extraction module, classifier learning module, and taking actions with the trained classifier to TPC. We attempt to extract the essential SS characteristics of texture images in texture feature extraction based on the local differential-based SS representation and SWDM-based feature extraction. Schematic ISM-TPI can be found in Figure5.Figure 5
Schematic ISM-TPI.


### 4.1. SWDM-Based TP Feature Extraction

SIGDF and OAGDF are used to obtain the omnidirectional SS characteristics of texture image at first. Then, the triparameter SWDM of the filter responses is approximated to obtain the image’s structure feature vector. Figure6 displays a texture image with its SWDM-based image statistical modeling result. As can be seen from Figure 6, the SWDM can well fit the LSD of the GDF-based local SS characteristics and is apparently superior to the Gaussian distribution approximation result.Figure 6
SWDM-based statistical modeling of image spatial structure characteristics. (a) The original texture image, D23 in the Brodatz texture database; (b) the corresponding SSs of D23 with IGDF and AGDF, respectively; (c) SWDM-based statistical modeling results and the fitting results are apparently superior to the normal distribution (ND) estimation results based on the model evaluation indices, statisticsχ2, and the Kullback-Leibler divergence.To achieve the omnidirectional SS characteristics, we discretize the omnidirectional interval0,2π into P orientations uniformly. Then the SWDM parameters of the P orientations are obtained to construct the texture image feature. Figure 7 displays the polar plot of the omnidirectional SWDM parameter features of three texture images convoluted with two kinds of OGDFs: one is a second-order IGDF and the other is the OAGDF. The SWDM parameter-based features are obvious different between different texture patterns, which makes the automatic texture pattern classification possible.Figure 7
SWDM-based omnidirectional texture image feature.Besides the directional information, the observation scale is another important factor affecting the SS characterization. Hence, spatial structure characteristics observed in different scales are also analyzed by changing the size ofσ in this paper.(38)FITP=FMOSIGDFs,FMOOGDFswhere FMOIGDFs collects the SWDM parameters of the SS convolved by the SIGDFs in fMOSIGDF and FMOOGDFs collects the SWDM parameters of the texture image’s SS characteristics convolved with OAGDFs in fMOSIGDF.(39)FMOIGDFs=μσi,θp,λσi,θp,βσi,θp∣1≤i≤S,1≤p≤P(40)FMOOGDFs=μσi,ρj,θp,λσi,ρj,θp,βσi,ρj,θp∣1≤i≤S,1≤j≤H,1≤p≤Pwhere S is the total number of GKF scales, H is the number of anisotropic factor for MOOAGDFs generation, and μσi,θp, λσi,θpβσi,θp are the three parameters of the estimated SWDM of texture image I filtered with a SIGDF, whose Gaussian kernel scale is σi and the orientation is θp; μσi,ρj,θp, λσi,ρj,θp, βσi,ρj,θp are the corresponding three SWDM parameters of ITP I filtered with an OAGDF, whose orientation is θp, anisotropic factor is ρi, and the scale parameter of the eGKF is σi.Hence, the statistical modeling-based SS characteristics of TP consist of two parts: one part is the multiscale and omnidirectional SIGDFs- (MOSIGDF-) based SWDM parameter feature, denoting asFMOSIGDFs, and the other is the multiscale and omnidirectional OAGDFs- (MOOAGDFs-) based SWDM parameter features, abbreviated as FMOOGDFs. To obtain the local SS characteristics of TP, we can divide the image into multiple subblocks of fixed size, then each subblock can be treated as an independent image. Given B subblocks are conducted detailed texture image feature analysis, the expanded TP feature can be consequently expressed by(41)FITP_Expand=FITP∪FMOIGDFs,FMOOGDFsb∣b=1,2,…,Bwhere FMOIGDFs,FMOOGDFsb denotes the SWDM parameter feature of the local subimage b.


### 4.2. PLS-DA Based Classifier

TPC is a nonlinear function approximation task, with the multivariate image feature variables as the independent variables and the textural pattern labels as the dependent variables. Currently, the commonly used methods [57] to solve this task include the minimum distance classification (MDC), the K-nearest neighbor method (KNN), the soft independent modeling of class analogy classification (SIMCA), artificial neural network (ANN), the support vector machine (SVM), the partial least square-discriminant analysis (PLS-DA) method, and so on.The PLS-DA [58] is a discriminant analysis method based on the PLS regression, which considers the class member information provided by the auxiliary matrix in the process of constructing the classification model. PLS-DA has the advantages of PLS method in feature extraction and noise reduction. It has been revealed that PLS-DA can effectively achieve the dimension reduction of input variables so as to obtain good classification performance [59]. It inevitably gives information redundant in the SWDM-based texture feature because of the introduction of the filter bank. So, the PLS-DA based classifier is adopted to build the TP classification models in this work.Suppose there areT labeled texture image samples xt,ytt=1,2,…,T, where xt is the image SS feature vector of each texture images, whose corresponding category label is yt, the effective mapping model from xt to yt can be established by the PLS-DA method to achieve TPC. Based on the PLS-DA method, xt and yt can be expressed by using the following model:(42)xt=ztPT+etyt=mtQT+ηtwhere P and Q are parameter matrixes, satisfying PTP=I, QTQ=I, et and ηt are the regression residual error, and zt and mt are the latent variables in the PLS model, which are mainly used to build the correlation model between xt and yt. The linear regression model between zt and mt can be expressed as mt=ztB, where the regression factor B can be expressed as B=ztTzt-1ztTmt.As a result, in the classification task, firstly, the score vectorztest of the sample xtest to be tested is obtained according to parameter matrix P, then the category label y^test corresponding to the sample can be acquired based on the following expression(43)y^test=ztestBQFor the two-category classification problem, supposeω0,ω1 is the class label, the obtained category label y^test for any test image is not necessarily limited to {-1,1}; the following rules can be chosen to distinguish the image classes:(44)y^testclass=1"belong  toω0"ify^test≥y--1"belong  toω1"ify^test≥y-where y- is the appropriate threshold value and the selection of the suitable y- can minimize the classification error.Concerning the multiclass TPC task, the decision-making tree method can be used to solve this problem withn-classes samples; n-1 PLS-DA classifiers are constructed in series to realize the accurate classification and recognition of n-classes texture image samples.


## 4.1. SWDM-Based TP Feature Extraction

SIGDF and OAGDF are used to obtain the omnidirectional SS characteristics of texture image at first. Then, the triparameter SWDM of the filter responses is approximated to obtain the image’s structure feature vector. Figure6 displays a texture image with its SWDM-based image statistical modeling result. As can be seen from Figure 6, the SWDM can well fit the LSD of the GDF-based local SS characteristics and is apparently superior to the Gaussian distribution approximation result.Figure 6
SWDM-based statistical modeling of image spatial structure characteristics. (a) The original texture image, D23 in the Brodatz texture database; (b) the corresponding SSs of D23 with IGDF and AGDF, respectively; (c) SWDM-based statistical modeling results and the fitting results are apparently superior to the normal distribution (ND) estimation results based on the model evaluation indices, statisticsχ2, and the Kullback-Leibler divergence.To achieve the omnidirectional SS characteristics, we discretize the omnidirectional interval0,2π into P orientations uniformly. Then the SWDM parameters of the P orientations are obtained to construct the texture image feature. Figure 7 displays the polar plot of the omnidirectional SWDM parameter features of three texture images convoluted with two kinds of OGDFs: one is a second-order IGDF and the other is the OAGDF. The SWDM parameter-based features are obvious different between different texture patterns, which makes the automatic texture pattern classification possible.Figure 7
SWDM-based omnidirectional texture image feature.Besides the directional information, the observation scale is another important factor affecting the SS characterization. Hence, spatial structure characteristics observed in different scales are also analyzed by changing the size ofσ in this paper.(38)FITP=FMOSIGDFs,FMOOGDFswhere FMOIGDFs collects the SWDM parameters of the SS convolved by the SIGDFs in fMOSIGDF and FMOOGDFs collects the SWDM parameters of the texture image’s SS characteristics convolved with OAGDFs in fMOSIGDF.(39)FMOIGDFs=μσi,θp,λσi,θp,βσi,θp∣1≤i≤S,1≤p≤P(40)FMOOGDFs=μσi,ρj,θp,λσi,ρj,θp,βσi,ρj,θp∣1≤i≤S,1≤j≤H,1≤p≤Pwhere S is the total number of GKF scales, H is the number of anisotropic factor for MOOAGDFs generation, and μσi,θp, λσi,θpβσi,θp are the three parameters of the estimated SWDM of texture image I filtered with a SIGDF, whose Gaussian kernel scale is σi and the orientation is θp; μσi,ρj,θp, λσi,ρj,θp, βσi,ρj,θp are the corresponding three SWDM parameters of ITP I filtered with an OAGDF, whose orientation is θp, anisotropic factor is ρi, and the scale parameter of the eGKF is σi.Hence, the statistical modeling-based SS characteristics of TP consist of two parts: one part is the multiscale and omnidirectional SIGDFs- (MOSIGDF-) based SWDM parameter feature, denoting asFMOSIGDFs, and the other is the multiscale and omnidirectional OAGDFs- (MOOAGDFs-) based SWDM parameter features, abbreviated as FMOOGDFs. To obtain the local SS characteristics of TP, we can divide the image into multiple subblocks of fixed size, then each subblock can be treated as an independent image. Given B subblocks are conducted detailed texture image feature analysis, the expanded TP feature can be consequently expressed by(41)FITP_Expand=FITP∪FMOIGDFs,FMOOGDFsb∣b=1,2,…,Bwhere FMOIGDFs,FMOOGDFsb denotes the SWDM parameter feature of the local subimage b.


## 4.2. PLS-DA Based Classifier

TPC is a nonlinear function approximation task, with the multivariate image feature variables as the independent variables and the textural pattern labels as the dependent variables. Currently, the commonly used methods [57] to solve this task include the minimum distance classification (MDC), the K-nearest neighbor method (KNN), the soft independent modeling of class analogy classification (SIMCA), artificial neural network (ANN), the support vector machine (SVM), the partial least square-discriminant analysis (PLS-DA) method, and so on.The PLS-DA [58] is a discriminant analysis method based on the PLS regression, which considers the class member information provided by the auxiliary matrix in the process of constructing the classification model. PLS-DA has the advantages of PLS method in feature extraction and noise reduction. It has been revealed that PLS-DA can effectively achieve the dimension reduction of input variables so as to obtain good classification performance [59]. It inevitably gives information redundant in the SWDM-based texture feature because of the introduction of the filter bank. So, the PLS-DA based classifier is adopted to build the TP classification models in this work.Suppose there areT labeled texture image samples xt,ytt=1,2,…,T, where xt is the image SS feature vector of each texture images, whose corresponding category label is yt, the effective mapping model from xt to yt can be established by the PLS-DA method to achieve TPC. Based on the PLS-DA method, xt and yt can be expressed by using the following model:(42)xt=ztPT+etyt=mtQT+ηtwhere P and Q are parameter matrixes, satisfying PTP=I, QTQ=I, et and ηt are the regression residual error, and zt and mt are the latent variables in the PLS model, which are mainly used to build the correlation model between xt and yt. The linear regression model between zt and mt can be expressed as mt=ztB, where the regression factor B can be expressed as B=ztTzt-1ztTmt.As a result, in the classification task, firstly, the score vectorztest of the sample xtest to be tested is obtained according to parameter matrix P, then the category label y^test corresponding to the sample can be acquired based on the following expression(43)y^test=ztestBQFor the two-category classification problem, supposeω0,ω1 is the class label, the obtained category label y^test for any test image is not necessarily limited to {-1,1}; the following rules can be chosen to distinguish the image classes:(44)y^testclass=1"belong  toω0"ify^test≥y--1"belong  toω1"ify^test≥y-where y- is the appropriate threshold value and the selection of the suitable y- can minimize the classification error.Concerning the multiclass TPC task, the decision-making tree method can be used to solve this problem withn-classes samples; n-1 PLS-DA classifiers are constructed in series to realize the accurate classification and recognition of n-classes texture image samples.


## 5. Experimental Verification

### 5.1. Texture Data Sets and Experiment Setup

Three commonly used texture data sets, the Brodatz texture database (Brodatz [60]), the Columbia-Utrecht-Reflectance and Texture Database (CURet [61]), and the texture database (TDB [62]), were used to validate the classification accuracy (CA) of the proposed ISM-TPI and compare it with the state of the art. Detailed training and testing data sets are explained as follows.Thesmall-scale Brodatz data setBSs(24 classes)was chosen and processed as reported in literature [3]. The 24 homogeneous texture patterns downloaded from the normalized Brodatz’s texture (NBT) database [63] can be seen in Figure 8., each of which was originally 640 × 640 pixels and was partitioned into 25 nonoverlapped subimages of size of 128 × 128 pixels, of which 13 samples were used for classifier training and the remaining 12 samples of each TP were used for testing.Figure 8
24-class small-scale Brodatz data setBSc.Thelarge-scale challenging Brodatz data setBLs(112 classes)contains the entire 112 TPs in the NBT database. We call BLcas a challenging data set due to the great visual similarity on different texture patterns. For instance, some different texture patterns in the data set are visually belonging to the same vision scene with different imaging distance, e.g., D6 and D14, D23 and D27, D25 and D26, and D79 and D78, while some texture patterns may incur visual ambiguities, e.g., D40 and D42, D103 and D104, etc. Some easily confusing TPs can be seen in Figure 9. The training and testing samples were processed the same as the BSc data set.Figure 9
Some confusing texture patterns inBLc.Thesmall-scale CURet data setCSs(20 classes)was chosen from the CURet data set, where the texture patterns were captured from the real-world surface (with different capturing scales, different surface height variations and different orientations, which was believed can be used to span a wide range of geometric and photometric properties. We collect 90 images of the size of 128 × 128 (trimmed and rescaled) for each of the texture pattern and half the samples are used for classifier training. We display the texture pattern in Figure 10.Figure 10
20-class texture patterns ofCSc.Thefull CURet data setCF(61 classes)was chosen from the CURet database, where all the 61 real-world surface texture patterns were used for performance comparison. The samples were prepared the same with database CSc.The originalTDB data setTD (25 classes) includes 25 texture classes and each texture pattern contains 40 samples with the size of 640 × 480 pixels. The samples from each texture pattern are collect from different orientations and scales; some texture patterns can be seen from Figure 11. This data set is also a challenging testing set due to the large variability in the imaged appearance of the same texture pattern captured under different illuminations, positions, etc.Figure 11
TDB data setTD.Table1 summarizes the related information of the data sets used in the experiments. Since the multichannel color content does not introduce discriminative information (because they are very similar for all the textural images), we just consider the gray values of the texture images for feature analysis and TPC. Each experiment is repeated 30 times by randomly choosing the samples with replacement; as a result, the effectiveness of the method can be evaluated by calculating the average classification accuracies of the experiments with the standard deviations.Table 1
Summary of the data sets used in the experiments.
data set viewpoint change scale difference illumination variation texture classes sample size samples per class TrainingNT Test B S s   √ 24 128 × 128 25 13 12 B L s  √ √ 112 128 × 128 25 13 12 C S s √ √ √ 20 128 × 128 90 45 45 C F √ √ √ 61 128 × 128 90 45 45 T D √ √ √ 25 640 × 480 40 20 20


### 5.2. Textural Feature Characterization Methods in the Comparison Study

The proposed ISM-TPI method was compared with the state of the art. A summary of these comparative methods is displayed in Table2 and the implementations details are described as follows.Table 2
Summary of the comparative methods.
Methods Feature characteristic Feature dimension per image Classifier GLCM [48] Statistical measures 14 × 16 × N x / 32 × N y / 32 PLS-DA  GWT [49] Statistical measures 70 × N x / 32 × N y / 32 PLS-DA  MRMRF [50] +PLS-DA Statistical model parameters 46 × N x / 32 × N y / 32 PLS-DA  MRMRF+NLC Statistical model parameters 46 × N x / 32 × N y / 32 NLC approach [50]  VZ algorithm [11] Frequency histogram of learnedtextons 10 × N T P Nearest neighbor classifier(NNC) with theχ2 statistic for histogram distance measurement  PATCH [28] Frequency histogram of learnedtextons 10 × N T P NNC, Joint classifier  PATCH-MRF [28] Two-dimensional histogram 200 × 10 × N T P NNC, Joint clssifier  WLD [12] One-dimensional histogram encoded from a differential excitation and gradient orientation-paired two-dimensional histogram 8 × 6 × 20 K-nearest neighbor classifier, K=3, with normalized histogram intersection as the similarity measurement of two histogram  PRICoLBP [16] +SVM Two-dimensional histogram: Pairwise rotation invariant cooccurrence LBP feature 2 × 3 × 590 SVM  PRICoLBP+PLS-DA Two-dimensional histogram: Pairwise rotation invariant cooccurrence LBP feature 2 × 3 × 590 PLS-DA  ISM-TPI Statistical model parameters P × 3 × F × N x / 32 × N y / 32 PLS-DAGLCM [48]. GLCM method characterizes the texture image by calculating how often pairs of pixels with specific values and in a specified spatial relationship occur in an image. In the comparative study, texture images were firstly quantized into 16, 32, 64, and 256 grayscale levels. Then, we obtain 16 GLCMs with different spatial relationship settings on each specific gray level image. The 16 GLCMs can be denoted as {GLCMx,yli,θj∣,li∈2,5,10,15,θj∈0∘,45∘,90∘,135∘}, where li denotes the pixel offset and θj denote the orientation of two cooccurrence pixel intensitiesx andy. Successively, 14 Haralick statistical measures [48] were calculated based on each GLCMli,θj, including energy, contrast ratio, correlation, entropy, and so on. Consequently, a 14 × 16-tuple vector for each texture image was obtained by concatenating the average value of each statistical measure on four different grayscale levels. However, to make fuller representation of the texture structure, texture image can be divided into nonoverlapped subimages with the size of no smaller than 32 × 32 pixels and then each subimage can be treated as an independent image to obtain the corresponding GLCM-based statistical measures. Finally, we can obtain a 14×16×Nx/32×Ny/32-tuple feature vector for texture pattern analysis, where Nx and Ny represent the image width and image height, respectively. For a fair comparison; the classifier was PLS-DA, the classifier structure with the training and test samples were used the same as the proposed ISM-TPI.GWT [49]. We applied the Gabor wavelet transform with 5 scales and 7 directions on texture image and we obtain 35 Gabor filter bank responses. The Gabor filter bank are designed based on the strategy that ensure the half-peak magnitude support of the filter response in the frequency spectrum touch each other as described in literature [64]. The statistics, mean and variance, of the energy spectrum of each Gabor filtering responses were computed to constitute the GWT texture feature. It was treated the same as the GLCM feature; GWT features based on the local subblock were also obtained. Hence, the ultimate GWT-based texture feature is a 70×Nx/32×Ny/32-tuple feature vector.MRMRF [50]. It applies a 3-layer Haar wavelet decomposition to construct the multiresolution image representation firstly, and then the MRF model in each decomposition subband is established. Each subband image is modeled with 4 nonzero MRF model parameters. For the 9 detailed subbands, a total of 36 parameters can obtained. At the same time 10 wavelet energy signatures (including the low-pass subimage energy signature) are collected. Hence, we can obtain 36+10 parameters to generate the MRMRF feature vector. If we divided the image into nonoverlapped subimages of 32×32 pixels and each subimage is processed the same as an independent image. We can obtain a 46×Nx/32×Ny/32-tuple feature vector. The classifier in the comparative study we used the PLS-DA classifier as well as the NLC approach as displayed in [50].VZ algorithm [11]. Texture pattern is represented by the joint probability distribution of image filtering responses (MR8 filter bank is used in the experiments). The distribution is represented by the frequency histogram of filter response cluster centers (textons). For each texture class 10 texton cluster centers (texton) are learned by the standard K-means algorithm as experimented in [11]. Hence an image can be represented by a 10×NTP -bin histogram, and NTP is the number of texture patterns.PATCH [28]. Row pixel intensities of local patches of size 7 × 7 from the selected training images in a texture class are aggregated and clustered. The set of cluster centers from all the classes comprises thetexton dictionary. Both the training and the testing are performed in the patch domain. The classification is still achieved by nearest neighbor classifier with the χ2 statistic.PATCH-MRF [28]. A texture image is represented by a two-dimensional histogram: one dimension for the quantized bins of the patch center pixel and the other dimension for the learnedtextons from the patch with the center pixel excluded. As reported in [28], it chooses 200 bins for the central pixel with 10×NTP textons, which can achieve relative better performance.WLD [12]. A texture image is represented by a two-dimensional histogram, one dimension for the pixel’s differential excitation and the other for the orientation. Hence, a size of T×C histogram is achieved, where C represents the number of cells in each orientation. Then, each column of subhistogram is evenly divided into M(C=S×M) segments and then the two-dimensional histogram is regrouped and encoded into a one-dimensional histogram H. In the literature [12], the parameters are set as T=8, M=6, S=20.PRICoLBP [16]. PRICoLBP is an extended local binary pattern feature. According to the processing method reported in literature [16], we extract six cooccurrence patterns from three scales and two orientations and thus achieve a 2 × 3 × 590-dimensional feature vector. Concerning the classifier, the SVM classifier with RBF kernel andone-versus-the-rest strategy for the multiclass classification task was used in the literature [16]. In this work, the PLS-DA classifier is also applied for the classification performance comparison.


### 5.3. Experimental Results

#### 5.3.1. Validation Experiments

The SWDM-based TP feature extraction results are directly related to the SS characterization results, which depend on the configuration of GDFs,FMOSIGDFs and FMOOGDFs. To facilitate describing, we denote Gi as the i-th omnidirectional SIGDFs (OSIGDFs) with five specific scales, which includes five SIGDF templates, Gi,σj, where 1≤j≤5 and σj∈2,22,4,42,8, and we denote ∇EGρ as the omnidirectional OAGDF (OOAGDF) with the fixed anisotropic coefficient ρ under five specific scales σj∈2,22,4,42,8. The multidirectional SS characteristics of texture image in the real application are implemented by discretizing the range of 0,2π into P orientations uniformly. In the experiments, P=30. And regarding the local SS expression, we divide the original image into the nonoverlapped subimages with size of 32 × 32; namely, the samples in the BSc, BLc, CSc, and CF are partitioned into 16 subimages, while the samples in the TD are divided into 300 subimages, each subimage is treated as an independent image and the expanded TP feature are extracted as described in the formula (41).We firstly consider the classification performance based on the OSIGDFs and OOAGDFs independently. Figure12 displays the OSIGDFs-based classification accuracies (CAs) of the proposed ISM-TPI method on the five prepared data sets.Figure 12
Classification accuracies based on OSIGDFs.
(a) (b) (c) (d) (e)As can be seen from Figure12, in the small-scale data sets, e.g., BSs, CSs, it yields apparently higher classification accuracies, whereas on the challenging data sets, e.g., BLs, CF,the classification performance declines for the great similarities existing in the visually confusing TPs.Regarding the influence on the CAs of the differential-order of the oriented GDFs, the selection of the oriented third-order GDF alone,G3, achieves the relative lower classification accuracies on whatever the data set is, whereas the classification accuracies (CAs) based on the G1 and G2 are relatively higher than the exclusive usage of the G3. However, G1 and G2 generally reach close CAs on all of the five prepared data sets except on the data set CSc,where CAs of G1 are apparently higher than that of G2.The combination of more than one differential-orders OSIGDFs can achieve much higher CAs, especially the full combination of the three OSIGDFs, denoted asG1+G2+G3 in Figure 12, which can achieve the CAs as high as 92% on average on the data sets BSs, CSs, and TD and can also achieve CAs nearly close to the 90% on the challenging data sets BLs and CF.Some microstructural difference on different TP can be reflected by the local anisotropic and directional features. Hence, effective exhibition and representation of anisotropic features of the SS characteristics exhibiting in different texture patterns can facilitate the visually similar TP analysis in theory. We conduct experiments with only OAGDFs to obtain the anisotropic characteristics for classification performance comparison. Figure13 displays the average classification accuracies with standard deviation values on the five prepared five data sets depending on OOAGDFs with different anisotropic factors ρ, ρ∈2,3,5,7,9,11.Figure 13
TP classification accuracies based on OOAGDFs.Comparing the OSIGDFs-based CAs in Figure12 and the OOAGDFs-based CAs in Figure 13, OOAGDFs can achieve higher average classification accuracies (the average classification accuracy on the five data sets is almost as high as 94%) than that of OSIGDFs. The main reason is that the SSs of many TP exhibit obvious anisotropic behavior, which can be obtained by OOAGDFs with different parameter setting.However, OOAGDFs-based TPC would yield higher variance value of CAs on these data sets, because many TPs exhibit obvious anisotropic organization behavior of the local homogeneous fragments, when the anisotropic factor of the OSIGDF is consistent to the anisotropy of the SS of TP; it will achieve good classification performance and vice versa. Hence, it is necessary to combine OSIGDFs and OOAGDFs with multiple scales and multiple potential anisotropic factors to obtain sufficient characteristic details of the SS of TPs to achieve better classification performance on various TPs. Table3 enumerates the CAs on five data sets based on the combined OGDFs, including the first-order, second-order, and third-order IGDFs and AGDFs with five specific anisotropic factors, ρ∈2,3,5,7,9. In Table 3, the remarks, Com1=G1+G2+G3, Com2=∇EG2+∇EG3+∇EG5+∇EG7+∇EG9, and Com3=Com1+Com2, and the Gaussian scale parameter in OSIGDFs and OOAGDFs both include 5 different values, σj∈{2,22,4,42,8}, and the omnidirectionality is implemented by discretizing the range 0,2π into P=30 directions uniformly.Table 3
SWDM-based classification accuracies based on the filter bank combination of OSIGDFs and OOAGDFs.
B S s B L s C S s C F T D C o m 1 97.34 ± 4.34 92.34 ± 4.34 95.34 ± 4.34 94.34 ± 4.34 92.34 ± 3.55 C o m 2 98.67 ± 8.34 93.67 ± 3.34 97.67 ± 5.34 96.07 ± 3.34 92.97 ± 4.34 C o m 3 98.70 ± 0.89 94.70 ± 1.39 98.82 ± 0.79 98.04 ± 1.89 94.01 ± 0.89The classification results in Table3 indicate that the proposed ISM-TPI method can achieve quite high and stable classification accuracies on the five data sets. Especially, based on the full combination of OSIGDFs and OOAGDFs (G1+G2+G3+∇EG2+∇EG3+∇EG5+∇EG7+∇EG9), it achieves high classification accuracies over 98% on data sets, BSs, CSs, and CF and over 94% on the challenging data sets BLs and TD. The classification performance is quite stable considering the low STD values of CAs of the repeated experiments on almost all the five data sets, which are no greater than 2%.


#### 5.3.2. Comparative Experiments

In terms of the characteristic of texture features extracted by these comparative methods, they are all statistics-related texture analysis methods, which can be categorized into three classes. One is the statistical measurement-based methods, including, GLCM, GWT, one is the statistical distribution modeling-related methods, MRMRF, and the other is histogram-based local descriptor method, including WLD and PRICoLBP. A summary description of these comparative TPC methods, including the feature description, feature dimensions and the classifier selection, is displayed in Table2. The classification accuracies of the comparative methods on these data sets are displayed in Table 4, where the results marked with an asterisk are directly taken from the recent literature and marked with references.Table 4
Classification accuracies of comparative methods. Results marked with asterisk are directly taken from the literature with references and N/A means the corresponding index is absent in the referred literature.
Methods data sets B S s B L s C S s C F T D GLCM 89.21 ± 4.52 86.45 ± 5.34 91.45 ± 2.34 90.78 ± 5.34 84.45 ± 1.46 GWT 90.12 ± 3.45 89.12 ± 3.45 92.12 ± 1.45 91.32 ± 2.69 86.12 ± 4.89 MRMRF+PLS-DA 92.87 ± 4.36 89.87 ± 6.36 92.87 ± 3.36 89.87 ± 6.36 85.87 ± 4.36 MRMRF+NLC 91.45 ± 5.23 90.45 ± 1.23 91.85 ± 4.23 90.47 ± 3.82 84.45 ± 5.23 VZ algorithm 95.34 ± 4.54 91.45±2.01 97.12 ± 3.78 (96.40 ± N/A)∗ [28] 90.34 ± 2.14 PATCH 96.12 ± 3.45 (91.80 ± 0.98)∗ [3] 96.80 ± 0.79 (96.19 ± N/A)∗[28] 91.04 ± 2.02 PATCH-MRF 96.12 ± 3.43 92.04±1.68 97.12 ± 1.34 (97.47 ± N/A)∗[28] 92.89 ± 1.89 WLD (97.5 ± 0.6)∗ [12] 94.56 ± 0.83 98.92 ± 1.20 96.98 ± 2.12 93.45 ± 0.89 PRICoLBP+SVM (96.9 ± N/A)∗[16] 93.28±0.89 98.68 ± 0.98 98.42 ± 1.78 92.34 ± 1.23 PRICoLBP+PLS-DA 97.62 ± 2.29 93.89±1.02 98.01 ± 1.25 97.89 ± 0.89 93.12 ± 1.56 ISM-TPI (Com3) 98.70 ± 0.89 94.70 ± 1.39 98.82 ± 0.79 98.04 ± 1.89 94.01 ± 0.89As can be seen from Table4, the proposed ISM-TPI can achieve comparable or superior classification performance with the state-of-the-art methods. In terms of CA results on data sets, BSs, BLs, and TD, the proposed ISM-TPI achieves higher CAs than all the comparative methods. Regarding the CAs results on data sets CSs, almost all the methods can achieve higher classification performance; where the classic statistical measurement-based methods achieve about 92% CA, the methods VZ algorithm and PATCH can achieve as high as 97% CAs and the methods WLD and PRICoLBP can achieve the CAs higher over 98%. On data set CSs, the best CA result is 98.92%, coming from the WLD method, which obtains the histogram information of the differential excitation and the gradient orientation, representing the salient local texture pattern, or in other words, the spatial structural characteristic as the proposed ISM-TPI method. On data set CF, the best classification result is 98.42%, coming from the PRICoLBP+SVM-based TPC method. However, the average CAs of ISM-TPI on the two data sets are only inferior as low as 0.1% and 0.38% than the best classification results from WLD and PRICoLBP+SVM, respectively. At the same time, CAs of the proposed ISM-TPI are higher than the other state-of-the-art methods, e.g., VZ algorithm, Path, and Path-MRF on the two data sets.


#### 5.3.3. Discriminant Power of the SWDM-Based Feature

The main reason for achieving high classification accuracies is due to the stronger discriminant power of the SWDM-based texture features, which are demonstrated to be directly related to the perceptual properties of HVP system as explained in Section2.2. Hence, the proposed ISM-TPI method can obtain the essential spatial structural characteristics of texture images and can be adaptive to various texture patterns. In addition, the introduced SIGDFs and OAGDFs are effective filter banks for multiscale and omnidirectional spatial structural characteristics expression. Hence, SWDM-based features can achieve comparative and even superior classification performance compared with the state-of-the-art TPC methods.The SWDM parameters of two groups of typical confusing TPs are displayed in Figure14. D23, D27, and D31 represent the rubble-like texture pattern with different visual scales; the visual difference of these TPs can be expressed by the scale parameter β and location parameter μ, which is plot visually in the Figure 14(a). Regarding the TPs, D80, D82, and D84, whose surface organization of LHFs is resembling that may even mislead the visual identification of human beings. However, the essential distinctive SS of these TPs can be expressed clearly by the SWDM parameters λ,β,μ on different observation directions as can be seen in Figure 14(b).Figure 14
SWDM feature comparison of two groups of confusing TPs in Brodatz database, whose distinctive texture pattern characteristics can be expressed by the SWDM parametersλ,β,μ. (a) Rubble-like TPs, D23, D27, and D31. (b) Different TPs with great surface resemblance, D80, D82, and D84.
(a) (b)Some commonly used methods, e.g., GLCM, GWT, and MRMRF, cannot achieve high CAs on the challenging data sets mainly because that they may achieve similar feature parameters for the confusing TP analysis with great visual similarity. Figure15 displays four kind of texture feature extraction results of the confusing TPs, D80, D82, and D84, by the GLCM-based texture analysis method, GWT texture analysis method, WLD method, and VZ algorithm. The GLCM features are computed as follows. Image pixel gray-value of the original images are quantized into 8 gray levels and the GLCMs are computed under four orientation, 0°, 45°, 90°, and 135°, each orientation with 8 distances, 1,2,4,6,8,10,12,15; namely, we achieve 32 GLCMs for each ITP, and then the four GLCM-based statistics, contrast, correlation, energy, and homogeneity, are extracted and displayed in Figure 15(a).Figure 15
Feature extraction of TPs, D80, D82, and D84, based on GLCM, GWT, and VZ algorithm and WLD methods. (a) Original TPs, D80, D82, and D84, apparently with great visual similarity; (b) GLCM-based texture features; (c) GWT features, which tend to be the same trend in all the forty subbands; (d) histogram-based texture feature representation based on VZ algorithm proposed by Varma and Zisserman [11], with the first 100textons learned from the training image with MR8 filter banks class by class. PATCH and PATCH-MRF methods are their more recently texture classification methods, which has the similarity processing scheme but the original pixel patches are used and which can achieve the similar results for TP representation. (e) WLD histogram of D80, D82, and D84.
(a) (b) (c)  (d) (e)Regarding the GWT features, we compute the Gabor wavelet filters by the following construction formula:(45)hσ,θx,y;γ,fv,φ=π2σ2exp⁡-12x′2+γ2y′2σ2·cos⁡2πfvx′+φwhere x′=xcos⁡θ+ysin⁡θ, y′=ycos⁡θ-xsin⁡θ, σ is the standard deviation of the Gaussian factor (in pixels), θ represents the preferred orientation (in radians), φ denotes the phase offset (in radians), γ represents the spatial aspect ratio (of the x- and y-axis of the Gaussian ellipse), and 2π/fv means the preferred wavelength (period of the cosine factor, united in pixels).In the experiments, we choose 8 orientations in0.π, and five different wavelengths, namely, 2π/fv∈2,4,6,8,10,12, and we specify the γ as 2. Convolved with the forty Gabor filters, we can obtain the statistics mean and variance the texture images in each Gabor subband, which is plotted in Figure 15(c). The WLD features and the VZ algorithm features are histogram-based texture representation. The feature extraction steps fully follow the protocol in the literature [11, 12]. As can be seen clearly, GWT features, GLCM-based statistics, WLD histogram, and VZ algorithm-based texture features of D80, D82, and D84 exhibit strong similarity, some elements of which are even overlapped or only have very small difference, which seriously affected the discriminative ability of these methods for TPC.In a word, the state-of-the-art methods either compute some statistical measurement, e.g., GLCM, GWT, or make joint statistics of the image local structure with one-dimensional or two-dimensional histograms, whose resultant models are not based on the visual perception mechanism; hence, the ultimate results of the texture feature would have low discriminant power for the visual similarity texture patterns or make mistake in identifying the true TP of texture images sampled from large intravariant TPs, whereas the proposed SWDM reveals the essential LSD behavior of LHFs of the SS of TPs and the SWDM parameters under different observation scales and orientations can effectively obtain the distinctive visual features of TPs; hence the proposed ISM-TPI achieve relative higher classification performance, especially with better identification ability to the classification of the TPs with great visual resemblance or TPs with large visual variants.


## 5.1. Texture Data Sets and Experiment Setup

Three commonly used texture data sets, the Brodatz texture database (Brodatz [60]), the Columbia-Utrecht-Reflectance and Texture Database (CURet [61]), and the texture database (TDB [62]), were used to validate the classification accuracy (CA) of the proposed ISM-TPI and compare it with the state of the art. Detailed training and testing data sets are explained as follows.Thesmall-scale Brodatz data setBSs(24 classes)was chosen and processed as reported in literature [3]. The 24 homogeneous texture patterns downloaded from the normalized Brodatz’s texture (NBT) database [63] can be seen in Figure 8., each of which was originally 640 × 640 pixels and was partitioned into 25 nonoverlapped subimages of size of 128 × 128 pixels, of which 13 samples were used for classifier training and the remaining 12 samples of each TP were used for testing.Figure 8
24-class small-scale Brodatz data setBSc.Thelarge-scale challenging Brodatz data setBLs(112 classes)contains the entire 112 TPs in the NBT database. We call BLcas a challenging data set due to the great visual similarity on different texture patterns. For instance, some different texture patterns in the data set are visually belonging to the same vision scene with different imaging distance, e.g., D6 and D14, D23 and D27, D25 and D26, and D79 and D78, while some texture patterns may incur visual ambiguities, e.g., D40 and D42, D103 and D104, etc. Some easily confusing TPs can be seen in Figure 9. The training and testing samples were processed the same as the BSc data set.Figure 9
Some confusing texture patterns inBLc.Thesmall-scale CURet data setCSs(20 classes)was chosen from the CURet data set, where the texture patterns were captured from the real-world surface (with different capturing scales, different surface height variations and different orientations, which was believed can be used to span a wide range of geometric and photometric properties. We collect 90 images of the size of 128 × 128 (trimmed and rescaled) for each of the texture pattern and half the samples are used for classifier training. We display the texture pattern in Figure 10.Figure 10
20-class texture patterns ofCSc.Thefull CURet data setCF(61 classes)was chosen from the CURet database, where all the 61 real-world surface texture patterns were used for performance comparison. The samples were prepared the same with database CSc.The originalTDB data setTD (25 classes) includes 25 texture classes and each texture pattern contains 40 samples with the size of 640 × 480 pixels. The samples from each texture pattern are collect from different orientations and scales; some texture patterns can be seen from Figure 11. This data set is also a challenging testing set due to the large variability in the imaged appearance of the same texture pattern captured under different illuminations, positions, etc.Figure 11
TDB data setTD.Table1 summarizes the related information of the data sets used in the experiments. Since the multichannel color content does not introduce discriminative information (because they are very similar for all the textural images), we just consider the gray values of the texture images for feature analysis and TPC. Each experiment is repeated 30 times by randomly choosing the samples with replacement; as a result, the effectiveness of the method can be evaluated by calculating the average classification accuracies of the experiments with the standard deviations.Table 1
Summary of the data sets used in the experiments.
data set viewpoint change scale difference illumination variation texture classes sample size samples per class TrainingNT Test B S s   √ 24 128 × 128 25 13 12 B L s  √ √ 112 128 × 128 25 13 12 C S s √ √ √ 20 128 × 128 90 45 45 C F √ √ √ 61 128 × 128 90 45 45 T D √ √ √ 25 640 × 480 40 20 20


## 5.2. Textural Feature Characterization Methods in the Comparison Study

The proposed ISM-TPI method was compared with the state of the art. A summary of these comparative methods is displayed in Table2 and the implementations details are described as follows.Table 2
Summary of the comparative methods.
Methods Feature characteristic Feature dimension per image Classifier GLCM [48] Statistical measures 14 × 16 × N x / 32 × N y / 32 PLS-DA  GWT [49] Statistical measures 70 × N x / 32 × N y / 32 PLS-DA  MRMRF [50] +PLS-DA Statistical model parameters 46 × N x / 32 × N y / 32 PLS-DA  MRMRF+NLC Statistical model parameters 46 × N x / 32 × N y / 32 NLC approach [50]  VZ algorithm [11] Frequency histogram of learnedtextons 10 × N T P Nearest neighbor classifier(NNC) with theχ2 statistic for histogram distance measurement  PATCH [28] Frequency histogram of learnedtextons 10 × N T P NNC, Joint classifier  PATCH-MRF [28] Two-dimensional histogram 200 × 10 × N T P NNC, Joint clssifier  WLD [12] One-dimensional histogram encoded from a differential excitation and gradient orientation-paired two-dimensional histogram 8 × 6 × 20 K-nearest neighbor classifier, K=3, with normalized histogram intersection as the similarity measurement of two histogram  PRICoLBP [16] +SVM Two-dimensional histogram: Pairwise rotation invariant cooccurrence LBP feature 2 × 3 × 590 SVM  PRICoLBP+PLS-DA Two-dimensional histogram: Pairwise rotation invariant cooccurrence LBP feature 2 × 3 × 590 PLS-DA  ISM-TPI Statistical model parameters P × 3 × F × N x / 32 × N y / 32 PLS-DAGLCM [48]. GLCM method characterizes the texture image by calculating how often pairs of pixels with specific values and in a specified spatial relationship occur in an image. In the comparative study, texture images were firstly quantized into 16, 32, 64, and 256 grayscale levels. Then, we obtain 16 GLCMs with different spatial relationship settings on each specific gray level image. The 16 GLCMs can be denoted as {GLCMx,yli,θj∣,li∈2,5,10,15,θj∈0∘,45∘,90∘,135∘}, where li denotes the pixel offset and θj denote the orientation of two cooccurrence pixel intensitiesx andy. Successively, 14 Haralick statistical measures [48] were calculated based on each GLCMli,θj, including energy, contrast ratio, correlation, entropy, and so on. Consequently, a 14 × 16-tuple vector for each texture image was obtained by concatenating the average value of each statistical measure on four different grayscale levels. However, to make fuller representation of the texture structure, texture image can be divided into nonoverlapped subimages with the size of no smaller than 32 × 32 pixels and then each subimage can be treated as an independent image to obtain the corresponding GLCM-based statistical measures. Finally, we can obtain a 14×16×Nx/32×Ny/32-tuple feature vector for texture pattern analysis, where Nx and Ny represent the image width and image height, respectively. For a fair comparison; the classifier was PLS-DA, the classifier structure with the training and test samples were used the same as the proposed ISM-TPI.GWT [49]. We applied the Gabor wavelet transform with 5 scales and 7 directions on texture image and we obtain 35 Gabor filter bank responses. The Gabor filter bank are designed based on the strategy that ensure the half-peak magnitude support of the filter response in the frequency spectrum touch each other as described in literature [64]. The statistics, mean and variance, of the energy spectrum of each Gabor filtering responses were computed to constitute the GWT texture feature. It was treated the same as the GLCM feature; GWT features based on the local subblock were also obtained. Hence, the ultimate GWT-based texture feature is a 70×Nx/32×Ny/32-tuple feature vector.MRMRF [50]. It applies a 3-layer Haar wavelet decomposition to construct the multiresolution image representation firstly, and then the MRF model in each decomposition subband is established. Each subband image is modeled with 4 nonzero MRF model parameters. For the 9 detailed subbands, a total of 36 parameters can obtained. At the same time 10 wavelet energy signatures (including the low-pass subimage energy signature) are collected. Hence, we can obtain 36+10 parameters to generate the MRMRF feature vector. If we divided the image into nonoverlapped subimages of 32×32 pixels and each subimage is processed the same as an independent image. We can obtain a 46×Nx/32×Ny/32-tuple feature vector. The classifier in the comparative study we used the PLS-DA classifier as well as the NLC approach as displayed in [50].VZ algorithm [11]. Texture pattern is represented by the joint probability distribution of image filtering responses (MR8 filter bank is used in the experiments). The distribution is represented by the frequency histogram of filter response cluster centers (textons). For each texture class 10 texton cluster centers (texton) are learned by the standard K-means algorithm as experimented in [11]. Hence an image can be represented by a 10×NTP -bin histogram, and NTP is the number of texture patterns.PATCH [28]. Row pixel intensities of local patches of size 7 × 7 from the selected training images in a texture class are aggregated and clustered. The set of cluster centers from all the classes comprises thetexton dictionary. Both the training and the testing are performed in the patch domain. The classification is still achieved by nearest neighbor classifier with the χ2 statistic.PATCH-MRF [28]. A texture image is represented by a two-dimensional histogram: one dimension for the quantized bins of the patch center pixel and the other dimension for the learnedtextons from the patch with the center pixel excluded. As reported in [28], it chooses 200 bins for the central pixel with 10×NTP textons, which can achieve relative better performance.WLD [12]. A texture image is represented by a two-dimensional histogram, one dimension for the pixel’s differential excitation and the other for the orientation. Hence, a size of T×C histogram is achieved, where C represents the number of cells in each orientation. Then, each column of subhistogram is evenly divided into M(C=S×M) segments and then the two-dimensional histogram is regrouped and encoded into a one-dimensional histogram H. In the literature [12], the parameters are set as T=8, M=6, S=20.PRICoLBP [16]. PRICoLBP is an extended local binary pattern feature. According to the processing method reported in literature [16], we extract six cooccurrence patterns from three scales and two orientations and thus achieve a 2 × 3 × 590-dimensional feature vector. Concerning the classifier, the SVM classifier with RBF kernel andone-versus-the-rest strategy for the multiclass classification task was used in the literature [16]. In this work, the PLS-DA classifier is also applied for the classification performance comparison.


## 5.3. Experimental Results

### 5.3.1. Validation Experiments

The SWDM-based TP feature extraction results are directly related to the SS characterization results, which depend on the configuration of GDFs,FMOSIGDFs and FMOOGDFs. To facilitate describing, we denote Gi as the i-th omnidirectional SIGDFs (OSIGDFs) with five specific scales, which includes five SIGDF templates, Gi,σj, where 1≤j≤5 and σj∈2,22,4,42,8, and we denote ∇EGρ as the omnidirectional OAGDF (OOAGDF) with the fixed anisotropic coefficient ρ under five specific scales σj∈2,22,4,42,8. The multidirectional SS characteristics of texture image in the real application are implemented by discretizing the range of 0,2π into P orientations uniformly. In the experiments, P=30. And regarding the local SS expression, we divide the original image into the nonoverlapped subimages with size of 32 × 32; namely, the samples in the BSc, BLc, CSc, and CF are partitioned into 16 subimages, while the samples in the TD are divided into 300 subimages, each subimage is treated as an independent image and the expanded TP feature are extracted as described in the formula (41).We firstly consider the classification performance based on the OSIGDFs and OOAGDFs independently. Figure12 displays the OSIGDFs-based classification accuracies (CAs) of the proposed ISM-TPI method on the five prepared data sets.Figure 12
Classification accuracies based on OSIGDFs.
(a) (b) (c) (d) (e)As can be seen from Figure12, in the small-scale data sets, e.g., BSs, CSs, it yields apparently higher classification accuracies, whereas on the challenging data sets, e.g., BLs, CF,the classification performance declines for the great similarities existing in the visually confusing TPs.Regarding the influence on the CAs of the differential-order of the oriented GDFs, the selection of the oriented third-order GDF alone,G3, achieves the relative lower classification accuracies on whatever the data set is, whereas the classification accuracies (CAs) based on the G1 and G2 are relatively higher than the exclusive usage of the G3. However, G1 and G2 generally reach close CAs on all of the five prepared data sets except on the data set CSc,where CAs of G1 are apparently higher than that of G2.The combination of more than one differential-orders OSIGDFs can achieve much higher CAs, especially the full combination of the three OSIGDFs, denoted asG1+G2+G3 in Figure 12, which can achieve the CAs as high as 92% on average on the data sets BSs, CSs, and TD and can also achieve CAs nearly close to the 90% on the challenging data sets BLs and CF.Some microstructural difference on different TP can be reflected by the local anisotropic and directional features. Hence, effective exhibition and representation of anisotropic features of the SS characteristics exhibiting in different texture patterns can facilitate the visually similar TP analysis in theory. We conduct experiments with only OAGDFs to obtain the anisotropic characteristics for classification performance comparison. Figure13 displays the average classification accuracies with standard deviation values on the five prepared five data sets depending on OOAGDFs with different anisotropic factors ρ, ρ∈2,3,5,7,9,11.Figure 13
TP classification accuracies based on OOAGDFs.Comparing the OSIGDFs-based CAs in Figure12 and the OOAGDFs-based CAs in Figure 13, OOAGDFs can achieve higher average classification accuracies (the average classification accuracy on the five data sets is almost as high as 94%) than that of OSIGDFs. The main reason is that the SSs of many TP exhibit obvious anisotropic behavior, which can be obtained by OOAGDFs with different parameter setting.However, OOAGDFs-based TPC would yield higher variance value of CAs on these data sets, because many TPs exhibit obvious anisotropic organization behavior of the local homogeneous fragments, when the anisotropic factor of the OSIGDF is consistent to the anisotropy of the SS of TP; it will achieve good classification performance and vice versa. Hence, it is necessary to combine OSIGDFs and OOAGDFs with multiple scales and multiple potential anisotropic factors to obtain sufficient characteristic details of the SS of TPs to achieve better classification performance on various TPs. Table3 enumerates the CAs on five data sets based on the combined OGDFs, including the first-order, second-order, and third-order IGDFs and AGDFs with five specific anisotropic factors, ρ∈2,3,5,7,9. In Table 3, the remarks, Com1=G1+G2+G3, Com2=∇EG2+∇EG3+∇EG5+∇EG7+∇EG9, and Com3=Com1+Com2, and the Gaussian scale parameter in OSIGDFs and OOAGDFs both include 5 different values, σj∈{2,22,4,42,8}, and the omnidirectionality is implemented by discretizing the range 0,2π into P=30 directions uniformly.Table 3
SWDM-based classification accuracies based on the filter bank combination of OSIGDFs and OOAGDFs.
B S s B L s C S s C F T D C o m 1 97.34 ± 4.34 92.34 ± 4.34 95.34 ± 4.34 94.34 ± 4.34 92.34 ± 3.55 C o m 2 98.67 ± 8.34 93.67 ± 3.34 97.67 ± 5.34 96.07 ± 3.34 92.97 ± 4.34 C o m 3 98.70 ± 0.89 94.70 ± 1.39 98.82 ± 0.79 98.04 ± 1.89 94.01 ± 0.89The classification results in Table3 indicate that the proposed ISM-TPI method can achieve quite high and stable classification accuracies on the five data sets. Especially, based on the full combination of OSIGDFs and OOAGDFs (G1+G2+G3+∇EG2+∇EG3+∇EG5+∇EG7+∇EG9), it achieves high classification accuracies over 98% on data sets, BSs, CSs, and CF and over 94% on the challenging data sets BLs and TD. The classification performance is quite stable considering the low STD values of CAs of the repeated experiments on almost all the five data sets, which are no greater than 2%.


### 5.3.2. Comparative Experiments

In terms of the characteristic of texture features extracted by these comparative methods, they are all statistics-related texture analysis methods, which can be categorized into three classes. One is the statistical measurement-based methods, including, GLCM, GWT, one is the statistical distribution modeling-related methods, MRMRF, and the other is histogram-based local descriptor method, including WLD and PRICoLBP. A summary description of these comparative TPC methods, including the feature description, feature dimensions and the classifier selection, is displayed in Table2. The classification accuracies of the comparative methods on these data sets are displayed in Table 4, where the results marked with an asterisk are directly taken from the recent literature and marked with references.Table 4
Classification accuracies of comparative methods. Results marked with asterisk are directly taken from the literature with references and N/A means the corresponding index is absent in the referred literature.
Methods data sets B S s B L s C S s C F T D GLCM 89.21 ± 4.52 86.45 ± 5.34 91.45 ± 2.34 90.78 ± 5.34 84.45 ± 1.46 GWT 90.12 ± 3.45 89.12 ± 3.45 92.12 ± 1.45 91.32 ± 2.69 86.12 ± 4.89 MRMRF+PLS-DA 92.87 ± 4.36 89.87 ± 6.36 92.87 ± 3.36 89.87 ± 6.36 85.87 ± 4.36 MRMRF+NLC 91.45 ± 5.23 90.45 ± 1.23 91.85 ± 4.23 90.47 ± 3.82 84.45 ± 5.23 VZ algorithm 95.34 ± 4.54 91.45±2.01 97.12 ± 3.78 (96.40 ± N/A)∗ [28] 90.34 ± 2.14 PATCH 96.12 ± 3.45 (91.80 ± 0.98)∗ [3] 96.80 ± 0.79 (96.19 ± N/A)∗[28] 91.04 ± 2.02 PATCH-MRF 96.12 ± 3.43 92.04±1.68 97.12 ± 1.34 (97.47 ± N/A)∗[28] 92.89 ± 1.89 WLD (97.5 ± 0.6)∗ [12] 94.56 ± 0.83 98.92 ± 1.20 96.98 ± 2.12 93.45 ± 0.89 PRICoLBP+SVM (96.9 ± N/A)∗[16] 93.28±0.89 98.68 ± 0.98 98.42 ± 1.78 92.34 ± 1.23 PRICoLBP+PLS-DA 97.62 ± 2.29 93.89±1.02 98.01 ± 1.25 97.89 ± 0.89 93.12 ± 1.56 ISM-TPI (Com3) 98.70 ± 0.89 94.70 ± 1.39 98.82 ± 0.79 98.04 ± 1.89 94.01 ± 0.89As can be seen from Table4, the proposed ISM-TPI can achieve comparable or superior classification performance with the state-of-the-art methods. In terms of CA results on data sets, BSs, BLs, and TD, the proposed ISM-TPI achieves higher CAs than all the comparative methods. Regarding the CAs results on data sets CSs, almost all the methods can achieve higher classification performance; where the classic statistical measurement-based methods achieve about 92% CA, the methods VZ algorithm and PATCH can achieve as high as 97% CAs and the methods WLD and PRICoLBP can achieve the CAs higher over 98%. On data set CSs, the best CA result is 98.92%, coming from the WLD method, which obtains the histogram information of the differential excitation and the gradient orientation, representing the salient local texture pattern, or in other words, the spatial structural characteristic as the proposed ISM-TPI method. On data set CF, the best classification result is 98.42%, coming from the PRICoLBP+SVM-based TPC method. However, the average CAs of ISM-TPI on the two data sets are only inferior as low as 0.1% and 0.38% than the best classification results from WLD and PRICoLBP+SVM, respectively. At the same time, CAs of the proposed ISM-TPI are higher than the other state-of-the-art methods, e.g., VZ algorithm, Path, and Path-MRF on the two data sets.


### 5.3.3. Discriminant Power of the SWDM-Based Feature

The main reason for achieving high classification accuracies is due to the stronger discriminant power of the SWDM-based texture features, which are demonstrated to be directly related to the perceptual properties of HVP system as explained in Section2.2. Hence, the proposed ISM-TPI method can obtain the essential spatial structural characteristics of texture images and can be adaptive to various texture patterns. In addition, the introduced SIGDFs and OAGDFs are effective filter banks for multiscale and omnidirectional spatial structural characteristics expression. Hence, SWDM-based features can achieve comparative and even superior classification performance compared with the state-of-the-art TPC methods.The SWDM parameters of two groups of typical confusing TPs are displayed in Figure14. D23, D27, and D31 represent the rubble-like texture pattern with different visual scales; the visual difference of these TPs can be expressed by the scale parameter β and location parameter μ, which is plot visually in the Figure 14(a). Regarding the TPs, D80, D82, and D84, whose surface organization of LHFs is resembling that may even mislead the visual identification of human beings. However, the essential distinctive SS of these TPs can be expressed clearly by the SWDM parameters λ,β,μ on different observation directions as can be seen in Figure 14(b).Figure 14
SWDM feature comparison of two groups of confusing TPs in Brodatz database, whose distinctive texture pattern characteristics can be expressed by the SWDM parametersλ,β,μ. (a) Rubble-like TPs, D23, D27, and D31. (b) Different TPs with great surface resemblance, D80, D82, and D84.
(a) (b)Some commonly used methods, e.g., GLCM, GWT, and MRMRF, cannot achieve high CAs on the challenging data sets mainly because that they may achieve similar feature parameters for the confusing TP analysis with great visual similarity. Figure15 displays four kind of texture feature extraction results of the confusing TPs, D80, D82, and D84, by the GLCM-based texture analysis method, GWT texture analysis method, WLD method, and VZ algorithm. The GLCM features are computed as follows. Image pixel gray-value of the original images are quantized into 8 gray levels and the GLCMs are computed under four orientation, 0°, 45°, 90°, and 135°, each orientation with 8 distances, 1,2,4,6,8,10,12,15; namely, we achieve 32 GLCMs for each ITP, and then the four GLCM-based statistics, contrast, correlation, energy, and homogeneity, are extracted and displayed in Figure 15(a).Figure 15
Feature extraction of TPs, D80, D82, and D84, based on GLCM, GWT, and VZ algorithm and WLD methods. (a) Original TPs, D80, D82, and D84, apparently with great visual similarity; (b) GLCM-based texture features; (c) GWT features, which tend to be the same trend in all the forty subbands; (d) histogram-based texture feature representation based on VZ algorithm proposed by Varma and Zisserman [11], with the first 100textons learned from the training image with MR8 filter banks class by class. PATCH and PATCH-MRF methods are their more recently texture classification methods, which has the similarity processing scheme but the original pixel patches are used and which can achieve the similar results for TP representation. (e) WLD histogram of D80, D82, and D84.
(a) (b) (c)  (d) (e)Regarding the GWT features, we compute the Gabor wavelet filters by the following construction formula:(45)hσ,θx,y;γ,fv,φ=π2σ2exp⁡-12x′2+γ2y′2σ2·cos⁡2πfvx′+φwhere x′=xcos⁡θ+ysin⁡θ, y′=ycos⁡θ-xsin⁡θ, σ is the standard deviation of the Gaussian factor (in pixels), θ represents the preferred orientation (in radians), φ denotes the phase offset (in radians), γ represents the spatial aspect ratio (of the x- and y-axis of the Gaussian ellipse), and 2π/fv means the preferred wavelength (period of the cosine factor, united in pixels).In the experiments, we choose 8 orientations in0.π, and five different wavelengths, namely, 2π/fv∈2,4,6,8,10,12, and we specify the γ as 2. Convolved with the forty Gabor filters, we can obtain the statistics mean and variance the texture images in each Gabor subband, which is plotted in Figure 15(c). The WLD features and the VZ algorithm features are histogram-based texture representation. The feature extraction steps fully follow the protocol in the literature [11, 12]. As can be seen clearly, GWT features, GLCM-based statistics, WLD histogram, and VZ algorithm-based texture features of D80, D82, and D84 exhibit strong similarity, some elements of which are even overlapped or only have very small difference, which seriously affected the discriminative ability of these methods for TPC.In a word, the state-of-the-art methods either compute some statistical measurement, e.g., GLCM, GWT, or make joint statistics of the image local structure with one-dimensional or two-dimensional histograms, whose resultant models are not based on the visual perception mechanism; hence, the ultimate results of the texture feature would have low discriminant power for the visual similarity texture patterns or make mistake in identifying the true TP of texture images sampled from large intravariant TPs, whereas the proposed SWDM reveals the essential LSD behavior of LHFs of the SS of TPs and the SWDM parameters under different observation scales and orientations can effectively obtain the distinctive visual features of TPs; hence the proposed ISM-TPI achieve relative higher classification performance, especially with better identification ability to the classification of the TPs with great visual resemblance or TPs with large visual variants.


## 5.3.1. Validation Experiments

The SWDM-based TP feature extraction results are directly related to the SS characterization results, which depend on the configuration of GDFs,FMOSIGDFs and FMOOGDFs. To facilitate describing, we denote Gi as the i-th omnidirectional SIGDFs (OSIGDFs) with five specific scales, which includes five SIGDF templates, Gi,σj, where 1≤j≤5 and σj∈2,22,4,42,8, and we denote ∇EGρ as the omnidirectional OAGDF (OOAGDF) with the fixed anisotropic coefficient ρ under five specific scales σj∈2,22,4,42,8. The multidirectional SS characteristics of texture image in the real application are implemented by discretizing the range of 0,2π into P orientations uniformly. In the experiments, P=30. And regarding the local SS expression, we divide the original image into the nonoverlapped subimages with size of 32 × 32; namely, the samples in the BSc, BLc, CSc, and CF are partitioned into 16 subimages, while the samples in the TD are divided into 300 subimages, each subimage is treated as an independent image and the expanded TP feature are extracted as described in the formula (41).We firstly consider the classification performance based on the OSIGDFs and OOAGDFs independently. Figure12 displays the OSIGDFs-based classification accuracies (CAs) of the proposed ISM-TPI method on the five prepared data sets.Figure 12
Classification accuracies based on OSIGDFs.
(a) (b) (c) (d) (e)As can be seen from Figure12, in the small-scale data sets, e.g., BSs, CSs, it yields apparently higher classification accuracies, whereas on the challenging data sets, e.g., BLs, CF,the classification performance declines for the great similarities existing in the visually confusing TPs.Regarding the influence on the CAs of the differential-order of the oriented GDFs, the selection of the oriented third-order GDF alone,G3, achieves the relative lower classification accuracies on whatever the data set is, whereas the classification accuracies (CAs) based on the G1 and G2 are relatively higher than the exclusive usage of the G3. However, G1 and G2 generally reach close CAs on all of the five prepared data sets except on the data set CSc,where CAs of G1 are apparently higher than that of G2.The combination of more than one differential-orders OSIGDFs can achieve much higher CAs, especially the full combination of the three OSIGDFs, denoted asG1+G2+G3 in Figure 12, which can achieve the CAs as high as 92% on average on the data sets BSs, CSs, and TD and can also achieve CAs nearly close to the 90% on the challenging data sets BLs and CF.Some microstructural difference on different TP can be reflected by the local anisotropic and directional features. Hence, effective exhibition and representation of anisotropic features of the SS characteristics exhibiting in different texture patterns can facilitate the visually similar TP analysis in theory. We conduct experiments with only OAGDFs to obtain the anisotropic characteristics for classification performance comparison. Figure13 displays the average classification accuracies with standard deviation values on the five prepared five data sets depending on OOAGDFs with different anisotropic factors ρ, ρ∈2,3,5,7,9,11.Figure 13
TP classification accuracies based on OOAGDFs.Comparing the OSIGDFs-based CAs in Figure12 and the OOAGDFs-based CAs in Figure 13, OOAGDFs can achieve higher average classification accuracies (the average classification accuracy on the five data sets is almost as high as 94%) than that of OSIGDFs. The main reason is that the SSs of many TP exhibit obvious anisotropic behavior, which can be obtained by OOAGDFs with different parameter setting.However, OOAGDFs-based TPC would yield higher variance value of CAs on these data sets, because many TPs exhibit obvious anisotropic organization behavior of the local homogeneous fragments, when the anisotropic factor of the OSIGDF is consistent to the anisotropy of the SS of TP; it will achieve good classification performance and vice versa. Hence, it is necessary to combine OSIGDFs and OOAGDFs with multiple scales and multiple potential anisotropic factors to obtain sufficient characteristic details of the SS of TPs to achieve better classification performance on various TPs. Table3 enumerates the CAs on five data sets based on the combined OGDFs, including the first-order, second-order, and third-order IGDFs and AGDFs with five specific anisotropic factors, ρ∈2,3,5,7,9. In Table 3, the remarks, Com1=G1+G2+G3, Com2=∇EG2+∇EG3+∇EG5+∇EG7+∇EG9, and Com3=Com1+Com2, and the Gaussian scale parameter in OSIGDFs and OOAGDFs both include 5 different values, σj∈{2,22,4,42,8}, and the omnidirectionality is implemented by discretizing the range 0,2π into P=30 directions uniformly.Table 3
SWDM-based classification accuracies based on the filter bank combination of OSIGDFs and OOAGDFs.
B S s B L s C S s C F T D C o m 1 97.34 ± 4.34 92.34 ± 4.34 95.34 ± 4.34 94.34 ± 4.34 92.34 ± 3.55 C o m 2 98.67 ± 8.34 93.67 ± 3.34 97.67 ± 5.34 96.07 ± 3.34 92.97 ± 4.34 C o m 3 98.70 ± 0.89 94.70 ± 1.39 98.82 ± 0.79 98.04 ± 1.89 94.01 ± 0.89The classification results in Table3 indicate that the proposed ISM-TPI method can achieve quite high and stable classification accuracies on the five data sets. Especially, based on the full combination of OSIGDFs and OOAGDFs (G1+G2+G3+∇EG2+∇EG3+∇EG5+∇EG7+∇EG9), it achieves high classification accuracies over 98% on data sets, BSs, CSs, and CF and over 94% on the challenging data sets BLs and TD. The classification performance is quite stable considering the low STD values of CAs of the repeated experiments on almost all the five data sets, which are no greater than 2%.


## 5.3.2. Comparative Experiments

In terms of the characteristic of texture features extracted by these comparative methods, they are all statistics-related texture analysis methods, which can be categorized into three classes. One is the statistical measurement-based methods, including, GLCM, GWT, one is the statistical distribution modeling-related methods, MRMRF, and the other is histogram-based local descriptor method, including WLD and PRICoLBP. A summary description of these comparative TPC methods, including the feature description, feature dimensions and the classifier selection, is displayed in Table2. The classification accuracies of the comparative methods on these data sets are displayed in Table 4, where the results marked with an asterisk are directly taken from the recent literature and marked with references.Table 4
Classification accuracies of comparative methods. Results marked with asterisk are directly taken from the literature with references and N/A means the corresponding index is absent in the referred literature.
Methods data sets B S s B L s C S s C F T D GLCM 89.21 ± 4.52 86.45 ± 5.34 91.45 ± 2.34 90.78 ± 5.34 84.45 ± 1.46 GWT 90.12 ± 3.45 89.12 ± 3.45 92.12 ± 1.45 91.32 ± 2.69 86.12 ± 4.89 MRMRF+PLS-DA 92.87 ± 4.36 89.87 ± 6.36 92.87 ± 3.36 89.87 ± 6.36 85.87 ± 4.36 MRMRF+NLC 91.45 ± 5.23 90.45 ± 1.23 91.85 ± 4.23 90.47 ± 3.82 84.45 ± 5.23 VZ algorithm 95.34 ± 4.54 91.45±2.01 97.12 ± 3.78 (96.40 ± N/A)∗ [28] 90.34 ± 2.14 PATCH 96.12 ± 3.45 (91.80 ± 0.98)∗ [3] 96.80 ± 0.79 (96.19 ± N/A)∗[28] 91.04 ± 2.02 PATCH-MRF 96.12 ± 3.43 92.04±1.68 97.12 ± 1.34 (97.47 ± N/A)∗[28] 92.89 ± 1.89 WLD (97.5 ± 0.6)∗ [12] 94.56 ± 0.83 98.92 ± 1.20 96.98 ± 2.12 93.45 ± 0.89 PRICoLBP+SVM (96.9 ± N/A)∗[16] 93.28±0.89 98.68 ± 0.98 98.42 ± 1.78 92.34 ± 1.23 PRICoLBP+PLS-DA 97.62 ± 2.29 93.89±1.02 98.01 ± 1.25 97.89 ± 0.89 93.12 ± 1.56 ISM-TPI (Com3) 98.70 ± 0.89 94.70 ± 1.39 98.82 ± 0.79 98.04 ± 1.89 94.01 ± 0.89As can be seen from Table4, the proposed ISM-TPI can achieve comparable or superior classification performance with the state-of-the-art methods. In terms of CA results on data sets, BSs, BLs, and TD, the proposed ISM-TPI achieves higher CAs than all the comparative methods. Regarding the CAs results on data sets CSs, almost all the methods can achieve higher classification performance; where the classic statistical measurement-based methods achieve about 92% CA, the methods VZ algorithm and PATCH can achieve as high as 97% CAs and the methods WLD and PRICoLBP can achieve the CAs higher over 98%. On data set CSs, the best CA result is 98.92%, coming from the WLD method, which obtains the histogram information of the differential excitation and the gradient orientation, representing the salient local texture pattern, or in other words, the spatial structural characteristic as the proposed ISM-TPI method. On data set CF, the best classification result is 98.42%, coming from the PRICoLBP+SVM-based TPC method. However, the average CAs of ISM-TPI on the two data sets are only inferior as low as 0.1% and 0.38% than the best classification results from WLD and PRICoLBP+SVM, respectively. At the same time, CAs of the proposed ISM-TPI are higher than the other state-of-the-art methods, e.g., VZ algorithm, Path, and Path-MRF on the two data sets.


## 5.3.3. Discriminant Power of the SWDM-Based Feature

The main reason for achieving high classification accuracies is due to the stronger discriminant power of the SWDM-based texture features, which are demonstrated to be directly related to the perceptual properties of HVP system as explained in Section2.2. Hence, the proposed ISM-TPI method can obtain the essential spatial structural characteristics of texture images and can be adaptive to various texture patterns. In addition, the introduced SIGDFs and OAGDFs are effective filter banks for multiscale and omnidirectional spatial structural characteristics expression. Hence, SWDM-based features can achieve comparative and even superior classification performance compared with the state-of-the-art TPC methods.The SWDM parameters of two groups of typical confusing TPs are displayed in Figure14. D23, D27, and D31 represent the rubble-like texture pattern with different visual scales; the visual difference of these TPs can be expressed by the scale parameter β and location parameter μ, which is plot visually in the Figure 14(a). Regarding the TPs, D80, D82, and D84, whose surface organization of LHFs is resembling that may even mislead the visual identification of human beings. However, the essential distinctive SS of these TPs can be expressed clearly by the SWDM parameters λ,β,μ on different observation directions as can be seen in Figure 14(b).Figure 14
SWDM feature comparison of two groups of confusing TPs in Brodatz database, whose distinctive texture pattern characteristics can be expressed by the SWDM parametersλ,β,μ. (a) Rubble-like TPs, D23, D27, and D31. (b) Different TPs with great surface resemblance, D80, D82, and D84.
(a) (b)Some commonly used methods, e.g., GLCM, GWT, and MRMRF, cannot achieve high CAs on the challenging data sets mainly because that they may achieve similar feature parameters for the confusing TP analysis with great visual similarity. Figure15 displays four kind of texture feature extraction results of the confusing TPs, D80, D82, and D84, by the GLCM-based texture analysis method, GWT texture analysis method, WLD method, and VZ algorithm. The GLCM features are computed as follows. Image pixel gray-value of the original images are quantized into 8 gray levels and the GLCMs are computed under four orientation, 0°, 45°, 90°, and 135°, each orientation with 8 distances, 1,2,4,6,8,10,12,15; namely, we achieve 32 GLCMs for each ITP, and then the four GLCM-based statistics, contrast, correlation, energy, and homogeneity, are extracted and displayed in Figure 15(a).Figure 15
Feature extraction of TPs, D80, D82, and D84, based on GLCM, GWT, and VZ algorithm and WLD methods. (a) Original TPs, D80, D82, and D84, apparently with great visual similarity; (b) GLCM-based texture features; (c) GWT features, which tend to be the same trend in all the forty subbands; (d) histogram-based texture feature representation based on VZ algorithm proposed by Varma and Zisserman [11], with the first 100textons learned from the training image with MR8 filter banks class by class. PATCH and PATCH-MRF methods are their more recently texture classification methods, which has the similarity processing scheme but the original pixel patches are used and which can achieve the similar results for TP representation. (e) WLD histogram of D80, D82, and D84.
(a) (b) (c)  (d) (e)Regarding the GWT features, we compute the Gabor wavelet filters by the following construction formula:(45)hσ,θx,y;γ,fv,φ=π2σ2exp⁡-12x′2+γ2y′2σ2·cos⁡2πfvx′+φwhere x′=xcos⁡θ+ysin⁡θ, y′=ycos⁡θ-xsin⁡θ, σ is the standard deviation of the Gaussian factor (in pixels), θ represents the preferred orientation (in radians), φ denotes the phase offset (in radians), γ represents the spatial aspect ratio (of the x- and y-axis of the Gaussian ellipse), and 2π/fv means the preferred wavelength (period of the cosine factor, united in pixels).In the experiments, we choose 8 orientations in0.π, and five different wavelengths, namely, 2π/fv∈2,4,6,8,10,12, and we specify the γ as 2. Convolved with the forty Gabor filters, we can obtain the statistics mean and variance the texture images in each Gabor subband, which is plotted in Figure 15(c). The WLD features and the VZ algorithm features are histogram-based texture representation. The feature extraction steps fully follow the protocol in the literature [11, 12]. As can be seen clearly, GWT features, GLCM-based statistics, WLD histogram, and VZ algorithm-based texture features of D80, D82, and D84 exhibit strong similarity, some elements of which are even overlapped or only have very small difference, which seriously affected the discriminative ability of these methods for TPC.In a word, the state-of-the-art methods either compute some statistical measurement, e.g., GLCM, GWT, or make joint statistics of the image local structure with one-dimensional or two-dimensional histograms, whose resultant models are not based on the visual perception mechanism; hence, the ultimate results of the texture feature would have low discriminant power for the visual similarity texture patterns or make mistake in identifying the true TP of texture images sampled from large intravariant TPs, whereas the proposed SWDM reveals the essential LSD behavior of LHFs of the SS of TPs and the SWDM parameters under different observation scales and orientations can effectively obtain the distinctive visual features of TPs; hence the proposed ISM-TPI achieve relative higher classification performance, especially with better identification ability to the classification of the TPs with great visual resemblance or TPs with large visual variants.


## 6. Conclusions

The visual appearance of a texture pattern is visually determined by the random layout of LHFs in the imaged surface. The effective characterization of the LSD of LHFs in texture image, determined by the CSS of texture pattern, is the foundation for accurate TPC. This paper presents a TPC method from the image statistical modeling aspect, denoted as ISM-TPI, which mainly concerns the LSD of the CSS of TP. It has made a theoretical explanation of the Weibull distribution behavior of the LHFs in the imaging process by introducing the SFT, which consequently derives a SWDM to extract the distinctive SS characteristics of texture images for TPC. Omnidirectional image SS details under different observation scales are achieved by the local differential operation with two kinds of OGDFs with high computational efficiency, including OSIGDF and OOAGDF. The WDM features are demonstrated to be directly related to the HVP system, which can achieve high classification accuracies combined with the PLS-DA classifier. Compared with the state-of-the-art methods, the proposed ISM-TPI method can obtain higher and stable classification accuracy especially has good identification ability to TPC for the structurally similar and visually confusing texture images.



---

*Source: 6163787-2018-12-02.xml*