Application of Unmanned Aerial Vehicles and Image Processing Techniques in Monitoring Underwater Coastal Protection Measures

: A prerequisite for solving issues associated with surf zone variability, which affect human activity in coastal zones, is an accurate estimation of the effects of coastal protection methods. There-fore, performing frequent monitoring activities, especially when applying new nature-friendly coastal defense methods, is a major challenge. In this manuscript, we propose a pipeline for performing low-cost monitoring using RGB images, accessed by an unmanned aerial vehicle (UAV) and a four-level analysis architecture of an underwater object detection methodology. First, several color-based pre-processing activities were applied. Second, contrast-limited adaptive histogram equalization and the Hough transform methodology were used to automatically detect the underwater, circle-shaped elements of a hybrid coastal defense construction. An alternative pipeline was used to detect holes in the circle-shaped elements with an adaptive thresholding method; this pipeline was subsequently applied to the normalized images. Finally, the concatenation of the results from both the methods and the validation processes were performed. The results indicate that our automated monitoring tool works for RGB images captured by a low-cost consumer UAV. The experimental results showed that our pipeline achieved an average error of four pixels in the test set.


Introduction
Nearshore morphology is continuously changing. Occasionally, this is due to drastic transformations during extreme storm events, and gradual change occurs because of dayto-day interactions with waves and currents [1][2][3][4]. Numerous coastal change models have been based on simplified global shoreline movement patterns [4][5][6] and a site-specific local and regional analysis that implements real coastal measurements-which are mostly surveyed as a cross-section of the coast [7][8][9][10][11][12]-to provide evidence for the vulnerability of sandy coasts to erosion. Therefore, the appropriate management of coastal zones is closely linked to the balance between natural processes and protection techniques [13]. For years, coastal protection measures have been dominated by hard structural solutions, such as seawalls, breakwaters, groins, and dikes [14]. However, traditional coastal defenses can no longer keep pace with climate change and continuous coastal population growth [15]. The support of new nature-friendly methods of coastal protection planning is becoming a major

Study Area and Materials
Artificial reefs, constructed from concrete habitat modules, are new to the Baltic Sea, and no monitoring of such underwater structures has been conducted before. They were constructed under the project "Protection of the sea shore at Łeba, Rowy, and Ustka" in 2014-2016. During the construction phase, aerial images using UAV were taken of all underwater thresholds in Ustka, Rowy, and Łeba ( Figure 1). The first orthophotos of all structures were created in 2016. Due to the apparent significant changes in the structure of the underwater thresholds between years and the high concentration of concrete elements in the structure, the thresholds in Rowy ( Figure 2) were selected for further study.
Remote Sens. 2022, 14, x FOR PEER REVIEW 3 of 18 regional significance. Furthermore, only a few studies have highlighted the usefulness of UAVs within integrated coastal zone management, especially for the documentation of post-storm coastal erosion, as well as monitoring damage to coastal infrastructures [33,40,41]. Thus, reliable information about the behavior of underwater structures benefits and enables numerous strategic decisions in coastal protection processes.

Study Area and Materials
Artificial reefs, constructed from concrete habitat modules, are new to the Baltic Sea, and no monitoring of such underwater structures has been conducted before. They were constructed under the project "Protection of the sea shore at Łeba, Rowy, and Ustka" in 2014-2016. During the construction phase, aerial images using UAV were taken of all underwater thresholds in Ustka, Rowy, and Łeba ( Figure 1). The first orthophotos of all structures were created in 2016. Due to the apparent significant changes in the structure of the underwater thresholds between years and the high concentration of concrete elements in the structure, the thresholds in Rowy ( Figure 2) were selected for further study.  This section of the Polish coast exhibits diverse characteristics, but primarily assumes the form of eroded cliffs. The most dynamic coastal zone is the beach, which experiences cyclical changes (i.e., increasing from spring to summer and decreasing due to marine erosion in the autumn and winter months) [42]. Remote Sens. 2022, 14, x FOR PEER REVIEW 3 of 18 regional significance. Furthermore, only a few studies have highlighted the usefulness of UAVs within integrated coastal zone management, especially for the documentation of post-storm coastal erosion, as well as monitoring damage to coastal infrastructures [33,40,41]. Thus, reliable information about the behavior of underwater structures benefits and enables numerous strategic decisions in coastal protection processes.

