3D Ocean Water Wave Surface Analysis on Airborne LiDAR Bathymetric Point Clouds

: Water wave monitoring is a vital issue for coastal research and plays a key role in geomorphological changes, erosion and sediment transportation, coastal hazards, risk assessment, and decision making. However, despite missing data and the difﬁculty of capturing the data of nearshore ﬁeldwork, the analysis of water wave surface parameters is still able to be discussed. In this paper, we propose a novel approach for accurate detection and analysis of water wave surface from Airborne LiDAR Bathymetry (ALB) large-scale point clouds data. In our proposed method we combined the modiﬁed Density-Based Spatial Clustering of Applications with Noise (DBSCAN) clustering method with a connectivity constraint and a multi-level analysis of ocean water surface. We adapted for most types of wave shape anatomies in shallow waters, nearshore, and onshore of the coastal zone. We used a wavelet analysis ﬁlter to detect the water wave surface. Then, through the Fourier Transformation Approach, we estimated the parameters of wave height, wavelength, and wave orientation. The comparison between the LiDAR measure estimation technique and available buoy data was then presented. We quantiﬁed the performance of the algorithm by measuring the precision and recall for the waves identiﬁcation without evaluating the degree of over-segmentation. The proposed method achieves 87% accuracy of wave identiﬁcation in the shallow water of coastal zones.


