sensors

Article
Change Analysis in Structural Laser Scanning Point
Clouds: The Baseline Method
Yueqian Shen 1,2, *, Roderik Lindenbergh 2 and Jinhu Wang 2,3
1 School of Earth Science and Engineering, Hohai University, No. 1, Xikang Road, Nanjing 210098, China
2 Department of Geoscience and Remote Sensing, Delft University of Technology, Stevinweg 1, 2628 CN Delft,
The Netherlands; [email protected] (R.L.); [email protected] (J.W.)
3 Key Laboratory of quantitative Remote Sensing Information Technology, Academy of Opto-Electronics,
Chinese Academy of Sciences, No. 9, Deng Zhuang South Road, Haidian District, Beijing 100094, China
* Correspondence: [email protected]; Tel.: +86-158-5067-9847

Academic Editor: Dale A. Quattrochi

Received: 13 August 2016; Accepted: 16 December 2016; Published: 24 December 2016
Abstract: A method is introduced for detecting changes from point clouds that avoids registration.
For many applications, changes are detected between two scans of the same scene obtained at
different times. Traditionally, these scans are aligned to a common coordinate system having the
disadvantage that this registration step introduces additional errors. In addition, registration requires
stable targets or features. To avoid these issues, we propose a change detection method based on
so-called baselines. Baselines connect feature points within one scan. To analyze changes, baselines
connecting corresponding points in two scans are compared. As feature points either targets or
virtual points corresponding to some reconstructable feature in the scene are used. The new method
is implemented on two scans sampling a masonry laboratory building before and after seismic testing,
that resulted in damages in the order of several centimeters. The centres of the bricks of the laboratory
building are automatically extracted to serve as virtual points. Baselines connecting virtual points
and/or target points are extracted and compared with respect to a suitable structural coordinate
system. Changes detected from the baseline analysis are compared to a traditional cloud to cloud
change analysis demonstrating the potential of the new method for structural analysis.

Keywords: terrestrial laser scanning; change detection; masonry buildings; baselines; structural
analysis

1. Introduction
Terrestrial Light Detection and Ranging (Lidar) generates 3D coordinates of an object point by
measuring the horizontal and vertical angle and the distance between the scanner’s centre and the
object point. Due to its ability to provide dense, fast and accurate measurements, the use of terrestrial
Lidar for various surveying applications such as deformation monitoring and damage detection has
increased rapidly [1]. Terrestrial Lidar is used in civil engineering applications in structural monitoring
of tunnels [2], road modelling [3], sign inventory [4], evaluating deformations and/or geometric
changes [5], road identification [6], drift detection [7], natural hazards and structural analyses in
cultural heritage [8,9]. The application of terrestrial Lidar data for accurately measuring deflection of
loaded beams is discussed by Park et al. [10].
The density of points on the object’s surface can be predefined by the user and is limited by the
minimum angle increment of the system. Depending on the distance from the scanner and the amount
of scans, a very high point density can be achieved. In addition, the reflectance of the surface may be
measured by recording the intensity of the reflected laser beam.
Combined with high quality static GNSS positioning and precise tachometry, terrestrial Lidar
was used for high precision monitoring of deformations during longer time intervals [11]. A terrestrial

Sensors 2017, 17, 26; doi:10.3390/s17010026 www.mdpi.com/journal/sensors

Sensors 2017, 17, 26 2 of 25

Lidar-based vehicle detection approach was proposed by using the so-called Probability Hypothesis
Density [12]. Automatic processing of point clouds was carried out to so-called build geometric
models suitable for structural analysis purposes [13]. A method used the triangulation, reflectance
and RGB triplets to obtain a proximity-based segmentation [14]. In recent decades, the advanced
analysis of masonry structures in large point clouds also received considerable attention [15]. A novel
segmentation algorithm was proposed to enable the automatic segmentation of masonry blocks from
3D point clouds acquired by terrestrial Lidar technology [16]. Despite all the improvements, the
challenges in the seismic assessment of these structures remains rather difficult [17]. What’s more,
these type of structures are quite large due to the numerous variations of masonry, the large scatter of
in situ material properties, and the impossibility of reducing it all in a specimen.
Change detection on point clouds is a rather new technique considering that it has been mainly
used by professional surveyors until now [18]. Comparison of the surface geometry in different
epochs is often performed by reconstructing the surface models in each epoch [19,20]. By comparing
the resulting surface models, useful change information is obtained. Combined with conventional
surveying devices such as total station and GPS, georeferenced topographic data is acquired by
terrestrial Lidar for further scene comparison [21]. Subtraction of a resampled data set was used
to detect changes on a hydropower station [22]. One common way to perform change detection is
to compare the 3D coordinates of corresponding points from two epochs. Before the comparison,
two point clouds should be aligned in the same coordinate system which is called registration.
Therefore, registration is a crucial step in change detection applications. Point clouds can be registered
by iteratively decreasing the distances to arbitrary closest points in overlapping areas, which is the
basis of Iterative Closest Point (ICP)-like algorithms [23,24], or by matching explicitly derived features
points [25–28]. However, even minor misalignments of two epochs may lead to erroneous results
when detecting changes. Another common technique, to perform change detection on a point cloud,
which also requires registration in advance, is to compute its distance to a 3D reference model. This can
be done either directly or by creating an intermediate model on top of the points. The reference
model can either be theoretical or also created from real data. Cloud-to cloud and cloud-to-mesh
distances have been very well studied based on above principles. The academic software, “Metro”
and “Mesh”, allows one to compare the difference between a pair of surfaces by adopting a surface
sampling approach [29]. An efficient method was proposed to estimate the distance between discrete
3D surfaces represented by triangular 3D meshes which was based on an approximation of the
Hausdorff distance [30]. Several simple cloud-to-cloud comparison techniques based on a specific
octree structure were presented and implemented [31]. Although the feasibility and effectiveness of
the cloud-to-cloud and cloud-to-mesh distances has been demonstrated in previous work, the absolute
value of displacement can’t be estimated.
To overcome these issues we propose instead to compare corresponding baselines extracted
independently in each of the two scans. Here a baseline is defined as a 3D line segment connecting
two points in one scan. Two baselines from different scans are corresponding if they connect the
same features. Feature identification and extraction have drawn many scholars’ attention in recent
years. A laser-based approach for door and handle identification was proposed and implemented [32].
Euclidean was applied for tabletop object detection which was efficient and enables real-time plane
segmentation [33]. The Fast Point Feature Histograms was used for labelling 3D points with different
geometric surface primitives [34]. As features we propose to use two type of points. First, target points
identified by spherical or planar targets as placed by the surveyor in the scene, and, second so-called
virtual points, which give the 3D location of a feature that is well-reconstructable from the 3D scan
data. What virtual points are suitable depends on the particular scene that is considered.
This proposed methodology is illustrated to detect changes on a masonry building during seismic
testing, mainly for monitoring purposes but also to be able to provide change information for further
structural analysis. In this case, as in most cases, target points are easily identified in the scan data.
The extraction of virtual points requires more work. In this case we choose as virtual points the centres
Sensors 2017, 17, 26 3 of 25

Sensors 2017, 17, 26 3 of 25

of bricks. In the building scan, one wall that is expected to be stable is available, while another wall is
actually cracked
centres during
of bricks. In thethe seismic
building testing.
scan, one wall that is expected to be stable is available, while another
The isextraction
wall of virtual
actually cracked points
during thestarts
seismicwith separation of mortar and bricks in the scan data using
testing.
k-meansThe clustering of of
extraction the TLS intensity
virtual points startsattribute. Next, point
with separation clouds
of mortar andare segmented
bricks in the scananddata3D virtual
using
point locations
k-means are estimated
clustering of the TLS at the brick attribute.
intensity centres. Given both clouds
Next, point the resulting target points
are segmented and 3Dand virtual
virtual
point
points, locations
baselines areconstructed.
are estimated atFinally,
the brickbycentres. Givencorresponding
comparing both the resulting target points
baselines from theandtwovirtual
epochs,
points,
changes baselines
in X, Y and Zare constructed.
direction Finally, structural
of a suitable by comparing corresponding
coordinate system arebaselines from the two
extracted.
epochs, changes in X,ofY the
As an evaluation and proposed
Z directionapproach,
of a suitableresults
structural
are coordinate
comparedsystem are extracted.
to results from traditional
As an evaluation of the proposed approach, results are compared
methods. That is, for a traditional approach, scan data is first registered in a common to results from traditional
coordinate
methods. That is, for a traditional approach, scan data is first registered in a common coordinate
system. Next, virtual point changes and cloud-to-cloud distances between the aligned point clouds
system. Next, virtual point changes and cloud-to-cloud distances between the aligned point clouds
are estimated. As a first benefit of the proposed method we were able to identify a target that
are estimated. As a first benefit of the proposed method we were able to identify a target that was
was apparently
apparentlymoved moved by by an external
an external agentagent
early early
in the in the reconstruction.
reconstruction. A analysis
A detailed detailedofanalysis
the changeof the
change
results will be presented in the results section. To summarize, the proposed method is a new,new,
results will be presented in the results section. To summarize, the proposed method is a
alternative approach
alternative approach forfor
change
change detection
detectionthatthat eliminates
eliminates ananoften
oftenunnecessary
unnecessary registration
registration stepstep
andand
its associated errors.
its associated errors.

2. Materials andand
2. Materials Methods
Methods
In order to detect
In order changes
to detect changesbybybaselines,
baselines, aa workflow
workflow isis introduced
introducedininFigure
Figure
1. 1.
TheThe
left left
box box
summarizes
summarizes thethe
newnew baseline
baseline methodwhile
method whilethe
theright
right box shows
showsthethecomparison
comparison approaches.
approaches.

Figure1.1.Change
Figure Change detection workflow.
detection workflow.

2.1. Seismic
2.1. Seismic Experiment
Experiment Descriptionand
Description andScan
ScanData
Data Acquisition
Acquisition
A seismic experiment has been carried out on a masonry house in the Stevin lab of TU Delft. The
A seismic experiment has been carried out on a masonry house in the Stevin lab of TU Delft.
purpose of the experiment was to explore how resilient houses are against earthquakes. The house
The purpose of the experiment was to explore how resilient houses are against earthquakes. The house
was built of bricklayers with mortar and calcium silicate. The experiment was performed from 3
was built of bricklayers with mortar and calcium silicate. The experiment was performed from
Sensors 2017, 17, 26 4 of 25

3 December 2015 to 16 December 2015. During the experiment the house was shaken back and forth
Sensors 2017, 17, 26 4 of 25
repeatedly and repeated scanning was carried out.
A December
C10 Scan2015 Station scanner (Leica,
to 16 December Heerbrugg,
2015. During Switzerland)
the experiment the house was used for
was shaken backthis experiment.
and forth
The Leica C10 Scan
repeatedly andStation
repeatedscanner
scanningiswasa time-of flight scanner with an effective range of 300 m at 90%
carried out.
reflectively. ASpecifications
C10 Scan Station scannerthat
indicate (Leica,
theHeerbrugg,
accuracySwitzerland) was used for thisisexperiment.
of a single measurement 6 mm (oneThe sigma) in
Leica C10 Scan Station scanner is a time-of flight
position and 4 mm (one sigma) in depth at ranges up to 50 m [35]. scanner with an effective range of 300 m at 90%
reflectively. Specifications indicate that the accuracy of a single measurement is 6 mm (one sigma) in
The location of the scan positions should be planned such that the number of scans is minimized
position and 4 mm (one sigma) in depth at ranges up to 50 m [35].
while avoiding occlusions
The location in order
of the scan to ensure
positions should be the full coverage
planned of number
such that the the whole monitoring
of scans is minimized area [36].
Actually, oneavoiding
while stationocclusions
is enough in to
ordercover the façade
to ensure the full of the moving
coverage masonry
of the whole housearea
monitoring in this
[36]. project.
However, we planned
Actually, one station twoisstations,
enough toone onthe
cover thefaçade
left side of moving
of the the house and the
masonry houseother on project.
in this the right side,
However, we planned two stations, one on the left side of the house
to get additional information about the changes. The positions of the moving wall, stable and the other on the right side,wall, the
to get additional information about the changes. The positions of the moving wall, stable wall, the
scanner and the targets are shown in a top view of the experiment scene, see Figure 2. The left side of
scanner and the targets are shown in a top view of the experiment scene, see Figure 2. The left side of
the scene contains a big door and lacks stable locations for positioning further targets. But, as only the
the scene contains a big door and lacks stable locations for positioning further targets. But, as only
comparison methodsmethods
the comparison need registration, the used
need registration, target
the used distribution
target distributionhas noeffect
has no effectonon
thethe results
results of of the
proposedthemethod.
proposed method.

