RECONSTRUCTING 2-D / 3-D BUILDING SHAPES FROM SPACEBORNE TOMOGRAPHIC SYNTHETIC APERTURE RADAR DATA

In this paper, we present an approach that allows automatic (parametric) reconstruction of building shapes in 2-D/3-D using TomoSAR point clouds. These point clouds are generated by processing radar image stacks via advanced interferometric technique, called SAR tomography. The proposed approach reconstructs the building outline by exploiting both the available roof and façade information. Roof points are extracted out by employing a surface normals based region growing procedure via selected seed points while the extraction of façade points is based on thresholding the point scatterer density SD estimated by robust M-estimator. Spatial clustering is then applied to the extracted roof points in a way such that each roof cluster represents an individual building. Extracted façade points are reconstructed and afterwards incorporated to the segmented roof cluster to reconstruct the complete building shape. Initial building footprints are derived by employing alpha shapes method that are later regularized. Finally, rectilinear constraints are added to yield better geometrically looking building shapes. The proposed approach is illustrated and validated by examples using TomoSAR point clouds generated from a stack of TerraSAR-X high-resolution spotlight images from ascending orbit only covering two different test areas with one containing relatively smaller buildings in densely populated regions and the other containing moderate sized buildings in the city of Las Vegas. * Corresponding author. This is useful to know for communication with the appropriate person in cases with more than one author.


INTRODUCTION
With data provided by modern meter-resolution SAR sensors and advanced multi-pass interferometric techniques such as tomographic SAR inversion (TomoSAR), it is now possible to generate 4-D (space-time) point clouds of the illuminated area with point density of approx. 1 million points/km 2 .However, due to side looking geometry, these point clouds exhibit much higher density of points on building façades in contrast to nadir looking LiDAR geometry (typically used for object reconstruction).Moreover, temporally incoherent objects such as trees cannot be reconstructed from multi-pass spaceborne SAR image stacks and provide moderate 3-D positioning accuracy in the order of 1m as compared to airborne LiDAR systems (around 0.1m).Despite of these special considerations, object reconstruction from these high quality point clouds can greatly support the reconstruction of dynamic city models that could be potentially used to monitor and visualize the dynamics of urban infrastructure in very high level of details.Motivated by these chances, earlier approaches have been proposed to reconstruct building façades from this class of data.E.g., experimental results provided in (Zhu, 2014) and (Shahzad, 2014) over smaller and larger areas demonstrate that façade reconstruction is an appropriate first step to detect and reconstruct building shape when dense points on the façade are available.In particular, when data from multiple views e.g., from both ascending and descending orbits, are available, the full shape of buildings can be reconstructed using extracted façade points.However, there are cases when no or only few façade points are available.This happens usually for lower height buildings and renders detection of façade points/regions very challenging.Moreover, problems related to the visibility of façades mainly pointing towards the azimuth direction can also cause difficulties in deriving the complete structure of an individual building.These problems motivate us to reconstruct 2-D/3-D building shape (footprint) via roof point analysis.In this paper, we propose solutions to the following two cases: 1) When only roof points are available, i.e., no or very few façade points exist; and/or 2) Data is acquired from one orbit, e.g., ascending orbit only.
In such a case, obvious occlusion due to side looking SAR geometry renders façade points to be available only from one side.To reconstruct the other side of the building, information related to roof points thus needs to be exploited.

APPROACH OVERVIEW
Figure 1 presents the complete processing chain of the proposed approach.Building points are sequentially extracted out by first extracting façade points followed by extraction of roof points.Façade points are utilized to reconstruct sides of the buildings that are visible in the data.Later roof points are spatially clustered and respective façade-roof pair is identified.Based on the availability of façade points, i.e., case 1 where no façaderoof pair is found and only roof points are available or case 2 where façade-roof pair is found, initial building footprints are reconstructed by alpha shapes method.Usually such reconstructed building outlines are quite rough and therefore cannot be directly used.Therefore they are regularized to yield better geometrically looking building shapes.The paper is organized as follows: In the following two sections 3 and 4, we explain the processing steps in detail.In sections 5 and 6, the experimental results on two different test areas obtained from the TomoSAR point cloud generated from a TerraSARX high-resolution spotlight data stack (ascending orbit only), are presented and discussed.Finally, in section 7, a conclusion about the proposed approach is drawn, and future perspectives are discussed.

