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 /**

    81  * Computes the 2D convex hull of a set of points using Graham's scanning

    82  * algorithm. The algorithm has been implemented as described in Cormen,

    83  * Leiserson, and Rivest's Introduction to Algorithms. The running time of

    84  * this algorithm is O(n log n), where n is the number of input points.

    85  *

    86  * @param vertices [[x1, y1], [x2, y2], …]

    87  * @returns polygon [[x1, y1], [x2, y2], …]

    88  */

    89 d3.geom.hull = function(vertices) {

    90   if (vertices.length < 3) return [];

    91 

    92   var len = vertices.length,

    93       plen = len - 1,

    94       points = [],

    95       stack = [],

    96       i, j, h = 0, x1, y1, x2, y2, u, v, a, sp;

    97 

    98   // find the starting ref point: leftmost point with the minimum y coord

    99   for (i=1; i<len; ++i) {

   100     if (vertices[i][1] < vertices[h][1]) {

   101       h = i;

   102     } else if (vertices[i][1] == vertices[h][1]) {

   103       h = (vertices[i][0] < vertices[h][0] ? i : h);

   104     }

   105   }

   106 

   107   // calculate polar angles from ref point and sort

   108   for (i=0; i<len; ++i) {

   109     if (i === h) continue;

   110     y1 = vertices[i][1] - vertices[h][1];

   111     x1 = vertices[i][0] - vertices[h][0];

   112     points.push({angle: Math.atan2(y1, x1), index: i});

   113   }

   114   points.sort(function(a, b) { return a.angle - b.angle; });

   115 

   116   // toss out duplicate angles

   117   a = points[0].angle;

   118   v = points[0].index;

   119   u = 0;

   120   for (i=1; i<plen; ++i) {

   121     j = points[i].index;

   122     if (a == points[i].angle) {

   123       // keep angle for point most distant from the reference

   124       x1 = vertices[v][0] - vertices[h][0];

   125       y1 = vertices[v][1] - vertices[h][1];

   126       x2 = vertices[j][0] - vertices[h][0];

   127       y2 = vertices[j][1] - vertices[h][1];

   128       if ((x1*x1 + y1*y1) >= (x2*x2 + y2*y2)) {

   129         points[i].index = -1;

   130       } else {

   131         points[u].index = -1;

   132         a = points[i].angle;

   133         u = i;

   134         v = j;

   135       }

   136     } else {

   137       a = points[i].angle;

   138       u = i;

   139       v = j;

   140     }

   141   }

   142 

   143   // initialize the stack

   144   stack.push(h);

   145   for (i=0, j=0; i<2; ++j) {

   146     if (points[j].index !== -1) {

   147       stack.push(points[j].index);

   148       i++;

   149     }

   150   }

   151   sp = stack.length;

   152 

   153   // do graham's scan

   154   for (; j<plen; ++j) {

   155     if (points[j].index === -1) continue; // skip tossed out points

   156     while (!d3_geom_hullCCW(stack[sp-2], stack[sp-1], points[j].index, vertices)) {

   157       --sp;

   158     }

   159     stack[sp++] = points[j].index;

   160   }

   161 

   162   // construct the hull

   163   var poly = [];

   164   for (i=0; i<sp; ++i) {

   165     poly.push(vertices[stack[i]]);

   166   }

   167   return poly;

   168 }

   169 

   170 // are three points in counter-clockwise order?

   171 function d3_geom_hullCCW(i1, i2, i3, v) {

   172   var t, a, b, c, d, e, f;

   173   t = v[i1]; a = t[0]; b = t[1];

   174   t = v[i2]; c = t[0]; d = t[1];

   175   t = v[i3]; e = t[0]; f = t[1];

   176   return ((f-b)*(c-a) - (d-b)*(e-a)) > 0;

   177 }

   178 // Note: requires coordinates to be counterclockwise and convex!

   179 d3.geom.polygon = function(coordinates) {

   180 

   181   coordinates.area = function() {

   182     var i = 0,

   183         n = coordinates.length,

   184         a = coordinates[n - 1][0] * coordinates[0][1],

   185         b = coordinates[n - 1][1] * coordinates[0][0];

   186     while (++i < n) {

   187       a += coordinates[i - 1][0] * coordinates[i][1];

   188       b += coordinates[i - 1][1] * coordinates[i][0];

   189     }

   190     return (b - a) * .5;

   191   };

   192 

   193   coordinates.centroid = function(k) {

   194     var i = -1,

   195         n = coordinates.length - 1,

   196         x = 0,

   197         y = 0,

   198         a,

   199         b,

   200         c;

   201     if (!arguments.length) k = -1 / (6 * coordinates.area());

   202     while (++i < n) {

   203       a = coordinates[i];

   204       b = coordinates[i + 1];

   205       c = a[0] * b[1] - b[0] * a[1];

   206       x += (a[0] + b[0]) * c;

   207       y += (a[1] + b[1]) * c;

   208     }

   209     return [x * k, y * k];

   210   };

   211 

   212   // The Sutherland-Hodgman clipping algorithm.

   213   coordinates.clip = function(subject) {

   214     var input,

   215         i = -1,

   216         n = coordinates.length,

   217         j,

   218         m,

   219         a = coordinates[n - 1],

   220         b,

   221         c,

   222         d;

   223     while (++i < n) {

   224       input = subject.slice();

   225       subject.length = 0;

   226       b = coordinates[i];

   227       c = input[(m = input.length) - 1];

   228       j = -1;

   229       while (++j < m) {

   230         d = input[j];

   231         if (d3_geom_polygonInside(d, a, b)) {

   232           if (!d3_geom_polygonInside(c, a, b)) {

   233             subject.push(d3_geom_polygonIntersect(c, d, a, b));

   234           }

   235           subject.push(d);

   236         } else if (d3_geom_polygonInside(c, a, b)) {

   237           subject.push(d3_geom_polygonIntersect(c, d, a, b));

   238         }

   239         c = d;

   240       }

   241       a = b;

   242     }

   243     return subject;

   244   };

   245 

   246   return coordinates;

   247 };

   248 

   249 function d3_geom_polygonInside(p, a, b) {

   250   return (b[0] - a[0]) * (p[1] - a[1]) < (b[1] - a[1]) * (p[0] - a[0]);

   251 }

   252 

   253 // Intersect two infinite lines cd and ab.

   254 function d3_geom_polygonIntersect(c, d, a, b) {

   255   var x1 = c[0], x2 = d[0], x3 = a[0], x4 = b[0],

   256       y1 = c[1], y2 = d[1], y3 = a[1], y4 = b[1],

   257       x13 = x1 - x3,

   258       x21 = x2 - x1,

   259       x43 = x4 - x3,

   260       y13 = y1 - y3,

   261       y21 = y2 - y1,

   262       y43 = y4 - y3,

   263       ua = (x43 * y13 - y43 * x13) / (y43 * x21 - x43 * y21);

   264   return [x1 + ua * x21, y1 + ua * y21];

   265 }

   266 // Adapted from Nicolas Garcia Belmonte's JIT implementation:

   267 // http://blog.thejit.org/2010/02/12/voronoi-tessellation/

   268 // http://blog.thejit.org/assets/voronoijs/voronoi.js

   269 // See lib/jit/LICENSE for details.

   270 

   271 // Notes:

   272 //

   273 // This implementation does not clip the returned polygons, so if you want to

   274 // clip them to a particular shape you will need to do that either in SVG or by

   275 // post-processing with d3.geom.polygon's clip method.

   276 //

   277 // If any vertices are coincident or have NaN positions, the behavior of this

   278 // method is undefined. Most likely invalid polygons will be returned. You

   279 // should filter invalid points, and consolidate coincident points, before

   280 // computing the tessellation.

   281 

   282 /**

   283  * @param vertices [[x1, y1], [x2, y2], …]

   284  * @returns polygons [[[x1, y1], [x2, y2], …], …]

   285  */

   286 d3.geom.voronoi = function(vertices) {

   287   var polygons = vertices.map(function() { return []; });

   288 

   289   d3_voronoi_tessellate(vertices, function(e) {

   290     var s1,

   291         s2,

   292         x1,

   293         x2,

   294         y1,

   295         y2;

   296     if (e.a === 1 && e.b >= 0) {

   297       s1 = e.ep.r;

   298       s2 = e.ep.l;

   299     } else {

   300       s1 = e.ep.l;

   301       s2 = e.ep.r;

   302     }

   303     if (e.a === 1) {

   304       y1 = s1 ? s1.y : -1e6;

   305       x1 = e.c - e.b * y1;

   306       y2 = s2 ? s2.y : 1e6;

   307       x2 = e.c - e.b * y2;

   308     } else {

   309       x1 = s1 ? s1.x : -1e6;

   310       y1 = e.c - e.a * x1;

   311       x2 = s2 ? s2.x : 1e6;

   312       y2 = e.c - e.a * x2;

   313     }

   314     var v1 = [x1, y1],

   315         v2 = [x2, y2];

   316     polygons[e.region.l.index].push(v1, v2);

   317     polygons[e.region.r.index].push(v1, v2);

   318   });

   319 

   320   // Reconnect the polygon segments into counterclockwise loops.

   321   return polygons.map(function(polygon, i) {

   322     var cx = vertices[i][0],

   323         cy = vertices[i][1];

   324     polygon.forEach(function(v) {

   325       v.angle = Math.atan2(v[0] - cx, v[1] - cy);

   326     });

   327     return polygon.sort(function(a, b) {

   328       return a.angle - b.angle;

   329     }).filter(function(d, i) {

   330       return !i || (d.angle - polygon[i - 1].angle > 1e-10);

   331     });

   332   });

   333 };

   334 

   335 var d3_voronoi_opposite = {"l": "r", "r": "l"};

   336 

   337 function d3_voronoi_tessellate(vertices, callback) {

   338 

   339   var Sites = {

   340     list: vertices

   341       .map(function(v, i) {

   342         return {

   343           index: i,

   344           x: v[0],

   345           y: v[1]

   346         };

   347       })

   348       .sort(function(a, b) {

   349         return a.y < b.y ? -1

   350           : a.y > b.y ? 1

   351           : a.x < b.x ? -1

   352           : a.x > b.x ? 1

   353           : 0;

   354       }),

   355     bottomSite: null

   356   };

   357 

   358   var EdgeList = {

   359     list: [],

   360     leftEnd: null,

   361     rightEnd: null,

   362 

   363     init: function() {

   364       EdgeList.leftEnd = EdgeList.createHalfEdge(null, "l");

   365       EdgeList.rightEnd = EdgeList.createHalfEdge(null, "l");

   366       EdgeList.leftEnd.r = EdgeList.rightEnd;

   367       EdgeList.rightEnd.l = EdgeList.leftEnd;

   368       EdgeList.list.unshift(EdgeList.leftEnd, EdgeList.rightEnd);

   369     },

   370 

   371     createHalfEdge: function(edge, side) {

   372       return {

   373         edge: edge,

   374         side: side,

   375         vertex: null,

   376         "l": null,

   377         "r": null

   378       };

   379     },

   380 

   381     insert: function(lb, he) {

   382       he.l = lb;

   383       he.r = lb.r;

   384       lb.r.l = he;

   385       lb.r = he;

   386     },

   387 

   388     leftBound: function(p) {

   389       var he = EdgeList.leftEnd;

   390       do {

   391         he = he.r;

   392       } while (he != EdgeList.rightEnd && Geom.rightOf(he, p));

   393       he = he.l;

   394       return he;

   395     },

   396 

   397     del: function(he) {

   398       he.l.r = he.r;

   399       he.r.l = he.l;

   400       he.edge = null;

   401     },

   402 

   403     right: function(he) {

   404       return he.r;

   405     },

   406 

   407     left: function(he) {

   408       return he.l;

   409     },

   410 

   411     leftRegion: function(he) {

   412       return he.edge == null

   413           ? Sites.bottomSite

   414           : he.edge.region[he.side];

   415     },

   416 

   417     rightRegion: function(he) {

   418       return he.edge == null

   419           ? Sites.bottomSite

   420           : he.edge.region[d3_voronoi_opposite[he.side]];

   421     }

   422   };

   423 

   424   var Geom = {

   425 

   426     bisect: function(s1, s2) {

   427       var newEdge = {

   428         region: {"l": s1, "r": s2},

   429         ep: {"l": null, "r": null}

   430       };

   431 

   432       var dx = s2.x - s1.x,

   433           dy = s2.y - s1.y,

   434           adx = dx > 0 ? dx : -dx,

   435           ady = dy > 0 ? dy : -dy;

   436 

   437       newEdge.c = s1.x * dx + s1.y * dy

   438           + (dx * dx + dy * dy) * .5;

   439 

   440       if (adx > ady) {

   441         newEdge.a = 1;

   442         newEdge.b = dy / dx;

   443         newEdge.c /= dx;

   444       } else {

   445         newEdge.b = 1;

   446         newEdge.a = dx / dy;

   447         newEdge.c /= dy;

   448       }

   449 

   450       return newEdge;

   451     },

   452 

   453     intersect: function(el1, el2) {

   454       var e1 = el1.edge,

   455           e2 = el2.edge;

   456       if (!e1 || !e2 || (e1.region.r == e2.region.r)) {

   457         return null;

   458       }

   459       var d = (e1.a * e2.b) - (e1.b * e2.a);

   460       if (Math.abs(d) < 1e-10) {

   461         return null;

   462       }

   463       var xint = (e1.c * e2.b - e2.c * e1.b) / d,

   464           yint = (e2.c * e1.a - e1.c * e2.a) / d,

   465           e1r = e1.region.r,

   466           e2r = e2.region.r,

   467           el,

   468           e;

   469       if ((e1r.y < e2r.y) ||

   470          (e1r.y == e2r.y && e1r.x < e2r.x)) {

   471         el = el1;

   472         e = e1;

   473       } else {

   474         el = el2;

   475         e = e2;

   476       }

   477       var rightOfSite = (xint >= e.region.r.x);

   478       if ((rightOfSite && (el.side === "l")) ||

   479         (!rightOfSite && (el.side === "r"))) {

   480         return null;

   481       }

   482       return {

   483         x: xint,

   484         y: yint

   485       };

   486     },

   487 

   488     rightOf: function(he, p) {

   489       var e = he.edge,

   490           topsite = e.region.r,

   491           rightOfSite = (p.x > topsite.x);

   492 

   493       if (rightOfSite && (he.side === "l")) {

   494         return 1;

   495       }

   496       if (!rightOfSite && (he.side === "r")) {

   497         return 0;

   498       }

   499       if (e.a === 1) {

   500         var dyp = p.y - topsite.y,

   501             dxp = p.x - topsite.x,

   502             fast = 0,

   503             above = 0;

   504 

   505         if ((!rightOfSite && (e.b < 0)) ||

   506           (rightOfSite && (e.b >= 0))) {

   507           above = fast = (dyp >= e.b * dxp);

   508         } else {

   509           above = ((p.x + p.y * e.b) > e.c);

   510           if (e.b < 0) {

   511             above = !above;

   512           }

   513           if (!above) {

   514             fast = 1;

   515           }

   516         }

   517         if (!fast) {

   518           var dxs = topsite.x - e.region.l.x;

   519           above = (e.b * (dxp * dxp - dyp * dyp)) <

   520             (dxs * dyp * (1 + 2 * dxp / dxs + e.b * e.b));

   521 

   522           if (e.b < 0) {

   523             above = !above;

   524           }

   525         }

   526       } else /* e.b == 1 */ {

   527         var yl = e.c - e.a * p.x,

   528             t1 = p.y - yl,

   529             t2 = p.x - topsite.x,

   530             t3 = yl - topsite.y;

   531 

   532         above = (t1 * t1) > (t2 * t2 + t3 * t3);

   533       }

   534       return he.side === "l" ? above : !above;

   535     },

   536 

   537     endPoint: function(edge, side, site) {

   538       edge.ep[side] = site;

   539       if (!edge.ep[d3_voronoi_opposite[side]]) return;

   540       callback(edge);

   541     },

   542 

   543     distance: function(s, t) {

   544       var dx = s.x - t.x,

   545           dy = s.y - t.y;

   546       return Math.sqrt(dx * dx + dy * dy);

   547     }

   548   };

   549 

   550   var EventQueue = {

   551     list: [],

   552 

   553     insert: function(he, site, offset) {

   554       he.vertex = site;

   555       he.ystar = site.y + offset;

   556       for (var i=0, list=EventQueue.list, l=list.length; i<l; i++) {

   557         var next = list[i];

   558         if (he.ystar > next.ystar ||

   559           (he.ystar == next.ystar &&

   560           site.x > next.vertex.x)) {

   561           continue;

   562         } else {

   563           break;

   564         }

   565       }

   566       list.splice(i, 0, he);

   567     },

   568 

   569     del: function(he) {

   570       for (var i=0, ls=EventQueue.list, l=ls.length; i<l && (ls[i] != he); ++i) {}

   571       ls.splice(i, 1);

   572     },

   573 

   574     empty: function() { return EventQueue.list.length === 0; },

   575 

   576     nextEvent: function(he) {

   577       for (var i=0, ls=EventQueue.list, l=ls.length; i<l; ++i) {

   578         if (ls[i] == he) return ls[i+1];

   579       }

   580       return null;

   581     },

   582 

   583     min: function() {

   584       var elem = EventQueue.list[0];

   585       return {

   586         x: elem.vertex.x,

   587         y: elem.ystar

   588       };

   589     },

   590 

   591     extractMin: function() {

   592       return EventQueue.list.shift();

   593     }

   594   };

   595 

   596   EdgeList.init();

   597   Sites.bottomSite = Sites.list.shift();

   598 

   599   var newSite = Sites.list.shift(), newIntStar;

   600   var lbnd, rbnd, llbnd, rrbnd, bisector;

   601   var bot, top, temp, p, v;

   602   var e, pm;

   603 

   604   while (true) {

   605     if (!EventQueue.empty()) {

   606       newIntStar = EventQueue.min();

   607     }

   608     if (newSite && (EventQueue.empty()

   609       || newSite.y < newIntStar.y

   610       || (newSite.y == newIntStar.y

   611       && newSite.x < newIntStar.x))) { //new site is smallest

   612       lbnd = EdgeList.leftBound(newSite);

   613       rbnd = EdgeList.right(lbnd);

   614       bot = EdgeList.rightRegion(lbnd);

   615       e = Geom.bisect(bot, newSite);

   616       bisector = EdgeList.createHalfEdge(e, "l");

   617       EdgeList.insert(lbnd, bisector);

   618       p = Geom.intersect(lbnd, bisector);

   619       if (p) {

   620         EventQueue.del(lbnd);

   621         EventQueue.insert(lbnd, p, Geom.distance(p, newSite));

   622       }

   623       lbnd = bisector;

   624       bisector = EdgeList.createHalfEdge(e, "r");

   625       EdgeList.insert(lbnd, bisector);

   626       p = Geom.intersect(bisector, rbnd);

   627       if (p) {

   628         EventQueue.insert(bisector, p, Geom.distance(p, newSite));

   629       }

   630       newSite = Sites.list.shift();

   631     } else if (!EventQueue.empty()) { //intersection is smallest

   632       lbnd = EventQueue.extractMin();

   633       llbnd = EdgeList.left(lbnd);

   634       rbnd = EdgeList.right(lbnd);

   635       rrbnd = EdgeList.right(rbnd);

   636       bot = EdgeList.leftRegion(lbnd);

   637       top = EdgeList.rightRegion(rbnd);

   638       v = lbnd.vertex;

   639       Geom.endPoint(lbnd.edge, lbnd.side, v);

   640       Geom.endPoint(rbnd.edge, rbnd.side, v);

   641       EdgeList.del(lbnd);

   642       EventQueue.del(rbnd);

   643       EdgeList.del(rbnd);

   644       pm = "l";

   645       if (bot.y > top.y) {

   646         temp = bot;

   647         bot = top;

   648         top = temp;

   649         pm = "r";

   650       }

   651       e = Geom.bisect(bot, top);

   652       bisector = EdgeList.createHalfEdge(e, pm);

   653       EdgeList.insert(llbnd, bisector);

   654       Geom.endPoint(e, d3_voronoi_opposite[pm], v);

   655       p = Geom.intersect(llbnd, bisector);

   656       if (p) {

   657         EventQueue.del(llbnd);

   658         EventQueue.insert(llbnd, p, Geom.distance(p, bot));

   659       }

   660       p = Geom.intersect(bisector, rrbnd);

   661       if (p) {

   662         EventQueue.insert(bisector, p, Geom.distance(p, bot));

   663       }

   664     } else {

   665       break;

   666     }

   667   }//end while

   668 

   669   for (lbnd = EdgeList.right(EdgeList.leftEnd);

   670       lbnd != EdgeList.rightEnd;

   671       lbnd = EdgeList.right(lbnd)) {

   672     callback(lbnd.edge);

   673   }

   674 }

   675 /**

   676 * @param vertices [[x1, y1], [x2, y2], …]

   677 * @returns triangles [[[x1, y1], [x2, y2], [x3, y3]], …]

   678  */

   679 d3.geom.delaunay = function(vertices) {

   680   var edges = vertices.map(function() { return []; }),

   681       triangles = [];

   682 

   683   // Use the Voronoi tessellation to determine Delaunay edges.

   684   d3_voronoi_tessellate(vertices, function(e) {

   685     edges[e.region.l.index].push(vertices[e.region.r.index]);

   686   });

   687 

   688   // Reconnect the edges into counterclockwise triangles.

   689   edges.forEach(function(edge, i) {

   690     var v = vertices[i],

   691         cx = v[0],

   692         cy = v[1];

   693     edge.forEach(function(v) {

   694       v.angle = Math.atan2(v[0] - cx, v[1] - cy);

   695     });

   696     edge.sort(function(a, b) {

   697       return a.angle - b.angle;

   698     });

   699     for (var j = 0, m = edge.length - 1; j < m; j++) {

   700       triangles.push([v, edge[j], edge[j + 1]]);

   701     }

   702   });

   703 

   704   return triangles;

   705 };

   706 // Constructs a new quadtree for the specified array of points. A quadtree is a

   707 // two-dimensional recursive spatial subdivision. This implementation uses

   708 // square partitions, dividing each square into four equally-sized squares. Each

   709 // point exists in a unique node; if multiple points are in the same position,

   710 // some points may be stored on internal nodes rather than leaf nodes. Quadtrees

   711 // can be used to accelerate various spatial operations, such as the Barnes-Hut

   712 // approximation for computing n-body forces, or collision detection.

   713 d3.geom.quadtree = function(points, x1, y1, x2, y2) {

   714   var p,

   715       i = -1,

   716       n = points.length;

   717 

   718   // Type conversion for deprecated API.

   719   if (n && isNaN(points[0].x)) points = points.map(d3_geom_quadtreePoint);

   720 

   721   // Allow bounds to be specified explicitly.

   722   if (arguments.length < 5) {

   723     if (arguments.length === 3) {

   724       y2 = x2 = y1;

   725       y1 = x1;

   726     } else {

   727       x1 = y1 = Infinity;

   728       x2 = y2 = -Infinity;

   729 

   730       // Compute bounds.

   731       while (++i < n) {

   732         p = points[i];

   733         if (p.x < x1) x1 = p.x;

   734         if (p.y < y1) y1 = p.y;

   735         if (p.x > x2) x2 = p.x;

   736         if (p.y > y2) y2 = p.y;

   737       }

   738 

   739       // Squarify the bounds.

   740       var dx = x2 - x1,

   741           dy = y2 - y1;

   742       if (dx > dy) y2 = y1 + dx;

   743       else x2 = x1 + dy;

   744     }

   745   }

   746 

   747   // Recursively inserts the specified point p at the node n or one of its

   748   // descendants. The bounds are defined by [x1, x2] and [y1, y2].

   749   function insert(n, p, x1, y1, x2, y2) {

   750     if (isNaN(p.x) || isNaN(p.y)) return; // ignore invalid points

   751     if (n.leaf) {

   752       var v = n.point;

   753       if (v) {

   754         // If the point at this leaf node is at the same position as the new

   755         // point we are adding, we leave the point associated with the

   756         // internal node while adding the new point to a child node. This

   757         // avoids infinite recursion.

   758         if ((Math.abs(v.x - p.x) + Math.abs(v.y - p.y)) < .01) {

   759           insertChild(n, p, x1, y1, x2, y2);

   760         } else {

   761           n.point = null;

   762           insertChild(n, v, x1, y1, x2, y2);

   763           insertChild(n, p, x1, y1, x2, y2);

   764         }

   765       } else {

   766         n.point = p;

   767       }

   768     } else {

   769       insertChild(n, p, x1, y1, x2, y2);

   770     }

   771   }

   772 

   773   // Recursively inserts the specified point p into a descendant of node n. The

   774   // bounds are defined by [x1, x2] and [y1, y2].

   775   function insertChild(n, p, x1, y1, x2, y2) {

   776     // Compute the split point, and the quadrant in which to insert p.

   777     var sx = (x1 + x2) * .5,

   778         sy = (y1 + y2) * .5,

   779         right = p.x >= sx,

   780         bottom = p.y >= sy,

   781         i = (bottom << 1) + right;

   782 

   783     // Recursively insert into the child node.

   784     n.leaf = false;

   785     n = n.nodes[i] || (n.nodes[i] = d3_geom_quadtreeNode());

   786 

   787     // Update the bounds as we recurse.

   788     if (right) x1 = sx; else x2 = sx;

   789     if (bottom) y1 = sy; else y2 = sy;

   790     insert(n, p, x1, y1, x2, y2);

   791   }

   792 

   793   // Create the root node.

   794   var root = d3_geom_quadtreeNode();

   795 

   796   root.add = function(p) {

   797     insert(root, p, x1, y1, x2, y2);

   798   };

   799 

   800   root.visit = function(f) {

   801     d3_geom_quadtreeVisit(f, root, x1, y1, x2, y2);

   802   };

   803 

   804   // Insert all points.

   805   points.forEach(root.add);

   806   return root;

   807 };

   808 

   809 function d3_geom_quadtreeNode() {

   810   return {

   811     leaf: true,

   812     nodes: [],

   813     point: null

   814   };

   815 }

   816 

   817 function d3_geom_quadtreeVisit(f, node, x1, y1, x2, y2) {

   818   if (!f(node, x1, y1, x2, y2)) {

   819     var sx = (x1 + x2) * .5,

   820         sy = (y1 + y2) * .5,

   821         children = node.nodes;

   822     if (children[0]) d3_geom_quadtreeVisit(f, children[0], x1, y1, sx, sy);

   823     if (children[1]) d3_geom_quadtreeVisit(f, children[1], sx, y1, x2, sy);

   824     if (children[2]) d3_geom_quadtreeVisit(f, children[2], x1, sy, sx, y2);

   825     if (children[3]) d3_geom_quadtreeVisit(f, children[3], sx, sy, x2, y2);

   826   }

   827 }

   828 

   829 function d3_geom_quadtreePoint(p) {

   830   return {

   831     x: p[0],

   832     y: p[1]

   833   };

   834 }

   835 })();