SEMI-AUTOMATIC CO-REGISTRATION OF PHOTOGRAMMETRIC AND LIDAR DATA USING BUILDINGS

In this work, the co-registration steps between LiDAR and photogrammetric DSM 3Ddata are analyzed and a solution based on automated plane matching is proposed and implemented. For a robust 3D geometric transformation both planes and points are used. Initially planes are chosen as the co-registration primitives. To confine the search space for the plane matching a sequential automatic building matching is performed first. For matching buildings from the LiDAR and the photogrammetric data, a similarity objective function is formed based on the roof height difference (RHD), the 3D histogram of the building attributes, and the building boundary area of a building. A region growing algorithm based on a Triangulated Irregular Network (TIN) is implemented to extract planes from both datasets. Next, an automatic successive process for identifying and matching corresponding planes from the two datasets has been developed and implemented. It is based on the building boundary region and determines plane pairs through a robust matching process thus eliminating outlier pairs. The selected correct plane pairs are the input data for the geometric transformation process. The 3D conformal transformation method in conjunction with the attitude quaternion is applied to obtain the transformation parameters using the normal vectors of the corresponding plane pairs. Following the mapping of one dataset onto the coordinate system of the other, the Iterative Closest Point (ICP) algorithm is then applied, using the corresponding building point clouds to further refine the transformation solution. The results indicate that the combination of planes and points improve the co-registration outcomes.


INTRODUCTION
Three dimensional building data are used extensively in numerous city applications such as modeling of urban environments, city planning, emergency response, environmental assessment and determining spatio-temporal urban changes.Building databases are generated by various data sources but mainly aerial and satellite images and airborne laser scanning data are the most common ones.The complimentary characteristics of LiDAR and photogrammetric data for data conflation or change detection required the proper coregistrations these two data sets.Therefore horizontal and vertical alignment of various 3D datasets is necessary by determining the rigid transformation that co-registers different 3D datasets in a common coordinate reference system using primitives, such as points, lines and planes.Considering the urban environment and the 3D building modeling applications we observe that linear features and planar surfaces are important invariant geometric and semantic elements.Planes and lines are typical primitives that can be found in man-made environments such as cities due to the nature of the city landscape.For example, buildings are formed by planar surfaces.Habib et al. (2004) used straight lines as the registration primitives, while planes for surface matching have been also used (Habib and Schenk, 1999;Schenk and Chatho, 2002;Sampath andShan, 2006, Brenner et al., 2008).The quality of the integration outcomes unquestionably depends on the registration process towards respective data (Postolov et al., 1999).Current approaches rely on the manual identification of corresponding spatial features (points, lines planes) between the two data sets to be used for determining the transformation parameters between the two coordinate systems.However, the volume of data makes the workload of manual registration quite tedious.Hence, automatic 3D datasets co-registration processes are on demand to efficiently align different datasets into a common reference system.
Given two arbitrary datasets, overlapping regions between the datasets must be determined and corresponding primitives must be defined within the common regions (Planitz et al., 2005).A pair-wise correspondence is required to ensure primitives matching between feature representations.The process of primitive matching allows an algorithm to identify potential correspondence between primitives, and more importantly discard grossly erroneous matches.This ensures that matches only occur between similar features.Therefore, it is essential to determine how to extract features/primitives from a dataset.Considerations such as which feature extraction algorithm to be used are necessary.Features must be represented in such a way that they are comparable to other feature representations.
In this paper emphasis is given in introducing automation in the co-registration process of LiDAR and photogrammetric DSM 3D points, particularly by automatically identify and generate conjugate geometric primitives in the different data sets to serve as anchor features for mapping one dataset onto another; and automatically localize and determine the correspondence of common primitives.This is followed by the actual mapping of one dataset onto the spatial reference system of the other through a 3D geometric transformation and by assessing the accuracy of the proposed method.For a robust 3D geometric transformation both planes and points are used.Initially planes are chosen as the co-registration primitives.A region growing algorithm based on a Triangulated Irregular Network (TIN) is implemented to extract planes from both datasets.Point clouds are subsequently used as another registration primitive to complement the plane-based registration.Next, an automatic iterative process for identifying and matching corresponding planes from the two datasets has been developed and implemented.The extracted planes are associated as plane pairs, initially by a building matching process which is then followed by the plane matching algorithm.After the correspondence between the buildings has been established, a plane matching process is implemented to create plane pairs through a robust matching process, which thus eliminates outlier pairs.The selected correct plane pairs are the input data for the geometric transformation process.Three different geometric registration algorithms are used to obtain accurate transformation parameters between the two datasets.Following the mapping of one dataset onto the coordinate system of the other, the corresponding building point clouds are used to further refine the transformation solution.

