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