Integrating Three-Dimensional Benthic Habitat Characterization Techniques into Ecological Monitoring of Coral Reefs

Long-term ecological monitoring of reef fish populations often requires the simultaneous collection of data on benthic habitats in order to account for the effects of these variables on fish assemblage structure. Here, we described an approach to benthic surveys that uses photogrammetric techniques to facilitate the extraction of quantitative metrics for characterization of benthic habitats from the resulting three-dimensional (3D) reconstruction of coral reefs. Out of 92 sites surveyed in the Northwestern Hawaiian Islands, photographs from 85 sites achieved complete alignment and successfully produced 3D reconstructions and digital elevation models (DEMs). Habitat metrics extracted from the DEMs were generally correlated with one another, with the exception of curvature measures, indicating that complexity and curvature measures should be treated separately when quantifying the habitat structure. Fractal dimension D64, calculated by changing resolutions of the DEMs from 1 cm to 64 cm, had the best correlations with other habitat metrics. Fractal dimension was also less affected by changes in orientations of the models compared to surface complexity or slope. These results showed that fractal dimension can be used as a single measure of complexity for the characterization of coral reef habitats. Further investigations into metrics for 3D characterization of habitats should consider relevant spatial scales and focus on obtaining variables that can complement fractal dimension in the characterization of reef habitats.


Introduction
Long-term ecological monitoring of reef fish populations generally aims to evaluate the current status of the ecosystem and to detect any trend or sudden changes in the ecosystem due to natural and/or anthropogenic disturbances.One of challenges in conducting fish monitoring is that biotic and abiotic factors, such as water motion, depth and habitat, can affect the structure of fish assemblages [1].Benthic composition can influence the distributions of fishes that rely on benthic organisms for food [2].Measures of a habitat's structural complexity, such as numbers and sizes of holes and topographic relief, also affect the diversity and abundance of fish assemblages by altering the amount of available refuge space, capacity for niche specialization and the encounter rate between competitors [3][4][5][6][7].In coral reef environments, disturbances such as tropical storms [8], bleaching events [9,10] and the decline of herbivores [11] all have the potential to alter benthic cover and structural complexity, which in turn affect the structure of fish assemblages.Therefore, variations in the structure of fish assemblages explained by these variables need to be considered when assessing temporal changes in the fish assemblages [12].
In recent years, technological advancements in image sensors and data-processing capacity has led to the successful use of photogrammetric techniques (e.g., structure-from-motion photogrammetry) as a cost-effective means to generate three-dimensional (3D) reconstructions of coral reef habitats [13][14][15][16].The resulting 3D models can be used to quantify the habitat structure by extracting metrics (e.g., fractal dimension, vector dispersion, rugosity, surface complexity, slope and curvature) directly from the digital surface models [15] or from the digital elevation models (DEMs) exported from the 3D models [13].Habitat metrics extracted from 3D models of coral colonies have been, in turn, linked to reef fish assemblages [17].Three-dimensional reconstructions can also be used for volumetric analyses to track temporal changes in the habitat structure and/or for obtaining data on the distribution of benthic species after exporting them as orthophotomosaics [14,18,19].While the use of photogrammetric techniques requires more time for post-survey data processing than in situ observations, it potentially generates high-resolution data [14] with relatively high precision and accuracy [20], and it is less prone to observer bias [21].Image acquisition for 3D reconstruction is also less time consuming in the field than in situ (i.e., underwater) methods that are commonly used to characterize either the habitat structure (e.g., chain-and-tape rugosity and Habitat Assessment Scores [5]) or benthic cover [22].
Here, we describe a workflow for obtaining habitat variables through benthic surveys, using structure-from-motion photogrammetric techniques, and subsequent processing of the resulting DEMs.We applied these methods to a reef monitoring program in the Northwestern Hawaiian Islands (NWHI) to quantify the 3D habitat structure of coral reefs throughout this region.The NWHI consist of ten major islands and atolls that span approximately 2000 km (Figure 1).They are part of the Papahānaumokuākea Marine National Monument, a marine protected area encompassing 1.5 million square kilometers that is one of the largest conservation areas in the world.While monitoring of shallow-water (≤30 m) coral reef fish and reef habitats has been conducted annually in the NWHI during the reef assessment and monitoring program (RAMP) expedition each summer, maintaining a consistent and sufficient sample size (i.e., number of survey sites) at each island/atoll can be challenging due the large spatial area and remote location of the monument.Nonetheless, a recent study investigating temporal changes in the fish assemblages in the NWHI revealed the importance of obtaining habitat variables when conducting fish surveys [23].The purpose of this study was, therefore, to establish a procedure to integrate a benthic survey using photogrammetric techniques into the reef fish monitoring in the NWHI without substantial increases in the time required for divers to complete their visual surveys at each site.Owing to a plethora of habitat metrics that can be extracted from DEMs using various geographic information system (GIS) software/packages, identifying the appropriate and relevant metrics for habitat characterization is also important.It streamlines the post-expedition data processing and allows scientists and managers to obtain data quickly.However, application of the procedures described here is not limited to fish monitoring surveys in the NWHI and can be used to accompany any other underwater visual surveys that benefit from high-resolution habitat data.use of photogrammetric techniques requires more time for post-survey data processing than in situ observations, it potentially generates high-resolution data [14] with relatively high precision and accuracy [20], and it is less prone to observer bias [21].Image acquisition for 3D reconstruction is also less time consuming in the field than in situ (i.e., underwater) methods that are commonly used to characterize either the habitat structure (e.g., chain-and-tape rugosity and Habitat Assessment Scores [5]) or benthic cover [22].Here, we describe a workflow for obtaining habitat variables through benthic surveys, using structure-from-motion photogrammetric techniques, and subsequent processing of the resulting DEMs.We applied these methods to a reef monitoring program in the Northwestern Hawaiian Islands (NWHI) to quantify the 3D habitat structure of coral reefs throughout this region.The NWHI consist of ten major islands and atolls that span approximately 2000 km (Figure 1).They are part of the Papahānaumokuākea Marine National Monument, a marine protected area encompassing 1.5 million square kilometers that is one of the largest conservation areas in the world.While monitoring of shallow-water (≤30 m) coral reef fish and reef habitats has been conducted annually in the NWHI during the reef assessment and monitoring program (RAMP) expedition each summer, maintaining a consistent and sufficient sample size (i.e., number of survey sites) at each island/atoll can be challenging due the large spatial area and remote location of the monument.Nonetheless, a recent