IDENTIFICATION OF PRIMITIVES
In our case we are interested in integrating photogrammetric and LiDAR data.The space of comparison was based on the minimum processing of the data sets.Thus, the framework for the comparison is based on two 3D point sets that cover the same area.However, there is no one-to-one point correspondence between the photogrammetric and LiDAR point-based 3D surfaces.In addition the 3D point data covering the same area usually are geo-referenced either to local coordinate systems or to ones defined by their acquired sensors.Thus, even in the case of having corresponding features their 3D positions will not coincide.The critical task is to automatically identify common features between the two 3D point data sets that will allow the determination of the parameters of the 3D transformation.
Traditionally comparison of digital surface models (DSM) data is based on establishing the point correspondence by establishing a common XY grid and derived the Z values of each DSM by interpolation.This allows for the estimation of the displacements between the two data sets but obviously is error prone to the interpolation, particularly in urban areas.Another approach used for deriving the alignment between two point clouds is the Iterative Closest Point (ICP) method (Besl and McKay, 1992).It is an iterative procedure requiring good initial transformation values and it is based on the minimization of the distance between a point belonging to one data set and its closest points in the second data set.In this work we propose a two prone approach.Initially use the normal vectors of plane primitives in corresponding buildings to determine the transformation parameters.Then use the ICP algorithm applied to the initial 3D point clouds of the corresponding building using these derived transformation parameters as the initial transformation parameters.There are numerous planar patches in urban environments, especially roofs of buildings, which make the selection process of common planes incredibly laborious.It will save significant time and effort if the search pool for plane correspondence can be narrowed down into building level.Therefore, in our case the question is to how to first automatically determine the corresponding buildings in the two data sets and second how to automatically identify the corresponding plane pairs within the area covered by each of the corresponding buildings.

BUILDING EXTRACTION AND MATCHING
The LiDAR and DSM points of buildings were extracted using the TerraScan software (Terrasolid, 2011).Initially low elevation and ground points were extracted.The building candidate points were extracted by height distance from ground points followed by manual editing to remove any non-building points.Buildings have various characteristics and identifying the same building between the photogrammetric DSM and LiDAR point cloud often requires a lot of information, time and effort.The buildings were much through a stepwise process, first by grouping them in single and multi-layered buildings, then by applying three matching tests.

Grouping of single buildings
First a building group algorithm is implemented based on a common building characteristic.All buildings are divided into two main categories: Single-layer roof buildings and multi-layer roof buildings using a grouping algorithm.The algorithm is based on a common characteristic of the building called Roof Height Difference (RHD).The Roof Height Difference is defined as the difference between the highest roof and the lowest roof patch of the same building.Since single-layer roof building has only one roof, the RHD of any single-layer roof buildings is 0 while that of a multi-layer roof building is a positive number.Therefore, the RHD parameter is introduced as a shift and rotation invariance factor to determine whether a building is a single-roof building or not.

