Tessellating shapes on top of Bing Maps in a WP7 app

Before I complete the 'How many pins can Bing Maps handle in a WP7 app...' set of posts. I want to take one last diversion and show how you can tessellate shapes on top of the Bing Maps control in a WP7 app. This builds on top of a previous post about drawing shapes on top of Bing Maps. If you're unsure what tessellation is, then the following definition should give you a good idea:

A tessellation is created when a shape is repeated over and over again covering a plane without any gaps or overlaps.

Only three regular polygons tessellate in the Euclidean plane: triangles, squares or hexagons. Shown below are some graphical examples of how these shapes tessellate:

I'm only going to show how you can tessellate these shapes on top of the map control. The effect will be rendered as a set of MapPolygon class instances. You could achieve a similar effect with a set of MapPolyline class instances and it would be probably a lot easier but, and it is an important 'but', if you use a MapPolyline you won't be able to use the polygon colour fill functionality. The MapPolygon instances are render as a layer on top of the map control, the control provides this functionality out of the box, this means the XAML is very clean and easy to understand:

As you can see the polygons are rendered as part of the MapItemsControl instance, the item source is bound to the Polygons property on the ViewModel. The item template is shown below and you can see this is also very simple. I've also included the brush definitions so you can see where the polygon colours are defined:

Before I dive into the code used to generate the polygons lets see the effect I'm going to generate, shown below are a set of tessellated shapes at different sizes. The density of polygons on the screen is dependent on the zoom level of the map control - here I'm using a default zoom level of 13 when the app starts, obviously you can zoom in and out as required and the polygons will automatically re size. This functionality is provided by the parent MapLayer class instance.

Here I'm using a set of pre-defined shapes, these are defined in the ViewModel constructor which is shown below. You can see each Shape class instance has a Name property and 2 function delegate properties - PolygonFunc and AdjacentPolygonFunc. This is the same style I used in the previous post.

public MapViewModel(ILog log)
{
    this.log = log;

    this.polygons = new ObservableCollection<Polygon>();
    this.datum = new GeoCoordinate(51.54590803376571, -0.10334014892578125);
    this.Centre = new GeoCoordinate(51.561811605968394, -0.0883626937866211);

    this.Zoom = 13;
            
    this.shapes = new ObservableCollection<Shape>
    {
        new Shape { Name = "No Shape", PolygonFunc = centre => new Polygon(), AdjacentPolygonFunc = centre => new List<Polygon>() },
        new Shape { Name = "Square (100 m)", PolygonFunc = MapFuncs.Square(100), AdjacentPolygonFunc = MapFuncs.AdjacentSquares(100) },
        new Shape { Name = "Square (200 m)", PolygonFunc = MapFuncs.Square(200), AdjacentPolygonFunc = MapFuncs.AdjacentSquares(200) },
        new Shape { Name = "Square (500 m)", PolygonFunc = MapFuncs.Square(500), AdjacentPolygonFunc = MapFuncs.AdjacentSquares(500) },
        new Shape { Name = "Square (1 km)", PolygonFunc = MapFuncs.Square(1000), AdjacentPolygonFunc = MapFuncs.AdjacentSquares(1000) },
        new Shape { Name = "Square (1 mile)", PolygonFunc = MapFuncs.Square(1609.347), AdjacentPolygonFunc = MapFuncs.AdjacentSquares(1609.347) },
        new Shape { Name = "Triangle (3 sides, 500 m)", PolygonFunc = MapFuncs.Triangle(500), AdjacentPolygonFunc = MapFuncs.AdjacentTriangles(500) },
        new Shape { Name = "Triangle (3 sides, 1 km)", PolygonFunc = MapFuncs.Triangle(1000), AdjacentPolygonFunc = MapFuncs.AdjacentTriangles(1000) },
        new Shape { Name = "Triangle (3 sides, 1 mile)", PolygonFunc = MapFuncs.Triangle(1609.347), AdjacentPolygonFunc = MapFuncs.AdjacentTriangles(1609.347) },
        new Shape { Name = "Hexagon (6 sides, 500 m)", PolygonFunc = MapFuncs.Hexagon(500), AdjacentPolygonFunc = MapFuncs.AdjacentHexagons(500) },
        new Shape { Name = "Hexagon (6 sides, 1 km)", PolygonFunc = MapFuncs.Hexagon(1000), AdjacentPolygonFunc = MapFuncs.AdjacentHexagons(1000) },
        new Shape { Name = "Hexagon (6 sides, 1 mile)", PolygonFunc = MapFuncs.Hexagon(1609.347), AdjacentPolygonFunc = MapFuncs.AdjacentHexagons(1609.347) },
    };
    this.SelectedShape = this.shapes[5];
}