Image Acquisition
Images of coral reefs were collected from islands/atolls of the NWHI during the 2017 RAMP expedition aboard NOAA ship Survey sites at each island or atoll were randomly chosen using the GIS software ArcMap (v.10.4) by generating random points (i.e., coordinates) on bathymetric and bottom composition maps, in which target survey sites were hard-bottom habitats at depths ≤30 m.
At each site, a pair of SCUBA divers first laid a 30-m transect line along a depth contour to survey the reef fish inside two adjacent cylinders that were 15 m in diameter using a stationary point count method [24].After the completion of the fish survey, the divers collected images along the transect line following a procedure that was similar to the one described in Burns et al. [13].A ground control point (GCP) unit was placed at one end of the transect line.The GCP unit, consisting of three polyvinyl reflectors, was used to obtain accurate x, y and z reference coordinates for scaling.After setting the white balance using a white slate, overlapping photographs were manually taken along the transect line using the single shooting setting, down on one side and back on the other side, while swimming slowly at approximately 3 m above the substrate aiming to cover a 90-m 2 (30-m length × 3-m width) area with 70-80% overlaps between images.While this simple approach of image acquisition from above allowed the divers to complete the process relatively quickly and, thus, to maintain the number of monitoring surveys that they could complete per day, certain structures such as overhangs associated with table corals would unlikely be captured in our 3D reconstructions.As our main focus was to characterize each reef habitat as a whole, and not the individual coral colonies, we considered this as an acceptable trade-off for maintaining the sample size for the monitoring program.The number of images captured at each site was, on average, 252 ± 50.48 SD (standard deviation).All photographs were taken using a Canon EOS Rebel SL1 digital SLR camera with an 18-55 mm lens in an Ikelite housing with a 6-inch dome port.The camera lens was set to the focal length of 18 mm with a shutter speed of 1/250, an aperture of f/10 and an auto ISO.

