RECONSTRUCTION OF BRIDGE SUPERSTRUCTURES FROM AIRBORNE LASER SCANNING POINT CLOUDS

The paper deals with the 3D reconstruction of bridges from Airborne Laser Scanning point clouds and cadastral footprints. The generated realistic 3D objects can be used to enhance city models. While other studies have focused on bridge decks to fill gaps in digital elevation models, this paper focuses on the decomposition of superstructures into construction elements such as pylons, cables and arches. For this purpose, the bridge type is classified, and a combination of model-based and data-based methods is used that are built on the detection of arcs, catenaries, and line segments in the point clouds. The described techniques were successfully applied to create 3D models of the Rhine bridges in the German state of North Rhine-Westphalia.


INTRODUCTION
In the last decade, work on automatic computation of 3D city models from Airborne Laser Scanning (ALS) point clouds and cadastral footprints has focused on building roofs. Only a few papers deal with the reconstruction of bridges, see the overview in (Wang et al., 2018). Most papers are concerned with the detection of bridges in point clouds but not with model reconstruction, see e.g. (Sithole , Vosselman, 2003). Another field of research is the use of Terrestrial Laser Scanning (TLS) for bridge inspection, see (Truong-Hong , Laefer, 2014). These methods focus on a detailed reconstruction, e.g., with a mesh to detect deformation, but not on a generalized representation of whole bridges for the use in city models. Recently, a deep learning algorithm was used in (Hu et al., 2021) to reconstruct two large suspension bridges from individually acquired UAV images. Whilst the previous study (Goebbels, 2021) of our research group focused on assembling smooth surfaces of bridge decks from planar polygons, it also dealt with simple but fast reconstruction of building elements that are above the deck level of the bridge. This was achieved with a general data-based approach not restricted to suspension bridges: first, planes were matched to laser scanning points, and then points were projected onto these planes, thereby generating 2D images with sparsely set pixels. These pixels were connected to nearest neighbors by line segments, contours were detected, and extruded in 3D. Although this algorithm makes bridges recognizable, it does not allow to distinguish between different construc- * Corresponding author  tion elements, and there are artifacts due to resolution and vegetation. Therefore, we enhance the proposed technique in the following sections and discuss a mixed data-and model-based approach to reconstruct what we call superstructures. Throughout this paper, we use the term "superstructure" only for large building elements that are above the deck level of the bridge. This excludes the deck and everything below. We also do not consider railings. These elements of the bridge have to be reconstructed with other algorithms, e.g., with the one described in (Goebbels, 2021).

METHODOLOGY
The first question that arises is whether sparse ALS point clouds (see Figure 1) contain enough information to identify structural elements such as cables, pylons, steel truss, and arches. These objects can be described as 3D curves. One algorithm capable of identifying 3D curves in point clouds is TriplClust (Dalitz et al., 2019). This algorithm returns clusters of points that belong to the same curve. The examples in Figure 2 show that large construction parts are recognizable and are correctly segmented. As the runtime of TriplClust is quadratic in the number of points, we ran it on only 20% of the data points. This thinning had the side effect that some structural elements, like vertical pillars, were less well represented in the data and were thus classified as "noise" by TriplClust. Another problem with TriplClust is its large number of parameters. In this case, it turned out, however, that the only critical parameter was the minimum distance at which parallel curves are split up, and we set it to 0.3 m (by specifying TriplClust variables s = 0.1, t = 3.0). TriplClust does not yield curves, but clusters of points, to which a curve might be subsequently fitted. It is thus desirable to have a faster algorithm that directly yields the fitted curves.
In our use case, line and circle segments as well as catenaries are embedded in planes, and the problem can therefore be reduced to 2D curve fitting, which is easier to treat than 3D curve fitting. Thus, we first apply the RANSAC algorithm to find planes with many inlier points, see (Fischler , Bolles, 1981). Within these planes, we then apply RANSAC again to detect planar curve segments. Therefore, we get not only points belonging to a curve, as in TriplClust, but also the type of the curve, from which we can classify the type of bridge. Detecting planes prior to detecting curves has some advantages. Planes can be adjusted to the cadastral footprint to avoid RANSAC inaccuracies, and it is easier to combine different curve segments if a common plane is known, see Section 3. Planes also help to detect barely visible structures such as vertical cables of suspension bridges, see Section 4. We even need adjacent planes to add construction elements between them in a model-based way, see Section 5. Thus, we propose the following algorithm: • Planes are iteratively detected with RANSAC and adjusted according to model knowledge, • 3D points are projected to each plane. Then for each plane, line, circle and catenary segments are fitted with RANSAC to an upper boundary curve of projections. The type of the first (largest) curve segment of the first detected plane (with highest number of inliers) defines the type of the bridge, see Section 3.
• For each plane, the frequency and position of vertical cables or struts is determined using a discrete Fourier transform, see Section 4.
• Based on model knowledge, bridge type specific rules are used to generate a 3D model from upper boundary curves and frequency information, see Section 5.
While most construction parts can typically be covered with few planes, there exist cable-stayed bridges where the cables run in individual directions. Therefore, additionally to the described algorithm, we also experimented with RANSAC on the 3D cloud to detect straight lines, see Section 6.
The method was applied to 28 large bridges spanning the Rhine River in the German state of North Rhine-Westphalia 1 and two bridges in the harbor of Duisburg. The used cadastral footprints and ALS point clouds with a resolution between four and ten points per square meter are freely available at GeoBasis NRW 2 . Information about instruments and scanning settings are not provided.

