A method for calculating building height based on 3D Gaussian technology

Abstract

The application provides a kind of aerial data building height calculation method based on 3D Gauss technology, comprising obtaining aerial photograph;Carrying out data pretreatment;Carrying out geographic coordinate conversion;Taking all aerial photographs, calculating its position in local coordinate system before conversion, applying conversion function f_best (x, y, z), make local coordinate system and target coordinate system align;The internal parameter, external parameter, sparse point cloud of converted aerial photograph are input into 3DGS algorithm, and scene reconstruction is carried out;Obtain user click position;Calculate the building height of the user click position.The aerial data building height calculation method based on 3D Gauss technology proposed by the application introduces 3DGS technology into city building three-dimensional modeling, for high-definition image data collected by unmanned aerial vehicle, while taking into account the advantages of 3DGS technology, realizes interactive browsing and automatic height measurement task.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of aerial building height calculation, and particularly relates to a method for calculating building height of aerial data based on 3D Gaussian technology. BACKGROUND

[0002] With the rise of smart city construction, three-dimensional modeling of city-level large scenes has wide application value in many fields, including city roaming, autonomous driving and navigation guide, etc. Traditional three-dimensional modeling methods, such as oblique photography modeling, express scenes by display mode, and describe scene information through elements such as geometry, texture and material. This method has a long modeling process, produces a large amount of data, and in order to achieve high fidelity, often requires manual participation in a large amount of editing work, and the data production efficiency is low.

[0003] 3D Gaussian is a new breakthrough in the field of multi-view synthesis, which represents the geometry and scene in 3D space as a large number of Gaussian distributions, and then projects and renders pixel colors on these Gaussian distributions through the parallel computing power of GPU to generate realistic 2D images. First, the colmap software is used to calculate the camera intrinsic and extrinsic parameters of each photo to generate a sparse point cloud under the local coordinate system. Then, the 3DGS technology starts from sparse points and represents the scene with a three-dimensional Gaussian distribution, which retains the ideal properties of continuous body radiation field for scene optimization, while avoiding unnecessary calculations in blank space; according to the camera intrinsic and extrinsic parameters, the 3D Gaussian distribution is projected onto a two-dimensional plane, and then the 3DGS uses a fast visibility-aware rendering algorithm to calculate the color of each pixel, which allows real-time rendering, supports anisotropic splatting, and accelerates training. According to the loss calculated from the inference image and the real image, the parameters of the 3DGS are optimized through back propagation; at the same time, the three-dimensional Gaussian function is optimized by interleaving or density control, especially the anisotropic covariance, to achieve accurate representation of the scene.

[0004] However, using 3D Gaussian technology for building height measurement of aerial data has the following problems: the previous height measurement method has the problems of low precision, low efficiency, high cost and great difficulty, which is difficult to promote; the model established by the 3DGS algorithm has high fidelity and simplifies the modeling steps, but the model established by the 3DGS algorithm is only a visual model similar to the real world, but the model cannot directly recover the real spatial distance. The 3DGS modeling algorithm based on depth map has the problems of unclear boundary under the current view and inaccurate depth acquisition at the boundary, which also directly affects the accuracy of spatial distance calculation.

[0005] The invention patent application with the application number CN119762983A discloses a building height change identification method, device and equipment based on remote sensing space-time information and a storage medium, comprising: obtaining a building bottom contour using a vector map; obtaining building top contours at different times using multi-temporal high-resolution optical remote sensing images; obtaining building top-bottom contour projection displacement according to the Euclidean distance between the geometric center point of the building bottom contour and the geometric center point of the building top contour; calculating multi-building relative projection displacement ratios and single-building absolute projection displacement ratios of different buildings according to the building top-bottom contour projection displacement at different times; calculating satellite inclination tangent ratios according to the satellite observation angles at different times; and identifying height changes of the buildings according to the multi-building relative projection displacement ratios, the single-building absolute projection displacement ratios and the satellite inclination tangent ratios. The disadvantage of this method is that it is difficult to accurately determine the geometric center point of the top contour of an irregularly shaped building with a complex top structure, which may cause projection displacement calculation errors and affect the accuracy of height change identification.

