3Dgs rendering method and system for hdr large scene

Abstract

The present application belongs to the field of scene reconstruction, and particularly relates to a 3DGS rendering method and system for HDR large scenes, comprising: acquiring a multi-exposure image sequence and camera pose parameters of a target scene; reconstructing a high dynamic range light field to obtain light field configuration data comprising a dynamic range compression curve, indoor and outdoor area segmentation masks, and a material reflectivity model; iteratively optimizing spherical harmonic coefficients of a Gaussian sphere in a three-dimensional Gaussian splatting model according to the light field configuration data until the area proportion difference between overexposed and underexposed regions is less than a first threshold value; real-time rendering and collecting pixel-level brightness distribution, calculating the average brightness ratio of overexposed and underexposed regions, and if the ratio is less than a second threshold value, directly outputting, otherwise adjusting the spherical harmonic coefficients and tone mapping operators and re-rendering until the requirements are met. Through the light field guided spherical harmonic coefficient optimization and real-time dynamic adjustment, the present application effectively solves the overexposure, underexposure and light and shadow abrupt change problems in HDR scenes, and significantly improves the visual realism and consistency of large scene virtual roaming.

Description

Technical Field

[0001] This invention belongs to the field of scene reconstruction, and in particular relates to a 3DGS rendering method and system for large HDR scenes. Background Technology

[0002] With the deep application of digital twin and virtual reality technologies in cultural heritage protection, smart agriculture monitoring, and intelligent city management, the lighting environment faced by 3D scene reconstruction and real-time rendering is becoming increasingly complex. In high dynamic range (HDR) scenarios such as virtual tours of ancient building complexes, full-cycle monitoring of farmland, and digital displays of large venues, rendering systems typically need to handle dramatic changes in lighting, from direct outdoor sunlight to deep indoor shadows, while simultaneously ensuring consistency of light and shadow under multiple time periods, angles, and light sources, as well as display adaptation for multiple terminals such as mobile devices and VR / AR. The system is prone to issues such as overexposure in bright areas leading to loss of detail, underexposure in shadow areas resulting in unrecognizable textures, and visual breaks caused by abrupt changes in lighting at the boundary between indoor and outdoor environments. Moreover, these problems may dynamically evolve with changes in viewpoint, time, and scene switching.

[0003] In existing technologies, a common approach is to perform lighting compensation based on traditional 3D reconstruction and rendering pipelines. This includes methods such as multi-exposure image fusion, tone mapping, or global illumination pre-calculation, combined with fixed exposure parameters or post-processing enhancement modules, attempting to alleviate uneven lighting. While this approach is effective in single-scene, static lighting conditions, it struggles in large scenes with dynamically changing lighting, significant differences in brightness across multiple regions, and high demands for real-time interaction. Fixed-parameter exposure strategies fail to balance global brightness, post-processing enhancement cannot recover lost texture details, and there is a lack of underlying optimization of the lighting response capabilities of the 3D Gaussian sphere itself. Consequently, the rendering results still exhibit noticeable visual flaws in areas of dramatic lighting changes.

[0004] Another approach employs Neural Radiation Field (NeRF) or its variants for scene representation and rendering, utilizing implicit neural representations to jointly model lighting and geometry. While this method enhances the realism of the rendering to some extent, existing NeRF models still have significant shortcomings in HDR scene adaptation: the finite-band encoding used to reduce computational complexity may weaken the fidelity of high dynamic brightness information, making it difficult to fully preserve texture details in extreme lighting areas while maintaining rendering efficiency; the non-editable nature of implicit representations prevents fine-tuning of parameters for specific lighting areas, and its rendering speed is insufficient to meet the real-time interactive needs of large scenes; furthermore, existing solutions lack differentiated processing for the brightness response characteristics of different display devices when adapting to multiple terminals, resulting in significant differences in the visual effects of the same HDR scene on different terminals.

[0005] Therefore, the technical problem that the existing technology urgently needs to solve is how to achieve vector-level high-fidelity extraction and expression of lighting information of large scenes with high dynamic range under the 3D Gaussian Splash (3DGS) rendering framework, dynamically adjust the lighting response parameters of the Gaussian sphere while identifying areas of dramatic lighting changes, achieve smooth transition of light and shadow through adaptive exposure fusion and indoor and outdoor lighting transition modeling, and combine multi-terminal display adaptation mechanism to ensure visual consistency of rendering effects on different devices, thereby fundamentally solving the problems of glare, dead black and sudden lighting changes in large-scene virtual roaming. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention proposes a 3DGS rendering method and system for large HDR scenes, comprising: acquiring multi-exposure image sequences and camera pose parameters of the target scene; reconstructing a high dynamic range illumination field to obtain illumination field configuration data including dynamic range compression curves, indoor and outdoor region segmentation masks, and material reflectivity models; iteratively optimizing the spherical harmonic coefficients of the Gaussian sphere in the 3D Gaussian splash model until the difference in the area ratio of overexposed and underexposed regions is less than a first threshold; rendering and acquiring pixel-level brightness distribution in real time, calculating the average brightness ratio of overexposed and underexposed regions, and directly outputting the value if it is less than a second threshold; otherwise, adjusting the spherical harmonic coefficients and tone mapping operators and re-rendering until the threshold is met. This invention effectively solves the problems of overexposure, underexposure, and abrupt changes in light and shadow in HDR scenes through illumination field-guided spherical harmonic coefficient optimization and real-time dynamic adjustment, significantly improving the visual realism and consistency of large-scene virtual roaming.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] 3DGS rendering methods for large HDR scenes include:

[0009] Acquire multi-exposure image sequences of the target scene and the corresponding camera pose parameters;

[0010] Based on the multi-exposure image sequence and the camera pose parameters, a high dynamic range illumination field is reconstructed to obtain illumination field configuration data. The illumination field configuration data includes at least a dynamic range compression curve, an indoor / outdoor region segmentation mask, and a material reflectivity model.

[0011] Based on the lighting field configuration data, the spherical harmonic coefficients of each Gaussian sphere in the three-dimensional Gaussian splash model are iteratively optimized until the difference between the area ratio of overexposed region and the area ratio of underexposed region in the current rendering frame is less than the first preset threshold, and a lighting adaptation completion signal is output.

[0012] In response to receiving the lighting adaptation completion signal, the three-dimensional Gaussian splash model is rendered in real time at a preset rendering frame rate, and pixel-level brightness distribution data of each rendering frame is collected during the rendering process. The pixel-level brightness distribution data includes the brightness value of each pixel and an overexposure or underexposure indicator determined based on a brightness threshold.

[0013] Specifically, the 3DGS rendering method also includes:

[0014] Calculate the brightness ratio between the average brightness of the overexposed areas and the average brightness of the underexposed areas in the current rendering frame, and determine whether the brightness ratio is less than a preset second threshold.

[0015] If the brightness ratio is less than the second threshold, the current rendering frame is output based on the current spherical harmonic coefficients and the preset global tone mapping operator.

[0016] If the brightness ratio is greater than or equal to the second threshold, then the exposure fusion weight adjustment and tone mapping curve correction are calculated based on the pixel-level brightness distribution data. The exposure fusion weight adjustment is applied to the spherical harmonic coefficients of the Gaussian sphere to update the spherical harmonic coefficients, and the tone mapping curve correction is applied to the global tone mapping operator to update the global tone mapping operator. Then, the step of performing real-time rendering of the three-dimensional Gaussian splash model at a preset rendering frame rate is returned until the brightness ratio is less than the second threshold.

[0017] Specifically, iterative optimization is performed on the spherical harmonic coefficients of each Gaussian sphere in the 3D Gaussian splash model, including:

[0018] Based on the indoor / outdoor region segmentation mask, the Gaussian spheres in the 3D Gaussian splash model are divided into an indoor region Gaussian sphere set, an outdoor region Gaussian sphere set, and a boundary Gaussian sphere. A first optimization step size is assigned to the indoor region Gaussian sphere set, a second optimization step size is assigned to the outdoor region Gaussian sphere set, and a third optimization step size, which is between the first and second optimization step sizes, is assigned to the boundary Gaussian spheres; wherein the first optimization step size is smaller than the second optimization step size.

[0019] The current viewpoint is selected sequentially from multiple preset viewpoints, and the current rendering frame of the current viewpoint is rendered based on the initial spherical harmonic coefficients obtained from initialization. The pixel radiance value of the current rendering frame is mapped to the display brightness value according to the dynamic range compression curve.

[0020] Based on preset overexposure and underexposure thresholds, each pixel of the current rendering frame is marked with an overexposure or underexposure flag. The first area ratio of the overexposure region and the second area ratio of the underexposure region in the current rendering frame are calculated, and the difference between the first area ratio and the second area ratio is calculated as the current difference.

[0021] Based on the diffuse reflection coefficient and specular reflection coefficient in the material reflectivity model, the spherical harmonic coefficient gradients of the Gaussian sphere in the indoor area, the Gaussian sphere in the outdoor area, and the boundary Gaussian sphere are calculated respectively, and the spherical harmonic coefficients of the corresponding Gaussian spheres are updated using the first optimization step size, the second optimization step size, and the third optimization step size respectively.

[0022] Repeat the step of sequentially selecting the current viewing direction from multiple preset viewing directions until the current difference is less than the first preset threshold in three or more consecutive viewing directions. Then stop the iterative optimization and output a lighting adaptation completion signal.

[0023] Specifically, the Gaussian spheres in the three-dimensional Gaussian splash model are divided into an indoor region Gaussian sphere set, an outdoor region Gaussian sphere set, and a boundary Gaussian sphere, including:

[0024] For each Gaussian sphere, obtain its center point coordinates and the spatial extension range defined by the covariance matrix, wherein the spatial extension range is an ellipsoidal region whose semi-axis length is determined by the eigenvalues ​​of the covariance matrix;

[0025] Based on the spatial extension range, calculate the set of voxels covered by the Gaussian sphere in the indoor and outdoor region segmentation mask, where the set of voxels is all voxels that intersect with the ellipsoidal region;

[0026] The number of indoor voxels and the number of outdoor voxels in the voxel set are counted, and the proportion of indoor voxels and the proportion of outdoor voxels are calculated respectively.

[0027] When the proportion of indoor voxels is greater than or equal to the first division threshold, the Gaussian sphere is assigned to the indoor region Gaussian sphere set; when the proportion of outdoor voxels is greater than or equal to the second division threshold, the Gaussian sphere is assigned to the outdoor region Gaussian sphere set.

[0028] When both the indoor and outdoor voxel proportions are less than the corresponding division thresholds, the Gaussian sphere is marked as a boundary Gaussian sphere. The third optimization step size is obtained by linear interpolation based on the ratio of the indoor voxel proportion to the outdoor voxel proportion.

[0029] Specifically, mapping the pixel radiance value of the current rendered frame to a display brightness value based on the dynamic range compression curve includes:

[0030] Obtain the camera intrinsic and extrinsic parameter matrices corresponding to the current viewpoint direction; for each pixel in the current rendering frame, query all Gaussian spheres that intersect with the light rays of each pixel through ray tracing or Gaussian splashing algorithm, calculate the radiance value of the corresponding pixel based on the initial spherical harmonic coefficients of each Gaussian sphere, and generate a radiance image in the current viewpoint direction.

[0031] The radiometric image is subjected to pixel-level validity verification to detect whether there are invalid radiometric values. Invalid radiometric values ​​include NaN, infinity, or values ​​that exceed the preset physical range. If they exist, the radiometric value of the pixel is replaced with the average radiometric value of the adjacent valid pixels or the preset default radiometric value.

[0032] Read the minimum radiance, maximum radiance, minimum displayable brightness, maximum displayable brightness, and contrast adjustment factor from the dynamic range compression curve;

[0033] For each pixel in the radiometric image, its radiometric value is compared with the minimum and maximum radiometric values. If it is less than the minimum radiometric value, it is set to the minimum radiometric value. If it is greater than the maximum radiometric value, it is set to the maximum radiometric value.

[0034] The cropped radiance value is logarithmically mapped according to the contrast adjustment factor to calculate the corresponding display brightness value, which is between the minimum and maximum displayable brightness. The mapped display brightness value is then output as the pixel brightness value of the current rendering frame.

[0035] Specifically, the exposure fusion weight adjustment and tone mapping curve correction are calculated based on the pixel-level brightness distribution data, including:

[0036] Obtain pixel-level brightness distribution data of the current rendering frame, and calculate the proportion of overexposed pixels to the total number of pixels as the first overexposure ratio, and calculate the proportion of underexposed pixels to the total number of pixels as the first underexposure ratio.

[0037] The ratio of the first overexposure ratio to the first underexposure ratio is calculated as the first scaling factor;

[0038] When the first overexposure ratio is greater than the first underexposure ratio, a first gain coefficient is calculated based on the first ratio factor. The first gain coefficient is equal to the first ratio factor multiplied by the first preset gain, and the first gain coefficient is used as the exposure fusion weight adjustment amount. The exposure fusion weight adjustment amount is used to enhance the weight of the high frequency component in the spherical harmonic coefficient of the Gaussian sphere corresponding to the overexposure area.

[0039] When the first underexposure ratio is greater than the first overexposure ratio, a second gain coefficient is calculated based on the reciprocal of the first ratio factor. The second gain coefficient is equal to the reciprocal of the first ratio factor multiplied by a second preset gain. The second gain coefficient is used as the exposure fusion weight adjustment amount. The exposure fusion weight adjustment amount is used to enhance the weight of the low-frequency components in the spherical harmonic coefficient of the Gaussian sphere corresponding to the underexposure region.

[0040] Calculate the deviation between the brightness ratio of the current rendered frame and the preset second threshold, multiply the deviation by the third preset gain to obtain the tone mapping curve correction amount, which is used to adjust the contrast adjustment factor of the dynamic range compression curve.

[0041] When both the first overexposure ratio and the first underexposure ratio exceed the third ratio threshold, the exposure fusion weight adjustment amount and the tone mapping curve correction amount are calculated simultaneously.

[0042] Specifically, applying the exposure fusion weight adjustment to the spherical harmonic coefficients of the Gaussian sphere includes:

[0043] For each Gaussian sphere in the set of Gaussian spheres in the indoor region, the weight of the 0th to 2nd order low-frequency components in its spherical harmonic coefficients is increased according to the exposure fusion weight adjustment amount, and the increase of the low-frequency components is the exposure fusion weight adjustment amount multiplied by the first region coefficient.

[0044] For each Gaussian sphere in the outdoor area Gaussian sphere set, the weight of the 3rd order and above high-frequency components in its spherical harmonic coefficients is increased according to the exposure fusion weight adjustment amount, and the increase of the high-frequency components is the exposure fusion weight adjustment amount multiplied by the second area coefficient.

[0045] For each Gaussian sphere in the set of boundary Gaussian spheres, the weights of the low-frequency and high-frequency components in its spherical harmonic coefficients are increased simultaneously according to the exposure fusion weight adjustment amount. The increase in the low-frequency component is the exposure fusion weight adjustment amount multiplied by the indoor voxel percentage of the Gaussian sphere, and the increase in the high-frequency component is the exposure fusion weight adjustment amount multiplied by the outdoor voxel percentage of the Gaussian sphere.

[0046] Specifically, applying the tone mapping curve correction to the global tone mapping operator includes:

[0047] Read the current contrast adjustment factor from the current global tone mapping operator, add the contrast adjustment factor to the tone mapping curve correction amount, and obtain the updated contrast adjustment factor;

[0048] The updated contrast adjustment factor is written into the global tone mapping operator to replace the original contrast adjustment factor, while keeping the minimum radiosity, maximum radiosity, minimum displayable brightness, and maximum displayable brightness in the dynamic range compression curve unchanged.

[0049] Return to the step of performing real-time rendering of the three-dimensional Gaussian splash model at a preset rendering frame rate, regenerate the current rendering frame based on the updated spherical harmonic coefficients and the updated global tone mapping operator, and recalculate the brightness ratio.

[0050] The maximum number of iterations is set to Q. If the brightness ratio is less than the second threshold within the maximum number of iterations, the adjustment is terminated and the current rendering frame is output. If the brightness ratio is still not less than the second threshold after the maximum number of iterations, the current rendering frame is forcibly output and alarm information is recorded.

[0051] A 3DGS rendering system designed for large HDR scenes includes:

[0052] The acquisition module is used to acquire multi-exposure image sequences of the target scene and the corresponding camera pose parameters;

[0053] The illumination field reconstruction module reconstructs a high dynamic range illumination field based on the multi-exposure image sequence and the camera pose parameters to obtain illumination field configuration data. The illumination field configuration data includes at least a dynamic range compression curve, an indoor / outdoor region segmentation mask, and a material reflectivity model.

[0054] The spherical harmonic coefficient pre-optimization module is used to perform iterative optimization of the spherical harmonic coefficients of each Gaussian sphere in the three-dimensional Gaussian splash model according to the illumination field configuration data, until the difference between the area ratio of overexposed area and the area ratio of underexposed area in the current rendering frame is less than a first preset threshold, and outputs an illumination adaptation completion signal.

