2bd06093c7
Applies thin film iridescence to metals in Metallic BSDF and Principled BSDF. To get the complex IOR values for each spectral band from F82 Tint colors, the code uses the parametrization from "Artist Friendly Metallic Fresnel", where the g parameter is set to F82. This IOR is used to find the phase shift, but reflectance is still calculated with the F82 Tint formula after adjusting F0 for the film's IOR. Co-authored-by: Lukas Stockner <lukas@lukasstockner.de> Co-authored-by: Weizhen Huang <weizhen@blender.org> Co-authored-by: RobertMoerland <rmoerlandrj@gmail.com> Pull Request: https://projects.blender.org/blender/blender/pulls/141131
568 lines
23 KiB
C++
568 lines
23 KiB
C++
/* SPDX-FileCopyrightText: 2009-2010 Sony Pictures Imageworks Inc., et al. All Rights Reserved.
|
|
* SPDX-FileCopyrightText: 2011-2022 Blender Foundation
|
|
*
|
|
* SPDX-License-Identifier: BSD-3-Clause
|
|
*
|
|
* Adapted code from Open Shading Language. */
|
|
|
|
#pragma once
|
|
|
|
#include "kernel/types.h"
|
|
|
|
#include "kernel/util/colorspace.h"
|
|
#include "kernel/util/lookup_table.h"
|
|
|
|
#include "util/color.h"
|
|
|
|
CCL_NAMESPACE_BEGIN
|
|
|
|
struct FresnelThinFilm {
|
|
float thickness;
|
|
float ior;
|
|
};
|
|
|
|
template<typename T> struct complex {
|
|
T re;
|
|
T im;
|
|
|
|
ccl_device_inline_method complex<T> operator*=(ccl_private const complex<T> &other)
|
|
{
|
|
const T im = this->re * other.im + this->im * other.re;
|
|
this->re = this->re * other.re - this->im * other.im;
|
|
this->im = im;
|
|
return *this;
|
|
}
|
|
|
|
ccl_device_inline_method complex<T> operator*(ccl_private const float &other)
|
|
{
|
|
return complex<T>{this->re * other, this->im * other};
|
|
}
|
|
};
|
|
|
|
/* Compute fresnel reflectance for perpendicular (aka S-) and parallel (aka P-) polarized light.
|
|
* If requested by the caller, r_cos_phi is set to the cosine of the phase shift on reflection.
|
|
* Also returns the dot product of the refracted ray and the normal as `cos_theta_t`, as it is
|
|
* used when computing the direction of the refracted ray. */
|
|
ccl_device float2 fresnel_dielectric_polarized(float cos_theta_i,
|
|
const float eta,
|
|
ccl_private float *r_cos_theta_t,
|
|
ccl_private float2 *r_cos_phi)
|
|
{
|
|
kernel_assert(!isnan_safe(cos_theta_i));
|
|
|
|
/* Using Snell's law, calculate the squared cosine of the angle between the surface normal and
|
|
* the transmitted ray. */
|
|
const float eta_cos_theta_t_sq = sqr(eta) - (1.0f - sqr(cos_theta_i));
|
|
if (eta_cos_theta_t_sq <= 0) {
|
|
/* Total internal reflection. */
|
|
if (r_cos_phi) {
|
|
/* The following code would compute the proper phase shift on TIR.
|
|
* However, for the current user of this computation (the iridescence code),
|
|
* this doesn't actually affect the result, so don't bother with the computation for now.
|
|
*
|
|
* `const float fac = sqrtf(1.0f - sqr(cosThetaI) - sqr(eta));`
|
|
* `r_phi->x = -2.0f * atanf(fac / cosThetaI);`
|
|
* `r_phi->y = -2.0f * atanf(fac / (cosThetaI * sqr(eta)));`
|
|
*/
|
|
*r_cos_phi = one_float2();
|
|
}
|
|
return one_float2();
|
|
}
|
|
|
|
cos_theta_i = fabsf(cos_theta_i);
|
|
/* Relative to the surface normal. */
|
|
const float cos_theta_t = -safe_sqrtf(eta_cos_theta_t_sq) / eta;
|
|
|
|
if (r_cos_theta_t) {
|
|
*r_cos_theta_t = cos_theta_t;
|
|
}
|
|
|
|
/* Amplitudes of reflected waves. */
|
|
const float r_s = (cos_theta_i + eta * cos_theta_t) / (cos_theta_i - eta * cos_theta_t);
|
|
const float r_p = (cos_theta_t + eta * cos_theta_i) / (eta * cos_theta_i - cos_theta_t);
|
|
|
|
if (r_cos_phi) {
|
|
*r_cos_phi = make_float2(2 * (r_s >= 0.0f) - 1, 2 * (r_p >= 0.0f) - 1);
|
|
}
|
|
|
|
/* Return squared amplitude to get the fraction of reflected energy. */
|
|
return make_float2(sqr(r_s), sqr(r_p));
|
|
}
|
|
|
|
/* Compute fresnel reflectance for unpolarized light. */
|
|
ccl_device_forceinline float fresnel_dielectric(const float cos_theta_i,
|
|
const float eta,
|
|
ccl_private float *r_cos_theta_t)
|
|
{
|
|
return average(fresnel_dielectric_polarized(cos_theta_i, eta, r_cos_theta_t, nullptr));
|
|
}
|
|
|
|
/* Refract the incident ray, given the cosine of the refraction angle and the relative refractive
|
|
* index of the incoming medium w.r.t. the outgoing medium. */
|
|
ccl_device_inline float3 refract_angle(const float3 incident,
|
|
const float3 normal,
|
|
const float cos_theta_t,
|
|
const float inv_eta)
|
|
{
|
|
return (inv_eta * dot(normal, incident) + cos_theta_t) * normal - inv_eta * incident;
|
|
}
|
|
|
|
ccl_device float fresnel_dielectric_cos(const float cosi, const float eta)
|
|
{
|
|
// compute fresnel reflectance without explicitly computing
|
|
// the refracted direction
|
|
const float c = fabsf(cosi);
|
|
float g = eta * eta - 1 + c * c;
|
|
if (g > 0) {
|
|
g = sqrtf(g);
|
|
const float A = (g - c) / (g + c);
|
|
const float B = (c * (g + c) - 1) / (c * (g - c) + 1);
|
|
return 0.5f * A * A * (1 + B * B);
|
|
}
|
|
return 1.0f; // TIR(no refracted component)
|
|
}
|
|
|
|
/* Approximates the average single-scattering Fresnel for a given IOR.
|
|
* This is defined as the integral over 0...1 of 2*cosI * F(cosI, eta) d_cosI, with F being
|
|
* the real dielectric Fresnel.
|
|
* The implementation here uses a numerical fit from "Revisiting Physically Based Shading
|
|
* at Imageworks" by Christopher Kulla and Alejandro Conty. */
|
|
ccl_device_inline float fresnel_dielectric_Fss(const float eta)
|
|
{
|
|
if (eta < 1.0f) {
|
|
return 0.997118f + eta * (0.1014f - eta * (0.965241f + eta * 0.130607f));
|
|
}
|
|
return (eta - 1.0f) / (4.08567f + 1.00071f * eta);
|
|
}
|
|
|
|
/* Evaluates the Fresnel equations at a dielectric-conductor interface, calculating reflectances
|
|
* and phase shifts due to reflection if requested. The phase shifts phi_s and phi_p are returned
|
|
* as phasor_s = exp(i * phi_s) and phasor_p = exp(i * phi_p).
|
|
* This code is based on equations from section 14.4.1 of Principles of Optics 7th ed. by Born and
|
|
* Wolf, but uses `n + ik` instead of `n(1 + ik)` for IOR. The phase shifts are calculated so that
|
|
* phi_p = phi_s at 90 degree incidence to match fresnel_dielectric_polarized. */
|
|
ccl_device void fresnel_conductor_polarized(const float cosi,
|
|
const float ambient_ior,
|
|
const complex<Spectrum> conductor_ior,
|
|
ccl_private Spectrum *r_R_s,
|
|
ccl_private Spectrum *r_R_p,
|
|
ccl_private complex<Spectrum> *r_phasor_s,
|
|
ccl_private complex<Spectrum> *r_phasor_p)
|
|
{
|
|
const float eta1 = ambient_ior;
|
|
const Spectrum eta2 = conductor_ior.re;
|
|
const Spectrum k2 = conductor_ior.im;
|
|
|
|
const float eta1_sq = sqr(eta1);
|
|
const Spectrum eta2_sq = sqr(eta2);
|
|
const Spectrum k2_sq = sqr(k2);
|
|
const Spectrum two_eta2_k2 = 2.0f * eta2 * k2;
|
|
|
|
const Spectrum t1 = eta2_sq - k2_sq - eta1_sq * (1.0f - sqr(cosi));
|
|
const Spectrum t2 = sqrt(sqr(t1) + sqr(two_eta2_k2));
|
|
|
|
const Spectrum u_sq = max(0.5f * (t2 + t1), zero_spectrum());
|
|
const Spectrum v_sq = max(0.5f * (t2 - t1), zero_spectrum());
|
|
const Spectrum u = sqrt(u_sq);
|
|
const Spectrum v = sqrt(v_sq);
|
|
|
|
if (r_R_s && r_R_p) {
|
|
*r_R_s = safe_divide(sqr(eta1 * cosi - u) + v_sq, sqr(eta1 * cosi + u) + v_sq);
|
|
|
|
const Spectrum t3 = (eta2_sq - k2_sq) * cosi;
|
|
const Spectrum t4 = two_eta2_k2 * cosi;
|
|
*r_R_p = safe_divide(sqr(t3 - eta1 * u) + sqr(t4 - eta1 * v),
|
|
sqr(t3 + eta1 * u) + sqr(t4 + eta1 * v));
|
|
}
|
|
|
|
if (r_phasor_s && r_phasor_p) {
|
|
const Spectrum re_s = -u_sq - v_sq + sqr(eta1 * cosi);
|
|
const Spectrum im_s = -2.0f * eta1 * cosi * v;
|
|
const Spectrum mag_s = sqrt(sqr(re_s) + sqr(im_s));
|
|
r_phasor_s->re = select(is_zero_mask(mag_s), one_spectrum(), re_s / mag_s);
|
|
r_phasor_s->im = select(is_zero_mask(mag_s), zero_spectrum(), im_s / mag_s);
|
|
|
|
const Spectrum re_p = sqr((eta2_sq + k2_sq) * cosi) - eta1_sq * (u_sq + v_sq);
|
|
const Spectrum im_p = 2.0f * eta1 * cosi * (two_eta2_k2 * u - (eta2_sq - k2_sq) * v);
|
|
const Spectrum mag_p = sqrt(sqr(re_p) + sqr(im_p));
|
|
r_phasor_p->re = select(is_zero_mask(mag_p), one_spectrum(), re_p / mag_p);
|
|
r_phasor_p->im = select(is_zero_mask(mag_p), zero_spectrum(), im_p / mag_p);
|
|
}
|
|
}
|
|
|
|
/* Calculates Fresnel reflectance at a dielectric-conductor interface given the relative IOR. */
|
|
ccl_device Spectrum fresnel_conductor(const float cosi, const complex<Spectrum> ior)
|
|
{
|
|
Spectrum R_s, R_p;
|
|
fresnel_conductor_polarized(cosi, 1.0f, ior, &R_s, &R_p, nullptr, nullptr);
|
|
return (R_s + R_p) * 0.5f;
|
|
}
|
|
|
|
/* Computes the average single-scattering Fresnel for the F82 metallic model. */
|
|
ccl_device_inline Spectrum fresnel_f82_Fss(const Spectrum F0, const Spectrum B)
|
|
{
|
|
return mix(F0, one_spectrum(), 1.0f / 21.0f) - B * (1.0f / 126.0f);
|
|
}
|
|
|
|
/* Precompute the B term for the F82 metallic model, given a tint factor. */
|
|
ccl_device_inline Spectrum fresnel_f82tint_B(const Spectrum F0, const Spectrum tint)
|
|
{
|
|
/* In the classic F82 model, the F82 input directly determines the value of the Fresnel
|
|
* model at ~82°, similar to F0 and F90.
|
|
* With F82-Tint, on the other hand, the value at 82° is the value of the classic Schlick
|
|
* model multiplied by the tint input.
|
|
* Therefore, the factor follows by setting F82Tint(cosI) = FSchlick(cosI) - b*cosI*(1-cosI)^6
|
|
* and F82Tint(acos(1/7)) = FSchlick(acos(1/7)) * f82_tint and solving for b. */
|
|
const float f = 6.0f / 7.0f;
|
|
const float f5 = sqr(sqr(f)) * f;
|
|
const Spectrum F_schlick = mix(F0, one_spectrum(), f5);
|
|
return F_schlick * (7.0f / (f5 * f)) * (one_spectrum() - tint);
|
|
}
|
|
|
|
/* Precompute the B term for the F82 metallic model, given the F82 value. */
|
|
ccl_device_inline Spectrum fresnel_f82_B(const Spectrum F0, const Spectrum F82)
|
|
{
|
|
const float f = 6.0f / 7.0f;
|
|
const float f5 = sqr(sqr(f)) * f;
|
|
const Spectrum F_schlick = mix(F0, one_spectrum(), f5);
|
|
return (7.0f / (f5 * f)) * (F_schlick - F82);
|
|
}
|
|
|
|
/* Evaluate the F82 metallic model for the given parameters. */
|
|
ccl_device_inline Spectrum fresnel_f82(const float cosi, const Spectrum F0, const Spectrum B)
|
|
{
|
|
const float s = saturatef(1.0f - cosi);
|
|
const float s5 = sqr(sqr(s)) * s;
|
|
const Spectrum F_schlick = mix(F0, one_spectrum(), s5);
|
|
return saturate(F_schlick - B * cosi * s5 * s);
|
|
}
|
|
|
|
/* Approximates the average single-scattering Fresnel for a physical conductor. */
|
|
ccl_device_inline Spectrum fresnel_conductor_Fss(const complex<Spectrum> ior)
|
|
{
|
|
/* In order to estimate Fss of the conductor, we fit the F82 model to it based on the
|
|
* value at 0° and ~82° and then use the analytic expression for its Fss. */
|
|
const Spectrum F0 = fresnel_conductor(1.0f, ior);
|
|
const Spectrum F82 = fresnel_conductor(1.0f / 7.0f, ior);
|
|
return saturate(fresnel_f82_Fss(F0, fresnel_f82_B(F0, F82)));
|
|
}
|
|
|
|
ccl_device float ior_from_F0(const float f0)
|
|
{
|
|
const float sqrt_f0 = sqrtf(clamp(f0, 0.0f, 0.99f));
|
|
return (1.0f + sqrt_f0) / (1.0f - sqrt_f0);
|
|
}
|
|
|
|
ccl_device float F0_from_ior(const float ior)
|
|
{
|
|
return sqr((ior - 1.0f) / (ior + 1.0f));
|
|
}
|
|
|
|
ccl_device float schlick_fresnel(const float u)
|
|
{
|
|
const float m = clamp(1.0f - u, 0.0f, 1.0f);
|
|
const float m2 = m * m;
|
|
return m2 * m2 * m; // pow(m, 5)
|
|
}
|
|
|
|
/* Calculate the fresnel color, which is a blend between white and the F0 color */
|
|
ccl_device_forceinline Spectrum interpolate_fresnel_color(const float3 L,
|
|
const float3 H,
|
|
const float ior,
|
|
Spectrum F0)
|
|
{
|
|
/* Compute the real Fresnel term and remap it from real_F0..1 to F0..1.
|
|
* The reason why we use this remapping instead of directly doing the
|
|
* Schlick approximation mix(F0, 1.0, (1.0-cosLH)^5) is that for cases
|
|
* with similar IORs (e.g. ice in water), the relative IOR can be close
|
|
* enough to 1.0 that the Schlick approximation becomes inaccurate. */
|
|
const float real_F = fresnel_dielectric_cos(dot(L, H), ior);
|
|
const float real_F0 = fresnel_dielectric_cos(1.0f, ior);
|
|
|
|
return mix(F0, one_spectrum(), inverse_lerp(real_F0, 1.0f, real_F));
|
|
}
|
|
|
|
/* If the shading normal results in specular reflection in the lower hemisphere, raise the shading
|
|
* normal towards the geometry normal so that the specular reflection is just above the surface.
|
|
* Only used for glossy materials. */
|
|
ccl_device float3 ensure_valid_specular_reflection(const float3 Ng, const float3 I, float3 N)
|
|
{
|
|
const float3 R = 2 * dot(N, I) * N - I;
|
|
|
|
const float Iz = dot(I, Ng);
|
|
kernel_assert(Iz >= 0);
|
|
|
|
/* Reflection rays may always be at least as shallow as the incoming ray. */
|
|
const float threshold = min(0.9f * Iz, 0.01f);
|
|
if (dot(Ng, R) >= threshold) {
|
|
return N;
|
|
}
|
|
|
|
/* Form coordinate system with Ng as the Z axis and N inside the X-Z-plane.
|
|
* The X axis is found by normalizing the component of N that's orthogonal to Ng.
|
|
* The Y axis isn't actually needed.
|
|
*/
|
|
const float3 X = safe_normalize_fallback(N - dot(N, Ng) * Ng, N);
|
|
|
|
/* Calculate N.z and N.x in the local coordinate system.
|
|
*
|
|
* The goal of this computation is to find a N' that is rotated towards Ng just enough
|
|
* to lift R' above the threshold (here called t), therefore dot(R', Ng) = t.
|
|
*
|
|
* According to the standard reflection equation,
|
|
* this means that we want dot(2*dot(N', I)*N' - I, Ng) = t.
|
|
*
|
|
* Since the Z axis of our local coordinate system is Ng, dot(x, Ng) is just x.z, so we get
|
|
* 2*dot(N', I)*N'.z - I.z = t.
|
|
*
|
|
* The rotation is simple to express in the coordinate system we formed -
|
|
* since N lies in the X-Z-plane, we know that N' will also lie in the X-Z-plane,
|
|
* so N'.y = 0 and therefore dot(N', I) = N'.x*I.x + N'.z*I.z .
|
|
*
|
|
* Furthermore, we want N' to be normalized, so N'.x = sqrt(1 - N'.z^2).
|
|
*
|
|
* With these simplifications, we get the equation
|
|
* 2*(sqrt(1 - N'.z^2)*I.x + N'.z*I.z)*N'.z - I.z = t,
|
|
* or
|
|
* 2*sqrt(1 - N'.z^2)*I.x*N'.z = t + I.z * (1 - 2*N'.z^2),
|
|
* after rearranging terms.
|
|
* Raise both sides to the power of two and substitute terms with
|
|
* a = I.x^2 + I.z^2,
|
|
* b = 2*(a + Iz*t),
|
|
* c = (Iz + t)^2,
|
|
* we obtain
|
|
* 4*a*N'.z^4 - 2*b*N'.z^2 + c = 0.
|
|
*
|
|
* The only unknown here is N'.z, so we can solve for that.
|
|
*
|
|
* The equation has four solutions in general, two can immediately be discarded because they're
|
|
* negative so N' would lie in the lower hemisphere; one solves
|
|
* 2*sqrt(1 - N'.z^2)*I.x*N'.z = -(t + I.z * (1 - 2*N'.z^2))
|
|
* instead of the original equation (before squaring both sides).
|
|
* Therefore only one root is valid.
|
|
*/
|
|
|
|
const float Ix = dot(I, X);
|
|
|
|
const float a = sqr(Ix) + sqr(Iz);
|
|
const float b = 2.0f * (a + Iz * threshold);
|
|
const float c = sqr(threshold + Iz);
|
|
|
|
/* In order that the root formula solves 2*sqrt(1 - N'.z^2)*I.x*N'.z = t + I.z - 2*I.z*N'.z^2,
|
|
* Ix and (t + I.z * (1 - 2*N'.z^2)) must have the same sign (the rest terms are non-negative by
|
|
* definition). */
|
|
const float Nz2 = (Ix < 0) ? 0.25f * (b + safe_sqrtf(sqr(b) - 4.0f * a * c)) / a :
|
|
0.25f * (b - safe_sqrtf(sqr(b) - 4.0f * a * c)) / a;
|
|
|
|
const float Nx = safe_sqrtf(1.0f - Nz2);
|
|
const float Nz = safe_sqrtf(Nz2);
|
|
|
|
return Nx * X + Nz * Ng;
|
|
}
|
|
|
|
/* Do not call #ensure_valid_specular_reflection if the primitive type is curve or if the geometry
|
|
* normal and the shading normal is the same. */
|
|
ccl_device float3 maybe_ensure_valid_specular_reflection(ccl_private ShaderData *sd,
|
|
const float3 N)
|
|
{
|
|
if ((sd->flag & SD_USE_BUMP_MAP_CORRECTION) == 0) {
|
|
return N;
|
|
}
|
|
if ((sd->type & PRIMITIVE_CURVE) || isequal(sd->Ng, N)) {
|
|
return N;
|
|
}
|
|
return ensure_valid_specular_reflection(sd->Ng, sd->wi, N);
|
|
}
|
|
|
|
/* Principled Hair albedo and absorption coefficients. */
|
|
ccl_device_inline float bsdf_principled_hair_albedo_roughness_scale(
|
|
const float azimuthal_roughness)
|
|
{
|
|
const float x = azimuthal_roughness;
|
|
return (((((0.245f * x) + 5.574f) * x - 10.73f) * x + 2.532f) * x - 0.215f) * x + 5.969f;
|
|
}
|
|
|
|
ccl_device_inline Spectrum
|
|
bsdf_principled_hair_sigma_from_reflectance(const Spectrum color, const float azimuthal_roughness)
|
|
{
|
|
const Spectrum sigma = log(max(color, zero_spectrum())) /
|
|
bsdf_principled_hair_albedo_roughness_scale(azimuthal_roughness);
|
|
return sigma * sigma;
|
|
}
|
|
|
|
ccl_device_inline Spectrum bsdf_principled_hair_sigma_from_concentration(const float eumelanin,
|
|
const float pheomelanin)
|
|
{
|
|
const float3 eumelanin_color = make_float3(0.506f, 0.841f, 1.653f);
|
|
const float3 pheomelanin_color = make_float3(0.343f, 0.733f, 1.924f);
|
|
|
|
return eumelanin * rgb_to_spectrum(eumelanin_color) +
|
|
pheomelanin * rgb_to_spectrum(pheomelanin_color);
|
|
}
|
|
|
|
/* Computes the weight for base closure(s) which are layered under another closure.
|
|
* layer_albedo is an estimate of the top layer's reflectivity, while weight is the closure weight
|
|
* of the entire base+top combination. */
|
|
ccl_device_inline Spectrum closure_layering_weight(const Spectrum layer_albedo,
|
|
const Spectrum weight)
|
|
{
|
|
return weight * saturatef(1.0f - reduce_max(safe_divide_color(layer_albedo, weight)));
|
|
}
|
|
|
|
/* ******** Thin-film iridescence implementation ********
|
|
*
|
|
* Based on "A Practical Extension to Microfacet Theory for the Modeling of Varying Iridescence"
|
|
* by Laurent Belcour and Pascal Barla.
|
|
* https://belcour.github.io/blog/research/publication/2017/05/01/brdf-thin-film.html.
|
|
*/
|
|
|
|
/**
|
|
* Evaluate the sensitivity functions for the Fourier-space spectral integration.
|
|
*/
|
|
ccl_device_inline complex<Spectrum> iridescence_lookup_sensitivity(KernelGlobals kg,
|
|
const float OPD)
|
|
{
|
|
/* The LUT covers 0 to 60 um. */
|
|
float x = M_2PI_F * OPD / 60000.0f;
|
|
const int size = THIN_FILM_TABLE_SIZE;
|
|
|
|
const float3 re = make_float3(
|
|
lookup_table_read(kg, x, kernel_data.tables.thin_film_table + 0 * size, size),
|
|
lookup_table_read(kg, x, kernel_data.tables.thin_film_table + 1 * size, size),
|
|
lookup_table_read(kg, x, kernel_data.tables.thin_film_table + 2 * size, size));
|
|
const float3 im = make_float3(
|
|
lookup_table_read(kg, x, kernel_data.tables.thin_film_table + 3 * size, size),
|
|
lookup_table_read(kg, x, kernel_data.tables.thin_film_table + 4 * size, size),
|
|
lookup_table_read(kg, x, kernel_data.tables.thin_film_table + 5 * size, size));
|
|
|
|
return {re, im};
|
|
}
|
|
|
|
template<typename SpectrumOrFloat>
|
|
ccl_device_inline float3 iridescence_airy_summation(KernelGlobals kg,
|
|
const float R12,
|
|
const SpectrumOrFloat R23,
|
|
const float OPD,
|
|
const complex<SpectrumOrFloat> phasor)
|
|
{
|
|
const float T121 = 1.0f - R12;
|
|
|
|
const SpectrumOrFloat R123 = R12 * R23;
|
|
const SpectrumOrFloat r123 = sqrt(R123);
|
|
const SpectrumOrFloat Rs = sqr(T121) * R23 / (1.0f - R123);
|
|
|
|
/* Initialize complex number for exp(i * phi)^m, equivalent to {cos(m * phi), sin(m * phi)} as
|
|
* used in equation 10. */
|
|
complex<SpectrumOrFloat> accumulator = phasor;
|
|
|
|
/* Perform summation over path order differences (equation 10). */
|
|
Spectrum R = make_spectrum(Rs + R12); /* C0 */
|
|
SpectrumOrFloat Cm = (Rs - T121);
|
|
|
|
/* Truncate after m=3, higher differences have barely any impact. */
|
|
for (int m = 1; m < 4; m++) {
|
|
Cm *= r123;
|
|
const complex<Spectrum> S = iridescence_lookup_sensitivity(kg, m * OPD);
|
|
R += Cm * 2.0f * (accumulator.re * S.re + accumulator.im * S.im);
|
|
accumulator *= phasor;
|
|
}
|
|
return R;
|
|
}
|
|
|
|
/* Template meta-programming helper to be able to have an if-constexpr expression
|
|
* to switch between conductive (for Spectrum) or dielectric (for float) Fresnel.
|
|
* Essentially std::is_same<T, Spectrum>, but also works on GPU. */
|
|
template<class T> struct fresnel_info;
|
|
template<> struct fresnel_info<float> {
|
|
ccl_static_constexpr bool conductive = false;
|
|
};
|
|
template<> struct fresnel_info<Spectrum> {
|
|
ccl_static_constexpr bool conductive = true;
|
|
};
|
|
|
|
template<typename SpectrumOrFloat>
|
|
ccl_device Spectrum fresnel_iridescence(KernelGlobals kg,
|
|
const float ambient_ior,
|
|
const FresnelThinFilm thin_film,
|
|
const complex<SpectrumOrFloat> substrate_ior,
|
|
ccl_private const Spectrum *F82,
|
|
const float cos_theta_1,
|
|
ccl_private float *r_cos_theta_3)
|
|
{
|
|
/* For films below 1nm, the wave-optic-based Airy summation approach no longer applies,
|
|
* so blend towards the case without coating. */
|
|
float film_ior = thin_film.ior;
|
|
if (thin_film.thickness < 1.0f) {
|
|
film_ior = mix(ambient_ior, film_ior, smoothstep(0.0f, 1.0f, thin_film.thickness));
|
|
}
|
|
|
|
float cos_theta_2;
|
|
/* The real component of exp(i * phi12), equivalent to cos(phi12). */
|
|
float2 phasor12_real;
|
|
|
|
/* Compute reflection at the top interface (ambient to film). */
|
|
const float2 R12 = fresnel_dielectric_polarized(
|
|
cos_theta_1, film_ior / ambient_ior, &cos_theta_2, &phasor12_real);
|
|
if (isequal(R12, one_float2())) {
|
|
/* TIR at the top interface. */
|
|
return one_spectrum();
|
|
}
|
|
|
|
/* Compute reflection at the bottom interface (film to substrate). */
|
|
SpectrumOrFloat R23_s, R23_p;
|
|
complex<SpectrumOrFloat> phasor23_s, phasor23_p;
|
|
if constexpr (fresnel_info<SpectrumOrFloat>::conductive) {
|
|
/* Material is a conductor. */
|
|
if (F82 != nullptr) {
|
|
/* Calculate reflectance using the F82 model if the caller requested it. */
|
|
|
|
/* Scale n and k by the film ior, and recompute F0. */
|
|
const Spectrum n = substrate_ior.re / film_ior;
|
|
const Spectrum k_sq = sqr(substrate_ior.im / film_ior);
|
|
const Spectrum F0 = (sqr(n - 1.0f) + k_sq) / (sqr(n + 1.0f) + k_sq);
|
|
R23_s = fresnel_f82(-cos_theta_2, F0, fresnel_f82_B(F0, *F82));
|
|
R23_p = R23_s;
|
|
|
|
fresnel_conductor_polarized(
|
|
-cos_theta_2, film_ior, substrate_ior, nullptr, nullptr, &phasor23_s, &phasor23_p);
|
|
}
|
|
else {
|
|
fresnel_conductor_polarized(
|
|
-cos_theta_2, film_ior, substrate_ior, &R23_s, &R23_p, &phasor23_s, &phasor23_p);
|
|
}
|
|
}
|
|
else {
|
|
/* Material is a dielectric. */
|
|
float2 phasor23_real;
|
|
const float2 R23 = fresnel_dielectric_polarized(
|
|
-cos_theta_2, substrate_ior.re / film_ior, r_cos_theta_3, &phasor23_real);
|
|
|
|
if (isequal(R23, one_float2())) {
|
|
/* TIR at the bottom interface.
|
|
* All the Airy summation math still simplifies to 1.0 in this case. */
|
|
return one_spectrum();
|
|
}
|
|
|
|
R23_s = R23.x;
|
|
R23_p = R23.y;
|
|
phasor23_s = {phasor23_real.x, 0.0f};
|
|
phasor23_p = {phasor23_real.y, 0.0f};
|
|
}
|
|
|
|
/* Compute optical path difference inside the thin film. */
|
|
const float OPD = -2.0f * film_ior * thin_film.thickness * cos_theta_2;
|
|
|
|
/* Compute full phase shifts due to reflection, as a complex number exp(i * (phi23 + phi21)).
|
|
* This complex form avoids the atan2 and cos calls needed to directly get the phase shift. */
|
|
const complex<SpectrumOrFloat> phasor_s = phasor23_s * -phasor12_real.x;
|
|
const complex<SpectrumOrFloat> phasor_p = phasor23_p * -phasor12_real.y;
|
|
|
|
/* Perform Airy summation and average the polarizations. */
|
|
const Spectrum R_s = iridescence_airy_summation(kg, R12.x, R23_s, OPD, phasor_s);
|
|
const Spectrum R_p = iridescence_airy_summation(kg, R12.y, R23_p, OPD, phasor_p);
|
|
|
|
return saturate(mix(R_s, R_p, 0.5f));
|
|
}
|
|
|
|
CCL_NAMESPACE_END
|