Cycles: svn merge -r41225:41232 ^/trunk/blender
[blender.git] / intern / cycles / kernel / osl / vol_subsurface.cpp
1 /*
2  * Adapted from Open Shading Language with this license:
3  *
4  * Copyright (c) 2009-2010 Sony Pictures Imageworks Inc., et al.
5  * All Rights Reserved.
6  *
7  * Modifications Copyright 2011, Blender Foundation.
8  * 
9  * Redistribution and use in source and binary forms, with or without
10  * modification, are permitted provided that the following conditions are
11  * met:
12  * * Redistributions of source code must retain the above copyright
13  *   notice, this list of conditions and the following disclaimer.
14  * * Redistributions in binary form must reproduce the above copyright
15  *   notice, this list of conditions and the following disclaimer in the
16  *   documentation and/or other materials provided with the distribution.
17  * * Neither the name of Sony Pictures Imageworks nor the names of its
18  *   contributors may be used to endorse or promote products derived from
19  *   this software without specific prior written permission.
20  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23  * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24  * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31  */
32
33 #include <OpenImageIO/fmath.h>
34
35 #include <OSL/genclosure.h>
36
37 #include "osl_closures.h"
38
39 CCL_NAMESPACE_BEGIN
40
41 using namespace OSL;
42
43 // Computes scattering properties based on Jensen's reparameterization
44 // described in:
45 //    http://graphics.ucsd.edu/~henrik/papers/fast_bssrdf/
46
47 class SubsurfaceClosure : public VolumeClosure {
48 public:
49     float  m_g;
50     float  m_eta;
51     Color3 m_mfp, m_albedo;
52     static float root_find_Rd(const float Rd0, const float A) {
53         // quick exit for trivial cases
54         if (Rd0 <= 0) return 0;
55         const float A43 = A * 4.0f / 3.0f;
56         // Find alpha such that f(alpha) = Rd (see eq.15). A simple bisection
57         // method can be used because this function is monotonicaly increasing.
58         float lo = 0, hi = 1;
59         for (int i = 0; i < 20; i++) { // 2^20 divisions should be sufficient
60             // eval function at midpoint
61             float alpha = 0.5f * (lo + hi);
62             float a1 = sqrtf(3 * (1 - alpha));
63             float e1 = expf(-a1);
64             float e2 = expf(-A43 * a1);
65             float Rd = 0.5f * alpha * (1 + e2) * e1 - Rd0;
66             if (fabsf(Rd) < 1e-6f)
67                 return alpha; // close enough
68             else if (Rd > 0)
69                 hi = alpha; // root is on left side
70             else
71                 lo = alpha; // root is on right side
72         }
73         // didn't quite converge, pick result in the middle of remaining interval
74         return 0.5f * (lo + hi);
75     }
76     SubsurfaceClosure() { }
77
78     void setup()
79     {
80         ior(m_eta);
81
82         if (m_g >=  0.99f) m_g =  0.99f;
83         if (m_g <= -0.99f) m_g = -0.99f;
84
85         // eq.10
86         float inv_eta = 1 / m_eta;
87         float Fdr = -1.440f * inv_eta * inv_eta + 0.710 * inv_eta + 0.668f + 0.0636 * m_eta;
88         float A = (1 + Fdr) / (1 - Fdr);
89         // compute sigma_s, sigma_a (eq.16)
90         Color3 alpha_prime = Color3 (root_find_Rd(m_albedo[0], A),
91                                      root_find_Rd(m_albedo[1], A),
92                                      root_find_Rd(m_albedo[2], A));
93         Color3 sigma_t_prime = Color3 (m_mfp.x > 0 ? 1.0f / (m_mfp[0] * sqrtf(3 * (1 - alpha_prime[0]))) : 0.0f,
94                                        m_mfp.y > 0 ? 1.0f / (m_mfp[1] * sqrtf(3 * (1 - alpha_prime[1]))) : 0.0f,
95                                        m_mfp.z > 0 ? 1.0f / (m_mfp[2] * sqrtf(3 * (1 - alpha_prime[2]))) : 0.0f);
96         Color3 sigma_s_prime = alpha_prime * sigma_t_prime;
97
98         sigma_s((1.0f / (1 - m_g)) * sigma_s_prime);
99         sigma_a(sigma_t_prime - sigma_s_prime);
100     }
101
102     bool mergeable (const ClosurePrimitive *other) const {
103         const SubsurfaceClosure *comp = (const SubsurfaceClosure *)other;
104         return m_g == comp->m_g && VolumeClosure::mergeable(other);
105     }
106
107     size_t memsize () const { return sizeof(*this); }
108
109     const char *name () const { return "subsurface"; }
110
111     void print_on (std::ostream &out) const {
112         out << name() << " ()";
113     }
114
115     virtual Color3 eval_phase(const Vec3 &omega_in, const Vec3 &omega_out) const {
116         float costheta = omega_in.dot(omega_out);
117         float ph = 0.25f * float(M_1_PI) * ((1 - m_g * m_g) / powf(1 + m_g * m_g - 2.0f * m_g * costheta, 1.5f));
118         return Color3 (ph, ph, ph);
119     }
120 };
121
122
123
124 ClosureParam closure_subsurface_params[] = {
125     CLOSURE_FLOAT_PARAM (SubsurfaceClosure, m_eta),
126     CLOSURE_FLOAT_PARAM (SubsurfaceClosure, m_g),
127     CLOSURE_COLOR_PARAM (SubsurfaceClosure, m_mfp),
128     CLOSURE_COLOR_PARAM (SubsurfaceClosure, m_albedo),
129     CLOSURE_STRING_KEYPARAM("label"),
130     CLOSURE_FINISH_PARAM(SubsurfaceClosure) };
131
132 CLOSURE_PREPARE(closure_subsurface_prepare, SubsurfaceClosure)
133
134 CCL_NAMESPACE_END
135