- Email: [email protected]
Textural analysis of disordered materials with multifractals
Physica A 267 (1999) 221–238
Textural analysis of disordered materials with multifractals a Centre
Antoine Saucier a , Jiri Muller b; ∗
de Recherche en Calcul AppliquÃe, 5160 boul. DÃecarie, bureau 400, MontrÃeal (QuÃebec) Canada H3X 2H9 b Institutt for energiteknikk, Instituttveien 18, P.O. Box 40, N-2007 Kjeller, Norway Received 1 July 1998; received in revised form 13 November 1998
Abstract We propose a systematic method to choose the scaling range for multifractal analysis, and we illustrate this method by examining the texture of paper formation and pore space of sedimentary chalks. To clarify the statistical meaning of a texture characterization based on the moments of a measure (such as multifractal analysis), we examine the connection between these moments and c 1999 Elsevier Science B.V. All rights the multipoint correlations of the underlying signal. reserved. PACS: G1.43-j Keywords: Multifractals; Scaling range; Texture; Paper formation; Porous rocks
1. Introduction Textural analysis [1] can be of great signi cance to characterise disordered materials and to predict their physical properties. Such predictions can be made, e.g. by searching correlations between the physical properties of interest and textural parameters that are extracted from an image analysis of microstructure [2]. Textural analysis can also be useful to distinguish and classify di erent zones in a signal or an image. In the oil industry, for instance, one is regularly confronted with huge banks of images and=or signals (well logs, images from microresistivity well logs, etc.). The shear size of these data banks makes it impractical, if not impossible, to use only visual inspection for classi cation purposes. Automatic segmentation tools can be a valuable aid to analyse this type of data. In this context, texture analysis is a way to quantify what can be ∗
Corresponding author. Fax: +47-63-81-11-68; e-mail: [email protected].
c 1999 Elsevier Science B.V. All rights reserved. 0378-4371/99/$ – see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 9 8 ) 0 0 6 5 5 - 4
222
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
perceived visually. In this paper, we examine the problem of texture characterisation from the multifractal standpoint. This approach is appropriate when suitable properties are satis ed, which is the case for many objects in nature that exhibit fractal and multifractal behaviour over more than one order of magnitude in length (or time) scale [3]. The use of multifractal analysis raises several diculties in practice. In particular, the choice of an appropriate scaling range is a crucial step [4,5]. Indeed, this range must be selected before the analysis can be performed and its choice has usually a signi cant e ect on the fractal dimensions. Moreover, the choice of a scaling range is also important to justify the legitimacy and the relevance of the fractal approach. The rst diculty in selecting a suitable scaling range is to choose unambiguously, e.g. by visual inspection, a linear region on a log–log plot (i.e. the logarithm of the generating function versus the logarithm of a length scale). Other considerations that are speci c to the data analysed, such as a lower (or upper) bound on the scaling range, must also be taken into account in this choice. The range selected should exhibit a good linearity for all values of q, otherwise general multifractality constraints, such as 0 (q)¿0 and 00 (q)60, may not be satis ed (these constraints are actually consequences of the linearity for all q). Judging this linearity requires a correct examination of several logarithmic (log–log) plots simultaneously. Choosing among all these plots an optimal range is not a trivial task. For this reason, we believe that de ning better what is meant by “a good scaling range”, and proposing a tool to help locating such a range is a worthwhile goal in the context of multifractal analysis. In this paper we propose a de nition and a method to nd appropriate scaling ranges for multifractal analysis. The method proposed is helpful to analyse individual samples, and can be particularly useful to analyse long time series with a gliding window technique. The application of our method is presently illustrated by analysing the texture of paper formation images, and pore space of North Sea sedimentary chalk. However, in future, it will be extended to the analysis of microresistivity and image well logs, incorporating a sliding window processing.
2. Multifractal analysis 2.1. Implementation In this section we recall brie y basic notions of multifractal analysis and expose our implementation of this method, that includes uncertainty estimates that play a signi cant role for the choice of a scaling range. We use in this paper the multifractal formalism that was originally introduced to characterise measures [6]. Consider a 2D random positive signal S(ri ) de ned on a set of points ri = (xi ; yi ) located on a regular grid. The measure of a square box Br () of size centred on a point r is de ned by X S(ri ) : (2.1) pr () = ri ∈Br ()
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
223
If S(ri ) is interpreted to be a mass localised at point ri , then pr () is simply the total mass contained in the box Br (). We say that the measure has multifractal properties if the moments h[pr ()]q i, where q is a real parameter and angle brackets denote an ensemble average, take approximately the form [7] h[pr ()]q i = h[pr (1 )]q i(=1 )D+ (q)
(2.2a)
in a nite interval [min ; max ] called scaling range and min 61 6max . The exponents (q) are called mass exponents, and D is the dimension of the embedding space. Strictly speaking, the measure is said to be multifractal if (q) depends non-linearly on q. The form (2.2a) is valid if pr () ¿ 0 for all r and for min ¡ ¡ max , since 0q is not de ned for non-integer q. The more general form of (2.2a), valid for all q and all , is hq (pr ())i = hq (pr (1 ))i(=1 ) D+ (q) ;
(2.2b)
where q (x) ≡ xq if x ¿ 0 and q (0) = 0. The scaling properties of a multifractal measure are usually described in terms of its generalised dimensions D(q), which are related to mass exponents by D(q) = (q)=(q − 1). We stress that D(q) is not singular at q = 1 because (q) vanishes as q = 1; that is, (q) ∼ (q − 1)D(1) as q → 1. One of the simplest way to estimate generalised dimensions for a measure is the method of moments [7], using a generating function. In this paper, we always consider a collection of M samples of equal size that are derived from the same statistical ensemble. If a single sample is available, then we split it into several disjoint sub-samples of the same size. For instance, a single image will be split into four disjoint quadrants of the same size that will be regarded as four independent samples. Having several samples allow us to de ne precisely the uncertainties on the generating function. Firstly, we de ne for each sample an “individual” generating function q () by q () = −D hq (pr ())iS , where D is the dimension of the embedding space, and the angle brackets h: : :iS denote a spatial average. q () is therefore simply proportional to the moment hq (pr ())iS . Secondly, we de ne the “global” generating function ˆq () obtained from all M samples by ˆq () =
q () ; ( 1 ())q
(2.3a)
q () =
M 1 X (i) q () M
(2.3b)
where
i=1
and q(i) () denotes the individual generating function of sample i. De nition (2.3b) insures that ˆq () remains proportional to the same spatial average, 1 whereas the denominator of (2.3a) guarantees that ˆq () is normalised exactly, i.e. ˆ1 () = 1 for 1
Indeed, if hf(r)iS(Di ) denotes a spatial average of a function f(r) over a domain Di , then the spa-
tial average of f on a union M −1
PM
i=1
SM
i=1
Di of disjoint domains of the same size satis es hf(r)iS(∪M
hf(r)iS(Di ) . In the text, the function that plays the role of f is q (pr ()).
D) i=1 i
=
224
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
all , which is important to obtain accurate mass exponents. For multifractals, ˆq () satis es the power-law property ˆq () ∼ (q)
(2.4)
in the scaling range, which is a consequence of (2.2a) and (2.2b). The exponents (q) are obtained for each q by tting ˆq () to the model (2.4) using linear least-squares regressions of Zq = ln( ˆq ()) versus ln(). D(q) is then calculated from (q) according to D(q) = (q)=(q − 1). We derive the uncertainties used in the linear regressions from de nition (2.3a) as follows. Using the notations Zq ≡ ln( ˆq ()) and Xq = q (), the logarithm of (2.3a) takes the form Zq =ln(Xq )−q ln(X1 ). The uncertainty on Zq is de ned to be the standard deviation of the uctuation Zq = Zq − hZq i. Generally, the uctuation of a random variable R(V1 ; : : : ; VN ) which is a function of N variables V1 ; : : : ; VN is estimated by the formula R =
N X @R (hV1 i; : : : ; hVN i)Vi ; @Vi
(2.5)
i=1
where Vi ≡ Vi − hVi i. Using (2.5) on Zq , we obtain Zq = Xq =hXq i − qX1 =hX1 i. Assuming that the M samples are independent and belong to the same statistical ensemble, it follows from (2.3b) that hXq i = hXq i; 2 (Xq ) = 2 (Xq )=M and cov(Xq ; X1 ) = cov(Xq ; X1 )=M . Squaring Zq and averaging it nally yields the variance of Zq , ! 2 2 (Xq ) cov(Xq ; X1 ) 1 (X1 ) 2 ; (2.6) + q − 2q h(Zq ) i = M hXq i hX1 i hXq ihX1 i where cov(Xq ; X1 )=hXq X1 i−hXq ihX1 i. In this paper, we use formula (2.6) to calculate the uncertainties on generating functions. 2.2. The meaning of texture from a multifractal standpoint: connection between moments and multipoint correlations One of our motivations for using the moments h(pr ())q i to characterise the underlying signal S(r) lies in their ability to capture multipoint correlations. In this section, we will rst show that this ability is not restricted to pairs of points, but also includes triplets, quadruplets, etc., and more generally n-tuplets. In a second step, we will propose a comparison of integer and non-integer moments with respect to their dependence on multipoint correlations, and derive some interesting implications for multifractals. The simplest case is the moment of order 1 of pr (), i.e. the expectation value. Averaging (2.1) yields hpr ()i =
N X
hS(ri )i = N ;
(2.7)
i=1
where we assumed that hS(ri )i = for all i (translation invariance), and where the ri0 s of summation (2.1) have been indexed from i = 1; : : : ; N . We see from Eq. (2.7) that
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
225
hpr ()i is simply determined by the mean of the signal. Consider next h(pr ())2 i, which is the moment of order 2 of the measure. Squaring (2.1), then expanding and averaging, we obtain h(pr ())2 i =
N X N X
hS(ri )S(rj )i ;
(2.8)
i=1 j=1
h(pr ())2 i is therefore expressed in terms of the two-point correlation coecients hS(ri )S(rj )i, that involve all the pairs of points (ri ; rj ) which are contained in the box Br (). More generally, the moment of order n takes the form h(pr ())n i =
N N X X i1 =1 i2 =1
···
N X
hS(ri1 )S(ri2 ) · · · S(rin )i :
(2.9)
in =1
It is therefore determined by the n-point correlation coecients of the form hS(ri1 )S(ri2 ) · · · S(rin )i. The above expansions show that the integer moments of pr () can be used to characterise multipoint correlations in a signal. The order of the moment controls the number of points involved in the correlation coecients, and the box size controls the maximum distance separating two points. The main problem with the direct characterisation of the n-point correlation coecients is that they are functions of n variables, and are therefore dicult to represent. In contrast, the generating function of order q is expressed as a sum of correlation coecients of di erent orders and spatial locations ri (Eq. (2.9)), but depends only on one variable . Knowing the generating function for all scales and several q’s is an indirect way to characterise higher-order correlation functions. The main simpli cation gained from this approach is that for each q the generating function is a function of a single variable . Hence, the problem of characterising several functions of n variables (n¿1), i.e. the correlation coecients, has been reduced to the problem of characterising several functions of a single variable, i.e. the generating function for several q’s. This is our main motivation for using generating functions to describe texture. The understanding of non-integer moments is more complicated because simple expansions of type (2.9) cannot be used to reveal the role of multipoint correlations. However, we will show in the following that non-integer moments can be expressed in terms of integer moments. Consider a positive random variable P and its moment hP q i, where q is a real number. We can always write q in the form q = n + , where n is an integer and is a real number such that 06 ¡ 1. On one hand hP q i can obviously be written in the form hP q i = hP n P i :
(2.10)
On the other hand, we can expand the function P in Taylor series at the point P = hPi as long as hPi¿0. Truncating this expansion to the second order, and using the notation hP q i ≡ M [q], we obtain P ∼ = M [1] + M [1]−1 (P − M [1]) + 12 M [1]−2 ( − 1)(P − M [1])2 :
(2.11)
226
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Multiplying (2.11) by P n , then expanding and taking the expectation value yields (after simpli cation) M [n + ] ∼ = 12 (2 − 3 + 2 )M [1] M [n] + (2 − )M [1]−1 M [n + 1] 1 2 (
− 1)M [1]−2 M [n + 2] ;
(2.12)
M [n + ] is now expressed in terms of the integer moments M [1]; M [n]; M [n + 1] and M [n + 2]. A more accurate estimate of M [n + ] can be obtained by expanding P to higher orders (Eq. (2.11)), which will introduce an additional dependence of M [n + ] on M [n + 3]; M [n + 4], etc. Non-integer moments of order n + are therefore determined by integer moments of order m¿n, and consequently they depend on correlation coecients involving m points, where m¿n. Relationship (2.12) has interesting implications for multifractals. Indeed, we will show in the following that this relationship between moments of di erent orders leads to a corresponding relationship between dimensions. Consider multifractals for which the moments h(pr ())q i = M [q] are related to a scale ratio 0 ¡ 61 by the formula M [q] = 1+ (q) ;
(2.13)
i.e. we assume a unit prefactor in Eq. (2.2a), as well as D = 1. Substituting (2.13) into (2.12), and then using the property (1) = 0, we obtain 1 1+ (q) ∼ = −1 ((2 − 3 + 2 )2+ (n) − 2( − 2)1+ (n+1) + ( − 1) (n+2) ) : 2 (2.14) In Appendix A, we show that in the limit → 1 expansion (2.14) leads to a relationship between (n + ) and the mass exponents { (n); (n + 1); (n + 2)}: (n ˆ + ) ∼ = 12 (2 − 3 + 2 ) (n) + (2 − ) (n + 1) + 12 ( − 1) (n + 2) :
(2.15)
It is stressed that (n+) ˆ is an estimate (i.e. approximation) of the exact value (n+). It is seen from Eq. (2.15) that (n+) ˆ is a linear combination of (n); (n+1); (n+2), where the coecients of each (m) are functions of only. Higher-order expansions of (2.11) would introduce a dependence of (n ˆ + ) on (n + 3); (n + 4), etc. Since mass exponents and generalised dimensions are linked by (q) = (q − 1)D(q), we can ˆ also obtain from (2.15) an estimate D(q) of D(q) ˆ + ) = D(n
1 [1=2(2 − 3 + 2 )(n − 1)D(n) + (2 − )nD(n + 1) n+−1 +1=2( − 1)(n + 1)D(n + 2)] :
(2.16)
The special case n = 1 and then = 0 yields in particular ˆ D(1) = 2D(2) − D(3) :
(2.17)
It is stressed that setting = 0 in Eq. (2.16) yields the trivial identity D(n) = D(n) ˆ for all n except for n = 1. Repeating this whole derivation of D(1) with higher-order
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
227
expansions (starting from Eq. (2.11)) yields successively the following higher-order estimates: Order 2 :
ˆ D(1) = 2D(2) − D(3) ;
Order 3 :
ˆ D(1) = 3D(2) − 3D(3) + D(4) ;
Order 4 :
ˆ D(1) = 4D(2) − 6D(3) + 4D(4) − D(5) ;
Order 5 :
ˆ D(1) = 5D(2) − 10D(3) + 10D(4) − 5D(5) + D(6) :
(2.18)
It is amusing to notice that the coecients of the dimensions in (2.18) form a Pascal triangle. The accuracy of these formulas can be examined with the binomial measure for which the entropy dimension is given exactly by D(1) = −(w1 log2 (w1 ) + w2 log2 (w2 )). Using, for instance, the weights w1 = 0:4 and w2 = 0:6, we obtain D(1) = 0:970951, whereas the approximations (2.18) yield successively the estimates {0:968582; 0:971306; 0:971045; 0:970934} that approach the correct value. Formulas (2.18) show that the information dimension D(1) depends on the other integer dimensions, and therefore depends on multipoint correlations. In the above example, the formula of order 3 gives an error smaller than 0.04% which implies that for this multifractal, D(1) is essentially determined by correlation coecients involving 2; 3, and 4 points.
3. Optimal choice of a scaling range The dimensions D(q) are extracted from the generating function with linear regressions done in the scaling range. However, this range is not known a priori, and it must be chosen before the analysis. The choice of a scaling range should be done carefully because it can signi cantly a ect D(q). Moreover, an inadequate choice can lead to an illegitimate function D(q), i.e. to a function that violates some basic properties of multifractals, such as D0 (q)60 for instance. In this section, we propose for the choice of a scaling range some guidelines on which an automatic selection procedure can be based. On one hand, it is desirable to use a wide scaling range to exploit the capacity of multifractal analysis to characterise multipoint statistics. On the other hand, the width of the scaling range is limited by the goodness of t. Indeed, it is normally easier to obtain a good linear t with a narrow range than with a broad range (a good linear t of the generating function can always be achieved at the price of shrinking the scaling range suciently, as long as the error bars do not vanish. Indeed, the generating function is smooth (di erentiable) and is therefore linear locally. From this standpoint, the broadest intervals giving an acceptable t are interesting candidates for the scaling range. Unfortunately, the goodness of t is not always a sucient condition to assess linearity in a given interval.
228
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Indeed, large error bars can result in an acceptable t even if the generating function is far from linear. In this situation, the dimension function D(q) derived from the regressions may not be legitimate in the sense that it could violate some characteristic properties of multifractals, such as the condition D0 (q)60 for instance. To avoid such unsatisfactory scaling ranges, it is necessary to check that some basic properties of mass exponents are indeed satis ed. For this purpose, we propose to use the well-known properties [7] 0 (q)¿0 ;
00 (q)60
(3.1)
as constraints on multifractality. We regard these conditions as minimal requirements of multifractality. In particular, the property D0 (q)60 can be derived from these conditions. The multifractality conditions (3.1) are direct consequences of linearity. Indeed, if the power-law property h(p(n ))q i = nD+ (q)
(3.2)
holds for a range of values of q and n , then taking the logarithm yields ln(h(p(n ))q i) = (D + (q))ln(n ) :
(3.3)
Taking the derivative of (3.3) with respect to q, we obtain h(p(n ))q ln(p(n ))i = 0 (q)ln(n ) : h(p(n ))q i
(3.4)
Because n ¡ 1 and 0¡p(n )¡1, by de nition, it follows from Eq. (3.4) that 0 (q)¿0. Taking another derivative of (3.4) with respect to q, it can be shown that 00 (q)60 is also satis ed. Eq. (3.3) expresses the linearity of the plot ln(h(p(n ))q i) versus ln(n ) for each q, and we have shown that inequalities (3.1) are the direct consequences of this linearity. Conversely, forcing the scaling range to respect these inequalities is an indirect way to obtain linearity. Experience has shown us that a linear regression performed on a function F(q; x) that is not quite close to the ane form F(q; x) = (q)x + c(q), e.g. a function with constant curvature, will typically not give a regression coecient (q) that satis es both the multifractality conditions (3.1). This scaling range selection criterion becomes more stringent when a broader range of q’s is used. Combined with the goodness of t, conditions (3.1) applied to a broad range of q’s typically select intervals that are quite linear. Following the above arguments, we de ne a fractal range to be an interval that gives an acceptable linear t to the generating function (i.e. a reduced chi-square smaller than 1) and that satis es the multifractality conditions (3.1). In general, there will be several fractal ranges of a given size n, and the one that gives the best t (i.e. the smallest value of the reduced chi-square) will be called the best-ÿt fractal range of size n. Recalling that wide scaling ranges (i.e. large n) yield a better characterisation, we nally arrive at the following de nition of an “optimal” scaling range, that embodies all the above argumentation: The optimal scaling range is the widest and best-ÿt fractal range. The automatic selection procedure that corresponds to this de nition is the following. For each interval width n, we consider all subintervals and do for each
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
229
one the linear regression ln( ˆq ()) versus ln(), get the reduced chi-square r2 and the mass exponents (q). Next, we compute numerically 0 (q) and 00 (q). Among the intervals that satisfy the multifractality conditions (3.1), if they exist, we select the one that has the smallest r2 . If this minimum r2 is smaller than 1, then we nd the best- t fractal range of size n. This procedure is repeated for each n, from the smallest (n = 3 consecutive points) to the largest. The widest interval obtained with this procedure is the optimal scaling range chosen for the nal analysis. In the following section, we will use this procedure to locate the optimal scaling range within a given interval. It is stressed that the range obtained depends on the interval chosen for the search. For instance, searching the interval [1; 100] could give the optimal scaling range [50; 100], for instance. In a second step, we can search the interval [1; 50] and obtain another optimal – but narrower – scaling range, e.g. [10; 20]. This iterative procedure allows to nd several valid scaling ranges. It should be stressed that in this paper we adopt a rather pragmatic attitude towards scaling: an interval quali es as a legitimate range if a good t can be obtained with model (3.3), and if inequalities (3.1) are satis ed for all q’s in this range (otherwise the function (q) does not make any sense in the context of multifractals). The procedure may therefore select rather narrow intervals as acceptable scaling ranges. We do not make any statement about the actual nature of the experimental data tted, i.e it could di er from a power law. All we say is that the intervals selected with our method allow to obtain a good t and produce a legitimate multifractal spectrum.
4. First example: analysis of paper formation Paper formation is a technical term that refers to the variations of optical density in a sheet of paper. These variations are easily observed by viewing a paper sheet with a light source behind it. It is the disordered spatial distribution of paper bres that produces the optical density variations and thus e ects the texture of the paper. In the paper industry, it is common practice to digitise the light intensity transmitted through a sheet of paper, and then to analyse this data (formation analysis). A good knowledge of paper formation is of practical signi cance because it can be related to paper’s physical properties such as strength and printability. In this section, we show an example where multifractal analysis can be appropriate for the characterisation of paper formation, and could therefore be a useful complement to the more conventional methods based on Fourier analysis [8] that are used in the paper industry. We have used three paper formation images obtained from three di erent types of papers which were given the symbolic names “sum1”, “ort1” and “kaj1”. The images were digitised with a resolution of 512 × 512, for a physical size of 6:4 cm × 6:4 cm. These paper formations are shown in Fig. 1. It can be seen that the illumination is inhomogeneous, i.e. it tends to be stronger at the centre than on the edges of the image. We partly compensated for this inhomogeneous illumination by subtracting from the image a smoothed version of the same image. The smoothed image was obtained by
230
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Fig. 1. Paper formations images. Top: “kaj1”; middle: “ort1”; bottom: “sum1”.
tting a polynomial of order 4 to each line of the image (Fig. 2a and b). In addition to these formation samples, we constructed an arti cial sample for which the light intensity g(x; y) was sampled for each pixel from a uniform random number generator ranging from 0 and 255.
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
231
Fig. 2. (a) Top: Plot showing the light intensity (irregular curve) together with the fourth-order polynomial tted through it (smooth dotted curve). The polynomial background is subtracted from the signal to correct for the inhomogeneous illumination. (b) bottom: On the left is the original formation “kaj1”, and on the right its transform corrected for the inhomogeneous illumination.
To focus on the variability of the formation image, we applied a transformation to the original images. Let g(x; y) denote the transmitted light intensity at point (x; y) (partly compensated for the inhomogeneous illumination). A transformed signal S(x; y) was obtained by computing the squared norm of the intensity gradient, i.e. S(x; y) = (@g=@x)2 + (@g=@y)2 , and S(x; y) was used to de ne the measure according to (2:1). A typical view of this transformed signal is shown in Fig. 3. The singular appearance of S(x; y) is striking (i.e. spikes having a broad range of amplitudes), which makes this image qualitatively suitable for multifractal analysis, which is a singularity analysis. The relevance of multifractal analysis as a characterisation method is linked to the degree of validity of the power-law property (2:2) (or (2:4) for the generating function).
232
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Fig. 3. A portion of the transformed signal S(x; y) obtained from the formation kaj1. The singular appearance of this signal is striking.
More precisely, the multifractal approach is most appropriate when the scaling range [min ; max ] is broad, thus allowing to characterise variability over a wide range of length scales with a one-variable dimension function D(q). We show in Fig. 4 the generating functions obtained from the three paper samples, plotted for q = 2. These three functions are similar (they look virtually superposed) and exhibit a power-law behaviour (i.e. linear region in log–log coordinates) that starts at the largest scale and extends down to about 10 pixels, which corresponds to about 1.5 decade. These linear segments are the optimal scaling ranges that were obtained with the selection procedure described in Section 3. For the three types of paper formations analysed, we can say that the multifractal model (2:2) is valid for a broad range of length scales. We display in Fig. 5 the generalised dimensions D(q) obtained from the optimal scaling ranges for the three paper samples, with 06q62: For this arti cial sample we obtained D(q) ≈ 2 for each q, which is the signature of a (statistically) uniform image, i.e. an image for which measure and area are proportional (i.e. “monofractal”). The rst obvious observation is that all samples have high values of D(1), i.e. values which are close to 2, thus indicating a rather homogeneous spatial distribution of the measure. Nevertheless, the di erences between the D(q)’s and their uncertainties (error bars) allow to distinguish the di erent formations unambiguously. In our future work we wish to correlate some features of the D(q) curve (such as D(1)) with the physical properties of the paper.
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
233
Fig. 4. Generating functions obtained from three di erent paper formations, displayed for q = 2. The optimal scaling ranges of each sample are marked by squares, circles and crosses. The linear region starts at the largest scales and extends down to about 10 pixels, which corresponds to approximately 1.5 decade.
Fig. 5. Generalised dimensions obtained from each of the paper samples, together with the dimensions obtained from an arti cial random (statistically uniform) sample.
5. Second example: analysis of pore space in sedimentary chalk The rocks analysed are chalk cores extracted from a geological formation in the Valhall eld in the North sea (oil reservoir). We analysed six core plugs labelled 1,3,4,7,9 and 10. The core plugs were used to prepare thin sections suitable for the
234
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Fig. 6. Samples of the thin sections of the rocks analysed.
analysis of the geometry of pore space with scanning electron microscope (SEM). The SEM images of the thin sections were digitised with a resolution of 512 × 512, for a physical size of about 100 m. We obtained four SEM images for each plug. These images were next thresholded in such a way that pore space would be given a black colour, while the rock matrix received a white colour. More details on this data acquisition procedure can be found in [9]. Samples of these images are shown in Fig. 6.
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
235
Fig. 7. Typical generating functions for the chalk rocks. The error bars are small for all samples. The dots mark the optimal scaling ranges. No acceptable scaling range was found at small scales, i.e. between 1 and 10 pixels. Top: The optimal scaling range obtained from the complete interval 1–512 is shown as “range1”, and corresponds to the homogeneous limit, i.e. D(q) ≈ 2 for all q. The optimal scaling range obtained from the interval 1– 40 is shown as “range2”. Bottom: Because of the tiny error bars, we obtained two optimal scaling ranges at large scales (“range1” and “range2”) that both correspond to rather homogeneous behaviour, i.e. D(q) ≈ 2. A third optimal scaling range (“range3”), rather narrow, was found by searching the interval 1–90.
For such geometrical multifractals [4,5], we de ne the measure of a square to be the area of the pore space contained in this square. In other words, the function S(r) used in (2:1) is simply the indicator function of the pore space, i.e. S(r) = 1 if r is the pore space, and S(r) = 0 otherwise. The four images corresponding to each core plug were used to calculate a single generating function according to formula (2:3b). Typical generating functions obtained from these rocks are shown in Fig. 7. We can summarise
236
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Fig. 8. Generalised dimensions obtained for each of the rock samples.
their behaviour as follows. Searching an optimal scaling range in the whole interval, i.e. 1–512 pixels, we obtained for all rocks a scaling range located in the large scales (e.g. 40 –512) that correspond to the homogeneous limit, 2 i.e. such that D(q) ≈ 2 for all q (monofractal behaviour). Searching the remaining interval (in this case 1– 40) for an optimal scaling range yields another range, narrower than the rst one, e.g. 10 – 40. For all rock samples, it was found that there was no acceptable scaling range between 1 and 10 pixels (the multifractality conditions (3:1) are never satis ed in this interval, even if a good t can be obtained), which are length scales of the order of the typical pore size (about 10 pixels). For several rock samples, we found three disjoint scaling ranges before reaching the length scale of 10 pixels (Fig. 7). The informative behaviour of the generating function, i.e. away from the large-scale homogeneous limit, occurs typically between 1 and 80 pixels, which corresponds to about two decades. The acceptable scaling ranges lie somewhere between 10 and 80 pixels, i.e. above the typical pore size (about 10 pixels) and below the homogeneous limit (about 80 pixels), and the width of these ranges reaches at most three quarters of a decade, i.e. about half of the informative range [1,8]. The D(q)’s obtained from the scaling ranges contained in the interval 10 – 40 were plotted together in Fig. 8. We notice a distinct multifractal behaviour in sedimentary rocks which we also observed in our earlier work [9,10]. Such behaviour can be attributed to the ow properties of these rocks [9], and this is a subject of our future investigations.
2
For a large box size L, the measure px (L) of a box centered on a point x is simply proportional to the box area and to the macroscopic porosity of the medium, and does not vary with x.
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
237
6. Conclusions We proposed guidelines for the choice of the scaling range for multifractal analysis. For paper formations, we discovered a convincing multifractal behaviour, i.e. a power law that extends over most of the range of scales available. Multifractal analysis could therefore be a valuable tool for the characterisation of paper formations. The analysis of the porous rock images showed that the multifractal approach allows to characterise the images variability over approximately half of the informative range of scales. Appendix A Eq. (2.14) (main text) is an expansion of the form M [q] = l+ (q) ∼ =
N X
a i i :
(A.1)
i=1
Using L’Hospital’s rule, we can show that the asymptotic property N X
ai i ∼ !
(A.2)
i=1
as → 1 holds for any N , where the exponent ! is given by !, , N N N X X X i ai ai i ai : ln() = ! = lim ln →1
i=1
i=1
(A.3)
i=1
PN Moreover, it can be shown that i=1 ai = 1 and consequently (A.1) gives for all q an exact estimate of M [q] for = 1, i.e. M [q] = 1. Comparing (A.1) and (A.2) leads us to de ne an estimate (q) ˆ of the mass exponent (q) by ! ≡ 1 + (q). ˆ Using result PN (A.3), we obtain (q) ˆ = −1 + i=1 ai i , which takes the explicit form given by (2:15) in the main text. References [1] R. Haralick et al., Textural features for image classi cation, IEEE Trans. Systems Man Cybernet. SMC 3 (6) (1973) 610 – 621. [2] J.C. Russ, The Image Processing Handbook, Boca Raton, CRC Press, 1995. [3] D. Avnir, D. Lidar, O. Malcai, On the abundance of fractals, in: M.M. Novak, T.G. Dewey (Eds.), Fractals in the Natural and Applied Sciences, World Scienti c, Singapore, 1997. [4] A. Saucier, J. Muller, Characterization of porous media with geometrical multifractals, Fractals (4) (1993) 894 –903. [5] A. Saucier, J. Muller, Remarks on some properties of geometrical multifractals, Physica A 199 (1993) 350–362. [6] T.C. Halsey, M.H. Jensen, L.P. Kadano , I. Procaccia, B.I. Shraiman, Fractal measures and their singularities; The characterization of strange sets, Phys. Rev. A 33 (1986) 1141–1151. [7] C. Meneveau, K.R. Sreenivasan, The multifractal nature of turbulent energy dissipation, J. Fluid. Mech. 224 (1991) 429– 484.
238
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
[8] J. Muller, G. Helgesen, IFE Report IFE=KR=F-90=121. [9] J. Muller, Characterization of the North Sea chalk by multifractal analysis, J. Geophys. Res. 99 (B4) (1994) 7275 –7280. [10] A. Saucier, O.K. Huseby, J. Muller, Electrical texture characterization of dipmeter microresistivity signals using multifractal analysis, J. Geophys. Res. 102 (B5) (1997) 10 327–10 337.
Textural analysis of disordered materials with multifractals a Centre
Antoine Saucier a , Jiri Muller b; ∗
de Recherche en Calcul AppliquÃe, 5160 boul. DÃecarie, bureau 400, MontrÃeal (QuÃebec) Canada H3X 2H9 b Institutt for energiteknikk, Instituttveien 18, P.O. Box 40, N-2007 Kjeller, Norway Received 1 July 1998; received in revised form 13 November 1998
Abstract We propose a systematic method to choose the scaling range for multifractal analysis, and we illustrate this method by examining the texture of paper formation and pore space of sedimentary chalks. To clarify the statistical meaning of a texture characterization based on the moments of a measure (such as multifractal analysis), we examine the connection between these moments and c 1999 Elsevier Science B.V. All rights the multipoint correlations of the underlying signal. reserved. PACS: G1.43-j Keywords: Multifractals; Scaling range; Texture; Paper formation; Porous rocks
1. Introduction Textural analysis [1] can be of great signi cance to characterise disordered materials and to predict their physical properties. Such predictions can be made, e.g. by searching correlations between the physical properties of interest and textural parameters that are extracted from an image analysis of microstructure [2]. Textural analysis can also be useful to distinguish and classify di erent zones in a signal or an image. In the oil industry, for instance, one is regularly confronted with huge banks of images and=or signals (well logs, images from microresistivity well logs, etc.). The shear size of these data banks makes it impractical, if not impossible, to use only visual inspection for classi cation purposes. Automatic segmentation tools can be a valuable aid to analyse this type of data. In this context, texture analysis is a way to quantify what can be ∗
Corresponding author. Fax: +47-63-81-11-68; e-mail: [email protected].
c 1999 Elsevier Science B.V. All rights reserved. 0378-4371/99/$ – see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 9 8 ) 0 0 6 5 5 - 4
222
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
perceived visually. In this paper, we examine the problem of texture characterisation from the multifractal standpoint. This approach is appropriate when suitable properties are satis ed, which is the case for many objects in nature that exhibit fractal and multifractal behaviour over more than one order of magnitude in length (or time) scale [3]. The use of multifractal analysis raises several diculties in practice. In particular, the choice of an appropriate scaling range is a crucial step [4,5]. Indeed, this range must be selected before the analysis can be performed and its choice has usually a signi cant e ect on the fractal dimensions. Moreover, the choice of a scaling range is also important to justify the legitimacy and the relevance of the fractal approach. The rst diculty in selecting a suitable scaling range is to choose unambiguously, e.g. by visual inspection, a linear region on a log–log plot (i.e. the logarithm of the generating function versus the logarithm of a length scale). Other considerations that are speci c to the data analysed, such as a lower (or upper) bound on the scaling range, must also be taken into account in this choice. The range selected should exhibit a good linearity for all values of q, otherwise general multifractality constraints, such as 0 (q)¿0 and 00 (q)60, may not be satis ed (these constraints are actually consequences of the linearity for all q). Judging this linearity requires a correct examination of several logarithmic (log–log) plots simultaneously. Choosing among all these plots an optimal range is not a trivial task. For this reason, we believe that de ning better what is meant by “a good scaling range”, and proposing a tool to help locating such a range is a worthwhile goal in the context of multifractal analysis. In this paper we propose a de nition and a method to nd appropriate scaling ranges for multifractal analysis. The method proposed is helpful to analyse individual samples, and can be particularly useful to analyse long time series with a gliding window technique. The application of our method is presently illustrated by analysing the texture of paper formation images, and pore space of North Sea sedimentary chalk. However, in future, it will be extended to the analysis of microresistivity and image well logs, incorporating a sliding window processing.
2. Multifractal analysis 2.1. Implementation In this section we recall brie y basic notions of multifractal analysis and expose our implementation of this method, that includes uncertainty estimates that play a signi cant role for the choice of a scaling range. We use in this paper the multifractal formalism that was originally introduced to characterise measures [6]. Consider a 2D random positive signal S(ri ) de ned on a set of points ri = (xi ; yi ) located on a regular grid. The measure of a square box Br () of size centred on a point r is de ned by X S(ri ) : (2.1) pr () = ri ∈Br ()
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
223
If S(ri ) is interpreted to be a mass localised at point ri , then pr () is simply the total mass contained in the box Br (). We say that the measure has multifractal properties if the moments h[pr ()]q i, where q is a real parameter and angle brackets denote an ensemble average, take approximately the form [7] h[pr ()]q i = h[pr (1 )]q i(=1 )D+ (q)
(2.2a)
in a nite interval [min ; max ] called scaling range and min 61 6max . The exponents (q) are called mass exponents, and D is the dimension of the embedding space. Strictly speaking, the measure is said to be multifractal if (q) depends non-linearly on q. The form (2.2a) is valid if pr () ¿ 0 for all r and for min ¡ ¡ max , since 0q is not de ned for non-integer q. The more general form of (2.2a), valid for all q and all , is hq (pr ())i = hq (pr (1 ))i(=1 ) D+ (q) ;
(2.2b)
where q (x) ≡ xq if x ¿ 0 and q (0) = 0. The scaling properties of a multifractal measure are usually described in terms of its generalised dimensions D(q), which are related to mass exponents by D(q) = (q)=(q − 1). We stress that D(q) is not singular at q = 1 because (q) vanishes as q = 1; that is, (q) ∼ (q − 1)D(1) as q → 1. One of the simplest way to estimate generalised dimensions for a measure is the method of moments [7], using a generating function. In this paper, we always consider a collection of M samples of equal size that are derived from the same statistical ensemble. If a single sample is available, then we split it into several disjoint sub-samples of the same size. For instance, a single image will be split into four disjoint quadrants of the same size that will be regarded as four independent samples. Having several samples allow us to de ne precisely the uncertainties on the generating function. Firstly, we de ne for each sample an “individual” generating function q () by q () = −D hq (pr ())iS , where D is the dimension of the embedding space, and the angle brackets h: : :iS denote a spatial average. q () is therefore simply proportional to the moment hq (pr ())iS . Secondly, we de ne the “global” generating function ˆq () obtained from all M samples by ˆq () =
q () ; ( 1 ())q
(2.3a)
q () =
M 1 X (i) q () M
(2.3b)
where
i=1
and q(i) () denotes the individual generating function of sample i. De nition (2.3b) insures that ˆq () remains proportional to the same spatial average, 1 whereas the denominator of (2.3a) guarantees that ˆq () is normalised exactly, i.e. ˆ1 () = 1 for 1
Indeed, if hf(r)iS(Di ) denotes a spatial average of a function f(r) over a domain Di , then the spa-
tial average of f on a union M −1
PM
i=1
SM
i=1
Di of disjoint domains of the same size satis es hf(r)iS(∪M
hf(r)iS(Di ) . In the text, the function that plays the role of f is q (pr ()).
D) i=1 i
=
224
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
all , which is important to obtain accurate mass exponents. For multifractals, ˆq () satis es the power-law property ˆq () ∼ (q)
(2.4)
in the scaling range, which is a consequence of (2.2a) and (2.2b). The exponents (q) are obtained for each q by tting ˆq () to the model (2.4) using linear least-squares regressions of Zq = ln( ˆq ()) versus ln(). D(q) is then calculated from (q) according to D(q) = (q)=(q − 1). We derive the uncertainties used in the linear regressions from de nition (2.3a) as follows. Using the notations Zq ≡ ln( ˆq ()) and Xq = q (), the logarithm of (2.3a) takes the form Zq =ln(Xq )−q ln(X1 ). The uncertainty on Zq is de ned to be the standard deviation of the uctuation Zq = Zq − hZq i. Generally, the uctuation of a random variable R(V1 ; : : : ; VN ) which is a function of N variables V1 ; : : : ; VN is estimated by the formula R =
N X @R (hV1 i; : : : ; hVN i)Vi ; @Vi
(2.5)
i=1
where Vi ≡ Vi − hVi i. Using (2.5) on Zq , we obtain Zq = Xq =hXq i − qX1 =hX1 i. Assuming that the M samples are independent and belong to the same statistical ensemble, it follows from (2.3b) that hXq i = hXq i; 2 (Xq ) = 2 (Xq )=M and cov(Xq ; X1 ) = cov(Xq ; X1 )=M . Squaring Zq and averaging it nally yields the variance of Zq , ! 2 2 (Xq ) cov(Xq ; X1 ) 1 (X1 ) 2 ; (2.6) + q − 2q h(Zq ) i = M hXq i hX1 i hXq ihX1 i where cov(Xq ; X1 )=hXq X1 i−hXq ihX1 i. In this paper, we use formula (2.6) to calculate the uncertainties on generating functions. 2.2. The meaning of texture from a multifractal standpoint: connection between moments and multipoint correlations One of our motivations for using the moments h(pr ())q i to characterise the underlying signal S(r) lies in their ability to capture multipoint correlations. In this section, we will rst show that this ability is not restricted to pairs of points, but also includes triplets, quadruplets, etc., and more generally n-tuplets. In a second step, we will propose a comparison of integer and non-integer moments with respect to their dependence on multipoint correlations, and derive some interesting implications for multifractals. The simplest case is the moment of order 1 of pr (), i.e. the expectation value. Averaging (2.1) yields hpr ()i =
N X
hS(ri )i = N ;
(2.7)
i=1
where we assumed that hS(ri )i = for all i (translation invariance), and where the ri0 s of summation (2.1) have been indexed from i = 1; : : : ; N . We see from Eq. (2.7) that
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
225
hpr ()i is simply determined by the mean of the signal. Consider next h(pr ())2 i, which is the moment of order 2 of the measure. Squaring (2.1), then expanding and averaging, we obtain h(pr ())2 i =
N X N X
hS(ri )S(rj )i ;
(2.8)
i=1 j=1
h(pr ())2 i is therefore expressed in terms of the two-point correlation coecients hS(ri )S(rj )i, that involve all the pairs of points (ri ; rj ) which are contained in the box Br (). More generally, the moment of order n takes the form h(pr ())n i =
N N X X i1 =1 i2 =1
···
N X
hS(ri1 )S(ri2 ) · · · S(rin )i :
(2.9)
in =1
It is therefore determined by the n-point correlation coecients of the form hS(ri1 )S(ri2 ) · · · S(rin )i. The above expansions show that the integer moments of pr () can be used to characterise multipoint correlations in a signal. The order of the moment controls the number of points involved in the correlation coecients, and the box size controls the maximum distance separating two points. The main problem with the direct characterisation of the n-point correlation coecients is that they are functions of n variables, and are therefore dicult to represent. In contrast, the generating function of order q is expressed as a sum of correlation coecients of di erent orders and spatial locations ri (Eq. (2.9)), but depends only on one variable . Knowing the generating function for all scales and several q’s is an indirect way to characterise higher-order correlation functions. The main simpli cation gained from this approach is that for each q the generating function is a function of a single variable . Hence, the problem of characterising several functions of n variables (n¿1), i.e. the correlation coecients, has been reduced to the problem of characterising several functions of a single variable, i.e. the generating function for several q’s. This is our main motivation for using generating functions to describe texture. The understanding of non-integer moments is more complicated because simple expansions of type (2.9) cannot be used to reveal the role of multipoint correlations. However, we will show in the following that non-integer moments can be expressed in terms of integer moments. Consider a positive random variable P and its moment hP q i, where q is a real number. We can always write q in the form q = n + , where n is an integer and is a real number such that 06 ¡ 1. On one hand hP q i can obviously be written in the form hP q i = hP n P i :
(2.10)
On the other hand, we can expand the function P in Taylor series at the point P = hPi as long as hPi¿0. Truncating this expansion to the second order, and using the notation hP q i ≡ M [q], we obtain P ∼ = M [1] + M [1]−1 (P − M [1]) + 12 M [1]−2 ( − 1)(P − M [1])2 :
(2.11)
226
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Multiplying (2.11) by P n , then expanding and taking the expectation value yields (after simpli cation) M [n + ] ∼ = 12 (2 − 3 + 2 )M [1] M [n] + (2 − )M [1]−1 M [n + 1] 1 2 (
− 1)M [1]−2 M [n + 2] ;
(2.12)
M [n + ] is now expressed in terms of the integer moments M [1]; M [n]; M [n + 1] and M [n + 2]. A more accurate estimate of M [n + ] can be obtained by expanding P to higher orders (Eq. (2.11)), which will introduce an additional dependence of M [n + ] on M [n + 3]; M [n + 4], etc. Non-integer moments of order n + are therefore determined by integer moments of order m¿n, and consequently they depend on correlation coecients involving m points, where m¿n. Relationship (2.12) has interesting implications for multifractals. Indeed, we will show in the following that this relationship between moments of di erent orders leads to a corresponding relationship between dimensions. Consider multifractals for which the moments h(pr ())q i = M [q] are related to a scale ratio 0 ¡ 61 by the formula M [q] = 1+ (q) ;
(2.13)
i.e. we assume a unit prefactor in Eq. (2.2a), as well as D = 1. Substituting (2.13) into (2.12), and then using the property (1) = 0, we obtain 1 1+ (q) ∼ = −1 ((2 − 3 + 2 )2+ (n) − 2( − 2)1+ (n+1) + ( − 1) (n+2) ) : 2 (2.14) In Appendix A, we show that in the limit → 1 expansion (2.14) leads to a relationship between (n + ) and the mass exponents { (n); (n + 1); (n + 2)}: (n ˆ + ) ∼ = 12 (2 − 3 + 2 ) (n) + (2 − ) (n + 1) + 12 ( − 1) (n + 2) :
(2.15)
It is stressed that (n+) ˆ is an estimate (i.e. approximation) of the exact value (n+). It is seen from Eq. (2.15) that (n+) ˆ is a linear combination of (n); (n+1); (n+2), where the coecients of each (m) are functions of only. Higher-order expansions of (2.11) would introduce a dependence of (n ˆ + ) on (n + 3); (n + 4), etc. Since mass exponents and generalised dimensions are linked by (q) = (q − 1)D(q), we can ˆ also obtain from (2.15) an estimate D(q) of D(q) ˆ + ) = D(n
1 [1=2(2 − 3 + 2 )(n − 1)D(n) + (2 − )nD(n + 1) n+−1 +1=2( − 1)(n + 1)D(n + 2)] :
(2.16)
The special case n = 1 and then = 0 yields in particular ˆ D(1) = 2D(2) − D(3) :
(2.17)
It is stressed that setting = 0 in Eq. (2.16) yields the trivial identity D(n) = D(n) ˆ for all n except for n = 1. Repeating this whole derivation of D(1) with higher-order
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
227
expansions (starting from Eq. (2.11)) yields successively the following higher-order estimates: Order 2 :
ˆ D(1) = 2D(2) − D(3) ;
Order 3 :
ˆ D(1) = 3D(2) − 3D(3) + D(4) ;
Order 4 :
ˆ D(1) = 4D(2) − 6D(3) + 4D(4) − D(5) ;
Order 5 :
ˆ D(1) = 5D(2) − 10D(3) + 10D(4) − 5D(5) + D(6) :
(2.18)
It is amusing to notice that the coecients of the dimensions in (2.18) form a Pascal triangle. The accuracy of these formulas can be examined with the binomial measure for which the entropy dimension is given exactly by D(1) = −(w1 log2 (w1 ) + w2 log2 (w2 )). Using, for instance, the weights w1 = 0:4 and w2 = 0:6, we obtain D(1) = 0:970951, whereas the approximations (2.18) yield successively the estimates {0:968582; 0:971306; 0:971045; 0:970934} that approach the correct value. Formulas (2.18) show that the information dimension D(1) depends on the other integer dimensions, and therefore depends on multipoint correlations. In the above example, the formula of order 3 gives an error smaller than 0.04% which implies that for this multifractal, D(1) is essentially determined by correlation coecients involving 2; 3, and 4 points.
3. Optimal choice of a scaling range The dimensions D(q) are extracted from the generating function with linear regressions done in the scaling range. However, this range is not known a priori, and it must be chosen before the analysis. The choice of a scaling range should be done carefully because it can signi cantly a ect D(q). Moreover, an inadequate choice can lead to an illegitimate function D(q), i.e. to a function that violates some basic properties of multifractals, such as D0 (q)60 for instance. In this section, we propose for the choice of a scaling range some guidelines on which an automatic selection procedure can be based. On one hand, it is desirable to use a wide scaling range to exploit the capacity of multifractal analysis to characterise multipoint statistics. On the other hand, the width of the scaling range is limited by the goodness of t. Indeed, it is normally easier to obtain a good linear t with a narrow range than with a broad range (a good linear t of the generating function can always be achieved at the price of shrinking the scaling range suciently, as long as the error bars do not vanish. Indeed, the generating function is smooth (di erentiable) and is therefore linear locally. From this standpoint, the broadest intervals giving an acceptable t are interesting candidates for the scaling range. Unfortunately, the goodness of t is not always a sucient condition to assess linearity in a given interval.
228
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Indeed, large error bars can result in an acceptable t even if the generating function is far from linear. In this situation, the dimension function D(q) derived from the regressions may not be legitimate in the sense that it could violate some characteristic properties of multifractals, such as the condition D0 (q)60 for instance. To avoid such unsatisfactory scaling ranges, it is necessary to check that some basic properties of mass exponents are indeed satis ed. For this purpose, we propose to use the well-known properties [7] 0 (q)¿0 ;
00 (q)60
(3.1)
as constraints on multifractality. We regard these conditions as minimal requirements of multifractality. In particular, the property D0 (q)60 can be derived from these conditions. The multifractality conditions (3.1) are direct consequences of linearity. Indeed, if the power-law property h(p(n ))q i = nD+ (q)
(3.2)
holds for a range of values of q and n , then taking the logarithm yields ln(h(p(n ))q i) = (D + (q))ln(n ) :
(3.3)
Taking the derivative of (3.3) with respect to q, we obtain h(p(n ))q ln(p(n ))i = 0 (q)ln(n ) : h(p(n ))q i
(3.4)
Because n ¡ 1 and 0¡p(n )¡1, by de nition, it follows from Eq. (3.4) that 0 (q)¿0. Taking another derivative of (3.4) with respect to q, it can be shown that 00 (q)60 is also satis ed. Eq. (3.3) expresses the linearity of the plot ln(h(p(n ))q i) versus ln(n ) for each q, and we have shown that inequalities (3.1) are the direct consequences of this linearity. Conversely, forcing the scaling range to respect these inequalities is an indirect way to obtain linearity. Experience has shown us that a linear regression performed on a function F(q; x) that is not quite close to the ane form F(q; x) = (q)x + c(q), e.g. a function with constant curvature, will typically not give a regression coecient (q) that satis es both the multifractality conditions (3.1). This scaling range selection criterion becomes more stringent when a broader range of q’s is used. Combined with the goodness of t, conditions (3.1) applied to a broad range of q’s typically select intervals that are quite linear. Following the above arguments, we de ne a fractal range to be an interval that gives an acceptable linear t to the generating function (i.e. a reduced chi-square smaller than 1) and that satis es the multifractality conditions (3.1). In general, there will be several fractal ranges of a given size n, and the one that gives the best t (i.e. the smallest value of the reduced chi-square) will be called the best-ÿt fractal range of size n. Recalling that wide scaling ranges (i.e. large n) yield a better characterisation, we nally arrive at the following de nition of an “optimal” scaling range, that embodies all the above argumentation: The optimal scaling range is the widest and best-ÿt fractal range. The automatic selection procedure that corresponds to this de nition is the following. For each interval width n, we consider all subintervals and do for each
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
229
one the linear regression ln( ˆq ()) versus ln(), get the reduced chi-square r2 and the mass exponents (q). Next, we compute numerically 0 (q) and 00 (q). Among the intervals that satisfy the multifractality conditions (3.1), if they exist, we select the one that has the smallest r2 . If this minimum r2 is smaller than 1, then we nd the best- t fractal range of size n. This procedure is repeated for each n, from the smallest (n = 3 consecutive points) to the largest. The widest interval obtained with this procedure is the optimal scaling range chosen for the nal analysis. In the following section, we will use this procedure to locate the optimal scaling range within a given interval. It is stressed that the range obtained depends on the interval chosen for the search. For instance, searching the interval [1; 100] could give the optimal scaling range [50; 100], for instance. In a second step, we can search the interval [1; 50] and obtain another optimal – but narrower – scaling range, e.g. [10; 20]. This iterative procedure allows to nd several valid scaling ranges. It should be stressed that in this paper we adopt a rather pragmatic attitude towards scaling: an interval quali es as a legitimate range if a good t can be obtained with model (3.3), and if inequalities (3.1) are satis ed for all q’s in this range (otherwise the function (q) does not make any sense in the context of multifractals). The procedure may therefore select rather narrow intervals as acceptable scaling ranges. We do not make any statement about the actual nature of the experimental data tted, i.e it could di er from a power law. All we say is that the intervals selected with our method allow to obtain a good t and produce a legitimate multifractal spectrum.
4. First example: analysis of paper formation Paper formation is a technical term that refers to the variations of optical density in a sheet of paper. These variations are easily observed by viewing a paper sheet with a light source behind it. It is the disordered spatial distribution of paper bres that produces the optical density variations and thus e ects the texture of the paper. In the paper industry, it is common practice to digitise the light intensity transmitted through a sheet of paper, and then to analyse this data (formation analysis). A good knowledge of paper formation is of practical signi cance because it can be related to paper’s physical properties such as strength and printability. In this section, we show an example where multifractal analysis can be appropriate for the characterisation of paper formation, and could therefore be a useful complement to the more conventional methods based on Fourier analysis [8] that are used in the paper industry. We have used three paper formation images obtained from three di erent types of papers which were given the symbolic names “sum1”, “ort1” and “kaj1”. The images were digitised with a resolution of 512 × 512, for a physical size of 6:4 cm × 6:4 cm. These paper formations are shown in Fig. 1. It can be seen that the illumination is inhomogeneous, i.e. it tends to be stronger at the centre than on the edges of the image. We partly compensated for this inhomogeneous illumination by subtracting from the image a smoothed version of the same image. The smoothed image was obtained by
230
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Fig. 1. Paper formations images. Top: “kaj1”; middle: “ort1”; bottom: “sum1”.
tting a polynomial of order 4 to each line of the image (Fig. 2a and b). In addition to these formation samples, we constructed an arti cial sample for which the light intensity g(x; y) was sampled for each pixel from a uniform random number generator ranging from 0 and 255.
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
231
Fig. 2. (a) Top: Plot showing the light intensity (irregular curve) together with the fourth-order polynomial tted through it (smooth dotted curve). The polynomial background is subtracted from the signal to correct for the inhomogeneous illumination. (b) bottom: On the left is the original formation “kaj1”, and on the right its transform corrected for the inhomogeneous illumination.
To focus on the variability of the formation image, we applied a transformation to the original images. Let g(x; y) denote the transmitted light intensity at point (x; y) (partly compensated for the inhomogeneous illumination). A transformed signal S(x; y) was obtained by computing the squared norm of the intensity gradient, i.e. S(x; y) = (@g=@x)2 + (@g=@y)2 , and S(x; y) was used to de ne the measure according to (2:1). A typical view of this transformed signal is shown in Fig. 3. The singular appearance of S(x; y) is striking (i.e. spikes having a broad range of amplitudes), which makes this image qualitatively suitable for multifractal analysis, which is a singularity analysis. The relevance of multifractal analysis as a characterisation method is linked to the degree of validity of the power-law property (2:2) (or (2:4) for the generating function).
232
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Fig. 3. A portion of the transformed signal S(x; y) obtained from the formation kaj1. The singular appearance of this signal is striking.
More precisely, the multifractal approach is most appropriate when the scaling range [min ; max ] is broad, thus allowing to characterise variability over a wide range of length scales with a one-variable dimension function D(q). We show in Fig. 4 the generating functions obtained from the three paper samples, plotted for q = 2. These three functions are similar (they look virtually superposed) and exhibit a power-law behaviour (i.e. linear region in log–log coordinates) that starts at the largest scale and extends down to about 10 pixels, which corresponds to about 1.5 decade. These linear segments are the optimal scaling ranges that were obtained with the selection procedure described in Section 3. For the three types of paper formations analysed, we can say that the multifractal model (2:2) is valid for a broad range of length scales. We display in Fig. 5 the generalised dimensions D(q) obtained from the optimal scaling ranges for the three paper samples, with 06q62: For this arti cial sample we obtained D(q) ≈ 2 for each q, which is the signature of a (statistically) uniform image, i.e. an image for which measure and area are proportional (i.e. “monofractal”). The rst obvious observation is that all samples have high values of D(1), i.e. values which are close to 2, thus indicating a rather homogeneous spatial distribution of the measure. Nevertheless, the di erences between the D(q)’s and their uncertainties (error bars) allow to distinguish the di erent formations unambiguously. In our future work we wish to correlate some features of the D(q) curve (such as D(1)) with the physical properties of the paper.
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
233
Fig. 4. Generating functions obtained from three di erent paper formations, displayed for q = 2. The optimal scaling ranges of each sample are marked by squares, circles and crosses. The linear region starts at the largest scales and extends down to about 10 pixels, which corresponds to approximately 1.5 decade.
Fig. 5. Generalised dimensions obtained from each of the paper samples, together with the dimensions obtained from an arti cial random (statistically uniform) sample.
5. Second example: analysis of pore space in sedimentary chalk The rocks analysed are chalk cores extracted from a geological formation in the Valhall eld in the North sea (oil reservoir). We analysed six core plugs labelled 1,3,4,7,9 and 10. The core plugs were used to prepare thin sections suitable for the
234
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Fig. 6. Samples of the thin sections of the rocks analysed.
analysis of the geometry of pore space with scanning electron microscope (SEM). The SEM images of the thin sections were digitised with a resolution of 512 × 512, for a physical size of about 100 m. We obtained four SEM images for each plug. These images were next thresholded in such a way that pore space would be given a black colour, while the rock matrix received a white colour. More details on this data acquisition procedure can be found in [9]. Samples of these images are shown in Fig. 6.
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
235
Fig. 7. Typical generating functions for the chalk rocks. The error bars are small for all samples. The dots mark the optimal scaling ranges. No acceptable scaling range was found at small scales, i.e. between 1 and 10 pixels. Top: The optimal scaling range obtained from the complete interval 1–512 is shown as “range1”, and corresponds to the homogeneous limit, i.e. D(q) ≈ 2 for all q. The optimal scaling range obtained from the interval 1– 40 is shown as “range2”. Bottom: Because of the tiny error bars, we obtained two optimal scaling ranges at large scales (“range1” and “range2”) that both correspond to rather homogeneous behaviour, i.e. D(q) ≈ 2. A third optimal scaling range (“range3”), rather narrow, was found by searching the interval 1–90.
For such geometrical multifractals [4,5], we de ne the measure of a square to be the area of the pore space contained in this square. In other words, the function S(r) used in (2:1) is simply the indicator function of the pore space, i.e. S(r) = 1 if r is the pore space, and S(r) = 0 otherwise. The four images corresponding to each core plug were used to calculate a single generating function according to formula (2:3b). Typical generating functions obtained from these rocks are shown in Fig. 7. We can summarise
236
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
Fig. 8. Generalised dimensions obtained for each of the rock samples.
their behaviour as follows. Searching an optimal scaling range in the whole interval, i.e. 1–512 pixels, we obtained for all rocks a scaling range located in the large scales (e.g. 40 –512) that correspond to the homogeneous limit, 2 i.e. such that D(q) ≈ 2 for all q (monofractal behaviour). Searching the remaining interval (in this case 1– 40) for an optimal scaling range yields another range, narrower than the rst one, e.g. 10 – 40. For all rock samples, it was found that there was no acceptable scaling range between 1 and 10 pixels (the multifractality conditions (3:1) are never satis ed in this interval, even if a good t can be obtained), which are length scales of the order of the typical pore size (about 10 pixels). For several rock samples, we found three disjoint scaling ranges before reaching the length scale of 10 pixels (Fig. 7). The informative behaviour of the generating function, i.e. away from the large-scale homogeneous limit, occurs typically between 1 and 80 pixels, which corresponds to about two decades. The acceptable scaling ranges lie somewhere between 10 and 80 pixels, i.e. above the typical pore size (about 10 pixels) and below the homogeneous limit (about 80 pixels), and the width of these ranges reaches at most three quarters of a decade, i.e. about half of the informative range [1,8]. The D(q)’s obtained from the scaling ranges contained in the interval 10 – 40 were plotted together in Fig. 8. We notice a distinct multifractal behaviour in sedimentary rocks which we also observed in our earlier work [9,10]. Such behaviour can be attributed to the ow properties of these rocks [9], and this is a subject of our future investigations.
2
For a large box size L, the measure px (L) of a box centered on a point x is simply proportional to the box area and to the macroscopic porosity of the medium, and does not vary with x.
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
237
6. Conclusions We proposed guidelines for the choice of the scaling range for multifractal analysis. For paper formations, we discovered a convincing multifractal behaviour, i.e. a power law that extends over most of the range of scales available. Multifractal analysis could therefore be a valuable tool for the characterisation of paper formations. The analysis of the porous rock images showed that the multifractal approach allows to characterise the images variability over approximately half of the informative range of scales. Appendix A Eq. (2.14) (main text) is an expansion of the form M [q] = l+ (q) ∼ =
N X
a i i :
(A.1)
i=1
Using L’Hospital’s rule, we can show that the asymptotic property N X
ai i ∼ !
(A.2)
i=1
as → 1 holds for any N , where the exponent ! is given by !, , N N N X X X i ai ai i ai : ln() = ! = lim ln →1
i=1
i=1
(A.3)
i=1
PN Moreover, it can be shown that i=1 ai = 1 and consequently (A.1) gives for all q an exact estimate of M [q] for = 1, i.e. M [q] = 1. Comparing (A.1) and (A.2) leads us to de ne an estimate (q) ˆ of the mass exponent (q) by ! ≡ 1 + (q). ˆ Using result PN (A.3), we obtain (q) ˆ = −1 + i=1 ai i , which takes the explicit form given by (2:15) in the main text. References [1] R. Haralick et al., Textural features for image classi cation, IEEE Trans. Systems Man Cybernet. SMC 3 (6) (1973) 610 – 621. [2] J.C. Russ, The Image Processing Handbook, Boca Raton, CRC Press, 1995. [3] D. Avnir, D. Lidar, O. Malcai, On the abundance of fractals, in: M.M. Novak, T.G. Dewey (Eds.), Fractals in the Natural and Applied Sciences, World Scienti c, Singapore, 1997. [4] A. Saucier, J. Muller, Characterization of porous media with geometrical multifractals, Fractals (4) (1993) 894 –903. [5] A. Saucier, J. Muller, Remarks on some properties of geometrical multifractals, Physica A 199 (1993) 350–362. [6] T.C. Halsey, M.H. Jensen, L.P. Kadano , I. Procaccia, B.I. Shraiman, Fractal measures and their singularities; The characterization of strange sets, Phys. Rev. A 33 (1986) 1141–1151. [7] C. Meneveau, K.R. Sreenivasan, The multifractal nature of turbulent energy dissipation, J. Fluid. Mech. 224 (1991) 429– 484.
238
A. Saucier, J. Muller / Physica A 267 (1999) 221–238
[8] J. Muller, G. Helgesen, IFE Report IFE=KR=F-90=121. [9] J. Muller, Characterization of the North Sea chalk by multifractal analysis, J. Geophys. Res. 99 (B4) (1994) 7275 –7280. [10] A. Saucier, O.K. Huseby, J. Muller, Electrical texture characterization of dipmeter microresistivity signals using multifractal analysis, J. Geophys. Res. 102 (B5) (1997) 10 327–10 337.