CN105336003B - The method for drawing out three-dimensional terrain model with reference to the real-time smoothness of GPU technologies - Google Patents

Abstract

A kind of method that real-time smoothness of combination GPU technologies draws out three-dimensional terrain model, belongs to image processing technique field.The purpose of the present invention is using global digital elevation model as data source, Programmable GPU technology based on current hot topic, caching multiplexing during multiple drafting is realized, efficiently reduces the method that the real-time smoothness of combination GPU technologies of space load draws out three-dimensional terrain model that calculates.The present invention step be：Multi-resolution pyramid model construction, the elimination of image noise, image is filtered, to the earth after plane projection according to etc. longitudes and latitudes carry out piecemeal, while according to top-down mode, build the pyramidal layer of different levels rank.The present invention is based on the acceleration and enhancing that Programmable GPU technology realizes terrain rendering, i.e. by using each stage of Shader Language control pattern rendering pipeline, the vertex information of altitude data and index information are generated into two texture deposit video memorys respectively, for entire landform draw call；Using tessellation and fractal technology in geometry stage opposite vertexes interpolation and offset, generating process details, when making up lack of resolution the phenomenon that terrain mesh corner angle.

Description

Method for smoothly drawing three-dimensional terrain model in real time by combining GPU technology

Technical Field

The invention belongs to the technical field of image processing.

Background

With the increasing maturity of aerospace science and technology and remote sensing technology, the spatial resolution of sensors carried by reconnaissance satellites and various aircrafts surrounding a near-earth orbit is continuously improved, and the digital photogrammetry technology is more and more precise and accurate. A data set formed by gathering various remote sensing data is an important information source of geospatial information, and the remote sensing images and digital elevations are the basis and the root of the whole spatial data frame. On the premise that the application demand of the basic geographic data is rapidly increased, the data volume of satellite remote sensing images and digital elevations is increasingly huge, how to effectively manage and utilize the data volume, and the problem becomes an urgent need to be solved in the field of remote sensing data processing and geographic information systems.

The development of modern military technology is continuously pushing the battlefield space to the direction of multi-dimension and multi-element. In order to meet complex, staggered and instantaneously changeable battlefield environments and situations, how to assist a commander to comprehensively deploy operation actions and accurately master battlefield information becomes an important target and content of military information. Meanwhile, the battle space of people is continuously expanded, the development of geospatial information is slightly lagged, the data acquisition and visualization capacity cannot meet the real-time and accurate requirements due to the development limitation of a sensor technology and a computer graphic technology in the early stage, the process of understanding and expressing the spatial information can be simplified by a two-dimensional plane processing mode, but the authenticity and the integrity of the spatial information are sacrificed, particularly the topological relation between the elevation information and the three-dimensional space is sacrificed. However, as the remote sensing science and technology and computer graphics are changing day by day, the efficient and vivid three-dimensional display of the obtained high-resolution elevation and image is no longer a technical challenge difficult to realize, and the development of the subject field also prompts military information to make up for the technical defect of insufficient three-dimensional expression capability.

In the military field, the development of war forms requires that the commander can not master the battlefield environment any more on a certain fixed area or information in a certain specific form, and the commander needs to control the battlefield timely, flexibly and omnidirectionally to acquire the victory of local war under modern informatization conditions, and whether a plurality of elements of the battlefield environment can be accurately, clearly and comprehensively displayed in real time, so as to carry out auxiliary decision-making for the commander, which is a serious challenge for military information.

Disclosure of Invention

The invention aims to realize cache multiplexing during multiple drawing by taking a global digital elevation model as a data source and based on the current popular programmable GPU technology, and effectively reduces a method for calculating space load and smoothly drawing a three-dimensional terrain model in real time by combining the GPU technology.

The method comprises the following steps:

constructing a multi-resolution pyramid model of the three-dimensional terrain model:

(1) preprocessing flow of topographic data: filtering the original image to obtain a noise-free image, and constructing a multi-resolution pyramid;

i, eliminating image noise:

the noise model conforms to the characteristics of additive noise:

wherein, I (x, y) is a gray value at the image (x, y), O (x, y) is a gray value without noise pollution, n (x, y) is a random noise value independent of O (x, y), and p is a probability of noise pollution;

the process of filtering the image comprises the following steps:

calculating the uniformity image _ HSV of the whole image, and traversing pixels of the whole image;

uniformity is a parameter that measures the contrast of noise signals in an image, image I _M×N And a window W _r×r The uniformity of (c) is defined as follows:

assuming that the size of the image I is M × N, the gray-scale value of any point (x, y) in the image is f (x, y), and the average gray-scale value of the whole image is average (I) _M×N ) The calculation formula of (2) is as follows:

the size of the sliding window is r × r, and if the coordinates of the center pixel of the window are (i, j) and the gray value is f (i, j), the average value of the pixels in the window is average (W) _r×r ) The calculation formula of (2) is as follows:

(ii) counting the maximum value point pixel value iMax in the current window, the number nMax of the maximum value point pixel values iMax, the number sumPixel of the pixels in the window and the block uniformity wnd _ HSV of the window;

(iii) judging whether the gray value I (x, y) of the current window center pixel point (x, y) is equal to the pixel maximum value iMax in the window or whether the window uniformity wnd _ HSV is greater than the image uniformity image _ HSV;

(iv) performing Lagrange interpolation on the four directions which are judged as the noise points; the interpolation method of Lagrange is as follows:

wherein, y _i Represents x _i Outputting a function value; the interpolation algorithm is used for processing one-dimensional data, when the two-dimensional image is expanded, adjacent pixel points in four directions including horizontal, vertical and two diagonal directions of an interpolation pixel point (x, y) are taken for interpolation, and the average value of the interpolation pixel points is taken as the final interpolation result;

(2) constructing a multi-resolution pyramid by the data after the noise points are eliminated:

the pyramid structure is a multi-resolution hierarchical model and is constructed by a multiplying power method, namely, for the same region, a series of grids are adopted for representation, and the sampling intervals of upper and lower adjacent grids have the same proportional coefficient m, so that the resolution is reduced layer by layer from bottom to top of the pyramid, and the grid size is correspondingly reduced;

let the size of the original grid be S ₀ Resolution is R ₀ If the constructed pyramid has l layers and the pyramid top is the 0 th layer, then the grid size s and resolution r of the k-th layer can be expressed as:

partitioning the earth after plane projection according to equal longitude and latitude, and simultaneously constructing pyramid layers of different levels according to a top-down mode, wherein the construction scheme is as follows:

the latitude span of the earth is-90 degrees to 90 degrees, the longitude span is-180 degrees to 180 degrees, and the width-height ratio on a projection plane is just 2;

(II) the ratio of the resolutions of LOD of each layer is 2, namely the resolution of the sub-block is twice that of the parent block, the pixel size of the sub-block is 2 times that of the previous level in the same geographic space, and the geographic space area covered by the parent block is 4 times that of the sub-block in the same screen resolution;

(III) the block sequence of each layer adopts simple row-column coding, the ratio of the number of blocks between the next layer and the previous layer and the ratio of the number of blocks between the horizontal layer and the longitudinal layer of each layer are both 2, according to the habit of a computer image coordinate system, the starting point is defined at the upper left corner, the line number is from top to bottom, and the column number is from left to right; the number of blocks of layer 0 is 4X 2, and the number of blocks of layer k-1 is 2 ^k+1 ×2 ^k ；

(IV) each partition contains DEM data for a corresponding geographic location, which must be odd in resolution size, i.e., (2n + I) x (2n + I);

(V) the bottom layer of the pyramid is the highest resolution, and is directly obtained from original data, and the LOD level where the pyramid is located is determined by the resolution and the range of the pyramid; the resolution of texture data and normal data corresponding to each terrain block must be higher than that of DEM (digital elevation model) so as to ensure smooth interpolation and avoid rendering distortion; assuming that the DEM resolution of a certain layer is r, the resolution of the DEM is defined as r' = nxx (r-l), N is larger than l, and N belongs to N; a global scale pyramid model is obtained.

And (3) based on a GPU (graphics processing unit) pipeline drawing enhancement technology, drawing the data of the global pyramid model:

storing the vertex information and the index information of the elevation separately by using shader language, and transmitting the vertex information and the index information as texture data into a video memory at one time;

a. data structure and organization of nodes:

a chunk quad-tree structure is constructed by adopting a Chunked LOD algorithm, namely, a single vertex or a triangle is not used as a processing unit, but a grid model including a texture mapping is used for packaging the triangle into strips to complete rendering at one time, and the batch processing capacity of a GPU is fully utilized;

b. multiplexing of the buffer: the vertex coordinates (x, y) in the same layer are all transformed in the same way by adopting the method used by Geo-Clipmap;

c. and (3) data of the pyramid model are drawn in real time: the method is divided into two types, one is a static node, and the other is a dynamic node;

i, static node data rendering: during rendering, the nodes with the corresponding resolution are transferred into a video memory to be directly rendered;

(1) bit operation of ring address:

using a grid cache with a fixed size to acquire elevation data of any node in the quadtree, and rendering the whole terrain space by setting an offset address and a geographical space span; before the visual pyramid clipping is added, a visual point is taken to be vertically downward, the visual range of the visual point is a pixel grid with n multiplied by n fixed size, and the geographic spatial resolution r represented by each pixel in the layer is obtained according to the pyramid level l where the visual point is located _l ＝2 ^-l-9 X π × R, where R represents the radius of the earth, in (gx) _diff ，gy _diff ) Representing the displacement of the viewpoint projected on the terrain plane, the distance of the pixels where the viewpoint moves on the plane is represented asAndthe pixel range covered by the field of view is (x) _min ，y _min ) And (x) _max ，y _max ) Representing;

the node which moves out of the view field is replaced by the node which enters the view field for the first time, and the rest nodes are fixed, so that increment updating is completed, and for any point (x, y) in the view field range, the position of the node in the node array is (x mod n, y mod n), wherein 'mod' represents modular operation;

(2) determination of an error metric:

the error measurement calculation mode of the round LOD is adopted, namely, the errors of the nodes are calculated layer by layer from the leaf nodes along the quadtree, and the error measurement epsilon of one node _c Epsilon of its four child nodes _a Is added to its own error metric epsilon _■ ；ε _c The calculation method of (a) is as follows:

therefore, when the screen space error is larger than the threshold value, the node with smaller error is selected to be continuously searched downwards until the requirement is met or the leaf node is reached;

(3) visibility elimination:

the elimination method when the viewpoint is closer:

when the viewpoint is closer to the terrain, an intersection test is carried out by adopting the bounding box and the view cone, namely, the feature points of the bounding box and 6 planes of the view cone are sequentially intersected to determine whether the feature points belong to the half space of the plane forward direction or are intersected with the plane; the 6 planes of the view cone are determined by the equation f (x, y, z) = ax + by + cz + dw =0, where the vector (a, b, c) represents the normal vector of the plane and the value dw represents the distance of the plane from the origin; the 9 characteristic points of the bounding box consist of 8 angular points and a central point, the visibility of the terrain grid can be judged by comparing the characteristic points of the bounding box, and the size is equal to half of the length of the diagonal line of the bounding box;

the elimination method when the viewpoint is far:

when the viewpoint is far away from the terrain, judging by using a distance coefficient between the node and the viewpoint, and defining the size of the distance coefficient as the actual distance represented by the node bounding box, wherein the ratio of the length L of a projection plane to the distance D between the bounding box and the projection plane is a = L/D;

taking the coordinates v (x) of the viewpoint ₀ ，y ₀ ，z ₀ ) The bounding box center coordinate is c (x) ₁ ，y ₁ ，z ₁ ) If the direction vector from the viewpoint to the ground is d (m, n, p) = v-c;

let 8 vertex coordinates of bounding box be p _i (x _i ，y _i ，z _i ) Passing point P _i The equation of a straight line parallel to d is:

f(x，y，z)＝m(x-x _i )+n(y-y _i )+p(z-z _i )＝0 (3.9)

the near shearing surface of the visual cone body is taken as a projection plane which is vertical to a direction vector d, and the plane equation is set as

mx+ny+pz＝0 (3.10)

By combining the above two equations, p can be obtained _i Coordinates p of the projected point on the plane _i ' is a

Calculate p for all 8 vertices _i Then, taking maximum and minimum values min X, max X, min Y, max Y, min Z and maxZ in 8 coordinate points; order to

The value of the distance coefficient a is:

when the viewpoint height is larger than a set threshold value, if the a value of the node bounding box is smaller than 1, namely the size of the node bounding box is less than one pixel on a projection plane, the node is considered invisible and is not rendered any more, otherwise, the node is rendered continuously;

(4) gap elimination and geometric transition:

before node screening, firstly, calculating the geographical range covered by each LOD level so as to generate the quantity of triangles with approximately consistent granularity on a window plane during rendering;

let the starting distance of LOD of each layer be start _i I represents the level of the current LOD, the range of the covered distance is range, the distance from the viewpoint to any vertex of the current LOD is dist, and then the current transition region [ mStart, mEnd ] is]And the transition coefficient k of the vertex is:

wherein clamp (x, a, b) is a clamping function, i.e.

After obtaining the transition coefficient, defining the size of the mesh where the vertex is located as n × n, the depth of the quadtree where the vertex is located as l, and the normalized coordinate of the vertex in the mesh as (x, y), and then the actual distance scaling value scale of the mesh size relative to the world coordinate and the coordinate (x ', y') after the transition are obtained by the following formula:

where frac (x) represents the fractional part of x;

finally, the vertex coordinates are transformed from (x, y) to the position of (x/2, y/2), and the corresponding elevation is obtained through bilinear interpolation;

II, dynamic node: when the quad-tree structure extends to the leaf nodes, the requirement of a viewpoint threshold value cannot be met, and the nodes are generated by GPU real-time calculation and added into the quad-tree;

and (3) generation of dynamic nodes:

(1) subdivision of terrain models: the subdivision is that face points and edge points are inserted into the original vertexes of the quadrilateral meshes according to subdivision masks and are connected with the original vertexes to form a new quadrilateral, the degrees of the generated non-boundary vertexes are 4, the degrees of the boundary vertexes are 3 and are all regular points;

the regular point subdivision process: to quadrangle V ₅ V ₆ V ₇ V ₈ When the subdivision is carried out, four edge points E are formed on four edges of the division ₀ E ₁ E ₂ E ₃ Meanwhile, a dough point F is also provided, so that the original dough sheet is split into four sub dough sheets; e ₀ ，E ₁ ，E ₂ ，E ₃ The value of (a) is as follows:

the value of the face point F is given by equation 4.2:

wherein the content of the first and second substances,

(2) random fractal: representing the coordinates of a point on the surface by a variable (X, y), wherein X (X, y) represents the height of the point on the surface; the FBM surface is defined as a random process with an index of H (0 < H < 1) in a certain probability space

A random variable X: r → R ² Satisfies the following conditions:

(1) With probability 1, X (0, 0) =0, i.e. the process starts from the origin; and X (X, y) is a continuous function with respect to (X, y);

(2) For any variable (x, y), (h, k) e R ² With a two-dimensional increment X (X + h, y + k) obeying the expectation of 0 and a variance of (h) ² -k ² ) ^H Is normally distributed with a probability satisfying

Two fractal algorithms were used to generate different patterns of noise maps, the Diamond-Square algorithm and Perlin noise:

Diamond-Square algorithm: in the iteration process, each subdivision can execute two steps of operation, namely a Diamond step and a Square step; the Diamond step is used for generating the center of a square, taking four vertexes to calculate an average value, and adding a random quantity to serve as a value of a central point; the Diamond step generates four triangular pyramids; the Square step is to take the four vertexes of the pyramid, and take the average value of the four vertexes and the random variable as the midpoints of the four sides of the original Square, and the result of each subdivision is consistent;

perlin noise: for a noise image, firstly, integer coordinate points are endowed with normal vectors and gray values in random directions, and besides, the pixel value of any point is determined by the normal vectors of four vertexes in the nearest adjacent cell where the point is located and the direction vector connecting the point and the four vertexes, and the method is provided withFor the vector of points found in the graph, n is the number of noise superpositions, and is generally taken to be n ∈ [3, 12 ]](ii) a Let the wavelength of the noise function be λ, the frequency ω = λ ^-0.5i The gray value of the pointComprises the following steps:

wherein the offset is taken as-0.1 of the offset,is a noise function;

the preprocessing method is adopted to generate various types of noise textures, the additional expense of increasing a GPU (graphics processing unit) in real-time calculation is avoided, in the rasterization stage, a texture graph containing noise is read by a pixel shader and is fused with an elevation graph, meanwhile, the noise graphs are scaled and overlapped in different proportions, infinite noise textures are expanded, and therefore the drawn periodicity and the drawn mode effect are eliminated.

The invention provides a new method and means for solving the above problems, namely, the efficient organization and management of spatial data are performed, and the spatial data are expressed in a digital and visual form, while the visualization of the terrain information researched herein is just the content which is the first and basic content for visualizing the battlefield environment: in application, any battle is carried out in a certain battlefield geographic space; technically, the visualization of the geospatial environment is the basis and the carrier for the visualization superposition display of other elements, and has important significance for researching the three-dimensional expression of military intelligence. The advantages are that:

1. the method analyzes the modeling method and the basic characteristics of the terrain elevation, discusses the current research situation and the hot spot problem in the field, establishes the technical route of the text, analyzes and solves the preprocessing problem of the elevation data, such as the transformation of a plane projection coordinate system, the elimination of noise points in an elevation map and the like, and lays a foundation for organizing and managing the elevation data.

2. The method comprises the steps of establishing a global multi-resolution virtual terrain environment by using a pyramid model, designing and realizing the pyramid model in a global range by combining a multi-resolution LOD model technology, including pyramid establishment rules and a hierarchical partitioning scheme, indexing data by adopting a quadtree structure, and realizing rapid positioning and updating of data nodes by using bit operation and a queue structure.

3. Analyzing and comparing several classical terrain LOD algorithms, and combining a quad-tree structure of a Chunked LOD algorithm and annular addressing of a Geometry Clipmap algorithm to realize global terrain drawing and incremental updating of a rendering model; a more accurate visibility elimination method and a continuous seamless detail transition scheme based on the combination of the bounding sphere and the bounding box are provided, and drawing distortion and model flicker are effectively avoided.

