ClangFormat: apply to source, most of intern
[blender.git] / intern / cycles / util / util_math.h
1 /*
2  * Copyright 2011-2013 Blender Foundation
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16
17 #ifndef __UTIL_MATH_H__
18 #define __UTIL_MATH_H__
19
20 /* Math
21  *
22  * Basic math functions on scalar and vector types. This header is used by
23  * both the kernel code when compiled as C++, and other C++ non-kernel code. */
24
25 #ifndef __KERNEL_GPU__
26 #  include <cmath>
27 #endif
28
29 #ifndef __KERNEL_OPENCL__
30 #  include <float.h>
31 #  include <math.h>
32 #  include <stdio.h>
33 #endif /* __KERNEL_OPENCL__ */
34
35 #include "util/util_types.h"
36
37 CCL_NAMESPACE_BEGIN
38
39 /* Float Pi variations */
40
41 /* Division */
42 #ifndef M_PI_F
43 #  define M_PI_F (3.1415926535897932f) /* pi */
44 #endif
45 #ifndef M_PI_2_F
46 #  define M_PI_2_F (1.5707963267948966f) /* pi/2 */
47 #endif
48 #ifndef M_PI_4_F
49 #  define M_PI_4_F (0.7853981633974830f) /* pi/4 */
50 #endif
51 #ifndef M_1_PI_F
52 #  define M_1_PI_F (0.3183098861837067f) /* 1/pi */
53 #endif
54 #ifndef M_2_PI_F
55 #  define M_2_PI_F (0.6366197723675813f) /* 2/pi */
56 #endif
57 #ifndef M_1_2PI_F
58 #  define M_1_2PI_F (0.1591549430918953f) /* 1/(2*pi) */
59 #endif
60 #ifndef M_SQRT_PI_8_F
61 #  define M_SQRT_PI_8_F (0.6266570686577501f) /* sqrt(pi/8) */
62 #endif
63 #ifndef M_LN_2PI_F
64 #  define M_LN_2PI_F (1.8378770664093454f) /* ln(2*pi) */
65 #endif
66
67 /* Multiplication */
68 #ifndef M_2PI_F
69 #  define M_2PI_F (6.2831853071795864f) /* 2*pi */
70 #endif
71 #ifndef M_4PI_F
72 #  define M_4PI_F (12.566370614359172f) /* 4*pi */
73 #endif
74
75 /* Float sqrt variations */
76 #ifndef M_SQRT2_F
77 #  define M_SQRT2_F (1.4142135623730950f) /* sqrt(2) */
78 #endif
79 #ifndef M_LN2_F
80 #  define M_LN2_F (0.6931471805599453f) /* ln(2) */
81 #endif
82 #ifndef M_LN10_F
83 #  define M_LN10_F (2.3025850929940457f) /* ln(10) */
84 #endif
85
86 /* Scalar */
87
88 #ifdef _WIN32
89 #  ifndef __KERNEL_OPENCL__
90 ccl_device_inline float fmaxf(float a, float b)
91 {
92   return (a > b) ? a : b;
93 }
94
95 ccl_device_inline float fminf(float a, float b)
96 {
97   return (a < b) ? a : b;
98 }
99 #  endif /* !__KERNEL_OPENCL__ */
100 #endif   /* _WIN32 */
101
102 #ifndef __KERNEL_GPU__
103 using std::isfinite;
104 using std::isnan;
105 using std::sqrt;
106
107 ccl_device_inline int abs(int x)
108 {
109   return (x > 0) ? x : -x;
110 }
111
112 ccl_device_inline int max(int a, int b)
113 {
114   return (a > b) ? a : b;
115 }
116
117 ccl_device_inline int min(int a, int b)
118 {
119   return (a < b) ? a : b;
120 }
121
122 ccl_device_inline float max(float a, float b)
123 {
124   return (a > b) ? a : b;
125 }
126
127 ccl_device_inline float min(float a, float b)
128 {
129   return (a < b) ? a : b;
130 }
131
132 ccl_device_inline double max(double a, double b)
133 {
134   return (a > b) ? a : b;
135 }
136
137 ccl_device_inline double min(double a, double b)
138 {
139   return (a < b) ? a : b;
140 }
141
142 /* These 2 guys are templated for usage with registers data.
143  *
144  * NOTE: Since this is CPU-only functions it is ok to use references here.
145  * But for other devices we'll need to be careful about this.
146  */
147
148 template<typename T> ccl_device_inline T min4(const T &a, const T &b, const T &c, const T &d)
149 {
150   return min(min(a, b), min(c, d));
151 }
152
153 template<typename T> ccl_device_inline T max4(const T &a, const T &b, const T &c, const T &d)
154 {
155   return max(max(a, b), max(c, d));
156 }
157 #endif /* __KERNEL_GPU__ */
158
159 ccl_device_inline float min4(float a, float b, float c, float d)
160 {
161   return min(min(a, b), min(c, d));
162 }
163
164 ccl_device_inline float max4(float a, float b, float c, float d)
165 {
166   return max(max(a, b), max(c, d));
167 }
168
169 #ifndef __KERNEL_OPENCL__
170 /* Int/Float conversion */
171
172 ccl_device_inline int as_int(uint i)
173 {
174   union {
175     uint ui;
176     int i;
177   } u;
178   u.ui = i;
179   return u.i;
180 }
181
182 ccl_device_inline uint as_uint(int i)
183 {
184   union {
185     uint ui;
186     int i;
187   } u;
188   u.i = i;
189   return u.ui;
190 }
191
192 ccl_device_inline uint as_uint(float f)
193 {
194   union {
195     uint i;
196     float f;
197   } u;
198   u.f = f;
199   return u.i;
200 }
201
202 ccl_device_inline int __float_as_int(float f)
203 {
204   union {
205     int i;
206     float f;
207   } u;
208   u.f = f;
209   return u.i;
210 }
211
212 ccl_device_inline float __int_as_float(int i)
213 {
214   union {
215     int i;
216     float f;
217   } u;
218   u.i = i;
219   return u.f;
220 }
221
222 ccl_device_inline uint __float_as_uint(float f)
223 {
224   union {
225     uint i;
226     float f;
227   } u;
228   u.f = f;
229   return u.i;
230 }
231
232 ccl_device_inline float __uint_as_float(uint i)
233 {
234   union {
235     uint i;
236     float f;
237   } u;
238   u.i = i;
239   return u.f;
240 }
241
242 ccl_device_inline int4 __float4_as_int4(float4 f)
243 {
244 #  ifdef __KERNEL_SSE__
245   return int4(_mm_castps_si128(f.m128));
246 #  else
247   return make_int4(
248       __float_as_int(f.x), __float_as_int(f.y), __float_as_int(f.z), __float_as_int(f.w));
249 #  endif
250 }
251
252 ccl_device_inline float4 __int4_as_float4(int4 i)
253 {
254 #  ifdef __KERNEL_SSE__
255   return float4(_mm_castsi128_ps(i.m128));
256 #  else
257   return make_float4(
258       __int_as_float(i.x), __int_as_float(i.y), __int_as_float(i.z), __int_as_float(i.w));
259 #  endif
260 }
261 #endif /* __KERNEL_OPENCL__ */
262
263 /* Versions of functions which are safe for fast math. */
264 ccl_device_inline bool isnan_safe(float f)
265 {
266   unsigned int x = __float_as_uint(f);
267   return (x << 1) > 0xff000000u;
268 }
269
270 ccl_device_inline bool isfinite_safe(float f)
271 {
272   /* By IEEE 754 rule, 2*Inf equals Inf */
273   unsigned int x = __float_as_uint(f);
274   return (f == f) && (x == 0 || x == (1u << 31) || (f != 2.0f * f)) && !((x << 1) > 0xff000000u);
275 }
276
277 ccl_device_inline float ensure_finite(float v)
278 {
279   return isfinite_safe(v) ? v : 0.0f;
280 }
281
282 #ifndef __KERNEL_OPENCL__
283 ccl_device_inline int clamp(int a, int mn, int mx)
284 {
285   return min(max(a, mn), mx);
286 }
287
288 ccl_device_inline float clamp(float a, float mn, float mx)
289 {
290   return min(max(a, mn), mx);
291 }
292
293 ccl_device_inline float mix(float a, float b, float t)
294 {
295   return a + t * (b - a);
296 }
297 #endif /* __KERNEL_OPENCL__ */
298
299 #ifndef __KERNEL_CUDA__
300 ccl_device_inline float saturate(float a)
301 {
302   return clamp(a, 0.0f, 1.0f);
303 }
304 #endif /* __KERNEL_CUDA__ */
305
306 ccl_device_inline int float_to_int(float f)
307 {
308   return (int)f;
309 }
310
311 ccl_device_inline int floor_to_int(float f)
312 {
313   return float_to_int(floorf(f));
314 }
315
316 ccl_device_inline int quick_floor_to_int(float x)
317 {
318   return float_to_int(x) - ((x < 0) ? 1 : 0);
319 }
320
321 ccl_device_inline int ceil_to_int(float f)
322 {
323   return float_to_int(ceilf(f));
324 }
325
326 ccl_device_inline float signf(float f)
327 {
328   return (f < 0.0f) ? -1.0f : 1.0f;
329 }
330
331 ccl_device_inline float nonzerof(float f, float eps)
332 {
333   if (fabsf(f) < eps)
334     return signf(f) * eps;
335   else
336     return f;
337 }
338
339 ccl_device_inline float smoothstepf(float f)
340 {
341   float ff = f * f;
342   return (3.0f * ff - 2.0f * ff * f);
343 }
344
345 ccl_device_inline int mod(int x, int m)
346 {
347   return (x % m + m) % m;
348 }
349
350 ccl_device_inline float3 float2_to_float3(const float2 a)
351 {
352   return make_float3(a.x, a.y, 0.0f);
353 }
354
355 ccl_device_inline float3 float4_to_float3(const float4 a)
356 {
357   return make_float3(a.x, a.y, a.z);
358 }
359
360 ccl_device_inline float4 float3_to_float4(const float3 a)
361 {
362   return make_float4(a.x, a.y, a.z, 1.0f);
363 }
364
365 ccl_device_inline float inverse_lerp(float a, float b, float x)
366 {
367   return (x - a) / (b - a);
368 }
369
370 /* Cubic interpolation between b and c, a and d are the previous and next point. */
371 ccl_device_inline float cubic_interp(float a, float b, float c, float d, float x)
372 {
373   return 0.5f *
374              (((d + 3.0f * (b - c) - a) * x + (2.0f * a - 5.0f * b + 4.0f * c - d)) * x +
375               (c - a)) *
376              x +
377          b;
378 }
379
380 CCL_NAMESPACE_END
381
382 #include "util/util_math_int2.h"
383 #include "util/util_math_int3.h"
384 #include "util/util_math_int4.h"
385
386 #include "util/util_math_float2.h"
387 #include "util/util_math_float3.h"
388 #include "util/util_math_float4.h"
389
390 #include "util/util_rect.h"
391
392 CCL_NAMESPACE_BEGIN
393
394 #ifndef __KERNEL_OPENCL__
395 /* Interpolation */
396
397 template<class A, class B> A lerp(const A &a, const A &b, const B &t)
398 {
399   return (A)(a * ((B)1 - t) + b * t);
400 }
401
402 #endif /* __KERNEL_OPENCL__ */
403
404 /* Triangle */
405
406 #ifndef __KERNEL_OPENCL__
407 ccl_device_inline float triangle_area(const float3 &v1, const float3 &v2, const float3 &v3)
408 #else
409 ccl_device_inline float triangle_area(const float3 v1, const float3 v2, const float3 v3)
410 #endif
411 {
412   return len(cross(v3 - v2, v1 - v2)) * 0.5f;
413 }
414
415 /* Orthonormal vectors */
416
417 ccl_device_inline void make_orthonormals(const float3 N, float3 *a, float3 *b)
418 {
419 #if 0
420   if(fabsf(N.y) >= 0.999f) {
421     *a = make_float3(1, 0, 0);
422     *b = make_float3(0, 0, 1);
423     return;
424   }
425   if(fabsf(N.z) >= 0.999f) {
426     *a = make_float3(1, 0, 0);
427     *b = make_float3(0, 1, 0);
428     return;
429   }
430 #endif
431
432   if (N.x != N.y || N.x != N.z)
433     *a = make_float3(N.z - N.y, N.x - N.z, N.y - N.x);  //(1,1,1)x N
434   else
435     *a = make_float3(N.z - N.y, N.x + N.z, -N.y - N.x);  //(-1,1,1)x N
436
437   *a = normalize(*a);
438   *b = cross(N, *a);
439 }
440
441 /* Color division */
442
443 ccl_device_inline float3 safe_invert_color(float3 a)
444 {
445   float x, y, z;
446
447   x = (a.x != 0.0f) ? 1.0f / a.x : 0.0f;
448   y = (a.y != 0.0f) ? 1.0f / a.y : 0.0f;
449   z = (a.z != 0.0f) ? 1.0f / a.z : 0.0f;
450
451   return make_float3(x, y, z);
452 }
453
454 ccl_device_inline float3 safe_divide_color(float3 a, float3 b)
455 {
456   float x, y, z;
457
458   x = (b.x != 0.0f) ? a.x / b.x : 0.0f;
459   y = (b.y != 0.0f) ? a.y / b.y : 0.0f;
460   z = (b.z != 0.0f) ? a.z / b.z : 0.0f;
461
462   return make_float3(x, y, z);
463 }
464
465 ccl_device_inline float3 safe_divide_even_color(float3 a, float3 b)
466 {
467   float x, y, z;
468
469   x = (b.x != 0.0f) ? a.x / b.x : 0.0f;
470   y = (b.y != 0.0f) ? a.y / b.y : 0.0f;
471   z = (b.z != 0.0f) ? a.z / b.z : 0.0f;
472
473   /* try to get gray even if b is zero */
474   if (b.x == 0.0f) {
475     if (b.y == 0.0f) {
476       x = z;
477       y = z;
478     }
479     else if (b.z == 0.0f) {
480       x = y;
481       z = y;
482     }
483     else
484       x = 0.5f * (y + z);
485   }
486   else if (b.y == 0.0f) {
487     if (b.z == 0.0f) {
488       y = x;
489       z = x;
490     }
491     else
492       y = 0.5f * (x + z);
493   }
494   else if (b.z == 0.0f) {
495     z = 0.5f * (x + y);
496   }
497
498   return make_float3(x, y, z);
499 }
500
501 /* Rotation of point around axis and angle */
502
503 ccl_device_inline float3 rotate_around_axis(float3 p, float3 axis, float angle)
504 {
505   float costheta = cosf(angle);
506   float sintheta = sinf(angle);
507   float3 r;
508
509   r.x = ((costheta + (1 - costheta) * axis.x * axis.x) * p.x) +
510         (((1 - costheta) * axis.x * axis.y - axis.z * sintheta) * p.y) +
511         (((1 - costheta) * axis.x * axis.z + axis.y * sintheta) * p.z);
512
513   r.y = (((1 - costheta) * axis.x * axis.y + axis.z * sintheta) * p.x) +
514         ((costheta + (1 - costheta) * axis.y * axis.y) * p.y) +
515         (((1 - costheta) * axis.y * axis.z - axis.x * sintheta) * p.z);
516
517   r.z = (((1 - costheta) * axis.x * axis.z - axis.y * sintheta) * p.x) +
518         (((1 - costheta) * axis.y * axis.z + axis.x * sintheta) * p.y) +
519         ((costheta + (1 - costheta) * axis.z * axis.z) * p.z);
520
521   return r;
522 }
523
524 /* NaN-safe math ops */
525
526 ccl_device_inline float safe_sqrtf(float f)
527 {
528   return sqrtf(max(f, 0.0f));
529 }
530
531 ccl_device float safe_asinf(float a)
532 {
533   return asinf(clamp(a, -1.0f, 1.0f));
534 }
535
536 ccl_device float safe_acosf(float a)
537 {
538   return acosf(clamp(a, -1.0f, 1.0f));
539 }
540
541 ccl_device float compatible_powf(float x, float y)
542 {
543 #ifdef __KERNEL_GPU__
544   if (y == 0.0f) /* x^0 -> 1, including 0^0 */
545     return 1.0f;
546
547   /* GPU pow doesn't accept negative x, do manual checks here */
548   if (x < 0.0f) {
549     if (fmodf(-y, 2.0f) == 0.0f)
550       return powf(-x, y);
551     else
552       return -powf(-x, y);
553   }
554   else if (x == 0.0f)
555     return 0.0f;
556 #endif
557   return powf(x, y);
558 }
559
560 ccl_device float safe_powf(float a, float b)
561 {
562   if (UNLIKELY(a < 0.0f && b != float_to_int(b)))
563     return 0.0f;
564
565   return compatible_powf(a, b);
566 }
567
568 ccl_device float safe_divide(float a, float b)
569 {
570   return (b != 0.0f) ? a / b : 0.0f;
571 }
572
573 ccl_device float safe_logf(float a, float b)
574 {
575   if (UNLIKELY(a <= 0.0f || b <= 0.0f))
576     return 0.0f;
577
578   return safe_divide(logf(a), logf(b));
579 }
580
581 ccl_device float safe_modulo(float a, float b)
582 {
583   return (b != 0.0f) ? fmodf(a, b) : 0.0f;
584 }
585
586 ccl_device_inline float sqr(float a)
587 {
588   return a * a;
589 }
590
591 ccl_device_inline float pow20(float a)
592 {
593   return sqr(sqr(sqr(sqr(a)) * a));
594 }
595
596 ccl_device_inline float pow22(float a)
597 {
598   return sqr(a * sqr(sqr(sqr(a)) * a));
599 }
600
601 ccl_device_inline float beta(float x, float y)
602 {
603 #ifndef __KERNEL_OPENCL__
604   return expf(lgammaf(x) + lgammaf(y) - lgammaf(x + y));
605 #else
606   return expf(lgamma(x) + lgamma(y) - lgamma(x + y));
607 #endif
608 }
609
610 ccl_device_inline float xor_signmask(float x, int y)
611 {
612   return __int_as_float(__float_as_int(x) ^ y);
613 }
614
615 ccl_device float bits_to_01(uint bits)
616 {
617   return bits * (1.0f / (float)0xFFFFFFFF);
618 }
619
620 /* projections */
621 ccl_device_inline float2 map_to_tube(const float3 co)
622 {
623   float len, u, v;
624   len = sqrtf(co.x * co.x + co.y * co.y);
625   if (len > 0.0f) {
626     u = (1.0f - (atan2f(co.x / len, co.y / len) / M_PI_F)) * 0.5f;
627     v = (co.z + 1.0f) * 0.5f;
628   }
629   else {
630     u = v = 0.0f;
631   }
632   return make_float2(u, v);
633 }
634
635 ccl_device_inline float2 map_to_sphere(const float3 co)
636 {
637   float l = len(co);
638   float u, v;
639   if (l > 0.0f) {
640     if (UNLIKELY(co.x == 0.0f && co.y == 0.0f)) {
641       u = 0.0f; /* othwise domain error */
642     }
643     else {
644       u = (1.0f - atan2f(co.x, co.y) / M_PI_F) / 2.0f;
645     }
646     v = 1.0f - safe_acosf(co.z / l) / M_PI_F;
647   }
648   else {
649     u = v = 0.0f;
650   }
651   return make_float2(u, v);
652 }
653
654 CCL_NAMESPACE_END
655
656 #endif /* __UTIL_MATH_H__ */