key: cord-0201850-7ur3kqno
authors: Banerjee, Debasmita; Jha, Amit Kumar; Iyengar, A.N.Sekar; Janaki, M. S.
title: Fractal dimension analysis of spatio-temporal patterns using image processing and nonlinear time-series analysis
date: 2020-10-31
journal: nan
DOI: nan
sha: 81222bb1d648958fce1b760e0457a75a053f43b6
doc_id: 201850
cord_uid: 7ur3kqno

This article deals with the estimation of fractal dimension of spatio-temporal patterns that are generated by numerically solving the Swift Hohenberg (SH) equation. The patterns were converted into a spatial series (analogous to time series) which were shown to be chaotic by evaluating the largest Lyapunov exponent. We have applied several nonlinear time-series analysis techniques like Detrended fluctuation and Rescaled range on these spatial data to obtain Hurst exponent values that reveal spatial series data to be long range correlated. We have estimated fractal dimension from the Hurst and power law exponent and found the value lying between 1 and 2. The novelty of our approach lies in estimating fractal dimension using image to data conversion and spatial series analysis techniques, crucial for experimentally obtained images.

The ubiquitous presence of fractals have long motivated the scientific community to quantify nature's geometrical complexity. According to historical records, the term was first introduced by Benoit Mandelbrot in 1975 to expand the concept of theoretical fractional dimensions to geometric patterns in nature [1] . A defining characteristic of fractals is its infinite number of scale lengths. Due to this, the structure of a fractal is identical in response to random magnification or scale variation. A mathematical characterizing quantity of fractals is called fractal dimension (FD) that measures complexity of patterns as a ratio of the change in structure to the change in scale. Our study involves estimation of the Fractal Dimension of the spatial data generated from the evolution patterns of Swift Hohenberg equation. This equation is a parabolic partial differential equation of fourth order known for modeling patterns in simple fluids (e.g. Rayleigh-Bénard convection) and in a diverse collection of complex fluids and biological formation. Previous theoretical measurements of fractal dimensions have been mainly obtained by box counting method [2] which involves analyzing complex patterns by scanning and breaking dataset ( or images ) in finer pieces. In our work, applying numerical solving techniques, we have obtained non-linear spatial series as solution of the SH equation. These non-linear spatial series are then converted into color plots that give patterned images.

We have further retrieved the same non-linear spatial series by converting the image intensity into a spatial series (analogous to time series) using image processing method and then performed non-linear spatial series analysis. The purpose was to demonstrate the usefulness of the method in case of experimental systems when the governing equation is not known and one only has images of the pattern but want to obtain non-linear spatial series for fractal dimension calculation. Analysis of the fractal dimension of the spatial series from the power law exponents reveals a non-integer value between 1 to 2 for our data. The merit of this alternative method also lies in its compatibility with other nonlinear time series analysis techniques such as Recurrence analysis etc. Having extensive applications on experimental images from varied sections of sciences including evolution biology and surface physics, makes it specifically appealing. The system shows long range correlation as the Hurst exponent generated by both the Detrended Fluctuation and Rescaled Range Analysis follows the range between 0.51 and 0.87. We have also calculated fractal dimension from the Hurst exponents which provides a comparable result with that of Hurst exponents, estimated from power law coefficient.

The presence of chaos has recently been reported in the lab-grown experiment to study pattern in two dimensions [3, 4] where the authors were able to detect system transition from long range correlation to anticorrelation between the initial and final state of the experimental patterns using a novel image processing technique. Our current study is closely based on the proposed image to data conversion technique of above mentioned work.

The theory of chaos in a finite-dimensional dynamical system such as in nonlinear ordinary differential equations [5, 6, 7, 8] has been extensively studied. However, exploration of deterministic chaos in partial differential equations (PDE) offer challenging problems [9, 10, 11] . An investigation on the simplest nonlinear PDE reveals that the Kuramoto-Sivashinsky is the simplest known PDE with quadratic or cubic nonlinearity and periodic boundary conditions to show chaos [12, 13] . We have used a model equation known as Swift Hohenberg (SH) equation of the fourth order [14] parabolic type that is famous for its use in pattern formation [15] due to finite length instability, such as in Rayleigh-Bénard convection [16, 17] . Due to the vast applicability of the SH equation, it is represented as a generic model for the spatiotemporal dynamics of spatially extended systems [18] . It is well known that ordinary complicated mixing show multifractality (for two-components or multi-components fluids) [19] . Being a model pattern forming equation for simple and complex fluids, the idea can also be applied to the patterns generated by SH equation [20] if the evolution dynamics show scale invariance with the existence of self-regulatory internal mechanisms that drive the system spontaneously to a statistical stationary state. The motivation for choosing SH model lies in the statements of Cross' and Hohenberg's work [21] where they have emphasized the omnipresent presence of self-similar systems over a wide range of parameter in natural nonequilibrium systems carefully tested by Bak et al. [22] . Also, spatio-temporal chaos rarely exhibits scaling behavior and that occurs only due to parameter tuning. Considering such occurrences of multi-fractal scaling properties conditioned to parameter tuning [23, 24, 25, 26] , we have considered SH equation to be a potential candidate for exhibiting multi-fractality in the ensuing steps of the investigation.

