diff --git a/src/core/integrators/photon_map/PhotonTracer.cpp b/src/core/integrators/photon_map/PhotonTracer.cpp
@@ -128,7 +128,7 @@ static bool evalBeam1D(const PhotonBeam &beam, PathSampleGenerator &sampler, con
         mediumQuery.setFarT(t);
         beamEstimate += medium->sigmaT(hitPoint)*invSinTheta/(2.0f*radius)
                 *medium->phaseFunction(hitPoint)->eval(beam.dir, -ray.dir())
-                *medium->transmittance(sampler, mediumQuery)*beam.power;
+                *medium->transmittance(sampler, mediumQuery, true, false)*beam.power;
 
         return true;
     }
@@ -149,7 +149,7 @@ static bool evalPlane0D(const PhotonPlane0D &p, PathSampleGenerator &sampler, co
             mediumQuery.setFarT(t);
             beamEstimate += sqr(medium->sigmaT(hitPoint))*std::abs(invDet)
                     *medium->phaseFunction(hitPoint)->eval(p.d1, -ray.dir())
-                    *medium->transmittance(sampler, mediumQuery)*p.power;
+                    *medium->transmittance(sampler, mediumQuery, true, false)*p.power;
 
             return true;
         }
@@ -187,7 +187,7 @@ bool evalPlane1D(const PhotonPlane1D &p, PathSampleGenerator &sampler, const Ray
             Ray mediumQuery(ray);
             mediumQuery.setFarT(t);
 
-            controlVariate -= medium->transmittance(sampler, mediumQuery)*(maxT - minT);
+            controlVariate -= medium->transmittance(sampler, mediumQuery, true, false)*(maxT - minT);
         }
 
         beamEstimate += sqr(medium->sigmaT(v2))*medium->phaseFunction(v2)->eval(p.d1, -ray.dir())*p.power*controlVariate;
@@ -289,7 +289,7 @@ Vec3f PhotonTracer::traceSensorPath(Vec2u pixel, const KdTree<Photon> &surfaceTr
                     mediumQuery.setFarT(t);
                     estimate += (3.0f*INV_PI*sqr(1.0f - distSq/p.radiusSq))/p.radiusSq
                             *medium->phaseFunction(p.pos)->eval(p.dir, -ray.dir())
-                            *medium->transmittance(sampler, mediumQuery)*p.power;
+                            *medium->transmittance(sampler, mediumQuery, true, false)*p.power;
                 };
                 auto beamContribution = [&](uint32 photonIndex, const Vec3pf *bounds, float tMin, float tMax) {
                     const PhotonBeam &beam = beams[photonIndex];
@@ -343,7 +343,7 @@ Vec3f PhotonTracer::traceSensorPath(Vec2u pixel, const KdTree<Photon> &surfaceTr
 
                 result += throughput*estimate;
             }
-            throughput *= medium->transmittance(sampler, ray);
+            throughput *= medium->transmittance(sampler, ray, true, true);
         }
         if (!didHit || !includeSurfaces)
             break;

diff --git a/src/core/io/Scene.cpp b/src/core/io/Scene.cpp
@@ -91,6 +91,11 @@ std::shared_ptr<T> fetchObject(const std::vector<std::shared_ptr<T>> &list, cons
     }
 }
 
+std::shared_ptr<Transmittance> Scene::fetchTransmittance(JsonPtr value) const
+{
+    return instantiate<Transmittance>(value, *this);
+}
+
 std::shared_ptr<PhaseFunction> Scene::fetchPhase(JsonPtr value) const
 {
     return instantiate<PhaseFunction>(value, *this);

diff --git a/src/core/io/Scene.hpp b/src/core/io/Scene.hpp
@@ -13,6 +13,8 @@
 #include "ImageIO.hpp"
 #include "Path.hpp"
 
+#include "transmittances/Transmittance.hpp"
+
 #include "phasefunctions/PhaseFunction.hpp"
 
 #include "integrators/Integrator.hpp"
@@ -66,6 +68,7 @@ class Scene : public JsonSerializable
     virtual void loadResources() override;
     virtual void saveResources() override;
 
+    std::shared_ptr<Transmittance> fetchTransmittance(JsonPtr value) const;
     std::shared_ptr<PhaseFunction> fetchPhase(JsonPtr value) const;
     std::shared_ptr<Medium> fetchMedium(JsonPtr value) const;
     std::shared_ptr<Grid> fetchGrid(JsonPtr value) const;

diff --git a/src/core/media/AtmosphericMedium.cpp b/src/core/media/AtmosphericMedium.cpp
@@ -110,10 +110,10 @@ inline float AtmosphericMedium::densityIntegral(float h, float t0, float t1) con
         return (SQRT_PI*0.5f/s)*std::exp((-h*h + _radius*_radius)*s*s)*Erf::erfDifference(s*t0, s*t1);
 }
 
-inline float AtmosphericMedium::inverseOpticalDepth(double h, double t0, double sigmaT, double xi) const
+inline float AtmosphericMedium::inverseOpticalDepth(double h, double t0, double tau) const
 {
     double s = _effectiveFalloffScale;
-    double inner = std::erf(s*t0) - 2.0*double(INV_SQRT_PI)*std::exp(s*s*(h - _radius)*(h + _radius))*s*std::log(xi)/sigmaT;
+    double inner = std::erf(s*t0) + 2.0*double(INV_SQRT_PI)*std::exp(s*s*(h - _radius)*(h + _radius))*s*tau;
 
     if (inner >= 1.0)
         return Ray::infinity();
@@ -134,24 +134,26 @@ bool AtmosphericMedium::sampleDistance(PathSampleGenerator &sampler, const Ray &
     float maxT = ray.farT() + t0;
     if (_absorptionOnly) {
         sample.t = ray.farT();
-        sample.weight = std::exp(-_sigmaT*densityIntegral(h, t0, maxT));
+        Vec3f tau = densityIntegral(h, t0, maxT)*_sigmaT;
+        sample.weight = _transmittance->eval(tau, state.firstScatter, true);
         sample.pdf = 1.0f;
         sample.exited = true;
     } else {
         int component = sampler.nextDiscrete(3);
         float sigmaTc = _sigmaT[component];
-        float xi = 1.0f - sampler.next1D();
+        float tauC = _transmittance->sample(sampler, state.firstScatter)/sigmaTc;
 
-        float t = inverseOpticalDepth(h, t0, sigmaTc, xi);
+        float t = inverseOpticalDepth(h, t0, tauC);
         sample.t = min(t, maxT);
-        sample.weight = std::exp(-_sigmaT*densityIntegral(h, t0, sample.t));
+        Vec3f tau = densityIntegral(h, t0, sample.t)*_sigmaT;
+        sample.weight = _transmittance->eval(tau, state.firstScatter, sample.exited);
         sample.exited = (t >= maxT);
         if (sample.exited) {
-            sample.pdf = sample.weight.avg();
+            sample.pdf = _transmittance->surfaceProbability(tau, state.firstScatter).avg();
         } else {
             float rho = density(h, sample.t);
-            sample.pdf = (rho*_sigmaT*sample.weight).avg();
-            sample.weight *= rho*_sigmaS;
+            sample.pdf = (rho*_sigmaT*_transmittance->mediumPdf(tau, state.firstScatter)).avg();
+            sample.weight *= rho*_sigmaS*_transmittance->sigmaBar();
         }
         sample.weight /= sample.pdf;
         sample.t -= t0;
@@ -164,17 +166,19 @@ bool AtmosphericMedium::sampleDistance(PathSampleGenerator &sampler, const Ray &
     return true;
 }
 
-Vec3f AtmosphericMedium::transmittance(PathSampleGenerator &/*sampler*/, const Ray &ray) const
+Vec3f AtmosphericMedium::transmittance(PathSampleGenerator &/*sampler*/, const Ray &ray, bool startOnSurface,
+        bool endOnSurface) const
 {
     Vec3f p = (ray.pos() - _center);
     float t0 = p.dot(ray.dir());
     float t1 = ray.farT() + t0;
     float  h = (p - t0*ray.dir()).length();
 
-    return std::exp(-_sigmaT*densityIntegral(h, t0, t1));
+    Vec3f tau = _sigmaT*densityIntegral(h, t0, t1);
+    return _transmittance->eval(tau, startOnSurface, endOnSurface);
 }
 
-float AtmosphericMedium::pdf(PathSampleGenerator &/*sampler*/, const Ray &ray, bool onSurface) const
+float AtmosphericMedium::pdf(PathSampleGenerator &/*sampler*/, const Ray &ray, bool startOnSurface, bool endOnSurface) const
 {
     if (_absorptionOnly) {
         return 1.0f;
@@ -184,33 +188,13 @@ float AtmosphericMedium::pdf(PathSampleGenerator &/*sampler*/, const Ray &ray, b
         float t1 = ray.farT() + t0;
         float  h = (p - t0*ray.dir()).length();
 
-        Vec3f transmittance = std::exp(-_sigmaT*densityIntegral(h, t0, t1));
-        if (onSurface) {
-            return transmittance.avg();
-        } else {
-            return (density(h, t0)*_sigmaT*transmittance).avg();
-        }
-    }
-}
 
-Vec3f AtmosphericMedium::transmittanceAndPdfs(PathSampleGenerator &/*sampler*/, const Ray &ray, bool startOnSurface,
-        bool endOnSurface, float &pdfForward, float &pdfBackward) const
-{
-    Vec3f p = (ray.pos() - _center);
-    float t0 = p.dot(ray.dir());
-    float t1 = ray.farT() + t0;
-    float  h = (p - t0*ray.dir()).length();
-
-    Vec3f transmittance = std::exp(-_sigmaT*densityIntegral(h, t0, t1));
-
-    if (_absorptionOnly) {
-        pdfForward = pdfBackward = 1.0f;
-    } else {
-        pdfForward  =   endOnSurface ? transmittance.avg() : (density(h, t1)*_sigmaT*transmittance).avg();
-        pdfBackward = startOnSurface ? transmittance.avg() : (density(h, t0)*_sigmaT*transmittance).avg();
+        Vec3f tau = _sigmaT*densityIntegral(h, t0, t1);
+        if (endOnSurface)
+            return _transmittance->surfaceProbability(tau, startOnSurface).avg();
+        else
+            return (density(h, t0)*_sigmaT*_transmittance->mediumPdf(tau, startOnSurface)).avg();
     }
-
-    return transmittance;
 }
 
 }
diff --git a/src/core/media/AtmosphericMedium.hpp b/src/core/media/AtmosphericMedium.hpp
@@ -28,7 +28,7 @@ class AtmosphericMedium : public Medium
     inline float density(Vec3f p) const;
     inline float density(float h, float t0) const;
     inline float densityIntegral(float h, float t0, float t1) const;
-    inline float inverseOpticalDepth(double h, double t0, double sigmaT, double xi) const;
+    inline float inverseOpticalDepth(double h, double t0, double tau) const;
 
 public:
     AtmosphericMedium();
@@ -46,10 +46,9 @@ class AtmosphericMedium : public Medium
 
     virtual bool sampleDistance(PathSampleGenerator &sampler, const Ray &ray,
             MediumState &state, MediumSample &sample) const override;
-    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray) const override;
-    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool onSurface) const override;
-    virtual Vec3f transmittanceAndPdfs(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
-            bool endOnSurface, float &pdfForward, float &pdfBackward) const override;
+    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
+            bool endOnSurface) const override;
+    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface, bool endOnSurface) const override;
 };
 
 }

