A method, medium, and system for topographic mapping based on the fusion of laser point clouds and optical images.

Abstract

This invention provides a topographic mapping method, medium, and system based on the fusion of laser point clouds and optical images, belonging to the field of surveying and mapping technology. This invention performs joint radiometric calibration of the reflection intensity of optical images and point clouds, outputting a radiometrically aligned image and a radiometrically calibrated point cloud. The above two types of data are input into a geometry-guided cross-modal dynamic graph attention fusion artificial intelligence model. Utilizing a differentiable K-nearest neighbor dynamic graph reconstruction module and an asymmetric bidirectional cross-attention mechanism, the geometric topological relationship of the point cloud actively guides image feature extraction, outputting dense point cloud semantic labels, surface normal vector estimation, occlusion confidence map, and texture quality assessment map. Based on the occlusion confidence map and texture quality assessment map, a dynamic task queue and work-stealing scheduling algorithm are used to generate topographic mapping results in parallel. This solves the problem of cross-modal feature alignment failure caused by the failure of geometric topological relationship to actively guide image feature extraction during cross-modal fusion of laser point clouds and optical images.

Description

Technical Field

[0001] This invention belongs to the field of surveying and mapping technology, and specifically relates to a topographic and geomorphological surveying method, medium, and system based on the fusion of laser point clouds and optical images. Background Technology

[0002] Topographic mapping is a core technology for acquiring three-dimensional morphological and semantic information of the Earth's surface, and it is widely used in fields such as land surveying, disaster assessment, and engineering site selection. The current mainstream approach is to fuse point clouds acquired by airborne lidar with synchronous optical imagery, leveraging the precise geometric depth information of the point clouds and the rich texture information of the imagery to achieve complementarity. In traditional fusion methods, spatial domain filtering is typically performed on the point clouds for noise reduction, followed by matching manually designed geometric descriptors with image features, and finally, data fusion is completed through interpolation or voxel projection.

[0003] However, traditional fusion methods have significant drawbacks when dealing with complex terrain. Spatial domain filtering struggles to distinguish between real terrain details and random noise on irregular point cloud manifolds, leading to excessive smoothing of terrain boundary details. More critically, manually designed geometric descriptors suffer from metric failure in cross-modal matching; that is, the 3D geometric topology of the point cloud cannot actively guide the extraction direction of texture features from optical images. The feature extraction processes of the two modalities are independent of each other, causing feature alignment errors in the common embedding space to accumulate with increasing terrain complexity.

[0004] In current laser point cloud and optical image fusion, due to the significant differences in the radiation physics models of the two modalities and the fact that the geometric topology of the point cloud is not involved in the image feature extraction process, the cross-modal feature alignment mechanism loses its geometric constraints when terrain occlusion is severe or texture degradation occurs, leading to systematic deviations in semantic annotation and normal vector estimation. In other words, existing technologies suffer from a technical problem where cross-modal feature alignment fails to function properly due to the lack of proactive guidance from geometric topology for image feature extraction during laser point cloud and optical image fusion. Summary of the Invention

[0005] In view of this, the present invention provides a topographic mapping method, medium and system based on the fusion of laser point clouds and optical images, which can solve the technical problem in the prior art that the cross-modal feature alignment fails due to the failure of geometric topological relationships to actively guide image feature extraction during cross-modal fusion of laser point clouds and optical images.

[0006] The present invention is implemented as follows: The first aspect of the present invention provides a topographic mapping method based on the fusion of laser point clouds and optical images, comprising the following steps:

[0007] The original point cloud acquired by lidar is subjected to point cloud noise separation and terrain frequency analysis based on spectral signal processing, and the denoised terrain point cloud is output.

[0008] The hollow regions in the denoised terrain point cloud are filled by an adaptive point cloud hole completion algorithm based on hydrodynamic potential flow theory, and the complete terrain point cloud is output.

[0009] Flat field correction and atmospheric radiative correction are performed on optical images. The incident angle is normalized and distance compensation is calibrated for the reflection intensity of complete terrain point clouds. A joint radiative transfer model is established, and radiative aligned images and radiatively calibrated point clouds are output.

[0010] The radiometrically calibrated point cloud and the radiometrically aligned image are input into the geometry-guided cross-modal dynamic graph attention fusion model. The geometry-guided cross-modal dynamic graph attention fusion model includes a point cloud branch and an image branch. The point cloud branch extracts geometric features using a differentiable K-nearest neighbor dynamic graph reconstruction module. The two branches interact across modalities through an asymmetric bidirectional cross attention module. The point cloud feature vector is used as a query to guide image feature aggregation, and the point cloud normal vector field is used as the geometric position encoding of the image branch. The output includes dense point cloud semantic labels, surface normal vector estimation, occlusion confidence map, and texture quality assessment map.

[0011] Based on the occlusion confidence map and texture quality assessment map, a dynamic task queue and work-stealing scheduling algorithm are used to distribute tasks to multiple CUDA streams and generate terrain mapping results in parallel.

[0012] Specifically, the point cloud noise separation and terrain frequency analysis based on spectral signal processing involves constructing a K-nearest neighbor weighted graph with each point in the original point cloud as a node. The edge weights are jointly determined by the Euclidean distance and the cosine of the angle between the normal vectors. The eigenvalue decomposition of the normalized graph Laplacian matrix is ​​calculated, the point cloud elevation signal is projected onto the spectral domain, and the spectral domain filter is adaptively designed according to the terrain fractal dimension. Compressed sensing reconstruction is performed through the sparse prior of the graph signal.

[0013] The cutoff frequency of the spectral domain filter is modulated by the power spectrum of the optical image texture. Low-frequency coefficients corresponding to the macroscopic undulations of the terrain are preserved, mid-frequency coefficients corresponding to the boundaries of ground features are selectively preserved, and high-frequency coefficients corresponding to random noise are attenuated. The computational core of the compressed sensing reconstruction is approximately accelerated by randomized singular value decomposition.

[0014] Specifically, the point cloud cavity adaptive completion algorithm based on the fluid dynamics potential flow theory treats the cavity boundary points of the denoised terrain point cloud as solid wall boundary conditions and the cavity region as a fluid domain. The Laplace equation is numerically solved using the finite difference method on an adaptive grid. The grid resolution is adaptively controlled by the curvature of the boundary points. Gravitational potential correction terms are superimposed along the slope direction in the slope region. The texture gradient of the radiation-aligned image is injected into the Laplace equation as an anisotropic diffusion correction term.

[0015] The numerical solution of the Laplace equation uses the multigrid preconditional conjugate gradient method to accelerate iterative convergence. The completion points are generated along the gradient direction of the potential field, and the density of the completion points is adaptively matched with the density of the surrounding denoised terrain point cloud.

[0016] The joint radiative transfer model aligns the two modal radiative distributions in a common semantic feature space by using contrastive learning loss to align the optical image radiative values ​​after flat field correction and atmospheric radiative correction with the complete terrain point cloud reflection intensity after incident angle normalization and distance compensation calibration.

[0017] In the geometrically guided cross-modal dynamic graph attention fusion model, the point cloud branch adopts an improved PointNet++ backbone network, and the image branch adopts a Swing Transformer backbone to extract multi-scale hierarchical features and introduces deformable convolutional modules at each scale. The features of the two branches interact across three scale levels, and an asymmetric bidirectional cross attention module is designed at each level.

[0018] The geometrically guided cross-modal dynamic graph attention fusion model introduces a geometric consistency gating mechanism to calculate the cosine similarity of the features of the two branches. When the cosine similarity is lower than the geometric consistency gating threshold, an internal iterative refinement loop is triggered. The internal iterative refinement loop is an embedded lightweight iterative registration sub-network that updates the projection matrix of Query and Key in each iteration.

[0019] The geometry-guided cross-modal dynamic graph attention fusion model designs cross-modal jump paths, directly injecting the original geometric features of the bottom layer of the point cloud branch into the high-level semantic decision layer; cross-layer jump connections within the same branch are retained simultaneously, and the feature fusion weights of the two types of jump paths are adaptively allocated by hierarchical attention gating.

[0020] The output layer of the geometry-guided cross-modal dynamic graph attention fusion model is jointly trained through a multi-task loss function, which consists of a weighted sum of four terms: semantic cross-entropy loss, normal vector cosine loss, confidence binary cross-entropy loss, and texture quality mean square error loss.

[0021] The geometry-guided cross-modal dynamic graph attention fusion model supports streaming block inference and adopts an overlapping sliding window strategy. During training, a course learning strategy is adopted, gradually transitioning from flat terrain samples to mountain samples, and using gradient norm stability as a criterion to determine the course advancement nodes.

[0022] The relationship between the mean confidence score of the occlusion confidence map and the batch size is determined by the modal fusion quality evaluation function. The modal fusion quality evaluation value is calculated by weighted summation of the mean of the occlusion confidence map, the mean of the texture quality evaluation map, and the mean of the cosine similarity of the two-branch features. The geometric consistency gating threshold, batch size, and number of CUDA streams are dynamically adjusted based on the modal fusion quality evaluation value.

[0023] The training dataset of the geometry-guided cross-modal dynamic graph attention fusion model covers four types of terrain: plains, hills, mountains, and canyons. The labeled categories include ground, vegetation, buildings, water bodies, and bare rocks. Data augmentation is performed through random rotation, point cloud density downsampling, and image brightness perturbation.