Introduction
Water wave surface monitoring at nearshore plays a key role in facilitating the updating of basic geographic information including geomorphological change, coastal hazards, erosion, and sediment transportation. It aids the management and maintenance of coastal areas. In addition to the function of damaging and warning system, there is a high probability of nearshore infrastructures being damaged in locations with waves. Hence, it is crucial to detect the waves before they reach coastal zones. The main challenge of ocean water wave detection is primarily due to the difficulty of acquiring high-quality data especially in the nearshore and breaking surf zone. Figure 1 shows several typical waves on nearshore shallow water.
Wave height estimation analysis has been ongoing research throughout the last century. Many researchers have studied the measurement of wave heights in coastal zones via Synthetic Aperture Radar (SAR) satellite imaging. Traditionally, a theoretical and parametric analysis was performed for the retrieval of the two-dimensional ocean waves spectra from SAR [1][2][3]. In the literature review, the C-band of Sentinel-1 SAR imagery has a wide range of products for ocean wave analysis (e.g., [4][5][6]). Parameters such as the peak of wavelengths and wave direction are calculated from C-band Scan SAR images spectra with a semi-empirical algorithm. Smit et al. [7] assimilate significant wave height from a deployed network of free-drifting satellite-connected surface weather buoys. Their assimilation strategy improves forecast accuracy with a 27% reduction in root-mean-square error in significant wave heights overall.  Furthermore, wave spectra derived from SAR images acquired by ENVISATs and SEASAT SAR were studied (e.g., [8][9][10][11]). The wave spectra inversion scheme is tailored for shallow water. The results were compared to in situ measurements by seven sensors, deployed in a field experiment. They found a reasonable agreement between in situ and SAR observations for wave heights and directions. The TanDEM-X SAR satellite also was used for ocean surface observations. There is a sensor capability within the contexts of tidal current mapping and wind-driven ocean currents as well as detection of areas with surface films [12,13].
Some other researchers employ Artificial Neural Networks (ANN) on SAR images. For instance, Ardhuin et al. [14] and Wang et al. [15] used ANN to understand the nonlinear relationship between the input SAR images parameters and output geophysical wave parameters with a data combination of altimeters and in situ buoys. Likewise, deep learning was applied to forecast time series of Significant Wave Height (SWH) [16]. They proposed a machine learning technique to create a model that relates the ocean image properties to geophysical wave parameters and predicts the sea state parameters. Their datasets comprise the SWH and wind speed data which were captured from six buoy stations in the Taiwan Strait and its adjacent waters.
Besides normal meteorological conditions, measurement of wave heights in ice zones can be implemented through SAR imaging Sentinel 1A [17,18]. They developed a nonlinear algorithm to assess the phase-resolved deterministic maps of wave-induced orbital velocities. The performance of their method for the estimation of wave parameters is expected to work best when the shortest wave components can be neglected. Despite some good results, ocean waves are not always visible in SAR images and there is a gap in detection criteria in terms of wave height, length, and direction [10]. Therefore, researchers investigated other approaches.
Microwave-radar based products were used for ocean surface observation (e.g., [19][20][21]). This research monitors the surface roughness, which is the source of the radar signal.
Hwang et al. [21] studied on surface roughness and foam which are the two main contributors to ocean surface microwave thermal emission. Other research used GNSS reflectometrybased on the signal-to-noise ratio (SNR) for ocean surface observation. The distance between an antenna and the water surface is considered by analyzing the oscillation of the SNR observation [22][23][24][25]. The relationship between the sea surface roughness and the attenuation of the SNR oscillation will help to estimate the numerous significant wave height (SWH) by analyzing the SNR signal. The lack of dimension of the surface roughness spectrum and wind speed is empirically parameterized and defined as the ratio of wind friction velocity and phase speed of the surface roughness wave component [26]. Regarding numerical mathematics models, a developed model was proposed for the elevation spectrum of wind-generated ocean waves. The proposed model is a modification of the spectrum. The result of these modifications is a model simplification compared with the Pierson model [27].
Some other researchers have investigated the wave breaking analysis of the nearshore zones in both field study and laboratory experiments. In the outdoor field study, research was conducted at a tropical microtidal intermediate sandy beach to investigate the shortterm swash-zone hydrodynamics and morphodynamics under variable wave conditions. The data was obtained from a 2D Light Detection and Ranging (LiDAR) scanner for wave height calculation at the lower foreshore and swash-induced topographic changed [28]. Another outdoor field experiment used 2D LiDAR scanners deployed from the pier at Saltburn-by-the-Sea, the UK for 6 days. The scanners monitored the surface elevation of nearshore waves from the breaking point to the run-up limit of temporal and spatial to capture the time-varying free surface throughout the surf and swash zones [29]. In the indoor laboratory experiment [30], the breakpoint location and time-varying water surface wave breaking were measured by mixing LiDAR and Sonar data. They used three SICK LMS511 commercial 2D LiDAR scanners along with multibeam sonar to obtain a complete surface elevation. Bryan et al. [30] obtained a good performance on breaking point location and changing wave height with an assessment of LiDAR data.
Additionally, some other methods were done on 2D imagery data such as infrared imagery and stereo frames. In this case, Deep Convolutional Neural Networks (DCNNs) perform to estimate wave breaking types from infrared imagery at the surf zone. Then, they classify the breaker type through logistic regression on the extracted image features. Their classification accuracy reached 89% and 93% [31,32]. In other research, the development of a remote wave gauging technique performs to estimate wave height and period from the imagery of waves in the surf zone. In this study, their implementation was done on the three image-based datasets. The suggested network is trained using integration of coincident images and in situ wave measurements. They achieved RMS errors of 0.14 m and 0.41 s for height and period [33]. In particular, [34][35][36] proposed a series of procedures for coastal wave-tracking using coastal video imagery with deep neural networks. The existing methods consist of stereo frames based methods. The speed of instant wave speed in the video domain is estimated through learning the behavior of propagated waves in the surf zone.
Recently, Airborne LiDAR widely used for coastal shallow water monitoring [37][38][39]. In this regard, the impact of upstream offshore wind farms on sea wave state has been studied with systematic flights deploying an airborne laser scanner. The analysis of the spectral energy distribution represented a re-distribution of the wave energy within the downstream area with enhanced energy at smaller wavelengths. Their results clearly show the effect of upstream wind at a distance of 55 km. Moreover, the wave profiles and coastal sea topography were studied with a commercial portable 2D laser scanner deployed from a fixed-wing aircraft. It presents sea surface topography and wave profiles on low altitude surveys (<300 m). Relative observed wave heights compared with a waverider buoy located in the same area during the LiDAR survey. The results show variations of the sea surface to within approximately 10 cm. In recent research, airborne LiDAR was used to measure ocean waves. The LiDAR data captured in the German Bight and processed for significant wave height which could be used to improve Wave Analysis Model (WAM). The results showed that SWH varied between 0.4 m and 2.3 m, and data is validated against buoy measurements during direct overpasses. They have below 15% deviation with the outputs of the third-generation numerical WAM [35,40].
Except for surface detection, the internal waves in shallow water are detected under the right conditions along with a plankton layer associated with a density gradient [41][42][43]. Likewise, two-dimensional measurements of ocean wave displacement were studied via an airborne LiDAR scanning system. The active ranging by the scanning optics can obtain passive measurements of surface emissivity and process them into a binary image. These measurements can provide information on the average statistics and the spatial distribution of breaking waves [44].
Few efforts used 3D LiDAR point clouds for ocean water wave detection and analysis. Our proposed method was done on large-scale 3D LiDAR point cloud dataset. Our method can adapt to different types of oceanic, river, and wetlands water wave shape. Figure 2 represents the flowchart of our method. The main contributions of this paper are summarized as follows: (1) Presenting a novel water wave detection pipeline-interpretable method-which can simultaneously target different types of water waves in Airborne LiDAR Bathymetry (ALB) point clouds.
(2) Estimation of wave parameters such as wave height, wave length, and wave orientation.
(3) Exploring a modified Density-Based Spatial Clustering of Applications with Noise (DBSCAN) method with a connectivity constraint and a multi-level 2.5D analysis of water surface. The paper continues below. The proposed method, experimental results, and discussion are described in Sections 2-4 respectively. Furthermore, the last section presents the conclusion and future directions.

