ec0c3114213578593c2e68afdb31cb05cf51afe6
[blender.git] / intern / cycles / kernel / kernel_volume.h
1 /*
2  * Copyright 2011-2013 Blender Foundation
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License
15  */
16
17 CCL_NAMESPACE_BEGIN
18
19 typedef enum VolumeIntegrateResult {
20         VOLUME_PATH_SCATTERED = 0,
21         VOLUME_PATH_ATTENUATED = 1,
22         VOLUME_PATH_MISSED = 2
23 } VolumeIntegrateResult;
24
25 /* Volume shader properties
26  *
27  * extinction coefficient = absorption coefficient + scattering coefficient
28  * sigma_t = sigma_a + sigma_s */
29
30 typedef struct VolumeShaderCoefficients {
31         float3 sigma_a;
32         float3 sigma_s;
33         float3 emission;
34 } VolumeShaderCoefficients;
35
36 /* evaluate shader to get extinction coefficient at P */
37 ccl_device bool volume_shader_extinction_sample(KernelGlobals *kg, ShaderData *sd, PathState *state, float3 P, float3 *extinction)
38 {
39         sd->P = P;
40         shader_eval_volume(kg, sd, state->volume_stack, PATH_RAY_SHADOW, SHADER_CONTEXT_SHADOW);
41
42         if(!(sd->flag & (SD_ABSORPTION|SD_SCATTER)))
43                 return false;
44
45         float3 sigma_t = make_float3(0.0f, 0.0f, 0.0f);
46
47         for(int i = 0; i < sd->num_closure; i++) {
48                 const ShaderClosure *sc = &sd->closure[i];
49
50                 if(CLOSURE_IS_VOLUME(sc->type))
51                         sigma_t += sc->weight;
52         }
53
54         *extinction = sigma_t;
55         return true;
56 }
57
58 /* evaluate shader to get absorption, scattering and emission at P */
59 ccl_device bool volume_shader_sample(KernelGlobals *kg, ShaderData *sd, PathState *state, float3 P, VolumeShaderCoefficients *coeff)
60 {
61         sd->P = P;
62         shader_eval_volume(kg, sd, state->volume_stack, state->flag, SHADER_CONTEXT_VOLUME);
63
64         if(!(sd->flag & (SD_ABSORPTION|SD_SCATTER|SD_EMISSION)))
65                 return false;
66         
67         coeff->sigma_a = make_float3(0.0f, 0.0f, 0.0f);
68         coeff->sigma_s = make_float3(0.0f, 0.0f, 0.0f);
69         coeff->emission = make_float3(0.0f, 0.0f, 0.0f);
70
71         for(int i = 0; i < sd->num_closure; i++) {
72                 const ShaderClosure *sc = &sd->closure[i];
73
74                 if(sc->type == CLOSURE_VOLUME_ABSORPTION_ID)
75                         coeff->sigma_a += sc->weight;
76                 else if(sc->type == CLOSURE_EMISSION_ID)
77                         coeff->emission += sc->weight;
78                 else if(CLOSURE_IS_VOLUME(sc->type))
79                         coeff->sigma_s += sc->weight;
80         }
81
82         /* when at the max number of bounces, treat scattering as absorption */
83         if(sd->flag & SD_SCATTER) {
84                 if(state->volume_bounce >= kernel_data.integrator.max_volume_bounce) {
85                         coeff->sigma_a += coeff->sigma_s;
86                         coeff->sigma_s = make_float3(0.0f, 0.0f, 0.0f);
87                         sd->flag &= ~SD_SCATTER;
88                         sd->flag |= SD_ABSORPTION;
89                 }
90         }
91
92         return true;
93 }
94
95 ccl_device float3 volume_color_attenuation(float3 sigma, float t)
96 {
97         return make_float3(expf(-sigma.x * t), expf(-sigma.y * t), expf(-sigma.z * t));
98 }
99
100 ccl_device bool volume_stack_is_heterogeneous(KernelGlobals *kg, VolumeStack *stack)
101 {
102         for(int i = 0; stack[i].shader != SHADER_NO_ID; i++) {
103                 int shader_flag = kernel_tex_fetch(__shader_flag, (stack[i].shader & SHADER_MASK)*2);
104
105                 if(shader_flag & SD_HETEROGENEOUS_VOLUME)
106                         return true;
107         }
108
109         return false;
110 }
111
112 /* Volume Shadows
113  *
114  * These functions are used to attenuate shadow rays to lights. Both absorption
115  * and scattering will block light, represented by the extinction coefficient. */
116
117 /* homogenous volume: assume shader evaluation at the starts gives
118  * the extinction coefficient for the entire line segment */
119 ccl_device void kernel_volume_shadow_homogeneous(KernelGlobals *kg, PathState *state, Ray *ray, ShaderData *sd, float3 *throughput)
120 {
121         float3 sigma_t;
122
123         if(volume_shader_extinction_sample(kg, sd, state, ray->P, &sigma_t))
124                 *throughput *= volume_color_attenuation(sigma_t, ray->t);
125 }
126
127 /* heterogeneous volume: integrate stepping through the volume until we
128  * reach the end, get absorbed entirely, or run out of iterations */
129 ccl_device void kernel_volume_shadow_heterogeneous(KernelGlobals *kg, PathState *state, Ray *ray, ShaderData *sd, float3 *throughput)
130 {
131         float3 tp = *throughput;
132         const float tp_eps = 1e-10f; /* todo: this is likely not the right value */
133
134         /* prepare for stepping */
135         int max_steps = kernel_data.integrator.volume_max_steps;
136         float step = kernel_data.integrator.volume_step_size;
137         float random_jitter_offset = lcg_step_float(&state->rng_congruential) * step;
138
139         /* compute extinction at the start */
140         float t = 0.0f;
141         float3 P = ray->P;
142         float3 sigma_t;
143
144         if(!volume_shader_extinction_sample(kg, sd, state, P, &sigma_t))
145                 sigma_t = make_float3(0.0f, 0.0f, 0.0f);
146
147         for(int i = 0; i < max_steps; i++) {
148                 /* advance to new position */
149                 float new_t = min(ray->t, t + random_jitter_offset + i * step);
150                 float3 new_P = ray->P + ray->D * new_t;
151                 float3 new_sigma_t;
152
153                 /* compute attenuation over segment */
154                 if(volume_shader_extinction_sample(kg, sd, state, new_P, &new_sigma_t)) {
155                         /* todo: we could avoid computing expf() for each step by summing,
156                          * because exp(a)*exp(b) = exp(a+b), but we still want a quick
157                          * tp_eps check too */
158                         tp *= volume_color_attenuation(0.5f*(sigma_t + new_sigma_t), new_t - t);
159
160                         /* stop if nearly all light blocked */
161                         if(tp.x < tp_eps && tp.y < tp_eps && tp.z < tp_eps)
162                                 break;
163
164                         sigma_t = new_sigma_t;
165                 }
166                 else {
167                         /* skip empty space */
168                         sigma_t = make_float3(0.0f, 0.0f, 0.0f);
169                 }
170
171                 /* stop if at the end of the volume */
172                 t = new_t;
173                 if(t == ray->t)
174                         break;
175         }
176
177         *throughput = tp;
178 }
179
180 /* get the volume attenuation over line segment defined by ray, with the
181  * assumption that there are no surfaces blocking light between the endpoints */
182 ccl_device_noinline void kernel_volume_shadow(KernelGlobals *kg, PathState *state, Ray *ray, float3 *throughput)
183 {
184         ShaderData sd;
185         shader_setup_from_volume(kg, &sd, ray, state->bounce);
186
187         if(volume_stack_is_heterogeneous(kg, state->volume_stack))
188                 kernel_volume_shadow_heterogeneous(kg, state, ray, &sd, throughput);
189         else
190                 kernel_volume_shadow_homogeneous(kg, state, ray, &sd, throughput);
191 }
192
193 /* Volume Path */
194
195 /* homogenous volume: assume shader evaluation at the starts gives
196  * the volume shading coefficient for the entire line segment */
197 ccl_device VolumeIntegrateResult kernel_volume_integrate_homogeneous(KernelGlobals *kg,
198         PathState *state, Ray *ray, ShaderData *sd, PathRadiance *L, float3 *throughput,
199         RNG *rng)
200 {
201         VolumeShaderCoefficients coeff;
202
203         if(!volume_shader_sample(kg, sd, state, ray->P, &coeff))
204                 return VOLUME_PATH_MISSED;
205
206         int closure_flag = sd->flag;
207         float t = ray->t;
208         float3 new_tp;
209         float3 transmittance;
210
211         /* randomly scatter, and if we do t is shortened */
212         if(closure_flag & SD_SCATTER) {
213                 float3 sigma_t = coeff.sigma_a + coeff.sigma_s;
214
215                 /* set up variables for sampling */
216                 float rphase = path_state_rng_1D(kg, rng, state, PRNG_PHASE);
217                 int channel = (int)(rphase*3.0f);
218                 sd->randb_closure = rphase*3.0f - channel;
219
220                 /* pick random color channel, we use the Veach one-sample
221                  * model with balance heuristic for the channels */
222                 float sample_sigma_t;
223
224                 if(channel == 0)
225                         sample_sigma_t = sigma_t.x;
226                 else if(channel == 1)
227                         sample_sigma_t = sigma_t.y;
228                 else
229                         sample_sigma_t = sigma_t.z;
230
231                 /* distance sampling */
232                 if(kernel_data.integrator.volume_homogeneous_sampling == 0 || !kernel_data.integrator.num_all_lights) { 
233                         /* xi is [0, 1[ so log(0) should never happen, division by zero is
234                          * avoided because sample_sigma_t > 0 when SD_SCATTER is set */
235                         float xi = path_state_rng_1D(kg, rng, state, PRNG_SCATTER_DISTANCE);
236                         float sample_t = min(t, -logf(1.0f - xi)/sample_sigma_t);
237
238                         transmittance = volume_color_attenuation(sigma_t, sample_t);
239
240                         if(sample_t < t) {
241                                 float pdf = dot(sigma_t, transmittance);
242                                 new_tp = *throughput * coeff.sigma_s * transmittance * (3.0f / pdf);
243                                 t = sample_t;
244                         }
245                         else {
246                                 float pdf = (transmittance.x + transmittance.y + transmittance.z);
247                                 new_tp = *throughput * transmittance * (3.0f / pdf);
248                         }
249                 }
250                 /* equi-angular sampling */
251                 else {
252                         /* decide if we are going to scatter or not, based on sigma_t. this
253                          * is not ideal, instead we should perhaps split the path here and
254                          * do both, and at least add multiple importance sampling */
255                         float xi = path_state_rng_1D(kg, rng, state, PRNG_SCATTER_DISTANCE);
256                         float sample_transmittance = expf(-sample_sigma_t * t);
257
258                         if(xi < sample_transmittance) {
259                                 /* no scattering */
260                                 float3 transmittance = volume_color_attenuation(sigma_t, t);
261                                 float pdf = (transmittance.x + transmittance.y + transmittance.z);
262                                 new_tp = *throughput * transmittance * (3.0f / pdf);
263                         }
264                         else {
265                                 /* rescale random number so we can reuse it */
266                                 xi = (xi - sample_transmittance)/(1.0f - sample_transmittance);
267
268                                 /* equi-angular scattering somewhere on segment 0..t */
269                                 /* see "Importance Sampling Techniques for Path Tracing in Participating Media" */
270
271                                 /* light RNGs */
272                                 float light_t = path_state_rng_1D(kg, rng, state, PRNG_LIGHT);
273                                 float light_u, light_v;
274                                 path_state_rng_2D(kg, rng, state, PRNG_LIGHT_U, &light_u, &light_v);
275
276                                 /* light sample */
277                                 LightSample ls;
278                                 light_sample(kg, light_t, light_u, light_v, ray->time, ray->P, &ls);
279                                 if(ls.pdf == 0.0f)
280                                         return VOLUME_PATH_MISSED;
281
282                                 /* sampling */
283                                 float delta = dot((ls.P - ray->P) , ray->D);
284                                 float D = sqrtf(len_squared(ls.P - ray->P) - delta * delta);
285                                 float theta_a = -atan2f(delta, D);
286                                 float theta_b = atan2f(t - delta, D);
287                                 float t_ = D * tan((xi * theta_b) + (1 - xi) * theta_a);
288
289                                 float pdf = D / ((theta_b - theta_a) * (D * D + t_ * t_));
290                                 float sample_t = min(t, delta + t_);
291
292                                 transmittance = volume_color_attenuation(sigma_t, sample_t);
293
294                                 new_tp = *throughput * coeff.sigma_s * transmittance / ((1.0f - sample_transmittance) * pdf);
295                                 t = sample_t;
296                         }
297                 }
298         }
299         else if(closure_flag & SD_ABSORPTION) {
300                 /* absorption only, no sampling needed */
301                 transmittance = volume_color_attenuation(coeff.sigma_a, t);
302                 new_tp = *throughput * transmittance;
303         }
304
305         /* integrate emission attenuated by extinction
306          * integral E * exp(-sigma_t * t) from 0 to t = E * (1 - exp(-sigma_t * t))/sigma_t
307          * this goes to E * t as sigma_t goes to zero
308          *
309          * todo: we should use an epsilon to avoid precision issues near zero sigma_t */
310         if(closure_flag & SD_EMISSION) {
311                 float3 emission = coeff.emission;
312
313                 if(closure_flag & SD_ABSORPTION) {
314                         float3 sigma_t = coeff.sigma_a + coeff.sigma_s;
315
316                         emission.x *= (sigma_t.x > 0.0f)? (1.0f - transmittance.x)/sigma_t.x: t;
317                         emission.y *= (sigma_t.y > 0.0f)? (1.0f - transmittance.y)/sigma_t.y: t;
318                         emission.z *= (sigma_t.z > 0.0f)? (1.0f - transmittance.z)/sigma_t.z: t;
319                 }
320                 else
321                         emission *= t;
322
323                 path_radiance_accum_emission(L, *throughput, emission, state->bounce);
324         }
325
326         /* modify throughput */
327         if(closure_flag & (SD_ABSORPTION|SD_SCATTER)) {
328                 *throughput = new_tp;
329
330                 /* prepare to scatter to new direction */
331                 if(t < ray->t) {
332                         /* adjust throughput and move to new location */
333                         sd->P = ray->P + t*ray->D;
334
335                         return VOLUME_PATH_SCATTERED;
336                 }
337         }
338
339         return VOLUME_PATH_ATTENUATED;
340 }
341
342 /* heterogeneous volume: integrate stepping through the volume until we
343  * reach the end, get absorbed entirely, or run out of iterations */
344 ccl_device VolumeIntegrateResult kernel_volume_integrate_heterogeneous(KernelGlobals *kg,
345         PathState *state, Ray *ray, ShaderData *sd, PathRadiance *L, float3 *throughput, RNG *rng)
346 {
347         VolumeShaderCoefficients coeff;
348         float3 tp = *throughput;
349         const float tp_eps = 1e-10f; /* todo: this is likely not the right value */
350
351         /* prepare for stepping */
352         int max_steps = kernel_data.integrator.volume_max_steps;
353         float step = kernel_data.integrator.volume_step_size;
354         float random_jitter_offset = lcg_step_float(&state->rng_congruential) * step;
355
356         /* compute coefficients at the start */
357         float t = 0.0f;
358         float3 P = ray->P;
359
360         if(!volume_shader_sample(kg, sd, state, P, &coeff)) {
361                 coeff.sigma_a = make_float3(0.0f, 0.0f, 0.0f);
362                 coeff.sigma_s = make_float3(0.0f, 0.0f, 0.0f);
363                 coeff.emission = make_float3(0.0f, 0.0f, 0.0f);
364         }
365
366         /* accumulate these values so we can use a single stratified number to sample */
367         float3 accum_transmittance = make_float3(1.0f, 1.0f, 1.0f);
368         float3 accum_sigma_t = make_float3(0.0f, 0.0f, 0.0f);
369         float3 accum_sigma_s = make_float3(0.0f, 0.0f, 0.0f);
370
371         /* cache some constant variables */
372         float nlogxi;
373         int channel = -1;
374         bool has_scatter = false;
375
376         for(int i = 0; i < max_steps; i++) {
377                 /* advance to new position */
378                 float new_t = min(ray->t, t + random_jitter_offset + i * step);
379                 float3 new_P = ray->P + ray->D * new_t;
380                 VolumeShaderCoefficients new_coeff;
381
382                 /* compute segment */
383                 if(volume_shader_sample(kg, sd, state, new_P, &new_coeff)) {
384                         int closure_flag = sd->flag;
385                         float dt = new_t - t;
386                         float3 new_tp;
387                         float3 transmittance;
388                         bool scatter = false;
389
390                         /* randomly scatter, and if we do dt and new_t are shortened */
391                         if((closure_flag & SD_SCATTER) || (has_scatter && (closure_flag & SD_ABSORPTION))) {
392                                 has_scatter = true;
393
394                                 /* average sigma_t and sigma_s over segment */
395                                 float3 last_sigma_t = coeff.sigma_a + coeff.sigma_s;
396                                 float3 new_sigma_t = new_coeff.sigma_a + new_coeff.sigma_s;
397                                 float3 sigma_t = 0.5f*(last_sigma_t + new_sigma_t);
398                                 float3 sigma_s = 0.5f*(coeff.sigma_s + new_coeff.sigma_s);
399
400                                 /* lazily set up variables for sampling */
401                                 if(channel == -1) {
402                                         float xi = path_state_rng_1D(kg, rng, state, PRNG_SCATTER_DISTANCE);
403                                         nlogxi = -logf(1.0f - xi);
404
405                                         float rphase = path_state_rng_1D(kg, rng, state, PRNG_PHASE);
406                                         channel = (int)(rphase*3.0f);
407                                         sd->randb_closure = rphase*3.0f - channel;
408                                 }
409
410                                 /* pick random color channel, we use the Veach one-sample
411                                  * model with balance heuristic for the channels */
412                                 float sample_sigma_t;
413
414                                 if(channel == 0)
415                                         sample_sigma_t = accum_sigma_t.x + dt*sigma_t.x;
416                                 else if(channel == 1)
417                                         sample_sigma_t = accum_sigma_t.y + dt*sigma_t.y;
418                                 else
419                                         sample_sigma_t = accum_sigma_t.z + dt*sigma_t.z;
420
421                                 if(nlogxi < sample_sigma_t) {
422                                         /* compute sampling distance */
423                                         sample_sigma_t /= new_t;
424                                         new_t = nlogxi/sample_sigma_t;
425                                         dt = new_t - t;
426
427                                         transmittance = volume_color_attenuation(sigma_t, dt);
428
429                                         accum_transmittance *= transmittance;
430                                         accum_sigma_t = (accum_sigma_t + dt*sigma_t)/new_t;
431                                         accum_sigma_s = (accum_sigma_s + dt*sigma_s)/new_t;
432
433                                         /* todo: it's not clear to me that this is correct if we move
434                                          * through a color volumed, needs verification */
435                                         float pdf = dot(accum_sigma_t, accum_transmittance);
436                                         new_tp = tp * accum_sigma_s * transmittance * (3.0f / pdf);
437
438                                         scatter = true;
439                                 }
440                                 else {
441                                         transmittance = volume_color_attenuation(sigma_t, dt);
442
443                                         accum_transmittance *= transmittance;
444                                         accum_sigma_t += dt*sigma_t;
445                                         accum_sigma_s += dt*sigma_s;
446
447                                         new_tp = tp * transmittance;
448                                 }
449                         }
450                         else if(closure_flag & SD_ABSORPTION) {
451                                 /* absorption only, no sampling needed */
452                                 float3 sigma_a = 0.5f*(coeff.sigma_a + new_coeff.sigma_a);
453                                 transmittance = volume_color_attenuation(sigma_a, dt);
454
455                                 accum_transmittance *= transmittance;
456                                 accum_sigma_t += dt*sigma_a;
457
458                                 new_tp = tp * transmittance;
459
460                                 /* todo: we could avoid computing expf() for each step by summing,
461                                  * because exp(a)*exp(b) = exp(a+b), but we still want a quick
462                                  * tp_eps check too */
463                         }
464
465                         /* integrate emission attenuated by absorption 
466                          * integral E * exp(-sigma_t * t) from 0 to t = E * (1 - exp(-sigma_t * t))/sigma_t
467                          * this goes to E * t as sigma_t goes to zero
468                          *
469                          * todo: we should use an epsilon to avoid precision issues near zero sigma_t */
470                         if(closure_flag & SD_EMISSION) {
471                                 float3 emission = 0.5f*(coeff.emission + new_coeff.emission);
472
473                                 if(closure_flag & SD_ABSORPTION) {
474                                         float3 sigma_t = 0.5f*(coeff.sigma_a + coeff.sigma_s + new_coeff.sigma_a + new_coeff.sigma_s);
475
476                                         emission.x *= (sigma_t.x > 0.0f)? (1.0f - transmittance.x)/sigma_t.x: dt;
477                                         emission.y *= (sigma_t.y > 0.0f)? (1.0f - transmittance.y)/sigma_t.y: dt;
478                                         emission.z *= (sigma_t.z > 0.0f)? (1.0f - transmittance.z)/sigma_t.z: dt;
479                                 }
480                                 else
481                                         emission *= dt;
482
483                                 path_radiance_accum_emission(L, tp, emission, state->bounce);
484                         }
485
486                         /* modify throughput */
487                         if(closure_flag & (SD_ABSORPTION|SD_SCATTER)) {
488                                 tp = new_tp;
489
490                                 /* stop if nearly all light blocked */
491                                 if(tp.x < tp_eps && tp.y < tp_eps && tp.z < tp_eps) {
492                                         tp = make_float3(0.0f, 0.0f, 0.0f);
493                                         break;
494                                 }
495
496                                 /* prepare to scatter to new direction */
497                                 if(scatter) {
498                                         /* adjust throughput and move to new location */
499                                         sd->P = ray->P + new_t*ray->D;
500                                         *throughput = tp;
501
502                                         return VOLUME_PATH_SCATTERED;
503                                 }
504                         }
505
506                         coeff = new_coeff;
507                 }
508                 else {
509                         /* skip empty space */
510                         coeff.sigma_a = make_float3(0.0f, 0.0f, 0.0f);
511                         coeff.sigma_s = make_float3(0.0f, 0.0f, 0.0f);
512                         coeff.emission = make_float3(0.0f, 0.0f, 0.0f);
513                 }
514
515                 /* stop if at the end of the volume */
516                 t = new_t;
517                 if(t == ray->t)
518                         break;
519         }
520
521         /* include pdf for volumes with scattering */
522         if(has_scatter) {
523                 float pdf = (accum_transmittance.x + accum_transmittance.y + accum_transmittance.z);
524                 if(pdf > 0.0f)
525                         tp *= (3.0f/pdf);
526         }
527
528         *throughput = tp;
529
530         return VOLUME_PATH_ATTENUATED;
531 }
532
533 /* get the volume attenuation and emission over line segment defined by
534  * ray, with the assumption that there are no surfaces blocking light
535  * between the endpoints */
536 ccl_device_noinline VolumeIntegrateResult kernel_volume_integrate(KernelGlobals *kg,
537         PathState *state, ShaderData *sd, Ray *ray, PathRadiance *L, float3 *throughput, RNG *rng)
538 {
539         /* workaround to fix correlation bug in T38710, can find better solution
540          * in random number generator later, for now this is done here to not impact
541          * performance of rendering without volumes */
542         RNG tmp_rng = cmj_hash(*rng, state->rng_offset);
543
544         shader_setup_from_volume(kg, sd, ray, state->bounce);
545
546         if(volume_stack_is_heterogeneous(kg, state->volume_stack))
547                 return kernel_volume_integrate_heterogeneous(kg, state, ray, sd, L, throughput, &tmp_rng);
548         else
549                 return kernel_volume_integrate_homogeneous(kg, state, ray, sd, L, throughput, &tmp_rng);
550 }
551
552 /* Volume Stack
553  *
554  * This is an array of object/shared ID's that the current segment of the path
555  * is inside of. */
556
557 ccl_device void kernel_volume_stack_init(KernelGlobals *kg, VolumeStack *stack)
558 {
559         /* todo: this assumes camera is always in air, need to detect when it isn't */
560         if(kernel_data.background.volume_shader == SHADER_NO_ID) {
561                 stack[0].shader = SHADER_NO_ID;
562         }
563         else {
564                 stack[0].shader = kernel_data.background.volume_shader;
565                 stack[0].object = ~0;
566                 stack[1].shader = SHADER_NO_ID;
567         }
568 }
569
570 ccl_device void kernel_volume_stack_enter_exit(KernelGlobals *kg, ShaderData *sd, VolumeStack *stack)
571 {
572         /* todo: we should have some way for objects to indicate if they want the
573          * world shader to work inside them. excluding it by default is problematic
574          * because non-volume objects can't be assumed to be closed manifolds */
575
576         if(!(sd->flag & SD_HAS_VOLUME))
577                 return;
578         
579         if(sd->flag & SD_BACKFACING) {
580                 /* exit volume object: remove from stack */
581                 for(int i = 0; stack[i].shader != SHADER_NO_ID; i++) {
582                         if(stack[i].object == sd->object) {
583                                 /* shift back next stack entries */
584                                 do {
585                                         stack[i] = stack[i+1];
586                                         i++;
587                                 }
588                                 while(stack[i].shader != SHADER_NO_ID);
589
590                                 return;
591                         }
592                 }
593         }
594         else {
595                 /* enter volume object: add to stack */
596                 int i;
597
598                 for(i = 0; stack[i].shader != SHADER_NO_ID; i++) {
599                         /* already in the stack? then we have nothing to do */
600                         if(stack[i].object == sd->object)
601                                 return;
602                 }
603
604                 /* if we exceed the stack limit, ignore */
605                 if(i >= VOLUME_STACK_SIZE-1)
606                         return;
607
608                 /* add to the end of the stack */
609                 stack[i].shader = sd->shader;
610                 stack[i].object = sd->object;
611                 stack[i+1].shader = SHADER_NO_ID;
612         }
613 }
614
615 CCL_NAMESPACE_END
616