[0024] Specifically, the dynamic task queue and work-stealing scheduling algorithm pre-estimates the computational cost of each task block based on the local density of the complete terrain point cloud, subdivides the task into subtasks with approximately equal computational cost, maintains an independent task queue for each CUDA stream, and steals subtasks from the tail of the queue of the most heavily loaded CUDA stream. This, combined with the asynchronous execution mechanism of the CUDA pipeline, allows data transmission and computation to overlap.

[0025] Specifically, the K-value of the K-nearest neighbor weighted graph ranges from 10 to 20; the weight ratio of the Euclidean distance to the cosine of the angle between the normal vector and the weights ranges from 6:4 to 7:3; the number of retained feature vectors k ranges from 50 to 200; the K-value of the differentiable K-nearest neighbor dynamic graph reconstruction module ranges from 8 to 32; the geometric consistency gating threshold ranges from 0.5 to 0.8; the ratio of the sliding window step size to the window size in the overlapping sliding window strategy ranges from 0.5 to 0.75; the grid spacing in high curvature regions ranges from 0.02 to 0.05 m, and the grid spacing in low curvature regions ranges from 0.1 to 0.3 m; and the terrain erosion coefficient ranges from 0.01 to 0.05. .

[0026] A second aspect of the present invention provides a computer-readable storage medium storing program instructions, which, when executed in a computer, are used to perform the above-described method for topographic mapping based on the fusion of laser point clouds and optical images.

[0027] A third aspect of the present invention provides a topographic mapping system based on the fusion of laser point clouds and optical images, comprising the aforementioned computer-readable storage medium, wherein the system is a computer, the computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor for executing program instructions stored in the computer-readable storage medium.

[0028] This invention solves the technical problem of cross-modal feature alignment failure by constructing a geometrically guided cross-modal dynamic graph attention fusion model and embedding a differentiable K-nearest neighbor dynamic graph reconstruction module into an asymmetric bidirectional cross-attention framework. This enables the geometric topological relationship of the radiometric calibration point cloud to actively guide the extraction of texture features from the radiometrically aligned image, thereby solving the technical problem of cross-modal feature alignment failure.

[0029] In traditional methods, the two modalities are extracted independently, and geometric information does not constrain image branches, leading to feature alignment failure under occlusion or texture degradation conditions. This invention uses the point cloud normal vector field as the geometric position encoding of image branches and uses the point cloud feature vector as a query to perform spatial attention queries on image features. This constrains the aggregation direction of image texture features to the geometric topology of the point cloud, fundamentally eliminating the alignment error caused by the independent extraction of the two modalities. A geometric consistency gating mechanism triggers an internal iterative refinement loop when the semantics of the two branch features are inconsistent, ensuring robustness of feature alignment in the cross-modal common embedding space.

[0030] In summary, the present invention solves the technical problem mentioned in the background art of cross-modal feature alignment failure caused by the failure of geometric topological relationships to actively guide image feature extraction during cross-modal fusion of laser point clouds and optical images. Attached Figure Description

[0031] Figure 1 This is a flowchart of the method of the present invention.

[0032] Figure 2 The structural diagram of a geometry-guided cross-modal dynamic graph attention fusion model.

[0033] Figure 3 This is a spatial distribution diagram of the cutoff frequency of the spectral domain filter as a function of the texture power spectrum.

[0034] Figure 4 Comparison of tangential continuity before and after filling in the voids in the point cloud.

[0035] Figure 5 The cosine similarity distribution of the two modal features before and after alignment of the joint radiative transfer model is shown in the figure.

[0036] Figure 6 Confidence maps for occlusion in different terrain zones.

[0037] Figure 7Output a semantic label map for dense point clouds.

[0038] Figure 8 Output the texture quality evaluation map.

[0039] Figure 9 A spatial distribution map for estimating the surface normal vectors.

[0040] Figure 10 The graph shows the spectral response curves before and after radiation alignment for each band. Detailed Implementation

[0041] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0042] like Figure 1 The diagram shown is a flowchart of a topographic mapping method based on the fusion of laser point clouds and optical images, provided by the first aspect of this invention. This method includes the following steps:

[0043] S01. Perform point cloud noise separation and terrain frequency analysis based on spectral signal processing on the original point cloud acquired by lidar, and output denoised terrain point cloud;

[0044] S02. The hollow areas in the denoised terrain point cloud are filled by the point cloud hole adaptive completion algorithm based on the hydrodynamic potential flow theory, and the complete terrain point cloud is output.

[0045] S03. Perform flat field correction and atmospheric radiation correction on the optical image, normalize the incident angle and calibrate the distance compensation for the reflection intensity of the complete terrain point cloud, establish a joint radiative transfer model, and output the radiation-aligned image and radiation-calibrated point cloud.

[0046] S04. Input the radiometrically calibrated point cloud and the radiometrically aligned image into the geometry-guided cross-modal dynamic graph attention fusion model, and output dense point cloud semantic labels, surface normal vector estimation, occlusion confidence map and texture quality assessment map.

[0047] S05. Based on the occlusion confidence map and texture quality assessment map, the dynamic task queue and work-stealing scheduling algorithm are used to distribute the tasks to multiple CUDA streams and generate terrain mapping results in parallel.

[0048] The specific implementation steps of point cloud noise separation and terrain frequency analysis based on spectral signal processing are as follows: A K-nearest neighbor weighted graph is constructed with each point in the original point cloud as a node, where the K value ranges from 10 to 20; the edge weights are determined by a combined weighting of the Euclidean distance and the cosine of the angle between the normal vectors, and the ratio of these weights is determined based on reconstruction accuracy after multiple rounds of comparative experiments on four types of terrain samples (plains, hills, mountains, and canyons), with a typical range of 6:4 to 7:3; the normalized graph Laplacian matrix is ​​calculated. eigenvalue decomposition ,in For the graph Fourier basis matrix, The matrix is ​​a diagonal eigenvalue matrix, with all units being dimensionless. Point cloud elevation signals are projected into the spectral domain, and spectral filters are adaptively designed based on the fractal dimension of the terrain. Low-frequency coefficients, corresponding to macroscopic terrain undulations, are preserved; mid-frequency coefficients, corresponding to feature boundaries, are selectively preserved; and high-frequency coefficients, corresponding to random noise, are attenuated. The cutoff frequency is modulated by the power spectrum of the optical image texture, with higher cutoff frequencies in texture-rich areas and lower cutoff frequencies in uniform areas. Finally, compressed sensing reconstruction is performed using sparse priors of the image signal. The computational core is approximately accelerated by randomized singular value decomposition, resulting in an overall computational complexity of O(log n). ,in The number of point cloud nodes. To preserve the number of feature vectors, The range is 50 to 200, determined by experimental verification after achieving a balance between terrain reconstruction accuracy and computational cost.

[0049] The proposed point cloud noise separation and terrain frequency analysis algorithm based on spectral signal processing transfers spectral analysis theory from signal processing to irregular point cloud manifolds, enabling frequency domain analysis on the point cloud structure. Traditional point cloud denoising methods process points one by one in the spatial domain, making it difficult to distinguish between real terrain details and random noise. This algorithm constructs a graph Laplacian operator to decompose the elevation signal into components with different oscillation frequencies on the graph topology, allowing macroscopic terrain undulations, feature boundaries, and random noise to be naturally separated in the spectral domain. The denoising process adaptively adjusts the retained frequency bands according to the complexity of the feature texture, avoiding excessive smoothing of terrain details. Simultaneously, under sparse point cloud conditions, it reconstructs complete terrain through compressed sensing, breaking through the dependence of traditional spatial domain filtering on point cloud density. This provides frequency-consistent and geometrically complete denoised terrain point clouds for subsequent hole completion and cross-modal fusion.

[0050] The specific implementation steps of the point cloud cavity adaptive completion algorithm based on fluid dynamics potential flow theory are as follows: the cavity boundary points of the denoised terrain point cloud are regarded as solid wall boundary conditions, and the cavity region is regarded as a fluid domain; the Laplace equation is numerically solved using the finite difference method on the adaptive mesh. ,in The velocity potential function is expressed in units of 1000 kJ / m². The mesh resolution is adaptively controlled by the curvature of the boundary points, with the mesh spacing ranging from 0.02 to 0.05 in high curvature regions. The grid spacing in the low curvature region ranges from 0.1 to 0.3. The curvature threshold was determined through statistical analysis of various terrain samples obtained from field measurements. A gravitational potential correction term was superimposed along the slope direction in the slope region; the correction term coefficient was determined by both the slope angle and the terrain erosion coefficient, which ranged from 0.01 to 0.05. The method was determined through iterative experiments comparing measured field terrain data and numerical simulations. The texture gradient of the radiatively aligned image was injected into the Laplace equation as an anisotropic diffusion correction term, guiding the completion points to retain more detail in directions with rich texture variations. Completion point positions were generated along the potential field gradient direction to ensure tangential continuity between the completed surface and boundary points. The iterative convergence was accelerated using a multi-grid preconditioned conjugate gradient method, achieving a computational complexity of O(n log n). The density of the completed points is adaptively matched with the density of the surrounding known denoised terrain point cloud to output a complete terrain point cloud.

