Application of multifractal analysis to segmentation of water bodies in optical satellite images

In this paper, we study the characteristics of multifractal spectra obtained for remote sensing images through the coarse theory and propose a method for segmentation of water bodies based on it. In the first place, spectra of self-similar images created with Iterated Function Systems are calculated and compared with their statistical spectra. Then, optical remote sensing images are studied, emphasizing the differences between their spectra and those obtained for synthetic images. Attention is focused on the concavity of real image spectra and on its interpretation in terms of the analyzed images. This led us to the proposition of a segmentation method for water bodies. The method is tested and the results are compared with water masks of the regions under study. Comments are made about the limitations of the proposed method.


I. Introduction
Fractal theory was initially proposed by Mandelbrot [1] and extensively developed up to these days [2][3][4][5].Multifractal theory can be seen as an extension of fractal formalism, centered in subsets with different scaling behavior coexisting simultaneously.Since its introduction, multifractal analysis has been developed and applied in a wide range of disciplines, like medicine [6], neuroscience [7], economy [8][9][10], ecology [11] and theoretical physics [12][13][14].Many methods have been proposed for approximating multifractal spectra.A summary of the most important ones can be found in [15].
Particularly, multifractal analysis has been used via different formalisms to extract information from images.
In a series of three articles, Arneodo et al. [16][17][18] develop a multifractal analysis method for images based on wavelet coefficients.Their approach leads to a multifractal spectra expressed in terms of Hölder exponent for functions, a solid way of studying singularities.However, the present work leads to using Hölder exponent for measures.Despite the fact that there are theoretical connections between those spectra, the passage from one to the other is not straightforward, making any comparison of numerical results quite complicated.Nevertheless, it is interesting to see how those singularity spectra are useful for characterizing captures with diverse textures.In [16], the authors propose performing a multifractal analysis of images through the socalled wavelet transform modulus maxima method and testing it in random self-affine surfaces.In [17], this group applies the proposed formalism to multifractal synthetic images simulating cloud struc-tures, while in [18], they utilize it on Landsat images with cloud cover.
In a similar line, Wendt et al. present the use of wavelet leaders for studying regularity and obtaining the multifractal spectra of images with different textures [19].
In a classical article, Lévy Véhel and Mignot [20] characterize different types of 2D singularities via Hölder exponents for the sum measure and various capacities.Once the response of these exponents to the singularities is studied, the authors use them to obtain the edges in both, synthetic and real images.
Multifractal analysis based on measures has also been applied to medical images.Nilsson [21] studies the possibility of using this methodology to distinguish healthy tissue from a potentially tumoral one based on the differences in the width of the corresponding spectra.Vasiljevic et al. [22] use particular parameters of multifractal analysis, like some generalized dimensions or the range of Hölder exponents, for classifying the primary cancer in cases of metastasis bone disease.
Millán et al. [23] also pay attention to several selected parameters estimated from multifractal spectra and use them together with a log-Levy stable model for characterizing soil surface moisture distributions.In geophysics, Posadas et al. [24] characterize the dynamics of preferential water flow in soils and porous media with some particular parameters of multifractal spectra of MRI data.In engineering, multifractal spectra has been proposed as a method for evaluating the integrity of concrete shear walls through the study of images of fracture patterns caused by earthquakes [25].They concluded that multifractal parameters move toward higher values as crack patterns extend and grow.
As it has been pointed out, many of the cited works oriented to image analysis do not deal with the multifractal spectrum as a whole, but rather focus only on the values of some parameters to characterize the processes under study.This is a clear difference with the present contribution, which seeks to study the shape of the spectrum for real images and how it is associated with different regions of the captures.
In recent decades, on the other hand, remote sensing has become an extremely useful tool for studying the complex dynamics of the Earth's surface.This technique gives us a huge amount of periodic and varied information about it.
Along with the rapid development of remote sensing, interest in the study of geometrical and self-similarity properties of different elements of the Earth's surface has grown.Some examples are the study of coastlines [26,27], river basins [28], urban sprawls [29,30] and forests [31].Researchers have also turned their attention to the characteristics of changes introduced in the ground by various environmental disasters such as wildfire [32,33], oil spills [34], seismic activity [35] and others.These studies, mostly based on fractal theory, show the intricate nature of the structures present in satellite images captured by different sensors.
Taking all this into account, it is reasonable to use multifractal approach on remote sensing images in order to develop segmentation tools that could allow us to differentiate various ground covers based on their textural features.
In this work, we use the so-called coarse multifractal theory [5] because it is computationally manageable and allows us to study, at the same time, regularity of each pixel of an image and the structures that points of similar regularity determine.This dual approach, local and global, becomes useful for classifying regions according to their textural properties.Such a task would be very hard to carry out by implementing a method that analyzes only globally the properties of an image.
Since this work was done in the framework of an interdisciplinary project devoted to plain rivers modeling and flood predictions, we decided to direct our attention towards segmentation of water bodies.Our motivation is to develop a reliable method for segmenting flooded areas accurately in order to use the obtained information as input in a prediction model.
This goal helped us devote our efforts in examining geometrical properties of multifractal spectra corresponding to real images and to interpret them in terms of geographical components of the scenes.In particular, concavity of spectra caught our attention because particularities on it seemed to be correlated with water covers in the regions of interest.
This article is organized as follows: Mathematical tools used in the analysis of the images are succinctly presented in section II, where central concepts of multifractal theory and Iterated Function Systems are defined.An interpretation of theoretical concepts, which will allow further considera-
tions on results, is also presented at the end of this section.
Section III is mainly devoted to explain how theoretical concepts are adapted for image analysis and segmentation.A procedure for testing multifractal spectra obtained for synthetic images is also explained here.
At the beginning of section IV, we present the results of testing multifractal method with theoretical images.Then, we show the first spectra of remote sensing images emphasizing the differences with the ones corresponding to theoretical images.After analyzing real images, we state a segmentation criterion for flooded areas and present comparisons between its results and water masks of the images.The conclusions of this research are presented in section V.

