cgv
Loading...
Searching...
No Matches
functions.h
1#pragma once
2
3#include <cmath>
4#include <cassert>
5#include <limits>
6#include <vector>
7
8#include "constants.h"
9
10#include "lib_begin.h"
11
12namespace cgv {
13namespace math {
14
16template<typename T>
17T sgn(const T& v) { return (T(0) < v) - (v < T(0)); }
19template<typename T>
20T sign(const T&v) { return (v >= T(0)) ? T(1) : T(-1); }
22template<typename T>
23T sign(const T& a, const T& b) { return (b < 0) ? -std::abs(a) : std::abs(a); }
25template<typename T>
26T step(const T& a, const T& b) { return (b < a) ? T(0) : T(1); }
28template<typename T>
29T plus(const T& v) { return v >= T(0) ? v : T(0); }
31template<typename T>
32T clamp(const T& v, const T& a, const T& b) { return v > b ? b : (v < a ? a : v); }
34template<typename T>
35T saturate(const T& v) { return v > T(1) ? T(1) : (v < T(0) ? T(0) : v); }
37template <typename T>
38T lerp(const T& a, const T& b, T t) { return (1 - t) * a + t * b; }
40template<typename T>
41T normalize(const T& v, const T& a, const T& b) { return (v - a) / (b - a); }
43template<typename T>
44T map(const T& v, const T& a, const T& b, const T& c, const T& d) { return c + (d - c) * ((v - a) / (b - a)); }
45
46namespace detail {
47
49template <typename T>
50T erfccheb(const T z)
51{
52 static const int ncof = 28;
53 static const double cof[28] = {
54 -1.3026537197817094, 6.4196979235649026e-1,
55 1.9476473204185836e-2,-9.561514786808631e-3,-9.46595344482036e-4,
56 3.66839497852761e-4,4.2523324806907e-5,-2.0278578112534e-5,
57 -1.624290004647e-6,1.303655835580e-6,1.5626441722e-8,-8.5238095915e-8,
58 6.529054439e-9,5.059343495e-9,-9.91364156e-10,-2.27365122e-10,
59 9.6467911e-11, 2.394038e-12,-6.886027e-12,8.94487e-13, 3.13092e-13,
60 -1.12708e-13,3.81e-16,7.106e-15,-1.523e-15,-9.4e-17,1.21e-16,-2.8e-17
61 };
62
63 int j;
64 double t, ty, tmp, d = 0.0, dd = 0.0;
65 assert(z >= 0.0);
66
67 t = 2.0 / (2.0 + z);
68 ty = 4.0*t - 2.0;
69 for (j = ncof - 1; j > 0; j--)
70 {
71 tmp = d;
72 d = ty * d - dd + cof[j];
73 dd = tmp;
74 }
75 return static_cast<T>(t*exp(-z * z + 0.5*(cof[0] + ty * d) - dd));
76}
77
78} // namespace detail
79
81template <typename T>
82T erf(const T x)
83{
84 if (x >= 0.)
85 return T(1.0) - detail::erfccheb(x);
86 else
87 return detail::erfccheb(-x) - T(1.0);
88}
89
91template <typename T>
92T erfc(const T x)
93{
94 if (x >= 0.)
95 return detail::erfccheb(x);
96 else
97 return 2.0 - detail::erfccheb(-x);
98}
99
101template <typename T>
102T erfcx(const T x)
103{
104 return exp(x*x)* erfccheb(x);
105}
106
108template <typename T>
109T erfc_inv(const T p)
110{
111 double x, err, t, pp;
112 if (p >= 2.0) return -100.;
113 if (p <= 0.0) return 100.;
114 pp = (p < 1.0) ? p : 2. - p;
115 t = sqrt(-2.*log(pp / 2.));
116 x = -0.70711*((2.30753 + t * 0.27061) / (1. + t * (0.99229 + t * 0.04481)) - t);
117 for (int j = 0; j < 2; j++) {
118 err = erfc(x) - pp;
119 x += err / (1.12837916709551257*exp(-(x*x)) - x * err);
120 }
121 return (p < 1.0 ? x : -x);
122}
123
125template <typename T>
126T erf_inv(const T p)
127{
128 return erfc_inv(1. - p);
129}
130
131
133template <typename T>
134T Phi(const T& x)
135{
136 return (T)(0.5*erfc(-x*sqrt(0.5)));
137}
138
140template <typename T>
141T inv_Phi(const T& p)
142{
143 static T a[] = {(T)-3.969683028665376e+01,(T) 2.209460984245205e+02,(T) -2.759285104469687e+02,(T) 1.383577518672690e+02,(T) -3.066479806614716e+01,(T) 2.506628277459239e+00};
144 static T b[] = {(T)-5.447609879822406e+01,(T) 1.615858368580409e+02,(T) -1.556989798598866e+02,(T) 6.680131188771972e+01,(T) -1.328068155288572e+01};
145 static T c[] = {(T)-7.784894002430293e-03,(T)-3.223964580411365e-01,(T)-2.400758277161838e+00,(T)-2.549732539343734e+00,(T) 4.374664141464968e+00,(T) 2.938163982698783e+00};
146 static T d[] = {(T)7.784695709041462e-03,(T)3.224671290700398e-01,(T)2.445134137142996e+00,(T)3.754408661907416e+00};
147
148 if (p <= 0)
149 return -std::numeric_limits<T>::infinity();
150 if (p >= 1)
151 return std::numeric_limits<T>::infinity();
152 T x;
153 if (p < (T)0.02425) {
154 T q = sqrt(-(T)2.0*log(p));
155 x = (((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /
156 ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1);
157 }
158 else if (p <= (T)0.97575) {
159 T q = p - (T)0.5;
160 T r = q*q;
161 x = (((((a(1)*r+a(2))*r+a(3))*r+a(4))*r+a(5))*r+a(6))*q /
162 (((((b(1)*r+b(2))*r+b(3))*r+b(4))*r+b(5))*r+1);
163 }
164 else {
165 T q = sqrt(-(T)2.0*log(1-p));
166 x = -(((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /
167 ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1);
168 }
169 return x;
170}
172extern CGV_API double compute_unit_ball_volume(unsigned n);
174extern CGV_API double compute_ball_volume(unsigned n, double R);
176extern CGV_API double compute_unit_sphere_area(unsigned n);
178extern CGV_API double compute_sphere_area(unsigned n, double R);
179
180//returns v*v
181template <typename T>
182T sqr(const T& v)
183{
184 return v*v;
185}
186
187//return the maximum of a and b
188template <typename T>
189T maximum(const T& a, const T&b)
190{
191 return a > b ? a : b;
192}
193
194//return minimum of a and b
195template<typename T>
196T minimum(const T &a, const T&b)
197{
198 return a < b ? a : b;
199}
200
201//return the maximum of a, b and c
202template <typename T>
203T maximum(const T& a, const T&b, const T& c)
204{
205 if(a > b)
206 {
207 return a > c ? a : c;
208 }
209 else
210 {
211 return b > c ? b : c;
212 }
213}
214
215//return the minimum of a, b and c
216template <typename T>
217T minimum(const T& a, const T&b, const T& c)
218{
219 if(a < b)
220 {
221 return a < c ? a : c;
222 }
223 else
224 {
225 return b < c ? b : c;
226 }
227}
228
229//convert angles measured in degrees to angles measured in radians
230template <typename T>
231T deg2rad(T deg)
232{
233 return static_cast<T>(deg * static_cast<T>(cgv::math::constants::pi / 180.0));
234}
235
236//convert angles measured in radians to angles measured in degrees
237template <typename T>
238T rad2deg(T rad)
239{
240 return static_cast<T>(rad * static_cast<T>(180.0 * cgv::math::constants::inv_pi));
241}
242
243//ln(gamma(x))
244extern CGV_API double gamma_ln(const double xx);
245
246//factorial of n (n!)
247extern CGV_API double fac(const int n);
248
249//returns ln(n!)
250extern CGV_API double fac_ln(const int n);
251
252//beta function
253extern CGV_API double beta(const double z, const double w);
254
255//binomial coefficient n over k
256extern CGV_API double nchoosek(const int n, const int k);
257
258} // namespace math
259} // namespace cgv
260
261#include <cgv/config/lib_end.h>
cgv
this header is dependency free
Definition print.h:11