Automatic tree detection and attribute characterization using portable terrestrial lidar

This study details a methodology to automatically detect the positions of and dasometric information about individual Eucalyptus trees from a point cloud acquired with a portable LiDAR system. Currently, the implementation of portable laser scanners (PLS) in forest inventories is being studied, since they allow for significantly reduced field-work time and costs when compared to the traditional inventory methods and other LiDAR systems. However, it has been shown that their operability and efficiency are dependent upon the species assessed, and therefore, there is a need for more research assessing different types of stands and species. Additionally, a few studies have been conducted in Eucalyptus stands, one of the tree genus that is most commonly planted around the world. In this study, a PLS system was tested in a Eucalyptus globulus stand to obtain different metrics of individual trees. An automatic methodology to obtain inventory data (individual tree positions, DBH, diameter at different heights, and height of individual trees) was developed using public domain software. The results were compared to results obtained with a static terrestrial laser scanner (TLS). The methodology was able to identify 100% of the trees present in the stand in both the PLS and TLS point clouds. For the PLS point cloud, the RMSE of the DBH obtained was 0.0716, and for the TLS point cloud, it was 0.176. The RMSE for height for the PLS point cloud was 3.415 m, while for the PLS point cloud, it was 10.712 m. This study demonstrates the applicability of PLS systems for the estimation of the metrics of individual trees in adult Eucalyptus globulus stands.


