CN105336003A - Three-dimensional terrain model real-time smooth drawing method with combination of GPU technology - Google Patents

Abstract

The invention provides a three-dimensional terrain model real-time smooth drawing method with combination of a GPU technology, and belongs to the technical field of image processing. The objective of the invention is to provide the three-dimensional terrain model real-time smooth drawing method with combination of the GPU technology so that cache reuse in multiple times of drawing can be realized based on the current popular programmable GPU technology with a global digital elevation model acting as a data source, and load of computation space is effectively reduced. The method comprises the steps of construction of a multi-resolution pyramid model, elimination of image noise points, filtering of images, partitioning of planar projection of the earth according to equal latitude and longitude, and construction of different hierarchical levels of pyramid layers according to a mode from the top to the bottom. Acceleration and enhancement of terrain rendering are realized based on the programmable GPU technology, i.e. all phases of a graphical drawing pipeline are controlled by using shader languages, two and textures are respectively generated by vertex information and index information of elevation data to be stored in video memory for scheduling of whole terrain drawing; and vertex interpolation and migration are performed in the geometric phase by utilizing a curved surface subdivision and fractal technology so that procedural details are generated and the phenomenon of edges and corners of the terrain mesh when resolution is insufficient can be compensated.

Description

The method of three-dimensional terrain model is drawn out in conjunction with the real-time smoothness of GPU technology

Technical field

The invention belongs to image processing technique field.

Background technology

Along with Aero-Space science and technology and remote sensing technology increasingly mature, the spatial resolution being looped around sensor entrained by the reconnaissance satellite of near-earth orbit and various aircraft is improving constantly, digital photogrammetry technology also more and more meticulousr accurately.Multiple remotely-sensed data the data set that converges be the important information source of geospatial intelligence, remote sensing image wherein and digital elevation are then basis and the foundation of whole digital city.To under the burgeoning prerequisite of the application demand of this geo-spatial data, the data volume of satellite remote-sensing image and digital elevation is also and then day by day huge, how effective management and use are carried out to it, become problem in the urgent need to address in Remote Sensing Data Processing and Geographic Information System field, the particularly sensor information of release management and organizational scheduling magnanimity, has extremely important realistic meaning for the various application demands such as territory planning, Mapping remote sensing technology, military traffic, hydraulic and electric engineering, agricultural environmental protection provide in geographic information services.

The development of modern military technology, constantly advances the direction of space, battlefield towards various dimensions, many key elements.In the face of complicated staggered, fast changing battlefield surroundings and situation, commander how is assisted to dispose operation all sidedly, grasp important goal and the content that battle field information has just become military information exactly.Meanwhile, our battle space is in continuous expansion, and the development of geospatial intelligence shows slightly delayed, in early days due to the development restriction of sensor technology and computer graphics techniques, the acquisition of data and visualization capability are difficult to reach requirement accurately in real time, the processing mode of two dimensional surface can simplify the process that spatial information is understood and expressed, but to sacrifice the authenticity and integrity of spatial information, and especially elevation information and three-dimensional topological relation are cost.But making rapid progress when remote sensing science and technology and computer graphics, make the high-resolution elevation of acquisition and image, 3-D display efficiently true to nature is no longer the technological challenge being difficult to realize, and the development of ambit also all impels military information must make up the technological deficiency of its three-dimensional expression scarce capacity.

In military field, the demand for development commander of the form of war can not be confined on the information of a certain fixed area or certain particular form the grasp of battlefield surroundings again, want to obtain the triumph of local war under Modern Information based condition, just must control battlefield flexibly, in all directions in time, and can the accurate clear lot of essential factors representing battlefield surroundings all sidedly in real time, for officers and men carry out aid decision making, it is the severe challenge that military information faces.

Summary of the invention

The object of the invention is to take global digital elevation model as data source, based on the Programmable GPU technology of current hot topic, the buffer memory achieved when repeatedly drawing is multiplexing, efficiently reduces the method drawing out three-dimensional terrain model in conjunction with the real-time smoothness of GPU technology of computer memory load.

Step of the present invention is:

The multi-resolution pyramid model construction of three-dimensional terrain model:

1. the pretreatment process of terrain data: filtering is carried out to raw video and obtains muting image, and build multi-resolution pyramid;

I, the elimination of image noise:

Noise model meets the feature of additive noise:

(2.12)

Wherein, for being positioned at image the gray-scale value at place, for the gray-scale value not by noise pollution, be independent of random noise value, for being subject to the probability of noise pollution;

Steps flow chart image being carried out to filtering is as follows:

(I), calculate the uniformity coefficient image_HSV of entire image, then the pixel of entire image is traveled through;

Uniformity coefficient weighs the parameter of noise signal contrast in image, image i _{m× n} and window w _{r× r} uniformity coefficient be defined as follows respectively:

(2.15)

Suppose image isize is m× n, any point in image gray-scale value be , the average gray value average of entire image ( i _{m× n} ) computing formula be:

(2.13)

The size of moving window is r× rif the coordinate of window center pixel is , gray-scale value is , then pixel average average in this window ( w _{r× r} ) computing formula be:

(2.14)

Maximum point pixel value iMax in (II), statistics current window, and number nMax, the number of pixels sumPixel in window, and the block uniformity wnd_HSV of window;

(III), judge current window central pixel point gray-scale value whether equal pixel maximum value iMax in window, or whether window uniformity coefficient wnd_HSV is greater than image uniformity coefficient image_HSV;

(IV), to being judged to be that the four direction of noise spot carries out Lagrange's interpolation; The interpolation method of Lagrange is:

(2.16)

Wherein, represent the functional value gone out; Herein, interpolation algorithm processes for one-dimensional data, when being extended on bidimensional image, gets interpolating pixel point level, vertical, two diagonals adjacent pixels point altogether on four directions carries out interpolation, and gets its mean value as last interpolation result;

2. by the data construct multi-resolution pyramid after above-mentioned removing noise:

Pyramid structure is a kind of multiresolution hierarchical model, and adopt multiplying power method to build, namely to areal, adopt a series of grid representation, the sampling interval of neighbouring grid differs identical scale-up factor m, bottom-up from pyramid like this, resolution successively reduces, and size of mesh opening is corresponding reduction also;

If original mesh is of a size of s ₀, resolution is r ₀, constructed pyramid has llayer, pyramid top is the 0th layer, so the kthe size of mesh opening of layer sand resolution rcan be expressed as:

(2.7)

Carry out piecemeal to the earth after plane projection according to waiting longitude and latitude, simultaneously according to top-down mode, build the pyramidal layer of different levels rank, the scheme of structure is as follows:

(I) the latitude span of the earth is-90 ° to 90 °, longitude span is-180 ° to 180 °, the ratio of wide height is on a projection plane just 2:1, therefore quaternary tree is divided into 8 root nodes the 0th grade time, laterally equal 90 ° of longitudinal span, with the satisfied requirement waiting longitude and latitude to divide, it is invalid that the coordinate figure except this scope is all considered as;

(II) ratio of the resolution of every layer of LOD is 2, and namely the resolution of sub-block is two times of his father's block, identical geographical space, and the Pixel Dimensions of sub-block is 2 times of upper level, and under identical screen resolution, the geographical space area that father's block covers is 4 times of sub-block;

(III) the block sequence of every layer adopts simple ranks to encode, ratio, every layer of transverse direction of lower block number between one deck and last layer are 2:1 with the ratio of the block number of longitudinal direction, and according to the custom of computer picture coordinate system, starting point is defined in the upper left corner, from top to bottom, column number from left to right for line number; The block count of the 0th layer is 4 × 2, and the block count of kth-1 layer is ;

(IV) contain the dem data in corresponding geographic position in each piecemeal, its resolution dimensions must be odd number, namely ;

(V) pyramidal bottom is highest resolution, directly obtains from raw data, and the LOD level at its place is determined by himself resolution and affiliated scope; The data texturing that each terrain block is corresponding and normal direction data resolution must, higher than the resolution of DEM, to ensure interpolation smoothing, be avoided playing up distortion; Suppose that certain one deck DEM resolution is r, their resolution definition be here ; Obtain the pyramid model of global scale.

The method drawing out three-dimensional terrain model in conjunction with the real-time smoothness of GPU technology according to claim 1, is characterized in that: enhancing technology drawn by the streamline based on GPU, and the data of the pyramid model of above-mentioned global scale are drawn:

Utilize Shader Language the vertex information of elevation and index information separately to be stored, and import video memory into as data texturing is disposable;

The data structure of a, node and tissue:

ChunkedLOD algorithm is adopted to construct block quad-tree structure, namely no longer with single summit or triangle for processing unit, but to comprise the grid model of texture mapping, triangle is packaged into disposable the completing of band and plays up, make full use of the batch processing ability of GPU;

B, buffer memory multiplexing: adopt the method that Geo-Clipmap uses, the apex coordinate in same layer all carry out same conversion;

The data real-time rendering of c, pyramid model: be divided into two classes, one is static node, and another kind is dynamic node;

I, static node data render: namely obtaining the same area different resolution node by sampling in preprocessing process, when playing up, directly playing up by the node of corresponding resolution is called in video memory;

1. ring-type gets the bit arithmetic of location:

Use the grid buffer memory of a fixed size to obtain the altitude figures of arbitrary node in quaternary tree, by setting offset address and geographical space span, play up whole topographic space; Before adding and looking vertebra cutting, get viewpoint vertically downward, the visual range of viewpoint is one n× nthe pixel grid of fixed size, according to the pyramid level at viewpoint place l, obtain the geographical space resolution in this layer representated by each pixel , wherein rrepresent earth radius, use represent that viewpoint is projected in displacement in landform plane, the pixel distance of so viewpoint movement is in the plane expressed as with , the pixel coverage that visual field covers by with represent;

The node shifting out visual field is replaced by the node entering visual field for the first time, and all the other nodes are motionless, this completes incremental update, for any point in field range , its position stored in node array is , wherein modular arithmetic is asked in " mod " expression;

2. the determination of error metrics:

Adopt the error metrics account form of ChunkedLOD, namely from leaf node, along the error of quaternary tree upwards successively computing node, the error metrics of a node for its four child nodes maximal value add itself error metrics ; computing method as follows:

(3.6)

Like this, when screen space error is greater than threshold value, select to continue to find the less node of error downwards, until meet the demands or arrive leaf node;

3. visibility culling:

Elimination method when viewpoint is nearer:

When view distance landform is nearer, adopts bounding box and carry out test for intersection depending on centrum, namely asking friendship successively by the unique point of bounding box and 6 planes depending on centrum, with determine these unique points whether belong to plane forward that half space or and Plane intersects; Depending on 6 planes of vertebra by equation determine, wherein vector represent the normal vector of this plane, numerical value dwrepresent the distance of plane and initial point; 9 unique points of bounding box are made up of 8 angle points and a central point, can be judged the observability of terrain mesh, get by the unique point comparing bounding box sizeequal the half of the diagonally opposing corner line length of bounding box;

Elimination method when viewpoint is far away:

When viewpoint is away from landform, utilize the distance coefficient of node and viewpoint to judge, the actual range of definition distance coefficient size representated by node bounding box, in projection plane length with the distance of bounding box on projection plane ratio, namely ;

Get eye coordinates , bounding box centre coordinate is , then viewpoint to the direction vector on ground is ;

If 8 of bounding box apex coordinates are , cross point with dparallel straight-line equation is:

(3.9)

Depending on the nearly shear surface of centrum as projection plane, with direction vector vertically, if its plane equation is

(3.10)

Above-mentioned two equations of simultaneous can be tried to achieve subpoint coordinate in above-mentioned plane for

(3.11)

Calculate all 8 summits after, get the Min-max in 8 coordinate points , , , , , ; Order

(3.12)

Distance coefficient avalue be:

(3.13)

During judgement, when viewpoint height is greater than setting threshold value, if node bounding box avalue is less than 1, when namely size is less than a pixel on a projection plane, can thinks invisible, no longer play up, otherwise continues to play up this node;

4. gap is eliminated and geometry transition:

Before carrying out node screening, first calculate the geographic range that each LOD level covers, to produce granularity when playing up being similar to consistent number of triangles on window plane;

If the initial distance of every layer of LOD is start _i , irepresent the level at current LOD place, its distance range covered is range, and viewpoint to the distance on current any one summit of LOD is dist, then current transitional region [ mStart, mEnd] and the transfer coefficient on this summit kfor:

(3.14)

Wherein, for clamper function, namely ;

After obtaining transfer coefficient, definition place, summit size of mesh opening is n× n, the quaternary tree degree of depth at place is l, summit normalized coordinate in grid is , then size of mesh opening is relative to the actual range scaling value under world coordinates scaleand the coordinate after transition tried to achieve by following formula:

(3.15)

Wherein, expression is got fraction part;

Finally, apex coordinate can be from transform to position, the elevation of its correspondence is obtained by bilinear interpolation;

II, dynamic node: be when quad-tree structure extends to leaf node, still can not meet the demand of viewpoint threshold value, is calculated in real time by GPU and generates and the node added in quaternary tree;

The generation of dynamic node:

1. the segmentation of relief block: segmentation inserts millet cake and edge point according to segmentation mask in the original vertices of quadrilateral mesh, and connect into new quadrilateral with original vertices, and the non-border vertices number of degrees of generation are 4, and the border vertices number of degrees are 3, are all regular points;

Regular point segmentation process: to quadrilateral v ₅ v ₆ v ₇ v ₈when segmenting, four edge points can be formed on its four edges e ₀ e ₁ e ₂ e ₃, also have a millet cake simultaneously f, make original dough sheet split into four sub-dough sheets; e ₀, e ₁, e ₂, e ₃the following formula of value:

(4.1)

Millet cake fvalue drawn by formula 4.2:

(4.2)

Wherein, ;

2. random fractal: use variable represent the coordinate that curved surface is put, represent this height on curved surface; Then FBM curve surface definition is on certain probability space, and index is a stochastic process

Stochastic variable meet:

(1) set up with probability 1 , namely process is from initial point; And be about continuous function;

(2) to aleatory variable , its two-dimentional increment obeying expectation is 0, and variance is normal distribution, its probability meet

(4.3)；

Adopting two kinds of fractal algorithms to generate the noise pattern of different pattern, is Diamond-Square algorithm and Perlin noise respectively:

Diamond-Square algorithm: segment each time in an iterative process and all can perform two step operations, Diamond step and Square step; Diamond step, for generating foursquare center, is got four summit calculating mean values, is added a random quantity, as the value of central point; Diamond step generates four triangular pyramids; Square step is then four summits of getting pyramid, and add that stochastic variable is as former foursquare four edges mid point with its mean value, the result segmented each time is consistent;

Perlin noise: for an amplitude and noise acoustic image, first normal vector and the gray-scale value of random direction is given to rounded coordinate point, the pixel value of any point is by the normal vector on four summits in this cell that place is the most contiguous in addition, decides with the direction vector being connected this point and four summits

If for the vector of required point in figure, nfor the stacking fold of noise, generally get ; If noise function wavelength is , then frequency , the gray-scale value of this point for:

(4.4)

Wherein, offsetfor skew measures-0.1, it is a noise function;

Pretreated method is taked to generate polytype Noise texture, avoid calculating the extra expenses increasing GPU in real time, in rasterization stage, pixel coloring device reads a texture maps comprising noise and elevation map merges, noise pattern carries out different proportion ground scaling and superposes simultaneously, expand infinitely-great Noise texture, eliminate with this periodicity and medelling effect drawn out.

Of the present invention " proposition of digital earth and " information graticule mesh "; provide new ways and means for overcoming the above problems; namely organization and administration are efficiently carried out to spatial data; and express with digitizing, visual form; and the terrain information studied herein is visual; just by content primary and basic for Data for Virtual Battle Space Visualization: in application, any operation is all carry out in certain battlefield geographical space; Technically, geographical space environment visual is basis and the carrier of the visual Overlapping display of other key elements, has great significance to the three-dimensional expression of research military information.Its advantage is:

1. analyze modeling method and the basic characteristics of landform altitude, present Research and the hot issue in this field are inquired into, establish technology path herein, the problem of pretreatment of analyze and solve altitude figures as the elimination etc. of noise in the conversion of plane projection coordinate system and elevation map, for organization and management altitude figures is laid a good foundation.

2. propose and use pyramid model to build Global Multi-resolution Virtual Terrain Environment, and in conjunction with multiresolution LOD modelling technique, design and Implement the pyramid model of global range, comprise pyramidal structure rule and hierarchical block scheme, adopt quad-tree structure to carry out index to data, utilize bit arithmetic and queue structure to realize back end quick position and renewal.

