diff --git a/.gitignore b/.gitignore
index acdb42f8..e051e093 100644
--- a/.gitignore
+++ b/.gitignore
@@ -31,7 +31,10 @@ __pycache__
 
 # Results
 results/.ipynb_checkpoints
-results/estimator_comparison_simple
+results/estimator_comparison_vmf_conductor
+results/estimator_comparison_vmf_dielectric
+results/estimator_comparison_envmap_conductor
+results/estimator_comparison_envmap_dielectric
 results/estimator_comparison_optimization
 results/estimator_comparison_attached_reparam
 results/estimator_comparison_veach
diff --git a/README.md b/README.md
index 89a5c6e8..3dd974e7 100644
--- a/README.md
+++ b/README.md
@@ -25,9 +25,11 @@ The `results` directory contains a set of Jupyter notebooks that reproduce some
       - Parameter-dependent discontinuities from moving geometry.
       - Parameter-dependent discontinuities from static geometry in presence of attached sampling strategies.
 
-3. **Estimator comparison (simple)**
+3. **Estimator comparison (vMF conductor)**
 
-   Forward and reverse mode differentiation of a simple microfacet material's roughness parameter. Compares a long list of differential estimators:
+   Forward and reverse mode differentiation of a rough conductor roughness
+   parameter in a simple test scene with a vMF emitter.
+   Compares a long list of differential estimators:
       - Detached emitter sampling
       - Detached BSDF sampling
       - Attached BSDF sampling
@@ -40,7 +42,19 @@ The `results` directory contains a set of Jupyter notebooks that reproduce some
       - Detached BSDF sampling using the differential microfacet sampling strategy
       - Detached MIS weights, detached differential microfacet BSDF sampling, detached emitter sampling
 
-4. **Estimator comparison (roughness optimization)**
+4. **Estimator comparison (vMF dielectric)**
+
+   Same as (3) but with a rough dielectric material instead.
+
+5. **Estimator comparison (envmap conductor)**
+
+   Same as (3) but with an environment emitter instead of the vMF.
+
+6. **Estimator comparison (envmap dielectric)**
+
+   Same as (5) but with a rough dielectric material instead.
+
+7. **Estimator comparison (roughness optimization)**
    
    A simple gradient based optimization recovering a roughness texture from a given reference rendering. Compares three different BSDF sampling techniques under two different illumination conditions. (Similar to Figure 12 from the paper):
    
@@ -48,14 +62,14 @@ The `results` directory contains a set of Jupyter notebooks that reproduce some
       - Attached BSDF sampling
       - Detached BSDF sampling using the differential microfacet sampling strategy
 
-5. **Estimator comparison (Veach)**
+8. **Estimator comparison (Veach)**
 
    Reproduces the (forward mode) gradient images of the classical Veach test scene from Figure 11 in the paper. Compared to the simpler test scenes above, this scene contains (static) visibility discontinuities that causes naïve attached estimators to be biased. In addition to the previously listed estimators, two additional ones involving reparameterized attached sampling are used:
    
       - Reparameterized attached BSDF sampling
       - Attached MIS weights, reparameterized attached BSDF sampling, detached emitter sampling
 
-6. **Estimator comparison (attached reparameterized)**
+9. **Estimator comparison (attached reparameterized)**
 
    Again illustrates the bias of attached BSDF sampling in presence of static visibility discontinuities, this time visualizing parameter gradient textures produced by reverse mode differentiation. Based on Figure 10 of the paper.
    
