BLI_math_rotation: properly name the quaternion power function.
[blender.git] / source / blender / blenlib / intern / math_geom_inline.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: some of this file.
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  * */
25
26 /** \file blender/blenlib/intern/math_geom_inline.c
27  *  \ingroup bli
28  */
29
30 #ifndef __MATH_GEOM_INLINE_C__
31 #define __MATH_GEOM_INLINE_C__
32
33 #include "BLI_math.h"
34
35 #include <string.h>
36
37 /* A few small defines. Keep'em local! */
38 #define SMALL_NUMBER  1.e-8f
39
40 /********************************** Polygons *********************************/
41
42 MINLINE float cross_tri_v2(const float v1[2], const float v2[2], const float v3[2])
43 {
44         return (v1[0] - v2[0]) * (v2[1] - v3[1]) + (v1[1] - v2[1]) * (v3[0] - v2[0]);
45 }
46
47 MINLINE float area_tri_signed_v2(const float v1[2], const float v2[2], const float v3[2])
48 {
49         return 0.5f * ((v1[0] - v2[0]) * (v2[1] - v3[1]) + (v1[1] - v2[1]) * (v3[0] - v2[0]));
50 }
51
52 MINLINE float area_tri_v2(const float v1[2], const float v2[2], const float v3[2])
53 {
54         return fabsf(area_tri_signed_v2(v1, v2, v3));
55 }
56
57 MINLINE float area_squared_tri_v2(const float v1[2], const float v2[2], const float v3[2])
58 {
59         float area = area_tri_signed_v2(v1, v2, v3);
60         return area * area;
61 }
62
63 /****************************** Spherical Harmonics **************************/
64
65 MINLINE void zero_sh(float r[9])
66 {
67         memset(r, 0, sizeof(float) * 9);
68 }
69
70 MINLINE void copy_sh_sh(float r[9], const float a[9])
71 {
72         memcpy(r, a, sizeof(float) * 9);
73 }
74
75 MINLINE void mul_sh_fl(float r[9], const float f)
76 {
77         int i;
78
79         for (i = 0; i < 9; i++)
80                 r[i] *= f;
81 }
82
83 MINLINE void add_sh_shsh(float r[9], const float a[9], const float b[9])
84 {
85         int i;
86
87         for (i = 0; i < 9; i++)
88                 r[i] = a[i] + b[i];
89 }
90
91 MINLINE float dot_shsh(const float a[9], const float b[9])
92 {
93         float r = 0.0f;
94         int i;
95
96         for (i = 0; i < 9; i++)
97                 r += a[i] * b[i];
98
99         return r;
100 }
101
102 MINLINE float diffuse_shv3(float sh[9], const float v[3])
103 {
104         /* See formula (13) in:
105          * "An Efficient Representation for Irradiance Environment Maps" */
106         static const float c1 = 0.429043f, c2 = 0.511664f, c3 = 0.743125f;
107         static const float c4 = 0.886227f, c5 = 0.247708f;
108         float x, y, z, sum;
109
110         x = v[0];
111         y = v[1];
112         z = v[2];
113
114         sum = c1 * sh[8] * (x * x - y * y);
115         sum += c3 * sh[6] * z * z;
116         sum += c4 * sh[0];
117         sum += -c5 * sh[6];
118         sum += 2.0f * c1 * (sh[4] * x * y + sh[7] * x * z + sh[5] * y * z);
119         sum += 2.0f * c2 * (sh[3] * x + sh[1] * y + sh[2] * z);
120
121         return sum;
122 }
123
124 MINLINE void vec_fac_to_sh(float r[9], const float v[3], const float f)
125 {
126         /* See formula (3) in:
127          * "An Efficient Representation for Irradiance Environment Maps" */
128         float sh[9], x, y, z;
129
130         x = v[0];
131         y = v[1];
132         z = v[2];
133
134         sh[0] = 0.282095f;
135
136         sh[1] = 0.488603f * y;
137         sh[2] = 0.488603f * z;
138         sh[3] = 0.488603f * x;
139
140         sh[4] = 1.092548f * x * y;
141         sh[5] = 1.092548f * y * z;
142         sh[6] = 0.315392f * (3.0f * z * z - 1.0f);
143         sh[7] = 1.092548f * x * z;
144         sh[8] = 0.546274f * (x * x - y * y);
145
146         mul_sh_fl(sh, f);
147         copy_sh_sh(r, sh);
148 }
149
150 MINLINE float eval_shv3(float sh[9], const float v[3])
151 {
152         float tmp[9];
153
154         vec_fac_to_sh(tmp, v, 1.0f);
155         return dot_shsh(tmp, sh);
156 }
157
158 MINLINE void madd_sh_shfl(float r[9], const float sh[9], const float f)
159 {
160         float tmp[9];
161
162         copy_sh_sh(tmp, sh);
163         mul_sh_fl(tmp, f);
164         add_sh_shsh(r, r, tmp);
165 }
166
167 /* get the 2 dominant axis values, 0==X, 1==Y, 2==Z */
168 MINLINE void axis_dominant_v3(int *r_axis_a, int *r_axis_b, const float axis[3])
169 {
170         const float xn = fabsf(axis[0]);
171         const float yn = fabsf(axis[1]);
172         const float zn = fabsf(axis[2]);
173
174         if      (zn >= xn && zn >= yn) { *r_axis_a = 0; *r_axis_b = 1; }
175         else if (yn >= xn && yn >= zn) { *r_axis_a = 0; *r_axis_b = 2; }
176         else                           { *r_axis_a = 1; *r_axis_b = 2; }
177 }
178
179 /* same as axis_dominant_v3 but return the max value */
180 MINLINE float axis_dominant_v3_max(int *r_axis_a, int *r_axis_b, const float axis[3])
181 {
182         const float xn = fabsf(axis[0]);
183         const float yn = fabsf(axis[1]);
184         const float zn = fabsf(axis[2]);
185
186         if      (zn >= xn && zn >= yn) { *r_axis_a = 0; *r_axis_b = 1; return zn; }
187         else if (yn >= xn && yn >= zn) { *r_axis_a = 0; *r_axis_b = 2; return yn; }
188         else                           { *r_axis_a = 1; *r_axis_b = 2; return xn; }
189 }
190
191 /* get the single dominant axis value, 0==X, 1==Y, 2==Z */
192 MINLINE int axis_dominant_v3_single(const float vec[3])
193 {
194         const float x = fabsf(vec[0]);
195         const float y = fabsf(vec[1]);
196         const float z = fabsf(vec[2]);
197         return ((x > y) ?
198                ((x > z) ? 0 : 2) :
199                ((y > z) ? 1 : 2));
200 }
201
202 /* the dominant axis of an orthogonal vector */
203 MINLINE int axis_dominant_v3_ortho_single(const float vec[3])
204 {
205         const float x = fabsf(vec[0]);
206         const float y = fabsf(vec[1]);
207         const float z = fabsf(vec[2]);
208         return ((x < y) ?
209                ((x < z) ? 0 : 2) :
210                ((y < z) ? 1 : 2));
211 }
212
213 MINLINE int max_axis_v3(const float vec[3])
214 {
215         const float x = vec[0];
216         const float y = vec[1];
217         const float z = vec[2];
218         return ((x > y) ?
219                ((x > z) ? 0 : 2) :
220                ((y > z) ? 1 : 2));
221 }
222
223 MINLINE int min_axis_v3(const float vec[3])
224 {
225         const float x = vec[0];
226         const float y = vec[1];
227         const float z = vec[2];
228         return ((x < y) ?
229                ((x < z) ? 0 : 2) :
230                ((y < z) ? 1 : 2));
231 }
232
233 /**
234  * Simple method to find how many tri's we need when we already know the corner+poly count.
235  *
236  * \param poly_count: The number of ngon's/tris (1-2 sided faces will give incorrect results)
237  * \param corner_count: also known as loops in BMesh/DNA
238  */
239 MINLINE int poly_to_tri_count(const int poly_count, const int corner_count)
240 {
241         BLI_assert(!poly_count || corner_count > poly_count * 2);
242         return corner_count - (poly_count * 2);
243 }
244
245 MINLINE float plane_point_side_v3(const float plane[4], const float co[3])
246 {
247         return dot_v3v3(co, plane) + plane[3];
248 }
249
250 /* useful to calculate an even width shell, by taking the angle between 2 planes.
251  * The return value is a scale on the offset.
252  * no angle between planes is 1.0, as the angle between the 2 planes approaches 180d
253  * the distance gets very high, 180d would be inf, but this case isn't valid */
254 MINLINE float shell_angle_to_dist(const float angle)
255 {
256         return (UNLIKELY(angle < SMALL_NUMBER)) ? 1.0f : fabsf(1.0f / cosf(angle));
257 }
258 /**
259  * equivalent to ``shell_angle_to_dist(angle_normalized_v3v3(a, b))``
260  */
261 MINLINE float shell_v3v3_normalized_to_dist(const float a[3], const float b[3])
262 {
263         const float angle_cos = fabsf(dot_v3v3(a, b));
264         BLI_ASSERT_UNIT_V3(a);
265         BLI_ASSERT_UNIT_V3(b);
266         return (UNLIKELY(angle_cos < SMALL_NUMBER)) ? 1.0f : (1.0f / angle_cos);
267 }
268 /**
269  * equivalent to ``shell_angle_to_dist(angle_normalized_v2v2(a, b))``
270  */
271 MINLINE float shell_v2v2_normalized_to_dist(const float a[2], const float b[2])
272 {
273         const float angle_cos = fabsf(dot_v2v2(a, b));
274         BLI_ASSERT_UNIT_V2(a);
275         BLI_ASSERT_UNIT_V2(b);
276         return (UNLIKELY(angle_cos < SMALL_NUMBER)) ? 1.0f : (1.0f / angle_cos);
277 }
278
279 /**
280  * equivalent to ``shell_angle_to_dist(angle_normalized_v3v3(a, b) / 2)``
281  */
282 MINLINE float shell_v3v3_mid_normalized_to_dist(const float a[3], const float b[3])
283 {
284         float angle_cos;
285         float ab[3];
286         BLI_ASSERT_UNIT_V3(a);
287         BLI_ASSERT_UNIT_V3(b);
288         add_v3_v3v3(ab, a, b);
289         angle_cos = (normalize_v3(ab) != 0.0f) ? fabsf(dot_v3v3(a, ab)) : 0.0f;
290         return (UNLIKELY(angle_cos < SMALL_NUMBER)) ? 1.0f : (1.0f / angle_cos);
291 }
292
293 /**
294  * equivalent to ``shell_angle_to_dist(angle_normalized_v2v2(a, b) / 2)``
295  */
296 MINLINE float shell_v2v2_mid_normalized_to_dist(const float a[2], const float b[2])
297 {
298         float angle_cos;
299         float ab[2];
300         BLI_ASSERT_UNIT_V2(a);
301         BLI_ASSERT_UNIT_V2(b);
302         add_v2_v2v2(ab, a, b);
303         angle_cos = (normalize_v2(ab) != 0.0f) ? fabsf(dot_v2v2(a, ab)) : 0.0f;
304         return (UNLIKELY(angle_cos < SMALL_NUMBER)) ? 1.0f : (1.0f / angle_cos);
305 }
306
307 #undef SMALL_NUMBER
308
309 #endif /* __MATH_GEOM_INLINE_C__ */