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