Structure-from-Motion Approach for Characterization of Bioerosion Patterns Using UAV Imagery

3.1. Data Collection

The 3D model was achieved using images taken with a consumer-grade camera (Lumix DMC-TZ10, 12 MP, Panasonic, Osaka, Japan) mounted on a UAV. The flight vehicle is a hexacopter (DJI F550 air frame, DJI, Shenzhen, China) commanded by an open-source ArduPilot Mega 2.6 with a generic 2D gimbal (Figure 2a,b). The UAV includes an on-board GPS/IMU, logging at 5 Hz. Its total weight is about 3 kg. The estimated price, including camera and hexacopter (complete kit), is around US$1250.

To ensure high degree of overlap, flight speed was set at 1 m s⁻¹ and the images were taken at 3 m intervals, resulting in nearly 400 images. The UAV followed a predefined flight path, covering a set of flight lines (five in total) which are parallel to the vertical surface (Figure 2c). Each flight line was duplicated, acquiring two views on the same surface, with an angle of approximately 25° (Figure 2c). The total flight time took about 30 min. The field work was conducted in the early morning, when the illumination of the whole area was sufficiently homogeneous.

3.2. 3D Model Reconstruction

After carrying out the selection of images to be analyzed, the first step was to identify and match features automatically in consecutive images through an object recognition system called scale invariant feature transform, using the open-source software package VisualSFM [24]. The approach transforms an image into a large set of local feature vectors. The features are invariant to the image scaling and rotation, but they are partially invariant to changes in illumination conditions and 3D camera viewpoint [25]. Within this process, all unstable points were filtered. A second step was to estimate both intrinsic (geometric and optical characteristics) and extrinsic (3D position and orientation of the camera frame in a relative coordinate system) camera parameters by determining optimal camera positions using the bundle adjustment algorithm into the VisualSFM software. A sparse point cloud is built, representing the most important features (Figure 2c).

In order to generate a dense point cloud, Clustering Views for Multi-view Stereo (CMVS) and Patch-based Multi-view Stereo (PMVS) algorithms were applied using the VisualSFM software. The software provides a convenient menu interface to run external tools like the aforementioned algorithms. The previously determined camera positions were used as input. CMVS decomposes input images into clusters (sub-models) and PMVS reconstructs 3D data for each cluster [26].

3.3. Georeferencing

The resulting dense point cloud (11,364,917 points) is generated in an arbitrary (image) coordinate system, which must be transformed into an absolute coordinate system. To achieve this transformation, firstly, eight well-distributed anthropogenic benchmarks along the cliff (four posts in the upper part; three safety cones and one sign board in the lower part) were considered to aid the location in the point cloud. The 3D position of these points was performed using a Sokkia Radian IS real-time kinematic differential Global Positioning System (DGPS-RTK) unit (Sokkia, Mississauga, Canada) in order to obtain high accuracy (horizontal accuracy = 10 mm, vertical accuracy = 20 mm). The location is given in UTM coordinate system (zone 20S) that then was converted into manipulable data. Finally, georeferencing process of the 3D point cloud was carried out in VisualSFM software. This software allows to combine all the resultant sub-models into a single, larger-scale model and, once georeferenced, permits to save the model in PLY format (Stanford Triangle Format).

3.4. Testing the Model

With the aim of testing the model, we collected 3D coordinates for a large number of points over a portion of the vertical surface by using a low cost do-it-yourself (DIY) terrestrial laser scanning (TLS). Therefore, both SFM photogrammetry and DIY-TLS point cloud datasets were finely aligned and then compared in terms of absolute distance by means of the open-source software CloudCompare [27].

The DIY-TLS (previously tested) consists of a commercial laser range finder (model UT390B, Uni-T, Shanghai, China; accuracy: ±2 mm) that rotates around two axes (x, y) using stepper motors; the equipment is mounted on a tripod. The sampling rate is 1 point every 1.1 s. A program was developed to control scanning motion and laser pulses using the Arduino open-source electronics prototyping platform (Mega 2560, Arduino, Torino, Italy).

Also, a radio telemetry system (3DR radiotelemetry Kit-915 MHz, 3DR, San Diego, CA, USA) with a range up to 500 m was used. Data extraction from the DIY-TLS is based on time-of-flight of the laser pulse to travel from its source to the target object and back; therefore, it computes the distance based on the travel speed of the pulse [28]. Coordinates of each point are acquired in a polar coordinate system (ρ, ϑ, φ). Subsequently, these coordinates were transformed into the Cartesian system (x, y, z) and corrected based on known points of reference (DGPS-RTK).