3. the com-parison and analysis Terrain LOD algorithm of several classics, gets location in conjunction with the quad-tree structure of ChunkedLOD algorithm and the ring-type of GeometryClipmap algorithm, realizes the terrain rendering in the whole world and the incremental update of rending model; Propose the details transition scheme of the method for the visibility culling more accurately and continuous seamless combined with bounding box based on encircle sphere, efficiently avoid the scintillation of drawing distortion and model.

4. acceleration and the enhancing of terrain rendering is achieved based on Programmable GPU technology, namely each stage by using Shader Language to control graphics rendering pipeline, the vertex information of altitude figures and index information are generated respectively two textures stored in video memory, call for whole terrain rendering; Utilize surface subdivision and fractal technology in geometry stage opposite vertexes interpolation and skew, generative process details, the phenomenon of terrain mesh corner angle when making up lack of resolution.

Accompanying drawing explanation

Fig. 1 is utm projection figure;

Fig. 2 is multi-resolution pyramid basic block diagram;

Fig. 3 is Global Multi-resolution landform hierarchical block conceptual scheme;

Fig. 4 is the process flow diagram of filtering noise;

Fig. 5 is the design sketch of filtering noise;

Fig. 6 is the different level of detail figure using different step values to obtain from identical summit buffering;

Fig. 7 is the generative process figure of dynamic node;

Fig. 8 is the schematic diagram of segmentation process;

Fig. 9 is the determination figure of pixel value;

Figure 10 is the composite diagram of noise pattern;

Figure 11 is the noise pattern that two kinds of algorithms generate;

Figure 12 is details synthetic effect figure;

Figure 13 is the elevation map metadata of two groups of data;

Figure 14 is the general view of two scenes;

Figure 15 is the seamless continuous drafting design sketch of LOD that the present invention realizes;

Figure 16 is that the real-time details played up strengthens design sketch;

Figure 17 utilizes GPU monitoring tool NVIDIAnTune and GPU-Z to carry out the runnability of system monitoring the record figure obtained;

Figure 18 is the triangle number of two groups of experiments when roam according to route and the fast variation diagram of corresponding frame;

Figure 19 is that geometric error is converted to screen error figure;

Figure 20 is the what comes into a driver's body figure in perspective projection;

Figure 21 is the result figure after looking vertebra cutting;

Figure 22 is the details allocative decision figure of 6 grades of LOD;

Figure 23 is the geometry transition figure of grid;

Figure 24 is multistage geometry transition figure.

Embodiment

1. have studied pretreatment process terrain data being carried out to organization and administration, from the projected coordinate system that many expression models and elevation adopt, discussing projective transformation successively, noise is eliminated and the implementation method of interpolation sampling, preparing for building multi-resolution Topographic Data pyramid.

2. have studied the scheme of geospace subdivision and the Method of Data Organization of multi-resolution pyramid hierarchical block, carried out organizing storing according to the file layout designed herein to Global Multi-resolution Topographic Data; Use quad-tree structure to each back end index building, and fast algorithm is inserted in the renewal proposing quaternary tree.

3. have studied the Drawing of Three-Dimensional Topography enhancing technology based on GPU.Analyze operation and the division of labor in each stage of graphics rendering pipeline, discuss the rendering algorithm that terrain data combines with GPU programming technique, take full advantage of the computation capability of GPU, realize the multiplexing and real-time details of buffer memory to strengthen, greatly reduce the load of CPU and internal memory, improve rendering performance and display effect.

4., based on research contents herein, give the Outline Design model of global terrain data manage and display system, divided the function of modules, and the workflow of discussing system, with being verified.

The present invention draws out the multi-resolution pyramid model construction of three-dimensional terrain model:

1. the pretreatment process of terrain data: by I, the conversion of projection coordinate, II, the elimination and III of image noise, sample interpolation build multi-resolution pyramid and unified data layout;

I, the conversion of projection coordinate: the anti-solution to utm projection:

Real-time three-dimensional display system can adopt a kind of Descartes's right-handed coordinate system to be unified in the model in world space.But consider and set up georef, it is the rectangular coordinate system in space of initial point that this coordinate system is comparatively close to earth centroid, and all position coordinateses are all in units of rice, uses represent the skew relative to barycenter.But this seat target value is usually very large, has mostly exceeded the expression scope of single knuckle type data, only can preserve the integral part (the highest as positive and negative identifier) of 23 bits in the form of single knuckle type data, change into the decimal system namely 2 ²³-1=8388607, and fraction part only has at most 7.In WGS84 reference system, the equatorial radius of the earth just reaches 6378137m, like this uses single knuckle type form, and most high spatial resolution also can not more than 0.8m(6378137/8388607 ≈ 0.8).Under such threshold value, the distortion that various floating data approximate solution causes will occur, as summit overlap and viewpoint shake etc.

Another one factor is, rectangular coordinate system in space cannot indicate the precision information of drawing intuitively, on the contrary, a lot of geographical or projected coordinate system but can represent more geography information, geographical co-ordinate system is relevant to the reference ellipsoid that it adopts, as longitude and latitude geographic coordinate; Earth model then projects on some simple Developing Surfaces by projection coordinate's system, what LambertConformalConic(LCC) adopt is positive axis conformal conic projection, UniversalTransverseMercator(UTM) what adopt is that transverse axis isogonism cuts conic projection.Thus geospace is projected in method in rectangular coordinate system and inadvisable, this position needing the positive inverse transformation between projected coordinate system and geographic coordinate system to determine a little equally.

What the DEM as terrain source data adopted is WGS84 is as with reference to spheroid, and the coordinate format adopted has two kinds, and one is UTM coordinate, and unit is rice or the km of long, and one is WGS84 latitude and longitude coordinates, and unit is the angle of type double precision.In order to Uniform data format, according to the actual demand of people, when organizing dem data, adopt latitude and longitude coordinates as our linear module.Therefore, need the dem data in units of UTM coordinate to be converted to latitude and longitude coordinates.

Separate UTM coordinate is counter, be namely by , its projection equation:

(2.1)

Wherein, represent the projection value of latitude, represent the projection value of longitude;

If the meridian circle longitude at central meridian place is , each latitude circle corresponding, except equatorial projection is straight line, all the other equal inwardly protrudings are right the intersection point of straight line and central meridian is referred to as end point , its latitude and longitude coordinates is , the latitude and longitude coordinates of this Points on Straight Line with point is that the expression formula of initial point is set to:

(2.2)

And then , former formula just transforms to solve , with the process of three parameters; Because the point on straight line is symmetrical projection relative to end point, and on same straight line value constant, then must be about even function, must be about odd function, according to fourier function launch can obtain:

(2.3)

Wherein, be about function, and , namely ;

By can obtain:

(2.4)

solve and can be drawn by the algorithm of Gauss projection, and then each term coefficient in formula can be obtained successively, finally can obtain:

(2.5)

Wherein being expressed as follows of parameters:

(2.6)；

II, the elimination of image noise:

By analyzing the noise characteristic in elevation map, can find, the noise mainly white noise of elevation, and noise area is comparatively large, in block in image, instead of simple salt noise, noise model meets the feature of additive noise:

(2.12)

Wherein, for being positioned at image the gray-scale value at place, for the gray-scale value not by noise pollution, be independent of random noise value, for being subject to the probability of noise pollution;

Therefore, when filtering, no longer consider that minimal value is as noise spot, and the original dimension of filter window set comparatively greatly, to avoid undetected block distortion.

Before treatment, first clear and definite several definition is once needed:

Suppose image isize is m× n, any point in image gray-scale value be , the average gray value average of entire image ( i _{m× n} ) computing formula be:

(2.13)

The size of moving window is r× rif the coordinate of window center pixel is , gray-scale value is , then pixel average average in this window ( w _{r× r} ) computing formula be:

(2.14)

Uniformity coefficient weighs the parameter of noise signal contrast in image, image i _{m× n} and window w _{r× r} uniformity coefficient be defined as follows respectively:

(2.15)

The interpolation method of Lagrange is:

(2.16)

Wherein, represent the functional value gone out; Herein, interpolation algorithm processes for one-dimensional data, and when being extended on bidimensional image, we get interpolating pixel point level, vertical, two diagonals adjacent pixels point altogether on four directions carries out interpolation, and gets its mean value as last interpolation result;

The experiment effect of algorithm as shown in Figure 5, can see that noise is all detected and successful filtering, although the elevation around noise affects by it improved certain height, relative to landform whole after medium filtering by level and smooth effect the method or acceptable.

