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