MULTI-DIMENSIONAL QUALITY ASSESSMENT OF PHOTOGRAMMETRIC AND LIDAR DATASETS BASED ON A VECTOR APPROACH

The production of 3D city models requires the reconstruction of individual 3D building models. As the performance of data acquisition methods improves, the quality evaluation of building models in 3D has become an important issue. The main objective of the presented work is to introduce a multi-dimensional approach for assessing the quality of 3D building vector models. This approach performs assessments in 1D, 2D and 3D by comparing calculated building models to their reference. For 1D assessment, homologous points in two buildings to be compared are analyzed. For 2D assessment, homologous planes enter in the evaluation process. Quality of the planes under study is assessed by calculating a set of indices in vector format. For 3D assessment, building models are considered as one object. Quality of the buildings is assessed by calculation of vector volumetric quality factors. These factors require the not trivial calculation of vector intersection volumes which calculation is presented in the paper. Intersection volume is defined by superimposing the building model to be tested with the reference one. The multi-dimensional vector assessment approach has been applied to evaluate the building models produced with three different reconstruction processes created from different types of datasets. The datasets are obtained by photogrammetry (UltraCam-X and Zeiss LMK cameras), by LiDAR, and also by integration of photogrammetric and LiDAR datasets. The 1D, 2D or 3D assessment approach allows highlighting the source of deviations in the tested buildings. The error budget affecting the final product is not only composed of errors due to the reconstruction algorithm. Also errors due to the quality of the raw data, the processing of LiDAR data, of aerial data and the shape of the produced buildings should be considered.


INTRODUCTION
Photogrammetric and LiDAR data are used since many years for the 3D reconstruction of objects such as buildings.These applications need accurate data and methods in order to produce good results.Quality assessment is critical for 3D data production and is important for several reasons.Firstly, it may give important information about deficiencies of an approach and may take place to help in focusing a further research activity.Secondly, accuracy evaluation is needed in order to compare the results of the different approaches and to convince a user (Schuster and Weidner, 2003).Several methods were presented in order to evaluate photogrammetric and LiDAR datasets.Calculation of Root Mean Square Errors (RMSE) for analyzing the precision of the complete 3D building is an interesting process.It shows the shifts between reference and test models in X, Y and Z directions.(Grussenmeyer et al., 1994) proposed statistical techniques in order to calculate RMSE by point and line based assessment.Several approaches using quality factors for quality evaluation of building models were introduced ((McKeown et al., 2000), (Ragia, 2000)).Based on these, a further approach was developed, which introduced alternative quality measurements (Schuster and Weidner, 2003).All accuracy assessments processes include three fundamental steps (Congalton and Kass, 2009).Firstly, the accuracy assessment sample should be designed.The sampling design issues are similar to those traditionally addressed by survey sampling methodology: how to choose sample in a cost-effective and at the same time statistically rigorous manner?Application of basic sampling designs such as simple random, stratified random, systematic and cluster have been summarized in (Congalton and Kass, 2009).Secondly, data must be collected for each sample; and finally, results must be analyzed.Because high quality reference data are difficult and expensive to obtain, reference model may be considered as error-free, more accurate than the test model, or with the same accuracy as the test model (Meidow and Schuster, 2005).Once the reference data are available, the assessment process can start.Related works were done by the Photogrammetry and Geomatics Group at INSA-Strasbourg.In the context of assessing the quality of planes detection in a 3D building reconstruction process based on LIDAR data, several solutions have been suggested (Tarsha-Kurdi et al., 2008).For evaluating the quality of geometric facade models reconstructed from TLS data, (Landes et al., 2012a) suggested the use of quality factors and RMSE calculations.Also, the evaluation of characteristic planes extracted from digital airborne sensors have been published in (Mohamed and Grussenmeyer, 2011).A new approach has been proposed for assessing 3D building models based on vector volumetric quality factors (Landes et al., 2012b).The main objective of the presented work is to introduce a multi-dimensional approach for evaluating the quality of 3D building models reconstruction.The approach used in this research is vector based.

1D assessment
Point accuracy assessment allows to evaluate fully 3D geometry datasets by comparing points to points (Rottensteiner and Briese, 2002).Calculations are performed to compare two different 3D building models.This is done by computing RMSE based on the deviations between both models (reference and test), in X, Y and Z directions.Deviations are not calculated between homologous nodes, but between centers of gravity of homologous planes that compose the tested and respectively the reference building.