diff --git a/results/Estimator comparison (envmap conductor).ipynb b/results/Estimator comparison (envmap conductor).ipynb
new file mode 100644
index 00000000..d5be83f8
--- /dev/null
+++ b/results/Estimator comparison (envmap conductor).ipynb	
@@ -0,0 +1,407 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15f8e794",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:24:40.951530Z",
+     "start_time": "2022-08-14T04:24:40.667076Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import mitsuba as mi\n",
+    "mi.set_variant(\"llvm_ad_rgb\")\n",
+    "import drjit as dr\n",
+    "\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+    "import cmap_diff\n",
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "%matplotlib inline\n",
+    "\n",
+    "import os\n",
+    "base_dir = 'estimator_comparison_envmap_conductor'\n",
+    "if not os.path.exists(base_dir):\n",
+    "    os.makedirs(base_dir)\n",
+    "    \n",
+    "mi.Thread.thread().logger().set_log_level(mi.LogLevel.Warn)\n",
+    "\n",
+    "def produce_plots(methods, method_names, directory, gradients=False):\n",
+    "    # Suitable values to compare gradients computed with all different methods (determined manually)\n",
+    "    scales = np.array([20, 10, 5, 2])\n",
+    "\n",
+    "    for method, method_name in zip(methods, method_names):\n",
+    "        fig, axes = plt.subplots(ncols=len(alphas), nrows=1, figsize=(12,4))\n",
+    "        for idx_alpha, _ in enumerate(alphas):\n",
+    "            name = \"a_{:02d}\".format(idx_alpha)\n",
+    "            path = '{}/{}/{}/{}.exr'.format(base_dir, directory, method, name)\n",
+    "            ax = axes[idx_alpha]\n",
+    "            ax.set_xticks([]); ax.set_yticks([])\n",
+    "            data = np.array(mi.Bitmap(path)).astype('float32')\n",
+    "            if gradients:\n",
+    "                vminmax = scales[idx_alpha]\n",
+    "                data_ = data[:,:,0] if len(data.shape) == 3 else data\n",
+    "                im = ax.imshow(data_, cmap='diff', vmin=-vminmax, vmax=+vminmax)\n",
+    "            else:\n",
+    "                data = np.clip(data**(1/2.2), 0.0, 1.0) # Crude gamma correction\n",
+    "                im = ax.imshow(data, cmap='gray')\n",
+    "\n",
+    "        fig.suptitle(method_name, y=0.98, size=13, weight='bold')\n",
+    "        fig.text(0.125, 0.85, 'Surface roughness --->',\n",
+    "                 ha='left', size=14, weight='bold')\n",
+    "        for idx_alpha, alpha in enumerate(alphas):\n",
+    "            axes[idx_alpha].set_title(r\"$\\alpha={}$\".format(alpha), size=14)\n",
+    "        outname = '{}/{}/{}.jpg'.format(base_dir, directory, method)\n",
+    "        plt.savefig(outname, dpi=150, pad_inches=0.1, bbox_inches='tight')\n",
+    "        plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f00c2237",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:24:40.974613Z",
+     "start_time": "2022-08-14T04:24:40.964613Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# We look at a the simple test scene of a rough conductor plane, lit by an environment map.\n",
+    "# This allows us to judge the effectiveness of the different MC estimators at various\n",
+    "# combinations of emitter concentration and surface roughness parameters.\n",
+    "# No discontinuities are present in this scene.\n",
+    "\n",
+    "# Scene parameters varied in each configuration\n",
+    "alpha_key = 'plane.bsdf.alpha.data'\n",
+    "\n",
+    "# Values to consider in the comparison\n",
+    "alphas = [0.02, 0.05, 0.1, 0.2]\n",
+    "\n",
+    "# Sampels per pixel for all renderings\n",
+    "spp = 128\n",
+    "\n",
+    "primal_methods = [\n",
+    "    'primal_bs',\n",
+    "    'primal_es',\n",
+    "    'primal_mis',\n",
+    "]\n",
+    "primal_method_names = [\n",
+    "    'Primal BSDF sampling',\n",
+    "    'Primal emitter sampling',\n",
+    "    'Primal MIS (BSDF + emitter sampling)',\n",
+    "]\n",
+    "\n",
+    "diff_methods = [\n",
+    "    'es_detached',\n",
+    "    \n",
+    "    'bs_detached',\n",
+    "    'bs_attached',\n",
+    "    \n",
+    "    'mis_detached_detached',\n",
+    "    'mis_attached_attached',\n",
+    "    'mis_detached_attached',\n",
+    "    'mis_attached_detached',\n",
+    "    \n",
+    "    'bs_detached_diff',\n",
+    "    'mis_detached_detached_diff',\n",
+    "]\n",
+    "diff_method_names = [\n",
+    "    'Detached emitter sampling',\n",
+    "    \n",
+    "    'Detached BSDF sampling',\n",
+    "    'Attached BSDF sampling',\n",
+    "    \n",
+    "    'Detached MIS weights, detached BSDF sampling, detached emitter sampling',\n",
+    "    'Attached MIS weights, attached BSDF sampling, detached emitter sampling',\n",
+    "    'Detached MIS weights, attached BSDF sampling, detached emitter sampling (BIASED!)',\n",
+    "    'Attached MIS weights, detached BSDF sampling, detached emitter sampling',\n",
+    "    \n",
+    "    'Detached diff. BSDF sampling',\n",
+    "    'Detached MIS weights, detached diff. BSDF sampling, detached emitter sampling',\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6ef9f2d6",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:24:46.121144Z",
+     "start_time": "2022-08-14T04:24:40.981562Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "# #####################\n",
+    "# # Primal estimators #\n",
+    "# #####################\n",
+    "\n",
+    "for method, method_name in zip(primal_methods, primal_method_names):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/primal/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "        \n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method, 'hide_emitters': True})\n",
+    "        \n",
+    "    scene = mi.load_file('scenes/conductor_plane_envmap.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_alpha, alpha in enumerate(alphas):\n",
+    "        print(\"  - {}/{}\".format(idx_alpha + 1, len(alphas)), end='\\r')\n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "        params.update()\n",
+    "\n",
+    "        image = mi.render(scene, params, integrator=integrator, seed=0, spp=spp)\n",
+    "        outname = \"{}/a_{:02d}.exr\".format(method_dir, idx_alpha)\n",
+    "        mi.util.convert_to_bitmap(image, uint8_srgb=False).write(outname)\n",
+    "    print(\"\")\n",
+    "    \n",
+    "produce_plots(primal_methods, primal_method_names, 'primal')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96bf3d5f",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:24:58.561975Z",
+     "start_time": "2022-08-14T04:24:46.122511Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "######################################\n",
+    "# Finite differences (FORWARD mode)  #\n",
+    "######################################\n",
+    "# Produce a gradient image as a result.\n",
+    "\n",
+    "eps = 1e-4\n",
+    "\n",
+    "for method, method_name in zip(['fd'], ['Finite differences']):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_forward/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "\n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': 'primal_mis', 'hide_emitters': True})\n",
+    "    \n",
+    "    scene = mi.load_file('scenes/conductor_plane_envmap.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_alpha, alpha in enumerate(alphas):\n",
+    "        print(\"  - {}/{}\".format(idx_alpha + 1, len(alphas)), end='\\r') \n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "        params.update()\n",
+    "\n",
+    "        # Render first image\n",
+    "        image_0 = mi.render(scene, params, integrator=integrator, seed=0, spp=512)\n",
+    "\n",
+    "        # Apply FD\n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*(alpha + eps))\n",
+    "        params.update()\n",
+    "\n",
+    "        # Render second image\n",
+    "        image_1 = mi.render(scene, params, integrator=integrator, seed=0, spp=512)\n",
+    "\n",
+    "        # Save output gradient image\n",
+    "        image_grad = (image_1 - image_0) / eps\n",
+    "        outname = \"{}/a_{:02d}.exr\".format(method_dir, idx_alpha)\n",
+    "        mi.util.convert_to_bitmap(image_grad, uint8_srgb=False).write(outname)\n",
+    "            \n",
+    "produce_plots(['fd'], ['Finite differences'], 'gradients_forward', gradients=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fbec4a0d",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:25:41.377382Z",
+     "start_time": "2022-08-14T04:24:58.562772Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "######################################\n",
+    "# Gradient estimators (FORWARD mode) #\n",
+    "######################################\n",
+    "# These propagate a scalar input gradient (in the surface roughness) to the output pixels,\n",
+    "# i.e. produce a gradient image as a result.\n",
+    "\n",
+    "for method, method_name in zip(diff_methods, diff_method_names):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_forward/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "        \n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method, 'hide_emitters': True})\n",
+    "        \n",
+    "    scene = mi.load_file('scenes/conductor_plane_envmap.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_alpha, alpha in enumerate(alphas):\n",
+    "        print(\"  - {}/{}\".format(idx_alpha + 1, len(alphas)), end='\\r') \n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "\n",
+    "        # Diff. input parameter\n",
+    "        pi = mi.Float(0.0)\n",
+    "        dr.enable_grad(pi)\n",
+    "        dr.set_grad(pi, 1.0)\n",
+    "        params[alpha_key] += pi\n",
+    "        params.update()\n",
+    "\n",
+    "        image_grad = mi.Float(0.0)\n",
+    "        if 'diff' in method:\n",
+    "            # Differential sampling strategy, use antithetic sampling.\n",
+    "            # Note that we use the same seed twice, and use half the number of samples for each pass.\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp//2, antithetic_pass=False)\n",
+    "            dr.forward(pi)\n",
+    "            image_grad = dr.grad(image)\n",
+    "\n",
+    "            params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "            params[alpha_key] += pi\n",
+    "            params.update()\n",
+    "\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp//2, antithetic_pass=True)\n",
+    "            dr.forward(pi)\n",
+    "            image_grad += dr.grad(image)\n",
+    "\n",
+    "            # Average both passes\n",
+    "            image_grad *= 0.5\n",
+    "        else:\n",
+    "            # Produce differentiable rendering\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp)\n",
+    "            # And propagate derivatives forwards through it\n",
+    "            dr.forward(pi)\n",
+    "            image_grad = dr.grad(image)\n",
+    "\n",
+    "        # Save output gradient image\n",
+    "        outname = \"{}/a_{:02d}.exr\".format(method_dir, idx_alpha)\n",
+    "        mi.util.convert_to_bitmap(image_grad, uint8_srgb=False).write(outname)\n",
+    "            \n",
+    "produce_plots(diff_methods, diff_method_names, 'gradients_forward', gradients=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "71be7f79",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:26:25.537236Z",
+     "start_time": "2022-08-14T04:25:41.378161Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "######################################\n",
+    "# Gradient estimators (REVERSE mode) #\n",
+    "######################################\n",
+    "# These propagate an output image gradient to the (textured) scene input parameters,\n",
+    "# i.e. the result are texture gradient images.\n",
+    "\n",
+    "for method, method_name in zip(diff_methods, diff_method_names):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_reverse/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "        \n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method, 'hide_emitters': True})\n",
+    "        \n",
+    "    scene = mi.load_file('scenes/conductor_plane_envmap.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    img_shape     = scene.sensors()[0].film().crop_size()\n",
+    "    \n",
+    "    for idx_alpha, alpha in enumerate(alphas):\n",
+    "        print(\"  - {}/{}\".format(idx_alpha + 1, len(alphas)), end='\\r') \n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "        params.update()\n",
+    "\n",
+    "        dr.enable_grad(params[alpha_key])\n",
+    "\n",
+    "        grad_backward = mi.Float(0.0)\n",
+    "        if 'diff' in method:\n",
+    "            # Differential sampling strategy, use antithetic sampling.\n",
+    "            # Note that we use the same seed twice, and use half the number of samples for each pass.\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp//2, antithetic_pass=False)\n",
+    "            dr.backward(image)\n",
+    "\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp//2, antithetic_pass=True)\n",
+    "            dr.backward(image)\n",
+    "\n",
+    "            # Average accumulated gradients from both passes\n",
+    "            grad_backward = 0.5 * dr.grad(params[alpha_key])\n",
+    "        else:\n",
+    "            # Produce differentiable rendering\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp)\n",
+    "\n",
+    "            # And propagate derivatives backwards through it\n",
+    "            dr.backward(image)\n",
+    "            grad_backward = dr.grad(params[alpha_key])\n",
+    "\n",
+    "        # Save output gradient image\n",
+    "        outname = \"{}/a_{:02d}.exr\".format(method_dir, idx_alpha)\n",
+    "        mi.util.convert_to_bitmap(grad_backward, uint8_srgb=False).write(outname)\n",
+    "    \n",
+    "produce_plots(diff_methods, diff_method_names, 'gradients_reverse', gradients=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a89d5169",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/results/Estimator comparison (envmap dielectric).ipynb b/results/Estimator comparison (envmap dielectric).ipynb
new file mode 100644
index 00000000..68365631
--- /dev/null
+++ b/results/Estimator comparison (envmap dielectric).ipynb	
@@ -0,0 +1,407 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15f8e794",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:45:06.837104Z",
+     "start_time": "2022-08-14T04:45:06.806412Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import mitsuba as mi\n",
+    "mi.set_variant(\"llvm_ad_rgb\")\n",
+    "import drjit as dr\n",
+    "\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+    "import cmap_diff\n",
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "%matplotlib inline\n",
+    "\n",
+    "import os\n",
+    "base_dir = 'estimator_comparison_envmap_dielectric'\n",
+    "if not os.path.exists(base_dir):\n",
+    "    os.makedirs(base_dir)\n",
+    "    \n",
+    "mi.Thread.thread().logger().set_log_level(mi.LogLevel.Warn)\n",
+    "\n",
+    "def produce_plots(methods, method_names, directory, gradients=False):\n",
+    "    # Suitable values to compare gradients computed with all different methods (determined manually)\n",
+    "    scales = np.array([1, 0.5, 0.2, 0.1])\n",
+    "\n",
+    "    for method, method_name in zip(methods, method_names):\n",
+    "        fig, axes = plt.subplots(ncols=len(alphas), nrows=1, figsize=(12,4))\n",
+    "        for idx_alpha, _ in enumerate(alphas):\n",
+    "            name = \"a_{:02d}\".format(idx_alpha)\n",
+    "            path = '{}/{}/{}/{}.exr'.format(base_dir, directory, method, name)\n",
+    "            ax = axes[idx_alpha]\n",
+    "            ax.set_xticks([]); ax.set_yticks([])\n",
+    "            data = np.array(mi.Bitmap(path)).astype('float32')\n",
+    "            if gradients:\n",
+    "                vminmax = scales[idx_alpha]\n",
+    "                data_ = data[:,:,0] if len(data.shape) == 3 else data\n",
+    "                im = ax.imshow(data_, cmap='diff', vmin=-vminmax, vmax=+vminmax)\n",
+    "            else:\n",
+    "                data = np.clip(data**(1/2.2), 0.0, 1.0) # Crude gamma correction\n",
+    "                im = ax.imshow(data, cmap='gray')\n",
+    "\n",
+    "        fig.suptitle(method_name, y=0.98, size=13, weight='bold')\n",
+    "        fig.text(0.125, 0.85, 'Surface roughness --->',\n",
+    "                 ha='left', size=14, weight='bold')\n",
+    "        for idx_alpha, alpha in enumerate(alphas):\n",
+    "            axes[idx_alpha].set_title(r\"$\\alpha={}$\".format(alpha), size=14)\n",
+    "        outname = '{}/{}/{}.jpg'.format(base_dir, directory, method)\n",
+    "        plt.savefig(outname, dpi=150, pad_inches=0.1, bbox_inches='tight')\n",
+    "        plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f00c2237",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:43:47.252023Z",
+     "start_time": "2022-08-14T04:43:47.241600Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# We look at a the simple test scene of a rough conductor plane, lit by an environment map.\n",
+    "# This allows us to judge the effectiveness of the different MC estimators at various\n",
+    "# combinations of emitter concentration and surface roughness parameters.\n",
+    "# No discontinuities are present in this scene.\n",
+    "\n",
+    "# Scene parameters varied in each configuration\n",
+    "alpha_key = 'plane.bsdf.alpha.data'\n",
+    "\n",
+    "# Values to consider in the comparison\n",
+    "alphas = [0.02, 0.05, 0.1, 0.2]\n",
+    "\n",
+    "# Sampels per pixel for all renderings\n",
+    "spp = 128\n",
+    "\n",
+    "primal_methods = [\n",
+    "    'primal_bs',\n",
+    "    'primal_es',\n",
+    "    'primal_mis',\n",
+    "]\n",
+    "primal_method_names = [\n",
+    "    'Primal BSDF sampling',\n",
+    "    'Primal emitter sampling',\n",
+    "    'Primal MIS (BSDF + emitter sampling)',\n",
+    "]\n",
+    "\n",
+    "diff_methods = [\n",
+    "    'es_detached',\n",
+    "    \n",
+    "    'bs_detached',\n",
+    "    'bs_attached',\n",
+    "    \n",
+    "    'mis_detached_detached',\n",
+    "    'mis_attached_attached',\n",
+    "    'mis_detached_attached',\n",
+    "    'mis_attached_detached',\n",
+    "    \n",
+    "    'bs_detached_diff',\n",
+    "    'mis_detached_detached_diff',\n",
+    "]\n",
+    "diff_method_names = [\n",
+    "    'Detached emitter sampling',\n",
+    "    \n",
+    "    'Detached BSDF sampling',\n",
+    "    'Attached BSDF sampling',\n",
+    "    \n",
+    "    'Detached MIS weights, detached BSDF sampling, detached emitter sampling',\n",
+    "    'Attached MIS weights, attached BSDF sampling, detached emitter sampling',\n",
+    "    'Detached MIS weights, attached BSDF sampling, detached emitter sampling (BIASED!)',\n",
+    "    'Attached MIS weights, detached BSDF sampling, detached emitter sampling',\n",
+    "    \n",
+    "    'Detached diff. BSDF sampling',\n",
+    "    'Detached MIS weights, detached diff. BSDF sampling, detached emitter sampling',\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6ef9f2d6",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:43:58.615363Z",
+     "start_time": "2022-08-14T04:43:53.307567Z"
+    },
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "# #####################\n",
+    "# # Primal estimators #\n",
+    "# #####################\n",
+    "\n",
+    "for method, method_name in zip(primal_methods, primal_method_names):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/primal/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "        \n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method, 'hide_emitters': True})\n",
+    "        \n",
+    "    scene = mi.load_file('scenes/dielectric_plane_envmap.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_alpha, alpha in enumerate(alphas):\n",
+    "        print(\"  - {}/{}\".format(idx_alpha + 1, len(alphas)), end='\\r')\n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "        params.update()\n",
+    "\n",
+    "        image = mi.render(scene, params, integrator=integrator, seed=0, spp=spp)\n",
+    "        outname = \"{}/a_{:02d}.exr\".format(method_dir, idx_alpha)\n",
+    "        mi.util.convert_to_bitmap(image, uint8_srgb=False).write(outname)\n",
+    "    print(\"\")\n",
+    "    \n",
+    "produce_plots(primal_methods, primal_method_names, 'primal')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96bf3d5f",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:46:10.852004Z",
+     "start_time": "2022-08-14T04:45:58.401590Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "######################################\n",
+    "# Finite differences (FORWARD mode)  #\n",
+    "######################################\n",
+    "# Produce a gradient image as a result.\n",
+    "\n",
+    "eps = 1e-4\n",
+    "\n",
+    "for method, method_name in zip(['fd'], ['Finite differences']):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_forward/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "\n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': 'primal_mis', 'hide_emitters': True})\n",
+    "    \n",
+    "    scene = mi.load_file('scenes/dielectric_plane_envmap.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_alpha, alpha in enumerate(alphas):\n",
+    "        print(\"  - {}/{}\".format(idx_alpha + 1, len(alphas)), end='\\r') \n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "        params.update()\n",
+    "\n",
+    "        # Render first image\n",
+    "        image_0 = mi.render(scene, params, integrator=integrator, seed=0, spp=512)\n",
+    "\n",
+    "        # Apply FD\n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*(alpha + eps))\n",
+    "        params.update()\n",
+    "\n",
+    "        # Render second image\n",
+    "        image_1 = mi.render(scene, params, integrator=integrator, seed=0, spp=512)\n",
+    "\n",
+    "        # Save output gradient image\n",
+    "        image_grad = (image_1 - image_0) / eps\n",
+    "        outname = \"{}/a_{:02d}.exr\".format(method_dir, idx_alpha)\n",
+    "        mi.util.convert_to_bitmap(image_grad, uint8_srgb=False).write(outname)\n",
+    "            \n",
+    "produce_plots(['fd'], ['Finite differences'], 'gradients_forward', gradients=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fbec4a0d",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:45:58.400654Z",
+     "start_time": "2022-08-14T04:45:16.191444Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "######################################\n",
+    "# Gradient estimators (FORWARD mode) #\n",
+    "######################################\n",
+    "# These propagate a scalar input gradient (in the surface roughness) to the output pixels,\n",
+    "# i.e. produce a gradient image as a result.\n",
+    "\n",
+    "for method, method_name in zip(diff_methods, diff_method_names):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_forward/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "        \n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method, 'hide_emitters': True})\n",
+    "        \n",
+    "    scene = mi.load_file('scenes/dielectric_plane_envmap.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_alpha, alpha in enumerate(alphas):\n",
+    "        print(\"  - {}/{}\".format(idx_alpha + 1, len(alphas)), end='\\r') \n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "\n",
+    "        # Diff. input parameter\n",
+    "        pi = mi.Float(0.0)\n",
+    "        dr.enable_grad(pi)\n",
+    "        dr.set_grad(pi, 1.0)\n",
+    "        params[alpha_key] += pi\n",
+    "        params.update()\n",
+    "\n",
+    "        image_grad = mi.Float(0.0)\n",
+    "        if 'diff' in method:\n",
+    "            # Differential sampling strategy, use antithetic sampling.\n",
+    "            # Note that we use the same seed twice, and use half the number of samples for each pass.\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp//2, antithetic_pass=False)\n",
+    "            dr.forward(pi)\n",
+    "            image_grad = dr.grad(image)\n",
+    "\n",
+    "            params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "            params[alpha_key] += pi\n",
+    "            params.update()\n",
+    "\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp//2, antithetic_pass=True)\n",
+    "            dr.forward(pi)\n",
+    "            image_grad += dr.grad(image)\n",
+    "\n",
+    "            # Average both passes\n",
+    "            image_grad *= 0.5\n",
+    "        else:\n",
+    "            # Produce differentiable rendering\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp)\n",
+    "            # And propagate derivatives forwards through it\n",
+    "            dr.forward(pi)\n",
+    "            image_grad = dr.grad(image)\n",
+    "\n",
+    "        # Save output gradient image\n",
+    "        outname = \"{}/a_{:02d}.exr\".format(method_dir, idx_alpha)\n",
+    "        mi.util.convert_to_bitmap(image_grad, uint8_srgb=False).write(outname)\n",
+    "            \n",
+    "produce_plots(diff_methods, diff_method_names, 'gradients_forward', gradients=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "71be7f79",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:48:25.111654Z",
+     "start_time": "2022-08-14T04:47:39.901103Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "######################################\n",
+    "# Gradient estimators (REVERSE mode) #\n",
+    "######################################\n",
+    "# These propagate an output image gradient to the (textured) scene input parameters,\n",
+    "# i.e. the result are texture gradient images.\n",
+    "\n",
+    "for method, method_name in zip(diff_methods, diff_method_names):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_reverse/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "        \n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method, 'hide_emitters': True})\n",
+    "        \n",
+    "    scene = mi.load_file('scenes/dielectric_plane_envmap.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    img_shape     = scene.sensors()[0].film().crop_size()\n",
+    "    \n",
+    "    for idx_alpha, alpha in enumerate(alphas):\n",
+    "        print(\"  - {}/{}\".format(idx_alpha + 1, len(alphas)), end='\\r') \n",
+    "        params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "        params.update()\n",
+    "\n",
+    "        dr.enable_grad(params[alpha_key])\n",
+    "\n",
+    "        grad_backward = mi.Float(0.0)\n",
+    "        if 'diff' in method:\n",
+    "            # Differential sampling strategy, use antithetic sampling.\n",
+    "            # Note that we use the same seed twice, and use half the number of samples for each pass.\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp//2, antithetic_pass=False)\n",
+    "            dr.backward(image)\n",
+    "\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp//2, antithetic_pass=True)\n",
+    "            dr.backward(image)\n",
+    "\n",
+    "            # Average accumulated gradients from both passes\n",
+    "            grad_backward = 0.5 * dr.grad(params[alpha_key])\n",
+    "        else:\n",
+    "            # Produce differentiable rendering\n",
+    "            image = mi.render(scene, params,\n",
+    "                              integrator=integrator, seed=0, spp=spp)\n",
+    "\n",
+    "            # And propagate derivatives backwards through it\n",
+    "            dr.backward(image)\n",
+    "            grad_backward = dr.grad(params[alpha_key])\n",
+    "\n",
+    "        # Save output gradient image\n",
+    "        outname = \"{}/a_{:02d}.exr\".format(method_dir, idx_alpha)\n",
+    "        mi.util.convert_to_bitmap(grad_backward, uint8_srgb=False).write(outname)\n",
+    "    \n",
+    "produce_plots(diff_methods, diff_method_names, 'gradients_reverse', gradients=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a8311bfb",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/results/Estimator comparison (simple).ipynb b/results/Estimator comparison (vMF conductor).ipynb
similarity index 82%
rename from results/Estimator comparison (simple).ipynb
rename to results/Estimator comparison (vMF conductor).ipynb
index 84d387f8..a290ea90 100644
--- a/results/Estimator comparison (simple).ipynb	
+++ b/results/Estimator comparison (vMF conductor).ipynb	
@@ -6,8 +6,8 @@
    "id": "15f8e794",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-07-27T08:38:36.252638Z",
