Better grid rasterization method for hair volumetric grids.

This is based on the paper
"Detail Preserving Continuum Simulation of Straight Hair"
(McAdams, Selle, Ward, 2009)

The main difference is that hair line segments are used rather than only
rasterizing velocity at the vertices. This gives a much better coverage
of the hair volume grid, otherwise gaps can be produced at smaller grid
cell sizes and the distribution is uneven along the hair curve.

The algorithm for rasterizing is a variation of Bresenham's algorithm
extended onto 3D grids.

Conflicts:
	source/blender/physics/intern/BPH_mass_spring.cpp

3 files changed, 373 insertions, 10 deletions

diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp
index 7a2d85f140d..16828d83f6e 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -247,6 +247,7 @@ static void cloth_setup_constraints(ClothModifierData *clmd, ColliderContacts *c
 //			BKE_sim_debug_data_add_dot(clmd->debug_data, collpair->pb, 1, 0, 0, "collision", hash_collpair(937, collpair));

 //			BKE_sim_debug_data_add_line(clmd->debug_data, collpair->pa, collpair->pb, 0.7, 0.7, 0.7, "collision", hash_collpair(938, collpair));

 			

+#if 0

 			{ /* DEBUG */

 				float nor[3];

 				mul_v3_v3fl(nor, collpair->normal, -collpair->distance);

@@ -254,6 +255,7 @@ static void cloth_setup_constraints(ClothModifierData *clmd, ColliderContacts *c
 //				BKE_sim_debug_data_add_vector(clmd->debug_data, collpair->pb, impulse, 1, 1, 0, "collision", hash_collpair(940, collpair));

 //				BKE_sim_debug_data_add_vector(clmd->debug_data, collpair->pb, collpair->normal, 1, 1, 0, "collision", hash_collpair(941, collpair));

 			}

+#endif

 		}

 	}

 }

@@ -419,6 +421,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
 		/* XXX assuming same restlen for ij and jk segments here, this can be done correctly for hair later */

 		BPH_mass_spring_force_spring_bending_angular(data, s->ij, s->kl, s->mn, s->target, kb, cb);

 		

+#if 0

 		{

 			float x_kl[3], x_mn[3], v[3], d[3];

 			

@@ -435,6 +438,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s,
 //			BKE_sim_debug_data_add_vector(clmd->debug_data, x, d, 1, 0.4, 0.4, "target", hash_vertex(7983, s->kl));

 		}

 #endif

+#endif

 	}

 }

 

@@ -584,8 +588,140 @@ static void cloth_calc_volume_force(ClothModifierData *clmd)
 }

 #endif

 

+/* returns active vertexes' motion state, or original location the vertex is disabled */

+BLI_INLINE bool cloth_get_grid_location(Implicit_Data *data, const float cell_scale[3], const float cell_offset[3],

+                                        ClothVertex *vert, float x[3], float v[3])

+{

+	bool is_motion_state;

+	if (vert->solver_index < 0) {

+		copy_v3_v3(x, vert->x);

+		copy_v3_v3(v, vert->v);

+		is_motion_state = false;

+	}

+	else {

+		BPH_mass_spring_get_position(data, vert->solver_index, x);

+		BPH_mass_spring_get_new_velocity(data, vert->solver_index, v);

+		is_motion_state = true;

+	}

+	

+	mul_v3_v3(x, cell_scale);

+	add_v3_v3(x, cell_offset);

+	

+	return is_motion_state;

+}

+

+/* returns next spring forming a continous hair sequence */

+BLI_INLINE LinkNode *hair_spring_next(LinkNode *spring_link)

+{

+	ClothSpring *spring = (ClothSpring *)spring_link->link;

+	LinkNode *next = spring_link->next;

+	if (next) {

+		ClothSpring *next_spring = (ClothSpring *)next->link;

+		if (next_spring->type == CLOTH_SPRING_TYPE_STRUCTURAL && next_spring->kl == spring->ij)

+			return next;

+	}

+	return NULL;

+}

+

+/* XXX this is nasty: cloth meshes do not explicitly store

+ * the order of hair segments!

+ * We have to rely on the spring build function for now,

+ * which adds structural springs in reverse order:

+ *   (3,4), (2,3), (1,2)

+ * This is currently the only way to figure out hair geometry inside this code ...

+ */

