WO2010126176A1

WO2010126176A1 - Method for partitioning region represented by contours into smaller polygonal zones and calculating digital elevation model data and geographic surface configuration data, and recording medium in which program for implementing method is recorded

Info

Publication number: WO2010126176A1
Application number: PCT/KR2009/002237
Authority: WO
Inventors: 최준수; 김천
Original assignee: Choi Joonsoo; Kim Choen
Priority date: 2009-04-28
Filing date: 2009-04-29
Publication date: 2010-11-04
Also published as: KR100916474B1

Abstract

Disclosed is a method for partitioning a region represented by contours into smaller polygonal zones which are able to express the geographic features of the region well and calculating digital elevation model data and geographic surface configuration data by using the same. Geographic partitioning is implemented by comprising the steps of: using contour sets to obtain a medial axis for each sub-domain that is a zone between adjacent contours; then partitioning the each sub-domain into smaller zones for the all segments of the medial axis. The method: calculates a geographic feature for each smaller zone after the geographic partitioning; and then calculates the elevation of each lattice point of a digital elevation model (DEM) by interpolating the elevations of the adjacent two contours based on an interpolation method which is predetermined according to the geographic feature of the smaller zone that includes the lattice points. Monotone rational Hermite spline or convex rational Hermite spline is used for the interpolation calculation, and the elevation of a desired point can be calculated from the contour sets. Surface reconfiguration data are obtained in the case of: partitioning the segments of the medial axis into the smaller zones; then straightening the segments of the medial axis; and then calculating three-dimensional coordinates for the polygons of the partitioned smaller zones. The DEM resulting from the aforementioned technology can be used for the orthorectification of high resolution satellite or aerial images. In addition, the DEM is applicable to various fields such as visualization of the surface of the earth, elevation calculation of a certain surface point, and the like.

Description

A recording medium that decomposes the area represented by the contour into polygonal detailed areas, calculates the digital elevation model data and the topographic surface composition data, and a recording medium on which a program for executing the method is recorded.

The present invention relates to the element technology required to visualize the topographic features of a region represented by contours, and more particularly, to decompose the region to be represented by the contour lines into detailed regions of the polygon mesh that can express the topographic features of the region. And how to generate three-dimensional digital elevation model (DEM) data for the area, or the data needed to construct the surface of the area, using the detail area obtained through the detailed area. .

In applications such as geographic information system (GIS), flight simulation, military training and operational simulation, three-dimensional computer games, virtual reality, etc., terrain data measured using satellite photographs, aerial photographs, ground surveys, etc. Requires skills to represent the terrain in the area. There are two general ways to represent the terrain. First, there is a method of expressing the terrain as a set of contour lines connecting points representing the same altitude, and second, there is a three-dimensional DEM that displays the altitude of the grid points divided into finely uniform intervals. In addition, in order to visualize the topographical surface of a region in three dimensions, a topographical surface construction technique for constructing the topographical surface of the region in close approximation to reality is also required.

In recent years, high resolution satellite images have been widely used in a variety of applications, including general mapping, photogrammetry, and geographic information system (GIS) data acquisition and visualization. Since satellite images are distorted by the tilt of the satellite sensors and topographical variations on their surface, the images must be orthorectified in order to be able to be used in each application. The satellite image represents an image distortion phenomenon in which the image pixel points are not placed at the correct position due to various factors such as the attitude of the satellite and the degree of tilt of the camera when the image is taken. Considering all the factors that cause such distortion, the process of removing the distortion of the image so that the position of each pixel point in the image coincides with the position on the topographical map by geometrically reconfiguring the environment as it was at the time of image shooting is called orthodontic correction. . Orthorectification of satellite images is usually obtained by resampling pixels of satellite images using rational polynomial coefficients (RPCs). However, when the glass polynomial coefficients are not given, the optimal glass polynomial coefficients are indirectly computed using the optimization technique of the glass polynomial coefficients using ground control points (GCPs) and DEMs that are well distributed on the surface.

DEM is a numerical representation of a surface topography or terrain surface in the form of a raster grid for a particular area, which consists of the topographic elevation for each point of the ground surface sampled at regular intervals. In addition to orthodontic correction of remotely detected images, DEMs are used in forms that help to systematically interrogate modeling, analysis, visualization, and topographical features. DEMs are created from various resources. The DEM can be generated by automatically extracting the DEM from stereo satellite images (scenes) or stereo digital aerial photographs. However, automatic extraction from stereo images is not always satisfactory in mountainous regions. In most cases, DEMs for mountainous regions are obtained through interpolation calculations of many existing digital contour maps made by direct survey of the terrain.

However, existing known methods for generating DEM from contours are difficult to reproduce the same result because the method of obtaining height data of a specific point is not clear and only outlined, and the approximate method of calculating height data between two contours is theoretical. Incorrect or weak grounds. Among the prior art, there is also known a technique for creating a new intermediate contour between existing contours and generating a DEM from the contours. However, this method is not only complicated, but also does not guarantee that it will always produce nice intermediate contours. In addition, a technique for calculating the height of an arbitrary point using a skeleton having a concept similar to the intermediate axis used in the present invention is also known. In this case, however, the interpolation method used is not applied consistently.

On the other hand, a general surface reconstruction of a three-dimensional object from a set of contours is to determine a surface, which approximates an unknown surface using geometric information within the contour. Contours are usually contours or contours obtained from cross-sectional data of computerized tomography (CT), magnetic resonance imager (MRI) or range sensors. General surface reconstruction from contours can be divided into three subproblems. One is the problem of correspondence, the other is the branching problem, and the last is the tiling problem. The correspondence problem is to determine which contours of a given level should be connected to which contours of adjacent levels. When a contour at a level corresponds to a single contour at an adjacent level, it is called a one-to-one correspondence. When a contour corresponds to two or more contours of adjacent levels, it is called a "one-to-many" response. Branching problems arise when there is a one-to-many or many-to-many response between adjacent levels. It is how to connect the corresponding contours of adjacent levels. The tiling problem is to construct a surface that connects the corresponding contour sets into a triangular mesh.

In reconstructing terrain surfaces from contours, only one-to-one or one-to-many responses occur because one contour encompasses one or more adjacent contours. More complex "many-to-many" responses or branching structures do not occur in this problem. The problem of extracting the DEM from the contours is simpler than the terrain reconstruction problem, since the approximate surface does not need to be constructed and the problem can be solved with correspondence and branching subproblems. However, the method of calculating the height value for any point used when extracting the DEM can be applied when calculating the height value of the vertex of the triangular mesh in the tiling problem of the surface reconstruction problem.

In the above-mentioned applications, the surface reconstruction of the terrain transforms the contour data or terrain data of the DEM into a polygon mesh that can represent the terrain as a set of three-dimensional triangles or rectangles, and then visualizes the terrain by rendering the polygons. do. As a method of converting a contour into a polygonal mesh, there is a method of connecting points on the contour to a triangle by the dilloni triangulation technique. (A) shows a sub-domain between contour c ₁ with vertices v ₁ , v ₂ , ..., v ₁₁ and contour c ₁ with vertices u ₁ , u ₂ , ..., u ₆ . This is an example of a triangular mesh expressed by the Deloni triangulation technique. When the height h ₁ to height h ₂ of the contours of the contour c ₂ c ₁ in other cases, such a triangular mesh is prevent good representation of the terrain. For example, a triangle such as a triangle v ₁ u ₁ v ₇ , v ₁ u ₁ v ₂ is connected to two contours of different heights to represent an inclined plane. However, the triangles v ₁ v ₇ v ₉ , v ₁ v ₉ v ₁₁ , v ₁₁ v ₉ v ₁₀ connect the points of the contour line c ₁ to represent the flat surface, which is a good representation of the topographical features with curvature between the two contours. can not do. In (b) of FIG. 1, one of three vertices of a triangle is a case where the sub-domain is triangulated to be located at another contour line. In this case, as illustrated in the drawings, an elongated unreal triangle may be formed, or in the case where an unusual contour is given, triangulation may not be possible.

Topographical features of some regions include peaks, pits, hillsides, ridges, valleys, and canyons. DEM data or terrain surface reconstruction data obtained by processing contours need to reflect these topographical features of the terrain as closely as possible. In order to obtain DEM data or terrain surface reconstruction data, it is known that it is effective to divide the terrain represented by the contour into smaller regions. However, according to the analysis of the present inventors, even if the existing DEM data generation method or terrain surface reconstruction method that approached, it is difficult to generate DEM data or terrain surface reconstruction data that properly express the actual characteristics of the terrain. The main reason for this is that it is not based on the terrain segmentation method that can express the topographic features well.

The present inventors have made various efforts and explorations to find an improvement of the terrain segmentation method, and as a result, in order to better express the topographical characteristics of the terrain (three-dimensional surface of the object) represented by the contour set, the terrain is divided into a plurality of small regions. In other words, it was found that it is necessary to decompose the polygon into planes, but to find the critical points that can express the topographical features and divide them using the critical points. Accordingly, a first object of the present invention is to provide a new terrain segmentation method for dividing a terrain represented by contour lines into a plurality of small regions capable of expressing actual geographical features.

It is a further object of the present invention to provide a new method that can generate DEM data using such a terrain segmentation method, which very closely reflects the actual topographical features of the terrain from a given set of contours.

It is still another object of the present invention to provide a method that can calculate the height for an arbitrary position very close to the actual height, given the set of contour lines, using the terrain segmentation method as described above.

It is still another object of the present invention to provide a method of reconstructing a three-dimensional surface that very closely represents the actual features of the terrain from a given set of contours using the terrain segmentation method as described above.

It is yet another object of the present invention to provide a method for more accurately calculating the elevation of a desired point in a region represented by contour lines.

Furthermore, another object of the present invention is to provide a computer readable recording medium recording such application programs that can be executed by a computer or other computing device.

According to an aspect of the present invention for achieving the above object, a method performed by a computer and decomposing a terrain represented by a contour set into detailed regions, by sub-domains, which are regions between adjacent contour lines using the contour set. Obtaining an intermediate axis; And decomposing each sub-domain into zones for all segments of its intermediate axis, wherein the segment of the intermediate axis is a collection of all center points of the largest circles that contact the same two active boundary elements, An active boundary element is provided that decomposes the terrain to be represented by a set of contour lines into subregions, characterized by being defined by corners and concave vertices of the sub-domain that abut the maximum circle.

The step of obtaining the intermediate axis comprises: creating a contour tree representing correspondences between contours from data relating to the contour set; Calculating all sub-domains that are regions between adjacent contour lines using the contour tree; And calculating an intermediate axis for each sub-domain.

The boundary of the subregion is the point of contact between the center point of two largest circles located at each of the important point where the new segment starts and the important point which ends each time while tracking the intermediate axis, and the two active boundary elements which each of the two maximum circle contacts. It is defined as the radius segment of the largest circle connecting. The key points include, among the points on the intermediate axis, end points that are points of the intermediate axis that intersect the corner boundaries of the polygon, junction points that are the center points of the largest circle in contact with three or more active boundary elements, and two points. There are three types of transition points, which are the center points of the maximum circle in contact with the active boundary element, the center point where only one of the two active boundary elements changes.

According to another aspect of the present invention for achieving the above object, a method of calculating an elevation of a desired point using a set of contour lines, which calculates an approximate inclination line passing through the point for which the elevation is to be obtained. step; And calculating the elevation of the point as the elevation of the point by interpolating the distance between two adjacent contours of the point along the approximate inclination line and the elevation of the two contours adjacent to the point using a predetermined interpolation method. There is provided an elevation calculation method comprising a.

