* Point Density texture

Replaced the previous KD-tree (for caching points) with a 
BVH-tree (thanks to Andre 'jaguarandi' Pinto for help here!).

The bvh is quite a bit faster and doesn't suffer some of the
artifacts that were apparent with the kd-tree.

I've also added a choice of falloff types: Standard, Smooth, and 
Sharp. Standard gives a harder edge, easier to see individual 
particles, and when used with a larger radius, Smooth and Sharp 
falloffs make a much cloudier appearance possible. See the image 
below (note the settings and render times too)

http://mke3.net/blender/devel/rendering/volumetrics/pointdensity_bvh.jpg

6 files changed, 183 insertions, 45 deletions

diff --git a/source/blender/blenkernel/intern/texture.c b/source/blender/blenkernel/intern/texture.c
index 45fd11faa32..ec34fcf4103 100644
--- a/source/blender/blenkernel/intern/texture.c
+++ b/source/blender/blenkernel/intern/texture.c
@@ -43,7 +43,7 @@
 #include "BLI_blenlib.h"
 #include "BLI_arithb.h"
 #include "BLI_rand.h"
-#include "BLI_kdtree.h"
+#include "BLI_kdopbvh.h"
 
 #include "DNA_texture_types.h"
 #include "DNA_key_types.h"
@@ -474,7 +474,7 @@ void default_tex(Tex *tex)
 
 	if (tex->pd) {
 		tex->pd->radius = 0.3f;
-		tex->pd->nearest = 5;
+		tex->pd->falloff_type = TEX_PD_FALLOFF_STD;
 	}
 
 	pit = tex->plugin;
@@ -874,9 +874,10 @@ PointDensity *BKE_add_pointdensity(void)
 	
 	pd= MEM_callocN(sizeof(PointDensity), "pointdensity");
 	pd->radius = 0.3f;
-	pd->nearest = 5;
+	pd->falloff_type = TEX_PD_FALLOFF_STD;
 	pd->source = TEX_PD_PSYS;
 	pd->point_tree = NULL;
+	//pd->point_data = NULL;
 	
 	return pd;
 } 
@@ -887,6 +888,7 @@ PointDensity *BKE_copy_pointdensity(PointDensity *pd)
 
 	pdn= MEM_dupallocN(pd);
 	pdn->point_tree = NULL;
+	//pdn->point_data = NULL;
 	
 	return pd;
 }
@@ -894,9 +896,15 @@ PointDensity *BKE_copy_pointdensity(PointDensity *pd)
 void BKE_free_pointdensitydata(PointDensity *pd)
 {
 	if (pd->point_tree) {
-		BLI_kdtree_free(pd->point_tree);
+		BLI_bvhtree_free(pd->point_tree);
 		pd->point_tree = NULL;
 	}
+	/*
+	if (pd->point_data) {
+		MEM_freeN(pd->point_data);
+		pd->point_data = NULL;
+	}
+	*/
 }
 
 void BKE_free_pointdensity(PointDensity *pd)
diff --git a/source/blender/blenlib/BLI_kdopbvh.h b/source/blender/blenlib/BLI_kdopbvh.h
index e3591a84e98..8912ac5bd13 100644
--- a/source/blender/blenlib/BLI_kdopbvh.h
+++ b/source/blender/blenlib/BLI_kdopbvh.h
@@ -71,6 +71,9 @@ typedef void (*BVHTree_NearestPointCallback) (void *userdata, int index, const f
 /* callback must update hit in case it finds a nearest successful hit */
 typedef void (*BVHTree_RayCastCallback) (void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit);
 
+/* callback to range search query */
+typedef void (*BVHTree_RangeQuery) (void *userdata, int index, float squared_dist, float radius);
+
 
 BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis);
 void BLI_bvhtree_free(BVHTree *tree);
