- Email: [email protected]
Synthesis of bidimensional α-stable models with long-range dependence
Signal Processing 82 (2002) 1927 – 1940 www.elsevier.com/locate/sigpro
Synthesis of bidimensional -stable models with long-range dependence B%eatrice Pesquet-Popescua; ∗ , Jean-Christophe Pesquetb a Dept.
Traitement du Signal et des Images, Ecole Nationale Superieure des Telecommunications, 46, rue Barrault, 75634 Paris Cedex 13, France b Institut Gaspard Monge and URA CNRS 820, Universit e de Marne-la-Vallee, 5, Boulevard Descartes, Champs sur Marne, 77454 Marne la Vallee Cedex 2, France Received 5 May 2001; received in revised form 9 February 2002
Abstract In the context of linear modeling, the main advantage of stable distributions is to allow the de2nition of non-Gaussian processes whose statistical properties are easy to characterize. In this work, we are interested in the design of a speci2c class of 2D discrete-space processes with stable distributions. A frequency domain method for the synthesis of these 2elds will be proposed which is similar to algorithms already used in the Gaussian case. The considered models represent 2D -stable extensions of the 1D fractional Gaussian noise. They exhibit long-range dependence properties and, consequently, they could provide interesting alternatives to image modeling techniques based on the 2D fractional Brownian motion. However, the conditions for the existence of such processes deserve special attention and they are derived in this paper. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Stable process; Long-range dependence; Bidimensional -stable models; Fractional stable noise; Texture modeling; Image synthesis
1. Introduction In the last years, signal processing using fractal models has represented an active research area stimulated by the plethora of applications (compression, computer graphics [16,25,24] geophysics [21], traf2c modeling in computer networks [2], turbulence phenomena [1,3], 2nance, etc.) where these concepts are potentially interesting. Numerous works on “fractals” have popularized the notion of “self-similarity”.
∗
Corresponding author. E-mail addresses: [email protected] (B. Pesquet-Popescu), [email protected] (J-C. Pesquet).
This concept is at the origin of relevant models for several natural phenomena. In the meantime, the non-stationary structure of self-similar processes put them at the core of the most recent preoccupations in signal and image processing. Stable processes have turned out to be good models for many impulsive signals and noises [17], when the probability distributions of the highly variable data have “heavy” tails. These distributions have in2nite variance and unde2ned higher-order moments, but it was pointed out in [36] that many signal processing algorithms based on second-order statistics can be transposed to fractional lower-order moments. In this paper, we present new approaches for generating fractal random 2elds and, by this way, we
0165-1684/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 2 ) 0 0 3 2 0 - 1
1928
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
generalize the results existing in the 1D case. The resulting textures can be further used in classi2cation or recognition tasks, as well as for numerical terrain models. In this later application, the fractal characteristic allows to realize more realistic processing. In the case of an interpolation for example, this leads to better results than those obtained by traditional methods relying on autoregressive or deterministic models [30]. The proposed models are interesting for modeling and synthesizing textures with impulsive and long-range dependence (LRD) behaviors. Related existing works concern mainly the Gaussian case, i.e., two-dimensional versions of the fractional Brownian motion (fBm) [8,11,14,15,20,26,27,34,38] and its discrete counterparts [28,31]. In works like [8,14,15] these models have been used for texture synthesis, while in [11,38], they have been exploited for segmentation of synthetic or satellite images (in the later case, the diGculty of computing texture metrics in high-speckle SAR imagery being one of the motivations for using fractal models). The strong relationship between LRD models and fractals (when the Hurst parameter H is such that 1 2 ¡ H ¡ 1) has been emphasized in several papers [12,22], and since then a lot of applications of these models have been considered, related to 1D and 2D phenomena. In [13], 2D isotropic FARIMA processes driven by one-sided exponential noise are used to model the sea clutter texture in SAR images collected by RADARSAT sensors. In [19], a 1=f -type spectral behavior is exhibited for ultrasound RF echoes, and the parameters of this model are related to tissues characterization features. The authors also point out the usefulness of a 3D model, but, for the sake of simplicity, they built a 1D model of the backscattered signal. Potential applications of 2D -stable models can be found in diHerent 2elds. It seems that three applications could be mainly concerned in image analysis, due to the impulsiveness of the underlying data: SAR imaging, ultrasound medical imaging and astronomical imaging. The impulsiveness of certain classes of images has been modeled by -stable 2elds [18,32] or by L%evy processes [33]. In this paper, we propose a 2D linear fractional process driven by a stable noise, that can be a useful model in the kind of applications cited above. In the context of linear modeling,
the main advantage of stable distributions is to allow us to de2ne non-Gaussian processes whose probability laws are easily deduced from those of their driving noise. In this work, we will be interested in the design of a speci2c class of 2D discrete-space processes with stable distributions. Furthermore, the considered models have LRD properties. Conditions for the existence of such long dependent -stable 2elds are provided, that lead to the de2nition of the validity domain for the fractional parameter. An FFT-based synthesis method is also proposed, which is close to other synthesis methods developed in the Gaussian case [14]. In order to illustrate the potential interest of the proposed model, we also provide preliminary results suggesting that they can be used to simulate speckle noise in SAR images and the cosmic background in astronomical images. The outline of this paper is as follows. In Section 2 we recall some useful facts on stable distributions. Conditions for the existence of linear 2D stable processes are then studied in Section 3. Section 4 introduces a new model for stable 2elds which, in some sense, is the bidimensional equivalent of 1D FARIMA(0; d; 0) stable processes. In Section 5, we describe the synthesis algorithm and provide simulation examples to illustrate the eHectiveness of our approach. Some concluding remarks are given in Section 6.
2. Stable distributions The family of stable distributions [35] is interesting, since linear transforms preserve the distribution of any linear combination of independent, identically distributed -stable random variables. This property is directly related to the de2nition of stable random variables: a random variable X has a stable law if for all positive numbers A and B, there exists a positive number C and a real number D such that d
AX1 + BX2 = CX + D; where X1 and X2 are independent random variables with the same distribution as X . Another fundamental result is the generalized central limit theorem [6], which represents an alternative
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
de2nition of a stable law: a random variable X has a stable law if and only if there exists a sequence of independently, identically distributed random variables (Yn )n∈N and sequences of positive numbers (dn )n∈Z and real numbers (an )n∈Z such that Y1 + Y2 + · · · + Yn d + an ⇒ X; dn d
where the symbol ⇒ denotes the convergence in distribution. There does not exist explicit forms for most of the probability densities of stable variables. Exceptions are the Gaussian distribution ( = 2), the Cauchy distribution ( = 1) and the L%evy distribution ( = 12 ). Hence, an important way to characterize an -stable law is by means of its characteristic function which has the following form, in the symmetric case: E{eX } = exp(− || );
E|X |p ¡ ∞ as E|X |p = ∞
0 ¡ p ¡ ;
as p ¿ :
3. Linear stable random elds In the sequel, we will be concerned with 2D linear stable 2elds given by un; m =
∞ ∞ k=−∞ l=−∞
hk; l wn−k; m−l ;
where (wn; m )(n; m)∈Z2 are i.i.d. SS random variables with scale parameter w ¿ 0. As in2nite summations are involved in the above expression, the existence of un; m needs to be studied more precisely. Furthermore, some caution should be taken when de2ning the convergence of a 2D series. This makes the construction of a 2D linear random 2eld somehow more intricate than the construction of a 1D linear process. It can be shown [35] that the series in (1) converges absolutely almost surely if and only if ∞ ∞
(1)
|hk; l | ¡ ∞;
when 0 ¡ ¡ 1;
(2)
k=−∞ l=−∞ ∞ ∞ k=−∞ l=−∞
1 ¡ ∞; |hk; l | ln min |hk; l |w ; 2
(3)
when = 1, ∞ ∞
where ∈ (0; 2] and ¿ 0 are the two parameters of the symmetric -stable (SS) law. We note X ∼ S (; 0; 0). The parameters and of this distribution are called, respectively, the index of stability and the scale parameter. For a Gaussian random variable ( = 2); is proportional to the standard deviation, while the mean is zero. When ¿ 1, the scale parameter induces a norm X = on a vector space of jointly SS random variables. One of the major diGculties encountered when considering stable random variables is the impossibility of de2ning the variance and, when 6 1, even the mean. More precisely, if we have a random variable X ∼ S (; 0; 0), where 0 ¡ ¡ 2, then
1929
|hk; l | ¡ ∞;
when ¿ 1:
(4)
k=−∞ l=−∞
Furthermore, when ¿ 1, a necessary and suGcient condition for the series de2ning un; m to be convergent in the sense of the norm · (or equivalently, in the Lp sense if 1 6 p ¡ ) is ∞ ∞
|hk; l | ¡ ∞:
(5)
k=−∞ l=−∞
This condition is derived in Appendix A. Note that (5) reduces to the classical condition [4] (hk; l )(k; l)∈Z2 ∈ ‘2 (Z2 ) when (wn; m )(n; m)∈Z2 is a zero-mean white Gaussian noise. When one of the above existence conditions is satis2ed, the 2ltered random 2eld un; m is SS with scale parameter ∞ 1= ∞ |hk; l | : (6) u = w k=−∞ l=−∞
Three remarks can be made at this point: • When the existence conditions are satis2ed, the resulting process is (strictly) stationary as (1)
1930
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
obviously shows that its 2nite probability distributions are invariant under space translation. • According to (4), when almost sure absolute convergence is obtained, (hk; l )(k; l)∈Z2 ∈ ‘1 (Z2 ). As, when ¿ 1, ‘1 (Z2 ) ⊂ ‘ (Z2 ), we deduce from condition (5) that convergence also holds in Lp sense with 1 6 p ¡ . This result is consistent with the facts that almost sure convergence entails convergence in probability and that, in a space of jointly SS random variables, convergence in probability is equivalent to convergence w.r.t. the covariation norm · [35]. • Another possible approach to establish existence conditions is to consider (1) as an integral with respect to a singular SS measure and apply the results in [35].
4. 2D fractional -stable processes Our objective is now to de2ne a 2D discrete-space process having LRD properties. We recall that, in the 1D Gaussian case, such a process un is provided by an FARIMA(0; d; 0) with d ∈ R+ [12], which satis2es (1 − q−1 )d un = wn ; where wn is a zero-mean white Gaussian noise with variance 22d "2w and q±1 un = un±1 . As (1 − q−1 ) is a discrete-time derivation operator, (1 − q−1 )d represents a discrete-time fractional derivation of order d(∈ (0; 1=2)). The power spectrum density of the process is given by S(!) =
"2w |sin !=2|2d
and consequently, this random process can be viewed as the output of a linear 2lter with frequency response |sin(!=2)|−d driven by a white Gaussian noise. Note that the LRD behavior is caused by the divergence of the frequency response at ! = 0. By analogy, we de2ne a 2D fractional stable process un; m as the output of a bidimensional 2lter with frequency response Hd (!x ; !y ) =
1 (sin2 !x =2 + sin2 !y =2)d
(7)
whose input is an SS i.i.d. sequence wn; m with scale parameter w . When d ¡ 1, Hd belongs to L1 ([ − &; &]2 ) and it is proved in Appendix B that the impulse response hk; l of the 2lter is such that hk; l = O((k 2 + l2 )d−1 ) 2
(8)
2
when k + l → ∞. We deduce that, for a ¿ 0, |hk; l |a ¡ ∞ (k;l)∈Z2
if d ¡ 1 − 1=a. As the long dependence property is desired (d ¿ 0), conditions (2), (3) and (5) stated in the previous section show that the existence of un; m is guaranteed only if ¿ 1. In this case, the domain of validity of d is the interval (0; 1 − 1=). Note that this result is consistent with both those existing for 1D stable FARIMA processes and those recently established for 2D Gaussian processes with LRD [26,29]. 5. Experimental results 5.1. Synthesis method The method used to generate the proposed fractional model is a frequency-domain synthesis. The involved bidimensional linear 2ltering is based on a FFT implementation, using the form of the frequency response given by Eq. (7). From this point of view, this synthesis method is close to the method proposed in [14] for the synthesis of a 2D isotropic fBm, and it is also similar to the algorithm we proposed in [28] to generate 2D anisotropic long-range-dependent 2elds with 2nite variance. The corresponding Matlab code is given in Appendix C (function LRD2Dalpha). The -stable driving noise has been generated using the method described in [35, p. 41]. In Figs. 1 and 2, we present two 512 × 512 realizations of the proposed 2D fractional model in the Gaussian case ( = 2; d = 0:2 and 0.4). For the same values (d = 0:2 and 0.4) of the fractional parameter, we show in Figs. 3 and 4 two realizations of an -stable model with = 1:8, while in Figs. 5 and 6 two realization are displayed for = 1:2 and d = 0:1, resp., d = 0:15.
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
1931
Fig. 1. Left: a realization of the proposed fractional model in the Gaussian case for = 2, d = 0:2. Right: the same realization, shown in log scales.
Fig. 2. Left: a realization of the proposed fractional model in the Gaussian case for = 2, d = 0:4. Right: the same realization, shown in log-scales.
Note that, as decreases, the constraint on the choice of d becomes stronger: the upper bound on d is indeed 1 − 1=. According to this constraint, the value of the parameter d cannot be chosen greater than 0.1667 in the case = 1:2. We remark the clear diHerence between the Gaussian case and the in2nite variance situation. This is illustrated by the spikyness of the synthesized 2elds.
5.2. Post-processing The impulsiveness may render the visual inspection of the synthesized images diGcult. The characteristics of the synthesized 2eld become more visible if we apply a nonlinear post-processing. The most common nonlinearities one can apply in order to reduce the dynamics of the synthesized 2eld are a logarithmic transform, a hyperbolic tangent
1932
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
Fig. 3. Left: a realization of the proposed fractional model in the case = 1:8, d = 0:2. Right: the same realization, shown in log-scales.
Fig. 4. Left: a realization of the proposed fractional model in the case = 1:8, d = 0:4. Right: the same realization, shown in log-scales.
and a power law. The last operation is quite common in image processing, since it corresponds to a “gamma-correction” of the display. We 2rst compare these kinds of nonlinear operations on synthetic textures. The representation of the above process in log-scales is illustrated in Figs. 3– 6. For comparison, the Gaussian processes with the same LRD parameter are also represented in Figs. 1 and 2 in log-scales. Note that the log-scaled representation always corresponds to a second order
random 2eld. Indeed, if X is a SS random variable, log |X | is a second order random variable [23] with mean 1 E{log |X |} = CE − 1 + log X and variance &2 var (log |X |) = 6
1 1 + 2 2
;
where CE 0:57721566 is the Euler’s constant.
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
1933
Fig. 5. Left: a realization of the proposed fractional model in the case = 1:2, d = 0:1. Right: the same realization, shown in log-scales.
Fig. 6. Left: a realization of the proposed fractional model in the case = 1:2, d = 0:15. Right: the same realization, shown in log-scales.
Another kind of nonlinearity which will be used subsequently is expressed by the following transformation: ) |fn; m | gn; m = tanh |fn; m |* ; (9) ( where (gn; m )(n; m)∈Z2 denotes the resulting 2eld and (, ) and * are positive constants. Note that when * ¡ =2, gn; m is a 2eld with 2nite variance. The statistics of the data after this nonlinearity are not as well studied
as for the log-transform and three parameters have to be set for this transform. However, it provides much more Oexibility to 2t real data. Note that the kind of composite model we are using, which includes a dynamical linear part, followed by a nonlinear static part is usual in nonlinear control [7] and neural networks [10], and it can also be found in statistics [9]. This step may be important in practical applications as it can be expected that many natural textures actually correspond to
1934
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
5.3. Modeling of natural images
Fig. 7. The original SAR image of the Aix area.
2nite variance processes. Of course, the nonlinearity has to be adapted to the application at hand. This point is illustrated in the following two application examples.
We applied our composite model to simulate speckle noise in SAR images and the cosmic background in astronomical images. An ERS-1 SAR image of the region of Aix en Provence, France, is presented in Fig. 7. The original data are coded on 16 bits per pixel. The noise is roughly separated from the main features in the image using an 8 × 8 median 2lter. We compare this noise image with a synthetic SS 2eld with parameters = 1:95, d = 0:35, on which we applied the nonlinear operation (9) with ( = 500, ) = 12 and * = 14 . Fig. 8 shows a typical realization generated in this way. The second example corresponds to an astronomical image taken by the Hubble space telescope. Fig. 10 represents an active star-forming region within galaxy NGC 6822. Still using a basic median 2ltering technique, we separated the astronomical background, and compared it with a synthetic 2eld, generated by a SS 2eld with parameters =1:6, d=0:2, transformed with the same kind of nonlinearities (with ( = 50, ) = 0:6 and * = 0:35) as for the previous example. The two images are shown in Fig. 11. In order to illustrate the usefulness of the two assumptions (impulsiveness and LRD) on which our
Fig. 8. Left: original speckle noise in the SAR image in Fig. 7. Right: generated noise, using a composite model with the parameters = 1:95, d = 0:35 and the nonlinear processing (9).
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
1935
Fig. 9. Generated noise, using the same model as in Fig. 8, in which we drop the assumption of: impulsiveness—at left; LRD—at right.
6. Conclusions
Fig. 10. Original astronomical image of the NGC 6822 nebula.
model relies, we also provide for the two considered examples the 2elds synthesized using a Gaussian model with the same LRD parameter and a SS 2eld with the same stability parameter, but no correlation between samples. These textures are presented in Figs. 9 and 12.
In this paper, a new stochastic model for -stable 2elds with LRD properties has been proposed. Some of the mathematical properties of these processes have been studied. At the same time, a frequency-domain method has been described to generate realizations of this random 2eld. At this point, it is important to emphasize that, although other non-Gaussian distributions can be considered [31], the -stable model is especially interesting as it leads to a both simple and elegant characterization of the statistics of the resulting model. It is now possible to envisage applications of this model for texture characterization=classi2cation and noise modeling in SAR and astronomical imaging. However, the preliminary results shown in the paper would need to be further investigated. In particular, is would be useful to study the statistical estimators for the parameters and d of the proposed model. However, this is not an easy task. The statistical study of the estimators for the LRD parameters of Gaussian 2elds is still an on-going research topic [5]. In our case, the problem is even more critical due to diGculty of de2ning a spectral representation for the considered stable processes. The study of parameter estimation methods could be the topic of a future work. Some preliminary work may be found in [32].
1936
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
Fig. 11. Left: original astronomical background. Right: generated noise, using a composite model with the parameters = 1:6, d = 0:2 and the nonlinear processing (9).
Fig. 12. Generated noise, using the same model as in Fig. 11, in which we drop the assumption of: impulsiveness—at left; LRD—at right.
Acknowledgements The authors would like to thank Dr. E. Kuruoglu for his numerous comments which were helpful in improving this paper. Appendix A. Condition for the existence of linear 2D stable processes Let ¿ 1 and S() be a vector space of jointly SS random variables de2ned on a probability space (+; T; P). When 1 6 p ¡ , S() is a vector sub-
space of the Banach space Lp (+; T; P) and the norms · and · Lp are equivalent on S() [35] as there exists Cp; ∈ R∗+ such that, for all X ∈ S(), X = Cp; X Lp : Consequently, S() (endowed with the norm · ) is a Banach space if and only if it is a closed set in Lp (+; T; P). Coming back to the notations of Section 3, let us now de2ne L() = span {wk; l ; (k; l) ∈ Z2 }:
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
According to the previous remark, L() endowed with the norm · is a Banach space. We will show that Proposition 1. The series (k; l)∈Z2 hk; l wn−k; m−l converges in L() (or; equivalently; Lp (+; T; P); 1 6 p ¡ ) if and only if (5) is satis
hk; l wn−k; m−l ¡ j:
(k;l)∈J
As L() is a Banach space; a necessary and suGcient condition for the convergence of the series is provided by the Cauchy criterion: for all j ¿ 0; there exists a 2nite subset Ij of Z2 such that; for all 2nite set J ⊂ Z2 with J ∩ Ij = ∅;
hk; l wn−k; m−l ¡ j:
(k;l)∈J
By continuity of the norm, K K hk; l wn−k; m−l un; m = lim K→∞ k=−K l=−K
1=
|hk; l | w ¡ j:
(k;l)∈J
As we have assumed that w = 0; the previous inequality reduces to the Cauchy criterion for the series (k; l)∈Z2 |hk; l | and condition (5) is obtained.
∞ ∞
=
un; m = lim
K→∞
k=−K l=−K
Appendix B. Asymptotic behaviour of the impulse response of the 2D synthesis fractional lter Let 0 ¡ d ¡ 1 and consider & & e−ix·! f(x) = d!; 2 2 d −& −& (sin !x + sin !y ) where ! = (!x ; !y )T , x = (x; y)T and · is the usual scalar product: x · ! = x!x + y!y . Let ’ ∈ C ∞ (R2 ), such that ’(!)=0 when |!|∞ =max{|!x |; |!y |} ¿ & and ’(!) = 1 when |!|∞ 6 &=2. We can write f(x) = f1 (x) + f2 (x) + f3 (x); where
&
−&
f2 (x) =
&
−&
R2
e−ix·!
1 − ’(!) d!; (sin !x + sin2 !y )d 2
e−ix·! ’(!)
1 1 × − 2 2 2 d (!x + !y2 )d (sin !x + sin !y ) ’(!) f3 (x) = e−ix·! 2 d!: 2 d (! 2 R x + !y )
d!;
As the kernel of f1 is a C ∞ function, it follows that f1 decays rapidly. For example, & & 2 |x| f1 (x) = − e−ix·! 2 −&
×
hk; l wn−k; m−l :
1= |hk; l | w
and relation (6) follows.
Note that we have, in particular, K K
k=−∞ l=−∞
f1 (x) =
As J is 2nite and (wn−k; m−l )(k; l)∈J are i.i.d. SS random variables; the previous inequality is equivalent to
1937
−&
1 − ’(!) (sin2 !x + sin2 !y )d
d! = O(1);
1938
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
where 2 denotes the Laplacian operator. This implies f1 (x) = O(|x|−2 ) at least (actually, it has a faster decay). By using a 2rst order Taylor expansion, we have, for ! = 0: 3(!) =
=
=
1 1 − 2 2 d (! + !y2 )d (sin !x + sin !y ) x (!x2 + !y2 )d − (sin2 !x + sin2 !y )d 2
(sin !x + sin !y )d (!x2 + !y2 )d !x2 + !y2 − sin2 !x − sin2 !y
(sin2 !x + sin2 !y )d (!x2 + !y2 )d
d d−1 ;
where ∈ ]sin2 !x +sin2 !y ; !x2 +!y2 [. We deduce that 0 6 3(!) 6 d
!x2 + !y2 − sin2 !x − sin2 !y
(sin2 !x + sin2 !y )d (!x2 + !y2 )
:
This allows us to conclude that 3(!) = O(|!|2−2d );
as ! → 0
since !x2 −sin2 !x =O(!x4 ) and O(!x4 +!y4 )=O(|!|4 ). Using the same kind of arguments, it can be shown that the 2rst order partial derivatives of 3 are O(|!|1−2d ) and the second-order partial derivatives are O(|!|−2d ). Therefore, 3 and its 2rst and second-partial derivatives are locally integrable since 2d−1 ¡ 1. Similarly to f1 , it follows that f2 (x) = O(|x|−2 ). Finally, f3 (x) =
R2
e−ix·!
+
R2
f(x) = O(|x|−2 ) + C|x|−2+2d = O(|x|−2+2d ): Eq. (8) follows by noting that hk; l = f(k; l)=(2&)2 . Appendix C. Matlab codes
2
2
In conclusion,
’(!) − 1 d! (!x2 + !y2 )d
e−ix·! (!x2 + !y2 )−d d!:
It can be easily veri2ed that 2((’(!) − 1)=|!|2d ) ∈ L1 (R2 ). It follows that the Fourier transform of ((’(!) − 1)=|!|2d ) is O(|x|−2 ). The Fourier transform of |!|−2d is well known [37] and it is equal to C|x|−2+2d , where C is an explicitly known constant. Therefore, we have f3 (x) = O(|x|−2 ) + C|x|−2+2d .
function [s,S] = LRD2Dalpha(a, N, d) % [s,S] = LRD2Dalpha(a, N, d) : generation of a fractional % symmetric alpha stable field % a: stability index (1 ¡ a ¡ = 2) % (the driving noise has unit scale parameter) % N : horizontal and vertical field dimension % d : fractional parameter (0 ¡ d ¡ 1-1=a) % s : generated field % S : frequency response of fractional synthesis filter if (rem(N,2)~=0) error(’the dimension must be even’) end % Generation of stable driving noise w = srstabl1(a,N); W = fft2(w); % Creation of the spectrum density of the process, % for the half frequency plane ind = 0:N=2; ind2 = 0:N-1; [Ind1,Ind2] = meshgrid(ind,ind2); Ind1 = Ind1’; Ind2 = Ind2’; S(ind+1,ind2+1)=1.=((sin(pi* Ind1=N).ˆ 2 +sin(pi * Ind2=N). ˆ 2+eps).ˆ d); S(1,1) = 0; % Symmetrization of the spectrum density in order to generate % a real field S = symfft2(S); % Computation of the field in the spatial domain s = real(ifft2(S.* W)); % % % end of LRD2Dalpha % % %
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
function st = srstabl1(a,N) % st = srstabl1(a,N) : generates a field of iid symmetric alpha-stable % random variables with unit scale parameter % a : stability index % N : horizontal and vertical dimension of generated field gamma = (rand(N)-1=2)* pi; w = exprnd(1,N,N); % Use of the statistical toolbox st = sin(a * gamma).=(cos(gamma)). ˆ (1=a).*(cos((1-a)* gamma).=w). ˆ ((1-a)=a); % % % end of srstabl1 % % % function Ss = symfft2(S) % Ss = symfft2(S) : symmetrization of the 2D DFT field S(1:N=2+1,1:N) % of a real 2D NxN field % values of S(1,N=2+2:N) and S(N=2+1,N=2+2:N) may be arbitrary % Ss(1,1), Ss(1,N=2+1), Ss(N=2+1,1) and Ss(N=2+1,N=2+1) are % constrained to be real N = size(S,2); if size(S,1)~ = N=2+1 error(’Incorrect dimensions of original field’) end Ss = S; % Checking for the symmetries in the original field ind = N=2+1:N-1; Ss(1,1) = real(Ss(1,1)); Ss(1,N=2+1) = real(Ss(1,N=2+1)); Ss(1,ind+1) = conj(Ss(1,N-ind+1)); Ss(N=2+1,1) = real(Ss(N=2+1,1)); Ss(N=2+1,N=2+1) = real(Ss(N=2+1,N=2+1)); Ss(N=2+1,ind+1) = conj(Ss(N=2+1,N-ind+1)); % Symmetrization Ss(ind+1,1) = conj(S(N-ind+1,1)); ind2 = 1:N-1; Ss(ind+1,ind2+1) = conj(S(N-ind+1,N ind2+1)); % % % end of symfft2 % % %
1939
References [1] P. Abry, Ondelettes et Turbulences—Multir%esolutions, algorithmes de d%ecompositions, invariance d’%echelle et signaux de pression, Diderot, Editeurs des Sciences et des Arts, Paris, 1997. [2] P. Abry, D. Veitch, Wavelet analysis of long-range-dependent traGc, IEEE Trans. Inform. Theory 44 (1) (January 1998) 2–15. [3] A. Arneodo, F. Argoul, E. Bacry, J. Elezgaray, J.-F. Muzy, Ondelettes, multifractales et turbulence—de l’ADN aux croissances cristallines, Diderot, Editeurs des Sciences et des Arts, Paris, 1995. [4] P.J. Brockwell, R.A. Davis, Time Series: Theory and Methods, Springer, New York, 1995. [5] O. Capp%e, E. Moulines, J.-C. Pesquet, A. Petropulu, X. Yang, Long-range dependence and heavy tail modeling for teletraGc data, IEEE Signal Process. Mag., 2002, 19, 14–27. [6] W. Feller, An Introduction to Probability Theory and Its Applications, Vol. II, Wiley, New York, 1971. [7] R.A. Freeman, P.V. Kokotovic, Robust Nonlinear Control Design, Birkhauser, Boston, MA, 1996. [8] R. Ghozi, B.C. Levy, Critical Markov random 2elds and fractional Brownian motion in texture synthesis, Proceedings of IEEE International Conference on Image Processing, Vol. 3, Austin, TX, 1994, pp. 426 – 430. [9] W. Greblicki, M. Pawlak, Nonparametric identi2cation of Hammerstein systems, IEEE Trans. Inform. Theory 45 (1989) 409–418. [10] S. Haykin, Neural Networks: A Comprehensive Foundation, IEEE Press=Macmillan College Publishing Company, New York, 1994. [11] S. Hoefer, F. Heil, M. Pandit, R. Kumaresan, Segmentation of textures with diHerent roughness using the model of isotropic two-dimensional fractional Brownian motion, Proceedings of ICASSP 5 (1993) 53–56. [12] J.R.M. Hosking, Fractional diHerencing, Biometrika 68 (1) (1981) 165–176. [13] J. Ilow, H. Leung, Self-similar modeling using FARIMA processes with applications to satellite images, IEEE Trans. Image Process. 10 (5) (May 2001) 792–797. [14] L.M. Kaplan, C.-C.J. Kuo, An improved method for 2D self-similar image synthesis, IEEE Trans. Image Process 5 (5) (May 1996) 754–761. [15] R. Kashyap, P. Lapsa, Synthesis and estimation of random 2elds using long-correlation models, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-6 (November 1994) 800–809. [16] J. Keller, R. Crownover, R. Chen, Characteristics of natural scenes related to the fractal dimension, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-9 (September 1987) 621–627. [17] E.E. Kuruoglu, W.J. Fitzgerald, P.J.W. Rayner, Near optimal detection of signals in impulsive noise modelled with a symmetric alpha-stable distribution, IEEE Commun. Lett. 2 (10) (October 1998) 282–284.
1940
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
[18] E.E. Kuruoglu, J. Zerubia, Modelling Images Using alpha-stable textures, Physics in Signal and Image Processing Conference, January 2001. [19] M. Kutay, A. Petropulu, C. Piccoli, On modeling biomedical ultrasound RF echoes using a power-law shot-noise model, IEEE Trans. Ultrasonics Ferroelectrics Frequency Control 48 (4) (July 2001) 953–968. [20] S. L%eger, M. Pontier, Drap Brownien Fractionnaire, C. R. Acad. Sci. 329 (1999) 893–898. [21] B.B. Mandelbrot, The Fractal Geometry of Nature, W.H. Freeman and Company Eds, New York, 1983. [22] B.B. Mandelbrot, J.W. Van Ness, Fractional Brownian motions, fractional noises and applications, SIAM Rev. 10 (4) (1968) 422– 437. [23] C. Nikias, M. Shao, Signal Processing with Alpha-stable Distributions and Applications, Wiley, New York, 1995. [24] H.O. Peitgen, H. Jurgens, D. Saupe, Chaos and Fractals, Springer, New York, 1993. [25] O. Peitgen, D. Saupe (Eds.), The Science of Fractal Images, Springer, New York, 1988. [26] B. Pesquet-Popescu, Mod%elisation bidimensionnelle de processus non stationnaires et application aW l’%etude du fond sous-marin, Doctorat en Sciences de l’Ecole Normale Sup%erieure de Cachan, 1998. [27] B. Pesquet-Popescu, Wavelet packet analysis of 2D processes with stationary fractional increments, IEEE Trans. Inf. Theory 45 (3) (April 1999) 1033–1038. [28] B. Pesquet-Popescu, ModWeles fractionnaires bidimensionnels aW espace discret, Tech. Sci. Inform. 20 (9) (November 2001) 1173–1200. [29] B. Pesquet-Popescu, P. Larzabal, Analyse de textures aW l’aide de modWeles anisotropes aW longue d%ependance, Colloque GRETSI, Grenoble, France, 15 –19 September 1997, pp. 635 – 638.
[30] B. Pesquet-Popescu, P. Larzabal, Interpolation of nonstationary 2elds with stationary increments, Proceedings of ICASSP, Vol. IV, Seattle, USA, 1998, pp. 2197–2200. [31] B. Pesquet-Popescu, J.-C. Pesquet, Non-Gaussian anisotropic 2D models with long-range dependence, Proceedings of IEEE Workshop on Non-linear Signal and Image Processing, Antalya, 1999, pp. 867–871. [32] B. Pesquet-Popescu, J.-C. Pesquet, Bidimensional -stable models with long-range dependence, Proceedings of IEEE Workshop on Non-linear Signal and Image Processing, Antalya, 1999, pp. 199 –203. [33] O.V. Poliannikov, Y. Bao, H. Krim, Levy processes for image modelling, Proceedings of IEEE Signal Processing Workshop on Higher-Order Statistics, Caesarea, Isreal, June 14 –16, 1999. [34] I.S. Reed, P.C. Lee, T.K. Truong, Spectral representation of fractional brownian motion in n dimensions and its properties, IEEE Trans. Inform. Theory 41 (5) (September 1995) 1439– 1451. [35] G. Samorodnitsky, M.S. Taqqu, Stable Non-Gaussian Random Processes: Stochastic Models with In2nite Variance, Chapman and Hall, New York, 1994. [36] M. Shao, C.L. Nikias, Signal processing with fractional lower order moments: stable processes and their applications, Proc. IEEE 8 (7) (July 1993) 984–1009. [37] E.M. Stein, Singular Integrals and DiHerentiability Properties of Functions, Princeton University Press, Princeton, 1970. [38] C. Stewart, B. Moghaddam, K. Hintz, L. Novak, Fractional Brownian motion models for synthetic aperture radar imagery scene segmentation, Proc. IEEE 81 (10) (October 1993) 1511–1522.
Synthesis of bidimensional -stable models with long-range dependence B%eatrice Pesquet-Popescua; ∗ , Jean-Christophe Pesquetb a Dept.
Traitement du Signal et des Images, Ecole Nationale Superieure des Telecommunications, 46, rue Barrault, 75634 Paris Cedex 13, France b Institut Gaspard Monge and URA CNRS 820, Universit e de Marne-la-Vallee, 5, Boulevard Descartes, Champs sur Marne, 77454 Marne la Vallee Cedex 2, France Received 5 May 2001; received in revised form 9 February 2002
Abstract In the context of linear modeling, the main advantage of stable distributions is to allow the de2nition of non-Gaussian processes whose statistical properties are easy to characterize. In this work, we are interested in the design of a speci2c class of 2D discrete-space processes with stable distributions. A frequency domain method for the synthesis of these 2elds will be proposed which is similar to algorithms already used in the Gaussian case. The considered models represent 2D -stable extensions of the 1D fractional Gaussian noise. They exhibit long-range dependence properties and, consequently, they could provide interesting alternatives to image modeling techniques based on the 2D fractional Brownian motion. However, the conditions for the existence of such processes deserve special attention and they are derived in this paper. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Stable process; Long-range dependence; Bidimensional -stable models; Fractional stable noise; Texture modeling; Image synthesis
1. Introduction In the last years, signal processing using fractal models has represented an active research area stimulated by the plethora of applications (compression, computer graphics [16,25,24] geophysics [21], traf2c modeling in computer networks [2], turbulence phenomena [1,3], 2nance, etc.) where these concepts are potentially interesting. Numerous works on “fractals” have popularized the notion of “self-similarity”.
∗
Corresponding author. E-mail addresses: [email protected] (B. Pesquet-Popescu), [email protected] (J-C. Pesquet).
This concept is at the origin of relevant models for several natural phenomena. In the meantime, the non-stationary structure of self-similar processes put them at the core of the most recent preoccupations in signal and image processing. Stable processes have turned out to be good models for many impulsive signals and noises [17], when the probability distributions of the highly variable data have “heavy” tails. These distributions have in2nite variance and unde2ned higher-order moments, but it was pointed out in [36] that many signal processing algorithms based on second-order statistics can be transposed to fractional lower-order moments. In this paper, we present new approaches for generating fractal random 2elds and, by this way, we
0165-1684/02/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 2 ) 0 0 3 2 0 - 1
1928
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
generalize the results existing in the 1D case. The resulting textures can be further used in classi2cation or recognition tasks, as well as for numerical terrain models. In this later application, the fractal characteristic allows to realize more realistic processing. In the case of an interpolation for example, this leads to better results than those obtained by traditional methods relying on autoregressive or deterministic models [30]. The proposed models are interesting for modeling and synthesizing textures with impulsive and long-range dependence (LRD) behaviors. Related existing works concern mainly the Gaussian case, i.e., two-dimensional versions of the fractional Brownian motion (fBm) [8,11,14,15,20,26,27,34,38] and its discrete counterparts [28,31]. In works like [8,14,15] these models have been used for texture synthesis, while in [11,38], they have been exploited for segmentation of synthetic or satellite images (in the later case, the diGculty of computing texture metrics in high-speckle SAR imagery being one of the motivations for using fractal models). The strong relationship between LRD models and fractals (when the Hurst parameter H is such that 1 2 ¡ H ¡ 1) has been emphasized in several papers [12,22], and since then a lot of applications of these models have been considered, related to 1D and 2D phenomena. In [13], 2D isotropic FARIMA processes driven by one-sided exponential noise are used to model the sea clutter texture in SAR images collected by RADARSAT sensors. In [19], a 1=f -type spectral behavior is exhibited for ultrasound RF echoes, and the parameters of this model are related to tissues characterization features. The authors also point out the usefulness of a 3D model, but, for the sake of simplicity, they built a 1D model of the backscattered signal. Potential applications of 2D -stable models can be found in diHerent 2elds. It seems that three applications could be mainly concerned in image analysis, due to the impulsiveness of the underlying data: SAR imaging, ultrasound medical imaging and astronomical imaging. The impulsiveness of certain classes of images has been modeled by -stable 2elds [18,32] or by L%evy processes [33]. In this paper, we propose a 2D linear fractional process driven by a stable noise, that can be a useful model in the kind of applications cited above. In the context of linear modeling,
the main advantage of stable distributions is to allow us to de2ne non-Gaussian processes whose probability laws are easily deduced from those of their driving noise. In this work, we will be interested in the design of a speci2c class of 2D discrete-space processes with stable distributions. Furthermore, the considered models have LRD properties. Conditions for the existence of such long dependent -stable 2elds are provided, that lead to the de2nition of the validity domain for the fractional parameter. An FFT-based synthesis method is also proposed, which is close to other synthesis methods developed in the Gaussian case [14]. In order to illustrate the potential interest of the proposed model, we also provide preliminary results suggesting that they can be used to simulate speckle noise in SAR images and the cosmic background in astronomical images. The outline of this paper is as follows. In Section 2 we recall some useful facts on stable distributions. Conditions for the existence of linear 2D stable processes are then studied in Section 3. Section 4 introduces a new model for stable 2elds which, in some sense, is the bidimensional equivalent of 1D FARIMA(0; d; 0) stable processes. In Section 5, we describe the synthesis algorithm and provide simulation examples to illustrate the eHectiveness of our approach. Some concluding remarks are given in Section 6.
2. Stable distributions The family of stable distributions [35] is interesting, since linear transforms preserve the distribution of any linear combination of independent, identically distributed -stable random variables. This property is directly related to the de2nition of stable random variables: a random variable X has a stable law if for all positive numbers A and B, there exists a positive number C and a real number D such that d
AX1 + BX2 = CX + D; where X1 and X2 are independent random variables with the same distribution as X . Another fundamental result is the generalized central limit theorem [6], which represents an alternative
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
de2nition of a stable law: a random variable X has a stable law if and only if there exists a sequence of independently, identically distributed random variables (Yn )n∈N and sequences of positive numbers (dn )n∈Z and real numbers (an )n∈Z such that Y1 + Y2 + · · · + Yn d + an ⇒ X; dn d
where the symbol ⇒ denotes the convergence in distribution. There does not exist explicit forms for most of the probability densities of stable variables. Exceptions are the Gaussian distribution ( = 2), the Cauchy distribution ( = 1) and the L%evy distribution ( = 12 ). Hence, an important way to characterize an -stable law is by means of its characteristic function which has the following form, in the symmetric case: E{eX } = exp(− || );
E|X |p ¡ ∞ as E|X |p = ∞
0 ¡ p ¡ ;
as p ¿ :
3. Linear stable random elds In the sequel, we will be concerned with 2D linear stable 2elds given by un; m =
∞ ∞ k=−∞ l=−∞
hk; l wn−k; m−l ;
where (wn; m )(n; m)∈Z2 are i.i.d. SS random variables with scale parameter w ¿ 0. As in2nite summations are involved in the above expression, the existence of un; m needs to be studied more precisely. Furthermore, some caution should be taken when de2ning the convergence of a 2D series. This makes the construction of a 2D linear random 2eld somehow more intricate than the construction of a 1D linear process. It can be shown [35] that the series in (1) converges absolutely almost surely if and only if ∞ ∞
(1)
|hk; l | ¡ ∞;
when 0 ¡ ¡ 1;
(2)
k=−∞ l=−∞ ∞ ∞ k=−∞ l=−∞
1 ¡ ∞; |hk; l | ln min |hk; l |w ; 2
(3)
when = 1, ∞ ∞
where ∈ (0; 2] and ¿ 0 are the two parameters of the symmetric -stable (SS) law. We note X ∼ S (; 0; 0). The parameters and of this distribution are called, respectively, the index of stability and the scale parameter. For a Gaussian random variable ( = 2); is proportional to the standard deviation, while the mean is zero. When ¿ 1, the scale parameter induces a norm X = on a vector space of jointly SS random variables. One of the major diGculties encountered when considering stable random variables is the impossibility of de2ning the variance and, when 6 1, even the mean. More precisely, if we have a random variable X ∼ S (; 0; 0), where 0 ¡ ¡ 2, then
1929
|hk; l | ¡ ∞;
when ¿ 1:
(4)
k=−∞ l=−∞
Furthermore, when ¿ 1, a necessary and suGcient condition for the series de2ning un; m to be convergent in the sense of the norm · (or equivalently, in the Lp sense if 1 6 p ¡ ) is ∞ ∞
|hk; l | ¡ ∞:
(5)
k=−∞ l=−∞
This condition is derived in Appendix A. Note that (5) reduces to the classical condition [4] (hk; l )(k; l)∈Z2 ∈ ‘2 (Z2 ) when (wn; m )(n; m)∈Z2 is a zero-mean white Gaussian noise. When one of the above existence conditions is satis2ed, the 2ltered random 2eld un; m is SS with scale parameter ∞ 1= ∞ |hk; l | : (6) u = w k=−∞ l=−∞
Three remarks can be made at this point: • When the existence conditions are satis2ed, the resulting process is (strictly) stationary as (1)
1930
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
obviously shows that its 2nite probability distributions are invariant under space translation. • According to (4), when almost sure absolute convergence is obtained, (hk; l )(k; l)∈Z2 ∈ ‘1 (Z2 ). As, when ¿ 1, ‘1 (Z2 ) ⊂ ‘ (Z2 ), we deduce from condition (5) that convergence also holds in Lp sense with 1 6 p ¡ . This result is consistent with the facts that almost sure convergence entails convergence in probability and that, in a space of jointly SS random variables, convergence in probability is equivalent to convergence w.r.t. the covariation norm · [35]. • Another possible approach to establish existence conditions is to consider (1) as an integral with respect to a singular SS measure and apply the results in [35].
4. 2D fractional -stable processes Our objective is now to de2ne a 2D discrete-space process having LRD properties. We recall that, in the 1D Gaussian case, such a process un is provided by an FARIMA(0; d; 0) with d ∈ R+ [12], which satis2es (1 − q−1 )d un = wn ; where wn is a zero-mean white Gaussian noise with variance 22d "2w and q±1 un = un±1 . As (1 − q−1 ) is a discrete-time derivation operator, (1 − q−1 )d represents a discrete-time fractional derivation of order d(∈ (0; 1=2)). The power spectrum density of the process is given by S(!) =
"2w |sin !=2|2d
and consequently, this random process can be viewed as the output of a linear 2lter with frequency response |sin(!=2)|−d driven by a white Gaussian noise. Note that the LRD behavior is caused by the divergence of the frequency response at ! = 0. By analogy, we de2ne a 2D fractional stable process un; m as the output of a bidimensional 2lter with frequency response Hd (!x ; !y ) =
1 (sin2 !x =2 + sin2 !y =2)d
(7)
whose input is an SS i.i.d. sequence wn; m with scale parameter w . When d ¡ 1, Hd belongs to L1 ([ − &; &]2 ) and it is proved in Appendix B that the impulse response hk; l of the 2lter is such that hk; l = O((k 2 + l2 )d−1 ) 2
(8)
2
when k + l → ∞. We deduce that, for a ¿ 0, |hk; l |a ¡ ∞ (k;l)∈Z2
if d ¡ 1 − 1=a. As the long dependence property is desired (d ¿ 0), conditions (2), (3) and (5) stated in the previous section show that the existence of un; m is guaranteed only if ¿ 1. In this case, the domain of validity of d is the interval (0; 1 − 1=). Note that this result is consistent with both those existing for 1D stable FARIMA processes and those recently established for 2D Gaussian processes with LRD [26,29]. 5. Experimental results 5.1. Synthesis method The method used to generate the proposed fractional model is a frequency-domain synthesis. The involved bidimensional linear 2ltering is based on a FFT implementation, using the form of the frequency response given by Eq. (7). From this point of view, this synthesis method is close to the method proposed in [14] for the synthesis of a 2D isotropic fBm, and it is also similar to the algorithm we proposed in [28] to generate 2D anisotropic long-range-dependent 2elds with 2nite variance. The corresponding Matlab code is given in Appendix C (function LRD2Dalpha). The -stable driving noise has been generated using the method described in [35, p. 41]. In Figs. 1 and 2, we present two 512 × 512 realizations of the proposed 2D fractional model in the Gaussian case ( = 2; d = 0:2 and 0.4). For the same values (d = 0:2 and 0.4) of the fractional parameter, we show in Figs. 3 and 4 two realizations of an -stable model with = 1:8, while in Figs. 5 and 6 two realization are displayed for = 1:2 and d = 0:1, resp., d = 0:15.
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
1931
Fig. 1. Left: a realization of the proposed fractional model in the Gaussian case for = 2, d = 0:2. Right: the same realization, shown in log scales.
Fig. 2. Left: a realization of the proposed fractional model in the Gaussian case for = 2, d = 0:4. Right: the same realization, shown in log-scales.
Note that, as decreases, the constraint on the choice of d becomes stronger: the upper bound on d is indeed 1 − 1=. According to this constraint, the value of the parameter d cannot be chosen greater than 0.1667 in the case = 1:2. We remark the clear diHerence between the Gaussian case and the in2nite variance situation. This is illustrated by the spikyness of the synthesized 2elds.
5.2. Post-processing The impulsiveness may render the visual inspection of the synthesized images diGcult. The characteristics of the synthesized 2eld become more visible if we apply a nonlinear post-processing. The most common nonlinearities one can apply in order to reduce the dynamics of the synthesized 2eld are a logarithmic transform, a hyperbolic tangent
1932
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
Fig. 3. Left: a realization of the proposed fractional model in the case = 1:8, d = 0:2. Right: the same realization, shown in log-scales.
Fig. 4. Left: a realization of the proposed fractional model in the case = 1:8, d = 0:4. Right: the same realization, shown in log-scales.
and a power law. The last operation is quite common in image processing, since it corresponds to a “gamma-correction” of the display. We 2rst compare these kinds of nonlinear operations on synthetic textures. The representation of the above process in log-scales is illustrated in Figs. 3– 6. For comparison, the Gaussian processes with the same LRD parameter are also represented in Figs. 1 and 2 in log-scales. Note that the log-scaled representation always corresponds to a second order
random 2eld. Indeed, if X is a SS random variable, log |X | is a second order random variable [23] with mean 1 E{log |X |} = CE − 1 + log X and variance &2 var (log |X |) = 6
1 1 + 2 2
;
where CE 0:57721566 is the Euler’s constant.
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
1933
Fig. 5. Left: a realization of the proposed fractional model in the case = 1:2, d = 0:1. Right: the same realization, shown in log-scales.
Fig. 6. Left: a realization of the proposed fractional model in the case = 1:2, d = 0:15. Right: the same realization, shown in log-scales.
Another kind of nonlinearity which will be used subsequently is expressed by the following transformation: ) |fn; m | gn; m = tanh |fn; m |* ; (9) ( where (gn; m )(n; m)∈Z2 denotes the resulting 2eld and (, ) and * are positive constants. Note that when * ¡ =2, gn; m is a 2eld with 2nite variance. The statistics of the data after this nonlinearity are not as well studied
as for the log-transform and three parameters have to be set for this transform. However, it provides much more Oexibility to 2t real data. Note that the kind of composite model we are using, which includes a dynamical linear part, followed by a nonlinear static part is usual in nonlinear control [7] and neural networks [10], and it can also be found in statistics [9]. This step may be important in practical applications as it can be expected that many natural textures actually correspond to
1934
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
5.3. Modeling of natural images
Fig. 7. The original SAR image of the Aix area.
2nite variance processes. Of course, the nonlinearity has to be adapted to the application at hand. This point is illustrated in the following two application examples.
We applied our composite model to simulate speckle noise in SAR images and the cosmic background in astronomical images. An ERS-1 SAR image of the region of Aix en Provence, France, is presented in Fig. 7. The original data are coded on 16 bits per pixel. The noise is roughly separated from the main features in the image using an 8 × 8 median 2lter. We compare this noise image with a synthetic SS 2eld with parameters = 1:95, d = 0:35, on which we applied the nonlinear operation (9) with ( = 500, ) = 12 and * = 14 . Fig. 8 shows a typical realization generated in this way. The second example corresponds to an astronomical image taken by the Hubble space telescope. Fig. 10 represents an active star-forming region within galaxy NGC 6822. Still using a basic median 2ltering technique, we separated the astronomical background, and compared it with a synthetic 2eld, generated by a SS 2eld with parameters =1:6, d=0:2, transformed with the same kind of nonlinearities (with ( = 50, ) = 0:6 and * = 0:35) as for the previous example. The two images are shown in Fig. 11. In order to illustrate the usefulness of the two assumptions (impulsiveness and LRD) on which our
Fig. 8. Left: original speckle noise in the SAR image in Fig. 7. Right: generated noise, using a composite model with the parameters = 1:95, d = 0:35 and the nonlinear processing (9).
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
1935
Fig. 9. Generated noise, using the same model as in Fig. 8, in which we drop the assumption of: impulsiveness—at left; LRD—at right.
6. Conclusions
Fig. 10. Original astronomical image of the NGC 6822 nebula.
model relies, we also provide for the two considered examples the 2elds synthesized using a Gaussian model with the same LRD parameter and a SS 2eld with the same stability parameter, but no correlation between samples. These textures are presented in Figs. 9 and 12.
In this paper, a new stochastic model for -stable 2elds with LRD properties has been proposed. Some of the mathematical properties of these processes have been studied. At the same time, a frequency-domain method has been described to generate realizations of this random 2eld. At this point, it is important to emphasize that, although other non-Gaussian distributions can be considered [31], the -stable model is especially interesting as it leads to a both simple and elegant characterization of the statistics of the resulting model. It is now possible to envisage applications of this model for texture characterization=classi2cation and noise modeling in SAR and astronomical imaging. However, the preliminary results shown in the paper would need to be further investigated. In particular, is would be useful to study the statistical estimators for the parameters and d of the proposed model. However, this is not an easy task. The statistical study of the estimators for the LRD parameters of Gaussian 2elds is still an on-going research topic [5]. In our case, the problem is even more critical due to diGculty of de2ning a spectral representation for the considered stable processes. The study of parameter estimation methods could be the topic of a future work. Some preliminary work may be found in [32].
1936
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
Fig. 11. Left: original astronomical background. Right: generated noise, using a composite model with the parameters = 1:6, d = 0:2 and the nonlinear processing (9).
Fig. 12. Generated noise, using the same model as in Fig. 11, in which we drop the assumption of: impulsiveness—at left; LRD—at right.
Acknowledgements The authors would like to thank Dr. E. Kuruoglu for his numerous comments which were helpful in improving this paper. Appendix A. Condition for the existence of linear 2D stable processes Let ¿ 1 and S() be a vector space of jointly SS random variables de2ned on a probability space (+; T; P). When 1 6 p ¡ , S() is a vector sub-
space of the Banach space Lp (+; T; P) and the norms · and · Lp are equivalent on S() [35] as there exists Cp; ∈ R∗+ such that, for all X ∈ S(), X = Cp; X Lp : Consequently, S() (endowed with the norm · ) is a Banach space if and only if it is a closed set in Lp (+; T; P). Coming back to the notations of Section 3, let us now de2ne L() = span {wk; l ; (k; l) ∈ Z2 }:
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
According to the previous remark, L() endowed with the norm · is a Banach space. We will show that Proposition 1. The series (k; l)∈Z2 hk; l wn−k; m−l converges in L() (or; equivalently; Lp (+; T; P); 1 6 p ¡ ) if and only if (5) is satis
hk; l wn−k; m−l ¡ j:
(k;l)∈J
As L() is a Banach space; a necessary and suGcient condition for the convergence of the series is provided by the Cauchy criterion: for all j ¿ 0; there exists a 2nite subset Ij of Z2 such that; for all 2nite set J ⊂ Z2 with J ∩ Ij = ∅;
hk; l wn−k; m−l ¡ j:
(k;l)∈J
By continuity of the norm, K K hk; l wn−k; m−l un; m = lim K→∞ k=−K l=−K
1=
|hk; l | w ¡ j:
(k;l)∈J
As we have assumed that w = 0; the previous inequality reduces to the Cauchy criterion for the series (k; l)∈Z2 |hk; l | and condition (5) is obtained.
∞ ∞
=
un; m = lim
K→∞
k=−K l=−K
Appendix B. Asymptotic behaviour of the impulse response of the 2D synthesis fractional lter Let 0 ¡ d ¡ 1 and consider & & e−ix·! f(x) = d!; 2 2 d −& −& (sin !x + sin !y ) where ! = (!x ; !y )T , x = (x; y)T and · is the usual scalar product: x · ! = x!x + y!y . Let ’ ∈ C ∞ (R2 ), such that ’(!)=0 when |!|∞ =max{|!x |; |!y |} ¿ & and ’(!) = 1 when |!|∞ 6 &=2. We can write f(x) = f1 (x) + f2 (x) + f3 (x); where
&
−&
f2 (x) =
&
−&
R2
e−ix·!
1 − ’(!) d!; (sin !x + sin2 !y )d 2
e−ix·! ’(!)
1 1 × − 2 2 2 d (!x + !y2 )d (sin !x + sin !y ) ’(!) f3 (x) = e−ix·! 2 d!: 2 d (! 2 R x + !y )
d!;
As the kernel of f1 is a C ∞ function, it follows that f1 decays rapidly. For example, & & 2 |x| f1 (x) = − e−ix·! 2 −&
×
hk; l wn−k; m−l :
1= |hk; l | w
and relation (6) follows.
Note that we have, in particular, K K
k=−∞ l=−∞
f1 (x) =
As J is 2nite and (wn−k; m−l )(k; l)∈J are i.i.d. SS random variables; the previous inequality is equivalent to
1937
−&
1 − ’(!) (sin2 !x + sin2 !y )d
d! = O(1);
1938
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
where 2 denotes the Laplacian operator. This implies f1 (x) = O(|x|−2 ) at least (actually, it has a faster decay). By using a 2rst order Taylor expansion, we have, for ! = 0: 3(!) =
=
=
1 1 − 2 2 d (! + !y2 )d (sin !x + sin !y ) x (!x2 + !y2 )d − (sin2 !x + sin2 !y )d 2
(sin !x + sin !y )d (!x2 + !y2 )d !x2 + !y2 − sin2 !x − sin2 !y
(sin2 !x + sin2 !y )d (!x2 + !y2 )d
d d−1 ;
where ∈ ]sin2 !x +sin2 !y ; !x2 +!y2 [. We deduce that 0 6 3(!) 6 d
!x2 + !y2 − sin2 !x − sin2 !y
(sin2 !x + sin2 !y )d (!x2 + !y2 )
:
This allows us to conclude that 3(!) = O(|!|2−2d );
as ! → 0
since !x2 −sin2 !x =O(!x4 ) and O(!x4 +!y4 )=O(|!|4 ). Using the same kind of arguments, it can be shown that the 2rst order partial derivatives of 3 are O(|!|1−2d ) and the second-order partial derivatives are O(|!|−2d ). Therefore, 3 and its 2rst and second-partial derivatives are locally integrable since 2d−1 ¡ 1. Similarly to f1 , it follows that f2 (x) = O(|x|−2 ). Finally, f3 (x) =
R2
e−ix·!
+
R2
f(x) = O(|x|−2 ) + C|x|−2+2d = O(|x|−2+2d ): Eq. (8) follows by noting that hk; l = f(k; l)=(2&)2 . Appendix C. Matlab codes
2
2
In conclusion,
’(!) − 1 d! (!x2 + !y2 )d
e−ix·! (!x2 + !y2 )−d d!:
It can be easily veri2ed that 2((’(!) − 1)=|!|2d ) ∈ L1 (R2 ). It follows that the Fourier transform of ((’(!) − 1)=|!|2d ) is O(|x|−2 ). The Fourier transform of |!|−2d is well known [37] and it is equal to C|x|−2+2d , where C is an explicitly known constant. Therefore, we have f3 (x) = O(|x|−2 ) + C|x|−2+2d .
function [s,S] = LRD2Dalpha(a, N, d) % [s,S] = LRD2Dalpha(a, N, d) : generation of a fractional % symmetric alpha stable field % a: stability index (1 ¡ a ¡ = 2) % (the driving noise has unit scale parameter) % N : horizontal and vertical field dimension % d : fractional parameter (0 ¡ d ¡ 1-1=a) % s : generated field % S : frequency response of fractional synthesis filter if (rem(N,2)~=0) error(’the dimension must be even’) end % Generation of stable driving noise w = srstabl1(a,N); W = fft2(w); % Creation of the spectrum density of the process, % for the half frequency plane ind = 0:N=2; ind2 = 0:N-1; [Ind1,Ind2] = meshgrid(ind,ind2); Ind1 = Ind1’; Ind2 = Ind2’; S(ind+1,ind2+1)=1.=((sin(pi* Ind1=N).ˆ 2 +sin(pi * Ind2=N). ˆ 2+eps).ˆ d); S(1,1) = 0; % Symmetrization of the spectrum density in order to generate % a real field S = symfft2(S); % Computation of the field in the spatial domain s = real(ifft2(S.* W)); % % % end of LRD2Dalpha % % %
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
function st = srstabl1(a,N) % st = srstabl1(a,N) : generates a field of iid symmetric alpha-stable % random variables with unit scale parameter % a : stability index % N : horizontal and vertical dimension of generated field gamma = (rand(N)-1=2)* pi; w = exprnd(1,N,N); % Use of the statistical toolbox st = sin(a * gamma).=(cos(gamma)). ˆ (1=a).*(cos((1-a)* gamma).=w). ˆ ((1-a)=a); % % % end of srstabl1 % % % function Ss = symfft2(S) % Ss = symfft2(S) : symmetrization of the 2D DFT field S(1:N=2+1,1:N) % of a real 2D NxN field % values of S(1,N=2+2:N) and S(N=2+1,N=2+2:N) may be arbitrary % Ss(1,1), Ss(1,N=2+1), Ss(N=2+1,1) and Ss(N=2+1,N=2+1) are % constrained to be real N = size(S,2); if size(S,1)~ = N=2+1 error(’Incorrect dimensions of original field’) end Ss = S; % Checking for the symmetries in the original field ind = N=2+1:N-1; Ss(1,1) = real(Ss(1,1)); Ss(1,N=2+1) = real(Ss(1,N=2+1)); Ss(1,ind+1) = conj(Ss(1,N-ind+1)); Ss(N=2+1,1) = real(Ss(N=2+1,1)); Ss(N=2+1,N=2+1) = real(Ss(N=2+1,N=2+1)); Ss(N=2+1,ind+1) = conj(Ss(N=2+1,N-ind+1)); % Symmetrization Ss(ind+1,1) = conj(S(N-ind+1,1)); ind2 = 1:N-1; Ss(ind+1,ind2+1) = conj(S(N-ind+1,N ind2+1)); % % % end of symfft2 % % %
1939
References [1] P. Abry, Ondelettes et Turbulences—Multir%esolutions, algorithmes de d%ecompositions, invariance d’%echelle et signaux de pression, Diderot, Editeurs des Sciences et des Arts, Paris, 1997. [2] P. Abry, D. Veitch, Wavelet analysis of long-range-dependent traGc, IEEE Trans. Inform. Theory 44 (1) (January 1998) 2–15. [3] A. Arneodo, F. Argoul, E. Bacry, J. Elezgaray, J.-F. Muzy, Ondelettes, multifractales et turbulence—de l’ADN aux croissances cristallines, Diderot, Editeurs des Sciences et des Arts, Paris, 1995. [4] P.J. Brockwell, R.A. Davis, Time Series: Theory and Methods, Springer, New York, 1995. [5] O. Capp%e, E. Moulines, J.-C. Pesquet, A. Petropulu, X. Yang, Long-range dependence and heavy tail modeling for teletraGc data, IEEE Signal Process. Mag., 2002, 19, 14–27. [6] W. Feller, An Introduction to Probability Theory and Its Applications, Vol. II, Wiley, New York, 1971. [7] R.A. Freeman, P.V. Kokotovic, Robust Nonlinear Control Design, Birkhauser, Boston, MA, 1996. [8] R. Ghozi, B.C. Levy, Critical Markov random 2elds and fractional Brownian motion in texture synthesis, Proceedings of IEEE International Conference on Image Processing, Vol. 3, Austin, TX, 1994, pp. 426 – 430. [9] W. Greblicki, M. Pawlak, Nonparametric identi2cation of Hammerstein systems, IEEE Trans. Inform. Theory 45 (1989) 409–418. [10] S. Haykin, Neural Networks: A Comprehensive Foundation, IEEE Press=Macmillan College Publishing Company, New York, 1994. [11] S. Hoefer, F. Heil, M. Pandit, R. Kumaresan, Segmentation of textures with diHerent roughness using the model of isotropic two-dimensional fractional Brownian motion, Proceedings of ICASSP 5 (1993) 53–56. [12] J.R.M. Hosking, Fractional diHerencing, Biometrika 68 (1) (1981) 165–176. [13] J. Ilow, H. Leung, Self-similar modeling using FARIMA processes with applications to satellite images, IEEE Trans. Image Process. 10 (5) (May 2001) 792–797. [14] L.M. Kaplan, C.-C.J. Kuo, An improved method for 2D self-similar image synthesis, IEEE Trans. Image Process 5 (5) (May 1996) 754–761. [15] R. Kashyap, P. Lapsa, Synthesis and estimation of random 2elds using long-correlation models, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-6 (November 1994) 800–809. [16] J. Keller, R. Crownover, R. Chen, Characteristics of natural scenes related to the fractal dimension, IEEE Trans. Pattern Anal. Mach. Intell. PAMI-9 (September 1987) 621–627. [17] E.E. Kuruoglu, W.J. Fitzgerald, P.J.W. Rayner, Near optimal detection of signals in impulsive noise modelled with a symmetric alpha-stable distribution, IEEE Commun. Lett. 2 (10) (October 1998) 282–284.
1940
B. Pesquet-Popescu, J-C. Pesquet / Signal Processing 82 (2002) 1927 – 1940
[18] E.E. Kuruoglu, J. Zerubia, Modelling Images Using alpha-stable textures, Physics in Signal and Image Processing Conference, January 2001. [19] M. Kutay, A. Petropulu, C. Piccoli, On modeling biomedical ultrasound RF echoes using a power-law shot-noise model, IEEE Trans. Ultrasonics Ferroelectrics Frequency Control 48 (4) (July 2001) 953–968. [20] S. L%eger, M. Pontier, Drap Brownien Fractionnaire, C. R. Acad. Sci. 329 (1999) 893–898. [21] B.B. Mandelbrot, The Fractal Geometry of Nature, W.H. Freeman and Company Eds, New York, 1983. [22] B.B. Mandelbrot, J.W. Van Ness, Fractional Brownian motions, fractional noises and applications, SIAM Rev. 10 (4) (1968) 422– 437. [23] C. Nikias, M. Shao, Signal Processing with Alpha-stable Distributions and Applications, Wiley, New York, 1995. [24] H.O. Peitgen, H. Jurgens, D. Saupe, Chaos and Fractals, Springer, New York, 1993. [25] O. Peitgen, D. Saupe (Eds.), The Science of Fractal Images, Springer, New York, 1988. [26] B. Pesquet-Popescu, Mod%elisation bidimensionnelle de processus non stationnaires et application aW l’%etude du fond sous-marin, Doctorat en Sciences de l’Ecole Normale Sup%erieure de Cachan, 1998. [27] B. Pesquet-Popescu, Wavelet packet analysis of 2D processes with stationary fractional increments, IEEE Trans. Inf. Theory 45 (3) (April 1999) 1033–1038. [28] B. Pesquet-Popescu, ModWeles fractionnaires bidimensionnels aW espace discret, Tech. Sci. Inform. 20 (9) (November 2001) 1173–1200. [29] B. Pesquet-Popescu, P. Larzabal, Analyse de textures aW l’aide de modWeles anisotropes aW longue d%ependance, Colloque GRETSI, Grenoble, France, 15 –19 September 1997, pp. 635 – 638.
[30] B. Pesquet-Popescu, P. Larzabal, Interpolation of nonstationary 2elds with stationary increments, Proceedings of ICASSP, Vol. IV, Seattle, USA, 1998, pp. 2197–2200. [31] B. Pesquet-Popescu, J.-C. Pesquet, Non-Gaussian anisotropic 2D models with long-range dependence, Proceedings of IEEE Workshop on Non-linear Signal and Image Processing, Antalya, 1999, pp. 867–871. [32] B. Pesquet-Popescu, J.-C. Pesquet, Bidimensional -stable models with long-range dependence, Proceedings of IEEE Workshop on Non-linear Signal and Image Processing, Antalya, 1999, pp. 199 –203. [33] O.V. Poliannikov, Y. Bao, H. Krim, Levy processes for image modelling, Proceedings of IEEE Signal Processing Workshop on Higher-Order Statistics, Caesarea, Isreal, June 14 –16, 1999. [34] I.S. Reed, P.C. Lee, T.K. Truong, Spectral representation of fractional brownian motion in n dimensions and its properties, IEEE Trans. Inform. Theory 41 (5) (September 1995) 1439– 1451. [35] G. Samorodnitsky, M.S. Taqqu, Stable Non-Gaussian Random Processes: Stochastic Models with In2nite Variance, Chapman and Hall, New York, 1994. [36] M. Shao, C.L. Nikias, Signal processing with fractional lower order moments: stable processes and their applications, Proc. IEEE 8 (7) (July 1993) 984–1009. [37] E.M. Stein, Singular Integrals and DiHerentiability Properties of Functions, Princeton University Press, Princeton, 1970. [38] C. Stewart, B. Moghaddam, K. Hintz, L. Novak, Fractional Brownian motion models for synthetic aperture radar imagery scene segmentation, Proc. IEEE 81 (10) (October 1993) 1511–1522.