CurveExtractorSquares.h 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470
  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 3 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 CurveExtractorSquares : 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. CurveExtractorSquares(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. for (int x = 0, xp = 1; xp < this->mXBound; ++x, ++xp)
  53. {
  54. // Get the image values at the corners of the square.
  55. int i00 = x + this->mXBound * y;
  56. int i10 = i00 + 1;
  57. int i01 = i00 + this->mXBound;
  58. int i11 = i10 + this->mXBound;
  59. int64_t f00 = this->mPixels[i00];
  60. int64_t f10 = this->mPixels[i10];
  61. int64_t f01 = this->mPixels[i01];
  62. int64_t f11 = this->mPixels[i11];
  63. // Construct the vertices and edges of the level curve in
  64. // the square. The x, xp, y and yp values are implicitly
  65. // converted from int to int64_t (which is guaranteed to
  66. // be correct).
  67. ProcessSquare(vertices, edges, x, xp, y, yp, f00, f10, f11, f01);
  68. }
  69. }
  70. }
  71. protected:
  72. void ProcessSquare(std::vector<Vertex>& vertices, std::vector<Edge>& edges,
  73. int64_t x, int64_t xp, int64_t y, int64_t yp,
  74. int64_t f00, int64_t f10, int64_t f11, int64_t f01)
  75. {
  76. int64_t xn0, yn0, xn1, yn1, d0, d1, d2, d3, det;
  77. if (f00 != 0)
  78. {
  79. // convert to case "+***"
  80. if (f00 < 0)
  81. {
  82. f00 = -f00;
  83. f10 = -f10;
  84. f11 = -f11;
  85. f01 = -f01;
  86. }
  87. if (f10 > 0)
  88. {
  89. if (f11 > 0)
  90. {
  91. if (f01 > 0)
  92. {
  93. // ++++
  94. return;
  95. }
  96. else if (f01 < 0)
  97. {
  98. // +++-
  99. d0 = f11 - f01;
  100. xn0 = f11 * x - f01 * xp;
  101. d1 = f00 - f01;
  102. yn1 = f00 * yp - f01 * y;
  103. this->AddEdge(vertices, edges, xn0, d0, yp, 1, x, 1, yn1, d1);
  104. }
  105. else
  106. {
  107. // +++0
  108. this->AddVertex(vertices, x, 1, yp, 1);
  109. }
  110. }
  111. else if (f11 < 0)
  112. {
  113. d0 = f10 - f11;
  114. yn0 = f10 * yp - f11 * y;
  115. if (f01 > 0)
  116. {
  117. // ++-+
  118. d1 = f01 - f11;
  119. xn1 = f01 * xp - f11 * x;
  120. this->AddEdge(vertices, edges, xp, 1, yn0, d0, xn1, d1, yp, 1);
  121. }
  122. else if (f01 < 0)
  123. {
  124. // ++--
  125. d1 = f01 - f00;
  126. yn1 = f01 * y - f00 * yp;
  127. this->AddEdge(vertices, edges, x, 1, yn1, d1, xp, 1, yn0, d0);
  128. }
  129. else
  130. {
  131. // ++-0
  132. this->AddEdge(vertices, edges, x, 1, yp, 1, xp, 1, yn0, d0);
  133. }
  134. }
  135. else
  136. {
  137. if (f01 > 0)
  138. {
  139. // ++0+
  140. this->AddVertex(vertices, xp, 1, yp, 1);
  141. }
  142. else if (f01 < 0)
  143. {
  144. // ++0-
  145. d0 = f01 - f00;
  146. yn0 = f01 * y - f00 * yp;
  147. this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yn0, d0);
  148. }
  149. else
  150. {
  151. // ++00
  152. this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yp, 1);
  153. }
  154. }
  155. }
  156. else if (f10 < 0)
  157. {
  158. d0 = f00 - f10;
  159. xn0 = f00 * xp - f10 * x;
  160. if (f11 > 0)
  161. {
  162. d1 = f11 - f10;
  163. yn1 = f11 * y - f10 * yp;
  164. if (f01 > 0)
  165. {
  166. // +-++
  167. this->AddEdge(vertices, edges, xn0, d0, y, 1, xp, 1, yn1, d1);
  168. }
  169. else if (f01 < 0)
  170. {
  171. // +-+-
  172. d3 = f11 - f01;
  173. xn1 = f11 * x - f01 * xp;
  174. d2 = f01 - f00;
  175. yn0 = f01 * y - f00 * yp;
  176. if (d0*d3 > 0)
  177. {
  178. det = xn1 * d0 - xn0 * d3;
  179. }
  180. else
  181. {
  182. det = xn0 * d3 - xn1 * d0;
  183. }
  184. if (det > 0)
  185. {
  186. this->AddEdge(vertices, edges, xn1, d3, yp, 1, xp, 1, yn1, d1);
  187. this->AddEdge(vertices, edges, xn0, d0, y, 1, x, 1, yn0, d2);
  188. }
  189. else if (det < 0)
  190. {
  191. this->AddEdge(vertices, edges, xn1, d3, yp, 1, x, 1, yn0, d2);
  192. this->AddEdge(vertices, edges, xn0, d0, y, 1, xp, 1, yn1, d1);
  193. }
  194. else
  195. {
  196. this->AddEdge(vertices, edges, xn0, d0, yn0, d2, xn0, d0, y, 1);
  197. this->AddEdge(vertices, edges, xn0, d0, yn0, d2, xn0, d0, yp, 1);
  198. this->AddEdge(vertices, edges, xn0, d0, yn0, d2, x, 1, yn0, d2);
  199. this->AddEdge(vertices, edges, xn0, d0, yn0, d2, xp, 1, yn0, d2);
  200. }
  201. }
  202. else
  203. {
  204. // +-+0
  205. this->AddEdge(vertices, edges, xn0, d0, y, 1, xp, 1, yn1, d1);
  206. this->AddVertex(vertices, x, 1, yp, 1);
  207. }
  208. }
  209. else if (f11 < 0)
  210. {
  211. if (f01 > 0)
  212. {
  213. // +--+
  214. d1 = f11 - f01;
  215. xn1 = f11 * x - f01 * xp;
  216. this->AddEdge(vertices, edges, xn0, d0, y, 1, xn1, d1, yp, 1);
  217. }
  218. else if (f01 < 0)
  219. {
  220. // +---
  221. d1 = f01 - f00;
  222. yn1 = f01 * y - f00 * yp;
  223. this->AddEdge(vertices, edges, x, 1, yn1, d1, xn0, d0, y, 1);
  224. }
  225. else
  226. {
  227. // +--0
  228. this->AddEdge(vertices, edges, x, 1, yp, 1, xn0, d0, y, 1);
  229. }
  230. }
  231. else
  232. {
  233. if (f01 > 0)
  234. {
  235. // +-0+
  236. this->AddEdge(vertices, edges, xp, 1, yp, 1, xn0, d0, y, 1);
  237. }
  238. else if (f01 < 0)
  239. {
  240. // +-0-
  241. d1 = f01 - f00;
  242. yn1 = f01 * y - f00 * yp;
  243. this->AddEdge(vertices, edges, x, 1, yn1, d1, xn0, d0, y, 1);
  244. this->AddVertex(vertices, xp, 1, yp, 1);
  245. }
  246. else
  247. {
  248. // +-00
  249. this->AddEdge(vertices, edges, xp, 1, yp, 1, xn0, d0, yp, 1);
  250. this->AddEdge(vertices, edges, xn0, d0, yp, 1, x, 1, yp, 1);
  251. this->AddEdge(vertices, edges, xn0, d0, yp, 1, xn0, d0, y, 1);
  252. }
  253. }
  254. }
  255. else
  256. {
  257. if (f11 > 0)
  258. {
  259. if (f01 > 0)
  260. {
  261. // +0++
  262. this->AddVertex(vertices, xp, 1, y, 1);
  263. }
  264. else if (f01 < 0)
  265. {
  266. // +0+-
  267. d0 = f11 - f01;
  268. xn0 = f11 * x - f01 * xp;
  269. d1 = f00 - f01;
  270. yn1 = f00 * yp - f01 * y;
  271. this->AddEdge(vertices, edges, xn0, d0, yp, 1, x, 1, yn1, d1);
  272. this->AddVertex(vertices, xp, 1, y, 1);
  273. }
  274. else
  275. {
  276. // +0+0
  277. this->AddVertex(vertices, xp, 1, y, 1);
  278. this->AddVertex(vertices, x, 1, yp, 1);
  279. }
  280. }
  281. else if (f11 < 0)
  282. {
  283. if (f01 > 0)
  284. {
  285. // +0-+
  286. d0 = f11 - f01;
  287. xn0 = f11 * x - f01 * xp;
  288. this->AddEdge(vertices, edges, xp, 1, y, 1, xn0, d0, yp, 1);
  289. }
  290. else if (f01 < 0)
  291. {
  292. // +0--
  293. d0 = f01 - f00;
  294. yn0 = f01 * y - f00 * yp;
  295. this->AddEdge(vertices, edges, xp, 1, y, 1, x, 1, yn0, d0);
  296. }
  297. else
  298. {
  299. // +0-0
  300. this->AddEdge(vertices, edges, xp, 1, y, 1, x, 1, yp, 1);
  301. }
  302. }
  303. else
  304. {
  305. if (f01 > 0)
  306. {
  307. // +00+
  308. this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yp, 1);
  309. }
  310. else if (f01 < 0)
  311. {
  312. // +00-
  313. d0 = f00 - f01;
  314. yn0 = f00 * yp - f01 * y;
  315. this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yn0, d0);
  316. this->AddEdge(vertices, edges, xp, 1, yn0, d0, xp, 1, yp, 1);
  317. this->AddEdge(vertices, edges, xp, 1, yn0, d0, x, 1, yn0, d0);
  318. }
  319. else
  320. {
  321. // +000
  322. this->AddEdge(vertices, edges, x, 1, yp, 1, x, 1, y, 1);
  323. this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, y, 1);
  324. }
  325. }
  326. }
  327. }
  328. else if (f10 != 0)
  329. {
  330. // convert to case 0+**
  331. if (f10 < 0)
  332. {
  333. f10 = -f10;
  334. f11 = -f11;
  335. f01 = -f01;
  336. }
  337. if (f11 > 0)
  338. {
  339. if (f01 > 0)
  340. {
  341. // 0+++
  342. this->AddVertex(vertices, x, 1, y, 1);
  343. }
  344. else if (f01 < 0)
  345. {
  346. // 0++-
  347. d0 = f11 - f01;
  348. xn0 = f11 * x - f01 * xp;
  349. this->AddEdge(vertices, edges, x, 1, y, 1, xn0, d0, yp, 1);
  350. }
  351. else
  352. {
  353. // 0++0
  354. this->AddEdge(vertices, edges, x, 1, yp, 1, x, 1, y, 1);
  355. }
  356. }
  357. else if (f11 < 0)
  358. {
  359. if (f01 > 0)
  360. {
  361. // 0+-+
  362. d0 = f10 - f11;
  363. yn0 = f10 * yp - f11 * y;
  364. d1 = f01 - f11;
  365. xn1 = f01 * xp - f11 * x;
  366. this->AddEdge(vertices, edges, xp, 1, yn0, d0, xn1, d1, yp, 1);
  367. this->AddVertex(vertices, x, 1, y, 1);
  368. }
  369. else if (f01 < 0)
  370. {
  371. // 0+--
  372. d0 = f10 - f11;
  373. yn0 = f10 * yp - f11 * y;
  374. this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, yn0, d0);
  375. }
  376. else
  377. {
  378. // 0+-0
  379. d0 = f10 - f11;
  380. yn0 = f10 * yp - f11 * y;
  381. this->AddEdge(vertices, edges, x, 1, y, 1, x, 1, yn0, d0);
  382. this->AddEdge(vertices, edges, x, 1, yn0, d0, x, 1, yp, 1);
  383. this->AddEdge(vertices, edges, x, 1, yn0, d0, xp, 1, yn0, d0);
  384. }
  385. }
  386. else
  387. {
  388. if (f01 > 0)
  389. {
  390. // 0+0+
  391. this->AddVertex(vertices, x, 1, y, 1);
  392. this->AddVertex(vertices, xp, 1, yp, 1);
  393. }
  394. else if (f01 < 0)
  395. {
  396. // 0+0-
  397. this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, yp, 1);
  398. }
  399. else
  400. {
  401. // 0+00
  402. this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yp, 1);
  403. this->AddEdge(vertices, edges, x, 1, yp, 1, x, 1, y, 1);
  404. }
  405. }
  406. }
  407. else if (f11 != 0)
  408. {
  409. // convert to case 00+*
  410. if (f11 < 0)
  411. {
  412. f11 = -f11;
  413. f01 = -f01;
  414. }
  415. if (f01 > 0)
  416. {
  417. // 00++
  418. this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, y, 1);
  419. }
  420. else if (f01 < 0)
  421. {
  422. // 00+-
  423. d0 = f01 - f11;
  424. xn0 = f01 * xp - f11 * x;
  425. this->AddEdge(vertices, edges, x, 1, y, 1, xn0, d0, y, 1);
  426. this->AddEdge(vertices, edges, xn0, d0, y, 1, xp, 1, y, 1);
  427. this->AddEdge(vertices, edges, xn0, d0, y, 1, xn0, d0, yp, 1);
  428. }
  429. else
  430. {
  431. // 00+0
  432. this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yp, 1);
  433. this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yp, 1);
  434. }
  435. }
  436. else if (f01 != 0)
  437. {
  438. // cases 000+ or 000-
  439. this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, y, 1);
  440. this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yp, 1);
  441. }
  442. else
  443. {
  444. // case 0000
  445. this->AddEdge(vertices, edges, x, 1, y, 1, xp, 1, y, 1);
  446. this->AddEdge(vertices, edges, xp, 1, y, 1, xp, 1, yp, 1);
  447. this->AddEdge(vertices, edges, xp, 1, yp, 1, x, 1, yp, 1);
  448. this->AddEdge(vertices, edges, x, 1, yp, 1, x, 1, y, 1);
  449. }
  450. }
  451. };
  452. }