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