[0055] A real-time rendering engine is used to respond to the received lighting adaptation completion signal, to render the three-dimensional Gaussian splash model in real time at a preset rendering frame rate, and to collect pixel-level brightness distribution data for each rendering frame during the rendering process. The pixel-level brightness distribution data includes the brightness value of each pixel and an overexposure or underexposure flag determined based on a brightness threshold.

[0056] Specifically, the 3DGS rendering system also includes:

[0057] The optimization module is configured as follows:

[0058] Calculate the brightness ratio between the average brightness of the overexposed areas and the average brightness of the underexposed areas in the current rendering frame, and determine whether the brightness ratio is less than a preset second threshold.

[0059] If the brightness ratio is less than the second threshold, the current rendering frame is output based on the current spherical harmonic coefficients and the preset global tone mapping operator.

[0060] If the brightness ratio is greater than or equal to the second threshold, then the exposure fusion weight adjustment and tone mapping curve correction are calculated based on the pixel-level brightness distribution data. The exposure fusion weight adjustment is applied to the spherical harmonic coefficients of the Gaussian sphere to update the spherical harmonic coefficients, and the tone mapping curve correction is applied to the global tone mapping operator to update the global tone mapping operator. Then, the step of performing real-time rendering of the three-dimensional Gaussian splash model at a preset rendering frame rate is returned until the brightness ratio is less than the second threshold.

[0061] Compared with the prior art, the beneficial effects of the present invention are:

[0062] This invention constructs a technical path for reconstructing high dynamic range (HDR) lighting fields from multi-exposure image sequences. It obtains lighting field configuration data including dynamic range compression curves, indoor / outdoor region segmentation masks, and material reflectivity models, thus solving the rendering starting point deviation problem caused by missing lighting information in traditional methods with accurate prior information. Based on this, it employs an iterative optimization strategy for spherical harmonic coefficients guided by the lighting field, ensuring that the spherical harmonic coefficients of the Gaussian sphere accurately carry HDR lighting information. Furthermore, it combines indoor / outdoor region segmentation masks to classify the Gaussian sphere and allocate differentiated step sizes, resolving the convergence speed imbalance problem caused by differences in lighting characteristics between indoor low-frequency regions and outdoor high-frequency regions. Finally, it introduces overexposed and underexposed regions. The joint constraint mechanism of area ratio and multi-view loop optimization overcome the problems of insufficient viewpoint generalization ability and difficulty in controlling overexposure and underexposure coupling caused by single-viewpoint optimization. In the real-time rendering stage, a closed-loop feedback mechanism based on brightness ratio is established to dynamically calculate the exposure fusion weight adjustment amount and tone mapping curve correction amount and apply them to the spherical harmonic coefficient and global tone mapping operator respectively, effectively solving the temporal coupling problem between dynamic range compression curve and spherical harmonic coefficient optimization as well as the inter-frame brightness jitter problem. This embodiment achieves a visual consistency effect of clear details in bright areas, visible textures in dark areas, and natural transition of indoor and outdoor light and shadow in HDR large-scene virtual roaming, significantly improving the rendering realism and user experience. Attached Figure Description

[0063] Figure 1 This is a flowchart of the 3DGS rendering method for large HDR scenes according to Embodiment 1 of the present invention;

[0064] Figure 2 This is the logic diagram for Gaussian sphere partitioning in Embodiment 1 of the present invention;

[0065] Figure 3 This is a block diagram of the 3DGS rendering system for HDR large scenes in Embodiment 2 of the present invention. Detailed Implementation

[0066] Example 1

[0067] Please see Figure 1The present invention provides an embodiment of a 3DGS rendering method for large HDR scenes, comprising the following steps:

[0068] S1. Obtain the multi-exposure image sequence of the target scene and the corresponding camera pose parameters; In this embodiment, the camera pose parameters include a camera intrinsic parameter matrix and a camera extrinsic parameter matrix. The camera intrinsic parameter matrix includes focal length, principal point coordinates and distortion coefficients. The camera extrinsic parameter matrix includes a rotation matrix and a translation vector. The rotation matrix describes the rotation relationship between the camera coordinate system and the world coordinate system. The translation vector describes the spatial position of the camera optical center in the world coordinate system.

[0069] S2. Based on the multi-exposure image sequence and the camera pose parameters, reconstruct the high dynamic range illumination field to obtain illumination field configuration data. The illumination field configuration data includes at least the dynamic range compression curve, indoor and outdoor region segmentation mask, and material reflectivity model.

[0070] S3. Based on the lighting field configuration data, perform iterative optimization on the spherical harmonic coefficients of each Gaussian sphere in the three-dimensional Gaussian splash model. If the difference between the area ratio of overexposed region and the area ratio of underexposed region in the current rendering frame is less than a first preset threshold, output a lighting adaptation completion signal. If the difference is greater than or equal to the first preset threshold, continue to perform iterative optimization steps, including calculating the gradient of the spherical harmonic coefficients based on the rendering results of the current view direction and updating the spherical harmonic coefficients, and then reselecting the next view direction from multiple preset view directions for rendering and difference calculation, until the difference is less than the first preset threshold in multiple consecutive view directions, stop iterative optimization and output a lighting adaptation completion signal.

[0071] S4. In response to receiving the lighting adaptation completion signal, the three-dimensional Gaussian splash model is rendered in real time at a preset rendering frame rate, and pixel-level brightness distribution data of each rendering frame is collected during the rendering process. The pixel-level brightness distribution data includes the brightness value of each pixel and an overexposure or underexposure indicator determined based on the brightness threshold.

[0072] S5. Calculate the brightness ratio between the average brightness of the overexposed area and the average brightness of the underexposed area in the current rendering frame, and determine whether the brightness ratio is less than a preset second threshold.

[0073] S6. If the brightness ratio is less than the second threshold, the current rendering frame is output based on the current spherical harmonic coefficients and the preset global tone mapping operator.

[0074] S7. If the brightness ratio is greater than or equal to the second threshold, then calculate the exposure fusion weight adjustment and tone mapping curve correction based on the pixel-level brightness distribution data. Apply the exposure fusion weight adjustment to the spherical harmonic coefficients of the Gaussian sphere to update the spherical harmonic coefficients, apply the tone mapping curve correction to the global tone mapping operator to update the global tone mapping operator, and return to the step of real-time rendering of the 3D Gaussian splash model at a preset rendering frame rate until the brightness ratio is less than the second threshold. It should be further noted that this embodiment combines the 3D Gaussian splash model with a specific implementation in the context of a digital farmland scenario in a smart agriculture demonstration zone:

[0075] The first step involves selecting a 1 / 50-second standard exposure frame from the multi-exposure image sequence as a reference image and extracting scale-invariant feature transformation (SMT) feature points. A contrast threshold of 0.04 is set, and nearest neighbor distance ratio matching is performed on adjacent exposure images with a matching threshold of 0.8. A random sampling consensus algorithm is used to eliminate mismatched point pairs, with 500 iterations and a projection error threshold of 1.5 pixels. Using the selected matching point pairs, bundle adjustment is performed to jointly optimize the camera pose parameters, including the camera intrinsic matrix, extrinsic matrix, and the 3D coordinates of the sparse feature points, reconstructing a sparse 3D point cloud of the target scene containing 2 million feature points.

[0076] For the sparse 3D point cloud, a patch-based multi-view stereo matching algorithm is used for densification reconstruction to generate a dense 3D point cloud. For each 3D point in the dense point cloud, a set of points within its 16-neighborhood is taken, and the covariance matrix of the neighborhood point set is calculated using principal component analysis. The eigenvector corresponding to the smallest eigenvalue is taken as the normal vector of that point, thus assigning normal vector information to each dense point. Simultaneously, the RGB values ​​of the corresponding pixels in the 1 / 50-second standard exposure frame are assigned to the dense point as initial color information, ultimately obtaining a scene dense 3D point cloud containing 500 million points, carrying normal vectors and initial colors, providing a geometric reference for the Gaussian sphere center point and pose initialization.

[0077] The second step involves traversing each 3D point in the dense 3D point cloud and using its spatial coordinates as the center coordinates of the Gaussian sphere to be generated. An initial covariance matrix is ​​constructed using the normal vector of this 3D point as the axis. This covariance matrix is ​​a diagonal matrix, with its diagonal elements consisting of the squares of the lengths of the three semi-axises. Here, the lengths of the three semi-axises are initialized to 0.01 meters, thus defining the initial spatial extension range of the Gaussian sphere in space. The opacity parameter of each Gaussian sphere is initialized to 1.0. The 0th order spherical harmonic coefficients of each Gaussian sphere are initialized by dividing the RGB three-channel color values ​​of the corresponding 3D point of the Gaussian sphere by the constant value of the 0th order real spherical harmonic basis function, which is 1 / The value is 0.282095; thus, three 0th-order spherical harmonic coefficients are obtained, corresponding to the three RGB color channels respectively. Through the above operations, a total of 500 million initial Gaussian spheres are generated. Each initial Gaussian sphere contains core optimizable parameters such as center point coordinates, covariance matrix, opacity, and 0th-order spherical harmonic coefficients.

[0078] The third step involves constructing a differentiable loss function using the pixel color values ​​of the 1 / 50-second standard exposure frame and the depth map generated by the multi-view stereo matching algorithm in the first step as supervised ground truth. This loss function is a weighted sum of the color reprojection L1 error term and the depth reprojection L1 error term, with the color error term having a weight of 0.8 and the depth error term having a weight of 0.2. Iterative optimization is performed using the Adam optimizer, with hyperparameters β1=0.9 and β2=0.999, an initial iteration step size of 0.005, and a total of 30,000 iterations. In each iteration, based on the current center point coordinates, covariance matrix, opacity, and 0th-order spherical harmonic coefficients of each Gaussian sphere, a rendered image of the current viewpoint is generated using the Gaussian splashing algorithm. The loss function is then calculated, and all parameters to be optimized are updated via backpropagation.

[0079] During optimization, an adaptive splitting and pruning operation of the Gaussian sphere is performed every 1000 iterations. For any Gaussian sphere, if the pixel rendering reprojection error is greater than a preset error threshold of 0.01, the Gaussian sphere is split: the eigenvalues ​​of the original covariance matrix are scaled by a factor of 0.5, resulting in two new covariance matrices for the child Gaussian spheres; the center point coordinates of the two child Gaussian spheres are offset by 0.5 times the original semi-axis length (0.01 meters in this case) along the positive and negative directions of the Gaussian sphere's normal vector, respectively, i.e., 0.005 meters; the child Gaussian spheres inherit the opacity parameter and spherical harmonic coefficients of the original Gaussian sphere. For any Gaussian sphere, if the opacity is lower than the preset opacity threshold of 0.1, it is directly removed from the model. After 30,000 iterations and adaptive density control every 1000 iterations, an initial 3D Gaussian splash model composed of 50 million effective Gaussian spheres is obtained. This model has accurately reconstructed the geometric structure of the farmland scene and has basic color expression capabilities.

[0080] The fourth step involves updating the parameters of each Gaussian sphere in the initial 3D Gaussian splash model using the illumination field configuration data, ensuring that its illumination response characteristics are consistent with the reconstructed high dynamic range illumination field. The specific update operations are as follows:

[0081] First, a pre-classification of Gaussian sphere regions is performed. The center point coordinates of each Gaussian sphere and its spatial extent (ellipsoidal region) determined by the eigenvalues ​​of its covariance matrix are obtained. Based on the indoor / outdoor region segmentation mask (i.e., crop-soil segmentation mask) in the illumination field configuration data, a voxel intersection query is performed. The number of indoor voxels (crop voxels) and outdoor voxels (soil voxels) marked as such in the voxel set covered by the Gaussian sphere are counted, and the proportions of indoor and outdoor voxels are calculated respectively. When the proportion of indoor voxels is greater than or equal to a first classification threshold, the Gaussian sphere is classified into the indoor region Gaussian sphere set; when the proportion of outdoor voxels is greater than or equal to a second classification threshold, it is classified into the outdoor region Gaussian sphere set; otherwise, it is marked as a boundary Gaussian sphere. Both the first and second classification thresholds are preset values ​​between 0.5 and 0.8.

[0082] Secondly, for each Gaussian sphere, its spatial location is determined by the coordinates of its center point in the voxel grid of the illumination field radiation distribution. A trilinear interpolation algorithm is used to calculate the radiance value corresponding to the location by weighted averaging the radiance values ​​of the eight corner points of the voxel where the center point is located, and this value is then updated as the reference radiance of the Gaussian sphere.

[0083] Then, based on the region type to which the Gaussian sphere belongs, the corresponding material parameters are obtained from the material reflectivity model: for Gaussian spheres belonging to the indoor region Gaussian sphere set, the diffuse reflection coefficient corresponding to the voxel is obtained; for Gaussian spheres belonging to the outdoor region Gaussian sphere set, the roughness and metallicity corresponding to the voxel are obtained; for boundary Gaussian spheres, the diffuse reflection coefficient, roughness, and metallicity are obtained simultaneously, and then weighted and fused according to the proportion of indoor and outdoor voxels in subsequent gradient calculations.

[0084] Finally, full-order spherical harmonic coefficient initialization is performed to replace the initial 0th-order spherical harmonic coefficients. Specifically, for each Gaussian sphere, 64 sampling directions are uniformly generated on the unit sphere surface using the Fibonacci grid sampling method. For each sampling direction, the theoretical radiance value is calculated using the corresponding material reflection model based on the region type of the Gaussian sphere: for indoor Gaussian spheres, the updated reference radiance is used as the incident light irradiance, and the Lambertian diffuse reflection model is used to calculate the radiance by combining the diffuse reflection coefficient and the dot product of the normal vector; for outdoor Gaussian spheres, the Cook-Torrance micro-plane specular reflection model is used to calculate the radiance by combining the incident light direction, sampling direction, normal vector, roughness, and metallicity. After traversing the 64 sampling directions, an observation vector containing 64 theoretical radiance values ​​is constructed. The highest spherical harmonic order is 5 (corresponding to (5+1)). 2=36 spherical harmonic coefficients), and the coefficient matrix is ​​solved by fitting the least squares method to minimize the sum of the squared errors between the output value of the spherical harmonic function in each sampling direction and the corresponding theoretical radiance value in the observation vector. The 36 coefficients obtained are the full-order spherical harmonic coefficients initialized for the Gaussian sphere.

[0085] After completing all the above update operations, a three-dimensional Gaussian splash model is obtained in which each Gaussian sphere has the same reference radiance, material parameters and full-order spherical harmonic coefficients as the illumination field. This model can be directly used as the initial state for the subsequent illumination field-guided iterative optimization steps of the spherical harmonic coefficients.

[0086] It should be further explained that one method for reconstructing the high dynamic range illumination field in this embodiment is as follows:

[0087] S201. Extract scale-invariant feature transformation points from each frame of the multi-exposure image sequence. Perform nearest neighbor distance ratio matching on image frames with adjacent exposure parameters, and use a random sampling consensus algorithm to remove mismatched point pairs. Calculate the inter-frame homography transformation matrix based on the filtered matching point pairs. Perform perspective transformation and resampling on the images according to the homography transformation matrix, so that the same spatial point in different exposure images is mapped to the same pixel coordinates, resulting in an aligned multi-exposure image sequence. Next, a specific complete example illustrates the entire process. The example is only to illustrate the feasibility at the computational level and does not represent actual values. The specific values ​​can be determined by those skilled in the art through simulation experiments or physical experiments. For example, it specifically includes:

[0088] The global average brightness value is calculated for each frame: 40 for a 1 / 200 second image (underexposed), 120 for a 1 / 50 second image (normal), 220 for a 1 / 10 second image (overexposed), and 250 for a 0.2 second image (overexposed). For underexposed images (1 / 200 second) with an average brightness less than 50, the contrast threshold in the scale-invariant feature transform algorithm is increased from the default 0.04 to 0.06 to suppress soil noise in dark areas; for overexposed images (1 / 10 second and 0.2 second) with an average brightness greater than 200, the contrast threshold is decreased to 0.02 to preserve the highlight edge features of corn leaves (such as vein reflections).

[0089] Calculate the brightness ratio of adjacent exposure images (the ratio of the average brightness of the high-exposure image to the average brightness of the low-exposure image): The brightness ratio of 1 / 200 second to 1 / 50 second is 120 / 40=3.0, which is between 2 and 4, so the threshold for nearest neighbor distance ratio matching is set to 0.8; the brightness ratio of 1 / 50 second to 1 / 10 second is 220 / 120≈1.83, which is less than 2, so the threshold is tightened to 0.7; the brightness ratio of 1 / 10 second to 0.2 second is 250 / 220≈1.14, which is also less than 2, so the threshold remains at 0.7. In this embodiment, the brightness ratio between adjacent exposed images reflects the degree of exposure difference and feature sharpness between the images. When the brightness ratio is between 2 and 4, it indicates that the brightness difference between the two images is large but not excessive. At this time, the overexposed or underexposed areas have not completely lost texture information, the contrast of feature points is still high, and the matching reliability is strong. Therefore, a higher threshold of 0.8 is used to retain more matching point pairs and enhance the stability of subsequent homography matrix estimation. When the brightness ratio is less than 2, the image brightness is similar and the feature discrimination decreases. The threshold needs to be tightened to 0.7 to ensure the matching quality. When the brightness ratio is greater than 4, the overexposed or underexposed areas are too large, causing feature point degradation. In this case, the number of iterations needs to be increased or a weighted strategy needs to be adopted to eliminate mismatches. This dynamic adjustment mechanism quantifies the illumination difference between images through the brightness ratio, thereby adaptively controlling the matching threshold to ensure that the feature matching of different exposure image sequences in high dynamic range large scenes such as farmland and ancient buildings is both sufficient and robust.