[0051] The adaptive point cloud cavity completion algorithm based on hydrodynamic potential flow theory introduces the smoothness and boundary adaptability of the potential flow field from continuum mechanics into the discrete point cloud completion problem. The potential flow field naturally satisfies the harmonicity described by the Laplace equation within the cavity, ensuring global geometric smoothness and tangential continuity at the boundary of the completed points, avoiding geometric distortions caused by interpolation methods in large cavity regions. A gravity potential correction term ensures that the completed results conform to the physical laws of gravity erosion landforms in slope regions, reducing topographically unreasonable completion shapes. The introduction of the radiometrically aligned image texture gradient as an anisotropic diffusion correction term allows the completion process to proceed collaboratively under geometric and radiometric texture constraints, ensuring geometric and aesthetic consistency of the completed region and providing a physically reasonable and geometrically continuous complete terrain point cloud for subsequent cross-modal fusion.

[0052] The flat-field correction is a process of calibrating and compensating for spatial radiation inhomogeneities caused by lens vignetting and CCD non-uniform response in optical cameras. Flat-field images acquired under uniform diffuse light sources are used to extract the gain coefficient for each pixel, and the acquired images are normalized pixel-by-pixel. Atmospheric radiation correction is a process of removing the influence of atmospheric scattering and absorption on image radiance based on an atmospheric radiative transfer model. Atmospheric optical path parameters are obtained through synchronously measured atmospheric parameters or estimated based on the dark pixel method. Incident angle normalization corrects the laser return intensity of complete terrain point clouds according to the cosine effect of the incident angle, ensuring that the intensity value reflects the true reflectivity of the target. Distance compensation calibration restores the laser return intensity according to the squared distance attenuation law; the calibration coefficients are obtained by fitting multi-distance measured data collected on a standard target with known reflectivity. The joint radiative transfer model is a model constructed by aligning the two-modal radiative distributions in a common semantic feature space using contrastive learning loss, with the radiative values ​​of the optical image after flat field correction and atmospheric radiative correction and the reflection intensity of the complete terrain point cloud after incident angle normalization and distance compensation calibration. The output results are denoted as radiative aligned image and radiative calibration point cloud, respectively.

[0053] The geometry-guided cross-modal dynamic graph attention fusion model is based on deep learning. Its specific structure is as follows: the model input includes two parallel branches. The first branch receives the 3D coordinates and normal vectors of the radiometrically calibrated point cloud, and the second branch receives multispectral image patches from the radiometrically aligned image. The first branch uses an improved PointNet++ backbone network, introducing a differentiable K-nearest neighbor dynamic graph reconstruction module at each level. The edge weights are determined by a weighted average of three geometric features: point distance, normal angle, and curvature difference. The weight coefficients are obtained through end-to-end learning on multi-scene training sets, with K ranging from 8 to 32. The K value is determined iteratively using feature consistency loss as an indicator through ablation experiments on terrain samples with different point cloud densities. The second branch uses a Swin Transformer backbone to extract multi-scale hierarchical features and introduces deformable convolution modules at each scale to adapt to feature space distortion caused by terrain deformation. The two branches interact across three scale levels, with each level employing an asymmetric bidirectional cross-attention module: the first branch point cloud feature vector serves as the query for spatial attention queries on the second branch image features; the second branch image texture gradient vector serves as the auxiliary key to guide the aggregation direction of the first branch point cloud local features; and the value is provided by the second branch image features. The second branch image features use the first branch point cloud normal field as geometric position encoding to enhance spatial perception. A geometric consistency gating mechanism is introduced: the cosine similarity between the two branches is calculated, and when the cosine similarity is lower than the geometric consistency gating threshold, an internal iterative refinement loop is triggered, iterating a maximum of three times until convergence. The geometric consistency gating threshold ranges from 0.5 to 0.7, determined through grid search experiments on the validation set targeting the average intersection-union ratio of semantic labels. The internal iterative refinement loop is an embedded lightweight iterative registration sub-network, updating the projection matrices of the query and key in each iteration. The network design employs cross-modal skipping paths, directly injecting the raw geometric features from the lower layers of the first branch into the higher-level semantic decision layer. This prevents fine-grained geometric information from being overly abstracted and diluted in deep networks. Within the same branch, cross-layer skipping connections are simultaneously preserved, and the feature fusion weights for both types of skipping paths are adaptively allocated by hierarchical attention gating. The output layer simultaneously generates dense point cloud semantic labels, surface normal vector estimates, occlusion confidence maps, and texture quality assessments. Figure 4The network employs a joint output, trained using a multi-task loss function. This function consists of a weighted sum of four terms: semantic cross-entropy loss, normal vector cosine loss, confidence binary cross-entropy loss, and texture quality mean square error loss. The weights of each term are determined iteratively on the validation set through multiple rounds of ablation experiments. The network supports streaming block inference and uses an overlapping sliding window strategy. The ratio of the sliding window step size to the window size ranges from 0.5 to 0.75, determined experimentally based on the minimum overlap rate required for seamless stitching. During training, a course-learning strategy is used, gradually transitioning from flat terrain samples to mountainous samples, with gradient norm stability determining the course progression nodes.

[0054] In the geometrically guided cross-modal dynamic graph attention fusion model, the connection weight allocation between neurons is jointly determined by the graph edge weights of the differentiable K-nearest neighbor dynamic graph reconstruction module and the attention weights of the asymmetric bidirectional cross-attention module; the feature transfer allocation between layers is adjusted by the cross-modal jump path weight gating; the batch allocation of data loops dynamically adjusts the batch size based on the mean confidence of the occlusion confidence map; the memory allocation is pre-estimated and segmented according to the product of the feature map size and the number of channels at each scale level; the memory allocation pre-allocates buffers according to the number of blocks in the complete terrain point cloud and the size of the sliding window overlap area; the CUDA stream allocation allocates the forward computation of the first and second branches to independent CUDA streams for asynchronous execution; the CUDA thread block allocation dynamically divides the thread block size after estimating the number of points per block based on the local density of the complete terrain point cloud; the CUDA shared memory allocation is allocated according to the neighborhood point coordinate cache requirements of the K-nearest neighbor dynamic graph reconstruction module; and the hierarchical allocation assigns the lower-level feature extraction layer and the higher-level semantic fusion layer to different computation priority queues.

[0055] The relationship between the mean confidence score of the occlusion confidence map and the batch size is determined by a modal fusion quality evaluation function, which uses the mean confidence score of the occlusion confidence map as the basis for evaluation. (Dimensionless) Mean of Texture Quality Assessment Map (Dimensionless) Mean cosine similarity of features between two branches (Dimensionless) Calculation of fusion quality assessment value (Dimensionless), the calculation formula is as follows ;in , , The three weights are dimensionless, with a sum of 1, and each coefficient ranging from 0.2 to 0.5. They are determined iteratively using the average intersection-union ratio of semantic labels as the evaluation metric, through comparative experiments on the validation set with different weight combinations. At this time, the geometric consistency gating threshold is lowered to 0.5-0.6, and the batch size is taken as the largest batch size. ;when At this time, the geometric consistency gating threshold is maintained between 0.6 and 0.7, and the batch size is set to... ;when At that time, the geometric consistency gating threshold is increased to 0.7–0.75, and the batch size is taken as... The number of CUDA streams increased from 2 to 4; when At that time, the geometric consistency gating threshold was increased to 0.75–0.8, and the batch size was set to [value missing]. The number of CUDA streams increased from 2 to 8; (Dimensionless, unit is individual) is determined by the ratio of video memory capacity to video memory usage per sample.

[0056] The steps for establishing the training dataset for the geometry-guided cross-modal dynamic graph attention fusion model specifically include: collecting raw point clouds and synchronous optical multispectral images from airborne lidar covering four types of terrain: plains, hills, mountains, and canyons, with a point cloud density ranging from 10 to 100 points / ... The ground resolution of the imagery ranges from 0.05 to 0.3. The original point cloud was manually semantically annotated, with annotation categories including ground, vegetation, buildings, water bodies, and bare rock; the optical multispectral image was radiometrically calibrated to generate texture quality assessment maps and occlusion confidence maps corresponding to the original point cloud; data augmentation was performed through random rotation, point cloud density downsampling, and image brightness perturbation, with an augmentation factor of 5 to 10 times; the dataset was divided into training, validation, and test sets in a 7:2:1 ratio.

[0057] The specific steps for training the geometry-guided cross-modal dynamic graph attention fusion model include: employing a course learning strategy; in the first stage, training is conducted using only flat terrain samples, with a gradient norm below 0.1 as the course progression criterion; in the second stage, hill and building samples are added; and in the third stage, mountain and canyon samples are added. The optimizer used is AdamW, with an initial learning rate range of [missing information]. The weight decay coefficient range is All parameters were determined through grid search experiments on the validation set; the learning rate adopted a cosine annealing strategy, with a total of 80-120 training rounds; the geometric consistency gating threshold was determined through grid search experiments on the validation set with the average intersection-union ratio of semantic labels as the target; the weights of each item in the multi-task loss function were determined iteratively on the validation set through multiple rounds of ablation experiments to determine the optimal combination.

[0058] The technical effects of the geometrically guided cross-modal dynamic graph attention fusion model are as follows: By embedding the differentiable K-nearest neighbor dynamic graph reconstruction module into an asymmetric bidirectional cross-attention framework, the model actively guides the extraction of texture features from the radiometrically aligned image by the geometric topological relationship of the radiometrically calibrated point cloud, overcoming the problem of descriptor metric failure in cross-modal matching caused by traditional handcrafted features. The geometric consistency gating mechanism triggers an internal iterative refinement loop when the semantics of the two branches are inconsistent, ensuring the robustness of feature alignment in the cross-modal common embedding space and fundamentally suppressing mismatches caused by differences in the radiometric physical model. The cross-modal jump path directly injects the original geometric information from the lower layer of the first branch into the higher-level decision, preserving the contribution of the fine-grained terrain structure in deep fusion. The multi-task loss function introduces Laplacian surface smoothness constraints and radiometric consistency constraints into the training process, constraining the fusion results in both geometric rationality and stability across illumination conditions, thus improving the physical reliability and semantic consistency of topographic mapping results.