Generation of 3D Models
Three-dimensional models were constructed from the acquired imagery using the software Agisoft PhotoScan Professional (v.1.4), following the procedure of Burns et al. [13] (Table 1).Briefly, a sparse 3D point cloud was generated through the photo alignment process for each survey site.The key point limit and tie point limit were initially set to 50,000 and 5000, respectively, with the generic preselection enabled.However, if all the photos did not align with these initial settings, they were increased to 70,000 and 8000, respectively, with the generic preselection disabled.The x, y and z values from the GCP unit were used as known reference coordinates to scale each model accurately.After optimizing the camera alignment using the reference coordinates, a dense point cloud was generated and, in turn, used to generate a triangulated mesh model.A spatially rectified orthophotomosaic and DEM were then exported from the Agisoft PhotoScan at a 1-cm resolution as GeoTIFF elevation data files with local coordinates for further processing.The process was completed on either an Intel Xeon workstation with a 32GB RAM and FirePro W5100 GPU or Intel Core i7 laptops with a 64GB RAM and GeForce GTX 980 GPU.

Process Settings
Align Photos High accuracy, generic preselection enabled, 50,000 key point limit, 5000 tie point limit.
Optimize Camera Alignment (Use all the ones selected by the software.) Build Dense Cloud Medium quality, mild depth filtering, reuse depth maps disabled.
Build Mesh Arbitrary surface type, high face count, interpolation enabled, calculate vertex colors enabled.
Build Texture Adaptive orthophoto mapping mode, mosaic blending mode, texture size/count 16,384, enable hole filling.

Quantification of the Habitat Structure
Quantification of the habitat structure was done using the software ArcMap v.10.4 with 3D Analyst and Spatial Analyst licenses and Benthic Terrain Modeler (BTM) v.3.0 [25], as well as the statistical software R v.3.4.4 [26] with the sp [27,28], raster [29] and rgeos [30] packages.The DEMs were imported into the programs as raster files for spatial analyses.All habitat metrics obtained in ArcMap and R (Table 2, also see Sections 2.3.1 and 2.3.2 for details) were examined for their relationships using Pearson correlation coefficients (ρ) obtained in R.

ArcMap Procedure
Each DEM and orthophotomosaic were imported as data layers into the ArcMap software.Linear rugosity, surface complexity and slope were obtained using the functional surface toolset available with the 3D Analyst license.For linear rugosity, five polylines (straight lines, each approximately 20-m in length) were haphazardly drawn on each DEM avoiding the edge of the DEM where the amount of image overlap tended to be less than that in the central area.ArcMap's add surface information tool was used to add surface information (e.g., shape distance, surface distance and the minimum, maximum and average slope) to each polyline and the ratio of the surface distance and shape distance was calculated.The mean of the ratios from the five polylines was used as a measure of linear rugosity for each DEM.For surface complexity and slope, five rectangular polygons of an approximately 8-m 2 area (for a total of 40-m 2 area) were haphazardly digitized on a DEM avoiding the edge of the DEM, and surface information was added to each polygon.Surface complexity for the DEM was calculated as the ratio of the total 3D surface area of the five polygons to the total 2D surface areas of the five polygons.Slope was calculated as the mean of the average slope values in the attribute table of the polygons.
Three-dimensional topographic features, including curvature and viewshed, were calculated using the surface toolset available with the Spatial Analyst license.Each DEM was first clipped with the five polygons that were used to obtain surface complexity and slope.Curvature values (planform curvature and profile curvature) were obtained with the curvature tool using the clipped DEM as the input raster and taking the mean values from the layer properties of two output raster layers.Viewshed values (viewshed visible area and non-visible area) were also obtained with the viewshed tool by using the clipped DEM as the input raster.Five observer points were digitized on the centers of each of the five polygons and then used to calculate the spatial area that could be seen from the observer points (visible) and the area that could not be seen from any points (non-visible).
The Benthic Terrain Modeler [25] was used to obtain the slope, terrain ruggedness (VRM) and the surface area to planar area ratio.For each of the three BTM variables, the clipped DEM was used as an input bathymetry raster and the mean value was taken from the layer property of the output raster created by the respective BTM function.