[0090] A progressive random sampling consensus algorithm is used to eliminate mismatches. For example, the initial homography matrix is ​​quickly estimated with 500 iterations, and the proportion of inliers with projection errors less than 1.5 pixels is counted. The proportion of inliers at 1 / 200 second and 1 / 50 second is 78%, and this is maintained for 500 iterations. The proportion of inliers at 1 / 50 second and 1 / 10 second is 52% (below 60%), and the number of iterations is increased to 1000. The proportion of inliers at 1 / 10 second and 0.2 second is 35% (below 40%), and the number of iterations is increased to 2000, until the inlier proportions stabilize.

[0091] For the finally selected correct matching point pairs, a weighted direct linear transformation is used to recalculate the inter-frame homography matrix. The weights are set according to the Euclidean distance of the feature points from the image center, with greater weights for points closer to the center (e.g., center point weight 1.0, edge point weight 0.6) to reduce the impact of perspective distortion. The calculated homography matrix parameters for 1 / 200 second and 1 / 50 second are: first row [0.999, -0.015, 1.2], second row [0.014, 0.999, 0.8], third row [0.00005, -0.0001, 1].

[0092] Using a 1 / 50-second exposure image as a baseline, the other three images are transformed using perspective based on their corresponding homography matrices. Bicubic interpolation is then used to resample pixels, replacing conventional bilinear interpolation, to preserve high-frequency details (such as corn leaf texture). Ultimately, the alignment error of feature points such as leaf tips and veins of the same corn plant across different exposure images is less than 0.3 pixels.

[0093] S202. For each pixel position in the aligned multi-exposure image sequence, obtain the original brightness value of that pixel in each frame image. Convert the original brightness value into a linear radiance value according to a pre-calibrated camera response function. Use a preset weighting function to weight and fuse all linear radiance values ​​at the same pixel position to generate a high dynamic range radiance value for that pixel position. After traversing all pixel positions, a high dynamic range radiance image is obtained. Next, a specific complete example will be used to illustrate the entire process. The example is only to illustrate the feasibility at the computational level and does not represent the actual values. The specific values ​​can be determined by those skilled in the art through simulation experiments or physical experiments. For example, it specifically includes:

[0094] The camera response function was calibrated using 20 neutral gray card images. A fifth-order polynomial was used to fit and obtain independent response curves for the red, green, and blue channels (e.g., the red channel response function is R1 = 0.82L0 + 0.09L0). 2 -0.018L0 3 (where L0 is the original brightness value) to eliminate color cast.

[0095] Taking a pixel coordinate (1500, 1200) of the top leaf of a corn plant as an example, the original brightness values ​​are obtained from the four aligned images: I1=40, I2=120, I3=220, I4=250. Based on the corresponding color channel response curves, these values ​​are converted to linear radiance values, resulting in E1=0.3, E2=2.2, E3=9.5, and E4=16.0 (unit: watts per steradian per square meter).

[0096] The confidence level for each raw brightness value was calculated using a Gaussian function with a mean of 128 and an adaptive standard deviation σ as the parameter. σ was determined based on the image's ISO sensitivity: σ was 30 for ISO 100 and 40 for ISO 200. In this example, all four frames had an ISO of 200, so σ = 40. The confidence levels were calculated as follows: w1 = 0.25, w2 = 0.92, w3 = 0.65, w4 = 0.15 (I4 = 250, close to overexposure, low confidence level).

[0097] Calculate the brightness variance within a 3×3 neighborhood of the pixel, with a preset variance threshold of 20. If the variance is less than 20 (smooth region), multiply the confidence score by 1.2; if the variance is between 20 and 80 (general texture), the confidence score remains unchanged; if the variance is greater than or equal to 80 (edge ​​or noise), multiply the confidence score by 0.8. In this example, the pixel is located in the middle of the leaf, with a neighborhood variance of 18, belonging to a smooth region. Therefore, the confidence scores are multiplied by 1.2 respectively, updating to w1=0.3, w2=1.10, w3=0.78, and w4=0.18. Simultaneously, since I4=250>250, it is determined to be an overexposed pixel, and its confidence score is set to zero.

[0098] The four radiance values ​​are sorted by confidence level and their cumulative weights are calculated: the sorted values ​​are E4 (confidence level 0), E1 (0.3), E3 (0.78), and E2 (1.10). E3 corresponds to a cumulative weight of 50%, so the initial radiance value is 9.5. Since this pixel is overexposed, bilinear interpolation is used to complete the measurement using the initial radiance values ​​of non-overexposed pixels in the neighborhood (e.g., neighborhood pixel radiance values ​​of 9.2, 9.6, and 9.3), resulting in a final radiance value of 9.45.

[0099] After traversing all pixels, a high dynamic range radiance image covering 100 acres of farmland is generated, with an image size of 8000×6000 pixels and a dynamic range from 0.2 to 50.0 watts per spheradeitude per square meter.

[0100] S203. Based on the camera intrinsic and extrinsic parameter matrices included in the camera pose parameters, calculate the corresponding normalized ray direction for each pixel in the high dynamic range radiosity image, transform the ray direction to the world coordinate system, and perform spatial point sampling along the ray direction at a preset voxel resolution. Assign the high dynamic range radiosity value of the pixel to the first voxel intersecting the pre-generated scene geometry surface along the ray direction. After performing the same operation on all images, eliminate assignment noise through multi-view consistency verification to generate a voxel grid of illumination field radiation distribution in three-dimensional space. Next, a specific complete example will be used to illustrate the entire process. The example is only to illustrate the feasibility at the computational level and does not represent the actual values. The specific values ​​can be determined by those skilled in the art through simulation experiments or physical experiments. For example, this embodiment takes a digital farmland in a smart agriculture demonstration area as an example. Density clustering processing is performed on the pre-generated farmland scene point cloud, which is acquired by UAV oblique photography and LiDAR, with a total number of approximately 500 million points. By calculating the average point density in the local neighborhood, the point cloud is divided into dense and sparse regions. For dense areas such as the corn plant canopy, the point density is greater than 1000 points per cubic meter, and the voxel resolution is set to 0.005 meters to preserve the detailed structure of plant leaves, ears, etc. For flat and sparse areas such as soil and field ridges, the point density is less than or equal to 1000 points per cubic meter, and the voxel resolution is set to 0.02 meters to control the overall number of voxels while ensuring reconstruction accuracy.

[0101] Taking the pixel corresponding to the base of a corn stalk as an example, the pixel coordinates are (2000, 3000), and its radiance value is 5.2. Based on the camera intrinsic parameter matrix, the normalized ray direction in the camera coordinate system is calculated, where the focal length fx is 3500 pixels, fy is 3500 pixels, the principal point cx is 2048 pixels, and cy is 1536 pixels. The calculated ray directions are 0.314, 0.285, and 1.0. Based on the extrinsic matrix of this frame, the ray direction is transformed to the world coordinate system. The elements of the rotation matrix R are R11=0.995, R12=-0.087, R13=0.052, R21=0.086, R22=0.996, R23=0.028, R31=-0.054, R32=-0.023, and R33=0.998. The translation vectors T0 are 50.0 meters, 30.0 meters, and 20.0 meters. After transformation, the direction vectors in the world coordinate system are 0.298, 0.271, and 0.915.

[0102] Equal-interval sampling is performed along the ray direction with a step size of half the current voxel resolution. For dense corn canopy areas, the step size is 0.0025 meters, and for flat soil areas, it is 0.01 meters. An octree data structure is used to query the Euclidean distance between each sampling point and the nearest point in the point cloud. When this distance is less than 0.005 meters, the ray is considered to intersect the scene's geometric surface. This distance threshold is set based on a combination of point cloud accuracy and reconstruction error to ensure the reliability of the intersection determination. The world coordinates of the intersection points are recorded as 52.3 meters, 32.5 meters, and 1.2 meters, and the corresponding voxel indices are 5230, 3250, and 120. Simultaneously, the angle between the ray direction and the surface normal is calculated. The calculated angle is 25°, less than the preset angle threshold of 60°. Therefore, the confidence weight of this observation is set to 1.0. This angle threshold is used to eliminate observations with low confidence at grazing angles to avoid introducing large errors.

[0103] The voxel is assigned a radiosity value of 5.2, and the spatial distance from the camera to the voxel is recorded as approximately 58.5 meters, along with the corresponding confidence weight. This process is repeated for each pixel in all four images, where each voxel may be observed from multiple viewpoints.

[0104] For all radiance observations collected for each voxel, a fusion weight is calculated for each observation. The fusion weight consists of the product of a confidence weight, a cosine attenuation factor, and a distance attenuation factor. The cosine attenuation factor is the cosine of the angle between the observation direction and the surface normal; in this example, the angle is 25°, and the cosine value is 0.906. The distance attenuation factor is a Gaussian function, specifically exp(-d²) divided by σ_d², where d represents the Euclidean distance from the spatial point corresponding to the pixel to the camera's optical center, used to quantify the impact of observation distance on the confidence of the radiance value, and σ_d represents the scale parameter of distance attenuation, set to 10.0 meters, used to control the rate of weight attenuation with increasing distance. This distance attenuation mechanism ensures that observations closer to the camera receive higher weights, while the weights of observations farther away decrease exponentially, effectively suppressing noise introduced by long-distance observations due to viewpoint deviation, atmospheric scattering, or insufficient resolution. In this example, when the distance is 58.5 meters, d squared is 3422 square meters and σ_d squared is 100 square meters. The calculated attenuation factor is approximately 0.00003. Due to its extremely low weight, this observation was discarded in actual fusion, indicating that effective observations usually come from multiple perspectives at relatively close distances.

[0105] For each voxel, the radiance values ​​and corresponding weights of all valid observations are collected. In this example, radiance values ​​of 5.1, 5.2, 5.3, and 5.0 and their corresponding fusion weights are collected. After removing outlier observations with fusion weights below 0.1, a weighted average is performed to obtain the final radiance value of 5.18 for this voxel. The weight removal threshold of 0.1 is used to exclude observations with excessively low confidence, ensuring the stability and accuracy of the voxel radiance values.

[0106] After traversing all images and voxels according to the above process, a three-dimensional light field radiation distribution voxel grid covering the farmland is finally generated, with a total effective voxel count of approximately 50 million. The voxel resolution in the dense corn canopy area is 0.005 meters, and the voxel resolution in the flat soil area is 0.02 meters, achieving a refined representation and efficient storage of light information in different areas.

[0107] S204. The voxel grid of the illumination field radiation distribution is converted into a multi-channel feature tensor. The multi-channel feature tensor contains the radiance value, normal vector, and spatial coordinates of each voxel. The multi-channel feature tensor is input into a pre-trained three-dimensional convolutional neural network to obtain the probability value of each voxel belonging to the indoor or outdoor region. Conditional random field post-processing is performed on the probability value to generate a binarized indoor / outdoor region segmentation mask. Next, a specific complete example is used to illustrate the entire process. The example is only to illustrate the feasibility at the computational level and does not represent the actual values. The specific values ​​can be determined by those skilled in the art through simulation experiments or physical experiments. For example, the pre-generated three-dimensional voxel grid is feature-converted to construct a seven-channel feature tensor. The seven channels are, in order, the radiance value, the three components of the normal vector, and the three components of the spatial coordinates. The feature tensor was input into a pre-trained 3D convolutional neural network (3DU-Net), which uses a 3D U-shaped network structure. The encoder consists of four convolutional layers, each with a kernel size of 3×3×3, and output channels of 32, 64, 128, and 256, respectively. The decoder consists of four deconvolutional layers. The network was trained on a labeled dataset containing 10 million voxels across various crop types. The dataset covers major crop varieties such as maize, wheat, and rice, along with their corresponding soil backgrounds. Voxel-level annotations were jointly performed manually using a combination of multispectral imagery and LiDAR data.

[0108] The network outputs the probability values ​​for each voxel belonging to the crop category and the probability values ​​for each voxel belonging to the soil category. Taking the corn leaf voxel as an example, the network outputs a crop probability of 0.95 and a soil probability of 0.05; taking the bare soil voxel as an example, the network outputs a crop probability of 0.02 and a soil probability of 0.98. The sum of the crop probability and the soil probability is 1, and both are normalized using the Softmax function.

[0109] Conditional random field (CRF) post-processing is applied to the probability tensor output by the network to optimize the edge consistency and spatial continuity of the segmentation results. A graph model with voxels as nodes is constructed, where neighboring voxels are defined as 26 neighborhoods in 3D space, meaning each voxel is connected to its 26 surrounding voxels. The connection weights between neighboring voxels are adaptively adjusted based on radiosity differences and spatial distances; voxels with smaller radiosity differences and closer spatial distances have larger connection weights. Specifically, a Gaussian kernel function is used to fuse these differences for calculation. The voxel class probabilities are updated using a CRF iterative optimization algorithm, with 20 iterations. This number of iterations is determined based on the convergence curve of segmentation accuracy versus iterations on the validation set. When the number of iterations exceeds 20, the improvement in segmentation accuracy is less than 0.5%, so 20 iterations are selected as the preset value to balance accuracy and efficiency.

[0110] After iterative optimization, a binary crop-soil segmentation mask is generated based on the optimized category probabilities. A probability threshold of 0.5 is set. For each voxel, if its crop probability is greater than or equal to 0.5, it is classified as a crop region; if the crop probability is less than 0.5, it is classified as a soil region. This probability threshold is determined based on the minimum error rate criterion in Bayesian decision theory to ensure that classification error is minimized. The generated segmentation mask data is approximately 500 million bits, and after run-length encoding compression, its size is approximately 0.6 GB, with a compression ratio of approximately 12:1, achieving efficient storage of large-scale voxel segmentation results.

[0111] S205. For all voxels marked as indoor areas in the indoor / outdoor area segmentation mask, the diffuse reflection coefficient of each voxel is solved by least squares fitting based on the Lambertian diffuse reflection model and the incident light direction and radiance value at the voxel location. For all voxels marked as outdoor areas, the specular reflection coefficient of each voxel is solved by nonlinear optimization based on the microplane theory specular reflection model and the voxel normal, observation direction, and radiance value. The diffuse reflection coefficient and the specular reflection coefficient are combined into a material reflectivity model. Next, a specific complete example is used to illustrate the entire process. The example is only to illustrate the feasibility at the computational level and does not represent the actual values. The specific values ​​can be determined by simulation experiments or physical experiments conducted by those skilled in the art. For example, it specifically includes:

[0112] For the crop area, a maize leaf voxel was used as an example for material parameter fitting. The spatial coordinates of this leaf voxel are X = 52.3 m, Y = 32.5 m, Z = 1.2 m, with normal vectors of 0.12, 0.15, and 0.98, incident light directions of 0.65, 0.48, and 0.59, and a radiance value of 5.2. Material parameter optimization was carried out in two stages. The first stage was based on the Lambertian diffuse reflection model, assuming the leaf surface is approximately an ideal diffuse reflector. The initial diffuse reflection coefficient was fitted using the radiance observations of this voxel at eight different viewing angles, and the fitting result was 0.68. This diffuse reflection coefficient describes the leaf's ability to diffusely reflect incident light, and its value ranges between 0 and 1, with a larger value indicating stronger diffuse reflection. The second stage introduced a specular reflection term based on the fixed diffuse reflection coefficient to characterize the weak specular reflection characteristics of the leaf surface. Through optimization, the roughness parameter was obtained as 0.45 and the metallicity parameter as 0.02. The roughness parameter describes the smoothness of the surface microstructure; the smaller the value, the smoother the surface and the more concentrated the specular reflection. Its value ranges from 0 to 1. The metallicity parameter describes the metallic properties of the material; the closer the value is to 1, the stronger the metallic texture. Its value also ranges from 0 to 1. Since the actual specular component of a leaf is very weak, a low metallicity value is chosen, consistent with the physical and optical characteristics of plant leaves.