[0059] The specific implementation of the dynamic task queue and work-stealing scheduling algorithm is as follows: the computational cost of each task block is estimated in advance based on the local density of the complete terrain point cloud, and the task is subdivided into subtasks with approximately equal computational cost; each CUDA stream maintains an independent task queue, and the idle CUDA stream steals subtasks from the tail of the queue of the most heavily loaded CUDA stream; the asynchronous execution mechanism of the CUDA pipeline is used to make data transmission and computation overlap; the number of concurrent CUDA streams ranges from 2 to 8, constrained by the memory capacity and task block size, and the optimal value is selected after experimentally measuring the GPU utilization under each number of concurrent streams.

[0060] The terrain fractal dimension is a dimensionless parameter describing the self-similarity of the terrain surface, estimated by the variance ratio of the point cloud elevation signal at different scales, and used to adaptively adjust the cutoff frequency of the spectral domain filter. Compressed sensing reconstruction utilizes the sparsity of the signal in the graph Fourier domain to recover the complete signal from undersampled observations through convex optimization. The Nyquist sampling constraint refers to the constraint in traditional signal reconstruction that the sampling frequency must be no less than twice the highest frequency of the signal. The multigrid preconditioned conjugate gradient method is a numerical method that alternately solves linear equations on multiple resolution grids to accelerate iterative convergence, used to accelerate the numerical solution of the potential flow equation. The anisotropic diffusion correction term introduces a direction-dependent diffusion coefficient into the Laplace equation, with the diffusion intensity determined by the gradient magnitude of the radiometrically aligned image texture, reducing the diffusion intensity at texture boundaries to preserve the geometric details of feature boundaries. The randomized singular value decomposition is a method that approximates the low-rank singular value decomposition of the matrix through random projection, used to accelerate the core computational step of the normalized graph Laplace matrix eigenvalue decomposition. The differentiable K-nearest neighbor dynamic graph reconstruction module refers to a module where the gradient of the graph edge weights with respect to the input geometric features exists and can be backpropagated, allowing the graph structure to be optimized end-to-end during network training. The asymmetric bidirectional cross-attention module refers to a module where the two branches have asymmetrical roles in attention calculation; the first branch provides the query based on point cloud features, while the second branch provides the key and value based on image features. These roles do not interchange, unlike the symmetric attention mechanism. The curriculum learning strategy refers to a training strategy that gradually increases the difficulty of training samples from easy to difficult during training to alleviate the gradient instability problem in the early stages of training. The cross-modal jump path refers to a cross-branch feature transfer path that directly connects the original geometric features at the bottom layer of the first branch to the high-level semantic decision layer, distinguishing it from cross-layer jump connections within the same branch. The geometric consistency gating mechanism refers to a gating structure that determines whether to trigger the internal iterative refinement loop based on the cosine similarity of the features of the two branches, used to adaptively adjust the refinement depth of cross-modal feature alignment.

[0061] The specific implementation of step S01 is as follows: Using each point in the original point cloud as a node, a K-nearest neighbor weighted graph is constructed based on Euclidean distance, with K ranging from 10 to 20. Edge weights are determined by a joint weighting of the Euclidean distance and the cosine of the angle between the normal vectors. The typical ratio of the Euclidean distance weight to the cosine of the angle between the normal vectors is 6:4 to 7:3. This ratio was determined based on reconstruction accuracy analysis after multiple rounds of comparative experiments on four types of terrain samples: plains, hills, mountains, and canyons. Subsequently, the normalized graph Laplacian matrix is ​​calculated. eigenvalue decomposition ,in For the graph Fourier basis matrix, The matrix is ​​a diagonal eigenvalue matrix. After projecting the point cloud elevation signal into the spectral domain, a spectral domain filter is adaptively designed based on the terrain fractal dimension, which is estimated by the variance ratio of the point cloud elevation signal at different scales. Low-frequency coefficients, corresponding to macroscopic terrain undulations, are preserved; mid-frequency coefficients, corresponding to feature boundaries, are selectively preserved; and high-frequency coefficients, corresponding to random noise, are attenuated. The cutoff frequency is modulated by the optical image texture power spectrum, with higher cutoff frequencies in textured areas and lower cutoff frequencies in uniform areas. Finally, compressed sensing reconstruction is performed using sparse priors of the map signal. The computational core is approximately accelerated by randomized singular value decomposition, and the overall computational complexity is O(log n). ,in The number of point cloud nodes. To preserve the number of feature vectors, k ranges from 50 to 200. The purpose of this step is to perform frequency domain analysis on irregular point cloud manifolds, so that the macroscopic topographic relief, feature boundaries and random noise are naturally separated in the spectral domain, and the output is a denoised topographic point cloud with consistent frequency and geometric integrity, providing high-quality input for subsequent steps.

[0062] The specific implementation of step S02 is as follows: the void boundary points of the denoised terrain point cloud are regarded as solid wall boundary conditions, the void region is regarded as a fluid domain, and the Laplace equation is numerically solved using the finite difference method on the adaptive grid. ,in The velocity potential function is used. The grid resolution is adaptively controlled by the curvature of the boundary points. The grid spacing ranges from 0.02 to 0.05 m in high-curvature regions and from 0.1 to 0.3 m in low-curvature regions. The curvature threshold is determined through statistical analysis of various terrain samples obtained from field measurements. A gravitational potential correction term is superimposed along the slope direction in the slope region. The correction term coefficient is determined by both the slope angle and the topographic erosion coefficient, which ranges from 0.01 to 0.05. The texture gradient of the radiatively aligned image is injected into the Laplacian equation as an anisotropic diffusion correction term, guiding the completion points to retain more detail in directions with rich texture variations. Completion point positions are generated along the potential gradient direction to ensure tangential continuity between the completed surface and boundary points. The numerical solution of the Laplacian equation employs a multigrid preconditioned conjugate gradient method to accelerate iterative convergence, achieving a computational complexity of O(log n). The density of the completed points is adaptively matched with the density of the surrounding known denoised terrain point cloud. The purpose of this step is to collaboratively complete the point cloud holes under geometric and radial texture constraints, ensuring the consistency of the geometry and appearance of the completed area, and outputting a physically reasonable and geometrically continuous complete terrain point cloud.

[0063] The specific implementation of step S03 is as follows: Flat-field correction uses flat-field images acquired under uniform diffuse light sources. Gain coefficients for each pixel are extracted, and the acquired image is normalized pixel-by-pixel to eliminate spatial radiation inhomogeneity caused by lens vignetting and CCD non-uniform response. Atmospheric radiation correction uses an atmospheric radiative transfer model to remove the influence of atmospheric scattering and absorption on image radiance. Atmospheric optical path parameters are obtained through synchronously measured atmospheric parameters or estimated based on the dark pixel method. The laser return intensity of the complete terrain point cloud is normalized according to the incident angle cosine effect and calibrated for distance compensation according to the squared distance attenuation law. The calibration coefficients are obtained by fitting multi-distance measured data collected on a standard target with known reflectivity. The joint radiative transfer model aligns the radiative values ​​of the processed optical image and the reflection intensity of the complete terrain point cloud in a common semantic feature space using a contrastive learning loss to align the two modal radiative distributions. The output results are denoted as the radiative aligned image and the radiative calibrated point cloud, respectively. The purpose of this step is to eliminate the interference of differences in the radiative physics model on subsequent cross-modal fusion.

[0064] The specific implementation of step S04 is as follows: The geometry-guided cross-modal dynamic graph attention fusion model comprises two parallel branches. The point cloud branch uses an improved PointNet++ backbone network, introducing a differentiable K-nearest neighbor dynamic graph reconstruction module at each level. The edge weights are determined by a weighted average of three geometric features: point distance, normal angle, and curvature difference, with K ranging from 8 to 32. The image branch uses a Swin Transformer backbone to extract multi-scale hierarchical features and introduces deformable convolution modules at each scale to adapt to feature space distortion caused by ground deformation. The features of the two branches interact across three scale levels. Each level designs an asymmetric bidirectional cross-attention module: the point cloud feature vector serves as the query for spatial attention queries on image features, the image texture gradient vector serves as the auxiliary key to guide the aggregation direction of local point cloud features, and the value is provided by the image features; the image branch uses the point cloud normal vector field as a geometric position encoding to enhance spatial perception capabilities. The geometric consistency gating mechanism calculates the cosine similarity of features between two branches. When the similarity falls below the geometric consistency gating threshold, an internal iterative refinement loop is triggered, iterating a maximum of 3 times. The geometric consistency gating threshold ranges from 0.5 to 0.8. A cross-modal jump path directly injects the original geometric features of the lower-level point cloud branches into the high-level semantic decision layer. The output layer is jointly trained using a multi-task loss function, consisting of a weighted sum of four terms: semantic cross-entropy loss, normal vector cosine loss, confidence binary cross-entropy loss, and texture quality mean square error loss. The network supports streaming block inference and employs an overlapping sliding window strategy, with the ratio of sliding window stride to window size ranging from 0.5 to 0.75.

