path: root/assets/js/d3-geo-voronoi.js

// 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 });

})));