Study Area and Materials
Artificial reefs, constructed from concrete habitat modules, are new to the Baltic Sea, and no monitoring of such underwater structures has been conducted before. They were constructed under the project "Protection of the sea shore at Łeba, Rowy, and Ustka" in 2014-2016. During the construction phase, aerial images using UAV were taken of all underwater thresholds in Ustka, Rowy, and Łeba ( Figure 1). The first orthophotos of all structures were created in 2016. Due to the apparent significant changes in the structure of the underwater thresholds between years and the high concentration of concrete elements in the structure, the thresholds in Rowy ( Figure 2) were selected for further study.  This section of the Polish coast exhibits diverse characteristics, but primarily assumes the form of eroded cliffs. The most dynamic coastal zone is the beach, which experiences cyclical changes (i.e., increasing from spring to summer and decreasing due to marine erosion in the autumn and winter months) [42]. This section of the Polish coast exhibits diverse characteristics, but primarily assumes the form of eroded cliffs. The most dynamic coastal zone is the beach, which experiences cyclical changes (i.e., increasing from spring to summer and decreasing due to marine erosion in the autumn and winter months) [42]. The monitored underwater construction in Rowy is a 750 m long structure that is located approximately 170 m offshore at depths ranging from 2 to 4 m ( Figure 3). It consists of 680 concrete coils, each 2.4 m in diameter, with a 0.8 m diameter void at the top. In 2016, the 1st module on the west side contained 121 elements, the 2nd and 3rd modules contained 143 elements each, the 4th module contained 150 elements, and the last module contained 123 elements. All groups were separated by four segments made of stone modules. The monitored underwater construction in Rowy is a 750 m long structure that is located approximately 170 m offshore at depths ranging from 2 to 4 m ( Figure 3). It consists of 680 concrete coils, each 2.4 m in diameter, with a 0.8 m diameter void at the top. In 2016, the 1st module on the west side contained 121 elements, the 2nd and 3rd modules contained 143 elements each, the 4th module contained 150 elements, and the last module contained 123 elements. All groups were separated by four segments made of stone modules. The data used in this study cover three separate UAV surveys conducted on 19 May 2016, 8 February 2018, and 10 September 2021. To collect aerial data from the first 2 flights, we used a low-budget DJI Inspire 1 unmanned aerial vehicle with a 12 MPix camera, 3.6 mm lens, F/2.8 aperture, ISO 200, and a shutter speed of 1/120 of a second. The last flight was made with a DJI Matrice 300 RTK drone, with a 12 MPix camera, 4.5 mm lens, wideangle lens, F/2.8 aperture, ISO 110, shutter from 1/80 to 1/200, depending on the light. The time of data acquisition was determined by the angle of incidence of the solar rays to minimize the reflections on the water surface in relation to the position of the camera lens. Finally, the measurements were started at 8:30 a.m. in 2016, at 11:30 a.m. in 2018, and at 8:00 a.m. in 2021. The first and third survey were conducted without additional filters on the lens, and the second survey used a polarizing filter to eliminate light reflections on the water surface.
Flight parameters, such as flight altitude, spacing between the axes of the series, and distance between the centers of subsequent photos (i.e., the length of the shooting base, which determines where the shutter is turned on and off) were automatically calculated after entering the expected output pixel size. The values included the camera parameters, longitudinal and lateral coverage, and vehicle speed. The mission plan was formulated using the DJI software. Sidelap and overlap during the mission were set at 75%.
The experimental assumptions for using a low-cost UAV required the carefully planning of all surveys based on the weather; therefore, wind and wave conditions in the two surveys were similar, to allow a comparison of the acquired data. Constant observation of meteorological conditions was conducted at the Interdisciplinary Centre for Mathematical and Computational Modelling of the Warsaw University website; and the state of the sea was monitored using the high-resolution operational wave model (WAM), which was validated for the Baltic Sea in the framework of the Hindcast of Dynamic Processes of the Ocean and Coastal Areas of Europe (HIPOCAS) project [43]. Moreover, before the final decision to initiate the survey, images from the web cameras in the port of Rowy were analyzed.
The weather conditions during the 2016 survey were extremely suitable; that is, southern winds with a maximum speed of 2 m/s were observed, and no precipitation, cloud cover, or waves were identified. During the second survey in 2018, the weather conditions were extremely similar; that is, no waves, clouds, or precipitation, were spotted, and the maximum speed of the southerly wind was 2 m/s. The temperature difference was significant (12 °C and −12 °C in the first and second surveys, respectively). During the last survey in 2021, the weather condition was similar, with wind at a maximum speed of 3 m/s, and no precipitation, cloud cover, or waves were identified. The temperature was about 15 °C. The data used in this study cover three separate UAV surveys conducted on 19 May 2016, 8 February 2018, and 10 September 2021. To collect aerial data from the first 2 flights, we used a low-budget DJI Inspire 1 unmanned aerial vehicle with a 12 MPix camera, 3.6 mm lens, F/2.8 aperture, ISO 200, and a shutter speed of 1/120 of a second. The last flight was made with a DJI Matrice 300 RTK drone, with a 12 MPix camera, 4.5 mm lens, wide-angle lens, F/2.8 aperture, ISO 110, shutter from 1/80 to 1/200, depending on the light. The time of data acquisition was determined by the angle of incidence of the solar rays to minimize the reflections on the water surface in relation to the position of the camera lens. Finally, the measurements were started at 8:30 a.m. in 2016, at 11:30 a.m. in 2018, and at 8:00 a.m. in 2021. The first and third survey were conducted without additional filters on the lens, and the second survey used a polarizing filter to eliminate light reflections on the water surface.
Flight parameters, such as flight altitude, spacing between the axes of the series, and distance between the centers of subsequent photos (i.e., the length of the shooting base, which determines where the shutter is turned on and off) were automatically calculated after entering the expected output pixel size. The values included the camera parameters, longitudinal and lateral coverage, and vehicle speed. The mission plan was formulated using the DJI software. Sidelap and overlap during the mission were set at 75%.
The experimental assumptions for using a low-cost UAV required the carefully planning of all surveys based on the weather; therefore, wind and wave conditions in the two surveys were similar, to allow a comparison of the acquired data. Constant observation of meteorological conditions was conducted at the Interdisciplinary Centre for Mathematical and Computational Modelling of the Warsaw University website; and the state of the sea was monitored using the high-resolution operational wave model (WAM), which was validated for the Baltic Sea in the framework of the Hindcast of Dynamic Processes of the Ocean and Coastal Areas of Europe (HIPOCAS) project [43]. Moreover, before the final decision to initiate the survey, images from the web cameras in the port of Rowy were analyzed.
The weather conditions during the 2016 survey were extremely suitable; that is, southern winds with a maximum speed of 2 m/s were observed, and no precipitation, cloud cover, or waves were identified. During the second survey in 2018, the weather conditions were extremely similar; that is, no waves, clouds, or precipitation, were spotted, and the maximum speed of the southerly wind was 2 m/s. The temperature difference was significant (12 • C and −12 • C in the first and second surveys, respectively). During the last survey in 2021, the weather condition was similar, with wind at a maximum speed of 3 m/s, and no precipitation, cloud cover, or waves were identified. The temperature was about 15 • C.
Obtaining aerial data required a perfect match to atmospheric conditions, while determining ground control points for the georeferencing process required careful consideration. To analyze and compare orthophotos from consecutive years, it was necessary to find Remote Sens. 2022, 14, 458 5 of 18 common points on all developed mosaics. For this purpose, using a Trimble RTK (real time kinematic) receiver, fixed GCP points were determined along the wooden spurs ( Figure 4). Obtaining aerial data required a perfect match to atmospheric conditions, while determining ground control points for the georeferencing process required careful consideration. To analyze and compare orthophotos from consecutive years, it was necessary to find common points on all developed mosaics. For this purpose, using a Trimble RTK (real time kinematic) receiver, fixed GCP points were determined along the wooden spurs ( Figure 4). In addition, vertical warning signs (yellow-painted metal pipes 20 cm in diameter with a cross on top), located in the axis of the entire structure and protruding about 1.5 m above the water surface, were used to align the orthophotos in the post-production process. The coordinates of these warning signs were measured with a Leica TS11 total station and the distance to the last wooden piles in the rows of piling spurs located along the entire structure was checked.
All orthophotos were compiled from the vertical images taken by UAV and using Agisoft's Metashape software. This software is used for automatic image processing and creates a dense point cloud using image matching and image georeferencing algorithms ( Figure 5). Alignment and matching of orthophotos from all measurements were performed using QGIS software. All orthophotos from 2016, 2018, and 2021 were cut into 5 sectors, presenting only the modules that contained concrete circles. The rest of the structure made of stones, located between these modules, was omitted from further analyses ( Figure 6). A preliminary validation process was also performed in QGIS software, using manual determination of the centers of all circles visible underwater. In addition, vertical warning signs (yellow-painted metal pipes 20 cm in diameter with a cross on top), located in the axis of the entire structure and protruding about 1.5 m above the water surface, were used to align the orthophotos in the post-production process. The coordinates of these warning signs were measured with a Leica TS11 total station and the distance to the last wooden piles in the rows of piling spurs located along the entire structure was checked.
All orthophotos were compiled from the vertical images taken by UAV and using Agisoft's Metashape software. This software is used for automatic image processing and creates a dense point cloud using image matching and image georeferencing algorithms ( Figure 5). Obtaining aerial data required a perfect match to atmospheric conditions, while determining ground control points for the georeferencing process required careful consideration. To analyze and compare orthophotos from consecutive years, it was necessary to find common points on all developed mosaics. For this purpose, using a Trimble RTK (real time kinematic) receiver, fixed GCP points were determined along the wooden spurs ( Figure 4). In addition, vertical warning signs (yellow-painted metal pipes 20 cm in diameter with a cross on top), located in the axis of the entire structure and protruding about 1.5 m above the water surface, were used to align the orthophotos in the post-production process. The coordinates of these warning signs were measured with a Leica TS11 total station and the distance to the last wooden piles in the rows of piling spurs located along the entire structure was checked.
All orthophotos were compiled from the vertical images taken by UAV and using Agisoft's Metashape software. This software is used for automatic image processing and creates a dense point cloud using image matching and image georeferencing algorithms ( Figure 5). Alignment and matching of orthophotos from all measurements were performed using QGIS software. All orthophotos from 2016, 2018, and 2021 were cut into 5 sectors, presenting only the modules that contained concrete circles. The rest of the structure made of stones, located between these modules, was omitted from further analyses ( Figure 6). A preliminary validation process was also performed in QGIS software, using manual determination of the centers of all circles visible underwater. Alignment and matching of orthophotos from all measurements were performed using QGIS software. All orthophotos from 2016, 2018, and 2021 were cut into 5 sectors, presenting only the modules that contained concrete circles. The rest of the structure made of stones, located between these modules, was omitted from further analyses ( Figure 6). A preliminary validation process was also performed in QGIS software, using manual determination of the centers of all circles visible underwater. Remote Sens. 2022, 14, x FOR PEER REVIEW 6 of 18

