RAIL TRACK DETECTION AND MODELLING IN MOBILE LASER SCANNER DATA

We present a method for detecting and modelling rails in mobile laser scanner data. The detection is based on the properties of the rail tracks and contact wires such as relative height, linearity and relative position with respect to other objects. Points classified as rail track are used in a 3D modelling algorithm. The modelling is done by first fitting a parametric model of a rail piece to the points along each track, and estimating the position and orientation parameters of each piece model. For each position and orientation parameter a smooth low-order Fourier curve is interpolated. Using all interpolated parameters a mesh model of the rail is reconstructed. The method is explained using two areas from a dataset acquired by a LYNX mobile mapping system in a mountainous area. Residuals between railway laser points and 3D models are in the range of 2 cm. It is concluded that a curve fitting algorithm is essential to reliably and accurately model the rail tracks by using the knowledge that railways are following a continuous and smooth path.


INTRODUCTION
Rail track irregularities have a large effect on railway safety and operation.To ensure a good maintenance of the rails, frequent measurements are needed, which are costly and require specific tools for different aspects of the rails geometry.Mobile laser scanning (MLS) offers the advantage of acquiring accurate 3D measurements of all the objects present in the railway environment in a short operational time.Another advantage of using mobile laser scanner data mounted on a train vehicle is that there is no need for surveyors to enter the rail track.
Producing 3D models of the rail environment is useful for many applications such as asset inventories, analysing the minimum free passage space, and determination of measurements such as platform position in relation to the rail track.Manual detection and modelling of the rail tracks in a point cloud is largely impractical for two main reasons.Firstly, recognizing and precisely delineating the tracks in sparse points is difficult for a human user.Secondly, the extensive length of the rails makes their detection and modelling a very tedious task.
The idea is to apply a knowledge based classification which takes advantage of the regularity in a railway environment to classify points on the objects of interest.The focus in this paper is on detection of points on railway tracks, followed by a 3D modelling step.Local properties are piecewise linearity of rail tracks.Two parallel tracks form a pair at a certain fixed distance, i.e. the gauge.Globally railways follow a continuous and smooth pattern.Our contribution is the integration of local and global geometric properties of the railway during both the detection and modelling steps.After describing the related work (section 2) our detection algorithm is explained in section 3. Section 4 handles the modelling steps, followed by an explanation and analyses of our results in section 5. Conclusions and future work are presented in section 6.

RELATED WORK
Classifying laser data has been a research topic for several years.In airborne laser scanner (ALS) data the challenges are on analyzing the influences of training sizes, feature selection (Pal and Foody, 2010) and the classifier itself (Pal and Foody, 2012).Often the aim of classifying ALS data is to find specific objects, such as buildings and vegetation (Xu et al., 2012), or to filter non ground points (Tovari and Pfeifer, 2005).Automation in finding the rail track centre lines using high resolution aerial imagery and lidar can be found in Beger et al., (2011).Jeon (2010) detects and models catenary wires from airborne laser scanner systems.Rutzinger et al (2011) describe the feasibility of building footprint extraction from MLS data.Several studies show that MLS data can be used for asset inventory of the rail side hardware and engineering design work, and to extract highly accurate spatial information for construction applications and maintenance (Leslar et al, 2010).Although these studies show the potential of using MLS data for applications, the measurements themselves are still based on manual interpretations.Chan and Lichti (2011) explain how to fit catenary curves to power cables, in order to calibrate a mobile mapping system.Arastounia (2012) describes two methods for detecting rail points, based on a template-based matching and a region growing approach.The template matching emphasizes the detection of points on parallel tracks, whereas the region growing approach uses the knowledge on the continuous shape of the tracks.Fitting parametric models to points has been described in many reverse engineering projects, for example curve fitting (Werman and Keren, 2001), for modelling 3D buildings (Verma et al, 2006), (Maas and Vosselman, 1999) and industrial installations (Rabbani and Van den Heuvel, 2004).
To our knowledge there is no approach that detects and models rail tracks by fitting parametric models, using both local and global properties of the rail way environment.

