Merged changes in the trunk up to revision 47700.
[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 post;
45 } MotionTransform;
46
47 /* transform decomposed in rotation/translation/scale. we use the same data
48  * structure as Transform, and tightly pack decomposition into it. first the
49  * rotation (4), then translation (3), then 3x3 scale matrix (9) */
50
51 /* Functions */
52
53 __device_inline float3 transform_perspective(const Transform *t, const float3 a)
54 {
55         float4 b = make_float4(a.x, a.y, a.z, 1.0f);
56         float3 c = make_float3(dot(t->x, b), dot(t->y, b), dot(t->z, b));
57         float w = dot(t->w, b);
58
59         return (w != 0.0f)? c/w: make_float3(0.0f, 0.0f, 0.0f);
60 }
61
62 __device_inline float3 transform_point(const Transform *t, const float3 a)
63 {
64         float3 c = make_float3(
65                 a.x*t->x.x + a.y*t->x.y + a.z*t->x.z + t->x.w,
66                 a.x*t->y.x + a.y*t->y.y + a.z*t->y.z + t->y.w,
67                 a.x*t->z.x + a.y*t->z.y + a.z*t->z.z + t->z.w);
68
69         return c;
70 }
71
72 __device_inline float3 transform_direction(const Transform *t, const float3 a)
73 {
74         float3 c = make_float3(
75                 a.x*t->x.x + a.y*t->x.y + a.z*t->x.z,
76                 a.x*t->y.x + a.y*t->y.y + a.z*t->y.z,
77                 a.x*t->z.x + a.y*t->z.y + a.z*t->z.z);
78
79         return c;
80 }
81
82 __device_inline float3 transform_direction_transposed(const Transform *t, const float3 a)
83 {
84         float3 x = make_float3(t->x.x, t->y.x, t->z.x);
85         float3 y = make_float3(t->x.y, t->y.y, t->z.y);
86         float3 z = make_float3(t->x.z, t->y.z, t->z.z);
87
88         return make_float3(dot(x, a), dot(y, a), dot(z, a));
89 }
90
91 #ifndef __KERNEL_GPU__
92
93 __device_inline void print_transform(const char *label, const Transform& t)
94 {
95         print_float4(label, t.x);
96         print_float4(label, t.y);
97         print_float4(label, t.z);
98         print_float4(label, t.w);
99         printf("\n");
100 }
101
102 __device_inline Transform transform_transpose(const Transform a)
103 {
104         Transform t;
105
106         t.x.x = a.x.x; t.x.y = a.y.x; t.x.z = a.z.x; t.x.w = a.w.x;
107         t.y.x = a.x.y; t.y.y = a.y.y; t.y.z = a.z.y; t.y.w = a.w.y;
108         t.z.x = a.x.z; t.z.y = a.y.z; t.z.z = a.z.z; t.z.w = a.w.z;
109         t.w.x = a.x.w; t.w.y = a.y.w; t.w.z = a.z.w; t.w.w = a.w.w;
110
111         return t;
112 }
113
114 __device_inline Transform operator*(const Transform a, const Transform b)
115 {
116         Transform c = transform_transpose(b);
117         Transform t;
118
119         t.x = make_float4(dot(a.x, c.x), dot(a.x, c.y), dot(a.x, c.z), dot(a.x, c.w));
120         t.y = make_float4(dot(a.y, c.x), dot(a.y, c.y), dot(a.y, c.z), dot(a.y, c.w));
121         t.z = make_float4(dot(a.z, c.x), dot(a.z, c.y), dot(a.z, c.z), dot(a.z, c.w));
122         t.w = make_float4(dot(a.w, c.x), dot(a.w, c.y), dot(a.w, c.z), dot(a.w, c.w));
123
124         return t;
125 }
126
127 __device_inline Transform make_transform(float a, float b, float c, float d,
128                                                                         float e, float f, float g, float h,
129                                                                         float i, float j, float k, float l,
130                                                                         float m, float n, float o, float p)
131 {
132         Transform t;
133
134         t.x.x = a; t.x.y = b; t.x.z = c; t.x.w = d;
135         t.y.x = e; t.y.y = f; t.y.z = g; t.y.w = h;
136         t.z.x = i; t.z.y = j; t.z.z = k; t.z.w = l;
137         t.w.x = m; t.w.y = n; t.w.z = o; t.w.w = p;
138
139         return t;
140 }
141
142 __device_inline Transform transform_translate(float3 t)
143 {
144         return make_transform(
145                 1, 0, 0, t.x,
146                 0, 1, 0, t.y,
147                 0, 0, 1, t.z,
148                 0, 0, 0, 1);
149 }
150
151 __device_inline Transform transform_translate(float x, float y, float z)
152 {
153         return transform_translate(make_float3(x, y, z));
154 }
155
156 __device_inline Transform transform_scale(float3 s)
157 {
158         return make_transform(
159                 s.x, 0, 0, 0,
160                 0, s.y, 0, 0,
161                 0, 0, s.z, 0,
162                 0, 0, 0, 1);
163 }
164
165 __device_inline Transform transform_scale(float x, float y, float z)
166 {
167         return transform_scale(make_float3(x, y, z));
168 }
169
170 __device_inline Transform transform_perspective(float fov, float n, float f)
171 {
172         Transform persp = make_transform(
173                 1, 0, 0, 0,
174                 0, 1, 0, 0,
175                 0, 0, f / (f - n), -f*n / (f - n),
176                 0, 0, 1, 0);
177
178         float inv_angle = 1.0f/tanf(0.5f*fov);
179
180         Transform scale = transform_scale(inv_angle, inv_angle, 1);
181
182         return scale * persp;
183 }
184
185 __device_inline Transform transform_rotate(float angle, float3 axis)
186 {
187         float s = sinf(angle);
188         float c = cosf(angle);
189         float t = 1.f - c;
190
191         axis = normalize(axis);
192
193         return make_transform(
194                 axis.x*axis.x*t + c,
195                 axis.x*axis.y*t - s*axis.z,
196                 axis.x*axis.z*t + s*axis.y,
197                 0.0f,
198
199                 axis.y*axis.x*t + s*axis.z,
200                 axis.y*axis.y*t + c,
201                 axis.y*axis.z*t - s*axis.x,
202                 0.0f,
203
204                 axis.z*axis.x*t - s*axis.y,
205                 axis.z*axis.y*t + s*axis.x,
206                 axis.z*axis.z*t + c,
207                 0.0f,
208
209                 0.0f, 0.0f, 0.0f, 1.0f);
210 }
211
212 __device_inline Transform transform_euler(float3 euler)
213 {
214         return
215                 transform_rotate(euler.x, make_float3(1.0f, 0.0f, 0.0f)) *
216                 transform_rotate(euler.y, make_float3(0.0f, 1.0f, 0.0f)) *
217                 transform_rotate(euler.z, make_float3(0.0f, 0.0f, 1.0f));
218 }
219
220 __device_inline Transform transform_orthographic(float znear, float zfar)
221 {
222         return transform_scale(1.0f, 1.0f, 1.0f / (zfar-znear)) *
223                 transform_translate(0.0f, 0.0f, -znear);
224 }
225
226 __device_inline Transform transform_identity()
227 {
228         return transform_scale(1.0f, 1.0f, 1.0f);
229 }
230
231 __device_inline bool operator==(const Transform& A, const Transform& B)
232 {
233         return memcmp(&A, &B, sizeof(Transform)) == 0;
234 }
235
236 __device_inline bool operator!=(const Transform& A, const Transform& B)
237 {
238         return !(A == B);
239 }
240
241 __device_inline float3 transform_get_column(const Transform *t, int column)
242 {
243         return make_float3(t->x[column], t->y[column], t->z[column]);
244 }
245
246 __device_inline void transform_set_column(Transform *t, int column, float3 value)
247 {
248         t->x[column] = value.x;
249         t->y[column] = value.y;
250         t->z[column] = value.z;
251 }
252
253 Transform transform_inverse(const Transform& a);
254
255 __device_inline bool transform_uniform_scale(const Transform& tfm, float& scale)
256 {
257         /* the epsilon here is quite arbitrary, but this function is only used for
258          * surface area and bump, where we except it to not be so sensitive */
259         Transform ttfm = transform_transpose(tfm);
260         float eps = 1e-6f; 
261         
262         float sx = len_squared(float4_to_float3(tfm.x));
263         float sy = len_squared(float4_to_float3(tfm.y));
264         float sz = len_squared(float4_to_float3(tfm.z));
265         float stx = len_squared(float4_to_float3(ttfm.x));
266         float sty = len_squared(float4_to_float3(ttfm.y));
267         float stz = len_squared(float4_to_float3(ttfm.z));
268
269         if(fabsf(sx - sy) < eps && fabsf(sx - sz) < eps &&
270            fabsf(sx - stx) < eps && fabsf(sx - sty) < eps &&
271            fabsf(sx - stz) < eps) {
272                 scale = sx;
273                 return true;
274         }
275    
276    return false;
277 }
278
279 __device_inline bool transform_negative_scale(const Transform& tfm)
280 {
281         float3 c0 = transform_get_column(&tfm, 0);
282         float3 c1 = transform_get_column(&tfm, 1);
283         float3 c2 = transform_get_column(&tfm, 2);
284
285         return (dot(cross(c0, c1), c2) < 0.0f);
286 }
287
288 __device_inline Transform transform_clear_scale(const Transform& tfm)
289 {
290         Transform ntfm = tfm;
291
292         transform_set_column(&ntfm, 0, normalize(transform_get_column(&ntfm, 0)));
293         transform_set_column(&ntfm, 1, normalize(transform_get_column(&ntfm, 1)));
294         transform_set_column(&ntfm, 2, normalize(transform_get_column(&ntfm, 2)));
295
296         return ntfm;
297 }
298
299 #endif
300
301 /* Motion Transform */
302
303 __device_inline float4 quat_interpolate(float4 q1, float4 q2, float t)
304 {
305         float costheta = dot(q1, q2);
306
307         if(costheta > 0.9995f) {
308                 return normalize((1.0f - t)*q1 + t*q2);
309         }
310         else  {
311                 float theta = acosf(clamp(costheta, -1.0f, 1.0f));
312                 float thetap = theta * t;
313                 float4 qperp = normalize(q2 - q1 * costheta);
314                 return q1 * cosf(thetap) + qperp * sinf(thetap);
315         }
316 }
317
318 __device_inline Transform transform_quick_inverse(Transform M)
319 {
320         Transform R;
321         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);
322
323         det = (det != 0.0f)? 1.0f/det: 0.0f;
324
325         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);
326         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);
327         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);
328         float3 T = -make_float3(M.x.w, M.y.w, M.z.w);
329
330         R.x = make_float4(Rx.x, Rx.y, Rx.z, dot(Rx, T));
331         R.y = make_float4(Ry.x, Ry.y, Ry.z, dot(Ry, T));
332         R.z = make_float4(Rz.x, Rz.y, Rz.z, dot(Rz, T));
333         R.w = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
334
335         return R;
336 }
337
338 __device_inline void transform_compose(Transform *tfm, const Transform *decomp)
339 {
340         /* rotation */
341         float q0, q1, q2, q3, qda, qdb, qdc, qaa, qab, qac, qbb, qbc, qcc;
342
343         q0 = M_SQRT2_F * decomp->x.w;
344         q1 = M_SQRT2_F * decomp->x.x;
345         q2 = M_SQRT2_F * decomp->x.y;
346         q3 = M_SQRT2_F * decomp->x.z;
347
348         qda = q0*q1;
349         qdb = q0*q2;
350         qdc = q0*q3;
351         qaa = q1*q1;
352         qab = q1*q2;
353         qac = q1*q3;
354         qbb = q2*q2;
355         qbc = q2*q3;
356         qcc = q3*q3;
357
358         float3 rotation_x = make_float3(1.0f-qbb-qcc, -qdc+qab, qdb+qac);
359         float3 rotation_y = make_float3(qdc+qab, 1.0f-qaa-qcc, -qda+qbc);
360         float3 rotation_z = make_float3(-qdb+qac, qda+qbc, 1.0f-qaa-qbb);
361
362         /* scale */
363         float3 scale_x = make_float3(decomp->y.w, decomp->z.z, decomp->w.y);
364         float3 scale_y = make_float3(decomp->z.x, decomp->z.w, decomp->w.z);
365         float3 scale_z = make_float3(decomp->z.y, decomp->w.x, decomp->w.w);
366
367         /* compose with translation */
368         tfm->x = make_float4(dot(rotation_x, scale_x), dot(rotation_x, scale_y), dot(rotation_x, scale_z), decomp->y.x);
369         tfm->y = make_float4(dot(rotation_y, scale_x), dot(rotation_y, scale_y), dot(rotation_y, scale_z), decomp->y.y);
370         tfm->z = make_float4(dot(rotation_z, scale_x), dot(rotation_z, scale_y), dot(rotation_z, scale_z), decomp->y.z);
371         tfm->w = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
372 }
373
374 __device void transform_motion_interpolate(Transform *tfm, const MotionTransform *motion, float t)
375 {
376         Transform decomp;
377
378         decomp.x = quat_interpolate(motion->pre.x, motion->post.x, t);
379         decomp.y = (1.0f - t)*motion->pre.y + t*motion->post.y;
380         decomp.z = (1.0f - t)*motion->pre.z + t*motion->post.z;
381         decomp.w = (1.0f - t)*motion->pre.w + t*motion->post.w;
382
383         transform_compose(tfm, &decomp);
384 }
385
386 #ifndef __KERNEL_GPU__
387
388 __device_inline bool operator==(const MotionTransform& A, const MotionTransform& B)
389 {
390         return (A.pre == B.pre && A.post == B.post);
391 }
392
393 void transform_motion_decompose(MotionTransform *decomp, const MotionTransform *motion);
394
395 #endif
396
397 CCL_NAMESPACE_END
398
399 #endif /* __UTIL_TRANSFORM_H__ */
400