A Method and System for Building 3D Modeling Based on Unmanned Aerial Vehicle (UAV)

Abstract

This invention relates to the field of 3D modeling technology, and discloses a method and system for 3D architectural modeling based on unmanned aerial vehicles (UAVs). The method includes acquiring a sequence of architectural images; obtaining a dense motion distribution map based on the sequence of architectural images; obtaining a dynamic region mask based on the dense motion distribution map; obtaining a clean image sequence based on the dynamic region mask; obtaining a set of 3D contour lines based on the clean image sequence; obtaining a refined set of structural points based on the set of 3D contour lines and a preset architectural geometry template; obtaining a precise 3D point cloud based on the refined set of structural points; obtaining a preliminary architectural model based on the precise 3D point cloud; and obtaining a complete architectural model based on the preliminary architectural model. This method improves the accuracy of 3D architectural reconstruction in complex environments.

Description

Technical Field

[0001] This invention relates to the field of 3D modeling technology, and in particular to a method and system for 3D architectural modeling based on unmanned aerial vehicle (UAV) assistance. Background Technology

[0002] Currently, in the fields of modern urban planning and building management, the accurate acquisition and presentation of architectural spatial information using advanced technologies is particularly crucial, as it directly relates to urban development efficiency and the safe and sustainable development of infrastructure. 3D architectural modeling requires accurate spatial data as its foundation; however, in dynamic and complex environments, data acquisition is often hampered by moving objects such as pedestrians and vehicles, making it difficult to distinguish between useful and irrelevant content during information processing, thus affecting the realism and reliability of the modeling results. Therefore, there is an urgent need to develop technologies capable of effectively identifying and retaining key spatial information while eliminating interfering factors in diverse scenarios, in order to improve the accuracy and application value of 3D modeling.

[0003] In one existing technology, a drone first orbits the target building, capturing overlapping sequences of images covering the building's surfaces using its onboard visible light camera. Then, software based on Structure of Motion (SfM) principles automatically extracts and matches a large number of two-dimensional feature points from the images, generating a three-dimensional point cloud through triangulation. In this process, the software primarily extracts and matches feature points based on pixel color, texture, and brightness variations. Next, simple filtering based on spatial distance or density is applied to the generated point cloud to remove obviously outlier noise points. Finally, the processed point cloud is used as input to generate a three-dimensional mesh model of the building using a surface reconstruction algorithm. In this process, all feature points captured and successfully matched by the camera, regardless of whether they belong to the building itself or are temporary pedestrians or vehicles, are treated equally as geometric data for calculation. However, existing technologies rely on color and texture for feature extraction, failing to distinguish between static building features and dynamic distractions. Moving objects are incorrectly matched as geometric data, resulting in a point cloud containing significant noise.

[0004] Therefore, existing technologies cannot improve the accuracy of 3D reconstruction of buildings in complex environments. Summary of the Invention

[0005] This invention provides a method and system for building 3D modeling based on unmanned aerial vehicles (UAVs) to improve the accuracy of building 3D reconstruction in complex environments.

[0006] Firstly, in order to solve the above-mentioned technical problems, the present invention provides a method for building 3D modeling based on unmanned aerial vehicle (UAV) assistance, comprising: Acquire building image sequences; Based on the building image sequence, the pixel displacement is calculated using the optical flow method, and the pixel displacement is analyzed into motion amplitude and motion direction to generate a dense motion distribution map. Based on the dense motion distribution map, candidate motion point sets are selected based on motion amplitude, and directional consistency analysis is performed on the candidate motion point sets to mark the dynamic region mask. The dynamic pixels marked by the dynamic region mask are removed from the building image sequence to obtain a static background pixel set. The missing regions in the static background pixel set are filled by bilinear interpolation to obtain a clean image sequence. Stable edge pixels with stable temporal sequence are selected from the clean image sequence to form a stable edge set. Based on the preset camera calibration parameters, the stable edge set is back-projected in three-dimensional space to construct a three-dimensional contour line set. Based on the set of three-dimensional contour lines and the preset architectural geometry template, a geometric conformity comparison is performed, and three-dimensional points that meet the preset geometric tolerance are selected to form a refined set of structural points. Incremental motion structure recovery is performed based on the refined structural point set to obtain an initial three-dimensional point cloud. Then, bundle adjustment optimization is performed on the three-dimensional coordinates in the initial three-dimensional point cloud that are greater than a preset reprojection error threshold to obtain an accurate three-dimensional point cloud. Based on the precise 3D point cloud and the preset global reference direction, surface normal vector estimation and consistency adjustment are performed to obtain an oriented normal vector field. Then, the candidate triangular facet set constructed based on the oriented normal vector field is screened for normal vector consistency to obtain a preliminary building model. Detect and repair mesh holes in the preliminary building model to obtain a closed surface mesh, and perform mesh subdivision and smoothing processing based on the closed surface mesh to obtain a complete building model.

[0007] Secondly, the present invention provides a UAV-assisted 3D building modeling system, comprising: The data acquisition module is used to acquire building image sequences; The optical flow module is used to calculate pixel displacement based on the building image sequence using the optical flow method, and to resolve the pixel displacement into motion amplitude and motion direction to generate a dense motion distribution map; The analysis module is used to filter out candidate motion point sets based on motion amplitude according to the dense motion distribution map, and to perform directional consistency analysis on the candidate motion point sets to mark dynamic region masks. The filling module is used to remove dynamic pixels marked by the dynamic region mask from the building image sequence to obtain a static background pixel set, and to fill the missing regions in the static background pixel set with bilinear interpolation to obtain a clean image sequence. The mapping module is used to select temporally stable edge pixels from the clean image sequence to form a stable edge set, and to perform three-dimensional spatial back projection mapping on the stable edge set based on preset camera calibration parameters to construct a three-dimensional contour line set. The comparison module is used to compare the geometric conformity of the three-dimensional contour line set with the preset architectural geometry template, and to select three-dimensional points that meet the preset geometric tolerance to form a refined structural point set. The optimization module is used to perform incremental motion structure recovery based on the refined structural point set to obtain an initial three-dimensional point cloud, and to perform bundle adjustment optimization on the three-dimensional coordinates in the initial three-dimensional point cloud that are greater than a preset reprojection error threshold to obtain an accurate three-dimensional point cloud. The filtering module is used to perform surface normal vector estimation and consistency adjustment based on the accurate 3D point cloud and the preset global reference direction to obtain an oriented normal vector field, and to perform normal vector consistency filtering on the candidate triangular facet set constructed based on the oriented normal vector field to obtain a preliminary building model. The output module is used to detect and repair mesh holes in the preliminary building model to obtain a closed surface mesh, and to perform mesh subdivision and smoothing processing based on the closed surface mesh to obtain a complete building model.

[0008] Compared with the prior art, the present invention has the following beneficial effects: (1) This invention obtains a dense motion distribution map through pixel motion analysis, and effectively identifies and removes dynamic interference areas in the image sequence based on deviation detection and dynamic region masking technology. Combined with background restoration, a clean image sequence is obtained, thereby reducing the impact of noise on three-dimensional reconstruction and improving the geometric realism and structural integrity of building models in complex environments.

[0009] (2) Based on the three-dimensional contour line set and the preset building geometric template, the present invention performs geometric conformity screening, extracts and refines the structural point set, and obtains accurate three-dimensional point cloud through observation data fusion and overall adjustment optimization, thereby realizing high-fidelity reconstruction of the main structure of the building and enhancing the accuracy and reliability of the model at the level of detail.

[0010] (3) This invention constructs a preliminary building model by estimating the surface normal and reconstructing the continuous surface, and further detects mesh voids, performs Poisson surface repair and Laplace smoothing, and finally outputs a complete and smooth three-dimensional building mesh model, which improves the continuity and visual effect of the model surface and enhances the direct usability of the three-dimensional reconstruction results. Attached Figure Description

[0011] Figure 1 This is a schematic diagram of the process of the UAV-assisted 3D building modeling method provided in the first embodiment of the present invention; Figure 2This is a schematic diagram of the structure of a UAV-assisted 3D building modeling system provided in the second embodiment of the present invention. Detailed Implementation

