*3.1. Multifractal Spectrum Based on MF-DFA*

The method of MF-DFA is widely applied in scaling analysis due to its high accuracy and easy implementation. For grayscale images, it is not appropriate to calculate the multifractal spectrum with a gray series by the approach of one-dimensional MF-DFA. Generalizing the one-dimensional method to two-dimensional one can better express the information of the surface with self-similar properties [29]. Specifically, the process of calculating the multifractal spectrum of the grayscale image by using the two-dimensional MF-DFA method is determined as follows.

Step 1. Regard a microscopic image as a self-similar surface with a size of *M* × *N*, which is represented by a matrix *X*(*i*, *j*), *i* = 1, 2, ..., *M* and *j* = 1, 2, ..., *N*. Partition the surface into *Ms* × *Ns* (*Ms* ≡ [*M*/*s*], *Ns* ≡ [*N*/*s*]) none-overlapping square subdomains of equal length *s*. Each subdomain is denoted by *Xm*,*<sup>n</sup>* = *Xm*,*n*(*i*, *j*) with *Xm*,*n*(*i*, *j*) = *X*(*r* + *i*, *t* + *j*) for 1 ≤ *i*, *j* ≤ *s* where *r* = (*m* − 1)*s*, *t* = (*n* − 1)*s*.

Step 2. For each subdomain *Xm*,*n*, the cumulative sum is constructed as follows

$$G\_{m,n}(i,j) = \sum\_{k\_1=1}^{i} \sum\_{k\_2=1}^{j} X\_{m,n}(k\_1, k\_2),\tag{1}$$

where 1 ≤ *i*, *j* ≤ *s*, *m* = 1, 2, ..., *Ms*, *n* = 1, 2, ..., *Ns*. Note that *Gm*,*<sup>n</sup>* = *Gm*,*n*(*i*, *j*)(*i*, *j* = 1, 2, ..,*s*) is a surface.

Step 3. The local trend *G*˜*m*,*<sup>n</sup>* for each surface *Gm*,*<sup>n</sup>* can be obtained by fitting it with a pre-chosen bivariate polynomial function. In this paper, we adopt the trending function as

$$G\_{m,n}(i,j) = ai + bj + c \tag{2}$$

where *a*, *b*, and *c* are free parameters to be estimated by the least-squares method. We can determine the residual matrix *ym*,*n*(*i*, *j*) with

$$g\_{m,n}(i,j) = G\_{m,n}(i,j) - G\_{m,n}(i,j). \tag{3}$$

Step 4. The detrended fluctuation *F*(*m*, *n*,*s*) for each subdomain *Xm*,*<sup>n</sup>* can be defined via the variance of *ym*,*n*(*i*, *j*) as follows

$$F^2(m, n, s) = \frac{1}{s^2} \sum\_{i=1}^s \sum\_{j=1}^s y^2\_{m, n}(i, j). \tag{4}$$

Step 5. The *q* − *th* order fluctuation is obtained by averaging over all the subdomain

$$F\_{\emptyset}(s) = \exp\left\{ \frac{1}{M\_s N\_s} \sum\_{m=1}^{M\_s} \sum\_{n=1}^{N\_s} \ln[F(m, n, s)] \right\}, q = 0 \tag{5}$$

$$F\_q(s) = [\frac{1}{M\_s N\_s} \sum\_{m=1}^{M\_s} \sum\_{n=1}^{N\_s} [F(m, n, s)]^q]^{1/q}, q \neq 0. \tag{6}$$

Step 6. The scaling relation of the fluctuation can be determined by analyzing the log-log *Fq*(*s*) versus the *s* for different values of *s* ranging from 6 to (*M*, *N*) /4, which reads

$$F\_q(\mathbf{s}) \propto \mathbf{s}^{h(q)}.\tag{7}$$

The scaling exponent *h*(*q*) can be obtained by the linear regression of ln *Fq*(*s*) to ln *s*, which is also called the generalized Hurst index. For each *q*, the corresponding traditional scaling exponent as *τ*(*q*)

$$
\pi(q) = qh(q) - D\_f. \tag{8}
$$

Note that, *Df* represents the fractal dimension of the geometric support. For the two-dimensional microscopic image of this paper, we take the value of *Df* = 2.

Step 7. The multifractal surface can be characterized by *Holder* ¨ exponent *α*(*q*) and singularity spectrum *f*(*α*), which are given by the Legendre transform [30].

$$a(q) = \pi'(q) = h(q) + qh'(q) \tag{9}$$

$$f(a) = qa(q) - \pi(q) = q[a - h(q)] + 2.\tag{10}$$

Here, the multifractal singularity spectrum *f*(*α*) is a continuous exponential spectrum used to characterize multiple fractal sets, which provides a complete statistical description of the internal inconsistencies of fractals.
