d3.geom.js 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836
  1. // See d3.LICENSE.txt for license.
  2. (function(){d3.geom = {};
  3. /**
  4. * Computes a contour for a given input grid function using the <a
  5. * href="https://en.wikipedia.org/wiki/Marching_squares">marching
  6. * squares</a> algorithm. Returns the contour polygon as an array of points.
  7. *
  8. * @param grid a two-input function(x, y) that returns true for values
  9. * inside the contour and false for values outside the contour.
  10. * @param start an optional starting point [x, y] on the grid.
  11. * @returns polygon [[x1, y1], [x2, y2], …]
  12. */
  13. d3.geom.contour = function(grid, start) {
  14. var s = start || d3_geom_contourStart(grid), // starting point
  15. c = [], // contour polygon
  16. x = s[0], // current x position
  17. y = s[1], // current y position
  18. dx = 0, // next x direction
  19. dy = 0, // next y direction
  20. pdx = NaN, // previous x direction
  21. pdy = NaN, // previous y direction
  22. i = 0;
  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. // determine next direction
  31. if (i === 6) {
  32. dx = pdy === -1 ? -1 : 1;
  33. dy = 0;
  34. } else if (i === 9) {
  35. dx = 0;
  36. dy = pdx === 1 ? -1 : 1;
  37. } else {
  38. dx = d3_geom_contourDx[i];
  39. dy = d3_geom_contourDy[i];
  40. }
  41. // update contour polygon
  42. if (dx != pdx && dy != pdy) {
  43. c.push([x, y]);
  44. pdx = dx;
  45. pdy = dy;
  46. }
  47. x += dx;
  48. y += dy;
  49. } while (s[0] != x || s[1] != y);
  50. return c;
  51. };
  52. // lookup tables for marching directions
  53. var d3_geom_contourDx = [1, 0, 1, 1,-1, 0,-1, 1,0, 0,0,0,-1, 0,-1,NaN],
  54. d3_geom_contourDy = [0,-1, 0, 0, 0,-1, 0, 0,1,-1,1,1, 0,-1, 0,NaN];
  55. function d3_geom_contourStart(grid) {
  56. var x = 0,
  57. y = 0;
  58. // search for a starting point; begin at origin
  59. // and proceed along outward-expanding diagonals
  60. while (true) {
  61. if (grid(x,y)) {
  62. return [x,y];
  63. }
  64. if (x === 0) {
  65. x = y + 1;
  66. y = 0;
  67. } else {
  68. x = x - 1;
  69. y = y + 1;
  70. }
  71. }
  72. }
  73. /**
  74. * Computes the 2D convex hull of a set of points using Graham's scanning
  75. * algorithm. The algorithm has been implemented as described in Cormen,
  76. * Leiserson, and Rivest's Introduction to Algorithms. The running time of
  77. * this algorithm is O(n log n), where n is the number of input points.
  78. *
  79. * @param vertices [[x1, y1], [x2, y2], …]
  80. * @returns polygon [[x1, y1], [x2, y2], …]
  81. */
  82. d3.geom.hull = function(vertices) {
  83. if (vertices.length < 3) return [];
  84. var len = vertices.length,
  85. plen = len - 1,
  86. points = [],
  87. stack = [],
  88. i, j, h = 0, x1, y1, x2, y2, u, v, a, sp;
  89. // find the starting ref point: leftmost point with the minimum y coord
  90. for (i=1; i<len; ++i) {
  91. if (vertices[i][1] < vertices[h][1]) {
  92. h = i;
  93. } else if (vertices[i][1] == vertices[h][1]) {
  94. h = (vertices[i][0] < vertices[h][0] ? i : h);
  95. }
  96. }
  97. // calculate polar angles from ref point and sort
  98. for (i=0; i<len; ++i) {
  99. if (i === h) continue;
  100. y1 = vertices[i][1] - vertices[h][1];
  101. x1 = vertices[i][0] - vertices[h][0];
  102. points.push({angle: Math.atan2(y1, x1), index: i});
  103. }
  104. points.sort(function(a, b) { return a.angle - b.angle; });
  105. // toss out duplicate angles
  106. a = points[0].angle;
  107. v = points[0].index;
  108. u = 0;
  109. for (i=1; i<plen; ++i) {
  110. j = points[i].index;
  111. if (a == points[i].angle) {
  112. // keep angle for point most distant from the reference
  113. x1 = vertices[v][0] - vertices[h][0];
  114. y1 = vertices[v][1] - vertices[h][1];
  115. x2 = vertices[j][0] - vertices[h][0];
  116. y2 = vertices[j][1] - vertices[h][1];
  117. if ((x1*x1 + y1*y1) >= (x2*x2 + y2*y2)) {
  118. points[i].index = -1;
  119. } else {
  120. points[u].index = -1;
  121. a = points[i].angle;
  122. u = i;
  123. v = j;
  124. }
  125. } else {
  126. a = points[i].angle;
  127. u = i;
  128. v = j;
  129. }
  130. }
  131. // initialize the stack
  132. stack.push(h);
  133. for (i=0, j=0; i<2; ++j) {
  134. if (points[j].index !== -1) {
  135. stack.push(points[j].index);
  136. i++;
  137. }
  138. }
  139. sp = stack.length;
  140. // do graham's scan
  141. for (; j<plen; ++j) {
  142. if (points[j].index === -1) continue; // skip tossed out points
  143. while (!d3_geom_hullCCW(stack[sp-2], stack[sp-1], points[j].index, vertices)) {
  144. --sp;
  145. }
  146. stack[sp++] = points[j].index;
  147. }
  148. // construct the hull
  149. var poly = [];
  150. for (i=0; i<sp; ++i) {
  151. poly.push(vertices[stack[i]]);
  152. }
  153. return poly;
  154. }
  155. // are three points in counter-clockwise order?
  156. function d3_geom_hullCCW(i1, i2, i3, v) {
  157. var t, a, b, c, d, e, f;
  158. t = v[i1]; a = t[0]; b = t[1];
  159. t = v[i2]; c = t[0]; d = t[1];
  160. t = v[i3]; e = t[0]; f = t[1];
  161. return ((f-b)*(c-a) - (d-b)*(e-a)) > 0;
  162. }
  163. // Note: requires coordinates to be counterclockwise and convex!
  164. d3.geom.polygon = function(coordinates) {
  165. coordinates.area = function() {
  166. var i = 0,
  167. n = coordinates.length,
  168. a = coordinates[n - 1][0] * coordinates[0][1],
  169. b = coordinates[n - 1][1] * coordinates[0][0];
  170. while (++i < n) {
  171. a += coordinates[i - 1][0] * coordinates[i][1];
  172. b += coordinates[i - 1][1] * coordinates[i][0];
  173. }
  174. return (b - a) * .5;
  175. };
  176. coordinates.centroid = function(k) {
  177. var i = -1,
  178. n = coordinates.length,
  179. x = 0,
  180. y = 0,
  181. a,
  182. b = coordinates[n - 1],
  183. c;
  184. if (!arguments.length) k = -1 / (6 * coordinates.area());
  185. while (++i < n) {
  186. a = b;
  187. b = coordinates[i];
  188. c = a[0] * b[1] - b[0] * a[1];
  189. x += (a[0] + b[0]) * c;
  190. y += (a[1] + b[1]) * c;
  191. }
  192. return [x * k, y * k];
  193. };
  194. // The Sutherland-Hodgman clipping algorithm.
  195. coordinates.clip = function(subject) {
  196. var input,
  197. i = -1,
  198. n = coordinates.length,
  199. j,
  200. m,
  201. a = coordinates[n - 1],
  202. b,
  203. c,
  204. d;
  205. while (++i < n) {
  206. input = subject.slice();
  207. subject.length = 0;
  208. b = coordinates[i];
  209. c = input[(m = input.length) - 1];
  210. j = -1;
  211. while (++j < m) {
  212. d = input[j];
  213. if (d3_geom_polygonInside(d, a, b)) {
  214. if (!d3_geom_polygonInside(c, a, b)) {
  215. subject.push(d3_geom_polygonIntersect(c, d, a, b));
  216. }
  217. subject.push(d);
  218. } else if (d3_geom_polygonInside(c, a, b)) {
  219. subject.push(d3_geom_polygonIntersect(c, d, a, b));
  220. }
  221. c = d;
  222. }
  223. a = b;
  224. }
  225. return subject;
  226. };
  227. return coordinates;
  228. };
  229. function d3_geom_polygonInside(p, a, b) {
  230. return (b[0] - a[0]) * (p[1] - a[1]) < (b[1] - a[1]) * (p[0] - a[0]);
  231. }
  232. // Intersect two infinite lines cd and ab.
  233. function d3_geom_polygonIntersect(c, d, a, b) {
  234. var x1 = c[0], x2 = d[0], x3 = a[0], x4 = b[0],
  235. y1 = c[1], y2 = d[1], y3 = a[1], y4 = b[1],
  236. x13 = x1 - x3,
  237. x21 = x2 - x1,
  238. x43 = x4 - x3,
  239. y13 = y1 - y3,
  240. y21 = y2 - y1,
  241. y43 = y4 - y3,
  242. ua = (x43 * y13 - y43 * x13) / (y43 * x21 - x43 * y21);
  243. return [x1 + ua * x21, y1 + ua * y21];
  244. }
  245. // Adapted from Nicolas Garcia Belmonte's JIT implementation:
  246. // http://blog.thejit.org/2010/02/12/voronoi-tessellation/
  247. // http://blog.thejit.org/assets/voronoijs/voronoi.js
  248. // See lib/jit/LICENSE for details.
  249. // Notes:
  250. //
  251. // This implementation does not clip the returned polygons, so if you want to
  252. // clip them to a particular shape you will need to do that either in SVG or by
  253. // post-processing with d3.geom.polygon's clip method.
  254. //
  255. // If any vertices are coincident or have NaN positions, the behavior of this
  256. // method is undefined. Most likely invalid polygons will be returned. You
  257. // should filter invalid points, and consolidate coincident points, before
  258. // computing the tessellation.
  259. /**
  260. * @param vertices [[x1, y1], [x2, y2], …]
  261. * @returns polygons [[[x1, y1], [x2, y2], …], …]
  262. */
  263. d3.geom.voronoi = function(vertices) {
  264. var polygons = vertices.map(function() { return []; });
  265. d3_voronoi_tessellate(vertices, function(e) {
  266. var s1,
  267. s2,
  268. x1,
  269. x2,
  270. y1,
  271. y2;
  272. if (e.a === 1 && e.b >= 0) {
  273. s1 = e.ep.r;
  274. s2 = e.ep.l;
  275. } else {
  276. s1 = e.ep.l;
  277. s2 = e.ep.r;
  278. }
  279. if (e.a === 1) {
  280. y1 = s1 ? s1.y : -1e6;
  281. x1 = e.c - e.b * y1;
  282. y2 = s2 ? s2.y : 1e6;
  283. x2 = e.c - e.b * y2;
  284. } else {
  285. x1 = s1 ? s1.x : -1e6;
  286. y1 = e.c - e.a * x1;
  287. x2 = s2 ? s2.x : 1e6;
  288. y2 = e.c - e.a * x2;
  289. }
  290. var v1 = [x1, y1],
  291. v2 = [x2, y2];
  292. polygons[e.region.l.index].push(v1, v2);
  293. polygons[e.region.r.index].push(v1, v2);
  294. });
  295. // Reconnect the polygon segments into counterclockwise loops.
  296. return polygons.map(function(polygon, i) {
  297. var cx = vertices[i][0],
  298. cy = vertices[i][1];
  299. polygon.forEach(function(v) {
  300. v.angle = Math.atan2(v[0] - cx, v[1] - cy);
  301. });
  302. return polygon.sort(function(a, b) {
  303. return a.angle - b.angle;
  304. }).filter(function(d, i) {
  305. return !i || (d.angle - polygon[i - 1].angle > 1e-10);
  306. });
  307. });
  308. };
  309. var d3_voronoi_opposite = {"l": "r", "r": "l"};
  310. function d3_voronoi_tessellate(vertices, callback) {
  311. var Sites = {
  312. list: vertices
  313. .map(function(v, i) {
  314. return {
  315. index: i,
  316. x: v[0],
  317. y: v[1]
  318. };
  319. })
  320. .sort(function(a, b) {
  321. return a.y < b.y ? -1
  322. : a.y > b.y ? 1
  323. : a.x < b.x ? -1
  324. : a.x > b.x ? 1
  325. : 0;
  326. }),
  327. bottomSite: null
  328. };
  329. var EdgeList = {
  330. list: [],
  331. leftEnd: null,
  332. rightEnd: null,
  333. init: function() {
  334. EdgeList.leftEnd = EdgeList.createHalfEdge(null, "l");
  335. EdgeList.rightEnd = EdgeList.createHalfEdge(null, "l");
  336. EdgeList.leftEnd.r = EdgeList.rightEnd;
  337. EdgeList.rightEnd.l = EdgeList.leftEnd;
  338. EdgeList.list.unshift(EdgeList.leftEnd, EdgeList.rightEnd);
  339. },
  340. createHalfEdge: function(edge, side) {
  341. return {
  342. edge: edge,
  343. side: side,
  344. vertex: null,
  345. "l": null,
  346. "r": null
  347. };
  348. },
  349. insert: function(lb, he) {
  350. he.l = lb;
  351. he.r = lb.r;
  352. lb.r.l = he;
  353. lb.r = he;
  354. },
  355. leftBound: function(p) {
  356. var he = EdgeList.leftEnd;
  357. do {
  358. he = he.r;
  359. } while (he != EdgeList.rightEnd && Geom.rightOf(he, p));
  360. he = he.l;
  361. return he;
  362. },
  363. del: function(he) {
  364. he.l.r = he.r;
  365. he.r.l = he.l;
  366. he.edge = null;
  367. },
  368. right: function(he) {
  369. return he.r;
  370. },
  371. left: function(he) {
  372. return he.l;
  373. },
  374. leftRegion: function(he) {
  375. return he.edge == null
  376. ? Sites.bottomSite
  377. : he.edge.region[he.side];
  378. },
  379. rightRegion: function(he) {
  380. return he.edge == null
  381. ? Sites.bottomSite
  382. : he.edge.region[d3_voronoi_opposite[he.side]];
  383. }
  384. };
  385. var Geom = {
  386. bisect: function(s1, s2) {
  387. var newEdge = {
  388. region: {"l": s1, "r": s2},
  389. ep: {"l": null, "r": null}
  390. };
  391. var dx = s2.x - s1.x,
  392. dy = s2.y - s1.y,
  393. adx = dx > 0 ? dx : -dx,
  394. ady = dy > 0 ? dy : -dy;
  395. newEdge.c = s1.x * dx + s1.y * dy
  396. + (dx * dx + dy * dy) * .5;
  397. if (adx > ady) {
  398. newEdge.a = 1;
  399. newEdge.b = dy / dx;
  400. newEdge.c /= dx;
  401. } else {
  402. newEdge.b = 1;
  403. newEdge.a = dx / dy;
  404. newEdge.c /= dy;
  405. }
  406. return newEdge;
  407. },
  408. intersect: function(el1, el2) {
  409. var e1 = el1.edge,
  410. e2 = el2.edge;
  411. if (!e1 || !e2 || (e1.region.r == e2.region.r)) {
  412. return null;
  413. }
  414. var d = (e1.a * e2.b) - (e1.b * e2.a);
  415. if (Math.abs(d) < 1e-10) {
  416. return null;
  417. }
  418. var xint = (e1.c * e2.b - e2.c * e1.b) / d,
  419. yint = (e2.c * e1.a - e1.c * e2.a) / d,
  420. e1r = e1.region.r,
  421. e2r = e2.region.r,
  422. el,
  423. e;
  424. if ((e1r.y < e2r.y) ||
  425. (e1r.y == e2r.y && e1r.x < e2r.x)) {
  426. el = el1;
  427. e = e1;
  428. } else {
  429. el = el2;
  430. e = e2;
  431. }
  432. var rightOfSite = (xint >= e.region.r.x);
  433. if ((rightOfSite && (el.side === "l")) ||
  434. (!rightOfSite && (el.side === "r"))) {
  435. return null;
  436. }
  437. return {
  438. x: xint,
  439. y: yint
  440. };
  441. },
  442. rightOf: function(he, p) {
  443. var e = he.edge,
  444. topsite = e.region.r,
  445. rightOfSite = (p.x > topsite.x);
  446. if (rightOfSite && (he.side === "l")) {
  447. return 1;
  448. }
  449. if (!rightOfSite && (he.side === "r")) {
  450. return 0;
  451. }
  452. if (e.a === 1) {
  453. var dyp = p.y - topsite.y,
  454. dxp = p.x - topsite.x,
  455. fast = 0,
  456. above = 0;
  457. if ((!rightOfSite && (e.b < 0)) ||
  458. (rightOfSite && (e.b >= 0))) {
  459. above = fast = (dyp >= e.b * dxp);
  460. } else {
  461. above = ((p.x + p.y * e.b) > e.c);
  462. if (e.b < 0) {
  463. above = !above;
  464. }
  465. if (!above) {
  466. fast = 1;
  467. }
  468. }
  469. if (!fast) {
  470. var dxs = topsite.x - e.region.l.x;
  471. above = (e.b * (dxp * dxp - dyp * dyp)) <
  472. (dxs * dyp * (1 + 2 * dxp / dxs + e.b * e.b));
  473. if (e.b < 0) {
  474. above = !above;
  475. }
  476. }
  477. } else /* e.b == 1 */ {
  478. var yl = e.c - e.a * p.x,
  479. t1 = p.y - yl,
  480. t2 = p.x - topsite.x,
  481. t3 = yl - topsite.y;
  482. above = (t1 * t1) > (t2 * t2 + t3 * t3);
  483. }
  484. return he.side === "l" ? above : !above;
  485. },
  486. endPoint: function(edge, side, site) {
  487. edge.ep[side] = site;
  488. if (!edge.ep[d3_voronoi_opposite[side]]) return;
  489. callback(edge);
  490. },
  491. distance: function(s, t) {
  492. var dx = s.x - t.x,
  493. dy = s.y - t.y;
  494. return Math.sqrt(dx * dx + dy * dy);
  495. }
  496. };
  497. var EventQueue = {
  498. list: [],
  499. insert: function(he, site, offset) {
  500. he.vertex = site;
  501. he.ystar = site.y + offset;
  502. for (var i=0, list=EventQueue.list, l=list.length; i<l; i++) {
  503. var next = list[i];
  504. if (he.ystar > next.ystar ||
  505. (he.ystar == next.ystar &&
  506. site.x > next.vertex.x)) {
  507. continue;
  508. } else {
  509. break;
  510. }
  511. }
  512. list.splice(i, 0, he);
  513. },
  514. del: function(he) {
  515. for (var i=0, ls=EventQueue.list, l=ls.length; i<l && (ls[i] != he); ++i) {}
  516. ls.splice(i, 1);
  517. },
  518. empty: function() { return EventQueue.list.length === 0; },
  519. nextEvent: function(he) {
  520. for (var i=0, ls=EventQueue.list, l=ls.length; i<l; ++i) {
  521. if (ls[i] == he) return ls[i+1];
  522. }
  523. return null;
  524. },
  525. min: function() {
  526. var elem = EventQueue.list[0];
  527. return {
  528. x: elem.vertex.x,
  529. y: elem.ystar
  530. };
  531. },
  532. extractMin: function() {
  533. return EventQueue.list.shift();
  534. }
  535. };
  536. EdgeList.init();
  537. Sites.bottomSite = Sites.list.shift();
  538. var newSite = Sites.list.shift(), newIntStar;
  539. var lbnd, rbnd, llbnd, rrbnd, bisector;
  540. var bot, top, temp, p, v;
  541. var e, pm;
  542. while (true) {
  543. if (!EventQueue.empty()) {
  544. newIntStar = EventQueue.min();
  545. }
  546. if (newSite && (EventQueue.empty()
  547. || newSite.y < newIntStar.y
  548. || (newSite.y == newIntStar.y
  549. && newSite.x < newIntStar.x))) { //new site is smallest
  550. lbnd = EdgeList.leftBound(newSite);
  551. rbnd = EdgeList.right(lbnd);
  552. bot = EdgeList.rightRegion(lbnd);
  553. e = Geom.bisect(bot, newSite);
  554. bisector = EdgeList.createHalfEdge(e, "l");
  555. EdgeList.insert(lbnd, bisector);
  556. p = Geom.intersect(lbnd, bisector);
  557. if (p) {
  558. EventQueue.del(lbnd);
  559. EventQueue.insert(lbnd, p, Geom.distance(p, newSite));
  560. }
  561. lbnd = bisector;
  562. bisector = EdgeList.createHalfEdge(e, "r");
  563. EdgeList.insert(lbnd, bisector);
  564. p = Geom.intersect(bisector, rbnd);
  565. if (p) {
  566. EventQueue.insert(bisector, p, Geom.distance(p, newSite));
  567. }
  568. newSite = Sites.list.shift();
  569. } else if (!EventQueue.empty()) { //intersection is smallest
  570. lbnd = EventQueue.extractMin();
  571. llbnd = EdgeList.left(lbnd);
  572. rbnd = EdgeList.right(lbnd);
  573. rrbnd = EdgeList.right(rbnd);
  574. bot = EdgeList.leftRegion(lbnd);
  575. top = EdgeList.rightRegion(rbnd);
  576. v = lbnd.vertex;
  577. Geom.endPoint(lbnd.edge, lbnd.side, v);
  578. Geom.endPoint(rbnd.edge, rbnd.side, v);
  579. EdgeList.del(lbnd);
  580. EventQueue.del(rbnd);
  581. EdgeList.del(rbnd);
  582. pm = "l";
  583. if (bot.y > top.y) {
  584. temp = bot;
  585. bot = top;
  586. top = temp;
  587. pm = "r";
  588. }
  589. e = Geom.bisect(bot, top);
  590. bisector = EdgeList.createHalfEdge(e, pm);
  591. EdgeList.insert(llbnd, bisector);
  592. Geom.endPoint(e, d3_voronoi_opposite[pm], v);
  593. p = Geom.intersect(llbnd, bisector);
  594. if (p) {
  595. EventQueue.del(llbnd);
  596. EventQueue.insert(llbnd, p, Geom.distance(p, bot));
  597. }
  598. p = Geom.intersect(bisector, rrbnd);
  599. if (p) {
  600. EventQueue.insert(bisector, p, Geom.distance(p, bot));
  601. }
  602. } else {
  603. break;
  604. }
  605. }//end while
  606. for (lbnd = EdgeList.right(EdgeList.leftEnd);
  607. lbnd != EdgeList.rightEnd;
  608. lbnd = EdgeList.right(lbnd)) {
  609. callback(lbnd.edge);
  610. }
  611. }
  612. /**
  613. * @param vertices [[x1, y1], [x2, y2], …]
  614. * @returns triangles [[[x1, y1], [x2, y2], [x3, y3]], …]
  615. */
  616. d3.geom.delaunay = function(vertices) {
  617. var edges = vertices.map(function() { return []; }),
  618. triangles = [];
  619. // Use the Voronoi tessellation to determine Delaunay edges.
  620. d3_voronoi_tessellate(vertices, function(e) {
  621. edges[e.region.l.index].push(vertices[e.region.r.index]);
  622. });
  623. // Reconnect the edges into counterclockwise triangles.
  624. edges.forEach(function(edge, i) {
  625. var v = vertices[i],
  626. cx = v[0],
  627. cy = v[1];
  628. edge.forEach(function(v) {
  629. v.angle = Math.atan2(v[0] - cx, v[1] - cy);
  630. });
  631. edge.sort(function(a, b) {
  632. return a.angle - b.angle;
  633. });
  634. for (var j = 0, m = edge.length - 1; j < m; j++) {
  635. triangles.push([v, edge[j], edge[j + 1]]);
  636. }
  637. });
  638. return triangles;
  639. };
  640. // Constructs a new quadtree for the specified array of points. A quadtree is a
  641. // two-dimensional recursive spatial subdivision. This implementation uses
  642. // square partitions, dividing each square into four equally-sized squares. Each
  643. // point exists in a unique node; if multiple points are in the same position,
  644. // some points may be stored on internal nodes rather than leaf nodes. Quadtrees
  645. // can be used to accelerate various spatial operations, such as the Barnes-Hut
  646. // approximation for computing n-body forces, or collision detection.
  647. d3.geom.quadtree = function(points, x1, y1, x2, y2) {
  648. var p,
  649. i = -1,
  650. n = points.length;
  651. // Type conversion for deprecated API.
  652. if (n && isNaN(points[0].x)) points = points.map(d3_geom_quadtreePoint);
  653. // Allow bounds to be specified explicitly.
  654. if (arguments.length < 5) {
  655. if (arguments.length === 3) {
  656. y2 = x2 = y1;
  657. y1 = x1;
  658. } else {
  659. x1 = y1 = Infinity;
  660. x2 = y2 = -Infinity;
  661. // Compute bounds.
  662. while (++i < n) {
  663. p = points[i];
  664. if (p.x < x1) x1 = p.x;
  665. if (p.y < y1) y1 = p.y;
  666. if (p.x > x2) x2 = p.x;
  667. if (p.y > y2) y2 = p.y;
  668. }
  669. // Squarify the bounds.
  670. var dx = x2 - x1,
  671. dy = y2 - y1;
  672. if (dx > dy) y2 = y1 + dx;
  673. else x2 = x1 + dy;
  674. }
  675. }
  676. // Recursively inserts the specified point p at the node n or one of its
  677. // descendants. The bounds are defined by [x1, x2] and [y1, y2].
  678. function insert(n, p, x1, y1, x2, y2) {
  679. if (isNaN(p.x) || isNaN(p.y)) return; // ignore invalid points
  680. if (n.leaf) {
  681. var v = n.point;
  682. if (v) {
  683. // If the point at this leaf node is at the same position as the new
  684. // point we are adding, we leave the point associated with the
  685. // internal node while adding the new point to a child node. This
  686. // avoids infinite recursion.
  687. if ((Math.abs(v.x - p.x) + Math.abs(v.y - p.y)) < .01) {
  688. insertChild(n, p, x1, y1, x2, y2);
  689. } else {
  690. n.point = null;
  691. insertChild(n, v, x1, y1, x2, y2);
  692. insertChild(n, p, x1, y1, x2, y2);
  693. }
  694. } else {
  695. n.point = p;
  696. }
  697. } else {
  698. insertChild(n, p, x1, y1, x2, y2);
  699. }
  700. }
  701. // Recursively inserts the specified point p into a descendant of node n. The
  702. // bounds are defined by [x1, x2] and [y1, y2].
  703. function insertChild(n, p, x1, y1, x2, y2) {
  704. // Compute the split point, and the quadrant in which to insert p.
  705. var sx = (x1 + x2) * .5,
  706. sy = (y1 + y2) * .5,
  707. right = p.x >= sx,
  708. bottom = p.y >= sy,
  709. i = (bottom << 1) + right;
  710. // Recursively insert into the child node.
  711. n.leaf = false;
  712. n = n.nodes[i] || (n.nodes[i] = d3_geom_quadtreeNode());
  713. // Update the bounds as we recurse.
  714. if (right) x1 = sx; else x2 = sx;
  715. if (bottom) y1 = sy; else y2 = sy;
  716. insert(n, p, x1, y1, x2, y2);
  717. }
  718. // Create the root node.
  719. var root = d3_geom_quadtreeNode();
  720. root.add = function(p) {
  721. insert(root, p, x1, y1, x2, y2);
  722. };
  723. root.visit = function(f) {
  724. d3_geom_quadtreeVisit(f, root, x1, y1, x2, y2);
  725. };
  726. // Insert all points.
  727. points.forEach(root.add);
  728. return root;
  729. };
  730. function d3_geom_quadtreeNode() {
  731. return {
  732. leaf: true,
  733. nodes: [],
  734. point: null
  735. };
  736. }
  737. function d3_geom_quadtreeVisit(f, node, x1, y1, x2, y2) {
  738. if (!f(node, x1, y1, x2, y2)) {
  739. var sx = (x1 + x2) * .5,
  740. sy = (y1 + y2) * .5,
  741. children = node.nodes;
  742. if (children[0]) d3_geom_quadtreeVisit(f, children[0], x1, y1, sx, sy);
  743. if (children[1]) d3_geom_quadtreeVisit(f, children[1], sx, y1, x2, sy);
  744. if (children[2]) d3_geom_quadtreeVisit(f, children[2], x1, sy, sx, y2);
  745. if (children[3]) d3_geom_quadtreeVisit(f, children[3], sx, sy, x2, y2);
  746. }
  747. }
  748. function d3_geom_quadtreePoint(p) {
  749. return {
  750. x: p[0],
  751. y: p[1]
  752. };
  753. }
  754. })();