Introduction
Forest inventories constitute an essential tool for sustainable forest management, since they provide relevant information from an ecological and economical perspective (Fankhauser et al. 2018). Forest inventories are usually based on the measurement of tree attributes, such as height, diameter at breast height (DBH), stem profiles, and first branch height, among others (Maas et al. 2008). These forest attributes are traditionally obtained through fieldwork using measuring instruments, such as calipers, hypsometers, and tape measures (Piermattei et al. 2019). However, conducting research in this way usually entails a great investment in terms of economic resources and time (Dalla Corte et al. 2020), thus limiting the periodicity, scale, and scope of these inventories (Liang et al. 2018).
In recent years, advances in remote-sensing technology have provided solutions and aided in overcoming the Communicated by Roetzer. drawbacks of traditional approaches to information acquisition for forest inventories; these advances have led to the enhancement of the provided inventory data (White et al. 2016). LiDAR technology is especially useful in this field. It provides 3D representations of target areas and allows for the characterization of the structure of the above ground biomass (Lechner et al. 2020;Shugart et al. 2010).
LiDAR devices can be mounted on different platforms, each with its associated strengths and weaknesses: aerial platforms (either manned or unmanned aerial vehicles), ground platforms (such as terrestrial laser scanners), mobile terrestrial platforms, and on-board satellite platforms (Michez et al. 2016). Airborne laser scanning (ALS) systems are the most widespread; they have proven to be very helpful in describing and characterizing forest variables at relatively large scales (Kangas et al. 2018). However, due to the fact that ALS point clouds hardly provide representations of tree trunks, they have limitations when it comes to providing tree diameter information, such as diameter at breast height (DBH) (Chen et al. 2019;Gollob et al. 2020).
Ground-based LiDAR systems allow for the creation of 3D representations of the vegetation structure at the lower canopy level with higher point densities than ALS systems (Donager et al. 2021;Hilker et al. 2010;Xia et al. 2021). The ground-based devices can be divided into terrestrial laser scanners (TLS) and the more-recently-developed portable systems, known as mobile laser scanners (MLS).
TLS are stationary systems that must be fixed on a tripod; they can provide high-quality diameter metrics for multiple tree species in different environments (Brede et al. 2022;Calders et al. 2020;Panagiotidis et al. 2022;Gao and Kan 2022). The main weakness of these systems, in terms of operability, stems from their static nature (Liang et al. 2016;Tremblay and Béland 2018), which generates occlusion problems (tree occlusion, treetop occlusion, or even stem occlusion due to the presence of shrubs or branches) .
MLS systems work mounted on moving platforms (e.g., vehicles) (Del Perugia et al. 2019). When it can be carried by a human, the system is known as a personal laser scanner (PLS), although there are a variety of terms that are used to refer to these systems: handheld mobile laser scanners (HMLS), wearable laser scanners (WLS), and handheld laser scanners (HLS) (Gollob et al. 2020). According to Bauwens et al. (2016), PLS are more suitable for forest environments than other MLS systems with integrated GNSS devices; high-quality GNSS is not usually feasible under forest cover. Furthermore, the other MLS systems cannot guarantee the desired continuous and nondestructive acquisition of data. Additionally, thanks to its portability, PLS systems significantly reduce field work time and costs when compared to the other laser scanning systems; they also improve point cloud density and avoid the occlusion problems that are characteristic of TLS systems Chiappini et al. 2022;Donager et al. 2021).
The applicability of PLS systems for forest inventories is being assessed in the scientific literature. An example of the implementation of PLS into forest inventories can be found in Hyyppä et al. (2020). They developed a methodology to assess stem curve and stem volume using PLS in boreal forests. Other examples are the studies by Chiappini et al. (2022) and Donager et al. (2021). They used PLS systems to detect single trees and obtain stand variables, such as DBH and tree height in Pinus nigra and Pinus ponderosa stands, respectively; both obtaining satisfactory results. Gollob et al. (2020) and Gollob et al. (2021) show a comparison of the efficiency of different ground-based LiDAR devices in the detection of individual trees and the estimation of DBH in different types of temperate-forest stands (i.e., species composition and basal area). They demonstrated the superiority of PLS over TLS. However, Gollob et al. (2020) found that DBH accuracies significantly differed among the tree species analyzed when using the PLS in temperate forests. Bienert et al. (2018) also reported that certain temperateforest tree species might compromise the PLS's operability. Given these observations, and considering the novelty of the technology, Donager et al. (2021) indicated the need for more research in different types of stands and tree species, as well as research on sampling methods and methodological variables, to achieve a better understanding of the capabilities and potential of PLS.
Some studies have focused on the appropriacy of PLS systems for Eucalyptus stands, Eucalyptus being one of the most common genera in forest plantations worldwide (FAO 2005;Messier et al. 2021). Marselis et al. (2016), for instance, used a PLS in a natural Eucalyptus forest to automatically detect individual trees and calculate their DBH, obtaining high-accuracy results. Conversely, Levick et al. (2021) reported difficulties in deriving Eucalyptus tree diameter metrics using PLS in the Australian tropical savanna. Additionally, Camaretta et al. (2021) declared that for their study site (Eucalyptus trees of less than 10 m in height), it was not possible to directly calculate DBH using 3D modeling techniques similar to the ones used in the previously mentioned studies; instead, they used allometric equations along with the height and the crown diameters calculated from the LiDAR point cloud. While these studies constitute valuable examples, research on PLS to capture Eucalyptus stands can, nevertheless, be considered scarce. Additionally, all of the studies mentioned were carried out in Eucalyptus stands in their natural ecosystem. Considering the importance of Eucalyptus trees in forestry, further exploration of PLS capabilities in Eucalyptus stands is necessary. Furthermore, exploration into methodologies incorporating openaccess software would contribute to operational applicability (Zeybek and Vatandaslar 2021).
This study focuses on the capability of a PLS system to detect individual trees and measure tree metrics in a Eucalyptus stand. Open-access-software-tools were used to develop the methodology. In particular, a Eucalyptus globulus plot was captured using an SLAM-based PLS system to obtain inventory data (individual tree positions, DBH, diameter at different heights, and height of individual trees). The results obtained were compared, in terms of accuracy and efficiency in data collection, to the results obtained using a TLS system and to the results obtained using traditional field methods.

Study area
The present study was performed in a Eucalyptus globulus stand, located on Illa do Covo, an Island in Pontevedra, Galicia (Northwestern Spain); see Fig. 1. The island has an area of 7.6 ha. The stand is free of understory and is composed of large individuals with a maximum height of approximately 50 m. A group of representative trees was selected. These individuals covered an area of 1500 m 2 . Figure 2 shows a

Materials
The TLS data were collected using an FARO ® Focus3D 120 laser scanner. This scanner weighs 5 kg, has a measurement range of 120 m, an accuracy range of up to 2 mm at 25 m, a 360º horizontal field of view, and a 305º vertical field of view, and captures up to 976,000 points/second. Six reference spheres were used to co-register the data collected.
The PLS data were collected using a GeoSLAM™ ZEB Horizon scanner. The scanner head of this device weighs 1.5 kg, has a measurement range of 100 m, an accuracy range of between 1 and 3 cm, a 270º horizontal field of view, and a 360º vertical field of view, and captures 300,000 points per second. This device uses a simultaneous localization and mapping (SLAM) algorithm to reference the laser distance measurements in three-dimensional space, while the device is in motion, without the need for a global navigation satellite system (GNSS) (Cabo et al. 2018;Gollob et al. 2020).
A diameter tape, a measuring tape, and a Vertex IV hypsometer were used to collect field data using the traditional methods.