+static LinkNode *cloth_continuum_add_hair_segments(HairVertexGrid *grid, const float cell_scale[3], const float cell_offset[3], Cloth *cloth, LinkNode *spring_link)

+{

+	Implicit_Data *data = cloth->implicit;

+	LinkNode *next_spring_link = NULL; /* return value */

+	ClothSpring *spring1, *spring2, *spring3;

+	ClothVertex *verts = cloth->verts;

+	ClothVertex *vert3, *vert4;

+	float x1[3], v1[3], x2[3], v2[3], x3[3], v3[3], x4[3], v4[3];

+	float dir1[3], dir2[3], dir3[3];

+	

+	spring1 = NULL;

+	spring2 = NULL;

+	spring3 = (ClothSpring *)spring_link->link;

+	

+	zero_v3(x1); zero_v3(v1);

+	zero_v3(dir1);

+	zero_v3(x2); zero_v3(v2);

+	zero_v3(dir2);

+	

+	vert3 = &verts[spring3->kl];

+	cloth_get_grid_location(data, cell_scale, cell_offset, vert3, x3, v3);

+	vert4 = &verts[spring3->ij];

+	cloth_get_grid_location(data, cell_scale, cell_offset, vert4, x4, v4);

+	sub_v3_v3v3(dir3, x4, x3);

+	normalize_v3(dir3);

+	

+	while (spring_link) {

+		/* move on */

+		spring1 = spring2;

+		spring2 = spring3;

+		

+		vert3 = vert4;

+		

+		copy_v3_v3(x1, x2); copy_v3_v3(v1, v2);

+		copy_v3_v3(x2, x3); copy_v3_v3(v2, v3);

+		copy_v3_v3(x3, x4); copy_v3_v3(v3, v4);

+		

+		copy_v3_v3(dir1, dir2);

+		copy_v3_v3(dir2, dir3);

+		

+		/* read next segment */

+		next_spring_link = spring_link->next;

+		spring_link = hair_spring_next(spring_link);

+		

+		if (spring_link) {

+			spring3 = (ClothSpring *)spring_link->link;

+			vert4 = &verts[spring3->ij];

+			cloth_get_grid_location(data, cell_scale, cell_offset, vert4, x4, v4);

+			sub_v3_v3v3(dir3, x4, x3);

+			normalize_v3(dir3);

+		}

+		else {

+			spring3 = NULL;

+			vert4 = NULL;

+			zero_v3(x4); zero_v3(v4);

+			zero_v3(dir3);

+		}

+		

+		BPH_hair_volume_add_segment(grid, x1, v1, x2, v2, x3, v3, x4, v4,

+		                            spring1 ? dir1 : NULL,

+		                            dir2,

+		                            spring3 ? dir3 : NULL);

+	}

+	

+	/* last segment */

+	spring1 = spring2;

+	spring2 = spring3;

+	spring3 = NULL;

+	

+	vert3 = vert4;

+	vert4 = NULL;

+	

+	copy_v3_v3(x1, x2); copy_v3_v3(v1, v2);

+	copy_v3_v3(x2, x3); copy_v3_v3(v2, v3);

+	copy_v3_v3(x3, x4); copy_v3_v3(v3, v4);

+	zero_v3(x4);        zero_v3(v4);

+	

+	copy_v3_v3(dir1, dir2);

+	copy_v3_v3(dir2, dir3);

+	zero_v3(dir3);

+	

+	BPH_hair_volume_add_segment(grid, x1, v1, x2, v2, x3, v3, x4, v4,

+	                            spring1 ? dir1 : NULL,

+	                            dir2,

+	                            NULL);

+	

+	return next_spring_link;

+}

