Cycles: change __device and similar qualifiers to ccl_device in kernel code.
[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_math.h"
25 #include "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  * For the DecompMotionTransform we drop scale from pre/post. */
45
46 typedef struct ccl_may_alias MotionTransform {
47         Transform pre;
48         Transform mid;
49         Transform post;
50 } MotionTransform;
51
52 typedef struct DecompMotionTransform {
53         Transform mid;
54         float4 pre_x, pre_y;
55         float4 post_x, post_y;
56 } DecompMotionTransform;
57
58 /* Functions */
59
60 ccl_device_inline float3 transform_perspective(const Transform *t, const float3 a)
61 {
62         float4 b = make_float4(a.x, a.y, a.z, 1.0f);
63         float3 c = make_float3(dot(t->x, b), dot(t->y, b), dot(t->z, b));
64         float w = dot(t->w, b);
65
66         return (w != 0.0f)? c/w: make_float3(0.0f, 0.0f, 0.0f);
67 }
68
69 ccl_device_inline float3 transform_point(const Transform *t, const float3 a)
70 {
71         float3 c = make_float3(
72                 a.x*t->x.x + a.y*t->x.y + a.z*t->x.z + t->x.w,
73                 a.x*t->y.x + a.y*t->y.y + a.z*t->y.z + t->y.w,
74                 a.x*t->z.x + a.y*t->z.y + a.z*t->z.z + t->z.w);
75
76         return c;
77 }
78
79 ccl_device_inline float3 transform_direction(const Transform *t, const float3 a)
80 {
81         float3 c = make_float3(
82                 a.x*t->x.x + a.y*t->x.y + a.z*t->x.z,
83                 a.x*t->y.x + a.y*t->y.y + a.z*t->y.z,
84                 a.x*t->z.x + a.y*t->z.y + a.z*t->z.z);
85
86         return c;
87 }
88
89 ccl_device_inline float3 transform_direction_transposed(const Transform *t, const float3 a)
90 {
91         float3 x = make_float3(t->x.x, t->y.x, t->z.x);
92         float3 y = make_float3(t->x.y, t->y.y, t->z.y);
93         float3 z = make_float3(t->x.z, t->y.z, t->z.z);
94
95         return make_float3(dot(x, a), dot(y, a), dot(z, a));
96 }
97
98 ccl_device_inline Transform transform_transpose(const Transform a)
99 {
100         Transform t;
101
102         t.x.x = a.x.x; t.x.y = a.y.x; t.x.z = a.z.x; t.x.w = a.w.x;
103         t.y.x = a.x.y; t.y.y = a.y.y; t.y.z = a.z.y; t.y.w = a.w.y;
104         t.z.x = a.x.z; t.z.y = a.y.z; t.z.z = a.z.z; t.z.w = a.w.z;
105         t.w.x = a.x.w; t.w.y = a.y.w; t.w.z = a.z.w; t.w.w = a.w.w;
106
107         return t;
108 }
109
110 ccl_device_inline Transform make_transform(float a, float b, float c, float d,
111                                                                         float e, float f, float g, float h,
112                                                                         float i, float j, float k, float l,
113                                                                         float m, float n, float o, float p)
114 {
115         Transform t;
116
117         t.x.x = a; t.x.y = b; t.x.z = c; t.x.w = d;
118         t.y.x = e; t.y.y = f; t.y.z = g; t.y.w = h;
119         t.z.x = i; t.z.y = j; t.z.z = k; t.z.w = l;
120         t.w.x = m; t.w.y = n; t.w.z = o; t.w.w = p;
121
122         return t;
123 }
124
125 #ifndef __KERNEL_GPU__
126
127 ccl_device_inline Transform operator*(const Transform a, const Transform b)
128 {
129         Transform c = transform_transpose(b);
130         Transform t;
131
132         t.x = make_float4(dot(a.x, c.x), dot(a.x, c.y), dot(a.x, c.z), dot(a.x, c.w));
133         t.y = make_float4(dot(a.y, c.x), dot(a.y, c.y), dot(a.y, c.z), dot(a.y, c.w));
134         t.z = make_float4(dot(a.z, c.x), dot(a.z, c.y), dot(a.z, c.z), dot(a.z, c.w));
135         t.w = make_float4(dot(a.w, c.x), dot(a.w, c.y), dot(a.w, c.z), dot(a.w, c.w));
136
137         return t;
138 }
139
140 ccl_device_inline void print_transform(const char *label, const Transform& t)
141 {
142         print_float4(label, t.x);
143         print_float4(label, t.y);
144         print_float4(label, t.z);
145         print_float4(label, t.w);
146         printf("\n");
147 }
148
149 ccl_device_inline Transform transform_translate(float3 t)
150 {
151         return make_transform(
152                 1, 0, 0, t.x,
153                 0, 1, 0, t.y,
154                 0, 0, 1, t.z,
155                 0, 0, 0, 1);
156 }
157
158 ccl_device_inline Transform transform_translate(float x, float y, float z)
159 {
160         return transform_translate(make_float3(x, y, z));
161 }
162
163 ccl_device_inline Transform transform_scale(float3 s)
164 {
165         return make_transform(
166                 s.x, 0, 0, 0,
167                 0, s.y, 0, 0,
168                 0, 0, s.z, 0,
169                 0, 0, 0, 1);
170 }
171
172 ccl_device_inline Transform transform_scale(float x, float y, float z)
173 {
174         return transform_scale(make_float3(x, y, z));
175 }
176
177 ccl_device_inline Transform transform_perspective(float fov, float n, float f)
178 {
179         Transform persp = make_transform(
180                 1, 0, 0, 0,
181                 0, 1, 0, 0,
182                 0, 0, f / (f - n), -f*n / (f - n),
183                 0, 0, 1, 0);
184
185         float inv_angle = 1.0f/tanf(0.5f*fov);
186
187         Transform scale = transform_scale(inv_angle, inv_angle, 1);
188
189         return scale * persp;
190 }
191
192 ccl_device_inline Transform transform_rotate(float angle, float3 axis)
193 {
194         float s = sinf(angle);
195         float c = cosf(angle);
196         float t = 1.0f - c;
197
198         axis = normalize(axis);
199
200         return make_transform(
201                 axis.x*axis.x*t + c,
202                 axis.x*axis.y*t - s*axis.z,
203                 axis.x*axis.z*t + s*axis.y,
204                 0.0f,
205
206                 axis.y*axis.x*t + s*axis.z,
207                 axis.y*axis.y*t + c,
208                 axis.y*axis.z*t - s*axis.x,
209                 0.0f,
210
211                 axis.z*axis.x*t - s*axis.y,
212                 axis.z*axis.y*t + s*axis.x,
213                 axis.z*axis.z*t + c,
214                 0.0f,
215
216                 0.0f, 0.0f, 0.0f, 1.0f);
217 }
218
219 ccl_device_inline Transform transform_euler(float3 euler)
220 {
221         return
222                 transform_rotate(euler.x, make_float3(1.0f, 0.0f, 0.0f)) *
223                 transform_rotate(euler.y, make_float3(0.0f, 1.0f, 0.0f)) *
224                 transform_rotate(euler.z, make_float3(0.0f, 0.0f, 1.0f));
225 }
226
227 ccl_device_inline Transform transform_orthographic(float znear, float zfar)
228 {
229         return transform_scale(1.0f, 1.0f, 1.0f / (zfar-znear)) *
230                 transform_translate(0.0f, 0.0f, -znear);
231 }
232
233 ccl_device_inline Transform transform_identity()
234 {
235         return transform_scale(1.0f, 1.0f, 1.0f);
236 }
237
238 ccl_device_inline bool operator==(const Transform& A, const Transform& B)
239 {
240         return memcmp(&A, &B, sizeof(Transform)) == 0;
241 }
242
243 ccl_device_inline bool operator!=(const Transform& A, const Transform& B)
244 {
245         return !(A == B);
246 }
247
248 ccl_device_inline float3 transform_get_column(const Transform *t, int column)
249 {
250         return make_float3(t->x[column], t->y[column], t->z[column]);
251 }
252
253 ccl_device_inline void transform_set_column(Transform *t, int column, float3 value)
254 {
255         t->x[column] = value.x;
256         t->y[column] = value.y;
257         t->z[column] = value.z;
258 }
259
260 Transform transform_inverse(const Transform& a);
261
262 ccl_device_inline bool transform_uniform_scale(const Transform& tfm, float& scale)
263 {
264         /* the epsilon here is quite arbitrary, but this function is only used for
265          * surface area and bump, where we except it to not be so sensitive */
266         Transform ttfm = transform_transpose(tfm);
267         float eps = 1e-6f; 
268         
269         float sx = len_squared(float4_to_float3(tfm.x));
270         float sy = len_squared(float4_to_float3(tfm.y));
271         float sz = len_squared(float4_to_float3(tfm.z));
272         float stx = len_squared(float4_to_float3(ttfm.x));
273         float sty = len_squared(float4_to_float3(ttfm.y));
274         float stz = len_squared(float4_to_float3(ttfm.z));
275
276         if(fabsf(sx - sy) < eps && fabsf(sx - sz) < eps &&
277            fabsf(sx - stx) < eps && fabsf(sx - sty) < eps &&
278            fabsf(sx - stz) < eps) {
279                 scale = sx;
280                 return true;
281         }
282    
283    return false;
284 }
285
286 ccl_device_inline bool transform_negative_scale(const Transform& tfm)
287 {
288         float3 c0 = transform_get_column(&tfm, 0);
289         float3 c1 = transform_get_column(&tfm, 1);
290         float3 c2 = transform_get_column(&tfm, 2);
291
292         return (dot(cross(c0, c1), c2) < 0.0f);
293 }
294
295 ccl_device_inline Transform transform_clear_scale(const Transform& tfm)
296 {
297         Transform ntfm = tfm;
298
299         transform_set_column(&ntfm, 0, normalize(transform_get_column(&ntfm, 0)));
300         transform_set_column(&ntfm, 1, normalize(transform_get_column(&ntfm, 1)));
301         transform_set_column(&ntfm, 2, normalize(transform_get_column(&ntfm, 2)));
302
303         return ntfm;
304 }
305
306 #endif
307
308 /* Motion Transform */
309
310 ccl_device_inline float4 quat_interpolate(float4 q1, float4 q2, float t)
311 {
312         /* use simpe nlerp instead of slerp. it's faster and almost the same */
313         return normalize((1.0f - t)*q1 + t*q2);
314
315 #if 0
316         /* note: this does not ensure rotation around shortest angle, q1 and q2
317          * are assumed to be matched already in transform_motion_decompose */
318         float costheta = dot(q1, q2);
319
320         /* possible optimization: it might be possible to precompute theta/qperp */
321
322         if(costheta > 0.9995f) {
323                 /* linear interpolation in degenerate case */
324                 return normalize((1.0f - t)*q1 + t*q2);
325         }
326         else  {
327                 /* slerp */
328                 float theta = acosf(clamp(costheta, -1.0f, 1.0f));
329                 float4 qperp = normalize(q2 - q1 * costheta);
330                 float thetap = theta * t;
331                 return q1 * cosf(thetap) + qperp * sinf(thetap);
332         }
333 #endif
334 }
335
336 ccl_device_inline Transform transform_quick_inverse(Transform M)
337 {
338         /* possible optimization: can we avoid doing this altogether and construct
339          * the inverse matrix directly from negated translation, transposed rotation,
340          * scale can be inverted but what about shearing? */
341         Transform R;
342         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);
343
344         det = (det != 0.0f)? 1.0f/det: 0.0f;
345
346         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);
347         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);
348         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);
349         float3 T = -make_float3(M.x.w, M.y.w, M.z.w);
350
351         R.x = make_float4(Rx.x, Rx.y, Rx.z, dot(Rx, T));
352         R.y = make_float4(Ry.x, Ry.y, Ry.z, dot(Ry, T));
353         R.z = make_float4(Rz.x, Rz.y, Rz.z, dot(Rz, T));
354         R.w = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
355
356         return R;
357 }
358
359 ccl_device_inline void transform_compose(Transform *tfm, const Transform *decomp)
360 {
361         /* rotation */
362         float q0, q1, q2, q3, qda, qdb, qdc, qaa, qab, qac, qbb, qbc, qcc;
363
364         q0 = M_SQRT2_F * decomp->x.w;
365         q1 = M_SQRT2_F * decomp->x.x;
366         q2 = M_SQRT2_F * decomp->x.y;
367         q3 = M_SQRT2_F * decomp->x.z;
368
369         qda = q0*q1;
370         qdb = q0*q2;
371         qdc = q0*q3;
372         qaa = q1*q1;
373         qab = q1*q2;
374         qac = q1*q3;
375         qbb = q2*q2;
376         qbc = q2*q3;
377         qcc = q3*q3;
378
379         float3 rotation_x = make_float3(1.0f-qbb-qcc, -qdc+qab, qdb+qac);
380         float3 rotation_y = make_float3(qdc+qab, 1.0f-qaa-qcc, -qda+qbc);
381         float3 rotation_z = make_float3(-qdb+qac, qda+qbc, 1.0f-qaa-qbb);
382
383         /* scale */
384         float3 scale_x = make_float3(decomp->y.w, decomp->z.z, decomp->w.y);
385         float3 scale_y = make_float3(decomp->z.x, decomp->z.w, decomp->w.z);
386         float3 scale_z = make_float3(decomp->z.y, decomp->w.x, decomp->w.w);
387
388         /* compose with translation */
389         tfm->x = make_float4(dot(rotation_x, scale_x), dot(rotation_x, scale_y), dot(rotation_x, scale_z), decomp->y.x);
390         tfm->y = make_float4(dot(rotation_y, scale_x), dot(rotation_y, scale_y), dot(rotation_y, scale_z), decomp->y.y);
391         tfm->z = make_float4(dot(rotation_z, scale_x), dot(rotation_z, scale_y), dot(rotation_z, scale_z), decomp->y.z);
392         tfm->w = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
393 }
394
395 /* Disabled for now, need arc-length parametrization for constant speed motion.
396  * #define CURVED_MOTION_INTERPOLATE */
397
398 ccl_device void transform_motion_interpolate(Transform *tfm, const DecompMotionTransform *motion, float t)
399 {
400         /* possible optimization: is it worth it adding a check to skip scaling?
401          * it's probably quite uncommon to have scaling objects. or can we skip
402          * just shearing perhaps? */
403         Transform decomp;
404
405 #ifdef CURVED_MOTION_INTERPOLATE
406         /* 3 point bezier curve interpolation for position */
407         float3 Ppre = float4_to_float3(motion->pre_y);
408         float3 Pmid = float4_to_float3(motion->mid.y);
409         float3 Ppost = float4_to_float3(motion->post_y);
410
411         float3 Pcontrol = 2.0f*Pmid - 0.5f*(Ppre + Ppost);
412         float3 P = Ppre*t*t + Pcontrol*2.0f*t*(1.0f - t) + Ppost*(1.0f - t)*(1.0f - t);
413
414         decomp.y.x = P.x;
415         decomp.y.y = P.y;
416         decomp.y.z = P.z;
417 #endif
418
419         /* linear interpolation for rotation and scale */
420         if(t < 0.5f) {
421                 t *= 2.0f;
422
423                 decomp.x = quat_interpolate(motion->pre_x, motion->mid.x, t);
424 #ifdef CURVED_MOTION_INTERPOLATE
425                 decomp.y.w = (1.0f - t)*motion->pre_y.w + t*motion->mid.y.w;
426 #else
427                 decomp.y = (1.0f - t)*motion->pre_y + t*motion->mid.y;
428 #endif
429         }
430         else {
431                 t = (t - 0.5f)*2.0f;
432
433                 decomp.x = quat_interpolate(motion->mid.x, motion->post_x, t);
434 #ifdef CURVED_MOTION_INTERPOLATE
435                 decomp.y.w = (1.0f - t)*motion->mid.y.w + t*motion->post_y.w;
436 #else
437                 decomp.y = (1.0f - t)*motion->mid.y + t*motion->post_y;
438 #endif
439         }
440
441         decomp.z = motion->mid.z;
442         decomp.w = motion->mid.w;
443
444         /* compose rotation, translation, scale into matrix */
445         transform_compose(tfm, &decomp);
446 }
447
448 #ifndef __KERNEL_GPU__
449
450 ccl_device_inline bool operator==(const MotionTransform& A, const MotionTransform& B)
451 {
452         return (A.pre == B.pre && A.post == B.post);
453 }
454
455 float4 transform_to_quat(const Transform& tfm);
456 void transform_motion_decompose(DecompMotionTransform *decomp, const MotionTransform *motion, const Transform *mid);
457
458 #endif
459
460 CCL_NAMESPACE_END
461
462 #endif /* __UTIL_TRANSFORM_H__ */
463