Steps flow chart image being carried out to filtering is as follows:

(I), to arrange initial filter window size basicWndSize be 11 × 11, maximum window size maxWndSize is 19 × 19, noise counter noiseCount=0 is set, and characterize the identifier isNoise=false whether pixel is noise, calculate the uniformity coefficient image_HSV of entire image, then the pixel of entire image is traveled through;

Maximum point pixel value iMax in (II), statistics current window, and number nMax, the number of pixels sumPixel in window, and the block uniformity wnd_HSV of window;

(III), judge current window central pixel point gray-scale value whether equal pixel maximum value iMax in window, or whether window uniformity coefficient wnd_HSV is greater than image uniformity coefficient image_HSV; If yes, then continue to judge whether the number nMax of maximum value pixel in current window is less than sumPixel/3, if yes, then mark center pixel is noise spot isNoise=TRUE, noiseCount++; If NO, then think that current extreme value density is excessive, and then have influence on the judgement to noise, now filter window is expanded basicWndSize+=2, jump to step 2; If filter window reaches maximal value, i.e. basicWndSize==maxWndSize, then mark center pixel is noise spot isNoise=TRUE, noiseCount++; Otherwise, medium filtering is carried out to this shop; If NO, then judge that current point is not image noise point, then jump out circulation, moving window is moved to next pixel, until complete traversal to entire image;

(IV), to being judged to be that the four direction of noise spot carries out Lagrange's interpolation;

III, sample interpolation: namely the interpolation of dem data adopts bilinear interpolation method; Get the pixel value of its mean value as noise, the identifier isNoise of this point is set to false, simultaneously noiseCount--.

2. multi-resolution pyramid is built:

Dimensional topography is when carrying out real-time rendering, in order to ensure the continuity of visual effect, need initial landform Grid Sampling or interpolation, form the grid that different resolution is expressed, calling or removing of corresponding resolution grid is controlled according to the error metrics that algorithm presets, but this multi-resolution grid is difficult to obtain by calculating in real time and be applied to playing up of graphic hardware, the less fashion of size of mesh opening is not enough to the decline causing frame speed, but just to seem on the impact of frame speed when viewpoint collecting scope is larger and can not be ignored.

In order to head it off, create the concept of multi-resolution pyramid, pyramid structure is a kind of multiresolution hierarchical model, from the principle of model, pyramid is a kind of resolution model continuously, but be difficult to when building the consecutive variations realizing resolution, the complexity of algorithm and the redundancy in space can improve greatly, therefore, during pyramid construction, the multiplying power method that adopts builds, namely to areal more, adopt a series of grid representation, the sampling interval of neighbouring grid differs identical scale-up factor m(usually get m=2), bottom-up from pyramid like this, resolution successively reduces, and size of mesh opening is corresponding reduction also.If original mesh is of a size of s ₀, resolution is r ₀, constructed pyramid has llayer (pyramid top is the 0th layer), so the kthe size of mesh opening of layer sand resolution rcan be expressed as:

(2.7)

Carry out piecemeal to the earth after plane projection according to waiting longitude and latitude, simultaneously according to top-down mode, build the pyramidal layer of different levels rank, see Fig. 3, the scheme of structure is as follows:

(I) the latitude span of the earth is-90 ° to 90 °, longitude span is-180 ° to 180 °, the ratio of wide height is on a projection plane just 2:1, therefore quaternary tree is divided into 8 root nodes the 0th grade time, laterally equal 90 ° of longitudinal span, with the satisfied requirement waiting longitude and latitude to divide, it is invalid that the coordinate figure except this scope is all considered as;

(II) ratio of the resolution of every layer of LOD is 2, and namely the resolution of sub-block is two times of his father's block, identical geographical space, and the Pixel Dimensions of sub-block is 2 times of upper level, and under identical screen resolution, the geographical space area that father's block covers is 4 times of sub-block;

(III) the block sequence of every layer adopts simple ranks to encode, ratio, every layer of transverse direction of lower block number between one deck and last layer are 2:1 with the ratio of the block number of longitudinal direction, and according to the custom of computer picture coordinate system, starting point is defined in the upper left corner, from top to bottom, column number from left to right for line number; The block count of the 0th layer is 4 × 2, and the block count of kth-1 layer is ;

(IV) contain the dem data in corresponding geographic position in each piecemeal, its resolution dimensions must be odd number, namely ; To ensure that block adjacent with it has the Duplication of a pixel, be so not only convenient to the appearance avoiding T-shaped crack, during to altitude figures interpolation, adjacent terrain block junction height can be made consistent, eliminate the generation of terrace with edge phenomenon;

(V) pyramidal bottom is highest resolution, directly obtains from raw data, and the LOD level at its place is determined by himself resolution and affiliated scope; The data texturing that each terrain block is corresponding and normal direction data resolution must, higher than the resolution of DEM, to ensure interpolation smoothing, be avoided playing up distortion; Suppose that certain one deck DEM resolution is r, their resolution definition be here ; Obtain the pyramid model of global scale.

Its property parameters is as shown in the table:

The hierarchical block result of table pyramid model

Wherein, R represents the earth major radius that WGS84 defines: 6378137.0m.It is the span difference of longitude and latitude to the foundation of global metadata piecemeal, and the unit of the 4th, 5 row represents with distance in table, to meet the use habit of people better, their value is taken from the longest latitude distance of place block leap and the ratio with pixel resolution thereof.

Enhancing technology drawn by the streamline that the present invention is based on GPU: dynamic image pro cess:

Utilize Shader Language the vertex information of elevation and index information separately to be stored, and import video memory into as data texturing is disposable;

The data structure of a, node and tissue:

The world is divided into many zonules, eachly represents with rectangular parallelepiped, be referred to as block (sector), these blocks formed a hierarchical tree, each block is made up of four less sub-block.Carry out ground scape when playing up, first judge which block is visible, if father's block is invisible, so child is also invisible; Otherwise then child nodes is all visible.And when father's block is crossing with looking centrum, just select its sub-block, it is continued be used as father's block to test its observability.Until all visible blocks are found or complete the traversal of all nodes.

ChunkedLOD algorithm is adopted to construct block quad-tree structure, namely no longer with single summit or triangle for processing unit, but to comprise the grid model of texture mapping, triangle is packaged into disposable the completing of band and plays up, make full use of the batch processing ability of GPU;

The elevation map getting a 2D is example, first uses Geomipmap to carry out filtering to it, builds the mipmap pyramid that contains L layer.Here, the data of every layer are all saved as 2D texture with suitable size by us, instead of the buffering by vertical array, it linearly being turned to 1D thought.Pyramidal every one deck contains the sampled point of n × n, and we store a single pass elevation map for it, a three-channel vertex normal figure and three-channel texture maps.

Like this, our node data structures that can be defined as follows in tinter:

structVSOUTPUT{

vectorposition:POSITION;

floatheight:TEXCOORD0;

float3normal:TEXCOORD1;

float2uv:TEXCOORD2;

}

The resolution of normal map and texture is allowed to be the twice of elevation, because each summit only has the effect of a normal very fuzzy, simultaneously also can be very coarse after adding illumination, and the normal map of two resolutions can create softer model by interpolation.The scope of the texture restriction viewpoint of low resolution, each surface type is produced large-area fuzzy, therefore fixture resolution also should not be less than the resolution of height map.

B, buffer memory multiplexing: adopt the method that Geo-Clipmap uses, the apex coordinate in same layer all carry out same conversion;

Set up two classes, one is cTerrain, and it represents whole landform, help the part landform cutting out actual needs, and this landform is put under block corresponding to grid.These blocks use second class to represent, cTerrainSection stores elevation and the normal of each block in this type of.Such as, the landform on 256 × 256 summits, establishment cTerrain is stored whole data set by us.Landform is subdivided into the summit block of 32 × 32 by cTerrain object further, creates 64 cTerrainSection objects altogether and represents each block.Since all cTerrainSection objects are all equidimensions, etc. quantity summit, we just can set up an independent indexed cache object and use it to play up each block of landform, for different levels, only corresponding step value need be moved right one, as shown in Figure 6.

The code of index building buffering is as follows:

boolcTerrain::buildIndexBuffer(

IntgridVerts, the number of vertices on // Grid Edge

Intstep) step value that // index buffering uses

