A study on the detection of bulging disease in ancient city walls based on fitted initial outer planes from 3D point cloud data

Bulging is a common disease suffered by ancient city walls, especially clay brick walls. Bulging affects or even threatens the stability of outer structures of ancient city walls. Therefore, accurate identification of bulges is one of the urgent issues to be addressed in ancient city wall conservation. In this study, a novel method is proposed for detecting and identifying bulge information using point cloud data collected from the outer surface of an ancient city wall. The non-bulging points are identified in the sliced point cloud of a city wall, and then the initial plane of the city wall is fitted more accurately with the aid of spatial constraints. Finally, the bulging damage is identified and extracted by comparing the initial plane with the original point cloud. This method was applied to the project of ancient city wall conservation and was compared with other existing methods in effectiveness and accuracy. The results show that the proposed method can detect efficiently budging disease in ancient city walls for architectural heritage conservation.


Introduction
Ancient city walls are architectural heritage sites that possess great historical value and reflect early urban planning and architectural techniques. One of the most well-preserved types of walls is the outer clad brickwork, which comprises rammed earth inside and bricks outside [1]. Bulging refers to a type of external deformation damage of ancient city walls in which protruding bricks on the surface of the city wall form a cavity [2], and the main manifestation of this damage is the outward displacement of the bricks on the outer surface and damage to the brick layer structure. In severe cases, the collapse of the outer clad brickwork can occur, resulting in irreversible wall damage. It is important to obtain the extent and degree of bulging in a timely manner when considering preventive conservation and repair of ancient city walls.
Currently, existing studies have focused on the identification of bulging diseases and the monitoring of changes. According to the difference in the dielectric constant between the bulging cavity and the wall, ground penetrating radar can detect the existence of bulging in city walls from the reflected wave image [3]. Changes in the bulging can be obtained from the spatial coordinates of the observation target on the outer surface of the wall as monitored by a total station [4,5], and they can also be obtained by sensors embedded in the wall, such as optical fibres [6]. However, the installation of equipment inside the wall is required in these existing methods, e.g., [6], which could cause unnecessary damage to the original city wall. or the coordinate data cannot cover the entire outer surface of city wall with sufficient density in some methods, e.g., [4,5],therefore, it is difficult to simultaneously obtain the accurate degree and range of bulging.
Terrestrial laser scanning (TLS) is a measurement technique that actively acquires high-precision threedimensional (3D) geometric information on the surface of an object [7,8]. Over the last few years, point cloud data obtained from TLS have been widely used in damage detection and deformation analysis of various types of buildings [9][10][11]. Processing point cloud data by slicing the point cloud to obtain a profile point cloud and then performing feature extraction analysis has become a common approach. For example, in tunnel boring projects, after obtaining the tunnel profile (tunnel contour lines) based on the sliced point cloud data, the tunnel deformation can be analyzed by counting the deviations in the measurement points [12][13][14]. To fit sliced point cloud data, Elberink used Markov chain Monte Carlo sampling and the Fourier interpolation method to obtain continuous smooth railway track curves [15], and Xu used the B-spline fitting method to extract the outlines of arched components [16].
In this paper, 3D TLS is adopted for the detection of bulging disease in ancient city walls. The main strategy is to acquire point cloud data through scanning the outer surface of ancient city walls and then compare them with the initial plane of ancient city walls to determine the bulging range and degree. Theoretically, the initial plane of an ancient city wall can be calculated in a very straight way if it is perpendicular to the ground surface. However, most surfaces of ancient city walls are not perpendicular to the ground surface and form a convergence angle with the plumb plane. Therefore, it is the key issue how to determine the initial plane of the city wall surface before the wall became bulging.
Currently, there are two methods that are employed to solve this issue. One method is to determine the inclination angle of the initial plane according to historical records. For example, Cao determined the height and the length of a wall according to a description of the Forbidden City wall in the Code of Great Ming Dynasty. This historic text states that "the height of the wall is three zhangs"("zhang" is a length measure unit in ancient China, 1 zhang = 3.33 m) and "the top shrinks inward by one chi, two cuns and five fens"(1 chi = 10 cuns = 100 fens = 0.33 m) [17]. These values were then used to calculate the inclination angle. Because the units used in historical records are not sufficiently accurate, the accuracy in the first measurement method is limited. Particularly, for the city wall with long history, the overall settlement, inclination, and other deformations make it difficult in determining the initial position of wall surface from historical records.
The other method is to measure and fit the spatial point coordinate data of the wall surface to obtain the initial plane. For example, 3D point cloud data are acquired by scanning the outer surface of an ancient city wall, and the feature plane method [18] or the random sampling consensus (RANSAC) algorithm [19] is then directly applied to fit the point cloud to obtain the initial plane. In this method, point cloud data are directly used to determine the initial plane without considering the different features of the point clouds in both the bulging area and the non-bulging, as a result, the initial plane calculated for the wall would be in poor agreement with the true initial plane. In addition, RANSAC method generally has an unsatisfactory fitting effect when the noise point exceeds 60% [20]. Therefore, the accuracy of the existing methods for extracting bulging disease information of city walls cannot meet the requirements of extracting bulging disease information for city wall conservation.
In this study, a new method is proposed to detect bulging disease in ancient city walls based on point clouds obtained from 3D TLS. In the new method, a 3D comparison between the fitting plane and the point cloud of the city wall is used to detect bulging areas. The key point is to determine the initial fitting plane of the city wall as a benchmark. An ideal solution is to use points in flat areas (non-bulging areas) only to fit the initial plane in order to reduce the effect from points in non-flat areas (bulging areas) such as bulging areas. This impacts the performance but also causes errors due to the noise points in non-flat areas.
In order to overcome the difficulties in the existing methods, the proposed method adopts a two-levels strategy to determine the plane from the huge size of 3D point clouds. At the first level, we cut the point cloud slices and resample the points to achieve rough noise elimination. At the next level, spatial constraints like the missing, spalling, and bulging points are introduced to further remove the noise and refine the final plane. The result of the initial plane will be derived only from point cloud data in non-bulging areas and avoid the effect of the noise data in possible building areas. Therefore, the fitted initial plane and extracted bulging areas will be more accurate.

