merge from master
[blender.git] / source / blender / blenkernel / intern / particle_distribute.c
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) 2007 by Janne Karhu.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Raul Fernandez Hernandez (Farsthary),
24  *                 Stephen Swhitehorn,
25  *                 Lukas Toenne
26  *
27  * ***** END GPL LICENSE BLOCK *****
28  */
29
30 /** \file blender/blenkernel/intern/particle_distribute.c
31  *  \ingroup bke
32  */
33
34 #include <string.h>
35
36 #include "MEM_guardedalloc.h"
37
38 #include "BLI_utildefines.h"
39 #include "BLI_jitter_2d.h"
40 #include "BLI_kdtree.h"
41 #include "BLI_math.h"
42 #include "BLI_math_geom.h"
43 #include "BLI_rand.h"
44 #include "BLI_sort.h"
45 #include "BLI_task.h"
46
47 #include "DNA_mesh_types.h"
48 #include "DNA_meshdata_types.h"
49 #include "DNA_modifier_types.h"
50 #include "DNA_particle_types.h"
51 #include "DNA_scene_types.h"
52
53 #include "BKE_cdderivedmesh.h"
54 #include "BKE_DerivedMesh.h"
55 #include "BKE_global.h"
56 #include "BKE_mesh.h"
57 #include "BKE_object.h"
58 #include "BKE_particle.h"
59
60 static int psys_render_simplify_distribution(ParticleThreadContext *ctx, int tot);
61
62 static void alloc_child_particles(ParticleSystem *psys, int tot)
63 {
64         if (psys->child) {
65                 /* only re-allocate if we have to */
66                 if (psys->part->childtype && psys->totchild == tot) {
67                         memset(psys->child, 0, tot*sizeof(ChildParticle));
68                         return;
69                 }
70
71                 MEM_freeN(psys->child);
72                 psys->child=NULL;
73                 psys->totchild=0;
74         }
75
76         if (psys->part->childtype) {
77                 psys->totchild= tot;
78                 if (psys->totchild)
79                         psys->child= MEM_callocN(psys->totchild*sizeof(ChildParticle), "child_particles");
80         }
81 }
82
83 static void distribute_simple_children(Scene *scene, Object *ob, DerivedMesh *finaldm, DerivedMesh *deformdm, ParticleSystem *psys)
84 {
85         ChildParticle *cpa = NULL;
86         int i, p;
87         int child_nbr= psys_get_child_number(scene, psys);
88         int totpart= psys_get_tot_child(scene, psys);
89
90         alloc_child_particles(psys, totpart);
91
92         cpa = psys->child;
93         for (i=0; i<child_nbr; i++) {
94                 for (p=0; p<psys->totpart; p++,cpa++) {
95                         float length=2.0;
96                         cpa->parent=p;
97                                         
98                         /* create even spherical distribution inside unit sphere */
99                         while (length>=1.0f) {
100                                 cpa->fuv[0]=2.0f*BLI_frand()-1.0f;
101                                 cpa->fuv[1]=2.0f*BLI_frand()-1.0f;
102                                 cpa->fuv[2]=2.0f*BLI_frand()-1.0f;
103                                 length=len_v3(cpa->fuv);
104                         }
105
106                         cpa->num=-1;
107                 }
108         }
109         /* dmcache must be updated for parent particles if children from faces is used */
110         psys_calc_dmcache(ob, finaldm, deformdm, psys);
111 }
112 static void distribute_grid(DerivedMesh *dm, ParticleSystem *psys)
113 {
114         ParticleData *pa=NULL;
115         float min[3], max[3], delta[3], d;
116         MVert *mv, *mvert = dm->getVertDataArray(dm,0);
117         int totvert=dm->getNumVerts(dm), from=psys->part->from;
118         int i, j, k, p, res=psys->part->grid_res, size[3], axis;
119
120         /* find bounding box of dm */
121         if (totvert > 0) {
122                 mv=mvert;
123                 copy_v3_v3(min, mv->co);
124                 copy_v3_v3(max, mv->co);
125                 mv++;
126                 for (i = 1; i < totvert; i++, mv++) {
127                         minmax_v3v3_v3(min, max, mv->co);
128                 }
129         }
130         else {
131                 zero_v3(min);
132                 zero_v3(max);
133         }
134
135         sub_v3_v3v3(delta, max, min);
136
137         /* determine major axis */
138         axis = axis_dominant_v3_single(delta);
139          
140         d = delta[axis]/(float)res;
141
142         size[axis] = res;
143         size[(axis+1)%3] = (int)ceil(delta[(axis+1)%3]/d);
144         size[(axis+2)%3] = (int)ceil(delta[(axis+2)%3]/d);
145
146         /* float errors grrr.. */
147         size[(axis+1)%3] = MIN2(size[(axis+1)%3],res);
148         size[(axis+2)%3] = MIN2(size[(axis+2)%3],res);
149
150         size[0] = MAX2(size[0], 1);
151         size[1] = MAX2(size[1], 1);
152         size[2] = MAX2(size[2], 1);
153
154         /* no full offset for flat/thin objects */
155         min[0]+= d < delta[0] ? d/2.f : delta[0]/2.f;
156         min[1]+= d < delta[1] ? d/2.f : delta[1]/2.f;
157         min[2]+= d < delta[2] ? d/2.f : delta[2]/2.f;
158
159         for (i=0,p=0,pa=psys->particles; i<res; i++) {
160                 for (j=0; j<res; j++) {
161                         for (k=0; k<res; k++,p++,pa++) {
162                                 pa->fuv[0] = min[0] + (float)i*d;
163                                 pa->fuv[1] = min[1] + (float)j*d;
164                                 pa->fuv[2] = min[2] + (float)k*d;
165                                 pa->flag |= PARS_UNEXIST;
166                                 pa->hair_index = 0; /* abused in volume calculation */
167                         }
168                 }
169         }
170
171         /* enable particles near verts/edges/faces/inside surface */
172         if (from==PART_FROM_VERT) {
173                 float vec[3];
174
175                 pa=psys->particles;
176
177                 min[0] -= d/2.0f;
178                 min[1] -= d/2.0f;
179                 min[2] -= d/2.0f;
180
181                 for (i=0,mv=mvert; i<totvert; i++,mv++) {
182                         sub_v3_v3v3(vec,mv->co,min);
183                         vec[0]/=delta[0];
184                         vec[1]/=delta[1];
185                         vec[2]/=delta[2];
186                         pa[((int)(vec[0] * (size[0] - 1))  * res +
187                             (int)(vec[1] * (size[1] - 1))) * res +
188                             (int)(vec[2] * (size[2] - 1))].flag &= ~PARS_UNEXIST;
189                 }
190         }
191         else if (ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
192                 float co1[3], co2[3];
193
194                 MFace *mface= NULL, *mface_array;
195                 float v1[3], v2[3], v3[3], v4[4], lambda;
196                 int a, a1, a2, a0mul, a1mul, a2mul, totface;
197                 int amax= from==PART_FROM_FACE ? 3 : 1;
198
199                 totface=dm->getNumTessFaces(dm);
200                 mface=mface_array=dm->getTessFaceDataArray(dm,CD_MFACE);
201                 
202                 for (a=0; a<amax; a++) {
203                         if (a==0) { a0mul=res*res; a1mul=res; a2mul=1; }
204                         else if (a==1) { a0mul=res; a1mul=1; a2mul=res*res; }
205                         else { a0mul=1; a1mul=res*res; a2mul=res; }
206
207                         for (a1=0; a1<size[(a+1)%3]; a1++) {
208                                 for (a2=0; a2<size[(a+2)%3]; a2++) {
209                                         mface= mface_array;
210
211                                         pa = psys->particles + a1*a1mul + a2*a2mul;
212                                         copy_v3_v3(co1, pa->fuv);
213                                         co1[a] -= d < delta[a] ? d/2.f : delta[a]/2.f;
214                                         copy_v3_v3(co2, co1);
215                                         co2[a] += delta[a] + 0.001f*d;
216                                         co1[a] -= 0.001f*d;
217
218                                         struct IsectRayPrecalc isect_precalc;
219                                         float ray_direction[3];
220                                         sub_v3_v3v3(ray_direction, co2, co1);
221                                         isect_ray_tri_watertight_v3_precalc(&isect_precalc, ray_direction);
222
223                                         /* lets intersect the faces */
224                                         for (i=0; i<totface; i++,mface++) {
225                                                 copy_v3_v3(v1, mvert[mface->v1].co);
226                                                 copy_v3_v3(v2, mvert[mface->v2].co);
227                                                 copy_v3_v3(v3, mvert[mface->v3].co);
228
229                                                 bool intersects_tri = isect_ray_tri_watertight_v3(co1,
230                                                                                                   &isect_precalc,
231                                                                                                   v1, v2, v3,
232                                                                                                   &lambda, NULL);
233                                                 if (intersects_tri) {
234                                                         if (from==PART_FROM_FACE)
235                                                                 (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
236                                                         else /* store number of intersections */
237                                                                 (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
238                                                 }
239
240                                                 if (mface->v4 && (!intersects_tri || from==PART_FROM_VOLUME)) {
241                                                         copy_v3_v3(v4, mvert[mface->v4].co);
242
243                                                         if (isect_ray_tri_watertight_v3(co1,
244                                                                                         &isect_precalc,
245                                                                                         v1, v3, v4,
246                                                                                         &lambda, NULL)) {
247                                                                 if (from==PART_FROM_FACE)
248                                                                         (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
249                                                                 else
250                                                                         (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
251                                                         }
252                                                 }
253                                         }
254
255                                         if (from==PART_FROM_VOLUME) {
256                                                 int in=pa->hair_index%2;
257                                                 if (in) pa->hair_index++;
258                                                 for (i=0; i<size[0]; i++) {
259                                                         if (in || (pa+i*a0mul)->hair_index%2)
260                                                                 (pa+i*a0mul)->flag &= ~PARS_UNEXIST;
261                                                         /* odd intersections == in->out / out->in */
262                                                         /* even intersections -> in stays same */
263                                                         in=(in + (pa+i*a0mul)->hair_index) % 2;
264                                                 }
265                                         }
266                                 }
267                         }
268                 }
269         }
270
271         if (psys->part->flag & PART_GRID_HEXAGONAL) {
272                 for (i=0,p=0,pa=psys->particles; i<res; i++) {
273                         for (j=0; j<res; j++) {
274                                 for (k=0; k<res; k++,p++,pa++) {
275                                         if (j%2)
276                                                 pa->fuv[0] += d/2.f;
277
278                                         if (k%2) {
279                                                 pa->fuv[0] += d/2.f;
280                                                 pa->fuv[1] += d/2.f;
281                                         }
282                                 }
283                         }
284                 }
285         }
286
287         if (psys->part->flag & PART_GRID_INVERT) {
288                 for (i=0; i<size[0]; i++) {
289                         for (j=0; j<size[1]; j++) {
290                                 pa=psys->particles + res*(i*res + j);
291                                 for (k=0; k<size[2]; k++, pa++) {
292                                         pa->flag ^= PARS_UNEXIST;
293                                 }
294                         }
295                 }
296         }
297
298         if (psys->part->grid_rand > 0.f) {
299                 float rfac = d * psys->part->grid_rand;
300                 for (p=0,pa=psys->particles; p<psys->totpart; p++,pa++) {
301                         if (pa->flag & PARS_UNEXIST)
302                                 continue;
303
304                         pa->fuv[0] += rfac * (psys_frand(psys, p + 31) - 0.5f);
305                         pa->fuv[1] += rfac * (psys_frand(psys, p + 32) - 0.5f);
306                         pa->fuv[2] += rfac * (psys_frand(psys, p + 33) - 0.5f);
307                 }
308         }
309 }
310
311 /* modified copy from rayshade.c */
312 static void hammersley_create(float *out, int n, int seed, float amount)
313 {
314         RNG *rng;
315
316         double offs[2], t;
317
318         rng = BLI_rng_new(31415926 + n + seed);
319         offs[0] = BLI_rng_get_double(rng) + (double)amount;
320         offs[1] = BLI_rng_get_double(rng) + (double)amount;
321         BLI_rng_free(rng);
322
323         for (int k = 0; k < n; k++) {
324                 BLI_hammersley_1D(k, &t);
325
326                 out[2*k + 0] = fmod((double)k/(double)n + offs[0], 1.0);
327                 out[2*k + 1] = fmod(t + offs[1], 1.0);
328         }
329 }
330
331 /* almost exact copy of BLI_jitter_init */
332 static void init_mv_jit(float *jit, int num, int seed2, float amount)
333 {
334         RNG *rng;
335         float *jit2, x, rad1, rad2, rad3;
336         int i, num2;
337
338         if (num==0) return;
339
340         rad1= (float)(1.0f/sqrtf((float)num));
341         rad2= (float)(1.0f/((float)num));
342         rad3= (float)sqrtf((float)num)/((float)num);
343
344         rng = BLI_rng_new(31415926 + num + seed2);
345         x= 0;
346         num2 = 2 * num;
347         for (i=0; i<num2; i+=2) {
348         
349                 jit[i] = x + amount*rad1*(0.5f - BLI_rng_get_float(rng));
350                 jit[i+1] = i/(2.0f*num) + amount*rad1*(0.5f - BLI_rng_get_float(rng));
351                 
352                 jit[i]-= (float)floor(jit[i]);
353                 jit[i+1]-= (float)floor(jit[i+1]);
354                 
355                 x+= rad3;
356                 x -= (float)floor(x);
357         }
358
359         jit2= MEM_mallocN(12 + 2*sizeof(float)*num, "initjit");
360
361         for (i=0 ; i<4 ; i++) {
362                 BLI_jitterate1((float (*)[2])jit, (float (*)[2])jit2, num, rad1);
363                 BLI_jitterate1((float (*)[2])jit, (float (*)[2])jit2, num, rad1);
364                 BLI_jitterate2((float (*)[2])jit, (float (*)[2])jit2, num, rad2);
365         }
366         MEM_freeN(jit2);
367         BLI_rng_free(rng);
368 }
369
370 static void psys_uv_to_w(float u, float v, int quad, float *w)
371 {
372         float vert[4][3], co[3];
373
374         if (!quad) {
375                 if (u+v > 1.0f)
376                         v= 1.0f-v;
377                 else
378                         u= 1.0f-u;
379         }
380
381         vert[0][0] = 0.0f; vert[0][1] = 0.0f; vert[0][2] = 0.0f;
382         vert[1][0] = 1.0f; vert[1][1] = 0.0f; vert[1][2] = 0.0f;
383         vert[2][0] = 1.0f; vert[2][1] = 1.0f; vert[2][2] = 0.0f;
384
385         co[0] = u;
386         co[1] = v;
387         co[2] = 0.0f;
388
389         if (quad) {
390                 vert[3][0] = 0.0f; vert[3][1] = 1.0f; vert[3][2] = 0.0f;
391                 interp_weights_poly_v3( w,vert, 4, co);
392         }
393         else {
394                 interp_weights_poly_v3( w,vert, 3, co);
395                 w[3] = 0.0f;
396         }
397 }
398
399 /* Find the index in "sum" array before "value" is crossed. */
400 static int distribute_binary_search(float *sum, int n, float value)
401 {
402         int mid, low = 0, high = n - 1;
403
404         if (high == low)
405                 return low;
406
407         if (sum[low] >= value)
408                 return low;
409
410         if (sum[high - 1] < value)
411                 return high;
412
413         while (low < high) {
414                 mid = (low + high) / 2;
415                 
416                 if ((sum[mid] >= value) && (sum[mid - 1] < value))
417                         return mid;
418                 
419                 if (sum[mid] > value) {
420                         high = mid - 1;
421                 }
422                 else {
423                         low = mid + 1;
424                 }
425         }
426
427         return low;
428 }
429
430 /* the max number if calls to rng_* funcs within psys_thread_distribute_particle
431  * be sure to keep up to date if this changes */
432 #define PSYS_RND_DIST_SKIP 2
433
434 /* note: this function must be thread safe, for from == PART_FROM_CHILD */
435 #define ONLY_WORKING_WITH_PA_VERTS 0
436 static void distribute_from_verts_exec(ParticleTask *thread, ParticleData *pa, int p)
437 {
438         ParticleThreadContext *ctx= thread->ctx;
439         MFace *mface;
440
441         mface = ctx->dm->getTessFaceDataArray(ctx->dm, CD_MFACE);
442
443         int rng_skip_tot = PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
444
445         /* TODO_PARTICLE - use original index */
446         pa->num = ctx->index[p];
447
448         zero_v4(pa->fuv);
449
450         if (pa->num != DMCACHE_NOTFOUND && pa->num < ctx->dm->getNumVerts(ctx->dm)) {
451
452                 /* This finds the first face to contain the emitting vertex,
453                  * this is not ideal, but is mostly fine as UV seams generally
454                  * map to equal-colored parts of a texture */
455                 for (int i = 0; i < ctx->dm->getNumTessFaces(ctx->dm); i++, mface++) {
456                         if (ELEM(pa->num, mface->v1, mface->v2, mface->v3, mface->v4)) {
457                                 unsigned int *vert = &mface->v1;
458
459                                 for (int j = 0; j < 4; j++, vert++) {
460                                         if (*vert == pa->num) {
461                                                 pa->fuv[j] = 1.0f;
462                                                 break;
463                                         }
464                                 }
465
466                                 break;
467                         }
468                 }
469         }
470         
471 #if ONLY_WORKING_WITH_PA_VERTS
472         if (ctx->tree) {
473                 KDTreeNearest ptn[3];
474                 int w, maxw;
475                 
476                 psys_particle_on_dm(ctx->dm,from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co1,0,0,0,orco1,0);
477                 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &orco1, 1, 1);
478                 maxw = BLI_kdtree_find_nearest_n(ctx->tree,orco1,ptn,3);
479                 
480                 for (w=0; w<maxw; w++) {
481                         pa->verts[w]=ptn->num;
482                 }
483         }
484 #endif
485         
486         if (rng_skip_tot > 0) /* should never be below zero */
487                 BLI_rng_skip(thread->rng, rng_skip_tot);
488 }
489
490 static void distribute_from_faces_exec(ParticleTask *thread, ParticleData *pa, int p) {
491         ParticleThreadContext *ctx= thread->ctx;
492         DerivedMesh *dm= ctx->dm;
493         float randu, randv;
494         int distr= ctx->distr;
495         int i;
496         int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
497
498         MFace *mface;
499         
500         pa->num = i = ctx->index[p];
501         mface = dm->getTessFaceData(dm,i,CD_MFACE);
502         
503         switch (distr) {
504                 case PART_DISTR_JIT:
505                         if (ctx->jitlevel == 1) {
506                                 if (mface->v4)
507                                         psys_uv_to_w(0.5f, 0.5f, mface->v4, pa->fuv);
508                                 else
509                                         psys_uv_to_w(1.0f / 3.0f, 1.0f / 3.0f, mface->v4, pa->fuv);
510                         }
511                         else {
512                                 float offset = fmod(ctx->jitoff[i] + (float)p, (float)ctx->jitlevel);
513                                 if (!isnan(offset)) {
514                                         psys_uv_to_w(ctx->jit[2*(int)offset], ctx->jit[2*(int)offset+1], mface->v4, pa->fuv);
515                                 }
516                         }
517                         break;
518                 case PART_DISTR_RAND:
519                         randu= BLI_rng_get_float(thread->rng);
520                         randv= BLI_rng_get_float(thread->rng);
521                         rng_skip_tot -= 2;
522                         
523                         psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
524                         break;
525         }
526         pa->foffset= 0.0f;
527         
528         if (rng_skip_tot > 0) /* should never be below zero */
529                 BLI_rng_skip(thread->rng, rng_skip_tot);
530 }
531
532 static void distribute_from_volume_exec(ParticleTask *thread, ParticleData *pa, int p) {
533         ParticleThreadContext *ctx= thread->ctx;
534         DerivedMesh *dm= ctx->dm;
535         float *v1, *v2, *v3, *v4, nor[3], co[3];
536         float cur_d, min_d, randu, randv;
537         int distr= ctx->distr;
538         int i, intersect, tot;
539         int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
540         
541         MFace *mface;
542         MVert *mvert=dm->getVertDataArray(dm,CD_MVERT);
543         
544         pa->num = i = ctx->index[p];
545         mface = dm->getTessFaceData(dm,i,CD_MFACE);
546         
547         switch (distr) {
548                 case PART_DISTR_JIT:
549                         if (ctx->jitlevel == 1) {
550                                 if (mface->v4)
551                                         psys_uv_to_w(0.5f, 0.5f, mface->v4, pa->fuv);
552                                 else
553                                         psys_uv_to_w(1.0f / 3.0f, 1.0f / 3.0f, mface->v4, pa->fuv);
554                         }
555                         else {
556                                 float offset = fmod(ctx->jitoff[i] + (float)p, (float)ctx->jitlevel);
557                                 if (!isnan(offset)) {
558                                         psys_uv_to_w(ctx->jit[2*(int)offset], ctx->jit[2*(int)offset+1], mface->v4, pa->fuv);
559                                 }
560                         }
561                         break;
562                 case PART_DISTR_RAND:
563                         randu= BLI_rng_get_float(thread->rng);
564                         randv= BLI_rng_get_float(thread->rng);
565                         rng_skip_tot -= 2;
566                         
567                         psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
568                         break;
569         }
570         pa->foffset= 0.0f;
571         
572         /* experimental */
573         tot=dm->getNumTessFaces(dm);
574         
575         psys_interpolate_face(mvert,mface,0,0,pa->fuv,co,nor,0,0,0,0);
576         
577         normalize_v3(nor);
578         negate_v3(nor);
579         
580         min_d=FLT_MAX;
581         intersect=0;
582         
583         for (i=0,mface=dm->getTessFaceDataArray(dm,CD_MFACE); i<tot; i++,mface++) {
584                 if (i==pa->num) continue;
585                 
586                 v1=mvert[mface->v1].co;
587                 v2=mvert[mface->v2].co;
588                 v3=mvert[mface->v3].co;
589                 
590                 if (isect_ray_tri_v3(co, nor, v2, v3, v1, &cur_d, NULL)) {
591                         if (cur_d<min_d) {
592                                 min_d=cur_d;
593                                 pa->foffset=cur_d*0.5f; /* to the middle of volume */
594                                 intersect=1;
595                         }
596                 }
597                 if (mface->v4) {
598                         v4=mvert[mface->v4].co;
599                         
600                         if (isect_ray_tri_v3(co, nor, v4, v1, v3, &cur_d, NULL)) {
601                                 if (cur_d<min_d) {
602                                         min_d=cur_d;
603                                         pa->foffset=cur_d*0.5f; /* to the middle of volume */
604                                         intersect=1;
605                                 }
606                         }
607                 }
608         }
609         if (intersect==0)
610                 pa->foffset=0.0;
611         else {
612                 switch (distr) {
613                         case PART_DISTR_JIT:
614                                 pa->foffset *= ctx->jit[p % (2 * ctx->jitlevel)];
615                                 break;
616                         case PART_DISTR_RAND:
617                                 pa->foffset *= BLI_frand();
618                                 break;
619                 }
620         }
621         
622         if (rng_skip_tot > 0) /* should never be below zero */
623                 BLI_rng_skip(thread->rng, rng_skip_tot);
624 }
625
626 static void distribute_children_exec(ParticleTask *thread, ChildParticle *cpa, int p) {
627         ParticleThreadContext *ctx= thread->ctx;
628         Object *ob= ctx->sim.ob;
629         DerivedMesh *dm= ctx->dm;
630         float orco1[3], co1[3], nor1[3];
631         float randu, randv;
632         int cfrom= ctx->cfrom;
633         int i;
634         int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
635         
636         MFace *mf;
637         
638         if (ctx->index[p] < 0) {
639                 cpa->num=0;
640                 cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]=0.0f;
641                 cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
642                 return;
643         }
644         
645         mf= dm->getTessFaceData(dm, ctx->index[p], CD_MFACE);
646         
647         randu= BLI_rng_get_float(thread->rng);
648         randv= BLI_rng_get_float(thread->rng);
649         rng_skip_tot -= 2;
650         
651         psys_uv_to_w(randu, randv, mf->v4, cpa->fuv);
652         
653         cpa->num = ctx->index[p];
654         
655         if (ctx->tree) {
656                 KDTreeNearest ptn[10];
657                 int w,maxw;//, do_seams;
658                 float maxd /*, mind,dd */, totw= 0.0f;
659                 int parent[10];
660                 float pweight[10];
661                 
662                 psys_particle_on_dm(dm,cfrom,cpa->num,DMCACHE_ISCHILD,cpa->fuv,cpa->foffset,co1,nor1,NULL,NULL,orco1,NULL);
663                 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &orco1, 1, 1);
664                 maxw = BLI_kdtree_find_nearest_n(ctx->tree,orco1,ptn,3);
665                 
666                 maxd=ptn[maxw-1].dist;
667                 /* mind=ptn[0].dist; */ /* UNUSED */
668                 
669                 /* the weights here could be done better */
670                 for (w=0; w<maxw; w++) {
671                         parent[w]=ptn[w].index;
672                         pweight[w]=(float)pow(2.0,(double)(-6.0f*ptn[w].dist/maxd));
673                 }
674                 for (;w<10; w++) {
675                         parent[w]=-1;
676                         pweight[w]=0.0f;
677                 }
678                 
679                 for (w=0,i=0; w<maxw && i<4; w++) {
680                         if (parent[w]>=0) {
681                                 cpa->pa[i]=parent[w];
682                                 cpa->w[i]=pweight[w];
683                                 totw+=pweight[w];
684                                 i++;
685                         }
686                 }
687                 for (;i<4; i++) {
688                         cpa->pa[i]=-1;
689                         cpa->w[i]=0.0f;
690                 }
691                 
692                 if (totw > 0.0f) {
693                         for (w = 0; w < 4; w++) {
694                                 cpa->w[w] /= totw;
695                         }
696                 }
697                 
698                 cpa->parent=cpa->pa[0];
699         }
700
701         if (rng_skip_tot > 0) /* should never be below zero */
702                 BLI_rng_skip(thread->rng, rng_skip_tot);
703 }
704
705 static void exec_distribute_parent(TaskPool * __restrict UNUSED(pool), void *taskdata, int UNUSED(threadid))
706 {
707         ParticleTask *task = taskdata;
708         ParticleSystem *psys= task->ctx->sim.psys;
709         ParticleData *pa;
710         int p;
711
712         BLI_rng_skip(task->rng, PSYS_RND_DIST_SKIP * task->begin);
713         
714         pa= psys->particles + task->begin;
715         switch (psys->part->from) {
716                 case PART_FROM_FACE:
717                         for (p = task->begin; p < task->end; ++p, ++pa)
718                                 distribute_from_faces_exec(task, pa, p);
719                         break;
720                 case PART_FROM_VOLUME:
721                         for (p = task->begin; p < task->end; ++p, ++pa)
722                                 distribute_from_volume_exec(task, pa, p);
723                         break;
724                 case PART_FROM_VERT:
725                         for (p = task->begin; p < task->end; ++p, ++pa)
726                                 distribute_from_verts_exec(task, pa, p);
727                         break;
728         }
729 }
730
731 static void exec_distribute_child(TaskPool * __restrict UNUSED(pool), void *taskdata, int UNUSED(threadid))
732 {
733         ParticleTask *task = taskdata;
734         ParticleSystem *psys = task->ctx->sim.psys;
735         ChildParticle *cpa;
736         int p;
737         
738         /* RNG skipping at the beginning */
739         cpa = psys->child;
740         for (p = 0; p < task->begin; ++p, ++cpa) {
741                 if (task->ctx->skip) /* simplification skip */
742                         BLI_rng_skip(task->rng, PSYS_RND_DIST_SKIP * task->ctx->skip[p]);
743                 
744                 BLI_rng_skip(task->rng, PSYS_RND_DIST_SKIP);
745         }
746                 
747         for (; p < task->end; ++p, ++cpa) {
748                 if (task->ctx->skip) /* simplification skip */
749                         BLI_rng_skip(task->rng, PSYS_RND_DIST_SKIP * task->ctx->skip[p]);
750                 
751                 distribute_children_exec(task, cpa, p);
752         }
753 }
754
755 static int distribute_compare_orig_index(const void *p1, const void *p2, void *user_data)
756 {
757         int *orig_index = (int *) user_data;
758         int index1 = orig_index[*(const int *)p1];
759         int index2 = orig_index[*(const int *)p2];
760
761         if (index1 < index2)
762                 return -1;
763         else if (index1 == index2) {
764                 /* this pointer comparison appears to make qsort stable for glibc,
765                  * and apparently on solaris too, makes the renders reproducible */
766                 if (p1 < p2)
767                         return -1;
768                 else if (p1 == p2)
769                         return 0;
770                 else
771                         return 1;
772         }
773         else
774                 return 1;
775 }
776
777 static void distribute_invalid(Scene *scene, ParticleSystem *psys, int from)
778 {
779         if (from == PART_FROM_CHILD) {
780                 ChildParticle *cpa;
781                 int p, totchild = psys_get_tot_child(scene, psys);
782
783                 if (psys->child && totchild) {
784                         for (p=0,cpa=psys->child; p<totchild; p++,cpa++) {
785                                 cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3] = 0.0;
786                                 cpa->foffset= 0.0f;
787                                 cpa->parent=0;
788                                 cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
789                                 cpa->num= -1;
790                         }
791                 }
792         }
793         else {
794                 PARTICLE_P;
795                 LOOP_PARTICLES {
796                         pa->fuv[0] = pa->fuv[1] = pa->fuv[2] = pa->fuv[3] = 0.0;
797                         pa->foffset= 0.0f;
798                         pa->num= -1;
799                 }
800         }
801 }
802
803 /* Creates a distribution of coordinates on a DerivedMesh       */
804 /* This is to denote functionality that does not yet work with mesh - only derived mesh */
805 static int psys_thread_context_init_distribute(ParticleThreadContext *ctx, ParticleSimulationData *sim, int from)
806 {
807         Scene *scene = sim->scene;
808         DerivedMesh *finaldm = sim->psmd->dm_final;
809         Object *ob = sim->ob;
810         ParticleSystem *psys= sim->psys;
811         ParticleData *pa=0, *tpars= 0;
812         ParticleSettings *part;
813         ParticleSeam *seams= 0;
814         KDTree *tree=0;
815         DerivedMesh *dm= NULL;
816         float *jit= NULL;
817         int i, p=0;
818         int cfrom=0;
819         int totelem=0, totpart, *particle_element=0, children=0, totseam=0;
820         int jitlevel= 1, distr;
821         float *element_weight=NULL,*jitter_offset=NULL, *vweight=NULL;
822         float cur, maxweight=0.0, tweight, totweight, inv_totweight, co[3], nor[3], orco[3];
823         
824         if (ELEM(NULL, ob, psys, psys->part))
825                 return 0;
826         
827         part=psys->part;
828         totpart=psys->totpart;
829         if (totpart==0)
830                 return 0;
831         
832         if (!finaldm->deformedOnly && !finaldm->getTessFaceDataArray(finaldm, CD_ORIGINDEX)) {
833                 printf("Can't create particles with the current modifier stack, disable destructive modifiers\n");
834 // XXX          error("Can't paint with the current modifier stack, disable destructive modifiers");
835                 return 0;
836         }
837         
838         /* XXX This distribution code is totally broken in case from == PART_FROM_CHILD, it's always using finaldm
839          *     even if use_modifier_stack is unset... But making things consistent here break all existing edited
840          *     hair systems, so better wait for complete rewrite.
841          */
842
843         psys_thread_context_init(ctx, sim);
844         
845         /* First handle special cases */
846         if (from == PART_FROM_CHILD) {
847                 /* Simple children */
848                 if (part->childtype != PART_CHILD_FACES) {
849                         BLI_srandom(31415926 + psys->seed + psys->child_seed);
850                         distribute_simple_children(scene, ob, finaldm, sim->psmd->dm_deformed, psys);
851                         return 0;
852                 }
853         }
854         else {
855                 /* Grid distribution */
856                 if (part->distr==PART_DISTR_GRID && from != PART_FROM_VERT) {
857                         BLI_srandom(31415926 + psys->seed);
858
859                         if (psys->part->use_modifier_stack) {
860                                 dm = finaldm;
861                         }
862                         else {
863                                 dm = CDDM_from_mesh((Mesh*)ob->data);
864                         }
865                         DM_ensure_tessface(dm);
866
867                         distribute_grid(dm,psys);
868
869                         if (dm != finaldm) {
870                                 dm->release(dm);
871                         }
872
873                         return 0;
874                 }
875         }
876         
877         /* Create trees and original coordinates if needed */
878         if (from == PART_FROM_CHILD) {
879                 distr=PART_DISTR_RAND;
880                 BLI_srandom(31415926 + psys->seed + psys->child_seed);
881                 dm= finaldm;
882
883                 /* BMESH ONLY */
884                 DM_ensure_tessface(dm);
885
886                 children=1;
887
888                 tree=BLI_kdtree_new(totpart);
889
890                 for (p=0,pa=psys->particles; p<totpart; p++,pa++) {
891                         psys_particle_on_dm(dm,part->from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co,nor,0,0,orco,NULL);
892                         BKE_mesh_orco_verts_transform((Mesh*)ob->data, &orco, 1, 1);
893                         BLI_kdtree_insert(tree, p, orco);
894                 }
895
896                 BLI_kdtree_balance(tree);
897
898                 totpart = psys_get_tot_child(scene, psys);
899                 cfrom = from = PART_FROM_FACE;
900         }
901         else {
902                 distr = part->distr;
903                 BLI_srandom(31415926 + psys->seed);
904                 
905                 if (psys->part->use_modifier_stack)
906                         dm = finaldm;
907                 else
908                         dm= CDDM_from_mesh((Mesh*)ob->data);
909
910                 DM_ensure_tessface(dm);
911
912                 /* we need orco for consistent distributions */
913                 if (!CustomData_has_layer(&dm->vertData, CD_ORCO))
914                         DM_add_vert_layer(dm, CD_ORCO, CD_ASSIGN, BKE_mesh_orco_verts_get(ob));
915
916                 if (from == PART_FROM_VERT) {
917                         MVert *mv= dm->getVertDataArray(dm, CD_MVERT);
918                         float (*orcodata)[3] = dm->getVertDataArray(dm, CD_ORCO);
919                         int totvert = dm->getNumVerts(dm);
920
921                         tree=BLI_kdtree_new(totvert);
922
923                         for (p=0; p<totvert; p++) {
924                                 if (orcodata) {
925                                         copy_v3_v3(co,orcodata[p]);
926                                         BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co, 1, 1);
927                                 }
928                                 else
929                                         copy_v3_v3(co,mv[p].co);
930                                 BLI_kdtree_insert(tree, p, co);
931                         }
932
933                         BLI_kdtree_balance(tree);
934                 }
935         }
936
937         /* Get total number of emission elements and allocate needed arrays */
938         totelem = (from == PART_FROM_VERT) ? dm->getNumVerts(dm) : dm->getNumTessFaces(dm);
939
940         if (totelem == 0) {
941                 distribute_invalid(scene, psys, children ? PART_FROM_CHILD : 0);
942
943                 if (G.debug & G_DEBUG)
944                         fprintf(stderr,"Particle distribution error: Nothing to emit from!\n");
945
946                 if (dm != finaldm) dm->release(dm);
947
948                 BLI_kdtree_free(tree);
949
950                 return 0;
951         }
952
953         element_weight  = MEM_callocN(sizeof(float)*totelem, "particle_distribution_weights");
954         particle_element= MEM_callocN(sizeof(int)*totpart, "particle_distribution_indexes");
955         jitter_offset   = MEM_callocN(sizeof(float)*totelem, "particle_distribution_jitoff");
956
957         /* Calculate weights from face areas */
958         if ((part->flag&PART_EDISTR || children) && from != PART_FROM_VERT) {
959                 MVert *v1, *v2, *v3, *v4;
960                 float totarea=0.f, co1[3], co2[3], co3[3], co4[3];
961                 float (*orcodata)[3];
962                 
963                 orcodata= dm->getVertDataArray(dm, CD_ORCO);
964
965                 for (i=0; i<totelem; i++) {
966                         MFace *mf=dm->getTessFaceData(dm,i,CD_MFACE);
967
968                         if (orcodata) {
969                                 copy_v3_v3(co1, orcodata[mf->v1]);
970                                 copy_v3_v3(co2, orcodata[mf->v2]);
971                                 copy_v3_v3(co3, orcodata[mf->v3]);
972                                 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co1, 1, 1);
973                                 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co2, 1, 1);
974                                 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co3, 1, 1);
975                                 if (mf->v4) {
976                                         copy_v3_v3(co4, orcodata[mf->v4]);
977                                         BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co4, 1, 1);
978                                 }
979                         }
980                         else {
981                                 v1= (MVert*)dm->getVertData(dm,mf->v1,CD_MVERT);
982                                 v2= (MVert*)dm->getVertData(dm,mf->v2,CD_MVERT);
983                                 v3= (MVert*)dm->getVertData(dm,mf->v3,CD_MVERT);
984                                 copy_v3_v3(co1, v1->co);
985                                 copy_v3_v3(co2, v2->co);
986                                 copy_v3_v3(co3, v3->co);
987                                 if (mf->v4) {
988                                         v4= (MVert*)dm->getVertData(dm,mf->v4,CD_MVERT);
989                                         copy_v3_v3(co4, v4->co);
990                                 }
991                         }
992
993                         cur = mf->v4 ? area_quad_v3(co1, co2, co3, co4) : area_tri_v3(co1, co2, co3);
994                         
995                         if (cur > maxweight)
996                                 maxweight = cur;
997
998                         element_weight[i] = cur;
999                         totarea += cur;
1000                 }
1001
1002                 for (i=0; i<totelem; i++)
1003                         element_weight[i] /= totarea;
1004
1005                 maxweight /= totarea;
1006         }
1007         else {
1008                 float min=1.0f/(float)(MIN2(totelem,totpart));
1009                 for (i=0; i<totelem; i++)
1010                         element_weight[i]=min;
1011                 maxweight=min;
1012         }
1013
1014         /* Calculate weights from vgroup */
1015         vweight = psys_cache_vgroup(dm,psys,PSYS_VG_DENSITY);
1016
1017         if (vweight) {
1018                 if (from==PART_FROM_VERT) {
1019                         for (i=0;i<totelem; i++)
1020                                 element_weight[i]*=vweight[i];
1021                 }
1022                 else { /* PART_FROM_FACE / PART_FROM_VOLUME */
1023                         for (i=0;i<totelem; i++) {
1024                                 MFace *mf=dm->getTessFaceData(dm,i,CD_MFACE);
1025                                 tweight = vweight[mf->v1] + vweight[mf->v2] + vweight[mf->v3];
1026                                 
1027                                 if (mf->v4) {
1028                                         tweight += vweight[mf->v4];
1029                                         tweight /= 4.0f;
1030                                 }
1031                                 else {
1032                                         tweight /= 3.0f;
1033                                 }
1034
1035                                 element_weight[i]*=tweight;
1036                         }
1037                 }
1038                 MEM_freeN(vweight);
1039         }
1040
1041         /* Calculate total weight of all elements */
1042         int totmapped = 0;
1043         totweight = 0.0f;
1044         for (i = 0; i < totelem; i++) {
1045                 if (element_weight[i] > 0.0f) {
1046                         totmapped++;
1047                         totweight += element_weight[i];
1048                 }
1049         }
1050
1051         if (totmapped == 0) {
1052                 /* We are not allowed to distribute particles anywhere... */
1053                 return 0;
1054         }
1055
1056         inv_totweight = 1.0f / totweight;
1057
1058         /* Calculate cumulative weights.
1059          * We remove all null-weighted elements from element_sum, and create a new mapping
1060          * 'activ'_elem_index -> orig_elem_index.
1061          * This simplifies greatly the filtering of zero-weighted items - and can be much more efficient
1062          * especially in random case (reducing a lot the size of binary-searched array)...
1063          */
1064         float *element_sum = MEM_mallocN(sizeof(*element_sum) * totmapped, __func__);
1065         int *element_map = MEM_mallocN(sizeof(*element_map) * totmapped, __func__);
1066         int i_mapped = 0;
1067
1068         for (i = 0; i < totelem && element_weight[i] == 0.0f; i++);
1069         element_sum[i_mapped] = element_weight[i] * inv_totweight;
1070         element_map[i_mapped] = i;
1071         i_mapped++;
1072         for (i++; i < totelem; i++) {
1073                 if (element_weight[i] > 0.0f) {
1074                         element_sum[i_mapped] = element_sum[i_mapped - 1] + element_weight[i] * inv_totweight;
1075                         /* Skip elements which weight is so small that it does not affect the sum. */
1076                         if (element_sum[i_mapped] > element_sum[i_mapped - 1]) {
1077                                 element_map[i_mapped] = i;
1078                                 i_mapped++;
1079                         }
1080                 }
1081         }
1082         totmapped = i_mapped;
1083
1084         /* Finally assign elements to particles */
1085         if ((part->flag & PART_TRAND) || (part->simplify_flag & PART_SIMPLIFY_ENABLE)) {
1086                 for (p = 0; p < totpart; p++) {
1087                         /* In theory element_sum[totmapped - 1] should be 1.0,
1088                          * but due to float errors this is not necessarily always true, so scale pos accordingly. */
1089                         const float pos = BLI_frand() * element_sum[totmapped - 1];
1090                         const int eidx = distribute_binary_search(element_sum, totmapped, pos);
1091                         particle_element[p] = element_map[eidx];
1092                         BLI_assert(pos <= element_sum[eidx]);
1093                         BLI_assert(eidx ? (pos > element_sum[eidx - 1]) : (pos >= 0.0f));
1094                         jitter_offset[particle_element[p]] = pos;
1095                 }
1096         }
1097         else {
1098                 double step, pos;
1099                 
1100                 step = (totpart < 2) ? 0.5 : 1.0 / (double)totpart;
1101                 /* This is to address tricky issues with vertex-emitting when user tries (and expects) exact 1-1 vert/part
1102                  * distribution (see T47983 and its two example files). It allows us to consider pos as
1103                  * 'midpoint between v and v+1' (or 'p and p+1', depending whether we have more vertices than particles or not),
1104                  * and avoid stumbling over float impression in element_sum.
1105                  * Note: moved face and volume distribution to this as well (instead of starting at zero),
1106                  * for the same reasons, see T52682. */
1107                 pos = (totpart < totmapped) ? 0.5 / (double)totmapped : step * 0.5;  /* We choose the smaller step. */
1108
1109                 for (i = 0, p = 0; p < totpart; p++, pos += step) {
1110                         for ( ; (i < totmapped - 1) && (pos > (double)element_sum[i]); i++);
1111
1112                         particle_element[p] = element_map[i];
1113
1114                         jitter_offset[particle_element[p]] = pos;
1115                 }
1116         }
1117
1118         MEM_freeN(element_sum);
1119         MEM_freeN(element_map);
1120
1121         /* For hair, sort by origindex (allows optimization's in rendering), */
1122         /* however with virtual parents the children need to be in random order. */
1123         if (part->type == PART_HAIR && !(part->childtype==PART_CHILD_FACES && part->parents != 0.0f)) {
1124                 int *orig_index = NULL;
1125
1126                 if (from == PART_FROM_VERT) {
1127                         if (dm->numVertData)
1128                                 orig_index = dm->getVertDataArray(dm, CD_ORIGINDEX);
1129                 }
1130                 else {
1131                         if (dm->numTessFaceData)
1132                                 orig_index = dm->getTessFaceDataArray(dm, CD_ORIGINDEX);
1133                 }
1134
1135                 if (orig_index) {
1136                         BLI_qsort_r(particle_element, totpart, sizeof(int), distribute_compare_orig_index, orig_index);
1137                 }
1138         }
1139
1140         /* Create jittering if needed */
1141         if (distr==PART_DISTR_JIT && ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
1142                 jitlevel= part->userjit;
1143                 
1144                 if (jitlevel == 0) {
1145                         jitlevel= totpart/totelem;
1146                         if (part->flag & PART_EDISTR) jitlevel*= 2;     /* looks better in general, not very scientific */
1147                         if (jitlevel<3) jitlevel= 3;
1148                 }
1149                 
1150                 jit= MEM_callocN((2+ jitlevel*2)*sizeof(float), "jit");
1151
1152                 /* for small amounts of particles we use regular jitter since it looks
1153                  * a bit better, for larger amounts we switch to hammersley sequence 
1154                  * because it is much faster */
1155                 if (jitlevel < 25)
1156                         init_mv_jit(jit, jitlevel, psys->seed, part->jitfac);
1157                 else
1158                         hammersley_create(jit, jitlevel+1, psys->seed, part->jitfac);
1159                 BLI_array_randomize(jit, 2*sizeof(float), jitlevel, psys->seed); /* for custom jit or even distribution */
1160         }
1161
1162         /* Setup things for threaded distribution */
1163         ctx->tree= tree;
1164         ctx->seams= seams;
1165         ctx->totseam= totseam;
1166         ctx->sim.psys= psys;
1167         ctx->index= particle_element;
1168         ctx->jit= jit;
1169         ctx->jitlevel= jitlevel;
1170         ctx->jitoff= jitter_offset;
1171         ctx->weight= element_weight;
1172         ctx->maxweight= maxweight;
1173         ctx->cfrom= cfrom;
1174         ctx->distr= distr;
1175         ctx->dm= dm;
1176         ctx->tpars= tpars;
1177
1178         if (children) {
1179                 totpart= psys_render_simplify_distribution(ctx, totpart);
1180                 alloc_child_particles(psys, totpart);
1181         }
1182
1183         return 1;
1184 }
1185
1186 static void psys_task_init_distribute(ParticleTask *task, ParticleSimulationData *sim)
1187 {
1188         /* init random number generator */
1189         int seed = 31415926 + sim->psys->seed;
1190         
1191         task->rng = BLI_rng_new(seed);
1192 }
1193
1194 static void distribute_particles_on_dm(ParticleSimulationData *sim, int from)
1195 {
1196         TaskScheduler *task_scheduler;
1197         TaskPool *task_pool;
1198         ParticleThreadContext ctx;
1199         ParticleTask *tasks;
1200         DerivedMesh *finaldm = sim->psmd->dm_final;
1201         int i, totpart, numtasks;
1202         
1203         /* create a task pool for distribution tasks */
1204         if (!psys_thread_context_init_distribute(&ctx, sim, from))
1205                 return;
1206         
1207         task_scheduler = BLI_task_scheduler_get();
1208         task_pool = BLI_task_pool_create(task_scheduler, &ctx);
1209         
1210         totpart = (from == PART_FROM_CHILD ? sim->psys->totchild : sim->psys->totpart);
1211         psys_tasks_create(&ctx, 0, totpart, &tasks, &numtasks);
1212         for (i = 0; i < numtasks; ++i) {
1213                 ParticleTask *task = &tasks[i];
1214                 
1215                 psys_task_init_distribute(task, sim);
1216                 if (from == PART_FROM_CHILD)
1217                         BLI_task_pool_push(task_pool, exec_distribute_child, task, false, TASK_PRIORITY_LOW);
1218                 else
1219                         BLI_task_pool_push(task_pool, exec_distribute_parent, task, false, TASK_PRIORITY_LOW);
1220         }
1221         BLI_task_pool_work_and_wait(task_pool);
1222         
1223         BLI_task_pool_free(task_pool);
1224         
1225         psys_calc_dmcache(sim->ob, finaldm, sim->psmd->dm_deformed, sim->psys);
1226         
1227         if (ctx.dm != finaldm)
1228                 ctx.dm->release(ctx.dm);
1229         
1230         psys_tasks_free(tasks, numtasks);
1231         
1232         psys_thread_context_free(&ctx);
1233 }
1234
1235 /* ready for future use, to emit particles without geometry */
1236 static void distribute_particles_on_shape(ParticleSimulationData *sim, int UNUSED(from))
1237 {
1238         distribute_invalid(sim->scene, sim->psys, 0);
1239
1240         fprintf(stderr,"Shape emission not yet possible!\n");
1241 }
1242
1243 void distribute_particles(ParticleSimulationData *sim, int from)
1244 {
1245         PARTICLE_PSMD;
1246         int distr_error=0;
1247
1248         if (psmd) {
1249                 if (psmd->dm_final)
1250                         distribute_particles_on_dm(sim, from);
1251                 else
1252                         distr_error=1;
1253         }
1254         else
1255                 distribute_particles_on_shape(sim, from);
1256
1257         if (distr_error) {
1258                 distribute_invalid(sim->scene, sim->psys, from);
1259
1260                 fprintf(stderr,"Particle distribution error!\n");
1261         }
1262 }
1263
1264 /* ======== Simplify ======== */
1265
1266 static float psys_render_viewport_falloff(double rate, float dist, float width)
1267 {
1268         return pow(rate, dist / width);
1269 }
1270
1271 static float psys_render_projected_area(ParticleSystem *psys, const float center[3], float area, double vprate, float *viewport)
1272 {
1273         ParticleRenderData *data = psys->renderdata;
1274         float co[4], view[3], ortho1[3], ortho2[3], w, dx, dy, radius;
1275         
1276         /* transform to view space */
1277         copy_v3_v3(co, center);
1278         co[3] = 1.0f;
1279         mul_m4_v4(data->viewmat, co);
1280         
1281         /* compute two vectors orthogonal to view vector */
1282         normalize_v3_v3(view, co);
1283         ortho_basis_v3v3_v3(ortho1, ortho2, view);
1284
1285         /* compute on screen minification */
1286         w = co[2] * data->winmat[2][3] + data->winmat[3][3];
1287         dx = data->winx * ortho2[0] * data->winmat[0][0];
1288         dy = data->winy * ortho2[1] * data->winmat[1][1];
1289         w = sqrtf(dx * dx + dy * dy) / w;
1290
1291         /* w squared because we are working with area */
1292         area = area * w * w;
1293
1294         /* viewport of the screen test */
1295
1296         /* project point on screen */
1297         mul_m4_v4(data->winmat, co);
1298         if (co[3] != 0.0f) {
1299                 co[0] = 0.5f * data->winx * (1.0f + co[0] / co[3]);
1300                 co[1] = 0.5f * data->winy * (1.0f + co[1] / co[3]);
1301         }
1302
1303         /* screen space radius */
1304         radius = sqrtf(area / (float)M_PI);
1305
1306         /* make smaller using fallof once over screen edge */
1307         *viewport = 1.0f;
1308
1309         if (co[0] + radius < 0.0f)
1310                 *viewport *= psys_render_viewport_falloff(vprate, -(co[0] + radius), data->winx);
1311         else if (co[0] - radius > data->winx)
1312                 *viewport *= psys_render_viewport_falloff(vprate, (co[0] - radius) - data->winx, data->winx);
1313
1314         if (co[1] + radius < 0.0f)
1315                 *viewport *= psys_render_viewport_falloff(vprate, -(co[1] + radius), data->winy);
1316         else if (co[1] - radius > data->winy)
1317                 *viewport *= psys_render_viewport_falloff(vprate, (co[1] - radius) - data->winy, data->winy);
1318         
1319         return area;
1320 }
1321
1322 /* BMESH_TODO, for orig face data, we need to use MPoly */
1323 static int psys_render_simplify_distribution(ParticleThreadContext *ctx, int tot)
1324 {
1325         DerivedMesh *dm = ctx->dm;
1326         Mesh *me = (Mesh *)(ctx->sim.ob->data);
1327         MFace *mf, *mface;
1328         MVert *mvert;
1329         ParticleRenderData *data;
1330         ParticleRenderElem *elems, *elem;
1331         ParticleSettings *part = ctx->sim.psys->part;
1332         float *facearea, (*facecenter)[3], size[3], fac, powrate, scaleclamp;
1333         float co1[3], co2[3], co3[3], co4[3], lambda, arearatio, t, area, viewport;
1334         double vprate;
1335         int *facetotvert;
1336         int a, b, totorigface, totface, newtot, skipped;
1337
1338         /* double lookup */
1339         const int *index_mf_to_mpoly;
1340         const int *index_mp_to_orig;
1341
1342         if (part->ren_as != PART_DRAW_PATH || !(part->draw & PART_DRAW_REN_STRAND))
1343                 return tot;
1344         if (!ctx->sim.psys->renderdata)
1345                 return tot;
1346
1347         data = ctx->sim.psys->renderdata;
1348         if (data->timeoffset)
1349                 return 0;
1350         if (!(part->simplify_flag & PART_SIMPLIFY_ENABLE))
1351                 return tot;
1352
1353         mvert = dm->getVertArray(dm);
1354         mface = dm->getTessFaceArray(dm);
1355         totface = dm->getNumTessFaces(dm);
1356         totorigface = me->totpoly;
1357
1358         if (totface == 0 || totorigface == 0)
1359                 return tot;
1360
1361         index_mf_to_mpoly = dm->getTessFaceDataArray(dm, CD_ORIGINDEX);
1362         index_mp_to_orig  = dm->getPolyDataArray(dm, CD_ORIGINDEX);
1363         if (index_mf_to_mpoly == NULL) {
1364                 index_mp_to_orig = NULL;
1365         }
1366
1367         facearea = MEM_callocN(sizeof(float) * totorigface, "SimplifyFaceArea");
1368         facecenter = MEM_callocN(sizeof(float[3]) * totorigface, "SimplifyFaceCenter");
1369         facetotvert = MEM_callocN(sizeof(int) * totorigface, "SimplifyFaceArea");
1370         elems = MEM_callocN(sizeof(ParticleRenderElem) * totorigface, "SimplifyFaceElem");
1371
1372         if (data->elems)
1373                 MEM_freeN(data->elems);
1374
1375         data->do_simplify = true;
1376         data->elems = elems;
1377         data->index_mf_to_mpoly = index_mf_to_mpoly;
1378         data->index_mp_to_orig  = index_mp_to_orig;
1379
1380         /* compute number of children per original face */
1381         for (a = 0; a < tot; a++) {
1382                 b = (index_mf_to_mpoly) ? DM_origindex_mface_mpoly(index_mf_to_mpoly, index_mp_to_orig, ctx->index[a]) : ctx->index[a];
1383                 if (b != ORIGINDEX_NONE) {
1384                         elems[b].totchild++;
1385                 }
1386         }
1387
1388         /* compute areas and centers of original faces */
1389         for (mf = mface, a = 0; a < totface; a++, mf++) {
1390                 b = (index_mf_to_mpoly) ? DM_origindex_mface_mpoly(index_mf_to_mpoly, index_mp_to_orig, a) : a;
1391
1392                 if (b != ORIGINDEX_NONE) {
1393                         copy_v3_v3(co1, mvert[mf->v1].co);
1394                         copy_v3_v3(co2, mvert[mf->v2].co);
1395                         copy_v3_v3(co3, mvert[mf->v3].co);
1396
1397                         add_v3_v3(facecenter[b], co1);
1398                         add_v3_v3(facecenter[b], co2);
1399                         add_v3_v3(facecenter[b], co3);
1400
1401                         if (mf->v4) {
1402                                 copy_v3_v3(co4, mvert[mf->v4].co);
1403                                 add_v3_v3(facecenter[b], co4);
1404                                 facearea[b] += area_quad_v3(co1, co2, co3, co4);
1405                                 facetotvert[b] += 4;
1406                         }
1407                         else {
1408                                 facearea[b] += area_tri_v3(co1, co2, co3);
1409                                 facetotvert[b] += 3;
1410                         }
1411                 }
1412         }
1413
1414         for (a = 0; a < totorigface; a++)
1415                 if (facetotvert[a] > 0)
1416                         mul_v3_fl(facecenter[a], 1.0f / facetotvert[a]);
1417
1418         /* for conversion from BU area / pixel area to reference screen size */
1419         BKE_mesh_texspace_get(me, 0, 0, size);
1420         fac = ((size[0] + size[1] + size[2]) / 3.0f) / part->simplify_refsize;
1421         fac = fac * fac;
1422
1423         powrate = log(0.5f) / log(part->simplify_rate * 0.5f);
1424         if (part->simplify_flag & PART_SIMPLIFY_VIEWPORT)
1425                 vprate = pow(1.0f - part->simplify_viewport, 5.0);
1426         else
1427                 vprate = 1.0;
1428
1429         /* set simplification parameters per original face */
1430         for (a = 0, elem = elems; a < totorigface; a++, elem++) {
1431                 area = psys_render_projected_area(ctx->sim.psys, facecenter[a], facearea[a], vprate, &viewport);
1432                 arearatio = fac * area / facearea[a];
1433
1434                 if ((arearatio < 1.0f || viewport < 1.0f) && elem->totchild) {
1435                         /* lambda is percentage of elements to keep */
1436                         lambda = (arearatio < 1.0f) ? powf(arearatio, powrate) : 1.0f;
1437                         lambda *= viewport;
1438
1439                         lambda = MAX2(lambda, 1.0f / elem->totchild);
1440
1441                         /* compute transition region */
1442                         t = part->simplify_transition;
1443                         elem->t = (lambda - t < 0.0f) ? lambda : (lambda + t > 1.0f) ? 1.0f - lambda : t;
1444                         elem->reduce = 1;
1445
1446                         /* scale at end and beginning of the transition region */
1447                         elem->scalemax = (lambda + t < 1.0f) ? 1.0f / lambda : 1.0f / (1.0f - elem->t * elem->t / t);
1448                         elem->scalemin = (lambda + t < 1.0f) ? 0.0f : elem->scalemax * (1.0f - elem->t / t);
1449
1450                         elem->scalemin = sqrtf(elem->scalemin);
1451                         elem->scalemax = sqrtf(elem->scalemax);
1452
1453                         /* clamp scaling */
1454                         scaleclamp = (float)min_ii(elem->totchild, 10);
1455                         elem->scalemin = MIN2(scaleclamp, elem->scalemin);
1456                         elem->scalemax = MIN2(scaleclamp, elem->scalemax);
1457
1458                         /* extend lambda to include transition */
1459                         lambda = lambda + elem->t;
1460                         if (lambda > 1.0f)
1461                                 lambda = 1.0f;
1462                 }
1463                 else {
1464                         lambda = arearatio;
1465
1466                         elem->scalemax = 1.0f; //sqrt(lambda);
1467                         elem->scalemin = 1.0f; //sqrt(lambda);
1468                         elem->reduce = 0;
1469                 }
1470
1471                 elem->lambda = lambda;
1472                 elem->scalemin = sqrtf(elem->scalemin);
1473                 elem->scalemax = sqrtf(elem->scalemax);
1474                 elem->curchild = 0;
1475         }
1476
1477         MEM_freeN(facearea);
1478         MEM_freeN(facecenter);
1479         MEM_freeN(facetotvert);
1480
1481         /* move indices and set random number skipping */
1482         ctx->skip = MEM_callocN(sizeof(int) * tot, "SimplificationSkip");
1483
1484         skipped = 0;
1485         for (a = 0, newtot = 0; a < tot; a++) {
1486                 b = (index_mf_to_mpoly) ? DM_origindex_mface_mpoly(index_mf_to_mpoly, index_mp_to_orig, ctx->index[a]) : ctx->index[a];
1487
1488                 if (b != ORIGINDEX_NONE) {
1489                         if (elems[b].curchild++ < ceil(elems[b].lambda * elems[b].totchild)) {
1490                                 ctx->index[newtot] = ctx->index[a];
1491                                 ctx->skip[newtot] = skipped;
1492                                 skipped = 0;
1493                                 newtot++;
1494                         }
1495                         else skipped++;
1496                 }
1497                 else skipped++;
1498         }
1499
1500         for (a = 0, elem = elems; a < totorigface; a++, elem++)
1501                 elem->curchild = 0;
1502
1503         return newtot;
1504 }