style cleanup, brackets in else/if, some indentation.
[blender.git] / source / blender / render / intern / source / volume_precache.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Matt Ebb, Ra˙l Fern·ndez Hern·ndez (Farsthary).
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/render/intern/source/volume_precache.c
29  *  \ingroup render
30  */
31
32
33 #include <math.h>
34 #include <stdlib.h>
35 #include <string.h>
36 #include <float.h>
37
38 #include "MEM_guardedalloc.h"
39
40 #include "BLI_blenlib.h"
41 #include "BLI_math.h"
42 #include "BLI_threads.h"
43 #include "BLI_voxel.h"
44 #include "BLI_utildefines.h"
45
46 #include "PIL_time.h"
47
48 #include "RE_shader_ext.h"
49
50 #include "DNA_material_types.h"
51
52 #include "rayintersection.h"
53 #include "rayobject.h"
54 #include "render_types.h"
55 #include "rendercore.h"
56 #include "renderdatabase.h"
57 #include "volumetric.h"
58 #include "volume_precache.h"
59
60 #if defined( _MSC_VER ) && !defined( __cplusplus )
61 # define inline __inline
62 #endif // defined( _MSC_VER ) && !defined( __cplusplus )
63
64 #include "BKE_global.h"
65
66 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
67 /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
68 /* only to be used here in this file, it's for speed */
69 extern struct Render R;
70 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
71
72 /* *** utility code to set up an individual raytree for objectinstance, for checking inside/outside *** */
73
74 /* Recursive test for intersections, from a point inside the mesh, to outside
75  * Number of intersections (depth) determine if a point is inside or outside the mesh */
76 static int intersect_outside_volume(RayObject *tree, Isect *isect, float *offset, int limit, int depth)
77 {
78         if (limit == 0) return depth;
79         
80         if (RE_rayobject_raycast(tree, isect)) {
81                 
82                 isect->start[0] = isect->start[0] + isect->dist*isect->dir[0];
83                 isect->start[1] = isect->start[1] + isect->dist*isect->dir[1];
84                 isect->start[2] = isect->start[2] + isect->dist*isect->dir[2];
85                 
86                 isect->dist = FLT_MAX;
87                 isect->skip = RE_SKIP_VLR_NEIGHBOUR;
88                 isect->orig.face= isect->hit.face;
89                 isect->orig.ob= isect->hit.ob;
90                 
91                 return intersect_outside_volume(tree, isect, offset, limit-1, depth+1);
92         } else {
93                 return depth;
94         }
95 }
96
97 /* Uses ray tracing to check if a point is inside or outside an ObjectInstanceRen */
98 static int point_inside_obi(RayObject *tree, ObjectInstanceRen *UNUSED(obi), float *co)
99 {
100         Isect isect= {{0}};
101         float dir[3] = {0.0f,0.0f,1.0f};
102         int final_depth=0, depth=0, limit=20;
103         
104         /* set up the isect */
105         copy_v3_v3(isect.start, co);
106         copy_v3_v3(isect.dir, dir);
107         isect.mode= RE_RAY_MIRROR;
108         isect.last_hit= NULL;
109         isect.lay= -1;
110         
111         isect.dist = FLT_MAX;
112         isect.orig.face= NULL;
113         isect.orig.ob = NULL;
114
115         final_depth = intersect_outside_volume(tree, &isect, dir, limit, depth);
116         
117         /* even number of intersections: point is outside
118          * odd number: point is inside */
119         if (final_depth % 2 == 0) return 0;
120         else return 1;
121 }
122
123 /* find the bounding box of an objectinstance in global space */
124 void global_bounds_obi(Render *re, ObjectInstanceRen *obi, float *bbmin, float *bbmax)
125 {
126         ObjectRen *obr = obi->obr;
127         VolumePrecache *vp = obi->volume_precache;
128         VertRen *ver= NULL;
129         float co[3];
130         int a;
131         
132         if (vp->bbmin != NULL && vp->bbmax != NULL) {
133                 copy_v3_v3(bbmin, vp->bbmin);
134                 copy_v3_v3(bbmax, vp->bbmax);
135                 return;
136         }
137         
138         vp->bbmin = MEM_callocN(sizeof(float)*3, "volume precache min boundbox corner");
139         vp->bbmax = MEM_callocN(sizeof(float)*3, "volume precache max boundbox corner");
140         
141         INIT_MINMAX(bbmin, bbmax);
142         
143         for(a=0; a<obr->totvert; a++) {
144                 if((a & 255)==0) ver= obr->vertnodes[a>>8].vert;
145                 else ver++;
146                 
147                 copy_v3_v3(co, ver->co);
148                 
149                 /* transformed object instance in camera space */
150                 if(obi->flag & R_TRANSFORMED)
151                         mul_m4_v3(obi->mat, co);
152                 
153                 /* convert to global space */
154                 mul_m4_v3(re->viewinv, co);
155                 
156                 DO_MINMAX(co, vp->bbmin, vp->bbmax);
157         }
158         
159         copy_v3_v3(bbmin, vp->bbmin);
160         copy_v3_v3(bbmax, vp->bbmax);
161         
162 }
163
164 /* *** light cache filtering *** */
165
166 static float get_avg_surrounds(float *cache, int *res, int xx, int yy, int zz)
167 {
168         int x, y, z, x_, y_, z_;
169         int added=0;
170         float tot=0.0f;
171         
172         for (z=-1; z <= 1; z++) {
173                 z_ = zz+z;
174                 if (z_ >= 0 && z_ <= res[2]-1) {
175                 
176                         for (y=-1; y <= 1; y++) {
177                                 y_ = yy+y;
178                                 if (y_ >= 0 && y_ <= res[1]-1) {
179                                 
180                                         for (x=-1; x <= 1; x++) {
181                                                 x_ = xx+x;
182                                                 if (x_ >= 0 && x_ <= res[0]-1) {
183                                                         const int i= V_I(x_, y_, z_, res);
184                                                         
185                                                         if (cache[i] > 0.0f) {
186                                                                 tot += cache[i];
187                                                                 added++;
188                                                         }
189                                                         
190                                                 }
191                                         }
192                                 }
193                         }
194                 }
195         }
196         
197         if (added > 0) tot /= added;
198         
199         return tot;
200 }
201
202 /* function to filter the edges of the light cache, where there was no volume originally.
203  * For each voxel which was originally external to the mesh, it finds the average values of
204  * the surrounding internal voxels and sets the original external voxel to that average amount.
205  * Works almost a bit like a 'dilate' filter */
206 static void lightcache_filter(VolumePrecache *vp)
207 {
208         int x, y, z;
209
210         for (z=0; z < vp->res[2]; z++) {
211                 for (y=0; y < vp->res[1]; y++) {
212                         for (x=0; x < vp->res[0]; x++) {
213                                 /* trigger for outside mesh */
214                                 const int i= V_I(x, y, z, vp->res);
215                                 
216                                 if (vp->data_r[i] < -0.f)
217                                         vp->data_r[i] = get_avg_surrounds(vp->data_r, vp->res, x, y, z);
218                                 if (vp->data_g[i] < -0.f)
219                                         vp->data_g[i] = get_avg_surrounds(vp->data_g, vp->res, x, y, z);
220                                 if (vp->data_b[i] < -0.f)
221                                         vp->data_b[i] = get_avg_surrounds(vp->data_b, vp->res, x, y, z);
222                         }
223                 }
224         }
225 }
226
227 #if 0
228 static void lightcache_filter2(VolumePrecache *vp)
229 {
230         int x, y, z;
231         float *new_r, *new_g, *new_b;
232         int field_size = vp->res[0]*vp->res[1]*vp->res[2]*sizeof(float);
233         
234         new_r = MEM_mallocN(field_size, "temp buffer for light cache filter r channel");
235         new_g = MEM_mallocN(field_size, "temp buffer for light cache filter g channel");
236         new_b = MEM_mallocN(field_size, "temp buffer for light cache filter b channel");
237         
238         memcpy(new_r, vp->data_r, field_size);
239         memcpy(new_g, vp->data_g, field_size);
240         memcpy(new_b, vp->data_b, field_size);
241         
242         for (z=0; z < vp->res[2]; z++) {
243                 for (y=0; y < vp->res[1]; y++) {
244                         for (x=0; x < vp->res[0]; x++) {
245                                 /* trigger for outside mesh */
246                                 const int i= V_I(x, y, z, vp->res);
247                                 if (vp->data_r[i] < -0.f)
248                                         new_r[i] = get_avg_surrounds(vp->data_r, vp->res, x, y, z);
249                                 if (vp->data_g[i] < -0.f)
250                                         new_g[i] = get_avg_surrounds(vp->data_g, vp->res, x, y, z);
251                                 if (vp->data_b[i] < -0.f)
252                                         new_b[i] = get_avg_surrounds(vp->data_b, vp->res, x, y, z);
253                         }
254                 }
255         }
256         
257         SWAP(float *, vp->data_r, new_r);
258         SWAP(float *, vp->data_g, new_g);
259         SWAP(float *, vp->data_b, new_b);
260         
261         if (new_r) { MEM_freeN(new_r); new_r=NULL; }
262         if (new_g) { MEM_freeN(new_g); new_g=NULL; }
263         if (new_b) { MEM_freeN(new_b); new_b=NULL; }
264 }
265 #endif
266
267 static inline int ms_I(int x, int y, int z, int *n) //has a pad of 1 voxel surrounding the core for boundary simulation
268
269         /* different ordering to light cache */
270         return x*(n[1]+2)*(n[2]+2) + y*(n[2]+2) + z;    
271 }
272
273 static inline int v_I_pad(int x, int y, int z, int *n) //has a pad of 1 voxel surrounding the core for boundary simulation
274
275         /* same ordering to light cache, with padding */
276         return z*(n[1]+2)*(n[0]+2) + y*(n[0]+2) + x;    
277 }
278
279 static inline int lc_to_ms_I(int x, int y, int z, int *n)
280
281         /* converting light cache index to multiple scattering index */
282         return (x-1)*(n[1]*n[2]) + (y-1)*(n[2]) + z-1;
283 }
284
285 /* *** multiple scattering approximation *** */
286
287 /* get the total amount of light energy in the light cache. used to normalize after multiple scattering */
288 static float total_ss_energy(Render *re, int do_test_break, VolumePrecache *vp)
289 {
290         int x, y, z;
291         int *res = vp->res;
292         float energy=0.f;
293         
294         for (z=0; z < res[2]; z++) {
295                 for (y=0; y < res[1]; y++) {
296                         for (x=0; x < res[0]; x++) {
297                                 const int i=V_I(x, y, z, res);
298                         
299                                 if (vp->data_r[i] > 0.f) energy += vp->data_r[i];
300                                 if (vp->data_g[i] > 0.f) energy += vp->data_g[i];
301                                 if (vp->data_b[i] > 0.f) energy += vp->data_b[i];
302                         }
303                 }
304
305                 if (do_test_break && re->test_break(re->tbh)) break;
306         }
307         
308         return energy;
309 }
310
311 static float total_ms_energy(Render *re, int do_test_break, float *sr, float *sg, float *sb, int *res)
312 {
313         int x, y, z;
314         float energy=0.f;
315         
316         for (z=1;z<=res[2];z++) {
317                 for (y=1;y<=res[1];y++) {
318                         for (x=1;x<=res[0];x++) {
319                                 const int i = ms_I(x,y,z,res);
320                                 
321                                 if (sr[i] > 0.f) energy += sr[i];
322                                 if (sg[i] > 0.f) energy += sg[i];
323                                 if (sb[i] > 0.f) energy += sb[i];
324                         }
325                 }
326
327                 if (do_test_break && re->test_break(re->tbh)) break;
328         }
329         
330         return energy;
331 }
332
333 static void ms_diffuse(Render *re, int do_test_break, float *x0, float *x, float diff, int *n) //n is the unpadded resolution
334 {
335         int i, j, k, l;
336         const float dt = VOL_MS_TIMESTEP;
337         size_t size = n[0]*n[1]*n[2];
338         const float a = dt*diff*size;
339         
340         for (l=0; l<20; l++)
341         {
342                 for (k=1; k<=n[2]; k++)
343                 {
344                         for (j=1; j<=n[1]; j++)
345                         {
346                                 for (i=1; i<=n[0]; i++)
347                                 {
348                                    x[v_I_pad(i,j,k,n)] = (x0[v_I_pad(i,j,k,n)]) + a*(   x0[v_I_pad(i-1,j,k,n)]+ x0[v_I_pad(i+1,j,k,n)]+ x0[v_I_pad(i,j-1,k,n)]+
349                                                                                                                                                 x0[v_I_pad(i,j+1,k,n)]+ x0[v_I_pad(i,j,k-1,n)]+x0[v_I_pad(i,j,k+1,n)]
350                                                                                                                                                 ) / (1+6*a);
351                                 }
352                         }
353
354                         if (do_test_break && re->test_break(re->tbh)) break;
355                 }
356
357                 if (re->test_break(re->tbh)) break;
358         }
359 }
360
361 static void multiple_scattering_diffusion(Render *re, VolumePrecache *vp, Material *ma)
362 {
363         const float diff = ma->vol.ms_diff * 0.001f;    /* compensate for scaling for a nicer UI range */
364         const int simframes = (int)(ma->vol.ms_spread * (float)MAX3(vp->res[0], vp->res[1], vp->res[2]));
365         const int shade_type = ma->vol.shade_type;
366         float fac = ma->vol.ms_intensity;
367         
368         int x, y, z, m;
369         int *n = vp->res;
370         const int size = (n[0]+2)*(n[1]+2)*(n[2]+2);
371         const int do_test_break = (size > 100000);
372         double time, lasttime= PIL_check_seconds_timer();
373         float total;
374         float c=1.0f;
375         float origf;    /* factor for blending in original light cache */
376         float energy_ss, energy_ms;
377
378         float *sr0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
379         float *sr=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
380         float *sg0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
381         float *sg=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
382         float *sb0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
383         float *sb=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
384
385         total = (float)(n[0]*n[1]*n[2]*simframes);
386         
387         energy_ss = total_ss_energy(re, do_test_break, vp);
388         
389         /* Scattering as diffusion pass */
390         for (m=0; m<simframes; m++)
391         {
392                 /* add sources */
393                 for (z=1; z<=n[2]; z++)
394                 {
395                         for (y=1; y<=n[1]; y++)
396                         {
397                                 for (x=1; x<=n[0]; x++)
398                                 {
399                                         const int i = lc_to_ms_I(x, y ,z, n);   //lc index                                      
400                                         const int j = ms_I(x, y, z, n);                 //ms index
401                                         
402                                         time= PIL_check_seconds_timer();
403                                         c++;                                                                            
404                                         if (vp->data_r[i] > 0.0f)
405                                                 sr[j] += vp->data_r[i];
406                                         if (vp->data_g[i] > 0.0f)
407                                                 sg[j] += vp->data_g[i];
408                                         if (vp->data_b[i] > 0.0f)
409                                                 sb[j] += vp->data_b[i];
410                                         
411                                         /* Displays progress every second */
412                                         if(time-lasttime>1.0) {
413                                                 char str[64];
414                                                 BLI_snprintf(str, sizeof(str), "Simulating multiple scattering: %d%%", (int)(100.0f * (c / total)));
415                                                 re->i.infostr= str;
416                                                 re->stats_draw(re->sdh, &re->i);
417                                                 re->i.infostr= NULL;
418                                                 lasttime= time;
419                                         }
420                                 }
421                         }
422
423                         if (do_test_break && re->test_break(re->tbh)) break;
424                 }
425
426                 if (re->test_break(re->tbh)) break;
427
428                 SWAP(float *,sr,sr0);
429                 SWAP(float *,sg,sg0);
430                 SWAP(float *,sb,sb0);
431
432                 /* main diffusion simulation */
433                 ms_diffuse(re, do_test_break, sr0, sr, diff, n);
434                 ms_diffuse(re, do_test_break, sg0, sg, diff, n);
435                 ms_diffuse(re, do_test_break, sb0, sb, diff, n);
436                 
437                 if (re->test_break(re->tbh)) break;
438         }
439         
440         /* normalisation factor to conserve energy */
441         energy_ms = total_ms_energy(re, do_test_break, sr, sg, sb, n);
442         fac *= (energy_ss / energy_ms);
443         
444         /* blend multiple scattering back in the light cache */
445         if (shade_type == MA_VOL_SHADE_SHADEDPLUSMULTIPLE) {
446                 /* conserve energy - half single, half multiple */
447                 origf = 0.5f;
448                 fac *= 0.5f;
449         } else {
450                 origf = 0.0f;
451         }
452
453         for (z=1;z<=n[2];z++)
454         {
455                 for (y=1;y<=n[1];y++)
456                 {
457                         for (x=1;x<=n[0];x++)
458                         {
459                                 const int i = lc_to_ms_I(x, y ,z, n);   //lc index                                      
460                                 const int j = ms_I(x, y, z, n);                 //ms index
461                                 
462                                 vp->data_r[i] = origf * vp->data_r[i] + fac * sr[j];
463                                 vp->data_g[i] = origf * vp->data_g[i] + fac * sg[j];
464                                 vp->data_b[i] = origf * vp->data_b[i] + fac * sb[j];
465                         }
466                 }
467
468                 if (do_test_break && re->test_break(re->tbh)) break;
469         }
470
471         MEM_freeN(sr0);
472         MEM_freeN(sr);
473         MEM_freeN(sg0);
474         MEM_freeN(sg);
475         MEM_freeN(sb0);
476         MEM_freeN(sb);
477 }
478
479
480
481 #if 0 // debug stuff
482 static void *vol_precache_part_test(void *data)
483 {
484         VolPrecachePart *pa = data;
485
486         printf("part number: %d \n", pa->num);
487         printf("done: %d \n", pa->done);
488         printf("x min: %d   x max: %d \n", pa->minx, pa->maxx);
489         printf("y min: %d   y max: %d \n", pa->miny, pa->maxy);
490         printf("z min: %d   z max: %d \n", pa->minz, pa->maxz);
491
492         return NULL;
493 }
494 #endif
495
496 typedef struct VolPrecacheQueue {
497         ThreadQueue *work;
498         ThreadQueue *done;
499 } VolPrecacheQueue;
500
501 /* Iterate over the 3d voxel grid, and fill the voxels with scattering information
502  *
503  * It's stored in memory as 3 big float grids next to each other, one for each RGB channel.
504  * I'm guessing the memory alignment may work out better this way for the purposes
505  * of doing linear interpolation, but I haven't actually tested this theory! :)
506  */
507 static void *vol_precache_part(void *data)
508 {
509         VolPrecacheQueue *queue = (VolPrecacheQueue*)data;
510         VolPrecachePart *pa;
511
512         while ((pa = BLI_thread_queue_pop(queue->work))) {
513                 ObjectInstanceRen *obi = pa->obi;
514                 RayObject *tree = pa->tree;
515                 ShadeInput *shi = pa->shi;
516                 float scatter_col[3] = {0.f, 0.f, 0.f};
517                 float co[3], cco[3], view[3];
518                 int x, y, z, i;
519                 int res[3];
520
521                 if (pa->re->test_break && pa->re->test_break(pa->re->tbh))
522                         break;
523
524                 res[0]= pa->res[0];
525                 res[1]= pa->res[1];
526                 res[2]= pa->res[2];
527
528                 for (z= pa->minz; z < pa->maxz; z++) {
529                         co[2] = pa->bbmin[2] + (pa->voxel[2] * (z + 0.5f));
530                         
531                         for (y= pa->miny; y < pa->maxy; y++) {
532                                 co[1] = pa->bbmin[1] + (pa->voxel[1] * (y + 0.5f));
533                                 
534                                 for (x=pa->minx; x < pa->maxx; x++) {
535                                         co[0] = pa->bbmin[0] + (pa->voxel[0] * (x + 0.5f));
536                                         
537                                         if (pa->re->test_break && pa->re->test_break(pa->re->tbh))
538                                                 break;
539                                         
540                                         /* convert from world->camera space for shading */
541                                         mul_v3_m4v3(cco, pa->viewmat, co);
542                                         
543                                         i= V_I(x, y, z, res);
544                                         
545                                         // don't bother if the point is not inside the volume mesh
546                                         if (!point_inside_obi(tree, obi, cco)) {
547                                                 obi->volume_precache->data_r[i] = -1.0f;
548                                                 obi->volume_precache->data_g[i] = -1.0f;
549                                                 obi->volume_precache->data_b[i] = -1.0f;
550                                                 continue;
551                                         }
552                                         
553                                         copy_v3_v3(view, cco);
554                                         normalize_v3(view);
555                                         vol_get_scattering(shi, scatter_col, cco, view);
556                                 
557                                         obi->volume_precache->data_r[i] = scatter_col[0];
558                                         obi->volume_precache->data_g[i] = scatter_col[1];
559                                         obi->volume_precache->data_b[i] = scatter_col[2];
560                                         
561                                 }
562                         }
563                 }
564
565                 BLI_thread_queue_push(queue->done, pa);
566         }
567         
568         return NULL;
569 }
570
571
572 static void precache_setup_shadeinput(Render *re, ObjectInstanceRen *obi, Material *ma, ShadeInput *shi)
573 {
574         memset(shi, 0, sizeof(ShadeInput)); 
575         shi->depth= 1;
576         shi->mask= 1;
577         shi->mat = ma;
578         shi->vlr = NULL;
579         memcpy(&shi->r, &shi->mat->r, 23*sizeof(float));        // note, keep this synced with render_types.h
580         shi->har= shi->mat->har;
581         shi->obi= obi;
582         shi->obr= obi->obr;
583         shi->lay = re->lay;
584 }
585
586 static void precache_init_parts(Render *re, RayObject *tree, ShadeInput *shi, ObjectInstanceRen *obi, int totthread, int *parts)
587 {
588         VolumePrecache *vp = obi->volume_precache;
589         int i=0, x, y, z;
590         float voxel[3];
591         int sizex, sizey, sizez;
592         float bbmin[3], bbmax[3];
593         int *res;
594         int minx, maxx;
595         int miny, maxy;
596         int minz, maxz;
597         
598         if (!vp) return;
599
600         BLI_freelistN(&re->volume_precache_parts);
601         
602         /* currently we just subdivide the box, number of threads per side */
603         parts[0] = parts[1] = parts[2] = totthread;
604         res = vp->res;
605         
606         /* using boundbox in worldspace */
607         global_bounds_obi(re, obi, bbmin, bbmax);
608         sub_v3_v3v3(voxel, bbmax, bbmin);
609         
610         voxel[0] /= (float)res[0];
611         voxel[1] /= (float)res[1];
612         voxel[2] /= (float)res[2];
613
614         for (x=0; x < parts[0]; x++) {
615                 sizex = ceil(res[0] / (float)parts[0]);
616                 minx = x * sizex;
617                 maxx = minx + sizex;
618                 maxx = (maxx>res[0])?res[0]:maxx;
619                 
620                 for (y=0; y < parts[1]; y++) {
621                         sizey = ceil(res[1] / (float)parts[1]);
622                         miny = y * sizey;
623                         maxy = miny + sizey;
624                         maxy = (maxy>res[1])?res[1]:maxy;
625                         
626                         for (z=0; z < parts[2]; z++) {
627                                 VolPrecachePart *pa= MEM_callocN(sizeof(VolPrecachePart), "new precache part");
628                                 
629                                 sizez = ceil(res[2] / (float)parts[2]);
630                                 minz = z * sizez;
631                                 maxz = minz + sizez;
632                                 maxz = (maxz>res[2])?res[2]:maxz;
633                                 
634                                 pa->re = re;
635                                 pa->num = i;
636                                 pa->tree = tree;
637                                 pa->shi = shi;
638                                 pa->obi = obi;
639                                 copy_m4_m4(pa->viewmat, re->viewmat);
640                                 
641                                 copy_v3_v3(pa->bbmin, bbmin);
642                                 copy_v3_v3(pa->voxel, voxel);
643                                 copy_v3_v3_int(pa->res, res);
644                                 
645                                 pa->minx = minx; pa->maxx = maxx;
646                                 pa->miny = miny; pa->maxy = maxy;
647                                 pa->minz = minz; pa->maxz = maxz;
648                                 
649                                 
650                                 BLI_addtail(&re->volume_precache_parts, pa);
651                                 
652                                 i++;
653                         }
654                 }
655         }
656 }
657
658 /* calculate resolution from bounding box in world space */
659 static int precache_resolution(Render *re, VolumePrecache *vp, ObjectInstanceRen *obi, int res)
660 {
661         float dim[3], div;
662         float bbmin[3], bbmax[3];
663         
664         /* bound box in global space */
665         global_bounds_obi(re, obi, bbmin, bbmax);
666         sub_v3_v3v3(dim, bbmax, bbmin);
667         
668         div = MAX3(dim[0], dim[1], dim[2]);
669         dim[0] /= div;
670         dim[1] /= div;
671         dim[2] /= div;
672         
673         vp->res[0] = ceil(dim[0] * res);
674         vp->res[1] = ceil(dim[1] * res);
675         vp->res[2] = ceil(dim[2] * res);
676         
677         if ((vp->res[0] < 1) || (vp->res[1] < 1) || (vp->res[2] < 1))
678                 return 0;
679         
680         return 1;
681 }
682
683 /* Precache a volume into a 3D voxel grid.
684  * The voxel grid is stored in the ObjectInstanceRen, 
685  * in camera space, aligned with the ObjectRen's bounding box.
686  * Resolution is defined by the user.
687  */
688 static void vol_precache_objectinstance_threads(Render *re, ObjectInstanceRen *obi, Material *ma)
689 {
690         VolumePrecache *vp;
691         VolPrecachePart *pa;
692         RayObject *tree;
693         ShadeInput shi;
694         ListBase threads;
695         VolPrecacheQueue queue;
696         int parts[3] = {1, 1, 1}, totparts;
697         
698         int counter=0;
699         int totthread = re->r.threads, thread;
700         
701         double time, lasttime= PIL_check_seconds_timer();
702         
703         R = *re;
704
705         /* create a raytree with just the faces of the instanced ObjectRen, 
706          * used for checking if the cached point is inside or outside. */
707         tree = makeraytree_object(&R, obi);
708         if (!tree) return;
709
710         vp = MEM_callocN(sizeof(VolumePrecache), "volume light cache");
711         obi->volume_precache = vp;
712         
713         if (!precache_resolution(re, vp, obi, ma->vol.precache_resolution)) {
714                 MEM_freeN(vp);
715                 vp = NULL;
716                 return;
717         }
718
719         vp->data_r = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data red channel");
720         vp->data_g = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data green channel");
721         vp->data_b = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data blue channel");
722         if (vp->data_r==NULL || vp->data_g==NULL || vp->data_b==NULL) {
723                 MEM_freeN(vp);
724                 return;
725         }
726
727         /* Need a shadeinput to calculate scattering */
728         precache_setup_shadeinput(re, obi, ma, &shi);
729         
730         precache_init_parts(re, tree, &shi, obi, totthread, parts);
731         totparts = parts[0] * parts[1] * parts[2];
732
733         /* setup work and done queues */
734         queue.work = BLI_thread_queue_init();
735         queue.done = BLI_thread_queue_init();
736         BLI_thread_queue_nowait(queue.work);
737
738         for(pa= re->volume_precache_parts.first; pa; pa= pa->next)
739                 BLI_thread_queue_push(queue.work, pa);
740         
741         /* launch threads */
742         BLI_init_threads(&threads, vol_precache_part, totthread);
743
744         for(thread= 0; thread<totthread; thread++)
745                 BLI_insert_thread(&threads, &queue);
746         
747         /* loop waiting for work to be done */
748         while(counter < totparts) {
749                 if(re->test_break && re->test_break(re->tbh))
750                         break;
751
752                 if(BLI_thread_queue_pop_timeout(queue.done, 50))
753                         counter++;
754
755                 time= PIL_check_seconds_timer();
756                 if(time-lasttime>1.0) {
757                         char str[64];
758                         BLI_snprintf(str, sizeof(str), "Precaching volume: %d%%", (int)(100.0f * ((float)counter / (float)totparts)));
759                         re->i.infostr= str;
760                         re->stats_draw(re->sdh, &re->i);
761                         re->i.infostr= NULL;
762                         lasttime= time;
763                 }
764         }
765         
766         /* free */
767         BLI_end_threads(&threads);
768         BLI_thread_queue_free(queue.work);
769         BLI_thread_queue_free(queue.done);
770         BLI_freelistN(&re->volume_precache_parts);
771         
772         if(tree) {
773                 //TODO: makeraytree_object creates a tree and saves it on OBI, if we free this tree we should also clear other pointers to it
774                 //RE_rayobject_free(tree);
775                 //tree= NULL;
776         }
777         
778         if (ELEM(ma->vol.shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SHADEDPLUSMULTIPLE)) {
779                 /* this should be before the filtering */
780                 multiple_scattering_diffusion(re, obi->volume_precache, ma);
781         }
782                 
783         lightcache_filter(obi->volume_precache);
784 }
785
786 static int using_lightcache(Material *ma)
787 {
788         return (((ma->vol.shadeflag & MA_VOL_PRECACHESHADING) && (ma->vol.shade_type == MA_VOL_SHADE_SHADED))
789                 || (ELEM(ma->vol.shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SHADEDPLUSMULTIPLE)));
790 }
791
792 /* loop through all objects (and their associated materials)
793  * marked for pre-caching in convertblender.c, and pre-cache them */
794 void volume_precache(Render *re)
795 {
796         ObjectInstanceRen *obi;
797         VolumeOb *vo;
798
799         re->i.infostr= "Volume preprocessing";
800         re->stats_draw(re->sdh, &re->i);
801
802         for(vo= re->volumes.first; vo; vo= vo->next) {
803                 if (using_lightcache(vo->ma)) {
804                         for(obi= re->instancetable.first; obi; obi= obi->next) {
805                                 if (obi->obr == vo->obr) {
806                                         vol_precache_objectinstance_threads(re, obi, vo->ma);
807
808                                         if(re->test_break && re->test_break(re->tbh))
809                                                 break;
810                                 }
811                         }
812
813                         if(re->test_break && re->test_break(re->tbh))
814                                 break;
815                 }
816         }
817         
818         re->i.infostr= NULL;
819         re->stats_draw(re->sdh, &re->i);
820 }
821
822 void free_volume_precache(Render *re)
823 {
824         ObjectInstanceRen *obi;
825         
826         for(obi= re->instancetable.first; obi; obi= obi->next) {
827                 if (obi->volume_precache != NULL) {
828                         MEM_freeN(obi->volume_precache->data_r);
829                         MEM_freeN(obi->volume_precache->data_g);
830                         MEM_freeN(obi->volume_precache->data_b);
831                         MEM_freeN(obi->volume_precache->bbmin);
832                         MEM_freeN(obi->volume_precache->bbmax);
833                         MEM_freeN(obi->volume_precache);
834                         obi->volume_precache = NULL;
835                 }
836         }
837         
838         BLI_freelistN(&re->volumes);
839 }
840
841 int point_inside_volume_objectinstance(Render *re, ObjectInstanceRen *obi, float *co)
842 {
843         RayObject *tree;
844         int inside=0;
845         
846         tree = makeraytree_object(re, obi);
847         if (!tree) return 0;
848         
849         inside = point_inside_obi(tree, obi, co);
850         
851         //TODO: makeraytree_object creates a tree and saves it on OBI, if we free this tree we should also clear other pointers to it
852         //RE_rayobject_free(tree);
853         //tree= NULL;
854         
855         return inside;
856 }
857