3DGS-based 3D reconstruction method for panoramic images

Abstract

This invention provides a method for 3D reconstruction of panoramic images based on 3DGS, which relates to the fields of computer vision and image processing. The method comprises: Step 1: Acquiring a sequence of panoramic images of the scene to be reconstructed and uniformly encoding them via equirectangular projection to obtain a coded dataset; Step 2: Performing panoramic structure-from-motion processing on the coded dataset to generate six-degree-of-freedom camera pose parameters for each image and a corresponding sparse 3D point cloud. This method achieves high-precision and high-efficiency 3D reconstruction of panoramic scenes, capable of outputting high-quality panoramic images or dense 3D models, improving both the quality and efficiency of 3D reconstruction.

Description

Technical Field

[0001] The present invention relates to the technical field of computer vision and image processing, and in particular to a method for three-dimensional reconstruction of panoramic images based on 3DGS. Background Art

[0002] With the development of VR, AR, and digital twin technologies, 3D reconstruction based on panoramic images has broad prospects in scenarios such as urban modeling. However, existing technologies still have some areas for improvement:

[0003] The feature point extraction and elimination process relies on traditional algorithms. The feature matching accuracy may be affected in partially occluded areas, with repeated or weak textures. There is also room for improvement in stability in scenes with large viewing angles.

[0004] When dealing with extreme viewing angles, irregularly structured scenes, or areas with large parallax, there may be room for improvement in point cloud density and structural clarity. Current technology may suffer from reduced reconstruction accuracy due to image distortion in polar regions. Resource efficiency and convergence speed could be improved in repetitive scene reconstruction. Traditional projection models are still not well-suited to the spatial consistency requirements of panoramic viewpoints. Gaussian modeling also faces challenges, such as the need for enhanced detail control and further optimization of redundant points. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a panoramic image three-dimensional reconstruction method based on 3DGS, which effectively improves the accuracy and efficiency of panoramic image three-dimensional reconstruction and reduces training resource consumption.

[0006] In order to solve the above technical problems, the technical solutions of the present invention are as follows:

[0007] In a first aspect, a method for 3D reconstruction of panoramic images based on 3DGS, the method comprising:

[0008] Step 1: Collect a panoramic image sequence of the scene to be reconstructed, perform unified encoding through equirectangular projection, and obtain a coded data set;

[0009] Step 2: Perform panoramic motion recovery structure processing on the encoded dataset to generate the six-degree-of-freedom camera pose parameters and corresponding sparse three-dimensional point cloud for each image;

[0010] Step 3: Initialize a 3D Gaussian primitive set based on the sparse 3D point cloud, assigning position, rotation parameters, initial scale value, color parameters, and opacity parameters to each 3D point. If there is a historical reconstruction result that coincides with the current scene, fuse the historical reconstruction result with the current Gaussian primitive set through pose registration to obtain a fused Gaussian primitive set.

[0011] Step 4: Perform distortion correction on the top or bottom polar regions of the panoramic image sequence, and reconstruct the corrected image set using local spherical remapping or cubic projection. Iterative training is performed based on the fused Gaussian basis set and the corrected image set. During the training process, the Gaussian basis set is densified based on the spatial gradient of the panoramic screen. Training is terminated when the iteration limit is reached, resulting in an optimized Gaussian basis set.

[0012] Step 5: Perform panoramic rendering on the optimized Gaussian primitive set to output a panoramic image or a dense three-dimensional model.

[0013] Furthermore, a panoramic image sequence of the scene to be reconstructed is collected and uniformly encoded through equirectangular projection to obtain a coded data set, including:

[0014] Capturing a multi-view spherical image sequence of a scene by a panoramic camera;

[0015] Each spherical image is converted into a two-dimensional plane image using equidistant cylindrical projection, and a mapping relationship between spherical view coordinates and plane pixel coordinates is established;

[0016] All converted images are resolution-normalized to generate a coded dataset in a unified format.

[0017] Furthermore, panoramic structure-from-motion processing is performed on the encoded dataset to generate the six-degree-of-freedom camera pose parameters and corresponding sparse three-dimensional point clouds for each image, including:

[0018] For the images in the unified format encoding dataset, feature points with panoramic invariance are extracted, and matching point pairs are filtered according to the spherical perspective similarity to obtain the filtered matching point pairs;

[0019] Based on the screened matching point pairs, the camera pose is solved using the incremental structure-from-motion algorithm to generate the initial 3D point coordinates.

[0020] The camera pose and 3D point coordinates are jointly optimized using bundle adjustment, the panoramic projection consistency of each 3D point is verified, and the optimized camera pose parameter set and a sparse 3D point cloud that meets the consistency constraints are output.

[0021] Furthermore, based on the filtered matching point pairs, the camera pose is solved using an incremental structure-from-motion algorithm to generate the initial 3D point coordinates, including:

[0022] Taking the camera position of the first image as the origin of the world coordinate system, initialize the first camera pose based on the filtered matching point pairs, and mark the corresponding view as a registered view;

[0023] From the remaining unregistered views, a new view that has common view feature points with the registered view is selected, and the rotation matrix and translation vector of the new view relative to the world coordinate system are solved using the PnP algorithm to obtain the six-degree-of-freedom camera pose;

[0024] Based on the new view pose and common view feature points, a triangulation operation is performed on the new view and the registered view to generate the initial 3D point coordinates bound to the world coordinate system.

[0025] Furthermore, a 3D Gaussian primitive set is initialized based on the sparse 3D point cloud, and each 3D point is assigned a position, rotation parameter, initial scale value, color parameter, and opacity parameter. When there is a historical reconstruction result that coincides with the current scene, the historical reconstruction result is fused with the current Gaussian primitive set through pose registration to obtain a fused Gaussian primitive set, including:

[0026] Convert the sparse 3D point cloud into an initial Gaussian primitive set, where each Gaussian primitive inherits the position and color parameters of the 3D point, and initializes the rotation parameters to no rotation, the scale to a microcube, and the opacity to 0.8;

[0027] When a historical reconstruction model is detected, feature matching is performed based on the initial Gaussian primitive set, and the reprojection error of the historical Gaussian primitive is calculated using the camera pose parameters, and invalid Gaussians with an error of more than 2 pixels are eliminated;

[0028] The verified historical Gaussian primitives are merged with the initial Gaussian primitive set, and the final fused Gaussian primitive set is output.

[0029] Furthermore, a distortion correction operation is performed on the top or bottom polar region of the panoramic image sequence, and a corrected image set is obtained by local spherical remapping or cubic projection reconstruction. Iterative training is performed based on the fused Gaussian basis set and the corrected image set. During the training process, the Gaussian basis set is densified based on the spatial gradient of the panoramic screen. When the iteration limit is reached, the training is terminated to obtain an optimized Gaussian basis set, including:

[0030] For the top or bottom polar region of the panoramic image sequence, local spherical remapping or cubic projection reconstruction is used to output a corrected image set;

[0031] The image set is corrected based on the fused Gaussian basis set and iterative training is performed. When the number of iterations reaches the upper limit, the training is terminated and the optimized Gaussian basis set is output.

[0032] Furthermore, the image set is corrected based on the fused Gaussian basis set and iterative training is performed. When the number of iterations reaches the upper limit, the training is terminated and the optimized Gaussian basis set is output, including:

[0033] Use the rectified image set as training input data to initialize the parameters of the fused Gaussian basis set;

[0034] Randomly select training views and corresponding camera poses from the rectified image set, and generate simulated views through the equirectangular projection rendering pipeline based on the camera pose and the current Gaussian primitive set;

[0035] Calculate the loss function between the simulated view and the real training view, and update the position, rotation, scale, color, and opacity parameters of the Gaussian primitives based on the backpropagation of the loss function to obtain the updated Gaussian primitive set;

[0036] Based on the updated Gaussian basis set, the gradient magnitude of each Gaussian basis in the panoramic screen space is calculated. If the gradient magnitude is greater than the splitting threshold, a splitting operation is performed to generate a sub-Gaussian; if the gradient magnitude is less than the cloning threshold, a cloning operation is performed; if the gradient magnitude is equal to the cloning threshold, the original state is maintained;

[0037] When the number of iterations reaches the preset upper limit, the training is terminated and the optimized Gaussian basis set is output.

[0038] Furthermore, the optimized Gaussian primitive set is subjected to panoramic rendering to output a panoramic image or a dense 3D model, including:

[0039] The optimized Gaussian primitive set is input into the equirectangular projection rendering pipeline, and an equirectangular projection grid is generated according to the output resolution of the panoramic image.

[0040] For each pixel in the grid, a ray is generated from the camera center according to the spherical view direction corresponding to the plane coordinate mapping;

[0041] Calculate the intersection weights of the ray with all Gaussian primitives in the scene, and generate the final color of the current pixel by weighted accumulation of the color values ​​and opacity contributed by each Gaussian.

[0042] Traverse all pixels to complete panoramic image rendering output, or generate a dense 3D point cloud model based on the spatial position and scale parameters of the optimized Gaussian primitive set.

