Merge remote-tracking branch 'origin/blender-v2.90-release' into soc-2020-production-ready-light-tree-2

* removed unused update light picking code
* switched all light tree function parameters to be in the same P, V, t order
* renamed N_pick to V_pick since it is not always a normal

build_files/build_environment/install_deps.sh

@@ -72,7 +72,7 @@ COMMANDLINE=$@
 
 DISTRO=""
 RPM=""
-SRC="$HOME/src/blender-deps"
+SRC="$HOME/Programming/blender-git/src/blender-deps"
 INST="/opt/lib"
 TMP="/tmp"
 CWD=$PWD

build_files/windows/configure_msbuild.cmd

@@ -77,4 +77,4 @@ echo msbuild ^
 	/verbosity:minimal ^
 	/p:platform=%MSBUILD_PLATFORM% ^
 	/flp:Summary;Verbosity=minimal;LogFile=%BUILD_DIR%\Build.log >> %BUILD_DIR%\rebuild.cmd
-echo echo %%TIME%% ^>^> buildtime.txt >> %BUILD_DIR%\rebuild.cmd
+echo echo %%TIME%% ^>^> buildtime.txt >> %BUILD_DIR%\rebuild.cmd

build_files/windows/configure_ninja.cmd

@@ -93,4 +93,4 @@ echo   call "%VCVARS%" %BUILD_ARCH% >> %BUILD_DIR%\rebuild.cmd
 echo ^) >> %BUILD_DIR%\rebuild.cmd
 echo echo %%TIME%% ^> buildtime.txt >> %BUILD_DIR%\rebuild.cmd
 echo ninja install >> %BUILD_DIR%\rebuild.cmd 
-echo echo %%TIME%% ^>^> buildtime.txt >> %BUILD_DIR%\rebuild.cmd
+echo echo %%TIME%% ^>^> buildtime.txt >> %BUILD_DIR%\rebuild.cmd

intern/cycles/blender/addon/presets.py

@@ -68,6 +68,8 @@ class AddPresetSampling(AddPresetBase, Operator):
         "cycles.subsurface_samples",
         "cycles.volume_samples",
         "cycles.use_square_samples",
+        "cycles.use_light_tree",
+        "cycles.splitting_threshold",
         "cycles.progressive",
         "cycles.seed",
         "cycles.sample_clamp_direct",

intern/cycles/blender/addon/properties.py

@@ -291,6 +291,19 @@ class CyclesRenderSettings(bpy.types.PropertyGroup):
         default=False,
     )
 
+    use_light_tree: BoolProperty(
+        name="Light Tree",
+        description="Samples many lights more efficiently",
+        default=False,
+    )
+
+    splitting_threshold: FloatProperty(
+        name="Splitting",
+        description="Amount of lights to sample at a time, from one light at 0.0, to adaptively more lights as needed, to all lights at 1.0",
+        min=0.0, max=1.0,
+        default=0.0,
+    )
+
     samples: IntProperty(
         name="Samples",
         description="Number of samples to render for each pixel",

intern/cycles/blender/addon/ui.py

@@ -326,6 +326,35 @@ class CYCLES_RENDER_PT_sampling_advanced(CyclesButtonsPanel, Panel):
                 layout.row().prop(cscene, "use_layer_samples")
                 break
 
+class CYCLES_RENDER_PT_sampling_light_tree(CyclesButtonsPanel, Panel):
+    bl_label = "Light Tree"
+    bl_parent_id = "CYCLES_RENDER_PT_sampling"
+    bl_options = {'DEFAULT_CLOSED'}
+
+    @classmethod
+    def poll(cls, context):
+        return (context.scene.cycles.feature_set == 'EXPERIMENTAL')
+
+    def draw_header(self, context):
+        layout = self.layout
+        scene = context.scene
+        cscene = scene.cycles
+
+        layout.prop(cscene, "use_light_tree", text="")
+
+    def draw(self, context):
+        layout = self.layout
+        layout.use_property_split = True
+        layout.use_property_decorate = False
+
+        scene = context.scene
+        cscene = scene.cycles
+
+        layout.active = cscene.use_light_tree
+        row = layout.row(align=True)
+        row.label(text="") # create empty column
+        row.prop(cscene, "splitting_threshold", text="Splitting")
+
 
 class CYCLES_RENDER_PT_sampling_total(CyclesButtonsPanel, Panel):
     bl_label = "Total Samples"
@@ -2272,6 +2301,7 @@ classes = (
     CYCLES_RENDER_PT_sampling_adaptive,
     CYCLES_RENDER_PT_sampling_denoising,
     CYCLES_RENDER_PT_sampling_advanced,
+    CYCLES_RENDER_PT_sampling_light_tree,
     CYCLES_RENDER_PT_light_paths,
     CYCLES_RENDER_PT_light_paths_max_bounces,
     CYCLES_RENDER_PT_light_paths_clamping,

intern/cycles/blender/blender_sync.cpp

@@ -304,6 +304,8 @@ void BlenderSync::sync_integrator()
   integrator->method = (Integrator::Method)get_enum(
       cscene, "progressive", Integrator::NUM_METHODS, Integrator::PATH);
 
+  integrator->use_light_tree = get_boolean(cscene, "use_light_tree");
+  integrator->splitting_threshold = get_float(cscene, "splitting_threshold");
   integrator->sample_all_lights_direct = get_boolean(cscene, "sample_all_lights_direct");
   integrator->sample_all_lights_indirect = get_boolean(cscene, "sample_all_lights_indirect");
   integrator->light_sampling_threshold = get_float(cscene, "light_sampling_threshold");
@@ -359,8 +361,13 @@ void BlenderSync::sync_integrator()
     integrator->ao_bounces = 0;
   }
 
-  if (integrator->modified(previntegrator))
+  if (integrator->modified(previntegrator)) {
     integrator->tag_update(scene);
+  }
+  if (integrator->use_light_tree != previntegrator.use_light_tree ||
+      integrator->splitting_threshold != previntegrator.splitting_threshold) {
+    scene->light_manager->tag_update(scene);
+  }
 }
 
 /* Film */

intern/cycles/kernel/kernel_bake.h

@@ -81,6 +81,9 @@ ccl_device_noinline void compute_light_pass(
               kg, sd, emission_sd, L, &state, &ray, &throughput, &ss_indirect)) {
         while (ss_indirect.num_rays) {
           kernel_path_subsurface_setup_indirect(kg, &ss_indirect, &state, &ray, L, &throughput);
+          indirect_sd.P_pick = sd->P_pick;
+          indirect_sd.V_pick = sd->V_pick;
+          indirect_sd.t_pick = sd->t_pick;
           kernel_path_indirect(kg, &indirect_sd, emission_sd, &ray, throughput, &state, L);
         }
         is_sss_sample = true;
@@ -97,6 +100,9 @@ ccl_device_noinline void compute_light_pass(
         state.ray_t = 0.0f;
 #  endif
         /* compute indirect light */
+        indirect_sd.P_pick = sd->P_pick;
+        indirect_sd.V_pick = sd->V_pick;
+        indirect_sd.t_pick = sd->t_pick;
         kernel_path_indirect(kg, &indirect_sd, emission_sd, &ray, throughput, &state, L);
 
         /* sum and reset indirect light pass variables for the next samples */

intern/cycles/kernel/kernel_emission.h

@@ -222,6 +222,7 @@ ccl_device_noinline_cpu float3 indirect_primitive_emission(
     /* multiple importance sampling, get triangle light pdf,
      * and compute weight with respect to BSDF pdf */
     float pdf = triangle_light_pdf(kg, sd, t);
+    pdf *= light_distribution_pdf(kg, sd->P_pick, sd->V_pick, sd->t_pick, sd->prim, sd->object);
     float mis_weight = power_heuristic(bsdf_pdf, pdf);
 
     return L * mis_weight;
@@ -270,11 +271,16 @@ ccl_device_noinline_cpu void indirect_lamp_emission(KernelGlobals *kg,
       kernel_volume_shadow(kg, emission_sd, state, &volume_ray, &volume_tp);
       lamp_L *= volume_tp;
     }
+
 #endif
 
     if (!(state->flag & PATH_RAY_MIS_SKIP)) {
       /* multiple importance sampling, get regular light pdf,
        * and compute weight with respect to BSDF pdf */
+
+      /* multiply with light picking probablity to pdf */
+      ls.pdf *= light_distribution_pdf(
+          kg, emission_sd->P_pick, emission_sd->V_pick, emission_sd->t_pick, ~ls.lamp, -1);
       float mis_weight = power_heuristic(state->ray_pdf, ls.pdf);
       lamp_L *= mis_weight;
     }
@@ -326,10 +332,14 @@ ccl_device_noinline_cpu float3 indirect_background(KernelGlobals *kg,
   /* Background MIS weights. */
 #  ifdef __BACKGROUND_MIS__
   /* Check if background light exists or if we should skip pdf. */
+
+  /* consider shading point at previous non-transparent bounce */
+  float3 P_pick = ray->P - state->ray_t * ray->D;
+
   if (!(state->flag & PATH_RAY_MIS_SKIP) && kernel_data.background.use_mis) {
     /* multiple importance sampling, get background light pdf for ray
      * direction, and compute weight with respect to BSDF pdf */
-    float pdf = background_light_pdf(kg, ray->P, ray->D);
+    float pdf = background_light_pdf(kg, P_pick, ray->D);
     float mis_weight = power_heuristic(state->ray_pdf, pdf);
 
     return L * mis_weight;

intern/cycles/kernel/kernel_light.h

@@ -35,8 +35,31 @@ typedef struct LightSample {
   LightType type; /* type of light */
 } LightSample;
 
+/* Updates the position and normal used to pick a light for direct lighting.
+ *
+ * The importance calculation for the light tree is different for scattering events on a surface
+ * and in a volume. For volume events we need to know the point where the ray started and the point
+ * it scattered. To do this we keep track of the ray direction and length. For surface events we
+ * need the normal at the surface. In this case we set the ray length to -1 to mark that surface
+ * importance sampling should be used.
+ */
+ccl_device void kernel_update_light_picking(ShaderData *sd, Ray *ray)
+{
+  if (ray) {
+    sd->P_pick = ray->P;
+    sd->V_pick = ray->D;
+    sd->t_pick = ray->t;
+  }
+  else {
+    sd->P_pick = sd->P;
+    sd->V_pick = sd->N;
+    sd->t_pick = -1.0f;
+  }
+}
+
 /* Regular Light */
 
+/* returns the PDF of sampling a point on this lamp */
 ccl_device_inline bool lamp_light_sample(
     KernelGlobals *kg, int lamp, float randu, float randv, float3 P, LightSample *ls)
 {
@@ -153,8 +176,6 @@ ccl_device_inline bool lamp_light_sample(
     }
   }
 
-  ls->pdf *= kernel_data.integrator.pdf_lights;
-
   return (ls->pdf > 0.0f);
 }
 
@@ -291,8 +312,6 @@ ccl_device bool lamp_light_eval(
     return false;
   }
 
-  ls->pdf *= kernel_data.integrator.pdf_lights;
-
   return true;
 }
 
@@ -328,18 +347,15 @@ ccl_device_inline bool triangle_world_space_vertices(
   return has_motion;
 }
 
-ccl_device_inline float triangle_light_pdf_area(KernelGlobals *kg,
-                                                const float3 Ng,
-                                                const float3 I,
-                                                float t)
+ccl_device_inline float triangle_light_pdf_area(
+    KernelGlobals *kg, const float3 Ng, const float3 I, float t, float triangle_area)
 {
-  float pdf = kernel_data.integrator.pdf_triangles;
   float cos_pi = fabsf(dot(Ng, I));
 
   if (cos_pi == 0.0f)
     return 0.0f;
 
-  return t * t * pdf / cos_pi;
+  return t * t / (cos_pi * triangle_area);
 }
 
 ccl_device_forceinline float triangle_light_pdf(KernelGlobals *kg, ShaderData *sd, float t)
@@ -349,7 +365,7 @@ ccl_device_forceinline float triangle_light_pdf(KernelGlobals *kg, ShaderData *s
    * to the length of the edges of the triangle. */
 
   float3 V[3];
-  bool has_motion = triangle_world_space_vertices(kg, sd->object, sd->prim, sd->time, V);
+  triangle_world_space_vertices(kg, sd->object, sd->prim, sd->time, V);
 
   const float3 e0 = V[1] - V[0];
   const float3 e1 = V[2] - V[0];
@@ -380,34 +396,12 @@ ccl_device_forceinline float triangle_light_pdf(KernelGlobals *kg, ShaderData *s
       return 0.0f;
     }
     else {
-      float area = 1.0f;
-      if (has_motion) {
-        /* get the center frame vertices, this is what the PDF was calculated from */
-        triangle_world_space_vertices(kg, sd->object, sd->prim, -1.0f, V);
-        area = triangle_area(V[0], V[1], V[2]);
-      }
-      else {
-        area = 0.5f * len(N);
-      }
-      const float pdf = area * kernel_data.integrator.pdf_triangles;
-      return pdf / solid_angle;
+      return 1.0f / solid_angle;
     }
   }
   else {
-    float pdf = triangle_light_pdf_area(kg, sd->Ng, sd->I, t);
-    if (has_motion) {
-      const float area = 0.5f * len(N);
-      if (UNLIKELY(area == 0.0f)) {
-        return 0.0f;
-      }
-      /* scale the PDF.
-       * area = the area the sample was taken from
-       * area_pre = the are from which pdf_triangles was calculated from */
-      triangle_world_space_vertices(kg, sd->object, sd->prim, -1.0f, V);
-      const float area_pre = triangle_area(V[0], V[1], V[2]);
-      pdf = pdf * area_pre / area;
-    }
-    return pdf;
+    const float area = 0.5f * len(N);
+    return triangle_light_pdf_area(kg, sd->Ng, sd->I, t, area);
   }
 }
 
@@ -425,7 +419,7 @@ ccl_device_forceinline void triangle_light_sample(KernelGlobals *kg,
    * to the length of the edges of the triangle. */
 
   float3 V[3];
-  bool has_motion = triangle_world_space_vertices(kg, object, prim, time, V);
+  triangle_world_space_vertices(kg, object, prim, time, V);
 
   const float3 e0 = V[1] - V[0];
   const float3 e1 = V[2] - V[0];
@@ -434,7 +428,6 @@ ccl_device_forceinline void triangle_light_sample(KernelGlobals *kg,
   const float3 N0 = cross(e0, e1);
   float Nl = 0.0f;
   ls->Ng = safe_normalize_len(N0, &Nl);
-  float area = 0.5f * Nl;
 
   /* flip normal if necessary */
   const int object_flag = kernel_tex_fetch(__object_flag, object);
@@ -537,13 +530,7 @@ ccl_device_forceinline void triangle_light_sample(KernelGlobals *kg,
       return;
     }
     else {
-      if (has_motion) {
-        /* get the center frame vertices, this is what the PDF was calculated from */
-        triangle_world_space_vertices(kg, object, prim, -1.0f, V);
-        area = triangle_area(V[0], V[1], V[2]);
-      }
-      const float pdf = area * kernel_data.integrator.pdf_triangles;
-      ls->pdf = pdf / solid_angle;
+      ls->pdf = 1.0f / solid_angle;
     }
   }
   else {
@@ -564,20 +551,53 @@ ccl_device_forceinline void triangle_light_sample(KernelGlobals *kg,
     ls->P = u * V[0] + v * V[1] + t * V[2];
     /* compute incoming direction, distance and pdf */
     ls->D = normalize_len(ls->P - P, &ls->t);
-    ls->pdf = triangle_light_pdf_area(kg, ls->Ng, -ls->D, ls->t);
-    if (has_motion && area != 0.0f) {
-      /* scale the PDF.
-       * area = the area the sample was taken from
-       * area_pre = the are from which pdf_triangles was calculated from */
-      triangle_world_space_vertices(kg, object, prim, -1.0f, V);
-      const float area_pre = triangle_area(V[0], V[1], V[2]);
-      ls->pdf = ls->pdf * area_pre / area;
-    }
+    float area = 0.5f * Nl;
+    ls->pdf = triangle_light_pdf_area(kg, ls->Ng, -ls->D, ls->t, area);
+
     ls->u = u;
     ls->v = v;
   }
 }
 
+/* chooses either to sample the light tree, distant or background lights by
+ * sampling a CDF based on energy */
+ccl_device int light_group_distribution_sample(KernelGlobals *kg, float *randu)
+{
+  /* This is basically std::upper_bound as used by pbrt, to find a point light or
+   * triangle to emit from, proportional to area. a good improvement would be to
+   * also sample proportional to power, though it's not so well defined with
+   * arbitrary shaders. */
+  const int num_groups = LIGHTGROUP_NUM;
+  int first = 0;
+  int len = num_groups + 1;
+  float r = *randu;
+  // todo: refactor this into its own function. It is used in several places
+  while (len > 0) {
+    int half_len = len >> 1;
+    int middle = first + half_len;
+
+    if (r < kernel_tex_fetch(__light_group_sample_cdf, middle)) {
+      len = half_len;
+    }
+    else {
+      first = middle + 1;
+      len = len - half_len - 1;
+    }
+  }
+
+  /* Clamping should not be needed but float rounding errors seem to
+   * make this fail on rare occasions. */
+  int index = clamp(first - 1, 0, num_groups - 1);
+
+  /* Rescale to reuse random number. this helps the 2D samples within
+   * each area light be stratified as well. */
+  float distr_min = kernel_tex_fetch(__light_group_sample_cdf, index);
+  float distr_max = kernel_tex_fetch(__light_group_sample_cdf, index + 1);
+  *randu = (r - distr_min) / (distr_max - distr_min);
+
+  return index;
+}
+
 /* Light Distribution */
 
 ccl_device int light_distribution_sample(KernelGlobals *kg, float *randu)
@@ -623,22 +643,691 @@ ccl_device_inline bool light_select_reached_max_bounces(KernelGlobals *kg, int i
   return (bounce > kernel_tex_fetch(__lights, index).max_bounces);
 }
 
-ccl_device_noinline bool light_sample(KernelGlobals *kg,
-                                      int lamp,
-                                      float randu,
-                                      float randv,
-                                      float time,
-                                      float3 P,
-                                      int bounce,
-                                      LightSample *ls)
+/*
+ * Finds the solid angle of the smallest cone with vertex P that contains the bounding box. If P is
+ * inside the bounding box light can be in any direction so use the entire sphere.
+ */
+ccl_device float calc_bbox_solid_angle(float3 P,
+                                       float3 centroid_to_P_dir,
+                                       float3 bboxMin,
+                                       float3 bboxMax)
+{
+  if (P.x < bboxMax.x && P.y < bboxMax.y && P.z < bboxMax.z && P.x > bboxMin.x &&
+      P.y > bboxMin.y && P.z > bboxMin.z) {
+    /* P is inside bounding box */
+    return M_PI_F;
+  }
+  else {
+    /* Find the smallest cone that contains the bounding box by checking which bbox vertex is
+     * farthest out. If we use a bounding sphere we get a too big cone. For example consider a long
+     * skinny bbox oriented with P next to one of the small sides. */
+    float theta_u = 0;
+    float3 corners[8];
+    corners[0] = bboxMin;
+    corners[1] = make_float3(bboxMin.x, bboxMin.y, bboxMax.z);
+    corners[2] = make_float3(bboxMin.x, bboxMax.y, bboxMin.z);
+    corners[3] = make_float3(bboxMin.x, bboxMax.y, bboxMax.z);
+    corners[4] = make_float3(bboxMax.x, bboxMin.y, bboxMin.z);
+    corners[5] = make_float3(bboxMax.x, bboxMin.y, bboxMax.z);
+    corners[6] = make_float3(bboxMax.x, bboxMax.y, bboxMin.z);
+    corners[7] = bboxMax;
+    for (int i = 0; i < 8; ++i) {
+      float3 P_to_corner = normalize(P - corners[i]);
+      const float cos_theta_u = dot(-centroid_to_P_dir, P_to_corner);
+      theta_u = fmaxf(fast_acosf(cos_theta_u), theta_u);
+    }
+    return theta_u;
+  }
+}
+
+/* calculates the importance metric for the given node and shading point P
+ * If t_max is negative assume that V is the normal vector and use surface importance calculation.
+ * Otherwise assume that V is the ray direction and use volume importance.
+ */
+ccl_device float calc_importance(KernelGlobals *kg,
+                                 float3 P,
+                                 float3 V,
+                                 float t_max,
+                                 float3 bboxMax,
+                                 float3 bboxMin,
+                                 float theta_o,
+                                 float theta_e,
+                                 float3 axis,
+                                 float energy,
+                                 float3 centroid)
+{
+  if (t_max < 0.0f) {
+    /* eq. 3 */
+
+    const float3 centroid_to_P = P - centroid;
+    const float3 centroid_to_P_dir = normalize(centroid_to_P);
+    const float r2 = len_squared(bboxMax - centroid);
+    float d2 = len_squared(centroid_to_P);
+
+    /* based on comment in the implementation details of the paper */
+    const bool splitting = kernel_data.integrator.splitting_threshold != 0.0f;
+    if (!splitting) {
+      d2 = fmaxf(d2, r2 * 0.25f);
+    }
+
+    /* "theta_u captures the solid angle of the entire box" */
+
+    float theta_u = calc_bbox_solid_angle(P, centroid_to_P_dir, bboxMin, bboxMax);
+
+    /* cos(theta') */
+    float cos_theta = dot(axis, centroid_to_P_dir);
+    const float theta = fast_acosf(cos_theta);
+    const float theta_prime = fmaxf(theta - theta_o - theta_u, 0.0f);
+    if (theta_prime >= theta_e) {
+      return 0.0f;
+    }
+    const float cos_theta_prime = fast_cosf(theta_prime);
+
+    /* f_a|cos(theta'_i)| -- diffuse approximation */
+    const float cos_theta_i = dot(V, -centroid_to_P_dir);
+    const float theta_i = fast_acosf(cos_theta_i);
+    const float theta_i_prime = fmaxf(theta_i - theta_u, 0.0f);
+    const float cos_theta_i_prime = fast_cosf(theta_i_prime);
+    const float abs_cos_theta_i_prime = fabsf(cos_theta_i_prime);
+    /* doing something similar to bsdf_diffuse_eval_reflect() */
+    /* TODO: Use theta_i or theta_i_prime here? */
+    const float f_a = fmaxf(cos_theta_i_prime, 0.0f) * M_1_PI_F;
+
+    return f_a * abs_cos_theta_i_prime * energy * cos_theta_prime / d2;
+  }
+  else {
+    const float3 p_to_c = centroid - P;
+    /* find closest point to centroid */
+    const float t = fminf(t_max, fmaxf(0.0f, dot(V, p_to_c)));
+    const float3 geometry_point = P + V * t;
+    const float d_min = len(centroid - geometry_point);
+
+    const float3 V0 = normalize(P - centroid);
+    const float3 V1 = normalize(P + V * fminf(t_max, 1e12f) - centroid);
+    const float3 O0 = V0;
+    float3 O1, O2;
+    make_orthonormals_tangent(O0, V1, &O1, &O2);
+
+    const float b_max = fmaxf(dot(V0, axis), dot(V1, axis));
+    const float O0_dot_a = dot(O0, axis);
+    const float O1_dot_a = dot(O1, axis);
+    const float cos_phi_o = O0_dot_a / sqrtf(O1_dot_a * O1_dot_a + O0_dot_a * O0_dot_a); /* Eq 5 */
+    const float sin_phi_o = sqrtf(1 - cos_phi_o * cos_phi_o);
+    const float3 v = O0 * cos_phi_o + O1 * sin_phi_o;
+
+    /* Eq 6 */
+    float cos_theta_min;
+    if (O1_dot_a < 0 || dot(V0, V1) > cos_phi_o) {
+      cos_theta_min = b_max;
+    }
+    else {
+      cos_theta_min = dot(axis, v);
+    }
+
+    float theta_min = fast_acosf(cos_theta_min);
+    float theta_u = calc_bbox_solid_angle(
+        geometry_point, normalize(geometry_point - centroid), bboxMin, bboxMax);
+    float theta_prime = fmaxf(theta_min - theta_o - theta_u, 0);
+
+    if (theta_prime >= theta_e) {
+      return 0;
+    }
+
+    return energy * fast_cosf(theta_prime) / d_min; /* Eq 7 */
+  }
+}
+
+/* the energy, spatial and orientation bounds for the light are loaded and decoded
+ * and then this information is used to calculate the importance for this light.
+ * This function is used to calculate the importance for a light in a node
+ * containing several lights. */
+ccl_device float calc_light_importance(
+    KernelGlobals *kg, float3 P, float3 V, float t_max, int node_offset, int light_offset)
+{
+  /* find offset into light_tree_leaf_emitters array */
+  int first_emitter = kernel_tex_fetch(__leaf_to_first_emitter, node_offset / 4);
+  kernel_assert(first_emitter != -1);
+  int offset = first_emitter + light_offset * 3;
+
+  /* get relevant information to be able to calculate the importance */
+  const float4 data0 = kernel_tex_fetch(__light_tree_leaf_emitters, offset + 0);
+  const float4 data1 = kernel_tex_fetch(__light_tree_leaf_emitters, offset + 1);
+  const float4 data2 = kernel_tex_fetch(__light_tree_leaf_emitters, offset + 2);
+
+  /* decode data for this light */
+  const float3 bbox_min = make_float3(data0.x, data0.y, data0.z);
+  const float3 bbox_max = make_float3(data0.w, data1.x, data1.y);
+  const float theta_o = data1.z;
+  const float theta_e = data1.w;
+  const float3 axis = make_float3(data2.x, data2.y, data2.z);
+  const float energy = data2.w;
+  const float3 centroid = 0.5f * (bbox_max + bbox_min);
+
+  return calc_importance(
+      kg, P, V, t_max, bbox_max, bbox_min, theta_o, theta_e, axis, energy, centroid);
+}
+
+/* the combined energy, spatial and orientation bounds for all the lights for the
+ * given node are loaded and decoded and then this information is used to
+ * calculate the importance for this node. */
+ccl_device float calc_node_importance(
+    KernelGlobals *kg, float3 P, float3 V, float t_max, int node_offset)
+{
+  /* load the data for this node */
+  const float4 node0 = kernel_tex_fetch(__light_tree_nodes, node_offset + 0);
+  const float4 node1 = kernel_tex_fetch(__light_tree_nodes, node_offset + 1);
+  const float4 node2 = kernel_tex_fetch(__light_tree_nodes, node_offset + 2);
+  const float4 node3 = kernel_tex_fetch(__light_tree_nodes, node_offset + 3);
+
+  /* decode the data so it can be used to calculate the importance */
+  const float energy = node0.x;
+  const float3 bbox_min = make_float3(node1.x, node1.y, node1.z);
+  const float3 bbox_max = make_float3(node1.w, node2.x, node2.y);
+  const float theta_o = node2.z;
+  const float theta_e = node2.w;
+  const float3 axis = make_float3(node3.x, node3.y, node3.z);
+  const float3 centroid = 0.5f * (bbox_max + bbox_min);
+
+  return calc_importance(
+      kg, P, V, t_max, bbox_max, bbox_min, theta_o, theta_e, axis, energy, centroid);
+}
+
+/* given a node offset, this function loads and decodes the minimum amount of
+ * data needed for a the given node to be able to only either identify if it is
+ * a leaf node or how to find its two children
+ *
+ * child_o  ffset is an offset into the nodes array to this nodes right child. the
+ * left child has index node_offset+4.
+ * distribution_id corresponds to an offset into the distribution array for the
+ * first light contained in this node. num_emitters is how many lights there are
+ * in this node. */
+ccl_device void update_node(
+    KernelGlobals *kg, int node_offset, int *child_offset, int *distribution_id, int *num_emitters)
+{
+  float4 node = kernel_tex_fetch(__light_tree_nodes, node_offset);
+  (*child_offset) = __float_as_int(node.y);
+  (*distribution_id) = __float_as_int(node.z);
+  (*num_emitters) = __float_as_int(node.w);
+}
+
+/* picks one of the distant lights and computes the probability of picking it */
+ccl_device void light_distant_sample(
+    KernelGlobals *kg, float3 P, float *randu, int *index, float *pdf)
+{
+  light_distribution_sample(kg, randu);  // rescale random number
+
+  /* choose one of the distant lights randomly */
+  int num_distant = kernel_data.integrator.num_distant_lights;
+  int light = min((int)(*randu * (float)num_distant), num_distant - 1);
+
+  /* This assumes the distant lights are next to each other in the
+   * distribution array starting at distant_lights_offset. */
+  int distant_lights_offset = kernel_data.integrator.distant_lights_offset;
+
+  *index = light + distant_lights_offset;
+  *pdf = kernel_data.integrator.inv_num_distant_lights;
+}
+
+/* picks the background light and sets the probability of picking it */
+ccl_device void light_background_sample(
+    KernelGlobals *kg, float3 P, float *randu, int *index, float *pdf)
+{
+  *index = kernel_tex_fetch(__lamp_to_distribution, kernel_data.integrator.background_light_index);
+  *pdf = 1.0f;
+}
+
+/* picks a light from the light tree and returns its index and the probability of
+ * picking this light. */
+ccl_device void light_tree_sample(KernelGlobals *kg,
+                                  float3 P,
+                                  float3 V,
+                                  float t_max,
+                                  float *randu,
+                                  int *index,
+                                  float *pdf_factor)
+{
+  int sampled_index = -1;
+  *pdf_factor = 1.0f;
+
+  int offset = 0;
+  int right_child_offset, distribution_id, num_emitters;
+  do {
+
+    /* read in first part of node of light tree */
+    update_node(kg, offset, &right_child_offset, &distribution_id, &num_emitters);
+
+    /* Found a leaf - Choose which light to use */
+    if (right_child_offset == -1) {  // Found a leaf
+      if (num_emitters == 1) {
+        sampled_index = distribution_id;
+      }
+      else {
+
+        /* At a leaf node containing several lights. Pick one of these
+         * by creating and sampling a CDF based on the importance metric.
+         *
+         * The number of lights in this leaf node is not known at compile
+         * time and dynamic allocation is not allowed on the GPU, so
+         * some more computations have to be done instead.
+         * (TODO: Could we allocate a fixed array of the same size as
+         * the maximum allowed number of lights per node in the
+         * construction algorithm? i.e. max_lights_in_node.)
+         *
+         * First, the total importance of all the lights are calculated.
+         * Then, a linear loop over the lights are done where the
+         * current CDF value is calculated. This loop can stop as soon
+         * as the random value used to sample the CDF is less than the
+         * current CDF value. The sampled light has index i-1 if i is
+         * the iteration counter of the loop over the lights. This is
+         * similar to for example light the old light_distribution_sample()
+         * except not having an array to store the CDF in. */
+        float sum = 0.0f;
+        for (int i = 0; i < num_emitters; ++i) {
+          sum += calc_light_importance(kg, P, V, t_max, offset, i);
+        }
+
+        if (sum == 0.0f) {
+          *pdf_factor = 0.0f;
+          return;
+        }
+
+        float sum_inv = 1.0f / sum;
+
+        float cdf_L = 0.0f;
+        float cdf_R = 0.0f;
+        float prob = 0.0f;
+        int light = 0;
+        for (int i = 1; i < num_emitters + 1; ++i) {
+          prob = calc_light_importance(kg, P, V, t_max, offset, i - 1) * sum_inv;
+          cdf_R = cdf_L + prob;
+          if (*randu < cdf_R) {
+            light = i - 1;
+            break;
+          }
+          cdf_L = cdf_R;
+        }
+
+        sampled_index = distribution_id + light;
+        *pdf_factor *= prob;
+        /* rescale random number */
+        *randu = (*randu - cdf_L) / (cdf_R - cdf_L);
+      }
+      break;
+    }
+    else {  // Interior node, pick left or right randomly
+
+      /* calculate probability of going down left node */
+      int child_offsetL = offset + 4;
+      int child_offsetR = 4 * right_child_offset;
+      float I_L = calc_node_importance(kg, P, V, t_max, child_offsetL);
+      float I_R = calc_node_importance(kg, P, V, t_max, child_offsetR);
+      if ((I_L == 0.0f) && (I_R == 0.0f)) {
+        *pdf_factor = 0.0f;
+        break;
+      }
+
+      float P_L = I_L / (I_L + I_R);
+
+      /* choose which node to go down */
+      if (*randu <= P_L) {  // Going down left node
+        /* rescale random number */
+        *randu = *randu / P_L;
+
+        offset = child_offsetL;
+        *pdf_factor *= P_L;
+      }
+      else {  // Going down right node
+        /* rescale random number */
+        *randu = (*randu * (I_L + I_R) - I_L) / I_R;
+
+        offset = child_offsetR;
+        *pdf_factor *= 1.0f - P_L;
+      }
+    }
+  } while (true);
+
+  *index = sampled_index;
+}
+
+/* converts from an emissive triangle index to the corresponding
+ * light distribution index. */
+ccl_device int triangle_to_distribution(KernelGlobals *kg, int triangle_id, int object_id)
+{
+  /* binary search to find triangle_id which then gives distribution_id */
+  /* equivalent to implementation of std::lower_bound */
+  /* todo: of complexity log(N) now. could be made constant with a hash table? */
+  /* __triangle_to_distribution is an array of uints of the format below:
+   * [triangle_id0, object_id0, distribution_id0, triangle_id1,... ]
+   * where e.g. [triangle_id0,object_id0] corresponds to distribution id
+   * distribution_id0
+   */
+  int first = 0;
+  int last = kernel_data.integrator.num_triangle_lights;
+  int count = last - first;
+  int middle, step;
+  while (count > 0) {
+    step = count / 2;
+    middle = first + step;
+    int triangle = kernel_tex_fetch(__triangle_to_distribution, middle * 3);
+    if (triangle < triangle_id) {
+      first = middle + 1;
+      count -= step + 1;
+    }
+    else
+      count = step;
+  }
+
+  /* If instancing then we can have several triangles with the same triangle_id
+   * so loop over object_id too. */
+  /* todo: do a binary search here too if many instances */
+  while (true) {
+    int object = kernel_tex_fetch(__triangle_to_distribution, first * 3 + 1);
+    if (object == object_id)
+      break;
+    ++first;
+  }
+
+  kernel_assert(kernel_tex_fetch(__triangle_to_distribution, first * 3) == triangle_id);
+
+  return kernel_tex_fetch(__triangle_to_distribution, first * 3 + 2);
+}
+
+/* Decides whether to go down both childen or only one in the tree traversal.
+ * The split heuristic is based on the variance of the lighting within the node.
+ * There are two types of variances that are considered: variance in energy and
+ * in the distance term 1/d^2. The variance in energy is pre-computed on the
+ * host but the distance term is calculated here. These variances are then
+ * combined and normalized to get the final splitting heuristic. High variance
+ * leads to a lower splitting heuristic which leads to more splits during the
+ * traversal. */
+ccl_device bool split(KernelGlobals *kg, float3 P, int node_offset)
+{
+  /* early exists if never/always splitting */
+  const float threshold = kernel_data.integrator.splitting_threshold;
+  if (threshold == 0.0f) {
+    return false;
+  }
+  else if (threshold == 1.0f) {
+    return true;
+  }
+
+  /* extract bounding box of cluster */
+  const float4 node1 = kernel_tex_fetch(__light_tree_nodes, node_offset + 1);
+  const float4 node2 = kernel_tex_fetch(__light_tree_nodes, node_offset + 2);
+  const float3 bboxMin = make_float3(node1.x, node1.y, node1.z);
+  const float3 bboxMax = make_float3(node1.w, node2.x, node2.y);
+
+  /* if P is inside bounding sphere then split */
+  const float3 centroid = 0.5f * (bboxMax + bboxMin);
+  const float radius_squared = len_squared(bboxMax - centroid);
+  const float dist_squared = len_squared(centroid - P);
+
+  if (dist_squared <= radius_squared) {
+    return true;
+  }
+
+  /* eq. 8 & 9 */
+
+  /* the integral in eq. 8 requires us to know the interval the distance can
+   * be in: [a,b]. This is found by considering a bounding sphere around the
+   * bounding box of the node and "a" then becomes the smallest distance to
+   * this sphere and "b" becomes the largest. */
+  const float radius = sqrt(radius_squared);
+  const float dist = sqrt(dist_squared);
+  const float a = dist - radius;
+  const float b = dist + radius;
+
+  const float g_mean = 1.0f / (a * b);
+  const float g_mean_squared = g_mean * g_mean;
+  const float a3 = a * a * a;
+  const float b3 = b * b * b;
+  const float g_variance = (b3 - a3) / (3.0f * (b - a) * a3 * b3) - g_mean_squared;
+
+  /* eq. 10 */
+  const float4 node0 = kernel_tex_fetch(__light_tree_nodes, node_offset);
+  const float4 node3 = kernel_tex_fetch(__light_tree_nodes, node_offset + 3);
+  const float energy = node0.x;
+  const float e_variance = node3.w;
+  const float num_emitters = (float)__float_as_int(node0.w);
+  const float num_emitters_squared = num_emitters * num_emitters;
+  const float e_mean = energy / num_emitters;
+  const float e_mean_squared = e_mean * e_mean;
+  const float variance = (e_variance * (g_variance + g_mean_squared) +
+                          e_mean_squared * g_variance) *
+                         num_emitters_squared;
+
+  /* normalize the variance heuristic to be within [0,1]. Note that high
+   * variance corresponds to a low normalized variance. To give an idea of
+   * how this normalization function looks like:
+   *            variance: 0    1   10  100  1000  10000  100000
+   * normalized variance: 1  0.8  0.7  0.5   0.4    0.3     0.2  */
+  const float variance_normalized = sqrt(sqrt(1.0f / (1.0f + sqrt(variance))));
+
+  return variance_normalized < threshold;
+}
+
+/* given a light in the form of a distribution id, this function computes the
+ * the probability of picking it using the light tree. this mimics the
+ * probability calculations in accum_light_tree_contribution()
+ *
+ * the nodes array contains all the nodes of the tree and each interior node
+ * has its left child directly after it in the nodes array and the right child
+ * is also after it at the second_child_offset.
+ *
+ * Given the structure of the nodes array we can find the path from the root
+ * to the leaf node the given light belongs to as follows:
+ *
+ * 1. Find the offset of the leaf node the given light belongs to.
+ * 2. Traverse the tree in a top-down manner where the decision to go down the
+ *    left or right child is determined as follows:
+ *    If the node we are looking for has a lower offset than the right child
+ *    then it belongs to a node between the current node and the right child.
+ *    That is, we should go down the left node. Otherwise, the right node.
+ *    This is done recursively until the leaf node is found.
+ * 3. Before going down the left or right child:
+ *    a) If we are splitting then the probability is not affected.
+ *    b) If we are not splitting then the probability is multiplied by the
+ *       probability of choosing this particular child node.
+ */
+ccl_device float light_tree_pdf(KernelGlobals *kg,
+                                float3 P,
+                                float3 V,
+                                float t_max,
+                                int distribution_id,
+                                int offset,
+                                float pdf,
+                                bool can_split)
+{
+
+  /* find mapping from distribution_id to node_id */
+  int node_id = kernel_tex_fetch(__light_distribution_to_node, distribution_id);
+
+  /* read in first part of node of light tree */
+  int right_child_offset, first_distribution_id, num_emitters;
+  update_node(kg, offset, &right_child_offset, &first_distribution_id, &num_emitters);
+
+  while (right_child_offset != -1) {
+    int child_offsetL = offset + 4;
+    int child_offsetR = 4 * right_child_offset;
+
+    /* choose whether to go down both(split) or only one of the children */
+    if (can_split && split(kg, P, offset)) {
+      /* go down to the child node that is an ancestor of this node_id
+       * without changing the probability since we split here */
+
+      if (node_id < child_offsetR) {
+        offset = child_offsetL;
+      }
+      else {
+        offset = child_offsetR;
+      }
+
+      return light_tree_pdf(kg, P, V, t_max, distribution_id, offset, pdf, true);
+    }
+    else {
+      /* go down one of the child nodes */
+
+      /* evaluate the importance of each of the child nodes */
+      float I_L = calc_node_importance(kg, P, V, t_max, child_offsetL);
+      float I_R = calc_node_importance(kg, P, V, t_max, child_offsetR);
+
+      if ((I_L == 0.0f) && (I_R == 0.0f)) {
+        return 0.0f;
+      }
+
+      /* calculate the probability of going down the left node */
+      float P_L = I_L / (I_L + I_R);
+
+      /* choose which node to go down */
+      if (node_id < child_offsetR) {
+        offset = child_offsetL;
+        pdf *= P_L;
+      }
+      else {
+        offset = child_offsetR;
+        pdf *= 1.0f - P_L;
+      }
+      update_node(kg, offset, &right_child_offset, &first_distribution_id, &num_emitters);
+    }
+  }
+
+  /* if there are several emitters in this leaf then pick one of them */
+  if (num_emitters > 1) {
+
+    /* the case of being a light inside a leaf node with several lights.
+     * during sampling, a CDF is created based on importance, so here
+     * the probability of sampling this light using the CDF has to be
+     * computed. This is done by dividing the importance of this light
+     * by the total sum of the importance of all lights in the leaf. */
+    float sum = 0.0f;
+    for (int i = 0; i < num_emitters; ++i) {
+      sum += calc_light_importance(kg, P, V, t_max, offset, i);
+    }
+
+    if (sum == 0.0f) {
+      return 0.0f;
+    }
+
+    pdf *= calc_light_importance(
+               kg, P, V, t_max, offset, distribution_id - first_distribution_id) /
+           sum;
+  }
+
+  return pdf;
+}
+
+/* computes the the probability of picking the given light out of all lights.
+ * this mimics the probability calculations in light_distribution_sample() */
+ccl_device float light_distribution_pdf(
+    KernelGlobals *kg, float3 P, float3 V, float t_max, int prim_id, int object_id)
+{
+  /* convert from triangle/lamp to light distribution */
+  int distribution_id;
+  if (prim_id >= 0) {  // Triangle_id = prim_id
+    distribution_id = triangle_to_distribution(kg, prim_id, object_id);
+  }
+  else {  // Lamp
+    int lamp_id = -prim_id - 1;
+    distribution_id = kernel_tex_fetch(__lamp_to_distribution, lamp_id);
+  }
+
+  kernel_assert((distribution_id >= 0) &&
+                (distribution_id < kernel_data.integrator.num_distribution));
+
+  /* compute picking pdf for this light */
+  if (kernel_data.integrator.use_light_tree) {
+    /* find out which group of lights to sample */
+    int group;
+    if (prim_id >= 0) {
+      group = LIGHTGROUP_TREE;
+    }
+    else {
+      int lamp = -prim_id - 1;
+      int light_type = kernel_tex_fetch(__lights, lamp).type;
+      if (light_type == LIGHT_DISTANT) {
+        group = LIGHTGROUP_DISTANT;
+      }
+      else if (light_type == LIGHT_BACKGROUND) {
+        group = LIGHTGROUP_BACKGROUND;
+      }
+      else {
+        group = LIGHTGROUP_TREE;
+      }
+    }
+
+    /* get probabilty to sample this group of lights */
+    float group_prob = kernel_tex_fetch(__light_group_sample_prob, group);
+    float pdf = group_prob;
+
+    if (group == LIGHTGROUP_TREE) {
+      pdf *= light_tree_pdf(kg, P, V, t_max, distribution_id, 0, 1.0f, true);
+    }
+    else if (group == LIGHTGROUP_DISTANT) {
+      pdf *= kernel_data.integrator.inv_num_distant_lights;
+    }
+    else if (group == LIGHTGROUP_BACKGROUND) {
+      /* there is only one background light so nothing to do here */
+    }
+    else {
+      kernel_assert(false);
+    }
+
+    return pdf;
+  }
+  else {
+    const ccl_global KernelLightDistribution *kdistribution = &kernel_tex_fetch(
+        __light_distribution, distribution_id);
+    return kdistribution->area * kernel_data.integrator.pdf_inv_totarea;
+  }
+}
+
+/* picks a light and returns its index and the probability of picking it */
+ccl_device void light_distribution_sample(
+    KernelGlobals *kg, float3 P, float3 V, float t_max, float *randu, int *index, float *pdf)
+{
+  if (kernel_data.integrator.use_light_tree) {
+    /* sample light type distribution */
+    int group = light_group_distribution_sample(kg, randu);
+    float group_prob = kernel_tex_fetch(__light_group_sample_prob, group);
+
+    if (group == LIGHTGROUP_TREE) {
+      light_tree_sample(kg, P, V, t_max, randu, index, pdf);
+    }
+    else if (group == LIGHTGROUP_DISTANT) {
+      light_distant_sample(kg, P, randu, index, pdf);
+    }
+    else if (group == LIGHTGROUP_BACKGROUND) {
+      light_background_sample(kg, P, randu, index, pdf);
+    }
+    else {
+      kernel_assert(false);
+    }
+
+    *pdf *= group_prob;
+  }
+  else {  // Sample light distribution CDF
+    *index = light_distribution_sample(kg, randu);
+    const ccl_global KernelLightDistribution *kdistribution = &kernel_tex_fetch(
+        __light_distribution, *index);
+    *pdf = kdistribution->area * kernel_data.integrator.pdf_inv_totarea;
+  }
+}
+
+/* picks a point on a given light and computes the probability of picking this point*/
+ccl_device void light_point_sample(KernelGlobals *kg,
+                                   int lamp,
+                                   float randu,
+                                   float randv,
+                                   float time,
+                                   float3 P,
+                                   int bounce,
+                                   int distribution_id,
+                                   LightSample *ls)
 {
   if (lamp < 0) {
-    /* sample index */
-    int index = light_distribution_sample(kg, &randu);
-
-    /* fetch light data */
+    /* fetch light data and compute rest of light pdf */
     const ccl_global KernelLightDistribution *kdistribution = &kernel_tex_fetch(
-        __light_distribution, index);
+        __light_distribution, distribution_id);
     int prim = kdistribution->prim;
 
     if (prim >= 0) {
@@ -647,17 +1336,54 @@ ccl_device_noinline bool light_sample(KernelGlobals *kg,
 
       triangle_light_sample(kg, prim, object, randu, randv, time, ls, P);
       ls->shader |= shader_flag;
-      return (ls->pdf > 0.0f);
+      return;
     }
-
     lamp = -prim - 1;
   }
 
   if (UNLIKELY(light_select_reached_max_bounces(kg, lamp, bounce))) {
+    ls->pdf = 0.0f;
+    return;
+  }
+
+  if (!lamp_light_sample(kg, lamp, randu, randv, P, ls)) {
+    ls->pdf = 0.0f;
+    return;
+  }
+}
+
+/* picks a light and then picks a point on the light and computes the
+ * probability of doing so. V and t_max are only used when using light tree for sampling. If t_max
+ * is negative V should be the normal vector at P. Otherwise V should be the direction of the light
+ * ray starting at P.
+ */
+ccl_device_noinline bool light_sample(KernelGlobals *kg,
+                                      int lamp,
+                                      float randu,
+                                      float randv,
+                                      float time,
+                                      float3 P,
+                                      float3 V,
+                                      float t_max,
+                                      int bounce,
+                                      LightSample *ls)
+{
+  /* pick a light and compute the probability of picking this light */
+  float pdf_factor = 0.0f;
+  int index = -1;
+  light_distribution_sample(kg, P, V, t_max, &randu, &index, &pdf_factor);
+
+  if (pdf_factor == 0.0f) {
     return false;
   }
 
-  return lamp_light_sample(kg, lamp, randu, randv, P, ls);
+  /* pick a point on the light and the probability of picking this point */
+  light_point_sample(kg, lamp, randu, randv, time, P, bounce, index, ls);
+
+  /* combine pdfs */
+  ls->pdf *= pdf_factor;
+
+  return (ls->pdf > 0.0f);
 }
 
 ccl_device_inline int light_select_num_samples(KernelGlobals *kg, int index)

intern/cycles/kernel/kernel_path.h

@@ -99,6 +99,7 @@ ccl_device_forceinline void kernel_path_lamp_emission(KernelGlobals *kg,
     Ray light_ray ccl_optional_struct_init;
 
     light_ray.P = ray->P - state->ray_t * ray->D;
+
     state->ray_t += isect->t;
     light_ray.D = ray->D;
     light_ray.t = state->ray_t;
@@ -185,6 +186,8 @@ ccl_device_forceinline VolumeIntegrateResult kernel_path_volume(KernelGlobals *k
     shader_setup_from_volume(kg, sd, &volume_ray);
     kernel_volume_decoupled_record(kg, state, &volume_ray, sd, &volume_segment, step_size);
 
+    kernel_update_light_picking(sd, &volume_ray);
+
     volume_segment.sampling_method = sampling_method;
 
     /* emission */
@@ -231,6 +234,8 @@ ccl_device_forceinline VolumeIntegrateResult kernel_path_volume(KernelGlobals *k
     VolumeIntegrateResult result = kernel_volume_integrate(
         kg, state, sd, &volume_ray, L, throughput, step_size);
 
+    kernel_update_light_picking(sd, NULL);
+
 #    ifdef __VOLUME_SCATTER__
     if (result == VOLUME_PATH_SCATTERED) {
       /* direct lighting */
@@ -390,6 +395,7 @@ ccl_device void kernel_path_indirect(KernelGlobals *kg,
 
     /* path iteration */
     for (;;) {
+
       /* Find intersection with objects in scene. */
       Intersection isect;
       bool hit = kernel_path_scene_intersect(kg, state, ray, &isect, L);
@@ -408,6 +414,7 @@ ccl_device void kernel_path_indirect(KernelGlobals *kg,
       else if (result == VOLUME_PATH_MISSED) {
         break;
       }
+
 #    endif /* __VOLUME__*/
 
       /* Shade background. */
@@ -453,6 +460,8 @@ ccl_device void kernel_path_indirect(KernelGlobals *kg,
           throughput /= probability;
         }
 
+        kernel_update_light_picking(sd, NULL);
+
 #    ifdef __DENOISING_FEATURES__
         kernel_update_denoising_features(kg, sd, state, L);
 #    endif
@@ -518,7 +527,6 @@ ccl_device_forceinline void kernel_path_integrate(KernelGlobals *kg,
 
   /* Shader data memory used for both volumes and surfaces, saves stack space. */
   ShaderData sd;
-
 #  ifdef __SUBSURFACE__
   SubsurfaceIndirectRays ss_indirect;
   kernel_path_subsurface_init_indirect(&ss_indirect);
@@ -528,6 +536,7 @@ ccl_device_forceinline void kernel_path_integrate(KernelGlobals *kg,
 
     /* path iteration */
     for (;;) {
+
       /* Find intersection with objects in scene. */
       Intersection isect;
       bool hit = kernel_path_scene_intersect(kg, state, ray, &isect, L);
@@ -590,6 +599,8 @@ ccl_device_forceinline void kernel_path_integrate(KernelGlobals *kg,
           throughput /= probability;
         }
 
+        kernel_update_light_picking(&sd, NULL);
+
 #  ifdef __DENOISING_FEATURES__
         kernel_update_denoising_features(kg, &sd, state, L);
 #  endif

intern/cycles/kernel/kernel_path_branched.h

@@ -102,6 +102,8 @@ ccl_device_forceinline void kernel_branched_path_volume(KernelGlobals *kg,
     shader_setup_from_volume(kg, sd, &volume_ray);
     kernel_volume_decoupled_record(kg, state, &volume_ray, sd, &volume_segment, step_size);
 
+    kernel_update_light_picking(sd, &volume_ray);
+
     /* direct light sampling */
     if (volume_segment.closure_flag & SD_SCATTER) {
       volume_segment.sampling_method = volume_stack_sampling_method(kg, state->volume_stack);
@@ -134,6 +136,9 @@ ccl_device_forceinline void kernel_branched_path_volume(KernelGlobals *kg,
 
         if (result == VOLUME_PATH_SCATTERED &&
             kernel_path_volume_bounce(kg, sd, &tp, &ps, &L->state, &pray)) {
+          indirect_sd->P_pick = sd->P_pick;
+          indirect_sd->V_pick = sd->V_pick;
+          indirect_sd->t_pick = sd->t_pick;
           kernel_path_indirect(kg, indirect_sd, emission_sd, &pray, tp * num_samples_inv, &ps, L);
 
           /* for render passes, sum and reset indirect light pass variables
@@ -173,6 +178,8 @@ ccl_device_forceinline void kernel_branched_path_volume(KernelGlobals *kg,
       VolumeIntegrateResult result = kernel_volume_integrate(
           kg, &ps, sd, &volume_ray, L, &tp, step_size);
 
+      kernel_update_light_picking(sd, &volume_ray);
+
 #      ifdef __VOLUME_SCATTER__
       if (result == VOLUME_PATH_SCATTERED) {
         /* todo: support equiangular, MIS and all light sampling.
@@ -180,6 +187,9 @@ ccl_device_forceinline void kernel_branched_path_volume(KernelGlobals *kg,
         kernel_path_volume_connect_light(kg, sd, emission_sd, tp, state, L);
 
         if (kernel_path_volume_bounce(kg, sd, &tp, &ps, &L->state, &pray)) {
+          indirect_sd->P_pick = sd->P_pick;
+          indirect_sd->V_pick = sd->V_pick;
+          indirect_sd->t_pick = sd->t_pick;
           kernel_path_indirect(kg, indirect_sd, emission_sd, &pray, tp, &ps, L);
 
           /* for render passes, sum and reset indirect light pass variables
@@ -265,7 +275,9 @@ ccl_device_noinline_cpu void kernel_branched_path_surface_indirect_light(KernelG
       }
 
       ps.rng_hash = state->rng_hash;
-
+      indirect_sd->P_pick = sd->P_pick;
+      indirect_sd->V_pick = sd->V_pick;
+      indirect_sd->t_pick = sd->t_pick;
       kernel_path_indirect(kg, indirect_sd, emission_sd, &bsdf_ray, tp * num_samples_inv, &ps, L);
 
       /* for render passes, sum and reset indirect light pass variables
@@ -393,6 +405,7 @@ ccl_device void kernel_branched_path_integrate(KernelGlobals *kg,
    * Indirect bounces are handled in kernel_branched_path_surface_indirect_light().
    */
   for (;;) {
+
     /* Find intersection with objects in scene. */
     Intersection isect;
     bool hit = kernel_path_scene_intersect(kg, &state, &ray, &isect, L);
@@ -445,6 +458,8 @@ ccl_device void kernel_branched_path_integrate(KernelGlobals *kg,
         }
       }
 
+      kernel_update_light_picking(&sd, NULL);
+
 #    ifdef __DENOISING_FEATURES__
       kernel_update_denoising_features(kg, &sd, &state, L);
 #    endif

intern/cycles/kernel/kernel_path_subsurface.h

@@ -68,6 +68,7 @@ ccl_device_inline
        * integration loop stops when this function returns true.
        */
       subsurface_scatter_multi_setup(kg, &ss_isect, hit, sd, state, bssrdf_type, bssrdf_roughness);
+      kernel_update_light_picking(sd, NULL);
 
       kernel_path_surface_connect_light(kg, sd, emission_sd, *throughput, state, L);
 

intern/cycles/kernel/kernel_path_surface.h

@@ -16,6 +16,220 @@
 
 CCL_NAMESPACE_BEGIN
 
+/* connect the given light sample with the shading point and calculate its
+ * contribution and accumulate it to L */
+ccl_device void accum_light_contribution(KernelGlobals *kg,
+                                         ShaderData *sd,
+                                         ShaderData *emission_sd,
+                                         LightSample *ls,
+                                         ccl_addr_space PathState *state,
+                                         Ray *light_ray,
+                                         BsdfEval *L_light,
+                                         PathRadiance *L,
+                                         bool *is_lamp,
+                                         float terminate,
+                                         float3 throughput,
+                                         float scale)
+{
+  if (direct_emission(kg, sd, emission_sd, ls, state, light_ray, L_light, is_lamp, terminate)) {
+    /* trace shadow ray */
+    float3 shadow;
+
+    if (!shadow_blocked(kg, sd, emission_sd, state, light_ray, &shadow)) {
+      /* accumulate */
+      path_radiance_accum_light(
+          kg, L, state, throughput * scale, L_light, shadow, scale, *is_lamp);
+    }
+    else {
+      path_radiance_accum_total_light(L, state, throughput * scale, L_light);
+    }
+  }
+}
+
+/* The accum_light_tree_contribution() function does the following:
+ * 1. Recursive tree traversal using splitting. This picks one or more lights.
+ * 2. For each picked light, a position on the light is also chosen.
+ * 3. The total contribution of all these light samples are evaluated and
+ *    accumulated to L. */
+ccl_device void accum_light_tree_contribution(KernelGlobals *kg,
+                                              float randu,
+                                              float randv,
+                                              int offset,
+                                              float pdf_factor,
+                                              bool can_split,
+                                              float3 throughput,
+                                              float scale_factor,
+                                              PathRadiance *L,
+                                              ccl_addr_space PathState *state,
+                                              ShaderData *sd,
+                                              ShaderData *emission_sd)
+{
+  float3 P = sd->P_pick;
+  float3 V = sd->V_pick;
+  float t = sd->t_pick;
+
+  float time = sd->time;
+  int bounce = state->bounce;
+
+  float randu_stack[64];
+  float randv_stack[64];
+  int offset_stack[64];
+  float pdf_stack[64];
+
+  randu_stack[0] = randu;
+  randv_stack[0] = randv;
+  offset_stack[0] = offset;
+  pdf_stack[0] = pdf_factor;
+
+  int stack_idx = 0;
+
+  while (stack_idx > -1) {
+    randu = randu_stack[stack_idx];
+    randv = randv_stack[stack_idx];
+    offset = offset_stack[stack_idx];
+    pdf_factor = pdf_stack[stack_idx];
+    /* read in first part of node of light tree */
+    int right_child_offset, distribution_id, num_emitters;
+    update_node(kg, offset, &right_child_offset, &distribution_id, &num_emitters);
+
+    /* found a leaf */
+    if (right_child_offset == -1) {
+
+      /* if there are several emitters in this leaf then pick one of them */
+      if (num_emitters > 1) {
+
+        /* create and sample CDF without dynamic allocation.
+         * see comment in light_tree_sample() for this piece of code */
+        float sum = 0.0f;
+        for (int i = 0; i < num_emitters; ++i) {
+          sum += calc_light_importance(kg, P, V, t, offset, i);
+        }
+
+        if (sum == 0.0f) {
+          --stack_idx;
+          continue;
+        }
+
+        float sum_inv = 1.0f / sum;
+        float cdf_L = 0.0f;
+        float cdf_R = 0.0f;
+        float prob = 0.0f;
+        int light = num_emitters - 1;
+        for (int i = 1; i < num_emitters + 1; ++i) {
+          prob = calc_light_importance(kg, P, V, t, offset, i - 1) * sum_inv;
+          cdf_R = cdf_L + prob;
+          if (randu < cdf_R) {
+            light = i - 1;
+            break;
+          }
+
+          cdf_L = cdf_R;
+        }
+        distribution_id += light;
+        pdf_factor *= prob;
+
+        /* rescale random number */
+        randu = (randu - cdf_L) / (cdf_R - cdf_L);
+      }
+
+      /* pick a point on the chosen light(distribution_id) and calculate the
+       * probability of picking this point */
+      LightSample ls;
+      light_point_sample(kg, -1, randu, randv, time, P, bounce, distribution_id, &ls);
+
+      /* combine pdfs */
+      ls.pdf *= pdf_factor;
+
+      if (ls.pdf <= 0.0f) {
+        --stack_idx;
+        continue;
+      }
+
+      /* compute and accumulate the total contribution of this light */
+      Ray light_ray;
+      light_ray.t = 0.0f;
+#ifdef __OBJECT_MOTION__
+      light_ray.time = sd->time;
+#endif
+      BsdfEval L_light;
+      bool is_lamp;
+      float terminate = path_state_rng_light_termination(kg, state);
+      accum_light_contribution(kg,
+                               sd,
+                               emission_sd,
+                               &ls,
+                               state,
+                               &light_ray,
+                               &L_light,
+                               L,
+                               &is_lamp,
+                               terminate,
+                               throughput,
+                               scale_factor);
+
+      --stack_idx;
+      can_split = true;
+      continue;
+    }
+    else {  // Interior node, choose which child(ren) to go down
+
+      int child_offsetL = offset + 4;
+      int child_offsetR = 4 * right_child_offset;
+
+      /* choose whether to go down both(split) or only one of the children */
+      if (can_split && split(kg, P, offset)) {
+        /* go down both child nodes */
+        randu_stack[stack_idx] = randu;
+        randv_stack[stack_idx] = randv;
+        offset_stack[stack_idx] = child_offsetL;
+        pdf_stack[stack_idx] = pdf_factor;
+
+        ++stack_idx;
+        randu_stack[stack_idx] = randu;
+        randv_stack[stack_idx] = randv;
+        offset_stack[stack_idx] = child_offsetR;
+        pdf_stack[stack_idx] = pdf_factor;
+      }
+      else {
+        /* go down one of the child nodes */
+
+        /* evaluate the importance of each of the child nodes */
+        float I_L = calc_node_importance(kg, P, V, t, child_offsetL);
+        float I_R = calc_node_importance(kg, P, V, t, child_offsetR);
+
+        if ((I_L == 0.0f) && (I_R == 0.0f)) {
+          return;
+        }
+
+        /* calculate the probability of going down the left node */
+        float P_L = I_L / (I_L + I_R);
+
+        /* choose which node to go down */
+        if (randu <= P_L) {  // Going down left node
+          /* rescale random number */
+          randu = randu / P_L;
+
+          offset = child_offsetL;
+          pdf_factor *= P_L;
+        }
+        else {  // Going down right node
+          /* rescale random number */
+          randu = (randu * (I_L + I_R) - I_L) / I_R;
+
+          offset = child_offsetR;
+          pdf_factor *= 1.0f - P_L;
+        }
+
+        can_split = false;
+        randu_stack[stack_idx] = randu;
+        randv_stack[stack_idx] = randv;
+        offset_stack[stack_idx] = offset;
+        pdf_stack[stack_idx] = pdf_factor;
+      }
+    }
+  }
+}
+
 #if defined(__BRANCHED_PATH__) || defined(__SUBSURFACE__) || defined(__SHADOW_TRICKS__) || \
     defined(__BAKING__)
 /* branched path tracing: connect path directly to position on one or more lights and add it to L
@@ -34,104 +248,174 @@ ccl_device_noinline_cpu void kernel_branched_path_surface_connect_light(
   /* sample illumination from lights to find path contribution */
   BsdfEval L_light ccl_optional_struct_init;
 
-  int num_lights = 0;
-  if (kernel_data.integrator.use_direct_light) {
-    if (sample_all_lights) {
-      num_lights = kernel_data.integrator.num_all_lights;
-      if (kernel_data.integrator.pdf_triangles != 0.0f) {
-        num_lights += 1;
-      }
+  bool use_light_tree = kernel_data.integrator.use_light_tree;
+  if (use_light_tree) {
+    Ray light_ray;
+    bool is_lamp;
+
+    light_ray.t = 0.0f;
+#    ifdef __OBJECT_MOTION__
+    light_ray.time = sd->time;
+#    endif
+
+    int index;
+    float randu, randv;
+    path_state_rng_2D(kg, state, PRNG_LIGHT_U, &randu, &randv);
+
+    /* sample light group distribution */
+    int group = light_group_distribution_sample(kg, &randu);
+    float group_prob = kernel_tex_fetch(__light_group_sample_prob, group);
+    float pdf = 1.0f;
+
+    if (group == LIGHTGROUP_TREE) {
+      /* accumulate contribution to L from potentially several lights */
+      accum_light_tree_contribution(kg,
+                                    randu,
+                                    randv,
+                                    0,
+                                    group_prob,
+                                    true,
+                                    throughput,
+                                    num_samples_adjust,
+                                    L,  // todo: is num_samples_adjust correct here?
+                                    state,
+                                    sd,
+                                    emission_sd);
+
+      /* have accumulated all the contributions so return */
+      return;
+    }
+    else if (group == LIGHTGROUP_DISTANT) {
+      /* pick a single distant light */
+      light_distant_sample(kg, sd->P, &randu, &index, &pdf);
+    }
+    else if (group == LIGHTGROUP_BACKGROUND) {
+      /* pick a single background light */
+      light_background_sample(kg, sd->P, &randu, &index, &pdf);
     }
     else {
-      num_lights = 1;
+      kernel_assert(false);
     }
+
+    /* sample a point on the given distant/background light */
+    LightSample ls;
+    light_point_sample(kg, -1, randu, randv, sd->time, sd->P, state->bounce, index, &ls);
+
+    /* combine pdfs */
+    ls.pdf *= group_prob;
+
+    if (ls.pdf <= 0.0f)
+      return;
+
+    /* accumulate the contribution of this distant/background light to L */
+    float terminate = path_state_rng_light_termination(kg, state);
+    accum_light_contribution(kg,
+                             sd,
+                             emission_sd,
+                             &ls,
+                             state,
+                             &light_ray,
+                             &L_light,
+                             L,
+                             &is_lamp,
+                             terminate,
+                             throughput,
+                             num_samples_adjust);
   }
-
-  for (int i = 0; i < num_lights; i++) {
-    /* sample one light at random */
-    int num_samples = 1;
-    int num_all_lights = 1;
-    uint lamp_rng_hash = state->rng_hash;
-    bool double_pdf = false;
-    bool is_mesh_light = false;
-    bool is_lamp = false;
-
-    if (sample_all_lights) {
-      /* lamp sampling */
-      is_lamp = i < kernel_data.integrator.num_all_lights;
-      if (is_lamp) {
-        if (UNLIKELY(light_select_reached_max_bounces(kg, i, state->bounce))) {
-          continue;
+  else {
+    int num_lights = 0;
+    if (kernel_data.integrator.use_direct_light) {
+      if (sample_all_lights) {
+        num_lights = kernel_data.integrator.num_all_lights;
+        if (kernel_data.integrator.pdf_triangles != 0.0f) {
+          num_lights += 1;
         }
-        num_samples = ceil_to_int(num_samples_adjust * light_select_num_samples(kg, i));
-        num_all_lights = kernel_data.integrator.num_all_lights;
-        lamp_rng_hash = cmj_hash(state->rng_hash, i);
-        double_pdf = kernel_data.integrator.pdf_triangles != 0.0f;
       }
-      /* mesh light sampling */
       else {
-        num_samples = ceil_to_int(num_samples_adjust * kernel_data.integrator.mesh_light_samples);
-        double_pdf = kernel_data.integrator.num_all_lights != 0;
-        is_mesh_light = true;
+        num_lights = 1;
       }
     }
 
-    float num_samples_inv = num_samples_adjust / (num_samples * num_all_lights);
+    for (int i = 0; i < num_lights; i++) {
+      /* sample one light at random */
+      int num_samples = 1;
+      uint lamp_rng_hash = state->rng_hash;
+      bool double_pdf = false;
+      bool is_mesh_light = false;
+      bool is_lamp = false;
 
-    for (int j = 0; j < num_samples; j++) {
-      Ray light_ray ccl_optional_struct_init;
-      light_ray.t = 0.0f; /* reset ray */
-#    ifdef __OBJECT_MOTION__
-      light_ray.time = sd->time;
-#    endif
-      bool has_emission = false;
-
-      if (kernel_data.integrator.use_direct_light && (sd->flag & SD_BSDF_HAS_EVAL)) {
-        float light_u, light_v;
-        path_branched_rng_2D(
-            kg, lamp_rng_hash, state, j, num_samples, PRNG_LIGHT_U, &light_u, &light_v);
-        float terminate = path_branched_rng_light_termination(
-            kg, lamp_rng_hash, state, j, num_samples);
-
-        /* only sample triangle lights */
-        if (is_mesh_light && double_pdf) {
-          light_u = 0.5f * light_u;
+      if (sample_all_lights) {
+        /* lamp sampling */
+        is_lamp = i < kernel_data.integrator.num_all_lights;
+        if (is_lamp) {
+          if (UNLIKELY(light_select_reached_max_bounces(kg, i, state->bounce))) {
+            continue;
+          }
+          num_samples = ceil_to_int(num_samples_adjust * light_select_num_samples(kg, i));
+          lamp_rng_hash = cmj_hash(state->rng_hash, i);
+          double_pdf = kernel_data.integrator.pdf_triangles != 0.0f;
         }
+        /* mesh light sampling */
+        else {
+          num_samples = ceil_to_int(num_samples_adjust *
+                                    kernel_data.integrator.mesh_light_samples);
+          double_pdf = kernel_data.integrator.num_all_lights != 0;
+          is_mesh_light = true;
+        }
+      }
 
-        LightSample ls ccl_optional_struct_init;
-        const int lamp = is_lamp ? i : -1;
-        if (light_sample(kg, lamp, light_u, light_v, sd->time, sd->P, state->bounce, &ls)) {
-          /* The sampling probability returned by lamp_light_sample assumes that all lights were
-           * sampled. However, this code only samples lamps, so if the scene also had mesh lights,
-           * the real probability is twice as high. */
-          if (double_pdf) {
-            ls.pdf *= 2.0f;
+      float num_samples_inv = num_samples_adjust / num_samples;
+
+      for (int j = 0; j < num_samples; j++) {
+        Ray light_ray ccl_optional_struct_init;
+        light_ray.t = 0.0f; /* reset ray */
+#    ifdef __OBJECT_MOTION__
+        light_ray.time = sd->time;
+#    endif
+
+        if (kernel_data.integrator.use_direct_light && (sd->flag & SD_BSDF_HAS_EVAL)) {
+          float light_u, light_v;
+          path_branched_rng_2D(
+              kg, lamp_rng_hash, state, j, num_samples, PRNG_LIGHT_U, &light_u, &light_v);
+          float terminate = path_branched_rng_light_termination(
+              kg, lamp_rng_hash, state, j, num_samples);
+
+          /* only sample triangle lights */
+          if (is_mesh_light && double_pdf) {
+            light_u = 0.5f * light_u;
           }
 
-          has_emission = direct_emission(
-              kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate);
-        }
-      }
-
-      /* trace shadow ray */
-      float3 shadow;
-
-      const bool blocked = shadow_blocked(kg, sd, emission_sd, state, &light_ray, &shadow);
-
-      if (has_emission) {
-        if (!blocked) {
-          /* accumulate */
-          path_radiance_accum_light(kg,
-                                    L,
-                                    state,
-                                    throughput * num_samples_inv,
-                                    &L_light,
-                                    shadow,
-                                    num_samples_inv,
-                                    is_lamp);
-        }
-        else {
-          path_radiance_accum_total_light(L, state, throughput * num_samples_inv, &L_light);
+          LightSample ls ccl_optional_struct_init;
+          const int lamp = is_lamp ? i : -1;
+          if (light_sample(kg,
+                           lamp,
+                           light_u,
+                           light_v,
+                           sd->time,
+                           sd->P_pick,
+                           sd->V_pick,
+                           sd->t_pick,
+                           state->bounce,
+                           &ls)) {
+            /* The sampling probability returned by lamp_light_sample assumes that all lights were
+             * sampled. However, this code only samples lamps, so if the scene also had mesh
+             * lights, the real probability is twice as high. */
+            if (double_pdf) {
+              ls.pdf *= 2.0f;
+            }
+            accum_light_contribution(kg,
+                                     sd,
+                                     emission_sd,
+                                     &ls,
+                                     state,
+                                     &light_ray,
+                                     &L_light,
+                                     L,
+                                     &is_lamp,
+                                     terminate,
+                                     throughput,
+                                     num_samples_inv);
+          }
         }
       }
     }
@@ -225,11 +509,11 @@ ccl_device_inline void kernel_path_surface_connect_light(KernelGlobals *kg,
   int all = (state->flag & PATH_RAY_SHADOW_CATCHER);
   kernel_branched_path_surface_connect_light(kg, sd, emission_sd, state, throughput, 1.0f, L, all);
 #  else
+
   /* sample illumination from lights to find path contribution */
   Ray light_ray ccl_optional_struct_init;
   BsdfEval L_light ccl_optional_struct_init;
   bool is_lamp = false;
-  bool has_emission = false;
 
   light_ray.t = 0.0f;
 #    ifdef __OBJECT_MOTION__
@@ -241,27 +525,32 @@ ccl_device_inline void kernel_path_surface_connect_light(KernelGlobals *kg,
     path_state_rng_2D(kg, state, PRNG_LIGHT_U, &light_u, &light_v);
 
     LightSample ls ccl_optional_struct_init;
-    if (light_sample(kg, -1, light_u, light_v, sd->time, sd->P, state->bounce, &ls)) {
+    if (light_sample(kg,
+                     -1,
+                     light_u,
+                     light_v,
+                     sd->time,
+                     sd->P_pick,
+                     sd->V_pick,
+                     sd->t_pick,
+                     state->bounce,
+                     &ls)) {
       float terminate = path_state_rng_light_termination(kg, state);
-      has_emission = direct_emission(
-          kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate);
+      accum_light_contribution(kg,
+                               sd,
+                               emission_sd,
+                               &ls,
+                               state,
+                               &light_ray,
+                               &L_light,
+                               L,
+                               &is_lamp,
+                               terminate,
+                               throughput,
+                               1.0f);
     }
   }
 
-  /* trace shadow ray */
-  float3 shadow;
-
-  const bool blocked = shadow_blocked(kg, sd, emission_sd, state, &light_ray, &shadow);
-
-  if (has_emission) {
-    if (!blocked) {
-      /* accumulate */
-      path_radiance_accum_light(kg, L, state, throughput, &L_light, shadow, 1.0f, is_lamp);
-    }
-    else {
-      path_radiance_accum_total_light(L, state, throughput, &L_light);
-    }
-  }
 #  endif
 #endif
 }
@@ -310,6 +599,7 @@ ccl_device bool kernel_path_surface_bounce(KernelGlobals *kg,
 
     /* setup ray */
     ray->P = ray_offset(sd->P, (label & LABEL_TRANSMIT) ? -sd->Ng : sd->Ng);
+    kernel_update_light_picking(sd, NULL);
     ray->D = normalize(bsdf_omega_in);
 
     if (state->bounce == 0)
@@ -342,6 +632,8 @@ ccl_device bool kernel_path_surface_bounce(KernelGlobals *kg,
 
     /* setup ray position, direction stays unchanged */
     ray->P = ray_offset(sd->P, -sd->Ng);
+    kernel_update_light_picking(sd, NULL);
+
 #  ifdef __RAY_DIFFERENTIALS__
     ray->dP = sd->dP;
 #  endif

intern/cycles/kernel/kernel_path_volume.h

@@ -43,7 +43,16 @@ ccl_device_inline void kernel_path_volume_connect_light(KernelGlobals *kg,
     path_state_rng_2D(kg, state, PRNG_LIGHT_U, &light_u, &light_v);
 
     LightSample ls ccl_optional_struct_init;
-    if (light_sample(kg, -1, light_u, light_v, sd->time, sd->P, state->bounce, &ls)) {
+    if (light_sample(kg,
+                     -1,
+                     light_u,
+                     light_v,
+                     sd->time,
+                     sd->P_pick,
+                     sd->V_pick,
+                     sd->t_pick,
+                     state->bounce,
+                     &ls)) {
       float terminate = path_state_rng_light_termination(kg, state);
       has_emission = direct_emission(
           kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate);
@@ -128,6 +137,210 @@ ccl_device_noinline_cpu bool kernel_path_volume_bounce(KernelGlobals *kg,
 }
 
 #  if !defined(__SPLIT_KERNEL__) && (defined(__BRANCHED_PATH__) || defined(__VOLUME_DECOUPLED__))
+ccl_device void accum_light_tree_contribution_volume(KernelGlobals *kg,
+                                                     float randu,
+                                                     float randv,
+                                                     int offset,
+                                                     float pdf_factor,
+                                                     bool can_split,
+                                                     float3 throughput,
+                                                     float scale_factor,
+                                                     PathRadiance *L,
+                                                     ccl_addr_space PathState *state,
+                                                     ShaderData *sd,
+                                                     ShaderData *emission_sd,
+                                                     Ray *ray,
+                                                     const VolumeSegment *segment)
+{
+  float3 P = sd->P_pick;
+  float3 V = sd->V_pick;
+  float t = sd->t_pick;
+
+  float time = sd->time;
+  int bounce = state->bounce;
+
+  float randu_stack[64];
+  float randv_stack[64];
+  int offset_stack[64];
+  float pdf_stack[64];
+
+  randu_stack[0] = randu;
+  randv_stack[0] = randv;
+  offset_stack[0] = offset;
+  pdf_stack[0] = pdf_factor;
+
+  int stack_idx = 0;
+
+  while (stack_idx > -1) {
+    randu = randu_stack[stack_idx];
+    randv = randv_stack[stack_idx];
+    offset = offset_stack[stack_idx];
+    pdf_factor = pdf_stack[stack_idx];
+    /* read in first part of node of light tree */
+    int right_child_offset, distribution_id, num_emitters;
+    update_node(kg, offset, &right_child_offset, &distribution_id, &num_emitters);
+
+    /* found a leaf */
+    if (right_child_offset == -1) {
+
+      /* if there are several emitters in this leaf then pick one of them */
+      if (num_emitters > 1) {
+
+        /* create and sample CDF without dynamic allocation.
+         * see comment in light_tree_sample() for this piece of code */
+        float sum = 0.0f;
+        for (int i = 0; i < num_emitters; ++i) {
+          sum += calc_light_importance(kg, P, V, t, offset, i);
+        }
+
+        if (sum == 0.0f) {
+          --stack_idx;
+          continue;
+        }
+
+        float sum_inv = 1.0f / sum;
+        float cdf_L = 0.0f;
+        float cdf_R = 0.0f;
+        float prob = 0.0f;
+        int light = num_emitters - 1;
+        for (int i = 1; i < num_emitters + 1; ++i) {
+          prob = calc_light_importance(kg, P, V, t, offset, i - 1) * sum_inv;
+          cdf_R = cdf_L + prob;
+          if (randu < cdf_R) {
+            light = i - 1;
+            break;
+          }
+
+          cdf_L = cdf_R;
+        }
+        distribution_id += light;
+        pdf_factor *= prob;
+
+        /* rescale random number */
+        randu = (randu - cdf_L) / (cdf_R - cdf_L);
+      }
+
+      /* pick a point on the chosen light(distribution_id) and calculate the
+       * probability of picking this point */
+      LightSample ls;
+      light_point_sample(kg, -1, randu, randv, time, P, bounce, distribution_id, &ls);
+
+      /* combine pdfs */
+      ls.pdf *= pdf_factor;
+
+      /* compute and accumulate the total contribution of this light */
+      Ray light_ray;
+      light_ray.t = 0.0f;
+#    ifdef __OBJECT_MOTION__
+      light_ray.time = sd->time;
+#    endif
+      float3 tp = throughput;
+      bool has_emission = false;
+      bool is_lamp = false;
+      BsdfEval L_light ccl_optional_struct_init;
+      /* sample position on volume segment */
+      float rphase = path_branched_rng_1D(
+          kg, state->rng_hash, state, distribution_id, 1.0f, PRNG_PHASE_CHANNEL);
+      float rscatter = path_branched_rng_1D(
+          kg, state->rng_hash, state, distribution_id, 1.0f, PRNG_SCATTER_DISTANCE);
+
+      VolumeIntegrateResult result = kernel_volume_decoupled_scatter(kg,
+                                                                     state,
+                                                                     ray,
+                                                                     sd,
+                                                                     &tp,
+                                                                     rphase,
+                                                                     rscatter,
+                                                                     segment,
+                                                                     (ls.t != FLT_MAX) ? &ls.P :
+                                                                                         NULL,
+                                                                     false);
+      if (result == VOLUME_PATH_SCATTERED) {
+        light_point_sample(kg, -1, randu, randv, time, sd->P_pick, bounce, distribution_id, &ls);
+        if (ls.pdf <= 0.0f) {
+          --stack_idx;
+          continue;
+        }
+        /* sample random light */
+        float terminate = path_branched_rng_light_termination(
+            kg, state->rng_hash, state, distribution_id, 1.0f);
+        has_emission = direct_emission(
+            kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate);
+      }
+
+      /* trace shadow ray */
+      float3 shadow;
+
+      const bool blocked = shadow_blocked(kg, sd, emission_sd, state, &light_ray, &shadow);
+
+      if (has_emission && !blocked) {
+        /* accumulate */
+        path_radiance_accum_light(kg, L, state, tp, &L_light, shadow, 1.0f, is_lamp);
+      }
+
+      --stack_idx;
+      can_split = true;
+      continue;
+    }
+    else {  // Interior node, choose which child(ren) to go down
+
+      int child_offsetL = offset + 4;
+      int child_offsetR = 4 * right_child_offset;
+
+      /* choose whether to go down both(split) or only one of the children */
+      if (can_split && split(kg, P, offset)) {
+        /* go down both child nodes */
+        randu_stack[stack_idx] = randu;
+        randv_stack[stack_idx] = randv;
+        offset_stack[stack_idx] = child_offsetL;
+        pdf_stack[stack_idx] = pdf_factor;
+
+        ++stack_idx;
+        randu_stack[stack_idx] = randu;
+        randv_stack[stack_idx] = randv;
+        offset_stack[stack_idx] = child_offsetR;
+        pdf_stack[stack_idx] = pdf_factor;
+      }
+      else {
+        /* go down one of the child nodes */
+
+        /* evaluate the importance of each of the child nodes */
+        float I_L = calc_node_importance(kg, P, V, t, child_offsetL);
+        float I_R = calc_node_importance(kg, P, V, t, child_offsetR);
+
+        if ((I_L == 0.0f) && (I_R == 0.0f)) {
+          return;
+        }
+
+        /* calculate the probability of going down the left node */
+        float P_L = I_L / (I_L + I_R);
+
+        /* choose which node to go down */
+        if (randu <= P_L) {  // Going down left node
+          /* rescale random number */
+          randu = randu / P_L;
+
+          offset = child_offsetL;
+          pdf_factor *= P_L;
+        }
+        else {  // Going down right node
+          /* rescale random number */
+          randu = (randu * (I_L + I_R) - I_L) / I_R;
+
+          offset = child_offsetR;
+          pdf_factor *= 1.0f - P_L;
+        }
+
+        can_split = false;
+        randu_stack[stack_idx] = randu;
+        randv_stack[stack_idx] = randv;
+        offset_stack[stack_idx] = offset;
+        pdf_stack[stack_idx] = pdf_factor;
+      }
+    }
+  }
+}
+
 ccl_device void kernel_branched_path_volume_connect_light(KernelGlobals *kg,
                                                           ShaderData *sd,
                                                           ShaderData *emission_sd,
@@ -141,113 +354,228 @@ ccl_device void kernel_branched_path_volume_connect_light(KernelGlobals *kg,
 #    ifdef __EMISSION__
   BsdfEval L_light ccl_optional_struct_init;
 
-  int num_lights = 1;
-  if (sample_all_lights) {
-    num_lights = kernel_data.integrator.num_all_lights;
-    if (kernel_data.integrator.pdf_triangles != 0.0f) {
-      num_lights += 1;
+  bool use_light_tree = kernel_data.integrator.use_light_tree;
+  if (use_light_tree) {
+    Ray light_ray;
+    bool is_lamp;
+
+    light_ray.t = 0.0f;
+#      ifdef __OBJECT_MOTION__
+    light_ray.time = sd->time;
+#      endif
+
+    int index;
+    float randu, randv;
+    path_state_rng_2D(kg, state, PRNG_LIGHT_U, &randu, &randv);
+
+    /* sample light group distribution */
+    int group = light_group_distribution_sample(kg, &randu);
+    float group_prob = kernel_tex_fetch(__light_group_sample_prob, group);
+    float pdf = 1.0f;
+
+    if (group == LIGHTGROUP_TREE) {
+      /* accumulate contribution to L from potentially several lights */
+      accum_light_tree_contribution_volume(kg,
+                                           randu,
+                                           randv,
+                                           0,
+                                           group_prob,
+                                           true,
+                                           throughput,
+                                           1.0f,
+                                           L,  // todo: is num_samples_adjust correct here?
+                                           state,
+                                           sd,
+                                           emission_sd,
+                                           ray,
+                                           segment);
+
+      /* have accumulated all the contributions so return */
+      return;
+    }
+    else if (group == LIGHTGROUP_DISTANT) {
+      /* pick a single distant light */
+      light_distant_sample(kg, sd->P, &randu, &index, &pdf);
+    }
+    else if (group == LIGHTGROUP_BACKGROUND) {
+      /* pick a single background light */
+      light_background_sample(kg, sd->P, &randu, &index, &pdf);
+    }
+    else {
+      kernel_assert(false);
+    }
+
+    /* sample a point on the given distant/background light */
+    LightSample ls;
+    light_point_sample(kg, -1, randu, randv, sd->time, sd->P, state->bounce, index, &ls);
+
+    /* combine pdfs */
+    ls.pdf *= group_prob;
+
+    if (ls.pdf <= 0.0f)
+      return;
+
+    float3 tp = throughput;
+    bool has_emission = false;
+    /* sample position on volume segment */
+    float rphase = path_branched_rng_1D(
+        kg, state->rng_hash, state, index, 1.0f, PRNG_PHASE_CHANNEL);
+    float rscatter = path_branched_rng_1D(
+        kg, state->rng_hash, state, index, 1.0f, PRNG_SCATTER_DISTANCE);
+
+    VolumeIntegrateResult result = kernel_volume_decoupled_scatter(kg,
+                                                                   state,
+                                                                   ray,
+                                                                   sd,
+                                                                   &tp,
+                                                                   rphase,
+                                                                   rscatter,
+                                                                   segment,
+                                                                   (ls.t != FLT_MAX) ? &ls.P :
+                                                                                       NULL,
+                                                                   false);
+    if (result == VOLUME_PATH_SCATTERED) {
+      light_point_sample(kg, -1, randu, randv, sd->time, sd->P_pick, state->bounce, index, &ls);
+      /* sample random light */
+      float terminate = path_branched_rng_light_termination(
+          kg, state->rng_hash, state, index, 1.0f);
+      has_emission = direct_emission(
+          kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate);
+    }
+
+    /* trace shadow ray */
+    float3 shadow;
+
+    const bool blocked = shadow_blocked(kg, sd, emission_sd, state, &light_ray, &shadow);
+
+    if (has_emission && !blocked) {
+      /* accumulate */
+      path_radiance_accum_light(kg, L, state, tp, &L_light, shadow, 1.0f, is_lamp);
     }
   }
-
-  for (int i = 0; i < num_lights; ++i) {
-    /* sample one light at random */
-    int num_samples = 1;
-    int num_all_lights = 1;
-    uint lamp_rng_hash = state->rng_hash;
-    bool double_pdf = false;
-    bool is_mesh_light = false;
-    bool is_lamp = false;
-
+  else {
+    int num_lights = 1;
     if (sample_all_lights) {
-      /* lamp sampling */
-      is_lamp = i < kernel_data.integrator.num_all_lights;
-      if (is_lamp) {
-        if (UNLIKELY(light_select_reached_max_bounces(kg, i, state->bounce))) {
-          continue;
-        }
-        num_samples = light_select_num_samples(kg, i);
-        num_all_lights = kernel_data.integrator.num_all_lights;
-        lamp_rng_hash = cmj_hash(state->rng_hash, i);
-        double_pdf = kernel_data.integrator.pdf_triangles != 0.0f;
-      }
-      /* mesh light sampling */
-      else {
-        num_samples = kernel_data.integrator.mesh_light_samples;
-        double_pdf = kernel_data.integrator.num_all_lights != 0;
-        is_mesh_light = true;
+      num_lights = kernel_data.integrator.num_all_lights;
+      if (kernel_data.integrator.pdf_triangles != 0.0f) {
+        num_lights += 1;
       }
     }
+    for (int i = 0; i < num_lights; ++i) {
+      /* sample one light at random */
+      int num_samples = 1;
+      uint lamp_rng_hash = state->rng_hash;
+      bool double_pdf = false;
+      bool is_mesh_light = false;
+      bool is_lamp = false;
 
-    float num_samples_inv = 1.0f / (num_samples * num_all_lights);
-
-    for (int j = 0; j < num_samples; j++) {
-      Ray light_ray ccl_optional_struct_init;
-      light_ray.t = 0.0f; /* reset ray */
-#      ifdef __OBJECT_MOTION__
-      light_ray.time = sd->time;
-#      endif
-      bool has_emission = false;
-
-      float3 tp = throughput;
-
-      if (kernel_data.integrator.use_direct_light) {
-        /* sample random position on random light/triangle */
-        float light_u, light_v;
-        path_branched_rng_2D(
-            kg, lamp_rng_hash, state, j, num_samples, PRNG_LIGHT_U, &light_u, &light_v);
-
-        /* only sample triangle lights */
-        if (is_mesh_light && double_pdf) {
-          light_u = 0.5f * light_u;
-        }
-
-        LightSample ls ccl_optional_struct_init;
-        const int lamp = is_lamp ? i : -1;
-        light_sample(kg, lamp, light_u, light_v, sd->time, ray->P, state->bounce, &ls);
-
-        /* sample position on volume segment */
-        float rphase = path_branched_rng_1D(
-            kg, state->rng_hash, state, j, num_samples, PRNG_PHASE_CHANNEL);
-        float rscatter = path_branched_rng_1D(
-            kg, state->rng_hash, state, j, num_samples, PRNG_SCATTER_DISTANCE);
-
-        VolumeIntegrateResult result = kernel_volume_decoupled_scatter(kg,
-                                                                       state,
-                                                                       ray,
-                                                                       sd,
-                                                                       &tp,
-                                                                       rphase,
-                                                                       rscatter,
-                                                                       segment,
-                                                                       (ls.t != FLT_MAX) ? &ls.P :
-                                                                                           NULL,
-                                                                       false);
-
-        if (result == VOLUME_PATH_SCATTERED) {
-          /* todo: split up light_sample so we don't have to call it again with new position */
-          if (light_sample(kg, lamp, light_u, light_v, sd->time, sd->P, state->bounce, &ls)) {
-            if (double_pdf) {
-              ls.pdf *= 2.0f;
-            }
-
-            /* sample random light */
-            float terminate = path_branched_rng_light_termination(
-                kg, state->rng_hash, state, j, num_samples);
-            has_emission = direct_emission(
-                kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate);
+      if (sample_all_lights) {
+        /* lamp sampling */
+        is_lamp = i < kernel_data.integrator.num_all_lights;
+        if (is_lamp) {
+          if (UNLIKELY(light_select_reached_max_bounces(kg, i, state->bounce))) {
+            continue;
           }
+          lamp_rng_hash = cmj_hash(state->rng_hash, i);
+          double_pdf = kernel_data.integrator.pdf_triangles != 0.0f;
+        }
+        /* mesh light sampling */
+        else {
+          num_samples = kernel_data.integrator.mesh_light_samples;
+          double_pdf = kernel_data.integrator.num_all_lights != 0;
+          is_mesh_light = true;
         }
       }
 
-      /* trace shadow ray */
-      float3 shadow;
+      float num_samples_inv = 1.0f / num_samples;
 
-      const bool blocked = shadow_blocked(kg, sd, emission_sd, state, &light_ray, &shadow);
+      for (int j = 0; j < num_samples; j++) {
+        Ray light_ray ccl_optional_struct_init;
+        light_ray.t = 0.0f; /* reset ray */
+#      ifdef __OBJECT_MOTION__
+        light_ray.time = sd->time;
+#      endif
+        bool has_emission = false;
 
-      if (has_emission && !blocked) {
-        /* accumulate */
-        path_radiance_accum_light(
-            kg, L, state, tp * num_samples_inv, &L_light, shadow, num_samples_inv, is_lamp);
+        float3 tp = throughput;
+
+        if (kernel_data.integrator.use_direct_light) {
+          /* sample random position on random light/triangle */
+          float light_u, light_v;
+          path_branched_rng_2D(
+              kg, lamp_rng_hash, state, j, num_samples, PRNG_LIGHT_U, &light_u, &light_v);
+
+          /* only sample triangle lights */
+          if (is_mesh_light && double_pdf) {
+            light_u = 0.5f * light_u;
+          }
+
+          LightSample ls ccl_optional_struct_init;
+          const int lamp = is_lamp ? i : -1;
+          light_sample(kg,
+                       lamp,
+                       light_u,
+                       light_v,
+                       sd->time,
+                       sd->P_pick,
+                       sd->V_pick,
+                       sd->t_pick,
+                       state->bounce,
+                       &ls);
+
+          /* sample position on volume segment */
+          float rphase = path_branched_rng_1D(
+              kg, state->rng_hash, state, j, num_samples, PRNG_PHASE_CHANNEL);
+          float rscatter = path_branched_rng_1D(
+              kg, state->rng_hash, state, j, num_samples, PRNG_SCATTER_DISTANCE);
+
+          VolumeIntegrateResult result = kernel_volume_decoupled_scatter(
+              kg,
+              state,
+              ray,
+              sd,
+              &tp,
+              rphase,
+              rscatter,
+              segment,
+              (ls.t != FLT_MAX) ? &ls.P : NULL,
+              false);
+
+          if (result == VOLUME_PATH_SCATTERED) {
+            /* todo: split up light_sample so we don't have to call it again with new position */
+            if (light_sample(kg,
+                             lamp,
+                             light_u,
+                             light_v,
+                             sd->time,
+                             sd->P_pick,
+                             sd->V_pick,
+                             sd->t_pick,
+                             state->bounce,
+                             &ls)) {
+              if (double_pdf) {
+                ls.pdf *= 2.0f;
+              }
+
+              /* sample random light */
+              float terminate = path_branched_rng_light_termination(
+                  kg, state->rng_hash, state, j, num_samples);
+              has_emission = direct_emission(
+                  kg, sd, emission_sd, &ls, state, &light_ray, &L_light, &is_lamp, terminate);
+            }
+          }
+        }
+
+        /* trace shadow ray */
+        float3 shadow;
+
+        const bool blocked = shadow_blocked(kg, sd, emission_sd, state, &light_ray, &shadow);
+
+        if (has_emission && !blocked) {
+          /* accumulate */
+          path_radiance_accum_light(
+              kg, L, state, tp * num_samples_inv, &L_light, shadow, num_samples_inv, is_lamp);
+        }
       }
     }
   }

intern/cycles/kernel/kernel_textures.h

@@ -66,6 +66,14 @@ KERNEL_TEX(KernelLightDistribution, __light_distribution)
 KERNEL_TEX(KernelLight, __lights)
 KERNEL_TEX(float2, __light_background_marginal_cdf)
 KERNEL_TEX(float2, __light_background_conditional_cdf)
+KERNEL_TEX(float4, __light_tree_nodes)
+KERNEL_TEX(uint, __light_distribution_to_node)
+KERNEL_TEX(uint, __lamp_to_distribution)
+KERNEL_TEX(uint, __triangle_to_distribution)
+KERNEL_TEX(float, __light_group_sample_cdf)
+KERNEL_TEX(float, __light_group_sample_prob)
+KERNEL_TEX(float4, __light_tree_leaf_emitters)
+KERNEL_TEX(int, __leaf_to_first_emitter)
 
 /* particles */
 KERNEL_TEX(KernelParticle, __particles)

intern/cycles/kernel/kernel_types.h

@@ -626,6 +626,16 @@ enum PanoramaType {
   PANORAMA_NUM_TYPES,
 };
 
+/* Light Sampling Group */
+
+enum LightGroup {
+  LIGHTGROUP_TREE,
+  LIGHTGROUP_DISTANT,
+  LIGHTGROUP_BACKGROUND,
+
+  LIGHTGROUP_NUM,
+};
+
 /* Differential */
 
 typedef struct differential3 {
@@ -939,6 +949,12 @@ typedef ccl_addr_space struct ccl_align(16) ShaderData
   float3 N;
   /* true geometric normal */
   float3 Ng;
+  /* position used in light picking */
+  float3 P_pick;
+  /* normal or ray direction used in light picking */
+  float3 V_pick;
+  /* ray dist used for light picking */
+  float t_pick;
   /* view/incoming direction */
   float3 I;
   /* shader id */
@@ -1320,13 +1336,22 @@ static_assert_align(KernelBackground, 16);
 
 typedef struct KernelIntegrator {
   /* emission */
+  int use_light_tree;
+  float splitting_threshold;
   int use_direct_light;
   int use_ambient_occlusion;
   int num_distribution;
   int num_all_lights;
+  int num_light_nodes;
+  int num_triangle_lights;
+  int num_distant_lights;
+  float inv_num_distant_lights;
   float pdf_triangles;
   float pdf_lights;
+  float pdf_inv_totarea;
   float light_inv_rr_threshold;
+  int distant_lights_offset;
+  int background_light_index;
 
   /* bounces */
   int min_bounce;
@@ -1389,7 +1414,7 @@ typedef struct KernelIntegrator {
 
   int max_closures;
 
-  int pad1, pad2;
+  int pad1;
 } KernelIntegrator;
 static_assert_align(KernelIntegrator, 16);
 
@@ -1533,8 +1558,10 @@ typedef struct KernelLight {
 static_assert_align(KernelLight, 16);
 
 typedef struct KernelLightDistribution {
+  float area;
   float totarea;
   int prim;
+  float pad1, pad2, pad3;
   union {
     struct {
       int shader_flag;

intern/cycles/kernel/kernel_volume.h

@@ -1125,6 +1125,11 @@ ccl_device VolumeIntegrateResult kernel_volume_decoupled_scatter(KernelGlobals *
   /* move to new position */
   sd->P = ray->P + sample_t * ray->D;
 
+  kernel_update_light_picking(sd, NULL);
+
+  sd->V_pick = ray->D;
+  sd->t_pick = ray->t;
+
   return VOLUME_PATH_SCATTERED;
 }
 #  endif /* __SPLIT_KERNEL */

intern/cycles/kernel/split/kernel_direct_lighting.h

@@ -86,7 +86,16 @@ ccl_device void kernel_direct_lighting(KernelGlobals *kg,
       float terminate = path_state_rng_light_termination(kg, state);
 
       LightSample ls;
-      if (light_sample(kg, -1, light_u, light_v, sd->time, sd->P, state->bounce, &ls)) {
+      if (light_sample(kg,
+                       -1,
+                       light_u,
+                       light_v,
+                       sd->time,
+                       sd->P_pick,
+                       sd->V_pick,
+                       sd->t_pick,
+                       state->bounce,
+                       &ls)) {
         Ray light_ray;
         light_ray.time = sd->time;
 

intern/cycles/kernel/split/kernel_do_volume.h

@@ -63,6 +63,8 @@ ccl_device_noinline bool kernel_split_branched_path_volume_indirect_light_iter(K
     VolumeIntegrateResult result = kernel_volume_integrate(
         kg, ps, sd, &volume_ray, L, tp, step_size);
 
+    kernel_update_light_picking(sd, NULL);
+
 #  ifdef __VOLUME_SCATTER__
     if (result == VOLUME_PATH_SCATTERED) {
       /* direct lighting */

intern/cycles/kernel/split/kernel_lamp_emission.h

@@ -59,7 +59,6 @@ ccl_device void kernel_lamp_emission(KernelGlobals *kg)
     Ray ray = kernel_split_state.ray[ray_index];
     ccl_global Intersection *isect = &kernel_split_state.isect[ray_index];
     ShaderData *sd = kernel_split_sd(sd, ray_index);
-
     kernel_path_lamp_emission(kg, state, &ray, throughput, isect, sd, L);
   }
 }

intern/cycles/render/CMakeLists.txt

@@ -30,6 +30,7 @@ set(SRC
   integrator.cpp
   jitter.cpp
   light.cpp
+  light_tree.cpp
   merge.cpp
   mesh.cpp
   mesh_displace.cpp

intern/cycles/render/integrator.cpp

@@ -77,6 +77,11 @@ NODE_DEFINE(Integrator)
   SOCKET_BOOLEAN(sample_all_lights_direct, "Sample All Lights Direct", true);
   SOCKET_BOOLEAN(sample_all_lights_indirect, "Sample All Lights Indirect", true);
   SOCKET_FLOAT(light_sampling_threshold, "Light Sampling Threshold", 0.05f);
+  SOCKET_BOOLEAN(use_light_tree, "Use light tree to optimize many light sampling", false);
+  SOCKET_FLOAT(splitting_threshold,
+               "Amount of lights to sample at a time, from one light at 0.0, to adaptively more "
+               "lights as needed, to all lights at 1.0",
+               0.0f);
 
   static NodeEnum method_enum;
   method_enum.insert("path", PATH);

intern/cycles/render/integrator.h

@@ -74,6 +74,8 @@ class Integrator : public Node {
   bool sample_all_lights_direct;
   bool sample_all_lights_indirect;
   float light_sampling_threshold;
+  bool use_light_tree;
+  float splitting_threshold;
 
   int adaptive_min_samples;
   float adaptive_threshold;

intern/cycles/render/light.cpp

@@ -20,6 +20,7 @@
 #include "render/film.h"
 #include "render/graph.h"
 #include "render/integrator.h"
+#include "render/light_tree.h"
 #include "render/mesh.h"
 #include "render/nodes.h"
 #include "render/object.h"
@@ -274,7 +275,89 @@ bool LightManager::object_usable_as_light(Object *object)
   return false;
 }
 
-void LightManager::device_update_distribution(Device *,
+float LightManager::distant_lights_energy(const Scene *scene, const vector<Primitive> &prims)
+{
+  float luminance = 0.0f;
+  float3 emission;
+  foreach (Primitive prim, prims) {
+
+    if (prim.prim_id >= 0)
+      continue;
+
+    const Light *lamp = scene->lights[prim.lamp_id];
+    if (lamp->type != LIGHT_DISTANT)
+      continue;
+
+    /* get emission from shader */
+    bool is_constant_emission = lamp->shader->is_constant_emission(&emission);
+    if (!is_constant_emission)
+      continue;  // TODO: Properly handle this case
+
+    luminance += scene->shader_manager->linear_rgb_to_gray(emission);
+  }
+
+  /* TODO: could project each bbox onto a disk outside the scene and sum up
+   * all the projected areas instead if this results in too high sampling */
+
+  /* get radius of bounding sphere of scene */
+  BoundBox scene_bbox = BoundBox::empty;
+  foreach (Object *object, scene->objects) {
+    // TODO: What about transforms?
+    scene_bbox.grow(object->bounds);
+  }
+
+  float radius_squared = len_squared(scene_bbox.max - scene_bbox.center());
+
+  return M_PI_F * radius_squared * luminance;
+}
+
+float LightManager::background_light_energy(Device *device,
+                                            DeviceScene *dscene,
+                                            Scene *scene,
+                                            Progress &progress,
+                                            const vector<Primitive> &prims)
+{
+  /* compute energy for all background lights */
+  float average_luminance = 0.0f;
+  size_t num_pixels = 0;
+  /* find background lights */
+  foreach (Primitive prim, prims) {
+
+    if (prim.prim_id >= 0)
+      continue;
+
+    const Light *lamp = scene->lights[prim.lamp_id];
+    if (lamp->type != LIGHT_BACKGROUND)
+      continue;
+
+    vector<float3> pixels;
+    int2 res = get_background_map_resolution(lamp, scene);
+    shade_background_pixels(device, dscene, res.x, res.y, pixels, progress);
+    num_pixels += pixels.size();
+    for (int i = 0; i < pixels.size(); ++i) {
+      average_luminance += scene->shader_manager->linear_rgb_to_gray(pixels[i]);
+    }
+
+    break;
+  }
+
+  if (num_pixels == 0)
+    return 0.0f;
+
+  average_luminance /= (float)num_pixels;
+
+  /* get radius of bounding sphere of scene */
+  BoundBox scene_bbox = BoundBox::empty;
+  foreach (Object *object, scene->objects) {
+    // TODO: What about transforms?
+    scene_bbox.grow(object->bounds);
+  }
+  float radius_squared = len_squared(scene_bbox.max - scene_bbox.center());
+
+  return M_PI_F * radius_squared * average_luminance;
+}
+
+void LightManager::device_update_distribution(Device *device,
                                               DeviceScene *dscene,
                                               Scene *scene,
                                               Progress &progress)
@@ -285,28 +368,35 @@ void LightManager::device_update_distribution(Device *,
   size_t num_lights = 0;
   size_t num_portals = 0;
   size_t num_background_lights = 0;
+  size_t num_distant_lights = 0;
   size_t num_triangles = 0;
 
   bool background_mis = false;
 
-  foreach (Light *light, scene->lights) {
-    if (light->is_enabled) {
-      num_lights++;
-    }
-    if (light->is_portal) {
-      num_portals++;
-    }
-  }
+  /* The emissive_prims vector contains all emissive primitives in the scene,
+   * i.e., all mesh light triangles and all lamps. The order of the primitives
+   * in the vector is important since it has the same order as the
+   * light_distribution array.
+   *
+   * If using the light tree then the order is important since the light tree
+   * reordered the lights so lights in the same node are next to each other
+   * in memory.
+   *
+   * If NOT using the light tree then the order is important since during
+   * sampling we assume all triangles are first in the array. */
+  vector<Primitive> emissive_prims;
+  emissive_prims.reserve(scene->lights.size());
 
+  int object_id = 0;
   foreach (Object *object, scene->objects) {
     if (progress.get_cancel())
       return;
 
     if (!object_usable_as_light(object)) {
+      object_id++;
       continue;
     }
-
-    /* Count triangles. */
+    /* Count emissive triangles. */
     Mesh *mesh = static_cast<Mesh *>(object->geometry);
     size_t mesh_num_triangles = mesh->num_triangles();
     for (size_t i = 0; i < mesh_num_triangles; i++) {
@@ -316,35 +406,272 @@ void LightManager::device_update_distribution(Device *,
                            scene->default_surface;
 
       if (shader->use_mis && shader->has_surface_emission) {
+        emissive_prims.push_back(Primitive(i + mesh->prim_offset, object_id));
         num_triangles++;
       }
     }
+
+    object_id++;
+  }
+
+  /* light index is the index of this lamp in the device lights array*/
+  int light_index = 0;
+
+  /* light_id is the index of this lamp in the scene lights array */
+  int light_id = 0;
+
+  /* the light index of the background light */
+  int background_index = -1;
+  foreach (Light *light, scene->lights) {
+    if (light->is_enabled) {
+      emissive_prims.push_back(Primitive(~light_index, light_id));
+      num_lights++;
+      if (light->type == LIGHT_BACKGROUND) {
+        background_index = light_index;
+      }
+      light_index++;
+    }
+    if (light->is_portal) {
+      num_portals++;
+    }
+    light_id++;
   }
 
   size_t num_distribution = num_triangles + num_lights;
   VLOG(1) << "Total " << num_distribution << " of light distribution primitives.";
 
+  if (scene->integrator->use_light_tree) {
+
+    /* create light tree */
+    double time_start = time_dt();
+    LightTree light_tree(emissive_prims, scene, 64);
+    VLOG(1) << "Light tree build time: " << time_dt() - time_start;
+
+    /* the light tree reorders the primitives so update emissive_prims */
+    const vector<Primitive> &ordered_prims = light_tree.get_primitives();
+    emissive_prims = ordered_prims;
+
+    if (progress.get_cancel())
+      return;
+
+    /* create the nodes to be used on the device */
+    const vector<CompactNode> &nodes = light_tree.get_nodes();
+    float4 *dnodes = dscene->light_tree_nodes.alloc(nodes.size() * LIGHT_TREE_NODE_SIZE);
+
+    /* convert each compact node into 4xfloat4
+     * 4 for energy, right_child_offset, prim_id, num_emitters
+     * 4 for bbox.min + bbox.max[0]
+     * 4 for bbox.max[1-2], theta_o, theta_e
+     * 4 for axis + energy variance */
+    size_t offset = 0;
+    size_t num_leaf_lights = 0;
+    foreach (CompactNode node, nodes) {
+      dnodes[offset].x = node.energy;
+      dnodes[offset].y = __int_as_float(node.right_child_offset);
+      dnodes[offset].z = __int_as_float(node.first_prim_offset);
+      dnodes[offset].w = __int_as_float(node.num_lights);
+
+      dnodes[offset + 1].x = node.bounds_s.min[0];
+      dnodes[offset + 1].y = node.bounds_s.min[1];
+      dnodes[offset + 1].z = node.bounds_s.min[2];
+      dnodes[offset + 1].w = node.bounds_s.max[0];
+
+      dnodes[offset + 2].x = node.bounds_s.max[1];
+      dnodes[offset + 2].y = node.bounds_s.max[2];
+      dnodes[offset + 2].z = node.bounds_o.theta_o;
+      dnodes[offset + 2].w = node.bounds_o.theta_e;
+
+      dnodes[offset + 3].x = node.bounds_o.axis[0];
+      dnodes[offset + 3].y = node.bounds_o.axis[1];
+      dnodes[offset + 3].z = node.bounds_o.axis[2];
+      dnodes[offset + 3].w = node.energy_variance;
+
+      offset += 4;
+
+      if ((node.right_child_offset == -1) && (node.num_lights > 1)) {
+        num_leaf_lights += node.num_lights;
+      }
+    }
+
+    /* store information needed for importance computations for each emitter
+     * in leaf nodes containing several emitters.
+     *
+     * each leaf node with several emitters stores relevant information about
+     * its emitters in the light_tree_leaf_emitters array. each such node
+     * also stores an offset into the light_tree_leaf_emitters array to where
+     * its first light is. this offset is stored in leaf_to_first_emitter.
+     */
+    float4 *leaf_emitters = dscene->light_tree_leaf_emitters.alloc(num_leaf_lights * 3);
+    int *leaf_to_first_emitter = dscene->leaf_to_first_emitter.alloc(nodes.size());
+
+    offset = 0;
+    for (int i = 0; i < nodes.size(); ++i) {
+      const CompactNode &node = nodes[i];
+
+      /* only store this information for leaf nodes with several emitters */
+      if (!((node.right_child_offset == -1) && (node.num_lights > 1))) {
+        leaf_to_first_emitter[i] = -1;
+        continue;
+      }
+
+      leaf_to_first_emitter[i] = offset;
+
+      int start = node.first_prim_offset;  // distribution id
+      int end = start + node.num_lights;
+      for (int j = start; j < end; ++j) {
+
+        /* todo: is there a better way than recalcing this? */
+        /* have getters for the light tree that just accesses build_data? */
+        BoundBox bbox = light_tree.compute_bbox(emissive_prims[j]);
+        Orientation bcone = light_tree.compute_bcone(emissive_prims[j]);
+        float energy = light_tree.compute_energy(emissive_prims[j]);
+
+        leaf_emitters[offset].x = bbox.min[0];
+        leaf_emitters[offset].y = bbox.min[1];
+        leaf_emitters[offset].z = bbox.min[2];
+        leaf_emitters[offset].w = bbox.max[0];
+
+        leaf_emitters[offset + 1].x = bbox.max[1];
+        leaf_emitters[offset + 1].y = bbox.max[2];
+        leaf_emitters[offset + 1].z = bcone.theta_o;
+        leaf_emitters[offset + 1].w = bcone.theta_e;
+
+        leaf_emitters[offset + 2].x = bcone.axis[0];
+        leaf_emitters[offset + 2].y = bcone.axis[1];
+        leaf_emitters[offset + 2].z = bcone.axis[2];
+        leaf_emitters[offset + 2].w = energy;
+        offset += 3;
+      }
+    }
+
+    if (progress.get_cancel())
+      return;
+
+    /* create CDF for distant lights, background lights and light tree */
+    float tree_energy = (nodes.size() > 0) ? nodes[0].energy : 0.0f;
+    float distant_energy = distant_lights_energy(scene, emissive_prims);
+    float background_energy = background_light_energy(
+        device, dscene, scene, progress, emissive_prims);
+
+    /* stores the function that the CDF will be generated from */
+    float3 func = make_float3(tree_energy, distant_energy, background_energy);
+
+    /* probs stores the probability of sampling each of the light groups.
+     * probs[0] corresponds to the probability to sample the tree, etc. */
+    float3 probs;
+    float4 cdf = make_float4(0.0f);
+    const int num_func_values = LIGHTGROUP_NUM;
+    const int num_cdf_values = num_func_values + 1;
+    for (int i = 1; i < num_cdf_values; ++i) {
+      cdf[i] = cdf[i - 1] + func[i - 1];
+    }
+    float func_integral = cdf[num_func_values];
+    if (func_integral == 0.0f) {  // Sample uniformly if no energy
+      for (int i = 1; i < num_cdf_values; ++i) {
+        cdf[i] = (float)i / num_func_values;
+      }
+      probs = make_float3(1.0f / (float)num_func_values);
+    }
+    else {
+      for (int i = 0; i < num_cdf_values; ++i) {
+        cdf[i] /= func_integral;
+      }
+      probs = func / func_integral;
+    }
+
+    /* create and fill device arrays for the light group probabilities and CDF */
+    float *type_cdf = dscene->light_group_sample_cdf.alloc(num_cdf_values);
+    for (int i = 0; i < num_cdf_values; ++i) {
+      type_cdf[i] = cdf[i];
+    }
+
+    float *type_prob = dscene->light_group_sample_prob.alloc(num_func_values);
+    for (int i = 0; i < num_func_values; ++i) {
+      type_prob[i] = probs[i];
+    }
+
+    if (progress.get_cancel())
+      return;
+
+    /* find mapping between distribution_id to node_id, used for MIS */
+    uint *distribution_to_node = dscene->light_distribution_to_node.alloc(num_distribution);
+
+    for (int i = 0; i < nodes.size(); ++i) {
+      const CompactNode &node = nodes[i];
+      if (node.right_child_offset != -1)
+        continue;  // Skip interior nodes
+
+      int start = node.first_prim_offset;  // distribution id
+      int end = start + node.num_lights;
+      for (int j = start; j < end; ++j) {
+        distribution_to_node[j] = 4 * i;
+      }
+    }
+  }
+
+  if (progress.get_cancel())
+    return;
+
+  /* find mapping between lamp_id to distribution_id, used for MIS */
+  uint *lamp_to_distribution = dscene->lamp_to_distribution.alloc(num_lights);
+  for (int i = 0; i < emissive_prims.size(); ++i) {
+    const Primitive &prim = emissive_prims[i];
+    if (prim.prim_id >= 0)
+      continue;                       // Skip triangles
+    int lamp_id = -prim.prim_id - 1;  // This should not use prim.lamp_id
+    lamp_to_distribution[lamp_id] = i;
+  }
+
+  /* find mapping between [triangle_id, object_id] to distribution_id, used for MIS */
+  /* tri_to_distr has the following format:
+   * [triangle_id0, object_id0, distrib_id0, triangle_id1,..]
+   * where [triangle_idX, object_idX] is mapped to distrib_idX. */
+  vector<std::tuple<uint, uint, uint>> tri_to_distr;
+  tri_to_distr.reserve(num_triangles);
+  for (int i = 0; i < emissive_prims.size(); ++i) {
+    const Primitive &prim = emissive_prims[i];
+    if (prim.prim_id < 0)
+      continue;  // Skip lamps
+    tri_to_distr.push_back(std::make_tuple(prim.prim_id, prim.object_id, i));
+  }
+
+  std::sort(tri_to_distr.begin(), tri_to_distr.end());
+
+  assert(num_triangles == tri_to_distr.size());
+  uint *triangle_to_distribution = dscene->triangle_to_distribution.alloc(num_triangles * 3);
+  for (int i = 0; i < tri_to_distr.size(); ++i) {
+    triangle_to_distribution[3 * i] = std::get<0>(tri_to_distr[i]);
+    triangle_to_distribution[3 * i + 1] = std::get<1>(tri_to_distr[i]);
+    triangle_to_distribution[3 * i + 2] = std::get<2>(tri_to_distr[i]);
+  }
+
+  if (progress.get_cancel())
+    return;
+
+  /* create light distribution in same order as the emissive_prims */
+
   /* emission area */
   KernelLightDistribution *distribution = dscene->light_distribution.alloc(num_distribution + 1);
   float totarea = 0.0f;
 
   /* triangles */
   size_t offset = 0;
-  int j = 0;
 
-  foreach (Object *object, scene->objects) {
+  assert(emissive_prims.size() == num_distribution);
+
+  /* create distributions for mesh lights */
+  foreach (Primitive prim, emissive_prims) {
     if (progress.get_cancel())
       return;
 
-    if (!object_usable_as_light(object)) {
-      j++;
+    if (prim.prim_id < 0) {  // Early exit for lights
+      offset++;
       continue;
     }
+
+    const Object *object = scene->objects[prim.object_id];
     /* Sum area. */
     Mesh *mesh = static_cast<Mesh *>(object->geometry);
-    bool transform_applied = mesh->transform_applied;
-    Transform tfm = object->tfm;
-    int object_id = j;
     int shader_flag = 0;
 
     if (!(object->visibility & PATH_RAY_DIFFUSE)) {
@@ -364,39 +691,16 @@ void LightManager::device_update_distribution(Device *,
       use_light_visibility = true;
     }
 
-    size_t mesh_num_triangles = mesh->num_triangles();
-    for (size_t i = 0; i < mesh_num_triangles; i++) {
-      int shader_index = mesh->shader[i];
-      Shader *shader = (shader_index < mesh->used_shaders.size()) ?
-                           mesh->used_shaders[shader_index] :
-                           scene->default_surface;
+    int triangle_id = prim.prim_id - mesh->prim_offset;
+    const float area = mesh->compute_triangle_area(triangle_id, object->tfm);
+    distribution[offset].area = area;
+    distribution[offset].totarea = totarea;
+    distribution[offset].prim = prim.prim_id;
+    distribution[offset].mesh_light.shader_flag = shader_flag;
+    distribution[offset].mesh_light.object_id = prim.object_id;
+    offset++;
 
-      if (shader->use_mis && shader->has_surface_emission) {
-        distribution[offset].totarea = totarea;
-        distribution[offset].prim = i + mesh->prim_offset;
-        distribution[offset].mesh_light.shader_flag = shader_flag;
-        distribution[offset].mesh_light.object_id = object_id;
-        offset++;
-
-        Mesh::Triangle t = mesh->get_triangle(i);
-        if (!t.valid(&mesh->verts[0])) {
-          continue;
-        }
-        float3 p1 = mesh->verts[t.v[0]];
-        float3 p2 = mesh->verts[t.v[1]];
-        float3 p3 = mesh->verts[t.v[2]];
-
-        if (!transform_applied) {
-          p1 = transform_point(&tfm, p1);
-          p2 = transform_point(&tfm, p2);
-          p3 = transform_point(&tfm, p3);
-        }
-
-        totarea += triangle_area(p1, p2, p3);
-      }
-    }
-
-    j++;
+    totarea += area;
   }
 
   float trianglearea = totarea;
@@ -404,19 +708,29 @@ void LightManager::device_update_distribution(Device *,
   /* point lights */
   float lightarea = (totarea > 0.0f) ? totarea / num_lights : 1.0f;
   bool use_lamp_mis = false;
+  offset = 0;
 
-  int light_index = 0;
-  foreach (Light *light, scene->lights) {
-    if (!light->is_enabled)
+  /* create distributions for lights */
+  foreach (Primitive prim, emissive_prims) {
+    if (progress.get_cancel())
+      return;
+
+    if (prim.prim_id >= 0) {  // Early exit for mesh lights
+      offset++;
       continue;
+    }
 
+    const Light *light = scene->lights[prim.lamp_id];
+
+    distribution[offset].area = lightarea;
     distribution[offset].totarea = totarea;
-    distribution[offset].prim = ~light_index;
+    distribution[offset].prim = prim.prim_id;
     distribution[offset].lamp.pad = 1.0f;
     distribution[offset].lamp.size = light->size;
     totarea += lightarea;
 
     if (light->type == LIGHT_DISTANT) {
+      num_distant_lights++;
       use_lamp_mis |= (light->angle > 0.0f && light->use_mis);
     }
     else if (light->type == LIGHT_POINT || light->type == LIGHT_SPOT) {
@@ -430,11 +744,11 @@ void LightManager::device_update_distribution(Device *,
       background_mis |= light->use_mis;
     }
 
-    light_index++;
     offset++;
   }
 
   /* normalize cumulative distribution functions */
+  distribution[num_distribution].area = 0.0f;
   distribution[num_distribution].totarea = totarea;
   distribution[num_distribution].prim = 0.0f;
   distribution[num_distribution].lamp.pad = 0.0f;
@@ -456,12 +770,34 @@ void LightManager::device_update_distribution(Device *,
   kintegrator->use_direct_light = (totarea > 0.0f);
 
   if (kintegrator->use_direct_light) {
+
+    /* update light tree */
+    kintegrator->use_light_tree = scene->integrator->use_light_tree;
+    kintegrator->splitting_threshold = scene->integrator->splitting_threshold;
+
+    dscene->light_tree_nodes.copy_to_device();
+    dscene->light_distribution_to_node.copy_to_device();
+    dscene->lamp_to_distribution.copy_to_device();
+    dscene->triangle_to_distribution.copy_to_device();
+    dscene->light_group_sample_cdf.copy_to_device();
+    dscene->light_group_sample_prob.copy_to_device();
+    dscene->leaf_to_first_emitter.copy_to_device();
+    dscene->light_tree_leaf_emitters.copy_to_device();
+    kintegrator->num_light_nodes = dscene->light_tree_nodes.size() / LIGHT_TREE_NODE_SIZE;
+    // TODO: Currently this is only the correct offset when using light tree
+    kintegrator->distant_lights_offset = num_distribution - num_distant_lights;
+    kintegrator->background_light_index = background_index;
+
     /* number of emissives */
     kintegrator->num_distribution = num_distribution;
+    kintegrator->num_triangle_lights = num_triangles;
+    kintegrator->num_distant_lights = num_distant_lights;
+    kintegrator->inv_num_distant_lights = 1.0f / (float)num_distant_lights;
 
     /* precompute pdfs */
     kintegrator->pdf_triangles = 0.0f;
     kintegrator->pdf_lights = 0.0f;
+    kintegrator->pdf_inv_totarea = 1.0f / totarea;
 
     /* sample one, with 0.5 probability of light or triangle */
     kintegrator->num_all_lights = num_lights;
@@ -509,10 +845,24 @@ void LightManager::device_update_distribution(Device *,
     kbackground->map_weight = background_mis ? 1.0f : 0.0f;
   }
   else {
+    dscene->light_group_sample_cdf.free();
+    dscene->light_group_sample_prob.free();
+    dscene->leaf_to_first_emitter.free();
+    dscene->light_tree_leaf_emitters.free();
     dscene->light_distribution.free();
+    dscene->light_tree_nodes.free();
+    dscene->light_distribution_to_node.free();
+    dscene->lamp_to_distribution.free();
+    dscene->triangle_to_distribution.free();
 
+    kintegrator->pdf_inv_totarea = 0.0f;
+    kintegrator->num_light_nodes = 0;
+    kintegrator->num_triangle_lights = 0;
+    kintegrator->use_light_tree = false;
     kintegrator->num_distribution = 0;
     kintegrator->num_all_lights = 0;
+    kintegrator->num_distant_lights = 0;
+    kintegrator->inv_num_distant_lights = 0.0f;
     kintegrator->pdf_triangles = 0.0f;
     kintegrator->pdf_lights = 0.0f;
     kintegrator->use_lamp_mis = false;
@@ -521,8 +871,9 @@ void LightManager::device_update_distribution(Device *,
     kbackground->portal_offset = 0;
     kbackground->portal_weight = 0.0f;
     kbackground->sun_weight = 0.0f;
+    kintegrator->distant_lights_offset = 0;
+    kintegrator->background_light_index = 0;
     kbackground->map_weight = 0.0f;
-
     kfilm->pass_shadow_scale = 1.0f;
   }
 }
@@ -662,14 +1013,6 @@ void LightManager::device_update_background(Device *device,
   /* If the resolution isn't set manually, try to find an environment texture. */
   if (res.x == 0) {
     res = environment_res;
-    if (res.x > 0 && res.y > 0) {
-      VLOG(2) << "Automatically set World MIS resolution to " << res.x << " by " << res.y << "\n";
-    }
-  }
-  /* If it's still unknown, just use the default. */
-  if (res.x == 0 || res.y == 0) {
-    res = make_int2(1024, 512);
-    VLOG(2) << "Setting World MIS resolution to default\n";
   }
   kbackground->map_res_x = res.x;
   kbackground->map_res_y = res.y;
@@ -1000,11 +1343,19 @@ void LightManager::device_update(Device *device,
 void LightManager::device_free(Device *, DeviceScene *dscene, const bool free_background)
 {
   dscene->light_distribution.free();
+  dscene->light_tree_nodes.free();
+  dscene->light_distribution_to_node.free();
+  dscene->lamp_to_distribution.free();
+  dscene->triangle_to_distribution.free();
   dscene->lights.free();
   if (free_background) {
     dscene->light_background_marginal_cdf.free();
     dscene->light_background_conditional_cdf.free();
   }
+  dscene->light_group_sample_prob.free();
+  dscene->light_group_sample_cdf.free();
+  dscene->leaf_to_first_emitter.free();
+  dscene->light_tree_leaf_emitters.free();
   dscene->ies_lights.free();
 }
 
@@ -1131,4 +1482,32 @@ void LightManager::device_update_ies(DeviceScene *dscene)
   }
 }
 
+int2 LightManager::get_background_map_resolution(const Light *background_light, const Scene *scene)
+{
+  int2 res = make_int2(background_light->map_resolution, background_light->map_resolution / 2);
+  /* If the resolution isn't set manually, try to find an environment texture. */
+  if (res.x == 0) {
+    Shader *shader = (scene->background->shader) ? scene->background->shader :
+                                                   scene->default_background;
+    foreach (ShaderNode *node, shader->graph->nodes) {
+      if (node->type == EnvironmentTextureNode::node_type) {
+        EnvironmentTextureNode *env = (EnvironmentTextureNode *)node;
+        ImageMetaData metadata = env->handle.metadata();
+        res.x = max(res.x, metadata.width);
+        res.y = max(res.y, metadata.height);
+      }
+    }
+    if (res.x > 0 && res.y > 0) {
+      VLOG(2) << "Automatically set World MIS resolution to " << res.x << " by " << res.y << "\n";
+    }
+  }
+  /* If it's still unknown, just use the default. */
+  if (res.x == 0 || res.y == 0) {
+    res = make_int2(1024, 512);
+    VLOG(2) << "Setting World MIS resolution to default\n";
+  }
+
+  return res;
+}
+
 CCL_NAMESPACE_END

intern/cycles/render/light.h

@@ -31,6 +31,7 @@ CCL_NAMESPACE_BEGIN
 class Device;
 class DeviceScene;
 class Object;
+class Primitive;
 class Progress;
 class Scene;
 class Shader;
@@ -128,6 +129,18 @@ class LightManager {
   /* Check whether light manager can use the object as a light-emissive. */
   bool object_usable_as_light(Object *object);
 
+  /* compute energy of background lights */
+  float background_light_energy(Device *device,
+                                DeviceScene *dscene,
+                                Scene *scene,
+                                Progress &progress,
+                                const vector<Primitive> &prims);
+
+  /* compute energy of distant lights */
+  float distant_lights_energy(const Scene *scene, const vector<Primitive> &prims);
+
+  int2 get_background_map_resolution(const Light *background_light, const Scene *scene);
+
   struct IESSlot {
     IESFile ies;
     uint hash;

intern/cycles/render/light_tree.cpp

@@ -0,0 +1,601 @@
+/*
+ * Copyright 2011-2018 Blender Foundation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "render/light_tree.h"
+#include "render/light.h"
+#include "render/mesh.h"
+#include "render/object.h"
+
+#include "util/util_foreach.h"
+#include "util/util_logging.h"
+
+CCL_NAMESPACE_BEGIN
+
+LightTree::LightTree(const vector<Primitive> &prims_,
+                     Scene *scene_,
+                     const uint max_lights_in_node_)
+    : max_lights_in_node(max_lights_in_node_), scene(scene_)
+{
+
+  if (prims_.empty())
+    return;
+
+  /* background and distant lights are not added to the light tree and are
+   * considered seperately. so here all primitives except background and
+   * distant lights are moved into a local primitives array */
+  primitives.reserve(prims_.size());
+  vector<Primitive> distant_lights;
+  vector<Primitive> background_lights;
+  foreach (Primitive prim, prims_) {
+
+    /* put background and distant lights into their own arrays */
+    if (prim.prim_id < 0) {
+      const Light *lamp = scene->lights[prim.lamp_id];
+      if (lamp->type == LIGHT_DISTANT) {
+        distant_lights.push_back(prim);
+        continue;
+      }
+      else if (lamp->type == LIGHT_BACKGROUND) {
+        background_lights.push_back(prim);
+        continue;
+      }
+    }
+
+    primitives.push_back(prim);
+  }
+
+  /* initialize build_data array that stores the energy and spatial and
+   * orientation bounds for each light. */
+  vector<BVHPrimitiveInfo> build_data;
+  build_data.reserve(primitives.size());
+  for (int i = 0; i < primitives.size(); ++i) {
+    BoundBox bbox = compute_bbox(primitives[i]);
+    Orientation bcone = compute_bcone(primitives[i]);
+    float energy = compute_energy(primitives[i]);
+
+    build_data.push_back(BVHPrimitiveInfo(i, bbox, bcone, energy));
+  }
+
+  /* recursively build BVH tree */
+  uint total_nodes = 0;
+  vector<Primitive> ordered_prims;
+  ordered_prims.reserve(primitives.size());
+  BVHBuildNode *root = recursive_build(
+      0, primitives.size(), build_data, total_nodes, ordered_prims);
+
+  /* order the primitives array so lights belonging to the same node are
+   * next to each other */
+  primitives.swap(ordered_prims);
+  ordered_prims.clear();
+  build_data.clear();
+
+  /* add background lights to the primitives array */
+  for (int i = 0; i < background_lights.size(); ++i) {
+    primitives.push_back(background_lights[i]);
+  }
+
+  /* add distant lights to the end of primitives array */
+  for (int i = 0; i < distant_lights.size(); ++i) {
+    primitives.push_back(distant_lights[i]);
+  }
+
+  VLOG(1) << "Total BVH nodes: " << total_nodes;
+
+  if (!root)
+    return;
+
+  /* convert to linear representation of the tree */
+  nodes.resize(total_nodes);
+  int offset = 0;
+  flattenBVHTree(*root, offset);
+
+  assert(offset == total_nodes);
+}
+
+int LightTree::flattenBVHTree(const BVHBuildNode &node, int &offset)
+{
+
+  CompactNode &compact_node = nodes[offset];
+  compact_node.bounds_s = node.bbox;
+  compact_node.bounds_o = node.bcone;
+
+  int my_offset = offset++;
+  if (node.is_leaf) {
+    /* create leaf node */
+    assert(!node.children[0] && !node.children[1]);
+    compact_node.energy = node.energy;
+    compact_node.energy_variance = node.energy_variance;
+    compact_node.first_prim_offset = node.first_prim_offset;
+    compact_node.num_lights = node.num_lights;
+  }
+  else {
+    /* create interior compact node */
+    compact_node.num_lights = node.num_lights;
+    compact_node.energy = node.energy;
+    compact_node.energy_variance = node.energy_variance;
+    assert(node.children[0] && node.children[1]);
+    flattenBVHTree(*node.children[0], offset);
+    compact_node.right_child_offset = flattenBVHTree(*node.children[1], offset);
+    compact_node.energy = node.energy;
+  }
+
+  return my_offset;
+}
+
+BoundBox LightTree::compute_bbox(const Primitive &prim)
+{
+  BoundBox bbox = BoundBox::empty;
+  if (prim.prim_id >= 0) {
+    /* extract bounding box from emissive triangle */
+    const Object *object = scene->objects[prim.object_id];
+    const Mesh *mesh = (Mesh *)object->geometry;
+    const int triangle_id = prim.prim_id - mesh->prim_offset;
+    const Mesh::Triangle triangle = mesh->get_triangle(triangle_id);
+
+    float3 p0 = mesh->verts[triangle.v[0]];
+    float3 p1 = mesh->verts[triangle.v[1]];
+    float3 p2 = mesh->verts[triangle.v[2]];
+
+    /* instanced mesh lights have not applied their transform at this point.
+     * in this case, these points have to be transformed to get the proper
+     * spatial bound. */
+    if (!mesh->transform_applied) {
+      const Transform &tfm = object->tfm;
+      p0 = transform_point(&tfm, p0);
+      p1 = transform_point(&tfm, p1);
+      p2 = transform_point(&tfm, p2);
+    }
+
+    bbox.grow(p0);
+    bbox.grow(p1);
+    bbox.grow(p2);
+  }
+  else {
+    /* extract bounding box from lamp based on light type */
+    Light *lamp = scene->lights[prim.lamp_id];
+    if (lamp->type == LIGHT_POINT || lamp->type == LIGHT_SPOT) {
+      float radius = lamp->size;
+      bbox.grow(lamp->co + make_float3(radius));
+      bbox.grow(lamp->co - make_float3(radius));
+    }
+    else if (lamp->type == LIGHT_AREA) {
+      const float3 center = lamp->co;
+      const float3 half_axisu = 0.5f * lamp->axisu * (lamp->sizeu * lamp->size);
+      const float3 half_axisv = 0.5f * lamp->axisv * (lamp->sizev * lamp->size);
+      const float3 p0 = center - half_axisu - half_axisv;
+      const float3 p1 = center - half_axisu + half_axisv;
+      const float3 p2 = center + half_axisu - half_axisv;
+      const float3 p3 = center + half_axisu + half_axisv;
+
+      bbox.grow(p0);
+      bbox.grow(p1);
+      bbox.grow(p2);
+      bbox.grow(p3);
+    }
+    else {
+      /* LIGHT_DISTANT and LIGHT_BACKGROUND are handled separately */
+      assert(false);
+    }
+  }
+
+  return bbox;
+}
+
+Orientation LightTree::compute_bcone(const Primitive &prim)
+{
+  Orientation bcone;
+  if (prim.prim_id >= 0) {
+    /* extract bounding cone from emissive triangle */
+    const Object *object = scene->objects[prim.object_id];
+    const Mesh *mesh = (Mesh *)object->geometry;
+    const int triangle_id = prim.prim_id - mesh->prim_offset;
+    const Mesh::Triangle triangle = mesh->get_triangle(triangle_id);
+
+    float3 p0 = mesh->verts[triangle.v[0]];
+    float3 p1 = mesh->verts[triangle.v[1]];
+    float3 p2 = mesh->verts[triangle.v[2]];
+
+    if (!mesh->transform_applied) {
+      const Transform &tfm = object->tfm;
+      p0 = transform_point(&tfm, p0);
+      p1 = transform_point(&tfm, p1);
+      p2 = transform_point(&tfm, p2);
+    }
+
+    float3 normal = make_float3(1.0f, 0.0f, 0.0f);
+    const float3 norm = cross(p1 - p0, p2 - p0);
+    const float normlen = len(norm);
+    if (normlen != 0.0f) {
+      normal = norm / normlen;
+    }
+
+    bcone.axis = normal;
+    bcone.theta_o = 0.0f;
+    bcone.theta_e = M_PI_2_F;
+  }
+  else {
+    Light *lamp = scene->lights[prim.lamp_id];
+    bcone.axis = lamp->dir / len(lamp->dir);
+    if (lamp->type == LIGHT_POINT) {
+      bcone.theta_o = M_PI_F;
+      bcone.theta_e = M_PI_2_F;
+    }
+    else if (lamp->type == LIGHT_SPOT) {
+      bcone.theta_o = 0;
+      bcone.theta_e = lamp->spot_angle * 0.5f;
+    }
+    else if (lamp->type == LIGHT_AREA) {
+      bcone.theta_o = 0;
+      bcone.theta_e = M_PI_2_F;
+    }
+  }
+
+  return bcone;
+}
+
+float LightTree::compute_energy(const Primitive &prim)
+{
+  float3 emission = make_float3(0.0f);
+  Shader *shader = NULL;
+
+  if (prim.prim_id >= 0) {
+    /* extract shader from emissive triangle */
+    const Object *object = scene->objects[prim.object_id];
+    const Mesh *mesh = (Mesh *)object->geometry;
+    const int triangle_id = prim.prim_id - mesh->prim_offset;
+
+    int shader_index = mesh->shader[triangle_id];
+    shader = mesh->used_shaders.at(shader_index);
+
+    /* get emission from shader */
+    bool is_constant_emission = shader->is_constant_emission(&emission);
+    if (!is_constant_emission) {
+      emission = make_float3(1.0f);
+    }
+
+    const Transform &tfm = scene->objects[prim.object_id]->tfm;
+    float area = mesh->compute_triangle_area(triangle_id, tfm);
+
+    emission *= area * 4;
+  }
+  else {
+    const Light *light = scene->lights[prim.lamp_id];
+
+    emission = light->strength;
+
+    /* calculate the max emission in a single direction. */
+    if (light->type == LIGHT_POINT) {
+      emission /= M_PI_F;
+    }
+    else if (light->type == LIGHT_SPOT) {
+      emission /= M_PI_F;
+    }
+    else if (light->type == LIGHT_AREA) {
+    }
+    else {
+      /* LIGHT_DISTANT and LIGHT_BACKGROUND are handled separately */
+      assert(false);
+    }
+  }
+
+  return scene->shader_manager->linear_rgb_to_gray(emission);
+}
+
+Orientation LightTree::combine_bounding_cones(const vector<Orientation> &bcones)
+{
+
+  if (bcones.size() == 0) {
+    return Orientation();
+  }
+  else if (bcones.size() == 1) {
+    return bcones[0];
+  }
+
+  Orientation cone = bcones[0];
+  for (int i = 1; i < bcones.size(); ++i) {
+    cone = cone_union(cone, bcones[i]);
+  }
+
+  return cone;
+}
+
+/* Algorithm 1 */
+Orientation LightTree::cone_union(const Orientation &cone1, const Orientation &cone2)
+{
+  const Orientation *a = &cone1;
+  const Orientation *b = &cone2;
+  if (b->theta_o > a->theta_o) {
+    a = &cone2;
+    b = &cone1;
+  }
+
+  float theta_d = safe_acosf(dot(a->axis, b->axis));
+
+  float theta_e = fmaxf(a->theta_e, b->theta_e);
+  if (fminf(theta_d + b->theta_o, M_PI_F) <= a->theta_o) {
+    return Orientation(a->axis, a->theta_o, theta_e);
+  }
+
+  float theta_o = (a->theta_o + theta_d + b->theta_o) * 0.5f;
+  if (M_PI_F <= theta_o) {
+    return Orientation(a->axis, M_PI_F, theta_e);
+  }
+
+  float theta_r = theta_o - a->theta_o;
+  float3 axis = rotate_around_axis(a->axis, cross(a->axis, b->axis), theta_r);
+  axis = normalize(axis);
+  return Orientation(axis, theta_o, theta_e);
+}
+
+float LightTree::calculate_cone_measure(const Orientation &bcone)
+{
+  /* eq. 1 */
+  float theta_w = fminf(bcone.theta_o + bcone.theta_e, M_PI_F);
+  return M_2PI_F *
+         (1.0f - cosf(bcone.theta_o) + 0.5f * (theta_w - bcone.theta_o) * sinf(bcone.theta_o) +
+          0.25f * cosf(bcone.theta_o) - 0.25f * cosf(bcone.theta_o - 2.0f * theta_w));
+}
+
+void LightTree::split_saoh(const BoundBox &centroid_bbox,
+                           const vector<BVHPrimitiveInfo> &build_data,
+                           const int start,
+                           const int end,
+                           const int num_buckets,
+                           const float node_M_Omega,
+                           const BoundBox &node_bbox,
+                           float &min_cost,
+                           int &min_dim,
+                           int &min_bucket)
+{
+
+  struct BucketInfo {
+    BucketInfo() : count(0), energy(0.0f)
+    {
+      bounds = BoundBox::empty;
+    }
+
+    int count;
+    float energy;     // total energy
+    BoundBox bounds;  // bounds of all primitives
+    Orientation bcone;
+  };
+
+  min_cost = std::numeric_limits<float>::max();
+  min_bucket = -1;
+
+  const float extent_max = max3(centroid_bbox.size());
+  for (int dim = 0; dim < 3; ++dim) {
+
+    vector<BucketInfo> buckets(num_buckets);
+    vector<vector<Orientation>> bucket_bcones(num_buckets);
+
+    /* calculate total energy in each bucket and a bbox of it */
+    const float extent = centroid_bbox.max[dim] - centroid_bbox.min[dim];
+    if (extent == 0.0f) {  // All dims cannot be zero
+      continue;
+    }
+
+    const float extent_inv = 1.0f / extent;
+    for (unsigned int i = start; i < end; ++i) {
+      int bucket_id = (int)((float)num_buckets *
+                            (build_data[i].centroid[dim] - centroid_bbox.min[dim]) * extent_inv);
+      if (bucket_id == num_buckets)
+        bucket_id = num_buckets - 1;
+      buckets[bucket_id].count++;
+      buckets[bucket_id].energy += build_data[i].energy;
+      buckets[bucket_id].bounds.grow(build_data[i].bbox);
+      bucket_bcones[bucket_id].push_back(build_data[i].bcone);
+    }
+
+    for (unsigned int i = 0; i < num_buckets; ++i) {
+      if (buckets[i].count != 0) {
+        buckets[i].bcone = combine_bounding_cones(bucket_bcones[i]);
+      }
+    }
+
+    /* compute costs for splitting at bucket boundaries */
+    vector<float> cost(num_buckets - 1);
+    BoundBox bbox_L, bbox_R;
+    float energy_L, energy_R;
+    vector<Orientation> bcones_L, bcones_R;
+
+    for (int i = 0; i < num_buckets - 1; ++i) {
+      bbox_L = BoundBox::empty;
+      bbox_R = BoundBox::empty;
+      energy_L = 0;
+      energy_R = 0;
+      bcones_L.clear();
+      bcones_R.clear();
+
+      /* L corresponds to all buckets up to and including i */
+      for (int j = 0; j <= i; ++j) {
+        if (buckets[j].count != 0) {
+          energy_L += buckets[j].energy;
+          bbox_L.grow(buckets[j].bounds);
+          bcones_L.push_back(buckets[j].bcone);
+        }
+      }
+
+      /* R corresponds to bucket i+1 and all after */
+      for (int j = i + 1; j < num_buckets; ++j) {
+        if (buckets[j].count != 0) {
+          energy_R += buckets[j].energy;
+          bbox_R.grow(buckets[j].bounds);
+          bcones_R.push_back(buckets[j].bcone);
+        }
+      }
+
+      /* eq. 2 */
+      const Orientation bcone_L = combine_bounding_cones(bcones_L);
+      const Orientation bcone_R = combine_bounding_cones(bcones_R);
+      const float M_Omega_L = calculate_cone_measure(bcone_L);
+      const float M_Omega_R = calculate_cone_measure(bcone_R);
+      const float K = extent_max * extent_inv;
+
+      cost[i] = K * (energy_L * M_Omega_L * bbox_L.area() + energy_R * M_Omega_R * bbox_R.area()) /
+                (node_M_Omega * node_bbox.area());
+    }
+
+    /* update minimum cost, dim and bucket */
+    for (int i = 0; i < num_buckets - 1; ++i) {
+      if (cost[i] < min_cost) {
+        min_cost = cost[i];
+        min_dim = dim;
+        min_bucket = i;
+      }
+    }
+  }
+}
+
+BVHBuildNode *LightTree::recursive_build(const uint start,
+                                         const uint end,
+                                         vector<BVHPrimitiveInfo> &build_data,
+                                         uint &total_nodes,
+                                         vector<Primitive> &ordered_prims)
+{
+  if (build_data.size() == 0)
+    return NULL;
+
+  total_nodes++;
+  BVHBuildNode *node = new BVHBuildNode();
+
+  /* compute bounds and energy for all emissive primitives in node */
+  BoundBox node_bbox = BoundBox::empty;
+  vector<Orientation> bcones;
+  bcones.reserve(end - start);
+  double node_energy = 0.0;
+  double node_energy_sum_squared = 0.0;
+  uint num_lights = end - start;
+
+  for (unsigned int i = start; i < end; ++i) {
+    const BVHPrimitiveInfo &light = build_data.at(i);
+    node_bbox.grow(light.bbox);
+    bcones.push_back(light.bcone);
+
+    double energy = (double)light.energy;
+    node_energy += energy;
+    node_energy_sum_squared += energy * energy;
+  }
+
+  /* pre-calculate energy variance for the splitting heuristic */
+  double node_energy_mean = node_energy / (double)num_lights;
+  double node_energy_variance = node_energy_sum_squared / (double)num_lights -  // E[e^2]
+                                node_energy_mean * node_energy_mean;            // E[e]^2
+  node_energy_variance = max(node_energy_variance, 0.0);
+
+  Orientation node_bcone = combine_bounding_cones(bcones);
+  bcones.clear();
+  const float node_M_Omega = calculate_cone_measure(node_bcone);
+
+  if (num_lights == 1) {
+    /* create leaf */
+    int first_prim_offset = ordered_prims.size();
+    int prim = build_data.at(start).primitive_offset;
+    ordered_prims.push_back(primitives.at(prim));
+
+    node->init_leaf(
+        first_prim_offset, num_lights, node_bbox, node_bcone, node_energy, node_energy_variance);
+    return node;
+  }
+  else {
+    /* compute spatial bound for primitive centroids */
+    BoundBox centroid_bbox = BoundBox::empty;
+    for (unsigned int i = start; i < end; ++i) {
+      centroid_bbox.grow(build_data.at(i).centroid);
+    }
+
+    /* find dimension of bounding box with maximum extent */
+    float3 diag = centroid_bbox.size();
+    int max_dim;
+    if (diag[0] > diag[1] && diag[0] > diag[2]) {
+      max_dim = 0;
+    }
+    else if (diag[1] > diag[2]) {
+      max_dim = 1;
+    }
+    else {
+      max_dim = 2;
+    }
+
+    /* checks special case if all lights are in the same place */
+    if (centroid_bbox.max[max_dim] == centroid_bbox.min[max_dim]) {
+      /* create leaf */
+      int first_prim_offset = ordered_prims.size();
+      for (int i = start; i < end; ++i) {
+        int prim = build_data.at(i).primitive_offset;
+        ordered_prims.push_back(primitives.at(prim));
+      }
+
+      node->init_leaf(
+          first_prim_offset, num_lights, node_bbox, node_bcone, node_energy, node_energy_variance);
+
+      return node;
+    }
+    else {
+
+      /* find dimension and bucket with smallest SAOH cost */
+      const int num_buckets = 12;
+      float min_cost;
+      int min_dim, min_bucket;
+      split_saoh(centroid_bbox,
+                 build_data,
+                 start,
+                 end,
+                 num_buckets,
+                 node_M_Omega,
+                 node_bbox,
+                 min_cost,
+                 min_dim,
+                 min_bucket);
+      assert(min_dim != -1);
+
+      int mid = 0;
+      if (num_lights > max_lights_in_node || min_cost < (float)node_energy) {
+        /* partition primitives */
+        BVHPrimitiveInfo *mid_ptr = std::partition(
+            &build_data[start],
+            &build_data[end - 1] + 1,
+            CompareToBucket(min_bucket, num_buckets, min_dim, centroid_bbox));
+        mid = mid_ptr - &build_data[0];
+      }
+      else {
+        /* create leaf */
+        int first_prim_offset = ordered_prims.size();
+        for (int i = start; i < end; ++i) {
+          int prim = build_data.at(i).primitive_offset;
+          ordered_prims.push_back(primitives.at(prim));
+        }
+
+        node->init_leaf(first_prim_offset,
+                        num_lights,
+                        node_bbox,
+                        node_bcone,
+                        node_energy,
+                        node_energy_variance);
+        return node;
+      }
+
+      /* build children */
+      BVHBuildNode *left = recursive_build(start, mid, build_data, total_nodes, ordered_prims);
+      BVHBuildNode *right = recursive_build(mid, end, build_data, total_nodes, ordered_prims);
+      node->init_interior(left, right, node_bcone, num_lights, node_energy, node_energy_variance);
+    }
+  }
+
+  return node;
+}
+
+CCL_NAMESPACE_END

intern/cycles/render/light_tree.h

@@ -0,0 +1,356 @@
+/*
+ * Copyright 2011-2018 Blender Foundation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef __LIGHT_TREE_H__
+#define __LIGHT_TREE_H__
+
+#include "util/util_boundbox.h"
+
+CCL_NAMESPACE_BEGIN
+
+class Light;
+class Object;
+class Scene;
+
+#define LIGHT_TREE_NODE_SIZE 4
+
+/* Data structure to represent orientation bounds. It consists of two bounding
+ * cones represented by a direction(axis) and two angles out from this axis.
+ * This can be thought of as two cones.
+ */
+struct Orientation {
+  Orientation()
+  {
+    axis = make_float3(0.0f, 0.0f, 0.0f);
+    theta_o = 0;
+    theta_e = 0;
+  }
+
+  Orientation(const float3 &a, float o, float e) : axis(a), theta_o(o), theta_e(e)
+  {
+  }
+
+  /* orientation/direction of the cones */
+  float3 axis;
+
+  /* angle bounding light orientations */
+  float theta_o;
+
+  /* angle bounding the directions light can be emitted in */
+  float theta_e;
+};
+
+/* Temporary data structure for nodes during construction.
+ * After the construction is complete a final step converts the tree consisting
+ * of these nodes into a tree consisting of CompactNode:s. */
+struct BVHBuildNode {
+
+  BVHBuildNode()
+  {
+    children[0] = children[1] = NULL;
+    bbox = BoundBox::empty;
+  }
+
+  /* initializes this node as a leaf node */
+  void init_leaf(
+      uint first, uint n, const BoundBox &b, const Orientation &c, double e, double e_var)
+  {
+    first_prim_offset = first;
+    num_lights = n;
+    bbox = b;
+    bcone = c;
+    energy = (float)e;
+    energy_variance = (float)e_var;
+    is_leaf = true;
+  }
+
+  /* initializes this node as an interior node */
+  void init_interior(
+      BVHBuildNode *c0, BVHBuildNode *c1, const Orientation &c, uint n, double e, double e_var)
+  {
+    bbox = merge(c0->bbox, c1->bbox);
+    bcone = c;
+
+    children[0] = c0;
+    children[1] = c1;
+
+    num_lights = n;
+    energy = (float)e;
+    energy_variance = (float)e_var;
+    is_leaf = false;
+  }
+
+  /* spatial and orientation bounds */
+  BoundBox bbox;
+  Orientation bcone;
+
+  /* total energy and energy variance for the lights in the node */
+  float energy, energy_variance;
+
+  /* pointers to the two children */
+  BVHBuildNode *children[2];
+
+  /* each leaf node contains one or more lights. lights that are contained in
+   * the same node are stored next to each other in the ordered primitives
+   * array. this offset points to the first of these lights. num_lights below
+   * can be used to find the last light for this node */
+  uint first_prim_offset;
+
+  /* total number of lights contained in this node */
+  uint num_lights;
+
+  /* if this node is a leaf or not */
+  bool is_leaf;
+};
+
+// TODO: Have this struct in kernel_types.h instead?
+/* A more memory efficient representation of BVHBuildNode above. This is the
+ * structure of the nodes on the device. */
+struct CompactNode {
+
+  CompactNode()
+      : right_child_offset(-1),
+        first_prim_offset(-1),
+        num_lights(-1),
+        bounds_s(BoundBox::empty),
+        energy(0.0f),
+        energy_variance(0.0f)
+  {
+    bounds_o.axis = make_float3(0.0f);
+    bounds_o.theta_o = 0.0f;
+    bounds_o.theta_e = 0.0f;
+  }
+
+  /* All compact nodes are stored in a single array. interior nodes can find
+   * their two child nodes as follows:
+   * - the left child node is directly after its parent in the nodes array
+   * - the right child node is at the offset below in the nodes array.
+   *
+   * This offset is default initialized to -1 and will only change if this is
+   * and interior node. this therefore used to see if a node is a leaf/interior
+   * node as well. */
+  int right_child_offset;
+
+  /* see comment in BVHBuildNode for the variable with the same name */
+  int first_prim_offset;
+
+  /* total number of lights contained in this node */
+  int num_lights;
+
+  /* spatial and orientation bounds */
+  BoundBox bounds_s;
+  Orientation bounds_o;
+
+  /* total energy and energy variance for the lights in the node */
+  float energy, energy_variance;
+};
+
+/* Helper struct that is only used during the construction of the tree */
+struct BVHPrimitiveInfo {
+
+  BVHPrimitiveInfo()
+  {
+    bbox = BoundBox::empty;
+  }
+
+  BVHPrimitiveInfo(uint offset, const BoundBox &bounds, const Orientation &oBounds, float e)
+      : primitive_offset(offset), bbox(bounds), centroid(bounds.center()), energy(e)
+  {
+    bcone.axis = oBounds.axis;
+    bcone.theta_o = oBounds.theta_o;
+    bcone.theta_e = oBounds.theta_e;
+  }
+
+  /* this primitives offset into the unordered primtives array */
+  uint primitive_offset;
+
+  /* spatial and orientation bounds */
+  BoundBox bbox;
+  float3 centroid;
+  Orientation bcone;
+
+  /* total energy of this emissive primitive */
+  float energy;
+};
+
+/* A custom pointer struct that points to an emissive triangle or a lamp. */
+class Primitive {
+ public:
+  /* If prim_id >= 0 then the primitive is a triangle and prim_id is a global
+   * triangle index.
+   * If prim_id < 0 then the primitive is a lamp and -prim_id-1 is an index
+   * into the lights array on the device. */
+  int prim_id;
+  union {
+    /* which object the triangle belongs to */
+    int object_id;
+    /* index for this lamp in the scene->lights array */
+    int lamp_id;
+  };
+  Primitive(int prim, int object) : prim_id(prim), object_id(object)
+  {
+  }
+};
+
+/* Compare operator that returns true if the given light is in a lower
+ * bucket than a given split_bucket. This is used to partition lights into lights
+ * for the left and right child during tree construction. */
+struct CompareToBucket {
+  CompareToBucket(int split, int num, int d, const BoundBox &b) : centroid_bbox(b)
+  {
+    split_bucket = split;
+    num_buckets = num;
+    dim = d;
+    inv_extent = 1.0f / (centroid_bbox.max[dim] - centroid_bbox.min[dim]);
+  }
+
+  bool operator()(const BVHPrimitiveInfo &p) const
+  {
+    int bucket_id = (int)((float)num_buckets * (p.centroid[dim] - centroid_bbox.min[dim]) *
+                          inv_extent);
+    if (bucket_id == num_buckets) {
+      bucket_id = num_buckets - 1;
+    }
+
+    return bucket_id <= split_bucket;
+  }
+
+  /* everything lower or equal to the split_bucket is considered to be in one
+   * child and everything above will be considered to belong to the other. */
+  int split_bucket;
+
+  /* the total number of buckets that are considered for this dimension(dim) */
+  int num_buckets;
+
+  /* the construction creates candidate splits along the three dimensions.
+   * this variable stores which dimension is currently being split along.*/
+  int dim;
+
+  /* storing the inverse extend of the bounding box along the current
+   * dimension to only have to do the division once instead of everytime the
+   * operator() is called. */
+  float inv_extent;
+
+  /* bound for the centroids of all lights of the current node being split */
+  const BoundBox &centroid_bbox;
+};
+
+/* This class takes a set of lights as input and organizes them into a light
+ * hierarchy. This hierarchy is represented as a Bounding Volume Hierarchy(BVH).
+ * This is the process to acheive this:
+ * 1. For each given light, important information is gathered
+ *    - Bounding box of the light
+ *    - Bounding cones of the light
+ *    - The energy of the light
+ *    This first calculated and then stored as BVHPrimitiveInfo for each light.
+ * 2. A top-down recursive build algorithm creates a BVH consisting of
+ *    BVHBuildNode:s which each are allocated randomly on the heap with new.
+ *    This step also reorders the given array of lights such that lights
+ *    belonging to the same node are next to each other in the primitives array.
+ * 3. A final step converts this BVH into a more memory efficient layout where
+ *    each BVHBuildNode is converted to a CompactNode and all of these nodes
+ *    are placed next to each other in memory in a single nodes array.
+ *
+ *  This structure is based on PBRTs geometry BVH implementation.
+ **/
+class LightTree {
+ public:
+  LightTree(const vector<Primitive> &prims_, Scene *scene_, const uint max_lights_in_node_);
+
+  /* returns the ordered emissive primitives */
+  const vector<Primitive> &get_primitives() const
+  {
+    return primitives;
+  }
+
+  /* returns the array of nodes */
+  const vector<CompactNode> &get_nodes() const
+  {
+    return nodes;
+  }
+
+  /* computes the bounding box for the given light */
+  BoundBox compute_bbox(const Primitive &prim);
+
+  /* computes the orientation bounds for the given light. */
+  Orientation compute_bcone(const Primitive &prim);
+
+  /* computes the emitted energy for the given light. this is done by
+   * integrating the constant emission over the angles of the sphere it emits
+   * light in. */
+  float compute_energy(const Primitive &prim);
+
+ private:
+  /* the top-down recursive build algorithm mentioned in step 2 above */
+  BVHBuildNode *recursive_build(const uint start,
+                                const uint end,
+                                vector<BVHPrimitiveInfo> &build_data,
+                                uint &total_nodes,
+                                vector<Primitive> &ordered_prims);
+
+  /* returns the union of two orientation bounds and returns the result */
+  Orientation cone_union(const Orientation &a, const Orientation &b);
+
+  /* returns the union of several orientation bounds */
+  Orientation combine_bounding_cones(const vector<Orientation> &bcones);
+
+  /* calculates the cone measure in the surface area orientation heuristic */
+  float calculate_cone_measure(const Orientation &bcone);
+
+  /* takes a node and the lights contained in it as input and returns a way to
+   * split the node into two child nodes. This is done as follows:
+   * 1. A bounding box of all lights centroid is constructed
+   * 2. A set of candidate splits(proposed left and right child nodes) are
+   *    created.
+   *    - This is done by partitioning the bounding box into two regions.
+   *      All lights in the same region belongs to the same child node. This
+   *      is done for several partions of the bounding box.
+   * 3. Each such candidate is evaluated using the Surface Area Orientation
+   *    Heuristic(SAOH).
+   * 4. The candidate split with the minimum cost(heuristic) is returned */
+  void split_saoh(const BoundBox &centroid_bbox,
+                  const vector<BVHPrimitiveInfo> &build_data,
+                  const int start,
+                  const int end,
+                  const int num_buckets,
+                  const float node_M_Omega,
+                  const BoundBox &node_bbox,
+                  float &min_cost,
+                  int &min_dim,
+                  int &min_bucket);
+
+  /* this method performs step 3 above. */
+  int flattenBVHTree(const BVHBuildNode &node, int &offset);
+
+  /* contains all the lights in the scene. when the constructor has finished,
+   * these will be ordered such that lights belonging to the same node will be
+   * next to each other in this array. */
+  vector<Primitive> primitives;
+
+  /* the maximum number of allowed lights in each leaf node */
+  uint max_lights_in_node;
+
+  /* pointer to the scene. this is used to access the objects, the lights and
+   * the shader manager of the scene.*/
+  Scene *scene;
+
+  /* the nodes of the light hierarchy */
+  vector<CompactNode> nodes;
+};
+
+CCL_NAMESPACE_END
+
+#endif  // __LIGHT_TREE_H__

intern/cycles/render/mesh.h

@@ -86,6 +86,25 @@ class Mesh : public Geometry {
     return tri;
   }
 
+  float compute_triangle_area(size_t i, const Transform &tfm) const
+  {
+    Mesh::Triangle t = get_triangle(i);
+    if (!t.valid(&verts[0])) {
+      return 0.0f;
+    }
+    float3 p1 = verts[t.v[0]];
+    float3 p2 = verts[t.v[1]];
+    float3 p3 = verts[t.v[2]];
+
+    if (!transform_applied) {
+      p1 = transform_point(&tfm, p1);
+      p2 = transform_point(&tfm, p2);
+      p3 = transform_point(&tfm, p3);
+    }
+
+    return triangle_area(p1, p2, p3);
+  }
+
   size_t num_triangles() const
   {
     return triangles.size() / 3;

intern/cycles/render/mesh_volume.cpp

@@ -159,7 +159,7 @@ class VolumeMeshBuilder {
   bool empty_grid() const;
 
 #ifdef WITH_OPENVDB
-  template <typename GridType>
+  template<typename GridType>
   void merge_grid(openvdb::GridBase::ConstPtr grid, bool do_clipping, float volume_clipping)
   {
     typename GridType::ConstPtr typed_grid = openvdb::gridConstPtrCast<GridType>(grid);
@@ -189,7 +189,9 @@ VolumeMeshBuilder::VolumeMeshBuilder()
 }
 
 #ifdef WITH_OPENVDB
-void VolumeMeshBuilder::add_grid(openvdb::GridBase::ConstPtr grid, bool do_clipping, float volume_clipping)
+void VolumeMeshBuilder::add_grid(openvdb::GridBase::ConstPtr grid,
+                                 bool do_clipping,
+                                 float volume_clipping)
 {
   /* set the transform of our grid from the first one */
   if (first_grid) {
@@ -415,22 +417,38 @@ static openvdb::GridBase::ConstPtr openvdb_grid_from_device_texture(device_textu
   typename GridType::Ptr sparse = GridType::create(ValueType(0.0f));
   openvdb::tools::copyFromDense(dense, *sparse, ValueType(volume_clipping));
 
-  /* copyFromDense will remove any leaf node that contains constant data and replace it with a tile,
-   * however, we need to preserve the leaves in order to generate the mesh, so revoxelize the leaves
-   * that were pruned. This should not affect areas that were skipped due to the volume_clipping parameter. */
+  /* copyFromDense will remove any leaf node that contains constant data and replace it with a
+   * tile, however, we need to preserve the leaves in order to generate the mesh, so revoxelize the
+   * leaves that were pruned. This should not affect areas that were skipped due to the
+   * volume_clipping parameter. */
   sparse->tree().voxelizeActiveTiles();
 
   /* Compute index to world matrix. */
-  float3 voxel_size = make_float3(1.0f / image_memory->data_width, 1.0f / image_memory->data_height, 1.0f / image_memory->data_depth);
+  float3 voxel_size = make_float3(1.0f / image_memory->data_width,
+                                  1.0f / image_memory->data_height,
+                                  1.0f / image_memory->data_depth);
 
   transform_3d = transform_inverse(transform_3d);
 
-  openvdb::Mat4R index_to_world_mat((double)(voxel_size.x * transform_3d[0][0]), 0.0, 0.0, 0.0,
-                           0.0, (double)(voxel_size.y * transform_3d[1][1]), 0.0, 0.0,
-                           0.0, 0.0, (double)(voxel_size.z * transform_3d[2][2]), 0.0,
-                           (double)transform_3d[0][3], (double)transform_3d[1][3], (double)transform_3d[2][3], 1.0);
+  openvdb::Mat4R index_to_world_mat((double)(voxel_size.x * transform_3d[0][0]),
+                                    0.0,
+                                    0.0,
+                                    0.0,
+                                    0.0,
+                                    (double)(voxel_size.y * transform_3d[1][1]),
+                                    0.0,
+                                    0.0,
+                                    0.0,
+                                    0.0,
+                                    (double)(voxel_size.z * transform_3d[2][2]),
+                                    0.0,
+                                    (double)transform_3d[0][3],
+                                    (double)transform_3d[1][3],
+                                    (double)transform_3d[2][3],
+                                    1.0);
 
-  openvdb::math::Transform::Ptr index_to_world_tfm = openvdb::math::Transform::createLinearTransform(index_to_world_mat);
+  openvdb::math::Transform::Ptr index_to_world_tfm =
+      openvdb::math::Transform::createLinearTransform(index_to_world_mat);
 
   sparse->setTransform(index_to_world_tfm);
 
@@ -472,19 +490,16 @@ void GeometryManager::create_volume_mesh(Mesh *mesh, Progress &progress)
       device_texture *image_memory = handle.image_memory();
 
       if (image_memory->data_elements == 1) {
-        grid = openvdb_grid_from_device_texture<openvdb::FloatGrid>(image_memory,
-                                                                    mesh->volume_clipping,
-                                                                    handle.metadata().transform_3d);
+        grid = openvdb_grid_from_device_texture<openvdb::FloatGrid>(
+            image_memory, mesh->volume_clipping, handle.metadata().transform_3d);
       }
       else if (image_memory->data_elements == 3) {
-        grid = openvdb_grid_from_device_texture<openvdb::Vec3fGrid>(image_memory,
-                                                                    mesh->volume_clipping,
-                                                                    handle.metadata().transform_3d);
+        grid = openvdb_grid_from_device_texture<openvdb::Vec3fGrid>(
+            image_memory, mesh->volume_clipping, handle.metadata().transform_3d);
       }
       else if (image_memory->data_elements == 4) {
-        grid = openvdb_grid_from_device_texture<openvdb::Vec4fGrid>(image_memory,
-                                                                    mesh->volume_clipping,
-                                                                    handle.metadata().transform_3d);
+        grid = openvdb_grid_from_device_texture<openvdb::Vec4fGrid>(
+            image_memory, mesh->volume_clipping, handle.metadata().transform_3d);
       }
     }
 

intern/cycles/render/scene.cpp

@@ -79,7 +79,15 @@ DeviceScene::DeviceScene(Device *device)
       shaders(device, "__shaders", MEM_GLOBAL),
       lookup_table(device, "__lookup_table", MEM_GLOBAL),
       sample_pattern_lut(device, "__sample_pattern_lut", MEM_GLOBAL),
-      ies_lights(device, "__ies", MEM_GLOBAL)
+      ies_lights(device, "__ies", MEM_GLOBAL),
+      light_tree_nodes(device, "__light_tree_nodes", MEM_GLOBAL),
+      light_distribution_to_node(device, "__light_distribution_to_node", MEM_GLOBAL),
+      lamp_to_distribution(device, "__lamp_to_distribution", MEM_GLOBAL),
+      triangle_to_distribution(device, "__triangle_to_distribution", MEM_GLOBAL),
+      light_group_sample_cdf(device, "__light_group_sample_cdf", MEM_GLOBAL),
+      light_group_sample_prob(device, "__light_group_sample_prob", MEM_GLOBAL),
+      leaf_to_first_emitter(device, "__leaf_to_first_emitter", MEM_GLOBAL),
+      light_tree_leaf_emitters(device, "__light_tree_leaf_emitters", MEM_GLOBAL)
 {
   memset((void *)&data, 0, sizeof(data));
 }

intern/cycles/render/scene.h

@@ -108,6 +108,14 @@ class DeviceScene {
   device_vector<KernelLight> lights;
   device_vector<float2> light_background_marginal_cdf;
   device_vector<float2> light_background_conditional_cdf;
+  device_vector<float4> light_tree_nodes;
+  device_vector<uint> light_distribution_to_node;
+  device_vector<uint> lamp_to_distribution;
+  device_vector<uint> triangle_to_distribution;
+  device_vector<float> light_group_sample_cdf;
+  device_vector<float> light_group_sample_prob;
+  device_vector<int> leaf_to_first_emitter;
+  device_vector<float4> light_tree_leaf_emitters;
 
   /* particles */
   device_vector<KernelParticle> particles;

source/blender/python/mathutils/mathutils_bvhtree.c

@@ -549,8 +549,7 @@ static bool py_bvhtree_overlap_cb(void *userdata, int index_a, int index_b, int
     }
   }
 
-  return (isect_tri_tri_v3(
-              UNPACK3(tri_a_co), UNPACK3(tri_b_co), ix_pair[0], ix_pair[1]) &&
+  return (isect_tri_tri_v3(UNPACK3(tri_a_co), UNPACK3(tri_b_co), ix_pair[0], ix_pair[1]) &&
           ((verts_shared == 0) || (len_squared_v3v3(ix_pair[0], ix_pair[1]) > data->epsilon)));
 }