source/blender/draw/CMakeLists.txt

@@ -153,6 +153,7 @@ set(SRC
   engines/eevee_next/eevee_shader.cc
   engines/eevee_next/eevee_shadow.cc
   engines/eevee_next/eevee_sync.cc
+  engines/eevee_next/eevee_subsurface.cc
   engines/eevee_next/eevee_velocity.cc
   engines/eevee_next/eevee_view.cc
   engines/eevee_next/eevee_world.cc
@@ -284,6 +285,7 @@ set(SRC
   engines/eevee_next/eevee_shader.hh
   engines/eevee_next/eevee_shadow.hh
   engines/eevee_next/eevee_sync.hh
+  engines/eevee_next/eevee_subsurface.hh
   engines/eevee_next/eevee_velocity.hh
   engines/eevee_next/eevee_view.hh
   engines/eevee_next/eevee_world.hh
@@ -489,6 +491,7 @@ set(GLSL_SRC
   engines/eevee_next/shaders/eevee_shadow_tilemap_init_comp.glsl
   engines/eevee_next/shaders/eevee_shadow_tilemap_lib.glsl
   engines/eevee_next/shaders/eevee_spherical_harmonics_lib.glsl
+  engines/eevee_next/shaders/eevee_subsurface_eval_frag.glsl
   engines/eevee_next/shaders/eevee_surf_deferred_frag.glsl
   engines/eevee_next/shaders/eevee_surf_depth_frag.glsl
   engines/eevee_next/shaders/eevee_surf_forward_frag.glsl

source/blender/draw/engines/eevee_next/eevee_instance.cc

@@ -204,6 +204,7 @@ void Instance::end_sync()
   shadows.end_sync(); /** \note: Needs to be before lights. */
   lights.end_sync();
   sampling.end_sync();
+  subsurface.end_sync();
   film.end_sync();
   cryptomatte.end_sync();
   pipelines.end_sync();

source/blender/draw/engines/eevee_next/eevee_instance.hh

@@ -30,6 +30,7 @@
 #include "eevee_sampling.hh"
 #include "eevee_shader.hh"
 #include "eevee_shadow.hh"
+#include "eevee_subsurface.hh"
 #include "eevee_sync.hh"
 #include "eevee_view.hh"
 #include "eevee_world.hh"
@@ -48,6 +49,7 @@ class Instance {
   ShaderModule &shaders;
   SyncModule sync;
   MaterialModule materials;
+  SubsurfaceModule subsurface;
   PipelineModule pipelines;
   ShadowModule shadows;
   LightModule lights;
@@ -94,6 +96,7 @@ class Instance {
       : shaders(*ShaderModule::module_get()),
         sync(*this),
         materials(*this),
+        subsurface(*this),
         pipelines(*this),
         shadows(*this),
         lights(*this),

source/blender/draw/engines/eevee_next/eevee_shader_shared.hh

@@ -965,6 +965,7 @@ using ShadowPageCacheBuf = draw::StorageArrayBuffer<uint2, SHADOW_MAX_PAGE, true
 using ShadowTileMapDataBuf = draw::StorageVectorBuffer<ShadowTileMapData, SHADOW_MAX_TILEMAP>;
 using ShadowTileMapClipBuf = draw::StorageArrayBuffer<ShadowTileMapClip, SHADOW_MAX_TILEMAP, true>;
 using ShadowTileDataBuf = draw::StorageArrayBuffer<ShadowTileDataPacked, SHADOW_MAX_TILE, true>;
+using SubsurfaceDataBuf = draw::UniformBuffer<SubsurfaceData>;
 using VelocityGeometryBuf = draw::StorageArrayBuffer<float4, 16, true>;
 using VelocityIndexBuf = draw::StorageArrayBuffer<VelocityIndex, 16>;
 using VelocityObjectBuf = draw::StorageArrayBuffer<float4x4, 16>;

source/blender/draw/engines/eevee_next/eevee_subsurface.cc

@@ -0,0 +1,195 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright 2021 Blender Foundation.
+ */
+
+/** \file
+ * \ingroup eevee
+ *
+ */
+
+#include "BLI_vector.hh"
+
+#include "eevee_instance.hh"
+#include "eevee_subsurface.hh"
+
+#include <iostream>
+
+namespace blender::eevee {
+
+/* -------------------------------------------------------------------- */
+/** \name Subsurface
+ *
+ * \{ */
+
+/* TODO(fclem) Only enable this module if there is any SSS object in the scene. */
+void SubsurfaceModule::end_sync()
+{
+  data_.jitter_threshold = inst_.scene->eevee.sss_jitter_threshold;
+  if (data_.sample_len != inst_.scene->eevee.sss_samples) {
+    /* Convert sample count from old implementation which was using a separable filter. */
+    /* TODO(fclem) better remapping. */
+    // data_.sample_len = square_f(1 + 2 * inst_.scene->eevee.sss_samples);
+    data_.sample_len = 55;
+  }
+
+  if (transmittance_tx == nullptr) {
+    precompute_transmittance_profile();
+  }
+
+  precompute_samples_location();
+
+  data_.push_update();
+}
+
+void SubsurfaceModule::precompute_samples_location()
+{
+  /* Precompute sample position with white albedo. */
+  float d = burley_setup(1.0f, 1.0f);
+
+  float rand_u = inst_.sampling.rng_get(SAMPLING_SSS_U);
+  float rand_v = inst_.sampling.rng_get(SAMPLING_SSS_V);
+
+  double golden_angle = M_PI * (3.0 - sqrt(5.0));
+  for (auto i : IndexRange(data_.sample_len)) {
+    float theta = golden_angle * i + M_PI * 2.0f * rand_u;
+    /* Scale using rand_v in order to keep first sample always at center. */
+    float x = (1.0f + (rand_v / data_.sample_len)) * (i / (float)data_.sample_len);
+    float r = burley_sample(d, x);
+    data_.samples[i].x = cosf(theta) * r;
+    data_.samples[i].y = sinf(theta) * r;
+    data_.samples[i].z = 1.0f / burley_pdf(d, r);
+  }
+}
+
+void SubsurfaceModule::precompute_transmittance_profile()
+{
+  Vector<float> profile(SSS_TRANSMIT_LUT_SIZE);
+
+  /* Precompute sample position with white albedo. */
+  float radius = 1.0f;
+  float d = burley_setup(radius, 1.0f);
+
+  /* For each distance d we compute the radiance incoming from an hypothetical parallel plane. */
+  for (auto i : IndexRange(SSS_TRANSMIT_LUT_SIZE)) {
+    /* Distance from the lit surface plane.
+     * Compute to a larger maximum distance to have a smoother falloff for all channels. */
+    float lut_radius = SSS_TRANSMIT_LUT_RADIUS * radius;
+    float distance = lut_radius * (i + 1e-5f) / profile.size();
+    /* Compute radius of the footprint on the hypothetical plane. */
+    float r_fp = sqrtf(square_f(lut_radius) - square_f(distance));
+
+    profile[i] = 0.0f;
+    float area_accum = 0.0f;
+    for (auto j : IndexRange(SSS_TRANSMIT_LUT_STEP_RES)) {
+      /* Compute distance to the "shading" point through the medium. */
+      float r = (r_fp * (j + 0.5f)) / SSS_TRANSMIT_LUT_STEP_RES;
+      float r_prev = (r_fp * (j + 0.0f)) / SSS_TRANSMIT_LUT_STEP_RES;
+      float r_next = (r_fp * (j + 1.0f)) / SSS_TRANSMIT_LUT_STEP_RES;
+      r = hypotf(r, distance);
+      float R = burley_eval(d, r);
+      /* Since the profile and configuration are radially symmetrical we
+       * can just evaluate it once and weight it accordingly */
+      float disk_area = square_f(r_next) - square_f(r_prev);
+
+      profile[i] += R * disk_area;
+      area_accum += disk_area;
+    }
+    /* Normalize over the disk. */
+    profile[i] /= area_accum;
+  }
+  /* Make a smooth gradient from 1 to 0. */
+  float range = profile.first() - profile.last();
+  float offset = profile.last();
+  for (float &value : profile) {
+    value = (value - offset) / range;
+  }
+  profile.first() = 1;
+  profile.last() = 0;
+
+  transmittance_tx = GPU_texture_create_1d("SSSTransmittanceProfile",
+                                           profile.size(),
+                                           1,
+                                           GPU_R16F,
+                                           GPU_TEXTURE_USAGE_SHADER_READ,
+                                           profile.data());
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Christensen-Burley SSS model
+ *
+ * Based on: "Approximate Reflectance Profiles for Efficient Subsurface Scattering"
+ * by Per Christensen
+ * https://graphics.pixar.com/library/ApproxBSSRDF/approxbssrdfslides.pdf
+ * \{ */
+
+float SubsurfaceModule::burley_setup(float radius, float albedo)
+{
+  float A = albedo;
+  /* Diffuse surface transmission, equation (6). */
+  float s = 1.9f - A + 3.5f * square_f(A - 0.8f);
+  /* Mean free path length adapted to fit ancient Cubic and Gaussian models. */
+  float l = 0.25 * M_1_PI * radius;
+
+  return l / s;
+}
+
+float SubsurfaceModule::burley_sample(float d, float x_rand)
+{
+  x_rand *= SSS_BURLEY_TRUNCATE_CDF;
+
+  const float tolerance = 1e-6;
+  const int max_iteration_count = 10;
+  /* Do initial guess based on manual curve fitting, this allows us to reduce
+   * number of iterations to maximum 4 across the [0..1] range. We keep maximum
+   * number of iteration higher just to be sure we didn't miss root in some
+   * corner case.
+   */
+  float r;
+  if (x_rand <= 0.9) {
+    r = exp(x_rand * x_rand * 2.4) - 1.0;
+  }
+  else {
+    /* TODO(sergey): Some nicer curve fit is possible here. */
+    r = 15.0;
+  }
+  /* Solve against scaled radius. */
+  for (int i = 0; i < max_iteration_count; i++) {
+    float exp_r_3 = exp(-r / 3.0);
+    float exp_r = exp_r_3 * exp_r_3 * exp_r_3;
+    float f = 1.0 - 0.25 * exp_r - 0.75 * exp_r_3 - x_rand;
+    float f_ = 0.25 * exp_r + 0.25 * exp_r_3;
+
+    if (abs(f) < tolerance || f_ == 0.0) {
+      break;
+    }
+
+    r = r - f / f_;
+    if (r < 0.0) {
+      r = 0.0;
+    }
+  }
+
+  return r * d;
+}
+
+float SubsurfaceModule::burley_eval(float d, float r)
+{
+  if (r >= SSS_BURLEY_TRUNCATE * d) {
+    return 0.0;
+  }
+  /* Slide 33. */
+  float exp_r_3_d = expf(-r / (3.0f * d));
+  float exp_r_d = exp_r_3_d * exp_r_3_d * exp_r_3_d;
+  return (exp_r_d + exp_r_3_d) / (8.0f * (float)M_PI * d);
+}
+
+float SubsurfaceModule::burley_pdf(float d, float r)
+{
+  return burley_eval(d, r) / SSS_BURLEY_TRUNCATE_CDF;
+}
+
+/** \} */
+
+}  // namespace blender::eevee

source/blender/draw/engines/eevee_next/eevee_subsurface.hh

@@ -0,0 +1,69 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright 2021 Blender Foundation.
+ */
+
+/** \file
+ * \ingroup eevee
+ *
+ */
+
+#pragma once
+
+#include "eevee_shader.hh"
+#include "eevee_shader_shared.hh"
+
+namespace blender::eevee {
+
+/* -------------------------------------------------------------------- */
+/** \name Subsurface
+ *
+ * \{ */
+
+class Instance;
+
+struct SubsurfaceModule {
+ private:
+  Instance &inst_;
+  /** Contains samples locations. */
+  SubsurfaceDataBuf data_;
+  /** Contains translucence profile for a single color channel. */
+  GPUTexture *transmittance_tx = nullptr;
+
+ public:
+  SubsurfaceModule(Instance &inst) : inst_(inst)
+  {
+    /* Force first update. */
+    data_.sample_len = -1;
+  };
+
+  ~SubsurfaceModule()
+  {
+    GPU_TEXTURE_FREE_SAFE(transmittance_tx);
+  };
+
+  void end_sync();
+
+  const GPUUniformBuf *ubo_get(void) const
+  {
+    return data_;
+  }
+
+  GPUTexture **transmittance_ref_get(void)
+  {
+    return &transmittance_tx;
+  }
+
+ private:
+  void precompute_samples_location();
+  void precompute_transmittance_profile();
+
+  /** Christensen-Burley implementation. */
+  static float burley_setup(float radius, float albedo);
+  static float burley_sample(float d, float x_rand);
+  static float burley_eval(float d, float r);
+  static float burley_pdf(float d, float r);
+};
+
+/** \} */
+
+}  // namespace blender::eevee

source/blender/draw/engines/eevee_next/shaders/eevee_subsurface_eval_frag.glsl

@@ -0,0 +1,132 @@
+
+/**
+ * Postprocess diffuse radiance output from the diffuse evaluation pass to mimic subsurface
+ * transmission.
+ *
+ * This implementation follows the technique described in the siggraph presentation:
+ * "Efficient screen space subsurface scattering Siggraph 2018"
+ * by Evgenii Golubev
+ *
+ * But, instead of having all the precomputed weights for all three color primaries,
+ * we precompute a weight profile texture to be able to support per pixel AND per channel radius.
+ **/
+
+#pragma BLENDER_REQUIRE(common_view_lib.glsl)
+#pragma BLENDER_REQUIRE(common_math_geom_lib.glsl)
+#pragma BLENDER_REQUIRE(eevee_sampling_lib.glsl)
+
+vec3 burley_setup(vec3 radius, vec3 albedo)
+{
+  /* Scale albedo because we can have HDR value caused by BSDF sampling. */
+  vec3 A = albedo / max(1e-6, max_v3(albedo));
+  /* Diffuse surface transmission, equation (6). */
+  vec3 s = 1.9 - A + 3.5 * sqr(A - 0.8);
+  /* Mean free path length adapted to fit ancient Cubic and Gaussian models. */
+  vec3 l = 0.25 * M_1_PI * radius;
+
+  return l / s;
+}
+
+vec3 burley_eval(vec3 d, float r)
+{
+  /* Slide 33. */
+  vec3 exp_r_3_d = exp(-r / (3.0 * d));
+  vec3 exp_r_d = exp_r_3_d * exp_r_3_d * exp_r_3_d;
+  /** NOTE:
+   * - Surface albedo is applied at the end.
+   * - This is normalized diffuse model, so the equation is multiplied
+   *   by 2*pi, which also matches cdf().
+   */
+  return (exp_r_d + exp_r_3_d) / (4.0 * d);
+}
+
+void main(void)
+{
+  vec2 center_uv = uvcoordsvar.xy;
+  ivec2 texel = ivec2(gl_FragCoord.xy);
+
+  float gbuffer_depth = texelFetch(hiz_tx, texel, 0).r;
+  vec3 vP = get_view_space_from_depth(center_uv, gbuffer_depth);
+  vec4 tra_col_in = texelFetch(transmit_color_tx, texel, 0);
+  vec4 tra_nor_in = texelFetch(transmit_normal_tx, texel, 0);
+  vec4 tra_dat_in = texelFetch(transmit_data_tx, texel, 0);
+
+  ClosureDiffuse diffuse = gbuffer_load_diffuse_data(tra_col_in, tra_nor_in, tra_dat_in);
+
+  if (diffuse.sss_id == 0u) {
+    /* Normal diffuse is already in combined pass. */
+    /* Refraction also go into this case. */
+    out_combined = vec4(0.0);
+    return;
+  }
+
+  float max_radius = max_v3(diffuse.sss_radius);
+
+  float homcoord = ProjectionMatrix[2][3] * vP.z + ProjectionMatrix[3][3];
+  vec2 sample_scale = vec2(ProjectionMatrix[0][0], ProjectionMatrix[1][1]) *
+                      (0.5 * max_radius / homcoord);
+
+  float pixel_footprint = sample_scale.x * drw_view.viewport_size.x;
+  if (pixel_footprint <= 1.0) {
+    /* Early out. */
+    out_combined = vec4(texture(radiance_tx, center_uv).rgb * diffuse.color, 0.0);
+    return;
+  }
+
+  diffuse.sss_radius = max(vec3(1e-4), diffuse.sss_radius / max_radius) * max_radius;
+  vec3 d = burley_setup(diffuse.sss_radius, diffuse.color);
+
+  /* Do not rotate too much to avoid too much cache misses. */
+  float golden_angle = M_PI * (3.0 - sqrt(5.0));
+  float theta = interlieved_gradient_noise(gl_FragCoord.xy, 0, 0.0) * golden_angle;
+  float cos_theta = cos(theta);
+  float sin_theta = sqrt(1.0 - sqr(cos_theta));
+  mat2 rot = mat2(cos_theta, sin_theta, -sin_theta, cos_theta);
+
+  mat2 scale = mat2(sample_scale.x, 0.0, 0.0, sample_scale.y);
+  mat2 sample_space = scale * rot;
+
+  vec3 accum_weight = vec3(0.0);
+  vec3 accum = vec3(0.0);
+
+  /* TODO/OPTI(fclem) Make separate sample set for lower radius. */
+
+  for (int i = 0; i < sss_buf.sample_len; i++) {
+    vec2 sample_uv = center_uv + sample_space * sss_buf.samples[i].xy;
+    float pdf_inv = sss_buf.samples[i].z;
+
+    float sample_depth = textureLod(hiz_tx, sample_uv * hiz_buf.uv_scale, 0.0).r;
+    vec3 sample_vP = get_view_space_from_depth(sample_uv, sample_depth);
+
+    vec4 sample_data = texture(radiance_tx, sample_uv);
+    vec3 sample_radiance = sample_data.rgb;
+    uint sample_sss_id = uint(sample_data.a);
+
+    if (sample_sss_id != diffuse.sss_id) {
+      continue;
+    }
+
+    /* Discard out of bounds samples. */
+    if (any(lessThan(sample_uv, vec2(0.0))) || any(greaterThan(sample_uv, vec2(1.0)))) {
+      continue;
+    }
+
+    /* Slide 34. */
+    float r = distance(sample_vP, vP);
+    vec3 weight = burley_eval(d, r) * pdf_inv;
+
+    accum += sample_radiance * weight;
+    accum_weight += weight;
+  }
+  /* Normalize the sum (slide 34). */
+  accum /= accum_weight;
+  /* Apply surface color on final radiance. */
+  accum *= diffuse.color;
+
+  /* Debug, detect NaNs. */
+  if (any(isnan(accum))) {
+    accum = vec3(1.0, 0.0, 1.0);
+  }
+
+  out_combined = vec4(accum, 0.0);
+}

source/blender/draw/engines/eevee_next/shaders/infos/eevee_material_info.hh

@@ -180,6 +180,25 @@ GPU_SHADER_CREATE_INFO(eevee_surf_shadow)
     .fragment_source("eevee_surf_shadow_frag.glsl")
     .additional_info("eevee_camera", "eevee_utility_texture", "eevee_sampling_data");
 
+GPU_SHADER_CREATE_INFO(eevee_transmittance_data)
+    .define("SSS_TRANSMITTANCE")
+    .sampler(0, ImageType::FLOAT_1D, "sss_transmittance_tx");
+
+GPU_SHADER_CREATE_INFO(eevee_subsurface_eval)
+    .do_static_compilation(true)
+    .additional_info("eevee_shared")
+    .uniform_buf(2, "SubsurfaceData", "sss_buf")
+    .uniform_buf(1, "HiZData", "hiz_buf")
+    .sampler(0, ImageType::FLOAT_2D, "hiz_tx")
+    .sampler(1, ImageType::FLOAT_2D, "radiance_tx")
+    .sampler(2, ImageType::FLOAT_2D, "transmit_color_tx")
+    .sampler(3, ImageType::FLOAT_2D, "transmit_normal_tx")
+    .sampler(4, ImageType::FLOAT_2D, "transmit_data_tx")
+    .fragment_out(0, Type::VEC4, "out_combined")
+    .fragment_source("eevee_subsurface_eval_frag.glsl")
+    /* TODO(fclem) Output to diffuse pass without feedback loop. */
+    .additional_info("draw_fullscreen", "draw_view");
+
 #undef image_out
 #undef image_array_out