Relevance of multifractal textures in static images

Electronic Letters on Computer Vision and Image Analysis 1(1):35-49, 2003 Relevance of multifractal textures in static images Antonio Turiel Air Project - INRIA. Domaine de Voluceau BP105 78153 Le Chesnay CEDEX. France Received 19 June 2002; accepted 28 February 2003 Abstract In the latest years, multifractal analysis has been applied to image analysis. The multifractal framework takes advantage of multiscaling properties of images to decompose them as a collection of different fractal components, each one associated to a singularity exponent (an exponent characterizing the way in which that part of the image evolves under changes in scale). One of those components, characterized by the least possible exponent, seems to be the most informative about the whole image. Very recently it has been proposed an algorithm to reconstruct the image from this component, just using physical information conveyed by it. In this paper, we will show that the same algorithm can be used to assess the relevance of the other fractal parts of the image. Keywords: Singularity analysis, wavelets, reconstruction. 1 Introduction Edge detection and texture classification are two main tasks in image processing, recognition and classification [1]. Extraction of edges provides information about the objects composing the scene, sometimes allowing segmentation; edges are thus the main source of information in the image and serve well also for classifying purposes. Texture information is more subtle, concerning the patterns and regularities inside the objects, light rendering and similar features. They also provide an important amount of information and they are specially useful in classification and segmentation tasks. One of the reasons to introduce the multifractal formalism in image processing was to provide a unified, reasonable way to deal with edges and textures at the same time [2]. The multifractal classification splits the image in edge-like and texture-like sets, which are arranged according to their properties under changes in scale (that is, under zooms). This approach is specially well adapted to certain types of images (for instance, those of turbulent or chaotic nature, as multifractality arose to explain the statistical properties of turbulent flows), but a great variety of real world scenes seem to be well described in this framework [3]. There is another reason to use the multifractal formalism: due to some statistical properties, one of the fractal components issued from the multifractal classification allows reconstructing the whole image. The implementation of the reconstruction algorithm has been recently proposed [4]. That reconstruction algorithm was designed to work over the most edge-like of the fractal components (reconstructing from edge-like structures has been explored in several contexts from scale-space theory [5] to wavelet analysis [6]). The key point is that the same algorithm can potentially be applied to the other components of the multifractal decomposition. Correspondence to: [email protected] Recommended for acceptance by Joachim Weickert ELCVIA ISSN:1577-5097 Published by Computer Vision Center / Universitat Aut` onoma de Barcelona, Barcelona, Spain 36 A. Turiel / Electronic Letters on Computer Vision and Image Analysis 1(1):35-49, 2003 The goal of this paper is to use this algorithm to evaluate the relative importance of each one of those fractal components. The paper is structured as follows: in Section 2, the theoretical fundations of the multifractal framework are briefly explained and the main implications discussed. Section 3shows how to apply the formalism in practice, in particular to produce the multifractal decomposition. In Section 4the reconstruction algorithm is presented and its properties are discussed; next, in Section 5we will make use of it to obtain an assessment about the relevance of each fractal component. Finally, in Section 6the conclusions of our work are presented. For the purposes of illustration, we will make use of Lena’s picture (Figure 1) and we will apply our techniques on it. The image presents remarkable deviations from the multifractal scheme (for instance, it has fuzzy edges in out of focus objects and numerous coding and processing artifacts), but however it is rather well described as a multifractal object. 2 Multifractal framework The multifractal formalism was developed first in the study of turbulent flows [7], as a way to explain the properties under changes of scale of very turbulent systems. It has been applied to the study of different types of images by several authors [8,2], as images have some properties which resemble to those of turbulent flows. We briefly sketch here the basic concepts in the approach we are going to use; for further details the reader is referred to [2]. We will denote any image by c(~x)where ~x denotes the vector coordinates of the referred pixel and it is normalized so that its average over the image vanishes, hc(~x)i~x∈image = 0. Acording to [2] we define a positive measure µas follows: for any subset Aof the image, its measure µ(A)is given by: µ(A) = ZA d~y |∇c|(~y)(1) that is, the measure assigns a weight to the set Aequal to the sum of the absolute variations of the image over it. Texturized areas will contribute with larger weights to the measure µthan flatly illuminated, smooth surfaces. In fact we will not be interested in the value of the measure over sets of fixed size, but in its evolution under changes in scale (resolution) around each point. Given a collection of balls Br(~x)of radii rand center ~x, we will say that the measure µis multifractal if: µ(Br(~x)) ≈α(~x)r2+h(~x)(2) for r’s small enough. The exponent h(~x)is called the local singularity exponent, and characterizes the way in which image behaves under changes in the size parameter rat the particular point ~x∗. As we consider small r’s, the largest values of the measures µ(Br(~x)) correspond to the smallest values of the exponents h(~x). For that reason, we will be specially interested in negative singularity exponents, which are found at pixels which contribute strongly to the measure by themselves (take into account that we consider very small radii). One of the advantages of this definition is that what determines the value of h(~x)is not the absolute variation of c(~x) at the point ~x, but its relative importance compared to the variations at the surrounding points: multiplying c(~x) by a constant modifies α(~x)in eq. (2), but lefts h(~x)unchanged. The classification of points accordingly is local, in opposition with global thresholding techniques. Natural images, that is, real word scenes of “natural” objects are of multifractal character [9,2], what has been tested for a large variety of scenes [3] and even with color images [10]. This property is far from trivial, and accounts for a special arrangement of edges and textures in images. In the following, we will only discuss on this type of images, although the same methods could be applied to other as well. Assessment of multifractality ∗The prefactor (2in our case) in the definition of the singularity exponent, eq. (2), is conventionally set to the dimension of the embedding space. This normalization allows to compare results from subspaces of different dimensions: the value of h(~x)becomes independent of the dimension of the space. A. Turiel / Electronic Letters on Computer Vision and Image Analysis 1(1):35-49, 2003 37 on real, digitized images can not be easily performed by a direct application of eq. (2) because of several technical reasons: some interpolation mechanism should be devised to take into account non-integer radii, for instance (there may be also undesiderable long-range effects which should be filtered; see [2] for a full discussion). In order to obtain a good evaluation of the singularity exponents, singularity analysis via wavelet analysis [11] should be performed. Wavelet analysis is a quite straightforward generalization of the scaling measurements in eq. (2): insted of applying the measure over finite size balls of radii r, a convolution of the measure µwith a scaled version of a wavelet Ψis computed. More precisely, the wavelet projection TΨµ(~x, r) of the measure µat the point ~x and the scale ris defined as: TΨµ(~x, r) = Zd~y |∇c|(~y)1 r2Ψ(~x −~y r)(3) The measure µis multifractal (in the sense of eq. (2)) if and only if: TΨµ(~x, r)≈αΨ(~x)rh(~x)(4) for small scale parameters r. Notice that αΨis in general dependent of the wavelet Ψand the measure µ, but the scaling exponent h(~x)has exactly the same value than in eq. (2) and does only depend on µ, that is, on the image c(~x)†. From the theoretical point of view, the choice of the particular wavelet Ψis irrelevant for the determination of the exponents h(~x); it can be even chosen as a positive function‡. However, in practical grounds there are wavelets which resolve better the finer structures than other. In Figure 2we show the representations of the multifractal classifications for four different wavelets. We will discuss further about the choice of the wavelet in Section 3. Multifractal classification of points is the first stage for multifractal decomposition of images (what justifies the name “multifractal” for the method). Points in the image can be arranged in fractal components, each one associated to a value for the singularity exponent. Namely, the fractal component Fh0associated to the exponent h0is given by: Fh0={~x ∈image |h(~x) = h0}(5) As the measure verifies to be multifractal, every point in the image can be associated a particular singularity exponent, so the image can be decomposed as the union of all its fractal components. They are indeed fractal sets [2], their dimensions being connected with statistical properties of images [12]. The most interesting of those fractal components is the Most Singular Manifold (MSM) [9], which is the fractal component associated to the least possible exponent. This set is usually related to the edges present in the image [2]. The least possible exponent is usually denoted h∞and its associated manifold Fh∞is generally denoted F∞in short. 3 Multifractal decomposition A correct determination of the MSM F∞implies a good multifractal decomposition, according to what was explained in Section 2. The main point concerns the choice of the analyzing wavelet Ψ. Once the wavelet is fixed, the singularity exponents are computed at every point in the image. The exponents are obtained by means of a log-log regression applied to eq. (4) at every point, in a range of scales typically going from 1 to 8 pixels non uniformly sampled (see [2]). Once every point is assigned a singularity exponent, the value of h∞ is estimated. A usual way to do this consists of taking the average between the values associated to the 1% and the 5% most singular points [2]. The dispersion around this value is conventionally fixed depending on the application. In Figure 2we present the functions h(~x)for four different wavelets. Let us define them. Let †Let us remark that the normalization in the wavelet Ψelliminates the prefactor 2in the exponent ‡Positive functions are not proper (admissible) wavelets (an admissible wavelet has zero mean [11]). Recall that the admissibility condition is necessary for signal representation (reconstruction), but not for signal analysis. 38 A. Turiel / Electronic Letters on Computer Vision and Image Analysis 1(1):35-49, 2003 ~x = (x1, x2)be the position vector, r=qx2 1+x2 2its modulus. We will make use of the following wavelets Ψi(~x),i= 1,2,3,4: 1. Lorentzian wavelet: Ψ1(~x) = 1 1 + r2 2. First radial derivative of Lorentzian wavelet: Ψ2(~x) = dΨ1 dr (~x) = −2r (1 + r2)2 3. Gaussian wavelet: Ψ3(~x) = e−1 2r2 4. Second radial derivative of gaussian wavelet: Ψ4(~x) = d2Ψ1 dr2= (r2−1)e−1 2r2 Each one of those wavelets fits the best for a particular application. Lorentzian wavelet (Ψ1) is a possitive wavelet of slow decay at infinity. It is very good to resolve sharp (negative) singularities in the measure µ (good spatial localization), but it has the backdraw of being unable to distinguish all the singularities beyond h= 0 (it returns the value h= 0 for all of them); besides, it cannot be used to analize the signal c(~x)directly (a certain number of vanishing moments would be required [13,2]). The gaussian wavelet (Ψ3) cannot be either used over the signal itself, as it is positive also, but having fast decaying tails it is able to resolve the whole range of singularities (typically between −1and 2, see [2]); the backdraw is a worse spatial localization, specially for the MSM. The second derivative of the gaussian (Ψ4) is, from the theoretical point of view, the best possible choice for analyzing signals: it resolves the whole range of values of h(~x)and it can be even used over the signal itself, without necessity of constructing a measure. However, in practice it has very poor spatial localization, associated to an inner minimum scale of several pixels, necessary to separate positive from negative extrema in wavelet projections. The best choice in practice is then the derivative of Lorentzian wavelet (Ψ2), which arrives to a compromise in range of detected singularities, localization and applicability over the full signal. It is not well adapted for any one of those tasks (it truncates the range of singularities above h= 1, it blurs localization, it has not enough number of vanishig moments), but it is able to provide meaningful results in every context. In Figure 3several different fractal manifolds for our image are represented, every column showing the sets associated to each one of the wavelets discussed above. The first step is to compute h∞as described at the beginning of this Section, obtaining the different values for the different wavelets: h∞=−0.47 for Ψ1, h∞=−0.32 for Ψ2,h∞=−0.43 for Ψ3and h∞=−0.68 for Ψ4. As a general remark, wavelets with higher orders of derivative are more imprecise in the determination of this value, while positive wavelets throw more similar results. Once the value of h∞has been obtained, we isolate the MSM, defining it as the set of points ~x for which h∞−∆h≤h(~x)< h∞+∆hwith a conventionally fixed value of the dispersion ∆h; in the following we take ∆h= 0.15. We represent also the other fractal manifolds according to the given dispersion, so the nth MSM will be the set of points ~x for which h∞+ (2n−3)∆h≤h(~x)< h∞+ (2n−1)∆h(the MSM itself is the first MSM). Finally, we define the manifold of excluded points or excluded manifod as the set of points ~x such that h(~x)< h∞−∆h, that is, which are more singular than expected. Those events are generally associated to the borders of the image and some particular events, which happen to have singularities close to −1, typical to isolated edges [2]. In Figure 3, we show all those manifolds. A. Turiel / Electronic Letters on Computer Vision and Image Analysis 1(1):35-49, 2003 39 4 Reconstructing from edges Recently, an algorithm to reconstruct the whole image from the most singular of its fractal components has been proposed [4]. We will not go into details about the reconstruction algorithm; we will just present the final formula and discuss it. The reader is referred to the original paper. The reconstruction formula intends to reproduce the whole image from the value of the gradient field over the MSM. First, let us define the essential gradient over a general set F. We define it as a vector function which is only different from zero over the set F, namely: ~vF(~x)≡ ∇c(~x)δF(~x)(6) where the symbol δFstands for a delta function on the set F. The reconstruction algorithm is given by the following expression: c(~x) = ~g ⊗~vF∞(~x)(7) where ⊗stands for the convolution and the reconstructing kernel~g is given in the Fourier space by the following expression: ˆ ~g(~ f) = i ~ f f2(8) In the above expression, the symbolˆstands for the Fourier transform, ~ fis the spatial frequency (the variable in the Fourier domain) and i≡√−1. The reconstruction formula states that it is possible to retrieve the image from the essential gradient associated to the MSM F∞. Note, however, that the formula could be applied to any set F; we will denote by cFthe image retrieved from the essential gradient associated to the set F; namely: cF(~x) = ~g ⊗~vF(~x)(9) We will call eq. (9) the generalized reconstruction formula. In this language, the reconstruction formula states that cF∞=c. The generalized reconstruction formula has some nice properties. •It is linear in the reconstructing data: If the set Fis the disjoint union of two sets F1and F2(i.e., F=F1∪F2, with F1∩F2=∅), then cF=cF1+cF2. This comes from the fact that ~vF1∪F2=~vF1+~vF2if the sets are disjoint, and the associativity of the convolution product. •It always exists a set from which reconstruction is perfect: If F=<2, that is, the whole image, ~vF=∇c, but as ˆ ∇c(~ f) = −i~ fˆc(~ f)and taking into account the definition of ~g, trivially cF=c. Taking into account both remarks, we conclude that if Fcis the complementary set of a set F,cF+cFc=c, which can also be expressed as c−cF=cFc, that is, the reconstruction from the complementary of Fis equivalent to the error image (the difference between the reconstruction and the actual image). The reconstruction formula states that there exists a rather sparse set F∞from which the reconstruction is perfect (equivalently, the reconstruction error is zero). In practice, however, a good determination of F∞is sometimes difficult. In such cases, the generalized reconstruction formula allows measuring how relevant the points not included in that set are, for instance just measuring the PSNR’s for the reconstructed images. Due to linearity, the same measure can be interpreted as the decrease in the error associated to the inclusion of those points in the estimate of F∞. We apply those ideas in the next section to interpret the importance of the different fractal components extracted according different wavelet projections. 40 A. Turiel / Electronic Letters on Computer Vision and Image Analysis 1(1):35-49, 2003 5 Relevance of the fractal manifolds We will make an assess about the relative importance of the fractal manifolds by means of the generalized reconstruction formula. In Figure 4we show the different images reconstructed from the manifolds presented in Figure 3using eq. (9); in Table 1the associated PSNRs can be found. We see that the MSM provides always the greatest amount of information about the image, which is reflected both by visual inspection and the values of the PSNR. However, the second manifold contains a significant amount of information, which reflects in the recogniscible structures which are reconstructed from it and still significant values of PSNR (in the case of the Lorentzian wavelet (first column), the second MSM contains all the other points, because it is not able to distinguish singularities above h= 0 and they are truncated to that value). The other manifolds (when they can be distinguished, that is, when considering wavelets other that Lorentzian) contain significantly very few information. The excluded manifold deserves a particular comment. It contains very sharp edges and it accounts for global illumination conditions (for instance, more light over the hat or the shoulder, the global focus on the right of the image,...). It should be included in any reasonable reconstructing set, even if their statistical properties (due to boundary conditions or strong fluctuations) may constitute a deviation from the multifractal model (as they more singular than what is predicted by the model). The relative importance of each fractal part can be better understood looking at Figure 5, in which images in Figure 4are progressively summed up from top the column to the bottom; due to linearity of eq. (9), the resulting images are equivalent to the result of reconstructing from the succesive union of manifolds in Figure 3. It is obvious from visual inspection that after the second manifold very few information is incorporated in the successive additions, a fact also evidenced by the associated PSNRs, Table 2. 6 Conclusions In this paper we have recalled the multifractal formalism, which stands to be a method for classifying points in images according to their singular character. We have seen that this rather mathematical characterization (the singularity exponent) has an interpretation in terms of relative informative relevance: the most singular points are the most informative about the scene. This characterization of the informational content is made by means of the reconstruction algorithm [4], which was proposed as a way to reconstruct images from edges, derived from simple, general assumptions. The properties of the reconstruction algorithm allow to isolate the contribution of every point in the final reconstructed image. We have made use of it to assess the qualities as reconstructing sets of the different fractal components spawned in the multifractal scheme. The method proposed here could be used to determine which properties (edges, textures) are important to keep in order to have a good visual performance in compressed images and which ones could be removed without affecting significantly the quality. It is a rather natural technique, as it is based on physical properties of images. It is important to notice that the reconstruction algorithm can be considered an edge-detection based coding scheme, much in the way of the modern techniques of ridgelets and curvelets [14], which have been shown to be very efficient for image coding. In order to implement compressing techniques using the reconstruction algorithm, high performance reconstructing sets should be extracted from images. The technique of singularity classification is a good first approach to obtain that set, but the multifractal model is just approximate for general real word images (it was derived for a subset of so-called natural scenes) and so the MSM is just an approximation to the best reconstructing set. Besides, singularity detection is a complicated technique, which requires fine tuning in the choice of the analyzing wavelet. In spite of all those backdraws, it is possible to obtain good performance just using singularity analysis (see Figure 6). From a more general perspective, however, other methods for the extraction of the reconstructing set need to be devised. A. Turiel / Electronic Letters on Computer Vision and Image Analysis 1(1):35-49, 2003 41 Acknowledgements A. Turiel is financially supported by a post-doctoral grant from INRIA. I am grateful to J. Grazzini for his comments and suggestions. References [1] D. Marr, Vision, W.H. Freeman and Co., New York, 1982. [2] A. Turiel and N. Parga, “The multi-fractal structure of contrast changes in natural images: from sharp edges to textures,” Neural Computation, vol. 12, pp. 763–793, 2000. [3] A. Nevado, A. Turiel, and N. Parga, “Scene dependence of the non-gaussian scaling properties of natural images,” Network, vol. 11, pp. 131–152, 2000. [4] A. Turiel and A. del Pozo, “Reconstructing images from their most singular fractal manifold,” IEEE Trans. on Im. Proc., vol. 11, pp. 345–350, 2002. [5] R. Hummel and R. Moniot, “Reconstructions from zero crossings in scale space,” IEEE Trans. on Acoustics, Speech and Signal Processing, vol. 37, no. 12, pp. 2111–2130, December 1989. [6] S. Mallat and S. Zhong, “Wavelet transform maxima and multiscale edges,” in Wavelets and their applications, Ruskai M. B. et al, Ed. Jones and Bartlett, Boston, 1991. [7] U. Frisch, Turbulence, Cambridge Univ. Press, Cambridge MA, 1995. [8] J. L´ evy-V´ ehel, “Introduction to the multifractal analysis of images,” in Fractal Image Encoding and Analysis, Y. Fisher, Ed. Springer Verlag, 1997. [9] A. Turiel, G. Mato, N. Parga, and J. P. Nadal, “The self-similarity properties of natural images resemble those of turbulent flows,” Physical Review Letters, vol. 80, pp. 1098–1101, 1998. [10] A. Turiel, N. Parga, D. Ruderman, and T. Cronin, “Multiscaling and information content of natural color images,” Physical Review E, vol. 62, pp. 1138–1148, 2000. [11] I. Daubechies, Ten lectures on wavelets, CBMS-NSF Series in Ap. Math. Capital City Press, Montpelier, Vermont, 1992. [12] G. Parisi and U. Frisch, “On the singularity structure of fully developed turbulence,” in Turbulence and Predictability in Geophysical Fluid Dynamics. Proc. Intl. School of Physics E. Fermi, M. Ghil, R. Benzi, and G. Parisi, Eds., Amsterdam, 1985, pp. 84–87, North Holland. [13] A. Arneodo, F. Argoul, E. Bacry, J. Elezgaray, and J. F. Muzy, Ondelettes, multifractales et turbulence, Diderot Editeur, Paris, France, 1995. [14] E. J. Cand` es and D. L. Donoho, “Curvelets, multiresolution representation, and scaling laws,” in Wavelet Applications in Signal and Image Processing VIII, A. Aldroubi, A. F. Laine, and M. A. Unser, Eds., 2000, vol. 4119 of Proc. SPIE. 42 A. Turiel / Electronic Letters on Computer Vision and Image Analysis 1(1):35-49, 2003 7 Lena’s picture Figure 1: Lena’s image A. Turiel / Electronic Letters on Computer Vision and Image Analysis 1(1):35-49, 2003 43 8 Singularity exponents with four wavelets Figure 2: Multifractal decompositions on Lena’s image for Lorentzian wavelet and its derivative (top) and Gaussian wavelet and its second derivative (bottom) (see Section 5). The smaller is the singularity exponent, the brighter is the point.