# Automatic Detection and Modeling of Underground Pipes Using a Portable 3D LiDAR System

Ahmad Aijazi, Laurent Malaterre, Laurent Trassoudaine, Thierry Chateau, Paul Checchin

► To cite this version:

Ahmad Aijazi, Laurent Malaterre, Laurent Trassoudaine, Thierry Chateau, Paul Checchin. Automatic Detection and Modeling of Underground Pipes Using a Portable 3D LiDAR System. Sensors, 2019, 19 (24), pp.5345. 10.3390/s19245345. [https://hal.science/hal-03089726](https://hal.science/hal-03089726)

HAL Id: hal-03089726

https://hal.science/hal-03089726

Submitted on 8 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Automatic Detection and Modeling of Underground Pipes Using a Portable 3D LiDAR System

Ahmad K. Aijazi *, Laurent Malaterre, Laurent Trassoudaine * and Thierry Chateau
and Paul Checchin *

Institut Pascal, UMR 6602, Université Clermont Auvergne, CNRS, SIGMA Clermont,
F-63000 Clermont-Ferrand, France; laurent.malaterre@uca.fr (L.M.); thierry.chateau@uca.fr (T.C.)

* Correspondence: kamalaijazi@gmail.com (A.K.A.); laurent.trassoudaine@uca.fr (L.T.); paul.checchin@uca.fr (P.C.)

Received: 8 November 2019; Accepted: 1 December 2019; Published: 4 December 2019

<figure><img src="image_3.png" /></figure>

check for
updates

**Abstract:** Automatic and accurate mapping and modeling of underground infrastructure has become indispensable for several important tasks ranging from urban planning and construction to safety and hazard mitigation. However, this offers several technical and operational challenges. The aim of this work is to develop a portable automated mapping solution for the 3D mapping and modeling of underground pipe networks during renovation and installation work when the infrastructure is being laid down in open trenches. The system is used to scan the trench and then the 3D scans obtained from the system are registered together to form a 3D point cloud of the trench containing the pipe network using a modified global ICP (iterative closest point) method. In the 3D point cloud, pipe-like structures are segmented using fuzzy C-means clustering and then modeled using a nested MSAC (M-estimator SAMpling Consensus) algorithm. The proposed method is evaluated on real data pertaining to three different sites, containing several different types of pipes. We report an overall registration error of less than 7%, an overall segmentation accuracy of 85% and an overall modeling error of less than 5%. The evaluated results not only demonstrate the efficacy but also the suitability of the proposed solution.

**Keywords:** 3D point cloud; LiDAR; pipes; automatic detection; segmentation; portable 3D scanning system

## 1. Introduction

The underground space in urban areas has played an important role in urban development over recent years. Underground utilities are critical infrastructure for all modern cities. They carry drinking water, storm water, sewage, natural gas, electric power, telecommunications, etc. In most cities, this underground infrastructure reflects the growth and history of the city.

Mapping underground infrastructure has become indispensable for many important tasks pertaining to urban planning, safety and hazard mitigation, construction and maintenance, environmental impact, resource efficiency, etc. However, this offers several technical and operational challenges. The different infrastructure elements are usually placed at different depths resulting in lots of occlusions and overlap. The pipes and cables in the underground networks may well be conceptually simple, but they are narrow objects of different sizes, spanned over large geographical areas and made of multiple segments. For these very reasons the use of 3D technology from different scanning angles is more suited for such a task, as it allows a more accurate 3D mapping and modeling of the complex infrastructure. Several new regulations introduced recently in major cities impose a class-A geo-referencing/mapping (accuracy of 10 cm on each axis) of all the sensitive underground networks in urban areas. This is essential, especially during construction and renovation tasks, where digging in
the vicinity of safety-critical assets requires a high degree of vigilance from the construction workers and machine operators.

Today, for this task, surveyors are involved to conduct a manual survey at various stages when different underground infrastructures are laid. This is not only costly, but it lacks the necessary accuracy and it is also very time-consuming, considering the tight schedule of construction sites. To reduce the impact of these factors, automation is needed. The aim of this work is to develop a portable low cost automated mapping solution for the 3D mapping and modeling of underground networks during renovation and installation work, when the infrastructure is being laid down in open trenches.

## 2. Related Work

Underground mapping of pipe networks have been a topic of interest in the past several years now. Most of the pipe detection and mapping is done over the ground using ground penetrating radar (GPR) such as presented in [1–3]. However, dealing with open ditches with visible pipe networks is a different issue. Some works [4,5] have employed 2D camera images to detect and model these visible pipes in open ditches; however, these methods usually do not meet the accuracy requirements. The use of 3D light detection and ranging (LiDAR) sensors for mapping and modeling such visible pipe networks is quite recent. However, the problem of detecting and modeling pipes in industrial environment has received some attention in the scientific community in recent years. For instance, Rabbani et al. [6] presented a segmentation technique from LiDAR data employing a smoothness constraint based segmentation technique [6]. In this bottom-up approach, a *k*-nearest neighbors (*k*-NN) method is used to fit a plane to the neighbors to estimate the surface normal for each 3D point. The local surface curvature is then approximated using the residual of the plane fitting. Whereas a large residual indicates a curved surface, a small residual indicates a planar surface. A region growing method is then employed, in which points with the smallest residual are used as seed points, with neighboring points having residuals below a predefined threshold added to the current region based on the smoothness constraint. One of the shortcomings of this method is the computational complexity of the algorithm as *k*-NN needs to be computed for every point.

Liu et al. [7] presented a hierarchical structure detection and decomposition method to detect pipe-like objects in LIDAR data in oil refineries. In this method, the ground plane was first extracted using a Gauss map, and then the point cloud was projected onto the ground plane. In this method, the ground plane was first extracted using a Gauss map, and then the point cloud was projected onto the ground plane reducing the problem to finding circles in $R^2$, rather than finding pipes in $R^3$, using a random sample consensus (RANSAC) based method. In the first step, the detected pipes that were perpendicular to the ground plane were removed, and then the remaining points were projected onto different planes perpendicular to the estimated ground plane. A circle fitting technique was used to extract pipes parallel to the ground. The approach, which is computationally more efficient, resembles, in many ways, the joint cylinder segmentation and cylinder fitting approach. Also, similar to our proposed approach, it also considers pipes predominantly parallel to the ground plane.

In [8], an adaptive segmentation and extraction approach for planar and cylindrical features in terrestrial LIDAR data is presented. The method used eigenvectors and eigenvalues, computed from each point’s local neighborhood, to determine whether the points belong to planar or cylindrical surfaces. The local density variation and noise were taken into consideration for robustness. The directional and positional parameters of the cylindrical features were used to generate final segmentation results through clustering.

Masuda [9] proposed a method to recognize planes and cylinders from 3D industrial data. However, the method is semi-automatic as the region of interest has to be manually selected in advance. Matsunuma et al. [10] implemented a similar algorithm but it requires additional information such as brightness/intensity which is not always available. In [11], a skeletonization algorithm is used to classify scanned 3D point clouds into several groups, each corresponding to a single skeleton.
Although the proposed method effectively extracts featured structures and also identifies connection relationships in the piping network, it somewhat struggles to find straight portions of the pipe in the point cloud. The authors of [12] used the skeleton extraction algorithm with Laplacian smoothing to determine the central axis of the pipe. This axis was then divided to estimate the curved region and the T shape of the cylinder. A multiple cylinder detection method that used a coarse-to-fine approach and a clustering algorithm was presented in [13]. As this method employed a 5D Hough transform to estimate the cylinder shape, it was computationally inefficient and so was generally computed separately in 2D and 3D. In [7], the data were first extracted and removed vertically and horizontally from the ground, and then the Hough transform was applied to project a plane on the normal vector followed by a cylinder shape estimation through circle fitting [14,15]. Huang and You [16] proposed a method to classify 3D point clouds using fast point feature histogram (FPFH) and support vector machine (SVM) and then segmented it using a classical flood-fill algorithm. In this, segmented data matching was done using RANSAC. However, the proposed method worked well mainly for relatively complex shapes having unique feature points. Similarly, in [17] primitive models were extracted in point clouds using FPFH, and Adaboost with library models was used to classify objects. The objects were merged and adjusted based on directionality and the start and end positions. Qiu et al. [18] computed a plane from normals (estimated through RANSAC) on the Gaussian sphere. The radius and length of the cylinder were estimated by projecting the segmented points on the plane. In order to estimate the cylindrical parameters of pipes, the authors in [19] employed RANSAC fitting to compute the center and radius of the sphere assuming to be the radius of the cylindrical axis of the pipe. Principal component analysis (PCA) was then employed to find the straight and curved portions of the pipes, using a discriminant based on eigenvalues. Where the curved portion was matched by a Catmull-Rom spline [20], the linear region was taken as the sum of the vectors of the linear candidates.

In addition to LiDARs, some studies have also employed RGB-D cameras, however, to little success. Nahangi et al. [21] proposed a radius estimation algorithm for accurate geometrical pipe feature detection, and using a radius feature metric, evaluated the performance of certain range cameras. The results showed an average error ranging from 10% to 18% at close distance measurements. Proença and Gao [22] proposed a method for fast plane and cylinder extraction to improve visual odometry performance on scenes made of cylindrical surfaces. In order to improve the efficiency and cater for sensor noise, image cells were used instead of 3D points. However, this resulted in surfaces smaller than the patch size to be filtered out. The authors demonstrated, in [23], a technique for the detection and pose estimation of arbitrary cylindrical structures in RGB-D data using a robust soft-voting scheme based on the generalized Hough transform. The rejection of outliers, mainly arising from planar surfaces that contaminate the orientation voting space, was improved by incorporating the curvature information in the voting scheme. However, the results were only evaluated on synthetic data.

The surface obtained after feature extraction needs to be represented and reconstructed. One of the most effective methods is using implicit functions for this task [24,25]. This helps to generate tightly bound 3D models from discrete data points [26], repair missing data points, fill holes, and filter data acquisition noise [27]. However, for larger and more complex point clouds, estimating an accurate implicit equation is difficult. In our work, we divide the point clouds of segmented pipe structures into further subsections to estimate the parameters of the implicit function more accurately, similar to [28], in which the 3D point cloud is first divided into patches, and then each patch is fitted by different implicit functions before the fitting results are combined [29]. This technique, not only results in low computational load and a small amount of data, but it also allows local smoothing [28].

Most of these works either used single scans or an already registered point cloud, minimizing the role of registration, and they also focused more on a structured industrial environment. To the best of our knowledge, no work has directly used terrestrial LiDAR data to detect, map, and model pipes in open ditches in the outdoor environment. The unstructured environment poses several other challenges such as sparsity of point cloud, limited access and visibility, undue clutter especially
due to vegetation and dirt/rubble, unforeseen occlusions, indiscriminating surfaces, etc. In this work, we use 3D data to detect, map, and model underground pipe networks during renovation and installation work when the infrastructure is being laid down/examined in open trenches, using a portable data acquisition system (explained in Section 3.1). Different scans obtained from the system are first registered together using a modified global iterative closest point (ICP) method to form a complete 3D point cloud of the scene (Section 3.2). This 3D point cloud is then analyzed to extract the ditch region, and then, within the ditch, pipe-like structures are segmented, employing fuzzy C-means clustering (Section 3.3). Different parameters of the segmented pipes are estimated using a nested MSAC (M-estimator SAMpling Consensus) method and then used to model the pipe network (Section 3.4). The results are evaluated on real data obtained from different sites to show the efficacy and viability of the method (Section 4). After a thorough analysis of the proposed method and relevant discussion, the conclusion is presented in Section 5.

# 3. Materials and Methods

## 3.1. Portable 3D Data Acquisition System

The portable system used to collect data consists of a VLP-16 [30] laser scanner, a low cost AHRS (Attitude Heading and Reference System) [31], and a simple low cost Global Positioning System (GPS) receiver [32] as shown in Figure 1.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** VLP-16 LiDAR (a), GPS receiver (b), and low cost attitude heading and reference system (AHRS) (c) used in the portable 3D data acquisition system.</figcaption></figure>

The VLP-16 is the third generation laser sensor produced by Velodyne. The salient features of the laser scanner are presented in [30]. Mounted on the same rigid assembly, the AHRS provides the rough orientation of the 3D scans, while the GPS provides its rough position, using rigid transformations of the later sensors in the VLP-16 frame of reference.
Due to the limited vertical FOV (field of view) of the VLP-16 sensor, the user walks around the ditch extending the sensor assembly over it to effectively capture/scan the pipe network inside as shown in Figure 2.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Scanning a ditch using the portable 3D data acquisition system.</figcaption></figure>

The data is acquired and stored on the computer (laptop) which is placed in a backpack along with the portable battery. The total weight of the sensor assembly is about 2 kg which makes it easy to carry and maneuver. The average total acquisition time for a ditch of about 5 m in length is about 4 min.

The rough orientation and position of the 3D scans help to provide an initial guess/alignment for multiple scan registration to form a single registered 3D point cloud, as explained in the next section.

## 3.2. Registration of 3D Scans

The rough alignment with the help of the AHRS and GPS data was not enough to form a coherent registered 3D point cloud. Thus, this initial information was used to finely register the acquired 3D scans to form a registered 3D point cloud using a modified global ICP technique as explained below.

In the global ICP the different 3D scans are finely registered using a type of bundle adjustment method as presented in [33]. The proposed method takes into account all the 3D scans at the same time and searches correspondences in the overlapping regions of different scans. Due to the fact that the LiDAR has a fast scan rate (10 Hz), that the sensor assembly is moved quite slowly over the ditch (normal walking speed), and that the FOV of the sensor remains focused on the ditch during the whole scanning process, there is a large amount of overlap (overlapping regions) found in the different scans and hence several correspondences. The method is summarized in Figure 3.

It is shown in [34] that the registration results are strongly affected by the choice of correspondences. A single correspondence is defined by two points from overlapping regions along with their normal vectors (estimated from the neighboring points using a principal component analysis method [35]). As a point and its normal vector define a tangent plane, thus, a correspondence represents two homologous tangent planes in object space.

The correspondences are established for each pair of overlapping regions in three distinct steps: the Selection, Matching, and Rejection steps.

In the selection step, a subset of points is selected within the overlap region using maximum leverage sampling (MLS). This technique helps select those points, which are best suited for the estimation of the parameters. For this, the effect of each point on the parameter estimation, i.e., its
leverage, is considered. The points with the maximum leverage, which correspond to the lowest redundancy, are selected. This strategy considers the coordinates and the normal vectors of the points.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Summary of the proposed registration method.</figcaption></figure>

The matching step establishes the correspondences. Each selected point from the selection step is paired to the nearest neighbor (the closest point) in the overlapping scans. The nearest neighbor search is efficiently realized using short for *k*-dimensional trees (*k*-d trees). As in the adjustment phase the point-to-plane distance is minimized for each correspondence, two associated points do not have to be identical in the object space, they only have to belong to the same (tangent) plane.
The rejection step aims to identify and subsequently reject unreliable or false correspondences. Based on our application and type of data, each correspondence is tested using the following criteria:

* Plane roughness of corresponding points: In order to minimize the point-to-plane distance, the reliability of both normal vectors must be ensured. This condition is not met if the scanned object cannot be appropriately modeled by a plane, e.g., in the case of vegetation.

* Angle between the normal vectors of corresponding points: To ensure that two corresponding points belong to the same plane, the angle $\alpha$ between the normal vectors of these points should not exceed an upper limit $\alpha_{max}$. The value of $\alpha_{max} = 5^{\circ}$ is used in this work.

* Point-to-plane distance between corresponding points: Apart from a few false correspondences, the a priori point-to-plane distances $d_1, d_2, \dots, d_n$ are assumed to be normally distributed for each individual pair of overlapping scans. A robust estimator for the standard deviation of this contaminated set of correspondences is given by: $\sigma_{med} = 1.4826 \times \text{med}$, where med is the median of the absolute differences (with respect to the median) [36]. In this work, all correspondences with a point-to-plane distance outside the range of $\tilde{d} \pm 3 \times \sigma_{med}$ are rejected, where $\tilde{d}$ denotes the median of the point-to-plane distances of all correspondences that passed the first two criteria.

These criteria ensure that only the most reliable correspondences are retained. As is evident from these criteria, most of the correspondences belonging to structured parts of the environment are used whereas unreliable correspondences belonging to unstructured environment such as vegetation are rejected. Even if some false correspondences still remain after these rejection criteria, a robust adjustment method used helps in their detection and removal.

In the adjustment phase a robust least squares adjustment is performed to minimize the weighted sum of squared point-to-plane distances (which is defined as the perpendicular distance from one point to the tangent plane of the other point):

$$ \sum_{i=1}^{n} (w_i d_i^2) \rightarrow \text{minimum} \qquad (1) $$

where $w_i$ is the weight, and $d_i$ is the point-to-plane distance of the i-th correspondence. The signed point-to-plane distance $d_i$ is defined as the perpendicular distance from a point to a plane. It is conveniently expressed by the Hessian normal form:

$$ d_i = (p_i - q_i)^T \cdot n_i \qquad (2) $$

where $p_i$ and $q_i$ are the corresponding points of the i-th correspondence, and $n_i$ is the normal vector associated to the point $p_i$.

The weight of the correspondences $w_i$ is estimated iteratively (i.e., before each adjustment) from the point clouds itself or, more specifically, from the previously established correspondences. Assuming that the correspondence i belongs to the scan pair k, then its weight is determined by:

$$ w_i = \frac{1}{\sigma_k^2} \qquad (3) $$

where $\sigma_k$ is the $\sigma_{med}$ value of all (non-rejected) point-to-plane distances belonging to the scan pair k.

We find that although we are mainly concerned with 3D points belonging to the ditch, the ditch points are not fully used in the registration phase as many of these are rejected due to the strict rejection criteria. This is why the choice of a longer range sensor such as VLP-16 (range > 100 m) is useful for our application, as it allows scanning parts of more structured environments (for example, building facades, sign boards, and roads), which are predominantly used for registration (see Figures 4 and 5 for some qualitative results).
<figure><img src="image_7.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b)</figcaption></figure>

<caption>**Figure 5.** (a) The unregistered 3D point clouds of our Site-3 in the city of Pérignat-lès-Sarliève (b) the registered 3D point cloud after the proposed registration process. Colors represent the elevation of the 3D points with blue representing the lowest and red the highest elevation.</caption>

### 3.3. Segmentation of 3D Point Clouds

Once the different scans are registered together to form a 3D point cloud, it is segmented to extract objects/points of interests. The segmentation process consists of two main steps:

#### 3.3.1. Semantic Filtering

In order to extract pipes, we are only interested in the 3D points belonging to the ditch. So in order to extract the part of the point cloud belonging to the ditch, we apply certain semantic filters. The filters are based on distance, height, and profile. As we know that the scanning was done around
the ditch itself, we only take into account the 3D points with horizontal and vertical distance of less than 5 m. As the scanning range of the VLP-16 scanner used is >100 m, this pre-filtering helps us to remove large amount of unwanted 3D points from our analysis. This reduced 3D point cloud still contains some unwanted 3D points belonging to the ground/vegetation surrounding the ditch; hence, a profile based filtering is introduced to further narrow down the desired 3D points in the point cloud. Using a simple sweep scan method, as employed in [37], we analyze the height profile (Z direction) of the 3D point clouds with respect to the plane (X,Y) as shown in Figure 6. A crater usually represents the points belonging to the ditch. These 3D points are then extracted and used for further analyses.

<figure><img src="image_12.png" /><figcaption>Figure 6. (a) The delimited zone around the ditch, (b) the extracted ditch points after analysis of height profile in the X-Y plane.</figcaption></figure>

### 3.3.2. Fuzzy C-Means Clustering

Once the 3D points belonging to the ditch are extracted, they are segmented/clustered to form potential objects. In such a non-structured environment, the clustering/segmentation task, especially if only based on 3D information, can be quite difficult. This type of 3D data is widely hampered due to jump edges [38], feature occlusions [39], varying point resolution [40], close objects of large dimensions, and measurement noise [40], etc. For these reasons, a more robust method of modified fuzzy C-means clustering is employed in this work. Compared to other conventional methods, it is more suitable for this type of data. The originally conceived method [41] has been widely implemented (both in original and modified forms) in 2D image segmentation, especially for medical analysis applications [42,43]. Recently it has also been employed in the 3D domain with success [44,45].

Evolved from hard C-means (HCM) [46], fuzzy C-means (FCM) is a clustering method, in which each data point can belong to multiple clusters with varying degrees of membership. HCM allows only binary classification (0 or 1 according to distance) which makes it simple and fast, however, the precision can be lacking. Based on HCM, the proposed FCM method modifies the class membership function from {0, 1} to [0, 1]. As a result, each sample point has a specific membership degree relative to each of the class, and for each point, the sum of these membership degrees (one membership degree for every class) is 1.

Fuzzy clustering aims to minimize the following cost function as presented in [46]:

$$J_m(U, V) = \sum_{j=1}^{n} \sum_{i=1}^{c} u_{ij}^{m} d_{ij}^{2}. \qquad (4)$$

The sample space $X = \{x_1, x_2, ..., x_n\} \subset \mathbb{R}^s$, where $s$ and $n$ are the dimension and the number of samples respectively. $c$ is the number of clustering classes, $m > 1$ is the fuzzy factor, $d_{ij} = \|x_j - v_i\|$ is the distance between sample point $x_j$ and $i$-th cluster center $v_i$. $V = [v_{ij}]_{c \times s}$ with $v_i \subset \mathbb{R}^s$. $U = [u_{ij}]_{c \times n}$
is the fuzzy membership matrix satisfying Equation (4) with $u_{ij}$ being the membership degree of $x_j$ belonging to the $i$-th class satisfying the following conditions:

$$ \left. \begin{array}{l} \sum_{i=1}^{c} u_{ij} > 1, \quad 1 \le j \le n \\ u_{ij} \ge 0, \quad 1 \le i \le c, 1 \le j \le n \\ \sum_{j=1}^{n} u_{ij} > 0, \quad 1 \le i \le c \end{array} \right\} . \qquad (5) $$