@@ -93,5 +96,9 @@ int BLI_bvhtree_find_nearest(BVHTree *tree, const float *co, BVHTreeNearest *nea
 
 int BLI_bvhtree_ray_cast(BVHTree *tree, const float *co, const float *dir, float radius, BVHTreeRayHit *hit, BVHTree_RayCastCallback callback, void *userdata);
 
+/* range query */
+int BLI_bvhtree_range_query(BVHTree *tree, const float *co, float radius, BVHTree_RangeQuery callback, void *userdata);
+
+
 #endif // BLI_KDOPBVH_H
 
diff --git a/source/blender/blenlib/intern/BLI_kdopbvh.c b/source/blender/blenlib/intern/BLI_kdopbvh.c
index 30472beb3e6..f82b6b7f162 100644
--- a/source/blender/blenlib/intern/BLI_kdopbvh.c
+++ b/source/blender/blenlib/intern/BLI_kdopbvh.c
@@ -1171,7 +1171,7 @@ static float squared_dist(const float *a, const float *b)
 }
 
 //Determines the nearest point of the given node BV. Returns the squared distance to that point.
-static float calc_nearest_point(BVHNearestData *data, BVHNode *node, float *nearest)
+static float calc_nearest_point(const float *proj, BVHNode *node, float *nearest)
 {
 	int i;
 	const float *bv = node->bv;
@@ -1179,12 +1179,12 @@ static float calc_nearest_point(BVHNearestData *data, BVHNode *node, float *near
 	//nearest on AABB hull
 	for(i=0; i != 3; i++, bv += 2)
 	{
-		if(bv[0] > data->proj[i])
+		if(bv[0] > proj[i])
 			nearest[i] = bv[0];
-		else if(bv[1] < data->proj[i])
+		else if(bv[1] < proj[i])
 			nearest[i] = bv[1];
 		else
-			nearest[i] = data->proj[i];
+			nearest[i] = proj[i];
 	}
 
 /*
@@ -1206,7 +1206,7 @@ static float calc_nearest_point(BVHNearestData *data, BVHNode *node, float *near
 		}
 	}
 */
-	return squared_dist(data->co, nearest);
+	return squared_dist(proj, nearest);
 }
 
 
@@ -1229,7 +1229,7 @@ static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node)
 		else
 		{
 			data->nearest.index	= node->index;
-			data->nearest.dist	= calc_nearest_point(data, node, data->nearest.co);
+			data->nearest.dist	= calc_nearest_point(data->proj, node, data->nearest.co);
 		}
 	}
 	else
@@ -1243,7 +1243,7 @@ static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node)
 
 			for(i=0; i != node->totnode; i++)
 			{
-				if( calc_nearest_point(data, node->children[i], nearest) >= data->nearest.dist) continue;
+				if( calc_nearest_point(data->proj, node->children[i], nearest) >= data->nearest.dist) continue;
 				dfs_find_nearest_dfs(data, node->children[i]);
 			}
 		}
@@ -1251,7 +1251,7 @@ static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node)
 		{
 			for(i=node->totnode-1; i >= 0 ; i--)
 			{
-				if( calc_nearest_point(data, node->children[i], nearest) >= data->nearest.dist) continue;
+				if( calc_nearest_point(data->proj, node->children[i], nearest) >= data->nearest.dist) continue;
 				dfs_find_nearest_dfs(data, node->children[i]);
 			}
 		}
@@ -1261,7 +1261,7 @@ static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node)
 static void dfs_find_nearest_begin(BVHNearestData *data, BVHNode *node)
 {
 	float nearest[3], sdist;
-	sdist = calc_nearest_point(data, node, nearest);
+	sdist = calc_nearest_point(data->proj, node, nearest);
 	if(sdist >= data->nearest.dist) return;
 	dfs_find_nearest_dfs(data, node);
 }
@@ -1298,7 +1298,7 @@ static void bfs_find_nearest(BVHNearestData *data, BVHNode *node)
 	}
 
 	current.node = node;
-	current.dist = calc_nearest_point(data, node, nearest);
+	current.dist = calc_nearest_point(data->proj, node, nearest);
 
 	while(current.dist < data->nearest.dist)
 	{
@@ -1326,7 +1326,7 @@ static void bfs_find_nearest(BVHNearestData *data, BVHNode *node)
 				}
 
 				heap[heap_size].node = current.node->children[i];
