cd869046cc1ea9232d9af32ede52288701d1932b
[blender.git] / source / blender / physics / intern / hair_volume.cpp
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) Blender Foundation
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Janne Karhu, Lukas Toenne
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/physics/intern/hair_volume.cpp
29  *  \ingroup bph
30  */
31
32 #include "MEM_guardedalloc.h"
33
34 extern "C" {
35 #include "BLI_math.h"
36 #include "BLI_utildefines.h"
37
38 #include "DNA_texture_types.h"
39
40 #include "BKE_effect.h"
41 }
42
43 #include "implicit.h"
44 #include "eigen_utils.h"
45
46 /* ================ Volumetric Hair Interaction ================
47  * adapted from
48  *
49  * Volumetric Methods for Simulation and Rendering of Hair
50  *     (Petrovic, Henne, Anderson, Pixar Technical Memo #06-08, Pixar Animation Studios)
51  *
52  * as well as
53  *
54  * "Detail Preserving Continuum Simulation of Straight Hair"
55  *     (McAdams, Selle 2009)
56  */
57
58 /* Note about array indexing:
59  * Generally the arrays here are one-dimensional.
60  * The relation between 3D indices and the array offset is
61  *   offset = x + res_x * y + res_x * res_y * z
62  */
63
64 static float I[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
65
66 BLI_INLINE int floor_int(float value)
67 {
68         return value > 0.0f ? (int)value : ((int)value) - 1;
69 }
70
71 BLI_INLINE float floor_mod(float value)
72 {
73         return value - floorf(value);
74 }
75
76 BLI_INLINE int hair_grid_size(const int res[3])
77 {
78         return res[0] * res[1] * res[2];
79 }
80
81 typedef struct HairGridVert {
82         int samples;
83         float velocity[3];
84         float density;
85
86         float velocity_smooth[3];
87 } HairGridVert;
88
89 typedef struct HairGrid {
90         HairGridVert *verts;
91         int res[3];
92         float gmin[3], gmax[3];
93         float cellsize, inv_cellsize;
94 } HairGrid;
95
96 #define HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, axis) ( min_ii( max_ii( (int)((vec[axis] - gmin[axis]) * scale), 0), res[axis]-2 ) )
97
98 BLI_INLINE int hair_grid_offset(const float vec[3], const int res[3], const float gmin[3], float scale)
99 {
100         int i, j, k;
101         i = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 0);
102         j = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 1);
103         k = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 2);
104         return i + (j + k*res[1])*res[0];
105 }
106
107 BLI_INLINE int hair_grid_interp_weights(const int res[3], const float gmin[3], float scale, const float vec[3], float uvw[3])
108 {
109         int i, j, k, offset;
110
111         i = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 0);
112         j = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 1);
113         k = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 2);
114         offset = i + (j + k*res[1])*res[0];
115
116         uvw[0] = (vec[0] - gmin[0]) * scale - (float)i;
117         uvw[1] = (vec[1] - gmin[1]) * scale - (float)j;
118         uvw[2] = (vec[2] - gmin[2]) * scale - (float)k;
119
120 //      BLI_assert(0.0f <= uvw[0] && uvw[0] <= 1.0001f);
121 //      BLI_assert(0.0f <= uvw[1] && uvw[1] <= 1.0001f);
122 //      BLI_assert(0.0f <= uvw[2] && uvw[2] <= 1.0001f);
123
124         return offset;
125 }
126
127 BLI_INLINE void hair_grid_interpolate(const HairGridVert *grid, const int res[3], const float gmin[3], float scale, const float vec[3],
128                                       float *density, float velocity[3], float vel_smooth[3], float density_gradient[3], float velocity_gradient[3][3])
129 {
130         HairGridVert data[8];
131         float uvw[3], muvw[3];
132         int res2 = res[1] * res[0];
133         int offset;
134
135         offset = hair_grid_interp_weights(res, gmin, scale, vec, uvw);
136         muvw[0] = 1.0f - uvw[0];
137         muvw[1] = 1.0f - uvw[1];
138         muvw[2] = 1.0f - uvw[2];
139
140         data[0] = grid[offset              ];
141         data[1] = grid[offset            +1];
142         data[2] = grid[offset     +res[0]  ];
143         data[3] = grid[offset     +res[0]+1];
144         data[4] = grid[offset+res2         ];
145         data[5] = grid[offset+res2       +1];
146         data[6] = grid[offset+res2+res[0]  ];
147         data[7] = grid[offset+res2+res[0]+1];
148
149         if (density) {
150                 *density = muvw[2]*( muvw[1]*( muvw[0]*data[0].density + uvw[0]*data[1].density )   +
151                                       uvw[1]*( muvw[0]*data[2].density + uvw[0]*data[3].density ) ) +
152                             uvw[2]*( muvw[1]*( muvw[0]*data[4].density + uvw[0]*data[5].density )   +
153                                       uvw[1]*( muvw[0]*data[6].density + uvw[0]*data[7].density ) );
154         }
155
156         if (velocity) {
157                 int k;
158                 for (k = 0; k < 3; ++k) {
159                         velocity[k] = muvw[2]*( muvw[1]*( muvw[0]*data[0].velocity[k] + uvw[0]*data[1].velocity[k] )   +
160                                                  uvw[1]*( muvw[0]*data[2].velocity[k] + uvw[0]*data[3].velocity[k] ) ) +
161                                        uvw[2]*( muvw[1]*( muvw[0]*data[4].velocity[k] + uvw[0]*data[5].velocity[k] )   +
162                                                  uvw[1]*( muvw[0]*data[6].velocity[k] + uvw[0]*data[7].velocity[k] ) );
163                 }
164         }
165
166         if (vel_smooth) {
167                 int k;
168                 for (k = 0; k < 3; ++k) {
169                         vel_smooth[k] = muvw[2]*( muvw[1]*( muvw[0]*data[0].velocity_smooth[k] + uvw[0]*data[1].velocity_smooth[k] )   +
170                                                    uvw[1]*( muvw[0]*data[2].velocity_smooth[k] + uvw[0]*data[3].velocity_smooth[k] ) ) +
171                                          uvw[2]*( muvw[1]*( muvw[0]*data[4].velocity_smooth[k] + uvw[0]*data[5].velocity_smooth[k] )   +
172                                                    uvw[1]*( muvw[0]*data[6].velocity_smooth[k] + uvw[0]*data[7].velocity_smooth[k] ) );
173                 }
174         }
175
176         if (density_gradient) {
177                 density_gradient[0] = muvw[1] * muvw[2] * ( data[0].density - data[1].density ) +
178                                        uvw[1] * muvw[2] * ( data[2].density - data[3].density ) +
179                                       muvw[1] *  uvw[2] * ( data[4].density - data[5].density ) +
180                                        uvw[1] *  uvw[2] * ( data[6].density - data[7].density );
181
182                 density_gradient[1] = muvw[2] * muvw[0] * ( data[0].density - data[2].density ) +
183                                        uvw[2] * muvw[0] * ( data[4].density - data[6].density ) +
184                                       muvw[2] *  uvw[0] * ( data[1].density - data[3].density ) +
185                                        uvw[2] *  uvw[0] * ( data[5].density - data[7].density );
186
187                 density_gradient[2] = muvw[2] * muvw[0] * ( data[0].density - data[4].density ) +
188                                        uvw[2] * muvw[0] * ( data[1].density - data[5].density ) +
189                                       muvw[2] *  uvw[0] * ( data[2].density - data[6].density ) +
190                                        uvw[2] *  uvw[0] * ( data[3].density - data[7].density );
191         }
192
193         if (velocity_gradient) {
194                 /* XXX TODO */
195                 zero_m3(velocity_gradient);
196         }
197 }
198
199 void BPH_hair_volume_vertex_grid_forces(HairGrid *grid, const float x[3], const float v[3],
200                                         float smoothfac, float pressurefac, float minpressure,
201                                         float f[3], float dfdx[3][3], float dfdv[3][3])
202 {
203         float gdensity, gvelocity[3], ggrad[3], gvelgrad[3][3], gradlen;
204
205         hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, &gdensity, gvelocity, NULL, ggrad, gvelgrad);
206
207         zero_v3(f);
208         sub_v3_v3(gvelocity, v);
209         mul_v3_v3fl(f, gvelocity, smoothfac);
210
211         gradlen = normalize_v3(ggrad) - minpressure;
212         if (gradlen > 0.0f) {
213                 mul_v3_fl(ggrad, gradlen);
214                 madd_v3_v3fl(f, ggrad, pressurefac);
215         }
216
217         zero_m3(dfdx);
218
219         sub_m3_m3m3(dfdv, gvelgrad, I);
220         mul_m3_fl(dfdv, smoothfac);
221 }
222
223 void BPH_hair_volume_grid_interpolate(HairGrid *grid, const float x[3],
224                                       float *density, float velocity[3], float velocity_smooth[3], float density_gradient[3], float velocity_gradient[3][3])
225 {
226         hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, density, velocity, velocity_smooth, density_gradient, velocity_gradient);
227 }
228
229 void BPH_hair_volume_grid_velocity(HairGrid *grid, const float x[3], const float v[3],
230                                    float fluid_factor,
231                                    float r_v[3])
232 {
233         float gdensity, gvelocity[3], gvel_smooth[3], ggrad[3], gvelgrad[3][3];
234         float v_pic[3], v_flip[3];
235
236         hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->inv_cellsize, x, &gdensity, gvelocity, gvel_smooth, ggrad, gvelgrad);
237
238         /* velocity according to PIC method (Particle-in-Cell) */
239         copy_v3_v3(v_pic, gvel_smooth);
240
241         /* velocity according to FLIP method (Fluid-Implicit-Particle) */
242         sub_v3_v3v3(v_flip, gvel_smooth, gvelocity);
243         add_v3_v3(v_flip, v);
244
245         interp_v3_v3v3(r_v, v_pic, v_flip, fluid_factor);
246 }
247
248 void BPH_hair_volume_grid_clear(HairGrid *grid)
249 {
250         const int size = hair_grid_size(grid->res);
251         int i;
252         for (i = 0; i < size; ++i) {
253                 zero_v3(grid->verts[i].velocity);
254                 zero_v3(grid->verts[i].velocity_smooth);
255                 grid->verts[i].density = 0.0f;
256                 grid->verts[i].samples = 0;
257         }
258 }
259
260 BLI_INLINE bool hair_grid_point_valid(const float vec[3], float gmin[3], float gmax[3])
261 {
262         return !(vec[0] < gmin[0] || vec[1] < gmin[1] || vec[2] < gmin[2] ||
263                  vec[0] > gmax[0] || vec[1] > gmax[1] || vec[2] > gmax[2]);
264 }
265
266 BLI_INLINE float dist_tent_v3f3(const float a[3], float x, float y, float z)
267 {
268         float w = (1.0f - fabsf(a[0] - x)) * (1.0f - fabsf(a[1] - y)) * (1.0f - fabsf(a[2] - z));
269         return w;
270 }
271
272 BLI_INLINE float weights_sum(const float weights[8])
273 {
274         float totweight = 0.0f;
275         int i;
276         for (i = 0; i < 8; ++i)
277                 totweight += weights[i];
278         return totweight;
279 }
280
281 /* returns the grid array offset as well to avoid redundant calculation */
282 BLI_INLINE int hair_grid_weights(const int res[3], const float gmin[3], float scale, const float vec[3], float weights[8])
283 {
284         int i, j, k, offset;
285         float uvw[3];
286
287         i = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 0);
288         j = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 1);
289         k = HAIR_GRID_INDEX_AXIS(vec, res, gmin, scale, 2);
290         offset = i + (j + k*res[1])*res[0];
291
292         uvw[0] = (vec[0] - gmin[0]) * scale;
293         uvw[1] = (vec[1] - gmin[1]) * scale;
294         uvw[2] = (vec[2] - gmin[2]) * scale;
295
296         weights[0] = dist_tent_v3f3(uvw, (float)i    , (float)j    , (float)k    );
297         weights[1] = dist_tent_v3f3(uvw, (float)(i+1), (float)j    , (float)k    );
298         weights[2] = dist_tent_v3f3(uvw, (float)i    , (float)(j+1), (float)k    );
299         weights[3] = dist_tent_v3f3(uvw, (float)(i+1), (float)(j+1), (float)k    );
300         weights[4] = dist_tent_v3f3(uvw, (float)i    , (float)j    , (float)(k+1));
301         weights[5] = dist_tent_v3f3(uvw, (float)(i+1), (float)j    , (float)(k+1));
302         weights[6] = dist_tent_v3f3(uvw, (float)i    , (float)(j+1), (float)(k+1));
303         weights[7] = dist_tent_v3f3(uvw, (float)(i+1), (float)(j+1), (float)(k+1));
304
305 //      BLI_assert(fabsf(weights_sum(weights) - 1.0f) < 0.0001f);
306
307         return offset;
308 }
309
310 BLI_INLINE void grid_to_world(HairGrid *grid, float vecw[3], const float vec[3])
311 {
312         copy_v3_v3(vecw, vec);
313         mul_v3_fl(vecw, grid->cellsize);
314         add_v3_v3(vecw, grid->gmin);
315 }
316
317 void BPH_hair_volume_add_vertex(HairGrid *grid, const float x[3], const float v[3])
318 {
319         const int res[3] = { grid->res[0], grid->res[1], grid->res[2] };
320         float weights[8];
321         int di, dj, dk;
322         int offset;
323
324         if (!hair_grid_point_valid(x, grid->gmin, grid->gmax))
325                 return;
326
327         offset = hair_grid_weights(res, grid->gmin, grid->inv_cellsize, x, weights);
328
329         for (di = 0; di < 2; ++di) {
330                 for (dj = 0; dj < 2; ++dj) {
331                         for (dk = 0; dk < 2; ++dk) {
332                                 int voffset = offset + di + (dj + dk*res[1])*res[0];
333                                 int iw = di + dj*2 + dk*4;
334
335                                 grid->verts[voffset].density += weights[iw];
336                                 madd_v3_v3fl(grid->verts[voffset].velocity, v, weights[iw]);
337                         }
338                 }
339         }
340 }
341
342 #if 0
343 BLI_INLINE void hair_volume_eval_grid_vertex(HairGridVert *vert, const float loc[3], float radius, float dist_scale,
344                                              const float x2[3], const float v2[3], const float x3[3], const float v3[3])
345 {
346         float closest[3], lambda, dist, weight;
347
348         lambda = closest_to_line_v3(closest, loc, x2, x3);
349         dist = len_v3v3(closest, loc);
350
351         weight = (radius - dist) * dist_scale;
352
353         if (weight > 0.0f) {
354                 float vel[3];
355
356                 interp_v3_v3v3(vel, v2, v3, lambda);
357                 madd_v3_v3fl(vert->velocity, vel, weight);
358                 vert->density += weight;
359                 vert->samples += 1;
360         }
361 }
362
363 BLI_INLINE int major_axis_v3(const float v[3])
364 {
365         const float a = fabsf(v[0]);
366         const float b = fabsf(v[1]);
367         const float c = fabsf(v[2]);
368         return a > b ? (a > c ? 0 : 2) : (b > c ? 1 : 2);
369 }
370
371 BLI_INLINE void hair_volume_add_segment_2D(HairGrid *grid,
372                                            const float UNUSED(x1[3]), const float UNUSED(v1[3]), const float x2[3], const float v2[3],
373                                            const float x3[3], const float v3[3], const float UNUSED(x4[3]), const float UNUSED(v4[3]),
374                                            const float UNUSED(dir1[3]), const float dir2[3], const float UNUSED(dir3[3]),
375                                            int resj, int resk, int jmin, int jmax, int kmin, int kmax,
376                                            HairGridVert *vert, int stride_j, int stride_k, const float loc[3], int axis_j, int axis_k,
377                                            int debug_i)
378 {
379         const float radius = 1.5f;
380         const float dist_scale = grid->inv_cellsize;
381
382         int j, k;
383
384         /* boundary checks to be safe */
385         CLAMP_MIN(jmin, 0);
386         CLAMP_MAX(jmax, resj-1);
387         CLAMP_MIN(kmin, 0);
388         CLAMP_MAX(kmax, resk-1);
389
390         HairGridVert *vert_j = vert + jmin * stride_j;
391         float loc_j[3] = { loc[0], loc[1], loc[2] };
392         loc_j[axis_j] += (float)jmin;
393         for (j = jmin; j <= jmax; ++j, vert_j += stride_j, loc_j[axis_j] += 1.0f) {
394
395                 HairGridVert *vert_k = vert_j + kmin * stride_k;
396                 float loc_k[3] = { loc_j[0], loc_j[1], loc_j[2] };
397                 loc_k[axis_k] += (float)kmin;
398                 for (k = kmin; k <= kmax; ++k, vert_k += stride_k, loc_k[axis_k] += 1.0f) {
399
400                         hair_volume_eval_grid_vertex(vert_k, loc_k, radius, dist_scale, x2, v2, x3, v3);
401
402 #if 0
403                         {
404                                 float wloc[3], x2w[3], x3w[3];
405                                 grid_to_world(grid, wloc, loc_k);
406                                 grid_to_world(grid, x2w, x2);
407                                 grid_to_world(grid, x3w, x3);
408
409                                 if (vert_k->samples > 0)
410                                         BKE_sim_debug_data_add_circle(wloc, 0.01f, 1.0, 1.0, 0.3, "grid", 2525, debug_i, j, k);
411
412                                 if (grid->debug_value) {
413                                         BKE_sim_debug_data_add_dot(wloc, 1, 0, 0, "grid", 93, debug_i, j, k);
414                                         BKE_sim_debug_data_add_dot(x2w, 0.1, 0.1, 0.7, "grid", 649, debug_i, j, k);
415                                         BKE_sim_debug_data_add_line(wloc, x2w, 0.3, 0.8, 0.3, "grid", 253, debug_i, j, k);
416                                         BKE_sim_debug_data_add_line(wloc, x3w, 0.8, 0.3, 0.3, "grid", 254, debug_i, j, k);
417 //                                      BKE_sim_debug_data_add_circle(x2w, len_v3v3(wloc, x2w), 0.2, 0.7, 0.2, "grid", 255, i, j, k);
418                                 }
419                         }
420 #endif
421                 }
422         }
423 }
424
425 /* Uses a variation of Bresenham's algorithm for rasterizing a 3D grid with a line segment.
426  *
427  * The radius of influence around a segment is assumed to be at most 2*cellsize,
428  * i.e. only cells containing the segment and their direct neighbors are examined.
429  *
430  *
431  */
432 void BPH_hair_volume_add_segment(HairGrid *grid,
433                                  const float x1[3], const float v1[3], const float x2[3], const float v2[3],
434                                  const float x3[3], const float v3[3], const float x4[3], const float v4[3],
435                                  const float dir1[3], const float dir2[3], const float dir3[3])
436 {
437         const int res[3] = { grid->res[0], grid->res[1], grid->res[2] };
438
439         /* find the primary direction from the major axis of the direction vector */
440         const int axis0 = major_axis_v3(dir2);
441         const int axis1 = (axis0 + 1) % 3;
442         const int axis2 = (axis0 + 2) % 3;
443
444         /* vertex buffer offset factors along cardinal axes */
445         const int strides[3] = { 1, res[0], res[0] * res[1] };
446         const int stride0 = strides[axis0];
447         const int stride1 = strides[axis1];
448         const int stride2 = strides[axis2];
449
450         /* increment of secondary directions per step in the primary direction
451          * note: we always go in the positive direction along axis0, so the sign can be inverted
452          */
453         const float inc1 = dir2[axis1] / dir2[axis0];
454         const float inc2 = dir2[axis2] / dir2[axis0];
455
456         /* start/end points, so increment along axis0 is always positive */
457         const float *start = x2[axis0] < x3[axis0] ? x2 : x3;
458         const float *end   = x2[axis0] < x3[axis0] ? x3 : x2;
459         const float start0 = start[axis0], start1 = start[axis1], start2 = start[axis2];
460         const float end0 = end[axis0];
461
462         /* range along primary direction */
463         const int imin = max_ii(floor_int(start[axis0]) - 1, 0);
464         const int imax = min_ii(floor_int(end[axis0]) + 2, res[axis0]-1);
465
466         float h = 0.0f;
467         HairGridVert *vert0;
468         float loc0[3];
469         int j0, k0, j0_prev, k0_prev;
470         int i;
471
472         for (i = imin; i <= imax; ++i) {
473                 float shift1, shift2; /* fraction of a full cell shift [0.0, 1.0) */
474                 int jmin, jmax, kmin, kmax;
475
476                 h = CLAMPIS((float)i, start0, end0);
477
478                 shift1 = start1 + (h - start0) * inc1;
479                 shift2 = start2 + (h - start0) * inc2;
480
481                 j0_prev = j0;
482                 j0 = floor_int(shift1);
483
484                 k0_prev = k0;
485                 k0 = floor_int(shift2);
486
487                 if (i > imin) {
488                         jmin = min_ii(j0, j0_prev);
489                         jmax = max_ii(j0, j0_prev);
490                         kmin = min_ii(k0, k0_prev);
491                         kmax = max_ii(k0, k0_prev);
492                 }
493                 else {
494                         jmin = jmax = j0;
495                         kmin = kmax = k0;
496                 }
497
498                 vert0 = grid->verts + i * stride0;
499                 loc0[axis0] = (float)i;
500                 loc0[axis1] = 0.0f;
501                 loc0[axis2] = 0.0f;
502
503                 hair_volume_add_segment_2D(grid, x1, v1, x2, v2, x3, v3, x4, v4, dir1, dir2, dir3,
504                                            res[axis1], res[axis2], jmin-1, jmax+2, kmin-1, kmax+2,
505                                            vert0, stride1, stride2, loc0, axis1, axis2,
506                                            i);
507         }
508 }
509 #else
510 BLI_INLINE void hair_volume_eval_grid_vertex_sample(HairGridVert *vert, const float loc[3], float radius, float dist_scale,
511                                                     const float x[3], const float v[3])
512 {
513         float dist, weight;
514
515         dist = len_v3v3(x, loc);
516
517         weight = (radius - dist) * dist_scale;
518
519         if (weight > 0.0f) {
520                 madd_v3_v3fl(vert->velocity, v, weight);
521                 vert->density += weight;
522                 vert->samples += 1;
523         }
524 }
525
526 /* XXX simplified test implementation using a series of discrete sample along the segment,
527  * instead of finding the closest point for all affected grid vertices.
528  */
529 void BPH_hair_volume_add_segment(HairGrid *grid,
530                                  const float UNUSED(x1[3]), const float UNUSED(v1[3]), const float x2[3], const float v2[3],
531                                  const float x3[3], const float v3[3], const float UNUSED(x4[3]), const float UNUSED(v4[3]),
532                                  const float UNUSED(dir1[3]), const float UNUSED(dir2[3]), const float UNUSED(dir3[3]))
533 {
534         const float radius = 1.5f;
535         const float dist_scale = grid->inv_cellsize;
536
537         const int res[3] = { grid->res[0], grid->res[1], grid->res[2] };
538         const int stride[3] = { 1, res[0], res[0] * res[1] };
539         const int num_samples = 10;
540
541         int s;
542
543         for (s = 0; s < num_samples; ++s) {
544                 float x[3], v[3];
545                 int i, j, k;
546
547                 float f = (float)s / (float)(num_samples-1);
548                 interp_v3_v3v3(x, x2, x3, f);
549                 interp_v3_v3v3(v, v2, v3, f);
550
551                 int imin = max_ii(floor_int(x[0]) - 2, 0);
552                 int imax = min_ii(floor_int(x[0]) + 2, res[0]-1);
553                 int jmin = max_ii(floor_int(x[1]) - 2, 0);
554                 int jmax = min_ii(floor_int(x[1]) + 2, res[1]-1);
555                 int kmin = max_ii(floor_int(x[2]) - 2, 0);
556                 int kmax = min_ii(floor_int(x[2]) + 2, res[2]-1);
557
558                 for (k = kmin; k <= kmax; ++k) {
559                         for (j = jmin; j <= jmax; ++j) {
560                                 for (i = imin; i <= imax; ++i) {
561                                         float loc[3] = { (float)i, (float)j, (float)k };
562                                         HairGridVert *vert = grid->verts + i * stride[0] + j * stride[1] + k * stride[2];
563
564                                         hair_volume_eval_grid_vertex_sample(vert, loc, radius, dist_scale, x, v);
565                                 }
566                         }
567                 }
568         }
569 }
570 #endif
571
572 void BPH_hair_volume_normalize_vertex_grid(HairGrid *grid)
573 {
574         int i, size = hair_grid_size(grid->res);
575         /* divide velocity with density */
576         for (i = 0; i < size; i++) {
577                 float density = grid->verts[i].density;
578                 if (density > 0.0f)
579                         mul_v3_fl(grid->verts[i].velocity, 1.0f/density);
580         }
581 }
582
583 static const float density_threshold = 0.001f; /* cells with density below this are considered empty */
584
585 /* Contribution of target density pressure to the laplacian in the pressure poisson equation.
586  * This is based on the model found in
587  * "Two-way Coupled SPH and Particle Level Set Fluid Simulation" (Losasso et al., 2008)
588  */
589 BLI_INLINE float hair_volume_density_divergence(float density, float target_density, float strength)
590 {
591         if (density > density_threshold && density > target_density)
592                 return strength * logf(target_density / density);
593         else
594                 return 0.0f;
595 }
596
597 bool BPH_hair_volume_solve_divergence(HairGrid *grid, float /*dt*/, float target_density, float target_strength)
598 {
599         const float flowfac = grid->cellsize;
600         const float inv_flowfac = 1.0f / grid->cellsize;
601
602         /*const int num_cells = hair_grid_size(grid->res);*/
603         const int res[3] = { grid->res[0], grid->res[1], grid->res[2] };
604         const int resA[3] = { grid->res[0] + 2, grid->res[1] + 2, grid->res[2] + 2 };
605
606         const int stride0 = 1;
607         const int stride1 = grid->res[0];
608         const int stride2 = grid->res[1] * grid->res[0];
609         const int strideA0 = 1;
610         const int strideA1 = grid->res[0] + 2;
611         const int strideA2 = (grid->res[1] + 2) * (grid->res[0] + 2);
612
613         const int num_cells = res[0] * res[1] * res[2];
614         const int num_cellsA = (res[0] + 2) * (res[1] + 2) * (res[2] + 2);
615
616         HairGridVert *vert_start = grid->verts - (stride0 + stride1 + stride2);
617         HairGridVert *vert;
618         int i, j, k;
619
620 #define MARGIN_i0 (i < 1)
621 #define MARGIN_j0 (j < 1)
622 #define MARGIN_k0 (k < 1)
623 #define MARGIN_i1 (i >= resA[0]-1)
624 #define MARGIN_j1 (j >= resA[1]-1)
625 #define MARGIN_k1 (k >= resA[2]-1)
626
627 #define NEIGHBOR_MARGIN_i0 (i < 2)
628 #define NEIGHBOR_MARGIN_j0 (j < 2)
629 #define NEIGHBOR_MARGIN_k0 (k < 2)
630 #define NEIGHBOR_MARGIN_i1 (i >= resA[0]-2)
631 #define NEIGHBOR_MARGIN_j1 (j >= resA[1]-2)
632 #define NEIGHBOR_MARGIN_k1 (k >= resA[2]-2)
633
634         BLI_assert(num_cells >= 1);
635
636         /* Calculate divergence */
637         lVector B(num_cellsA);
638         for (k = 0; k < resA[2]; ++k) {
639                 for (j = 0; j < resA[1]; ++j) {
640                         for (i = 0; i < resA[0]; ++i) {
641                                 int u = i * strideA0 + j * strideA1 + k * strideA2;
642                                 bool is_margin = MARGIN_i0 || MARGIN_i1 || MARGIN_j0 || MARGIN_j1 || MARGIN_k0 || MARGIN_k1;
643
644                                 if (is_margin) {
645                                         B[u] = 0.0f;
646                                         continue;
647                                 }
648
649                                 vert = vert_start + i * stride0 + j * stride1 + k * stride2;
650
651                                 const float *v0 = vert->velocity;
652                                 float dx = 0.0f, dy = 0.0f, dz = 0.0f;
653                                 if (!NEIGHBOR_MARGIN_i0)
654                                         dx += v0[0] - (vert - stride0)->velocity[0];
655                                 if (!NEIGHBOR_MARGIN_i1)
656                                         dx += (vert + stride0)->velocity[0] - v0[0];
657                                 if (!NEIGHBOR_MARGIN_j0)
658                                         dy += v0[1] - (vert - stride1)->velocity[1];
659                                 if (!NEIGHBOR_MARGIN_j1)
660                                         dy += (vert + stride1)->velocity[1] - v0[1];
661                                 if (!NEIGHBOR_MARGIN_k0)
662                                         dz += v0[2] - (vert - stride2)->velocity[2];
663                                 if (!NEIGHBOR_MARGIN_k1)
664                                         dz += (vert + stride2)->velocity[2] - v0[2];
665
666                                 float divergence = -0.5f * flowfac * (dx + dy + dz);
667
668                                 /* adjustment term for target density */
669                                 float target = hair_volume_density_divergence(vert->density, target_density, target_strength);
670
671                                 /* B vector contains the finite difference approximation of the velocity divergence.
672                                  * Note: according to the discretized Navier-Stokes equation the rhs vector
673                                  * and resulting pressure gradient should be multiplied by the (inverse) density;
674                                  * however, this is already included in the weighting of hair velocities on the grid!
675                                  */
676                                 B[u] = divergence - target;
677
678 #if 0
679                                 {
680                                         float wloc[3], loc[3];
681                                         float col0[3] = {0.0, 0.0, 0.0};
682                                         float colp[3] = {0.0, 1.0, 1.0};
683                                         float coln[3] = {1.0, 0.0, 1.0};
684                                         float col[3];
685                                         float fac;
686
687                                         loc[0] = (float)(i - 1);
688                                         loc[1] = (float)(j - 1);
689                                         loc[2] = (float)(k - 1);
690                                         grid_to_world(grid, wloc, loc);
691
692                                         if (divergence > 0.0f) {
693                                                 fac = CLAMPIS(divergence * target_strength, 0.0, 1.0);
694                                                 interp_v3_v3v3(col, col0, colp, fac);
695                                         }
696                                         else {
697                                                 fac = CLAMPIS(-divergence * target_strength, 0.0, 1.0);
698                                                 interp_v3_v3v3(col, col0, coln, fac);
699                                         }
700                                         if (fac > 0.05f)
701                                                 BKE_sim_debug_data_add_circle(grid->debug_data, wloc, 0.01f, col[0], col[1], col[2], "grid", 5522, i, j, k);
702                                 }
703 #endif
704                         }
705                 }
706         }
707
708         /* Main Poisson equation system:
709          * This is derived from the discretezation of the Poisson equation
710          *   div(grad(p)) = div(v)
711          *
712          * The finite difference approximation yields the linear equation system described here:
713          * https://en.wikipedia.org/wiki/Discrete_Poisson_equation
714          */
715         lMatrix A(num_cellsA, num_cellsA);
716         /* Reserve space for the base equation system (without boundary conditions).
717          * Each column contains a factor 6 on the diagonal
718          * and up to 6 factors -1 on other places.
719          */
720         A.reserve(Eigen::VectorXi::Constant(num_cellsA, 7));
721
722         for (k = 0; k < resA[2]; ++k) {
723                 for (j = 0; j < resA[1]; ++j) {
724                         for (i = 0; i < resA[0]; ++i) {
725                                 int u = i * strideA0 + j * strideA1 + k * strideA2;
726                                 bool is_margin = MARGIN_i0 || MARGIN_i1 || MARGIN_j0 || MARGIN_j1 || MARGIN_k0 || MARGIN_k1;
727
728                                 vert = vert_start + i * stride0 + j * stride1 + k * stride2;
729                                 if (!is_margin && vert->density > density_threshold) {
730                                         int neighbors_lo = 0;
731                                         int neighbors_hi = 0;
732                                         int non_solid_neighbors = 0;
733                                         int neighbor_lo_index[3];
734                                         int neighbor_hi_index[3];
735                                         int n;
736
737                                         /* check for upper bounds in advance
738                                          * to get the correct number of neighbors,
739                                          * needed for the diagonal element
740                                          */
741                                         if (!NEIGHBOR_MARGIN_k0 && (vert - stride2)->density > density_threshold)
742                                                 neighbor_lo_index[neighbors_lo++] = u - strideA2;
743                                         if (!NEIGHBOR_MARGIN_j0 && (vert - stride1)->density > density_threshold)
744                                                 neighbor_lo_index[neighbors_lo++] = u - strideA1;
745                                         if (!NEIGHBOR_MARGIN_i0 && (vert - stride0)->density > density_threshold)
746                                                 neighbor_lo_index[neighbors_lo++] = u - strideA0;
747                                         if (!NEIGHBOR_MARGIN_i1 && (vert + stride0)->density > density_threshold)
748                                                 neighbor_hi_index[neighbors_hi++] = u + strideA0;
749                                         if (!NEIGHBOR_MARGIN_j1 && (vert + stride1)->density > density_threshold)
750                                                 neighbor_hi_index[neighbors_hi++] = u + strideA1;
751                                         if (!NEIGHBOR_MARGIN_k1 && (vert + stride2)->density > density_threshold)
752                                                 neighbor_hi_index[neighbors_hi++] = u + strideA2;
753
754                                         /*int liquid_neighbors = neighbors_lo + neighbors_hi;*/
755                                         non_solid_neighbors = 6;
756
757                                         for (n = 0; n < neighbors_lo; ++n)
758                                                 A.insert(neighbor_lo_index[n], u) = -1.0f;
759                                         A.insert(u, u) = (float)non_solid_neighbors;
760                                         for (n = 0; n < neighbors_hi; ++n)
761                                                 A.insert(neighbor_hi_index[n], u) = -1.0f;
762                                 }
763                                 else {
764                                         A.insert(u, u) = 1.0f;
765                                 }
766                         }
767                 }
768         }
769
770         ConjugateGradient cg;
771         cg.setMaxIterations(100);
772         cg.setTolerance(0.01f);
773
774         cg.compute(A);
775
776         lVector p = cg.solve(B);
777
778         if (cg.info() == Eigen::Success) {
779                 /* Calculate velocity = grad(p) */
780                 for (k = 0; k < resA[2]; ++k) {
781                         for (j = 0; j < resA[1]; ++j) {
782                                 for (i = 0; i < resA[0]; ++i) {
783                                         int u = i * strideA0 + j * strideA1 + k * strideA2;
784                                         bool is_margin = MARGIN_i0 || MARGIN_i1 || MARGIN_j0 || MARGIN_j1 || MARGIN_k0 || MARGIN_k1;
785                                         if (is_margin)
786                                                 continue;
787
788                                         vert = vert_start + i * stride0 + j * stride1 + k * stride2;
789                                         if (vert->density > density_threshold) {
790                                                 float p_left   = p[u - strideA0];
791                                                 float p_right  = p[u + strideA0];
792                                                 float p_down   = p[u - strideA1];
793                                                 float p_up     = p[u + strideA1];
794                                                 float p_bottom = p[u - strideA2];
795                                                 float p_top    = p[u + strideA2];
796
797                                                 /* finite difference estimate of pressure gradient */
798                                                 float dvel[3];
799                                                 dvel[0] = p_right - p_left;
800                                                 dvel[1] = p_up - p_down;
801                                                 dvel[2] = p_top - p_bottom;
802                                                 mul_v3_fl(dvel, -0.5f * inv_flowfac);
803
804                                                 /* pressure gradient describes velocity delta */
805                                                 add_v3_v3v3(vert->velocity_smooth, vert->velocity, dvel);
806                                         }
807                                         else {
808                                                 zero_v3(vert->velocity_smooth);
809                                         }
810                                 }
811                         }
812                 }
813
814 #if 0
815                 {
816                         int axis = 0;
817                         float offset = 0.0f;
818
819                         int slice = (offset - grid->gmin[axis]) / grid->cellsize;
820
821                         for (k = 0; k < resA[2]; ++k) {
822                                 for (j = 0; j < resA[1]; ++j) {
823                                         for (i = 0; i < resA[0]; ++i) {
824                                                 int u = i * strideA0 + j * strideA1 + k * strideA2;
825                                                 bool is_margin = MARGIN_i0 || MARGIN_i1 || MARGIN_j0 || MARGIN_j1 || MARGIN_k0 || MARGIN_k1;
826                                                 if (i != slice)
827                                                         continue;
828
829                                                 vert = vert_start + i * stride0 + j * stride1 + k * stride2;
830
831                                                 float wloc[3], loc[3];
832                                                 float col0[3] = {0.0, 0.0, 0.0};
833                                                 float colp[3] = {0.0, 1.0, 1.0};
834                                                 float coln[3] = {1.0, 0.0, 1.0};
835                                                 float col[3];
836                                                 float fac;
837
838                                                 loc[0] = (float)(i - 1);
839                                                 loc[1] = (float)(j - 1);
840                                                 loc[2] = (float)(k - 1);
841                                                 grid_to_world(grid, wloc, loc);
842
843                                                 float pressure = p[u];
844                                                 if (pressure > 0.0f) {
845                                                         fac = CLAMPIS(pressure * grid->debug1, 0.0, 1.0);
846                                                         interp_v3_v3v3(col, col0, colp, fac);
847                                                 }
848                                                 else {
849                                                         fac = CLAMPIS(-pressure * grid->debug1, 0.0, 1.0);
850                                                         interp_v3_v3v3(col, col0, coln, fac);
851                                                 }
852                                                 if (fac > 0.05f)
853                                                         BKE_sim_debug_data_add_circle(grid->debug_data, wloc, 0.01f, col[0], col[1], col[2], "grid", 5533, i, j, k);
854
855                                                 if (!is_margin) {
856                                                         float dvel[3];
857                                                         sub_v3_v3v3(dvel, vert->velocity_smooth, vert->velocity);
858 //                                                      BKE_sim_debug_data_add_vector(grid->debug_data, wloc, dvel, 1, 1, 1, "grid", 5566, i, j, k);
859                                                 }
860
861                                                 if (!is_margin) {
862                                                         float d = CLAMPIS(vert->density * grid->debug2, 0.0f, 1.0f);
863                                                         float col0[3] = {0.3, 0.3, 0.3};
864                                                         float colp[3] = {0.0, 0.0, 1.0};
865                                                         float col[3];
866
867                                                         interp_v3_v3v3(col, col0, colp, d);
868 //                                                      if (d > 0.05f)
869 //                                                              BKE_sim_debug_data_add_dot(grid->debug_data, wloc, col[0], col[1], col[2], "grid", 5544, i, j, k);
870                                                 }
871                                         }
872                                 }
873                         }
874                 }
875 #endif
876
877                 return true;
878         }
879         else {
880                 /* Clear result in case of error */
881                 for (i = 0, vert = grid->verts; i < num_cells; ++i, ++vert) {
882                         zero_v3(vert->velocity_smooth);
883                 }
884
885                 return false;
886         }
887 }
888
889 #if 0 /* XXX weighting is incorrect, disabled for now */
890 /* Velocity filter kernel
891  * See https://en.wikipedia.org/wiki/Filter_%28large_eddy_simulation%29
892  */
893
894 BLI_INLINE void hair_volume_filter_box_convolute(HairVertexGrid *grid, float invD, const int kernel_size[3], int i, int j, int k)
895 {
896         int res = grid->res;
897         int p, q, r;
898         int minp = max_ii(i - kernel_size[0], 0), maxp = min_ii(i + kernel_size[0], res-1);
899         int minq = max_ii(j - kernel_size[1], 0), maxq = min_ii(j + kernel_size[1], res-1);
900         int minr = max_ii(k - kernel_size[2], 0), maxr = min_ii(k + kernel_size[2], res-1);
901         int offset, kernel_offset, kernel_dq, kernel_dr;
902         HairGridVert *verts;
903         float *vel_smooth;
904
905         offset = i + (j + k*res)*res;
906         verts = grid->verts;
907         vel_smooth = verts[offset].velocity_smooth;
908
909         kernel_offset = minp + (minq + minr*res)*res;
910         kernel_dq = res;
911         kernel_dr = res * res;
912         for (r = minr; r <= maxr; ++r) {
913                 for (q = minq; q <= maxq; ++q) {
914                         for (p = minp; p <= maxp; ++p) {
915
916                                 madd_v3_v3fl(vel_smooth, verts[kernel_offset].velocity, invD);
917
918                                 kernel_offset += 1;
919                         }
920                         kernel_offset += kernel_dq;
921                 }
922                 kernel_offset += kernel_dr;
923         }
924 }
925
926 void BPH_hair_volume_vertex_grid_filter_box(HairVertexGrid *grid, int kernel_size)
927 {
928         int size = hair_grid_size(grid->res);
929         int kernel_sizev[3] = {kernel_size, kernel_size, kernel_size};
930         int tot;
931         float invD;
932         int i, j, k;
933
934         if (kernel_size <= 0)
935                 return;
936
937         tot = kernel_size * 2 + 1;
938         invD = 1.0f / (float)(tot*tot*tot);
939
940         /* clear values for convolution */
941         for (i = 0; i < size; ++i) {
942                 zero_v3(grid->verts[i].velocity_smooth);
943         }
944
945         for (i = 0; i < grid->res; ++i) {
946                 for (j = 0; j < grid->res; ++j) {
947                         for (k = 0; k < grid->res; ++k) {
948                                 hair_volume_filter_box_convolute(grid, invD, kernel_sizev, i, j, k);
949                         }
950                 }
951         }
952
953         /* apply as new velocity */
954         for (i = 0; i < size; ++i) {
955                 copy_v3_v3(grid->verts[i].velocity, grid->verts[i].velocity_smooth);
956         }
957 }
958 #endif
959
960 HairGrid *BPH_hair_volume_create_vertex_grid(float cellsize, const float gmin[3], const float gmax[3])
961 {
962         float scale;
963         float extent[3];
964         int resmin[3], resmax[3], res[3];
965         float gmin_margin[3], gmax_margin[3];
966         int size;
967         HairGrid *grid;
968         int i;
969
970         /* sanity check */
971         if (cellsize <= 0.0f)
972                 cellsize = 1.0f;
973         scale = 1.0f / cellsize;
974
975         sub_v3_v3v3(extent, gmax, gmin);
976         for (i = 0; i < 3; ++i) {
977                 resmin[i] = floor_int(gmin[i] * scale);
978                 resmax[i] = floor_int(gmax[i] * scale) + 1;
979
980                 /* add margin of 1 cell */
981                 resmin[i] -= 1;
982                 resmax[i] += 1;
983
984                 res[i] = resmax[i] - resmin[i] + 1;
985                 /* sanity check: avoid null-sized grid */
986                 if (res[i] < 4) {
987                         res[i] = 4;
988                         resmax[i] = resmin[i] + 4;
989                 }
990                 /* sanity check: avoid too large grid size */
991                 if (res[i] > MAX_HAIR_GRID_RES) {
992                         res[i] = MAX_HAIR_GRID_RES;
993                         resmax[i] = resmin[i] + MAX_HAIR_GRID_RES;
994                 }
995
996                 gmin_margin[i] = (float)resmin[i] * cellsize;
997                 gmax_margin[i] = (float)resmax[i] * cellsize;
998         }
999         size = hair_grid_size(res);
1000
1001         grid = (HairGrid *)MEM_callocN(sizeof(HairGrid), "hair grid");
1002         grid->res[0] = res[0];
1003         grid->res[1] = res[1];
1004         grid->res[2] = res[2];
1005         copy_v3_v3(grid->gmin, gmin_margin);
1006         copy_v3_v3(grid->gmax, gmax_margin);
1007         grid->cellsize = cellsize;
1008         grid->inv_cellsize = scale;
1009         grid->verts = (HairGridVert *)MEM_callocN(sizeof(HairGridVert) * size, "hair voxel data");
1010
1011         return grid;
1012 }
1013
1014 void BPH_hair_volume_free_vertex_grid(HairGrid *grid)
1015 {
1016         if (grid) {
1017                 if (grid->verts)
1018                         MEM_freeN(grid->verts);
1019                 MEM_freeN(grid);
1020         }
1021 }
1022
1023 void BPH_hair_volume_grid_geometry(HairGrid *grid, float *cellsize, int res[3], float gmin[3], float gmax[3])
1024 {
1025         if (cellsize) *cellsize = grid->cellsize;
1026         if (res) copy_v3_v3_int(res, grid->res);
1027         if (gmin) copy_v3_v3(gmin, grid->gmin);
1028         if (gmax) copy_v3_v3(gmax, grid->gmax);
1029 }
1030
1031 #if 0
1032 static HairGridVert *hair_volume_create_collision_grid(ClothModifierData *clmd, lfVector *lX, unsigned int numverts)
1033 {
1034         int res = hair_grid_res;
1035         int size = hair_grid_size(res);
1036         HairGridVert *collgrid;
1037         ListBase *colliders;
1038         ColliderCache *col = NULL;
1039         float gmin[3], gmax[3], scale[3];
1040         /* 2.0f is an experimental value that seems to give good results */
1041         float collfac = 2.0f * clmd->sim_parms->collider_friction;
1042         unsigned int    v = 0;
1043         int                 i = 0;
1044
1045         hair_volume_get_boundbox(lX, numverts, gmin, gmax);
1046         hair_grid_get_scale(res, gmin, gmax, scale);
1047
1048         collgrid = MEM_mallocN(sizeof(HairGridVert) * size, "hair collider voxel data");
1049
1050         /* initialize grid */
1051         for (i = 0; i < size; ++i) {
1052                 zero_v3(collgrid[i].velocity);
1053                 collgrid[i].density = 0.0f;
1054         }
1055
1056         /* gather colliders */
1057         colliders = get_collider_cache(clmd->scene, NULL, NULL);
1058         if (colliders && collfac > 0.0f) {
1059                 for (col = colliders->first; col; col = col->next) {
1060                         MVert *loc0 = col->collmd->x;
1061                         MVert *loc1 = col->collmd->xnew;
1062                         float vel[3];
1063                         float weights[8];
1064                         int di, dj, dk;
1065
1066                         for (v=0; v < col->collmd->numverts; v++, loc0++, loc1++) {
1067                                 int offset;
1068
1069                                 if (!hair_grid_point_valid(loc1->co, gmin, gmax))
1070                                         continue;
1071
1072                                 offset = hair_grid_weights(res, gmin, scale, lX[v], weights);
1073
1074                                 sub_v3_v3v3(vel, loc1->co, loc0->co);
1075
1076                                 for (di = 0; di < 2; ++di) {
1077                                         for (dj = 0; dj < 2; ++dj) {
1078                                                 for (dk = 0; dk < 2; ++dk) {
1079                                                         int voffset = offset + di + (dj + dk*res)*res;
1080                                                         int iw = di + dj*2 + dk*4;
1081
1082                                                         collgrid[voffset].density += weights[iw];
1083                                                         madd_v3_v3fl(collgrid[voffset].velocity, vel, weights[iw]);
1084                                                 }
1085                                         }
1086                                 }
1087                         }
1088                 }
1089         }
1090         free_collider_cache(&colliders);
1091
1092         /* divide velocity with density */
1093         for (i = 0; i < size; i++) {
1094                 float density = collgrid[i].density;
1095                 if (density > 0.0f)
1096                         mul_v3_fl(collgrid[i].velocity, 1.0f/density);
1097         }
1098
1099         return collgrid;
1100 }
1101 #endif
1102
1103 bool BPH_hair_volume_get_texture_data(HairGrid *grid, VoxelData *vd)
1104 {
1105         int totres, i;
1106         int depth;
1107
1108         vd->resol[0] = grid->res[0];
1109         vd->resol[1] = grid->res[1];
1110         vd->resol[2] = grid->res[2];
1111
1112         totres = hair_grid_size(grid->res);
1113
1114         if (vd->hair_type == TEX_VD_HAIRVELOCITY) {
1115                 depth = 4;
1116                 vd->data_type = TEX_VD_RGBA_PREMUL;
1117         }
1118         else {
1119                 depth = 1;
1120                 vd->data_type = TEX_VD_INTENSITY;
1121         }
1122
1123         if (totres > 0) {
1124                 vd->dataset = (float *)MEM_mapallocN(sizeof(float) * depth * (totres), "hair volume texture data");
1125
1126                 for (i = 0; i < totres; ++i) {
1127                         switch (vd->hair_type) {
1128                                 case TEX_VD_HAIRDENSITY:
1129                                         vd->dataset[i] = grid->verts[i].density;
1130                                         break;
1131
1132                                 case TEX_VD_HAIRRESTDENSITY:
1133                                         vd->dataset[i] = 0.0f; // TODO
1134                                         break;
1135
1136                                 case TEX_VD_HAIRVELOCITY: {
1137                                         vd->dataset[i + 0*totres] = grid->verts[i].velocity[0];
1138                                         vd->dataset[i + 1*totres] = grid->verts[i].velocity[1];
1139                                         vd->dataset[i + 2*totres] = grid->verts[i].velocity[2];
1140                                         vd->dataset[i + 3*totres] = len_v3(grid->verts[i].velocity);
1141                                         break;
1142                                 }
1143                                 case TEX_VD_HAIRENERGY:
1144                                         vd->dataset[i] = 0.0f; // TODO
1145                                         break;
1146                         }
1147                 }
1148         }
1149         else {
1150                 vd->dataset = NULL;
1151         }
1152
1153         return true;
1154 }