Merged revision(s) 58452-58584 from trunk/blender into soc-2013-dingto.
[blender-staging.git] / intern / cycles / util / util_transform.h
index c43736fb2e402b50802f050537f93ef4af345c22..66801e90b56d601e6759063da71d3efd924a6463 100644 (file)
@@ -28,6 +28,8 @@
 
 CCL_NAMESPACE_BEGIN
 
+/* Data Types */
+
 typedef struct Transform {
        float4 x, y, z, w; /* rows */
 
@@ -37,31 +39,62 @@ typedef struct Transform {
 #endif
 } Transform;
 
-__device_inline float3 transform(const Transform *t, const float3 a)
+/* 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;
+
+typedef struct DecompMotionTransform {
+       Transform mid;
+       float4 pre_x, pre_y;
+       float4 post_x, post_y;
+} DecompMotionTransform;
+
+/* Functions */
+
+__device_inline float3 transform_perspective(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));
+       float w = dot(t->w, b);
 
-       return c/dot(t->w, b);
+       return (w != 0.0f)? c/w: make_float3(0.0f, 0.0f, 0.0f);
 }
 
-__device_inline float3 transform_direction(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, 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 + 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;
 }
 
-#ifndef __KERNEL_GPU__
+__device_inline float3 transform_direction(const Transform *t, const float3 a)
+{
+       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);
 
-__device_inline void print_transform(const char *label, const Transform& t)
+       return c;
+}
+
+__device_inline float3 transform_direction_transposed(const Transform *t, const float3 a)
 {
-       print_float4(label, t.x);
-       print_float4(label, t.y);
-       print_float4(label, t.z);
-       print_float4(label, t.w);
-       printf("\n");
+       float3 x = make_float3(t->x.x, t->y.x, t->z.x);
+       float3 y = make_float3(t->x.y, t->y.y, t->z.y);
+       float3 z = make_float3(t->x.z, t->y.z, t->z.z);
+
+       return make_float3(dot(x, a), dot(y, a), dot(z, a));
 }
 
 __device_inline Transform transform_transpose(const Transform a)
@@ -76,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);
@@ -89,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)
@@ -151,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);
 
@@ -208,22 +252,29 @@ __device_inline float3 transform_get_column(const Transform *t, int column)
        return make_float3(t->x[column], t->y[column], t->z[column]);
 }
 
+__device_inline void transform_set_column(Transform *t, int column, float3 value)
+{
+       t->x[column] = value.x;
+       t->y[column] = value.y;
+       t->z[column] = value.z;
+}
+
 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-7f; 
-       
-       float sx = len(float4_to_float3(tfm.x));
-       float sy = len(float4_to_float3(tfm.y));
-       float sz = len(float4_to_float3(tfm.z));
-       float stx = len(float4_to_float3(ttfm.x));
-       float sty = len(float4_to_float3(ttfm.y));
-       float stz = len(float4_to_float3(ttfm.z));
+       float eps = 1e-6f; 
        
+       float sx = len_squared(float4_to_float3(tfm.x));
+       float sy = len_squared(float4_to_float3(tfm.y));
+       float sz = len_squared(float4_to_float3(tfm.z));
+       float stx = len_squared(float4_to_float3(ttfm.x));
+       float sty = len_squared(float4_to_float3(ttfm.y));
+       float stz = len_squared(float4_to_float3(ttfm.z));
+
        if(fabsf(sx - sy) < eps && fabsf(sx - sz) < eps &&
           fabsf(sx - stx) < eps && fabsf(sx - sty) < eps &&
           fabsf(sx - stz) < eps) {
@@ -234,6 +285,178 @@ __device_inline bool transform_uniform_scale(const Transform& tfm, float& scale)
    return false;
 }
 