RANSAC-BASED DETECTION AND CLASSIFICATION
Our algorithm reconstructs superstructures of four main types of bridges: real and self-anchored suspension bridges, tied-arch bridges, cable-stayed bridges, and steel truss bridges (without arches). To this end, it looks for arcs, lines, and catenaries that coincide with many points.
Pylons and cables may lie outside the footprint polygon from the cadastre. Therefore, the algorithm considers points within the dilated footprint of a bridge that lie above deck level. To this end, it estimates the local elevation of the deck at a point (x, y) by calculating the 0.25 quartile of the z-coordinates of all points with x-y-coordinates in a surrounding of (x, y).
Arcs, lines, and catenaries lie within one or more 2D planes. One can directly detect these curve segments by applying RAN-SAC to the subset of the 3D point cloud. However, as reasoned before, it is advantageous to first select planes with maximum numbers of inlier points and then detect the curves within the planes. In general, relevant planes are orientated in the longitudinal direction of a bridge. This direction corresponds with a longest footprint edge from cadastral data. With the exception of tied-arch bridges, relevant planes are also often exactly vertical. By adjusting planes accordingly, they can be estimated with higher precision than by only using somewhat randomly distributed cloud points. Once a plane is selected, points within a threshold distance can be projected orthogonally onto the plane, and a binary image of isolated points showing a layer of the superstructure can be derived (cf. (Goebbels, 2021)). Each pixel represents 0.1 m × 0.1 m. This resolution matches the accuracy of the point clouds that were available to us. Then 2D computer vision techniques can be applied to the binary image. It is also easy to analyze the image from top to bottom and restrict RANSAC to detect curves only in relevant image areas. This is described in the following subsections. Building elements have some width. To avoid detecting objects twice, there has to be a minimal distance d between adjacent vertical planes. We chose d > min{7, w/2} meters, where w is the width of the bridge and seven meters is taken as a small default lane width.
If no plane with a sufficient number of inlier points is found then the bridge is classified to have no significant structures above deck level, and the following steps are not performed.

Upper boundary curve in a projection image
A curve through the highest projected points (upper boundary curve) characterizes the type of bridge. This curve is represented by the graph of a real-valued function.
• For a suspension bridge, this function consists approximately of parts that are a solution (catenary, cf. Figure 3) Figure 3. Three function segments of type 1 k cosh(kx + c) + h fitted to a main cable of the largest German suspension bridge in Emmerich; projected points are black.
of the differential equation • For an arch bridge, the function is defined piecewise by (upper) circle segments y2(x) = r 2 − (x + c) 2 + h.
To identify one of the previously discussed bridge types, only the longest segment of the upper boundary curve on the longitudinal plane with the largest number of points is considered. If, for example, this segment can be described with a function of type y1 then the class "suspension bridge" is chosen. The curve y(x) on the plane is represented by sample points (xj, y(xj)), j ∈ {1, 2, . . . , n}. To obtain the sample points and cope with sparsity, the highest pixels of each consecutive group of ten image columns (range of one meter) is determined. For each of the ten columns, the algorithm checks if its highest pixel/point is closer than 3 m to the highest pixel of the group. If this is the case, the point is added to a list sorted in ascending order by the x-coordinates. Then outlier points are removed that differ from their neighbors in height by more than one meter.
To detect the parts of the function y(x), a RANSAC algorithm is used iteratively that estimates the type y k , k ∈ {1, 2, 3, 4}, and the parameters of the largest remaining part. RANSAC finds a best model by iteratively estimating model parameters of all four model types and counting inliers. These are pixels closer to the graph of y k (x) than a threshold distance (0.3 m).
Within each internal RANSAC iteration, three different indices j l ∈ {1, . . . , n − 10} belonging to remaining sample points (xj l , y(xj l )), l ∈ {1, 2, 3}, are randomly selected. From these points, the model parameters of the different model types are calculated and the inliers are counted. Also, the first and the last inlier points are noted to determine the position of the estimated segment. RANSAC might find models that describe interrupted curves. If interruptions are wider than 10 m, only one connected segment of the curve is counted. The segment is chosen randomly.
Points belonging to a curve segment estimated with RANSAC (due to a highest inlier count) are removed from the list of sample points prior to the detection of the next segment.