[0065] The specific implementation of step S05 is as follows: based on the mean of the occlusion confidence map Mean of texture quality assessment map Mean cosine similarity to the features of the two branches According to the modal fusion quality evaluation function Calculate the fusion quality assessment value Weighting coefficient , , The sum is 1, and the coefficients range from 0.2 to 0.5. (Based on...) The values ​​dynamically adjust the geometric consistency gating threshold, batch size, and number of CUDA streams: when Take the largest batch size The number of CUDA streams is 2; when Batch size ;when Batch size The number of CUDA streams increases to 4; when Batch size The number of CUDA streams has been increased to 8. The dynamic task queue and work-stealing scheduling algorithm pre-estimates the computational cost of each task block based on the local density of the complete terrain point cloud, subdivides the task into subtasks with approximately equal computational cost, and each CUDA stream maintains an independent task queue. Idle CUDA streams steal subtasks from the tail of the queue of the most loaded CUDA stream. Combined with the asynchronous execution mechanism of the CUDA pipeline, data transmission and computation overlap, and the number of concurrent CUDA streams ranges from 2 to 8.

[0066] It should be noted that the key technologies of this invention include: a point cloud noise separation technology based on spectral signal processing transfers the spectral analysis theory in signal processing to irregular point cloud manifolds, enabling the natural separation of macroscopic terrain undulations, ground feature boundaries, and random noise in the spectral domain, overcoming the inherent defect of traditional spatial domain point-by-point filtering that cannot distinguish between terrain details and random noise; a cavity completion technology based on hydrodynamic potential flow theory utilizes the harmonicity of the Laplace equation to ensure the tangential continuity of the completed surface, and the gravity potential correction term makes the completion result conform to the physical laws of gravity erosion landforms, fundamentally avoiding the problem of geometric distortion caused by interpolation methods in large cavity areas; and a geometry-guided cross-modal dynamic graph attention fusion technology, through a differentiable K-nearest neighbor dynamic graph reconstruction module and an asymmetric bidirectional cross-attention mechanism, enables the point cloud geometric topology relationship to actively guide the extraction direction of image texture features, and the geometric consistency gating mechanism triggers an iterative refinement loop when feature semantics are inconsistent, ensuring the robustness of feature alignment in the cross-modal common embedding space. The synergistic effect of the three key technologies is reflected in the following aspects: spectral denoising provides point clouds with accurate geometric topological relationships for dynamic graph construction; potential flow completion provides geometrically continuous complete input for attention fusion; and radiometric calibration eliminates the radiometric differences between the two modes. Together, these three technologies ensure the robustness of cross-modal feature alignment in three dimensions: data quality, geometric integrity, and radiometric consistency, thus constraining the final topographic mapping results in terms of semantic consistency and physical reliability.

[0067] It should be noted that in large-scale mountainous and canyon topographic mapping, terrain occlusion leads to significant visual differences between lidar point clouds and optical images. This results in systematic deviations in the feature descriptors of the same feature in both modalities. Fixed-threshold feature matching strategies struggle to adaptively handle these differences in feature alignment quality under varying degrees of occlusion, leading to regional biases in semantic annotation and normal vector estimation. The root cause of this technical problem lies in the different imaging mechanisms of lidar and optical cameras. LiDAR actively emits laser pulses to acquire distance information, while optical cameras passively receive solar radiation reflected from terrain features to acquire texture information. In steep slopes or canyons, their effective observation areas do not overlap, resulting in significantly lower cosine similarity of feature descriptors in the common embedding space compared to flat terrain. Traditional fixed-threshold matching strategies fail to detect these feature quality changes caused by terrain occlusion, applying the same alignment constraint to all regions and causing accumulated feature alignment errors in occluded areas. Common solutions to this problem include introducing an occlusion mask to mark occluded areas and skip matching, or using single-modal point cloud features for semantic annotation of occluded areas. However, the generation of occlusion masks relies on accurate geometric registration. In large-scale complex terrain, occlusion mask errors can propagate to subsequent semantic annotations. Single-modal annotation completely abandons the potential contribution of image texture information to the semantic judgment of occluded areas, resulting in a significant decrease in classification accuracy in areas with mixed vegetation and bare rock. This invention effectively solves this technical problem. The geometric consistency gating mechanism calculates the cosine similarity of the two-branch features in real time during the inference phase. When the cosine similarity is lower than the geometric consistency gating threshold, an internal iterative refinement loop is automatically triggered. By updating the projection matrices of the Query and Key, feature alignment is gradually tightened, iterating a maximum of 3 times until convergence. This eliminates the need for pre-generating occlusion masks and achieves adaptive alignment depth adjustment for regions with different degrees of occlusion. The modal fusion quality evaluation function further integrates the mean of the occlusion confidence map, the mean of the texture quality evaluation map, and the mean of the cosine similarity of the two-branch features into a fusion quality evaluation value. Based on the fusion quality evaluation value, the batch size and the number of CUDA streams are dynamically adjusted, tilting computational resources towards areas with severe occlusion and high fusion difficulty, thus improving the fusion quality of difficult areas while ensuring overall processing efficiency. The cross-modal skip path directly injects the original geometric features of the low-level point cloud branch into the high-level semantic decision layer, so that even in the occluded areas with degraded texture, geometric constraints continue to participate in the final semantic output, avoiding the dilution of geometric information caused by occlusion of high-level fused features.

[0068] A second aspect of the present invention provides a computer-readable storage medium storing program instructions, which, when executed in a computer, are used to perform the above-described method for topographic mapping based on the fusion of laser point clouds and optical images.

[0069] A third aspect of the present invention provides a topographic mapping system based on the fusion of laser point cloud and optical image, comprising the aforementioned computer-readable storage medium. The system is any one of a computer, a server, or a microcontroller. The computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor that executes the program instructions stored in the computer-readable storage medium.

[0070] Specifically, the principle of this invention is:

[0071] The core reason why this invention can solve the above-mentioned technical problems is that the geometric topological relationship is explicitly introduced into the image feature extraction process, so that the feature extraction of the two modes is no longer independent of each other.

[0072] In traditional methods, point cloud branches and image branches extract features independently before matching is performed, and geometric information does not constrain the direction of image texture aggregation. When terrain occlusion is severe, the image branch lacks geometric guidance, and the spatial aggregation area of ​​texture features is misaligned with the boundaries of real-world features, leading to a decrease in the cosine similarity of feature vectors in the cross-modal common embedding space, ultimately causing semantic annotation bias. The logical rationality of this invention is reflected in the following aspects.

[0073] First, the differentiable K-nearest neighbor dynamic graph reconstruction module enables the gradient of graph edge weights with respect to the input geometric features to propagate backward. The graph structure is optimized end-to-end during network training, ensuring that geometric topological relationships continue to participate in feature learning throughout the entire training process, rather than being fixed and used only in the preprocessing stage.

[0074] Second, the asymmetric bidirectional cross-attention module uses point cloud feature vectors as queries and image texture gradient vectors as auxiliary keys to guide the aggregation direction of local point cloud features. At the same time, it uses point cloud normal vector fields as geometric position encoding of image branches, which geometrically constrains the spatial perception capability of image texture features and establishes cross-modal geometric coupling from two directions, rather than unidirectional information transmission.

[0075] Third, the geometric consistency gating mechanism judges the alignment quality by calculating the cosine similarity of the two branches of features. When the cosine similarity is lower than the threshold, it triggers an internal iterative refinement loop, which iterates a maximum of 3 times. By updating the projection matrix of Query and Key, the feature alignment is gradually tightened, enabling the model to adaptively adjust the alignment depth under difficult conditions such as terrain occlusion or texture degradation, rather than processing all samples with the same depth.

[0076] Fourth, the cross-modal skip path directly injects the original geometric features of the low-level point cloud branch into the high-level semantic decision layer, preventing fine-grained geometric information from being over-abstracted and diluted in the deep network, and ensuring the continuous contribution of the fine-grained terrain structure to the final semantic output.

[0077] Fifth, at the data level, in step S01, this invention decomposes the elevation signal into components with different oscillation frequencies using a spectral signal processing method. This naturally separates the macroscopic topographic undulations, feature boundaries, and random noise in the spectral domain, providing a frequency-consistent and geometrically complete denoised point cloud for subsequent cross-modal fusion, reducing noise interference with geometric topological relationships at the source. In step S02, a hole-filling algorithm based on hydrodynamic potential flow theory ensures the tangential continuity of the filled surface, providing a geometrically continuous complete point cloud for the fusion model and avoiding interference from geometric discontinuities at hole boundaries in dynamic graph construction. In step S03, a joint radiative transfer model aligns the radiative distributions of the two modes through contrastive learning loss, eliminating the interference of differences in radiative physics models on feature matching and providing radiatively consistent input data for cross-modal dynamic graph attention fusion. The above steps lay the foundation for cross-modal fusion in step S04 from three dimensions: data quality, geometric integrity, and radiative consistency. The overall scheme forms a closed loop in the logical link.

[0078] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0079] The specific implementation method of step S01 is as follows.

[0080] A K-nearest neighbor weighted graph is constructed using each point in the original point cloud as a node, where K ranges from 10 to 20. Edge weights... The value is determined by a weighted sum of the Euclidean distance and the cosine of the angle between the normal vectors, and the calculation formula is as follows:

[0081] ;

