Mesh.cs 70 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768
  1. // -----------------------------------------------------------------------
  2. // <copyright file="Mesh.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. {
  9. using System;
  10. using System.Collections.Generic;
  11. using Animation.TriangleNet.Geometry;
  12. using Animation.TriangleNet.Logging;
  13. using Animation.TriangleNet.Meshing;
  14. using Animation.TriangleNet.Meshing.Data;
  15. using Animation.TriangleNet.Meshing.Iterators;
  16. using Animation.TriangleNet.Tools;
  17. using Animation.TriangleNet.Topology;
  18. /// <summary>
  19. /// Mesh data structure.
  20. /// </summary>
  21. internal class Mesh : IMesh
  22. {
  23. #region Variables
  24. IPredicates predicates;
  25. ILog<LogItem> logger;
  26. QualityMesher qualityMesher;
  27. // Stack that maintains a list of recently flipped triangles.
  28. Stack<Otri> flipstack;
  29. // TODO: Check if custom hashmap implementation could be faster.
  30. // Using hashsets for memory management should quite fast.
  31. internal TrianglePool triangles;
  32. internal Dictionary<int, SubSegment> subsegs;
  33. internal Dictionary<int, Vertex> vertices;
  34. // Hash seeds (should belong to mesh instance)
  35. internal int hash_vtx = 0;
  36. internal int hash_seg = 0;
  37. internal int hash_tri = 0;
  38. internal List<Point> holes;
  39. internal List<RegionPointer> regions;
  40. // TODO: remove mesh_dim, invertices and insegments
  41. // Other variables.
  42. internal Rectangle bounds; // x and y bounds.
  43. internal int invertices; // Number of input vertices.
  44. internal int insegments; // Number of input segments.
  45. internal int undeads; // Number of input vertices that don't appear in the mesh.
  46. internal int mesh_dim; // Dimension (ought to be 2).
  47. internal int nextras = 0; // Number of attributes per vertex.
  48. //internal int eextras; // Number of attributes per triangle.
  49. internal int hullsize; // Number of edges in convex hull.
  50. internal int steinerleft; // Number of Steiner points not yet used.
  51. internal bool checksegments; // Are there segments in the triangulation yet?
  52. internal bool checkquality; // Has quality triangulation begun yet?
  53. // Triangular bounding box vertices.
  54. internal Vertex infvertex1, infvertex2, infvertex3;
  55. internal TriangleLocator locator;
  56. // Controls the behavior of the mesh instance.
  57. internal Behavior behavior;
  58. // The current node numbering
  59. internal NodeNumbering numbering;
  60. #endregion
  61. #region Public properties
  62. /// <summary>
  63. /// Gets the mesh bounding box.
  64. /// </summary>
  65. public Rectangle Bounds
  66. {
  67. get { return this.bounds; }
  68. }
  69. /// <summary>
  70. /// Gets the mesh vertices.
  71. /// </summary>
  72. public ICollection<Vertex> Vertices
  73. {
  74. get { return this.vertices.Values; }
  75. }
  76. /// <summary>
  77. /// Gets the mesh holes.
  78. /// </summary>
  79. public IList<Point> Holes
  80. {
  81. get { return this.holes; }
  82. }
  83. /// <summary>
  84. /// Gets the mesh triangles.
  85. /// </summary>
  86. public ICollection<Triangle> Triangles
  87. {
  88. get { return this.triangles; }
  89. }
  90. /// <summary>
  91. /// Gets the mesh segments.
  92. /// </summary>
  93. public ICollection<SubSegment> Segments
  94. {
  95. get { return this.subsegs.Values; }
  96. }
  97. /// <summary>
  98. /// Gets the mesh edges.
  99. /// </summary>
  100. public IEnumerable<Edge> Edges
  101. {
  102. get
  103. {
  104. var e = new EdgeIterator(this);
  105. while (e.MoveNext())
  106. {
  107. yield return e.Current;
  108. }
  109. }
  110. }
  111. /// <summary>
  112. /// Gets the number of input vertices.
  113. /// </summary>
  114. public int NumberOfInputPoints
  115. {
  116. get { return invertices; }
  117. }
  118. /// <summary>
  119. /// Gets the number of mesh edges.
  120. /// </summary>
  121. public int NumberOfEdges
  122. {
  123. get { return (3 * triangles.Count + hullsize) / 2; }
  124. }
  125. /// <summary>
  126. /// Indicates whether the input is a PSLG or a point set.
  127. /// </summary>
  128. public bool IsPolygon
  129. {
  130. get { return this.insegments > 0; }
  131. }
  132. /// <summary>
  133. /// Gets the current node numbering.
  134. /// </summary>
  135. public NodeNumbering CurrentNumbering
  136. {
  137. get { return numbering; }
  138. }
  139. #endregion
  140. #region "Outer space" variables
  141. internal const int DUMMY = -1;
  142. // The triangle that fills "outer space," called 'dummytri', is pointed to
  143. // by every triangle and subsegment on a boundary (be it outer or inner) of
  144. // the triangulation. Also, 'dummytri' points to one of the triangles on
  145. // the convex hull (until the holes and concavities are carved), making it
  146. // possible to find a starting triangle for point location.
  147. // 'dummytri' and 'dummysub' are generally required to fulfill only a few
  148. // invariants: their vertices must remain NULL and 'dummytri' must always
  149. // be bonded (at offset zero) to some triangle on the convex hull of the
  150. // mesh, via a boundary edge. Otherwise, the connections of 'dummytri' and
  151. // 'dummysub' may change willy-nilly. This makes it possible to avoid
  152. // writing a good deal of special-case code (in the edge flip, for example)
  153. // for dealing with the boundary of the mesh, places where no subsegment is
  154. // present, and so forth. Other entities are frequently bonded to
  155. // 'dummytri' and 'dummysub' as if they were real mesh entities, with no
  156. // harm done.
  157. internal Triangle dummytri;
  158. // Set up 'dummysub', the omnipresent subsegment pointed to by any
  159. // triangle side or subsegment end that isn't attached to a real
  160. // subsegment.
  161. internal SubSegment dummysub;
  162. private void Initialize()
  163. {
  164. dummysub = new SubSegment();
  165. dummysub.hash = DUMMY;
  166. // Initialize the two adjoining subsegments to be the omnipresent
  167. // subsegment. These will eventually be changed by various bonding
  168. // operations, but their values don't really matter, as long as they
  169. // can legally be dereferenced.
  170. dummysub.subsegs[0].seg = dummysub;
  171. dummysub.subsegs[1].seg = dummysub;
  172. // Set up 'dummytri', the 'triangle' that occupies "outer space."
  173. dummytri = new Triangle();
  174. dummytri.hash = dummytri.id = DUMMY;
  175. // Initialize the three adjoining triangles to be "outer space." These
  176. // will eventually be changed by various bonding operations, but their
  177. // values don't really matter, as long as they can legally be
  178. // dereferenced.
  179. dummytri.neighbors[0].tri = dummytri;
  180. dummytri.neighbors[1].tri = dummytri;
  181. dummytri.neighbors[2].tri = dummytri;
  182. // Initialize the three adjoining subsegments of 'dummytri' to be
  183. // the omnipresent subsegment.
  184. dummytri.subsegs[0].seg = dummysub;
  185. dummytri.subsegs[1].seg = dummysub;
  186. dummytri.subsegs[2].seg = dummysub;
  187. }
  188. #endregion
  189. /// <summary>
  190. /// Initializes a new instance of the <see cref="Mesh" /> class.
  191. /// </summary>
  192. public Mesh(Configuration config)
  193. {
  194. Initialize();
  195. logger = Log.Instance;
  196. behavior = new Behavior();
  197. vertices = new Dictionary<int, Vertex>();
  198. subsegs = new Dictionary<int, SubSegment>();
  199. triangles = config.TrianglePool();
  200. flipstack = new Stack<Otri>();
  201. holes = new List<Point>();
  202. regions = new List<RegionPointer>();
  203. steinerleft = -1;
  204. this.predicates = config.Predicates();
  205. this.locator = new TriangleLocator(this, predicates);
  206. }
  207. public void Refine(QualityOptions quality, bool delaunay = false)
  208. {
  209. invertices = vertices.Count;
  210. if (behavior.Poly)
  211. {
  212. insegments = behavior.useSegments ? subsegs.Count : hullsize;
  213. }
  214. Reset();
  215. if (qualityMesher == null)
  216. {
  217. qualityMesher = new QualityMesher(this, new Configuration());
  218. }
  219. // Enforce angle and area constraints.
  220. qualityMesher.Apply(quality, delaunay);
  221. }
  222. /// <summary>
  223. /// Renumber vertex and triangle id's.
  224. /// </summary>
  225. public void Renumber()
  226. {
  227. this.Renumber(NodeNumbering.Linear);
  228. }
  229. /// <summary>
  230. /// Renumber vertex and triangle id's.
  231. /// </summary>
  232. public void Renumber(NodeNumbering num)
  233. {
  234. // Don't need to do anything if the nodes are already numbered.
  235. if (num == this.numbering)
  236. {
  237. return;
  238. }
  239. int id;
  240. if (num == NodeNumbering.Linear)
  241. {
  242. id = 0;
  243. foreach (var node in this.vertices.Values)
  244. {
  245. node.id = id++;
  246. }
  247. }
  248. else if (num == NodeNumbering.CuthillMcKee)
  249. {
  250. var rcm = new CuthillMcKee();
  251. var iperm = rcm.Renumber(this);
  252. // Permute the node indices.
  253. foreach (var node in this.vertices.Values)
  254. {
  255. node.id = iperm[node.id];
  256. }
  257. }
  258. // Remember the current numbering.
  259. numbering = num;
  260. // Triangles will always be numbered from 0 to n-1
  261. id = 0;
  262. foreach (var item in this.triangles)
  263. {
  264. item.id = id++;
  265. }
  266. }
  267. #region Misc
  268. /// <summary>
  269. /// Set QualityMesher for mesh refinement.
  270. /// </summary>
  271. /// <param name="qmesher"></param>
  272. internal void SetQualityMesher(QualityMesher qmesher)
  273. {
  274. qualityMesher = qmesher;
  275. }
  276. internal void CopyTo(Mesh target)
  277. {
  278. target.vertices = this.vertices;
  279. target.triangles = this.triangles;
  280. target.subsegs = this.subsegs;
  281. target.holes = this.holes;
  282. target.regions = this.regions;
  283. target.hash_vtx = this.hash_vtx;
  284. target.hash_seg = this.hash_seg;
  285. target.hash_tri = this.hash_tri;
  286. target.numbering = this.numbering;
  287. target.hullsize = this.hullsize;
  288. }
  289. /// <summary>
  290. /// Reset all the mesh data. This method will also wipe
  291. /// out all mesh data.
  292. /// </summary>
  293. private void ResetData()
  294. {
  295. vertices.Clear();
  296. triangles.Restart();
  297. subsegs.Clear();
  298. holes.Clear();
  299. regions.Clear();
  300. this.hash_vtx = 0;
  301. this.hash_seg = 0;
  302. this.hash_tri = 0;
  303. flipstack.Clear();
  304. hullsize = 0;
  305. Reset();
  306. locator.Reset();
  307. }
  308. /// <summary>
  309. /// Reset the mesh triangulation state.
  310. /// </summary>
  311. private void Reset()
  312. {
  313. numbering = NodeNumbering.None;
  314. undeads = 0; // No eliminated input vertices yet.
  315. checksegments = false; // There are no segments in the triangulation yet.
  316. checkquality = false; // The quality triangulation stage has not begun.
  317. Statistic.InCircleCount = 0;
  318. Statistic.CounterClockwiseCount = 0;
  319. Statistic.InCircleAdaptCount = 0;
  320. Statistic.CounterClockwiseAdaptCount = 0;
  321. Statistic.Orient3dCount = 0;
  322. Statistic.HyperbolaCount = 0;
  323. Statistic.CircleTopCount = 0;
  324. Statistic.CircumcenterCount = 0;
  325. }
  326. /// <summary>
  327. /// Read the vertices from memory.
  328. /// </summary>
  329. /// <param name="data">The input data.</param>
  330. internal void TransferNodes(IList<Vertex> points)
  331. {
  332. this.invertices = points.Count;
  333. this.mesh_dim = 2;
  334. this.bounds = new Rectangle();
  335. if (this.invertices < 3)
  336. {
  337. logger.Error("Input must have at least three input vertices.", "Mesh.TransferNodes()");
  338. throw new Exception("Input must have at least three input vertices.");
  339. }
  340. var v = points[0];
  341. #if USE_ATTRIBS
  342. // Check attributes.
  343. this.nextras = v.attributes == null ? 0 : v.attributes.Length;
  344. #endif
  345. // Simple heuristic to check if ids are already set. We assume that if the
  346. // first two vertex ids are distinct, then all input vertices have pairwise
  347. // distinct ids.
  348. bool userId = (v.id != points[1].id);
  349. foreach (var p in points)
  350. {
  351. if (userId)
  352. {
  353. p.hash = p.id;
  354. // Make sure the hash counter gets updated.
  355. hash_vtx = Math.Max(p.hash + 1, hash_vtx);
  356. }
  357. else
  358. {
  359. p.hash = p.id = hash_vtx++;
  360. }
  361. this.vertices.Add(p.hash, p);
  362. this.bounds.Expand(p);
  363. }
  364. }
  365. /// <summary>
  366. /// Construct a mapping from vertices to triangles to improve the speed of
  367. /// point location for segment insertion.
  368. /// </summary>
  369. /// <remarks>
  370. /// Traverses all the triangles, and provides each corner of each triangle
  371. /// with a pointer to that triangle. Of course, pointers will be overwritten
  372. /// by other pointers because (almost) each vertex is a corner of several
  373. /// triangles, but in the end every vertex will point to some triangle
  374. /// that contains it.
  375. /// </remarks>
  376. internal void MakeVertexMap()
  377. {
  378. Otri tri = default(Otri);
  379. Vertex triorg;
  380. foreach (var t in this.triangles)
  381. {
  382. tri.tri = t;
  383. // Check all three vertices of the triangle.
  384. for (tri.orient = 0; tri.orient < 3; tri.orient++)
  385. {
  386. triorg = tri.Org();
  387. triorg.tri = tri;
  388. }
  389. }
  390. }
  391. #endregion
  392. #region Factory
  393. /// <summary>
  394. /// Create a new triangle with orientation zero.
  395. /// </summary>
  396. /// <param name="newotri">Reference to the new triangle.</param>
  397. internal void MakeTriangle(ref Otri newotri)
  398. {
  399. Triangle tri = triangles.Get();
  400. //tri.id = tri.hash;
  401. tri.subsegs[0].seg = dummysub;
  402. tri.subsegs[1].seg = dummysub;
  403. tri.subsegs[2].seg = dummysub;
  404. tri.neighbors[0].tri = dummytri;
  405. tri.neighbors[1].tri = dummytri;
  406. tri.neighbors[2].tri = dummytri;
  407. newotri.tri = tri;
  408. newotri.orient = 0;
  409. }
  410. /// <summary>
  411. /// Create a new subsegment with orientation zero.
  412. /// </summary>
  413. /// <param name="newsubseg">Reference to the new subseg.</param>
  414. internal void MakeSegment(ref Osub newsubseg)
  415. {
  416. var seg = new SubSegment();
  417. seg.hash = this.hash_seg++;
  418. seg.subsegs[0].seg = dummysub;
  419. seg.subsegs[1].seg = dummysub;
  420. seg.triangles[0].tri = dummytri;
  421. seg.triangles[1].tri = dummytri;
  422. newsubseg.seg = seg;
  423. newsubseg.orient = 0;
  424. subsegs.Add(seg.hash, seg);
  425. }
  426. #endregion
  427. #region Manipulation
  428. /// <summary>
  429. /// Insert a vertex into a Delaunay triangulation, performing flips as necessary
  430. /// to maintain the Delaunay property.
  431. /// </summary>
  432. /// <param name="newvertex">The point to be inserted.</param>
  433. /// <param name="searchtri">The triangle to start the search.</param>
  434. /// <param name="splitseg">Segment to split.</param>
  435. /// <param name="segmentflaws">Check for creation of encroached subsegments.</param>
  436. /// <param name="triflaws">Check for creation of bad quality triangles.</param>
  437. /// <returns>If a duplicate vertex or violated segment does not prevent the
  438. /// vertex from being inserted, the return value will be ENCROACHINGVERTEX if
  439. /// the vertex encroaches upon a subsegment (and checking is enabled), or
  440. /// SUCCESSFULVERTEX otherwise. In either case, 'searchtri' is set to a handle
  441. /// whose origin is the newly inserted vertex.</returns>
  442. /// <remarks>
  443. /// The point 'newvertex' is located. If 'searchtri.triangle' is not NULL,
  444. /// the search for the containing triangle begins from 'searchtri'. If
  445. /// 'searchtri.triangle' is NULL, a full point location procedure is called.
  446. /// If 'insertvertex' is found inside a triangle, the triangle is split into
  447. /// three; if 'insertvertex' lies on an edge, the edge is split in two,
  448. /// thereby splitting the two adjacent triangles into four. Edge flips are
  449. /// used to restore the Delaunay property. If 'insertvertex' lies on an
  450. /// existing vertex, no action is taken, and the value DUPLICATEVERTEX is
  451. /// returned. On return, 'searchtri' is set to a handle whose origin is the
  452. /// existing vertex.
  453. ///
  454. /// InsertVertex() does not use flip() for reasons of speed; some
  455. /// information can be reused from edge flip to edge flip, like the
  456. /// locations of subsegments.
  457. ///
  458. /// Param 'splitseg': Normally, the parameter 'splitseg' is set to NULL,
  459. /// implying that no subsegment should be split. In this case, if 'insertvertex'
  460. /// is found to lie on a segment, no action is taken, and the value VIOLATINGVERTEX
  461. /// is returned. On return, 'searchtri' is set to a handle whose primary edge is the
  462. /// violated subsegment.
  463. /// If the calling routine wishes to split a subsegment by inserting a vertex in it,
  464. /// the parameter 'splitseg' should be that subsegment. In this case, 'searchtri'
  465. /// MUST be the triangle handle reached by pivoting from that subsegment; no point
  466. /// location is done.
  467. ///
  468. /// Param 'segmentflaws': Flags that indicate whether or not there should
  469. /// be checks for the creation of encroached subsegments. If a newly inserted
  470. /// vertex encroaches upon subsegments, these subsegments are added to the list
  471. /// of subsegments to be split if 'segmentflaws' is set.
  472. ///
  473. /// Param 'triflaws': Flags that indicate whether or not there should be
  474. /// checks for the creation of bad quality triangles. If bad triangles are
  475. /// created, these are added to the queue if 'triflaws' is set.
  476. /// </remarks>
  477. internal InsertVertexResult InsertVertex(Vertex newvertex, ref Otri searchtri,
  478. ref Osub splitseg, bool segmentflaws, bool triflaws)
  479. {
  480. Otri horiz = default(Otri);
  481. Otri top = default(Otri);
  482. Otri botleft = default(Otri), botright = default(Otri);
  483. Otri topleft = default(Otri), topright = default(Otri);
  484. Otri newbotleft = default(Otri), newbotright = default(Otri);
  485. Otri newtopright = default(Otri);
  486. Otri botlcasing = default(Otri), botrcasing = default(Otri);
  487. Otri toplcasing = default(Otri), toprcasing = default(Otri);
  488. Otri testtri = default(Otri);
  489. Osub botlsubseg = default(Osub), botrsubseg = default(Osub);
  490. Osub toplsubseg = default(Osub), toprsubseg = default(Osub);
  491. Osub brokensubseg = default(Osub);
  492. Osub checksubseg = default(Osub);
  493. Osub rightsubseg = default(Osub);
  494. Osub newsubseg = default(Osub);
  495. BadSubseg encroached;
  496. //FlipStacker newflip;
  497. Vertex first;
  498. Vertex leftvertex, rightvertex, botvertex, topvertex, farvertex;
  499. Vertex segmentorg, segmentdest;
  500. int region;
  501. double area;
  502. InsertVertexResult success;
  503. LocateResult intersect;
  504. bool doflip;
  505. bool mirrorflag;
  506. bool enq;
  507. if (splitseg.seg == null)
  508. {
  509. // Find the location of the vertex to be inserted. Check if a good
  510. // starting triangle has already been provided by the caller.
  511. if (searchtri.tri.id == DUMMY)
  512. {
  513. // Find a boundary triangle.
  514. horiz.tri = dummytri;
  515. horiz.orient = 0;
  516. horiz.Sym();
  517. // Search for a triangle containing 'newvertex'.
  518. intersect = locator.Locate(newvertex, ref horiz);
  519. }
  520. else
  521. {
  522. // Start searching from the triangle provided by the caller.
  523. searchtri.Copy(ref horiz);
  524. intersect = locator.PreciseLocate(newvertex, ref horiz, true);
  525. }
  526. }
  527. else
  528. {
  529. // The calling routine provides the subsegment in which
  530. // the vertex is inserted.
  531. searchtri.Copy(ref horiz);
  532. intersect = LocateResult.OnEdge;
  533. }
  534. if (intersect == LocateResult.OnVertex)
  535. {
  536. // There's already a vertex there. Return in 'searchtri' a triangle
  537. // whose origin is the existing vertex.
  538. horiz.Copy(ref searchtri);
  539. locator.Update(ref horiz);
  540. return InsertVertexResult.Duplicate;
  541. }
  542. if ((intersect == LocateResult.OnEdge) || (intersect == LocateResult.Outside))
  543. {
  544. // The vertex falls on an edge or boundary.
  545. if (checksegments && (splitseg.seg == null))
  546. {
  547. // Check whether the vertex falls on a subsegment.
  548. horiz.Pivot(ref brokensubseg);
  549. if (brokensubseg.seg.hash != DUMMY)
  550. {
  551. // The vertex falls on a subsegment, and hence will not be inserted.
  552. if (segmentflaws)
  553. {
  554. enq = behavior.NoBisect != 2;
  555. if (enq && (behavior.NoBisect == 1))
  556. {
  557. // This subsegment may be split only if it is an
  558. // internal boundary.
  559. horiz.Sym(ref testtri);
  560. enq = testtri.tri.id != DUMMY;
  561. }
  562. if (enq)
  563. {
  564. // Add the subsegment to the list of encroached subsegments.
  565. encroached = new BadSubseg();
  566. encroached.subseg = brokensubseg;
  567. encroached.org = brokensubseg.Org();
  568. encroached.dest = brokensubseg.Dest();
  569. qualityMesher.AddBadSubseg(encroached);
  570. }
  571. }
  572. // Return a handle whose primary edge contains the vertex,
  573. // which has not been inserted.
  574. horiz.Copy(ref searchtri);
  575. locator.Update(ref horiz);
  576. return InsertVertexResult.Violating;
  577. }
  578. }
  579. // Insert the vertex on an edge, dividing one triangle into two (if
  580. // the edge lies on a boundary) or two triangles into four.
  581. horiz.Lprev(ref botright);
  582. botright.Sym(ref botrcasing);
  583. horiz.Sym(ref topright);
  584. // Is there a second triangle? (Or does this edge lie on a boundary?)
  585. mirrorflag = topright.tri.id != DUMMY;
  586. if (mirrorflag)
  587. {
  588. topright.Lnext();
  589. topright.Sym(ref toprcasing);
  590. MakeTriangle(ref newtopright);
  591. }
  592. else
  593. {
  594. // Splitting a boundary edge increases the number of boundary edges.
  595. hullsize++;
  596. }
  597. MakeTriangle(ref newbotright);
  598. // Set the vertices of changed and new triangles.
  599. rightvertex = horiz.Org();
  600. leftvertex = horiz.Dest();
  601. botvertex = horiz.Apex();
  602. newbotright.SetOrg(botvertex);
  603. newbotright.SetDest(rightvertex);
  604. newbotright.SetApex(newvertex);
  605. horiz.SetOrg(newvertex);
  606. // Set the region of a new triangle.
  607. newbotright.tri.label = botright.tri.label;
  608. if (behavior.VarArea)
  609. {
  610. // Set the area constraint of a new triangle.
  611. newbotright.tri.area = botright.tri.area;
  612. }
  613. if (mirrorflag)
  614. {
  615. topvertex = topright.Dest();
  616. newtopright.SetOrg(rightvertex);
  617. newtopright.SetDest(topvertex);
  618. newtopright.SetApex(newvertex);
  619. topright.SetOrg(newvertex);
  620. // Set the region of another new triangle.
  621. newtopright.tri.label = topright.tri.label;
  622. if (behavior.VarArea)
  623. {
  624. // Set the area constraint of another new triangle.
  625. newtopright.tri.area = topright.tri.area;
  626. }
  627. }
  628. // There may be subsegments that need to be bonded
  629. // to the new triangle(s).
  630. if (checksegments)
  631. {
  632. botright.Pivot(ref botrsubseg);
  633. if (botrsubseg.seg.hash != DUMMY)
  634. {
  635. botright.SegDissolve(dummysub);
  636. newbotright.SegBond(ref botrsubseg);
  637. }
  638. if (mirrorflag)
  639. {
  640. topright.Pivot(ref toprsubseg);
  641. if (toprsubseg.seg.hash != DUMMY)
  642. {
  643. topright.SegDissolve(dummysub);
  644. newtopright.SegBond(ref toprsubseg);
  645. }
  646. }
  647. }
  648. // Bond the new triangle(s) to the surrounding triangles.
  649. newbotright.Bond(ref botrcasing);
  650. newbotright.Lprev();
  651. newbotright.Bond(ref botright);
  652. newbotright.Lprev();
  653. if (mirrorflag)
  654. {
  655. newtopright.Bond(ref toprcasing);
  656. newtopright.Lnext();
  657. newtopright.Bond(ref topright);
  658. newtopright.Lnext();
  659. newtopright.Bond(ref newbotright);
  660. }
  661. if (splitseg.seg != null)
  662. {
  663. // Split the subsegment into two.
  664. splitseg.SetDest(newvertex);
  665. segmentorg = splitseg.SegOrg();
  666. segmentdest = splitseg.SegDest();
  667. splitseg.Sym();
  668. splitseg.Pivot(ref rightsubseg);
  669. InsertSubseg(ref newbotright, splitseg.seg.boundary);
  670. newbotright.Pivot(ref newsubseg);
  671. newsubseg.SetSegOrg(segmentorg);
  672. newsubseg.SetSegDest(segmentdest);
  673. splitseg.Bond(ref newsubseg);
  674. newsubseg.Sym();
  675. newsubseg.Bond(ref rightsubseg);
  676. splitseg.Sym();
  677. // Transfer the subsegment's boundary marker to the vertex if required.
  678. if (newvertex.label == 0)
  679. {
  680. newvertex.label = splitseg.seg.boundary;
  681. }
  682. }
  683. if (checkquality)
  684. {
  685. flipstack.Clear();
  686. flipstack.Push(default(Otri)); // Dummy flip (see UndoVertex)
  687. flipstack.Push(horiz);
  688. }
  689. // Position 'horiz' on the first edge to check for
  690. // the Delaunay property.
  691. horiz.Lnext();
  692. }
  693. else
  694. {
  695. // Insert the vertex in a triangle, splitting it into three.
  696. horiz.Lnext(ref botleft);
  697. horiz.Lprev(ref botright);
  698. botleft.Sym(ref botlcasing);
  699. botright.Sym(ref botrcasing);
  700. MakeTriangle(ref newbotleft);
  701. MakeTriangle(ref newbotright);
  702. // Set the vertices of changed and new triangles.
  703. rightvertex = horiz.Org();
  704. leftvertex = horiz.Dest();
  705. botvertex = horiz.Apex();
  706. newbotleft.SetOrg(leftvertex);
  707. newbotleft.SetDest(botvertex);
  708. newbotleft.SetApex(newvertex);
  709. newbotright.SetOrg(botvertex);
  710. newbotright.SetDest(rightvertex);
  711. newbotright.SetApex(newvertex);
  712. horiz.SetApex(newvertex);
  713. // Set the region of the new triangles.
  714. newbotleft.tri.label = horiz.tri.label;
  715. newbotright.tri.label = horiz.tri.label;
  716. if (behavior.VarArea)
  717. {
  718. // Set the area constraint of the new triangles.
  719. area = horiz.tri.area;
  720. newbotleft.tri.area = area;
  721. newbotright.tri.area = area;
  722. }
  723. // There may be subsegments that need to be bonded
  724. // to the new triangles.
  725. if (checksegments)
  726. {
  727. botleft.Pivot(ref botlsubseg);
  728. if (botlsubseg.seg.hash != DUMMY)
  729. {
  730. botleft.SegDissolve(dummysub);
  731. newbotleft.SegBond(ref botlsubseg);
  732. }
  733. botright.Pivot(ref botrsubseg);
  734. if (botrsubseg.seg.hash != DUMMY)
  735. {
  736. botright.SegDissolve(dummysub);
  737. newbotright.SegBond(ref botrsubseg);
  738. }
  739. }
  740. // Bond the new triangles to the surrounding triangles.
  741. newbotleft.Bond(ref botlcasing);
  742. newbotright.Bond(ref botrcasing);
  743. newbotleft.Lnext();
  744. newbotright.Lprev();
  745. newbotleft.Bond(ref newbotright);
  746. newbotleft.Lnext();
  747. botleft.Bond(ref newbotleft);
  748. newbotright.Lprev();
  749. botright.Bond(ref newbotright);
  750. if (checkquality)
  751. {
  752. flipstack.Clear();
  753. flipstack.Push(horiz);
  754. }
  755. }
  756. // The insertion is successful by default, unless an encroached
  757. // subsegment is found.
  758. success = InsertVertexResult.Successful;
  759. if (newvertex.tri.tri != null)
  760. {
  761. // Store the coordinates of the triangle that contains newvertex.
  762. newvertex.tri.SetOrg(rightvertex);
  763. newvertex.tri.SetDest(leftvertex);
  764. newvertex.tri.SetApex(botvertex);
  765. }
  766. // Circle around the newly inserted vertex, checking each edge opposite it
  767. // for the Delaunay property. Non-Delaunay edges are flipped. 'horiz' is
  768. // always the edge being checked. 'first' marks where to stop circling.
  769. first = horiz.Org();
  770. rightvertex = first;
  771. leftvertex = horiz.Dest();
  772. // Circle until finished.
  773. while (true)
  774. {
  775. // By default, the edge will be flipped.
  776. doflip = true;
  777. if (checksegments)
  778. {
  779. // Check for a subsegment, which cannot be flipped.
  780. horiz.Pivot(ref checksubseg);
  781. if (checksubseg.seg.hash != DUMMY)
  782. {
  783. // The edge is a subsegment and cannot be flipped.
  784. doflip = false;
  785. if (segmentflaws)
  786. {
  787. // Does the new vertex encroach upon this subsegment?
  788. if (qualityMesher.CheckSeg4Encroach(ref checksubseg) > 0)
  789. {
  790. success = InsertVertexResult.Encroaching;
  791. }
  792. }
  793. }
  794. }
  795. if (doflip)
  796. {
  797. // Check if the edge is a boundary edge.
  798. horiz.Sym(ref top);
  799. if (top.tri.id == DUMMY)
  800. {
  801. // The edge is a boundary edge and cannot be flipped.
  802. doflip = false;
  803. }
  804. else
  805. {
  806. // Find the vertex on the other side of the edge.
  807. farvertex = top.Apex();
  808. // In the incremental Delaunay triangulation algorithm, any of
  809. // 'leftvertex', 'rightvertex', and 'farvertex' could be vertices
  810. // of the triangular bounding box. These vertices must be
  811. // treated as if they are infinitely distant, even though their
  812. // "coordinates" are not.
  813. if ((leftvertex == infvertex1) || (leftvertex == infvertex2) ||
  814. (leftvertex == infvertex3))
  815. {
  816. // 'leftvertex' is infinitely distant. Check the convexity of
  817. // the boundary of the triangulation. 'farvertex' might be
  818. // infinite as well, but trust me, this same condition should
  819. // be applied.
  820. doflip = predicates.CounterClockwise(newvertex, rightvertex, farvertex) > 0.0;
  821. }
  822. else if ((rightvertex == infvertex1) ||
  823. (rightvertex == infvertex2) ||
  824. (rightvertex == infvertex3))
  825. {
  826. // 'rightvertex' is infinitely distant. Check the convexity of
  827. // the boundary of the triangulation. 'farvertex' might be
  828. // infinite as well, but trust me, this same condition should
  829. // be applied.
  830. doflip = predicates.CounterClockwise(farvertex, leftvertex, newvertex) > 0.0;
  831. }
  832. else if ((farvertex == infvertex1) ||
  833. (farvertex == infvertex2) ||
  834. (farvertex == infvertex3))
  835. {
  836. // 'farvertex' is infinitely distant and cannot be inside
  837. // the circumcircle of the triangle 'horiz'.
  838. doflip = false;
  839. }
  840. else
  841. {
  842. // Test whether the edge is locally Delaunay.
  843. doflip = predicates.InCircle(leftvertex, newvertex, rightvertex, farvertex) > 0.0;
  844. }
  845. if (doflip)
  846. {
  847. // We made it! Flip the edge 'horiz' by rotating its containing
  848. // quadrilateral (the two triangles adjacent to 'horiz').
  849. // Identify the casing of the quadrilateral.
  850. top.Lprev(ref topleft);
  851. topleft.Sym(ref toplcasing);
  852. top.Lnext(ref topright);
  853. topright.Sym(ref toprcasing);
  854. horiz.Lnext(ref botleft);
  855. botleft.Sym(ref botlcasing);
  856. horiz.Lprev(ref botright);
  857. botright.Sym(ref botrcasing);
  858. // Rotate the quadrilateral one-quarter turn counterclockwise.
  859. topleft.Bond(ref botlcasing);
  860. botleft.Bond(ref botrcasing);
  861. botright.Bond(ref toprcasing);
  862. topright.Bond(ref toplcasing);
  863. if (checksegments)
  864. {
  865. // Check for subsegments and rebond them to the quadrilateral.
  866. topleft.Pivot(ref toplsubseg);
  867. botleft.Pivot(ref botlsubseg);
  868. botright.Pivot(ref botrsubseg);
  869. topright.Pivot(ref toprsubseg);
  870. if (toplsubseg.seg.hash == DUMMY)
  871. {
  872. topright.SegDissolve(dummysub);
  873. }
  874. else
  875. {
  876. topright.SegBond(ref toplsubseg);
  877. }
  878. if (botlsubseg.seg.hash == DUMMY)
  879. {
  880. topleft.SegDissolve(dummysub);
  881. }
  882. else
  883. {
  884. topleft.SegBond(ref botlsubseg);
  885. }
  886. if (botrsubseg.seg.hash == DUMMY)
  887. {
  888. botleft.SegDissolve(dummysub);
  889. }
  890. else
  891. {
  892. botleft.SegBond(ref botrsubseg);
  893. }
  894. if (toprsubseg.seg.hash == DUMMY)
  895. {
  896. botright.SegDissolve(dummysub);
  897. }
  898. else
  899. {
  900. botright.SegBond(ref toprsubseg);
  901. }
  902. }
  903. // New vertex assignments for the rotated quadrilateral.
  904. horiz.SetOrg(farvertex);
  905. horiz.SetDest(newvertex);
  906. horiz.SetApex(rightvertex);
  907. top.SetOrg(newvertex);
  908. top.SetDest(farvertex);
  909. top.SetApex(leftvertex);
  910. // Assign region.
  911. // TODO: check region ok (no Math.Min necessary)
  912. region = Math.Min(top.tri.label, horiz.tri.label);
  913. top.tri.label = region;
  914. horiz.tri.label = region;
  915. if (behavior.VarArea)
  916. {
  917. if ((top.tri.area <= 0.0) || (horiz.tri.area <= 0.0))
  918. {
  919. area = -1.0;
  920. }
  921. else
  922. {
  923. // Take the average of the two triangles' area constraints.
  924. // This prevents small area constraints from migrating a
  925. // long, long way from their original location due to flips.
  926. area = 0.5 * (top.tri.area + horiz.tri.area);
  927. }
  928. top.tri.area = area;
  929. horiz.tri.area = area;
  930. }
  931. if (checkquality)
  932. {
  933. flipstack.Push(horiz);
  934. }
  935. // On the next iterations, consider the two edges that were exposed (this
  936. // is, are now visible to the newly inserted vertex) by the edge flip.
  937. horiz.Lprev();
  938. leftvertex = farvertex;
  939. }
  940. }
  941. }
  942. if (!doflip)
  943. {
  944. // The handle 'horiz' is accepted as locally Delaunay.
  945. if (triflaws)
  946. {
  947. // Check the triangle 'horiz' for quality.
  948. qualityMesher.TestTriangle(ref horiz);
  949. }
  950. // Look for the next edge around the newly inserted vertex.
  951. horiz.Lnext();
  952. horiz.Sym(ref testtri);
  953. // Check for finishing a complete revolution about the new vertex, or
  954. // falling outside of the triangulation. The latter will happen when
  955. // a vertex is inserted at a boundary.
  956. if ((leftvertex == first) || (testtri.tri.id == DUMMY))
  957. {
  958. // We're done. Return a triangle whose origin is the new vertex.
  959. horiz.Lnext(ref searchtri);
  960. Otri recenttri = default(Otri);
  961. horiz.Lnext(ref recenttri);
  962. locator.Update(ref recenttri);
  963. return success;
  964. }
  965. // Finish finding the next edge around the newly inserted vertex.
  966. testtri.Lnext(ref horiz);
  967. rightvertex = leftvertex;
  968. leftvertex = horiz.Dest();
  969. }
  970. }
  971. }
  972. /// <summary>
  973. /// Create a new subsegment and inserts it between two triangles. Its
  974. /// vertices are properly initialized.
  975. /// </summary>
  976. /// <param name="tri">The new subsegment is inserted at the edge
  977. /// described by this handle.</param>
  978. /// <param name="subsegmark">The marker 'subsegmark' is applied to the
  979. /// subsegment and, if appropriate, its vertices.</param>
  980. internal void InsertSubseg(ref Otri tri, int subsegmark)
  981. {
  982. Otri oppotri = default(Otri);
  983. Osub newsubseg = default(Osub);
  984. Vertex triorg, tridest;
  985. triorg = tri.Org();
  986. tridest = tri.Dest();
  987. // Mark vertices if possible.
  988. if (triorg.label == 0)
  989. {
  990. triorg.label = subsegmark;
  991. }
  992. if (tridest.label == 0)
  993. {
  994. tridest.label = subsegmark;
  995. }
  996. // Check if there's already a subsegment here.
  997. tri.Pivot(ref newsubseg);
  998. if (newsubseg.seg.hash == DUMMY)
  999. {
  1000. // Make new subsegment and initialize its vertices.
  1001. MakeSegment(ref newsubseg);
  1002. newsubseg.SetOrg(tridest);
  1003. newsubseg.SetDest(triorg);
  1004. newsubseg.SetSegOrg(tridest);
  1005. newsubseg.SetSegDest(triorg);
  1006. // Bond new subsegment to the two triangles it is sandwiched between.
  1007. // Note that the facing triangle 'oppotri' might be equal to 'dummytri'
  1008. // (outer space), but the new subsegment is bonded to it all the same.
  1009. tri.SegBond(ref newsubseg);
  1010. tri.Sym(ref oppotri);
  1011. newsubseg.Sym();
  1012. oppotri.SegBond(ref newsubseg);
  1013. newsubseg.seg.boundary = subsegmark;
  1014. }
  1015. else if (newsubseg.seg.boundary == 0)
  1016. {
  1017. newsubseg.seg.boundary = subsegmark;
  1018. }
  1019. }
  1020. /// <summary>
  1021. /// Transform two triangles to two different triangles by flipping an edge
  1022. /// counterclockwise within a quadrilateral.
  1023. /// </summary>
  1024. /// <param name="flipedge">Handle to the edge that will be flipped.</param>
  1025. /// <remarks>Imagine the original triangles, abc and bad, oriented so that the
  1026. /// shared edge ab lies in a horizontal plane, with the vertex b on the left
  1027. /// and the vertex a on the right. The vertex c lies below the edge, and
  1028. /// the vertex d lies above the edge. The 'flipedge' handle holds the edge
  1029. /// ab of triangle abc, and is directed left, from vertex a to vertex b.
  1030. ///
  1031. /// The triangles abc and bad are deleted and replaced by the triangles cdb
  1032. /// and dca. The triangles that represent abc and bad are NOT deallocated;
  1033. /// they are reused for dca and cdb, respectively. Hence, any handles that
  1034. /// may have held the original triangles are still valid, although not
  1035. /// directed as they were before.
  1036. ///
  1037. /// Upon completion of this routine, the 'flipedge' handle holds the edge
  1038. /// dc of triangle dca, and is directed down, from vertex d to vertex c.
  1039. /// (Hence, the two triangles have rotated counterclockwise.)
  1040. ///
  1041. /// WARNING: This transformation is geometrically valid only if the
  1042. /// quadrilateral adbc is convex. Furthermore, this transformation is
  1043. /// valid only if there is not a subsegment between the triangles abc and
  1044. /// bad. This routine does not check either of these preconditions, and
  1045. /// it is the responsibility of the calling routine to ensure that they are
  1046. /// met. If they are not, the streets shall be filled with wailing and
  1047. /// gnashing of teeth.
  1048. ///
  1049. /// Terminology
  1050. ///
  1051. /// A "local transformation" replaces a small set of triangles with another
  1052. /// set of triangles. This may or may not involve inserting or deleting a
  1053. /// vertex.
  1054. ///
  1055. /// The term "casing" is used to describe the set of triangles that are
  1056. /// attached to the triangles being transformed, but are not transformed
  1057. /// themselves. Think of the casing as a fixed hollow structure inside
  1058. /// which all the action happens. A "casing" is only defined relative to
  1059. /// a single transformation; each occurrence of a transformation will
  1060. /// involve a different casing.
  1061. /// </remarks>
  1062. internal void Flip(ref Otri flipedge)
  1063. {
  1064. Otri botleft = default(Otri), botright = default(Otri);
  1065. Otri topleft = default(Otri), topright = default(Otri);
  1066. Otri top = default(Otri);
  1067. Otri botlcasing = default(Otri), botrcasing = default(Otri);
  1068. Otri toplcasing = default(Otri), toprcasing = default(Otri);
  1069. Osub botlsubseg = default(Osub), botrsubseg = default(Osub);
  1070. Osub toplsubseg = default(Osub), toprsubseg = default(Osub);
  1071. Vertex leftvertex, rightvertex, botvertex;
  1072. Vertex farvertex;
  1073. // Identify the vertices of the quadrilateral.
  1074. rightvertex = flipedge.Org();
  1075. leftvertex = flipedge.Dest();
  1076. botvertex = flipedge.Apex();
  1077. flipedge.Sym(ref top);
  1078. // SELF CHECK
  1079. //if (top.triangle.id == DUMMY)
  1080. //{
  1081. // logger.Error("Attempt to flip on boundary.", "Mesh.Flip()");
  1082. // flipedge.LnextSelf();
  1083. // return;
  1084. //}
  1085. //if (checksegments)
  1086. //{
  1087. // flipedge.SegPivot(ref toplsubseg);
  1088. // if (toplsubseg.ss != Segment.Empty)
  1089. // {
  1090. // logger.Error("Attempt to flip a segment.", "Mesh.Flip()");
  1091. // flipedge.LnextSelf();
  1092. // return;
  1093. // }
  1094. //}
  1095. farvertex = top.Apex();
  1096. // Identify the casing of the quadrilateral.
  1097. top.Lprev(ref topleft);
  1098. topleft.Sym(ref toplcasing);
  1099. top.Lnext(ref topright);
  1100. topright.Sym(ref toprcasing);
  1101. flipedge.Lnext(ref botleft);
  1102. botleft.Sym(ref botlcasing);
  1103. flipedge.Lprev(ref botright);
  1104. botright.Sym(ref botrcasing);
  1105. // Rotate the quadrilateral one-quarter turn counterclockwise.
  1106. topleft.Bond(ref botlcasing);
  1107. botleft.Bond(ref botrcasing);
  1108. botright.Bond(ref toprcasing);
  1109. topright.Bond(ref toplcasing);
  1110. if (checksegments)
  1111. {
  1112. // Check for subsegments and rebond them to the quadrilateral.
  1113. topleft.Pivot(ref toplsubseg);
  1114. botleft.Pivot(ref botlsubseg);
  1115. botright.Pivot(ref botrsubseg);
  1116. topright.Pivot(ref toprsubseg);
  1117. if (toplsubseg.seg.hash == DUMMY)
  1118. {
  1119. topright.SegDissolve(dummysub);
  1120. }
  1121. else
  1122. {
  1123. topright.SegBond(ref toplsubseg);
  1124. }
  1125. if (botlsubseg.seg.hash == DUMMY)
  1126. {
  1127. topleft.SegDissolve(dummysub);
  1128. }
  1129. else
  1130. {
  1131. topleft.SegBond(ref botlsubseg);
  1132. }
  1133. if (botrsubseg.seg.hash == DUMMY)
  1134. {
  1135. botleft.SegDissolve(dummysub);
  1136. }
  1137. else
  1138. {
  1139. botleft.SegBond(ref botrsubseg);
  1140. }
  1141. if (toprsubseg.seg.hash == DUMMY)
  1142. {
  1143. botright.SegDissolve(dummysub);
  1144. }
  1145. else
  1146. {
  1147. botright.SegBond(ref toprsubseg);
  1148. }
  1149. }
  1150. // New vertex assignments for the rotated quadrilateral.
  1151. flipedge.SetOrg(farvertex);
  1152. flipedge.SetDest(botvertex);
  1153. flipedge.SetApex(rightvertex);
  1154. top.SetOrg(botvertex);
  1155. top.SetDest(farvertex);
  1156. top.SetApex(leftvertex);
  1157. }
  1158. /// <summary>
  1159. /// Transform two triangles to two different triangles by flipping an edge
  1160. /// clockwise within a quadrilateral. Reverses the flip() operation so that
  1161. /// the data structures representing the triangles are back where they were
  1162. /// before the flip().
  1163. /// </summary>
  1164. /// <param name="flipedge"></param>
  1165. /// <remarks>
  1166. /// See above Flip() remarks for more information.
  1167. ///
  1168. /// Upon completion of this routine, the 'flipedge' handle holds the edge
  1169. /// cd of triangle cdb, and is directed up, from vertex c to vertex d.
  1170. /// (Hence, the two triangles have rotated clockwise.)
  1171. /// </remarks>
  1172. internal void Unflip(ref Otri flipedge)
  1173. {
  1174. Otri botleft = default(Otri), botright = default(Otri);
  1175. Otri topleft = default(Otri), topright = default(Otri);
  1176. Otri top = default(Otri);
  1177. Otri botlcasing = default(Otri), botrcasing = default(Otri);
  1178. Otri toplcasing = default(Otri), toprcasing = default(Otri);
  1179. Osub botlsubseg = default(Osub), botrsubseg = default(Osub);
  1180. Osub toplsubseg = default(Osub), toprsubseg = default(Osub);
  1181. Vertex leftvertex, rightvertex, botvertex;
  1182. Vertex farvertex;
  1183. // Identify the vertices of the quadrilateral.
  1184. rightvertex = flipedge.Org();
  1185. leftvertex = flipedge.Dest();
  1186. botvertex = flipedge.Apex();
  1187. flipedge.Sym(ref top);
  1188. farvertex = top.Apex();
  1189. // Identify the casing of the quadrilateral.
  1190. top.Lprev(ref topleft);
  1191. topleft.Sym(ref toplcasing);
  1192. top.Lnext(ref topright);
  1193. topright.Sym(ref toprcasing);
  1194. flipedge.Lnext(ref botleft);
  1195. botleft.Sym(ref botlcasing);
  1196. flipedge.Lprev(ref botright);
  1197. botright.Sym(ref botrcasing);
  1198. // Rotate the quadrilateral one-quarter turn clockwise.
  1199. topleft.Bond(ref toprcasing);
  1200. botleft.Bond(ref toplcasing);
  1201. botright.Bond(ref botlcasing);
  1202. topright.Bond(ref botrcasing);
  1203. if (checksegments)
  1204. {
  1205. // Check for subsegments and rebond them to the quadrilateral.
  1206. topleft.Pivot(ref toplsubseg);
  1207. botleft.Pivot(ref botlsubseg);
  1208. botright.Pivot(ref botrsubseg);
  1209. topright.Pivot(ref toprsubseg);
  1210. if (toplsubseg.seg.hash == DUMMY)
  1211. {
  1212. botleft.SegDissolve(dummysub);
  1213. }
  1214. else
  1215. {
  1216. botleft.SegBond(ref toplsubseg);
  1217. }
  1218. if (botlsubseg.seg.hash == DUMMY)
  1219. {
  1220. botright.SegDissolve(dummysub);
  1221. }
  1222. else
  1223. {
  1224. botright.SegBond(ref botlsubseg);
  1225. }
  1226. if (botrsubseg.seg.hash == DUMMY)
  1227. {
  1228. topright.SegDissolve(dummysub);
  1229. }
  1230. else
  1231. {
  1232. topright.SegBond(ref botrsubseg);
  1233. }
  1234. if (toprsubseg.seg.hash == DUMMY)
  1235. {
  1236. topleft.SegDissolve(dummysub);
  1237. }
  1238. else
  1239. {
  1240. topleft.SegBond(ref toprsubseg);
  1241. }
  1242. }
  1243. // New vertex assignments for the rotated quadrilateral.
  1244. flipedge.SetOrg(botvertex);
  1245. flipedge.SetDest(farvertex);
  1246. flipedge.SetApex(leftvertex);
  1247. top.SetOrg(farvertex);
  1248. top.SetDest(botvertex);
  1249. top.SetApex(rightvertex);
  1250. }
  1251. /// <summary>
  1252. /// Find the Delaunay triangulation of a polygon that has a certain "nice" shape.
  1253. /// This includes the polygons that result from deletion of a vertex or insertion
  1254. /// of a segment.
  1255. /// </summary>
  1256. /// <param name="firstedge">The primary edge of the first triangle.</param>
  1257. /// <param name="lastedge">The primary edge of the last triangle.</param>
  1258. /// <param name="edgecount">The number of sides of the polygon, including its
  1259. /// base.</param>
  1260. /// <param name="doflip">A flag, wether to perform the last flip.</param>
  1261. /// <param name="triflaws">A flag that determines whether the new triangles should
  1262. /// be tested for quality, and enqueued if they are bad.</param>
  1263. /// <remarks>
  1264. // This is a conceptually difficult routine. The starting assumption is
  1265. // that we have a polygon with n sides. n - 1 of these sides are currently
  1266. // represented as edges in the mesh. One side, called the "base", need not
  1267. // be.
  1268. //
  1269. // Inside the polygon is a structure I call a "fan", consisting of n - 1
  1270. // triangles that share a common origin. For each of these triangles, the
  1271. // edge opposite the origin is one of the sides of the polygon. The
  1272. // primary edge of each triangle is the edge directed from the origin to
  1273. // the destination; note that this is not the same edge that is a side of
  1274. // the polygon. 'firstedge' is the primary edge of the first triangle.
  1275. // From there, the triangles follow in counterclockwise order about the
  1276. // polygon, until 'lastedge', the primary edge of the last triangle.
  1277. // 'firstedge' and 'lastedge' are probably connected to other triangles
  1278. // beyond the extremes of the fan, but their identity is not important, as
  1279. // long as the fan remains connected to them.
  1280. //
  1281. // Imagine the polygon oriented so that its base is at the bottom. This
  1282. // puts 'firstedge' on the far right, and 'lastedge' on the far left.
  1283. // The right vertex of the base is the destination of 'firstedge', and the
  1284. // left vertex of the base is the apex of 'lastedge'.
  1285. //
  1286. // The challenge now is to find the right sequence of edge flips to
  1287. // transform the fan into a Delaunay triangulation of the polygon. Each
  1288. // edge flip effectively removes one triangle from the fan, committing it
  1289. // to the polygon. The resulting polygon has one fewer edge. If 'doflip'
  1290. // is set, the final flip will be performed, resulting in a fan of one
  1291. // (useless?) triangle. If 'doflip' is not set, the final flip is not
  1292. // performed, resulting in a fan of two triangles, and an unfinished
  1293. // triangular polygon that is not yet filled out with a single triangle.
  1294. // On completion of the routine, 'lastedge' is the last remaining triangle,
  1295. // or the leftmost of the last two.
  1296. //
  1297. // Although the flips are performed in the order described above, the
  1298. // decisions about what flips to perform are made in precisely the reverse
  1299. // order. The recursive triangulatepolygon() procedure makes a decision,
  1300. // uses up to two recursive calls to triangulate the "subproblems"
  1301. // (polygons with fewer edges), and then performs an edge flip.
  1302. //
  1303. // The "decision" it makes is which vertex of the polygon should be
  1304. // connected to the base. This decision is made by testing every possible
  1305. // vertex. Once the best vertex is found, the two edges that connect this
  1306. // vertex to the base become the bases for two smaller polygons. These
  1307. // are triangulated recursively. Unfortunately, this approach can take
  1308. // O(n^2) time not only in the worst case, but in many common cases. It's
  1309. // rarely a big deal for vertex deletion, where n is rarely larger than
  1310. // ten, but it could be a big deal for segment insertion, especially if
  1311. // there's a lot of long segments that each cut many triangles. I ought to
  1312. // code a faster algorithm some day.
  1313. /// </remarks>
  1314. private void TriangulatePolygon(Otri firstedge, Otri lastedge,
  1315. int edgecount, bool doflip, bool triflaws)
  1316. {
  1317. Otri testtri = default(Otri);
  1318. Otri besttri = default(Otri);
  1319. Otri tempedge = default(Otri);
  1320. Vertex leftbasevertex, rightbasevertex;
  1321. Vertex testvertex;
  1322. Vertex bestvertex;
  1323. int bestnumber = 1;
  1324. // Identify the base vertices.
  1325. leftbasevertex = lastedge.Apex();
  1326. rightbasevertex = firstedge.Dest();
  1327. // Find the best vertex to connect the base to.
  1328. firstedge.Onext(ref besttri);
  1329. bestvertex = besttri.Dest();
  1330. besttri.Copy(ref testtri);
  1331. for (int i = 2; i <= edgecount - 2; i++)
  1332. {
  1333. testtri.Onext();
  1334. testvertex = testtri.Dest();
  1335. // Is this a better vertex?
  1336. if (predicates.InCircle(leftbasevertex, rightbasevertex, bestvertex, testvertex) > 0.0)
  1337. {
  1338. testtri.Copy(ref besttri);
  1339. bestvertex = testvertex;
  1340. bestnumber = i;
  1341. }
  1342. }
  1343. if (bestnumber > 1)
  1344. {
  1345. // Recursively triangulate the smaller polygon on the right.
  1346. besttri.Oprev(ref tempedge);
  1347. TriangulatePolygon(firstedge, tempedge, bestnumber + 1, true, triflaws);
  1348. }
  1349. if (bestnumber < edgecount - 2)
  1350. {
  1351. // Recursively triangulate the smaller polygon on the left.
  1352. besttri.Sym(ref tempedge);
  1353. TriangulatePolygon(besttri, lastedge, edgecount - bestnumber, true, triflaws);
  1354. // Find 'besttri' again; it may have been lost to edge flips.
  1355. tempedge.Sym(ref besttri);
  1356. }
  1357. if (doflip)
  1358. {
  1359. // Do one final edge flip.
  1360. Flip(ref besttri);
  1361. if (triflaws)
  1362. {
  1363. // Check the quality of the newly committed triangle.
  1364. besttri.Sym(ref testtri);
  1365. qualityMesher.TestTriangle(ref testtri);
  1366. }
  1367. }
  1368. // Return the base triangle.
  1369. besttri.Copy(ref lastedge);
  1370. }
  1371. /// <summary>
  1372. /// Delete a vertex from a Delaunay triangulation, ensuring that the
  1373. /// triangulation remains Delaunay.
  1374. /// </summary>
  1375. /// <param name="deltri"></param>
  1376. /// <remarks>The origin of 'deltri' is deleted. The union of the triangles
  1377. /// adjacent to this vertex is a polygon, for which the Delaunay triangulation
  1378. /// is found. Two triangles are removed from the mesh.
  1379. ///
  1380. /// Only interior vertices that do not lie on segments or boundaries
  1381. /// may be deleted.
  1382. /// </remarks>
  1383. internal void DeleteVertex(ref Otri deltri)
  1384. {
  1385. Otri countingtri = default(Otri);
  1386. Otri firstedge = default(Otri), lastedge = default(Otri);
  1387. Otri deltriright = default(Otri);
  1388. Otri lefttri = default(Otri), righttri = default(Otri);
  1389. Otri leftcasing = default(Otri), rightcasing = default(Otri);
  1390. Osub leftsubseg = default(Osub), rightsubseg = default(Osub);
  1391. Vertex delvertex;
  1392. Vertex neworg;
  1393. int edgecount;
  1394. delvertex = deltri.Org();
  1395. VertexDealloc(delvertex);
  1396. // Count the degree of the vertex being deleted.
  1397. deltri.Onext(ref countingtri);
  1398. edgecount = 1;
  1399. while (!deltri.Equals(countingtri))
  1400. {
  1401. edgecount++;
  1402. countingtri.Onext();
  1403. }
  1404. if (edgecount > 3)
  1405. {
  1406. // Triangulate the polygon defined by the union of all triangles
  1407. // adjacent to the vertex being deleted. Check the quality of
  1408. // the resulting triangles.
  1409. deltri.Onext(ref firstedge);
  1410. deltri.Oprev(ref lastedge);
  1411. TriangulatePolygon(firstedge, lastedge, edgecount, false, behavior.NoBisect == 0);
  1412. }
  1413. // Splice out two triangles.
  1414. deltri.Lprev(ref deltriright);
  1415. deltri.Dnext(ref lefttri);
  1416. lefttri.Sym(ref leftcasing);
  1417. deltriright.Oprev(ref righttri);
  1418. righttri.Sym(ref rightcasing);
  1419. deltri.Bond(ref leftcasing);
  1420. deltriright.Bond(ref rightcasing);
  1421. lefttri.Pivot(ref leftsubseg);
  1422. if (leftsubseg.seg.hash != DUMMY)
  1423. {
  1424. deltri.SegBond(ref leftsubseg);
  1425. }
  1426. righttri.Pivot(ref rightsubseg);
  1427. if (rightsubseg.seg.hash != DUMMY)
  1428. {
  1429. deltriright.SegBond(ref rightsubseg);
  1430. }
  1431. // Set the new origin of 'deltri' and check its quality.
  1432. neworg = lefttri.Org();
  1433. deltri.SetOrg(neworg);
  1434. if (behavior.NoBisect == 0)
  1435. {
  1436. qualityMesher.TestTriangle(ref deltri);
  1437. }
  1438. // Delete the two spliced-out triangles.
  1439. TriangleDealloc(lefttri.tri);
  1440. TriangleDealloc(righttri.tri);
  1441. }
  1442. /// <summary>
  1443. /// Undo the most recent vertex insertion.
  1444. /// </summary>
  1445. /// <remarks>
  1446. /// Walks through the list of transformations (flips and a vertex insertion)
  1447. /// in the reverse of the order in which they were done, and undoes them.
  1448. /// The inserted vertex is removed from the triangulation and deallocated.
  1449. /// Two triangles (possibly just one) are also deallocated.
  1450. /// </remarks>
  1451. internal void UndoVertex()
  1452. {
  1453. Otri fliptri;
  1454. Otri botleft = default(Otri), botright = default(Otri), topright = default(Otri);
  1455. Otri botlcasing = default(Otri), botrcasing = default(Otri), toprcasing = default(Otri);
  1456. Otri gluetri = default(Otri);
  1457. Osub botlsubseg = default(Osub), botrsubseg = default(Osub), toprsubseg = default(Osub);
  1458. Vertex botvertex, rightvertex;
  1459. // Walk through the list of transformations (flips and a vertex insertion)
  1460. // in the reverse of the order in which they were done, and undo them.
  1461. while (flipstack.Count > 0)
  1462. {
  1463. // Find a triangle involved in the last unreversed transformation.
  1464. fliptri = flipstack.Pop();
  1465. // We are reversing one of three transformations: a trisection of one
  1466. // triangle into three (by inserting a vertex in the triangle), a
  1467. // bisection of two triangles into four (by inserting a vertex in an
  1468. // edge), or an edge flip.
  1469. if (flipstack.Count == 0)
  1470. {
  1471. // Restore a triangle that was split into three triangles,
  1472. // so it is again one triangle.
  1473. fliptri.Dprev(ref botleft);
  1474. botleft.Lnext();
  1475. fliptri.Onext(ref botright);
  1476. botright.Lprev();
  1477. botleft.Sym(ref botlcasing);
  1478. botright.Sym(ref botrcasing);
  1479. botvertex = botleft.Dest();
  1480. fliptri.SetApex(botvertex);
  1481. fliptri.Lnext();
  1482. fliptri.Bond(ref botlcasing);
  1483. botleft.Pivot(ref botlsubseg);
  1484. fliptri.SegBond(ref botlsubseg);
  1485. fliptri.Lnext();
  1486. fliptri.Bond(ref botrcasing);
  1487. botright.Pivot(ref botrsubseg);
  1488. fliptri.SegBond(ref botrsubseg);
  1489. // Delete the two spliced-out triangles.
  1490. TriangleDealloc(botleft.tri);
  1491. TriangleDealloc(botright.tri);
  1492. }
  1493. else if (flipstack.Peek().tri == null) // Dummy flip
  1494. {
  1495. // Restore two triangles that were split into four triangles,
  1496. // so they are again two triangles.
  1497. fliptri.Lprev(ref gluetri);
  1498. gluetri.Sym(ref botright);
  1499. botright.Lnext();
  1500. botright.Sym(ref botrcasing);
  1501. rightvertex = botright.Dest();
  1502. fliptri.SetOrg(rightvertex);
  1503. gluetri.Bond(ref botrcasing);
  1504. botright.Pivot(ref botrsubseg);
  1505. gluetri.SegBond(ref botrsubseg);
  1506. // Delete the spliced-out triangle.
  1507. TriangleDealloc(botright.tri);
  1508. fliptri.Sym(ref gluetri);
  1509. if (gluetri.tri.id != DUMMY)
  1510. {
  1511. gluetri.Lnext();
  1512. gluetri.Dnext(ref topright);
  1513. topright.Sym(ref toprcasing);
  1514. gluetri.SetOrg(rightvertex);
  1515. gluetri.Bond(ref toprcasing);
  1516. topright.Pivot(ref toprsubseg);
  1517. gluetri.SegBond(ref toprsubseg);
  1518. // Delete the spliced-out triangle.
  1519. TriangleDealloc(topright.tri);
  1520. }
  1521. flipstack.Clear();
  1522. }
  1523. else
  1524. {
  1525. // Undo an edge flip.
  1526. Unflip(ref fliptri);
  1527. }
  1528. }
  1529. }
  1530. #endregion
  1531. #region Dealloc
  1532. /// <summary>
  1533. /// Deallocate space for a triangle, marking it dead.
  1534. /// </summary>
  1535. /// <param name="dyingtriangle"></param>
  1536. internal void TriangleDealloc(Triangle dyingtriangle)
  1537. {
  1538. // Mark the triangle as dead. This makes it possible to detect dead
  1539. // triangles when traversing the list of all triangles.
  1540. Otri.Kill(dyingtriangle);
  1541. triangles.Release(dyingtriangle);
  1542. }
  1543. /// <summary>
  1544. /// Deallocate space for a vertex, marking it dead.
  1545. /// </summary>
  1546. /// <param name="dyingvertex"></param>
  1547. internal void VertexDealloc(Vertex dyingvertex)
  1548. {
  1549. // Mark the vertex as dead. This makes it possible to detect dead
  1550. // vertices when traversing the list of all vertices.
  1551. dyingvertex.type = VertexType.DeadVertex;
  1552. vertices.Remove(dyingvertex.hash);
  1553. }
  1554. /// <summary>
  1555. /// Deallocate space for a subsegment, marking it dead.
  1556. /// </summary>
  1557. /// <param name="dyingsubseg"></param>
  1558. internal void SubsegDealloc(SubSegment dyingsubseg)
  1559. {
  1560. // Mark the subsegment as dead. This makes it possible to detect dead
  1561. // subsegments when traversing the list of all subsegments.
  1562. Osub.Kill(dyingsubseg);
  1563. subsegs.Remove(dyingsubseg.hash);
  1564. }
  1565. #endregion
  1566. }
  1567. }