Hopefully the delegate method names are obvious, PolygonFunc calculates the geo-location points for the shape and the AdjacentPolygonFunc. The shapes are divided into 2 distinct groupings - squares and everything else. Squares are easy to calculate because you only have to calculate geo-locations for 4 points using constant bearings  - 0, 90, 180 & 270. Everything else is where it starts to get complicated, you have to use a circle and calculate the position of each point on the circumference. The diameter for the circle comes from the defined size of the shape.

Shown below is the PolygonFunc method for creating a hexagon, as you can see the hexagon is generated using the Polygon method we also set the Height & Width of the generated polygon.

public static Func<GeoCoordinate, Polygon> Hexagon(double diameter)
{
    Func<GeoCoordinate, Polygon> func = location =>
    {
        var hexagon = Polygon(location, 6, diameter, 0, 0, 0);

        var north = hexagon.Locations[0];
        var south = hexagon.Locations[3];
        var east = hexagon.Locations[1];
        var west = hexagon.Locations[5];
              
        hexagon.Height = north.GetDistanceTo(south);
        hexagon.Width = east.GetDistanceTo(west);
                
        return hexagon;
    };

    return func;
}

The Polygon method is shown below, this is ultimately the same method used to create triangles as well as hexagons:

private static Polygon Polygon(GeoCoordinate location, int sides, double diameter, double offset, double offsetBearing, double startAngle)
{
    var centre = location;
    if (offset != 0)
    {
        centre = CalculateUsingHaversine(centre, offset, offsetBearing);
    }

    var polygon = new Polygon();

    var radius = diameter / 2;
    var angle = CircleInDegresses / sides;

    for (var i = 0; i < sides; i++)
    {
        var aggregatedAngle = (i * angle) + startAngle;
        polygon.Locations.Add(CalculateUsingHaversine(centre, radius, aggregatedAngle));
        polygon.Centre = centre;
        polygon.TileBoundingRectangle = Square(diameter)(centre).TileBoundingRectangle;
    }

    return polygon;
}

For this demo I'm using the Haversine method for calculating distance between two geo-locations. There are more accurate ways to calculate geo-locations, I'm aware of the error factor (0.3 %) and happy to negate this for Haversine formula. The reason is simple - in a real-world app I won't be displaying a large number of polygons over a large areas (see the 'square 100 metre' in the above screenshot) so therefore the total distances calculated would not be large and the error factor would be small. Shown below is a screenshot of a single hexagon polygon rendered onto the map control.

The purpose of the AdjacentPolygonFunc method should be obvious - to calculate the adjacent polygons for an existing polygon. The implementation of this function depends on the type of polygon, the number of adjacent polygons for a square is different to a triangle or hexagon. Shown below are the rendered adjacents; I've added the order in which the adjacent polygons are calculated. An understand of the order of calculation is important when it comes to rendering the rest of the polygons for the visible map control.

As I said the adjacent methods are dependent on the type of polygon and calculating the adjacents for a square is the easiest. Shown below is the method used for adjacent hexagons. It requires not only the ability to calculate the polygon points but the ability to calculate the points in the correct geo-locations - the points have to be 'massaged' into the correction geo-locations.