-				heap[heap_size].dist = calc_nearest_point(data, current.node->children[i], nearest);
+				heap[heap_size].dist = calc_nearest_point(data->proj, current.node->children[i], nearest);
 
 				if(heap[heap_size].dist >= data->nearest.dist) continue;
 				heap_size++;
@@ -1524,3 +1524,90 @@ int BLI_bvhtree_ray_cast(BVHTree *tree, const float *co, const float *dir, float
 	return data.hit.index;
 }
 
+/*
+ * Range Query - as request by broken :P
+ *
+ * Allocs and fills an array with the indexs of node that are on the given spherical range (center, radius) 
+ * Returns the size of the array.
+ */
+typedef struct RangeQueryData
+{
+	BVHTree *tree;
+	const float *center;
+	float radius;			//squared radius
+
+	int hits;
+
+	BVHTree_RangeQuery callback;
+	void *userdata;
+
+
+} RangeQueryData;
+
+
+static void dfs_range_query(RangeQueryData *data, BVHNode *node)
+{
+	if(node->totnode == 0)
+	{
+
+		//Calculate the node min-coords (if the node was a point then this is the point coordinates)
+		float co[3];
+		co[0] = node->bv[0];
+		co[1] = node->bv[2];
+		co[2] = node->bv[4];
+
+	}
+	else
+	{
+		int i;
+		for(i=0; i != node->totnode; i++)
+		{
+			float nearest[3];
+			float dist = calc_nearest_point(data->center, node->children[i], nearest);
+			if(dist < data->radius)
+			{
+				//Its a leaf.. call the callback
+				if(node->children[i]->totnode == 0)
+				{
+					data->hits++;
+					data->callback( data->userdata, node->children[i]->index, dist, data->radius );
+				}
+				else
+					dfs_range_query( data, node->children[i] );
+			}
+		}
+	}
+}
+
+int BLI_bvhtree_range_query(BVHTree *tree, const float *co, float radius, BVHTree_RangeQuery callback, void *userdata)
+{
+	BVHNode * root = tree->nodes[tree->totleaf];
+
+	RangeQueryData data;
+	data.tree = tree;
+	data.center = co;
+	data.radius = radius*radius;
+	data.hits = 0;
+
+	data.callback = callback;
+	data.userdata = userdata;
+
+	if(root != NULL)
+	{
+		float nearest[3];
+		float dist = calc_nearest_point(data.center, root, nearest);
+		if(dist < data.radius)
+		{
+			//Its a leaf.. call the callback
+			if(root->totnode == 0)
+			{
+				data.hits++;
+				data.callback( data.userdata, root->index, dist, data.radius );
+			}
+			else
+				dfs_range_query( &data, root );
+		}
+	}
+
+	return data.hits;
+}
diff --git a/source/blender/makesdna/DNA_texture_types.h b/source/blender/makesdna/DNA_texture_types.h
index 52055f24f82..d23672f9284 100644
--- a/source/blender/makesdna/DNA_texture_types.h
+++ b/source/blender/makesdna/DNA_texture_types.h
@@ -130,7 +130,7 @@ typedef struct EnvMap {
 typedef struct PointDensity {
 	short flag;
 
-	short nearest;
+	short falloff_type;
 	float radius;
 
 	short source;
@@ -144,7 +144,10 @@ typedef struct PointDensity {
 	
 	short pdpad2;
 	
-	void *point_tree;		/* the kd-tree containing points */
+	void *point_tree;		/* the acceleration tree containing points */
+	//void *point_data;		/* dynamically allocated extra for extra information, like particle age */
+	//int pdpad3;
+	
 } PointDensity;
 
 typedef struct Tex {
@@ -415,10 +418,16 @@ typedef struct TexMapping {
 #define TEX_PD_OBJECT		1
 #define TEX_PD_FILE			2
 
+/* falloff_type */
+#define TEX_PD_FALLOFF_STD		0
+#define TEX_PD_FALLOFF_SMOOTH	1
+#define TEX_PD_FALLOFF_SHARP	2
+
 /* psys_cache_space */
 #define TEX_PD_OBJECTLOC	0
 #define TEX_PD_OBJECTSPACE	1
 #define TEX_PD_WORLDSPACE	2
 
+
 #endif
 
diff --git a/source/blender/render/intern/source/pointdensity.c b/source/blender/render/intern/source/pointdensity.c
index 17bc6bb58cf..4422b9fbbdd 100644
--- a/source/blender/render/intern/source/pointdensity.c
+++ b/source/blender/render/intern/source/pointdensity.c
@@ -27,7 +27,7 @@
 #include <stdio.h>
 
 #include "BLI_arithb.h"
-#include "BLI_kdtree.h"
+#include "BLI_kdopbvh.h"
 
 #include "BKE_DerivedMesh.h"
 #include "BKE_global.h"
@@ -71,7 +71,7 @@ static void pointdensity_cache_psys(Render *re, PointDensity *pd, Object *ob, Pa
 	/* in case ob->imat isn't up-to-date */
 	Mat4Invert(ob->imat, ob->obmat);
 	
-	pd->point_tree = BLI_kdtree_new(psys->totpart+psys->totchild);
+	pd->point_tree = BLI_bvhtree_new(psys->totpart+psys->totchild, 0.0, 2, 6);
 	
 	if (psys->totchild > 0 && !(psys->part->draw & PART_DRAW_PARENT))
 		childexists = 1;
@@ -93,11 +93,12 @@ static void pointdensity_cache_psys(Render *re, PointDensity *pd, Object *ob, Pa
 				/* TEX_PD_WORLDSPACE */
 			}
 			
-			BLI_kdtree_insert(pd->point_tree, i, partco, NULL);
+			BLI_bvhtree_insert(pd->point_tree, i, partco, 1);
 		}
 	}
 	
-	BLI_kdtree_balance(pd->point_tree);
+	BLI_bvhtree_balance(pd->point_tree);
+	
 	psys_render_restore(ob, psys);
 }
 
@@ -112,7 +113,7 @@ static void pointdensity_cache_object(Render *re, PointDensity *pd, ObjectRen *o
 	/* in case ob->imat isn't up-to-date */
 	Mat4Invert(obr->ob->imat, obr->ob->obmat);
 	
-	pd->point_tree = BLI_kdtree_new(obr->totvert);
+	pd->point_tree = BLI_bvhtree_new(obr->totvert, 0.0, 2, 6);
 	
 	for(i=0; i<obr->totvert; i++) {
 		float ver_co[3];
@@ -128,10 +129,10 @@ static void pointdensity_cache_object(Render *re, PointDensity *pd, ObjectRen *o
 			Mat4MulVecfl(re->viewinv, ver_co);
 		}
 		
-		BLI_kdtree_insert(pd->point_tree, i, ver_co, NULL);
+		BLI_bvhtree_insert(pd->point_tree, i, ver_co, 1);
 	}
 	
-	BLI_kdtree_balance(pd->point_tree);
+	BLI_bvhtree_balance(pd->point_tree);
 
 }
 static void cache_pointdensity(Render *re, Tex *tex)
@@ -139,7 +140,7 @@ static void cache_pointdensity(Render *re, Tex *tex)
 	PointDensity *pd = tex->pd;
 
 	if (pd->point_tree) {
-		BLI_kdtree_free(pd->point_tree);
+		BLI_bvhtree_free(pd->point_tree);
 		pd->point_tree = NULL;
 	}
 	
@@ -178,9 +179,16 @@ static void free_pointdensity(Render *re, Tex *tex)
 	PointDensity *pd = tex->pd;
 
 	if (pd->point_tree) {
-		BLI_kdtree_free(pd->point_tree);
+		BLI_bvhtree_free(pd->point_tree);
 		pd->point_tree = NULL;
 	}
+
+	/*
+	if (pd->point_data) {
+		MEM_freeN(pd->point_data);
+		pd->point_data = NULL;
+	}
+	*/
 }
 
 
@@ -216,33 +224,49 @@ void free_pointdensities(Render *re)
 	}
 }
 
