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