public static Func<Polygon, List<Polygon>> AdjacentHexagons(double length)
{
    Func<Polygon, List<Polygon>> func = hexagon =>
    {
        var adjacentPolygons = new List<Polygon>();
        GeoCoordinate deltaNorthWest;
        GeoCoordinate deltaNorthEast;
        GeoCoordinate deltaSouthWest;
        GeoCoordinate deltaSouthEast;

        // Calculate the horizontal distance between the edge of the hexagon and the edge of the bounding rectangle
        var eastEdge = new GeoCoordinate(hexagon.Locations.OrderByDescending(l => l.Latitude).Select(l => l.Latitude).First(),
                                            hexagon.Locations.OrderByDescending(l => l.Longitude).Select(l => l.Longitude).First());

        var horizontalDelta = eastEdge.GetDistanceTo(hexagon.TileBoundingRectangle.Northeast) * 2;

        // Calculate the vertical distance between the top of the hexagon and the top of the east edge
        var northEdge1 = hexagon.Locations.OrderByDescending(l => l.Latitude).First();
        var northEdge2 = hexagon.Locations.OrderByDescending(l => l.Longitude).ThenByDescending(l => l.Latitude).First();

        var verticalDelta = Math.Cos(DegreeToRadian(60.00)) * northEdge1.GetDistanceTo(northEdge2);

        // Calculate top centre hexagon position....
        deltaNorthWest = CalculateUsingHaversine(hexagon.TileBoundingRectangle.Northwest, verticalDelta, 180);
        deltaNorthEast = CalculateUsingHaversine(hexagon.TileBoundingRectangle.Northeast, verticalDelta, 180);

        var topCentre = AdjacentSquare(deltaNorthWest, deltaNorthEast, 0, 0, length);
        topCentre.Locations = Hexagon(length)(topCentre.Centre).Locations;

        // Calculate bottom centre hexagon position....
        deltaSouthWest = CalculateUsingHaversine(hexagon.TileBoundingRectangle.Southwest, verticalDelta, 0);
        deltaSouthEast = CalculateUsingHaversine(hexagon.TileBoundingRectangle.Southeast, verticalDelta, 0);

        var bottomCentre = AdjacentSquare(deltaSouthWest, deltaSouthEast, 180, 180, length);
        bottomCentre.Locations = Hexagon(length)(bottomCentre.Centre).Locations;

        // Calculate centre right hexagon position...
        deltaNorthEast = CalculateUsingHaversine(hexagon.TileBoundingRectangle.Northeast, horizontalDelta, 270);
        deltaSouthEast = CalculateUsingHaversine(hexagon.TileBoundingRectangle.Southeast, horizontalDelta, 270);

        var centreRight = AdjacentSquare(deltaNorthEast, deltaSouthEast, 90, 90, length);
        centreRight.Locations = Hexagon(length)(centreRight.Centre).Locations;

        // Calculate centre left hexagon position...
        deltaNorthWest = CalculateUsingHaversine(hexagon.TileBoundingRectangle.Northwest, horizontalDelta, 90);
        deltaSouthWest = CalculateUsingHaversine(hexagon.TileBoundingRectangle.Southwest, horizontalDelta, 90);

        var centreLeft = AdjacentSquare(deltaNorthWest, deltaSouthWest, 270, 270, length);
        centreLeft.Locations = Hexagon(length)(centreLeft.Centre).Locations;

        // Calculate top right hexagon position...
        deltaNorthEast = CalculateUsingHaversine(topCentre.TileBoundingRectangle.Northeast, horizontalDelta * 4.225, 270);
        deltaSouthEast = CalculateUsingHaversine(topCentre.TileBoundingRectangle.Southeast, horizontalDelta * 4.225, 270);

        var topRight = AdjacentSquare(deltaNorthEast, deltaSouthEast, 90, 90, length);
        topRight.Locations = Hexagon(length)(topRight.Centre).Locations;

        // Calculate top left hexagon position...
        deltaNorthWest = CalculateUsingHaversine(topCentre.TileBoundingRectangle.Northwest, horizontalDelta * 4.225, 90);
        deltaSouthWest = CalculateUsingHaversine(topCentre.TileBoundingRectangle.Southwest, horizontalDelta * 4.225, 90);

        var topLeft = AdjacentSquare(deltaNorthWest, deltaSouthWest, 270, 270, length);
        topLeft.Locations = Hexagon(length)(topLeft.Centre).Locations;

        // Calculate bottom right hexagon position...
        deltaNorthEast = CalculateUsingHaversine(bottomCentre.TileBoundingRectangle.Northeast, horizontalDelta * 4.225, 270);
        deltaSouthEast = CalculateUsingHaversine(bottomCentre.TileBoundingRectangle.Southeast, horizontalDelta * 4.225, 270);

        var bottomRight = AdjacentSquare(deltaNorthEast, deltaSouthEast, 90, 90, length);
        bottomRight.Locations = Hexagon(length)(bottomRight.Centre).Locations;

        // Calculate bottom left hexagon position...
        deltaNorthWest = CalculateUsingHaversine(bottomCentre.TileBoundingRectangle.Northwest, horizontalDelta * 4.225, 90);
        deltaSouthWest = CalculateUsingHaversine(bottomCentre.TileBoundingRectangle.Southwest, horizontalDelta * 4.225, 90);

        var bottomLeft = AdjacentSquare(deltaNorthWest, deltaSouthWest, 270, 270, length);
        bottomLeft.Locations = Hexagon(length)(bottomLeft.Centre).Locations;

        adjacentPolygons.Add(topRight);
        adjacentPolygons.Add(centreRight);
        adjacentPolygons.Add(bottomRight);
        adjacentPolygons.Add(bottomLeft);
        adjacentPolygons.Add(centreLeft);
        adjacentPolygons.Add(topLeft);

        return adjacentPolygons;
    };

    return func;
}