Methodology and Results
Owing to artificial reef construction, which includes a stone mound in addition to the habitat modules, further analysis included only the five groups of habitat modules, without considering the stone mound separating these modules in further work. The concrete blocks were arranged in fairly equal rows (the 1st module on the west side consists of 2 segments in which the blocks are arranged in 5 and 6 rows, the middle 3 modules have 9 rows each, and the eastern last module has 2 segments containing 8 and 9 rows). In several rows, concrete blocks were broken or overturned during the construction phase in 2016 and these elements were excluded from further analysis at the start. Only those elements where a circular shape could be recognized were included in the study. In 2018 and 2021, structural damage and displacement of the circles in each row are evident.
Each concrete block had a diameter of 2.4 m and different heights (i.e., 1.5, 2, and 2.5 m) depending on the depth of the bottom. The walls of the blocks were 15 cm thick and comprised 12 conical holes distributed evenly along the circumference and one conical hole on the upper surface. The characteristic shape and construction enabled the implementation of an automatic algorithm to detect submerged concrete circle-shaped elements in the UAV images acquired during the monitoring campaigns.
The proposed solution was based on image processing and analysis techniques. A four-step algorithm ( Figure 7) was proposed, and an object detection methodology provided a list of circle-shaped origins stored in a floating point format (coordinates in the image plane) for each analyzed series of UAV images. The algorithm was implemented and validated in Python using the OpenCV, NumPy, Pillow, and Matplotlib standard libraries.