Thus, FCM can be considered as a constrained optimization problem with the following main steps [46]:

### Step I Initialization

* Cluster center $V^{(0)}$.

* Iterative index $k = 0$.

* Maximum iteration number $K$.

* Threshold $\epsilon > 0$.

### Step II Calculation of the membership degree matrix $U(k)$

* Given $\forall j, r, i \neq r, d_{ij}(k) > 0$, then

$$ u_{ij} = \frac{1}{\sum_{r=1}^{c} \left[ \frac{d_{ij}(k)}{d_{rj}(k)} \right]^{\frac{2}{m-1}}} \qquad (6) $$

* If there exist $j$ and $r$ values that make $d_{ij}(k) = 0$, then

$$ u_{rj}(k) = 1, \quad u_{ij}(k) = 0 \quad (i \neq r). \qquad (7) $$

### Step III Computation of the cluster center $V_i(k+1)$

$$ V_i(k+1) = \frac{\sum_{j=1}^{n} u_{ij}^{m}(k)x_j}{\sum_{j=1}^{n} u_{ij}^{m}(k)}. \qquad (8) $$

### Step IV Iteration or termination

* If $\|V_i(k+1) - V_i(k)\| < \epsilon$ or $k > K$, iteration terminates at this instant.

* Otherwise, $k = k + 1$ and back to Step II.