[0113] For the soil region, material parameters were optimized using a bare soil voxel as an example. The spatial coordinates of this soil voxel are X = 53.1 m, Y = 33.2 m, and Z = 0.8 m, with normal vectors of 0.05, 0.10, and 0.99, observation directions of 0.58, 0.62, and 0.53, incident light directions of 0.65, 0.48, and 0.59, and a radiance value of 3.5. Radiance observations of this voxel were collected from eight different viewpoints, yielding values ​​of 3.5, 3.4, 3.6, 3.5, 3.3, 3.7, 3.4, and 3.6. The objective function was constructed as the sum of squares of the radiance residuals from each viewpoint plus a neighborhood regularization term. The weight coefficient of the neighborhood regularization term was set to 0.1 to constrain the smoothness of material parameters of adjacent voxels, avoiding drastic parameter fluctuations due to observation noise. The Levenberg-Marquardt nonlinear optimization algorithm was used for the solution, with initial roughness values ​​of 0.4 and metallicity values ​​of 0.1, and an upper limit of 30 iterations. After convergence, the roughness was 0.38 and the metallicity was 0.12. An early stopping mechanism was implemented during the optimization process. If the residual decrease was less than 0.01 for five consecutive iterations, the optimization was considered to have entered a local minimum or stalled. In this case, a random restart strategy was executed, randomly perturbing the current parameters within a range of ±10% before re-optimizing to escape the local minimum region.

[0114] After optimizing the material parameters of all voxels, a bilateral filtering post-processing is performed on the generated material parameter field to achieve a smooth transition between adjacent voxels while maintaining sharp edges. The bilateral filtering has a standard deviation of 0.1 meters in the spatial domain. This parameter is set based on a combination of voxel mesh resolution and scene scale to ensure that the filtering window covers a sufficient number of adjacent voxels to achieve effective smoothing. In terms of value range, the roughness parameter has a standard deviation of 0.05, and the metallicity parameter has a standard deviation of 0.1. These two standard deviations are determined based on the natural variation range of the material parameters, ensuring that voxels with parameter differences less than the standard deviation smooth each other during filtering, while voxels with differences greater than the standard deviation maintain sharp edges, thus achieving boundary preservation at the junction of crops and soil. Finally, a material reflectivity model covering 50 million voxels is generated, where the crop region stores the diffuse reflectance coefficient, and the soil region stores the specular reflectance coefficient, i.e., the roughness and metallicity parameters, providing physically accurate material property support for subsequent lighting rendering.

[0115] S206. Statistically analyze the radiance values ​​of all voxels in the voxel grid of the illumination field radiation distribution, determine the minimum and maximum radiance values, obtain the maximum and minimum displayable brightness of the target display device, design a smooth mapping function to map the minimum radiance value to the minimum displayable brightness and the maximum radiance value to the maximum displayable brightness, and save the parameters of the mapping function as a dynamic range compression curve. Next, a specific complete example will illustrate the entire process. This example is only to illustrate the feasibility at the computational level and does not represent actual values. Specific values ​​can be determined by those skilled in the art through simulation experiments or physical experiments. For example, it specifically includes:

[0116] The radiance values ​​of all voxels in the voxel grid of the illumination field radiation distribution are statistically analyzed, and a one-dimensional frequency distribution histogram of the radiance values ​​is generated. The horizontal axis of the histogram is the range of radiance values, and the vertical axis is the number of voxels falling into each range.

[0117] Based on the frequency distribution histogram, the cumulative distribution function of the radiance value is determined, and the first radiance value corresponding to the preset first percentile and the second radiance value corresponding to the preset second percentile in the cumulative distribution function are calculated. The first percentile is a preset value between 1 and 5, and the second percentile is a preset value between 95 and 99.

[0118] Using the first radiance value as the lower limit of the mapping and the second radiance value as the upper limit of the mapping, the maximum and minimum displayable brightness of the target display device are obtained.

[0119] The radiance value range is divided into three segments: dark area, middle area, and bright area. The dark area is the part where the radiance value is less than the first radiance value, the middle area is the part where the radiance value is between the first radiance value and the second radiance value, and the bright area is the part where the radiance value is greater than the second radiance value.

[0120] A first mapping function is designed for the dark area. The first mapping function is a linear function that maps the minimum radiance value to the minimum displayable brightness and maps the first radiance value to a first intermediate brightness. The first intermediate brightness is a preset percentage between 10% and 20% of the maximum displayable brightness.

[0121] A second mapping function is designed for the intermediate range. This second mapping function is a logarithmic mapping function, used to map the clipped radiance values ​​between the first and second radiance values ​​to the target display brightness range. The complete expression of the logarithmic mapping function is as follows:

[0122] ,

[0123] Among them, L d L is the output display brightness value after mapping the current pixel. mind L is the minimum displayable brightness of the target display device. maxd L represents the maximum displayable brightness of the target display device. rad L is the radiosity value of the current pixel after cropping. rad_low L is the first radiance value (the lower limit of the mapping in the middle interval). rad_high Let α be the second radiance value (the upper limit of the mapping in the middle interval), and let α be the contrast adjustment factor. This contrast adjustment factor is embedded in the independent variable of the logarithmic function as a linear multiplier to control the slope of the middle segment of the logarithmic mapping curve. A larger α value results in a steeper slope in the middle segment and a stronger contrast in the mapped image; a smaller α value results in a shallower slope in the middle segment and a smoother contrast in the mapped image. The contrast adjustment factor α is adaptively calculated based on the standard deviation of the radiance values ​​within the middle interval, and its adaptive calculation function is: , where σ rad σ is the standard deviation of all radiance values ​​within the intermediate interval. globalThis represents the global standard deviation of the radiance values ​​of all voxels in the voxel grid of the target scene's illumination field radiation distribution. The above function implements an adaptive adjustment logic where the larger the standard deviation of the radiance values ​​in the middle interval, the larger the contrast adjustment factor α becomes. At the same time, through the fixed coefficient constraint of the function, the calculated result of α is always within the preset range of 1.0 to 10.0. If the calculated α exceeds this range, it is clamped to the corresponding range boundary value to ensure the calculation stability of the mapping function and the consistency of the display effect.

[0124] The larger the standard deviation, the larger the contrast adjustment factor; wherein the contrast adjustment factor is a dimensionless positive number, and its value ranges from 1.0 to 10.0, and is used to control the middle slope of the logarithmic mapping function in the dynamic range compression curve;

[0125] A third mapping function is designed for the bright area. The third mapping function is a linear function or a compression function. The second radiance value is mapped to the second intermediate brightness, and the maximum radiance value is mapped to the maximum displayable brightness. The second intermediate brightness is a preset percentage between 80% and 90% of the maximum displayable brightness.

[0126] Obtain the brightness response curve of the target display device. The brightness response curve describes the nonlinear relationship between the input brightness value and the actual display brightness. Based on the brightness response curve, pre-correct the first mapping function, the second mapping function, and the third mapping function so that the corrected mapping function restores the designed brightness relationship when outputting on the display device.

[0127] The parameters of the first, second, and third mapping functions, along with the first and second radiometric values, are saved together as the configuration parameters for the dynamic range compression curve.

[0128] S207. Output the indoor and outdoor area segmentation mask, the material reflectivity model, and the dynamic range compression curve as the illumination field configuration data.

[0129] It should be further noted that, before performing iterative optimization on the spherical harmonic coefficients of each Gaussian sphere in the three-dimensional Gaussian splash model, this embodiment also includes a spherical harmonic coefficient initialization step, including:

[0130] For each Gaussian sphere in the 3D Gaussian splash model, obtain the 3D spatial coordinates corresponding to its center position;

[0131] The radiance value at this coordinate is read from the voxel grid of the illumination field radiation distribution and used as the reference radiance of the Gaussian sphere;

[0132] The type of the Gaussian sphere is determined based on the indoor and outdoor area segmentation mask. If it is an indoor area, the corresponding diffuse reflection coefficient is obtained from the material reflectivity model. If it is an outdoor area, the roughness and metallicity in the corresponding specular reflection coefficient are obtained.

[0133] Based on the reference emissivity and material parameters of the Gaussian sphere, a spherical sampling method is used to uniformly generate K sampling directions on a unit spherical surface, where K is a preset integer greater than or equal to 32;

[0134] For each sampling direction, the theoretical radiance value under that direction is calculated by combining the incident light direction and the material reflection model, and a radiance value vector is constructed.

[0135] Using the radiance value vector as the observed value, the spherical harmonic coefficients are fitted by the least squares method to solve for the coefficient matrix that minimizes the error between the output of the spherical harmonic function in each sampling direction and the theoretical radiance value.

[0136] The obtained spherical harmonic coefficients are used as the initial spherical harmonic coefficients of the Gaussian sphere for subsequent iterative optimization.

[0137] Next, a specific complete example will be used to illustrate the whole process. The example is only to illustrate the feasibility of the calculation and does not represent the actual values. The specific values ​​can be determined by those skilled in the art through simulation experiments or physical experiments. For example, this embodiment takes the digital farmland scene of a smart agriculture demonstration area as an example to initialize the spherical harmonic coefficients of the Gaussian sphere in the three-dimensional Gaussian splash model corresponding to a certain leaf at the top of the corn plant canopy.

[0138] First, the three-dimensional spatial coordinates of the center of the Gaussian sphere were obtained as 12.5 meters, 8.3 meters, and 2.1 meters. The radiance value at these coordinates was read from the voxel grid of the three-dimensional illumination field radiation distribution covering 100 acres of farmland as the reference radiance, which was measured to be 8.5 watts per steradian per square meter. If there was no corresponding radiance value in the voxel grid at these coordinates, it was calculated from neighboring voxels using trilinear interpolation to ensure the continuity of the radiance value.

[0139] Based on the crop soil segmentation mask, the region to which the Gaussian sphere belongs is determined to be a crop region. The corresponding leaf material parameters are obtained from the material reflectance model, including a roughness of 0.25 and a metallicity of 0.05. The roughness parameter describes the smoothness of the leaf surface microstructure, while the metallicity parameter characterizes the strength of the leaf's metallic texture. Both values ​​range from 0 to 1. Crop leaves primarily exhibit diffuse reflection, and the low metallicity value aligns with their physical and optical characteristics.

[0140] Based on the reference radiometry of 8.5 and the material parameters of the Gaussian sphere, 64 sampling directions were uniformly generated on the unit sphere surface using the Fibonacci grid sampling method. The selection of the number of sampling directions was determined based on the order of the spherical harmonic function and the sampling theorem. When the order of the spherical harmonic function was set to 5, the minimum number of sampling points required was twice the square of the order, i.e., 50. The actual selection of 64 sampling directions satisfied the requirements of the Nyquist sampling theorem, ensuring the numerical stability of the spherical harmonic coefficient fitting.

[0141] For each sampling direction, the theoretical radiance value in that direction is calculated by combining the solar incidence direction with the material reflection function based on the Cook-Torrance model. The solar incidence direction is obtained during the illumination field reconstruction stage; in this example, a zenith angle of 30 degrees and an azimuth angle of 120 degrees are used, representing the incidence angle of the main light source in the scene. The Cook-Torrance model comprehensively considers diffuse and specular reflection components. The diffuse reflection component is determined by the reference radiance and roughness parameters, while the specular reflection component is calculated by the metallicity parameter and the geometric relationship between the incidence and observation directions. By traversing 64 sampling directions, an observation vector containing 64 theoretical radiance values ​​is constructed.

[0142] Using the theoretical radiance vector as the observed value, the spherical harmonic coefficients are fitted using the least squares method. The spherical harmonic function is set to order 5, corresponding to 36 coefficients to be solved. This order is determined based on the complexity of the illumination field's radiation distribution; a 5th-order spherical harmonic function can accurately represent the illumination variations caused by a mixture of diffuse and specular reflections, while avoiding excessive computational overhead introduced by higher-order coefficients. The objective function of the least squares method is to minimize the sum of the squares of the differences between the spherical harmonic function output value and the theoretical radiance value in each sampling direction. This results in a set of 36 5th-order spherical harmonic coefficients, which serve as the initial spherical harmonic coefficients for the Gaussian sphere.

[0143] Through the initialization process, the spherical harmonic coefficients of each Gaussian sphere already carry the true radiometric information and material physical properties of that location before the iterative optimization begins. This ensures that the corn leaves can accurately reproduce their true reflective characteristics under different light directions, providing an initial state close to the ideal solution for subsequent iterative optimization, effectively reducing the number of iterations and avoiding getting trapped in local optima.

[0144] This embodiment introduces an initialization step based on the radiation distribution of the illumination field and the material reflectivity model before iterative optimization of the spherical harmonic coefficients. By obtaining the reference radiance value corresponding to the center coordinates of each Gaussian sphere, and combining this with the region type determined by the indoor / outdoor region segmentation mask, the roughness and metallicity in the corresponding diffuse reflection coefficient or specular reflection coefficient are obtained. This allows the initial spherical harmonic coefficients to directly carry the real radiance information and material physical properties of the scene. A spherical sampling method is used to uniformly generate multiple sampling directions on a unit sphere. The theoretical radiance values ​​in each direction are calculated by combining the incident illumination direction and the corresponding material reflection model, and a radiance value vector is constructed. Then, the least squares method is used to fit and solve the sphere... The coefficient matrix with the smallest error between the harmonic function output and the theoretical radiance value ensures that the initial spherical harmonic coefficients can accurately represent the illumination response characteristics of the Gaussian sphere in multiple viewing directions. This initialization method fundamentally solves the problems caused by conventional random initialization or simple assignment, where the high-frequency components of the spherical harmonic coefficients are insufficient to support the detail restoration of bright areas and the low-frequency components are too weak, resulting in the loss of texture in dark areas. It makes the starting point of iterative optimization closer to the ideal state, effectively reducing the number of subsequent iterations and avoiding getting trapped in local optima. At the same time, it provides accurate initial values ​​for the step size and gradient calculation of the differential optimization of indoor and outdoor areas, ultimately improving the convergence efficiency of the spherical harmonic coefficient iterative optimization and the visual realism of HDR large scene rendering.

[0145] In the iterative optimization of spherical harmonic coefficients for a 3D Gaussian splash model for large HDR scenes, the radiance values ​​of the high dynamic range illumination field can span two orders of magnitude. Conventional initialization methods result in insufficient high-frequency components of the spherical harmonic coefficients to support detail restoration in bright areas, while weak low-frequency components lead to loss of texture in dark areas, causing the optimization starting point to deviate from the ideal state. Simultaneously, there is a coupling relationship between the area ratio control of overexposed and underexposed regions; enhancing high-frequency components may trigger new overexposure, while increasing low-frequency components may lead to new underexposure. The lack of a joint constraint mechanism makes it difficult for both to converge synchronously. Furthermore, optimization based solely on a single viewpoint easily falls into local optima, causing the optimized spherical harmonic coefficients to exhibit drastic changes in overexposure or underexposure when the viewpoint changes, resulting in insufficient viewpoint generalization ability. In large scenes, the lighting characteristics of indoor and outdoor areas differ significantly; therefore, a unified iteration step size and optimization... The strategy results in slow convergence in indoor areas and oscillation and divergence in outdoor areas; a temporal coupling occurs between the dynamic range compression curve and the spherical harmonic coefficient optimization. Adjusting the spherical harmonic coefficient changes the rendering brightness distribution, thus affecting the applicability of the compression curve. Changes in the compression curve parameters, in turn, change the overexposure and underexposure judgment results. The lack of a collaborative optimization mechanism in the existing technology leads to repeated adjustments; moreover, the use of a fixed value for the first preset threshold cannot adapt to the brightness distribution characteristics of different scenes. If the threshold is too loose, the overexposure and underexposure problems will not be fully resolved; if the threshold is too tight, over-optimization or even failure to converge will occur. The above problems together result in low convergence efficiency and poor stability in the iterative optimization process of the spherical harmonic coefficient, and the final rendering effect is difficult to meet the visual consistency requirements of HDR large scenes. Therefore, it needs to be further explained that this embodiment performs iterative optimization on the spherical harmonic coefficients of each Gaussian sphere in the 3D Gaussian splash model, including:

[0146] S301. Based on the indoor and outdoor area segmentation mask, the Gaussian spheres in the three-dimensional Gaussian splash model are divided into an indoor area Gaussian sphere set, an outdoor area Gaussian sphere set, and a boundary Gaussian sphere. A first optimization step size is assigned to the indoor area Gaussian sphere set, and a second optimization step size is assigned to the outdoor area Gaussian sphere set, wherein the first optimization step size is smaller than the second optimization step size.

[0147] Further explanation is needed; please refer to [link / reference]. Figure 2 In this embodiment, the Gaussian spheres in the three-dimensional Gaussian splash model are divided into an indoor region Gaussian sphere set, an outdoor region Gaussian sphere set, and a boundary Gaussian sphere, including:

[0148] For each Gaussian sphere, obtain its center point coordinates and the spatial extension range defined by the covariance matrix, wherein the spatial extension range is an ellipsoidal region whose semi-axis length is determined by the eigenvalues ​​of the covariance matrix;

[0149] Based on the spatial extension range, calculate the set of voxels covered by the Gaussian sphere in the indoor and outdoor region segmentation mask, where the set of voxels is all voxels that intersect with the ellipsoidal region;

[0150] The number of indoor voxels and the number of outdoor voxels in the voxel set are counted, and the proportion of indoor voxels and the proportion of outdoor voxels are calculated respectively.