Parameters of catenaries
To find the parameters of a catenary y1(x) = 1 k cosh(kx + c) + h, we do not solve the non-linear equations y(xj l ) = 1 k cosh(kxj l + c) + h but obtain the parameters k and c from the derivative y ′ 1 (x) = sinh(kx + c). To this end, we estimate y ′ 1 (xj l ), l ∈ {1, 2}, via difference quotients We calculate the difference quotient with points that are at least one meter apart (index plus ten) to prevent extinctions in the numerator. Then the system of two linear equations (l ∈ {1, 2}) has the unique solution After obtaining k and c, parameter h is determined by If parameter k is negative, then the estimated function is concave and does not represent a cable. It might represent the deck of the bridge. The inlier count is set to zero.

Parameters of circles
If the three points (xj l , y(xj l )) lie on a straight line, i.e.
the inlier count is set to zero. Otherwise, the center (x0, y0) of the circle is calculated with Cramer's rule as the unique solution of a system of linear equations arising from the condition that the three points must have the same distance to the center: Thus, the parameters of y2(x) = r 2 − (x + c) 2 + h are given via c = −x0, h = y0 and Inliers can be counted for arguments x with r 2 − (x + c) 2 ≥ 0.

Parameters of lines
Parameters of y3(x) = mx + c are estimated with m := y(xj 2 ) − y(xj 1 ) c = y(xj 1 ) − mxj 1 . Note that xj 2 and xj 1 indicate different image columns. If |m| < 1 4 , the value of m is replaced by 0 to test for y4(x) = c. If |m| > 3, the inlier count is set to zero.
Straight lines can be approximated very well by circle and catenary segments. We prefer simple straight lines over the other curves. If a bridge is not identified as a suspension or arch bridge, a line is selected as the best model if it has at least 75% of the inlier count of the other models. To avoid confusing lines with arcs, the radius r of a circle is restricted to 0 < r < 10000 meters.

Assembling curve segments
Typically, RANSAC finds multiple curves that have to be connected to a single upper boundary curve. Curves are discarded if their points are completely below the 1/3 quantile of superstructure height values. Then it is very likely that they represent deck structures. The remaining curve segments are sorted by the x-coordinates of their starting positions. Then each curve segment is extended until a function value is approximately equal to the function value of the next curve at the considered position (if possible). This gives the end point of the current curve segment and the start point of the next segment. The start of the first curve and the end of the last curve is adjusted to meet the deck level if this is possible within the footprint.

VERTICAL CABLES OF SUSPENSION BRIDGES, STRUTS OF ARCH BRIDGES
Vertical cables or struts connecting main cables or arches to the deck have a much smaller diameter than other bridge members. Therefore, only a few points belong to these connectors. For example, they were classified as noise by TriplClust, see Figure 2. Some, but not all connectors can be detected with RANSAC or with a Hough transform. An example is the Uerdingen Bridge in Figure 4 that connects the cities of Krefeld and Duisburg. For some reason, the density of the point cloud is much lower on the Duisburg side, so the vertical struts are not visible. The pattern needs to be continued from the Krefeld side. Thus, we propose a different technique. To perform a model-based reconstruction, a histogram function is used, see Figures 4, 5, 6, and 7 for results. To avoid the influence of pylons, each segment of the upper boundary curve is analyzed with a separate histogram function. Let n be the number of image columns in this segment. We define function values g x n 2π of a function g(x), g : [0, 2π) → R, by a weighted number of projected points in columnx ∈ {0, 1, . . . , n − 1}. Only points at least three meters above the deck level and no closer than three meters to the main cable or arch are counted. The distance is needed to exclude points belonging to railings, vehicles on the deck, truss elements of arches, etc. The count is weighted by dividing by the distance between upper boundary curve and deck. The "three meters" restriction is necessary to focus on cable points. Weighting compensates for the fact that cables or struts are of different lengths.
We are looking for a higher frequency j with dominant amplitude 2|g ∧ (j)|. Low frequencies describe the basic shape of the histogram curve, higher frequencies relate to details. Thus, low frequencies j are ignored that imply distances n j between adjacent cables exceeding a threshold of 20 m. A residual frequency j0 with largest amplitude is determined. Then vertical connectors should be placed according to the maxima of hj 0 (x). Since −ψj 0 /j0 defines a horizontal shift to the right, the vertical lines are shifted by −nψ j 0 2πj 0 image columns. Typically, vertical struts or cables are adjusted to the center of the arch or catenary, i.e., they have to be symmetrical about x = −c for an arch or x = − c k for a catenary. Either a strut or cable is placed exactly at this center position, or the center position is in the middle between two of them. The algorithm selects one of the two options that best matches the calculated displacement.
It turned out that the described method is not very stable for sparse point clouds with less than five points per square meter. To some extent, it can be stabilized by drawing line segments (with a Hough transform) into the point cloud image before computing the histogram. Figure 5 shows that distances between cables can be estimated slightly differently for different planes. Therefore, the frequency computation is performed only for all curve segments of the plane with the largest number of inliers. Then the smallest of all computed distances is chosen. This distance and the corresponding type of alignment to the center position is applied to all curve segments on all detected planes. Unfortunately, there are bridges that connect two cities so that their footprint is divided in two. In such a case, different calculated frequencies may still occur if the visibility of the cables varies. If the distance between cables is not within a threshold interval, default values are applied instead.