Methodology
This section is divided into four subsections: the process of obtaining the field measurements, the scanning, and preprocessing methodologies for each of the two devices, and the process followed to estimate dasometric variables from the point clouds obtained. An overview of the method is presented in Fig. 3.

Acquisition of tree metrics using the traditional methods
All of the trees within the study plot were measured using the traditional tools. The first metric obtained was the diameters at breast height (DBH); measured with the diameter tape. The measurement was performed in the direction orthogonal to the vertical axis of the tree trunks. The second variable was the height of the trees, which was measured using a Vertex IV hypsometer. To do this, the optimum position to view the treetop had to be selected. The trees lying on or near the contour of the plot could be easily observed, but the trees in the plot's interior proved more troublesome as it was difficult to view their treetops given the proximity and height of the neighboring trees.

Point cloud acquisition and post-processing for the terrestrial laser scanner
Given the range of the device and the stand area in this study, two scans were taken with the TLS. This device is stationary, so it was necessary to acquire several scans to optimize the 3D measurement of the trees' understory with minimum occlusions. Most TLS systems can be configured to satisfy the study requirements. In this case, the parameters were defined to strike a balance between maximizing the resolution of the point clouds and minimizing the amount of time invested in field work. They are described in Table 1. The time invested in data collection was also measured.
Once the single point clouds are obtained, they can be combined and transformed into a common coordinate system. Several procedures can be used to do this (Mora et al. 2021). In open environments where there are no regular surfaces, and vegetation may be moved slightly by the wind, the most highly recommended procedure is to use fixed references like spheres, cylinders, or planar targets (Wilkes et al. 2017). In this case, six reference spheres were distributed within the plot and scanned together with the trees in the scene. The co-registering of point clouds was performed using FAROScene software; the result was the generation of one single point cloud.
The whole point cloud was classified into two classes: ground points and non-ground points. This step is needed for further normalization. It was done with the "lasground" tool from the LAStools software. The tool settings were adjusted to obtain an optimized performance.
The normalization of the classified point cloud was performed to convert the altitude of points into their heights above ground. This step allows for the automation of the Canopy Height Model generation. These models represent vegetation heights as raster maps. The normalization algorithm applied consisted of a k-nearest neighborhood (KNN) interpolation with an inverse distance weighting (IDW). It was implemented using the "lidR" package from the Rstudio software.
The final processing step was noise removal. This was done by defining the admissible values for the normalized heights of the points. These threshold values were 0 and 60 m, and all points below the minimum threshold value or above the maximum threshold value were removed from the point cloud.

Point cloud acquisition and post-processing with the portable laser scanner
The PLS acquisition started with the design of the scanning trajectory. This was carried out in such a way as to minimize occlusions, while optimizing acquisition time and ensuring a closed loop to minimize SLAM drift (Gollob et al. 2020). The operator walked along the designed route while holding the scanning device. The raw LiDAR data were processed with GeoSLAM Hub software and one single point cloud was obtained. To streamline the subsequent processing, point cloud density was reduced using an Octree decimation algorithm with ten subdivisions. This final point cloud was classified and normalized like the previous one; the noise was removed as well.

Estimation of dasometric variables
The procedure to estimate the dasometric variables from the point clouds obtained was designed to be fully automated. It is based on the identification of individual trees, with each subset being subsequently further processed to estimate the DBH and heights of those trees.