-     "start_time": "2022-07-27T08:38:35.952633Z"
+     "end_time": "2022-08-14T04:26:21.156855Z",
+     "start_time": "2022-08-14T04:26:21.120135Z"
     }
    },
    "outputs": [],
@@ -24,7 +24,7 @@
     "%matplotlib inline\n",
     "\n",
     "import os\n",
-    "base_dir = 'estimator_comparison_simple'\n",
+    "base_dir = 'estimator_comparison_vmf_conductor'\n",
     "if not os.path.exists(base_dir):\n",
     "    os.makedirs(base_dir)\n",
     "    \n",
@@ -76,8 +76,8 @@
    "id": "f00c2237",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-07-27T08:38:39.036056Z",
-     "start_time": "2022-07-27T08:38:39.026138Z"
+     "end_time": "2022-08-14T04:26:24.976503Z",
+     "start_time": "2022-08-14T04:26:24.973718Z"
     }
    },
    "outputs": [],
@@ -145,9 +145,10 @@
    "id": "6ef9f2d6",
    "metadata": {
     "ExecuteTime": {
-     "start_time": "2022-07-27T08:19:06.642Z"
+     "end_time": "2022-08-14T04:26:44.187785Z",
+     "start_time": "2022-08-14T04:26:26.504582Z"
     },
-    "scrolled": false
+    "scrolled": true
    },
    "outputs": [],
    "source": [
@@ -163,7 +164,7 @@
     "        \n",
     "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method})\n",
     "        \n",
-    "    scene = mi.load_file('scenes/microfacet_plane.xml')\n",
+    "    scene = mi.load_file('scenes/conductor_plane_vmf.xml')\n",
     "    params = mi.traverse(scene)\n",
     "    params.keep([alpha_key, kappa_key])\n",
     "    texture_shape = params[alpha_key].shape\n",
@@ -187,9 +188,59 @@
    "cell_type": "code",
    "execution_count": null,
    "id": "96bf3d5f",
-   "metadata": {},
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:27:45.277800Z",
+     "start_time": "2022-08-14T04:27:02.782197Z"
+    },
+    "scrolled": true
+   },
    "outputs": [],