[0151] When the proportion of indoor voxels is greater than or equal to the first classification threshold, the Gaussian sphere is classified into the indoor region Gaussian sphere set; when the proportion of outdoor voxels is greater than or equal to the second classification threshold, the Gaussian sphere is classified into the outdoor region Gaussian sphere set. In this embodiment, the first classification threshold and the second classification threshold are both preset values ​​between 0.5 and 0.8. This is to ensure that the Gaussian sphere is classified as indoor or outdoor only when the proportion of voxels in a certain type of region exceeds half and reaches sufficient confidence, avoiding misclassification due to insufficient proportion of a single type. At the same time, a buffer zone is reserved to identify Gaussian spheres that cover both indoor and outdoor regions as boundary Gaussian spheres for special processing, thereby achieving a smooth transition between indoor and outdoor lighting.

[0152] When both the indoor and outdoor voxel proportions are less than the corresponding division thresholds, the Gaussian sphere is marked as a boundary Gaussian sphere, and a third optimization step size between the first and second optimization step sizes is assigned to the boundary Gaussian sphere. The third optimization step size is obtained by linear interpolation based on the ratio of the indoor voxel proportion to the outdoor voxel proportion.

[0153] During the iterative optimization process, the above division steps are re-executed every preset number of rounds to update the category assignment of the Gaussian sphere. The preset number of rounds is an integer between 10 and 20, which aims to balance the timeliness of updating the category assignment of the Gaussian sphere and the computational overhead during the iterative optimization process. This ensures that the region assignment can be re-divided in a timely manner when there are significant changes in the indoor and outdoor lighting boundaries, while avoiding unnecessary computational redundancy caused by overly frequent updates.

[0154] S302. Select the current viewing direction sequentially from multiple preset viewing directions, render the current rendering frame of the current viewing direction based on the initial spherical harmonic coefficients obtained from initialization, and map the pixel radiance value of the current rendering frame to the display brightness value according to the dynamic range compression curve.

[0155] Furthermore, in this embodiment, if the preset multiple viewing directions are not distributed reasonably, the optimization result may only meet the requirements for overexposure and underexposure areas in some of the viewing directions involved in the optimization, while the rendering effect fails in other viewing directions not involved in the optimization. Therefore, this embodiment selects the current viewing direction sequentially from the preset multiple viewing directions, specifically including:

[0156] Obtain the 3D bounding box of the target scene, and calculate the coordinates of the center point and the length of the longest diagonal of the bounding box; construct an observation sphere with the center point as the center and a radius of 1.5 to 2 times the length of the longest diagonal; uniformly sample and generate N candidate view directions on the observation sphere, where N is a preset integer greater than or equal to 12; for each candidate view direction, calculate the angle between the vector pointing to the center point and the normal vector of the main facade of the scene; select candidate view directions with an angle between 30 degrees and 150 degrees as effective view directions; select M view directions from the effective view directions at equal intervals as a preset set of multiple view directions, where M is a preset integer greater than or equal to 6; during the iterative optimization process, select the current view direction in sequence according to the preset order of the M view directions, and repeat the selection cycle after completing one round of M view directions.

[0157] S303. Based on preset overexposure and underexposure thresholds, mark each pixel of the current rendering frame with an overexposure or underexposure flag, calculate the first area ratio of the overexposure region and the second area ratio of the underexposure region in the current rendering frame, and calculate the difference between the first area ratio and the second area ratio as the current difference.

[0158] Furthermore, in this embodiment, the area ratio statistics of overexposed and underexposed areas need to distinguish between areas of different importance. For example, the overexposed and underexposed control of the main scene area (such as the gate of the Hall of Supreme Harmony) should be more stringent than that of the background area (such as the sky); otherwise, it will lead to the loss of detail in the main scene area and over-optimization in the background area. Therefore, this embodiment calculates the first area ratio of overexposed areas and the second area ratio of underexposed areas in the current rendering frame, specifically including:

[0159] Based on the indoor / outdoor region segmentation mask, the pixels of the current rendered frame are divided into an indoor pixel set and an outdoor pixel set. The number of overexposed pixels and underexposed pixels are counted for each set. The proportion of overexposed pixels in the indoor region to the total number of indoor pixels is calculated as the first indoor overexposure proportion, and the proportion of underexposed pixels in the indoor region to the total number of indoor pixels is calculated as the first indoor underexposure proportion. Similarly, the proportion of overexposed pixels in the outdoor region to the total number of outdoor pixels is calculated as the first outdoor overexposure proportion, and the proportion of underexposed pixels in the outdoor region to the total number of outdoor pixels is calculated as the first outdoor underexposure proportion. The first indoor overexposure proportion and the first outdoor overexposure proportion are weighted and summed to obtain the first area proportion, where the weight coefficient for the indoor region is greater than that for the outdoor region. The first indoor underexposure proportion and the first outdoor underexposure proportion are weighted and summed to obtain the second area proportion, where the weight coefficient for the indoor region is greater than that for the outdoor region.

[0160] It should be further explained that this embodiment marks each pixel of the current rendered frame with an overexposure or underexposure flag, and adopts a staged dynamic threshold adjustment strategy to adapt to the optimization needs of different iteration stages, specifically including:

[0161] Initialize the overexposure threshold to 240 and the underexposure threshold to 30 for the first iteration stage; initialize the overexposure threshold to 230 and the underexposure threshold to 25 for the second iteration stage; initialize the overexposure threshold to 220 and the underexposure threshold to 20 for the third iteration stage; select the corresponding overexposure threshold and underexposure threshold according to the iteration stage of the current iteration number k; compare the display brightness value of each pixel in the current rendering frame with the overexposure threshold of the current stage, and mark pixels with a value greater than the threshold as overexposed pixels; compare the display brightness value of each pixel in the current rendering frame with the underexposure threshold of the current stage, and mark pixels with a value less than the threshold as underexposed pixels; record the overexposure or underexposure flag of each pixel for subsequent area proportion statistics.

[0162] S304. Based on the diffuse reflection coefficient and specular reflection coefficient in the material reflectivity model, calculate the spherical harmonic coefficient gradients of the Gaussian sphere in the indoor area, the Gaussian sphere in the outdoor area, and the boundary Gaussian sphere, respectively, and update the spherical harmonic coefficients of the corresponding Gaussian spheres using the first optimization step size, the second optimization step size, and the third optimization step size.

[0163] It should be further explained that this embodiment calculates the spherical harmonic coefficient gradients of the Gaussian sphere in the indoor region, the Gaussian sphere in the outdoor region, and the boundary Gaussian sphere, specifically including:

[0164] For each Gaussian sphere in the set of Gaussian spheres in the indoor region, its corresponding diffuse reflection coefficient is obtained; based on the Lambertian diffuse reflection model, the theoretical radiance value of the Gaussian sphere under the current viewing direction is calculated, which is equal to the diffuse reflection coefficient multiplied by the dot product of the incident light direction and the normal vector of the Gaussian sphere; the difference between the actual radiance value of the corresponding pixel of the Gaussian sphere in the current rendering frame and the theoretical radiance value is calculated as the first radiance residual; the first radiance residual is backpropagated to the spherical harmonic coefficients of the Gaussian sphere to obtain the indoor spherical harmonic coefficient gradient dominated by low-frequency components; in a specific embodiment of the present invention, the process of backpropagating the first radiance residual to the spherical harmonic coefficients of the Gaussian sphere in the indoor region is as follows: for the current Gaussian sphere in the set of Gaussian spheres in the indoor region... A Gaussian sphere is used to obtain its actual rendered radiance value under the current viewpoint direction. Simultaneously, the theoretical radiance value of the Gaussian sphere under the current viewpoint direction is calculated based on the Lambert diffuse reflection model. The calculation of the Lambert diffuse reflection model is a well-known technique in the field. This theoretical radiance value is determined by those skilled in the art based on the diffuse reflection coefficient, incident radiance, normal vector, and incident light direction of the Gaussian sphere, according to the standard calculation method of Lambert's cosine law. That is, the theoretical radiance value is equal to the diffuse reflection coefficient of the Gaussian sphere multiplied by the incident radiance, and then multiplied by the larger of the dot product of the Gaussian sphere's normal vector and the incident light direction, and 0. The actual rendered radiance value is subtracted from the theoretical radiance value to obtain the first radiance residual, and this first radiance residual is used as the basis for further calculation. The square of the radiance residual is used as the loss function value. Using the chain rule for multivariate functions, the actual rendered radiance value is expressed as the sum of the product of each order of spherical harmonic coefficient and the corresponding order of real spherical harmonic basis function in the current observation direction. Each order is determined by the order l and degree m of the spherical harmonic function, and the range of l and m is determined by the preset highest spherical harmonic order. The partial derivative of any l-th order m-th spherical harmonic coefficient is calculated, yielding a gradient value that is twice the first radiance residual multiplied by the value of the corresponding l-th order m-th real spherical harmonic basis function in the current observation direction. Subsequently, the calculated full-order spherical harmonic coefficient gradient is weighted according to indoor lighting characteristics, and the gradient values ​​of all low-frequency components with orders l of 0, 1, and 2 are... Multiply the gradient values ​​of all high-frequency components with a spherical harmonic coefficient order l of 3 or higher by a preset first low-frequency weighting coefficient, which is exemplarily set to 1.5. This yields the indoor region spherical harmonic coefficient gradient dominated by low-frequency components. Finally, perform numerical validity checks and boundary clamping on each weighted gradient component. If a gradient component is not a NaN or is infinite, set its value to 0. If a gradient component is greater than a preset upper gradient threshold of 1.0, clamp its value to 1.0. If a gradient component is less than a preset lower gradient threshold of -1.0, clamp its value to -1.The gradient is set to 0, while the remaining valid gradient components remain unchanged. All gradient components, after the above verification and clamping processes, are used as the spherical harmonic coefficient update gradient for the Gaussian sphere in the current iteration step for this indoor region, providing a low-frequency-dominated gradient direction and amplitude basis for the iterative update of the spherical harmonic coefficients.

[0165] For each Gaussian sphere in the outdoor area Gaussian sphere set, its corresponding roughness and metallicity are obtained; based on the Cook-Torrance specular reflection model, combined with the current viewing direction and the incident light direction, the theoretical radiosity value of the Gaussian sphere under the current viewing direction is calculated; the difference between the actual radiosity value of the corresponding pixel of the Gaussian sphere in the current rendering frame and the theoretical radiosity value is calculated as the second radiosity residual; the second radiosity residual is backpropagated to the spherical harmonic coefficients of the Gaussian sphere to obtain the outdoor spherical harmonic coefficient gradient containing high-frequency components; in a specific embodiment of the present invention, the process of backpropagating the second radiosity residual to the spherical harmonic coefficients of the outdoor area Gaussian sphere is as follows: for the outdoor area Gaussian sphere For the current Gaussian sphere in the set, obtain its actual rendered radiance value under the current viewpoint direction. Simultaneously, calculate the theoretical radiance value of the Gaussian sphere under the current viewpoint direction based on the Cook-Torrance microplane specular reflection model. This theoretical radiance value is composed of the diffuse reflection component and the specular reflection component. The calculation of the specular reflection component is a well-known technique in the art, determined by those skilled in the art based on the roughness parameters, metallicity parameters, normal vector, incident light direction, and current viewpoint direction of the Gaussian sphere, according to the standard calculation method of the normal distribution function, geometric occlusion function, and Fresnel term in the Cook-Torrance model. Subtract the theoretical radiance value from the actual rendered radiance value. The radiance value is used to obtain the second radiance residual, and the square of this second radiance residual is used as the loss function value. Using the chain rule of multivariate functions, the actual rendered radiance value is expressed as the sum of the product of each order of spherical harmonic coefficient and the corresponding order of real spherical harmonic basis function in the current observation direction. The partial derivative is taken for any l-th order m-th spherical harmonic coefficient, and the gradient value of this spherical harmonic coefficient is twice the second radiance residual multiplied by the value of the real spherical harmonic basis function corresponding to the l-th order m-th order in the current observation direction. Then, the calculated full-order spherical harmonic coefficient gradient is subjected to outdoor illumination characteristic adaptation weighting processing, and the gradient values ​​of all low-frequency components with spherical harmonic coefficient order l of 0, 1, and 2 are multiplied by the preset second low-frequency weight. The second low-frequency weighting coefficient is exemplarily set to 0.5. The gradient values ​​of all high-frequency components with a spherical harmonic coefficient order l of 3 or higher are multiplied by the preset second high-frequency weighting coefficient, which is exemplarily set to 2.0. This yields the outdoor area spherical harmonic coefficient gradient dominated by high-frequency components. Finally, each weighted gradient component undergoes numerical validity verification and boundary clamping. If a gradient component is not a NaN or is infinite, its value is set to 0. If a gradient component is greater than the preset upper gradient threshold of 1.0, its value is clamped to 1.0. If a gradient component is less than the preset lower gradient threshold of -1.0, its value is clamped to -1.The gradient is set to 0, while the remaining valid gradient components remain unchanged. All gradient components, after the above verification and clamping process, are used as the spherical harmonic coefficient update gradient for the Gaussian sphere in the current iteration step for this outdoor region. This provides a basis for the gradient direction and amplitude, focusing on high-frequency details, for the iterative update of the spherical harmonic coefficients.

[0166] For each Gaussian sphere in the boundary Gaussian sphere set, obtain its corresponding diffuse reflection coefficient, roughness, and metallicity, as well as the indoor and outdoor voxel proportions of the Gaussian sphere; calculate the first theoretical radiosity value of the Gaussian sphere in the current viewing direction based on the Lambertian diffuse reflection model, and calculate the second theoretical radiosity value of the Gaussian sphere in the current viewing direction based on the Cook-Torrance specular reflection model; weight and fuse the first and second theoretical radiosity values ​​according to the indoor and outdoor voxel proportions to obtain the comprehensive theoretical radiosity value of the Gaussian sphere; calculate the actual radiosity value of the corresponding pixel of the Gaussian sphere in the current rendering frame. The difference between the radiance value and the comprehensive theoretical radiance value is used as the third radiance residual; the third radiance residual is backpropagated to the spherical harmonic coefficients of the Gaussian sphere to obtain the boundary spherical harmonic coefficient gradient; in a specific embodiment of the present invention, the process of backpropagating the third radiance residual to the spherical harmonic coefficients of the boundary Gaussian sphere is as follows: for the current Gaussian sphere marked as the boundary Gaussian sphere, firstly, its actual rendered radiance value under the current viewing direction is obtained, and simultaneously, the first theoretical radiance value of the boundary Gaussian sphere under the current viewing direction is calculated based on the Lambertian diffuse reflection model. The calculation of the Lambertian diffuse reflection model is a well-known technique in the art, and the first theoretical radiance value is obtained by... Those skilled in the art determine the first theoretical radiance value based on the diffuse reflection coefficient, incident radiance, normal vector, and incident illumination direction of the boundary Gaussian sphere, using the standard calculation method of Lambert's cosine law. That is, the first theoretical radiance value is equal to the larger of the diffuse reflection coefficient of the boundary Gaussian sphere multiplied by the incident radiance, and then multiplied by the larger of the dot product of the normal vector and the incident illumination direction of the boundary Gaussian sphere and 0. Furthermore, based on the Cook-Torrance microplane specular reflection model, the second theoretical radiance value of the boundary Gaussian sphere at the current viewing angle is calculated. This second theoretical radiance value is composed of the addition of the diffuse reflection component and the specular reflection component, where the specular reflection component... The calculation is a well-known technique in the field, and is determined by those skilled in the art based on the roughness parameters, metallicity parameters, normal vector, incident light direction, and current viewing direction of the boundary Gaussian sphere, according to the standard calculation method of the normal distribution function, geometric occlusion function, and Fresnel term in the Cook-Torrance model; subsequently, the indoor voxel ratio and outdoor voxel ratio calculated during the region division of the boundary Gaussian sphere are obtained, and the indoor voxel ratio is used as the first fusion weight and the outdoor voxel ratio is used as the second fusion weight to perform a weighted sum of the first theoretical radiance value and the second theoretical radiance value to obtain the comprehensive theoretical radiance value.Subtracting the comprehensive theoretical radiance value from the actual rendered radiance value yields the third radiance residual, and the square of this third radiance residual is used as the loss function value. Using the chain rule for multivariate functions, the actual rendered radiance value is expressed as the sum of the products of each order of spherical harmonic coefficient and the corresponding order of real spherical harmonic basis function in the current observation direction. Taking the partial derivative of any l-th order m-th spherical harmonic coefficient yields a gradient value that is twice the third radiance residual multiplied by the value of the l-th order m-th corresponding real spherical harmonic basis function in the current observation direction. Subsequently, a weighted processing for boundary lighting transition adaptation is performed on the calculated full-order spherical harmonic coefficient gradient. The gradient values ​​of all low-frequency components with spherical harmonic coefficient orders l of 0, 1, and 2 are multiplied by the indoor voxel percentage of the boundary Gaussian sphere, and the gradient values ​​of all high-frequency components with spherical harmonic coefficient orders l of 3 and above are multiplied by the outdoor voxel percentage of the boundary Gaussian sphere. This yields a boundary region spherical harmonic coefficient gradient that balances smooth transitions between indoor and outdoor lighting. Finally, each weighted gradient component undergoes numerical validity verification and boundary clamping. If a gradient component is non-numeric (NaN) or infinite, its value is set to 0. If a gradient component is greater than the preset upper gradient threshold of 1.0, its value is clamped to 1.0. If a gradient component is less than the preset lower gradient threshold of -1.0, its value is clamped to -1.0. The remaining valid gradient components remain unchanged. All gradient components after the above verification and clamping are used as the spherical harmonic coefficient update gradient for the current iteration step of the boundary Gaussian sphere, providing a gradient direction and amplitude basis for the iterative update of the spherical harmonic coefficient to adapt to the transition characteristics of indoor and outdoor lighting.