In the above elevation calculation method, an approximate slope line passing through the point is the maximum circle contacting the contact point p ₁ which meets when the vertical line is lowered from the two adjacent contour lines c ₁ and c ₂ to the contour line c ₁ close to the point. When contacting the other contour line c ₂ at the point p ₂ , the first line segment connecting the contact point p ₁ and the center point of the largest circle and the second line segment connecting the center point of the maximum circle and the contact point p ₂ It is determined by the line connecting. If there is no other contour line inside the contour line including the point, the approximate slope line passing through the point is different from the maximum circle in contact with the contact point p ₁ that meets when the vertical line is lowered to the point on the contour line closest to the point point. Also at point p ₂ , when contacted with the contour line, the first line segment connecting the contact point p ₁ and the center point of the largest circle and the second line segment connecting the center point of the maximum circle and the contact point p ₂ are determined. . Further, among the two adjacent contour lines c ₁ and c ₂ , the first maximum circle MC ₁ contacting the contact point p ₁ that meets when the vertical line is lowered to the contour line c ₁ close to the point is the other contour line c ₂ . If not, the approximate inclination line passing through the point is i) the first largest circle MC located on the intermediate axis MX of the contact p ₁ and the two adjacent contour lines c ₁ , c ₂ . ₁ ) a first line connecting the center point p ₅ of the center, ii) from the center point p ₅ of the first largest circle MC ₁ along the segment of the dangling subtree in the cycle direction of the intermediate axis MX, to the rest of the contour (c ₂₎ a second maximum circle (MC ₂₎ second line segments and, iii consisting of segments of dangling subtree to the center point p ₃ of the tangent at the contact point p ₂ in) Finally the second 2 A line connecting the center point p ₃ of the maximum circle MC ₂ and the third line connecting the contact point p ₂ . It is decided.

The elevation calculation method further includes the step of obtaining the topographical feature of the detailed area to which the point to which the elevation is to belong. If the terrain features correspond to hills, valleys, or ridges and the elevations of two adjacent contours are different, it is preferable to apply the forged glass Hermite spline interpolation method during the elevation interpolation. If the feature of the terrain is the same as the canyon, or if the heights of two adjacent contours are equal, or if the height of the two adjacent contours defining the ridge is the same, or the feature of the terrain may be represented as convex (or concave). In the case of terrain such as a summit or a pond, it is desirable to apply a convex hullite spline interpolation method that reflects such terrain features well in the elevation interpolation.

According to another aspect of the present invention for achieving the above object, a method for generating a DEM using a contour set, which is performed by a computer, by using the contour set to form a sub-domain intermediate axis that is an area between adjacent contour lines Obtaining; Decomposing each sub-domain into zones for all segments of its intermediate axis; Calculating a topographic feature of each detail area; And calculating the elevation of each grid point of the DEM by interpolating the elevations of two adjacent contour lines using a predetermined interpolation method according to the terrain feature of the detailed area to which the grid point belongs. And a set of all the center points of the largest circles in contact with the two active boundary elements, wherein the active boundary elements are defined as corners and concave vertices of the sub-domain that abut the maximum circle.

In the DEM generating method, the terrain feature calculation step of the detail region inspects the structure of the intermediate axis, and (1) when the segment belongs to the dangling subtree, the detail region of the segment is a relative elevation of adjacent boundary contours. Determined as part of any one of the ridges, valleys, and canyons, and (2) if the segment belongs to a cycle, the detail area of that segment shall be determined as hilly, ridge, or canyon, (3 ) If a segment belongs to a route, the detail region of that segment is determined by one of ridges, canyons, and hills, and (4) if the sub-domain does not have internal boundary contours, the region of that sub-domain is Decide on a summit or pool as a whole. Wherein the cycle is a subset of intermediate axis edges with a closed shape, the dangling subtree is a maximum connected intermediate axis edge subset that does not belong to any cycle of the intermediate axis or does not belong to a path connecting the cycle, and the path is a cycle Subset of mid-axis edges connecting them together.

In the DEM generation method, the step of calculating the elevation of the grid point, the step of calculating the approximate slope line passing through the grid point to obtain the elevation; And interpolating a distance ratio between two adjacent contour lines along the approximate inclination line based on the corresponding grid point and an elevation of the two adjacent contour lines using the interpolation method. When the approximate slope line is drawn at the maximum circle at any point on the intermediate axis, the approximate inclination line falls to two contacts (p ₁ , p ₂ ) with the two adjacent contour lines (c ₁ , c ₂ ) at the center of the maximum circle. As a line connecting line segments, it is defined as a line passing through the grid point for which the elevation is to be obtained.

In the DEM generation method, if the topographical feature of the detailed region to which the lattice point belongs is hilly, valley, or ridge, and the elevations of two adjacent contours are different from each other, a forged glass Hermite spline interpolation method is applied during interpolation. It is preferable. In addition, even when the topographic features of the detailed region to which the grid points belong are the same as the canyon, or the heights of two adjacent contour lines are the same, or the ridges, the heights of two adjacent contour lines defining the ridges are the same, or the topographic features are convex (or In the case of terrain, such as a summit or a pond, which may be represented by concaveness, it is desirable to apply a convex hullite spline interpolation method that reflects such terrain features well in the elevation interpolation.

According to another aspect of the present invention for achieving the above object, a computer-implemented method of reconstructing a three-dimensional surface of a terrain represented by a contour set, using the contour set is a sub-domain that is an area between adjacent contour lines. Obtaining a star intermediate axis; Decomposing each sub-domain into zones for all segments of its intermediate axis; If a segment of the intermediate axis is represented by a parabola, it straightens by connecting the main points at both ends of the parabolic segment in a straight line; And calculating three-dimensional coordinates of the decomposed subregion polygons, that is, two-dimensional plane coordinates and an elevation, wherein the segment of the intermediate axis is a set of all center points of the maximum circles contacting the same two active boundary elements. A boundary element is provided, wherein the three-dimensional surface reconstruction method of the contour representation terrain is defined by the edges and concave vertices of the sub-domain that abut the maximum circle.

In the surface reconstruction method, it is preferable that the quadrangles of the disassembled subregions further include subdividing into two triangles.

Also in the surface reconstruction method, the elevation calculation step includes calculating an approximate slope line passing through the elevation calculation point; And calculating the elevation of the point as the elevation of the point by interpolating the distance between two adjacent contours of the point along the approximate inclination line and the elevation of the two contours adjacent to the point using a predetermined interpolation method. It is preferable to include.

According to another aspect of the present invention for achieving the above object, a computer-implemented method of reconstructing a three-dimensional surface of a terrain represented by a contour set, using the contour set is a sub-domain that is an area between adjacent contour lines. Obtaining a star intermediate axis; Decomposing each sub-domain into zones for all segments of its intermediate axis; Calculating a topographic feature of each detail area; Calculating the elevation of each grid point of the DEM by interpolating the elevations of two adjacent contour lines using a predetermined interpolation method according to the terrain feature of the detailed area to which the grid point belongs; Decomposing the entire grid point of the DEM result data into a mesh of polygons such as a plurality of triangles or squares by connecting three or four adjacent grid points; And obtaining and allocating planar coordinates and elevations, which are three-dimensional coordinate values of vertices of each polygon, from the DEM result data, wherein the segment of the intermediate axis is a set of all center points of the maximum circles contacting the same two active boundary elements, The active boundary element is provided with a three-dimensional surface reconstruction method of contour representation terrain, characterized in that it is defined by the corners and concave vertices of the sub-domain that abut the maximum circle.

According to another aspect of the present invention for achieving the above object, there is provided a computer-readable recording medium recording a program for executing each step of each method as described above.

In order to extract the terrain-related data from the contour line, such as obtaining the height data of a specific position from the contour line data, it can be simply extracted directly from the contour data. However, it is more effective to first analyze the features of the terrain that the contour data means or to represent and then indirectly extract the desired data from the features of the terrain. For example, if you find critical points that represent topographical features well from the contours, and find a break line that connects them, the topographic features that are represented by lines such as mountain peaks, ridges, or valleys of the terrain. Can be expressed well. In previous studies, the line elements were extracted from the contours and the terrain data were indirectly extracted. However, the terrain segmentation method proposed by the present invention is an extension of such a concept, and the polygons are represented as planes by decomposing the contour region from the contour data into polygons connecting the critical point on the dividing line and the points on the contour line. It can be used to express geographic features such as hills, ridges, valleys, and canyons. Therefore, compared to the existing technology, it is possible to express various topographical features closer to reality.

In the elevation calculation technique at a particular location, interpolation was previously used that is proportional to the vertical distance from that particular location to two contours. This method does not take into account the inclination lines present on the terrain and also has the disadvantage of being inconsistent in all positions. On the other hand, the elevation calculation technique proposed by the present invention more realistically calculates the slope line of the terrain in the region between the contour lines and uses an interpolation method proportional to the distance on the slope line. Elevation can be calculated in a conventional way. Using this method, contours can be used to calculate approximate elevation data for a specific location more accurately than conventional methods, without actually going to the field and using techniques such as triangulation or Global Positioning System (GPS). .

On the other hand, the DEM generation method proposed by the present invention divides the area represented by the contour into small sub-regions, each sub-region to have such features as hilly, ridge, valley, canyon, summit or pool Based on the decomposition method, it is possible to calculate DEM data (elevation data for grid points) that reflects the actual topographical features of the area. In addition, it is possible to calculate more accurate height data for a specific position by applying an interpolation method for the terrain feature when calculating the elevation. For example, for hills, ridges, and valleys, forged glass Hermite interpolation can be used to calculate height data, and in the case of tops or puddles, convex glass Hermite interpolation can be used to calculate more realistic height data. . By calculating the height of grid points such as the checkerboard extracted in this way, the DEM can be made to be a more realistic surface by changing the height smoothly rather than changing the height rapidly.

Finally, in order to reconstruct the three-dimensional surface of the terrain using the contour data, first, in order to create a three-dimensional surface from the contour data, the contour is converted into a polygon mesh that can represent the terrain as a set of three-dimensional triangles or squares. Then, the terrain is visualized by rendering the polygons. Existing methods of converting contour lines into polygon meshes have disadvantages such as unrealistic polygon meshes or unrealistic shapes of squares or triangles. However, if the polygon mesh is made by the method proposed in the present invention, it is possible to reconstruct a three-dimensional surface of a shape that conforms to the features of the terrain a little more than the conventional method.

(A) of FIG. 1 is an example in which a sub-domain between two contours is represented by a triangular mesh by a dilloni triangulation technique, and FIG. 1 (b) shows that one of three vertices of a triangle is located at another contour. Example of triangulation of a domain.

2 shows a contour tree and sub-domains of an exemplary set of contours.

3 and 4 are views for explaining the concept of the intermediate axis and the detail region of a polygon having a simple polygon and a hole.

5 and 6 are views for explaining a method for calculating the intermediate axis of the polygon.

FIG. 7 shows a method for determining the end point of the intermediate axis segment in the case of FIG. 6 (a).

FIG. 8 shows a method for determining the end point of the intermediate axis segment in the case of FIG. 6 (b).

FIG. 9 shows a method for determining the end point of the intermediate axis segment in the case of FIG.

FIG. 10 illustrates the decomposition of a particular sub-domain into subregions based on a segment of a typical intermediate axis with cycles and dangling subtrees.

FIG. 11 is an example plotting a planar graph showing the mid-axis of a sub-domain with four contours c ₂ , c ₃ , c ₄ , c ₅ inside contour c ₁ .

FIG. 12 is a diagram for explaining a conventional linear interpolation method for calculating an elevation of an arbitrary point P between two contour lines.

FIG. 13 is a diagram for describing a method of obtaining an inclination line and linearly interpolating an elevation of a desired point using the gradient line.

FIG. 14 is a diagram for explaining a method of performing elevation calculation of a hillside region by the forged hermite interpolation method.

15 is a diagram for explaining a tilt calculation method at an intersection point between an intermediate slope line and a contour line.

16 is a diagram for explaining an interpolation method relating to a top and a pool.

17 and 19 are diagrams for explaining elevation calculations when the terrain feature is a ridge.

18A and 18B are views for explaining elevation calculations when the terrain feature is a valley.

18A, 18C, and 20 are views for explaining elevation calculations when the terrain feature is a canyon.

21 illustrates a configuration of an application program implementing the DEM automatic generation method according to the present invention and a configuration of an apparatus for executing the same.

Fig. 22 is a flow chart showing the flow of a procedure for generating DEM result data by this application program processing DEM input data.

FIG. 23 is a diagram for explaining a calculation method of grid point coordinates for a set of contour lines.

FIG. 24A is a graphic representation of an exemplary contour line used as the source data of the present invention, and FIG. 24B is a graphic representation of DEM result data generated by the method proposed by the present invention.

FIG. 25 illustrates a straightening of all segments of the intermediate axis illustrated in FIG. 4, and FIG. 26 illustrates a state in which all quadrangles of FIG. 25 are divided into two triangles.

27 schematically shows the data structure of a set of contour lines.

Hereinafter, with reference to the accompanying drawings will be described in detail for the practice of the present invention.

1. Terrain segmentation method

(1.1) source data

Contours can be obtained from stereo satellite images, stereo aerial photographs, and other direct surveys. The source data of the present invention is such contour data. FIG. 24A is a graphical representation of exemplary contour lines used as source data of the present invention. For reference, Figure 24 (b) is a graphical representation of the DEM result data generated by the method proposed by the present invention.

Contours as source data are approximated as simple polygons that connect adjacent vertices with line segments, and all points on the contour have the same height value. In order to define the contour set, an integer n representing the total number of contour lines constituting the contour set and a definition of each of the n contour lines are required. Each contour is defined by the number k ₁ of the total vertices to which the contours are to be expressed, the height h ₁ of the contours, and the two-dimensional coordinates of k ₁ vertices. 27 illustrates a data structure of such a contour set.

Since the contours are represented by polygons, the vertices of the contours can be known by connecting all the vertices of the polygons. Therefore, knowing the total number of vertices and coordinates of a contour can trace the contour of the contour. Terrain modeling contours have certain limitations with regard to the geometry allowed. For example, two contours should not overlap each other. Terrain is modeled in the form of a complete nesting hierarchy of contours. That is, some contours encompass any number of other contours, but are completely contained within only one other contour of the next higher level, and no contours intersect the other contours.

(1.2) Method of Terrain Segmentation

(1.2.1) In order to create a polygon mesh that can accurately represent the terrain to be represented by contours, it is necessary to find the topographic features that are indirectly represented in the contour lines and use them to make a polygon mesh. For example, if you find critical points that represent topographical features well from the contours, and find a break line that connects them, the topographic features that are represented by lines such as mountain peaks, ridges, or valleys of the terrain. Can be expressed well. Also, by creating polygons connecting the critical points on the dividing line and the points on the contour line, it is possible to express the topographical features of hills, ridges, valleys, and canyons where these polygons are represented by faces.

The terrain segmentation approach proposed by the present invention is a polygonal mesh, that is, a region having geographic features such as hills, ridges, valleys, and canyons represented by planes of sub-domains, which are divided by contour lines. Decompose into rectangular polygons, but find the critical points from the contours that represent the topographic features of the sub-domains, and break them up using break lines connecting the critical points. to be. Each polygonal area thus broken down is called "zones." After all, the detail region is a polygon region connecting the critical points on the dividing line connecting the critical points and the points on the contour line. This will be described in more detail below.

(1.2.2) In order to find the critical point, the first thing to do is to build a contour tree representing the correspondence between contours from the data on the set of contours entered. Contour trees are data structures that describe hierarchical relationships between contour lines. The corresponding problem of determining adjacent contours to be connected together can be solved by utilizing the contour tree. 2 shows a contour tree and sub-domains of an exemplary set of contours. Total contour lines c ₁ , c ₂ ,. For example, we can determine the correspondence between contour lines by considering the inclusion relation between c and ₁₄ . The correspondence can then be used to build a contour tree for all contours. Inclusion relationships between the contour lines can be determined using the trajectory (coordinate) and height information of each contour line. For example, n contour sets c ₁ , c ₂ ,... Given the source data (the total number of contours, the number and coordinates of the total vertices of each contour, and the height of each contour) for c _n , the trajectories of the contours are calculated using the number of vertices and the coordinates of each contour. . The traces of all the contour lines obtained can be expressed as shown in FIG. Using the obtained contour trajectory and the height of each contour line, it is possible to grasp the containment relationship between the contour lines, that is, the correspondence relationship between the contour lines, and create a contour tree according to the correspondence relationship. The tree structure of FIG. 2B is a contour tree with respect to the contour line of (a).

(1.2.3) Once the contour tree is found, it is used to calculate all the sub-domains. The sub-domain Ω _i is defined as the area between each adjacent contour in the contour tree. Thus, a sub-domain is represented as a polygon with no holes inside or with one or more holes. That is, the outer boundary contours defining the sub-domains can be viewed as a polygon, and the inner boundary contours within them can be viewed as holes in the polygon. Generalizing this, suppose that contour c _i has an elevation value h _i and c _i lies inside some c _j (where j <i for i ≥ 2). The outermost contour c ₁ encompasses all other contours and it is a domain defined by two open sets of planes, the unqualified domains Ω ₀ and c ₁ , which in turn are n mutually open sets by different contours. (disjoint open sets) are separated by Ω _i . Each sub-domain Ω _i has an outer boundary contour c _i and k _i immediately adjacent inner boundary contours c _j (where j> i). For example, in FIG. 2B, the sub-domain Ω ₁ is an approximate polygon defined by the outer boundary contour c ₁ and is a polygonal region having one polygon hole defined by the inner boundary contour c ₂ therein, Sub-domain Ω ₇ is an approximate polygon defined by the outer boundary contour c ₇ and is a polygonal region with three polygonal holes defined therein by the inner boundary contours c ₈ , c ₉ and c ₁₂ .

For simplicity, assume that all k _i inner boundary contours have the same elevation and that the elevation is different from the elevation of the outer boundary contour c _i .

If k _i = 0, the sub-domain Ω _i has a local pole and the sub-domain is normal or a pool. In FIG. 2B, for example, the sub-domains Ω ₈ and Ω ₁₄ correspond to this. If c _i is the inner boundary of the sub-domain Ω _j , consider them by considering the elevation relationship between h _i and h _j . If h _i > h _j , Ω _i is normal and else it is a pool.

If k _i = 1, contour c _i corresponds to one contour c _{i + 1} . In FIG. 2B, for example, the sub-domain Ω ₁ corresponds to this. In this case, the sub-domain Ω _i is a mixture of hills, ridges, valleys or canyons. If h _{i + 1} > h _i , as in the summit and pool cases, the sub-domain Ω _i is a mixture of outbound hills and ridges; in other cases, inbound hills, valleys and canyons .

If k _i > 1, contour c _i corresponds to k _i contours. This is equivalent domain Ω ₇ - In Fig. 2 (b), for example sub. This case is called one-to-many branching and it forms a complex terrain. In that case the sub-domain Ω _i is a mixture of hills, ridges, valleys and canyons.

(1.2.4) After all sub-domains have been found, calculate the intermediate axis for each sub-domain. First, the method of calculating the intermediate axis from the outer boundary contours and the inner boundary contours that define the sub-domain Ω _k considers the outer boundary contour as a polygon, the inner boundary contour as the hole of the polygon, and the polygon Calculate the intermediate axis of. A sub-domain with only outer boundary contours can be viewed as a simple polygon and its mid-axis calculated.

3 exemplarily shows the shape of the intermediate axis of the simple polygon P. FIG. The medial axis (MX) of a simple polygon P is the trace of the center points of the circles (called 'Maximal Circles (MC)') contained within the polygon P and tangent to the polygon P at two or more points. The relationship between these maximum circles is such that no maximum circle is a complete subset of the other maximum circles in the polygon P. The middle axis MX is a simple tree graph consisting of vertices and edges connecting the vertices, where the edges are parabolic arcs or straight line segments. In the case of a polygon including at least one hole therein, the graph of the intermediate axis may further include cycles surrounding each hole, and in some cases, may further include a path connecting each cycle. The middle axis MX is closely related to the generalized Voronoi diagram of the corner segments of polygon P and their endpoints. A generalized Voronoi diagram of a polygon's edges and concave vertices (RV, a vertex larger than 180 degrees in size) differs from the polygon's mid-axis at a portion near the concave vertex of the polygon. . The intermediate axis does not have any edges connected to the concave vertices RV, whereas in a typical Voronoi diagram there is always an edge connected to the concave vertex RV. There is a generalized Voronoi diagram of n line segments and an optimal O (nlogn) algorithm that computes the intermediate axis. In the case of a simple polygon with holes inside the polygon, the middle axis of the polygon can be calculated in O (n (logn + h)) time. Where n is the number of edges of the polygon and the hole, and h is the number of holes. The intermediate axis of a polygon with holes is a general planar graph with parabolic arcs or straight line segments. Computing the intermediate axis of a simple polygon or a polygon having several holes therein is called a medial axis transformation (MAT).

Regarding the server-domain represented by the approximate polygon and the intermediate axis (MX) creation of each sub-domain, the concept of key points and segments of the intermediate axis is introduced as follows.

Key points (KP): The corners and concave vertices of polygons that contact the maximum circles (MCs) are called active boundary elements, and are located on the middle axis (MX) depending on the number of active boundary elements that contact the maximum circle. Classify the types of points. Among the points on the intermediate axis, the points on the intermediate axis that intersect the edge boundary of the polygon are called end points (EP), and the center points of the maximum circles MC which contact three or more active boundary elements are the junction points (JP). ), The center point of the maximum circle (MC) in contact with the two boundary elements is the regular point (RP), and in particular, the normal point at which one of the two active boundary elements changes is called the transition point (TP). do. And these three points are called the important point (KP) of the intermediate axis (MX).

Segment of the intermediate axis: The segment Sg of the intermediate axis MX is the maximum subset of the intermediate axis that is uniquely associated with two distinct active boundary elements. That is, each segment of the intermediate axis is a set of all the center points of the maximum circles MC that contact the same two active boundary elements. In Figure 3 a segment of a transition point TP ₁ and TP ₂ in both end points is a set of the largest circle center to the edges e and the RV apex of the polygon in the active boundary element. Both endpoints of each segment Sg are important points.

Based on these concepts, a corresponding intermediate axis can be generated for each sub-domain Ω _i defined in the hierarchy of terrain contours. In this case, the outer contour of each sub-domain Ω _i is a polygon defining the intermediate axis, and the inner contours of each sub-domain Ω _i are holes in that polygon. 5 and 6 are views for explaining a method for calculating the intermediate axis of the polygon. A method of calculating the intermediate axis of the polygon will be described in detail with reference to this. The middle axis of the polygon is made in the form of a normal graph. As illustrated in FIG. 3, the middle axis of the polygon without a hole is in the form of a tree graph without a cycle, and the intermediate axis of the polygon with a hole is in the form of a general graph in which a cycle exists as illustrated in FIG. 4. The vertices of this graph are the midpoints of the middle axis, and the edges of the graph are segments that connect the two important points, with the active circle element of the largest circle centered on the segment. For example, FIG. 5 (b) shows a graph showing the middle axis MX of a given polygon.

Suppose that a given circle is given a maximum circle with a center point and an active boundary element of that maximum circle. In the first case, one of the convex vertices (v ₁ , v ₂ , v ₃ , v ₄ , v ₆ , v ₇ ) of the polygon begins with the key points. Starting from this center point, the intermediate axis is calculated by tracing the center point of the maximum circle that is in contact with the active boundary element and the size of the circle changes continuously. When tracking a center point, we find one segment that connects the key point that we started with by calculating the next key point. If the newly calculated next critical point is the end point, we finish tracking the center of the maximum circle. For example, in the case of tracking from the important point j ₂ to the end point v ₁ in FIG. 5 (b), a segment connecting j ₂ and v ₁ is created and the tracking ends. If the newly calculated next critical point is the transition point (t ₁ , t ₂ ), one of the active boundary elements of the largest circle is changed and another new segment starts from this transition point. For example, in the case of tracking along the transition point t ₁ at the critical point j ₂ in FIG. 5 (b), the active boundary element changes from e ₁ , e ₅ to e ₁ , v ₅ , and a new segment Begins. If the newly calculated next critical point is the connection point (j ₁ , j ₂ , j ₃ , j ₄ ), the maximum circle is in contact with three or more active boundary elements, in which case one or more segments are connected at the connection point. We start to create a new segment using the two active boundary elements that make up each segment. The remaining segments start with the segment you just started tracking, then track all the new segments, then start tracking. For example, in the case of tracking from the critical point t ₂ to the connection point j ₃ in FIG. 5B, two additional segments are connected at the connection point j ₃ . Either first track the segment associated with v ₄ or the segment associated with j ₄ , then track all segments that can be tracked from that segment, then track the other segment. With a slight modification of the above method, you can calculate the intermediate axis for polygons with holes in the same way. In other words, a cycle surrounding a hole can be created in the mid-axis graph for a polygon with a hole. In this case, while tracking a segment, you no longer need to track the next segment if the endpoint of that segment returns to any connection point it has ever tracked.

There are three types of segments created by two active boundary elements, as shown in FIG. As illustrated in (a) of FIG. 6, when the active boundary element is two vertices, the segment is in the form of a straight line. As illustrated in (b), one active boundary element is a vertex and the other active boundary element is In the case of polygon edges, it becomes parabolic, and as illustrated in (c), the segments created by two polygon edge active boundary elements are in the form of straight lines.

In creating the intermediate axis, the end point of each segment is i) the point where the extension line of the segment meets a line extending perpendicular to the polygon edge at the concave vertex of adjacent active boundary elements, or ii) at the convex vertex of the adjacent active boundary elements. The point that meets the extension line that bisects the vertex angle is determined to meet first. FIG. 7 shows a method for determining the end point of the intermediate axis segment in the case of FIG. 6 (a). In FIG. 7A, when the active boundary element is two concave vertices v _i and v _j , the end point of the middle axis segment MX _m is the extension line of the segment i) the active boundary element edge e _{i + 1} at the concave vertex v _i . in that the lines and meets vertical extension to m _i and ii) a concave vertex v point for the active boundary element edge e _{j + 1} from _j to see the vertically extending line and determines a point of m _j before meeting the end of the intermediate shaft MX _m . In the case of FIG. 7B, the end point of the intermediate axis segment MX _m is determined to be m _i in the same manner. FIG. 8 shows a method for determining the end point of the intermediate axis segment in the case of FIG. 6 (b). In (a) of FIG. 8, since the middle axis segment MX _m meets the line extending vertically with respect to the edge e _{i + 1} at the concave vertex v _i , the point m _i _precedes the point m _j that meets the extension line of the bisector of vertex v _j . Determine the point m _i as the end point of the intermediate axis segment MX _m . Similarly, the end point of the intermediate axis segment MX _m in FIG. 8 (b) is determined as m _i , and the starting point of the intermediate axis segment MX _m in FIG. 8 (c) is m _i and the end point is determined as m _j , FIG. 8. In (d), the end point of the intermediate axis segment MX _m is determined by m _j . FIG. 9 shows a method for determining the end point of the intermediate axis segment in the case of FIG. An intermediate shaft that is an extension of the segment MX _m convex vertex v in _i points m _i is the active boundary element concave vertex v _j of intersection with each bisector h _i1 of the meeting and the vertical extension h _j2 for the edge e _i connected thereto so ahead of m _j Point m _i is determined as the end point of the intermediate axis segment MX _m .

