Merged revision(s) 58452-58584 from trunk/blender into soc-2013-dingto.
[blender-staging.git] / intern / cycles / util / util_transform.h
index e4897ee6787e9aa5ca80b47387fe21c938f74127..66801e90b56d601e6759063da71d3efd924a6463 100644 (file)
@@ -39,14 +39,23 @@ typedef struct Transform {
 #endif
 } Transform;
 
-typedef struct MotionTransform {
+/* transform decomposed in rotation/translation/scale. we use the same data
+ * structure as Transform, and tightly pack decomposition into it. first the
+ * rotation (4), then translation (3), then 3x3 scale matrix (9).
+ *
+ * For the DecompMotionTransform we drop scale from pre/post. */
+
+typedef struct __may_alias MotionTransform {
        Transform pre;
+       Transform mid;
        Transform post;
 } MotionTransform;
 
-/* transform decomposed in rotation/translation/scale. we use the same data
- * structure as Transform, and tightly pack decomposition into it. first the
- * rotation (4), then translation (3), then 3x3 scale matrix (9) */
+typedef struct DecompMotionTransform {
+       Transform mid;
+       float4 pre_x, pre_y;
+       float4 post_x, post_y;
+} DecompMotionTransform;
 
 /* Functions */
 
@@ -61,16 +70,20 @@ __device_inline float3 transform_perspective(const Transform *t, const float3 a)
 
 __device_inline float3 transform_point(const Transform *t, const float3 a)
 {
-       float4 b = make_float4(a.x, a.y, a.z, 1.0f);
-       float3 c = make_float3(dot(t->x, b), dot(t->y, b), dot(t->z, b));
+       float3 c = make_float3(
+               a.x*t->x.x + a.y*t->x.y + a.z*t->x.z + t->x.w,
+               a.x*t->y.x + a.y*t->y.y + a.z*t->y.z + t->y.w,
+               a.x*t->z.x + a.y*t->z.y + a.z*t->z.z + t->z.w);
 
        return c;
 }
 
 __device_inline float3 transform_direction(const Transform *t, const float3 a)
 {
-       float4 b = make_float4(a.x, a.y, a.z, 0.0f);
-       float3 c = make_float3(dot(t->x, b), dot(t->y, b), dot(t->z, b));
+       float3 c = make_float3(
+               a.x*t->x.x + a.y*t->x.y + a.z*t->x.z,
+               a.x*t->y.x + a.y*t->y.y + a.z*t->y.z,
+               a.x*t->z.x + a.y*t->z.y + a.z*t->z.z);
 
        return c;
 }
@@ -84,17 +97,6 @@ __device_inline float3 transform_direction_transposed(const Transform *t, const
        return make_float3(dot(x, a), dot(y, a), dot(z, a));
 }
 
-#ifndef __KERNEL_GPU__
-
-__device_inline void print_transform(const char *label, const Transform& t)
-{
-       print_float4(label, t.x);
-       print_float4(label, t.y);
-       print_float4(label, t.z);
-       print_float4(label, t.w);
-       printf("\n");
-}
-
 __device_inline Transform transform_transpose(const Transform a)
 {
        Transform t;
@@ -107,6 +109,23 @@ __device_inline Transform transform_transpose(const Transform a)
        return t;
 }
 
+__device_inline Transform make_transform(float a, float b, float c, float d,
+                                                                       float e, float f, float g, float h,
+                                                                       float i, float j, float k, float l,
+                                                                       float m, float n, float o, float p)
+{
+       Transform t;
+
+       t.x.x = a; t.x.y = b; t.x.z = c; t.x.w = d;
+       t.y.x = e; t.y.y = f; t.y.z = g; t.y.w = h;
+       t.z.x = i; t.z.y = j; t.z.z = k; t.z.w = l;
+       t.w.x = m; t.w.y = n; t.w.z = o; t.w.w = p;
+
+       return t;
+}
+
+#ifndef __KERNEL_GPU__
+
 __device_inline Transform operator*(const Transform a, const Transform b)
 {
        Transform c = transform_transpose(b);
@@ -120,19 +139,13 @@ __device_inline Transform operator*(const Transform a, const Transform b)
        return t;
 }
 
-__device_inline Transform make_transform(float a, float b, float c, float d,
-                                                                       float e, float f, float g, float h,
-                                                                       float i, float j, float k, float l,
-                                                                       float m, float n, float o, float p)
+__device_inline void print_transform(const char *label, const Transform& t)
 {
-       Transform t;
-
-       t.x.x = a; t.x.y = b; t.x.z = c; t.x.w = d;
-       t.y.x = e; t.y.y = f; t.y.z = g; t.y.w = h;
-       t.z.x = i; t.z.y = j; t.z.z = k; t.z.w = l;
-       t.w.x = m; t.w.y = n; t.w.z = o; t.w.w = p;
-
-       return t;
+       print_float4(label, t.x);
+       print_float4(label, t.y);
+       print_float4(label, t.z);
+       print_float4(label, t.w);
+       printf("\n");
 }
 
 __device_inline Transform transform_translate(float3 t)
@@ -182,7 +195,7 @@ __device_inline Transform transform_rotate(float angle, float3 axis)
 {
        float s = sinf(angle);
        float c = cosf(angle);
-       float t = 1.f - c;
+       float t = 1.0f - c;
 
        axis = normalize(axis);
 
@@ -251,7 +264,7 @@ Transform transform_inverse(const Transform& a);
 __device_inline bool transform_uniform_scale(const Transform& tfm, float& scale)
 {
        /* the epsilon here is quite arbitrary, but this function is only used for
-          surface area and bump, where we except it to not be so sensitive */
+        * surface area and bump, where we except it to not be so sensitive */
        Transform ttfm = transform_transpose(tfm);
        float eps = 1e-6f; 
        
@@ -298,21 +311,35 @@ __device_inline Transform transform_clear_scale(const Transform& tfm)
 
 __device_inline float4 quat_interpolate(float4 q1, float4 q2, float t)
 {
+       /* use simpe nlerp instead of slerp. it's faster and almost the same */
+       return normalize((1.0f - t)*q1 + t*q2);
+
+#if 0
+       /* note: this does not ensure rotation around shortest angle, q1 and q2
+        * are assumed to be matched already in transform_motion_decompose */
        float costheta = dot(q1, q2);
 
+       /* possible optimization: it might be possible to precompute theta/qperp */
+
        if(costheta > 0.9995f) {
+               /* linear interpolation in degenerate case */
                return normalize((1.0f - t)*q1 + t*q2);
        }
        else  {
+               /* slerp */
                float theta = acosf(clamp(costheta, -1.0f, 1.0f));
-               float thetap = theta * t;
                float4 qperp = normalize(q2 - q1 * costheta);
+               float thetap = theta * t;
                return q1 * cosf(thetap) + qperp * sinf(thetap);
        }
+#endif
 }
 
 __device_inline Transform transform_quick_inverse(Transform M)
 {
+       /* possible optimization: can we avoid doing this altogether and construct
+        * the inverse matrix directly from negated translation, transposed rotation,
+        * scale can be inverted but what about shearing? */
        Transform R;
        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);
 
@@ -367,15 +394,56 @@ __device_inline void transform_compose(Transform *tfm, const Transform *decomp)
        tfm->w = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
 }
 
-__device void transform_motion_interpolate(Transform *tfm, const MotionTransform *motion, float t)
+/* Disabled for now, need arc-length parametrization for constant speed motion.
+ * #define CURVED_MOTION_INTERPOLATE */
+
+__device void transform_motion_interpolate(Transform *tfm, const DecompMotionTransform *motion, float t)
 {
+       /* possible optimization: is it worth it adding a check to skip scaling?
+        * it's probably quite uncommon to have scaling objects. or can we skip
+        * just shearing perhaps? */
        Transform decomp;
 
-       decomp.x = quat_interpolate(motion->pre.x, motion->post.x, t);
-       decomp.y = (1.0f - t)*motion->pre.y + t*motion->post.y;
-       decomp.z = (1.0f - t)*motion->pre.z + t*motion->post.z;
-       decomp.w = (1.0f - t)*motion->pre.w + t*motion->post.w;
+#ifdef CURVED_MOTION_INTERPOLATE
+       /* 3 point bezier curve interpolation for position */
+       float3 Ppre = float4_to_float3(motion->pre_y);
+       float3 Pmid = float4_to_float3(motion->mid.y);
+       float3 Ppost = float4_to_float3(motion->post_y);
+
+       float3 Pcontrol = 2.0f*Pmid - 0.5f*(Ppre + Ppost);
+       float3 P = Ppre*t*t + Pcontrol*2.0f*t*(1.0f - t) + Ppost*(1.0f - t)*(1.0f - t);
+
+       decomp.y.x = P.x;
+       decomp.y.y = P.y;
+       decomp.y.z = P.z;
+#endif
+
+       /* linear interpolation for rotation and scale */
+       if(t < 0.5f) {
+               t *= 2.0f;
+
+               decomp.x = quat_interpolate(motion->pre_x, motion->mid.x, t);
+#ifdef CURVED_MOTION_INTERPOLATE
+               decomp.y.w = (1.0f - t)*motion->pre_y.w + t*motion->mid.y.w;
+#else
+               decomp.y = (1.0f - t)*motion->pre_y + t*motion->mid.y;
+#endif
+       }
+       else {
+               t = (t - 0.5f)*2.0f;
+
+               decomp.x = quat_interpolate(motion->mid.x, motion->post_x, t);
+#ifdef CURVED_MOTION_INTERPOLATE
+               decomp.y.w = (1.0f - t)*motion->mid.y.w + t*motion->post_y.w;
+#else
+               decomp.y = (1.0f - t)*motion->mid.y + t*motion->post_y;
+#endif
+       }
+
+       decomp.z = motion->mid.z;
+       decomp.w = motion->mid.w;
 
+       /* compose rotation, translation, scale into matrix */
        transform_compose(tfm, &decomp);
 }
 
@@ -386,7 +454,8 @@ __device_inline bool operator==(const MotionTransform& A, const MotionTransform&
        return (A.pre == B.pre && A.post == B.post);
 }
 
-void transform_motion_decompose(MotionTransform *decomp, const MotionTransform *motion);
+float4 transform_to_quat(const Transform& tfm);
+void transform_motion_decompose(DecompMotionTransform *decomp, const MotionTransform *motion, const Transform *mid);
 
 #endif