[0167] The indoor spherical harmonic coefficient gradient, the outdoor spherical harmonic coefficient gradient, and the boundary spherical harmonic coefficient gradient are stored respectively.

[0168] It should be further explained that, in this embodiment, when updating the spherical harmonic coefficients using the first, second, and third optimization step sizes, a spatial smoothing regularization process is introduced to maintain the continuity of the spherical harmonic coefficients between adjacent Gaussian spheres. Specifically, this includes:

[0169] For each Gaussian sphere, construct its neighborhood Gaussian sphere set. The neighborhood is defined as other Gaussian spheres within a sphere centered at the center of the current Gaussian sphere and with a preset neighborhood radius of 0.05 meters to 0.1 meters. Calculate the difference in spherical harmonic coefficients between the current Gaussian sphere and each of its neighboring Gaussian spheres. This difference is the sum of the squares of the differences between the corresponding components of each order of spherical harmonic coefficient. Sum the difference values ​​of all neighboring Gaussian spheres to obtain a spatial inconsistency measure for the Gaussian sphere. Calculate a regularization gradient based on this spatial inconsistency measure, with the gradient direction pointing towards reducing the difference with neighboring Gaussian spheres. The direction is as follows: Multiply the regularization gradient by a preset regularization weight coefficient to obtain the regularization correction amount. The regularization weight coefficient is a preset value between 0.01 and 0.1, which aims to balance the spatial smoothing effect and the fidelity of the original lighting data. While effectively suppressing abrupt changes in spherical harmonic coefficients between adjacent Gaussian spheres and avoiding visual discontinuities in the rendering results, it also prevents excessive smoothing from causing the loss of high-frequency lighting details. Add the regularization correction amount to the spherical harmonic coefficient gradient calculated from the radiosity residual to obtain the comprehensive gradient. Use the comprehensive gradient to replace the original gradient and update the spherical harmonic coefficients in combination with the corresponding optimization step size.

[0170] It should be further explained that, in this embodiment, when updating the spherical harmonic coefficients of the indoor region, outdoor region, and boundary Gaussian sphere using the first, second, and third optimization step sizes, the step size decreases with the number of iterations to facilitate convergence, specifically including:

[0171] Initialize the current iteration count k=1, and set the maximum iteration count K_max to a preset value between 100 and 200; calculate the first decay factor, the second decay factor, and the third decay factor based on the current iteration count k, wherein the decay factors decrease as k increases; specifically, an exponential decay strategy is adopted: first decay factor = α1 k Second attenuation factor = α2 k The third attenuation factor = α3 k α1, α2, and α3 are the preset decay base rates of the indoor region, outdoor region, and boundary Gaussian sphere, respectively, and all are preset values ​​between 0.95 and 0.99. The aim is to make the optimization step size decrease exponentially with the number of iterations, maintain a large step size in the early stage of iteration to quickly approach the optimal solution, and make fine adjustments with a small step size in the later stage of iteration, thereby balancing convergence speed and stability.

[0172] Multiply the first optimization step size by the first decay factor to obtain the first actual step size of the current iteration; multiply the second optimization step size by the second decay factor to obtain the second actual step size of the current iteration; multiply the third optimization step size by the third decay factor to obtain the third actual step size of the current iteration; multiply the indoor spherical harmonic coefficient gradient by the first actual step size to obtain the indoor spherical harmonic coefficient update amount, and add it to the spherical harmonic coefficient of the current indoor region Gaussian sphere to obtain the updated indoor spherical harmonic coefficient; multiply the outdoor spherical harmonic coefficient gradient by the second actual step size to obtain the outdoor spherical harmonic coefficient update amount, and add it to the spherical harmonic coefficient of the current outdoor region Gaussian sphere to obtain the updated outdoor spherical harmonic coefficient; multiply the boundary spherical harmonic coefficient gradient by the third actual step size to obtain the boundary spherical harmonic coefficient update amount, and add it to the spherical harmonic coefficient of the current boundary Gaussian sphere to obtain the updated boundary spherical harmonic coefficient; record the various spherical harmonic coefficients after this iteration for rendering in the next iteration cycle.

[0173] It should be further explained that, in this embodiment, when updating the spherical harmonic coefficients using the first, second, and third optimization step sizes, a differentiated step size adjustment is also performed based on the order of the spherical harmonic coefficients to balance the convergence speed of low-frequency and high-frequency components. Specifically, this includes:

[0174] The spherical harmonic coefficients are divided into low-frequency and high-frequency groups according to their order l, where l = 0 to 2 is the low-frequency group and l = 3 and above is the high-frequency group. A first-order step size factor is assigned to the low-frequency group, and a second-order step size factor is assigned to the high-frequency group. The first-order step size factor is greater than the second-order step size factor. In each iteration, the first, second, and third actual step sizes are multiplied by the first-order step size factor to obtain the actual step size of the low-frequency group. The first, second, and third actual step sizes are multiplied by the second-order step size factor to obtain the actual step size of the high-frequency group. The gradients of the spherical harmonic coefficients of the low-frequency group and the high-frequency group are calculated for each type of Gaussian sphere. The low-frequency group spherical harmonic coefficients are updated using the corresponding low-frequency group actual step size, and the high-frequency group spherical harmonic coefficients are updated using the corresponding high-frequency group actual step size.

[0175] S305. Repeat the step of sequentially selecting the current viewing direction from multiple preset viewing directions until the current difference is less than the first preset threshold in three or more consecutive viewing directions. Then stop the iterative optimization and output the illumination adaptation completion signal.

[0176] It should be further noted that, to ensure the stability of convergence and avoid premature termination due to accidental fluctuations, this embodiment has refined the stopping conditions, specifically including:

[0177] A sliding window of length L is set, where L is a preset integer greater than or equal to 3. After each calculation of a current difference, it is stored in the sliding window, and the oldest difference record in the window is removed. When the sliding window is full, it is determined whether all L current differences in the window are less than the first preset threshold. If they are all less than the first preset threshold, the ratio of the mean to the standard deviation of the L differences in the window is further calculated. When the ratio of the mean to the standard deviation is less than the first stability threshold, it is determined that the current difference has stabilized and converged. At the same time, it is determined whether the mean of the differences in the most recent two consecutive rounds of M view directions shows a downward trend. If it shows an upward trend, even if the current differences are all less than the first preset threshold, the iteration does not stop, and the optimization continues. Only when all differences in the window are less than the first preset threshold, the ratio of the mean to the standard deviation is less than the first stability threshold, and the mean of the differences shows a downward trend or remains flat, is the illumination adaptation completion signal output.

[0178] It should be further noted that, in order to verify the effectiveness of the optimization results in all preset viewing directions, this embodiment introduces a periodic full-view verification mechanism, specifically including:

[0179] Every T iterations (T being a preset integer between 5 and 10), a full-view verification is performed. During the full-view verification, M verification frames are rendered sequentially from M preset viewpoints. The overexposed and underexposed area percentages of each verification frame are calculated, and the difference between each frame is determined. It is then determined whether the differences between the M verification frames are all less than the first preset threshold. If the differences between the M verification frames are all less than the first preset threshold, a full-view verification pass flag is recorded. If the difference between any verification frame is greater than or equal to the first preset threshold, iterative optimization continues, and the count for three consecutive viewpoints meeting the standard is reset. Only when three or more consecutive viewpoints meet the standard simultaneously and the most recent full-view verification passes is a lighting adaptation completion signal output.

[0180] Next, a detailed description of the above-described iterative optimization process for spherical harmonic coefficients will be provided using a specific and complete embodiment. This embodiment is only used to illustrate the feasibility of the technical solution of the present invention at the computational level and does not represent a limitation on the scope of protection of the present invention. Those skilled in the art can determine the specific parameter values ​​based on the actual application scenario through simulation experiments or physical experiments. For example, this embodiment uses a digital farmland scenario in a smart agriculture demonstration zone as the application object to illustrate the specific implementation of the above-described iterative optimization process:

[0181] Based on the crop-soil segmentation mask corresponding to the indoor and outdoor region segmentation mask, the 50 million Gaussian spheres in the 3D Gaussian splash model are divided into a set of Gaussian spheres for the crop region corresponding to the indoor region, a set of Gaussian spheres for the soil region corresponding to the outdoor region, and a set of Gaussian spheres for the boundary region. The set of Gaussian spheres for the crop region contains approximately 20 million Gaussian spheres, the set of Gaussian spheres for the soil region contains approximately 25 million Gaussian spheres, and the set of Gaussian spheres for the boundary region contains approximately 5 million Gaussian spheres. A first optimization step size of 0.001 is assigned to the crop region, a second optimization step size of 0.005 is assigned to the soil region, and a third optimization step size for the boundary Gaussian spheres is calculated by linear interpolation between the first and second optimization step sizes based on their indoor and outdoor voxel proportions. The formula for calculating the third optimization step size is: Third optimization step size = First optimization step size × Indoor voxel proportion + Second optimization step size × Outdoor voxel proportion. For example, if the proportion of indoor voxels in a certain boundary Gaussian sphere is 0.4 and the proportion of outdoor voxels is 0.6, then the third optimization step size for this boundary Gaussian sphere is 0.001 × 0.4 + 0.005 × 0.6 = 0.0032. The above Gaussian sphere classification and optimization step size allocation operation is re-executed every 15 iterations to adapt to the dynamic changes in the spatial distribution of Gaussian spheres during the optimization process.

[0182] The current viewpoint is selected sequentially from 12 preset viewpoint directions. These 12 viewpoint directions are constructed as follows: First, the 3D bounding box of the target farmland scene is obtained, and the coordinates of the center point of the bounding box are calculated to be 0 meters, 0 meters, and 1.5 meters, with the longest diagonal length being 100 meters. An observation sphere is constructed with this center point as the center and a radius of 1.8 times the longest diagonal length. 24 candidate viewpoint directions are generated by uniformly sampling on the observation sphere. For each candidate viewpoint direction, the angle between the vector pointing to the center point of the scene bounding box and the normal vector of the scene's main facade is calculated. Candidate viewpoint directions with an angle between 30 degrees and 150 degrees are selected as effective viewpoint directions. Finally, 12 viewpoint directions are selected from the effective viewpoint directions at equal intervals as preset multiple viewpoint directions for iterative optimization.

[0183] Based on the initial spherical harmonic coefficients obtained from initialization, a rendering frame is obtained in the current viewing direction. The pixel radiance values ​​of the current rendering frame are mapped to display brightness values ​​according to the dynamic range compression curve. The minimum radiance value of this dynamic range compression curve is 0.2, the maximum radiance value is 45.0, and the contrast adjustment factor is 2.5. Based on a phased dynamic threshold adjustment strategy, each pixel in the current rendering frame is marked with an overexposure or underexposure flag: In the early stages of iterative optimization (iterations 1-50), an overexposure brightness threshold of 240 and an underexposure brightness threshold of 30 are used; in the middle stages (iterations 51-100), the overexposure brightness threshold is tightened to 230, and the underexposure brightness threshold is tightened to 25; in the later stages (iterations 101 and beyond), the overexposure brightness threshold is further tightened to 220, and the underexposure brightness threshold is tightened to 20. The optimization accuracy is gradually improved through phased threshold adjustments. Based on the crop-soil segmentation mask, the pixels of the current rendered frame are divided into a crop pixel set and a soil pixel set. The percentages of overexposed pixels in the crop region are calculated as 0.08 and underexposed pixels as 0.02, while the percentages of overexposed pixels in the soil region are 0.12 and underexposed pixels as 0.04. A weighted summation method is used to calculate the overall area percentage, with a weight coefficient of 0.7 for the crop region and 0.3 for the soil region. The first area percentage of the overexposed region is calculated to be 0.092, and the second area percentage of the underexposed region is calculated to be 0.026. The difference between the first and second area percentages, 0.066, is taken as the current difference.

[0184] Based on the material reflectivity model, the spherical harmonic coefficient gradients of different types of Gaussian spheres are calculated: For Gaussian spheres in crop areas, the theoretical radiance value is calculated based on the Lambertian diffuse reflection model, which corresponds to a diffuse reflection coefficient of 0.65. For example, the actual rendered radiance value of a Gaussian sphere for a corn leaf is 8.5, and the theoretical radiance value is 8.2. The first radiance residual is calculated to be 0.3. Through the backpropagation process of the aforementioned first radiance residual, the spherical harmonic coefficient gradient of the crop area, dominated by low-frequency components, is obtained. For Gaussian spheres in soil areas, the theoretical radiance value is calculated based on the Cook-Torrance specular reflection model, which corresponds to a roughness parameter of 0.3 and a metallicity parameter of 0.1. For example, the actual rendered radiance value of a Gaussian sphere for soil is 11.8, and the theoretical radiance value is 12. The second radiance residual is calculated to be -0.2. Through the backpropagation process of the aforementioned second radiance residual, the spherical harmonic coefficient gradient of the soil region containing high-frequency components is obtained. For the boundary Gaussian sphere, the first theoretical radiance value is calculated to be 7.5 based on the Lambert diffuse reflection model, and the second theoretical radiance value is calculated to be 11.2 based on the Cook-Torrance specular reflection model. According to the indoor voxel ratio of 0.4 and the outdoor voxel ratio of 0.6 for the boundary Gaussian sphere, the first theoretical radiance value and the second theoretical radiance value are weighted and fused to obtain a comprehensive theoretical radiance value of 10.2. The actual rendered radiance value of the boundary Gaussian sphere is 10.5, and the third radiance residual is calculated to be 0.3. Through the backpropagation process of the aforementioned third radiance residual, the spherical harmonic coefficient gradient of the boundary region is obtained. In this embodiment, the weighted fusion weights of the boundary Gaussian sphere are directly determined by the proportion of indoor and outdoor voxels covered by its space. The weights based on the Lambertian model are equal to the indoor voxel proportions, and the weights based on the Cook-Torrance model are equal to the outdoor voxel proportions; the sum of these two weights is 1. This weight determination method is based on the continuity assumption of physical space affiliation. That is, the higher the proportion of indoor voxels in the area covered by the Gaussian sphere, the more its lighting response tends to reflect the diffuse characteristics of indoor materials; conversely, the higher the proportion of outdoor voxels, the more it tends to reflect the high-reflectivity characteristics of outdoor materials. This method of assigning weights based on voxel proportions is consistent with the linear interpolation strategy of the boundary Gaussian sphere optimization step size, achieving a smooth transition in lighting response from indoor to outdoor areas and avoiding the problem of abrupt changes in light and shadow caused by abrupt changes in material models at the boundary.

[0185] In the process of updating the spherical harmonic coefficients using the corresponding optimization step size, spatial smoothing regularization is introduced. Specifically, a neighborhood Gaussian sphere set is constructed for each Gaussian sphere with a neighborhood radius of 0.05 meters. The difference in spherical harmonic coefficients between the current Gaussian sphere and each of its neighboring Gaussian spheres is calculated. This difference is the sum of the squares of the differences between the corresponding components of each order of spherical harmonic coefficients. The summation of the difference values ​​for all neighboring Gaussian spheres yields a spatial inconsistency measure of 0.5 for the Gaussian sphere. A regularization gradient is calculated based on the spatial inconsistency measure, pointing in the direction of reducing the difference with neighboring Gaussian spheres. The regularization gradient is multiplied by a preset regularization weight coefficient of 0.05 to obtain a regularization correction of 0.025. This regularization correction is added to the spherical harmonic coefficient gradient calculated from the radiometric residual to obtain a comprehensive gradient. This comprehensive gradient replaces the original gradient, and the spherical harmonic coefficients are updated using the corresponding optimization step size. The regularization weight coefficient ranges from 0.01 to 0.1. In this embodiment, 0.05 is selected to balance the spatial smoothing effect and the fidelity of the original lighting data. While effectively suppressing the abrupt change in the spherical harmonic coefficients between adjacent Gaussian spheres and avoiding visual discontinuities in the rendering results, it also prevents excessive smoothing from causing the loss of high-frequency lighting details.