FCM increases the optimizing capability of HCM, but it decreases the convergence rate [47]. Choosing a higher value of $c$ results in over-segmentation/clustering of the 3D point cloud. This ultimately helps segment small and very close objects as well.

Once the point cloud is over-segmented into smaller segments, they are grouped together to form larger objects. The smaller segments are grouped together along the longer dimension of the ditch. The direction of the longer dimension is simply found by analyzing the bounding box of the 3D point cloud of the ditch. As we are searching for pipe-like structures, supposedly along the longer dimension of the ditch, subsequent colinear segments are grouped together to form larger segments as shown in Figure 7. These larger segments represent potential pipe-like objects.
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Segmentation results. (a) The over-segmented point cloud after fuzzy clustering, (b) the results after the agglomerative clustering step.</figcaption></figure>

### 3.4. Modeling of Detected Pipes

Once potential pipe-like objects are segmented in the point cloud, they are first analyzed to determine whether they belong to a cylindrical pipe and, if so, their parameters are estimated for modeling.

The 3D points belonging to a segmented object are first divided into a finite number of equal segments along the object's length (i.e., $\frac{Segment\_Length}{l_n}$ where $l_n$ is taken as 25). This helps to better estimate cylindrical parameters for curved pipe segments. A cylinder model is then used with MSAC (M-estimator SAMpling Consensus) [48] to parameterize the potential pipe surface for each segment which best fits the sampled measurements. This is implemented with two nested MSAC subroutines: First is line fitting, and second is the circle fitting of the projection of points onto the plane orthogonal to the line.