2D assessment
Our method to evaluate 2D data is based on the comparison of 3D planes of two building models (reference and test) by calculating a set of indices.This approach uses well known quantities, as mentioned in (McKeown et al., 2000), (McGlone and Shufelt, 1994), (Ragia, 2000), (Schuster and Weidner, 2003) and (Landes et al., 2012b).These indices are namely the detection rate (ρd), the branch factor (ρb), the quality rate (ρq), miss factor (ρm), and false alarm rate (ρf ) as summarized in Table 1.The principle idea of these indices is based on the relation between the reference surface area (Ar) and the tested surface area (At) as shown in Figure 1(a).

Quality factor
Explanation Detection rate: it is the ratio between the ρd = Ar∩At Ar intersection area between two planes and the ρd ∈ [0 : 1] reference plane.ρd = 1 means that the calculated polygon is perfectly superposed Quality rate: it is the ratio between the ρq = Ar∩At Ar∪At intersection area between two planes and the ρq ∈ [0 : 1] union of two planes.ρq = 1 means that both polygons are perfectly superposed.Branch factor: it is the ratio of the part of the ρb = At\Ar Ar∩At reference polygon which is not included pb ≥ 0 in the polygon under study and the intersection of the two polygons.Miss factor: Ratio of the part of the polygon ρm = Ar\At Ar∩At being evaluated which is not included in a ρm ≥ 0 reference polygon and the intersection of the two polygons.False alarm rate: it is the ratio of the part of ρf = At\Ar Ar reference polygon which is not included ρf ≥ 0 in the polygon under study compared to the area of the reference polygon.
Table 1: Quality indices based on surface areas ratios.

3D assessment
Considering a building as one object to be evaluated, volume comparisons are rather appropriate than surface areas comparisons.Therefore, the quality of buildings is assessed by the calculation of volumetric quality factors.The quality factors in 3D are deduced from the 2D quality factors and depend on ratios of volumes.These factors take into account the volume of the intersection as well as the union volume of two vector buildings.
A first experiment has been introduced in (Landes et al., 2012b).
The principle idea of these indices is based on the relation between the reference volume (V r) and the tested volume (V t) as shown in Figure 1(b).
Equations in Table 2 detail the volumetric quality factors.Satisfying results are reached when the value of V ρd and V ρq is close to 1, and the three others are close to 0.
Quality factor Explanation Volumetric detection rate: it is the ratio the two buildings and the volume of the reference building.Volumetric quality rate: it is the ratio V ρq = V r∩V t V r∪V t between the parts which are common to V ρq ∈ [0 : 1] both volumes and the union of two volumes.Volumetric branch factor: it is the ratio of the part of reference building volume V ρb ≥ 0 which is not included in the building volume under study and the intersection of two volumes.Volumetric miss factor: it is the ratio of V ρm = V r\V t V r∩V t the part of the volume being evaluated V ρm ≥ 0 which is not included in reference volume and the intersection of two volumes.Volumetric false alarm rate: it is the ratio of the part of reference volume which V ρf ≥ 0 is not included in volume under study compared to the volume of reference polygon.
Table 2: Quality indices based on volume ratios.

GEOMETRIC COMPUTATIONS
The operation of determining the union and intersection of areas (or volumes) when two models must be compared is a topic of geometric computation.This section describes the algorithms leading to calculate union and intersection of areas/volumes which come into consideration when a model to be assessed is compared to a reference model.First of all the areas and volumes of the models under study must be calculated.

Surface area computation
For the area computation in vector form, we use the formula based on Green's theorem (O'Rourke, J, 1998) and given in Equation (1).
Where n is the number of vertices, Xn+1 = X1 and Yn+1 = Y1.X and Y are the coordinates of vertices points of the polygon numbered in ascending order.