[0186] The optimization step size is dynamically adjusted using an exponential decay strategy with each iteration, and the decay base rate is set to 0.98, meaning that the actual step size in the k-th iteration is the initial step size multiplied by 0.98. k In the first iteration, the first actual step size was calculated to be 0.001 × 0.98 = 0.00098, the second actual step size was 0.005 × 0.98 = 0.0049, and the third actual step size was 0.0032 × 0.98 = 0.003136. Simultaneously, differentiated step size adjustments were made based on the order of the spherical harmonic coefficients. The spherical harmonic coefficients were divided into low-frequency and high-frequency groups, with orders 0 to 2 being the low-frequency group and orders 3 and above being the high-frequency group. A first-order step size factor of 1.0 was assigned to the low-frequency group, and a second-order step size factor of 0.5 was assigned to the high-frequency group. This ensured that the low-frequency components converged at a faster rate, while the high-frequency components underwent fine-tuning at a slower rate, avoiding numerical oscillations during the optimization process.

[0187] The iterative optimization process is repeated until a preset convergence stopping condition is met. This condition is as follows: the current difference in five consecutive viewing directions is less than a first preset threshold of 0.02; the ratio of the mean to the standard deviation of all differences within a sliding window of length 5 is less than a first stability threshold of 0.2; the mean difference in the last two consecutive rounds of twelve viewing directions shows a decreasing trend; and in every 10 rounds of full-view verification, the difference in all twelve verification rendering frames is less than the first preset threshold of 0.02. The first preset threshold of 0.02 is set based on the scene brightness distribution characteristics to ensure that the deviation in the area ratio of overexposed and underexposed regions is controlled within an acceptable range. The first stability threshold of 0.2 is used to determine the fluctuation of the difference sequence to ensure a stable and reliable convergence state. When the above convergence stopping condition is met, the iterative optimization stops and a lighting adaptation completion signal is output. At this point, the optimized spherical harmonic coefficients can accurately express the highlight reflection details of corn leaves and the diffuse reflection texture of the soil, achieving a high-quality rendering effect for large HDR scenes.

[0188] This embodiment addresses the convergence speed imbalance caused by significant differences in indoor and outdoor lighting characteristics due to the introduction of a voxel coverage statistical method based on spatial extension range. This method finely divides the Gaussian sphere into indoor, outdoor, and boundary regions and assigns differentiated optimization step sizes, enabling stable convergence in low-frequency indoor regions with smaller step sizes, rapid detail capture in high-frequency outdoor regions with larger step sizes, and smooth transition in boundary regions using interpolation step sizes. By constructing a uniform sampling multi-viewpoint structure on the observation sphere and iteratively optimizing, it solves the problem of insufficient viewpoint generalization ability caused by single-viewpoint optimization, ensuring that the spherical harmonic coefficients maintain a balance between overexposed and underexposed areas under all preset viewpoints. Furthermore, by using region-weighted area ratio statistics based on indoor and outdoor segmentation masks and dynamically adjusted overexposed and underexposed thresholds in stages, it solves the detail loss problem caused by differences in importance between the subject and background regions, focusing optimization efforts on key areas of the scene. Finally, it utilizes a Lambertian volume model and Cook-Torr... The ance model calculates the spherical harmonic coefficient gradients of indoor and outdoor Gaussian spheres separately and performs weighted fusion on the boundary Gaussian spheres, solving the gradient direction deviation problem caused by differences in material reflection characteristics, making the gradient calculation more consistent with the actual physical laws of lighting. By introducing spatial smoothing regularization, the visual discontinuity problem caused by abrupt changes in the spherical harmonic coefficients of adjacent Gaussian spheres is solved, ensuring a smooth transition in the rendering results in space. Through the step size decay mechanism and differential step size adjustment of the spherical harmonic coefficient order, the problems of oscillation in the later stage of iteration and the mismatch between low-frequency and high-frequency convergence speeds are solved, enabling rapid approximation in the early stage and fine optimization in the later stage. By judging the ratio of the mean to the standard deviation of the sliding window and periodic full-view verification, the problem of premature termination due to accidental fluctuations or failure of views not involved in optimization is solved, ensuring that the stopping conditions are both stable and reliable and fully cover all key views. Finally, it achieves efficient convergence, stable control and visually consistent rendering effects for iterative optimization of spherical harmonic coefficients in HDR large scenes.

[0189] It should be further explained that this embodiment outputs the current rendered frame based on the current spherical harmonic coefficients and a preset global tone mapping operator, specifically including:

[0190] S601. Obtain the camera intrinsic and extrinsic parameter matrices corresponding to the current viewpoint direction; for each pixel in the current rendering frame, query all Gaussian spheres that intersect with the light rays of that pixel through ray tracing or Gaussian splashing algorithm, calculate the color contribution of each Gaussian sphere in the current viewpoint direction based on the current spherical harmonic coefficients, accumulate them to obtain the radiometric value of that pixel, and generate the radiometric image in the current viewpoint direction.

[0191] S602. Read all configuration parameters of the dynamic range compression curve from the current global tone mapping operator, including minimum radiance, maximum radiance, minimum displayable brightness, maximum displayable brightness, first radiance value, second radiance value, and contrast adjustment factor.

[0192] S603. For each pixel in the radiometric image, compare its radiometric value with the minimum radiometric value and the maximum radiometric value. If it is less than the minimum radiometric value, set it to be equal to the minimum radiometric value. If it is greater than the maximum radiometric value, set it to be equal to the maximum radiometric value.

[0193] S604. Determine the brightness range to which the pixel belongs based on the first radiance value and the second radiance value. If the radiance value is less than the first radiance value, use the first mapping function to calculate the display brightness value. If the radiance value is within the range of the first radiance value and the second radiance value, use the second mapping function to calculate the display brightness value. If the radiance value is greater than the second radiance value, use the third mapping function to calculate the display brightness value. Use the calculated display brightness value as the final output brightness value of the pixel, combine all pixels to generate the current rendering frame, and output it to the display device or storage medium.

[0194] During real-time rendering, due to the dynamic changes in lighting conditions with viewpoint switching and scene content in high dynamic range scenes, the global tone mapping operator based on statically optimized spherical harmonic coefficients and fixed parameters is difficult to adapt in real time to the distribution changes of overexposed and underexposed areas in each rendering frame, resulting in frequent brightness imbalances in the rendering frames: highlight details are lost in overexposed areas, and dark textures are invisible in underexposed areas. Moreover, the control of overexposed and underexposed areas is coupled, and adjusting the spherical harmonic coefficients or tone mapping operator alone can easily cause adjustment oscillations or slow convergence. At the same time, the response characteristics of different areas (indoor, outdoor, and boundary) to brightness adjustment vary significantly, and there is a lack of differentiated adjustment mechanisms for area type and spherical harmonic coefficient frequency bands (low frequency and high frequency), making it difficult to achieve fine control of brightness balance, which seriously affects the visual continuity and realism of large-scene virtual roaming. Therefore, it needs to be further explained that this embodiment calculates the exposure fusion weight adjustment amount and tone mapping curve correction amount based on the pixel-level brightness distribution data, including:

[0195] S701. Obtain the pixel-level brightness distribution data of the current rendering frame, and calculate the proportion of overexposed pixels to the total number of pixels as the first overexposure ratio, and calculate the proportion of underexposed pixels to the total number of pixels as the first underexposure ratio.

[0196] S702. Calculate the ratio of the first overexposure ratio to the first underexposure ratio as the first scaling factor;

[0197] S703. When the first overexposure ratio is greater than the first underexposure ratio, a first gain coefficient is calculated based on the first ratio factor. The first gain coefficient is equal to the first ratio factor multiplied by the first preset gain. The first gain coefficient is used as the exposure fusion weight adjustment amount. The exposure fusion weight adjustment amount is used to enhance the weight of the high-frequency components in the spherical harmonic coefficient of the Gaussian sphere corresponding to the overexposure area.

[0198] S704. When the first underexposure ratio is greater than the first overexposure ratio, a second gain coefficient is calculated based on the reciprocal of the first ratio factor. The second gain coefficient is equal to the reciprocal of the first ratio factor multiplied by a second preset gain. The second gain coefficient is used as the exposure fusion weight adjustment amount. The exposure fusion weight adjustment amount is used to enhance the weight of the low-frequency components in the spherical harmonic coefficients of the Gaussian sphere corresponding to the underexposure area. The first preset gain is a preset constant between 0.2 and 0.5, used to control the enhancement amplitude of the high-frequency components. The second preset gain is a preset constant between 0.2 and 0.5, used to control the enhancement amplitude of the low-frequency components. In this embodiment, the spherical harmonic coefficients are divided into low-frequency components and high-frequency components according to the order l. The spherical harmonic coefficients corresponding to orders l=0 to 2 are low-frequency components, used to express the information of smooth areas where the light intensity changes slowly with the viewing angle, including overall brightness, diffuse hue, and shadow transition. The spherical harmonic coefficients corresponding to orders l=3 and above are high-frequency components, used to express the information of detailed areas where the light intensity changes rapidly with the viewing angle, including specular reflection, specular gloss, and texture edges. The determination of the low-frequency and high-frequency components is based on the definition of the order of the spherical harmonic function. The lower the order, the lower the spatial frequency, and the higher the order, the higher the spatial frequency.

[0199] S705. Calculate the deviation between the brightness ratio of the current rendered frame and the preset second threshold, multiply the deviation by the third preset gain to obtain the tone mapping curve correction amount, and use the tone mapping curve correction amount to adjust the contrast adjustment factor of the dynamic range compression curve.

[0200] S706. When both the first overexposure ratio and the first underexposure ratio exceed the third ratio threshold, the exposure fusion weight adjustment amount and the tone mapping curve correction amount are calculated simultaneously, and the third ratio threshold is a preset value between 0.1 and 0.2.

[0201] It should be further explained that in this embodiment, the exposure blending weight adjustment is applied to the spherical harmonic coefficient of the Gaussian sphere, specifically including:

[0202] S707. For each Gaussian sphere in the set of Gaussian spheres in the indoor area, increase the weight of the 0th to 2nd order low-frequency components in its spherical harmonic coefficients according to the exposure fusion weight adjustment amount. The increase of the low-frequency components is the exposure fusion weight adjustment amount multiplied by the first region coefficient, where the first region coefficient is a preset value between 0.8 and 1.0. The increase method is to multiply the 0th to 2nd order low-frequency components in the current spherical harmonic coefficients by (1 + the increase of the low-frequency components) to obtain the updated low-frequency components. If the updated low-frequency components exceed the preset upper limit of the amplitude of the spherical harmonic coefficients, then truncate them to the upper limit of the amplitude.

[0203] S708. For each Gaussian sphere in the outdoor area Gaussian sphere set, increase the weight of the 3rd order and above high-frequency components in its spherical harmonic coefficients according to the exposure fusion weight adjustment amount. The increase of the high-frequency components is the exposure fusion weight adjustment amount multiplied by the second region coefficient, where the second region coefficient is a preset value between 0.8 and 1.0. The increase method is to multiply the 3rd order and above high-frequency components in the current spherical harmonic coefficients by (1 + the increase of the high-frequency components) to obtain the updated high-frequency components. If the updated high-frequency components exceed the preset upper limit of the spherical harmonic coefficient amplitude, then truncate them to the upper limit of the amplitude.

[0204] S709. For each Gaussian sphere in the boundary Gaussian sphere set, the weights of the low-frequency and high-frequency components in its spherical harmonic coefficients are increased simultaneously according to the exposure fusion weight adjustment amount. The increase in the low-frequency component is the exposure fusion weight adjustment amount multiplied by the indoor voxel ratio of the Gaussian sphere, and the increase in the high-frequency component is the exposure fusion weight adjustment amount multiplied by the outdoor voxel ratio of the Gaussian sphere. The increase method is to multiply the 0th to 2nd order low-frequency components in the current spherical harmonic coefficients by (1 + the increase in the low-frequency component), and multiply the 3rd order and above high-frequency components in the current spherical harmonic coefficients by (1 + the increase in the high-frequency component), respectively, to obtain the updated low-frequency and high-frequency components. If any of the updated components exceeds the preset upper limit of the spherical harmonic coefficient amplitude, it is truncated to the upper limit of the amplitude. In this embodiment, the increase in the low-frequency component and the increase in the high-frequency component are both dimensionless scaling factors.

[0205] It should be further explained that in this embodiment, the tone mapping curve correction is applied to the global tone mapping operator, specifically including:

[0206] S710. Read the current contrast adjustment factor from the current global tone mapping operator, add the contrast adjustment factor to the tone mapping curve correction amount, and obtain the updated contrast adjustment factor.

[0207] S711. Write the updated contrast adjustment factor into the global tone mapping operator to replace the original contrast adjustment factor, while keeping the minimum radiosity, maximum radiosity, minimum displayable brightness, maximum displayable brightness, first radiosity value, and second radiosity value unchanged in the dynamic range compression curve.

[0208] S712. Return to the step of performing real-time rendering of the three-dimensional Gaussian splash model at a preset rendering frame rate, regenerate the current rendering frame based on the updated spherical harmonic coefficients and the updated global tone mapping operator, and recalculate the brightness ratio; set the maximum number of iterations to Q times, if the brightness ratio is less than the second threshold within the maximum number of iterations, terminate the adjustment and output the current rendering frame, if the brightness ratio is still not less than the second threshold after exceeding the maximum number of iterations, force the output of the current rendering frame and record alarm information.

[0209] Next, a specific and complete example will be used to illustrate the entire process. This example is only to illustrate the feasibility at the computational level and does not represent actual values. The specific values ​​can be determined by those skilled in the art through simulation experiments or physical experiments. For example, this embodiment takes a digital farmland in a smart agriculture demonstration area as an example. After completing the iterative optimization of the spherical harmonic coefficients and outputting the illumination adaptation completion signal, the system renders the three-dimensional Gaussian splash model of the corn farmland in real time at a preset frame rate of 60 frames per second. When the rendering perspective quickly switches from high altitude to close-up of the corn canopy, the pixel-level brightness distribution data of the current rendering frame shows that the proportion of overexposed pixels is 0.18 and the proportion of underexposed pixels is 0.04. The first overexposed proportion is greater than the first underexposed proportion, and the first proportion factor is calculated to be 4.5. Multiplying this first proportion factor by the first preset gain of 0.3, the exposure fusion weight adjustment amount of 1.35 is obtained. The first preset gain of 0.3 is used to control the enhancement amplitude of high-frequency components to avoid over-adjustment that may cause new brightness imbalances. Simultaneously, the deviation of 1.0 between the current rendering frame's brightness ratio of 2.8 and the second threshold of 1.8 is calculated. This deviation is multiplied by the third preset gain of 0.5 to obtain a tone mapping curve correction amount of 0.5. The second threshold of 1.8 is a preset brightness balance judgment benchmark, and the third preset gain of 0.5 is used to control the adjustment step size of the tone mapping curve to prevent excessive correction from causing abrupt contrast changes. Since both the overexposure ratio of 0.18 and the underexposure ratio of 0.04 exceed the third ratio threshold of 0.15, both adjustment amounts are calculated simultaneously. The third ratio threshold of 0.15 is used to determine whether it is necessary to jointly adjust the spherical harmonic coefficient and the tone mapping operator to cope with complex scenes where overexposure and underexposure coexist.

[0210] An exposure blending weight adjustment of 1.35 is applied to the spherical harmonic coefficients of the Gaussian sphere. For the soil region Gaussian sphere (outdoor region), the increase in high-frequency components is 1.35 multiplied by the second region coefficient 0.9, which equals 1.215, or 2.215 times the high-frequency components. The second region coefficient 0.9 controls the enhancement intensity of high-frequency components in the outdoor region to highlight the highlight texture of the soil under direct sunlight. For the crop region Gaussian sphere (indoor region), the increase in low-frequency components is 1.35 multiplied by the first region coefficient 0.8, which equals 1.08, or 2.08 times the low-frequency components. The first region coefficient 0.8 controls the enhancement intensity of low-frequency components in the indoor region to brighten the dark details of the crop canopy. For the boundary Gaussian sphere at the junction of corn and soil, the indoor voxel ratio is 0.3 and the outdoor voxel ratio is 0.7. The increase of the low-frequency component is 1.35 multiplied by 0.3 equals 0.405, and the increase of the high-frequency component is 1.35 multiplied by 0.7 equals 0.945. These are multiplied by the corresponding coefficients to update the spherical harmonic coefficients, so as to achieve a smooth connection between the indoor and outdoor lighting transition areas.

[0211] A tone mapping curve correction of 0.5 is applied to the global tone mapping operator. The contrast adjustment factor of 2.5 is read from the current operator and updated to 3.0. Simultaneously, the minimum radiosity value of 0.2, the maximum radiosity value of 45.0, the minimum displayable brightness of 0.05, the maximum displayable brightness of 1000, the first radiosity value of 0.8, and the second radiosity value of 40.0 remain unchanged. This ensures that dynamic range compression is adjusted only in the contrast adjustment dimension, avoiding changes to the overall brightness mapping range. Returning to the real-time rendering step, a new rendering frame for the current viewpoint is generated based on the updated spherical harmonics and the global tone mapping operator. The brightness ratio is recalculated and reduced to 1.6, which is less than the second threshold of 1.8. The adjustment is terminated, and the rendering frame is output. This achieves a balanced brightness effect where the highlights of the corn leaves are clear and the textures in the darker parts of the soil are visible, ensuring the visual continuity and realism of the virtual farmland walkthrough.