4. The acceleration and the enhancement of terrain rendering are realized on the basis of a programmable GPU technology, namely, vertex information and index information of elevation data are respectively generated into two textures to be stored in a display memory for the whole terrain rendering and calling by using a shader language to control each stage of a graphics rendering pipeline; and interpolating and offsetting the vertex in a geometric stage by utilizing a surface subdivision and fractal technology to generate procedural details, thereby making up the phenomenon of the diagonalization of the time-grid when the resolution is insufficient.

Drawings

FIG. 1 is a UTM projection;

FIG. 2 is a basic block diagram of a multi-resolution pyramid;

FIG. 3 is a diagram of a global multi-resolution terrain layering chunking scheme;

FIG. 4 is a flow chart for filtering noise;

FIG. 5 is a graph of the effect of noise filtering;

FIG. 6 is a diagram of different levels of detail obtained from the same vertex buffer using different step values;

FIG. 7 is a diagram of a dynamic node generation process;

FIG. 8 is a schematic diagram of a subdivision process;

fig. 9 is a determination map of pixel values;

FIG. 10 is a composite of noise maps;

FIG. 11 is a noise map generated by two algorithms;

FIG. 12 is a diagram of the effect of the detail synthesis;

FIG. 13 is elevation primitive data for two sets of data;

FIG. 14 is a bird's eye view of two scenarios;

FIG. 15 is a LOD seamless continuous drawing effect diagram implemented by the present invention;

FIG. 16 is a diagram of real-time detail enhancement effects of rendering;

FIG. 17 is a record chart of the performance of the system monitored by the GPU monitoring tools NVIDIA nTune and GPU-Z;

FIG. 18 is a graph of the number of triangles and corresponding frame rate changes for two sets of experiments when roaming along a fixed route;

FIG. 19 is a graph of geometric error conversion to screen error;

FIG. 20 is a perspective view in perspective projection;

FIG. 21 is a diagram showing the result after vertebral cutting;

FIG. 22 is a diagram of a detail allocation scheme for a level 6 LOD;

FIG. 23 is a geometric transition diagram of a mesh;

FIG. 24 is a multi-level geometric transition diagram.

Detailed Description

1. A preprocessing flow for organizing and managing the topographic data is researched, starting from a plurality of projection coordinate systems adopted for representing models and elevations, implementation methods of projection transformation, noise elimination and interpolation sampling are discussed in sequence, and preparation is made for constructing a multi-resolution topographic data pyramid.

2. The scheme of earth space subdivision and the data organization mode of multi-resolution pyramid layering and blocking are researched, and global multi-resolution terrain data are organized and stored according to a file format designed herein; and (3) constructing indexes for each data node by using a quadtree structure, and providing an update insertion fast algorithm of the quadtree.

3. GPU-based three-dimensional terrain rendering enhancement techniques are studied. The operation and the division of labor of each stage of the graphic drawing assembly line are analyzed, the drawing algorithm combining the terrain data and the GPU programming technology is discussed, the parallel computing capability of the GPU is fully utilized, the cache multiplexing and the real-time detail enhancement are realized, the loads of a CPU and an internal memory are greatly reduced, and the drawing performance and the display effect are improved.

4. Based on the research content, a general design model of the global terrain data management and display system is provided, the functions of all modules are divided, the working process of the system is discussed, and the system is verified by experiments.

The invention discloses a multi-resolution pyramid model construction method for drawing a three-dimensional terrain model, which comprises the following steps:

(1) preprocessing flow of topographic data: the data format is unified through the transformation of projection coordinates, the elimination of image noise points, the construction of a multi-resolution pyramid through sampling interpolation and the like;

i, transformation of projection coordinates: inverse solution to UTM projection:

real-time three-dimensional display systems employ a cartesian right-hand coordinate system to unify models in world space. But considering the establishment of a global geographic reference system, this coordinate system is relatively close to a spatial rectangular coordinate system with the earth's centroid as the origin, and all positional coordinates are in meters, with (x, y, z) representing the offset from the centroid. However, the value of such coordinates is usually very large, and mostly exceeds the representation range of single floating point type data, and the format of the single floating point type data can only store the integer part of 23 bits (the highest bit is used as a positive and negative identifier) and convert the integer part into decimal, namely 2 ²³ -1=8388607, while the fractional part has only 7 bits at most. In the WGS84 reference frame, the equatorial radius of the Earth reaches 6378137m, so with the single floating point format, the highest spatial resolution does not exceed 0.8m (6378137/8388607 ≈ 0.8). At such a threshold, distortion due to various approximate solutions of floating point data, such as vertex overlap and view point jitter, may occur.

Another factor is that the rectangular spatial coordinate system cannot visually represent drawing accuracy information, but on the contrary, many geographical or projection coordinate systems can represent more geographical information, and the geographical coordinate system is related to the adopted reference ellipsoid, such as longitude and latitude geographical coordinates; the projection coordinate system projects the earth model onto some simple expansion curved surfaces, and Lambert structural common (LCC) uses a positive axis equiangular cone projection, and Universal Transform Mercator (UTM) uses a horizontal axis equiangular cylinder projection. The method of projecting geospatial space in a rectangular coordinate system is therefore not advisable, which also requires a positive and negative transformation between the projected coordinate system and the geographical coordinate system to determine the position of the point.

The DEM as a terrain data source adopts WGS84 as a reference ellipsoid, and the adopted coordinate formats include two types, namely UTM coordinates in the unit of long and integral meters or kilometers, WGS84 longitude and latitude coordinates in the unit of double-precision type angles. In order to unify data formats, according to the actual requirements of people, when DEM data are organized, longitude and latitude coordinates are adopted as measurement units of people. Therefore, the DEM data in units of UTM coordinates needs to be converted into longitude and latitude coordinates.

Inverse solution of UTM coordinates, i.e. fromThe projection equation is as follows:

wherein x represents a projected value of latitude, and y represents a projected value of longitude;

let the meridian value of the meridian of the central meridian be L ₀ Corresponding to each weft circle, except for the equator projection which is a straight line, the other weft circles are convex inwards, the intersection point of the X = X straight line and the central meridian is called a bottom point F, and the longitude and latitude coordinates of the bottom point F are (B) _F ，L ₀ ) The longitude and latitude coordinates of the point on the straight line are set as the expression with the point F as the origin:

and further B = B _F +B′，L＝L ₀ + L', the original equation is converted to solve B _F The processes of three parameters of B 'and L'; since the point on the straight line is projected symmetrically relative to the bottom point and the value of x on the same straight line is not changed, B 'is necessarily an even function relative to y, L' is necessarily an odd function relative to y, and the expansion is obtained according to the fourier function:

wherein n is _i (i =0,1,2..) is a function of x, andthat is to say

ByThe following results were obtained:

n ₀ the solution of (2) can be obtained by the algorithm of Gaussian projection, and then each coefficient in the formula can be sequentially solved, and finally, the following results can be obtained:

wherein the respective parameters are represented as follows:

II, eliminating image noise:

by analyzing the noise characteristics in the high-level diagram, it can be found that the noise of the elevation is mainly white noise, the noise area is large, the noise is blocky in the image instead of pure salt noise, and the noise model conforms to the characteristics of additive noise:

wherein, I (x, y) is a gray value at the image (x, y), O (x, y) is a gray value without noise pollution, n (x, y) is a random noise value independent of O (x, y), and p is a probability of noise pollution;

therefore, during filtering, the minimum value is not considered as a noise point, and the initial size of the filtering window is set to be larger so as to avoid missing detection of block noise.

Before processing, the following definitions need to be first specified:

assuming that the size of the image I is M × N, the gray-scale value of any point (x, y) in the image is f (x, y), and the average gray-scale value of the whole image is average (I) _M×N ) The calculation formula of (2) is as follows:

the size of the sliding window is r × r, and if the coordinates of the center pixel of the window are (i, j) and the gray value is f (i, j), the average value of the pixels in the window is average (W) _r×r ) The calculation formula of (c) is:

uniformity is a parameter that measures the contrast of noise signals in an image, image I _M×N And a window W _r×r The uniformity of (c) is defined as follows:

the interpolation method of Lagrange is as follows:

wherein, y _i Represents x _i Outputting a function value; the interpolation algorithm is used for processing one-dimensional data, when the two-dimensional image is expanded, adjacent pixel points in four directions including horizontal, vertical and two diagonal directions of an interpolation pixel point (x, y) are taken for interpolation, and the average value of the interpolation pixel points is taken as the final interpolation result;

the experimental effect of the algorithm is shown in fig. 5, and it can be seen that noise is completely detected and successfully filtered, and although the elevation around the noise is improved by a certain height due to the influence of the noise, the method is acceptable compared with the effect that the whole terrain is smoothed after median filtering.

The process flow of the step of filtering the image is as follows:

setting an initial filtering window size basicWndSize to be 11 × 11, setting a maximum window size maxWndSize to be 19 × 19, setting a noise counter noiseCount =0, and representing whether a pixel point is a noise identifier isNoise = false, calculating uniformity image _ HSV of the whole image, and traversing pixels of the whole image;

