INTRODUCTION
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.
METHODS
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 gridbygrid 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 gridbygrid, 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
gridbygrid. 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 nonclosed 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, b_{1}~b_{12}: 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 nonclosed contours with the boundaries
(b_{1}~b_{12} 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 nonclosed contours with the boundaries: Once
point
has been removed, the contour line does not intersect with side b_{7},
b_{8}, which affects the contour filling. Therefore, it is necessary
to find the intersection of line
and side b_{7}, b_{8}. 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 nonclosed
region which is formed by nonclosed 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 nonclosed region, since the start and end points of nonclosed 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
nonclosed contour has two intersections with rectangular boundary, each intersection
point will participate in the construction of nonclosed area twice, while a
boundary vertex only participates in equivalent area construction once.

Fig. 4: 
Sorting boundary points of nonclosed contours counterclockwise 

Fig. 5: 
Numbering boundary points of nonclosed 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).
Nonclosed areas do not overlap, but may contain closed areas. Closed areas
may also contain closed areas. Therefore the nonclosed 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(ad): 
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 nonclosed 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 nonclosed region. If there
are two different elevations, then assign their average as the value of the
nonclosed area.
RESULTS
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.
CONCLUSION
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.
ACKNOWLEDGMENTS
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).