Multifractal wavelet filter of natural images

Natural images are characterized by the multiscaling properties of their contrast gradient, in addition to their power spectrum. In this work we show that those properties uniquely define an {\em intrinsic wavelet} and present a suitable technique to obtain it from an ensemble of images. Once this wavelet is known, images can be represented as expansions in the associated wavelet basis. The resulting code has the remarkable properties that it separates independent features at different resolution level, reducing the redundancy, and remains essentially unchanged under changes in the power spectrum. The possible generalization of this representation to other systems is discussed.


I. INTRODUCTION
Given the complexity and degree of redundancy of natural images, the early visual system had to find good coding procedures to represent the visual stimuli internally. To achieve this goal, the visual system must have learnt the regularities present in the environment where the organism lived [1]. If there are image features that tend to appear together, a cell responding quasi-optimally to them is rather likely to exist. To find such a representation one has first to understand the statistical properties of visual scenes common in the environment. In particular, the relevance of the second order statistics has been pointed out some time ago [2], and internal representations that eliminate these correlations have been discussed [3]. However, even if whitening represents an improvement of the code, it still leaves much geometrical structure that should be dealt with more properly [4]. A more systematic study of statistical regularities that go beyond the two-point correlations has began rather recently [5][6][7]. A novel approach to understand the statistical properties of natural images has been proposed in [6] where the non-gaussian statistics of changes in contrast has been characterized and explained by means of a stochastic multiplicative process: Contrast changes at a given scale are obtained from those at a coarser scale by multiplication with an independent random variable. This implies a linear relation between the logarithms of the variables at two different scales (this property has also been discussed in [9], although the existence of a multiplicative process was not noticed). A very rich geometrical structure has emerged from those studies: contrast changes are organized in such a way that pixels in the image can be classified according to the strength of the singularities of the contrast gradient [7]. It was also checked that the multiplicative process is present in very different sets of images and in color natural images [8].
In this work we show that when those findings are taken into account in an image model there appears a very compact representation of the visual world. In particular, a wavelet filter that guarantees that variations in contrast at different scales are related by a multiplicative process can be derived experimentally from a dataset of natural images, with the only assumption of the existence of such a filter. Furthermore, the code that it defines has eliminated a great deal of redundancy.

II. WAVELET BASIS AND MULTIRESOLUTION ANALYSIS
We start by considering the projection of the contrast c( x) on a dyadic wavelet setΨ j k ( x) =Ψ(2 j x − k). The wavelet Ψ j k focuses on the details of the image at the scale r, where r = 2 −j , j ∈ Z, at the sampling points x 0 = 2 −j k. Here k ≡ (k 1 , k 2 ), with k 1 , k 2 ∈ Z [10]. IfΨ defines a wavelet basis, the wavelet projections T r Ψ c( x 0 ) of the field of luminosity contrast, characterize completely the image. This is the reason of the name "multiresolution analysis": the wavelet projection is a description of c( x) at the point x 0 when the image is observed at a variable scale r. This type of analysis reaches a compromise between localization in position and in spatial frequency [10]. If the discrete basis is complete, then the contrast can be expanded in a wavelet basis orthogonal toΨ using the wavelet projections as coefficients. The dual basis Ψ j k verifies that: and the contrast c( x) can then be expressed as: where α j k = T 2 −j Ψ c(2 −j k) = 2 2j Ψ j k |c . The wavelet basis generated by Ψ will be called the representation basis. Its dual basis, defined byΨ will be referred to as the analizer basis. Once one of these basis is known, the other is completely defined by eq. (2).