II. Mathematical tools
In this section, we define the mathematical tools that we are going to use.This is basically done to clarify the notation and check the meaning of some theoretical concepts, which will help interpret results later.For a deeper explanation of the concepts of measure, fine and coarse multifractal theory, we refer the reader to [3,5,36].

i. Measures on R n
A measure over a set X ⊂ R n is defined as a function which associates a value in [0, ∞) to every subset of χ [5].The support of a measure (spt(µ)) is defined as the largest closed subset of the set χ for which every open neighborhood of every point has positive measure.A finite measure µ can be considered as a mass distribution over the support of the measure.By definition, a measure is additive, that is: the measure of a union of two disjoint subsets is the sum of the measures of the subsets.
The class of Borel sets is defined as the smallest collection of subsets of X ⊂ R n which satisfy: 1. every open set and every closed set is a Borel set; 2. the union of every finite or countable collection of Borel sets is a Borel set, and the intersection of every finite or countable collection of Borel sets is a Borel set.
If the Borel sets of a set X can be measured with a measure µ, then that measure is called a Borel measure.
In this article, we are going to use "sum measure": let g(k, l) be the intensity in a grey scale of a point (k, l), and let Ω be the set of all points in a neighborhood of width i centered in a point (m, n), then the measure of Ω is defined: ii. Hölder exponent and fine theory Given a Borel measure µ and a point x in the support of µ, its local dimension [5] in x is if the limit exists, where B(x, r) is a ball of radius r centered in x.This limit is usually called the Hölder exponent (α).For α ≥ 0 we define We define the Hausdorff fine multifractal spectrum f H (α) as the function that maps α to the Hausdorff dimension dim H of the set F α , that is: Although the Hausdorff dimension is more convenient for theoretical issues, in practice, it is almost impossible to estimate from experimental data.That is the reason why the so-called 'Coarse Theory' is used to obtain the multifractal spectrum.

