diff options
Diffstat (limited to 'assets/js/d3-geo-voronoi.js')
-rw-r--r-- | assets/js/d3-geo-voronoi.js | 1753 |
1 files changed, 1753 insertions, 0 deletions
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 }); + +}))); |