Cycles: motion blur is now curved and passes exactly through the midpoint.
[blender-staging.git] / intern / cycles / util / util_transform.h
1 /*
2  * Copyright 2011, Blender Foundation.
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18
19 #ifndef __UTIL_TRANSFORM_H__
20 #define __UTIL_TRANSFORM_H__
21
22 #ifndef __KERNEL_GPU__
23 #include <string.h>
24 #endif
25
26 #include "util_math.h"
27 #include "util_types.h"
28
29 CCL_NAMESPACE_BEGIN
30
31 /* Data Types */
32
33 typedef struct Transform {
34         float4 x, y, z, w; /* rows */
35
36 #ifndef __KERNEL_GPU__
37         float4 operator[](int i) const { return *(&x + i); }
38         float4& operator[](int i) { return *(&x + i); }
39 #endif
40 } Transform;
41
42 typedef struct MotionTransform {
43         Transform pre;
44         Transform mid;
45         Transform post;
46 } MotionTransform;
47
48 /* transform decomposed in rotation/translation/scale. we use the same data
49  * structure as Transform, and tightly pack decomposition into it. first the
50  * rotation (4), then translation (3), then 3x3 scale matrix (9) */
51
52 /* Functions */
53
54 __device_inline float3 transform_perspective(const Transform *t, const float3 a)
55 {
56         float4 b = make_float4(a.x, a.y, a.z, 1.0f);
57         float3 c = make_float3(dot(t->x, b), dot(t->y, b), dot(t->z, b));
58         float w = dot(t->w, b);
59
60         return (w != 0.0f)? c/w: make_float3(0.0f, 0.0f, 0.0f);
61 }
62
63 __device_inline float3 transform_point(const Transform *t, const float3 a)
64 {
65         float3 c = make_float3(
66                 a.x*t->x.x + a.y*t->x.y + a.z*t->x.z + t->x.w,
67                 a.x*t->y.x + a.y*t->y.y + a.z*t->y.z + t->y.w,
68                 a.x*t->z.x + a.y*t->z.y + a.z*t->z.z + t->z.w);
69
70         return c;
71 }
72
73 __device_inline float3 transform_direction(const Transform *t, const float3 a)
74 {
75         float3 c = make_float3(
76                 a.x*t->x.x + a.y*t->x.y + a.z*t->x.z,
77                 a.x*t->y.x + a.y*t->y.y + a.z*t->y.z,
78                 a.x*t->z.x + a.y*t->z.y + a.z*t->z.z);
79
80         return c;
81 }
82
83 __device_inline float3 transform_direction_transposed(const Transform *t, const float3 a)
84 {
85         float3 x = make_float3(t->x.x, t->y.x, t->z.x);
86         float3 y = make_float3(t->x.y, t->y.y, t->z.y);
87         float3 z = make_float3(t->x.z, t->y.z, t->z.z);
88
89         return make_float3(dot(x, a), dot(y, a), dot(z, a));
90 }
91
92 #ifndef __KERNEL_GPU__
93
94 __device_inline void print_transform(const char *label, const Transform& t)
95 {
96         print_float4(label, t.x);
97         print_float4(label, t.y);
98         print_float4(label, t.z);
99         print_float4(label, t.w);
100         printf("\n");
101 }
102
103 __device_inline Transform transform_transpose(const Transform a)
104 {
105         Transform t;
106
107         t.x.x = a.x.x; t.x.y = a.y.x; t.x.z = a.z.x; t.x.w = a.w.x;
108         t.y.x = a.x.y; t.y.y = a.y.y; t.y.z = a.z.y; t.y.w = a.w.y;
109         t.z.x = a.x.z; t.z.y = a.y.z; t.z.z = a.z.z; t.z.w = a.w.z;
110         t.w.x = a.x.w; t.w.y = a.y.w; t.w.z = a.z.w; t.w.w = a.w.w;
111
112         return t;
113 }
114
115 __device_inline Transform operator*(const Transform a, const Transform b)
116 {
117         Transform c = transform_transpose(b);
118         Transform t;
119
120         t.x = make_float4(dot(a.x, c.x), dot(a.x, c.y), dot(a.x, c.z), dot(a.x, c.w));
121         t.y = make_float4(dot(a.y, c.x), dot(a.y, c.y), dot(a.y, c.z), dot(a.y, c.w));
122         t.z = make_float4(dot(a.z, c.x), dot(a.z, c.y), dot(a.z, c.z), dot(a.z, c.w));
123         t.w = make_float4(dot(a.w, c.x), dot(a.w, c.y), dot(a.w, c.z), dot(a.w, c.w));
124
125         return t;
126 }
127
128 __device_inline Transform make_transform(float a, float b, float c, float d,
129                                                                         float e, float f, float g, float h,
130                                                                         float i, float j, float k, float l,
131                                                                         float m, float n, float o, float p)
132 {
133         Transform t;
134
135         t.x.x = a; t.x.y = b; t.x.z = c; t.x.w = d;
136         t.y.x = e; t.y.y = f; t.y.z = g; t.y.w = h;
137         t.z.x = i; t.z.y = j; t.z.z = k; t.z.w = l;
138         t.w.x = m; t.w.y = n; t.w.z = o; t.w.w = p;
139
140         return t;
141 }
142
143 __device_inline Transform transform_translate(float3 t)
144 {
145         return make_transform(
146                 1, 0, 0, t.x,
147                 0, 1, 0, t.y,
148                 0, 0, 1, t.z,
149                 0, 0, 0, 1);
150 }
151
152 __device_inline Transform transform_translate(float x, float y, float z)
153 {
154         return transform_translate(make_float3(x, y, z));
155 }
156
157 __device_inline Transform transform_scale(float3 s)
158 {
159         return make_transform(
160                 s.x, 0, 0, 0,
161                 0, s.y, 0, 0,
162                 0, 0, s.z, 0,
163                 0, 0, 0, 1);
164 }
165
166 __device_inline Transform transform_scale(float x, float y, float z)
167 {
168         return transform_scale(make_float3(x, y, z));
169 }
170
171 __device_inline Transform transform_perspective(float fov, float n, float f)
172 {
173         Transform persp = make_transform(
174                 1, 0, 0, 0,
175                 0, 1, 0, 0,
176                 0, 0, f / (f - n), -f*n / (f - n),
177                 0, 0, 1, 0);
178
179         float inv_angle = 1.0f/tanf(0.5f*fov);
180
181         Transform scale = transform_scale(inv_angle, inv_angle, 1);
182
183         return scale * persp;
184 }
185
186 __device_inline Transform transform_rotate(float angle, float3 axis)
187 {
188         float s = sinf(angle);
189         float c = cosf(angle);
190         float t = 1.f - c;
191
192         axis = normalize(axis);
193
194         return make_transform(
195                 axis.x*axis.x*t + c,
196                 axis.x*axis.y*t - s*axis.z,
197                 axis.x*axis.z*t + s*axis.y,
198                 0.0f,
199
200                 axis.y*axis.x*t + s*axis.z,
201                 axis.y*axis.y*t + c,
202                 axis.y*axis.z*t - s*axis.x,
203                 0.0f,
204
205                 axis.z*axis.x*t - s*axis.y,
206                 axis.z*axis.y*t + s*axis.x,
207                 axis.z*axis.z*t + c,
208                 0.0f,
209
210                 0.0f, 0.0f, 0.0f, 1.0f);
211 }
212
213 __device_inline Transform transform_euler(float3 euler)
214 {
215         return
216                 transform_rotate(euler.x, make_float3(1.0f, 0.0f, 0.0f)) *
217                 transform_rotate(euler.y, make_float3(0.0f, 1.0f, 0.0f)) *
218                 transform_rotate(euler.z, make_float3(0.0f, 0.0f, 1.0f));
219 }
220
221 __device_inline Transform transform_orthographic(float znear, float zfar)
222 {
223         return transform_scale(1.0f, 1.0f, 1.0f / (zfar-znear)) *
224                 transform_translate(0.0f, 0.0f, -znear);
225 }
226
227 __device_inline Transform transform_identity()
228 {
229         return transform_scale(1.0f, 1.0f, 1.0f);
230 }
231
232 __device_inline bool operator==(const Transform& A, const Transform& B)
233 {
234         return memcmp(&A, &B, sizeof(Transform)) == 0;
235 }
236
237 __device_inline bool operator!=(const Transform& A, const Transform& B)
238 {
239         return !(A == B);
240 }
241
242 __device_inline float3 transform_get_column(const Transform *t, int column)
243 {
244         return make_float3(t->x[column], t->y[column], t->z[column]);
245 }
246
247 __device_inline void transform_set_column(Transform *t, int column, float3 value)
248 {
249         t->x[column] = value.x;
250         t->y[column] = value.y;
251         t->z[column] = value.z;
252 }
253
254 Transform transform_inverse(const Transform& a);
255
256 __device_inline bool transform_uniform_scale(const Transform& tfm, float& scale)
257 {
258         /* the epsilon here is quite arbitrary, but this function is only used for
259          * surface area and bump, where we except it to not be so sensitive */
260         Transform ttfm = transform_transpose(tfm);
261         float eps = 1e-6f; 
262         
263         float sx = len_squared(float4_to_float3(tfm.x));
264         float sy = len_squared(float4_to_float3(tfm.y));
265         float sz = len_squared(float4_to_float3(tfm.z));
266         float stx = len_squared(float4_to_float3(ttfm.x));
267         float sty = len_squared(float4_to_float3(ttfm.y));
268         float stz = len_squared(float4_to_float3(ttfm.z));
269
270         if(fabsf(sx - sy) < eps && fabsf(sx - sz) < eps &&
271            fabsf(sx - stx) < eps && fabsf(sx - sty) < eps &&
272            fabsf(sx - stz) < eps) {
273                 scale = sx;
274                 return true;
275         }
276    
277    return false;
278 }
279
280 __device_inline bool transform_negative_scale(const Transform& tfm)
281 {
282         float3 c0 = transform_get_column(&tfm, 0);
283         float3 c1 = transform_get_column(&tfm, 1);
284         float3 c2 = transform_get_column(&tfm, 2);
285
286         return (dot(cross(c0, c1), c2) < 0.0f);
287 }
288
289 __device_inline Transform transform_clear_scale(const Transform& tfm)
290 {
291         Transform ntfm = tfm;
292
293         transform_set_column(&ntfm, 0, normalize(transform_get_column(&ntfm, 0)));
294         transform_set_column(&ntfm, 1, normalize(transform_get_column(&ntfm, 1)));
295         transform_set_column(&ntfm, 2, normalize(transform_get_column(&ntfm, 2)));
296
297         return ntfm;
298 }
299
300 #endif
301
302 /* Motion Transform */
303
304 __device_inline float4 quat_interpolate(float4 q1, float4 q2, float t)
305 {
306         float costheta = dot(q1, q2);
307
308         /* rotate around shortest angle */
309         if(costheta < 0.0f) {
310                 costheta = -costheta;
311                 q1 = -q1;
312         }
313
314         if(costheta > 0.9995f) {
315                 /* linear interpolation in degenerate case */
316                 return normalize((1.0f - t)*q1 + t*q2);
317         }
318         else  {
319                 /* slerp */
320                 float theta = acosf(clamp(costheta, -1.0f, 1.0f));
321                 float thetap = theta * t;
322                 float4 qperp = normalize(q2 - q1 * costheta);
323                 return q1 * cosf(thetap) + qperp * sinf(thetap);
324         }
325 }
326
327 __device_inline Transform transform_quick_inverse(Transform M)
328 {
329         Transform R;
330         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);
331
332         det = (det != 0.0f)? 1.0f/det: 0.0f;
333
334         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);
335         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);
336         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);
337         float3 T = -make_float3(M.x.w, M.y.w, M.z.w);
338
339         R.x = make_float4(Rx.x, Rx.y, Rx.z, dot(Rx, T));
340         R.y = make_float4(Ry.x, Ry.y, Ry.z, dot(Ry, T));
341         R.z = make_float4(Rz.x, Rz.y, Rz.z, dot(Rz, T));
342         R.w = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
343
344         return R;
345 }
346
347 __device_inline void transform_compose(Transform *tfm, const Transform *decomp)
348 {
349         /* rotation */
350         float q0, q1, q2, q3, qda, qdb, qdc, qaa, qab, qac, qbb, qbc, qcc;
351
352         q0 = M_SQRT2_F * decomp->x.w;
353         q1 = M_SQRT2_F * decomp->x.x;
354         q2 = M_SQRT2_F * decomp->x.y;
355         q3 = M_SQRT2_F * decomp->x.z;
356
357         qda = q0*q1;
358         qdb = q0*q2;
359         qdc = q0*q3;
360         qaa = q1*q1;
361         qab = q1*q2;
362         qac = q1*q3;
363         qbb = q2*q2;
364         qbc = q2*q3;
365         qcc = q3*q3;
366
367         float3 rotation_x = make_float3(1.0f-qbb-qcc, -qdc+qab, qdb+qac);
368         float3 rotation_y = make_float3(qdc+qab, 1.0f-qaa-qcc, -qda+qbc);
369         float3 rotation_z = make_float3(-qdb+qac, qda+qbc, 1.0f-qaa-qbb);
370
371         /* scale */
372         float3 scale_x = make_float3(decomp->y.w, decomp->z.z, decomp->w.y);
373         float3 scale_y = make_float3(decomp->z.x, decomp->z.w, decomp->w.z);
374         float3 scale_z = make_float3(decomp->z.y, decomp->w.x, decomp->w.w);
375
376         /* compose with translation */
377         tfm->x = make_float4(dot(rotation_x, scale_x), dot(rotation_x, scale_y), dot(rotation_x, scale_z), decomp->y.x);
378         tfm->y = make_float4(dot(rotation_y, scale_x), dot(rotation_y, scale_y), dot(rotation_y, scale_z), decomp->y.y);
379         tfm->z = make_float4(dot(rotation_z, scale_x), dot(rotation_z, scale_y), dot(rotation_z, scale_z), decomp->y.z);
380         tfm->w = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
381 }
382
383 __device void transform_motion_interpolate(Transform *tfm, const MotionTransform *motion, float t)
384 {
385         Transform decomp;
386
387         /* 3 point bezier curve interpolation for position */
388         float3 Ppre = float4_to_float3(motion->pre.y);
389         float3 Pmid = float4_to_float3(motion->mid.y);
390         float3 Ppost = float4_to_float3(motion->post.y);
391
392         float3 Pcontrol = 2.0f*Pmid - 0.5f*(Ppre + Ppost);
393         float3 P = Ppre*t*t + Pcontrol*2.0f*t*(1.0f - t) + Ppost*(1.0f - t)*(1.0f - t);
394
395         decomp.y.x = P.x;
396         decomp.y.y = P.y;
397         decomp.y.z = P.z;
398
399         /* linear interpolation for rotation and scale */
400         if(t < 0.5f) {
401                 t *= 2.0f;
402
403                 decomp.x = quat_interpolate(motion->pre.x, motion->mid.x, t);
404                 decomp.y.w = (1.0f - t)*motion->pre.y.w + t*motion->mid.y.w;
405                 decomp.z = (1.0f - t)*motion->pre.z + t*motion->mid.z;
406                 decomp.w = (1.0f - t)*motion->pre.w + t*motion->mid.w;
407         }
408         else {
409                 t = (t - 0.5f)*2.0f;
410
411                 decomp.x = quat_interpolate(motion->mid.x, motion->post.x, t);
412                 decomp.y.w = (1.0f - t)*motion->mid.y.w + t*motion->post.y.w;
413                 decomp.z = (1.0f - t)*motion->mid.z + t*motion->post.z;
414                 decomp.w = (1.0f - t)*motion->mid.w + t*motion->post.w;
415         }
416
417         /* compose rotation, translation, scale into matrix */
418         transform_compose(tfm, &decomp);
419 }
420
421 #ifndef __KERNEL_GPU__
422
423 __device_inline bool operator==(const MotionTransform& A, const MotionTransform& B)
424 {
425         return (A.pre == B.pre && A.post == B.post);
426 }
427
428 void transform_motion_decompose(MotionTransform *decomp, const MotionTransform *motion, const Transform *mid);
429
430 #endif
431
432 CCL_NAMESPACE_END
433
434 #endif /* __UTIL_TRANSFORM_H__ */
435