Gridding Geometries (or, “Creating patchwork animals in SQL Server”)

This is one of those things that I can’t imagine anybody would ever really want to do but, seeing as I haven’t posted anything for a while, I thought I’d write about it just in case it’s of use to someone…

First, suppose you had an arbitrary geometry. Here’s an example:

DECLARE @snake geometry = geometry::Parse('CIRCULARSTRING(0 0, 5 -5, 10 0, 15 5, 20 0, 25 -5, 29 -1)').STBuffer(2).STDifference(geometry::Parse('MULTIPOINT(28.5 -0.5, 29.75 -0.8)').STBuffer(0.6)).STUnion(geometry::Parse('MULTIPOINT(28.5 -0.5, 29.75 -0.8)').STBuffer(0.2)).STUnion(geometry::Parse('POLYGON((29.2 0.5, 29.5 2.7, 29.79 2.2, 30.22 2.75, 30 0.5, 29.2 0.5))'));
SELECT @snake;

In in the immortal words of Rolf Harris, “can you guess what it is yet?”.

Yes, that’s right – it’s a snake and it looks like this:

Now suppose that you wanted to divide that geometry up into a series of individual polygonal tiles on a regular-spaced grid. You could do so using the following T-SQL:

-- Determine the extent of the grid required to cover the geometry
DECLARE @grid geometry = @snake.STEnvelope();
DECLARE @gridwidth float = @grid.STPointN(3).STX - @grid.STPointN(1).STX;
DECLARE @gridheight float = @grid.STPointN(3).STY - @grid.STPointN(2).STY;

-- Work out the lower-left corner of the grid
DECLARE @xoffset float = @grid.STPointN(1).STX,
        @yoffset float = @grid.STPointN(1).STY;

-- Specify the number of columns and rows in the grid
DECLARE @numcolumns float = 80, @numrows float = 35;

-- Work out the height and width of each tile
DECLARE @cellwidth float = @gridwidth / @numcolumns;
DECLARE @cellheight float = @gridheight / @numrows;

-- Create a table variable to hold the tiles
DECLARE @gridcells table (id int, geom geometry)