(1.2.5) After finding the intermediate axis for the entire sub-domain Ω _i , decompose each sub-domain into zones for all segments of its intermediate axis.

The radius of the maximum circle connecting the center of the two largest circles located at each of the starting point and ending point of each new segment and the contact of the two active boundary elements that each of the two maximum circles are in contact with, along the mid-axis. It becomes the boundary line of this subregion. According to the terrain segmentation method of the present invention, the sub-domain, which is the object to be divided, is divided into small regions (detailed regions) corresponding to each segment of the intermediate axis by finding its intermediate axis, thereby matching the face (detailed region) to a line (segment). Region division in the form of a form is called quasi-dual medial axis decomposition. Referring to FIG. 3, if the transition point TP ₂ is encountered while tracking the intermediate axis MX toward the transition point TP ₁ past the connection point JP ₂ , the transition point TP ₁ , which is the starting point of the new segment Sg, is met. The transition point TP is the maximum circle MC ₁ in the _first two active boundary element e ₁ and RV ₁ two radii of the largest circle MC ₁ to connect the contact point of the line R ₃ and R ₄ are one of detail area 30 It becomes a boundary line. Also keep track of the intermediate axis (MX) past the transition point TP ₁ to meet the new transition point TP ₂ . Across the transition point TP ₂ maximum circle MC ₂ are two active boundary element e ₁ and RV ₁ two radial lines R ₁ of the maximum circle MC ₂ connecting the point of contact with, and R ₂ is a detail area 30 of the It is the border of the side. In other words, the four radial segments R ₁ , R ₂ , R ₃ , R ₄ and the two active boundary elements e ₁ which fall on two active boundary elements e ₁ and RV ₁ at both ends TP ₁ and TP ₂ of the segment Sg And the region 30 surrounded by RV ₁ is referred to as a detailed zone corresponding to the segment Sg. As such, the closed area by the subregion boundary line and the two active boundary elements created at the critical point of both ends of the segment becomes the subregion corresponding to the segment.

