diff --git a/primitive/triangle.cpp b/primitive/triangle.cpp
index 34da50a..3632b6b 100644
--- a/primitive/triangle.cpp
+++ b/primitive/triangle.cpp
@@ -100,8 +100,13 @@ bool Triangle::intersect(Ray &ray) const {
     return false;
 
   // Calculate the normal
-  // IMPLEMENT smooth triangles, if available
-  ray.normal = normalized(crossProduct(edge1, edge2));
+  if (length(this->normal[0]) * length(this->normal[1]) * length(this->normal[2]) > EPSILON)
+    ray.normal = normalized(u * this->normal[1] + v * this->normal[2] + (1 - u - v) * this->normal[0]);
+  else
+    ray.normal = normalized(crossProduct(edge1, edge2));
+  // calculate the tangent and bitangent vectors as well
+  ray.tangent = normalized(u * this->tangent[1] + v * this->tangent[2] + (1 - u - v) * this->tangent[0]);
+  ray.bitangent = normalized(u * this->bitangent[1] + v * this->bitangent[2] + (1 - u - v) * this->bitangent[0]);
 
   // Calculate the surface position
   ray.surface = u * this->surface[1] + v * this->surface[2] + (1 - u - v) * this->surface[0];
diff --git a/shader/brdfshader.cpp b/shader/brdfshader.cpp
index ab77d37..0c40a96 100644
--- a/shader/brdfshader.cpp
+++ b/shader/brdfshader.cpp
@@ -6,8 +6,42 @@ BrdfShader::BrdfShader(char const *fileName, Color const &scale)
     : scale(scale), brdf(std::make_unique<BRDFRead>(fileName)) {}
 
 Color BrdfShader::shade(Scene const &scene, Ray const &ray) const {
-  Color illuminationColor;
+  // Calculate theta and phi
+  float thetaIn = std::acos(dotProduct(-ray.normal, ray.direction));
+  float phiIn = 0.0f;
 
-  // IMPLEMENT ME
+  // Derive local coordinate system
+  Vector3d const x = crossProduct(-ray.direction, ray.normal);
+  Vector3d const y = crossProduct(ray.normal, x);
+
+  // Accumulate the light over all light sources
+  Color illuminationColor;
+  for (const auto &light : scene.lights()) {
+    Light::Illumination illum;
+    illum = light->illuminate(scene, ray);
+
+    // Diffuse term
+    float const cosine = dotProduct(-illum.direction, ray.normal);
+    if (cosine > 0) {
+      Color color;
+
+      // Avoid numeric instability
+      if (cosine < 1) {
+        float const thetaOut = std::acos(cosine);
+
+        // Project outgoing vector into local coordinate system
+        Vector3d const c = crossProduct(-illum.direction, ray.normal);
+        float const phiOut = std::atan2(dotProduct(c, y), dotProduct(c, x));
+
+        color = Color(brdf->lookupBrdfValues(thetaIn, phiIn, thetaOut, phiOut));
+      } else {
+        color = Color(brdf->lookupBrdfValues(thetaIn, phiIn, 0, 0));
+      }
+
+      // Calculate colors
+      Color const diffuseColor = scale * color * cosine;
+      illuminationColor += diffuseColor * illum.color;
+    }
+  }
   return illuminationColor;
 }
diff --git a/shader/cooktorranceshader.cpp b/shader/cooktorranceshader.cpp
index 3f9cff5..ccfc01d 100644
--- a/shader/cooktorranceshader.cpp
+++ b/shader/cooktorranceshader.cpp
@@ -4,10 +4,51 @@
 
 CookTorranceShader::CookTorranceShader(Color const &diffCol, Color const &ctCol, float IOR, float roughness, float diffCoeff, float ctCoeff) : diffuseColor(diffCol * diffCoeff), ctColor(ctCol * ctCoeff), F0(IOR), m(roughness) {}
 
+float CookTorranceShader::D(float NdotH) const {
+  // Beckmann distribution
+  float const r2 = m * m;
+  float const NdotH2 = NdotH * NdotH;
+  return expf((NdotH2 - 1.0f) / (r2 * NdotH2)) / (4.0f * r2 * powf(NdotH, 4.0f));
+}
+
+float CookTorranceShader::F(float VdotH) const {
+  // Schlicks approximation
+  return F0 + (1.0f - F0) * powf(1.0f - VdotH, 5);
+}
+
+float CookTorranceShader::G(float NdotH, float NdotV, float VdotH, float NdotL) const { return std::min(1.0f, std::min(2.0f * NdotH * NdotV / VdotH, 2.0f * NdotH * NdotL / VdotH)); }
+
 Color CookTorranceShader::shade(Scene const &scene, Ray const &ray) const {
   Color fragmentColor;
 
-  // IMPLEMENT ME
+  if (m >= 0.0f) {
+    // Accumulate the light over all light sources
+    for (const auto &light : scene.lights()) {
+      Light::Illumination illum;
+      illum = light->illuminate(scene, ray);
+
+      float const NdotL = std::max(0.0f, dotProduct(-illum.direction, ray.normal));
+      if (NdotL <= 0.0f)
+        continue;
+
+      // Diffuse term
+      Color const diffuse = this->diffuseColor / float(PI);
+      fragmentColor += diffuse * NdotL * illum.color;
+
+      // Cook-Torrance term
+      // half angle vector
+      Vector3d const H = normalized(-illum.direction - ray.direction);
+      float const NdotH = std::max(0.0f, dotProduct(ray.normal, H));
+      float const NdotV = std::max(0.0f, dotProduct(ray.normal, -ray.direction));
+      float const VdotH = std::max(0.0f, dotProduct(-ray.direction, H));
+
+      if (NdotV * NdotL > EPSILON) {
+        Color const specular = this->ctColor * (F(VdotH) * D(NdotH) * G(NdotH, NdotV, VdotH, NdotL)) / (float(PI) * NdotV * NdotL);
+
+        fragmentColor += specular * NdotL * illum.color;
+      }
+    }
+  }
 
   return fragmentColor;
 }
\ No newline at end of file
diff --git a/shader/lambertshader.cpp b/shader/lambertshader.cpp
index 3aebacd..0e2c4a6 100644
--- a/shader/lambertshader.cpp
+++ b/shader/lambertshader.cpp
@@ -7,7 +7,13 @@ LambertShader::LambertShader(Color const &diffuseColor) : diffuseColor(diffuseCo
 Color LambertShader::shade(Scene const &scene, Ray const &ray) const {
   Color fragmentColor;
 
-  // IMPLEMENT ME
+  // Accumulate the light over all light sources
+  for (const auto &light : scene.lights()) {
+    Light::Illumination const illum = light->illuminate(scene, ray);
+    // Diffuse term
+    Color const diffuse = this->diffuseColor * std::max(dotProduct(-illum.direction, ray.normal), 0.0f);
+    fragmentColor += diffuse * illum.color;
+  }
 
   return fragmentColor;
 }
diff --git a/shader/phongshader.cpp b/shader/phongshader.cpp
index da6d3e2..ac9c3d2 100644
--- a/shader/phongshader.cpp
+++ b/shader/phongshader.cpp
@@ -10,7 +10,27 @@ PhongShader::PhongShader(Color const &diffuseColor, float diffuseCoefficient, Co
 Color PhongShader::shade(Scene const &scene, Ray const &ray) const {
   Color fragmentColor;
 
-  // IMPLEMENT ME
+  // Calculate the reflection vector
+  Vector3d const reflection = ray.direction - 2 * dotProduct(ray.normal, ray.direction) * ray.normal;
+
+  // Accumulate the light over all light sources
+  for (const auto &light : scene.lights()) {
+    Light::Illumination illum;
+    illum = light->illuminate(scene, ray);
+
+    // Diffuse term
+    Color const diffuse =
+        this->diffuseCoefficient * this->diffuseColor * std::max(dotProduct(-illum.direction, ray.normal), 0.0f);
+    fragmentColor += diffuse * illum.color;
+
+    // Specular term
+    float const cosine = dotProduct(-illum.direction, reflection);
+    if (cosine > 0) {
+      Color const specular = this->specularCoefficient * this->specularColor // highlight
+                             * powf(cosine, this->shininessExponent);        // shininess factor
+      fragmentColor += specular * illum.color;
+    }
+  }
 
   return fragmentColor;
 }