From 3a2a2e452d4df6ab5975f59b10418ddaa243d8c5 Mon Sep 17 00:00:00 2001
From: Jon Atkins <github@jonatkins.com>
Date: Mon, 20 Jan 2014 08:24:01 +0000
Subject: [PATCH] more s2geometry code we have enough now to calculate the cell
 containing a lat/lng, and to calculate the lat/lng for a cell center/boundary
 not certain the hibbert curve stuff is right yet...

---
 external/s2geometry.js | 201 +++++++++++++++++++++++++++++++++++------
 1 file changed, 174 insertions(+), 27 deletions(-)

diff --git a/external/s2geometry.js b/external/s2geometry.js
index 78f106f8..d7fc661d 100644
--- a/external/s2geometry.js
+++ b/external/s2geometry.js
@@ -10,18 +10,24 @@
 // to convert a lat,lng into a cell id:
 // - convert lat,lng to x,y,z
 // - convert x,y,z into face,u,v
-// - u,v scaled to s,t (i,j?) with quadratic formula
-// - s,t (i,j?) converted to a position along a Hubbert space-filling curve
+// - u,v scaled to s,t with quadratic formula
+// - s,t converted to integer i,j offsets
+// - i,j converted to a position along a Hubbert space-filling curve
 // - combine face,position to get the cell id
 
