GEOREFERENCING UAS DERIVATIVES THROUGH POINT CLOUD REGISTRATION WITH ARCHIVED LIDAR DATASETS

Georeferencing gathered images is a common step before performing spatial analysis and other processes on acquired datasets using unmanned aerial systems (UAS). Methods of applying spatial information to aerial images or their derivatives is through onboard GPS (Global Positioning Systems) geotagging, or through tying of models through GCPs (Ground Control Points) acquired in the field. Currently, UAS (Unmanned Aerial System) derivatives are limited to meter-levels of accuracy when their generation is unaided with points of known position on the ground. The use of ground control points established using survey-grade GPS or GNSS receivers can greatly reduce model errors to centimeter levels. However, this comes with additional costs not only with instrument acquisition and survey operations, but also in actual time spent in the field. This study uses a workflow for cloud-based post-processing of UAS data in combination with already existing LiDAR data. The georeferencing of the UAV point cloud is executed using the Iterative Closest Point algorithm (ICP). It is applied through the open-source CloudCompare software (Girardeau-Montaut, 2006) on a ‘skeleton point cloud’. This skeleton point cloud consists of manually extracted features consistent on both LiDAR and UAV data. For this cloud, roads and buildings with minimal deviations given their differing dates of acquisition are considered consistent. Transformation parameters are computed for the skeleton cloud which could then be applied to the whole UAS dataset. In addition, a separate cloud consisting of non-vegetation features automatically derived using CANUPO classification algorithm (Brodu and Lague, 2012) was used to generate a separate set of parameters. Ground survey is done to validate the transformed cloud. An RMSE value of around 16 centimeters was found when comparing validation data to the models georeferenced using the CANUPO cloud and the manual skeleton cloud. Cloud-to-cloud distance computations of CANUPO and manual skeleton clouds were obtained with values for both equal to around 0.67 meters at 1.73 standard deviation.


INTRODUCTION
The significant advancement of photogrammetric processing software has transformed UAVs (Unmanned Aerial Vehicles) from being recreational devices into important tools for the geospatial field.UAVs or pilotless remotely controlled aircraft capture images from an aerial perspective.These images can further be processed to generate Digital Elevation Models (DEMs), and Orthomosaics for mapping.Without using GCPs, GPS tracking devices that can be carried by these small instruments are limited to meter-levels of accuracy such as the popular Sensefly eBee (Sensefly, 2015).Adding GCPs (Ground Control Points) in processing can greatly reduce model errors to centimeter levels but these can only be acquired by doing ground surveys in the field.For survey grade mapping of large areas, such as whole cadasters, UAS manufacturers have developed a way of using RTK with the UAV simultaneously.As the extra equipment brings accuracies to centimeter level, the cost of procuring the system as a whole also increases.Moreover, using any of these extra steps in the field would require more manpower and more time.Eliminating on-site methods, this research looks at post-processing steps instead.
One alternative method is by using already existing data of relatively higher levels of accuracy and use them as reference for the UAV derivatives.Such reference data must also have high availability to be of actual use or at least have overlaps with the study areas to be mapped.For the case of the Philippines, a high accuracy collection of data is provided by its nationwide Aerial LiDAR Surveys (ALS).From as early as 2011, parts of the Philippines have been mapped using LiDAR (Light Detection and Ranging) with derivatives being used for flood modeling and hazard mapping (DREAM Program, 2012).Aside from disasters and hazards, programs have also been conceptualized with aims for using those same point clouds in providing a detailed inventory of the country's natural resources (Blanco et al, 2015).In 2016, an online web portal was launched that serves as a data distribution center of these various LiDAR derivatives such as DEMs, maps, orthophotos, and classified point clouds.These LiDAR data have accuracies of 0.5 meters on the horizontal plane and 0.2 meters along the vertical axis.
Given this collection of reference data, automated georeferencing becomes a viable approach for UAV data.A similar concept applied to TLS or Terrestrial Laser Scanning can be used for this study.TLS is a grounded version of ALS opting instead for a terrestrial perspective in scanning objects and environments.In place of an aerial platform, an instrument set-up is placed to perform the scan and if need be, multiple set-ups to cover the whole scene of interest.A simple step called Point Cloud Registration could get a whole environment in one coordinate system by fusing datasets from multiple views (Rajendra et al, 2014).Cloud registration is performed along overlapping areas of different scans and as those areas of overlap become consistent with each other, so do the rest of the dataset.The parameters used for this process include translation and rotation.This ensures that the cloud transforms only with location and distance.