R Procedure
The aggregate function in the raster package was utilized to remove the edge of a DEM by aggregating the DEM to either a 64-or 128-cm resolution.Specifically, the aggregate function with an aggregation factor n would produce a cell containing a null value if any of the n × n cells being aggregated were outside the extremity of the DEM.After removing these cells with null values from the aggregated raster layer and converting the raster to a polygon, it was used to clip the original DEM (1-cm resolution), thereby removing the edge.As clipping DEMs with polygons created at a 128-cm aggregation generally removed too much area (see the Result section for details), the surface complexity, slope and curvature were calculated using DEMs clipped with polygons created at a 64-cm aggregation (Text S1).
Surface complexity was obtained as the ratio of the 3D surface area along reef contours to the 2D surface (planar) area.The three-dimensional surface area was calculated using the surfaceArea function in the sp package.Slope and curvature were calculated using the surface properties of the raster cells in 3 × 3 windows, as in the ArcMap's surface toolset available with the Spatial Analyst license (Text S2, also see the ArcMap tool Reference [31] and the References therein [32,33] for details).For measures of curvature, the profile curvature and planform curvature were calculated as described in Zevenbergen and Throne [33].
Fractal dimension has been previously suggested as a matric for quantifying the shapes of corals or structural complexity of coral reefs [15,[34][35][36].Fractal dimension (D) was calculated as D = 2 − slope logS(δ)/log(δ), in which δ is a resolution of the raster and S(δ) is the 3D surface area at the given resolution [15].Therefore, calculating the fractal dimension required changing the resolution of the DEMs using the aggregate function.Two separate calculations were done for each DEM, one changing the resolutions from 1 cm to 2, 4, 8, 16, 32 and 64 cm (D 64 ), and the other changing the resolutions from 1 cm to 2, 4, 8, 16, 32, 64 and 128 cm (D 128 ).The two maximum aggregations of 64 cm and 128 cm were chosen based on preliminary investigations where we examined the decreases in 3D surface areas with increases in the aggregation factors (i.e., decrease in resolutions).In theory, fractal dimension remains the same across different spatial scales as the relationship between log(δ) and logS(δ) is linear [34].(i.e., The slope of log(δ) and logS(δ) is a constant.)In our case where we investigate habitat metrics within a limited spatial area, working from small to large scales, however, the 3D surface area starts approaching an asymptote as it approaches the 2D planar area because the 3D surface area cannot be smaller than the 2D planar area.The maximum aggregations, thus, needed to be chosen to minimize the effect of the change in the slope.DEMs clipped with polygons created at a 64-cm aggregation and a 128-cm aggregation were used to obtain D 64 and D 128 , respectively.

Results
In total, the divers collected overlapping photographs from 92 sites in the NWHI, with images from 85 sites achieving complete alignment and producing 3D models and DEMs for analyses of the habitat structure.Ground sampling distances (resolution/pixel) of the 3D models generated in Agisoft Photoscan ranged from 0.00026 to 0.00084 m/pixel, with errors of 0.5-2.3pixels, confirming that the designated resolution of 0.01 m (1 cm) for DEMs was within the range of model accuracy.The coverage areas of the 85 DEMs ranged from 57.5 to 284.2 m 2 , with a mean coverage of 142.3 m 2 (Table 3); 78% of the DEMs had coverage between 100 and 200 m 2 .Imagery from the seven sites that did not have complete alignment contained photographs that had a relatively low overlap (visual estimation of less than 50%) with others.These alignment issues produced gaps in the resulting 3D point clouds, so these sites were not included in the present study.
DEMs clipped with polygons to remove the edges retained, on average, 76% of the planar areas of the original DEMs for the 64-cm aggregation, and 54% of the planar areas for the 128-cm aggregation (Table 3).For DEMs clipped with a 64-cm aggregation, retention of the planar areas ranged from 54% to 86%, with 14 DEMs retaining less than 70%.For DEMs clipped with a 128-cm aggregation, retention of the planar areas ranged from 19% to 73%, with 54 DEMs retaining less than 60% of the original planar areas.
Table 3. Results of clipping digital elevation models (DEMs) with polygons, showing the mean, SD (standard deviation), median, minimum and maximum values of the planar areas of the original DEMs (original planar area), planar areas after the 64-cm aggregation (Planar area ×64), planar areas after the 128-cm aggregation (Planar area ×128), percent retention of planar areas after the 64-cm aggregation (area retained ×64), and percent retention of planar areas after the 128-cm aggregation (area retained ×128).Strong correlations were detected between many of the habitat metrics quantified using ArcMap and R (Figure 2).For example, the linear rugosity obtained using the ArcMap's functional surface toolset had strong correlations (|ρ| ≥ 0.75) with the surface complexity obtained using the functional surface toolset and viewshed (non-visible) obtained using the surface toolset, as well as fractal dimension, surface complexity and slope calculated in R. Similarly, the VRM and the surface area to planar area ratio obtained using the BTM in ArcMap were strongly correlated with each other and they also had relatively strong correlations (|ρ| ≥ 0.65) with the fractal dimension, surface complexity and slope calculated in R and the viewshed obtained using the surface toolset.Planform curvature and profile curvature had the correlation of ≈ 1, whether they were obtained in ArcMap or R, and the curvatures obtained in ArcMap had strong negative correlations (ρ = −0.92)with curvatures obtained in R.These curvature measures did not have strong correlations with other complexity variables, with the exception of the slope obtained using the functional surface toolset in ArcMap.
Values of the fractal dimension ranged from 2.029 to 2.249 for D 64 and 2.026 to 2.230 for D 128 , with larger values corresponding to higher structural complexity of the habitat (Figure 3).The generally lower values of D 128 in comparison to D 64 were due to the inclusion of data at the 128-cm aggregation, resulting in a gentler slope in some cases as values of the 3D surface areas approached the values of their 2D planar areas (e.g., FFS-4297 and LIS-4199 in Figure 3).However, the two fractal dimensions, D 64 and D 128 , had a high correlation of ρ = 0.99.They also had relatively strong correlations (|ρ| ≥ 0.7) with all other habitat metrics, with the exception of the curvature measures and slope obtained using the functional surface toolset.While D 128 had slightly higher correlations with rugosity, surface complexity and slope than D 64 , D 64 had higher correlations with the VRM and the surface area to planar area ratio obtained using the BTM than D 128 . in R.These curvature measures did not have strong correlations with other complexity variables, with the exception of the slope obtained using the functional surface toolset in ArcMap.
Values of the fractal dimension ranged from 2.029 to 2.249 for D64 and 2.026 to 2.230 for with larger values corresponding to higher structural complexity of the habitat (Figure 3).The generally lower values of D128 in comparison to D64 were due to the inclusion of data at the 128-cm aggregation, resulting in a gentler slope in some cases as values of the 3D surface areas approached the values of their 2D planar areas (e.g., FFS-4297 and LIS-4199 in Figure 3).However, the two fractal dimensions, D64 and D128, had a high correlation of ρ = 0.99.They also had relatively strong correlations (|ρ| ≥ 0.7) with all other habitat metrics, with the exception of the curvature measures and slope obtained using the functional surface toolset.While D128 had slightly higher correlations with rugosity, surface complexity and slope than D64, D64 had higher correlations with the VRM and the surface area to planar area ratio obtained using the BTM than D128.

Discussion
This study utilized photogrammetric techniques to create 3D reconstructions of coral reef habitats in the NWHI.Multiple metrics of habitat structure were extracted from the resulting DEMs to establish a procedure for time-efficient characterization of benthic habitats that can accompany the Monument's reef fish monitoring survey.The amount of time required for divers to complete image acquisition at each site was approximately 10 minutes, although the time requirement varied slightly depending on the presence/absence of current and surge, with surge being more disruptive than current.Given the amount of information we could extract from the resulting 3D models (e.g., habitat metrics as described in the present study, distribution of benthic species and 3D volume), we considered the 10-minute timeframe cost-effective, although the exact amount of time required to

Discussion
This study utilized photogrammetric techniques to create 3D reconstructions of coral reef habitats in the NWHI.Multiple metrics of habitat structure were extracted from the resulting DEMs to establish a procedure for time-efficient characterization of benthic habitats that can accompany the Monument's reef fish monitoring survey.The amount of time required for divers to complete image acquisition at each site was approximately 10 min, although the time requirement varied slightly depending on the presence/absence of current and surge, with surge being more disruptive than current.Given the amount of information we could extract from the resulting 3D models (e.g., habitat metrics as described in the present study, distribution of benthic species and 3D volume), we considered the 10-min timeframe cost-effective, although the exact amount of time required to obtain the same information in situ would depend on individual surveyors and habitat types.Nonetheless, the relatively quick timeframe required to complete the image acquisition provided the divers with ample time to complete all fish surveys, highlighting that the workflow for 3D mapping described here is a viable method to be integrated into reef fish surveys during annual monitoring expeditions.
The metrics for the 3D habitat structure obtained in the present study generally correlated well with one another, with the exception of the curvature measures, indicating that complexity and curvature measures should be treated separately when quantifying benthic habitat.Correlations between viewshed and other measures of structural complexity were somewhat unexpected but explicable, as the amount of area that is visible or not visible from observer points would be affected by terrain heterogeneity.Among all the metrics examined, fractal dimension, especially D 64 , had the best correlations with other measures of habitat complexity computed using both R and ArcMap.Unlike all other habitat metrics calculated using DEMs at a 1-cm resolution, fractal dimension combined information obtained at various spatial scales (i.e., 3D surface areas of DEMs at different resolutions in our case), allowing for the quantitative characterization of the habitat structure as a whole [36].While the calculation of fractal dimension required changing resolutions of each DEM from 1 to 2, 4, 8, 16, 32, 64 and 128 cm, this was achieved relatively quickly in R using the aggregate function in the raster package.Overall, processing DEMs in R required time to write all the scripts, but once they were written, it automated the whole process, from importing image files to changing resolutions, obtaining different habitat metrics and exporting data as a data frame all at once.In ArcMap, the same analysis would have required, for each survey site, manually processing each of the DEMs at different resolutions to obtain 3D surface areas and manually compiling the data from different resolutions to calculate fractal dimension.This advantage of speed and automation, as well as the high correlations between fractal dimension obtained in R and the metrics obtained in ArcMap, demonstrated the high efficiency of the fractal analysis in R relative to the manual approach in ArcMap.
Measures of rugosity (i.e., linear rugosity and 3D surface complexity) had a correlation of >0.85 when either the functional surface toolset in ArcMap or R was used.On the other hand, the surface area to planar area ratio and the VRM obtained using the BTM had much lower correlations with these values.Although these metrics obtained in the BTM measured terrain heterogeneity, they were calculated for each raster cell using the 3 × 3 neighborhood [25], instead of examining the total 3D distance/area versus the total planar distance/area of a DEM.Similarly, the slope calculated using the BTM (BTM slope) and R (R slope) had a strong correlation of 0.85, but they both had correlations of ~0.5 with the slope calculated using the functional surface toolset in ArcMap.This was also likely due to differences in how these values were calculated; BTM slope and R slope were calculated in the same way using a moving 3 × 3 window, while the slope from the functional surface toolset was measured in grade (the tangent of the angle of a surface).It is also interesting to note that the slope in the functional surface toolset had stronger correlations with measures of rugosity, as well as the planform and profile curvatures, than with either the BTM slope or R slope (Figure 2).This emphasized the importance of understanding the underlying mathematical equations for each function when obtaining habitat metrics using a software package.The variability in strength of correlations among these metrics warrants future investigation into how each metric, or combination of metrics, is associated with benthic cover, and ultimately, with the structure of reef fish assemblages.
Profile curvature and planform curvature measure curvature in two different directions: profile curvature in the direction of slope and planform curvature in the direction that is transverse to the slope [33].Therefore, it was surprising that these two curvature measures had a correlation of ≈1 in both ArcMap and R. Profile curvature and planform curvature values were obtained for every raster cell of a DEM and were subsequently averaged.Therefore, the correlation of ≈1 between these variables means that curvature parallel to slope and curvature perpendicular to slope were the same when they were averaged across the entire coral reef area at each survey site.As curvature measures the rate of change in slope gradient and direction [33], averaging all the cell values to represent each survey area may be an oversimplification for such non-static variables.Rather, examining the range and/or distribution of the curvature values for each DEM may be more useful for assessing structural dynamics of coral reef habitats.Curvature is also strongly affected by resolutions of DEMs [37].While obtaining curvature values at a 1-cm resolution is particularly useful when comparing structural complexity of coral species [13], using lower resolutions (larger cell sizes) should also be considered depending on the spatial scale of interest.Finally, the negative correlations observed between the profile and planform curvature values obtained in ArcMap and R were due to changes in the plus/minus signs in front of the respective equations in ArcMap [38,39].As there are many different equations for curvature that are implemented in different software packages, with some inconsistent usage of the plus/minus signs as observed in the present study [38,39], one should use caution when interpreting curvature values.
Placing the GCP unit on a completely flat surface prior to image acquisition in the field was difficult due to the complex structure of the benthos.This often resulted in the placement of GCPs where the axes between the GCPs were not parallel along the depth contour of the substrate.Even in these cases, scaling of the 3D models was still correct, but the models became tilted in relation to the contour of the benthic substrate.This tilting of the 3D models required some caution when exporting the DEMs, and in some cases, the models needed to be rotated to ensure the exported DEMs had x and y values corresponding to the horizontal distances along the survey isobath, and z-values corresponding to the planar (overhead "birds-eye-view") angle from which the images were taken.In order to investigate the effects of this model "tilting" on habitat metrics, we performed a post-hoc analysis after randomly choosing 3D models from five sites and haphazardly rotating each of the five models in six different ways in Agisoft PhotoScan.After each rotation, a DEM was exported using the planar projection with "current view" as a projection plane and it was processed in R for fractal dimension, surface complexity, slope, profile curvature and planform curvature.The rotations were done in the extreme in some cases, as observed by the changes in the range of minimum and maximum heights of DEMs (Figure 4g).While it would be unlikely for models with such levels of "tilting" to go unnoticed and be processed further without any corrections, we wanted to test the degree to which tilting of the models affected the resulting metrics of the habitat structure.The rotations did affect the habitat metrics, with sites with higher complexity generally being more affected than sites with lower complexity (Figure 4).Among fractal dimension, surface complexity and slope, fractal dimension had the least amount of changes in the values after rotations (within each site) relative to the ranges of values observed from the 85 sites, although the change was still relatively high for the manually tilted site with the highest habitat complexity (Figure 4a-d).Future studies employing underwater 3D reconstruction techniques on coral reefs should ensure that an adequate number of GCPs are used and that the relative depths of the GCPs are recorded to validate the correct model orientation.As the quality of habitat metrics obtained from DEMs depends on the accuracy of the DEMs, the benefits of having ample GCPs dispersed throughout the study plot and accurately recording the depth of each GCP greatly outweigh the cost of the additional time and effort for divers.The benthic survey using photogrammetric techniques and the subsequent extraction of habitat metrics described here can be used with any ecological monitoring programs that require highresolution data on reef habitat structures.Image collection simply using two swaths along a transect line with two or more GCP units is time-efficient, and the resulting models have sufficient accuracy for extraction of habitat metrics.The use of fractal dimension as a single measure of habitat complexity can simplify and streamline post-survey data processing.This is particularly important if a large number of sites are periodically surveyed and image processing to obtain data needs to keep up with an influx of photographs from survey sites.It also frees up some time if the 3D models are used to obtain other benthic characteristics, such as 3D volume or distributions of benthic species.Further investigations into curvature measures should consider relevant spatial scales and focus on obtaining variables that can complement fractal dimension in characterizing a reef habitat, using DEMs with appropriate resolutions if necessary.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1,Text S1: R code for processing DEMs, including calculation of the fractal dimension, Text S2: R code for a custom function to obtain the slope and curvature.The benthic survey using photogrammetric techniques and the subsequent extraction of habitat metrics described here can be used with any ecological monitoring programs that require high-resolution data on reef habitat structures.Image collection simply using two swaths along a transect line with two or more GCP units is time-efficient, and the resulting models have sufficient accuracy for extraction of habitat metrics.The use of fractal dimension as a single measure of habitat complexity can simplify and streamline post-survey data processing.This is particularly important if a large number of sites are periodically surveyed and image processing to obtain data needs to keep up with an influx of photographs from survey sites.It also frees up some time if the 3D models are used to obtain other benthic characteristics, such as 3D volume or distributions of benthic species.Further investigations into curvature measures should consider relevant spatial scales and focus on obtaining variables that can complement fractal dimension in characterizing a reef habitat, using DEMs with appropriate resolutions if necessary.

Figure 1 .
Figure 1.Map of the Hawaiian Archipelago, including the Northwestern Hawaiian Islands (NWHI).The boundary of the Papahānaumokuākea Marine National Monument is shown by the dotted line.

Figure 1 .
Figure 1.Map of the Hawaiian Archipelago, including the Northwestern Hawaiian Islands (NWHI).The boundary of the Papahānaumokuākea Marine National Monument is shown by the dotted line.

Figure 2 .
Figure 2. Pearson correlation coefficients (ρ) between habitat metrics obtained using the ArcMap path distance toolset, ArcMap surface toolset, ArcMap BTM and R.

Figure 2 .
Figure 2. Pearson correlation coefficients (ρ) between habitat metrics obtained using the ArcMap path distance toolset, ArcMap surface toolset, ArcMap BTM and R.

Figure 4 .
Figure 4. Habitat metrics calculated in R from five randomly chosen sites, showing (a) D64, (b) D128, (c) 3D surface complexity, (d) slope, (e) profile curvature and (f) planform curvature.The 3D models were haphazardly tilted in six different ways, generating six DEMs for each site.Habitat metrics were then calculated from each DEM.The blue circles show values before tilting.The bars show minimum and maximum values from the seven DEMs (one before titling and six tilted), with white circles showing the mean values.Red dotted lines show minimum and maximum values from the 85 sites, without tilting.Differences in the maximum heights and minimum heights, (Max Height)-(Min Height), of the DEMs are also shown (g) in reference to the extent of model tilting used in this posthoc analysis.

Figure 4 .
Figure 4. Habitat metrics calculated in R from five randomly chosen sites, showing (a) D 64 , (b) D 128 , (c) 3D surface complexity, (d) slope, (e) profile curvature and (f) planform curvature.The 3D models were haphazardly tilted in six different ways, generating six DEMs for each site.Habitat metrics were then calculated from each DEM.The blue circles show values before tilting.The bars show minimum and maximum values from the seven DEMs (one before titling and six tilted), with white circles showing the mean values.Red dotted lines show minimum and maximum values from the 85 sites, without tilting.Differences in the maximum heights and minimum heights, (Max Height)-(Min Height), of the DEMs are also shown (g) in reference to the extent of model tilting used in this post-hoc analysis.

Table 2 .
List of habitat metrics obtained and the variable names used in the present study.Software and a toolset/license for ArcMap used for each metric are also listed.