KNOWLEDGE BASED DETECTION OF RAIL TRACKS
The aim of this step is to detect and select laser points on a railway track.These points will be used to fit 3D models of the railway.As the rail modelling stage (described in section 4) includes a curve fitting stage, it is considered to be less harmful to miss a few railway points (false negatives), than to include too many false railway points (false positives).Our detection algorithm is based on object specific properties of the railway environment, which are listed in  1. Properties of points on rail tracks and contact wires.
These properties are implemented in a four step approach.In the first step the height distribution of all laser points are analyzed per 1 by 1 m grid cell.The reason for this grid size is that the expectation is that there is maximum one rail track within a grid cell.In figure 1, points can be seen from 3 neighboring grid cells, in a horizontal viewing direction.The purpose of this initialization step is to roughly indicate whether in this grid cell there are points on rails and/or wires.Basic assumption for the rough detection of rail points is that the grid cell contains points on the terrain, there are several points slightly above the terrain (potentially the points on the railway track), there are almost no points between the rail track and the wires, and there may be some points on the wires at a certain height above the terrain.Starting point is the determination of DTM height per grid cell.As an initial guess the 10%-ile height of all points within the grid cell.If there are less than 10 % of the points within 0.5 and 4.5 m above DTM height, there may be rail and or wire points.For roughly detecting points on wires the assumption is that points on wires are between 5.5 and 6.5 m above terrain level, only the lowest 5 cm of those points are taken as potential contact wire points.All points within 0.5 meter above DTM height (called "the terrain points") are further analyzed for rail point detection, see figure 2. If the difference between the 98%ile height of the terrain points and the 10%-ile point is larger than 10 cm, there may be a rail track inside the grid cell.All points within 10 cm of the 98%-ile point are potentially rail track points, but only if this is not the majority within the grid cell, see figure 2. These criteria come from the knowledge that in most situations the railway is about 10-15 cm elevated above the terrain.For the exceptions, e.g. at road crossings, our modeling strategy (see section 4) is designed to bridge the gaps in the detected railway points.The second step is to keep only points that represent linear structures.Within a grid cell of 1 by 1 m there is maximum one piece of rail track and/or 1 piece of wire.For the line fitting a RANSAC algorithm (Fischler and Bolles, 1981) is used, where the assumption is that the majority of the roughly detected rail track points within a grid cell, actually fits to one line within a certain buffer (say 0.05 m).Two parameters are used here: the percentage of inliers and the maximum distance of inliers to the line.Results of a successful line fitting are the inlier points which represent the detected rail tracks and wires, plus a 3D fitted line through the inliers.When closely looking at this step, it is obvious that the detection algorithm is designed for detecting points at the top part of the rail track.The RANSAC line fitting selects points within a buffer of 0.05 m on points that are slightly elevated above the terrain.For reliably fitting a complete rail track model, it is desired to also include the points at the foot of the rail track.Only after a successful inlier detection of the rail points, other laser points are added if there are within 0.12m distance, i.e. the height of the body of the track, to the line.For rail tracks it always implies the addition of the points directly beneath the top of the rail track.Until this point, the processing has been done grid cell wise, analyzing all points within a square meter.The result is a rough classification of the point cloud based on histogram analysis in combination with linearity restrictions.From now, the processing will be done on the roughly classified point cloud.
The third step is optional and contains an extra filtering step by keeping rail points only if there are wire points within a certain 2D and 3D distance, and vice versa.The assumption is that there is a wire somewhere near the rail track, for example within 2 m in the horizontal plane and between 5 and 6 meter in vertical direction.As a result one can immediately assign a wire ID to the rail points on a pair of rails, making it possible to determine the number of tracks in a certain area.Assumption is that there is one contact wire for each pair of rails.
The fourth step is to check whether two rail tracks are parallel and at a certain constant distance from each other.For each RANSAC line it is checked whether there is another parallel line at a certain distance (the gauge in this dataset is about 1.45 ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey m).This is done by projecting mid points of the lines to other parallel lines and see whether the distance is between 1.3 and 1.6 m.If this is the case, it is a strong indication that this pair of lines, and thus two pairs of inlier points, are actually representing rail track pieces.If a corresponding line cannot be found, there can be two reasons: firstly the processed line can represent some false positives, or secondly, at the location of the corresponding rail track there is no line detected by the RANSAC algorithm.So, at that location the roughly classified points are checked whether they fit to the railway track, but for some reasons were not yet classified as such or not considered to fit to a RANSAC line.These potential points are added to the detected rail track points.A connected component analyses is performed to group points on each railway track.Only points at large components are finally classified as rail track points.The main purpose of this fourth step is to decrease both the number of false negatives (by adding potential points) and the number of false positives (by removing the small components).After finding correspondences between two tracks, it is determined for each pair whether points belong to the left or to the right track.At this stage it is known per laser point whether it belongs to a rail track and to which of the two rail track within a pair.This will be the input for the railway modeling stage.