Methodology
As described above, the existing second method does not consider the features of different area and utilizes all point cloud data in the process of determining the initial plane. It causes poor accuracy of initial plane, but also needs extra workload to meet the requirements of preprocessing including the removal of redundant point clouds, such as parapet walls, drainage outlets and ground.
To reduce or eliminate the effect of building areas (noise data) on the initial plane and improve the performance, one strategy of point cloud slicing is firstly introduced in this study. Slicing the point cloud of the city wall can sample the points at the specified height of the city wall, and fitting these points can obtain ordered discrete points representing the bulging change of the city wall at the corresponding height. Then, the discrete points in the non-bulging region are identified based on the curved spatial geometric features in the sliced data. Thirdly, plane fitting was performed on the points of the non-bulging parts to obtain the initial plane of the city wall surface, the wall bulging disease information is extracted by comparing the initial plane with the original point cloud of the wall, and an accuracy evaluation is performed. Finally, the disease information is visualized. The main workflow is shown in Fig. 1.
In general, the method includes four main parts: point cloud slicing and fitting, feature points identification, fitting of the city wall initial plane and bulging disease information extraction, and accuracy evaluation.

Point cloud data slicing
Usually, the coordinate system of the 3D point cloud data collected in the field is the instrument coordinate system, which cannot be directly input into the algorithm. For the convenience, an independent coordinate system is necessary to be established according to the structural characteristics of city walls, illustrated in Fig. 2, where the X-axis is the intersection line between the base of the wall and the ground as determined through both endpoints of the wall using the right direction (facing the wall) as the positive direction. the Z-axis trends upwards vertically, the Y-axis is perpendicular to the XOZ plane and points towards the wall. The origin of the coordinates is selected at any point on the X-axis. The 3D point cloud data are converted into an independent coordinate system to facilitate point cloud data processing and algorithm description.
When a plane passes through the point cloud data, all points in a certain neighbourhood of the plane location are extracted and projected to the plane, these projection points are called the slice of point cloud data, and the width of neighbourhood is called the slice width. The slice of the point cloud data can be viewed as the point set located on a certain plane in the point cloud data. Assuming that several planes with equal intervals that are parallel to the horizontal plane pass through the wall point cloud data, several corresponding slices are obtained (Fig. 2). Let the slice width be w, the slice interval be d 1 , and the upper and lower limits of the slice height be h max and h min , respectively. Then, a total of L sections can be obtained. The Z-coordinates range of the J-th section containing discrete points is expressed as follows: The 2D discrete points of the J-th section are resampled at an interval of d 2 to obtain N J discrete points, and the discrete points are sorted in descending order according to the X-axis values. Each section is resampled and sorted, and a point set Q is formed as follows: where K represents the K-th discrete point in the J-th section. The schematic diagram of point cloud slice is shown in Fig. 2.