Fig. 4 shows the decomposition of a sub-domain between the polygon with holes, i.e., the inner contour 42 and the outer contour 44, into subregions, in which case the subregions are divided in the same way. That is, the radial line made in both end JP ₁ and JP ₂ of the segment Sg to the maximum circle MC ₁ and MC contact ₂ is in contact with each of the two active boundary element e ₁ and e ₂ of the respective end R _5, R _6, R ₇ , R _{8 together} with the two active boundary elements e ₁ and e _{2 become the} subregion 40 of the segment Sg. FIG. 4 illustrates an example in which sub-domains with holes are decomposed into subregions according to such subregion determination criteria, and both end points of all segments belonging to the intermediate axis of the polygon with holes and the two active boundaries of the segments. Shows all radial segments connecting points that touch the element. The radius segment (the detail region boundary line) and the active boundary element divide the inside of the polygon into detailed regions of the segments.

2. Determination of topographic features of subregions

Within the sub-domain Ω _i may include various kinds of terrain features. Therefore, the sub-topographic characteristics such as to determine the geographical characteristics of the domain Ω _i are each subdivision ridge (ridges), the valley (valleys), hills (hillsides), Canyon (canyons), pool (pits), the top (peaks) It is necessary to classify which of these is applicable. Therefore, after sub-domain Ω _i is divided into multiple subregions using intermediate axis segments, the sub-domain OMEGA _i is matched with the corresponding terrain feature for each subregion. In the method of the present invention, identification of the topographical feature is accomplished by examining the structure of the intermediate axis MX by pseudo-dual intermediate axis division of the target terrain. As illustrated in FIG. 10, the intermediate axis MX of the sub-domain Ω _i , defined as the area between the outer contour line 56b and the inner contour line 56a, is a general planar graph with cycles 52.

Here, cycle 52 refers to a subset of intermediate axis edges having a closed shape, and the path is a subset of intermediate axis edges connecting the cycles to each other. The maximum connected intermediate axis edge subset that does not belong to any cycle of the intermediate axis MX or does not belong to the path connecting the cycles is also called dangling sub-tree 50. FIG. 11 is an example of plotting a planar graph showing the mid-axis of a sub-domain with four contours c ₂ , c ₃ , c ₄ , c ₅ inside contour c ₁ . As such, the segments of the intermediate axis belong to the cycle of the intermediate axis graph, belong to the path connecting the cycles, or belong to the dangling sub-tree. 11 has a four-cycle _{_{_{CL 1, CL 2, CL 3}}} , CL 4 to the intermediate shaft graph, P ₁ is a path (Path) to connect the cycles CL ₂ and CL _4. The detailed areas corresponding to the segments of the intermediate axis are classified as follows according to whether the segments belong to a cycle, a dangling subtree, or a path connecting the cycles.

(2.1) If a segment belongs to a dangling subtree, the detail region of the segment may be part of a ridge, part of a valley, or part of a canyon, depending on the relative elevation of the bounding contours. In addition, the sum of the corresponding detail regions of all segments in one dangling subtree becomes a continuous ridge or a continuous valley or canyon. For example, in FIG. 10, when the elevation of the inner contour 56a is higher than the elevation of the outer contour 56b, the detail region 50 of the dangling subtree is ridge, whereas in the opposite case, the dangling subtree Detail area 50 is a valley. In FIG. 10, when the elevation of the inner contour line 56a is the same as the elevation of the outer contour line 56b, the detail region 50 of the dangling subtree becomes a ridge or a canyon. The distinction between these two cases cannot be judged by the boundary contours alone, but can be judged by the difference between the altitude of the contour line just outside the contour boundary and the elevation of the inner contour line just inside the contour boundary. For example, if the elevation of the contour just outside of the outer boundary contour is high, the detail area of the segment is a canyon, otherwise it is a ridge.

(2.2) When a segment belongs to a cycle CL, the subregions of that segment are hillside, ridge, or canyon. In FIG. 10, if there is a difference between the elevations of the inner contour line 56a and the outer contour line 56b, the detailed region 54 is an area representing a hill. In FIG. 11, if there is a difference between the elevations of the inner contour line c ₂ and the outer contour line c ₁ , the detailed area Z ₁ is an area representing a hill. Among these segments, the relative elevations of the boundary contours are the same, in which case the area of the segment becomes a ridge or canyon. The distinction between these two cases cannot be judged by the boundary contours alone, but can be judged by the difference between the altitude of the contour line just outside the contour boundary and the elevation of the inner contour line just inside the contour boundary. For example, if the elevation of the contour just outside of the outer boundary contour is high, the detail area of the segment is a canyon, otherwise it is a ridge. And the internal contour c ₂ in Figure 11 another if there is no internal contours difference in elevation of c _3, subdivision Z ₂ is c _2, the height is c _2, c face is higher than the _third subdivision of c ₃ the outer contours c ₁ of contours Z ₂ is a ridge, otherwise a Cebu area Z ₂ is a canyon.

(2.3) If a segment belongs to a cycle linking cycles, the subdivision of that segment will be a ridge, canyon, or hillside for the same reasons as described in (2).

(2.4) In addition, if there are no inner boundary contours in the sub-domain Ω _i , the region of the subdomain may be peaks or pits as a whole. The distinction between these two topologies results in Ω _i being normal if the contour just outside the outer boundary of the subdomain sub-domain Ω _i is lower than the outer boundary, or otherwise the pool.