-
+//NOTE: compared to the google S2 geometry library, we vary from their code in the following ways
+// - cell IDs: they combine face and the hilbert curve position into a single 64 bit number. this gives efficient space
+//             and speed. javascript doesn't have appropriate data types, and speed is not cricical, so we use
+//             as [face,[bitpair,bitpair,...]] instead
+// - i,j: they always use 30 bits, adjusting as needed. we use 0 to (1<<level)-1 instead
+//        (so GetSizeIJ for a cell is always 1)
 
 (function() {
 
 window.S2 = {};
 
 
-S2.LatLngToXYZ = function(latLng) {
+var LatLngToXYZ = function(latLng) {
   var d2r = L.LatLng.DEG_TO_RAD;
 
   var phi = latLng.lat*d2r;
@@ -32,6 +38,15 @@ S2.LatLngToXYZ = function(latLng) {
   return [Math.cos(theta)*cosphi, Math.sin(theta)*cosphi, Math.sin(phi)];
 };
 
+var XYZToLatLng = function(xyz) {
+  var r2d = L.LatLng.RAD_TO_DEG;
+
+  var lat = Math.atan2(xyz[2], Math.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]));
+  var lng = Math.atan2(xyz[1], xyz[0]);
+
+  return L.latLng(lat*r2d, lng*r2d);
+};
+
 var largestAbsComponent = function(xyz) {
   var temp = [Math.abs(xyz[0]), Math.abs(xyz[1]), Math.abs(xyz[2])];
 
@@ -39,7 +54,7 @@ var largestAbsComponent = function(xyz) {
     if (temp[0] > temp[2]) {
       return 0;
     } else {
-      return 1;
+      return 2;
     }
   } else {
     if (temp[1] > temp[2]) {
@@ -67,25 +82,10 @@ var faceXYZToUV = function(face,xyz) {
   return [u,v];
 }
 
-var singleSTtoUV = function(st) {
-
-  if (st >= 0.5) {
-    return (1/3.0) * (4*st*st - 1);
-  } else {
-    return (1/3.0) * (1 - (4*(1-st)*(1-st)));
-  }
-}
-
-var singleUVtoST = function(uv) {
-  if (uv >= 0) {
-    return 0.5 * Math.sqrt (1 + 3*uv);
-  } else {
-    return 1 - 0.5 * Math.sqrt (1 - 3*uv);
-  }
-}
 
 
-S2.XYZToFaceUV = function(xyz) {
+
+var XYZToFaceUV = function(xyz) {
   var face = largestAbsComponent(xyz);
 
   if (xyz[face] < 0) {
@@ -97,15 +97,162 @@ S2.XYZToFaceUV = function(xyz) {
   return [face, uv];
 };
 
-S2.UVToST = function(uv) {
-  return [singleUVToST(uv[0]), singleUVToST(uv[1])];
+var FaceUVToXYZ = function(face,uv) {
+  var u = uv[0];
+  var v = uv[1];
+
+  switch (face) {
+    case 0: return [ 1, u, v];
+    case 1: return [-u, 1, v];
+    case 2: return [-u,-v, 1];
+    case 3: return [-1,-v,-u];
+    case 4: return [ v,-1,-u];
+    case 5: return [ v, u,-1];
+    default: throw {error: 'Invalid face'};
+  }
 };
 
 
-var singleSTtoIJ = function(st) {
-  var ij = Math.floor(st * kMaxSize);
-  return Math.max(0, Math.min(kMaxSize-1, ij));
+var STToUV = function(st) {
+  var singleSTtoUV = function(st) {
+    if (st >= 0.5) {
+      return (1/3.0) * (4*st*st - 1);
+    } else {
+      return (1/3.0) * (1 - (4*(1-st)*(1-st)));
+    }
+  }
+
+  return [singleSTtoUV(st[0]), singleSTtoUV(st[1])];
 };
 
 
+
+var UVToST = function(uv) {
+  var singleUVtoST = function(uv) {
+    if (uv >= 0) {
+      return 0.5 * Math.sqrt (1 + 3*uv);
+    } else {
+      return 1 - 0.5 * Math.sqrt (1 - 3*uv);
+    }
+  }
+
+  return [singleUVtoST(uv[0]), singleUVtoST(uv[1])];
+};
+
+
+var STToIJ = function(st,order) {
+  var maxSize = (1<<order);
+
+  var singleSTtoIJ = function(st) {
+    var ij = Math.floor(st * maxSize);
+    return Math.max(0, Math.min(maxSize-1, ij));
+  };
+
+  return [singleSTtoIJ(st[0]), singleSTtoIJ(st[1])];
+};
+
+
+var IJToST = function(ij,order,offsets) {
+  var maxSize = (1<<order);
+
+  return [
+    (ij[0]+offsets[0])/maxSize,
+    (ij[1]+offsets[1])/maxSize
+  ];
+}
+
+// hilbert space-filling curve
+// based on http://blog.notdot.net/2009/11/Damn-Cool-Algorithms-Spatial-indexing-with-Quadtrees-and-Hilbert-Curves
+// note: rather then calculating the final integer hilbert position, we just return the list of quads
+// this ensures no precision issues whth large orders (S3 cell IDs use up to 30), and is more
+// convenient for pulling out the individual bits as needed later
+var pointToHilbertQuadList = function(x,y,order) {
+  var hilbertMap = {
+    'a': [ [0,'d'], [1,'a'], [3,'b'], [2,'a'] ],
+    'b': [ [2,'b'], [1,'b'], [3,'a'], [0,'c'] ],
+    'c': [ [2,'c'], [3,'d'], [1,'c'], [0,'b'] ],
+    'd': [ [0,'a'], [3,'c'], [1,'d'], [2,'d'] ]  
+  };
+
+  var currentSquare='a';
+  var positions = [];
+
+  for (var i=order-1; i>=0; i--) {
+
+    var mask = 1<<i;
+
+    var quad_x = x&mask ? 1 : 0;
+    var quad_y = y&mask ? 1 : 0;
+
+    var t = hilbertMap[currentSquare][quad_x*2+quad_y];
+
+    positions.push(t[0]);
+
+    currentSquare = t[1];
+  }
+
+  return positions;
+};
+
+
+
+
+// S2Cell class
+
+S2.S2Cell = function(){};
+
+//static method to construct
+S2.S2Cell.FromLatLng = function(latLng,level) {
+
+  var xyz = LatLngToXYZ(latLng);
+
+  var faceuv = XYZToFaceUV(xyz);
+  var st = UVToST(faceuv[1]);
+
+  var ij = STToIJ(st,level);
+
+
+  var result = new S2.S2Cell();
+  result.face = faceuv[0];
+  result.ij = ij;
+  result.level = level;
+
+  return result;
+};
+
+
+S2.S2Cell.prototype.getLatLng = function() {
+  var st = IJToST(this.ij,this.level, [0.5,0.5]);
+  var uv = STToUV(st);
+  var xyz = FaceUVToXYZ(this.face, uv);
+
+  return XYZToLatLng(xyz);  
+};
+
+S2.S2Cell.prototype.getCornerLatLngs = function() {
+  var result = [];
+  var offsets = [
+    [ 0.0, 0.0 ],
+    [ 0.0, 1.0 ],
+    [ 1.0, 1.0 ],
+    [ 1.0, 0.0 ]
+  ];
+
+  for (var i=0; i<4; i++) {
+    var st = IJToST(this.ij,this.level, offsets[i]);
+    var uv = STToUV(st);
+    var xyz = FaceUVToXYZ(this.face, uv);
+
+    result.push ( XYZToLatLng(xyz) );
+  }
+  return result;
+};
+
+
+S2.S2Cell.prototype.getFaceAndQuads = function() {
+  var quads = pointToHilbertQuadList(this.ij[0], this.ij[1], this.level);
+
+  return [this.face,quads];
+};
+
 })();