Piecewise fitting of slice contours
If there is no bulging on the outer surface of the city wall, then the discrete points in the successively connected (1) (2) Q = q JK q JK = x JK , y JK ,   slices can be approximated as a straight line. otherwise, they can be represented as a curve, and the curved section can be considered a representation of the bulging area on the slice [20]. Therefore, the key to bulge extraction is to identify the curved part of the line that connects the discrete points. Due to the presence of weathering or damage to the exterior clad bricks of ancient city walls (Fig. 3), the discrete points in the connected slices can also have many small local depressions. Therefore, to obtain the part of the curve that represents the bulge in the sliced data, it is necessary to fit the discrete points in the slice to a continuous smooth curve and, at the same time, remove the small local depressions as noise points. The actual masonry spalling and missing of the wall is shown in Fig. 3. For curve fitting of the slice discrete points, only an appropriate algorithm can fit the curve to reflect the curved characteristics. otherwise, too much detail could be retained or necessary detail could be removed [21]. Here, a piecewise fitting strategy is adopted. The discrete points are first segmented, and then a polynomial is combined with the cubic Hermite interpolation function for fitting. At the same time, the improved K-means clustering algorithm is used to eliminate the noise points to ensure that the curved characteristics of the curve can be expressed, and reliable fitting accuracy can be obtained.
The discrete points of the section J in the point set Q form a point set Q 1 , which is expressed as follows: The discrete curvature of each point in the point set is calculated, and K-means clustering is performed on the discrete points according to the curvature value. According to the clustering results, the noise points are filtered through a neighbourhood search. The specific method is as follows.

Calculation of the discrete curvature of a point set
The three consecutive points in the point set Q 1 are q JK-1 , q JK , q JK+1 , and t 1 , t 2 represent the lengths of the vectors between q JK and two nearby points, as shown in Fig. 4.
The curve parameter equation composed of three points is obtained as follows: The length of two segments is: If the following conditions are met, then the expression becomes: Then a 1 , a 2 , a 3 , b 1 , b 2 , b 3 can be solved according to the following equation: Discrete curvature κ at point q JK is expressed [22] as follows:

Noise filtering
Noise filtering mainly removes the small local depressions in the slice caused by damage to the city wall. Considering that these small depressions have a small range and large curvature value at the bottommost point of the depression, this area can be searched by setting the curvature threshold of the discrete points, and at the same time, the width and depth range of the small depressions can be limited to find other points inside the depressions, thus facilitating their identification. Because the curvature of the discrete points at the bottom of these small depressions is similar, it can be obtained by the K-means clustering method. K-means clustering is a classic unsupervised clustering algorithm, and its main idea is to group similar samples together to form clusters according to the distance between samples [23]. This approach is performed on the point set according to the discrete curvature, and the cluster with the largest discrete curvature is obtained, that is, the set of the bottom points of all small depressions. Then, the two thresholds, the depth and width, of the small depression are set and defined by the difference S of the Y-axis between adjacent points and the range coefficient R ′ to determine the adjacent area. The specific noise filtering algorithm is obtained as follows: (1) The discrete points contained in the cluster with the largest discrete curvature are extracted as seed points to form a seed point set O. (2) Centerd at the seed point, within the range R ′ , a search is performed on both sides of the X-axis.
If the absolute value of the difference between the Y-coordinate of the seed point and the neighbouring point on one side is less than S, the neighbouring point is considered a noise point, and the search continues until the edge point of one side (at distance R ′ ) is determined. If the absolute value of the difference between the Y-coordinate of the seed point and the neighbouring point on one side is greater than S, the iteration is stopped on that side, as shown in Fig. 5. The search on both sides is complete. (3) The seed point set O and obtained noise points are removed, and the remaining points forms the point set Q 2 , which is expressed as follows

Piecewise curve fitting of point sets
To accurately represent the curved characteristics of the slice contours, it is necessary to adopt a reasonable (9) fitting strategy and function model [24]. Piecewise fitting is a method that divides a segment of complex discrete points into several segments of simple discrete points for curve fitting. this strategy has high fitting accuracy but maintains the curved characteristics of the curve. In this study, a cubic polynomial function is used to achieve piecewise fitting of the sliced data after the noise is removed, but the fitting curve of each segment is not continuous at the breakpoints. For a segment of discrete points, the endpoint condition of the cubic Hermite interpolation can satisfy the first-order continuity of the interpolation function at the endpoints [25,26]. In this study, cubic Hermite interpolation is used to achieve connection between the polynomial fitting curves.
Let the function expression of the cubic polynomial be expressed as follows: In Eq. (10), b is the number of curve fittings, and a 0 , a 1 , a 2 , … a b are the undetermined regression coefficients of polynomial fitting.
Let the function expression of cubic Hermite interpolation be expressed as follows: In Eq. (11),α 1 , α 2 , α 3 , α 4 are the undetermined regression coefficients of the interpolation function. If the coordinates of two endpoints are (x L ,y L ), (x R ,y R ) in the interpolation function, then the expression of constraint (12) exists and is unique [27]:

Fig. 5 Schematic diagram of seed point extension
Piecewise fitting is performed on the discrete points in point set Q 2 , and the fitting method is as follows: (1) The number of segments u for piecewise fitting of point set Q 2 is determined, the points are evenly distributed to each segment, and the number of discrete points in each segment is recorded as n 1 ,n 2 ,n 3 ,…,n u . (2) Cubic polynomial fitting is performed on the first segment of discrete points (x i ,y i ), i = 1,2,3,…n 1 , and the obtained polynomial equation is ing on the point ( x n 1 -4 , y n 1 -4 ) is performed, and the number of times the point is weighted is generally 0.2n 1 . Both ( x l 1 , y l 1 ), l 1 = n 1 − 4, n 1 − 3, n 1 − 2,…n 1 and (x i , y i ), i = n 1 , n 1 + 1, n 1 + 2,…2n are fitted with cubic polynomials, and the obtained polynomial equation is denoted as f 2 (x). Discretization f 2 (x) at x = x i , x = n 1 − 4, n 1 − 3, n 1 − 2,…2n 1 obtains discrete points ( x l 2 , y l 2 ), l 2 = 1, 2, 3,…n 2 . (4) For ( x l 1 , y l 1 ), l 1 = 1, 2, 3,…n 1 and ( x l 2 , y l 2 ), l 2 = 1, 2, 3,…n 2 , ( x l 1 , y l 1 ), l 1 = n 1 − 4, n 1 − 3, n 1 − 2,…n 1 and ( x l 2 , y l 2 ), l 2 = 1, 2, 3,…5 represent overlapping sections of two discrete data segments in the X-direction. This overlapping part is represented by an interpolation curve that employs ( x n 1 -4 , y n 1 -4 ) and ( x 5 , y 5 ) as the two endpoints of the interpolation curve, and the equation of the interpolation curve is obtained through the constraint in Eq. (12). (5) Steps (3) to (4) are repeated to obtain a complete, continuous and derivable curve for the discrete points of one section.
At the X-coordinate of each discrete point in point set Q 2 , the fitted curve is discretized, and the corrected coordinates obtained from the discretization form the overall point set P, which is expressed as follows:

Extraction of non-bulging points
To extract the non-bulging points more accurately from the sliced data, it is necessary to consider the spatial geometric characteristics of either the bulging area or the non-bulging area on the slice contour line. If there is no bulging of the city wall, the slice contour line may be a straight line, and the points in the sliced data can be directly extracted as non-bulging points. When there is only one bulging point, there may be a local bending and three inflection points (the two ends and the apex of the curve) on the slice contour line. If there are multiple bulging points in the section, there may be multiple bendings in the slice contour line, and the corresponding inflection points may also increase. For simplicity, the discrete points at the inflection points (both ends of the bending) and in a certain neighbourhood range can be treated as non-bulging points when the slope of the line formed between adjacent discrete points is small, as shown in Fig. 6. Therefore, the identification of non-bulging points can be achieved by identifying the inflection points and its type in the sliced data and determining the bending of the contour lines at the inflection points within a certain range on both sides. The identification of the type of inflection point is mainly based on the direction of the Y-component of the tangent vectors of the discrete points on the left and right sides of the inflection point. For an inflection point at one endpoint of the bending, the direction of the Y-component of the tangent vector of the discrete point on the left side is towards the wall, and the opposite direction corresponds to the discrete point on the right side. After the inflection points at the two ends of the bending are identified, the non-bulging points is determined by the slope of the connecting line between the two adjacent points. This study sets the threshold of the line slope between two adjacent points as δ and the specific algorithm for identifying non-bulging points is obtained as follows, the form of identifying non-bulging areas is shown in Fig. 6.

Inflection point identification
For the point p JF (x JF ,y JF ) in point set P, the slope of tangent at x = x JF is ε JF . then, the Y-component of the unit tangent vector at this point is: If point p JF satisfies condition (15), then point p JF can be considered the inflection point at the endpoint of a bending. This relationship is expressed as follows: (14) T x JF = sin arctan ε JF . Although some inflection points after the above treatment are closer to the initial plane than neighborhood points, but they are at the intersection of two bulging areas and still on the bulging area. In order to reduce the influence of these points on the fitting result, when the absolute value of the slope of an inflection point and other inflection points is greater than threshold δ , this point is no longer regarded as an inflection point.

