[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [54648] trunk/blender/source/blender/ blenkernel: Fix for [#34099] Particles leaking from moving meshes

[Alex Fraser]

Revision: 54648
          http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=54648
Author:   z0r
Date:     2013-02-19 02:24:52 +0000 (Tue, 19 Feb 2013)
Log Message:
-----------
Fix for [#34099] Particles leaking from moving meshes

This was caused by a floating point precision error. During collision detection, Newton-Raphson iteration is used to find the exact time of the collision. But when using subframes, the initial Newton step was too small. Now the initial step is given in absolute units. When subframes = 0, this should behave almost the same as before.

Thanks to Janne Karhu, Lukas Toenne and Ton Roosendaal for their help with this patch, and to AutoCRC for funding.

Modified Paths:
--------------
    trunk/blender/source/blender/blenkernel/BKE_particle.h
    trunk/blender/source/blender/blenkernel/intern/particle_system.c

Modified: trunk/blender/source/blender/blenkernel/BKE_particle.h
===================================================================
--- trunk/blender/source/blender/blenkernel/BKE_particle.h	2013-02-19 00:41:13 UTC (rev 54647)
+++ trunk/blender/source/blender/blenkernel/BKE_particle.h	2013-02-19 02:24:52 UTC (rev 54648)
@@ -220,7 +220,10 @@
 
 	ParticleCollisionElement pce;
 
-	float total_time, inv_timestep;
+	/* total_time is the amount of time in this subframe
+	 * inv_total_time is the opposite
+	 * inv_timestep is the inverse of the amount of time in this frame */
+	float total_time, inv_total_time, inv_timestep;
 
 	float radius;
 	float co1[3], co2[3];

Modified: trunk/blender/source/blender/blenkernel/intern/particle_system.c
===================================================================
--- trunk/blender/source/blender/blenkernel/intern/particle_system.c	2013-02-19 00:41:13 UTC (rev 54647)
+++ trunk/blender/source/blender/blenkernel/intern/particle_system.c	2013-02-19 02:24:52 UTC (rev 54648)
@@ -3028,13 +3028,20 @@
 	normalize_qt(pa->state.rot);
 }
 
-/************************************************/
-/*			Collisions							*/
-/************************************************/
+/************************************************
+ *			Collisions
+ *
+ * The algorithm is roughly:
+ *  1. Use a BVH tree to search for faces that a particle may collide with.
+ *  2. Use Newton's method to find the exact time at which the collision occurs.
+ *     http://en.wikipedia.org/wiki/Newton's_method
+ *
+ ************************************************/
 #define COLLISION_MAX_COLLISIONS	10
 #define COLLISION_MIN_RADIUS 0.001f
 #define COLLISION_MIN_DISTANCE 0.0001f
 #define COLLISION_ZERO 0.00001f
+#define COLLISION_INIT_STEP 0.00008f
 typedef float (*NRDistanceFunc)(float *p, float radius, ParticleCollisionElement *pce, float *nor);
 static float nr_signed_distance_to_plane(float *p, float radius, ParticleCollisionElement *pce, float *nor)
 {
@@ -3189,16 +3196,20 @@
 /* find first root in range [0-1] starting from 0 */
 static float collision_newton_rhapson(ParticleCollision *col, float radius, ParticleCollisionElement *pce, NRDistanceFunc distance_func)
 {
-	float t0, t1, d0, d1, dd, n[3];
+	float t0, t1, dt_init, d0, d1, dd, n[3];
 	int iter;
 
 	pce->inv_nor = -1;
 
+	/* Initial step size should be small, but not too small or floating point
+	 * precision errors will appear. - z0r */
+	dt_init = COLLISION_INIT_STEP * col->inv_total_time;
+
 	/* start from the beginning */
 	t0 = 0.f;
 	collision_interpolate_element(pce, t0, col->f, col);
 	d0 = distance_func(col->co1, radius, pce, n);
-	t1 = 0.001f;
+	t1 = dt_init;
 	d1 = 0.f;
 
 	for (iter=0; iter<10; iter++) {//, itersum++) {
@@ -3208,11 +3219,6 @@
 
 		d1 = distance_func(pce->p, radius, pce, n);
 
-		/* no movement, so no collision */
-		if (d1 == d0) {
-			return -1.f;
-		}
-
 		/* particle already inside face, so report collision */
 		if (iter == 0 && d0 < 0.f && d0 > -radius) {
 			copy_v3_v3(pce->p, col->co1);
@@ -3220,7 +3226,24 @@
 			pce->inside = 1;
 			return 0.f;
 		}
-		
+
+		/* Zero gradient (no movement relative to element). Can't step from
+		 * here. */
+		if (d1 == d0) {
+			/* If first iteration, try from other end where the gradient may be
+			 * greater. Note: code duplicated below. */
+			if (iter == 0) {
+				t0 = 1.f;
+				collision_interpolate_element(pce, t0, col->f, col);
+				d0 = distance_func(col->co2, radius, pce, n);
+				t1 = 1.0f - dt_init;
+				d1 = 0.f;
+				continue;
+			}
+			else
+				return -1.f;
+		}
+
 		dd = (t1-t0)/(d1-d0);
 
 		t0 = t1;
@@ -3228,14 +3251,14 @@
 
 		t1 -= d1*dd;
 
-		/* particle movin away from plane could also mean a strangely rotating face, so check from end */
+		/* Particle moving away from plane could also mean a strangely rotating
+		 * face, so check from end. Note: code duplicated above. */
 		if (iter == 0 && t1 < 0.f) {
 			t0 = 1.f;
 			collision_interpolate_element(pce, t0, col->f, col);
 			d0 = distance_func(col->co2, radius, pce, n);
-			t1 = 0.999f;
+			t1 = 1.0f - dt_init;
 			d1 = 0.f;
-
 			continue;
 		}
 		else if (iter == 1 && (t1 < -COLLISION_ZERO || t1 > 1.f))
@@ -3683,6 +3706,7 @@
 	memset(&col, 0, sizeof(ParticleCollision));
 
 	col.total_time = timestep * dfra;
+	col.inv_total_time = 1.0f/col.total_time;
 	col.inv_timestep = 1.0f/timestep;
 
 	col.cfra = cfra;