intern/cycles/kernel/closure/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 #ifndef __VOLUME_H__
  18 #define __VOLUME_H__
  19
  20 CCL_NAMESPACE_BEGIN
  21
  22 /* VOLUME EXTINCTION */
  23
  24 ccl_device void volume_extinction_setup(ShaderData *sd, float3 weight)
  25 {
  26         if(sd->flag & SD_EXTINCTION) {
  27                 sd->closure_transparent_extinction += weight;
  28         }
  29         else {
  30                 sd->flag |= SD_EXTINCTION;
  31                 sd->closure_transparent_extinction = weight;
  32         }
  33 }
  34
  35 /* HENYEY-GREENSTEIN CLOSURE */
  36
  37 typedef ccl_addr_space struct HenyeyGreensteinVolume {
  38         SHADER_CLOSURE_BASE;
  39
  40         float g;
  41 } HenyeyGreensteinVolume;
  42
  43 /* Given cosine between rays, return probability density that a photon bounces
  44  * to that direction. The g parameter controls how different it is from the
  45  * uniform sphere. g=0 uniform diffuse-like, g=1 close to sharp single ray. */
  46 ccl_device float single_peaked_henyey_greenstein(float cos_theta, float g)
  47 {
  48         return ((1.0f - g * g) / safe_powf(1.0f + g * g - 2.0f * g * cos_theta, 1.5f)) * (M_1_PI_F * 0.25f);
  49 };
  50
  51 ccl_device int volume_henyey_greenstein_setup(HenyeyGreensteinVolume *volume)
  52 {
  53         volume->type = CLOSURE_VOLUME_HENYEY_GREENSTEIN_ID;
  54
  55         /* clamp anisotropy to avoid delta function */
  56         volume->g = signf(volume->g) * min(fabsf(volume->g), 1.0f - 1e-3f);
  57
  58         return SD_SCATTER;
  59 }
  60
  61 ccl_device bool volume_henyey_greenstein_merge(const ShaderClosure *a, const ShaderClosure *b)
  62 {
  63         const HenyeyGreensteinVolume *volume_a = (const HenyeyGreensteinVolume*)a;
  64         const HenyeyGreensteinVolume *volume_b = (const HenyeyGreensteinVolume*)b;
  65
  66         return (volume_a->g == volume_b->g);
  67 }
  68
  69 ccl_device float3 volume_henyey_greenstein_eval_phase(const ShaderClosure *sc, const float3 I, float3 omega_in, float *pdf)
  70 {
  71         const HenyeyGreensteinVolume *volume = (const HenyeyGreensteinVolume*)sc;
  72         float g = volume->g;
  73
  74         /* note that I points towards the viewer */
  75         if(fabsf(g) < 1e-3f) {
  76                 *pdf = M_1_PI_F * 0.25f;
  77         }
  78         else {
  79                 float cos_theta = dot(-I, omega_in);
  80                 *pdf = single_peaked_henyey_greenstein(cos_theta, g);
  81         }
  82
  83         return make_float3(*pdf, *pdf, *pdf);
  84 }
  85
  86 ccl_device float3 henyey_greenstrein_sample(float3 D, float g, float randu, float randv, float *pdf)
  87 {
  88         /* match pdf for small g */
  89         float cos_theta;
  90         bool isotropic = fabsf(g) < 1e-3f;
  91
  92         if(isotropic) {
  93                 cos_theta = (1.0f - 2.0f * randu);
  94                 if(pdf) {
  95                         *pdf = M_1_PI_F * 0.25f;
  96                 }
  97         }
  98         else {
  99                 float k = (1.0f - g * g) / (1.0f - g + 2.0f * g * randu);
 100                 cos_theta = (1.0f + g * g - k * k) / (2.0f * g);
 101                 if(pdf) {
 102                         *pdf = single_peaked_henyey_greenstein(cos_theta, g);
 103                 }
 104         }
 105
 106         float sin_theta = safe_sqrtf(1.0f - cos_theta * cos_theta);
 107         float phi = M_2PI_F * randv;
 108         float3 dir = make_float3(sin_theta * cosf(phi), sin_theta * sinf(phi), cos_theta);
 109
 110         float3 T, B;
 111         make_orthonormals(D, &T, &B);
 112         dir = dir.x * T + dir.y * B + dir.z * D;
 113
 114         return dir;
 115 }
 116
 117 ccl_device int volume_henyey_greenstein_sample(const ShaderClosure *sc, float3 I, float3 dIdx, float3 dIdy, float randu, float randv,
 118         float3 *eval, float3 *omega_in, float3 *domega_in_dx, float3 *domega_in_dy, float *pdf)
 119 {
 120         const HenyeyGreensteinVolume *volume = (const HenyeyGreensteinVolume*)sc;
 121         float g = volume->g;
 122
 123         /* note that I points towards the viewer and so is used negated */
 124         *omega_in = henyey_greenstrein_sample(-I, g, randu, randv, pdf);
 125         *eval = make_float3(*pdf, *pdf, *pdf); /* perfect importance sampling */
 126
 127 #ifdef __RAY_DIFFERENTIALS__
 128         /* todo: implement ray differential estimation */
 129         *domega_in_dx = make_float3(0.0f, 0.0f, 0.0f);
 130         *domega_in_dy = make_float3(0.0f, 0.0f, 0.0f);
 131 #endif
 132
 133         return LABEL_VOLUME_SCATTER;
 134 }
 135
 136 /* VOLUME CLOSURE */
 137
 138 ccl_device float3 volume_phase_eval(const ShaderData *sd, const ShaderClosure *sc, float3 omega_in, float *pdf)
 139 {
 140         kernel_assert(sc->type == CLOSURE_VOLUME_HENYEY_GREENSTEIN_ID);
 141
 142         return volume_henyey_greenstein_eval_phase(sc, sd->I, omega_in, pdf);
 143 }
 144
 145 ccl_device int volume_phase_sample(const ShaderData *sd, const ShaderClosure *sc, float randu,
 146         float randv, float3 *eval, float3 *omega_in, differential3 *domega_in, float *pdf)
 147 {
 148         int label;
 149
 150         switch(sc->type) {
 151                 case CLOSURE_VOLUME_HENYEY_GREENSTEIN_ID:
 152                         label = volume_henyey_greenstein_sample(sc, sd->I, sd->dI.dx, sd->dI.dy, randu, randv, eval, omega_in, &domega_in->dx, &domega_in->dy, pdf);
 153                         break;
 154                 default:
 155                         *eval = make_float3(0.0f, 0.0f, 0.0f);
 156                         label = LABEL_NONE;
 157                         break;
 158         }
 159
 160         return label;
 161 }
 162
 163 CCL_NAMESPACE_END
 164
 165 #endif