Rework of volume shading
[blender.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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, 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.
25  *
26  * ***** END GPL LICENSE BLOCK *****
27  */
28
29 #include <math.h>
30 #include <stdlib.h>
31 #include <string.h>
32 #include <float.h>
33
34 #include "MEM_guardedalloc.h"
35
36 #include "BLI_blenlib.h"
37 #include "BLI_arithb.h"
38 #include "BLI_threads.h"
39 #include "BLI_voxel.h"
40
41 #include "PIL_time.h"
42
43 #include "RE_shader_ext.h"
44 #include "RE_raytrace.h"
45
46 #include "DNA_material_types.h"
47
48 #include "render_types.h"
49 #include "renderdatabase.h"
50 #include "volumetric.h"
51 #include "volume_precache.h"
52
53 #if defined( _MSC_VER ) && !defined( __cplusplus )
54 # define inline __inline
55 #endif // defined( _MSC_VER ) && !defined( __cplusplus )
56
57 #include "BKE_global.h"
58
59 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
60 /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
61 /* only to be used here in this file, it's for speed */
62 extern struct Render R;
63 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
64
65 /* *** utility code to set up an individual raytree for objectinstance, for checking inside/outside *** */
66
67 /* Recursive test for intersections, from a point inside the mesh, to outside
68  * Number of intersections (depth) determine if a point is inside or outside the mesh */
69 int intersect_outside_volume(RayTree *tree, Isect *isect, float *offset, int limit, int depth)
70 {
71         if (limit == 0) return depth;
72         
73         if (RE_ray_tree_intersect(tree, isect)) {
74                 float hitco[3];
75                 
76                 hitco[0] = isect->start[0] + isect->labda*isect->vec[0];
77                 hitco[1] = isect->start[1] + isect->labda*isect->vec[1];
78                 hitco[2] = isect->start[2] + isect->labda*isect->vec[2];
79                 VecAddf(isect->start, hitco, offset);
80
81                 return intersect_outside_volume(tree, isect, offset, limit-1, depth+1);
82         } else {
83                 return depth;
84         }
85 }
86
87 /* Uses ray tracing to check if a point is inside or outside an ObjectInstanceRen */
88 int point_inside_obi(RayTree *tree, ObjectInstanceRen *obi, float *co)
89 {
90         float maxsize = RE_ray_tree_max_size(tree);
91         Isect isect;
92         float vec[3] = {0.0f,0.0f,1.0f};
93         int final_depth=0, depth=0, limit=20;
94         
95         /* set up the isect */
96         memset(&isect, 0, sizeof(isect));
97         VECCOPY(isect.start, co);
98         isect.end[0] = co[0] + vec[0] * maxsize;
99         isect.end[1] = co[1] + vec[1] * maxsize;
100         isect.end[2] = co[2] + vec[2] * maxsize;
101         
102         /* and give it a little offset to prevent self-intersections */
103         VecMulf(vec, 1e-5);
104         VecAddf(isect.start, isect.start, vec);
105         
106         isect.mode= RE_RAY_MIRROR;
107         isect.face_last= NULL;
108         isect.lay= -1;
109         
110         final_depth = intersect_outside_volume(tree, &isect, vec, limit, depth);
111         
112         /* even number of intersections: point is outside
113          * odd number: point is inside */
114         if (final_depth % 2 == 0) return 0;
115         else return 1;
116 }
117
118 static int inside_check_func(Isect *is, int ob, RayFace *face)
119 {
120         return 1;
121 }
122 static void vlr_face_coords(RayFace *face, float **v1, float **v2, float **v3, float **v4)
123 {
124         VlakRen *vlr= (VlakRen*)face;
125
126         *v1 = (vlr->v1)? vlr->v1->co: NULL;
127         *v2 = (vlr->v2)? vlr->v2->co: NULL;
128         *v3 = (vlr->v3)? vlr->v3->co: NULL;
129         *v4 = (vlr->v4)? vlr->v4->co: NULL;
130 }
131
132 RayTree *create_raytree_obi(ObjectInstanceRen *obi, float *bbmin, float *bbmax)
133 {
134         int v;
135         VlakRen *vlr= NULL;
136         
137         /* create empty raytree */
138         RayTree *tree = RE_ray_tree_create(64, obi->obr->totvlak, bbmin, bbmax,
139                 vlr_face_coords, inside_check_func, NULL, NULL);
140         
141         /* fill it with faces */
142         for(v=0; v<obi->obr->totvlak; v++) {
143                 if((v & 255)==0)
144                         vlr= obi->obr->vlaknodes[v>>8].vlak;
145                 else
146                         vlr++;
147         
148                 RE_ray_tree_add_face(tree, 0, vlr);
149         }
150         
151         RE_ray_tree_done(tree);
152         
153         return tree;
154 }
155
156 /* *** light cache filtering *** */
157
158 static float get_avg_surrounds(float *cache, int *res, int xx, int yy, int zz)
159 {
160         int x, y, z, x_, y_, z_;
161         int added=0;
162         float tot=0.0f;
163         
164         for (z=-1; z <= 1; z++) {
165                 z_ = zz+z;
166                 if (z_ >= 0 && z_ <= res[2]-1) {
167                 
168                         for (y=-1; y <= 1; y++) {
169                                 y_ = yy+y;
170                                 if (y_ >= 0 && y_ <= res[1]-1) {
171                                 
172                                         for (x=-1; x <= 1; x++) {
173                                                 x_ = xx+x;
174                                                 if (x_ >= 0 && x_ <= res[0]-1) {
175                                                 
176                                                         if (cache[ V_I(x_, y_, z_, res) ] > 0.0f) {
177                                                                 tot += cache[ V_I(x_, y_, z_, res) ];
178                                                                 added++;
179                                                         }
180                                                         
181                                                 }
182                                         }
183                                 }
184                         }
185                 }
186         }
187         
188         if (added > 0) tot /= added;
189         
190         return tot;
191 }
192
193 /* function to filter the edges of the light cache, where there was no volume originally.
194  * For each voxel which was originally external to the mesh, it finds the average values of
195  * the surrounding internal voxels and sets the original external voxel to that average amount.
196  * Works almost a bit like a 'dilate' filter */
197 static void lightcache_filter(VolumePrecache *vp)
198 {
199         int x, y, z;
200
201         for (z=0; z < vp->res[2]; z++) {
202                 for (y=0; y < vp->res[1]; y++) {
203                         for (x=0; x < vp->res[0]; x++) {
204                                 /* trigger for outside mesh */
205                                 if (vp->data_r[ V_I(x, y, z, vp->res) ] < -0.f)
206                                         vp->data_r[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_r, vp->res, x, y, z);
207                                 if (vp->data_g[ V_I(x, y, z, vp->res) ] < -0.f)
208                                         vp->data_g[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_g, vp->res, x, y, z);
209                                 if (vp->data_b[ V_I(x, y, z, vp->res) ] < -0.f)
210                                         vp->data_b[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_b, vp->res, x, y, z);
211                         }
212                 }
213         }
214 }
215
216 static void lightcache_filter2(VolumePrecache *vp)
217 {
218         int x, y, z;
219         float *new_r, *new_g, *new_b;
220         int field_size = vp->res[0]*vp->res[1]*vp->res[2]*sizeof(float);
221         
222         new_r = MEM_mallocN(field_size, "temp buffer for light cache filter r channel");
223         new_g = MEM_mallocN(field_size, "temp buffer for light cache filter g channel");
224         new_b = MEM_mallocN(field_size, "temp buffer for light cache filter b channel");
225         
226         memcpy(new_r, vp->data_r, field_size);
227         memcpy(new_g, vp->data_g, field_size);
228         memcpy(new_b, vp->data_b, field_size);
229         
230         for (z=0; z < vp->res[2]; z++) {
231                 for (y=0; y < vp->res[1]; y++) {
232                         for (x=0; x < vp->res[0]; x++) {
233                                 /* trigger for outside mesh */
234                                 if (vp->data_r[ V_I(x, y, z, vp->res) ] < -0.f)
235                                         new_r[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_r, vp->res, x, y, z);
236                                 if (vp->data_g[ V_I(x, y, z, vp->res) ] < -0.f)
237                                         new_g[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_g, vp->res, x, y, z);
238                                 if (vp->data_b[ V_I(x, y, z, vp->res) ] < -0.f)
239                                         new_b[ V_I(x, y, z, vp->res) ] = get_avg_surrounds(vp->data_b, vp->res, x, y, z);
240                         }
241                 }
242         }
243         
244         SWAP(float *, vp->data_r, new_r);
245         SWAP(float *, vp->data_g, new_g);
246         SWAP(float *, vp->data_b, new_b);
247         
248         if (new_r) { MEM_freeN(new_r); new_r=NULL; }
249         if (new_g) { MEM_freeN(new_g); new_g=NULL; }
250         if (new_b) { MEM_freeN(new_b); new_b=NULL; }
251 }
252
253 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
254
255         return z*(n[1]+2)*(n[0]+2) + y*(n[0]+2) + x;
256 }
257
258
259 /* *** multiple scattering approximation *** */
260
261 /* get the total amount of light energy in the light cache. used to normalise after multiple scattering */
262 static float total_ss_energy(VolumePrecache *vp)
263 {
264         int x, y, z;
265         int *res = vp->res;
266         float energy=0.f;
267         
268         for (z=0; z < res[2]; z++) {
269                 for (y=0; y < res[1]; y++) {
270                         for (x=0; x < res[0]; x++) {
271                                 if (vp->data_r[ V_I(x, y, z, res) ] > 0.f) energy += vp->data_r[ V_I(x, y, z, res) ];
272                                 if (vp->data_g[ V_I(x, y, z, res) ] > 0.f) energy += vp->data_g[ V_I(x, y, z, res) ];
273                                 if (vp->data_b[ V_I(x, y, z, res) ] > 0.f) energy += vp->data_b[ V_I(x, y, z, res) ];
274                         }
275                 }
276         }
277         
278         return energy;
279 }
280
281 static float total_ms_energy(float *sr, float *sg, float *sb, int *res)
282 {
283         int x, y, z, i;
284         float energy=0.f;
285         
286         for (z=1;z<=res[2];z++) {
287                 for (y=1;y<=res[1];y++) {
288                         for (x=1;x<=res[0];x++) {
289                         
290                                 i = ms_I(x,y,z,res);
291                                 if (sr[i] > 0.f) energy += sr[i];
292                                 if (sg[i] > 0.f) energy += sg[i];
293                                 if (sb[i] > 0.f) energy += sb[i];
294                         }
295                 }
296         }
297         
298         return energy;
299 }
300
301 static void ms_diffuse(int b, float* x0, float* x, float diff, int *n)
302 {
303         int i, j, k, l;
304         const float dt = VOL_MS_TIMESTEP;
305         const float a = dt*diff*n[0]*n[1]*n[2];
306         
307         for (l=0; l<20; l++)
308         {
309                 for (k=1; k<=n[2]; k++)
310                 {
311                         for (j=1; j<=n[1]; j++)
312                         {
313                                 for (i=1; i<=n[0]; i++)
314                                 {
315                                         x[ms_I(i,j,k,n)] = (x0[ms_I(i,j,k,n)] + a*(
316                                                  x[ms_I(i-1,j,k,n)]+x[ms_I(i+1,j,k,n)]+
317                                                  x[ms_I(i,j-1,k,n)]+x[ms_I(i,j+1,k,n)]+
318                                                  x[ms_I(i,j,k-1,n)]+x[ms_I(i,j,k+1,n)]))/(1+6*a);
319                                 }
320                         }
321                 }
322         }
323 }
324
325 void multiple_scattering_diffusion(Render *re, VolumePrecache *vp, Material *ma)
326 {
327         const float diff = ma->vol.ms_diff * 0.001f;    /* compensate for scaling for a nicer UI range */
328         const float simframes = ma->vol.ms_steps;
329         const int shade_type = ma->vol.shade_type;
330         float fac = ma->vol.ms_intensity;
331         
332         int x, y, z, m;
333         int *n = vp->res;
334         const int size = (n[0]+2)*(n[1]+2)*(n[2]+2);
335         double time, lasttime= PIL_check_seconds_timer();
336         float total;
337         float c=1.0f;
338         int i;
339         float origf;    /* factor for blending in original light cache */
340         float energy_ss, energy_ms;
341
342         float *sr0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
343         float *sr=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
344         float *sg0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
345         float *sg=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
346         float *sb0=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
347         float *sb=(float *)MEM_callocN(size*sizeof(float), "temporary multiple scattering buffer");
348
349         total = (float)(n[0]*n[1]*n[2]*simframes);
350         
351         energy_ss = total_ss_energy(vp);
352         
353         /* Scattering as diffusion pass */
354         for (m=0; m<simframes; m++)
355         {
356                 /* add sources */
357                 for (z=1; z<=n[2]; z++)
358                 {
359                         for (y=1; y<=n[1]; y++)
360                         {
361                                 for (x=1; x<=n[0]; x++)
362                                 {
363                                         i = V_I((x-1), (y-1), (z-1), n);
364                                         time= PIL_check_seconds_timer();
365                                         c++;
366                                                                                 
367                                         if (vp->data_r[i] > 0.f)
368                                                 sr[ms_I(x,y,z,n)] += vp->data_r[i];
369                                         if (vp->data_g[i] > 0.f)
370                                                 sg[ms_I(x,y,z,n)] += vp->data_g[i];
371                                         if (vp->data_b[i] > 0.f)
372                                                 sb[ms_I(x,y,z,n)] += vp->data_b[i];
373                                         
374                                         /* Displays progress every second */
375                                         if(time-lasttime>1.0f) {
376                                                 char str[64];
377                                                 sprintf(str, "Simulating multiple scattering: %d%%", (int)
378                                                                 (100.0f * (c / total)));
379                                                 re->i.infostr= str;
380                                                 re->stats_draw(re->sdh, &re->i);
381                                                 re->i.infostr= NULL;
382                                                 lasttime= time;
383                                         }
384                                 }
385                         }
386                 }
387                 SWAP(float *, sr, sr0);
388                 SWAP(float *, sg, sg0);
389                 SWAP(float *, sb, sb0);
390
391                 /* main diffusion simulation */
392                 ms_diffuse(0, sr0, sr, diff, n);
393                 ms_diffuse(0, sg0, sg, diff, n);
394                 ms_diffuse(0, sb0, sb, diff, n);
395                 
396                 if (re->test_break(re->tbh)) break;
397         }
398         
399         /* normalisation factor to conserve energy */
400         energy_ms = total_ms_energy(sr, sg, sb, n);
401         fac *= (energy_ss / energy_ms);
402         
403         /* blend multiple scattering back in the light cache */
404         if (shade_type == MA_VOL_SHADE_SINGLEPLUSMULTIPLE) {
405                 /* conserve energy - half single, half multiple */
406                 origf = 0.5f;
407                 fac *= 0.5f;
408         } else {
409                 origf = 0.0f;
410         }
411
412         for (z=1;z<=n[2];z++)
413         {
414                 for (y=1;y<=n[1];y++)
415                 {
416                         for (x=1;x<=n[0];x++)
417                         {
418                                 int index=(x-1)*n[1]*n[2] + (y-1)*n[2] + z-1;
419                                 vp->data_r[index] = origf * vp->data_r[index] + fac * sr[ms_I(x,y,z,n)];
420                                 vp->data_g[index] = origf * vp->data_g[index] + fac * sg[ms_I(x,y,z,n)];
421                                 vp->data_b[index] = origf * vp->data_b[index] + fac * sb[ms_I(x,y,z,n)];
422                         }
423                 }
424         }
425
426         MEM_freeN(sr0);
427         MEM_freeN(sr);
428         MEM_freeN(sg0);
429         MEM_freeN(sg);
430         MEM_freeN(sb0);
431         MEM_freeN(sb);
432 }
433
434
435
436 #if 0 // debug stuff
437 static void *vol_precache_part_test(void *data)
438 {
439         VolPrecachePart *pa = data;
440
441         printf("part number: %d \n", pa->num);
442         printf("done: %d \n", pa->done);
443         printf("x min: %d   x max: %d \n", pa->minx, pa->maxx);
444         printf("y min: %d   y max: %d \n", pa->miny, pa->maxy);
445         printf("z min: %d   z max: %d \n", pa->minz, pa->maxz);
446
447         return NULL;
448 }
449 #endif
450
451 /* Iterate over the 3d voxel grid, and fill the voxels with scattering information
452  *
453  * It's stored in memory as 3 big float grids next to each other, one for each RGB channel.
454  * I'm guessing the memory alignment may work out better this way for the purposes
455  * of doing linear interpolation, but I haven't actually tested this theory! :)
456  */
457 static void *vol_precache_part(void *data)
458 {
459         VolPrecachePart *pa =  (VolPrecachePart *)data;
460         ObjectInstanceRen *obi = pa->obi;
461         RayTree *tree = pa->tree;
462         ShadeInput *shi = pa->shi;
463         float scatter_col[3] = {0.f, 0.f, 0.f};
464         float co[3];
465         int x, y, z;
466         const int res[3]= {pa->res[0], pa->res[1], pa->res[2]};
467
468         for (z= pa->minz; z < pa->maxz; z++) {
469                 co[2] = pa->bbmin[2] + (pa->voxel[2] * (z + 0.5f));
470                 
471                 for (y= pa->miny; y < pa->maxy; y++) {
472                         co[1] = pa->bbmin[1] + (pa->voxel[1] * (y + 0.5f));
473                         
474                         for (x=pa->minx; x < pa->maxx; x++) {
475                                 co[0] = pa->bbmin[0] + (pa->voxel[0] * (x + 0.5f));
476                                 
477                                 // don't bother if the point is not inside the volume mesh
478                                 if (!point_inside_obi(tree, obi, co)) {
479                                         obi->volume_precache->data_r[ V_I(x, y, z, res) ] = -1.0f;
480                                         obi->volume_precache->data_g[ V_I(x, y, z, res) ] = -1.0f;
481                                         obi->volume_precache->data_b[ V_I(x, y, z, res) ] = -1.0f;
482                                         continue;
483                                 }
484                                 
485                                 VecCopyf(shi->view, co);
486                                 Normalize(shi->view);
487                                 vol_get_scattering(shi, scatter_col, co);
488                         
489                                 obi->volume_precache->data_r[ V_I(x, y, z, res) ] = scatter_col[0];
490                                 obi->volume_precache->data_g[ V_I(x, y, z, res) ] = scatter_col[1];
491                                 obi->volume_precache->data_b[ V_I(x, y, z, res) ] = scatter_col[2];
492                         }
493                 }
494         }
495         
496         pa->done = 1;
497         
498         return 0;
499 }
500
501
502 static void precache_setup_shadeinput(Render *re, ObjectInstanceRen *obi, Material *ma, ShadeInput *shi)
503 {
504         memset(shi, 0, sizeof(ShadeInput)); 
505         shi->depth= 1;
506         shi->mask= 1;
507         shi->mat = ma;
508         shi->vlr = NULL;
509         memcpy(&shi->r, &shi->mat->r, 23*sizeof(float));        // note, keep this synced with render_types.h
510         shi->har= shi->mat->har;
511         shi->obi= obi;
512         shi->obr= obi->obr;
513         shi->lay = re->scene->lay;
514 }
515
516 static void precache_init_parts(Render *re, RayTree *tree, ShadeInput *shi, ObjectInstanceRen *obi, int totthread, int *parts)
517 {
518         VolumePrecache *vp = obi->volume_precache;
519         int i=0, x, y, z;
520         float voxel[3];
521         int sizex, sizey, sizez;
522         float *bbmin=obi->obr->boundbox[0], *bbmax=obi->obr->boundbox[1];
523         int *res;
524         int minx, maxx;
525         int miny, maxy;
526         int minz, maxz;
527         
528         if (!vp) return;
529
530         BLI_freelistN(&re->volume_precache_parts);
531         
532         /* currently we just subdivide the box, number of threads per side */
533         parts[0] = parts[1] = parts[2] = totthread;
534         res = vp->res;
535         
536         VecSubf(voxel, bbmax, bbmin);
537         
538         voxel[0] /= res[0];
539         voxel[1] /= res[1];
540         voxel[2] /= res[2];
541
542         for (x=0; x < parts[0]; x++) {
543                 sizex = ceil(res[0] / (float)parts[0]);
544                 minx = x * sizex;
545                 maxx = minx + sizex;
546                 maxx = (maxx>res[0])?res[0]:maxx;
547                 
548                 for (y=0; y < parts[1]; y++) {
549                         sizey = ceil(res[1] / (float)parts[1]);
550                         miny = y * sizey;
551                         maxy = miny + sizey;
552                         maxy = (maxy>res[1])?res[1]:maxy;
553                         
554                         for (z=0; z < parts[2]; z++) {
555                                 VolPrecachePart *pa= MEM_callocN(sizeof(VolPrecachePart), "new precache part");
556                                 
557                                 sizez = ceil(res[2] / (float)parts[2]);
558                                 minz = z * sizez;
559                                 maxz = minz + sizez;
560                                 maxz = (maxz>res[2])?res[2]:maxz;
561                                                 
562                                 pa->done = 0;
563                                 pa->working = 0;
564                                 
565                                 pa->num = i;
566                                 pa->tree = tree;
567                                 pa->shi = shi;
568                                 pa->obi = obi;
569                                 VECCOPY(pa->bbmin, bbmin);
570                                 VECCOPY(pa->voxel, voxel);
571                                 VECCOPY(pa->res, res);
572                                 
573                                 pa->minx = minx; pa->maxx = maxx;
574                                 pa->miny = miny; pa->maxy = maxy;
575                                 pa->minz = minz; pa->maxz = maxz;
576                                 
577                                 
578                                 BLI_addtail(&re->volume_precache_parts, pa);
579                                 
580                                 i++;
581                         }
582                 }
583         }
584 }
585
586 static VolPrecachePart *precache_get_new_part(Render *re)
587 {
588         VolPrecachePart *pa, *nextpa=NULL;
589         
590         for (pa = re->volume_precache_parts.first; pa; pa=pa->next)
591         {
592                 if (pa->done==0 && pa->working==0) {
593                         nextpa = pa;
594                         break;
595                 }
596         }
597
598         return nextpa;
599 }
600
601 static int precache_resolution(VolumePrecache *vp, float *bbmin, float *bbmax, int res)
602 {
603         float dim[3], div;
604         
605         VecSubf(dim, bbmax, bbmin);
606         
607         div = MAX3(dim[0], dim[1], dim[2]);
608         dim[0] /= div;
609         dim[1] /= div;
610         dim[2] /= div;
611         
612         vp->res[0] = dim[0] * (float)res;
613         vp->res[1] = dim[1] * (float)res;
614         vp->res[2] = dim[2] * (float)res;
615         
616         if ((vp->res[0] < 1) || (vp->res[1] < 1) || (vp->res[2] < 1))
617                 return 0;
618         
619         return 1;
620 }
621
622 /* Precache a volume into a 3D voxel grid.
623  * The voxel grid is stored in the ObjectInstanceRen, 
624  * in camera space, aligned with the ObjectRen's bounding box.
625  * Resolution is defined by the user.
626  */
627 void vol_precache_objectinstance_threads(Render *re, ObjectInstanceRen *obi, Material *ma)
628 {
629         VolumePrecache *vp;
630         VolPrecachePart *nextpa, *pa;
631         RayTree *tree;
632         ShadeInput shi;
633         ListBase threads;
634         float *bbmin=obi->obr->boundbox[0], *bbmax=obi->obr->boundbox[1];
635         int parts[3], totparts;
636         
637         int caching=1, counter=0;
638         int totthread = re->r.threads;
639         
640         double time, lasttime= PIL_check_seconds_timer();
641         
642         R = *re;
643
644         /* create a raytree with just the faces of the instanced ObjectRen, 
645          * used for checking if the cached point is inside or outside. */
646         tree = create_raytree_obi(obi, bbmin, bbmax);
647         if (!tree) return;
648
649         vp = MEM_callocN(sizeof(VolumePrecache), "volume light cache");
650         
651         if (!precache_resolution(vp, bbmin, bbmax, ma->vol.precache_resolution)) {
652                 MEM_freeN(vp);
653                 vp = NULL;
654                 return;
655         }
656
657         vp->data_r = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data red channel");
658         vp->data_g = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data green channel");
659         vp->data_b = MEM_callocN(sizeof(float)*vp->res[0]*vp->res[1]*vp->res[2], "volume light cache data blue channel");
660         obi->volume_precache = vp;
661
662         /* Need a shadeinput to calculate scattering */
663         precache_setup_shadeinput(re, obi, ma, &shi);
664         
665         precache_init_parts(re, tree, &shi, obi, totthread, parts);
666         totparts = parts[0] * parts[1] * parts[2];
667         
668         BLI_init_threads(&threads, vol_precache_part, totthread);
669         
670         while(caching) {
671
672                 if(BLI_available_threads(&threads) && !(re->test_break(re->tbh))) {
673                         nextpa = precache_get_new_part(re);
674                         if (nextpa) {
675                                 nextpa->working = 1;
676                                 BLI_insert_thread(&threads, nextpa);
677                         }
678                 }
679                 else PIL_sleep_ms(50);
680
681                 caching=0;
682                 counter=0;
683                 for(pa= re->volume_precache_parts.first; pa; pa= pa->next) {
684                         
685                         if(pa->done) {
686                                 counter++;
687                                 BLI_remove_thread(&threads, pa);
688                         } else
689                                 caching = 1;
690                 }
691                 
692                 if (re->test_break(re->tbh) && BLI_available_threads(&threads)==totthread)
693                         caching=0;
694                 
695                 time= PIL_check_seconds_timer();
696                 if(time-lasttime>1.0f) {
697                         char str[64];
698                         sprintf(str, "Precaching volume: %d%%", (int)(100.0f * ((float)counter / (float)totparts)));
699                         re->i.infostr= str;
700                         re->stats_draw(re->sdh, &re->i);
701                         re->i.infostr= NULL;
702                         lasttime= time;
703                 }
704         }
705         
706         BLI_end_threads(&threads);
707         BLI_freelistN(&re->volume_precache_parts);
708         
709         if(tree) {
710                 RE_ray_tree_free(tree);
711                 tree= NULL;
712         }
713         
714         lightcache_filter(obi->volume_precache);
715         
716         if (ELEM(ma->vol.shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SINGLEPLUSMULTIPLE))
717         {
718                 multiple_scattering_diffusion(re, vp, ma);
719         }
720 }
721
722 static int using_lightcache(Material *ma)
723 {
724         return (((ma->vol.shadeflag & MA_VOL_PRECACHESHADING) && (ma->vol.shade_type == MA_VOL_SHADE_SINGLE))
725                 || (ELEM(ma->vol.shade_type, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SINGLEPLUSMULTIPLE)));
726 }
727
728 /* loop through all objects (and their associated materials)
729  * marked for pre-caching in convertblender.c, and pre-cache them */
730 void volume_precache(Render *re)
731 {
732         ObjectInstanceRen *obi;
733         VolumeOb *vo;
734
735         for(vo= re->volumes.first; vo; vo= vo->next) {
736                 if (using_lightcache(vo->ma)) {
737                         for(obi= re->instancetable.first; obi; obi= obi->next) {
738                                 if (obi->obr == vo->obr) {
739                                         vol_precache_objectinstance_threads(re, obi, vo->ma);
740                                 }
741                         }
742                 }
743         }
744         
745         re->i.infostr= NULL;
746         re->stats_draw(re->sdh, &re->i);
747 }
748
749 void free_volume_precache(Render *re)
750 {
751         ObjectInstanceRen *obi;
752         
753         for(obi= re->instancetable.first; obi; obi= obi->next) {
754                 if (obi->volume_precache != NULL) {
755                         MEM_freeN(obi->volume_precache->data_r);
756                         MEM_freeN(obi->volume_precache->data_g);
757                         MEM_freeN(obi->volume_precache->data_b);
758                         MEM_freeN(obi->volume_precache);
759                         obi->volume_precache = NULL;
760                 }
761         }
762         
763         BLI_freelistN(&re->volumes);
764 }
765
766 int point_inside_volume_objectinstance(ObjectInstanceRen *obi, float *co)
767 {
768         RayTree *tree;
769         int inside=0;
770         
771         tree = create_raytree_obi(obi, obi->obr->boundbox[0], obi->obr->boundbox[1]);
772         if (!tree) return 0;
773         
774         inside = point_inside_obi(tree, obi, co);
775         
776         RE_ray_tree_free(tree);
777         tree= NULL;
778         
779         return inside;
780 }
781