(ii) counting the maximum value point pixel value iMax in the current window, the number nMax of the maximum value point pixel values iMax, the number sumPixel of the pixels in the window and the block uniformity wnd _ HSV of the window;

(iii) judging whether the gray value I (x, y) of the current window center pixel point (x, y) is equal to the pixel maximum value iMax in the window or whether the window uniformity wnd _ HSV is greater than the image uniformity image _ HSV; if so, continuously judging whether the number nMax of the maximum value pixels in the current window is smaller than sumPixel/3, if so, marking the central pixel point as a noise point isNoise = TRUE and noiseCount + +; if not, the current extreme value density is considered to be too large, so that the judgment on the noise is influenced, at the moment, the filtering window is expanded by basicWndSize + =2, and the step 2 is skipped; if the filtering window has reached the maximum value, i.e. basicWndSize = = maxwdsisze, then mark the center pixel as noise point isNoise = TRUE, noiseCount + +; otherwise, performing median filtering on the shop; if not, judging that the current point is not the image noise point, jumping out of the loop, and moving the sliding window to the next pixel until the whole image is traversed;

(iv) performing Lagrange interpolation on the four directions determined as the noise points;

III, sampling interpolation: namely, a bilinear interpolation method is adopted for the interpolation of DEM data; the average value is taken as the pixel value of the noise point, and the identifier isNoise of the point is set to false, and meanwhile, noise count-.

(2) Constructing a multi-resolution pyramid:

when the three-dimensional terrain is drawn in real time, in order to ensure the continuity of visual effect, the original terrain grids need to be sampled or interpolated to form grids expressed by different resolutions, and the calling or removal of the corresponding resolution grids is controlled according to the error degree preset by the algorithm.

In order to solve the problem, a multi-resolution pyramid concept is generated, a pyramid structure is a multi-resolution hierarchical model, and in principle, a pyramid is a continuous resolution model, but continuous change of resolution is difficult to realize during construction, and complexity of an algorithm and spatial redundancy are greatly improved, so that when the pyramid is constructed, a multiplying power method is adopted for construction, namely a series of grids are adopted for representing the same region, sampling intervals of upper and lower adjacent grids are different by a same proportion coefficient m (usually m = 2), so that from bottom to top of the pyramid, resolution is reduced layer by layer, and grid size is correspondingly reduced. Let the size of the original grid be S ₀ Resolution is R ₀ The constructed pyramid has l layers (the pyramid top is the 0 th layer), then the grid size s and resolution r of the k-th layer can be expressed as:

partitioning the earth after plane projection according to equal longitude and latitude, and simultaneously constructing pyramid layers with different levels according to a top-down mode, wherein the constructed scheme is as follows as shown in FIG. 3:

the latitude span of the earth is-90 degrees to 90 degrees, the longitude span is-180 degrees to 180 degrees, the width-to-height ratio on a projection plane is just 2;

(II) the ratio of the resolutions of LODs in each layer is 2, namely the resolutions of the sub-blocks are twice of the resolutions of the parent blocks, the sub-blocks have the same geographic space, the pixel size of the sub-blocks is 2 times of the previous level, and the geographic space area covered by the parent blocks is 4 times of that of the sub-blocks under the same screen resolution;

(III) the block sequence of each layer adopts simple row-column coding, the ratio of the number of blocks between the next layer and the previous layer and the ratio of the number of blocks between the horizontal layer and the longitudinal layer of each layer are both 2, according to the habit of a computer image coordinate system, the starting point is defined at the upper left corner, the line number is from top to bottom, and the column number is from left to right; the number of blocks of layer 0 is 4X 2, and the number of blocks of layer k-1 is 2 ^k+1 ×2 ^k ；

(IV) each partition contains DEM data for a corresponding geographic location, which must be odd in resolution size, i.e., (2n + 1) × (2n + 1); the overlapping rate of one pixel with the adjacent block is ensured, so that the occurrence of T-shaped cracks is avoided, the heights of the connecting parts of the adjacent terrain blocks can be consistent when the elevation data is interpolated, and the occurrence of the prism phenomenon is eliminated;

(V) the bottom layer of the pyramid is the highest resolution, and is directly obtained from original data, and the LOD level where the pyramid is located is determined by the resolution and the range of the pyramid; the resolution of texture data and normal data corresponding to each terrain block must be higher than that of DEM (digital elevation model) so as to ensure smooth interpolation and avoid rendering distortion; assuming that the DEM resolution of a certain layer is r, the resolution of the DEM is defined as r' = N x (r-1), N is larger than 1, N belongs to N; a global scale pyramid model is obtained.

The attribute parameters are shown in the following table:

hierarchical chunking results for table pyramid models

Wherein, R represents the earth's major radius defined by WGS 84: 6378137.0m. The global data block is based on the difference of latitude and longitude, and the units in the 4 th and 5 th columns in the table are represented by distance to better conform to the use habit of people, and the values are taken from the longest latitude distance spanned by the block and the ratio of the longest latitude distance to the pixel resolution.

The invention relates to a pipeline drawing enhancement technology based on a GPU: and (3) dynamic image processing:

storing the vertex information and the index information of the elevation separately by using shader language, and transmitting the vertex information and the index information as texture data into a video memory at one time;

a. data structure and organization of nodes:

the world is divided into a number of small regions, each represented by a rectangular solid called a sector, which are formed into a hierarchical tree, each of which is made up of four smaller sub-sectors. When rendering the ground scene, firstly judging which blocks are visible, and if a parent block is invisible, then a child is invisible; otherwise, the child nodes are all visible. When the parent block intersects the view pyramid, its child block is selected and continues to be the parent block to test its visibility. Until all visible blocks are found or the traversal of all nodes is completed.

A blocky quadtree structure is constructed by adopting a Chunked LOD algorithm, namely, a single vertex or a triangle is not used as a processing unit, but a mesh model including a texture map is used for packaging the triangle into strips to complete rendering at one time, and the batch processing capacity of the GPU is fully utilized;

taking a 2D elevation map as an example, filtering the 2D elevation map by using Geomipmap to construct a mipmap pyramid containing L layers. Here we save the data of each layer as a 2D texture at the appropriate size, rather than the buffer that we consider linearizing it to 1D through a vertex array. Each layer of the pyramid comprises n multiplied by n sampling points, and a single-channel elevation map, a three-channel vertex normal map and a three-channel texture map are stored for the sampling points.

Thus, we can define the following node data structure in the shader:

the resolution of the normal map and the texture is twice as high, because the effect of only one normal line at each vertex is fuzzy, and the normal map with twice resolution can be rough after illumination is added, and a softer model can be created through interpolation. The low resolution texture limits the range of the viewpoint, creating large area blurring for each surface type, and the texture resolution should not be less than the resolution of the height map.

b. Multiplexing of the buffer memory: the vertex coordinates (x, y) in the same layer are all transformed in the same way by adopting the method used by Geo-Clipmap;

two classes are established, one is cTerrain, which represents the whole terrain to help cut out the part of the terrain actually needed, and the terrain is cut into the corresponding blocks of the grid. These tiles are represented using a second class, cTerrannSection, which stores the elevation and normal for each tile. For example, for a 256 x 256 vertex terrain, we will create a cTerrain to store the entire data set. The cTerranain object further subdivides the terrain into 32 x 32 vertex tiles, creating a total of 64 cTerranSection objects to represent each tile. Since all the cTerranSection objects are of equal size and equal number of vertices, we can create a single index cache object and use it to render each tile of terrain by shifting the corresponding step value to the right by one for different levels, as shown in FIG. 6.

The code to build the index buffer is as follows:

bool cTerrain::buildIndexBuffer(

int gridVerts,// number of vertices on grid edge

int step)// step value used for index buffering

{ int totalscripts = gridVerts-1; // number of stripes

int index _ tripwise = gridVerts < <1; // each stripe occupies two rows of vertices

// the number of vertices increases due to the presence of 0 triangles

int total_indice＝(indice_stripwise*totalStirps)+(totalStrips<<1)-2；

int*index＝new int[total_indice]；

int vert,start_vert＝0；

for(int j＝0；j<total_strips；++j)

{// generating an index buffer for a line

vert＝start_vert；

for(int k＝0；k<xVerts；++k)

{*(index++)＝vert；

*(index++)＝vert+gridVerts*step；

vert+＝xStep；

}

start_vert+＝lineStep；

if(j+1<total_strips)

{// add 0 triangle at the end of each line

*(index++)＝(vert-xStep)+lineStep；

*(index++)＝start_vert；

}}}

The specific way to construct the vertex buffer is to shrink the position and texture coordinates (x, y) of each block to within (0, 1). During rendering, the coordinates are enlarged to a required scale, and the coordinates can be moved to the corresponding world position by adding a proper offset value. While the individual terrain blocks use a series of similar xy-positions and uv-texture data, we store them once and refer to them when rendering each cTerranAinSection object within the cTerranain object. Thus we build two sets of vertex data streams: the shared vertex data stored in the cTerrain and the independent vertex data of each cTerrain section greatly reduce the occupied memory space.

The vertex buffer code is constructed as follows:

bool cTerrain::create(

cTexture*heightMap,

const cRect3d&worldExtents,

int shift)