[0043] In a second aspect, a computing device includes:

[0044] one or more processors;

[0045] The storage device is used to store one or more programs, and when the one or more programs are executed by the one or more processors, the one or more processors implement the method.

[0046] According to a third aspect, a computer-readable storage medium stores a program, which implements the method described above when executed by a processor.

[0047] The above solution of the present invention includes at least the following beneficial effects:

[0048] Panoramic images are encoded using equirectangular projection, and panoramic SfM is used to generate camera poses and sparse point clouds, ensuring data format consistency and camera parameter accuracy. When historical reconstruction results exist, Gaussian primitives are fused through pose registration to avoid repeated training and reduce computing resource consumption. For example, in the same or overlapping scenes, historical high-confidence Gaussians can be directly used as initial parameters to accelerate optimization convergence and improve reconstruction efficiency. To address the severe distortion problem of equirectangular projection in the polar regions (±90° viewing angle), local spherical remapping or cube map stitching is used to correct the spatial distortion of polar images, avoid misleading reconstruction caused by projection errors, and improve modeling stability and accuracy in edge areas (such as the top / bottom of the image).

[0049] Dynamically adjusting the Gaussian distribution based on the spatial gradient of the panoramic screen enables adaptive modeling of spatially non-uniform structures. In resource-constrained scenarios, it prioritizes optimization of key areas, balancing reconstruction efficiency and detail. The rendering framework based on equirectangular projection supports 360° omnidirectional viewing angles. Combining a joint loss function of photometric error and structural similarity error (MSE+SSIM) ensures lighting consistency and visual continuity in the rendered image, ultimately outputting a high-precision panoramic image or dense 3D model. Through mechanisms such as distortion correction and gradient-driven densification, it overcomes the matching error issues of traditional methods in scenes such as polar regions and weak textures. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 It is a flowchart of a method for 3D reconstruction of panoramic images based on 3DGS provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0051] Exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although exemplary embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. Rather, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0052] like Figure 1 As shown, an embodiment of the present invention proposes a method for 3D reconstruction of panoramic images based on 3DGS, the method comprising the following steps:

[0053] Step 1: Collect a panoramic image sequence of the scene to be reconstructed, perform unified encoding through equirectangular projection, and obtain a coded data set;

[0054] Step 2: Perform panoramic motion recovery structure processing on the encoded dataset to generate the six-degree-of-freedom camera pose parameters and corresponding sparse three-dimensional point cloud for each image;

[0055] Step 3: Initialize a 3D Gaussian primitive set based on the sparse 3D point cloud, assigning position, rotation parameters, initial scale value, color parameters, and opacity parameters to each 3D point. If there is a historical reconstruction result that coincides with the current scene, fuse the historical reconstruction result with the current Gaussian primitive set through pose registration to obtain a fused Gaussian primitive set.

[0056] Step 4: Perform distortion correction on the top or bottom polar regions of the panoramic image sequence, and reconstruct the corrected image set using local spherical remapping or cubic projection. Iterative training is performed based on the fused Gaussian basis set and the corrected image set. During the training process, the Gaussian basis set is densified based on the spatial gradient of the panoramic screen. Training is terminated when the iteration limit is reached, resulting in an optimized Gaussian basis set.

[0057] Step 5: Perform panoramic rendering on the optimized Gaussian primitive set to output a panoramic image or a dense three-dimensional model.

[0058] In an embodiment of the present invention, a panoramic image sequence is uniformly encoded through equidistant cylindrical projection to ensure data format consistency; six-degree-of-freedom camera poses and sparse point clouds are generated to accurately locate the camera view and scene structure, laying the spatial coordinate foundation for Gaussian model initialization. Gaussian primitives are initialized based on sparse point clouds, and pose registration and fusion of historical reconstruction results are combined to avoid repeated training of repeated scenes, accelerate model convergence, and reduce computing resource consumption. Polar projection distortion is corrected through spherical remapping or cube mapping to improve the modeling accuracy of edge areas; Gaussian distribution is dynamically optimized based on the densification control (splitting, cloning, pruning) of panoramic gradients to achieve encryption of detail areas and elimination of redundant points, balancing reconstruction efficiency and accuracy. Based on the optimized Gaussian set, omnidirectional continuous panoramic images or dense three-dimensional models are generated to meet the high-precision visual requirements of scenes such as VR / AR, ensuring light consistency and structural integrity.

[0059] In a preferred embodiment of the present invention, the above step 1 of acquiring a panoramic image sequence of the scene to be reconstructed and uniformly encoding it through equirectangular projection to obtain a coded data set may include:

[0060] Step 100, capturing a multi-view spherical image sequence of a scene by a panoramic camera;

[0061] Step 101: convert each spherical image into a two-dimensional plane image using equirectangular projection, and establish a mapping relationship between spherical view coordinates and plane pixel coordinates;

[0062] Step 102: Perform resolution standardization on all converted images to generate a coded data set in a unified format.

[0063] In an embodiment of the present invention, a panoramic camera (such as a fisheye lens stitching camera or a multi-lens array camera) is used, and the lens field of view must cover a horizontal 360° and a vertical 180° range to ensure that there are no blind spots. Acquisition parameter settings: the exposure time is dynamically adjusted according to the scene lighting (such as ISO100-800), and the resolution is usually set to 4096×2048 or higher to retain sufficient details; the frame rate must meet the static or dynamic shooting requirements of the scene (such as 1-5fps for static scenes and ≥30fps for dynamic scenes). If fixed-position acquisition is used, the camera remains in the center position and takes a set of images at a preset angle (such as every 10°) by rotating the pan-tilt head to ensure that the overlap rate of adjacent images is ≥50%; if mobile acquisition is used, the camera trajectory must be recorded (such as through real-time positioning through SLAM) to ensure position consistency between consecutive frames.

[0064] Image sequence acquisition process:

[0065] Shoot sequentially in a clockwise or counterclockwise direction, covering the entire viewing angle to form a closed image sequence. For example, horizontally, starting at 0°, capture one image every 10°, for a total of 36. Vertically, capture can be divided into multiple layers (such as zenith, horizontal, and nadir) to ensure that all areas of the sphere are covered. During the capture process, image clarity, exposure consistency, and view overlap are checked in real time to eliminate blurry, overexposed, or insufficiently overlapped images.

[0066] Step 101: Each pixel of the spherical image corresponds to a direction vector in three-dimensional space, expressed as longitude θ (range 0~2π) and latitude φ (range - ~ ) is represented by the equidistant cylindrical projection. The spherical surface is "unfolded" into a rectangular plane. The horizontal coordinate x of the plane image corresponds to θ, and the vertical coordinate y corresponds to φ. The mapping relationship is: the horizontal coordinate x of the plane = ( )× plane image width W, where θ=0 corresponds to the left edge of the plane, and θ=2π corresponds to the right edge (needs to be spliced ​​cyclically); the plane vertical coordinate y=( +0.5)×plane image height H, where φ=- Corresponding to the bottom of the plane (nadir), φ= Corresponds to the top (zenith).

[0067] Pixel value interpolation calculation process:

[0068] For the target pixel coordinates (x, y) in the plane image, according to the mapping relationship of the equidistant cylindrical projection, its direction coordinates in the original spherical image are reversed:

[0069] Longitude θ = (x ÷ W) × 2π, where W is the width of the plane image and θ is in the range [0, 2π);

[0070] Latitude φ=(y÷H-0.5)×π, where H is the plane image height, and the range of φ is [- , ]. Since x and y are integer pixel coordinates, the calculated θ and φ are non-integer (such as θ = 3.1416 radians), while the pixel coordinates of the original spherical image are integers (such as θ1 = 3, θ2 = 4), so it is necessary to determine the pixel values ​​at non-integer coordinates through interpolation.

[0071] Separate the integer and fractional parts of θ and φ:

[0072] Let the integer part of θ be θ_int = floor(θ), and the fractional part be θ_frac = θ - θ_int (in the range [0, 1));

[0073] Let the integer part of φ be φ_int = floor(φ), and the fractional part be φ_frac = φ - φ_int (in the range [0, 1)).

[0074] Thus, the four pixel points adjacent to (θ, φ) in the original spherical image are determined:

[0075] Upper left point: (θ1, φ1) = (θ_int, φ_int);

[0076] Upper right point: (θ2, φ1) = (θ_int+1, φ_int);

[0077] Lower left point: (θ1, φ2) = (θ_int, φ_int+1);

[0078] Lower right point: (θ2, φ2)=(θ_int+1, φ_int+1);

[0079] (If θ_int = 2π-1, then θ2 needs to loop back to 0 to ensure the continuity of the left and right boundaries of the spherical image).

[0080] The weight coefficients are based on the fractional parts of θ and φ, θ_frac and φ_frac, following the principle of "the closer the distance, the higher the weight":

[0081] Horizontal weight (θ direction, controlling the influence of left and right pixels):