Methodology and Results
Owing to artificial reef construction, which includes a stone mound in addition to the habitat modules, further analysis included only the five groups of habitat modules, without considering the stone mound separating these modules in further work. The concrete blocks were arranged in fairly equal rows (the 1st module on the west side consists of 2 segments in which the blocks are arranged in 5 and 6 rows, the middle 3 modules have 9 rows each, and the eastern last module has 2 segments containing 8 and 9 rows). In several rows, concrete blocks were broken or overturned during the construction phase in 2016 and these elements were excluded from further analysis at the start. Only those elements where a circular shape could be recognized were included in the study. In 2018 and 2021, structural damage and displacement of the circles in each row are evident.
Each concrete block had a diameter of 2.4 m and different heights (i.e., 1.5, 2, and 2.5 m) depending on the depth of the bottom. The walls of the blocks were 15 cm thick and comprised 12 conical holes distributed evenly along the circumference and one conical hole on the upper surface. The characteristic shape and construction enabled the implementation of an automatic algorithm to detect submerged concrete circle-shaped elements in the UAV images acquired during the monitoring campaigns.
The proposed solution was based on image processing and analysis techniques. A fourstep algorithm (Figure 7) was proposed, and an object detection methodology provided a list of circle-shaped origins stored in a floating point format (coordinates in the image plane) for each analyzed series of UAV images. The algorithm was implemented and validated in Python using the OpenCV, NumPy, Pillow, and Matplotlib standard libraries.

Image Processing
The original input images were delivered from UAV surveys in RGB color format with a possibly high spatial resolution. Because the color characteristics of individual images differed significantly (due to various environmental conditions, such as lighting, time of day, and season), several color-based pre-processing activities were required. Thus, each input image was converted into five dedicated representations: grayscale, CIE Lab, HSV, YCbCr, and CMYK. In this study, different channels from each of the five representations were tested and observed to provide a higher accuracy than that of a perpixel-based approach [44]. Finally, the five most promising channels were chosen for further analysis: cyan from CMYK, Cr from YCbCr, hue from HSV, a channel from CIELab, and a grayscale image ( Figure 8).

Image Processing
The original input images were delivered from UAV surveys in RGB color format with a possibly high spatial resolution. Because the color characteristics of individual images differed significantly (due to various environmental conditions, such as lighting, time of day, and season), several color-based pre-processing activities were required. Thus, each input image was converted into five dedicated representations: grayscale, CIE Lab, HSV, YCbCr, and CMYK. In this study, different channels from each of the five representations were tested and observed to provide a higher accuracy than that of a per-pixel-based approach [44]. Finally, the five most promising channels were chosen for further analysis: cyan from CMYK, Cr from YCbCr, hue from HSV, a channel from CIELab, and a grayscale image ( Figure 8).

Image Processing
The original input images were delivered from UAV surveys in RGB color format with a possibly high spatial resolution. Because the color characteristics of individual images differed significantly (due to various environmental conditions, such as lighting, time of day, and season), several color-based pre-processing activities were required. Thus, each input image was converted into five dedicated representations: grayscale, CIE Lab, HSV, YCbCr, and CMYK. In this study, different channels from each of the five representations were tested and observed to provide a higher accuracy than that of a perpixel-based approach [44]. Finally, the five most promising channels were chosen for further analysis: cyan from CMYK, Cr from YCbCr, hue from HSV, a channel from CIELab, and a grayscale image ( Figure 8).

Circle Detection Pipeline
The second step, which performs the detection of the circular elements of underwater construction, is the core of this approach. As the objects we are looking for are circles, we have proposed several methods to perform automatic circle recognition for each singlechannel image [45]. First, we normalized all the extracted matrices to a 0-255 range and stored them in an unsigned int 8-bit format. Second, contrast-limited adaptive histogram equalization (CLAHE) [46] was applied to all the normalized images, with the following parameters: clip limit of 2.0 and tile grid size of 8 × 8 pixels. In each image, the detection of circular shapes with Hough transform (Figure 9) [47] was realized by assuming a very narrow (pre-defined) range of radii. For each individual list of detected candidates saved in the earlier step, we integrated circles that were closer than the distance equal to the radius by averaging their positions. Finally, we concatenated the five resulting lists (one for each pre-defined channel) and removed duplicates ( Figure 9) (candidates that were close to each other) by using radius as a criterion.

Circle Detection Pipeline
The second step, which performs the detection of the circular elements of underwater construction, is the core of this approach. As the objects we are looking for are circles, we have proposed several methods to perform automatic circle recognition for each singlechannel image [45]. First, we normalized all the extracted matrices to a 0-255 range and stored them in an unsigned int 8-bit format. Second, contrast-limited adaptive histogram equalization (CLAHE) [46] was applied to all the normalized images, with the following parameters: clip limit of 2.0 and tile grid size of 8 × 8 pixels. In each image, the detection of circular shapes with Hough transform (Figure 9) [47] was realized by assuming a very narrow (pre-defined) range of radii. For each individual list of detected candidates saved in the earlier step, we integrated circles that were closer than the distance equal to the radius by averaging their positions. Finally, we concatenated the five resulting lists (one for each pre-defined channel) and removed duplicates (Figure 9) (candidates that were close to each other) by using radius as a criterion.

Hole Detection Pipeline
In the next stage, we performed hole detection for the grayscale, cyan, and hue channels ( Figure 10). In this process, we detected the centers of circles that were significantly darker than their surroundings; the images were normalized to 0-255 and 8-bit representations and histogram equalization was performed according to global equalization criteria. The prepared images were segmented using adaptive thresholding. However, the images had to

Hole Detection Pipeline
In the next stage, we performed hole detection for the grayscale, cyan, and hue channels ( Figure 10).

Circle Detection Pipeline
The second step, which performs the detection of the circular elements of underwater construction, is the core of this approach. As the objects we are looking for are circles, we have proposed several methods to perform automatic circle recognition for each singlechannel image [45]. First, we normalized all the extracted matrices to a 0-255 range and stored them in an unsigned int 8-bit format. Second, contrast-limited adaptive histogram equalization (CLAHE) [46] was applied to all the normalized images, with the following parameters: clip limit of 2.0 and tile grid size of 8 × 8 pixels. In each image, the detection of circular shapes with Hough transform (Figure 9) [47] was realized by assuming a very narrow (pre-defined) range of radii. For each individual list of detected candidates saved in the earlier step, we integrated circles that were closer than the distance equal to the radius by averaging their positions. Finally, we concatenated the five resulting lists (one for each pre-defined channel) and removed duplicates (Figure 9) (candidates that were close to each other) by using radius as a criterion.

