Visual reconstruction method of architectural space under laser point cloud big data

2.1.1. Laser point cloud big data collection in building space

In order to obtain the complete 3D laser point cloud data, the following layout principles of 3D laser scanners should be followed: First, through on-the-spot reconnaissance, the operation scope should be arranged according to the topographic features and features of the survey area, and scanning stations should be arranged, and the positions of scanning instruments should be reasonably arranged to ensure that the data of adjacent stations have overlapping parts [13]. At the same time, the number of stations should be reduced as much as possible, so as to reduce the amount of data and the splicing error before the stations. Secondly, ensure less occlusion and avoid interference from other people, so that the data contains complete building space information.

When collecting the big data of laser point cloud in building space, the building is scanned by Rigel VZ-600 3D laser scanner, the scanning speed is 127,000 pts/s, and the point spacing of the scanner is < 1 mm; The distance accuracy is ±6 mm, and the scanning field of view is 360°×100° (horizontal× vertical). Set the scanning mode of collecting laser point cloud big data in building space as long distance mode. The three-dimensional point cloud data of building space is obtained at different sites, so that the building can be completely scanned and the building target can not be missed.

2.1.2. Laser point cloud big data organization method based on point cloud pyramid

The big data of building laser power supply collected by three-dimensional laser scanner contains massive point cloud data. The point cloud pyramid method is used to organize and process the collected building laser point cloud big data, which simplifies the visual reconstruction process of building space visualization. Point cloud pyramid is a kind of point cloud data organization structure similar to image pyramid type. Its basic principle is the same as that of the image pyramid, that is, the idea of layering and blocking is used to organize point cloud data [14], but the difference is that the image data adopts the method of image resampling and two-dimensional blocking, while the point cloud pyramid uses point cloud data thinning and spatial blocking to organize point cloud data. The pyramid structure with a factor of 2 is adopted. During the construction, the laser point cloud data of the original building are thinned out according to the different thinning ratios [15]. The point cloud data size of the upper layer in the pyramid structure is 1/4 of that of the lower layer, and the data range is 4 times that of the lower layer. The construction method of point cloud pyramid is similar to that of image pyramid, which is mainly divided into two steps: data thinning and data blocking. The specific steps are as follows.

2.2.1. Feature extraction of building space lines based on improved mean-shift

The amount of three-dimensional point cloud data in building space is huge, so it takes a long time to directly use point cloud data for plane segmentation, although the method is intuitive. At the same time, when the three-dimensional laser point cloud data of building space contains curved walls, building curved walls by plane or piecewise linear fitting will cause the model to be imprecise, and surface modeling can be transformed into curve extraction by building slices, so an effective curve feature extraction method is needed. Using the improved Mean-shift point cloud line feature extraction method, the building space line feature is extracted. Firstly, the concepts of straight line and curve are defined. Line features include straight line segments and non-straight curve segments. A straight line is a series of points with the same direction. On the contrary, the direction of a curve will change from one point to another. The curve is represented by a series of uniform sampling points, while the straight line can be represented by both uniform sampling points and two points, depending on the specific situation. In order to simplify the problem of feature extraction, the original unordered point cloud is transformed into an orderly organization form of hierarchical structure [20], that is, it is processed by slicing. The basic idea of the point cloud slicing method is that the original point cloud is segmented by a set of parallel planes along a given direction and at a given interval, and the layered point cloud is projected onto the corresponding reference plane to form an outline point cloud slice. Point cloud slicing technology can segment scattered point clouds in layers, and the segmentation direction can be user-defined direction or along the $X$ axis, $Y$ axis or $Z$ axis.

Based on the assumption of Manhattan space, the reconstruction direction is generally select $Z$ axis direction when building space modeling, the number of division layers is user-defined. Suppose there are scattered point clouds $P=\left\{p_1,p_2,\cdots ,p_n\right\}$, $p_i=\left\{x,y,z\right\}\in R$, point cloud coordinate range can be recorded as $(x_{\mathrm{m}\mathrm{i}\mathrm{n}},y_{\mathrm{m}\mathrm{i}\mathrm{n}},z_{\mathrm{m}\mathrm{i}\mathrm{n}})~(x_{\mathrm{m}\mathrm{a}\mathrm{x}},y_{\mathrm{m}\mathrm{a}\mathrm{x}},z_{\mathrm{m}\mathrm{a}\mathrm{x}})$. Suppose along $Z$ axis direction for hierarchical segmentation of building space point cloud data, the number of layers is $n$, the thickness of each layer is:

where, $z_{thickness}$ is the layered thickness.

The layered point cloud is sliced by the projection plane method. For the layered point cloud, the middle position is defined as the projection reference plane of the layer, and all the data points in the layer are projected to the reference plane to form a point cloud slice.

The calculation formula of each reference plane is as follows:

where, $z_i$ represents the projection reference plane of the layer $i$, and the normal vector of the plane points in the positive direction of the $Z$ axis, $z_{\mathrm{m}\mathrm{i}\mathrm{n}}$ represents$z$value of the lowest projected reference plane.

2.2.2. Line tracing with double radius threshold

The double radius threshold line tracing method is used to extract the continuous polyline from the point cloud data of building space. Using distance threshold $d_1$ speed up the search and ensure the continuity and smoothness of the building polyline. The basic steps of line tracing algorithm with double radius threshold are as follows:

(1) Construct a $KD$ tree of the point cloud slice;

(2) Randomly select a point and search for the nearest neighbor within its $d_2$ radius, computing the $PCA$ main direction of the point cloud according to the nearest neighbor set $Q$. Search along the main direction, and then search along the reverse direction of the main direction.

(3) Perthrough the nearest neighbor within the radius $d_2$ to obtain all points whose distance to the current point is greater than $d_1$, the point with the smallest angle $\mathrm{\Delta }\theta$ was added to the curve branch as new points and marked as a node. All points that less than $d_1$ are marked as visited.

(4) Until there are no neighboring points that meet the conditions, merge two curves searched in opposite directions and return the curves.

(5) Repeat steps (2)-(4) until all slicing points are marked as visited and branching points.

2.2.3. Straight line discrimination of building space

The curve obtained by tracing consists of a series of nodes, and the distance between adjacent nodes is relatively uniform. After the online tracking step, the detected building wall lines need to be divided into straight lines and curves. It is assumed that the detected building curve can be defined by defining a pair $(p,v)$ to represent, in which $p$ is a point on a straight line, $v$ is the direction vector of a straight line. The value of $p$ takes the median of all vertex positions in the polyline. $v$ is a normalized direction vector, and its value is determined by the median of the direction vector of each component straight line segment of a polyline:

The calculation formula of the distance from each vertex of the building to the estimated straight line $d_i$ is as follows:

The angle ${\theta }_i$ between the vector from each point $p_i$ to the center point $p$ and the linear direction vector is calculated as follows:

When $\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_i\right)>\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Delta }\theta \text{'}\right)‖\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_i\right)<\mathrm{c}\mathrm{o}\mathrm{s}(\pi -\mathrm{\Delta }\theta \text{'})$ and $d_i<\mathrm{\Delta }d$ established simultaneously at the same time, $L_i$ is a straight line; Otherwise, $L_i$ is a curve.

2.2.4. Regularization of building space

Walls, doors, windows and other structures in architectural space usually present plane or straight-line geometric characteristics, mostly parallel or orthogonal. There is noise in 3D laser point cloud data, and the line feature extraction results will deviate from these rules. Therefore, it is necessary to regularize 3D laser point cloud data of buildings. The initial values of parameters of straight line segment are obtained by using the median of line point coordinates and the median of line direction. The straight line parameter is expressed as $(p_i,{\overrightarrow{v}}_i)$, record the included angle between other straight lines and reference straight lines as $\theta$, $\mathrm{\Delta }\theta$ is the angle threshold, and the expression can be obtained as follows:

If the first line refers to a straight line containing the maximum number of point clouds $L_{ref}$, then the conditional rules of the other orthogonal and parallel lines $L_i$ are defined as: When $\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)>\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Delta }\theta \right)‖\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)<\mathrm{c}\mathrm{o}\mathrm{s}(\pi -\mathrm{\Delta }\theta )\text{,}$$L_i$ is parallel; When $\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)<\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{\Delta }\theta \right)\text{,}$$L_i$is orthogonal.

