Texture stack influences are now all separate values, and negative
[blender-staging.git] / source / blender / blenkernel / intern / particle.c
1 /* particle.c
2  *
3  *
4  * $Id$
5  *
6  * ***** BEGIN GPL LICENSE BLOCK *****
7  *
8  * This program is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License
10  * as published by the Free Software Foundation; either version 2
11  * of the License, or (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software Foundation,
20  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
21  *
22  * The Original Code is Copyright (C) 2007 by Janne Karhu.
23  * All rights reserved.
24  *
25  * The Original Code is: all of this file.
26  *
27  * Contributor(s): none yet.
28  *
29  * ***** END GPL LICENSE BLOCK *****
30  */
31
32 #include <stdlib.h>
33 #include <math.h>
34 #include <string.h>
35
36 #include "MEM_guardedalloc.h"
37
38 #include "DNA_scene_types.h"
39 #include "DNA_boid_types.h"
40 #include "DNA_particle_types.h"
41 #include "DNA_mesh_types.h"
42 #include "DNA_meshdata_types.h"
43 #include "DNA_modifier_types.h"
44 #include "DNA_object_force.h"
45 #include "DNA_texture_types.h"
46 #include "DNA_material_types.h"
47 #include "DNA_object_types.h"
48 #include "DNA_curve_types.h"
49 #include "DNA_key_types.h"
50 #include "DNA_ipo_types.h"      // XXX old animation system stuff to remove!
51
52 #include "BLI_arithb.h"
53 #include "BLI_blenlib.h"
54 #include "BLI_dynstr.h"
55 #include "BLI_kdtree.h"
56 #include "BLI_listbase.h"
57 #include "BLI_rand.h"
58 #include "BLI_threads.h"
59
60 #include "BKE_anim.h"
61
62 #include "BKE_boids.h"
63 #include "BKE_cloth.h"
64 #include "BKE_effect.h"
65 #include "BKE_global.h"
66 #include "BKE_main.h"
67 #include "BKE_lattice.h"
68 #include "BKE_utildefines.h"
69 #include "BKE_displist.h"
70 #include "BKE_particle.h"
71 #include "BKE_DerivedMesh.h"
72 #include "BKE_object.h"
73 #include "BKE_cloth.h"
74 #include "BKE_material.h"
75 #include "BKE_key.h"
76 #include "BKE_library.h"
77 #include "BKE_depsgraph.h"
78 #include "BKE_modifier.h"
79 #include "BKE_mesh.h"
80 #include "BKE_cdderivedmesh.h"
81 #include "BKE_pointcache.h"
82
83 #include "RE_render_ext.h"
84
85 static void get_cpa_texture(DerivedMesh *dm, Material *ma, int face_index,
86                                 float *fuv, float *orco, ParticleTexture *ptex, int event);
87 static void get_child_modifier_parameters(ParticleSettings *part, ParticleThreadContext *ctx,
88                                 ChildParticle *cpa, short cpa_from, int cpa_num, float *cpa_fuv, float *orco, ParticleTexture *ptex);
89 static void do_child_modifiers(ParticleSimulationData *sim,
90                                 ParticleTexture *ptex, ParticleKey *par, float *par_rot, ChildParticle *cpa,
91                                 float *orco, float mat[4][4], ParticleKey *state, float t);
92
93 /* few helpers for countall etc. */
94 int count_particles(ParticleSystem *psys){
95         ParticleSettings *part=psys->part;
96         PARTICLE_P;
97         int tot=0;
98
99         LOOP_SHOWN_PARTICLES {
100                 if(pa->alive == PARS_UNBORN && (part->flag & PART_UNBORN)==0);
101                 else if(pa->alive == PARS_DEAD && (part->flag & PART_DIED)==0);
102                 else tot++;
103         }
104         return tot;
105 }
106 int count_particles_mod(ParticleSystem *psys, int totgr, int cur){
107         ParticleSettings *part=psys->part;
108         PARTICLE_P;
109         int tot=0;
110
111         LOOP_SHOWN_PARTICLES {
112                 if(pa->alive == PARS_UNBORN && (part->flag & PART_UNBORN)==0);
113                 else if(pa->alive == PARS_DEAD && (part->flag & PART_DIED)==0);
114                 else if(p%totgr==cur) tot++;
115         }
116         return tot;
117 }
118 /* we allocate path cache memory in chunks instead of a big continguous
119  * chunk, windows' memory allocater fails to find big blocks of memory often */
120
121 #define PATH_CACHE_BUF_SIZE 1024
122
123 static ParticleCacheKey **psys_alloc_path_cache_buffers(ListBase *bufs, int tot, int steps)
124 {
125         LinkData *buf;
126         ParticleCacheKey **cache;
127         int i, totkey, totbufkey;
128
129         tot= MAX2(tot, 1);
130         totkey = 0;
131         cache = MEM_callocN(tot*sizeof(void*), "PathCacheArray");
132
133         while(totkey < tot) {
134                 totbufkey= MIN2(tot-totkey, PATH_CACHE_BUF_SIZE);
135                 buf= MEM_callocN(sizeof(LinkData), "PathCacheLinkData");
136                 buf->data= MEM_callocN(sizeof(ParticleCacheKey)*totbufkey*steps, "ParticleCacheKey");
137
138                 for(i=0; i<totbufkey; i++)
139                         cache[totkey+i] = ((ParticleCacheKey*)buf->data) + i*steps;
140
141                 totkey += totbufkey;
142                 BLI_addtail(bufs, buf);
143         }
144
145         return cache;
146 }
147
148 static void psys_free_path_cache_buffers(ParticleCacheKey **cache, ListBase *bufs)
149 {
150         LinkData *buf;
151
152         if(cache)
153                 MEM_freeN(cache);
154
155         for(buf= bufs->first; buf; buf=buf->next)
156                 MEM_freeN(buf->data);
157         BLI_freelistN(bufs);
158 }
159
160 /************************************************/
161 /*                      Getting stuff                                           */
162 /************************************************/
163 /* get object's active particle system safely */
164 ParticleSystem *psys_get_current(Object *ob)
165 {
166         ParticleSystem *psys;
167         if(ob==0) return 0;
168
169         for(psys=ob->particlesystem.first; psys; psys=psys->next){
170                 if(psys->flag & PSYS_CURRENT)
171                         return psys;
172         }
173         
174         return 0;
175 }
176 short psys_get_current_num(Object *ob)
177 {
178         ParticleSystem *psys;
179         short i;
180
181         if(ob==0) return 0;
182
183         for(psys=ob->particlesystem.first, i=0; psys; psys=psys->next, i++)
184                 if(psys->flag & PSYS_CURRENT)
185                         return i;
186         
187         return i;
188 }
189 void psys_set_current_num(Object *ob, int index)
190 {
191         ParticleSystem *psys;
192         short i;
193
194         if(ob==0) return;
195
196         for(psys=ob->particlesystem.first, i=0; psys; psys=psys->next, i++) {
197                 if(i == index)
198                         psys->flag |= PSYS_CURRENT;
199                 else
200                         psys->flag &= ~PSYS_CURRENT;
201         }
202 }
203 Object *psys_find_object(Scene *scene, ParticleSystem *psys)
204 {
205         Base *base = scene->base.first;
206         ParticleSystem *tpsys;
207
208         for(base = scene->base.first; base; base = base->next) {
209                 for(tpsys = base->object->particlesystem.first; psys; psys=psys->next) {
210                         if(tpsys == psys)
211                                 return base->object;
212                 }
213         }
214
215         return NULL;
216 }
217 Object *psys_get_lattice(ParticleSimulationData *sim)
218 {
219         Object *lattice=0;
220         
221         if(psys_in_edit_mode(sim->scene, sim->psys)==0){
222
223                 ModifierData *md = (ModifierData*)psys_get_modifier(sim->ob, sim->psys);
224
225                 for(; md; md=md->next){
226                         if(md->type==eModifierType_Lattice){
227                                 LatticeModifierData *lmd = (LatticeModifierData *)md;
228                                 lattice=lmd->object;
229                                 break;
230                         }
231                 }
232                 if(lattice)
233                         init_latt_deform(lattice,0);
234         }
235
236         return lattice;
237 }
238 void psys_disable_all(Object *ob)
239 {
240         ParticleSystem *psys=ob->particlesystem.first;
241
242         for(; psys; psys=psys->next)
243                 psys->flag |= PSYS_DISABLED;
244 }
245 void psys_enable_all(Object *ob)
246 {
247         ParticleSystem *psys=ob->particlesystem.first;
248
249         for(; psys; psys=psys->next)
250                 psys->flag &= ~PSYS_DISABLED;
251 }
252 int psys_in_edit_mode(Scene *scene, ParticleSystem *psys)
253 {
254         return (scene->basact && (scene->basact->object->mode & OB_MODE_PARTICLE_EDIT) && psys==psys_get_current((scene->basact)->object) && (psys->edit || psys->pointcache->edit));
255 }
256 static void psys_create_frand(ParticleSystem *psys)
257 {
258         int i;
259         float *rand = psys->frand = MEM_callocN(PSYS_FRAND_COUNT * sizeof(float), "particle randoms");
260
261         BLI_srandom(psys->seed);
262
263         for(i=0; i<1024; i++, rand++)
264                 *rand = BLI_frand();
265 }
266 int psys_check_enabled(Object *ob, ParticleSystem *psys)
267 {
268         ParticleSystemModifierData *psmd;
269         Mesh *me;
270
271         if(psys->flag & PSYS_DISABLED || psys->flag & PSYS_DELETE || !psys->part)
272                 return 0;
273
274         if(ob->type == OB_MESH) {
275                 me= (Mesh*)ob->data;
276                 if(me->mr && me->mr->current != 1)
277                         return 0;
278         }
279
280         psmd= psys_get_modifier(ob, psys);
281         if(psys->renderdata) {
282                 if(!(psmd->modifier.mode & eModifierMode_Render))
283                         return 0;
284         }
285         else if(!(psmd->modifier.mode & eModifierMode_Realtime))
286                 return 0;
287
288         /* perhaps not the perfect place, but we have to be sure the rands are there before usage */
289         if(!psys->frand)
290                 psys_create_frand(psys);
291         else if(psys->recalc & PSYS_RECALC_RESET) {
292                 MEM_freeN(psys->frand);
293                 psys_create_frand(psys);
294         }
295         
296         return 1;
297 }
298
299 /************************************************/
300 /*                      Freeing stuff                                           */
301 /************************************************/
302 void psys_free_settings(ParticleSettings *part)
303 {
304         free_partdeflect(part->pd);
305         free_partdeflect(part->pd2);
306
307         if(part->effector_weights)
308                 MEM_freeN(part->effector_weights);
309
310         boid_free_settings(part->boids);
311 }
312
313 void free_hair(Object *ob, ParticleSystem *psys, int dynamics)
314 {
315         PARTICLE_P;
316
317         if(psys->part->type != PART_HAIR)
318                 return;
319
320         LOOP_PARTICLES {
321                 if(pa->hair)
322                         MEM_freeN(pa->hair);
323                 pa->hair = NULL;
324                 pa->totkey = 0;
325         }
326
327         psys->flag &= ~PSYS_HAIR_DONE;
328
329         if(psys->clmd) {
330                 if(dynamics) {
331                         BKE_ptcache_free_list(&psys->ptcaches);
332                         psys->clmd->point_cache = psys->pointcache = NULL;
333                         psys->clmd->ptcaches.first = psys->clmd->ptcaches.first = NULL;
334
335                         modifier_free((ModifierData*)psys->clmd);
336                         
337                         psys->clmd = NULL;
338                 }
339                 else {
340                         cloth_free_modifier(ob, psys->clmd);
341                 }
342         }
343
344         if(psys->hair_in_dm)
345                 psys->hair_in_dm->release(psys->hair_in_dm);
346         psys->hair_in_dm = NULL;
347
348         if(psys->hair_out_dm)
349                 psys->hair_out_dm->release(psys->hair_out_dm);
350         psys->hair_out_dm = NULL;
351 }
352 void free_keyed_keys(ParticleSystem *psys)
353 {
354         PARTICLE_P;
355
356         if(psys->part->type == PART_HAIR)
357                 return;
358
359         if(psys->particles && psys->particles->keys) {
360                 MEM_freeN(psys->particles->keys);
361
362                 LOOP_PARTICLES {
363                         if(pa->keys) {
364                                 pa->keys= NULL;
365                                 pa->totkey= 0;
366                         }
367                 }
368         }
369 }
370 static void free_child_path_cache(ParticleSystem *psys)
371 {
372         psys_free_path_cache_buffers(psys->childcache, &psys->childcachebufs);
373         psys->childcache = NULL;
374         psys->totchildcache = 0;
375 }
376 void psys_free_path_cache(ParticleSystem *psys, PTCacheEdit *edit)
377 {
378         if(edit) {
379                 psys_free_path_cache_buffers(edit->pathcache, &edit->pathcachebufs);
380                 edit->pathcache= NULL;
381                 edit->totcached= 0;
382         }
383         if(psys) {
384                 psys_free_path_cache_buffers(psys->pathcache, &psys->pathcachebufs);
385                 psys->pathcache= NULL;
386                 psys->totcached= 0;
387
388                 free_child_path_cache(psys);
389         }
390 }
391 void psys_free_children(ParticleSystem *psys)
392 {
393         if(psys->child) {
394                 MEM_freeN(psys->child);
395                 psys->child=0;
396                 psys->totchild=0;
397         }
398
399         free_child_path_cache(psys);
400 }
401 void psys_free_particles(ParticleSystem *psys)
402 {
403         PARTICLE_P;
404
405         if(psys->particles) {
406                 if(psys->part->type==PART_HAIR) {
407                         LOOP_PARTICLES {
408                                 if(pa->hair)
409                                         MEM_freeN(pa->hair);
410                         }
411                 }
412                 
413                 if(psys->particles->keys)
414                         MEM_freeN(psys->particles->keys);
415                 
416                 if(psys->particles->boid)
417                         MEM_freeN(psys->particles->boid);
418
419                 MEM_freeN(psys->particles);
420                 psys->particles= NULL;
421                 psys->totpart= 0;
422         }
423 }
424 void psys_free_pdd(ParticleSystem *psys)
425 {
426         if(psys->pdd) {
427                 if(psys->pdd->cdata)
428                         MEM_freeN(psys->pdd->cdata);
429                 psys->pdd->cdata = NULL;
430
431                 if(psys->pdd->vdata)
432                         MEM_freeN(psys->pdd->vdata);
433                 psys->pdd->vdata = NULL;
434
435                 if(psys->pdd->ndata)
436                         MEM_freeN(psys->pdd->ndata);
437                 psys->pdd->ndata = NULL;
438
439                 if(psys->pdd->vedata)
440                         MEM_freeN(psys->pdd->vedata);
441                 psys->pdd->vedata = NULL;
442         }
443 }
444 /* free everything */
445 void psys_free(Object *ob, ParticleSystem * psys)
446 {       
447         if(psys){
448                 int nr = 0;
449                 ParticleSystem * tpsys;
450                 
451                 psys_free_path_cache(psys, NULL);
452
453                 free_hair(ob, psys, 1);
454
455                 psys_free_particles(psys);
456
457                 if(psys->edit && psys->free_edit)
458                         psys->free_edit(psys->edit);
459
460                 if(psys->child){
461                         MEM_freeN(psys->child);
462                         psys->child = 0;
463                         psys->totchild = 0;
464                 }
465                 
466                 // check if we are last non-visible particle system
467                 for(tpsys=ob->particlesystem.first; tpsys; tpsys=tpsys->next){
468                         if(tpsys->part)
469                         {
470                                 if(ELEM(tpsys->part->ren_as,PART_DRAW_OB,PART_DRAW_GR))
471                                 {
472                                         nr++;
473                                         break;
474                                 }
475                         }
476                 }
477                 // clear do-not-draw-flag
478                 if(!nr)
479                         ob->transflag &= ~OB_DUPLIPARTS;
480
481                 if(psys->part){
482                         psys->part->id.us--;            
483                         psys->part=0;
484                 }
485
486                 BKE_ptcache_free_list(&psys->ptcaches);
487                 psys->pointcache = NULL;
488                 
489                 BLI_freelistN(&psys->targets);
490
491                 BLI_kdtree_free(psys->tree);
492
493                 pdEndEffectors(&psys->effectors);
494
495                 if(psys->frand)
496                         MEM_freeN(psys->frand);
497
498                 if(psys->pdd) {
499                         psys_free_pdd(psys);
500                         MEM_freeN(psys->pdd);
501                 }
502
503                 MEM_freeN(psys);
504         }
505 }
506
507 /************************************************/
508 /*                      Rendering                                                       */
509 /************************************************/
510 /* these functions move away particle data and bring it back after
511  * rendering, to make different render settings possible without
512  * removing the previous data. this should be solved properly once */
513
514 typedef struct ParticleRenderElem {
515         int curchild, totchild, reduce;
516         float lambda, t, scalemin, scalemax;
517 } ParticleRenderElem;
518
519 typedef struct ParticleRenderData {
520         ChildParticle *child;
521         ParticleCacheKey **pathcache;
522         ParticleCacheKey **childcache;
523         ListBase pathcachebufs, childcachebufs;
524         int totchild, totcached, totchildcache;
525         DerivedMesh *dm;
526         int totdmvert, totdmedge, totdmface;
527
528         float mat[4][4];
529         float viewmat[4][4], winmat[4][4];
530         int winx, winy;
531
532         int dosimplify;
533         int timeoffset;
534         ParticleRenderElem *elems;
535         int *origindex;
536 } ParticleRenderData;
537
538 static float psys_render_viewport_falloff(double rate, float dist, float width)
539 {
540         return pow(rate, dist/width);
541 }
542
543 static float psys_render_projected_area(ParticleSystem *psys, float *center, float area, double vprate, float *viewport)
544 {
545         ParticleRenderData *data= psys->renderdata;
546         float co[4], view[3], ortho1[3], ortho2[3], w, dx, dy, radius;
547         
548         /* transform to view space */
549         VECCOPY(co, center);
550         co[3]= 1.0f;
551         Mat4MulVec4fl(data->viewmat, co);
552         
553         /* compute two vectors orthogonal to view vector */
554         VECCOPY(view, co);
555         Normalize(view);
556         VecOrthoBasisf(view, ortho1, ortho2);
557
558         /* compute on screen minification */
559         w= co[2]*data->winmat[2][3] + data->winmat[3][3];
560         dx= data->winx*ortho2[0]*data->winmat[0][0];
561         dy= data->winy*ortho2[1]*data->winmat[1][1];
562         w= sqrt(dx*dx + dy*dy)/w;
563
564         /* w squared because we are working with area */
565         area= area*w*w;
566
567         /* viewport of the screen test */
568
569         /* project point on screen */
570         Mat4MulVec4fl(data->winmat, co);
571         if(co[3] != 0.0f) {
572                 co[0]= 0.5f*data->winx*(1.0f + co[0]/co[3]);
573                 co[1]= 0.5f*data->winy*(1.0f + co[1]/co[3]);
574         }
575
576         /* screen space radius */
577         radius= sqrt(area/M_PI);
578
579         /* make smaller using fallof once over screen edge */
580         *viewport= 1.0f;
581
582         if(co[0]+radius < 0.0f)
583                 *viewport *= psys_render_viewport_falloff(vprate, -(co[0]+radius), data->winx);
584         else if(co[0]-radius > data->winx)
585                 *viewport *= psys_render_viewport_falloff(vprate, (co[0]-radius) - data->winx, data->winx);
586
587         if(co[1]+radius < 0.0f)
588                 *viewport *= psys_render_viewport_falloff(vprate, -(co[1]+radius), data->winy);
589         else if(co[1]-radius > data->winy)
590                 *viewport *= psys_render_viewport_falloff(vprate, (co[1]-radius) - data->winy, data->winy);
591         
592         return area;
593 }
594
595 void psys_render_set(Object *ob, ParticleSystem *psys, float viewmat[][4], float winmat[][4], int winx, int winy, int timeoffset)
596 {
597         ParticleRenderData*data;
598         ParticleSystemModifierData *psmd= psys_get_modifier(ob, psys);
599
600         if(!G.rendering)
601                 return;
602         if(psys->renderdata)
603                 return;
604
605         data= MEM_callocN(sizeof(ParticleRenderData), "ParticleRenderData");
606
607         data->child= psys->child;
608         data->totchild= psys->totchild;
609         data->pathcache= psys->pathcache;
610         data->pathcachebufs.first = psys->pathcachebufs.first;
611         data->pathcachebufs.last = psys->pathcachebufs.last;
612         data->totcached= psys->totcached;
613         data->childcache= psys->childcache;
614         data->childcachebufs.first = psys->childcachebufs.first;
615         data->childcachebufs.last = psys->childcachebufs.last;
616         data->totchildcache= psys->totchildcache;
617
618         if(psmd->dm)
619                 data->dm= CDDM_copy(psmd->dm);
620         data->totdmvert= psmd->totdmvert;
621         data->totdmedge= psmd->totdmedge;
622         data->totdmface= psmd->totdmface;
623
624         psys->child= NULL;
625         psys->pathcache= NULL;
626         psys->childcache= NULL;
627         psys->totchild= psys->totcached= psys->totchildcache= 0;
628         psys->pathcachebufs.first = psys->pathcachebufs.last = NULL;
629         psys->childcachebufs.first = psys->childcachebufs.last = NULL;
630
631         Mat4CpyMat4(data->winmat, winmat);
632         Mat4MulMat4(data->viewmat, ob->obmat, viewmat);
633         Mat4MulMat4(data->mat, data->viewmat, winmat);
634         data->winx= winx;
635         data->winy= winy;
636
637         data->timeoffset= timeoffset;
638
639         psys->renderdata= data;
640 }
641
642 void psys_render_restore(Object *ob, ParticleSystem *psys)
643 {
644         ParticleRenderData*data;
645         ParticleSystemModifierData *psmd= psys_get_modifier(ob, psys);
646
647         data= psys->renderdata;
648         if(!data)
649                 return;
650         
651         if(data->elems)
652                 MEM_freeN(data->elems);
653
654         if(psmd->dm) {
655                 psmd->dm->needsFree= 1;
656                 psmd->dm->release(psmd->dm);
657         }
658
659         psys_free_path_cache(psys, NULL);
660
661         if(psys->child){
662                 MEM_freeN(psys->child);
663                 psys->child= 0;
664                 psys->totchild= 0;
665         }
666
667         psys->child= data->child;
668         psys->totchild= data->totchild;
669         psys->pathcache= data->pathcache;
670         psys->pathcachebufs.first = data->pathcachebufs.first;
671         psys->pathcachebufs.last = data->pathcachebufs.last;
672         psys->totcached= data->totcached;
673         psys->childcache= data->childcache;
674         psys->childcachebufs.first = data->childcachebufs.first;
675         psys->childcachebufs.last = data->childcachebufs.last;
676         psys->totchildcache= data->totchildcache;
677
678         psmd->dm= data->dm;
679         psmd->totdmvert= data->totdmvert;
680         psmd->totdmedge= data->totdmedge;
681         psmd->totdmface= data->totdmface;
682         psmd->flag &= ~eParticleSystemFlag_psys_updated;
683
684         if(psmd->dm)
685                 psys_calc_dmcache(ob, psmd->dm, psys);
686
687         MEM_freeN(data);
688         psys->renderdata= NULL;
689 }
690
691 int psys_render_simplify_distribution(ParticleThreadContext *ctx, int tot)
692 {
693         DerivedMesh *dm= ctx->dm;
694         Mesh *me= (Mesh*)(ctx->sim.ob->data);
695         MFace *mf, *mface;
696         MVert *mvert;
697         ParticleRenderData *data;
698         ParticleRenderElem *elems, *elem;
699         ParticleSettings *part= ctx->sim.psys->part;
700         float *facearea, (*facecenter)[3], size[3], fac, powrate, scaleclamp;
701         float co1[3], co2[3], co3[3], co4[3], lambda, arearatio, t, area, viewport;
702         double vprate;
703         int *origindex, *facetotvert;
704         int a, b, totorigface, totface, newtot, skipped;
705
706         if(part->ren_as!=PART_DRAW_PATH || !(part->draw & PART_DRAW_REN_STRAND))
707                 return tot;
708         if(!ctx->sim.psys->renderdata)
709                 return tot;
710
711         data= ctx->sim.psys->renderdata;
712         if(data->timeoffset)
713                 return 0;
714         if(!(part->simplify_flag & PART_SIMPLIFY_ENABLE))
715                 return tot;
716
717         mvert= dm->getVertArray(dm);
718         mface= dm->getFaceArray(dm);
719         origindex= dm->getFaceDataArray(dm, CD_ORIGINDEX);
720         totface= dm->getNumFaces(dm);
721         totorigface= me->totface;
722
723         if(totface == 0 || totorigface == 0 || origindex == NULL)
724                 return tot;
725
726         facearea= MEM_callocN(sizeof(float)*totorigface, "SimplifyFaceArea");
727         facecenter= MEM_callocN(sizeof(float[3])*totorigface, "SimplifyFaceCenter");
728         facetotvert= MEM_callocN(sizeof(int)*totorigface, "SimplifyFaceArea");
729         elems= MEM_callocN(sizeof(ParticleRenderElem)*totorigface, "SimplifyFaceElem");
730
731         if(data->elems)
732                 MEM_freeN(data->elems);
733
734         data->dosimplify= 1;
735         data->elems= elems;
736         data->origindex= origindex;
737
738         /* compute number of children per original face */
739         for(a=0; a<tot; a++) {
740                 b= origindex[ctx->index[a]];
741                 if(b != -1)
742                         elems[b].totchild++;
743         }
744
745         /* compute areas and centers of original faces */
746         for(mf=mface, a=0; a<totface; a++, mf++) {
747                 b= origindex[a];
748
749                 if(b != -1) {
750                         VECCOPY(co1, mvert[mf->v1].co);
751                         VECCOPY(co2, mvert[mf->v2].co);
752                         VECCOPY(co3, mvert[mf->v3].co);
753
754                         VECADD(facecenter[b], facecenter[b], co1);
755                         VECADD(facecenter[b], facecenter[b], co2);
756                         VECADD(facecenter[b], facecenter[b], co3);
757
758                         if(mf->v4) {
759                                 VECCOPY(co4, mvert[mf->v4].co);
760                                 VECADD(facecenter[b], facecenter[b], co4);
761                                 facearea[b] += AreaQ3Dfl(co1, co2, co3, co4);
762                                 facetotvert[b] += 4;
763                         }
764                         else {
765                                 facearea[b] += AreaT3Dfl(co1, co2, co3);
766                                 facetotvert[b] += 3;
767                         }
768                 }
769         }
770
771         for(a=0; a<totorigface; a++)
772                 if(facetotvert[a] > 0)
773                         VecMulf(facecenter[a], 1.0f/facetotvert[a]);
774
775         /* for conversion from BU area / pixel area to reference screen size */
776         mesh_get_texspace(me, 0, 0, size);
777         fac= ((size[0] + size[1] + size[2])/3.0f)/part->simplify_refsize;
778         fac= fac*fac;
779
780         powrate= log(0.5f)/log(part->simplify_rate*0.5f);
781         if(part->simplify_flag & PART_SIMPLIFY_VIEWPORT)
782                 vprate= pow(1.0 - part->simplify_viewport, 5.0);
783         else
784                 vprate= 1.0;
785
786         /* set simplification parameters per original face */
787         for(a=0, elem=elems; a<totorigface; a++, elem++) {
788                 area = psys_render_projected_area(ctx->sim.psys, facecenter[a], facearea[a], vprate, &viewport);
789                 arearatio= fac*area/facearea[a];
790
791                 if((arearatio < 1.0f || viewport < 1.0f) && elem->totchild) {
792                         /* lambda is percentage of elements to keep */
793                         lambda= (arearatio < 1.0f)? pow(arearatio, powrate): 1.0f;
794                         lambda *= viewport;
795
796                         lambda= MAX2(lambda, 1.0f/elem->totchild);
797
798                         /* compute transition region */
799                         t= part->simplify_transition;
800                         elem->t= (lambda-t < 0.0f)? lambda: (lambda+t > 1.0f)? 1.0f-lambda: t;
801                         elem->reduce= 1;
802
803                         /* scale at end and beginning of the transition region */
804                         elem->scalemax= (lambda+t < 1.0f)? 1.0f/lambda: 1.0f/(1.0f - elem->t*elem->t/t);
805                         elem->scalemin= (lambda+t < 1.0f)? 0.0f: elem->scalemax*(1.0f-elem->t/t);
806
807                         elem->scalemin= sqrt(elem->scalemin);
808                         elem->scalemax= sqrt(elem->scalemax);
809
810                         /* clamp scaling */
811                         scaleclamp= MIN2(elem->totchild, 10.0f);
812                         elem->scalemin= MIN2(scaleclamp, elem->scalemin);
813                         elem->scalemax= MIN2(scaleclamp, elem->scalemax);
814
815                         /* extend lambda to include transition */
816                         lambda= lambda + elem->t;
817                         if(lambda > 1.0f)
818                                 lambda= 1.0f;
819                 }
820                 else {
821                         lambda= arearatio;
822
823                         elem->scalemax= 1.0f; //sqrt(lambda);
824                         elem->scalemin= 1.0f; //sqrt(lambda);
825                         elem->reduce= 0;
826                 }
827
828                 elem->lambda= lambda;
829                 elem->scalemin= sqrt(elem->scalemin);
830                 elem->scalemax= sqrt(elem->scalemax);
831                 elem->curchild= 0;
832         }
833
834         MEM_freeN(facearea);
835         MEM_freeN(facecenter);
836         MEM_freeN(facetotvert);
837
838         /* move indices and set random number skipping */
839         ctx->skip= MEM_callocN(sizeof(int)*tot, "SimplificationSkip");
840
841         skipped= 0;
842         for(a=0, newtot=0; a<tot; a++) {
843                 b= origindex[ctx->index[a]];
844                 if(b != -1) {
845                         if(elems[b].curchild++ < ceil(elems[b].lambda*elems[b].totchild)) {
846                                 ctx->index[newtot]= ctx->index[a];
847                                 ctx->skip[newtot]= skipped;
848                                 skipped= 0;
849                                 newtot++;
850                         }
851                         else skipped++;
852                 }
853                 else skipped++;
854         }
855
856         for(a=0, elem=elems; a<totorigface; a++, elem++)
857                 elem->curchild= 0;
858
859         return newtot;
860 }
861
862 int psys_render_simplify_params(ParticleSystem *psys, ChildParticle *cpa, float *params)
863 {
864         ParticleRenderData *data;
865         ParticleRenderElem *elem;
866         float x, w, scale, alpha, lambda, t, scalemin, scalemax;
867         int b;
868
869         if(!(psys->renderdata && (psys->part->simplify_flag & PART_SIMPLIFY_ENABLE)))
870                 return 0;
871         
872         data= psys->renderdata;
873         if(!data->dosimplify)
874                 return 0;
875         
876         b= data->origindex[cpa->num];
877         if(b == -1)
878                 return 0;
879
880         elem= &data->elems[b];
881
882         lambda= elem->lambda;
883         t= elem->t;
884         scalemin= elem->scalemin;
885         scalemax= elem->scalemax;
886
887         if(!elem->reduce) {
888                 scale= scalemin;
889                 alpha= 1.0f;
890         }
891         else {
892                 x= (elem->curchild+0.5f)/elem->totchild;
893                 if(x < lambda-t) {
894                         scale= scalemax;
895                         alpha= 1.0f;
896                 }
897                 else if(x >= lambda+t) {
898                         scale= scalemin;
899                         alpha= 0.0f;
900                 }
901                 else {
902                         w= (lambda+t - x)/(2.0f*t);
903                         scale= scalemin + (scalemax - scalemin)*w;
904                         alpha= w;
905                 }
906         }
907
908         params[0]= scale;
909         params[1]= alpha;
910
911         elem->curchild++;
912
913         return 1;
914 }
915
916 /************************************************/
917 /*                      Interpolation                                           */
918 /************************************************/
919 static float interpolate_particle_value(float v1, float v2, float v3, float v4, float *w, int four)
920 {
921         float value;
922
923         value= w[0]*v1 + w[1]*v2 + w[2]*v3;
924         if(four)
925                 value += w[3]*v4;
926         
927         return value;
928 }
929 static void weighted_particle_vector(float *v1, float *v2, float *v3, float *v4, float *weights, float *vec)
930 {
931         vec[0]= weights[0]*v1[0] + weights[1]*v2[0] + weights[2]*v3[0] + weights[3]*v4[0];
932         vec[1]= weights[0]*v1[1] + weights[1]*v2[1] + weights[2]*v3[1] + weights[3]*v4[1];
933         vec[2]= weights[0]*v1[2] + weights[1]*v2[2] + weights[2]*v3[2] + weights[3]*v4[2];
934 }
935 void psys_interpolate_particle(short type, ParticleKey keys[4], float dt, ParticleKey *result, int velocity)
936 {
937         float t[4];
938
939         if(type<0) {
940                 VecfCubicInterpol(keys[1].co, keys[1].vel, keys[2].co, keys[2].vel, dt, result->co, result->vel);
941         }
942         else {
943                 key_curve_position_weights(dt, t, type);
944
945                 weighted_particle_vector(keys[0].co, keys[1].co, keys[2].co, keys[3].co, t, result->co);
946
947                 if(velocity){
948                         float temp[3];
949
950                         if(dt>0.999f){
951                                 key_curve_position_weights(dt-0.001f, t, type);
952                                 weighted_particle_vector(keys[0].co, keys[1].co, keys[2].co, keys[3].co, t, temp);
953                                 VECSUB(result->vel, result->co, temp);
954                         }
955                         else{
956                                 key_curve_position_weights(dt+0.001f, t, type);
957                                 weighted_particle_vector(keys[0].co, keys[1].co, keys[2].co, keys[3].co, t, temp);
958                                 VECSUB(result->vel, temp, result->co);
959                         }
960                 }
961         }
962 }
963
964
965
966 typedef struct ParticleInterpolationData {
967         HairKey *hkey[2];
968
969         DerivedMesh *dm;
970         MVert *mvert[2];
971
972         int keyed;
973         ParticleKey *kkey[2];
974
975         PointCache *cache;
976
977         PTCacheEditPoint *epoint;
978         PTCacheEditKey *ekey[2];
979
980         float birthtime, dietime;
981         int bspline;
982 } ParticleInterpolationData;
983 /* Assumes pointcache->mem_cache exists, so for disk cached particles call psys_make_temp_pointcache() before use */
984 static void get_pointcache_keys_for_time(Object *ob, PointCache *cache, int index, float t, ParticleKey *key1, ParticleKey *key2)
985 {
986         static PTCacheMem *pm = NULL; /* not thread safe */
987
988         if(index < 0) { /* initialize */
989                 pm = cache->mem_cache.first;
990
991                 if(pm)
992                         pm = pm->next;
993         }
994         else {
995                 if(pm) {
996                         while(pm && pm->next && (float)pm->frame < t)
997                                 pm = pm->next;
998
999                         BKE_ptcache_make_particle_key(key2, pm->index_array ? pm->index_array[index] : index, pm->data, (float)pm->frame);
1000                         BKE_ptcache_make_particle_key(key1, pm->prev->index_array ? pm->prev->index_array[index] : index, pm->prev->data, (float)pm->prev->frame);
1001                 }
1002                 else if(cache->mem_cache.first) {
1003                         PTCacheMem *pm2 = cache->mem_cache.first;
1004                         BKE_ptcache_make_particle_key(key2, pm2->index_array ? pm2->index_array[index] : index, pm2->data, (float)pm2->frame);
1005                         copy_particle_key(key1, key2, 1);
1006                 }
1007         }
1008 }
1009 static void init_particle_interpolation(Object *ob, ParticleSystem *psys, ParticleData *pa, ParticleInterpolationData *pind)
1010 {
1011
1012         if(pind->epoint) {
1013                 PTCacheEditPoint *point = pind->epoint;
1014
1015                 pind->ekey[0] = point->keys;
1016                 pind->ekey[1] = point->totkey > 1 ? point->keys + 1 : NULL;
1017
1018                 pind->birthtime = *(point->keys->time);
1019                 pind->dietime = *((point->keys + point->totkey - 1)->time);
1020         }
1021         else if(pind->keyed) {
1022                 ParticleKey *key = pa->keys;
1023                 pind->kkey[0] = key;
1024                 pind->kkey[1] = pa->totkey > 1 ? key + 1 : NULL;
1025
1026                 pind->birthtime = key->time;
1027                 pind->dietime = (key + pa->totkey - 1)->time;
1028         }
1029         else if(pind->cache) {
1030                 get_pointcache_keys_for_time(ob, pind->cache, -1, 0.0f, NULL, NULL);
1031
1032                 pind->birthtime = pa ? pa->time : pind->cache->startframe;
1033                 pind->dietime = pa ? pa->dietime : pind->cache->endframe;
1034         }
1035         else {
1036                 HairKey *key = pa->hair;
1037                 pind->hkey[0] = key;
1038                 pind->hkey[1] = key + 1;
1039
1040                 pind->birthtime = key->time;
1041                 pind->dietime = (key + pa->totkey - 1)->time;
1042
1043                 if(pind->dm) {
1044                         pind->mvert[0] = CDDM_get_vert(pind->dm, pa->hair_index);
1045                         pind->mvert[1] = pind->mvert[0] + 1;
1046                 }
1047         }
1048 }
1049 static void edit_to_particle(ParticleKey *key, PTCacheEditKey *ekey)
1050 {
1051         VECCOPY(key->co, ekey->co);
1052         if(ekey->vel) {
1053                 VECCOPY(key->vel, ekey->vel);
1054         }
1055         key->time = *(ekey->time);
1056 }
1057 static void hair_to_particle(ParticleKey *key, HairKey *hkey)
1058 {
1059         VECCOPY(key->co, hkey->co);
1060         key->time = hkey->time;
1061 }
1062
1063 static void mvert_to_particle(ParticleKey *key, MVert *mvert, HairKey *hkey)
1064 {
1065         VECCOPY(key->co, mvert->co);
1066         key->time = hkey->time;
1067 }
1068
1069 static void do_particle_interpolation(ParticleSystem *psys, int p, ParticleData *pa, float t, float frs_sec, ParticleInterpolationData *pind, ParticleKey *result)
1070 {
1071         PTCacheEditPoint *point = pind->epoint;
1072         ParticleKey keys[4];
1073         int point_vel = (point && point->keys->vel);
1074         float real_t, dfra, keytime;
1075
1076         /* interpret timing and find keys */
1077         if(point) {
1078                 if(result->time < 0.0f)
1079                         real_t = -result->time;
1080                 else
1081                         real_t = *(pind->ekey[0]->time) + t * (*(pind->ekey[0][point->totkey-1].time) - *(pind->ekey[0]->time));
1082
1083                 while(*(pind->ekey[1]->time) < real_t)
1084                         pind->ekey[1]++;
1085
1086                 pind->ekey[0] = pind->ekey[1] - 1;
1087         }
1088         else if(pind->keyed) {
1089                 /* we have only one key, so let's use that */
1090                 if(pind->kkey[1]==NULL) {
1091                         copy_particle_key(result, pind->kkey[0], 1);
1092                         return;
1093                 }
1094
1095                 if(result->time < 0.0f)
1096                         real_t = -result->time;
1097                 else
1098                         real_t = pind->kkey[0]->time + t * (pind->kkey[0][pa->totkey-1].time - pind->kkey[0]->time);
1099
1100                 if(psys->part->phystype==PART_PHYS_KEYED && psys->flag & PSYS_KEYED_TIMING) {
1101                         ParticleTarget *pt = psys->targets.first;
1102
1103                         pt=pt->next;
1104
1105                         while(pt && pa->time + pt->time < real_t)
1106                                 pt= pt->next;
1107
1108                         if(pt) {
1109                                 pt=pt->prev;
1110
1111                                 if(pa->time + pt->time + pt->duration > real_t)
1112                                         real_t = pa->time + pt->time;
1113                         }
1114                         else
1115                                 real_t = pa->time + ((ParticleTarget*)psys->targets.last)->time;
1116                 }
1117
1118                 CLAMP(real_t, pa->time, pa->dietime);
1119
1120                 while(pind->kkey[1]->time < real_t)
1121                         pind->kkey[1]++;
1122                 
1123                 pind->kkey[0] = pind->kkey[1] - 1;
1124         }
1125         else if(pind->cache) {
1126                 if(result->time < 0.0f) /* flag for time in frames */
1127                         real_t = -result->time;
1128                 else
1129                         real_t = pa->time + t * (pa->dietime - pa->time);
1130         }
1131         else {
1132                 if(result->time < 0.0f)
1133                         real_t = -result->time;
1134                 else
1135                         real_t = pind->hkey[0]->time + t * (pind->hkey[0][pa->totkey-1].time - pind->hkey[0]->time);
1136
1137                 while(pind->hkey[1]->time < real_t) {
1138                         pind->hkey[1]++;
1139                         pind->mvert[1]++;
1140                 }
1141
1142                 pind->hkey[0] = pind->hkey[1] - 1;
1143         }
1144
1145         /* set actual interpolation keys */
1146         if(point) {
1147                 edit_to_particle(keys + 1, pind->ekey[0]);
1148                 edit_to_particle(keys + 2, pind->ekey[1]);
1149         }
1150         else if(pind->dm) {
1151                 pind->mvert[0] = pind->mvert[1] - 1;
1152                 mvert_to_particle(keys + 1, pind->mvert[0], pind->hkey[0]);
1153                 mvert_to_particle(keys + 2, pind->mvert[1], pind->hkey[1]);
1154         }
1155         else if(pind->keyed) {
1156                 memcpy(keys + 1, pind->kkey[0], sizeof(ParticleKey));
1157                 memcpy(keys + 2, pind->kkey[1], sizeof(ParticleKey));
1158         }
1159         else if(pind->cache) {
1160                 get_pointcache_keys_for_time(NULL, pind->cache, p, real_t, keys+1, keys+2);
1161         }
1162         else {
1163                 hair_to_particle(keys + 1, pind->hkey[0]);
1164                 hair_to_particle(keys + 2, pind->hkey[1]);
1165         }
1166
1167         /* set secondary interpolation keys for hair */
1168         if(!pind->keyed && !pind->cache && !point_vel) {
1169                 if(point) {
1170                         if(pind->ekey[0] != point->keys)
1171                                 edit_to_particle(keys, pind->ekey[0] - 1);
1172                         else
1173                                 edit_to_particle(keys, pind->ekey[0]);
1174                 }
1175                 else if(pind->dm) {
1176                         if(pind->hkey[0] != pa->hair)
1177                                 mvert_to_particle(keys, pind->mvert[0] - 1, pind->hkey[0] - 1);
1178                         else
1179                                 mvert_to_particle(keys, pind->mvert[0], pind->hkey[0]);
1180                 }
1181                 else {
1182                         if(pind->hkey[0] != pa->hair)
1183                                 hair_to_particle(keys, pind->hkey[0] - 1);
1184                         else
1185                                 hair_to_particle(keys, pind->hkey[0]);
1186                 }
1187
1188                 if(point) {
1189                         if(pind->ekey[1] != point->keys + point->totkey - 1)
1190                                 edit_to_particle(keys + 3, pind->ekey[1] + 1);
1191                         else
1192                                 edit_to_particle(keys + 3, pind->ekey[1]);
1193                 }
1194                 else if(pind->dm) {
1195                         if(pind->hkey[1] != pa->hair + pa->totkey - 1)
1196                                 mvert_to_particle(keys + 3, pind->mvert[1] + 1, pind->hkey[1] + 1);
1197                         else
1198                                 mvert_to_particle(keys + 3, pind->mvert[1], pind->hkey[1]);
1199                 }
1200                 else {
1201                         if(pind->hkey[1] != pa->hair + pa->totkey - 1)
1202                                 hair_to_particle(keys + 3, pind->hkey[1] + 1);
1203                         else
1204                                 hair_to_particle(keys + 3, pind->hkey[1]);
1205                 }
1206         }
1207
1208         dfra = keys[2].time - keys[1].time;
1209         keytime = (real_t - keys[1].time) / dfra;
1210
1211         /* convert velocity to timestep size */
1212         if(pind->keyed || pind->cache || point_vel){
1213                 VecMulf(keys[1].vel, dfra / frs_sec);
1214                 VecMulf(keys[2].vel, dfra / frs_sec);
1215                 QuatInterpol(result->rot,keys[1].rot,keys[2].rot,keytime);
1216         }
1217
1218         /* now we should have in chronologiacl order k1<=k2<=t<=k3<=k4 with keytime between [0,1]->[k2,k3] (k1 & k4 used for cardinal & bspline interpolation)*/
1219         psys_interpolate_particle((pind->keyed || pind->cache || point_vel) ? -1 /* signal for cubic interpolation */
1220                 : (pind->bspline ? KEY_BSPLINE : KEY_CARDINAL)
1221                 ,keys, keytime, result, 1);
1222
1223         /* the velocity needs to be converted back from cubic interpolation */
1224         if(pind->keyed || pind->cache || point_vel)
1225                 VecMulf(result->vel, frs_sec / dfra);
1226 }
1227 /************************************************/
1228 /*                      Particles on a dm                                       */
1229 /************************************************/
1230 /* interpolate a location on a face based on face coordinates */
1231 void psys_interpolate_face(MVert *mvert, MFace *mface, MTFace *tface, float (*orcodata)[3], float *w, float *vec, float *nor, float *utan, float *vtan, float *orco,float *ornor){
1232         float *v1=0, *v2=0, *v3=0, *v4=0;
1233         float e1[3],e2[3],s1,s2,t1,t2;
1234         float *uv1, *uv2, *uv3, *uv4;
1235         float n1[3], n2[3], n3[3], n4[3];
1236         float tuv[4][2];
1237         float *o1, *o2, *o3, *o4;
1238
1239         v1= (mvert+mface->v1)->co;
1240         v2= (mvert+mface->v2)->co;
1241         v3= (mvert+mface->v3)->co;
1242         VECCOPY(n1,(mvert+mface->v1)->no);
1243         VECCOPY(n2,(mvert+mface->v2)->no);
1244         VECCOPY(n3,(mvert+mface->v3)->no);
1245         Normalize(n1);
1246         Normalize(n2);
1247         Normalize(n3);
1248
1249         if(mface->v4) {
1250                 v4= (mvert+mface->v4)->co;
1251                 VECCOPY(n4,(mvert+mface->v4)->no);
1252                 Normalize(n4);
1253                 
1254                 vec[0]= w[0]*v1[0] + w[1]*v2[0] + w[2]*v3[0] + w[3]*v4[0];
1255                 vec[1]= w[0]*v1[1] + w[1]*v2[1] + w[2]*v3[1] + w[3]*v4[1];
1256                 vec[2]= w[0]*v1[2] + w[1]*v2[2] + w[2]*v3[2] + w[3]*v4[2];
1257
1258                 if(nor){
1259                         if(mface->flag & ME_SMOOTH){
1260                                 nor[0]= w[0]*n1[0] + w[1]*n2[0] + w[2]*n3[0] + w[3]*n4[0];
1261                                 nor[1]= w[0]*n1[1] + w[1]*n2[1] + w[2]*n3[1] + w[3]*n4[1];
1262                                 nor[2]= w[0]*n1[2] + w[1]*n2[2] + w[2]*n3[2] + w[3]*n4[2];
1263                         }
1264                         else
1265                                 CalcNormFloat4(v1,v2,v3,v4,nor);
1266                 }
1267         }
1268         else {
1269                 vec[0]= w[0]*v1[0] + w[1]*v2[0] + w[2]*v3[0];
1270                 vec[1]= w[0]*v1[1] + w[1]*v2[1] + w[2]*v3[1];
1271                 vec[2]= w[0]*v1[2] + w[1]*v2[2] + w[2]*v3[2];
1272                 
1273                 if(nor){
1274                         if(mface->flag & ME_SMOOTH){
1275                                 nor[0]= w[0]*n1[0] + w[1]*n2[0] + w[2]*n3[0];
1276                                 nor[1]= w[0]*n1[1] + w[1]*n2[1] + w[2]*n3[1];
1277                                 nor[2]= w[0]*n1[2] + w[1]*n2[2] + w[2]*n3[2];
1278                         }
1279                         else
1280                                 CalcNormFloat(v1,v2,v3,nor);
1281                 }
1282         }
1283         
1284         /* calculate tangent vectors */
1285         if(utan && vtan){
1286                 if(tface){
1287                         uv1= tface->uv[0];
1288                         uv2= tface->uv[1];
1289                         uv3= tface->uv[2];
1290                         uv4= tface->uv[3];
1291                 }
1292                 else{
1293                         uv1= tuv[0]; uv2= tuv[1]; uv3= tuv[2]; uv4= tuv[3];
1294                         spheremap(v1[0], v1[1], v1[2], uv1, uv1+1);
1295                         spheremap(v2[0], v2[1], v2[2], uv2, uv2+1);
1296                         spheremap(v3[0], v3[1], v3[2], uv3, uv3+1);
1297                         if(v4)
1298                                 spheremap(v4[0], v4[1], v4[2], uv4, uv4+1);
1299                 }
1300
1301                 if(v4){
1302                         s1= uv3[0] - uv1[0];
1303                         s2= uv4[0] - uv1[0];
1304
1305                         t1= uv3[1] - uv1[1];
1306                         t2= uv4[1] - uv1[1];
1307
1308                         VecSubf(e1, v3, v1);
1309                         VecSubf(e2, v4, v1);
1310                 }
1311                 else{
1312                         s1= uv2[0] - uv1[0];
1313                         s2= uv3[0] - uv1[0];
1314
1315                         t1= uv2[1] - uv1[1];
1316                         t2= uv3[1] - uv1[1];
1317
1318                         VecSubf(e1, v2, v1);
1319                         VecSubf(e2, v3, v1);
1320                 }
1321
1322                 vtan[0] = (s1*e2[0] - s2*e1[0]);
1323                 vtan[1] = (s1*e2[1] - s2*e1[1]);
1324                 vtan[2] = (s1*e2[2] - s2*e1[2]);
1325
1326                 utan[0] = (t1*e2[0] - t2*e1[0]);
1327                 utan[1] = (t1*e2[1] - t2*e1[1]);
1328                 utan[2] = (t1*e2[2] - t2*e1[2]);
1329         }
1330
1331         if(orco) {
1332                 if(orcodata) {
1333                         o1= orcodata[mface->v1];
1334                         o2= orcodata[mface->v2];
1335                         o3= orcodata[mface->v3];
1336
1337                         if(mface->v4) {
1338                                 o4= orcodata[mface->v4];
1339                                 orco[0]= w[0]*o1[0] + w[1]*o2[0] + w[2]*o3[0] + w[3]*o4[0];
1340                                 orco[1]= w[0]*o1[1] + w[1]*o2[1] + w[2]*o3[1] + w[3]*o4[1];
1341                                 orco[2]= w[0]*o1[2] + w[1]*o2[2] + w[2]*o3[2] + w[3]*o4[2];
1342
1343                                 if(ornor)
1344                                         CalcNormFloat4(o1, o2, o3, o4, ornor);
1345                         }
1346                         else {
1347                                 orco[0]= w[0]*o1[0] + w[1]*o2[0] + w[2]*o3[0];
1348                                 orco[1]= w[0]*o1[1] + w[1]*o2[1] + w[2]*o3[1];
1349                                 orco[2]= w[0]*o1[2] + w[1]*o2[2] + w[2]*o3[2];
1350
1351                                 if(ornor)
1352                                         CalcNormFloat(o1, o2, o3, ornor);
1353                         }
1354                 }
1355                 else {
1356                         VECCOPY(orco, vec);
1357                         if(ornor)
1358                                 VECCOPY(ornor, nor);
1359                 }
1360         }
1361 }
1362 void psys_interpolate_uvs(MTFace *tface, int quad, float *w, float *uvco)
1363 {
1364         float v10= tface->uv[0][0];
1365         float v11= tface->uv[0][1];
1366         float v20= tface->uv[1][0];
1367         float v21= tface->uv[1][1];
1368         float v30= tface->uv[2][0];
1369         float v31= tface->uv[2][1];
1370         float v40,v41;
1371
1372         if(quad) {
1373                 v40= tface->uv[3][0];
1374                 v41= tface->uv[3][1];
1375
1376                 uvco[0]= w[0]*v10 + w[1]*v20 + w[2]*v30 + w[3]*v40;
1377                 uvco[1]= w[0]*v11 + w[1]*v21 + w[2]*v31 + w[3]*v41;
1378         }
1379         else {
1380                 uvco[0]= w[0]*v10 + w[1]*v20 + w[2]*v30;
1381                 uvco[1]= w[0]*v11 + w[1]*v21 + w[2]*v31;
1382         }
1383 }
1384
1385 void psys_interpolate_mcol(MCol *mcol, int quad, float *w, MCol *mc)
1386 {
1387         char *cp, *cp1, *cp2, *cp3, *cp4;
1388
1389         cp= (char *)mc;
1390         cp1= (char *)&mcol[0];
1391         cp2= (char *)&mcol[1];
1392         cp3= (char *)&mcol[2];
1393         
1394         if(quad) {
1395                 cp4= (char *)&mcol[3];
1396
1397                 cp[0]= (int)(w[0]*cp1[0] + w[1]*cp2[0] + w[2]*cp3[0] + w[3]*cp4[0]);
1398                 cp[1]= (int)(w[0]*cp1[1] + w[1]*cp2[1] + w[2]*cp3[1] + w[3]*cp4[1]);
1399                 cp[2]= (int)(w[0]*cp1[2] + w[1]*cp2[2] + w[2]*cp3[2] + w[3]*cp4[2]);
1400                 cp[3]= (int)(w[0]*cp1[3] + w[1]*cp2[3] + w[2]*cp3[3] + w[3]*cp4[3]);
1401         }
1402         else {
1403                 cp[0]= (int)(w[0]*cp1[0] + w[1]*cp2[0] + w[2]*cp3[0]);
1404                 cp[1]= (int)(w[0]*cp1[1] + w[1]*cp2[1] + w[2]*cp3[1]);
1405                 cp[2]= (int)(w[0]*cp1[2] + w[1]*cp2[2] + w[2]*cp3[2]);
1406                 cp[3]= (int)(w[0]*cp1[3] + w[1]*cp2[3] + w[2]*cp3[3]);
1407         }
1408 }
1409
1410 static float psys_interpolate_value_from_verts(DerivedMesh *dm, short from, int index, float *fw, float *values)
1411 {
1412         if(values==0 || index==-1)
1413                 return 0.0;
1414
1415         switch(from){
1416                 case PART_FROM_VERT:
1417                         return values[index];
1418                 case PART_FROM_FACE:
1419                 case PART_FROM_VOLUME:
1420                 {
1421                         MFace *mf=dm->getFaceData(dm,index,CD_MFACE);
1422                         return interpolate_particle_value(values[mf->v1],values[mf->v2],values[mf->v3],values[mf->v4],fw,mf->v4);
1423                 }
1424                         
1425         }
1426         return 0.0;
1427 }
1428
1429 /* conversion of pa->fw to origspace layer coordinates */
1430 static void psys_w_to_origspace(float *w, float *uv)
1431 {
1432         uv[0]= w[1] + w[2];
1433         uv[1]= w[2] + w[3];
1434 }
1435
1436 /* conversion of pa->fw to weights in face from origspace */
1437 static void psys_origspace_to_w(OrigSpaceFace *osface, int quad, float *w, float *neww)
1438 {
1439         float v[4][3], co[3];
1440
1441         v[0][0]= osface->uv[0][0]; v[0][1]= osface->uv[0][1]; v[0][2]= 0.0f;
1442         v[1][0]= osface->uv[1][0]; v[1][1]= osface->uv[1][1]; v[1][2]= 0.0f;
1443         v[2][0]= osface->uv[2][0]; v[2][1]= osface->uv[2][1]; v[2][2]= 0.0f;
1444
1445         psys_w_to_origspace(w, co);
1446         co[2]= 0.0f;
1447         
1448         if(quad) {
1449                 v[3][0]= osface->uv[3][0]; v[3][1]= osface->uv[3][1]; v[3][2]= 0.0f;
1450                 MeanValueWeights(v, 4, co, neww);
1451         }
1452         else {
1453                 MeanValueWeights(v, 3, co, neww);
1454                 neww[3]= 0.0f;
1455         }
1456 }
1457
1458 /* find the derived mesh face for a particle, set the mf passed. this is slow
1459  * and can be optimized but only for many lookups. returns the face index. */
1460 int psys_particle_dm_face_lookup(Object *ob, DerivedMesh *dm, int index, float *fw, struct LinkNode *node)
1461 {
1462         Mesh *me= (Mesh*)ob->data;
1463         MFace *mface;
1464         OrigSpaceFace *osface;
1465         int *origindex;
1466         int quad, findex, totface;
1467         float uv[2], (*faceuv)[2];
1468
1469         mface = dm->getFaceDataArray(dm, CD_MFACE);
1470         origindex = dm->getFaceDataArray(dm, CD_ORIGINDEX);
1471         osface = dm->getFaceDataArray(dm, CD_ORIGSPACE);
1472
1473         totface = dm->getNumFaces(dm);
1474         
1475         if(osface==NULL || origindex==NULL) {
1476                 /* Assume we dont need osface data */
1477                 if (index <totface) {
1478                         //printf("\tNO CD_ORIGSPACE, assuming not needed\n");
1479                         return index;
1480                 } else {
1481                         printf("\tNO CD_ORIGSPACE, error out of range\n");
1482                         return DMCACHE_NOTFOUND;
1483                 }
1484         }
1485         else if(index >= me->totface)
1486                 return DMCACHE_NOTFOUND; /* index not in the original mesh */
1487
1488         psys_w_to_origspace(fw, uv);
1489         
1490         if(node) { /* we have a linked list of faces that we use, faster! */
1491                 for(;node; node=node->next) {
1492                         findex= GET_INT_FROM_POINTER(node->link);
1493                         faceuv= osface[findex].uv;
1494                         quad= mface[findex].v4;
1495
1496                         /* check that this intersects - Its possible this misses :/ -
1497                          * could also check its not between */
1498                         if(quad) {
1499                                 if(IsectPQ2Df(uv, faceuv[0], faceuv[1], faceuv[2], faceuv[3]))
1500                                         return findex;
1501                         }
1502                         else if(IsectPT2Df(uv, faceuv[0], faceuv[1], faceuv[2]))
1503                                 return findex;
1504                 }
1505         }
1506         else { /* if we have no node, try every face */
1507                 for(findex=0; findex<totface; findex++) {
1508                         if(origindex[findex] == index) {
1509                                 faceuv= osface[findex].uv;
1510                                 quad= mface[findex].v4;
1511
1512                                 /* check that this intersects - Its possible this misses :/ -
1513                                  * could also check its not between */
1514                                 if(quad) {
1515                                         if(IsectPQ2Df(uv, faceuv[0], faceuv[1], faceuv[2], faceuv[3]))
1516                                                 return findex;
1517                                 }
1518                                 else if(IsectPT2Df(uv, faceuv[0], faceuv[1], faceuv[2]))
1519                                         return findex;
1520                         }
1521                 }
1522         }
1523
1524         return DMCACHE_NOTFOUND;
1525 }
1526
1527 static int psys_map_index_on_dm(DerivedMesh *dm, int from, int index, int index_dmcache, float *fw, float foffset, int *mapindex, float *mapfw)
1528 {
1529         if(index < 0)
1530                 return 0;
1531
1532         if (dm->deformedOnly || index_dmcache == DMCACHE_ISCHILD) {
1533                 /* for meshes that are either only defined or for child particles, the
1534                  * index and fw do not require any mapping, so we can directly use it */
1535                 if(from == PART_FROM_VERT) {
1536                         if(index >= dm->getNumVerts(dm))
1537                                 return 0;
1538
1539                         *mapindex = index;
1540                 }
1541                 else  { /* FROM_FACE/FROM_VOLUME */
1542                         if(index >= dm->getNumFaces(dm))
1543                                 return 0;
1544
1545                         *mapindex = index;
1546                         QUATCOPY(mapfw, fw);
1547                 }
1548         } else {
1549                 /* for other meshes that have been modified, we try to map the particle
1550                  * to their new location, which means a different index, and for faces
1551                  * also a new face interpolation weights */
1552                 if(from == PART_FROM_VERT) {
1553                         if (index_dmcache == DMCACHE_NOTFOUND || index_dmcache > dm->getNumVerts(dm))
1554                                 return 0;
1555
1556                         *mapindex = index_dmcache;
1557                 }
1558                 else  { /* FROM_FACE/FROM_VOLUME */
1559                         /* find a face on the derived mesh that uses this face */
1560                         MFace *mface;
1561                         OrigSpaceFace *osface;
1562                         int i;
1563
1564                         i = index_dmcache;
1565
1566                         if(i== DMCACHE_NOTFOUND || i >= dm->getNumFaces(dm))
1567                                 return 0;
1568
1569                         *mapindex = i;
1570
1571                         /* modify the original weights to become
1572                          * weights for the derived mesh face */
1573                         osface= dm->getFaceDataArray(dm, CD_ORIGSPACE);
1574                         mface= dm->getFaceData(dm, i, CD_MFACE);
1575
1576                         if(osface == NULL)
1577                                 mapfw[0]= mapfw[1]= mapfw[2]= mapfw[3]= 0.0f;
1578                         else
1579                                 psys_origspace_to_w(&osface[i], mface->v4, fw, mapfw);
1580                 }
1581         }
1582
1583         return 1;
1584 }
1585
1586 /* interprets particle data to get a point on a mesh in object space */
1587 void psys_particle_on_dm(DerivedMesh *dm, int from, int index, int index_dmcache, float *fw, float foffset, float *vec, float *nor, float *utan, float *vtan, float *orco, float *ornor)
1588 {
1589         float tmpnor[3], mapfw[4];
1590         float (*orcodata)[3];
1591         int mapindex;
1592
1593         if(!psys_map_index_on_dm(dm, from, index, index_dmcache, fw, foffset, &mapindex, mapfw)) {
1594                 if(vec) { vec[0]=vec[1]=vec[2]=0.0; }
1595                 if(nor) { nor[0]=nor[1]=0.0; nor[2]=1.0; }
1596                 if(orco) { orco[0]=orco[1]=orco[2]=0.0; }
1597                 if(ornor) { ornor[0]=ornor[1]=0.0; ornor[2]=1.0; }
1598                 if(utan) { utan[0]=utan[1]=utan[2]=0.0; }
1599                 if(vtan) { vtan[0]=vtan[1]=vtan[2]=0.0; }
1600
1601                 return;
1602         }
1603
1604         orcodata= dm->getVertDataArray(dm, CD_ORCO);
1605
1606         if(from == PART_FROM_VERT) {
1607                 dm->getVertCo(dm,mapindex,vec);
1608
1609                 if(nor) {
1610                         dm->getVertNo(dm,mapindex,nor);
1611                         Normalize(nor);
1612                 }
1613
1614                 if(orco)
1615                         VECCOPY(orco, orcodata[mapindex])
1616
1617                 if(ornor) {
1618                         dm->getVertNo(dm,mapindex,nor);
1619                         Normalize(nor);
1620                 }
1621
1622                 if(utan && vtan) {
1623                         utan[0]= utan[1]= utan[2]= 0.0f;
1624                         vtan[0]= vtan[1]= vtan[2]= 0.0f;
1625                 }
1626         }
1627         else { /* PART_FROM_FACE / PART_FROM_VOLUME */
1628                 MFace *mface;
1629                 MTFace *mtface;
1630                 MVert *mvert;
1631
1632                 mface=dm->getFaceData(dm,mapindex,CD_MFACE);
1633                 mvert=dm->getVertDataArray(dm,CD_MVERT);
1634                 mtface=CustomData_get_layer(&dm->faceData,CD_MTFACE);
1635
1636                 if(mtface)
1637                         mtface += mapindex;
1638
1639                 if(from==PART_FROM_VOLUME) {
1640                         psys_interpolate_face(mvert,mface,mtface,orcodata,mapfw,vec,tmpnor,utan,vtan,orco,ornor);
1641                         if(nor)
1642                                 VECCOPY(nor,tmpnor);
1643
1644                         Normalize(tmpnor);
1645                         VecMulf(tmpnor,-foffset);
1646                         VECADD(vec,vec,tmpnor);
1647                 }
1648                 else
1649                         psys_interpolate_face(mvert,mface,mtface,orcodata,mapfw,vec,nor,utan,vtan,orco,ornor);
1650         }
1651 }
1652
1653 float psys_particle_value_from_verts(DerivedMesh *dm, short from, ParticleData *pa, float *values)
1654 {
1655         float mapfw[4];
1656         int mapindex;
1657
1658         if(!psys_map_index_on_dm(dm, from, pa->num, pa->num_dmcache, pa->fuv, pa->foffset, &mapindex, mapfw))
1659                 return 0.0f;
1660         
1661         return psys_interpolate_value_from_verts(dm, from, mapindex, mapfw, values);
1662 }
1663
1664 ParticleSystemModifierData *psys_get_modifier(Object *ob, ParticleSystem *psys)
1665 {
1666         ModifierData *md;
1667         ParticleSystemModifierData *psmd;
1668
1669         for(md=ob->modifiers.first; md; md=md->next){
1670                 if(md->type==eModifierType_ParticleSystem){
1671                         psmd= (ParticleSystemModifierData*) md;
1672                         if(psmd->psys==psys){
1673                                 return psmd;
1674                         }
1675                 }
1676         }
1677         return NULL;
1678 }
1679 /************************************************/
1680 /*                      Particles on a shape                            */
1681 /************************************************/
1682 /* ready for future use */
1683 static void psys_particle_on_shape(int distr, int index, float *fuv, float *vec, float *nor, float *utan, float *vtan, float *orco, float *ornor)
1684 {
1685         /* TODO */
1686         float zerovec[3]={0.0f,0.0f,0.0f};
1687         if(vec){
1688                 VECCOPY(vec,zerovec);
1689         }
1690         if(nor){
1691                 VECCOPY(nor,zerovec);
1692         }
1693         if(utan){
1694                 VECCOPY(utan,zerovec);
1695         }
1696         if(vtan){
1697                 VECCOPY(vtan,zerovec);
1698         }
1699         if(orco){
1700                 VECCOPY(orco,zerovec);
1701         }
1702         if(ornor){
1703                 VECCOPY(ornor,zerovec);
1704         }
1705 }
1706 /************************************************/
1707 /*                      Particles on emitter                            */
1708 /************************************************/
1709 void psys_particle_on_emitter(ParticleSystemModifierData *psmd, int from, int index, int index_dmcache, float *fuv, float foffset, float *vec, float *nor, float *utan, float *vtan, float *orco, float *ornor){
1710         if(psmd){
1711                 if(psmd->psys->part->distr==PART_DISTR_GRID && psmd->psys->part->from != PART_FROM_VERT){
1712                         if(vec){
1713                                 VECCOPY(vec,fuv);
1714                         }
1715                         return;
1716                 }
1717                 /* we cant use the num_dmcache */
1718                 psys_particle_on_dm(psmd->dm,from,index,index_dmcache,fuv,foffset,vec,nor,utan,vtan,orco,ornor);
1719         }
1720         else
1721                 psys_particle_on_shape(from,index,fuv,vec,nor,utan,vtan,orco,ornor);
1722
1723 }
1724 /************************************************/
1725 /*                      Path Cache                                                      */
1726 /************************************************/
1727 static float vert_weight(MDeformVert *dvert, int group)
1728 {
1729         MDeformWeight *dw;
1730         int i;
1731         
1732         if(dvert) {
1733                 dw= dvert->dw;
1734                 for(i= dvert->totweight; i>0; i--, dw++) {
1735                         if(dw->def_nr == group) return dw->weight;
1736                         if(i==1) break; /*otherwise dw will point to somewhere it shouldn't*/
1737                 }
1738         }
1739         return 0.0;
1740 }
1741
1742 static void do_prekink(ParticleKey *state, ParticleKey *par, float *par_rot, float time, float freq, float shape, float amplitude, short type, short axis, float obmat[][4])
1743 {
1744         float vec[3]={0.0,0.0,0.0}, q1[4]={1,0,0,0},q2[4];
1745         float t;
1746
1747         CLAMP(time,0.0,1.0);
1748
1749         if(shape!=0.0f && type!=PART_KINK_BRAID) {
1750                 if(shape<0.0f)
1751                         time= (float)pow(time, 1.0+shape);
1752                 else
1753                         time= (float)pow(time, 1.0/(1.0-shape));
1754         }
1755
1756         t=time;
1757
1758         t*=(float)M_PI*freq;
1759
1760         if(par==0) return;
1761
1762         switch(type){
1763                 case PART_KINK_CURL:
1764                         vec[axis]=1.0;
1765                         if(par_rot)
1766                                 QUATCOPY(q2,par_rot)
1767                         else
1768                                 vectoquat(par->vel,axis,(axis+1)%3, q2);
1769                         QuatMulVecf(q2,vec);
1770                         VecMulf(vec,amplitude);
1771                         VECADD(state->co,state->co,vec);
1772
1773                         VECSUB(vec,state->co,par->co);
1774
1775                         if(t!=0.0)
1776                                 VecRotToQuat(par->vel,t,q1);
1777                         
1778                         QuatMulVecf(q1,vec);
1779                         
1780                         VECADD(state->co,par->co,vec);
1781                         break;
1782                 case PART_KINK_RADIAL:
1783                         VECSUB(vec,state->co,par->co);
1784
1785                         Normalize(vec);
1786                         VecMulf(vec,amplitude*(float)sin(t));
1787
1788                         VECADD(state->co,state->co,vec);
1789                         break;
1790                 case PART_KINK_WAVE:
1791                         vec[axis]=1.0;
1792                         if(obmat)
1793                                 Mat4Mul3Vecfl(obmat,vec);
1794
1795                         if(par_rot)
1796                                 QuatMulVecf(par_rot,vec);
1797
1798                         Projf(q1,vec,par->vel);
1799                         
1800                         VECSUB(vec,vec,q1);
1801                         Normalize(vec);
1802
1803                         VecMulf(vec,amplitude*(float)sin(t));
1804
1805                         VECADD(state->co,state->co,vec);
1806                         break;
1807                 case PART_KINK_BRAID:
1808                         if(par){
1809                                 float y_vec[3]={0.0,1.0,0.0};
1810                                 float z_vec[3]={0.0,0.0,1.0};
1811                                 float vec_from_par[3], vec_one[3], radius, state_co[3];
1812                                 float inp_y,inp_z,length;
1813                                 
1814                                 if(par_rot)
1815                                         QUATCOPY(q2,par_rot)
1816                                 else
1817                                         vectoquat(par->vel,axis,(axis+1)%3,q2);
1818                                 QuatMulVecf(q2,y_vec);
1819                                 QuatMulVecf(q2,z_vec);
1820                                 
1821                                 VECSUB(vec_from_par,state->co,par->co);
1822                                 VECCOPY(vec_one,vec_from_par);
1823                                 radius=Normalize(vec_one);
1824
1825                                 inp_y=Inpf(y_vec,vec_one);
1826                                 inp_z=Inpf(z_vec,vec_one);
1827
1828                                 if(inp_y>0.5){
1829                                         VECCOPY(state_co,y_vec);
1830
1831                                         VecMulf(y_vec,amplitude*(float)cos(t));
1832                                         VecMulf(z_vec,amplitude/2.0f*(float)sin(2.0f*t));
1833                                 }
1834                                 else if(inp_z>0.0){
1835                                         VECCOPY(state_co,z_vec);
1836                                         VecMulf(state_co,(float)sin(M_PI/3.0f));
1837                                         VECADDFAC(state_co,state_co,y_vec,-0.5f);
1838
1839                                         VecMulf(y_vec,-amplitude*(float)cos(t + M_PI/3.0f));
1840                                         VecMulf(z_vec,amplitude/2.0f*(float)cos(2.0f*t + M_PI/6.0f));
1841                                 }
1842                                 else{
1843                                         VECCOPY(state_co,z_vec);
1844                                         VecMulf(state_co,-(float)sin(M_PI/3.0f));
1845                                         VECADDFAC(state_co,state_co,y_vec,-0.5f);
1846
1847                                         VecMulf(y_vec,amplitude*(float)-sin(t+M_PI/6.0f));
1848                                         VecMulf(z_vec,amplitude/2.0f*(float)-sin(2.0f*t+M_PI/3.0f));
1849                                 }
1850
1851                                 VecMulf(state_co,amplitude);
1852                                 VECADD(state_co,state_co,par->co);
1853                                 VECSUB(vec_from_par,state->co,state_co);
1854
1855                                 length=Normalize(vec_from_par);
1856                                 VecMulf(vec_from_par,MIN2(length,amplitude/2.0f));
1857
1858                                 VECADD(state_co,par->co,y_vec);
1859                                 VECADD(state_co,state_co,z_vec);
1860                                 VECADD(state_co,state_co,vec_from_par);
1861
1862                                 shape=(2.0f*(float)M_PI)*(1.0f+shape);
1863
1864                                 if(t<shape){
1865                                         shape=t/shape;
1866                                         shape=(float)sqrt((double)shape);
1867                                         VecLerpf(state->co,state->co,state_co,shape);
1868                                 }
1869                                 else{
1870                                         VECCOPY(state->co,state_co);
1871                                 }
1872                         }
1873                         break;
1874         }
1875 }
1876
1877 static void do_clump(ParticleKey *state, ParticleKey *par, float time, float clumpfac, float clumppow, float pa_clump)
1878 {
1879         if(par && clumpfac!=0.0){
1880                 float clump, cpow;
1881
1882                 if(clumppow<0.0)
1883                         cpow=1.0f+clumppow;
1884                 else
1885                         cpow=1.0f+9.0f*clumppow;
1886
1887                 if(clumpfac<0.0) /* clump roots instead of tips */
1888                         clump = -clumpfac*pa_clump*(float)pow(1.0-(double)time,(double)cpow);
1889                 else
1890                         clump = clumpfac*pa_clump*(float)pow((double)time,(double)cpow);
1891                 VecLerpf(state->co,state->co,par->co,clump);
1892         }
1893 }
1894 void precalc_guides(ParticleSimulationData *sim, ListBase *effectors)
1895 {
1896         EffectedPoint point;
1897         ParticleKey state;
1898         EffectorData efd;
1899         EffectorCache *eff;
1900         ParticleSystem *psys = sim->psys;
1901         EffectorWeights *weights = sim->psys->part->effector_weights;
1902         GuideEffectorData *data;
1903         PARTICLE_P;
1904
1905         if(!effectors)
1906                 return;
1907
1908         LOOP_PARTICLES {
1909                 psys_particle_on_emitter(sim->psmd,sim->psys->part->from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,state.co,0,0,0,0,0);
1910                 pd_point_from_particle(sim, pa, &state, &point);
1911
1912                 for(eff = effectors->first; eff; eff=eff->next) {
1913                         if(eff->pd->forcefield != PFIELD_GUIDE)
1914                                 continue;
1915
1916                         if(!eff->guide_data)
1917                                 eff->guide_data = MEM_callocN(sizeof(GuideEffectorData)*psys->totpart, "GuideEffectorData");
1918
1919                         data = eff->guide_data + p;
1920
1921                         VECSUB(efd.vec_to_point, state.co, eff->guide_loc);
1922                         VECCOPY(efd.nor, eff->guide_dir);
1923                         efd.distance = VecLength(efd.vec_to_point);
1924
1925                         VECCOPY(data->vec_to_point, efd.vec_to_point);
1926                         data->strength = effector_falloff(eff, &efd, &point, weights);
1927                 }
1928         }
1929 }
1930 int do_guides(ListBase *effectors, ParticleKey *state, int index, float time)
1931 {
1932         EffectorCache *eff;
1933         PartDeflect *pd;
1934         Curve *cu;
1935         ParticleKey key, par;
1936         GuideEffectorData *data;
1937
1938         float effect[3] = {0.0f, 0.0f, 0.0f}, veffect[3] = {0.0f, 0.0f, 0.0f};
1939         float guidevec[4], guidedir[3], rot2[4], temp[3];
1940         float guidetime, radius, angle, totstrength = 0.0f;
1941         float vec_to_point[3];
1942
1943         if(effectors) for(eff = effectors->first; eff; eff=eff->next) {
1944                 pd = eff->pd;
1945
1946                 if(pd->forcefield != PFIELD_GUIDE)
1947                         continue;
1948
1949                 data = eff->guide_data + index;
1950
1951                 if(data->strength <= 0.0f)
1952                         continue;
1953
1954                 guidetime = time / (1.0 - pd->free_end);
1955
1956                 if(guidetime>1.0f)
1957                         continue;
1958
1959                 cu = (Curve*)eff->ob->data;
1960
1961                 if(pd->flag & PFIELD_GUIDE_PATH_ADD) {
1962                         if(where_on_path(eff->ob, data->strength * guidetime, guidevec, guidedir, NULL, &radius)==0)
1963                                 return 0;
1964                 }
1965                 else {
1966                         if(where_on_path(eff->ob, guidetime, guidevec, guidedir, NULL, &radius)==0)
1967                                 return 0;
1968                 }
1969
1970                 Mat4MulVecfl(eff->ob->obmat, guidevec);
1971                 Mat4Mul3Vecfl(eff->ob->obmat, guidedir);
1972
1973                 Normalize(guidedir);
1974
1975                 VECCOPY(vec_to_point, data->vec_to_point);
1976
1977                 if(guidetime != 0.0){
1978                         /* curve direction */
1979                         Crossf(temp, eff->guide_dir, guidedir);
1980                         angle = Inpf(eff->guide_dir, guidedir)/(VecLength(eff->guide_dir));
1981                         angle = saacos(angle);
1982                         VecRotToQuat(temp, angle, rot2);
1983                         QuatMulVecf(rot2, vec_to_point);
1984
1985                         /* curve tilt */
1986                         VecRotToQuat(guidedir, guidevec[3] - eff->guide_loc[3], rot2);
1987                         QuatMulVecf(rot2, vec_to_point);
1988                 }
1989
1990                 /* curve taper */
1991                 if(cu->taperobj)
1992                         VecMulf(vec_to_point, calc_taper(eff->scene, cu->taperobj, (int)(data->strength*guidetime*100.0), 100));
1993
1994                 else{ /* curve size*/
1995                         if(cu->flag & CU_PATH_RADIUS) {
1996                                 VecMulf(vec_to_point, radius);
1997                         }
1998                 }
1999                 par.co[0] = par.co[1] = par.co[2] = 0.0f;
2000                 VECCOPY(key.co, vec_to_point);
2001                 do_prekink(&key, &par, 0, guidetime, pd->kink_freq, pd->kink_shape, pd->kink_amp, pd->kink, pd->kink_axis, 0);
2002                 do_clump(&key, &par, guidetime, pd->clump_fac, pd->clump_pow, 1.0f);
2003                 VECCOPY(vec_to_point, key.co);
2004
2005                 VECADD(vec_to_point, vec_to_point, guidevec);
2006                 //VECSUB(pa_loc,pa_loc,pa_zero);
2007                 VECADDFAC(effect, effect, vec_to_point, data->strength);
2008                 VECADDFAC(veffect, veffect, guidedir, data->strength);
2009                 totstrength += data->strength;
2010         }
2011
2012         if(totstrength != 0.0){
2013                 if(totstrength > 1.0)
2014                         VecMulf(effect, 1.0f / totstrength);
2015                 CLAMP(totstrength, 0.0, 1.0);
2016                 //VECADD(effect,effect,pa_zero);
2017                 VecLerpf(state->co, state->co, effect, totstrength);
2018
2019                 Normalize(veffect);
2020                 VecMulf(veffect, VecLength(state->vel));
2021                 VECCOPY(state->vel, veffect);
2022                 return 1;
2023         }
2024         return 0;
2025 }
2026 static void do_rough(float *loc, float mat[4][4], float t, float fac, float size, float thres, ParticleKey *state)
2027 {
2028         float rough[3];
2029         float rco[3];
2030
2031         if(thres!=0.0)
2032                 if((float)fabs((float)(-1.5+loc[0]+loc[1]+loc[2]))<1.5f*thres) return;
2033
2034         VECCOPY(rco,loc);
2035         VecMulf(rco,t);
2036         rough[0]=-1.0f+2.0f*BLI_gTurbulence(size, rco[0], rco[1], rco[2], 2,0,2);
2037         rough[1]=-1.0f+2.0f*BLI_gTurbulence(size, rco[1], rco[2], rco[0], 2,0,2);
2038         rough[2]=-1.0f+2.0f*BLI_gTurbulence(size, rco[2], rco[0], rco[1], 2,0,2);
2039
2040         VECADDFAC(state->co,state->co,mat[0],fac*rough[0]);
2041         VECADDFAC(state->co,state->co,mat[1],fac*rough[1]);
2042         VECADDFAC(state->co,state->co,mat[2],fac*rough[2]);
2043 }
2044 static void do_rough_end(float *loc, float mat[4][4], float t, float fac, float shape, ParticleKey *state)
2045 {
2046         float rough[2];
2047         float roughfac;
2048
2049         roughfac=fac*(float)pow((double)t,shape);
2050         VECCOPY(rough,loc);
2051         rough[0]=-1.0f+2.0f*rough[0];
2052         rough[1]=-1.0f+2.0f*rough[1];
2053         VecMulf(rough,roughfac);
2054
2055         VECADDFAC(state->co,state->co,mat[0],rough[0]);
2056         VECADDFAC(state->co,state->co,mat[1],rough[1]);
2057 }
2058 static void do_path_effectors(ParticleSimulationData *sim, int i, ParticleCacheKey *ca, int k, int steps, float *rootco, float effector, float dfra, float cfra, float *length, float *vec)
2059 {
2060         float force[3] = {0.0f,0.0f,0.0f};
2061         ParticleKey eff_key;
2062         EffectedPoint epoint;
2063
2064         /* Don't apply effectors for dynamic hair, otherwise the effectors don't get applied twice. */
2065         if(sim->psys->flag & PSYS_HAIR_DYNAMICS)
2066                 return;
2067
2068         VECCOPY(eff_key.co,(ca-1)->co);
2069         VECCOPY(eff_key.vel,(ca-1)->vel);
2070         QUATCOPY(eff_key.rot,(ca-1)->rot);
2071
2072         pd_point_from_particle(sim, sim->psys->particles+i, &eff_key, &epoint);
2073         pdDoEffectors(sim->psys->effectors, sim->colliders, sim->psys->part->effector_weights, &epoint, force, NULL);
2074
2075         VecMulf(force, effector*pow((float)k / (float)steps, 100.0f * sim->psys->part->eff_hair) / (float)steps);
2076
2077         VecAddf(force, force, vec);
2078
2079         Normalize(force);
2080
2081         VECADDFAC(ca->co, (ca-1)->co, force, *length);
2082
2083         if(k < steps) {
2084                 VecSubf(vec, (ca+1)->co, ca->co);
2085                 *length = VecLength(vec);
2086         }
2087 }
2088 static int check_path_length(int k, ParticleCacheKey *keys, ParticleCacheKey *state, float max_length, float *cur_length, float length, float *dvec)
2089 {
2090         if(*cur_length + length > max_length){
2091                 VecMulf(dvec, (max_length - *cur_length) / length);
2092                 VECADD(state->co, (state - 1)->co, dvec);
2093                 keys->steps = k;
2094                 /* something over the maximum step value */
2095                 return k=100000;
2096         }
2097         else {
2098                 *cur_length+=length;
2099                 return k;
2100         }
2101 }
2102 static void offset_child(ChildParticle *cpa, ParticleKey *par, ParticleKey *child, float flat, float radius)
2103 {
2104         VECCOPY(child->co,cpa->fuv);
2105         VecMulf(child->co,radius);
2106
2107         child->co[0]*=flat;
2108
2109         VECCOPY(child->vel,par->vel);
2110
2111         QuatMulVecf(par->rot,child->co);
2112
2113         QUATCOPY(child->rot,par->rot);
2114
2115         VECADD(child->co,child->co,par->co);
2116 }
2117 float *psys_cache_vgroup(DerivedMesh *dm, ParticleSystem *psys, int vgroup)
2118 {
2119         float *vg=0;
2120
2121         if(vgroup < 0) {
2122                 /* hair dynamics pinning vgroup */
2123
2124         }
2125         else if(psys->vgroup[vgroup]){
2126                 MDeformVert *dvert = dm->getVertDataArray(dm, CD_MDEFORMVERT);
2127                 if(dvert){
2128                         int totvert=dm->getNumVerts(dm), i;
2129                         vg=MEM_callocN(sizeof(float)*totvert, "vg_cache");
2130                         if(psys->vg_neg&(1<<vgroup)){
2131                                 for(i=0; i<totvert; i++)
2132                                         vg[i]=1.0f-vert_weight(dvert+i,psys->vgroup[vgroup]-1);
2133                         }
2134                         else{
2135                                 for(i=0; i<totvert; i++)
2136                                         vg[i]=vert_weight(dvert+i,psys->vgroup[vgroup]-1);
2137                         }
2138                 }
2139         }
2140         return vg;
2141 }
2142 void psys_find_parents(ParticleSimulationData *sim)
2143 {
2144         ParticleSettings *part=sim->psys->part;
2145         KDTree *tree;
2146         ChildParticle *cpa;
2147         int p, totparent,totchild=sim->psys->totchild;
2148         float co[3], orco[3];
2149         int from=PART_FROM_FACE;
2150         totparent=(int)(totchild*part->parents*0.3);
2151
2152         if(G.rendering && part->child_nbr && part->ren_child_nbr)
2153                 totparent*=(float)part->child_nbr/(float)part->ren_child_nbr;
2154
2155         tree=BLI_kdtree_new(totparent);
2156
2157         for(p=0,cpa=sim->psys->child; p<totparent; p++,cpa++){
2158                 psys_particle_on_emitter(sim->psmd,from,cpa->num,DMCACHE_ISCHILD,cpa->fuv,cpa->foffset,co,0,0,0,orco,0);
2159                 BLI_kdtree_insert(tree, p, orco, NULL);
2160         }
2161
2162         BLI_kdtree_balance(tree);
2163
2164         for(; p<totchild; p++,cpa++){
2165                 psys_particle_on_emitter(sim->psmd,from,cpa->num,DMCACHE_ISCHILD,cpa->fuv,cpa->foffset,co,0,0,0,orco,0);
2166                 cpa->parent=BLI_kdtree_find_nearest(tree, orco, NULL, NULL);
2167         }
2168
2169         BLI_kdtree_free(tree);
2170 }
2171
2172 static void get_strand_normal(Material *ma, float *surfnor, float surfdist, float *nor)
2173 {
2174         float cross[3], nstrand[3], vnor[3], blend;
2175
2176         if(!((ma->mode & MA_STR_SURFDIFF) || (ma->strand_surfnor > 0.0f)))
2177                 return;
2178
2179         if(ma->mode & MA_STR_SURFDIFF) {
2180                 Crossf(cross, surfnor, nor);
2181                 Crossf(nstrand, nor, cross);
2182
2183                 blend= INPR(nstrand, surfnor);
2184                 CLAMP(blend, 0.0f, 1.0f);
2185
2186                 VecLerpf(vnor, nstrand, surfnor, blend);
2187                 Normalize(vnor);
2188         }
2189         else
2190                 VECCOPY(vnor, nor)
2191         
2192         if(ma->strand_surfnor > 0.0f) {
2193                 if(ma->strand_surfnor > surfdist) {
2194                         blend= (ma->strand_surfnor - surfdist)/ma->strand_surfnor;
2195                         VecLerpf(vnor, vnor, surfnor, blend);
2196                         Normalize(vnor);
2197                 }
2198         }
2199
2200         VECCOPY(nor, vnor);
2201 }
2202
2203 static int psys_threads_init_path(ParticleThread *threads, Scene *scene, float cfra, int editupdate)
2204 {
2205         ParticleThreadContext *ctx= threads[0].ctx;
2206 /*      Object *ob= ctx->sim.ob; */
2207         ParticleSystem *psys= ctx->sim.psys;
2208         ParticleSettings *part = psys->part;
2209         ParticleEditSettings *pset = &scene->toolsettings->particle;
2210         int totparent=0, between=0;
2211         int steps = (int)pow(2.0,(double)part->draw_step);
2212         int totchild = psys->totchild;
2213         int i, seed, totthread= threads[0].tot;
2214
2215         /*---start figuring out what is actually wanted---*/
2216         if(psys_in_edit_mode(scene, psys))
2217                 if(psys->renderdata==0 && (psys->edit==NULL || pset->flag & PE_DRAW_PART)==0)
2218                         totchild=0;
2219
2220         if(totchild && part->from!=PART_FROM_PARTICLE && part->childtype==PART_CHILD_FACES){
2221                 totparent=(int)(totchild*part->parents*0.3);
2222                 
2223                 if(G.rendering && part->child_nbr && part->ren_child_nbr)
2224                         totparent*=(float)part->child_nbr/(float)part->ren_child_nbr;
2225
2226                 /* part->parents could still be 0 so we can't test with totparent */
2227                 between=1;
2228         }
2229
2230         if(psys->renderdata)
2231                 steps=(int)pow(2.0,(double)part->ren_step);
2232         else{
2233                 totchild=(int)((float)totchild*(float)part->disp/100.0f);
2234                 totparent=MIN2(totparent,totchild);
2235         }
2236
2237         if(totchild==0) return 0;
2238
2239         /* init random number generator */
2240         if(ctx->sim.psys->part->flag & PART_ANIM_BRANCHING)
2241                 seed= 31415926 + ctx->sim.psys->seed + (int)cfra;
2242         else
2243                 seed= 31415926 + ctx->sim.psys->seed;
2244         
2245         if(part->flag & PART_BRANCHING || ctx->editupdate || totchild < 10000)
2246                 totthread= 1;
2247         
2248         for(i=0; i<totthread; i++) {
2249                 threads[i].rng_path= rng_new(seed);
2250                 threads[i].tot= totthread;
2251         }
2252
2253         /* fill context values */
2254         ctx->between= between;
2255         ctx->steps= steps;
2256         ctx->totchild= totchild;
2257         ctx->totparent= totparent;
2258         ctx->parent_pass= 0;
2259         ctx->cfra= cfra;
2260
2261         psys->lattice = psys_get_lattice(&ctx->sim);
2262
2263         /* cache all relevant vertex groups if they exist */
2264         if(part->from!=PART_FROM_PARTICLE){
2265                 ctx->vg_length = psys_cache_vgroup(ctx->dm,psys,PSYS_VG_LENGTH);
2266                 ctx->vg_clump = psys_cache_vgroup(ctx->dm,psys,PSYS_VG_CLUMP);
2267                 ctx->vg_kink = psys_cache_vgroup(ctx->dm,psys,PSYS_VG_KINK);
2268                 ctx->vg_rough1 = psys_cache_vgroup(ctx->dm,psys,PSYS_VG_ROUGH1);
2269                 ctx->vg_rough2 = psys_cache_vgroup(ctx->dm,psys,PSYS_VG_ROUGH2);
2270                 ctx->vg_roughe = psys_cache_vgroup(ctx->dm,psys,PSYS_VG_ROUGHE);
2271                 if(psys->part->flag & PART_CHILD_EFFECT)
2272                         ctx->vg_effector = psys_cache_vgroup(ctx->dm,psys,PSYS_VG_EFFECTOR);
2273         }
2274
2275         /* set correct ipo timing */
2276 #if 0 // XXX old animation system
2277         if(part->flag&PART_ABS_TIME && part->ipo){
2278                 calc_ipo(part->ipo, cfra);
2279                 execute_ipo((ID *)part, part->ipo);
2280         }
2281 #endif // XXX old animation system
2282
2283         return 1;
2284 }
2285
2286 /* note: this function must be thread safe, except for branching! */
2287 static void psys_thread_create_path(ParticleThread *thread, struct ChildParticle *cpa, ParticleCacheKey *keys, int i)
2288 {
2289         ParticleThreadContext *ctx= thread->ctx;
2290         Object *ob= ctx->sim.ob;
2291         ParticleSystem *psys = ctx->sim.psys;
2292         ParticleSettings *part = psys->part;
2293         ParticleCacheKey **cache= psys->childcache;
2294         ParticleCacheKey **pcache= psys->pathcache;
2295         ParticleCacheKey *state, *par = NULL, *key[4];
2296         ParticleData *pa=NULL;
2297         ParticleTexture ptex;
2298         float *cpa_fuv=0, *par_rot=0;
2299         float co[3], orco[3], ornor[3], hairmat[4][4], t, cpa_1st[3], dvec[3];
2300         float branch_begin, branch_end, branch_prob, rough_rand;
2301         float length, max_length = 1.0f, cur_length = 0.0f;
2302         float eff_length, eff_vec[3];
2303         int k, cpa_num;
2304         short cpa_from;
2305
2306         if(part->flag & PART_BRANCHING) {
2307                 branch_begin=rng_getFloat(thread->rng_path);
2308                 branch_end=branch_begin+(1.0f-branch_begin)*rng_getFloat(thread->rng_path);
2309                 branch_prob=rng_getFloat(thread->rng_path);
2310                 rough_rand=rng_getFloat(thread->rng_path);
2311         }
2312         else {
2313                 branch_begin= 0.0f;
2314                 branch_end= 0.0f;
2315                 branch_prob= 0.0f;
2316                 rough_rand= 0.0f;
2317         }
2318
2319         if(i<psys->totpart){
2320                 branch_begin=0.0f;
2321                 branch_end=1.0f;
2322                 branch_prob=0.0f;
2323         }
2324
2325         if(ctx->between){
2326                 int w, needupdate;
2327                 float foffset;
2328
2329                 if(ctx->editupdate && !(part->flag & PART_BRANCHING)) {
2330                         needupdate= 0;
2331                         w= 0;
2332                         while(w<4 && cpa->pa[w]>=0) {
2333                                 if(psys->edit->points[cpa->pa[w]].flag & PEP_EDIT_RECALC) {
2334                                         needupdate= 1;
2335                                         break;
2336                                 }
2337                                 w++;
2338                         }
2339
2340                         if(!needupdate)
2341                                 return;
2342                         else
2343                                 memset(keys, 0, sizeof(*keys)*(ctx->steps+1));
2344                 }
2345
2346                 /* get parent paths */
2347                 w= 0;
2348                 while(w<4 && cpa->pa[w]>=0){
2349                         key[w] = pcache[cpa->pa[w]];
2350                         w++;
2351                 }
2352
2353                 /* get the original coordinates (orco) for texture usage */
2354                 cpa_num = cpa->num;
2355                 
2356                 foffset= cpa->foffset;
2357                 cpa_fuv = cpa->fuv;
2358                 cpa_from = PART_FROM_FACE;
2359
2360                 psys_particle_on_emitter(ctx->sim.psmd,cpa_from,cpa_num,DMCACHE_ISCHILD,cpa->fuv,foffset,co,ornor,0,0,orco,0);
2361
2362                 if(part->path_start==0.0f) {
2363                         /* we need to save the actual root position of the child for positioning it accurately to the surface of the emitter */
2364                         VECCOPY(cpa_1st,co);
2365                         Mat4MulVecfl(ob->obmat,cpa_1st);
2366                 }
2367
2368                 pa = psys->particles + cpa->parent;
2369
2370                 psys_mat_hair_to_global(ob, ctx->sim.psmd->dm, psys->part->from, pa, hairmat);
2371
2372                 pa=0;
2373         }
2374         else{
2375                 if(ctx->editupdate && !(part->flag & PART_BRANCHING)) {
2376                         if(!(psys->edit->points[cpa->parent].flag & PEP_EDIT_RECALC))
2377                                 return;
2378
2379                         memset(keys, 0, sizeof(*keys)*(ctx->steps+1));
2380                 }
2381
2382                 /* get the parent path */
2383                 key[0]=pcache[cpa->parent];
2384
2385                 /* get the original coordinates (orco) for texture usage */
2386                 pa=psys->particles+cpa->parent;
2387
2388                 cpa_from=part->from;
2389                 cpa_num=pa->num;
2390                 cpa_fuv=pa->fuv;
2391
2392                 psys_particle_on_emitter(ctx->sim.psmd,cpa_from,cpa_num,DMCACHE_ISCHILD,cpa_fuv,pa->foffset,co,ornor,0,0,orco,0);
2393
2394                 psys_mat_hair_to_global(ob, ctx->sim.psmd->dm, psys->part->from, pa, hairmat);
2395         }
2396
2397         keys->steps = ctx->steps;
2398
2399         /* correct child ipo timing */
2400 #if 0 // XXX old animation system
2401         if((part->flag&PART_ABS_TIME)==0 && part->ipo){
2402                 float dsta=part->end-part->sta;
2403                 calc_ipo(part->ipo, 100.0f*(ctx->cfra-(part->sta+dsta*cpa->rand[1]))/(part->lifetime*(1.0f - part->randlife*cpa->rand[0])));
2404                 execute_ipo((ID *)part, part->ipo);
2405         }
2406 #endif // XXX old animation system
2407
2408         /* get different child parameters from textures & vgroups */
2409         get_child_modifier_parameters(part, ctx, cpa, cpa_from, cpa_num, cpa_fuv, orco, &ptex);
2410
2411         if(ptex.exist < PSYS_FRAND(i + 24)) {
2412                 keys->steps = -1;
2413                 return;
2414         }
2415
2416         /* create the child path */
2417         for(k=0,state=keys; k<=ctx->steps; k++,state++){
2418                 if(ctx->between){
2419                         int w=0;
2420
2421                         state->co[0] = state->co[1] = state->co[2] = 0.0f;
2422                         state->vel[0] = state->vel[1] = state->vel[2] = 0.0f;
2423                         state->rot[0] = state->rot[1] = state->rot[2] = state->rot[3] = 0.0f;
2424
2425                         //QUATCOPY(state->rot,key[0]->rot);
2426
2427                         /* child position is the weighted sum of parent positions */
2428                         while(w<4 && cpa->pa[w]>=0){
2429                                 state->co[0] += cpa->w[w] * key[w]->co[0];
2430                                 state->co[1] += cpa->w[w] * key[w]->co[1];
2431                                 state->co[2] += cpa->w[w] * key[w]->co[2];
2432
2433                                 state->vel[0] += cpa->w[w] * key[w]->vel[0];
2434                                 state->vel[1] += cpa->w[w] * key[w]->vel[1];
2435                                 state->vel[2] += cpa->w[w] * key[w]->vel[2];
2436                                 key[w]++;
2437                                 w++;
2438                         }
2439                         if(part->path_start==0.0f) {
2440                                 if(k==0){
2441                                         /* calculate the offset between actual child root position and first position interpolated from parents */
2442                                         VECSUB(cpa_1st,cpa_1st,state->co);
2443                                 }
2444                                 /* apply offset for correct positioning */
2445                                 VECADD(state->co,state->co,cpa_1st);
2446                         }
2447                 }
2448                 else{
2449                         /* offset the child from the parent position */
2450                         offset_child(cpa, (ParticleKey*)key[0], (ParticleKey*)state, part->childflat, part->childrad);
2451
2452                         key[0]++;
2453                 }
2454         }
2455
2456         /* apply effectors */
2457         if(part->flag & PART_CHILD_EFFECT) {
2458                 for(k=0,state=keys; k<=ctx->steps; k++,state++) {
2459                         if(k) {
2460                                 do_path_effectors(&ctx->sim, cpa->pa[0], state, k, ctx->steps, keys->co, ptex.effector, 0.0f, ctx->cfra, &eff_length, eff_vec);
2461                         }
2462                         else {
2463                                 VecSubf(eff_vec,(state+1)->co,state->co);
2464                                 eff_length= VecLength(eff_vec);
2465                         }
2466                 }
2467         }
2468
2469         for(k=0,state=keys; k<=ctx->steps; k++,state++){
2470                 t=(float)k/(float)ctx->steps;
2471
2472                 if(ctx->totparent){
2473                         if(i>=ctx->totparent) {
2474                                 /* this is now threadsafe, virtual parents are calculated before rest of children */
2475                                 par = cache[cpa->parent] + k;
2476                         }
2477                         else
2478                                 par=0;
2479                 }
2480                 else if(cpa->parent>=0){
2481                         par=pcache[cpa->parent]+k;
2482                         par_rot = par->rot;
2483                 }
2484
2485                 /* apply different deformations to the child path */
2486                 do_child_modifiers(&ctx->sim, &ptex, (ParticleKey *)par, par_rot, cpa, orco, hairmat, (ParticleKey *)state, t);
2487
2488                 /* TODO: better branching */
2489                 //if(part->flag & PART_BRANCHING && ctx->between == 0 && part->flag & PART_ANIM_BRANCHING)
2490                 //      rough_t = t * rough_rand;
2491                 //else
2492                 //      rough_t = t;
2493
2494                 /* TODO: better branching */
2495                 //if(part->flag & PART_BRANCHING && ctx->between==0){
2496                 //      if(branch_prob > part->branch_thres){
2497                 //              branchfac=0.0f;
2498                 //      }
2499                 //      else{
2500                 //              if(part->flag & PART_SYMM_BRANCHING){
2501                 //                      if(t < branch_begin || t > branch_end)
2502                 //                              branchfac=0.0f;
2503                 //                      else{
2504                 //                              if((t-branch_begin)/(branch_end-branch_begin)<0.5)
2505                 //                                      branchfac=2.0f*(t-branch_begin)/(branch_end-branch_begin);
2506                 //                              else
2507                 //                                      branchfac=2.0f*(branch_end-t)/(branch_end-branch_begin);
2508
2509                 //                              CLAMP(branchfac,0.0f,1.0f);
2510                 //                      }
2511                 //              }
2512                 //              else{
2513                 //                      if(t < branch_begin){
2514                 //                              branchfac=0.0f;
2515                 //                      }
2516                 //                      else{
2517                 //                              branchfac=(t-branch_begin)/((1.0f-branch_begin)*0.5f);
2518                 //                              CLAMP(branchfac,0.0f,1.0f);
2519                 //                      }
2520                 //              }
2521                 //      }
2522
2523                 //      if(i<psys->totpart)
2524                 //              VecLerpf(state->co, (pcache[i] + k)->co, state->co, branchfac);
2525                 //      else
2526                 //              /* this is not threadsafe, but should only happen for
2527                 //               * branching particles particles, which are not threaded */
2528                 //              VecLerpf(state->co, (cache[i - psys->totpart] + k)->co, state->co, branchfac);
2529                 //}
2530
2531                 /* we have to correct velocity because of kink & clump */
2532                 if(k>1){
2533                         VECSUB((state-1)->vel,state->co,(state-2)->co);
2534                         VecMulf((state-1)->vel,0.5);
2535
2536                         if(ctx->ma && (part->draw & PART_DRAW_MAT_COL))
2537                                 get_strand_normal(ctx->ma, ornor, cur_length, (state-1)->vel);
2538                 }
2539
2540                 /* check if path needs to be cut before actual end of data points */
2541                 if(k){
2542                         VECSUB(dvec,state->co,(state-1)->co);
2543                         length=1.0f/(float)ctx->steps;
2544                         k=check_path_length(k,keys,state,max_length,&cur_length,length,dvec);
2545                 }
2546                 else{
2547                         /* initialize length calculation */
2548                         max_length= ptex.length;
2549                         cur_length= 0.0f;
2550                 }
2551
2552                 if(ctx->ma && (part->draw & PART_DRAW_MAT_COL)) {
2553                         VECCOPY(state->col, &ctx->ma->r)
2554                         get_strand_normal(ctx->ma, ornor, cur_length, state->vel);
2555                 }
2556         }
2557 }
2558
2559 static void *exec_child_path_cache(void *data)
2560 {
2561         ParticleThread *thread= (ParticleThread*)data;
2562         ParticleThreadContext *ctx= thread->ctx;
2563         ParticleSystem *psys= ctx->sim.psys;
2564         ParticleCacheKey **cache= psys->childcache;
2565         ChildParticle *cpa;
2566         int i, totchild= ctx->totchild, first= 0;
2567
2568         if(thread->tot > 1){
2569                 first= ctx->parent_pass? 0 : ctx->totparent;
2570                 totchild= ctx->parent_pass? ctx->totparent : ctx->totchild;
2571         }
2572         
2573         cpa= psys->child + first + thread->num;
2574         for(i=first+thread->num; i<totchild; i+=thread->tot, cpa+=thread->tot)
2575                 psys_thread_create_path(thread, cpa, cache[i], i);
2576
2577         return 0;
2578 }
2579
2580 void psys_cache_child_paths(ParticleSimulationData *sim, float cfra, int editupdate)
2581 {
2582         ParticleSettings *part = sim->psys->part;
2583         ParticleThread *pthreads;
2584         ParticleThreadContext *ctx;
2585         ParticleCacheKey **cache;
2586         ListBase threads;
2587         int i, totchild, totparent, totthread;
2588
2589         if(sim->psys->flag & PSYS_GLOBAL_HAIR)
2590                 return;
2591
2592         pthreads= psys_threads_create(sim);
2593
2594         if(!psys_threads_init_path(pthreads, sim->scene, cfra, editupdate)) {
2595                 psys_threads_free(pthreads);
2596                 return;
2597         }
2598
2599         ctx= pthreads[0].ctx;
2600         totchild= ctx->totchild;
2601         totparent= ctx->totparent;
2602
2603         if(editupdate && sim->psys->childcache && !(part->flag & PART_BRANCHING) && totchild == sim->psys->totchildcache) {
2604                 cache = sim->psys->childcache;
2605         }
2606         else {
2607                 /* clear out old and create new empty path cache */
2608                 free_child_path_cache(sim->psys);
2609                 sim->psys->childcache= psys_alloc_path_cache_buffers(&sim->psys->childcachebufs, totchild, ctx->steps+1);
2610                 sim->psys->totchildcache = totchild;
2611         }
2612
2613         totthread= pthreads[0].tot;
2614
2615         if(totthread > 1) {
2616
2617                 /* make virtual child parents thread safe by calculating them first */
2618                 if(totparent) {
2619                         BLI_init_threads(&threads, exec_child_path_cache, totthread);
2620                         
2621                         for(i=0; i<totthread; i++) {
2622                                 pthreads[i].ctx->parent_pass = 1;
2623                                 BLI_insert_thread(&threads, &pthreads[i]);
2624                         }
2625
2626                         BLI_end_threads(&threads);
2627
2628                         for(i=0; i<totthread; i++)
2629                                 pthreads[i].ctx->parent_pass = 0;
2630                 }
2631
2632                 BLI_init_threads(&threads, exec_child_path_cache, totthread);
2633
2634                 for(i=0; i<totthread; i++)
2635                         BLI_insert_thread(&threads, &pthreads[i]);
2636
2637                 BLI_end_threads(&threads);
2638         }
2639         else
2640                 exec_child_path_cache(&pthreads[0]);
2641
2642         psys_threads_free(pthreads);
2643 }
2644 /* Calculates paths ready for drawing/rendering.                                                                        */
2645 /* -Usefull for making use of opengl vertex arrays for super fast strand drawing.       */
2646 /* -Makes child strands possible and creates them too into the cache.                           */
2647 /* -Cached path data is also used to determine cut position for the editmode tool.      */
2648 void psys_cache_paths(ParticleSimulationData *sim, float cfra)
2649 {
2650         PARTICLE_PSMD;
2651         ParticleEditSettings *pset = &sim->scene->toolsettings->particle;
2652         ParticleSystem *psys = sim->psys;
2653         ParticleSettings *part = psys->part;
2654         ParticleCacheKey *ca, **cache= psys->pathcache;
2655
2656         DerivedMesh *hair_dm = psys->hair_out_dm;
2657         
2658         ParticleKey result;
2659         
2660         Material *ma;
2661         ParticleInterpolationData pind;
2662
2663         PARTICLE_P;
2664         
2665         float birthtime = 0.0, dietime = 0.0;
2666         float t, time = 0.0, dfra = 1.0, frs_sec = sim->scene->r.frs_sec;
2667         float col[4] = {0.5f, 0.5f, 0.5f, 1.0f};
2668         float prev_tangent[3], hairmat[4][4];
2669         float rotmat[3][3];
2670         int k;
2671         int steps = (int)pow(2.0, (double)(psys->renderdata ? part->ren_step : part->draw_step));
2672         int totpart = psys->totpart;
2673         float length, vec[3];
2674         float *vg_effector= NULL, effector=0.0f;
2675         float *vg_length= NULL, pa_length=1.0f;
2676         int keyed, baked;
2677
2678         /* we don't have anything valid to create paths from so let's quit here */
2679         if((psys->flag & PSYS_HAIR_DONE || psys->flag & PSYS_KEYED || psys->pointcache->flag & PTCACHE_BAKED)==0)
2680                 return;
2681
2682         if(psys_in_edit_mode(sim->scene, psys))
2683                 if(psys->renderdata==0 && (psys->edit==NULL || pset->flag & PE_DRAW_PART)==0)
2684                         return;
2685         
2686         BLI_srandom(psys->seed);
2687
2688         keyed = psys->flag & PSYS_KEYED;
2689         baked = !hair_dm && psys->pointcache->flag & PTCACHE_BAKED;
2690
2691         /* clear out old and create new empty path cache */
2692         psys_free_path_cache(psys, psys->edit);
2693         cache= psys->pathcache= psys_alloc_path_cache_buffers(&psys->pathcachebufs, totpart, steps+1);
2694
2695         psys->lattice = psys_get_lattice(sim);
2696         ma= give_current_material(sim->ob, psys->part->omat);
2697         if(ma && (psys->part->draw & PART_DRAW_MAT_COL))
2698                 VECCOPY(col, &ma->r)
2699
2700         if(psys->part->from!=PART_FROM_PARTICLE && !(psys->flag & PSYS_GLOBAL_HAIR)) {
2701                 if(!(psys->part->flag & PART_CHILD_EFFECT))
2702                         vg_effector = psys_cache_vgroup(psmd->dm, psys, PSYS_VG_EFFECTOR);
2703                 
2704                 if(!psys->totchild)
2705                         vg_length = psys_cache_vgroup(psmd->dm, psys, PSYS_VG_LENGTH);
2706         }
2707
2708         /*---first main loop: create all actual particles' paths---*/
2709         LOOP_SHOWN_PARTICLES {
2710                 if(!psys->totchild) {
2711                         BLI_srandom(psys->seed + p);
2712                         pa_length = 1.0f - part->randlength * BLI_frand();
2713                         if(vg_length)
2714                                 pa_length *= psys_particle_value_from_verts(psmd->dm,part->from,pa,vg_length);
2715                 }
2716
2717                 pind.keyed = keyed;
2718                 pind.cache = baked ? psys->pointcache : NULL;
2719                 pind.epoint = NULL;
2720                 pind.bspline = (psys->part->flag & PART_HAIR_BSPLINE);
2721                 pind.dm = hair_dm;
2722
2723                 memset(cache[p], 0, sizeof(*cache[p])*(steps+1));
2724
2725                 cache[p]->steps = steps;
2726
2727                 /*--get the first data points--*/
2728                 init_particle_interpolation(sim->ob, sim->psys, pa, &pind);
2729
2730                 /* hairmat is needed for for non-hair particle too so we get proper rotations */
2731                 psys_mat_hair_to_global(sim->ob, psmd->dm, psys->part->from, pa, hairmat);
2732                 VECCOPY(rotmat[0], hairmat[2]);
2733                 VECCOPY(rotmat[1], hairmat[1]);
2734                 VECCOPY(rotmat[2], hairmat[0]);
2735
2736                 if(part->draw & PART_ABS_PATH_TIME) {
2737                         birthtime = MAX2(pind.birthtime, part->path_start);
2738                         dietime = MIN2(pind.dietime, part->path_end);
2739                 }
2740                 else {
2741                         float tb = pind.birthtime;
2742                         birthtime = tb + part->path_start * (pind.dietime - tb);
2743                         dietime = tb + part->path_end * (pind.dietime - tb);
2744                 }
2745
2746                 if(birthtime >= dietime) {
2747                         cache[p]->steps = -1;
2748                         continue;
2749                 }
2750
2751                 dietime = birthtime + pa_length * (dietime - birthtime);
2752
2753                 /*--interpolate actual path from data points--*/
2754                 for(k=0, ca=cache[p]; k<=steps; k++, ca++){
2755                         time = (float)k / (float)steps;
2756
2757                         t = birthtime + time * (dietime - birthtime);
2758
2759                         result.time = -t;
2760
2761                         do_particle_interpolation(psys, p, pa, t, frs_sec, &pind, &result);
2762
2763                         /* dynamic hair is in object space */
2764                         /* keyed and baked are allready in global space */
2765                         if(hair_dm)
2766                                 Mat4MulVecfl(sim->ob->obmat, result.co);
2767                         else if(!keyed && !baked && !(psys->flag & PSYS_GLOBAL_HAIR))
2768                                 Mat4MulVecfl(hairmat, result.co);
2769
2770                         VECCOPY(ca->co, result.co);
2771                         VECCOPY(ca->col, col);
2772                 }
2773                 
2774                 /*--modify paths and calculate rotation & velocity--*/
2775
2776                 VecSubf(vec,(cache[p]+1)->co,cache[p]->co);
2777                 length = VecLength(vec);
2778
2779                 effector= 1.0f;
2780                 if(vg_effector)
2781                         effector*= psys_particle_value_from_verts(psmd->dm,psys->part->from,pa,vg_effector);
2782
2783                 for(k=0, ca=cache[p]; k<=steps; k++, ca++) {
2784                         if(!(psys->flag & PSYS_GLOBAL_HAIR)) {
2785                         /* apply effectors */
2786                                 if(!(psys->part->flag & PART_CHILD_EFFECT) && k)
2787                                         do_path_effectors(sim, p, ca, k, steps, cache[p]->co, effector, dfra, cfra, &length, vec);
2788
2789                                 /* apply guide curves to path data */
2790                                 if(sim->psys->effectors && (psys->part->flag & PART_CHILD_EFFECT)==0)
2791                                         /* ca is safe to cast, since only co and vel are used */
2792                                         do_guides(sim->psys->effectors, (ParticleKey*)ca, p, (float)k/(float)steps);
2793
2794                                 /* apply lattice */
2795                                 if(psys->lattice)
2796                                         calc_latt_deform(psys->lattice, ca->co, 1.0f);
2797
2798                                 /* figure out rotation */
2799                                 
2800                                 if(k) {
2801                                         float cosangle, angle, tangent[3], normal[3], q[4];
2802
2803                                         if(k == 1) {
2804                                                 /* calculate initial tangent for incremental rotations */
2805                                                 VECSUB(tangent, ca->co, (ca - 1)->co);
2806                                                 VECCOPY(prev_tangent, tangent);
2807                                                 Normalize(prev_tangent);
2808
2809                                                 /* First rotation is based on emitting face orientation.                */
2810                                                 /* This is way better than having flipping rotations resulting  */
2811                                                 /* from using a global axis as a rotation pole (vec_to_quat()). */
2812                                                 /* It's not an ideal solution though since it disregards the    */
2813                                                 /* initial tangent, but taking that in to account will allow    */
2814                                                 /* the possibility of flipping again. -jahka                                    */
2815                                                 Mat3ToQuat_is_ok(rotmat, (ca-1)->rot);
2816                                         }
2817                                         else {
2818                                                 VECSUB(tangent, ca->co, (ca - 1)->co);
2819                                                 Normalize(tangent);
2820
2821                                                 cosangle= Inpf(tangent, prev_tangent);
2822
2823                                                 /* note we do the comparison on cosangle instead of
2824                                                 * angle, since floating point accuracy makes it give
2825                                                 * different results across platforms */
2826                                                 if(cosangle > 0.999999f) {
2827                                                         QUATCOPY((ca - 1)->rot, (ca - 2)->rot);
2828                                                 }
2829                                                 else {
2830                                                         angle= saacos(cosangle);
2831                                                         Crossf(normal, prev_tangent, tangent);
2832                                                         VecRotToQuat(normal, angle, q);
2833                                                         QuatMul((ca - 1)->rot, q, (ca - 2)->rot);
2834                                                 }
2835
2836                                                 VECCOPY(prev_tangent, tangent);
2837                                         }
2838
2839                                         if(k == steps)
2840                                                 QUATCOPY(ca->rot, (ca - 1)->rot);
2841                                 }
2842
2843                         }
2844                         
2845                         /* set velocity */
2846
2847                         if(k){
2848                                 VECSUB(ca->vel, ca->co, (ca-1)->co);
2849
2850                                 if(k==1) {
2851                                         VECCOPY((ca-1)->vel, ca->vel);
2852                                 }
2853
2854                         }
2855                 }
2856         }
2857
2858         psys->totcached = totpart;
2859
2860         if(psys && psys->lattice){
2861                 end_latt_deform(psys->lattice);
2862                 psys->lattice= NULL;
2863         }
2864
2865         if(vg_effector)
2866                 MEM_freeN(vg_effector);
2867
2868         if(vg_length)
2869                 MEM_freeN(vg_length);
2870 }
2871 void psys_cache_edit_paths(Scene *scene, Object *ob, PTCacheEdit *edit, float cfra)
2872 {
2873         ParticleCacheKey *ca, **cache= edit->pathcache;
2874         ParticleEditSettings *pset = &scene->toolsettings->particle;
2875         
2876         PTCacheEditPoint *point = edit->points;
2877         PTCacheEditKey *ekey = NULL;
2878
2879         ParticleSystem *psys = edit->psys;
2880         ParticleSystemModifierData *psmd = psys_get_modifier(ob, psys);
2881         ParticleData *pa = psys ? psys->particles : NULL;
2882
2883         ParticleInterpolationData pind;
2884         ParticleKey result;
2885         
2886         float birthtime = 0.0, dietime = 0.0;
2887         float t, time = 0.0, keytime = 0.0, frs_sec;
2888         float hairmat[4][4];
2889         int k,i;
2890         int steps = (int)pow(2.0, (double)pset->draw_step);
2891         int totpart = edit->totpoint;
2892         float sel_col[3];
2893         float nosel_col[3];
2894
2895         steps = MAX2(steps, 4);
2896
2897         if(!cache || edit->totpoint != edit->totcached) {
2898                 /* clear out old and create new empty path cache */
2899                 psys_free_path_cache(edit->psys, edit);
2900                 cache= edit->pathcache= psys_alloc_path_cache_buffers(&edit->pathcachebufs, totpart, steps+1);
2901         }
2902
2903         frs_sec = (psys || edit->pid.flag & PTCACHE_VEL_PER_SEC) ? 25.0f : 1.0f;
2904
2905         sel_col[0] = (float)edit->sel_col[0] / 255.0f;
2906         sel_col[1] = (float)edit->sel_col[1] / 255.0f;
2907         sel_col[2] = (float)edit->sel_col[2] / 255.0f;
2908         nosel_col[0] = (float)edit->nosel_col[0] / 255.0f;
2909         nosel_col[1] = (float)edit->nosel_col[1] / 255.0f;
2910         nosel_col[2] = (float)edit->nosel_col[2] / 255.0f;
2911
2912         /*---first main loop: create all actual particles' paths---*/
2913         for(i=0; i<totpart; i++, pa+=pa?1:0, point++){
2914                 if(edit->totcached && !(point->flag & PEP_EDIT_RECALC))
2915                         continue;
2916
2917                 ekey = point->keys;
2918
2919                 pind.keyed = 0;
2920                 pind.cache = NULL;
2921                 pind.epoint = point;
2922                 pind.bspline = psys ? (psys->part->flag & PART_HAIR_BSPLINE) : 0;
2923                 pind.dm = NULL;
2924
2925                 memset(cache[i], 0, sizeof(*cache[i])*(steps+1));
2926
2927                 cache[i]->steps = steps;
2928
2929                 /*--get the first data points--*/
2930                 init_particle_interpolation(ob, psys, pa, &pind);
2931
2932                 if(psys)
2933                         psys_mat_hair_to_global(ob, psmd->dm, psys->part->from, pa, hairmat);
2934
2935                 birthtime = pind.birthtime;
2936                 dietime = pind.dietime;
2937
2938                 if(birthtime >= dietime) {
2939                         cache[i]->steps = -1;
2940                         continue;
2941                 }
2942
2943                 /*--interpolate actual path from data points--*/
2944                 for(k=0, ca=cache[i]; k<=steps; k++, ca++){
2945                         time = (float)k / (float)steps;
2946
2947                         t = birthtime + time * (dietime - birthtime);
2948
2949                         result.time = -t;
2950
2951                         do_particle_interpolation(psys, i, pa, t, frs_sec, &pind, &result);
2952
2953                          /* non-hair points are allready in global space */
2954                         if(psys && !(psys->flag & PSYS_GLOBAL_HAIR))
2955                                 Mat4MulVecfl(hairmat, result.co);
2956
2957                         VECCOPY(ca->co, result.co);
2958
2959                         ca->vel[0] = ca->vel[1] = 0.0f;
2960                         ca->vel[1] = 1.0f;
2961
2962                         /* selection coloring in edit mode */
2963                         if((ekey + (pind.ekey[0] - point->keys))->flag & PEK_SELECT){
2964                                 if((ekey + (pind.ekey[1] - point->keys))->flag & PEK_SELECT){
2965                                         VECCOPY(ca->col, sel_col);
2966                                 }
2967                                 else{
2968                                         keytime = (t - (*pind.ekey[0]->time))/((*pind.ekey[1]->time) - (*pind.ekey[0]->time));
2969                                         VecLerpf(ca->col, sel_col, nosel_col, keytime);
2970                                 }
2971                         }
2972                         else{
2973                                 if((ekey + (pind.ekey[1] - point->keys))->flag & PEK_SELECT){
2974                                         keytime = (t - (*pind.ekey[0]->time))/((*pind.ekey[1]->time) - (*pind.ekey[0]->time));
2975                                         VecLerpf(ca->col, nosel_col, sel_col, keytime);
2976                                 }
2977                                 else{
2978                                         VECCOPY(ca->col, nosel_col);
2979                                 }
2980                         }
2981
2982                         ca->time = t;
2983                 }
2984         }
2985
2986         edit->totcached = totpart;
2987 }
2988 /************************************************/
2989 /*                      Particle Key handling                           */
2990 /************************************************/
2991 void copy_particle_key(ParticleKey *to, ParticleKey *from, int time){
2992         if(time){
2993                 memcpy(to,from,sizeof(ParticleKey));
2994         }
2995         else{
2996                 float to_time=to->time;
2997                 memcpy(to,from,sizeof(ParticleKey));
2998                 to->time=to_time;
2999         }
3000 }
3001 void psys_get_from_key(ParticleKey *key, float *loc, float *vel, float *rot, float *time){
3002         if(loc) VECCOPY(loc,key->co);
3003         if(vel) VECCOPY(vel,key->vel);
3004         if(rot) QUATCOPY(rot,key->rot);
3005         if(time) *time=key->time;
3006 }
3007 /*-------changing particle keys from space to another-------*/
3008 #if 0
3009 static void key_from_object(Object *ob, ParticleKey *key){
3010         float q[4];
3011
3012         VECADD(key->vel,key->vel,key->co);
3013
3014         Mat4MulVecfl(ob->obmat,key->co);
3015         Mat4MulVecfl(ob->obmat,key->vel);
3016         Mat4ToQuat(ob->obmat,q);
3017
3018         VECSUB(key->vel,key->vel,key->co);
3019         QuatMul(key->rot,q,key->rot);
3020 }
3021 #endif
3022
3023 static void triatomat(float *v1, float *v2, float *v3, float (*uv)[2], float mat[][4])
3024 {
3025         float det, w1, w2, d1[2], d2[2];
3026
3027         memset(mat, 0, sizeof(float)*4*4);
3028         mat[3][3]= 1.0f;
3029
3030         /* first axis is the normal */
3031         CalcNormFloat(v1, v2, v3, mat[2]);
3032
3033         /* second axis along (1, 0) in uv space */
3034         if(uv) {
3035                 d1[0]= uv[1][0] - uv[0][0];
3036                 d1[1]= uv[1][1] - uv[0][1];
3037                 d2[0]= uv[2][0] - uv[0][0];
3038                 d2[1]= uv[2][1] - uv[0][1];
3039
3040                 det = d2[0]*d1[1] - d2[1]*d1[0];
3041
3042                 if(det != 0.0f) {
3043                         det= 1.0f/det;
3044                         w1= -d2[1]*det;
3045                         w2= d1[1]*det;
3046
3047                         mat[1][0]= w1*(v2[0] - v1[0]) + w2*(v3[0] - v1[0]);
3048                         mat[1][1]= w1*(v2[1] - v1[1]) + w2*(v3[1] - v1[1]);
3049                         mat[1][2]= w1*(v2[2] - v1[2]) + w2*(v3[2] - v1[2]);
3050                         Normalize(mat[1]);
3051                 }
3052                 else
3053                         mat[1][0]= mat[1][1]= mat[1][2]= 0.0f;
3054         }
3055         else {
3056                 VecSubf(mat[1], v2, v1);
3057                 Normalize(mat[1]);
3058         }
3059         
3060         /* third as a cross product */
3061         Crossf(mat[0], mat[1], mat[2]);
3062 }
3063
3064 static void psys_face_mat(Object *ob, DerivedMesh *dm, ParticleData *pa, float mat[][4], int orco)
3065 {
3066         float v[3][3];
3067         MFace *mface;
3068         OrigSpaceFace *osface;
3069         float (*orcodata)[3];
3070
3071         int i = pa->num_dmcache==DMCACHE_NOTFOUND ? pa->num : pa->num_dmcache;
3072         
3073         if (i==-1 || i >= dm->getNumFaces(dm)) { Mat4One(mat); return; }
3074
3075         mface=dm->getFaceData(dm,i,CD_MFACE);
3076         osface=dm->getFaceData(dm,i,CD_ORIGSPACE);
3077         
3078         if(orco && (orcodata=dm->getVertDataArray(dm, CD_ORCO))) {
3079                 VECCOPY(v[0], orcodata[mface->v1]);
3080                 VECCOPY(v[1], orcodata[mface->v2]);
3081                 VECCOPY(v[2], orcodata[mface->v3]);
3082
3083                 /* ugly hack to use non-transformed orcos, since only those
3084                  * give symmetric results for mirroring in particle mode */
3085                 transform_mesh_orco_verts(ob->data, v, 3, 1);
3086         }
3087         else {
3088                 dm->getVertCo(dm,mface->v1,v[0]);