diff --git a/src/libbsdl/CMakeLists.txt b/src/libbsdl/CMakeLists.txt
index 5e551b2bf..aa1d68e8e 100644
--- a/src/libbsdl/CMakeLists.txt
+++ b/src/libbsdl/CMakeLists.txt
@@ -7,4 +7,13 @@ include(bsdl.cmake)
 # BSDL has small tables for doing spectral render in sRGB
 # but for other color spaces tell this function which ones
 # and it will bake Jakob-Hanika coefficient tables.
-add_bsdl_library(BSDL) # SPECTRAL_COLOR_SPACES "ACEScg")
\ No newline at end of file
+add_bsdl_library(BSDL) # SPECTRAL_COLOR_SPACES "ACEScg")
+
+
+if (OSL_BUILD_TESTS AND BUILD_TESTING)
+    add_executable (bsdl_test bsdl_test.cpp)
+    set_target_properties (bsdl_test PROPERTIES FOLDER "Unit Tests")
+    target_include_directories (bsdl_test BEFORE PRIVATE ${OpenImageIO_INCLUDES})
+    target_link_libraries (bsdl_test PRIVATE BSDL OpenImageIO::OpenImageIO)
+    add_test (unit_bsdl ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/bsdl_test)
+endif()
diff --git a/src/libbsdl/bsdl_test.cpp b/src/libbsdl/bsdl_test.cpp
new file mode 100644
index 000000000..9d0f54962
--- /dev/null
+++ b/src/libbsdl/bsdl_test.cpp
@@ -0,0 +1,93 @@
+// Copyright Contributors to the Open Shading Language project.
+// SPDX-License-Identifier: BSD-3-Clause
+// https://github.com/AcademySoftwareFoundation/OpenShadingLanguage
+
+#define BSDL_UNROLL()
+
+#include <BSDL/config.h>
+
+using BSDLConfig = bsdl::BSDLDefaultConfig;
+
+#include <BSDL/MTX/bsdf_schlick_impl.h>
+
+#include <OpenImageIO/unittest.h>
+
+using namespace bsdl;
+
+
+
+struct TestDielectricBSDF : mtx::DielectricBSDF<mtx::SchlickFresnel> {
+    using Base = mtx::DielectricBSDF<mtx::SchlickFresnel>;
+    using Base::Base;
+    using Base::reflection_probability;
+};
+
+
+
+static void
+test_colored_schlick_sampling()
+{
+    const float rgb[] = { 0.1f, 0.2f, 0.8f };
+    const Power F([&](int i) { return rgb[i]; }, 0.0f);
+    const mtx::SchlickFresnel fresnel(F, F, 5.0f, 1.5f, false);
+    const GGXDist dist(0.25f, 0.0f);
+    const Imath::V3f wo(0.0f, 0.0f, 1.0f);
+    const TestDielectricBSDF bsdf(dist, fresnel, wo.z, 0.25f, true, 0.0f);
+
+    const float probability = (rgb[0] + rgb[1] + rgb[2]) / 3.0f;
+    OIIO_CHECK_EQUAL_THRESH(bsdf.reflection_probability(F), probability, 1e-6f);
+    OIIO_CHECK_EQUAL_THRESH(bsdf.reflection_probability(Power::UNIT()), 1.0f,
+                            1e-6f);
+
+    const Sample reflected   = bsdf.sample(wo, 0.5f, 0.5f, probability - 0.01f);
+    const Sample transmitted = bsdf.sample(wo, 0.5f, 0.5f, probability + 0.01f);
+    OIIO_CHECK_ASSERT(reflected.wi.z > 0.0f);
+    OIIO_CHECK_ASSERT(transmitted.wi.z < 0.0f);
+
+    const float rgb2[] = { 0.2f, 0.4f, 0.6f };
+    const Power F2([&](int i) { return rgb2[i]; }, 0.0f);
+    const mtx::SchlickFresnel fresnel2(F2, F2, 5.0f, 1.5f, false);
+    const TestDielectricBSDF bsdf2(dist, fresnel2, wo.z, 0.25f, true, 0.0f);
+    const float probability2 = (rgb2[0] + rgb2[1] + rgb2[2]) / 3.0f;
+
+    const Sample reflected2   = bsdf2.sample(wo, 0.5f, 0.5f, 0.3f);
+    const Sample transmitted2 = bsdf2.sample(wo, 0.5f, 0.5f, 0.8f);
+    const Sample reflected1   = bsdf.sample(wo, 0.5f, 0.5f, 0.3f);
+    const Sample transmitted1 = bsdf.sample(wo, 0.5f, 0.5f, 0.8f);
+    OIIO_CHECK_EQUAL_THRESH(reflected1.pdf / reflected2.pdf,
+                            probability / probability2, 1e-6f);
+    const float transmission_probability  = 1.0f - probability;
+    const float transmission_probability2 = 1.0f - probability2;
+    OIIO_CHECK_EQUAL_THRESH(transmitted1.pdf / transmitted2.pdf,
+                            transmission_probability
+                                / transmission_probability2,
+                            1e-6f);
+    for (int i = 0; i < 3; ++i) {
+        OIIO_CHECK_EQUAL_THRESH(reflected1.weight[i] / reflected2.weight[i],
+                                (rgb[i] * probability2)
+                                    / (rgb2[i] * probability),
+                                1e-6f);
+        OIIO_CHECK_EQUAL_THRESH(transmitted1.weight[i] / transmitted2.weight[i],
+                                ((1.0f - rgb[i]) * transmission_probability2)
+                                    / ((1.0f - rgb2[i])
+                                       * transmission_probability),
+                                1e-6f);
+    }
+
+    const float spectral[] = { 0.2f, 0.4f, 0.6f, 0.8f };
+    const Power Fs([&](int i) { return spectral[i]; }, 500.0f);
+    const mtx::SchlickFresnel spectral_fresnel(Fs, Fs, 5.0f, 1.5f, false);
+    const TestDielectricBSDF spectral_bsdf(dist, spectral_fresnel, wo.z, 0.25f,
+                                           true, 500.0f);
+    OIIO_CHECK_EQUAL_THRESH(spectral_bsdf.reflection_probability(Fs), 0.5f,
+                            1e-6f);
+}
+
+
+
+int
+main(int /*argc*/, char* /*argv*/[])
+{
+    test_colored_schlick_sampling();
+    return unit_test_failures;
+}
diff --git a/src/libbsdl/include/BSDL/MTX/bsdf_dielectric_decl.h b/src/libbsdl/include/BSDL/MTX/bsdf_dielectric_decl.h
index 19a052e57..37da18b55 100644
--- a/src/libbsdl/include/BSDL/MTX/bsdf_dielectric_decl.h
+++ b/src/libbsdl/include/BSDL/MTX/bsdf_dielectric_decl.h
@@ -42,7 +42,7 @@ template<typename Fresnel> struct DielectricBSDF {
 
     BSDL_INLINE_METHOD
     DielectricBSDF(const GGXDist& dist, const Fresnel& fresnel, float cosNO,
-                   float roughness, bool dorefr);
+                   float roughness, bool dorefr, float lambda_0);
 
     DielectricBSDF() = default;
 
