diff --git a/include/color.h b/include/color.h
index 3f1d731..5e3a958 100644
--- a/include/color.h
+++ b/include/color.h
@@ -70,6 +70,16 @@ inline ColorF32 operator*(ColorF32 c, float f)
     return result;
 }
 
+// f * c
+inline ColorF32 operator*(float f, ColorF32 c)
+{
+    ColorF32 result;
+
+    result = c * f;
+
+    return result;
+}
+
 #define COLOR_H
 #endif
 
diff --git a/include/engine.h b/include/engine.h
index bb9d214..8221e4c 100644
--- a/include/engine.h
+++ b/include/engine.h
@@ -19,7 +19,8 @@ enum Engine_Input
     ROTATE_Z_POS,
     ROTATE_Z_NEG,
     SCALE_UP,
-    SCALE_DOWN
+    SCALE_DOWN,
+    SHADING_TOGGLE
 };
 
 
diff --git a/include/geometry.h b/include/geometry.h
index cf42e63..77253bd 100644
--- a/include/geometry.h
+++ b/include/geometry.h
@@ -17,7 +17,22 @@ struct Material
 struct Face
 {
     unsigned int vertIndex[3];
-    ColorU32 color;
+    Vector normal;
+    ColorF32 color;
+};
+
+struct Vertex
+{
+    Point point;
+    Vector normal;
+    ColorF32 color;
+};
+
+struct MeshRenderData
+{
+    std::vector<Vertex> vertsTransformed;
+    std::vector<Face> culledFaces;
+    bool smooth;
 };
 
 struct Mesh
@@ -37,7 +52,7 @@ struct Mesh
 
     inline void CullBackfaces(Point camPosition)
     {
-        culledFaces.clear();
+        renderData.culledFaces.clear();
 
         for (size_t f = 0; f < faces.size(); ++f)
         {
@@ -45,41 +60,51 @@ struct Mesh
             unsigned int v1 = faces[f].vertIndex[1];
             unsigned int v2 = faces[f].vertIndex[2];
 
-            Vector v01 = vertsTransformed[v1] - vertsTransformed[v0];
-            Vector v02 = vertsTransformed[v2] - vertsTransformed[v0];
+            Vector v01 =
+                renderData.vertsTransformed[v1].point
+                - renderData.vertsTransformed[v0].point;
+
+            Vector v02 =
+                renderData.vertsTransformed[v2].point -
+                renderData.vertsTransformed[v0].point;
 
             Vector normal = Vector::Cross(v01, v02);
 
+            // Store normal for flat shading
+            faces[f].normal = normal;
+            faces[f].normal.Normalize();
+
             // Invert for Blender-compatibility
             normal = -normal;
 
             // Eye vector to viewport
-            Vector view = camPosition - vertsTransformed[v0];
+            Vector view = camPosition - renderData.vertsTransformed[v0].point;
 
             float dot = Vector::Dot(normal, view);
 
             if (dot < EPSILON_E3)
             {
-                culledFaces.push_back(faces[f]);
+                renderData.culledFaces.push_back(faces[f]);
             }
         }
     }
 
+
     Point position;
     float rotation[3];
     float scale;
-    std::vector<Point> verts;
-    std::vector<Point> vertsTransformed;
+
+    std::vector<Vertex> verts;
     std::vector<Face> faces;
-    std::vector<Face> culledFaces;
     Material material;
+
+    MeshRenderData renderData;
 };
 
 
 // PUBLIC FUNCTIONS
-void FillTriangle(
-        Engine_Buffer &buffer, ColorU32 &color,
-        Point &p0, Point &p1, Point &p2);
+void RenderMesh(
+        Engine_Buffer &buffer, MeshRenderData &mesh);
 
 
 #define GEOMETRY_H
diff --git a/include/matrix.h b/include/matrix.h
index 62b53d9..8d12db7 100644
--- a/include/matrix.h
+++ b/include/matrix.h
@@ -74,6 +74,7 @@ inline Point operator*(Point v, Matrix m)
     return result;
 }
 