[0082] The weights of the upper left point and the lower left point are wtheta_left=1-θfrac, where θfrac is the fractional part of the longitude θ, and its value range is [0, 1). When θ is an integer, θfrac=0, and wtheta_left=1.

[0083] The weights of the upper right point and the lower right point are wtheta_right = θfrac. When θ approaches the next integer coordinate, θfrac approaches 1 and wtheta_right approaches 1.

[0084] Longitudinal weight (φ direction, controlling the influence of upper and lower pixels):

[0085] The weights of the upper left point and the upper right point are wphi_top = 1 - φfrac, where φfrac is the fractional part of the latitude φ, and the value range is [0, 1). When φ is an integer, φfrac = 0, and wphi_top = 1.

[0086] The weights of the lower left point and the lower right point are wphi_bottom = φfrac. When φ approaches the next integer coordinate, φfrac approaches 1 and wphi_bottom approaches 1.

[0087] Four-point combination weight (spatial position weighted):

[0088] The weight of the upper left point is w1 = wtheta_left × wphi_top; it represents the joint influence of the left and upper pixels. When θfrac = 0 and φfrac = 0, w1 = 1, and the upper left pixel value is directly taken.

[0089] The weight of the top right point is w2 = wtheta_right × wphi_top; it represents the joint influence of the right and top pixels. When θfrac = 1 and φfrac = 0 (the limit case), w2 = 1.

[0090] The weight of the lower left point is w3 = wtheta_left × wphi_bottom; it represents the joint influence of the left and bottom pixels. When θfrac = 0 and φfrac = 1 (the extreme case), w3 = 1.

[0091] The weight of the bottom right point is w4 = wtheta_right × wphi_bottom, representing the combined influence of the right and bottom pixels. When θfrac = 1 and φfrac = 1 (the limiting case), w4 = 1. Regardless of the values ​​of θfrac and φfrac, w1 + w2 + w3 + w4 = (wtheta_left + wtheta_right) × (wphi_top + wphi_bottom) = 1 × 1 = 1.

[0092] Extract the RGB values ​​of four adjacent pixels:

[0093] Upper left pixel value: RGB1 = (R1, G1, B1);

[0094] Upper right pixel value: RGB2 = (R2, G2, B2);

[0095] Lower left pixel value: RGB3 = (R3, G3, B3);

[0096] Lower right pixel value: RGB4 = (R4, G4, B4);

[0097] The RGB value of the target pixel is obtained by weighted summation:

[0098] R=w1×R1+w2×R2+w3×R3+w4×R4;

[0099] G=w1×G1+w2×G2+w3×G3+w4×G4;

[0100] B=w1×B1+w2×B2+w3×B3+w4×B4.

[0101] This process smoothly transitions the colors of adjacent pixels to avoid jagged edges or color mutations caused by direct rounding, ensuring the visual continuity of the projected flat image.

[0102] Polar distortion preprocessing:

[0103] The range of latitude φ is [- , ], corresponding to the vertical coordinate of the plane image (nadir to zenith). When |φ|>80° (about 1.396 radians), it is defined as the polar region (top φ>80° or bottom φ<-80°), and this region has severe vertical compression in the equirectangular projection:

[0104] Normal area (e.g., |φ| ≤ 60°): The number of plane pixels corresponding to a unit latitude is approximately H / 6 (H is the image height);

[0105] Polar regions (|φ|>80°): The number of plane pixels corresponding to a unit latitude exceeds H / 2, resulting in a pixel density in the polar regions that is more than three times that of normal areas, forming a "pixel pile-up" phenomenon.

[0106] The dense pixels in the polar regions lead to large differences in pixel spans of the same physical size in the three-dimensional space, misleading feature matching and depth estimation in subsequent SfM. The uneven spatial distribution of pixels in the polar regions can cause deviations in gradient calculations during rendering, resulting in an imbalance in the density distribution of Gaussian modeling.

[0107] Convert the vertical pixel coordinate y of the plane image to the latitude φ:

[0108] Top polar region: when φ>80°, y>H×(0.5+80° / 180°)≈H×0.944;

[0109] Bottom polar region: When φ<-80°, y <H×(0.5-80° / 180°)≈H×0.056。

[0110] The polar region is divided into N×N grids (e.g., N=16), where each grid corresponds to a small area of ​​the spherical polar region, to facilitate local scaling operations.

[0111] Compress the polar latitude φ from [80°, 90°] to [70°, 80°] to alleviate the problem of dense vertical pixels. The scaling factor k(φ) is defined as: , ( is the target latitude, such as =80°). For each pixel in the polar region, vertical scaling is performed according to k(φ):

