Fixed non_recursive BVHbuild with openmp

CHanged the BENCH functions to use: gettimeofday (wall time) instead of clock (cpu time)
This was to test if the openmp was working right.

4 files changed, 102 insertions, 32 deletions

diff --git a/source/blender/blenkernel/BKE_deform.h b/source/blender/blenkernel/BKE_deform.h
index b0570e8bd45..73a9b2b5d4e 100644
--- a/source/blender/blenkernel/BKE_deform.h
+++ b/source/blender/blenkernel/BKE_deform.h
@@ -46,8 +46,8 @@ int get_defgroup_num (struct Object *ob, struct bDeformGroup *dg);
 int get_named_vertexgroup_num (Object *ob, char *name);
 void unique_vertexgroup_name (struct bDeformGroup *dg, struct Object *ob);
 
-float deformvert_get_weight(struct MDeformVert *dvert, int group_num);
-float vertexgroup_get_vertex_weight(struct MDeformVert *dvert, int index, int group_num);
+float deformvert_get_weight(const struct MDeformVert *dvert, int group_num);
+float vertexgroup_get_vertex_weight(const struct MDeformVert *dvert, int index, int group_num);
 
 #endif
 
diff --git a/source/blender/blenkernel/intern/deform.c b/source/blender/blenkernel/intern/deform.c
index 888432005d0..3143c5e4df2 100644
--- a/source/blender/blenkernel/intern/deform.c
+++ b/source/blender/blenkernel/intern/deform.c
@@ -221,11 +221,11 @@ void unique_vertexgroup_name (bDeformGroup *dg, Object *ob)
 	}	
 }
 