Parametric model of a rail piece
A rail track is modelled as a set of smaller rail pieces.A rail piece is defined by seven shape parameters, as shown in Figure 4, and six orientation parameters that specify the position and rotations of the local coordinate system of the piece with respect to the global coordinate system of the point cloud.To fit the rail pieces, first the points detected on each track are partitioned into straight segments using a planimetric grid of 50 m cell size.For the points within each cell the eigen vector corresponding to the largest eigen value is calculated as the main axis of the segment.The points are then again subdivided into segments of size L (piece length) along the main axis.
In the fitting of the piece model to the segmented points, the shape parameters are considered fixed and only the orientation parameters of the model are estimated.The least-squares method for the estimation of the orientation parameters often fails because of the sparsity of the points.We therefore use a Markov Chain Monte Carlo (MCMC) algorithm to obtain an estimate by sampling the joint probability distribution of the orientation parameters.Formally, for a given point segment D the aim is to find a model M i that maximizes the probability P(M i |D) = ηP(D|M i )P(M i ), where P(D|M i ) is a measure of how well the model fits the data points D, P(M i ) is the model prior, and η is a normalization factor that is independent of M i .We define: where d is the mean distance from the points to the model M i .The point-model distance is defined as the smallest distance between the point and each of the planar patches of the model (see Fig. 4).However, distance calculation based on the mathematical equation of a plane may lead to small distances for points that are outside the patch boundary but lie on the extension of the patch plane.To exclude such incorrect distances we introduce the following additional condition: for a point-plane distance to be accepted the point should orthogonally project inside the polygon that encloses the patch.Distances not fulfilling this condition are excluded from the calculation of d in Eq. ( 1).The prior P(M i ) is used to incorporate our prior knowledge of the model parameters.We expect with a high probability that the rotation parameters of each rail piece are only slightly different from the previous piece, and that the position of a piece is close to the center of its corresponding point segment.To include these we model the prior with a normal distribution centered around the expected rotation and position parameters of the piece model: where u = [x o , y o , z o , v, f, k] T is the vector of six orientation parameters of model M i , µ is the mean vector and ∑ is the covariance of the parameters in u.The mean vector contains the expected rotation (i.e. the rotation of the previous piece) and position parameters (i.e. the center of the point segment).For the first piece in each track we set v=0 and f=0, while k is obtained as the orientation of the main axis of the point segment (eigen vector corresponding to the largest eigen value of the segment).The covariance of the prior distribution ∑ is chosen by assuming large variances for the expected orientation parameters and no correlation between them.The large variances ensure that the Markov chain does not get stuck at the mean of the prior distribution.
For the MCMC sampling we use the Metropolis-Hastings algorithm (Hastings, 1970), which starts at a random initial point and recursively provides samples µ t , t=1,…,n from the target distribution P(M i |D).The samples are drawn from a proposal distribution, which is chosen as a Gaussian with variances larger than the prior distribution.An estimate of the expectation of model parameters is then obtained by the ergodic mean of the samples: θ ∑ θ , where m specifies the number of so-called burn-in samples.Further details about the Metropolis-Hastings algorithm can be found in (Gilks et al., 1995).