Recognition of individual trees in the point cloud
The recognition of individual trees in the point cloud was done by slicing the point cloud and performing a cluster analysis of the slices. This is outlined in a diagram in Fig. 4. The first step in this workflow consists of extracting a thin slice of the point cloud at breast height (DBH-1.3 m). The X, Y coordinates of the points contained in the slice were used to perform a cluster analysis with the aim of grouping into a cluster all of the points that belong to a certain tree stem. The cluster algorithm that was used is the density-based cluster algorithm (DBSCAN) (Ester et al. 1996). This algorithm is a density-based clustering algorithm which clusters points that are in areas of high point density, separating them from points lying in areas with low point density values. Consequently, points that are in low density regions can be considered outliers. This is an important step, because it avoids mistakenly clustering as tree stems points that correspond to different features. For example, scattered points that belong to the upper part of shrubs can be mistakenly included in a tree stem cluster, and this step serves to avoid this problem. The configuration of two parameters is needed: the maximum Euclidean distance permitted to consider two points Third, the horizontal projection of the slice is subjected to a cluster analysis. Fourth, each obtained cluster is adjusted to a circumference as neighbors (eps) and the minimum number of points that a cluster can contain (minPts). This process was performed for both the TLS and the PLS point clouds. Different slice thicknesses were tested for each of the clouds to find the most suitable size considering the noise, occlusions, and irregular spatial distribution of points. The TLS and the PLS point clouds required different eps and minPts to efficiently obtain the tree stem clusters.
Once the cluster analysis is performed, the group of points within a cluster is adjusted to a circumference line using Umbach and Jones's (2003) full least-squares method. This function minimizes the sum of the squares of the distances from the points to the circle (Umbach and Jones 2003). The X, Y central coordinates of the fitted circles were extracted and considered to be the central coordinates for each tree.

Diameter estimation
The DBH estimations were obtained by extracting the radius of each circle generated in the previous step (see "Recognition of individual trees in the point cloud"). These were compared with the measured diameters collected in field work. Linear regressions were performed and the Root-Mean-Error Square (RMSE) was calculated.
Diameters at higher heights are commonly used to obtain stem curves (Holopainen et al. 2011) and for volumetric estimation purposes. In this case, diameters at 4 m above the ground were measured in both point clouds to compare the performance of the two different scanners on this task. Point clustering, circle fitting, and circle radius extraction were performed analogously to the DBH measurement procedure. Since the resolution of the point cloud decreases with increased distance from the scanner to the object (Liang et al. 2014), the estimation of diameters at this height was done using the original point clouds, prior to decimation.

Height estimation
To estimate the heights of individual trees, the CHM (Canopy height model) was created from the normalized point cloud. This model consists of a raster containing information about the height of the vegetation in each pixel. The CHM was obtained using the "LidR" package from the R software. A 0.5 m resolution was selected, in an attempt to strike a balance between height generalization and excessive height detail. Subsequently, the height of each tree was automatically estimated by extracting the CHM value for its previously calculated central coordinates. This procedure was performed for both point clouds. The heights obtained were compared with the data obtained in the field. Linear regressions were performed and the root-mean-error square (RMSE) was calculated. Nevertheless, as problems often arise when using the Vertex to measure tall trees (Fernandes da Silva et al. 2012) or due to a lack of visibility of the tree apex (Larjavaara and Muller-Landau 2013), a new set of ground truth values was obtained. Tree heights were measured directly from the PLShh point cloud as this was the point cloud that had the least occlusion in the upper part of the canopy. An example is presented in Fig. 5. Table 2 shows the time invested in collecting data for each of the two scanning devices. The TLS point cloud acquisition took 37 min, while only 8 min were needed for the PLS.

Point cloud acquisition and post-processing
The registration error in co-registering the two TLS point clouds acquired was an RMSE of 14.6 mm. The integrated TLS point cloud had 11,589,772 points. For the PLS, the point cloud was directly downloaded from the scanning device. It had 41,902,931 points. In this case, the density reduction step mentioned in the methodology was applied to obtain a decimated point cloud; it contained 12,942,132 points. Figure 6 shows the plane and perspective views of the final normalized point clouds. This figure illustrates the differences between the two different devices' point clouds in terms of the completeness of the objects. This can be appreciated both in vertical and horizontal structure. Due to the occlusions typical of the TLS, there are few points that reach the tops of the trees in the upper part of the canopy. The PLS

Recognition of individual trees in the point cloud
The point cloud slices extracted at DBH from each point cloud are presented in Fig. 7. The final slice thicknesses selected for each point cloud were 2 cm in the case of the TLS and 1 cm in the case of the PLS.
The final minPts and eps selected differed between the two point clouds. For the TLS point cloud, the eps was 0.4 and the minPts was 10. In the case of the PLS, the eps was 0.5 and the minPts was also 10. As an example of the clustering step, the cluster result for the PLS point cloud is shown in Fig. 8.