+__device_inline bool transform_negative_scale(const Transform& tfm)
+{
+       float3 c0 = transform_get_column(&tfm, 0);
+       float3 c1 = transform_get_column(&tfm, 1);
+       float3 c2 = transform_get_column(&tfm, 2);
+
+       return (dot(cross(c0, c1), c2) < 0.0f);
+}
+
+__device_inline Transform transform_clear_scale(const Transform& tfm)
+{
+       Transform ntfm = tfm;
+
+       transform_set_column(&ntfm, 0, normalize(transform_get_column(&ntfm, 0)));
+       transform_set_column(&ntfm, 1, normalize(transform_get_column(&ntfm, 1)));
+       transform_set_column(&ntfm, 2, normalize(transform_get_column(&ntfm, 2)));
+
+       return ntfm;
+}
+
+#endif
+
+/* Motion Transform */
+
+__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));
+               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);
+
+       det = (det != 0.0f)? 1.0f/det: 0.0f;
+
+       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);
+       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);
+       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);
+       float3 T = -make_float3(M.x.w, M.y.w, M.z.w);
+
+       R.x = make_float4(Rx.x, Rx.y, Rx.z, dot(Rx, T));
+       R.y = make_float4(Ry.x, Ry.y, Ry.z, dot(Ry, T));
+       R.z = make_float4(Rz.x, Rz.y, Rz.z, dot(Rz, T));
+       R.w = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
+
+       return R;
+}
+
+__device_inline void transform_compose(Transform *tfm, const Transform *decomp)
+{
+       /* rotation */
+       float q0, q1, q2, q3, qda, qdb, qdc, qaa, qab, qac, qbb, qbc, qcc;
+
+       q0 = M_SQRT2_F * decomp->x.w;
+       q1 = M_SQRT2_F * decomp->x.x;
+       q2 = M_SQRT2_F * decomp->x.y;
+       q3 = M_SQRT2_F * decomp->x.z;
+
+       qda = q0*q1;
+       qdb = q0*q2;
+       qdc = q0*q3;
+       qaa = q1*q1;
+       qab = q1*q2;
+       qac = q1*q3;
+       qbb = q2*q2;
+       qbc = q2*q3;
+       qcc = q3*q3;
+
+       float3 rotation_x = make_float3(1.0f-qbb-qcc, -qdc+qab, qdb+qac);
+       float3 rotation_y = make_float3(qdc+qab, 1.0f-qaa-qcc, -qda+qbc);
+       float3 rotation_z = make_float3(-qdb+qac, qda+qbc, 1.0f-qaa-qbb);
+
+       /* scale */
+       float3 scale_x = make_float3(decomp->y.w, decomp->z.z, decomp->w.y);
+       float3 scale_y = make_float3(decomp->z.x, decomp->z.w, decomp->w.z);
+       float3 scale_z = make_float3(decomp->z.y, decomp->w.x, decomp->w.w);
+
+       /* compose with translation */
+       tfm->x = make_float4(dot(rotation_x, scale_x), dot(rotation_x, scale_y), dot(rotation_x, scale_z), decomp->y.x);
+       tfm->y = make_float4(dot(rotation_y, scale_x), dot(rotation_y, scale_y), dot(rotation_y, scale_z), decomp->y.y);
+       tfm->z = make_float4(dot(rotation_z, scale_x), dot(rotation_z, scale_y), dot(rotation_z, scale_z), decomp->y.z);
+       tfm->w = make_float4(0.0f, 0.0f, 0.0f, 1.0f);
+}
+
+/* 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;
+
+#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);
+}
+
+#ifndef __KERNEL_GPU__
+
+__device_inline bool operator==(const MotionTransform& A, const MotionTransform& B)
+{
+       return (A.pre == B.pre && A.post == B.post);
+}
+
+float4 transform_to_quat(const Transform& tfm);
+void transform_motion_decompose(DecompMotionTransform *decomp, const MotionTransform *motion, const Transform *mid);
+
 #endif
 
 CCL_NAMESPACE_END