diff --git a/src/core/media/ExponentialMedium.cpp b/src/core/media/ExponentialMedium.cpp
@@ -92,14 +92,13 @@ inline float ExponentialMedium::densityIntegral(float x, float dx, float tMax) c
         return (std::exp(-x) - std::exp(-dx*tMax - x))/dx;
 }
 
-inline float ExponentialMedium::inverseOpticalDepth(float x, float dx, float sigmaT, float logXi) const
+inline float ExponentialMedium::inverseOpticalDepth(float x, float dx, float tau) const
 {
     if (dx == 0.0f) {
-        float effectiveSigmaTc = sigmaT*std::exp(-x);
-        return -logXi/effectiveSigmaTc;
+        return tau/std::exp(-x);
     } else {
-        float denom = sigmaT + dx*std::exp(x)*logXi;
-        return denom <= 0.0f ? Ray::infinity() : std::log(sigmaT/denom)/dx;
+        float denom = 1.0f - dx*std::exp(x)*tau;
+        return denom <= 0.0f ? Ray::infinity() : -std::log(denom)/dx;
     }
 }
 
@@ -117,25 +116,26 @@ bool ExponentialMedium::sampleDistance(PathSampleGenerator &sampler, const Ray &
         if (maxT == Ray::infinity() && dx <= 0.0f)
             return false;
         sample.t = maxT;
-        sample.weight = std::exp(-_sigmaT*densityIntegral(x, dx, ray.farT()));
+        Vec3f tau = densityIntegral(x, dx, ray.farT())*_sigmaT;
+        sample.weight = _transmittance->eval(tau, state.firstScatter, true);
         sample.pdf = 1.0f;
         sample.exited = true;
     } else {
         int component = sampler.nextDiscrete(3);
         float sigmaTc = _sigmaT[component];
-        float xi = 1.0f - sampler.next1D();
-        float logXi = std::log(xi);
+        float tauC = _transmittance->sample(sampler, state.firstScatter)/sigmaTc;
 
-        float t = inverseOpticalDepth(x, dx, sigmaTc, logXi);
+        float t = inverseOpticalDepth(x, dx, tauC);
         sample.t = min(t, maxT);
-        sample.weight = std::exp(-_sigmaT*densityIntegral(x, dx, sample.t));
+        Vec3f tau = densityIntegral(x, dx, sample.t)*_sigmaT;
+        sample.weight = _transmittance->eval(tau, state.firstScatter, sample.exited);
         sample.exited = (t >= maxT);
         if (sample.exited) {
-            sample.pdf = sample.weight.avg();
+            sample.pdf = _transmittance->surfaceProbability(tau, state.firstScatter).avg();
         } else {
             float rho = density(x, dx, sample.t);
-            sample.pdf = (rho*_sigmaT*sample.weight).avg();
-            sample.weight *= rho*_sigmaS;
+            sample.pdf = (rho*_sigmaT*_transmittance->mediumPdf(tau, state.firstScatter)).avg();
+            sample.weight *= rho*_sigmaS*_transmittance->sigmaBar();
         }
         sample.weight /= sample.pdf;
 
@@ -146,55 +146,34 @@ bool ExponentialMedium::sampleDistance(PathSampleGenerator &sampler, const Ray &
 
     return true;
 }
-Vec3f ExponentialMedium::transmittance(PathSampleGenerator &/*sampler*/, const Ray &ray) const
+Vec3f ExponentialMedium::transmittance(PathSampleGenerator &/*sampler*/, const Ray &ray, bool startOnSurface,
+        bool endOnSurface) const
 {
     float  x = _falloffScale*(ray.pos() - _unitPoint).dot(_unitFalloffDirection);
     float dx = _falloffScale*ray.dir().dot(_unitFalloffDirection);
 
-    if (ray.farT() == Ray::infinity() && dx <= 0.0f)
+    if (ray.farT() == Ray::infinity() && dx <= 0.0f) {
         return Vec3f(0.0f);
-    else
-        return std::exp(-_sigmaT*densityIntegral(x, dx, ray.farT()));
+    } else {
+        Vec3f tau = _sigmaT*densityIntegral(x, dx, ray.farT());
+        return _transmittance->eval(tau, startOnSurface, endOnSurface);
+    }
 }
 
-float ExponentialMedium::pdf(PathSampleGenerator &/*sampler*/, const Ray &ray, bool onSurface) const
+float ExponentialMedium::pdf(PathSampleGenerator &/*sampler*/, const Ray &ray, bool startOnSurface, bool endOnSurface) const
 {
     if (_absorptionOnly) {
         return 1.0f;
     } else {
         float  x = _falloffScale*(ray.pos() - _unitPoint).dot(_unitFalloffDirection);
         float dx = _falloffScale*ray.dir().dot(_unitFalloffDirection);
 
-        Vec3f transmittance = std::exp(-_sigmaT*densityIntegral(x, dx, ray.farT()));
-        if (onSurface) {
-            return transmittance.avg();
-        } else {
-            return (density(x, dx, ray.farT())*_sigmaT*transmittance).avg();
-        }
+        Vec3f tau = _sigmaT*densityIntegral(x, dx, ray.farT());
+        if (endOnSurface)
+            return _transmittance->surfaceProbability(tau, startOnSurface).avg();
+        else
+            return (density(x, dx, ray.farT())*_sigmaT*_transmittance->mediumPdf(tau, startOnSurface)).avg();
     }
 }
 
-Vec3f ExponentialMedium::transmittanceAndPdfs(PathSampleGenerator &/*sampler*/, const Ray &ray, bool startOnSurface,
-        bool endOnSurface, float &pdfForward, float &pdfBackward) const
-{
-    float  x = _falloffScale*(ray.pos() - _unitPoint).dot(_unitFalloffDirection);
-    float dx = _falloffScale*ray.dir().dot(_unitFalloffDirection);
-
-    if (ray.farT() == Ray::infinity() && dx <= 0.0f) {
-        pdfForward = pdfBackward = 0.0f;
-        return Vec3f(0.0f);
-    }
-
-    Vec3f transmittance = std::exp(-_sigmaT*densityIntegral(x, dx, ray.farT()));
-
-    if (_absorptionOnly) {
-        pdfForward = pdfBackward = 1.0f;
-    } else {
-        pdfForward  =   endOnSurface ? transmittance.avg() : (density(x, dx, ray.farT())*_sigmaT*transmittance).avg();
-        pdfBackward = startOnSurface ? transmittance.avg() : (density(x, dx,       0.0f)*_sigmaT*transmittance).avg();
-    }
-
-    return transmittance;
-}
-
 }
diff --git a/src/core/media/ExponentialMedium.hpp b/src/core/media/ExponentialMedium.hpp
@@ -21,7 +21,7 @@ class ExponentialMedium : public Medium
     inline float density(Vec3f p) const;
     inline float density(float x, float dx, float t) const;
     inline float densityIntegral(float x, float dx, float tMax) const;
-    inline float inverseOpticalDepth(float x, float dx, float sigmaT, float logXi) const;
+    inline float inverseOpticalDepth(float x, float dx, float tau) const;
 
 public:
     ExponentialMedium();
@@ -39,10 +39,9 @@ class ExponentialMedium : public Medium
 
     virtual bool sampleDistance(PathSampleGenerator &sampler, const Ray &ray,
             MediumState &state, MediumSample &sample) const override;
-    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray) const override;
-    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool onSurface) const override;
-    virtual Vec3f transmittanceAndPdfs(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
-            bool endOnSurface, float &pdfForward, float &pdfBackward) const override;
+    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
+            bool endOnSurface) const override;
+    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface, bool endOnSurface) const override;
 };
 
 }

diff --git a/src/core/media/HomogeneousMedium.cpp b/src/core/media/HomogeneousMedium.cpp
@@ -3,7 +3,6 @@
 #include "sampling/PathSampleGenerator.hpp"
 
 #include "math/TangentFrame.hpp"