{ inttotalStrips=gridVerts-1; The number of // band

Intindice_tripwise=gridVerts<< 1; // each band appears on two lines summit

// due to the existence of 0 triangle, number of vertices increases

inttotal_indice=(indice_stripwise*totalStirps)+(totalStrips<<1)-2;

int*index=newint[total_indice];

intvert,start_vert=0;

for(intj=0;j<total_strips;++j)

The index buffering of { // generate a line

vert=start_vert;

for(intk=0;k<xVerts;++k)

{*(index++)=vert;

*(index++)=vert+gridVerts*step;

vert+=xStep;

}

start_vert+=lineStep;

if(j+1<total_strips)

{ // increase by 0 triangle at the end of every a line

*(index++)=(vert-xStep)+lineStep;

*(index++)=start_vert;

}}}

Build the concrete way of summit buffering, by position and the texture coordinate of each piece be contracted in scope in.When playing up, coordinate is amplified to the ratio of needs by us, and in addition namely suitable off-set value is movable to corresponding world locations.And each terrain block uses a series of similar xy position and uv data texturing, we quote them again by when playing up each cTerrainSection object in their disposable storages and cTerrain object.We construct two groups of vertex data streams like this: be stored in the shared vertex data of cTerrain and the independent vertices data of each cTerrainSection, greatly reduce the memory headroom that it takies.

Build summit buffering code as follows:

boolcTerrain::create(

cTexture*heightMap,

constcRect3d&worldExtents,

intshift)

{

// scaling ratio is set

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;

// initial offset is set

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=

newsLocalVertex[m_sectorVerts*m_sectorVerts];

// summit of sharing and texture coordinate write and flow

for(inty=0;y<m_sectorVerts;++y)

{

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

for(intx=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;

}}}

The real-time rendering of c, three-dimensional terrain model: be divided into two classes, one is static node, and another is dynamic node;

I, large-scale terrain data render: namely obtaining the same area different resolution node by sampling in preprocessing process, when playing up, directly playing up by the node of corresponding resolution is called in video memory;

Adopt the quad-tree structure of ChunkedLOD to organize terrain data, draw in conjunction with the nested grid pattern in Geo-Clipmap simultaneously, and use the transition of a kind of continuous print geometry to adjust the fission process of grid when LOD changes, avoid the generation in sudden change and crack.

1. ring-type gets the bit arithmetic of location:

Use the grid buffer memory of a fixed size to obtain the altitude figures of arbitrary node in quaternary tree, by setting offset address and geographical space span, play up whole topographic space; Before adding and looking vertebra cutting, get viewpoint vertically downward, the visual range of viewpoint is one n× nthe pixel grid of fixed size, according to the pyramid level at viewpoint place l, obtain the geographical space resolution in this layer representated by each pixel , wherein rrepresent earth radius, use represent that viewpoint is projected in displacement in landform plane, the pixel distance of so viewpoint movement is in the plane expressed as with , the pixel coverage that visual field covers by with represent;

Pixel coverage after viewpoint moves becomes with , when upgrading node queue, from node to queue, map below demand fulfillment three conditions:

(1) same element all the time the same element in corresponding queue , namely from iarrive meet the mapping relations of 1:1;

(2) to arbitrary element , necessarily exist correspond, ensure the continuity mapped;

(3) if , so , ensure the nonredundancy mapped.

First condition ensure that the possibility of incremental update, has joined in queue once node, if node is still retained in visual field after viewpoint moves, so does not just need it to shift out queue; Second and the 3rd condition prevent data from being rewritten, make queue minimally space.Realize ring-type and get location, utilize modular arithmetic to realize above-mentioned mapping.

The node shifting out visual field is replaced by the node entering visual field for the first time, and all the other nodes are motionless, this completes incremental update, for any point in field range , its position stored in node array is , wherein modular arithmetic is asked in " mod " expression;

Code is as follows:

intCoordWrap(intval,intlimit)

{

If (val<0) returnlimit+val; When // displacement is negative direction, shift out behind border around value

If (val>limit-1) returnval-limit; When // displacement is positive dirction, shift out behind border around value

Returnval; If // node is also in scope, its value is constant

}

for(intx=0;x<=n-1;x++)

{

for(inty=0;y<=n-1;y++)

{

// obtain the coordinate position of node in queue

x0=CoordWrap(x+xdiff/blockDim,n);

y0=coordWrap(y+ydiff/blockDim,n);

for(inti=0;i<=blockDim-1;i++)

{

for(intj=0;j<=blockDim-1;j++)

The pixel value on all summits in // acquisition node

GetPixel(i+xdiff+xmin,i+ydiff+ymin,x0,y0&pixel);

}

}

}

Wherein, represent the actual coordinate of node, represent the coordinate of node in queue.So, be no matter viewpoint to any direction motion time, all can be read the coordinate of correspondence position by a function, eliminate the tedious steps repeatedly judging the positive negative direction of displacement in the circulating cycle.

When viewpoint position to move the width of a node to upper right side, node ABCDEFG enters activeregion, and " L " shape node in the lower left corner correspondingly can be removed, and in nodal cache, base node and side node are substituted respectively.When continuing to move the width of a node to upper right, node HIJKLMN continues to enter activeregion, except the second row and the original node of secondary series are removed, also note that node A and G has also been removed simultaneously, queue array is now different according to shading, four regions are divided into, except original node in the upper right corner remains unchanged, EFMN is a region, JKDL is a region, HIBC also constitutes a region, the mode that ring-type gets location makes " L " shape update area in queue, become "+" shape, thus read more new node time can also read according to these three Region dividing.

2. the determination of error metrics:

Adopt the error metrics account form of ChunkedLOD, namely from leaf node, along the error of quaternary tree upwards successively computing node, ensure that the error metrics of father node is greater than child node, all the time for the error metrics of each node , value in , wherein for being removed the sum on summit in the grid that node comprises, represent in simplification process the height difference that after removing a summit, landform produces.And in order to ensure the nested of error, the error metrics of a node for its four child nodes maximal value add itself error metrics ; computing method as follows:

(3.6)

Error be the geometric error based on three-dimensional article space measurement, when utilizing perspective projection geometric error to be projected in the screen space of two dimension, real-time LOD selects the screen space error threshold that can utilize two dimension better.After determining viewpoint parameter, according to current geometric error determine corresponding screen error , as shown in figure 19, represent the width of viewing volume, represent the distance of viewpoint to current layer, for field angle size, for screen differentiates (in units of pixel).Can be drawn by similar triangles:

(3.7)；

Like this, when screen space error is greater than threshold value, select to continue to find the less node of error downwards, until meet the demands or arrive leaf node;

3. visibility culling:

Choosing of terrain mesh is the process of an iteration, can constantly repeat to calculate distance and error to each node grid, so just relates to a large amount of caps.In order to reduce pending number of grid, the judgment mode adopting bounding box to combine with encircle sphere improves arithmetic speed.The enclosure body of an object this object can be wrapped in interior a kind of shape, by getting rid of rapidly and looking the disjoint enclosure body of centrum, can reduce much unnecessary computing.If not crossing with enclosure body depending on vertebra, so it must be crossing with the object of enclosure body inside, and the simple enclosure body of selected shape is than judging that the object that its inner profile is complicated will save a lot of time.

Elimination method when viewpoint is nearer:

When view distance landform is nearer, adopts bounding box and carry out test for intersection depending on centrum, namely asking friendship successively by the unique point of bounding box and 6 planes depending on centrum, with determine these unique points whether belong to plane forward that half space or and Plane intersects; This just means and to be made up of the positive semispace of these planes depending on centrum space.Therefore, if rectangle is present in the negative semispace of any cone plane, so this object is just outside the cone.Depending on 6 planes of vertebra by equation determine, wherein vector represent the normal vector of this plane, numerical value dwrepresent the distance of plane and initial point; 9 unique points of bounding box are made up of 8 angle points and a central point, can be judged the observability of terrain mesh, get by the unique point comparing bounding box sizeequal the half of the diagonally opposing corner line length of bounding box;

For any point P in space with plane relation, computing formula:

(3.8)

Namely the absolute value of result represents a little to the vertical range of plane, compare absolute value with sizesize, just can obtain node and whether look vertebra inside.

Elimination method when viewpoint is far away:

When viewpoint is away from landform, utilize the distance coefficient of node and viewpoint to judge, the actual range of definition distance coefficient size representated by node bounding box, in projection plane length with the distance of bounding box on projection plane ratio, namely ;

Get eye coordinates , bounding box centre coordinate is , then viewpoint to the direction vector on ground is ;

