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