Extraction of non-bulging points
Using the retained inflection points as the center, calculates the slope of the retained inflection points and discrete points on both sides. If the slope is less than δ , the discrete point is retained, and the search continues for the next discrete point on the same side. If the slope is greater than δ , the calculation on this side is complete. If two points p JF 1 , p JF 2 (x JF 1 <x JF 2 ) on both sides of the inflection point p JF satisfy Eq. (17) as follows: Then the discrete points in the interval [ x JF 1 , x JF 2 ] are denoted as non-bulging points. Each inflection points and the discrete points in its neighbourhood are identified to complete the identification of the non-bulging points in the sliced data.
Steps 1-2 are repeated to extract the non-bulging points in each section. All non-bulging points form a 3D point set P*(expression (18)), and the Z-coordinates correspond to the Z-coordinates of the original point cloud associated with each point.

Plane fitting and bulging extraction
The robust eigenvalue method is used to fit the plane comprising the non-bulging points and is based on the eigenvalue method. The calculated unit normal vector of the fitting plane and the distance from the origin to the plane are used as parameters, all the points are divided into inner points and outliers, and the final fitted spatial plane equation is obtained by iteration [28]. This method is widely used to eliminate noise points during the fitting process [29][30][31], can reduce the extraction error of nonbulging points, and can remove the few discrete points that lie on the bulging area but are arranged in parallel with the X-axis.
The specific calculation method for fitting initial planes is described as follows: (1) Plane fitting is performed for all points in point set P* through the eigenvalue method, and the obtained equation that approximates the city wall initial plane is Ax + By + Cz + D = 0. The orthogonal distance from each point along the normal line to the plane is expressed as follows: (2) The calculation of the standard deviation (SD) of the orthogonal distance H j is obtained as follows: where: (3) When the H j of a certain point is less than 2σ , the point is retained. otherwise, the point is considered an abnormal point and eliminated. (4) All of the retained points are fitted using the eigenvalue method, the orthogonal distance from each point to the fitting plane is calculated, and Steps 1 to 2 are repeated. If the orthogonal distances from all retained points to the plane are within the threshold, iteration ends, and the final fitting plane equation is obtained.
After the initial plane is obtained, the proceeded operation is to detect possible bulging areas. As described in "Piecewise fitting of slice contours", contour lines can be obtained within the slices at different heights of the city (19)  wall. However, the contour lines can only represent the horizontal bulge variation at the corresponding height of the city wall. Moreover, there are slice intervals between the contour lines, the "bulges" obtained only by the contour lines is discrete and cannot accurately reflect the degree and extent of the bulge in the vertical direction. Moreover, unlike initial plane, bulging surfaces vary from one to another without predefined well-known mathematical forms to express them; therefore, it is not able or difficult to use the same fitting method to determine all bulging areas only from the sampled discrete bulging points. Therefore, in order to obtain accurately a complete bulge situation, the initial plane of the wall is taken as the reference and then the comparative analysis is performed between the initial plane and the whole point cloud of the wall. During the process of comparison, those points which are above the initial plane will be marked as bulging. Simultaneously, the offset distance between each original wall point and the initial plane is calculated. The offset distance is defined as the orthogonal distance from each point of the wall to the initial plane along the normal line. The offset distance is also recorded as the relative height of each wall point to the initial plane.

Error evaluation
The accuracy of the bulging extraction method was evaluated from two aspects: the accuracy of curve fitting and the angle and distance deviation between fitting planes. The accuracy of curve fitting is usually evaluated using posterior variance, small error probability, and the mean absolute percentage error (MAPE) [32,33]. The accuracy analysis method for the curve fitting used in this paper is as follows: Point set P is obtained by curve fitting and discretization as follows: Point set‾P is obtained by correlating each point in set P to the point before fitting as follows: According to the residual error e JF = y JF − y JF at any point in point set P, the SD of the fitted original sequence is as follows: The SD of the absolute error is expressed as follows:  where Then, the posterior error ratio is as follows: The small error probability is expressed by: The MAPE is calculated as follows: Based on the above calculation results, the piecewise curve fitting grade and the fitting accuracy can be determined.
The algorithm robustness is verified by calculating the inclination of the fitting plane, which is obtained by processing the point clouds of each city wall segment in the same global coordinate system. The calculation of the plane inclination is shown in Fig. 7.
In the figure, A x, y, z is the unit normal vector of the fitting plane, and B(0, 0, 1) is the unit normal vector in the Z-direction. Then, the inclination angle α formed by the fitting plane and the XOZ plane is obtained as follows: (29) p e JF − e JF < 0.6745S 1 . In the same way, the inclination angle formed by the fitting plane and the YOZ plane can be obtained.

