Subscribe Now Subscribe Today
Research Article

Contour Generation and Filling Algorithm of Irregular Region Based on Rectangular Grid Method

Chi Chen, Mei Li, Jinghe Zhao, Tianshu Yi and Dongqi Fan
Facebook Twitter Digg Reddit Linkedin StumbleUpon E-mail

An approach for contour generation and filling of irregular regions based on rectangular grid method is proposed. This approach entails contour generation in a rectangular area containing the irregular region of interest using rectangular grid method, followed by removal of the equivalent points outside the irregular region. An improved algorithm of ray method to determine the spatial relation between the equivalent points and irregular region is also reported. The above method has been successfully used in a geological disaster monitoring platform for drawing the isopleth maps of rainfall and soil moisture content in various regions.

Related Articles in ASCI
Similar Articles in this Journal
Search in Google Scholar
View Citation
Report Citation

  How to cite this article:

Chi Chen, Mei Li, Jinghe Zhao, Tianshu Yi and Dongqi Fan, 2014. Contour Generation and Filling Algorithm of Irregular Region Based on Rectangular Grid Method. Journal of Applied Sciences, 14: 368-373.

DOI: 10.3923/jas.2014.368.373

Received: September 25, 2013; Accepted: November 20, 2013; Published: February 08, 2014


Contours can be drawn with either rectangular grid or triangular mesh method: Rectangular grid method is very effective in drawing contours of a rectangular area, but most practical applications deal with irregular areas, such as the geological fault, contour map and isothermal chart of area with irregular boundaries. Triangular mesh method is often used instead in these applications to draw the contours of irregular region (Fu et al., 2006), which allows faithful representation of original data with a relatively simple program but lack capability for extrapolation. Wang (1991) proposed a rectangular grid method searching directly inside the irregular region, which however required large amount of computation. In this study, an improved method for contour generation for irregular regions based on established rectangular grid method is devised which has the characteristics of less calculation. At first, a rectangular area is created which contains the irregular region of interest; contours are then drawn for the rectangular area by use of the rectangular grid. The contours outside the irregular region are then erased, with the end results being the desired contours for the irregular area.

When removing the contours outside the irregular region, the ray method is used to determine the positional relationship between points and irregular region (Chen et al., 2007). In this study, an improved algorithm of the ray method with higher efficiency is also reported.


Rectangular grid method : The rectangular grid method is composed of three sequential steps, i.e., discrete data gridding, the contours generation and filling.

Discrete data gridding is to constitute a rectangle by the minimum and maximum values of the original data points horizontal and vertical coordinates; then, according to the requirement of the curve smoothness, the rectangle could be divided into m*n small rectangles, which constitute the rectangular grids; finally, rectangular grid point values can be estimated by the raw data using interpolation algorithm. Commonly used interpolation methods include inverse distance weighting interpolation method, Kriging interpolation method and natural neighbor interpolation method, etc.

The gridding of discrete data is followed by contour generation. This is accomplished by obtaining the coordinates of the equivalent points of each contour across the grid edge by tracking each contour grid-by-grid starting from the boundaries via linear interpolation method. The results (the new coordinates of points) are then stored in an array or a list and these points are eventually connected to form a continuous smooth contour.

The filling of the contour map entails fill the area between two consecutive contours with different colors, which allows better visualization of the horizontal and vertical changes of physical quantitative value.

Contour generation is the most important step for plotting of a contour map. This involves the calculation of equivalent points, contour tracking and contour drawing (Tang et al., 2008).

Calculation of equivalent points: The calculation of equivalent points is to find the coordinates of intersection between contours and grid boundaries. It should first be decided whether there is an intersection based on the coordinates of grid points. Then the coordinates of intersections are calculated by the linear interpolation algorithm.

Contour tracking: After all equivalent points are found, the contours are formed by connecting the points one by one in order by certain rules (Wang, 2008). This contour tracking process includes the following steps:

Determination of the starting point: The contours are divided into closed and open curves. Both the head and tail of an open curve are located on the boundary, therefore search for the starting point should start from the outer boundary; if no starting point is found, search the equivalent points from left to right and from the bottom to the top for the starting point. As for closed curve, any equivalent point can serve as the head or the tail
Equivalent point tracking: After the starting point of the contour is determined, track the equivalent points grid-by-grid, which are then connected sequentially to form the contour. The trends of contours entering or exiting the grid need to be considered when tracking equivalent points grid-by-grid. The contours have four possible trends to go into the grid: from bottom to top, from top to bottom, from left to right or from right to left. After a contour enters the grid, it exits from the other three sides of the grid. When a contour leaves one grid to track the next equivalent point, its direction and distance from the current equivalent point should be consider comprehensively to avoid the intersection of the contours and other unexpected outcomes
Determination of the end point: The starting point and the end point of an open curve are on the boundary of the rectangular area. As for closed curve, the end point is the same as the starting point. When tracking contours, the starting point of the contour should be labeled so that the contour could be terminated when the next intersection point becomes the starting point
Other problems: Confusion may arise when one equivalent point happens to be a grid point, as the latter is also the intersection of four adjacent grids. This could be resolved by adding a small value to the equivalent point