[0012] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0013] Reference Figure 1 The first embodiment of the present invention provides a method for building 3D modeling based on unmanned aerial vehicle (UAV) assistance, including the following steps: S11, acquire building image sequence; S12, Based on the building image sequence, the pixel displacement is calculated using the optical flow method, and the pixel displacement is analyzed into motion amplitude and motion direction to generate a dense motion distribution map; S13, Based on the dense motion distribution map, candidate motion point sets are selected based on motion amplitude, and directional consistency analysis is performed on the candidate motion point sets to mark the dynamic region mask. S14, remove the dynamic pixels marked by the dynamic region mask from the building image sequence to obtain a static background pixel set, and fill the missing regions in the static background pixel set with bilinear interpolation to obtain a clean image sequence. S15, Select temporally stable edge pixels from the clean image sequence to form a stable edge set, and perform three-dimensional spatial back projection mapping on the stable edge set based on preset camera calibration parameters to construct a three-dimensional contour line set; S16. Based on the set of three-dimensional contour lines and the preset architectural geometry template, a geometric conformity comparison is performed, and three-dimensional points that meet the preset geometric tolerance are selected to form a refined set of structural points. S17, Incremental motion structure recovery is performed based on the refined structural point set to obtain an initial three-dimensional point cloud, and bundle adjustment optimization is performed on the three-dimensional coordinates in the initial three-dimensional point cloud that are greater than a preset reprojection error threshold to obtain an accurate three-dimensional point cloud. S18. Based on the precise 3D point cloud and the preset global reference direction, surface normal vector estimation and consistency adjustment are performed to obtain an oriented normal vector field. The candidate triangular facet set constructed based on the oriented normal vector field is then screened for normal vector consistency to obtain a preliminary building model. S19, detect and repair mesh holes in the preliminary building model to obtain a closed surface mesh, and perform mesh subdivision and smoothing processing based on the closed surface mesh to obtain a complete building model.

[0014] In step S11, a sequence of building images is acquired.

[0015] Specifically, an architectural image sequence refers to a collection of digital images arranged in time-stamp order, captured by a visible light camera mounted on a drone during its orbital flight around a building. Each frame in this architectural image sequence contains color pixel information of the architectural scene, and there is a high degree of overlap between adjacent frames. The image data also includes intrinsic and extrinsic parameters of the camera. The intrinsic parameters include focal length and principal point coordinates, while the extrinsic parameters include the camera's position and orientation in space at the time each frame was captured.

[0016] The continuous inter-frame overlap and temporal sequence of the building image sequence provide the necessary data foundation for pixel motion analysis in step S12. The color pixel information in the sequence serves as input data for subsequent grayscale processing, while the accompanying camera parameters provide the geometric projection relationship for the reconstruction from the two-dimensional edge feature set to the three-dimensional contour line set in step S15. The acquisition of the image sequence ensures that subsequent steps can be processed based on a consistent observation data source.

[0017] In step S12, the pixel displacement is calculated using the optical flow method based on the building image sequence, and the pixel displacement is analyzed into motion amplitude and motion direction to generate a dense motion distribution map.

[0018] In one implementation, a grayscale image set is obtained by performing a color space conversion through grayscale processing based on the architectural image sequence. The difference in grayscale value between each corresponding pixel in two adjacent frames of the grayscale image set is calculated to obtain the grayscale gradient field; Based on the grayscale gradient field, an optical flow calculation equation is constructed by a preset constant brightness constraint and a preset local smoothing constraint. The pixel displacement field is obtained by numerically solving the optical flow calculation equation using the Lucas-Kanade method. The displacement vector of each pixel in the pixel displacement field is transformed into polar coordinates to obtain a displacement analytical field that includes the motion amplitude and motion direction; The motion amplitude and direction of each pixel in the displacement analysis field are mapped to color depth and direction arrow, respectively, to obtain a dense motion distribution map.

[0019] Specifically, in step S12, pixel motion analysis is performed based on the building image sequence. The input data is the building image sequence containing color pixel information and camera parameters obtained in step S11. First, the sequence is converted to grayscale. The brightness values ​​of the red, green, and blue channels of each pixel in each frame of the color image are weighted and summed according to preset weighting coefficients (e.g., 0.299, 0.587, 0.114) to calculate a grayscale value, thereby converting the entire building image sequence into a grayscale image set containing only single-channel brightness information. This processing aims to reduce the data dimensionality and complexity of subsequent calculations.

[0020] The difference in grayscale values ​​between corresponding pixels in two adjacent frames of the grayscale image set is calculated. Specifically, the grayscale value of a pixel in the later frame is subtracted from the grayscale value of the same pixel in the previous frame. This difference forms a two-dimensional matrix of the same size as the image, called the grayscale gradient field. This grayscale gradient field characterizes the intensity of pixel brightness changes between adjacent time points.

[0021] The optical flow calculation equation is constructed based on the gray-level gradient field. This equation is based on two physical constraints. The first is a constant brightness constraint, which assumes that the gray-level value of the projection position of the same object point remains unchanged in adjacent frames. The second is a local smoothness constraint, which assumes that adjacent pixels have similar motion vectors. These two constraints are expressed mathematically as follows: the constant brightness constraint is represented by the dot product of the spatial gradient of the pixel gray level and the displacement vector, plus the partial derivative of the gray level with respect to time, the sum of which is zero. The partial derivative of the gray level with respect to time is characterized by the gray-level gradient field calculated above. The local smoothness constraint is measured by the magnitude of the change in the displacement vector within its local neighborhood (usually using the sum of the squares of its spatial gradients). These two constraint terms are linearly combined using a preset weighting coefficient λ to form a least-squares optimization objective function for each pixel displacement vector, which is the optical flow calculation equation. The weighting coefficient λ is used to balance the relative importance of the two constraints. Its value is determined by analyzing historical image sequences. Specifically, multiple drone flight image sequences containing typical building backgrounds but without persistent large-scale dynamic interference are collected as a calibration set. For each sequence in the calibration set, the median reprojection error that only satisfies the constant brightness constraint and the median displacement field change that only satisfies the local smoothness constraint are calculated. The initial value of λ is set to the average of the ratio of the median reprojection error to the median displacement field change of all calibration sequences. This method aims to balance the relative importance of the two constraints in typical building shooting scenarios.

[0022] The constructed optical flow calculation equations are numerically solved using the Lucas-Kanade method. This method assumes that all pixels share the same displacement vector within a small neighborhood window of each pixel. The constant brightness constraint equations corresponding to all pixels within the window are simultaneously solved to form an overdetermined linear system of equations. This system is then solved using the least squares method to obtain the optimal displacement vector estimate for the center pixel. By traversing all pixels in the entire image, a complete pixel displacement field is finally obtained, where each pixel corresponds to a two-dimensional displacement vector.

[0023] The displacement vector of each pixel in the pixel displacement field is transformed into polar coordinates. For each displacement vector, its magnitude is calculated to obtain the motion amplitude, and its angle with the positive horizontal axis is calculated to obtain the motion direction. The motion amplitudes and directions of all pixels together constitute the analytical displacement field.

[0024] The displacement analysis field is visualized and encoded to generate a dense motion distribution map. The encoding rule is to linearly map the motion amplitude of each pixel to the depth value of a color (e.g., grayscale or thermal spectrum), with a larger amplitude resulting in a darker or warmer color. Simultaneously, the motion direction of each pixel is represented by a short line segment or arrow superimposed on the pixel location, with the arrow pointing in the same direction as the motion. This dense motion distribution map serves as input data for step S13. Its significance lies in visually and quantitatively revealing the motion of all pixels in the entire image sequence, providing a data foundation for subsequently distinguishing between the global background flow field caused by the UAV's own motion and local abnormal motion caused by the motion of independent objects in the scene. A motion amplitude threshold (e.g., selecting the 95th percentile of motion amplitude values ​​in the historical displacement analysis field) can be used in step S13 to initially screen candidate pixels with significant motion.

[0025] In step S13, based on the dense motion distribution map, a set of candidate motion points is selected based on the motion amplitude, and a directional consistency analysis is performed on the set of candidate motion points to mark the dynamic region mask.