Hole Detection Pipeline
In the next stage, we performed hole detection for the grayscale, cyan, and hue channels ( Figure 10). In this process, we detected the centers of circles that were significantly darker than their surroundings; the images were normalized to 0-255 and 8-bit representations and histogram equalization was performed according to global equalization criteria. The prepared images were segmented using adaptive thresholding. However, the images had to In this process, we detected the centers of circles that were significantly darker than their surroundings; the images were normalized to 0-255 and 8-bit representations and histogram equalization was performed according to global equalization criteria. The prepared images were segmented using adaptive thresholding. However, the images had to be blurred with a Gaussian mask of 15 × 15 pixels to smooth small objects that interrupted the threshold process. Prepared images were analyzed using a sequence of morphological opening and closing operations, with an ellipse of 5 × 5 pixels being used as a structuring element. This allowed us to find isolated contours (blobs) that satisfy the following conditions: a blob area with a larger range (600, 1400) and an aspect ratio of width to height in a specific range (0.8, 1.2). For each individual list of candidates, we integrated holes (candidates) that were closer than the distance equal to the radius by averaging their positions and removing duplicates (holes close to each other were removed by using radius as a criterion: Figure 11-left side).
Remote Sens. 2022, 14, x FOR PEER REVIEW 9 of 18 be blurred with a Gaussian mask of 15 × 15 pixels to smooth small objects that interrupted the threshold process. Prepared images were analyzed using a sequence of morphological opening and closing operations, with an ellipse of 5 × 5 pixels being used as a structuring element. This allowed us to find isolated contours (blobs) that satisfy the following conditions: a blob area with a larger range (600, 1400) and an aspect ratio of width to height in a specific range (0.8, 1.2). For each individual list of candidates, we integrated holes (candidates) that were closer than the distance equal to the radius by averaging their positions and removing duplicates (holes close to each other were removed by using radius as a criterion: Figure 11-left side). Figure 11. Detection of circle centers (left) and detected objects of specific geometrical properties in respective color channels: green objects-gray channel; blue objects-cyan channel; red objectshue channel.

Concatenation of the Results
The last analysis step included the concatenation of the results from the processing pipelines. First, the duplicates were removed by averaging the origins. Second, isolated candidates far from their neighbors were removed. We considered the distance between origins (based on the experiments, a value of 200 pixels was used). After removing all the duplicates and isolated candidates, we generated a list of circle origins for all test campaigns.

Validation Process
In the process of validating the algorithm, two qualified operators (specialists in the interpretation of remote sensing data of the coastal zone) independently identified the origin of the circles visible at the selected test site. In addition, a specially developed webbased program, based on the Leaflet Javascript library and accessible via a web browser, was used to benchmark the location recognition of individual infrastructure elements. This program works based on an algorithm and manual classification. A portal was created to allow multiple users (mostly senior students from various marine, geography, and computer science majors) to simultaneously enter identification points for the center of concrete circles. Each user, after registering, was given a random set of images with habitat modules from different years, and their task was to mark the centers of all circles visible in each image ( Figure 12A). In this way, results were collected from approximately 40 additional operators who collectively marked over 25,000 points ( Figure 12C). The resulting point coordinates were then stored in a PostgreSQL database ( Figure 12B). The database allowed the rejection of significant outliers and averaged the data to centroids in the WGS84 system (EPSG:4326). Figure 11. Detection of circle centers (left) and detected objects of specific geometrical properties in respective color channels: green objects-gray channel; blue objects-cyan channel; red objectshue channel.

Concatenation of the Results
The last analysis step included the concatenation of the results from the processing pipelines. First, the duplicates were removed by averaging the origins. Second, isolated candidates far from their neighbors were removed. We considered the distance between origins (based on the experiments, a value of 200 pixels was used). After removing all the duplicates and isolated candidates, we generated a list of circle origins for all test campaigns.