Case study and results
The Forbidden City is located in Beijing and was built in 1420. This ancient city wall type is outer clad brickwork, its height is 9.27 m, its bottom width is 8.55 m, and its top width is 6.63 m. Although the Forbidden City wall is one of the most intact ancient city walls, it has experienced bulging disease. In particular, the inner part of western section of the wall has suffered more serious bulging damage, and temporary reinforcement measures have been taken in some areas, as shown in Fig. 8.
To extract information about the bulging disease experienced by this section of the city wall, point cloud data were collected from the outer surface over a length of approximately 190 m using a 3D laser scanner. Data were acquired using a FARO Focus3D X130 scanner with a measurement rate of 976,000 points/SEC and ranging accuracy of ± 2 mm at 25 m. During the scanning process, a total of 13 stations were scanned, each station is evenly distributed in the section of the experimental city wall, about 15 m away from the city wall, and the resolution of the obtained point cloud is 5 mm. The point cloud is acquired by frontal scanning. The material of each wall brick is the same. The wall bricks are laid by mortar, but the mortar is not exposed to the outer surface of the wall, so there is no difference in reflectivity caused by the material of the wall.
We use Geomagic Studio software for point cloud registration. Considering that the city wall scaffolds and the large amount of data could lead to a decrease in the calculation speed, the city wall was divided into three segments of 29 m, 54 m and 59 m in length. The preprocessed point cloud is shown in Fig. 9.

City wall point cloud slicing
Using the third segment of the city wall as an example, the point cloud of the preprocessed city wall was sliced. The lower and upper limits of the slice height (h min and h max ) are 0.2 m and 8.6 m, respectively, the slice width w is 5 mm, the slice interval d 1 is 30 cm, and a total of 29 sections were obtained, as shown in Fig. 10a. The slices of each section were sampled with horizontal interval d 2 of 50 cm, the discrete points in each section were sorted from small to large according to the X-axis value, and the results after sampling are shown in Fig. 10b.

Fitting contour lines
Using the slice at a height of 5 m (the 17th section) as an example, the Z-coordinate change of the discrete points in the same section is within 1 mm, as shown in Fig. 11a. The discrete points contained in each section were clustered according to the discrete curvature, and the   Fig. 11b. According to the length of the city wall and the sampling interval, the slice was evenly segmented according to the number of discrete points. In this study, the slice was divided into nine segments. Among them, the first to the eigth segments each of the fitted curve, each contain 11 discrete points, and the remaining segment contained 10 discrete points. The results after piecewise fitting are shown in Fig. 12, i.e., a first-order continuous and derivable smooth curve.

Identification of non-bulging points
The fitted curve was discretized again at the X-coordinate of each discrete point to obtain a new set of points. After the slope of the tangent line and the Y-component of the tangent vector of each point in the new point set were calculated, the identification of bulging points and extraction of non-bulging points were performed. The direction of the Y-component of the tangent vector of each point is shown in Fig. 13.
The slope threshold δ was set to 0.004 and the extracted non-bulging points are shown in Fig. 14a. The identified non-bulging points in both the 7th section (at 2 m) and the 27th section (at 8 m) are shown in Fig. 14b, c.

Extraction and visualization of wall bulging points
The non-bulging points of each section were extracted as shown in Fig. 15a. All these points were fitted using the robust eigenvalue method, and the results of robust plane fitting were obtained through iterative operations and are shown in Fig. 15b.
The city wall initial point cloud was compared with the fitting plane, and the bulging points of this segment of the city wall are shown in Fig. 16c. Figure 16a, b, respectively show the bulging of first and second segments. The bulging can be clearly presented in the colour map after visualization. For the city wall studied in this paper, there are four bulges with a horizontal length greater than 10 m, the bulging height usually reaches 15-20 cm, and for the most severe area, the bulging height can reach 30 cm.

Slicing Parameter selection and its effect on results
It is a crucial step to select the appropriate parameters for the performance of the proposed method. The main parameters interval point cloud slicing, piecewise fitting and non-bulging point extraction. The requirements for determining the parameters and the correlations between parameters are discussed in this section.