In a set of points $\Lambda$, we choose randomly two points $\lambda_a$ and $\lambda_b$ such that $\lambda_a \neq \lambda_b$. A line $u$ is fitted through $\lambda_a$ and $\lambda_b$ and all the points in $\Lambda$ are then projected onto the plane orthogonal to $u$ (i.e., $\text{orth}_u$) to obtain a set $V \in \mathbb{R}^2$. The density of $V$ is computed and a subset $W$ is chosen to separate points of high density. This helps reduce the effect of different types of noise, for example due to branches, during the estimation/parameterization of the cylinder radius. In order to fit a circle to the points in $W$, a MSAC subroutine is then used. Three unique points $w_c$, $w_d$, and $w_e$ are chosen randomly as the minimal sample set (MSS) and a circle (defined by its center $c \in \mathbb{R}^2$) is estimated/parameterized based on the Kasa Method [49]. Using Equation (9), error is computed for all points in $W$ and, based on $||\epsilon_i|| < \epsilon_0$, the inlier set $I$ is determined. The MSAC model is relatively insensitive to the choice of $\epsilon_0$.

$$ \epsilon_i = \sqrt{(w_i - c)^T (w_i - c) - r.} \qquad (9) $$

The outer/top-level MSAC cylinder fit subroutine continues until a stopping criterion is reached. It returns the cylinder parameters which instantiated $I$ with the highest likelihood as the best model. The algorithm is summarized in Algorithm 1.

Once the cylindrical parameters (line and diameter) for each segment is estimated, the different segment lines are concatenated and second order polynomial estimation is used to obtain a smooth center line. The radius $R$ of the curved cylinder pipe is considered as the mean of the radii of all the different segments.

Even though the length of the cylinder pipe is equal to the length of the curve line estimated using the 3D points, the smooth center line is then extrapolated up to the full length of the ditch (limited by the length and width of the ditch) as the pipe curvature and cross-sectional diameter usually remain constant right through its length. This helps to obtain the complete pipe in the ditch and also helps cater for missing pipe segments due to lack of 3D points, occlusions, etc.
**Algorithm 1:** A summary of the proposed pipe modeling method.

**Data:** 3D points of segmented object

**Result:** model parameters *I* and *O*

