GENERALIZED EUCLIDIAN
LINEOFSIGHT
ON REGULAR HEXAGONAL GRIDS
by
Ron C. Boyd
B.S., Texas A&M University, 1997
A thesis submitted to the
University of Colorado Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Computer Science
2008
2008 by Ron C Boyd
All rights reserved
This thesis for the Master of Science
Degree by
Ron C. Boyd
Has been approved
By
Tom Altman
Date
Mike Radenkovic
Boyd, Ron C. (M.S., Computer Science)
Generalized Euclidean Lineofsight on Regular Hexagonal Grids
Thesis Directed by Professor Boris Stilman
ABSTRACT
While rasterization on hexagonal grids is a problem that has been solved by an
adaptation of Bresenham's line rasterization algorithm, lineofsight in hexagonal
space has received far less study. In a lineofsight situation, it is desirable that all
hexagons that intersect a line segment on a hexagonal grid be checked and in some
situations that the hexagons that intersect the segment be checked in order from the
beginning to the end of the segment. The algorithm presented solves the lineofsight
problem regardless of what indexing scheme is used on the hexagonal grid and is at
least one order of magnitude faster than previously published algorithms.
This abstract accurately represents the content of the candidates thesis. I recommend
its publication.
Signed
Boris Stilman
CONTENTS
Figures.....................................................................vii
Chapters
1. Introduction.............................................................1
2. Hexagonal Grids..........................................................2
2.1 Orientation and Indexing.............................................2
2.2 The Tex Grid.........................................................3
2.3 Stepping.............................................................5
2.4 Vancouver Distance...................................................6
2.4.1 Vancouver Rings.................................................6
2.4.2 Vancouver Squares...............................................6
2.5 Average Distance.....................................................7
2.5.1 Number of Cells in a Given Vancouver Ring.......................7
2.5.2 Total Hexagons Inside a Vancouver Ring..........................9
2.5.3 Average Distance to All Hexagons in a Vancouver Square..........9
2.6 Primary Axis.......................................................11
2.7 Parallel Lines......................................................12
2.7.1 Point Case.....................................................13
2.7.2 Flat Case......................................................14
2.8 Stepping Diameter...................................................15
3. Lineofsight Algorithm.................................................16
3.1 Determining Which Case..............................................16
3.2 Intersection Distances..............................................17
3.3 Point Case..........................................................19
3.3.1 Stepping Determination.........................................19
3.3.2 Special Cases..................................................21
3.3.2.1 Ray Crossing a Parallel Line in a Hexagonal Wall..............21
3.3.2.2 Relative Slope Equal to 0.....................................21
3.4 Flat Case...........................................................22
3.4.1 Stepping Determination.........................................23
4. Modifications of the Algorithm..........................................25
4.1 NonCentered Start and End Locations................................25
4.2 Extension to Threedimensions.......................................25
v
4.3 Mirroring Quadrant Four Onto Quadrant One.........................29
5. Empirical Results......................................................31
5.1 Performance.......................................................31
5.2 Steps per Calculation.............................................32
5.2.1 Theoretical.................................................33
5.2.1.1 Point Case..................................................33
5.2.1.2 Flat Case...................................................34
5.2.2 Experimental................................................34
5.3 Comparison........................................................35
6. Application............................................................36
7. Conclusion.............................................................38
Appendix
A. Hexagonal LineofSight Algorithm......................................39
Bibliography
BIBLIOGRAPHY..............................................................55
vi
LIST OF FIGURES
Figures
Figure 2.1 Hexagonal direction and vertex indexing............................2
Figure 2.2 Hexagonal distances................................................3
Figure 2.3 Hexagonal axes.....................................................4
Figure 2.4 Tex grid axes......................................................4
Figure 2.5 Hexagonal neighbor indexing........................................5
Figure 2.6 Vancouver rings....................................................7
Figure 2.7 Equilateral triangle from hexagonal neighbors.......................8
Figure 2.8 Extension of equilateral legs......................................8
Figure 2.9 Comer of Vancouver square with odd iindex on center...............10
Figure 2.10 Comer of Vancouver square with even iindex on center.............10
Figure 2.11 Primary axis directions...........................................12
Figure 2.12 Point case parallel lines.........................................12
Figure 2.13 Flat case parallel lines..........................................13
Figure 2.14 Point case example................................................14
Figure 2.15 Flat case example.................................................14
Figure 3.1 LOS case indices...................................................16
Figure 3.2 Point case primary distances.......................................18
Figure 3.3 Flat case primary distances........................................19
Figure 3.4 Point case stepping values.........................................20
Figure 3.5 r crossing a parallel line and a hexagonal boundary simultaneously.21
Figure 3.6 Relative slope equal to zero.......................................22
Figure 3.7 Flat case stepping example.........................................23
Figure 3.8 Flat case stepping values..........................................24
Figure 0.1 150 degree multiplier determination................................26
Figure 0.2 30 degree multiplier determination.................................28
Figure 0.3 60 degree multiplier determination.................................29
Figure 0.4 120 degree multiplier determination................................29
Figure 5.1 Average hexagons per second by length..............................32
Figure 6.1 LG Architecture....................................................37
vii
1. Introduction
Hexagonal grids have very appealing qualities in several areas of computation. In
path finding, the fact that all neighbors of a hexagon share an edge with the hexagon
and form a closed neighborhood around the hexagon [2] is very appealing for
reducing graphical artifacts associated with squarebased grids [7] and for increasing
the fidelity of paths produced without increasing the branching factor of search
algorithms[l]. In computer graphics, Rogers [9] showed that hexagonal grids offer
better solutions than squarebased grids for sampling a plane with a gridbased
structure.
As increasingly more software begins to work with underlying structures based on
hexagonal grids, efficient algorithms for rasterization of lines and curves have been
developed and studied [5]. The line rasterization techniques have mostly been based
on Bresenhams rasterization algorithm and are very efficient for providing a path
between cells of a hexagonal grid; however, the algorithm does not return every
hexagon that the line being rasterized intersects.
For lineofsight operations, every hexagon that intersects a lineofsight ray needs to
be checked. Very little research regarding hexagonal lineofsight has been
published. In fact, all hexagonal lineofsight references that were found referred to
work done by Verbrugge [4] and a suggested algorithm for finding all hexes in a grid
that intersect a given line. The basic operation of Verbrugges algorithm is to test the
comers of the starting hexagon to see which points lay on which side of the lineof
sight ray and use that information to step to the next hexagon along the path, where
the process is repeated. While the algorithm can be made to find all hexagons that
intersect a line of sight ray, its inefficiency lies in the fact that it fails to use the
symmetry of the hexagonal grid to its advantage. The symmetry and regularity of a
hexagonal grid is what allows the algorithm presented in this paper to reduce all cases
to two cases and supply exact intersection distances with minimal calculation.
1
2. Hexagonal Grids
2.1 Orientation and Indexing
For the purposes of this paper, hexagons will always be oriented with a flat side up.
Sides will be indexed starting with zero pointing up and increment in a clockwise
direction around the hexagon. Vertex indexing will be determined such that the
vertex at the clockwise end of an edge shares the index of that edge.
o
Figure 2.1 Hexagonal direction and vertex indexing
To more easily describe distances that are relevant to hexagons, let us define two
variables that will correspond to the inside and outside radii of the hexagon. Let hexa
be the outside radius of the hexagons, and hex}, be the inside radius.
2
Figure 2.2 Hexagonal distances
hexh = hexa
2
2.2 The Tex Grid
Standard axes on hexagonal grids are usually separated by either sixty or one hundred
and twenty degrees. These types of axes make many things about hexagonal grids
quite convenient. Steps along a given axis are always in the same hexagonal
direction [3, 6],
3
Figure 2.3 Hexagonal axes
The problem with those types of axes is that representing a square area becomes a
difficulty. Solving the square area representation was the initial purpose of the tex
grid [1]. The axes of the tex grid are orthogonal. The iaxis of the tex grid alternates
up and down as it goes to the side.
4
While this allows the tex grid to represent square areas with the most fidelity that is
possible with a hexagonal grid, stepping along the iaxis is no longer as simple as
incrementing or decrementing the i value in an (i, j) pair. Whether or not the current i
value is even or odd determines how stepping should occur.
2.3 Stepping
When stepping is discussed in this paper, it denotes the process of incrementing the (i,
j) indices of the current cell. The algorithm presented was developed to work on a tex
grid and the GetDirection() function included is specific to that type of grid. For
other indexing schemes, the only significant changes to any of the code provided
should occur in the GetDirection() function. In Fig. 2.5 the indexing changes
resulting from stepping in any of the six primary directions is shown.
Figure 2.5 Hexagonal neighbor indexing
5
2.4 Vancouver Distance
The Vancouver Distance [1] on a tex grid is defined as follows:
Given two nodes (xl,yl) and (x2,y2),
jc =*1*2 
y=\\yly2\\
f xl(mod 2), ify\< y2 and x(mod 2) = 1
correction <
x2(mod 2), if yl > y2 and x(mod 2) = 1
VancouverDistance = max
0,T
x
2
+ x correction
A simple visual interpretation of the Vancouver Distance is to consider concentric
rings of hexagons around an origin.
2.4.1 Vancouver Rings
Let a Vancouver ring of size n be the set of all hexagons in a hexagonal grid that are
Vancouver distance n from a given origin. Vancouver rings are shown in Fig. 2.6.
2.4.2 Vancouver Squares
Let a Vancouver square of size n be a square region of hexagonal cells that just
contains a Vancouver ring of size n.
6
Figure 2.6 Vancouver rings
2.5 Average Distance
Finding the average distance from the origin to all of the cells in a Vancouver square
of size n is not trivial. Finding the average distance from the origin to every cell
within a Vancouver ring of radius n is trivial. The only thing that need be known is
the number of cells in each Vancouver ring.
2.5.1 Number of Cells in a Given Vancouver Ring
It is easy to prove that there are 6n cells in a Vancouver ring of radius n. First note
that connecting the centers of three adjacent hexagons in a hexagonal grid yields an
equilateral triangle.
7
Figure 2.7 Equilateral triangle from hexagonal neighbors
If two of the legs of the equilateral triangle are then extended by the diameter of a
hexagon, 2 hex},, each of those legs is extended along one of the primary directions
for a hexagon and will therefore end up at the center of another hexagon. Since two
legs of the triangle have been extended by one hexagonal diameter without changing
angles, the third leg must also be extended by one hexagonal diameter to maintain the
triangle. Since the third leg is also in a primary hexagonal direction, adding 2 hexh
along that leg will add one more cell to that leg.
8
Since the triangle given represents one sixth of the Vancouver ring, each Vancouver
ring has 6n hexagons.
2.5.2 Total Hexagons Inside a Vancouver Ring
The total number of hexagons that are within a Vancouver ring of size n is:
n
1 + ^ 6/ = 1 + 3 (n)(n +1)
2.5.3 Average Distance to All Hexagons in a Vancouver Square
The average distance from the origin of all cells out to and including the nth
Vancouver ring is:
This does not include the hexagons in the comers of the Vancouver square that are
above and below the sides of the nth Vancouver ring. In order for those hexagons to
be accounted for, we must first understand those areas. If the center of the Vancouver
rings has an odd iindex value, the lower left comer area will begin with n hexagons
in a row and decrement by two through successive rows until the comer is reached.
n
2()(rt + l)(2 + l)
l + 30)( + l)
i=i
n
9
Figure 2.9 Corner of Vancouver square with odd iindex on center
If the center has an even iindex value, the rows will start with n1 hexagons in a row
and decrement by two until the comer is reached.
Figure 2.10 Corner of Vancouver square with even iindex on center
The gaps at the top of the square area have their behavior reversed from those at the
bottom based on whether or not the center has an even or odd iindex value.
Therefore, the overall number of extra hexagons in the square area just containing the
nth Vancouver ring is:
10
I fl+1 I
L 2 J L2J
2* Â£ (77 2/ + 2) + 2 ^ (77 2/ +1)
i=l 7=1
To find the average distance from the center of a square area containing the nth
Vancouver ring to every hexagon in that square area is then:
n+1
L 2 J L2J
77(77 +1)(2/7 +1)+2 ^ [(77 2i + 2)(n + /)] + 2 ^ [(77 2/ + 1)(t7 + /)]
_________________________1=]_____________________________M________________________
1 + 3(77)(77 + 1)+2 Â£ (772/ + 2) + 2*Â£(772/ + l)
/=1 /=1
2.6 Primary Axis
For any given lineofsight ray r, there is a reference ray ro that is used to determine
which case will be used in the lineofsight algorithm. Reference rays, hereby known
as primary axes, occur every thirty degrees in the quadratic plane with zero degrees
having the standard mathematical definition of being down the xaxis in the positive
direction. Reference rays originate from the center of the hexagonal cell that contains
the beginning of r.
11
Figure 2.11 Primary axis directions
2.7 Parallel Lines
The base concept of this algorithm is that there are two sets of parallel lines that
dictate which stepping behavior should be used when moving through the hexagonal
grid. Given a ray from the center of a hexagonal cell, the ray will leave the cell
within fifteen degrees of either a comer of the hexagonal cell, or within fifteen
degrees of the middle of a side of the hexagonal cell. If the ray is within fifteen
degrees of a comer, it will be referred to as a point case and otherwise will be
referred to as a flat case. For each case, there is an associated set of parallel lines
that constitute behavioral boundaries for the algorithm.
12
Figure 2.13 Flat case parallel lines
2.7.1 Point Case
For rays that leave the starting cell within fifteen degrees of comer c, the set of
parallel lines is determined by edges (c + 5) mod 6 and (c + 2) mod 6 of every
hexagon in the grid. In this case, all of the parallel lines are spaced evenly and are
one half of the inside diameter of the hexagon apart.
13
Figure 2.14 Point case example
2.7.2 Flat Case
For rays that leave within fifteen degrees of the middle of the edge of a hexagon, the
parallel lines are the set of lines that pass through a comer of a hexagon in the
hexagonal grid, but are orthogonal to the edge that the ray exited through. Unlike the
point case, the parallel lines of the flat case are not evenly spaced. The spacing
between the lines alternates between one half of the outside diameter of the hexagons
and one quarter of the outside diameter of the hexagons.
Figure 2.15 Flat case example
14
2.8 Stepping Diameter
In order to determine how many steps should be taken between parallel line crossings,
the distance along r must be normalized by some value related to the effective
diameter of the hexagons. Let that diameter be known as the stepping diameter. For
the point case, the stepping diameter is 1.5 hexa, and for the flat case it is 2 hext,.
15
3. Lineofsight Algorithm
3.1 Determining Which Case
While a lineofsight ray r is associated with either a flat or point case, each case must
be further broken down to being either clockwise or counterclockwise from the
primary axis ro of that specific case as the rotational orientation will affect the
stepping of the algorithm. For the purposes of this paper, the individual cases will
have the indices indicated in Fig. 3.1. Cases in quadrants two and three of the
quadratic plane are mirror images of the cases in quadrants one and four and can be
computed by switching the iindex incrementing value from 1 to 1. It is also
possible to collapse quadrant four into quadrant one through the symmetry exhibited
across the xaxis, but for clarity of the algorithm presented, that step is not shown.
For the purposes of this paper, let it be assumed that all r have Ay >= 0.
Figure 3.1 LOS case indices
16
Point Cases Flat Cases
0 1
3 2
4 5
6 7
9 8
10 11
Table 3.1 Point/Flat case breakdown by index
To determine which case r belongs to, a comparison of the slope of r to known slopes
of primary axes is used. This step can be done very efficiently since the slopes of the
primary axes are constant.
3.2 Intersection Distances
Since the algorithm steps through the cells in a discrete manner, it needs to know how
many steps to make before changing its behavior. Determination of the distance
along the primary axis is straightforward. It is necessary to first determine the slope
of the lineofsight ray relative to the primary axis. Let us call this relative slope mrei.
Using the trigonometric identity of Tan(A B), we find:
Given:
mo = slope of r
mj = slope of ro
If the lineofsight ray is counterclockwise from the primary axis:
Else:
(mo m\)
mrei =
(1 + mom\)
17
(mi mo)
nirei =
(1 + momi)
Since spacing between parallel lines is known, the distance from ro to the first parallel
line is a constant for each case. Let us call that height heightprimary Finding the
distance to the first parallel line projected to ro is achieved by dividing the
heightprimary by nirei
disti
heightpn/
primary
primary
nirei
Figure 3.2 Point case primary distances
For the point case, heightprimary is one half of the inside diameter of a hexagonal cell.
For the flat case, heightprimary is one quarter of the outside diameter of a hexagonal
cell.
18
Figure 3.3 Flat case primary distances
distprimary is divided by the stepping diameter to find the number of discrete hexagonal
steps along the primary axis between parallel line intersections, which shall be known
as distx.
distx =
disti
primary
stepping diameter
3.3 Point Case
Given a lineofsight ray r that passes out of the originating hexagon near comer c,
the algorithm will need to step through the hexagonal array by alternating moving in
directions c and n,n = (c + 1) mod 6. If the r exits clockwise from vertex c, then the
algorithm will start with direction n, otherwise it will start with direction c. When r
crosses one of the parallel lines, the alternation between c and n reverses.
3.3.1 Stepping Determination
For the purposes of stepping, the ro has value zero at the comer c that the lineofsight
ray is exiting the originating hexagon near, not at the origin of r. Integer increments
along ro are spaced at the distance of the outside diameter of the hexagons. This
leaves the center of the originating cell at negative twothirds on ro. The origin of r is
not used as the origin of ro to solve the problem of cell boundaries along parallel lines
19
alternating between being onehalf of the outside diameter of a hexagon apart and a
full outside diameter of a hexagon apart.
Let distancehex be the value that will be used to determine how many steps should be
taken until one of the parallel lines is crossed. Initially, distancehex will be negative
twothirds plus distx. Each time a parallel line is crossed, distancehex is incremented
by distx.
If the lineofsight ray is crossing the first parallel line in Fig 3.1 at any point in the
range:
the algorithm continues in a normal manner. If the crossing is in the range:
AND the algorithm's next step is going to be in the n direction, then a special case has
arisen and the algorithm must make a special step.
Figure 3.4 Point case stepping values
20
3.3.2 Special Cases
3.3.2.1 Ray Crossing a Parallel Line in a Hexagonal Wall
A special case arises when r crosses one of the parallel lines while simultaneously
passing through the side of a cell. The point circled in Fig 3.5 shows an occurrence
of this special case. The algorithm will need to move perpendicular to the parallel
lines, and not do the normal n, c swap for the direction of the next step that would
occur when a parallel line is crossed.
Figure 3.5 r crossing a parallel line and a hexagonal boundary simultaneously
3.3.2.2 Relative Slope Equal to 0
In cases where mrei is zero, every other crossing of a hexagonal boundary is going to
have a pair of hexagons to check. The algorithm will simply alternate between c and
n, using direction c first, and do a separate check in the direction perpendicular to and
across ro after each step in direction c.
21
Figure 3.6 Relative slope equal to zero
3.4 Flat Case
In contrast to the behaviors of the point case, the flat case has two different behaviors
For a lineofsight ray r exiting the originating cell in direction c, while r stays
between a larger spaced set of parallel lines, the algorithm simply steps through the
hexagons in the direction c. Once the ray passes into the region between two closely
spaced parallel lines, it alternates between directions q, q = (c + 5) mod 6 and n, n =
(c + 1) mod 6.
22
Figure 3.7 Flat case stepping example
3.4.1 Stepping Determination
For the purposes of stepping, r0 has value zero at the middle of the edge c that r
passes through. Integer increments along the primary axis are spaced by the inside
diameter of the hexagons.
23
Figure 3.8 Flat case stepping values
Let distancehex be the distance along the primary axis that will be used to determine
how many cellular steps to make, distancehex will initially be negative distx. The
negative offset allows the algorithm to treat the initial case as a regular case where 2 *
distx is added to distancehex for stepping between the more widely spaced parallel
lines, distancehex is incremented by only distx between the more closely spaced
parallel lines. Additionally, when the lineofsight ray passes from a narrow parallel
band to a wide parallel band distancehex is decremented by one half because of the
relative offset in the direction of the primary axis of the perpendicular walls of the
next set of hexagonal cells.
24
4. Modifications of the Algorithm
4.1 NonCentered Start and End Locations
Extension of this algorithm to handle start and end locations that are not cellcentered
is relatively easy. Case determination, stepping, and distance between parallel line
crossings is handled in the same manner as with cellcentered endpoints. Only the
initial value of distancehex needs to be modified for the algorithm to run normally.
distance^ is modified by the difference between the first parallel line intersection
points of r and a line with a slope equal to that of r that is cellcentered.
For the flat cases, it is also possible for the origin of a lineofsight ray to fall within
one of the more closely spaced sets of parallel lines. In this situation it is necessary to
swap the type of checking that the algorithm begins with. Also, distancehex would
begin with an initial value of zero instead of negative distx. Finally, the origin of the
cellcentered reference line used to find the intersection point would need to be
translated in a direction perpendicular to the primary axis to the parallel line that is
closer to the cell center than the origin of r. Once the initial conditions have been set
for distancehex, the algorithm can continue with no modification.
4.2 Extension to Threedimensions
Extending the algorithm to include layers of hexagonal prisms is not hard if you rely
on the symmetry of the hexagonal grid. For any given lineofsight ray r, there are a
limited number of different angle lines r will cross.
For a flat case with cellcentered begin and end points, r starts in a more widely
spaced band of parallel lines and crosses only cell boundaries that are orthogonal to
the primary axis ro. The distance along r between each of those crossings is a
constant value. Once r crosses into one of the more closely spaced parallel line sets is
where symmetry shines through to save time on calculation. At first glance it appears
that there are many line crossings that would need to be checked to find the distance
25
to intersections along r. Since those lines all have constant slopes of either thirty or
one hundred and fifty degrees relative to ro, we need only calculate a single value per
slope that can be multiplied by the distance from where r crosses into the band and
the next stepping diameter increment.
In order to find the intersection points along r, we need to calculate the length of the
hypotenuse of the triangle formed by distanceprimary and heightprimary This calculation
need only be done once per lineofsight check. Let us call this hypotenuse value
hypotenuseprimaiy. As the algorithm is executed on r, the distance along r at each
intersection point with a parallel line is easily found by adding hypotenusePnmary onto
a tracking variable each time r passes out of a narrower parallel band and adding
double hypotenuseprimary each time r passes out of a wider parallel band.
Figure 0.1 150 degree multiplier determination
In Fig. 4.1, r is passing into and through a narrower parallel line set. r enters at A and
intersects both of the hexagonal boundaries originating at B. No matter how near or
far A is from B, the angles will never change for a given r. Using the law of sines:
K Xl50
sin(30$) sin(150)
Solving for x and simplifying constants yields a value that can be multiplied by the
appropriate K value for any a ABC to find x since all a ABC will be similar.
26
*150 = *K
sin(30$)
Since we have only the slope of r and do not actually know the angle 6, we can
substitute
Xl50
2
tan(30#)
^1 + tan2 (300)
K
For simplicity of reading, let:
ox tan(30)tan(6?)
mx so = tan(30 G) =
l + tan(30) tan(#)
We now have x expressed solely in terms that we have already calculated for the line
ofsight algorithm. Mi so is composed of constants and the value tan(d) which is mrej.
Xl50 =K
mi so
yj\ + mi50 777150
The other case from the narrow band is identical except for the constants.
Figure 0.2 30 degree multiplier determination
X30 =
1
2
sin(15O0)
K
Which leads us to:
And finally:
mao = tan(150#) =
tan(150)tan(#)
l + tan(150) tan(#)
*30 =2K
mm
Vl+ 777 30 77730
The same process is used to find the multiplicative values relevant to the point cases.
28
120H
Figure 0.3 60 degree multiplier determination
Figure 0.4 120 degree multiplier determination
In order for these distances to be useful, we need to know the amount of total
elevation change between the start and end of r and normalize that difference over the
length of r so that elevations can be determined by simply multiplying the current
distance along r by that normalized value.
For the point case special case of r crossing a parallel line at the same time as it
crosses a hexagonal boundary, the length along the ray is trivial since it is at a parallel
border and is therefore an integer multiple of hypotenusePnmary
4.3 Mirroring Quadrant Four Onto Quadrant One
29
Collapsing cases six through eleven into cases zero through five would not be a
challenging task. In the algorithm presented, quadrants two and three have already
been collapsed into quadrants one and four through the use of the variable m_istep,
which is one for rays heading to the right of the yaxis and negative one for rays
heading to the left of the yaxis. If a similar variable m jstep were used, only
modifications to the function GoDirection() would need be made.
Through the process of reduction of cases to result in the two cases that are presented
here, it was found that the reduction in complexity results in a decrease in efficiency.
Depending upon the length of the lineofsight check to be made, versions of this
algorithm that have been fully expanded run anywhere from thirtysix to sixty percent
faster. As such, further reductions were deemed unnecessary. Table 5.1 shows
relative performance of several versions of the algorithm.
30
5. Empirical Results
5.1 Performance
For Vancouver squares of increasing size, the algorithm given approaches a
maximum number of hexagons processed per second. This implies that the runtime
of the algorithm is O(n) with respect to the length of the lineofsight check.
Studying the content of the algorithm confirms that position. There are a few
calculations that need be made at the beginning of a lineofsight check no matter the
length and then the algorithm steps through the cells with a minimum of calculation.
As n increases, the cost of those few starting calculations and the price of calling the
function itself become overshadowed by the speed of the stepping. Table 5.1 shows
the average hexagons per second of three versions of the algorithm. Algorithm
refers to the algorithm as it is presented here, with function calls in the body. No
Funcs refers to the algorithm presented here with all of the work being done by the
function calls in Algorithm inplace. Expanded refers to the algorithm in its
largest form with each of the directional cases having its own section in code with no
function calls. The relative performance of Verbrugges lineofsight algorithm is
presented as Verbrugge to show the relative performance of the two algorithms.
The higher startup costs of the algorithm presented here versus Verbrugges
algorithm are shown by the curved appearance of the datasets associated with
Algorithm, No Funcs, and Expanded.
31
Figure 5.1 Average hexagons per second by length
5.2 Steps per Calculation
Since each case deals with relative slopes between zero and fifteen degrees, let us
assume that the average relative slope dealt with by the algorithm is seven and one
half degrees. We should expect to see an average number of steps between parallel
lines of distx.
, distprimary
distx =
stepping diameter
_ heightp, imury
CtlStprimary 
ITlrel
32
5.2.1 Theoretical
5.2.1.1 Point Case
For a point case, heightprimary is hexh
_ hexh
CtlSlprimary 
Wrel
>/3
The identity hexh = ^hexa then gives us:
hexa
dish
primary '
YYlrel
hexa
distx 
ttlrel
stepping diameter
Point cases have a stepping diameter of
hexa
2
disL = mn'
hexa
2
distx =
3 mrel
33
5.2.1.2 Flat Case
For a flat case, heightprimary is
hexa
hexa
dish
primary
2
m,el
Flat cases have a stepping diameter of 2 hex/, which leads us to:
hexa
2
distx =
Mrel
2 hexh
Using the hexf, identity mentioned above:
hexa
~Y~
distx = mrxl = i
hexayJ3 2J 3 mrei
The flat case is always either in a narrow band, where it makes two steps per stepping
diameter or its in a wide band where distx is doubled, so the number of steps is
always double the value calculated.
distx 
s
mrel
5.2.2 Experimental
By modifying the algorithm to count the number of times that distancehex is
incremented and the number of steps that are taken, an average number of steps per
34
parallel line crossing can be found. For a Vancouver square of size one thousand, the
algorithm averaged 4.278939 steps per parallel line crossing. Based on the results of
sections 5.2.1 and 5.2.2, we expect to see:
= 4.38541
V3 tan(7.5)
There is an error of approximately 2.43 percent between the expected number of steps
per parallel line crossing and what experiments showed to be the actual number of
steps per parallel line crossing. The difference corresponds to using an average angle
of 7.68 degrees instead of the assumed 7.5 degree average. Due to discretization, an
even distribution of slopes is only possible at infinite range, so an average angle that
is slightly higher than expected is not alarming.
5.3 Comparison
When comparing the runtimes of the algorithm given here to Verbrugges hexagonal
lineofsight algorithm, this algorithm is at least one order of magnitude faster.
Additionally, extension of this algorithm to three dimensions is easier due to
controlling the angles of intersection and knowing exactly what type of intersection is
going to occur next.
35
6. Application
The lineofsight algorithm presented here was developed for use in the Stilman
Advanced Strategies LG Package software. When Stilman won DARPA program
contract NBCHC040152, this algorithm was essential for being able to compute a
great deal of lineofsight information in a reasonable amount of time. Areas up to
six hundred and twentyfive square kilometers were rasterized onto a hexagonal grid
with the hexagons having an inside diameter of twentyfive meters. The complete
size of the resulting array of hexagonal prisms is 32,925,000 cells. Lineofsight
information was needed for the entire area.
In Fig 6.1, the LG LOS component, which is based on the algorithm presented here,
is shown to be a very lowlevel component of the overall LG system. Extensive
testing and verification of the lineofsight algorithm has been done over the course of
several years and fault has yet to be found with the correctness of the algorithm.
36
Figure 6.1 LG Architecture
37
7. Conclusion
The lineofsight algorithm introduced here is an O(n) solution to checking complete
lineofsight. Reliance of this algorithm on the symmetry and regularity of a
hexagonal grid via the dual sets of parallel lines give it an edge over other lineof
sight algorithms. Efficient operation of this algorithm allowed Stilman to not only
fulfill the RAID program requirements, but to exceed them and best human
commanders in battlefield simulations. It is even feasible that this algorithm could be
used in the world of image processing, where the benefits of hexagonal grids are well
known.
38
APPENDIX
A. Hexagonal LineofSight Algorithm
void CHexRasterizer::GoDirection(int dir, int &curi, int &curj)
{
GetDirection(dir, curi, curj, curi, curj);
}
void CHexRasterizer::GetDirection(int dir, int &curi, int &curj, int &tempi, int &tempj)
{
if(m_istep < 0)
{
if(dir = 1)
{
dir = 5;
}
else if(dir = 2)
{
dir = 4;
}
else if(dir = 5)
{
dir = 1;
}
else if(dir = 4)
{
dir = 2;
}
}
switch(dir)
{
case(O):
{
tempi = curi;
tempj = curj + 1;
break;
}
case(l):
{
tempi = curi + 1;
tempj = curj + !(tempi % 2);
break;
}
case(2):
39
{
tempi = curi + 1;
tempj = curj (tempi % 2);
break;
}
case(3):
{
tempi = curi;
tempj = curj 1;
break;
}
case(4):
{
tempi = curi 1;
tempj = curj (tempi % 2);
break;
}
case(5):
{
tempi = curi 1;
tempj = curj + !(tempi % 2);
break;
}
}
}
void CHexRasterizer::LOSCheckAlgorithm(int bi, int bj, int ei, int ej, CString filename)
{
m_curi = bi;
mcurj = bj;
m_ei = ei;
m_ej = ej;
//get rectangular coordinates
GetRectangularCoordinates(m_curi, m_curj, m_xl, m_yl);
GetRectangularCoordinates(m_ei, m_ej, m_x2, m_y2);
//vector of the line segment
m ray.x = m_x2 m_xl;
m ray.y = m_y2 m y 1;
m_ray.z = 0;
m_rayLength = m_ray.length();
#ifdef EXPORTMATLABVISIBILITYCHECK
CMatlabOutput mout;
CString tempf;
tempf.Format("%s\\out%d_%d_to_%d_%d.m", filename, bi, bj, ei, ej);
mout.SetFile(tempf);
mout.OpenFile();
40
mout.RenewFile(FALSE);
#endif
//find slope of ray
if((m_ray.x > SLOPETOL)  (m_ray.x < SLOPETOL))
{
m_slope = m_ray.y / m_ray.x;
}
else if(m_ray.y > 0)
{
m_slope = 2000;
}
else
{
m_slope = 2000;
}
//now find direction
m_direction = 1;
if(m_ray.x >= 0) // quads 1 and 4
{
m_istep = 1;
}
else
{
//force the slope to first and fourth quads for direction finding
mslope *= 1;
m_istep = 1;
}
if(m_slope >= 0) // quad 1
{
myinc = 1;
if(m_slope >= SQ3INV)
{
if(m_slope >= SQ3)
{
if(m_slope >= TAN75)
{
m_relativeSlope = 1.0 / m_slope;
indirection = 5; // [75, 90]
{
m_subtype = FALSE;
m_relativeSlope = (m_slope SQ3) / (1.0 + m slope SQ3)
indirection = 4; // [60, 75)
}
}
else if(m_slope >= 1.0)
41
{
msubtype = TRUE;
m_relativeSlope = (SQ3 m_slope) / (1.0 + SQ3 m_slope);
in direction = 3; // [45, 60)
m_relativeSlope = (m_slope SQ31NV) / (1.0 + m_slope SQ3INV);
indirection = 2; // [30,45)
}
}
else if(m_slope >= TAN15)
{
mjelativeSlope = (SQ3INV m slope) / (1 + SQ3INV m_slope);
m_direction = 1; // [ 15, 30)
else
msubtype = FALSE;
mrelativeSlope = m_slope;
indirection = 0; // [0, 15)
}
}
else // quad 4
myinc = 1;
if(m_slope <= SQ3INV)
{
if(m_slope <= SQ3)
{
if(m_slope <= TAN75)
{
m_relativeSlope = 1.0/ m_slope;
indirection =11;// [75, 90]
}
else
m_subtype = FALSE;
m_relativeSlope = (SQ3 + m slope) / (SQ3 m slope 1.0);
m_direction = 10; // [60, 75)
}
}
else if(m_slope <= 1.0)
{
m_subtype = TRUE;
m relativeSlope = (SQ3 + m_slope) / (1.0 SQ3 m slope);
m direction = 9; // [45, 60)
else
{
42
mrelativeSlope = (SQ3INV m slope) / (1.0 SQ3INV m_slope);
in direction = 8; // [30, 45)
}
}
else if(m_slope <= TAN 15)
{
m relativeSlope = (m_slope + SQ31NV) / (1.0 SQ3INV m slope);
m_direction = 7; // [15, 30)
}
else
msubtype = TRUE;
m_relativeSlope = mslope;
myinc = 1;
m_direction = 6; // (0, 15)
}
}
#ifdef EXPORT_MATLAB_VISIBILITY_CHECK
mout.BeginArray("x");
mout. AddArrayPoint(m_x 1);
mout.AddArrayPoint(m_x2);
mout.EndArray();
mout.BeginArray("y");
mout. Add ArrayPoint(m_y 1);
mout.AddArrayPoint(m_y2);
mout.EndArray();
mout. Plot Arrays(2);
mout.PrintHexagonFromCoordinates(m_curi, m curj, 0);
#endif
CHECKCELL(m_curi, m curj);
switch(mdirection)
{
//FLAT CASE
case(l): //15..30 degrees
case(2): // 30..45 degrees
case(5): // 75..90 degrees
case(7): // 30..15 degrees
case(8): // 45..30 degrees
case(l 1): // 90..75 degrees
{
if(m_relativeSlope < SLOPETOL)
{
//straight out the row of cells
#ifdef FAILS AFEMODE
m_farout = m_ray Length + m_hex_h2;
m_intpoint = m_hex_h;
while(m_intpoint < m_farout)
#else
43
#endif
while(l)
switch(m_direction)
{
case(l):
case(2):
{
GoDirection(l, m_curi, mcurj);
break;
}
case(5):
{
GoDirection(0, m_curi, m curj);
break;
}
case(7):
case(8):
{
GoDirection(2, m_curi, m_curj);
break;
}
case(l 1):
{
GoDirection(3, m_curi, m_curj);
break;
}
CHECKCELL(m_curi, m_curj);
#ifdef FAILSAFE_MODE
mintpoint += m_hex_h2;
#endif
}
#ifdef FAILSAFE_MODE
//if it reached this, there's trouble
CString temp;
temp.Format("Failed to reach target for (%d, %d) to (%d, %d)", bi, bj, ei, ej);
AfxMessageBox(temp);
#endif
}
else
// distance along primary axis to hitting the next parallel line
musedist = m_hex_a02 / mrelativeSlope;
//number of hexagons along primary axis between line crossings
m_distx = m_usedist / m_hex_h2;
//start at m_distx since it starts at a hex center
m_distance = m_distx;
m_type = TRUE;
44
m_stePPer = 0;
#ifdef FAILSAFEMODE
m_usedlen = sqrt(m_usedist m usedist + m_hex_a02 m_hex_a02)
m_totallen = musedlen;
mfarout = mray Length + m_hex_h2;
while(m_totallen < m_farout)
#else
#endif
while(l)
if(m_type)
{
m_distance += m_distx + mdistx;
//half is subtracted because the LOS
//ray starts halfway through a wide band
m dval = m distance 0.5;
for(;m_stepper
{
switch(mdirection)
{
case(l):
case(2):
{
GoDirection(l, mcuri, m_curj);
break;
}
case(5):
{
GoDirection(0, m_curi, mcurj);
break;
}
case(7):
case(8):
{
GoDirection(2, m_curi, m curj);
break;
}
case(l 1):
{
GoDirection(3, m curi, m_curj);
break;
}
}
CHECKCELL(m_curi, m curj);
}
#ifdef F AILS AFEMODE
m totallen += m_usedlen + m_usedlen;
#endif
45
}
else
//in a halfheight space between parallel lines
//m distx only added once
m_distance += m_distx;
//make one less step than called for so the last
//step can be in only one direction
m dval = m distance 1;
for(;m_stepper
{
switch( m_direction)
{
case(l):
{
GoDirection(2, m_curi, m_curj);
break;
}
case(2):
{
GoDirection(0, m_curi, m_curj);
break;
}
case(5):
{
GoDirection( 1, mcuri, m_curj);
break;
}
case(7):
{
GoDirection(l, m curi, mcurj);
break;
}
case(8):
{
GoDirection(3, m_curi, mcurj);
break;
}
case(l 1):
GoDirection(2, m curi, m_curj);
break;
}
}
CHECKCELL(m_curi, m curj);
switch( m_direction)
{
case(l):
46
{
GoDirection(0, mcuri, m_curj);
break;
}
case(2):
{
GoDirection(2, m_curi, mcurj);
break;
}
case(5):
{
GoDirection(5, m_curi, m_curj);
break;
}
case(7):
{
GoDirection(3, m_curi, m curj);
break;
}
case(8):
{
GoDirection(l, m curi, m_curj);
break;
}
case(l 1):
{
GoDirection(4, m curi, m_curj);
break;
}
CHECKCELL(m_curi, m curj);
}
//Make the last half of a step
switch(mdirection)
{
case(l):
{
GoDirection(2, m curi, m curj);
break;
}
case(2):
{
GoDirection(0, m_curi, m curj);
break;
}
case(5):
GoDirection(l, m_curi, m curj);
47
break;
}
case(7):
{
GoDirection( 1, m_curi, mcurj);
break;
}
case(8):
GoDirection(3, m_curi, m_curj);
break;
}
case(l 1):
{
GoDirection(2, m_curi, m_curj);
break;
}
CHECKCELL(m_curi, m curj);
//this is done because you've stepped into a new row
//and there is now a 1/2 cell offset
m_distance = 0.5;
#ifdef FAILSAFE MODE
#endif
m_totallen += musedlen;
}
m_type = !m_type;
}
#ifdef FAILS AFE_MODE
//if it reached this, there's trouble
CString temp;
temp.Format("Failed to reach target for (%d, %d) to (%d, %d)", bi, bj, ei, ej);
AfxMessageBox(temp);
#endif
}
break;
//POINT CASE
case(0): // 0..15 degrees
case(3): // 45..60 degrees
case(4): // 60..75 degrees
case(6): // 15..0 degrees
case(9): // 45..60 degrees
case(10): // 60..75 degrees
{
if(m_relativeSlope < SLOPETOL)
{
48
//have to alternate checking one then two cells
#ifdef FAILS AFE_MODE
m_farout = m_ray Length + m_hex_a2;
//the location of the cell boundary intersection
m_intpoint = mhexa;
while(m_intpoint < m_farout)
#else
#endif
while(l)
switch(m_direction)
{
case(O):
case(6):
{
GoDirection(l, m_curi, m_curj);
CHECKCELL(m_curi, mcurj);
GetDirection(3, m_curi, m_curj, m_ti, m_tj);
CHECKCELL(m_ti, m_tj);
GoDirection(2, m_curi, m_curj);
CHECKCELL(m_curi, m_curj);
break;
}
case(3):
case(4):
{
GoDirection(0, mcuri, m curj);
CHECKCELL(m_curi, m_curj);
GetDirection(2, m curi, m_curj, m_ti, m_tj);
CHECKCELL(m_ti, m_tj);
GoDirection(l, m_curi, m_curj);
CHECKCELL(m_curi, m curj);
break;
}
case(9):
case) 10):
{
GoDirection(3, m curi, m_curj);
CHECKCELL(m_curi, m_curj);
GetDirection) 1, m_curi, m cuij, m_ti, m_tj);
CHECKCELL(m_ti, m_tj);
GoDirection(2, m_curi, m curj);
CHECKCELL(m_curi, m_curj);
break;
}
#ifdef FAILSAFE MODE
#endif
m_intpoint += m_hex_a3;
49
}
#ifdef FAILSAFEMODE
//if it reached this, there's trouble
CString temp;
temp.Format("Failed to reach target for (%d, %d) to (%d, %d)", bi, bj, ei, ej);
AfxMessageBox(temp);
#endif
}
else
{
//project LOS ray up to first intersection onto the primary axis
m_usedist = mhexh / mrelativeSlope;
//the number of steps out the axis to make
m_distx = m_usedist / mhexal5 3.0;
//this starts at 2 so that last step can by independently
//analyzed
mdistance = 2.0;
niStepper = 0;
mtype = TRUE;
#ifdef FAILSAFEMODE
mtotallen = 0;
//the hypotenuse of the triangle, distance along the LOS ray
//between intersection points
m_usedlen = sqrt(m_usedist m_usedist + m_hex_h m_hex_h);
m_farout = mrayLength + musedlen;
while(m_totallen < m farout)
#else
while(l)
#endif
{
m_distance += m_distx;
for(;m_stepper
{
if(m_subtype = m_type)
{
switch( redirection)
case(0):
case(6):
case(9):
case(10):
{
GoDirection(2, mcuri, mcurj);
break;
}
case(3):
case(4):
GoDirection(l, m curi, m curj);
50
}
break;
else
switch( redirection)
{
case(O):
case(6):
{
GoDirection( 1, m_curi, m_curj);
break;
}
case(3):
case(4):
{
GoDirection(0, m_curi, mcurj);
break;
}
case(9):
case(10):
GoDirection(3, m curi, m_curj);
break;
m_subtype = !m_subtype;
CHECKCELL(m_curi, m curj);
if(((m_direction = 3)  (m_direction == 6) 
(indirection = 9)) ?
(m_type == m subtype): (mtype != m_subtype))
{
// add 1 back on to see if you fall on the far side of
// the parallel cell boundary
mdval = m_distance + 1.0;
if(m_stepper < m_dval)
{
if(m_subtype = mtype)
{
switch( indirection)
{
case(0):
case(3):
case(4):
{
51
GoDirection(l, mcuri, mcurj);
break;
case(6):
case(9):
case(10):
GoDirection(2, m_curi, m_curj);
break;
}
}
else
switch(m_direction)
{
case(O):
{
GoDirection(l, m_curi, m_curj);
break;
}
case(3):
case(4):
{
GoDirection(0, m_curi, m curj);
break;
}
case(6):
{
GoDirection(2, m_curi, m curj);
break;
}
case(9):
case(10):
{
GoDirection(3, mcuri, m_curj);
break;
}
m_subtype = !m_subtype;
CHECKCELL(m_curi, m curj);
mstepper += 3;
}
else // have to make the perpedicular step
{
switch(mdirection)
{
case(O):
52
{
GoDirection(0, mcuri, m_curj);
break;
}
case(3):
{
GoDirection(2, m curi, m_curj);
break;
}
case(4):
{
GoDirection(5, m curi, mcurj);
break;
}
case(6):
{
GoDirection(3, m_curi, m_curj);
break;
}
case(9):
{
GoDirection(l, m_curi, m_curj);
break;
}
case(10):
{
GoDirection(4, m_curi, m_curj);
break;
}
CHECKCELL(m_curi, m curj);
}
}
#ifdef FAILS AFEMODE
mtotallen += musedlen;
#endif
mtype = !m_type;
}
#ifdef FAILS AFEMODE
//if it reached this, there's trouble
CString temp;
temp.Format("Failed to reach target for (%d, %d) to (%d, %d)", bi, bj, ei, ej);
AfxMessageBox(temp);
#endif
}
break;
}
}
53
#ifdef EXPORTMATLABVISIBILITYCHECK
mout.CloseFile();
#endif
}
54
BIBLIOGRAPHY
[1] Yap, P. 2002. GridBased PathFinding. In Proceedings of the 15th Conference of
the Canadian Society For Computational Studies of intelligence on Advances in
Artificial intelligence (May 27 29, 2002). R. Cohen and B. Spencer, Eds. Lecture
Notes In Computer Science, vol. 2338. SpringerVerlag, London, 4455.
[2] Wesley E. Snyder, Hairong Qi, William Sander, "A coordinate system for
hexagonal pixels" Proc. of SPIE, Vol. 3661, pp. 716727, 1999.
[3] Luczak, E. and Rosenfeld, A. 1976. Distance on a Hexagonal Grid. IEEE Trans.
Comput. 25, 5 (May. 1976), 532533. DOI=
http://dx.doi.Org/l 0.1109/TC. 1976.1674642
[4] Verbrugge, C. (1997) Hex Grids. Unpublished manuscript.
[5] Wiithrich, C. A. and Stucki, P. 1991. An algorithmic comparison between square
and hexagonalbased grids. CVGIP: Graph. Models Image Process. 53, 4 (Jun. 1991),
324339. DOI= http://dx.doi.org/10.1016/10499652r9n90Q36J
[6] Umanskiy, Oleg. LG Hypergames for RealWorld Defense Systems Masters
thesis, University of Colorado at Denver, 2001
[7] Brimkov, V. E. and Bameva, R. P. 2005. Analytical Honeycomb Geometry for
Raster and Volume Graphics. Comput. J. 48, 2 (Mar. 2005), 180199. DOI=
http://dx.doi.org/10.1093/cominl/bxh07 5
55