If 8 of bounding box apex coordinates are , cross point with dparallel straight-line equation is:

(3.9)

Depending on the nearly shear surface of centrum as projection plane, with direction vector vertically, if its plane equation is

(3.10)

Above-mentioned two equations of simultaneous can be tried to achieve subpoint coordinate in above-mentioned plane for

(3.11)

Calculate all 8 summits after, get the Min-max in 8 coordinate points , , , , , ; Order

(3.12)

Distance coefficient avalue be:

(3.13)

During judgement, when viewpoint height is greater than setting threshold value, if node bounding box avalue is less than 1, when namely size is less than a pixel on a projection plane, can thinks invisible, no longer play up, otherwise continues to play up this node;

4. gap is eliminated and geometry transition:

Propose the control technology of a kind of new LOD level conversion, leave out the degeneration triangle update between different resolution grid simultaneously, avoid unnecessary play up burden and cause thus play up distortion, to realize the LOD transition of seamless process.

In order to control to choose suitable level of detail under certain distance, before carrying out node screening, first calculate the geographic range that each LOD level covers, to produce granularity when playing up being similar to consistent number of triangles on window plane; Due to the definition of quad-tree structure, every layer of details is fixed as 4 with the ratio of the vertex number of its upper level, is also thus just 2 at one-dimensional square distance multiple upwards, and therefore the transfer coefficient of different resolution is set in about 2 to ensure that triangle number is approximate consistent.Figure 22 is the distance range that different LOD levels covers, and along with resolution reduces, the geographic range shared by every layer of terrain block is 2 times of upper level, and last one deck (apart from viewpoint farthest, level of detail is minimum) namely represents the visual range of whole viewpoint.

In order to ensure can to seamlessly transit between adjacent LOD level, retain landform precision as much as possible simultaneously, we adopt the transition scheme consistent with Geo-Clipmap algorithm, namely in the region that low resolution grid is adjacent with fine-resolution meshes, a part is marked off as transitional region, transition from low to high, chooses 33% of each LOD scope as transitional region in experiment.The spinoff done like this to produce more triangle, but be different from leg-of-mutton generation method of degenerating, and the triangle generated by top displacement method wants simple many in complexity.

The summit of all transitional regions is all convert according to self unique error metrics, and namely transition is carried out in pointwise, and two kinds of mapping modes are all supported on each summit: by itself to more high-resolution transition with by transition to more low resolution itself.Getting the second mapping mode is example, namely blending operation is the mesh of vertices of each 3 × 3, is made up of 8 triangles, by mobile 4 limit mid points step by step and 1 central point position to lower-left angle point, make it the mesh of vertices being reduced to 2 × 2, be made up of 2 triangles.Whole process seems, two triangles as the upper right corner expand gradually, and remaining 6 triangles little by little to left down direction be collapsed into degeneration triangle and final no longer grid turns to the process of pixel.This method can eliminate the gap between adjacent LOD level effectively, and produces level and smooth transition, there will not be T-shaped connection.

First, calculate displacement coefficient is decided on summit size to the approximate distance of viewpoint, the vertex position being used for calculating distance obtains according to the position proportional of grid in world coordinate system is approximate, because grid is arranged continuously in coordinate, the vertex position then obtained also should be fixed, the edge vertices of adjacent mesh must be consistent, so just can ensure can not produce gap between grid.In an experiment, with take from latitude and longitude value respectively, height value is then drawn by elevation map continuous sampling.By Figure 22, if the initial distance of every layer of LOD is start _i , irepresent the level at current LOD place, its distance range covered is range, and viewpoint to the distance on current any one summit of LOD is dist, then current transitional region [ mStart, mEnd] and the transfer coefficient on this summit kfor:

(3.14)

Wherein, for clamper function, namely ;

So, ksize be namely limited at 0(and do not converted completely, rest on low resolution position) and 1(convert completely, triangle number becomes original twice) between, and along with kvalue become large gradually from 0 to 1, grid also achieves by low resolution to high-resolution transition.After acquisition transfer coefficient, definition place, summit size of mesh opening is n× n, the quaternary tree degree of depth at place is l, summit normalized coordinate in grid is , then size of mesh opening is relative to the actual range scaling value under world coordinates scaleand the coordinate after transition tried to achieve by following formula:

(3.15)

Wherein, expression is got fraction part;

Finally, apex coordinate can be from transform to position, the elevation of its correspondence is obtained by bilinear interpolation; After the transition of all summits completes, current grid will overlap completely with the fine-resolution meshes of its lower one deck, and the process of period can not produce any gap.Figure 23 is the mesh shape that different transfer coefficient is corresponding.

In fact, because viewpoint is continually varying, and in whole visual field, selected node resolution reduces gradually along with the intensification of what comes into a driver's, thus all can there is the geometry transition of more than one deck in each frame, as shown in figure 24, multilayer transitional zone can closely connect adjacent resolution grid, avoids the generation in crack.Experiment shows, when carrying out terrain walkthrough, " the kick sense " that cause when between grid, level and smooth continuous print transition method effectively eliminates node segmentation, frame speed is also stabilized in more than 60 frames, can not produce visual stagnation.

II, dynamic node, it is when quad-tree structure extends to leaf node, still can not meet the demand of viewpoint threshold value, is calculated in real time generate and the node added in quaternary tree by GPU.

The generation of dynamic node:

The segmentation of relief block: segmentation inserts millet cake and edge point according to segmentation mask in the original vertices of quadrilateral mesh, and connect into new quadrilateral with original vertices, and the non-border vertices number of degrees of generation are 4, and the border vertices number of degrees are 3, are all regular points;

As shown in Figure 8, regular point segmentation process: to quadrilateral v ₅ v ₆ v ₇ v ₈when segmenting, four edge points can be formed on its four edges e ₀ e ₁ e ₂ e ₃, also have a millet cake simultaneously f, make original dough sheet split into four sub-dough sheets; e ₀, e ₁, e ₂, e ₃the following formula of value:

(4.1)

Millet cake fvalue drawn by formula 4.2:

(4.2)

Wherein, ;

This algorithm of subdivision can only be applied on the regular grid of regular point place, does not consider for singular point.Due to the feature of ChunkedLOD algorithm, thinning process independently carries out in each piece, grid in adjacent block is had no way of finding out about it, when mask moves to the edge of block, extrapolate and replace angle point and all can produce false edge in traditional limit, in order to ensure the visual sense of reality, we store the limit sampled point of adjacent block again in block structure, only just use it as reference mark when opposite side carries out refinement, segmentation in block is not then considered it, exchanges the continuity at edge for increase low volume data for cost.

Random fractal: use variable represent the coordinate that curved surface is put, represent this height on curved surface; Then FBM curve surface definition is on certain probability space, and index is a stochastic process

Stochastic variable meet:

(1) set up with probability 1 , namely process is from initial point; And be about continuous function;

(2) to aleatory variable , its two-dimentional increment obeying expectation is 0, and variance is normal distribution, its probability meet

(4.3)；

The dimensional topography using fractal technology to draw can produce the scene true to nature infinitely with self-similarity, is easy to smoothing and interpolation.Here adopting two kinds of fractal algorithms to generate the noise pattern of different pattern, is Diamond-Square algorithm and Perlin noise respectively.

1.Diamond-Square algorithm

The essence of algorithm is the Random Mid-point Displacement algorithm on two dimensional surface.The square that it is made up of Seed Points starts, and after several times Diamond-Square algorithm iteration, constantly segmentation seed square is to obtain dimensional topography true to nature.Segment each time in an iterative process and all can perform two step operations, Diamond step and Square step.

Diamond step, for generating foursquare center, is got four summit calculating mean values, is added a random quantity, as the value of central point.Diamond step generates four triangular pyramids.Square step is then four summits of getting pyramid, and add that stochastic variable (obeying identical probability distribution with the stochastic variable that Diamond walks) is as former foursquare four edges mid point with its mean value, the result segmented each time is consistent with Fig. 8.

In each stage of process, random value is always created in fixed range, and we claim present scope to be delta, total size be [- delta, delta].When upper segmentation once, deltaa scaling value can be multiplied by, be referred to as roughness (roughness), between 0 to 1, be used for reducing deltaspan.High roughness value can produce more unordered landform, and low roughness can make landform more and more smooth.