Diagonal cables of cable-stayed bridges
If the analysis of the upper boundary curve in a projection image indicates a cable-stayed bridge, the cables are detected iteratively with RANSAC for diagonal straight lines. In each iteration, points of the current upper boundary curve y(x) are removed. Then a new upper boundary curve is computed as above (Section 3.1), and RANSAC detects its line segments, which are then combined, see Figures 8,9,10,13,and 24. Lines are extended until they meet each other or until they reach deck level within a threshold distance.

Steel truss bridges
A dominant horizontal line, that is at least tree meters above deck level, characterizes steel truss bridges. We model such bridges generically with post-less truss, see Karl Lehr Bridges in Figure 14. Often, truss elements run diagonally so that their number can be calculated by dividing the length of a bridge segment by its height above the deck. Additionally, the frequency of truss elements is determined by applying a Discrete Fourier Transform to a histogram function, similar to Section 4. In this case, only points in the lower half of the superstructure are counted. Otherwise, the histogram function would be almost constant. If the DFT-based frequency is significantly (a) One half of Theodor Heuss Bridge, cf. Figure 1 (b) Oberkassel Bridge Figure 9. Cable-stayed bridges in Düsseldorf.
(a) Points projected to the first RANSAC plane (b) Model-based reconstruction Figure 10. Düsseldorf Airport Bridge: triangular pylons do not match with our model assumptions. Figure 11. The cables of the Severin Bridge in Cologne were estimated directly on the 3D point cloud using RANSAC. The pylon was reconstructed data-based, see Section 6.2.
higher than the frequency obtained from length and height then additional truss elements are added, see Figure 16.

Pylons
Typically, the pylons are located at intersection points of the model curves when the points represent local maxima. We add them model-based by estimating their width from the point cloud. The heights are taken from the upper boundary curves. By averaging different values, all pylons are extended to the same height, see Figure 21. In case of suspension bridges, adjacent pylons belonging to different planes are connected below their tops.

Model-based enrichment of arches
We use model knowledge to show tied-arch bridges in a beautiful way. If, unlike the Bridge of Solidarity in Figure 22, arches do not reach the deck level at their end points within a threshold distance of three meters, a vertical end line is added. In this case, a second arch is also added below each detected one. The additional arch connects the end points on deck level. Diagonal connections are added between the two arches and between consecutive vertical strut elements. Also generic connections between adjacent planes are added, see Figures 14, 17, 22, and 23. At least, this fits with a certain type of tied-arch bridges. Examples can be found in Cologne and Duisburg. Some of the Elbe bridges in Hamburg are also of this type. However, the railway bridge between Moers and Duisburg in Figure 17 has a different truss pattern in reality. In contrast to Hohenzollern Bridge, the arches of South Bridge in Figure 23 continue below the deck. But structures below deck level are typically occluded in ALS point clouds.

3D model generation
When curves and vertical cable lines are detected in RANSAC planes, they are drawn into an image from which contours are detected and extruded in 3D. However, depending on the bridge type classification, some curve segments may be discarded. We avoid mixtures of different bridge types.
Truss elements between planes, model based pylons and cables detected with 3D-RANSAC (see next section) are modeled as cuboids.
All 3D objects are encoded in CityGML as separate bridge construction elements, so that semantics can be assigned to them, see Figure 15. CityGML is the standard format for 3D city models, see (Gröger et al., 2012), and (Kutzner et al., 2020). Figure 26 shows an integration of some computed bridges into a CityGML-based model of Cologne.

