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