Here, the patterns are nothing but the solution of the SH equation which has been represented using MATLAB colorbar. Time evolution is generated by obtaining snapshots at different time instants whereas space evolution is determined by varying the parameter values. The color difference between each unit pattern points to a different intensity value which in turn serves as the spatial series data. The data exhibits a robust power law behavior where the log-log plot shows two scaling regions with two distinct slopes. The range of power law scaling factor, β varies from 1.3 to 2.8.

In Section II, we describe the SH equation and the method used to numerically solve it. Section III explains the method used to extract the data points for the spatial series generated from the patterns formed by the SH equation. Section IV shows the spatial series and its corresponding phase space plots. In Section V, we present results of the long range correlations of the chaotic data obtained by Detrended fluctuation analysis (DFA) and Rescaled range analysis (R/S). In section VI, power spectrum by fast Fourier transform and its log power-log frequency plots are shown. The presence of more than one slope in power law diagram confirms multifractality in the SH equation. In Section VII, we have presented the summary of the results [27, 28] .

2 Patterns generated from the numerically approximated solutions of the SH Equation

We have used the reduced third order form of the SH equation where z variable has been omitted. The Swift Hohenberg equation [29] is described as

In this equation, is the real bifurcation parameter whereas v is the strength of the quadratic equation. Together with ψ 3 the term vψ 2 − ψ 3 makes the nonlinear term in the equation. Here, the ∇ 2 term in the equation is called a spatial derivative complex. Using the Fourier spectral method [30] and Exponential Time Differencing method [31] of order 2 the numerically approximated solutions of the SH equation have been obtained. The Fourier spectral method has been used to transform the PDEs into ODEs by calculating the Fourier modes. After that, the ODEs are solved using the Exponential Time Differencing method (ETD2). On applying these methods to our equation, a matrix containing the intensity for the x-y plane is generated.

The entire simulation has been carried out using MATLAB [32] . At t=0, the initial condition is chosen to be a random matrix.

Ψ(X(t)) = Ψ(X(0)) (2)

X signifies a two dimensional space i.e. x and y space plane. Here Ψ(X(0)) is a random number from the interval − √ ≤ Ψ(X(0)) ≤ √ . Here, v and values are the only control parameters which can be varied to observe spatial changes. To notice the time evolution, we have recorded the image changes at t=5 seconds, t=10 seconds, t=15 seconds, and t=20 seconds. For our simulation we have chosen 256 x 256 grid points with a time step of 0.1.

Images 1-16 of Figure 1 are the spatio-temporal evolution of SH equation with a defined initial condition and a set of parameters i.e. v and and t. The red disks represent the highest image intensity value (portrayed by jet colorbar, MATLAB) whereas blue represents the lowest value. The next section has been dedicated to explaining the conversion of image intensity to spatial data.

Initially, we have obtained numerically approximated solution from the SH equation and have reconstructed the solutions to images. This makes the spatio-temporal evolution and its pattern forming characteristics to become visually perceptive. One of the basic characteristics of image processing in MATLAB is that an image can be converted into a matrix and the opposite also holds true. This unique convertibility has been used here to extract back the spatial data from the images where each pixel (color grain) corresponds to a lower or a higher value of the matrix. Here, the solution of the equation is the intensity values of images represented by ψ. For this simulation, "Jet" colorbar has been chosen to represent the patterns in Matlab figures. The key significance of colorbar is its usability in defining color scheme for many types of visualizations, such as surfaces and patches. The intensity of the images is represented in such a way that a redder appearance shows a higher intensity value whereas a bluer appearance shows a lower intensity value. Accompanying colorbars with Figure. 1 describe the entire range.