A. Multiscaling in natural images
In [7] it was experimentally proved that there exists a certain class of wavelets, such that some of their low order moments are zero, for which: (the angular brackets denote the average over an ensemble of images), a property known as Self-Similarity (SS). It was also found [6] that the SS exponents τ p have a non-trivial dependence on p which can be explained by means of a stochastic multiplicative process (see e.g. [11]). In terms of the wavelet coefficients α j k , this process is expressed as: where [ κ] denotes the vector with components given by the integer part (rouding down) of those of κ. In [6,7] eq. (5) was understood in the statistical sense, where the variables η j k are statistically independent of the α j−1, k 2 . The SS exponents τ p define completely the distribution of the random variables η j k and vice-versa. Besides, the distribution of the η's depends only on the ratio between the two scales of the wavelet projections [6]. Here this ratio has been fixed at 2 for any pair of consecutive scales. This implies that all the η j k 's are identically distributed.
In this work we go beyond those results by showing that natural images posses an intrinsic wavelet for which eq. (5) is fulfilled point by point. This means that the equality holds for any image, at any scale and position and that the variables η j k can be extracted directly from: In order to obtain these variables and to verify their statistical properties, one has first to determine this wavelet. Under the assumptions that the η's obtained from eq. (6) are scale independent and equally distributed variables, the representation wavelet Ψ can be experimentally obtained from a statistical analysis of the image ensemble. This will be shown in the next section. The validity of these two hypothesis on the η's has to be verified a posteriori, once the wavelet is known. This is done in Section IV as follows: first the the analizer waveletΨ is obtained from the representation wavelet Ψ. In turn,Ψ can be used to evaluate the coefficients α j k 's and from them the η j k 's. Once these are known it is finally checked that they are indeed scale-independent, identically distributed random variables, so checking the self-consistency of the multiplicative process model.

III. THE REPRESENTATION MOTHER WAVELET
We suppose that the contrast field obtained from our image dataset can be expanded as a superposition of wavelets, eq.(3), and that the α j k 's verify eq. (5), where the η j k 's are scale independent equally distributed variables. Since the images have finite size we will take j ≥ 0 where Ψ 0 0 covers the whole image and it represents the mother wavelet, Ψ 0 0 ≡ Ψ. For the same reason, the range of k = (k 1 , k 2 ) at the scale j is bounded as: k 1 , k 2 = 0, 1, ..., 2 j − 1.
Averaging c( x) at each point x, it is found that the average contrast can be represented as a simple wavelet superposition: Here |η| is the first order moment of the distribution of the |η|'s, and we have used the assumptions that all the η j k 's have the same marginal distribution and are independent across the scales. By Fourier transforming this field, C( f ), one easily obtains the Fourier transform of the representation wavelet,Ψ( f ), that reads: where Λ( f ) = (1 − e −2πif1 )(1 − e −2πif2 ). This expression is very appealing. The right hand side compares the average contrast at two consecutive scales (related by a factor 2), and it expresses that the wavelet is obtained as an observation of the scale transformation properties of the images. The average |η| has an a priori known value: |η| = 1 2 [7]. This allows to use eq. (7) with no a priori knowledge about the data.
To obtain experimentally the representation wavelet, eq. (7) was applied to a large ensemble of natural images. The data were 200 1024 × 1024 images selected at random from van Hateren's dataset [12]. Figure 1 shows the representation wavelet obtained with this procedure. It exhibits two clear features: it is right-left symmetric and roughly up-down antisymmetric with a sharp central discontinuity (see Figure 2). The first property is expected: our world remains statistically unchanged when it is reflected from left to right. The rough up-down antisymmetry and the related discontinuity are probably due to the sharp contrast between the sky and the ground.

