diff --git a/internal/ceres/autodiff_benchmarks/CMakeLists.txt b/internal/ceres/autodiff_benchmarks/CMakeLists.txt
index e42e584..29a238c 100644
--- a/internal/ceres/autodiff_benchmarks/CMakeLists.txt
+++ b/internal/ceres/autodiff_benchmarks/CMakeLists.txt
@@ -13,6 +13,7 @@
   generate_benchmark_functor(SnavelyReprojectionError snavely_reprojection_error.h)
   generate_benchmark_functor(Linear1CostFunction linear_cost_functions.h)
   generate_benchmark_functor(Linear10CostFunction linear_cost_functions.h)
+  generate_benchmark_functor(Brdf brdf_cost_function.h)
 
   list(APPEND CERES_BENCHMARK_FLAGS "-DWITH_CODE_GENERATION")
 endif()
@@ -46,10 +47,20 @@
 add_dependencies_to_benchmark(linear_benchmark_fast_math)
 target_compile_options(linear_benchmark_fast_math PRIVATE ${CERES_BENCHMARK_FAST_MATH_FLAGS})
 
+add_executable(brdf_benchmark brdf_benchmark.cc)
+add_dependencies_to_benchmark(brdf_benchmark)
+target_compile_options(brdf_benchmark PRIVATE ${CERES_BENCHMARK_FLAGS})
+
+add_executable(brdf_benchmark_fast_math brdf_benchmark.cc)
+add_dependencies_to_benchmark(brdf_benchmark_fast_math)
+target_compile_options(brdf_benchmark_fast_math PRIVATE ${CERES_BENCHMARK_FAST_MATH_FLAGS})
+
 if(CODE_GENERATION)
   target_link_libraries(snavely_reprojection_error_benchmark PUBLIC SnavelyReprojectionError)
   target_link_libraries(snavely_reprojection_error_benchmark_fast_math PUBLIC SnavelyReprojectionError)
   target_link_libraries(linear_benchmark PUBLIC Linear1CostFunction Linear10CostFunction)
   target_link_libraries(linear_benchmark_fast_math PUBLIC Linear1CostFunction Linear10CostFunction)
+  target_link_libraries(brdf_benchmark PUBLIC Brdf)
+  target_link_libraries(brdf_benchmark_fast_math PUBLIC Brdf)
 endif()
 