Façade points extraction
As already mentioned, side-looking SAR geometry enables rich number of points on building façades.TomoSAR points therefore when projected onto ground (xy plane) exhibit higher scatterer (point) density SD in vertical façade regions as compared to nonfaçade regions.It is mostly true due to the existence of strong corner reflectors, e.g., window frames on the building façades.Taking this fact into account, we earlier proposed an approach that robustly estimates the SD while incorporating the façade geometry (Shahzad, 2014).The basic idea of the approach is as follows:


For each 3-D TomoSAR point p, points within its local neighborhood c v are used for SD estimation.c v includes all of those points that lie inside a vertical cylinder centered at p.


To emphasize the building façades, we incorporate façade geometry in estimating SD, i.e., we estimate the direction of the local neighborhood via line fitting using robust Mestimator.


The estimated line describes the main principal axis of the cylindrical footprint of the local neighbourhood.
Orthogonal distance for every point in c v is then calculated from the principal axis (shifted to the point p) and the subset of points d v having distances less than d are taken as "inliers" and used in SD estimation.
SD for each point is thus defined as the number of points within a directional (cylindrical) neighborhood window divided by the area of the window: number of points in Area of where dc vv  but includes only those points that lie close to the principal axis of points in c v .
Façade points are then extracted out by applying a two step procedure:  Detection of probable façade points by applying soft thresholding to the SD estimated via aforementioned Mestimator based directional method;  Rejection/removal of false positives by utilizing 3-D surface normals estimates.
As shown in (Shahzad, 2014), the two step approach allows us to robustly extract façade points over a large area where both high and low buildings are present.

Roof points extraction
Buildings are first extracted directly from the unstructured TomoSAR point cloud.The rationale behind the approach is the assumption that buildings and other man-made structures are elevated objects within the vicinity of their surrounding region.The procedure thus first determines the transition regions/points in the whole point cloud.The above procedure results in too many points occurring on transition regions although most of these points lie on building boundaries.Therefore instead of testing neighbours of each transition point, density based clustering procedure is adopted which spatially cluster these transition points.Robust 3-D surface normals, based on minimum covariance determinant (MCD) method, are then estimated for the remaining (non-transition) points and later 3-D region growing procedure is adopted by choosing appropriate seed points.Neighbouring points are added into the region based on the similarity of their surface normals.A minimum height constraint, computed from height statistics of the neighbouring points belonging to each transition cluster, is also imposed in the growing procedure to cope for inaccurate positioning accuracies of 3-D points.The procedure is adopted for all the transition clusters and union of all the points in the grown regions thus results in extraction of building roof points.

Selection of seed points
If we denote the transition clusters as 1,...,  Seed points extracted from the above procedure can be used in the region growing procedure.However since we are extracting them based on the object elevation, there might also be points belonging to transition clusters that are not part/ near buildings but rather occur from other originating sources e.g., lamp post etc.Thus growth of the seed points extracted from such transitions clusters need to be restricted in the region growing procedure.This is resolved by introducing two constraints, minimum height constraint

Computation of surface normals
The similarity criterion used by the region growing procedure for adding points in the cluster is the surface normal vector locally computed at the point of interest.Use of surface normals allow extraction of points belonging to flat or polyhedral roof  (Hoppe, 1992).This implies that the surface normals can be directly estimated for each 3-D point via eigenvalue/eigenvector analysis of 3-D (i.e., 3x3) covariance matrix Thus starting from a seed, all points that share similar normal orientation within its neighbourhood are added into the region provided they satisfy already mentioned two constraints min h and min  .The procedure repeats itself until all the seed points have been utilized resulting in the set of points belonging to building roofs.

