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