[0082] In the formula, , They are nodes , The three-dimensional coordinates, in meters. The Euclidean distance between two points is expressed in meters. The distance attenuation scale parameter, in meters, is estimated from the average distance between local points in the point cloud. The Gaussian distance decay term is dimensionless. For nodes With nodes The angle between the normal vectors, in rad. The normal vector consistency term is dimensionless. The distance weighting coefficient is dimensionless. The normal vector weight coefficients are dimensionless; after weighting the two terms... Dimensionless boundary weight The typical range is 6:4 to 7:3, determined by multiple rounds of comparative experiments on four types of terrain samples: plains, hills, mountains, and canyons, based on the reconstruction accuracy index.

[0083] Constructing the normalized graph Laplacian matrix Its eigenvalues ​​are decomposed into:

[0084] ;

[0085] In the formula, Total number of point cloud nodes The graph Fourier basis matrix has column vectors that are orthogonal eigenvectors and are dimensionless. It is a diagonal eigenvalue matrix, and all elements are dimensionless. for The transpose of the point cloud elevation signal vector. (Unit: m) Projected onto the spectral domain, the spectral coefficient vector is obtained. , The unit is meters (m), which represents the amplitude of each frequency component of the elevation signal under the Fourier basis of the graph.

[0086] Spectral domain filter response value Based on the fractal dimension of the terrain Adaptive design, fractal dimension Elevation signals at different scales Estimation of variance ratios under the following conditions:

[0087] ;

[0088] In the formula, This is a scale parameter, in meters. For scale Lower elevation signal variance, in units of and All are dimensionless logarithmic values, ratios This is a dimensionless parameter, typically ranging from 2.0 to 2.8, describing the degree of self-similarity of the terrain surface. Cutoff frequency. (Dimensionless) Derived from the mean of the power spectrum of optical image texture (Dimensionless) modulation, the modulation relationship is:

[0089] ;

