merge with 2.5 (not trunk, last merge message said that on accident) at r22252
[blender.git] / source / blender / blenkernel / intern / smoke.c
1 /**
2  * BKE_cloth.h
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) Blender Foundation.
23  * All rights reserved.
24  *
25  * The Original Code is: all of this file.
26  *
27  * Contributor(s): Daniel Genrich
28  *
29  * ***** END GPL LICENSE BLOCK *****
30  */
31
32 /* Part of the code copied from elbeem fluid library, copyright by Nils Thuerey */
33
34 #include <GL/glew.h>
35
36 #include "MEM_guardedalloc.h"
37
38 #include <float.h>
39 #include <math.h>
40 #include "stdio.h"
41
42 #include "BLI_linklist.h"
43 #include "BLI_rand.h"
44 #include "BLI_jitter.h"
45 #include "BLI_blenlib.h"
46 #include "BLI_arithb.h"
47 #include "BLI_edgehash.h"
48 #include "BLI_kdtree.h"
49 #include "BLI_kdopbvh.h"
50
51 #include "BKE_bvhutils.h"
52 #include "BKE_cdderivedmesh.h"
53 #include "BKE_customdata.h"
54 #include "BKE_DerivedMesh.h"
55 #include "BKE_modifier.h"
56 #include "BKE_particle.h"
57 #include "BKE_utildefines.h"
58
59 #include "DNA_customdata_types.h"
60 #include "DNA_group_types.h"
61 #include "DNA_mesh_types.h"
62 #include "DNA_meshdata_types.h"
63 #include "DNA_modifier_types.h"
64 #include "DNA_object_types.h"
65 #include "DNA_particle_types.h"
66 #include "DNA_scene_types.h"
67 #include "DNA_smoke_types.h"
68
69 #include "smoke_API.h"
70
71 #include "BKE_smoke.h"
72
73 #ifdef _WIN32
74 #include <time.h>
75 #include <stdio.h>
76 #include <conio.h>
77 #include <windows.h>
78
79 static LARGE_INTEGER liFrequency;
80 static LARGE_INTEGER liStartTime;
81 static LARGE_INTEGER liCurrentTime;
82
83 static void tstart ( void )
84 {
85         QueryPerformanceFrequency ( &liFrequency );
86         QueryPerformanceCounter ( &liStartTime );
87 }
88 static void tend ( void )
89 {
90         QueryPerformanceCounter ( &liCurrentTime );
91 }
92 static double tval()
93 {
94         return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart)* (double)1000.0/(double)liFrequency.QuadPart ));
95 }
96 #else
97 #include <sys/time.h>
98 static struct timeval _tstart, _tend;
99 static struct timezone tz;
100 static void tstart ( void )
101 {
102         gettimeofday ( &_tstart, &tz );
103 }
104 static void tend ( void )
105 {
106         gettimeofday ( &_tend,&tz );
107 }
108 static double tval()
109 {
110         double t1, t2;
111         t1 = ( double ) _tstart.tv_sec*1000 + ( double ) _tstart.tv_usec/ ( 1000 );
112         t2 = ( double ) _tend.tv_sec*1000 + ( double ) _tend.tv_usec/ ( 1000 );
113         return t2-t1;
114 }
115 #endif
116
117 struct Object;
118 struct Scene;
119 struct DerivedMesh;
120 struct SmokeModifierData;
121
122 // forward declerations
123 static void get_cell(struct SmokeModifierData *smd, float *pos, int *cell, int correct);
124 static void get_bigcell(struct SmokeModifierData *smd, float *pos, int *cell, int correct);
125 void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *tris, int numfaces, int numtris, int **tridivs, float cell_len);
126
127 #define TRI_UVOFFSET (1./4.)
128
129
130 BVHTree *bvhtree_build_from_smoke ( float mat[4][4], MFace *mfaces, unsigned int numfaces, MVert *x, unsigned int numverts, float epsilon )
131 {
132         BVHTree *tree;
133         float co[12];
134         int i;
135         MFace *tface = mfaces;
136
137         // calc quads 
138         // todo
139
140         tree = BLI_bvhtree_new ( numfaces, epsilon, 2, 6 );
141
142         // fill tree
143         for ( i = 0; i < numfaces; i++, tface++ )
144         {
145                 VECCOPY ( &co[0*3], x[tface->v1].co );
146                 Mat4MulVecfl (mat, &co[0*3]);
147                 VECCOPY ( &co[1*3], x[tface->v2].co );
148                 Mat4MulVecfl (mat, &co[1*3]);
149                 VECCOPY ( &co[2*3], x[tface->v3].co );
150                 Mat4MulVecfl (mat, &co[2*3]);
151
152                 if ( tface->v4 )
153                 {
154                         VECCOPY ( &co[3*3], x[tface->v4].co );
155                         Mat4MulVecfl (mat, &co[3*3]);
156                 }
157
158                 BLI_bvhtree_insert ( tree, i, co, ( mfaces->v4 ? 4 : 3 ) );
159         }
160
161         // balance tree
162         BLI_bvhtree_balance ( tree );
163
164         return tree;
165 }
166
167 void bvhtree_update_from_smoke ( float mat[4][4], BVHTree * bvhtree, MFace *faces, int numfaces, MVert *x, int numverts)
168 {
169         int i;
170         MFace *mfaces = faces;
171         float co[12];
172         int ret = 0;
173
174         if ( !bvhtree )
175                 return;
176
177         if ( x )
178         {
179                 for ( i = 0; i < numfaces; i++, mfaces++ )
180                 {
181                         VECCOPY ( &co[0*3], x[mfaces->v1].co );
182                         Mat4MulVecfl (mat, &co[0*3]);
183                         VECCOPY ( &co[1*3], x[mfaces->v2].co );
184                         Mat4MulVecfl (mat, &co[1*3]);
185                         VECCOPY ( &co[2*3], x[mfaces->v3].co );
186                         Mat4MulVecfl (mat, &co[2*3]);
187
188                         if ( mfaces->v4 )
189                         {
190                                 VECCOPY ( &co[3*3], x[mfaces->v4].co );
191                                 Mat4MulVecfl (mat, &co[3*3]);
192                         }
193
194                         ret = BLI_bvhtree_update_node ( bvhtree, i, co, NULL, ( mfaces->v4 ? 4 : 3 ) );
195
196                         // check if tree is already full
197                         if ( !ret )
198                                 break;
199                 }
200
201                 BLI_bvhtree_update_tree ( bvhtree );
202         }
203 }
204
205 int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, DerivedMesh *dm)
206 {
207         if((smd->type & MOD_SMOKE_TYPE_DOMAIN) && smd->domain && !smd->domain->fluid)
208         {
209                 size_t i;
210                 float min[3] = {FLT_MAX, FLT_MAX, FLT_MAX}, max[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
211                 float size[3];
212                 MVert *verts = dm->getVertArray(dm);
213                 float scale = 0.0;
214                 int res;                
215
216                 res = smd->domain->maxres;
217
218                 // get BB of domain
219                 for(i = 0; i < dm->getNumVerts(dm); i++)
220                 {
221                         float tmp[3];
222
223                         VECCOPY(tmp, verts[i].co);
224                         Mat4MulVecfl(ob->obmat, tmp);
225
226                         // min BB
227                         min[0] = MIN2(min[0], tmp[0]);
228                         min[1] = MIN2(min[1], tmp[1]);
229                         min[2] = MIN2(min[2], tmp[2]);
230
231                         // max BB
232                         max[0] = MAX2(max[0], tmp[0]);
233                         max[1] = MAX2(max[1], tmp[1]);
234                         max[2] = MAX2(max[2], tmp[2]);
235                 }
236
237                 VECCOPY(smd->domain->p0, min);
238                 VECCOPY(smd->domain->p1, max);
239
240                 // calc other res with max_res provided
241                 VECSUB(size, max, min);
242                 if(size[0] > size[1])
243                 {
244                         if(size[0] > size[1])
245                         {
246                                 scale = res / size[0];
247                                 smd->domain->dx = size[0] / res;
248                                 smd->domain->res[0] = res;
249                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
250                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
251                         }
252                         else
253                         {
254                                 scale = res / size[1];
255                                 smd->domain->dx = size[1] / res;
256                                 smd->domain->res[1] = res;
257                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
258                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
259                         }
260                 }
261                 else
262                 {
263                         if(size[1] > size[2])
264                         {
265                                 scale = res / size[1];
266                                 smd->domain->dx = size[1] / res;
267                                 smd->domain->res[1] = res;
268                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
269                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
270                         }
271                         else
272                         {
273                                 scale = res / size[2];
274                                 smd->domain->dx = size[2] / res;
275                                 smd->domain->res[2] = res;
276                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
277                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
278                         }
279                 }
280
281                 // TODO: put in failsafe if res<=0 - dg
282
283                 // printf("res[0]: %d, res[1]: %d, res[2]: %d\n", smd->domain->res[0], smd->domain->res[1], smd->domain->res[2]);
284
285                 // dt max is 0.1
286                 smd->domain->fluid = smoke_init(smd->domain->res, smd->domain->amplify, smd->domain->p0, smd->domain->p1, 2.5 / FPS);
287                 smd->time = scene->r.cfra;
288                 smd->domain->firstframe = smd->time;
289                 
290                 smoke_initBlenderRNA(smd->domain->fluid, &(smd->domain->alpha), &(smd->domain->beta));
291
292                 return 1;
293         }
294         else if((smd->type & MOD_SMOKE_TYPE_FLOW) && smd->flow)
295         {
296                 // handle flow object here
297                 // XXX TODO
298
299                 smd->time = scene->r.cfra;
300
301                 // update particle lifetime to be one frame
302                 // smd->flow->psys->part->lifetime = scene->r.efra + 1;
303
304                 return 1;
305         }
306         else if((smd->type & MOD_SMOKE_TYPE_COLL))
307         {
308                 smd->time = scene->r.cfra;
309
310                 // todo: delete this when loading colls work -dg
311                 if(!smd->coll)
312                         smokeModifier_createType(smd);
313
314                 if(!smd->coll->points)
315                 {
316                         // init collision points
317                         SmokeCollSettings *scs = smd->coll;
318                         MVert *mvert = dm->getVertArray(dm);
319                         MFace *mface = dm->getTessFaceArray(dm);
320                         size_t i = 0, divs = 0;
321                         int *tridivs = NULL;
322                         float cell_len = 1.0 / 50.0; // for res = 50
323                         size_t newdivs = 0;
324                         size_t max_points = 0;
325                         size_t quads = 0, facecounter = 0;
326
327                         // count quads
328                         for(i = 0; i < dm->getNumFaces(dm); i++)
329                         {
330                                 if(mface[i].v4)
331                                         quads++;
332                         }
333
334                         calcTriangleDivs(ob, mvert, dm->getNumVerts(dm), mface,  dm->getNumFaces(dm), dm->getNumFaces(dm) + quads, &tridivs, cell_len);
335
336                         // count triangle divisions
337                         for(i = 0; i < dm->getNumFaces(dm) + quads; i++)
338                         {
339                                 divs += (tridivs[3 * i] + 1) * (tridivs[3 * i + 1] + 1) * (tridivs[3 * i + 2] + 1);
340                         }
341
342                         // printf("divs: %d\n", divs);
343
344                         scs->points = MEM_callocN(sizeof(float) * (dm->getNumVerts(dm) + divs) * 3, "SmokeCollPoints");
345
346                         for(i = 0; i < dm->getNumVerts(dm); i++)
347                         {
348                                 float tmpvec[3];
349                                 VECCOPY(tmpvec, mvert[i].co);
350                                 Mat4MulVecfl (ob->obmat, tmpvec);
351                                 VECCOPY(&scs->points[i * 3], tmpvec);
352                         }
353                         
354                         for(i = 0, facecounter = 0; i < dm->getNumFaces(dm); i++)
355                         {
356                                 int again = 0;
357                                 do
358                                 {
359                                         size_t j, k;
360                                         int divs1 = tridivs[3 * facecounter + 0];
361                                         int divs2 = tridivs[3 * facecounter + 1];
362                                         int divs3 = tridivs[3 * facecounter + 2];
363                                         float side1[3], side2[3], trinormorg[3], trinorm[3];
364                                         
365                                         if(again == 1 && mface[i].v4)
366                                         {
367                                                 VECSUB(side1,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
368                                                 VECSUB(side2,  mvert[ mface[i].v4 ].co, mvert[ mface[i].v1 ].co);
369                                         }
370                                         else
371                                         {
372                                                 VECSUB(side1,  mvert[ mface[i].v2 ].co, mvert[ mface[i].v1 ].co);
373                                                 VECSUB(side2,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
374                                         }
375
376                                         Crossf(trinormorg, side1, side2);
377                                         Normalize(trinormorg);
378                                         VECCOPY(trinorm, trinormorg);
379                                         VecMulf(trinorm, 0.25 * cell_len);
380
381                                         for(j = 0; j <= divs1; j++)
382                                         {
383                                                 for(k = 0; k <= divs2; k++)
384                                                 {
385                                                         float p1[3], p2[3], p3[3], p[3]={0,0,0}; 
386                                                         const float uf = (float)(j + TRI_UVOFFSET) / (float)(divs1 + 0.0);
387                                                         const float vf = (float)(k + TRI_UVOFFSET) / (float)(divs2 + 0.0);
388                                                         float tmpvec[3];
389                                                         
390                                                         if(uf+vf > 1.0) 
391                                                         {
392                                                                 // printf("bigger - divs1: %d, divs2: %d\n", divs1, divs2);
393                                                                 continue;
394                                                         }
395
396                                                         VECCOPY(p1, mvert[ mface[i].v1 ].co);
397                                                         if(again == 1 && mface[i].v4)
398                                                         {
399                                                                 VECCOPY(p2, mvert[ mface[i].v3 ].co);
400                                                                 VECCOPY(p3, mvert[ mface[i].v4 ].co);
401                                                         }
402                                                         else
403                                                         {
404                                                                 VECCOPY(p2, mvert[ mface[i].v2 ].co);
405                                                                 VECCOPY(p3, mvert[ mface[i].v3 ].co);
406                                                         }
407
408                                                         VecMulf(p1, (1.0-uf-vf));
409                                                         VecMulf(p2, uf);
410                                                         VecMulf(p3, vf);
411                                                         
412                                                         VECADD(p, p1, p2);
413                                                         VECADD(p, p, p3);
414
415                                                         if(newdivs > divs)
416                                                                 printf("mem problem\n");
417
418                                                         // mMovPoints.push_back(p + trinorm);
419                                                         VECCOPY(tmpvec, p);
420                                                         VECADD(tmpvec, tmpvec, trinorm);
421                                                         Mat4MulVecfl (ob->obmat, tmpvec);
422                                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
423                                                         newdivs++;
424
425                                                         if(newdivs > divs)
426                                                                 printf("mem problem\n");
427
428                                                         // mMovPoints.push_back(p - trinorm);
429                                                         VECCOPY(tmpvec, p);
430                                                         VECSUB(tmpvec, tmpvec, trinorm);
431                                                         Mat4MulVecfl (ob->obmat, tmpvec);
432                                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
433                                                         newdivs++;
434                                                 }
435                                         }
436
437                                         if(again == 0 && mface[i].v4)
438                                                 again++;
439                                         else
440                                                 again = 0;
441
442                                         facecounter++;
443
444                                 } while(again!=0);
445                         }
446
447                         scs->numpoints = dm->getNumVerts(dm) + newdivs;
448
449                         MEM_freeN(tridivs);
450                 }
451
452                 if(!smd->coll->bvhtree)
453                 {
454                         smd->coll->bvhtree = bvhtree_build_from_smoke ( ob->obmat, dm->getTessFaceArray(dm), dm->getNumFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm), 0.0 );
455                 }
456
457         }
458
459         return 0;
460 }
461
462 /*! init triangle divisions */
463 void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *faces, int numfaces, int numtris, int **tridivs, float cell_len) 
464 {
465         // mTriangleDivs1.resize( faces.size() );
466         // mTriangleDivs2.resize( faces.size() );
467         // mTriangleDivs3.resize( faces.size() );
468
469         size_t i = 0, facecounter = 0;
470         float maxscale[3] = {1,1,1}; // = channelFindMaxVf(mcScale);
471         float maxpart = ABS(maxscale[0]);
472         float scaleFac = 0;
473         float fsTri = 0;
474         if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
475         if(ABS(maxscale[2])>maxpart) maxpart = ABS(maxscale[2]);
476         scaleFac = 1.0 / maxpart;
477         // featureSize = mLevel[mMaxRefine].nodeSize
478         fsTri = cell_len * 0.5 * scaleFac;
479
480         if(*tridivs)
481                 MEM_freeN(*tridivs);
482
483         *tridivs = MEM_callocN(sizeof(int) * numtris * 3, "Smoke_Tridivs");
484
485         for(i = 0, facecounter = 0; i < numfaces; i++) 
486         {
487                 float p0[3], p1[3], p2[3];
488                 float side1[3];
489                 float side2[3];
490                 float side3[3];
491                 int divs1=0, divs2=0, divs3=0;
492
493                 VECCOPY(p0, verts[faces[i].v1].co);
494                 Mat4MulVecfl (ob->obmat, p0);
495                 VECCOPY(p1, verts[faces[i].v2].co);
496                 Mat4MulVecfl (ob->obmat, p1);
497                 VECCOPY(p2, verts[faces[i].v3].co);
498                 Mat4MulVecfl (ob->obmat, p2);
499
500                 VECSUB(side1, p1, p0);
501                 VECSUB(side2, p2, p0);
502                 VECSUB(side3, p1, p2);
503
504                 if(INPR(side1, side1) > fsTri*fsTri) 
505                 { 
506                         float tmp = Normalize(side1);
507                         divs1 = (int)(tmp/fsTri); 
508                 }
509                 if(INPR(side2, side2) > fsTri*fsTri) 
510                 { 
511                         float tmp = Normalize(side2);
512                         divs2 = (int)(tmp/fsTri); 
513                         
514                         /*
515                         // debug
516                         if(i==0)
517                                 printf("b tmp: %f, fsTri: %f, divs2: %d\n", tmp, fsTri, divs2);
518                         */
519                 }
520
521                 (*tridivs)[3 * facecounter + 0] = divs1;
522                 (*tridivs)[3 * facecounter + 1] = divs2;
523                 (*tridivs)[3 * facecounter + 2] = divs3;
524
525                 // TODO quad case
526                 if(faces[i].v4)
527                 {
528                         divs1=0, divs2=0, divs3=0;
529
530                         facecounter++;
531                         
532                         VECCOPY(p1, verts[faces[i].v3].co);
533                         Mat4MulVecfl (ob->obmat, p1);
534                         VECCOPY(p2, verts[faces[i].v4].co);
535                         Mat4MulVecfl (ob->obmat, p2);
536
537                         VECSUB(side1, p1, p0);
538                         VECSUB(side2, p2, p0);
539                         VECSUB(side3, p1, p2);
540
541                         if(INPR(side1, side1) > fsTri*fsTri) 
542                         { 
543                                 float tmp = Normalize(side1);
544                                 divs1 = (int)(tmp/fsTri); 
545                         }
546                         if(INPR(side2, side2) > fsTri*fsTri) 
547                         { 
548                                 float tmp = Normalize(side2);
549                                 divs2 = (int)(tmp/fsTri); 
550                         }
551
552                         (*tridivs)[3 * facecounter + 0] = divs1;
553                         (*tridivs)[3 * facecounter + 1] = divs2;
554                         (*tridivs)[3 * facecounter + 2] = divs3;
555                 }
556                 facecounter++;
557         }
558 }
559
560 void smokeModifier_freeDomain(SmokeModifierData *smd)
561 {
562         if(smd->domain)
563         {
564                 // free visualisation buffers
565                 if(smd->domain->bind)
566                 {
567                         glDeleteTextures(smd->domain->max_textures, (GLuint *)smd->domain->bind);
568                         MEM_freeN(smd->domain->bind);
569                 }
570                 smd->domain->max_textures = 0; // unnecessary but let's be sure
571
572                 if(smd->domain->tray)
573                         MEM_freeN(smd->domain->tray);
574                 if(smd->domain->tvox)
575                         MEM_freeN(smd->domain->tvox);
576                 if(smd->domain->traybig)
577                         MEM_freeN(smd->domain->traybig);
578                 if(smd->domain->tvoxbig)
579                         MEM_freeN(smd->domain->tvoxbig);
580
581                 if(smd->domain->fluid)
582                 {
583                         smoke_free(smd->domain->fluid);
584                 }
585                 MEM_freeN(smd->domain);
586                 smd->domain = NULL;
587         }
588 }
589
590 void smokeModifier_freeFlow(SmokeModifierData *smd)
591 {
592         if(smd->flow)
593         {
594                 MEM_freeN(smd->flow);
595                 smd->flow = NULL;
596         }
597 }
598
599 void smokeModifier_freeCollision(SmokeModifierData *smd)
600 {
601         if(smd->coll)
602         {
603                 if(smd->coll->points)
604                 {
605                         MEM_freeN(smd->coll->points);
606                         smd->coll->points = NULL;
607                 }
608
609                 if(smd->coll->bvhtree)
610                 {
611                         BLI_bvhtree_free(smd->coll->bvhtree);
612                         smd->coll->bvhtree = NULL;
613                 }
614
615                 MEM_freeN(smd->coll);
616                 smd->coll = NULL;
617         }
618 }
619
620 void smokeModifier_reset(struct SmokeModifierData *smd)
621 {
622         if(smd)
623         {
624                 if(smd->domain)
625                 {
626                         // free visualisation buffers
627                         if(smd->domain->bind)
628                         {
629                                 glDeleteTextures(smd->domain->max_textures, (GLuint *)smd->domain->bind);
630                                 MEM_freeN(smd->domain->bind);
631                                 smd->domain->bind = NULL;
632                         }
633                         smd->domain->max_textures = 0;
634                         smd->domain->viewsettings = 0; // reset view for new frame
635
636                         if(smd->domain->tray)
637                                 MEM_freeN(smd->domain->tray);
638                         if(smd->domain->tvox)
639                                 MEM_freeN(smd->domain->tvox);
640                         if(smd->domain->traybig)
641                                 MEM_freeN(smd->domain->traybig);
642                         if(smd->domain->tvoxbig)
643                                 MEM_freeN(smd->domain->tvoxbig);
644
645                         smd->domain->tvox = NULL;
646                         smd->domain->tray = NULL;
647                         smd->domain->tvoxbig = NULL;
648                         smd->domain->traybig = NULL;
649
650                         if(smd->domain->fluid)
651                         {
652                                 smoke_free(smd->domain->fluid);
653                                 smd->domain->fluid = NULL;
654                         }
655                 }
656                 else if(smd->flow)
657                 {
658                                                 
659                 }
660                 else if(smd->coll)
661                 {
662                         if(smd->coll->points)
663                         {
664                                 MEM_freeN(smd->coll->points);
665                                 smd->coll->points = NULL;
666                         }
667
668                         if(smd->coll->bvhtree)
669                         {
670                                 BLI_bvhtree_free(smd->coll->bvhtree);
671                                 smd->coll->bvhtree = NULL;
672                         }
673                 }
674         }
675 }
676
677 void smokeModifier_free (SmokeModifierData *smd)
678 {
679         if(smd)
680         {
681                 smokeModifier_freeDomain(smd);
682                 smokeModifier_freeFlow(smd);
683                 smokeModifier_freeCollision(smd);
684         }
685 }
686
687 void smokeModifier_createType(struct SmokeModifierData *smd)
688 {
689         if(smd)
690         {
691                 if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
692                 {
693                         if(smd->domain)
694                                 smokeModifier_freeDomain(smd);
695
696                         smd->domain = MEM_callocN(sizeof(SmokeDomainSettings), "SmokeDomain");
697
698                         smd->domain->smd = smd;
699
700                         /* set some standard values */
701                         smd->domain->fluid = NULL;
702                         smd->domain->eff_group = NULL;
703                         smd->domain->fluid_group = NULL;
704                         smd->domain->coll_group = NULL;
705                         smd->domain->maxres = 32;
706                         smd->domain->amplify = 2;
707                         smd->domain->omega = 1.0;
708                         smd->domain->alpha = -0.001;
709                         smd->domain->beta = 0.1;
710                         smd->domain->flags = 0;
711                         smd->domain->noise = MOD_SMOKE_NOISEWAVE;
712                         smd->domain->visibility = 1;
713
714                         // init 3dview buffer
715                         smd->domain->tvox = NULL;
716                         smd->domain->tray = NULL;
717                         smd->domain->tvoxbig = NULL;
718                         smd->domain->traybig = NULL;
719                         smd->domain->viewsettings = 0;
720                         smd->domain->bind = NULL;
721                         smd->domain->max_textures = 0;
722                 }
723                 else if(smd->type & MOD_SMOKE_TYPE_FLOW)
724                 {
725                         if(smd->flow)
726                                 smokeModifier_freeFlow(smd);
727
728                         smd->flow = MEM_callocN(sizeof(SmokeFlowSettings), "SmokeFlow");
729
730                         smd->flow->smd = smd;
731
732                         /* set some standard values */
733                         smd->flow->density = 1.0;
734                         smd->flow->temp = 1.0;
735
736                         smd->flow->psys = NULL;
737
738                 }
739                 else if(smd->type & MOD_SMOKE_TYPE_COLL)
740                 {
741                         if(smd->coll)
742                                 smokeModifier_freeCollision(smd);
743
744                         smd->coll = MEM_callocN(sizeof(SmokeCollSettings), "SmokeColl");
745
746                         smd->coll->smd = smd;
747                         smd->coll->points = NULL;
748                         smd->coll->numpoints = 0;
749                         smd->coll->bvhtree = NULL;
750                 }
751         }
752 }
753
754 // forward declaration
755 void smoke_calc_transparency(struct SmokeModifierData *smd, float *light, int big);
756
757 void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc)
758 {       
759         if(scene->r.cfra >= smd->time)
760                 smokeModifier_init(smd, ob, scene, dm);
761
762         if((smd->type & MOD_SMOKE_TYPE_FLOW))
763         {
764                 if(scene->r.cfra > smd->time)
765                 {
766                         // XXX TODO
767                         smd->time = scene->r.cfra;
768                 }
769                 else if(scene->r.cfra < smd->time)
770                 {
771                         smd->time = scene->r.cfra;
772                         smokeModifier_reset(smd);
773                 }
774         }
775         else if(smd->type & MOD_SMOKE_TYPE_COLL)
776         {
777                 if(scene->r.cfra > smd->time)
778                 {
779                         // XXX TODO
780                         smd->time = scene->r.cfra;
781                         
782                         if(smd->coll->bvhtree)
783                                 bvhtree_update_from_smoke ( ob->obmat, smd->coll->bvhtree, dm->getTessFaceArray(dm), dm->getNumFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm));
784                         else
785                                 printf("smoke coll with no bvh\n");
786                 }
787                 else if(scene->r.cfra < smd->time)
788                 {
789                         smd->time = scene->r.cfra;
790                         smokeModifier_reset(smd);
791                 }
792         }
793         else if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
794         {
795                 SmokeDomainSettings *sds = smd->domain;
796                 
797                 if(scene->r.cfra > smd->time)
798                 {
799                         GroupObject *go = NULL;
800                         Base *base = NULL;
801                         int cnt_domain = 0;
802                         
803                         tstart();
804
805                         sds->viewsettings = 0; // reset view for new frame
806
807                         // check for 2nd domain, if not there -> no groups are necessary
808                         for(base = scene->base.first; base; base= base->next) 
809                         {
810                                 Object *ob1= base->object;
811                                 
812                                 if(ob1 && ob1 != ob)
813                                 {
814                                         ModifierData *tmd = modifiers_findByType(ob1, eModifierType_Smoke);
815
816                                         if(tmd && tmd->mode & (eModifierMode_Realtime | eModifierMode_Render))
817                                         {
818                                                 SmokeModifierData *tsmd = (SmokeModifierData *)tmd;
819
820                                                 if((tsmd->type & MOD_SMOKE_TYPE_DOMAIN))
821                                                 {
822                                                         cnt_domain++;
823                                                 }
824                                         }
825                                 }
826                         }
827
828                         // do flows and fluids
829                         if(sds->fluid_group || !cnt_domain)
830                         {
831                                 Object *otherobj = NULL;
832                                 ModifierData *md = NULL;
833
834                                 if(cnt_domain && !sds->fluid_group) // we use groups since we have 2 domains
835                                         go = sds->fluid_group->gobject.first;
836                                 else
837                                         base = scene->base.first;
838
839                                 while(base || go)
840                                 {
841                                         otherobj = NULL;
842
843                                         if(cnt_domain && !sds->fluid_group) 
844                                         {
845                                                 if(go->ob)
846                                                         otherobj = go->ob;
847                                         }
848                                         else
849                                                 otherobj = base->object;
850
851                                         if(!otherobj)
852                                         {
853                                                 if(cnt_domain && !sds->fluid_group)
854                                                         go = go->next;
855                                                 else
856                                                         base= base->next;
857
858                                                 continue;
859                                         }
860
861                                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
862                                         
863                                         // check for active smoke modifier
864                                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
865                                         {
866                                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
867                                                 
868                                                 // check for initialized smoke object
869                                                 if((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow)
870                                                 {
871                                                         // we got nice flow object
872                                                         SmokeFlowSettings *sfs = smd2->flow;
873                                                         
874                                                         if(sfs->psys && sfs->psys->part && sfs->psys->part->type==PART_EMITTER) // is particle system selected
875                                                         {
876                                                                 ParticleSystem *psys = sfs->psys;
877                                                                 ParticleSettings *part=psys->part;
878                                                                 ParticleData *pa = NULL;
879                                                                 int p = 0;
880                                                                 float *density = smoke_get_density(sds->fluid);
881                                                                 float *bigdensity = smoke_get_bigdensity(sds->fluid);
882                                                                 float *heat = smoke_get_heat(sds->fluid);
883                                                                 float *velocity_x = smoke_get_velocity_x(sds->fluid);
884                                                                 float *velocity_y = smoke_get_velocity_y(sds->fluid);
885                                                                 float *velocity_z = smoke_get_velocity_z(sds->fluid);
886                                                                 int bigres[3];
887
888                                                                 smoke_get_bigres(smd->domain->fluid, bigres);
889                                                                 
890                                                                 // mostly copied from particle code
891                                                                 for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++)
892                                                                 {
893                                                                         int cell[3];
894                                                                         size_t i = 0;
895                                                                         size_t index = 0;
896                                                                         int badcell = 0;
897                                                                         
898                                                                         if(pa->alive == PARS_KILLED) continue;
899                                                                         else if(pa->alive == PARS_UNBORN && (part->flag & PART_UNBORN)==0) continue;
900                                                                         else if(pa->alive == PARS_DEAD && (part->flag & PART_DIED)==0) continue;
901                                                                         else if(pa->flag & (PARS_UNEXIST+PARS_NO_DISP)) continue;
902                                                                         
903                                                                         // VECCOPY(pos, pa->state.co);
904                                                                         // Mat4MulVecfl (ob->imat, pos);
905                                                                         
906                                                                         // 1. get corresponding cell
907                                                                         get_cell(smd, pa->state.co, cell, 0);
908                                                                 
909                                                                         // check if cell is valid (in the domain boundary)
910                                                                         for(i = 0; i < 3; i++)
911                                                                         {
912                                                                                 if((cell[i] > sds->res[i] - 1) || (cell[i] < 0))
913                                                                                 {
914                                                                                         badcell = 1;
915                                                                                         break;
916                                                                                 }
917                                                                         }
918                                                                                 
919                                                                         if(badcell)
920                                                                                 continue;
921                                                                         
922                                                                         // 2. set cell values (heat, density and velocity)
923                                                                         index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
924                                                                         
925                                                                         if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW)) // this is inflow
926                                                                         {
927                                                                                 heat[index] = sfs->temp;
928                                                                                 density[index] = sfs->density;
929                                                                                 velocity_x[index] = pa->state.vel[0];
930                                                                                 velocity_y[index] = pa->state.vel[1];
931                                                                                 velocity_z[index] = pa->state.vel[2];
932
933                                                                                 // we need different handling for the high-res feature
934                                                                                 if(bigdensity)
935                                                                                 {
936                                                                                         // init all surrounding cells according to amplification, too
937                                                                                         int i, j, k;
938                                                                                         for(i = 0; i < smd->domain->amplify; i++)
939                                                                                                 for(j = 0; j < smd->domain->amplify; j++)
940                                                                                                         for(k = 0; k < smd->domain->amplify; k++)
941                                                                                                         {
942                                                                                                                 index = smoke_get_index(smd->domain->amplify * cell[0] + i, bigres[0], smd->domain->amplify * cell[1] + j, bigres[1], smd->domain->amplify * cell[2] + k);
943                                                                                                                 bigdensity[index] = sfs->density;
944                                                                                                         }
945                                                                                 }
946                                                                         }
947                                                                         else // outflow
948                                                                         {
949                                                                                 heat[index] = 0.f;
950                                                                                 density[index] = 0.f;
951                                                                                 velocity_x[index] = 0.f;
952                                                                                 velocity_y[index] = 0.f;
953                                                                                 velocity_z[index] = 0.f;
954
955                                                                                 // we need different handling for the high-res feature
956                                                                                 if(bigdensity)
957                                                                                 {
958                                                                                         // init all surrounding cells according to amplification, too
959                                                                                         int i, j, k;
960                                                                                         for(i = 0; i < smd->domain->amplify; i++)
961                                                                                                 for(j = 0; j < smd->domain->amplify; j++)
962                                                                                                         for(k = 0; k < smd->domain->amplify; k++)
963                                                                                                         {
964                                                                                                                 index = smoke_get_index(smd->domain->amplify * cell[0] + i, bigres[0], smd->domain->amplify * cell[1] + j, bigres[1], smd->domain->amplify * cell[2] + k);
965                                                                                                                 bigdensity[index] = 0.f;
966                                                                                                         }
967                                                                                 }
968                                                                         }
969                                                                 }
970                                                         }       
971                                                 }       
972                                         }
973
974                                         if(cnt_domain && !sds->fluid_group)
975                                                 go = go->next;
976                                         else
977                                                 base= base->next;
978                                 }
979                         }
980
981                         // do effectors
982                         /*
983                         if(sds->eff_group)
984                         {
985                                 for(go = sds->eff_group->gobject.first; go; go = go->next) 
986                                 {
987                                         if(go->ob)
988                                         {
989                                                 if(ob->pd)
990                                                 {
991                                                         
992                                                 }
993                                         }
994                                 }
995                         }
996                         */
997
998                         // do collisions        
999                         if(sds->coll_group || !cnt_domain)
1000                         {
1001                                 Object *otherobj = NULL;
1002                                 ModifierData *md = NULL;
1003
1004                                 if(cnt_domain && !sds->coll_group) // we use groups since we have 2 domains
1005                                         go = sds->coll_group->gobject.first;
1006                                 else
1007                                         base = scene->base.first;
1008
1009                                 while(base || go)
1010                                 {
1011                                         otherobj = NULL;
1012
1013                                         if(cnt_domain && !sds->coll_group) 
1014                                         {
1015                                                 if(go->ob)
1016                                                         otherobj = go->ob;
1017                                         }
1018                                         else
1019                                                 otherobj = base->object;
1020
1021                                         if(!otherobj)
1022                                         {
1023                                                 if(cnt_domain && !sds->coll_group)
1024                                                         go = go->next;
1025                                                 else
1026                                                         base= base->next;
1027
1028                                                 continue;
1029                                         }
1030                         
1031                                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
1032                                         
1033                                         // check for active smoke modifier
1034                                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
1035                                         {
1036                                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
1037
1038                                                 if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll)
1039                                                 {
1040                                                         // we got nice collision object
1041                                                         SmokeCollSettings *scs = smd2->coll;
1042                                                         int cell[3];
1043                                                         size_t index = 0;
1044                                                         size_t i, j;
1045                                                         unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid);
1046
1047                                                         // int BLI_bvhtree_find_nearest(BVHTree *tree, const float *co, BVHTreeNearest *nearest, BVHTree_NearestPointCallback callback, void *userdata);
1048
1049                                                         for(i = 0; i < scs->numpoints; i++)
1050                                                         {
1051                                                                 int badcell = 0;
1052
1053                                                                 // 1. get corresponding cell
1054                                                                 get_cell(smd, &scs->points[3 * i], cell, 0);
1055                                                         
1056                                                                 // check if cell is valid (in the domain boundary)
1057                                                                 for(j = 0; j < 3; j++)
1058                                                                         if((cell[j] > sds->res[j] - 1) || (cell[j] < 0))
1059                                                                         {
1060                                                                                 badcell = 1;
1061                                                                                 break;
1062                                                                         }
1063                                                                                 
1064                                                                 if(badcell)
1065                                                                         continue;
1066
1067                                                                 // 2. set cell values (heat, density and velocity)
1068                                                                 index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
1069                                                                 
1070                                                                 // printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]);
1071                                                                 // printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index);
1072                                                                         
1073                                                                 obstacles[index] = 1;
1074
1075                                                                 // for moving gobstacles
1076                                                                 /*
1077                                                                 const LbmFloat maxVelVal = 0.1666;
1078                                                                 const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
1079
1080                                                                 LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { 
1081                                                                 const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; 
1082                                                                 USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); 
1083                                                                 if(usqr>maxusqr) { 
1084                                                                         // cutoff at maxVelVal 
1085                                                                         for(int jj=0; jj<3; jj++) { 
1086                                                                                 if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  
1087                                                                                 if(objvel[jj]<0.) objvel[jj] = -maxVelVal; 
1088                                                                         } 
1089                                                                 } } 
1090
1091                                                                 const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); 
1092                                                                 const LbmVec oldov=objvel; // debug
1093                                                                 objvel = vec2L((*pNormals)[n]) *dp;
1094                                                                 */
1095                                                         }
1096                                                 }
1097                                         }
1098
1099                                         if(cnt_domain && !sds->coll_group)
1100                                                 go = go->next;
1101                                         else
1102                                                 base= base->next;
1103                                 }
1104                         }
1105                         
1106                         // set new time
1107                         smd->time = scene->r.cfra;
1108
1109                         // simulate the actual smoke (c++ code in intern/smoke)
1110                         smoke_step(sds->fluid);
1111
1112                         tend();
1113                         printf ( "Frame: %d, Time: %f\n", (int)smd->time, ( float ) tval() );
1114                 }
1115                 else if(scene->r.cfra < smd->time)
1116                 {
1117                         // we got back in time, reset smoke in this case (TODO: use cache later)
1118                         smd->time = scene->r.cfra;
1119                         smokeModifier_reset(smd);
1120                 }
1121         }
1122 }
1123
1124 // update necessary information for 3dview
1125 void smoke_prepare_View(SmokeModifierData *smd, float *light)
1126 {
1127         float *density = NULL;
1128         int x, y, z;
1129
1130         if(!smd->domain->tray)
1131         {
1132                 // TRay is for self shadowing
1133                 smd->domain->tray = MEM_callocN(sizeof(float)*smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2], "Smoke_tRay");
1134         }
1135         if(!smd->domain->tvox)
1136         {
1137                 // TVox is for tranaparency
1138                 smd->domain->tvox = MEM_callocN(sizeof(float)*smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2], "Smoke_tVox");
1139         }
1140
1141         // update 3dview
1142         density = smoke_get_density(smd->domain->fluid);
1143         for(x = 0; x < smd->domain->res[0]; x++)
1144                         for(y = 0; y < smd->domain->res[1]; y++)
1145                                 for(z = 0; z < smd->domain->res[2]; z++)
1146                                 {
1147                                         size_t index;
1148
1149                                         index = smoke_get_index(x, smd->domain->res[0], y, smd->domain->res[1], z);
1150                                         // Transparency computation
1151                                         // formula taken from "Visual Simulation of Smoke" / Fedkiw et al. pg. 4
1152                                         // T_vox = exp(-C_ext * h)
1153                                         // C_ext/sigma_t = density * C_ext
1154                                         smoke_set_tvox(smd, index, exp(-density[index] * 7.0 * smd->domain->dx));
1155         }
1156         smoke_calc_transparency(smd, light, 0);
1157 }
1158
1159 // update necessary information for 3dview ("high res" option)
1160 void smoke_prepare_bigView(SmokeModifierData *smd, float *light)
1161 {
1162         float *density = NULL;
1163         size_t i = 0;
1164         int bigres[3];
1165
1166         smoke_get_bigres(smd->domain->fluid, bigres);
1167
1168         if(!smd->domain->traybig)
1169         {
1170                 // TRay is for self shadowing
1171                 smd->domain->traybig = MEM_callocN(sizeof(float)*bigres[0]*bigres[1]*bigres[2], "Smoke_tRayBig");
1172         }
1173         if(!smd->domain->tvoxbig)
1174         {
1175                 // TVox is for tranaparency
1176                 smd->domain->tvoxbig = MEM_callocN(sizeof(float)*bigres[0]*bigres[1]*bigres[2], "Smoke_tVoxBig");
1177         }
1178
1179         density = smoke_get_bigdensity(smd->domain->fluid);
1180         for (i = 0; i < bigres[0] * bigres[1] * bigres[2]; i++)
1181         {
1182                 // Transparency computation
1183                 // formula taken from "Visual Simulation of Smoke" / Fedkiw et al. pg. 4
1184                 // T_vox = exp(-C_ext * h)
1185                 // C_ext/sigma_t = density * C_ext
1186                 smoke_set_bigtvox(smd, i, exp(-density[i] * 7.0 * smd->domain->dx / smd->domain->amplify) );
1187         }
1188         smoke_calc_transparency(smd, light, 1);
1189 }
1190
1191
1192 float smoke_get_tvox(SmokeModifierData *smd, size_t index)
1193 {
1194         return smd->domain->tvox[index];
1195 }
1196
1197 void smoke_set_tvox(SmokeModifierData *smd, size_t index, float tvox)
1198 {
1199         smd->domain->tvox[index] = tvox;
1200 }
1201
1202 float smoke_get_tray(SmokeModifierData *smd, size_t index)
1203 {
1204         return smd->domain->tray[index];
1205 }
1206
1207 void smoke_set_tray(SmokeModifierData *smd, size_t index, float transparency)
1208 {
1209         smd->domain->tray[index] = transparency;
1210 }
1211
1212 float smoke_get_bigtvox(SmokeModifierData *smd, size_t index)
1213 {
1214         return smd->domain->tvoxbig[index];
1215 }
1216
1217 void smoke_set_bigtvox(SmokeModifierData *smd, size_t index, float tvox)
1218 {
1219         smd->domain->tvoxbig[index] = tvox;
1220 }
1221
1222 float smoke_get_bigtray(SmokeModifierData *smd, size_t index)
1223 {
1224         return smd->domain->traybig[index];
1225 }
1226
1227 void smoke_set_bigtray(SmokeModifierData *smd, size_t index, float transparency)
1228 {
1229         smd->domain->traybig[index] = transparency;
1230 }
1231
1232 long long smoke_get_mem_req(int xres, int yres, int zres, int amplify)
1233 {
1234           int totalCells = xres * yres * zres;
1235           int amplifiedCells = totalCells * amplify * amplify * amplify;
1236
1237           // print out memory requirements
1238           long long int coarseSize = sizeof(float) * totalCells * 22 +
1239                            sizeof(unsigned char) * totalCells;
1240
1241           long long int fineSize = sizeof(float) * amplifiedCells * 7 + // big grids
1242                          sizeof(float) * totalCells * 8 +     // small grids
1243                          sizeof(float) * 128 * 128 * 128;     // noise tile
1244
1245           long long int totalMB = (coarseSize + fineSize) / (1024 * 1024);
1246
1247           return totalMB;
1248 }
1249
1250
1251 static void calc_voxel_transp(SmokeModifierData *smd, int *pixel, float *tRay)
1252 {
1253         // printf("Pixel(%d, %d, %d)\n", pixel[0], pixel[1], pixel[2]);
1254         const size_t index = smoke_get_index(pixel[0], smd->domain->res[0], pixel[1], smd->domain->res[1], pixel[2]);
1255
1256         // T_ray *= T_vox
1257         *tRay *= smoke_get_tvox(smd, index);
1258 }
1259
1260 static void calc_voxel_transp_big(SmokeModifierData *smd, int *pixel, float *tRay)
1261 {
1262         int bigres[3];
1263         size_t index;
1264
1265         smoke_get_bigres(smd->domain->fluid, bigres);
1266         index = smoke_get_index(pixel[0], bigres[0], pixel[1], bigres[1], pixel[2]);
1267
1268         /*
1269         if(index > bigres[0]*bigres[1]*bigres[2])
1270                 printf("pixel[0]: %d, [1]: %d, [2]: %d\n", pixel[0], pixel[1], pixel[2]);
1271         */
1272
1273         // T_ray *= T_vox
1274         *tRay *= smoke_get_bigtvox(smd, index);
1275 }
1276
1277 static void bresenham_linie_3D(SmokeModifierData *smd, int x1, int y1, int z1, int x2, int y2, int z2, float *tRay, int big)
1278 {
1279     int dx, dy, dz, i, l, m, n, x_inc, y_inc, z_inc, err_1, err_2, dx2, dy2, dz2;
1280     int pixel[3];
1281
1282     pixel[0] = x1;
1283     pixel[1] = y1;
1284     pixel[2] = z1;
1285
1286     dx = x2 - x1;
1287     dy = y2 - y1;
1288     dz = z2 - z1;
1289
1290     x_inc = (dx < 0) ? -1 : 1;
1291     l = abs(dx);
1292     y_inc = (dy < 0) ? -1 : 1;
1293     m = abs(dy);
1294     z_inc = (dz < 0) ? -1 : 1;
1295     n = abs(dz);
1296     dx2 = l << 1;
1297     dy2 = m << 1;
1298     dz2 = n << 1;
1299
1300     if ((l >= m) && (l >= n)) {
1301         err_1 = dy2 - l;
1302         err_2 = dz2 - l;
1303         for (i = 0; i < l; i++) {
1304                 if(!big)
1305                                 calc_voxel_transp(smd, pixel, tRay);
1306                         else
1307                                 calc_voxel_transp_big(smd, pixel, tRay);
1308                 if(*tRay < 0.0f)
1309                         return;
1310             if (err_1 > 0) {
1311                 pixel[1] += y_inc;
1312                 err_1 -= dx2;
1313             }
1314             if (err_2 > 0) {
1315                 pixel[2] += z_inc;
1316                 err_2 -= dx2;
1317             }
1318             err_1 += dy2;
1319             err_2 += dz2;
1320             pixel[0] += x_inc;
1321         }
1322     } else if ((m >= l) && (m >= n)) {
1323         err_1 = dx2 - m;
1324         err_2 = dz2 - m;
1325         for (i = 0; i < m; i++) {
1326                 if(!big)
1327                                 calc_voxel_transp(smd, pixel, tRay);
1328                         else
1329                                 calc_voxel_transp_big(smd, pixel, tRay);
1330                 if(*tRay < 0.0f)
1331                         return;
1332             if (err_1 > 0) {
1333                 pixel[0] += x_inc;
1334                 err_1 -= dy2;
1335             }
1336             if (err_2 > 0) {
1337                 pixel[2] += z_inc;
1338                 err_2 -= dy2;
1339             }
1340             err_1 += dx2;
1341             err_2 += dz2;
1342             pixel[1] += y_inc;
1343         }
1344     } else {
1345         err_1 = dy2 - n;
1346         err_2 = dx2 - n;
1347         for (i = 0; i < n; i++) {
1348                 if(!big)
1349                                 calc_voxel_transp(smd, pixel, tRay);
1350                         else
1351                                 calc_voxel_transp_big(smd, pixel, tRay);
1352                 if(*tRay < 0.0f)
1353                         return;
1354             if (err_1 > 0) {
1355                 pixel[1] += y_inc;
1356                 err_1 -= dz2;
1357             }
1358             if (err_2 > 0) {
1359                 pixel[0] += x_inc;
1360                 err_2 -= dz2;
1361             }
1362             err_1 += dy2;
1363             err_2 += dx2;
1364             pixel[2] += z_inc;
1365         }
1366     }
1367     if(!big)
1368         calc_voxel_transp(smd, pixel, tRay);
1369     else
1370         calc_voxel_transp_big(smd, pixel, tRay);
1371 }
1372
1373 static void get_cell(struct SmokeModifierData *smd, float *pos, int *cell, int correct)
1374 {
1375         float tmp[3];
1376
1377         VECSUB(tmp, pos, smd->domain->p0);
1378         VecMulf(tmp, 1.0 / smd->domain->dx);
1379
1380         if(correct)
1381         {
1382                 cell[0] = MIN2(smd->domain->res[0] - 1, MAX2(0, (int)floor(tmp[0])));
1383                 cell[1] = MIN2(smd->domain->res[1] - 1, MAX2(0, (int)floor(tmp[1])));
1384                 cell[2] = MIN2(smd->domain->res[2] - 1, MAX2(0, (int)floor(tmp[2])));
1385         }
1386         else
1387         {
1388                 cell[0] = (int)floor(tmp[0]);
1389                 cell[1] = (int)floor(tmp[1]);
1390                 cell[2] = (int)floor(tmp[2]);
1391         }
1392 }
1393 static void get_bigcell(struct SmokeModifierData *smd, float *pos, int *cell, int correct)
1394 {
1395         float tmp[3];
1396         int res[3];
1397         smoke_get_bigres(smd->domain->fluid, res);
1398
1399         VECSUB(tmp, pos, smd->domain->p0);
1400
1401         VecMulf(tmp, smd->domain->amplify / smd->domain->dx );
1402
1403         if(correct)
1404         {
1405                 cell[0] = MIN2(res[0] - 1, MAX2(0, (int)floor(tmp[0])));
1406                 cell[1] = MIN2(res[1] - 1, MAX2(0, (int)floor(tmp[1])));
1407                 cell[2] = MIN2(res[2] - 1, MAX2(0, (int)floor(tmp[2])));
1408         }
1409         else
1410         {
1411                 cell[0] = (int)floor(tmp[0]);
1412                 cell[1] = (int)floor(tmp[1]);
1413                 cell[2] = (int)floor(tmp[2]);
1414         }
1415 }
1416
1417
1418 void smoke_calc_transparency(struct SmokeModifierData *smd, float *light, int big)
1419 {
1420         int x, y, z;
1421         float bv[6];
1422         int res[3];
1423         float bigfactor = 1.0;
1424
1425         // x
1426         bv[0] = smd->domain->p0[0];
1427         bv[1] = smd->domain->p1[0];
1428         // y
1429         bv[2] = smd->domain->p0[1];
1430         bv[3] = smd->domain->p1[1];
1431         // z
1432         bv[4] = smd->domain->p0[2];
1433         bv[5] = smd->domain->p1[2];
1434 /*
1435         printf("bv[0]: %f, [1]: %f, [2]: %f, [3]: %f, [4]: %f, [5]: %f\n", bv[0], bv[1], bv[2], bv[3], bv[4], bv[5]);
1436
1437         printf("p0[0]: %f, p0[1]: %f, p0[2]: %f\n", smd->domain->p0[0], smd->domain->p0[1], smd->domain->p0[2]);
1438         printf("p1[0]: %f, p1[1]: %f, p1[2]: %f\n", smd->domain->p1[0], smd->domain->p1[1], smd->domain->p1[2]);
1439         printf("dx: %f, amp: %d\n", smd->domain->dx, smd->domain->amplify);
1440 */
1441         if(!big)
1442         {
1443                 res[0] = smd->domain->res[0];
1444                 res[1] = smd->domain->res[1];
1445                 res[2] = smd->domain->res[2];
1446         }
1447         else
1448         {
1449                 smoke_get_bigres(smd->domain->fluid, res);
1450                 bigfactor = 1.0 / smd->domain->amplify;
1451         }
1452
1453 #pragma omp parallel for schedule(static) private(y, z) shared(big, smd, light, res, bigfactor)
1454         for(x = 0; x < res[0]; x++)
1455                 for(y = 0; y < res[1]; y++)
1456                         for(z = 0; z < res[2]; z++)
1457                         {
1458                                 float voxelCenter[3];
1459                                 size_t index;
1460                                 float pos[3];
1461                                 int cell[3];
1462                                 float tRay = 1.0;
1463
1464                                 index = smoke_get_index(x, res[0], y, res[1], z);
1465
1466                                 // voxelCenter = m_voxelarray[i].GetCenter();
1467                                 voxelCenter[0] = smd->domain->p0[0] + smd->domain->dx * bigfactor * x + smd->domain->dx * bigfactor * 0.5;
1468                                 voxelCenter[1] = smd->domain->p0[1] + smd->domain->dx * bigfactor * y + smd->domain->dx * bigfactor * 0.5;
1469                                 voxelCenter[2] = smd->domain->p0[2] + smd->domain->dx * bigfactor * z + smd->domain->dx * bigfactor * 0.5;
1470
1471                                 // printf("vc[0]: %f, vc[1]: %f, vc[2]: %f\n", voxelCenter[0], voxelCenter[1], voxelCenter[2]);
1472                                 // printf("light[0]: %f, light[1]: %f, light[2]: %f\n", light[0], light[1], light[2]);
1473
1474                                 // get starting position (in voxel coords)
1475                                 if(BLI_bvhtree_bb_raycast(bv, light, voxelCenter, pos) > FLT_EPSILON)
1476                                 {
1477                                         // we're ouside
1478                                         // printf("out: pos[0]: %f, pos[1]: %f, pos[2]: %f\n", pos[0], pos[1], pos[2]);
1479                                         if(!big)
1480                                                 get_cell(smd, pos, cell, 1);
1481                                         else
1482                                                 get_bigcell(smd, pos, cell, 1);
1483                                 }
1484                                 else
1485                                 {
1486                                         // printf("in: pos[0]: %f, pos[1]: %f, pos[2]: %f\n", light[0], light[1], light[2]);
1487                                         // we're inside
1488                                         if(!big)
1489                                                 get_cell(smd, light, cell, 1);
1490                                         else
1491                                                 get_bigcell(smd, light, cell, 1);
1492                                 }
1493
1494                                 // printf("cell - [0]: %d, [1]: %d, [2]: %d\n", cell[0], cell[1], cell[2]);
1495                                 bresenham_linie_3D(smd, cell[0], cell[1], cell[2], x, y, z, &tRay, big);
1496
1497                                 if(!big)
1498                                         smoke_set_tray(smd, index, tRay);
1499                                 else
1500                                         smoke_set_bigtray(smd, index, tRay);
1501                         }
1502 }
1503