Multi-layer buildings -first matching test
All multi-layer roof buildings are further divided into several smaller groups according to their RHD.To avoid a building pair being divided into two neighbouring groups, a buffer zone is defined by the standard deviation value  assuming a normal distribution of RHD for all buildings to enlarge the space of grouping.Following the building grouping, the buildings that have similar attribute are further grouped together.
Additionally to the RHD, we introduce a novel metric named 3D histogram of individual buildings as another element to better represent the complexity of a building and generate links for corresponding buildings in two datasets.The 3D building histogram can represent a random multi-layer building based on certain attributes.Let the x-axis represent the different directions of normal vector of each roof patch in a building by calculating the zenith angles of the normal vectors; the y-axis records elevation gaps between each roof patch; and the z-axis is the sum of areas of the roof patches which are located in the specific space defined by direction (x) and elevation (y).The building 3D histogram then contains both direct and inherent information about each building.An example can be seen in Figure 1 where (1a) displays the aerial image of a building; (1b) is the 3D histogram of the building from the LiDAR data; and (1c) is the 3D histogram of the same building from the DSM data.The higher class in  -The block distance between the two histograms, SD.SD 12 is used to calculate the distance between two columns x 1 and x 2 in the 3D building histogram: where, k is the number of dimensions of a histogram (k = 3 in this example).
-The z-value of the largest class in the 3D histogram, ZLB: The z-value of the largest class in 3D histogram (ZLB) is used as a criterion to identify if two random buildings are alike.The ZLB represents the area of the largest layer in one building.Small layers may vary due to noises, system errors and/or segmentation precision, but large and flat areas are always more resistant to those kinds of disturbance.
Besides these two indices of 3D building histogram, two other indices are also combined into the building matching equation.
-The area of the individual building, Area; -Roof Height Difference of an individual building, RHD.
Thus, the objective function SM of the similarity measure for building matching using the 3D building histograms is defined as: (2) where The similarity function is applied between one building B1 from dataset 1 and with all the buildings in the same group from dataset 2. If building B2 from dataset 2 is the corresponding building of building B1, then in the ideal case, the similarity measure between them will be 0 (perfect match).In reality however, when a building B2 found that minimizes the objective similarity function with building B1 this will be the corresponding building of building B1.Therefore: where n is the number of buildings in the same group from system 2.This pair (B1, B2) of buildings will be used for finding corresponding plane pairs within these matched buildings.

Multi-layer buildings -second matching test
We observe that the two histograms that represent corresponding buildings in two different datasets they are different in shape and can vary considerably from each other, due to the different point density of two different datasets and the goodness of segmentation.This is the major reason why some buildings that should be corresponding failed during the first building matching test.Therefore, a second building matching process is applied to reduce these failed building cases from elimination during the plane matching phase.The failed building cases are put together as a new group and a new building matching scheme is applied on this new group.Since the 3D histogram cannot correctly represent these buildings a new similarity matching objective function SM2, is defined as following: After the second test, most buildings should pass with their corresponding buildings.If some buildings still failed the matching test after the second test, then these buildings participate in a third test this time together with the all the single-layer buildings grouped during the initial grouping of all buildings.

Single layer and unclassified multi-layer buildingsthird matching test
From previous tests, we have observed that the coverage area of building is relatively more stable than the Roof Height Difference (RHD).What is more, single-layer buildings were not included in the first two matching tests because these buildings lack the RHD information.Therefore, the remaining multi-layer buildings from the first two matching tests and all single-layer buildings are grouped together and building area will be the only criterion used to calculate the similarity objective function SM3, where:

PLANE EXTRACTION AND MATCHING
First, a planar segmentation is applied to both LiDAR point clouds and photogrammetric DSM to generate plane patches for the buildings based on TIN generation.A region growing algorithm based on the similarity of the normal vectors of the TIN triangles was used for the extraction of planes.
Following the determination of building correspondence, plane matching is addressed at the building level within the coverage of the building, thus greatly reducing the spatial search window.Since the potential plane primitives should be easy to represent, a rectangularity test is applied initially before finding the corresponding plane pairs in order to select as much as possible regular shapes.Three parameters were used to characterize a plane.Its area and two shape measures, the rectangularity and the compactness.The two shape parameters describe how 'rectangular' and 'square' a plane patch is.They are defined as: If a plane has area, compactness and rectangularity values that fall into certain threshold zone, then it will be labelled as one of the candidate planes for the plane matching process.The values of the thresholds are determined experimentally based on the footprint size of the buildings and of known polygon shapes.
Then all the rectangular planes in one building pair will be matched according to their boundary box area.However, this plane matching based only on the plane area can introduce some matching errors because some planes in the same building may have similar area.In this case, a plane in one building in system 1 can find more than one corresponding planes in the corresponding building system 2. Therefore spatial distribution of the planes within the building was ap[plied to eliminate potential plane matching errors.This was done by comparing distance between the plane centroids in one dataset to that in the other dataset.

CO-REGISTRATION TRANSFORMATIONS
The final steps to the co-registration process are the actual mapping of one dataset onto the spatial reference system of the other through a geometric transformation function, and the accuracy assessment of the process.The proposed methodology for establishing the mapping parameters from one system to the other is based on a rigid 3D surface registration based on the common primitives between the two data sets.The registration process is performed in two steps to obtain accurate transformation parameters between the LiDAR and the DSM datasets.First using the corresponding plane pairs, a 3D conformal transformation based on their centroids and normal vectors at the centroids is applied to obtain the transformation parameters.Second using these estimated parameters the 3D point datasets is mapped into the other coordinate system is performed followed by the Iterative Closest Point (ICP) algorithm.In the latter the corresponding building point clouds are used to further refine the transformation solution.
The 3D conformal transformation parameters are estimated in two steps.The rotation matrix is estimated by least squares adjustment from the k pairs of conjugated normal vectors n as: Non-accurately derived planes and thus normal vectors lead to quite different angular parameters.In our case, and generally in urban environments, most plane primitives are horizontal planes and their normal unit vectors are very close to (0, 0, 1) values, that is they are almost parallel.Hence the values of the three rotational angles are numerically very sensitive to the quality of planes.Another way to represent attitude relationship between two coordinate systems is the rotation quaternion (or attitude quaternion).A quaternion q as either a 4x1 vector q= mathematical notation for representing orientations and rotations of objects in three dimensions.By using the characteristics of quaternion, we know that the relationship between rotation matrix R and q is as follows (Zhao, 2009): The R is orthogonal matrix, and contains only three unknowns according to the characteristics of unit quaternions.If (X) (Y) (Z) be the estimates of normal vector components of X Y Z, in system 2 then the normal equation of least squares of rotation quaternion is: The average translation of the k pairs of centroids is utilized for the estimation of the translation parameters.Finally, the 3D conformal transformation of data set A to the coordinate system of data set B is computes as: The transformation results were assessed in the object space as the differences of normal vectors

APPLICATION
The Keele Campus of York University was used as the study area.It has 66 building.The DSM was generated from UltraCam digital aerial images taken in 2005 from a flying height of about 1660m.It has a resolution of 0.78m and therefore a point density of 5.8 points/m 2 .Its X and Y coordinates are in the NAD83 UTM coordinate system and the Z elevations are orthometric heights.The lidar data were generated in 2009 from a flying height of about 2000m and have point density of about 3 points/m 2 .Its X and Y are in the NAD83 UTM coordinate system and its Z are ellipsoidal heights.The point density of photogrammetric DSM is approximately 2 times of that of LiDAR data clouds.
Among the 59 buildings that passed the building matching tests there are 53 buildings that contain at least one rectangular plane based on the rectangularity and compactness tests.The number of the rectangular planes from these buildings is 138 in LiDAR dataset and 212 in DSM system, respectively.That is to say, some buildings may have more than one plane that passed the rectangularity test while some others may not have any rectangular plane.Thus, 138 rectangular planes following the rectangularity test are the input of our plane matching process.After comparing the areas of two planes in the two different systems and eliminating outlier pairs using the linear fitting of the plane distances, 25 pairs of planes that satisfied the conditions set were matched and used as the co-registration corresponding primitives.These 25 pairs are from 21 buildings which contain 46 rectangular planes as candidate planes for matching.The low rate of generating these 25 plane pairs from 138 rectangular planes is attributed to: 1) the quality of planes extracted from the two systems that have different point densities; and 2) the strict constrains (thresholds of area and plane distances) to make sure the correspondence between planes are rightly connected.
From the 25 corresponding pairs, 15 pairs of the normal vectors of planes were randomly chosen as the input to the 3D conformal transformation.The other 10 pairs were used as "check normal vectors" for evaluation purpose.The differences at these check planes are given in Tables 1 and 2.  After the transformation parameters were calculated, they were applied to the LiDAR data points to transform them to the coordinate system of the DSM.Then the overlap area between the convex hulls of the transformed LiDAR buildings and their corresponding DSM buildings was calculated.Their estimated mean value is 98.60%.
The boundary of a building was selected as an example to display the transformation performance.In Figure 2, the blue line is the building boundary of the photogrammetric DSM and the red line is the building boundary of the LiDAR data.After applying the 3D-conformal transformation, it can be seen that there is still a slight  rotational bias between the transformed LiDAR building boundary and the corresponding photogrammetric building boundary.Since now the point clouds for each building for both data sets are closely located the ICP algorithm is applied to get the final adjusted transformation parameters.The ''iterative'' aspect of the ICP comes from the fact that the correspondences are reconsidered as the solution comes closer to the local minimum of the error.As any gradient descent method, the ICP is applicable when there is prior a relatively good starting point.The previously derived conformal transformation parameters were used as initial values for the ICP algorithm.After applying the ICP algorithm to the transformed LiDAR data shown now with blue line in Figure 2, the small bias has been almost eliminated.The overlap area between the convex hulls of the transformed LiDAR buildings and its corresponding DSM buildings was calculated.Their mean value was estimated to be 98.78% which is slightly higher than the overlap ratio (98.60%) that calculated after the 3D conformal transformation is performed.
Figure 2. Results from 3D conformal and ICP transformations

