7, 2004
2
Numerical models of binary random elds on the basis of thresholds of Gaussian functions S.M. Prigarin, A. Martin, G. Winkler 519.6
. ., ., . -
// " . $. . %% / &'(. " . %- . | ( , 2004. | . 7, 2. | ". 165{175. " # "$ . % & & '" $ $ . & &
" , $ " " # , ") *
, " "
.
: , , , .
Prigarin S.M., Martin A., Winkler G. Numerical models of binary random .elds on the basis of thresholds of Gaussian functions // Siberian J. of Numer. Mathematics / Sib. Branch of Russ. Acad. of Sci. | Novosibirsk, 2004. | Vol. 7, 2. | P. 165{175. We present a few numerical models of binary stochastic +elds based on the thresholds of Gaussian functions and discuss the results of numerical experiments on estimating the models' parameters and simulation of the observed data. The considered models can be used, in particular, for texture analysis and synthesis, for simulation of stochastic structure of clouds in the atmosphere, as well as for solving other problems when statistical analysis and construction of binary random +elds are a part of research.
Key words: stochastic simulation, numerical models of random elds, binary elds, texture analysis and
synthesis.
1. Description of the models We study speci.c models of binary random .elds "(x), where x is a (multi-dimensional) parameter. Assume that the .elds values can attain only two numbers, "(x) 2 f0 1g, and the model can be described in the following way. Let u(x) be a real-valued Gaussian random .eld and A0 A1 = (;1 1) (A0 \ A1 = ) is a splitting of the real axes into two disjoint sets. Let us set ( u(x) 2 A0 "(x) = 01 for (1) for u(x) 2 A : 1
Supported by INTAS under Grant 01-0239, RFBR under Grant 03-05-64655, SB RAS under Grant 2003-2, and Grant \Leading Scienti+c School" HIII-1271.2003.1
166
S.M. Prigarin, A. Martin, G. Winkler
In this paper, we make two other assumptions to simplify the consideration. First, the Gaussian .eld u is supposed to be homogeneous with mean zero and the correlation function Ku(x) = Eu(x + y)u(x) and, second, the sets A0 , A1 are supposed to be of the form
A0 = (;1 d) A1 = 3d +1) d 2 (;1 +1) for Model A and
A0 = (;d d) A1 = (;1 d] 3d +1) d 2 (0 +1) for Model B. The parameter d will be called a threshold level and the models of binary .elds will be called threshold models. Remark 1. Evidently, more general threshold models can be considered, but at the .rst stage we restrict the study to these two models because of their simplicity and adaptability to a considerably large number of applied problems. In particular, these models occur in simulation of a stochastic structure of broken clouds (see, e.g., 31, 2]) and atmospheric precipitation .elds 33, 4]. In addition, we expect that the threshold models can be a promising tool in texture synthesis and analysis (cf. 35]). For a threshold model, the binary .eld " (just as the Gaussian .eld u) is a homogeneous one. Its average is m" = P (" = 1) = P (u 2 A1 ), which is equal to m" = 1 ; 6(d)
(2)
for Model A and to
m" = 1 ; erf(d) = 2(1 ; 6(d)) = 26(;d) (3) for Model B. Here 6 is the function of the standard normal distribution and erf() is the error function: erf(a) = 2(6(a) ; 0:5) a 0: For the covariance function K" (x) = E"(x + y)"(y) of the .eld "(x) we have K" (x) = R(Ku (x)) where and
R() =
Z f(): "()"()=1g
' ( ) d d
(4) (5)
q 2 2 ;1 ' ( ) = 2 1 ; 2 exp +2(1 ;;22)
is the probability density of a two-dimensional Gaussian random vector with zero mean, the unit variance of components and the correlation coe8cient between the components.
2. Distortion of correlations The covariance distortion properties for the general point-wise nonlinear transformations of the Gaussian functions were studied, for example, in 36, 7]. Here we concentrate our attention on relation (5) for Model A (cf. 33, 4, 6]) and Model B.
Numerical models of binary random elds on the basis of thresholds of Gaussian functions
167
First, consider Model A. In this case, formula (5) can be simpli.ed as follows + Z 1 +Z 1 R() = ' ( ) d d: d
d
For numerical calculations it is reasonable to turn to Owen's function (see, e.g., 38] and the references there): R() = 6(;d) ; 2T (d a) (5A) where s (6) a = 11 ; +
and
Z T (d a) = 21 e;d2 (1+u2 )=2 1 +duu2 0 a
is Owen's function that satis.es the following relations 38] Zh Zk '( ) d d = 21 6(h) + 12 6(k) ; T (h a1 ) ; T (k a2 ) ;1 ;1 for hk > 0,
Zh Zk ;1 ;1
for hk 0, where and
' ( ) d d = 21 6(h) + 12 6(k) ; T (h a1 ) ; T (k a2 ) ; 21 k ; h a1 = p h 1 ; 2
h ; k a2 = p k 1 ; 2
T (h a) = ;T (h ;a) T (h a) = T (;h a) T (h a) + T (ah 1=a) = 41 ; erf(h) erf(ah) (a 0 h 0) a) T (h 0) = 0 T (0 a) = arctg( 2 2 (h) T (h 1) = 1 ; 24erf(h) T (h 1) = 1 ; 4 erf (h 0): 8
As the average m" of the binary .eld " is equal to p1 = Pr(" = 1) = 1 ; 6(d) = 6(;d) and the variance of the .eld is equal to p1 p0 , where p0 = Pr(" = 0) = 6(d), then for the correlation and the normalized correlation functions we can write down E("(x + y) ; m")("(y) ; m") = 6(d)6(;d) ; 2T (d a) E("(x + y) ; m")("(y) ; m") = 1 ; 2T (d a) 2 ; min(p0 p1 ) 1 V" p0p1 max(p0 p1 ) where a is de.ned by (6) and = Ku (x) in (6). The functions above are odd with respect to .
168
S.M. Prigarin, A. Martin, G. Winkler
For Model B we have
R() = 2
Z Z
Z Z;d
+1 +1
+1
d
d ;1
d
' ( ) d d + 2
' ( ) d d
and using general relations for Owen's function one can obtain
R() = 46(;d) ; 4fT (d a) + T (d 1=a)g (5B) V" = 26(;d)(1 ; 26(;d)) d > 0: E("(x + y) ; m")("(y) ; m") = 46(d)6(;d) ; 4fT (d a) + T (d 1=a)g: q Here a = 11+; and = Ku (x). Note, that for Model B the function R is even, R(;) = R().
3. Numerical modeling of the Gaussian elds To numerically simulate a binary random .eld "(x) according to the threshold model (1) it is necessary to construct realizations of the homogeneous Gaussian .eld u(x). For that we used two well-known general methods: the spectral models and the moving average scheme (see, e.g., 36, 9]). As we performed numerical experiments only with random .elds on the plane, all the formulas below are given for the two-dimensional parameter x = (x1 x2 ).
3.1. The moving average scheme. The moving average scheme for a homogeneous Gaussian .eld with zero mean on the grid of integers can be described by formulas 36, 9, 10] 1 X 1 1 X 1 X X unm = Akj wn;km;j = wkj An;km;j (7) k=;1 j =;1
Z Z
k=;1 j =;1
a( )ei(k+ j) d d Akj = (21 )2 ; ; X X ;i(k+ j) q 2 a( ) = Akj e = (2) f ( ) k j X X 2 XX 1 f ( ) = (2)2 Knme;i(n+ m) = (21 )2 Akj e;i(k+ j ) n m k j Z Z Knm = f ( )ei(n+ m) d d :
(8) (9) (10)
; ;
Here wnm are orthogonal standard random variables (Ewnm = 0, Ejwnm j2 = 1, Ewnm w
f (; ; ) = f ( ) a(; ; ) = a( ) K;n;m = Knm A;k;j = Akj : For the numerical simulation of a moving average .eld the series in the formulas above should be truncated:
169
Numerical models of binary random elds on the basis of thresholds of Gaussian functions
unm =
N X M X k=;N j =;M
Akj wn;km;j :
(11)
Note, that for model (11) the correlations Knm are equal to zero if n > 2N or m > 2M .
3.2. Spectral models. A general theory of the spectral models of the Gaussian random .elds and additional references can be found in 39]. Below we present a brief information about the spectral models of isotropic homogeneous Gaussian .elds on the plane (exactly these models were used in the numerical experiment described below). A correlation function of any homogeneous isotropic .eld u(x1 x2 ) on the plane can be written down in the form Z1 2 Ku (x1 x2 ) = J0 (z (x2 + y2 )1=2 )G(dz ) (12) 0
where 2 is variance of the .eld, J0 is the Bessel function of the .rst kind, G(dz ) is a \radial" spectral measure on 30 1), G30 1) = 1. Further we assume that the measure G has a density: G(dz ) = g(z )dz . (Table of spectral densities and the corresponding correlation functions of isotropic .elds is presented in 39].) Assume that 0 = Z0 < Z1 < : : : < ZN ;1 < ZN = 1. For an approximate simulation of the isotropic homogeneous Gaussian .eld u(x1 x2 ) on the plane let us consider the following spectral model: N M X X uNM (x1 x2 ) = cn M ;1=2 (;2 ln nm)1=2 n=1
m=1
cos3x1 zn cos !nm + x2 zn sin !nm + 2nm ] where
c2n = G3Zn;1 Zn ) =
ZZn Zn;1
g(z ) dz
(13)
; nm ) !nm = (m M
zn are random variables distributed in 3Zn;1 Zn ) according to the density g(z )=c2n and nm , nm, nm are independent random variables uniformly distributed in 30 1]. Thus, the .eld u is approximated by the spectral model uNM which is the sum of N M harmonics with random amplitudes and random frequencies. The randomness of frequencies ensures the coincidence of correlation functions of the .eld u and its spectral model uNM (this technique is known as the spectrum randomization principle, see, e.g., 311] { one of the basic papers on spectral models). On the other hand, random frequencies make the distributions of the .eld uNM nonGaussian. But a sequence of the .elds uNM is asymptotically Gaussian if NM ! 1 and maxnN (c2n =M ) ! 0 (for convergence of spectral models see 39, 11]). The algorithm of the modeling is in the creation of the following arrays: A(n m) = cn (;2(ln nm )=M )1=2 D(n m) = 2nm B (n m) = zn cos !nm C (n m) = zn sin !nm and the value of the .eld at the desired point (x1 x2 ) is calculated by the formula N X M X uNM (x1 x2 ) = A(n m) cos3B (n m)x1 + C (n m)x2 + D(n m)]: n=1 m=1
170
S.M. Prigarin, A. Martin, G. Winkler
Di?erent modi.cations of the spectral model are used with the purpose to simplify the algorithm or to reproduce some properties of the .eld. In particular, the following changes in (12) (individually and in combination) are acceptable 1) zn = znm 2) nm = n 3) nm = where znm , n , are independent random variables with the corresponding distributions. Algorithmically simpler, but less @exible is Model (13), where cn = N ;1=2 and zn are independent and distributed on the whole semi-axis 30 1) with the density g(z ). The choice of a version of the simulation algorithm and its parameters (particularly, the number of harmonics) is speci.ed by the objective to represent more or less in detail the corresponding parts of the spectrum. It can strongly a?ect the form of realizations uNM (x1 x2 ).
4. Estimation of parameters and tuning of the model The parameters of the threshold model for a binary .eld "(x) are: the threshold level d and the correlation function Ku (x) of the homogeneous Gaussian .eld u(x). If a set of realizations of the binary .eld " is available, the parameters of the model can be found in the following way. First, we compute an estimation m" of the mean value m" for the .eld ". Then the approximation d for the threshold level can be found from equation m" = 1 ; 6(d ) for Model A and from the equation m" = 26(;d ) for Model B. Second, we compute an estimate K" (x) of the covariance function K" (x) and .nd an approximation Ku (x) for the correlation function Ku (x) of the Gaussian .eld u(x) on the basis of the equation K"(x) = R(Ku(x)) (14) where the function R is de.ned from (5), (5A), (5B).
Remark 2. The direct inversion of (14) can fail because the result Ku may not be posi-
tive de.nite and the inversion can be ill-de.ned or multi-valued (especially, in the case of Model B). Thus, to construct the estimate Ku except equation (14), that can be accepted approximately, it is reasonable to take into account additional considerations, like positive de.niteness, assumptions about a sign of the correlation function for the Gaussian .eld (for Model B), smoothness, the rate of convergence to zero, and etc. The next step in constructing a simulation algorithm is to set a numerical model of the Gaussian .eld u(x). First, it is necessary to choose a model to be most appropriate for the behavior of the obtained correlation function Ku . Then the parameters of the numerical model should be determined. For the moving average scheme (11) the parameters are values of orders N , M and the coe8cients Akj (the latter can be de.ned by formulas (8){(10)). The spectral models have a \richer" system of parameters as compared to the moving average scheme, and these models seem to be more suitable for the simulation of the .elds with longrange correlations. An additional operation here is estimating the spectrum. To calculate redial spectral density g(z ) of the .eld for the isotropic spectral models on the plane we used numerical approximation of the following representation + Z1 q g(z) = z rJ0 (zr)K (r) dr K ( x21 + x22 ) = Ku (x1 x2 ): 0
Numerical models of binary random elds on the basis of thresholds of Gaussian functions
171
Finally, we should mention an important stage of the construction of the threshold models of binary .elds to be performed at the very beginning on the basis of realization characteristics: selection of the type of the threshold model (Models A, B, the \negative" versions when zero and one swap places, or another model).
5. Results of numerical experiments The results of numerical experiments are presented in Figures 1{7. On the left, there is a realization of the observed binary .eld, and on the right { there is a realization of the threshold model constructed on the basis of the observed realization by the method described in the previous section. Some of the results can be interpreted like reasonable, but some of them are imperfect. The following conclusions were made during our studies: 1. Model B is more complicated and unstable in tuning than Model A. 2. Spectral models are more promising for random .elds with the long-distance correlations than the moving average scheme because the procedure to de.ne coe8cients for a high order moving average scheme is very time-consuming. 3. The proposed threshold models are not appropriate for simulation of the .elds with a regular structure (see Figures 3, 4, where the observed textures are borrowed from 312]). The explanation of the fact is the following. The Gaussian distribution has the highest entropy among all distributions with .xed moments of the .rst and the second orders. The threshold models inherit the property of high entropy, which is incompatible with a regular structure. All realizations on the right of Figures 1{7 were simulated by using spectral models of the Guassian .elds (with 100 80 harmonics) and threshold Model A (except Figure 3, where we tried to apply Model B). For simplicity, in the numerical experiment we assumed the observed .elds to be homogeneous and isotropic, while for some of the .elds this hypothesis evidently fails. Let us mention here that the spectral method and the moving average scheme can be applied to simulate homogeneous anisotropic Gaussian random .elds and even nonhomogeneous .elds as well (see 39]).
Figure 1. The observed binary eld simulated by the threshold method and a second order moving average scheme (the left image) and the simulated realization (the right image)
172
S.M. Prigarin, A. Martin, G. Winkler
Figure 2. The observed binary elds (on the left) (by K. Rodenacker), and the simulated realizations of the threshold model (on the right)
Numerical models of binary random elds on the basis of thresholds of Gaussian functions
173
Figure 3. The observed binary eld on the left is a binary version of a part of D75 from 12], and the simulated realization of the threshold model is on the right
Figure 4. The observed binary eld on the left is a part of D101 from 12], and the simulated realization of the threshold model is on the right
Figure 5. The observed binary eld on the left is a binary version of a part of D4 from 12], and the simulated realization of the threshold model is on the right
174
S.M. Prigarin, A. Martin, G. Winkler
Figure 6. The observed binary eld on the left is a binary version of a part of D100 from 12], and the simulated realization of the threshold model is on the right
Figure 7. The observed binary eld on the left is a binary version of a part of D40 from 12], and the simulated realization of the threshold model is on the right
Acknowledgements. We are grateful to Karsten Rodenacker who provided us with images that we used in our numerical experiments.
References
1] Kargin B.A., Prigarin S.M. Simulation of cumulus clouds for investigation of solar radiation transfer processes by Monte Carlo method // Optika Atm. i Okeana. | 1994. | Vol. 7, 9. | P. 1275{1287. 2] Prigarin S.M., Kargin B.A., Oppel U.G. Random elds of broken clouds and their associated direct solar radiation, scattered transmission, and albedo // Pure and Applied Optics, A. | 1998. | Vol. 7, . 6. | P. 1389{1402. 3] Ogorodnikov V.A. Testing of stochastic models for series of dry and rainy days // Proceedings ZapSibNIGMI Goskomgidrometa. | 1989. | Vol. 86. | P. 74{79 (in Russian). 4] Anisimova A.V. Numerical modeling of random indicator elds for rainy precipitation // Proceedings of Young Scientists Conference. | Novosibirsk: ICMMG, 1997. | P. 3{15 (in Russian).
Numerical models of binary random elds on the basis of thresholds of Gaussian functions
175
5] Winkler G. Image Analysis, Random Fields and Markov Chain Monte Carlo Methods. | Berlin: Springer, 2003. 6] Ogorodnikov V.A., Prigarin S.M. Numerical Modelling of Random Processes and Fields: Algorithms and Applications. | The Netherlands, Utrecht: VSP, 1996. 7] Prigarin S.M. Introduction to Numerical Modeling of Random Processes and Fields. | Novosibirsk: Univ. Press, 1999 (in Russian). 8] Smirnov N.V., Bolshev L.N. Tables for Calculating the Normal Distribution Functions / Izv. Akad. Nauk SSSR. | 1962 (in Russian). 9] Prigarin S.M. Spectral Models of Random Fields in Monte Carlo Methods. | The Netherlands, Utrecht: VSP, 2001. 10] Tovstik T.M. Simulation of homogeneous Gaussian eld // Proceedings of Symposium \Methods of Representation and Apparat. Anal. of Random Processes and Fields". | Leningrad, 1978. | P. 75{77 (in Russian). 11] Mikhailov G.A. Numerical construction of a random eld with given spectral density // Dokl. Akad. Nauk SSSR. | 1978. | Vol. 238, 4. | P. 793{795 (in Russian). 12] Brodatz P. Textures: a Photographic Album for Artists and Designers. | New York: Dover Publications, 1966.
Sergei M. Prigarin
Institute of Computational Mathematics and Mathematical Geophysics, Siberian Branch, Russian Academy of Sciences, Prosp. Akad. Lavrentyeva, 6, Novosibirsk, 630090, Russia E-mail:
[email protected]
Andreas Martin
Institute of Biomathematics and Biometry GSF-National Research Centre for Environment and Health 85764 Neuherberg/M"unchen, Germany E-mail:
[email protected]
Gerhard Winkler
Institute of Biomathematics and Biometry GSF-National Research Centre for Environment and Health 85764 Neuherberg/M"unchen, Germany E-mail:
[email protected]
The article submitted May 14, 2003 Revision submitted October 14, 2003
176
