BoundedVoronoiLegacy.cs 24 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693
  1. // -----------------------------------------------------------------------
  2. // <copyright file="BoundedVoronoi.cs" company="">
  3. // Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
  4. // </copyright>
  5. // -----------------------------------------------------------------------
  6. namespace UnityEngine.U2D.Animation.TriangleNet
  7. .Voronoi.Legacy
  8. {
  9. using System;
  10. using System.Collections.Generic;
  11. using Animation.TriangleNet.Topology;
  12. using Animation.TriangleNet.Geometry;
  13. /// <summary>
  14. /// The Bounded Voronoi Diagram is the dual of a PSLG triangulation.
  15. /// </summary>
  16. /// <remarks>
  17. /// 2D Centroidal Voronoi Tessellations with Constraints, 2010,
  18. /// Jane Tournois, Pierre Alliez and Olivier Devillers
  19. /// </remarks>
  20. [Obsolete("Use TriangleNet.Voronoi.BoundedVoronoi class instead.")]
  21. internal class BoundedVoronoiLegacy : IVoronoi
  22. {
  23. IPredicates predicates = RobustPredicates.Default;
  24. Mesh mesh;
  25. Point[] points;
  26. List<VoronoiRegion> regions;
  27. // Used for new points on segments.
  28. List<Point> segPoints;
  29. int segIndex;
  30. Dictionary<int, SubSegment> subsegMap;
  31. bool includeBoundary = true;
  32. /// <summary>
  33. /// Initializes a new instance of the <see cref="BoundedVoronoiLegacy" /> class.
  34. /// </summary>
  35. /// <param name="mesh">Mesh instance.</param>
  36. public BoundedVoronoiLegacy(Mesh mesh)
  37. : this(mesh, true)
  38. {
  39. }
  40. /// <summary>
  41. /// Initializes a new instance of the <see cref="BoundedVoronoiLegacy" /> class.
  42. /// </summary>
  43. /// <param name="mesh">Mesh instance.</param>
  44. public BoundedVoronoiLegacy(Mesh mesh, bool includeBoundary)
  45. {
  46. this.mesh = mesh;
  47. this.includeBoundary = includeBoundary;
  48. Generate();
  49. }
  50. /// <summary>
  51. /// Gets the list of Voronoi vertices.
  52. /// </summary>
  53. public Point[] Points
  54. {
  55. get { return points; }
  56. }
  57. /// <summary>
  58. /// Gets the list of Voronoi regions.
  59. /// </summary>
  60. public ICollection<VoronoiRegion> Regions
  61. {
  62. get { return regions; }
  63. }
  64. public IEnumerable<IEdge> Edges
  65. {
  66. get { return EnumerateEdges(); }
  67. }
  68. /// <summary>
  69. /// Computes the bounded voronoi diagram.
  70. /// </summary>
  71. private void Generate()
  72. {
  73. mesh.Renumber();
  74. mesh.MakeVertexMap();
  75. // Allocate space for voronoi diagram
  76. this.regions = new List<VoronoiRegion>(mesh.vertices.Count);
  77. this.points = new Point[mesh.triangles.Count];
  78. this.segPoints = new List<Point>(mesh.subsegs.Count * 4);
  79. ComputeCircumCenters();
  80. TagBlindTriangles();
  81. foreach (var v in mesh.vertices.Values)
  82. {
  83. // TODO: Need a reliable way to check if a vertex is on a segment
  84. if (v.type == VertexType.FreeVertex || v.label == 0)
  85. {
  86. ConstructCell(v);
  87. }
  88. else if (includeBoundary)
  89. {
  90. ConstructBoundaryCell(v);
  91. }
  92. }
  93. // Add the new points on segments to the point array.
  94. int length = points.Length;
  95. Array.Resize<Point>(ref points, length + segPoints.Count);
  96. for (int i = 0; i < segPoints.Count; i++)
  97. {
  98. points[length + i] = segPoints[i];
  99. }
  100. this.segPoints.Clear();
  101. this.segPoints = null;
  102. }
  103. private void ComputeCircumCenters()
  104. {
  105. Otri tri = default(Otri);
  106. double xi = 0, eta = 0;
  107. Point pt;
  108. // Compue triangle circumcenters
  109. foreach (var item in mesh.triangles)
  110. {
  111. tri.tri = item;
  112. pt = predicates.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
  113. pt.id = item.id;
  114. points[item.id] = pt;
  115. }
  116. }
  117. /// <summary>
  118. /// Tag all blind triangles.
  119. /// </summary>
  120. /// <remarks>
  121. /// A triangle is said to be blind if the triangle and its circumcenter
  122. /// lie on two different sides of a constrained edge.
  123. /// </remarks>
  124. private void TagBlindTriangles()
  125. {
  126. int blinded = 0;
  127. Stack<Triangle> triangles;
  128. subsegMap = new Dictionary<int, SubSegment>();
  129. Otri f = default(Otri);
  130. Otri f0 = default(Otri);
  131. Osub e = default(Osub);
  132. Osub sub1 = default(Osub);
  133. // Tag all triangles non-blind
  134. foreach (var t in mesh.triangles)
  135. {
  136. // Use the infected flag for 'blinded' attribute.
  137. t.infected = false;
  138. }
  139. // for each constrained edge e of cdt do
  140. foreach (var ss in mesh.subsegs.Values)
  141. {
  142. // Create a stack: triangles
  143. triangles = new Stack<Triangle>();
  144. // for both adjacent triangles fe to e tagged non-blind do
  145. // Push fe into triangles
  146. e.seg = ss;
  147. e.orient = 0;
  148. e.Pivot(ref f);
  149. if (f.tri.id != Mesh.DUMMY && !f.tri.infected)
  150. {
  151. triangles.Push(f.tri);
  152. }
  153. e.Sym();
  154. e.Pivot(ref f);
  155. if (f.tri.id != Mesh.DUMMY && !f.tri.infected)
  156. {
  157. triangles.Push(f.tri);
  158. }
  159. // while triangles is non-empty
  160. while (triangles.Count > 0)
  161. {
  162. // Pop f from stack triangles
  163. f.tri = triangles.Pop();
  164. f.orient = 0;
  165. // if f is blinded by e (use P) then
  166. if (TriangleIsBlinded(ref f, ref e))
  167. {
  168. // Tag f as blinded by e
  169. f.tri.infected = true;
  170. blinded++;
  171. // Store association triangle -> subseg
  172. subsegMap.Add(f.tri.hash, e.seg);
  173. // for each adjacent triangle f0 to f do
  174. for (f.orient = 0; f.orient < 3; f.orient++)
  175. {
  176. f.Sym(ref f0);
  177. f0.Pivot(ref sub1);
  178. // if f0 is finite and tagged non-blind & the common edge
  179. // between f and f0 is unconstrained then
  180. if (f0.tri.id != Mesh.DUMMY && !f0.tri.infected && sub1.seg.hash == Mesh.DUMMY)
  181. {
  182. // Push f0 into triangles.
  183. triangles.Push(f0.tri);
  184. }
  185. }
  186. }
  187. }
  188. }
  189. blinded = 0;
  190. }
  191. /// <summary>
  192. /// Check if given triangle is blinded by given segment.
  193. /// </summary>
  194. /// <param name="tri">Triangle.</param>
  195. /// <param name="seg">Segments</param>
  196. /// <returns>Returns true, if the triangle is blinded.</returns>
  197. private bool TriangleIsBlinded(ref Otri tri, ref Osub seg)
  198. {
  199. Point c, pt;
  200. Vertex torg = tri.Org();
  201. Vertex tdest = tri.Dest();
  202. Vertex tapex = tri.Apex();
  203. Vertex sorg = seg.Org();
  204. Vertex sdest = seg.Dest();
  205. c = this.points[tri.tri.id];
  206. if (SegmentsIntersect(sorg, sdest, c, torg, out pt, true))
  207. {
  208. return true;
  209. }
  210. if (SegmentsIntersect(sorg, sdest, c, tdest, out pt, true))
  211. {
  212. return true;
  213. }
  214. if (SegmentsIntersect(sorg, sdest, c, tapex, out pt, true))
  215. {
  216. return true;
  217. }
  218. return false;
  219. }
  220. private void ConstructCell(Vertex vertex)
  221. {
  222. VoronoiRegion region = new VoronoiRegion(vertex);
  223. regions.Add(region);
  224. Otri f = default(Otri);
  225. Otri f_init = default(Otri);
  226. Otri f_next = default(Otri);
  227. Osub sf = default(Osub);
  228. Osub sfn = default(Osub);
  229. Point cc_f, cc_f_next, p;
  230. int n = mesh.triangles.Count;
  231. // Call P the polygon (cell) in construction
  232. List<Point> vpoints = new List<Point>();
  233. // Call f_init a triangle incident to x
  234. vertex.tri.Copy(ref f_init);
  235. if (f_init.Org() != vertex)
  236. {
  237. throw new Exception("ConstructCell: inconsistent topology.");
  238. }
  239. // Let f be initialized to f_init
  240. f_init.Copy(ref f);
  241. // Call f_next the next triangle counterclockwise around x
  242. f_init.Onext(ref f_next);
  243. // repeat ... until f = f_init
  244. do
  245. {
  246. // Call Lffnext the line going through the circumcenters of f and f_next
  247. cc_f = this.points[f.tri.id];
  248. cc_f_next = this.points[f_next.tri.id];
  249. // if f is tagged non-blind then
  250. if (!f.tri.infected)
  251. {
  252. // Insert the circumcenter of f into P
  253. vpoints.Add(cc_f);
  254. if (f_next.tri.infected)
  255. {
  256. // Call S_fnext the constrained edge blinding f_next
  257. sfn.seg = subsegMap[f_next.tri.hash];
  258. // Insert point Lf,f_next /\ Sf_next into P
  259. if (SegmentsIntersect(sfn.Org(), sfn.Dest(), cc_f, cc_f_next, out p, true))
  260. {
  261. p.id = n + segIndex++;
  262. segPoints.Add(p);
  263. vpoints.Add(p);
  264. }
  265. }
  266. }
  267. else
  268. {
  269. // Call Sf the constrained edge blinding f
  270. sf.seg = subsegMap[f.tri.hash];
  271. // if f_next is tagged non-blind then
  272. if (!f_next.tri.infected)
  273. {
  274. // Insert point Lf,f_next /\ Sf into P
  275. if (SegmentsIntersect(sf.Org(), sf.Dest(), cc_f, cc_f_next, out p, true))
  276. {
  277. p.id = n + segIndex++;
  278. segPoints.Add(p);
  279. vpoints.Add(p);
  280. }
  281. }
  282. else
  283. {
  284. // Call Sf_next the constrained edge blinding f_next
  285. sfn.seg = subsegMap[f_next.tri.hash];
  286. // if Sf != Sf_next then
  287. if (!sf.Equal(sfn))
  288. {
  289. // Insert Lf,fnext /\ Sf and Lf,fnext /\ Sfnext into P
  290. if (SegmentsIntersect(sf.Org(), sf.Dest(), cc_f, cc_f_next, out p, true))
  291. {
  292. p.id = n + segIndex++;
  293. segPoints.Add(p);
  294. vpoints.Add(p);
  295. }
  296. if (SegmentsIntersect(sfn.Org(), sfn.Dest(), cc_f, cc_f_next, out p, true))
  297. {
  298. p.id = n + segIndex++;
  299. segPoints.Add(p);
  300. vpoints.Add(p);
  301. }
  302. }
  303. }
  304. }
  305. // f <- f_next
  306. f_next.Copy(ref f);
  307. // Call f_next the next triangle counterclockwise around x
  308. f_next.Onext();
  309. }
  310. while (!f.Equals(f_init));
  311. // Output: Bounded Voronoi cell of x in counterclockwise order.
  312. region.Add(vpoints);
  313. }
  314. private void ConstructBoundaryCell(Vertex vertex)
  315. {
  316. VoronoiRegion region = new VoronoiRegion(vertex);
  317. regions.Add(region);
  318. Otri f = default(Otri);
  319. Otri f_init = default(Otri);
  320. Otri f_next = default(Otri);
  321. Otri f_prev = default(Otri);
  322. Osub sf = default(Osub);
  323. Osub sfn = default(Osub);
  324. Vertex torg, tdest, tapex, sorg, sdest;
  325. Point cc_f, cc_f_next, p;
  326. int n = mesh.triangles.Count;
  327. // Call P the polygon (cell) in construction
  328. List<Point> vpoints = new List<Point>();
  329. // Call f_init a triangle incident to x
  330. vertex.tri.Copy(ref f_init);
  331. if (f_init.Org() != vertex)
  332. {
  333. throw new Exception("ConstructBoundaryCell: inconsistent topology.");
  334. }
  335. // Let f be initialized to f_init
  336. f_init.Copy(ref f);
  337. // Call f_next the next triangle counterclockwise around x
  338. f_init.Onext(ref f_next);
  339. f_init.Oprev(ref f_prev);
  340. // Is the border to the left?
  341. if (f_prev.tri.id != Mesh.DUMMY)
  342. {
  343. // Go clockwise until we reach the border (or the initial triangle)
  344. while (f_prev.tri.id != Mesh.DUMMY && !f_prev.Equals(f_init))
  345. {
  346. f_prev.Copy(ref f);
  347. f_prev.Oprev();
  348. }
  349. f.Copy(ref f_init);
  350. f.Onext(ref f_next);
  351. }
  352. if (f_prev.tri.id == Mesh.DUMMY)
  353. {
  354. // For vertices on the domain boundaray, add the vertex. For
  355. // internal boundaries don't add it.
  356. p = new Point(vertex.x, vertex.y);
  357. p.id = n + segIndex++;
  358. segPoints.Add(p);
  359. vpoints.Add(p);
  360. }
  361. // Add midpoint of start triangles' edge.
  362. torg = f.Org();
  363. tdest = f.Dest();
  364. p = new Point((torg.x + tdest.x) / 2, (torg.y + tdest.y) / 2);
  365. p.id = n + segIndex++;
  366. segPoints.Add(p);
  367. vpoints.Add(p);
  368. // repeat ... until f = f_init
  369. do
  370. {
  371. // Call Lffnext the line going through the circumcenters of f and f_next
  372. cc_f = this.points[f.tri.id];
  373. if (f_next.tri.id == Mesh.DUMMY)
  374. {
  375. if (!f.tri.infected)
  376. {
  377. // Add last circumcenter
  378. vpoints.Add(cc_f);
  379. }
  380. // Add midpoint of last triangles' edge (chances are it has already
  381. // been added, so post process cell to remove duplicates???)
  382. torg = f.Org();
  383. tapex = f.Apex();
  384. p = new Point((torg.x + tapex.x) / 2, (torg.y + tapex.y) / 2);
  385. p.id = n + segIndex++;
  386. segPoints.Add(p);
  387. vpoints.Add(p);
  388. break;
  389. }
  390. cc_f_next = this.points[f_next.tri.id];
  391. // if f is tagged non-blind then
  392. if (!f.tri.infected)
  393. {
  394. // Insert the circumcenter of f into P
  395. vpoints.Add(cc_f);
  396. if (f_next.tri.infected)
  397. {
  398. // Call S_fnext the constrained edge blinding f_next
  399. sfn.seg = subsegMap[f_next.tri.hash];
  400. // Insert point Lf,f_next /\ Sf_next into P
  401. if (SegmentsIntersect(sfn.Org(), sfn.Dest(), cc_f, cc_f_next, out p, true))
  402. {
  403. p.id = n + segIndex++;
  404. segPoints.Add(p);
  405. vpoints.Add(p);
  406. }
  407. }
  408. }
  409. else
  410. {
  411. // Call Sf the constrained edge blinding f
  412. sf.seg = subsegMap[f.tri.hash];
  413. sorg = sf.Org();
  414. sdest = sf.Dest();
  415. // if f_next is tagged non-blind then
  416. if (!f_next.tri.infected)
  417. {
  418. tdest = f.Dest();
  419. tapex = f.Apex();
  420. // Both circumcenters lie on the blinded side, but we
  421. // have to add the intersection with the segment.
  422. // Center of f edge dest->apex
  423. Point bisec = new Point((tdest.x + tapex.x) / 2, (tdest.y + tapex.y) / 2);
  424. // Find intersection of seg with line through f's bisector and circumcenter
  425. if (SegmentsIntersect(sorg, sdest, bisec, cc_f, out p, false))
  426. {
  427. p.id = n + segIndex++;
  428. segPoints.Add(p);
  429. vpoints.Add(p);
  430. }
  431. // Insert point Lf,f_next /\ Sf into P
  432. if (SegmentsIntersect(sorg, sdest, cc_f, cc_f_next, out p, true))
  433. {
  434. p.id = n + segIndex++;
  435. segPoints.Add(p);
  436. vpoints.Add(p);
  437. }
  438. }
  439. else
  440. {
  441. // Call Sf_next the constrained edge blinding f_next
  442. sfn.seg = subsegMap[f_next.tri.hash];
  443. // if Sf != Sf_next then
  444. if (!sf.Equal(sfn))
  445. {
  446. // Insert Lf,fnext /\ Sf and Lf,fnext /\ Sfnext into P
  447. if (SegmentsIntersect(sorg, sdest, cc_f, cc_f_next, out p, true))
  448. {
  449. p.id = n + segIndex++;
  450. segPoints.Add(p);
  451. vpoints.Add(p);
  452. }
  453. if (SegmentsIntersect(sfn.Org(), sfn.Dest(), cc_f, cc_f_next, out p, true))
  454. {
  455. p.id = n + segIndex++;
  456. segPoints.Add(p);
  457. vpoints.Add(p);
  458. }
  459. }
  460. else
  461. {
  462. // Both circumcenters lie on the blinded side, but we
  463. // have to add the intersection with the segment.
  464. // Center of f_next edge org->dest
  465. Point bisec = new Point((torg.x + tdest.x) / 2, (torg.y + tdest.y) / 2);
  466. // Find intersection of seg with line through f_next's bisector and circumcenter
  467. if (SegmentsIntersect(sorg, sdest, bisec, cc_f_next, out p, false))
  468. {
  469. p.id = n + segIndex++;
  470. segPoints.Add(p);
  471. vpoints.Add(p);
  472. }
  473. }
  474. }
  475. }
  476. // f <- f_next
  477. f_next.Copy(ref f);
  478. // Call f_next the next triangle counterclockwise around x
  479. f_next.Onext();
  480. }
  481. while (!f.Equals(f_init));
  482. // Output: Bounded Voronoi cell of x in counterclockwise order.
  483. region.Add(vpoints);
  484. }
  485. /// <summary>
  486. /// Determines the intersection point of the line segment defined by points A and B with the
  487. /// line segment defined by points C and D.
  488. /// </summary>
  489. /// <param name="seg">The first segment AB.</param>
  490. /// <param name="pc">Endpoint C of second segment.</param>
  491. /// <param name="pd">Endpoint D of second segment.</param>
  492. /// <param name="p">Reference to the intersection point.</param>
  493. /// <param name="strictIntersect">If false, pa and pb represent a line.</param>
  494. /// <returns>Returns true if the intersection point was found, and stores that point in X,Y.
  495. /// Returns false if there is no determinable intersection point, in which case X,Y will
  496. /// be unmodified.
  497. /// </returns>
  498. private bool SegmentsIntersect(Point p1, Point p2, Point p3, Point p4, out Point p, bool strictIntersect)
  499. {
  500. p = null; // TODO: Voronoi SegmentsIntersect
  501. double Ax = p1.x, Ay = p1.y;
  502. double Bx = p2.x, By = p2.y;
  503. double Cx = p3.x, Cy = p3.y;
  504. double Dx = p4.x, Dy = p4.y;
  505. double distAB, theCos, theSin, newX, ABpos;
  506. // Fail if either line segment is zero-length.
  507. if (Ax == Bx && Ay == By || Cx == Dx && Cy == Dy) return false;
  508. // Fail if the segments share an end-point.
  509. if (Ax == Cx && Ay == Cy || Bx == Cx && By == Cy
  510. || Ax == Dx && Ay == Dy || Bx == Dx && By == Dy)
  511. {
  512. return false;
  513. }
  514. // (1) Translate the system so that point A is on the origin.
  515. Bx -= Ax; By -= Ay;
  516. Cx -= Ax; Cy -= Ay;
  517. Dx -= Ax; Dy -= Ay;
  518. // Discover the length of segment A-B.
  519. distAB = Math.Sqrt(Bx * Bx + By * By);
  520. // (2) Rotate the system so that point B is on the positive X axis.
  521. theCos = Bx / distAB;
  522. theSin = By / distAB;
  523. newX = Cx * theCos + Cy * theSin;
  524. Cy = Cy * theCos - Cx * theSin; Cx = newX;
  525. newX = Dx * theCos + Dy * theSin;
  526. Dy = Dy * theCos - Dx * theSin; Dx = newX;
  527. // Fail if segment C-D doesn't cross line A-B.
  528. if (Cy < 0 && Dy < 0 || Cy >= 0 && Dy >= 0 && strictIntersect) return false;
  529. // (3) Discover the position of the intersection point along line A-B.
  530. ABpos = Dx + (Cx - Dx) * Dy / (Dy - Cy);
  531. // Fail if segment C-D crosses line A-B outside of segment A-B.
  532. if (ABpos < 0 || ABpos > distAB && strictIntersect) return false;
  533. // (4) Apply the discovered position to line A-B in the original coordinate system.
  534. p = new Point(Ax + ABpos * theCos, Ay + ABpos * theSin);
  535. // Success.
  536. return true;
  537. }
  538. // TODO: Voronoi enumerate edges
  539. private IEnumerable<IEdge> EnumerateEdges()
  540. {
  541. // Copy edges
  542. Point first, last;
  543. var edges = new List<IEdge>(this.Regions.Count * 2);
  544. foreach (var region in this.Regions)
  545. {
  546. first = null;
  547. last = null;
  548. foreach (var pt in region.Vertices)
  549. {
  550. if (first == null)
  551. {
  552. first = pt;
  553. last = pt;
  554. }
  555. else
  556. {
  557. edges.Add(new Edge(last.id, pt.id));
  558. last = pt;
  559. }
  560. }
  561. if (region.Bounded && first != null)
  562. {
  563. edges.Add(new Edge(last.id, first.id));
  564. }
  565. }
  566. return edges;
  567. }
  568. }
  569. }