[0026] In one implementation, pixels in the dense motion distribution map whose motion amplitude is greater than a preset motion amplitude threshold are selected to obtain a candidate motion point set. Based on the dense motion distribution map, a global motion model is obtained by estimating the motion model using a random sampling consensus algorithm. Based on the candidate motion point set and the global motion model, background motion compensation is performed using a model prediction control algorithm to obtain a compensation vector set; Calculate the angle between the motion direction of each pixel in the candidate motion point set and the corresponding vector direction in the compensation vector set to obtain the direction difference value; The regions containing pixels whose directional difference values ​​are greater than a preset difference threshold are marked to obtain a dynamic region mask.

[0027] Specifically, the input data for step S13 is the dense motion distribution map obtained in step S12. This map encodes the motion amplitude and direction information of each location in pixels. First, the dense motion distribution map is filtered based on a preset motion amplitude threshold, and all pixels with motion amplitude values ​​greater than the threshold are extracted to form a candidate motion point set. The motion amplitude threshold is determined by analyzing dense motion distribution maps generated from image sequences of historical clean architectural scenes (i.e., without dynamic interference objects). Specifically, the set of motion amplitude values ​​of each pixel in all historical distribution maps is calculated, and the 95th percentile of the values ​​in this set, arranged in ascending order, is taken as the threshold.

[0028] Based on the dense motion distribution map, a random sampling consensus algorithm is used to estimate the global motion model. This algorithm first randomly selects a minimum subset (e.g., for an affine motion model, the minimum subset consists of three sets of non-collinear displacement vectors) from all pixel displacement data represented by the dense motion distribution map. An initial global motion model parameter is then fitted using this subset. This initial model is then used to test all displacement data, calculating the deviation between each pixel displacement data and the model's predicted displacement. Data with deviations less than a preset tolerance error are marked as inliers. This process of random sampling, model fitting, and inlier counting is repeated multiple times (e.g., one thousand times), and finally, the model parameter fitted in the iteration with the largest number of inliers is selected as the final global motion model. This model describes the background motion patterns affecting the entire image caused by the UAV's own flight or lens translation and scaling.

[0029] Based on the candidate motion point set and the global motion model, background motion compensation is performed using a model prediction control algorithm. For each pixel in the candidate motion point set, the displacement vector that the point should produce under pure background motion is predicted using the estimated global motion model. This predicted vector is the background motion compensation vector for that point. The compensation vectors of all candidate points constitute the compensation vector set.

[0030] The directional difference between the motion direction of each pixel in the candidate motion point set and the corresponding vector direction in the compensation vector set is calculated. For each candidate point, the actual motion direction angle recorded in the dense motion distribution map is subtracted from the calculated direction angle of the background motion compensation vector to obtain the absolute value of the angle difference. Since the angle is periodic, if the absolute value is greater than 180 degrees, it is subtracted from 360 degrees, and the final directional difference value is between 0 and 180 degrees.

[0031] Based on a preset directional difference threshold, all candidate points are judged, and pixels with directional difference values ​​greater than the threshold are marked. The difference threshold is also determined based on historical data analysis. The set of directional difference values ​​of all pixels in a historical clean scene after background compensation is calculated, and the 90th percentile of the values ​​in this set, sorted in ascending order, is taken as the threshold. The connected regions containing all marked pixels are extracted and merged to form a binary image mask, where the marked region has a value of one and the remaining regions have a value of zero. This mask is the dynamic region mask. This mask serves as the input for step S14, and its significance lies in accurately identifying pixel regions in the image sequence that are significantly inconsistent with the global background motion pattern. These regions correspond to independent moving objects in the scene, providing direct spatial location information for subsequent removal of dynamic interference from the image.

[0032] In step S14, dynamic pixels marked by the dynamic region mask are removed from the building image sequence to obtain a static background pixel set, and bilinear interpolation is performed to fill the missing regions in the static background pixel set to obtain a clean image sequence.

[0033] In one implementation, the dynamic region mask is spatially aligned with the building image sequence to obtain a dynamic pixel index; Extract the pixels in the building image sequence that are not marked by the dynamic pixel index to obtain a static background pixel set; The static background pixel set is subjected to temporal median filtering to obtain a stable background pixel set; Based on the stable background pixel set, the missing regions are filled using a bilinear interpolation algorithm to obtain a clean image sequence.

[0034] Specifically, the input data for step S14 consists of the dynamic region mask generated in step S13 and the building image sequence obtained in step S11. The dynamic region mask is a binary image, where a pixel value of one corresponds to a region determined to be in abnormal motion, and a pixel value of zero corresponds to a static background region. The building image sequence is raw data containing multiple frames of color images.

[0035] The dynamic region mask is spatially aligned with the architectural image sequence. Since the mask is derived from image sequence analysis and its size is identical to each frame of the image, the alignment operation is a direct pixel position mapping. For each frame in the image sequence, the pixel coordinates with a value of one in the dynamic region mask are recorded to form the dynamic pixel index for that frame. This index is a list that lists the coordinates of all pixels in the image that need to be removed.

[0036] Extract the static background pixel set based on the dynamic pixel index. For each frame in the architectural image sequence, iterate through all pixel positions. If the coordinates of a position are not in the dynamic pixel index list corresponding to that frame, the pixel is considered a static background pixel. Collect the RGB color value and coordinate information of that pixel. Perform this operation on all frames in the sequence. The final result is a set containing data from multiple frames, but consisting only of the static pixels of each frame, i.e., the static background pixel set. In this set, the same spatial coordinate point may have multiple pixel values ​​in different frames (when it is not marked as dynamic), and may also have no value in some frames (when it is marked as dynamic).

[0037] Temporal median filtering is performed on the static background pixel set to obtain a stable background pixel set. For each spatial coordinate location on the image plane, all pixel values ​​(from different frames) it possesses in the static background pixel set are checked. This is done separately for each color channel (red, green, blue), sorting these pixel values ​​by numerical value and taking the median. This median is assigned to this coordinate location as its stable background color value. If a coordinate location has no pixel value in the static background pixel set (i.e., this location is marked as dynamic in all frames), then the value for this location remains empty in this step. After this processing, each spatial coordinate location obtains an RGB color value that best represents the static background, obtained through temporal statistics. These values ​​constitute the stable background pixel set, which is essentially a complete background image that eliminates transient dynamic interference, but may contain empty areas.

[0038] Based on a stable background pixel set, a bilinear interpolation algorithm is used to fill in missing regions to obtain a clean image sequence. For the coordinates of a missing value in the stable background pixel set, the four nearest valid pixels (i.e., points whose values ​​are not missing) are found in its surrounding neighborhood. Using the coordinates and color values ​​of these four points, interpolation calculations are performed on the red, green, and blue channels respectively. The interpolation weights are determined by the horizontal and vertical distances between the point to be filled and its four neighbors. The color value of the point to be filled is equal to the weighted sum of the color values ​​of the four neighbors, weighted according to the reciprocal of their relative distances. After filling all missing positions, a complete background image without any gaps is obtained. This image is output as the content of each frame of the clean image sequence. This means that the clean image sequence is a set of visually consistent images containing only static architectural backgrounds, and each frame originates from the temporal statistics and spatial restoration of static background information from multiple original images.

[0039] This sequence serves as the input for step S15. Its significance lies in completely removing visual interference and erroneous textures caused by dynamic objects such as pedestrians and vehicles, providing a stable and consistent image foundation for subsequent geometric feature extraction based on edges and structures, thereby ensuring that the reconstructed 3D contours only reflect the geometric characteristics of the building itself.

[0040] In step S15, stable edge pixels with stable temporal sequence are selected from the clean image sequence to form a stable edge set, and the stable edge set is back-projected into three-dimensional space based on preset camera calibration parameters to construct a three-dimensional contour line set.