Data Description
The ALB LiDAR data in this paper are captured by a RIEGL VQ-880-G scanning system, as shown in Figure 3. The RIEGL VQ-880-G scanning system integrates laser scanners, RGB camera, and inertial navigation devices which are installed on traveling aircraft to capture accurate, high density, and precise 3D point clouds. As opposed to airborne topographic LiDAR, which uses an infrared wavelength of 1064 nm, bathymetric LiDAR systems use a green wavelength of 532 nm to penetrate the water column for measuring the seafloor. The infrared laser beam pattern is linear and the green laser beam acts in a circular pattern. The dataset recorded by 3D LiDAR sensor is mounted on a moving air-vehicle passing through various coastal-urban scenes. The data was acquired in the beachside areas of Tampa Bay, Florida, United States, on 15 October 2015, during 1 flight in the afternoon (Time zone-GMT-05:00). Figure 4 shows the flight trajectory of Airborne LiDAR. The study area covers a large geospatial area (latitude 27 • 51 35 N, longitude 82 • 50 54 W), as a part of Florida beachside. The total number of ALB point clouds at shallow water are around 5 million points by 1.2 km × 1.8 km around.
The raw captured point cloud is adjusted with the National Oceanic and Atmospheric Administration (NOAA) benchmarks as ground control points. This transformation results in having geo-referencing for further calculation. Figure 5 shows the wave in deep and shallow water through an Airborne LiDAR Bathymetry (ALB) system at coastal urban zone. Likewise, it represents the original and geo-referenced bathymetry point clouds.
The proposed method combines a modified DBSCAN clustering method with a connectivity constraint and a multi-level 2.5 D analysis to detect water wave surface. Our approach is adaptive for most types of wave shape anatomies in deep and shallow waters, nearshore and onshore. Figures 6 and 7 show the global and detailed pipeline of water wave detection. The input ALB point clouds are shifted onto the refined coordinates system, by selection planar surface inlier and discretization to voxel indices. The resulting voxel indices set will be used for the 2.5D projection. In detail, the description of each step of our pipeline is shown in Algorithm 1. Algorithm 1 waves detection from point clouds() Input: A point cloud P Input: n x , n y , n z the number of range subdivisions by axis Input: α, the standard deviation factor, used for plane inliers selection Input: θ, the rotation angle step used for the minimum bounding box computation (typically 1 degree) Input: The approximate size of the 2.5D height image: m × n Input: The threshold δ s.t. any fluctuation with a magnitude lower than δ is considered stationary Result: A copy of P with and additional field label. The label is 0 for stationary points, ≥0 for waves clusters, and −1 otherwise (e.g., underwater points) 1: (u 1 , u 2 , u 3 ), I plane = LRF(P, n x , n y , n z , α, θ), (LRF is LocalReferenceFrame) 2: Compute P local the projection of the point cloud onto the local reference frame 3: P bins , p min , r = Discretize(P local , n x , n y , n z ) 4: Compute the height map: I, M, m 0 , n 0 = ComputeHeghtMap(P bin , m, n) 5: Extract I c the sub-image of I without empty pixels 6: Decompose I c in low frequency and high frequency coefficients using 2D wavelets analysis 7: Compute I low the reconstruction of I using the low frequency components exclusively 8: I d = I − I low 9: I stationary , I + , I − = SeparateFluctuations(I d , δ) 10: I +,labels = DBSCAN(I + , √ 2, connectivetyDistance), we use √ 2 as a growing radius to exclusively allow merging pixels within 1 pixel of distance the current cluster 11: I −,lables = DBSCAN(I − , √ 2, connetivityDistance) 12: P labels = AssignLabelsToPointCloud(P, I stationary,I +,labels ,I −,labels ) return: P labels