METHODOLOGY
Imagery was acquired for this study using a Sensefly eBee drone flown over the College of Science Complex of the University of the Philippines -Diliman.It is a fixed-wing type UAV with a 96 cm wingspan weighing approximately 700 grams.It performs data acquisition by means of an 18 MP RGB camera loaded inside the drone with a circular opening at its belly for the lens.The UAS included a flight planning tool that automatically designs a mission for the flight using the desired parameters such as image overlap, ground sampling distance, and survey area.The option of doing a crisscrossed flight was implemented meaning that the whole area was surveyed using two sets of flightlines with the second set perpendicular in direction to the first one.This ensures a higher overlap and redundancy for the data and minimizes any effect of deformation in the generated cloud.With a less deformed cloud, a rigid transformation can be implemented on the dataset.A crisscrossing flight also minimizes shadows in the created cloud by increasing viewing angles with increased camera stations.The same software is also used in geotagging the collected images using data from an onboard GPS (Global Positioning System) and IMU (Inertial Measurement Unit).
The 3D modeling itself was done in Pix4D by processing the geotagged images.An initial sparse point cloud is first created after matching recurring key points in multiple images that were detected using computer vision algorithms.This process also refines the spatial coordinates from the image geotags.Next, a denser point cloud is produced by using camera coordinates refined from the previous step as well as the camera calibration parameters automatically derived.The results for point cloud densification were saved and exported in LAZ format.Normal Pix4D pipeline finishes with mesh and orthomosaic generation, but for this study the raster files produced by Pix4D were used for feature detection to be used in registration.The final UAS derivatives will be generated using LAStools (Rapidlasso, GmbH).
Due to the immense size of the obtained point cloud, it is preferable to split it first into smaller tiles and then 'sort' the data.Pix4D sequences each point by the number of matches from the obtained images, the tool lassort will sequence the points instead by location.This step should reduce the processing times for the succeeding tools in the pipeline such as DEM generation in blast2dem which uses Delaunay triangulation before generating a raster.A separate lighter version of the sorted cloud was used in obtaining transformation parameters for faster computations.A thinning algorithm was used to thin the sorted point cloud by picking central points on a grid with dimensions defined by the user as a step size parameter.Generation of point clouds by photogrammetry is also prone to outliers.This thinning step lessens high and low noise in order to prevent them from being included in registration computations, One problem of performing comparative studies between 3D data from different acquisition dates is the natural tendency of features to change.These differences could erroneously contribute to the calculated errors or worse alter the registration results if not noted correctly during processing.Obvious elevation differences due to vegetation and structural developments can easily be detected using DSMs or archived imagery such as Google Earth.The orthomosaic produced by Pix4d showed well that the majority of change can be attributed to vegetation growth and at certain locations, the lack thereof.
Upon detection, features of significant change was masked from processing to maintain the integrity of both georeferencing and validation steps.Consistent features between the two 3D datasets was then selected such as roads and buildings.A drawn polyline shapefile following non-occluded road centerlines can be buffered to estimate whole roads.Additionally, consistent roofs of buildings can be extracted by following their edges from the orthomosaic.These features were extracted from the thinned cloud using lasclip as the 'skeleton' point cloud.Using the thinned cloud eliminated noise along building facades ensured the planarity of smooth surfaces for better registration with their LiDAR counterparts.
For larger datasets, manual extraction of possible consistent features is a tedious task.A classification algorithm used for LiDAR point clouds was tested for this process but yielded less than favorable results.For this study, a point cloud classification algorithm implemented in CANUPO software suite (Brodu and Lague, 2012) was used to separate vegetation from nonvegetation classes as the first class contributes largest to the changes between the two datasets.The algorithm uses a multiscale approach wherein classes are distinguished through dimensionality in varying scales.For example, as discussed by its authors, CANUPO observes vegetation in a centimeter scale as objects comprised of 1-dimensional (sticks) and 2-dimensional objects (leaves), while as a 3-dimensional complex shrub object when viewed in a larger scale (decimeters and above).A smooth ground area however, consistently appears as a 2-dimensional object in both scales.Figure 3    Given the 'skeleton' point cloud and a reference point cloud, registration using the Iterative Closest Point algorithm could then be executed.ICP algorithm is an iterative process that can be divided into three steps.Considering two point clouds A and B, the algorithm first performs matching by identifying for each point in cloud A its nearest counterpart in cloud B. The distance between all point pairs is considered as the error metric.The goal of the second step is to minimize that error metric.Transformation would then constitute the final step as the data points are transformed depending on the spatial deviations of both clouds.These three are done repeatedly until the maximum number of iterations has been reached or the algorithm itself has converged (Kjer and Wilm, 2010).
For ICP to converge to good results, the skeleton point cloud was roughly translated along the horizontal plane.The transformation parameters on the skeleton cloud was extracted from CloudCompare's console and was easily applied on the entirety of the dataset.The resulting georeferenced point cloud, as well as the georeferenced skeleton cloud was saved and compared with the LiDAR dataset.This was done for the two skeleton point clouds: manual identification and CANUPO classification.
Cloud-to-cloud distance computations were executed as well for the two registered clouds versus the reference LiDAR cloud.The values and their standard deviations were recorded.In addition, a ground validation survey was conducted using a Total Station.An observation point was established and the distances relative to this point were measured along the XY-plane and Z-axis.A total of six points were observed in the field and identified in the registered orthomosaic.A comparison of model distances and ground distances yielded error values used for accuracy assessment of the georeferenced models in this study.
With a georeferenced point cloud, LAStools was used in generating a new set of DSMs and a new set of orthomosaics using the RGB statistics of points in the registered dataset.