[0006] The invention patent application with the application number CN119131119A discloses a building height determination method and system, a storage medium and a computer program product, comprising: determining a top prediction box of a building to be identified and a plurality of bottom candidate boxes of the building to be identified in a remote sensing image; determining a bottom prediction box of the building to be identified among the plurality of bottom candidate boxes according to the geometric features of the top prediction box and the geometric features of the bottom candidate boxes; and analyzing the pixel offset between the top prediction box and the bottom prediction box to determine the height of the building to be identified, wherein the pixel offset represents the pixel dimension distance from the center point of the top prediction box to the center point of the bottom prediction box, and the pixel offset is obtained from the remote sensing image. The disadvantage of this method is that it is difficult to handle occluded or overlapping buildings by analyzing a single perspective remote sensing image; weather conditions affect the stability of top / bottom detection. SUMMARY

[0007] To solve the above technical problems, the present application proposes a method for calculating the height of buildings based on 3D Gaussian technology aerial data, which introduces 3DGS technology into city building three-dimensional modeling, and realizes interactive browsing and automatic height measurement task while taking into account the advantages of 3DGS technology for high-definition image data collected by unmanned aerial vehicles.

[0008] The object of the present application is to provide a method for calculating the height of buildings based on 3D Gaussian technology aerial data, comprising obtaining aerial photographs, and further comprising the following steps:

[0009] Step 1: data preprocessing;

[0010] Step 2: geographic coordinate conversion;

[0011] Step 3: Take all aerial photos, calculate their positions in the local coordinate system before conversion, and apply the conversion function f_best (x, y, z) to align the local coordinate system with the target coordinate system;

[0012] Step 4: Input the converted aerial photos' intrinsic parameters, extrinsic parameters, and sparse point cloud into the 3DGS algorithm for scene reconstruction;

[0013] Step 5: Obtain the user's click position;

[0014] Step 6: Calculate the building height at the user's click position.

[0015] Preferably, the step 1 includes downsampling the aerial photos, estimating the camera intrinsic parameters and extrinsic parameters.

[0016] In any of the above schemes, preferably, the step 2 includes the following sub-steps:

[0017] Step 21: Extract the latitude, longitude, and altitude from the EXIF information in the aerial photos and convert them to the target coordinate system;

[0018] Step 22: Randomly take several aerial photos, calculate their positions in the local coordinate system before conversion, then extract the corresponding converted coordinates, and calculate the similarity transformation function f(x, y, z);

[0019] Step 23: Take all aerial photos, calculate their positions in the local coordinate system before conversion, and input them into the similarity transformation function f(x, y, z) to obtain the converted coordinates X', Y', Z';

[0020] Step 24: For the same aerial photo, calculate the distance between the converted coordinates X', Y', Z' and X, Y, Z. If the distance is less than the threshold, it is recorded as an inlier, if the distance is greater than or equal to the threshold, it is recorded as an outlier, and the number of inliers is counted;

[0021] Step 25: When the number of inliers no longer increases, record the conversion function f_best(x, y, z), otherwise go to step 22.

[0022] In any of the above schemes, preferably, the conversion method includes obtaining the converted coordinates, the formula is:

[0023] X=(N+alt)*cos_lat*cos_lon

[0024] Y=(N+alt)*cos_lat*sin_lon

[0025] Z=(N*(1- e 2 )+alt)*sin_lat

[0026] wherein X, Y, Z are the transformed coordinates, alt is the altitude, cos_lat is the cosine of the latitude of the aerial photograph, cos_lon is the cosine of the longitude of the aerial photograph, sin_lat is the sine of the latitude of the aerial photograph, sin_lon is the sine of the longitude of the aerial photograph, e 2 is the square of the first eccentricity of the earth ellipsoid, N is the radius of curvature of the prime vertical, and the formula is:

[0027]

[0028] wherein a is the major axis of the earth, approximately 6,378,137 meters.

[0029] In any of the above solutions, preferably, the method of transformation comprises continuing the transformation to obtain the result coordinates, with the formula:

[0030] R * (point_i - point_0)

[0031] wherein point_i is the coordinates Xi, Yi, Zi of each point after transformation, point_0 is the coordinates X0, Y0, Z0 of the point corresponding to the first aerial photograph after transformation, and R is the rotation matrix, with the formula:

[0032]

[0033] wherein sin_lat_point_0 is the sine function applied to the latitude of the first aerial photograph, sin_lon_point_0 is the sine function applied to the longitude of the first aerial photograph, cos_lat_point_0 is the cosine function applied to the latitude of the first aerial photograph, and cos_lon_point_0 is the cosine function applied to the longitude of the first aerial photograph.

[0034] In any of the above solutions, preferably, the step 22 comprises estimating the transformation matrix between the two point sets by the umeyama algorithm, i.e. solving R and t to optimize the following function:

[0035]

[0036] The transformation matrix is solved as:

[0037]

[0038] wherein n is the total number of points, i is the i-th point, p i is the i-th point of the source point set, p i is the i-th point of the target point set, c is the scaling factor, R is the rotation matrix, t is the translation matrix, and T is the transformation matrix between the two point sets.

[0039] In any of the above solutions, preferably, the step 6 comprises the following sub-steps:

[0040] Step 61: Take the user click position as the tracking point, rotate the view angle around the tracking point for a circle, and synthesize the depth maps of multiple roof view angles;

[0041] Step 62: For each depth map, filter the depth using the confidence function and discard the wrong depth;

[0042] Step 63: Convert all pixels with correct depth into a point cloud;

[0043] Step 64: Take the maximum value of z coordinates in the point cloud, subtract the height of the building bottom, and get the height of the building.

[0044] In any of the above solutions, preferably, the step 61 comprises the following sub-steps:

[0045] Step 611: Obtain the synthesized view angle information;

[0046] Step 612: Construct a rotation matrix;

[0047] Step 613: Calculate the camera position, i.e., the translation vector;

[0048] Step 614: Construct a complete camera_to_world matrix, including the rotation matrix and the translation vector.

[0049] In any of the above solutions, preferably, the view angle information includes forward, right and down vectors, the forward vector is a unit vector pointing from the camera to the target point, the right vector is obtained by the cross product of the world_up vector and the forward vector, and the correct down vector is obtained by the cross product of the forward vector and the right vector.

[0050] In any of the above solutions, preferably, the right vector, the down vector and the forward vector are arranged in columns to obtain the rotation matrix.

[0051] In any of the above solutions, preferably, the formula of the translation vector is

[0052] T=[ camera_x, camera_y, camera_z]

[0053] camera_x = target_x + distance * sin(angle) * cos(pitch_angle)

[0054] camera_y = target_y + distance * sin(pitch_angle)

[0055] camera_z = target_z + distance * cos(angle) * cos(pitch_angle)

[0056] where camera_x, camera_y, camera_z are the x, y, z coordinates of the camera position, target_x, target_y, target_z are the x, y, z coordinates of the tracking point, distance is the distance set by the user from the camera position to the tracking point, angle is the horizontal angle, and pitch_angle is the pitch angle.

[0057] In any of the above solutions, preferably, the confidence calculation formula is:

[0058]

[0059] where N_confidence is the ordered sequence of all Gaussian points sorted by depth, α m_confidence is the opacity of the m_confidence-th Gaussian point, α n_confidence is the opacity of the n_confidence-th Gaussian point, m_confidence is the m_confidence-th Gaussian point in the ordered sequence, and n_confidence is the n_confidence-th Gaussian point in the ordered sequence.

[0060] In any of the above solutions, preferably, the confidence function formula is:

[0061]

[0062] where threshold = 0.999.

[0063] In any of the above solutions, preferably, when the confidence function calculation result of a certain pixel is 1, the depth synthesis is accurate, and the depth corresponding to the pixel is retained; when the confidence function calculation result of a certain pixel is 0, the depth synthesis is inaccurate, and the depth corresponding to the pixel is discarded.

[0064] In any of the above solutions, preferably, the step 63 includes the following sub-steps:

[0065] Step 631: For each pixel (u, v) in the depth map, whose corresponding depth value is d, calculate the three-dimensional coordinates (x c ,y c ,zc

[0066] Step 632: convert the point (x c ,y c ,z c ) in the camera coordinates into a point (x w ,y w ,z w ) in the world coordinates.

[0067] In any of the above solutions, preferably, the calculation formula of the three-dimensional coordinates (x c ,y c ,z c ) is

[0068]

[0069] wherein K is the camera intrinsic parameter, and d is the depth value corresponding to the pixel (u, v).

[0070] In any of the above solutions, preferably, the conversion formula of the point (x w ,y w ,z w ) is

[0071]

[0072] wherein c2w is the transformation matrix from the camera coordinate system to the world coordinate system.

[0073] The present application provides a kind of based on 3D Gaussian technique's aerial data building height calculation method, multi-view height finding (highest point) process, traverse each view angle near roof, without artificial special finding, strong universality, realize automatic finding highest point. BRIEF DESCRIPTION OF DRAWINGS

[0074] Figure 1 It is the flow chart of a preferred embodiment of the building height calculation method based on 3D Gaussian technique's aerial data according to the present application.

[0075] Figure 2 It is the flow chart of another preferred embodiment of the building height calculation method based on 3D Gaussian technique's aerial data according to the present application.

[0076] Figure 3 It is the schematic diagram of an embodiment of user selection building of the building height calculation method based on 3D Gaussian technique's aerial data according to the present application.

[0077] Figure 4 It is the schematic diagram of an embodiment of view angle rotation of the building height calculation method based on 3D Gaussian technique's aerial data according to the present application. DETAILED DESCRIPTION​

[0078] The application will be further described below in connection with the accompanying drawings and specific embodiments.

[0079] Embodiment one

[0080] As shown in the figure, step 1000 is performed to obtain aerial photos. Figure 1

[0081] Step 1100 is performed to perform data preprocessing, including downsampling the aerial photos, estimating camera intrinsic parameters and extrinsic parameters.

[0082] Step 1200 is performed to perform geographic coordinate conversion, including the following sub-steps:

[0083] Step 1210 is performed to extract the latitude, longitude and altitude in the EXIF information in the aerial photos, and convert them to the target coordinate system, the conversion method including obtaining the converted coordinates, the formula being:

[0084] X=(N+alt)*cos_lat*cos_lon

[0085] Y=(N+alt)*cos_lat*sin_lon

[0086] Z=(N*(1- e 2 )+alt)*sin_lat

[0087] Wherein, X, Y, Z are the converted coordinates, alt is the altitude, cos_lat is the cosine value of the latitude of the aerial photo, cos_lon is the cosine value of the longitude of the aerial photo, sin_lat is the sine value of the latitude of the aerial photo, sin_lon is the sine value of the longitude of the aerial photo, e 2 is the square of the first eccentricity of the earth ellipsoid, N is the radius of curvature of the prime vertical circle, the formula being:

[0088]

[0089] Wherein, α is the long semi-axis of the earth, about 6,378,137 meters.

[0090] The conversion method includes continuing to convert to obtain the result coordinates, the formula being:

[0091] R*(point_i-point_0)

[0092] Wherein, point_i is the coordinates Xi, Yi, Zi of each point after conversion, point_0 is the point coordinates X0, Y0, Z0 corresponding to the first aerial photo after conversion, R is the rotation matrix, the formula being:

[0093] ​

[0094] wherein sin_lat_point_0 is a sine function applied to the latitude of the first aerial photograph, sin_lon_point_0 is a sine function applied to the longitude of the first aerial photograph, cos_lat_point_0 is a cosine function applied to the latitude of the first aerial photograph, and cos_lon_point_0 is a cosine function applied to the longitude of the first aerial photograph.

[0095] Step 1220 is performed, a few aerial photographs are randomly taken, their positions in the local coordinate system before conversion are calculated, and then the corresponding converted coordinates are extracted, and a similarity transformation function f(x, y, z) is calculated.

[0096] The transformation matrix between the two point sets is estimated by the umeyama algorithm, that is, R and t are solved to make the following function optimal:

[0097]

[0098] The transformation matrix is solved:

[0099]

[0100] wherein n is the total number of points, i is the i-th point, p i is the i-th point of the source point set, p i is the i-th point of the target point set, c is a scaling factor, R is a rotation matrix, t is a translation matrix, and T is the transformation matrix between the two point sets.

[0101] Step 1230 is performed, all aerial photographs are taken, their positions in the local coordinate system before conversion are calculated, and they are input into the similarity transformation function f(x, y, z) to obtain the converted coordinates X', Y', and Z'.

[0102] Step 1240 is performed, for the same aerial photograph, the distance between the converted coordinates X', Y', and Z' and X, Y, and Z is calculated, if the distance is less than a threshold value, it is recorded as an inlier, if the distance is greater than or equal to the threshold value, it is recorded as an outlier, and the number of inliers is counted.

[0103] Step 1250 is performed, when the number of inliers no longer increases, the conversion function f_best(x, y, z) is recorded, otherwise step 1220 is performed.

[0104] Step 1300 is performed, all aerial photographs are taken, their positions in the local coordinate system before conversion are calculated, and the conversion function f_best(x, y, z) is applied to align the local coordinate system with the target coordinate system.

[0105] Perform step 1400, input the converted inner parameter, outer parameter and sparse point cloud of aerial photograph into 3DGS algorithm to perform scene reconstruction.

[0106] Perform step 1500 to obtain the user click position.

[0107] Perform step 1600 to calculate the building height of the user click position, including the following sub-steps:

[0108] Perform step 1610 to take the user click position as a tracking point, rotate the view angle around the tracking point for a circle, and synthesize a plurality of depth maps of roof view angles, including the following sub-steps:

[0109] Perform step 1611 to obtain the synthesized view angle information, the view angle information including forward vector, right vector and down vector, the forward vector being a unit vector from the camera to the target point, the right vector being obtained by the cross product of the world_up vector and the forward vector, and the correct down vector being obtained by the cross product of the forward vector and the right vector.

[0110] Perform step 1612 to construct a rotation matrix, and arrange the right vector, the down vector and the forward vector in columns to obtain the rotation matrix.

[0111] Perform step 1613 to calculate the camera position, i.e. the translation vector, the formula of the translation vector being

[0112] T = [ camera_x, camera_y, camera_z]

[0113] camera_x = target_x + distance * sin(angle) * cos(pitch_angle)

[0114] camera_y = target_y + distance * sin(pitch_angle)

[0115] camera_z = target_z + distance * cos(angle) * cos(pitch_angle)

[0116] Wherein, camera_x, camera_y, camera_z are the x, y, z coordinates of the camera position, target_x, target_y, target_z are the x, y, z coordinates of the tracking point, distance is the distance set by the user from the camera position to the tracking point, angle is the horizontal angle, and its value is equal to 360 / user specified interval number; pitch_angle is the pitch angle, which is 45 degrees.

[0117] Step 1614 is performed to build a complete camera_to_world matrix, including a rotation matrix and a translation vector.

[0118] Step 1620 is performed to filter the depth using a confidence function for each depth map, and discard the wrong depth, and the confidence calculation formula is:

[0119]

[0120] Wherein, N_confidence is the ordered sequence of all Gaussian points sorted by depth, α m_confidence is the opacity of the m_confidence Gaussian point, α n_confidence is the opacity of the n_confidence Gaussian point, m_confidence is the m_confidence Gaussian point in the ordered sequence, and n_confidence is the n_confidence Gaussian point in the ordered sequence.

[0121] The confidence function formula is:

[0122]

[0123] Wherein, threshold=0.999.

[0124] When the confidence function calculation result of a certain pixel is 1, the depth synthesis is accurate, and the depth corresponding to the pixel is retained; when the confidence function calculation result of a certain pixel is 0, the depth synthesis is inaccurate, and the depth corresponding to the pixel is discarded.

[0125] Step 1630 is performed to convert all pixels with correct depth into point cloud, including the following sub-steps:

[0126] Step 1631 is performed to calculate the three-dimensional coordinates (x c ,y c ,z c ) of the point in the camera coordinate system for each pixel (u, v) in the depth map, wherein the depth value corresponding to the pixel is d. c ,y c ,zc The calculation formula of the camera coordinate system is

[0127]

[0128] wherein K is the camera internal parameter, and d is the depth value corresponding to the pixel (u, v).

[0129] In step 1632, the point (x c ,y c ,z c ) in the camera coordinate system is converted into the point (x w ,y w ,z w ) in the world coordinate system, and the conversion formula of the point (x w ,y w ,z w ) is

[0130]

[0131] wherein c2w is the transformation matrix from the camera coordinate system to the world coordinate system.

[0132] In step 1640, the maximum value of the z coordinate in the point cloud is taken, and the building bottom height is subtracted to obtain the building height.

[0133] Embodiment two

[0134] The application discloses a kind of based on 3D Gauss technique's aerial data building height calculation method, scene reconstruction is carried out based on 3DGS algorithm, only need to provide unmanned aerial vehicle aerial picture;Then, develop multi-view high (highest point) process, traverse a circle around roof, find the highest point by the information of multiple perspectives, realize high precision, low cost and small difficulty.3DGS model directly restores real spatial distance.Using coordinate system conversion function, the aerial picture shooting position (local coordinate system) calculated by 3DGS algorithm is converted to global unified coordinate system (unit meter), and using outlier rejection algorithm, the robustness of conversion function is improved, and real distance is restored.In depth map synthesis, carry out outlier rejection: based on the confidence of the depth value of each pixel based on the cumulative opacity of pixel, design confidence filtering function to filter depth, and remove depth value with larger error.

[0135] As shown in Figure 2 , a kind of based on 3D Gauss technique's aerial data building height calculation method includes the following steps:

[0136] Step 1: data preprocessing

[0137] Data preprocessing includes downsampling (according to actual hardware equipment and rendering demand) to aerial photograph, estimating camera internal parameter and external parameter and the like.

[0138] Step 2: Coordinate system conversion

[0139] Step 2.1: Extract the latitude, longitude, and altitude from the EXIF information in the aerial photograph, and convert them directly to the target coordinate system using the following formulas:

[0140] X = (N + alt) * cos_lat * cos_lon

[0141] Y = (N + alt) * cos_lat * sin_lon

[0142] Z = (N * (1 - e 2 ) + alt) * sin_lat

[0143] where X, Y, and Z are the converted coordinates; suffix lat is the latitude, suffix lon is the longitude, prefix sin is the sine function, prefix cos is the cosine function, alt is the altitude in meters; e 2 is the first eccentricity of the Earth's ellipsoid, approximately 0.00669438002290, and N is the radius of the prime vertical circle, with the formula:

[0144]

[0145] where α is the length of the Earth's semi-major axis, approximately 6,378,137 meters.

[0146] Then, continue the conversion to get the result coordinates:

[0147] R * (point_i - point_0)

[0148] where point_i is the coordinate Xi, Yi, Zi of each point after the above conversion, point_0 is the coordinate X0, Y0, Z0 of the first aerial photograph after the above conversion, and R is the rotation matrix, with the formula:

[0149] [-sin_lon, cos_lon, 0., -sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat, cos_lat * cos_lon, cos_lat * sin_lon, sin_lat;]

[0150] where suffix lat is the latitude of the first aerial photograph, suffix lon is the longitude of the first aerial photograph, prefix sin is the sine function, and prefix cos is the cosine function.

[0151] Step 2.2: Randomly take several aerial photos, calculate their positions in the local coordinate system before transformation, then extract the corresponding coordinates after transformation (result of step 2.1), and calculate the similarity transformation function f(x, y, z).

[0152] Step 2.3: Take all aerial photos, calculate their positions in the local coordinate system before transformation, and input them into the similarity transformation function f(x, y, z) to obtain the transformed coordinates X', Y', Z'.

[0153] Step 2.4: For the same aerial photo, calculate the distance between the transformed coordinates X', Y', Z' and X, Y, Z. If the distance is less than the threshold, it is considered an "inlier", if the distance is greater than or equal to the threshold, it is considered an "outlier", and count the number of "inliers".

[0154] Step 2.5: If the number of "inliers" no longer increases, record the transformation function f_best(x, y, z) and go to step 3, otherwise go to step 2.2.

[0155] Step 3: Take all aerial photos, calculate their positions in the local coordinate system before transformation, and apply the transformation function f_best(x, y, z) to align the local coordinate system with the target coordinate system.

[0156] Step 4: Input the intrinsic parameters, extrinsic parameters, and sparse point cloud of the transformed aerial photos into the 3DGS algorithm to perform scene reconstruction.

[0157] Step 5: The system reads the user's click position, as shown in Figure 3 .

[0158] Step 6: Height measurement process

[0159] Step 6.1: As shown in Figure 4 , use the user's click position as the tracking point, and rotate the view around the tracking point for a full circle to synthesize multiple depth maps of the roof view.

[0160] Where, the synthesized view information needs the following parameters:

[0161] forward vector: a unit vector pointing from the camera to the target point

[0162] right vector: obtained by the cross product of the world_up vector and the forward vector

[0163] down vector: obtained by the cross product of the forward vector and the right vector

[0164] Construct a rotation matrix: arrange the right, down, and forward vectors in columns

[0165] Compute camera position, i.e. translation vector

[0166] Construct full camera_to_world matrix: contains rotation matrix and translation vector

[0167] Step 6.2: For each depth image, filter the depth using a confidence function, discard the wrong depth.

[0168] The confidence formula is:

[0169]

[0170] where N_confidence is the ordered sequence of all the gauss points by depth, a m_confidence is the opacity of the m_confidence-th gauss point, a n_confidence is the opacity of the n_confidence-th gauss point, m_confidence is the m_confidence-th gauss point in the ordered sequence, n_confidence is the n_confidence-th gauss point in the ordered sequence, when the number of gauss points is enough, mathematically, the upper limit of the whole sum is 1. Therefore, if the cumulative opacity reaches 0.999, the confidence of the depth value reaches 0.999.

[0171] The confidence function formula is:

[0172]

[0173] threshold is generally taken as 0.999, when the confidence of a certain pixel is input into the formula, if the result is 1, the depth synthesis is accurate, and the depth corresponding to the pixel is retained; if the result is 0, the depth synthesis is inaccurate, and the depth corresponding to the pixel is discarded.

[0174] Step 6.3: Convert all pixels with correct depth into point cloud.

[0175] For each pixel (u, v) in the depth image, its corresponding depth value is d, then the three-dimensional coordinates (x c ,y c ,z c ) of the point in the camera coordinate system can be calculated by the following formula:

[0176]

[0177] where K is the camera intrinsic parameter, (x c ,y c ,z c ) is the three-dimensional coordinates in the camera coordinate system.

[0178] Then, the point (x c ,y c ,z c ) in the camera coordinate is converted into the point (x w ,y w ,z w ) in the world coordinate.

[0179]

[0180] wherein c2w is the transformation matrix from the camera coordinate system to the world coordinate system.

[0181] Step 6.4: Take the maximum value of the z coordinate in the point cloud, subtract the height of the building bottom, and obtain the height of the building.

[0182] Step 6.5: The user selects whether to continue, if yes, go to step 6.1; otherwise, end.

[0183] In order to better understand the present application, the above is described in detail in combination with the specific embodiments of the present application, but is not a limitation on the present application. Any simple modification made to the above embodiments according to the technical essence of the present application still belongs to the scope of the technical solution of the present application. In the specification, each embodiment focuses on the difference from other embodiments, and the same or similar parts of each embodiment can be referred to. For system embodiments, since they basically correspond to method embodiments, the description is relatively simple, and the relevant parts can be referred to the part of the method embodiment.

Claims

1. A method for calculating building height based on aerial photography data using 3D Gaussian technology, comprising acquiring aerial photographs, characterized in that, It also includes the following steps: Step 1: Perform data preprocessing; Step 2: Perform geographic coordinate transformation; Step 3: Take all aerial photographs, calculate their positions in the local coordinate system before transformation, and apply the transformation function f_best(x,y,z) to align the local coordinate system with the target coordinate system; Step 4: Input the intrinsic parameters, extrinsic parameters, and sparse point cloud data of the converted aerial photographs into the 3DGS algorithm for scene reconstruction; Step 5: Obtain the user's click location; Step 6: Calculate the building height at the user-clicked location, including the following sub-steps: Step 61: Using the user's click location as the tracking point, rotate the viewpoint around the tracking point to synthesize multiple depth maps from the rooftop viewpoint; Step 62: For each depth map, use a confidence function to filter the depths, discarding erroneous depths. The confidence calculation formula is: ， Where N_confidence is the ordered sequence of all Gaussian points sorted by depth, and α m_confidence Let α be the opacity at the m-th Gaussian point. n_confidence Let be the opacity of the nth Gaussian point, m_confidence be the mth Gaussian point in the ordered sequence, and n_confidence be the nth Gaussian point in the ordered sequence. The confidence function formula is: ， Where threshold = 0.999; Step 63: Convert all pixels with the correct depth into point clouds; Step 64: Take the maximum z-coordinate in the point cloud, subtract the building's bottom height, and obtain the building height.

2. The method for calculating building height based on aerial photography data using 3D Gaussian technology as described in claim 1, characterized in that, Step 1 includes downsampling the aerial photographs and estimating the camera's intrinsic and extrinsic parameters.

3. The method for calculating building height based on aerial photography data using 3D Gaussian technology as described in claim 2, characterized in that, Step 2 includes the following sub-steps: Step 21: Extract the latitude, longitude, and altitude from the EXIF ​​information of the aerial photograph and convert them to the target coordinate system; Step 22: Randomly select several aerial photographs, calculate their positions in the local coordinate system before transformation, then extract the corresponding transformed coordinates, and calculate the similarity transformation function f(x,y,z); Step 23: Take all aerial photographs, calculate their positions in the local coordinate system before the transformation, and input them into the similarity transformation function f(x,y,z) to obtain the transformed coordinates X', Y', Z'; Step 24: For the same aerial photograph, calculate the distance between the converted coordinates X', Y', Z' and X, Y, Z. If the distance is less than the threshold, it is recorded as an interior point; if the distance is greater than or equal to the threshold, it is recorded as an exterior point. Count the number of interior points. Step 25: When the number of interior points no longer increases, record the transformation function f_best(x,y,z); otherwise, go to step 22.

4. The method for calculating building height based on aerial photography data using 3D Gaussian technology as described in claim 3, characterized in that, The transformation method includes obtaining the transformed coordinates, using the following formula: X=(N+alt)*cos_lat*cos_lon Y=(N+alt)*cos_lat*sin_lon Z=(N*(1-e 2 )+alt)*sin_lat Where X, Y, and Z are the transformed coordinates, alt is the altitude, cos_lat is the cosine of the latitude of the aerial photograph, cos_lon is the cosine of the longitude of the aerial photograph, sin_lat is the sine of the latitude of the aerial photograph, sin_lon is the sine of the longitude of the aerial photograph, and e 2 Let N be the square of the first eccentricity of the Earth's ellipsoid, and N be the radius of curvature of the prime mover. The formula is: ， Where α is the Earth's semi-major axis.

5. The method for calculating building height based on aerial photography data using 3D Gaussian technology as described in claim 4, characterized in that, The transformation method includes further transformation to obtain the resulting coordinates, as shown in the formula: R*(point_i-point_0) Where point_i are the coordinates Xi, Yi, Zi of each point after transformation, point_0 are the coordinates X0, Y0, Z0 of the point corresponding to the first aerial photograph after transformation, and R is the rotation matrix, with the formula: ， Here, sin_lat_point_0 applies a sine function to the latitude of the first aerial photograph, sin_lon_point_0 applies a sine function to the longitude of the first aerial photograph, cos_lat_point_0 applies a cosine function to the latitude of the first aerial photograph, and cos_lon_point_0 applies a cosine function to the longitude of the first aerial photograph.

6. The method for calculating building height based on aerial photography data using 3D Gaussian technology as described in claim 5, characterized in that, Step 61 includes the following sub-steps: Step 611: Obtain the synthesized viewpoint information; Step 612: Construct the rotation matrix; Step 613: Calculate the camera position, i.e., the translation vector; Step 614: Construct the complete camera_to_world matrix, including the rotation matrix and translation vector.

7. The method for calculating building height based on aerial photography data using 3D Gaussian technology as described in claim 6, characterized in that, The viewpoint information includes a forward vector, a right vector, and a down vector. The forward vector is a unit vector pointing from the camera to the target point. The right vector is obtained by the cross product of the world_up vector and the forward vector. The correct down vector is obtained by the cross product of the forward vector and the right vector.

8. The method for calculating building height based on aerial photography data using 3D Gaussian technology as described in claim 7, characterized in that, The rotation matrix is ​​obtained by arranging the right vector, down vector, and forward vector column by column.

Images