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