Intersection area computation
The operation of determining the intersection area of two vector polygons is performed in two steps: the detection of the points located inside the polygons and the detection of lines intersections.
Step1: classification of the points located inside the polygons As shown in Figure 2(a), points 3 and 5 are "inside points".In Figure 2(b), all points of the red polygon are classified as inside points.Matlab built-in function uses a simple and commonly used technique for point-in-polygon detection, and works as follows.Assuming the polygon is defined by n points in an array P , this algorithm computes the summation of angles between the query point and every pair of points defining each edge of the polygon (i.e. the angle is formed by the P [n] point, query point, and P [n + 1]).If this summation computes to 2π (or near 2π within some tolerance), then the point is inside the polygon.If the summation computes to zero (or near zero) then the point is outside the polygon.
Step2: detection of points at the intersection between lines For example, in Figure 2(a), the point of intersection between line (2, 3) and line (5, 6) is a "point of intersection" to consider, as well as the intersection point of lines (4, 3) and (5, 8).However, the intersection between lines (4, 3) and (6, 7) is not an intersection point to consider.That's why it must be checked if the intersection point lies on the edges of both polygons.This check is done by computing the distance between the intersection point and the two points of the line.If the maximum of the two distances is shorter than the edge length, the point of intersection belongs to the edge.Then the resulting intersection shape can be calculated based on the coordinates of its vertices (see Equation ( 1)).
Figure 2: Intersection area calculation for vector polygons in 2D.

Volume calculation
In order to compute the volume of 3D objects, one has to solve two tasks; firstly determine the convex hull of the given boundary, then, calculate the volume of the resulting 3D polyhedron.Convex hull is the boundary of a closed convex surface generated by applying Delaunay triangulation on the corner points.In three dimensions, the convex hull corresponds to a closed polyhedron.Convex hull calculation is a hard process in computational geometry (Barber et al., 1996).A large class of algorithms that compute the exact volume of a convex object is based on triangulation methods (Büeler et al., 1997).The result of convex hull calculation for a gable roof building is shown in Figure 3(a).The convex hull of a set of points in two or three dimensions is given by Matlab built-in function (convhull in 2D or convhulln in 3D) as presented in Equation ( 2).These functions use meshed objects for storing and displaying polyhedra.The faces of such polyhedra are triangles.
In 3D, the boundary of the convex hull, K, is represented by a triangulated 3D object.It is a set of triangular facets in face-vertex format that is indexed with respect to the point array.Each row of the matrix K represents a triangle.The volume, V, bounded by the 3D convex hull can optionally be returned by convhulln.
In a first step, the meshed model is computed.It is defined by the input points.The sum of the volumes composing the meshed model equals to the volume of the convex object.

Intersection volume calculation
The calculation of volume of intersection between two 3D models in vector format is more accurate than in raster format, but also much more complicated.We propose an algorithm allowing to simplify this process.Our method consists of extracting the vertices of 3D intersection volumes.The flowchart in Figure 4 shows the process of intersection volume calculation.The following steps are performed: Step 1: detection of inside points Extraction of the first group of points is defined by searching about vertices of the reference points located inside the model to assess.This can be done by applying the function "convhulln" to the points of the model to assess.The points that are inside the reference model are located on the positive side of the plane normals of all of the faces.The result of this process is shown in Figure 3(a) where inside points are in red and outside points in green color.
Step 2: creation of boundary lines and their intersection with planes The second group of points describing the intersection shape can be determined by calculating the intersection of the lines composing the reference model with all planes that are composing the model to assess.This process can be achieved firstly by separating the edge lines of each plane of the reference model.Then, the duplicated lines are cancelled in order to avoid repeating the same process.After that, the intersections of all lines with all planes of the model to assess are calculated.In order to check if the resulting point is located on the edge line (and not on the extension of the edge line) and simultaneously belongs to the face of model to assess, two tests should be made.Firstly, we test if the point is placed on the edge line by distance computation as shown in section 3.2.A second test is achieved by looking for the "points inside a polygon" in 2D (in the frame of the intersected plane), as explained in section 3.2.
Step 3: repetition of steps 1 and 2 Steps 1 and 2 are repeated by replacing the reference model by the model to assess for the process leading to edge line creation.
Then the intersection between lines of the model to assess with all the planes of reference model is performed.Finally, as a result of these steps, the coordinates of vertices of the intersection shape are determined and the volume of this shape can be calculated.Figure 3(b) shows the intersection shape (filled with red color) obtained by the intersection of a reference model and a test model.