[0041] In one implementation, edge features are extracted from the clean image sequence using the Canny edge detection algorithm to obtain an edge feature set. Based on the edge feature set, edge key points are extracted using the Harris corner detection algorithm, and the Lucas-Kanade method is used to perform temporal position tracking on the edge key points to obtain a set of motion trajectories. The trajectory points in the motion trajectory set whose standard deviation of the position sequence is less than a preset position standard deviation threshold are marked to obtain a stable point set, and stable edge pixels associated with the stable point set are extracted from the edge feature set to form a stable edge set; Based on the stable edge set and preset camera calibration parameters, a three-dimensional spatial coordinate system mapping is performed through back projection calculation to obtain a three-dimensional line segment set. Based on the three-dimensional line segment set, a three-dimensional structural line set is obtained by structural integration through preset collinearity constraints and preset coplanarity constraints. Based on the three-dimensional structural line set, a spatial topology connection graph is constructed, and based on the spatial topology connection graph, loop segments are extracted to obtain a three-dimensional contour line set.

[0042] Specifically, the input data for step S15 is the clean image sequence output from step S14. This sequence consists of multiple identical frames containing only static architectural backgrounds. First, the Canny edge detection algorithm is applied to each frame in the clean image sequence. This algorithm first converts the input color image to a grayscale image, then uses a two-dimensional Gaussian filter to convolve the grayscale image to smooth noise. The filter size and standard deviation parameters are preset based on common image resolutions. Next, the magnitude and direction of the brightness gradient at each pixel in the image are calculated in the horizontal and vertical directions. The horizontal gradient is obtained by convolving the image with a Sobel operator, and the vertical gradient is obtained by convolving with another Sobel operator. The gradient magnitude of each pixel is calculated by taking the square root of the sum of the squares of the horizontal and vertical gradients, and the gradient direction is calculated by the two-parameter arctangent function of the ratio of the vertical to the horizontal gradient.

[0043] Non-maximum suppression is performed. For each pixel, interpolation is performed along its gradient direction to obtain the gradient magnitude of two adjacent pixels in that direction. If the gradient magnitude of the current pixel is not less than the gradient magnitudes of these two adjacent pixels, it is retained; otherwise, the suppression is zero. Finally, a dual-threshold hysteresis connection is performed, setting a high threshold and a low threshold. The high threshold is determined by analyzing the set of edge gradient magnitudes from historical clean images and taking the 70th percentile of the values ​​in the set after ascending order. The low threshold is set to 0.4 times the high threshold. Pixels with gradient magnitudes greater than the high threshold are marked as strong edge points, points with gradient magnitudes lower than the low threshold are suppressed, and points in between are weak edge points. Finally, only those weak edge points connected to the strong edge points are retained. The binary image formed by all the retained pixels is the edge feature set.

[0044] Edge keypoints are extracted from the edge feature set. The Harris corner detection algorithm is used to calculate the corner response function value for each pixel in the image. This calculation is based on the sum of squares of the image grayscale changes within a small window when the pixel makes a tiny translation. This sum of squares can be approximated as a quadratic function whose eigenvalues ​​are related to the corner response. Specifically, the gradients of the image in the x and y directions are calculated, and the autocorrelation matrix of the gradients is constructed. Then, the determinant of this matrix is ​​subtracted from the square of its trace multiplied by an empirical coefficient (e.g., 0.04), and the result is used as the corner response value. Non-maximum suppression and a threshold are used to select the points with the largest local response values ​​as edge keypoints. This threshold is determined by analyzing the Harris response value set in historical building images and taking the 80th percentile.

[0045] The Lucas-Kanade optical flow method is used to perform temporal position tracking of the extracted edge keypoints across each frame of the clean image sequence. Since the content is identical across frames of the clean image sequence, ideally, the keypoint positions should remain unchanged. For each keypoint, its position in the next frame is solved using an optical flow equation similar to that in step S12 within its neighborhood window, thus obtaining the position sequence of the keypoint across all frames, i.e., a motion trajectory. The trajectories of all keypoints constitute a set of motion trajectories.

[0046] The stability of each trajectory in the motion trajectory set is evaluated. The standard deviation of the coordinates of all points on each trajectory is calculated. This is done by first calculating the mean of the x and y coordinates of each point, then calculating the sum of the squares of the differences between each coordinate value and the mean, dividing by the number of coordinates, and taking the square root to obtain the standard deviation of the trajectory position sequence. A standard deviation threshold is set, determined by analyzing the set of standard deviations of keypoint trajectory positions in historical clean image sequences and taking the 75th percentile. Trajectory points with a standard deviation less than this threshold are marked as stable points. From the initially extracted edge feature set, all pixel line segments connecting these stable points are extracted; this set of line segments is called the stable edge set.

[0047] 3D spatial location reconstruction is performed based on a stable edge set and preset camera calibration parameters. For each 2D pixel line segment in the stable edge set, its two endpoints are first normalized using the camera intrinsic matrix, and then its coordinates in 3D space are calculated using triangulation methods in conjunction with the camera extrinsic matrix corresponding to that frame of image. Specifically, using the observed positions of the line segment in at least two frames of imagery, its corresponding projection ray equation is constructed, and the closest point between the two rays in 3D space is solved, with this closest point being taken as the endpoint of the 3D line segment. This process is repeated for all edge line segments, resulting in a set consisting of line segment endpoint pairs in 3D space, called the 3D line segment set.

[0048] Geometric constraints are applied to a set of 3D line segments for structural integration. Collinearity constraint merges multiple short line segments that are close to each other and nearly parallel in direction into a longer straight line segment. The criteria are that the Euclidean distance between the endpoints of the line segments is less than a merging threshold and the angle between the direction vectors of the line segments is less than an angular threshold. Coplanarity constraint adjusts multiple line segments located in the same plane to be strictly coplanar. Principal component analysis is used to fit the plane containing the endpoints of these line segments, and the endpoints of all line segments are projected onto this fitted plane. The resulting set of 3D line segments after integration is called a 3D structural line set.

[0049] A spatial topological connection graph is constructed based on the 3D structural line set. The nodes of this graph are the endpoints of 3D line segments. If the Euclidean distance between two endpoints is less than a connection threshold (determined by analyzing the 85th percentile of historical 3D line segment endpoint distances), an edge is established between these two nodes. Based on this graph, a graph search algorithm is used to find all closed loop paths, i.e., paths that start from a node, follow edges, and ultimately return to that node without revisiting the edges. The 3D line segments corresponding to these loops are extracted to form the final 3D contour line set.

[0050] From the clean imagery after removing dynamic interference, the geometric representation of the building's main structural edge contour in three-dimensional space is accurately recovered by combining temporal stability verification and three-dimensional geometric constraints. The three-dimensional contour set eliminates unreliable features caused by image noise or texture, focusing on the building's inherent structural edges. It provides high-fidelity input data for the refinement and screening based on preset geometric templates in step S16, and is the key geometric foundation for achieving high-precision three-dimensional building reconstruction.

[0051] In step S16, the geometric conformity is compared with the three-dimensional contour line set and the preset building geometry template, and three-dimensional points that meet the preset geometric tolerance are selected to form a refined structural point set.

[0052] In one implementation, a geometric conformity matrix is ​​generated by calculating the Euclidean distance deviation and included angle deviation between each line segment of the three-dimensional contour line set and the corresponding line segment of the preset architectural geometry template. Mark the three-dimensional contour line segment regions corresponding to the elements in the geometric conformity matrix whose Euclidean distance deviation is less than a preset distance deviation and whose included angle deviation is less than a preset included angle deviation threshold, and obtain the effective structural region; Extract the points corresponding to the effective structural regions from the three-dimensional contour set to obtain a refined structural point set.

[0053] Specifically, the input data for step S16 consists of the three-dimensional contour line set obtained in step S15 and a preset architectural geometry template. The architectural geometry template is a predefined data structure used to describe the spatial geometric patterns of typical architectural structures, such as common planes, straight edges, right-angle connections, and symmetry.

[0054] This template is built upon a database of 3D models of historical buildings. The models in the database are sourced from publicly available standard architectural model libraries or historical project data that has been precisely surveyed. The construction process begins by normalizing the orientation of all models in the database to ensure alignment of their main directions. Then, an octree-based spatial partitioning method is used to divide the model space into voxel meshes of different levels. At each voxel level, the types of geometric primitives (such as planar pieces and line segments) and their frequencies are statistically analyzed across all models. Finally, combinations of geometric primitives that appear frequently across multiple levels and have stable spatial relationships are extracted, parameterized (e.g., the normal vector and distance from the origin of a plane, the start and end coordinates and direction vector of a line segment), and stored to form the architectural geometric template. When used, the parameterized geometric primitives of this template are instantiated as specific 3D reference line segments.