In the process of point cloud data regularization, in order to obtain the optimal straight line parameters for each straight line, the constrained least square method is used to solve this problem. For each straight line, it can be expressed by a straight line equation as follows:

The straight line equation corresponding to the straight line segment projected by the reference wall line can be expressed as:

If the line is parallel to the reference line, there are:

If the line is perpendicular to the reference line, there are:

Using the above process, a linear equation group is obtained. $Ax=b$, in which $A$ is a coefficient matrix representing a system of linear equations, $x$ represents an unknown parameter vector, because $b=\varphi$, can get $x=(c_0,c_1,\cdots ,c_i,a_0,b_0{)}^T$. QR decomposition is made with coefficient matrix $A$, make the singular value SVD decomposition of the triangular array $R$ in the decomposition result to obtain the final solution of $x$. After the wall lines of buildings are regularized vertically and parallelly, the collinearity rule is used to merge the parallel straight line segments whose mutual distance is less than a given threshold, and the collinearity line segments are merged by the average method.

2.3.1. Translation matching of characteristic points of building laser point cloud big data

The initial value of translation matrix of laser point cloud big data in building space is determined, which provides the basis for feature point matching. For the selection of the initial value of the translation matrix, the centroid of two reconstructed point clouds is obtained. Let an arbitrary set of three-dimensional point clouds $P$, which should contain complete three-dimensional structural information of the target, and solve the number of point clouds in each point $P_i$ domain, determines the weight of a single point, which is inversely proportional to the number of point clouds, that is, when the neighborhood of a certain point $U(P_i,\delta )$ its number of inner point clouds is $n$, its weight is $k/n$ ($k$ is constant). When the number of point clouds in the neighborhood is 0, this point is likely to be a wrong matching point, so its weight is defined as 0, and then the weighted average coordinates of all points are calculated:

Similarly, it can be calculated $\stackrel{-}{y}$, $\stackrel{-}{z}$ and finally get the point $O_p{\left[\stackrel{-}{x},\stackrel{-}{y},\stackrel{-}{z}\right]}^T$ which is the geometric center coordinate of the point cloud, and then, the weighted average distance from all point clouds to $O_p$ is calculated as follows:

Similarly, the weighted average distance ${\stackrel{-}{s}}_q$ from point to $O_q$ within point set $Q$ can be calculated, and the initial scale factor $\lambda$ can be calculated by formula $\lambda ={\stackrel{-}{s}}_q/{\stackrel{-}{s}}_p$.

The coordinate difference between two groups of point clouds is used as the initial value of translation matrix and applied to plane registration of point cloud data:

Finally get $\lambda T$ is the initial value of registration.

2.3.2. Spatial matching of big data features of building laser point clouds

Through multi-view model registration, the visual reconstruction of architectural space is completed. The process of model registration is to match and align the point clouds obtained from different viewpoints, which is an important step in the process of obtaining the object surface contour and 3D reconstruction of the scene. Find the two closest points in two point clouds on the object, and use them as corresponding points in the process of processing, find out all such point pairs, and then calculate the rotation and translation transformation between the two point clouds according to these point pairs. The purpose of model registration is to find a rigid transformation between two local point clouds, so that one point cloud can overlap with the other point cloud best after transformation. Iterative Corresponded Point, ICP method (ICP) is used to reconstruct the architectural space visually. This method gives two three-dimensional data points set in different coordinate systems, and solves the spatial transformation of these two three-dimensional data points sets, so that they can accurately match in space.

If the first three-dimensional data point set is represented by $\left\{P_i\left|P_i\in \right.R^3,i=\mathrm{1,2},3,\cdots ,N\right\}$, the second three-dimensional data point set is represented by $\left\{Q_i\left|Q_i\in \right.R^3,i=\mathrm{1,2},3,\cdots ,N\right\}$, the objective function of two point sets alignment matching is as follows:

ICP algorithm is actually an optimal matching algorithm based on the least square method, which repeatedly executes the process of “determining the corresponding relation point set-calculating the optimal rigid body transformation” until it meets a convergence criterion representing correct matching. The purpose of ICP algorithm is to find the rotation matrix between the reference point set and the target point set $R$ and translation vector $T$ the transformation between two sets of matching data can meet the optimal matching under a certain metric criterion. Make the coordinate of the target point set $P$ is to be $\left\{P_i\left|P_i\in \right.R^3\right\}$, the coordinates of reference point set $Q$ is $\left\{Q_i\left|Q_i\in \right.R^3\right\}$, in the $k$ iteration, the point coordinate corresponding to the coordinate of the point set $P$ is $\left\{Q_i\left|Q_i\in \right.R^3\right\}$, solve the transformation matrix between $P$ and $Q$, the original transformation is updated until the average value of the distance between the data is less than the given threshold. The detailed steps are as follows:

(1) Take a subset $P_i^k\in P$ contained in the target point set $P$.

(2) Solve a reference point set $Q$ with its corresponding point $Q_i^k\in Q$, make $\left|P_i^k-Q_i^k\right|=\mathrm{m}\mathrm{i}\mathrm{n}$.

(3) Determine the rotation matrix $R_i^k$ and translation vector $T_i^k$, to satisfy the value of its smallest ${\sum }_{i=1}^N(Q_i^k-(R{P}_i^k+T){)}^2$.

(4) Calculate $P^{k+1}=\left\{P_i^{k+1}\left|P_i^{k+1}=R^k{P}_i^k+P^k,P_i^k\in P\right.\right\}$.

(5) Calculate $d^{k+1}=\frac{1}{n}{\sum }_{i=1}^n(P_i^{k+1}-Q_i^k{)}^2$.

(6) If $d^{k+1}$ is greater than or equal to a given value $\tau$, return (2) until $d^{k+1}<\tau$ or the number of iterations $k$ is more than the preset maximum number of iterations.

The ICP algorithm is calculated by the optimal analytical method $P$ and $Q$ with its corresponding points in the two-point set are matched by the transformation between them, and the registration parameters between the data sets are optimized by iteration. The singular value decomposition method is used to optimize the analysis. If the motion parameters are represented by matrix $T=(R,t)$, point set of $P$ and $Q$ are given at the same time. At any point $p_i$, after transformation, it is a point such as $p_{i+t}$. The error of that corresponds to the corresponding point $q_i$ is $(q_i-R{p}_{i+t}{)}^2$. The total error is:

In order to convert the rotation matrix $R$ and displacement vector $T$ we subtract the coordinates of all points from the coordinates of the center of mass of these points. Set $p$ is the center of mass of point set $P$, set $q$ is the center of mass of point set $Q$, then the displaced point is: ${p\text{'}}_i=p_i-p$, ${q\text{'}}_i=q_i-q$. Use a unit quaternion $q$ represent rotation matrix $R:R=R\left(q\right)$, where the unit quaternion $q=(q_0,q_x,q_y,q_z{)}^T$ meet the conditions: $q_0\ge 0$, $q_0^2+q_x^2+q_y^2+q_z^2=1$.

The cross covariance matrix of point set $P$ and point set $Q$ is:

Among them, $p$ and $q$ respectively represent the center of mass point sets $P$ and $Q$. Using anti-symmetric matrix $A=\mathrm{\Sigma }-{\mathrm{\Sigma }}^T$ and column vector $\mathrm{\Delta }=(A_{23},A_{31},A_{12}{)}^T$ construct a 4×4 symmetric matrix:

Among them, $I_3$ is a 3×3 identity matrix. The unit eigenvector corresponding to the maximum eigenvalue of matrix $Q\left({\mathrm{\Sigma }}_{pq}\right)$ is the quaternion that satisfies the rotation matrix, and the rotation matrix is calculated. After calculating the rotation matrix $R$, combined with the translation matrix $T$, each point in the 3D point cloud $A$ reconstructed by one viewpoint using $R$ and $T$, is transformed according to equation $B=AR+T$ to obtain the coordinate value of this point cloud under the coordinate system of the 3D point cloud $B$ reconstructed from another perspective, the two groups of point clouds are placed in the same coordinate system to obtain the final visual reconstruction result of architectural space visualization. The results of visual reconstruction of architectural space are visually displayed by Lumoin3D software, which is a three-dimensional visualization software, including a rich model base, and intuitively displays the results of three-dimensional visual reconstruction of architectural space.