ConstraintMesher.cs 52 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232
  1. // -----------------------------------------------------------------------
  2. // <copyright file="ConstraintMesher.cs">
  3. // Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html
  4. // Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
  5. // </copyright>
  6. // -----------------------------------------------------------------------
  7. namespace UnityEngine.U2D.Animation.TriangleNet
  8. .Meshing
  9. {
  10. using System;
  11. using System.Collections.Generic;
  12. using Animation.TriangleNet.Geometry;
  13. using Animation.TriangleNet.Logging;
  14. using Animation.TriangleNet.Meshing.Iterators;
  15. using Animation.TriangleNet.Topology;
  16. internal class ConstraintMesher
  17. {
  18. IPredicates predicates;
  19. Mesh mesh;
  20. Behavior behavior;
  21. TriangleLocator locator;
  22. List<Triangle> viri;
  23. ILog<LogItem> logger;
  24. public ConstraintMesher(Mesh mesh, Configuration config)
  25. {
  26. this.mesh = mesh;
  27. this.predicates = config.Predicates();
  28. this.behavior = mesh.behavior;
  29. this.locator = mesh.locator;
  30. this.viri = new List<Triangle>();
  31. logger = Log.Instance;
  32. }
  33. /// <summary>
  34. /// Insert segments into the mesh.
  35. /// </summary>
  36. /// <param name="input">The polygon.</param>
  37. /// <param name="options">Constraint options.</param>
  38. public void Apply(IPolygon input, ConstraintOptions options)
  39. {
  40. behavior.Poly = input.Segments.Count > 0;
  41. // Copy constraint options
  42. if (options != null)
  43. {
  44. behavior.ConformingDelaunay = options.ConformingDelaunay;
  45. behavior.Convex = options.Convex;
  46. behavior.NoBisect = options.SegmentSplitting;
  47. if (behavior.ConformingDelaunay)
  48. {
  49. behavior.Quality = true;
  50. }
  51. }
  52. //if (input.EdgeMarkers != null)
  53. //{
  54. // behavior.UseBoundaryMarkers = true;
  55. //}
  56. behavior.useRegions = input.Regions.Count > 0;
  57. // Ensure that no vertex can be mistaken for a triangular bounding
  58. // box vertex in insertvertex().
  59. mesh.infvertex1 = null;
  60. mesh.infvertex2 = null;
  61. mesh.infvertex3 = null;
  62. if (behavior.useSegments)
  63. {
  64. // Segments will be introduced next.
  65. mesh.checksegments = true;
  66. // Insert PSLG segments and/or convex hull segments.
  67. FormSkeleton(input);
  68. }
  69. if (behavior.Poly && (mesh.triangles.Count > 0))
  70. {
  71. // Copy holes and regions
  72. mesh.holes.AddRange(input.Holes);
  73. mesh.regions.AddRange(input.Regions);
  74. // Carve out holes and concavities.
  75. CarveHoles();
  76. }
  77. }
  78. /// <summary>
  79. /// Find the holes and infect them. Find the area constraints and infect
  80. /// them. Infect the convex hull. Spread the infection and kill triangles.
  81. /// Spread the area constraints.
  82. /// </summary>
  83. private void CarveHoles()
  84. {
  85. Otri searchtri = default(Otri);
  86. Vertex searchorg, searchdest;
  87. LocateResult intersect;
  88. Triangle[] regionTris = null;
  89. var dummytri = mesh.dummytri;
  90. if (!mesh.behavior.Convex)
  91. {
  92. // Mark as infected any unprotected triangles on the boundary.
  93. // This is one way by which concavities are created.
  94. InfectHull();
  95. }
  96. if (!mesh.behavior.NoHoles)
  97. {
  98. // Infect each triangle in which a hole lies.
  99. foreach (var hole in mesh.holes)
  100. {
  101. // Ignore holes that aren't within the bounds of the mesh.
  102. if (mesh.bounds.Contains(hole))
  103. {
  104. // Start searching from some triangle on the outer boundary.
  105. searchtri.tri = dummytri;
  106. searchtri.orient = 0;
  107. searchtri.Sym();
  108. // Ensure that the hole is to the left of this boundary edge;
  109. // otherwise, locate() will falsely report that the hole
  110. // falls within the starting triangle.
  111. searchorg = searchtri.Org();
  112. searchdest = searchtri.Dest();
  113. if (predicates.CounterClockwise(searchorg, searchdest, hole) > 0.0)
  114. {
  115. // Find a triangle that contains the hole.
  116. intersect = mesh.locator.Locate(hole, ref searchtri);
  117. if ((intersect != LocateResult.Outside) && (!searchtri.IsInfected()))
  118. {
  119. // Infect the triangle. This is done by marking the triangle
  120. // as infected and including the triangle in the virus pool.
  121. searchtri.Infect();
  122. viri.Add(searchtri.tri);
  123. }
  124. }
  125. }
  126. }
  127. }
  128. // Now, we have to find all the regions BEFORE we carve the holes, because locate() won't
  129. // work when the triangulation is no longer convex. (Incidentally, this is the reason why
  130. // regional attributes and area constraints can't be used when refining a preexisting mesh,
  131. // which might not be convex; they can only be used with a freshly triangulated PSLG.)
  132. if (mesh.regions.Count > 0)
  133. {
  134. int i = 0;
  135. regionTris = new Triangle[mesh.regions.Count];
  136. // Find the starting triangle for each region.
  137. foreach (var region in mesh.regions)
  138. {
  139. regionTris[i] = dummytri;
  140. // Ignore region points that aren't within the bounds of the mesh.
  141. if (mesh.bounds.Contains(region.point))
  142. {
  143. // Start searching from some triangle on the outer boundary.
  144. searchtri.tri = dummytri;
  145. searchtri.orient = 0;
  146. searchtri.Sym();
  147. // Ensure that the region point is to the left of this boundary
  148. // edge; otherwise, locate() will falsely report that the
  149. // region point falls within the starting triangle.
  150. searchorg = searchtri.Org();
  151. searchdest = searchtri.Dest();
  152. if (predicates.CounterClockwise(searchorg, searchdest, region.point) > 0.0)
  153. {
  154. // Find a triangle that contains the region point.
  155. intersect = mesh.locator.Locate(region.point, ref searchtri);
  156. if ((intersect != LocateResult.Outside) && (!searchtri.IsInfected()))
  157. {
  158. // Record the triangle for processing after the
  159. // holes have been carved.
  160. regionTris[i] = searchtri.tri;
  161. regionTris[i].label = region.id;
  162. regionTris[i].area = region.area;
  163. }
  164. }
  165. }
  166. i++;
  167. }
  168. }
  169. if (viri.Count > 0)
  170. {
  171. // Carve the holes and concavities.
  172. Plague();
  173. }
  174. if (regionTris != null)
  175. {
  176. var iterator = new RegionIterator(mesh);
  177. for (int i = 0; i < regionTris.Length; i++)
  178. {
  179. if (regionTris[i].id != Mesh.DUMMY)
  180. {
  181. // Make sure the triangle under consideration still exists.
  182. // It may have been eaten by the virus.
  183. if (!Otri.IsDead(regionTris[i]))
  184. {
  185. // Apply one region's attribute and/or area constraint.
  186. iterator.Process(regionTris[i]);
  187. }
  188. }
  189. }
  190. }
  191. // Free up memory (virus pool should be empty anyway).
  192. viri.Clear();
  193. }
  194. /// <summary>
  195. /// Create the segments of a triangulation, including PSLG segments and edges
  196. /// on the convex hull.
  197. /// </summary>
  198. private void FormSkeleton(IPolygon input)
  199. {
  200. // The segment endpoints.
  201. Vertex p, q;
  202. mesh.insegments = 0;
  203. if (behavior.Poly)
  204. {
  205. // If the input vertices are collinear, there is no triangulation,
  206. // so don't try to insert segments.
  207. if (mesh.triangles.Count == 0)
  208. {
  209. return;
  210. }
  211. // If segments are to be inserted, compute a mapping
  212. // from vertices to triangles.
  213. if (input.Segments.Count > 0)
  214. {
  215. mesh.MakeVertexMap();
  216. }
  217. // Read and insert the segments.
  218. foreach (var seg in input.Segments)
  219. {
  220. mesh.insegments++;
  221. p = seg.GetVertex(0);
  222. q = seg.GetVertex(1);
  223. if ((p.x == q.x) && (p.y == q.y))
  224. {
  225. if (Log.Verbose)
  226. {
  227. logger.Warning("Endpoints of segment (IDs " + p.id + "/" + q.id + ") are coincident.",
  228. "Mesh.FormSkeleton()");
  229. }
  230. }
  231. else
  232. {
  233. InsertSegment(p, q, seg.Label);
  234. }
  235. }
  236. }
  237. if (behavior.Convex || !behavior.Poly)
  238. {
  239. // Enclose the convex hull with subsegments.
  240. MarkHull();
  241. }
  242. }
  243. #region Carving holes
  244. /// <summary>
  245. /// Virally infect all of the triangles of the convex hull that are not
  246. /// protected by subsegments. Where there are subsegments, set boundary
  247. /// markers as appropriate.
  248. /// </summary>
  249. private void InfectHull()
  250. {
  251. Otri hulltri = default(Otri);
  252. Otri nexttri = default(Otri);
  253. Otri starttri = default(Otri);
  254. Osub hullsubseg = default(Osub);
  255. Vertex horg, hdest;
  256. var dummytri = mesh.dummytri;
  257. // Find a triangle handle on the hull.
  258. hulltri.tri = dummytri;
  259. hulltri.orient = 0;
  260. hulltri.Sym();
  261. // Remember where we started so we know when to stop.
  262. hulltri.Copy(ref starttri);
  263. // Go once counterclockwise around the convex hull.
  264. do
  265. {
  266. // Ignore triangles that are already infected.
  267. if (!hulltri.IsInfected())
  268. {
  269. // Is the triangle protected by a subsegment?
  270. hulltri.Pivot(ref hullsubseg);
  271. if (hullsubseg.seg.hash == Mesh.DUMMY)
  272. {
  273. // The triangle is not protected; infect it.
  274. if (!hulltri.IsInfected())
  275. {
  276. hulltri.Infect();
  277. viri.Add(hulltri.tri);
  278. }
  279. }
  280. else
  281. {
  282. // The triangle is protected; set boundary markers if appropriate.
  283. if (hullsubseg.seg.boundary == 0)
  284. {
  285. hullsubseg.seg.boundary = 1;
  286. horg = hulltri.Org();
  287. hdest = hulltri.Dest();
  288. if (horg.label == 0)
  289. {
  290. horg.label = 1;
  291. }
  292. if (hdest.label == 0)
  293. {
  294. hdest.label = 1;
  295. }
  296. }
  297. }
  298. }
  299. // To find the next hull edge, go clockwise around the next vertex.
  300. hulltri.Lnext();
  301. hulltri.Oprev(ref nexttri);
  302. while (nexttri.tri.id != Mesh.DUMMY)
  303. {
  304. nexttri.Copy(ref hulltri);
  305. hulltri.Oprev(ref nexttri);
  306. }
  307. }
  308. while (!hulltri.Equals(starttri));
  309. }
  310. /// <summary>
  311. /// Spread the virus from all infected triangles to any neighbors not
  312. /// protected by subsegments. Delete all infected triangles.
  313. /// </summary>
  314. /// <remarks>
  315. /// This is the procedure that actually creates holes and concavities.
  316. ///
  317. /// This procedure operates in two phases. The first phase identifies all
  318. /// the triangles that will die, and marks them as infected. They are
  319. /// marked to ensure that each triangle is added to the virus pool only
  320. /// once, so the procedure will terminate.
  321. ///
  322. /// The second phase actually eliminates the infected triangles. It also
  323. /// eliminates orphaned vertices.
  324. /// </remarks>
  325. void Plague()
  326. {
  327. Otri testtri = default(Otri);
  328. Otri neighbor = default(Otri);
  329. Osub neighborsubseg = default(Osub);
  330. Vertex testvertex;
  331. Vertex norg, ndest;
  332. var dummysub = mesh.dummysub;
  333. var dummytri = mesh.dummytri;
  334. bool killorg;
  335. // Loop through all the infected triangles, spreading the virus to
  336. // their neighbors, then to their neighbors' neighbors.
  337. for (int i = 0; i < viri.Count; i++)
  338. {
  339. // WARNING: Don't use foreach, mesh.viri list may get modified.
  340. testtri.tri = viri[i];
  341. // A triangle is marked as infected by messing with one of its pointers
  342. // to subsegments, setting it to an illegal value. Hence, we have to
  343. // temporarily uninfect this triangle so that we can examine its
  344. // adjacent subsegments.
  345. // TODO: Not true in the C# version (so we could skip this).
  346. testtri.Uninfect();
  347. // Check each of the triangle's three neighbors.
  348. for (testtri.orient = 0; testtri.orient < 3; testtri.orient++)
  349. {
  350. // Find the neighbor.
  351. testtri.Sym(ref neighbor);
  352. // Check for a subsegment between the triangle and its neighbor.
  353. testtri.Pivot(ref neighborsubseg);
  354. // Check if the neighbor is nonexistent or already infected.
  355. if ((neighbor.tri.id == Mesh.DUMMY) || neighbor.IsInfected())
  356. {
  357. if (neighborsubseg.seg.hash != Mesh.DUMMY)
  358. {
  359. // There is a subsegment separating the triangle from its
  360. // neighbor, but both triangles are dying, so the subsegment
  361. // dies too.
  362. mesh.SubsegDealloc(neighborsubseg.seg);
  363. if (neighbor.tri.id != Mesh.DUMMY)
  364. {
  365. // Make sure the subsegment doesn't get deallocated again
  366. // later when the infected neighbor is visited.
  367. neighbor.Uninfect();
  368. neighbor.SegDissolve(dummysub);
  369. neighbor.Infect();
  370. }
  371. }
  372. }
  373. else
  374. { // The neighbor exists and is not infected.
  375. if (neighborsubseg.seg.hash == Mesh.DUMMY)
  376. {
  377. // There is no subsegment protecting the neighbor, so
  378. // the neighbor becomes infected.
  379. neighbor.Infect();
  380. // Ensure that the neighbor's neighbors will be infected.
  381. viri.Add(neighbor.tri);
  382. }
  383. else
  384. {
  385. // The neighbor is protected by a subsegment.
  386. // Remove this triangle from the subsegment.
  387. neighborsubseg.TriDissolve(dummytri);
  388. // The subsegment becomes a boundary. Set markers accordingly.
  389. if (neighborsubseg.seg.boundary == 0)
  390. {
  391. neighborsubseg.seg.boundary = 1;
  392. }
  393. norg = neighbor.Org();
  394. ndest = neighbor.Dest();
  395. if (norg.label == 0)
  396. {
  397. norg.label = 1;
  398. }
  399. if (ndest.label == 0)
  400. {
  401. ndest.label = 1;
  402. }
  403. }
  404. }
  405. }
  406. // Remark the triangle as infected, so it doesn't get added to the
  407. // virus pool again.
  408. testtri.Infect();
  409. }
  410. foreach (var virus in viri)
  411. {
  412. testtri.tri = virus;
  413. // Check each of the three corners of the triangle for elimination.
  414. // This is done by walking around each vertex, checking if it is
  415. // still connected to at least one live triangle.
  416. for (testtri.orient = 0; testtri.orient < 3; testtri.orient++)
  417. {
  418. testvertex = testtri.Org();
  419. // Check if the vertex has already been tested.
  420. if (testvertex != null)
  421. {
  422. killorg = true;
  423. // Mark the corner of the triangle as having been tested.
  424. testtri.SetOrg(null);
  425. // Walk counterclockwise about the vertex.
  426. testtri.Onext(ref neighbor);
  427. // Stop upon reaching a boundary or the starting triangle.
  428. while ((neighbor.tri.id != Mesh.DUMMY) &&
  429. (!neighbor.Equals(testtri)))
  430. {
  431. if (neighbor.IsInfected())
  432. {
  433. // Mark the corner of this triangle as having been tested.
  434. neighbor.SetOrg(null);
  435. }
  436. else
  437. {
  438. // A live triangle. The vertex survives.
  439. killorg = false;
  440. }
  441. // Walk counterclockwise about the vertex.
  442. neighbor.Onext();
  443. }
  444. // If we reached a boundary, we must walk clockwise as well.
  445. if (neighbor.tri.id == Mesh.DUMMY)
  446. {
  447. // Walk clockwise about the vertex.
  448. testtri.Oprev(ref neighbor);
  449. // Stop upon reaching a boundary.
  450. while (neighbor.tri.id != Mesh.DUMMY)
  451. {
  452. if (neighbor.IsInfected())
  453. {
  454. // Mark the corner of this triangle as having been tested.
  455. neighbor.SetOrg(null);
  456. }
  457. else
  458. {
  459. // A live triangle. The vertex survives.
  460. killorg = false;
  461. }
  462. // Walk clockwise about the vertex.
  463. neighbor.Oprev();
  464. }
  465. }
  466. if (killorg)
  467. {
  468. // Deleting vertex
  469. testvertex.type = VertexType.UndeadVertex;
  470. mesh.undeads++;
  471. }
  472. }
  473. }
  474. // Record changes in the number of boundary edges, and disconnect
  475. // dead triangles from their neighbors.
  476. for (testtri.orient = 0; testtri.orient < 3; testtri.orient++)
  477. {
  478. testtri.Sym(ref neighbor);
  479. if (neighbor.tri.id == Mesh.DUMMY)
  480. {
  481. // There is no neighboring triangle on this edge, so this edge
  482. // is a boundary edge. This triangle is being deleted, so this
  483. // boundary edge is deleted.
  484. mesh.hullsize--;
  485. }
  486. else
  487. {
  488. // Disconnect the triangle from its neighbor.
  489. neighbor.Dissolve(dummytri);
  490. // There is a neighboring triangle on this edge, so this edge
  491. // becomes a boundary edge when this triangle is deleted.
  492. mesh.hullsize++;
  493. }
  494. }
  495. // Return the dead triangle to the pool of triangles.
  496. mesh.TriangleDealloc(testtri.tri);
  497. }
  498. // Empty the virus pool.
  499. viri.Clear();
  500. }
  501. #endregion
  502. #region Segment insertion
  503. /// <summary>
  504. /// Find the first triangle on the path from one point to another.
  505. /// </summary>
  506. /// <param name="searchtri"></param>
  507. /// <param name="searchpoint"></param>
  508. /// <returns>
  509. /// The return value notes whether the destination or apex of the found
  510. /// triangle is collinear with the two points in question.</returns>
  511. /// <remarks>
  512. /// Finds the triangle that intersects a line segment drawn from the
  513. /// origin of 'searchtri' to the point 'searchpoint', and returns the result
  514. /// in 'searchtri'. The origin of 'searchtri' does not change, even though
  515. /// the triangle returned may differ from the one passed in. This routine
  516. /// is used to find the direction to move in to get from one point to
  517. /// another.
  518. /// </remarks>
  519. private FindDirectionResult FindDirection(ref Otri searchtri, Vertex searchpoint)
  520. {
  521. Otri checktri = default(Otri);
  522. Vertex startvertex;
  523. Vertex leftvertex, rightvertex;
  524. double leftccw, rightccw;
  525. bool leftflag, rightflag;
  526. startvertex = searchtri.Org();
  527. rightvertex = searchtri.Dest();
  528. leftvertex = searchtri.Apex();
  529. // Is 'searchpoint' to the left?
  530. leftccw = predicates.CounterClockwise(searchpoint, startvertex, leftvertex);
  531. leftflag = leftccw > 0.0;
  532. // Is 'searchpoint' to the right?
  533. rightccw = predicates.CounterClockwise(startvertex, searchpoint, rightvertex);
  534. rightflag = rightccw > 0.0;
  535. if (leftflag && rightflag)
  536. {
  537. // 'searchtri' faces directly away from 'searchpoint'. We could go left
  538. // or right. Ask whether it's a triangle or a boundary on the left.
  539. searchtri.Onext(ref checktri);
  540. if (checktri.tri.id == Mesh.DUMMY)
  541. {
  542. leftflag = false;
  543. }
  544. else
  545. {
  546. rightflag = false;
  547. }
  548. }
  549. while (leftflag)
  550. {
  551. // Turn left until satisfied.
  552. searchtri.Onext();
  553. if (searchtri.tri.id == Mesh.DUMMY)
  554. {
  555. logger.Error("Unable to find a triangle on path.", "Mesh.FindDirection().1");
  556. throw new Exception("Unable to find a triangle on path.");
  557. }
  558. leftvertex = searchtri.Apex();
  559. rightccw = leftccw;
  560. leftccw = predicates.CounterClockwise(searchpoint, startvertex, leftvertex);
  561. leftflag = leftccw > 0.0;
  562. }
  563. while (rightflag)
  564. {
  565. // Turn right until satisfied.
  566. searchtri.Oprev();
  567. if (searchtri.tri.id == Mesh.DUMMY)
  568. {
  569. logger.Error("Unable to find a triangle on path.", "Mesh.FindDirection().2");
  570. throw new Exception("Unable to find a triangle on path.");
  571. }
  572. rightvertex = searchtri.Dest();
  573. leftccw = rightccw;
  574. rightccw = predicates.CounterClockwise(startvertex, searchpoint, rightvertex);
  575. rightflag = rightccw > 0.0;
  576. }
  577. if (leftccw == 0.0)
  578. {
  579. return FindDirectionResult.Leftcollinear;
  580. }
  581. else if (rightccw == 0.0)
  582. {
  583. return FindDirectionResult.Rightcollinear;
  584. }
  585. else
  586. {
  587. return FindDirectionResult.Within;
  588. }
  589. }
  590. /// <summary>
  591. /// Find the intersection of an existing segment and a segment that is being
  592. /// inserted. Insert a vertex at the intersection, splitting an existing subsegment.
  593. /// </summary>
  594. /// <param name="splittri"></param>
  595. /// <param name="splitsubseg"></param>
  596. /// <param name="endpoint2"></param>
  597. /// <remarks>
  598. /// The segment being inserted connects the apex of splittri to endpoint2.
  599. /// splitsubseg is the subsegment being split, and MUST adjoin splittri.
  600. /// Hence, endpoints of the subsegment being split are the origin and
  601. /// destination of splittri.
  602. ///
  603. /// On completion, splittri is a handle having the newly inserted
  604. /// intersection point as its origin, and endpoint1 as its destination.
  605. /// </remarks>
  606. private void SegmentIntersection(ref Otri splittri, ref Osub splitsubseg, Vertex endpoint2)
  607. {
  608. Osub opposubseg = default(Osub);
  609. Vertex endpoint1;
  610. Vertex torg, tdest;
  611. Vertex leftvertex, rightvertex;
  612. Vertex newvertex;
  613. InsertVertexResult success;
  614. var dummysub = mesh.dummysub;
  615. double ex, ey;
  616. double tx, ty;
  617. double etx, ety;
  618. double split, denom;
  619. // Find the other three segment endpoints.
  620. endpoint1 = splittri.Apex();
  621. torg = splittri.Org();
  622. tdest = splittri.Dest();
  623. // Segment intersection formulae; see the Antonio reference.
  624. tx = tdest.x - torg.x;
  625. ty = tdest.y - torg.y;
  626. ex = endpoint2.x - endpoint1.x;
  627. ey = endpoint2.y - endpoint1.y;
  628. etx = torg.x - endpoint2.x;
  629. ety = torg.y - endpoint2.y;
  630. denom = ty * ex - tx * ey;
  631. if (denom == 0.0)
  632. {
  633. logger.Error("Attempt to find intersection of parallel segments.",
  634. "Mesh.SegmentIntersection()");
  635. throw new Exception("Attempt to find intersection of parallel segments.");
  636. }
  637. split = (ey * etx - ex * ety) / denom;
  638. // Create the new vertex.
  639. newvertex = new Vertex(
  640. torg.x + split * (tdest.x - torg.x),
  641. torg.y + split * (tdest.y - torg.y),
  642. splitsubseg.seg.boundary
  643. #if USE_ATTRIBS
  644. , mesh.nextras
  645. #endif
  646. );
  647. newvertex.hash = mesh.hash_vtx++;
  648. newvertex.id = newvertex.hash;
  649. #if USE_ATTRIBS
  650. // Interpolate its attributes.
  651. for (int i = 0; i < mesh.nextras; i++)
  652. {
  653. newvertex.attributes[i] = torg.attributes[i] + split * (tdest.attributes[i] - torg.attributes[i]);
  654. }
  655. #endif
  656. #if USE_Z
  657. newvertex.z = torg.z + split * (tdest.z - torg.z);
  658. #endif
  659. mesh.vertices.Add(newvertex.hash, newvertex);
  660. // Insert the intersection vertex. This should always succeed.
  661. success = mesh.InsertVertex(newvertex, ref splittri, ref splitsubseg, false, false);
  662. if (success != InsertVertexResult.Successful)
  663. {
  664. logger.Error("Failure to split a segment.", "Mesh.SegmentIntersection()");
  665. throw new Exception("Failure to split a segment.");
  666. }
  667. // Record a triangle whose origin is the new vertex.
  668. newvertex.tri = splittri;
  669. if (mesh.steinerleft > 0)
  670. {
  671. mesh.steinerleft--;
  672. }
  673. // Divide the segment into two, and correct the segment endpoints.
  674. splitsubseg.Sym();
  675. splitsubseg.Pivot(ref opposubseg);
  676. splitsubseg.Dissolve(dummysub);
  677. opposubseg.Dissolve(dummysub);
  678. do
  679. {
  680. splitsubseg.SetSegOrg(newvertex);
  681. splitsubseg.Next();
  682. }
  683. while (splitsubseg.seg.hash != Mesh.DUMMY);
  684. do
  685. {
  686. opposubseg.SetSegOrg(newvertex);
  687. opposubseg.Next();
  688. }
  689. while (opposubseg.seg.hash != Mesh.DUMMY);
  690. // Inserting the vertex may have caused edge flips. We wish to rediscover
  691. // the edge connecting endpoint1 to the new intersection vertex.
  692. FindDirection(ref splittri, endpoint1);
  693. rightvertex = splittri.Dest();
  694. leftvertex = splittri.Apex();
  695. if ((leftvertex.x == endpoint1.x) && (leftvertex.y == endpoint1.y))
  696. {
  697. splittri.Onext();
  698. }
  699. else if ((rightvertex.x != endpoint1.x) || (rightvertex.y != endpoint1.y))
  700. {
  701. logger.Error("Topological inconsistency after splitting a segment.", "Mesh.SegmentIntersection()");
  702. throw new Exception("Topological inconsistency after splitting a segment.");
  703. }
  704. // 'splittri' should have destination endpoint1.
  705. }
  706. /// <summary>
  707. /// Scout the first triangle on the path from one endpoint to another, and check
  708. /// for completion (reaching the second endpoint), a collinear vertex, or the
  709. /// intersection of two segments.
  710. /// </summary>
  711. /// <param name="searchtri"></param>
  712. /// <param name="endpoint2"></param>
  713. /// <param name="newmark"></param>
  714. /// <returns>Returns true if the entire segment is successfully inserted, and false
  715. /// if the job must be finished by ConstrainedEdge().</returns>
  716. /// <remarks>
  717. /// If the first triangle on the path has the second endpoint as its
  718. /// destination or apex, a subsegment is inserted and the job is done.
  719. ///
  720. /// If the first triangle on the path has a destination or apex that lies on
  721. /// the segment, a subsegment is inserted connecting the first endpoint to
  722. /// the collinear vertex, and the search is continued from the collinear
  723. /// vertex.
  724. ///
  725. /// If the first triangle on the path has a subsegment opposite its origin,
  726. /// then there is a segment that intersects the segment being inserted.
  727. /// Their intersection vertex is inserted, splitting the subsegment.
  728. /// </remarks>
  729. private bool ScoutSegment(ref Otri searchtri, Vertex endpoint2, int newmark)
  730. {
  731. Otri crosstri = default(Otri);
  732. Osub crosssubseg = default(Osub);
  733. Vertex leftvertex, rightvertex;
  734. FindDirectionResult collinear;
  735. collinear = FindDirection(ref searchtri, endpoint2);
  736. rightvertex = searchtri.Dest();
  737. leftvertex = searchtri.Apex();
  738. if (((leftvertex.x == endpoint2.x) && (leftvertex.y == endpoint2.y)) ||
  739. ((rightvertex.x == endpoint2.x) && (rightvertex.y == endpoint2.y)))
  740. {
  741. // The segment is already an edge in the mesh.
  742. if ((leftvertex.x == endpoint2.x) && (leftvertex.y == endpoint2.y))
  743. {
  744. searchtri.Lprev();
  745. }
  746. // Insert a subsegment, if there isn't already one there.
  747. mesh.InsertSubseg(ref searchtri, newmark);
  748. return true;
  749. }
  750. else if (collinear == FindDirectionResult.Leftcollinear)
  751. {
  752. // We've collided with a vertex between the segment's endpoints.
  753. // Make the collinear vertex be the triangle's origin.
  754. searchtri.Lprev();
  755. mesh.InsertSubseg(ref searchtri, newmark);
  756. // Insert the remainder of the segment.
  757. return ScoutSegment(ref searchtri, endpoint2, newmark);
  758. }
  759. else if (collinear == FindDirectionResult.Rightcollinear)
  760. {
  761. // We've collided with a vertex between the segment's endpoints.
  762. mesh.InsertSubseg(ref searchtri, newmark);
  763. // Make the collinear vertex be the triangle's origin.
  764. searchtri.Lnext();
  765. // Insert the remainder of the segment.
  766. return ScoutSegment(ref searchtri, endpoint2, newmark);
  767. }
  768. else
  769. {
  770. searchtri.Lnext(ref crosstri);
  771. crosstri.Pivot(ref crosssubseg);
  772. // Check for a crossing segment.
  773. if (crosssubseg.seg.hash == Mesh.DUMMY)
  774. {
  775. return false;
  776. }
  777. else
  778. {
  779. // Insert a vertex at the intersection.
  780. SegmentIntersection(ref crosstri, ref crosssubseg, endpoint2);
  781. crosstri.Copy(ref searchtri);
  782. mesh.InsertSubseg(ref searchtri, newmark);
  783. // Insert the remainder of the segment.
  784. return ScoutSegment(ref searchtri, endpoint2, newmark);
  785. }
  786. }
  787. }
  788. /// <summary>
  789. /// Enforce the Delaunay condition at an edge, fanning out recursively from
  790. /// an existing vertex. Pay special attention to stacking inverted triangles.
  791. /// </summary>
  792. /// <param name="fixuptri"></param>
  793. /// <param name="leftside">Indicates whether or not fixuptri is to the left of
  794. /// the segment being inserted. (Imagine that the segment is pointing up from
  795. /// endpoint1 to endpoint2.)</param>
  796. /// <remarks>
  797. /// This is a support routine for inserting segments into a constrained
  798. /// Delaunay triangulation.
  799. ///
  800. /// The origin of fixuptri is treated as if it has just been inserted, and
  801. /// the local Delaunay condition needs to be enforced. It is only enforced
  802. /// in one sector, however, that being the angular range defined by
  803. /// fixuptri.
  804. ///
  805. /// This routine also needs to make decisions regarding the "stacking" of
  806. /// triangles. (Read the description of ConstrainedEdge() below before
  807. /// reading on here, so you understand the algorithm.) If the position of
  808. /// the new vertex (the origin of fixuptri) indicates that the vertex before
  809. /// it on the polygon is a reflex vertex, then "stack" the triangle by
  810. /// doing nothing. (fixuptri is an inverted triangle, which is how stacked
  811. /// triangles are identified.)
  812. ///
  813. /// Otherwise, check whether the vertex before that was a reflex vertex.
  814. /// If so, perform an edge flip, thereby eliminating an inverted triangle
  815. /// (popping it off the stack). The edge flip may result in the creation
  816. /// of a new inverted triangle, depending on whether or not the new vertex
  817. /// is visible to the vertex three edges behind on the polygon.
  818. ///
  819. /// If neither of the two vertices behind the new vertex are reflex
  820. /// vertices, fixuptri and fartri, the triangle opposite it, are not
  821. /// inverted; hence, ensure that the edge between them is locally Delaunay.
  822. /// </remarks>
  823. private void DelaunayFixup(ref Otri fixuptri, bool leftside)
  824. {
  825. Otri neartri = default(Otri);
  826. Otri fartri = default(Otri);
  827. Osub faredge = default(Osub);
  828. Vertex nearvertex, leftvertex, rightvertex, farvertex;
  829. fixuptri.Lnext(ref neartri);
  830. neartri.Sym(ref fartri);
  831. // Check if the edge opposite the origin of fixuptri can be flipped.
  832. if (fartri.tri.id == Mesh.DUMMY)
  833. {
  834. return;
  835. }
  836. neartri.Pivot(ref faredge);
  837. if (faredge.seg.hash != Mesh.DUMMY)
  838. {
  839. return;
  840. }
  841. // Find all the relevant vertices.
  842. nearvertex = neartri.Apex();
  843. leftvertex = neartri.Org();
  844. rightvertex = neartri.Dest();
  845. farvertex = fartri.Apex();
  846. // Check whether the previous polygon vertex is a reflex vertex.
  847. if (leftside)
  848. {
  849. if (predicates.CounterClockwise(nearvertex, leftvertex, farvertex) <= 0.0)
  850. {
  851. // leftvertex is a reflex vertex too. Nothing can
  852. // be done until a convex section is found.
  853. return;
  854. }
  855. }
  856. else
  857. {
  858. if (predicates.CounterClockwise(farvertex, rightvertex, nearvertex) <= 0.0)
  859. {
  860. // rightvertex is a reflex vertex too. Nothing can
  861. // be done until a convex section is found.
  862. return;
  863. }
  864. }
  865. if (predicates.CounterClockwise(rightvertex, leftvertex, farvertex) > 0.0)
  866. {
  867. // fartri is not an inverted triangle, and farvertex is not a reflex
  868. // vertex. As there are no reflex vertices, fixuptri isn't an
  869. // inverted triangle, either. Hence, test the edge between the
  870. // triangles to ensure it is locally Delaunay.
  871. if (predicates.InCircle(leftvertex, farvertex, rightvertex, nearvertex) <= 0.0)
  872. {
  873. return;
  874. }
  875. // Not locally Delaunay; go on to an edge flip.
  876. }
  877. // else fartri is inverted; remove it from the stack by flipping.
  878. mesh.Flip(ref neartri);
  879. fixuptri.Lprev(); // Restore the origin of fixuptri after the flip.
  880. // Recursively process the two triangles that result from the flip.
  881. DelaunayFixup(ref fixuptri, leftside);
  882. DelaunayFixup(ref fartri, leftside);
  883. }
  884. /// <summary>
  885. /// Force a segment into a constrained Delaunay triangulation by deleting the
  886. /// triangles it intersects, and triangulating the polygons that form on each
  887. /// side of it.
  888. /// </summary>
  889. /// <param name="starttri"></param>
  890. /// <param name="endpoint2"></param>
  891. /// <param name="newmark"></param>
  892. /// <remarks>
  893. /// Generates a single subsegment connecting 'endpoint1' to 'endpoint2'.
  894. /// The triangle 'starttri' has 'endpoint1' as its origin. 'newmark' is the
  895. /// boundary marker of the segment.
  896. ///
  897. /// To insert a segment, every triangle whose interior intersects the
  898. /// segment is deleted. The union of these deleted triangles is a polygon
  899. /// (which is not necessarily monotone, but is close enough), which is
  900. /// divided into two polygons by the new segment. This routine's task is
  901. /// to generate the Delaunay triangulation of these two polygons.
  902. ///
  903. /// You might think of this routine's behavior as a two-step process. The
  904. /// first step is to walk from endpoint1 to endpoint2, flipping each edge
  905. /// encountered. This step creates a fan of edges connected to endpoint1,
  906. /// including the desired edge to endpoint2. The second step enforces the
  907. /// Delaunay condition on each side of the segment in an incremental manner:
  908. /// proceeding along the polygon from endpoint1 to endpoint2 (this is done
  909. /// independently on each side of the segment), each vertex is "enforced"
  910. /// as if it had just been inserted, but affecting only the previous
  911. /// vertices. The result is the same as if the vertices had been inserted
  912. /// in the order they appear on the polygon, so the result is Delaunay.
  913. ///
  914. /// In truth, ConstrainedEdge() interleaves these two steps. The procedure
  915. /// walks from endpoint1 to endpoint2, and each time an edge is encountered
  916. /// and flipped, the newly exposed vertex (at the far end of the flipped
  917. /// edge) is "enforced" upon the previously flipped edges, usually affecting
  918. /// only one side of the polygon (depending upon which side of the segment
  919. /// the vertex falls on).
  920. ///
  921. /// The algorithm is complicated by the need to handle polygons that are not
  922. /// convex. Although the polygon is not necessarily monotone, it can be
  923. /// triangulated in a manner similar to the stack-based algorithms for
  924. /// monotone polygons. For each reflex vertex (local concavity) of the
  925. /// polygon, there will be an inverted triangle formed by one of the edge
  926. /// flips. (An inverted triangle is one with negative area - that is, its
  927. /// vertices are arranged in clockwise order - and is best thought of as a
  928. /// wrinkle in the fabric of the mesh.) Each inverted triangle can be
  929. /// thought of as a reflex vertex pushed on the stack, waiting to be fixed
  930. /// later.
  931. ///
  932. /// A reflex vertex is popped from the stack when a vertex is inserted that
  933. /// is visible to the reflex vertex. (However, if the vertex behind the
  934. /// reflex vertex is not visible to the reflex vertex, a new inverted
  935. /// triangle will take its place on the stack.) These details are handled
  936. /// by the DelaunayFixup() routine above.
  937. /// </remarks>
  938. private void ConstrainedEdge(ref Otri starttri, Vertex endpoint2, int newmark)
  939. {
  940. Otri fixuptri = default(Otri), fixuptri2 = default(Otri);
  941. Osub crosssubseg = default(Osub);
  942. Vertex endpoint1;
  943. Vertex farvertex;
  944. double area;
  945. bool collision;
  946. bool done;
  947. endpoint1 = starttri.Org();
  948. starttri.Lnext(ref fixuptri);
  949. mesh.Flip(ref fixuptri);
  950. // 'collision' indicates whether we have found a vertex directly
  951. // between endpoint1 and endpoint2.
  952. collision = false;
  953. done = false;
  954. do
  955. {
  956. farvertex = fixuptri.Org();
  957. // 'farvertex' is the extreme point of the polygon we are "digging"
  958. // to get from endpoint1 to endpoint2.
  959. if ((farvertex.x == endpoint2.x) && (farvertex.y == endpoint2.y))
  960. {
  961. fixuptri.Oprev(ref fixuptri2);
  962. // Enforce the Delaunay condition around endpoint2.
  963. DelaunayFixup(ref fixuptri, false);
  964. DelaunayFixup(ref fixuptri2, true);
  965. done = true;
  966. }
  967. else
  968. {
  969. // Check whether farvertex is to the left or right of the segment being
  970. // inserted, to decide which edge of fixuptri to dig through next.
  971. area = predicates.CounterClockwise(endpoint1, endpoint2, farvertex);
  972. if (area == 0.0)
  973. {
  974. // We've collided with a vertex between endpoint1 and endpoint2.
  975. collision = true;
  976. fixuptri.Oprev(ref fixuptri2);
  977. // Enforce the Delaunay condition around farvertex.
  978. DelaunayFixup(ref fixuptri, false);
  979. DelaunayFixup(ref fixuptri2, true);
  980. done = true;
  981. }
  982. else
  983. {
  984. if (area > 0.0)
  985. {
  986. // farvertex is to the left of the segment.
  987. fixuptri.Oprev(ref fixuptri2);
  988. // Enforce the Delaunay condition around farvertex, on the
  989. // left side of the segment only.
  990. DelaunayFixup(ref fixuptri2, true);
  991. // Flip the edge that crosses the segment. After the edge is
  992. // flipped, one of its endpoints is the fan vertex, and the
  993. // destination of fixuptri is the fan vertex.
  994. fixuptri.Lprev();
  995. }
  996. else
  997. {
  998. // farvertex is to the right of the segment.
  999. DelaunayFixup(ref fixuptri, false);
  1000. // Flip the edge that crosses the segment. After the edge is
  1001. // flipped, one of its endpoints is the fan vertex, and the
  1002. // destination of fixuptri is the fan vertex.
  1003. fixuptri.Oprev();
  1004. }
  1005. // Check for two intersecting segments.
  1006. fixuptri.Pivot(ref crosssubseg);
  1007. if (crosssubseg.seg.hash == Mesh.DUMMY)
  1008. {
  1009. mesh.Flip(ref fixuptri); // May create inverted triangle at left.
  1010. }
  1011. else
  1012. {
  1013. // We've collided with a segment between endpoint1 and endpoint2.
  1014. collision = true;
  1015. // Insert a vertex at the intersection.
  1016. SegmentIntersection(ref fixuptri, ref crosssubseg, endpoint2);
  1017. done = true;
  1018. }
  1019. }
  1020. }
  1021. }
  1022. while (!done);
  1023. // Insert a subsegment to make the segment permanent.
  1024. mesh.InsertSubseg(ref fixuptri, newmark);
  1025. // If there was a collision with an interceding vertex, install another
  1026. // segment connecting that vertex with endpoint2.
  1027. if (collision)
  1028. {
  1029. // Insert the remainder of the segment.
  1030. if (!ScoutSegment(ref fixuptri, endpoint2, newmark))
  1031. {
  1032. ConstrainedEdge(ref fixuptri, endpoint2, newmark);
  1033. }
  1034. }
  1035. }
  1036. /// <summary>
  1037. /// Insert a PSLG segment into a triangulation.
  1038. /// </summary>
  1039. /// <param name="endpoint1"></param>
  1040. /// <param name="endpoint2"></param>
  1041. /// <param name="newmark"></param>
  1042. private void InsertSegment(Vertex endpoint1, Vertex endpoint2, int newmark)
  1043. {
  1044. Otri searchtri1 = default(Otri), searchtri2 = default(Otri);
  1045. Vertex checkvertex = null;
  1046. var dummytri = mesh.dummytri;
  1047. // Find a triangle whose origin is the segment's first endpoint.
  1048. searchtri1 = endpoint1.tri;
  1049. if (searchtri1.tri != null)
  1050. {
  1051. checkvertex = searchtri1.Org();
  1052. }
  1053. if (checkvertex != endpoint1)
  1054. {
  1055. // Find a boundary triangle to search from.
  1056. searchtri1.tri = dummytri;
  1057. searchtri1.orient = 0;
  1058. searchtri1.Sym();
  1059. // Search for the segment's first endpoint by point location.
  1060. if (locator.Locate(endpoint1, ref searchtri1) != LocateResult.OnVertex)
  1061. {
  1062. logger.Error("Unable to locate PSLG vertex in triangulation.", "Mesh.InsertSegment().1");
  1063. throw new Exception("Unable to locate PSLG vertex in triangulation.");
  1064. }
  1065. }
  1066. // Remember this triangle to improve subsequent point location.
  1067. locator.Update(ref searchtri1);
  1068. // Scout the beginnings of a path from the first endpoint
  1069. // toward the second.
  1070. if (ScoutSegment(ref searchtri1, endpoint2, newmark))
  1071. {
  1072. // The segment was easily inserted.
  1073. return;
  1074. }
  1075. // The first endpoint may have changed if a collision with an intervening
  1076. // vertex on the segment occurred.
  1077. endpoint1 = searchtri1.Org();
  1078. // Find a triangle whose origin is the segment's second endpoint.
  1079. checkvertex = null;
  1080. searchtri2 = endpoint2.tri;
  1081. if (searchtri2.tri != null)
  1082. {
  1083. checkvertex = searchtri2.Org();
  1084. }
  1085. if (checkvertex != endpoint2)
  1086. {
  1087. // Find a boundary triangle to search from.
  1088. searchtri2.tri = dummytri;
  1089. searchtri2.orient = 0;
  1090. searchtri2.Sym();
  1091. // Search for the segment's second endpoint by point location.
  1092. if (locator.Locate(endpoint2, ref searchtri2) != LocateResult.OnVertex)
  1093. {
  1094. logger.Error("Unable to locate PSLG vertex in triangulation.", "Mesh.InsertSegment().2");
  1095. throw new Exception("Unable to locate PSLG vertex in triangulation.");
  1096. }
  1097. }
  1098. // Remember this triangle to improve subsequent point location.
  1099. locator.Update(ref searchtri2);
  1100. // Scout the beginnings of a path from the second endpoint
  1101. // toward the first.
  1102. if (ScoutSegment(ref searchtri2, endpoint1, newmark))
  1103. {
  1104. // The segment was easily inserted.
  1105. return;
  1106. }
  1107. // The second endpoint may have changed if a collision with an intervening
  1108. // vertex on the segment occurred.
  1109. endpoint2 = searchtri2.Org();
  1110. // Insert the segment directly into the triangulation.
  1111. ConstrainedEdge(ref searchtri1, endpoint2, newmark);
  1112. }
  1113. /// <summary>
  1114. /// Cover the convex hull of a triangulation with subsegments.
  1115. /// </summary>
  1116. private void MarkHull()
  1117. {
  1118. Otri hulltri = default(Otri);
  1119. Otri nexttri = default(Otri);
  1120. Otri starttri = default(Otri);
  1121. // Find a triangle handle on the hull.
  1122. hulltri.tri = mesh.dummytri;
  1123. hulltri.orient = 0;
  1124. hulltri.Sym();
  1125. // Remember where we started so we know when to stop.
  1126. hulltri.Copy(ref starttri);
  1127. // Go once counterclockwise around the convex hull.
  1128. do
  1129. {
  1130. // Create a subsegment if there isn't already one here.
  1131. mesh.InsertSubseg(ref hulltri, 1);
  1132. // To find the next hull edge, go clockwise around the next vertex.
  1133. hulltri.Lnext();
  1134. hulltri.Oprev(ref nexttri);
  1135. while (nexttri.tri.id != Mesh.DUMMY)
  1136. {
  1137. nexttri.Copy(ref hulltri);
  1138. hulltri.Oprev(ref nexttri);
  1139. }
  1140. }
  1141. while (!hulltri.Equals(starttri));
  1142. }
  1143. #endregion
  1144. }
  1145. }