[0055] Specifically, the preset architectural geometry template is constructed by quantitatively analyzing a diverse database of standard architectural 3D models. After posture normalization, an octree is used to divide the model space. For each voxel unit, all geometric primitives passing through or located in that unit in the model are counted, such as line segments and planes. A lower limit for the frequency of occurrence is set, for example, all of them appear in 30% of the database models, and primitives that meet this frequency are marked as candidates. For candidate primitives, their spatial relationships are further analyzed. For line primitives, it is checked whether they meet the parallel or perpendicular relationship with other nearby candidate primitives (the deviation of the included angle from the standard angle is less than 5 degrees). For plane primitives, it is checked whether the angle between their normal and the vertical direction is less than 10 degrees (representing the wall surface). The combination of primitives that simultaneously meet the frequency condition and the spatial relationship condition, along with their parameters (direction, position), are stored to form a template. The distance deviation threshold and included angle deviation threshold in the geometric conformity comparison are set according to the statistical distribution of the average fitting error between all template primitives in the template library and their corresponding parts in the original model, such as taking the average value plus twice the standard deviation.

[0056] Geometric conformity screening is performed based on the set of 3D contour lines and the preset architectural geometry template. First, two types of deviations are calculated between each 3D line segment in the 3D contour line set and each reference line segment in the architectural geometry template. The first is the Euclidean distance deviation, calculated by taking the midpoints of both the 3D line segment and the reference line segment and calculating the straight-line distance between these two points. The second is the angle deviation, calculated by obtaining the direction vectors (unit vectors pointing from the starting point to the ending point) of both the 3D line segment and the reference line segment, calculating the cosine of the angle between these two direction vectors, and then using the inverse cosine function to obtain the angle value.

[0057] After performing the above calculations on all line segment pairs, a feature matching degree matrix is ​​obtained. The rows of this matrix correspond to the line segments in the 3D contour line set, and the columns correspond to the reference line segments in the architectural geometry template. Each element in the matrix is ​​a tuple containing the Euclidean distance deviation and the included angle deviation between the line segment in that row and the line segment in that column.

[0058] The feature matching matrix is ​​filtered based on preset distance deviation thresholds and angle deviation thresholds. The distance deviation threshold is determined by analyzing the set of distance deviation values ​​of correctly matched line segment pairs in historical successful reconstruction cases, and taking the 90th percentile of this set of values ​​after sorting in ascending order. The angle deviation threshold is determined by analyzing the angle deviation values ​​in the same set and taking the 90th percentile. In the feature matching matrix, matrix elements that simultaneously satisfy both the Euclidean distance deviation and the angle deviation threshold are marked. Each marked matrix element's corresponding row index points to a specific line segment in the 3D contour line set.

[0059] The set of three-dimensional contour lines corresponding to all marked elements, and the three-dimensional spatial region covered by it, is defined as the effective structural region. Finally, from the original set of three-dimensional contour lines, the endpoints (i.e., start and end points) of all line segments located within the effective structural region are extracted, and the coordinates of these three-dimensional points are collected to form the refined structural point set.

[0060] This step compares the reconstructed 3D contour with prior architectural geometry knowledge to filter out geometric elements that conform to typical architectural structural rules. This effectively filters out contour segments caused by reconstruction errors, residual noise, or atypical structures, ensuring that the refined structural point set is highly concentrated on the main and regular structural features of the building itself. This step provides cleaner and more reliable 3D spatial point data for subsequent observation data fusion and overall adjustment optimization, and is a key step in improving the geometric regularity and accuracy of the final architectural model.

[0061] In step S17, incremental motion structure recovery is performed based on the refined structure point set to obtain an initial three-dimensional point cloud, and bundle adjustment optimization is performed on the three-dimensional coordinates in the initial three-dimensional point cloud that are greater than a preset reprojection error threshold to obtain an accurate three-dimensional point cloud.

[0062] In one implementation, based on the refined structure point set, the camera pose and 3D point cloud are estimated by an incremental motion recovery structure algorithm to obtain the initial camera pose and an initial 3D point cloud containing multiple frames of eponymous feature points. Calculate the reprojection error of the same feature points in multiple frames of the initial 3D point cloud to obtain the reprojection error set; The three-dimensional coordinates corresponding to the reprojection errors in the reprojection error set that are greater than a preset reprojection error threshold are selected, and the three-dimensional coordinate deviation is calculated based on the selection results to obtain the coordinate deviation vector. Based on the coordinate deviation vector and the initial camera pose, global optimization is performed using the bundle adjustment algorithm to obtain the optimized camera pose; Based on the optimized camera pose, forward intersection calculation is performed on the refined structural point set to obtain an accurate 3D point cloud.

[0063] Specifically, the input data for step S17 is the refined structure point set obtained in step S16. Each point in this point set is associated with its two-dimensional pixel coordinate observation in multiple frames of the clean image sequence, that is, one three-dimensional point corresponds to multiple two-dimensional observation points.

[0064] Based on the refined structure point set, camera pose and 3D point cloud estimation are performed using an incremental structure-of-motion (SOMO) algorithm. The process begins by selecting two initial frames from a clean image sequence. Based on the 2D feature matching point pairs corresponding to the refined structure points in these two frames, the fundamental matrix is ​​calculated and decomposed to obtain the initial relative pose of the camera in these two frames. Subsequently, the initial 3D point coordinates corresponding to these matching point pairs are calculated using triangulation. Next, the algorithm incrementally adds new image frames. For each new image frame, a perspective N-point problem is solved, and the new camera pose is estimated using the reconstructed 3D points and their 2D observations in the new image. Then, based on the newly estimated camera pose, the existing 3D points are triangulated and expanded, and new 3D points are reconstructed using the new feature matching pairs appearing in the new image. This process is iterated until the camera poses of all image frames are estimated, generating an initial 3D point cloud supported by feature points of the same name from multiple frames and the corresponding initial camera pose.

[0065] Calculate the reprojection error for each 3D point in the initial 3D point cloud. For each 3D point, using its associated initial camera pose, reproject the 3D point coordinates onto each image plane that observed it, obtaining a 2D projected point coordinate. Calculate the Euclidean distance between this projected point coordinate and its actual 2D observation point coordinate in the image; this distance is the reprojection error of that point in that frame. The reprojection errors of a 3D point across all observation frames constitute a set, and the sum of all 3D point error sets forms the reprojection error set.

[0066] Error terms exceeding a preset reprojection error threshold are filtered from the reprojection error set. This reprojection error threshold is determined by analyzing historical successful 3D reconstruction project datasets. Specifically, the final optimized reprojection error values ​​of all 3D points from these projects are collected, and the 95th percentile of this set of values, sorted in ascending order, is used as the threshold. For each filtered reprojection error exceeding the threshold, it corresponds to an observation of a specific 3D point in a specific frame of an image. Based on these anomalous observations, the deviation estimate of the associated 3D point coordinates, i.e., the coordinate deviation vector, is calculated. The calculation method typically involves back-calculating based on the magnitude and direction of the error from the observation, combined with the camera pose, to obtain a position offset vector in 3D space.

[0067] Based on the coordinate deviation vector and the initial camera pose, global optimization is performed using the bundle adjustment algorithm. Bundle adjustment constructs a nonlinear least squares optimization problem that includes all camera pose parameters and all 3D point coordinate parameters. Its objective function is the sum of squares of the reprojection errors of all 3D points across all image frames, while the coordinate deviation vector calculated in the previous steps is added to the objective function as part of the prior constraints on the positions of certain 3D points. The Levenberg-Marquardt iterative algorithm is used to solve this optimization problem, minimizing the overall objective function value by adjusting the camera pose parameters and 3D point coordinate parameters. After the optimization process converges, a more accurate set of camera pose parameters is obtained, referred to as the optimized camera pose.

