PROGRAMMING  TIPS

An April Fool's Hoax

Some subtle ideas from computational geometry can be used to form the adjacency graph of a map, which in turn can be used to resolve a famous four-coloring hoax.

by Stan Wagon

In my column four issues ago, Joan Hutchinson and I discussed an algorithm that could be used to 4-color planar graphs. But coloring maps is more difficult since the map must first be turned into a graph, and if some countries are nonconvex then that could be tricky. In this column I discuss how the triangulation of a country can be used as the basis of a scheme to create a planar graph corresponding to any map. My motivation is a map presented as a hoax by Martin Gardner in the April 1, 1975 issue of Scientific American: he asserted that the map required five colors. By the end of this column we will see how Mathematica can be used to color maps with nonconvex countries. We will conclude by taking a trip, via the dual map construction, to the planet Bizarro.

.  Triangulating a Polygon

The following seems obvious but, after digesting the statement, you might pause and think about how it would be proved: Every polygon has two vertices such that the line connecting them lies inside the polygon. To highlight the difficulty of relying on intution, I'll point out that the analogous statement for polyhedra in 3-space is false! There are polyhedra that have no interior diagonal. But for plane polygons the result is true and is quite important. For example, we can use it to prove that every polygon can be triangulated by repeatedly drawing lines from vertices to other vertices. Just choose an interior diagonal and use recursion on the two polygons that it leads to (they each have fewer sides than the original). This leads to a workable algorithm that we will implement.