Seemingly redundant, this double conversion (data to image and image to data) has its use in experimental systems, where the spatial data is not known. The technique can be very convenient for obtaining information directly from experimentally generated 3-dimensional images such as pattern formation in developmental biology where intensity of a picture is of great importance along with 2 space dimensions. Once converted from patterns,the data has multitude of applications including their use for traditional nonlinear time series analysis. We are proposing this method for medical imaging, vegetation patterns, Turing patterns, reactiondiffusion systems and any other pattern formation variants where space, time or space-time based evolution is experienced. A similar but experimental approach i.e. image to data conversion for studying chaotic dynamics of microscale patterns formed due to demixing of two fluids has been reported in [33] where they have noticed a transition from long range correlation to long range anti-correlation as the system progresses towards the final state.

On the next step we have applied the nonlinear time series analysis on the extracted spatial series data and have estimated the fractal dimension of the system. The two dimensions of the patterns i.e. the abscissa X and ordinate Y of the patterns are the two space variables x and y. Each unit intensity of an image is the solution of SH equation for two space variables namely x and y at a certain time where the two varying parameter namely v and is assigned to a specific value. The method proposed in this work is a simple image to matrix conversion method, where these color intensity values represent specific colors as given in the colorbar. Now, each row of the image is taken and then concatenated with the next immediate row. Once the entire 2D image is converted into a single row data, it is plotted as a spatial series where y-axis represents the intensity values and x axis represents the data length. A flow chart (Figure 4 ), explaining the method is also attached for convenience.

Since we have obtained a set of pattern forming images depending on varying time and parameter values, we have taken a methodical approach to choose three images out of that set to be used further for our analysis. The process of selection is made in such a way, that the change of pattern is noticed due to the change of and time. For this experiment, we have taken Image 1 as our initial sample which has v = 0.1, = 0.3, t = 5 seconds. Along with it we have taken Image 4 (v = 0.1, = 0.3, t = 20 seconds -time increased), and Image 12 (v = 0.1, = −0.3, t = 20 seconds -and time increased) as a representation of spatio-temporal evolution due to the difference in time and one of the varying control parameter . After acquiring the entire spatial series, we have looked at the 50th row data and 350th row data of each of the aforementioned images.

The intensity profile of the solution depicted in images are converted to data points that generated a specific row vs column matrix. In this paper, we have used a direct approach to convert the image-generated matrix to a spatial series by simply translating data points into a 2D spatial plot. If the matrix has m columns and n rows then it is first converted into a single row containing n × m columns i.e. 1D array by simply removing the rows from second row onward, until the n th row and concatenating them with the first row. When plotted ( Figure. 2), the result is a spatial series with image intensity solution in its Y axis and data length in its X axis. From the plots it is apparent that the data do not hold a periodic nature and gives a typical random appearance. However, apparent randomness is not enough to justify the data and therefore the data has been used for further analysis.

In dynamical systems, phase space projection plays a decisive role in determining all possible states of a system. We have obtained a normalized differential phase space plot via plotting a graph between spatial series data points versus its derivative. The phase space projections represented in the Figure 3 show trajectories of a dynamical system for different row data. The data used is filtrated from noise by making a use of low pass Butterworth filter. The noise is introduced when the spatial series data is differentiated. This happens because of the sharp changes in the spatial series data where the derivative values can not be computed. The low pass Butterworth filter removes those high frequency fluctuation which are redundant for our calculation method.

One of the primary points to be noted from the phase space results is the position of the periodical trajectories or orbits. The trajectory encircles around one basin of attraction [34] and suddenly it starts to circle another. The trajectories represent a very dense pathway which does not follow the same route again however co-exist with its neighbor in a close distance making the area thicker in appearance. This corroborates with the concept of an attractor which by definition means the set of states, invariant under dynamics, towards which neighboring states in a given basin of attraction asymptotically approach in the course of dynamic evolution. [35] As commonly known, this behavior is natural to chaotic data sets. Although it does not offer any specific information about the dataset but significantly affirms the nature of the data being chaotic. One of many other popular approaches to detect chaos has been assessing whether the largest Lyapunov exponent (λ) gives a positive value. The Lyapunov exponent of a dynamical system is a quantity that characterizes the rate of separation of infinitesimally close trajectories. If λ > 0, the nearby trajectories separate exponentially, determining chaos. Using the TISEAN package [36] 

