From fdcb6c8829a460bcc612de68676f822c045bd8b7 Mon Sep 17 00:00:00 2001
From: Your Name <you@example.com>
Date: Tue, 26 Dec 2017 15:54:59 +0100
Subject: Update assets

---
 assets/js/d3-geo-voronoi.js | 1753 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 1753 insertions(+)
 create mode 100644 assets/js/d3-geo-voronoi.js

(limited to 'assets/js/d3-geo-voronoi.js')

diff --git a/assets/js/d3-geo-voronoi.js b/assets/js/d3-geo-voronoi.js
new file mode 100644
index 0000000..3d0711c
--- /dev/null
+++ b/assets/js/d3-geo-voronoi.js
@@ -0,0 +1,1753 @@
+// https://github.com/Fil/d3-geo-voronoi Version 0.0.5. Copyright 2017 Philippe Riviere.
+(function (global, factory) {
+	typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports, require('d3-array'), require('d3-collection'), require('d3-geo'), require('d3-voronoi')) :
+	typeof define === 'function' && define.amd ? define(['exports', 'd3-array', 'd3-collection', 'd3-geo', 'd3-voronoi'], factory) :
+	(factory((global.d3 = global.d3 || {}),global.d3,global.d3,global.d3,global.d3));
+}(this, (function (exports,d3Array,d3Collection,d3Geo,d3Voronoi) { 'use strict';
+
+//
+//  Copyright (c) 2016, by Loren Petrich
+//
+// Distributed under the terms of the MIT License
+//
+//  http://lpetrich.org/
+//
+
+// Calculates the spherical Delaunay triangulation of a set of points
+// These points are entered as an array of arrays of coordinates: 0, 1, 2
+// Any extra members are ignored
+
+// FindDelaunayTriangulation(Positions) and
+// FindDelaunayTriangulationIndexed(Positions, Indices)
+// work from an array of points as specified above,
+// the second one working from a set of indices into the array,
+// and return an object with these members:
+
+// positions -- vectors on a unit sphere
+// indices -- of all the vertices
+// triangles -- array of TriangleObject
+// edges -- array of EdgeObject
+// hull -- array of vertex indices -- the convex hull
+// vor_positions -- positions of triangles' circumcircle centers (Voronoi vertices)
+// vor_edges -- pair of indices in vor_positions (empty one: [-1,-1])
+// vor_polygons -- object indexed by vertex index,
+//    and containing edges (EdgeObject), triangles (TriangleObject),
+//    and boundary (vertices in vor_positions)
+//    Open ones have a -1 at each end.
+
+// Warning: ImproveTriangulation() is mysteriously buggy
+// and is effectively disabled for now
+
+function dotprd(x,y)
+{
+	var sum = 0.0;
+	for (var ic=0; ic<3; ic++)
+		sum += x[ic]*y[ic];
+	return sum;
+}
+
+function crsprd(x,y)
+{
+	var prod = new Array(3);
+	for (var ic=0; ic<3; ic++)
+	{
+		var ic1 = ic + 1;
+		if (ic1 >= 3) ic1 -= 3;
+		var ic2 = ic + 2;
+		if (ic2 >= 3) ic2 -= 3;
+		prod[ic] = x[ic1]*y[ic2] - x[ic2]*y[ic1];
+	}
+	return prod;
+}
+
+function triple_prd(x,y,z)
+{
+	return dotprd(crsprd(x,y),z);
+}
+
+// This distance formula has some nice properties:
+// distance and not square of distance;
+// the square roots give better numerical resolution
+// distance of antipode(p) to p' = - (distance of p to p')
+// Range: -2 to +2
+function ptdist(x,y)
+{
+	var dst1 = 0.0;
+	var dst2 = 0.0;
+	for (var ic=0; ic<3; ic++)
+	{
+		var diff1 = y[ic] - x[ic];
+		dst1 += diff1*diff1;
+		var diff2 = y[ic] + x[ic];
+		dst2 += diff2*diff2;
+	}
+	return Math.sqrt(dst1) - Math.sqrt(dst2);
+}
+
+function Normalize(vec)
+{
+	var vecres = new Array(3);
+	var sum = 0.0, nrmult;
+	for (var ic=0; ic<3; ic++)
+	{
+		var val = vec[ic];
+		sum += val*val;
+	}
+	if (sum > 0)
+		nrmult = 1/Math.sqrt(sum);
+	else
+		nrmult = 0;
+	for (var ic=0; ic<3; ic++)
+	{
+		vecres[ic] = nrmult*vec[ic];
+	}
+	return vecres;
+}
+
+function crsprd_ix(Positions,ix,iy)
+{
+	return crsprd(Positions[ix],Positions[iy]);
+}
+
+function triple_prd_ix(Positions,ix,iy,iz)
+{
+	return triple_prd(Positions[ix],Positions[iy],Positions[iz]);
+}
+
+function ptdist_ix(Positions,ix,iy)
+{
+	return ptdist(Positions[ix],Positions[iy]);
+}
+
+// Returns a zero 3-vector
+function zerovec()
+{
+	var vec = new Array(3);
+	for (var ic=0; ic<3; ic++)
+		vec[ic] = 0.0;
+	return vec;
+}
+
+// Implements copying
+function vec_copy(x)
+{
+	var vec = new Array(3);
+	for (var ic=0; ic<3; ic++)
+		vec[ic] = x[ic];
+	return vec;
+}
+
+// Implements x += y
+function vec_add_to(x,y)
+{
+	for (var ic=0; ic<3; ic++)
+		x[ic] += y[ic];
+}
+
+// Implements x *= y
+function vec_mult_scalar_to(x,y)
+{
+	for (var ic=0; ic<3; ic++)
+		x[ic] *= y;
+}
+
+// Implements x - y
+function vec_difference(x,y)
+{
+	var diff = zerovec();
+	for (var ic=0; ic<3; ic++)
+		diff[ic] = x[ic] - y[ic];
+	return diff;
+}
+
+// JavaScript's counterpart of "null" / "None":
+function IsNull(x)
+{
+	return (typeof(x) == 'undefined');
+}
+
+function TrianglesEqual(tr1, tr2)
+{
+	if (IsNull(tr1)) return false;
+	if (IsNull(tr2)) return false;
+	
+	for (var iv=0; iv<3; iv++)
+		if (tr1.verts[iv] != tr2.verts[iv])
+			return false;
+	
+	return true;
+}
+
+function EdgesEqual(ed1, ed2)
+{
+	if (IsNull(ed1)) return false;
+	if (IsNull(ed2)) return false;
+	
+	for (var iv=0; iv<2; iv++)
+		if (ed1.verts[iv] != ed2.verts[iv])
+			return false;
+	
+	return true;
+}
+
+function max(x,y)
+{
+	return (y > x) ? y : x;
+}
+
+function TriangleObject(Positions, verts)
+{
+	this.verts = verts;
+	this.edges = new Array(3);
+	
+	// Find directions for testing whether a point is inside
+	this.dirs = new Array(3);
+	for (var ic=0; ic<3; ic++)
+	{
+		var ic1 = ic + 1;
+		if (ic1 >= 3) ic1 -= 3;
+		var ic2 = ic + 2;
+		if (ic2 >= 3) ic2 -= 3;
+		this.dirs[ic] = crsprd_ix(Positions,verts[ic1],verts[ic2]);
+	}
+	
+	// Tetrahedral volume factor
+	this.vol = triple_prd_ix(Positions,verts[0],verts[1],verts[2]);
+	
+	// Adjust to get the signs correct for the point-inside test;
+	// the vertex opposite the edges' vertices ought to give a dot product of 1
+	for (var ic=0; ic<3; ic++)
+		vec_mult_scalar_to(this.dirs[ic],1/this.vol);
+	
+	// Circumcircle test
+	var ccdir = zerovec();
+	for (var ic=0; ic<3; ic++)
+		vec_add_to(ccdir,this.dirs[ic]);
+	this.ccdir = Normalize(ccdir);
+	
+	var ccdsq = 0;
+	for (var ic=0; ic<3; ic++)
+		ccdsq += ptdist(this.ccdir,Positions[verts[ic]]);
+	ccdsq /= 3;
+	this.ccdsq = ccdsq;
+}
+
+// For copying in vertex info from another triangle
+TriangleObject.prototype.copy_vert_info = function(src)
+{
+	this.verts = src.verts;
+	this.dirs = src.dirs;
+	this.vol = src.vol;
+	this.ccdir = src.ccdir;
+	this.ccdsq = src.ccdsq;
+};
+
+TriangleObject.prototype.IsVertexOrderCorrect = function()
+{
+	return this.vol >= 0;
+};
+
+TriangleObject.prototype.IsPointInside = function(p)
+{
+	 for (var ic=0; ic<3; ic++)
+	 	if (dotprd(p,this.dirs[ic]) < 0) return false;
+	 
+	 return true;
+};
+
+TriangleObject.prototype.IsPointInCircumcircle = function(p)
+{
+	return (ptdist(this.ccdir,p) < this.ccdsq);
+};
+
+TriangleObject.prototype.IsVertex = function(ix)
+{
+	for (var ic=0; ic<3; ic++)
+		if (ix == this.verts[ic]) return true;
+	return false;
+};
+
+TriangleObject.prototype.VertexIndexIn = function(ix)
+{
+	for (var ic=0; ic<3; ic++)
+		if (ix == this.verts[ic]) return ic;
+	return -1;
+};
+
+TriangleObject.prototype.EdgeIndexIn = function(ed)
+{
+	for (var ic=0; ic<3; ic++)
+		if (EdgesEqual(this.edges[ic], ed)) return ic;
+	return -1;
+};
+
+function EdgeObject(verts)
+{
+	this.verts = verts;
+	this.polys = new Array(2);
+}
+
+EdgeObject.prototype.IsVertex = function(ix)
+{
+	for (var ic=0; ic<2; ic++)
+		if (ix == this.verts[ic]) return true;
+	return false;
+};
+
+EdgeObject.prototype.VertexIndexIn = function(ix)
+{
+	for (var ic=0; ic<2; ic++)
+		if (ix == this.verts[ic]) return ic;
+	return -1;
+};
+
+EdgeObject.prototype.PolyIndexIn = function(pl)
+{
+	for (var ic=0; ic<2; ic++)
+		if (TrianglesEqual(this.polys[ic],pl)) return ic;
+	return -1;
+};
+
+function EdgeCheckObject(Positions,verts)
+{
+	this.verts = verts;
+	this.pdst = ptdist_ix(Positions,verts[0],verts[1]);
+	this.direc = Normalize(crsprd_ix(Positions,verts[0],verts[1]));
+	var midpnt = zerovec();
+	vec_add_to(midpnt,Positions[verts[0]]);
+	vec_add_to(midpnt,Positions[verts[1]]);
+	this.midpnt = Normalize(midpnt);
+}
+
+// Check on the possible intersection with another edge-check object
+// return a boolean of whether or not it does
+EdgeCheckObject.prototype.intersects = function(Positions,other)
+{
+	// Assume that sharing a vertex means non-intersecting
+	for (var ic=0; ic<2; ic++)
+		for (var ict=0; ict<2; ict++)
+			if (this.verts[ic] == other.verts[ict]) return false;
+	
+	// Find intersection point; will test it and its antipode
+	var itsc = Normalize(crsprd(this.direc, other.direc));
+	
+	// Find dot product with midpoints to test if the intersection
+	// is in the near hemispheres of the lines' midpoints.
+	// If it is in both near hemispheres or both far hemispheres, it's OK
+	// In both far hemispheres: antipode is in both near hemispheres
+	var near0 = dotprd(itsc,this.midpnt) > 0;
+	var near1 = dotprd(itsc,other.midpnt) > 0;
+	if (near1 != near0) return false;
+	
+	var pd0 = [];
+	for (var ic=0; ic<2; ic++)
+	{
+		var pd = ptdist(itsc, Positions[this.verts[ic]]);
+		pd0.push(pd);
+	}
+	var pd1 = [];
+	for (var ic=0; ic<2; ic++)
+	{
+		var pd = ptdist(itsc, Positions[other.verts[ic]]);
+		pd1.push(pd);
+	}
+	var mxpd0 = max(pd0[0],pd0[1]);
+	var mxpd1 = max(pd1[0],pd1[1]);
+	if ((mxpd0 <= this.pdst) && (mxpd1 <= other.pdst) && near0) return true;
+	
+	// Handle its antipode; use antipode-related shortcuts
+	// like reversing the distance value and the hemisphere-presence value
+	vec_mult_scalar_to(itsc, -1);
+	near0 = !near0;
+	for (var ic=0; ic<2; ic++)
+	{
+		pd0[ic] = - pd0[ic];
+		pd1[ic] = - pd1[ic];
+	}
+	mxpd0 = max(pd0[0],pd0[1]);
+	mxpd1 = max(pd1[0],pd1[1]);
+	if ((mxpd0 <= this.pdst) && (mxpd1 <= other.pdst) && near0) return true;
+	
+	return false;
+};
+
+// Adds to an array if it was not already present;
+// Must resort to this kludge because JavaScript doesn't have sets
+function AddUnique(arr, x)
+{
+	for (var i=0; i<arr.length; i++)
+		if (x == arr[i]) return;
+	arr.push(x);
+}
+
+// Version for edges, since testing equality of objects
+// doesn't work that well in JavaScript
+function AddUniqueEdge(arr, ed)
+{
+	for (var i=0; i<arr.length; i++)
+		if (EdgesEqual(arr[i],ed)) return;
+	arr.push(ed);
+}
+
+// Find the set intersection
+function FindShared(arr1, arr2)
+{
+	var resarr = [];
+	for (var i1=0; i1<arr1.length; i1++)
+	{
+		var x1 = arr1[i1];
+		for (var i2=0; i2<arr2.length; i2++)
+		{
+			var x2 = arr2[i2];
+			if (x1 == x2)
+			{
+				resarr.push(x1);
+				break;
+			}
+		}
+	}
+	
+	return resarr;
+}
+
+// Version for edges
+function FindSharedEdges(arr1, arr2)
+{
+	var resarr = [];
+	for (var i1=0; i1<arr1.length; i1++)
+	{
+		var x1 = arr1[i1];
+		for (var i2=0; i2<arr2.length; i2++)
+		{
+			var x2 = arr2[i2];
+			if (EdgesEqual(x1,x2))
+			{
+				resarr.push(x1);
+				break;
+			}
+		}
+	}
+	
+	return resarr;
+}
+
+// Takes all the members of of arr2 out of arr1
+// and ignores the arr2 members not present in arr1
+function FindSetDifference(arr1, arr2)
+{
+	if (arr2.length == 0) return;
+	var diffarr = [];
+	for (var i1=0; i1<arr1.length; i1++)
+	{
+		var x1 = arr1[i1];
+		var AddThisOne = true;
+		for (var i2=0; i2<arr2.length; i2++)
+		{
+			var x2 = arr2[i2];
+			if (x2 == x1)
+			{
+				AddThisOne = false;
+				break;
+			}
+		}
+		if (AddThisOne) diffarr.push(x1);
+	}
+	
+	// Clear the array
+	arr1.splice(0,arr1.length);
+	
+	for (var i=0; i<diffarr.length; i++)
+		arr1.push(diffarr[i]);
+}
+
+// Version for edges
+function FindSetDifferenceEdges(arr1, arr2)
+{
+	if (arr2.length == 0) return;
+	var diffarr = [];
+	for (var i1=0; i1<arr1.length; i1++)
+	{
+		var x1 = arr1[i1];
+		var AddThisOne = true;
+		for (var i2=0; i2<arr2.length; i2++)
+		{
+			var x2 = arr2[i2];
+			if (EdgesEqual(x1,x2))
+			{
+				AddThisOne = false;
+				break;
+			}
+		}
+		if (AddThisOne) diffarr.push(x1);
+	}
+	
+	// Clear the array
+	arr1.splice(0,arr1.length);
+	
+	for (var i=0; i<diffarr.length; i++)
+		arr1.push(diffarr[i]);
+}
+
+
+// Specified by index ix; returns whether it was possible to do so
+function AddPointInside(TriSet, ix)
+{
+	var Positions = TriSet.positions;
+	var p = Positions[ix];
+	
+	var NumTris = TriSet.triangles.length;
+	for (var j=0; j<NumTris; j++)
+	{
+		var tri = TriSet.triangles[j];
+		if (tri.IsPointInside(p))
+		{
+			// Create three new triangles and their edges
+			var eds = tri.edges;
+			var trixs = [];
+			for (var ic=0; ic<3; ic++)
+				trixs.push(eds[ic].PolyIndexIn(tri));
+			
+			var newtris = Array(3);
+			var neweds = Array(3);
+			for (var ic=0; ic<3; ic++)
+			{
+				var ic1 = ic + 1;
+				if (ic1 >= 3) ic1 -= 3;
+				newtris[ic] = new TriangleObject(Positions,[tri.verts[ic],tri.verts[ic1],ix]);
+				neweds[ic] = new EdgeObject([tri.verts[ic],ix]);
+			}
+			
+			// Connect those triangles and edges
+			for (var ic=0; ic<3; ic++)
+			{
+				var ic1 = ic + 1;
+				if (ic1 >= 3) ic1 -= 3;
+				newtris[ic].edges[0] = neweds[ic1];
+				newtris[ic].edges[1] = neweds[ic];
+				neweds[ic].polys[0] = newtris[ic];
+				neweds[ic1].polys[1] = newtris[ic];
+			}
+		
+			// Find which external edges go with which triangles
+			for (var ic=0; ic<3; ic++)
+			{
+				var ed = eds[ic];
+				var trix = trixs[ic];
+				for (var ict=0; ict<3; ict++)
+				{
+					var newtri = newtris[ict];
+					var numverts = 0;
+					for (var iv=0; iv<2; iv++)
+					{
+						if (newtri.IsVertex(ed.verts[iv]))
+							numverts++;
+						if (numverts == 2)
+						{
+							ed.polys[trix] = newtri;
+							newtri.edges[2] = ed;
+							break;
+						}
+					}
+				}
+			}
+			
+			// Insert those triangles and edges into the lists
+			TriSet.triangles[j] = newtris[0];
+			for (var ic=1; ic<3; ic++)
+				TriSet.triangles.push(newtris[ic]);
+			for (var ic=0; ic<3; ic++)
+				TriSet.edges.push(neweds[ic]);
+			
+			// All done; indicate that the point was added
+			return true;
+		}
+	}
+	
+	// The point was inside no triangle, and thus was not added
+	return false;
+}
+
+function ImproveTriangulation(TriSet)
+{
+	var Positions = TriSet.positions;
+	var quad_verts = new Array(4);
+	for (var itr=0; itr<100; itr++)
+	{
+		var numflips = 0;
+		for (var i=0; i<TriSet.edges.length; i++)
+		{
+			var edge = TriSet.edges[i];
+			var tris = edge.polys;
+			
+			// Skip over external edges
+			if (IsNull(tris[0])) continue;
+			if (IsNull(tris[1])) continue;
+			
+			// Find the containing quadrangle's vertices
+			for (var ic=0; ic<3; ic++)
+			{
+				var ix = tris[0].verts[ic];
+				if (!edge.IsVertex(ix)) break;
+			}
+			var ic1 = ic + 1;
+			if (ic1 >= 3) ic1 -= 3;
+			var ic2 = ic + 2;
+			if (ic2 >= 3) ic2 -= 3;
+			quad_verts[0] = ix;
+			quad_verts[1] = tris[0].verts[ic1];
+			quad_verts[3] = tris[0].verts[ic2];
+			
+			for (var ic=0; ic<3; ic++)
+			{
+				var ix = tris[1].verts[ic];
+				if (!edge.IsVertex(ix)) break;
+			}
+			quad_verts[2] = ix;
+			
+			// Are the non-edge points in the other triangles' circumcircles?
+			var incc0 = tris[0].IsPointInCircumcircle(Positions[quad_verts[2]]);
+			var incc1 = tris[1].IsPointInCircumcircle(Positions[quad_verts[0]]);
+			if ((!incc0) && (!incc1)) continue;
+			
+			// Are the would-be triangles properly oriented?
+			var newtri0 = new TriangleObject(Positions, [quad_verts[0],quad_verts[1],quad_verts[2]]);
+			if (!newtri0.IsVertexOrderCorrect()) continue;
+			var newtri1 = new TriangleObject(Positions, [quad_verts[0],quad_verts[2],quad_verts[3]]);
+			if (!newtri1.IsVertexOrderCorrect()) continue;
+			
+			// If so, then flip
+			numflips++;
+			
+			// Adjust the edge and triangle memberships:
+			// 0-3 goes from 0 to 1, 1-2 goes from 1 to 0
+			for (var ic=0; ic<3; ic++)
+			{
+				var ed = tris[0].edges[ic];
+				if (EdgesEqual(ed,edge)) continue;
+				else if (ed.IsVertex(quad_verts[3]))
+				{
+					var ed03 = ed;
+					var edix03 = ic;
+					break;
+				}
+			}
+			for (var ic=0; ic<3; ic++)
+			{
+				var ed = tris[1].edges[ic];
+				if (EdgesEqual(ed,edge)) continue;
+				else if (ed.IsVertex(quad_verts[1]))
+				{
+					var ed12 = ed;
+					var edix12 = ic;
+					break;
+				}
+			}
+						
+			var trix0 = ed03.PolyIndexIn(tris[0]);
+			var trix1 = ed12.PolyIndexIn(tris[1]);
+			
+			ed03.polys[trix0] = tris[1];
+			ed12.polys[trix1] = tris[0];
+			tris[0].edges[edix03] = ed12;
+			tris[1].edges[edix12] = ed03;
+			
+			// Add the vertices
+			tris[0].copy_vert_info(newtri0);
+			tris[1].copy_vert_info(newtri1);
+			edge.verts = [quad_verts[0], quad_verts[2]];
+		}
+		if (numflips == 0) break;
+	}
+}
+
+
+function FindConvexHull(TriSet)
+{
+	// var Positions = TriSet.positions;
+	
+	// Find boundary loop -- use as convex hull
+	var NextVertex = new Object;
+	var VtxStart = -1;
+	for (var i=0; i<TriSet.edges.length; i++)
+	{
+		var edge = TriSet.edges[i];
+		
+		// Find a boundary one -- look for the triangle that it contains
+		if (IsNull(edge.polys[0]))
+		{
+			if (IsNull(edge.polys[1]))
+				continue;
+			else
+				var tri = edge.polys[1];
+		}
+		else
+		{
+			if (IsNull(edge.polys[1]))
+				var tri = edge.polys[0];
+			else
+				continue;
+		}
+		
+		// Ensure that the hull is in the same direction as the triangles
+		var ix0 = edge.verts[0];
+		var ix1 = edge.verts[1];
+		var posdiff = tri.VertexIndexIn(ix1) - tri.VertexIndexIn(ix0);
+		if (posdiff < 0) posdiff += 3;
+		if (posdiff != 1)
+		{
+			var ixs = ix0;
+			ix0 = ix1;
+			ix1 = ixs;
+		}
+		
+		NextVertex[ix0] = ix1;
+		VtxStart = ix0;
+	}
+	
+	if (VtxStart >= 0)
+	{
+		var ix = VtxStart;
+		var hull = [ix];
+		while(true)
+		{
+			var ixnext = NextVertex[ix];
+			if (ixnext == VtxStart) break;
+			hull.push(ixnext);
+			ix = ixnext;
+		}
+		
+		TriSet.hull = hull;
+	}
+}
+
+
+// Finds the dual of the Delaunay triangulation
+// Won't bother to do the sort of connectivity
+// that was necessary for the Delaunay triangulation
+function FindVoronoiDiagram(TriSet)
+{
+	// Special cases: 3 or fewer points
+	if (TriSet.triangles.length == 1)
+	{
+		// A single triangle
+		if (TriSet.hull.length == 3)
+		{
+			var tri = TriSet.triangles[0];
+			TriSet.vor_positions.push(tri.ccdir);
+			for (var k=0; k<3; k++)
+			{
+				var kx = k + 1;
+				if (kx >= 3) kx = 0;
+				var ky = k - 1;
+				if (ky < 0) ky = 2;
+				
+				var v1 = TriSet.positions[TriSet.hull[k]];
+				var v2 = TriSet.positions[TriSet.hull[kx]];
+				var posdiff = vec_difference(v2,v1);
+				TriSet.vor_positions.push(Normalize(crsprd(posdiff,tri.ccdir)));
+				TriSet.vor_edges.push([0, k+1, 4]);
+				
+				var ix = TriSet.hull[k];
+				TriSet.vor_polygons[ix] = new Object;
+				var vor_poly = 	TriSet.vor_polygons[ix];
+				var iy = TriSet.hull[ky];
+				for (var l=0; l<3; l++)
+				{
+					var edge = TriSet.edges[l];
+					var shrd = FindShared([iy,ix],edge.verts);
+					if (shrd.length == 2) break;
+				}
+				
+				vor_poly.edges = [edge];
+				vor_poly.triangles = [tri];
+				vor_poly.boundary = [0, ky+1, 4, k+1];
+			}
+			var ept = vec_copy(tri.ccdir);
+			vec_mult_scalar_to(ept,-1);
+			TriSet.vor_positions.push(ept);
+		}
+		return;
+	}
+	else if (TriSet.triangles.length == 0)
+	{
+		// A biangle
+		if (TriSet.hull.length == 2)
+		{
+			var v0 = TriSet.positions[TriSet.hull[0]];
+			var v1 = TriSet.positions[TriSet.hull[1]];
+			
+			var vt0 = zerovec();
+			vec_add_to(vt0,v0);
+			vec_add_to(vt0,v1);
+			vt0 = Normalize(vt0);
+			TriSet.vor_positions.push(vt0);
+			
+			var vt1 = Normalize(crsprd(v0,v1));
+			TriSet.vor_positions.push(vt1);
+			
+			var vt2 = vec_copy(vt0);
+			vec_mult_scalar_to(vt2,-1);
+			TriSet.vor_positions.push(vt2);
+			
+			var vt3 = vec_copy(vt1);
+			vec_mult_scalar_to(vt3,-1);
+			TriSet.vor_positions.push(vt3);
+			
+			TriSet.vor_edges.push([0, 1, 2, 3, 0]);
+			
+			edge = TriSet.edges[0];
+			for (var k=0; k<2; k++)
+			{
+				var ix = TriSet.hull[k];
+				TriSet.vor_polygons[ix] = new Object;
+				var vor_poly = 	TriSet.vor_polygons[ix];
+				vor_poly.edges = [edge];
+				vor_poly.triangles = [0];
+				if (k == 0)
+					vor_poly.boundary = [0, 1, 2, 3];
+				else if (k == 1)
+					vor_poly.boundary = [0, 3, 2, 1];
+			}
+		}
+		return;
+	}
+	
+	// Create the array of Voronoi-vertex positions:
+	// Add indices to the triangle objects for convenience
+	for (var i=0; i<TriSet.triangles.length; i++)
+	{
+		var tri = TriSet.triangles[i];
+		tri.index = i;
+		TriSet.vor_positions.push(tri.ccdir);
+	}
+	
+	// Voronoi edges: a cinch
+	// Voronoi edges parallel original edges
+	for (var i=0; i<TriSet.edges.length; i++)
+	{
+		var edge = TriSet.edges[i];
+		if (!IsNull(edge.polys[0]) && !IsNull(edge.polys[1]))
+		{
+			var vor_edge = [edge.polys[0].index, edge.polys[1].index];
+			TriSet.vor_edges.push(vor_edge);
+		}
+	}
+	
+	// Voronoi polygons: -1 at ends means an open one
+	
+	// First, collect the edges and triangles at each vertex
+	// Put them into vor_polygons, because each one
+	// is for each original vertex
+	for (var i=0; i<TriSet.indices.length; i++)
+	{
+		var ix = TriSet.indices[i];
+		TriSet.vor_polygons[ix] = new Object;
+		
+		var vor_poly = 	TriSet.vor_polygons[ix];
+		vor_poly.edges = [];
+		vor_poly.triangles = [];
+		vor_poly.boundary = [];
+	}
+	
+	for (var i=0; i<TriSet.edges.length; i++)
+	{
+		var edge = TriSet.edges[i];
+		for (var ic=0; ic<2; ic++)
+			TriSet.vor_polygons[edge.verts[ic]].edges.push(edge);
+	}
+	
+	for (var i=0; i<TriSet.triangles.length; i++)
+	{
+		var tri = TriSet.triangles[i];
+		for (var ic=0; ic<3; ic++)
+			TriSet.vor_polygons[tri.verts[ic]].triangles.push(tri);
+	}
+	
+	for (var i=0; i<TriSet.indices.length; i++)
+	{
+		var ix = TriSet.indices[i];
+		var vor_poly = 	TriSet.vor_polygons[ix];
+		
+		// First triangle
+		var init_tri = vor_poly.triangles[0];
+		var tri = init_tri;
+		vor_poly.boundary.push(tri.index);
+		
+		// First edge
+		for (var ic=0; ic<3; ic++)
+		{
+			var edge = tri.edges[ic];
+			if (edge.IsVertex(ix))
+				break;
+		}
+		var init_edge = edge;
+		
+		// The next triangle and edge
+		var IsInside = false;
+		while(true)
+		{
+			var iv = edge.PolyIndexIn(tri);
+			tri = edge.polys[1-iv];
+			if (IsNull(tri)) break;
+			if (TrianglesEqual(tri,init_tri))
+			{
+				IsInside = true;
+				break;
+			}
+			
+			vor_poly.boundary.push(tri.index);
+			
+			for (var ic=0; ic<3; ic++)
+			{
+				var next_edge = tri.edges[ic];
+				if (EdgesEqual(next_edge,edge)) continue;
+				if (next_edge.IsVertex(ix))
+				{
+					edge = next_edge;
+					break;
+				}
+			}
+		}
+		
+		if (!IsInside)
+		{
+			vor_poly.boundary.reverse();
+			tri = init_tri;
+			
+			// First edge the other way
+			for (var ic=0; ic<3; ic++)
+			{
+				edge = tri.edges[ic];
+				if (EdgesEqual(edge,init_edge)) continue;
+				if (edge.IsVertex(ix))
+					break;
+			}
+			
+			while(true)
+			{
+				var iv = edge.PolyIndexIn(tri);
+				tri = edge.polys[1-iv];
+				if (IsNull(tri)) break;
+				
+				vor_poly.boundary.push(tri.index);
+				
+				for (var ic=0; ic<3; ic++)
+				{
+					var next_edge = tri.edges[ic];
+					if (EdgesEqual(next_edge,edge)) continue;
+					if (next_edge.IsVertex(ix))
+					{
+						edge = next_edge;
+						break;
+					}
+				}
+			}
+		}
+		
+		// Add -1 on ends for open polygon:
+		if (!IsInside)
+		{
+			vor_poly.boundary.reverse();
+			vor_poly.boundary.push(-1);
+			vor_poly.boundary.reverse();
+			vor_poly.boundary.push(-1);
+		}
+	}
+	
+	// Handle the area outside of the convex hull
+	if (TriSet.hull.length >= 3)
+	{
+		// Set up the initial boundary lines
+		// The boundary lines contain:
+		// Index of Voronoi vertex / triangle center / intersection (in VorPos)
+		// Indices of original vertices on each side of the line
+		var VorBdLns = new Array();
+		var Positions = TriSet.positions;
+		var hlen = TriSet.hull.length;
+		for (var ic=0; ic<hlen; ic++)
+		{
+			var ix = TriSet.hull[ic];
+			var icx = ic + 1;
+			if (icx >= hlen) icx = 0;
+			var ixa = TriSet.hull[icx];
+			var edset1 = TriSet.vor_polygons[ix].edges;
+			var edset2 = TriSet.vor_polygons[ixa].edges;
+			var edsetshr = FindSharedEdges(edset1,edset2);
+			var edge = edsetshr[0];
+			var tvrt = edge.polys[0].index;
+			var vt0 = Positions[ix];
+			var vt1 = Positions[ixa];
+			var vtdf = vec_difference(vt1,vt0);
+			// Contains: triangle index (Voronoi vertex),
+			// vertex index 1 (Voronoi region), position
+			// vertex index 2 (Voronoi region), position,
+			// great-circle normal
+			var VorBdLn = [tvrt, TriSet.vor_positions[tvrt], ix, vt0, ixa, vt1, vtdf];
+			VorBdLns.push(VorBdLn);
+		}
+		
+		// Find intersections
+		
+		while (VorBdLns.length > 3)
+		{
+			// Check all combinations of neighbors
+			var n = VorBdLns.length;
+			var itscpts = new Array();
+			var ptitscs = new Array();
+			for (var k=0; k<n; k++)
+				ptitscs.push(new Array());
+			for (var k=0; k<n; k++)
+			{
+				// Find the intersection point; use the convex hull's direction
+				var kx = k + 1;
+				if (kx >= n) kx = 0;
+				var itscpt = Normalize(crsprd(VorBdLns[k][6],VorBdLns[kx][6]));
+				vec_mult_scalar_to(itscpt,-1);
+				ptitscs[k].push(itscpts.length);
+				ptitscs[kx].push(itscpts.length);
+				itscpts.push(itscpt);
+			}
+			// Find the intersection points that are the closest to their parent points
+			for (var k=0; k<n; k++)
+			{
+				var ptitsc = ptitscs[k];
+				if (ptitsc.length >= 2)
+				{
+					var dists = new Array();
+					for (var kp=0; kp<ptitsc.length; kp++)
+						dists.push(ptdist(itscpts[ptitsc[kp]],VorBdLns[k][1]));
+					var dx = 0;
+					var dmin = dists[dx];
+					for (var dxi=0; dxi<dists.length; dxi++)
+					{
+						var dst = dists[dxi];
+						if (dst < dmin)
+						{
+							dx = dxi; dmin = dst;
+						}
+					}
+					var ptitscrd = ptitsc[dx];
+				}
+				else if (ptitsc.length == 1)
+					var ptitscrd = ptitsc[0];
+				else
+					var ptitscrd = -1;
+				ptitscs[k] = ptitscrd;
+			}
+			
+			var NewVorBdLns = new Array();
+			for (var k=0; k<n; k++)
+			{
+				// Find all matched intersection points and add them
+				var kx = k + 1;
+				if (kx >= n) kx = 0;
+				var ky = k - 1;
+				if (ky < 0) ky = n - 1;
+				// 0 is lone, 1 is leading, 2 is trailing
+				// vorvtidx is the index of the Voronoi vertex
+				var pstat = 0;
+				var ptitsc = ptitscs[k], ptitsc_next;
+				if (ptitsc != -1)
+				{
+					var ptitsc_prev = ptitscs[ky];
+					if (ptitsc == ptitsc_prev)
+						pstat = 2;
+					else
+					{
+						ptitsc_next = ptitscs[kx];
+						if (ptitsc == ptitsc_next)
+							pstat = 1;
+					}
+				}
+				
+				if (pstat == 0)
+				{
+					// Keep the Voronoi line without merging
+					NewVorBdLns.push(VorBdLns[k]);
+				}
+				else if (pstat == 1)
+				{
+					// Merge the Voronoi lines and create a new one
+					var VorBdLn0 = VorBdLns[k];
+					var VorBdLn1 = VorBdLns[kx];
+					var itscpt = itscpts[k];
+					var tvrt0 = VorBdLn0[0];
+					var tvrt1 = VorBdLn1[0];
+					var PointOK = (tvrt1 != tvrt0);
+					if (PointOK)
+					{
+						var nitx = TriSet.vor_positions.length;
+						var ix0 = VorBdLn0[2];
+						var vt0 = VorBdLn0[3];
+						var ix1 = VorBdLn1[4];
+						var vt1 = VorBdLn1[5];
+						var dst_in = undefined;
+						var dst_out = undefined;
+						for (var m=0; m<n; m++)
+						{
+							var ptstm = ptdist(VorBdLns[m][3],itscpt);
+							var mrl = m - k;
+							while (mrl < 0) mrl += n;
+							while (mrl >= n) mrl -= n;
+							if (mrl <= 2)
+							{
+								if (dst_in == undefined)
+									dst_in = ptstm;
+								else if (ptstm < dst_in)
+									dst_in = ptstm;
+							}
+							else
+							{
+								if (dst_out == undefined)
+									dst_out = ptstm;
+								else if (ptstm < dst_out)
+									dst_out = ptstm;
+							}
+						}
+						PointOK = (dst_in < dst_out);
+					}
+					if (PointOK)
+					{
+						var vtdf = vec_difference(vt1,vt0);
+						var VorBdLn = [nitx, itscpt, ix0, vt0, ix1, vt1, vtdf];
+						NewVorBdLns.push(VorBdLn);
+						TriSet.vor_positions.push(itscpt);
+						var ixi = VorBdLn0[4];
+						// Should be equal:
+						// ixi = VorBdLn2[2];
+						TriSet.vor_edges.push([tvrt0, nitx]);
+						TriSet.vor_edges.push([tvrt1, nitx]);
+						// Add the point to the center Voronoi region and close it
+						TriSet.vor_polygons[ixi].boundary.shift();
+						var vpln = TriSet.vor_polygons[ixi].boundary.length;
+						TriSet.vor_polygons[ixi].boundary[vpln-1] = nitx;
+						// Add the point to the left Voronoi region
+						if (TriSet.vor_polygons[ix0].boundary[1] == tvrt0)
+						{
+							TriSet.vor_polygons[ix0].boundary.unshift(-1);
+							TriSet.vor_polygons[ix0].boundary[1] = nitx;
+						}
+						else
+						{
+							vpln = TriSet.vor_polygons[ix0].boundary.length;
+							if (TriSet.vor_polygons[ix0].boundary[vpln-2] == tvrt0)
+							{
+								TriSet.vor_polygons[ix0].boundary.push(-1);
+								vpln = TriSet.vor_polygons[ix0].boundary.length;
+								TriSet.vor_polygons[ix0].boundary[vpln-2] = nitx;
+							}
+						}
+						// Add the point to the right Voronoi region
+						if (TriSet.vor_polygons[ix1].boundary[1] == tvrt1)
+						{
+							TriSet.vor_polygons[ix1].boundary.unshift(-1);
+							TriSet.vor_polygons[ix1].boundary[1] = nitx;
+						}
+						else
+						{
+							vpln = TriSet.vor_polygons[ix1].boundary.length;
+							if (TriSet.vor_polygons[ix1].boundary[vpln-2] == tvrt1)
+							{
+								TriSet.vor_polygons[ix1].boundary.push(-1);
+								vpln = TriSet.vor_polygons[ix1].boundary.length;
+								TriSet.vor_polygons[ix1].boundary[vpln-2] = nitx;
+							}
+						}
+					}
+					else
+					{
+						NewVorBdLns.push(VorBdLn0);
+						NewVorBdLns.push(VorBdLn1);
+					}
+				}
+				/*
+				else if (pstat == 2)
+				{
+					// Do nothing
+				}
+				*/
+			}
+			
+			if (NewVorBdLns.length == VorBdLns.length) break;
+			VorBdLns = NewVorBdLns;
+		}
+		
+		// Special cases: only two or three points left
+		if (VorBdLns.length == 2)
+		{
+			if (VorBdLns[0][0] != VorBdLns[1][0])
+			{
+				var VorLn = [];
+				for (var k=0; k<2; k++)
+				{
+					// Connecting line
+					var kx = VorBdLns[k][0];
+					VorLn.push(kx);
+					// Close the Voronoi region by deleting the end -1's
+					kx = VorBdLns[k][2];
+					TriSet.vor_polygons[kx].boundary.shift();
+					TriSet.vor_polygons[kx].boundary.pop();
+				}
+				TriSet.vor_edges.push(VorLn);
+			}
+		}
+		else if (VorBdLns.length == 3)
+		{
+			var ic0 = VorBdLns[0][0];
+			var ic1 = VorBdLns[1][0];
+			var ic2 = VorBdLns[2][0];
+			if (ic0 != ic1 && ic0 != ic2 && ic1 != ic2)
+			{
+				var nitx = TriSet.vor_positions.length;
+				var v0 = VorBdLns[0][3];
+				var v1 = VorBdLns[1][3];
+				var v2 = VorBdLns[2][3];
+				var itscpt = zerovec();
+				vec_add_to(itscpt,crsprd(v0,v1));
+				vec_add_to(itscpt,crsprd(v1,v2));
+				vec_add_to(itscpt,crsprd(v2,v0));
+				itscpt = Normalize(itscpt);
+				vec_mult_scalar_to(itscpt,-1);
+				TriSet.vor_positions.push(itscpt);
+				for (var k=0; k<3; k++)
+				{
+					// Connecting line
+					var VorBdLn = VorBdLns[k];
+					TriSet.vor_edges.push([VorBdLn[0], nitx]);
+					// Add the point to the Voronoi region and close it
+					var ix = VorBdLn[2];
+					TriSet.vor_polygons[ix].boundary.shift();
+					var vpln = TriSet.vor_polygons[ix].boundary.length;
+					TriSet.vor_polygons[ix].boundary[vpln-1] = nitx;
+				}
+			}
+		}
+	}
+	
+	// Adjust the orientations
+	for (var k=0; k<TriSet.vor_polygons.length; k++)
+	{
+		vor_poly = TriSet.vor_polygons[k];
+		if (vor_poly.boundary.length >= 3 && vor_poly.boundary[0] >= 0)
+		{
+			tri = new TriangleObject(TriSet.vor_positions,vor_poly.boundary.slice(0,3));
+			if (!tri.IsVertexOrderCorrect())
+				vor_poly.boundary.reverse();
+		}
+	}
+}
+
+
+function FindDelaunayTriangulationIndexed(Positions, Indices)
+{
+	// Create the triangle-set object
+	var TriSet = new Object;
+	TriSet.positions = Positions;
+	TriSet.indices = Indices;
+	TriSet.triangles = [];
+	TriSet.edges = [];
+	TriSet.hull = [];
+	TriSet.vor_positions = [];
+	TriSet.vor_edges = [];
+	TriSet.vor_polygons = new Object;
+		
+	// Create the first triangle, if it is possible to create any
+	if (Indices.length < 3)
+	{
+		if (Indices.length == 2)
+		{
+			TriSet.edges.push(new EdgeObject(Indices));
+			TriSet.hull = Indices;
+		}
+		FindVoronoiDiagram(TriSet);
+		return TriSet;
+	}
+	
+	var tri = new TriangleObject(Positions, Indices.slice(0,3));
+	if (!tri.IsVertexOrderCorrect())
+		tri = new TriangleObject(Positions, [Indices[0],Indices[2],Indices[1]]);
+	TriSet.triangles.push(tri);
+	
+	var echs = new Array(3);
+	
+	for (var ic=0; ic<3; ic++)
+	{
+		var ic1 = ic + 1;
+		if (ic1 >= 3) ic1 -= 3;
+		var ix = Indices[ic];
+		var ix1 = Indices[ic1];
+		var vts = [ix, ix1];
+		var edge = new EdgeObject(vts);
+		var echeck = new EdgeCheckObject(Positions,vts);
+		echeck.edge = edge;
+		echs[ic] = echeck;
+		tri.edges[ic] = edge;
+		edge.polys[0] = tri;
+		TriSet.edges.push(edge);
+	}
+	
+	// Place those crossing checkers in a boundary object;
+	// will have to use various kludges since JavaScript doesn't have sets
+	var BoundaryVerts = Indices.slice(0,3);
+	var BoundaryEdges = echs;
+	var Verts = Object;
+	for (var ic=0; ic<3; ic++)
+	{
+		var ic1 = ic + 2;
+		if (ic1 >= 3) ic1 -= 3;
+		var ix = Indices[ic];
+		Verts[ix] = [echs[ic],echs[ic+1]];
+	}
+	
+	// Add points until it is no longer possible
+	for (var i=3; i<Indices.length; i++)
+	{
+		var ix = Indices[i];
+		
+		// If possible, add the point inside
+		if (AddPointInside(TriSet, ix)) continue;
+		
+		// Point was not inside
+		Verts[ix] = [];
+		var NewEdges = [];
+		var VertsAddedTo = [];
+		var EdgesToDelete = [];
+		
+		// Find all the non-intersecting edges
+		for (var j=0; j<BoundaryVerts.length; j++)
+		{
+			var ix1 = BoundaryVerts[j];
+			var echk = new EdgeCheckObject(Positions, [ix, ix1]);
+			var DoesIntersect = false;
+			for (var k=0; k<BoundaryEdges.length; k++)
+			{
+				var echk1 = BoundaryEdges[k];
+				DoesIntersect = echk.intersects(Positions,echk1);
+				if (DoesIntersect) break;
+			}
+			if (DoesIntersect) continue;
+			
+			var edge = new EdgeObject(echk.verts);
+			echk.edge = edge;
+			AddUniqueEdge(NewEdges,echk);
+			AddUniqueEdge(Verts[ix],echk);
+			AddUnique(VertsAddedTo,ix);
+			AddUniqueEdge(Verts[ix1],echk);
+			AddUnique(VertsAddedTo,ix1);
+		}
+		
+		// Add the new vertex itself
+		AddUnique(BoundaryVerts,ix);
+		
+		// Find all the triangles
+		for (var j=0; j<BoundaryEdges.length; j++)
+		{
+			var echk = BoundaryEdges[j];
+			var echks = [];
+			
+			for (var iv=0; iv<2; iv++)
+			{
+				var vset = FindSharedEdges(Verts[ix],Verts[echk.verts[iv]]);
+				if (vset.length == 0) continue;
+				echks.push(vset[0]);
+			}
+			if (echks.length < 2) continue;
+			
+			var empt_indx = -1;
+			for (var iv=0; iv<2; iv++)
+			{
+				if (IsNull(echk.edge.polys[iv]))
+				{
+					empt_indx = iv;
+					break;
+				}
+			}
+			if (empt_indx < 0) continue;
+			
+			var oldtri = echk.edge.polys[1-empt_indx];
+			var v0 = echk.verts[0];
+			var i0 = oldtri.VertexIndexIn(v0);
+			var v1 = echk.verts[1];
+			var i1 = oldtri.VertexIndexIn(v1);
+			var i01 = i1 - i0;
+			if (i01 < 0) i01 += 3;
+			if (i01 == 1)
+			{
+				// Order in original: other, v0, v1
+				var NewTriVerts = [ix, v1, v0];
+			}
+			else if (i01 == 2)
+			{
+				// Order in original: other, v1, v0
+				var NewTriVerts = [ix, v0, v1];
+			}
+			var tri = new TriangleObject(Positions, NewTriVerts);
+			if (!tri.IsVertexOrderCorrect()) continue;
+			
+			// Add the new triangle
+			// Also, add the new edges,
+			// or remove them from the lists if necessary
+			TriSet.triangles.push(tri);
+			echk.edge.polys[empt_indx] = tri;
+			tri.edges[0] = echk.edge;
+			tri.edges[1] = echks[0].edge;
+			tri.edges[2] = echks[1].edge;
+			AddUniqueEdge(EdgesToDelete, echk);
+			
+			for (var iv=0; iv<2; iv++)
+			{
+				var echki = echks[iv];
+				if (IsNull(echki.edge.polys[0]))
+				{
+					echki.edge.polys[0] = tri;
+					TriSet.edges.push(echki.edge);
+				}
+				else
+				{
+					echki.edge.polys[1] = tri;
+					AddUniqueEdge(EdgesToDelete,echki);
+				}
+			}
+		}
+		
+		// Add the new edges and remove the edges and vertices
+		// that are now in the interior
+		for (var j=0; j<NewEdges.length; j++)
+			AddUniqueEdge(BoundaryEdges,NewEdges[j]);
+		FindSetDifferenceEdges(BoundaryEdges, EdgesToDelete);
+		
+		var BoundaryVertsToRemove = [];
+		for (var j=0; j<VertsAddedTo.length; j++)
+		{
+			var ixa = VertsAddedTo[j];
+			FindSetDifferenceEdges(Verts[ixa],EdgesToDelete);
+			if (Verts[ixa].length == 0)
+				BoundaryVertsToRemove.push(ixa);
+		}
+		FindSetDifference(BoundaryVerts, BoundaryVertsToRemove);
+	}
+	
+	// Improve it iteratively
+	ImproveTriangulation(TriSet);
+	
+	// Find the boundary line of this region
+	FindConvexHull(TriSet);
+	
+	// Find the regions around each point:
+	FindVoronoiDiagram(TriSet);
+	
+	return TriSet;
+}
+
+function FindDelaunayTriangulation(Positions)
+{
+	var Indices = new Array(Positions.length);
+	for (var i=0; i<Indices.length; i++)
+		Indices[i] = i;
+	
+	return FindDelaunayTriangulationIndexed(Positions, Indices); 
+}
+
+//
+// (c) 2016 Philippe Riviere
+//
+// https://github.com/Fil/
+//
+// This software is distributed under the terms of the MIT License
+
+var geoVoronoi = function () {
+    var radians = Math.PI / 180;
+
+    var cartesian = function (spherical) {
+        var lambda = spherical[0] * radians,
+            phi = spherical[1] * radians,
+            cosphi = Math.cos(phi);
+        return [
+            cosphi * Math.cos(lambda),
+            cosphi * Math.sin(lambda),
+            Math.sin(phi)
+        ];
+    };
+
+    var spherical = function (cartesian) {
+        var r = Math.sqrt(cartesian[0] * cartesian[0] + cartesian[1] * cartesian[1]),
+            lat = Math.atan2(cartesian[2], r),
+            lng = Math.atan2(cartesian[1], cartesian[0]);
+        return [lng / radians, lat / radians];
+    };
+
+    var mapline = function (Positions, Verts) {
+        return Verts
+            .map(function (v) {
+                return spherical(Positions[v]);
+            });
+    };
+
+
+    var diagram = d3Voronoi.voronoi()([]);
+
+    var DT = diagram.DT = null,
+        sites = diagram.sites = [],
+        pos = diagram.pos = [],
+        x = function (d) {
+            if (typeof d == 'object' && 'type' in d) {
+                return d3Geo.geoCentroid(d)[0];
+            }
+            if (0 in d) return d[0];
+        },
+        y = function (d) {
+            if (typeof d == 'object' && 'type' in d) {
+                return d3Geo.geoCentroid(d)[1];
+            }
+            if (0 in d) return d[1];
+        };
+
+
+    var voro = function (data) {
+        diagram._hull = diagram._polygons = diagram._links = diagram._triangles = null;
+
+        if (typeof data == 'object' && data.type == 'FeatureCollection') {
+            data = data.features;
+        }
+        sites = data.map(function (site, i) {
+            site.index = i;
+            return site;
+        });
+        pos = data.map(function (site) {
+            return [x(site), y(site)];
+        });
+        DT = FindDelaunayTriangulation(pos.map(cartesian));
+        return diagram;
+    };
+
+    diagram.links = voro.links = function (s) {
+        if (s) voro(s);
+        if (diagram._links) return diagram._links;
+
+        var _index = d3Collection.map();
+
+        var features = DT.edges.map(function (i, n) {
+
+            _index.set(d3Array.extent(i.verts), n);
+
+            var properties = {
+                source: sites[i.verts[0]],
+                target: sites[i.verts[1]],
+                urquhart: true, // will be changed to false later
+                length: d3Geo.geoDistance(pos[i.verts[0]], pos[i.verts[1]])
+            };
+
+            // add left and right sites (?)
+
+            // make GeoJSON
+            return {
+                type: "Feature",
+                geometry: {
+                    type: 'LineString',
+                    coordinates: [spherical(DT.positions[i.verts[0]]), spherical(DT.positions[i.verts[1]])]
+                },
+                properties: properties
+            };
+        });
+
+        // Urquhart Graph? tag longer link from each triangle
+        DT.triangles.forEach(function (t) {
+            var l = 0,
+                length = 0,
+                remove, v;
+            for (var j = 0; j < 3; j++) {
+                v = d3Array.extent([t.verts[j], t.verts[(j + 1) % 3]]);
+                var n = _index.get(v);
+                length = features[n].properties.length;
+                if (length > l) {
+                    l = length;
+                    remove = n;
+                }
+            }
+            features[remove].properties.urquhart = false;
+        });
+
+        return diagram._links = {
+            type: "FeatureCollection",
+            features: features
+        };
+    };
+
+    diagram.triangles = voro.triangles = function (s) {
+        if (s) voro(s);
+        if (diagram._triangles) return diagram._triangles;
+
+        var features = DT.triangles
+            .map(function (t) {
+                t.spherical = t.verts.map(function (v) {
+                    return DT.positions[v];
+                })
+                    .map(spherical);
+
+                // correct winding order
+                if (t.ccdsq < 0) {
+                    t.spherical = t.spherical.reverse();
+                    t.ccdsq *= -1;
+                }
+
+                return t;
+            })
+
+            // make geojson
+            .map(function (t) {
+                return {
+                    type: "Feature",
+                    geometry: {
+                        type: "Polygon",
+                        coordinates: [t.spherical.concat([t.spherical[0]])]
+                    },
+                    properties: {
+                        sites: t.verts.map(function (i) {
+                            return sites[i];
+                        }),
+                        area: t.vol, // steradians
+                        circumcenter: spherical(t.ccdir),
+                        // ccdsq is *not* the geodesic distance
+                        /* circumradius: (2-t.ccdsq) * 53 */
+                    }
+                }
+            });
+
+        return diagram._triangles = {
+            type: "FeatureCollection",
+            features: features
+        };
+
+    };
+
+    diagram.polygons = voro.polygons = function (s) {
+        if (s) voro(s);
+        if (diagram._polygons) return diagram._polygons;
+
+        var features = DT.indices.map(function (i, n) {
+            var vor_poly = DT.vor_polygons[DT.indices[i]];
+            var geometry = {};
+
+            if (vor_poly == undefined) {
+                geometry.type = "Sphere";
+            } else {
+                var line = mapline(DT.vor_positions,
+                    vor_poly.boundary.concat([vor_poly.boundary[0]])
+                );
+
+                // correct winding order
+                var b = {
+                    type: "Polygon",
+                    coordinates: [[pos[i], line[0], line[1], pos[i]]]
+                };
+                if (d3Geo.geoArea(b) > 2 * Math.PI + 1e-10) {
+                    line = line.reverse();
+                }
+
+                geometry.type = "Polygon";
+                geometry.coordinates = [line];
+            }
+
+            return {
+                type: "Feature",
+                geometry: geometry,
+                properties: {
+                    site: sites[i],
+                    sitecoordinates: pos[i],
+                    neighbours: vor_poly.edges.map(function (e) {
+                        return e.verts.filter(function (j) {
+                            return j !== i;
+                        })[0];
+                    })
+                }
+            };
+        });
+
+        return diagram._polygons = {
+            type: "FeatureCollection",
+            features: features
+        };
+    };
+
+    diagram.hull = voro.hull = function (s) {
+        if (s) voro(s);
+        if (diagram._hull) return diagram._hull;
+
+        if (!DT.hull.length) {
+            return null; // What is a null GeoJSON?
+        }
+
+        // seems that DT.hull is always clockwise
+        var hull = DT.hull.reverse();
+
+        // make GeoJSON
+        return diagram._hull = {
+            type: "Feature",
+            geometry: {
+                type: "Polygon",
+                coordinates: [hull.concat([hull[0]]).map(function (i) {
+                    return pos[i];
+                })]
+            },
+            properties: {
+                sites: hull.map(function (i) {
+                    return sites[i];
+                })
+            }
+        };
+    };
+
+
+    diagram.find = function (x, y, radius) {
+        var features = diagram.polygons().features;
+
+        // optimization: start from most recent result
+        var i, next = diagram.find.found || 0;
+        var cell = features[next] || features[next = 0];
+
+        var dist = d3Geo.geoDistance([x, y], cell.properties.sitecoordinates);
+        do {
+            cell = features[i = next];
+            next = null;
+            cell.properties.neighbours.forEach(function (e) {
+                var ndist = d3Geo.geoDistance([x, y], features[e].properties.sitecoordinates);
+                if (ndist < dist) {
+                    dist = ndist;
+                    next = e;
+                    return;
+                }
+            });
+
+        } while (next !== null);
+        diagram.find.found = i;
+        if (!radius || dist < radius * radius) return cell.properties.site;
+    };
+
+    voro.x = function (f) {
+        if (!f) return x;
+        x = f;
+        return voro;
+    };
+    voro.y = function (f) {
+        if (!f) return y;
+        y = f;
+        return voro;
+    };
+    voro.extent = function (f) {
+        if (!f) return null;
+        return voro;
+    };
+    voro.size = function (f) {
+        if (!f) return null;
+        return voro;
+    };
+
+    return voro;
+
+};
+
+exports.geoVoronoi = geoVoronoi;
+
+Object.defineProperty(exports, '__esModule', { value: true });
+
+})));
-- 
cgit v1.2.3