RESULTS AND DISCUSSION
A dense point cloud consisting of around 106 million points was generated using Pix4d for an area of 16 hectares.This translates to a density of almost 72 points per cubic meter with a ground sampling distance of around 4 cm per pixel.However, the processed point cloud has to be clipped as the decrease in image overlap along its edges affects the integrity of modeling for these areas.What remained is a point cloud spanning 16 hectares on the ground with dimensions of 400 meters along the x and y axes and composed of around 106 million points.The camera alignment step of the modeling process reports RMS errors of 0.38, 0.41, and 0.43 m errors for x, y, and z respectively.This means that locally, the model has an accuracy at the sub-meter level.The total processing time from camera alignment to point cloud densification was clocked at 11 hours.Figure 5 shows the output of image processing: the generated DSM and RGB orthomosaic.
Road centerlines nearly 1 kilometer in length were digitized and buffered three meters on either side to generate road polygons.Consistent buildings were also added to this polygon shapefile which makes the total area of the skeleton point cloud equal to 2.5 hectares.This cloud was translated giving the translational matrix followed by the ICP matrix.Figure 6 delineates the original position of the skeleton point cloud over the georeferenced mosaic after performing registration.Overlaying the two shows a significant change in position after performing the transformation.For the second cloud, the non-registered orthomosaic generated by Pix4D was segmented using eCognition instead of manually delineating image features for training data collection.Selection of various objects that represent the two classes, vegetation and non-vegetation, was done and fed to the CANUPO algorithm for classification of the whole cloud.For this branch of processing, the non-vegetated classes are used in registration instead of a manually selected skeleton cloud.The workflow proceeds with registration of the whole UAS dataset and the generation of derivatives.
Figure 7 shows a 3-dimensional view of the UAS point cloud overlaid on the LiDAR point cloud after one of the registrations.
The UAS data has blended near seamlessly with the position of the LiDAR data while also providing an 'updated look' on the LiDAR cloud as features that weren't present in the LiDAR cloud before but are present now can be appended using the registered UAS data.
Raster subtraction of the generated DSM was done to illustrate areas of high and low elevation differences.As expected, vegetation growth constitutes negative values while the higher ones are from newly constructed buildings.Cyan colored areas in Figure 8 show pixels with elevation differences between positive and negative 20 cm.The red arrow point to buildings under construction that were finished in between the acquisition dates of LiDAR and UAS data.
The results of the field validation with the use of a Total Station is summarized in Tables 1 and 2 was observed that registered derivatives using the Manually Identified and CANUPO classified skeleton cloud had nearly similar RMSE values in the XYZ-axes.Both workflows had relative distance differences of around 16 centimeters from the field data.When compared to UAS derivatives before registration, this is a significant improvement considering that no GCPs were used in aiding the photogrammetric processing of the obtained images.
Using cloud-to-cloud distance computations for further comparison of the results of the two pipelines yielded a similar observation.For the Manual Identification cloud, a distance of 1.339 meters was computed with a standard deviation of 2.154, while the CANUPO Classification cloud computed a distance of 1.342 meters at 2.153 standard deviation.The meter level distances and standard deviations can be attributed to the fact that the entirety of both results were used in cloud-to-cloud distance computation.This included vegetation that experienced meterlevels of elevation changes as well as new buildings that changed between the acquisition dates of the two datasets.Performing the distance algorithm only on non-vegetation classes, the CANUPO cloud had a distance from the reference cloud of 0.68 meters at a smaller standard deviation of 1.73.These values are similar to the manually derived cloud with cloud-to-cloud distance of 0.67 meters with a standard deviation of 1.73 meters as well.The remaining source of deviations could be caused by elevation differences by new buildings which were still included in the non-vegetation class of the CANUPO cloud.The movement of the cloud before and after registration can further be illustrated by analyzing the cloud distances of both Manual and CANUPO cloud subsets.Before georeferencing, the CANUPO cloud had a cloud distance of 1.20 meters while the Manual cloud had a distance of 1.40 meters.This is nearly a 50% reduction in distance from the LiDAR data.
The similarity in accuracy assessment values for both pipelines, manual identification and CANUPO classification, suggests that for this dataset, CANUPO classification produces similar results while reducing the amount of time for pre-processing of the cloud.Instead of relying solely on the judgement of the processor for the entirety of the cloud, subsets could be used to permit automation.This could be helpful when performing a similar georeferencing method for larger datasets.