Figure Figure
2. A 2. A sketch
sketch showingthe
showing the position
position ofof
thethe
testtest
building (moving
building wall), a stable
(moving wall),wall and scan
a stable wall and
positions.
scan positions.
Figure 3 is a picture of the moving wall taken from the perspective of Station 1. Considering the
Figure 3 is aofpicture
structure ofareas
the wall, the moving
of interestwall taken
should from
cover thethe
wallperspective of Station
in a representative way.1. Therefore,
Considering the
structurepatches
of theA,wall,
B, Careas
and Dofatinterest
the top-left,
shouldtop-right,
cover bottom-left
the wall inand bottom-right ofway.
a representative the moving wall patches
Therefore,
respectively, are chosen as a source for virtual points, see Figure 3. Here we should
A, B, C and D at the top-left, top-right, bottom-left and bottom-right of the moving wall respectively, demonstrate that
if we choose the whole wall which contains too much outliers as well as the crakes and it will
are chosen as a source for virtual points, see Figure 3. Here we should demonstrate that if we choose
definitely complicate the pre-processing. Moreover, four patches are enough to show the structural
the whole wallofwhich
change contains
the wall. Actually,tooourmuch
paperoutliers as well as thethe
aim at demonstrating crakes andbaseline
proposed it will definitely
method, socomplicate
we
the pre-processing.
only selected the Moreover,
above four four patches are enough to show the structural change of the wall.
patches.
Actually, ourTwopaper
scansaim atmade
were demonstrating
from each scan theposition.
proposed Thebaseline
first scan method, so we only
was at a minimal selected
resolution, whichthe above
correspond
four patches. to 0.2 m in horizontal and vertical spacing when the range is 100 m [37], with an entire
field of vision (360 degree with respect to Z-axis and 270 degree for vertical amplitude) which is used
Two scans were made from each scan position. The first scan was at a minimal resolution, which
for obtaining a general frame of the whole scanner field of view, where smaller areas can be selected
correspond to 0.2 m in horizontal and vertical spacing when the range is 100 m [37], with an entire
for subsequent scans. The second high-resolution scan, which corresponds to 0.05 m in horizontal
field of and
vision (360spacing
vertical degreeatwith respect
a range of 100tomZ-axis and
[37], was of 270 degree
an area for vertical
including amplitude)
the moving which
wall, stable wall is used
for obtaining a general
and targets. Due toframe of theranges
the different whole to scanner field
the scanner, theof view, where
horizontal smaller
and vertical areas
spacing at can
patchbeA,selected
patch B, patch C and patch D is 0.0055 m, 0.0035 m, 0.0047 m and 0.0028 mm,
for subsequent scans. The second high-resolution scan, which corresponds to 0.05 m in horizontal andrespectively.
vertical spacing at a range of 100 m [37], was of an area including the moving wall, stable wall and
targets. Due to the different ranges to the scanner, the horizontal and vertical spacing at patch A, patch
B, patch C and patch D is 0.0055 m, 0.0035 m, 0.0047 m and 0.0028 mm, respectively.
Sensors 2017, 17, 26 5 of 25
Sensors 2017, 17, 26 5 of 25
Sensors 2017, 17, 26 5 of 25

Figure
Figure 3.3.Station
Stationview
view of
of the
the seismic
seismictesting
testinghouse.
house.
Figure 3. Station view of the seismic testing house.
Two epochs of TLS data (on 3 December 2015 and 16 December 2015) were obtained before and
Two epochs of TLS data (on 3 December 2015 and 16 December 2015) were obtained before and
afterTwo
the seismic
epochs oftesting,
TLS data see (on
Figure 4, in which
3 December 2015the
and induced damage
16 December is apparent.
2015) Each before
were obtained epoch and
was
afterafter
the the
scannedseismic
from testing,
almost
seismic thesee
testing, seeFigure
same place.4,4,
Figure ininwhich
Together
which the
with induced
thethe damage
3D coordinates
induced damage isis
of theapparent.
points in
apparent. Each
the
Each epoch
waswas
monitoring
epoch
scanned from
surface,
scanned thealmost
from almostthethe
coordinates same
of the
same place.
place.Together
reference with
targets
Together were
with thealso
the 3Dcoordinates
3D coordinates
obtained. ofof
They the
are
the points in the
distributed
points in the monitoring
away from
monitoring
surface, the the
the moving
surface, coordinates
coordinatesof the
wall so that reference
ofthere
the reference targets
are at least were
four
targets also
common
were alsoobtained.
targets They
Theyare
between
obtained. distributed
scans
are away
of different
distributed from the
epochs,
away from
moving wall Figure
compare
the moving so wall
that2.there
so thatare at least
there four
are at common
least targetstargets
four common between scans of
between different
scans epochs,
of different compare
epochs,
Figure 2.
compare Figure 2.

Figure 4. On the left: the masonry wall before the seismic testing; on the right: the damaged masonry
wall
Figure 4.after
Figure 4. On
On thethe
the seismic
left:
left: the testing.
the masonrywall
masonry wallbefore
before the seismic
seismictesting;
testing;ononthe right:
the thethe
right: damaged masonry
damaged masonry
wall after the seismic
wall after the seismic testing.testing.
2.2. Change Detection Based on Baselines
2.2. Change Detection Based on Baselines
2.2. Change Detection Based on Baselines
2.2.1. Coordinate Transformation
2.2.1. Coordinate Transformation
2.2.1. Coordinate
In order toTransformation
obtain more intuitive change information, a structural coordinate system is designed.
In order
The initial to obtain more
coordinates, X ' ,intuitive
Y ' and change
Z ' values
information,
obtained abystructural
the TLS, coordinate
are determinedsystemby is
thedesigned.
location
In order to obtain more 'intuitive ' change
' information, a structural coordinate system is designed.
of the
The initial (i.e., distance Xof ,TLS
TLS coordinates, Y to theZobject
and andobtained
values its orientation
by theangle),
TLS, are see the left figure
determined by thein Figure
location5.
The initial coordinates, X 0', Y 0 'and Z 0 values obtained by the TLS, are determined by the location of
of
Displacement the X of
the TLS (i.e.,indistance , YTLSand Z ' object
to the anddoes
direction its orientation angle),tosee
not correspond thethe left figure
expected in Figure
changes 5.
of the
the TLS (i.e., distance of TLS to the object and its orientation angle), see the left figure in Figure 5.
Displacement in the X ' , Y ' and Z ' direction does not correspond to the expected changes of the
Displacement in the X 0 , Y 0 and Z 0 direction does not correspond to the expected changes of the
structure (i.e., along or perpendicular to the building walls) which should be measured in the X, Y and
Sensors2017,
Sensors 2017,17,
17,26
26 66of
of25
25

structure (i.e., along or perpendicular to the building walls) which should be measured in the X , Y
Z
and Z direction
direction of the structure. Therefore,
of the structure. we define
Therefore, a structural
we define coordinate
a structural system
coordinate as follows,
system see the
as follows, see
the right figure in Figure 5. Here, the Y and Z axes are parallel to the building walls. The Z axis
right figure in Figure 5. Here, the Y and Z axes are parallel to the building walls. The Z axis is pointing
upwards
is pointingand the X axis
upwards andisthe
pointing outisthe
X axis building.
pointing out the building.

Figure 5.
Figure 5. TLS
TLS (left)
(left) and
and structural
structural (right)
(right) coordinate
coordinatesystem.
system.

The 3D transformation from the TLS coordinate system to the structural coordinate system
The 3D transformation from the TLS coordinate system to the structural coordinate system consists
consists of six parameters, three translations (i.e., three parameters in three orthogonal directions)
of six parameters, three translations (i.e., three parameters in three orthogonal directions) and three
and three rotational parameters. The translation moves 0 the 0 origin
0 0 O '  x0 ' , y0 ' , z0 '  in the TLS
rotational parameters. The translation moves the origin O ( x0 , y0 , z0 ) in the TLS coordinate system
coordinate
to the originsystem
O( x0 , yto the origin O  x0 , y0 , coordinate
0 , z0 ) in the structural z 0  in the structural
system, see coordinate
Figure 5. system, see Figure 5.
The
The testing
testing building
building is is part
part of
of aa larger
largerscene.
scene. The
The structural
structural coordinate
coordinate system
system is is established
established onon
the
the stable
stable place
place such
such as as laboratory
laboratory ground,
ground, ceiling
ceiling surface,
surface, the
the experimental
experimental steel frame, fixed wall,
and
andetetal.
al.In
Inour
ourprocessing,
processing,we weselect
selectthe fixed
the fixedwall, which
wall, which is parallel to the
is parallel testing
to the testing wallwall surface, to set
surface, to
up the structural coordinate system, compare Figure 2. The procedure
set up the structural coordinate system, compare Figure 2. The procedure is as follows: is as follows:

1.
1. Select the
Select the point
point cloud
cloud ofof the
the fixed
fixed wall.
wall.
2.
2. Fit a plane to the point cloud using Principle Component
Fit a plane to the point cloud using Principle Component Analysis
Analysis (PCA).
(PCA).
The PCA method is used to estimate the point clouds normal vector in our work [38–40].
3. The PCA
Project method
the originisofused
the to estimate
TLS the point
coordinate clouds
system normal
to the vector
plane and in our work
calculate [38–40].
the translational
3. parameter.
Project the origin of the TLS coordinate system to the plane and calculate the translational parameter.
'
The
The plane
plane equation
equation of
of the
the fixed
fixed wall
wall is
is obtained
obtained by by the
the last
last step.
step. Afterwards,
Afterwards,the origin O
theorigin O0 of
of
the
the TLS
TLS coordinate
coordinate system
system isis projected
projected toto the
the fixed
fixed wall,
wall,being
beingthe origin O
theorigin of the
O of the structural
structural
coordinate system. Assuming that the plane equation
coordinate system. Assuming that the plane equation of the of the fixed wall and the origin
origin ofof TLS
TLS
coordinate
coordinate system are expressed as:
Ax  By  Cz  D  0
Ax'
+ By + Cz + D = 0 (1)
O (1)
O0 =(0,
(0,0,0,0)0)
The origin of the structural coordinate system is computed as:
The origin of the structural coordinate system is computed as:
  AD  BD CD 
O 2 2 2
, 2 2
, 2   (2)
 A −BAD  C A  B− BD  C A  B 2−
2
C2 
 CD
O , , (2)
4. A2 +and
Calculate the rotational moves B2 + C2 A2 +the
transform B2 coordinates.
+ C 2 A2 + B2 + C 2

4. The station
Calculate theisrotational
always level
moves during scanning.the
and transform Therefore, the structural coordinate system is
coordinates.

established by rotating in the X 'OY ' '
plane. The normal vector n of the stable wall is regarded
The station is always level during scanning. Therefore, the structural coordinate system  is
as the positive X axis. The rotation angle  is the angle between
0 O0 Y 0 plane. → the normal vector n and
established by rotating
 in the X The normal vector n of the stable wall is regarded as
the X ' direction n X  [1, 0, 0] . Therefore, the 3D coordinate transformation is defined by:
'
Sensors 2017, 17, 26 7 of 25

→
the positive X axis. The rotation angle α is the angle between the normal vector n and the X 0
→
direction n X 0 = [1, 0, 0]. Therefore, the 3D coordinate transformation is defined by:
      
X cos α − sin α 0 X 0 ∆X
Y  =  sin α cos α 0 Y 0  +  ∆Y  (3)
      
Z 0 0 1 Z0 ∆Z

Here, α represents the rotation angle between X 0 O0 Y 0 and XOY; [∆X∆Y∆Z ]0 denotes the
 0
translation vector, and in this paper, [∆X∆Y∆Z ]0 = A2 +AD
2 2 , 2
BD
2 2 , 2
CD
2
B +C A + B +C A + B +C 2 .

2.2.2. Automatic Virtual Point Extraction

As introduced above, a baseline is defined by connecting two feature points, thus it is necessary
to select and extract obvious feature points (e.g., brick centres, door knobs, corner of the wall) from
different areas, including moving areas, stable area and targets.
In this research, we choose brick centres as our virtual points. To implement the proposed method,
bricks without damage are selected and three assumptions are made about the bricks after the seismic
test: (1) a brick is always rectangular; (2) no crack emerges; (3) brick movement was translational, not
rotational. These assumptions are easily met in this case, but particular assumptions are case study
dependent. In this study, each brick is sampled by several random points so we have to reconstruct
corresponding points in corresponding bricks. To do so we perform three steps to be discussed in detail
(1) separation of mortar and bricks; (2) construction of a histogram of mortar points; (3) estimating
brick centres. These steps will enable us to extract 3D brick centres from the brick point clouds for
detecting changes:

1. Separation of Mortar and Bricks Using K-means Clustering

Laser scanners not only provide information about the geometric position of a surface, but also
information about the portion of energy reflected by the surface, which depends on its reflectance
characteristic. The backscatter generated after collision of the laser beam with the object surface is
recorded by most terrestrial Lidar instruments as a function of time [41]. In our work, the mortar
and the brick are composed of different materials, so the signal intensity attribute is considered
a possibility to segment bricks from mortar using k-means clustering, a commonly used data
clustering technique for performing unsupervised tasks. It is used to cluster n objects of the input
dataset into k homogeneous partitions, k < n [42,43]. In our study, the wall is composed of brick
and mortar so we use this technique in a classic way with k = 2 to separate brick points from
mortar points based on their intensity.

2. Histogram of Y and Z Axes Using Mortar Points

The brick point clouds are obtained by the above steps but with no accurate boundaries between
different bricks. The wall is coursed masonry which means the masonry and the mortar are almost
level to the wall surface. In addition, the 3D coordinates of the points have been transformed to
the structural coordinate system as introduced in Section 2.2.1, which indicates that the X axis is
perpendicular to the plane of the masonry wall. Therefore, a histogram along the Y and Z axis is
used to estimate boundary lines between different bricks. The possible horizontal and vertical
boundary lines of the bricks are determined by the peaks of the histogram along the Y and Z axis
using the mortar points extracted in step 1. The procedure of determining the lines is as follows:
As shown in Figure 6a, a window width Lwindow and a step width Lstep is defined along the Y
or the Z axis which satisfies the equation on the width of the mortar Lmortar between two bricks
as follows:
Lmortar ≈ Lwindow + 2Lstep (4)
Sensors 2017, 17, 26 8 of 25

Sensors 2017, 17, 26 8 of 25

In our work, the number of window positions along Y-direction and Z direction is given by:
In our work, the number of window positions along Y-direction and Z direction
 is given by:
ymax − ymin − Lwindow zmax − zmin − Lwindow
  
ny = + 1, nz = +1 (5)
 ymax Lstep
ymin  Lwindow   zmax  zmin  Lwindow
step 
ny     1, nz   1 (5)
 Lstep   Lstep 
In which the sign [] means that the value is rounded up toward the nearest integer.
In which the sign   means that the value
is rounded up toward the nearest integer.
The number of points nyi i = 1, 2, · · · , ny /nzi (i = 1, 2, · · · , nz ) in each window is computed and
The(number
the of points nofyi the
x, y, z) coordinates i  1,2, y  / stored.
 , nare
points n zi  i  1,2,  , n z  in each window is computed and the
( x, y“line
The , z ) coordinates of the
density” along thepoints are stored. is expressed as:
Y- or Z-direction
The “line density” alongthe Y- or Z-direction isexpressed as:


Density_y = n + n +nny(i+1/) 3/L(3Lwindowi)i2,=3,2,
Density _ y y(ni−y (1i )1)  nyi
3,, (·n· ·,1)(ny − 1)
yi y ( i 1)   window    y  (6)
Density_z = nz(i−1) + nzi + nz(i+1) /(3Lwindow )(i = 2, 3, · · · , (nz − 1)) (6)
Density _ z   nz ( i 1)  nzi  nz ( i 1)  /  3 Lwindow   i  2, 3,  , ( nz  1) 

For each
For each window,
window, the
the line
line density
density gradient
gradient along
along the
the Y-direction
Y-directionor
orZ-direction
Z-directionisiscalculated
calculatedby:
by:
Grad  i ,1  Density _ y (i )  Density _ y (i  1)
Grad(i, 1) = Density_y(i ) − Density_y(i − 1) (7)
Grad  i , 2   Density _ y (i  1)  Density _ y (i ) (7)
Grad(i, 2) = Density_y(i + 1) − Density_y(i )
When Grad  i,1  0 and Grad  i,2  0 , the window of the mortar is determined. For the
When
analyzed Grad (i, 1) > 0itand
directions, Gradhas
always < 0, the
(i, 2a) higher window
density when of the
the mortar
mortar is is determined.
perpendicularFor to the
this
analyzed directions, it always has a higher density when the mortar
direction. Therefore, a density threshold is considered when the mortar window is estimated. is perpendicular to this
In
direction. Therefore, a density threshold is considered when the mortar
our work, the average number in each window width is taken as the threshold. That means that window is estimated.
In our work, the average number in each window width is takennas the threshold. That means
total n
along
that the the
along Y-direction
Y-directionandand
Z-direction,
Z-direction, thethethresholds
thresholdsare areεyy = ntotal and  ntotal
and εzz = total , which
in which
ny n y nnz z
ntotal represents the number of points belonging to this patch.
ntotal represents the number of points belonging to this patch.
3. Estimating the Brick Centre Position
3. Estimating the Brick Centre Position
The centre line position is estimated by calculating the average of the points in the mortar
The centre
window as line
shownposition is estimated
in Figure 6b. Fourbylines calculating the average
surrounding one brick of theare points
estimatedin the by mortar
above
window as shown in Figure 6b. Four lines surrounding one brick
procedure. Given the equations of the four lines, four intersection points are estimated. are estimated by above
Next,
procedure. Given the equations of the four lines, four intersection points
the 3D coordinates of the brick centre are calculated by averaging the four intersection points. are estimated. Next,
the 3D coordinates of the brick centre are calculated by averaging the four
Finally, the extracted virtual points are sorted to facilitate automatic identification of virtual intersection points.
Finally,from
points the extracted virtual
the same brick in points
differentareepochs.
sorted to
The facilitate
regular automatic
configuration identification
of the bricks of virtual
in the
points from the same brick in different epochs. The regular
masonry wall facilitates the automatic handling of the brick centre locations. configuration of the bricks in the
masonry wall facilitates the automatic handling of the brick centre locations.

(a) A sketch illustrating the parameters used (b) Extracting brick centres from mortar lines
in the extraction of the mortar lines
Figure6.6.Automatic
Figure Automaticvirtual
virtualpoint
pointextraction.
extraction.

2.2.3. Baseline Establishment

In our work, a baseline is defined as a 3D line segment connecting two points in one scan. The
points represent features such as points identified by spherical or planar targets placed by a surveyor
 Sensors 2017, 17, 26 9 of 25

2.2.3. Baseline Establishment

Sensors In
2017,
our17, 26
work, a baseline is defined as a 3D line segment connecting two points in one9 of 25
scan.
The points represent features such as points identified by spherical or planar targets placed by a
insurveyor
the sceneinandthe so-called
scene andvirtual points
so-called extracted
virtual pointsfrom the 3D
extracted scanthe
from data.
3D We
scanpropose
data. We topropose
use two totypeuse
oftwo
baselines for change detection which are established by connecting points that are
type of baselines for change detection which are established by connecting points that are expected expected to
indicate structural change of the moving wall. First, points subject to change are connected
to indicate structural change of the moving wall. First, points subject to change are connected to points to points
that
thatare
areexpected
expectedtotobe bestable,
stable,and,
and,second
secondpoints
pointssubject
subjecttotochange
changeareareconnected
connectedto toeach
eachother.
other.
By connecting two different points, a baseline is established. In Figure 7a,
By connecting two different points, a baseline is established. In Figure 7a, assuming assuming that Athat
1 and A1
Band
1 represent
B 1
the
represent points
the in
points epoch
in I while
epoch I A
while
2 and
A 2
B
and 2 Bare
2
the
are corresponding
the corresponding points
points ininepoch
epoch II, the
II, the
baselines A
baselines A1B1 1and
1B and A A22BB22 are
are the
thecorresponding
correspondingbaselines
baselineswithin
within one
one scan
scan in in different
different epochs.
epochs. ThisThis
step
is performed in an automated way, which is possible
step is performed in an automated way, which is possible due to the due to the previous sorting of the extracted
sorting of the extracted
virtualpoints.
virtual points.

(a) Baseline establishment (b) Baseline decomposition and comparison

Figure
Figure7.7.Baseline
Baselineestablishment,
establishment,decomposition
decompositionand
andcomparison.
comparison.

2.2.4.
2.2.4.Baseline
BaselineDecomposing
Decomposingand
andComparison
Comparison
For
Foreach
eachpair
pairofofcorresponding
correspondingbaselines
baselinesfrom
fromtwotwoepochs,
epochs,a adisplacement
displacementvector
vectorisiscomputed.
computed.
Next,
Next,the
the length changeand
length change andthe
the change
change projected
projected to X,
to the the X,direction,
Y, Z Y, Z direction, of the baselines
of the baselines are
are calculated.
calculated. Actually, the start of every 3D vector can be translated to the origin O of the
Actually, the start of every 3D vector can be translated to the origin O of the coordinate system. coordinate
system. Therefore,
Therefore, corresponding
corresponding baselinebaseline vectors
vectors from two from twohave
epochs epochs
thehave
samethe same beginning
beginning (i.e., O
(i.e., O represents
→ → 
A1 and A2 Ain
represents and A2I and
1 epochs in epochs I and II). illustration
II). A graphical A graphicalisillustration is shown
shown in Figure 7b: in
OBFigure 7b: OB1 and
1 and OB2 represent
 → 
OB 2 represent corresponding
corresponding baseline
baseline vectors vectors
in epochs inII,epochs
I and I 1and
thus, B B2 isII,the
thus, B1B2 is the
displacement displacement
vector while θ is
the rotation
vector while angle between
is the rotationthe twobetween
angle baselinethe
vectors;∆X 1 , ∆Y
two baseline 1 and ∆Z
vectors; X1 1are Y1 and
, the Z1 areofthe
coordinates the
baseline vector in epoch I while ∆X2 , ∆Y2 and ∆Z2 are the coordinates of the baseline vector in epoch
coordinates of the baseline vector in epoch I while X 2 , Y2 and Z 2 are the coordinates of the
II. Actually, the baseline change includes the length change and the rotation change. The rotation
baseline
change vector in transformed
is finally epoch II. Actually, the baseline
to the three changeof
axes direction includes the length
the structural changesystem.
coordinate and theTherefore,
rotation
change. The could
the changes rotation change is as
be estimated finally transformed
long as determiningtothe
thedirection
three axes direction
of the of the
axes which structural
is significantly
coordinate system. Therefore,
different from registration. the changes could be estimated as long as determining the direction of
the axes which is significantly different from registration.
2.3. Traditional Change Detection Methods
2.3. Traditional Change Detection Methods
In order to validate the effectiveness of the proposed approach, the baseline method is compared
In order
to two to validate
traditional the effectiveness
methods: direct virtualofpoints
the proposed approach,
comparison the baseline method
and cloud-to-cloud is compared
distances.
to two traditional methods: direct virtual points comparison and cloud-to-cloud distances.
2.3.1. Registration of Two Epochs
2.3.1. Registration of Two
In general, the Epochsof scanning stations are not the same in different epochs, so the
locations
coordinates of identical
In general, pointsofsampled
the locations scanninginstations
consecutive epochs
are not the are
samenotinexpected
differenttoepochs,
be equalsoeither.
the
coordinates of identical points sampled in consecutive epochs are not expected to be equal either. In
traditional change detection methods, two point clouds from two epochs are expected to be in a
common coordinate system. Therefore, registration of two epochs is a vital pre-processing step
required before detecting changes except for those monitoring applications where scan locations as
well as targets positions are fixed by forced centering. Such setup guarantees that the point cloud
Sensors 2017, 17, 26 10 of 25