-#include "math/FastMath.hpp"
 #include "math/Ray.hpp"
 
 #include "io/JsonObject.hpp"
@@ -74,27 +73,29 @@ bool HomogeneousMedium::sampleDistance(PathSampleGenerator &sampler, const Ray &
         if (maxT == Ray::infinity())
             return false;
         sample.t = maxT;
-        sample.weight = FastMath::exp(-_sigmaT*maxT);
+        sample.weight = _transmittance->eval(sample.t*_sigmaT, state.firstScatter, true);
         sample.pdf = 1.0f;
         sample.exited = true;
     } else {
         int component = sampler.nextDiscrete(3);
         float sigmaTc = _sigmaT[component];
 
-        float t = -std::log(1.0f - sampler.next1D())/sigmaTc;
+        float t = _transmittance->sample(sampler, state.firstScatter)/sigmaTc;
         sample.t = min(t, maxT);
         sample.continuedT = t;
-        sample.weight = FastMath::exp(-sample.t*_sigmaT);
-        sample.continuedWeight = FastMath::exp(-sample.continuedT*_sigmaT);
         sample.exited = (t >= maxT);
+        Vec3f tau = sample.t*_sigmaT;
+        Vec3f continuedTau = sample.continuedT*_sigmaT;
+        sample.weight = _transmittance->eval(tau, state.firstScatter, sample.exited);
+        sample.continuedWeight = _transmittance->eval(continuedTau, state.firstScatter, sample.exited);
         if (sample.exited) {
-            sample.pdf = sample.weight.avg();
+            sample.pdf = _transmittance->surfaceProbability(tau, state.firstScatter).avg();
         } else {
-            sample.pdf = (_sigmaT*sample.weight).avg();
-            sample.weight *= _sigmaS;
+            sample.pdf = (_sigmaT*_transmittance->mediumPdf(tau, state.firstScatter)).avg();
+            sample.weight *= _sigmaS*_transmittance->sigmaBar();
         }
         sample.weight /= sample.pdf;
-        sample.continuedWeight = _sigmaS*sample.continuedWeight/(_sigmaT*sample.continuedWeight).avg();
+        sample.continuedWeight = _sigmaS*_transmittance->sigmaBar()*sample.continuedWeight/(_sigmaT*_transmittance->mediumPdf(continuedTau, state.firstScatter)).avg();
 
         state.advance();
     }
@@ -104,87 +105,26 @@ bool HomogeneousMedium::sampleDistance(PathSampleGenerator &sampler, const Ray &
     return true;
 }
 
-bool HomogeneousMedium::invertDistance(WritablePathSampleGenerator &sampler, const Ray &ray, bool onSurface) const
-{
-    float maxT = ray.farT();
-    if (_absorptionOnly) {
-        if (maxT == Ray::infinity())
-            return false;
-        return true;
-    } else {
-        Vec3f Tr = std::exp(-_sigmaT*maxT);
-        Vec3f pdfs = onSurface ? Tr : _sigmaT*Tr;
-        Vec3f cdf(pdfs.x(), pdfs.x() + pdfs.y(), pdfs.sum());
-
-        float target = sampler.untracked1D()*cdf.z();
-        int component = target < cdf.x() ? 0 : (target < cdf.y() ? 1 : 2);
-        float Trc = Tr[component];
-
-        float xi = 1.0f - Trc;
-        if (onSurface)
-            xi -= (1.0f - Trc)*sampler.next1D();
-
-        sampler.putDiscrete(3, component);
-        sampler.put1D(xi);
-    }
-    return true;
-}
-
-Vec3f HomogeneousMedium::transmittance(PathSampleGenerator &/*sampler*/, const Ray &ray) const
+Vec3f HomogeneousMedium::transmittance(PathSampleGenerator &/*sampler*/, const Ray &ray, bool startOnSurface,
+        bool endOnSurface) const
 {
     if (ray.farT() == Ray::infinity())
         return Vec3f(0.0f);
-    else {
-        return FastMath::exp(-_sigmaT*ray.farT());
-    }
+    else
+        return _transmittance->eval(_sigmaT*ray.farT(), startOnSurface, endOnSurface);
 }
 
-float HomogeneousMedium::pdf(PathSampleGenerator &/*sampler*/, const Ray &ray, bool onSurface) const
+float HomogeneousMedium::pdf(PathSampleGenerator &/*sampler*/, const Ray &ray, bool startOnSurface, bool endOnSurface) const
 {
     if (_absorptionOnly) {
         return 1.0f;
     } else {
-        if (onSurface)
-            return FastMath::exp(-ray.farT()*_sigmaT).avg();
+        Vec3f tau = ray.farT()*_sigmaT;
+        if (endOnSurface)
+            return _transmittance->surfaceProbability(tau, startOnSurface).avg();
         else
-            return (_sigmaT*FastMath::exp(-ray.farT()*_sigmaT)).avg();
+            return (_sigmaT*_transmittance->mediumPdf(tau, startOnSurface)).avg();
     }
 }
 
-Vec3f HomogeneousMedium::transmittanceAndPdfs(PathSampleGenerator &/*sampler*/, const Ray &ray, bool startOnSurface,
-        bool endOnSurface, float &pdfForward, float &pdfBackward) const
-{
-    if (ray.farT() == Ray::infinity()) {
-        pdfForward = pdfBackward = 0.0f;
-        return Vec3f(0.0f);
-    } else if (_absorptionOnly) {
-        pdfForward = pdfBackward = 1.0f;
-        return FastMath::exp(-_sigmaT*ray.farT());
-    } else {
-        Vec3f weight = FastMath::exp(-_sigmaT*ray.farT());
-        pdfForward  =   endOnSurface ? weight.avg() : (_sigmaT*weight).avg();
-        pdfBackward = startOnSurface ? weight.avg() : (_sigmaT*weight).avg();
-        return weight;
-    }
-}
-
-bool HomogeneousMedium::invert(WritablePathSampleGenerator &sampler, const Ray &ray, bool onSurface) const
-{
-    if (_absorptionOnly)
-        return true;
-
-    Vec3f transmittance = std::exp(-ray.farT()*_sigmaT);
-    Vec3f pdfs = _sigmaT*transmittance;
-    float target = sampler.untracked1D()*pdfs.sum();
-    int component = (target < pdfs.x() ? 0 : (target < pdfs.x() + pdfs.y() ? 1 : 2));
-
-    float xi = 1.0f - transmittance[component];
-    if (onSurface)
-        xi += (1.0f - xi)*sampler.untracked1D();
-    sampler.putDiscrete(3, component);
-    sampler.put1D(xi);
-
-    return true;
-}
-
 }
diff --git a/src/core/media/HomogeneousMedium.hpp b/src/core/media/HomogeneousMedium.hpp
@@ -30,12 +30,8 @@ class HomogeneousMedium : public Medium
 
     virtual bool sampleDistance(PathSampleGenerator &sampler, const Ray &ray,
             MediumState &state, MediumSample &sample) const override;
-    virtual bool invertDistance(WritablePathSampleGenerator &sampler, const Ray &ray, bool onSurface) const;
-    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray) const override;
-    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool onSurface) const override;
-    virtual Vec3f transmittanceAndPdfs(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
-            bool endOnSurface, float &pdfForward, float &pdfBackward) const override;
-    virtual bool invert(WritablePathSampleGenerator &sampler, const Ray &ray, bool onSurface) const;
+    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface, bool endOnSurface) const override;
+    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface, bool endOnSurface) const override;
 
     Vec3f sigmaA() const { return _sigmaA; }
     Vec3f sigmaS() const { return _sigmaS; }

diff --git a/src/core/media/Medium.cpp b/src/core/media/Medium.cpp
@@ -1,5 +1,7 @@
 #include "Medium.hpp"
 
+#include "transmittances/ExponentialTransmittance.hpp"
+
 #include "phasefunctions/IsotropicPhaseFunction.hpp"
 
 #include "io/JsonObject.hpp"
@@ -8,7 +10,8 @@
 namespace Tungsten {
 
 Medium::Medium()
-: _phaseFunction(std::make_shared<IsotropicPhaseFunction>()),
+: _transmittance(std::make_shared<ExponentialTransmittance>()),
+  _phaseFunction(std::make_shared<IsotropicPhaseFunction>()),
   _maxBounce(1024)
 {
 }
@@ -19,6 +22,8 @@ void Medium::fromJson(JsonPtr value, const Scene &scene)
 
     if (auto phase = value["phase_function"])
         _phaseFunction = scene.fetchPhase(phase);
+    if (auto trans = value["transmittance"])
+        _transmittance = scene.fetchTransmittance(trans);
 
     value.getField("max_bounces", _maxBounce);
 }
@@ -27,6 +32,7 @@ rapidjson::Value Medium::toJson(Allocator &allocator) const
 {
     return JsonObject{JsonSerializable::toJson(allocator), allocator,
         "phase_function", *_phaseFunction,
+        "transmittance", *_transmittance,
         "max_bounces", _maxBounce
     };
 }
@@ -39,14 +45,19 @@ bool Medium::invertDistance(WritablePathSampleGenerator &/*sampler*/, const Ray
 Vec3f Medium::transmittanceAndPdfs(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
         bool endOnSurface, float &pdfForward, float &pdfBackward) const
 {
-    pdfForward = pdf(sampler, ray, endOnSurface);
-    pdfBackward = pdf(sampler, ray.scatter(ray.hitpoint(), -ray.dir(), 0.0f, ray.farT()), startOnSurface);
-    return transmittance(sampler, ray);
+    pdfForward  = pdf(sampler, ray, startOnSurface, endOnSurface);
+    pdfBackward = pdf(sampler, ray.scatter(ray.hitpoint(), -ray.dir(), 0.0f, ray.farT()), endOnSurface, startOnSurface);
+    return transmittance(sampler, ray, startOnSurface, endOnSurface);
 }
 
 const PhaseFunction *Medium::phaseFunction(const Vec3f &/*p*/) const
 {
     return _phaseFunction.get();
 }
 
+bool Medium::isDirac() const
+{
+    return _transmittance->isDirac();
+}
+
 }
