Fixed remaining syntax errors in OSL files. node_sepcomb_rgb.osl is split into 2...
[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
79         void setup()
80         {
81                 ior(m_eta);
82
83                 if (m_g >=  0.99f) m_g =  0.99f;
84                 if (m_g <= -0.99f) m_g = -0.99f;
85
86                 // eq.10
87                 float inv_eta = 1 / m_eta;
88                 float Fdr = -1.440f * inv_eta * inv_eta + 0.710 * inv_eta + 0.668f + 0.0636 * m_eta;
89                 float A = (1 + Fdr) / (1 - Fdr);
90                 // compute sigma_s, sigma_a (eq.16)
91                 Color3 alpha_prime = Color3(root_find_Rd(m_albedo[0], A),
92                                             root_find_Rd(m_albedo[1], A),
93                                             root_find_Rd(m_albedo[2], A));
94                 Color3 sigma_t_prime = Color3(m_mfp.x > 0 ? 1.0f / (m_mfp[0] * sqrtf(3 * (1 - alpha_prime[0]))) : 0.0f,
95                                               m_mfp.y > 0 ? 1.0f / (m_mfp[1] * sqrtf(3 * (1 - alpha_prime[1]))) : 0.0f,
96                                               m_mfp.z > 0 ? 1.0f / (m_mfp[2] * sqrtf(3 * (1 - alpha_prime[2]))) : 0.0f);
97                 Color3 sigma_s_prime = alpha_prime * sigma_t_prime;
98
99                 sigma_s((1.0f / (1 - m_g)) * sigma_s_prime);
100                 sigma_a(sigma_t_prime - sigma_s_prime);
101         }
102
103         bool mergeable(const ClosurePrimitive *other) const {
104                 const SubsurfaceClosure *comp = (const SubsurfaceClosure *)other;
105                 return m_g == comp->m_g && VolumeClosure::mergeable(other);
106         }
107
108         size_t memsize() const { return sizeof(*this); }
109
110         const char *name() const { return "subsurface"; }
111
112         void print_on(std::ostream &out) const {
113                 out << name() << " ()";
114         }
115
116         virtual Color3 eval_phase(const Vec3 &omega_in, const Vec3 &omega_out) const {
117                 float costheta = omega_in.dot(omega_out);
118                 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));
119                 return Color3(ph, ph, ph);
120         }
121 };
122
123
124
125 ClosureParam closure_subsurface_params[] = {
126         CLOSURE_FLOAT_PARAM(SubsurfaceClosure, m_eta),
127         CLOSURE_FLOAT_PARAM(SubsurfaceClosure, m_g),
128         CLOSURE_COLOR_PARAM(SubsurfaceClosure, m_mfp),
129         CLOSURE_COLOR_PARAM(SubsurfaceClosure, m_albedo),
130         CLOSURE_STRING_KEYPARAM("label"),
131         CLOSURE_FINISH_PARAM(SubsurfaceClosure)
132 };
133
134 CLOSURE_PREPARE(closure_subsurface_prepare, SubsurfaceClosure)
135
136 CCL_NAMESPACE_END
137