## 2 Methodology
### Surface Generation  
Once a cleaned point cloud has been obtained, the next step is to determine the optimal plane for the triangulation, which is obtained by solving an eigenvalue problem to determine the best fit plane through the dataset (e.g., ([PERSON] _et al._, 2007)). The plane is defined by a normal vector (\(N=\{N_{x}\), \(N_{y}\), \(N_{z}\}\)) and the centroid coordinates of the original dataset (\(P_{c}\)). (Note that for computational efficiency, one can use a fraction of the points in the dataset to estimate the normal of the plane, if desired).  
Next, the necessary transformation to align the best fit plane with the XY plane can be represented as rotation about a vector or a quaternion ([PERSON] _et al._, 2007; [PERSON] _et al._, 2007). The orientation is established such that the Z\({}^{\prime}\) axis is positive away from the rock slope so that erosion values will be negative. The vector \(\rho\) is first calculated by the cross product:  
\[\rho=\overline{N}\times\overline{\tau}z \tag{2}\]where \(\tau_{z}=\{\)0,0,1\(\}\), represents a normalized vector for the Z axis (_i.e._, normal of the XY plane). The angle of rotation (\(\alpha\)) can be found with the cross product, which can be simplified with \(\tau\) as:  
\[\alpha=\sin^{-1}\left(\left|\rho\right|\right) \tag{3}\]  
The resulting rotation vector and angle can then be represented as a quaternion [46] producing a rotation matrix, \(R\), such that the transformed point coordinates can be calculated by:  
\[\vec{p_{i}}^{\prime}=R\vec{P}_{i}-\vec{P}_{C} \tag{4}\]  
Next, the minimum and maximum values for all points (\(p^{\prime}\)) are found. The number of columns and rows for the grid structure with a cell size dimension of \(\sigma\) are calculated as:  
\[n_{cols}=\left[\frac{\left|max(p^{\prime}_{\chi})\right|-\left|min(p^{\prime} _{\chi})\right|}{\sigma}\right],\ n_{rows}=\left[\frac{\left|max(p^{\prime}_{y })\right|-\left|min(p^{\prime}_{y})\right|}{\sigma}\right] \tag{5}\]  
Each point (\(p^{\prime}\)) is then subdivided into the square grid cells based on a spatial index _IND_:  
\[IND_{i}=id_{x_{i}}+id_{y_{i}}\times n_{cols} \tag{6}\]  
where:  
\[id_{x_{i}}=\left|\frac{p^{\prime}_{x_{i}}-\left|min(p^{\prime}_{x})\right|}{ \sigma}\right|,\ id_{y_{i}}=\left|\frac{p^{\prime}_{y_{i}}-\left|min(p^{\prime }_{y})\right|}{\sigma}\right| \tag{7}\]  
Then, the centroid coordinates for each cell (\(j\)) are calculated from all \(m\) points within the cell:  
\[p^{\prime}_{c}{}_{j}=\frac{1}{m}\sum_{k=1}^{m}p^{\prime}_{k} \tag{8}\]  
Once centroids have been found for all cells, a hole filling operation can be completed to create interpolated points in cells with insufficient data. The specifics of the hole filling approach will be described in more detail in the next section. Following hole filling, the triangulation is then formed using the rules established in [PERSON] _et al._[40] (Figure 4). In this approach, the positions of the centroids of four adjacent cells are evaluated to determine the appropriate connections to avoid intersections of triangles. Although there are five potential scenarios, only two triangulation patterns are needed, which can quickly be determined with a point in triangle test. This method proves to be very efficient because the cell centroids are structured and simple to connect with these rules.  
Figure 4: Triangulation scheme adopted from [PERSON] _et al._[40].  
After triangulation, the centroids can then be transformed back to the original coordinate system:  
\[\vec{P}_{\textit{c}}=R^{-1}\vec{p}_{\textit{i}}{}^{{}^{\prime}}+P_{\textit{c}} \tag{9}\]  
Similar to averaging the X, Y, Z coordinates to determine the centroid coordinates, RGB color or Intensity values can be averaged for all points in the cell for texture mapping and further analysis.