-   "source": []
+   "source": [
+    "######################################\n",
+    "# Finite differences (FORWARD mode)  #\n",
+    "######################################\n",
+    "# Produce a gradient image as a result.\n",
+    "\n",
+    "eps = 1e-4\n",
+    "\n",
+    "for method, method_name in zip(['fd'], ['Finite differences']):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_forward/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "\n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': 'primal_mis'})\n",
+    "    \n",
+    "    scene = mi.load_file('scenes/conductor_plane_vmf.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key, kappa_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_kappa, kappa in enumerate(kappas):\n",
+    "        for idx_alpha, alpha in enumerate(alphas):\n",
+    "            print(\"  - {}/{}\".format(idx_kappa*len(alphas) + idx_alpha + 1, len(alphas)*len(kappas)), end='\\r') \n",
+    "            params[kappa_key] = kappa\n",
+    "            params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "            params.update()\n",
+    "            \n",
+    "            # Render first image\n",
+    "            image_0 = mi.render(scene, params, integrator=integrator, seed=0, spp=512)\n",
+    "            \n",
+    "            # Apply FD\n",
+    "            params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*(alpha + eps))\n",
+    "            params.update()\n",
+    "            \n",
+    "            # Render second image\n",
+    "            image_1 = mi.render(scene, params, integrator=integrator, seed=0, spp=512)\n",
+    "            \n",
+    "            # Save output gradient image\n",
+    "            image_grad = (image_1 - image_0) / eps\n",
+    "            outname = \"{}/k_{:02d}_a_{:02d}.exr\".format(method_dir, idx_kappa, idx_alpha)\n",
+    "            mi.util.convert_to_bitmap(image_grad, uint8_srgb=False).write(outname)\n",
+    "            \n",
+    "produce_plots(['fd'], ['Finite differences'], 'gradients_forward', gradients=True)"
+   ]
   },
   {
    "cell_type": "code",
@@ -197,10 +248,10 @@
    "id": "fbec4a0d",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-07-27T08:42:24.872980Z",
-     "start_time": "2022-07-27T08:42:09.617042Z"
+     "end_time": "2022-08-14T04:30:08.101692Z",
+     "start_time": "2022-08-14T04:27:48.850431Z"
     },