Planar Surface Extraction
We extract the water surface jointly with a coordinate system (local reference frame) in which the maximum and minimum 3D points draw a minimum volume bounding box around the water surface inliers. In this regard, first, we extract 3 principal components using Principal Component Analysis (PCA). The component with the lowest variance approximates the water surface's normal vector since the point cloud is extended horizontally. Then, to separate the water from under-water surface, we conduct a discretization of the point cloud that acts like a standardization method. In summary, we divide the range of coordinates values in each axis to N bins intervals. Each interval is then indexed from 0 to N bins . Afterwards, the 3D space is partitioned into a set of 3D boxes, or voxels, where each box is referred by 3 bins indices. Every 3D point is then replaced by the coordinates of the box where it belongs. We then use a linear regression coupled with an outlier detector as Random Sample Consensus (RANSAC) [45] to find the planar surface with maximum inliers along with the inlier's indices (a list of indices from 1 to N points ). Subsequently, the inliers are then used exclusively to re-estimate the local reference frame using PCA followed by an estimation of the minimum bounding box of the planar surface. The axes of the bounding box form the refined local reference frame. We finally project the input point cloud onto the refined coordinates system, select planar surface inliers, and discretize them to voxel indices. In Figure 8, the results of voxel indices set will be used for the 2.5D projection.

2.5D Projection
This step aims to project the planar point cloud onto an image of signed distances to the planar surface. We first create a fictional planar grid parallel to the planar surface (or share the same normal vector of the detected plane). The 3D plane inliers are then projected onto the grid. In other terms, we assign to each 3D point the row and column indices of the grid cell where this point falls into it. Finally, we compute the average height of points in each cell and store the value in the pixel value corresponding to cell indices. Figure 9 shows the 2.5D projection. Note that we keep track of points belonging to each cell by storing their indices and cell indices in a map. This allows the recovery of 3D point clouds from pixel clusters. The grid size must be carefully chosen to minimize the number of empty cells while keeping the intrinsic details of the shape of the wave. Since the tile is not guaranteed to be square-shaped, we expect to have boundary empty cells in the projection. For this reason, we extract the maximum size sub-image without empty pixels from the original projection.

Wavelets Analysis
We use a combination of two wavelet analysis techniques, which are described as follows.

Wavelet Denoising
Decompose the image into a global shape (low frequency) and details (high frequency). Then cancel reconstruct again. The resulting image should look like a trend formed by the wavelet kernel [46]. Here we used a Daubechy kernel that results in locally linear trends (the neighboring pixels' values are approximated by a linear function of pixels' coordinates).