Curve fitting
The rail pieces modelled by the MCMC sampling do not necessarily form a continuous and smooth rail track.A typical rail track is a combination of linear segments and smooth circular curves.Therefore, to obtain a continuous and smooth model of the entire rail track we interpolate the rail pieces with a Fourier series: where u j ,j=1,..,6 is jth orientation parameter, i=1,..,m is the rail piece number, a k , b k , w are unknown coefficients of the interpolation function and n is its order.By evaluating Eq. ( 3) with u j of all pieces a system of equations is obtained, which is then solved for the unknown coefficients.The estimated coefficients minimize the sum of squared differences between the piece parameter and the interpolated parameter.The interpolation coefficients are estimated for each parameter separately, and can be used to evaluate the parameter at any point along the rail track.The order of the Fourier series defines the flexibility in the fitted curve.Note that we want to avoid higher order Fourier series to prevent the occurrence of oscillation and waviness in the final model.When analysing the geometry of the rail pieces, the rotation parameter around the Y-axis of the rail model (see Figure 4), i.e. the φ angle, is the most sensitive to the random error and scarcity of the laser points.This can also be seen in the evaluation of the φ parameter.Therefore, it has been decided to determine φ by calculating the slope between two parallel rail pieces.This can easily be done as the position of the rail pieces is very well determined.This newly determined φ is assigned to both parallel rail pieces.In Díaz Benito (2012) a Bezier curve fitting approach has been described.This works well for almost straight rails giving continuity and smoothness to the piece model, however Bezier curves are tangent only to the initial and end sections, and not bounded at all to the local parameters of intermediate pieces.Therefore, such a model behaves poorly in intermediate sections of curved tracks, as it cannot closely follow the railpoints.

RESULTS AND ANALYSIS
The proposed methods were evaluated by processing two point clouds obtained by a mobile laser scanner on the Austrian railway.

Mobile laser scanner dataset used
The data has been acquired by Topscan, using an Optech Lynx V1 mobile mapping system.The system, containing two scanners, was mounted on a car which was placed on a train waggon.Point densities near the rail track are about 700 p/m 2 , containing points from both scanners.The target object is a curvy rail track in a mountainous area.Two areas are selected from this dataset.Area 1 is a track of 200 m, containing one curve.The rail track from area 2 follows a S-curve with a total length of 400 m.
Figure 5. Oblique view on the point cloud from Area 1 (left, 2.2M points) and Area 2 (right, 4.2M points).

Detection of points on rail tracks
In figure 6

Curve fitting
Figure 8 highlights the working of the curve fitting to individual fitted rail pieces.The rail track now forms a continuous and smooth curve.We have analysed the RMS residual between piece parameters and interpolated parameters for different number of terms in the Fourier series.It was found that for area 1 three Fourier terms, for area 2 five Fourier terms, for the position parameters, and one Fourier term for rotation parameters in both areas are a suitable choice as we do not see a substantial improvement of the residuals with increasing number of terms.Nearly 5% of the individual pieces did not have enough laser points to accurately fit a single piece, however the other 95% was more than enough to bridge the gaps after the curve fitting.A triangular mesh is built from the curve fitting result to conveniently visualize the final model.Figure 9 shows the curve fitting to each orientation parameter of the pieces.Each dark circle represents an orientation parameter of a piece and each red curve represents the interpolation function (Fourier series) for that parameter.

Final 3D rail model results
The meshes of the rail models are shown in figure 10 and 11.
When overlaying the models to the point cloud, one can globally see that the models fits nicely, although for accuracy analyses one needs to rely on quantitative measures.To analyze how well the models fit to the laser data, point-tomodel distances are calculated and shown in figure 12 and table 2. The two figures below show the distances between the points and the models visualized by color.It can be seen that the residual distances are more or less evenly distributed along the rails, except for a few outlier points (3.7-5.1%), which are seen in dark red.While the influence of outliers is evident from the large mean values, the median distances provide a reliable measure of the accuracy of the final models.It can be concluded from the median distances that the accuracy of modeling the railway tracks is below 2 cm.

