[grits] Add support for Mercator projections in tiles

[Andy Spencer <andys src gnome org>]

commit 04171adb63df1084005faadd94419a64db072596
Author: Andy Spencer <andy753421 gmail com>
Date:   Sat Jan 21 04:16:25 2012 +0000

    Add support for Mercator projections in tiles
    
    This split tiles in the correct

 src/data/grits-tms.c     |   62 +++++++----------------------
 src/grits-viewer.h       |    6 +-
 src/objects/grits-tile.c |   96 +++++++++++++++++++++++++++++++---------------
 src/objects/grits-tile.h |    3 +
 4 files changed, 86 insertions(+), 81 deletions(-)
---
diff --git a/src/data/grits-tms.c b/src/data/grits-tms.c
index cfe84f5..48f5b0c 100644
--- a/src/data/grits-tms.c
+++ b/src/data/grits-tms.c
@@ -25,63 +25,31 @@
 
 static gchar *_make_uri(GritsTms *tms, GritsTile *tile)
 {
-
-#if 0
-	/* This doesn't make any sense.. */
-	gdouble lon_rad = deg2rad(tile->edge.n + tile->edge.s)/2;
-	gdouble lat_rad = deg2rad(tile->edge.e + tile->edge.w)/2;
-	g_message("%lf,%lf", lat_rad, lon_rad);
-
-	/* Reproject the coordinates to the Mercator projection: */
-	gdouble x = lon_rad;
-	gdouble y = log(tan(lat_rad) + 1.0/cos(lat_rad));
-
-	/* Transform range of x and y to 0 - 1 and shift origin to top left */
-	x = (1.0 + (x / G_PI)) / 2.0;
-	y = (1.0 - (y / G_PI)) / 2.0;
-
-	/* Calculate the number of tiles across the map, n, using 2^zoom */
-	gint zoom = 0;
-	for (GritsTile *tmp = tile->parent; tmp; tmp = tmp->parent)
-		zoom++;
-	gint n = pow(2, zoom);
-
-	/* Multiply x and y by n. Round results down to give tilex and tiley. */
-	gint xtile = x * n;
-	gint ytile = y * n;
-
-	g_message("xy=%f,%f  zoom=%d  n=%d  xy_tiles=%d,%d",
-			x, y, zoom, n, xtile, ytile);
-#endif
-
-#if 1
-	/* This is broken */
 	gint zoom = 0;
 	for (GritsTile *tmp = tile->parent; tmp; tmp = tmp->parent)
 		zoom++;
 	gint breath = pow(2,zoom);
 
-	gdouble lon_pos =  (tile->edge.e+tile->edge.w)/2 + 180;
-	gdouble lat_pos = -(tile->edge.n+tile->edge.s)/2 +  90 - 4.9489;
+	gdouble lat_top = asinh(tan(deg2rad(tile->edge.n)));
+	gdouble lat_bot = asinh(tan(deg2rad(tile->edge.s)));
 
-	gdouble lon_total = 360;
-	gdouble lat_total = 85.0511*2;
+	gdouble lat_mid = (lat_top + lat_bot)/2.0;
+	gdouble lon_mid = (tile->edge.e + tile->edge.w)/2.0;
 
-	gdouble lon_pct = lon_pos / lon_total;
-	gdouble lat_pct = lat_pos / lat_total;
+	gdouble lat_pos = 1.0 - (lat_mid + G_PI) / (2.0*G_PI);
+	gdouble lon_pos = (lon_mid + 180.0) / 360.0;
 
-	gint xtile = lon_pct * breath;
-	gint ytile = lat_pct * breath;
+	gint xtile = lon_pos * breath;
+	gint ytile = lat_pos * breath;
 
-	//g_message("bbok=%f,%f,%f,%f",
+	//g_message("tile=%f,%f,%f,%f t=%p p=%p",
 	//		tile->edge.n, tile->edge.s,
-	//		tile->edge.e, tile->edge.w);
-	//g_message("pos=%f,%f total=%f,%f pct=%f,%f tile=%d,%d",
-	//		lon_pos,   lat_pos,
-	//		lon_total, lat_total,
-	//		lon_pct,   lat_pct,
-	//		xtile,     ytile);
-#endif
+	//		tile->edge.e, tile->edge.w, tile, tile->parent);
+	//g_message("top=%lf->%lf bot=%lf->%lf pos=%lf,%lf tile=%d,%d,%d",
+	//		tile->edge.n, lat_top,
+	//		tile->edge.s, lat_bot,
+	//		lat_pos, lon_pos,
+	//		zoom, xtile, ytile);
 
 	// http://tile.openstreetmap.org/<zoom>/<xtile>/<ytile>.png
 	return g_strdup_printf("%s/%d/%d/%d.%s",
diff --git a/src/grits-viewer.h b/src/grits-viewer.h
index b7ff216..15d9d70 100644
--- a/src/grits-viewer.h
+++ b/src/grits-viewer.h
@@ -39,10 +39,10 @@
 #include <glib-object.h>
 
 /* Projections */
-enum {
+typedef enum {
+	GRITS_PROJ_LATLON,
 	GRITS_PROJ_MERCATOR,
-	GRITS_PROJ_EQUIRECTANGULAR,
-} GritsProjection;
+} GritsProj;
 
 /* Rendering levels */
 /**
diff --git a/src/objects/grits-tile.c b/src/objects/grits-tile.c
index c8326da..68e71cc 100644
--- a/src/objects/grits-tile.c
+++ b/src/objects/grits-tile.c
@@ -33,6 +33,7 @@
  */
 
 #include <config.h>
+#include <math.h>
 #include "gtkgl.h"
 #include "grits-tile.h"
 
@@ -61,6 +62,8 @@ GritsTile *grits_tile_new(GritsTile *parent,
 	tile->atime  = time(NULL);
 	grits_bounds_set_bounds(&tile->coords, 0, 1, 1, 0);
 	grits_bounds_set_bounds(&tile->edge, n, s, e, w);
+	if (parent)
+		tile->proj = parent->proj;
 	return tile;
 }
 
@@ -133,6 +136,43 @@ static gboolean _grits_tile_precise(GritsPoint *eye, GritsBounds *bounds,
 	       tile_res < view_res;
 }
 
+static void _grits_tile_split_latlon(GritsTile *tile)
+{
+	const gdouble rows = G_N_ELEMENTS(tile->children);
+	const gdouble cols = G_N_ELEMENTS(tile->children[0]);
+	const gdouble lat_dist = tile->edge.n - tile->edge.s;
+	const gdouble lon_dist = tile->edge.e - tile->edge.w;
+	const gdouble lat_step = lat_dist / rows;
+	const gdouble lon_step = lon_dist / cols;
+
+	int row, col;
+	grits_tile_foreach_index(tile, row, col)
+		tile->children[row][col] = grits_tile_new(tile,
+			tile->edge.n - lat_step*(row+0),  // north
+			tile->edge.n - lat_step*(row+1),  // south
+			tile->edge.w + lon_step*(col+1),  // east
+			tile->edge.w + lon_step*(col+0)); // west
+}
+
+static void _grits_tile_split_mercator(GritsTile *tile)
+{
+	GritsTile *child = NULL;
+	GritsBounds tmp = tile->edge;
+
+	/* Project */
+	tile->edge.n = asinh(tan(deg2rad(tile->edge.n)));
+	tile->edge.s = asinh(tan(deg2rad(tile->edge.s)));
+
+	_grits_tile_split_latlon(tile);
+
+	/* Convert back to lat-lon */
+	tile->edge = tmp;
+	grits_tile_foreach(tile, child) {
+		child->edge.n = rad2deg(atan(sinh(child->edge.n)));
+		child->edge.s = rad2deg(atan(sinh(child->edge.s)));
+	}
+}
+
 /**
  * grits_tile_update:
  * @root:      the root tile to split
@@ -149,43 +189,37 @@ static gboolean _grits_tile_precise(GritsPoint *eye, GritsBounds *bounds,
  * the tile is recursively subdivided until a sufficient resolution is
  * achieved.
  */
-void grits_tile_update(GritsTile *root, GritsPoint *eye,
+void grits_tile_update(GritsTile *tile, GritsPoint *eye,
 		gdouble res, gint width, gint height,
 		GritsTileLoadFunc load_func, gpointer user_data)
 {
-	root->atime = time(NULL);
 	//g_debug("GritsTile: update - %p->atime = %u",
-	//		root, (guint)root->atime);
-	const gdouble rows = G_N_ELEMENTS(root->children);
-	const gdouble cols = G_N_ELEMENTS(root->children[0]);
-	const gdouble lat_dist = root->edge.n - root->edge.s;
-	const gdouble lon_dist = root->edge.e - root->edge.w;
-	const gdouble lat_step = lat_dist / rows;
-	const gdouble lon_step = lon_dist / cols;
-	int row, col;
-	grits_tile_foreach_index(root, row, col) {
-		GritsBounds edge;
-		edge.n = root->edge.n-(lat_step*(row+0));
-		edge.s = root->edge.n-(lat_step*(row+1));
-		edge.e = root->edge.w+(lon_step*(col+1));
-		edge.w = root->edge.w+(lon_step*(col+0));
-
-		GritsTile **child = &root->children[row][col];
-		if (!_grits_tile_precise(eye, &edge, res,
-				width/cols, height/rows)) {
-			if (!*child) {
-				*child = grits_tile_new(root, edge.n, edge.s,
-						edge.e, edge.w);
-				load_func(*child, user_data);
-			}
-			grits_tile_update(*child, eye,
-					res, width, height,
-					load_func, user_data);
-			GRITS_OBJECT(*child)->hidden = FALSE;
-		} else if (*child) {
-			GRITS_OBJECT(*child)->hidden = TRUE;
+	//		tile, (guint)tile->atime);
+
+	if (tile == NULL)
+		return;
+
+	GRITS_OBJECT(tile)->hidden = TRUE;
+	if (_grits_tile_precise(eye, &tile->edge, res, width, height))
+		return;
+	tile->atime = time(NULL);
+	GRITS_OBJECT(tile)->hidden = FALSE;
+
+	if (!tile->data)
+		load_func(tile, user_data);
+
+	if (!tile->children[0][0]) {
+		switch (tile->proj) {
+		case GRITS_PROJ_LATLON:   _grits_tile_split_latlon(tile);   break;
+		case GRITS_PROJ_MERCATOR: _grits_tile_split_mercator(tile); break;
 		}
 	}
+
+	GritsTile *child;
+	grits_tile_foreach(tile, child)
+		grits_tile_update(child, eye, res, width, height,
+				load_func, user_data);
+
 }
 
 /**
diff --git a/src/objects/grits-tile.h b/src/objects/grits-tile.h
index 309f4bc..f2b706f 100644
--- a/src/objects/grits-tile.h
+++ b/src/objects/grits-tile.h
@@ -53,6 +53,9 @@ struct _GritsTile {
 
 	/* Last access time (for garbage collection) */
 	time_t atime;
+
+	/* Projection used by tile data */
+	GritsProj proj;
 };
 
 struct _GritsTileClass {