@@ -58,9 +58,12 @@ template<typename Fresnel> struct DielectricBSDF {
     static constexpr const char* NS = "mtx";
 
 protected:
+    BSDL_INLINE_METHOD float reflection_probability(const Power& F) const;
+
     GGXDist d;
     Fresnel f;
     float E_ms;
+    float lambda_0;
     bool dorefr;
 };
 
diff --git a/src/libbsdl/include/BSDL/MTX/bsdf_dielectric_impl.h b/src/libbsdl/include/BSDL/MTX/bsdf_dielectric_impl.h
index c79c65fb5..6f8e183e1 100644
--- a/src/libbsdl/include/BSDL/MTX/bsdf_dielectric_impl.h
+++ b/src/libbsdl/include/BSDL/MTX/bsdf_dielectric_impl.h
@@ -67,8 +67,9 @@ template<typename Fresnel>
 BSDL_INLINE_METHOD
 DielectricBSDF<Fresnel>::DielectricBSDF(const GGXDist& dist,
                                         const Fresnel& fresnel, float cosNO,
-                                        float roughness, bool dorefr)
-    : d(dist), f(fresnel), dorefr(dorefr)
+                                        float roughness, bool dorefr,
+                                        float lambda_0)
+    : d(dist), f(fresnel), lambda_0(lambda_0), dorefr(dorefr)
 {
     if (!dorefr) {
         TabulatedEnergyCurve<spi::MiniMicrofacetGGX> curve(roughness, 0.0f);
@@ -83,7 +84,7 @@ DielectricReflFront::DielectricReflFront(float cosNO, float roughness_index,
     : DielectricBSDF<DielectricFresnel>(
         GGXDist(roughness_index, 0),
         DielectricFresnel::from_table_index(fresnel_index, false), cosNO,
-        roughness_index, false)
+        roughness_index, false, 1)
 {
 }
 
@@ -93,7 +94,7 @@ DielectricBothFront::DielectricBothFront(float cosNO, float roughness_index,
     : DielectricBSDF<DielectricFresnel>(
         GGXDist(roughness_index, 0),
         DielectricFresnel::from_table_index(fresnel_index, false), cosNO,
-        roughness_index, true)
+        roughness_index, true, 1)
 {
 }
 
@@ -103,10 +104,19 @@ DielectricBothBack::DielectricBothBack(float cosNO, float roughness_index,
     : DielectricBSDF<DielectricFresnel>(
         GGXDist(roughness_index, 0),
         DielectricFresnel::from_table_index(fresnel_index, true), cosNO,
-        roughness_index, true)
+        roughness_index, true, 1)
 {
 }
 
+template<typename Fresnel>
+BSDL_INLINE_METHOD float
+DielectricBSDF<Fresnel>::reflection_probability(const Power& F) const
+{
+    // Fresnel may return Power::UNIT(), so mask the inactive RGB lane before
+    // averaging.
+    return F.cliped_rgb(lambda_0).avg(lambda_0);
+}
+
 template<typename Fresnel>
 BSDL_INLINE_METHOD Sample
 DielectricBSDF<Fresnel>::eval(Imath::V3f wo, Imath::V3f wi) const
@@ -125,18 +135,19 @@ DielectricBSDF<Fresnel>::eval(Imath::V3f wo, Imath::V3f wi) const
         const float D  = d.D(m);
         const float G1 = d.G1(wo);
         const Power F  = f.eval(cosMO);
-        if (F.max() <= 0)
+        const float P  = reflection_probability(F);
+        if (P <= 0)
             return {};
         if constexpr (BSDLConfig::use_bvn_refraction) {
             // Reflection optimized density
             const float D_refl_D = d.D_refl_D(wo, m);
             const float D_refl   = D_refl_D * D;
-            const Power out = F * (d.G2_G1(wi, wo) * G1 / (D_refl_D * F.max()));
-            const float pdf = D_refl / (4.0f * cosNO) * F.max();
+            const Power out      = F * (d.G2_G1(wi, wo) * G1 / (D_refl_D * P));
+            const float pdf      = D_refl / (4.0f * cosNO) * P;
             return { wi, out, pdf, 0 };
         } else {
-            const Power out = F * d.G2_G1(wi, wo);
-            const float pdf = (G1 * D * F.max()) / (4.0f * cosNO);
+            const Power out = F * (d.G2_G1(wi, wo) / P);
+            const float pdf = (G1 * D * P) / (4.0f * cosNO);
             return { wi, out, pdf, 0 };
         }
     } else if (cosNI < 0) {
@@ -148,8 +159,10 @@ DielectricBSDF<Fresnel>::eval(Imath::V3f wo, Imath::V3f wi) const
         const float cosHI = Ht.dot(wi);
         if (cosHO <= 0 || cosHI >= 0)
             return {};
-        const Power Ft = Power::UNIT() - f.eval(cosHO);
-        if (Ht.z <= 0 || Ft.max() <= 0)
+        const Power F  = f.eval(cosHO);
+        const Power Ft = Power::UNIT() - F;
+        const float Pt = 1 - reflection_probability(F);
+        if (Ht.z <= 0 || Pt <= 0)
             return {};
         const float D  = d.D(Ht);
         const float G1 = d.G1(wo);
@@ -159,16 +172,15 @@ DielectricBSDF<Fresnel>::eval(Imath::V3f wo, Imath::V3f wi) const
             // Reflection optimized density
             const float D_refl_D = d.D_refl_D(wo, Ht);
             const float D_refl   = D_refl_D * D;
-            float pdf            = D_refl * J * Ft.max();
+            float pdf            = D_refl * J * Pt;
             const Power out      = Ft
                               * (d.G2_G1({ wi.x, wi.y, -wi.z }, wo) * G1
-                                 / (D_refl_D * Ft.max()));
+                                 / (D_refl_D * Pt));
             return { wi, out, pdf, 0 };
         } else {
-            const Power out = Ft
-                              * (d.G2_G1({ wi.x, wi.y, -wi.z }, wo) / Ft.max());
+            const Power out = Ft * (d.G2_G1({ wi.x, wi.y, -wi.z }, wo) / Pt);
 
-            float pdf = J * G1 * D * Ft.max();
+            float pdf = J * G1 * D * Pt;
             return { wi, out, pdf, 0 };
         }
 