[0068] Based on the optimized camera pose, forward intersection calculations are performed on the refined structure point set. Using the optimized, more accurate camera pose parameters, triangulation calculations are re-performed for all 2D observations associated with each point in the refined structure point set. This process solves for the optimal intersection point of multiple rays (from different images) back-projected into 3D space, thus obtaining a more precise coordinate position of the point in 3D space. After performing this calculation on all points, the resulting 3D point set is the accurate 3D point cloud.

[0069] This step significantly improves the spatial geometric accuracy of the 3D point cloud by fusing observation data from multiple image frames and performing rigorous overall adjustment optimization. It optimizes not only the coordinates of the 3D points but also the position and attitude parameters of the camera, achieving optimal geometric consistency across the entire reconstruction system. The accurate 3D point cloud eliminates errors introduced by inaccurate initial matching and pose estimation, providing a high-precision geometric foundation for subsequent surface reconstruction steps and serving as a core guarantee for achieving high-fidelity 3D architectural modeling.

[0070] In step S18, surface normal vector estimation and consistency adjustment are performed based on the accurate 3D point cloud and the preset global reference direction to obtain an oriented normal vector field. Then, the candidate triangular facet set constructed based on the oriented normal vector field is screened for normal vector consistency to obtain a preliminary building model.

[0071] In one implementation, based on the precise 3D point cloud, the normal vector is estimated using principal component analysis to obtain an initial normal vector set. Based on the initial normal vector set, the normal vector direction is adjusted for consistency through a preset global reference direction to obtain a directional normal vector field; Based on the directional normal vector field, a triangular mesh is constructed using the Delaunay triangulation algorithm to obtain a triangular mesh set; The triangular mesh set is back-projected into three-dimensional space to obtain a set of candidate triangular facets; From the set of candidate triangular faces, select all candidate triangular faces whose angle with the normal vector of the corresponding center point is less than a preset angle threshold to obtain a preliminary building model.

[0072] Specifically, the input data for step S18 is the precise 3D point cloud obtained in step S17, which consists of a series of points with precise 3D coordinates. First, based on the precise 3D point cloud, normal vector estimation is performed using principal component analysis. For each point in the point cloud, all other points within a certain radius of its neighborhood are searched. This neighborhood radius is determined by analyzing the distribution of local point density in historical point cloud data; specifically, the 80th percentile of the distance set between all points and their nearest neighbors, arranged in ascending order, is taken as the radius value. Then, a covariance matrix is ​​constructed using the 3D coordinates of all points within this neighborhood, and the eigenvalues ​​and eigenvectors of this covariance matrix are calculated. The eigenvector corresponding to the smallest eigenvalue represents the opposite direction of the most dispersed distribution of the neighborhood point set, and is initially estimated as the normal vector direction of the center point. The normal vectors of all points constitute the initial normal vector set.

[0073] In one implementation, the neighborhood radius is determined using an adaptive method. For each point in the point cloud, the average distance to its K nearest neighbors is first calculated, for example, K=10. This average distance is used as an initial estimate of the neighborhood radius of that point. To ensure that the neighborhood contains enough points for stable principal component analysis, a minimum number of points is set, such as 15 points. If the number of neighborhood points found according to the initial radius is less than this minimum, the radius is gradually expanded, for example, by 20% each time, until the number of points is met or a preset maximum radius limit is reached. This method can dynamically adjust the neighborhood size according to the local point density.

[0074] Based on the initial normal vector set, the normal vector directions are adjusted for consistency using a preset global reference direction. This preset global reference direction is determined as follows: After obtaining an accurate 3D point cloud, the distribution range of the entire point cloud along all three coordinate axes is calculated, and the coordinate axis direction with the largest distribution range (typically perpendicular to the ground) is defined as the global reference direction. For each normal vector in the initial normal vector set, the cosine of its angle with the global reference direction is calculated; if this value is negative, the normal vector is inverted. After this adjustment, all normal vectors point to approximately the same hemispherical space, thus obtaining a directional normal vector field.

[0075] Based on the oriented normal vector field, a triangular mesh is constructed using the Delaunay triangulation algorithm. This operation is performed in a two-dimensional parameter space. First, the three-dimensional point cloud is projected onto a two-dimensional plane, which is spanned by the eigenvectors corresponding to the first two largest eigenvalues ​​from the principal component analysis of the point cloud, to minimize projection distortion. On the resulting two-dimensional projected point set, Delaunay triangulation is performed, ensuring that the circumcircle of all triangles does not contain any other projected points. The final result is a triangular mesh set composed of triangles, where each triangle is defined by three two-dimensional projected point indices.

[0076] The triangular mesh set is back-projected into three-dimensional space. Based on the three vertex indices of each triangle in the two-dimensional triangular mesh, the coordinates of its three corresponding three-dimensional points in the original three-dimensional point cloud are found. These three three-dimensional points constitute a candidate triangular facet in three-dimensional space. After all triangles undergo this back-projection operation, the candidate triangular facet set is obtained.

[0077] The candidate triangular facets are filtered from the set. For each triangular facet in the set, the three-dimensional coordinates of its geometric center point are calculated, and the normal vector corresponding to that center point in the oriented normal vector field is found. Simultaneously, the normal vector of the triangular facet itself is calculated (obtained by the cross product of the two side vectors of the triangle and normalized). The cosine of the angle between these two normal vectors is calculated, and then the inverse cosine function is used to obtain the angle value. A preset angle threshold is determined by analyzing the set of angles between correct triangular facets and their center point normal vectors in the historical reconstruction model, and the 85th percentile of this set of values, sorted in ascending order, is taken as the threshold. All candidate triangular facets with angle values ​​less than this threshold are selected, and this set of facets constitutes the preliminary architectural model.

[0078] This step, based on high-precision 3D point clouds, estimates surface normal vectors with consistent orientations and constructs a continuous surface mesh using triangulation techniques. The oriented normal vector field ensures the local orientation consistency of the reconstructed surface, while filtering based on the angle between the normal vectors removes erroneous triangular faces that clearly deviate from the local surface geometry and may be due to point cloud noise or projection distortion. The preliminary architectural model is thus a 2D manifold surface expressed in triangular meshes, maintaining the continuity of the building's main surfaces, providing an accurate geometric and topological foundation for void detection and repair in step S19.

[0079] In step S19, mesh holes in the preliminary building model are detected and repaired to obtain a closed surface mesh. Based on the closed surface mesh, mesh subdivision and smoothing processing is performed to obtain a complete building model.

[0080] Based on the preliminary building model, the edge with one adjacent facet is marked as the boundary, and the boundary loop is closed according to the connection order of the vertices on the boundary to obtain a closed boundary loop; Mark the regions in the closed boundary loop whose perimeter is greater than a preset perimeter threshold and whose area is greater than a preset area threshold to obtain the voids to be filled; Based on the voids to be filled, a closed surface mesh is obtained by mesh repair using the Poisson surface reconstruction algorithm. Based on the closed surface mesh, the mesh is subdivided using the Loop subdivision surface algorithm to obtain a high-density subdivision mesh; Based on the high-density subdivided mesh, the mesh is smoothed using the Laplace smoothing algorithm to obtain a complete building model.

[0081] Specifically, the input data for step S19 is the preliminary building model obtained in step S18. This model is a three-dimensional mesh structure composed of triangular facets, with boundaries but not completely closed. First, a mesh integrity check is performed on the preliminary building model. All mesh edges in the model are traversed, and for each edge, the number of triangular facets sharing this edge is counted. If an edge is shared by only one triangular facet, then the edge is marked as a boundary. Subsequently, based on the connection relationship between vertices, the boundaries that are connected end to end are combined according to their connection order in three-dimensional space to form one or more closed, non-self-intersecting spatial polygonal loops, each loop being called a closed boundary loop.

[0082] Based on preset perimeter and area thresholds, holes requiring repair are identified from all closed boundary loops. For each closed boundary loop, its perimeter and the area of ​​the enclosed region are calculated. The perimeter is obtained by summing the three-dimensional spatial lengths of each edge on the loop. The area is calculated by first projecting the vertices of the loop onto a two-dimensional plane defined by its normal vector, and then calculating the area of ​​the projected two-dimensional polygon. Both the perimeter and area thresholds are determined by analyzing the size distribution of effective structural holes (i.e., real holes caused by missing data, not mesh noise) in the historical building mesh model. Specifically, the perimeter and area data of all manually confirmed real holes in the historical model are collected, and the 75th percentile of the perimeter data set in ascending order is taken as the perimeter threshold, and the 75th percentile of the area data set in ascending order is taken as the area threshold. The region defined by a closed boundary loop that simultaneously satisfies both the perimeter and area thresholds is marked as a hole to be filled.

