I Introduction
Ia Overview
Contentbased image retrieval (CBIR) is a special case of image classification. It can be viewed as the process of assigning a query image to a set of image classes, where each class represents a database image. Contentbased hereby refers to the mode of formulating the search query. As opposed to metadatabased image search, which relies on prelabeling the images beforehand, a CBIR system retrieves the best matches within the database with respect to the visual similarity to the query image.
In CBIR, the concepts of feature extraction (FE) and similarity measure (SM) play crucial roles. FE converts an image into a feature vector
of numerical values with the aim to produce a lowdimensional, yet sensible representation of the input image in the context of some particular application. An SM assigns a numerical value to a pair of two feature vectors. In this work, we assume that the SM is nonnegative and that a smaller SM value indicates a higher similarity and vice versa. A typical CBIR system is depicted in Fig.
1. The images to be extracted are processed via an FE algorithm and the extracted feature vectors (signatures) are stored in a database. The query image is fed to the same FE algorithm and the output is compared with all the feature vectors in the database by means of the defined SM. As a result, the images with the lowest SM values are returned. Consequently, designing a CBIR system boils down to the construction of an FE/SM framework.Nowadays, CBIR is employed to deal with big databases and vast amounts of data. Thus, in additional to high retrieval efficiency, execution speed is of great concern. This requires sufficient dimensionality reduction within the FE and a computationally efficient SM.
IB Related Work and Our Contribution
Extracting sensible information from images has been a field of research in digital signal processing for roughly half a century now. Naturally, the choice of the FE method depends on the specific application. For instance, when comparing images of landscapes, pictures taken in a forest will have different color distributions from those taken in a desert, so an FE based on the color properties of the image is a reasonable approach. On the other hand, images of objects are strongly characterized by shapes, so in this case the FE could rely on extracted contours. This work focuses on CBIR for textures. The earliest approaches for texture FE, such as the graylevel cooccurrence matrices (GLCM) [1] and the biologically motivated proposed Tamura features [2] have been presented in the 1970s and are often based on numerical measures for how nearby pixels relate to each other.
More recent approaches include incorporating spatialfrequency representations of the images [3], such as the Fast Wavelet Transform (FWT) [4], Gabor frames [5] or Curvelet frames [6]. A conventional approach is designing the features based on the energy of the respective transform coefficients, but it can be beneficial to look at their statistical properties [7, 8]. Thus, FE approaches based on the histograms of filterbank transforms have become more prominent, cf. [9, 10, 11, 12, 13, 14, 15]. Techniques of this kind are referred to as Subband Histogram methods in the work.
Linear filterbank transforms feed the input image to a bank of frequencyselective filters, yielding a set of bandpass signals as a representation. One problem of constructing FE based on such a decomposition is that the higher frequency subbands are prone to deformations in the spatial domain. Thereby, even slight deformations of an image will yield significantly different transform coefficients. However, simply neglecting the higherfrequency subbands is not desirable since they carry important information about such distinctive features as edges and corners. As a remedy, Mallat introduced the Scattering transform [16], a nonlinear signal representation involving filterbanks and the modulus operation which transforms highfrequency components into lowpass representations. The Scattering transform appears to perform well in texture classification tasks [17, 18, 19]. The key contribution of this work is to define a statistical model in order to design subband histogram FE for texture images based on their Scattering transform coefficients. Furthermore, we introduce a similarity measure based on the extracted features, with the property of being a Kernel and thus can be interpreted as an inner product of Textures transformed to some Hilbert space.
The paper is organized as follows. Section II provides necessary mathematical preliminaries of the problem of CBIR, subband histogram methods, and the basics of the Scattering transform. In Section III we develop a statistical framework of CBIR based on the Scattering transform. In Section IV several numerical experiments are presented to investigate the performance of the proposed approach in comparison with the state of the art methods.
Ii Mathematical Preliminaries
Iia Notations
Unless stated otherwise, a signal denotes an element of the Lebesgue space . The terms almost everywhere and almost all refer to conditions which hold for all , where can be any null set. For the sake of simplicity, we refer to as a Hilbert space of functions, even though it would be more precise to call them equivalence classes of functions that are equal almost everywhere. Boldfaced lowercase letters or describe vectors, while regular lowercase letters like or describe scalar values. Depending on the context, uppercase letters stand either for scalar values or for matrices. For a complex value , denotes its conjugate, and the real and imaginary part, respectively. Angle brackets denote the scalar product of and the canonical norm induced by it. An asterisk denotes the convolution of two signals .
IiB Subband Histogram Methods
The coefficients produced by a transform of the query image can be viewed as a set of realizations of a random experiment. Let describe the probability density function (PDF) of the transform coefficients of some database image.
The likelihood of a set of random realizations for a PDF is defined as
(1) 
and is a measure for the probability that is subjected to . That is to say, can be viewed as a measure of similarity between the query image and the database images. Since the natural logarithm is a monotone function, this is also true for the loglikelihood .
Let denote the set of the PDFs of a number of database images. The best match for a query can thus be determined via
(2) 
Expression (2) is called a Maximum Likelihood (ML) solution.
Assume the query is subjected to a PDF
. Then, by the law of large numbers, the logLikelihood
can be approximated by the negative crossentropy(3) 
Assuming a generative model, the crossentropy defines a similarity measure between the respective query and database image.
One of the most wellstudied transforms in image processing is the multiresolution decomposition produced by the FWT. In particular, it is known that the histograms of its bandpass subbands can be modeled by the Generalized Gaussian Distribution (GGD), c.f.[20], defined as
(4) 
where denotes the GammaFunction,
(5) 
and the parameter determines the scale of the distribution while describes the shape. The crossentropy between two GGDs ist given by
(6) 
Equation (6
) provides a parametrized SM for texture images. The GGD parameters can be estimated from each image using FE based on ML which together with the SM in (
6) defines a complete framework for CBIR. This approach was explored and discussed in [9], though it was equivalently formulated in terms of the KullbackLeibler divergence (KLD) rather than the crossentropy. It can be viewed as a blueprint for other subband histogram methods for which the motivation is twofold. First, the FWT is not the only filterbank transform for images and is not necessarily the best basis for texture FE. Furthermore, even though the crossentropy can be rigorously motivated by ML optimization, it lacks any geometrical interpretation, which poses the question, whether there are more suitable SMs for parametrized probability models.
IiC The Scattering transform
IiC1 Definition
Let be a rotationally symmetric window function with lowpass characteristics. Let and be fixed and a finite group of rotation matrices. With , we define the wavelet as
(7) 
Further, let us assume a lowpass and rotationally symmetric scaling function [21] and define
(8) 
such that for the respective Fourier transforms
,(9) 
holds, for almost all . Then the set
(10) 
constitutes a Parsevaltight frame [21]. It will be assumed to span in the following. Note that for two signals , the equation
(11) 
holds.
At the core of the Windowed Scattering transform (WST) [16] is a dyadic wavelet decomposition of the input signal with the complex modulus performed on the bandpass components, defined as
(12) 
The modulus operation traverses some of the energy of the bandpass signals towards lower frequencies. Therefore, can be applied to the output signals again. This yields a tree structure like the one depicted in Fig. 2. Note that (12) yields an infinite (but countable) number of output signals. However, in practice we deal with bandlimited input signals . Hence, there exists an integer with , such that Eq. (12) needs to be evaluated for only, which corresponds to a tree with finite breadth.
Basically, the idea of the WST is to apply successively to the input signal and only keep the lowpass signals, i.e. to neglect signals represented by the black nodes in Fig. 2.
The WST along the path of scaling factors and rotations is defined as
(13) 
Let us denote by the (ordered) set of all possible paths of size to . We can then define the Finitepath WST as
(14) 
IiC2 Properties
The output signals of the Finitepath WST are the elements of (14), also referred to as WST subbands in the following, which are all lowpass signals due to the filtering with . With increasing , captures more and more energy of the input signal and the WST representation becomes more expressive as the maximum path length grows. Let the Scattering norm be defined as
(15) 
As a consequence of the tight frame property (9) of , the infinitepath WST is unitary, i.e.
(16) 
In practice, maximum path depths of are expected to capture all essential information [22, 23]. Actually reconstructing from its WST involves a phase recovery problem. However, it is known [24] that digital realizations of wavelet representations similar to those employed in the WST are uniquely determined by their complex modulus.
The WST is nonexpansive [16], i.e. for two signals , it can be shown [16] that for any positive integer ,
(17) 
holds, implying that the WST is robust with respect to additive noise.
Since WST is a representation consisting solely of lowpass signals, it is stable with respect to small spatial translations and deformations. Specifically, let be a deformed or translated version of . Under certain mild assumptions about the underlying wavelet frame, it has been proven that the error is bounded asymptotically, cf. [16]. Note, that we are given some freedom of choice in the wavelet frame which allows for considerable flexibility in terms of parameters such as frequency selectivity or directionality.
Due to its stability and flexibility, the WST is a convenient signal representation and can be used as the foundation for the FE from images.
IiC3 Normalized WST
In order to reduce redundancy and to increase invariance to distortions, the Normalized WST (NWST) [23] was introduced. Let be the predecessor of , i.e. implies .
Let denote a very narrowband lowpass blurring filter. For the layers , the NWST is defined as
(18) 
In words, each subband of the WST is normalized by the respective parent subband, except for the subbands in the first layer which are normalized by the mean of the modulus of the input signal. In practice, a small constant is added to the denominator in order to avoid division by zero.
Iii Statistical Scattering CBIR
Iiia Subband Modeling
We propose to model the grayvalue distributions of the different WST subbands with parametrized PDFs and describe the images in terms of their respective parameters to obtain a complete FE mechanism on top of the WST.
The most distinctive features in textures are those of higher frequencies and are thus carried by the layers . These layers contain signals of the form
(19) 
Clearly, the modulus eliminates all negative values which is why a symmetric model such as the GGD is not appropriate anymore. The histograms of signals from descendant layers suggest a PDF that occupies only positive values and can describe a skew to the left, which makes the
Weibull Distribution (WD) a nearby choice. In fact, the WD was already successfully employed in the modeling of complex wavelet coefficients [15]. Fig. 3 exemplarily shows histograms of WST subbands of different textures with a fitted WD curve. Despite variations in skewness and spread, the WD fittings describe the actual histograms fairly well. The PDF of the WD is defined for as(20) 
Analogous to the GGD, is the scale parameter, whereas determines the shape.
The above argument is purely empirical. The following discussion aims to justify further the choice of the WD as a model for WST subbands (19) in the layers : Like for the layer , we neglect the impact of the lowpass filter . The approximately analytical bandpass filters produce signals with similar zeromean distributions in their real and imaginary parts. Moreover, it is known that bandpass components of a natural image can be well modeled by a GGD. This is also a reasonable assumption for bandpass components of Scattering subbands since they resemble natural images under dim lif´ghting conditions. It is therefore natural to assume the GGD with the same parameters as a model for the distribution of the real and imaginary parts of the bandpass signals involved in the WST. By neglecting statistical dependence between real and imaginary components, let us define a “complex GGD” with statistically independent real and imaginary parts as
(21) 
In order to determine the PDF for the modulus , we need to fix and integrate over the circles centered at 0 in the complex plane:
(22) 
Evaluating the integral analytically is demanding, if not impossible. However, it is known that the modulus of a complex variable with statistically independent and Gaussian distributed real and imaginary parts abides the Rayleigh distribution, which can be easily verified by substituting and into Eq. (22), i.e.
(23) 
The usage of the GGD for modeling FWT subbands was motivated by the need to model histograms which are more or less peaked in shape than it is the case for the Gaussian distribution. For this purpose, the additional shape parameter was introduced. This corresponds to introducing an additional shape parameter to the Rayleigh distribution. Altering the peakedness of the real and imaginary part corresponds roughly to altering the skewness of the respective modulus distribution. From Eq. (22), we can observe that , so a good alternative to for modeling the WST subbands would be a twoparameter generalization of the Rayleigh distribution which is controllable in its skewness without changing its value at 0. For , these requirements are fulfilled by the WD.
Even though it is sensible to model the WST with its Weibull coefficients, it is still questionable if this decision is justifiable for the NWST. It is certainly true for the first layer since it only involves an overall scaling. Unfortunately, this can not be assumed for the other descendant layers. Nevertheless, experiments show, that in practice this assumption still holds. However, for the most important ranges of , the multiplicative inverses of WD distributed samples exhibit histograms which again can be well modeled by the WD as can be seen in Figure 4. Thus, the normalized coefficients in (18) for are products of values whicht are close to be Weibull distributed and statistically independent. Thus the WD will again dominate the subband histograms.
We employ a numerical approach to estimate the Weibull parameters via minimizing the corresponding ML function. The derivation of the algorithm can be found in Appendix A.
IiiB Scattering Similarity
IiiB1 The Bhattacharyya Kernel
Both the crossentropy and the KLD can be derived for the WD in order to define an SM similar to the on discussed in Subsection IIB, cf. [25]. However, the derivation leads to very complex expressions.
As an alternative, we propose an SM based on the Bhattacharyya coefficient. For a pair of PDFs , it is defined as
(24) 
Expression (24
) is a popular choice for Kernels in the machine learning community
[26] and as such imposes a Hilbert space structure on probability based features.IiiB2 Weibull Similarity
Our aim is to derive a simple approximation of (24) for a pair of Weibull distributions . To this end, let us assume that . This is a justifiable assumption since the distributions are always compared for each subband individually. Let us further define
(25) 
This ultimately leads to
(26) 
For the sake of convenience, let us write the arithmetical and geometrical mean of two values
as(27) 
respectively. From the last equation of (26) we define our similarity measure for pairs of Weibull PDFs as
(28) 
In the derivation of , we replace and by each time it appears as the exponent of and , respectively. Even though this step is not necessary to eveluate the integral in (26), it ensures a straightforward understanding of (28), with the first factor measuring the similarity of the scale and the second factor measuring the similarity of the shape. Moreover, just like (24), the expression (28) defines a Kernel for Weibull PDFs, as it is shown in the following proposition.
Proposition 1.
Expression (28) is a Kernel for Weibull parameters.
Proof.
Note that for the approximation in (26) becomes an equality. I.e.,
(29) 
is a Kernel, since it is a Bhattacharyya coefficient of two PDFs. So is
(30) 
because both the shape and the scale parameter of the WD are constrained to be from the same domain . Therefore, Expression (28) is the product of two Kernels and as such a Kernel itself. ∎
The Kernel property of (28) provides us with some geometrical intuition, hence makes it a good choice as the fundamental building block of an SM. This includes the fact that it naturally induces a metric, which incorporates the notion of similarity as a geometric distance as is for instance the motivation behind multidimensional scaling [27].
It also enables the (N)WST to be conveniently subjected to many machine learning techniques such as Support Vector Machines.
IiiB3 Similarity for Multiple Subbands
Consider two PDFs
of two statistically independent random variables
,(31) 
It can be easily seen that
(32) 
Carrying over this insight to the derivation (26) yields, for a pair of sets of independent WDs with the parameter vectors ,
(33) 
For a pair of WST transforms with subbands, it is therefore straightforward to derive an SM. Since by our definition a low SM indicates a high similarity we apply the logarithm and reverse the sign which finally leads us to
(34) 
IiiC Implementing the Framework
The two previous sections provide the theoretical basis for building a complete CBIR based on WST Subband histograms. In order to extract the features, the WST (Section IIC) is performed on each image. Each subband of the layers is then subjected to ML in order to extract the WD parameters. If the number of subbands is , this amounts to feature vectors of size . In order to compare the query feature vector to the feature vectors in the database, the SM (34) us employed. As most of these steps can be performed independently of each other, the whole system is well suited to be implemented in a parallel manner, for example in applications for Big Data.
Iv Numerical Experiments
Iva Experimental Settings
IvA1 Implementation
All experiments were performed in a Matlab 2014a environment. The code reproducing the key results is available online^{1}^{1}1https://www.ldv.ei.tum.de/uploads/media/texture_retrieval_scattering_14.zip. The WST was performed by means of the ScatNet^{2}^{2}2http://www.di.ens.fr/data/software/scatnet/ toolbox, version 0.2. The Weibull parameters were extracted via wblfit from the Statistical Toolbox. Since it is the most standard subband histogram technique, we also implemented the FWT+GGD+KLD method according to [9] or comparison. Additionally, we implemented a method based the DualTree Complex Wavelet Transform (DTCWT), inspired by [15]. The DTCWT was performed by the DTCWT Transform Pack^{3}^{3}3http://wwwsigproc.eng.cam.ac.uk/Main/NGK, version 4.3. Similarly the subband histograms were modeled by the WD, but the KLD was replaced by the proposed SM in (34) for the sake of comparability.
IvA2 Choice of parameters
Along with the regular WST, all experiments were also conducted with its normalized version as defined in (18). The directionality parameter, corresponding to the cardinality of the rotation group , is fixed to be . All computations were performed with double precision and the subbands of the WST are critically sampled according to the ShannonNyquist sampling theorem. The bandwidth of the lowpass filter is chosen in a way such that it produces signals of size , in accordance with [9]. The WST is evaluated for the maximum path lengths and . In theory, a linear increase in leads to an exponential increase in the number of subbands. However, since the modulus operation hardly traverses any energy toward higher frequencies, they become negligible with increasing path length. Apart from the mentioned, default settings of ScatNet were used. In particular, the Morlet Wavelet was used as the bandpass filter on each layer of the WST.
The FWT relies on the D2 4tap Daubechies wavelet [21] as the underlying multiresolution representation. The DTCWT was run with the options antonini and qshift_a. In both cases, the size of the smallest subband was , which corresponds to three levels of decomposition.
IvA3 Data
The database D1 is the same as used in [9] and is based on the VisTex database^{4}^{4}4http://vismod.media.mit.edu/vismod/imagery/VisionTexture/distribution.html. It consists of 40 texture images of the size , each divided into 16 nonoverlapping patches. This amounts to 640 image samples of 40 different textures depicted in Fig. 5. The database D2 was generated from images of the following subset of the UIUC texture database^{5}^{5}5http://wwwcvr.ai.uiuc.edu/ponce_grp/data/: Bark1, Bark2, Wood2, Wood3, Water, Marble, Floor1, Floor2, Pebbles, Wall, Brick1, Glass1, Glass2, Carpet1, Carpet2, Wallpaper, Fur, Knit, Curdoroy, Plaid. Two images from each class are used to create five overlapping patches which are then scaled down to half the edge size. Hence, we get a database of 20 different texture classes each containing ten patches. Fig. 6 depicts one patch from each utilized image. All image patches are normalized to zero mean and unit energy, in order to avoid any bias caused by the overall lighting condition of each original texture image. The set of all patches generated from the same texture is considered a class. Its cardinality will be denoted by in the following. Consequently, for D1 and for D2. For each image patch, the most similar patches were retrieved. The retrieval rate for each patch is defined as the ratio of the number of retrieved patches from the same class to . The overall retrieval rate is the average of the retrieval rates for all the images in the database.
IvB Results
IvB1 Overview
The retrieval rates are summarized in Table I.








