Cycles: Add polarized Fresnel function for conductors

This PR adds a new `fresnel_conductor_polarized` function, which calculates reflectance and phase shift (if requested) for both parallel and perpendicular polarized light. This is needed for applying thin film iridescence to conductors (see !141131).

For consistency, this PR also makes `fresnel_conductor` call `fresnel_conductor_polarized` instead of using a fast approximation of the Fresnel equations that is inaccurate at lower n and k values. This will change the output of some Metallic BSDF renders using Physical Conductor and prevent discrepancies when enabling thin film iridescence.

I didn't do any rigorous performance testing, but from timing the functions outside of Blender, `fresnel_conductor_polarized` is significantly slower than the approximation, between 1.5-3x depending on the compiler. This makes sense because it has three square roots and the approximation has none. In some informal tests with metallic_multiggx_physical.blend modified to have more spheres, the new renders took around 1-2% longer on both CPU and GPU.

There are some avoidable inefficiencies in this approach of just calling `fresnel_conductor_polarized`:

- one of the three square roots could be saved since `fresnel_conductor` never needs the phase shift and there are simplifications possible when only calculating the reflectance
- there are several unnecessary multiplications by 1.0 since `fresnel_conductor` uses relative IOR and `fresnel_conductor_polarized` doesn't, though those could get optimized out if inlined

Pull Request: https://projects.blender.org/blender/blender/pulls/143903

intern/cycles/kernel/closure/bsdf_microfacet.h

@@ -42,7 +42,7 @@ struct FresnelDielectricTint {
 };
 
 struct FresnelConductor {
-  Spectrum n, k;
+  ComplexIOR<Spectrum> ior;
 };
 
 struct FresnelGeneralizedSchlick {
@@ -250,7 +250,7 @@ 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->n, fresnel->k);
+    *r_reflectance = fresnel_conductor(cos_theta_i, fresnel->ior);
     *r_transmittance = zero_spectrum();
   }
   else if (bsdf->fresnel_type == MicrofacetFresnel::F82_TINT) {
@@ -808,7 +808,7 @@ ccl_device void bsdf_microfacet_setup_fresnel_conductor(KernelGlobals kg,
   bsdf->sample_weight *= average(bsdf_microfacet_estimate_albedo(kg, sd, bsdf, true, true));
 
   if (preserve_energy) {
-    microfacet_ggx_preserve_energy(kg, bsdf, sd, fresnel_conductor_Fss(fresnel->n, fresnel->k));
+    microfacet_ggx_preserve_energy(kg, bsdf, sd, fresnel_conductor_Fss(fresnel->ior));
   }
 }
 

intern/cycles/kernel/closure/bsdf_util.h

@@ -16,6 +16,11 @@
 
 CCL_NAMESPACE_BEGIN
 
+template<typename T> struct ComplexIOR {
+  T eta;
+  T k;
+};
+
 /* 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.
  * Also returns the dot product of the refracted ray and the normal as `cos_theta_t`, as it is
@@ -112,16 +117,65 @@ ccl_device_inline float fresnel_dielectric_Fss(const float eta)
   return (eta - 1.0f) / (4.08567f + 1.00071f * eta);
 }
 
-ccl_device Spectrum fresnel_conductor(const float cosi, const Spectrum eta, const Spectrum k)
+/* 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.
+ * 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,
+                                            ccl_private Spectrum *r_R_s,
+                                            ccl_private Spectrum *r_R_p,
+                                            ccl_private Spectrum *r_phi_s,
+                                            ccl_private Spectrum *r_phi_p)
 {
-  const Spectrum cosi2 = make_spectrum(cosi * cosi);
-  const Spectrum one = make_spectrum(1.0f);
-  const Spectrum tmp_f = eta * eta + k * k;
-  const Spectrum tmp = tmp_f * cosi2;
-  const Spectrum Rparl2 = (tmp - (2.0f * eta * cosi) + one) / (tmp + (2.0f * eta * cosi) + one);
-  const Spectrum Rperp2 = (tmp_f - (2.0f * eta * cosi) + cosi2) /
-                          (tmp_f + (2.0f * eta * cosi) + cosi2);
-  return (Rparl2 + Rperp2) * 0.5f;
+  const float eta1 = ambient_ior;
+  const Spectrum eta2 = conductor_ior.eta;
+  const Spectrum k2 = conductor_ior.k;
+
+  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_float3());
+  const Spectrum v_sq = max(0.5f * (t2 - t1), zero_float3());
+  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);
+
+    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);
+  }
+
+  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);
+
+    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);
+  }
+}
+
+/* Calculates Fresnel reflectance at a dielectric-conductor interface given the relative IOR. */
+ccl_device Spectrum fresnel_conductor(const float cosi, const ComplexIOR<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. */
@@ -164,12 +218,12 @@ 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 Spectrum eta, const Spectrum k)
+ccl_device_inline Spectrum fresnel_conductor_Fss(const ComplexIOR<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, eta, k);
-  const Spectrum F82 = fresnel_conductor(1.0f / 7.0f, eta, k);
+  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)));
 }
 

