* Volume Rendering: Voxel data
[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., 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
40 #include "RE_shader_ext.h"
41 #include "RE_raytrace.h"
42
43 #include "DNA_material_types.h"
44 #include "DNA_group_types.h"
45 #include "DNA_lamp_types.h"
46
47 #include "BKE_global.h"
48
49 #include "render_types.h"
50 #include "pixelshading.h"
51 #include "shading.h"
52 #include "texture.h"
53 #include "volumetric.h"
54
55 #if defined( _MSC_VER ) && !defined( __cplusplus )
56 # define inline __inline
57 #endif // defined( _MSC_VER ) && !defined( __cplusplus )
58
59 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
60 /* defined in pipeline.c, is hardcopy of active dynamic allocated Render */
61 /* only to be used here in this file, it's for speed */
62 extern struct Render R;
63 /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
64
65 static int vol_backface_intersect_check(Isect *is, int ob, RayFace *face)
66 {
67         VlakRen *vlr = (VlakRen *)face;
68         
69         /* only consider faces away, so overlapping layers
70          * of foward facing geometry don't cause the ray to stop */
71         return (INPR(is->vec, vlr->n) < 0.0f);
72 }
73
74 /* TODO: Box or sphere intersection types could speed things up */
75 static int vol_get_bounds(ShadeInput *shi, float *co, float *vec, float *hitco, Isect *isect, int intersect_type, int checkfunc)
76 {
77         float maxsize = RE_ray_tree_max_size(R.raytree);
78         int intersected=0;
79
80         /* TODO: use object's bounding box to calculate max size */
81         VECCOPY(isect->start, co);
82         isect->end[0] = co[0] + vec[0] * maxsize;
83         isect->end[1] = co[1] + vec[1] * maxsize;
84         isect->end[2] = co[2] + vec[2] * maxsize;
85         
86         isect->mode= RE_RAY_MIRROR;
87         isect->oborig= RAY_OBJECT_SET(&R, shi->obi);
88         isect->face_last= NULL;
89         isect->ob_last= 0;
90         isect->lay= -1;
91         
92         if (intersect_type == VOL_BOUNDS_DEPTH) isect->faceorig= (RayFace*)shi->vlr;
93         else if (intersect_type == VOL_BOUNDS_SS) isect->faceorig= NULL;
94         
95         if (checkfunc==VOL_IS_BACKFACE)
96                 intersected = RE_ray_tree_intersect_check(R.raytree, isect, vol_backface_intersect_check);
97         else
98                 intersected = RE_ray_tree_intersect(R.raytree, isect);
99         
100         if(intersected)
101         {
102                 float isvec[3];
103
104                 VECCOPY(isvec, isect->vec);
105                 hitco[0] = isect->start[0] + isect->labda*isvec[0];
106                 hitco[1] = isect->start[1] + isect->labda*isvec[1];
107                 hitco[2] = isect->start[2] + isect->labda*isvec[2];
108                 
109                 return 1;
110         } else {
111                 return 0;
112         }
113 }
114
115 float vol_get_stepsize(struct ShadeInput *shi, int context)
116 {
117         if (shi->mat->vol_stepsize_type == MA_VOL_STEP_RANDOMIZED) {
118                 /* range between 0.75 and 1.25 */
119                 const float rnd = 0.5f * BLI_thread_frand(shi->thread) + 0.75f;
120         
121                 if (context == STEPSIZE_VIEW)
122                         return shi->mat->vol_stepsize * rnd;
123                 else if (context == STEPSIZE_SHADE)
124                         return shi->mat->vol_shade_stepsize * rnd;
125         }
126         else {  // MA_VOL_STEP_CONSTANT
127                 
128                 if (context == STEPSIZE_VIEW)
129                         return shi->mat->vol_stepsize;
130                 else if (context == STEPSIZE_SHADE)
131                         return shi->mat->vol_shade_stepsize;
132         }
133         
134         return shi->mat->vol_stepsize;
135 }
136
137 static float vol_get_depth_cutoff(struct ShadeInput *shi)
138 {
139         return shi->mat->vol_depth_cutoff;
140 }
141
142 /* SHADING */
143
144 static float D(ShadeInput *shi, int rgb, int x, int y, int z)
145 {
146         const int res = shi->mat->vol_precache_resolution;
147         CLAMP(x, 0, res-1);
148         CLAMP(y, 0, res-1);
149         CLAMP(z, 0, res-1);
150         return shi->obi->volume_precache[rgb*res*res*res + x*res*res + y*res + z];
151 }
152
153 inline float lerp(float t, float v1, float v2) {
154         return (1.f - t) * v1 + t * v2;
155 }
156
157 /* trilinear interpolation */
158 static void vol_get_precached_scattering(ShadeInput *shi, float *scatter_col, float *co)
159 {
160         const int res = shi->mat->vol_precache_resolution;
161         float voxx, voxy, voxz;
162         int vx, vy, vz;
163         float dx, dy, dz;
164         float d00, d10, d01, d11, d0, d1, d_final;
165         float bbmin[3], bbmax[3], dim[3];
166         int rgb;
167         
168         if (!shi->obi->volume_precache) return;
169         
170         VECCOPY(bbmin, shi->obi->obr->boundbox[0]);
171         VECCOPY(bbmax, shi->obi->obr->boundbox[1]);
172         VecSubf(dim, bbmax, bbmin);
173         
174         voxx = ((co[0] - bbmin[0]) / dim[0]) * res - 0.5f;
175         voxy = ((co[1] - bbmin[1]) / dim[1]) * res - 0.5f;
176         voxz = ((co[2] - bbmin[2]) / dim[2]) * res - 0.5f;
177         
178         vx = (int)voxx; vy = (int)voxy; vz = (int)voxz;
179         
180         dx = voxx - vx; dy = voxy - vy; dz = voxz - vz;
181         
182         for (rgb=0; rgb < 3; rgb++) {
183                 d00 = lerp(dx, D(shi, rgb, vx, vy, vz),                 D(shi, rgb, vx+1, vy, vz));
184                 d10 = lerp(dx, D(shi, rgb, vx, vy+1, vz),               D(shi, rgb, vx+1, vy+1, vz));
185                 d01 = lerp(dx, D(shi, rgb, vx, vy, vz+1),               D(shi, rgb, vx+1, vy, vz+1));
186                 d11 = lerp(dx, D(shi, rgb, vx, vy+1, vz+1),     D(shi, rgb, vx+1, vy+1, vz+1));
187                 d0 = lerp(dy, d00, d10);
188                 d1 = lerp(dy, d01, d11);
189                 d_final = lerp(dz, d0, d1);
190                 
191                 scatter_col[rgb] = d_final;
192         }
193 }
194
195 /* no interpolation */
196 static void vol_get_precached_scattering_nearest(ShadeInput *shi, float *scatter_col, float *co)
197 {
198         const int res = shi->mat->vol_precache_resolution;
199         int x,y,z;
200         float bbmin[3], bbmax[3], dim[3];
201
202         if (!shi->obi->volume_precache) return;
203         
204         VECCOPY(bbmin, shi->obi->obr->boundbox[0]);
205         VECCOPY(bbmax, shi->obi->obr->boundbox[1]);
206         VecSubf(dim, bbmax, bbmin);
207         
208         x = (int)(((co[0] - bbmin[0]) / dim[0]) * res);
209         y = (int)(((co[1] - bbmin[1]) / dim[1]) * res);
210         z = (int)(((co[2] - bbmin[2]) / dim[2]) * res);
211         
212         scatter_col[0] = shi->obi->volume_precache[0*res*res*res + x*res*res + y*res + z];
213         scatter_col[1] = shi->obi->volume_precache[1*res*res*res + x*res*res + y*res + z];
214         scatter_col[2] = shi->obi->volume_precache[2*res*res*res + x*res*res + y*res + z];
215 }
216
217 float vol_get_density(struct ShadeInput *shi, float *co)
218 {
219         float density = shi->mat->alpha;
220         float density_scale = shi->mat->vol_density_scale;
221         float col[3] = {0.0, 0.0, 0.0};
222         
223         if (shi->mat->flag & MA_IS_TEXTURED) {
224                 do_volume_tex(shi, co, MAP_ALPHA, col, &density);
225         }
226         
227         return density * density_scale;
228 }
229
230 /* compute emission component, amount of radiance to add per segment
231  * can be textured with 'emit' */
232 void vol_get_emission(ShadeInput *shi, float *em, float *co, float density)
233 {
234         float emission = shi->mat->emit;
235         float col[3] = {0.0, 0.0, 0.0};
236         
237         VECCOPY(col, &shi->mat->r);
238         
239         do_volume_tex(shi, co, MAP_EMIT+MAP_COL, col, &emission);
240         
241         em[0] = em[1] = em[2] = emission * density;
242         VecMulVecf(em, em, col);
243 }
244
245 /* scattering multiplier, values above 1.0 are non-physical, 
246  * but can be useful to tweak lighting */
247 void vol_get_scattering_fac(ShadeInput *shi, float *scatter_fac, float *co, float density)
248 {
249         *scatter_fac = shi->mat->vol_scattering;
250 }
251
252 /* phase function - determines in which directions the light 
253  * is scattered in the volume relative to incoming direction 
254  * and view direction */
255 float vol_get_phasefunc(ShadeInput *shi, short phasefunc_type, float g, float *w, float *wp)
256 {
257         const float costheta = Inpf(w, wp);
258         
259         if (phasefunc_type == MA_VOL_PH_ISOTROPIC) {
260                 return 1.f / (4.f * M_PI);
261         }
262         else if (phasefunc_type == MA_VOL_PH_MIEHAZY) {
263                 return (0.5f + 4.5f * powf(0.5 * (1.f + costheta), 8.f)) / (4.f*M_PI);
264         }
265         else if (phasefunc_type == MA_VOL_PH_MIEMURKY) {
266                 return (0.5f + 16.5f * powf(0.5 * (1.f + costheta), 32.f)) / (4.f*M_PI);
267         }
268         else if (phasefunc_type == MA_VOL_PH_RAYLEIGH) {
269                 return 3.f/(16.f*M_PI) * (1 + costheta * costheta);
270         }
271         else if (phasefunc_type == MA_VOL_PH_HG) {
272                 return 1.f / (4.f * M_PI) * (1.f - g*g) / powf(1.f + g*g - 2.f * g * costheta, 1.5f);
273         }
274         else if (phasefunc_type == MA_VOL_PH_SCHLICK) {
275                 const float k = 1.55f * g - .55f * g * g * g;
276                 const float kcostheta = k * costheta;
277                 return 1.f / (4.f * M_PI) * (1.f - k*k) / ((1.f - kcostheta) * (1.f - kcostheta));
278         } else {
279                 return 1.0f;
280         }
281 }
282
283 void vol_get_absorption(ShadeInput *shi, float *absorb_col, float *co)
284 {
285         float dummy = 1.0f;
286         const float absorption = shi->mat->vol_absorption;
287         
288         VECCOPY(absorb_col, shi->mat->vol_absorption_col);
289         
290         if (shi->mat->flag & MA_IS_TEXTURED)
291                 do_volume_tex(shi, co, MAP_COLMIR, absorb_col, &dummy);
292         
293         absorb_col[0] = (1.0f - absorb_col[0]) * absorption;
294         absorb_col[1] = (1.0f - absorb_col[1]) * absorption;
295         absorb_col[2] = (1.0f - absorb_col[2]) * absorption;
296 }
297
298 /* Compute attenuation, otherwise known as 'optical thickness', extinction, or tau.
299  * Used in the relationship Transmittance = e^(-attenuation)
300  */
301 void vol_get_attenuation(ShadeInput *shi, float *tau, float *co, float *endco, float density, float stepsize)
302 {
303         /* input density = density at co */
304         float absorb_col[3];
305         int s, nsteps;
306         float step_vec[3], step_sta[3], step_end[3];
307         const float dist = VecLenf(co, endco);
308
309         vol_get_absorption(shi, absorb_col, co);
310
311         nsteps = (int)((dist / stepsize) + 0.5);
312         
313         /* trigger for recalculating density */
314         if (density < -0.001f) density = vol_get_density(shi, co);
315         
316         if (nsteps == 1) {
317                 /* homogenous volume within the sampled distance */
318                 tau[0] = tau[1] = tau[2] = dist * density;
319                 
320                 VecMulVecf(tau, tau, absorb_col);
321                 return;
322         } else {
323                 tau[0] = tau[1] = tau[2] = 0.0;
324         }
325         
326         VecSubf(step_vec, endco, co);
327         VecMulf(step_vec, 1.0f / nsteps);
328         
329         VecCopyf(step_sta, co);
330         VecAddf(step_end, step_sta, step_vec);
331         
332         for (s = 0;  s < nsteps; s++) {
333                 if (s > 0)
334                         density = vol_get_density(shi, step_sta);
335                 
336                 tau[0] += stepsize * density;
337                 tau[1] += stepsize * density;
338                 tau[2] += stepsize * density;
339                 
340                 if (s < nsteps-1) {
341                         VECCOPY(step_sta, step_end);
342                         VecAddf(step_end, step_end, step_vec);
343                 }
344         }
345         VecMulVecf(tau, tau, absorb_col);
346 }
347
348 void vol_shade_one_lamp(struct ShadeInput *shi, float *co, LampRen *lar, float *lacol, float stepsize, float density)
349 {
350         float visifac, lv[3], lampdist;
351         float tau[3], tr[3]={1.0,1.0,1.0};
352         float hitco[3], *atten_co;
353         float p;
354         float scatter_fac;
355         float shade_stepsize = vol_get_stepsize(shi, STEPSIZE_SHADE);
356         
357         if (lar->mode & LA_LAYER) if((lar->lay & shi->obi->lay)==0) return;
358         if ((lar->lay & shi->lay)==0) return;
359         if (lar->energy == 0.0) return;
360         
361         visifac= lamp_get_visibility(lar, co, lv, &lampdist);
362         if(visifac==0.0f) return;
363
364         lacol[0] = lar->r;
365         lacol[1] = lar->g;
366         lacol[2] = lar->b;
367         
368         if(lar->mode & LA_TEXTURE) {
369                 shi->osatex= 0;
370                 do_lamp_tex(lar, lv, shi, lacol, LA_TEXTURE);
371         }
372
373         VecMulf(lacol, visifac*lar->energy);
374
375         if (ELEM(lar->type, LA_SUN, LA_HEMI))
376                 VECCOPY(lv, lar->vec);
377         VecMulf(lv, -1.0f);
378         
379         p = vol_get_phasefunc(shi, shi->mat->vol_phasefunc_type, shi->mat->vol_phasefunc_g, shi->view, lv);
380         VecMulf(lacol, p);
381         
382         if (shi->mat->vol_shadeflag & MA_VOL_ATTENUATED) {
383                 Isect is;
384                 
385                 /* find minimum of volume bounds, or lamp coord */
386                 if (vol_get_bounds(shi, co, lv, hitco, &is, VOL_BOUNDS_SS, 0)) {
387                         float dist = VecLenf(co, hitco);
388                         VlakRen *vlr = (VlakRen *)is.face;
389                         
390                         /* simple internal shadowing */
391                         if (vlr->mat->material_type == MA_SOLID) {
392                                 lacol[0] = lacol[1] = lacol[2] = 0.0f;
393                                 return;
394                         }
395
396                         if (ELEM(lar->type, LA_SUN, LA_HEMI))
397                                 atten_co = hitco;
398                         else if ( lampdist < dist ) {
399                                 atten_co = lar->co;
400                         } else
401                                 atten_co = hitco;
402                         
403                         vol_get_attenuation(shi, tau, co, atten_co, density, shade_stepsize);
404                         tr[0] = exp(-tau[0]);
405                         tr[1] = exp(-tau[1]);
406                         tr[2] = exp(-tau[2]);
407                         
408                         VecMulVecf(lacol, lacol, tr);
409                 }
410                 else {
411                         /* Point is on the outside edge of the volume,
412                          * therefore no attenuation, full transmission.
413                          * Radiance from lamp remains unchanged */
414                 }
415         }
416         
417         vol_get_scattering_fac(shi, &scatter_fac, co, density);
418         VecMulf(lacol, scatter_fac);
419 }
420
421 /* single scattering only for now */
422 void vol_get_scattering(ShadeInput *shi, float *scatter, float *co, float stepsize, float density)
423 {
424         ListBase *lights;
425         GroupObject *go;
426         LampRen *lar;
427         float col[3] = {0.f, 0.f, 0.f};
428         int i=0;
429
430         lights= get_lights(shi);
431         for(go=lights->first; go; go= go->next)
432         {
433                 float lacol[3] = {0.f, 0.f, 0.f};
434         
435                 i++;
436         
437                 lar= go->lampren;
438                 if (lar) {
439                         vol_shade_one_lamp(shi, co, lar, lacol, stepsize, density);
440                         VecAddf(col, col, lacol);
441                 }
442         }
443         
444         VECCOPY(scatter, col);
445 }
446
447         
448 /*
449 The main volumetric integrator, using an emission/absorption/scattering model.
450
451 Incoming radiance = 
452
453 outgoing radiance from behind surface * beam transmittance/attenuation
454 + added radiance from all points along the ray due to participating media
455         --> radiance for each segment = 
456                 (radiance added by scattering + radiance added by emission) * beam transmittance/attenuation
457
458 -- To find transmittance:
459         compute optical thickness with tau (perhaps involving monte carlo integration)
460         transmittance = exp(-tau)
461         
462 -- To find radiance from segments along the way:
463         find radiance for one step: 
464         - loop over lights and weight by phase function
465 */
466 static void volumeintegrate(struct ShadeInput *shi, float *col, float *co, float *endco)
467 {
468         float tr[3] = {1.0f, 1.0f, 1.0f};
469         float radiance[3] = {0.f, 0.f, 0.f}, d_radiance[3] = {0.f, 0.f, 0.f};
470         float stepsize = vol_get_stepsize(shi, STEPSIZE_VIEW);
471         int nsteps, s;
472         float tau[3], emit_col[3], scatter_col[3] = {0.0, 0.0, 0.0};
473         float stepvec[3], step_sta[3], step_end[3], step_mid[3];
474         float density = vol_get_density(shi, co);
475         const float depth_cutoff = vol_get_depth_cutoff(shi);
476         
477         /* multiply col_behind with beam transmittance over entire distance */
478         vol_get_attenuation(shi, tau, co, endco, density, stepsize);
479         tr[0] *= exp(-tau[0]);
480         tr[1] *= exp(-tau[1]);
481         tr[2] *= exp(-tau[2]);
482         VecMulVecf(radiance, tr, col);  
483         tr[0] = tr[1] = tr[2] = 1.0f;
484         
485         /* ray marching */
486         nsteps = (int)((VecLenf(co, endco) / stepsize) + 0.5);
487         
488         VecSubf(stepvec, endco, co);
489         VecMulf(stepvec, 1.0f / nsteps);
490         VecCopyf(step_sta, co);
491         VecAddf(step_end, step_sta, stepvec);
492         
493         /* get radiance from all points along the ray due to participating media */
494         for (s = 0; s < nsteps; s++) {
495
496                 if (s > 0) density = vol_get_density(shi, step_sta);
497                 
498                 /* there's only any use in shading here if there's actually some density to shade! */
499                 if (density > 0.01f) {
500                 
501                         /* transmittance component (alpha) */
502                         vol_get_attenuation(shi, tau, step_sta, step_end, density, stepsize);
503                         tr[0] *= exp(-tau[0]);
504                         tr[1] *= exp(-tau[1]);
505                         tr[2] *= exp(-tau[2]);
506                         
507                         step_mid[0] = step_sta[0] + (stepvec[0] * 0.5);
508                         step_mid[1] = step_sta[1] + (stepvec[1] * 0.5);
509                         step_mid[2] = step_sta[2] + (stepvec[2] * 0.5);
510                 
511                         /* incoming light via emission or scattering (additive) */
512                         vol_get_emission(shi, emit_col, step_mid, density);
513                         
514                         if ((shi->mat->vol_shadeflag & MA_VOL_PRECACHESHADING) &&
515                                 (shi->mat->vol_shadeflag & MA_VOL_ATTENUATED)) {
516                                 vol_get_precached_scattering(shi, scatter_col, step_mid);
517                         } else
518                                 vol_get_scattering(shi, scatter_col, step_mid, stepsize, density);
519                                                 
520                         VecMulf(scatter_col, density);
521                         VecAddf(d_radiance, emit_col, scatter_col);
522                         
523                         /*   Lv += Tr * (Lve() + Ld) */
524                         VecMulVecf(d_radiance, tr, d_radiance);
525                         VecMulf(d_radiance, stepsize);
526                         
527                         VecAddf(radiance, radiance, d_radiance);        
528                 }
529
530                 VecCopyf(step_sta, step_end);
531                 VecAddf(step_end, step_end, stepvec);
532                 
533                 /* luminance rec. 709 */
534                 if ((0.2126*tr[0] + 0.7152*tr[1] + 0.0722*tr[2]) < depth_cutoff) break; 
535         }
536         
537         VecCopyf(col, radiance);
538         col[3] = 1.0f -(tr[0] + tr[1] + tr[2]) * 0.333f;
539 }
540
541 static void shade_intersection(ShadeInput *shi, float *col, Isect *is)
542 {
543         ShadeInput shi_new;
544         ShadeResult shr_new;
545         
546         memset(&shi_new, 0, sizeof(ShadeInput)); 
547         
548         shi_new.mask= shi->mask;
549         shi_new.osatex= shi->osatex;
550         shi_new.thread= shi->thread;
551         shi_new.depth= 1;
552         shi_new.volume_depth= shi->volume_depth + 1;
553         shi_new.xs= shi->xs;
554         shi_new.ys= shi->ys;
555         shi_new.lay= shi->lay;
556         shi_new.passflag= SCE_PASS_COMBINED; /* result of tracing needs no pass info */
557         shi_new.combinedflag= 0xFFFFFF;          /* ray trace does all options */
558         shi_new.light_override= shi->light_override;
559         shi_new.mat_override= shi->mat_override;
560         
561         VECCOPY(shi_new.camera_co, is->start);
562         
563         memset(&shr_new, 0, sizeof(ShadeResult));
564
565         /* hardcoded limit of 100 for now - prevents problems in weird geometry */
566         if (shi->volume_depth < 100) {
567                 shade_ray(is, &shi_new, &shr_new);
568         }
569         
570         col[0] = shr_new.combined[0];
571         col[1] = shr_new.combined[1];
572         col[2] = shr_new.combined[2];
573         col[3] = shr_new.alpha;
574 }
575
576 static void vol_trace_behind(ShadeInput *shi, VlakRen *vlr, float *co, float *col)
577 {
578         Isect isect;
579         float maxsize = RE_ray_tree_max_size(R.raytree);
580
581         VECCOPY(isect.start, co);
582         isect.end[0] = isect.start[0] + shi->view[0] * maxsize;
583         isect.end[1] = isect.start[1] + shi->view[1] * maxsize;
584         isect.end[2] = isect.start[2] + shi->view[2] * maxsize;
585
586         isect.faceorig= (RayFace *)vlr;
587         
588         isect.mode= RE_RAY_MIRROR;
589         isect.oborig= RAY_OBJECT_SET(&R, shi->obi);
590         isect.face_last= NULL;
591         isect.ob_last= 0;
592         isect.lay= -1;
593         
594         /* check to see if there's anything behind the volume, otherwise shade the sky */
595         if(RE_ray_tree_intersect(R.raytree, &isect)) {
596                 shade_intersection(shi, col, &isect);
597         } else {
598                 shadeSkyView(col, co, shi->view, NULL, shi->thread);
599                 shadeSunView(col, shi->view);
600         }
601 }
602
603 /* the main entry point for volume shading */
604 void volume_trace(struct ShadeInput *shi, struct ShadeResult *shr)
605 {
606         float hitco[3], col[4] = {0.f,0.f,0.f,0.f};
607         Isect is;
608
609         memset(shr, 0, sizeof(ShadeResult));
610
611         /* if 1st hit normal is facing away from the camera, 
612          * then we're inside the volume already. */
613         if (shi->flippednor) {
614                 /* trace behind the 1st hit point */
615                 vol_trace_behind(shi, shi->vlr, shi->co, col);
616                 
617                 /* shade volume from 'camera' to 1st hit point */
618                 volumeintegrate(shi, col, shi->camera_co, shi->co);
619                 
620                 shr->combined[0] = col[0];
621                 shr->combined[1] = col[1];
622                 shr->combined[2] = col[2];
623                 
624                 if (shi->mat->vol_shadeflag & MA_VOL_USEALPHA) {
625                         if (col[3] > 1.0f)
626                                 col[3] = 1.0f;
627                 }
628                 else
629                         col[3] = 1.0f;
630                 shr->combined[3] = col[3];
631                 shr->alpha = col[3];
632                 
633                 VECCOPY(shr->diff, shr->combined);
634         }
635         /* trace to find a backface, the other side bounds of the volume */
636         /* (ray intersect ignores front faces here) */
637         else if (vol_get_bounds(shi, shi->co, shi->view, hitco, &is, VOL_BOUNDS_DEPTH, 0)) {
638                 VlakRen *vlr = (VlakRen *)is.face;
639                 
640                 /* if it's another face in the same material */
641                 if (vlr->mat == shi->mat) {
642                         /* trace behind the 2nd (raytrace) hit point */
643                         vol_trace_behind(shi, (VlakRen *)is.face, hitco, col);
644                 } else {
645                         shade_intersection(shi, col, &is);
646                 }
647         
648                 /* shade volume from 1st hit point to 2nd hit point */
649                 volumeintegrate(shi, col, shi->co, hitco);
650                 
651                 shr->combined[0] = col[0];
652                 shr->combined[1] = col[1];
653                 shr->combined[2] = col[2];
654                 
655                 if (shi->mat->vol_shadeflag & MA_VOL_USEALPHA) {
656                         if (col[3] > 1.0f)
657                                 col[3] = 1.0f;
658                 }
659                 else
660                         col[3] = 1.0f;
661                 shr->combined[3] = col[3];
662                 shr->alpha = col[3];
663                 
664                 VECCOPY(shr->diff, shr->combined);
665         }
666         else {
667                 shr->combined[0] = 0.0f;
668                 shr->combined[1] = 0.0f;
669                 shr->combined[2] = 0.0f;
670                 shr->combined[3] = shr->alpha =  1.0f;
671         }
672 }
673
674 /* Traces a shadow through the object, 
675  * pretty much gets the transmission over a ray path */
676 void volume_trace_shadow(struct ShadeInput *shi, struct ShadeResult *shr, struct Isect *last_is)
677 {
678         float hitco[3];
679         float tr[3] = {1.0,1.0,1.0};
680         float tau[3] = {0.0,0.0,0.0};
681         Isect is;
682         float shade_stepsize = vol_get_stepsize(shi, STEPSIZE_SHADE);
683
684         memset(shr, 0, sizeof(ShadeResult));
685         
686         /* if 1st hit normal is facing away from the camera, 
687          * then we're inside the volume already. */
688         if (shi->flippednor) {
689         
690                 vol_get_attenuation(shi, tau, last_is->start, shi->co, -1.0f, shade_stepsize);
691                 tr[0] = exp(-tau[0]);
692                 tr[1] = exp(-tau[1]);
693                 tr[2] = exp(-tau[2]);
694                 
695                 shr->combined[0] = tr[0];
696                 shr->combined[1] = tr[1];
697                 shr->combined[2] = tr[2];
698                 
699                 shr->combined[3] = 1.0f -(tr[0] + tr[1] + tr[2]) * 0.333f;
700                 shr->alpha = shr->combined[3];
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, 0)) {
705                 
706                 vol_get_attenuation(shi, tau, shi->co, hitco, -1.0f, shade_stepsize);
707                 tr[0] = exp(-tau[0]);
708                 tr[1] = exp(-tau[1]);
709                 tr[2] = exp(-tau[2]);
710                 
711                 shr->combined[0] = tr[0];
712                 shr->combined[1] = tr[1];
713                 shr->combined[2] = tr[2];
714                 
715                 shr->combined[3] = 1.0f -(tr[0] + tr[1] + tr[2]) * 0.333f;
716                 shr->alpha = shr->combined[3];
717
718         }
719         else {
720                 shr->combined[0] = 0.0f;
721                 shr->combined[1] = 0.0f;
722                 shr->combined[2] = 0.0f;
723                 shr->combined[3] = shr->alpha =  0.0f;
724         }
725 }
726