-float deformvert_get_weight(struct MDeformVert *dvert, int group_num)
+float deformvert_get_weight(const struct MDeformVert *dvert, int group_num)
 {
 	if(dvert)
 	{
-		MDeformWeight *dw = dvert->dw;
+		const MDeformWeight *dw = dvert->dw;
 		int i;
 
 		for(i=dvert->totweight; i>0; i--, dw++)
@@ -237,7 +237,7 @@ float deformvert_get_weight(struct MDeformVert *dvert, int group_num)
 	return 0.0;
 }
 
-float vertexgroup_get_vertex_weight(struct MDeformVert *dvert, int index, int group_num)
+float vertexgroup_get_vertex_weight(const struct MDeformVert *dvert, int index, int group_num)
 {
 	if(group_num == -1)
 		return 1.0;
diff --git a/source/blender/blenkernel/intern/shrinkwrap.c b/source/blender/blenkernel/intern/shrinkwrap.c
index f53c7c30390..292bdbebcc6 100644
--- a/source/blender/blenkernel/intern/shrinkwrap.c
+++ b/source/blender/blenkernel/intern/shrinkwrap.c
@@ -62,6 +62,8 @@
 /* Benchmark macros */
 #if 1
 
+
+#if 0
 #define BENCH(a)	\
 	do {			\
 		clock_t _clock_init = clock();	\
@@ -76,6 +78,30 @@
 #define BENCH_REPORT(name)	printf("%s: %fms\n", TO_STR(name), JOIN(_bench_total,name)*1000.0f/CLOCKS_PER_SEC)
 
 #else
+#include <sys/time.h>
+
+#define BENCH(a)	\
+	do {			\
+		double _t1, _t2;				\
+		struct timeval _tstart, _tend;	\
+		gettimeofday ( &_tstart, NULL);	\
+		(a);							\
+		gettimeofday ( &_tend, NULL);	\
+		_t1 = ( double ) _tstart.tv_sec + ( double ) _tstart.tv_usec/ ( 1000*1000 );	\
+		_t2 = ( double )   _tend.tv_sec + ( double )   _tend.tv_usec/ ( 1000*1000 );	\
+		printf("%s: %fms\n", #a, _t2-_t1);\
+	} while(0)
+
+#define BENCH_VAR(name)		clock_t JOIN(_bench_step,name) = 0, JOIN(_bench_total,name) = 0
+#define BENCH_BEGIN(name)	JOIN(_bench_step, name) = clock()
+#define BENCH_END(name)		JOIN(_bench_total,name) += clock() - JOIN(_bench_step,name)
+#define BENCH_RESET(name)	JOIN(_bench_total, name) = 0
+#define BENCH_REPORT(name)	printf("%s: %fms\n", TO_STR(name), JOIN(_bench_total,name)*1000.0f/CLOCKS_PER_SEC)
+
+
+#endif
+
+#else
 
 #define BENCH(a)	(a)
 #define BENCH_VAR(name)
@@ -1019,7 +1045,7 @@ void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, Object *ob, DerivedM
 		switch(smd->shrinkType)
 		{
 			case MOD_SHRINKWRAP_NEAREST_SURFACE:
-				shrinkwrap_calc_nearest_surface_point(&calc);
+				BENCH(shrinkwrap_calc_nearest_surface_point(&calc));
 			break;
 
 			case MOD_SHRINKWRAP_NORMAL:
@@ -1027,7 +1053,7 @@ void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, Object *ob, DerivedM
 			break;
 
 			case MOD_SHRINKWRAP_NEAREST_VERTEX:
-				shrinkwrap_calc_nearest_vertex(&calc);
+				BENCH(shrinkwrap_calc_nearest_vertex(&calc));
 			break;
 		}
 
@@ -1047,7 +1073,6 @@ void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
 {
 	int i;
 	int vgroup		= get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
-	float *co;
 
 	BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
 	BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
@@ -1055,15 +1080,17 @@ void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
 	MDeformVert *dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
 
 
-	bvhtree_from_mesh_verts(&treeData, calc->target, 0.0, 2, 6);
+	BENCH(bvhtree_from_mesh_verts(&treeData, calc->target, 0.0, 2, 6));
 	if(treeData.tree == NULL) return OUT_OF_MEMORY();
 
 	//Setup nearest
 	nearest.index = -1;
 	nearest.dist = FLT_MAX;
 
-	for(co = calc->vertexCos[i=0]; i<calc->numVerts; co = calc->vertexCos[++i])
+//#pragma omp parallel for private(i) private(nearest) schedule(static)
+	for(i = 0; i<calc->numVerts; ++i)
 	{
+		float *co = calc->vertexCos[i];
 		int index;
 		float tmp_co[3];
 		float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
@@ -1165,7 +1192,6 @@ void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
 	int i;
 	int vgroup		= get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
 	char use_normal = calc->smd->shrinkOpts;
-	float *co;
 
 	//setup raytracing
 	BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
@@ -1189,7 +1215,7 @@ void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
 
 	CDDM_calc_normals(calc->original);	//Normals maybe arent yet calculated
 
-	bvhtree_from_mesh_faces(&treeData, calc->target, calc->keptDist, 4, 6);
+	BENCH(bvhtree_from_mesh_faces(&treeData, calc->target, calc->keptDist, 4, 6));
 	if(treeData.tree == NULL) return OUT_OF_MEMORY();
 
 
@@ -1200,13 +1226,15 @@ void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
 		//TODO currently we need a copy in case object_get_derived_final returns an emDM that does not defines getVertArray or getFace array
 		limit_mesh = CDDM_copy( object_get_derived_final(calc->smd->cutPlane, CD_MASK_BAREMESH) );
 		if(limit_mesh)
-			bvhtree_from_mesh_faces(&limitData, limit_mesh, 0.0, 4, 6);
+			BENCH(bvhtree_from_mesh_faces(&limitData, limit_mesh, 0.0, 4, 6));
 		else
 			printf("CutPlane finalDerived mesh is null\n");
 	}
 
-	for(co = calc->vertexCos[i=0]; i<calc->numVerts; co = calc->vertexCos[++i])
+//#pragma omp parallel for private(i) private(hit) schedule(static)
+	for(i = 0; i<calc->numVerts; ++i)
 	{
+		float *co = calc->vertexCos[i];
 		float tmp_co[3], tmp_no[3];
 		float lim = 1000;		//TODO: we should use FLT_MAX here, but sweepsphere code isnt prepared for that
 		float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
@@ -1274,17 +1302,17 @@ void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
 void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
 {
 	int i;
-	int vgroup		= get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
-	float *co;
+
+	const int vgroup		 = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
+	const MDeformVert *dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
 
 	BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
 	BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
 
-	MDeformVert *dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
 
 
 	//Create a bvh-tree of the given target
-	bvhtree_from_mesh_faces( &treeData, calc->target, 0.0, 2, 6);
+	BENCH(bvhtree_from_mesh_faces( &treeData, calc->target, 0.0, 2, 6));
 	if(treeData.tree == NULL) return OUT_OF_MEMORY();
 
 	//Setup nearest
@@ -1293,8 +1321,10 @@ void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
 
 
 	//Find the nearest vertex 
-	for(co = calc->vertexCos[i=0]; i<calc->numVerts; co = calc->vertexCos[++i])
+//#pragma omp parallel for private(i) private(nearest) schedule(static)
+	for(i = 0; i<calc->numVerts; ++i)
 	{
+		float *co = calc->vertexCos[i];
 		int index;
 		float tmp_co[3];
 		float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
diff --git a/source/blender/blenlib/intern/BLI_kdopbvh.c b/source/blender/blenlib/intern/BLI_kdopbvh.c
index 9b36051e59f..fcac5df934e 100644
--- a/source/blender/blenlib/intern/BLI_kdopbvh.c
+++ b/source/blender/blenlib/intern/BLI_kdopbvh.c
@@ -632,7 +632,22 @@ static void omp_bvh_div_nodes(BVHTree *tree, BVHNode *node, int start, int end,
 }
 
 
+static void print_tree(BVHTree *tree, BVHNode *node, int depth)
+{
+	int i;
+	for(i=0; i<depth; i++) printf(" ");
+	printf(" - %d (%d): ", node->index, node - tree->nodearray);
+	for(i=2*tree->start_axis; i<2*tree->stop_axis; i++)
+		printf("%.3f ", node->bv[i]);
+	printf("\n");
+
+	for(i=0; i<tree->tree_type; i++)
+		if(node->children[i])
+			print_tree(tree, node->children[i], depth+1);
+}
+
 #if 0
+
 static void verify_tree(BVHTree *tree)
 {
 	int i, j, check = 0;
@@ -711,10 +726,8 @@ static void build_implicit_tree_helper(BVHTree *tree, BVHBuildHelper *data)
 		data->leafs_per_child[0] *= data->tree_type
 	);
 
-
 	data->branches_on_level[0] = 1;
 
-
 	//We could stop the loop first (but I am lazy to find out when)
 	for(depth = 1; depth < 32; depth++)
 	{
@@ -731,13 +744,15 @@ static void build_implicit_tree_helper(BVHTree *tree, BVHBuildHelper *data)
 static int implicit_leafs_index(BVHBuildHelper *data, int depth, int child_index)
 {
 	int min_leaf_index = child_index * data->leafs_per_child[depth-1];
-	if(min_leaf_index < data->remain_leafs)
+	if(min_leaf_index <= data->remain_leafs)
 		return min_leaf_index;
-	else
+	else if(data->leafs_per_child[depth])
 		return data->totleafs - (data->branches_on_level[depth-1] - child_index) * data->leafs_per_child[depth];
+	else
+		return data->remain_leafs;
 }
 
-//WARNING: Beautifull/tricky code starts here :P
+//WARNING: Beautiful/tricky code starts here :P
 //Generalized implicit trees
 static void non_recursive_bvh_div_nodes(BVHTree *tree)
 {
@@ -748,7 +763,7 @@ static void non_recursive_bvh_div_nodes(BVHTree *tree)
 	const int num_leafs   = tree->totleaf;
 	const int num_branches= MAX2(1, (num_leafs + tree_type - 3) / (tree_type-1) );
 
-	BVHNode*  branches_array = tree->nodearray + tree->totleaf - 1; // This ocde uses 1 index arrays
+	BVHNode*  branches_array = tree->nodearray + tree->totleaf - 1; // This code uses 1 index arrays
 	BVHNode** leafs_array    = tree->nodes;
 
 	BVHBuildHelper data;
@@ -765,7 +780,7 @@ static void non_recursive_bvh_div_nodes(BVHTree *tree)
 
 		depth++;
 
-//#pragma omp parallel for private(j) schedule(static)
+#pragma omp parallel for private(j) schedule(static)
 		for(j = i; j < end_j; j++)
 		{
 			int k;
@@ -776,8 +791,11 @@ static void non_recursive_bvh_div_nodes(BVHTree *tree)
 			int parent_leafs_begin = implicit_leafs_index(&data, depth, parent_level_index);
 			int parent_leafs_end   = implicit_leafs_index(&data, depth, parent_level_index+1);
 
+			//split_axis = (depth*2 % 6); //use this instead of the 2 following lines for XYZ splitting
+
 			refit_kdop_hull(tree, parent, parent_leafs_begin, parent_leafs_end);
 			split_axis = get_largest_axis(parent->bv);
+
 			parent->main_axis = split_axis / 2;
 
 			for(k = 0; k < tree_type; k++)
@@ -792,33 +810,55 @@ static void non_recursive_bvh_div_nodes(BVHTree *tree)
 
 				if(child_leafs_end - child_leafs_begin > 1)
 				{
-					parent->children[k] = branches_array  + child_index;
+					parent->children[k] = branches_array + child_index;
+
+/*
+					printf("Add child %d (%d) to branch %d\n",
+								branches_array  + child_index - tree->nodearray,
+								branches_array[ child_index ].index,
+								parent - tree->nodearray
+						);
+*/
+
 					partition_nth_element(leafs_array, child_leafs_begin, parent_leafs_end, child_leafs_end, split_axis);
 				}
 				else if(child_leafs_end - child_leafs_begin == 1)
 				{
+/*
+					printf("Add child %d (%d) to branch %d\n",
+								leafs_array[ child_leafs_begin ] - tree->nodearray,
+								leafs_array[ child_leafs_begin ]->index,
+								parent - tree->nodearray
+						);
+*/
 					parent->children[k] = leafs_array[ child_leafs_begin ];
 				}
 				else
 				{
 					parent->children[k] = NULL;
+					break;
 				}
+				parent->totnode = k+1;
 			}
 		}
 	}
-	tree->nodes[tree->totleaf] = branches_array+0;
-	tree->totbranch = num_branches;
 
 
+	for(i = 0; i<num_branches; i++)
+		tree->nodes[tree->totleaf + i] = branches_array + 1 + i;
+
+	tree->totbranch = num_branches;
+
+//	BLI_bvhtree_update_tree(tree); //Uncoment this for XYZ splitting
 }
 
 void BLI_bvhtree_balance(BVHTree *tree)
 {
-	assert(tree->totbranch == 0);
-	assert(tree->totleaf   != 0);
+	if(tree->totleaf == 0) return;
 
+	assert(tree->totbranch == 0);
 	non_recursive_bvh_div_nodes(tree);
-	printf("here\n");
+
 /*
 	if(tree->totleaf != 0)
 	{