Contour generation of irregular region: A regular rectangular region that contains the irregular region can always be found. Rectangular grid method will allow us to generate contours of this rectangular region and then remove the parts outside irregular region. The key steps of irregular region contour generation include regular region contour generation, contour processing outside the irregular region and intersections processing of non-closed contours with the boundaries.

Regular region contour generation: In order to reduce calculation, the selected regular region is the smallest rectangle that contains the irregular region. It is defined by maximal and minimal of the abscissa and ordinate values of the irregular region.

The contours in the regular region can be obtained by use of rectangular grid method, as shown in Fig. 1.

Contour processing out of the irregular region: Contour processing outside the irregular region allows us to obtain the contours of irregular region. This is achieved by first determination of the positional relationship between equivalent points and the irregular region using the ray method, followed by elimination of all of the equivalent points outside the irregular region.

Fig. 1: Contour generation in the regular region that contains the irregular region (The black line indicates the boundary of the irregular region, red dots the discrete data, the blue line the smallest rectangular region that contains the irregular region and the gray curves the contours generated)

Fig. 2: Positional relationship between point and irregular region of ray method (~: Determination points, ~: Rays from the determination points, b1~b12: Polygon vertices)

The method of determining the positional relationship between points and irregular region by ray method is first to draw a ray from the point in any direction, then calculating the number of the intersections between the ray and irregular borders. If the intersection number is even, it can be deduced that the point is outside the polygons. If the intersection number is odd, then the point is within the polygon. As shown in Fig. 2, - are the 5 points to be determined, which serves as the starting points of 5 downward rays. These rays intersect with the polygon in different ways: Ray - intersects with the sides of the polygon; ray - and ray - intersect with the vertices of the polygon; ray - and ray - coincide with some of the polygon sides. In different cases, the calculation methods of intersection number are different.

For situations which the ray intersects with the sides of the polygon such as ray -, the positional relationship can be directly determined according to above rules.

When the rays coincide with a side of the polygon, as are the cases of ray - and ray -, there will be numerous intersections. The intersection number counts as two if the abscissa value of the point equals both abscissae values of two endpoints of the side. Otherwise, the intersection number counts as zero. Then the positional relationship of the point and the irregular region can then be determined according to above rules by the following program:

Int flag = 0; // The number of the intersections of the ray and irregular borders
For(Travel all boundary points){
If (the point to be determined is on the top of line section of two adjacent boundary points){
Flag ++ ; // then there is a intersection of the ray and the line section

Fig. 3: Intersections processing of non-closed contours with the boundaries (b1~b12 are polygon vertices, point are the equivalent points of a contour line)

The afore mentioned algorithm applies when the ray intersects with the polygon vertices just like ray and ray .

The existing algorithms in literatures require at least three or more conditions to determine the relationship between the point and the irregular region. The proposed algorithm in this study needs only two, because when it compares the coordinates of the point with that of the coordinates of two end points, it also compares the coordinates of two endpoints at the same time. This has greatly simplified the determination of the positional relationship when the ray intersects with the polygon vertices or coincides with irregular boundaries.

Intersections processing of non-closed contours with the boundaries: Once point has been removed, the contour line does not intersect with side b7, b8, which affects the contour filling. Therefore, it is necessary to find the intersection of line and side b7, b8. As shown in Fig. 3, line intersect with all unparalleled sides, so the algorithm needs to guarantee that the intersection is point which is on line , rather than point or other points.

Contour map filling algorithm: The contour filling can be divided into two steps: construction and filling of the equivalence region.

Construction of the equivalence region: Equivalence regions are classified into two categories: closed region which is formed by a closed contour and non-closed region which is formed by non-closed contours and boundaries. They are constructed differently. A closed region can be easily constructed because it is formed by a closed contour. When a closed region is included in another closed region, they will need to be differentiated. This could be easily achieved by finding the minimal abscissa value of each closed contour and sorting them in the order from low to high. This way the inner region is always behind the outer closed region.

For a non-closed region, since the start and end points of non-closed contours are on the boundaries, all the start and end points and the region vertices constitute a collection of boundary points. These boundary points are sorted and numbered clockwise or counterclockwise, as illustrated in Fig. 4 and 5. For convenience a regular region is chosen to explain the processes.

An equivalent area is constructed for each boundary point as a starting point and the number of times each point is involved in the structure of equivalent area is counted. For example, the construction process for the equivalent area from point 1 is as follows: Search from starting point 1 to point 2; since point 2 is the intersection of the contour lines and the boundary, point 6 can be found which is on the same contour from boundary point set; in this fashion, search from 6 to 7, 7 to 10, 10 to 11, 11 to 18 and eventually back to 1. The search process ends after tracking back to the starting point and an equivalent area is formed.

This process can be repeated to construct all equivalent areas for each boundary point. If the boundary point is a boundary vertex, the process is skipped when the number for the boundary point to participate in the structure of equivalent area is 1. If the boundary point is an intersection of a contour and rectangular boundary, the process is skipped too when the number is 2. That is because a non-closed contour has two intersections with rectangular boundary, each intersection point will participate in the construction of non-closed area twice, while a boundary vertex only participates in equivalent area construction once.

Fig. 4: Sorting boundary points of non-closed contours counterclockwise

Fig. 5: Numbering boundary points of non-closed contours after sorting

Equivalent area filling: Each equivalent area needs to be assigned a value according to surrounded contours after all equivalent areas are formed. Different values correspond to a set of gradual color. These areas are then filled with corresponding colors (Zheng et al., 2010).

Non-closed areas do not overlap, but may contain closed areas. Closed areas may also contain closed areas. Therefore the non-closed areas are filled first, followed by the closed ones. Similarly outer closed areas are filled before inner closed ones and big areas before smaller ones.

For any equivalent area, the surrounding contours have two elevation values at most. If the surrounding contour has only one elevation value, then it will be assigned this value. If its surrounding contours have two elevation values, then the average value will be assigned.

Fig. 6(a-d): Contour generation and filling of irregular region (a) The rainfall data used as discrete (b) Contour drawn by standard data for drawing contour (the red rectangular grid method in the dots indicate the locations of the regular region rainfall stations) (c) Removing the parts outside irregular and (d) Contour map filling region

The following principles are used for value assignment of equivalent areas (Han et al., 2006):

For closed regions, the closed contours are already lined up according to their minimal abscissa values and their inclusion relations to each other can be determined. This is achieved by analysis of two adjacent closed contours sequentially to decide if one is inside another. If so, when their elevation values are equal, this value is then assigned to the closed region. Otherwise, their average is assigned instead. If not, it indicates that the closed region is just formed by the closed contour without enclosing other closed contours, therefore the elevation of the closed contour will be assigned to the closed region.

As for a non-closed region, if it is composed of rectangular boundary with only closed region in it, its value is 0. If there is a nonzero elevation of its two surrounded contours, assign the value to the non-closed region. If there are two different elevations, then assign their average as the value of the non-closed area.


As a test example the rainfall data of Dalian City in a month is analyzed. After discrete data gridding and finishing contour tracking, a Google's Kml file is created using the resulting data and the contours are drawn on the map by calling the Kml file with Openlayers API interface. Figure 6 illustrates the drawing and filling process of Dalian City’s rainfall map. The effect of irregular region contours generation and filling based on the proposed approach is highly satisfying.


In this study, a new approach for contour drawing of irregular regions is proposed. This approach generates contours in the rectangular area which just contains the irregular regions of interest. Because it is based on established rectangular grid method, it needs less computation than the previously published studies. An improved ray method to determine the spatial relation between the equivalent points and irregular region is also presented which has efficient algorithm. Experimental results show that the proposed method could satisfactorily implement the contour search, drawing and filling of any irregular regions.


This study was financially supported by the key support discipline construction Project "Communication and Information Engineering" of China University of Geosciences (Beijing) and the Fundamental Research Funds for the Central Universities (Grant No. 2652013112).

1:  Chen, R., J. Zhou and L. Yu, 2007. Fast method to determine spatial relationship between point and polygon. J. Xi'an Jiaotong Univ., 41: 59-63.
Direct Link  |  

2:  Fu, X.D., J.S. Zhao and S.X. Ma, 2006. Isoline making algorithm based on D-TIN in the irregular regions. J. Kunming Univ. Sci. Technol., 31: 46-50.
Direct Link  |  

3:  Han, L., H. Shi and Q. Zhang, 2006. An algorithm for the region filling of isograms based on border point tracing. Comput. Eng. Sci., 28: 66-68.
Direct Link  |  

4:  Tang, L., A.J. Xu and L.M. Fang, 2008. The algorithm of creating contour lines based on DEM. Proceedings of the SPIE Geoinformatics 2008 and Joint Conference on GIS and Built Environment: Advanced Spatial Data Models and Analyses, Volume 7146, June 28-29, 2008, Guangzhou, China -.

5:  Wang, T., 2008. An algorithm for extracting contour lines based on interval tree from grid DEM. Geo-Spatial Inform. Sci., 11: 103-106.
Direct Link  |  

6:  Wang, T., 1991. Irregular area isoplenth drawing. J. Fuxin Mining Inst., 10: 71-76.

7:  Zheng, Y., C. Yao, C. Zhang and B. Liu, 2010. Fast area filling algorithm based on topology direction of contour line. Oil Geophys. Prospect., 45: 899-908.
Direct Link  |  

©  2021 Science Alert. All Rights Reserved