CONCLUSIONS
In this study, a post-fieldwork methodology was designed to georeference derivatives generated from UAS images.With this workflow, a survey party could proceed with the UAV survey itself without the need for increasing the party's size, period in the field, or procuring specialized instruments in the hopes of getting data of higher spatial accuracy.
One of the limitations of this method of model production would be on the changes inherent in the area of study.High amounts of changes between the two datasets would limit the UAS data from having a reliable basis from which to derive the appropriate transformation parameters.
It is worth noting that the reference data used, even though LiDAR data, was still derived from a raster dataset.With LiDAR's multiple number of returns and and non-grid spacing, it is recommended to do further research on whether its use as reference could produce better results.Nevertheless, readily available LiDAR data in the Philippines are often archived and distributed as generated raster files for ease of distribution.For users especially those not in the research field with no access to these huge LAS files restricted point clouds, this research has shown that the use of archived LiDAR DSMs would suffice for georeferencing of UAV derived datasets.
Figure 2. The skeleton point cloud consisting of roads and roof of consistent buildings.

Figure 3 .
Figure 3. Brodu and Lague's dimensionality diagrams for four classes in a mountain river environment.

Figure 4 .
Figure 4. Equivalent point cloud of the LiDAR data used as a reference.

Figure 5 .
Figure 5.The georeferenced DSM and Orthomosaic of the study area.

Figure 6 .
Figure 6.Clipping polygons overlaid with the orthomosaic after registration.

Figure 7 .
Figure 7.The registered UAV point cloud overlaid with the reference LiDAR point cloud.

Figure 8 .
Figure 8. Raster difference of the the registered DSM and LiDAR DSM (units are in meters).
. From the two validation sets, it

Table 1 .
Manually Identified model to ground differences of measured distances from reference point

Table 2 .
CANUPO classified model to ground differences of measured distances from reference point