TriangleQuadTree.cs 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427
  1. // -----------------------------------------------------------------------
  2. // <copyright file="TriangleQuadTree.cs" company="">
  3. // Original code by Frank Dockhorn, [not available anymore: http://sourceforge.net/projects/quadtreesim/]
  4. // Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
  5. // </copyright>
  6. // -----------------------------------------------------------------------
  7. namespace UnityEngine.U2D.Animation.TriangleNet
  8. .Tools
  9. {
  10. using System.Collections.Generic;
  11. using System.Linq;
  12. using Animation.TriangleNet.Geometry;
  13. /// <summary>
  14. /// A Quadtree implementation optimized for triangles.
  15. /// </summary>
  16. internal class TriangleQuadTree
  17. {
  18. QuadNode root;
  19. internal ITriangle[] triangles;
  20. internal int sizeBound;
  21. internal int maxDepth;
  22. /// <summary>
  23. /// Initializes a new instance of the <see cref="TriangleQuadTree" /> class.
  24. /// </summary>
  25. /// <param name="mesh">Mesh containing triangles.</param>
  26. /// <param name="maxDepth">The maximum depth of the tree.</param>
  27. /// <param name="sizeBound">The maximum number of triangles contained in a leaf.</param>
  28. /// <remarks>
  29. /// The quadtree does not track changes of the mesh. If a mesh is refined or
  30. /// changed in any other way, a new quadtree has to be built to make the point
  31. /// location work.
  32. ///
  33. /// A node of the tree will be split, if its level if less than the max depth parameter
  34. /// AND the number of triangles in the node is greater than the size bound.
  35. /// </remarks>
  36. public TriangleQuadTree(Mesh mesh, int maxDepth = 10, int sizeBound = 10)
  37. {
  38. this.maxDepth = maxDepth;
  39. this.sizeBound = sizeBound;
  40. triangles = mesh.Triangles.ToArray();
  41. int currentDepth = 0;
  42. root = new QuadNode(mesh.Bounds, this, true);
  43. root.CreateSubRegion(++currentDepth);
  44. }
  45. public ITriangle Query(double x, double y)
  46. {
  47. var point = new Point(x, y);
  48. var indices = root.FindTriangles(point);
  49. foreach (var i in indices)
  50. {
  51. var tri = this.triangles[i];
  52. if (IsPointInTriangle(point, tri.GetVertex(0), tri.GetVertex(1), tri.GetVertex(2)))
  53. {
  54. return tri;
  55. }
  56. }
  57. return null;
  58. }
  59. /// <summary>
  60. /// Test, if a given point lies inside a triangle.
  61. /// </summary>
  62. /// <param name="p">Point to locate.</param>
  63. /// <param name="t0">Corner point of triangle.</param>
  64. /// <param name="t1">Corner point of triangle.</param>
  65. /// <param name="t2">Corner point of triangle.</param>
  66. /// <returns>True, if point is inside or on the edge of this triangle.</returns>
  67. internal static bool IsPointInTriangle(Point p, Point t0, Point t1, Point t2)
  68. {
  69. // TODO: no need to create new Point instances here
  70. Point d0 = new Point(t1.x - t0.x, t1.y - t0.y);
  71. Point d1 = new Point(t2.x - t0.x, t2.y - t0.y);
  72. Point d2 = new Point(p.x - t0.x, p.y - t0.y);
  73. // crossproduct of (0, 0, 1) and d0
  74. Point c0 = new Point(-d0.y, d0.x);
  75. // crossproduct of (0, 0, 1) and d1
  76. Point c1 = new Point(-d1.y, d1.x);
  77. // Linear combination d2 = s * d0 + v * d1.
  78. //
  79. // Multiply both sides of the equation with c0 and c1
  80. // and solve for s and v respectively
  81. //
  82. // s = d2 * c1 / d0 * c1
  83. // v = d2 * c0 / d1 * c0
  84. double s = DotProduct(d2, c1) / DotProduct(d0, c1);
  85. double v = DotProduct(d2, c0) / DotProduct(d1, c0);
  86. if (s >= 0 && v >= 0 && ((s + v) <= 1))
  87. {
  88. // Point is inside or on the edge of this triangle.
  89. return true;
  90. }
  91. return false;
  92. }
  93. internal static double DotProduct(Point p, Point q)
  94. {
  95. return p.x * q.x + p.y * q.y;
  96. }
  97. /// <summary>
  98. /// A node of the quadtree.
  99. /// </summary>
  100. class QuadNode
  101. {
  102. const int SW = 0;
  103. const int SE = 1;
  104. const int NW = 2;
  105. const int NE = 3;
  106. const double EPS = 1e-6;
  107. static readonly byte[] BITVECTOR = { 0x1, 0x2, 0x4, 0x8 };
  108. Rectangle bounds;
  109. Point pivot;
  110. TriangleQuadTree tree;
  111. QuadNode[] regions;
  112. List<int> triangles;
  113. byte bitRegions;
  114. public QuadNode(Rectangle box, TriangleQuadTree tree)
  115. : this(box, tree, false)
  116. {
  117. }
  118. public QuadNode(Rectangle box, TriangleQuadTree tree, bool init)
  119. {
  120. this.tree = tree;
  121. this.bounds = new Rectangle(box.Left, box.Bottom, box.Width, box.Height);
  122. this.pivot = new Point((box.Left + box.Right) / 2, (box.Bottom + box.Top) / 2);
  123. this.bitRegions = 0;
  124. this.regions = new QuadNode[4];
  125. this.triangles = new List<int>();
  126. if (init)
  127. {
  128. int count = tree.triangles.Length;
  129. // Allocate memory upfront
  130. triangles.Capacity = count;
  131. for (int i = 0; i < count; i++)
  132. {
  133. triangles.Add(i);
  134. }
  135. }
  136. }
  137. public List<int> FindTriangles(Point searchPoint)
  138. {
  139. int region = FindRegion(searchPoint);
  140. if (regions[region] == null)
  141. {
  142. return triangles;
  143. }
  144. return regions[region].FindTriangles(searchPoint);
  145. }
  146. public void CreateSubRegion(int currentDepth)
  147. {
  148. // The four sub regions of the quad tree
  149. // +--------------+
  150. // | nw 2 | ne 3 |
  151. // |------+pivot--|
  152. // | sw 0 | se 1 |
  153. // +--------------+
  154. Rectangle box;
  155. var width = bounds.Right - pivot.x;
  156. var height = bounds.Top - pivot.y;
  157. // 1. region south west
  158. box = new Rectangle(bounds.Left, bounds.Bottom, width, height);
  159. regions[0] = new QuadNode(box, tree);
  160. // 2. region south east
  161. box = new Rectangle(pivot.x, bounds.Bottom, width, height);
  162. regions[1] = new QuadNode(box, tree);
  163. // 3. region north west
  164. box = new Rectangle(bounds.Left, pivot.y, width, height);
  165. regions[2] = new QuadNode(box, tree);
  166. // 4. region north east
  167. box = new Rectangle(pivot.x, pivot.y, width, height);
  168. regions[3] = new QuadNode(box, tree);
  169. Point[] triangle = new Point[3];
  170. // Find region for every triangle vertex
  171. foreach (var index in triangles)
  172. {
  173. ITriangle tri = tree.triangles[index];
  174. triangle[0] = tri.GetVertex(0);
  175. triangle[1] = tri.GetVertex(1);
  176. triangle[2] = tri.GetVertex(2);
  177. AddTriangleToRegion(triangle, index);
  178. }
  179. for (int i = 0; i < 4; i++)
  180. {
  181. if (regions[i].triangles.Count > tree.sizeBound && currentDepth < tree.maxDepth)
  182. {
  183. regions[i].CreateSubRegion(currentDepth + 1);
  184. }
  185. }
  186. }
  187. void AddTriangleToRegion(Point[] triangle, int index)
  188. {
  189. bitRegions = 0;
  190. if (TriangleQuadTree.IsPointInTriangle(pivot, triangle[0], triangle[1], triangle[2]))
  191. {
  192. AddToRegion(index, SW);
  193. AddToRegion(index, SE);
  194. AddToRegion(index, NW);
  195. AddToRegion(index, NE);
  196. return;
  197. }
  198. FindTriangleIntersections(triangle, index);
  199. if (bitRegions == 0)
  200. {
  201. // we didn't find any intersection so we add this triangle to a point's region
  202. int region = FindRegion(triangle[0]);
  203. regions[region].triangles.Add(index);
  204. }
  205. }
  206. void FindTriangleIntersections(Point[] triangle, int index)
  207. {
  208. // PLEASE NOTE:
  209. // Handling of component comparison is tightly associated with the implementation
  210. // of the findRegion() function. That means when the point to be compared equals
  211. // the pivot point the triangle must be put at least into region 2.
  212. //
  213. // Linear equations are in parametric form.
  214. // pivot.x = triangle[0].x + t * (triangle[1].x - triangle[0].x)
  215. // pivot.y = triangle[0].y + t * (triangle[1].y - triangle[0].y)
  216. int k = 2;
  217. double dx, dy;
  218. // Iterate through all triangle laterals and find bounding box intersections
  219. for (int i = 0; i < 3; k = i++)
  220. {
  221. dx = triangle[i].x - triangle[k].x;
  222. dy = triangle[i].y - triangle[k].y;
  223. if (dx != 0.0)
  224. {
  225. FindIntersectionsWithX(dx, dy, triangle, index, k);
  226. }
  227. if (dy != 0.0)
  228. {
  229. FindIntersectionsWithY(dx, dy, triangle, index, k);
  230. }
  231. }
  232. }
  233. void FindIntersectionsWithX(double dx, double dy, Point[] triangle, int index, int k)
  234. {
  235. double t;
  236. // find intersection with plane x = m_pivot.dX
  237. t = (pivot.x - triangle[k].x) / dx;
  238. if (t < (1 + EPS) && t > -EPS)
  239. {
  240. // we have an intersection
  241. double yComponent = triangle[k].y + t * dy;
  242. if (yComponent < pivot.y && yComponent >= bounds.Bottom)
  243. {
  244. AddToRegion(index, SW);
  245. AddToRegion(index, SE);
  246. }
  247. else if (yComponent <= bounds.Top)
  248. {
  249. AddToRegion(index, NW);
  250. AddToRegion(index, NE);
  251. }
  252. }
  253. // find intersection with plane x = m_boundingBox[0].dX
  254. t = (bounds.Left - triangle[k].x) / dx;
  255. if (t < (1 + EPS) && t > -EPS)
  256. {
  257. // we have an intersection
  258. double yComponent = triangle[k].y + t * dy;
  259. if (yComponent < pivot.y && yComponent >= bounds.Bottom)
  260. {
  261. AddToRegion(index, SW);
  262. }
  263. else if (yComponent <= bounds.Top) // TODO: check && yComponent >= pivot.Y
  264. {
  265. AddToRegion(index, NW);
  266. }
  267. }
  268. // find intersection with plane x = m_boundingBox[1].dX
  269. t = (bounds.Right - triangle[k].x) / dx;
  270. if (t < (1 + EPS) && t > -EPS)
  271. {
  272. // we have an intersection
  273. double yComponent = triangle[k].y + t * dy;
  274. if (yComponent < pivot.y && yComponent >= bounds.Bottom)
  275. {
  276. AddToRegion(index, SE);
  277. }
  278. else if (yComponent <= bounds.Top)
  279. {
  280. AddToRegion(index, NE);
  281. }
  282. }
  283. }
  284. void FindIntersectionsWithY(double dx, double dy, Point[] triangle, int index, int k)
  285. {
  286. double t, xComponent;
  287. // find intersection with plane y = m_pivot.dY
  288. t = (pivot.y - triangle[k].y) / dy;
  289. if (t < (1 + EPS) && t > -EPS)
  290. {
  291. // we have an intersection
  292. xComponent = triangle[k].x + t * dx;
  293. if (xComponent > pivot.x && xComponent <= bounds.Right)
  294. {
  295. AddToRegion(index, SE);
  296. AddToRegion(index, NE);
  297. }
  298. else if (xComponent >= bounds.Left)
  299. {
  300. AddToRegion(index, SW);
  301. AddToRegion(index, NW);
  302. }
  303. }
  304. // find intersection with plane y = m_boundingBox[0].dY
  305. t = (bounds.Bottom - triangle[k].y) / dy;
  306. if (t < (1 + EPS) && t > -EPS)
  307. {
  308. // we have an intersection
  309. xComponent = triangle[k].x + t * dx;
  310. if (xComponent > pivot.x && xComponent <= bounds.Right)
  311. {
  312. AddToRegion(index, SE);
  313. }
  314. else if (xComponent >= bounds.Left)
  315. {
  316. AddToRegion(index, SW);
  317. }
  318. }
  319. // find intersection with plane y = m_boundingBox[1].dY
  320. t = (bounds.Top - triangle[k].y) / dy;
  321. if (t < (1 + EPS) && t > -EPS)
  322. {
  323. // we have an intersection
  324. xComponent = triangle[k].x + t * dx;
  325. if (xComponent > pivot.x && xComponent <= bounds.Right)
  326. {
  327. AddToRegion(index, NE);
  328. }
  329. else if (xComponent >= bounds.Left)
  330. {
  331. AddToRegion(index, NW);
  332. }
  333. }
  334. }
  335. int FindRegion(Point point)
  336. {
  337. int b = 2;
  338. if (point.y < pivot.y)
  339. {
  340. b = 0;
  341. }
  342. if (point.x > pivot.x)
  343. {
  344. b++;
  345. }
  346. return b;
  347. }
  348. void AddToRegion(int index, int region)
  349. {
  350. //if (!(m_bitRegions & BITVECTOR[region]))
  351. if ((bitRegions & BITVECTOR[region]) == 0)
  352. {
  353. regions[region].triangles.Add(index);
  354. bitRegions |= BITVECTOR[region];
  355. }
  356. }
  357. }
  358. }
  359. }