Signal Detrending Using Wavelets
The difference between the original image and the image of local trends will give an image of residuals. The residuals are chunks of the original 2D signal that separates it from a linear one. Figure 10 shows the wavelet analysis procedure.

Clustering
We perform clustering of each fluctuation image using a custom distance that: (1) Penalizes empty pixels (connectivity constraint) by returning and infinite as the distance between an empty pixel and a valid pixel. (2) Returns the Euclidean distance between the pixels' indices otherwise. Note that we do not use pixel values for the clustering as the shape of convex and concave waves are already separated by finding a suitable threshold for static points, as shown in Figure 11.

Recover 3D Points with Clustering Labels
After the clustering with DBSCAN, each wave pixel is assigned to a label and to an entry in the height map storing the indices of 3D points belonging to the pixel. To recover the clusters' point clouds, we add a label column for the 3D points, and we proceed as follows for each image (2 fluctuation images and a static pixels' image): -Initialize a label counter = 0; -For each label in the cluster, use the height map to recover all the indices of 3D points belonging to the cluster; · Replace the column label by the label counter value; · Increase the label counter and continue; -Return the point cloud with clusters' labels.

Wave Parameters Measurement
We estimate the wave parameters such as wave height measurement through Fourier Transform (FT) approach. Below Equations (1)-(3), show the sinusoidal wave Fourier equations: e j2π(ux+vy) = cos 2π (ux + vy) + J sin 2π (ux + vy) where u and v are spatial frequencies. In addition, one can write FT pairs as f (x, y) = F(u, v). The magnitude of the vector (u, v) gives a frequency, and its direction gives an orientation. The function is a sinusoid with this frequency along the direction and constant perpendicular to the direction [47]. The frequency vector is starting from the center to the white dot (max magnitude Fourier). The orientation is aligned with the wave orientation which means that the angle of the vector equals the orientation of the wave and its norm is related to wavelength. We calculate the angle by code and rotated the wave image to make it horizontal. Then we got the cross-section that travels the whole wave. Figures 12 and 13 show the Fast Fourier Transfer (FFT) approach [48]. We analyzed the spectrum to extract wave parameters measurement locally and stored them in additional values. By repeating this to every patch of 21 pixels × 21 pixels around a central pixel, the result is a feature map of 3 channels. 1st layer: wave angle, 2nd layer: wavelength, and 3rd layer: wave height.

Results
In this section, first the water wave detection results are present. Then through the FFT approach, we show the results of wave parameters analysis.

Water Wave Detection Result
The adjustment of the parameters for the proposed algorithm is intuitive. The experiment results prove the effectiveness of the algorithm, as shown in Figures 14 and 15. Figure 14a shows the top view of original point cloud at nearshore and shore and Figure 14b shows related water wave cluster detection result on the near-shore and on-shore zones. Figure 15 shows the results as 3D point cloud. We illustrate the situations in which the nearshore waves and breaking waves are placed differently: the near-shore waves are placed in the approaching shore-waves and breaking waves are placed in the breaking surf zone.
The wave cluster detection is based on wavelet analysis. The projection of the point cloud using the applied grid scheme is seen as a 2D signal and decomposed into low frequency and high frequency wavelet components. In our experiment, we used the Daubechie kernel with 3 levels of decomposition since it matches the wave patterns we want to capture. The high frequency components are then removed (overridden with zeros) and the signal is reconstructed. The reconstruction shows local trends in the point cloud projection. We then compute the residual image as the difference between the heights in the projection of the point cloud and the heights in the reconstructed projection. This difference indicates how far the signal variation is from the local trend. Typically, the residual image shows locally concave and convex crests. We then apply a threshold to remove small fluctuations and obtain a binary mask where high enough fluctuation pixels are associated with white pixels and small fluctuations are associated with black ones. The mask is finally applied to the original projection and produces Figure 14b. In the mentioned figure, the rejected fluctuations are shown with a blue mask. Selected points are shown in grayscale.