In traditional change detection methods, two point clouds from two epochs are expected to be in
a common coordinate system. Therefore, registration of two epochs is a vital pre-processing step
Sensors 2017, 17, 26 10 of 25
required before detecting changes except for those monitoring applications where scan locations as
well as targets
acquired at each positions
epoch isare fixed by
referred forced
to the same centering. Such
coordinate setupand
system guarantees
thereforethat the registration.
avoids point cloud
acquired at each epoch is referred to the same coordinate system and therefore avoids
However, it is often difficult in practice to fix both scan locations and targets by forced centering as registration.
However,
this requires it is often difficult
additional in practice
preparation. to fix both
Registration scan and
aligns locations and targets
combines multiple bydata
forced
intocentering as
a single set
this requires additional preparation. Registration aligns and combines multiple
of range data. Registration may introduce small misalignment errors. Even minor misalignments of data into a single set
of range
two data.may
epochs Registration
lead to may introduce
erroneous small
results misalignment
during detecting errors. Even In
changes. minor
this misalignments
case study, the of
two epochs may lead to erroneous
registration is based on control points. results during detecting changes. In this case study, the registration
is based
The on control
Leica points.
Cyclone software provides an automatic registration method, with targets made by
special materials. In thissoftware
The Leica Cyclone seismicprovides an automaticwe
testing experiment, registration method,
set four plane withintargets
targets fixed made
areas byas
special materials. In this seismic
control/tie points for registration. testing experiment, we set four plane targets in fixed areas as
control/tie points for registration.
2.3.2. Virtual Points Extraction and Comparison
2.3.2. Virtual Points Extraction and Comparison
Comparing corresponding points from different epochs [44] is a common way to obtain change
Comparing corresponding points from different epochs [44] is a common way to obtain change
information. Here we consider two ways, to extract virtual points automatically by the procedure
information. Here we consider two ways, to extract virtual points automatically by the procedure
discussed in Section 2.2.2 or manually. After determining virtual points and their correspondence
discussed in Section 2.2.2 or manually. After determining virtual points and their correspondence from
from two epochs, the changes in X, Y and Z direction are computed.
two epochs, the changes in X, Y and Z direction are computed.
2.3.3. Cloud-To-Cloud Distances
2.3.3. Cloud-To-Cloud Distances
It’s generally not so easy to get a clean and proper global model of a surface. Therefore, the idea
It’s generally not so easy to get a clean and proper global model of a surface. Therefore, the idea
of cloud to cloud distance is proposed which is the classical way to detect changes. The principle of
of cloud to cloud distance is proposed which is the classical way to detect changes. The principle of
nearest neighbour distance is used to compute distances between two points: for each point in the
nearest neighbour distance is used to compute distances between two points: for each point in the
compared cloud, the nearest point in the reference cloud is searched and their Euclidean distance is
compared cloud, the nearest point in the reference cloud is searched and their Euclidean distance is
computed [45], see Figure 8.
computed [45], see Figure 8.

Figure8.8.True
Figure Trueglobal
globalmodel
modelfor
forcloud-to-cloud
cloud-to-clouddistances.
distances.

Another
Another immediate
immediate way,
way, getting a better approximation
approximation of of the true distance to the reference
surface,
surface, is to get a local surface model.
model. WhenWhen thethe nearest
nearest point
point inin the
the reference
reference cloud
cloud isis determined,
determined,
the idea is first
first to
to locally
locallymodel
modelthethesurface
surfaceofofthe
thereference
reference cloud by fitting a mathematical
cloud by fitting a mathematical primitive primitiveon
on
thethe ‘nearest’
‘nearest’ point
point andand several
several of itsofneighbours,
its neighbours, see Figure
see Figure 9. The9.distance
The distance
to thistolocal
thismodel
local model is
is finally
finally reported.
reported. Common Common
ways to ways
locallytomodel
locally model are
a surface a surface are bycompare
by triangles, triangles, compare
Figure 9, byFigure
planes 9,or by
by
planes or by
otherwise otherwise
smooth smooth
patches. patches. The of
The effectiveness effectiveness
this methodofisthis method is
statistically statistically
more more or less
or less dependent on
dependent on the cloud
the cloud sampling sampling
and on and on how
how appropriate theappropriate theapproximation
local surface local surface approximation
is. is.
Sensors 2017, 17, 26 11 of 25
Sensors 2017, 17, 26 11 of 25

Figure 9. Local surface model for cloud-to-cloud distances.

Figure 9. Local surface model for cloud-to-cloud distances.

To
Toprofit
profitfrom
fromobservation
observationredundancy,
redundancy,least
leastsquare
squarefitting
fittingisisused
usedtotofitfitaaplane
planeor
oraahigher
higherorder
order
model to the points of the reference cloud.
model to the points of the reference cloud.
3. Results
3. Results
The results of brick centre extraction are presented below. Next, the resulting virtual points are
The results of brick centre extraction are presented below. Next, the resulting virtual points are
used to establish baselines. Then, differences in corresponding baselines are presented as change
used to establish baselines. Then, differences in corresponding baselines are presented as change
detection results and compared to direct virtual point comparison and cloud-to-cloud distances results.
detection results and compared to direct virtual point comparison and cloud-to-cloud distances
results.
3.1. Structural Coordinate System Establishment and Coordinate Transformation
In order toCoordinate
3.1. Structural facilitate interpretation of the results,
System Establishment we introduce
and Coordinate a so-called structural coordinate
Transformation
system. Within this coordinate system, obtained 3D change vectors are roughly decomposed into
In orderand
wall-parallel to facilitate interpretation
wall-perpendicular of the results,
change. As shownwe introduce
in Figure a2,so-called structural
the points coordinate
of the stable wall
system. Within this coordinate system, obtained 3D change vectors are roughly
were selected to establish a structural coordinate system according to the principle introduced in decomposed
into wall-parallel
Section and wall-perpendicular
2.2.1. The plane equations in epochchange. As shown
I and epoch in follows:
II are as Figure 2, the points of the stable wall
were selected to establish a structural coordinate system according to the principle introduced in
Section 2.2.1. The plane Epoch
equations in epoch
I : 0.9348x + I0.3552y
and epoch II are −
+ 0.0015z as 8.8534
follows:
= 0;
Epoch I: 0.9348 x  0.3552 y  0.0015 z  8.8534  0 ;
Epoch II : 0.9974x − 0.0727y − 0.0002z − 8.8518 = 0;

The translation vectorEpoch

and the 0.9974 x angle
II: rotation 0.0727are 0.0002 z  8.8518
y calculated  0 TLS
between ; and reference plane.
0
Because the scanner is always level, the Z axis is the same as the Z axis, and the normal vector of the
The translation vector and the rotation angle are calculated between TLS and reference plane.
stable wall is regarded as the positive the X 0 axis. Therefore, the 3D coordinate transformation can be
Because the scanner is always level, the Z axis is the same as the Z ' axis, and the normal vector
defined by:
of the stable wall is regarded as the positive the X ' axis. Therefore, the 3D coordinate
transformation can be defined X 0 by:
      
0.9348 −0.3552 0 X −8.2760
Epoch I :  Y 0  =  0.3552 0.0094 0  Y  +  −3.1451 ;
      
Z0  '

X 0.0000 0.9348 0.0000
 0.35521 0  ZX   8.2760 
−0.0132
 '     
 I: Y  0.3552
Epoch 0.0094 0 Y   3.1451 ; 
X  '   0.0000
0 0.9974 0.07270.00000 1  ZX   0.0132 
−8.8284
 0  Z       
Epoch II :  Y  =  −0.0727 0.0094 0  Y  +  0.6432 ;

Z0 0.0000 0.0000 1 Z 0.0020
 X '   0.9974 0.0727 0   X   8.8284 
Finally, point clouds from  ' epochs    the TLS coordinate

Epoch II: both are transformed from
Y     0.0727 0.0094 0 Y    0.6432 ;
system to the
 Z '   0.0000 0.0000 1  Z   0.0020
structural coordinate system, respectively.
      

Finally, point clouds from both epochs are transformed from the TLS coordinate system to the
structural coordinate system, respectively.
Sensors 2017, 17, 26 12 of 25

Sensors 2017,
3.2. Virtual Point17, Extraction
26 Results 12 of 25

3.2.1.3.2. Virtual
Automatic
Sensors
Point Extraction Results
2017, 17,Method
26 12 of 25
3.2.1.results
The Automatic Method virtual points automatically are shown in this section. Here, we take
of extracting
3.2. Virtual Point Extraction Results
patch A as Theexample
results to
of illustrate
extracting the method.
virtual points automatically are shown in this section. Here, we take
3.2.1. Automatic Method
patch
1. A as
First, theexample to illustrate
brick and the method.
mortar points are separated by k-means clustering based on their intensity.
1.
TheFirst, the
resultsbrick
of and mortar
extracting
Figure 10a shows the extracted mortar points points
virtual are separated
points epoch Iby
automatically
in k-means
areFigure
while shownclustering
in this
10b based onHere,
section.
shows both their
brickintensity.
we take
and mortar
Figure
patch A10aasshows
examplethe extracted
to illustratemortar
the points
method. in epoch I while Figure 10b shows both
points in different colors after separation in epoch II. Other researchers, using photogrammetry, pasted brick and mortar
smallpoints
black 1.inFirst,
different
targets at colors
the brick
the and after separation
mortar
centre in epoch
points are
of each brick. As aII.result,
separated Other researchers,
by k-means
brick using based
clustering
centre pointsphotogrammetry,
to havepasted
on their intensity.
appear the same
Figure
small 10a shows
black targetsthe extracted
at the centremortar
of eachpoints
brick.inAsepoch I while
a result, Figure
brick 10b
centre shows
points both brick
appear andthe
to have mortar
same
color as mortar points. However, we estimated the brick centre locations from lines passing through
points
color asinmortar
different colorsHowever,
points. after separation in epochthe
we estimated II. Other researchers,
brick centre using
locations photogrammetry,
from pasted
lines passing through
mortar point.
small
mortar blackTherefore,
point.targets
thethe
at the
Therefore,
appearance
centre of each of
appearance
the brick
brick.
of theAsbrick
centres
a result, brick
centres
will notpoints
centre
will
affectour
not affect
our
appear results.
to have the same
results.
color as mortar points. However, we estimated the brick centre locations from lines passing through
mortar point. Therefore, the appearance of the brick centres will not affect our results.

(a) Epoch I (b) Epoch II

Figure 10. Results of (a)
separating
Epoch I brick points from notably mortar (b) points in patch A using k-means
Figure 10. Results of separating brick points from notably mortar points Epoch
inII patch A using k-means
clustering, with k = 2. Blue points correspond to mortar points, point reflecting from small
clustering,
Figure with k = of
10. Results 2. separating
Blue points
brick correspond to mortar
points from notably points,
mortar point
points in patchreflecting from small
A using k-means
photogrammetric targets attached to the brick centres, or to some wires in front of the wall. The yellow
clustering, with
photogrammetric k = attached
targets 2. Blue points
to the correspond
brick toormortar
centres, to points,
some wires point
in reflecting
front of the from The
wall. small
yellow
points in general belong to brick surfaces.
pointsphotogrammetric
in general belongtargets attached
to brick to the brick centres, or to some wires in front of the wall. The yellow
surfaces.
points in general belong to brick surfaces.
2. Resulting mortar points are projected to the Z and Y axis respectively. Through in situ
measuring
2. by steel
Resulting ruler,
mortar the widthprojected
points of the mortar is Zestimated at Lmortar  0.012m Therefore, the
2. Resulting mortar pointsare to the
are projected to the Z andandYYaxis
axis respectively.
respectively. Through
Through in situ
in situ
window
measuring width
by steel and step
ruler, width
the the
width are defined
of of
thethe as
mortar L
mortaris  0.008m
isestimated
estimated at and L
at LLmortar  0.002m based on the
step = 0.012 m. Therefore, the
measuring by steel ruler, width mortar  0.012m Therefore, the
window

window width
equation
window Lmortar
andand
width Lwindow
step  2 width
width
step Lstepare defined
introduced
are as Lwindow
in Section
defined as =
2.2.4.
Lwindow 0.008 mand
andLthe
Afterwards,
 0.008m L = 0.002
brick
stepstep
m based
boundary
0.002m based lines
on the on the
are
equation Lmortar
estimated
equation by ≈
Lmortar Lwindow +22L
theLprinciple
window introduced
Lstepstep introduced
in Section
introduced in 2.2.4.
Section
in Section 2.2.4.2.2.4. Afterwards,
Afterwards, theboundary
the brick brick boundary
lines are lines
are estimated
estimatedby
bythe
the principle introduced
principle introduced in in Section
Section 2.2.4.
2.2.4.

(a) Epoch I Y axis (b) Epoch II Z axis

(a) Epoch I Y axis (b) Epoch II Z axis
Figure 11. Histogram of mortar point frequencies along the Z and the Y axes.
Figure
Figure 11. 11. Histogram
Histogram ofofmortar
mortarpoint
point frequencies
frequenciesalong thethe
along Z and the the
Z and Y axes.
Y axes.
Sensors 2017, 17, 26 13 of 25
Sensors 2017, 17, 26 13 of 25