Plane-independent reconstruction of cables
In addition to the general algorithm proposed in this paper, we also tested an alternative approach for bridges classified as cable-stayed bridges. Some cable-stayed bridges may have cables that do not lie within common planes. Although this was not the case for the bridges that we reconstructed, we tested with RANSAC to find cables of the Severin Bridge, the Rees Bridge, and the Wesel Bridge (see Figures 11,12,20) within the 3D cloud instead of planes. In this case, we computed all intersection points between cables. The highest intersection points   (that still are close to some inliers and have sufficient mutual distance) define the probable positions of pylons. These positions can be used to clip cables if necessary. Pylons may be mistaken as cables. To avoid this, all cables with lower end points near pylons are discarded. This also applies to cables with upper end points that are too close to the deck.

Data-based pylons
The method of Section 3 can be used to find the outer contour of pylons in a data-based fashion. Examples are shown in Figures 11 and 20. To this end, only RANSAC planes perpendicular to the longest footprint direction are considered. The upper boundary curve of projected points is constructed from line segments. In contrast to Section 3, lines are also recognized if their slope |m| exceeds 3. Then, the curve is expanded into a 3D model of the pylon.

RESULTS
As described in Sections 3-5, our algorithm assembles most bridges from curves detected on planes, see Table 1. Due to the nature of the algorithm, bundles of cables are modeled as one cable. This is not counted as an error.
The bridge type was correctly recognized for 27 out of a total of 30 bridges. Reasons for wrong detections were power lines and a combination of steel truss and arch bridge. Although power lines should look like catenaries, bridges without significant superstructures but with light rail tracks or electrified railroad lines can be mistaken for steel truss bridges, see Figure 19.
Whereas the detected arches, catenaries and line segments of upper boundary curves matched the real structures within the threshold resolution of 30 cm, see Figure 25, we faced some other problems: The frequency and shift of the vertical construction elements matched only approximately, see Figure 25, where green struts and cables can be compared with black line segments drawn on projected points based on a Hough transform. The more vertical lines visible in a projection image, the better the Fourier analysis results. Figure 18. Ebert Bridge in Duisburg is composed from two cadastral bridges that are drawn with different colors. One main pillar was not available from data.  (Figure 9) recognized. are missing.  Pillars and counter bearings were often missing or incomplete in cadastral data. Three bridges had incomplete cadastral footprints, and many footprints were divided into two or more separate bridges. This had an impact on the quality of the recognition. For example, the Rheinknie Bridge in Düsseldorf (see Figure 24) is composed from several parts so one does not find a medial axis that connects counter bearings. Yet, this is a prerequisite for high-quality tesselation of the bridge deck. The decomposition also separates the cables into smaller pieces that become too short to be detected. Another problem occurs when two cadastral footprints exactly meet at a pylon. Then the pylon is not recognized, see the Ebert Bridge in Figure 18.
Cable-stayed bridges with many cables were also a challenge due to the sparsity of the available point clouds. For example, not all cables of one pylon of the Theodor Heuss Bridge in Figure 9 could be found, and the second pylon with its cables is completely missing. The model of the Oberkassel Bridge shows that longer cables are easier to detect than short ones.
Most bridges do not have characteristic structures above deck level. Therefore, the small number of large suspension, arch, cable-stayed and steel truss bridges with scenic superstructures allows for manual interaction. The results were corrected by  choosing between the following options: • Disabling superstructure generation so that power lines or trees cannot be confused with truss elements • Adjustment of RANSAC thresholds • Plane-independent cable detection with RANSAC on the 3D point cloud, with or without clipping cables at estimated pylons (Section 6.1) • Data-based reconstruction of pylons based on projection images instead of model based pylon reconstruction (Section 6.2).

CONCLUSION AND FUTURE WORK
Even if manual corrections are necessary in some cases, the superstructures generated by the presented methods are much more realistic than the fully automatically computed structures in (Goebbels, 2021) on which we have built. The methods thus can fill the gap in most published city models, which lack the presence of bridges. In future work, it might be interesting to set up a database of rules and attributes similar to 3D building reconstruction, cf. (Wonka et al., 2003), to improve models based on grammars. Instead of sparse ALS clouds, it would be interesting to consider dense TLS clouds as they allow for a higher level of detail. We focused on bridges, but similar procedures can be applied to similar structures, e.g., power lines.