It is very useful to have on hand a library of routines to perform simple geometrical constructions. We will start with a few that will be needed for triangulation. See [O'Rourke 19xx] for proofs of various facts we will assume.

A key notion in working with polygons is that of signed area. If the points are given in the positive (counterclockwise) direction, the signed area is just the area; otherwise it is the negative of the area. The signed area of a triangle is given by one-half of the determinant of (({{(Subscript[x, 1]), (Subscript[y, 1]), 1}, {(Subscript[x, 2]), (Subscript[y, 2]), 1}, {(Subscript[x, 3]), (Subscript[y, 3]), 1}})). For speed we implement this determinant as a compiled function and use it in our SignedArea function.

\(SACom = 
      Compile[{x1, y1, x2, y2, x3, y3}, 
        0.5 \((\(-x2\)\ y1 + x3\ y1 + x1\ y2 - x3\ y2 - x1\ y3 + x2\ y3)\)]; 
    

SignedArea[{p_, \ q_, \ r_}] := SACom@@Flatten[{p, q, r}]\)Area[p_] := Abs[SignedArea[p]]

An important application is the orientation of a polygon, which we take to be +1 if the vertices are given in counterclockwise order, –1 otherwise. We could simply look at the sign of the signed area, but there is a slightly trickier way. For a triangle, we do follow the obvious approach.

Orientation[{p_,  q_,  r_}] := Sign[SignedArea[{p, q, r}]]

But the orientation of a polygon is determined by the orientation at a convex vertex (a vertex for which the interior angle is less than 180˚). We can find a convex vertex by using the first vertex in the list returned by Sort. While it seems inefficient to sort the whole list, the fact that Sort is built-in makes this adequately fast. We work here with indices rather than the points themselves, so that ConvexVertex returns the position of a convex vertex in a polygon.

ConvexVertex[p_] := Position[p, First[Sort[p]]][[1, 1]]

Orientation[poly_] := 
      Orientation[
          poly[[ConvexVertex[poly] + {-1, 0, 1} /. {
                0 -> Length[poly], Length[poly] + 1 -> 1}
                ]]] /; Length[poly] > 3

A less obvious use of orientation is to tell whether a point lies to the left of the directed infinite line determined by two other points. LeftOf returns True if the third point is strictly left of the line determined by the first two.

LeftOf[p_, q_, r_] := (Orientation[{p, q, r}] == 1)

With these few simple routines in hand, we can attack the triangulation problem. Any n-gon can be triangulated by drawing a certain number — n – 2 — of internal diagonals. We will want our routine to return the indices of the triangles in the triangulation. The proof/algorithm can be by recursion: we simply choose an interior diagonal and use it to break the polygon into two pieces, which can be triangulated recursively.  Here's a proof of the fact that every polygon has an interior diagonal. First choose a convex vertex v, and let u and w be the neighboring vertices. Consider the triangle formed by u, v, and w. If no other vertex lies inside or on this triangle, then uw is the desired interior diagonal. Otherwise, choose from those vertices lying inside the triangle the one that is farthest from the segment uw (see Figure 1). It is easy to see that the segment determined by v and this optimal vertex is the desired interior diagonal.

Code for Figure 1

Show[Graphics[{AbsolutePointSize[5],
Line[{aa,w,v,u,bb}={{3,-2},{2,0}, {0,-2}, {-2,0.2},{-3,-.31}}],
Line[pts = {{2.8, -.2},{1.8, 1.35}, {1.22678, 1.07989}, 
{1.6, 0.55}, {0.860153, 0.501752}, 
{1.22678, 0.0223221}, {0.733245, -0.358402}, 
{0.507631, -0.104586}, 
{0.0705044, -0.259695}, {0.282017, -0.49941}, 
{0.549934, -0.584016}, {0.16921, -0.950638}, 
{-0.930658, -0.0481823}, {0.112807, 0.233835}, 
{-1.42419, 1.07989}, {-3, -.2}}], Point /@ pts,
 Point /@ {u,v,w},
 AbsoluteDashing[{3,3}], Line[{v, pts[[-5]]}],
 Line[{u,w}],Line[pts[[{2,4}]]],
 Text["v", v-{.15,.15}],
  Text["u", u-{.0,.15}],
   Text["w", w-{0,.22}],  
   Text["x", pts[[-5]]-{.15,.15}],
   Text[Subscript["v",1], pts[[3]]-{.15,.15}],
  Text[Subscript["u",1], pts[[2]]-{.15,-.13}],
   Text[Subscript["w",1], pts[[4]]-{-.18,.15}] 
}], AspectRatio->Automatic, 
DefaultFont -> {"Times", 9},
Frame->True,FrameTicks->None,
PlotRange -> {{-2.5, 2.5}, {-2.3,1.6}}];

We can simply follow this proof to get an algorithm. First we get an (interior) polygon diagonal, which we return as a sorted pair. Inside works only for positively oriented triangles, so care is necessary. We use a sort based on Area (since altitude is proportional to area for a fixed base) to find the extreme vertex in the case that Δuvw contains other vertices.

\(Inside[{p_, q_, r_}, pt_] := 
			
      \((¬ LeftOf[q, p, pt] ∧ 
          ¬ LeftOf[r, q, pt] ∧ 
          ¬ LeftOf[p, r, pt])\)
		\)\(PolygonDiagonal[poly_] := 
      Module[
			{v = ConvexVertex[poly], n\  = \ Length[poly], 
          invertices, u, w}, {u, w}\  = \ 
          \((v\  + \ {\(-1\), 1}\  /. {0 -> n, n + 1 -> 1})\); 

        			invertices = 
          Select[Delete[poly, {{u}, {v}, {w}}], 
					\ \ \ \ 
            Inside[poly[[{u, v, w}]], #]
                \ &]; 
			
        Sort[If[invertices == {}, {u, w}, 

            					{
              \(Position[poly, 
                  Last[Sort[invertices, 
                      Area[
								{
                              poly[[u
                                ]], #1, 
                              poly[[w
                                ]]}] ≤ 
                          Area[{poly[[u
                                ]], #2, 
                              poly[[w
                                ]]}]&]]]
                  \)[[1, 1]], v}]]]\)

And now we can use recursion to triangulate an arbitrary polygon. Because we wish to work with indices, we define the recursive part of the routine to take two arguments: the main polygon (which won't change), and the list of indices defining the smaller polygon. This case takes two arguments, and so is distinguished from the main case, which uses the same name, Triangulate, but takes only one argument. It seems reasonable to return the triangles in the same orientation as the original polygon, and this is what flip does. The main recursive case (inds as an argument) needs no condition (/; Length[inds] > 3); it will be checked last because it is more general than the other case ({i1_, i2_, i3_}). This sort of thing can be checked with ??Triangulate). Sometimes recursive code is wasteful of memory and slower than a nonrecursive implementation. The reader might try a different implementation of this triangulation algorithm and see if it is faster. Indeed, there are other algorithms that are faster (see [O'Rourke 199x]). In any event, it is hard to beat the conciseness of the recursive code: It finds a diagonal, splits the polygon, triangulates the two pieces, and combines the triangles so obtained. Because we use the original indices, no work is needed when the triangles are combined.

\(Triangulate[polygon_List] := 
      \((flip = If[Orientation[polygon] == \(-1\), Reverse, Identity]; 
        flip/@Triangulate[flip[polygon], Range[Length[polygon]]])\)
\)\(Triangulate[poly_, inds_] := 
      Module[{d = 
            PolygonDiagonal[
              poly[[inds]]]}, 
        Join[Triangulate[poly, Take[inds, d]], 
          Triangulate[poly, Drop[inds, d + {1, \(-1\)}]]]]
\)Triangulate[_, {i_, j_, k_}] := {{i, j, k}}

Here is a nonconvex polygon that we will use as an example.

poly = {{2.03996, -1.62395}, {2.10014, 0.497341}, {0.174424, 
          0.91859}, {1.79924, -0.194712}, {0.746119, -0.164622}, {
          1.61871, -0.886764}, {-0.577808, -0.811541}, {
          -0.863655, 0.482296}, {-1.75129, 0.542475}, {-1.66102, 
          -1.2027}, {-0.653031, -1.62395}, {-1.88669, 
          -1.95493}, {1.46826, -2.12042}, {0.866476, -1.60891}, {
          2.22049, -2.0452}, {2.43112, -1.45846}};

Triangulate[poly]

\(\(Show[
      Graphics[{
		
          \({{GrayLevel[\ \((1 + Random[])\)/2], 
                  Polygon[poly[[#1]]]}, 
                
				
                Line[poly
                    [[Append[#1, First[#1]]
                      ]]]}&\)/@
				\ \ \ \ \ 
            Triangulate[poly], 

          		{AbsoluteThickness[2], 
					
            Line[Append[poly, \ First[poly]]]}}], 
	\ 
      AspectRatio -> Automatic]; \)\)

There are several improvements one could make to this triangulation algorithm. But the simple approach taken here is adequate for small polygons.

.  An Application

We now need the GraphColoring package, so we remove the duplicated names and load it.

\(Remove[ConvexVertex, \ SignedArea, \ Orientation, \ Triangulate]; 

    Needs["GraphColoring`"]\)

A triangulation of a polygon can be used as a sort of coordinate system to help one wind one's way through a complicated polygon. Suppose P is a polygon and c a point inside P. For the application we have in mind we seek disjoint, piecewise-linear paths from c (the "capital") to the midpoints of all the edges of P. We can use the abstract graph (in fact, it is a tree) formed by the triangles to find a sequence of triangles that will take us from the capital to each midpoint. Then we can carefully cross the sides of the triangles so that the paths are disjoint. Figure xx shows the labeled triangles from our example.

Code for Figure

\(tris\  = \ 	Triangulate[poly]; 


    Show[Graphics[{
          MapIndexed[{{GrayLevel[\ \((1 + Random[])\)/2], 
                  Polygon[poly[[#1]]]}, 
                
					Text[#2[\([1]\)], \ Centroid[poly[\([#1]\)]]], 

                					\ 
                Line[poly
                    [[Append[#1, First[#1]]
                      ]]]}&, 
				tris], 

          			{AbsolutePointSize[5], \ Point\ /@\ poly}, 

          		{AbsoluteThickness[2], 
					
            Line[Append[poly, \ First[poly]]]}}], 
	
      DefaultFont -> {"Times", 8}, AspectRatio -> Automatic]; \)

Now we can form the graph — in fact, it is a tree — corresponding to the triangles and their adjacencies. Then the ShortestPath routine from the Combinatorica package tells us how to get from any triangle to any other.

\(tris\  = \ Triangulate[poly]; 
			
\ 
    adj[i_, j_]\  := \ 
      If[Length[Intersection[tris[\([i]\)], tris[\([j]\)]]] == 2, 1, 0]
	
    \)\(tree\  = \ 
      Graph[Array[adj, {Length[tris], \ Length[tris]}], 
				
        \(Centroid[poly[\([#]\)]]&\)\ /@\ tris]; 


    ShortestPath[tree, 7, 1]\)

Here is a picture of the triangle tree.

ShowLabeledGraph[tree, LabelOffset -> 0, 
      VertexColor -> White, VertexSize -> AbsolutePointSize[11.5], 
      Background -> GrayLevel[0.5], 
      DefaultFont -> {"Times", 8}, Frame -> True, 
      FrameTicks -> None];

Now we come to the geometric crux of the dual graph construction: the PathDecomposition routine. This function takes a polygon and returns a list of disjoint paths from a suitably chosen capital to the midpoints of the edges. The first entry in the output is the capital and the second is the list of paths, each one excluding the capital and given in the natural order of the edges.

pd = PathDecomposition[poly];

Here is a view of a typical path decomposition.

Show[
      Graphics[{
          ({{GrayLevel[0.3 Random[] + 0.7], 
                    Polygon[poly[[#]]]}, 
                  GrayLevel[0.5], 
                  Line[poly
                      [[Append[#, First[#]]
                        ]]]}&)/@tris, {
            AbsolutePointSize[5], 
            Point[pd[[1]]]}, {
            AbsoluteThickness[1.5], Line[Append[poly, First[poly]]]}, {
            AbsolutePointSize[4.5], 
            Point/@((poly + RotateLeft[poly])/2)}, 
          (Line[Prepend[#, pd[[1]]]]&)
              /@pd[[2]]}], 
      AspectRatio -> Automatic];

In practice the country-to-country adjacencies in a map will not involve all the edges. So the PathDecomposition function accepts as a second argument a list of edges to be reached. Moreover it tries all the centroids to see if oone of them has the property that the edge-centers can all be reached by straight lines. If so, that one along with the stright lines are used; if not, the capital is taken to be the centroid of the largest triangle and piecewise-linear paths are used.

borders = {6, 9, 10, 12, 15}; pd = PathDecomposition[poly, borders, ShowTriangles -> True];

Code for Figure

borders = {6, 9, 10, 12, 15}; pd = PathDecomposition[poly, borders,  ShowTriangles -> True];

\(\(Show[
      Graphics[{

          			{AbsolutePointSize[5], 
            Point[pd[[1]]]}, {
            AbsoluteThickness[1.5], Line[Append[poly, First[poly]]]}, {
            AbsolutePointSize[4.5], 
            Point/@\((\((poly + RotateLeft[poly])\)/2)
                  \)[[borders]]}, 
          \((Line[Prepend[#, pd[[1]]]]&)
              \)\ /@\ pd[[2]]}], 
      AspectRatio -> Automatic]; \)\)

borders = {1,2,6,9,10,12};
pd = PathDecomposition[poly,borders];

Show[Graphics[{
  ({{GrayLevel[0.3 Random[]+0.7],
          Polygon[poly[[#]]]},
        GrayLevel[0.5],
        Line[poly
            [[Append[#,First[#]]]]]}&)/@
  tris,
  {AbsolutePointSize[5], Point[pd[[1]]]},
  {AbsoluteThickness[1.5], Line[Append[poly, First[poly]]]},
  {AbsolutePointSize[4.5],Point /@ ((poly + RotateLeft[poly])/2)[[borders]]},
 Line[Prepend[#, pd[[1]]]]& /@ pd[[2]]

}], AspectRatio->Automatic];

If one subdivides the edges so that each edge is really, say, 10 pieces in a row, then the preceding routine will yield paths that come closer to filling the polygon. This can be viewed as an approximation to a homeomorphism from the polygon and its interior to a closed disk, and I imagine that this might be useful in other geometric situations.

.  The Graph of a Map

The PathDecomposition function returns lists of points corresponding to the paths. So it is a simple programming problem to take this data and use it to go from a map to the graph that codes the adjacencies. The definition of the true dual graph is more subtle, and is addressed in section xx. For coloring problems, we need only a graph, which we call the AdjacencyGraph, that captures the adjacency information from a map. For this we will want to place a point inside each country and connect two such points if the countries are adjacent. But straight-line edges are inadequate, for it is critical that the graph we obtain in this way be planar; moreover, we will require that the edges in the graph stay inside the corresponding adjacent countries.

Thus a key point in building the graph is allowing piecewise lienar edges. For this we introduce PlanarGraphPL objects, which are similar to Graph objects in Combinatorica in that they have two arguments, an adjacency matrix and a list of points. However, the matrix can have a 0, a 1, or a list of points in each position. A list is interpreted as an adjacency, with the points of the list defining the path between the vertices. A 1 gets interpreted as a straight line, as in Graph.

As an example, we will use some eastern states of the U.S., together with Lake Michigan considered as a separate entity. The inclusion of the lake is inessential for the generation of a dual, but is critical in order that Euler's formula hold for the map, and that, in turn, is important in four-coloring algorithms. The omission of Lake Michigan from MapOfUSA was an oversight in the presentation in [Hutchinson and Wagon, 1997].

Sort[
      (Position[StateNames, #][[1, 1
              ]]&)/@{"MD", "VA", "PA", 
          "NY", "NJ", "MI", "MN", "IL", "IN", 
          "OH", "WV", "DE", "IA", "WI", "NC", 
          "KY", "LakeMichigan"}]

map = 
      InducedSubmap[
        MapOfUSA, {7, 10, 12, 13, 15, 18, 20, 21, 25, 29, 32, 33, 36, 43, 46, 
          47, 49}];

ShowGraph[AdjacencyGraph[map], 
      VertexSize -> AbsolutePointSize[4]];

But we can straighten many of the edges, provided we check that the new edges stay within the two adjacent countries and do not cross any other edges. This checking takes a while, but the result is a much more pleasing graph.This feature is activated via the StraightenEdges option to AdjacencyGraph and ShowMap.

ShowMap[map, ShowAdjacencyGraph -> True, 
      StraightenEdges -> True, BorderStyle -> GrayLevel[0.5], 
      CountryColors -> Array[GrayLevel[Random[] 0.45 + 0.5]&, 17], 
      BorderPoints -> None, VertexSize -> AbsolutePointSize[4], 
      Background -> None, Frame -> False];

Finally, here is the April Fool's map of Martin Gardner, with some variations from the original to make the coloring harder. The aspiring map-colorer might try to 4-color this map by hand. We also add the exterior face, which is a region with one edge doing double duty.

ShowMap[map = AddExteriorFace[MapAprilFoolsHarder], 
      CountryColors -> White, Background -> White, 
      BorderPoints -> False, Frame -> False, 
      AspectRatio -> 1/GoldenRatio];

Figure  shows the dual graph for this map, in both raw and straightened form. In this case the straightened graph is a vast improvement; however, from the point of view of a coloring algorithm, the raw graph contains enough information and that is the one that is used when a 4-coloring of the map is requested.

ShowGraph[AdjacencyGraph[map]];

ShowMap[map, ShowAdjacencyGraph -> True, 
      RepeatedStraightenings -> True];

Figure 10(a)

Figure 10(b)

Code For Figure (b)

\(\(ShowMap[map, ShowAdjacencyGraph -> True, 
      RepeatedStraightenings -> True, 
      CountryColors -> 
        Array[GrayLevel[Random[]\ 0.45 + 0.5]&, 
          Length[map[[1]]]], 
      Background -> White, 
	BorderStyle -> GrayLevel[0.5], 
      BorderPoints -> False, Frame -> False]; \)\)

.  Coloring the April Fool’s Map

One can ask ShowMap to 4-color a map. If we do that for the April Fool's map, then the dual graph construction and its coloring are transparent to the user. The algorithm used to 4-color a graph is based on Kempe chains, as described in my earlier column [Wagon 1997].

ShowMap[map, FourColorCountries -> True];

Code for Figure

\(\(ShowMap[map, FourColorCountries -> True, 

      CountryColors -> GrayLevel/@{0.9, \ 0.65, \ 0.4, 0}, 
	
      Background -> None, Frame -> None, BorderPoints -> None, 
      
AspectRatio -> 1/GoldenRatio]; \)\)

The algorithm used to 4-color a graph is based on Kempe chains, as described in my earlier column [Hutchinson and Wagon, 1997]. But I have made some important improvements. The algorithm as described in my previous column was based on Kempe's strategy of using Kempe color-switches, and this could fail on certain graphs, notably the Errera graph. But I. Kittell, in 1935, figured out a workaround: he allowed Kempe switches of adjacent vertices in the neighbor ring and this can often untangle a situation where Kempe chains get tangled. It turns out that allowing random Kempe switches, including those of the Kittell type, is a good strategy — indeed, perhaps it never fails (see [Hutchinson and Wagon, 1998]) — and my 4-coloring algorithm now invokes Kittell's strategy whenever Kempe's original method fails.

.  Continents on the Planet Bizarro

warning: can take an hour or more for some of these.

The adjacency graph is but a pale imitation of the true dual graph. The package has the capability of getting the true dual of a map (or graph), and so one can illustrate the fact that the dual map of the dual map of a map is the original map. To get the dual map of a map one must, given a map in the plane, form a graph by placing a vertex in each country and connecting two vertices by several edges, one for every edge segment of the map that witnesses their adjacency. This requires setting up MultiGraph objects. The multigraph must then be turned into a graph by adding vertices at all the bends in all the edges. But these vertices are not really there, which adds a complication (actually, a simplification) later; see step 5 below. Briefly, here is how to see the Bizarro world of western Europe that arises as the double dual.

1. Let map1 be the map of western Europe with exterior face added; use AddExteriorFace.

2. Let graph1 be dual multigraph of map1, with straightened edges; use DualMultiGraph.

3. Let graph1a be graph1 with added vertices at all the bends; use MultiGraphToGraph.

4. Turn graph1a into a PlanarMap object, map2, by extracting the face data and adding an exterior face. This is done by PlanarGraphToPlanarMap, which uses some functions from the Combinatorica package (FindCycle) and some geometric programming. map2 is the true dual map of map1; the number of faces in map2 equals the number of vertices in map 1, and vice versa.

5. Now get the adjacency graph of map1 (with exterior face added) as described in section 3; this places only one piecewise-linear edge between capitals of adjacent countries, and is carried out by AdjacencyGraph, again with the RepeatedStraightenings option turned on. This use of the adjacency graph, as opposed to the dual multigraph, is how the question of the extraneous vertices added in step 3 is dealt with. The faces in this graph form a map isomorphic to Europe! Use PlanarGraphPLToGraph to get a graph with vertices at all bends amd then PlanarGraphToPlanarMap to turn it into a PlanarMap object. This final map is isomorphic to the original map! The process is not really mysterious once when sees that the dual construction turns vertices into faces and faces into vertices; hence the double dual returns the original.

map1  =  MapOfWesternEurope;

ShowMap[AddExteriorFace[map1], Background -> None, 
      Frame -> False];

Code for Figure

\(\(ShowMap[AddExteriorFace[map1], Background -> None, 
	
      Frame -> False, 
      CountryColors -> 
        GrayLevel\ /@
          \ { .4,  .2,  .7,  .5,  .2,  .9,  .6,  .2,  .4, 1,  .2,  .8}, 
      BorderPointStyle -> AbsolutePointSize[2], 

      BorderStyle -> AbsoluteThickness[0.4]]; \)\)

Figure . Western Europe with an exterior face added. This map, viewed as a map on the sphere (so that the 4 extremal vertices and corresponding edges do not count) has 12 faces, 61 vertices, and 71 edges.

graph1 = DualMultiGraph[map1, RepeatedStraightenings -> True];

graph1a = MultiGraphToGraph[graph1];

\(\(ShowGraph[graph1a, VertexSize -> AbsolutePointSize[3], 

      Epilog\  -> 
        \ {AbsolutePointSize[6.5], \ Point\ /@\ graph1[\([2]\)]}]; \)\)

map2 = PlanarGraphToPlanarMap[graph1a];

\(\(ShowMap[map2, FourColorCountries -> True, 

      Epilog\  -> 
        \ {AbsolutePointSize[6.5], \ Point\ /@\ graph1[\([2]\)]}]; \)\)

graph2 = 
      AdjacencyGraph[AddExteriorFace[map2], 
        RepeatedStraightenings -> True];

map3 = PlanarGraphToPlanarMap[PlanarGraphPLToGraph[graph2]];

ShowMap[map3];

Code for Figure (a)

\(\(ShowGraph[graph1a, VertexSize -> AbsolutePointSize[3], 

      Epilog\  -> 
        \ {AbsolutePointSize[6.5], \ Point\ /@\ graph1[\([2]\)]}]; \)\)

Code for Figure (b)

\(rad\  = \ 25; 


    ShowMap[map2, Background -> None, Frame -> None, 
      FourColorCountries -> True, 
      CountryColors -> GrayLevel/@{0.13, 0.32, 0.65, 1}, 
      BorderPoints -> None, 

      Epilog\  -> 
        \ {AbsolutePointSize[6.5], \ White, 
          \(Disk[#, \ rad]&\)\ \ /@\ graph1[\([2]\)], 
		Black, 
          AbsoluteThickness[0.5], \ \ 
          \(Circle[#, \ rad]&\)\ /@\ graph1[\([2]\)]}]; \)

Figure 12(a)

Figure 12(b)

Code for Figure

ShowMap[map3, BorderPoints->None, 
CountryColors->GrayLevel /@ {0.2, 0.4, 0.7, 0.2, 0.9, 0.4, 1, 0.5, 0.2, 0.2, 0.6},
Background -> GrayLevel[0.8]];

What's the point? Well, it is certainly cool to see all this complicated programming come together and get a distorted version of Europe. These tools surely have some educational value, both in computational geometry, where they illustrate the power of triangulation, and in graph theory, where the images can help students learn about duals and double duals.

References

[O'Rourke, 1993]   J. O'Rourke, Computational Geometry in C, Cambridge Univ. Pr., New York.

[Hutchinson and Wagon, 1997]  Joan Hutchinson and Stan Wagon, The four-color theorem, Mathematica in Education and Research 6:1,
42–51.

[Hutchinson and Wagon, 1998] Kempe revisited, American Mathematical Monthly, to appear.

ABOUT THE AUTHOR

Stan Wagon uses Mathematica extensively in his research, teaching, and exposition and is especially appreciative of the new ways in which Mathematica allows us to look at mathematical objects. He is author or coauthor of several books, including Mathematica in Action, Animating Calculus, VisualDSolve, and Which Way Did the Bicycle Go?. He also teaches a course at Rocky Mountain Mathematica, a high-altitude workshop held in July in the mountains of Colorado.

Stan Wagon
Department of Mathematics and Computer Science
Macalester College
St. Paul, MN 55105
wagon@macalester.edu

Electronic Supplement:  GraphColoring.m