+// v *=m
 inline Point &operator*=(Point &v, Matrix m)
 {
     v = v * m;
@@ -81,6 +82,26 @@ inline Point &operator*=(Point &v, Matrix m)
     return v;
 }
 
+// v * m
+inline Vector operator*(Vector v, Matrix m)
+{
+    Vector result;
+
+    for (int col = 0; col < 4; ++col)
+    {
+        float sum = 0.0;
+
+        for (int row = 0; row < 4; ++row)
+        {
+            sum += v.e[row] * m.e[row][col];
+        }
+
+        result.e[col] = sum;
+    }
+
+    return result;
+}
+
 
 #define MATRIX_H
 #endif
diff --git a/include/vec.h b/include/vec.h
index 652e8c8..b63e41d 100644
--- a/include/vec.h
+++ b/include/vec.h
@@ -76,6 +76,48 @@ inline Vector operator-(Vector v)
     return result;
 }
 
+// v1 + v2
+inline Vector operator+(Vector v1, Vector v2)
+{
+    Vector result;
+
+    result.x = v1.x + v2.x;
+    result.y = v1.y + v2.y;
+    result.z = v1.z + v2.z;
+
+    return result;
+}
+
+// v1 += v2
+inline Vector &operator+=(Vector &v1, Vector v2)
+{
+    v1 = v1 + v2;
+
+    return v1;
+}
+
+// v / f
+inline Vector operator/(Vector v, float f)
+{
+    Vector result;
+
+    float inverse = 1.0f / f;
+
+    result.x = v.x * inverse;
+    result.y = v.y * inverse;
+    result.z = v.z * inverse;
+
+    return result;
+}
+
+// v /= f
+inline Vector &operator/=(Vector &v, float f)
+{
+    v = v / f;
+
+    return v;
+}
+
 
 #define VEC_H
 #endif
diff --git a/src/engine.cpp b/src/engine.cpp
index 6281e31..2a7ff32 100644
--- a/src/engine.cpp
+++ b/src/engine.cpp
@@ -67,7 +67,8 @@ void Engine_Render(Engine_Buffer &buffer, uint32_t input)
 
     for (size_t v = 0; v < mesh.verts.size(); ++v)
     {
-        mesh.vertsTransformed[v] = mesh.verts[v] * tScale * tRotate * tTranslate;
+        mesh.renderData.vertsTransformed[v].point = mesh.verts[v].point * tScale * tRotate * tTranslate;
+        mesh.renderData.vertsTransformed[v].normal = mesh.verts[v].normal * tScale * tRotate * tTranslate;
     }
 
 
@@ -75,28 +76,43 @@ void Engine_Render(Engine_Buffer &buffer, uint32_t input)
     mesh.CullBackfaces(camera.position);
 
 
-    // Color the faces (flat shading)
-    for (size_t f = 0; f < mesh.culledFaces.size(); ++f)
+    // Color the vertices for Gouraud shading
+    if (mesh.renderData.smooth)
     {
-        unsigned int v0 = mesh.culledFaces[f].vertIndex[0];
-        unsigned int v1 = mesh.culledFaces[f].vertIndex[1];
-        unsigned int v2 = mesh.culledFaces[f].vertIndex[2];
-
-        Vector v01 = mesh.vertsTransformed[v1] - mesh.vertsTransformed[v0];
-        Vector v02 = mesh.vertsTransformed[v2] - mesh.vertsTransformed[v0];
-
-        Vector normal = Vector::Cross(v01, v02);
-        normal.Normalize();
-
-        ColorF32 totalColor = lights.ambient.ComputeColor(mesh.material.kAmbient);
-
-        for (int c = 0; c < lights.diffuseCount; ++c)
+        for (size_t f = 0; f < mesh.renderData.culledFaces.size(); ++f)
         {
-            totalColor += lights.diffuse[c].ComputeColor(
-                    mesh.material.kDiffuse, normal);
-        }
+            for (int i = 0; i < 3; ++i)
+            {
+                unsigned int v = mesh.renderData.culledFaces[f].vertIndex[i];
 
-        mesh.faces[f].color = ColorF32::ConvertToU32(totalColor);
+                ColorF32 totalColor = lights.ambient.ComputeColor(mesh.material.kAmbient);
+
+                for (int c = 0; c < lights.diffuseCount; ++c)
+                {
+                    totalColor += lights.diffuse[c].ComputeColor(
+                            mesh.material.kDiffuse, mesh.renderData.vertsTransformed[v].normal);
+                }
+
+                mesh.renderData.vertsTransformed[v].color = totalColor;
+            }
+        }
+    }
+
+    // Color the face for flat shading
+    else
+    {
+        for (size_t f = 0; f < mesh.renderData.culledFaces.size(); ++f)
+        {
+            ColorF32 totalColor = lights.ambient.ComputeColor(mesh.material.kAmbient);
+
+            for (int c = 0; c < lights.diffuseCount; ++c)
+            {
+                totalColor += lights.diffuse[c].ComputeColor(
+                        mesh.material.kDiffuse, mesh.renderData.culledFaces[f].normal);
+            }
+
+            mesh.renderData.culledFaces[f].color = totalColor;
+        }
     }
 
 