3D BUILDINGS RECONSTRUCTION
In this paper, three methods for 3D buildings reconstruction are presented.

Test site and datasets
The study area is located in the city of Strasbourg, France.Digital aerial images from UltraCam-X (4 images) and frame Zeiss LMK (6 images) of the same area were available.Table 3 shows characteristics of the photogrammetric data.Moreover, LiDAR data of the same are has been captured in 2004 (Table 4).14 ground control points were measured by GNSS.The digital aerial photographs have unknown imaging orientation parameters.Using bundle block adjustment, the exterior orientation parameters of the images have been calculated by KLT software.In our application, the camera information was taken from the calibration sheet given by the camera owner.3D reference buildings models have been created based on the photogrammetric processing of images acquired with UltraCam-X stereo-pairs.After relative and absolute orientation of the images, an accuracy of about 16 cm in X,Y and 25 cm in Z can be estimated for a point digitized in the stereo-pairs.It is satisfactory, considering that the accuracy of the LIDAR point clouds used here is lower (around 30 to 40 cm in X, Y, Z).So, it has been decided that the 3D buildings from UltraCam-X stereo-pairs will be used as references for assessing the 3D buildings reconstructed a) from Zeiss LMK stereo-pairs (75 buildings), b) from LiDAR data (8 buildings), and c) from integration of LiDAR and UltraCam-X stereo-pairs (26 buildings).For illustration purposes, only 8 buildings are shown in Figure 7.The reconstructed buildings in the test site have three types of roofs: flat, hip and gable roofs.

3D models from aerial images
The geometry of objects (roofs, walls, and footprints) are extracted from multiple images.The flowchart of the semi-automatic approach is depicted in Figure 5.The first step of building reconstruction is roof digitizing.Then, the projection of these points onto the ground is done in order to obtain footprints and thus to create the walls.Finally, planes of faces and footprints are created.
The reconstruction approach is based on the assumption that: (a) every solid object can be described by a decomposition of its boundaries; (b) the walls are vertical and reach either the ground or another surface of the constructed model.The wall faces can be constructed using the outlines of the roofs (no facade details).Therefore it is not necessary to digitize the footprints of the buildings.In this work, we restrict our study to simple polyhedral models.Figure 7(b) presents 8 of the 75 samples of building models reconstructed in the test site (in yellow colour) as well as the corresponding reference (in red).

3D models from LiDAR and aerial images
The semi-automatic method proposed by (Zhang et al., 2011) is based on the complementarities of airborne LiDAR and optical imagery (UltraCam-X).It consists of 4 steps (Figure 6).
Firstly, the building is decomposed into several primitives.The primitives parameters are measured manually on LiDAR and aerial imagery, such as length, width, height, orientation and translation of the primitive.These measurements can be used as constraints or initial values in the following optimization procedure.Secondly, features primitives are selected on the imagery.Corners are detected on the optical imagery, and planes are selected in the LiDAR point cloud.These features are used as observations in the following optimization procedure.Thirdly, the algorithm optimizes primitives parameters by the constraints of LiDAR point cloud and imagery.Based on the type and parameters of primitives, the 3D coordinates of the features primitives, such as corners, can be calculated.These 3D coordinates will be used as computed values in the next iteration.Finally, 3D building models are produced.
Figure 6: Flowchart of the reconstruction process using LiDAR data and aerial images (Zhang et al., 2011).
This method has been applied to 26 building models of the test site.Figure 7(c) shows 8 of 26 samples of reconstructed buildings (in green) and their reference buildings (in red).

3D models from LiDAR datasets
In this part, a model-driven building reconstruction method using airborne LiDAR data is presented.This method has been carried out by Yong Xiao from the Chinese Academy of Sciences (China).This semi-automatic reconstruction process comprises 3 steps.At first, the point cloud covering the building is segmented to isolate building roofs.Then, a topological graph is constructed to recognize the shape of the buildings.Finally, once simple roof types are determined, building models are reconstructed with predefined models.Outlines of the buildings are first estimated with the minimum area bounding rectangle while the other key vertices and segments are obtained through the roof-topology graph.More details are given in (Verma et al., 2008).Figure 7(d) shows the results.The buildings reconstructed from LiDAR datasets are in cyan colour and their reference in red.