Parameters of point cloud slicing
The main parameters are listed in the Table 1.
The selection of h max , h min should avoid the influence of the redundant point cloud data from non-wall surfaces such as the drainage outlet and the parapet wall. therefore, h max should be large enough to avoid the   interference of the ground point clouds, and h min should be a certain distance away from the parapet wall. The vertical interval d 1 and the horizontal sampling interval d 2 directly determine the number of discrete points after sampling the wall point cloud. If d 1 and d 2 are too small, a large amount of redundant data will degrade the performance of the algorithm. For example, when two slices are on the same layer of brick, the contours obtained by each slice are almost the same, so the extracted non-bulging points are repeated and too dense. At the same time, too large values of d 1 and d 2 will cause distortion of the bulging characteristics of the sampling results. In terms of Shannon-Nyquist sampling theorem [34,35], the sampling frequency should be greater than twice the highest frequency in the signal. Therefore, the upper limit of d 1 and d 2 depends on the size of the bulging area in the actual experimental wall segment. According to the preliminary survey, the smallest bulging area has a horizontal length of about 3 m and a vertical length of about 1 m. Figure 17 shows the two smallest bulging areas in the experimental wall section. Therefore, in this example, the value of d 1 should not be greater than 0.5 m, and the value of d 2 should not be greater than 1.5 m. Considering that the cladding bricks outside the wall are laid in the pattern of gardon wall bond and the length is about 0.45 m, the width is half of the length, and the thickness is about 0.1 m. In order to avoid multiple sampling on the same brick, the proper value range of d 1 is from 0.1 m to 0.5 m, and the value range of d 2 is from 0.225 m to 1.5 m. The actual situation of two bulging is shown in Fig. 17.
To select a proper slice width w, the density of the actual point cloud data should be considered. Simultaneously, it should also ensure that the sliced point clouds contain sufficient discrete points with even distribution. In this study, the resolution of the point cloud was 5 mm, for the reason of computing speed, w is the most appropriate point cloud resolution, so w is set to 5 mm.

Piecewise fitting parameters
The main parameters for piecewise fitting are listed in Table 2.
The number of cluster centers V could interfere to some extent with the filtering of points at the damage area of brick walls. When V is smaller, the number of discrete points in each cluster could increase, which may lead to severe local data missing. When V is larger, for the cluster with the greatest discrete curvature, the included discrete points on obvious damage area of walls might be incomplete. Besides a reasonable   setting for V, the integrity of the identified original seed points after clustering should be maintained. In this study, when V is set to 4 or 5, the clustering result is better. When V is set to 3 or 6, the result is shown in Fig. 18. The adjacent area range coefficient R ′ and the Y-coordinate difference coefficient S can be used to identify one or multiple adjacent defect points on the basis of the cluster identification seed points, and the two need to be determined synchronously according to the actual brick damage length and depth of the city wall. R ′ should ensure that a certain seed point is centered, usually between 1 and 3 times the resampling interval, so as to avoid false identification of defect points using K-means clustering when there is no obvious defect point in a slice of a layer. S is an empirical value, usually considering the depth of the brick defect, it is more appropriate to take 1 cm to 2 cm.
To determine u, the number of discrete points in each slice needs to be considered. Using the method in "Piecewise fitting of slice contours", there are at least five discrete points in each segment after segmentation. When there are too many discrete points in each segment, loworder polynomial fitting may not be able to obtain the proper bulging curve of the encasement brick. The cubic curve contains two inflection points, and in this study, the smallest bulging area has a horizontal length of about 3 m. Combined with the actual situation, the length of segments are set about 4.5 m or slightly larger, which is enough to ensure the accuracy and curve characteristics. Table 3 shows the required parameter for non-bulging point extraction.

Parameters for non-bulging point extraction
Threshold δ plays a key role in the process of identifying non-bulging points. When δ is too large, the identified non-bulging points are inaccurate, and δ is too small, it could be difficult to completely identify the non-bulging points. δ is positively correlated with the bulge height of city wall. Generally, when the maximum bulge height of   the wall is about 30 cm, δ can be set to 0.003, which can be enlarged or reduced in the same proportion with the maximum bulge height.

Precision evaluation
The accuracy of the curve fitting was verified. In total, there are 2962 discrete points used to fit each curve, and 29 fitted curves are obtained. The model accuracy evaluation criteria for curve fitting are shown in Table 4, and the classification criteria for fitting accuracy are shown in Table 5. The calculated posterior error ratio C was 0.1886, the small error probability p was 0.9956, and the MAPE was 0.6876. In summary, the piecewise curve fitting method used to fit the contour line of the city wall has high accuracy.
Bulging extraction was performed on the three segments of the city wall, and the final fitting planes were obtained in the same global coordinate system, as shown in Fig. 19.
The inclinations of the three planes with both the XOY and YOZ planes are shown in Table 6 (unit: degree).

Comparison with existing methods
The existing two methods mentioned above were compared with the proposed method in this paper. The RANSAC algorithm [36,37] was directly used for plane fitting on the point cloud data of the city wall, and the corresponding bulging disease information was extracted. In addition, the initial inclination angle of this segment of the wall was determined based on the historical recorded receding angle, and the initial plane of the wall was obtained and compared with the original point cloud to extract the bulging disease information. Taking the first segment of the city wall as an example, the visualization results of the three methods for extracting bulging diseases are shown in Fig. 20.
In order to verify the proposed method and compare the accuracy of extracted bulging areas in different methods, the point clouds of repaired wall are used as the actual city wall plane to identify the bulging areas in the city wall before repaired.
The overall bulging area and maximum bulging value of the wall obtained by using the repaired wall as actual  bulging areas and three different methods are shown in Table 7. The actual photo and point cloud of the repaired wall are shown in Fig. 21. Figure 22 illustrates the results of comparing bulging areas extracted with different methods.
Although the repaired wall plane cannot be absolutely considered as the initial plane, the repairment does not change the shape of the wall. Therefore, the repaired wall plane is highly close to the initial plane. As shown in Table 7, our method is far more accurate than RANSAC method and historical record method. Compared with RANSAC method, the accuracy in the maximum bulging value and bulging area is improved by 23% and 49% respectively; while compared with historical record method, the accuracy is improved by 32% and 68%, respectively.
At present, many studies have shown that bulge of city wall is mainly caused by rainwater seepage from the top of the city wall and freezing-thawing [2,5,19,38,39], therefore, the bulges generally develop downward from the top of the city wall, and the bulge development trend is obvious. According to the above visualization results, the bulges extracted using the two existing methods originates from the middle part of the city wall and extends around, which is different from the actual bulging pattern. However, the bulging areas obtained by the method in this paper are mainly concentrated on the top of the city wall, and all bulges spreads downwards from the top. This also indicates the new method can reflect the actual situation more accurate than the other two methods.
At the same time, a small number of points located in spalling or missing bricks show negative bulging values in the visualization. However, the visualization results of the existing methods all contains a large area of negative bulging values, and the maximum may even reach negative 15-20 cm, as shown in Fig. 20. This result means that there is a large area of wall bricks falling off, which does not conform to the actual preservation of the city wall. Therefore, the initial plane obtained by existing methods is inaccurate.
The quantitative comparative analysis above shows that the methods for directly obtains the initial wall plane can be largely affected by the discrete points of the bulging area, and the fitted initial wall plane is offset to the outside of the wall, which makes the extracted bulging areas and height inaccurate. Due to the poor measurement precision of the historical data, and the city wall deformations, such as settlement and inclination, there is a certain gap between the recorded receding angle of the city wall and the actual receding angle, so the repaired actual initial plane of the city wall has a large error. Overall, the bulging results extracted by the proposed method in this paper are close to the actual bulging pattern, and the extraction of the bulging value and area has been significantly and reliably enhanced.

Conclusions
This study provides a practical and feasible research idea for bulging extraction of regular planar walls. Because our method depends on the accuracy of the fitting plane, we design two levels strategy to eliminate the noise at the bulge of the city wall. At the first level, we cut the point  cloud slices and resample the points to filter the rough noise. At the next level, on the basis of a method of fitting the contour line of the city wall, the spatial constraints like the missing, spalling, and bulging diseases of the city wall are introduced to filter the detailed smaller noise points. Because RANSAC method only obtains the fitting plane through one-step continuous iteration on the basis of random selection of points, it is not suitable for the case that the overall bulging range of the wall is large. Therefore, this method can't be applied to severe bulge disease detection. Compared to RANSAC method, our method can handle huge size data more efficiently and accurately, breaking through the degradation and limitation of RANSAC method that the error points of plane objects generally cannot exceed 60% of the total target. Therefore, the proposed method could be generalized to determine any geometry planes from the 3D point clouds.
In this paper, an actual city wall with bulging disease was taken as an example to verify the accuracy and reliability of the proposed method. The piecewise curve fitting satisfies the high-precision fitting requirement, avoids overfitting, and accurately reflects the bulging condition. The results obtained using the proposed method in this paper were also compared with the results obtained by direct fitting of the original point cloud data. The results showed that the proposed methods can avoid the effect from noise data and obtain more accurate and reliable results of building extraction than existing methods.
The method proposed in this study can be applied to the bulging extraction of walls with other shapes, such as brick walls, rammed earth walls, and concrete walls. However, the parameters used in the method need to be adjusted according to actual conditions. Realizing the adaptive selection of parameters could further improve the operational efficiency. At the same time, bulging at the bottom of city walls is usually rare, and the points at the bottom of a city wall can be weighted according to a certain rule in the process of initial plane fitting and obtain more accurate results.
However, it should be also noted that the proposed method in this study is focused on the ancient city wall with flat brick structure extending in a straight line. In general, most of the ancient city walls in China are buildings with regular flat shapes such as Xi 'an Ancient City Wall, Pingyao Ancient City Wall, and the ancient city wall of the Forbidden City in our case analysis [40]. The bulge of irregular walls will be taken into consideration in our further research.