Merge with trunk r41342
[blender-staging.git] / source / blender / render / intern / source / volumetric.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Matt Ebb, Raul Fernandez Hernandez (Farsthary)
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/render/intern/source/volumetric.c
29  *  \ingroup render
30  */
31
32
33 #include <math.h>
34 #include <stdlib.h>
35 #include <string.h>
36 #include <float.h>
37
38 #include "MEM_guardedalloc.h"
39
40 #include "BLI_blenlib.h"
41 #include "BLI_math.h"
42 #include "BLI_rand.h"
43 #include "BLI_voxel.h"
44 #include "BLI_utildefines.h"
45
46 #include "RE_shader_ext.h"
47
48 #include "DNA_material_types.h"
49 #include "DNA_group_types.h"
50 #include "DNA_lamp_types.h"
51 #include "DNA_meta_types.h"
52
53 #include "BKE_global.h"
54
55 #include "render_types.h"
56 #include "pixelshading.h"
57 #include "rayintersection.h"
58 #include "rayobject.h"
59 #include "shading.h"
60 #include "shadbuf.h"
61 #include "texture.h"
62 #include "volumetric.h"
63 #include "volume_precache.h"
64
65 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
66 /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
67 /* only to be used here in this file, it's for speed */
68 extern struct Render R;
69 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
70
71 /* luminance rec. 709 */
72 BM_INLINE float luminance(const float col[3])
73 {
74         return (0.212671f*col[0] + 0.71516f*col[1] + 0.072169f*col[2]);
75 }
76
77 /* tracing */
78 static float vol_get_shadow(ShadeInput *shi, LampRen *lar, const float co[3])
79 {
80         float visibility = 1.f;
81         
82         if(lar->shb) {
83                 float dxco[3]={0.f, 0.f, 0.f}, dyco[3]={0.f, 0.f, 0.f};
84                 
85                 visibility = testshadowbuf(&R, lar->shb, co, dxco, dyco, 1.0, 0.0);             
86         } else if (lar->mode & LA_SHAD_RAY) {
87                 /* trace shadow manually, no good lamp api atm */
88                 Isect is;
89                 
90                 copy_v3_v3(is.start, co);
91                 if(lar->type==LA_SUN || lar->type==LA_HEMI) {
92                         is.dir[0] = -lar->vec[0];
93                         is.dir[1] = -lar->vec[1];
94                         is.dir[2] = -lar->vec[2];
95                         is.dist = R.maxdist;
96                 } else {
97                         sub_v3_v3v3(is.dir, lar->co, is.start);
98                         is.dist = normalize_v3( is.dir );
99                 }
100
101                 is.mode = RE_RAY_MIRROR;
102                 is.check = RE_CHECK_VLR_NON_SOLID_MATERIAL;
103                 is.skip = 0;
104                 
105                 if(lar->mode & (LA_LAYER|LA_LAYER_SHADOW))
106                         is.lay= lar->lay;       
107                 else
108                         is.lay= -1;
109                         
110                 is.orig.ob = NULL;
111                 is.orig.face = NULL;
112                 is.last_hit = lar->last_hit[shi->thread];
113                 
114                 if(RE_rayobject_raycast(R.raytree,&is)) {
115                         visibility = 0.f;
116                 }
117                 
118                 lar->last_hit[shi->thread]= is.last_hit;
119         }
120         return visibility;
121 }
122
123 static int vol_get_bounds(ShadeInput *shi, const float co[3], const float vec[3], float hitco[3], Isect *isect, int intersect_type)
124 {
125         
126         copy_v3_v3(isect->start, co);
127         copy_v3_v3(isect->dir, vec);
128         isect->dist = FLT_MAX;
129         isect->mode= RE_RAY_MIRROR;
130         isect->last_hit = NULL;
131         isect->lay= -1;
132         isect->check= RE_CHECK_VLR_NONE;
133         
134         if (intersect_type == VOL_BOUNDS_DEPTH) {
135                 isect->skip = RE_SKIP_VLR_NEIGHBOUR;
136                 isect->orig.face = (void*)shi->vlr;
137                 isect->orig.ob = (void*)shi->obi;
138         } else { // if (intersect_type == VOL_BOUNDS_SS) {
139                 isect->skip= 0;
140                 isect->orig.face= NULL;
141                 isect->orig.ob = NULL;
142         }
143         
144         if(RE_rayobject_raycast(R.raytree, isect))
145         {
146                 hitco[0] = isect->start[0] + isect->dist*isect->dir[0];
147                 hitco[1] = isect->start[1] + isect->dist*isect->dir[1];
148                 hitco[2] = isect->start[2] + isect->dist*isect->dir[2];
149                 return 1;
150         } else {
151                 return 0;
152         }
153 }
154
155 static void shade_intersection(ShadeInput *shi, float col_r[4], Isect *is)
156 {
157         ShadeInput shi_new;
158         ShadeResult shr_new;
159         
160         memset(&shi_new, 0, sizeof(ShadeInput)); 
161         
162         shi_new.mask= shi->mask;
163         shi_new.osatex= shi->osatex;
164         shi_new.thread= shi->thread;
165         shi_new.depth = shi->depth + 1;
166         shi_new.volume_depth= shi->volume_depth + 1;
167         shi_new.xs= shi->xs;
168         shi_new.ys= shi->ys;
169         shi_new.lay= shi->lay;
170         shi_new.passflag= SCE_PASS_COMBINED; /* result of tracing needs no pass info */
171         shi_new.combinedflag= 0xFFFFFF;          /* ray trace does all options */
172         shi_new.light_override= shi->light_override;
173         shi_new.mat_override= shi->mat_override;
174         
175         copy_v3_v3(shi_new.camera_co, is->start);
176         
177         memset(&shr_new, 0, sizeof(ShadeResult));
178         
179         /* hardcoded limit of 100 for now - prevents problems in weird geometry */
180         if (shi->volume_depth < 100) {
181                 shade_ray(is, &shi_new, &shr_new);
182         }
183         
184         copy_v3_v3(col_r, shr_new.combined);
185         col_r[3] = shr_new.alpha;
186 }
187
188 static void vol_trace_behind(ShadeInput *shi, VlakRen *vlr, const float co[3], float col_r[4])
189 {
190         Isect isect;
191         
192         copy_v3_v3(isect.start, co);
193         copy_v3_v3(isect.dir, shi->view);
194         isect.dist = FLT_MAX;
195         
196         isect.mode= RE_RAY_MIRROR;
197         isect.check = RE_CHECK_VLR_NONE;
198         isect.skip = RE_SKIP_VLR_NEIGHBOUR;
199         isect.orig.ob = (void*) shi->obi;
200         isect.orig.face = (void*)vlr;
201         isect.last_hit = NULL;
202         isect.lay= -1;
203         
204         /* check to see if there's anything behind the volume, otherwise shade the sky */
205         if(RE_rayobject_raycast(R.raytree, &isect)) {
206                 shade_intersection(shi, col_r, &isect);
207         } else {
208                 shadeSkyView(col_r, co, shi->view, NULL, shi->thread);
209                 shadeSunView(col_r, shi->view);
210         } 
211 }
212
213
214 /* trilinear interpolation */
215 static void vol_get_precached_scattering(Render *re, ShadeInput *shi, float scatter_col[3], const float co[3])
216 {
217         VolumePrecache *vp = shi->obi->volume_precache;
218         float bbmin[3], bbmax[3], dim[3];
219         float world_co[3], sample_co[3];
220         
221         if (!vp) return;
222         
223         /* find sample point in global space bounding box 0.0-1.0 */
224         global_bounds_obi(re, shi->obi, bbmin, bbmax);
225         sub_v3_v3v3(dim, bbmax, bbmin);
226         mul_v3_m4v3(world_co, re->viewinv, co); 
227
228         /* sample_co in 0.0-1.0 */
229         sample_co[0] = (world_co[0] - bbmin[0]) / dim[0];
230         sample_co[1] = (world_co[1] - bbmin[1]) / dim[1];
231         sample_co[2] = (world_co[2] - bbmin[2]) / dim[2];
232
233         scatter_col[0] = voxel_sample_triquadratic(vp->data_r, vp->res, sample_co);
234         scatter_col[1] = voxel_sample_triquadratic(vp->data_g, vp->res, sample_co);
235         scatter_col[2] = voxel_sample_triquadratic(vp->data_b, vp->res, sample_co);
236 }
237
238 /* Meta object density, brute force for now 
239  * (might be good enough anyway, don't need huge number of metaobs to model volumetric objects */
240 static float metadensity(Object* ob, const float co[3])
241 {
242         float mat[4][4], imat[4][4], dens = 0.f;
243         MetaBall* mb = (MetaBall*)ob->data;
244         MetaElem* ml;
245         
246         /* transform co to meta-element */
247         float tco[3] = {co[0], co[1], co[2]};
248         mul_m4_m4m4(mat, ob->obmat, R.viewmat);
249         invert_m4_m4(imat, mat);
250         mul_m4_v3(imat, tco);
251         
252         for (ml = mb->elems.first; ml; ml=ml->next) {
253                 float bmat[3][3], dist2;
254                 
255                 /* element rotation transform */
256                 float tp[3] = {ml->x - tco[0], ml->y - tco[1], ml->z - tco[2]};
257                 quat_to_mat3( bmat,ml->quat);
258                 transpose_m3(bmat);     // rot.only, so inverse == transpose
259                 mul_m3_v3(bmat, tp);
260                 
261                 /* MB_BALL default */
262                 switch (ml->type) {
263                         case MB_ELIPSOID:
264                                 tp[0] /= ml->expx, tp[1] /= ml->expy, tp[2] /= ml->expz;
265                                 break;
266                         case MB_CUBE:
267                                 tp[2] = (tp[2] > ml->expz) ? (tp[2] - ml->expz) : ((tp[2] < -ml->expz) ? (tp[2] + ml->expz) : 0.f);
268                                 // no break, xy as plane
269                         case MB_PLANE:
270                                 tp[1] = (tp[1] > ml->expy) ? (tp[1] - ml->expy) : ((tp[1] < -ml->expy) ? (tp[1] + ml->expy) : 0.f);
271                                 // no break, x as tube
272                         case MB_TUBE:
273                                 tp[0] = (tp[0] > ml->expx) ? (tp[0] - ml->expx) : ((tp[0] < -ml->expx) ? (tp[0] + ml->expx) : 0.f);
274                 }
275                 
276                 /* ml->rad2 is not set */
277                 dist2 = 1.f - ((tp[0]*tp[0] + tp[1]*tp[1] + tp[2]*tp[2]) / (ml->rad*ml->rad));
278                 if (dist2 > 0.f)
279                         dens += (ml->flag & MB_NEGATIVE) ? -ml->s*dist2*dist2*dist2 : ml->s*dist2*dist2*dist2;
280         }
281         
282         dens -= mb->thresh;
283         return (dens < 0.f) ? 0.f : dens;
284 }
285
286 float vol_get_density(struct ShadeInput *shi, const float co[3])
287 {
288         float density = shi->mat->vol.density;
289         float density_scale = shi->mat->vol.density_scale;
290                 
291         if (shi->mat->mapto_textured & MAP_DENSITY)
292                 do_volume_tex(shi, co, MAP_DENSITY, NULL, &density, &R);
293         
294         // if meta-object, modulate by metadensity without increasing it
295         if (shi->obi->obr->ob->type == OB_MBALL) {
296                 const float md = metadensity(shi->obi->obr->ob, co);
297                 if (md < 1.f) density *= md;
298          }
299         
300         return density * density_scale;
301 }
302
303
304 /* Color of light that gets scattered out by the volume */
305 /* Uses same physically based scattering parameter as in transmission calculations, 
306  * along with artificial reflection scale/reflection color tint */
307 static void vol_get_reflection_color(ShadeInput *shi, float ref_col[3], const float co[3])
308 {
309         float scatter = shi->mat->vol.scattering;
310         float reflection= shi->mat->vol.reflection;
311         copy_v3_v3(ref_col, shi->mat->vol.reflection_col);
312         
313         if (shi->mat->mapto_textured & (MAP_SCATTERING+MAP_REFLECTION_COL))
314                 do_volume_tex(shi, co, MAP_SCATTERING+MAP_REFLECTION_COL, ref_col, &scatter, &R);
315         
316         /* only one single float parameter at a time... :s */
317         if (shi->mat->mapto_textured & (MAP_REFLECTION))
318                 do_volume_tex(shi, co, MAP_REFLECTION, NULL, &reflection, &R);
319         
320         ref_col[0] = reflection * ref_col[0] * scatter;
321         ref_col[1] = reflection * ref_col[1] * scatter;
322         ref_col[2] = reflection * ref_col[2] * scatter;
323 }
324
325 /* compute emission component, amount of radiance to add per segment
326  * can be textured with 'emit' */
327 static void vol_get_emission(ShadeInput *shi, float emission_col[3], const float co[3])
328 {
329         float emission = shi->mat->vol.emission;
330         copy_v3_v3(emission_col, shi->mat->vol.emission_col);
331         
332         if (shi->mat->mapto_textured & (MAP_EMISSION+MAP_EMISSION_COL))
333                 do_volume_tex(shi, co, MAP_EMISSION+MAP_EMISSION_COL, emission_col, &emission, &R);
334         
335         emission_col[0] = emission_col[0] * emission;
336         emission_col[1] = emission_col[1] * emission;
337         emission_col[2] = emission_col[2] * emission;
338 }
339
340
341 /* A combination of scattering and absorption -> known as sigma T.
342  * This can possibly use a specific scattering color, 
343  * and absorption multiplier factor too, but these parameters are left out for simplicity.
344  * It's easy enough to get a good wide range of results with just these two parameters. */
345 static void vol_get_sigma_t(ShadeInput *shi, float sigma_t[3], const float co[3])
346 {
347         /* technically absorption, but named transmission color 
348          * since it describes the effect of the coloring *after* absorption */
349         float transmission_col[3] = {shi->mat->vol.transmission_col[0], shi->mat->vol.transmission_col[1], shi->mat->vol.transmission_col[2]};
350         float scattering = shi->mat->vol.scattering;
351         
352         if (shi->mat->mapto_textured & (MAP_SCATTERING+MAP_TRANSMISSION_COL))
353                 do_volume_tex(shi, co, MAP_SCATTERING+MAP_TRANSMISSION_COL, transmission_col, &scattering, &R);
354         
355         sigma_t[0] = (1.0f - transmission_col[0]) + scattering;
356         sigma_t[1] = (1.0f - transmission_col[1]) + scattering;
357         sigma_t[2] = (1.0f - transmission_col[2]) + scattering;
358 }
359
360 /* phase function - determines in which directions the light 
361  * is scattered in the volume relative to incoming direction 
362  * and view direction */
363 static float vol_get_phasefunc(ShadeInput *UNUSED(shi), float g, const float w[3], const float wp[3])
364 {
365         const float normalize = 0.25f; // = 1.f/4.f = M_PI/(4.f*M_PI)
366         
367         /* normalization constant is 1/4 rather than 1/4pi, since
368          * Blender's shading system doesn't normalise for
369          * energy conservation - eg. multiplying by pdf ( 1/pi for a lambert brdf ).
370          * This means that lambert surfaces in Blender are pi times brighter than they 'should be'
371          * and therefore, with correct energy conservation, volumes will darker than other solid objects,
372          * for the same lighting intensity.
373          * To correct this, scale up the phase function values by pi
374          * until Blender's shading system supports this better. --matt
375          */
376         
377         if (g == 0.f) { /* isotropic */
378                 return normalize * 1.f;
379         } else {                /* schlick */
380                 const float k = 1.55f * g - .55f * g * g * g;
381                 const float kcostheta = k * dot_v3v3(w, wp);
382                 return normalize * (1.f - k*k) / ((1.f - kcostheta) * (1.f - kcostheta));
383         }
384         
385         /*
386          * not used, but here for reference:
387         switch (phasefunc_type) {
388                 case MA_VOL_PH_MIEHAZY:
389                         return normalize * (0.5f + 4.5f * powf(0.5 * (1.f + costheta), 8.f));
390                 case MA_VOL_PH_MIEMURKY:
391                         return normalize * (0.5f + 16.5f * powf(0.5 * (1.f + costheta), 32.f));
392                 case MA_VOL_PH_RAYLEIGH:
393                         return normalize * 3.f/4.f * (1 + costheta * costheta);
394                 case MA_VOL_PH_HG:
395                         return normalize * (1.f - g*g) / powf(1.f + g*g - 2.f * g * costheta, 1.5f));
396                 case MA_VOL_PH_SCHLICK:
397                 {
398                         const float k = 1.55f * g - .55f * g * g * g;
399                         const float kcostheta = k * costheta;
400                         return normalize * (1.f - k*k) / ((1.f - kcostheta) * (1.f - kcostheta));
401                 }
402                 case MA_VOL_PH_ISOTROPIC:
403                 default:
404                         return normalize * 1.f;
405         }
406         */
407 }
408
409 /* Compute transmittance = e^(-attenuation) */
410 static void vol_get_transmittance_seg(ShadeInput *shi, float tr[3], float stepsize, const float co[3], float density)
411 {
412         /* input density = density at co */
413         float tau[3] = {0.f, 0.f, 0.f};
414         const float stepd = density * stepsize;
415         float sigma_t[3];
416         
417         vol_get_sigma_t(shi, sigma_t, co);
418         
419         /* homogenous volume within the sampled distance */
420         tau[0] += stepd * sigma_t[0];
421         tau[1] += stepd * sigma_t[1];
422         tau[2] += stepd * sigma_t[2];
423         
424         tr[0] *= expf(-tau[0]);
425         tr[1] *= expf(-tau[1]);
426         tr[2] *= expf(-tau[2]);
427 }
428
429 /* Compute transmittance = e^(-attenuation) */
430 static void vol_get_transmittance(ShadeInput *shi, float tr[3], const float co[3], const float endco[3])
431 {
432         float p[3] = {co[0], co[1], co[2]};
433         float step_vec[3] = {endco[0] - co[0], endco[1] - co[1], endco[2] - co[2]};
434         float tau[3] = {0.f, 0.f, 0.f};
435
436         float t0 = 0.f;
437         float t1 = normalize_v3(step_vec);
438         float pt0 = t0;
439         
440         t0 += shi->mat->vol.stepsize * ((shi->mat->vol.stepsize_type == MA_VOL_STEP_CONSTANT) ? 0.5f : BLI_thread_frand(shi->thread));
441         p[0] += t0 * step_vec[0];
442         p[1] += t0 * step_vec[1];
443         p[2] += t0 * step_vec[2];
444         mul_v3_fl(step_vec, shi->mat->vol.stepsize);
445
446         for (; t0 < t1; pt0 = t0, t0 += shi->mat->vol.stepsize) {
447                 const float d = vol_get_density(shi, p);
448                 const float stepd = (t0 - pt0) * d;
449                 float sigma_t[3];
450                 
451                 vol_get_sigma_t(shi, sigma_t, p);
452                 
453                 tau[0] += stepd * sigma_t[0];
454                 tau[1] += stepd * sigma_t[1];
455                 tau[2] += stepd * sigma_t[2];
456                 
457                 add_v3_v3(p, step_vec);
458         }
459         
460         /* return transmittance */
461         tr[0] = expf(-tau[0]);
462         tr[1] = expf(-tau[1]);
463         tr[2] = expf(-tau[2]);
464 }
465
466 static void vol_shade_one_lamp(struct ShadeInput *shi, const float co[3], const float view[3], LampRen *lar, float lacol[3])
467 {
468         float visifac, lv[3], lampdist;
469         float tr[3]={1.0,1.0,1.0};
470         float hitco[3], *atten_co;
471         float p, ref_col[3];
472         
473         if (lar->mode & LA_LAYER) if((lar->lay & shi->obi->lay)==0) return;
474         if ((lar->lay & shi->lay)==0) return;
475         if (lar->energy == 0.0f) return;
476         
477         if ((visifac= lamp_get_visibility(lar, co, lv, &lampdist)) == 0.f) return;
478         
479         copy_v3_v3(lacol, &lar->r);
480         
481         if(lar->mode & LA_TEXTURE) {
482                 shi->osatex= 0;
483                 do_lamp_tex(lar, lv, shi, lacol, LA_TEXTURE);
484         }
485
486         mul_v3_fl(lacol, visifac);
487
488         if (ELEM(lar->type, LA_SUN, LA_HEMI))
489                 copy_v3_v3(lv, lar->vec);
490         negate_v3(lv);
491         
492         if (shi->mat->vol.shade_type == MA_VOL_SHADE_SHADOWED) {
493                 mul_v3_fl(lacol, vol_get_shadow(shi, lar, co));
494         }
495         else if (ELEM3(shi->mat->vol.shade_type, MA_VOL_SHADE_SHADED, MA_VOL_SHADE_MULTIPLE, MA_VOL_SHADE_SHADEDPLUSMULTIPLE))
496         {
497                 Isect is;
498                 
499                 if (shi->mat->vol.shadeflag & MA_VOL_RECV_EXT_SHADOW) {
500                         mul_v3_fl(lacol, vol_get_shadow(shi, lar, co));
501                         if (luminance(lacol) < 0.001f) return;
502                 }
503                 
504                 /* find minimum of volume bounds, or lamp coord */
505                 if (vol_get_bounds(shi, co, lv, hitco, &is, VOL_BOUNDS_SS)) {
506                         float dist = len_v3v3(co, hitco);
507                         VlakRen *vlr = (VlakRen *)is.hit.face;
508                         
509                         /* simple internal shadowing */
510                         if (vlr->mat->material_type == MA_TYPE_SURFACE) {
511                                 lacol[0] = lacol[1] = lacol[2] = 0.0f;
512                                 return;
513                         }
514
515                         if (ELEM(lar->type, LA_SUN, LA_HEMI))
516                                 /* infinite lights, can never be inside volume */
517                                 atten_co = hitco;
518                         else if ( lampdist < dist ) {
519                                 atten_co = lar->co;
520                         } else
521                                 atten_co = hitco;
522                         
523                         vol_get_transmittance(shi, tr, co, atten_co);
524                         
525                         mul_v3_v3v3(lacol, lacol, tr);
526                 }
527                 else {
528                         /* Point is on the outside edge of the volume,
529                          * therefore no attenuation, full transmission.
530                          * Radiance from lamp remains unchanged */
531                 }
532         }
533         
534         if (luminance(lacol) < 0.001f) return;
535         
536         normalize_v3(lv);
537         p = vol_get_phasefunc(shi, shi->mat->vol.asymmetry, view, lv);
538         
539         /* physically based scattering with non-physically based RGB gain */
540         vol_get_reflection_color(shi, ref_col, co);
541         
542         lacol[0] *= p * ref_col[0];
543         lacol[1] *= p * ref_col[1];
544         lacol[2] *= p * ref_col[2];
545 }
546
547 /* single scattering only for now */
548 void vol_get_scattering(ShadeInput *shi, float scatter_col[3], const float co[3], const float view[3])
549 {
550         ListBase *lights;
551         GroupObject *go;
552         LampRen *lar;
553
554         zero_v3(scatter_col);
555
556         lights= get_lights(shi);
557         for(go=lights->first; go; go= go->next)
558         {
559                 float lacol[3] = {0.f, 0.f, 0.f};
560                 lar= go->lampren;
561                 
562                 if (lar) {
563                         vol_shade_one_lamp(shi, co, view, lar, lacol);
564                         add_v3_v3(scatter_col, lacol);
565                 }
566         }
567 }
568
569         
570 /*
571 The main volumetric integrator, using an emission/absorption/scattering model.
572
573 Incoming radiance = 
574
575 outgoing radiance from behind surface * beam transmittance/attenuation
576 + added radiance from all points along the ray due to participating media
577         --> radiance for each segment = 
578                 (radiance added by scattering + radiance added by emission) * beam transmittance/attenuation
579 */
580
581 /* For ease of use, I've also introduced a 'reflection' and 'reflection color' parameter, which isn't 
582  * physically correct. This works as an RGB tint/gain on out-scattered light, but doesn't affect the light 
583  * that is transmitted through the volume. While having wavelength dependent absorption/scattering is more correct,
584  * it also makes it harder to control the overall look of the volume since coloring the outscattered light results
585  * in the inverse color being transmitted through the rest of the volume.
586  */
587 static void volumeintegrate(struct ShadeInput *shi, float col[4], const float co[3], const float endco[3])
588 {
589         float radiance[3] = {0.f, 0.f, 0.f};
590         float tr[3] = {1.f, 1.f, 1.f};
591         float p[3] = {co[0], co[1], co[2]};
592         float step_vec[3] = {endco[0] - co[0], endco[1] - co[1], endco[2] - co[2]};
593         const float stepsize = shi->mat->vol.stepsize;
594         
595         float t0 = 0.f;
596         float pt0 = t0;
597         float t1 = normalize_v3(step_vec);      /* returns vector length */
598         
599         t0 += stepsize * ((shi->mat->vol.stepsize_type == MA_VOL_STEP_CONSTANT) ? 0.5f : BLI_thread_frand(shi->thread));
600         p[0] += t0 * step_vec[0];
601         p[1] += t0 * step_vec[1];
602         p[2] += t0 * step_vec[2];
603         mul_v3_fl(step_vec, stepsize);
604         
605         for (; t0 < t1; pt0 = t0, t0 += stepsize) {
606                 const float density = vol_get_density(shi, p);
607                 
608                 if (density > 0.00001f) {
609                         float scatter_col[3] = {0.f, 0.f, 0.f}, emit_col[3];
610                         const float stepd = (t0 - pt0) * density;
611                         
612                         /* transmittance component (alpha) */
613                         vol_get_transmittance_seg(shi, tr, stepsize, co, density);
614                         
615                         if (t0 > t1 * 0.25f) {
616                                 /* only use depth cutoff after we've traced a little way into the volume */
617                                 if (luminance(tr) < shi->mat->vol.depth_cutoff) break;
618                         }
619                         
620                         vol_get_emission(shi, emit_col, p);
621                         
622                         if (shi->obi->volume_precache) {
623                                 float p2[3];
624                                 
625                                 p2[0] = p[0] + (step_vec[0] * 0.5f);
626                                 p2[1] = p[1] + (step_vec[1] * 0.5f);
627                                 p2[2] = p[2] + (step_vec[2] * 0.5f);
628                                 
629                                 vol_get_precached_scattering(&R, shi, scatter_col, p2);
630                         } else
631                                 vol_get_scattering(shi, scatter_col, p, shi->view);
632                         
633                         radiance[0] += stepd * tr[0] * (emit_col[0] + scatter_col[0]);
634                         radiance[1] += stepd * tr[1] * (emit_col[1] + scatter_col[1]);
635                         radiance[2] += stepd * tr[2] * (emit_col[2] + scatter_col[2]);
636                 }
637                 add_v3_v3(p, step_vec);
638         }
639         
640         /* multiply original color (from behind volume) with transmittance over entire distance */
641         mul_v3_v3v3(col, tr, col);
642         add_v3_v3(col, radiance);
643         
644         /* alpha <-- transmission luminance */
645         col[3] = 1.0f - luminance(tr);
646 }
647
648 /* the main entry point for volume shading */
649 static void volume_trace(struct ShadeInput *shi, struct ShadeResult *shr, int inside_volume)
650 {
651         float hitco[3], col[4] = {0.f,0.f,0.f,0.f};
652         float *startco, *endco;
653         int trace_behind = 1;
654         const int ztransp= ((shi->depth==0) && (shi->mat->mode & MA_TRANSP) && (shi->mat->mode & MA_ZTRANSP));
655         Isect is;
656
657         /* check for shading an internal face a volume object directly */
658         if (inside_volume == VOL_SHADE_INSIDE)
659                 trace_behind = 0;
660         else if (inside_volume == VOL_SHADE_OUTSIDE) {
661                 if (shi->flippednor)
662                         inside_volume = VOL_SHADE_INSIDE;
663         }
664         
665         if (ztransp && inside_volume == VOL_SHADE_INSIDE) {
666                 MatInside *mi;
667                 int render_this=0;
668                 
669                 /* don't render the backfaces of ztransp volume materials.
670                  
671                  * volume shading renders the internal volume from between the
672                  * ' view intersection of the solid volume to the
673                  * intersection on the other side, as part of the shading of
674                  * the front face.
675                  
676                  * Because ztransp renders both front and back faces independently
677                  * this will double up, so here we prevent rendering the backface as well, 
678                  * which would otherwise render the volume in between the camera and the backface
679                  * --matt */
680                 
681                 for (mi=R.render_volumes_inside.first; mi; mi=mi->next) {
682                         /* weak... */
683                         if (mi->ma == shi->mat) render_this=1;
684                 }
685                 if (!render_this) return;
686         }
687         
688
689         if (inside_volume == VOL_SHADE_INSIDE)
690         {
691                 startco = shi->camera_co;
692                 endco = shi->co;
693                 
694                 if (trace_behind) {
695                         if (!ztransp)
696                                 /* trace behind the volume object */
697                                 vol_trace_behind(shi, shi->vlr, endco, col);
698                 } else {
699                         /* we're tracing through the volume between the camera 
700                          * and a solid surface, so use that pre-shaded radiance */
701                         copy_v4_v4(col, shr->combined);
702                 }
703                 
704                 /* shade volume from 'camera' to 1st hit point */
705                 volumeintegrate(shi, col, startco, endco);
706         }
707         /* trace to find a backface, the other side bounds of the volume */
708         /* (ray intersect ignores front faces here) */
709         else if (vol_get_bounds(shi, shi->co, shi->view, hitco, &is, VOL_BOUNDS_DEPTH))
710         {
711                 VlakRen *vlr = (VlakRen *)is.hit.face;
712                 
713                 startco = shi->co;
714                 endco = hitco;
715                 
716                 if (!ztransp) {
717                         /* if it's another face in the same material */
718                         if (vlr->mat == shi->mat) {
719                                 /* trace behind the 2nd (raytrace) hit point */
720                                 vol_trace_behind(shi, (VlakRen *)is.hit.face, endco, col);
721                         } else {
722                                 shade_intersection(shi, col, &is);
723                         }
724                 }
725                 
726                 /* shade volume from 1st hit point to 2nd hit point */
727                 volumeintegrate(shi, col, startco, endco);
728         }
729         
730         if (ztransp)
731                 col[3] = col[3]>1.f?1.f:col[3];
732         else
733                 col[3] = 1.f;
734         
735         copy_v3_v3(shr->combined, col);
736         shr->alpha = col[3];
737         
738         copy_v3_v3(shr->diff, shr->combined);
739 }
740
741 /* Traces a shadow through the object, 
742  * pretty much gets the transmission over a ray path */
743 void shade_volume_shadow(struct ShadeInput *shi, struct ShadeResult *shr, struct Isect *last_is)
744 {
745         float hitco[3];
746         float tr[3] = {1.0,1.0,1.0};
747         Isect is= {{0}};
748         float *startco, *endco;
749         int intersect_type = VOL_BOUNDS_DEPTH;
750
751         memset(shr, 0, sizeof(ShadeResult));
752         
753         /* if 1st hit normal is facing away from the camera, 
754          * then we're inside the volume already. */
755         if (shi->flippednor) {
756                 startco = last_is->start;
757                 endco = shi->co;
758                 intersect_type = VOL_BOUNDS_SS;
759         }
760         
761         /* trace to find a backface, the other side bounds of the volume */
762         /* (ray intersect ignores front faces here) */
763         else if (vol_get_bounds(shi, shi->co, shi->view, hitco, &is, intersect_type)) {
764                 startco = shi->co;
765                 endco = hitco;
766         }
767         else {
768                 shr->combined[0] = shr->combined[1] = shr->combined[2] = 0.f;
769                 shr->alpha = shr->combined[3] = 1.f;
770                 return;
771         }
772
773         vol_get_transmittance(shi, tr, startco, endco);
774
775         
776         /* if we hit another face in the same volume bounds */
777         /* shift raytrace coordinates to the hit point, to avoid shading volume twice */
778         /* due to idiosyncracy in ray_trace_shadow_tra() */
779         if (is.hit.ob == shi->obi) {
780                 copy_v3_v3(shi->co, hitco);
781                 last_is->dist -= is.dist;
782                 shi->vlr = (VlakRen *)is.hit.face;
783         }
784
785         
786         copy_v3_v3(shr->combined, tr);
787         shr->combined[3] = 1.0f - luminance(tr);
788         shr->alpha = shr->combined[3];
789 }
790
791
792 /* delivers a fully filled in ShadeResult, for all passes */
793 void shade_volume_outside(ShadeInput *shi, ShadeResult *shr)
794 {
795         memset(shr, 0, sizeof(ShadeResult));
796         volume_trace(shi, shr, VOL_SHADE_OUTSIDE);
797 }
798
799
800 void shade_volume_inside(ShadeInput *shi, ShadeResult *shr)
801 {
802         MatInside *m;
803         Material *mat_backup;
804         ObjectInstanceRen *obi_backup;
805         float prev_alpha = shr->alpha;
806
807         /* XXX: extend to multiple volumes perhaps later */
808         mat_backup = shi->mat;
809         obi_backup = shi->obi;
810         
811         m = R.render_volumes_inside.first;
812         shi->mat = m->ma;
813         shi->obi = m->obi;
814         shi->obr = m->obi->obr;
815         
816         volume_trace(shi, shr, VOL_SHADE_INSIDE);
817         
818         shr->alpha = shr->alpha + prev_alpha;
819         CLAMP(shr->alpha, 0.0f, 1.0f);
820
821         shi->mat = mat_backup;
822         shi->obi = obi_backup;
823         shi->obr = obi_backup->obr;
824 }
825
826