1 (function(){d3.geom = {};
3 * Computes a contour for a given input grid function using the <a
4 * href="http://en.wikipedia.org/wiki/Marching_squares">marching
5 * squares</a> algorithm. Returns the contour polygon as an array of points.
7 * @param grid a two-input function(x, y) that returns true for values
8 * inside the contour and false for values outside the contour.
9 * @param start an optional starting point [x, y] on the grid.
10 * @returns polygon [[x1, y1], [x2, y2], …]
12 d3.geom.contour = function(grid, start) {
13 var s = start || d3_geom_contourStart(grid), // starting point
14 c = [], // contour polygon
15 x = s[0], // current x position
16 y = s[1], // current y position
17 dx = 0, // next x direction
18 dy = 0, // next y direction
19 pdx = NaN, // previous x direction
20 pdy = NaN, // previous y direction
24 // determine marching squares index
26 if (grid(x-1, y-1)) i += 1;
27 if (grid(x, y-1)) i += 2;
28 if (grid(x-1, y )) i += 4;
29 if (grid(x, y )) i += 8;
31 // determine next direction
33 dx = pdy == -1 ? -1 : 1;
37 dy = pdx == 1 ? -1 : 1;
39 dx = d3_geom_contourDx[i];
40 dy = d3_geom_contourDy[i];
43 // update contour polygon
44 if (dx != pdx && dy != pdy) {
52 } while (s[0] != x || s[1] != y);
57 // lookup tables for marching directions
58 var d3_geom_contourDx = [1, 0, 1, 1,-1, 0,-1, 1,0, 0,0,0,-1, 0,-1,NaN],
59 d3_geom_contourDy = [0,-1, 0, 0, 0,-1, 0, 0,1,-1,1,1, 0,-1, 0,NaN];
61 function d3_geom_contourStart(grid) {
65 // search for a starting point; begin at origin
66 // and proceed along outward-expanding diagonals
80 * Computes the 2D convex hull of a set of points using Graham's scanning
81 * algorithm. The algorithm has been implemented as described in Cormen,
82 * Leiserson, and Rivest's Introduction to Algorithms. The running time of
83 * this algorithm is O(n log n), where n is the number of input points.
85 * @param vertices [[x1, y1], [x2, y2], …]
86 * @returns polygon [[x1, y1], [x2, y2], …]
88 d3.geom.hull = function(vertices) {
89 if (vertices.length < 3) return [];
91 var len = vertices.length,
95 i, j, h = 0, x1, y1, x2, y2, u, v, a, sp;
97 // find the starting ref point: leftmost point with the minimum y coord
98 for (i=1; i<len; ++i) {
99 if (vertices[i][1] < vertices[h][1]) {
101 } else if (vertices[i][1] == vertices[h][1]) {
102 h = (vertices[i][0] < vertices[h][0] ? i : h);
106 // calculate polar angles from ref point and sort
107 for (i=0; i<len; ++i) {
108 if (i == h) continue;
109 y1 = vertices[i][1] - vertices[h][1];
110 x1 = vertices[i][0] - vertices[h][0];
111 points.push({angle: Math.atan2(y1, x1), index: i});
113 points.sort(function(a, b) { return a.angle - b.angle; });
115 // toss out duplicate angles
119 for (i=1; i<plen; ++i) {
121 if (a == points[i].angle) {
122 // keep angle for point most distant from the reference
123 x1 = vertices[v][0] - vertices[h][0];
124 y1 = vertices[v][1] - vertices[h][1];
125 x2 = vertices[j][0] - vertices[h][0];
126 y2 = vertices[j][1] - vertices[h][1];
127 if ((x1*x1 + y1*y1) >= (x2*x2 + y2*y2)) {
128 points[i].index = -1;
130 points[u].index = -1;
142 // initialize the stack
144 for (i=0, j=0; i<2; ++j) {
145 if (points[j].index != -1) {
146 stack.push(points[j].index);
153 for (; j<plen; ++j) {
154 if (points[j].index == -1) continue; // skip tossed out points
155 while (!d3_geom_hullCCW(stack[sp-2], stack[sp-1], points[j].index, vertices)) {
158 stack[sp++] = points[j].index;
161 // construct the hull
163 for (i=0; i<sp; ++i) {
164 poly.push(vertices[stack[i]]);
169 // are three points in counter-clockwise order?
170 function d3_geom_hullCCW(i1, i2, i3, v) {
171 var t, a, b, c, d, e, f;
172 t = v[i1]; a = t[0]; b = t[1];
173 t = v[i2]; c = t[0]; d = t[1];
174 t = v[i3]; e = t[0]; f = t[1];
175 return ((f-b)*(c-a) - (d-b)*(e-a)) > 0;
176 }// Note: requires coordinates to be counterclockwise and convex!
177 d3.geom.polygon = function(coordinates) {
179 coordinates.area = function() {
181 n = coordinates.length,
182 a = coordinates[n - 1][0] * coordinates[0][1],
183 b = coordinates[n - 1][1] * coordinates[0][0];
185 a += coordinates[i - 1][0] * coordinates[i][1];
186 b += coordinates[i - 1][1] * coordinates[i][0];
191 coordinates.centroid = function(k) {
193 n = coordinates.length - 1,
199 if (!arguments.length) k = 1 / (6 * coordinates.area());
202 b = coordinates[i + 1];
203 c = a[0] * b[1] - b[0] * a[1];
204 x += (a[0] + b[0]) * c;
205 y += (a[1] + b[1]) * c;
207 return [x * k, y * k];
210 // The Sutherland-Hodgman clipping algorithm.
211 coordinates.clip = function(subject) {
214 n = coordinates.length,
217 a = coordinates[n - 1],
222 input = subject.slice();
225 c = input[(m = input.length) - 1];
229 if (d3_geom_polygonInside(d, a, b)) {
230 if (!d3_geom_polygonInside(c, a, b)) {
231 subject.push(d3_geom_polygonIntersect(c, d, a, b));
234 } else if (d3_geom_polygonInside(c, a, b)) {
235 subject.push(d3_geom_polygonIntersect(c, d, a, b));
247 function d3_geom_polygonInside(p, a, b) {
248 return (b[0] - a[0]) * (p[1] - a[1]) < (b[1] - a[1]) * (p[0] - a[0]);
251 // Intersect two infinite lines cd and ab.
252 function d3_geom_polygonIntersect(c, d, a, b) {
253 var x1 = c[0], x2 = d[0], x3 = a[0], x4 = b[0],
254 y1 = c[1], y2 = d[1], y3 = a[1], y4 = b[1],
261 ua = (x43 * y13 - y43 * x13) / (y43 * x21 - x43 * y21);
262 return [x1 + ua * x21, y1 + ua * y21];
264 // Adapted from Nicolas Garcia Belmonte's JIT implementation:
265 // http://blog.thejit.org/2010/02/12/voronoi-tessellation/
266 // http://blog.thejit.org/assets/voronoijs/voronoi.js
267 // See lib/jit/LICENSE for details.
270 * @param vertices [[x1, y1], [x2, y2], …]
271 * @returns polygons [[[x1, y1], [x2, y2], …], …]
273 d3.geom.voronoi = function(vertices) {
274 var polygons = vertices.map(function() { return []; });
276 // Note: we expect the caller to clip the polygons, if needed.
277 d3_voronoi_tessellate(vertices, function(e) {
284 if (e.a == 1 && e.b >= 0) {
292 y1 = s1 ? s1.y : -1e6;
294 y2 = s2 ? s2.y : 1e6;
297 x1 = s1 ? s1.x : -1e6;
299 x2 = s2 ? s2.x : 1e6;
304 polygons[e.region.l.index].push(v1, v2);
305 polygons[e.region.r.index].push(v1, v2);
308 // Reconnect the polygon segments into counterclockwise loops.
309 return polygons.map(function(polygon, i) {
310 var cx = vertices[i][0],
312 polygon.forEach(function(v) {
313 v.angle = Math.atan2(v[0] - cx, v[1] - cy);
315 return polygon.sort(function(a, b) {
316 return a.angle - b.angle;
317 }).filter(function(d, i) {
318 return !i || (d.angle - polygon[i - 1].angle > 1e-10);
323 var d3_voronoi_opposite = {"l": "r", "r": "l"};
325 function d3_voronoi_tessellate(vertices, callback) {
329 .map(function(v, i) {
336 .sort(function(a, b) {
337 return a.y < b.y ? -1
352 EdgeList.leftEnd = EdgeList.createHalfEdge(null, "l");
353 EdgeList.rightEnd = EdgeList.createHalfEdge(null, "l");
354 EdgeList.leftEnd.r = EdgeList.rightEnd;
355 EdgeList.rightEnd.l = EdgeList.leftEnd;
356 EdgeList.list.unshift(EdgeList.leftEnd, EdgeList.rightEnd);
359 createHalfEdge: function(edge, side) {
369 insert: function(lb, he) {
376 leftBound: function(p) {
377 var he = EdgeList.leftEnd;
380 } while (he != EdgeList.rightEnd && Geom.rightOf(he, p));
391 right: function(he) {
399 leftRegion: function(he) {
400 return he.edge == null
402 : he.edge.region[he.side];
405 rightRegion: function(he) {
406 return he.edge == null
408 : he.edge.region[d3_voronoi_opposite[he.side]];
414 bisect: function(s1, s2) {
416 region: {"l": s1, "r": s2},
417 ep: {"l": null, "r": null}
420 var dx = s2.x - s1.x,
422 adx = dx > 0 ? dx : -dx,
423 ady = dy > 0 ? dy : -dy;
425 newEdge.c = s1.x * dx + s1.y * dy
426 + (dx * dx + dy * dy) * .5;
441 intersect: function(el1, el2) {
444 if (!e1 || !e2 || (e1.region.r == e2.region.r)) {
447 var d = (e1.a * e2.b) - (e1.b * e2.a);
448 if (Math.abs(d) < 1e-10) {
451 var xint = (e1.c * e2.b - e2.c * e1.b) / d,
452 yint = (e2.c * e1.a - e1.c * e2.a) / d,
457 if ((e1r.y < e2r.y) ||
458 (e1r.y == e2r.y && e1r.x < e2r.x)) {
465 var rightOfSite = (xint >= e.region.r.x);
466 if ((rightOfSite && (el.side == "l")) ||
467 (!rightOfSite && (el.side == "r"))) {
476 rightOf: function(he, p) {
478 topsite = e.region.r,
479 rightOfSite = (p.x > topsite.x);
481 if (rightOfSite && (he.side == "l")) {
484 if (!rightOfSite && (he.side == "r")) {
488 var dyp = p.y - topsite.y,
489 dxp = p.x - topsite.x,
493 if ((!rightOfSite && (e.b < 0)) ||
494 (rightOfSite && (e.b >= 0))) {
495 above = fast = (dyp >= e.b * dxp);
497 above = ((p.x + p.y * e.b) > e.c);
506 var dxs = topsite.x - e.region.l.x;
507 above = (e.b * (dxp * dxp - dyp * dyp)) <
508 (dxs * dyp * (1 + 2 * dxp / dxs + e.b * e.b));
514 } else /* e.b == 1 */ {
515 var yl = e.c - e.a * p.x,
517 t2 = p.x - topsite.x,
520 above = (t1 * t1) > (t2 * t2 + t3 * t3);
522 return he.side == "l" ? above : !above;
525 endPoint: function(edge, side, site) {
526 edge.ep[side] = site;
527 if (!edge.ep[d3_voronoi_opposite[side]]) return;
531 distance: function(s, t) {
534 return Math.sqrt(dx * dx + dy * dy);
541 insert: function(he, site, offset) {
543 he.ystar = site.y + offset;
544 for (var i=0, list=EventQueue.list, l=list.length; i<l; i++) {
546 if (he.ystar > next.ystar ||
547 (he.ystar == next.ystar &&
548 site.x > next.vertex.x)) {
554 list.splice(i, 0, he);
558 for (var i=0, ls=EventQueue.list, l=ls.length; i<l && (ls[i] != he); ++i) {}
562 empty: function() { return EventQueue.list.length == 0; },
564 nextEvent: function(he) {
565 for (var i=0, ls=EventQueue.list, l=ls.length; i<l; ++i) {
566 if (ls[i] == he) return ls[i+1];
572 var elem = EventQueue.list[0];
579 extractMin: function() {
580 return EventQueue.list.shift();
585 Sites.bottomSite = Sites.list.shift();
587 var newSite = Sites.list.shift(), newIntStar;
588 var lbnd, rbnd, llbnd, rrbnd, bisector;
589 var bot, top, temp, p, v;
593 if (!EventQueue.empty()) {
594 newIntStar = EventQueue.min();
596 if (newSite && (EventQueue.empty()
597 || newSite.y < newIntStar.y
598 || (newSite.y == newIntStar.y
599 && newSite.x < newIntStar.x))) { //new site is smallest
600 lbnd = EdgeList.leftBound(newSite);
601 rbnd = EdgeList.right(lbnd);
602 bot = EdgeList.rightRegion(lbnd);
603 e = Geom.bisect(bot, newSite);
604 bisector = EdgeList.createHalfEdge(e, "l");
605 EdgeList.insert(lbnd, bisector);
606 p = Geom.intersect(lbnd, bisector);
608 EventQueue.del(lbnd);
609 EventQueue.insert(lbnd, p, Geom.distance(p, newSite));
612 bisector = EdgeList.createHalfEdge(e, "r");
613 EdgeList.insert(lbnd, bisector);
614 p = Geom.intersect(bisector, rbnd);
616 EventQueue.insert(bisector, p, Geom.distance(p, newSite));
618 newSite = Sites.list.shift();
619 } else if (!EventQueue.empty()) { //intersection is smallest
620 lbnd = EventQueue.extractMin();
621 llbnd = EdgeList.left(lbnd);
622 rbnd = EdgeList.right(lbnd);
623 rrbnd = EdgeList.right(rbnd);
624 bot = EdgeList.leftRegion(lbnd);
625 top = EdgeList.rightRegion(rbnd);
627 Geom.endPoint(lbnd.edge, lbnd.side, v);
628 Geom.endPoint(rbnd.edge, rbnd.side, v);
630 EventQueue.del(rbnd);
639 e = Geom.bisect(bot, top);
640 bisector = EdgeList.createHalfEdge(e, pm);
641 EdgeList.insert(llbnd, bisector);
642 Geom.endPoint(e, d3_voronoi_opposite[pm], v);
643 p = Geom.intersect(llbnd, bisector);
645 EventQueue.del(llbnd);
646 EventQueue.insert(llbnd, p, Geom.distance(p, bot));
648 p = Geom.intersect(bisector, rrbnd);
650 EventQueue.insert(bisector, p, Geom.distance(p, bot));
657 for (lbnd = EdgeList.right(EdgeList.leftEnd);
658 lbnd != EdgeList.rightEnd;
659 lbnd = EdgeList.right(lbnd)) {
664 * @param vertices [[x1, y1], [x2, y2], …]
665 * @returns triangles [[[x1, y1], [x2, y2], [x3, y3]], …]
667 d3.geom.delaunay = function(vertices) {
668 var edges = vertices.map(function() { return []; }),
671 // Use the Voronoi tessellation to determine Delaunay edges.
672 d3_voronoi_tessellate(vertices, function(e) {
673 edges[e.region.l.index].push(vertices[e.region.r.index]);
676 // Reconnect the edges into counterclockwise triangles.
677 edges.forEach(function(edge, i) {
681 edge.forEach(function(v) {
682 v.angle = Math.atan2(v[0] - cx, v[1] - cy);
684 edge.sort(function(a, b) {
685 return a.angle - b.angle;
687 for (var j = 0, m = edge.length - 1; j < m; j++) {
688 triangles.push([v, edge[j], edge[j + 1]]);
695 * Constructs a new quadtree for the specified array of points. A quadtree is a
696 * two-dimensional recursive spatial subdivision. This implementation uses
697 * square partitions, dividing each square into four equally-sized squares. Each
698 * point exists in a unique node; if multiple points are in the same position,
699 * some points may be stored on internal nodes rather than leaf nodes. Quadtrees
700 * can be used to accelerate various spatial operations, such as the Barnes-Hut
701 * approximation for computing n-body forces, or collision detection.
703 * @param points [[x1, y1], [x2, y2], …]
704 * @return quadtree root {left: boolean, nodes: […], point: [x, y]}
706 d3.geom.quadtree = function(points) {
711 /* Compute bounds. */
712 var x1 = Number.POSITIVE_INFINITY, y1 = x1,
713 x2 = Number.NEGATIVE_INFINITY, y2 = x2;
716 if (p[0] < x1) x1 = p[0];
717 if (p[1] < y1) y1 = p[1];
718 if (p[0] > x2) x2 = p[0];
719 if (p[1] > y2) y2 = p[1];
722 /* Squarify the bounds. */
725 if (dx > dy) y2 = y1 + dx;
729 * @ignore Recursively inserts the specified point <i>p</i> at the node
730 * <i>n</i> or one of its descendants. The bounds are defined by [<i>x1</i>,
731 * <i>x2</i>] and [<i>y1</i>, <i>y2</i>].
733 function insert(n, p, x1, y1, x2, y2) {
734 if (isNaN(p[0]) || isNaN(p[1])) return; // ignore invalid points
739 * If the point at this leaf node is at the same position as the new
740 * point we are adding, we leave the point associated with the
741 * internal node while adding the new point to a child node. This
742 * avoids infinite recursion.
744 if ((Math.abs(v[0] - p[0]) + Math.abs(v[1] - p[1])) < .01) {
745 insertChild(n, p, x1, y1, x2, y2);
748 insertChild(n, v, x1, y1, x2, y2);
749 insertChild(n, p, x1, y1, x2, y2);
755 insertChild(n, p, x1, y1, x2, y2);
760 * @ignore Recursively inserts the specified point <i>p</i> into a
761 * descendant of node <i>n</i>. The bounds are defined by [<i>x1</i>,
762 * <i>x2</i>] and [<i>y1</i>, <i>y2</i>].
764 function insertChild(n, p, x1, y1, x2, y2) {
765 /* Compute the split point, and the quadrant in which to insert p. */
766 var sx = (x1 + x2) * .5,
770 i = (bottom << 1) + right;
772 /* Recursively insert into the child node. */
774 n = n.nodes[i] || (n.nodes[i] = d3_geom_quadtreeNode());
776 /* Update the bounds as we recurse. */
777 if (right) x1 = sx; else x2 = sx;
778 if (bottom) y1 = sy; else y2 = sy;
779 insert(n, p, x1, y1, x2, y2);
782 /* Create the root node. */
783 var root = d3_geom_quadtreeNode();
785 /* Insert all points. */
787 while (++i < n) insert(root, points[i], x1, y1, x2, y2);
789 /* Register a visitor function for the root. */
790 root.visit = function(f) {
791 d3_geom_quadtreeVisit(f, root, x1, y1, x2, y2);
797 function d3_geom_quadtreeNode() {
805 function d3_geom_quadtreeVisit(f, node, x1, y1, x2, y2) {
806 if (!f(node, x1, y1, x2, y2)) {
807 var sx = (x1 + x2) * .5,
809 children = node.nodes;
810 if (children[0]) d3_geom_quadtreeVisit(f, children[0], x1, y1, sx, sy);
811 if (children[1]) d3_geom_quadtreeVisit(f, children[1], sx, y1, x2, sy);
812 if (children[2]) d3_geom_quadtreeVisit(f, children[2], x1, sy, sx, y2);
813 if (children[3]) d3_geom_quadtreeVisit(f, children[3], sx, sy, x2, y2);