Water Waves Analysis
In our study, we analyze the wave height, wavelength, and wave direction through the ability of the Fast Fourier Transform (FFT) approach. Since FFT outputs are a spectrum, we analyze the spectrum to extract parameter measurements locally and stored them in additional values (Figure 12). We used Fast Fourier Transformation and expressed it as a linear sum of sine waves for each local point cloud (Z is a function of x and y). We picked dominant frequency assuming that locally, we should have 1 sine wave. Although it is 2D Fourier, however, pixel values store height. So we are decomposing height as a sine wave. Figure 16 shows the results of wave height, wavelength, and wave angle. In the height-map, each 3D point is associated with the wave height as a scalar field. Notice that red regions are associated with high fluctuation, whereas green and cold regions represent small ones. The blue points at the boundary are missing data due to convolution. The wave height does not show clear clusters of waves since several connected wave parts might share the same height. However, it shows regions of almost the same color. These indicate the extent of some frequencies: i.e., the extent of wave amplitude in space, or simply, the regions where the same height information can be used to describe the variation (Figure 16a). Similarly, in the wavelength map each 3D point is associated with the wavelength as a scalar field. All values seem to be relatively similar (cold range of colors), which makes sense visually from the length of the crests in the 3D render.
The wavelength information coupled with cell positions shows well defined wave parts (Figure 16b). In angle-map, each 3D point is associated with the wave angle as a scalar field as well. The selection of the points corresponding to the dominant angle shows large connected components (Figure 16c). The sea state parameters derived from ALB data provide information over a large area, including the coastal zone with similar products to the other altimetry products. In comparison, the low-resolution altimetry wave analysis and open ocean SAR wave mode products are not sufficient for local and regional applications in the complex coastal environment. Therefore, the ALB wave data would provide added value dealing with coastal processes. However, this is just one period scan time and we did an estimation of wave parameters by available data. Figure 16, Table 1, Table2 and Table3 show the results of wave height, wavelength, and wave angle at the nearshore zone.

Significant Wave Height Estimation
The natural definition of wave height (H) can be as the vertical distance between the highest and the lowest surface elevation in a wave shape [49]. The wave condition in a stationary record can be characterized with average wave parameters, such as the significant wave height and the significant wave period [50]. The significant wave height (SWH) is in proper correlation with the wave height as estimated 'visually' by experienced observers. However, this is not true for the significant wave period. Significant wave height (SWH) is defined as the mean of the highest one-third of waves in the wave record: where j is not the sequence number in the record (i.e., the sequence in time) but the rank number of the wave, based on wave height (i.e., j = 1 is the highest wave, j = 2 is the second-highest wave, etc.). As shown in Table 1, the agreement between the estimated Significant wave height (SWH) and the measured buoy data is reasonable. To this end, regarding Tables 2 and 3 and Equation (4), the SWH after geo-reference transformation is Hs = 0.49 m, which confirms with available buoys measurement data (station 3).

Discussion
In this section, in order to evaluate the performance of our method, ALB data was evaluated with available buoy data. Additionally, we calculated the assessment metrics with an annotated version of projection of the point cloud that we created.

In Situ Measurements
The sea surface is measured by directional floating buoys and the processed data is available every 30 min. In this study, three buoy stations observe the water level around the study area (station 1 'Clearwater Beach, FL', station 2 'Fred Howard Park, FL', and station 3 'Egmont Channel Entrance, FL') along with wind speed, air temperate, significant wave height, and wind direction, Figure 4. However, station 1 and 2 did not record any wave height data during flight missions in 2015. Therefore, there is missing wave height data. The only available wave height data belongs to station 3, which was far (around 40 km) and the average wave height in station 3 is between 0.4 m and 0.7 m in 15 October 2015. This correctly confirms the algorithm output. Moreover, this new observation with ALB compensates for the missed data and shows the high capabilities of Airborne LiDAR when in situ measurement data is missed.

