## 2 Materials and Methods
### Two-Dimensional Fractional Brownian Motion  
Falconer [7] called the process of 2D FBM an index-\(H\) Brownian function; [PERSON] et al. [30] called it isotropic 2D FBM, whose corresponding discrete version was called isotropic 2-dimensional discrete FBM; however, [PERSON] and [PERSON] [31] called it as a 2-variable FBM, whose corresponding discrete version was called a 2-variable discrete fractional Brownian surface. For simplicity, we call it 2D FBM, and an image of 2D FBM is a 2D FBM image. Suppose that an image of 2D FBM of size \(M\times N\) is denoted as  
\[I_{B}= \begin{bmatrix}B[0,0]&B[0,1]&\cdots&B[0,N-1]\\ B[1,0]&B[1,1]&\cdots&B[1,N-1]\\ \vdots&\vdots&\ddots&\vdots\\ B[M-1,0]&B[M-1,1]&\cdots&B[M-1,N-1]\end{bmatrix}  \tag{1}\]For 2 coordinates or points \(B[x_{1},\ y_{1}]\) and \(B[x_{2},\ y_{2}]\), their ACF is as follows:  
\[ \begin{split}& r_{BB}\bigg{[}\begin{pmatrix}x_{1}\\ y_{1}\end{pmatrix},\ \begin{pmatrix}x_{2}\\ y_{2}\end{pmatrix}\bigg{]}=E\{B[x_{1},\ y_{1}]B[x_{2},\ y_{2}]\}\\ &=\frac{\sigma^{2}}{2}\bigg{\|}\begin{pmatrix}x_{1}\\ y_{1}\end{pmatrix}\bigg{\|}^{2H}+\bigg{\|}\begin{pmatrix}x_{2}\\ y_{2}\end{pmatrix}\bigg{\|}^{2H}-\bigg{\|}\begin{pmatrix}x_{1}-x_{2}\\ y_{1}-y_{2}\end{pmatrix}\bigg{\|}^{2H},\end{split}  \tag{2}\]  
where \(H\) is the Hurst exponent. Obviously, the ACF is directly related to the distance between 2 coordinates as well as 2 respective distances from the origin, and hence it is non-stationary. For estimation, stationary ACFs are needed. Two approaches were proposed. The first was proposed by [PERSON] et al. [30], who used the second increment process of 2D FBM to obtain a stationary covariance matrix. The resulting process is called 2-dimensional Gaussian noise (FGN2) [32] or simply called the first 2D FGN. The second was proposed by [PERSON] and [PERSON] [31], who introduced another 2D FGN, simply called the second 2D FGN. [PERSON] [22] detailed these 2 2D FGNs and their related stationary covariances and symbols. Through the MLE, the mean-squared errors of the second 2D FGN are smaller than those of the first 2D FGN, but the first 2D FGN gives us more physical meanings because real surfaces or images may exist in the patterns of the first 2D FGN instead of the second 2D FGN. Therefore, in the paper, we adopted the first 2D FGN to establish a probability density function (PDF) in order to estimate the Hurst exponent via the MLE. To obtain this 2D FGN, we first represent \(I_{B}\) as a column vector of size \(MN\times 1\) as follows:  
\[\overset{\rightarrow}{\boldsymbol{B}}^{T}= \begin{bmatrix}\boldsymbol{B}_{0}^ {T}&\boldsymbol{B}_{1}^{T}&\cdots&\boldsymbol{B}_{M-1}^{T}\end{bmatrix} , \tag{3}\]  
where  
\[\boldsymbol{B}_{k}^{T}= \begin{bmatrix}B[k,0]&B[k,1]&\cdots&B[k,N-1]\end{bmatrix} ,\ k=0,\ 1,\ldots,\ M-1.\]  
Hence, its nonstationary covariance matrix is as follows:  
\[\boldsymbol{R}_{\overset{\rightarrow}{\boldsymbol{B}}\overset{\rightarrow} {\boldsymbol{B}}}=E\bigg{\{}\overset{\rightarrow}{\boldsymbol{B}}\overset{ \rightarrow}{\boldsymbol{B}}^{T}\bigg{\}} \tag{4}\]  
Next, [PERSON] et al. [30] introduced the following second increment process of 2D FBM  
\[ \begin{split} X[i,j]&=B[i,j]-B[i,j-1]-B[i-1,j]+B[i-1,j-1],\\ &\ i=1,\ 2,\ldots,\ M-1;\ j=1,\ 2,\ldots,\ N-1\end{split}  \tag{5}\]  
Likewise, we represent this 2D FGN image as a column vector of size \(N_{1}\times 1\), \(N_{1}=(M-1)(N-1)\):  
\[\overset{\rightarrow}{\boldsymbol{X}}^{T}= \begin{bmatrix}\boldsymbol{X}_{1}^ {T}&\boldsymbol{X}_{2}^{T}&\cdots&\boldsymbol{X}_{M-1}^{T}\end{bmatrix} , \tag{6}\]  
where  
\[\boldsymbol{X}_{k}^{T}= \begin{bmatrix}X[k,1]&X[k,2]&\cdots&X[k,N-1]\end{bmatrix} ,\ k=1,\ 2,\ldots,\ M-1.\]  
From the inner structure of the vector \(\overset{\rightarrow}{\boldsymbol{X}}^{T}\), it can also be viewed as a vector time series; that is, a time series consists of vectors.