iii. Coarse theory
Once α is calculated for all x ∈ spt(µ) as explained in the previous section , spt(µ) is covered with a regular lattice of n-dimensional boxes of width r (r-mesh).For r > 0 and α ≥ 0 we define N r (α) = #{ boxes of r-mesh such that µ(box) ≥ r α }. ( The coarse multifractal spectrum is defined in [5] as if the double limit exists.

iv. Moment sum and Legendre transformation
Another method to estimate the multifractal spectrum is based on the definition of the so-called partition function.Given a measure µ, its support is covered with an r-mesh.For q ∈ R and r > 0 partition function is defined as the sum of the q powers of the measures µ i of each box in the mesh, where N r is the total number of boxes in the rmesh.Consistently with our previous assumptions, χ q (r) satisfies a power law behavior with respect to its scale [5] at least for small r χ q (r) ∼ r −τ (q) , ( If we discretize the range of α in intervals of width (with sufficiently small), the sum in the definition of χ q (r) can be decomposed into a sum over regions where α is approximatively constant and of size proportional to r −f (α) .It is then easy to see that in the limit r → 0, the sum is dominated by the α(q) that minimizes the exponent qα−f (α).This leads to a relationship between f (α) and τ (q) in terms of Legendre transform, defined for α ≥ 0 by If τ (q) is a convex function, which ensures that the Legendre transform is uniquely determined, it is possible to calculate f (α) inverting the Legendre transform.The spectrum obtained with this procedure is usually known as Legendre Spectrum and denoted f L (α).This method is usually more numerically manageable than the one derived from coarse theory.However, its application is limited to cases in which τ (q) is a convex function of q, otherwise the result is meaningless.

v. Iterated function systems (IFS)
Let D be a closed subset of R n .A mapping S : An iterated function system is defined as a finite family of contractions S m with m ≥ 2. A nonempty compact subset F ⊆ D is called attractor for the IFS if F = m i=1 S i (F ).An IFS determines a unique attractor, which is usually a fractal.
Let us suppose we associate a probability to each contraction of an IFS.If we apply the contractions according to their probabilities on a point in the IFS attractor, we will obtain a self-similar measure.This method is known as Random Iteration Algorithm.
As shown in [5], multifractal spectrum of selfsimilar measures is strictly concave.In addition, proposition 11.7 of [5] states that f C (α) = f L (α) for that type of measures.

vi. Interpretation of α values
A typical multifractal spectrum of a set D constructed by a Random Iteration Algorithm is sketched in Fig. 1.An interpretation of α values arises from the spectra of these IFS.In [5], it is shown that for a self-similar measure created by an IFS {S i }, i = 1, . . ., m, with contraction ratios {c i } and probabilities {p i } If we assume that c i = c ∀i, α max is related to the minimum p i .Assuming also that p i = p j when i = j, there is a single point in D, let us say x αmax , for which α = α max .Every point in D has a probability greater than that point.Applying similar reasoning, the point in D for which α = α min , x αmin is associated with the highest probability.
From Eq. ( 2), balls of radius r centered in x αmax have lower measures as r becomes lower, but that reduction is not proportional to r 2 .Since the set is multifractal, discontinuities in the measures of the neighborhoods occur as r varies, resulting in a reduction of the measure proportional to r αmax , with α max > 2.

III. Description of methods
In this section, we describe the use of mathematical tools presented in section II to carry out the study of images.In the first subsection, we describe the algorithm used for synthesizing multifractal sets of well known properties.We also describe how statistical multifractal spectra of these sets are obtained through the method of moments.In subsection ii, we describe the application of coarse theory to synthetic images and how its results are compared with the statistical spectra.Finally, in subsection iii, the analysis of satellite images through the coarse theory is described.

i. Algorithm for creating synthetic multifractal images
The algorithm used to obtain multifractal sets is explained in a more detailed way in [21].This al-gorithm emulates deterministically the result of a Random Iteration Algorithm over a subset of R 2 with contractions and corresponding probabilities The contraction factor is 0.5 in all cases.
The algorithm begins with a square array of 2 N × 2 N pixels of value 1.The initial square is divided into four squares and each one is multiplied by the probability corresponding to the contraction that each square represents.This process is iterated over the resulting squares up to N iterations (resolution limit).
Images of 256 × 256 pixels were synthesized with this algorithm based on random probabilities vectors.
Because of the way in which these images are created, τ (q) must be monotonically concave in Eq. ( 9), allowing the use of Legendre transform for determining their multifractal spectra.Different meshes with boxes 64, 128 and 256 pixels wide were applied.Given a mesh, the sum of the values of the pixels within each box is calculated.These measures are raised to an exponent q and added, obtaining the value of the partition function for that mesh width and the exponent q.After making this calculation for each mesh, a linear fit of the logarithm of the partition functions versus the negative of the logarithms of the corresponding box widths is performed.This gives us τ (q) of Eq. ( 8) for each q.After applying Legendre transform to τ (q), a statistical estimate of the multifractal spectrum is obtained.
As we said in section II.v,Legendre spectra of self-similar synthetic sets must agree with coarse multifractal spectra.That is the reason for using the synthetic images as a test for the multifractal method based on coarse theory.

ii. Analysis of synthetic images through coarse theory
We choose B(x, r) in Eq. ( 2) to be square neighborhoods of (2k − 1) × (2k − 1) pixels, with k = 2, . . ., 8.Those neighborhoods are centered in the pixel whose Hölder exponent we are trying to calculate.The measure of all neighborhoods is calculated.A linear fit is then performed on the logarithms of the measures versus the logarithms of the neighborhood widths.The slope of the adjusted line is our estimate of α for that pixel.A padding of copies of the image is built for performing calculations on pixels next to the edges.Once α values for every pixel are determined, the range of α is divided into R = 10 classes bounded by each pair of successive elements of the list (α min , α min + δα 2 + s • δα, α max ), with δα = (α max − α min )/R and s = 0, . . ., R − 2. Those classes determine subsets of pixels with similar local regularity.Then, we calculate the box counting dimension of each subset.Different meshes are applied to the image and the number of boxes necessary to completely cover the selected subset is counted.For each subset, we make a linear fit of the logarithm of the number of boxes versus minus the logarithm of the width of the boxes.The slope of the fit is our estimate of f (α m ), where α m is the mean value of each class.The boxes of the meshes used in the images have 4, 8, 16, 32, 64, 128 and 256 pixels of width.
It was checked that the multifractal spectra obtained had uniform concavity and that they fulfilled proposition 11.7 of [5].Although proposition 11.7 states that f L (α) = f C (α) for a self-similar measure, given the discretizations made in the calculation of the spectrum through the coarse theory, it is reasonable that f C values obtained are lower than those of f L .Particularly, it was observed that by reducing R in calculation of f C , a better match between the spectra is obtained.
Errors in calculations of α for each pixel are obtained.Errors in α exponents will affect the errors in multifractal spectrum values.This is primarily to the possible change in the full range of α values found in the picture.In second place, if the pixels have associated different exponents, when the total range is divided into sub-ranges, this could change the set of pixels included in a class, altering the fractal dimension of the set.To account for this, multifractal spectra for three sets of α were calculated: the set of α obtained, and the sets resulting from adding and subtracting the corresponding error to the α value of each pixel.Spectra of these three configurations are obtained with their respec-tive errors.
The largest α error between the pixels of a class is taken as the error bound of α for that class.For each class, the error bound of f (α) is given by the distance between that point and the farthest error bound for that class.
Concavity is analyzed by calculating an approximation of the second derivative of f (α).Thus, for each point, it was calculated If ∆ 2 f (α i ) is negative for all i, then the function is concave downward.

iii. Analysis of satellite images
Optical images taken by OLI (Operational Land Imager) and TIRS (Thermal Infrared Sensor) sensors of Landsat 8 satellite, provided by the U.S. Geological Survey [37] were studied.Since the work was oriented to the detection of flooded regions, the images were selected for having remarkable water bodies and a low percentage of cloud cover.Analysis was performed particularly in fragments of band 5 (wavelength of 0.851-0.879µm) of those images.
Images were calibrated before being studied so that their pixels represent values of reflectance measured by the sensor in each position.For more information on the Landsat 8 platform we refer to [38].
Since it was found that τ (q) did not have a uniform concavity for satellite images, these images were not studied using the method of moments, and their spectra were obtained only through the coarse theory.
As for the calculations of multifractal spectra in synthetic images, we chose the B(x, r) in Eq. ( 2) to be square neighborhoods of (2k−1)×(2k−1) pixels, with k = 2, . . ., k end .Calculations on the pixels near the border of the analyzed region are made using adjacent regions of the image as padding.
The range of α is also divided into R classes bounded by each pair of successive elements of the list (α min , α min + δα 2 + s • δα, α max ), with δα = (α max − α min )/R and s = 0, . . ., R − 2. Box counting dimension is obtained in each subset with the procedure explained early.
Selection of k end , R and the size of boxes in the meshes used in box counting will be discussed in the next section.
The classification of pixels of images as flooded or not flooded is achieved through the application of thresholds on the multifractal spectra.Selection of those thresholds arose during the exploratory study of the graphs of f (α) obtained and will be detailed in section IV.

IV. Results and Discussion
Multifractal spectra of 600 synthetic images were obtained through coarse theory and compared with their statistical spectra.591 spectra showed uniform concavity and fulfilled f C ≤ f L .Two examples are presented in Fig. 2.
Taking into account the results of calculations on synthetic images, we decided to analyze some remote sensing images through coarse theory.Initially, we took fragments of 256 × 256 pixels and analyzed them with the same configuration as the theoretical multifractals.Four examples are presented in Fig. 4. We noticed that concavity was not uniform in the spectra of those images.Particularly, it seemed that for images where there were large water bodies such as lakes, this change in concavity occurred more markedly.
Calculation of multifractal spectra in case of images built with IFS requires taking as maximum neighborhood width the one that most closely matches the full width of the image because, otherwise, the result does not fit the statistical spectra.We suppose that this is due to the infinite nature of fractals.Therefore, to take all possible intervals within the limits of scale means to collect as much information as possible.
As evidence of this, we have performed some calculations on fractals created by IFS with 2048 × 2048 pixels, which are not showed here.When the points used in the calculus of α were plotted, we observed that, to a certain scale, the points seemed to have a slope.Then, it changed abruptly as a result of discontinuity in the measurement.This new slope was retained during some scales and then a new change occurred.However, when a larger number of scales was observed, a linear behavior was seen in the points, despite these apparent local changes in slope.Comparisons between multifractal spectra obtained for two synthetic images through the method of moments (red dots) and the coarse theory (blue dots).Black dots correspond to multifractal spectra resulting from α sets calculated by adding and subtracting the corresponding error to the α value of each pixel.They are used in the determination of error bonds.
The case of real images is quite different since they do not have an infinite number of possible scales.Therefore, there is no need to take all possible neighborhoods for calculating the Hölder exponent.
However, given that this work is oriented to take advantage of multifractal formalism for the detection and segmentation of water bodies, the choice of k end is not trivial.
Before addressing the issue of the selection of k end , we should explain the criterion for segmentation of water bodies proposed in this work.For that, note the analysis presented in Fig. 5.It was done on an image of 1024 × 1024 pixels with Returning to the selection of k end , we have noticed that if the largest neighborhood for each pixel in a flooded region is large enough to cover an edge of the water body, the pixels will have α ≥ 2. On the contrary, if for example the largest neighborhood centered in a pixel in the middle of a lagoon does not touch the edge of it, then that pixel will have α ≈ 2. Figure 6 presents an analysis of the same region as in Fig. 5, with k end = 8 and R = 30.It is clear that the pixels from the center of the water body have an α ≈ 2, as opposed to what happens in Fig. 5.In that situation, those pixels would not be classified as water.For this reason, it is necessary to insist on the fact that k end should be chosen to avoid this situation.Clearly, the proposed criterion is only applicable in cases where the multifractal spectrum possesses that particular form.For example, a water body such as a river will have a high α but the fractal dimension of this region will be approximately 1. Therefore, there may be no change in concavity in that case.On the contrary, if besides a river we also have a lagoon or a water body with a fractal dimension greater than one and the greatest neighborhoods are large enough to cover any of its edges, then we will find a region with a high α and a fractal dimension greater than 1, causing a minimum in the center of the spectrum.With all that in mind, we now present segmentations performed on some Landsat images suitably selected through visual inspection.Images are presented and described in Fig. 7.
Figures 8, 9, 10 and 11 show the results of the analysis.Values of α center found in spectra are       As a way of evaluating the efficiency of the proposed segmentation method, the results were compared with water masks of the regions, which were provided by our colleagues of The National Water Institute (INA): L. Giordano and J. Bianchi.Comparisons were made through confusion matrices and the calculation of Cohen's kappa coefficient for the correlation between segmentations and masks.
Confusion matrices are presented in Tables 1, 2, 3 and 4, while table 6 summarizes kappa values obtained.

V. Conclusions
In this work, we have studied properties of coarse multifractal spectra calculated for optical satellite images and we have compared them with spectra corresponding to self-similar images.We have analyzed the concavity of natural images spectra and interpreted it in terms of structures present in the images.
We have seen that, in presence of suitable water bodies, a local minimum exists in the multifractal spectra.Thanks to the dual approach of coarse multifractal analysis, we could associate Hölder exponents higher than the local minimum with regions covered by water.
Based on those results, we have introduced a new criterion for segmentation of water bodies.Of course, this criterion is far from universal.As we mentioned early, the particular shape of spectra used in segmentation depends on the presence of given structures in the studied scenes.Without those conditions, we are not able to perform an accurate segmentation.
We must remark that the whole analysis could not have been done by means of a statistical or global approach because we would have needed also a local characterization of pixels in order to detect the correlation that allowed us to state the segmentation criterion.
After discussing the applicability of the method, we tested the accuracy of the results by comparing them with water masks.The segmentations performed showed a high coincidence with the 'true' flooded regions, situation evidenced by Cohen's kappa values close to 1 presented in section IV.
It is interesting to note that the proposed method
is only applied in this work to one band of spectral information, but there is no limitation for its application to other bands or even sensors of different nature.It is reasonable then to assume that the study of the diverse remote sensing information currently available could lead to a more complete characterization of a region of interest.
Another possibility that could further extend the usefulness of this method is the implementation of different measures or even capabilities in addition to the one presented here.This could allow the detection of various structures or regions in an image, leading to more detailed segmentations.

Figure 1 :
Figure 1: Typical multifractal spectrum of a set constructed by a Random Iteration Algorithm.

Figure 2 :
Figure2: Comparisons between multifractal spectra obtained for two synthetic images through the method of moments (red dots) and the coarse theory (blue dots).Black dots correspond to multifractal spectra resulting from α sets calculated by adding and subtracting the corresponding error to the α value of each pixel.They are used in the determination of error bonds.

Figure 5 :
Figure 5: Analysis performed on an image of 1024× 1024 pixels with k end = 10 and R = 30.

Figure 6 :
Figure 6: Analysis performed on an image of 1024× 1024 pixels with k end = 8 and R = 30.

Table 5 :
Values of α center found in multifractal spectra of analyzed images.