Diamond-Square algorithm realization is simple, and arithmetic speed is very fast, is widely adopted.But owing to there is no information transmission between adjacent square, thus cause producing " creasing problem " in noise pattern, namely there is factitious linear track in image, and do not remove by local smoothing method; In addition, its grey iterative generation mode calculated amount is comparatively large, expends too much time and underaction.

2.Perlin noise

Perlin noise is a function generating random noise, Perlin noise is by obtaining net result to the synthesis of multifractal, namely the successively increase of finite bandwidth frequency is utilized to synthesize, and the calculating of each frequency range is separate in building-up process, therefore it has stronger control ability to the frequency of Fractal Terrain, dimension and other corresponding statistical property.

For an amplitude and noise acoustic image, first normal vector and the gray-scale value of random direction is given to rounded coordinate point, the pixel value of any point is by the normal vector on four summits in this cell that place is the most contiguous in addition, decides with the direction vector being connected this point and four summits.

For an amplitude and noise acoustic image, first normal vector and the gray-scale value of random direction is given to rounded coordinate point, the pixel value of any point is by the normal vector on four summits in this cell that place is the most contiguous in addition, decides with the direction vector being connected this point and four summits.

If for the vector of required point in figure, nfor the stacking fold of noise, generally get .If noise function wavelength is , then frequency , the gray-scale value of this point for:

(4.4)

Wherein, offsetfor skew measures-0.1, it is a noise function.Bottom panel show two independent noise add and after appearance.

Take pretreated method to generate polytype Noise texture, avoid calculating in real time the extra expenses increasing GPU, in rasterization stage, pixel coloring device reads a texture maps comprising noise and elevation map merges, and such implementation effect is quicker.Noise pattern carries out different proportion ground scaling and superposes simultaneously, just can expand infinitely-great Noise texture, eliminate with this periodicity and medelling effect drawn out.That Figure 11 shows is the Noise texture figure that two kinds of Noise Algorithms generate, Figure 12 is then the effect produced after being superposed with elevation respectively by two kinds of noises.As seen from the figure, the two can strengthen proceduring details in real time, simulates different terrain effects, and Diamond-Square noise is more close to sand ground effect, and Perlin noise is then close to rock effect.

Realize effect:

Two groups of experimental datas, construct the dimensional topography of world's elevation of 8 grades of LOD and the PugetSound of 12 grades of LOD respectively, wherein the resolution of 8 grades of world's elevations is 152m, cover the entirely spherical space within the whole north and south latitude of geographic range 83 °, the original mesh built is by 268,435,456 triangle compositions, and the elevation resolution of the PugetSound of 12 grades is 10m, the geographic range covered is 163km × 163km, the original mesh built is by 463, and 601,664 triangles form.Because the data layout of elevation information does not adopt compression, the LOD file size of 12 grades reaches 2.58GB, and the normal direction figure of its correspondence is then close to there being 4GB.

Following table is open the change that details strengthens front and back system performance.Data show, and along with the reduction of viewpoint, also gradually deeply, triangle number constantly increases LOD level, but the block crossing with viewing volume decreases, so the triangle number of the 8th grade of LOD reduces suddenly when highly lower simultaneously.9th grade and the 10th grade is the result after dynamic node adds, and triangle number starts to increase, but loading speed also slows down to some extent.Equally, first 8 grades is the static node generated by preprocessing process, and only need read corresponding terrain block from store when browsing, frame speed is comparatively stable, and the insertion of dynamic node brings certain burden to real-time rendering, terrain complexity and triangle number are depended in the change of frame speed.

Table static node compared with the working time of dynamic node

As can be seen from Figure 17 two groups of Monitoring Data, GPU assume responsibility for most computing load, the occupancy of CPU and internal memory is in lower level always, the quaternary tree degree of depth used due to two groups of data is different, the node number held in field range is also not quite similar, this causes video memory to take inconsistent main cause, but compare the size of source document GB level, and native system has achieved the function utilizing a small amount of buffer memory to browse large-scale terrain.

The data of statistics gained, the system triangle number that operationally each frame average is drawn reaches 265.6k, and frame speed is then close to 194s ^-1it can thus be appreciated that utilize GPU obvious to the efficiency enhancement effect drawn herein, and found out by Figure 18, frame speed is not strict negative correlativing relation with triangle number, along with the depth down of quaternary tree, the details of landform is more and more abundanter, can cause the obvious decline of frame speed, namely frame speed affects also comparatively large by the complexity of landform, and world's elevation of 8 grades mostly is smooth landform, therefore its frame speed is less by triangle number object image, and the PugetSound of 12 grades then fluctuates larger.Following table lists the frame speed ratio of some algorithms comparatively, but the hardware platform adopted due to different algorithms is different, and algorithm is herein comparatively strong to the dependence of GPU, is therefore difficult to carry out lateral comparison, the rendering hardware load that algorithm causes should be reduced as far as possible in next step research.

The drafting efficiency comparison of table algorithms of different

The present invention carries out the design and implimentation of landform multi-resolution display verification system.First the design proposal of system is described respectively from function two aspects of framed structure and modules, then the major function of system is demonstrated by two groups of experimental datas, and shown with chart, the inspection software eventually passing specialty has carried out performance test, proves that this system meets in real time and the rendering request of GPU close friend.

Claims (2)

1. draw out a method for three-dimensional terrain model in conjunction with the real-time smoothness of GPU technology, it is characterized in that:

The multi-resolution pyramid model construction of three-dimensional terrain model:

1. the pretreatment process of terrain data: filtering is carried out to raw video and obtains muting image, and build multi-resolution pyramid;

I, the elimination of image noise:

Noise model meets the feature of additive noise:

(2.12)

Wherein, for being positioned at image the gray-scale value at place, for the gray-scale value not by noise pollution, be independent of random noise value, for being subject to the probability of noise pollution;

Steps flow chart image being carried out to filtering is as follows:

(I), calculate the uniformity coefficient image_HSV of entire image, then the pixel of entire image is traveled through;

Uniformity coefficient weighs the parameter of noise signal contrast in image, image i _{m× n} and window w _{r× r} uniformity coefficient be defined as follows respectively:

(2.15)

Suppose image isize is m× n, any point in image gray-scale value be , the average gray value average of entire image ( i _{m× n} ) computing formula be:

(2.13)

The size of moving window is r× rif the coordinate of window center pixel is , gray-scale value is , then pixel average average in this window ( w _{r× r} ) computing formula be:

(2.14)

Maximum point pixel value iMax in (II), statistics current window, and number nMax, the number of pixels sumPixel in window, and the block uniformity wnd_HSV of window;

(III), judge current window central pixel point gray-scale value whether equal pixel maximum value iMax in window, or whether window uniformity coefficient wnd_HSV is greater than image uniformity coefficient image_HSV;

(IV), to being judged to be that the four direction of noise spot carries out Lagrange's interpolation; The interpolation method of Lagrange is:

(2.16)

Wherein, represent the functional value gone out; Herein, interpolation algorithm processes for one-dimensional data, when being extended on bidimensional image, gets interpolating pixel point level, vertical, two diagonals adjacent pixels point altogether on four directions carries out interpolation, and gets its mean value as last interpolation result;

2. by the data construct multi-resolution pyramid after above-mentioned removing noise:

Pyramid structure is a kind of multiresolution hierarchical model, and adopt multiplying power method to build, namely to areal, adopt a series of grid representation, the sampling interval of neighbouring grid differs identical scale-up factor m, bottom-up from pyramid like this, resolution successively reduces, and size of mesh opening is corresponding reduction also;

If original mesh is of a size of s ₀, resolution is r ₀, constructed pyramid has llayer, pyramid top is the 0th layer, so the kthe size of mesh opening of layer sand resolution rcan be expressed as:

(2.7)

Carry out piecemeal to the earth after plane projection according to waiting longitude and latitude, simultaneously according to top-down mode, build the pyramidal layer of different levels rank, the scheme of structure is as follows:

(I) the latitude span of the earth is-90 ° to 90 °, longitude span is-180 ° to 180 °, the ratio of wide height is on a projection plane just 2:1, therefore quaternary tree is divided into 8 root nodes the 0th grade time, laterally equal 90 ° of longitudinal span, with the satisfied requirement waiting longitude and latitude to divide, it is invalid that the coordinate figure except this scope is all considered as;

(II) ratio of the resolution of every layer of LOD is 2, and namely the resolution of sub-block is two times of his father's block, identical geographical space, and the Pixel Dimensions of sub-block is 2 times of upper level, and under identical screen resolution, the geographical space area that father's block covers is 4 times of sub-block;