Façade reconstruction
Extracted façade points are further segmented to points belonging to the same individual façade as follows:  Extracted points are coarsely clustered by a density based connectivity approach as proposed by (Ester, 1996);  Surface normals are computed locally for each point and the mean shift algorithm is used for clustering points having smaller angular difference in feature space (Gaussian image GI) into one cluster (Liu, 2008) (Cheng, 1995);  Previous step results in clusters of points that have similar normal directions but may be spatially far from each other.
To cope this, spatial connectivity is used for further clustering of points.
Each cluster is further classified into flat or curved surface by analyzing derivatives of the local orientation angle θ (= azimuthal angle of the surface normal).Identified façade clusters in xy plane are then modeled using the general polynomial equation (Zhu, 2014) (Shahzad, 2014): Here i and j are permuted accordingly, p is the order of polynomial, the number of terms in the above polynomial is equal to (p + 1)(p + 2)/2.The coefficients q a are estimated using weighted total least squares (WTLS) method where total least squares is utilized to cope for localization errors of TomoSAR points in both xy directions and the weight of each point is assigned equal to its corresponding SD.
After estimation of model parameters, the next step is to describe the overall shape of the building footprint by further identifying adjacent façades pairs and determining the intersection of the façade surfaces.The adjacency of façades is usually described by an adjacency matrix AM that is built up via connectivity analysis (Zhu, 2014).Identified adjacent façade segments are used to determine the vertex points (i.e., façade intersection lines in 3D) by computing the intersection points between any adjacent façade pair.Determination of these intersection points can sometimes become difficult if the transition points (i.e., points occurring at the transition region of two adjacent façades) are segmented as isolated small clusters rather than part of the corresponding adjacent façade segments.
As a consequence, it gets complicated to find a legitimate adjacent façade pair from which intersection points should be computed.To resolve this issue, such cases are first identified and then the intersection point is computed from the two largest segments only.The computed vertex points (i.e., the intersection point of the two adjacent façade pair and the other "open" ends of these façades) along with their estimated model parameters completes the façades reconstruction procedure.

Identification of legitimate façade-roof pair
Prior to reconstruct the footprint by tracing the boundary of roof points, the reconstructed façades information is incorporated.This is done by first identifying the legitimate (reconstructed) façade-roof pair by performing following steps:


Apply density based spatial clustering to the extracted roof points to cluster spatially connected roof points.Each cluster is thus taken as an individual building structure.


Compute the midpoint of each reconstructed façade and then search in orthogonal direction for roof clusters.


If the distance to the nearest roof cluster is within come reasonable limit d, merge the reconstructed façade points into the corresponding roof cluster.
Thus if there are no façade points available i.e., Case 1, we reconstruct the building shape based on roof points only.For the other case, we incorporate the façade points together with the roof points to determine the overall shape of the building footprint.Figure 3 shows the procedure for incorporating façades in the reconstruction procedure.

Coarse building footprint
Reconstruction of building shape is then initially obtained by applying alpha shapes (generalization of convex hull) around each segmented building (Edelsbrunner,1994).The output of alpha shape (or α-shape) algorithm results in vertices describing the coarse 2-D polygonal footprint of the building.The shape of the reconstructed building footprint however depends on the particular value of α.Larger values of α tends to describe the convex hull around the points.Thus setting α large enough makes it difficult or even impossible for the algorithm to determine the building boundary having concave shape e.g., an L-shaped building.An appropriate value of α thus need to be empirically found or estimated from the data.A good value of α that produce reliable building shape, including smaller structures, may be chosen as the twice of the mean Euclidean point distance among building points (Dorninger, 2008).

