* Point Density texture
authorMatt Ebb <matt@mke3.net>
Wed, 1 Oct 2008 03:35:53 +0000 (03:35 +0000)
committerMatt Ebb <matt@mke3.net>
Wed, 1 Oct 2008 03:35:53 +0000 (03:35 +0000)
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

source/blender/blenkernel/intern/texture.c
source/blender/blenlib/BLI_kdopbvh.h
source/blender/blenlib/intern/BLI_kdopbvh.c
source/blender/makesdna/DNA_texture_types.h
source/blender/render/intern/source/pointdensity.c
source/blender/src/buttons_shading.c

index 45fd11faa3227485462e3772c9093e70d2601575..ec34fcf4103f6e24bc9b59c84ab2994cf98e5e30 100644 (file)
@@ -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)
index e3591a84e98c5fbef717da7793825b63a03e5998..8912ac5bd13927f24b9c00989963c7e93678fee1 100644 (file)
@@ -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
 
index 30472beb3e638ff04eb6222f865ca7eeda9a0d31..f82b6b7f162b04761d6e43881cfa411c3666ae40 100644 (file)
@@ -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;
+}
index 52055f24f82018329a834408fbe2b7119d1407ec..d23672f9284678221d6673793c33d8540c581943 100644 (file)
@@ -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
 
index 17bc6bb58cf3247dde658a5a8bdf2b138dba603f..4422b9fbbdd95ef4fc2b93bd46b504a5254486d6 100644 (file)
@@ -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;
 
index 85154682819ca49a77cb84f2dd23f3332641a6c3..050dae6026ed515bea540e612b5d0c178f824d77 100644 (file)
@@ -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;