diff --git a/src/core/media/Medium.hpp b/src/core/media/Medium.hpp
@@ -1,6 +1,8 @@
 #ifndef MEDIUM_HPP_
 #define MEDIUM_HPP_
 
+#include "transmittances/Transmittance.hpp"
+
 #include "phasefunctions/PhaseFunction.hpp"
 
 #include "samplerecords/MediumSample.hpp"
@@ -20,6 +22,7 @@ class Scene;
 class Medium : public JsonSerializable
 {
 protected:
+    std::shared_ptr<Transmittance> _transmittance;
     std::shared_ptr<PhaseFunction> _phaseFunction;
     int _maxBounce;
 
@@ -60,11 +63,14 @@ class Medium : public JsonSerializable
     virtual bool sampleDistance(PathSampleGenerator &sampler, const Ray &ray,
             MediumState &state, MediumSample &sample) const = 0;
     virtual bool invertDistance(WritablePathSampleGenerator &sampler, const Ray &ray, bool onSurface) const;
-    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray) const = 0;
-    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool onSurface) const = 0;
+    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
+            bool endOnSurface) const = 0;
+    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface, bool endOnSurface) const = 0;
     virtual Vec3f transmittanceAndPdfs(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
             bool endOnSurface, float &pdfForward, float &pdfBackward) const;
     virtual const PhaseFunction *phaseFunction(const Vec3f &p) const;
+
+    bool isDirac() const;
 };
 
 }

diff --git a/src/core/media/VoxelMedium.cpp b/src/core/media/VoxelMedium.cpp
@@ -121,26 +121,28 @@ bool VoxelMedium::sampleDistance(PathSampleGenerator &sampler, const Ray &ray,
 
     if (_absorptionOnly) {
         sample.t = maxT;
-        sample.weight = _grid->transmittance(sampler, p, w, t0, t1, _sigmaT/wPrime);
+        Vec3f tau = _grid->opticalDepth(sampler, p, w, t0, t1)*(_sigmaT/wPrime);
+        sample.weight = _transmittance->eval(tau, state.firstScatter, true);
         sample.pdf = 1.0f;
         sample.exited = true;
     } else {
         int component = sampler.nextDiscrete(3);
         float sigmaTc = _sigmaT[component];
-        float xi = 1.0f - sampler.next1D();
-        float logXi = -std::log(xi);
+        float tauC = _transmittance->sample(sampler, state.firstScatter)/(sigmaTc/wPrime);
 
-        Vec2f tAndDensity = _grid->inverseOpticalDepth(sampler, p, w, t0, t1, sigmaTc/wPrime, logXi);
+        Vec2f tAndDensity = _grid->inverseOpticalDepth(sampler, p, w, t0, t1, tauC);
         sample.t = tAndDensity.x();
         sample.exited = (sample.t >= t1);
-        float opticalDepth = sample.exited ? tAndDensity.y()/sigmaTc : logXi/sigmaTc;
-        sample.weight = std::exp(-_sigmaT*opticalDepth);
+        if (sample.exited)
+            tauC = tAndDensity.y();
+        Vec3f tau = tauC*(_sigmaT/wPrime);
+        sample.weight = _transmittance->eval(tau, state.firstScatter, sample.exited);
         if (sample.exited) {
-            sample.pdf = sample.weight.avg();
+            sample.pdf = _transmittance->surfaceProbability(tau, state.firstScatter).avg();
         } else {
             float rho = tAndDensity.y();
-            sample.pdf = (rho*_sigmaT*sample.weight).avg();
-            sample.weight *= rho*_sigmaS;
+            sample.pdf = (rho*_sigmaT*_transmittance->mediumPdf(tau, state.firstScatter)).avg();
+            sample.weight *= rho*_sigmaS*_transmittance->sigmaBar();
         }
         sample.weight /= sample.pdf;
         sample.t /= wPrime;
@@ -153,7 +155,8 @@ bool VoxelMedium::sampleDistance(PathSampleGenerator &sampler, const Ray &ray,
     return true;
 }
 
-Vec3f VoxelMedium::transmittance(PathSampleGenerator &sampler, const Ray &ray) const
+Vec3f VoxelMedium::transmittance(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
+        bool endOnSurface) const
 {
     Vec3f p = _worldToGrid*ray.pos();
     Vec3f w = _worldToGrid.transformVector(ray.dir());
@@ -163,10 +166,11 @@ Vec3f VoxelMedium::transmittance(PathSampleGenerator &sampler, const Ray &ray) c
     if (!bboxIntersection(_gridBounds, p, w, t0, t1))
         return Vec3f(1.0f);
 
-    return _grid->transmittance(sampler, p, w, t0, t1, _sigmaT/wPrime);
+    Vec3f tau = _grid->opticalDepth(sampler, p, w, t0, t1)*(_sigmaT/wPrime);
+    return _transmittance->eval(tau, startOnSurface, endOnSurface);
 }
 
-float VoxelMedium::pdf(PathSampleGenerator &sampler, const Ray &ray, bool onSurface) const
+float VoxelMedium::pdf(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface, bool endOnSurface) const
 {
     if (_absorptionOnly) {
         return 1.0f;
@@ -179,38 +183,12 @@ float VoxelMedium::pdf(PathSampleGenerator &sampler, const Ray &ray, bool onSurf
         if (!bboxIntersection(_gridBounds, p, w, t0, t1))
             return 1.0f;
 
-        Vec3f transmittance = _grid->transmittance(sampler, p, w, t0, t1, _sigmaT/wPrime);
-        if (onSurface) {
-            return transmittance.avg();
-        } else {
-            return (_grid->density(p)*_sigmaT*transmittance).avg();
-        }
+        Vec3f tau = _grid->opticalDepth(sampler, p, w, t0, t1)*(_sigmaT/wPrime);
+        if (endOnSurface)
+            return _transmittance->surfaceProbability(tau, startOnSurface).avg();
+        else
+            return (_grid->density(p)*_sigmaT*_transmittance->mediumPdf(tau, startOnSurface)).avg();
     }
 }
 
-Vec3f VoxelMedium::transmittanceAndPdfs(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
-        bool endOnSurface, float &pdfForward, float &pdfBackward) const
-{
-    Vec3f p = _worldToGrid*ray.pos();
-    Vec3f w = _worldToGrid.transformVector(ray.dir());
-    float wPrime = w.length();
-    w /= wPrime;
-    float t0 = 0.0f, t1 = ray.farT()*wPrime;
-    if (!bboxIntersection(_gridBounds, p, w, t0, t1)) {
-        pdfForward = pdfBackward = 1.0f;
-        return Vec3f(1.0f);
-    }
-
-    Vec3f transmittance = _grid->transmittance(sampler, p, w, t0, t1, _sigmaT/wPrime);
-
-    if (_absorptionOnly) {
-        pdfForward = pdfBackward = 1.0f;
-    } else {
-        pdfForward  =   endOnSurface ? transmittance.avg() : (_grid->density(p)*_sigmaT*transmittance).avg();
-        pdfBackward = startOnSurface ? transmittance.avg() : (_grid->density(p)*_sigmaT*transmittance).avg();
-    }
-
-    return transmittance;
-}
-
 }
diff --git a/src/core/media/VoxelMedium.hpp b/src/core/media/VoxelMedium.hpp
@@ -36,10 +36,9 @@ class VoxelMedium : public Medium
 
     virtual bool sampleDistance(PathSampleGenerator &sampler, const Ray &ray,
             MediumState &state, MediumSample &sample) const override;
-    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray) const override;
-    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool onSurface) const override;
-    virtual Vec3f transmittanceAndPdfs(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
-            bool endOnSurface, float &pdfForward, float &pdfBackward) const override;
+    virtual Vec3f transmittance(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface,
+            bool endOnSurface) const override;
+    virtual float pdf(PathSampleGenerator &sampler, const Ray &ray, bool startOnSurface, bool endOnSurface) const override;
 };
 
 }