Refinement of the building footprint
Due to lesser point density of TomoSAR points, alpha shapes only define the coarse outline of an individual building and therefore the resulting polygons are irregular and contain shorter line segments that need to be regularized.
To resolve this, a two step regularization procedure is adopted:

1) Refinement of alpha shapes vertices via mean orientation estimation
The coarse reconstructed building footprint via alpha shapes algorithm is refined (or regularized) by computing mean orientation at each vertex point (Dorninger, 2008).If we denote 1,..., in z  as the 2-D vertices of the initial alpha polygon (with n equal to the total number of vertices of the polygon), mean orientation θ at each vertex i z is computed as where ij zz  is the slope computed using vertices i z and j z .
Starting from a vertex point, θ for the next consecutive vertex is checked.If the difference is less than certain value denoted as m  , it is removed.The removal procedure continues till a vertex is found whose mean orientation is greater than m  .This vertex is retained and the procedure for removing vertices is again started but this time the mean orientation of the next consecutive vertex is tested with the previously retained vertex.
The procedure finishes when the algorithm reaches at the same vertex where the refinement/regularization procedure began.

2) Addition of rectilinear constraints
Subsequently, rectilinear constraints are added to derive correct and better looking geometric building shapes.This step is based on the assumption that buildings are mostly composed of two dominant directions that are orthogonal to each other.Following steps are performed to obtain rectilinear building footprint (see Figure 5):

Estimation of principal orientation
Minimum bound rectangle (MBR) can be employed to determine the dominant/principal orientations of the building footprint.This can be done by building up MBR around 2-D building boundary points.The two orthogonal axes of the MBR then provides the direct estimate of the desired dominant orientations.(Sun, 2013) and (Arefi, 2013) have adopted MBR to estimate the building footprint orientation.Commonly, MBR is computed by a method known as "rotating calipers" (Toussaint, 1983) which is based on the theorem, proved in (Freeman, 1975), that any minimum area bounded rectangle is collinear with at least one of the sides of the convex hull.The convex hull is therefore first computed and later the bounding rectangle is sequentially computed by rotating the convex hull polygons in a way such that each side of the convex hull becomes parallel to x-axis.In each rotation, the area of the minimum bounding box around 2-D points is computed and the rotation angle that provides the minimum bounding area is used to determine the vertices of the desired MBR.
Use of MBR provides reasonable estimates for the principal orientation of the buildings.However, in some cases, it is possible that MBR suffers in accurate determination of building dominant directions.An example is shown on L-shaped buildings with noisy points in Figure 4. Due to less 3-D positioning accuracy of TomoSAR points, such situations often arise and other methods for estimating the principal orientation are therefore needed to be explored.
An alternative to MBR, several researches e.g., (Sampath, 2007) (Dorninger, 2008) and (Jarzabek-Rychard, 2012) have estimated the principal orientation by first refining the initial boundary polygons and later considered the longest polygonal edge as an estimate for the building footprint orientation.This idea seems to perform better for irregularly derived building outlines.In our work, we adopted the similar idea but instead of directly taking the longest edge of the refined polygons (i.e., after smoothing the coarse reconstructed building footprint), we propose modifications that improve the robustness of the estimated dominant directions.The proposed idea of computing the principal orientation of the building footprint is as follows: Define a vector q and compute the angular difference γ of each edge i e with respect to q.All the edges having 45

