Fix #32974: cycles curved motion blur is not working well combined with rotation,

problem is that the curved interpolation is not constant speed which leads to
mismatches. Turns out this is really hard to solve and implement efficiently, so
curved motion blur is disabled for now, until I can find a solution.

2 files changed, 34 insertions, 10 deletions

diff --git a/intern/cycles/util/util_transform.cpp b/intern/cycles/util/util_transform.cpp
index 70ee13d96d7..4eee024990f 100644
--- a/intern/cycles/util/util_transform.cpp
+++ b/intern/cycles/util/util_transform.cpp
@@ -251,6 +251,13 @@ void transform_motion_decompose(MotionTransform *decomp, const MotionTransform *
 	transform_decompose(&decomp->pre, &motion->pre);
 	transform_decompose(&decomp->mid, mid);
 	transform_decompose(&decomp->post, &motion->post);
+
+	/* ensure rotation around shortest angle, negated quaternions are the same
+	 * but this means we don't have to do the check in quat_interpolate */
+	if(dot(decomp->mid.x, decomp->post.x) < 0.0f)
+		decomp->mid.x = -decomp->mid.x;
+	if(dot(decomp->pre.x, decomp->mid.x) < 0.0f)
+		decomp->pre.x = -decomp->pre.x;
 }
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/util/util_transform.h b/intern/cycles/util/util_transform.h
index df525542207..65162ebf4e6 100644
--- a/intern/cycles/util/util_transform.h
+++ b/intern/cycles/util/util_transform.h
@@ -39,16 +39,16 @@ typedef struct Transform {
 #endif
 } Transform;
 
+/* 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 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) */
-
 /* Functions */
 
 __device_inline float3 transform_perspective(const Transform *t, const float3 a)
@@ -303,13 +303,11 @@ __device_inline Transform transform_clear_scale(const Transform& tfm)
 
 __device_inline float4 quat_interpolate(float4 q1, float4 q2, float t)
 {
+	/* 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);
 
-	/* rotate around shortest angle */
-	if(costheta < 0.0f) {
-		costheta = -costheta;
-		q1 = -q1;
-	}
+	/* possible optimization: it might be possible to precompute theta/qperp */
 
 	if(costheta > 0.9995f) {
 		/* linear interpolation in degenerate case */
@@ -318,14 +316,17 @@ __device_inline float4 quat_interpolate(float4 q1, float4 q2, float t)
 	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);
 	}
 }
 
 __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);
 
@@ -380,10 +381,17 @@ __device_inline void transform_compose(Transform *tfm, const Transform *decomp)
 	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 MotionTransform *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);
@@ -395,13 +403,18 @@ __device void transform_motion_interpolate(Transform *tfm, const MotionTransform
 	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
 		decomp.z = (1.0f - t)*motion->pre.z + t*motion->mid.z;
 		decomp.w = (1.0f - t)*motion->pre.w + t*motion->mid.w;
 	}
@@ -409,7 +422,11 @@ __device void transform_motion_interpolate(Transform *tfm, const MotionTransform
 		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 = (1.0f - t)*motion->mid.z + t*motion->post.z;
 		decomp.w = (1.0f - t)*motion->mid.w + t*motion->post.w;
 	}