[0083] For each void to be filled, a Poisson surface reconstruction algorithm is used for mesh repair. This algorithm takes the vertices and their normal vectors on the void boundary as input. First, an octree spatial structure covering the void region is constructed, and a spatial basis function is defined at each octree node. The core of the algorithm is solving a Poisson equation, whose source terms are defined by the divergence of the normal vector field at the boundary vertices. By discretizing and solving the Poisson equation, a scalar indicator function in three-dimensional space is obtained. The isosurface of this function (usually the zero-value surface) defines the smooth surface for filling the void. This isosurface is extracted to generate a triangular mesh that seamlessly connects to the watertight closed boundary. This newly generated mesh is then partially merged with the original mesh to form a closed surface mesh.

[0084] The closed surface mesh is subdivided to increase its capacity to accommodate geometric details. A Loop subdivision surface algorithm is employed, which recursively subdivides each triangular facet in the mesh. In each subdivision iteration, a new vertex is inserted at the midpoint of each edge, and the coordinates of the new vertex are calculated by weighting the positions of its adjacent old vertices; simultaneously, the coordinates of the existing old vertices are adjusted based on the positions of the newly inserted adjacent vertices. After one or more such subdivision iterations, a high-density subdivided mesh with a significantly increased number of triangular facets and a higher vertex density is obtained.

[0085] The high-density subdivided mesh is smoothed using the Laplacian smoothing algorithm. For each internal vertex (i.e., non-boundary vertex) in the mesh, the average coordinates of all its first-order adjacent vertices are calculated. Then, the current position of that vertex is shifted in the direction of this average value by a distance determined by a preset coefficient (e.g., 0.5). This process is iterated multiple times (e.g., 10 times) to make the distribution of mesh vertices in space more uniform, soften sharp edges, and make the shapes of triangular facets more regular, thus obtaining a complete architectural model with a smooth surface and well-preserved geometry.

[0086] This step automatically detects and fills in missing parts of the model surface caused by blind spots in data acquisition or interference areas removed in the early stages, and optimizes the topological quality and visual appearance of the mesh through subdivision and smoothing techniques. This makes the final output complete building model not only a geometrically closed entity, but also has a uniform mesh distribution and smooth surface transitions, significantly improving the model's direct usability and fidelity in visualization, simulation analysis, or further digital applications.

[0087] Reference Figure 2 The second embodiment of the present invention provides a UAV-assisted 3D building modeling system, comprising: The data acquisition module is used to acquire building image sequences; The optical flow module is used to calculate pixel displacement based on the building image sequence using the optical flow method, and to resolve the pixel displacement into motion amplitude and motion direction to generate a dense motion distribution map; The analysis module is used to filter out candidate motion point sets based on motion amplitude according to the dense motion distribution map, and to perform directional consistency analysis on the candidate motion point sets to mark dynamic region masks. The filling module is used to remove dynamic pixels marked by the dynamic region mask from the building image sequence to obtain a static background pixel set, and to fill the missing regions in the static background pixel set with bilinear interpolation to obtain a clean image sequence. The mapping module is used to select temporally stable edge pixels from the clean image sequence to form a stable edge set, and to perform three-dimensional spatial back projection mapping on the stable edge set based on preset camera calibration parameters to construct a three-dimensional contour line set. The comparison module is used to compare the geometric conformity of the three-dimensional contour line set with the preset architectural geometry template, and to select three-dimensional points that meet the preset geometric tolerance to form a refined structural point set. The optimization module is used to perform incremental motion structure recovery based on the refined structural point set to obtain an initial three-dimensional point cloud, and to perform bundle adjustment optimization on the three-dimensional coordinates in the initial three-dimensional point cloud that are greater than a preset reprojection error threshold to obtain an accurate three-dimensional point cloud. The filtering module is used to perform surface normal vector estimation and consistency adjustment based on the accurate 3D point cloud and the preset global reference direction to obtain an oriented normal vector field, and to perform normal vector consistency filtering on the candidate triangular facet set constructed based on the oriented normal vector field to obtain a preliminary building model. The output module is used to detect and repair mesh holes in the preliminary building model to obtain a closed surface mesh, and to perform mesh subdivision and smoothing processing based on the closed surface mesh to obtain a complete building model.

[0088] It should be noted that the UAV-assisted 3D building modeling system provided in this embodiment of the invention is used to execute all the process steps of the UAV-assisted 3D building modeling method in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0089] It should be noted that the system embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Furthermore, in the accompanying drawings of the system embodiments provided by this invention, the connection relationships between modules indicate that they have communication connections, which can be specifically implemented as one or more communication buses or signal lines. Those skilled in the art can understand and implement this without any creative effort.

[0090] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A method for building 3D modeling based on unmanned aerial vehicle assistance, characterized in that, include: Acquire architectural image sequences; Based on the building image sequence, the pixel displacement is calculated using the optical flow method, and the pixel displacement is analyzed into motion amplitude and motion direction to generate a dense motion distribution map. Based on the dense motion distribution map, candidate motion point sets are selected based on motion amplitude, and directional consistency analysis is performed on the candidate motion point sets to mark the dynamic region mask. The dynamic pixels marked by the dynamic region mask are removed from the building image sequence to obtain a static background pixel set. The missing regions in the static background pixel set are filled by bilinear interpolation to obtain a clean image sequence. Stable edge pixels with stable temporal sequence are selected from the clean image sequence to form a stable edge set. Based on the preset camera calibration parameters, the stable edge set is back-projected in three-dimensional space to construct a three-dimensional contour line set. Based on the set of three-dimensional contour lines and the preset architectural geometry template, a geometric conformity comparison is performed, and three-dimensional points that meet the preset geometric tolerance are selected to form a refined set of structural points. Incremental motion structure recovery is performed based on the refined structural point set to obtain an initial three-dimensional point cloud. Then, bundle adjustment optimization is performed on the three-dimensional coordinates in the initial three-dimensional point cloud that are greater than a preset reprojection error threshold to obtain an accurate three-dimensional point cloud. Based on the precise 3D point cloud and the preset global reference direction, surface normal vector estimation and consistency adjustment are performed to obtain an oriented normal vector field. Then, the candidate triangular facet set constructed based on the oriented normal vector field is screened for normal vector consistency to obtain a preliminary building model. Detect and repair mesh holes in the preliminary building model to obtain a closed surface mesh, and perform mesh subdivision and smoothing processing based on the closed surface mesh to obtain a complete building model.

2. The UAV-assisted 3D building modeling method according to claim 1, characterized in that, The step of calculating pixel displacement using optical flow based on the building image sequence, and resolving the pixel displacement into motion amplitude and direction to generate a dense motion distribution map includes: Based on the architectural image sequence, a grayscale image set is obtained by performing color space conversion through grayscale processing. The difference in grayscale value between each corresponding pixel in two adjacent frames of the grayscale image set is calculated to obtain the grayscale gradient field; Based on the grayscale gradient field, an optical flow calculation equation is constructed by a preset constant brightness constraint and a preset local smoothing constraint. The pixel displacement field is obtained by numerically solving the optical flow calculation equation using the Lucas-Kanade method. The displacement vector of each pixel in the pixel displacement field is transformed into polar coordinates to obtain a displacement analytical field that includes the motion amplitude and motion direction; The motion amplitude and direction of each pixel in the displacement analysis field are mapped to color depth and direction arrow, respectively, to obtain a dense motion distribution map.

