## 3 Materials and Methods
### WAI Reconstruction Algorithm Based on Finite Difference
#### 3.2.2 Reconstruction of the WAI  
Assume that \(h\)\((x,y)\) is the height of an arbitrary sampled point \((x,y)\) of the water surface, and \(\mathbf{Z}\) is a 2D height field. The gradient of \(\mathbf{Z}\) can be given by  
\[\
abla\mathbf{Z}=\frac{\partial\;h\;(x,y)}{\partial\;x}\;\hat{\mathbf{i}}+ \frac{\partial\;h\;(x,y)}{\partial\;y}\;\hat{\mathbf{j}}=\mathbf{Z}_{x}(x,y) \;\hat{\mathbf{i}}+\mathbf{Z}_{y}(x,y)\;\hat{\mathbf{j}}, \tag{8}\]  
where \(\mathbf{Z}_{x}(x,y)\), \(\mathbf{Z}_{y}(x,y)\) are the \(x\) and \(y\) components of the two-dimensional numerical gradient of the arbitrary sampled point \((x,y)\), respectively. \(\hat{\mathbf{i}}\) and \(\hat{\mathbf{j}}\) are the unit vectors in the \(x\), \(y\) axis directions, respectively. Therefore, we can obtain the first-order partial differential equation as follows:  
\[\begin{cases}\frac{\partial\;h\;(x,y)}{\partial x}=\mathbf{Z}_{x}(x,y)\\ \frac{\partial\;h\;(x,y)}{\partial y}=\mathbf{Z}_{y}(x,y)\end{cases}, \tag{9}\]  
The sampled normals \(\{\hat{\mathbf{N}}_{k}\}\) corresponding to the sampled points \(\{\mathbf{q}_{k}\}\), estimated in Section 3.2.1, are known. Moreover, the normal vector of an arbitrary sampled point on the wave surface can be expressed as \(\left[\mathbf{Z}_{x}(x,y)\;,\;\mathbf{Z}_{y}(x,y)\;,\;-1\right]\). For the gradient operator of the 2D discrete function \(h\)\((x,y)\), according to the finite difference theory [36], we use the central difference formula to approximate the first derivative. Therefore, Equation (9) becomes  
\[\begin{cases}\frac{\partial\;h\;(x,y)}{\partial x}=\frac{h\;(x_{i+1}y_{i})-h\; (x_{i-1}y_{i})}{x_{i+1}x_{i-1}}&i=1,2\ldots n;\;\Delta x_{i}=x_{i+1}-x_{i};\; \Delta x_{i-1}=x_{i}-x_{i-1}\\ \frac{\partial\;h\;(x,y)}{\partial y}=\frac{h\;(x_{i}y_{i+1})-h\;(x_{i}y_{i-1} )}{\Delta y_{i}+\Delta y_{i-1}}&j=1,2\ldots m;\;\Delta y_{i}=y_{i+1}-y_{i};\; \Delta y_{i-1}=y_{i}-y_{i-1}\end{cases}, \tag{10}\]where \(m\), \(n\) represent the dimensions of the grid of sampled points \(\{\hat{\mathbf{N}}_{k}\}\). Noting that \(\mathbf{H}\) is the vector of length \(\mathrm{mn}\) representing the height field \(h\)\((x,y)\) sampled on a \(m\times n\) grid, the two vectors of length \(\mathrm{mn}\) representing each component of the gradient field can be written as  
\[\begin{array}{c}\mathbf{G}_{x}\cdot\mathbf{H}=\mathbf{Z}_{x}\\ \mathbf{G}_{y}\cdot\mathbf{H}=\mathbf{Z}_{y}\end{array} \tag{11}\]  
where \(\mathbf{G}_{x}\) and \(\mathbf{G}_{y}\) are two sparse matrices, of size \(\mathrm{mn}\times\mathrm{mn}\), defining the linear combinations of the elements of \(\mathbf{H}\) to produce each gradient. Equation (11) can be merged into a single linear system,  
\[\mathbf{G}\cdot\mathbf{H}=[\mathbf{Z}_{x},\mathbf{Z}_{y}]^{\mathrm{T}}= \boldsymbol{\Xi}, \tag{12}\]  
where \(\mathbf{G}=[\mathbf{G}_{x},\mathbf{G}_{y}]^{\mathrm{T}}\) is a rectangular sparse matrix of size \(2\mathrm{mn}\times\mathrm{mn}\), and \(\boldsymbol{\Xi}\) is a vector of length \(2\mathrm{mn}\). This system thus gives \(2\mathrm{mn}\) equations with \(\mathrm{mn}\) unknowns. It is over-determined, so a direct inversion is not possible. However, an estimate of \(\mathbf{H}\) may be obtained, by minimizing the residual [37],  
\[\|\mathbf{G}\cdot\mathbf{H}\text{-}\boldsymbol{\Xi}\|^{2}\rightarrow\mathrm{ min}, \tag{13}\]  
where \(\|\cdot\|\) represents the Euclidean norm. The WAI shape estimated by the numerical integration of the inverse gradient operator is sparse and discontinuous, whereas the WAI is typically smooth and integrable. Thus, we further perform the bicubic interpolation algorithm [38] for \(\mathbf{H}\) to estimate the WAI shape.