The histograms of the number of points along Z and Y axes are shown in Figure 11. Figure 11a
presents the histogram
The histograms of the
of the number
number of points
of points along
along the YY axes
Z and axis are
in epoch
shownI while Figure
in Figure 11b shows
11. Figure 11a
the histogram
presents of the number
the histogram of pointsofalong
of the number the
points Z axis
along theinYepoch
axis inII.epoch I while Figure 11b shows the
3. Finally,
histogram of thethe 3D coordinates
number of allthe
of points along brick centres
Z axis are estimated
in epoch II. by implementing the method
proposed in Section
3. Finally, 2.2.2.
the 3D Patch A contains
coordinates 16 bricks;
of all brick centresPatch B containsby
are estimated three bricks; Patchthe
implementing C contains
method
15 bricks; in
proposed Patch D contains
Section six bricks.
2.2.2. Patch Figure
A contains 16 12 shows
bricks; the B
Patch extraction
contains results of patch
three bricks; A in
Patch epochs I
C contains
and
15 II. Patch D contains six bricks. Figure 12 shows the extraction results of patch A in epochs I and II.
bricks;

(a) Epoch I (b) Epoch II

Figure 12. Virtual point extraction results (patch A). The red dots are the final virtual points.
Figure 12. Virtual point extraction results (patch A). The red dots are the final virtual points.

3.2.2. Comparison to Manual Virtual Point Extraction

3.2.2. Comparison to Manual Virtual Point Extraction
The automatic extraction is evaluated against manual picking the centre points of each brick. A
The automatic extraction is evaluated against manual picking the centre points of each brick.
human operator imported the two aligned point clouds to the “Cloudcompare” software and picked
A human operator imported the two aligned point clouds to the “Cloudcompare” software and picked
points representing the brick centres. By considering the study cases patch A, B, C and D, which
points representing the brick centres. By considering the study cases patch A, B, C and D, which
contain a total of 40 bricks, in the moving wall, the average difference between manually picking and
contain a total of 40 bricks, in the moving wall, the average difference between manually picking and
automatic extracting are [−0.23,1.09,−1.27] mm in the X, Y and Z direction respectively. The standard
automatic extracting are [−0.23,1.09,−1.27] mm in the X, Y and Z direction respectively. The standard
deviation along the X, Y and Z direction is 0.92, 2.48 and 1.44 mm, respectively. This gives an
deviation along the X, Y and Z direction is 0.92, 2.48 and 1.44 mm, respectively. This gives an indication
indication for the effectiveness and feasibility of the proposed automatic extraction method.
for the effectiveness and feasibility of the proposed automatic extraction method.
3.3. Baseline
3.3. Baseline Establishment
Establishment andand Decomposing
Decomposing
To establish
To establish the
the baselines,
baselines, we
we first
first select
select the
the four
four patches
patches A,
A, B,
B, C
C and
and D.D. In
In addition,
addition, target points
target points
are identified. Through either linking different virtual brick centre points as extracted
are identified. Through either linking different virtual brick centre points as extracted from the patches, from the
patches,
or linkingorvirtual
linking virtual
brick centrebrick centre
points pointspoints,
to target to target points,are
baselines baselines are established.
established. Figure 18
Figure 13 contains 13
contains 18 such baselines. The 12 baselines shown as continuous lines, link brick centres
such baselines. The 12 baselines shown as continuous lines, link brick centres in patch A, B, C and D to in patch A,
B, Cthree
the and targets
D to theTa3,
three
Ta6targets Ta3,
and Ta7. Ta6
The sixand Ta7. The
baselines six baselines
shown as dashed shown
lines,as dashed
link lines, link
brick centres brick
to other
centres to other
brick centres. brick centres.
By comparing
By comparing corresponding
correspondingbaselines
baselinesinintwo twoepochs,
epochs,change
change vectors
vectorsareare
estimated.
estimated.Then, the
Then,
change
the vectors
change are are
vectors projected to the
projected X, YX,and
to the Z axes
Y and within
Z axes the the
within structural
structuralcoordinate system.
coordinate system.
SensorsSensors 17, 2617, 26
2017, 2017, 14 of 2514 of 25

Figure
Figure 13. 13.Baselines
Baselines from
from virtual
virtualpoints to target
points pointspoints
to target are indicated by continuous
are indicated lines, while lines,
by continuous
baselines connecting virtual points are shown as dashed lines.
while baselines connecting virtual points are shown as dashed lines.

3.4. Baseline Changes

3.4. Baseline Changes
3.4.1. Baseline Changes between Patches and Targets
3.4.1. Baseline Changes between Patches and Targets
Figure 14 shows the length change for a total of 120 pairs of corresponding baselines between
Figurepoints
target 14 shows the length
and virtual brick change for a total
centre points, of 120
in which pairs
“Avg” of corresponding
represents baselines
the mean baseline between
length
target points and virtual brick centre points, in which “Avg” represents
change and “SD” represents the standard deviation of the baseline length change. Here mean andthe mean baseline length
change and “SD”
standard represents
deviation the
are taken standard
from all pairsdeviation of the baseline
of corresponding baselines length change.
connecting Here
the same mean and
patches
standard
and/or deviation
targets. Thearemaximum
taken from all change,
length pairs offrom
corresponding
target “Ta6” tobaselines
the virtual connecting
brick centrethe same
point “A16”patches
in patch
and/or A, is
targets. The−67.79 mm in which
maximum lengththe minus from
change, sign indicates that the
target “Ta6” length
to the of baseline
virtual “Ta6-A16”
brick centre pointin“A16”
epoch
in patch A,IIisbecomes
−67.79shorter
mm incompared
which the to epoch
minusI.sign
All the baselinesthat
indicates fromthe“Ta6” and of
length “Ta7” to virtual
baseline brick
“Ta6-A16” in
centre points became shorter. This indicates that the moving wall as a whole tilted
epoch II becomes shorter compared to epoch I. All the baselines from “Ta6” and “Ta7” to virtual brick to the left during
the experiment. The baselines from “Ta3” to virtual brick centre points in patch A and patch C
centre points became shorter. This indicates that the moving wall as a whole tilted to the left during
become longer, while baselines from “Ta3” to virtual brick centre points in patch B and patch D
the experiment. The baselines from “Ta3” to virtual brick centre points in patch A and patch C become
became shorter. This contrasting pattern could match the previous pattern of the moving building
longer, while
tilting baselines
to the from
left: patch “Ta3”
A and to virtual
patch C movebrick
awaycentre pointswhile
from “Ta3”, in patch B and
patches patch
B and D became
D move towards shorter.
This “Ta3”.
contrasting pattern could match the previous pattern of the moving building tilting to the left:
patch A and For patch C move away the
better understanding from “Ta3”,inwhile
changes patches
different areas,B we
andanalyzed
D movethe towards
changes“Ta3”.
per group,
For better understanding
summarized in Table 1. As shownthe changes in different
in Figure 14, areas, we
the mean changes from analyzed the changes
patch A, patch B, patch Cper
and group,
patch D to “Ta6”, are −66.35, −31.08, −31.86 and −4.96 mm, respectively. The
summarized in Table 1. As shown in Figure 14, the mean changes from patch A, patch B, patch C mean changes from patch
A, patch
and patch D B,topatch
“Ta6”, C and
are patch
−66.35,D to−“Ta7”,
31.08, are −58.65,
−31.86 −30.07,
and −4.96−37.87
mm, and −6.58 mm, respectively.
respectively. The mean The changes
results indicate that the four patches moved towards “Ta6” and “Ta7”. However, the absolute values
from patch A, patch B, patch C and patch D to “Ta7”, are −58.65, −30.07, −37.87 and −6.58 mm,
are quite different and show that patch A, B and C have a more obvious tendency towards “Ta6” and
respectively. The results indicate that the four patches moved towards “Ta6” and “Ta7”. However,
“Ta7” than patch D.
the absolute values are quite different and show that patch A, B and C have a more obvious tendency
towards “Ta6” and “Ta7” than patch D.
Sensors 2017,
Sensors 2017, 17,
17, 26
26 15 of
15 of 25
25

Figure 14. Baseline length changes (from target points to virtual brick centre points). The quantitative
Figure
results 14.
are Baseline length
summarized in changes
Table 1. (from target points to virtual brick centre points). The quantitative
results are summarized in Table 1.

The mean change from patch A, patch B, patch C and patch D to “Ta3”, are 54.11 mm, −9.57 mm,
The mean change from patch A, patch B, patch C and patch D to “Ta3”, are 54.11 mm, −9.57 mm,
29.25 mm and 0.99 mm, respectively. “Ta3” is at the middle of the four patches along the X axis which
29.25 mm and 0.99 mm, , respectively. “Ta3” is at the middle of the four patches along the X axis
lead to a different change compared to “Ta6” and “Ta7”: patch A and patch C moved away from
which lead to a different change compared to “Ta6” and “Ta7”: patch A and patch C moved away
the “Ta3” after the seismic testing while patch B moved towards “Ta3”. Patch D has a small change
from the “Ta3” after the seismic testing while patch B moved towards “Ta3”. Patch D has a small
w.r.t. “Ta3” compared to the standard deviation of the change. Therefore we cannot draw a significant
change w.r.t. “Ta3” compared to the standard deviation of the change. Therefore we cannot draw a
conclusion on its movement relative to “Ta3”.
significant conclusion on its movement relative to “Ta3”.
Table 1. Baseline change between patches and targets.
Table 1. Baseline change between patches and targets.
Maximum Minimum Mean Change Standard Deviation
Target Patch Maximum Minimum Mean Change (mm)Standard Deviation
Target Patch Change (mm) Change (mm) (mm)
Change (mm) Change (mm) (mm) (mm)
A 57.86 50.70 54.11 2.05
A B −57.86
10.09 50.70
−8.98 −9.57 54.11 0.56 2.05
Ta3
B C 41.36
−10.09 17.04
−8.98 29.25 −9.57 7.32 0.56
Ta3 D −4.81 −1.10 0.99 2.92
C 41.36 17.04 29.25 7.32
D A −−4.81
67.79 −64.14
−1.10 −66.35 0.99 1.07 2.92
B −31.28 −30.95 −31.08 0.17
Ta6 A C −−67.79
41.32 −64.14
−21.02 −31.86 −66.35 6.14 1.07
B D −
−31.28
7.42 −3.32
−30.95 −4.96 −31.08 1.64 0.17
Ta6
C A −−41.32
59.88 −21.02
−56.48 −58.65 −31.86 1.00 6.14
D B −−7.42
30.77 −30.46
−3.32 −30.57 −4.96 0.17 1.64
Ta7
C −51.81 −22.18 −37.87 9.10
A D −−59.88
10.54 −56.48
−3.79 −6.58 −58.65 2.72 1.00
B −30.77 −30.46 −30.57 0.17
Ta7
C −51.81 −22.18 −37.87 9.10
3.4.2. Baseline Changes within a Patch
D −10.54 −3.79 −6.58 2.72
Obviously, the standard deviations of the mean change involving patch C are relatively large
which indicate Changes
3.4.2. Baseline that damage happened
within a Patchwithin patch C. To validate this hypothesis, we select six virtual
brick centre points, three in the top (C1, C2 and C3) and the other three in the bottom (C13, C14 and
Obviously, the standard deviations of the mean change involving patch C are relatively large
C15) of patch C. By linking virtual points in the top area to virtual points in the bottom area, 9 baselines
which indicate that damage happened within patch C. To validate this hypothesis, we select six
for patch C are established and decomposed. Here the baseline direction is taken from the bottom area
virtual brick centre points, three in the top (C1, C2 and C3) and the other three in the bottom (C13,
to the top area, which means the bottom area is taken as reference. The mean change of the baselines
C14 and C15) of patch C. By linking virtual points in the top area to virtual points in the bottom area,
in patch C is −35.1, 0.0 and 0.6 mm in the X, Y and Z direction, respectively, in which the minus
9 baselines for patch C are established and decomposed. Here the baseline direction is taken from the
sign indicates that the change in the X direction is pointing into the moving building. Their standard
bottom area to the top area, which means the bottom area is taken as reference. The mean change of
the baselines in patch C is −35.1, 0.0 and 0.6 mm in the X, Y and Z direction, respectively, in which
Sensors 2017, 17, 26 16 of 25
Sensors 2017, 17, 26 16 of 25
the minus sign indicates that the change in the X direction is pointing into the moving building. Their
standard deviations are 4.7, 0.9 and 0.1 mm, respectively. Because of the big change in the X direction,
deviations
Figure 15 isare 4.7, 0.9which
plotted, and 0.1 mm, the
shows respectively. Because
baseline (in patch of
C) the big change
change in the Xindirection
the X direction, Figure
for a total 15
of nine
is plotted, which shows the baseline (in patch C) change in the X direction for a total of nine
pairs of corresponding baselines between different brick centre points. It is concluded that damage pairs of
corresponding
happened to patchbaselines between
C in the different brick centre points. It is concluded that damage happened
X direction.
to patch C in the X direction.