EXPERIMENTAL RESULTS
To validate our building extraction method, we tested the algorithm on two different sites in the city of Las Vegas with one containing relatively smaller buildings in densely populated regions, Testsite 1, while the other containing moderate sized buildings, Testsite 2. TomoSAR point clouds for both these sites are generated from a stack of 25 TerraSARX high spotlight images from ascending orbit only using the Tomo-GENESIS software developed at the German Aerospace Center (Zhu, 2013).The number of TomoSAR points in both the area of interests are about 0.45 million and 0.4 million respectively.The area of each scene is around 1 km 2 each.Figure 6(a) and Figure 7(a) shows the optical images of our test areas.v as the one whose height is maximum.Surface normals based region growing begins the growing procedure based on similarity of surface normals (the angular deviation used for adding a point is set to 10 degrees).The procedure is followed for every seed point with addition of neighboring points in the region subject to two constraints min h (adaptively computed as explained in the section 2.2) and min  set to 5. Roof points are then extracted by taking union of all region grown points extracted from all seeds.The results of complete reconstruction procedure are depicted in Figure 7. Figure 7   angular difference from the ground surface.Finally, appropriately handling both cases 1 and 2 leads to the overall reconstruction of the building outlines.Figure 7(d) and Figure 7(e) shows the final reconstructed building shape in 2-D and 3-D.

DISCUSSION
The algorithm for extracting façade points has already shown promising results as depicted in (Zhu, 2014) and (Shahzad, 2014).The approach of extracting roof points seemingly provides promising results over the two different area of interests with different scales of the buildings.The overall approach, in general, is automatic but still requires tuning of parameters for extraction and reconstruction modules.Results shown on both test sites are obtained by same parameter settings.Use of 5m radius setting for spatially clustering and neighborhood selection renders the algorithm to separate two buildings only if they are at least farther than 5m from each other.Otherwise, the algorithm will merge them into one single cluster.
Another critical parameter used in roof points extraction procedure is the minimum height constraint min h that restrict addition of smaller height points during region growing process.Setting a low min h may cause the algorithm to fail as in that case many non building points will also be added into the region.This can happen for cases where the seed point is surrounded by flat terrain e.g., a parking lot, roads etc.Such situation can be avoided either by reducing the surface normal threshold used as similarity measure or setting a maximum size limit to each grown cluster.Alpha shapes method provides good initial estimates of building outlines.The value of α effects the shape of this initial polygon.Since with alpha shapes, it is not only possible to extract the outer outline of the building but instead with lower values of α, the method can also be employed to extract inner polygons of certain shaped buildings that are closed but possess non building parts in between e.g., a doughnut shaped building.Thus for buildings with no inner polygonal region, lower values of α can result in more than one polygon with one outer and the rest inner polygons.Moreover, there also might be case where the outer and inner polygons share one common vertex.To avoid these situations, value of α should not be set too small.After empirical testing, we have determined a good value of α = 5m (for our data) provides relatively good initial coarse outlines.Refinement of the initial alpha vertices is done by computing orientation at each vertex point.The mean orientation threshold m  used for merging two vertices into one group is set to 10 degrees.m  = 0 results in the original alpha polygons i.e., no refinement or Setting too high value for m  may however result in over refinement or smoothing.

CONCLUDING REMARKS AND FUTURE OUTLOOK
We have presented an approach that only utilized unstructured TomoSAR point clouds to reconstruct 2-D/3-D building shapes.The proposed approach analyze both façade and roof points to reconstruct the overall shape of the building footprint.The approach allows for a robust reconstruction of both higher façades and lower height buildings, and hence is well suited for urban monitoring of larger areas from space.The depicted results validate the effectiveness of the proposed approach.The proposed approach is automatic but parametric.The free parameters are set empirically in this work.A further detailed sensitivity analysis of these parameters is therefore necessary.Moreover, we have only presented visual results of the approach.A more detailed evaluation of the algorithm is needed to test its qualitative and quantitative performance.

Figure 1 :
Figure 1: Block diagram of the proposed approach.
from the subset of non transition points and the point having highest height value is chosen as the seed point.The motive behind this step is based on the assumption that the neighbouring points i k v includes both ground and building (roof) points.Thus setting a higher height value point as initial seed ensures that the region growing procedure adds points in the correct direction.Figure2pictorially depicts the procedure.