diff --git a/src/core/transmittances/DavisTransmittance.cpp b/src/core/transmittances/DavisTransmittance.cpp
@@ -0,0 +1,62 @@
+#include "DavisTransmittance.hpp"
+
+#include "io/JsonObject.hpp"
+
+namespace Tungsten {
+
+DavisTransmittance::DavisTransmittance()
+: _alpha(1.1f)
+{
+}
+
+void DavisTransmittance::fromJson(JsonPtr value, const Scene &scene)
+{
+    Transmittance::fromJson(value, scene);
+    value.getField("alpha", _alpha);
+
+    if (_alpha < 1 + 1e-6f) {
+        std::cout << "Warning: alpha parameter of Davis transmittance has to be > 1. Clamping the current value (" << _alpha << ")" << std::endl;
+        _alpha = 1 + 1e-6f;
+    }
+}
+
+rapidjson::Value DavisTransmittance::toJson(Allocator &allocator) const
+{
+    return JsonObject{Transmittance::toJson(allocator), allocator,
+        "type", "davis",
+        "alpha", _alpha
+    };
+}
+
+Vec3f DavisTransmittance::surfaceSurface(const Vec3f &tau) const
+{
+    return std::pow(1.0f + tau/_alpha, -_alpha);
+}
+Vec3f DavisTransmittance::surfaceMedium(const Vec3f &tau) const
+{
+    return std::pow(1.0f + tau/_alpha, -(_alpha + 1.0f));
+}
+Vec3f DavisTransmittance::mediumSurface(const Vec3f &tau) const
+{
+    return surfaceMedium(tau);
+}
+Vec3f DavisTransmittance::mediumMedium(const Vec3f &tau) const
+{
+    return (1.0f + 1.0f/_alpha)*std::pow(1.0f + tau/_alpha, -(_alpha + 2.0f));
+}
+
+float DavisTransmittance::sigmaBar() const
+{
+    return 1.0f;
+}
+
+float DavisTransmittance::sampleSurface(PathSampleGenerator &sampler) const
+{
+    return _alpha*(std::pow(1.0f - sampler.next1D(), -1.0f/_alpha) - 1.0f);
+}
+float DavisTransmittance::sampleMedium(PathSampleGenerator &sampler) const
+{
+    return _alpha*(std::pow(1.0f - sampler.next1D(), -1.0f/(1.0f + _alpha)) - 1.0f);
+}
+
+}
diff --git a/src/core/transmittances/DavisTransmittance.hpp b/src/core/transmittances/DavisTransmittance.hpp
@@ -0,0 +1,31 @@
+#ifndef DAVISTRANSMITTANCE_HPP_
+#define DAVISTRANSMITTANCE_HPP_
+
+#include "Transmittance.hpp"
+
+namespace Tungsten {
+
+class DavisTransmittance : public Transmittance
+{
+    float _alpha;
+
+public:
+    DavisTransmittance();
+
+    virtual void fromJson(JsonPtr value, const Scene &scene) override;
+    virtual rapidjson::Value toJson(Allocator &allocator) const override;
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const override final;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const override final;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const override final;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const override final;
+
+    virtual float sigmaBar() const override final;
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const override final;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const override final;
+};
+
+}
+
+#endif /* DAVISTRANSMITTANCE_HPP_ */
diff --git a/src/core/transmittances/DavisWeinsteinTransmittance.cpp b/src/core/transmittances/DavisWeinsteinTransmittance.cpp
@@ -0,0 +1,119 @@
+#include "DavisWeinsteinTransmittance.hpp"
+
+#include "math/MathUtil.hpp"
+
+#include "io/JsonObject.hpp"
+
+namespace Tungsten {
+
+DavisWeinsteinTransmittance::DavisWeinsteinTransmittance()
+: _h(0.75f),
+  _c(1.0f)
+{
+}
+
+void DavisWeinsteinTransmittance::fromJson(JsonPtr value, const Scene &scene)
+{
+    Transmittance::fromJson(value, scene);
+    value.getField("h", _h);
+    value.getField("c", _c);
+    if (_h < 0.5f || _h > 1.0f)
+        std::cout << "Warning: Valid range of Davis Weinstein transmittance is [0.5, 1.0]. Clamping current value (" << _h << ") to within range" << std::endl;
+    _h = clamp(_h, 0.5f, 1.0f);
+}
+
+rapidjson::Value DavisWeinsteinTransmittance::toJson(Allocator &allocator) const
+{
+    return JsonObject{Transmittance::toJson(allocator), allocator,
+        "type", "davis",
+        "h", _h,
+        "c", _c
+    };
+}
+
+float DavisWeinsteinTransmittance::computeAlpha(float tau) const
+{
+    float beta = 2.0f*_h - 1.0f;
+    return std::pow(tau, 1 - beta)/(std::pow(_c, 1 + beta));
+}
+
+Vec3f DavisWeinsteinTransmittance::surfaceSurface(const Vec3f &tau) const
+{
+    float alpha = computeAlpha(tau[0]);
+    float Tr = std::pow(1.0f + tau[0]/alpha, -alpha);
+
+    return Vec3f(std::isnan(Tr) ? 0 : Tr);
+}
+Vec3f DavisWeinsteinTransmittance::surfaceMedium(const Vec3f &tau) const
+{
+    float beta = 2.0f*_h - 1.0f;
+    float t = tau[0];
+    float alpha = computeAlpha(t);
+    float base = 1.0f + t/alpha;
+    float trSurface = std::pow(base, -alpha);
+
+    float Tr = trSurface*(beta/base - (beta - 1.0f)*alpha/t*std::log(base));
+
+    return Vec3f(std::isnan(Tr) ? 0 : Tr);
+}
+Vec3f DavisWeinsteinTransmittance::mediumSurface(const Vec3f &tau) const
+{
+    return surfaceMedium(tau);
+}
+Vec3f DavisWeinsteinTransmittance::mediumMedium(const Vec3f &tau) const
+{
+    float beta = 2.0f*_h - 1.0f;
+    float t = tau[0];
+    float alpha = computeAlpha(t);
+    float base = 1.0f + t/alpha;
+    float logBase = std::log(base);
+    float trSurface = std::pow(base, -alpha);
+    float term1 = beta*(-1.0f + beta*(1.0f + t) + (-1.0f + 2.0f*beta)*t/alpha)/(t*base*base);
+    float term2 = ((-1.0f + beta)*beta*alpha/(t*t)*(2.0f*t + base)*logBase)/base;
+    float term3 = (beta - 1.0f)*alpha/t*logBase;
+
+    float Tr = trSurface*(term1 - term2 + term3*term3);
+
+    return Vec3f(std::isnan(Tr) ? 0 : Tr);
+}
+
+float DavisWeinsteinTransmittance::sigmaBar() const
+{
+    return 1.0f;
+}
+
+// We have no real way of analytically sampling this function, so a simple bisection will have to do
+float DavisWeinsteinTransmittance::sampleSurface(PathSampleGenerator &sampler) const
+{
+    float xi = sampler.next1D();
+    auto cdf = [this](float tau) { return 1.0f - surfaceSurface(Vec3f(tau))[0]; };
+    float step = 1e6;
+    float result = step*2;
+    while (step > 1e-6) {
+        if (cdf(result) > xi)
+            result -= step;
+        else
+            result += step;
+        step /= 2;
+    }
+
+    return result;
+}
+float DavisWeinsteinTransmittance::sampleMedium(PathSampleGenerator &sampler) const
+{
+    float xi = sampler.next1D();
+    auto cdf = [this](float tau) { return 1.0f - mediumSurface(Vec3f(tau))[0]; };
+    float step = 1e6;
+    float result = step * 2;
+    while (step > 1e-6) {
+        if (cdf(result) > xi)
+            result -= step;
+        else
+            result += step;
+        step /= 2;
+    }
+
+    return result;
+}
+
+}
diff --git a/src/core/transmittances/DavisWeinsteinTransmittance.hpp b/src/core/transmittances/DavisWeinsteinTransmittance.hpp
@@ -0,0 +1,34 @@
+#ifndef DAVISWEINSTEINTRANSMITTANCE_HPP_
+#define DAVISWEINSTEINTRANSMITTANCE_HPP_
+
+#include "Transmittance.hpp"
+
+namespace Tungsten {
+
+class DavisWeinsteinTransmittance : public Transmittance
+{
+    float _h;
+    float _c;
+
+    float computeAlpha(float tau) const;
+
+public:
+    DavisWeinsteinTransmittance();
+
+    virtual void fromJson(JsonPtr value, const Scene &scene) override;
+    virtual rapidjson::Value toJson(Allocator &allocator) const override;
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const override final;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const override final;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const override final;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const override final;
+
+    virtual float sigmaBar() const override final;
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const override final;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const override final;
+};
+
+}
+
+#endif /* DAVISWEINSTEINTRANSMITTANCE_HPP_ */
diff --git a/src/core/transmittances/DoubleExponentialTransmittance.cpp b/src/core/transmittances/DoubleExponentialTransmittance.cpp
@@ -0,0 +1,68 @@
+#include "DoubleExponentialTransmittance.hpp"
+
+#include "sampling/UniformPathSampler.hpp"
+
+#include "math/MathUtil.hpp"
+
+#include "io/JsonObject.hpp"
+
+#include "Memory.hpp"
+
+namespace Tungsten {
+
+DoubleExponentialTransmittance::DoubleExponentialTransmittance()
+: _sigmaA(0.5f),
+  _sigmaB(10.0f)
+{
+}
+
+void DoubleExponentialTransmittance::fromJson(JsonPtr value, const Scene &scene)
+{
+    Transmittance::fromJson(value, scene);
+    value.getField("sigma_a", _sigmaA);
+    value.getField("sigma_b", _sigmaB);
+}
+
+rapidjson::Value DoubleExponentialTransmittance::toJson(Allocator &allocator) const
+{
+    return JsonObject{Transmittance::toJson(allocator), allocator,
+        "type", "double_exponential",
+        "sigma_a", _sigmaA,
+        "sigma_b", _sigmaB,
+    };
+}
+
+Vec3f DoubleExponentialTransmittance::surfaceSurface(const Vec3f &tau) const
+{
+    return 0.5f*(std::exp(-_sigmaA*tau) + std::exp(-_sigmaB*tau));
+}
+Vec3f DoubleExponentialTransmittance::surfaceMedium(const Vec3f &tau) const
+{
+    return 0.5f*(_sigmaA*std::exp(-_sigmaA*tau) + _sigmaB*std::exp(-_sigmaB*tau));
+}
+Vec3f DoubleExponentialTransmittance::mediumSurface(const Vec3f &tau) const
+{
+    return (_sigmaA*std::exp(-_sigmaA*tau) + _sigmaB*std::exp(-_sigmaB*tau))/(_sigmaA + _sigmaB);
+}
+Vec3f DoubleExponentialTransmittance::mediumMedium(const Vec3f &tau) const
+{
+    return (sqr(_sigmaA)*std::exp(-_sigmaA*tau) + sqr(_sigmaB)*std::exp(-_sigmaB*tau))/(_sigmaA + _sigmaB);
+}
+
+float DoubleExponentialTransmittance::sigmaBar() const
+{
+    return 0.5f*(_sigmaA + _sigmaB);
+}
+
+float DoubleExponentialTransmittance::sampleSurface(PathSampleGenerator &sampler) const
+{
+    float t = -std::log(1.0f - sampler.next1D());
+    return sampler.nextBoolean(0.5f) ? t/_sigmaA : t/_sigmaB;
+}
+float DoubleExponentialTransmittance::sampleMedium(PathSampleGenerator &sampler) const
+{
+    float t = -std::log(1.0f - sampler.next1D());
+    return sampler.nextBoolean(_sigmaA/(_sigmaA + _sigmaB)) ? t/_sigmaA : t/_sigmaB;
+}
+
+}
diff --git a/src/core/transmittances/DoubleExponentialTransmittance.hpp b/src/core/transmittances/DoubleExponentialTransmittance.hpp
@@ -0,0 +1,35 @@
+#ifndef DoubleExponentialTransmittance_HPP_
+#define DoubleExponentialTransmittance_HPP_
+
+#include "Transmittance.hpp"
+
+namespace Tungsten {
+
+class DoubleExponentialTransmittance : public Transmittance
+{
+    float _sigmaA;
+    float _sigmaB;
+
+public:
+    DoubleExponentialTransmittance();
+
+    virtual void fromJson(JsonPtr value, const Scene &scene) override;
+    virtual rapidjson::Value toJson(Allocator &allocator) const override;
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const override final;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const override final;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const override final;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const override final;
+
+    virtual float sigmaBar() const override final;
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const override final;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const override final;
+};
+
+
+}
+
+
+
+#endif /* DoubleExponentialTransmittance_HPP_ */
diff --git a/src/core/transmittances/ErlangTransmittance.cpp b/src/core/transmittances/ErlangTransmittance.cpp
@@ -0,0 +1,69 @@
+#include "ErlangTransmittance.hpp"
+
+#include "sampling/UniformPathSampler.hpp"
+
+#include "math/MathUtil.hpp"
+
+#include "io/JsonObject.hpp"
+
+#include "Memory.hpp"
+
+namespace Tungsten {
+
+ErlangTransmittance::ErlangTransmittance()
+: _lambda(5.0f)
+{
+}
+
+void ErlangTransmittance::fromJson(JsonPtr value, const Scene &scene)
+{
+    Transmittance::fromJson(value, scene);
+    value.getField("rate", _lambda);
+}
+
+rapidjson::Value ErlangTransmittance::toJson(Allocator &allocator) const
+{
+    return JsonObject{Transmittance::toJson(allocator), allocator,
+        "type", "erlang",
+        "rate", _lambda,
+    };
+}
+
+Vec3f ErlangTransmittance::surfaceSurface(const Vec3f &tau) const
+{
+    return 0.5f*std::exp(-_lambda*tau)*(2.0f + _lambda*tau);
+}
+Vec3f ErlangTransmittance::surfaceMedium(const Vec3f &tau) const
+{
+    return mediumSurface(tau)*_lambda*0.5f;
+}
+Vec3f ErlangTransmittance::mediumSurface(const Vec3f &tau) const
+{
+    return std::exp(-_lambda*tau)*(1.0f + _lambda*tau);
+}
+Vec3f ErlangTransmittance::mediumMedium(const Vec3f &tau) const
+{
+    return sqr(_lambda)*tau*std::exp(-_lambda*tau);
+}
+
+float ErlangTransmittance::sigmaBar() const
+{
+    return _lambda*0.5f;
+}
+
+float ErlangTransmittance::sampleSurface(PathSampleGenerator &sampler) const
+{
+    float xi = sampler.next1D();
+    float x = 0.5f;
+    for (int i = 0; i < 10; ++i) {
+        x += (xi - (1.0f - surfaceSurface(Vec3f(x))[0]))/surfaceMedium(Vec3f(x))[0];
+        x = max(x, 0.0f);
+    }
+    return x;
+}
+float ErlangTransmittance::sampleMedium(PathSampleGenerator &sampler) const
+{
+    return -1.0f/_lambda*std::log(sampler.next1D()*sampler.next1D());
+}
+
+}
diff --git a/src/core/transmittances/ErlangTransmittance.hpp b/src/core/transmittances/ErlangTransmittance.hpp
@@ -0,0 +1,34 @@
+#ifndef ErlangTransmittance_HPP_
+#define ErlangTransmittance_HPP_
+
+#include "Transmittance.hpp"
+
+namespace Tungsten {
+
+class ErlangTransmittance : public Transmittance
+{
+    float _lambda;
+
+public:
+    ErlangTransmittance();
+
+    virtual void fromJson(JsonPtr value, const Scene &scene) override;
+    virtual rapidjson::Value toJson(Allocator &allocator) const override;
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const override final;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const override final;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const override final;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const override final;
+
+    virtual float sigmaBar() const override final;
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const override final;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const override final;
+};
+
+
+}
+
+
+
+#endif /* ErlangTransmittance_HPP_ */
diff --git a/src/core/transmittances/ExponentialTransmittance.cpp b/src/core/transmittances/ExponentialTransmittance.cpp
@@ -0,0 +1,57 @@
+#include "ExponentialTransmittance.hpp"
+
+#include "sampling/UniformPathSampler.hpp"
+
+#include "math/MathUtil.hpp"
+#include "math/FastMath.hpp"
+
+#include "io/JsonObject.hpp"
+
+#include "Memory.hpp"
+
+namespace Tungsten {
+
+void ExponentialTransmittance::fromJson(JsonPtr value, const Scene &scene)
+{
+    Transmittance::fromJson(value, scene);
+}
+
+rapidjson::Value ExponentialTransmittance::toJson(Allocator &allocator) const
+{
+    return JsonObject{Transmittance::toJson(allocator), allocator,
+        "type", "exponential",
+    };
+}
+
+Vec3f ExponentialTransmittance::surfaceSurface(const Vec3f &tau) const
+{
+    return FastMath::exp(-tau);
+}
+Vec3f ExponentialTransmittance::surfaceMedium(const Vec3f &tau) const
+{
+    return FastMath::exp(-tau);
+}
+Vec3f ExponentialTransmittance::mediumSurface(const Vec3f &tau) const
+{
+    return FastMath::exp(-tau);
+}
+Vec3f ExponentialTransmittance::mediumMedium(const Vec3f &tau) const
+{
+    return FastMath::exp(-tau);
+}
+
+float ExponentialTransmittance::sigmaBar() const
+{
+    return 1.0f;
+}
+
+float ExponentialTransmittance::sampleSurface(PathSampleGenerator &sampler) const
+{
+    return -std::log(1.0f - sampler.next1D());
+}
+float ExponentialTransmittance::sampleMedium(PathSampleGenerator &sampler) const
+{
+    return -std::log(1.0f - sampler.next1D());
+}
+
+}
diff --git a/src/core/transmittances/ExponentialTransmittance.hpp b/src/core/transmittances/ExponentialTransmittance.hpp
@@ -0,0 +1,32 @@
+#ifndef EXPONENTIALTRANSMITTANCE_HPP_
+#define EXPONENTIALTRANSMITTANCE_HPP_
+
+#include "Transmittance.hpp"
+
+namespace Tungsten {
+
+class ExponentialTransmittance : public Transmittance
+{
+public:
+    ExponentialTransmittance() = default;
+
+    virtual void fromJson(JsonPtr value, const Scene &scene) override;
+    virtual rapidjson::Value toJson(Allocator &allocator) const override;
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const override final;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const override final;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const override final;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const override final;
+
+    virtual float sigmaBar() const override final;
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const override final;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const override final;
+};
+
+
+}
+
+
+
+#endif /* EXPONENTIALTRANSMITTANCE_HPP_ */
diff --git a/src/core/transmittances/InterpolatedTransmittance.cpp b/src/core/transmittances/InterpolatedTransmittance.cpp
@@ -0,0 +1,85 @@
+#include "InterpolatedTransmittance.hpp"
+
+#include "ExponentialTransmittance.hpp"
+#include "LinearTransmittance.hpp"
+#include "ErlangTransmittance.hpp"
+
+#include "sampling/UniformPathSampler.hpp"
+
+#include "math/MathUtil.hpp"
+
+#include "io/JsonObject.hpp"
+#include "io/Scene.hpp"
+
+#include "Memory.hpp"
+
+namespace Tungsten {
+
+InterpolatedTransmittance::InterpolatedTransmittance()
+: _trA(std::make_shared<LinearTransmittance>()),
+  _trB(std::make_shared<ErlangTransmittance>()),
+  _u(0.5f)
+{
+}
+
+void InterpolatedTransmittance::fromJson(JsonPtr value, const Scene &scene)
+{
+    Transmittance::fromJson(value, scene);
+    if (auto m = value["tr_a"]) _trA = scene.fetchTransmittance(m);
+    if (auto m = value["tr_b"]) _trB = scene.fetchTransmittance(m);
+    value.getField("ratio", _u);
+}
+
+rapidjson::Value InterpolatedTransmittance::toJson(Allocator &allocator) const
+{
+    return JsonObject{Transmittance::toJson(allocator), allocator,
+        "type", "interpolated",
+        "tr_a", _trA->toJson(allocator),
+        "tr_b", _trB->toJson(allocator),
+        "ratio", _u
+    };
+}
+
+Vec3f InterpolatedTransmittance::surfaceSurface(const Vec3f &tau) const
+{
+    return sigmaBar()*lerp(_trA->surfaceSurface(tau)/_trA->sigmaBar(), _trB->surfaceSurface(tau)/_trB->sigmaBar(), _u);
+}
+Vec3f InterpolatedTransmittance::surfaceMedium(const Vec3f &tau) const
+{
+    return mediumSurface(tau)*sigmaBar();
+}
+Vec3f InterpolatedTransmittance::mediumSurface(const Vec3f &tau) const
+{
+    return lerp(_trA->mediumSurface(tau), _trB->mediumSurface(tau), _u);
+}
+Vec3f InterpolatedTransmittance::mediumMedium(const Vec3f &tau) const
+{
+    Vec3f pdfA = _trA->mediumMedium(tau);
+    Vec3f pdfB = _trB->mediumMedium(tau);
+    Vec3f result;
+    for (int i = 0; i < 3; ++i) {
+        bool diracA = _trA->isDirac() && pdfA[i] > 0.0f;
+        bool diracB = _trB->isDirac() && pdfB[i] > 0.0f;
+        if (diracA ^ diracB)
+            result[i] = diracA ? pdfA[i] : pdfB[i];
+        else
+            result[i] = lerp(pdfA[i], pdfB[i], _u);
+    }
+    return result;
+}
+
+float InterpolatedTransmittance::sigmaBar() const
+{
+    return 1.0f/lerp(1.0f/_trA->sigmaBar(), 1.0f/_trB->sigmaBar(), _u);
+}
+
+float InterpolatedTransmittance::sampleSurface(PathSampleGenerator &sampler) const
+{
+    return sampler.nextBoolean(_u) ? _trB->sampleSurface(sampler) : _trA->sampleSurface(sampler);
+}
+float InterpolatedTransmittance::sampleMedium(PathSampleGenerator &sampler) const
+{
+    return sampler.nextBoolean(_u) ? _trB->sampleMedium(sampler) : _trA->sampleMedium(sampler);
+}
+
+}
diff --git a/src/core/transmittances/InterpolatedTransmittance.hpp b/src/core/transmittances/InterpolatedTransmittance.hpp
@@ -0,0 +1,39 @@
+#ifndef INTERPOLATEDTRANSMITTANCE_HPP_
+#define INTERPOLATEDTRANSMITTANCE_HPP_
+
+#include "Transmittance.hpp"
+
+#include <vector>
+#include <memory>
+
+namespace Tungsten {
+
+class InterpolatedTransmittance : public Transmittance
+{
+    std::shared_ptr<Transmittance> _trA;
+    std::shared_ptr<Transmittance> _trB;
+    float _u;
+
+public:
+    InterpolatedTransmittance();
+
+    virtual void fromJson(JsonPtr value, const Scene &scene) override;
+    virtual rapidjson::Value toJson(Allocator &allocator) const override;
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const override final;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const override final;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const override final;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const override final;
+
+    virtual float sigmaBar() const override final;
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const override final;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const override final;
+};
+
+
+}
+
+
+
+#endif /* INTERPOLATEDTRANSMITTANCE_HPP_ */
diff --git a/src/core/transmittances/LinearTransmittance.cpp b/src/core/transmittances/LinearTransmittance.cpp
@@ -0,0 +1,79 @@
+#include "LinearTransmittance.hpp"
+
+#include "sampling/UniformPathSampler.hpp"
+
+#include "math/MathUtil.hpp"
+
+#include "io/JsonObject.hpp"
+
+#include "Memory.hpp"
+
+namespace Tungsten {
+
+LinearTransmittance::LinearTransmittance()
+: _maxT(1.0f)
+{
+}
+
+void LinearTransmittance::fromJson(JsonPtr value, const Scene &scene)
+{
+    Transmittance::fromJson(value, scene);
+    value.getField("max_t", _maxT);
+}
+
+rapidjson::Value LinearTransmittance::toJson(Allocator &allocator) const
+{
+    return JsonObject{Transmittance::toJson(allocator), allocator,
+        "type", "linear",
+        "max_t", _maxT
+    };
+}
+
+Vec3f LinearTransmittance::surfaceSurface(const Vec3f &tau) const
+{
+    return 1.0f - min(tau/_maxT, Vec3f(1.0f));
+}
+Vec3f LinearTransmittance::surfaceMedium(const Vec3f &tau) const
+{
+    Vec3f result(1.0f/_maxT);
+    for (int i = 0; i < 3; ++i)
+        if (tau[i] > _maxT)
+            result[i] = 0.0f;
+    return result;
+}
+Vec3f LinearTransmittance::mediumSurface(const Vec3f &tau) const
+{
+    Vec3f result(1.0f);
+    for (int i = 0; i < 3; ++i)
+        if (tau[i] > _maxT)
+            result[i] = 0.0f;
+    return result;
+}
+Vec3f LinearTransmittance::mediumMedium(const Vec3f &tau) const
+{
+    Vec3f result;
+    for (int i = 0; i < 3; ++i)
+        result[i] = std::abs(tau[i] - _maxT) < 1e-3f ? 1.0f : 0.0f;
+    return result;
+}
+
+float LinearTransmittance::sigmaBar() const
+{
+    return 1.0f/_maxT;
+}
+
+bool LinearTransmittance::isDirac() const
+{
+    return true;
+}
+
+float LinearTransmittance::sampleSurface(PathSampleGenerator &sampler) const
+{
+    return _maxT*sampler.next1D();
+}
+float LinearTransmittance::sampleMedium(PathSampleGenerator &/*sampler*/) const
+{
+    return _maxT;
+}
+
+}
diff --git a/src/core/transmittances/LinearTransmittance.hpp b/src/core/transmittances/LinearTransmittance.hpp
@@ -0,0 +1,33 @@
+#ifndef LINEARTRANSMITTANCE_HPP_
+#define LINEARTRANSMITTANCE_HPP_
+
+#include "Transmittance.hpp"
+
+namespace Tungsten {
+
+class LinearTransmittance : public Transmittance
+{
+    float _maxT;
+
+public:
+    LinearTransmittance();
+
+    virtual void fromJson(JsonPtr value, const Scene &scene) override;
+    virtual rapidjson::Value toJson(Allocator &allocator) const override;
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const override final;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const override final;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const override final;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const override final;
+
+    virtual bool isDirac() const override final;
+
+    virtual float sigmaBar() const override final;
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const override final;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const override final;
+};
+
+}
+
+#endif /* LINEARTRANSMITTANCE_HPP_ */
diff --git a/src/core/transmittances/PulseTransmittance.cpp b/src/core/transmittances/PulseTransmittance.cpp
@@ -0,0 +1,108 @@
+#include "PulseTransmittance.hpp"
+
+#include "sampling/UniformPathSampler.hpp"
+
+#include "math/MathUtil.hpp"
+
+#include "io/JsonObject.hpp"
+#include "io/FileUtils.hpp"
+
+#include "Memory.hpp"
+
+namespace Tungsten {
+
+PulseTransmittance::PulseTransmittance()
+: _a(0.0f),
+  _b(1.0f),
+  _numPulses(4)
+{
+}
+
+void PulseTransmittance::fromJson(JsonPtr value, const Scene &scene)
+{
+    Transmittance::fromJson(value, scene);
+    value.getField("min", _a);
+    value.getField("max", _b);
+    value.getField("num_pulses", _numPulses);
+}
+
+rapidjson::Value PulseTransmittance::toJson(Allocator &allocator) const
+{
+    return JsonObject{Transmittance::toJson(allocator), allocator,
+        "type", "pulse",
+        "min", _a,
+        "max", _b,
+        "num_pulses", _numPulses,
+    };
+}
+
+bool PulseTransmittance::isDirac() const
+{
+    return true;
+}
+
+Vec3f PulseTransmittance::surfaceSurface(const Vec3f &tau) const
+{
+    Vec3f idxF = clamp(float(_numPulses)*(tau - _a)/(_b - _a) + 0.5f, Vec3f(0.0f), Vec3f(_numPulses));
+    Vec3i idx = Vec3i(idxF);
+
+    Vec3f height = Vec3f(_numPulses - idx)/float(_numPulses);
+
+    Vec3f cellIntegral = height*(idxF - Vec3f(idx));
+    for (int i = 0; i < 3; ++i) {
+        if (idx[i] > 0)
+            cellIntegral[i] += (idx[i] - 0.5f) - (idx[i]*(idx[i] - 1))/(2.0f*_numPulses);
+        else
+            cellIntegral[i] -= 0.5f;
+    }
+
+    return 1.0f - (2.0f/_numPulses)*cellIntegral;
+}
+Vec3f PulseTransmittance::surfaceMedium(const Vec3f &tau) const
+{
+    return 2.0f/(_b - _a)*mediumSurface(tau);
+}
+Vec3f PulseTransmittance::mediumSurface(const Vec3f &tau) const
+{
+    Vec3i idx = clamp(Vec3i(float(_numPulses)*(tau - _a)/(_b - _a) + 0.5f), Vec3i(0), Vec3i(_numPulses));
+    return 1.0f - Vec3f(idx)/float(_numPulses);
+}
+Vec3f PulseTransmittance::mediumMedium(const Vec3f &tau) const
+{
+    Vec3f idxF = clamp(float(_numPulses)*(tau - _a)/(_b - _a), Vec3f(0.0f), Vec3f(_numPulses));
+    Vec3i idx = Vec3i(idxF);
+    return (1.0f/_numPulses)*Vec3f(
+        std::abs(idxF[0] - idx[0] - 0.5f) < 1e-3f ? 1.0f : 0.0f,
+        std::abs(idxF[1] - idx[1] - 0.5f) < 1e-3f ? 1.0f : 0.0f,
+        std::abs(idxF[2] - idx[2] - 0.5f) < 1e-3f ? 1.0f : 0.0f
+    );
+}
+
+float PulseTransmittance::sigmaBar() const
+{
+    return 2.0f/(_b - _a);
+}
+
+float PulseTransmittance::sampleSurface(PathSampleGenerator &sampler) const
+{
+    float xi = sampler.next1D()*_numPulses*0.5f;
+    float delta = 1.0f/_numPulses;
+
+    for (int i = 0; i < _numPulses; ++i) {
+        float h0 = 1.0f - (i + 0.0f)*delta;
+        float h1 = 1.0f - (i + 1.0f)*delta;
+        xi -= h0*0.5f;
+        if (xi < 0.0f)
+            return _a + (i + 0.0f + 0.5f*sampler.next1D())*(_b - _a)*delta;
+        xi -= h1*0.5f;
+        if (xi < 0.0f)
+            return _a + (i + 0.5f + 0.5f*sampler.next1D())*(_b - _a)*delta;
+    }
+    return 0.0f;
+}
+float PulseTransmittance::sampleMedium(PathSampleGenerator &sampler) const
+{
+    return _a + (0.5f + int(sampler.next1D()*_numPulses))/float(_numPulses)*(_b - _a);
+}
+
+}
diff --git a/src/core/transmittances/PulseTransmittance.hpp b/src/core/transmittances/PulseTransmittance.hpp
@@ -0,0 +1,37 @@
+#ifndef PULSETRANSMITTANCE_HPP_
+#define PULSETRANSMITTANCE_HPP_
+
+#include "Transmittance.hpp"
+
+namespace Tungsten {
+
+class PulseTransmittance : public Transmittance
+{
+    float _a, _b;
+    int _numPulses;
+
+public:
+    PulseTransmittance();
+
+    virtual void fromJson(JsonPtr value, const Scene &scene) override;
+    virtual rapidjson::Value toJson(Allocator &allocator) const override;
+
+    virtual bool isDirac() const override final;
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const override final;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const override final;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const override final;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const override final;
+
+    virtual float sigmaBar() const override final;
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const override final;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const override final;
+};
+
+
+}
+
+
+
+#endif /* PULSETRANSMITTANCE_HPP_ */
diff --git a/src/core/transmittances/QuadraticTransmittance.cpp b/src/core/transmittances/QuadraticTransmittance.cpp
@@ -0,0 +1,68 @@
+#include "QuadraticTransmittance.hpp"
+
+#include "sampling/UniformPathSampler.hpp"
+
+#include "math/MathUtil.hpp"
+
+#include "io/JsonObject.hpp"
+
+#include "Memory.hpp"
+
+namespace Tungsten {
+
+QuadraticTransmittance::QuadraticTransmittance()
+: _maxT(0.75f)
+{
+}
+
+void QuadraticTransmittance::fromJson(JsonPtr value, const Scene &scene)
+{
+    Transmittance::fromJson(value, scene);
+    value.getField("max_t", _maxT);
+}
+
+rapidjson::Value QuadraticTransmittance::toJson(Allocator &allocator) const
+{
+    return JsonObject{Transmittance::toJson(allocator), allocator,
+        "type", "quadratic",
+        "max_t", _maxT
+    };
+}
+
+Vec3f QuadraticTransmittance::surfaceSurface(const Vec3f &tau) const
+{
+    Vec3f t = min(tau/_maxT, Vec3f(1.0f));
+    return 1.0f - 2.0f*t + t*t;
+}
+Vec3f QuadraticTransmittance::surfaceMedium(const Vec3f &tau) const
+{
+    return (2.0f/_maxT)*(1.0f - min(tau/_maxT, Vec3f(1.0f)));
+}
+Vec3f QuadraticTransmittance::mediumSurface(const Vec3f &tau) const
+{
+    return 1.0f - min(tau/_maxT, Vec3f(1.0f));
+}
+Vec3f QuadraticTransmittance::mediumMedium(const Vec3f &tau) const
+{
+    Vec3f result(1.0f/_maxT);
+    for (int i = 0; i < 3; ++i)
+        if (tau[i] > _maxT)
+            result[i] = 0.0f;
+    return result;
+}
+
+float QuadraticTransmittance::sigmaBar() const
+{
+    return 2.0f/_maxT;
+}
+
+float QuadraticTransmittance::sampleSurface(PathSampleGenerator &sampler) const
+{
+    return _maxT*(1.0f - std::sqrt(1.0f - sampler.next1D()));
+}
+float QuadraticTransmittance::sampleMedium(PathSampleGenerator &sampler) const
+{
+    return _maxT*sampler.next1D();
+}
+
+}
diff --git a/src/core/transmittances/QuadraticTransmittance.hpp b/src/core/transmittances/QuadraticTransmittance.hpp
@@ -0,0 +1,31 @@
+#ifndef QuadraticTransmittance_HPP_
+#define QuadraticTransmittance_HPP_
+
+#include "Transmittance.hpp"
+
+namespace Tungsten {
+
+class QuadraticTransmittance : public Transmittance
+{
+    float _maxT;
+
+public:
+    QuadraticTransmittance();
+
+    virtual void fromJson(JsonPtr value, const Scene &scene) override;
+    virtual rapidjson::Value toJson(Allocator &allocator) const override;
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const override final;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const override final;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const override final;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const override final;
+
+    virtual float sigmaBar() const override final;
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const override final;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const override final;
+};
+
+}
+
+#endif /* QuadraticTransmittance_HPP_ */
diff --git a/src/core/transmittances/Transmittance.hpp b/src/core/transmittances/Transmittance.hpp
@@ -0,0 +1,62 @@
+#ifndef TRANSMITTANCE_HPP_
+#define TRANSMITTANCE_HPP_
+
+#include "sampling/PathSampleGenerator.hpp"
+
+#include "math/Vec.hpp"
+
+#include "io/JsonSerializable.hpp"
+
+#include "StringableEnum.hpp"
+
+#include <algorithm>
+#include <array>
+
+namespace Tungsten {
+
+class Transmittance : public JsonSerializable
+{
+public:
+    virtual ~Transmittance() {}
+
+    inline Vec3f eval(const Vec3f &tau, bool startOnSurface, bool endOnSurface) const
+    {
+        if (startOnSurface && endOnSurface)
+            return surfaceSurface(tau);
+        else if (!startOnSurface && !endOnSurface)
+            return mediumMedium(tau)/sigmaBar();
+        else
+            return mediumSurface(tau);
+    }
+    inline float sample(PathSampleGenerator &sampler, bool startOnSurface) const
+    {
+        return startOnSurface ? sampleSurface(sampler) : sampleMedium(sampler);
+    }
+    inline Vec3f surfaceProbability(const Vec3f &tau, bool startOnSurface) const
+    {
+        return startOnSurface ? surfaceSurface(tau) : mediumSurface(tau);
+    }
+    inline Vec3f mediumPdf(const Vec3f &tau, bool startOnSurface) const
+    {
+        return startOnSurface ? surfaceMedium(tau) : mediumMedium(tau);
+    }
+
+    virtual bool isDirac() const // Returns true if mediumMedium is a dirac/sum of dirac deltas
+    {
+        return false;
+    }
+
+    virtual Vec3f surfaceSurface(const Vec3f &tau) const = 0;
+    virtual Vec3f surfaceMedium(const Vec3f &tau) const = 0;
+    virtual Vec3f mediumSurface(const Vec3f &tau) const = 0;
+    virtual Vec3f mediumMedium(const Vec3f &tau) const = 0;
+
+    virtual float sigmaBar() const = 0; // Returns surfaceMedium(x)/mediumSurface(x) (i.e. identical to surfaceMedium(0))
+
+    virtual float sampleSurface(PathSampleGenerator &sampler) const = 0;
+    virtual float sampleMedium(PathSampleGenerator &sampler) const = 0;
+};
+
+}
+
+#endif /* TRANSMITTANCE_HPP_ */
diff --git a/src/core/transmittances/TransmittanceFactory.cpp b/src/core/transmittances/TransmittanceFactory.cpp
@@ -0,0 +1,26 @@
+#include "TransmittanceFactory.hpp"
+#include "DoubleExponentialTransmittance.hpp"
+#include "DavisWeinsteinTransmittance.hpp"
+#include "InterpolatedTransmittance.hpp"
+#include "ExponentialTransmittance.hpp"
+#include "QuadraticTransmittance.hpp"
+#include "ErlangTransmittance.hpp"
+#include "LinearTransmittance.hpp"
+#include "DavisTransmittance.hpp"
+#include "PulseTransmittance.hpp"
+
+namespace Tungsten {
+
+DEFINE_STRINGABLE_ENUM(TransmittanceFactory, "transmittance", ({
+    {"double_exponential", std::make_shared<DoubleExponentialTransmittance>},
+    {"exponential", std::make_shared<ExponentialTransmittance>},
+    {"quadratic", std::make_shared<QuadraticTransmittance>},
+    {"linear", std::make_shared<LinearTransmittance>},
+    {"erlang", std::make_shared<ErlangTransmittance>},
+    {"davis", std::make_shared<DavisTransmittance>},
+    {"davis_weinstein", std::make_shared<DavisWeinsteinTransmittance>},
+    {"pulse", std::make_shared<PulseTransmittance>},
+    {"interpolated", std::make_shared<InterpolatedTransmittance>},
+}))
+
+}
diff --git a/src/core/transmittances/TransmittanceFactory.hpp b/src/core/transmittances/TransmittanceFactory.hpp
@@ -0,0 +1,17 @@
+#ifndef TRANSMITTANCEFACTORY_HPP_
+#define TRANSMITTANCEFACTORY_HPP_
+
+#include "StringableEnum.hpp"
+
+#include <functional>
+#include <memory>
+
+namespace Tungsten {
+
+class Transmittance;
+
+typedef StringableEnum<std::function<std::shared_ptr<Transmittance>()>> TransmittanceFactory;
+
+}
+
+#endif /* TRANSMITTANCEFACTORY_HPP_ */