Figure 15. Baselines change in the X direction in patch C.

Figure 15. Baselines change in the X direction in patch C.

3.4.3.
3.4.3. Baseline
Baseline Changes
Changes between
between Patches
Patches
The
The variations
variationsinin baselines
baselinesconnecting
connectingtarget points
target to virtual
points brick centre
to virtual brick points
centreindicate
points absolute
indicate
change. While variation in baselines between different virtual points indicate
absolute change. While variation in baselines between different virtual points indicate relative relative change from one
change
brick
from to
onethe other.
brick Figure
to the 16 Figure
other. shows the lengththe
16 shows change
lengthforchange
a total for
of 537 pairs
a total of of
537corresponding baselines
pairs of corresponding
between different virtual brick centre points, in which “Avg” represents
baselines between different virtual brick centre points, in which “Avg” represents the mean the mean corresponding
baselines length
corresponding change length
baselines and “SD” represents
change and “SD” the standardthe
represents deviation
standardofdeviation
corresponding baselines
of corresponding
length changes.
baselines length changes.
The
The change
change results
results are
are also
also summarized
summarized in in Table
Table 2.2. We
We selected
selected the
the baselines
baselines from
from virtual
virtual brick
brick
centre points in patch B, C and D to virtual brick centre points in patch A, as examples
centre points in patch B, C and D to virtual brick centre points in patch A, as examples to illustrate to illustrate the
change pattern.
the change Patch
pattern. D shows
Patch maximal
D shows maximalchange
changewithwith
regard to patch
regard A. ItsA.
to patch mean change
Its mean is 53.94
change mm.
is 53.94
Patch C shows minimal change with regard to patch A. The mean change
mm. Patch C shows minimal change with regard to patch A. The mean change is 2.14 mm with a is 2.14 mm with a standard
deviation of 2.20 mm.
standard deviation The mm.
of 2.20 meanThechange
meanofchange
patch Bofwithpatch regard
B withtoregard
patch A to is 44.37Amm
patch withmm
is 44.37 a much
with
smaller standard deviation. The results indicate that patch B and D moved
a much smaller standard deviation. The results indicate that patch B and D moved away from away from patch Apatch
after
the seismic
A after testing. testing.
the seismic
3.4.4. Displacement Vectors
3.4.4. Displacement Vectors
For easier interpretation of the detected changes, difference vectors between corresponding
For easier interpretation of the detected changes, difference vectors between corresponding
baselines are plotted relative to an arbitrary reference point “A1”. Figure 17a shows the position of
baselines are plotted relative to an arbitrary reference point “A1”. Figure 17a shows the position of
point “A1” and the other three patches. The difference vectors from the perspective of X axis between
point “A1” and the other three patches. The difference vectors from the perspective of X axis between
corresponding baselines incident to “A1” are plotted in Figure 17b–d. However, these figures do not
corresponding baselines incident to “A1” are plotted in Figure 17b–d. However, these figures do not
illustrate possible displacements in the X direction. Therefore, additional difference vectors in 3D
illustrate possible displacements in the X direction. Therefore, additional difference vectors in 3D
between corresponding baselines incident to “A1” are plotted in Figure 17e–g.
between corresponding baselines incident to “A1” are plotted in Figure 17e–g.
Sensors 2017, 17, 26 17 of 25
Sensors 2017, 17, 26 17 of 25

Figure 16. Baseline

Figure 16. lengthchanges
Baseline length changesbetween
betweenvirtual
virtual brick
brick centre
centre points.
points.

Table 2.
Table Baselinechange
2. Baseline changebetween
betweenpatches.
patches.

Maximum Change
Maximum Minimum Change
Minimum Mean
Mean Change Standard
Change Standard Deviation
Deviation
Patch Patch
Patch Patch(mm)
Change (mm) (mm)
Change (mm) (mm)
(mm) (mm)(mm)
B B 44.63 44.63 44.07
44.07 44.37
44.37 0.17 0.17
A CA C 6.93 6.93 0.06
0.06 2.14
2.14 2.20 2.20
D D 56.50 56.50 50.87
50.87 53.94
53.94 1.21 1.21
C C −9.60 −9.60 −7.66
−7.66 −8.56
−8.56 0.52 0.52
B B
D D 3.61 3.61 2.05
2.05 2.93
2.93 0.58 0.58
C DC D 15.49 15.49 13.33
13.33 14.56
14.56 0.64 0.64

By comparing the baselines projected to the X, Y and Z axes within the structural coordinate
By comparing the baselines projected to the X, Y and Z axes within the structural coordinate
system, relative change information is obtained. Change is considered relative if it corresponds to
system, relative change information is obtained. Change is considered relative if it corresponds to
change between different areas in the moving building, and therefore the change value here is
change between different areas in the moving building, and therefore the change value here is absolute
absolute not relative. As regards the structural coordinate system, we have checked the difference in
not relative. As regards the structural coordinate system, we have checked the difference in two epochs
two epochs which is little enough to be ignored, namely this will not affect the results of our proposed
which is little
methods. Forenough
easier toto be ignored,
implement, wenamely
projectedthis
thewill not affect
baselines theepochs
in two resultsrespectively.
of our proposed methods.
Full change
For easier to implement,
information we projected
is given in Table 3. Here wethetake
baselines in two epochs
the baselines respectively.
connected Fullpoints
to the virtual changeininformation
patch A
isas
given
examples. The maximum change direction is along the Y axis in the structural coordinateassystem.
in Table 3. Here we take the baselines connected to the virtual points in patch A examples.
The
Themaximum
mean change change
valuesdirection is along
are −44.38, −58.97the
andY−73.33
axis inmmthefor
structural
patch B, coordinate system.inThe
C, D, respectively, whichmean
change values
the minus signare −44.38,
denotes the− 58.97 and
direction in − 73.33
the mm for
opposite patch
Y axis B, C,
or left. D, respectively,
Corresponding in which
standard the minus
deviations
sign
are denotes the direction
all well below in the opposite
1 mm. Along the X axis, Y the
axismean
or left. Corresponding
change has a valuestandard deviations
of 2.82, −24.45 are mm
and 5.06 all well
below 1 mm.
for patch B, Along the X axis, the
C, D, respectively. Themeanmeanchange
changehasvalue
a value of 2.82,mm
of −24.45 −24.45 and 5.06
in patch mm for patch
C indicates that, B,
C,compared to otherThe
D, respectively. patches,
meanvirtual
change brick
value of −points
centre in patch
24.45 mm C have
in patch a relatively
C indicates large
that, change along
compared to other
the X axis,
patches, which
virtual is pointing
brick into the
centre points moving
in patch building.
C have Compared
a relatively largetochange
the other two
along directions,
the the is
X axis, which
pointing into the moving building. Compared to the other two directions, the mean change along the
Sensors 2017, 17, 26 18 of 25

Sensors 2017, 17, 26 18 of 25

Z axis is less significant: Its values are −0.39, 2.45, 2.82 mm for patch B, C, D, respectively with in all
mean change along the Z axis is less significant: Its values are −0.39, 2.45, 2.82 mm for patch B, C, D,
three cases a standard deviation of ~0.8 mm. The above structural change results correspond well to
respectively with in all three cases a standard deviation of ~0.8 mm. The above structural change
the emerging cracks/damage shown in Figure 3.
results correspond well to the emerging cracks/damage shown in Figure 3.

(a) Positions of patches and

(b) Patch B (c) Patch C (d) Patch D
reference point “A1”

(e) Patch B in 3D (f) Patch C in 3D

(g) Patch D in 3D

Figure 17. Displacement vectors. In the Figures (b–d) 2D vectors are given, indicating relative changes
Figure 17. Displacement vectors. In the Figures (b–d) 2D vectors are given, indicating relative changes
parallel to the affected wall. Full 3D vectors are given in Figures (e–g).
parallel to the affected wall. Full 3D vectors are given in Figures (e–g).
Table 3. Displacement of baselines connected to patch A.
Table 3. Displacement of baselines connected to patch A.
Maximum Change Minimum Change Mean Change Standard Deviation
Patch Patch Axis
(mm) (mm) (mm) (mm)
Maximum Minimum Mean Change Standard
Patch X
Patch Axis 4.40 0.56 2.82 1.00
Change (mm) Change (mm) (mm) Deviation (mm)
B Y −44.60 −44.21 −44.38 0.16
Z X 1.61 4.40 0.56
0.01 2.82
−0.39 1.000.81
BX Y −46.36 −44.60 −44.21
−1.33 −−24.45
44.38 0.1612.60
Z 1.61 0.01 −0.39 0.81
A C Y −60.00 −58.23 −58.97 0.61
Z X 4.77 − 46.36 −
1.381.33 −24.45
2.45 12.600.83
A Y − 60.00 −58.23 −58.97 0.614.98
CX 12.11 0.07 5.06
D Y Z −73.71 4.77 1.38
−73.00 2.45
−73.33 0.830.23
Z X 5.24 12.11 0.07
1.67 5.06
2.82 4.980.87
D Y −73.71 −73.00 −73.33 0.23
Z 5.24 1.67 2.82 0.87
3.5. Comparison to Traditional Change Detection Methods

3.5.1. Registration
3.5. Comparison of Two Change
to Traditional Epochs Detection Methods
The control point registration is performed by the Leica Cyclone software. As introduced before,
3.5.1. Registration of Two Epochs
four plane targets were set. Control/tie points were additionally measured by the scanner, to
reference
The thepoint
control scans registration
more accurately to a common
is performed bycoordinate
the Leicasystem. However,
Cyclone theAs
software. registration
introduced error
before,
of “Tag 5” was 0.010 m which was relatively high. Therefore, we have reason to assume that “Tag
four plane targets were set. Control/tie points were additionally measured by the scanner, to reference 5”
the scans more accurately to a common coordinate system. However, the registration error of “Tag
5” was 0.010 m which was relatively high. Therefore, we have reason to assume that “Tag 5” was
Sensors 2017, 17, 26 19 of 25

apparently moved
Sensors 2017, 17, 26 by an external agent during the experiment. Considering that the scanner is 19 always
of 25
level, three targets are enough for registration. Thus, using the remaining three targets, the two point
was are
clouds apparently
alignedmoved
whichby an to
lead external agent
a better during
result. Thethe experiment.errors
registration Considering that theare
of the targets scanner
shown is in
always
Table 4. level, three targets are enough for registration. Thus, using the remaining three targets, the
two point clouds are aligned which lead to a better result. The registration errors of the targets are
shown in Table 4. Table 4. Automatic registration error by Leica Cyclone software.

TableName
Target 4. Automatic registrationError
Weight error(m)
by Leica Cyclone software.
Error Vector (m)
Ta3
Target Name 1
Weight 0.001Error (m)
(0.001,0.000,0.000)
Error Vector (m)
(−0.001,0.000,0.000)
Ta3 Ta6 1
1 0.000
0.001 (0.001,0.000,0.000)
Ta7 1 0.001 (0.000,0.000,0.000)
Ta6 1 0.000 (−0.001,0.000,0.000)
Ta7 1 0.001 (0.000,0.000,0.000)
3.5.2. Virtual Point Comparison
3.5.2.
AfterVirtual Point Comparison
registration of the two epochs of data, all virtual brick centre points in each patch were
extracted. Weregistration
After compared theof the3Dtwo
coordinates
epochs ofofdata,
the corresponding
all virtual brickvirtual
centre points in
from two
each epochs
patch wereand
extracted.the
determined Wechange
compared
in X,theY3D
andcoordinates of the
Z direction. corresponding
Figure 18a–d show virtual points arrows
the change from two forepochs and
four patches
determined the change in X, Y and Z direction. Figure 18a–d show the change arrows
projected on the YOZ plane. In patches C and D, the length of the change arrows is shorter than in for four patches
projected
patches onBthe
A and YOZindicates
which plane. In that
patches C and D,in
the changes theY length
and Z of the change
direction arrowscompared
are small is shorter to
than in C
patch
andpatches A and
patch D. B which
In order indicates
to reveal that the
changes changes
in the in Y andchange
X direction, Z direction arefor
arrows small
fourcompared to patch
patches are plotted
in 3D in Figure 18e–h. From the above figures, we draw the conclusion that patch C moved inare
C and patch D. In order to reveal changes in the X direction, change arrows for four patches the X
plotted in 3D in Figure 18e–h. From the above figures, we draw the conclusion that patch C moved
direction. Furthermore, there have a big difference of every virtual points in patch C which prove that
in the X direction. Furthermore, there have a big difference of every virtual points in patch C which
damage happened in patch C not just shift movement.
prove that damage happened in patch C not just shift movement.