Database D1  75.50%  78.18%  78.93%  78.14%  84.90%  85.30%  
Database D2  52.39%  59.61%  66.50%  63.78%  66.94%  65.17% 


WDHMM  Rotated Wavelets  

76.93%  78.40%  80.05%  82.81% 
While WST+WD+ produces a similar result as DTCWT+WD+ on Database D1, NWST+WD+ is able to outperform all of the competing frameworks by for and for . Database D1 is widely used as a benchmark for CBIR retrieval. In order to provide a sense for the state of the art, Table II
summarizes the results from recent publications on comparable approaches: DWT + Generalized Gamma Distribution (G
D) [13], Wavelet Domain Hidden Markov Models (WDHMM)
[10] and Rotated Complex Wavelets [28]. Also, the original result for FWT+GGD+KLD from [9] was included, since we were not able to reproduce it. To the authors’ best knowledge, the best published result for this experiment so far was produced by Rotated Complex Wavelets with a retrieval rate of which is still outperformed by 2.09% by NWST+WD+ with and 2.49% with .Database D2 is considerably smaller than D1, but involves more variation within the classes, for instance, in terms of camera angle and deformation. Again, the WST greatly improves the retrieval performance. However, this time the regular WST does not fall behind the NWST. Also, increasing the maximum path length up to harms the performance. However, we can conclude that NWST+WD+ performs comparably well and produces the best results in both our test settings.
IvB2 Impact of Directionality
In general, natural textures contain strongly directional features. Hence, it is often beneficial to incorporate a frame of higher directionality as it is done for the Rotated Wavelets approach, for example. The DTCWT+WD+ method which mostly draws its benefits from increased directional selectivity in comparison with standard FWT based techniques is consistently and significantly outperformed by NWST+WD+, despite the same directionality properties. That is to say, there are reasons to believe that the success of the WST based methods is not exclusively due to increased directionality.
IvB3 Impact of Textural Structure
To further analyze the capabilities of the involved methods, consider Fig. 7. It shows the retrieval rates for each class in Database D1 produced by DTCWT as well as by WST and NWST with . DTCWT + GGD performed well in comparison to the WST based techniques for the classes Brick5, Buildings9, Metal0 and Wood2 and badly for Bark9, Fabric4, Leaves10 and Stone1. In general, it appears that DTCWT is most suitable for flat and (piecewise) smooth surfaces, while the two methods involving Scattering seem to be preferable for surfaces that are rather uneven, be it because they are rough or because they are bent or dented. This correlates with the premise that Scattering coefficients are less sensitive to spatial deformations, since on the one hand, images of rough surfaces carry significant fractions of highfrequency components and on the other hand spatial deformation in the image is often caused by locally bending or denting the depicted material. More specifically, NWST performed comparably well for the classes Bark0, Bark6, Flowers5, Terrain10 and regular WST for the classes Bark8, Fabric11, Grass1, Stone4. This suggests that WST works better on fine structures while NWST yields better results for coarse structures.
This assumption is corroborated by the following experiment. Fig. 8 depicts the retrieval rates for the Database D2 after blurring it with a Gaussian filter of different bandwidths. Unsurprisingly, the retrieval performance declines with stronger blurring for all methods. However, it is evident that the retrieval rate for WST declines significantly faster than the other two. This can be interpreted as a consequence of the distortion invariance of normalized Scattering coefficients. More importantly, this confirms that the regular WST method relies on fine features since they are most heavily affected by blurring. On the other hand, images of coarse patterns and flat and piecewise smooth surfaces are less affected by blurring such that the other two methods appear to be relatively robust.
V Conclusion and Outlook
In this work, we propose a subband histogram FE method based on the statistics of the WST of a texture image. Our derivation and analysis demonstrate how to extract statistical features from a WST representation of a texture image by means of matching a Weibull model via ML and how to measure the similarity of the respective feature vectors via the KLD. The proposed techniques come in handy when it comes to reducing the enormous degree of redundancy introduced by the WST since they represent each subband by as much as two real numbered parameters. The experiments show that the proposed method can outperform comparable algorithms based on filter bank transforms in terms of retrieval accuracy when designed as a CBIR system. Regular WST coefficients seem to work better on fine structures and Normalized WST coefficients on coarse structures.
Since approaches for rotation invariant Scattering representations are already available [17, 18], it appears feasible to extend the presented ideas towards a rotation invariant CBIR framework which would allow for more general problem settings.
The summation in (34) was partly motivated by assuming statistical independence of the subbands. However, neither does statistical independence hold for the WST subbands, nor for the subbands of other filterbank transforms, in practice. Accounting for statistical dependence in the model could lead to significant improvements in efficiency [12, 10].
a Estimation of Weibull parameters
In this section, we employ a numerical optimization approach to estimate the Weibull parameters. A more detailed derivation can be found in [29]. The loglikelihood amounts to
(35) 
The optimal parameters are identified by solving
(36) 
The first derivatives of with respect to and are computed as
(37) 
and
(38) 
Setting Eq. (38) to be equal to uniquely determines for . Solving for yields
(39) 
Substituting (39) in Eq. (37) and equating it with characterize the critical points of in terms of for , i.e.
(40) 
Since the function as defined in Eq. (35) is concave with respect to , a numerical algorithm such as the NewtonRaphson method can be used to obtain . Substituting this solution into Eq. (39) yields .
References
 [1] R. M. Haralick, K. Shanmugam, and I. Dinstein, “Textural features for image classification,” Systems, Man and Cybernetics, IEEE Transactions on, no. 6, pp. 610–621, 1973.
 [2] H. Tamura, S. Mori, and T. Yamawaki, “Textural features corresponding to visual perception,” Systems, Man and Cybernetics, IEEE Transactions on, vol. 8, no. 6, pp. 460–473, 1978.
 [3] T. Randen and J. Husoy, “Filtering for texture classification: a comparative study,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 21, no. 4, pp. 291–310, Apr 1999.
 [4] J. Smith and S.F. Chang, “Transform features for texture classification and discrimination in large image databases,” in Image Processing, 1994. Proceedings. ICIP94., IEEE International Conference, vol. 3, Nov 1994, pp. 407–411 vol.3.
 [5] B. S. Manjunath and W.Y. Ma, “Texture features for browsing and retrieval of image data,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 18, no. 8, pp. 837–842, 1996.
 [6] I. Sumana, M. Islam, D. Zhang, and G. Lu, “Content based image retrieval using curvelet transform,” in Multimedia Signal Processing, 2008 IEEE 10th Workshop on, Oct 2008, pp. 11–16.
 [7] G. Van de Wouwer, P. Scheunders, and D. Van Dyck, “Statistical texture characterization from discrete wavelet representations,” Image Processing, IEEE Transactions on, vol. 8, no. 4, pp. 592–598, 1999.
 [8] N. Vasconcelos, “On the efficient evaluation of probabilistic similarity functions for image retrieval,” Information Theory, IEEE Transactions on, vol. 50, no. 7, pp. 1482–1496, July 2004.
 [9] M. N. Do and M. Vetterli, “Waveletbased texture retrieval using generalized gaussian density and kullbackleibler distance,” IEEE Trans. Image Processing, vol. 11, no. 2, pp. 146–158, 2002.
 [10] ——, “Rotation invariant texture characterization and retrieval using steerable waveletdomain hidden markov models,” Multimedia, IEEE Transactions on, vol. 4, no. 4, pp. 517–527, 2002.
 [11] D.Y. Po and M. Do, “Directional multiscale modeling of images using the contourlet transform,” in Statistical Signal Processing, 2003 IEEE Workshop on, Sept 2003, pp. 262–265.
 [12] G. Tzagkarakis, B. BeferullLozano, and P. Tsakalides, “Rotationinvariant texture retrieval with gaussianized steerable pyramids,” Image Processing, IEEE Transactions on, vol. 15, no. 9, pp. 2702–2718, Sept 2006.
 [13] S. Choy and C. Tong, “Statistical wavelet subband characterization based on generalized gamma density and its application in texture retrieval,” Image Processing, IEEE Transactions on, vol. 19, no. 2, pp. 281–289, Feb 2010.
 [14] M. Allili, “Wavelet modeling using finite mixtures of generalized gaussian distributions: Application to texture discrimination and retrieval,” Image Processing, IEEE Transactions on, vol. 21, no. 4, pp. 1452–1464, April 2012.
 [15] R. Kwitt and A. Uhl, “Image similarity measurement by kullbackleibler divergences between complex wavelet subband statistics for texture retrieval,” in Image Processing, 2008. ICIP 2008. 15th IEEE International Conference on. IEEE, 2008, pp. 933–936.
 [16] S. Mallat, “Group invariant scattering,” Communications on Pure and Applied Mathematics, vol. 65, no. 10, pp. 1331–1398, 2012.