@@ -112,24 +128,11 @@ void Engine_Render(Engine_Buffer &buffer, uint32_t input)
 
     for (size_t v = 0; v < mesh.verts.size(); ++v)
     {
-        mesh.vertsTransformed[v] *= tView * tPersp * tScreen;
-        mesh.vertsTransformed[v] /= mesh.vertsTransformed[v].w;
+        mesh.renderData.vertsTransformed[v].point *= tView * tPersp * tScreen;
+        mesh.renderData.vertsTransformed[v].point /= mesh.renderData.vertsTransformed[v].point.w;
     }
 
-
-    // Fill each face with its respective color
-    for(size_t f = 0; f < mesh.culledFaces.size(); ++f)
-    {
-        unsigned int v0 = mesh.culledFaces[f].vertIndex[0];
-        unsigned int v1 = mesh.culledFaces[f].vertIndex[1];
-        unsigned int v2 = mesh.culledFaces[f].vertIndex[2];
-
-        FillTriangle(
-                buffer, mesh.faces[f].color,
-                mesh.vertsTransformed[v0],
-                mesh.vertsTransformed[v1],
-                mesh.vertsTransformed[v2]);
-    }
+    RenderMesh(buffer, mesh.renderData);
 }
 
 void Engine_Shutdown(void)
@@ -202,4 +205,13 @@ static void CheckInputs(uint32_t input)
     {
         mesh.scale -= 0.1f;
     }
+
+    if (CHECK_BIT(input, SHADING_TOGGLE))
+    {
+        mesh.renderData.smooth = true;
+    }
+    else
+    {
+        mesh.renderData.smooth = false;
+    }
 }
diff --git a/src/geometry.cpp b/src/geometry.cpp
index c3c8955..a24324d 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -29,60 +29,87 @@ struct BoundingBox
 };
 
 
-static void ComputeBarycenter(float *w, Point &v0, Point &v1, Point &v2, Point &p)
-{
-    float area = (v2.x - v0.x) * (v1.y - v0.y) - (v2.y - v0.y) * (v1.x - v0.x);
-
-    float result[3];
-    result[0] = (p.x - v1.x) * (v2.y - v1.y) - (p.y - v1.y) * (v2.x - v1.x);
-    result[1] = (p.x - v2.x) * (v0.y - v2.y) - (p.y - v2.y) * (v0.x - v2.x);
-    result[2] = (p.x - v0.x) * (v1.y - v0.y) - (p.y - v0.y) * (v1.x - v0.x);
-
-    result[0] /= area;
-    result[1] /= area;
-    result[2] /= area;
-
-    w[0] = result[0];
-    w[1] = result[1];
-    w[2] = result[2];
-}
-
-
 // PUBLIC FUNCTIONS