-    "scrolled": false
+    "scrolled": true
    },
    "outputs": [],
    "source": [
@@ -218,7 +269,7 @@
     "        \n",
     "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method})\n",
     "        \n",
-    "    scene = mi.load_file('scenes/microfacet_plane.xml')\n",
+    "    scene = mi.load_file('scenes/conductor_plane_vmf.xml')\n",
     "    params = mi.traverse(scene)\n",
     "    params.keep([alpha_key, kappa_key])\n",
     "    texture_shape = params[alpha_key].shape\n",
@@ -278,10 +329,10 @@
    "id": "71be7f79",
    "metadata": {
     "ExecuteTime": {
-     "end_time": "2022-07-27T08:46:15.541746Z",
-     "start_time": "2022-07-27T08:45:59.858925Z"
+     "end_time": "2022-08-14T04:32:37.787971Z",
+     "start_time": "2022-08-14T04:30:08.102545Z"
     },
-    "scrolled": false
+    "scrolled": true
    },
    "outputs": [],
    "source": [
@@ -299,7 +350,7 @@
     "        \n",
     "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method})\n",
     "        \n",
-    "    scene = mi.load_file('scenes/microfacet_plane.xml')\n",
+    "    scene = mi.load_file('scenes/conductor_plane_vmf.xml')\n",
     "    params = mi.traverse(scene)\n",
     "    params.keep([alpha_key, kappa_key])\n",
     "    texture_shape = params[alpha_key].shape\n",
@@ -369,7 +420,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.12"
+   "version": "3.9.13"
   }
  },
  "nbformat": 4,