{

// setting the scaling ratio

m_tableWidth＝(uint16)heightMap->width()；

m_tableHeight＝(uint16)heightMap->height()；

m_mapScale.x＝m_worldSize.x/m_tableWidth；

m_mapScale.y＝m_worldSize.y/m_tableHeight；

m_mapScale.z＝m_worldSize.z/255.0f；

// set initial offset

m_sectorCountX＝m_tableWidth>>m_sectorShift；

m_sectorCountY＝m_tableHeight>>m_sectorShift；

m_sectorSize.set(m_worldSize.x/m_sectorCountX,

m_worldSize.y/m_sectorCountY)；

cVector2vert(0.0f,0.0f)；

sLocalVertex*pVerts＝

new sLocalVertex[m_sectorVerts*m_sectorVerts]；

Writing shared vertex and texture coordinates into a stream

for(int y＝0；y<m_sectorVerts；++y)

{

vert.set(0.0f,y*cellSize.y)；

for(int x＝0；x<m_sectorVerts；++x)

{

pVerts[(y*m_sectorVerts)+x].xyPosition＝vert；

pVerts[(y*m_sectorVerts)+x].localUV.set(

(float)x/(float)(m_sectorVerts-1),

(float)y/(float)(m_sectorVerts-1))；

vert.x+＝cellSize.x；

}}}

c. Real-time drawing of the three-dimensional terrain model: the method is divided into two types, one is a static node, and the other is a dynamic node;

i, large-scale terrain data rendering: the method comprises the steps that nodes with different resolutions in the same area are obtained through sampling in the preprocessing process, and during rendering, the nodes with the corresponding resolutions are called into a video memory to be directly rendered;

the terrain data of the quad-tree structure tissue of the Chunked LOD is adopted, the nested grid mode in the Geo-Clipmap is combined for drawing, the splitting process of the grid when the LOD changes is adjusted by using continuous geometric transition, and mutation and cracks are avoided.

(1) Bit operation of ring address:

using a grid cache with a fixed size to acquire elevation data of any node in the quadtree, and rendering the whole terrain space by setting an offset address and a geographical space span; before adding visual pyramid clipping, a visual point is taken to be vertically downward, the visual range of the visual point is a pixel grid with n multiplied by n fixed size, and the geographic spatial resolution r represented by each pixel in the layer is obtained according to the pyramid level l of the visual point _l ＝2 ^-l-9 X π × R, where R represents the radius of the earth, in (gx) _diff ，gy _diff ) Representing the displacement of the viewpoint projected on the terrain plane, the distance of the pixel over which the viewpoint moves on the plane is represented asAndthe pixel range covered by the field of view is (x) _min ，y _min ) And (x) _max ，y _max ) Representing;

the pixel range after the viewpoint is moved becomes (x) _min +x _diff ，y _min +y _diff ) And (x) _max +x _diff ，y _max +y _diff ) When updating the node queue, the node-to-queue mapping needs to satisfy the following three conditions:

(1) The same element i always corresponds to the same element f (i) in the queue, namely the mapping relation from i to f (i) meets 1;

(2) F (x) of any element x is corresponding to the element x, so that the continuity of mapping is ensured;

(3) If it isThenThe non-redundancy of the mapping is guaranteed.

The first condition guarantees the possibility of incremental updates, once a node is added to the queue, it does not need to be removed from the queue if it remains in the field of view after the point of view has moved; the second and third conditions prevent data from being overwritten, keeping the queue at a minimum space. And realizing annular address fetching, and realizing the mapping by using modular operation.

The node which moves out of the view field is replaced by the node which enters the view field for the first time, and the rest nodes are fixed, so that increment updating is completed, and for any point (x, y) in the view field range, the position of the node in the node array is (x mod n, y mod n), wherein 'mod' represents modular operation;

the implementation code is as follows:

wherein, (i + x) _min +x _diff ，i+y _min +y _diff ) Actual coordinates of the nodes are represented, (x) ₀ ，y ₀ ) Representing the coordinates of the node in the queue. Therefore, no matter the viewpoint moves towards any direction, the coordinates of the corresponding position can be read through one function, and the complicated step of judging the positive and negative directions of displacement for many times in circulation is omitted.

When the viewpoint position moves by the width of one node to the upper right, the node ABCDEFG enters the active region, the corresponding L-shaped node at the lower left corner is removed, and the bottom node and the side node are respectively replaced in the node cache. When the node is moved upwards to the right by the width of one node, the node HIJKLMN continues to enter the active region, except that the original nodes of the second row and the second column are removed, it should be noted that the nodes a and G are also removed at the same time, the queue array at this time is divided into four regions according to different shading, except that the original node at the upper right corner remains unchanged, the EFMN is a region, the JKDL is a region, and the HIBC also forms a region, and the "L" -shaped update region is changed into a "+" shape in the queue by means of annular addressing, so that reading can be performed according to the three region divisions when reading the update node.

(2) Determination of an error metric:

adopting an error metric calculation mode of Chunked LOD, namely calculating the errors of the nodes layer by layer from leaf nodes and upwards along the quadtree, ensuring that the error metric of a father node is always larger than that of a child node, and determining the error metric e of each node _m Value is max { epsilon ₁ ，ε ₂ ，…，ε _n N is the total number of removed vertices, ε, in the mesh contained by the node _i Representing the difference in height of the terrain after removal of one vertex in the simplification process. And to ensure error nesting, the error metric ε of a node _c Epsilon of its four child nodes _a Is added to its own error metric epsilon _■ ；ε _c The calculation method of (a) is as follows:

the error ε is the geometric error based on three-dimensional object space measurements, and real-time LOD selection may better utilize the two-dimensional screen space error threshold when projecting the geometric error in two-dimensional screen space using perspective projection. After the viewpoint parameters are determined, the corresponding screen error ρ is determined according to the current geometric error ∈, as shown in fig. 19, ω represents the width of the view volume, d represents the distance from the viewpoint to the current layer, θ represents the field angle size, and x represents the screen resolution (in pixels). From similar triangles we can derive:

therefore, when the screen space error is larger than the threshold value, the node with smaller error is selected to be continuously searched downwards until the requirement is met or the leaf node is reached;

(3) visibility elimination:

the selection of the terrain grids is a repeated iteration process, and the distance and the error of each node grid can be continuously and repeatedly calculated, so that a large amount of intersection calculation is involved. In order to reduce the number of the grids to be processed, a judgment mode of combining the bounding box and the bounding sphere is adopted to improve the operation speed. The bounding volume of an object is a shape that can wrap the object inside, and by quickly excluding bounding volumes that do not intersect the view volume, many unnecessary operations can be reduced. If the visual cone is not intersected with the enclosure, the visual cone is not intersected with the object inside the enclosure necessarily, and much time is saved by selecting the enclosure with a simple shape compared with judging the object with a complex internal appearance.

The elimination method when the viewpoint is closer:

when the viewpoint is closer to the terrain, an intersection test is carried out by adopting the bounding box and the view cone, namely, the feature points of the bounding box and 6 planes of the view cone are sequentially intersected to determine whether the feature points belong to the half space of the plane forward direction or are intersected with the planes; this means that the view volume is made up of the positive half of these planes. Thus, if a rectangle exists in the negative half of any viewing cone plane, then the object is outside the viewing cone. The 6 planes of the vertebra are determined by the equation f (x, y, z) = ax + by + cz + dw =0, where the vector (a, b, c) represents the normal vector of the plane and the value dw represents the distance of the plane from the origin; the 9 characteristic points of the bounding box consist of 8 angular points and a central point, the visibility of the terrain grid can be judged by comparing the characteristic points of the bounding box, and the size is equal to half of the length of the diagonal line of the bounding box;

for any point P (x) in space ₀ ，y ₀ ，z ₀ ) In relation to the plane f (x, y, z) = ax + by + cz + dw =0, the formula is calculated:

the absolute value of the result is the vertical distance from the point to the plane, and the absolute value is compared with the size to determine whether the node is inside the view frustum.

The elimination method when the viewpoint is far:

when the viewpoint is far away from the terrain, judging by using a distance coefficient between the node and the viewpoint, and defining the size of the distance coefficient as the actual distance represented by the node bounding box, wherein the ratio of the length L of a projection plane to the distance D of the bounding box from the projection plane is defined, namely a = L/D;

taking the coordinates v (x) of the viewpoint ₀ ，y ₀ ，z ₀ ) The bounding box center coordinate is c (x) ₁ ，y ₁ ，z ₁ ) If the direction vector from the viewpoint to the ground is d (m, n, p) = v-c;

let 8 vertex coordinates of bounding box be p _i (x _i ，y _i ，z _i ) Passing point p _i The equation of a straight line parallel to d is:

f(x，y，z)＝m(x-x _i )+n(y-y _i )+p(z-Z _i )＝0 (3.9)

the near shearing surface of the visual cone body is taken as a projection plane, is vertical to the direction vector d, and has a plane equation of

mx+ny+pz＝0 (3.10)

By combining the above two equations to obtain p _i Coordinates p of the projected point on the plane _i ' is

Calculate p for all 8 vertices _i ' after, 8 are takenMaximum and minimum values min X, max X, min Y, max Y, min Z, max Z in the coordinate points; order to

The value of the distance coefficient a is:

when the viewpoint height is larger than a set threshold value, if the a value of the node bounding box is smaller than 1, namely the size of the node bounding box is less than one pixel on a projection plane, the node is considered invisible and is not rendered any more, otherwise, the node is rendered continuously;

(4) gap elimination and geometric transition:

a new LOD level conversion control technology is provided, and meanwhile, degraded triangle insertion operation among grids with different resolutions is eliminated, so that unnecessary rendering burden and rendering distortion caused by the unnecessary rendering burden are avoided, and seamless LOD transition is realized.

In order to control the selection of a proper level of detail at a certain distance, before node screening, the geographical range covered by each LOD level is calculated, so that the quantity of triangles with approximately consistent granularity is generated on a window plane during rendering; due to the definition of the quadtree structure, the ratio of the detail of each layer to the number of the top points of the previous layer is fixed to be 4, so the distance multiple in the one-dimensional direction is also 2, and therefore the transition coefficients of different resolutions are set to be about 2 to ensure that the number of triangles is approximately consistent. Fig. 22 shows the distance ranges covered by different LOD levels, as the resolution is reduced, the geographic range occupied by each terrain block is 2 times that of the upper level, and the last level (farthest from the viewpoint and lowest level of detail) represents the visible range of the whole viewpoint.

In order to ensure smooth transition between adjacent LOD levels and keep the terrain precision as much as possible, a transition scheme consistent with a Geo-Clipmap algorithm is adopted, namely, a part of a region adjacent to a low-resolution grid and a high-resolution grid is marked out to be used as a transition region, transition is carried out from low to high, and 33% of the range of each LOD is selected as the transition region in an experiment. The side effect of this is that more triangles are produced, but unlike the method of generating a degenerate triangle, the triangle generated by the vertex displacement method is much simpler in complexity.

The vertexes of all transition areas are transformed according to the unique error measurement, namely, the transition is carried out point by point, and each vertex supports two transformation modes: the transition from itself to higher resolution and from itself to lower resolution. Taking the second transformation as an example, the transition operation is that each vertex mesh of 3 × 3 is composed of 8 triangles, and the vertex mesh of 2 × 2 is simplified by gradually moving the positions of 4 side midpoints and 1 center point to the lower left corner point, and is composed of 2 triangles. The whole process looks like a process in which two triangles in the upper right corner are gradually enlarged, and the remaining 6 triangles are compressed little by little in the lower left direction into degenerate triangles and are finally no longer rasterized into pixels. This method can effectively eliminate gaps between adjacent LOD levels and produce smooth transitions without T-junctions.

Firstly, the approximate distance from a vertex to a viewpoint is calculated to determine the size of a displacement coefficient, and the vertex position for calculating the distance is approximately obtained according to the position proportion of the grids in the world coordinate system. In the experiment, x and y are taken from longitude and latitude values, respectively, and elevation values are obtained by continuously sampling an elevation map. From FIG. 22, let the starting distance of LOD of each layer be start _i I represents the level of the current LOD, the range of the covered distance is range, the distance from the viewpoint to any vertex of the current LOD is dist, and then the current transition region [ mStart, mEnd ] is]And the transition coefficient k of the vertex is:

wherein clamp (x, a, b) is a clamping function, i.e.

Thus, the size of k is limited to between 0 (no transformation at all, staying at the low resolution position) and 1 (transformation at all, the number of triangles becomes twice as large), and the mesh also achieves a transition from low resolution to high resolution as the value of k gradually increases from 0 to 1. After obtaining the transition coefficient, defining the size of the mesh where the vertex is located as n × n, the depth of the quadtree where the vertex is located as l, and the normalized coordinate of the vertex in the mesh as (x, y), and then the actual distance scaling value scale of the mesh size relative to the world coordinate and the coordinate (x ', y') after the transition are obtained by the following formula:

wherein frac (x) represents the fractional part of x;

finally, the vertex coordinates are transformed from (x, y) to the position of (x/2, y/2), and the corresponding elevation is obtained through bilinear interpolation; when all the vertex transitions are completed, the current mesh is completely overlapped with the high-resolution mesh of the next layer, and the process in the period does not generate any gaps. Fig. 23 shows the grid shapes corresponding to different transition coefficients.

In fact, since the viewpoint is continuously changed, and the resolution of the selected nodes gradually decreases as the view deepens in the whole view field, more than one layer of geometric transition occurs in each frame, as shown in fig. 24, a multi-layer transition zone can tightly connect adjacent resolution grids, and the generation of cracks is avoided. Experiments show that when terrain roaming is carried out, the smooth and continuous transition method among grids effectively eliminates 'sudden jump feeling' caused by node details, the frame rate is stabilized above 60 frames, and visual stagnation is avoided.

And II, dynamic nodes which are nodes generated by the GPU in real time through calculation and added into the quadtree and still cannot meet the requirement of a viewpoint threshold when the quadtree structure extends to leaf nodes.

And (3) generation of dynamic nodes:

subdivision of the terrain model: the subdivision is that face points and edge points are inserted on original vertexes of the quadrilateral grids according to subdivision masks, and the face points and the edge points are connected with the original vertexes to form a new quadrilateral, the degree of the generated non-boundary vertexes is 4, the degree of the boundary vertexes is 3, and the non-boundary vertexes are regular points; as shown in fig. 8, the regular point subdivision process: to quadrangle V ₅ V ₆ V ₇ V ₈ When the fine division is performed, four edge points E are formed on four edges of the fine division ₀ E ₁ E ₂ E ₃ Meanwhile, a dough point F is also provided, so that the original dough sheet is split into four sub-dough sheets; e ₀ ，E ₁ ，E ₂ ，E ₃ The value of (c) is as follows:

the value of the face point F is given by equation 4.2:

wherein the content of the first and second substances,

the subdivision algorithm can only be applied to a regular grid where regular points are located, and singular points are not considered. Due to the characteristics of the round LOD algorithm, the thinning process is independently carried out in each block, the grids in the adjacent blocks are not known, when a mask moves to the edge of each block, the traditional edge extrapolation and diagonal point replacement can generate false edges, in order to guarantee the visual reality, edge sampling points of the adjacent blocks are stored in the block structure, the edge sampling points are used as control points only when the opposite edges are thinned, the subdivision in the blocks is not considered, and the continuity of the edges is obtained at the cost of adding a small amount of data.

Random fractal: representing the coordinates of a point on the surface by a variable (X, y), wherein X (X, y) represents the height of the point on the surface; the FBM surface is defined as a random process with an index of H (0 < H < 1) in a certain probability space

A random variable X: r → R ² Satisfies the following conditions:

(1) With probability 1, X (0, 0) =0, i.e. the process starts from the origin; and X (X, y) is a continuous function with respect to (X, y);

(2) For any variable (x, y), (h, k) e R ² With a two-dimensional increment X (X + h, y + k) obeying the expectation of 0 and a variance of (h) ² -k ² ) ^H Normal distribution of which the probability satisfies

The three-dimensional terrain drawn by using the fractal technology can generate infinite vivid scenes with self-similarity, and is easy to smooth and interpolate. Two fractal algorithms are used here to generate different patterns of noise maps, the Diamond-Square algorithm and Perlin noise, respectively.

Diamond-Square algorithm

The essence of the algorithm is a random midpoint displacement algorithm on a two-dimensional plane. The method starts from a square formed by seed points, and after a plurality of iterations of a midpoint displacement method, the seed square is continuously subdivided to obtain a vivid three-dimensional terrain. Two steps, diamond and Square, are performed for each subdivision in the iterative process.

And the Diamond step is used for generating the center of the square, and four vertexes are taken to calculate the average value, and a random quantity is added to serve as the value of the central point. Four triangular pyramids are generated by the Diamond step. The Square step is to take the four vertices of the pyramid and use the average value plus the random variable (subject to the same probability distribution as the random variable in the Diamond step) as the midpoints of the four sides of the original Square, and the result of each subdivision is the same as that in fig. 8.

At each stage of the process, random values are always generated within a fixed range, we call the current range delta and the total range [ -delta, delta ]. In the next subdivision, delta is multiplied by a scaling value, called roughness (roughness), which is between 0 and 1, to reduce the range of delta values. High roughness values produce a more disordered topography, while low roughness results in an increasingly flat topography.

The Diamond-Square algorithm is simple to implement, high in operation speed and widely adopted. But because there is no information transmission between adjacent squares, it causes "crease problem" in the noise image, i.e. there is an unnatural linear track in the image and it cannot be removed by local smoothing; in addition, the iterative generation mode has large calculation amount, consumes excessive time and is not flexible.

Perlin noise

Perlin noise is a function for generating random noise, the Perlin noise obtains a final result through synthesis of multi-fractal, namely successive increase synthesis by using limited bandwidth frequency, and calculation of each frequency band is independent in the synthesis process, so that the Perlin noise has strong control capability on frequency, dimension and other corresponding statistical characteristics of terrain fractal.

For a noise image, firstly, a normal vector and a gray value in random directions are given to an integer coordinate point, and besides, the pixel value of any point is determined by the normal vectors of four vertexes in a nearest cell where the point is located and a direction vector connecting the point and the four vertexes.