diff --git a/internal/ceres/autodiff_benchmarks/brdf_benchmark.cc b/internal/ceres/autodiff_benchmarks/brdf_benchmark.cc
new file mode 100644
index 0000000..b66ecf3
--- /dev/null
+++ b/internal/ceres/autodiff_benchmarks/brdf_benchmark.cc
@@ -0,0 +1,116 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2020 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: darius.rueckert@fau.de (Darius Rueckert)
+
+#include <memory>
+
+#include "benchmark/benchmark.h"
+#include "brdf_cost_function.h"
+#include "ceres/ceres.h"
+#include "codegen/test_utils.h"
+
+namespace ceres {
+
+#ifdef WITH_CODE_GENERATION
+static void BM_BrdfCodeGen(benchmark::State& state) {
+  using FunctorType = ceres::internal::CostFunctionToFunctor<Brdf>;
+
+  double material[] = {1., 2., 3., 4., 5., 6., 7., 8., 9., 10.};
+  auto c = Eigen::Vector3d(0.1, 0.2, 0.3);
+  auto n = Eigen::Vector3d(-0.1, 0.5, 0.2).normalized();
+  auto v = Eigen::Vector3d(0.5, -0.2, 0.9).normalized();
+  auto l = Eigen::Vector3d(-0.3, 0.4, -0.3).normalized();
+  auto x = Eigen::Vector3d(0.5, 0.7, -0.1).normalized();
+  auto y = Eigen::Vector3d(0.2, -0.2, -0.2).normalized();
+
+  double* parameters[7] = {
+      material, c.data(), n.data(), v.data(), l.data(), x.data(), y.data()};
+
+  double jacobian[(10 + 6 * 3) * 3];
+  double residuals[3];
+  double* jacobians[7] = {
+      jacobian + 0,
+      jacobian + 10 * 3,
+      jacobian + 13 * 3,
+      jacobian + 16 * 3,
+      jacobian + 19 * 3,
+      jacobian + 22 * 3,
+      jacobian + 25 * 3,
+  };
+
+  std::unique_ptr<ceres::CostFunction> cost_function(new Brdf());
+
+  while (state.KeepRunning()) {
+    cost_function->Evaluate(parameters, residuals, jacobians);
+  }
+}
+
+BENCHMARK(BM_BrdfCodeGen);
+#endif
+
+static void BM_BrdfAutoDiff(benchmark::State& state) {
+  using FunctorType = ceres::internal::CostFunctionToFunctor<Brdf>;
+
+  double material[] = {1., 2., 3., 4., 5., 6., 7., 8., 9., 10.};
+  auto c = Eigen::Vector3d(0.1, 0.2, 0.3);
+  auto n = Eigen::Vector3d(-0.1, 0.5, 0.2).normalized();
+  auto v = Eigen::Vector3d(0.5, -0.2, 0.9).normalized();
+  auto l = Eigen::Vector3d(-0.3, 0.4, -0.3).normalized();
+  auto x = Eigen::Vector3d(0.5, 0.7, -0.1).normalized();
+  auto y = Eigen::Vector3d(0.2, -0.2, -0.2).normalized();
+
+  double* parameters[7] = {
+      material, c.data(), n.data(), v.data(), l.data(), x.data(), y.data()};
+
+  double jacobian[(10 + 6 * 3) * 3];
+  double residuals[3];
+  double* jacobians[7] = {
+      jacobian + 0,
+      jacobian + 10 * 3,
+      jacobian + 13 * 3,
+      jacobian + 16 * 3,
+      jacobian + 19 * 3,
+      jacobian + 22 * 3,
+      jacobian + 25 * 3,
+  };
+
+  std::unique_ptr<ceres::CostFunction> cost_function(
+      new ceres::AutoDiffCostFunction<FunctorType, 3, 10, 3, 3, 3, 3, 3, 3>(
+          new FunctorType));
+
+  while (state.KeepRunning()) {
+    cost_function->Evaluate(parameters, residuals, jacobians);
+  }
+}
+
+BENCHMARK(BM_BrdfAutoDiff);
+
+}  // namespace ceres
+
+BENCHMARK_MAIN();
diff --git a/internal/ceres/autodiff_benchmarks/brdf_cost_function.h b/internal/ceres/autodiff_benchmarks/brdf_cost_function.h
new file mode 100644
index 0000000..287c409
--- /dev/null
+++ b/internal/ceres/autodiff_benchmarks/brdf_cost_function.h
@@ -0,0 +1,235 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2020 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: darius.rueckert@fau.de (Darius Rueckert)
+//
+//
+#ifndef CERES_INTERNAL_AUTODIFF_BENCHMARK_BRDF_COST_FUNCTION_H_
+#define CERES_INTERNAL_AUTODIFF_BENCHMARK_BRDF_COST_FUNCTION_H_
+
+#include <Eigen/Core>
+#include <cmath>
+
+#include "ceres/codegen/codegen_cost_function.h"
+namespace ceres {
+
+// The brdf is based on:
+// Burley, Brent, and Walt Disney Animation Studios. "Physically-based shading
+// at disney." ACM SIGGRAPH. Vol. 2012. 2012.
+//
+// The implementation is based on:
+// https://github.com/wdas/brdf/blob/master/src/brdfs/disney.brdf
+struct Brdf : public ceres::CodegenCostFunction<3, 10, 3, 3, 3, 3, 3, 3> {
+ public:
+  Brdf() {}
+
+  template <typename T>
+  bool operator()(const T* const material,
+                  const T* const c_ptr,
+                  const T* const n_ptr,
+                  const T* const v_ptr,
+                  const T* const l_ptr,
+                  const T* const x_ptr,
+                  const T* const y_ptr,
+                  T* residual) const {
+    using ceres::Ternary;
+    using Vec2 = Eigen::Matrix<T, 2, 1>;
+    using Vec3 = Eigen::Matrix<T, 3, 1>;
+
+    T metallic = material[0];
+    T subsurface = material[1];
+    T specular = material[2];
+    T roughness = material[3];
+    T specular_tint = material[4];
+    T anisotropic = material[5];
+    T sheen = material[6];
+    T sheen_tint = material[7];
+    T clearcoat = material[8];
+    T clearcoat_gloss = material[9];
+
+    Eigen::Map<const Vec3> c(c_ptr);
+    Eigen::Map<const Vec3> n(n_ptr);
+    Eigen::Map<const Vec3> v(v_ptr);
+    Eigen::Map<const Vec3> l(l_ptr);
+    Eigen::Map<const Vec3> x(x_ptr);
+    Eigen::Map<const Vec3> y(y_ptr);
+
+    const T n_dot_l = n.dot(l);
+    const T n_dot_v = n.dot(v);
+
+    const Vec3 l_p_v = l + v;
+    const Vec3 h = l_p_v / l_p_v.norm();
+
+    const T n_dot_h = n.dot(h);
+    const T l_dot_h = l.dot(h);
+
+    const T h_dot_x = h.dot(x);
+    const T h_dot_y = h.dot(y);
+
+    const T c_dlum = T(0.3) * c[0] + T(0.6) * c[1] + T(0.1) * c[2];
+
+    const Vec3 c_tint = c / c_dlum;
+
+    const Vec3 c_spec0 =
+        Lerp(specular * T(0.08) *
+                 Lerp(Vec3(T(1), T(1), T(1)), c_tint, specular_tint),
+             c,
+             metallic);
+    const Vec3 c_sheen = Lerp(Vec3(T(1), T(1), T(1)), c_tint, sheen_tint);
+
+    // Diffuse fresnel - go from 1 at normal incidence to .5 at grazing
+    // and mix in diffuse retro-reflection based on roughness
+    const T fl = SchlickFresnel(n_dot_l);
+    const T fv = SchlickFresnel(n_dot_v);
+    const T fd_90 = T(0.5) + T(2) * l_dot_h * l_dot_h * roughness;
+    const T fd = Lerp(T(1), fd_90, fl) * Lerp(T(1), fd_90, fv);
+
+    // Based on Hanrahan-Krueger brdf approximation of isotropic bssrdf
+    // 1.25 scale is used to (roughly) preserve albedo
+    // Fss90 used to "flatten" retroreflection based on roughness
+    const T fss_90 = l_dot_h * l_dot_h * roughness;
+    const T fss = Lerp(T(1), fss_90, fl) * Lerp(T(1), fss_90, fv);
+    const T ss =
+        T(1.25) * (fss * (T(1) / (n_dot_l + n_dot_v) - T(0.5)) + T(0.5));
+
+    // specular
+    const T eps = T(0.001);
+    const T aspct = Aspect(anisotropic);
+    const T ax_temp = Square(roughness) / aspct;
+    const T ay_temp = Square(roughness) * aspct;
+    const T ax = Ternary(ax_temp < eps, eps, ax_temp);
+    const T ay = Ternary(ay_temp < eps, eps, ay_temp);
+    const T ds = GTR2Aniso(n_dot_h, h_dot_x, h_dot_y, ax, ay);
+    const T fh = SchlickFresnel(l_dot_h);
+    const Vec3 fs = Lerp(c_spec0, Vec3(T(1), T(1), T(1)), fh);
+    const T roughg = Square(roughness * T(0.5) + T(0.5));
+    const T ggxn_dot_l = SmithG_GGX(n_dot_l, roughg);
+    const T ggxn_dot_v = SmithG_GGX(n_dot_v, roughg);
+    const T gs = ggxn_dot_l * ggxn_dot_v;
+
+    // sheen
+    const Vec3 f_sheen = fh * sheen * c_sheen;
+
+    // clearcoat (ior = 1.5 -> F0 = 0.04)
+    const T a = Lerp(T(0.1), T(0.001), clearcoat_gloss);
+    const T dr = GTR1(n_dot_h, a);
+    const T fr = Lerp(T(0.04), T(1), fh);
+    const T cggxn_dot_l = SmithG_GGX(n_dot_l, T(0.25));
+    const T cggxn_dot_v = SmithG_GGX(n_dot_v, T(0.25));
+    const T gr = cggxn_dot_l * cggxn_dot_v;
+
+    const Vec3 result_no_cosine =
+        (T(1.0 / M_PI) * Lerp(fd, ss, subsurface) * c + f_sheen) *
+            (T(1) - metallic) +
+        gs * fs * ds +
+        Vec3(T(0.25), T(0.25), T(0.25)) * clearcoat * gr * fr * dr;
+    const Vec3 result = n_dot_l * result_no_cosine;
+    residual[0] = result(0);
+    residual[1] = result(1);
+    residual[2] = result(2);
+
+    return true;
+  }
+
+  template <typename T>
+  T SchlickFresnel(const T& u) const {
+    T m = T(1) - u;
+    const T m2 = m * m;
+    return m2 * m2 * m;  // (1-u)^5
+  }
+
+  template <typename T>
+  T Aspect(const T& anisotropic) const {
+    return T(sqrt(T(1) - anisotropic * T(0.9)));
+  }
+
+  template <typename T>
+  T SmithG_GGX(const T& n_dot_v, const T& alpha_g) const {
+    const T a = alpha_g * alpha_g;
+    const T b = n_dot_v * n_dot_v;
+    return T(1) / (n_dot_v + T(sqrt(a + b - a * b)));
+  }
+
+  // Generalized-Trowbridge-Reitz (GTR) Microfacet Distribution
+  // See paper, Appendix B
+  template <typename T>
+  T GTR1(const T& n_dot_h, const T& a) const {
+    T result = T(0);
+
+    CERES_IF(a >= T(1)) { result = T(1 / M_PI); }
+    CERES_ELSE {
+      const T a2 = a * a;
+      const T t = T(1) + (a2 - T(1)) * n_dot_h * n_dot_h;
+      result = (a2 - T(1)) / (T(M_PI) * T(log(a2) * t));
+    }
+    CERES_ENDIF;
+    return result;
+  }
+
+  template <typename T>
+  T GTR2Aniso(const T& n_dot_h,
+              const T& h_dot_x,
+              const T& h_dot_y,
+              const T& ax,
+              const T& ay) const {
+    return T(1) / (T(M_PI) * ax * ay *
+                   Square(Square(h_dot_x / ax) + Square(h_dot_y / ay) +
+                          n_dot_h * n_dot_h));
+  }
+
+  template <typename T>
+  inline T Lerp(const T& a, const T& b, const T& u) const {
+    return a + u * (b - a);
+  }
+
+  template <typename Derived1, typename Derived2>
+  typename Derived1::PlainObject Lerp(const Eigen::MatrixBase<Derived1>& a,
+                                      const Eigen::MatrixBase<Derived2>& b,
+                                      typename Derived1::Scalar alpha) const {
+    return (typename Derived1::Scalar(1) - alpha) * a + alpha * b;
+  }
+
+  template <typename T>
+  inline T Square(const T& x) const {
+    return x * x;
+  }
+
+#ifdef WITH_CODE_GENERATION
+#include "benchmarks/brdf.h"
+#else
+  virtual bool Evaluate(double const* const* parameters,
+                        double* residuals,
+                        double** jacobians) const {
+    return false;
+  }
+#endif
+};
+
+}  // namespace ceres
+
+#endif  // CERES_INTERNAL_AUTODIFF_BENCHMARK_BRDF_COST_FUNCTION_H_