1 **begin**
2 &nbsp;&nbsp;&nbsp;Divide object into *segments* of equal lengths
3 &nbsp;&nbsp;&nbsp;Use MSAC to fit cylinder
4 &nbsp;&nbsp;&nbsp;**forall** *segments* **do**
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Use MSAC to fit line
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**begin**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Select 2 points $\lambda_a$ and $\lambda_b$ randomly such that $\lambda_a \neq \lambda_b$
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;u = $\lambda_a - \lambda_b$
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;V = {$v_1, v_2, \dots, v_{|A|}$} | V $\in \mathbb{R}^2$, $v_i = \text{orth}_u(\lambda_a - \lambda_b)$
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Estimate the density of V
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;W = {$w_1, w_2, \dots, w_{|W|}$} $\subseteq$ V
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Use MSAC to fit circle
13 &nbsp;&nbsp;&nbsp;**begin**
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Select 3 points $w_c, w_d$ and $w_e$ such that $w_c \neq w_d$ and $w_c \neq w_e$
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Use these 3 points employing Kasa method to parameterize a circle with center
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$c \in \mathbb{R}^2$ and radius $r$
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\epsilon_i = \sqrt{(w_i - c)^T(w_i - c)} - r$
18 &nbsp;&nbsp;&nbsp;**end**
19 &nbsp;&nbsp;&nbsp;Compute inliers *I*: $\lambda_i \in I$ if $|\epsilon_i| < \epsilon_0$
20 &nbsp;&nbsp;&nbsp;Compute outliers *O*: $\lambda_i \in O$ if $|\epsilon_i| \ge \epsilon_0$

## 4. Experiments, Results, and Discussion

In order to evaluate the proposed method, real data were obtained using the portable LiDAR acquisition system explained in Section 3.1 at three different sites located in the city of Grenoble and Pérignat-lès-Sarlièves in France (location shown in Figure 8) at the time when underground pipes were being laid. The sites contained a number of ditches with different underground pipe networks. Qualitative results from some of these ditches are presented in Figures 9 and 10. As seen in the figures, the pipes are generally well segmented and modeled by the proposed method.

For quantitative analyses ground truth was obtained, for each site, using a Leica P-20 total scan station [50]. The scan station consists of an accurate rotating 3D Laser scanner with a calibrated colored camera. The colored 3D scans were obtained from different positions and viewing angles around the ditch and then registered together to form a high resolution 3D colored point cloud [51] for each site as shown in Figure 11. The ground truth was obtained at a high accuracy of less than 5 mm which is sufficient for evaluation purposes. With the help of the ground truth, different quantitative analysis were conducted as explained in the following sections.
<figure><img src="image_14.png" /><figcaption>(a) Site-1 and Site-2 in Grenoble city (FRANCE)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Site-3 in Pérignat-lès-Sarliève (FRANCE)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 9.** *Cont.*</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 9.** Some segmentation results from the three sites. Pipes are segmented in the 3D point clouds (a–f). The point clouds, apart from the segmented pipes, are colored based on the elevation of the 3D points with blue representing the lowest and red, the highest elevation.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 10.** (a–f) Modeling results from the three sites. The modeled pipes are presented in blue in the 3D point clouds that are colored based on the elevation of the 3D points with blue representing the lowest and red, the highest elevation.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 11.** *Cont.*</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 11.** (a,b) The global 3D scan (ground truth) of Site-1 and Site-2 respectively, (c–h) shows the 3D scan (ground truth) of some of the different types of underground pipe networks found on the three sites.</figcaption></figure>

4.1. Registration Accuracy

Some qualitative results for the registration process are presented in Figures 4 and 5. For the quantitative results, as there is no standard method to measure the registration accuracy, we measured it in an indirect manner. In order to assess the registration accuracy, dimensions of different objects and features along the same X,Y, and Z axes in the registered 3D point cloud were analyzed and compared with those in the ground truth. The difference in the dimensions is considered primarily due to the registration errors. The different objects and features 'n' selected for the analysis mainly include large objects from the structured environment present around the ditch, and on the site, such as perpendicularity, length and height of nearby buildings and walls, dimensions of several windows, height and width of piping canals and lighting poles, etc. The dimensions of the same selected features from both the ground truth and the registered point cloud were then compared using the mean absolute error ratio (MAER) calculated for each site, along each of the three axes, as follows:

$$MAER = \frac{1}{n} \sum_{i=1}^{n} \frac{\| Dim_{GT}^{i} - Dim_{M}^{i} \|}{Dim_{GT}^{i}} \qquad (10)$$
where $Dim_{GT}^i$ and $Dim_M^i$ are the i-th dimension measured in the ground truth and the registered point cloud respectively. The results are presented in Table 1.

<table><caption>**Table 1.** The mean absolute error ratio (MAER) along the X, Y, and Z axes for *n* number of reference dimensions chosen in the respective datasets.</caption><thead><tr><th rowspan="2"></th><th colspan="4">MAER</th></tr><tr><th>X Axis</th><th>Y Axis</th><th>Z Axis</th><th>s = √X<sup>2</sup> + Y<sup>2</sup> + Z<sup>2</sup></th></tr></thead><tbody><tr><td>Site-1 (n = 30)</td><td>0.025</td><td>0.031</td><td>0.043</td><td>0.058</td></tr><tr><td>Site-2 (n = 15)</td><td>0.035</td><td>0.048</td><td>0.059</td><td>0.084</td></tr><tr><td>Site-3 (n = 25)</td><td>0.027</td><td>0.039</td><td>0.046</td><td>0.066</td></tr><tr><td></td><td colspan="2"><strong>Average</strong></td><td></td><td>0.069 (= 6.9%)</td></tr></tbody></table>

The higher value of MAER for the registered point cloud of Site-2 is due to the fact that there were not a lot of buildings or other built infrastructure (structured environment) around the ditch area and so as a result the registration process was somewhat hampered due to the rejection criterion explained in Section 3.1.

This shows the importance of having more structured environment around the ditch area to improve the registration quality.

## 4.2. Segmentation and Classification Results

The segmentation and classification results were evaluated at 3D point level. With the help of the available ground truth, the 3D points in the point cloud were manually labeled as belonging to pipes or otherwise. The results expressed in Table 2 are evaluated using different standard evaluation metrics as presented in [52]. Even though all these metrics are commonly used to evaluate such algorithms, the MCC (Matthews correlation coefficient) is insensitive to different class sizes (as is the case with our application, the number of pipe segment points are generally quite inferior as compared to points belonging to the surrounding environment, i.e., non-pipe points) and is considered as the most balanced measure.

Like other metrics, MCC is calculated using the values counts of the true positives (i.e., correct detection of a 3D point belonging to a pipe), false positives, true negatives and false negatives. A coefficient value of +1 signifies a perfect prediction, 0 is no better than random prediction, and -1 implies total disagreement. The results, including an overall accuracy (ACC) and Matthews Correlation Coefficient (MCC) greater than 85% and +0.6, respectively, clearly show the efficiency and usefulness of the proposed method.