-void FillTriangle(
-        Engine_Buffer &buffer, ColorU32 &color,
-        Point &v0, Point &v1, Point &v2)
+void RenderMesh(Engine_Buffer &buffer, MeshRenderData &mesh)
 {
-    BoundingBox box(v0, v1, v2);
-    int yMin = (int)ceilf(box.yMin);
-    int yMax = (int)ceilf(box.yMax) - 1;
-    int xMin = (int)ceilf(box.xMin);
-    int xMax = (int)ceilf(box.xMax) - 1;
-
-    for (int y = yMin; y <= yMax; ++y)
+    for(size_t f = 0; f < mesh.culledFaces.size(); ++f)
     {
-        for (int x = xMin; x <= xMax; ++x)
+        unsigned int vIndex0 = mesh.culledFaces[f].vertIndex[0];
+        unsigned int vIndex1 = mesh.culledFaces[f].vertIndex[1];
+        unsigned int vIndex2 = mesh.culledFaces[f].vertIndex[2];
+
+        Vertex v0 = mesh.vertsTransformed[vIndex0];
+        Vertex v1 = mesh.vertsTransformed[vIndex1];
+        Vertex v2 = mesh.vertsTransformed[vIndex2];
+
+        // Bounding box used to for iterating over possible pixels of this triangle
+        BoundingBox box(v0.point, v1.point, v2.point);
+        int yMin = (int)ceilf(box.yMin);
+        int yMax = (int)ceilf(box.yMax) - 1;
+        int xMin = (int)ceilf(box.xMin);
+        int xMax = (int)ceilf(box.xMax) - 1;
+
+        // Constants for this triangle used for barycentric calculations
+        Vector v01 = v1.point - v0.point;
+        Vector v02 = v2.point - v0.point;
+        float dot0101 = Vector::Dot(v01, v01);
+        float dot0102 = Vector::Dot(v01, v02);
+        float dot0202 = Vector::Dot(v02, v02);
+
+        // Iterate bounding box and determine if each point is in the triangle
+        for (int y = yMin; y <= yMax; ++y)
         {
-            Point p(x, y, 1.0f);
-
-            float barycenter[3];
-
-            ComputeBarycenter(barycenter, v0, v1, v2, p);
-
-            if (barycenter[0] >= 0 && barycenter[1] >= 0 && barycenter[2] >= 0)
+            for (int x = xMin; x <= xMax; ++x)
             {
-                float zInv =
-                    (barycenter[0] * (1.0f/v0.z))
-                    + (barycenter[1] * (1.0f/v1.z))
-                    + (barycenter[2] * (1.0f/v2.z));
+                Point p(x, y, 1.0f);
 
-                int pixel = (y * buffer.width + x);
+                Vector v0P = p - v0.point;
+                float dot0P01 = Vector::Dot(v0P, v01);
+                float dot0P02 = Vector::Dot(v0P, v02);
+                float denomInv = 1.0f / ((dot0101 * dot0202) - (dot0102 * dot0102));
 
-                if (zInv > buffer.zbuffer[pixel])
+                float barycenter[3];
+                barycenter[1] = (dot0202 * dot0P01 - dot0102 * dot0P02) * denomInv;
+                barycenter[2] = (dot0101 * dot0P02 - dot0102 * dot0P01) * denomInv;
+                barycenter[0] = 1.0f - barycenter[1] - barycenter[2];
+
+                // Point is inside the triangle
+                if ( (barycenter[0] >= 0.0f)
+                        && (barycenter[1] >= 0.0f)
+                        && (barycenter[2] >= 0.0f))
                 {
-                    DrawPixel(buffer.buffer, buffer.width, color.u32, x, y);
+                    ColorF32 totalColor;
 
-                    buffer.zbuffer[pixel] = zInv;
+                    if (mesh.smooth)
+                    {
+                        totalColor =
+                            (barycenter[0] * v0.color)
+                            + (barycenter[1] * v1.color)
+                            + (barycenter[2] * v2.color);
+                    }
+                    else
+                    {
+                        totalColor = mesh.culledFaces[f].color;
+                    }
+
+                    ColorU32 color = ColorF32::ConvertToU32(totalColor);
+
+
+                    float z =
+                        (barycenter[0] * v0.point.z)
+                        + (barycenter[1] * v1.point.z)
+                        + (barycenter[2] * v2.point.z);
+
+                    float zInv = 1.0f / z;
+
+                    int pixel = (y * buffer.width + x);
+
+                    if (zInv > buffer.zbuffer[pixel])
+                    {
+                        DrawPixel(buffer.buffer, buffer.width, color.u32, x, y);
+
+                        buffer.zbuffer[pixel] = zInv;
+                    }
                 }
             }
         }
diff --git a/src/loader.cpp b/src/loader.cpp
index 0184086..4e249a9 100644
--- a/src/loader.cpp
+++ b/src/loader.cpp
@@ -6,6 +6,7 @@
 
 // STATIC PROTOTYPES
 static char *GetLine(char *buffer, int maxLength, FILE *fp);
+static void ComputeNormals(Mesh &mesh);
 
 
 // PUBLIC FUNCTIONS
@@ -31,14 +32,14 @@ int LoadMesh(char *filename, Mesh &mesh)
     {
         if (token[0] == 'v')
         {
-            Point v;
+            Vertex v;
 
             sscanf(
                     token, "%s %f %f %f",
                     garbage,
-                    &v.x,
-                    &v.y,
-                    &v.z);
+                    &v.point.x,
+                    &v.point.y,
+                    &v.point.z);
 
             mesh.verts.push_back(v);
         }
@@ -63,11 +64,13 @@ int LoadMesh(char *filename, Mesh &mesh)
         token = GetLine(buffer, sizeof(buffer), fp);
     }
 
+    ComputeNormals(mesh);
+
     printf("OBJ: %s\n", filename);
     printf("Verts: %lu\n", mesh.verts.size());
     printf("Faces: %lu\n", mesh.faces.size());
 
-    mesh.vertsTransformed.resize(mesh.verts.size());
+    mesh.renderData.vertsTransformed.resize(mesh.verts.size());
 
     fclose(fp);
 
@@ -93,3 +96,40 @@ static char *GetLine(char *buffer, int maxLength, FILE *fp)
         return buffer;
     }
 }