Detrended fluctuation analysis [40] tool, applied over the same spatial data, is a method to investigate the statistical self-similarity of a signal and has been proven to be useful in revealing the extent of long-range correlations in time series [41, 42] . A general feature experienced in all self-similar systems is that the temporal and/or spatial power-law correlations extend over several decades where intuitively one may anticipate that the physical laws would change drastically. A simplified and general definition characterizes a time series as stationary (in this case spatial series) when the mean, standard deviation and higher moments, as well as the correlation functions are invariant under time translation. Signals that do not obey these conditions are non stationary. However the upper mentioned Detrended fluctuation analysis allows the discovery of intrinsic self-similarity hidden inside an apparently non-stationary time series or spatial series. The method is a technique for measuring the same power law scaling observed through R/S analysis, explained later but it has been used specifically to address non stationaries. DFA is commonly known as enabling correct estimation of the power law scaling (Hurst exponent) of a systems' signal in the presence of (extrinsic) non stationaries while eliminating spurious detection of long-range dependence.

The algorithm can be explained in four steps as follows:-1. Integration: The spatial series data x(i) is shifted by the meanx and cumulatively summed.

2. Segmentation and Fitting: The obtained sum y(i) is then segmented into boxes of n equal length. Then a local trend is calculated by fitting polynomial y n (i) to each box.

3. Detrending: The detrended series is then calculated by the following equation.

4. Deriving DFA fluctuation function: The root mean square deviation from the trend, i.e. the fluctuation function is calculated following the equation below.

Where N = entire length of the series Following the power law correlation, self affinity of the spatial series can be represented as

where α, the scaling constant is obtained via the slope of a straight line fit to the log-log graph of n vs F (n). For our case, α = H + 1 [43, 44] .