+

 static void cloth_continuum_fill_grid(HairVertexGrid *grid, Cloth *cloth)

 {

+#if 0

 	Implicit_Data *data = cloth->implicit;

 	int numverts = cloth->numverts;

 	ClothVertex *vert;

@@ -594,10 +730,32 @@ static void cloth_continuum_fill_grid(HairVertexGrid *grid, Cloth *cloth)
 	for (i = 0, vert = cloth->verts; i < numverts; i++, vert++) {

 		float x[3], v[3];

 		

-		BPH_mass_spring_get_position(data, i, x);

-		BPH_mass_spring_get_new_velocity(data, i, v);

+		cloth_get_vertex_motion_state(data, vert, x, v);

 		BPH_hair_volume_add_vertex(grid, x, v);

 	}

+#else

+	LinkNode *link;

+	float cellsize[3], gmin[3], cell_scale[3], cell_offset[3];

+	

+	/* scale and offset for transforming vertex locations into grid space

+	 * (cell size is 0..1, gmin becomes origin)

+	 */

+	BPH_hair_volume_grid_geometry(grid, cellsize, NULL, gmin, NULL);

+	cell_scale[0] = cellsize[0] > 0.0f ? 1.0f / cellsize[0] : 0.0f;

+	cell_scale[1] = cellsize[1] > 0.0f ? 1.0f / cellsize[1] : 0.0f;

+	cell_scale[2] = cellsize[2] > 0.0f ? 1.0f / cellsize[2] : 0.0f;

+	mul_v3_v3v3(cell_offset, gmin, cell_scale);

+	negate_v3(cell_offset);

+	

+	link = cloth->springs;

+	while (link) {

+		ClothSpring *spring = (ClothSpring *)link->link;

+		if (spring->type == CLOTH_SPRING_TYPE_STRUCTURAL)

+			link = cloth_continuum_add_hair_segments(grid, cell_scale, cell_offset, cloth, link);

+		else

+			link = link->next;

+	}

+#endif

 	BPH_hair_volume_normalize_vertex_grid(grid);

 }

 

@@ -614,6 +772,7 @@ static void cloth_continuum_step(ClothModifierData *clmd)
 	float pressfac = parms->pressure;

 	float minpress = parms->pressure_threshold;

 	float gmin[3], gmax[3];

+	float cellsize[3];

 	int i;

 	

 	/* clear grid info */

@@ -626,6 +785,10 @@ static void cloth_continuum_step(ClothModifierData *clmd)
 	/* gather velocities & density */

 	if (smoothfac > 0.0f || pressfac > 0.0f) {

 		HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax);

+		BPH_hair_volume_set_debug_data(vertex_grid, clmd->debug_data);

+		

+		BPH_hair_volume_grid_geometry(vertex_grid, cellsize, NULL, NULL, NULL);

+		

 		cloth_continuum_fill_grid(vertex_grid, cloth);

 		

 #if 0

@@ -667,16 +830,20 @@ static void cloth_continuum_step(ClothModifierData *clmd)
 			a[(axis+1) % 3] = clmd->hair_grid_max[(axis+1) % 3] - clmd->hair_grid_min[(axis+1) % 3];

 			b[(axis+2) % 3] = clmd->hair_grid_max[(axis+2) % 3] - clmd->hair_grid_min[(axis+2) % 3];

 			

+			BKE_sim_debug_data_clear_category(clmd->debug_data, "grid velocity");

 			for (j = 0; j < size; ++j) {

 				for (i = 0; i < size; ++i) {

-					float x[3], v[3], nv[3];

+					float x[3], v[3], gvel[3], gdensity;

+					

 					madd_v3_v3v3fl(x, offset, a, (float)i / (float)(size-1));

 					madd_v3_v3fl(x, b, (float)j / (float)(size-1));

 					zero_v3(v);

 					

-					BPH_hair_volume_grid_velocity(vertex_grid, x, v, 0.0f, nv);

+					BPH_hair_volume_grid_interpolate(vertex_grid, x, &gdensity, gvel, NULL, NULL);

 					

-					BKE_sim_debug_data_add_vector(clmd->debug_data, x, nv, 0.4, 0, 1, "grid velocity", hash_int_2d(hash_int_2d(i, j), 3112));

+//					BKE_sim_debug_data_add_circle(clmd->debug_data, x, gdensity, 0.7, 0.3, 1, "grid density", hash_int_2d(hash_int_2d(i, j), 3111));

+					if (!is_zero_v3(gvel))

+						BKE_sim_debug_data_add_vector(clmd->debug_data, x, gvel, 0.4, 0, 1, "grid velocity", hash_int_2d(hash_int_2d(i, j), 3112));

 				}

 			}

 		}