Diameter estimation
The linear regressions comparing the diameters measured using the traditional field methods with the diameters measured using the TLS and the portable LiDAR are shown in Figs. 9 and 10, respectively. The R 2 obtained when comparing the TLS to the field data was R 2 = 0.5782. The R 2 obtained when comparing the PLS to the field data was R 2 = 0.8735. When using the TLS for DBH estimation, an RMSE of 0.176 was obtained. A higher RMSE was obtained when using the PLS for DBH estimation (RMSE = 0.071). An outlier was identified in the TLS data set; it was tree 12 (see Figs. 9, 10).
For the diameter estimation at 4 m, the thickness used in the extraction of the slice was 2 cm for the TLS point cloud and 1 cm for the PLS point cloud. Figure 11 shows a horizontal projection of the slices. In the case of the TLS slice, it is noticeable that some of the stems are not represented, while in the case of the PLS slices, some of the stems have a great deal of surrounding points that most likely correspond to noise.
The final parameters for the cluster analysis at 4 m were the following: for the TLS point cloud: minPts = 10 and eps = 0.4; for the PLS point cloud: minPts = 10, and eps = 0.5. Table 3 presents the final diameters obtained at a height of 4 m from the three point clouds (TLS, reduced PLS, and complete PLS). It should be noted that some of the estimated values are overestimations, since they make no physical sense: diameters should be smaller at greater heights. In fact, of the diameters at 4 m that were estimated from the reduced PLS point cloud, 58.3% can be deemed overestimates. However, when the complete PLS point cloud is used to estimate the diameters, the percentage of overestimates is only 23.08% .

Height estimation
The CHM raster layers obtained from each point cloud (TLS and PLS) are shown in Fig. 12. The cluster analyses are overlaid with the central coordinates detected for each tree. The CHM's statistics, the maximum and mean height of each CHM raster layer, are shown in Table 4. Figures 13 and 14 show the linear regression between the tree heights measured in the field using the traditional methods and the tree heights obtained from the TLS and PLS point clouds, respectively. The R 2 obtained from the comparison of the ground truth height with the TLS is lower than the R 2 obtained when comparing the ground truth height with the PLS (TLS height R 2 = 0.0534, PLS height R 2 = 0.1149). When using the TLS for height estimation, the calculated root-mean-square error (RMSE) was 19.593. The RMSE was 1.4 m lower when using the PLS for height estimation (RMSE = 18.189). These results are detailed in Table 5. The correlation between ground truth height values measured directly from the PLS point cloud, and height values measured automatically from each of the two different point clouds are shown in Fig. 15 (TLS) and Fig. 16 (PLS). In both cases, the R 2 values were greater than in the case of the the R 2 values when comparing the height measurement using automatic methods vs. traditional methods. Furthermore, the R 2 value remains greater for the PLS's height measurement, (R 2 0.867), than for the TLS's height measurement (R 2 0.505). An RMSE of 10.712 was obtained when using the TLS for height estimation, while the RMSE was lower (RMSE 3.415 m) when using the PLS (see Table 5).