Satellite Observations
The launched operational altimetry satellites, e.g., Sentinel-3A and sentinel-6 were references for sea-surface height measurements until at least 2030. Satellite altimetry is radar technique-based for ocean surface height measurements with a precision of a few centimeters. However, due to the lack of available altimetry data during the study project (15 October 2015), we assess our algorithm performance with the available buoy dataset (station 3) and annotation point cloud (Section 4.3).

Quantitative Assessment Measures
We calculate the assessment metrics based on each wave of nearshore point clouds, where the precision and recall are respectively defined as Equations (5) and (6) the F1measure is adopted to comprehensively evaluate the performance as Equation (7).
N TP is the number for the correct waves points extracted; N FP is the number for the wrong points extracted; N FN is the number for the waves points that failed to extract. The quantitative results are shown in Table 4. To evaluate the performance of our method, we created an annotated version of the projection of the point cloud as an orthoimage ( Figure 17). A human labeler defines manually the most prominent waves, in other terms, the waves with the largest support. Note that in the original point cloud, waves are invisible without any processing. The water surface appears to be flat. This is one of the reasons why we chose to label the 2.5D projection. Moreover, small fluctuations are nearly invisible to the human eye. Each wave segment is associated with an index. The predictions are also associated with indices that do not necessarily match true segments. The number of segments and the ordering of predicted labels are completely independent from the annotation. As a future work, we can investigate the creation of a complete labeling framework.
Visually, we notice that the predicted wave presents small clusters that collectively overlap with most of the true clusters. We conclude that our approach provides an oversegmentation of the water surface. As a perspective, we can continue the processing by grouping segments into consistent wave parts. This requires the definition of wave consistency and the implementation of an appropriate merging policy that achieves a trade-off between over-segmentation and under-segmentation.
We quantify the performance of the algorithm by measuring the precision and recall for the wave's identification without evaluating the degree of over-segmentation. We label all wave parts with +1 and the remaining pixels with −1. The predictions labels are changed accordingly. Note that rejected pixels (−1 as prediction) result from the small fluctuations filtering followed by DBSCAN outlier detections. The small fluctuations on the bottom right of the image are hard to identify by a human labeler. These are, however, detected by our system. The percentage of these pixels is, hence, expected to be small relative to true −1 pixels. In our evaluation, the recall of class −1 is 52% which suggests that almost 48% of 'not wave' pixels were predicted as wave, such that the percentage of precision class −1, relative to all rejected pixels is 22%. Visually, an important part of these pixels is located at the "invisible" fluctuations. For the wave class, the precision is relatively important and associated with 64% recall. Similarly, the predicted wave class matches the shape of true waves except for small fluctuations. However, the support thickness of the wave segment was not predicted accurately. The missed pixels on the borders contribute to the decrease of the recall.

Complexity Analysis
The algorithm takes 2 min: 4 s ± 10 s on Ubuntu 20.04 operating system. The time was measured while having all applications and internet connections shut down. The desktop GUI and visualizations saving were deactivated. The algorithm was executed on an i7 9th generation laptop on a single processor and without use of the Graphics Processing Unit (GPU).

Conclusions
In this paper, we proposed an algorithm for ocean wave detection and wave analysis parameters measurement through the Airborne LiDAR Bathymetry (ALB) point cloud. Our experiments show that different types of wave shape anatomies of the ocean surface can be detected and quantified by using ALB system. By leveraging an addressing strategy, the proposed method rapidly handled large volumes of ALB point cloud toward 3D water wave detection. Our algorithm was successfully applied to the 3D ALB dataset and achieves precision of 87% for wave identification and associated with 64% recall at the shallow water of coastal zones. Visually, the predicted wave class matches with the shape of true wave except for small fluctuations. Moreover, the only available buoy station data confirms the algorithm output. The ALB wave analysis through the FFT approach provides more detailed information about spatial parameters of the wave field in the coastal zone compared with a model forecast, especially when in situ measurement data are missed. Indeed, the integration of time-varying surface data and field study data with 3D LiDAR point clouds data will improve the wave analysis accuracy measurement for real-time processing. Further research is required for additional information, such as video imagery and adding a more accurate field study detailed dataset.