By applying these criteria, it is possible to determine that the subregion corresponding to all segments of the middle axis (MX) of sub-domain Ω _i corresponds to the topographical features of any one of the ridges, valleys, hills, canyons, summits, and ponds. have.

3. How to generate DEM data

(3.1) A method of generating DEM data from contour set data as source data using the terrain segmentation method described above will be described. The next task for generating the DEM is to calculate the heights for the surface points within each subdivision. The elevation of the desired ground points is calculated by the interpolation method of the elevations of two adjacent contours. As shown in FIG. 12, a conventional representative interpolation method for calculating the elevation of any point P between two contour lines 61, 62 is a vertical line, i.e., the shortest, from that point p to both contour lines 61, 62. Determined by the ratio of the distances d ₁ and d ₂ . This is represented by the following equation (1). In equation (1), h _p is the elevation at point P, h ₁ and h ₂ represent the elevations of two contours (61, 62), and d ₁ and d ₂ are the two contours (61,62) at point p, respectively. The shortest distance to) is shown, respectively.

......(One)

The problem with this shortest-range interpolation method, however, does not show consistency in elevation calculations. The elevation at point q in FIG. 12 is determined according to the ratio of the lengths of the two vertical lines, respectively, down to the two contour lines 61 and 62 at that point. Elevation has a property that varies somewhat continuously. Since the q point and the p point exist on almost the same slope line, interpolation using the same proportional expression based on the same slope line (path) can ensure the continuity of the elevation. Nevertheless, the existing interpolation method calculates the elevation by using different proportional equations when calculating the elevations of two points p and q, and thus the inconsistency is relatively high. Therefore, there is a lack as a model of the elevation calculation.

(3.2) The present invention proposes a new interpolation method that can calculate the elevation more precisely by using a medial gradient line model using a flux line model. Similar to the electric flux lines representing the electric field distribution between the high and low potential planes, the contours are viewed as equipotential lines and are directed from the high points (contour line c ₂ ) to the low points (contour line c ₁ ). It has been known previously that gradient line models can be created to describe the paths of the slopes of the statue. The slopes (Γ) have the properties of orthogonality, non-intersection, and monotonicity, as well as the nature of the electric flux in the electric field. Orthogonality is the property that slopes and contours are perpendicular to one another everywhere. The comparative difference is that the slope lines do not meet and do not prune each other. Monotonic variability is the property that elevation varies monotonically along the slope line between two neighboring contour lines with different elevations.

In interpolating the elevation of an arbitrary point between two contours, the present invention proposes a method of using a flux line between two contours, i.e., an intermediate gradient line (also called an "approximate gradient line"). That is, a method of calculating an approximate inclination line Γ passing through an arbitrary point p to obtain an elevation and interpolating using it. The approximate slope line through point p is considered to be the same as the flux line through point p. The elevation of point p can be obtained by finally calculating the flux line passing through that point and using it to establish a proportional equation for interpolation.