CONCLUDING REMARKS
In this work an automated method for the integration of photogrammetric and LiDAR data in urban 3D environments is being proposed.The alignment of heterogeneous data for their integration in urban environments serves numerous applications from database updating and mapping, to 3D modelling and city planning, to sensor calibration.Both planar and point data have been used to achieve a stable co-registration of LiDAR and photogrammetric DSM point data.Initially planes are selected as the main registration primitives because they are easy to be represented mathematically, they are stable spatial features, and building roofs are major features in urban environment.The approach is based on automatically deriving the conjugate planes between LiDAR and photogrammetric datasets.A building matching algorithm is first introduced to narrow down the spatial search window for the planes to be matched and speed up the process since there is no direct spatial relationship between the extracted planes in each dataset and the number of planes is large.Having defined the search space, this is followed by plane matching within the region covered by each building pair.After obtaining the plane pairs as registration primitives, two transformation methods, 3D conformal transformation along with quaternion rotation and the Iterative Closest Point (ICP) transformation were applied in order to obtain accurate and correct transformation parameters between LiDAR and photogrammetric data coordinate systems.The 3D conformal transformation applied first for the determination of the 6 transformation parameters.The quaternion rotation is also applied as we observed that the Euler rotations solution is sensitive to the directions and accuracy of the plane normal vectors.Then the ICP transformation was applied to already transformed building points to further refine the transformation parameters.Thus, the normal vectors and the centroids of the planes is the input for 3D conformal transformation, while the input for ICP is the 3D points of each building buildings and the rotations and translations from the conformal transformation as initial values.The ICP algorithm further reduces the mismatch of the two data sets and thus improved the accuracy of the initial transformation based on data uncertainty.Therefore accurate transformation parameters between LiDAR and photogrammetric DSM point clouds in urban environments can be obtained by combining 3D conformal transformation based on planar surfaces followed by the application of the ICP algorithm on the 3D building points using a combination of building / plane matching approach for the determination of the corresponding planes.Future work will explore the different roof (plane) orientations and the inclusion of linear features.
Figure (1b) and (1c) represents the bigger lower plane roof in Figure (1a) while the lower class represents the smaller upper plane roof in Figure (1a).As one can see, the 3D histogram of the building from both the LiDAR and the photogrammetric data are very much alike.Hence the 3D histogram of building is a good measurement for building matching between LiDAR buildings and DSM buildings.Each building can be represented by a 3D histogram.Therefore the building matching problem can be transferred to a 3D histogram matching problem.
Figure 1.3D histogram of a building Since the 3D histogram matching heavily depends on the segmentation results, we define two indices of the 3D histogram that are relatively invariant to over-segmentation and for the cases of over smoothness.The two indices that are used in our building matching algorithm are: the scalar part that equal to the real number d .Unit quaternions ||q||=1 provide a convenient is estimated from the conjugate plane centroids m as:

Table 1 .
Check planes normal vectors differences