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