(function(){d3.geom = {}; /** * Computes a contour for a given input grid function using the marching * squares algorithm. Returns the contour polygon as an array of points. * * @param grid a two-input function(x, y) that returns true for values * inside the contour and false for values outside the contour. * @param start an optional starting point [x, y] on the grid. * @returns polygon [[x1, y1], [x2, y2], …] */ d3.geom.contour = function(grid, start) { var s = start || d3_geom_contourStart(grid), // starting point c = [], // contour polygon x = s[0], // current x position y = s[1], // current y position dx = 0, // next x direction dy = 0, // next y direction pdx = NaN, // previous x direction pdy = NaN, // previous y direction i = 0; do { // determine marching squares index i = 0; if (grid(x-1, y-1)) i += 1; if (grid(x, y-1)) i += 2; if (grid(x-1, y )) i += 4; if (grid(x, y )) i += 8; // determine next direction if (i === 6) { dx = pdy === -1 ? -1 : 1; dy = 0; } else if (i === 9) { dx = 0; dy = pdx === 1 ? -1 : 1; } else { dx = d3_geom_contourDx[i]; dy = d3_geom_contourDy[i]; } // update contour polygon if (dx != pdx && dy != pdy) { c.push([x, y]); pdx = dx; pdy = dy; } x += dx; y += dy; } while (s[0] != x || s[1] != y); return c; }; // lookup tables for marching directions var d3_geom_contourDx = [1, 0, 1, 1,-1, 0,-1, 1,0, 0,0,0,-1, 0,-1,NaN], d3_geom_contourDy = [0,-1, 0, 0, 0,-1, 0, 0,1,-1,1,1, 0,-1, 0,NaN]; function d3_geom_contourStart(grid) { var x = 0, y = 0; // search for a starting point; begin at origin // and proceed along outward-expanding diagonals while (true) { if (grid(x,y)) { return [x,y]; } if (x === 0) { x = y + 1; y = 0; } else { x = x - 1; y = y + 1; } } } /** * Computes the 2D convex hull of a set of points using Graham's scanning * algorithm. The algorithm has been implemented as described in Cormen, * Leiserson, and Rivest's Introduction to Algorithms. The running time of * this algorithm is O(n log n), where n is the number of input points. * * @param vertices [[x1, y1], [x2, y2], …] * @returns polygon [[x1, y1], [x2, y2], …] */ d3.geom.hull = function(vertices) { if (vertices.length < 3) return []; var len = vertices.length, plen = len - 1, points = [], stack = [], i, j, h = 0, x1, y1, x2, y2, u, v, a, sp; // find the starting ref point: leftmost point with the minimum y coord for (i=1; i= (x2*x2 + y2*y2)) { points[i].index = -1; } else { points[u].index = -1; a = points[i].angle; u = i; v = j; } } else { a = points[i].angle; u = i; v = j; } } // initialize the stack stack.push(h); for (i=0, j=0; i<2; ++j) { if (points[j].index !== -1) { stack.push(points[j].index); i++; } } sp = stack.length; // do graham's scan for (; j 0; } // Note: requires coordinates to be counterclockwise and convex! d3.geom.polygon = function(coordinates) { coordinates.area = function() { var i = 0, n = coordinates.length, a = coordinates[n - 1][0] * coordinates[0][1], b = coordinates[n - 1][1] * coordinates[0][0]; while (++i < n) { a += coordinates[i - 1][0] * coordinates[i][1]; b += coordinates[i - 1][1] * coordinates[i][0]; } return (b - a) * .5; }; coordinates.centroid = function(k) { var i = -1, n = coordinates.length - 1, x = 0, y = 0, a, b, c; if (!arguments.length) k = -1 / (6 * coordinates.area()); while (++i < n) { a = coordinates[i]; b = coordinates[i + 1]; c = a[0] * b[1] - b[0] * a[1]; x += (a[0] + b[0]) * c; y += (a[1] + b[1]) * c; } return [x * k, y * k]; }; // The Sutherland-Hodgman clipping algorithm. coordinates.clip = function(subject) { var input, i = -1, n = coordinates.length, j, m, a = coordinates[n - 1], b, c, d; while (++i < n) { input = subject.slice(); subject.length = 0; b = coordinates[i]; c = input[(m = input.length) - 1]; j = -1; while (++j < m) { d = input[j]; if (d3_geom_polygonInside(d, a, b)) { if (!d3_geom_polygonInside(c, a, b)) { subject.push(d3_geom_polygonIntersect(c, d, a, b)); } subject.push(d); } else if (d3_geom_polygonInside(c, a, b)) { subject.push(d3_geom_polygonIntersect(c, d, a, b)); } c = d; } a = b; } return subject; }; return coordinates; }; function d3_geom_polygonInside(p, a, b) { return (b[0] - a[0]) * (p[1] - a[1]) < (b[1] - a[1]) * (p[0] - a[0]); } // Intersect two infinite lines cd and ab. function d3_geom_polygonIntersect(c, d, a, b) { var x1 = c[0], x2 = d[0], x3 = a[0], x4 = b[0], y1 = c[1], y2 = d[1], y3 = a[1], y4 = b[1], x13 = x1 - x3, x21 = x2 - x1, x43 = x4 - x3, y13 = y1 - y3, y21 = y2 - y1, y43 = y4 - y3, ua = (x43 * y13 - y43 * x13) / (y43 * x21 - x43 * y21); return [x1 + ua * x21, y1 + ua * y21]; } // Adapted from Nicolas Garcia Belmonte's JIT implementation: // http://blog.thejit.org/2010/02/12/voronoi-tessellation/ // http://blog.thejit.org/assets/voronoijs/voronoi.js // See lib/jit/LICENSE for details. // Notes: // // This implementation does not clip the returned polygons, so if you want to // clip them to a particular shape you will need to do that either in SVG or by // post-processing with d3.geom.polygon's clip method. // // If any vertices are coincident or have NaN positions, the behavior of this // method is undefined. Most likely invalid polygons will be returned. You // should filter invalid points, and consolidate coincident points, before // computing the tessellation. /** * @param vertices [[x1, y1], [x2, y2], …] * @returns polygons [[[x1, y1], [x2, y2], …], …] */ d3.geom.voronoi = function(vertices) { var polygons = vertices.map(function() { return []; }); d3_voronoi_tessellate(vertices, function(e) { var s1, s2, x1, x2, y1, y2; if (e.a === 1 && e.b >= 0) { s1 = e.ep.r; s2 = e.ep.l; } else { s1 = e.ep.l; s2 = e.ep.r; } if (e.a === 1) { y1 = s1 ? s1.y : -1e6; x1 = e.c - e.b * y1; y2 = s2 ? s2.y : 1e6; x2 = e.c - e.b * y2; } else { x1 = s1 ? s1.x : -1e6; y1 = e.c - e.a * x1; x2 = s2 ? s2.x : 1e6; y2 = e.c - e.a * x2; } var v1 = [x1, y1], v2 = [x2, y2]; polygons[e.region.l.index].push(v1, v2); polygons[e.region.r.index].push(v1, v2); }); // Reconnect the polygon segments into counterclockwise loops. return polygons.map(function(polygon, i) { var cx = vertices[i][0], cy = vertices[i][1]; polygon.forEach(function(v) { v.angle = Math.atan2(v[0] - cx, v[1] - cy); }); return polygon.sort(function(a, b) { return a.angle - b.angle; }).filter(function(d, i) { return !i || (d.angle - polygon[i - 1].angle > 1e-10); }); }); }; var d3_voronoi_opposite = {"l": "r", "r": "l"}; function d3_voronoi_tessellate(vertices, callback) { var Sites = { list: vertices .map(function(v, i) { return { index: i, x: v[0], y: v[1] }; }) .sort(function(a, b) { return a.y < b.y ? -1 : a.y > b.y ? 1 : a.x < b.x ? -1 : a.x > b.x ? 1 : 0; }), bottomSite: null }; var EdgeList = { list: [], leftEnd: null, rightEnd: null, init: function() { EdgeList.leftEnd = EdgeList.createHalfEdge(null, "l"); EdgeList.rightEnd = EdgeList.createHalfEdge(null, "l"); EdgeList.leftEnd.r = EdgeList.rightEnd; EdgeList.rightEnd.l = EdgeList.leftEnd; EdgeList.list.unshift(EdgeList.leftEnd, EdgeList.rightEnd); }, createHalfEdge: function(edge, side) { return { edge: edge, side: side, vertex: null, "l": null, "r": null }; }, insert: function(lb, he) { he.l = lb; he.r = lb.r; lb.r.l = he; lb.r = he; }, leftBound: function(p) { var he = EdgeList.leftEnd; do { he = he.r; } while (he != EdgeList.rightEnd && Geom.rightOf(he, p)); he = he.l; return he; }, del: function(he) { he.l.r = he.r; he.r.l = he.l; he.edge = null; }, right: function(he) { return he.r; }, left: function(he) { return he.l; }, leftRegion: function(he) { return he.edge == null ? Sites.bottomSite : he.edge.region[he.side]; }, rightRegion: function(he) { return he.edge == null ? Sites.bottomSite : he.edge.region[d3_voronoi_opposite[he.side]]; } }; var Geom = { bisect: function(s1, s2) { var newEdge = { region: {"l": s1, "r": s2}, ep: {"l": null, "r": null} }; var dx = s2.x - s1.x, dy = s2.y - s1.y, adx = dx > 0 ? dx : -dx, ady = dy > 0 ? dy : -dy; newEdge.c = s1.x * dx + s1.y * dy + (dx * dx + dy * dy) * .5; if (adx > ady) { newEdge.a = 1; newEdge.b = dy / dx; newEdge.c /= dx; } else { newEdge.b = 1; newEdge.a = dx / dy; newEdge.c /= dy; } return newEdge; }, intersect: function(el1, el2) { var e1 = el1.edge, e2 = el2.edge; if (!e1 || !e2 || (e1.region.r == e2.region.r)) { return null; } var d = (e1.a * e2.b) - (e1.b * e2.a); if (Math.abs(d) < 1e-10) { return null; } var xint = (e1.c * e2.b - e2.c * e1.b) / d, yint = (e2.c * e1.a - e1.c * e2.a) / d, e1r = e1.region.r, e2r = e2.region.r, el, e; if ((e1r.y < e2r.y) || (e1r.y == e2r.y && e1r.x < e2r.x)) { el = el1; e = e1; } else { el = el2; e = e2; } var rightOfSite = (xint >= e.region.r.x); if ((rightOfSite && (el.side === "l")) || (!rightOfSite && (el.side === "r"))) { return null; } return { x: xint, y: yint }; }, rightOf: function(he, p) { var e = he.edge, topsite = e.region.r, rightOfSite = (p.x > topsite.x); if (rightOfSite && (he.side === "l")) { return 1; } if (!rightOfSite && (he.side === "r")) { return 0; } if (e.a === 1) { var dyp = p.y - topsite.y, dxp = p.x - topsite.x, fast = 0, above = 0; if ((!rightOfSite && (e.b < 0)) || (rightOfSite && (e.b >= 0))) { above = fast = (dyp >= e.b * dxp); } else { above = ((p.x + p.y * e.b) > e.c); if (e.b < 0) { above = !above; } if (!above) { fast = 1; } } if (!fast) { var dxs = topsite.x - e.region.l.x; above = (e.b * (dxp * dxp - dyp * dyp)) < (dxs * dyp * (1 + 2 * dxp / dxs + e.b * e.b)); if (e.b < 0) { above = !above; } } } else /* e.b == 1 */ { var yl = e.c - e.a * p.x, t1 = p.y - yl, t2 = p.x - topsite.x, t3 = yl - topsite.y; above = (t1 * t1) > (t2 * t2 + t3 * t3); } return he.side === "l" ? above : !above; }, endPoint: function(edge, side, site) { edge.ep[side] = site; if (!edge.ep[d3_voronoi_opposite[side]]) return; callback(edge); }, distance: function(s, t) { var dx = s.x - t.x, dy = s.y - t.y; return Math.sqrt(dx * dx + dy * dy); } }; var EventQueue = { list: [], insert: function(he, site, offset) { he.vertex = site; he.ystar = site.y + offset; for (var i=0, list=EventQueue.list, l=list.length; i next.ystar || (he.ystar == next.ystar && site.x > next.vertex.x)) { continue; } else { break; } } list.splice(i, 0, he); }, del: function(he) { for (var i=0, ls=EventQueue.list, l=ls.length; i top.y) { temp = bot; bot = top; top = temp; pm = "r"; } e = Geom.bisect(bot, top); bisector = EdgeList.createHalfEdge(e, pm); EdgeList.insert(llbnd, bisector); Geom.endPoint(e, d3_voronoi_opposite[pm], v); p = Geom.intersect(llbnd, bisector); if (p) { EventQueue.del(llbnd); EventQueue.insert(llbnd, p, Geom.distance(p, bot)); } p = Geom.intersect(bisector, rrbnd); if (p) { EventQueue.insert(bisector, p, Geom.distance(p, bot)); } } else { break; } }//end while for (lbnd = EdgeList.right(EdgeList.leftEnd); lbnd != EdgeList.rightEnd; lbnd = EdgeList.right(lbnd)) { callback(lbnd.edge); } } /** * @param vertices [[x1, y1], [x2, y2], …] * @returns triangles [[[x1, y1], [x2, y2], [x3, y3]], …] */ d3.geom.delaunay = function(vertices) { var edges = vertices.map(function() { return []; }), triangles = []; // Use the Voronoi tessellation to determine Delaunay edges. d3_voronoi_tessellate(vertices, function(e) { edges[e.region.l.index].push(vertices[e.region.r.index]); }); // Reconnect the edges into counterclockwise triangles. edges.forEach(function(edge, i) { var v = vertices[i], cx = v[0], cy = v[1]; edge.forEach(function(v) { v.angle = Math.atan2(v[0] - cx, v[1] - cy); }); edge.sort(function(a, b) { return a.angle - b.angle; }); for (var j = 0, m = edge.length - 1; j < m; j++) { triangles.push([v, edge[j], edge[j + 1]]); } }); return triangles; }; // Constructs a new quadtree for the specified array of points. A quadtree is a // two-dimensional recursive spatial subdivision. This implementation uses // square partitions, dividing each square into four equally-sized squares. Each // point exists in a unique node; if multiple points are in the same position, // some points may be stored on internal nodes rather than leaf nodes. Quadtrees // can be used to accelerate various spatial operations, such as the Barnes-Hut // approximation for computing n-body forces, or collision detection. d3.geom.quadtree = function(points, x1, y1, x2, y2) { var p, i = -1, n = points.length; // Type conversion for deprecated API. if (n && isNaN(points[0].x)) points = points.map(d3_geom_quadtreePoint); // Allow bounds to be specified explicitly. if (arguments.length < 5) { if (arguments.length === 3) { y2 = x2 = y1; y1 = x1; } else { x1 = y1 = Infinity; x2 = y2 = -Infinity; // Compute bounds. while (++i < n) { p = points[i]; if (p.x < x1) x1 = p.x; if (p.y < y1) y1 = p.y; if (p.x > x2) x2 = p.x; if (p.y > y2) y2 = p.y; } // Squarify the bounds. var dx = x2 - x1, dy = y2 - y1; if (dx > dy) y2 = y1 + dx; else x2 = x1 + dy; } } // Recursively inserts the specified point p at the node n or one of its // descendants. The bounds are defined by [x1, x2] and [y1, y2]. function insert(n, p, x1, y1, x2, y2) { if (isNaN(p.x) || isNaN(p.y)) return; // ignore invalid points if (n.leaf) { var v = n.point; if (v) { // If the point at this leaf node is at the same position as the new // point we are adding, we leave the point associated with the // internal node while adding the new point to a child node. This // avoids infinite recursion. if ((Math.abs(v.x - p.x) + Math.abs(v.y - p.y)) < .01) { insertChild(n, p, x1, y1, x2, y2); } else { n.point = null; insertChild(n, v, x1, y1, x2, y2); insertChild(n, p, x1, y1, x2, y2); } } else { n.point = p; } } else { insertChild(n, p, x1, y1, x2, y2); } } // Recursively inserts the specified point p into a descendant of node n. The // bounds are defined by [x1, x2] and [y1, y2]. function insertChild(n, p, x1, y1, x2, y2) { // Compute the split point, and the quadrant in which to insert p. var sx = (x1 + x2) * .5, sy = (y1 + y2) * .5, right = p.x >= sx, bottom = p.y >= sy, i = (bottom << 1) + right; // Recursively insert into the child node. n.leaf = false; n = n.nodes[i] || (n.nodes[i] = d3_geom_quadtreeNode()); // Update the bounds as we recurse. if (right) x1 = sx; else x2 = sx; if (bottom) y1 = sy; else y2 = sy; insert(n, p, x1, y1, x2, y2); } // Create the root node. var root = d3_geom_quadtreeNode(); root.add = function(p) { insert(root, p, x1, y1, x2, y2); }; root.visit = function(f) { d3_geom_quadtreeVisit(f, root, x1, y1, x2, y2); }; // Insert all points. points.forEach(root.add); return root; }; function d3_geom_quadtreeNode() { return { leaf: true, nodes: [], point: null }; } function d3_geom_quadtreeVisit(f, node, x1, y1, x2, y2) { if (!f(node, x1, y1, x2, y2)) { var sx = (x1 + x2) * .5, sy = (y1 + y2) * .5, children = node.nodes; if (children[0]) d3_geom_quadtreeVisit(f, children[0], x1, y1, sx, sy); if (children[1]) d3_geom_quadtreeVisit(f, children[1], sx, y1, x2, sy); if (children[2]) d3_geom_quadtreeVisit(f, children[2], x1, sy, sx, y2); if (children[3]) d3_geom_quadtreeVisit(f, children[3], sx, sy, x2, y2); } } function d3_geom_quadtreePoint(p) { return { x: p[0], y: p[1] }; } })();