Validation Process
In the process of validating the algorithm, two qualified operators (specialists in the interpretation of remote sensing data of the coastal zone) independently identified the origin of the circles visible at the selected test site. In addition, a specially developed webbased program, based on the Leaflet Javascript library and accessible via a web browser, was used to benchmark the location recognition of individual infrastructure elements. This program works based on an algorithm and manual classification. A portal was created to allow multiple users (mostly senior students from various marine, geography, and computer science majors) to simultaneously enter identification points for the center of concrete circles. Each user, after registering, was given a random set of images with habitat modules from different years, and their task was to mark the centers of all circles visible in each image ( Figure 12A). In this way, results were collected from approximately 40 additional operators who collectively marked over 25,000 points ( Figure 12C). The resulting point coordinates were then stored in a PostgreSQL database ( Figure 12B). The database allowed the rejection of significant outliers and averaged the data to centroids in the WGS84 system (EPSG:4326). In the next step, the obtained data was transformed into a cartesian grid format corresponding to the pixel layout of the individual images in the X and Y axes system. The upper-left corner of the raster was taken as point 0,0. The program was written in the Python language. For analysis, the identified circle origins were averaged and taken as reference points, which were compared with the output of the origins obtained using the automatic method.
The origins of circles undetected by our algorithm were not considered during the accuracy estimation stage. The comparison utilized several stages. For each point in the reference list, the closest point in the list was searched for based on the coordinates in the image plane. The Euclidean distance between the points was calculated and stored as an error value. However, if several points were found within that distance, the closest one was chosen. Finally, if no points were found in the predefined tolerance area, the algorithm skipped to the next point in the reference list.
The study of variation in the position of points (geometric centers of circles) was based on the rejection of constancy of the position of any of the points. In this method, it was assumed that any number of test points could move in time in any direction. However, it was assumed that there would be a certain number of fixed points, and the constancy of these points could be confirmed by computational methods. At the same time, the variability of the position of the moving points will be determined in relation to the position of the points considered as a fixed. Therefore, research questions arose, as follows: How could we find fixed points? How could we determine the constancy of these points (reference datum)? How could we determine the displacements of the moving In the next step, the obtained data was transformed into a cartesian grid format corresponding to the pixel layout of the individual images in the X and Y axes system. The upper-left corner of the raster was taken as point 0,0. The program was written in the Python language. For analysis, the identified circle origins were averaged and taken as reference points, which were compared with the output of the origins obtained using the automatic method.
The origins of circles undetected by our algorithm were not considered during the accuracy estimation stage. The comparison utilized several stages. For each point in the reference list, the closest point in the list was searched for based on the coordinates in the image plane. The Euclidean distance between the points was calculated and stored as an error value. However, if several points were found within that distance, the closest one was chosen. Finally, if no points were found in the predefined tolerance area, the algorithm skipped to the next point in the reference list.
The study of variation in the position of points (geometric centers of circles) was based on the rejection of constancy of the position of any of the points. In this method, it was assumed that any number of test points could move in time in any direction. However, it was assumed that there would be a certain number of fixed points, and the constancy of these points could be confirmed by computational methods. At the same time, the variability of the position of the moving points will be determined in relation to the position of the points considered as a fixed. Therefore, research questions arose, as follows: How could we find fixed points? How could we determine the constancy of these points (reference datum)? How could we determine the displacements of the moving points relative to the fixed points? A method for answering these questions using coordinate transformation was proposed, as shown in Figure 13.
In the first stage of the 2D transformation (Figure 13), the criterion of points stability (reference datum) was determined. This criterion is defined as the uncertainty (mean error) of the point position after the transformation, not exceeding the double value of the standard deviation, which is described by Formula (1). This approach is applicable to the study of the deformation of engineering structures. The coordinates of the outlier points were not used to determine the coordinate transformation parameters, and these points were treated as moving and the displacement (deformation) values are determined for them [48,49].
Criterion (1) was applied iteratively, with the rejection of points that exceeded twice the standard uncertainty, until the set of points stabilized, i.e., the position uncertainties of all points met condition (1). These points were considered fixed, and the displacements of the moving points were determined relative to these points. For the transformation of coordinates, a conformal (Helemrt) transformation model with a center of gravity, in the form of a Q-ST (quasi similarity transformation), was adopted in Relation (2) [50], as follows: where: • -coordinates in the target reference frame (TRF); • -(quasi) rotation matrix; • -scale factor; • -coordinates in the transformed reference frame (trf); • -coordinates of the center of gravity in the transformed reference frame (centroid of scaled transformed reference frame); • -coordinates of the center of gravity in the target reference system (centroid of target reference frame); • -vector length in the target reference frame (TRF); In the first stage of the 2D transformation (Figure 13), the criterion of points stability (reference datum) was determined. This criterion is defined as the uncertainty (mean error) of the point position after the transformation, not exceeding the double value of the standard deviation, which is described by Formula (1). This approach is applicable to the study of the deformation of engineering structures. The coordinates of the outlier points were not used to determine the coordinate transformation parameters, and these points were treated as moving and the displacement (deformation) values are determined for them [48,49].
Criterion (1) was applied iteratively, with the rejection of points that exceeded twice the standard uncertainty, until the set of points stabilized, i.e., the position uncertainties of all points met condition (1). These points were considered fixed, and the displacements of the moving points were determined relative to these points. For the transformation of coordinates, a conformal (Helemrt) transformation model with a center of gravity, in the form of a Q-ST (quasi similarity transformation), was adopted in Relation (2) [50], as follows: where: • X TRF -coordinates in the target reference frame (TRF); • Q -(quasi) rotation matrix; • λ-scale factor; • X tr f -coordinates in the transformed reference frame (trf); • X ctr f -coordinates of the center of gravity in the transformed reference frame (centroid of scaled transformed reference frame); • X CTRF -coordinates of the center of gravity in the target reference system (centroid of target reference frame); • D TRF -vector length in the target reference frame (TRF); • D tr f -vector length in the transformed reference frame (trf).
The rotation matrix Q (quasi) in Relation (2), in 2D space, is given in Matrix (3), as follows: The coefficients of Matrix (3) were determined using the least squares method, according to Formula (4), as follows: where As a result of the iterative approach in the transformation in Relation (2), fixed points were selected to satisfy Condition (1). In the next step, the transformation in Relation (2) was performed for the moving points, that is, those that did not participate in the determination of the transformation coefficients. The mean errors after the Q-ST transformation in Relation (2), for the fixed points and the displacements for the moving points, are summarized in Table 1. The data presented in Table 1 show that the average transformation error on points determined as fixed is within the values of determination of geometric centers of circles (Figure 14), which indicates high stability of the position of points-the constancy of the reference system. At the same time, it can be noted that the average displacement of points in the period 2018-2021 for moving points is 34 pix, which for a ground sampling distance (GSD) of 2 cm is about 0.7 m.
The overall detection accuracy (understood as a number of correctly detected objects in comparison to the actual ones) is equal to 0.987, varying for particular years (see Figure 13, left chart). According to the cumulative histogram of error distributions (Figure 13, right chart), the detection accuracy is quite high. The processing of all channels ("all") instead of individual ones ("c"-cyan; "gr"-grayscale; "h"-HSV hue; "cr"-YCbCr chroma; "a"-Lab) gives the highest detection accuracy. The errors distribution for the cr channel is slightly better than for all the concatenated channels. Detection accuracy for all the channels is higher, which makes such an approach more sensible. The majority of the errors did not exceed 4 pixels, with an average of 4.06 and a median of 4.12 pixels. Pixel value has been presented here as a field pixel representation-GSD (ground sampling distance)-a value of 0.02 m. This is primarily because the operators marked the circle origins by considering the shadow inside the central part of the objects, whereas the algorithm detected the center on the two-dimensional projection of the three-dimensional structure. It can be observed on the spatial distribution plot (Figure 15), where the majority of errors are within offset (2.2). The closer look unveils that the accuracy for images taken in 2016 and 2018 is quite high. Moreover, the errors are kept in significantly low range. On the other hand, images taken in 2021 characterize lower accuracy and larger spread of errors distribution, since in these images have lower contrast coming from the fact that the underwater structures have been covered by sand. The data presented in Table 1 show that the average transformation error on points determined as fixed is within the values of determination of geometric centers of circles (Figure 14), which indicates high stability of the position of points-the constancy of the reference system. At the same time, it can be noted that the average displacement of points in the period 2018-2021 for moving points is 34 pix, which for a ground sampling distance (GSD) of 2 cm is about 0.7 m. The overall detection accuracy (understood as a number of correctly detected objects in comparison to the actual ones) is equal to 0.987, varying for particular years (see Figure  13, left chart). According to the cumulative histogram of error distributions (Figure 13, right chart), the detection accuracy is quite high. The processing of all channels ("all") instead of individual ones ("c"-cyan; "gr"-grayscale; "h"-HSV hue; "cr"-YCbCr chroma; "a"-Lab) gives the highest detection accuracy. The errors distribution for the cr channel is slightly better than for all the concatenated channels. Detection accuracy for all the channels is higher, which makes such an approach more sensible. The majority of the errors did not exceed 4 pixels, with an average of 4.06 and a median of 4.12 pixels. Pixel value has been presented here as a field pixel representation-GSD (ground sampling distance)-a value of 0,02 m. This is primarily because the operators marked the circle origins by considering the shadow inside the central part of the objects, whereas the algorithm detected the center on the two-dimensional projection of the three-dimensional structure. It can be observed on the spatial distribution plot (Figure 15), where the majority of errors are within offset (2.2). The closer look unveils that the accuracy for images taken in 2016 and 2018 is quite high. Moreover, the errors are kept in significantly low range. On the other hand, images taken in 2021 characterize lower accuracy and larger spread of errors distribution, since in these images have lower contrast coming from the fact that the underwater structures have been covered by sand. The additional verification of the employed approach involved the same set of experiments, yet with a focus on elementary algorithm stages (i.e., circle detection and hole detection only stages). We have calculated accuracy (as a percentage of detected objects in comparison to the actual ones) and error distribution as a distance between detected object position and the actual one. The detailed results are presented in Figures 16 and 17. It is worth noticing that the detecting holes in the hue channel for images taken in 2016 gave a very low accuracy, while for images taken in 2018 and 2021 it worked better. On the other hand, detecting holes in cyan channel for images taken in 2021 was also not efficient. As it can been seen, an approach founded on a single detection stage (no matter which one) cannot result in high accuracy and low detection errors. The simple solution, that makes the algorithm insensitive to color characteristics, is concatenation of the results from different channels. The additional verification of the employed approach involved the same set of experiments, yet with a focus on elementary algorithm stages (i.e., circle detection and hole detection only stages). We have calculated accuracy (as a percentage of detected objects in comparison to the actual ones) and error distribution as a distance between detected object position and the actual one. The detailed results are presented in Figures 16 and 17. It is worth noticing that the detecting holes in the hue channel for images taken in 2016 gave a very low accuracy, while for images taken in 2018 and 2021 it worked better. On the other hand, detecting holes in cyan channel for images taken in 2021 was also not efficient. As it can been seen, an approach founded on a single detection stage (no matter which one) cannot result in high accuracy and low detection errors. The simple solution, that makes the algorithm insensitive to color characteristics, is concatenation of the results from different channels.

Discussion
In this paper, we present an innovative work, where an off-the-shelf UAV with a regular RGB camera was used to monitor shallow hybrid coastal protection measures. Moreover, the remarkable potential of UAV image-based analysis has been proven by the implementation of bathymetry estimation methods [3,36,51] and monitoring techniques for submerged objects [37,52,53]; however, none of the presented research have been used to identify underwater coastal protection measures. They are focused on military purposes [52] or water organisms and waste [52,53]. Furthermore, the employed image processing methods have been successfully applied to circular objects detection in other fields of science [54,55]; hence, the application of Hough transform and related image processing-based methods to analyze underwater objects is clearly a new approach.
The principal part of the analyzed underwater artificial reef in our experiment contained multiple regular circle textures. Based on these elements, the implemented algorithm was optimized and customized to provide good identification results with the proposed methodology. Inevitably, some elements require more attention; however, the sections identified according to the construction documentation are crucial for reef stability. Circular concrete blocks are the most exposed to hydrometeorological conditions; thus, monitoring these elements is obligatory for construction efficiency.
The proposed methodology was used for five different color separation stages in order to be invariant to as many factors as possible. Nevertheless, several thresholds to analyze the data have been implemented. The algorithm was tuned to specific characteristics of input data, while using thresholds and predefined values that can be easily adopted; this enables the proposed method to be used in different applications and in various environmental conditions. The pixel size of circles or holes depends on their physical dimensions and is related to the imaging altitude and spatial resolution of images. It primarily

Discussion
In this paper, we present an innovative work, where an off-the-shelf UAV with a regular RGB camera was used to monitor shallow hybrid coastal protection measures. Moreover, the remarkable potential of UAV image-based analysis has been proven by the implementation of bathymetry estimation methods [3,36,51] and monitoring techniques for submerged objects [37,52,53]; however, none of the presented research have been used to identify underwater coastal protection measures. They are focused on military purposes [52] or water organisms and waste [52,53]. Furthermore, the employed image processing methods have been successfully applied to circular objects detection in other fields of science [54,55]; hence, the application of Hough transform and related image processing-based methods to analyze underwater objects is clearly a new approach.
The principal part of the analyzed underwater artificial reef in our experiment contained multiple regular circle textures. Based on these elements, the implemented algorithm was optimized and customized to provide good identification results with the proposed methodology. Inevitably, some elements require more attention; however, the sections identified according to the construction documentation are crucial for reef stability. Circular concrete blocks are the most exposed to hydrometeorological conditions; thus, monitoring these elements is obligatory for construction efficiency.
The proposed methodology was used for five different color separation stages in order to be invariant to as many factors as possible. Nevertheless, several thresholds to analyze the data have been implemented. The algorithm was tuned to specific characteristics of input data, while using thresholds and predefined values that can be easily adopted; this enables the proposed method to be used in different applications and in various environmental conditions. The pixel size of circles or holes depends on their physical dimensions and is related to the imaging altitude and spatial resolution of images. It primarily