<table><caption>**Table 2.** The results evaluated using different standard evaluation metrics.</caption><thead><tr><th></th><th></th><th>Site-1</th><th>Site-2</th><th>Site-3</th></tr></thead><tbody><tr><td>ACC</td><td>Accuracy</td><td>0.874</td><td>0.850</td><td>0.861</td></tr><tr><td>PPV</td><td>Positive Predictive Value</td><td>0.846</td><td>0.810</td><td>0.833</td></tr><tr><td>NPV</td><td>Negative Predictive Value</td><td>0.851</td><td>0.825</td><td>0.841</td></tr><tr><td>FDR</td><td>False Discovery Rate</td><td>0.154</td><td>0.190</td><td>0.167</td></tr><tr><td>F<sub>1</sub></td><td>F<sub>1</sub> measure</td><td>0.710</td><td>0.610</td><td>0.690</td></tr><tr><td>MCC</td><td>Matthews Correlation Coefficient</td><td>+0.631</td><td>+0.602</td><td>+0.629</td></tr></tbody></table>

## 4.3. Modeling Accuracy

The segmented pipes are modeled as cylindrical tubes/pipes. Two main parameters of these modeled pipes are evaluated (i.e., the position of center line and the radius) and compared with the corresponding values obtained from the ground truth. The position of the center line was measured
along the X, Y, and Z axes, using a common point as a measurement reference. The results are evaluated using the standard deviation σ as a measure to evaluate how much the center line of the modeled pipes deviates from the one in the ground truth, and Equation (10) is used to calculate error in the diameter estimation. The results are presented in Table 3.

<table><caption>**Table 3.** The evaluation of pipe modeling accuracy for the three sites.</caption><thead><tr><th></th><th>σ<sub>X</sub> (m)</th><th>σ<sub>Y</sub> (m)</th><th>σ<sub>Z</sub> (m)</th><th>Diameter (MAER)</th></tr></thead><tbody><tr><td>Site-1</td><td>0.017</td><td>0.018</td><td>0.051</td><td>0.051</td></tr><tr><td>Site-2</td><td>0.016</td><td>0.021</td><td>0.045</td><td>0.043</td></tr><tr><td>Site-3</td><td>0.014</td><td>0.018</td><td>0.041</td><td>0.045</td></tr><tr><td></td><td></td><td></td><td>Average</td><td>0.046 (= 4.6%)</td></tr></tbody></table>

A higher standard deviation is observed in the case of the Z-axis due to a lack of 3D points belonging to the bottom part of the pipe, as it is usually occluded during scanning, and also some lower parts of the pipe were extracted out as part of ground points (see Figure 12 for example) and so as a result there is a larger deviation from the true center line along the Z-axis. It is for this reason that the dimensions along the X-Y axis are used to estimate the pipe diameter rather than the height axis Z. The error in the diameter measurement is strongly affected by the segmentation errors, especially due to jump edges between two adjacent objects in the unstructured environment as shown in Figure 12, and also due to occlusions as shown in Figure 13 and discussed in Section 3.3.2. In Figure 13, it can be clearly seen that the vertical pipe structures (white in Figure 13b) are not detected/segmented as they are not considered as pipes by the proposed method, and also due to sand covering parts of the pipe (marked by black circle marks).

<figure><img src="image_21.png" /><figcaption>**Figure 12.** Black circle shows the missing 3D points from the bottom part of the segmented pipe due to occlusion and extraction as ground points whereas blue circle marks the incorrect segmentation with close objects due to jump edges.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 13.** In (a) and (b), black circle marks the portion of the pipe incompletely segmented due to occlusions (covered by sand).</figcaption></figure>

**5. Conclusions**

In this paper, we present a portable, automated mapping solution for the 3D mapping and modeling of underground pipe networks during renovation and installation work, when the infrastructure is being laid down in open trenches. The system is used to scan the trench and the 3D scans obtained from the system are registered together to form a 3D point cloud of the trench containing the pipe network using a modified global ICP method. The pipes are segmented in the resulting 3D point cloud using fuzzy C-means clustering and then modeled using a nested MSAC algorithm. The proposed method is evaluated on real data from three different sites, containing several different types of pipes. An overall registration error of less than 7% is obtained, and it is observed that the proposed registration method gives higher registration accuracy (i.e., lower registration errors) in the presence of structured environment around the site, as compared to less or no structured environment. An overall segmentation accuracy and modeling error of 85% and less than 5% respectively, are also reported. These results are sometimes hampered due to occlusions and jump edges. In order to improve the segmentation and consequently the modeling results, 2D cameras could be added to the portable system in future works. The 2D colored images would be registered with the 3D scans to obtain a colored 3D point cloud. This extra information would provide additional discriminating factors enabling better segmentation [53].

The evaluated results demonstrate the efficacy and the technical prowess of the proposed solution. This solution is also an important step in reducing operational costs and improving mapping accuracy of underground pipe networks.

**Author Contributions:** A.K.A., L.M. and P.C. conceived and designed the experiments; A.K.A. and L.M. performed the experiments; A.K.A. and P.C. analyzed the data; A.K.A. and L.M. contributed materials/analysis tools; A.K.A., P.C. and L.T. wrote the paper; supervision, project administration and funding acquisition, P.C. and T.C.

**Funding:** This work was sponsored by a public grant overseen by the Auvergne-Rhône-Alpes region as part of the Pack Ambition Recherche 2017.

**Acknowledgments:** The authors would also like to thank CEA Tech Grenoble, France, for their cooperation during the data acquisition.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Li, J.; Guo, T.; Leung, H.; Xu, H.; Liu, L.; Wang, B.; Liu, Y. Locating Underground Pipe Using Wideband Chaotic Ground Penetrating Radar. *Sensors* **2019**, *19*, 2913. doi:10.3390/s19132913. [CrossRef]

2. Lai, W.W.L.; Dérobert, X.; Annan, P. A review of Ground Penetrating Radar application in civil engineering: A 30-year journey from Locating and Testing to Imaging and Diagnosis. *NDT E Int.* **2018**, *96*, 58–78. doi:10.1016/j.ndteint.2017.04.002. [CrossRef]
3. Bai, X.; An, W.; Wang, B.; Jiang, J.; Zhang, Y.; Zhang, J. Automatic Identification of Underground Pipeline Based on Ground Penetrating Radar. In Proceedings of the International Conference on Wireless and Satellite Systems, WiSATS 2019: Wireless and Satellite Systems, Harbin, China, 12–13 January 2019; pp. 70–78.

4. FITESIC. RAPH. Available online: http://90.85.0.186/index.php/fr/innovations-technologies/fit-esic-innove/60-innovations-technologies-fr-fr/fit-esic-innove-fr-fr/158-raph (accessed on 16 August 2019).

5. reso3d. RESO 3D. Available online: http://reso3d.com/ (accessed on 16 August 2019).

