1 /*
2  * Copyright (c) 2011-2021, The DART development contributors
3  * All rights reserved.
4  *
5  * The list of contributors can be found at:
6  *   https://github.com/dartsim/dart/blob/master/LICENSE
7  *
8  * This file is provided under the following "BSD-style" License:
9  *   Redistribution and use in source and binary forms, with or
10  *   without modification, are permitted provided that the following
11  *   conditions are met:
12  *   * Redistributions of source code must retain the above copyright
13  *     notice, this list of conditions and the following disclaimer.
14  *   * Redistributions in binary form must reproduce the above
15  *     copyright notice, this list of conditions and the following
16  *     disclaimer in the documentation and/or other materials provided
17  *     with the distribution.
18  *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
19  *   CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
20  *   INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21  *   MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
22  *   DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
23  *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
24  *   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
25  *   LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
26  *   USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
27  *   AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
28  *   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
29  *   ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
30  *   POSSIBILITY OF SUCH DAMAGE.
31  */
32 
33 #ifndef DART_MATH_HELPERS_HPP_
34 #define DART_MATH_HELPERS_HPP_
35 
36 // Standard Libraries
37 #include <cfloat>
38 #include <climits>
39 #include <cmath>
40 #include <cstdlib>
41 #include <ctime>
42 #include <iomanip>
43 #include <iostream>
44 #include <random>
45 
46 // External Libraries
47 #include <Eigen/Dense>
48 // Local Headers
49 #include "dart/math/Constants.hpp"
50 #include "dart/math/MathTypes.hpp"
51 #include "dart/math/Random.hpp"
52 
53 namespace dart {
54 namespace math {
55 
56 //==============================================================================
57 template <typename T>
toRadian(const T & degree)58 constexpr T toRadian(const T& degree)
59 {
60   return degree * constants<T>::pi() / 180.0;
61 }
62 
63 //==============================================================================
64 template <typename T>
toDegree(const T & radian)65 constexpr T toDegree(const T& radian)
66 {
67   return radian * 180.0 / constants<T>::pi();
68 }
69 
70 /// \brief a cross b = (CR*a) dot b
71 /// const Matd CR(2,2,0.0,-1.0,1.0,0.0);
72 const Eigen::Matrix2d CR((Eigen::Matrix2d() << 0.0, -1.0, 1.0, 0.0).finished());
73 
delta(int _i,int _j)74 inline int delta(int _i, int _j)
75 {
76   if (_i == _j)
77     return 1;
78   return 0;
79 }
80 
81 template <typename T>
sign(T x,std::false_type)82 inline constexpr int sign(T x, std::false_type)
83 {
84   return static_cast<T>(0) < x;
85 }
86 
87 template <typename T>
sign(T x,std::true_type)88 inline constexpr int sign(T x, std::true_type)
89 {
90   return (static_cast<T>(0) < x) - (x < static_cast<T>(0));
91 }
92 
93 template <typename T>
sign(T x)94 inline constexpr int sign(T x)
95 {
96   return sign(x, std::is_signed<T>());
97 }
98 
sqr(double _x)99 inline double sqr(double _x)
100 {
101   return _x * _x;
102 }
103 
Tsinc(double _theta)104 inline double Tsinc(double _theta)
105 {
106   return 0.5 - sqrt(_theta) / 48;
107 }
108 
isZero(double _theta)109 inline bool isZero(double _theta)
110 {
111   return (std::abs(_theta) < 1e-6);
112 }
113 
asinh(double _X)114 inline double asinh(double _X)
115 {
116   return log(_X + sqrt(_X * _X + 1));
117 }
118 
acosh(double _X)119 inline double acosh(double _X)
120 {
121   return log(_X + sqrt(_X * _X - 1));
122 }
123 
atanh(double _X)124 inline double atanh(double _X)
125 {
126   return log((1 + _X) / (1 - _X)) / 2;
127 }
128 
asech(double _X)129 inline double asech(double _X)
130 {
131   return log((sqrt(-_X * _X + 1) + 1) / _X);
132 }
133 
acosech(double _X)134 inline double acosech(double _X)
135 {
136   return log((sign(_X) * sqrt(_X * _X + 1) + 1) / _X);
137 }
138 
acotanh(double _X)139 inline double acotanh(double _X)
140 {
141   return log((_X + 1) / (_X - 1)) / 2;
142 }
143 
round(double _x)144 inline double round(double _x)
145 {
146   return floor(_x + 0.5);
147 }
148 
round2(double _x)149 inline double round2(double _x)
150 {
151   int gintx = static_cast<int>(std::floor(_x));
152   if (_x - gintx < 0.5)
153     return static_cast<double>(gintx);
154   else
155     return static_cast<double>(gintx + 1.0);
156 }
157 
158 template <typename T>
clip(const T & val,const T & lower,const T & upper)159 inline T clip(const T& val, const T& lower, const T& upper)
160 {
161   return std::max(lower, std::min(val, upper));
162 }
163 
164 template <typename DerivedA, typename DerivedB>
clip(const Eigen::MatrixBase<DerivedA> & val,const Eigen::MatrixBase<DerivedB> & lower,const Eigen::MatrixBase<DerivedB> & upper)165 inline typename DerivedA::PlainObject clip(
166     const Eigen::MatrixBase<DerivedA>& val,
167     const Eigen::MatrixBase<DerivedB>& lower,
168     const Eigen::MatrixBase<DerivedB>& upper)
169 {
170   return lower.cwiseMax(val.cwiseMin(upper));
171 }
172 
isEqual(double _x,double _y)173 inline bool isEqual(double _x, double _y)
174 {
175   return (std::abs(_x - _y) < 1e-6);
176 }
177 
178 // check if it is an integer
isInt(double _x)179 inline bool isInt(double _x)
180 {
181   if (isEqual(round(_x), _x))
182     return true;
183   return false;
184 }
185 
186 /// \brief Returns whether _v is a NaN (Not-A-Number) value
isNan(double _v)187 inline bool isNan(double _v)
188 {
189 #ifdef _WIN32
190   return _isnan(_v) != 0;
191 #else
192   return std::isnan(_v);
193 #endif
194 }
195 
196 /// \brief Returns whether _m is a NaN (Not-A-Number) matrix
isNan(const Eigen::MatrixXd & _m)197 inline bool isNan(const Eigen::MatrixXd& _m)
198 {
199   for (int i = 0; i < _m.rows(); ++i)
200     for (int j = 0; j < _m.cols(); ++j)
201       if (isNan(_m(i, j)))
202         return true;
203 
204   return false;
205 }
206 
207 /// \brief Returns whether _v is an infinity value (either positive infinity or
208 /// negative infinity).
isInf(double _v)209 inline bool isInf(double _v)
210 {
211 #ifdef _WIN32
212   return !_finite(_v);
213 #else
214   return std::isinf(_v);
215 #endif
216 }
217 
218 /// \brief Returns whether _m is an infinity matrix (either positive infinity or
219 /// negative infinity).
isInf(const Eigen::MatrixXd & _m)220 inline bool isInf(const Eigen::MatrixXd& _m)
221 {
222   for (int i = 0; i < _m.rows(); ++i)
223     for (int j = 0; j < _m.cols(); ++j)
224       if (isInf(_m(i, j)))
225         return true;
226 
227   return false;
228 }
229 
230 /// \brief Returns whether _m is symmetric or not
isSymmetric(const Eigen::MatrixXd & _m,double _tol=1e-6)231 inline bool isSymmetric(const Eigen::MatrixXd& _m, double _tol = 1e-6)
232 {
233   std::size_t rows = _m.rows();
234   std::size_t cols = _m.cols();
235 
236   if (rows != cols)
237     return false;
238 
239   for (std::size_t i = 0; i < rows; ++i)
240   {
241     for (std::size_t j = i + 1; j < cols; ++j)
242     {
243       if (std::abs(_m(i, j) - _m(j, i)) > _tol)
244       {
245         std::cout << "A: " << std::endl;
246         for (std::size_t k = 0; k < rows; ++k)
247         {
248           for (std::size_t l = 0; l < cols; ++l)
249             std::cout << std::setprecision(4) << _m(k, l) << " ";
250           std::cout << std::endl;
251         }
252 
253         std::cout << "A(" << i << ", " << j << "): " << _m(i, j) << std::endl;
254         std::cout << "A(" << j << ", " << i << "): " << _m(i, j) << std::endl;
255         return false;
256       }
257     }
258   }
259 
260   return true;
261 }
262 
seedRand()263 inline unsigned seedRand()
264 {
265   time_t now = time(0);
266   unsigned char* p = reinterpret_cast<unsigned char*>(&now);
267   unsigned seed = 0;
268   std::size_t i;
269 
270   for (i = 0; i < sizeof(now); i++)
271     seed = seed * (UCHAR_MAX + 2U) + p[i];
272 
273   srand(seed);
274   return seed;
275 }
276 
277 /// \deprecated Please use Random::uniform() instead.
278 DART_DEPRECATED(6.7)
random(double _min,double _max)279 inline double random(double _min, double _max)
280 {
281   return _min
282          + ((static_cast<double>(rand()) / (RAND_MAX + 1.0)) * (_max - _min));
283 }
284 
285 /// \deprecated Please use Random::uniform() instead.
286 template <int N>
287 DART_DEPRECATED(6.7)
randomVector(double _min,double _max)288 Eigen::Matrix<double, N, 1> randomVector(double _min, double _max)
289 {
290   Eigen::Matrix<double, N, 1> v;
291   DART_SUPPRESS_DEPRECATED_BEGIN
292   for (std::size_t i = 0; i < N; ++i)
293     v[i] = random(_min, _max);
294   DART_SUPPRESS_DEPRECATED_END
295 
296   return v;
297 }
298 
299 /// \deprecated Please use Random::uniform() instead.
300 template <int N>
301 DART_DEPRECATED(6.7)
randomVector(double _limit)302 Eigen::Matrix<double, N, 1> randomVector(double _limit)
303 {
304   DART_SUPPRESS_DEPRECATED_BEGIN
305   return randomVector<N>(-std::abs(_limit), std::abs(_limit));
306   DART_SUPPRESS_DEPRECATED_END
307 }
308 
309 //==============================================================================
310 /// \deprecated Please use Random::uniform() instead.
311 DART_DEPRECATED(6.7)
randomVectorXd(std::size_t size,double min,double max)312 inline Eigen::VectorXd randomVectorXd(std::size_t size, double min, double max)
313 {
314   Eigen::VectorXd v = Eigen::VectorXd::Zero(size);
315 
316   DART_SUPPRESS_DEPRECATED_BEGIN
317   for (std::size_t i = 0; i < size; ++i)
318     v[i] = random(min, max);
319   DART_SUPPRESS_DEPRECATED_END
320 
321   return v;
322 }
323 
324 //==============================================================================
325 /// \deprecated Please use Random::uniform() instead.
326 DART_DEPRECATED(6.7)
randomVectorXd(std::size_t size,double limit)327 inline Eigen::VectorXd randomVectorXd(std::size_t size, double limit)
328 {
329   DART_SUPPRESS_DEPRECATED_BEGIN
330   return randomVectorXd(size, -std::abs(limit), std::abs(limit));
331   DART_SUPPRESS_DEPRECATED_END
332 }
333 
334 namespace suffixes {
335 
336 //==============================================================================
operator ""_pi(long double x)337 constexpr double operator"" _pi(long double x)
338 {
339   return x * constants<double>::pi();
340 }
341 
342 //==============================================================================
operator ""_pi(unsigned long long int x)343 constexpr double operator"" _pi(unsigned long long int x)
344 {
345   return operator"" _pi(static_cast<long double>(x));
346 }
347 
348 //==============================================================================
operator ""_rad(long double angle)349 constexpr double operator"" _rad(long double angle)
350 {
351   return angle;
352 }
353 
354 //==============================================================================
operator ""_rad(unsigned long long int angle)355 constexpr double operator"" _rad(unsigned long long int angle)
356 {
357   return operator"" _rad(static_cast<long double>(angle));
358 }
359 
360 //==============================================================================
operator ""_deg(long double angle)361 constexpr double operator"" _deg(long double angle)
362 {
363   return toRadian(angle);
364 }
365 
366 //==============================================================================
operator ""_deg(unsigned long long int angle)367 constexpr double operator"" _deg(unsigned long long int angle)
368 {
369   return operator"" _deg(static_cast<long double>(angle));
370 }
371 
372 } // namespace suffixes
373 
374 } // namespace math
375 
376 namespace Color {
377 
Red(double alpha)378 inline Eigen::Vector4d Red(double alpha)
379 {
380   return Eigen::Vector4d(0.9, 0.1, 0.1, alpha);
381 }
382 
Red()383 inline Eigen::Vector3d Red()
384 {
385   return Eigen::Vector3d(0.9, 0.1, 0.1);
386 }
387 
Fuchsia()388 inline Eigen::Vector3d Fuchsia()
389 {
390   return Eigen::Vector3d(1.0, 0.0, 0.5);
391 }
392 
Fuchsia(double alpha)393 inline Eigen::Vector4d Fuchsia(double alpha)
394 {
395   return Eigen::Vector4d(1.0, 0.0, 0.5, alpha);
396 }
397 
Orange(double alpha)398 inline Eigen::Vector4d Orange(double alpha)
399 {
400   return Eigen::Vector4d(1.0, 0.63, 0.0, alpha);
401 }
402 
Orange()403 inline Eigen::Vector3d Orange()
404 {
405   return Eigen::Vector3d(1.0, 0.63, 0.0);
406 }
407 
Green(double alpha)408 inline Eigen::Vector4d Green(double alpha)
409 {
410   return Eigen::Vector4d(0.1, 0.9, 0.1, alpha);
411 }
412 
Green()413 inline Eigen::Vector3d Green()
414 {
415   return Eigen::Vector3d(0.1, 0.9, 0.1);
416 }
417 
Blue(double alpha)418 inline Eigen::Vector4d Blue(double alpha)
419 {
420   return Eigen::Vector4d(0.1, 0.1, 0.9, alpha);
421 }
422 
Blue()423 inline Eigen::Vector3d Blue()
424 {
425   return Eigen::Vector3d(0.1, 0.1, 0.9);
426 }
427 
White(double alpha)428 inline Eigen::Vector4d White(double alpha)
429 {
430   return Eigen::Vector4d(1.0, 1.0, 1.0, alpha);
431 }
432 
White()433 inline Eigen::Vector3d White()
434 {
435   return Eigen::Vector3d(1.0, 1.0, 1.0);
436 }
437 
Black(double alpha)438 inline Eigen::Vector4d Black(double alpha)
439 {
440   return Eigen::Vector4d(0.05, 0.05, 0.05, alpha);
441 }
442 
Black()443 inline Eigen::Vector3d Black()
444 {
445   return Eigen::Vector3d(0.05, 0.05, 0.05);
446 }
447 
LightGray(double alpha)448 inline Eigen::Vector4d LightGray(double alpha)
449 {
450   return Eigen::Vector4d(0.9, 0.9, 0.9, alpha);
451 }
452 
LightGray()453 inline Eigen::Vector3d LightGray()
454 {
455   return Eigen::Vector3d(0.9, 0.9, 0.9);
456 }
457 
Gray(double alpha)458 inline Eigen::Vector4d Gray(double alpha)
459 {
460   return Eigen::Vector4d(0.6, 0.6, 0.6, alpha);
461 }
462 
Gray()463 inline Eigen::Vector3d Gray()
464 {
465   return Eigen::Vector3d(0.6, 0.6, 0.6);
466 }
467 
Random(double alpha)468 inline Eigen::Vector4d Random(double alpha)
469 {
470   return Eigen::Vector4d(
471       math::Random::uniform(0.0, 1.0),
472       math::Random::uniform(0.0, 1.0),
473       math::Random::uniform(0.0, 1.0),
474       alpha);
475 }
476 
Random()477 inline Eigen::Vector3d Random()
478 {
479   return Eigen::Vector3d(
480       math::Random::uniform(0.0, 1.0),
481       math::Random::uniform(0.0, 1.0),
482       math::Random::uniform(0.0, 1.0));
483 }
484 
485 } // namespace Color
486 
487 } // namespace dart
488 
489 #endif // DART_MATH_HELPERS_HPP_
490