IV. THE PROJECTION INTO THE EXPERIMENTAL BASIS
Once the analizer wavelet is obtained, it can be used to evaluate the coefficients α j k and from them, using eq. (6), the coefficients η j k . Then, it is possible to check whether the η j k are scale-independent, identically distributed variables. In the affirmative case, this fact selfconsistently demonstrates that eq. (7) is the intrinsic wavelet we are looking for. It is immediate to check that the signs of the η's are independent of their absolute values, |η|, and also scale-independent, so it is only necessary to verify the independence of the |η|'s.
The mutual dependence among the |η j k |'s was estimated by computing the correlation coefficients between two |η|'s. As the number of possible pairs is very large, only two types of correlations were considered, which should be maximal by construction: those between consecutive scales (j − 1 and j) at equivalent positions, and those between consecutive spatial positions ( k and k + d, where the components of d take only the values 0 and 1) at a fixed scale j. Both correlation coefficients, denoted as ρ j and ρ j d respectively, take values in [−1, 1] and give a dimensionless measure of the degree of statistical dependence. The observed values of ρ j are very close to 0 (|ρ j | < 10 −2 , j > 2), what confirms rather well the independence of the η j k 's under changes of scale. On the contrary, the correlation coefficients ρ j d are by no means negligible, although they are extremely short-ranged: the variables η j k and η j k ′ are independent when they are more than one pixel apart; after that distance the two-point correlation ρ j d decays dramatically. It is important to remark that there is no need of spatial independence of the η j k 's in our wavelet model. For this reason, although short-ranged, the observed dependence is a significant source of information about the remaining statistical structure at a fixed resolution layer. On the other hand, the η's define a system of almost completely independent variables.
To confirm the validity of the wavelet representation, eqs. (3) and (7), one has still to check that the η's are identically distributed. Figure 3 exhibits the distribution of ln |η j | for different scales j. The correspondence between them is really very good.

V. WAVELET TRANSPARENCY AND DECORRELATION
It has been frequently argued that the the signal that arrives to the primary visual cortex has been already decorrelated in previous stages of the visual system [3,13]. In this case, V1 would take care of coding more complex aspects of images. In this regard, the way in which the mother wavelet is constructed guarantees a remarkable property of transparency to the power spectrum: it defines a code that is somewhat independent of the second order statistics of the images. The power spectrum of natural images exhibits a power law behaviour, S( f ) ∼ f −(2−ǫ) [2]. Under the assumption of translational invariance, the application of the decorrelating filter to the contrast is equivalent to multiplication in the Fourier domain by f 1− ǫ 2 . Denoting the decorrelated contrast as Dc( x), it is immediate to see that it has a representation similar to eq. (3): where DΨ indicates the application of the decorrelating operator to the representation wavelet and α ′ j k = 2 j(1− ǫ 2 ) α j k . Defining now η ′ j k = 2 1− ǫ 2 η j k it is concluded that the decorrelated images also posses a random multiplicative process. The new representation wavelet DΨ can be obtained from an ensemble of decorrelated images by means of eq. (7). It is known that the value of the exponent ǫ has fluctuations from image to image [14]. The relevance of the wavelet transparency is that if the visual stimulus has been decorrelated before it reaches V1, this area can use the type of filters described in this paper regardless of the precise exponent of the power spectrum.

VI. CONCLUSIONS
We have shown that images can be represented as multiresolution objects in terms of an appropriate wavelet basis Ψ, in which each resolution level is an independent image. One of the advantages of this representation is that it is based on the observed properties of the contrast gradient [6,7], what in turn leads to an automatic reduction of the redundancy. At the same time the spatial correlations at a given scale are short-ranged, but still informative.
Once the wavelet is known, it is possible to devise a compression technique based on the properties of this representation. This would have two important advantages with respect to other wavelet representations (see e.g. [15]): first that the scale layers are independent, and second that there exists a simple model for the distribution of coeficcients [6][7][8] that can be used in the coding.
These results have been obtained with the simplest wavelet expansion, but the tools presented in this work can be taken as the starting point to look for more realistic visual filters. This search should be directed by the experimental observation that cells in V1 are edge detectors. From this perspective, the expansion in eq. (3) should be generalized to include an orientational degree of freedom. Filters of this type have been proposed in [16] and found from an independent component analysis of natural images [17] . These studies should be combined with the use of overcomplete basis, this introduces redundancy but the representation becomes stable under small changes in the images [18] This analysis could also be carried out over any physical system with multiscaling properties, eq. (4) (for instance, turbulent flows in fully developed turbulence [11]). For all these systems eq. (7) could be applied to obtain the associated representation wavelet Ψ. However, it would be necessary to check the statements of scale independence and identical distribution for the projection coefficients η j k . If both properties hold, the compact code so obtained would be a valuable tool, and the interpretation of the possible remaining structure at each scale very meaningful.