For a noise image, firstly, integer coordinate points are endowed with normal vectors and gray values in random directions, and besides, the pixel value of any point is determined by the normal vectors of four vertexes in a nearest adjacent cell where the point is located and a direction vector connecting the point and the four vertexes.

Is provided withFor the vector of points found in the figure, n is the number of noise superpositions, and is generally taken to be n ∈ [3, 12 ]]. Let the noise function wavelength be λ, then frequencyRate ω = λ ^-0.5i The gray value of the pointComprises the following steps:

wherein the offset is taken as-0.1 of the offset,is a noise function. The lower graph shows the appearance of two independent noise sums.

The preprocessing method is adopted to generate various types of noise textures, the extra expense of increasing a GPU (graphics processing unit) in real-time calculation is avoided, in the rasterization stage, a texture map containing noise is read by a pixel shader and is fused with an elevation map, and the execution effect is quicker. Meanwhile, the noise graph is scaled and superposed in different proportions, so that infinite noise textures can be expanded, and the drawn periodicity and modeling effects are eliminated. Fig. 11 shows the texture maps of the noise generated by two noise algorithms, and fig. 12 shows the effect of superimposing the two noises with elevation, respectively. As can be seen from the figure, the two can enhance the process details in real time and simulate different terrain effects, the Diamond-Square noise is closer to the sand effect, and the Perlin noise is closer to the rock effect.

The method has the following effects:

two sets of experimental data respectively construct three-dimensional terrains of world elevations of 8-level LOD and Puget Sound of 12-level LOD, wherein the resolution of the 8-level world elevations is 152m, the whole earth space within 83 degrees of south-north latitude in a geographic range is covered, the constructed original grid consists of 268,435 and 456 triangles, the elevation resolution of the 12-level Puget Sound is 10m, the geographic range is 163km multiplied by 163km, and the constructed original grid consists of 463,601 and 664 triangles. Because the data format of the elevation information is not compressed, the size of the 12-level LOD file reaches 2.58GB, and the corresponding normal map approaches 4GB.

The following table shows the change in system performance before and after the start-up detail enhancement. The data shows that as the viewpoint is lowered, the LOD level is also progressively deeper and the number of triangles is increasing, but at the same time the blocks that intersect the view volume are also decreasing, so the number of triangles for a level 8 LOD is suddenly reduced at a lower height. Levels 9 and 10 are the result of dynamic node joining, and the number of triangles starts to increase, but the loading speed also slows down. Similarly, the first 8 levels are static nodes generated by a preprocessing process, and when browsing, only the corresponding terrain blocks need to be read from the storage, the frame rate is stable, while the insertion of the dynamic nodes brings a certain burden to real-time rendering, and the change of the frame rate depends on the terrain complexity and the number of triangles.

Runtime comparison of static and dynamic nodes of a table

Viewpoint height (m)	LOD level where viewpoint is located	Number of triangles	Loading time(s)	Roaming frame rate(s) -1 )
					1846.2	5	45056	8.96	145.22
1282.8	6	40960	7.82	147.63
					1015.5	7	32768	7.46	146.97
882.2	8	6144	5.88	159.46
					424.7	9	30720	7.52	133.23
301.1	10	26624	7.33	140.60

As can be seen from the two sets of monitoring data in fig. 17, the GPU bears most of the operation load, the occupancy rates of the CPU and the memory are always at a lower level, and the numbers of nodes accommodated in the field of view range are also different due to the different depths of the quadtrees used by the two sets of data, which is a main cause of inconsistent video memory occupancy, but compared with the GB-level size of the original file, the system has implemented a function of browsing a large-scale terrain using a small amount of cache.

Counting the obtained data, the system averages the number of triangles drawn by each frame to 265.6k at the running time, and the frame rate approaches 194s ^-1 It can be seen from fig. 18 that the efficiency of drawing is improved significantly by using the GPU, and the frame rate and the number of triangles are not strictly negative correlation, and as the depth of the quadtree is deepened, the details of the terrain are more and more abundant, which leads to a significant decrease in the frame rate, i.e., the frame rate is greatly affected by the complexity of the terrain, and the 8-level world elevation is mostly flat terrain, so the frame rate is less affected by the number of images of the triangles, and the 12-level Puget Sound fluctuates greatly. The following table lists frame rate comparisons of some algorithms, but because different algorithms adopt different hardware platforms, and the algorithms in the text have strong dependency on the GPU, the horizontal comparison is difficult to perform, and the rendering hardware load caused by the algorithms should be reduced as much as possible in the next step of research.

Rendering efficiency comparison of different algorithms of a table

Algorithm	Average frame rate(s) -1 )
		Classical algorithm	10-20
ChunkedLOD	21-50
		Geo-Clipmap(2004)	95-120
Geo-Clipmap(2005)	130-185
		CDLOD	56-1800
Text algorithm	82-280

The invention designs and realizes a terrain multi-resolution display verification system. Firstly, a design scheme of the system is introduced from two aspects of a frame structure and functions of each module, then, two groups of experimental data are used for verifying main functions of the system and are displayed by using a chart, and finally, performance tests are performed by professional detection software, so that the system is proved to meet the real-time and GPU-friendly drawing requirements.

Claims (1)

1. A method for smoothly drawing a three-dimensional terrain model in real time by combining a GPU technology is characterized by comprising the following steps:

constructing a multi-resolution pyramid model of the three-dimensional terrain model:

(1) preprocessing flow of topographic data: filtering the original image to obtain a noise-free image, and constructing a multi-resolution pyramid;

i, eliminating image noise:

the noise model conforms to the characteristics of additive noise:

wherein, I (x, y) is a gray value at the image (x, y), O (x, y) is a gray value without noise pollution, n (x, y) is a random noise value independent of O (x, y), and p is a probability of noise pollution;

the process of filtering the image comprises the following steps:

calculating the uniformity image _ HSV of the whole image, and traversing pixels of the whole image;

uniformity is a parameter that measures the contrast of noise signals in an image, image I _M×N And a window W _r×r The uniformity of (c) is defined as follows:

assuming that the size of the image I is M × N, the gray-scale value of any point (x, y) in the image is f (x, y), and the average gray-scale value of the whole image is average (I) _M×N ) The calculation formula of (2) is as follows:

the size of the sliding window is r × r, and if the coordinates of the center pixel of the window are (i, j) and the gray value is f (i, j), the average value of the pixels in the window is average (W) _r×r ) The calculation formula of (2) is as follows:

(ii) counting the pixel value iMax of the maximum value point in the current window, the number nMax of the pixel value iMax, the number sumPexel of the pixels in the window and the block uniformity wnd _ HSV of the window;

(iii) judging whether the gray value I (x, y) of the current window center pixel point (x, y) is equal to the pixel maximum value iMax in the window or whether the window uniformity wnd _ HSV is greater than the image uniformity image _ HSV;

(iv) performing Lagrange interpolation on the four directions which are judged as the noise points; the Lagrange interpolation method comprises the following steps:

wherein, y _i Represents x _i The function value of (c); the interpolation algorithm is used for processing one-dimensional data, when the interpolation algorithm is expanded to a two-dimensional image, adjacent pixel points in four directions including horizontal, vertical and two diagonal directions of an interpolation pixel point (x, y) are taken for interpolation, and the average value of the adjacent pixel points is taken as the final interpolation result;

(2) constructing a multi-resolution pyramid by the data after the noise points are eliminated:

the pyramid structure is a multi-resolution hierarchical model and is constructed by adopting a multiplying power method, namely for the same region, a series of grids are adopted for representation, and the sampling intervals of upper and lower adjacent grids have the same proportional coefficient m, so that the resolution is reduced layer by layer from the bottom to the top of the pyramid, and the grid size is correspondingly reduced;

let the size of the original grid be S ₀ Resolution is R ₀ If the constructed pyramid has l layers and the pyramid top is the 0 th layer, then the grid size s and resolution r of the k-th layer can be expressed as:

wherein m is more than 1;

partitioning the earth after plane projection according to equal longitude and latitude, and simultaneously constructing pyramid layers of different levels according to a top-down mode, wherein the construction scheme is as follows:

the latitude span of the earth is-90 degrees to 90 degrees, the longitude span is-180 degrees to 180 degrees, the width-to-height ratio on a projection plane is just 2;

(II) the ratio of the resolutions of LOD of each layer is 2, namely the resolution of the sub-block is twice that of the parent block, the pixel size of the sub-block is 2 times that of the previous level in the same geographic space, and the geographic space area covered by the parent block is 4 times that of the sub-block in the same screen resolution;

(III) the block sequence of each layer adopts simple row-column coding, the ratio of the number of blocks between the next layer and the previous layer and the ratio of the number of blocks between the horizontal layer and the longitudinal layer of each layer are both 2, according to the habit of a computer image coordinate system, the starting point is defined at the upper left corner, the line number is from top to bottom, and the column number is from left to right; the number of blocks of layer 0 is 4X 2, and the number of blocks of layer k-1 is 2 ^k+1 ×2 ^k ；

(IV) each tile contains DEM data for the corresponding geographic location, which must be odd in resolution size, i.e., (2n + 1) × (2n + 1);