3.5. Calculations Derived from the Model

Three properties (for each vertex) were stored in PLY file format: spatial coordinates (x, y, z), normal (nx, ny, nz) and color (RGB). To carry out post-run analysis, the georeferenced model (PLY format) was saved in ASCII format using the CloudCompare software. Then, all the header field values in the text file (“.txt”) were removed. The analysis was run under the GNU Octave software [29]. In general terms, two sets of calculations were derived from the model: one corresponding to the topographic properties of the vertical surface, and the other corresponding to the detection of the BF. Finally, associations between the results were carried out using bivariate statistical analysis.

4.1. 3D Point Cloud Model

A perspective view of the dense point cloud model covering the total area is presented in Figure 3a. As was mentioned before, the sparse point cloud comprised 323,847 points, while the dense point cloud produced 11,364,917 points. A map of local point density of the 3D model corresponding to the study area (i.e., 3,714,511 points) is shown in Figure 3c. The map was created from a sphere of radius 0.5 m by using the CloudCompare software. Clearly, the point density is more concentrated over the most prominent portion of the vertical surface, due to the fact that the point density is inversely proportional to the distance from the equipment. In this latter case, the values of density were greater than 12,500. Low values of point density (from 1000 to 2500) in specific portions of the study area are attributed to the presence of complex topography, which might cause an inefficient data capture (Figure 3c). In general terms, although there is a slightly lower local point density inside of the BF compared to the rest of the surface (Figure 3d), it is noticeable that these points tend to lose detail in y spatial coordinate from a certain depth.

Also, a histogram of the local point density is presented in Figure 3e. It showed a normal distribution, with a mean of 9749 and a standard deviation of ±3237. The local point density ranged from 5 to 20,015 points. Very low values of density at the lower end is due to occasional presence of vegetation cover.

4.2. Test of the Model

The test of the model was carried out on a surface of about 6.7 m² containing a complex topography. A comparison of the two georeferenced point cloud datasets, that is, the SFM photogrammetry model (50,944 points) and the DIY-TLS (1761 points), was made (Figure 4a). The comparison consisted of measuring the absolute distance between each point in the compared dataset with its closest point (i.e., nearest neighbor distance). In general terms, a good accuracy between the two datasets is noticeable (Figure 4b). The absolute distance variability appears not to have a clear pattern in the tested surface.

The Table 1 contains the statistical calculations of the absolute distance between the two dataset for the three spatial coordinates. Results concerning x, y, z coordinates jointly, indicated a mean distance of ∼0.07 m with a standard deviation of ±0.035 m. The y coordinate presented the greatest mean distance, with a value of 0.04 m. The remaining two coordinates showed similar values of the order of 0.03 m. It must be taken into consideration that the mean absolute (three-dimensional) distance is within the limits of the root-mean square error (0.19 m) involved in the georeferencing process of the whole SFM model.

4.3. Topographic Properties of the Vertical Surface

As a result of the achieved level of density, it was possible to obtain a high resolution model (0.013 m) for resolving fine spatial details in topography. Figure 5a–d shows all the topographic properties that were considered in this study. Calculations of slope and surface roughness properties give a local measure, independent of small-scale topographic variability, having the potential to differentiate individual surface properties to which it is most sensitive. Our results show that both local slope and surface roughness (e.g., Laplacian operator method) properties, which ranged from 0.1° to 89.3° and from 198 to 3364, respectively (Figure 5b,d), were able to demonstrate topographic patterns visually consistent with the real environment at a high level of detail.

The mean slope angle of the study area is 47°; the highest values were found in several parts of the cliff, but more significantly in the top portion (Figure 5b). In the case of surface roughness, whichever method is used, it also increases near the top of the cliff (Figure 5c,d). The roughness based on the Laplacian operator method showed a mean variability of 22%. According to the map of roughness (Figure 5d), warm colors denotes high variability, showing convex and/or concave curvature of the surface. An important pattern to note is, for example, the presence of a narrow convex-concave surface over the upper portion of the cliff.

4.4. Bioerosion Pattern Characterization

As a result of processing different kinds of information provided through the use of SFM approach, it was possible to detect BF, and then to characterize their geometry. In this work, an innovative method for the detection of the BF was implemented, which is given by the addition of color data to spatial coordinates (particularly the y coordinate). A total of 858 BF were detected in the study area (i.e., vertical surface = 480 m²) (Figure 5e). These features comprise a bioeroded surface of the order of 5.7%, and a frequency of occurrence of 1.8 cases per square meter. However, the BF are concentrated close to the top of the cliff, at a mean height of about 6.5 m above the ground (Figure 5e).