+
+void accum_density_std(void *userdata, int index, float squared_dist, float squared_radius)
+{
+	float *density = userdata;
+	const float dist = squared_radius - squared_dist;
+	
+	*density+= dist;
+}
+
+void accum_density_smooth(void *userdata, int index, float squared_dist, float squared_radius)
+{
+	float *density = userdata;
+	const float dist = squared_radius - squared_dist;
+	
+	*density+= 3.0f*dist*dist - 2.0f*dist*dist*dist;
+}
+
+void accum_density_sharp(void *userdata, int index, float squared_dist, float squared_radius)
+{
+	float *density = userdata;
+	const float dist = squared_radius - squared_dist;
+	
+	*density+= dist*dist;
+}
+
 #define MAX_POINTS_NEAREST	25
 int pointdensitytex(Tex *tex, float *texvec, TexResult *texres)
 {
 	int rv = TEX_INT;
-	
 	PointDensity *pd = tex->pd;
-	KDTreeNearest nearest[MAX_POINTS_NEAREST];
 	float density=0.0f;
-	int n, neighbours=0;
 	
 	if ((!pd) || (!pd->point_tree)) {
 		texres->tin = 0.0f;
 		return 0;
 	}
 	
-	neighbours = BLI_kdtree_find_n_nearest(pd->point_tree, pd->nearest, texvec, NULL, nearest);
-	
-	for(n=1; n<neighbours; n++) {
-		if ( nearest[n].dist < pd->radius) {
-			float dist = 1.0 - (nearest[n].dist / pd->radius);
-			
-			density += 3.0f*dist*dist - 2.0f*dist*dist*dist;
-		}
-	}
-	
-	density /= neighbours;
-	density *= 1.0 / pd->radius;
+	if (pd->falloff_type == TEX_PD_FALLOFF_STD)
+		BLI_bvhtree_range_query(pd->point_tree, texvec, pd->radius, accum_density_std, &density);
+	else if (pd->falloff_type == TEX_PD_FALLOFF_SMOOTH)
+		BLI_bvhtree_range_query(pd->point_tree, texvec, pd->radius, accum_density_smooth, &density);
+	else if (pd->falloff_type == TEX_PD_FALLOFF_SHARP)
+		BLI_bvhtree_range_query(pd->point_tree, texvec, pd->radius, accum_density_sharp, &density);
 
 	texres->tin = density;
 
diff --git a/source/blender/src/buttons_shading.c b/source/blender/src/buttons_shading.c
index 85154682819..050dae6026e 100644
--- a/source/blender/src/buttons_shading.c
+++ b/source/blender/src/buttons_shading.c
@@ -751,12 +751,15 @@ static void texture_panel_pointdensity(Tex *tex)
 		uiDefBut(block, LABEL, B_NOP, "Density estimation:",
 			X2CLM1, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");
 
-		uiBlockBeginAlign(block);
 		uiDefButF(block, NUM, B_REDR, "Radius: ",
 			X2CLM1, yco-=BUTH, BUTW2, BUTH, &(pd->radius), 0.001, 100.0, 10, 2, "Radius to look for nearby particles within");
-		uiDefButS(block, NUM, B_REDR, "Nearby: ",
-			X2CLM1, yco-=BUTH, BUTW2, BUTH, &(pd->nearest), 2.0, 25.0, 10, 2, "The number of nearby particles to check for density (higher is more accurate, but slower)");
-		uiBlockEndAlign(block);
+		
+		yco -= YSPACE;
+		
+		uiDefBut(block, LABEL, B_NOP, "Falloff:",
+			X2CLM1, yco-=BUTH, BUTW2, BUTH, 0, 0, 0, 0, 0, "");	
+		uiDefButS(block, MENU, B_REDR, "Standard %x0|Smooth %x1|Sharp %x2",
+				X2CLM1, yco-=BUTH, BUTW2, BUTH, &pd->falloff_type, 0.0, 0.0, 0, 0, "Falloff type");
 		
 		yco = PANEL_YMAX;