diff --git a/results/Estimator comparison (vMF dielectric).ipynb b/results/Estimator comparison (vMF dielectric).ipynb
new file mode 100644
index 00000000..9f4d89a5
--- /dev/null
+++ b/results/Estimator comparison (vMF dielectric).ipynb	
@@ -0,0 +1,428 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "15f8e794",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:33:02.981409Z",
+     "start_time": "2022-08-14T04:33:02.661114Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import mitsuba as mi\n",
+    "mi.set_variant(\"llvm_ad_rgb\")\n",
+    "import drjit as dr\n",
+    "\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+    "import cmap_diff\n",
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "%matplotlib inline\n",
+    "\n",
+    "import os\n",
+    "base_dir = 'estimator_comparison_vmf_dielectric'\n",
+    "if not os.path.exists(base_dir):\n",
+    "    os.makedirs(base_dir)\n",
+    "    \n",
+    "mi.Thread.thread().logger().set_log_level(mi.LogLevel.Warn)\n",
+    "\n",
+    "def produce_plots(methods, method_names, directory, gradients=False):\n",
+    "    # Suitable values to compare gradients computed with all different methods (determined manually)\n",
+    "    scales = np.array([ \n",
+    "        [300, 20, 2, 0.5],\n",
+    "        [50, 15, 3, 0.5],\n",
+    "        [5, 4, 2, 0.5],\n",
+    "        [0.2, 0.2, 0.2, 0.2],\n",
+    "    ])\n",
+    "\n",
+    "    for method, method_name in zip(methods, method_names):\n",
+    "        fig, axes = plt.subplots(ncols=len(alphas), nrows=len(kappas), figsize=(12,12))\n",
+    "        for idx_kappa, _ in enumerate(kappas):\n",
+    "            for idx_alpha, _ in enumerate(alphas):\n",
+    "                name = \"k_{:02d}_a_{:02d}\".format(idx_kappa, idx_alpha)\n",
+    "                path = '{}/{}/{}/{}.exr'.format(base_dir, directory, method, name)\n",
+    "                ax = axes[idx_kappa, idx_alpha]\n",
+    "                ax.set_xticks([]); ax.set_yticks([])\n",
+    "                data = np.array(mi.Bitmap(path)).astype('float32')\n",
+    "                if gradients:\n",
+    "                    vminmax = scales[idx_kappa, idx_alpha]\n",
+    "                    data_ = data[:,:,0] if len(data.shape) == 3 else data\n",
+    "                    im = ax.imshow(data_, cmap='diff', vmin=-vminmax, vmax=+vminmax)\n",
+    "                else:\n",
+    "                    data = np.clip(data**(1/2.2), 0.0, 1.0) # Crude gamma correction\n",
+    "                    im = ax.imshow(data, cmap='gray')\n",
+    "                \n",
+    "            fig.suptitle(method_name, y=0.98, size=14, weight='bold')\n",
+    "            fig.text(0.125, 0.92, 'Surface roughness --->',\n",
+    "                     ha='left', size=14, weight='bold')\n",
+    "            fig.text(0.075, 0.88, '<--- Illumination concentration',\n",
+    "                     va='top', rotation='vertical', size=14, weight='bold')\n",
+    "            for idx_kappa, kappa in enumerate(kappas):\n",
+    "                axes[idx_kappa, 0].set_ylabel(r\"$\\kappa={}$\".format(kappa), size=14)\n",
+    "            for idx_alpha, alpha in enumerate(alphas):\n",
+    "                axes[0, idx_alpha].set_title(r\"$\\alpha={}$\".format(alpha), size=14)\n",
+    "        outname = '{}/{}/{}.jpg'.format(base_dir, directory, method)\n",
+    "        plt.savefig(outname, dpi=150, pad_inches=0.1, bbox_inches='tight')\n",
+    "        plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f00c2237",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:33:03.002026Z",
+     "start_time": "2022-08-14T04:33:02.993859Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# We look at a the simple test scene of a rough dielectric plane, lit by a smooth emitter based on a \n",
+    "# vMF distribution. This allows us to judge the effectiveness of the different MC estimators at various\n",
+    "# combinations of emitter concentration and surface roughness parameters.\n",
+    "# No discontinuities are present in this scene.\n",
+    "\n",
+    "# Scene parameters varied in each configuration\n",
+    "alpha_key = 'plane.bsdf.alpha.data'\n",
+    "kappa_key = 'emitter.kappa'\n",
+    "\n",
+    "# Values to consider in the comparison\n",
+    "kappas = [100000, 1000, 100, 5]\n",
+    "alphas = [0.02, 0.05, 0.1, 0.2]\n",
+    "\n",
+    "# Sampels per pixel for all renderings\n",
+    "spp = 128\n",
+    "\n",
+    "primal_methods = [\n",
+    "    'primal_bs',\n",
+    "    'primal_es',\n",
+    "    'primal_mis',\n",
+    "]\n",
+    "primal_method_names = [\n",
+    "    'Primal BSDF sampling',\n",
+    "    'Primal emitter sampling',\n",
+    "    'Primal MIS (BSDF + emitter sampling)',\n",
+    "]\n",
+    "\n",
+    "diff_methods = [\n",
+    "    'es_detached',\n",
+    "    \n",
+    "    'bs_detached',\n",
+    "    'bs_attached',\n",
+    "    \n",
+    "    'mis_detached_detached',\n",
+    "    'mis_attached_attached',\n",
+    "    'mis_detached_attached',\n",
+    "    'mis_attached_detached',\n",
+    "    \n",
+    "    'bs_detached_diff',\n",
+    "    'mis_detached_detached_diff',\n",
+    "]\n",
+    "diff_method_names = [\n",
+    "    'Detached emitter sampling',\n",
+    "    \n",
+    "    'Detached BSDF sampling',\n",
+    "    'Attached BSDF sampling',\n",
+    "    \n",
+    "    'Detached MIS weights, detached BSDF sampling, detached emitter sampling',\n",
+    "    'Attached MIS weights, attached BSDF sampling, detached emitter sampling',\n",
+    "    'Detached MIS weights, attached BSDF sampling, detached emitter sampling (BIASED!)',\n",
+    "    'Attached MIS weights, detached BSDF sampling, detached emitter sampling',\n",
+    "    \n",
+    "    'Detached diff. BSDF sampling',\n",
+    "    'Detached MIS weights, detached diff. BSDF sampling, detached emitter sampling',\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6ef9f2d6",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:33:21.229558Z",
+     "start_time": "2022-08-14T04:33:03.003558Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "# #####################\n",
+    "# # Primal estimators #\n",
+    "# #####################\n",
+    "\n",
+    "for method, method_name in zip(primal_methods, primal_method_names):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/primal/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "        \n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method})\n",
+    "        \n",
+    "    scene = mi.load_file('scenes/dielectric_plane_vmf.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key, kappa_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_kappa, kappa in enumerate(kappas):\n",
+    "        for idx_alpha, alpha in enumerate(alphas):\n",
+    "            print(\"  - {}/{}\".format(idx_kappa*len(alphas) + idx_alpha + 1, len(alphas)*len(kappas)), end='\\r')\n",
+    "            params[kappa_key] = kappa\n",
+    "            params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "            params.update()\n",
+    "            \n",
+    "            image = mi.render(scene, params, integrator=integrator, seed=0, spp=spp)\n",
+    "            outname = \"{}/k_{:02d}_a_{:02d}.exr\".format(method_dir, idx_kappa, idx_alpha)\n",
+    "            mi.util.convert_to_bitmap(image, uint8_srgb=False).write(outname)\n",
+    "    print(\"\")\n",
+    "    \n",
+    "produce_plots(primal_methods, primal_method_names, 'primal')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96bf3d5f",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:36:55.461027Z",
+     "start_time": "2022-08-14T04:33:21.230657Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "######################################\n",
+    "# Finite differences (FORWARD mode)  #\n",
+    "######################################\n",
+    "# Produce a gradient image as a result.\n",
+    "\n",
+    "eps = 1e-4\n",
+    "\n",
+    "for method, method_name in zip(['fd'], ['Finite differences']):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_forward/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "\n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': 'primal_mis'})\n",
+    "    \n",
+    "    scene = mi.load_file('scenes/dielectric_plane_vmf.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key, kappa_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_kappa, kappa in enumerate(kappas):\n",
+    "        for idx_alpha, alpha in enumerate(alphas):\n",
+    "            print(\"  - {}/{}\".format(idx_kappa*len(alphas) + idx_alpha + 1, len(alphas)*len(kappas)), end='\\r') \n",
+    "            params[kappa_key] = kappa\n",
+    "            params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "            params.update()\n",
+    "            \n",
+    "            # Render first image\n",
+    "            image_0 = mi.render(scene, params, integrator=integrator, seed=0, spp=2048)\n",
+    "            \n",
+    "            # Apply FD\n",
+    "            params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*(alpha + eps))\n",
+    "            params.update()\n",
+    "            \n",
+    "            # Render second image\n",
+    "            image_1 = mi.render(scene, params, integrator=integrator, seed=0, spp=2048)\n",
+    "            \n",
+    "            # Save output gradient image\n",
+    "            image_grad = (image_1 - image_0) / eps\n",
+    "            outname = \"{}/k_{:02d}_a_{:02d}.exr\".format(method_dir, idx_kappa, idx_alpha)\n",
+    "            mi.util.convert_to_bitmap(image_grad, uint8_srgb=False).write(outname)\n",
+    "            \n",
+    "produce_plots(['fd'], ['Finite differences'], 'gradients_forward', gradients=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fbec4a0d",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:39:13.114526Z",
+     "start_time": "2022-08-14T04:36:55.461737Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "######################################\n",
+    "# Gradient estimators (FORWARD mode) #\n",
+    "######################################\n",
+    "# These propagate a scalar input gradient (in the surface roughness) to the output pixels,\n",
+    "# i.e. produce a gradient image as a result.\n",
+    "\n",
+    "for method, method_name in zip(diff_methods, diff_method_names):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_forward/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "        \n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method})\n",
+    "        \n",
+    "    scene = mi.load_file('scenes/dielectric_plane_vmf.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key, kappa_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    \n",
+    "    for idx_kappa, kappa in enumerate(kappas):\n",
+    "        for idx_alpha, alpha in enumerate(alphas):\n",
+    "            print(\"  - {}/{}\".format(idx_kappa*len(alphas) + idx_alpha + 1, len(alphas)*len(kappas)), end='\\r') \n",
+    "            params[kappa_key] = kappa\n",
+    "            params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "            \n",
+    "            # Diff. input parameter\n",
+    "            pi = mi.Float(0.0)\n",
+    "            dr.enable_grad(pi)\n",
+    "            dr.set_grad(pi, 1.0)\n",
+    "            params[alpha_key] += pi\n",
+    "            params.update()\n",
+    "            \n",
+    "            image_grad = mi.Float(0.0)\n",
+    "            if 'diff' in method:\n",
+    "                # Differential sampling strategy, use antithetic sampling.\n",
+    "                # Note that we use the same seed twice, and use half the number of samples for each pass.\n",
+    "                image = mi.render(scene, params,\n",
+    "                                  integrator=integrator, seed=0, spp=spp//2, antithetic_pass=False)\n",
+    "                dr.forward(pi)\n",
+    "                image_grad = dr.grad(image)\n",
+    "                \n",
+    "                params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "                params[alpha_key] += pi\n",
+    "                params.update()\n",
+    "                \n",
+    "                image = mi.render(scene, params,\n",
+    "                                  integrator=integrator, seed=0, spp=spp//2, antithetic_pass=True)\n",
+    "                dr.forward(pi)\n",
+    "                image_grad += dr.grad(image)\n",
+    "                \n",
+    "                # Average both passes\n",
+    "                image_grad *= 0.5\n",
+    "            else:\n",
+    "                # Produce differentiable rendering\n",
+    "                image = mi.render(scene, params,\n",
+    "                                  integrator=integrator, seed=0, spp=spp)\n",
+    "                # And propagate derivatives forwards through it\n",
+    "                dr.forward(pi)\n",
+    "                image_grad = dr.grad(image)\n",
+    "                \n",
+    "            # Save output gradient image\n",
+    "            outname = \"{}/k_{:02d}_a_{:02d}.exr\".format(method_dir, idx_kappa, idx_alpha)\n",
+    "            mi.util.convert_to_bitmap(image_grad, uint8_srgb=False).write(outname)\n",
+    "\n",
+    "            \n",
+    "produce_plots(diff_methods, diff_method_names, 'gradients_forward', gradients=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "71be7f79",
+   "metadata": {
+    "ExecuteTime": {
+     "end_time": "2022-08-14T04:41:42.409267Z",
+     "start_time": "2022-08-14T04:39:13.115223Z"
+    },
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "######################################\n",
+    "# Gradient estimators (REVERSE mode) #\n",
+    "######################################\n",
+    "# These propagate an output image gradient to the (textured) scene input parameters,\n",
+    "# i.e. the result are texture gradient images.\n",
+    "\n",
+    "for method, method_name in zip(diff_methods, diff_method_names):\n",
+    "    print(\"* {}\".format(method_name))\n",
+    "    method_dir = \"{}/gradients_reverse/{}\".format(base_dir, method)\n",
+    "    if not os.path.exists(method_dir):\n",
+    "        os.makedirs(method_dir)\n",
+    "        \n",
+    "    integrator = mi.load_dict({'type': 'estimator_comparison', 'method': method})\n",
+    "        \n",
+    "    scene = mi.load_file('scenes/dielectric_plane_vmf.xml')\n",
+    "    params = mi.traverse(scene)\n",
+    "    params.keep([alpha_key, kappa_key])\n",
+    "    texture_shape = params[alpha_key].shape\n",
+    "    img_shape     = scene.sensors()[0].film().crop_size()\n",
+    "    \n",
+    "    for idx_kappa, kappa in enumerate(kappas):\n",
+    "        for idx_alpha, alpha in enumerate(alphas):\n",
+    "            print(\"  - {}/{}\".format(idx_kappa*len(alphas) + idx_alpha + 1, len(alphas)*len(kappas)), end='\\r') \n",
+    "            params[kappa_key] = kappa\n",
+    "            params[alpha_key] = mi.TensorXf(np.ones(texture_shape)*alpha)\n",
+    "            params.update()\n",
+    "            \n",
+    "            dr.enable_grad(params[alpha_key])\n",
+    "            \n",
+    "            grad_backward = mi.Float(0.0)\n",
+    "            if 'diff' in method:\n",
+    "                # Differential sampling strategy, use antithetic sampling.\n",
+    "                # Note that we use the same seed twice, and use half the number of samples for each pass.\n",
+    "                image = mi.render(scene, params,\n",
+    "                                  integrator=integrator, seed=0, spp=spp//2, antithetic_pass=False)\n",
+    "                dr.backward(image)\n",
+    "                \n",
+    "                image = mi.render(scene, params,\n",
+    "                                  integrator=integrator, seed=0, spp=spp//2, antithetic_pass=True)\n",
+    "                dr.backward(image)\n",
+    "                \n",
+    "                # Average accumulated gradients from both passes\n",
+    "                grad_backward = 0.5 * dr.grad(params[alpha_key])\n",
+    "            else:\n",
+    "                # Produce differentiable rendering\n",
+    "                image = mi.render(scene, params,\n",
+    "                                  integrator=integrator, seed=0, spp=spp)\n",
+    "                \n",
+    "                # And propagate derivatives backwards through it\n",
+    "                dr.backward(image)\n",
+    "                grad_backward = dr.grad(params[alpha_key])\n",
+    "                \n",
+    "            # Save output gradient image\n",
+    "            outname = \"{}/k_{:02d}_a_{:02d}.exr\".format(method_dir, idx_kappa, idx_alpha)\n",
+    "            mi.util.convert_to_bitmap(grad_backward, uint8_srgb=False).write(outname)\n",
+    "    \n",
+    "produce_plots(diff_methods, diff_method_names, 'gradients_reverse', gradients=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a93cc8d2",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/results/scenes/conductor_plane_envmap.xml b/results/scenes/conductor_plane_envmap.xml
new file mode 100644
index 00000000..a550a90f
--- /dev/null
+++ b/results/scenes/conductor_plane_envmap.xml
@@ -0,0 +1,56 @@
+<scene version="2.0.0">
+    <default name="spp" value="32"/>
+    <default name="res" value="256"/>
+
+    <integrator type="path"/>
+
+    <sensor type="perspective">
+        <string name="fov_axis" value="smaller"/>
+        <float name="near_clip" value="0.001"/>
+        <float name="far_clip" value="1000"/>
+        <float name="focus_distance" value="1000"/>
+        <float name="fov" value="60"/>
+        <transform name="to_world">
+            <lookat origin="0, 0, 2"
+                    target="0, 0, 0"
+                    up    ="0, 1, 0"/>
+            <rotate x="1" angle="-70"/>
+        </transform>
+        <sampler type="independent">
+            <integer name="sample_count" value="$spp"/>
+        </sampler>
+        <film type="hdrfilm">
+            <integer name="width" value="$res"/>
+            <integer name="height" value="$res"/>
+            <rfilter type="gaussian"/>
+            <string name="pixel_format" value="rgb"/>
+        </film>
+    </sensor>
+
+    <shape type="obj" id="plane">
+        <string name="filename" value="meshes/rounded_plane.obj"/>
+
+        <boolean name="face_normals" value="true"/>
+
+        <bsdf type="roughconductor">
+            <texture name="alpha" type="bitmap">
+                <string name="filename" value="textures/constant_roughness.jpg"/>
+                <boolean name="raw" value="true"/>
+            </texture>
+
+            <string name="material" value="Au"/>
+            <boolean name="differential_sampling" value="true"/>
+            <boolean name="principled_roughness_mapping" value="false"/>
+            <boolean name="sample_visible" value="true"/>
+        </bsdf>
+
+        <transform name="to_world">
+            <scale x="1" y="1" z="1"/>
+            <translate z="-0.2"/>
+        </transform>
+    </shape>
+
+    <emitter type="envmap" id="emitter">
+        <string name="filename" value="textures/envmap.exr"/>
+    </emitter>
+</scene>
diff --git a/results/scenes/microfacet_plane.xml b/results/scenes/conductor_plane_vmf.xml
similarity index 100%
rename from results/scenes/microfacet_plane.xml
rename to results/scenes/conductor_plane_vmf.xml
diff --git a/results/scenes/dielectric_plane_envmap.xml b/results/scenes/dielectric_plane_envmap.xml
new file mode 100644
index 00000000..73fb7d7a
--- /dev/null
+++ b/results/scenes/dielectric_plane_envmap.xml
@@ -0,0 +1,55 @@
+<scene version="2.0.0">
+    <default name="spp" value="32"/>
+    <default name="res" value="256"/>
+
+    <integrator type="path"/>
+
+    <sensor type="perspective">
+        <string name="fov_axis" value="smaller"/>
+        <float name="near_clip" value="0.001"/>
+        <float name="far_clip" value="1000"/>
+        <float name="focus_distance" value="1000"/>
+        <float name="fov" value="60"/>
+        <transform name="to_world">
+            <lookat origin="0, 0, 2"
+                    target="0, 0, 0"
+                    up    ="0, 1, 0"/>
+            <rotate x="1" angle="-70"/>
+        </transform>
+        <sampler type="independent">
+            <integer name="sample_count" value="$spp"/>
+        </sampler>
+        <film type="hdrfilm">
+            <integer name="width" value="$res"/>
+            <integer name="height" value="$res"/>
+            <rfilter type="gaussian"/>
+            <string name="pixel_format" value="rgb"/>
+        </film>
+    </sensor>
+
+    <shape type="obj" id="plane">
+        <string name="filename" value="meshes/rounded_plane.obj"/>
+
+        <boolean name="face_normals" value="true"/>
+
+        <bsdf type="roughdielectric">
+            <texture name="alpha" type="bitmap">
+                <string name="filename" value="textures/constant_roughness.jpg"/>
+                <boolean name="raw" value="true"/>
+            </texture>
+
+            <boolean name="differential_sampling" value="true"/>
+            <boolean name="principled_roughness_mapping" value="false"/>
+            <boolean name="sample_visible" value="true"/>
+        </bsdf>
+
+        <transform name="to_world">
+            <scale x="1" y="1" z="1"/>
+            <translate z="-0.2"/>
+        </transform>
+    </shape>
+
+    <emitter type="envmap" id="emitter">
+        <string name="filename" value="textures/envmap.exr"/>
+    </emitter>
+</scene>
diff --git a/results/scenes/dielectric_plane_vmf.xml b/results/scenes/dielectric_plane_vmf.xml
new file mode 100644
index 00000000..6ea44c3f
--- /dev/null
+++ b/results/scenes/dielectric_plane_vmf.xml
@@ -0,0 +1,58 @@
+<scene version="2.0.0">
+    <default name="spp" value="32"/>
+    <default name="res" value="256"/>
+
+    <integrator type="path"/>
+
+    <sensor type="perspective">
+        <string name="fov_axis" value="smaller"/>
+        <float name="near_clip" value="0.001"/>
+        <float name="far_clip" value="1000"/>
+        <float name="focus_distance" value="1000"/>
+        <float name="fov" value="60"/>
+        <transform name="to_world">
+            <lookat origin="0, 0, 2"
+                    target="0, 0, 0"
+                    up    ="0, 1, 0"/>
+            <rotate x="1" angle="-70"/>
+        </transform>
+        <sampler type="independent">
+            <integer name="sample_count" value="$spp"/>
+        </sampler>
+        <film type="hdrfilm">
+            <integer name="width" value="$res"/>
+            <integer name="height" value="$res"/>
+            <rfilter type="gaussian"/>
+            <string name="pixel_format" value="rgb"/>
+        </film>
+    </sensor>
+
+    <shape type="obj" id="plane">
+        <string name="filename" value="meshes/rounded_plane.obj"/>
+
+        <boolean name="face_normals" value="true"/>
+
+        <bsdf type="roughdielectric">
+            <texture name="alpha" type="bitmap">
+                <string name="filename" value="textures/constant_roughness.jpg"/>
+                <boolean name="raw" value="true"/>
+            </texture>
+
+            <boolean name="differential_sampling" value="true"/>
+            <boolean name="principled_roughness_mapping" value="false"/>
+            <boolean name="sample_visible" value="true"/>
+        </bsdf>
+
+        <transform name="to_world">
+            <scale x="1" y="1" z="1"/>
+            <translate z="-0.2"/>
+        </transform>
+    </shape>
+
+    <emitter type="vmf" id="emitter">
+        <float name="kappa" value="100"/>
+        <transform name="to_world">
+            <rotate x="1" angle="-120"/>
+        </transform>
+    </emitter>
+</scene>