(V) the bottom layer of the pyramid is the highest resolution, and is directly obtained from original data, and the LOD level where the pyramid is located is determined by the resolution and the range of the pyramid; the resolution of texture data and normal data corresponding to each terrain block must be higher than that of DEM (digital elevation model) so as to ensure smooth interpolation and avoid rendering distortion; assuming that the DEM resolution of a certain layer is r, the resolution of the DEM is defined as r' = N x (r-1), N is larger than l, and N belongs to N; obtaining a global scale pyramid model;

and (3) based on a GPU pipeline drawing enhancement technology, drawing the data of the global pyramid model:

storing the vertex information and the index information of the elevation separately by using shader language, and transmitting the vertex information and the index information as texture data into a video memory at one time;

a. data structure and organization of nodes:

a blocky quadtree structure is constructed by adopting a ChunkedLOD algorithm, namely, a single vertex or a triangle is not used as a processing unit, but a mesh model including a texture map is used for packaging the triangle into strips to complete rendering at one time, and the batch processing capacity of the GPU is fully utilized;

b. multiplexing of the buffer memory: the vertex coordinates (x, y) in the same layer are all transformed in the same way by adopting the method used by Geo-Clipmap;

c. and (3) data of the pyramid model are drawn in real time: the method is divided into two types, one is a static node, and the other is a dynamic node;

i, static node data rendering: during rendering, the nodes with the corresponding resolution are transferred into a video memory to be directly rendered;

(1) bit operation of ring address:

using a grid cache with a fixed size to acquire elevation data of any node in the quadtree, and rendering the whole terrain space by setting an offset address and a geographical space span; before adding visual pyramid clipping, a visual point is taken to be vertically downward, the visual range of the visual point is a pixel grid with n multiplied by n fixed size, and the geographic spatial resolution r represented by each pixel in the layer is obtained according to the pyramid level l of the visual point _l ＝2 ^-l-9 X π × R, where R represents the radius of the earth, in (gx) _diff ，gy _diff ) Representing the displacement of the viewpoint projected on the terrain plane, the distance of the pixels where the viewpoint moves on the plane is represented asAndthe pixel range covered by the field of view is (x) _min ，y _min ) And (x) _max ，y _max ) Representing;

the node which is moved out of the view field is replaced by the node which enters the view field for the first time, and other nodes are still, so that incremental updating is completed, and for any point (x, y) in the view field range, the position of the node in the node array is (xmodn, ymodn), wherein 'mod' represents modulo operation;

(2) determination of the error metric:

adopting an error metric calculation mode of ChunkedLOD, namely calculating the errors of the nodes layer by layer from leaf nodes along the quadtree, and calculating the error metric epsilon of one node _c Epsilon of its four child nodes _ci Plus its own error metric epsilon _m ；ε _c The calculation of (c) is as follows:

therefore, when the screen space error is larger than the threshold value, the node with smaller error is selected to be continuously searched downwards until the requirement is met or the leaf node is reached;

(3) visibility elimination:

the elimination method when the viewpoint is closer:

when the viewpoint is closer to the terrain, an intersection test is carried out by adopting the bounding box and the view cone, namely, the feature points of the bounding box and 6 planes of the view cone are sequentially intersected to determine whether the feature points belong to the half space of the plane forward direction or are intersected with the planes; the 6 planes of the view cone are determined by the equation f (x, y, z) = ax + by + cz + dw =0, where the vector (a, b, c) represents the normal vector of the plane and the value dw represents the distance of the plane from the origin; the 9 characteristic points of the bounding box consist of 8 angular points and a central point, the visibility of the terrain grid can be judged by comparing the characteristic points of the bounding box, and the size is equal to half of the length of the diagonal line of the bounding box;

the elimination method when the viewpoint is far:

when the viewpoint is far away from the terrain, judging by using a distance coefficient between the node and the viewpoint, and defining the size of the distance coefficient as the actual distance represented by the node bounding box, wherein the ratio of the length L of a projection plane to the distance D of the bounding box from the projection plane is defined, namely a = L/D;

taking the coordinates v (x) of the viewpoint ₀ ，y ₀ ，z ₀ ) The bounding box center coordinate is c (x) ₁ ，y ₁ ，z ₁ ) If the direction vector from the viewpoint to the ground is d (m, n, p) = v-c;

let 8 vertex coordinates of bounding box be p _i (x _i ，y _i ，z _i ) Passing point P _i The equation of a straight line parallel to d is:

f(x，y，z)＝m(x-x _i )+n(y-y _i )+p(z-z _i )＝0 (3.9)

the near shearing surface of the visual cone body is taken as a projection plane which is vertical to a direction vector d, and the plane equation is set as

mx+ny+pz＝0 (3.10)

By combining the above two equations, P can be obtained _i Coordinates P of the projected point on the plane _i ' is

Calculate P 'for all 8 vertices' _i Taking maximum and minimum values minX, maxX, minY, maxY, minZ and maxZ in the 8 coordinate points; order to

x＝maxX-minX

y＝maxY-minY

z＝maxZ-minZ (3.12)

The value of the distance coefficient a is:

when the viewpoint height is larger than a set threshold value, if the a value of the node bounding box is smaller than 1, namely the size of the node bounding box is less than one pixel on a projection plane, the node is considered invisible and is not rendered any more, otherwise, the node is rendered continuously;

(4) gap elimination and geometric transition:

before node screening, firstly, calculating the geographical range covered by each LOD level so as to generate the quantity of triangles with approximately consistent granularity on a window plane during rendering;

let the starting distance of LOD of each layer be start _i I represents the level of the current LOD, the range of the covered distance is range, and the distance from the viewpoint to any vertex of the current LOD is dist, then the current transition is realizedRegion [ mStart, mEnd]And the transition coefficient k of the vertex is:

wherein clamp (x, a, b) is a clamping function, i.e.

After obtaining the transition coefficient, defining the size of the grid where the vertex is located as n × n, the depth of the quadtree where the vertex is located as l, and the normalized coordinate of the vertex in the grid as (x, y), then the actual distance scaling value scale of the grid size relative to the world coordinate and the coordinate after transition (x ', y') are obtained by the following formula:

where frac (x) represents the fractional part of x;

finally, the vertex coordinates are transformed from (x, y) to the position of (x/2, y/2), and the corresponding elevation is obtained by bilinear interpolation;

II, dynamic node: when the quad-tree structure extends to the leaf nodes, the requirement of the viewpoint threshold value cannot be met, and the nodes are generated by the GPU in real time calculation and added into the quad-tree;

and (3) generation of dynamic nodes:

(1) subdivision of the terrain model: the subdivision is that face points and edge points are inserted on original vertexes of the quadrilateral grids according to subdivision masks, and the face points and the edge points are connected with the original vertexes to form a new quadrilateral, the degree of the generated non-boundary vertexes is 4, the degree of the boundary vertexes is 3, and the non-boundary vertexes are regular points;

the regular point subdivision process: to quadrangle V ₅ V ₆ V ₉ V ₁₀ When the fine division is performed, four edge points E are formed on four edges of the fine division ₀ E ₁ E ₂ E ₃ Meanwhile, a dough point F is also provided, so that the original dough sheet is split into four sub-dough sheets;E ₀ ，E ₁ ，E ₂ ，E ₃ the value of (a) is as follows:

the value of the face point F is given by equation 4.2:

wherein the content of the first and second substances,

(2) random fractal: the coordinates of a point on the surface are expressed by variable (X, y), and X (X, y) represents the height of the point on the surface; the FBM curved surface is defined as a random process with index H on a certain probability space, and 0-H-s-1

A random variable X: r → R ² Satisfies the following conditions:

(1) With probability 1, X (0, 0) =0, i.e. the process starts from the origin; and X (X, y) is a continuous function with respect to (X, y);

(2) For any variable (x, y), (h, k) e R ² With a two-dimensional increment X (X + h, y + k) obeying the expectation of 0 and a variance of (h) ² -k ² ) ^H Is normally distributed with a probability satisfying

Two fractal algorithms were used to generate different patterns of noise maps, the Diamond-Square algorithm and Perlin noise:

Diamond-Square algorithm: in the iteration process, each subdivision can execute two steps of operation, namely a Diamond step and a Square step; the Diamond step is used for generating the center of a square, taking four vertexes to calculate an average value, and adding a random quantity to serve as a value of a central point; four triangular pyramids are generated by the Diamond step; the step of Square is to take four vertexes of the pyramid, use the average value of the four vertexes plus random variables as the middle points of the four edges of the original Square, and the result of each subdivision is consistent;

perlin noise: for a noise image, firstly, a normal vector and a gray value in random directions are given to an integer coordinate point, except that the pixel value of any point is determined by the normal vectors of four vertexes in a nearest adjacent cell where the point is located and a direction vector connecting the point and the four vertexes,

is provided withIs the vector of the points in the graph, n is the superposition times of the noise, and n is taken as [3, 12 ]](ii) a Let the wavelength of the noise function be λ, the frequency ω = λ ^-0.5i The gray value of the pointComprises the following steps:

wherein the offset is offset by-0.1,is a noise function;

in the rasterization stage, a pixel shader reads a texture graph containing noise and an elevation graph to be fused, and meanwhile, the noise graph is scaled and superposed in different proportions to expand infinite noise textures, so that the drawn periodicity and the drawn mode effect are eliminated.