Cycles: Thin film iridescence for metals

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

intern/cycles/kernel/closure/bsdf_microfacet.h

@@ -29,11 +29,6 @@ enum MicrofacetFresnel {
   F82_TINT,
 };
 
-struct FresnelThinFilm {
-  float thickness;
-  float ior;
-};
-
 struct FresnelDielectricTint {
   FresnelThinFilm thin_film;
 
@@ -42,7 +37,8 @@ struct FresnelDielectricTint {
 };
 
 struct FresnelConductor {
-  ComplexIOR<Spectrum> ior;
+  FresnelThinFilm thin_film;
+  complex<Spectrum> ior;
 };
 
 struct FresnelGeneralizedSchlick {
@@ -57,6 +53,8 @@ struct FresnelGeneralizedSchlick {
 };
 
 struct FresnelF82Tint {
+  FresnelThinFilm thin_film;
+
   /* Perpendicular reflectivity. */
   Spectrum f0;
   /* Precomputed (1-cos)^6 factor for edge tint. */
@@ -250,7 +248,15 @@ ccl_device_forceinline void microfacet_fresnel(KernelGlobals kg,
   }
   else if (bsdf->fresnel_type == MicrofacetFresnel::CONDUCTOR) {
     ccl_private FresnelConductor *fresnel = (ccl_private FresnelConductor *)bsdf->fresnel;
-    *r_reflectance = fresnel_conductor(cos_theta_i, fresnel->ior);
+
+    if (fresnel->thin_film.thickness > THINFILM_THICKNESS_CUTOFF) {
+      *r_reflectance = fresnel_iridescence<Spectrum>(
+          kg, 1.0f, fresnel->thin_film, fresnel->ior, nullptr, cos_theta_i, r_cos_theta_t);
+    }
+    else {
+      *r_reflectance = fresnel_conductor(cos_theta_i, fresnel->ior);
+    }
+
     *r_transmittance = zero_spectrum();
   }
   else if (bsdf->fresnel_type == MicrofacetFresnel::F82_TINT) {
@@ -258,7 +264,24 @@ ccl_device_forceinline void microfacet_fresnel(KernelGlobals kg,
      * Essentially, this is the usual Schlick Fresnel with an additional cosI*(1-cosI)^6
      * term which modulates the reflectivity around acos(1/7) degrees (ca. 82°). */
     ccl_private FresnelF82Tint *fresnel = (ccl_private FresnelF82Tint *)bsdf->fresnel;
-    *r_reflectance = fresnel_f82(cos_theta_i, fresnel->f0, fresnel->b);
+
+    if (fresnel->thin_film.thickness > THINFILM_THICKNESS_CUTOFF) {
+      /* Estimate n and k by reinterpreting F0 and F82 as r and g from "Artist Friendly Metallic
+       * Fresnel" by Ole Gulbrandsen. */
+      const Spectrum r = min(fresnel->f0, make_float3(0.999f));
+      const Spectrum g = fresnel_f82(1.0f / 7.0f, fresnel->f0, fresnel->b);
+
+      const Spectrum sqrt_r = sqrt(r);
+      const Spectrum n = mix((1.0f + sqrt_r) / (1.0f - sqrt_r), (1.0f - r) / (1.0f + r), g);
+      const Spectrum k = safe_sqrt((r * sqr(n + 1) - sqr(n - 1)) / (1.0f - r));
+
+      *r_reflectance = fresnel_iridescence<Spectrum>(
+          kg, 1.0f, fresnel->thin_film, {n, k}, &g, cos_theta_i, r_cos_theta_t);
+    }
+    else {
+      *r_reflectance = fresnel_f82(cos_theta_i, fresnel->f0, fresnel->b);
+    }
+
     *r_transmittance = zero_spectrum();
   }
   else if (bsdf->fresnel_type == MicrofacetFresnel::GENERALIZED_SCHLICK) {
@@ -270,13 +293,8 @@ ccl_device_forceinline void microfacet_fresnel(KernelGlobals kg,
        * Principled BSDF for now, so it's fine to not support custom exponents and F90. */
       kernel_assert(fresnel->exponent < 0.0f);
       kernel_assert(fresnel->f90 == one_spectrum());
-      F = fresnel_iridescence(kg,
-                              1.0f,
-                              fresnel->thin_film.ior,
-                              bsdf->ior,
-                              cos_theta_i,
-                              fresnel->thin_film.thickness,
-                              r_cos_theta_t);
+      F = fresnel_iridescence<float>(
+          kg, 1.0f, fresnel->thin_film, {bsdf->ior, 0.0f}, nullptr, cos_theta_i, r_cos_theta_t);
       /* Apply F0 scaling (here per-channel, since iridescence produces colored output).
        * Note that the usual approach (as used below) cannot be used here, since F may be below
        * F0_real. Therefore, use a different approach: Scale the result by (F0 / F0_real), with
@@ -446,11 +464,18 @@ ccl_device Spectrum bsdf_microfacet_estimate_albedo(KernelGlobals kg,
   }
   else if (bsdf->fresnel_type == MicrofacetFresnel::F82_TINT) {
     ccl_private FresnelF82Tint *fresnel = (ccl_private FresnelF82Tint *)bsdf->fresnel;
-    const float rough = sqrtf(sqrtf(bsdf->alpha_x * bsdf->alpha_y));
-    const float s = lookup_table_read_3D(
-        kg, rough, cos_NI, 0.5f, kernel_data.tables.ggx_gen_schlick_s, 16, 16, 16);
-    /* TODO: Precompute B factor term and account for it here. */
-    reflectance = mix(fresnel->f0, one_spectrum(), s);
+
+    if (fresnel->thin_film.thickness > THINFILM_THICKNESS_CUTOFF) {
+      /* Precomputing LUTs for thin-film iridescence isn't viable, so fall back to the specular
+       * reflection approximation from the microfacet_fresnel call above in that case. */
+    }
+    else {
+      const float rough = sqrtf(sqrtf(bsdf->alpha_x * bsdf->alpha_y));
+      const float s = lookup_table_read_3D(
+          kg, rough, cos_NI, 0.5f, kernel_data.tables.ggx_gen_schlick_s, 16, 16, 16);
+      /* TODO: Precompute B factor term and account for it here. */
+      reflectance = mix(fresnel->f0, one_spectrum(), s);
+    }
   }
   else if ((bsdf->fresnel_type == MicrofacetFresnel::DIELECTRIC ||
             bsdf->fresnel_type == MicrofacetFresnel::DIELECTRIC_TINT) &&

intern/cycles/kernel/closure/bsdf_util.h

@@ -16,19 +16,37 @@
 
 CCL_NAMESPACE_BEGIN
 
-template<typename T> struct ComplexIOR {
-  T eta;
-  T k;
+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_phi is set to the phase shift on reflection.
+ * 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_phi)
+                                               ccl_private float2 *r_cos_phi)
 {
   kernel_assert(!isnan_safe(cos_theta_i));
 
@@ -37,7 +55,7 @@ ccl_device float2 fresnel_dielectric_polarized(float cos_theta_i,
   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_phi) {
+    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.
@@ -46,7 +64,7 @@ ccl_device float2 fresnel_dielectric_polarized(float cos_theta_i,
        * `r_phi->x = -2.0f * atanf(fac / cosThetaI);`
        * `r_phi->y = -2.0f * atanf(fac / (cosThetaI * sqr(eta)));`
        */
-      *r_phi = zero_float2();
+      *r_cos_phi = one_float2();
     }
     return one_float2();
   }
@@ -63,8 +81,8 @@ ccl_device float2 fresnel_dielectric_polarized(float cos_theta_i,
   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_phi) {
-    *r_phi = make_float2(r_s < 0.0f, r_p < 0.0f) * M_PI_F;
+  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. */
@@ -117,23 +135,23 @@ ccl_device_inline float fresnel_dielectric_Fss(const float eta)
   return (eta - 1.0f) / (4.08567f + 1.00071f * eta);
 }
 
-/* Evaluates the Fresnel equations at a dielectric-conductor interface. If requested by the caller,
- * sets r_R_s and r_R_p to the reflectances for perpendicular and parallel polarized light, and
- * sets r_phi_s and r_phi_p to the phase shifts due to reflection.
+/* 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 ComplexIOR<Spectrum> conductor_ior,
+                                            const complex<Spectrum> conductor_ior,
                                             ccl_private Spectrum *r_R_s,
                                             ccl_private Spectrum *r_R_p,
-                                            ccl_private Spectrum *r_phi_s,
-                                            ccl_private Spectrum *r_phi_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.eta;
-  const Spectrum k2 = conductor_ior.k;
+  const Spectrum eta2 = conductor_ior.re;
+  const Spectrum k2 = conductor_ior.im;
 
   const float eta1_sq = sqr(eta1);
   const Spectrum eta2_sq = sqr(eta2);
@@ -143,35 +161,37 @@ ccl_device void fresnel_conductor_polarized(const float cosi,
   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_float3());
-  const Spectrum v_sq = max(0.5f * (t2 - t1), zero_float3());
+  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 = (sqr(eta1 * cosi - u) + v_sq) / (sqr(eta1 * cosi + u) + v_sq);
+    *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;
-    const Spectrum R_p = (sqr(t3 - eta1 * u) + sqr(t4 - eta1 * v)) /
-                         (sqr(t3 + eta1 * u) + sqr(t4 + eta1 * v));
-    const auto mask = isequal_mask(eta2, zero_spectrum()) & isequal_mask(k2, zero_spectrum());
-    *r_R_p = select(mask, one_spectrum(), R_p);
+    *r_R_p = safe_divide(sqr(t3 - eta1 * u) + sqr(t4 - eta1 * v),
+                         sqr(t3 + eta1 * u) + sqr(t4 + eta1 * v));
   }
 
-  if (r_phi_s && r_phi_p) {
-    const Spectrum s_numerator = 2.0f * eta1 * cosi * v;
-    const Spectrum s_denominator = u_sq + v_sq - sqr(eta1 * cosi);
-    *r_phi_s = atan2(-s_numerator, -s_denominator);
+  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 p_numerator = 2.0f * eta1 * cosi * (two_eta2_k2 * u - (eta2_sq - k2_sq) * v);
-    const Spectrum p_denominator = sqr((eta2_sq + k2_sq) * cosi) - eta1_sq * (u_sq + v_sq);
-    *r_phi_p = atan2(p_numerator, p_denominator);
+    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 ComplexIOR<Spectrum> 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);
@@ -218,7 +238,7 @@ ccl_device_inline Spectrum fresnel_f82(const float cosi, const Spectrum F0, cons
 }
 
 /* Approximates the average single-scattering Fresnel for a physical conductor. */
-ccl_device_inline Spectrum fresnel_conductor_Fss(const ComplexIOR<Spectrum> ior)
+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. */
@@ -398,108 +418,150 @@ ccl_device_inline Spectrum closure_layering_weight(const Spectrum layer_albedo,
 
 /**
  * Evaluate the sensitivity functions for the Fourier-space spectral integration.
- * The code here uses the Gaussian fit for the CIE XYZ curves that is provided
- * in the reference implementation.
- * For details on what this actually represents, see the paper.
- * In theory we should pre-compute the sensitivity functions for the working RGB
- * color-space, remap them to be functions of (light) frequency, take their Fourier
- * transform and store them as a LUT that gets looked up here.
- * In practice, using the XYZ fit and converting the result from XYZ to RGB is easier.
  */
-ccl_device_inline Spectrum iridescence_lookup_sensitivity(KernelGlobals kg,
-                                                          const float OPD,
-                                                          const float shift)
+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 mag = make_float3(
+  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 phase = make_float3(
+  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 mag * cos(phase - shift);
+  return {re, im};
 }
 
+template<typename SpectrumOrFloat>
 ccl_device_inline float3 iridescence_airy_summation(KernelGlobals kg,
-                                                    const float T121,
                                                     const float R12,
-                                                    const float R23,
+                                                    const SpectrumOrFloat R23,
                                                     const float OPD,
-                                                    const float phi)
+                                                    const complex<SpectrumOrFloat> phasor)
 {
-  if (R23 == 1.0f) {
-    /* Shortcut for TIR on the bottom interface. */
-    return one_float3();
-  }
+  const float T121 = 1.0f - R12;
 
-  const float R123 = R12 * R23;
-  const float r123 = sqrtf(R123);
-  const float Rs = sqr(T121) * R23 / (1.0f - R123);
+  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). */
-  float3 R = make_float3(R12 + Rs); /* C0 */
-  float Cm = (Rs - T121);
+  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;
-    R += Cm * 2.0f * iridescence_lookup_sensitivity(kg, m * OPD, m * phi);
+    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,
-                                        float eta1,
-                                        float eta2,
-                                        float eta3,
-                                        float cos_theta_1,
-                                        const float thickness,
+                                        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 30nm, the wave-optic-based Airy summation approach no longer applies,
+  /* For films below 1nm, the wave-optic-based Airy summation approach no longer applies,
    * so blend towards the case without coating. */
-  if (thickness < 30.0f) {
-    eta2 = mix(eta1, eta2, smoothstep(0.0f, 30.0f, thickness));
+  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;
-  float2 phi12;
-  float2 phi23;
+  /* 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, eta2 / eta1, &cos_theta_2, &phi12);
+  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 optical path difference inside the thin film. */
-  const float OPD = -2.0f * eta2 * thickness * cos_theta_2;
+  /* 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. */
 
-  /* Compute reflection at the bottom interface (film to medium). */
-  const float2 R23 = fresnel_dielectric_polarized(
-      -cos_theta_2, eta3 / eta2, r_cos_theta_3, &phi23);
-  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();
+      /* 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 helper parameters. */
-  const float2 T121 = one_float2() - R12;
-  const float2 phi = make_float2(M_PI_F, M_PI_F) - phi12 + phi23;
+  /* 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. */
-  float3 R = mix(iridescence_airy_summation(kg, T121.x, R12.x, R23.x, OPD, phi.x),
-                 iridescence_airy_summation(kg, T121.y, R12.y, R23.y, OPD, phi.y),
-                 0.5f);
+  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(R);
+  return saturate(mix(R_s, R_p, 0.5f));
 }
 
 CCL_NAMESPACE_END

intern/cycles/kernel/osl/closures_setup.h

@@ -381,6 +381,9 @@ ccl_device void osl_closure_conductor_bsdf_setup(KernelGlobals kg,
     preserve_energy = (closure->distribution == make_string("multi_ggx", 16842698693386468366ull));
   }
 
+  fresnel->thin_film.thickness = closure->thinfilm_thickness;
+  fresnel->thin_film.ior = closure->thinfilm_ior;
+
   fresnel->ior = {rgb_to_spectrum(closure->ior), rgb_to_spectrum(closure->extinction)};
   bsdf_microfacet_setup_fresnel_conductor(kg, bsdf, sd, fresnel, preserve_energy);
 }
@@ -607,6 +610,9 @@ ccl_device void osl_closure_microfacet_f82_tint_setup(
   }
 
   fresnel->f0 = rgb_to_spectrum(closure->f0);
+  fresnel->thin_film.thickness = closure->thinfilm_thickness;
+  fresnel->thin_film.ior = closure->thinfilm_ior;
+
   bsdf_microfacet_setup_fresnel_f82_tint(
       kg, bsdf, sd, fresnel, rgb_to_spectrum(closure->f82), preserve_energy);
 }

intern/cycles/kernel/osl/closures_template.h

@@ -79,6 +79,8 @@ OSL_CLOSURE_STRUCT_BEGIN(ConductorBSDF, conductor_bsdf)
   OSL_CLOSURE_STRUCT_MEMBER(ConductorBSDF, VECTOR, packed_float3, ior, nullptr)
   OSL_CLOSURE_STRUCT_MEMBER(ConductorBSDF, VECTOR, packed_float3, extinction, nullptr)
   OSL_CLOSURE_STRUCT_MEMBER(ConductorBSDF, STRING, DeviceString, distribution, nullptr)
+  OSL_CLOSURE_STRUCT_MEMBER(ConductorBSDF, FLOAT, float, thinfilm_thickness, "thinfilm_thickness")
+  OSL_CLOSURE_STRUCT_MEMBER(ConductorBSDF, FLOAT, float, thinfilm_ior, "thinfilm_ior")
 OSL_CLOSURE_STRUCT_END(ConductorBSDF, conductor_bsdf)
 
 OSL_CLOSURE_STRUCT_BEGIN(GeneralizedSchlickBSDF, generalized_schlick_bsdf)
@@ -117,6 +119,9 @@ OSL_CLOSURE_STRUCT_BEGIN(MicrofacetF82Tint, microfacet_f82_tint)
   OSL_CLOSURE_STRUCT_MEMBER(MicrofacetF82Tint, FLOAT, float, alpha_y, nullptr)
   OSL_CLOSURE_STRUCT_MEMBER(MicrofacetF82Tint, VECTOR, packed_float3, f0, nullptr)
   OSL_CLOSURE_STRUCT_MEMBER(MicrofacetF82Tint, VECTOR, packed_float3, f82, nullptr)
+  OSL_CLOSURE_STRUCT_MEMBER(
+      MicrofacetF82Tint, FLOAT, float, thinfilm_thickness, "thinfilm_thickness")
+  OSL_CLOSURE_STRUCT_MEMBER(MicrofacetF82Tint, FLOAT, float, thinfilm_ior, "thinfilm_ior")
 OSL_CLOSURE_STRUCT_END(MicrofacetF82Tint, microfacet)
 
 OSL_CLOSURE_STRUCT_BEGIN(MicrofacetMultiGGXGlass, microfacet_multi_ggx_glass)

intern/cycles/kernel/osl/shaders/node_metallic_bsdf.osl

@@ -14,6 +14,8 @@ shader node_metallic_bsdf(color BaseColor = color(0.617, 0.577, 0.540),
                           float Roughness = 0.5,
                           float Anisotropy = 0.0,
                           float Rotation = 0.0,
+                          float ThinFilmThickness = 0.0,
+                          float ThinFilmIOR = 1.33,
                           normal Normal = N,
                           normal Tangent = 0.0,
                           output closure color BSDF = 0)
@@ -35,10 +37,29 @@ shader node_metallic_bsdf(color BaseColor = color(0.617, 0.577, 0.540),
   if (fresnel_type == "f82") {
     color F0 = clamp(BaseColor, color(0.0), color(1.0));
     color F82 = clamp(EdgeTint, color(0.0), color(1.0));
-    BSDF = microfacet_f82_tint(distribution, Normal, T, alpha_x, alpha_y, F0, F82);
+    BSDF = microfacet_f82_tint(distribution,
+                               Normal,
+                               T,
+                               alpha_x,
+                               alpha_y,
+                               F0,
+                               F82,
+                               "thinfilm_thickness",
+                               ThinFilmThickness,
+                               "thinfilm_ior",
+                               ThinFilmIOR);
   }
   else {
-    BSDF = conductor_bsdf(
-        Normal, T, alpha_x, alpha_y, max(IOR, 0.0), max(Extinction, 0.0), distribution);
+    BSDF = conductor_bsdf(Normal,
+                          T,
+                          alpha_x,
+                          alpha_y,
+                          max(IOR, 0.0),
+                          max(Extinction, 0.0),
+                          distribution,
+                          "thinfilm_thickness",
+                          ThinFilmThickness,
+                          "thinfilm_ior",
+                          ThinFilmIOR);
   }
 }

intern/cycles/kernel/osl/shaders/node_principled_bsdf.osl

@@ -155,7 +155,17 @@ shader node_principled_bsdf(string distribution = "multi_ggx",
   if (metallic > CLOSURE_WEIGHT_CUTOFF) {
     color F0 = clamped_base_color;
     color F82 = min(specular_tint, color(1.0));
-    MetallicBSDF = microfacet_f82_tint(distribution, Normal, T, alpha_x, alpha_y, F0, F82);
+    MetallicBSDF = microfacet_f82_tint(distribution,
+                                       Normal,
+                                       T,
+                                       alpha_x,
+                                       alpha_y,
+                                       F0,
+                                       F82,
+                                       "thinfilm_thickness",
+                                       ThinFilmThickness,
+                                       "thinfilm_ior",
+                                       ThinFilmIOR);
     BSDF = mix(BSDF, MetallicBSDF, metallic);
   }
 

intern/cycles/kernel/svm/closure.h

@@ -354,6 +354,9 @@ ccl_device
             fresnel->f0 = rgb_to_spectrum(clamped_base_color);
             const Spectrum f82 = min(specular_tint, one_spectrum());
 
+            fresnel->thin_film.thickness = thinfilm_thickness;
+            fresnel->thin_film.ior = thinfilm_ior;
+
             /* setup bsdf */
             sd->flag |= bsdf_microfacet_ggx_setup(bsdf);
             const bool is_multiggx = (distribution == CLOSURE_BSDF_MICROFACET_MULTI_GGX_GLASS_ID);
@@ -576,7 +579,21 @@ ccl_device
         bsdf->N = valid_reflection_N;
         bsdf->ior = 1.0f;
 
-        const ClosureType distribution = (ClosureType)data_node.z;
+        uint distribution_int;
+        uint thin_film_thickness_offset;
+        uint thin_film_ior_offset;
+        uint unused;
+        svm_unpack_node_uchar4(data_node.z,
+                               &distribution_int,
+                               &thin_film_thickness_offset,
+                               &thin_film_ior_offset,
+                               &unused);
+
+        const float thin_film_thickness = fmaxf(
+            stack_load_float(stack, thin_film_thickness_offset), 1e-5f);
+        const float thin_film_ior = fmaxf(stack_load_float(stack, thin_film_ior_offset), 1e-5f);
+
+        const ClosureType distribution = (ClosureType)distribution_int;
         /* Setup BSDF */
         if (distribution == CLOSURE_BSDF_MICROFACET_BECKMANN_ID) {
           sd->flag |= bsdf_microfacet_beckmann_setup(bsdf);
@@ -591,6 +608,13 @@ ccl_device
           ccl_private FresnelConductor *fresnel = (ccl_private FresnelConductor *)
               closure_alloc_extra(sd, sizeof(FresnelConductor));
 
+          if (!fresnel) {
+            break;
+          }
+
+          fresnel->thin_film.thickness = thin_film_thickness;
+          fresnel->thin_film.ior = thin_film_ior;
+
           const float3 n = max(stack_load_float3(stack, base_ior_offset), zero_float3());
           const float3 k = max(stack_load_float3(stack, edge_tint_k_offset), zero_float3());
 
@@ -601,6 +625,13 @@ ccl_device
           ccl_private FresnelF82Tint *fresnel = (ccl_private FresnelF82Tint *)closure_alloc_extra(
               sd, sizeof(FresnelF82Tint));
 
+          if (!fresnel) {
+            break;
+          }
+
+          fresnel->thin_film.thickness = thin_film_thickness;
+          fresnel->thin_film.ior = thin_film_ior;
+
           const float3 color = saturate(stack_load_float3(stack, base_ior_offset));
           const float3 tint = saturate(stack_load_float3(stack, edge_tint_k_offset));
 

intern/cycles/scene/shader.cpp

@@ -945,15 +945,13 @@ void ShaderManager::compute_thin_film_table(const Transform &xyz_to_rgb)
    * the XYZ-to-RGB matrix to get the RGB LUT.
    *
    * That's what this function does: We load the precomputed values, convert to RGB, normalize
-   * the result to make the DC term equal to 1, convert from real/imaginary to magnitude/phase
-   * since that form is smoother and therefore interpolates more nicely, and then store that
-   * into the final table that's used by the kernel.
+   * the result to make the DC term equal to 1, and then store that into the final table that's
+   * used by the kernel.
    */
   assert(sizeof(table_thin_film_cmf) == 6 * THIN_FILM_TABLE_SIZE * sizeof(float));
   thin_film_table.resize(6 * THIN_FILM_TABLE_SIZE);
 
   float3 normalization;
-  float3 prevPhase = zero_float3();
   for (int i = 0; i < THIN_FILM_TABLE_SIZE; i++) {
     const float *table_row = table_thin_film_cmf[i];
     /* Load precomputed resampled Fourier-transformed XYZ CMFs. */
@@ -971,21 +969,13 @@ void ShaderManager::compute_thin_film_table(const Transform &xyz_to_rgb)
       normalization = 1.0f / rgbReal;
     }
 
-    /* Convert the complex value into magnitude/phase representation. */
-    const float3 rgbMag = sqrt(sqr(rgbReal) + sqr(rgbImag));
-    float3 rgbPhase = atan2(rgbImag, rgbReal);
-
-    /* Unwrap phase to avoid jumps. */
-    rgbPhase -= M_2PI_F * round((rgbPhase - prevPhase) * M_1_2PI_F);
-    prevPhase = rgbPhase;
-
     /* Store in lookup table. */
-    thin_film_table[i + 0 * THIN_FILM_TABLE_SIZE] = rgbMag.x * normalization.x;
-    thin_film_table[i + 1 * THIN_FILM_TABLE_SIZE] = rgbMag.y * normalization.y;
-    thin_film_table[i + 2 * THIN_FILM_TABLE_SIZE] = rgbMag.z * normalization.z;
-    thin_film_table[i + 3 * THIN_FILM_TABLE_SIZE] = rgbPhase.x;
-    thin_film_table[i + 4 * THIN_FILM_TABLE_SIZE] = rgbPhase.y;
-    thin_film_table[i + 5 * THIN_FILM_TABLE_SIZE] = rgbPhase.z;
+    thin_film_table[i + 0 * THIN_FILM_TABLE_SIZE] = rgbReal.x * normalization.x;
+    thin_film_table[i + 1 * THIN_FILM_TABLE_SIZE] = rgbReal.y * normalization.y;
+    thin_film_table[i + 2 * THIN_FILM_TABLE_SIZE] = rgbReal.z * normalization.z;
+    thin_film_table[i + 3 * THIN_FILM_TABLE_SIZE] = rgbImag.x * normalization.x;
+    thin_film_table[i + 4 * THIN_FILM_TABLE_SIZE] = rgbImag.y * normalization.y;
+    thin_film_table[i + 5 * THIN_FILM_TABLE_SIZE] = rgbImag.z * normalization.z;
   }
 }
 

intern/cycles/scene/shader_nodes.cpp

@@ -2339,6 +2339,9 @@ NODE_DEFINE(MetallicBsdfNode)
   SOCKET_IN_FLOAT(anisotropy, "Anisotropy", 0.0f);
   SOCKET_IN_FLOAT(rotation, "Rotation", 0.0f);
 
+  SOCKET_IN_FLOAT(thin_film_thickness, "Thin Film Thickness", 0.0f);
+  SOCKET_IN_FLOAT(thin_film_ior, "Thin Film IOR", 1.33f);
+
   SOCKET_OUT_CLOSURE(BSDF, "BSDF");
 
   return type;
@@ -2386,6 +2389,9 @@ void MetallicBsdfNode::compile(SVMCompiler &compiler)
                                      compiler.stack_assign(input("Extinction")) :
                                      compiler.stack_assign(input("Edge Tint"));
 
+  const int thin_film_thickness_offset = compiler.stack_assign(input("Thin Film Thickness"));
+  const int thin_film_ior_offset = compiler.stack_assign(input("Thin Film IOR"));
+
   ShaderInput *roughness_in = input("Roughness");
   ShaderInput *anisotropy_in = input("Anisotropy");
 
@@ -2404,7 +2410,7 @@ void MetallicBsdfNode::compile(SVMCompiler &compiler)
       normal_offset,
       compiler.encode_uchar4(
           base_color_ior_offset, edge_tint_k_offset, rotation_offset, tangent_offset),
-      distribution);
+      compiler.encode_uchar4(distribution, thin_film_thickness_offset, thin_film_ior_offset));
 }
 
 void MetallicBsdfNode::compile(OSLCompiler &compiler)

intern/cycles/scene/shader_nodes.h

@@ -622,6 +622,8 @@ class MetallicBsdfNode : public BsdfNode {
   NODE_SOCKET_API(float, roughness)
   NODE_SOCKET_API(float, anisotropy)
   NODE_SOCKET_API(float, rotation)
+  NODE_SOCKET_API(float, thin_film_thickness)
+  NODE_SOCKET_API(float, thin_film_ior)
   NODE_SOCKET_API(ClosureType, distribution)
   NODE_SOCKET_API(ClosureType, fresnel_type)
 

intern/cycles/util/math_float3.h

@@ -515,6 +515,11 @@ ccl_device_inline float3 faceforward(const float3 vector,
 }
 #endif
 
+ccl_device_inline float3 safe_sqrt(const float3 a)
+{
+  return sqrt(max(a, zero_float3()));
+}
+
 ccl_device_inline float3 project(const float3 v, const float3 v_proj)
 {
   const float len_squared = dot(v_proj, v_proj);

intern/cycles/util/types_float3.h

@@ -117,6 +117,11 @@ ccl_device_inline float3 make_float3(const int3 i)
 #endif
 }
 
+ccl_device_inline float3 make_float3(const float3 a)
+{
+  return a;
+}
+
 ccl_device_inline void print_float3(const ccl_private char *label, const float3 a)
 {
 #ifdef __KERNEL_PRINTF__

source/blender/gpu/shaders/material/gpu_shader_material_metallic.glsl

@@ -30,6 +30,8 @@ void node_bsdf_metallic(float4 base_color,
                         float3 N,
                         float3 T,
                         float weight,
+                        float thin_film_thickness,
+                        float thin_film_ior,
                         const float do_multiscatter,
                         const float use_complex_ior,
                         out Closure result)

source/blender/nodes/shader/nodes/node_shader_bsdf_metallic.cc

@@ -11,6 +11,11 @@ namespace blender::nodes::node_shader_bsdf_metallic_cc {
 
 static void node_declare(NodeDeclarationBuilder &b)
 {
+  b.use_custom_socket_order();
+
+  b.add_output<decl::Shader>("BSDF");
+  b.add_default_layout();
+
   b.add_input<decl::Color>("Base Color")
       .default_value({0.617f, 0.577f, 0.540f, 1.0f})
       .description("Color of the material");
@@ -36,7 +41,6 @@ static void node_declare(NodeDeclarationBuilder &b)
       .description(
           "Microfacet roughness of the surface (0.0 is a perfect mirror reflection, 1.0 is "
           "completely rough)");
-  ;
   b.add_input<decl::Float>("Anisotropy")
       .default_value(0.0f)
       .min(0.0f)
@@ -54,7 +58,19 @@ static void node_declare(NodeDeclarationBuilder &b)
   b.add_input<decl::Vector>("Normal").hide_value();
   b.add_input<decl::Vector>("Tangent").hide_value();
   b.add_input<decl::Float>("Weight").available(false);
-  b.add_output<decl::Shader>("BSDF");
+
+  PanelDeclarationBuilder &film = b.add_panel("Thin Film").default_closed(true);
+  film.add_input<decl::Float>("Thin Film Thickness")
+      .default_value(0.0)
+      .min(0.0f)
+      .max(100000.0f)
+      .subtype(PROP_WAVELENGTH)
+      .description("Thickness of the film in nanometers");
+  film.add_input<decl::Float>("Thin Film IOR")
+      .default_value(1.33f)
+      .min(1.0f)
+      .max(1000.0f)
+      .description("Index of refraction (IOR) of the thin film");
 }
 
 static void node_shader_buts_metallic(uiLayout *layout, bContext * /*C*/, PointerRNA *ptr)
@@ -123,6 +139,8 @@ NODE_SHADER_MATERIALX_BEGIN
   NodeItem anisotropy = get_input_value("Anisotropy", NodeItem::Type::Color3);
   NodeItem normal = get_input_link("Normal", NodeItem::Type::Vector3);
   NodeItem tangent = get_input_link("Tangent", NodeItem::Type::Vector3);
+  NodeItem thin_film_thickness = get_input_value("Thin Film Thickness", NodeItem::Type::Float);
+  NodeItem thin_film_ior = get_input_value("Thin Film IOR", NodeItem::Type::Float);
 
   NodeItem ior_out, extinction_out;
   if (node_->custom2 == SHD_PHYSICAL_CONDUCTOR) {
@@ -143,7 +161,9 @@ NODE_SHADER_MATERIALX_BEGIN
                       {"tangent", tangent},
                       {"ior", ior_out},
                       {"extinction", extinction_out},
-                      {"roughness", roughness}});
+                      {"roughness", roughness},
+                      {"thinfilm_thickness", thin_film_thickness},
+                      {"thinfilm_ior", thin_film_ior}});
 }
 #endif
 NODE_SHADER_MATERIALX_END

tests/python/storm_render_tests.py

@@ -34,6 +34,7 @@ BLOCKLIST_USD = [
 BLOCKLIST_METAL = [
     # Thinfilm
     "principled.*thinfilm.*.blend",
+    "metallic.*thinfilm.*.blend",
     # Transparency
     "transparent.blend",
     "transparent_shadow.blend",