[0212] This embodiment acquires pixel-level brightness distribution data of the rendered frame in real time, calculates the proportion of overexposed and underexposed pixels and the deviation of the brightness ratio from the preset threshold, and dynamically generates exposure blending weight adjustment and tone mapping curve correction. It establishes a closed-loop feedback mechanism from brightness imbalance detection to parameter updates, solving the problems of overexposure highlight detail loss, underexposure shadow texture invisibility, and difficulty in maintaining brightness balance in real time caused by static optimization parameters failing to adapt to dynamic scene lighting changes. By applying exposure blending weight adjustment differentially according to indoor, outdoor, and boundary area types and spherical harmonic frequency bands, it enhances low-frequency components indoors to brighten shadows, enhances high-frequency components outdoors to restore highlight details, and adjusts boundary areas according to the indoor / outdoor voxel ratio. Weighted blending achieves a smooth transition, resolving the issues of adjustment oscillation and slow convergence caused by significant differences in the response characteristics of different regions to brightness adjustment. By applying the tone mapping curve correction amount to the contrast adjustment factor while keeping the radiometric range parameter unchanged, adaptive fine-tuning of the dynamic range compression curve is achieved, solving the problem of repetitive adjustments caused by the coupling of spherical harmonic coefficient adjustment and tone mapping. Finally, through multiple iterations until the brightness ratio converges, it ensures that each rendering frame can achieve a precise balance between overexposed and underexposed areas under different viewing angles and lighting conditions. This significantly improves the detail clarity of highlight areas, the texture visibility of shadow areas, and the visual consistency of indoor and outdoor light and shadow transitions in large-scene virtual roaming, effectively enhancing the user experience and scene realism.

[0213] Example 2

[0214] Please see Figure 3 Another embodiment of the present invention provides a 3DGS rendering system for large HDR scenes, comprising:

[0215] The acquisition module is used to acquire multi-exposure image sequences of the target scene and the corresponding camera pose parameters;

[0216] The illumination field reconstruction module reconstructs a high dynamic range illumination field based on the multi-exposure image sequence and the camera pose parameters to obtain illumination field configuration data. The illumination field configuration data includes at least a dynamic range compression curve, an indoor / outdoor region segmentation mask, and a material reflectivity model.

[0217] The spherical harmonic coefficient pre-optimization module is used to perform iterative optimization of the spherical harmonic coefficients of each Gaussian sphere in the three-dimensional Gaussian splash model according to the illumination field configuration data, until the difference between the area ratio of overexposed area and the area ratio of underexposed area in the current rendering frame is less than a first preset threshold, and outputs an illumination adaptation completion signal.

[0218] A real-time rendering engine is used to respond to the received lighting adaptation completion signal, to render the three-dimensional Gaussian splash model in real time at a preset rendering frame rate, and to collect pixel-level brightness distribution data of each rendering frame during the rendering process. The pixel-level brightness distribution data includes the brightness value of each pixel and an overexposure or underexposure flag determined based on a brightness threshold.

[0219] The optimization module is configured as follows:

[0220] Calculate the brightness ratio between the average brightness of the overexposed areas and the average brightness of the underexposed areas in the current rendering frame, and determine whether the brightness ratio is less than a preset second threshold.

[0221] If the brightness ratio is less than the second threshold, the current rendering frame is output based on the current spherical harmonic coefficients and the preset global tone mapping operator.

[0222] If the brightness ratio is greater than or equal to the second threshold, then the exposure fusion weight adjustment and tone mapping curve correction are calculated based on the pixel-level brightness distribution data. The exposure fusion weight adjustment is applied to the spherical harmonic coefficients of the Gaussian sphere to update the spherical harmonic coefficients, and the tone mapping curve correction is applied to the global tone mapping operator to update the global tone mapping operator. Then, the step of performing real-time rendering of the three-dimensional Gaussian splash model at a preset rendering frame rate is returned until the brightness ratio is less than the second threshold.

[0223] The embodiments of the present invention have been described above with reference to the accompanying drawings. However, the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments under the guidance of the present invention without departing from the spirit and scope of the present invention. All of these variations are within the protection scope of the present invention.

Claims

1. A 3DGS rendering method for large HDR scenes, characterized in that, include: Acquire multi-exposure image sequences of the target scene and the corresponding camera pose parameters; Based on the multi-exposure image sequence and the camera pose parameters, a high dynamic range illumination field is reconstructed to obtain illumination field configuration data. The illumination field configuration data includes at least a dynamic range compression curve, an indoor / outdoor region segmentation mask, and a material reflectivity model. Based on the lighting field configuration data, the spherical harmonic coefficients of each Gaussian sphere in the three-dimensional Gaussian splash model are iteratively optimized until the difference between the area ratio of overexposed region and the area ratio of underexposed region in the current rendering frame is less than the first preset threshold, and a lighting adaptation completion signal is output. In response to receiving the lighting adaptation completion signal, the three-dimensional Gaussian splash model is rendered in real time at a preset rendering frame rate, and pixel-level brightness distribution data of each rendering frame is collected during the rendering process. The pixel-level brightness distribution data includes the brightness value of each pixel and an overexposure or underexposure indicator determined based on a brightness threshold.

2. The 3DGS rendering method for large HDR scenes as described in claim 1, characterized in that, The 3DGS rendering method also includes: Calculate the brightness ratio between the average brightness of the overexposed areas and the average brightness of the underexposed areas in the current rendering frame, and determine whether the brightness ratio is less than a preset second threshold. If the brightness ratio is less than the second threshold, the current rendering frame is output based on the current spherical harmonic coefficients and the preset global tone mapping operator. If the brightness ratio is greater than or equal to the second threshold, then the exposure fusion weight adjustment and tone mapping curve correction are calculated based on the pixel-level brightness distribution data. The exposure fusion weight adjustment is applied to the spherical harmonic coefficients of the Gaussian sphere to update the spherical harmonic coefficients, and the tone mapping curve correction is applied to the global tone mapping operator to update the global tone mapping operator. Then, the step of performing real-time rendering of the three-dimensional Gaussian splash model at a preset rendering frame rate is returned until the brightness ratio is less than the second threshold.

3. The 3DGS rendering method for large HDR scenes as described in claim 2, characterized in that, The iterative optimization of the spherical harmonic coefficients of each Gaussian sphere in the three-dimensional Gaussian splash model includes: Based on the indoor / outdoor region segmentation mask, the Gaussian spheres in the 3D Gaussian splash model are divided into an indoor region Gaussian sphere set, an outdoor region Gaussian sphere set, and a boundary Gaussian sphere. A first optimization step size is assigned to the indoor region Gaussian sphere set, a second optimization step size is assigned to the outdoor region Gaussian sphere set, and a third optimization step size, which is between the first and second optimization step sizes, is assigned to the boundary Gaussian spheres; wherein the first optimization step size is smaller than the second optimization step size. The current viewpoint is selected sequentially from multiple preset viewpoints, and the current rendering frame of the current viewpoint is rendered based on the initial spherical harmonic coefficients obtained from initialization. The pixel radiance value of the current rendering frame is mapped to the display brightness value according to the dynamic range compression curve. Based on preset overexposure and underexposure thresholds, each pixel of the current rendering frame is marked with an overexposure or underexposure flag. The first area ratio of the overexposure region and the second area ratio of the underexposure region in the current rendering frame are calculated, and the difference between the first area ratio and the second area ratio is calculated as the current difference. Based on the diffuse reflection coefficient and specular reflection coefficient in the material reflectivity model, the spherical harmonic coefficient gradients of the Gaussian sphere in the indoor area, the Gaussian sphere in the outdoor area, and the boundary Gaussian sphere are calculated respectively, and the spherical harmonic coefficients of the corresponding Gaussian spheres are updated using the first optimization step size, the second optimization step size, and the third optimization step size respectively. Repeat the step of sequentially selecting the current viewing direction from multiple preset viewing directions until the current difference is less than the first preset threshold in three or more consecutive viewing directions. Then stop the iterative optimization and output a lighting adaptation completion signal.

4. The 3DGS rendering method for large HDR scenes as described in claim 3, characterized in that, The Gaussian spheres in the three-dimensional Gaussian splash model are divided into an indoor region Gaussian sphere set, an outdoor region Gaussian sphere set, and a boundary Gaussian sphere set, including: For each Gaussian sphere, obtain its center point coordinates and the spatial extension range defined by the covariance matrix, wherein the spatial extension range is an ellipsoidal region whose semi-axis length is determined by the eigenvalues ​​of the covariance matrix; Based on the spatial extension range, calculate the set of voxels covered by the Gaussian sphere in the indoor and outdoor region segmentation mask, where the set of voxels is all voxels that intersect with the ellipsoidal region; The number of indoor voxels and the number of outdoor voxels in the voxel set are counted, and the proportion of indoor voxels and the proportion of outdoor voxels are calculated respectively. When the proportion of indoor voxels is greater than or equal to the first division threshold, the Gaussian sphere is assigned to the indoor region Gaussian sphere set; when the proportion of outdoor voxels is greater than or equal to the second division threshold, the Gaussian sphere is assigned to the outdoor region Gaussian sphere set. When both the indoor and outdoor voxel proportions are less than the corresponding division thresholds, the Gaussian sphere is marked as a boundary Gaussian sphere. The third optimization step size is obtained by linear interpolation based on the ratio of the indoor voxel proportion to the outdoor voxel proportion.

5. The 3DGS rendering method for large HDR scenes as described in claim 4, characterized in that, Mapping the pixel radiance value of the current rendered frame to a display brightness value based on the dynamic range compression curve includes: For each pixel in the current rendering frame, all Gaussian spheres that intersect with the light rays of each pixel are queried by ray tracing or Gaussian splashing algorithm. The radiance value of the corresponding pixel is calculated based on the initial spherical harmonic coefficients of each Gaussian sphere, and a radiance image in the current view direction is generated. The radiometric image is subjected to pixel-level validity verification to detect whether there are invalid radiometric values. Invalid radiometric values ​​include NaN, infinity, or values ​​that exceed the preset physical range. If they exist, the radiometric value of the pixel is replaced with the average radiometric value of the adjacent valid pixels or the preset default radiometric value. Read the minimum radiance, maximum radiance, minimum displayable brightness, maximum displayable brightness, and contrast adjustment factor from the dynamic range compression curve; For each pixel in the radiometric image, its radiometric value is compared with the minimum and maximum radiometric values. If it is less than the minimum radiometric value, it is set to the minimum radiometric value. If it is greater than the maximum radiometric value, it is set to the maximum radiometric value. The cropped radiance value is logarithmically mapped according to the contrast adjustment factor to calculate the corresponding display brightness value, which is between the minimum and maximum displayable brightness. The mapped display brightness value is then output as the pixel brightness value of the current rendering frame.

6. The 3DGS rendering method for large HDR scenes as described in claim 5, characterized in that, The exposure fusion weight adjustment and tone mapping curve correction are calculated based on the pixel-level brightness distribution data, including: Obtain pixel-level brightness distribution data of the current rendering frame, and calculate the proportion of overexposed pixels to the total number of pixels as the first overexposure ratio, and calculate the proportion of underexposed pixels to the total number of pixels as the first underexposure ratio. The ratio of the first overexposure ratio to the first underexposure ratio is calculated as the first scaling factor; When the first overexposure ratio is greater than the first underexposure ratio, a first gain coefficient is calculated based on the first ratio factor. The first gain coefficient is equal to the first ratio factor multiplied by the first preset gain, and the first gain coefficient is used as the exposure fusion weight adjustment amount. The exposure fusion weight adjustment amount is used to enhance the weight of the high frequency component in the spherical harmonic coefficient of the Gaussian sphere corresponding to the overexposure area. When the first underexposure ratio is greater than the first overexposure ratio, a second gain coefficient is calculated based on the reciprocal of the first ratio factor. The second gain coefficient is equal to the reciprocal of the first ratio factor multiplied by a second preset gain. The second gain coefficient is used as the exposure fusion weight adjustment amount. The exposure fusion weight adjustment amount is used to enhance the weight of the low-frequency components in the spherical harmonic coefficient of the Gaussian sphere corresponding to the underexposure region. Calculate the deviation between the brightness ratio of the current rendered frame and the preset second threshold, multiply the deviation by the third preset gain to obtain the tone mapping curve correction amount, which is used to adjust the contrast adjustment factor of the dynamic range compression curve. When both the first overexposure ratio and the first underexposure ratio exceed the third ratio threshold, the exposure fusion weight adjustment amount and the tone mapping curve correction amount are calculated simultaneously.

7. The 3DGS rendering method for large HDR scenes as described in claim 6, characterized in that, Applying the exposure blending weight adjustment to the spherical harmonics of the Gaussian sphere specifically includes: For each Gaussian sphere in the set of Gaussian spheres in the indoor region, the weight of the 0th to 2nd order low-frequency components in its spherical harmonic coefficients is increased according to the exposure fusion weight adjustment amount, and the increase of the low-frequency components is the exposure fusion weight adjustment amount multiplied by the first region coefficient. For each Gaussian sphere in the outdoor area Gaussian sphere set, the weight of the 3rd order and above high-frequency components in its spherical harmonic coefficients is increased according to the exposure fusion weight adjustment amount, and the increase of the high-frequency components is the exposure fusion weight adjustment amount multiplied by the second area coefficient. For each Gaussian sphere in the set of boundary Gaussian spheres, the weights of the low-frequency and high-frequency components in its spherical harmonic coefficients are increased simultaneously according to the exposure fusion weight adjustment amount. The increase in the low-frequency component is the exposure fusion weight adjustment amount multiplied by the indoor voxel percentage of the Gaussian sphere, and the increase in the high-frequency component is the exposure fusion weight adjustment amount multiplied by the outdoor voxel percentage of the Gaussian sphere.

8. The 3DGS rendering method for large HDR scenes as described in claim 7, characterized in that, Applying the tone mapping curve correction to the global tone mapping operator specifically includes: Read the current contrast adjustment factor from the current global tone mapping operator, add the contrast adjustment factor to the tone mapping curve correction amount, and obtain the updated contrast adjustment factor; The updated contrast adjustment factor is written into the global tone mapping operator to replace the original contrast adjustment factor, while keeping the minimum radiosity, maximum radiosity, minimum displayable brightness, and maximum displayable brightness in the dynamic range compression curve unchanged. Return to the step of performing real-time rendering of the three-dimensional Gaussian splash model at a preset rendering frame rate, regenerate the current rendering frame based on the updated spherical harmonic coefficients and the updated global tone mapping operator, and recalculate the brightness ratio. The maximum number of iterations is set to Q. If the brightness ratio is less than the second threshold within the maximum number of iterations, the adjustment is terminated and the current rendering frame is output. If the brightness ratio is still not less than the second threshold after the maximum number of iterations, the current rendering frame is forcibly output and alarm information is recorded.

9. A 3DGS rendering system for HDR large scenes, used to implement the 3DGS rendering method for HDR large scenes as described in any one of claims 1-8, characterized in that, include: The acquisition module is used to acquire multi-exposure image sequences of the target scene and the corresponding camera pose parameters; The illumination field reconstruction module reconstructs a high dynamic range illumination field based on the multi-exposure image sequence and the camera pose parameters to obtain illumination field configuration data. The illumination field configuration data includes at least a dynamic range compression curve, an indoor / outdoor region segmentation mask, and a material reflectivity model. The spherical harmonic coefficient pre-optimization module is used to perform iterative optimization of the spherical harmonic coefficients of each Gaussian sphere in the three-dimensional Gaussian splash model according to the illumination field configuration data, until the difference between the area ratio of overexposed area and the area ratio of underexposed area in the current rendering frame is less than a first preset threshold, and outputs an illumination adaptation completion signal. A real-time rendering engine is used to respond to the received lighting adaptation completion signal, to render the three-dimensional Gaussian splash model in real time at a preset rendering frame rate, and to collect pixel-level brightness distribution data for each rendering frame during the rendering process. The pixel-level brightness distribution data includes the brightness value of each pixel and an overexposure or underexposure flag determined based on a brightness threshold.

10. The 3DGS rendering system for large HDR scenes as described in claim 9, characterized in that, The 3DGS rendering system also includes: The optimization module is configured as follows: Calculate the brightness ratio between the average brightness of the overexposed areas and the average brightness of the underexposed areas in the current rendering frame, and determine whether the brightness ratio is less than a preset second threshold. If the brightness ratio is less than the second threshold, the current rendering frame is output based on the current spherical harmonic coefficients and the preset global tone mapping operator. If the brightness ratio is greater than or equal to the second threshold, then the exposure fusion weight adjustment and tone mapping curve correction are calculated based on the pixel-level brightness distribution data. The exposure fusion weight adjustment is applied to the spherical harmonic coefficients of the Gaussian sphere to update the spherical harmonic coefficients, and the tone mapping curve correction is applied to the global tone mapping operator to update the global tone mapping operator. Then, the step of performing real-time rendering of the three-dimensional Gaussian splash model at a preset rendering frame rate is returned until the brightness ratio is less than the second threshold.

Images