diff options
Diffstat (limited to 'intern/cycles/kernel/closure/bssrdf.h')
-rw-r--r-- | intern/cycles/kernel/closure/bssrdf.h | 217 |
1 files changed, 137 insertions, 80 deletions
diff --git a/intern/cycles/kernel/closure/bssrdf.h b/intern/cycles/kernel/closure/bssrdf.h index 486de4ca65f..23b932a91c6 100644 --- a/intern/cycles/kernel/closure/bssrdf.h +++ b/intern/cycles/kernel/closure/bssrdf.h @@ -21,130 +21,187 @@ CCL_NAMESPACE_BEGIN -__device int bssrdf_setup(ShaderClosure *sc) +__device int bssrdf_setup(ShaderClosure *sc, ClosureType type) { if(sc->data0 < BSSRDF_MIN_RADIUS) { /* revert to diffuse BSDF if radius too small */ sc->data0 = 0.0f; sc->data1 = 0.0f; - return bsdf_diffuse_setup(sc); + int flag = bsdf_diffuse_setup(sc); + sc->type = CLOSURE_BSDF_BSSRDF_ID; + return flag; } else { - /* IOR param */ - sc->data1 = max(sc->data1, 1.0f); - sc->type = CLOSURE_BSSRDF_ID; + sc->data1 = clamp(sc->data1, 0.0f, 1.0f); /* texture blur */ + sc->type = type; return SD_BSDF|SD_BSDF_HAS_EVAL|SD_BSSRDF; } } -/* Simple Cubic BSSRDF falloff */ +/* Planar Truncated Gaussian + * + * Note how this is different from the typical gaussian, this one integrates + * to 1 over the plane (where you get an extra 2*pi*x factor). We are lucky + * that integrating x*exp(-x) gives a nice closed form solution. */ + +/* paper suggests 1/12.46 which is much too small, suspect it's *12.46 */ +#define GAUSS_TRUNCATE 12.46f -__device float bssrdf_cubic(float ld, float r) +__device float bssrdf_gaussian_eval(ShaderClosure *sc, float r) { - if(ld == 0.0f) - return (r == 0.0f)? 1.0f: 0.0f; + /* integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) from 0 to Rm + * = 1 - exp(-Rm*Rm/(2*v)) */ + const float v = sc->data0; + const float Rm = sqrtf(v*GAUSS_TRUNCATE); + + if(r >= Rm) + return 0.0f; - return powf(ld - min(r, ld), 3.0f) * 4.0f/powf(ld, 4.0f); + return expf(-r*r/(2.0f*v))/(2.0f*M_PI_F*v); } -/* Original BSSRDF fallof function */ - -typedef struct BSSRDFParams { - float eta; /* index of refraction */ - float sigma_t_; /* reduced extinction coefficient */ - float sigma_tr; /* effective extinction coefficient */ - float Fdr; /* diffuse fresnel reflectance */ - float D; /* diffusion constant */ - float A; - float alpha_; /* reduced albedo */ - float zr; /* distance of virtual lightsource above surface */ - float zv; /* distance of virtual lightsource below surface */ - float ld; /* mean free path */ - float ro; /* diffuse reflectance */ -} BSSRDFParams; - -__device float bssrdf_reduced_albedo_Rd(float alpha_, float A, float ro) +__device float bssrdf_gaussian_pdf(ShaderClosure *sc, float r) { - float sq; + /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */ + const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE); + + return bssrdf_gaussian_eval(sc, r) * (1.0f/(area_truncated)); +} + +__device void bssrdf_gaussian_sample(ShaderClosure *sc, float xi, float *r, float *h) +{ + /* xi = integrate (2*pi*r * exp(-r*r/(2*v)))/(2*pi*v)) = -exp(-r^2/(2*v)) + * r = sqrt(-2*v*logf(xi)) */ + + const float v = sc->data0; + const float Rm = sqrtf(v*GAUSS_TRUNCATE); + + /* 1.0 - expf(-Rm*Rm/(2*v)) simplified */ + const float area_truncated = 1.0f - expf(-0.5f*GAUSS_TRUNCATE); + + /* r(xi) */ + const float r_squared = -2.0f*v*logf(1.0f - xi*area_truncated); + *r = sqrtf(r_squared); + + /* h^2 + r^2 = Rm^2 */ + *h = sqrtf(Rm*Rm - r_squared); +} + +/* Planar Cubic BSSRDF falloff + * + * This is basically (Rm - x)^3, with some factors to normalize it. For sampling + * we integrate 2*pi*x * (Rm - x)^3, which gives us a quintic equation that as + * far as I can tell has no closed form solution. So we get an iterative solution + * instead with newton-raphson. */ + +__device float bssrdf_cubic_eval(ShaderClosure *sc, float r) +{ + const float Rm = sc->data0; + + if(r >= Rm) + return 0.0f; + + /* integrate (2*pi*r * 10*(R - r)^3)/(pi * R^5) from 0 to R = 1 */ + const float Rm5 = (Rm*Rm) * (Rm*Rm) * Rm; + const float f = Rm - min(r, Rm); + const float f3 = f*f*f; - sq = sqrtf(3.0f*(1.0f - alpha_)); - return (alpha_/2.0f)*(1.0f + expf((-4.0f/3.0f)*A*sq))*expf(-sq) - ro; + return (f3 * 10.0f) / (Rm5 * M_PI_F); } -__device float bssrdf_compute_reduced_albedo(float A, float ro) +__device float bssrdf_cubic_pdf(ShaderClosure *sc, float r) { - const float tolerance = 1e-8f; - const int max_iteration_count = 20; - float d, fsub, xn_1 = 0.0f, xn = 1.0f, fxn, fxn_1; + return bssrdf_cubic_eval(sc, r); +} + +/* solve 10x^2 - 20x^3 + 15x^4 - 4x^5 - xi == 0 */ +__device float bssrdf_cubic_quintic_root_find(float xi) +{ + /* newton-raphson iteration, usually succeeds in 2-4 iterations, except + * outside 0.02 ... 0.98 where it can go up to 10, so overall performance + * should not be too bad */ + const float tolerance = 1e-6f; + const int max_iteration_count = 10; + float x = 0.25f; int i; - /* use secant method to compute reduced albedo using Rd function inverse - * with a given reflectance */ - fxn = bssrdf_reduced_albedo_Rd(xn, A, ro); - fxn_1 = bssrdf_reduced_albedo_Rd(xn_1, A, ro); + for (i = 0; i < max_iteration_count; i++) { + float x2 = x*x; + float x3 = x2*x; + float nx = (1.0f - x); - for (i= 0; i < max_iteration_count; i++) { - fsub = (fxn - fxn_1); - if (fabsf(fsub) < tolerance) - break; - d = ((xn - xn_1)/fsub)*fxn; - if (fabsf(d) < tolerance) - break; + float f = 10.0f*x2 - 20.0f*x3 + 15.0f*x2*x2 - 4.0f*x2*x3 - xi; + float f_ = 20.0f*(x*nx)*(nx*nx); - xn_1 = xn; - fxn_1 = fxn; - xn = xn - d; + if(fabsf(f) < tolerance || f_ == 0.0f) + break; - if (xn > 1.0f) xn = 1.0f; - if (xn_1 > 1.0f) xn_1 = 1.0f; - - fxn = bssrdf_reduced_albedo_Rd(xn, A, ro); + x = clamp(x - f/f_, 0.0f, 1.0f); } - /* avoid division by zero later */ - if (xn <= 0.0f) - xn = 0.00001f; - - return xn; + return x; } -__device void bssrdf_setup_params(BSSRDFParams *ss, float refl, float radius, float ior) +__device void bssrdf_cubic_sample(ShaderClosure *sc, float xi, float *r, float *h) { - ss->eta = ior; - ss->Fdr = -1.440f/ior*ior + 0.710f/ior + 0.668f + 0.0636f*ior; - ss->A = (1.0f + ss->Fdr)/(1.0f - ss->Fdr); - ss->ld = radius; - ss->ro = min(refl, 0.999f); + const float Rm = sc->data0; + const float r_ = bssrdf_cubic_quintic_root_find(xi) * Rm; - ss->alpha_ = bssrdf_compute_reduced_albedo(ss->A, ss->ro); + *r = r_; - ss->sigma_tr = 1.0f/ss->ld; - ss->sigma_t_ = ss->sigma_tr/sqrtf(3.0f*(1.0f - ss->alpha_)); + /* h^2 + r^2 = Rm^2 */ + *h = sqrtf(Rm*Rm - r_*r_); +} - ss->D = 1.0f/(3.0f*ss->sigma_t_); +/* None BSSRDF falloff + * + * Samples distributed over disk with no falloff, for reference. */ - ss->zr = 1.0f/ss->sigma_t_; - ss->zv = ss->zr + 4.0f*ss->A*ss->D; +__device float bssrdf_none_eval(ShaderClosure *sc, float r) +{ + const float Rm = sc->data0; + return (r < Rm)? 1.0f: 0.0f; } -/* exponential falloff function */ +__device float bssrdf_none_pdf(ShaderClosure *sc, float r) +{ + /* integrate (2*pi*r)/(pi*Rm*Rm) from 0 to Rm = 1 */ + const float Rm = sc->data0; + const float area = (M_PI_F*Rm*Rm); + + return bssrdf_none_eval(sc, r) / area; +} -__device float bssrdf_original(const BSSRDFParams *ss, float r) +__device void bssrdf_none_sample(ShaderClosure *sc, float xi, float *r, float *h) { - if(ss->ld == 0.0f) - return (r == 0.0f)? 1.0f: 0.0f; + /* xi = integrate (2*pi*r)/(pi*Rm*Rm) = r^2/Rm^2 + * r = sqrt(xi)*Rm */ + const float Rm = sc->data0; + const float r_ = sqrtf(xi)*Rm; + + *r = r_; - float rr = r*r; - float sr, sv, Rdr, Rdv; + /* h^2 + r^2 = Rm^2 */ + *h = sqrtf(Rm*Rm - r_*r_); +} - sr = sqrtf(rr + ss->zr*ss->zr); - sv = sqrtf(rr + ss->zv*ss->zv); +/* Generic */ - Rdr = ss->zr*(1.0f + ss->sigma_tr*sr)*expf(-ss->sigma_tr*sr)/(sr*sr*sr); - Rdv = ss->zv*(1.0f + ss->sigma_tr*sv)*expf(-ss->sigma_tr*sv)/(sv*sv*sv); +__device void bssrdf_sample(ShaderClosure *sc, float xi, float *r, float *h) +{ + if(sc->type == CLOSURE_BSSRDF_CUBIC_ID) + bssrdf_cubic_sample(sc, xi, r, h); + else + bssrdf_gaussian_sample(sc, xi, r, h); +} - return ss->alpha_*(1.0f/M_4PI_F)*(Rdr + Rdv); +__device float bssrdf_pdf(ShaderClosure *sc, float r) +{ + if(sc->type == CLOSURE_BSSRDF_CUBIC_ID) + return bssrdf_cubic_pdf(sc, r); + else + return bssrdf_gaussian_pdf(sc, r); } CCL_NAMESPACE_END |