[0090] In the formula, The reference cutoff frequency is dimensionless and was obtained through experimental calibration on standard terrain samples. The texture modulation coefficient is dimensionless and obtained through experimental calibration. The mean power spectrum of local texture in the radiometrically aligned image is dimensionless and is calculated from the normalized power spectrum after the two-dimensional discrete Fourier transform of the image. Dimensionless, textured areas High, uniform area Low. Spectral coefficients after filtering. (Unit: m) Processed by frequency band: Low frequency band ( The cutoff frequency (dimensionless, determined through experimental calibration) marks the boundary between low and mid frequencies, and corresponding macroscopic fluctuations are fully preserved; the mid-frequency band... The corresponding ground feature boundaries are selectively preserved according to the following attenuation function:

[0091] ;

[0092] In the formula, This is the mid-frequency attenuation coefficient, dimensionless, with a value range of 0 to 1. For the first 1 eigenvalue, dimensionless This is the attenuation steepness coefficient, dimensionless, with an empirical value of 5–20, determined experimentally; high-frequency band. For random noise, the attenuation coefficient ranges from 0.05 to 0.2, and was determined experimentally.

[0093] Finally, compressed sensing reconstruction is performed using the sparse prior of the graph signal. Compressed sensing reconstruction utilizes the sparsity of the signal in the graph Fourier domain to recover the complete signal from undersampled observations through convex optimization. The following convex optimization problem is solved:

[0094] ;

[0095] In the formula, To reconstruct the elevation signal vector, the unit is m. The sparse regularization term for the graph Fourier domain is in m units. The observation sampling matrix is ​​dimensionless. Number of observation points To observe the elevation vector, the unit is meters. For data fidelity items, the unit is... This is the regularization coefficient, in units of... Make the regularization term (unit m) and the fidelity term (unit m) equal. After weighting, the units are unified to m. The empirical range was determined through grid search experiments on the validation set, using the root mean square error of reconstruction as the metric. ~ m^{-1}$ for norm for The square of the norm; the core calculation is accelerated by randomized singular value decomposition, which is a method of approximating the low-rank singular value decomposition of a matrix through random projection, retaining only the first... 1 eigenvector The range is 50 to 200, and the overall computational complexity is... Output denoised terrain point clouds.

[0096] The specific implementation method of step S02 is as follows.

[0097] The void boundary points of the denoised terrain point cloud are treated as solid wall boundary conditions, and the void region is treated as a fluid domain. The Laplace equation is numerically solved using the finite difference method on an adaptive mesh.

[0098] ;

[0099] In the formula, The velocity potential function is expressed in units of 1000 kJ / m². For the Laplace operator, the unit is . This unifies the dimensions of all terms in the equation. The mesh resolution is determined by the curvature of the boundary points. (Unit is) Adaptive control is employed, with grid spacing ranging from 0.02 to 0.05 m in high-curvature regions and from 0.1 to 0.3 m in low-curvature regions. The curvature threshold is determined through statistical analysis based on field measurements of various terrain samples. A gravitational potential correction term is superimposed on the slope region, resulting in the following equation:

[0100] ;

[0101] In the formula, Elevation, in meters (m) The elevation gradient is dimensionless. The slope angle is expressed in rad. Dimensionless slope factor This is the topographic erosion coefficient, in units of... ,make The dimensions are Right now ,and dimension There are differences, therefore a time scale parameter is introduced. (Unit: s) Write the coefficient of the correction term as (Unit: m), to unify the dimensions of all terms in the equation. The range is 0.01 to 0.05. , The timescale of topographic erosion was estimated and determined through iterative experiments comparing measured field topographic data with numerical simulations. A radiometrically aligned image texture gradient was also incorporated. (Unit: grayscale value / m) As an anisotropic diffusion correction term, the anisotropic diffusion coefficient matrix... (Unit is) The diffusion intensity is determined by the texture gradient magnitude. At the texture boundary, the diffusion intensity is reduced to preserve the geometric details of the ground feature boundary. The complete equation is expressed as:

[0102] ;

[0103] In the formula, Units are ,through The unit after action is Dimensionally unified with the second term Calculated by the following formula:

[0104] ;

[0105] In the formula, The isotropic diffusion coefficient is expressed in units of 1000 ppm. The empirical value range is 0.01 to 0.1. Determined through experimental fitting on standard terrain samples These are the anisotropic component coefficients, in units of... The empirical value range is 0.001 to 0.01. Determined through experimental fitting This represents the texture gradient magnitude, expressed in grayscale values ​​per mille (m). The texture gradient attenuation scale, measured in grayscale values ​​per mille (m), was determined through experimental calibration. Dimensionless Dimensionless attenuation factor for Identity matrix, dimensionless This is the texture gradient vector, with units of grayscale values ​​per m. This is the outer product matrix, with units of 1. To prevent division by zero, the unit is... Experience points are taken (grayscale value / m)^{2}$ Units are The product matrix becomes dimensionless after dividing by the outer product matrix of the numerator, and multiplied by... The following unit is ,and The dimensions are unified; the positions of the completion points are generated along the gradient direction of the potential field to ensure the tangential continuity between the completed surface and the boundary points; the iterative convergence is accelerated by the multigrid preconditioning conjugate gradient method, a numerical method that alternately solves the linear equation system on multiple resolution grids to accelerate iterative convergence, with a computational complexity of O(n log n). The point density is adaptively matched with the density of the surrounding known point cloud to output a complete terrain point cloud.

[0106] The specific implementation method of step S03 is as follows.

[0107] Flat field correction is achieved by extracting the gain coefficient per pixel. (Dimensionless) The acquired image is normalized pixel by pixel. Pixels in a flat-field image acquired under uniform diffuse light source The ratio of the response value at a given location to the mean value of the entire graph is calculated, where... , These represent the row and column coordinates of the pixels, in pixels. Atmospheric radiative correction is performed by removing the effects of scattering and absorption based on the atmospheric radiative transfer model. Atmospheric optical path parameters are obtained through synchronously measured atmospheric parameters or estimated based on the dark pixel method. Normalization of the incident angle affects the laser return intensity. (Unit: W / sr) According to the angle of incidence (Unit: rad) Cosine effect correction, distance compensation restored according to the squared distance attenuation law, calibration coefficients obtained by fitting multi-distance measured data collected on a standard target with known reflectivity. The joint radiative transfer model aligns the two-modal radiative distribution in the common semantic feature space through contrastive learning loss, the contrastive learning loss function is:

[0108] ;

[0109] In the formula, The number of samples in the batch is dimensionless. For the first A sample optical image radiation-normalized feature vector, dimensionless For the first Normalized feature vectors of reflectance intensity of each sample point cloud, dimensionless. For feature vector dimension The dot product of two vectors, dimensionless. This is a temperature coefficient, dimensionless, and typically takes a value of 0.07. The value represents the dimensionless contrast loss; the output results are denoted as the radiometric aligned image and the radiometric calibration point cloud, respectively.

[0110] The specific implementation method of step S04 is as follows.

[0111] The first branch of the improved multi-scale point network backbone introduces a differentiable K-nearest neighbor dynamic graph reconstruction module at each layer. This module refers to a module where the gradient of the graph edge weights with respect to the input geometric features exists and can propagate backwards, allowing the graph structure to be optimized end-to-end during network training. It is determined by a weighted average of three geometric features: the distance between points, the angle between normals, and the difference in curvature.

[0112] ;

[0113] In the formula, for Normalize the activation function so that Dimensionless The local neighborhood radius, in meters, is estimated from the local density of the point cloud. Dimensionless , They are nodes , Curvature, in units of The local maximum curvature is expressed in units of . ,make Dimensionless , , The weights are learnable, dimensionless, and obtained through end-to-end training; the K value ranges from 8 to 32, and is determined iteratively using feature consistency loss as an indicator through ablation experiments on terrain samples with different point cloud densities. The second branch's hierarchical visual transformation backbone extracts multi-scale hierarchical features, with deformable convolutional modules introduced at each scale, and offsets... (Unit: pixels) Obtained through convolutional learning, where the subscripts are... It serves as an index for sampling points of the convolution kernel, adapting to the distortion of the feature space caused by ground feature deformation.

[0114] The two branches interact across three scale levels, with each level employing an asymmetric bidirectional cross-attention module. This module means the two branches have asymmetrical roles in attention computation: the first branch provides the query matrix based on point cloud features, while the second branch provides the key and value matrices based on image features. These roles are not interchangeable. The attention weight matrix... The (dimensionless) calculation is as follows:

[0115] ;

[0116] In the formula, For point cloud feature query matrix, The feature matrix of the point cloud is dimensionless. The projective matrix is ​​a learnable query matrix, dimensionless. The image feature key matrix, The image feature matrix is ​​dimensionless. The learnable bond projection matrix is ​​dimensionless. The texture gradient auxiliary key matrix is ​​obtained by linear mapping of the image texture gradient and is dimensionless. The dimension of the key vector is dimensionless. Number of feature points in a point cloud Number of image feature points For a dimensionless inner product matrix, divide by (Dimensionless) after Output dimensionless weight matrix Value matrix Provided by image features, The projection matrix is ​​a learnable value and is dimensionless. The second branch image features are derived from the normal vector field of the first branch point cloud. (Dimensionless after normalization) as a geometric position encoding, encoding vector Injecting image features to enhance spatial perception The learnable mapping matrix is ​​dimensionless. Dimensionless.

[0117] Geometric consistency gating mechanisms are based on the cosine similarity of features in the two branches. (Dimensionless) Determine if the internal iterative refining loop is triggered:

[0118] ;

[0119] In the formula, The inner product of eigenvectors, dimensionless. and These are the two eigenvectors respectively. Norm, dimensionless Dimensionless; when Below the gate threshold When the value (range 0.5–0.7, determined through grid search experiments using the average intersection-union ratio of semantic labels on the validation set) is reached, the embedded lightweight iterative registration sub-network is triggered, updating the projection matrix in each iteration. , ,in For the number of iterations, and The incremental update matrix output by the subnetwork is dimensionless and iterates at most 3 times until convergence.

[0120] The network design employs cross-modal skipping paths, directly injecting the original geometric features from the lower layers of the first branch into the higher-level semantic decision layer. This prevents fine-grained geometric information from being overly abstracted and diluted in deep networks. Within the same branch, cross-layer skipping connections are simultaneously preserved, and the feature fusion weights for both types of skipping paths are adaptively allocated by hierarchical attention gating. Normal vector cosine loss is used. The formula used to constrain the accuracy of surface normal vector estimation is as follows:

[0121] ;

[0122] In the formula, Number of valid sampling points The network estimated normal vector is dimensionless. The true value normal vector is dimensionless. The inner product is dimensionless. and For the corresponding vector Norm, dimensionless This is a dimensionless loss value, ranging from 0 to 2. The multi-task loss function consists of a weighted sum of four terms:

[0123] ;

[0124] In the formula, Semantic cross-entropy loss, dimensionless The confidence level is a binary cross-entropy loss, dimensionless. The mean square error loss for texture quality is dimensionless. , , , These are the weighting coefficients, dimensionless, determined iteratively on the validation set through multiple rounds of ablation experiments. The sum of the four coefficients is typically 1. The total loss value is dimensionless. The network supports streaming block inference and adopts an overlapping sliding window strategy. The ratio of the sliding window step size to the window size ranges from 0.5 to 0.75, and is determined after experimental verification of the minimum overlap rate required for seamless stitching. During training, a curriculum learning strategy is adopted. The curriculum learning strategy refers to the training strategy of gradually increasing the difficulty of training samples from easy to difficult to alleviate the gradient instability problem in the early stage of training. It gradually transitions from flat terrain samples to mountainous samples, and the gradient norm stability is used as the criterion to determine the curriculum advancement node.

[0125] The specific implementation method of step S05 is as follows.

[0126] Modal fusion quality evaluation function based on occlusion confidence map mean (Dimensionless) Mean of Texture Quality Assessment Map (Dimensionless) Mean cosine similarity of features between two branches (Dimensionless) Calculation of fusion quality assessment value (Dimensionless):

[0127] ;

[0128] In the formula, , , The three weights are dimensionless, with a sum of 1. Each coefficient ranges from 0.2 to 0.5. They are iteratively determined using the average intersection-union ratio of semantic tags as the evaluation metric through validation set comparison experiments. Maximum batch size. Based on video memory capacity (Unit: bytes) and single-sample memory usage The quotient (in bytes) is determined as follows:

[0129] ;

[0130] In the formula, For floor operations, It is a dimensionless integer. Based on... Dynamically adjust the gate threshold value Batch size and number of concurrent streams. The computational load of each task block is estimated in advance based on the local density of the complete terrain point cloud. The task is subdivided into subtasks with approximately equal computational load. Each concurrent stream maintains an independent task queue. Idle streams steal subtasks from the tail of the queue of the most loaded stream. With the help of the pipeline asynchronous execution mechanism, data transmission and computation overlap. The number of concurrent streams ranges from 2 to 8. After experimentally measuring the processor utilization under each concurrency level, the optimal value is selected, and the terrain mapping results are generated in parallel.

[0131] To better understand and implement this invention, the following is a specific application scenario of this invention, Example 2:

[0132] To verify the effectiveness of the invention, technicians constructed a test environment and used an airborne lidar and a synchronous optical multispectral camera to collect data on a complex terrain area containing mountains, canyons, and vegetation cover. The test area was approximately 8... The terrain has an elevation difference of approximately 620 m, and the vegetation coverage is approximately 63%. The average density of the lidar point cloud is 42 points / ... The optical multispectral image has a ground resolution of 0.12 m and includes four bands: blue, green, red, and near-infrared. Multiple point cloud holes exist within the survey area due to ridge obstruction, accounting for approximately 7% of the total area. This represents a typical scenario with high cross-modal fusion difficulty.

[0133] After data acquisition, technicians first performed step S01 on the original point cloud, conducting point cloud noise separation and terrain frequency analysis based on spectral signal processing. A K-nearest neighbor weighted graph was constructed with each point in the original point cloud as a node, with K set to 15 and the weight ratio of Euclidean distance to the cosine of the angle between the Euclidean distance and the normal vector set to 7:3. The point cloud elevation signal was projected into the spectral domain through eigenvalue decomposition of the normalized graph Laplacian matrix. The terrain fractal dimension estimation results are shown in Table 1. A spectral domain filter was adaptively designed based on the fractal dimension, with the number of retained eigenvectors k set to 128. High-frequency noise components were effectively attenuated, while mid-frequency components of feature boundaries were selectively preserved, resulting in a denoised terrain point cloud. The spatial distribution of the spectral domain filter cutoff frequency with the texture power spectrum is shown in Table 1. Figure 3 As shown, the cutoff frequency is significantly higher in areas with rich texture and lower in areas with uniform texture, indicating that the adaptive adjustment of the filter is as expected.

[0134] Table 1. Estimation results of topographic fractal dimension for various terrain zones in the survey area.

[0135]

[0136] Then, step S02 is executed, employing an adaptive point cloud cavity completion algorithm based on hydrodynamic potential flow theory to complete the cavity regions in the denoised terrain point cloud. The cavity boundary points are treated as solid wall boundary conditions, and the Laplace equation is solved on an adaptive mesh. The mesh spacing is 0.03 m for high curvature regions and 0.2 m for low curvature regions. A gravitational potential correction term is superimposed on the slope region, and the terrain erosion coefficient is set to 0.03. .like Figure 4 As shown, the comparison of tangential continuity of the point cloud cavity region before and after completion indicates that the potential flow completion method achieves a smooth transition at the cavity boundary, without the geometric distortion phenomenon common in interpolation methods, and outputs a complete terrain point cloud.

[0137] Step S03 involves performing flat-field correction and atmospheric radiative correction on the optical multispectral image, and normalizing the incident angle and calibrating the distance compensation for the reflection intensity of the complete terrain point cloud. Flat-field correction uses the gain coefficient per pixel extracted from flat-field images acquired under uniform diffuse light, and atmospheric optical path parameters are estimated using the dark pixel method. Distance compensation calibration coefficients are obtained by fitting multi-distance measured data collected on a standard target with known reflectivity; the fitting results are shown in Table 2. The joint radiative transfer model aligns the two-mode radiative distributions using a contrastive learning loss; the feature cosine similarity distributions before and after alignment are shown in Table 2. Figure 5 As shown, the mean cosine similarity is significantly improved after alignment, and the output is a radiometric aligned image and a radiometric calibration point cloud.

[0138] Table 2. Fitting Results of Distance Compensation Calibration Coefficients

[0139]

[0140] Execute step S04, inputting the radiometric calibration point cloud and radiometrically aligned image into the geometry-guided cross-modal dynamic graph attention fusion model. For example... Figure 2 As shown, the model includes a point cloud branch and an image branch. The point cloud branch extracts geometric features using a differentiable K-nearest neighbor dynamic graph reconstruction module, with K set to 16. The two branches interact across three scale levels through an asymmetric bidirectional cross-attention module. The geometric consistency gating threshold is set to 0.62. The proportion of internal iterative refinement loops triggered in the canyon sidewall occlusion area is approximately 34%, with an average of approximately 2.1 iterations, indicating that the gating mechanism effectively identifies the feature alignment difficulties in occlusion areas. The average intersection-union (IUCN) ratios of semantic annotation results for each terrain region are shown in Table 3. The average IUCN ratio for bare rock and vegetation classification in the canyon sidewall area is significantly improved compared to the single-modal method, indicating that cross-modal geometric guidance has a substantial improvement effect on the semantic annotation of occlusion areas. Occlusion confidence maps for different terrain regions are shown in Table 3. Figure 6As shown, the confidence level of shading in the canyon sidewall area is significantly lower than that in the main ridge area of ​​the mountain, which is consistent with the actual shading distribution.

[0141] Table 3. Average Intersection over Union (IoU) Ratio Results of Semantic Labels for Various Topographic Zones

[0142]

[0143] In step S05, based on the occlusion confidence map and texture quality evaluation map, the fusion quality evaluation value is calculated using the modal fusion quality evaluation function. Weighting coefficient , , We take values ​​of 0.4, 0.35, and 0.25 respectively. In the area shaded by the canyon sidewall, The value is approximately 0.41, triggering concurrent processing of 4 CUDA streams, with a batch size of [value missing]. In the flat valley area, The value is approximately 0.83, using 2-way CUDA streams, and the batch size is... The dynamic task queue and work-stealing scheduling algorithm estimates the computational load of task blocks based on local point cloud density. Idle CUDA streams steal subtasks from the tail of the heaviest queue in real time, thus balancing the load across CUDA streams. The final output is terrain mapping results, including a dense point cloud semantic label map, a surface normal field map, and a texture quality assessment map, such as... Figure 7 , Figure 8 As shown in Table 4, the GPU utilization statistics for each processing stage show that the work-stealing scheduling algorithm keeps GPU utilization at a high level in complex occlusion areas.

[0144] Table 4. GPU utilization statistics for each processing stage

[0145]

[0146] Spatial distribution of surface normal vector estimation as follows Figure 9 As shown, the continuity of the normal vector field between the main ridge region and the canyon sidewall region is good, and no jump in the normal vector is caused by point cloud voids. This indicates that the synergistic effect of potential flow completion and cross-modal fusion effectively ensures the geometric continuity of the normal vector estimation. The spectral response curves before and after radiation alignment for each band are shown in the figure. Figure 10 As shown, after alignment, the distribution of radiation values ​​in each band tends to be consistent, and the difference in radiation between the two modes is effectively eliminated.

[0147] This invention brings the following technological advancements compared to traditional methods. In point cloud denoising, traditional spatial domain filtering attenuates all high-frequency components on irregular point cloud manifolds, failing to distinguish between real terrain details and random noise. This invention decomposes the elevation signal into components with different oscillation frequencies using the graph Laplacian operator, with the cutoff frequency modulated in real-time by the optical image texture power spectrum. This enables the denoising process to possess ground feature perception capabilities, preventing excessive smoothing of terrain details. In hole completion, traditional interpolation methods produce geometric distortions in large hole areas due to a lack of physical constraints. This invention introduces the harmonicity of the Laplacian equation to ensure tangential continuity of the completed surface, a gravitational potential correction term to make the completion result conform to the physical laws of gravity erosion landforms, and an anisotropic diffusion correction term to introduce radiation texture constraints, ensuring consistency of the completed area from both geometric and aesthetic dimensions. In cross-modal fusion, traditional methods rely on manually designed descriptors for fixed threshold matching, and geometric information does not constrain image feature extraction. This invention establishes bidirectional coupling between geometry and texture through a differentiable K-nearest neighbor dynamic graph reconstruction module and an asymmetric bidirectional cross-attention mechanism. The geometric consistency gating mechanism adaptively adjusts the alignment depth for regions with different degrees of occlusion, thus eliminating the alignment error caused by the independent extraction of features from the two modalities.

[0148] It should be noted that the variables involved in this invention are explained in detail in Tables 5 and 6.

[0149] Table 5. Variable Explanation Table (Part 1)

[0150]

[0151] Table 6. Variable Explanation Table (Part Two)

[0152]

[0153] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for topographic mapping based on fusion of laser point cloud and optical image, characterized in that, Includes the following steps: The original point cloud acquired by lidar is subjected to point cloud noise separation and terrain frequency analysis based on spectral signal processing, and the denoised terrain point cloud is output. The hollow regions in the denoised terrain point cloud are filled by an adaptive point cloud hole completion algorithm based on hydrodynamic potential flow theory, and the complete terrain point cloud is output. Flat field correction and atmospheric radiative correction are performed on optical images. The incident angle is normalized and distance compensation is calibrated for the reflection intensity of complete terrain point clouds. A joint radiative transfer model is established, and radiative aligned images and radiatively calibrated point clouds are output. The radiometrically calibrated point cloud and the radiometrically aligned image are input into the geometry-guided cross-modal dynamic graph attention fusion model. The geometry-guided cross-modal dynamic graph attention fusion model includes a point cloud branch and an image branch. The point cloud branch extracts geometric features using a differentiable K-nearest neighbor dynamic graph reconstruction module. The two branches interact across modalities through an asymmetric bidirectional cross attention module. The point cloud feature vector is used as a query to guide image feature aggregation, and the point cloud normal vector field is used as the geometric position encoding of the image branch. The output includes dense point cloud semantic labels, surface normal vector estimation, occlusion confidence map, and texture quality assessment map. Based on the occlusion confidence map and texture quality assessment map, a dynamic task queue and work-stealing scheduling algorithm are used to distribute tasks to multiple CUDA streams and generate terrain mapping results in parallel.

2. The method according to claim 1, wherein, The point cloud noise separation and terrain frequency analysis based on spectral signal processing specifically involves constructing a K-nearest neighbor weighted graph with each point in the original point cloud as a node. The edge weights are jointly determined by the Euclidean distance and the cosine of the angle between the normal vectors. The eigenvalue decomposition of the normalized graph Laplacian matrix is ​​calculated, the point cloud elevation signal is projected onto the spectral domain, and a spectral domain filter is adaptively designed according to the terrain fractal dimension. Compressed sensing reconstruction is performed through the sparse prior of the graph signal.

3. The topographic mapping method based on the fusion of laser point clouds and optical images according to claim 2, characterized in that, The cutoff frequency of the spectral domain filter is modulated by the power spectrum of the optical image texture. Low-frequency coefficients corresponding to the macroscopic undulations of the terrain are preserved, mid-frequency coefficients corresponding to the boundaries of ground features are selectively preserved, and high-frequency coefficients corresponding to random noise are attenuated. The computational core of the compressed sensing reconstruction is approximately accelerated by randomized singular value decomposition.

4. The topographic mapping method based on the fusion of laser point clouds and optical images according to claim 3, characterized in that, The adaptive point cloud cavity completion algorithm based on hydrodynamic potential flow theory specifically treats the cavity boundary points of the denoised terrain point cloud as solid wall boundary conditions and the cavity region as a fluid domain. The Laplace equation is numerically solved using the finite difference method on an adaptive grid. The grid resolution is adaptively controlled by the curvature of the boundary points. Gravitational potential correction terms are superimposed along the slope direction in the slope region. The texture gradient of the radiation-aligned image is injected into the Laplace equation as an anisotropic diffusion correction term.

5. The topographic mapping method based on the fusion of laser point clouds and optical images according to claim 4, characterized in that, The numerical solution of the Laplace equation uses the multigrid preconditional conjugate gradient method to accelerate iterative convergence. The completion points are generated along the gradient direction of the potential field, and the density of the completion points is adaptively matched with the density of the surrounding denoised terrain point cloud.

6. The topographic mapping method based on the fusion of laser point clouds and optical images according to claim 5, characterized in that, The joint radiative transfer model aligns the two modal radiative distributions in a common semantic feature space by using contrastive learning loss to match the radiative values ​​of the optical image after flat field correction and atmospheric radiative correction with the reflection intensity of the complete terrain point cloud after incident angle normalization and distance compensation calibration.

7. The topographic mapping method based on the fusion of laser point clouds and optical images according to claim 6, characterized in that, In the geometrically guided cross-modal dynamic graph attention fusion model, the point cloud branch adopts an improved PointNet++ backbone network, and the image branch adopts a Swing Transformer backbone to extract multi-scale hierarchical features and introduces deformable convolutional modules at each scale. The features of the two branches interact across three scale levels, and an asymmetric bidirectional cross attention module is designed at each level.

8. The topographic mapping method based on the fusion of laser point clouds and optical images according to claim 7, characterized in that, The geometry-guided cross-modal dynamic graph attention fusion model introduces a geometric consistency gating mechanism to calculate the cosine similarity of features between two branches. When the cosine similarity is lower than the geometric consistency gating threshold, an internal iterative refinement loop is triggered. The internal iterative refinement loop is an embedded lightweight iterative registration sub-network that updates the projection matrices of Query and Key in each iteration.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores program instructions, which, when executed in a computer, are used to perform the topographic mapping method based on the fusion of laser point clouds and optical images as described in any one of claims 1-8.

10. A topographic mapping system based on the fusion of laser point clouds and optical images, characterized in that, The system comprises the computer-readable storage medium of claim 9, wherein the system is a computer, the computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor that executes program instructions stored in the computer-readable storage medium.

Images