Discussion
In this paper, we present an innovative work, where an off-the-shelf UAV with a regular RGB camera was used to monitor shallow hybrid coastal protection measures. Moreover, the remarkable potential of UAV image-based analysis has been proven by the implementation of bathymetry estimation methods [3,36,51] and monitoring techniques for submerged objects [37,52,53]; however, none of the presented research have been used to identify underwater coastal protection measures. They are focused on military purposes [52] or water organisms and waste [52,53]. Furthermore, the employed image processing methods have been successfully applied to circular objects detection in other fields of science [54,55]; hence, the application of Hough transform and related image processing-based methods to analyze underwater objects is clearly a new approach.
The principal part of the analyzed underwater artificial reef in our experiment contained multiple regular circle textures. Based on these elements, the implemented algorithm was optimized and customized to provide good identification results with the proposed methodology. Inevitably, some elements require more attention; however, the sections identified according to the construction documentation are crucial for reef stability. Circular concrete blocks are the most exposed to hydrometeorological conditions; thus, monitoring these elements is obligatory for construction efficiency.
The proposed methodology was used for five different color separation stages in order to be invariant to as many factors as possible. Nevertheless, several thresholds to analyze the data have been implemented. The algorithm was tuned to specific characteristics of input data, while using thresholds and predefined values that can be easily adopted; this enables the proposed method to be used in different applications and in various environmental conditions. The pixel size of circles or holes depends on their physical dimensions and is related to the imaging altitude and spatial resolution of images. It primarily influences the averaging filter mask, morphological element size, and blob area. The thresholds used at the intensity segmentation stages are universal to some extent, since they have been calculated for several images acquired at different times and various seasons.
UAV, combined with GNSS measurements, has been shown in many works to be an efficient tool for high-resolution analysis in coastal zones [56][57][58]. It has been proved that UAV photogrammetry is useful when decimeter-level accuracy is required. Different surveys report photogrammetric model accuracy values on sub-centimeter values, respectively, for the horizontal component [59,60]. The accuracy we found in identification of underwater objects by UAV was very similar. The majority of errors did not exceed 4 pixels (8 cm) for the 3 surveys. Such a result is very similar to that found in [33,59,61,62], even though, in those surveys, UAV was used to identified elements on land instead of under water.
The realized experimentation demonstrates that the UAV is a suitable tool to monitor underwater structures withing accepted range of the error propagation. However, this is only valid if the accuracy is assessed and monitored throughout the entire process. Notably, both the circle and hole detection process documented here, and the considerations presented in this manuscript, were made possible because the weather conditions were thoroughly analyzed and tracked throughout the entire process. Although this research feature, in particular, made it possible to evaluate the usage of UAV monitoring of underwater structures, it is suggested that uncertainty is taken into account while processing UAV surveys for coastal monitoring of underwater objects.
The results reveal that the movement complexity of the identified objects is visible over time. This leads to a need for additional monitoring procedures and analyses of changes in reef construction.

Conclusions
This manuscript proposed a monitoring solution for underwater coastal protection measures by automated identification of crucial construction elements in RGB UAV photographs, which was performed using image processing algorithms. We applied two different pipelines and concatenated the outcomes. The experimental results showed that the accuracy of our algorithm meets the requirements of the monitoring process and can be used with readily available UAV equipment to save manpower and material costs. Furthermore, we will continue to conduct an analysis of various underwater protection measures to extend this research to include more application scenarios. Additionally, differences in the identified objects, which have been observed using UAV image time series, will be analyzed and compared against hydrometeorological conditions to determine the most significant parameters that influence the protection measures in coastal zones.
The UAV monitoring campaign turned out to be suitably accurate for studying coastal defense underwater structures, for two reasons especially: firstly, the practicality of these devices-this is due to their wide availability and operating velocity, especially in comparison with traditional expensive methods of mapping underwater objects with echosounders. Secondly, spectral RGB data represents a huge advantage when conducting 2D detection of underwater objects and their relationship to topography.
The realized UAV surveys in the underwater environment experience different types of limitations. They are represented by the complexity and, at a certain level, the impossibility of complete detection of single elements of analyzed construction, which are partly covered with sand. This factor strongly complicates the detection process, and final results can be assumed to be realistic when elements such as sand or vegetation have little or no influence.
Despite this, the fully automatic object detection algorithm presented here, which is able to work on RGB data, represents a rapid and useful way to monitor the underwater coastal structures. The high resolution of UAV data allowed rapid elaborations that, until today, were much more time consuming, and only in the best-case scenario (multibeam echosounders) did these attain a similar resolution. Here, we confirm the conclusions made by Fabbri et al., 2021, that the availability of small, high-resolution sensors that UAVs can transport, such as RGB cameras or multi-spectral cameras, is increasing, presenting a new approach for high-resolution spatial studies. Thus, the implementation of the presented methodology can deliver data on underwater coastal protection on demand, measuring comparably high spatial, and much higher temporal, resolutions than traditional measurements. These attributes are particularly useful for a range of coastal applications, such as monitoring and other activities within integrated coastal zone management.
Funding: This research was funded in part by the National Science Centre, Poland, through the project "Influence of surfzone and beach morphology on coastal cliff retreat-INSUMOR", project number 2020/39/B/ST10/01122.