The rescaled range analysis [45] is a statistical tool to measure the variability of a time series (or spatial series) introduced by the British hydrologist Harold Edwin Hurst (1880 Hurst ( -1978 [46] . The R/S analysis determines how the apparent variability of a series changes with the length of the time-period being considered.

Introduced by Hurst, the Rescaled range (R/S) analysis [47] splits a time series (or spatial series) into adjacent windows and inspects the range R of the integrated fluctuations, rescaling by the standard deviation S as a function of window size [48] . Compared to the DFA analysis examined above, the R/S analysis gives attention to the range of signal instead of finding fluctuations around trend. As a result, DFA is suitable to be applicable for nonstationary time series analysis in contrast to R/S. Here we applied both of the methods to the spatial series data we have received from the n × m matrix to draw both the DFA and R/S plots.

The R/S analysis can be performed by following the mentioned algorithm below.

1. Shifting: The spatial series data is shifted by the mean for X = X 1 , X 2 , X 3 ...X n Y t = X t −X

whereX = 1 n n i=1 X t and t = 1, 2, ..., n 2. Series formation: A cumulative deviate series Z has been calculated

3. Calculation of Range:

The standard devitation can be calculated as follows -

where m = mean of (X 1 , X 2 , ...X t ) 5. Rescaled Range series estimation: The rescaled range analysis can be done as follows -

where t = 1, 2, ...n The calculated ratio of R/S follows the relation R/S = d H where d is a constant and H is the Hurst Exponent. The slope of log(R/S) vs log(n) gives the value of H.

From the plot of logF (n) vs log(n) and log(R/S) and log(n), presented in Figure 6 , we have performed a straight line fitting from which the estimated Hurst Exponent (α) turns out to be greater than 0.5. It has been noticed that a larger Hurst exponent value shows a stronger trend. According to our result, the H value obtained from DFA stays between 0.5458 to 0.8108 whereas the same obtained from RS Analysis stays between 0.5109 to 0.8118. Both the results corroborate with the result obtained from power spectral density analysis and further strengthens the appearance of long range correlation [49, 50] . 6 Self-similarity and multifractality in spatio-temporal patterns A system going through evolution due to time exposure and varying control parameter values is famously known to produce a certain behavior which is not unique to a specific system, rather shared by a large class of systems [51] . Bak, Tang, Wiesenfeld [52] in 1987 have discovered that dynamical systems with many interacting constituents indeed may exhibit a general characteristic behavior in which the dynamical systems organize themselves in a complex method such that no single time or length exist to control the temporal evolution of the system producing scale-invariance. Surprisingly however, the associated statistical characteristics of the complex response can be explained using a simple power law.

Spectral analysis is a reversible technique via which a function is decomposed into a superposition of components, each with a specific temporal (or spatial) frequency. [53] In other words, it samples or chops down a signal over a period of time (or space) and divides it into its frequency components. One of the processes to study the patterns underlying in the frequency distribution is Fourier Transform, which is a linear method. With the help of Discrete Fourier Transform (DFT) technique, achieved efficiently using FFT algorithm, we have obtained Fourier Periodgram or power versus frequency plot of the available spatial data.

The purpose of this estimation is to find the periodicities in the data, by observing peaks at the frequencies corresponding to the periodicities. The divided components are single sinusoidal oscillations at distinct frequencies having their own amplitude and phase information. For a particular dataset, if the power versus frequency plot gives two sharp peaks then that means two frequency components are dominant in the data set.

It is well-known that a power-frequency curve of chaotic signal shows a continuous spectra over a limited range and the energy is spread over a wider bandwidth. On the other hand, power vs frequency plot of a periodic signal presents discrete spectra, where a finite number of frequencies contribute to the response.

In our investigation, Fourier Periodgram plot shows a continues broadband spectra for a specific range which in turn suggests the data being chaotic in nature, for the set of parameters chosen for this analysis [54, 55, 56] . Although Fourier analysis is known to be useful for measuring the frequency content of stationary or transient signals, they fail to distinguish one chaotic data from the other. As per the Fourier periodgram plots ( Figure 5 ) received in our case, we see a broadband appearance of frequencies. In this case frequency means the number of times a pattern has been repeated over a specific data length.

The power law can be described as a functional relationship between two quantities, where a relative change in one quantity results in a proportional relative change in the other quantity. Here, one quantity is related to power of another quantity. The main reason to bring power spectral density along with Detrended fluctuation and Rescaled range analyses explained immediately after, is to understand the multifractal nature of the data.

Mathematically, power spectral density [57] represented by s(f ) satisfies the following relation -

where s(f) represents power spectral density and β represents the power law scaling factor [58] . In case of our data β = 2H + 1 where H is the Hurst exponent [59] . This behavior produces the linear relationship when logarithms are taken of both s(f ) and f , and the straight-line on the log-log plot is often called the signature of a power law. The power law scaling factor β can be expressed as the roughness of a spatial series (or time series) with smoother self-similar series having larger values of β. 0 ≤ β ≤ 1 translates to weak long range correlation whereas 1 ≤ β ≤ 2 signifies strong long range correlation. Another way to quantify the fractal nature of the data is given by fractal dimension [60] (D), that comments about the density of fractals. It is expressed by the following equation -

where γ = |β| The power law scaling factor (β) of our data in FIG. 7 shows two different values for two specific regions of a log-log plot. The Beta values for slope 1 and slope 2 of 50th row and 350th row of Image 1, Image 4 , Image 12 come out to be -1. 48 

The central idea behind this work is calculating fractal dimension of the spatial series (analogous to time series) generated from SH equation for identifying possible presence of fractal behavior. The spatio-temporal patterns were converted into a spatial series data. The chaotic behaviour of this data has been established by estimating the positive Lyapunov exponents, phase space and power spectrum analyses. With the help of Detrended fluctuation and Rescaled range analyses, we have calculated the Hurst exponent and found the value to be greater than 0.5. This indicates that the data is long range correlated [61] . Estimation of the Hurst exponent from both DFA and R/S analyses lead to the calculation of Fractal Dimensions that lies ∈ (1, 2). Recently, both long range correlation and chaotic behavior [62] has been observed by researchers while analysing the genomes of SARS-CoV (2002) and SARS-CoV-2 (2019), responsible for the global pandemic . Our method can aid in analysing similar nonlinear time series and spatial series.

Further, we have carried out an estimation of Fourier power law analysis which confirms that the exponent value (β) lies between −1.04 to −2.57 for the entire data set. We have also calculated the fractal dimension from the power law exponent values and found them comparable to the data generated from DFA and R/S analysis. The power law (log power-log frequency plot) scaling factor highlights the presence of more than one exponents. Whenever the fractal dimension, as a single exponent is not enough to describe the complexity of the system, multifractality can generalize the complexity by using several exponent, capturing its dynamical essence. As seen in our case, different scaling exponents were essential for different segments of the same spatial series, indicating occurrence of the scaling behaviour and characterizing multifractality in our data set. In the patterns, this may indicate interlinked fractal subsets in the spatial series [63] and hence, our work can be used for analysing and studying systems with interlinked fractal subsets and multifractality behavior. 

The Fractal Geometry of Nature

Chaos in Ordinary Differential Equations

Chaos in Partial Differential Equations

Asian Pacific Confederation of Chemical Engineering congress program and abstracts

Environmental Data Analysis with Matlab

Proceedings of the Second IASTED International Conference on Financial Engineering and Applications

Self-Organized Criticality: Emergent Complex Behavior in Physical and Biological Systems