3. The UAV-assisted 3D building modeling method according to claim 1, characterized in that, The step of selecting candidate motion point sets based on the motion amplitude according to the dense motion distribution map, and performing directional consistency analysis on the candidate motion point sets to mark dynamic region masks includes: Pixels whose motion amplitude is greater than a preset motion amplitude threshold are selected from the dense motion distribution map to obtain a candidate motion point set; Based on the dense motion distribution map, a global motion model is obtained by estimating the motion model using a random sampling consensus algorithm. Based on the candidate motion point set and the global motion model, background motion compensation is performed using a model prediction control algorithm to obtain a compensation vector set; Calculate the angle between the motion direction of each pixel in the candidate motion point set and the corresponding vector direction in the compensation vector set to obtain the direction difference value; The regions containing pixels whose directional difference values ​​are greater than a preset difference threshold are marked to obtain a dynamic region mask.

4. The UAV-assisted 3D building modeling method according to claim 1, characterized in that, The process of removing dynamic pixels marked by the dynamic region mask from the architectural image sequence to obtain a static background pixel set, and then performing bilinear interpolation to fill the missing regions in the static background pixel set to obtain a clean image sequence, includes: The dynamic region mask is spatially aligned with the building image sequence to obtain a dynamic pixel index; Extract the pixels in the building image sequence that are not marked by the dynamic pixel index to obtain a static background pixel set; The static background pixel set is subjected to temporal median filtering to obtain a stable background pixel set; Based on the stable background pixel set, the missing regions are filled using a bilinear interpolation algorithm to obtain a clean image sequence.

5. The UAV-assisted 3D building modeling method according to claim 1, characterized in that, The process of selecting temporally stable edge pixels from the clean image sequence to form a stable edge set, and then performing three-dimensional spatial back-projection mapping on the stable edge set based on preset camera calibration parameters to construct a three-dimensional contour line set, includes: Based on the clean image sequence, edge features are extracted using the Canny edge detection algorithm to obtain an edge feature set; Based on the edge feature set, edge key points are extracted using the Harris corner detection algorithm, and the Lucas-Kanade method is used to perform temporal position tracking on the edge key points to obtain a set of motion trajectories. The trajectory points in the motion trajectory set whose standard deviation of the position sequence is less than a preset position standard deviation threshold are marked to obtain a stable point set, and stable edge pixels associated with the stable point set are extracted from the edge feature set to form a stable edge set; Based on the stable edge set and preset camera calibration parameters, a three-dimensional spatial coordinate system mapping is performed through back projection calculation to obtain a three-dimensional line segment set. Based on the three-dimensional line segment set, a three-dimensional structural line set is obtained by structural integration through preset collinearity constraints and preset coplanarity constraints. Based on the three-dimensional structural line set, a spatial topology connection graph is constructed, and based on the spatial topology connection graph, loop segments are extracted to obtain a three-dimensional contour line set.

6. The UAV-assisted 3D building modeling method according to claim 1, characterized in that, The step of comparing the geometric conformity of the three-dimensional contour line set with the preset architectural geometry template, and selecting three-dimensional points that meet the preset geometric tolerance to form a refined structural point set includes: A geometric conformity matrix is ​​generated by calculating the Euclidean distance deviation and included angle deviation between each line segment of the three-dimensional contour line set and the corresponding line segment of the preset architectural geometry template. Mark the three-dimensional contour line segment regions corresponding to the elements in the geometric conformity matrix whose Euclidean distance deviation is less than a preset distance deviation and whose included angle deviation is less than a preset included angle deviation threshold, and obtain the effective structural region; Extract the points corresponding to the effective structural regions from the three-dimensional contour set to obtain a refined structural point set.

7. The UAV-assisted 3D building modeling method according to claim 1, characterized in that, The incremental motion structure recovery based on the refined structural point set yields an initial 3D point cloud. Bundle adjustment optimization is then performed on the 3D coordinates in the initial 3D point cloud that exceed a preset reprojection error threshold to obtain an accurate 3D point cloud. This includes: Based on the refined structure point set, the camera pose and 3D point cloud are estimated by incremental structure motion recovery algorithm to obtain the initial camera pose and the initial 3D point cloud containing feature points with the same name in multiple frames. Calculate the reprojection error of the same feature points in multiple frames of the initial 3D point cloud to obtain the reprojection error set; The three-dimensional coordinates corresponding to the reprojection errors in the reprojection error set that are greater than a preset reprojection error threshold are selected, and the three-dimensional coordinate deviation is calculated based on the selection results to obtain the coordinate deviation vector. Based on the coordinate deviation vector and the initial camera pose, global optimization is performed using the bundle adjustment algorithm to obtain the optimized camera pose; Based on the optimized camera pose, forward intersection calculation is performed on the refined structural point set to obtain an accurate 3D point cloud.

8. The UAV-assisted 3D building modeling method according to claim 1, characterized in that, The process involves estimating and adjusting the surface normal vectors based on the precise 3D point cloud and a preset global reference direction to obtain an oriented normal vector field. Then, a set of candidate triangular facets constructed based on this oriented normal vector field is screened for normal vector consistency to obtain a preliminary architectural model, including: Based on the precise 3D point cloud, the normal vectors are estimated using principal component analysis to obtain the initial normal vector set; Based on the initial normal vector set, the normal vector direction is adjusted for consistency through a preset global reference direction to obtain a directional normal vector field; Based on the directional normal vector field, a triangular mesh is constructed using the Delaunay triangulation algorithm to obtain a triangular mesh set; The triangular mesh set is back-projected into three-dimensional space to obtain a set of candidate triangular facets; From the set of candidate triangular faces, select all candidate triangular faces whose angle with the normal vector of the corresponding center point is less than a preset angle threshold to obtain a preliminary building model.

9. The UAV-assisted 3D building modeling method according to claim 1, characterized in that, The process of detecting and repairing mesh holes in the preliminary building model to obtain a closed surface mesh, and then performing mesh subdivision and smoothing based on the closed surface mesh to obtain a complete building model, includes: Based on the preliminary building model, the edge with one adjacent facet is marked as the boundary, and the boundary is closed according to the connection order of the vertices on the boundary to obtain a closed boundary loop; Mark the regions in the closed boundary loop whose perimeter is greater than a preset perimeter threshold and whose area is greater than a preset area threshold to obtain the voids to be filled; Based on the voids to be filled, a closed surface mesh is obtained by mesh repair using the Poisson surface reconstruction algorithm. Based on the closed surface mesh, the mesh is subdivided using the Loop subdivision surface algorithm to obtain a high-density subdivision mesh; Based on the high-density subdivided mesh, the mesh is smoothed using the Laplace smoothing algorithm to obtain a complete building model.

10. A UAV-assisted 3D architectural modeling system, characterized in that, include: The data acquisition module is used to acquire building image sequences; The optical flow module is used to calculate pixel displacement based on the building image sequence using the optical flow method, and to resolve the pixel displacement into motion amplitude and motion direction to generate a dense motion distribution map; The analysis module is used to filter out candidate motion point sets based on motion amplitude according to the dense motion distribution map, and to perform directional consistency analysis on the candidate motion point sets to mark dynamic region masks. The filling module is used to remove dynamic pixels marked by the dynamic region mask from the building image sequence to obtain a static background pixel set, and to fill the missing regions in the static background pixel set with bilinear interpolation to obtain a clean image sequence. The mapping module is used to select temporally stable edge pixels from the clean image sequence to form a stable edge set, and to perform three-dimensional spatial back projection mapping on the stable edge set based on preset camera calibration parameters to construct a three-dimensional contour line set. The comparison module is used to compare the geometric conformity of the three-dimensional contour line set with the preset architectural geometry template, and to select three-dimensional points that meet the preset geometric tolerance to form a refined structural point set. The optimization module is used to perform incremental motion structure recovery based on the refined structural point set to obtain an initial three-dimensional point cloud, and to perform bundle adjustment optimization on the three-dimensional coordinates in the initial three-dimensional point cloud that are greater than a preset reprojection error threshold to obtain an accurate three-dimensional point cloud. The filtering module is used to perform surface normal vector estimation and consistency adjustment based on the accurate 3D point cloud and the preset global reference direction to obtain an oriented normal vector field, and to perform normal vector consistency filtering on the candidate triangular facet set constructed based on the oriented normal vector field to obtain a preliminary building model. The output module is used to detect and repair mesh holes in the preliminary building model to obtain a closed surface mesh, and to perform mesh subdivision and smoothing processing based on the closed surface mesh to obtain a complete building model.

Images