6. Rabbani, T.; van den Heuvel, F.; Vosselman, G. Segmentation of point clouds using smoothness constraints. In Proceedings of the ISPRS Commission V Symposium Vol. 35, Part 6: Image Engineering and Vision Metrology, International Society for Photogrammetry and Remote Sensing (ISPRS), Dresden, Germany, 25–27 September 2006; pp. 248–253.

7. Liu, Y.; Zhang, J.; Hou, J.; Ren, J.; Tang, W. Cylinder Detection in Large-Scale Point Cloud of Pipeline Plant. *IEEE Trans. Vis. Comput. Graph.* **2013**, *19*, 1700–1707. doi:10.1109/TVCG.2013.74. [CrossRef] [PubMed]

8. Lari, Z.; Habib, A. An adaptive approach for the segmentation and extraction of planar and linear/cylindrical features from laser scanning data. *ISPRS J. Photogramm. Remote Sens.* **2014**, *93*, 192–212. doi:10.1016/j.isprsjprs.2013.12.001. [CrossRef]

9. Masuda, H. Extraction of Surface Primitives using Geometric Constraints of Objects. *Proc. JSPE Semest. Meet.* **2009**, *2009A*, 899–900. doi:10.11522/pscjspe.2009A.0.899.0. [CrossRef]

10. Matsunuma, C.; Abe, Y.; Masuda, H. Reconstruction of Cylinder and Rectangle Faces by Detecting Edges in Large-Scale Point-Cloud. *Proc. JSPE Semest. Meet.* **2010**, *2010A*, 63–64. doi:10.11522/pscjspe.2010A.0.63.0. [CrossRef]

11. Bucksch, A.; Lindenbergh, R.; Menenti, M. SkelTre. *Vis. Comput.* **2010**, *26*, 1283–1300. doi:10.1007/s00371-010-0520-4. [CrossRef]

12. Hyojoo, S.; Changmin, K.; Changwan, K. Automatic 3D Reconstruction of As-built Pipeline Based on Curvature Computations from Laser-Scanned Data. In Proceedings of the Construction Research Congress 2014, Atlanta, GA, USA, 19–21 May 2014; pp. 925–934, doi:10.1061/9780784413517.095. [CrossRef]

13. Su, Y.T.; Bethel, J. Detection and robust estimation of cylinder features in point clouds. In Proceedings of the ASPRS Annual Conference on Opportunities for Emerging Geospatial Technologies, San Diego, California, USA, 26–30 April 2010; p. 7.

14. Patil, A.K.; Holi, P.; Lee, S.K.; Chai, Y.H. An adaptive approach for the reconstruction and modeling of as-built 3D pipelines from point clouds. *Autom. Constr.* **2017**, *75*, 65–78. doi:10.1016/j.autcon.2016.12.002. [CrossRef]

15. Ahmed, M.F.; Haas, C.T.; Haas, R. Automatic Detection of Cylindrical Objects in Built Facilities. *J. Comput. Civ. Eng.* **2014**, *28*, 04014009. doi:10.1061/(ASCE)CP.1943-5487.0000329. [CrossRef]

16. Huang, J.; You, S. Detecting Objects in Scene Point Cloud: A Combinational Approach. In Proceedings of the 2013 International Conference on 3D Vision-3DV 2013, Seattle, WA, USA, 29 June–1 July 2013; pp. 175–182, doi:10.1109/3DV.2013.31. [CrossRef]

17. Pang, G.; Qiu, R.; Huang, J.; You, S.; Neumann, U. Automatic 3D industrial point cloud modeling and recognition. In Proceedings of the 2015 14th IAPR International Conference on Machine Vision Applications (MVA), Tokyo, Japan, 18–22 May 2015; pp. 22–25, doi:10.1109/MVA.2015.7153124. [CrossRef]

18. Qiu, R.; Zhou, Q.Y.; Neumann, U. Pipe-Run Extraction and Reconstruction from Point Clouds. In Proceedings of the Computer Vision (ECCV 2014), Zurich, Switzerland, 6–12 September 2014; pp. 17–30.

19. Jin, Y.H.; Lee, W.H. Fast Cylinder Shape Matching Using Random Sample Consensus in Large Scale Point Cloud. *Appl. Sci.* **2019**, *9*, 974. doi:10.3390/app9050974. [CrossRef]

20. Twigg, C. Catmull-rom splines. *Computer* **2003**, *41*, 4–6.

21. Nahangi, M.; Czerniawski, T.; Haas, C.T.; Walbridge, S. Pipe radius estimation using Kinect range cameras. *Autom. Constr.* **2019**, *99*, 197–205. doi:10.1016/j.autcon.2018.12.015. [CrossRef]

22. Proença, P.F.; Gao, Y. Fast Cylinder and Plane Extraction from Depth Cameras for Visual Odometry. In Proceedings of the 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Madrid, Spain, 1–5 October 2018; pp. 6813–6820. doi:10.1109/IROS.2018.8593516. [CrossRef]

23. Figueiredo, R.; Moreno, P.; Bernardino, A. Robust cylinder detection and pose estimation using 3D point cloud information. In Proceedings of the 2017 IEEE International Conference on Autonomous Robot Systems and Competitions (ICARSC), Coimbra, Portugal, 26–28 April 2017; pp. 234–239, doi:10.1109/ICARSC.2017.7964081. [CrossRef]
24. Fuhrmann, S.; Goesele, M. Fusion of Depth Maps with Multiple Scales. Proceedings of the 2011 SIGGRAPH Asia Conference, New York, NY, USA, 12–15 December 2011; p. 148, doi:10.1145/2024156.2024182. [CrossRef]

25. Tagliasacchi, A.; Olson, M.; Zhang, H.; Hamarneh, G.; Cohen-Or, D. VASE: Volume-Aware Surface Evolution for Surface Reconstruction from Incomplete Point Clouds. *Comput. Graph. Forum* **2011**, *30*, 1563–1571. doi:10.1111/j.1467-8659.2011.02030.x. [CrossRef]

26. Wu, G.; Zhang, Y. A Novel Fractional Implicit Polynomial Approach for Stable Representation of Complex Shapes. *J. Math. Imaging Vis.* **2016**, *55*, 89–104. doi:10.1007/s10851-015-0615-7. [CrossRef]

27. Mizoguchi, T.; Date, H.; Kanai, S.; Kishinami, T. Segmentation of Scanned Mesh into Analytic Surfaces Based on Robust Curvature Estimation and Region Growing. In Proceedings of the Geometric Modeling and Processing (GMP 2006), Pittsburgh, PA, USA, 26–28 July 2006; Kim, M.S., Shimada, K., Eds.; Springer: Berlin/Heidelberg, Germany, 2006; pp. 644–654.

28. Liu, Y.; Song, Y.; Yang, Z.; Deng, J. Implicit surface reconstruction with total variation regularization. *Comput. Aided Geom. Des.* **2017**, *52*, 135–153. [CrossRef]