intern/cycles/kernel/osl/closures_setup.h

@@ -381,8 +381,7 @@ ccl_device void osl_closure_conductor_bsdf_setup(KernelGlobals kg,
     preserve_energy = (closure->distribution == make_string("multi_ggx", 16842698693386468366ull));
   }
 
-  fresnel->n = rgb_to_spectrum(closure->ior);
-  fresnel->k = rgb_to_spectrum(closure->extinction);
+  fresnel->ior = {rgb_to_spectrum(closure->ior), rgb_to_spectrum(closure->extinction)};
   bsdf_microfacet_setup_fresnel_conductor(kg, bsdf, sd, fresnel, preserve_energy);
 }
 

intern/cycles/kernel/svm/closure.h

@@ -573,8 +573,7 @@ ccl_device
           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());
 
-          fresnel->n = rgb_to_spectrum(n);
-          fresnel->k = rgb_to_spectrum(k);
+          fresnel->ior = {rgb_to_spectrum(n), rgb_to_spectrum(k)};
           bsdf_microfacet_setup_fresnel_conductor(kg, bsdf, sd, fresnel, is_multiggx);
         }
         else {

intern/cycles/util/math_float3.h

@@ -117,9 +117,14 @@ ccl_device_inline float3 operator+(const float3 a, const float3 b)
 #  endif
 }
 
-ccl_device_inline float3 operator+(const float3 a, const float f)
+ccl_device_inline float3 operator+(const float3 a, const float b)
 {
-  return a + make_float3(f);
+  return a + make_float3(b);
+}
+
+ccl_device_inline float3 operator+(const float a, const float3 b)
+{
+  return make_float3(a) + b;
 }
 
 ccl_device_inline float3 operator-(const float3 a, const float3 b)
@@ -131,9 +136,14 @@ ccl_device_inline float3 operator-(const float3 a, const float3 b)
 #  endif
 }
 
-ccl_device_inline float3 operator-(const float3 a, const float f)
+ccl_device_inline float3 operator-(const float3 a, const float b)
 {
-  return a - make_float3(f);
+  return a - make_float3(b);
+}
+
+ccl_device_inline float3 operator-(const float a, const float3 b)
+{
+  return make_float3(a) - b;
 }
 
 ccl_device_inline float3 operator+=(float3 &a, const float3 b)
@@ -576,7 +586,7 @@ ccl_device_inline float3 safe_pow(const float3 a, const float3 b)
   return make_float3(safe_powf(a.x, b.x), safe_powf(a.y, b.y), safe_powf(a.z, b.z));
 }
 
-ccl_device_inline auto component_wise_equal(const float3 a, const float3 b)
+ccl_device_inline auto isequal_mask(const float3 a, const float3 b)
 {
 #if defined(__KERNEL_METAL__)
   return a == b;
@@ -589,14 +599,14 @@ ccl_device_inline auto component_wise_equal(const float3 a, const float3 b)
 #endif
 }
 
-ccl_device_inline auto component_is_zero(const float3 a)
+ccl_device_inline auto is_zero_mask(const float3 a)
 {
-  return component_wise_equal(a, zero_float3());
+  return isequal_mask(a, zero_float3());
 }
 
 ccl_device_inline float3 safe_floored_fmod(const float3 a, const float3 b)
 {
-  return select(component_is_zero(b), zero_float3(), a - floor(a / b) * b);
+  return select(is_zero_mask(b), zero_float3(), a - floor(a / b) * b);
 }
 
 ccl_device_inline float3 wrap(const float3 value, const float3 max, const float3 min)
@@ -606,7 +616,7 @@ ccl_device_inline float3 wrap(const float3 value, const float3 max, const float3
 
 ccl_device_inline float3 safe_fmod(const float3 a, const float3 b)
 {
-  return select(component_is_zero(b), zero_float3(), fmod(a, b));
+  return select(is_zero_mask(b), zero_float3(), fmod(a, b));
 }
 
 ccl_device_inline float3 compatible_sign(const float3 v)