(a) Patch A (b) Patch B (c) Patch C (d) Patch D

(e) Patch A in 3D (f) Patch B in 3D

(g) Patch C in 3D (h) Patch D in 3D

Figure 18. Direct virtual point comparison results.

Figure 18. Direct virtual point comparison results.
Sensors 2017, 17, 26 20 of 25

In the top area (patch A and patch B), the maximum change direction is along the Y axis which
is horizontal and parallel to the moving wall. While in the bottom area (patch C and patch D), the
maximum change direction is in the X direction which is pointing into the moving wall. Along the
Z direction (vertical direction), all four patches have a similar tendency: the change is marginal with
a value of less than 5 mm. The change at the top area is slightly larger than at the bottom area.
Sensors 2017, 17, 26 20 of 25
This means that after the seismic testing, no obvious vertical change of the moving wall occurred.
In the Y direction, In the top area (patch A and patch B), the maximum change direction is along the Y axis which the bottom
the top area (patch A and patch B) has a significantly greater change than
is horizontal
area: the absolute displacement and parallel tovalues
the moving wall.four
of the Whilepatches
in the bottom
(A, area
B, C(patch
and CD)and inpatch
the D), the
Y direction are in
maximum change direction is in the X direction which is pointing into the moving wall. Along the Z
[75.75 76.98] direction
mm, [31.55 31.55] mm, [17.75 18.13] mm and [1.80 2.49] mm, respectively
(vertical direction), all four patches have a similar tendency: the change is marginal with a
(In the symbol [],
the left valuevalue
is the minimum value in a patch while the right value is the maximum
of less than 5 mm. The change at the top area is slightly larger than at the bottom area. This value in the same
means that
patch). Compared withafter∆Y, ∆Z, the
the seismic testing,
pattern of ∆X
no obvious vertical
has achange of the moving The
big difference. wall occurred.
changes In in
the patch A and
Y direction, the top area (patch A and patch B) has a significantly greater change than the bottom
patch B are smooth which indicate the patches shift as a whole, however, in patch C and patch D the
area: the absolute displacement values of the four patches (A, B, C and D) in the Y direction are in
changes vary[75.75
prominently which
76.98] mm, [31.55 indicates
31.55] that
mm, [17.75 apparent
18.13] damage
mm and [1.80 happened
2.49] mm, within
respectively patch C and patch
(In the symbol
[], the left value is the minimum value in a patch while the right value
D. The above virtual points comparison results provide information similar to the proposed is the maximum value in the baseline
same patch). Compared with ΔY, ΔZ, the pattern of ΔX has a big difference. The changes in patch A
method and and meanwhile verifies its effectiveness.
patch B are smooth which indicate the patches shift as a whole, however, in patch C and patch
D the changes vary prominently which indicates that apparent damage happened within patch C
3.5.3. Cloud-to-Cloud
and patch D.Distances
The above virtual points comparison results provide information similar to the
proposed baseline method and meanwhile verifies its effectiveness.
First, we plot the cloud-to-cloud distances for the entire moving wall, see Figure 19. As the wall is
approximately 3.5.3. Cloud-to-Cloud
planar, a LeastDistances
Squares best fitting plane is used while searching KNN points, where
each leaf node contains a number of predefined
First, we plot the cloud-to-cloud distances target points.
for the entire In this
moving wall, research, six
see Figure 19. Asnearby
the wall points were
is approximately planar, a Least Squares best fitting plane is used while searching KNN points, where
used. The maximum cloud-to-cloud distance is 0.0846 m in the bottom-left area. Actually, we can
each leaf node contains a number of predefined target points. In this research, six nearby points were
easily draw the
used.conclusion
The maximumfrom Figure 19
cloud-to-cloud that is
distance just using
0.0846 m incloud-to-cloud
the bottom-left area.distances
Actually, wegives
can incomplete
information on those
easily drawareas that got
the conclusion internally
from damaged.
Figure 19 that just usingIndeed, changes
cloud-to-cloud parallel
distances to a wall are difficult
gives incomplete
information on those areas that got internally damaged. Indeed, changes parallel to a wall are difficult
to identify. Figure 19 heights change parallel to the wall.
to identify. Figure 19 heights change parallel to the wall.

Figure 19. Scatter diagram of the cloud-to-cloud distances (entire wall).

Figure 19. Scatter diagram of the cloud-to-cloud distances (entire wall).
Figure 20 shows the cloud-to-cloud distances for the four patches (A, B, C, D). As shown in
Section 2.3.3, cloud-to-cloud distances are nearest neighbor distances not the distances between the
Figure 20 shows the
corresponding cloud-to-cloud
points distances
or features, therefore the changefor thearefour
values patches
shown in
different (A,
from the B, C,method
baseline D). As
Section 2.3.3, cloud-to-cloud distances are nearest neighbor distances not the distances between
the corresponding points or features, therefore the change values are different from the baseline
method and virtual point comparison. Please note that the cloud-to-cloud distances are absolute
values, and therefore have no sign, in contrast to the change vectors obtained by the baseline method.
The maximum cloud-to-cloud distances of patch A (top-left), B (top-right), C (bottom-left), and D
Sensors 2017, 17, 26 21 of 25

Sensors 2017, 17, 26 21 of 25

and virtual point comparison. Please note that the cloud-to-cloud distances are absolute values, and
therefore have no sign, in contrast to the change vectors obtained by the baseline method. The
maximum cloud-to-cloud
(bottom-right) are about 0.076,distances
0.018, of patch
0.061 and A (top-left),
0.050 B (top-right),
m, respectively. C (bottom-left),
In patch andare
A, the distances D
(bottom-right)
consistent fromareleftabout 0.076,The
to right. 0.018, 0.061change
biggest and 0.050is onm,the
respectively.
left side ofInthe
patch
patchA,which
the distances
shows theare
consistent
uniform from
move to left to right.
the left, namelyTheinbiggest change
the opposite is on the In
Y direction. leftpatch
side B,
ofthere
the patch which
is a hole shows
without the
points
uniform
which move
is the to the
results ofleft, namely inby
an occlusion the opposite
a piece Y direction.
of wood. In patch B, there
The cloud-to-cloud is a hole
distance without
of patch B ispoints
small
which is the
compared toresults
the otherof an occlusion
patches. In by a piece
patch of wood.
C, the distancesTheon cloud-to-cloud distance
the top are greater of patch
than B bottom
on the is small
compared
which to thethe
indicates other patches.
damage that In patch C,inthe
happened distances
this on the topwhich
area is significant are greater than on
corresponds tothe
the bottom
results
which indicates
discussed the damage
in Section 3.4.2. In that happened
the middle in thisD,area
of patch theisdistances
significant
arewhich corresponds
not consistent to the results
(an obvious fault
discussed in Section 3.4.2. In the middle of patch D, the distances are not consistent
exists). However, the absolute values of the distances are relatively low which indicates that only small(an obvious fault
exists). However,
damage happenedthe absolute
in this area. values of the distances are relatively low which indicates that only
small damage happened in this area.

(a) Patch A (b) Patch B

(c) Patch C (d) Patch D

Figure 20. Scatter diagram of cloud-to-cloud distances.

Figure 20. Scatter diagram of cloud-to-cloud distances.

3.5.4. Comparative Analysis

3.5.4. Comparative Analysis
From the above analysis it follows that the three methods give similar results in terms of
From the above analysis it follows that the three methods give similar results in terms of maximum
maximum change, mean change and change tendency which directly indicates the feasibility of our
change, mean change and change tendency which directly indicates the feasibility of our proposed
proposed method. However, in our experiment, we could make use of available stable targets to
method. However, in our experiment, we could make use of available stable targets to obtain a
obtain a good registration. If, for some practical reason, a scene does not have proper locations for
good registration. If, for some practical reason, a scene does not have proper locations for targets,
targets, misalignments will affect the results of traditional methods while our proposed method will
misalignments will affect the results of traditional methods while our proposed method will not
not suffer. In addition, the results clearly demonstrate some shortcomings of traditional change
suffer. In addition, the results clearly demonstrate some shortcomings of traditional change detection
detection methods. Cloud-to-cloud distances are sensitive to the configuration of the compared point
methods. Cloud-to-cloud distances are sensitive to the configuration of the compared point clouds.
clouds. Insertion or removal of small objects, like the cables visible in Figure 12, will locally result in
Insertion or removal of small objects, like the cables visible in Figure 12, will locally result in large cloud
large cloud to cloud distances. In addition, cloud to cloud comparison is not suitable for revealing
Sensors 2017, 17, 26 22 of 25

to cloud distances. In addition, cloud to cloud comparison is not suitable for revealing motion parallel
to for example a wall. Comparing virtual points extracted from registered point clouds overcomes
these issues, as virtual point comparison is fully 3D, while virtual point extraction implicitly ignores
unwanted objects like the discussed cables. Still, comparing coordinates from different point clouds
requires a reliable registration, which is not always achievable in practice. It should be noted though
that also virtual point extraction introduces one or even two additional steps, the preferably automatic
extraction of reconstructable points, and the identification of corresponding reconstructable points
between epochs. The possibility to do so will strongly depend on the scene considered, but many
scenes do have candidate virtual points in practice. It should be noted that the methods illustrated in
this paper are largely automated. Only selecting the corresponding patches required manual work,
but very precise selection was not required.

4. Conclusions and Recommendations

In this work we have proposed two new approaches to change detection in scenes sampled by
laser scanning. Both of these approaches use so-called virtual points. For us, virtual points are 3D
locations of features that are reconstructable from the available point cloud data. To enable comparison
of virtual points, it is required that correspondences are established between the same virtual points in
different scenes. One way of using virtual points is by directly comparing 3D coordinates between
corresponding virtual points in different point clouds after registration. Advantage of this approach
is that changes between well-established features are considered. Alternatively, baselines, i.e., lines
connecting virtual points in one scene, are extracted and corresponding baselines from different epochs
are compared to identify 3D change. This approach has the additional advantage that registration is
no longer required.
For interpretation of the results in the latter case, we defined a so-called structural coordinate
system. This coordinate system is used to decompose 3D change vectors into orthogonal components
parallel and perpendicular to the object under consideration. In our case this object is the wall
subject to change Definition and use of such coordinate system is essential to properly understand
detected movement, which in general is a combination of rotation and displacement. As the baseline
method identifies relative motion, the use of such structural coordinate system is also essential for the
identification of absolute movement, as for this an external point of view is required outside the area
that is subject to change. One could therefore argue that adding this structural coordinate system is a
kind of back-door registration. Indeed, it does add additional information on the relative alignment of
the scene in different epochs. Still, adding a structural coordinate system has less requirements and is
less labor intensive then performing a fine registration.
Both new methods were demonstrated on a masonry wall located at the TU Delft Stevin laboratory
subject to seismic testing. The wall was scanned before and also after the seismic tests which caused
obvious damage. The bricks constituting the masonry wall are clear, recognizable objects that are in
addition of a suitable size for a structural damage assessment. Therefore, an algorithm was developed
that automatically estimates the brick centre locations from the laser scan intensity data. The intensity
response allowed us to separate mortar points from brick points using k-means clustering. Next,
brick outlines and brick centre locations are estimated from projected point densities of the mortar
points. Baselines connecting brick centres and targets in different scans were established, decomposed
and analyzed.
The results clearly indicate how different parts of the damaged wall move in 3D with respect to
each other. By identifying a larger spread of differences in lengths of baselines all connected to one
part of the wall, very local damage, almost at brick level, could be identified in a next, more local
analysis. An obvious disadvantage of traditional cloud to cloud based change detection is that changes
parallel to a surface are difficult to detect. In that sense cloud to cloud based change detection is a
2.5 D method while both the virtual point method and the baseline method are truly 3D as it reports
Sensors 2017, 17, 26 23 of 25

3D change vectors as a result. Designing a suitable structural coordinate system with axes parallel and
perpendicular to the wall under study in this case helped the interpretation of the results.
The two proposed methods are general in the sense that they could be applied to any situation
where repeated scan data is available, provided that it is possible in the scene under study to identify
the corresponding features. We believe that identification of features is indeed possible in many
man-made objects, which makes also the proposed method widely applicable, e.g., [46,47]. It should
also be possible to extract features in e.g., trees, outcrops and terrain covered by boulders, so the
method should also work for geomorphological application. In addition, we have to demonstrate that
for different scenarios it is indeed possible to establish a suitable structural coordinate system, which
may be difficult in the open area with ragged topography. Still, to make the method powerful would
require to have a supporting method that automatically extracts and identify suitable features in the
given application, like in this paper the brick centres extraction algorithm. What should also be studied
in future are (i) methods to estimate the quality of the feature points, which could be; (ii) propagated
to the quality of the extracted baselined, to be finally used to (iii), test whether a change in baseline is
significant [19].

Acknowledgments: We would like to thank the China Scholarship Council (CSC) for enabling Yueqian Shen to
study for one-year at the Department of Geoscience and Remote Sensing, TU Delft. We also would like to thanks
ir. H.R. Schipper from TU Delft for organizing the seismic experiment.
Author Contributions: All authors have contributed to the manuscript. The main work was carried out by
Yueqian Shen and supervised by Roderik Lindenbergh. Yueqian Shen wrote the paper with contributions from the
other authors, all of whom proofread the paper. Jinhu Wang provided helpful feedback on the experiments design.
Conflicts of Interest: The authors declare no conflicts of interests.

References
1. Vosselman, G.; Maas, H.-G. Airborne and Terrestrial Laser Scanning; Whittles Publishing: Dunbeath, UK, 2010.
2. Lindenbergh, R.; Uchanski, L.; Bucksch, A.; Van Gosliga, R. Structural monitoring of tunnels using terrestrial
laser scanning. Rep. Geod. 2009, 2, 231–238.
3. Cai, H.; Rasdorf, W. Modeling Road Centerlines and Predicting Lengths in 3-D Using LIDAR Point Cloud
and Planimetric Road Centerline Data. Comput.-Aided Civ. Infrastruct. Eng. 2008, 23, 157–173. [CrossRef]
4. Wang, K.C.P.; Hou, Z.; Gong, W. Automated Road Sign Inventory System Based on Stereo Vision and
Tracking. Comput.-Aided Civ. Infrastruct. Eng. 2010, 25, 468–477. [CrossRef]
5. Pesci, A.; Teza, G.; Bonali, E. Terrestrial Laser Scanner Resolution: Numerical Simulations and Experiments
on Spatial Sampling Optimization. Remote Sens. 2011, 3, 167–184. [CrossRef]
6. Zhang, C.; Elaksher, A. An Unmanned Aerial Vehicle-Based Imaging System for 3D Measurement of
Unpaved Road Surface Distresses. Comput.-Aided Civ. Infrastruct. Eng. 2012, 27, 118–129. [CrossRef]
7. Gil, E.; Llorens, J.; Llop, J.; Fabregas, X.; Gallart, M. Use of a terrestrial LIDAR sensor for drift detection in
vineyard spraying. Sensors 2013, 13, 516–534. [CrossRef] [PubMed]
8. Pirotti, F. State of the Art of Ground and Aerial Laser Scanning Technologies for High-Resolution Topography
of the Earth Surface. Eur. J. Remote Sens. 2013, 46, 66–78. [CrossRef]
9. Capra, A.; Bertacchini, E.; Castagnetti, C.; Rivola, R.; Dubbini, M. Recent approaches in geodesy and
geomatics for structures monitoring. Rendiconti Lincei 2015, 26, 53–61. [CrossRef]
10. Park, H.S.; Lee, H.M.; Adeli, H. A New Approach for Health Monitoring of Structures Terrestrial Laser.
Comput.-Aided Civ. Infrastruct. Eng. 2007, 22, 19–30. [CrossRef]
11. Vezocnik, R.; Ambrozic, T.; Sterle, O.; Bilban, G.; Pfeifer, N.; Stopar, B. Use of terrestrial laser scanning
technology for long term high precision deformation monitoring. Sensors 2009, 9, 9873–9895. [CrossRef]
[PubMed]
12. Zhang, F.; Knoll, A. Vehicle Detection Based on Probability Hypothesis Density Filter. Sensors 2016, 16, 510.
[CrossRef] [PubMed]
13. Truong-Hong, L.; Laefer, D.F.; Hinks, T.; Carr, H. Combining an Angle Criterion with Voxelization and
the Flying Voxel Method in Reconstructing Building Models from LiDAR Data. Comput. Aided Civ.
Infrastruct. Eng. 2013, 28, 112–129. [CrossRef]
Sensors 2017, 17, 26 24 of 25

14. Sithole, G. Detection of bricks in a masonry wall. ISPRS Arch. 2008, XXXVII, 567–572.
15. Lourenço, P.B. Computations on historic masonry structures. Progress Struct. Eng. Mater. 2002, 4, 301–319.
[CrossRef]
16. Riveiro, B.; Lourenço, P.B.; Oliveira, D.V.; González-Jorge, H.; Arias, P. Automatic Morphologic Analysis of
Quasi-Periodic Masonry Walls from LiDAR. Comput.-Aided Civ. Infrastruct. Eng. 2016, 31, 305–319. [CrossRef]
17. Lourenço, P.B.; Mendes, N.; Ramos, L.F.; Oliveira, D.V. Analysis of masonry structures without box behavior.
Int. J. Archit. Herit. 2011, 5, 369–382. [CrossRef]
18. Lindenbergh, R.; Pietrzyk, P. Change detection and deformation analysis using static and mobile laser
scanning. Appl. Geomat. 2015, 7, 65–74. [CrossRef]
19. Van Gosliga, R.; Lindenbergh, R.; Pfeifer, N. Deformation analysis of a bored tunnel by means of terrestrial
laser scanning. ISPRS Arch. 2006, XXXVI, 167–172.
20. Lindenbergh, R.; Pfeifer, N. A statistical deformation analysis of two epochs of terrestrial laser data of a lock.
In Proceedings of the 7th Conference on Optical 3D Measurement, Vienna, Austria, 29–30 October 2005;
pp. 61–70.
21. Hofton, M.; Blair, J.B. Laser pulse correlation: A method for detecting subtle topographic change using
LiDAR return waveforms. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2001, 34, 181–184.
22. Schäfer, T.; Weber, T.; Kyrinovič, P.; Zámečniková, M. Deformation Measurement Using Terrestrial
Laser Scanning at the Hydropower Station of Gabčíkovo. In Proceedings of the INGEO 2004 and FIG
Regional Central and Eastern European Conference on Engineering Surveying, Bratislava, Slovakia,
10–13 November 2004.
23. Men, H.; Gebre, B.; Pochiraju, K. Color point cloud registration with 4D ICP algorithm. In Proceedings of the
2011 IEEE International Conference on Robotics and Automation (ICRA), Shanghai, China, 9–13 May 2011;
pp. 1511–1516.
24. Jost, T.; Hügli, H. Fast ICP Algorithms for Shape Registration, Joint Pattern Recognition Symposium; Springer:
Berlin, Heidelberg, 2002; pp. 91–99.
25. Rusu, R.B.; Blodow, N.; Marton, Z.C.; Beetz, M. Aligning point cloud views using persistent feature
histograms. In Proceedings of the 2008 IEEE/RSJ International Conference on Intelligent Robots and
Systems, Nice, France, 22–26 September 2008; pp. 3384–3391.
26. Jiang, J.; Cheng, J.; Chen, X. Registration for 3-D point cloud using angular-invariant feature. Neurocomputing
2009, 72, 3839–3844. [CrossRef]
27. Rabbani, T.; van den Heuvel, F. Automatic point cloud registration using constrained search for
corresponding objects. In Proceedings of the 7th Conference on Optical 3D Measurement Technique,
Vienna, Austria, 29–30 October 2005; pp. 177–186.
28. Gressin, A.; Mallet, C.; Demantké, J.; David, N. Towards 3D lidar point cloud registration improvement
using optimal neighborhood knowledge. ISPRS J. Photogramm. Remote Sens. 2013, 79, 240–251. [CrossRef]
29. Cignoni, P.; Rocchini, C.; Scopigno, R. Metro: Measuring Error on Simplified Surfaces. Comput. Graph. Forum
1998, 17, 167–174. [CrossRef]
30. Aspert, N.; Santa Cruz, D.; Ebrahimi, T. MESH: Measuring Errors between Surfaces Using the Hausdorff
Distance. In Proceedings of the IEEE International Conference in Multimedia and Expo, Lausanne,
Switzerland, 26–29 August 2002; pp. 705–708.
31. Girardeau-Montaut, D.; Roux, M.; Marc, R.; Thibault, G. Change detection on points cloud data acquired
with a ground laser scanner. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2005, 36, W19.
32. Rusu, R.B.; Meeussen, W.; Chitta, S.; Beetz, M. Laser-based perception for door and handle identification.
In Proceedings of the 2009 International Conference on Advanced Robotics, Munich, Germany,
22–26 June 2009; pp. 1–8.
33. Trevor, A.J.; Gedikli, S.; Rusu, R.B.; Christensen, H.I. Efficient organized point cloud segmentation with
connected components. In Proceedings of the 3rd Workshop on Semantic Perception Mapping and
Exploration (SPME), Karlsruhe, Germany, 22–29 March 2013.
34. Rusu, R.B.; Holzbach, A.; Blodow, N.; Beetz, M. Fast geometric point labeling using conditional random
fields. In Proceedings of the 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems,
St. Louis, MO, USA, 10–15 October 2009; pp. 7–12.
35. Geosystems, L. Leica ScanStation C10; Product Specifications: Heerbrugg, Switzerland, 2012.
Sensors 2017, 17, 26 25 of 25

36. Soudarissanane, S.; Lindenbergh, R. Optimizing terrestrial laser scanning measurement set-up.
In Proceedings of the 2011 International Society for Photogrammetry and Remote Sensing (ISPRS) Workshop
Laser Scanning, Calgary, AB, Canada, 29–31 August 2011.
37. Abbas, M.A.; Setan, H.; Majid, Z.; Chong, A.K.; Idris, K.M.; Aspuri, A. Calibration and accuracy assessment
of leica scanstation c10 terrestrial laser scanner. In Developments in Multidimensional Spatial Data Models;
Springer: Berlin, Germany, 2013; pp. 33–47.
38. Sanchez, V.; Zakhor, A. Planar 3D modeling of building interiors from point cloud data. In Proceedings of
the 2012 19th IEEE International Conference on Image, Orlando, FL, USA, 30 September–3 October 2012;
pp. 1777–1780.
39. Jolliffe, I. Principal Component Analysis; John Wiley & Sons: Hoboken, NJ, USA, 2014.
40. Van Goor, B.; Lindenbergh, R.; Soudarissanane, S. Identifying corresponding segments from repeated scan
data. In Proceedings of the 2011 International Society for Photogrammetry and Remote Sensing (ISPRS)
Workshop Laser Scanning, Calgary, AB, Canada, 29–31 August 2011.
41. Höfle, B.; Pfeifer, N. Correction of laser scanning intensity data: Data and model-driven approaches. ISPRS J.
Photogramm. Remote Sens. 2007, 62, 415–433. [CrossRef]
42. Forgey, E. Cluster analysis of multivariate data: Efficiency vs. interpretability of classification. Biometrics
1965, 21, 768–769.
43. MacQueen, J. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability: Statistics;
University of California Press: Berkeley, CA, USA, 1967.
44. Mémoli, F.; Sapiro, G. Comparing point clouds. In Proceedings of the 2004 Eurographics/ACM SIGGRAPH
Symposium on Geometry Processing, Nice, France, 8–10 July 2004; pp. 32–40.
45. Girardeau-Montaut, D. CloudCompare. 2014. Available online: http://www.danielgm.net/cc (accessed on
19 December 2016).
46. Steder, B.; Grisetti, G.; Burgard, W. Robust place recognition for 3D range data based on point features.
In Proceedings of the 2010 IEEE International Conference on Robotics and Automation (ICRA), Raleigh, NC,
USA, 3–8 May 2010; pp. 1400–1405.
47. Hänsch, R.; Weber, T.; Hellwich, O. Comparison of 3D interest point detectors and descriptors for point
cloud fusion. ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci. 2014, 2, 57. [CrossRef]

© 2016 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/).