A method of obtaining a flux line will be described with reference to FIG. 13A. To find the flux lines between two consecutive contour lines (c ₁ , c ₂ ), first find the intermediate axis of the two contour lines (c ₁ , c ₂ ). The middle axis (MX) is the trace of the center points of all the maximum circles between the two contour lines (c ₁ , c ₂ ), as mentioned above. And that the draw up any source in the sense of c on the intermediate shaft (MX) are obtained two contours (c _1, c ₂₎ and the two contact points (p _1, p _2). If you connect two segments that fall from the center of the maximum circle c to the two contacts (p ₁ , p ₂ ) with two contours (c ₁ , c ₂ ), that is the flux line for the center point c.

In (a) of Figure 13, the two contour lines at any point p to obtain the elevation (c _1, c ₂₎ the closest point on either of the above p _i is the center c _i on the intermediate shaft (MX) A maximum circle 72 is defined and the circle 72 abuts the other contours at point p _{i + 1} . The flux line with the path p _{i + 1-} > c _i- > p _i passing through the point p can be seen as an approximate slope that reflects the continuously varying elevation of the points between two adjacent contours. . This flux line, that is, the approximate inclination line, is called an intermediate inclination line Γ _i , and this intermediate inclination line Γ _i is a reference line for establishing a proportional formula for elevation calculation in the present invention. The middle slope line Γ _i is defined by two line segments [p _i , c _i ] and [c _i , p _{i + 1} ]. Where c _i is the center of the circle on its intermediate axis MX, and p _i and p _{i + 1} are the contacts of the maximum circle 72 centered on c _i and the two contour lines c ₁ , c ₂ . Therefore, any intermediate slope (Γ _i ) satisfies the orthogonality, the non-difference, and the monotonic variability. That is, they reach the contour lines c ₁ and c ₂ at right angles everywhere (orthogonality). The middle slope line between two contours of different elevations does not intersect each other (non-crossing). This is because the trajectories of the two tangent points on the two contour lines c ₁ , c ₂ of the largest circle moving in one direction along the intermediate axis segment move in the same direction along the contour line. The height of a point between two contour lines (c ₁ , c ₂ ) is calculated to be monotonous along the arbitrary intermediate slope line (Γ _i ) (monovariance).

The intermediate axis MX is a continuous line and therefore there is an intermediate slope line Γ _i passing through it at any point within the two contour lines c ₁ , c ₂ . Therefore, by using the intermediate gradient line (Γ _i ) model, it is possible to obtain the elevation of all the points on the line (Γ _i ). The elevation calculation of any point located on the same middle slope line Γ _i is interpolated using the same proportional expression based on the middle slope line Γ _i . Therefore, unlike the conventional interpolation method described with reference to FIG. 12, an elevation calculation that reflects a characteristic in which elevation is continuously changed along the same flux line may be performed, thereby maintaining consistency of the elevation calculation.

(3.3) Once the median slope line (Γ _i ) for the point at which the elevation is to be obtained is obtained, use it to calculate the elevation at that point using a suitable known interpolation method. The most representative interpolation method available is a linear interpolation method represented by equation (1). However, this linear interpolation guarantees the monotonicity of the middle slope lines, but the slope of the middle slope line suddenly changes at both ends of the middle slope line, so that the overall slope line is not smooth. Monotone Hermite interpolation is known as an interpolation method that ensures the monotonicity of intermediate slope lines as the overall intermediate slope line becomes smooth. It is more preferable to use this interpolation method as an interpolation method for calculating the elevation of the points on the intermediate gradient line Γ. It is preferable to apply the forged hullite interpolation method when the elevations of two adjacent contour lines are different, such as hillsides, valleys, ridges, and the like. However, when the topographic features have the same height of two adjacent contours, such as a canyon, it is more preferable to apply convex Hermite interpolation, which is an interpolation method that well reflects the topographical features. Even if the feature of the terrain is a ridge, the convex hull is also confined if the heights of two adjacent contours defining the ridge are the same, or if the feature of the terrain is a convex (or concave) terrain such as a summit or a pond. It is more preferable to apply the mite interpolation method.

(3.3.1) FIG. 14 is a diagram for explaining a method of performing elevation calculation in a hillside region by the forged hullite interpolation method (FIG. 14 is shown assuming i = 1 for simplification of notation. ). In order to interpolate the elevation of a point in a sub-domain of a set of contours, what terrain features the sub-domain areas between them (exactly the subregions of the segment of the intermediate axis dividing the sub-domain areas) have? Assume that it is predetermined in the manner described above. Intermediate slope line (Γ) is started from a point p _{i + 1} on the contour c _{i + 1} and it is assumed that ends at a point on the contour p _i c _i. The elevation h of any point p on the middle slope line Γ can be calculated by interpolation along the middle slope line Γ. Let d _i and d _{i + 1} (indicated by d ₁ and d ₂ in the figure) are the distances from point p to p _i and p _{i + 1} , respectively, along the middle slope line Γ. The simplest kind of interpolation is linear interpolation defined as in Equation (2) above (see FIG. 13B). Formula (2) is substantially the same as Formula (1).

......(2)

However, this linear interpolation method causes a problem that the intermediate gradient line is not smooth at the point on the contour line because the first derivative is discontinuous in the contour line. To avoid this problem, a cube interpolation method can be used to create a smoother intermediate slope. The reason is that the middle slope provides the continuity of the first derivative in the contour. In other words, the middle slope line is smooth in the contour line. However, this method has the disadvantage of interpolating elevations outside the effective range [h _i , h _{i + 1} ] of the elevations on its intermediate slope line (Γ). Forged Hermite interpolation utilizing cubic Hermite spline or rational Hermite spline can be used to eliminate this drawback. In addition, a monotone rational Hermite spline is used to satisfy this property, since the height must be monotone for the point between two points with different elevations.

In FIG. 14, the elevation h of any point p in the two contours (c _i , c _{i + 1} ) is determined according to the monotone rational Hermite spline. In this case, the spline on the middle slope line Γ is as shown in Equation (3) below. Computing the elevation h of any point p between two contours belonging to a hilly area using Equation (3) always guarantees the monotonic variability of the middle slope line (already proven. Figure 14 (b) shows the middle slope line (Γ) ) Shows a monotonous change in a slightly curved downward direction compared to the straight line.

...... (3)

Here, Δi = (h _{i + 1} -h _i ) / (d _i + d _{i + 1} ) and t = d _i / (d _i + d _{i + 1} ). g _i and g _{i + 1} are the slopes (or derivatives) at two points (p ₁ , p ₂ ) where the intermediate slope (Γ) meets two contours (c ₁ , c ₂ ). Where i = 1, 2, ...

To calculate the elevation in hills using Equation (3), it is necessary to calculate the slope g _j at the point p _j where the intermediate slope line (Γ) intersects the contour line c _j (where j = 1, 2, ... ). 15 is a diagram for explaining this. Gradient g _j is the inner contour of the like, the point of the contour passing through the p _j c _j shown in Figure 15 c _j ⁺ and the external contour c _j ^- two intermediate slope line Γ ⁺ (inner contour at the point p _j c _j ^A linear slope 92 between two points p _j ⁺ and p _j ^{− that} meets Γ ⁻ (which is an intermediate slope extending from the point p _j to the outer contour line c _j ⁻ ), respectively. More specifically, if the contour c _j is non nor even the outermost contour innermost contour slope line Γ (Γ ⁺ and ^Γ-sum of) the two sub-domain Ω _j ⁺ and Ω _^j-extended to Can be. Where c _j is the outer boundary of sub-domain Ω _j ⁺ or the inner boundary of sub-domain Ω _j ⁻ . and h _j ⁺ Ω _j ⁺ the elevation of the internal contours of, h _j ^- a _j Ω ^- is the elevation of the outer contour. Further, _j ⁺ l is the sub-domain and Ω _j ⁺ a medium slope of the line length extending in the Γ ^+, l _^j-sub-a length of ^- the intermediate slope line extending in the ^Γ-domain Ω _j. The slope (derivative) at the point p _j g _j can be approximated by the following equation (4).

......(4)

At the boundary points on the outermost contour c ₁ or the innermost contour c _n , the slope (derivative) is approximated by the following equation (5).

...... (5)

Where h ₁ is the height of the outermost contour c ₁ and h ₁ ⁺ is the height of the inner contour of c ₁ . h _n is the height of the innermost contour c _n and h _n ^- is the height of the outer contour of c _n . l ₁ ⁺ and l _n ⁻ are also the lengths of the extended intermediate slope lines Γ ⁺ and Γ ^{− in} the sub-domain Ω ₁ ⁺ and in the sub-domain Ω _n ⁻ as described above.

(3.3.2) Next, FIG. 16 is a diagram for explaining an interpolation method relating to a summit and a pond. If a contour does not contain any other contours, that contour area corresponds to either a top or a pond. The topographical feature of the innermost contour c _n inner region in FIG. 16 compares the height of the contour with the contour c _n-1 immediately adjacent to it and is normal if the former height h _n is higher than the latter height h _n-1 and vice versa. If it is a pool.

When the innermost contour c _n inner region is normal, the inner region is divided into a plurality of detailed regions by using the intermediate axis MX and important points in the inner region as shown in FIG. For example, to find the elevation h _p of a point _p , first obtain the flux line passing through that point, that is, the approximate (or intermediate) slope line Γ _p , as described previously. In FIG. 16, when there is no other contour line inside the contour line c _n including the point p, it is a normal or puddle terrain. The approximate inclination line passing through the point p is in contact with the contour line c _n at another point p _2, where the maximum circle MC contacting the contact point p ₁ meets when the vertical line is lowered to the point on the contour line c _n closest to the point p. At that time, the contact point p ₁ and the maximum circle (MC) connected to the second line segment connecting the center point p ₃ and the contact point p ₂ of the first line segments and the maximum circle (MC) connected to the center point p ₃ lines from just Approximate (or intermediate) slope line Γ _p .

The middle slope line Γ _p passing through the point p in FIG. 16A may be approximately modeled as the convex line 110 in FIG. 16B. Find the middle slope Γ _p and then use the ratio of the distances d ₁ and d ₂ from the two points p ₁ and p ₂ where the contour line c _n meets the middle slope Γ _q to the point p at which the elevation is to be obtained. Can be obtained by interpolation. Various known interpolation methods such as the linear interpolation method and the forged hermite interpolation method mentioned above may be applied here. In particular, terrain such as a summit or a pond is preferably interpolated while ensuring convexity (convex lines 110 and 112 represent the normal convexities in Figs. 16 (b) and (c), respectively). . The known convex rational Hermite spline has the advantage of always ensuring convexity. Equation (6) below is an equation for calculating an elevation of a desired point p by interpolating elevations of two points p ₁ and p ₂ according to the convex glass hermite spline. The elevation calculation of point q can be done in the same way as above.

(6)

here,

Is as described in equation (3).

Even if the contour c _n inner region is a pit, the same method can be used for the elevation calculation because the pond is a symmetric topography feature. When the point p in FIG. 16 (a) is any point in the pool area, the middle slope line Γ _p passing through the point is roughly concave as modeled by the concave line 114 in FIG. 16 (d). Sex needs to be guaranteed. Therefore, it is advisable to use interpolation to ensure concavity in elevation calculations for pond areas. Equation (6) can be applied as it is an interpolation equation that always guarantees concaveness.

(3.3.3) Next, FIG. 17 is a diagram for explaining elevation calculations when the terrain feature is a ridge. The following conditions may apply to ridges:

i) When two adjacent contours (c ₁ , c ₂ ) are different contours (i.e. c ₁ ≠ c ₂ and h ₁ ≠ h ₂ ), as in FIG. 17 (b), gp ₁ > 0 and gp ₂ > If 0, the regions within both contours (c ₁ , c ₂ ) are ridged.

ii) as in Figure 17 (c), if the two contours (c _1, c ₂₎ have different internal contours (i.e., c _₁ ≠ c ₁ and h ₁ = h ₂ and gp _1> 0 while gp _2> 0, The regions within the two contours (c ₁ , c ₂ ) are also ridges.

iii) As shown in FIG. 19, if two contours (c ₁ , c ₂ ) are different internal contours (ie c ₁ ≠ c ₂ and h ₁ = h ₂ ) or two contours (c ₁ , c ₂ ) If this is contour (ie c ₁ = c ₂ and h ₁ = h ₂ ), and gp ₁ > 0 and gp ₂ > 0, then the regions within the two contours (c ₁ , c ₂ ) correspond to ridges.

Elevation calculations in ridge sub-domains are similarly performed by first obtaining the intermediate slope and then interpolating the elevations of the two contours using proportional expressions. Consider the case where the elevation h of the point p is obtained from (a) of FIG. 17. Two adjacent contours (c _1, c ₂₎ are not in contact in the first maximum circle (MC _1, and the other contour lines (c ₂₎ in contact with the contact point p ₁ to see when got off the vertical line to the nearest contour (c ₁₎ to the point p This is the case of obtaining an approximate slope line passing through the point p. To find an approximate (or intermediate) slope line passing through the point p, first find a detail region including the point p for which the elevation is to be obtained, and the detail We look for the intermediate shaft (MX) segment corresponding to the region. points p closest adjacent contour line c to find the contact p ₁ vertical line meets down to _1. Further, the contact point p1 in contact with the center point of the two adjacent to the contour line in (c _1, Find the center point p ₅ of the first maximum circle MC ₁ located on the intermediate axis of c ₂ ), and find the line (first line segment) from the contact point p ₁ to the center point p ₅ . p from the center point ₅ of the MC ₁₎ cycle (MX-1) of said intermediate shaft (MX) Toward the dangling subtree while following a segment of (MX-2), the first and the other contour lines _{(c. 2)} second dangling subtree to the center point p ₃ of the largest circle (MC ₂₎ in contact at the contact point p ₂ in Finally, the segments (second line segments) of are obtained, and finally, the line (third line segment) connecting the center point p ₃ of the second largest circle MC ₂ and the contact point p ₂ is obtained. The line connecting the second line segment and the third line segment is an approximate inclination line 116 to be obtained, that is, the path of the approximate inclination line 116 is' p _1- >p-> p _5- > p _4- > p _3- > p ₂ '.

The obtained middle slope line 116 may be modeled as a curve having a steeper slope from the high point p ₂ to the low point p ₁ as shown in FIG. 17B. The point p is located on the mid-slope line 116, the elevation h of the point p is the basis of the intermediate slope line 116, the distance from p ₁ to p the distance from d ₁ and p to p ₂ d Depending on the ratio of ₂ , the elevations h ₁ and h ₂ of the two contours can be obtained by interpolating using the proportional equations presented by the linear interpolation method or the monotonic Hulmite interpolation method described above.

(3.3.4) FIG. 18 is a diagram for explaining the elevation calculation when the terrain feature is a valley. A condition that may correspond to a valley is that two adjacent contours (c ₁ , c ₂ ) are different contours (that is, c ₁ ≠ c ₂ and h ₁ ≠ h ₂ ), as shown in (a) and (b) of FIG. 18. When gp ₁ <0 and gp ₂ <0, the regions within the two contours (c ₁ , c ₂ ) are valleys.

Since valleys and ridges are symmetrical topographical features, elevation calculations for them can be done in the same way. Naturally, a case where the contour lines of h ₂ c ₂ adjacent contour c is lower than the elevation h ₁ of the _first (i.e., h ₂ <h _1). Elevation calculation of valley area is same as calculation of elevation of ridge area. After the middle slope line 120 is obtained, the altitudes h ₁ and h ₂ of the two contours are based on the linear interpolation method or forging hermite interpolation method. Interpolation is done using proportional expressions. However, the characteristic of the middle slope line 120 in the valley region is different from that in which the slope is modeled as a gentler curve from the high point p ₂ to the low point p ₁ , as shown in FIG. 18B.

(3.3.5) Finally, as shown in Fig. 20, if two contours (c ₁ , c ₂ ) are different internal contours (ie c ₁ ≠ c ₂ and h ₁ = h ₂ ) or c ₁ and c _2, yet the same contour, gp ₁ <0 and gp _2> 0, the internal contours of the sub-domain is the canyon. Also, as shown in FIGS. 18A and 18C, the two contours (c ₁ , c ₂ ) are different internal contours (ie, c ₁ ≠ c ₂ and h ₁ = h ₂ ), gp ₁ <0 and gp _{If 2} > 0, in this case the area within the two contours (c ₁ , c ₂ ) corresponds to the canyon. Even in the case of the canyon, the elevation calculation is symmetrical with the ridge shown in FIG. 19, and thus, the elevation calculation may be performed in the same manner as described above.

Elevation calculations in the saddle area do not need special consideration. An eye point is a point where several valleys and ridges meet each other and is calculated in the same way as the elevation calculation of a point belonging to a valley or ridge.

(3.4) DEM data generation system

On the other hand, the DEM generation method described above may be implemented as an application program that can be executed in a general-purpose computer or a dedicated computing device. If you give the input data needed to create the DEM on the computer and run the application, the desired DEM result data will be automatically generated. 21 illustrates a configuration of an application program 230 implementing the DEM automatic generation method and a configuration of an apparatus for executing the same. Fig. 22 is a flow chart showing the flow of a procedure for generating DEM result data by this application program processing DEM input data. The hardware configuration of the DEM calculation apparatus 200 includes a storage device 220 storing input data required for generating a DEM, result data obtained through a DEM generation operation, and an application for automatic DEM generation, and reading an application program from the storage device. A CPU 210 for processing the DEM input data to generate and store the DEM result data and storing the DEM result data in the storage device 220, and a data processing space for the processing of the CPU 210. It provides a memory 250. If necessary, a screen output unit 260 is further provided to display the DEM result data provided by the CPU 210 on the screen. The application 230 is largely a data input unit 232, contour tree configuration unit 234, intermediate axis calculation unit 236, pseudo-dual intermediate axis decomposition calculation unit 238, geographic feature determination unit 240 of the detail region ), A height calculator 242, and a DEM output unit 244.

The data input unit 232 has a function of inputting data related to several contours and data for generating a DEM. When the CPU 210 executes the application program 230, the data input unit 232 reads input data necessary for generating the DEM stored in the storage device 220 (step S10). As the input data, the contour set data to be generated by the DEM and grid point information for generating the DEM are required. Contours are represented by approximate polygons. The contour set c ₁ , c ₂ , ..., c _n data includes the number n of total contours, the number and coordinate values of the total vertices of each contour, and the data of the heights of the contours. The DEM result data is expressed as height values for grid points that are parallel to the coordinate axis and spaced evenly. Therefore, the coordinate data of the grid point must be given for the DEM generation. Since the number of grid points is numerous, it is common and reasonable to provide reference point coordinates and calculation conditions for calculating grid point coordinates. FIG. 23 is a diagram for explaining a calculation method of grid point coordinates for a set of contour lines by way of example. When the grid point P ₁₁ located at the lower left as the reference point 160 is illustrated in this figure, the reference point ( Given the coordinates P ₁₁ (x ₁ , y ₁ ) at 160) and the spacing d between the grid points, and the number of grid points u in the x-axis direction, and the number v of grid points in the y-axis direction, The coordinate value of the point can be calculated. In practice, the spacing d between grid points will be a much smaller value than that shown in FIG.

The CPU reads the input data for DEM generation and uses it to perform the operation according to the following procedure. First, the contour tree constructing unit 234 is executed to generate a contour tree representing correspondences of the contour lines by using input data for the entire contour sets c ₁ , c ₂ , ..., c _n (step S12). . In other words, as described above, the trajectories and elevations of all the contours are used to calculate the inclusion relations between the entire contours, to determine the neighboring contours of each contour line, and to construct the neighbor relations in a tree form.

Once the contour tree is created, all sub-domains Ω _i (where i = 1, 2, ..., n) are calculated (step S14). Then, for each of the obtained sub-domains, each sub-domain Ω _i , in order from the first sub-domain Ω ₁ to the last sub-domain Ω _k , where i = 1, 2, ..., n ), Using the data of the sub-domain (step S16), the following operation is performed (steps S16-S26).

(a.1) the intermediate shaft calculation unit 236, all the sub-intermediate shaft is calculated in the domain Ω _i, respectively (Step S18). The intermediate axis calculator 236 considers the outer boundary contours that define the sub-domain Ω _i , which is the area between each adjacent contours in the contour tree, as a polygon, and considers the inner boundary contours as the holes of the polygon, Has the function of calculating the intermediate axis of the polygon.

(a.2) After the intermediate axis is obtained, the pseudo-dual intermediate axis decomposition calculation unit 238 divides the calculated intermediate axis for each sub-domain by segment, calculates subregions corresponding to all segments, and divides the domain. Decompose into several zones (step S20).

(a.3) Then, for each subregion obtained through the decomposition of each sub-domain, the geographic feature determination unit 240 of the subregion is the topographical feature of the peaks and the pits. It is determined whether the hills (hillsides), ridges (valleys) or valleys (canyons) (step S22).

After determining the detail region and its topographical features for all sub-domains, the height calculator 242 performs elevation calculations for all grid points for which DEM data will be calculated. The height calculator 242 first calculates coordinates of the positions of all grid points for calculating the height by using the input data for generating the DEM. For every grid point, the coordinates are used to determine the subregion of the sub-domain to which each grid point belongs. The elevation of the grid point is calculated by applying an appropriate interpolation method according to the determined terrain feature (step S28).

A calculation procedure for this purpose will be described with reference to FIG. 23. Among all the grid points, the calculation procedure is performed on the vertical line of the i-th (where i = 1, 2,. The elevation of all grid points is calculated by performing the following process on the grid points.

(b.1) Calculate the x-coordinate of the i-th vertical line using the equation x _i = x ₁ + d × (i-1). However, the reference point coordinate p ₁₁ (x ₁ , y ₁ grid point spacing d is given as input data.

(b.2) From the bottom of the lattice points lying on the i-th vertical line, moving from the bottom to the y coordinate value of the lattice point of the j-th (where j = 1, 2, ...., v) _{Find j} = y ₁ + d × (j-1).

(b.3) When the coordinates (x _i , y _j ) of a specific grid point p _ij are obtained as described above, a calculation is performed using the contour tree to find the sub-domain Ω _ij to which the grid point p _ij belongs. Further, a calculation is performed to find out which zone in the sub-domain Ω _ij belongs to the grid point p _ij .

(b.4) When the subzone to which the lattice point p _ij belongs is found, a calculation is made to determine whether the terrain features of the subdivision are ridges, valleys, hills, canyons, ponds, or summits. .

(b.5) Calculate the height of the grid point p _ij by interpolating the elevations of two adjacent contour lines, using the interpolation method described above, according to the determined topographical features of the subregion.

As such, in a loop (b.1) that increases the value of i by 1 from 1 to u, a loop that increases the value of j by 1 by 1, that is, includes (b.2) through (b.5) By performing a double loop, elevation can be calculated according to all coordinates of the entire grid point array and the topographical features of the detailed area where each grid point is located.

When the elevation calculation is completed, the DEM output unit 244 outputs the coordinates of all grid points for generating the DEM and the elevation data of the grid points to the storage device 220 to be stored. If necessary, the calculated coordinates and elevation data (DEM result data) are also output to the screen output unit 260 to be displayed on the screen. Shown in (b) of FIG. 24 shows a screen in which DEM data is generated and rendered according to the method described above using the contour set of FIG. 24 (a) as source data.

4. Stereoscopic surface reconstruction method of terrain

Using the terrain segmentation method described above, the surface of the terrain can be reconstructed more realistically. To this end, it is necessary to create a triangular mesh for surface reconstruction using subregions obtained by decomposing sub-domains between contours. To make a triangular mesh, we first find the intermediate axis for the sub-domain as described above, and then get the subregions corresponding to all segments of this intermediate axis. This intermediate axis serves as a dividing line connecting the critical points that represent the features of the terrain. Some segments of the intermediate axis are represented by parabolas, which straighten those parabolic segments. That is, the parabolic segment is transformed into a straight line connecting the main points at both ends thereof. FIG. 25 is an example of straightening all segments of the intermediate axis illustrated in FIG. 4 in this manner. Sub-domains are divided into triangles or squares by the boundaries of all straight segments and all subregions. If necessary (ie, if the rectangle's vertices are not on a plane, etc.), all rectangles should be subdivided into two triangles. This is because triangles are better suited to representing planes than rectangles. FIG. 26 shows an example in which all the rectangles of FIG. 25A are divided into two triangles.

The sub-domain is decomposed into triangular or rectangular subregions, and then three-dimensional coordinates (ie, two-dimensional plane coordinates and height values) of the vertices of the divided triangular or rectangular subregions are calculated. The plane coordinates of the vertices of the detail region can be easily obtained in the decomposition process into the detail region. The elevations of the vertices of the detail region can be obtained directly from the data of the corresponding contours for the vertices located on two

adjacent contours

42 and 44, while the heights of the vertices not located on the two

contours

42 and 44 are In the above-described DEM data generation, the height of the grid points is calculated in the same manner as that of the grid points, that is, the elevation of two

adjacent contour lines

42 and 44 is interpolated by a suitable interpolation method.

The triangular (or triangular or rectangular) mesh data created in this way, that is, the three-dimensional coordinates of the detail region, is the surface reconstruction data of the terrain to be obtained immediately. Since this method is based on the terrain segmentation method, which can express the topographical features very closely to the real world, the representation reconstruction data obtained by this method expresses the topographical features of the sub-domain more realistically than the data obtained by the conventional method. Done.

Rendering with the three-dimensional coordinates of the obtained detail region makes it possible to visualize the surface of the detail region. When you connect the rendered surfaces of the entire detail area, that is the three-dimensional surface of the terrain that you want to represent as contours.

The three-dimensional surface reconstruction of the terrain can also be made using the DEM result data described earlier. The DEM result data consists of coordinates and elevation data for each grid point. The entire grid point of the DEM result data can be decomposed into a mesh of polygons such as a plurality of triangles or squares by connecting three or four adjacent grid points. The plane coordinates and elevations of the vertices of each polygon can be known from the DEM result data. The three-dimensional vertex coordinate values of the polygon meshes thus obtained form desired surface reconstruction data, and when the data is rendered, the three-dimensional surface of the terrain is obtained.

In describing the preferred embodiment of the present invention, several new and important ideas have been proposed: (Idea 1) finding the intermediate axis of the sub-domain of contours; (Idea 2) dividing the intermediate axis into segments based on the significant points and decomposing the sub-domain into detailed regions corresponding to each segment; (Idea 3) matching topographic features by decomposed subregions; (Idea 4) Find the median slope for adjacent contours; (Idea 5) selecting appropriate interpolation methods by terrain feature; (Idea 6) Calculate the elevation of the desired point using the middle slope line. Those skilled in the art will be able to derive various variations of the present invention by combining these ideas as appropriate. For example, it is possible to calculate the elevation of a desired point in a region represented by a set of contour lines. Elevation calculations at specific points may be performed by simply combining

ideas

1, 4, and 6, or alternatively, using ideas 1 through 6 in their entirety.

In addition, through the above description, those skilled in the art will understand the methods of the present invention, that is, a method of decomposing a terrain represented by a contour set into detailed regions, a method of calculating an elevation of a desired point using the contour set, and a contour set. The method of generating DEM by subdivision decomposition method and the method of reconstructing three-dimensional surface of contour display terrain can be made into application program that can be executed by general purpose personal computer or computing device such as microcomputer, microprocessor, DSP, etc. . Such application programs may be recorded and provided in the data storage means as well as data storage means that are known in the art such as a hard disk, a magnetic tape, a CD, a flash memory, and the like. A computing device such as a computer can implement the methods of the present invention by reading and executing the program from a recording medium on which such a program relating to the method of the present invention is recorded.

The height calculation technique of the specific position of the present invention can be used to calculate the approximate height data for the specific position more accurately using the contour line without actually going to the field and using a technique such as triangulation or Global Positioning System (GPS). have.

DEM automatic generation technology of the present invention is used to orthodontic high resolution satellite image or aerial image, and is also used in various fields such as terrain visualization, three-dimensional map, model production. It is also used to represent terrain more realistically in applications such as geographic information system (GIS), flight simulation, military training and operational simulation, 3D computer games, virtual reality, and more.

The present invention described above is merely a preferred embodiment of the present invention, and those skilled in the art can vary the present invention without departing from the spirit and scope of the present invention as set forth in the claims below. Can be modified and changed. Accordingly, all changes that come within the meaning or range of equivalency of the claims are to be embraced within their scope.

Claims

A computer-implemented method that decomposes the terrain represented by a set of contours into detailed regions

Obtaining an intermediate axis for each sub-domain which is an area between adjacent contour lines using the contour set; And

Decomposing each sub-domain into zones for all segments of its intermediate axis,

The segment of the intermediate axis is a set of contours defined as the set of all the center points of the maximum circles that contact the same two active boundary elements, and the active boundary elements are defined by the corners and concave vertices of the sub-domain that abut the maximum circle. To decompose the terrain represented by subdivisions.
The method of claim 1, wherein the obtaining of the intermediate axis comprises:

Creating a contour tree representing correspondences between the contour lines from data relating to the contour set;

Calculating all sub-domains that are regions between adjacent contour lines using the contour tree; And

Calculating the intermediate axis for each sub-domain, wherein the terrain to be represented by the contour set is decomposed into detailed regions.
The method of claim 1, wherein the intermediate axis is a trace of the center points of the maximum circles inscribed in the sub-domain at two or more points, and the relationship between the maximum circles is such that any maximum circle is a complete subset of another maximum circle in the polygon. A method of decomposing a terrain to be represented by a set of contours, which is characterized by having a non-relational relationship into detailed regions.
2. The method of claim 1, wherein the intermediate axis is a simple tree graph composed of vertices of the polygon and parabolic or straight edges connecting the vertices when the sub-domain is a simple polygon. In the case where the domain is a polygon including one or more holes therein, the contour line further includes cycles surrounding each hole in addition to the simple tree, and in some cases, the graph further includes a path connecting each cycle. How to decompose the terrain you want the set to represent into detail regions.
The boundary line of claim 1, wherein the boundary line of the subregion is a center point of two maximum circles positioned at each of a critical point at which a new segment starts and an important point at which each new segment starts while tracking the intermediate axis. Radial segment of the largest circle connecting the contact with the active boundary element,

The key points include, among the points on the intermediate axis, end points that are points of the intermediate axis that intersect the corner boundaries of the polygon, junction points that are the center points of the largest circle in contact with three or more active boundary elements, and two points. It is a center point of the maximum circle contacting the active boundary element, and there are three kinds of transition points, which are the center points where only one of the two active boundary elements changes. Way.
The method of claim 5, wherein the obtaining of the intermediate axis,

Calculating an intermediate axis by tracing the center point of the maximum circle in which the size of the circle continuously changes and contacts the active boundary element starting from one important point;

When tracking the center point, finding one segment that connects the key point that started by calculating the next key point and the next point that is calculated; And

If the newly calculated next critical point is the endpoint, ending the tracking of the center point of the largest circle; decomposing the terrain to be represented by the contour set into detailed regions.
7. The segment of claim 6, wherein the segment's endpoint is the point where the segment's extension meets i) a line extending perpendicular to the polygonal edge at the concave vertex of adjacent active boundary elements, or ii) at the convex vertex of the adjacent active boundary elements. A method of decomposing a terrain to be represented by a set of contour lines, which is defined as the point where the vertex angle meets an extension line that bisects the vertex angle first.
The method according to claim 5, wherein the center point of the two largest circles located at each of the starting point and the ending point of the new segment each time while tracking the intermediate axis, and a contact between the two active boundary elements which each of the two maximum circles are in contact with. The radius of the maximum circle connecting becomes the boundary of the subregion, and the area closed by the subregion boundary and two active boundary elements is the subregion corresponding to the segment. How to decompose it into regions.
A computer-generated method of calculating elevations at desired points using a set of contours,

Calculating an approximate slope line passing through the point for which the elevation is to be obtained; And

Calculating the elevation of the point as the elevation of the point by interpolating the distance between two adjacent contours of the point along the approximate inclination line and the elevation of the two contours adjacent to the point using a predetermined interpolation method. Elevation calculation method comprising the.
10. The method of claim 9, wherein the approximate inclination line passing through the point is the maximum circle in contact with the contact point p 1 which meets when the vertical line is lowered from the two adjacent contour lines (c 1 , c 2 ) to the contour line (c 1 ) close to the point. When contacting the other contour (c 2 ) at point p 2 , the first line segment connecting the contact point p 1 and the center point of the largest circle and the second line segment connecting the center point of the maximum circle and the contact point p 2 Elevation calculation method, characterized in that determined by the connecting line.
10. The method of claim 9, wherein the approximate slope line passing through the point when there is no other contour inside the contour including the point is the maximum contacting the contact point p 1 encountered when the vertical line falls to the point on the contour that is closest to the point. a source is connected to the center point and a second line segment connecting the contact point p 2 of another point p 2 when in the access and the contour lines, the contact points p 1 and the maximum source a first line segment and the maximum source connected to the center point of the line Elevation calculation method, characterized in that determined by.
The first maximum circle (MC 1 ) of the two adjacent contour lines (c 1 , c 2 ) contacting the contact point p 1 when the vertical line is lowered to the contour line (c 1 ) close to the point is the other contour line. When not in contact with (c 2 ), the approximate inclination line passing through the point is i) the first maximum located on the intermediate axis MX of the contact p 1 and the two adjacent contour lines c 1 , c 2 . A first line connecting the center point p 5 of the circle MC 1 , ii) from the center point p 5 of the first largest circle MC 1 along the segment of the dangling subtree in the cycle direction of the intermediate axis MX going, first the second line segments and, iii consisting of segments of dangling subtree to the center point p 3 of the second largest circle (MC 2) in contact at the contact point p 2 to the rest of the contour (c 2)) last to connect the third line segment connecting the contact point p 2 and p 3 and the second center point of the largest circle (MC 2) It is determined by,

The intermediate axis MX is a set of all the center points of the maximum circles abutting the same two active boundary elements of the two adjacent contour lines c 1 , c 2 , and the active boundary element is the two adjacent contour lines (abutting the maximum circle). c 1 , c 2 ), defined as the corners of the sub-domain and concave vertices, wherein the cycle of the intermediate axis MX is a subset of the intermediate axis edges having a closed shape, and the dangling subtree is Elevation calculation method, characterized in that the maximum connection intermediate axis edge subset that does not belong to any cycle of the intermediate axis (MX) or do not belong to the path connecting the cycle.
10. The method of claim 9, further comprising obtaining a topographical feature of the detailed area to which the point belongs,

And the forged glass Hermite spline interpolation method is applied to the elevation interpolation when the terrain feature corresponds to a hilly, valley, or ridge and the elevation of two adjacent contours is different.
10. The method of claim 9, further comprising obtaining a topographical feature of the detailed area to which the point belongs,

If the feature of the terrain is the same as the canyon, or if the heights of two adjacent contours are equal, or if the height of the two adjacent contours defining the ridge is the same, or the feature of the terrain may be represented as convex (or concave). In the case of terrain such as a summit or a pond, an elevation calculation method characterized by applying a convex hullite spline interpolation method that reflects such terrain features well during the elevation interpolation.
A computer-generated method of generating a DEM using a set of contours,

Obtaining an intermediate axis for each sub-domain which is an area between adjacent contour lines using the contour set;

Decomposing each sub-domain into zones for all segments of its intermediate axis;

Calculating a topographic feature of each detail area; And

Calculating the elevation of each grid point of the DEM by interpolating the elevations of two adjacent contour lines using a predetermined interpolation method according to the terrain feature of the detailed area to which the grid point belongs.

The segment of the intermediate axis is a set of all the center points of the largest circles in contact with the same two active boundary elements, wherein the active boundary elements are defined by the corners and concave vertices of the sub-domain that abut the maximum circle. Way.
The method of claim 15, wherein the calculating of the feature of the terrain in the subregions examines the structure of the intermediate axis.

(1) If a segment belongs to a dangling subtree, the detail region of that segment is determined as part of one of the ridges, valleys, and canyons, according to the relative elevation of adjacent boundary contours,

(2) If a segment belongs to a cycle, the subdivision of that segment is determined by either hilly, ridge, or canyon,

(3) if a segment belongs to a route, the subdivision of that segment is determined by one of ridges, canyons, and hills,

(4) if the sub-domain does not have internal boundary contours, the area of that sub-domain is determined as a normal or a pool as a whole;

The cycle is a subset of intermediate axis edges having a closed shape, the dangling subtree is a maximum connected intermediate axis edge subset that does not belong to any cycle of the intermediate axis or does not belong to a path connecting the cycles, and the path is a cycle of cycles DEM generation method characterized in that the subset of the connecting intermediate axis edges.
The method of claim 15, wherein the step of calculating the elevation of the grid point,

Calculating an approximate gradient line passing through the grid point for which the elevation is to be obtained; And interpolating the distance ratio between two adjacent contour lines (c 1 , c 2 ) along the approximate inclination line with respect to the corresponding grid point and the elevation of the two adjacent contour lines (c 1 , c 2 ) using the interpolation method. DEM generation method comprising the step of.
The approximate inclination line passing through the grid point is the maximum contacting the contact point p 1 which meets when the vertical line is lowered from the two adjacent contour lines (c 1 , c 2 ) to the contour line (c 1 ) close to the grid point. a circle when encountered in point p 2 to the other one of the contours (c 2), connecting the contact points p 1 and the center point of the largest circle which the first line segment and the maximum circle connecting the central point with the said contact point p 2 2 DEM generation method characterized in that determined by the line connecting the line segment.
Of claim 17 wherein, the contact in the case to the inside of the contour including the grid points do not have other contours meet when got off a vertical line to a point approximate the slope line is on the nearest contour in the lattice points that passes through the grid point p 1 to the When the maximum circle in contact with the contour line is in contact with the contour line at another point p 2 , the first line segment connecting the contact point p 1 and the center point of the maximum circle and the second line segment connecting the center point of the maximum circle and the contact point p 2 DEM generation method characterized in that determined by the connected line.
18. The method of claim 17, wherein the first largest circle MC 1 in contact with the contact point p 1 which meets when the vertical line is lowered from the two adjacent contour lines c 1 and c 2 to the contour line c 1 close to the lattice point is different. When not in contact with the contour line c 2 , the approximate gradient line passing through the grid point is i) the second point located on the intermediate axis MX of the contact p 1 and the two adjacent contour lines c 1 , c 2 . 1 first line connecting the center point p 5 of the largest circle MC 1 , ii) a segment of the dangling subtree in the cycle direction of the intermediate axis MX from the center point p 5 of the first largest circle MC 1 Second line segments consisting of segments of the dangling subtree up to the center point p 3 of the second largest circle MC 2 first contacting the other other contour line c 2 at a contact p 2 , and iii ) finally opening the third line segment connecting the contact point p 2 and p 3 and the second center point of the largest circle (MC 2) Is determined by lines,

