source: Dev/branches/Cartis/d3/d3.geom.js @ 297

Last change on this file since 297 was 76, checked in by fpvanagthoven, 14 years ago

d3

File size: 21.2 KB
Line 
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 */
12d3.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
58var 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
61function 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 */
89d3.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?
171function 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!
179d3.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
249function 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.
254function 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/**
272 * @param vertices [[x1, y1], [x2, y2], 
]
273 * @returns polygons [[[x1, y1], [x2, y2], 
], 
]
274 */
275d3.geom.voronoi = function(vertices) {
276  var polygons = vertices.map(function() { return []; });
277
278  // Note: we expect the caller to clip the polygons, if needed.
279  d3_voronoi_tessellate(vertices, function(e) {
280    var s1,
281        s2,
282        x1,
283        x2,
284        y1,
285        y2;
286    if (e.a === 1 && e.b >= 0) {
287      s1 = e.ep.r;
288      s2 = e.ep.l;
289    } else {
290      s1 = e.ep.l;
291      s2 = e.ep.r;
292    }
293    if (e.a === 1) {
294      y1 = s1 ? s1.y : -1e6;
295      x1 = e.c - e.b * y1;
296      y2 = s2 ? s2.y : 1e6;
297      x2 = e.c - e.b * y2;
298    } else {
299      x1 = s1 ? s1.x : -1e6;
300      y1 = e.c - e.a * x1;
301      x2 = s2 ? s2.x : 1e6;
302      y2 = e.c - e.a * x2;
303    }
304    var v1 = [x1, y1],
305        v2 = [x2, y2];
306    polygons[e.region.l.index].push(v1, v2);
307    polygons[e.region.r.index].push(v1, v2);
308  });
309
310  // Reconnect the polygon segments into counterclockwise loops.
311  return polygons.map(function(polygon, i) {
312    var cx = vertices[i][0],
313        cy = vertices[i][1];
314    polygon.forEach(function(v) {
315      v.angle = Math.atan2(v[0] - cx, v[1] - cy);
316    });
317    return polygon.sort(function(a, b) {
318      return a.angle - b.angle;
319    }).filter(function(d, i) {
320      return !i || (d.angle - polygon[i - 1].angle > 1e-10);
321    });
322  });
323};
324
325var d3_voronoi_opposite = {"l": "r", "r": "l"};
326
327function d3_voronoi_tessellate(vertices, callback) {
328
329  var Sites = {
330    list: vertices
331      .map(function(v, i) {
332        return {
333          index: i,
334          x: v[0],
335          y: v[1]
336        };
337      })
338      .sort(function(a, b) {
339        return a.y < b.y ? -1
340          : a.y > b.y ? 1
341          : a.x < b.x ? -1
342          : a.x > b.x ? 1
343          : 0;
344      }),
345    bottomSite: null
346  };
347
348  var EdgeList = {
349    list: [],
350    leftEnd: null,
351    rightEnd: null,
352
353    init: function() {
354      EdgeList.leftEnd = EdgeList.createHalfEdge(null, "l");
355      EdgeList.rightEnd = EdgeList.createHalfEdge(null, "l");
356      EdgeList.leftEnd.r = EdgeList.rightEnd;
357      EdgeList.rightEnd.l = EdgeList.leftEnd;
358      EdgeList.list.unshift(EdgeList.leftEnd, EdgeList.rightEnd);
359    },
360
361    createHalfEdge: function(edge, side) {
362      return {
363        edge: edge,
364        side: side,
365        vertex: null,
366        "l": null,
367        "r": null
368      };
369    },
370
371    insert: function(lb, he) {
372      he.l = lb;
373      he.r = lb.r;
374      lb.r.l = he;
375      lb.r = he;
376    },
377
378    leftBound: function(p) {
379      var he = EdgeList.leftEnd;
380      do {
381        he = he.r;
382      } while (he != EdgeList.rightEnd && Geom.rightOf(he, p));
383      he = he.l;
384      return he;
385    },
386
387    del: function(he) {
388      he.l.r = he.r;
389      he.r.l = he.l;
390      he.edge = null;
391    },
392
393    right: function(he) {
394      return he.r;
395    },
396
397    left: function(he) {
398      return he.l;
399    },
400
401    leftRegion: function(he) {
402      return he.edge == null
403          ? Sites.bottomSite
404          : he.edge.region[he.side];
405    },
406
407    rightRegion: function(he) {
408      return he.edge == null
409          ? Sites.bottomSite
410          : he.edge.region[d3_voronoi_opposite[he.side]];
411    }
412  };
413
414  var Geom = {
415
416    bisect: function(s1, s2) {
417      var newEdge = {
418        region: {"l": s1, "r": s2},
419        ep: {"l": null, "r": null}
420      };
421
422      var dx = s2.x - s1.x,
423          dy = s2.y - s1.y,
424          adx = dx > 0 ? dx : -dx,
425          ady = dy > 0 ? dy : -dy;
426
427      newEdge.c = s1.x * dx + s1.y * dy
428          + (dx * dx + dy * dy) * .5;
429
430      if (adx > ady) {
431        newEdge.a = 1;
432        newEdge.b = dy / dx;
433        newEdge.c /= dx;
434      } else {
435        newEdge.b = 1;
436        newEdge.a = dx / dy;
437        newEdge.c /= dy;
438      }
439
440      return newEdge;
441    },
442
443    intersect: function(el1, el2) {
444      var e1 = el1.edge,
445          e2 = el2.edge;
446      if (!e1 || !e2 || (e1.region.r == e2.region.r)) {
447        return null;
448      }
449      var d = (e1.a * e2.b) - (e1.b * e2.a);
450      if (Math.abs(d) < 1e-10) {
451        return null;
452      }
453      var xint = (e1.c * e2.b - e2.c * e1.b) / d,
454          yint = (e2.c * e1.a - e1.c * e2.a) / d,
455          e1r = e1.region.r,
456          e2r = e2.region.r,
457          el,
458          e;
459      if ((e1r.y < e2r.y) ||
460         (e1r.y == e2r.y && e1r.x < e2r.x)) {
461        el = el1;
462        e = e1;
463      } else {
464        el = el2;
465        e = e2;
466      }
467      var rightOfSite = (xint >= e.region.r.x);
468      if ((rightOfSite && (el.side === "l")) ||
469        (!rightOfSite && (el.side === "r"))) {
470        return null;
471      }
472      return {
473        x: xint,
474        y: yint
475      };
476    },
477
478    rightOf: function(he, p) {
479      var e = he.edge,
480          topsite = e.region.r,
481          rightOfSite = (p.x > topsite.x);
482
483      if (rightOfSite && (he.side === "l")) {
484        return 1;
485      }
486      if (!rightOfSite && (he.side === "r")) {
487        return 0;
488      }
489      if (e.a === 1) {
490        var dyp = p.y - topsite.y,
491            dxp = p.x - topsite.x,
492            fast = 0,
493            above = 0;
494
495        if ((!rightOfSite && (e.b < 0)) ||
496          (rightOfSite && (e.b >= 0))) {
497          above = fast = (dyp >= e.b * dxp);
498        } else {
499          above = ((p.x + p.y * e.b) > e.c);
500          if (e.b < 0) {
501            above = !above;
502          }
503          if (!above) {
504            fast = 1;
505          }
506        }
507        if (!fast) {
508          var dxs = topsite.x - e.region.l.x;
509          above = (e.b * (dxp * dxp - dyp * dyp)) <
510            (dxs * dyp * (1 + 2 * dxp / dxs + e.b * e.b));
511
512          if (e.b < 0) {
513            above = !above;
514          }
515        }
516      } else /* e.b == 1 */ {
517        var yl = e.c - e.a * p.x,
518            t1 = p.y - yl,
519            t2 = p.x - topsite.x,
520            t3 = yl - topsite.y;
521
522        above = (t1 * t1) > (t2 * t2 + t3 * t3);
523      }
524      return he.side === "l" ? above : !above;
525    },
526
527    endPoint: function(edge, side, site) {
528      edge.ep[side] = site;
529      if (!edge.ep[d3_voronoi_opposite[side]]) return;
530      callback(edge);
531    },
532
533    distance: function(s, t) {
534      var dx = s.x - t.x,
535          dy = s.y - t.y;
536      return Math.sqrt(dx * dx + dy * dy);
537    }
538  };
539
540  var EventQueue = {
541    list: [],
542
543    insert: function(he, site, offset) {
544      he.vertex = site;
545      he.ystar = site.y + offset;
546      for (var i=0, list=EventQueue.list, l=list.length; i<l; i++) {
547        var next = list[i];
548        if (he.ystar > next.ystar ||
549          (he.ystar == next.ystar &&
550          site.x > next.vertex.x)) {
551          continue;
552        } else {
553          break;
554        }
555      }
556      list.splice(i, 0, he);
557    },
558
559    del: function(he) {
560      for (var i=0, ls=EventQueue.list, l=ls.length; i<l && (ls[i] != he); ++i) {}
561      ls.splice(i, 1);
562    },
563
564    empty: function() { return EventQueue.list.length === 0; },
565
566    nextEvent: function(he) {
567      for (var i=0, ls=EventQueue.list, l=ls.length; i<l; ++i) {
568        if (ls[i] == he) return ls[i+1];
569      }
570      return null;
571    },
572
573    min: function() {
574      var elem = EventQueue.list[0];
575      return {
576        x: elem.vertex.x,
577        y: elem.ystar
578      };
579    },
580
581    extractMin: function() {
582      return EventQueue.list.shift();
583    }
584  };
585
586  EdgeList.init();
587  Sites.bottomSite = Sites.list.shift();
588
589  var newSite = Sites.list.shift(), newIntStar;
590  var lbnd, rbnd, llbnd, rrbnd, bisector;
591  var bot, top, temp, p, v;
592  var e, pm;
593
594  while (true) {
595    if (!EventQueue.empty()) {
596      newIntStar = EventQueue.min();
597    }
598    if (newSite && (EventQueue.empty()
599      || newSite.y < newIntStar.y
600      || (newSite.y == newIntStar.y
601      && newSite.x < newIntStar.x))) { //new site is smallest
602      lbnd = EdgeList.leftBound(newSite);
603      rbnd = EdgeList.right(lbnd);
604      bot = EdgeList.rightRegion(lbnd);
605      e = Geom.bisect(bot, newSite);
606      bisector = EdgeList.createHalfEdge(e, "l");
607      EdgeList.insert(lbnd, bisector);
608      p = Geom.intersect(lbnd, bisector);
609      if (p) {
610        EventQueue.del(lbnd);
611        EventQueue.insert(lbnd, p, Geom.distance(p, newSite));
612      }
613      lbnd = bisector;
614      bisector = EdgeList.createHalfEdge(e, "r");
615      EdgeList.insert(lbnd, bisector);
616      p = Geom.intersect(bisector, rbnd);
617      if (p) {
618        EventQueue.insert(bisector, p, Geom.distance(p, newSite));
619      }
620      newSite = Sites.list.shift();
621    } else if (!EventQueue.empty()) { //intersection is smallest
622      lbnd = EventQueue.extractMin();
623      llbnd = EdgeList.left(lbnd);
624      rbnd = EdgeList.right(lbnd);
625      rrbnd = EdgeList.right(rbnd);
626      bot = EdgeList.leftRegion(lbnd);
627      top = EdgeList.rightRegion(rbnd);
628      v = lbnd.vertex;
629      Geom.endPoint(lbnd.edge, lbnd.side, v);
630      Geom.endPoint(rbnd.edge, rbnd.side, v);
631      EdgeList.del(lbnd);
632      EventQueue.del(rbnd);
633      EdgeList.del(rbnd);
634      pm = "l";
635      if (bot.y > top.y) {
636        temp = bot;
637        bot = top;
638        top = temp;
639        pm = "r";
640      }
641      e = Geom.bisect(bot, top);
642      bisector = EdgeList.createHalfEdge(e, pm);
643      EdgeList.insert(llbnd, bisector);
644      Geom.endPoint(e, d3_voronoi_opposite[pm], v);
645      p = Geom.intersect(llbnd, bisector);
646      if (p) {
647        EventQueue.del(llbnd);
648        EventQueue.insert(llbnd, p, Geom.distance(p, bot));
649      }
650      p = Geom.intersect(bisector, rrbnd);
651      if (p) {
652        EventQueue.insert(bisector, p, Geom.distance(p, bot));
653      }
654    } else {
655      break;
656    }
657  }//end while
658
659  for (lbnd = EdgeList.right(EdgeList.leftEnd);
660      lbnd != EdgeList.rightEnd;
661      lbnd = EdgeList.right(lbnd)) {
662    callback(lbnd.edge);
663  }
664}
665/**
666* @param vertices [[x1, y1], [x2, y2], 
]
667* @returns triangles [[[x1, y1], [x2, y2], [x3, y3]], 
]
668 */
669d3.geom.delaunay = function(vertices) {
670  var edges = vertices.map(function() { return []; }),
671      triangles = [];
672
673  // Use the Voronoi tessellation to determine Delaunay edges.
674  d3_voronoi_tessellate(vertices, function(e) {
675    edges[e.region.l.index].push(vertices[e.region.r.index]);
676  });
677
678  // Reconnect the edges into counterclockwise triangles.
679  edges.forEach(function(edge, i) {
680    var v = vertices[i],
681        cx = v[0],
682        cy = v[1];
683    edge.forEach(function(v) {
684      v.angle = Math.atan2(v[0] - cx, v[1] - cy);
685    });
686    edge.sort(function(a, b) {
687      return a.angle - b.angle;
688    });
689    for (var j = 0, m = edge.length - 1; j < m; j++) {
690      triangles.push([v, edge[j], edge[j + 1]]);
691    }
692  });
693
694  return triangles;
695};
696// Constructs a new quadtree for the specified array of points. A quadtree is a
697// two-dimensional recursive spatial subdivision. This implementation uses
698// square partitions, dividing each square into four equally-sized squares. Each
699// point exists in a unique node; if multiple points are in the same position,
700// some points may be stored on internal nodes rather than leaf nodes. Quadtrees
701// can be used to accelerate various spatial operations, such as the Barnes-Hut
702// approximation for computing n-body forces, or collision detection.
703d3.geom.quadtree = function(points, x1, y1, x2, y2) {
704  var p,
705      i = -1,
706      n = points.length;
707
708  // Type conversion for deprecated API.
709  if (n && isNaN(points[0].x)) points = points.map(d3_geom_quadtreePoint);
710
711  // Allow bounds to be specified explicitly.
712  if (arguments.length < 5) {
713    if (arguments.length === 3) {
714      y2 = x2 = y1;
715      y1 = x1;
716    } else {
717      x1 = y1 = Infinity;
718      x2 = y2 = -Infinity;
719
720      // Compute bounds.
721      while (++i < n) {
722        p = points[i];
723        if (p.x < x1) x1 = p.x;
724        if (p.y < y1) y1 = p.y;
725        if (p.x > x2) x2 = p.x;
726        if (p.y > y2) y2 = p.y;
727      }
728
729      // Squarify the bounds.
730      var dx = x2 - x1,
731          dy = y2 - y1;
732      if (dx > dy) y2 = y1 + dx;
733      else x2 = x1 + dy;
734    }
735  }
736
737  // Recursively inserts the specified point p at the node n or one of its
738  // descendants. The bounds are defined by [x1, x2] and [y1, y2].
739  function insert(n, p, x1, y1, x2, y2) {
740    if (isNaN(p.x) || isNaN(p.y)) return; // ignore invalid points
741    if (n.leaf) {
742      var v = n.point;
743      if (v) {
744        // If the point at this leaf node is at the same position as the new
745        // point we are adding, we leave the point associated with the
746        // internal node while adding the new point to a child node. This
747        // avoids infinite recursion.
748        if ((Math.abs(v.x - p.x) + Math.abs(v.y - p.y)) < .01) {
749          insertChild(n, p, x1, y1, x2, y2);
750        } else {
751          n.point = null;
752          insertChild(n, v, x1, y1, x2, y2);
753          insertChild(n, p, x1, y1, x2, y2);
754        }
755      } else {
756        n.point = p;
757      }
758    } else {
759      insertChild(n, p, x1, y1, x2, y2);
760    }
761  }
762
763  // Recursively inserts the specified point p into a descendant of node n. The
764  // bounds are defined by [x1, x2] and [y1, y2].
765  function insertChild(n, p, x1, y1, x2, y2) {
766    // Compute the split point, and the quadrant in which to insert p.
767    var sx = (x1 + x2) * .5,
768        sy = (y1 + y2) * .5,
769        right = p.x >= sx,
770        bottom = p.y >= sy,
771        i = (bottom << 1) + right;
772
773    // Recursively insert into the child node.
774    n.leaf = false;
775    n = n.nodes[i] || (n.nodes[i] = d3_geom_quadtreeNode());
776
777    // Update the bounds as we recurse.
778    if (right) x1 = sx; else x2 = sx;
779    if (bottom) y1 = sy; else y2 = sy;
780    insert(n, p, x1, y1, x2, y2);
781  }
782
783  // Create the root node.
784  var root = d3_geom_quadtreeNode();
785
786  root.add = function(p) {
787    insert(root, p, x1, y1, x2, y2);
788  };
789
790  root.visit = function(f) {
791    d3_geom_quadtreeVisit(f, root, x1, y1, x2, y2);
792  };
793
794  // Insert all points.
795  points.forEach(root.add);
796  return root;
797};
798
799function d3_geom_quadtreeNode() {
800  return {
801    leaf: true,
802    nodes: [],
803    point: null
804  };
805}
806
807function d3_geom_quadtreeVisit(f, node, x1, y1, x2, y2) {
808  if (!f(node, x1, y1, x2, y2)) {
809    var sx = (x1 + x2) * .5,
810        sy = (y1 + y2) * .5,
811        children = node.nodes;
812    if (children[0]) d3_geom_quadtreeVisit(f, children[0], x1, y1, sx, sy);
813    if (children[1]) d3_geom_quadtreeVisit(f, children[1], sx, y1, x2, sy);
814    if (children[2]) d3_geom_quadtreeVisit(f, children[2], x1, sy, sx, y2);
815    if (children[3]) d3_geom_quadtreeVisit(f, children[3], sx, sy, x2, y2);
816  }
817}
818
819function d3_geom_quadtreePoint(p) {
820  return {
821    x: p[0],
822    y: p[1]
823  };
824}
825})();
Note: See TracBrowser for help on using the repository browser.