GVector.h 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510
  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/Logger.h>
  9. #include <Mathematics/Math.h>
  10. #include <vector>
  11. namespace WwiseGTE
  12. {
  13. template <typename Real>
  14. class GVector
  15. {
  16. public:
  17. // The tuple is length zero (uninitialized).
  18. GVector() = default;
  19. // The tuple is length 'size' and the elements are uninitialized.
  20. GVector(int size)
  21. {
  22. SetSize(size);
  23. }
  24. // For 0 <= d <= size, element d is 1 and all others are zero. If d
  25. // is invalid, the zero vector is created. This is a convenience for
  26. // creating the standard Euclidean basis vectors; see also
  27. // MakeUnit(int,int) and Unit(int,int).
  28. GVector(int size, int d)
  29. {
  30. SetSize(size);
  31. MakeUnit(d);
  32. }
  33. // The copy constructor, destructor, and assignment operator are
  34. // generated by the compiler.
  35. // Member access. SetSize(int) does not initialize the tuple. The
  36. // first operator[] returns a const reference rather than a Real
  37. // value. This supports writing via standard file operations that
  38. // require a const pointer to data.
  39. void SetSize(int size)
  40. {
  41. LogAssert(size >= 0, "Invalid size.");
  42. mTuple.resize(size);
  43. }
  44. inline int GetSize() const
  45. {
  46. return static_cast<int>(mTuple.size());
  47. }
  48. inline Real const& operator[](int i) const
  49. {
  50. return mTuple[i];
  51. }
  52. inline Real& operator[](int i)
  53. {
  54. return mTuple[i];
  55. }
  56. // Comparison (for use by STL containers).
  57. inline bool operator==(GVector const& vec) const
  58. {
  59. return mTuple == vec.mTuple;
  60. }
  61. inline bool operator!=(GVector const& vec) const
  62. {
  63. return mTuple != vec.mTuple;
  64. }
  65. inline bool operator< (GVector const& vec) const
  66. {
  67. return mTuple < vec.mTuple;
  68. }
  69. inline bool operator<=(GVector const& vec) const
  70. {
  71. return mTuple <= vec.mTuple;
  72. }
  73. inline bool operator> (GVector const& vec) const
  74. {
  75. return mTuple > vec.mTuple;
  76. }
  77. inline bool operator>=(GVector const& vec) const
  78. {
  79. return mTuple >= vec.mTuple;
  80. }
  81. // Special vectors.
  82. // All components are 0.
  83. void MakeZero()
  84. {
  85. std::fill(mTuple.begin(), mTuple.end(), (Real)0);
  86. }
  87. // Component d is 1, all others are zero.
  88. void MakeUnit(int d)
  89. {
  90. std::fill(mTuple.begin(), mTuple.end(), (Real)0);
  91. if (0 <= d && d < (int)mTuple.size())
  92. {
  93. mTuple[d] = (Real)1;
  94. }
  95. }
  96. static GVector Zero(int size)
  97. {
  98. GVector<Real> v(size);
  99. v.MakeZero();
  100. return v;
  101. }
  102. static GVector Unit(int size, int d)
  103. {
  104. GVector<Real> v(size);
  105. v.MakeUnit(d);
  106. return v;
  107. }
  108. protected:
  109. // This data structure takes advantage of the built-in operator[],
  110. // range checking and visualizers in MSVS.
  111. std::vector<Real> mTuple;
  112. };
  113. // Unary operations.
  114. template <typename Real>
  115. GVector<Real> operator+(GVector<Real> const& v)
  116. {
  117. return v;
  118. }
  119. template <typename Real>
  120. GVector<Real> operator-(GVector<Real> const& v)
  121. {
  122. GVector<Real> result(v.GetSize());
  123. for (int i = 0; i < v.GetSize(); ++i)
  124. {
  125. result[i] = -v[i];
  126. }
  127. return result;
  128. }
  129. // Linear-algebraic operations.
  130. template <typename Real>
  131. GVector<Real> operator+(GVector<Real> const& v0, GVector<Real> const& v1)
  132. {
  133. GVector<Real> result = v0;
  134. return result += v1;
  135. }
  136. template <typename Real>
  137. GVector<Real> operator-(GVector<Real> const& v0, GVector<Real> const& v1)
  138. {
  139. GVector<Real> result = v0;
  140. return result -= v1;
  141. }
  142. template <typename Real>
  143. GVector<Real> operator*(GVector<Real> const& v, Real scalar)
  144. {
  145. GVector<Real> result = v;
  146. return result *= scalar;
  147. }
  148. template <typename Real>
  149. GVector<Real> operator*(Real scalar, GVector<Real> const& v)
  150. {
  151. GVector<Real> result = v;
  152. return result *= scalar;
  153. }
  154. template <typename Real>
  155. GVector<Real> operator/(GVector<Real> const& v, Real scalar)
  156. {
  157. GVector<Real> result = v;
  158. return result /= scalar;
  159. }
  160. template <typename Real>
  161. GVector<Real>& operator+=(GVector<Real>& v0, GVector<Real> const& v1)
  162. {
  163. if (v0.GetSize() == v1.GetSize())
  164. {
  165. for (int i = 0; i < v0.GetSize(); ++i)
  166. {
  167. v0[i] += v1[i];
  168. }
  169. return v0;
  170. }
  171. LogError("Mismatched sizes.");
  172. }
  173. template <typename Real>
  174. GVector<Real>& operator-=(GVector<Real>& v0, GVector<Real> const& v1)
  175. {
  176. if (v0.GetSize() == v1.GetSize())
  177. {
  178. for (int i = 0; i < v0.GetSize(); ++i)
  179. {
  180. v0[i] -= v1[i];
  181. }
  182. return v0;
  183. }
  184. LogError("Mismatched sizes.");
  185. }
  186. template <typename Real>
  187. GVector<Real>& operator*=(GVector<Real>& v, Real scalar)
  188. {
  189. for (int i = 0; i < v.GetSize(); ++i)
  190. {
  191. v[i] *= scalar;
  192. }
  193. return v;
  194. }
  195. template <typename Real>
  196. GVector<Real>& operator/=(GVector<Real>& v, Real scalar)
  197. {
  198. if (scalar != (Real)0)
  199. {
  200. Real invScalar = (Real)1 / scalar;
  201. for (int i = 0; i < v.GetSize(); ++i)
  202. {
  203. v[i] *= invScalar;
  204. }
  205. return v;
  206. }
  207. LogError("Division by zero.");
  208. }
  209. // Geometric operations. The functions with 'robust' set to 'false' use
  210. // the standard algorithm for normalizing a vector by computing the length
  211. // as a square root of the squared length and dividing by it. The results
  212. // can be infinite (or NaN) if the length is zero. When 'robust' is set
  213. // to 'true', the algorithm is designed to avoid floating-point overflow
  214. // and sets the normalized vector to zero when the length is zero.
  215. template <typename Real>
  216. Real Dot(GVector<Real> const& v0, GVector<Real> const& v1)
  217. {
  218. if (v0.GetSize() == v1.GetSize())
  219. {
  220. Real dot(0);
  221. for (int i = 0; i < v0.GetSize(); ++i)
  222. {
  223. dot += v0[i] * v1[i];
  224. }
  225. return dot;
  226. }
  227. LogError("Mismatched sizes.");
  228. }
  229. template <typename Real>
  230. Real Length(GVector<Real> const& v, bool robust = false)
  231. {
  232. if (robust)
  233. {
  234. Real maxAbsComp = std::fabs(v[0]);
  235. for (int i = 1; i < v.GetSize(); ++i)
  236. {
  237. Real absComp = std::fabs(v[i]);
  238. if (absComp > maxAbsComp)
  239. {
  240. maxAbsComp = absComp;
  241. }
  242. }
  243. Real length;
  244. if (maxAbsComp > (Real)0)
  245. {
  246. GVector<Real> scaled = v / maxAbsComp;
  247. length = maxAbsComp * std::sqrt(Dot(scaled, scaled));
  248. }
  249. else
  250. {
  251. length = (Real)0;
  252. }
  253. return length;
  254. }
  255. else
  256. {
  257. return std::sqrt(Dot(v, v));
  258. }
  259. }
  260. template <typename Real>
  261. Real Normalize(GVector<Real>& v, bool robust = false)
  262. {
  263. if (robust)
  264. {
  265. Real maxAbsComp = std::fabs(v[0]);
  266. for (int i = 1; i < v.GetSize(); ++i)
  267. {
  268. Real absComp = std::fabs(v[i]);
  269. if (absComp > maxAbsComp)
  270. {
  271. maxAbsComp = absComp;
  272. }
  273. }
  274. Real length;
  275. if (maxAbsComp > (Real)0)
  276. {
  277. v /= maxAbsComp;
  278. length = std::sqrt(Dot(v, v));
  279. v /= length;
  280. length *= maxAbsComp;
  281. }
  282. else
  283. {
  284. length = (Real)0;
  285. for (int i = 0; i < v.GetSize(); ++i)
  286. {
  287. v[i] = (Real)0;
  288. }
  289. }
  290. return length;
  291. }
  292. else
  293. {
  294. Real length = std::sqrt(Dot(v, v));
  295. if (length > (Real)0)
  296. {
  297. v /= length;
  298. }
  299. else
  300. {
  301. for (int i = 0; i < v.GetSize(); ++i)
  302. {
  303. v[i] = (Real)0;
  304. }
  305. }
  306. return length;
  307. }
  308. }
  309. // Gram-Schmidt orthonormalization to generate orthonormal vectors from
  310. // the linearly independent inputs. The function returns the smallest
  311. // length of the unnormalized vectors computed during the process. If
  312. // this value is nearly zero, it is possible that the inputs are linearly
  313. // dependent (within numerical round-off errors). On input,
  314. // 1 <= numElements <= N and v[0] through v[numElements-1] must be
  315. // initialized. On output, the vectors v[0] through v[numElements-1]
  316. // form an orthonormal set.
  317. template <typename Real>
  318. Real Orthonormalize(int numInputs, GVector<Real>* v, bool robust = false)
  319. {
  320. if (v && 1 <= numInputs && numInputs <= v[0].GetSize())
  321. {
  322. for (int i = 1; i < numInputs; ++i)
  323. {
  324. if (v[0].GetSize() != v[i].GetSize())
  325. {
  326. LogError("Mismatched sizes.");
  327. }
  328. }
  329. Real minLength = Normalize(v[0], robust);
  330. for (int i = 1; i < numInputs; ++i)
  331. {
  332. for (int j = 0; j < i; ++j)
  333. {
  334. Real dot = Dot(v[i], v[j]);
  335. v[i] -= v[j] * dot;
  336. }
  337. Real length = Normalize(v[i], robust);
  338. if (length < minLength)
  339. {
  340. minLength = length;
  341. }
  342. }
  343. return minLength;
  344. }
  345. LogError("Invalid input.");
  346. }
  347. // Compute the axis-aligned bounding box of the vectors. The return value is
  348. // 'true' iff the inputs are valid, in which case vmin and vmax have valid
  349. // values.
  350. template <typename Real>
  351. bool ComputeExtremes(int numVectors, GVector<Real> const* v,
  352. GVector<Real>& vmin, GVector<Real>& vmax)
  353. {
  354. if (v && numVectors > 0)
  355. {
  356. for (int i = 1; i < numVectors; ++i)
  357. {
  358. if (v[0].GetSize() != v[i].GetSize())
  359. {
  360. LogError("Mismatched sizes.");
  361. }
  362. }
  363. int const size = v[0].GetSize();
  364. vmin = v[0];
  365. vmax = vmin;
  366. for (int j = 1; j < numVectors; ++j)
  367. {
  368. GVector<Real> const& vec = v[j];
  369. for (int i = 0; i < size; ++i)
  370. {
  371. if (vec[i] < vmin[i])
  372. {
  373. vmin[i] = vec[i];
  374. }
  375. else if (vec[i] > vmax[i])
  376. {
  377. vmax[i] = vec[i];
  378. }
  379. }
  380. }
  381. return true;
  382. }
  383. LogError("Invalid input.");
  384. }
  385. // Lift n-tuple v to homogeneous (n+1)-tuple (v,last).
  386. template <typename Real>
  387. GVector<Real> HLift(GVector<Real> const& v, Real last)
  388. {
  389. int const size = v.GetSize();
  390. GVector<Real> result(size + 1);
  391. for (int i = 0; i < size; ++i)
  392. {
  393. result[i] = v[i];
  394. }
  395. result[size] = last;
  396. return result;
  397. }
  398. // Project homogeneous n-tuple v = (u,v[n-1]) to (n-1)-tuple u.
  399. template <typename Real>
  400. GVector<Real> HProject(GVector<Real> const& v)
  401. {
  402. int const size = v.GetSize();
  403. if (size > 1)
  404. {
  405. GVector<Real> result(size - 1);
  406. for (int i = 0; i < size - 1; ++i)
  407. {
  408. result[i] = v[i];
  409. }
  410. return result;
  411. }
  412. else
  413. {
  414. return GVector<Real>();
  415. }
  416. }
  417. // Lift n-tuple v = (w0,w1) to (n+1)-tuple u = (w0,u[inject],w1). By
  418. // inference, w0 is a (inject)-tuple [nonexistent when inject=0] and w1
  419. // is a (n-inject)-tuple [nonexistent when inject=n].
  420. template <typename Real>
  421. GVector<Real> Lift(GVector<Real> const& v, int inject, Real value)
  422. {
  423. int const size = v.GetSize();
  424. GVector<Real> result(size + 1);
  425. int i;
  426. for (i = 0; i < inject; ++i)
  427. {
  428. result[i] = v[i];
  429. }
  430. result[i] = value;
  431. int j = i;
  432. for (++j; i < size; ++i, ++j)
  433. {
  434. result[j] = v[i];
  435. }
  436. return result;
  437. }
  438. // Project n-tuple v = (w0,v[reject],w1) to (n-1)-tuple u = (w0,w1). By
  439. // inference, w0 is a (reject)-tuple [nonexistent when reject=0] and w1
  440. // is a (n-1-reject)-tuple [nonexistent when reject=n-1].
  441. template <typename Real>
  442. GVector<Real> Project(GVector<Real> const& v, int reject)
  443. {
  444. int const size = v.GetSize();
  445. if (size > 1)
  446. {
  447. GVector<Real> result(size - 1);
  448. for (int i = 0, j = 0; i < size - 1; ++i, ++j)
  449. {
  450. if (j == reject)
  451. {
  452. ++j;
  453. }
  454. result[i] = v[j];
  455. }
  456. return result;
  457. }
  458. else
  459. {
  460. return GVector<Real>();
  461. }
  462. }
  463. }