[0112] Original pixel coordinate (y, x) → new coordinate (y', x), where y' = y + (y - y_p) × (1 - k(φ)), y_p is the starting pixel coordinate of the polar region;

[0113] After scaling, bicubic interpolation is used to fill blank pixels to maintain image smoothness. A transition band of 50-100 pixels is set at the boundary between polar and non-polar regions (e.g., y ≈ 0.944H and y ≈ 0.056H). A weighted blending approach (e.g., weights varying linearly with y) is used to avoid noticeable boundary discontinuities after scaling. For the projected planar image, the corresponding latitude φ is calculated for each pixel. If |φ| > 80°, the pixel is marked as a polar region (e.g., mask = 1, non-polar regions are mask = 0). Based on the projection mapping, for each pixel (y, x), φ = (y / H − 0.5) × 180° (converted to degrees) is calculated. If |φ| > 80°, the pixel is marked as a polar region. A binary mask matrix of the same size as the image is generated, or the pixel coordinate range of the polar regions is stored (e.g., y∈[y_top_start, H] for the top polar region and y∈[0, y_bottom_end] for the bottom polar region).

[0114] Step 102: Define the standard resolution (e.g., width W0 = 2048, height H0 = 1024), and calculate the scaling ratio based on the resolution of the original projected image (W1, H1):

[0115] Horizontal scaling ratio s_w= , vertical scaling ratio s_h= ;like ≠ , you need to choose either proportional scaling (maintaining aspect ratio) or cropping to fill:

[0116] Proportional scaling: Take s = min(s_w, s_h), first scale the image by s, then fill the edges with black borders or repeat pixels to the standard size;

[0117] Crop and fill: Scale by s_w and s_h respectively, crop the excess part or fill the blank, and ensure that the final size is (W0, H0).

[0118] Assume the original image has a width of W1 and a height of H1, and the target image has a width of W2 and a height of H2. Calculate the horizontal scaling factor s_w = W1 / W2 and the vertical scaling factor s_h = H1 / H2. For any pixel (x0, y0) in the target image (x0 ranging from 0 to W2-1, y0 ranging from 0 to H2-1), calculate its corresponding position in the original image: horizontal coordinate x1 = x0 × s_w. The resulting x1 is a floating-point value with a decimal. For example, if x0 = 10 and s_w = 1.5, then x1 = 15.2. The vertical coordinate y1 = y0 × s_h. Similarly, y1 is also a floating-point value. For example, if y0 = 5 and s_h = 1.2, then y1 = 6.4.

[0119] Find the integer part of x1 and y1 and determine the 4×4 pixel range around it:

[0120] Take the integer part of x1 as ix = floor(x1), then the four pixel positions in the horizontal direction are ix-1, ix, ix+1, and ix+2; take the integer part of y1 as iy = floor(y1), and the four pixel positions in the vertical direction are iy-1, iy, iy+1, and iy+2; if ix-1 is less than 0, set it to 0 (or use mirror symmetry, for example, map it to 1 when ix-1 = -1); if ix+2 exceeds the original image width W1-1, set it to W1-1 to ensure that the neighborhood does not exceed the image range.

[0121] For each neighborhood pixel point (i, j) (i is the horizontal position, j is the vertical position), calculate its distance to the target point (x1, y1):

[0122] Horizontal distance dx = |x1-i|, for example, if x1 = 15.2 and i = 15, then dx = 0.2; if i = 16, dx = 0.8;

[0123] Vertical distance dy = |y1-j|, ​​for example, y1 = 6.4, j = 6, dy = 0.4; j = 7, dy = 0.6;

[0124] According to the logic of the bicubic kernel function, the closer the distance, the greater the weight of the pixel, and the weight decays with the distance:

[0125] When the distance is between 0 and 1, the weight decreases rapidly as the distance increases;

[0126] When the distance is between 1 and 2, the weight decreases more slowly, but still contributes;

[0127] Pixels with a distance greater than 2 are considered to have a weight of 0 and are not included in the calculation (actually only the 4×4 neighborhood is considered, and the distance is at most 2).

[0128] For each neighborhood pixel (i, j), a weight is calculated based on its dx and dy values, then multiplied by the pixel's RGB value (or RGBA). The weighted values ​​of all 16 pixels are then summed. For example, if a pixel in the neighborhood has an R value of 200 and a weight of 0.3, and another pixel has an R value of 150 and a weight of 0.2, the sum is 200 × 0.3 + 150 × 0.2. This calculation is performed separately for the three RGB channels to obtain the R, G, and B values ​​of the target pixel. If the sum of all weights is not 1 (possibly due to edge processing), the accumulated result is divided by the total weight to ensure that the pixel value is within a reasonable range (e.g., 0-255).

[0129] Regardless of whether the original image is in JPG, BMP or other formats, it needs to be converted using image processing tools:

[0130] Read the original image data and analyze information such as width, height, and color channels;

[0131] Create a file structure in PNG format and set the compression method (such as lossless compression);

[0132] Write the pixel data of the original image to a PNG file, ensuring that pixel values ​​are not lost.

[0133] The unified color space is RGB or RGBA.

[0134] Color space detection and conversion:

[0135] If the original image is in CMYK format (common in printed images), it needs to be converted to RGB: using a specific conversion matrix, the four-channel values ​​of cyan, magenta, yellow, and black are converted to three channels of red, green, and blue. If it is in YCbCr format (common in video encoding), the RGB values ​​are calculated based on the brightness and color difference channels. If the image contains an transparency channel (Alpha), it is retained and converted to RGBA format. If not, a fully transparent Alpha channel (value 255) is added. Check the bit depth of each color channel (such as 16-bit, 32-bit, etc.) and convert it to 8-bit:

[0136] For 16-bit images, the value range of each channel is 0-65535, which is divided by 256 and rounded to 0-255 (for example, 65535÷256≈256, so 255 is taken);

[0137] For a 4-bit image, each channel has only 16 possible values ​​(0-15), which are multiplied by 17 to expand to 0-255 (e.g. 15×17=255);

[0138] Ensure that the bit depth of each channel is consistent at 8 bits after conversion.

[0139] Check whether the image contains a color profile (such as AdobeRGB, ProPhotoRGB, etc.). If a non-sRGB profile exists, use color management tools to convert the image data to the sRGB color gamut. For example, AdobeRGB has a wider color gamut, so colors outside the sRGB range need to be mapped to the nearest displayable color during conversion. After the conversion is complete, remove the original profile and embed a standard sRGB profile to ensure that all images display consistently across different devices.

[0140] Multi-view capture with a panoramic camera ensures that all 360° spatial information of the scene is fully captured, improving the accuracy and completeness of 3D point cloud reconstruction. Equirectangular projection converts spherical images into a two-dimensional plane, establishing a linear mapping between spherical perspective and plane pixels. Polar distortion preprocessing lays the foundation for subsequent distortion correction, reducing the impact of projection errors on reconstruction. Resolution normalization eliminates computational bias caused by resolution differences between images, ensuring consistency of input data during model training and improving reconstruction efficiency and accuracy.

[0141] In a preferred embodiment of the present invention, step 2, performing structure-from-motion processing on the encoded dataset to generate six-degree-of-freedom camera pose parameters and corresponding sparse three-dimensional point clouds for each image, may include:

[0142] Step 200: extracting panorama-invariant feature points from the images in the unified format encoding dataset, and filtering matching point pairs based on spherical view similarity to obtain filtered matching point pairs;

[0143] Step 201, based on the screened matching point pairs, solve the camera pose using an incremental structure-from-motion algorithm to generate initial 3D point coordinates, specifically including:

[0144] Step 2010: Using the camera position of the first image as the origin of the world coordinate system, initialize the first camera pose based on the filtered matching point pairs, and mark the corresponding view as a registered view;

[0145] Step 211: Select a new view that has common view feature points with the registered view from the remaining unregistered views, and use the PnP algorithm to solve the rotation matrix and translation vector of the new view relative to the world coordinate system to obtain the six-degree-of-freedom camera pose;

[0146] Step 2022 , based on the new view pose and the common view feature points, perform a triangulation operation on the new view and the registered view to generate initial 3D point coordinates bound to the world coordinate system;

[0147] Step 202 , using bundle adjustment to jointly optimize camera pose and 3D point coordinates, verifying the panoramic projection consistency of each 3D point, and outputting an optimized camera pose parameter set and a sparse 3D point cloud that satisfies consistency constraints.

[0148] In this embodiment of the present invention, each planar image is converted to spherical coordinates, and a three-dimensional spatial mapping is established using longitude θ (horizontally 0-360°) and latitude φ (vertically -90°-90°) to ensure panoramic perspective continuity during feature extraction. For spherical coordinate images, four scale spaces are constructed horizontally and vertically (e.g., original scale, 1 / 2, 1 / 4, and 1 / 8). Each scale space is smoothed using a Gaussian convolution kernel (standard deviation σ = 1.6) and then subtracted from the image below to generate a Difference of Gaussians (DoG) image. Within the three-dimensional neighborhood of the DoG image (a 3×3×3 window of the current layer and the layers above and below), pixel values ​​are compared to identify local maxima and minima, and the locations of feature points are determined. Each extreme point is accurately fitted to sub-pixel precision using a quadratic function. Edge artifacts and low-contrast points are also removed (e.g., curvature is determined using the Hessian matrix, eliminating points with a principal curvature ratio greater than 10).

[0149] A 16×16 local region is delineated in spherical coordinates, centered around the feature point (accounting for longitude cycles, any excess on the right is joined from the left). This region is then divided into 4×4 sub-blocks. For each sub-block, eight gradient histograms are calculated (based on pixel gradients in the θ and φ directions). This ultimately generates a 128-dimensional descriptor. This descriptor generation incorporates the principal direction of spherical coordinates (the principal direction of the gradient of the surrounding area, centered around the feature point, serves as the descriptor's orientation reference) to ensure rotation invariance under panoramic viewing angles. For all feature points between each pair in the image, the Euclidean distance between the descriptors is calculated, and candidate point pairs with a distance less than 40% of the maximum descriptor distance are retained (for example, if a pair's descriptor distance is 0.8 and the maximum distance is 2.0, it is retained).

[0150] The plane coordinates of the candidate point pairs are back-calculated to spherical coordinates (θ1, φ1) and (θ2, φ2). The angular distance between the two points on the sphere is calculated. If the angular distance is less than 15° (approximately 1 / 3 of the camera's field of view), the two points are considered to be in close proximity within the spherical view and retained as valid candidates. Four candidate point pairs are randomly selected and assumed to be correct matches. The fundamental matrix is ​​calculated using the 8-point method and then solved using SVD decomposition. This fundamental matrix is ​​used to verify all candidate point pairs. The epipolar distance (i.e., the perpendicular distance from a 2D point to the epipolar line) is calculated for each point pair. The number of inliers with an epipolar distance less than 1 pixel is counted. The fundamental matrix with the largest number of inliers is selected, and all inliers (epiperipheral distance ≤ 1 pixel) are retained as the final matching point pairs.

[0151] In step 210, the encoded dataset is traversed, and the sharpness (e.g., using the Laplacian operator variance) and exposure uniformity (e.g., pixel value standard deviation) of each image is calculated. The image with the highest sharpness and most uniform exposure is selected as the first reference view. (If multiple images have similar metrics, the one with the highest sequence number is selected.) The camera optical center of the reference view is set to the world coordinate system origin (0, 0, 0). The optical axis direction is defined as the positive z-axis, with the x-axis pointing horizontally to the right of the image and the y-axis pointing vertically upward, forming a right-handed coordinate system.

[0152] Camera intrinsic parameter initialization:

[0153] Assume that the camera intrinsic parameters are: focal length f = image width / 2 (for example, when the width is 2048 pixels, f = 1024), principal point (cx, cy) = (image width / 2, image height / 2), and distortion coefficient is set to 0 (subsequently optimized by bundle adjustment).

[0154] Pose markers:

[0155] Assign pose parameters to the reference view: the rotation matrix R is a 3×3 identity matrix, the translation vector t is (0, 0, 0), and it is marked as a “registered view” in the data structure, and its pose parameters and feature point coordinates are recorded.

[0156] Step 211: For all feature points in the registered view, an index structure (such as a KD tree or hash table) is created based on the descriptor to facilitate fast query. Each feature point stores its descriptor vector, spherical coordinates (θ, φ), and plane coordinates (u, v) in the image. For each unregistered view, all its feature points (assuming the number is N) are traversed. For each feature point F_unreg:

[0157] Extract its descriptor vector D_unreg and find the most similar feature point in the feature point index of the registered view (such as through K-nearest neighbor search, K=2).

[0158] Calculate the distance (such as Euclidean distance) between D_unreg and the candidate feature point descriptor. If the ratio of the closest distance to the second closest distance is less than 0.7 (empirical threshold), the matching point F_reg is considered to be found.

[0159] Establish a set of matching point pairs to avoid duplicate counting of the same 3D point (e.g., by using the hash value of the feature point's local descriptor to remove duplicates). Count the number of valid matching point pairs between each unregistered view and the registered view, recording this as the number of common view feature points. Traverse all unregistered views and select those with ≥20 common view feature point pairs (20 pairs is an empirical threshold to ensure pose solution stability). If multiple views meet these criteria, select the one with the largest number of common view feature points as the new view. If all views have fewer than 20 common view feature point pairs, terminate the incremental SfM process.

[0160] PnP algorithm pose solving process:

[0161] If there are no triangulated 3D points yet, virtual 3D points are generated based on the feature points of the registered view, assuming that the depth of all matching points is 1 meter. For example, the plane coordinates (u, v) of the feature point F_reg are back-projected into a ray in the camera coordinate system through the intrinsic parameter K, and the point on the ray 1 meter away from the optical center is taken as the virtual 3D point P_virtual = (X, Y, 1). Use the 3D points generated by the previous triangulation to ensure that each 3D point is observed by at least two views and the reprojection error is less than 2 pixels. Select at least 6 pairs of 3D-2D matching points (P_i, g_i) from the common view feature point pairs, where P_i is the 3D point coordinate and g_i is the 2D feature point coordinate in the new view. Convert the 3D point P_i from the world coordinate system to the camera coordinate system of the registered view:

[0162] The pose of the registered view is (R_reg, t_reg), then the coordinates in the camera coordinate system are P_cam = R_reg × P_i + t_reg (assuming that the origin of the world coordinate system is the optical center of the registered view, so t_reg = 0, simplified to P_cam = R_reg × P_i). For each matching point pair (P_cam, g_i), g_i = (u, v), representing the projection position of the 3D point on the image plane, the internal parameter K = (f_x, f_y, c_x, c_y), the projection equation is established as ; Arranged into a linear equation form as A×x=0, where x is the camera pose parameter vector; A is the coefficient matrix (each row corresponds to the constraint of a matching point); and is the direction vector component of the three-dimensional point in the camera coordinate system (i.e., the normalized horizontal and vertical directions); c_x and c_y are the coordinates of the principal point in the camera intrinsic parameters. Perform singular value decomposition on the coefficient matrix A to obtain the singular vector corresponding to the minimum singular value, which is the initial solution of the pose parameters. This solution contains the rotation matrix R and the initial estimate of the translation vector t (R may not be orthogonal at this time, and t is the scaled vector). Since the rotation matrix R obtained by SVD decomposition may not be an orthogonal matrix (i.e., ×R≠I), it is necessary to perform orthogonalization through QR decomposition, that is, decomposing R into an orthogonal matrix Q and an upper triangular matrix R', taking Q as the orthogonalized rotation matrix to ensure the geometric correctness of the rotation. is the transpose of the rotation matrix; I is the identity matrix. The reprojection error is the Euclidean distance between the projected coordinates of the 3D point in the new view and the actual feature point coordinates. The total error is the sum of the squares of the errors of all matching points.

[0163] Iterative optimization process:

[0164] The initial pose is (R0, t0) obtained by SVD decomposition. Calculate the initial total error E0, add a small perturbation to each pose parameter (rotation quaternion and translation vector), and calculate the Jacobian matrix J (partial derivatives of the error with respect to the parameter). Solve the normal equation ,in is the parameter update amount; represents the Jacobian matrix; λ is the damping factor. During the iteration process, λ can be dynamically changed according to the following strategy:

[0165] When the iteration reduces the error, λ is reduced by a factor of 0.1 to 0.5 (e.g., λ←λ×0.1), approaching the Gauss-Newton method and accelerating convergence;

[0166] When the error increases, λ is enlarged by a factor of 10 to 20 (e.g., λ←λ×10) to enhance the stability of the gradient descent and avoid divergence.

[0167] λ must be a positive number (λ>0), with an upper bound not exceeding (A too large λ will cause the optimization to degenerate into pure gradient descent, which will converge very slowly);

[0168] Update pose parameters: R = R + ΔR (through quaternion interpolation), t = t + Δt;

[0169] Calculate the new total error E1. If |E1-E0|<0.1 pixel or the number of iterations reaches 10, terminate the iteration.

[0170] Step 2012: For each pair of matching feature points (u1, v1) and (u2, v2) between the new view and the registered view:

[0171] Through the internal parameter K1 and pose (R1, t1) of the registered view, (u1, v1) is back-projected into a ray in the camera coordinate system: the ray direction is K1 -1 ×(u1, v1, 1), the starting point is the camera optical center (0, 0, 0). Through the intrinsic parameter K2 and pose (R2, t2) of the new view, (u2, v2) is back-projected into a ray in the camera coordinate system: the ray direction is K2 -1 ×(u2, v2, 1), the starting point is the camera optical center (0, 0, 0), and then converted to the world coordinate system through R2 and t2.

[0172] Find the closest point between two rays in the world coordinate system:

[0173] Construct a least squares problem and find the point P such that the sum of the squares of the distances from P to the two rays is minimized. Solve the coordinates of the point through matrix operations. The specific steps are as follows:

[0174] The two rays are represented as L1(t) = O1 + t × D1 and L2(s) = O2 + s × D2, where O1 = 0 (the optical center of the registered view) and O2 = t2 (the optical center of the new view). D1 and D2 are unit direction vectors. Calculate the projection of vector O2 onto L1 and the cross product of the two ray direction vectors. Solve the linear equations for t and s to obtain the closest point P. Calculate the reprojection error of P in the two views. If the error in any view is greater than 2 pixels, or the depth (the distance from P to the origin) is less than 1 meter (near point interference) or greater than 100 meters (far point noise), the point is discarded. Valid 3D points are added to the sparse point cloud collection, their coordinates recorded, timestamped, and associated with the corresponding indexes of the two views (for error calculation during subsequent bundle adjustment).

[0175] Step 202: For each registered view, store the rotation matrix R (represented as a quaternion to avoid gimbal lock) and translation vector t, for a total of 6 degrees of freedom. Each point stores the world coordinates (x, y, z) and the corresponding camera intrinsic parameters (shared global intrinsic parameters or independent intrinsic parameters for each view). For each 3D point P and its corresponding view i, calculate the reprojected coordinates of P (u_i, vi_i) = Ki × (Ri × P + ti) using the pose (Ri, ti) and intrinsic parameters Ki of view i. Calculate the pixel distance error from the actual feature point coordinates (u_i', vi_i'). The total error is the sum of the squared errors of all point-view pairs.

[0176] Using the sparse BundleAdjustment algorithm (e.g. using an incremental solver):

[0177] For each error term, calculate its partial derivative with respect to the camera pose and 3D point coordinates to form a sparse matrix (only non-zero elements are retained).

[0178] Solving normal equations ×J×Δx= ×b, where Δx is the variable update, solved using Cholesky decomposition or preconditioned conjugate gradient method. The camera pose is updated using quaternion interpolation, and the update is directly added to the 3D point coordinates. The process of constructing the Jacobian matrix, solving the normal equations to obtain the update, and updating the camera pose and 3D point coordinates is repeated until 20 iterations have been performed or the total error change is less than 0.01 pixel. At this point, the mean reprojection error can be reduced to less than 1 pixel.

[0179] For each optimized 3D point P:

[0180] Calculate the reprojection error for all relevant views, convert the pixel error to spherical angular error: angular error ≈ pixel error / focal length (radians), and then convert it to degrees. Calculate the mean and standard deviation of the angular error for all views. If the mean is greater than 5° or the standard deviation is greater than 3°, it indicates that the point is inconsistently projected under the panoramic view and may be a mismatched point, so it should be removed. The optimized camera pose parameter set (quaternion and translation vector for each view) and sparse 3D point cloud (coordinates and corresponding view list) are generated.

[0181] Extracting panoramic invariant feature points and screening matching point pairs can reduce the mismatch rate caused by polar distortion and large viewing angle spans in panoramic images, improve the accuracy of subsequent pose calculations, and provide a reliable feature matching foundation for 3D reconstruction. The incremental SfM algorithm adaptively constructs the 3D framework of a panoramic scene by gradually adding new views and solving poses, avoiding the computational complexity explosion caused by processing large amounts of data at once; the initial 3D point cloud generated by triangulation provides the basic spatial structure of the scene. Bundle adjustment jointly optimizes camera pose and 3D points to eliminate accumulated errors and improve global consistency; panoramic projection consistency verification ensures the projection accuracy of 3D points at all viewing angles, avoids reconstruction deviations in polar and edge regions, and generates a more accurate sparse point cloud structure.

[0182] In a preferred embodiment of the present invention, the above step 3 initializes a 3D Gaussian primitive set based on the sparse 3D point cloud, assigns a position, rotation parameter, initial scale value, color parameter, and opacity parameter to each 3D point, and when there is a historical reconstruction result that coincides with the current scene, fuses the historical reconstruction result with the current Gaussian primitive set through pose registration to obtain a fused Gaussian primitive set, which may include:

[0183] Step 300 , converting the sparse 3D point cloud into an initial Gaussian primitive set, where each Gaussian primitive inherits the position and color parameters of the 3D point, and initializes the rotation parameters to no rotation, the scale to a microcube, and the opacity to 0.8;

[0184] Step 301: When a historical reconstruction model is detected, feature matching is performed based on the initial Gaussian primitive set, and the reprojection error of the historical Gaussian primitives is calculated using the camera pose parameters, and invalid Gaussian primitives with an error of more than 2 pixels are eliminated;

[0185] Step 302 : Merge the verified valid historical Gaussian basis elements with the initial Gaussian basis element set, and output a final fused Gaussian basis element set.

[0186] In an embodiment of the present invention, sparse 3D point cloud data is first read from a storage device (such as a hard drive or memory). This data is stored as a collection of points, each containing 3D coordinates (X, Y, Z) and color information (RGB values). Statistical methods are then used to identify and remove outliers. For example, the average distance from all points to their K nearest neighbors is calculated. If a point's distance exceeds three times the global average distance, it is considered an outlier and removed. Alternatively, the RANSAC algorithm is used to randomly select points to fit a planar or spatial model, and points that do not conform to the model are considered outliers. For each 3D point, a Gaussian primitive is created, and the 3D coordinates (X, Y, Z) of the point are directly assigned to the center position of the corresponding Gaussian primitive, ensuring that the Gaussian primitive's spatial location is consistent with the original point. The RGB color value of the 3D point is directly assigned to the Gaussian primitive. For example, if the point's color is (255, 0, 0), the color of the corresponding Gaussian primitive is also set to red, preserving the color characteristics of the original point cloud. Initialize to an unrotated state, specifically represented by the unit quaternion (1, 0, 0, 0) or a 3×3 identity matrix. This means that the Gaussian primitives are initially rotated and maintain an isotropic spherical or cubic shape. Initialize the Gaussian primitives to a tiny cube with a side length of 0.01 in each dimension (the unit can be adjusted based on the scene scale, such as 0.01 meters), so that they resemble tiny point lights when visualized. Set the opacity of all Gaussian primitives to 0.8, giving them a semi-transparent effect when rendered. Encapsulate the parameters of each Gaussian primitive (position, rotation, scale, color, and opacity) into a separate data unit (such as a structure or object), which contains storage and access interfaces for the five parameter sets. Sequentially store all Gaussian primitive data units into a collection (such as an array, linked list, or hash table) to form the initial Gaussian primitive set.

[0187] Step 301 extracts key features of the current scene (such as landmark objects and geometric structures) and matches them with features in the historical scene library. If the matching degree exceeds a certain threshold (such as 80%), it is considered that there is overlap. If the camera pose changes slightly (for example, the translation distance is less than 0.5 meters and the rotation angle is less than 15 degrees), and the time interval between the current frame and the historical frame is short, it is judged to be a continuation of the same scene. If overlap is detected, the corresponding historical Gaussian primitive set is loaded from the database or file storing the historical data. This set contains the historically reconstructed Gaussian primitives and their parameters.

[0188] For the current Gaussian primitive set, the center position of each Gaussian primitive is selected as a feature point, and its local geometric features (such as surface normals and curvature) and color features (such as RGB histograms) are extracted. For the historical Gaussian primitive set, feature points and their feature information are similarly extracted. A descriptor matching algorithm (such as SIFT descriptor matching based on Euclidean distance) is used to establish a correspondence between the feature points of the current and historical Gaussian primitives. For example, the distance between the descriptors of two feature points is calculated. If the distance is less than a threshold (such as 0.6 times the descriptor length), it is considered a matching point pair. Using the matching point pairs, the iterative closest point algorithm (ICP) is used to calculate the transformation relationship between the current camera pose and the historical scene. The specific operation is to transform the points in the historical Gaussian primitive set to the current coordinate system using the transformation matrix (rotation matrix R and translation vector t) so that the two are spatially aligned. For each historical Gaussian primitive, its center position is projected onto the current image plane based on the current camera's intrinsic parameter matrix (focal length, principal point) and extrinsic parameter matrix (pose parameters). The resulting projected point coordinates (u, v) are then calculated. The pixel distance between the projected point (u, v) and the current observation point (e.g., the location of the corresponding feature point) is calculated. An error threshold (e.g., 2 pixels) is set. If the reprojection error of a historical Gaussian primitive exceeds this threshold, it is considered mismatched with the current scene and removed from the set. Only valid Gaussian primitives with an error within the threshold are retained to ensure geometric consistency between historical information and current data.

[0189] In step 302, the current Gaussian basis set and the retained historical Gaussian basis set are traversed. For any two Gaussian basis sets (one from the current set and one from the historical set), the spatial distance between their center positions is calculated. If the distance is less than a threshold (e.g., 0.05 meters), they are considered to represent the same spatial location and are merged. A weighted average is used to calculate the new location, with the weight determined by the number of observations or opacity of the Gaussian basis set. For example, if the opacity of the current Gaussian basis set is 0.8 and the opacity of the historical Gaussian basis set is 0.7, a weighted average is applied to the RGB color values ​​to retain more reliable color information. The rotation parameter can be selected to retain the one with the smaller variance (i.e., the more stable rotation state). The scale parameter is taken as the average or larger value based on the relative size of the two. The opacity is weighted averaged to ensure reasonable opacity after fusion. Gaussian basis sets whose center distance exceeds the threshold are directly added to the final set without merging to avoid accidentally deleting useful information.

[0190] Gaussian primitives can represent complex scenes in a compact parametric form, offering better continuity and visual quality than traditional point clouds. By integrating historical reconstruction results, duplication of work is avoided and reconstruction efficiency is improved. This approach is particularly suitable for long-term monitoring or incremental mapping scenarios. The reprojection error filtering mechanism ensures that only reliable historical information is retained, improving the accuracy and consistency of the final model. The multi-parameter (position, rotation, scale, color, opacity) nature of Gaussian primitives enables them to adapt to different scenario requirements, such as dynamic scenes or semantic enhancement. This method can be combined with other 3D reconstruction techniques (such as neural radiance fields) to further enhance the quality and flexibility of scene representation.

[0191] In a preferred embodiment of the present invention, step 4 above, performing a distortion correction operation on the top or bottom polar region of the panoramic image sequence, and reconstructing the image set by local spherical remapping or cubic projection, is performed. Iterative training is performed based on the fused Gaussian basis set and the corrected image set. During the training process, densification control is performed on the Gaussian basis set based on the spatial gradient of the panoramic screen. When the iteration limit is reached, the training is terminated to obtain the optimized Gaussian basis set. This may include:

[0192] Step 400 , reconstructing the top or bottom polar region of the panoramic image sequence using local spherical remapping or cubic projection to output a corrected image set;

[0193] Step 401, based on the fused Gaussian basis set, corrects the image set and performs iterative training. When the number of iterations reaches an upper limit, the training is terminated and the optimized Gaussian basis set is output. Specifically, the following steps are performed:

[0194] Step 4010: using the rectified image set as training input data to initialize the parameters of the fused Gaussian basis set;

[0195] Step 4011: randomly select a training view and a corresponding camera pose from the rectified image set, and generate a simulated view through an equirectangular projection rendering pipeline based on the camera pose and the current Gaussian basis set;

[0196] Step 4012: Calculate the loss function between the simulated view and the real training view, and update the position, rotation, scale, color, and opacity parameters of the Gaussian primitives by backpropagation according to the loss function to obtain an updated Gaussian primitive set.

[0197] Step 4013: Based on the updated Gaussian basis set, the gradient magnitude of each Gaussian basis in the panoramic screen space is calculated. If the gradient magnitude is greater than a splitting threshold, a splitting operation is performed to generate a sub-Gaussian. If the gradient magnitude is less than a cloning threshold, a cloning operation is performed. If the gradient magnitude is equal to the cloning threshold, the original state is maintained.

[0198] Step 4014: When the number of iterations reaches a preset upper limit, the training is terminated and the optimized Gaussian basis set is output.

[0199] In this embodiment of the present invention, for each frame of a panoramic image sequence, latitude mapping is used to determine the polar regions (e.g., the region with φ > 80° at the top or φ < -80° at the bottom). This region exhibits a high density of vertical pixels under equirectangular projection, resulting in distorted object shapes (e.g., straight lines becoming curved). The spatial distribution of polar region pixels is characteristic of this: per unit angle, the number of pixels occupied by the spherical region on a planar image is more than three times that of the central region, necessitating a coordinate transformation to redistribute pixel density.

[0200] Local spherical remapping process:

[0201] The polar regions are divided into M×N grids (e.g., 64×64), each grid corresponding to a local region of the spherical polar region. The original latitude φ and longitude θ range of each grid is recorded. For each pixel within the grid, the latitude φ is compressed from [80°, 90°] to [70°, 80°] by setting φ' = 70° + (φ - 80°) × 0.5, reducing the number of vertical pixels in the polar region by 50% while keeping the longitude θ unchanged. This transformation creates blank areas between grids. Bicubic interpolation is used to calculate the RGB values ​​of these blank areas based on the pixel values ​​of adjacent grids to ensure a smooth and seamless image. The spherical region corresponding to the polar region is projected onto a face of a cube (e.g., the top or bottom face). Each face of the cube is aligned with the longitude and latitude of the sphere. For example, the center of the top face corresponds to the zenith (φ = 90°), and the four edges correspond to the intersections of φ = 80° and θ = 0°, 90°, 180°, and 270°. The top surface of the cube is expanded into a square plane, with each side being half the height of the original image. This ensures that polar pixels are evenly distributed within the square and eliminates vertical compression. When splicing the cube-projected image with the non-polar image, a 100-pixel-wide transition band is set at the junction. Linear weighted fusion is used (for example, the first 50 pixels are the cube-projected pixel values, and the last 50 pixels are the original image pixel values) to avoid splicing artifacts.

[0202] In step 4010, the rectified image set is arranged in view order. The camera pose parameters (rotation matrix R, translation vector t) and intrinsic parameters (focal length f, principal point (cx, cy)) of each image are extracted and stored as training data pairs (image, pose). The position, rotation, scale, color, and opacity parameters of the fused Gaussian basis set are used as initial values, where:

[0203] The position remains the fused coordinates, the rotation is reset to the unit quaternion, the scale is adjusted to a uniform value of 0.02m×0.02m×0.02m, the color is the fused RGB value, and the opacity is set to 0.7 (reducing the initial opacity facilitates iterative optimization). A preset iteration limit (e.g., 500) is set, and the learning rate is initialized to 0.01, which is then exponentially decayed with the number of iterations (e.g., the learning rate is multiplied by 0.5 every 100 iterations).

[0204] Step 4011: Randomly select a frame from the rectified image set as the current training view, and obtain its camera pose (R_cam, t_cam) and intrinsic parameter K. For each Gaussian basis, calculate its coordinates in the camera coordinate system: P_cam = R_cam × P_world + t_cam (P_world is the world coordinate).

[0205] Convert the 3D Gaussian to a 2D image using the equirectangular projection pipeline:

[0206] Calculate the longitude θ and latitude φ of the Gaussian center: θ = atan2(P_cam.y, P_cam.x), φ = asin( ),in, The coordinates of a 3D point in the camera coordinate system are represented by a 3D vector (P_cam.x, P_cam.y, P_cam.z), which is converted to plane coordinates x = (θ / 2π) × W, y = (φ / π+0.5) × H (W, H are the image width and height). Based on the scale, rotation, and opacity of the Gaussian, an elliptical area is rendered at the plane coordinates. The color is the RGB value of the Gaussian, and the opacity is superimposed according to the parameters (such as alpha blending). The rendering results of all Gaussian primitives are superimposed to generate a simulated view, which has the same size as the training view. Figure 1 To.

[0207] In step 4012, the mean squared error (MSE) between the simulated and trained views is calculated, which is the average of the sum of squared RGB differences for each pixel, focusing on the polar pixels (with a weight of 1.5). A pre-trained CNN is used to extract feature maps from the two images. The cosine similarity of the feature maps is calculated to ensure structural consistency. The total loss is 0.8 × MSE + 0.2 × structural loss. The gradient of the loss with respect to the parameters of the Gaussian basis is calculated as follows:

[0208] Position gradient: Adjust along the camera optical axis to reduce the projection error;

[0209] Rotation gradient: Calculate the effect of rotation change on projection through quaternion perturbation;

[0210] Scale gradient: Increase or decrease the scale to make the Gaussian coverage area match the real object;

[0211] Color gradient: adjust the RGB value to be close to the training view color;

[0212] Opacity gradient: adjust the alpha value to make the rendered brightness consistent with the real image;

[0213] Update parameters according to the gradient direction and learning rate: P_new = P_old - learning rate × gradient, where the learning rate is 0.5 and P_old is the original value of the Gaussian basis parameter before this iterative update. For the rotation parameters, quaternion interpolation is used to update them. After the update, the quaternion is normalized to ensure that its modulus is 1 and maintain the orthogonality of the rotation matrix.

[0214] In step 4013, for each Gaussian basis element, the pixel gradient magnitude of its projected area in the latest simulated view is calculated (e.g., using the Sobel operator to calculate the gradient in the x and y directions, taking the L2 norm). The gradient magnitude reflects the detail richness of the Gaussian surface. Large gradients indicate rich edges or textures, requiring a denser Gaussian representation.

[0215] Split operation (gradient magnitude > split threshold):

[0216] Split the Gaussian into two sub-Gaussians, offset by ±δ (δ is 0.01m, along the gradient) from the center of the original Gaussian. Scale them down to 0.8 times their original size, inheriting the color and opacity of the original Gaussian. Add a small random perturbation to the rotation (e.g., quaternion multiplication (1, ε, ε, ε), where ε = 0.1). Remove the original Gaussian and add the two sub-Gaussians to the set, ensuring the total opacity remains unchanged (e.g., original α = 0.7, sub-Gaussian α = 0.35).

[0217] Clone operation (gradient amplitude < cloning threshold):

[0218] Duplicate the current Gaussian, keeping its position unchanged but scaling it up by 1.2x. Add random noise to the color (e.g., a ±5 RGB offset), and reduce its opacity by 0.1. This will be used to fill sparse areas. The cloned Gaussian complements the original, improving coverage in low-gradient areas. Set the split threshold to 15 (pixel gradient units) and the clone threshold to 5. Leave the intermediate values ​​(5-15) unchanged.

[0219] In step 4014, training is terminated when the number of iterations reaches a preset upper limit (e.g., 500) or the total loss change is less than 0.01 (for 10 consecutive iterations). Gaussians with opacity < 0.1 are removed, and Gaussians with a position distance < 0.02m and a similarity > 0.8 are merged (similarity is calculated based on color and scale). The optimized Gaussian basis set is saved. Each Gaussian basis contains the final position, rotation (quaternion), scale (3D vector), color (RGB), and opacity (α) parameters.

[0220] Polar region distortion correction eliminates the pixel density problem at the top / bottom of the panoramic image, restores the normal shape of polar region objects, makes the pixel distribution uniform, and the iterative training optimizes through the comparison between the simulated view and the real view, improving the position accuracy of the Gaussian basis elements. Screen space gradient control dynamically adjusts the Gaussian density, automatically splits the Gaussian in texture-rich regions (such as edges) to enhance the detail representation ability; clones the Gaussian in smooth regions to avoid holes.

[0221] In a preferred embodiment of the present invention, in step 5 above, performing panoramic rendering on the optimized Gaussian basis element set and outputting a panoramic image or a dense three-dimensional model may include:

[0222] Step 500, inputting the optimized Gaussian basis element set into an equirectangular projection rendering pipeline, and generating an equirectangular projection grid according to the output resolution of the panoramic image;

[0223] Step 501, for each pixel point in the grid, generating a ray emitted from the camera center according to the spherical line-of-sight direction mapped by the planar coordinates;

[0224] Step 502, calculating the intersection weights of the ray with all Gaussian basis elements in the scene, and generating the final color of the current pixel by weighted accumulation of the color values and opacities contributed by each Gaussian;

[0225] Step 503, traversing all pixels to complete the rendering output of the panoramic image, or generating a dense three-dimensional point cloud model based on the spatial position and scale parameters of the optimized Gaussian basis element set.

[0226] In an embodiment of the present invention, an optimized Gaussian basis element set is obtained. Each Gaussian basis element includes a position (X, Y, Z), a rotation quaternion (w, x, y, z), a three-dimensional scale (sx, sy, sz), a color (R, G, B), and an opacity α. The resolution of the output panoramic image (such as 4096×2048 pixels) is determined, and this resolution determines the density of the projection grid.

[0227] Horizontal direction: The longitude range [0°, 360°) is evenly mapped to the image width W pixels, and each pixel corresponds to a longitude increment of 360° / W;

[0228] Vertical direction: The latitude range from 90° south latitude (90°S) to 90° north latitude (90°N) (total span 180°) is evenly mapped to the image height H pixels, and the latitude increment corresponding to each pixel is 180° / H.

[0229] For each pixel (i, j) (0 ≤ i < W, 0 ≤ j < H), calculate its corresponding spherical coordinates:

[0230] Longitude θ = i × ;

[0231] Latitude φ = 90° - j × .

[0232] Create a two-dimensional array to store the spherical coordinates (θ, φ) of each grid point and record its pixel index (i, j) in the image to form an equidistant cylindrical projection grid.

[0233] Step 501, for each grid point (θ, φ), convert its spherical coordinates into a three-dimensional direction vector (vx, vy, vz) on the unit sphere: vx = cos(φ) × cos(θ); vy = cos(φ) × sin(θ); vz = sin(φ).

[0234] Starting from the camera center C, a ray r(e)=C+e×(vx, vy, vz) is generated along the direction vector (vx, vy, vz), where e is a parameter on the ray (e≥0). The ray represents the line of sight from the camera perspective through the pixel (i, j) on the image plane.

[0235] Step 502: transform the ray r(e) from the world coordinate system to the local coordinate system of each Gaussian primitive:

[0236] Apply an inverse rotation (the conjugate of the Gaussian basis rotation quaternion) and a translation (the negative Gaussian basis position) to transform the ray's origin and direction vector to local space. In the local coordinate system, the Gaussian basis is a scaled ellipsoid centered at the origin (with scale parameters (sx, sy, sz)). Calculate the intersection parameters umin and umax of the ray with this ellipsoid; if umin ≤ umax and umin ≥ 0, the ray intersects the Gaussian basis, and record the intersection interval [umin, umax]. For each point u in the intersection interval, calculate its weight in the Gaussian distribution:

[0237] Weight w(u)= , where d(u) is the distance from point u to the center of the Gaussian basis element, and σ is the standard deviation of the Gaussian distribution (related to the scale of the Gaussian basis element). Discretely sample the intersection interval [umin, umax] (e.g., 10 sample points), accumulate the weights of each sample point, and obtain the total intersection weight W of the ray with the Gaussian basis element.

[0238] Calculate the contribution of each Gaussian primitive to the final color: Contribution Color ;

[0239] Cumulative Opacity: ;

[0240] Update final color: .

[0241] In step 503, we traverse each pixel (i, j) in the equirectangular projection grid, repeating steps 501 and 502 to calculate the final color of each pixel. The color values ​​of all pixels are then applied to the corresponding positions in the output image to form a complete panoramic image. For each Gaussian primitive, N points (e.g., N = 100) are uniformly sampled on its surface. The locations of the sampling points are determined by the position, rotation, and scale of the Gaussian primitive, and each sampling point inherits the color (R, G, B) and opacity α of the Gaussian primitive. The sampling points of all Gaussian primitives are merged, removing duplicate points or points that are too close together (e.g., a distance threshold of 0.01 meters) to form the final dense 3D point cloud model.

[0242] Equirectangular projection ensures complete coverage of the 360°×180° viewing angle without any stitching. Weighted rendering based on Gaussian primitives preserves the fine structure of the scene, especially the polar details obtained through iterative optimization. By pre-calculating the projection grid and ray intersection weights, rendering speed is improved compared to traditional ray tracing. Each Gaussian primitive generates multiple sampling points, and the point cloud density is improved compared to the original sparse point cloud. Based on the optimized Gaussian primitive parameters, the geometric accuracy of the point cloud reaches the centimeter level and can be directly used for 3D modeling and measurement. From polar correction to Gaussian primitive optimization to rendering, a complete quality improvement chain is formed. The average reprojection error of panoramic images is <1 pixel, and the same Gaussian primitive set can support both panoramic image and 3D point cloud output.

[0243] An embodiment of the present invention further provides a computing device comprising: a processor and a memory storing a computer program, wherein the computer program, when executed by the processor, performs the above-described method. All implementations in the above-described method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0244] The embodiment of the present invention further provides a computer-readable storage medium storing instructions, which, when executed on a computer, causes the computer to execute the above-described method. All implementations in the above-described method embodiment are applicable to this embodiment and can achieve the same technical effects.

[0245] The above is a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as within the scope of protection of the present invention.

Claims

1. A 3DGS-based panoramic image 3D reconstruction method, characterized in that: The method comprises: Step 1: Collect a panoramic image sequence of the scene to be reconstructed, perform unified encoding through equirectangular projection, and obtain a coded data set; Step 2: Perform panoramic structure-from-motion processing on the encoded dataset to generate six-degree-of-freedom camera pose parameters and corresponding sparse three-dimensional point clouds for each image. This includes extracting panoramically invariant feature points from the images in the unified format encoded dataset, screening matching point pairs based on spherical view similarity, and obtaining screened matching point pairs. Based on the screened matching point pairs, the camera pose is solved using an incremental structure-from-motion algorithm to generate initial three-dimensional point coordinates. The camera pose and three-dimensional point coordinates are jointly optimized using bundle adjustment, and the panoramic projection consistency of each three-dimensional point is verified. The optimized camera pose parameter set and a sparse three-dimensional point cloud that meets the consistency constraint are output. Step 3: Initialize a 3D Gaussian primitive set based on the sparse 3D point cloud, assign position, rotation parameters, initial scale, color parameters and opacity parameters to each 3D point, and when there is a historical reconstruction result that coincides with the current scene, fuse the historical reconstruction result with the current Gaussian primitive set through pose registration to obtain a fused Gaussian primitive set, including converting the sparse 3D point cloud into an initial Gaussian primitive set, where each Gaussian primitive inherits the position and color parameters of the 3D point, and initializes the rotation parameters to a non-rotational state, a scale to a microscopic cube, and an opacity to 0.8; when a historical reconstruction model is detected, perform feature matching based on the initial Gaussian primitive set, and use the camera pose parameters to calculate the reprojection error of the historical Gaussian primitives, and eliminate invalid Gaussians with an error of more than 2 pixels; merge the verified valid historical Gaussian primitives with the initial Gaussian primitive set, and output the final fused Gaussian primitive set; Step 4: perform distortion correction on the top or bottom polar region of the panoramic image sequence, and reconstruct the corrected image set using local spherical remapping or cubic projection; perform iterative training based on the fused Gaussian basis set and the corrected image set. During the training process, perform densification control on the Gaussian basis set based on the panoramic screen spatial gradient. When the iteration limit is reached, the training is terminated to obtain an optimized Gaussian basis set, including reconstructing the top or bottom polar region of the panoramic image sequence using local spherical remapping or cubic projection, and outputting a corrected image set; use the corrected image set as training input data to initialize the parameters of the fused Gaussian basis set; randomly select training views and corresponding camera poses from the corrected image set, and perform the training on the Gaussian basis set based on the training set. The camera pose and the current Gaussian primitive set are used to generate a simulated view through an equirectangular projection rendering pipeline. The loss function between the simulated view and the real training view is calculated, and the position, rotation, scale, color, and opacity parameters of the Gaussian primitives are updated by backpropagation according to the loss function to obtain an updated Gaussian primitive set. Based on the updated Gaussian primitive set, the gradient magnitude of each Gaussian primitive in the panoramic screen space is calculated. If the gradient magnitude is greater than the splitting threshold, a splitting operation is performed to generate a sub-Gaussian. If the gradient magnitude is less than the cloning threshold, a cloning operation is performed. If the gradient magnitude is equal to the cloning threshold, the original state is maintained. When the number of iterations reaches the preset upper limit, the training is terminated and the optimized Gaussian primitive set is output. Step 5: Perform panoramic rendering on the optimized Gaussian primitive set to output a panoramic image or a dense three-dimensional model.

2. The method for 3D reconstruction of panoramic images based on 3DGS according to claim 1, characterized in that: Collect a panoramic image sequence of the scene to be reconstructed and uniformly encode it through equirectangular projection to obtain an encoded data set, including: Capturing a multi-view spherical image sequence of a scene by a panoramic camera; Each spherical image is converted into a two-dimensional plane image using equidistant cylindrical projection, and a mapping relationship between spherical view coordinates and plane pixel coordinates is established; All converted images are resolution-normalized to generate a coded dataset in a unified format.

3. The method for 3D reconstruction of panoramic images based on 3DGS according to claim 2, characterized in that: Based on the filtered matching point pairs, the camera pose is solved using the incremental structure-from-motion algorithm to generate the initial 3D point coordinates, including: Taking the camera position of the first image as the origin of the world coordinate system, initialize the first camera pose based on the filtered matching point pairs, and mark the corresponding view as a registered view; From the remaining unregistered views, a new view that has common view feature points with the registered view is selected, and the rotation matrix and translation vector of the new view relative to the world coordinate system are solved using the PnP algorithm to obtain the six-degree-of-freedom camera pose; Based on the new view pose and common view feature points, a triangulation operation is performed on the new view and the registered view to generate the initial 3D point coordinates bound to the world coordinate system.

4. The method for 3D reconstruction of panoramic images based on 3DGS according to claim 3, characterized in that: Perform panoramic rendering on an optimized Gaussian primitive set, outputting panoramic images or dense 3D models, including: The optimized Gaussian primitive set is input into the equirectangular projection rendering pipeline, and an equirectangular projection grid is generated according to the output resolution of the panoramic image. For each pixel in the grid, a ray is generated from the camera center according to the spherical view direction corresponding to the plane coordinate mapping; Calculate the intersection weights of the ray with all Gaussian primitives in the scene, and generate the final color of the current pixel by weighted accumulation of the color values ​​and opacity contributed by each Gaussian. Traverse all pixels to complete panoramic image rendering output, or generate a dense 3D point cloud model based on the spatial position and scale parameters of the optimized Gaussian primitive set.

5. A computing device, characterized in that include: one or more processors; A storage device for storing one or more programs, wherein when the one or more programs are executed by the one or more processors, the one or more processors implement the method according to any one of claims 1 to 4.

6. A computer-readable storage medium, characterized in that The computer-readable storage medium stores a program, which implements the method according to any one of claims 1 to 4 when executed by a processor.

Images