+
+
+static void ComputeNormals(Mesh &mesh)
+{
+    int *vertexNormalCount = (int*)calloc((size_t)(mesh.verts.size()), sizeof(int));
+
+    for (size_t f = 0; f < mesh.faces.size(); ++f)
+    {
+        unsigned int v0 = mesh.faces[f].vertIndex[0];
+        unsigned int v1 = mesh.faces[f].vertIndex[1];
+        unsigned int v2 = mesh.faces[f].vertIndex[2];
+
+        Vector v01 = mesh.verts[v1].point - mesh.verts[v0].point;
+        Vector v02 = mesh.verts[v2].point - mesh.verts[v0].point;
+
+        Vector normal = Vector::Cross(v01, v02);
+
+        // Add each vertex's normal to the sum for future averaging
+        mesh.verts[v0].normal += normal;
+        mesh.verts[v1].normal += normal;
+        mesh.verts[v2].normal += normal;
+
+        ++vertexNormalCount[v0];
+        ++vertexNormalCount[v1];
+        ++vertexNormalCount[v2];
+    }
+
+    for (size_t v = 0; v < mesh.verts.size(); ++v)
+    {
+        if (vertexNormalCount[v] > 0)
+        {
+            // Compute the average normal for this vertex
+            mesh.verts[v].normal /= vertexNormalCount[v];
+            mesh.verts[v].normal.Normalize();
+        }
+    }
+}
diff --git a/src/platform.cpp b/src/platform.cpp
index bd1bf8f..61339a0 100644
--- a/src/platform.cpp
+++ b/src/platform.cpp
@@ -169,6 +169,10 @@ static void HandleEvent(
                 {
                     SET_BIT(platform.input, ROTATE_Y_NEG);
                 } break;
+                case SDLK_g:
+                {
+                    SET_BIT(platform.input, SHADING_TOGGLE);
+                } break;
                 case SDLK_UP:
                 {
                     SET_BIT(platform.input, SCALE_UP);
@@ -231,6 +235,10 @@ static void HandleEvent(
                 {
                     CLEAR_BIT(platform.input, ROTATE_Y_NEG);
                 } break;
+                case SDLK_g:
+                {
+                    CLEAR_BIT(platform.input, SHADING_TOGGLE);
+                } break;
                 case SDLK_UP:
                 {
                     CLEAR_BIT(platform.input, SCALE_UP);