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