ASSESSMENT RESULTS
In this section, we assess the 3D vector building models reconstructed previously, by applying the multi-dimensional assessment approach suggested in section 2.

1D assessment
The reference building models are the models reconstructed from UltraCam-X (see section 4.2).Deviations are calculated between centers of gravity of homologous planes that compose the tested and respectively the reference building.

2D assessment
The statistics of the 2D quality indices calculated for assessing the faces of the building models are summarized in Table 6.For models created from aerial images, the mean values of ρd and ρq are about 0.9 and the other three indices are close to 0. The worse values obtained for the other models are explained by the high RMSE in Y and/or Z-directions (Table 5).This means that the 3D surface of building models extracted from stereo-pairs are close from each other.Also, the models reconstructed from LiDAR or integration of LiDAR and aerial images are less accurate than the models reconstructed from aerial images alone.However, the mean values of quality indices can not be considered alone.In order to evaluate the building reconstruction quality in

3D assessment
The results of the quality analysis based on volumetric quality factors are shown in Table 7.The mean values of V ρd and V ρq are about 0.9, while the other three indices are close to 0 for 75 models of aerial images.The mean values of volume quality indices, V ρd and V ρq are about 0.8.The other three indices are about 0.1 for 8 models created from LiDAR datasets and 26 models of integration of optical and LiDAR datasets.This means that the models reconstructed from aerial images are more accurate than the other models.The 3D assessment provides more information than the 2D assessment, because it takes into account the shift which might exist in the third dimension between the two 2D surfaces.This occurs for the quality indices obtained for the models reconstructed from the LiDAR dataset.All tests applied on the building models reveal that a systematic error affects the Z coordinates of the LI-DAR data used here.This vertical shift has already been observed in a study where ALS and TLS data were combined (Boulaassal et al., 2011).

CONCLUSIONS AND FUTURE WORK
In this paper, three semi-automatic methods for 3D building reconstruction in vector format have been mentioned and carried out on the same test site.The quality evaluation of these models has been achieved by applying the proposed multi-dimensional quality assessment approach.This approach considers the accuracy of the 3D building models based on the comparison of points in 1D, of surfaces in 2D, and of volumes in 3D.1D assessment gives an overall idea about the reliability of the reconstructed models.2D assessment checks the superimposition of faces, despite its dependency on the size of the polygons to be compared.3D assessment compares the buildings in 3D through the comparison of their volumes intersection.It is appropriate for detecting the direction of errors (shifts in X, Y, and Z or rotations).This multi-dimensional approach is suitable for 3D building vector models created from aerial images and/or LiDAR datasets.Future researches will focus on the extension of this approach to more complex building models.In this context, it will be focused on the benefits of using raster assessment approaches.

Figure 1 :
Figure 1: Relationship between reference and tested models in 2D (a) and 3D (b).

Figure 4 :
Figure 4: Flowchart of the method of intersection volume calculation.

Figure 5 :
Figure 5: Flowchart of 3D building modeling from aerial images.
(a) Part of the models under study (b) Models obtained from aerial images (8 of 75 buildings).(c) Models obtained from aerial images and LiDAR datasets (8 of 26).(d) Models obtained from LiDAR datasets.

Figure 7 :
Figure 7: Reference and tested building models.
detail and to analyze the values of the 2D quality factors, one should check the quality indices for each building separately.Moreover, surface metrics are affected by the building size.Small buildings generally lead to bad results regarding the 2D quality indices.Reconstruction ρd ρq ρb ρm ρf from Aerial images 0.938 0.891 0.089 0.085 0.062 Images & LiDAR 0.867 0.788 0.177 0.154 0.120 LiDAR only 0.840 0.711 0.219 0.250 0.189Table 6: 2D quality indices.

Table 4 :
Characteristics of the LiDAR datasets.

Table 5 :
Table5presents the RMSE results obtained.The models reconstructed from aerial images give better results than the other methods.The models reconstructed from integration of LiDAR and aerial images give high RMSE in Z-direction.Worse results are obtained for models reconstructed with LiDAR data only.The error budget is not only composed of errors due to the reconstruction algorithm, but also of errors coming from the raw data.For instance, low point cloud density and errors due to the georeferencing of the LIDAR and the optical data affect the final results.RMSE calculated on gravity centers of homologous planes.