The mean area of the BF is 320 cm² with a similar value of standard deviation, which is indicating a great variability in their size (Figure 6a). The frequency distribution of this variable follows a lognormal distribution as is shown in Figure 6a. There is a predominance of BF ranging from 50 to 200 cm². Regarding their perimeter, it reached a mean value of 76 cm with a relatively high standard deviation (±47 cm). The frequency distribution of the perimeter also follows a lognormal distribution (Figure 6b). Circularity index showed a mean value of 0.67, which indicates BF shapes tending to be more irregular or elongated.

The standard deviation of this index was low (±0.18), implying that the variability was insignificant. The frequency distribution of the circularity index follows a normal curve (Figure 6c); also, it can be seen that the histogram shows a second and less significant peak close to values of 0.9.

There is a high positive (potential) correlation (R² = 0.97) between area and perimeter (Figure 7a), demonstrating that the variables perform equivalently, as is expected. Plotting the circularity index against either area or perimeter, it shows a negative slope, that is, decreasing circularity with increasing area or perimeter (Figure 7b,c). A relatively high (logarithmic) correlation of 0.62 and 0.77 was found between circularity and both area and perimeter, respectively; however, certain degree of correlation is unavoidable since the same data is used in both variables.

The study area was partitioned into 125 subareas based on equal sized cells of 4.8 × 0.8 m, in order to examine trends and interrelationships between BF and surface topography. According to the plots of Figure 8, we could explain the behaviour of the geometric (BF) and topographic (vertical surface) properties against the frequency of BF. It was possible to observe the predominance of subareas with low and, to a lesser extent, middle frequency of the BF per square meter. In the case of mean values of both area and perimeter, a clear trend to increase with an increase in the frequency of BF was observed, indicating that the largest BF are situated closer to each other (Figure 8a,b). In contrast, the standard deviation values of area and perimeter show a trend toward decrease with an increase in the frequency of BF, implying a decrease in the variability of the dimension according as they are more concentrated. The mean and standard deviation of the circularity index also presents a mild decrease with increasing frequency of BF (Figure 8c).

When increasing the frequency of BF, an increase of the mean slope angle values of the vertical surface is observed (Figure 8d). The latter indicates that the largest and most concentrated BP tend to be situated on steeper surfaces, ranging from 50° to 60° (mean values). Figure 8e,f show the surface roughness property calculated by using circular standard deviation and Laplacian operator methods, respectively. Regarding the former, there is no defined trend. However, when we use the Laplacian operator method, an apparent trend can be observed, indicating that the mean roughness (convex and/or concave curvature) increase with an increase in the frequency of BF as is expected.

4.5. Applications and Limitations of SFM-UAV

As proven in several studies [10,12,13] and in this work, the topography generated by SFM photogrammetry offers similar performance in resolution and precision to that of alternative approaches (e.g., TLS, aerial laser scanning (ALS), or real-time kinematic GPS) but with the advantage that it is much cheaper and easier to implement. In addition, UAV technology is fully suitable for SFM photogrammetry approach, due to its wide range of operational possibilities (e.g., position and speed control of flight) for capturing photos. Although the use of SFM-UAV technology implies costs of equipment (in this case, ∼US$1,250) and sometimes personnel to make the assessment, it is much cheaper than, for example, TLS or ALS, usually by at least one order of magnitude. These advantages are of great importance in multitemporal studies, in order to evaluate the spatial dynamics of a given surface.

Unlike laser scanning systems such as TLS/ALS, SFM photogrammetry provides at the same time 3D spatial (x, y, z) normal (nx, ny, nz) and color (RGB) datasets. As a result of this large amount of data, SFM photogrammetry has become more demanding regarding to post-processing tasks. Particularly, the availability of color data represents a clear advantage in terms of practicability, since by means of alternative methodologies this information is only possible to obtain by adding a digital camera mounted independently of the measurement device [33–35]. In this work, the combined use of both 3D spatial and color datasets have contributed to the determination of the BF.

Despite the aforementioned advantages, there are still limitations in the field of the photogrammetry referred to the sunlight exposure over topographically complex surfaces (e.g., cavities, depressions, etc.), which might be well resolved by laser scanning systems. This limitation causes pixel values near zero, difficulting the process of feature matching. In the current study, when considering the y spatial coordinate, the point cloud well inside the BF begins to generate erroneous data at a given depth where the sunlight never penetrates.