The cycle of the intermediate axis MX is a subset of the intermediate axis edges having a closed shape, and the dangling subtree does not belong to any cycle of the intermediate axis MX or to a path connecting the cycles. And a subset of edges.
16. The method of claim 15, wherein if the topographical feature of the detailed region to which the lattice point belongs is a hill, a valley, or a ridge, and the elevations of two adjacent contours are different from each other, applying the forged glass Hermite spline interpolation method during the elevation interpolation. Characterized by the DEM generation method.
16. The terrain feature according to claim 15, wherein the topographic feature of the subregion to which the grid point belongs is the same as the canyon, or the height of the two adjacent contours defining the ridge is the same, or even if the height is equal to the ridge. In the case of terrain, such as a summit or a pond, which can be represented by this convexity (or concaveness), DEM generation, characterized in that the elevation interpolation is applied to the convex hullite spline interpolation method that reflects such terrain features well Way.
A computer-implemented method of reconstructing the three-dimensional surface of the terrain represented by a set of contours,

Obtaining an intermediate axis for each sub-domain which is an area between adjacent contour lines using the contour set;

Decomposing each sub-domain into zones for all segments of its intermediate axis;

If a segment of the intermediate axis is represented by a parabola, it straightens by connecting the main points at both ends of the parabolic segment in a straight line; And

