# Code is poetry
Critical infrastructure such as bridges and tunnels are surveyed on regular basis to detect fractures, deformations and other structural weaknesses. While this work is currently performed manually by surveyors, there is an opportunity to use 3D twins, LiDAR data and advanced machine learning models to improve a manual process which has remained little changed in past 100 years.
At a recent client project, our team was presented with the challenge to analyse point cloud data of UK railway tunnels, and to design, build and assess mathematical models which can automated the detection and categorisation of safety defects in rail tunnels.
The existing inspection process requires surveyors to physically walk through a tunnel and annotate defects on a tablet, necessitating the closure of tunnels for up to 12 hours, while automating data capture can potentially reduce that timeframe to 1-2 hours, unlocking 5x to 10x time savings.
The current blog article presents one part of our research, specifically related to the Alignment of Tunnel Point Cloud Data using Central Line Extraction and Region-Based Rotation.
Our team started by exploring large datasets, in which a single datapoint, scan of an individual tunnel, can be more than 30GB. We attempted produce a general enough algorithm that will work with most tunnel profiles. The geometry of tunnels can vary significantly in size, shape and texture. Some are "more circular" than others, some are more "curvy", yet others have have unique tunnel features such as open shafts or tunnel refuges. The construction material can also vary significantly – from stone, to masonry to concrete.
Before any defect detection algorithm can be run on tunnel point cloud data, a number of sequential preprocessing steps need to take place such as data cleaning (removal of unnecessary artefacts) and alignment. We created these preprocessing steps with the goal of taking any 3D tunnel scan and outputting a straight tube-like tunnel running parallel to, say, the x-axis. This lets the defect detection models generalise well by assuming that all tunnels start at 0 on the x-axis and run along it.
Essentially, we want to ‘straighten’ the curvy tunnel and align all cross sections. In collaboration with an MSc student from the School of Mathematics at the University of Edinburgh and two academics in the span a year we have been developing a new algorithm which transforms a point cloud scan (LiDAR 3D data) of any tube-like data (pipe, tunnel), with any profile shape and curvature, to a form where all data cross-sections are aligned.
The proposed algorithm transforms the tunnel by reducing this observation error between cross sections, which allows us to study more easily any tunnel lining deformation and degradation. Few previous algorithms have been developed to reduce the misalignment of cross-sections and the previous focus has focussed on tunnels with a circular profile, the proposed algorithm therefore usefully extends the range of tunnels which can be considered.
This article goes over the basics of the proposed algorithm using a noiseless toy example of a tunnel to show how well the model works and it will quantify the error.
For the purpose of this article we created a synthetic tunnel data with no noise to be able to exactly quantify the accuracy of the algorithm. The idea is loosely based on the following StackOverflow post.
The StackOverflow article goes in enough detail, so I won’t focus on explaining the logic behind the code. We want to create a cylinder shape which will follow a specific predefined curve. We will work with sin(\frac{x}{16}) to make the data look like a more realistic tunnel which won’t curve as often. We define a few helper functions.
d = 1/16
def f_(x):
return np.sin(x*d)
def normalize(vec):
return (lambda norm: [x/norm for x in vec])(sum(x**2 for x in vec) **0.5)
def w_(x):
return normalize([1, d*np.cos(x), 0])
def u_(x):
return normalize(np.cross(w_(x), [0, 0, 1]))
def v_(x):
return np.cross(w_(x), u_(x))
Then we find our XYZ’s using the parametrised 3D surface.
def surface(Q, T):
x = [q + np.cos(t)*np.array(u_(q)[0]) + np.sin(t)*v_(q)[0] for q, t in zip(Q, T)]
y = [f_(q) + np.cos(t)*np.array(u_(q)[1]) + np.sin(t)*v_(q)[1] for q, t in zip(Q, T)]
z = [np.cos(t)*np.array(u_(q)[2]) + np.sin(t)*v_(q)[2] for q, t in zip(Q, T)]
return x, y, z
ux, vx = np.meshgrid(np.linspace(0, (1/d)*2*np.pi, 500),
np.linspace(0, 2*np.pi, 500))
x, y, z = surface(ux.reshape(-1, ), vx.reshape(-1, ))
Finally our tube data looks like this (visuals made with Open3D):
Given the PCD, the first part of the algorithm finds a set of points (each with a given step-size ‘h’ apart) within the tunnel that aim to describe the tunnel's centre line.
Points along the tunnel are then selected along the central line by solving a series of optimization problems. For some choice of constant integer K , we proceed by defining an objective function as follows:
f(x) = median (\textbf{d}_{i}), \; \textrm{where} \\ \;\; \{ \; \textbf{d}_{i} = \lVert x - x_{i} \rVert _{2}, \; \textrm{and} \\ \;\;\;\;\; \textbf{x}_{i} \; \textrm{are the K nearest points in the point cloud to location \textbf{x}} \; \}
We intend to maximise this function, finding the point maximally away from the point cloud. To ensure we continue to find points on the centre line we add the following additional constraints:
h^2 - \lVert x - o \rVert^2_2 = 0, \\ n(x - o) \ge 0
Here \textbf{o} is the previously found centre point, and \textbf{n} is the vector between the previous two determined centre points (at the initial step these are the manually chosen down-tunnel direction). These two conditions ensure, respectively, that the next proposed centre point is exactly on a distance \textbf{h} from the last, and that it is further down the tunnel.
We then fit a cubic spline to the synthetic data with smoothness parameter s = 0.9 and split the tunnel into 100 cross sections when fitting using piecewise cubic splines. Again, to better understand how well this line fits the expected centre we can examine the projections of the line to the x-y, and x-z planes to compare the calculated centre line and the true centre line.
For the synthetic data we calculated the RMSE (Root Mean Square Error) of the determined points normalised by the length of the curve which was found to be 20.05 (length of the sine wave from 0 to 100.
In the \textbf{y} direction this is 0.00015 and in the \textbf{z} direction 0 (however our synthetic data is straight in the you example).
As a next step, we want to transform the PCD to a new axis, where, the \textbf{x}-axis will represent travelling down the tunnel, the \textbf{y}-axis will represent travelling across the tunnel and the \textbf{z}-axis will represent travelling vertically within the tunnel. Additionally, the transformation will result in the centre line of the tunnel laying on the \textbf{x}-axis.
To achieve this, the vectors between points sampled from the spline will be calculated and then used to define regions along the tunnel. Each region can be defined as between two planes, each sharing a normal equal to the vector joining the sample points mentioned before, and containing the sample point. A 2D visualisation of such a region can be seen in the figure below, where the vertical lines represent the boundary planes and the horizontal line represents the region's associated vector. After the algorithm assigns each datapoint to a region, each region is transformed such that the beginning point of the region's associated vector is mapped to the origin. Then, a further transformation is applied to the region's points that rotate all the points such that the region's associated vector lies on the \textbf{x}-axis.
With the noiseless data, this step is not really necessary, so if you are working with data which is quite smooth the solution might be good enough already. However, locally, cross sections can still be misaligned especially due to sensor errors while gathering data. To overcome this, for a given section, its data points can be projected onto the y-z plane where a smooth function may be constructed using Kernel Density Estimation (KDE). After the smooth function is built for each cross-section, neighbouring sections can be aligned by simply scoring the points of one section using the function from their neighbour. The required translation to align the \textbf{y} and \textbf{z} components of the data points can be found using an optimisation algorithm to maximise this score. This alignment can then be repeated for all the defined cross-sections of the tunnel.
The left image shows points from a cross-section of tunnel in blue that is centred on the origin point \textbf{(0,0)}, and points from a cross-section which has been adjusted such that it is centred on the point \textbf{(1,1)} shown in orange, with both cross-sections having 100 data points. The output of the proposed alignment algorithm is then shown in the right image of the figure.