29. Rouhani, M.; Sappa, A.D.; Boyer, E. Implicit B-Spline Surface Reconstruction. *IEEE Trans. Image Process.* **2015**, *24*, 22–32. doi:10.1109/TIP.2014.2366374. [CrossRef]

30. Velodyne. Velodyne LiDAR PUCK. Available online: https://velodynelidar.com/vlp-16.html (accessed on 16 August 2019).

31. WithRobot. myAHRS+. Available online: http://withrobot.com/en/sensor/myahrsplus/ (accessed on 16 August 2019).

32. ublox. C94-M8P—u-blox RTK Application Board Package. Available online: https://www.u-blox.com/en/product/c94-m8p (accessed on 16 August 2019).

33. Fioraio, N.; Konolige, K. Realtime visual and point cloud slam. In Proceedings of the RGB-D Workshop on Advanced Reasoning with Depth Cameras at Robotics: Science and Systems (RSS), Los Angeles, CA, USA, 27–30 June 2011; p. 3.

34. Glira, P.; Pfeifer, N.; Briese, C.; Ressl, C. A Correspondence Framework for ALS Strip Adjustments based on Variants of the ICP Algorithm Korrespondenzen für die ALS-Streifenausgleichung auf Basis von ICP. *Photogramm. Fernerkund. Geoinf.* **2015**, *2015*, 275–289. [CrossRef]

35. Abdi, H.; Williams, L.J. Principal component analysis. *Wiley Interdiscip. Rev. Comput. Stat.* **2010**, *2*, 433–459. doi:10.1002/wics.101. [CrossRef]

36. Hampel, F.R. The Influence Curve and its Role in Robust Estimation. *J. Am. Stat. Assoc.* **1974**, *69*, 383–393. doi: 10.1080/01621459.1974.10482962. [CrossRef]

37. Aijazi, A.K.; Checchin, P.; Trassoudaine, L. Automatic detection and feature estimation of windows in 3D urban point clouds exploiting façade symmetry and temporal correspondences. *Int. J. Remote Sens.* **2014**, *35*, 7726–7748. doi:10.1080/01431161.2014.975420. [CrossRef]

38. Ghorpade, V.K.; Checchin, P.; Trassoudaine, L. Line-of-sight-based ToF camera’s range image filtering for precise 3D scene reconstruction. In Proceedings of the 2015 European Conference on Mobile Robots (ECMR), Lincoln, UK, 2–4 September 2015; pp. 1–6, doi:10.1109/ECMR.2015.7324208. [CrossRef]

39. Aijazi, A.K.; Checchin, P.; Trassoudaine, L. Handling Occlusions for Accurate 3D Urban Cartography: A New Approach Based on Characterization and Multiple Passages. In Proceedings of the 2012 Second International Conference on 3D Imaging, Modeling, Processing, Visualization Transmission, Zurich, Switzerland, 13–15 October 2012; pp. 112–119, doi:10.1109/3DIMPVT.2012.26. [CrossRef]

40. Aijazi, A.K.; Checchin, P.; Trassoudaine, L. Automatic Removal of Imperfections and Change Detection for Accurate 3D Urban Cartography by Classification and Incremental Updating. *Remote Sens.* **2013**, *5*, 3701–3728. doi:10.3390/rs5083701. [CrossRef]

41. Bezdek, J.C.; Ehrlich, R.; Full, W. FCM: The fuzzy c-means clustering algorithm. *Comput. Geosci.* **1984**, *10*, 191–203. doi:10.1016/0098-3004(84)90020-7. [CrossRef]

42. Halder, A.; Maity, A.; Sarkar, A.; Das, A. A Dynamic Spatial Fuzzy C-Means Clustering-Based Medical Image Segmentation. In *Emerging Technologies in Data Mining and Information Security*; Abraham, A., Dutta, P., Mandal, J.K., Bhattacharya, A., Dutta, S., Eds.; Springer: Singapore, 2019; Volume 813, pp. 829–836.

43. Le, T. Fuzzy C-Means Clustering Interval Type-2 Cerebellar Model Articulation Neural Network for Medical Data Classification. *IEEE Access* **2019**, *7*, 20967–20973. doi:10.1109/ACCESS.2019.2895636. [CrossRef]
44. Cai, Z.; Ma, H.; Zhang, L. A Building Detection Method Based on Semi-Suppressed Fuzzy C-Means and Restricted Region Growing Using Airborne LiDAR. *Remote Sens.* **2019**, *11*, 848. doi:10.3390/rs11070848. [CrossRef]

45. Loutfia, E.; Mahmoud, H.; Amr, A.; Mahmoud, S. 3D model reconstruction from aerial ortho-imagery and LiDAR data. *J. Geomat.* **2017**, *11*, 9.

46. Bezdek, J.C. *Pattern Recognition with Fuzzy Objective Function Algorithms*; Springer Science & Business Media: Berlin, Germany, 2013.

47. Fan, J.L.; Zhen, W.Z.; Xie, W.X. Suppressed fuzzy c-means clustering algorithm. *Pattern Recognit. Lett.* **2003**, *24*, 1607–1612. doi:10.1016/S0167-8655(02)00401-4. [CrossRef]

48. Torr, P.; Zisserman, A. MLESAC: A New Robust Estimator with Application to Estimating Image Geometry. *Comput. Vis. Image Underst.* **2000**, *78*, 138–156. doi:10.1006/cviu.1999.0832. [CrossRef]

49. Kâsa, I. A circle fitting procedure and its error analysis. *IEEE Trans. Instrum. Meas.* **1976**, *IM-25*, 8–14. doi:10.1109/TIM.1976.6312298. [CrossRef]

50. Leica-Geosystems. Leica ScanStation P20. Available online: https://w3.leica-geosystems.com/downloads123/hds/hds/ScanStation_P20/brochures-datasheet/Leica_ScanStation_P20_DAT_us.pdf (accessed on 3 August 2019).

51. Aijazi, A.K.; Malaterre, L.; Tazir, M.L.; Trassoudaine, L.; Checchin, P. Detecting and Analyzing Corrosion Spots on the Hull of Large Marine Vessels using Colored 3D LiDAR Point Clouds. *ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci.* **2016**, *III-3*, 153–160. doi:10.5194/isprs-annals-III-3-153-2016. [CrossRef]

52. Vihinen, M. How to evaluate performance of prediction methods? Measures and their interpretation in variation effect analysis. *BMC Genom.* **2012**, *13*, S2. [CrossRef]

53. Gong, Z.; Lin, H.; Zhang, D.; Luo, Z.; Zelek, J.; Chen, Y.; Nurunnabi, A.; Wang, C.; Li, J. A Frustum-based probabilistic framework for 3D object detection by fusion of LiDAR and camera data. *ISPRS J. Photogramm. Remote Sens.* **2020**, *159*, 90–100, doi:10.1016/j.isprsjprs.2019.10.015. [CrossRef]

<figure><img src="image_23.png" /></figure>

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).