Discussion
According to the results presented above, the efficiency in data acquisition of PLS systems as compared to TLS systems or traditional fieldwork seems clear, an observation that has also arisen in the other recent studies focused on this topic (Balenović et al. 2021;Bauwens et al. 2016;Chen et al. 2019;Gollob et al. 2021;Chiappini et al. 2022). This is a key point when it comes to ensuring the possibility of incorporating LiDAR data into the forest inventory reference data.
Another key feature of the results presented here is the possibility of the automated detection of the trees within a stand. The methodology presented here appears to be successful at this task. Previous studies also rely on DBSCAN clustering to identify individual trees (Hyyppä et al. 2020). The use of this methodology conveys an advantage when compared to other methodologies, since it avoids constraints such as the need for training data (Zeybek and Vatandaslar 2021) or a reliance on commercial software (Levick et al. 2021). Such constraints can hinder the viability, from an operational standpoint, of the implementation of the PLS in forest inventories. Additionally, it was observed that after the tree detection and individualization, it was possible to obtain DBH and height metrics from the PLS point clouds. In fact, higher accuracy metrics for both DBH and height were obtained with the PLS than with the TLS. Nevertheless, PLS systems do not always yield more accurate DBH estimations; this depends on many factors such as the accuracy and range of the device used, the number of scanning positionings used for the TLS and the trajectory followed when scanning with the PLS. An example of this is analyzed in the study conducted by Balenovic et al. (2021). In the case of the present study, it should be noted that large trees are being studied (roughly 40 m tall). Some studies have previously reported differences in the accuracy of the data obtained depending on the size of the trees. For example, Gollob et al. (2020) found that DBH was generally underestimated when the DBH was greater than 10 cm and that the magnitude of this underestimation became greater with increased DBH. Xie et al. (2022) stated that when trees are big, the shapes of their stems' cross-sections tend to be irregular, thus, increasing deviations between different measuring methods. Some other authors, such as Ryding et al. (2015), Bauwens et al. (2016), andHartley et al. (2022), also noted that irregularities (bark roughness and textures, or non-regular cross-section shapes) have an influence on the measurement of diameters. In fact, the greatest deviation coincides with the largest trees (for example, trees number 7 and 12). The trees in the case of this study also present quite noticeable bark roughness and irregular shapes (especially in the case of tree number  12). These irregularities can lead to errors when automatically fitting the scanned points to circumferences, and furthermore, this effect may be greater if there are trunk occlusions. In fact, when tree number 12 is removed from the analysis, the R 2 of the TLS measurements increases to 0.95, exceeding the value obtained for the PLS, 0.87. Although this circle-fitting approach can lead to errors, as shown here, it is an approach which is efficient and easy to implement from an operational point of view and it has been followed in different studies Hyyppä et al. 2020).
Greater height accuracies are to be expected with the PLS as there are fewer occlusions than with the static TLS systems (Chiappini et al. 2022;Hyyppä et al. 2020). This might not only be crucial for obtaining higher height accuracy metrics but may also be of interest in terms of going beyond the traditional measurements. For example, PLS systems can be used effectively to obtain diameters at different heights, as shown in this study, and later extract stem curves and estimate volumes, as shown in Hyyppä et al. (2020).
In general, the accuracy metrics obtained in this study are lower than in the previously mentioned studies. For DBH estimations, Giannetti et al. (2018) reported an RMSE of 0.0113 m with a TLS point cloud, and an RMSE of 0.0128 m with a PLS point cloud. The study of Oveland et al. (2018) reported DBH estimations with an RMSE of 0.062 m using a TLS point cloud and an RMSE of 0.031 using a PLS point cloud. Gollob et al. (2020) According to the results of this study, DBH measurements can be efficiently obtained from PLS point clouds in Eucalyptus stands. In this line, Marselis et al. (2016) calculated DBH and obtained similar RMSE and R 2 values (3.79 cm and 0.72). However, these results contrast the conclusions of some other previous studies (Levick et al. 2021;Camarretta et al. 2021). This discrepancy may be due to the varying complexity of the stands analyzed.
The previously mentioned studies do not test the PLS's ability to calculate individual tree height, while this study highlights the superiority of the PLS in comparison to the TLS in this regard. Therefore, the results presented here show that PLS systems can be used efficiently for inventory purposes in Eucalyptus stands similar to the stand analyzed here. Considering the importance of the Eucalyptus genus for the forestry sector worldwide (FAO 2005), an important next step for the future of this research will be to test the capacity of the proposed methods in more complex Eucalyptus stands.

Conclusion
In conclusion, this study demonstrates the utility and suitability of the PLS with the SLAM algorithm for estimating forestry parameters of individual trees in Eucalyptus globulus stands. The potential of these systems is undeniable, and an effort is required on the part of the scientific community to develop different methodologies for their application. Aiding in the extraction of dasometric variables can facilitate the use of these types of systems that have the potential to become a powerful tool for monitoring the evolution of forest masses of interest and even to help optimize production masses.
The most notable contribution of this study is that it developed a methodology for automatically processing LiDAR data acquired from a PLS system to obtain the dasometric variables of a Eucalyptus globulus stand. This is noteworthy due to the importance of this species for economic purposes worldwide. Furthermore, the entire methodology is based on the use of public domain software, which is crucial in terms of its replicability in different stand conditions.

Data availability
All the authors declare that the data presented in the paper belong to the authors.

Conflict of interest
The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.