diff --git a/source/blender/physics/intern/hair_volume.c b/source/blender/physics/intern/hair_volume.c
index f083cdfd479..9e9051a3a92 100644
--- a/source/blender/physics/intern/hair_volume.c
+++ b/source/blender/physics/intern/hair_volume.c
@@ -36,6 +36,8 @@
 
 #include "DNA_texture_types.h"
 
+#include "BKE_effect.h"
+
 #include "implicit.h"
 
 /* ================ Volumetric Hair Interaction ================
@@ -87,6 +89,8 @@ typedef struct HairVertexGrid {
 	int res;
 	float gmin[3], gmax[3];
 	float scale[3];
+	
+	struct SimDebugData *debug_data;
 } HairVertexGrid;
 
 typedef struct HairColliderGrid {
@@ -231,6 +235,12 @@ void BPH_hair_volume_vertex_grid_forces(HairVertexGrid *grid, const float x[3],
 	mul_m3_fl(dfdv, smoothfac);
 }
 
+void BPH_hair_volume_grid_interpolate(HairVertexGrid *grid, const float x[3],
+                                      float *density, float velocity[3], float density_gradient[3], float velocity_gradient[3][3])
+{
+	hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->scale, x, density, velocity, density_gradient, velocity_gradient);
+}
+
 void BPH_hair_volume_grid_velocity(HairVertexGrid *grid, const float x[3], const float v[3],
                                    float fluid_factor,
                                    float r_v[3])
@@ -318,6 +328,168 @@ void BPH_hair_volume_add_vertex(HairVertexGrid *grid, const float x[3], const fl
 	}
 }
 
+BLI_INLINE void hair_volume_eval_grid_vertex(HairGridVert *vert, const float loc[3], float radius, float dist_scale,
+                                             const float x2[3], const float v2[3], const float x3[3], const float v3[3])
+{
+	float closest[3], lambda, dist, weight;
+	
+	lambda = closest_to_line_v3(closest, loc, x2, x3);
+	dist = len_v3v3(closest, loc);
+	
+	weight = (radius - dist) * dist_scale;
+	
+	if (weight > 0.0f) {
+		float vel[3];
+		
+		interp_v3_v3v3(vel, v2, v3, lambda);
+		madd_v3_v3fl(vert->velocity, vel, weight);
+		vert->density += weight;
+	}
+}
+
+BLI_INLINE int major_axis_v3(const float v[3])
+{
+	return v[0] > v[1] ? (v[0] > v[2] ? 0 : 2) : (v[1] > v[2] ? 1 : 2);
+}
+
+BLI_INLINE void grid_to_world(HairVertexGrid *grid, float vecw[3], const float vec[3])
+{
+	copy_v3_v3(vecw, vec);
+	mul_v3_v3(vecw, grid->scale);
+	add_v3_v3(vecw, grid->gmin);
+}
+
+/* Uses a variation of Bresenham's algorithm for rasterizing a 3D grid with a line segment.
+ * 
+ * The radius of influence around a segment is assumed to be at most 2*cellsize,
+ * i.e. only cells containing the segment and their direct neighbors are examined.
+ * 
+ * 
+ */
+void BPH_hair_volume_add_segment(HairVertexGrid *grid,
+                                 const float UNUSED(x1[3]), const float UNUSED(v1[3]), const float x2[3], const float v2[3],
+                                 const float x3[3], const float v3[3], const float UNUSED(x4[3]), const float UNUSED(v4[3]),
+                                 const float UNUSED(dir1[3]), const float dir2[3], const float UNUSED(dir3[3]))
+{
+	SimDebugData *debug_data = grid->debug_data;
+	
+	const int res[3] = { grid->res, grid->res, grid->res };
+	
+	/* find the primary direction from the major axis of the direction vector */
+	const int axis0 = major_axis_v3(dir2);
+	const int axis1 = (axis0 + 1) % 3;
+	const int axis2 = (axis0 + 2) % 3;
+	
+	/* range along primary direction */
+	const float h2 = x2[axis0], h3 = x3[axis0];
+	const float hmin = min_ff(h2, h3);
+	const float hmax = max_ff(h2, h3);
+	const int imin = max_ii((int)hmin, 0);
+	const int imax = min_ii((int)hmax + 1, res[axis0]);
+	
+	const float inc[2] = { dir2[axis1], dir2[axis2] };          /* increment of secondary directions per step in the primary direction */
+	const int grid_start1 = (int)x2[axis1];                     /* offset of cells on minor axes */
+	const int grid_start2 = (int)x2[axis2];                     /* offset of cells on minor axes */
+	
+	const float cellsize[3] = { grid->scale[axis0], grid->scale[axis1], grid->scale[axis2] };
+	float shift[2] = { x2[axis1] - floorf(x2[axis1]),           /* fraction of a full cell shift [0.0, 1.0) */
+	                   x2[axis2] - floorf(x2[axis2]) };
+	
+	/* vertex buffer offset factors along cardinal axes */
+	const int strides[3] = { 1, res[0], res[0] * res[1] };
+	/* change in offset when incrementing one of the axes */
+	const int stride0 = strides[axis0];
+	const int stride1 = strides[axis1];
+	const int stride2 = strides[axis2];
+	
+	const float radius = 1.5f;
+	/* XXX cell size should be fixed and uniform! */
+	const float dist_scale = 1.0f / cellsize[0];
+	
+	HairGridVert *vert0;
+	float loc0[3];
+	int j0, k0;
+	int i;
+	
+	(void)debug_data;
+	
+	j0 = grid_start1 - 1;
+	k0 = grid_start2 - 1;
+	vert0 = grid->verts + stride0 * imin + stride1 * j0 + stride2 * k0;
+	loc0[axis0] = (float)imin;
+	loc0[axis1] = (float)j0;
+	loc0[axis2] = (float)k0;
+	
+	/* loop over all planes crossed along the primary direction */
+	for (i = imin; i < imax; ++i, vert0 += stride0, loc0[axis0] += cellsize[0]) {
+		const int jmin = max_ii(j0, 0);
+		const int jmax = min_ii(j0 + 5, res[axis1]);
+		const int kmin = max_ii(k0, 0);
+		const int kmax = min_ii(k0 + 5, res[axis2]);
+		
+		/* XXX problem: this can be offset beyond range of this plane when jmin/kmin gets clamped,
+		 * for now simply calculate in outer loop with multiplication once
+		 */
+//		HairGridVert *vert1 = vert0;
+//		float loc1[3] = { loc0[0], loc0[1], loc0[2] };
+		HairGridVert *vert1 = grid->verts + stride0 * i + stride1 * jmin + stride2 * kmin;
+		float loc1[3];
+		int j, k;
+		
+		/* note: loc is in grid cell units,
+		 * distances are be scaled by cell size for weighting
+		 */
+		loc1[axis0] = (float)i;
+		loc1[axis1] = (float)jmin;
+		loc1[axis2] = (float)kmin;
+		
+		/* 2x2 cells can be hit directly by the segment between two planes,
+		 * margin is 1 cell, i.e. 4x4 cells are influenced at most,
+		 * -> evaluate 5x5 grid vertices on cell borders
+		 */
+		
+		for (j = jmin; j < jmax; ++j, vert1 += stride1, loc1[axis1] += 1.0f) {
+			HairGridVert *vert2 = vert1;
+			float loc2[3] = { loc1[0], loc1[1], loc1[2] };
+			
+			for (k = kmin; k < kmax; ++k, vert2 += stride2, loc2[axis2] += 1.0f) {
+				hair_volume_eval_grid_vertex(vert2, loc2, radius, dist_scale, x2, v2, x3, v3);
+			}
+		}
+		
+		/* increment */
+		add_v2_v2(shift, inc);
+		if (shift[0] > 1.0f) {
+			shift[0] -= 1.0f;
+			
+			j0 += 1;
+			vert0 += stride1;
+			loc0[axis1] += 1.0f;
+		}
+		else if (shift[0] < -1.0f) {
+			shift[0] += 1.0f;
+			
+			j0 -= 1;
+			vert0 -= stride1;
+			loc0[axis1] -= 1.0f;
+		}
+		if (shift[1] > 1.0f) {
+			shift[1] -= 1.0f;
+			
+			k0 += 1;
+			vert0 += stride2;
+			loc0[axis2] += 1.0f;
+		}
+		else if (shift[1] < -1.0f) {
+			shift[1] += 1.0f;
+			
+			k0 -= 1;
+			vert0 -= stride2;
+			loc0[axis2] -= 1.0f;
+		}
+	}
+}
+
 void BPH_hair_volume_normalize_vertex_grid(HairVertexGrid *grid)
 {
 	int i, size = hair_grid_size(grid->res);
@@ -402,15 +574,25 @@ void BPH_hair_volume_vertex_grid_filter_box(HairVertexGrid *grid, int kernel_siz
 
 HairVertexGrid *BPH_hair_volume_create_vertex_grid(int res, const float gmin[3], const float gmax[3])
 {
-	int size = hair_grid_size(res);
+	float cellsize[3], gmin_margin[3], gmax_margin[3];
+	int size;
 	HairVertexGrid *grid;
-	int	            i = 0;
+	int i;
+	
+	/* original cell size, before adding margin */
+	hair_grid_get_scale(res, gmin, gmax, cellsize);
+	
+	/* add margin of 1 cell */
+	res += 2;
+	size = hair_grid_size(res);
+	sub_v3_v3v3(gmin_margin, gmin, cellsize);
+	add_v3_v3v3(gmax_margin, gmax, cellsize);
 	
 	grid = MEM_callocN(sizeof(HairVertexGrid), "hair vertex grid");
 	grid->res = res;
-	copy_v3_v3(grid->gmin, gmin);
-	copy_v3_v3(grid->gmax, gmax);
-	hair_grid_get_scale(res, gmin, gmax, grid->scale);
+	copy_v3_v3(grid->gmin, gmin_margin);
+	copy_v3_v3(grid->gmax, gmax_margin);
+	copy_v3_v3(grid->scale, cellsize);
 	grid->verts = MEM_mallocN(sizeof(HairGridVert) * size, "hair voxel data");
 
 	/* initialize grid */
@@ -431,6 +613,11 @@ void BPH_hair_volume_free_vertex_grid(HairVertexGrid *grid)
 	}
 }
 
+void BPH_hair_volume_set_debug_data(HairVertexGrid *grid, SimDebugData *debug_data)
+{
+	grid->debug_data = debug_data;
+}
+
 void BPH_hair_volume_grid_geometry(HairVertexGrid *grid, float cellsize[3], int res[3], float gmin[3], float gmax[3])
 {
 	if (cellsize) copy_v3_v3(cellsize, grid->scale);
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index b7eeea146bc..c62d26b508f 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -173,14 +173,23 @@ struct VoxelData;
 
 struct HairVertexGrid *BPH_hair_volume_create_vertex_grid(int res, const float gmin[3], const float gmax[3]);
 void BPH_hair_volume_free_vertex_grid(struct HairVertexGrid *grid);
+void BPH_hair_volume_set_debug_data(struct HairVertexGrid *grid, struct SimDebugData *debug_data);
 void BPH_hair_volume_grid_geometry(struct HairVertexGrid *grid, float cellsize[3], int res[3], float gmin[3], float gmax[3]);
 
 void BPH_hair_volume_add_vertex(struct HairVertexGrid *grid, const float x[3], const float v[3]);
+void BPH_hair_volume_add_segment(struct HairVertexGrid *grid,
+                                 const float x1[3], const float v1[3], const float x2[3], const float v2[3],
+                                 const float x3[3], const float v3[3], const float x4[3], const float v4[3],
+                                 const float dir1[3], const float dir2[3], const float dir3[3]);
+
 void BPH_hair_volume_normalize_vertex_grid(struct HairVertexGrid *grid);
 #if 0 /* XXX weighting is incorrect, disabled for now */
 void BPH_hair_volume_vertex_grid_filter_box(struct HairVertexGrid *grid, int kernel_size);
 #endif
 
+void BPH_hair_volume_grid_interpolate(struct HairVertexGrid *grid, const float x[3],
+                                      float *density, float velocity[3], float density_gradient[3], float velocity_gradient[3][3]);
+
 /* Effect of fluid simulation grid on velocities.
  * fluid_factor controls blending between PIC (Particle-in-Cell)
  *     and FLIP (Fluid-Implicit-Particle) methods (0 = only PIC, 1 = only FLIP)