-- Now, loop through and cut up the geometry into tiles
DECLARE @x int = 0, @y int = 0;
WHILE @y < @numrows
BEGIN
WHILE @x < @numcolumns
BEGIN
INSERT INTO @gridcells VALUES(
@y * @numcolumns + @x,
'POLYGON((' + cast(@xoffset + (@x * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + (@y * @cellheight) AS varchar(32)) + ','
+ cast(@xoffset + ((@x + 1)  * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + (@y * @cellheight) AS varchar(32)) + ','
+ cast(@xoffset + ((@x + 1)  * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + ((@y + 1) * @cellheight) AS varchar(32)) + ','
+ cast(@xoffset + (@x * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + ((@y + 1) * @cellheight) AS varchar(32)) + ','
+ cast(@xoffset + (@x * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + (@y * @cellheight) AS varchar(32)) + '))'
)
SET @x = @x + 1;
END
SET @x = 0;
SET @y = @y + 1;
END

-- Finally, select the tiles from the table variable
SELECT geom.STIntersection(@snake) FROM @gridcells ORDER BY id;

My snake has now been chopped up into 800 little tiles, and here’s what the result looks like:

 

I wonder if this is how Elmer the Elephant was created?

Google’s April Fools’ Day 8-bit Map

Today, Google launched its annual April Fool’s joke – go to http://maps.google.com right now and you’ll see a new “Quest” map style in addition to the regular satellite and aerial map styles.

Looks a little…. familiar? Here’s the 8-bit tileset for Bing Maps I described just over two weeks ago:

You’re welcome, Google, the invoice is in the post…

Creating “The Legend of Zelda” map tiles from Bing Maps aerial imagery

Time for something light-hearted, I feel. Following on from last year’s Bing Maps tribute to Grand Theft Auto, I thought I’d do another videogame-themed map hack inspired by perhaps my most-loved game of all time – the Legend of Zelda.

Almost the entire year of my 10-year old life was spent playing the Legend of Zelda on the Nintendo Entertainment System. In those days, games rarely had their own in-game maps, and it was long before GameFAQs or Prima strategy guides existed, so the only way to keep track of where you’d been in an adventure game was to create your own maps. My mum would bring me home sheets of graph paper from work, and I would spend hours laboriously drawing the layout of each room in every dungeon, and how each of the screens in the overworld connected to each other. (If you’ve ever tried to do this yourself, you’ll appreciate the ingenuity of the “lost woods”, which are very hard to map on paper!).

As a result, I’ve always personally found the overworld maps of Zelda games to be a thing of beauty. Here’s the world map for the original Legend Of Zelda:

This one’s from Zelda II: The Adventure of Link:

and here’s The Legend of Zelda: A Link to the Past:

The world of Hyrule is undeniably beautiful – its deserts, rivers, and forests becoming ever more intricate in each successive game in the Zelda series. But wouldn’t it be cool to create an 8bit-style tileset of the real world that you navigate around?

Now, I’m no graphic artist and, even if I was, the thought of drawing the 275 billion tiles required to cover the world at zoom level 19 would probably put me off the idea anyway. So, I need some way to automatically determine whether a given area should be represented as a water tile, or a grass tile, say. And I wondered if it was possible to do this based on a simple image analysis of the corresponding aerial tile from Bing Maps.

Converting Bing Maps Aerial Imagery to 8bit tiles

So, here’s the plan: I build a Bing Maps application that uses a custom tile layer. The handler for the tile layer requests aerial imagery from Bing Maps. But, instead of returning the tile image direct to the client, it converts the features on that tile into appropriate 8bit sprites and returns a new tile composed of those sprites to the client instead.

To choose which tile images to use, I’ll divide each 256px x 256 tile retrieved from Bing Maps into smaller 32px x 32px subtiles. I’ll then do a simple RGB analysis to determine the “average” colour of each subtile.

If the average colour of that subtile is mostly green (that is, when comparing the RGB component values, G > B > R) then I’ll insert a grass tile:
If the average colour is mostly blue (that is, B > R and B > G) then insert a water tile:
If the average colour is mostly dark brown (R >= G and R >= B and B <= 90) then insert a rock tile:
If the average colour is mostly light brown (R >= G and R >= B and B > 90) then insert a dirt/sand tile:
To add a bit of variety, I’ll also make it so that every grass tile has a 1/50 chance of being replaced with a flower tile:

You can obviously tweak these, or add more rules to fit the tileset as you see fit.

Here’s the code to implement these rules in a C# handler:

<%@ WebHandler Language="C#" Class="Handler" %>

using System;
using System.Web;
using System.Drawing;
using System.Data.SQLite;
using System.IO;
using System.Drawing.Imaging;

public class Handler : IHttpHandler
{

  public void ProcessRequest(HttpContext context)
  {
    string quadKey = context.Request.QueryString["q"];

    // Create a random number generator to vary grass tiles
    Random random = new Random();
    
    // Retrieve the original aerial tile image
    System.Net.WebRequest request = System.Net.WebRequest.Create("http://ecn.t" + quadKey[quadKey.Length - 1] + ".tiles.virtualearth.net/tiles/a" + quadKey + "?g=774&mkt=en-gb&lbl=l1&stl=h&n=z");
    System.Net.WebResponse response = request.GetResponse();
    System.IO.Stream responseStream = response.GetResponseStream();
    Bitmap mapImage = new Bitmap(responseStream);

    // Set up a graphics class based on the tile bitmap
    using (Graphics graphics = System.Drawing.Graphics.FromImage(mapImage))
    {
      int tileWidth = 32, tileHeight=32;
      int numRows = 8, numCols = 8;
      
      // Loop through the aerial tile in 16x16 squares
      for (int row = 0; row < numRows; row++)
      {
        for (int col = 0; col < numCols; col++)
        {
          // Get the bitmap data for this square
          Rectangle tileRect = new Rectangle(col * tileWidth, row * tileHeight, tileWidth, tileHeight);
          System.Drawing.Imaging.BitmapData bmData = mapImage.LockBits(
            tileRect,
            System.Drawing.Imaging.ImageLockMode.ReadOnly,
            mapImage.PixelFormat);

          // Get the average R,G,B values of pixels in this square
          long[] totals = new long[] { 0, 0, 0 };
          unsafe
          {
            // 24bit image so 3 bytes per pixel (PNG + transparency would be 4)
            int PixelSize = 3;
            for (int y = 0; y < tileHeight; y++)
            {
              byte* p = (byte*)bmData.Scan0 + (y * bmData.Stride);
              for (int x = 0; x < tileWidth; x++)
              {
                totals[0] += p[x * PixelSize]; // Blue
                totals[1] += p[x * PixelSize + 1]; // Green
                totals[2] += p[x * PixelSize + 2]; // Red
              }
            }
          }
          mapImage.UnlockBits(bmData);

          // Work out the average RGB colour value
          int avgB = (int)(totals[0] / (bmData.Width * bmData.Height));
          int avgG = (int)(totals[1] / (bmData.Width * bmData.Height));
          int avgR = (int)(totals[2] / (bmData.Width * bmData.Height));

          // Determine average colour of this tile
          Color fillColour = Color.FromArgb(avgR, avgG, avgB);

          // Snow
          if (fillColour.G >= 225 && fillColour.B >= 225 && fillColour.R >= 225)
          {
            graphics.FillRectangle(Brushes.White, tileRect);
          }
          
          // Grass
          else if (fillColour.G >= fillColour.B && fillColour.G >= fillColour.R)
          {
            // Random chance of being a flower tile
            if (random.Next(0, 50) < 1)
            {
              Bitmap grassImage = (Bitmap)Image.FromFile(@"C:\Users\Alastair\Documents\Visual Studio 2008\WebSites\BingMaps\v7\ZeldaMap\flowers.png");
              graphics.DrawImage(grassImage, tileRect);
            }
            else
            {
              Bitmap grassImage = (Bitmap)Image.FromFile(@"C:\Users\Alastair\Documents\Visual Studio 2008\WebSites\BingMaps\v7\ZeldaMap\grass.png");
              graphics.DrawImage(grassImage, tileRect);
            }
          }
            
          // Water
          else if (fillColour.B >= fillColour.G && fillColour.B >= fillColour.R)
          {
            graphics.FillRectangle(new SolidBrush(Color.FromArgb(255,48,89,166)), tileRect);
          }
            
          // Dirt/Sand
          else if (fillColour.R >= fillColour.G && fillColour.R >= fillColour.B && fillColour.B <= 90)
          {
            Bitmap dirtImage = (Bitmap)Image.FromFile(@"C:\Users\Alastair\Documents\Visual Studio 2008\WebSites\BingMaps\v7\ZeldaMap\rock.png");
            graphics.DrawImage(dirtImage, tileRect);
          }
            
          // Rock
          else if (fillColour.R >= fillColour.G && fillColour.R >= fillColour.B && fillColour.B > 90)
          {
            Bitmap rockImage = (Bitmap)Image.FromFile(@"C:\Users\Alastair\Documents\Visual Studio 2008\WebSites\BingMaps\v7\ZeldaMap\dirt.png");
            graphics.DrawImage(rockImage, tileRect);
          }
            
          // Default - just fill tile in single average colour
          else {
            graphics.FillRectangle(new SolidBrush(fillColour), tileRect);
          }
        }
      }
    }

    // Send the resulting image back to the client
    context.Response.ContentType = "image/png";
    System.IO.MemoryStream memStream = new System.IO.MemoryStream();
    mapImage.Save(memStream, System.Drawing.Imaging.ImageFormat.Png);
    memStream.WriteTo(context.Response.OutputStream);
  }


  public bool IsReusable
  {
    get
    {
      return false;
    }
  }

}

(Note that, because I’m accessing aerial tiles directly from the Bing Maps tile servers rather than through the API, this method is unsupported and probably not allowed in a production environment. But then, why on earth would you ever want to do this in a production environment? It’s just for fun!)

Based on the rules above, let’s see how certain tiles are rendered:

0331	0331
1
01	01

That’s pretty much the effect I was after, and here’s what it looks like then you create a custom tile layer based on the above handler displayed on a Bing Maps control at zoom level 1. Pretty neat, huh?:

Adding a Player Sprite

To finish this example off, I want to add a player to the world as well. He’ll stay in the middle of the map, and the tiles will scroll under his feet as you press the cursor keys for him to navigate around. I’ll create a few simple frames of animation to show him walking in different directions, as shown in the sprite strip to the left.

The player sprite will appear in a div placed in the centre of the map, and I’ll use javascript to change the background of the div to the appropriate frame depending on whether the player is moving up, right, down, or left. Since there are three frames of animation for each direction, I’ll create a simple frame matrix that will loop through the relevant frames.

var frameMatrix = [1, 2, 0, 4, 5, 3, 7, 8, 6, 10, 11, 9];

In addition to animating the background of the player sprite div, pressing a cursor key will also cause the map to scroll in the appropriate direction. For this, I’ve added a custom Pan() method that extends the Microsoft.Maps.Map prototype.

Here’s the code listing for the HTML page:

<!DOCTYPE HTML PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html  xmlns="http://www.w3.org/1999/xhtml">
<head>
<title></title>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<script type="text/javascript" src="http://ecn.dev.virtualearth.net/mapcontrol/mapcontrol.ashx?v=7.0"></script>
<script type="text/javascript">

  // Define some globals for handling key input
  var keys = { UP: 38, DOWN: 40, LEFT: 37, RIGHT: 39 };
  var keysPressed = [];
  // These used to track player animation frames
  var frameMatrix = [1, 2, 0, 4, 5, 3, 7, 8, 6, 10, 11, 9];
  var currentFrame = 0;
  // The map. Obviously.
  var map = null;
  // Add a convenient method to pan the map  
  Microsoft.Maps.Map.prototype.Pan = function(x, y) {
    var pos = this.tryPixelToLocation(new Microsoft.Maps.Point(this.getWidth() / 2 + x, this.getHeight() / 2 - y), Microsoft.Maps.PixelReference.control);
    this.setView({ center: pos });
  };
  // Set the frame rate of how often the map moves/player animates
  var int = window.setInterval("animate()", 150);

  // This function gets called once per animation loop
  function animate() {
    // Only look at the first key pressed
    switch (keysPressed[0]) {
      case keys.UP:
        // If we've reached the end of the "up" animation frames, start again
        if (currentFrame >= 2) { currentFrame = 0; }
        // otherwise, go on to the next "up" animation frame
        else { currentFrame++ };
        // Pan the map up
        map.Pan(0, 10);
        break;
      case keys.DOWN:
        // If we've reached the end of the "down" animation frames, start again
        if (currentFrame < 3 || currentFrame >= 5) { currentFrame = 3; }
        // otherwise, go on to the next "down" animation frame
        else { currentFrame++ };
        // Pan the map down
        map.Pan(0, -10);
        break;
      case keys.LEFT:
        // If we've reached the end of the "left" animation frames, start again
        if (currentFrame < 6 || currentFrame >= 8) { currentFrame = 6; }
        // otherwise, go on to the next "left" animation frame
        else { currentFrame++ };
        // Pan the map left
        map.Pan(-10, 0);
        break;
      case keys.RIGHT:
        // If we've reached the end of the "right" animation frames, start again
        if (currentFrame < 9 || currentFrame >= 11) { currentFrame = 9; }
        // otherwise, go on to the next "right" animation frame
        else { currentFrame++ };
        // Pan the map right
        map.Pan(10, 0);
        break;
    }
    
    // Update the offset of the background of the player div
    document.getElementById("player").style.backgroundPosition = "0 " + -currentFrame * 32 + "px";
  }


  // This function keeps track of the key(s) currently being pressed.
  // Uses an array so that, for example, can use diagonal movement when
  // two keys are pressed.
  function handleKeyDown(e) {
    switch (e.keyCode) {
      case keys.UP:
        if (!inArray(keys.UP, keysPressed) && !inArray(keys.DOWN, keysPressed)) {
          keysPressed.push(keys.UP);  
        }
        break;
      case keys.DOWN:
        if (!inArray(keys.DOWN, keysPressed) && !inArray(keys.UP, keysPressed)) {
          keysPressed.push(keys.DOWN);
        }
        break;
      case keys.LEFT:
        if (!inArray(keys.LEFT, keysPressed) && !inArray(keys.RIGHT, keysPressed)) {
          keysPressed.push(keys.LEFT);
        }
        break;
      case keys.RIGHT: // Right:
        if (!inArray(keys.RIGHT, keysPressed) && !inArray(keys.LEFT, keysPressed)) {
          keysPressed.push(keys.RIGHT);
        }
        break;
    }
    e.handled = true;
  }

  // Remove key from the array when they are released
  function handleKeyUp(e) {
        keysPressed = removeFromArray(e.keyCode, keysPressed);
  }

  // Utility methods
  function inArray(element, arr) {
    for (var i = 0; i < arr.length; i++) {
      if (arr[i] === element) {
        return true;
      }
    }
    return false;
  }
  function removeFromArray(element, arr) {
    for (var i = 0; i < arr.length; i++) {
      if (arr[i] == element)
        arr.splice(i, 1);
    }
    return arr;
  }

  // Create the map
  function GetMap() {
    map = new Microsoft.Maps.Map(document.getElementById("mapDiv"),
      { credentials: "GETYEROWNBINGMAPSKEY!",
        center: new Microsoft.Maps.Location(50, 0),
        mapTypeId: Microsoft.Maps.MapTypeId.mercator,
        zoom: 5,
        inertiaIntensity: 0.00001,
        tileBuffer:3,
        disableKeyboardInput: true
      });

    // Create a custom tile source that points to the .NET handler
    var tileSource = new Microsoft.Maps.TileSource({ uriConstructor: location.href.substring(0, location.href.lastIndexOf('/')) + "/" + "handler.ashx?q={quadkey}" });

    // Construct the layer using the tile source
    var tilelayer = new Microsoft.Maps.TileLayer({ mercator: tileSource, opacity: 1 });

    // Push the tile layer to the map
    map.entities.push(tilelayer);

    // Position the player sprite in the middle of the map
    document.getElementById("player").style.backgroundImage = "url('player.png')";
    document.getElementById("player").style.left = map.getWidth() / 2 + "px";
    document.getElementById("player").style.top = map.getHeight() / 2 + "px";

    // Attach event handlers
    window.addEventListener('keydown', handleKeyDown, false);
    window.addEventListener('keyup', handleKeyUp, false);
  }  
</script>
</head>
<body onload="GetMap();">
  <div id='mapDiv' style="position:relative; width:640px; height:480px;"></div>
  <div id='player' style="position: absolute; z-index:10; width:20px; height:32px; background-repeat:no-repeat;"></div>
</body>
</html>

Finally, here’s a video of the finished product. Watch me walk from somewhere near Paris to Valencia and then onto Naples, all in glorious dynamically-generated 8-bit graphics!

Drawing Fractals with SQL Server Spatial

I can’t think of any practical purpose for this code, but here’s a recursive SQLCLR procedure to draw a Sierpinski triangle fractal (or, at least, an approximation of one to a given level of detail) in SQL Server just because you can.

Adapted from http://www.codeproject.com/KB/silverlight/SierpinskiTriangle.aspx

using System;
using System.Collections.Generic;
using System.Text;
using Microsoft.SqlServer.Types;
using System.Data.SqlTypes;
using System.Data.SqlClient;
using Microsoft.SqlServer.Server;
using System.Data;

namespace ProSpatial
{
  public partial class StoredProcedures
  {
    [Microsoft.SqlServer.Server.SqlProcedure]
    public static void SierpinskiTriangle(int size)
    {

      // Set properties of exterior equilateral triangle
      int w = size;  // Width. e.g. 512
      int h = (int)(w * Math.Sqrt(3) / 2);  // Height
      int[] x = { 0, w, w/2 };  // x vertices
      int[] y = { h, h, 0 };  // y vertices

      // Create a new SqlGeometry Instance to hold the output
      SqlGeometry Triangles = new SqlGeometry();
      Triangles.STSrid = 0;

      // Start recursion
      Triangles = drawSierpinskiTriangle(x, y, w/2, 2, Triangles);

      // Send the results back to SQL Server
      SendResults(Triangles);
    }


    private static SqlGeometry drawSierpinskiTriangle(int[] x, int[] y, int d, int dMin, SqlGeometry Triangles)
    {

      // If triangles are too small to render then make this the last recursion
      if (d <= dMin)
      {
        // Create a new triangle and add it to the collection
        SqlGeometry Polygon = TriangleFromPoints(x[0], y[0], x[1], y[1],  x[2], y[2]);
        Triangles = Triangles.STUnion(Polygon);
      }
      else
      {
        // Calculate centre of each side
        int xMc = (x[0] + x[1]) / 2, yMc = (y[0] + y[1]) / 2;
        int xMb = (x[0] + x[2]) / 2, yMb = (y[0] + y[2]) / 2;
        int xMa = (x[1] + x[2]) / 2, yMa = (y[1] + y[2]) / 2;

        // Subdivide into three new triangles
        int[] xNew1 = { x[0], xMc, xMb };
        int[] yNew1 = { y[0], yMc, yMb };
        Triangles = drawSierpinskiTriangle(xNew1, yNew1, d / 2, dMin, Triangles);

        int[] xNew2 = { x[1], xMc, xMa };
        int[] yNew2 = { y[1], yMc, yMa };
        Triangles = drawSierpinskiTriangle(xNew2, yNew2, d / 2, dMin, Triangles);

        int[] xNew3 = { x[2], xMb, xMa };
        int[] yNew3 = { y[2], yMb, yMa };
        Triangles = drawSierpinskiTriangle(xNew3, yNew3, d / 2, dMin, Triangles);
      }

      // Recursion finished - return the result
      return Triangles;
    }

    // Send the results back to the client
    private static void SendResults(SqlGeometry Triangles)
    {
      // Define the metadata of the results column
      SqlMetaData metadata = new SqlMetaData("Triangle", SqlDbType.Udt, typeof(SqlGeometry));

      // Create a record based on this metadata
      SqlDataRecord record = new SqlDataRecord(metadata);
      record.SetValue(0, Triangles);

      // Send the results back to the client
      SqlContext.Pipe.Send(record);
    }

    // Construct a triangle from 3 vertices
    private static SqlGeometry TriangleFromPoints(double x0, double y0, double x1, double y1, double x2, double y2)
    {
      SqlGeometryBuilder TriangleBuilder = new SqlGeometryBuilder();
      TriangleBuilder.SetSrid(0);
      TriangleBuilder.BeginGeometry(OpenGisGeometryType.Polygon);
      TriangleBuilder.BeginFigure(x0, y0);
      TriangleBuilder.AddLine(x1, y1);
      TriangleBuilder.AddLine(x2, y2);
      TriangleBuilder.AddLine(x0, y0);
      TriangleBuilder.EndFigure();
      TriangleBuilder.EndGeometry();
      return TriangleBuilder.ConstructedGeometry;
    }

  }
}

Import the assembly, register the procedure, and then execute in SQL Server. It requires a single parameter, size, which determines the overall scale of the result.

/* Import Assembly */
CREATE ASSEMBLY SierpinskiTriangle
FROM 'c:\users\alastair\documents\SierpinskiTriangle.dll'
WITH PERMISSION_SET = SAFE;
GO

/* Register function */
CREATE PROCEDURE dbo.SierpinskiTriangle(@size int)
AS EXTERNAL NAME SierpinskiTriangle.[ProSpatial.StoredProcedures].SierpinskiTriangle;
GO

/* Execute Procedure */
EXEC dbo.SierpinskiTriangle 512;

When executed with a size of 512, the result is a single MultiPolygon instance containing a set of Sierpinski triangles as follows: