Cycles: Improved robustness of hair motion blur.motion_curve_fix

In some instances, the number of control vertices of a hair could change mid-frame.
Cycles would then be unable to calculate proper motion blur for those hairs. This adds
interpolated CVs to fill in for the missing data. While this will not necessarily result in
a fully accurate reconstruction of the guide hair, it preserves motion blur instead of disabling it.

Reviewers: #cycles, sergey

Reviewed By: #cycles, sergey

Subscribers: sergey, brecht, #cycles

Tags: #cycles

Differential Revision: https://developer.blender.org/D3695

1 files changed, 83 insertions, 16 deletions

diff --git a/intern/cycles/kernel/kernel_montecarlo.h b/intern/cycles/kernel/kernel_montecarlo.h
index 9b96bb80c32..dde93844dd3 100644
--- a/intern/cycles/kernel/kernel_montecarlo.h
+++ b/intern/cycles/kernel/kernel_montecarlo.h
@@ -187,7 +187,10 @@ ccl_device float2 regular_polygon_sample(float corners, float rotation, float u,
 ccl_device float3 ensure_valid_reflection(float3 Ng, float3 I, float3 N)
 {
 	float3 R = 2*dot(N, I)*N - I;
-	if(dot(Ng, R) >= 0.05f) {
+
+	/* Reflection rays may always be at least as shallow as the incoming ray. */
+	float threshold = min(0.9f*dot(Ng, I), 0.01f);
+	if(dot(Ng, R) >= threshold) {
 		return N;
 	}
 
@@ -195,24 +198,88 @@ ccl_device float3 ensure_valid_reflection(float3 Ng, float3 I, float3 N)
 	 * The X axis is found by normalizing the component of N that's orthogonal to Ng.
 	 * The Y axis isn't actually needed.
 	 */
-	float3 X = normalize(N - dot(N, Ng)*Ng);
-
-	/* Calculate N.z and N.x in the local coordinate system. */
-	float Iz = dot(I, Ng);
-	float Ix2 = sqr(dot(I, X)), Iz2 = sqr(Iz);
-	float Ix2Iz2 = Ix2 + Iz2;
-
-	float a = safe_sqrtf(Ix2*(Ix2Iz2 - sqr(0.05f)));
-	float b = Iz*0.05f + Ix2Iz2;
-	float c = (a + b > 0.0f)? (a + b) : (-a + b);
+	float NdotNg = dot(N, Ng);
+	float3 X = normalize(N - NdotNg*Ng);
+
+	/* Calculate N.z and N.x in the local coordinate system.
+	 *
+	 * The goal of this computation is to find a N' that is rotated towards Ng just enough
+	 * to lift R' above the threshold (here called t), therefore dot(R', Ng) = t.
+	 *
+	 * According to the standard reflection equation, this means that we want dot(2*dot(N', I)*N' - I, Ng) = t.
+	 *
+	 * Since the Z axis of our local coordinate system is Ng, dot(x, Ng) is just x.z, so we get 2*dot(N', I)*N'.z - I.z = t.
+	 *
+	 * The rotation is simple to express in the coordinate system we formed - since N lies in the X-Z-plane, we know that
+	 * N' will also lie in the X-Z-plane, so N'.y = 0 and therefore dot(N', I) = N'.x*I.x + N'.z*I.z .
+	 *
+	 * Furthermore, we want N' to be normalized, so N'.x = sqrt(1 - N'.z^2).
+	 *
+	 * With these simplifications, we get the final equation 2*(sqrt(1 - N'.z^2)*I.x + N'.z*I.z)*N'.z - I.z = t.
+	 *
+	 * The only unknown here is N'.z, so we can solve for that.
+	 *
+	 * The equation has four solutions in general:
+	 *
+	 * N'.z = +-sqrt(0.5*(+-sqrt(I.x^2*(I.x^2 + I.z^2 - t^2)) + t*I.z + I.x^2 + I.z^2)/(I.x^2 + I.z^2))
+	 * We can simplify this expression a bit by grouping terms:
+	 *
+	 * a = I.x^2 + I.z^2
+	 * b = sqrt(I.x^2 * (a - t^2))
+	 * c = I.z*t + a
+	 * N'.z = +-sqrt(0.5*(+-b + c)/a)
+	 *
+	 * Two solutions can immediately be discarded because they're negative so N' would lie in the lower hemisphere.
+	 */
+	float Ix = dot(I, X), Iz = dot(I, Ng);
+	float Ix2 = sqr(Ix), Iz2 = sqr(Iz);
+	float a = Ix2 + Iz2;
+
+	float b = safe_sqrtf(Ix2*(a - sqr(threshold)));
+	float c = Iz*threshold + a;
+
+	/* Evaluate both solutions.
+	 * In many cases one can be immediately discarded (if N'.z would be imaginary or larger than one), so check for that first.
+	 * If no option is viable (might happen in extreme cases like N being in the wrong hemisphere), give up and return Ng. */
+	float fac = 0.5f/a;
+	float N1_z2 = fac*(b+c), N2_z2 = fac*(-b+c);
+	bool valid1 = (N1_z2 > 1e-5f) && (N1_z2 <= (1.0f + 1e-5f));
+	bool valid2 = (N2_z2 > 1e-5f) && (N2_z2 <= (1.0f + 1e-5f));
+
+	float2 N_new;
+	if(valid1 && valid2) {
+		/* If both are possible, do the expensive reflection-based check. */
+		float2 N1 = make_float2(safe_sqrtf(1.0f - N1_z2), safe_sqrtf(N1_z2));
+		float2 N2 = make_float2(safe_sqrtf(1.0f - N2_z2), safe_sqrtf(N2_z2));
+
+		float R1 = 2*(N1.x*Ix + N1.y*Iz)*N1.y - Iz;
+		float R2 = 2*(N2.x*Ix + N2.y*Iz)*N2.y - Iz;
+
+		valid1 = (R1 >= 1e-5f);
+		valid2 = (R2 >= 1e-5f);
+		if(valid1 && valid2) {
+			/* If both solutions are valid, return the one with the shallower reflection since it will be closer to the input
+			 * (if the original reflection wasn't shallow, we would not be in this part of the function). */
+			N_new = (R1 < R2)? N1 : N2;
+		}
+		else {
+			/* If only one reflection is valid (= positive), pick that one. */
+			N_new = (R1 > R2)? N1 : N2;
+		}
 
-	float Nz = safe_sqrtf(0.5f * c * (1.0f / Ix2Iz2));
-	float Nx = safe_sqrtf(1.0f - sqr(Nz));
+	}
+	else if(valid1 || valid2) {
+		/* Only one solution passes the N'.z criterium, so pick that one. */
+		float Nz2 = valid1? N1_z2 : N2_z2;
+		N_new = make_float2(safe_sqrtf(1.0f - Nz2), safe_sqrtf(Nz2));
+	}
+	else {
+		return Ng;
+	}
 
-	/* Transform back into global coordinates. */
-	return Nx*X + Nz*Ng;
+	return N_new.x*X + N_new.y*Ng;
 }
 
 CCL_NAMESPACE_END
 
-#endif /* __KERNEL_MONTECARLO_CL__ */
+#endif  /* __KERNEL_MONTECARLO_CL__ */