(III) the block sequence of every layer adopts simple ranks to encode, ratio, every layer of transverse direction of lower block number between one deck and last layer are 2:1 with the ratio of the block number of longitudinal direction, and according to the custom of computer picture coordinate system, starting point is defined in the upper left corner, from top to bottom, column number from left to right for line number; The block count of the 0th layer is 4 × 2, and the block count of kth-1 layer is ;

(IV) contain the dem data in corresponding geographic position in each piecemeal, its resolution dimensions must be odd number, namely ;

(V) pyramidal bottom is highest resolution, directly obtains from raw data, and the LOD level at its place is determined by himself resolution and affiliated scope; The data texturing that each terrain block is corresponding and normal direction data resolution must, higher than the resolution of DEM, to ensure interpolation smoothing, be avoided playing up distortion; Suppose that certain one deck DEM resolution is r, their resolution definition be here ; Obtain the pyramid model of global scale.

2. the method drawing out three-dimensional terrain model in conjunction with the real-time smoothness of GPU technology according to claim 1, is characterized in that: enhancing technology drawn by the streamline based on GPU, and the data of the pyramid model of above-mentioned global scale are drawn:

Utilize Shader Language the vertex information of elevation and index information separately to be stored, and import video memory into as data texturing is disposable;

The data structure of a, node and tissue:

ChunkedLOD algorithm is adopted to construct block quad-tree structure, namely no longer with single summit or triangle for processing unit, but to comprise the grid model of texture mapping, triangle is packaged into disposable the completing of band and plays up, make full use of the batch processing ability of GPU;

B, buffer memory multiplexing: adopt the method that Geo-Clipmap uses, the apex coordinate in same layer all carry out same conversion;

The data real-time rendering of c, pyramid model: be divided into two classes, one is static node, and another kind is dynamic node;

I, static node data render: namely obtaining the same area different resolution node by sampling in preprocessing process, when playing up, directly playing up by the node of corresponding resolution is called in video memory;

1. ring-type gets the bit arithmetic of location:

Use the grid buffer memory of a fixed size to obtain the altitude figures of arbitrary node in quaternary tree, by setting offset address and geographical space span, play up whole topographic space; Before adding and looking vertebra cutting, get viewpoint vertically downward, the visual range of viewpoint is one n× nthe pixel grid of fixed size, according to the pyramid level at viewpoint place l, obtain the geographical space resolution in this layer representated by each pixel , wherein rrepresent earth radius, use represent that viewpoint is projected in displacement in landform plane, the pixel distance of so viewpoint movement is in the plane expressed as with , the pixel coverage that visual field covers by with represent;

The node shifting out visual field is replaced by the node entering visual field for the first time, and all the other nodes are motionless, this completes incremental update, for any point in field range , its position stored in node array is , wherein modular arithmetic is asked in " mod " expression;

2. the determination of error metrics:

Adopt the error metrics account form of ChunkedLOD, namely from leaf node, along the error of quaternary tree upwards successively computing node, the error metrics of a node for its four child nodes maximal value add itself error metrics ; computing method as follows:

(3.6)

Like this, when screen space error is greater than threshold value, select to continue to find the less node of error downwards, until meet the demands or arrive leaf node;

3. visibility culling:

Elimination method when viewpoint is nearer:

When view distance landform is nearer, adopts bounding box and carry out test for intersection depending on centrum, namely asking friendship successively by the unique point of bounding box and 6 planes depending on centrum, with determine these unique points whether belong to plane forward that half space or and Plane intersects; Depending on 6 planes of vertebra by equation determine, wherein vector represent the normal vector of this plane, numerical value dwrepresent the distance of plane and initial point; 9 unique points of bounding box are made up of 8 angle points and a central point, can be judged the observability of terrain mesh, get by the unique point comparing bounding box sizeequal the half of the diagonally opposing corner line length of bounding box;

Elimination method when viewpoint is far away:

When viewpoint is away from landform, utilize the distance coefficient of node and viewpoint to judge, the actual range of definition distance coefficient size representated by node bounding box, in projection plane length with the distance of bounding box on projection plane ratio, namely ;

Get eye coordinates , bounding box centre coordinate is , then viewpoint to the direction vector on ground is ;

If 8 of bounding box apex coordinates are , cross point with dparallel straight-line equation is:

(3.9)

Depending on the nearly shear surface of centrum as projection plane, with direction vector vertically, if its plane equation is

(3.10)

Above-mentioned two equations of simultaneous can be tried to achieve subpoint coordinate in above-mentioned plane for

(3.11)

Calculate all 8 summits after, get the Min-max in 8 coordinate points , , , , , ; Order

(3.12)

Distance coefficient avalue be:

(3.13)

During judgement, when viewpoint height is greater than setting threshold value, if node bounding box avalue is less than 1, when namely size is less than a pixel on a projection plane, can thinks invisible, no longer play up, otherwise continues to play up this node;

4. gap is eliminated and geometry transition:

Before carrying out node screening, first calculate the geographic range that each LOD level covers, to produce granularity when playing up being similar to consistent number of triangles on window plane;

If the initial distance of every layer of LOD is start _i , irepresent the level at current LOD place, its distance range covered is range, and viewpoint to the distance on current any one summit of LOD is dist, then current transitional region [ mStart, mEnd] and the transfer coefficient on this summit kfor:

(3.14)

Wherein, for clamper function, namely ;

After obtaining transfer coefficient, definition place, summit size of mesh opening is n× n, the quaternary tree degree of depth at place is l, summit normalized coordinate in grid is , then size of mesh opening is relative to the actual range scaling value under world coordinates scaleand the coordinate after transition tried to achieve by following formula:

(3.15)

Wherein, expression is got fraction part;

Finally, apex coordinate can be from transform to position, the elevation of its correspondence is obtained by bilinear interpolation;

II, dynamic node: be when quad-tree structure extends to leaf node, still can not meet the demand of viewpoint threshold value, is calculated in real time by GPU and generates and the node added in quaternary tree;

The generation of dynamic node:

1. the segmentation of relief block: segmentation inserts millet cake and edge point according to segmentation mask in the original vertices of quadrilateral mesh, and connect into new quadrilateral with original vertices, and the non-border vertices number of degrees of generation are 4, and the border vertices number of degrees are 3, are all regular points;

Regular point segmentation process: to quadrilateral v ₅ v ₆ v ₇ v ₈when segmenting, four edge points can be formed on its four edges e ₀ e ₁ e ₂ e ₃, also have a millet cake simultaneously f, make original dough sheet split into four sub-dough sheets; e ₀, e ₁, e ₂, e ₃the following formula of value:

(4.1)

Millet cake fvalue drawn by formula 4.2:

(4.2)

Wherein, ;

2. random fractal: use variable represent the coordinate that curved surface is put, represent this height on curved surface; Then FBM curve surface definition is on certain probability space, and index is a stochastic process

Stochastic variable meet:

(1) set up with probability 1 , namely process is from initial point; And be about continuous function;

(2) to aleatory variable , its two-dimentional increment obeying expectation is 0, and variance is normal distribution, its probability meet

(4.3)；

Adopting two kinds of fractal algorithms to generate the noise pattern of different pattern, is Diamond-Square algorithm and Perlin noise respectively:

Diamond-Square algorithm: segment each time in an iterative process and all can perform two step operations, Diamond step and Square step; Diamond step, for generating foursquare center, is got four summit calculating mean values, is added a random quantity, as the value of central point; Diamond step generates four triangular pyramids; Square step is then four summits of getting pyramid, and add that stochastic variable is as former foursquare four edges mid point with its mean value, the result segmented each time is consistent;

Perlin noise: for an amplitude and noise acoustic image, first normal vector and the gray-scale value of random direction is given to rounded coordinate point, the pixel value of any point is by the normal vector on four summits in this cell that place is the most contiguous in addition, decides with the direction vector being connected this point and four summits

If for the vector of required point in figure, nfor the stacking fold of noise, generally get ; If noise function wavelength is , then frequency , the gray-scale value of this point for:

(4.4)

Wherein, offsetfor skew measures-0.1, it is a noise function;

Pretreated method is taked to generate polytype Noise texture, avoid calculating the extra expenses increasing GPU in real time, in rasterization stage, pixel coloring device reads a texture maps comprising noise and elevation map merges, noise pattern carries out different proportion ground scaling and superposes simultaneously, expand infinitely-great Noise texture, eliminate with this periodicity and medelling effect drawn out.