ClangFormat: apply to source, most of intern
[blender.git] / intern / cycles / util / util_transform.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_TRANSFORM_H__
18 #define __UTIL_TRANSFORM_H__
19
20 #ifndef __KERNEL_GPU__
21 #  include <string.h>
22 #endif
23
24 #include "util/util_math.h"
25 #include "util/util_types.h"
26
27 CCL_NAMESPACE_BEGIN
28
29 /* Affine transformation, stored as 4x3 matrix. */
30
31 typedef struct Transform {
32   float4 x, y, z;
33
34 #ifndef __KERNEL_GPU__
35   float4 operator[](int i) const
36   {
37     return *(&x + i);
38   }
39   float4 &operator[](int i)
40   {
41     return *(&x + i);
42   }
43 #endif
44 } Transform;
45
46 /* Transform decomposed in rotation/translation/scale. we use the same data
47  * structure as Transform, and tightly pack decomposition into it. first the
48  * rotation (4), then translation (3), then 3x3 scale matrix (9). */
49
50 typedef struct DecomposedTransform {
51   float4 x, y, z, w;
52 } DecomposedTransform;
53
54 /* Functions */
55
56 ccl_device_inline float3 transform_point(const Transform *t, const float3 a)
57 {
58   /* TODO(sergey): Disabled for now, causes crashes in certain cases. */
59 #if defined(__KERNEL_SSE__) && defined(__KERNEL_SSE2__)
60   ssef x, y, z, w, aa;
61   aa = a.m128;
62
63   x = _mm_loadu_ps(&t->x.x);
64   y = _mm_loadu_ps(&t->y.x);
65   z = _mm_loadu_ps(&t->z.x);
66   w = _mm_set_ps(1.0f, 0.0f, 0.0f, 0.0f);
67
68   _MM_TRANSPOSE4_PS(x, y, z, w);
69
70   ssef tmp = shuffle<0>(aa) * x;
71   tmp = madd(shuffle<1>(aa), y, tmp);
72   tmp = madd(shuffle<2>(aa), z, tmp);
73   tmp += w;
74
75   return float3(tmp.m128);
76 #else
77   float3 c = make_float3(a.x * t->x.x + a.y * t->x.y + a.z * t->x.z + t->x.w,
78                          a.x * t->y.x + a.y * t->y.y + a.z * t->y.z + t->y.w,
79                          a.x * t->z.x + a.y * t->z.y + a.z * t->z.z + t->z.w);
80
81   return c;
82 #endif
83 }
84
85 ccl_device_inline float3 transform_direction(const Transform *t, const float3 a)
86 {
87 #if defined(__KERNEL_SSE__) && defined(__KERNEL_SSE2__)
88   ssef x, y, z, w, aa;
89   aa = a.m128;
90   x = _mm_loadu_ps(&t->x.x);
91   y = _mm_loadu_ps(&t->y.x);
92   z = _mm_loadu_ps(&t->z.x);
93   w = _mm_setzero_ps();
94
95   _MM_TRANSPOSE4_PS(x, y, z, w);
96
97   ssef tmp = shuffle<0>(aa) * x;
98   tmp = madd(shuffle<1>(aa), y, tmp);
99   tmp = madd(shuffle<2>(aa), z, tmp);
100
101   return float3(tmp.m128);
102 #else
103   float3 c = make_float3(a.x * t->x.x + a.y * t->x.y + a.z * t->x.z,
104                          a.x * t->y.x + a.y * t->y.y + a.z * t->y.z,
105                          a.x * t->z.x + a.y * t->z.y + a.z * t->z.z);
106
107   return c;
108 #endif
109 }
110
111 ccl_device_inline float3 transform_direction_transposed(const Transform *t, const float3 a)
112 {
113   float3 x = make_float3(t->x.x, t->y.x, t->z.x);
114   float3 y = make_float3(t->x.y, t->y.y, t->z.y);
115   float3 z = make_float3(t->x.z, t->y.z, t->z.z);
116
117   return make_float3(dot(x, a), dot(y, a), dot(z, a));
118 }
119
120 ccl_device_inline Transform make_transform(float a,
121                                            float b,
122                                            float c,
123                                            float d,
124                                            float e,
125                                            float f,
126                                            float g,
127                                            float h,
128                                            float i,
129                                            float j,
130                                            float k,
131                                            float l)
132 {
133   Transform t;
134
135   t.x.x = a;
136   t.x.y = b;
137   t.x.z = c;
138   t.x.w = d;
139   t.y.x = e;
140   t.y.y = f;
141   t.y.z = g;
142   t.y.w = h;
143   t.z.x = i;
144   t.z.y = j;
145   t.z.z = k;
146   t.z.w = l;
147
148   return t;
149 }
150
151 /* Constructs a coordinate frame from a normalized normal. */
152 ccl_device_inline Transform make_transform_frame(float3 N)
153 {
154   const float3 dx0 = cross(make_float3(1.0f, 0.0f, 0.0f), N);
155   const float3 dx1 = cross(make_float3(0.0f, 1.0f, 0.0f), N);
156   const float3 dx = normalize((dot(dx0, dx0) > dot(dx1, dx1)) ? dx0 : dx1);
157   const float3 dy = normalize(cross(N, dx));
158   return make_transform(dx.x, dx.y, dx.z, 0.0f, dy.x, dy.y, dy.z, 0.0f, N.x, N.y, N.z, 0.0f);
159 }
160
161 #ifndef __KERNEL_GPU__
162
163 ccl_device_inline Transform operator*(const Transform a, const Transform b)
164 {
165   float4 c_x = make_float4(b.x.x, b.y.x, b.z.x, 0.0f);
166   float4 c_y = make_float4(b.x.y, b.y.y, b.z.y, 0.0f);
167   float4 c_z = make_float4(b.x.z, b.y.z, b.z.z, 0.0f);
168   float4 c_w = make_float4(b.x.w, b.y.w, b.z.w, 1.0f);
169
170   Transform t;
171   t.x = make_float4(dot(a.x, c_x), dot(a.x, c_y), dot(a.x, c_z), dot(a.x, c_w));
172   t.y = make_float4(dot(a.y, c_x), dot(a.y, c_y), dot(a.y, c_z), dot(a.y, c_w));
173   t.z = make_float4(dot(a.z, c_x), dot(a.z, c_y), dot(a.z, c_z), dot(a.z, c_w));
174
175   return t;
176 }
177
178 ccl_device_inline void print_transform(const char *label, const Transform &t)
179 {
180   print_float4(label, t.x);
181   print_float4(label, t.y);
182   print_float4(label, t.z);
183   printf("\n");
184 }
185
186 ccl_device_inline Transform transform_translate(float3 t)
187 {
188   return make_transform(1, 0, 0, t.x, 0, 1, 0, t.y, 0, 0, 1, t.z);
189 }
190
191 ccl_device_inline Transform transform_translate(float x, float y, float z)
192 {
193   return transform_translate(make_float3(x, y, z));
194 }
195
196 ccl_device_inline Transform transform_scale(float3 s)
197 {
198   return make_transform(s.x, 0, 0, 0, 0, s.y, 0, 0, 0, 0, s.z, 0);
199 }
200
201 ccl_device_inline Transform transform_scale(float x, float y, float z)
202 {
203   return transform_scale(make_float3(x, y, z));
204 }
205
206 ccl_device_inline Transform transform_rotate(float angle, float3 axis)
207 {
208   float s = sinf(angle);
209   float c = cosf(angle);
210   float t = 1.0f - c;
211
212   axis = normalize(axis);
213
214   return make_transform(axis.x * axis.x * t + c,
215                         axis.x * axis.y * t - s * axis.z,
216                         axis.x * axis.z * t + s * axis.y,
217                         0.0f,
218
219                         axis.y * axis.x * t + s * axis.z,
220                         axis.y * axis.y * t + c,
221                         axis.y * axis.z * t - s * axis.x,
222                         0.0f,
223
224                         axis.z * axis.x * t - s * axis.y,
225                         axis.z * axis.y * t + s * axis.x,
226                         axis.z * axis.z * t + c,
227                         0.0f);
228 }
229
230 /* Euler is assumed to be in XYZ order. */
231 ccl_device_inline Transform transform_euler(float3 euler)
232 {
233   return transform_rotate(euler.z, make_float3(0.0f, 0.0f, 1.0f)) *
234          transform_rotate(euler.y, make_float3(0.0f, 1.0f, 0.0f)) *
235          transform_rotate(euler.x, make_float3(1.0f, 0.0f, 0.0f));
236 }
237
238 ccl_device_inline Transform transform_identity()
239 {
240   return transform_scale(1.0f, 1.0f, 1.0f);
241 }
242
243 ccl_device_inline bool operator==(const Transform &A, const Transform &B)
244 {
245   return memcmp(&A, &B, sizeof(Transform)) == 0;
246 }
247
248 ccl_device_inline bool operator!=(const Transform &A, const Transform &B)
249 {
250   return !(A == B);
251 }
252
253 ccl_device_inline float3 transform_get_column(const Transform *t, int column)
254 {
255   return make_float3(t->x[column], t->y[column], t->z[column]);
256 }
257
258 ccl_device_inline void transform_set_column(Transform *t, int column, float3 value)
259 {
260   t->x[column] = value.x;
261   t->y[column] = value.y;
262   t->z[column] = value.z;
263 }
264
265 Transform transform_inverse(const Transform &a);
266 Transform transform_transposed_inverse(const Transform &a);
267
268 ccl_device_inline bool transform_uniform_scale(const Transform &tfm, float &scale)
269 {
270   /* the epsilon here is quite arbitrary, but this function is only used for
271    * surface area and bump, where we expect it to not be so sensitive */
272   float eps = 1e-6f;
273
274   float sx = len_squared(float4_to_float3(tfm.x));
275   float sy = len_squared(float4_to_float3(tfm.y));
276   float sz = len_squared(float4_to_float3(tfm.z));
277   float stx = len_squared(transform_get_column(&tfm, 0));
278   float sty = len_squared(transform_get_column(&tfm, 1));
279   float stz = len_squared(transform_get_column(&tfm, 2));
280
281   if (fabsf(sx - sy) < eps && fabsf(sx - sz) < eps && fabsf(sx - stx) < eps &&
282       fabsf(sx - sty) < eps && fabsf(sx - stz) < eps) {
283     scale = sx;
284     return true;
285   }
286
287   return false;
288 }
289
290 ccl_device_inline bool transform_negative_scale(const Transform &tfm)
291 {
292   float3 c0 = transform_get_column(&tfm, 0);
293   float3 c1 = transform_get_column(&tfm, 1);
294   float3 c2 = transform_get_column(&tfm, 2);
295
296   return (dot(cross(c0, c1), c2) < 0.0f);
297 }
298
299 ccl_device_inline Transform transform_clear_scale(const Transform &tfm)
300 {
301   Transform ntfm = tfm;
302
303   transform_set_column(&ntfm, 0, normalize(transform_get_column(&ntfm, 0)));
304   transform_set_column(&ntfm, 1, normalize(transform_get_column(&ntfm, 1)));
305   transform_set_column(&ntfm, 2, normalize(transform_get_column(&ntfm, 2)));
306
307   return ntfm;
308 }
309
310 ccl_device_inline Transform transform_empty()
311 {
312   return make_transform(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
313 }
314
315 #endif
316
317 /* Motion Transform */
318
319 ccl_device_inline float4 quat_interpolate(float4 q1, float4 q2, float t)
320 {
321   /* use simpe nlerp instead of slerp. it's faster and almost the same */
322   return normalize((1.0f - t) * q1 + t * q2);
323
324 #if 0
325   /* note: this does not ensure rotation around shortest angle, q1 and q2
326    * are assumed to be matched already in transform_motion_decompose */
327   float costheta = dot(q1, q2);
328
329   /* possible optimization: it might be possible to precompute theta/qperp */
330
331   if(costheta > 0.9995f) {
332     /* linear interpolation in degenerate case */
333     return normalize((1.0f - t)*q1 + t*q2);
334   }
335   else  {
336     /* slerp */
337     float theta = acosf(clamp(costheta, -1.0f, 1.0f));
338     float4 qperp = normalize(q2 - q1 * costheta);
339     float thetap = theta * t;
340     return q1 * cosf(thetap) + qperp * sinf(thetap);
341   }
342 #endif
343 }
344
345 ccl_device_inline Transform transform_quick_inverse(Transform M)
346 {
347   /* possible optimization: can we avoid doing this altogether and construct
348    * the inverse matrix directly from negated translation, transposed rotation,
349    * scale can be inverted but what about shearing? */
350   Transform R;
351   float det = M.x.x * (M.z.z * M.y.y - M.z.y * M.y.z) - M.y.x * (M.z.z * M.x.y - M.z.y * M.x.z) +
352               M.z.x * (M.y.z * M.x.y - M.y.y * M.x.z);
353   if (det == 0.0f) {
354     M.x.x += 1e-8f;
355     M.y.y += 1e-8f;
356     M.z.z += 1e-8f;
357     det = M.x.x * (M.z.z * M.y.y - M.z.y * M.y.z) - M.y.x * (M.z.z * M.x.y - M.z.y * M.x.z) +
358           M.z.x * (M.y.z * M.x.y - M.y.y * M.x.z);
359   }
360   det = (det != 0.0f) ? 1.0f / det : 0.0f;
361
362   float3 Rx = det * make_float3(M.z.z * M.y.y - M.z.y * M.y.z,
363                                 M.z.y * M.x.z - M.z.z * M.x.y,
364                                 M.y.z * M.x.y - M.y.y * M.x.z);
365   float3 Ry = det * make_float3(M.z.x * M.y.z - M.z.z * M.y.x,
366                                 M.z.z * M.x.x - M.z.x * M.x.z,
367                                 M.y.x * M.x.z - M.y.z * M.x.x);
368   float3 Rz = det * make_float3(M.z.y * M.y.x - M.z.x * M.y.y,
369                                 M.z.x * M.x.y - M.z.y * M.x.x,
370                                 M.y.y * M.x.x - M.y.x * M.x.y);
371   float3 T = -make_float3(M.x.w, M.y.w, M.z.w);
372
373   R.x = make_float4(Rx.x, Rx.y, Rx.z, dot(Rx, T));
374   R.y = make_float4(Ry.x, Ry.y, Ry.z, dot(Ry, T));
375   R.z = make_float4(Rz.x, Rz.y, Rz.z, dot(Rz, T));
376
377   return R;
378 }
379
380 ccl_device_inline void transform_compose(Transform *tfm, const DecomposedTransform *decomp)
381 {
382   /* rotation */
383   float q0, q1, q2, q3, qda, qdb, qdc, qaa, qab, qac, qbb, qbc, qcc;
384
385   q0 = M_SQRT2_F * decomp->x.w;
386   q1 = M_SQRT2_F * decomp->x.x;
387   q2 = M_SQRT2_F * decomp->x.y;
388   q3 = M_SQRT2_F * decomp->x.z;
389
390   qda = q0 * q1;
391   qdb = q0 * q2;
392   qdc = q0 * q3;
393   qaa = q1 * q1;
394   qab = q1 * q2;
395   qac = q1 * q3;
396   qbb = q2 * q2;
397   qbc = q2 * q3;
398   qcc = q3 * q3;
399
400   float3 rotation_x = make_float3(1.0f - qbb - qcc, -qdc + qab, qdb + qac);
401   float3 rotation_y = make_float3(qdc + qab, 1.0f - qaa - qcc, -qda + qbc);
402   float3 rotation_z = make_float3(-qdb + qac, qda + qbc, 1.0f - qaa - qbb);
403
404   /* scale */
405   float3 scale_x = make_float3(decomp->y.w, decomp->z.z, decomp->w.y);
406   float3 scale_y = make_float3(decomp->z.x, decomp->z.w, decomp->w.z);
407   float3 scale_z = make_float3(decomp->z.y, decomp->w.x, decomp->w.w);
408
409   /* compose with translation */
410   tfm->x = make_float4(
411       dot(rotation_x, scale_x), dot(rotation_x, scale_y), dot(rotation_x, scale_z), decomp->y.x);
412   tfm->y = make_float4(
413       dot(rotation_y, scale_x), dot(rotation_y, scale_y), dot(rotation_y, scale_z), decomp->y.y);
414   tfm->z = make_float4(
415       dot(rotation_z, scale_x), dot(rotation_z, scale_y), dot(rotation_z, scale_z), decomp->y.z);
416 }
417
418 /* Interpolate from array of decomposed transforms. */
419 ccl_device void transform_motion_array_interpolate(Transform *tfm,
420                                                    const ccl_global DecomposedTransform *motion,
421                                                    uint numsteps,
422                                                    float time)
423 {
424   /* Figure out which steps we need to interpolate. */
425   int maxstep = numsteps - 1;
426   int step = min((int)(time * maxstep), maxstep - 1);
427   float t = time * maxstep - step;
428
429   const ccl_global DecomposedTransform *a = motion + step;
430   const ccl_global DecomposedTransform *b = motion + step + 1;
431
432   /* Interpolate rotation, translation and scale. */
433   DecomposedTransform decomp;
434   decomp.x = quat_interpolate(a->x, b->x, t);
435   decomp.y = (1.0f - t) * a->y + t * b->y;
436   decomp.z = (1.0f - t) * a->z + t * b->z;
437   decomp.w = (1.0f - t) * a->w + t * b->w;
438
439   /* Compose rotation, translation, scale into matrix. */
440   transform_compose(tfm, &decomp);
441 }
442
443 #ifndef __KERNEL_GPU__
444
445 #  ifdef WITH_EMBREE
446 ccl_device void transform_motion_array_interpolate_straight(
447     Transform *tfm, const ccl_global DecomposedTransform *motion, uint numsteps, float time)
448 {
449   /* Figure out which steps we need to interpolate. */
450   int maxstep = numsteps - 1;
451   int step = min((int)(time * maxstep), maxstep - 1);
452   float t = time * maxstep - step;
453
454   const ccl_global DecomposedTransform *a = motion + step;
455   const ccl_global DecomposedTransform *b = motion + step + 1;
456   Transform step1, step2;
457
458   transform_compose(&step1, a);
459   transform_compose(&step2, b);
460
461   /* matrix lerp */
462   tfm->x = (1.0f - t) * step1.x + t * step2.x;
463   tfm->y = (1.0f - t) * step1.y + t * step2.y;
464   tfm->z = (1.0f - t) * step1.z + t * step2.z;
465 }
466 #  endif
467
468 class BoundBox2D;
469
470 ccl_device_inline bool operator==(const DecomposedTransform &A, const DecomposedTransform &B)
471 {
472   return memcmp(&A, &B, sizeof(DecomposedTransform)) == 0;
473 }
474
475 float4 transform_to_quat(const Transform &tfm);
476 void transform_motion_decompose(DecomposedTransform *decomp, const Transform *motion, size_t size);
477 Transform transform_from_viewplane(BoundBox2D &viewplane);
478
479 #endif
480
481 /* TODO(sergey): This is only for until we've got OpenCL 2.0
482  * on all devices we consider supported. It'll be replaced with
483  * generic address space.
484  */
485
486 #ifdef __KERNEL_OPENCL__
487
488 #  define OPENCL_TRANSFORM_ADDRSPACE_GLUE(a, b) a##b
489 #  define OPENCL_TRANSFORM_ADDRSPACE_DECLARE(function) \
490     ccl_device_inline float3 OPENCL_TRANSFORM_ADDRSPACE_GLUE(function, _addrspace)( \
491         ccl_addr_space const Transform *t, const float3 a) \
492     { \
493       Transform private_tfm = *t; \
494       return function(&private_tfm, a); \
495     }
496
497 OPENCL_TRANSFORM_ADDRSPACE_DECLARE(transform_point)
498 OPENCL_TRANSFORM_ADDRSPACE_DECLARE(transform_direction)
499 OPENCL_TRANSFORM_ADDRSPACE_DECLARE(transform_direction_transposed)
500
501 #  undef OPENCL_TRANSFORM_ADDRSPACE_DECLARE
502 #  undef OPENCL_TRANSFORM_ADDRSPACE_GLUE
503 #  define transform_point_auto transform_point_addrspace
504 #  define transform_direction_auto transform_direction_addrspace
505 #  define transform_direction_transposed_auto transform_direction_transposed_addrspace
506 #else
507 #  define transform_point_auto transform_point
508 #  define transform_direction_auto transform_direction
509 #  define transform_direction_transposed_auto transform_direction_transposed
510 #endif
511
512 CCL_NAMESPACE_END
513
514 #endif /* __UTIL_TRANSFORM_H__ */