CurveExtractorTriangles.h 8.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227
  1. // David Eberly, Geometric Tools, Redmond WA 98052
  2. // Copyright (c) 1998-2020
  3. // Distributed under the Boost Software License, Version 1.0.
  4. // https://www.boost.org/LICENSE_1_0.txt
  5. // https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
  6. // Version: 4.0.2019.08.13
  7. #pragma once
  8. #include <Mathematics/CurveExtractor.h>
  9. // The level set extraction algorithm implemented here is described
  10. // in Section 2 of the document
  11. // https://www.geometrictools.com/Documentation/ExtractLevelCurves.pdf
  12. namespace WwiseGTE
  13. {
  14. // The image type T must be one of the integer types: int8_t, int16_t,
  15. // int32_t, uint8_t, uint16_t or uint32_t. Internal integer computations
  16. // are performed using int64_t. The type Real is for extraction to
  17. // floating-point vertices.
  18. template <typename T, typename Real>
  19. class CurveExtractorTriangles : public CurveExtractor<T, Real>
  20. {
  21. public:
  22. // Convenience type definitions.
  23. typedef typename CurveExtractor<T, Real>::Vertex Vertex;
  24. typedef typename CurveExtractor<T, Real>::Edge Edge;
  25. // The input is a 2D image with lexicographically ordered pixels (x,y)
  26. // stored in a linear array. Pixel (x,y) is stored in the array at
  27. // location index = x + xBound * y. The inputs xBound and yBound must
  28. // each be 2 or larger so that there is at least one image square to
  29. // process. The inputPixels must be nonnull and point to contiguous
  30. // storage that contains at least xBound * yBound elements.
  31. CurveExtractorTriangles(int xBound, int yBound, T const* inputPixels)
  32. :
  33. CurveExtractor<T, Real>(xBound, yBound, inputPixels)
  34. {
  35. }
  36. // Extract level curves and return rational vertices. Use the
  37. // base-class Extract if you want real-valued vertices.
  38. virtual void Extract(T level, std::vector<Vertex>& vertices,
  39. std::vector<Edge>& edges) override
  40. {
  41. // Adjust the image so that the level set is F(x,y) = 0.
  42. int64_t levelI64 = static_cast<int64_t>(level);
  43. for (size_t i = 0; i < this->mPixels.size(); ++i)
  44. {
  45. int64_t inputI64 = static_cast<int64_t>(this->mInputPixels[i]);
  46. this->mPixels[i] = inputI64 - levelI64;
  47. }
  48. vertices.clear();
  49. edges.clear();
  50. for (int y = 0, yp = 1; yp < this->mYBound; ++y, ++yp)
  51. {
  52. int yParity = (y & 1);
  53. for (int x = 0, xp = 1; xp < this->mXBound; ++x, ++xp)
  54. {
  55. int xParity = (x & 1);
  56. // Get the image values at the corners of the square.
  57. int i00 = x + this->mXBound * y;
  58. int i10 = i00 + 1;
  59. int i01 = i00 + this->mXBound;
  60. int i11 = i10 + this->mXBound;
  61. int64_t f00 = this->mPixels[i00];
  62. int64_t f10 = this->mPixels[i10];
  63. int64_t f01 = this->mPixels[i01];
  64. int64_t f11 = this->mPixels[i11];
  65. // Construct the vertices and edges of the level curve in
  66. // the square. The x, xp, y and yp values are implicitly
  67. // converted from int to int64_t (which is guaranteed to
  68. // be correct).
  69. if (xParity == yParity)
  70. {
  71. ProcessTriangle(vertices, edges, x, y, f00, x, yp, f01, xp, y, f10);
  72. ProcessTriangle(vertices, edges, xp, yp, f11, xp, y, f10, x, yp, f01);
  73. }
  74. else
  75. {
  76. ProcessTriangle(vertices, edges, x, yp, f01, xp, yp, f11, x, y, f00);
  77. ProcessTriangle(vertices, edges, xp, y, f10, x, y, f00, xp, yp, f11);
  78. }
  79. }
  80. }
  81. }
  82. protected:
  83. void ProcessTriangle(std::vector<Vertex>& vertices, std::vector<Edge>& edges,
  84. int64_t x0, int64_t y0, int64_t f0,
  85. int64_t x1, int64_t y1, int64_t f1,
  86. int64_t x2, int64_t y2, int64_t f2)
  87. {
  88. int64_t xn0, yn0, xn1, yn1, d0, d1;
  89. if (f0 != 0)
  90. {
  91. // convert to case "+**"
  92. if (f0 < 0)
  93. {
  94. f0 = -f0;
  95. f1 = -f1;
  96. f2 = -f2;
  97. }
  98. if (f1 > 0)
  99. {
  100. if (f2 > 0)
  101. {
  102. // +++
  103. return;
  104. }
  105. else if (f2 < 0)
  106. {
  107. // ++-
  108. d0 = f0 - f2;
  109. xn0 = f0 * x2 - f2 * x0;
  110. yn0 = f0 * y2 - f2 * y0;
  111. d1 = f1 - f2;
  112. xn1 = f1 * x2 - f2 * x1;
  113. yn1 = f1 * y2 - f2 * y1;
  114. this->AddEdge(vertices, edges, xn0, d0, yn0, d0, xn1, d1, yn1, d1);
  115. }
  116. else
  117. {
  118. // ++0
  119. this->AddVertex(vertices, x2, 1, y2, 1);
  120. }
  121. }
  122. else if (f1 < 0)
  123. {
  124. d0 = f0 - f1;
  125. xn0 = f0 * x1 - f1 * x0;
  126. yn0 = f0 * y1 - f1 * y0;
  127. if (f2 > 0)
  128. {
  129. // +-+
  130. d1 = f2 - f1;
  131. xn1 = f2 * x1 - f1 * x2;
  132. yn1 = f2 * y1 - f1 * y2;
  133. this->AddEdge(vertices, edges, xn0, d0, yn0, d0, xn1, d1, yn1, d1);
  134. }
  135. else if (f2 < 0)
  136. {
  137. // +--
  138. d1 = f2 - f0;
  139. xn1 = f2 * x0 - f0 * x2;
  140. yn1 = f2 * y0 - f0 * y2;
  141. this->AddEdge(vertices, edges, xn0, d0, yn0, d0, xn1, d1, yn1, d1);
  142. }
  143. else
  144. {
  145. // +-0
  146. this->AddEdge(vertices, edges, x2, 1, y2, 1, xn0, d0, yn0, d0);
  147. }
  148. }
  149. else
  150. {
  151. if (f2 > 0)
  152. {
  153. // +0+
  154. this->AddVertex(vertices, x1, 1, y1, 1);
  155. }
  156. else if (f2 < 0)
  157. {
  158. // +0-
  159. d0 = f2 - f0;
  160. xn0 = f2 * x0 - f0 * x2;
  161. yn0 = f2 * y0 - f0 * y2;
  162. this->AddEdge(vertices, edges, x1, 1, y1, 1, xn0, d0, yn0, d0);
  163. }
  164. else
  165. {
  166. // +00
  167. this->AddEdge(vertices, edges, x1, 1, y1, 1, x2, 1, y2, 1);
  168. }
  169. }
  170. }
  171. else if (f1 != 0)
  172. {
  173. // convert to case 0+*
  174. if (f1 < 0)
  175. {
  176. f1 = -f1;
  177. f2 = -f2;
  178. }
  179. if (f2 > 0)
  180. {
  181. // 0++
  182. this->AddVertex(vertices, x0, 1, y0, 1);
  183. }
  184. else if (f2 < 0)
  185. {
  186. // 0+-
  187. d0 = f1 - f2;
  188. xn0 = f1 * x2 - f2 * x1;
  189. yn0 = f1 * y2 - f2 * y1;
  190. this->AddEdge(vertices, edges, x0, 1, y0, 1, xn0, d0, yn0, d0);
  191. }
  192. else
  193. {
  194. // 0+0
  195. this->AddEdge(vertices, edges, x0, 1, y0, 1, x2, 1, y2, 1);
  196. }
  197. }
  198. else if (f2 != 0)
  199. {
  200. // cases 00+ or 00-
  201. this->AddEdge(vertices, edges, x0, 1, y0, 1, x1, 1, y1, 1);
  202. }
  203. else
  204. {
  205. // case 000
  206. this->AddEdge(vertices, edges, x0, 1, y0, 1, x1, 1, y1, 1);
  207. this->AddEdge(vertices, edges, x1, 1, y1, 1, x2, 1, y2, 1);
  208. this->AddEdge(vertices, edges, x2, 1, y2, 1, x0, 1, y0, 1);
  209. }
  210. }
  211. };
  212. }