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