Figure 2 :
Figure 2: Illustration of seed point selection procedure.Seed point is selected as the highest point among the neighboring points of the transition cluster i k v .
computed via fitting best plane in least sense which is equivalent to performing principal component analysis (PCA) of the points in c v


. However, analysis of eigenvalue/eigenvector via classical PCA may fail to give precise estimate of the 3-D surface normal using TomoSAR point cloud due to presence of outliers and localization errors.Robust estimation of the covariance matrix c v  is therefore needed.To this end, we estimated c v  using robust minimum covariance determinant (MCD) method (Hubert, 2005).The method finds a subset (fraction)  of the data points ic pv  whose covariance matrix has the lowest determinant.The covariance matrix c v  estimated using MCD method is then used to determine the local 3-D surface normal at o p .If we denote a plane which robustly fits the neighbouring points n depicts the local 3-D surface normal at


Estimate dominant directions of the building footprint;  Assigning each smoothed polygonal line (i.e., line between two adjacent vertex points) to one of the principal direction;  Apply rectilinear transformation to each polygonal line by projecting it onto its corresponding dominant axis;  Computing intersection points between adjacent vertices.

Figure 3 :Figure 4 :
Figure 3: Reconstruction sequence for case 2 where both roof points and façade points are available: (a) Red points are the extracted roof points where green and blue points are the extracted façade points for the same building; (b) The red polygon represents the coarse boundary obtained by alpha shapes method after incorporating the reconstructed façades (shown as black lines); (c) Result of refining the coarse building outline by computing mean orientation of each vertex point and removing those consecutive vertices where there is no or very little change in the orientation direction; (d) Final result after adding rectilinear constraints the previous refined building boundary.

Figure 5 :
Figure 5: Illustration of building footprint reconstruction.The estimated dominant orientations, computed by the method described above, are shown by black arrows.Each edge of the coarsely reconstructed footprint is segmented (shown in green and red color) according to their parallelism with respect to the two dominant directions.Later, the final rectilinear shape (shown in blue color) of the building footprint is reconstructed by rectangular projection of each segmented edge onto its corresponding dominant axis.
   are grouped into one category while the edges having 45    into other.Then the length of adjacently connected series of edges (i.e., consecutive edges grouped into one category) is computed.Let us denote the largest length of adjacently connected series in any of the group as e l .The vector q is varied over the then used to determine the principal orientation.This is done by fitting a least squares line among those vertices.The orientation of the fitted line thus describe the main orientation of the 2-D building points.

Figure 6
(c) and Figure 7(b) shows the result of applying the roof extraction procedure on Testsite 1 and Testsite 2 respectively.

Figure 6 :
Figure 6: Roof points extraction results: (a) Optical image of the Testsite 1 containing relatively smaller buildings in densely populated regions in the city of Las Vegas © Google Earth; (b) TomoSAR points in UTM coordinates of the corresponding test image.The height of TomoSAR points is color-coded [unit: meter]; (c) Extracted roof points overlaid onto the optical image.

Cluster roof points by density based clustering approach Identification of façade-roof pair by computing the midpoint of each reconstructed façade and searching for appropriate roof cluster in orthogonal direction Pair found? Fusion of points by combining corresponding façades and roof points Determine building boundary points by employing alpha shapes method Apply line simplification algorithm by computing mean orientation at each vertex point Add rectilinear constraints by determining the principal orientation of buildings 2-D/3-D reconstructed building shape Yes (i.e., Case 2) No (i.e., Case 1)
(c)shows the spatially segmented roof points such that each segment represent an individual building.Black points depict the extracted façade points that are utilized prior to reconstruct the complete building shape.To estimate SD, 5m radius is used to determine the local (cylindrical) neighborhood Probable façade points are then extracted by setting a soft threshold to the maximum of SD histogram value.Soft threshold also results in many false positives which are rejected/removed by retaining only those thresholded points whose surface normals are parallel i.e., 15 c v around each point p while d is set to 1m.