CONCLUSION AND FUTURE WORK
We have described a method that can detect and model railways in an automated way.Combining three types of properties of the railway track, relative height, linearity and relative position to other objects, resulted in a highly accurate detection of rail points.The accuracy of the end results is in the order of 2 cm, which is acceptable for many applications that deal with the rail environment, such as asset inventories and visualisations.However, for highly detailed measurements with mm precision, e.g.calculation of wear of rail tracks, it is recommended to have more accurate measurements.
The piece wise determination of rail models can only be used if it is followed by a fitting algorithm that analyses the more global shape of the rail track.The reason is that for small pieces there are too few points to accurately determine all model parameters.Using a curve fitting algorithm is essential to fully use the knowledge that rail tracks are following a smooth path.This can visually be seen, and is grounded by the correction of the parameters without increase of the residuals between points and model.After back-projecting the 3D model of the rail tracks to the original point cloud it is possible to better analyse the detection strategy, and even calculate the number of false positives.This is useful to know the sensitivity of some of the parameters, such as the determination of DTM height per grid cell and setting a certain minimum number of inliers for RANSAC line fitting.Future work will focus on optimizing the detection strategy.Also, the influence of the length of individual rail pieces on the accuracy of the model parameters and the use of clothoids for transitions between straight lines and curves will be analysed further.
The algorithms have shown their performances on relative simple rail environments.In practice, switches, bridges and train stations may ask for other parameter settings or even a modified workflow.Future work is to detect points near special objects such as switches, followed by a parametric model fitting.This is feasible as we can accurately determine where two tracks would meet, and the variety of switch types is limited, at least within one country.

Figure 1 .
Figure 1.Starting point is determination of DTM height, followed by checking empty space between 0.5 and 4.5 m above DTM height.Distance between two ticks at the axis is one meter.

Figure 2 .
Figure 2. Potential rail track points are detected by histogram analyses of terrain points.

Figure 3
Figure 3 Detection steps for rail tracks: rough classification of terrain points after step 1 (upper left, green are potential rail track points, cyan are points lower than the DTM height, yellow are other points within a grid cell with rail points, orange are points in a grid cell without rails or wires).Middle left: result after RANSAC line fitting (step 2).Lower left: including points (yellow) within 0.12 m of fitted line.Upper right: filtering step by keeping only rail track points (green) within a certain distance of wire points (orange) (step 3).Lower right: red lines indicate which fitted lines are parallel with a perpendicular distance between 1.4 and 1.6 m (step 4).Red points are on one side of the track, green at the other side.

Figure 4
Figure 4 Rail piece coordinate system and model parameters 4.2 Parameter estimation by Markov Chain Monte Carlo the (intermediate) results of the detection algorithm are shown.It can be seen that the result of the histogram analyses is giving a rough idea on the potential locations of railways (b) and wires (c).However, information on linearity and relative distances to other objects are needed to remove many false positives.The detected wire points are used to keep only nearby rail track points and thus to filter points that were falsely classified as rail track points.The final detection results are shown in figure 6f.As there was no reference data and the selection of test data would be a tedious work, we have estimated the number of false positives and false negatives by analysing the data gaps and outliers during the modelling steps, see section 5.4 and 5.5.

Figure 6 .
Figure 6.Detection of rail points per step.Input point cloud (a), rough detection of rail points (b) and contact wires (c), inliers of RANSAC fitting (d), large segments of inliers on wires (e), final rail points (f) after checking parallelism and checking nearby wire points.

Figure 8 .
Figure 8. Individual rail pieces before (left) and after (right) Fourier curve fitting for area 1.

Figure 9 .
Figure 9. Curve fitting parameters from top to bottom X, Y, Z (left) and omega, phi and kappa (right).

Figure 11 .
Figure 11.3D models overlaid on points from area 2.

Figure 12 .
Figure 12.Distances between final model and railway points for dataset 1 (top) and 2 (below).

Table 1 .
Although in this paper the focus is particularly on rail tracks, the presence of contact

Table 2 .
Statistics on point-to-model distances for both areas.