[17]
L. Sifre and S. Mallat, “Combined scattering for rotation invariant texture
analysis,” in
European Symposium on Artificial Neural Networks
, 2012.  [18] ——, “Rotation, scaling and deformation invariant scattering for texture discrimination,” in Computer Vision and Pattern Recognition (CVPR), 2013 IEEE Conference on. IEEE, 2013, pp. 1233–1240.
 [19] J. Bruna, “Scattering representations for recognition,” Ph.D. dissertation, CMAP, Ecole Polytechnique, 2012.
 [20] S. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 11, no. 7, pp. 674–693, Jul 1989.
 [21] ——, A Wavelet Tour of Signal Processing  The Sparse Way. Elsevier Science, 2008, pp. 89–376.
 [22] J. Bruna and S. Mallat, “Invariant scattering convolution networks,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 35, no. 8, pp. 1872–1886, 2013.
 [23] J. Andèn and S. Mallat, “Deep scattering spectrum,” IEEE Transactions on Signal Processing, 2013.
 [24] I. Waldspurger, A. d’Aspremont, and S. Mallat, “Phase recovery, maxcut and complex semidefinite programming,” Mathematical Programming, pp. 1–35, 2012.
 [25] C. Bauckhage, “Computing the KullbackLeibler Divergence between two Weibull Distributions,” ArXiv eprints, 2013.
 [26] T. Jebara, R. Kondor, and A. Howard, “Probability product kernels,” The Journal of Machine Learning Research, vol. 5, pp. 819–844, 2004.
 [27] W. Torgerson, Theory and methods of scaling. Wiley, 1958.
 [28] M. Kokare, P. K. Biswas, and B. N. Chatterji, “Texture image retrieval using new rotated complex wavelet filters,” Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on, vol. 35, no. 6, pp. 1168–1178, 2005.
 [29] D. Sornette, Critical Phenomena in Natural Sciences: Chaos, Fractals, Selforganization and Disorder: Concepts and Tools, ser. Springer Series in Synergetics. Springer, 2006, pp. 185–194.