Calculating three-dimensional coordinates of the exploded detail polygons, that is, two-dimensional plane coordinates and elevation,

The segment of the intermediate axis is a set of all the center points of the maximum circles that are in contact with the same two active boundary elements, wherein the active boundary elements are defined by the corners and concave vertices of the sub-domain that abut the maximum circle. 3D surface reconstruction method of terrain.
24. The method of claim 23, further comprising: subdividing all of the decomposed subregions into two triangles.
24. The method of claim 23, wherein the elevations of the vertices of the subregions can be obtained directly from the elevation data of the corresponding contours for vertices located on two adjacent contours, and for vertices not located on the two contours. A three-dimensional surface reconstruction method of contour display terrain, which is obtained by interpolating elevations of two adjacent contours.
The three-dimensional surface of the contour display terrain according to claim 23, wherein the elevation calculation step calculates by interpolating the elevations of two adjacent contour lines using a predetermined interpolation method according to the terrain feature of the detailed area to which the elevation calculation point belongs. Reconstruction method.
24. The method of claim 23, wherein the elevation calculation step comprises: calculating an approximate slope line passing through the elevation calculation point; And calculating the elevation of the point as the elevation of the point by interpolating the distance ratio between two adjacent contours of the point along the approximate inclination line and the elevation of the two contours adjacent to the point by using a predetermined interpolation method. Three-dimensional surface reconstruction method of the contour display terrain comprising a.
24. The method of claim 23, further comprising: rendering with three-dimensional coordinates of all detail regions; And connecting the rendered surfaces of the entire subregions to represent the three-dimensional surface of the terrain.
A computer-implemented method of reconstructing the three-dimensional surface of the terrain represented by a set of contours,

Obtaining an intermediate axis for each sub-domain which is an area between adjacent contour lines using the contour set;

Decomposing each sub-domain into zones for all segments of its intermediate axis;

Calculating a topographic feature of each detail area;

Calculating the elevation of each grid point of the DEM by interpolating the elevations of two adjacent contour lines using a predetermined interpolation method according to the terrain feature of the detailed area to which the grid point belongs;

Decomposing the entire grid point of the DEM result data into a mesh of polygons such as a plurality of triangles or squares by connecting three or four adjacent grid points; And

Obtaining and assigning planar coordinates and elevations, which are three-dimensional coordinate values of vertices of each polygon, from the DEM result data,

The segment of the intermediate axis is a set of all the center points of the maximum circles that are in contact with the same two active boundary elements, wherein the active boundary elements are defined by the corners and concave vertices of the sub-domain that abut the maximum circle. 3D surface reconstruction method of terrain.
30. The method of claim 29, further comprising: rendering using three-dimensional vertex coordinate values of the polygon meshes; And connecting the rendered surfaces of the entire polygon meshes to represent the three-dimensional surface of the terrain.
A computer-readable recording medium having a computer recorded thereon a program for executing each step of the method according to any one of claims 1 to 30.