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