ecomp-sdk-app/src/main/webapp/static/fusion/d3/js/d3.geom.js

   1 (function(){d3.geom = {};
   2 /**
   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.
   6  *
   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], …]
  11  */
  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
  21       i = 0;
  22
  23   do {
  24     // determine marching squares index
  25     i = 0;
  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;
  30
  31     // determine next direction
  32     if (i == 6) {
  33       dx = pdy == -1 ? -1 : 1;
  34       dy = 0;
  35     } else if (i == 9) {
  36       dx = 0;
  37       dy = pdx == 1 ? -1 : 1;
  38     } else {
  39       dx = d3_geom_contourDx[i];
  40       dy = d3_geom_contourDy[i];
  41     }
  42
  43     // update contour polygon
  44     if (dx != pdx && dy != pdy) {
  45       c.push([x, y]);
  46       pdx = dx;
  47       pdy = dy;
  48     }
  49
  50     x += dx;
  51     y += dy;
  52   } while (s[0] != x || s[1] != y);
  53
  54   return c;
  55 };
  56
  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];
  60
  61 function d3_geom_contourStart(grid) {
  62   var x = 0,
  63       y = 0;
  64
  65   // search for a starting point; begin at origin
  66   // and proceed along outward-expanding diagonals
  67   while (true) {
  68     if (grid(x,y)) {
  69       return [x,y];
  70     }
  71     if (x == 0) {
  72       x = y + 1;
  73       y = 0;
  74     } else {
  75       x = x - 1;
  76       y = y + 1;
  77     }
  78   }
  79 }/**
  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.
  84  *
  85  * @param vertices [[x1, y1], [x2, y2], …]
  86  * @returns polygon [[x1, y1], [x2, y2], …]
  87  */
  88 d3.geom.hull = function(vertices) {
  89   if (vertices.length < 3) return [];
  90
  91   var len = vertices.length,
  92       plen = len - 1,
  93       points = [],
  94       stack = [],
  95       i, j, h = 0, x1, y1, x2, y2, u, v, a, sp;
  96
  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]) {
 100       h = i;
 101     } else if (vertices[i][1] == vertices[h][1]) {
 102       h = (vertices[i][0] < vertices[h][0] ? i : h);
 103     }
 104   }
 105
 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});
 112   }
 113   points.sort(function(a, b) { return a.angle - b.angle; });
 114
 115   // toss out duplicate angles
 116   a = points[0].angle;
 117   v = points[0].index;
 118   u = 0;
 119   for (i=1; i<plen; ++i) {
 120     j = points[i].index;
 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;
 129       } else {
 130         points[u].index = -1;
 131         a = points[i].angle;
 132         u = i;
 133         v = j;
 134       }
 135     } else {
 136       a = points[i].angle;
 137       u = i;
 138       v = j;
 139     }
 140   }
 141
 142   // initialize the stack
 143   stack.push(h);
 144   for (i=0, j=0; i<2; ++j) {
 145     if (points[j].index != -1) {
 146       stack.push(points[j].index);
 147       i++;
 148     }
 149   }
 150   sp = stack.length;
 151
 152   // do graham's scan
 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)) {
 156       --sp;
 157     }
 158     stack[sp++] = points[j].index;
 159   }
 160
 161   // construct the hull
 162   var poly = [];
 163   for (i=0; i<sp; ++i) {
 164     poly.push(vertices[stack[i]]);
 165   }
 166   return poly;
 167 }
 168
 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) {
 178
 179   coordinates.area = function() {
 180     var i = 0,
 181         n = coordinates.length,
 182         a = coordinates[n - 1][0] * coordinates[0][1],
 183         b = coordinates[n - 1][1] * coordinates[0][0];
 184     while (++i < n) {
 185       a += coordinates[i - 1][0] * coordinates[i][1];
 186       b += coordinates[i - 1][1] * coordinates[i][0];
 187     }
 188     return (b - a) * .5;
 189   };
 190
 191   coordinates.centroid = function(k) {
 192     var i = -1,
 193         n = coordinates.length - 1,
 194         x = 0,
 195         y = 0,
 196         a,
 197         b,
 198         c;
 199     if (!arguments.length) k = 1 / (6 * coordinates.area());
 200     while (++i < n) {
 201       a = coordinates[i];
 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;
 206     }
 207     return [x * k, y * k];
 208   };
 209
 210   // The Sutherland-Hodgman clipping algorithm.
 211   coordinates.clip = function(subject) {
 212     var input,
 213         i = -1,
 214         n = coordinates.length,
 215         j,
 216         m,
 217         a = coordinates[n - 1],
 218         b,
 219         c,
 220         d;
 221     while (++i < n) {
 222       input = subject.slice();
 223       subject.length = 0;
 224       b = coordinates[i];
 225       c = input[(m = input.length) - 1];
 226       j = -1;
 227       while (++j < m) {
 228         d = input[j];
 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));
 232           }
 233           subject.push(d);
 234         } else if (d3_geom_polygonInside(c, a, b)) {
 235           subject.push(d3_geom_polygonIntersect(c, d, a, b));
 236         }
 237         c = d;
 238       }
 239       a = b;
 240     }
 241     return subject;
 242   };
 243
 244   return coordinates;
 245 };
 246
 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]);
 249 }
 250
 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],
 255       x13 = x1 - x3,
 256       x21 = x2 - x1,
 257       x43 = x4 - x3,
 258       y13 = y1 - y3,
 259       y21 = y2 - y1,
 260       y43 = y4 - y3,
 261       ua = (x43 * y13 - y43 * x13) / (y43 * x21 - x43 * y21);
 262   return [x1 + ua * x21, y1 + ua * y21];
 263 }
 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.
 268
 269 /**
 270  * @param vertices [[x1, y1], [x2, y2], …]
 271  * @returns polygons [[[x1, y1], [x2, y2], …], …]
 272  */
 273 d3.geom.voronoi = function(vertices) {
 274   var polygons = vertices.map(function() { return []; });
 275
 276   // Note: we expect the caller to clip the polygons, if needed.
 277   d3_voronoi_tessellate(vertices, function(e) {
 278     var s1,
 279         s2,
 280         x1,
 281         x2,
 282         y1,
 283         y2;
 284     if (e.a == 1 && e.b >= 0) {
 285       s1 = e.ep.r;
 286       s2 = e.ep.l;
 287     } else {
 288       s1 = e.ep.l;
 289       s2 = e.ep.r;
 290     }
 291     if (e.a == 1) {
 292       y1 = s1 ? s1.y : -1e6;
 293       x1 = e.c - e.b * y1;
 294       y2 = s2 ? s2.y : 1e6;
 295       x2 = e.c - e.b * y2;
 296     } else {
 297       x1 = s1 ? s1.x : -1e6;
 298       y1 = e.c - e.a * x1;
 299       x2 = s2 ? s2.x : 1e6;
 300       y2 = e.c - e.a * x2;
 301     }
 302     var v1 = [x1, y1],
 303         v2 = [x2, y2];
 304     polygons[e.region.l.index].push(v1, v2);
 305     polygons[e.region.r.index].push(v1, v2);
 306   });
 307
 308   // Reconnect the polygon segments into counterclockwise loops.
 309   return polygons.map(function(polygon, i) {
 310     var cx = vertices[i][0],
 311         cy = vertices[i][1];
 312     polygon.forEach(function(v) {
 313       v.angle = Math.atan2(v[0] - cx, v[1] - cy);
 314     });
 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);
 319     });
 320   });
 321 };
 322
 323 var d3_voronoi_opposite = {"l": "r", "r": "l"};
 324
 325 function d3_voronoi_tessellate(vertices, callback) {
 326
 327   var Sites = {
 328     list: vertices
 329       .map(function(v, i) {
 330         return {
 331           index: i,
 332           x: v[0],
 333           y: v[1]
 334         };
 335       })
 336       .sort(function(a, b) {
 337         return a.y < b.y ? -1
 338           : a.y > b.y ? 1
 339           : a.x < b.x ? -1
 340           : a.x > b.x ? 1
 341           : 0;
 342       }),
 343     bottomSite: null
 344   };
 345
 346   var EdgeList = {
 347     list: [],
 348     leftEnd: null,
 349     rightEnd: null,
 350
 351     init: function() {
 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);
 357     },
 358
 359     createHalfEdge: function(edge, side) {
 360       return {
 361         edge: edge,
 362         side: side,
 363         vertex: null,
 364         "l": null,
 365         "r": null
 366       };
 367     },
 368
 369     insert: function(lb, he) {
 370       he.l = lb;
 371       he.r = lb.r;
 372       lb.r.l = he;
 373       lb.r = he;
 374     },
 375
 376     leftBound: function(p) {
 377       var he = EdgeList.leftEnd;
 378       do {
 379         he = he.r;
 380       } while (he != EdgeList.rightEnd && Geom.rightOf(he, p));
 381       he = he.l;
 382       return he;
 383     },
 384
 385     del: function(he) {
 386       he.l.r = he.r;
 387       he.r.l = he.l;
 388       he.edge = null;
 389     },
 390
 391     right: function(he) {
 392       return he.r;
 393     },
 394
 395     left: function(he) {
 396       return he.l;
 397     },
 398
 399     leftRegion: function(he) {
 400       return he.edge == null
 401           ? Sites.bottomSite
 402           : he.edge.region[he.side];
 403     },
 404
 405     rightRegion: function(he) {
 406       return he.edge == null
 407           ? Sites.bottomSite
 408           : he.edge.region[d3_voronoi_opposite[he.side]];
 409     }
 410   };
 411
 412   var Geom = {
 413
 414     bisect: function(s1, s2) {
 415       var newEdge = {
 416         region: {"l": s1, "r": s2},
 417         ep: {"l": null, "r": null}
 418       };
 419
 420       var dx = s2.x - s1.x,
 421           dy = s2.y - s1.y,
 422           adx = dx > 0 ? dx : -dx,
 423           ady = dy > 0 ? dy : -dy;
 424
 425       newEdge.c = s1.x * dx + s1.y * dy
 426           + (dx * dx + dy * dy) * .5;
 427
 428       if (adx > ady) {
 429         newEdge.a = 1;
 430         newEdge.b = dy / dx;
 431         newEdge.c /= dx;
 432       } else {
 433         newEdge.b = 1;
 434         newEdge.a = dx / dy;
 435         newEdge.c /= dy;
 436       }
 437
 438       return newEdge;
 439     },
 440
 441     intersect: function(el1, el2) {
 442       var e1 = el1.edge,
 443           e2 = el2.edge;
 444       if (!e1 || !e2 || (e1.region.r == e2.region.r)) {
 445         return null;
 446       }
 447       var d = (e1.a * e2.b) - (e1.b * e2.a);
 448       if (Math.abs(d) < 1e-10) {
 449         return null;
 450       }
 451       var xint = (e1.c * e2.b - e2.c * e1.b) / d,
 452           yint = (e2.c * e1.a - e1.c * e2.a) / d,
 453           e1r = e1.region.r,
 454           e2r = e2.region.r,
 455           el,
 456           e;
 457       if ((e1r.y < e2r.y) ||
 458          (e1r.y == e2r.y && e1r.x < e2r.x)) {
 459         el = el1;
 460         e = e1;
 461       } else {
 462         el = el2;
 463         e = e2;
 464       }
 465       var rightOfSite = (xint >= e.region.r.x);
 466       if ((rightOfSite && (el.side == "l")) ||
 467         (!rightOfSite && (el.side == "r"))) {
 468         return null;
 469       }
 470       return {
 471         x: xint,
 472         y: yint
 473       };
 474     },
 475
 476     rightOf: function(he, p) {
 477       var e = he.edge,
 478           topsite = e.region.r,
 479           rightOfSite = (p.x > topsite.x);
 480
 481       if (rightOfSite && (he.side == "l")) {
 482         return 1;
 483       }
 484       if (!rightOfSite && (he.side == "r")) {
 485         return 0;
 486       }
 487       if (e.a == 1) {
 488         var dyp = p.y - topsite.y,
 489             dxp = p.x - topsite.x,
 490             fast = 0,
 491             above = 0;
 492
 493         if ((!rightOfSite && (e.b < 0)) ||
 494           (rightOfSite && (e.b >= 0))) {
 495           above = fast = (dyp >= e.b * dxp);
 496         } else {
 497           above = ((p.x + p.y * e.b) > e.c);
 498           if (e.b < 0) {
 499             above = !above;
 500           }
 501           if (!above) {
 502             fast = 1;
 503           }
 504         }
 505         if (!fast) {
 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));
 509
 510           if (e.b < 0) {
 511             above = !above;
 512           }
 513         }
 514       } else /* e.b == 1 */ {
 515         var yl = e.c - e.a * p.x,
 516             t1 = p.y - yl,
 517             t2 = p.x - topsite.x,
 518             t3 = yl - topsite.y;
 519
 520         above = (t1 * t1) > (t2 * t2 + t3 * t3);
 521       }
 522       return he.side == "l" ? above : !above;
 523     },
 524
 525     endPoint: function(edge, side, site) {
 526       edge.ep[side] = site;
 527       if (!edge.ep[d3_voronoi_opposite[side]]) return;
 528       callback(edge);
 529     },
 530
 531     distance: function(s, t) {
 532       var dx = s.x - t.x,
 533           dy = s.y - t.y;
 534       return Math.sqrt(dx * dx + dy * dy);
 535     }
 536   };
 537
 538   var EventQueue = {
 539     list: [],
 540
 541     insert: function(he, site, offset) {
 542       he.vertex = site;
 543       he.ystar = site.y + offset;
 544       for (var i=0, list=EventQueue.list, l=list.length; i<l; i++) {
 545         var next = list[i];
 546         if (he.ystar > next.ystar ||
 547           (he.ystar == next.ystar &&
 548           site.x > next.vertex.x)) {
 549           continue;
 550         } else {
 551           break;
 552         }
 553       }
 554       list.splice(i, 0, he);
 555     },
 556
 557     del: function(he) {
 558       for (var i=0, ls=EventQueue.list, l=ls.length; i<l && (ls[i] != he); ++i) {}
 559       ls.splice(i, 1);
 560     },
 561
 562     empty: function() { return EventQueue.list.length == 0; },
 563
 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];
 567       }
 568       return null;
 569     },
 570
 571     min: function() {
 572       var elem = EventQueue.list[0];
 573       return {
 574         x: elem.vertex.x,
 575         y: elem.ystar
 576       };
 577     },
 578
 579     extractMin: function() {
 580       return EventQueue.list.shift();
 581     }
 582   };
 583
 584   EdgeList.init();
 585   Sites.bottomSite = Sites.list.shift();
 586
 587   var newSite = Sites.list.shift(), newIntStar;
 588   var lbnd, rbnd, llbnd, rrbnd, bisector;
 589   var bot, top, temp, p, v;
 590   var e, pm;
 591
 592   while (true) {
 593     if (!EventQueue.empty()) {
 594       newIntStar = EventQueue.min();
 595     }
 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);
 607       if (p) {
 608         EventQueue.del(lbnd);
 609         EventQueue.insert(lbnd, p, Geom.distance(p, newSite));
 610       }
 611       lbnd = bisector;
 612       bisector = EdgeList.createHalfEdge(e, "r");
 613       EdgeList.insert(lbnd, bisector);
 614       p = Geom.intersect(bisector, rbnd);
 615       if (p) {
 616         EventQueue.insert(bisector, p, Geom.distance(p, newSite));
 617       }
 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);
 626       v = lbnd.vertex;
 627       Geom.endPoint(lbnd.edge, lbnd.side, v);
 628       Geom.endPoint(rbnd.edge, rbnd.side, v);
 629       EdgeList.del(lbnd);
 630       EventQueue.del(rbnd);
 631       EdgeList.del(rbnd);
 632       pm = "l";
 633       if (bot.y > top.y) {
 634         temp = bot;
 635         bot = top;
 636         top = temp;
 637         pm = "r";
 638       }
 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);
 644       if (p) {
 645         EventQueue.del(llbnd);
 646         EventQueue.insert(llbnd, p, Geom.distance(p, bot));
 647       }
 648       p = Geom.intersect(bisector, rrbnd);
 649       if (p) {
 650         EventQueue.insert(bisector, p, Geom.distance(p, bot));
 651       }
 652     } else {
 653       break;
 654     }
 655   }//end while
 656
 657   for (lbnd = EdgeList.right(EdgeList.leftEnd);
 658       lbnd != EdgeList.rightEnd;
 659       lbnd = EdgeList.right(lbnd)) {
 660     callback(lbnd.edge);
 661   }
 662 }
 663 /**
 664 * @param vertices [[x1, y1], [x2, y2], …]
 665 * @returns triangles [[[x1, y1], [x2, y2], [x3, y3]], …]
 666  */
 667 d3.geom.delaunay = function(vertices) {
 668   var edges = vertices.map(function() { return []; }),
 669       triangles = [];
 670
 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]);
 674   });
 675
 676   // Reconnect the edges into counterclockwise triangles.
 677   edges.forEach(function(edge, i) {
 678     var v = vertices[i],
 679         cx = v[0],
 680         cy = v[1];
 681     edge.forEach(function(v) {
 682       v.angle = Math.atan2(v[0] - cx, v[1] - cy);
 683     });
 684     edge.sort(function(a, b) {
 685       return a.angle - b.angle;
 686     });
 687     for (var j = 0, m = edge.length - 1; j < m; j++) {
 688       triangles.push([v, edge[j], edge[j + 1]]);
 689     }
 690   });
 691
 692   return triangles;
 693 };
 694 /**
 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.
 702  *
 703  * @param points [[x1, y1], [x2, y2], …]
 704  * @return quadtree root {left: boolean, nodes: […], point: [x, y]}
 705  */
 706 d3.geom.quadtree = function(points) {
 707   var p,
 708       i = -1,
 709       n = points.length;
 710
 711   /* Compute bounds. */
 712   var x1 = Number.POSITIVE_INFINITY, y1 = x1,
 713       x2 = Number.NEGATIVE_INFINITY, y2 = x2;
 714   while (++i < n) {
 715     p = points[i];
 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];
 720   }
 721
 722   /* Squarify the bounds. */
 723   var dx = x2 - x1,
 724       dy = y2 - y1;
 725   if (dx > dy) y2 = y1 + dx;
 726   else x2 = x1 + dy;
 727
 728   /**
 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>].
 732    */
 733   function insert(n, p, x1, y1, x2, y2) {
 734     if (isNaN(p[0]) || isNaN(p[1])) return; // ignore invalid points
 735     if (n.leaf) {
 736       var v = n.point;
 737       if (v) {
 738         /*
 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.
 743          */
 744         if ((Math.abs(v[0] - p[0]) + Math.abs(v[1] - p[1])) < .01) {
 745           insertChild(n, p, x1, y1, x2, y2);
 746         } else {
 747           n.point = null;
 748           insertChild(n, v, x1, y1, x2, y2);
 749           insertChild(n, p, x1, y1, x2, y2);
 750         }
 751       } else {
 752         n.point = p;
 753       }
 754     } else {
 755       insertChild(n, p, x1, y1, x2, y2);
 756     }
 757   }
 758
 759   /**
 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>].
 763    */
 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,
 767         sy = (y1 + y2) * .5,
 768         right = p[0] >= sx,
 769         bottom = p[1] >= sy,
 770         i = (bottom << 1) + right;
 771
 772     /* Recursively insert into the child node. */
 773     n.leaf = false;
 774     n = n.nodes[i] || (n.nodes[i] = d3_geom_quadtreeNode());
 775
 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);
 780   }
 781
 782   /* Create the root node. */
 783   var root = d3_geom_quadtreeNode();
 784
 785   /* Insert all points. */
 786   i = -1;
 787   while (++i < n) insert(root, points[i], x1, y1, x2, y2);
 788
 789   /* Register a visitor function for the root. */
 790   root.visit = function(f) {
 791     d3_geom_quadtreeVisit(f, root, x1, y1, x2, y2);
 792   };
 793
 794   return root;
 795 };
 796
 797 function d3_geom_quadtreeNode() {
 798   return {
 799     leaf: true,
 800     nodes: [],
 801     point: null
 802   };
 803 }
 804
 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,
 808         sy = (y1 + y2) * .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);
 814   }
 815 }
 816 })()