@@ -196,8 +208,8 @@ DielectricBSDF<Fresnel>::sample(Imath::V3f wo, float randu, float randv,
     const float cosMO = wo.dot(m);
     if (cosMO <= 0)
         return {};
-    const float F       = f.eval(cosMO).max();
-    bool choose_reflect = randw < F;
+    const float P       = reflection_probability(f.eval(cosMO));
+    bool choose_reflect = randw < P;
     const Imath::V3f wi = choose_reflect ? reflect(wo, m)
                                          : refract(wo, m, f.refraction_eta());
     if ((choose_reflect && wi.z <= 0) || (!choose_reflect && wi.z >= 0))
@@ -245,7 +257,8 @@ DielectricLobe<BSDF_ROOT>::DielectricLobe(T* lobe, const BsdfGlobals& globals,
     DielectricFresnel fresnel(globals.relative_eta(IOR), globals.backfacing);
     E_ms = 0;
     spec = DielectricBSDF<DielectricFresnel>(GGXDist(roughness, aniso, rx < ry),
-                                             fresnel, cosNO, roughness, dorefr);
+                                             fresnel, cosNO, roughness, dorefr,
+                                             globals.lambda_0);
     if (dorefl && !dorefr) {
         E_ms = TabulatedEnergyCurve<DielectricReflFront>(roughness,
                                                          fresnel.table_index())
diff --git a/src/libbsdl/include/BSDL/MTX/bsdf_schlick_impl.h b/src/libbsdl/include/BSDL/MTX/bsdf_schlick_impl.h
index e58ec46b7..81570b508 100644
--- a/src/libbsdl/include/BSDL/MTX/bsdf_schlick_impl.h
+++ b/src/libbsdl/include/BSDL/MTX/bsdf_schlick_impl.h
@@ -73,7 +73,8 @@ SchlickLobe<BSDF_ROOT>::SchlickLobe(T* lobe, const BsdfGlobals& globals,
 
     E_ms = 0;
     spec = DielectricBSDF<SchlickFresnel>(GGXDist(roughness, aniso, rx < ry),
-                                          fresnel, cosNO, roughness, dorefr);
+                                          fresnel, cosNO, roughness, dorefr,
+                                          globals.lambda_0);
     if (dorefl && !dorefr) {
         // Energy compensation reuses the dielectric Fresnel albedo tables,
         // which assumes the Schlick curve matches the true dielectric Fresnel.