So how I do generate all the adjacent polygons for the visible area of the map control?

Simply using one of favourite programming techniques - recursion!

This is were the order of creating the adjacent polygons becomes important. The pattern follows the principle of always creating the adjacent polygons for the current polygon and then moving to the next adjacent polygon and repeating ad infinitum. There is a break clause for the recursion to prevent an infinite loop - if the current polygon is not visible it does not calculate the adjacent polygons. This means only the visible polygons are calculated for the visible area of the map control.

Shown below is the visible order the polygons are calculated for all supported shapes, I've shown the first 12 polygons being calculated and rendered.

The first 12 squares:

The first 12 triangles:

The first 12 hexagons:

Shown below is the recursive method I use to calculate all the visible adjacent polygons. As you can see this method is short and concise - all the logic for creating the adjacent polygons is abstracted away in the AdjacentPolygonFunc. The only other logic is to work out if an adjacent polygon is visible on the screen, this uses the Intersects method on the bounding rectangle of the map control. To prevent an infinite loop occurring we also check to see if an adjacent polygon has already been created and added to the existing polygons list.

private void AddVisibleAdjacentPolygons(Polygon polygon, ICollection<Polygon> existingPolygons)
{
    var adjacentPolygons = this.selectedShape.AdjacentPolygonFunc(polygon);

    foreach (var adjacentPolygon in adjacentPolygons)
    {
        if (!this.boundingRectangle.Intersects(adjacentPolygon.TileBoundingRectangle))
        {
            continue;
        }

        if (!existingPolygons.Contains(adjacentPolygon))
        {
            existingPolygons.Add(adjacentPolygon);
            this.AddVisibleAdjacentPolygons(adjacentPolygon, existingPolygons);
        }
    }
}

The only remaining visibility issue is the use of a 'datum' for tessellation. Wikipedia defines a 'geodetic datum is a reference from which measurements are made'.

In the examples above I'm using 'Highbury Corner' in central London as the datum point for tessellating shapes, this is why the polygons do not start in the centre off the visible bounding rectangle. The idea being where every I scroll the map control to, all the rendered polygons should start from this point.

Since it's possible I'm not rendering all the polygons between the datum and the centre of the visible bounding rectangle I need to be able to work out the closes point to the centre of the visible bounding rectangle when the datum point is not visible. This is achieved using geometry and Pythagoras theorem.

I need to be able to work out the shortest hypotenuse between the datum (start) and the centre of the visible bounding rectangle (end). The horizontal and vertical components of this hypotenuses has to be an absolute factor of the height and width of the polygon. Once this value is calculated I then use the Haversine formula to calculate the exact geo-location of the centre polygon.

public static GeoCoordinate CalculateClosestPolygon(GeoCoordinate start, GeoCoordinate end, Polygon polygon)
{
    var hypotenuse = start.GetDistanceTo(end);
    if (hypotenuse == 0)
    {
        return start;
    }

    var bearing = CalculateBearing(start, end);
    var absoluteBearing = Math.Abs(bearing);

    var vertical = Math.Sin(DegreeToRadian(absoluteBearing)) * hypotenuse;
    var horizontal = Math.Cos(DegreeToRadian(absoluteBearing)) * hypotenuse;

    var rectifiedHorizontalLength = RectifiedLength(horizontal, polygon.Width);
    var rectifiedVerticalLength = RectifiedLength(vertical, polygon.Height);

    var rectifiedHypotenuseLength = Math.Sqrt((rectifiedHorizontalLength * rectifiedHorizontalLength) + (rectifiedVerticalLength * rectifiedVerticalLength));

    var closestEnd = CalculateUsingHaversine(start, rectifiedHypotenuseLength, bearing);
    return closestEnd;
}

The final issue is really an observation about memory consumption - the more polygons I add to the map control the closer I get to the 90 Mb limit. As you can see from the screenshot above I managed to get the peak memory usage over 180 Mb! You can see from the output window in visual studio the reason why the right hand screenshot above peaks at such a high value, the number of visible polygons is over 3700!

This is another demonstration of the problems you can have with a data set with a large if not infinite size. To get round this issue I would limit the minimum size of the polygons depending on the zoom level of the map control.

If you've reached this far and are still wondering what's the use of tessellating polygons over the map control?

The simple answer is to cluster pins or other visual information that has a geo-location value. In the final 'How many pins can Bing Maps handle in a WP7 app...'  I will show how you can conflate (cluster) large numbers of map pins for a relatively small visible area.

The code makes use of the WP7Contrib for the base class for the ViewModel & Model classes and is referenced as an NuGet packages. The code for this demo app is available on SkyDrive.