diff --git a/3rdparty/OBJ_Loader/CMakeLists.txt b/3rdparty/OBJ_Loader/CMakeLists.txt
new file mode 100644
index 0000000..79150f5
--- /dev/null
+++ b/3rdparty/OBJ_Loader/CMakeLists.txt
@@ -0,0 +1,7 @@
+cmake_minimum_required(VERSION 3.14)
+project(OBJ_Loader)
+
+set(CMAKE_CXX_STANDARD 11)
+
+add_library(${PROJECT_NAME} INTERFACE)
+target_include_directories(${PROJECT_NAME} INTERFACE ${CMAKE_CURRENT_SOURCE_DIR}/include/)
diff --git a/3rdparty/OBJ_Loader/include/OBJ_Loader.h b/3rdparty/OBJ_Loader/include/OBJ_Loader.h
new file mode 100644
index 0000000..e64ae5f
--- /dev/null
+++ b/3rdparty/OBJ_Loader/include/OBJ_Loader.h
@@ -0,0 +1,1015 @@
+// OBJ_Loader.h - A Single Header OBJ Model Loader
+
+#pragma once
+
+// Iostream - STD I/O Library
+#include <iostream>
+
+// Vector - STD Vector/Array Library
+#include <vector>
+
+// String - STD String Library
+#include <string>
+
+// fStream - STD File I/O Library
+#include <fstream>
+
+// Math.h - STD math Library
+#include <math.h>
+
+// Print progress to console while loading (large models)
+#define OBJL_CONSOLE_OUTPUT
+
+// Namespace: OBJL
+//
+// Description: The namespace that holds eveyrthing that
+//	is needed and used for the OBJ Model Loader
+namespace objl {
+// Structure: Vector2
+//
+// Description: A 2D Vector that Holds Positional Data
+struct Vector2 {
+  // Default Constructor
+  Vector2() {
+    X = 0.0f;
+    Y = 0.0f;
+  }
+  // Variable Set Constructor
+  Vector2(float X_, float Y_) {
+    X = X_;
+    Y = Y_;
+  }
+  // Bool Equals Operator Overload
+  bool operator==(const Vector2 &other) const {
+    return (this->X == other.X && this->Y == other.Y);
+  }
+  // Bool Not Equals Operator Overload
+  bool operator!=(const Vector2 &other) const {
+    return !(this->X == other.X && this->Y == other.Y);
+  }
+  // Addition Operator Overload
+  Vector2 operator+(const Vector2 &right) const {
+    return Vector2(this->X + right.X, this->Y + right.Y);
+  }
+  // Subtraction Operator Overload
+  Vector2 operator-(const Vector2 &right) const {
+    return Vector2(this->X - right.X, this->Y - right.Y);
+  }
+  // Float Multiplication Operator Overload
+  Vector2 operator*(const float &other) const {
+    return Vector2(this->X * other, this->Y * other);
+  }
+
+  // Positional Variables
+  float X;
+  float Y;
+};
+
+// Structure: Vector3
+//
+// Description: A 3D Vector that Holds Positional Data
+struct Vector3 {
+  // Default Constructor
+  Vector3() {
+    X = 0.0f;
+    Y = 0.0f;
+    Z = 0.0f;
+  }
+  // Variable Set Constructor
+  Vector3(float X_, float Y_, float Z_) {
+    X = X_;
+    Y = Y_;
+    Z = Z_;
+  }
+  // Bool Equals Operator Overload
+  bool operator==(const Vector3 &other) const {
+    return (this->X == other.X && this->Y == other.Y && this->Z == other.Z);
+  }
+  // Bool Not Equals Operator Overload
+  bool operator!=(const Vector3 &other) const {
+    return !(this->X == other.X && this->Y == other.Y && this->Z == other.Z);
+  }
+  // Addition Operator Overload
+  Vector3 operator+(const Vector3 &right) const {
+    return Vector3(this->X + right.X, this->Y + right.Y, this->Z + right.Z);
+  }
+  // Subtraction Operator Overload
+  Vector3 operator-(const Vector3 &right) const {
+    return Vector3(this->X - right.X, this->Y - right.Y, this->Z - right.Z);
+  }
+  // Float Multiplication Operator Overload
+  Vector3 operator*(const float &other) const {
+    return Vector3(this->X * other, this->Y * other, this->Z * other);
+  }
+  // Float Division Operator Overload
+  Vector3 operator/(const float &other) const {
+    return Vector3(this->X / other, this->Y / other, this->Z / other);
+  }
+
+  // Positional Variables
+  float X;
+  float Y;
+  float Z;
+};
+
+// Structure: Vertex
+//
+// Description: Model Vertex object that holds
+//	a Position, Normal, and Texture Coordinate
+struct Vertex {
+  // Position Vector
+  Vector3 Position;
+
+  // Normal Vector
+  Vector3 Normal;
+
+  // Texture Coordinate Vector
+  Vector2 TextureCoordinate;
+};
+
+struct Material {
+  Material() {
+    Ns = 0.0f;
+    Ni = 0.0f;
+    d = 0.0f;
+    illum = 0;
+  }
+
+  // Material Name
+  std::string name;
+  // Ambient Color
+  Vector3 Ka;
+  // Diffuse Color
+  Vector3 Kd;
+  // Specular Color
+  Vector3 Ks;
+  // Specular Exponent
+  float Ns;
+  // Optical Density
+  float Ni;
+  // Dissolve
+  float d;
+  // Illumination
+  int illum;
+  // Ambient Texture Map
+  std::string map_Ka;
+  // Diffuse Texture Map
+  std::string map_Kd;
+  // Specular Texture Map
+  std::string map_Ks;
+  // Specular Hightlight Map
+  std::string map_Ns;
+  // Alpha Texture Map
+  std::string map_d;
+  // Bump Map
+  std::string map_bump;
+};
+
+// Structure: Mesh
+//
+// Description: A Simple Mesh Object that holds
+//	a name, a vertex list, and an index list
+struct Mesh {
+  // Default Constructor
+  Mesh() {}
+  // Variable Set Constructor
+  Mesh(std::vector<Vertex> &_Vertices, std::vector<unsigned int> &_Indices) {
+    Vertices = _Vertices;
+    Indices = _Indices;
+  }
+  // Mesh Name
+  std::string MeshName;
+  // Vertex List
+  std::vector<Vertex> Vertices;
+  // Index List
+  std::vector<unsigned int> Indices;
+
+  // Material
+  Material MeshMaterial;
+};
+
+// Namespace: Math
+//
+// Description: The namespace that holds all of the math
+//	functions need for OBJL
+namespace math {
+// Vector3 Cross Product
+Vector3 CrossV3(const Vector3 a, const Vector3 b) {
+  return Vector3(a.Y * b.Z - a.Z * b.Y, a.Z * b.X - a.X * b.Z,
+                 a.X * b.Y - a.Y * b.X);
+}
+
+// Vector3 Magnitude Calculation
+float MagnitudeV3(const Vector3 in) {
+  return (sqrtf(powf(in.X, 2) + powf(in.Y, 2) + powf(in.Z, 2)));
+}
+
+// Vector3 DotProduct
+float DotV3(const Vector3 a, const Vector3 b) {
+  return (a.X * b.X) + (a.Y * b.Y) + (a.Z * b.Z);
+}
+
+// Angle between 2 Vector3 Objects
+float AngleBetweenV3(const Vector3 a, const Vector3 b) {
+  float angle = DotV3(a, b);
+  angle /= (MagnitudeV3(a) * MagnitudeV3(b));
+  return angle = acosf(angle);
+}
+
+// Projection Calculation of a onto b
+Vector3 ProjV3(const Vector3 a, const Vector3 b) {
+  Vector3 bn = b / MagnitudeV3(b);
+  return bn * DotV3(a, bn);
+}
+} // namespace math
+
+// Namespace: Algorithm
+//
+// Description: The namespace that holds all of the
+// Algorithms needed for OBJL
+namespace algorithm {
+// Vector3 Multiplication Opertor Overload
+Vector3 operator*(const float &left, const Vector3 &right) {
+  return Vector3(right.X * left, right.Y * left, right.Z * left);
+}
+
+// A test to see if P1 is on the same side as P2 of a line segment ab
+bool SameSide(Vector3 p1, Vector3 p2, Vector3 a, Vector3 b) {
+  Vector3 cp1 = math::CrossV3(b - a, p1 - a);
+  Vector3 cp2 = math::CrossV3(b - a, p2 - a);
+
+  if (math::DotV3(cp1, cp2) >= 0)
+    return true;
+  else
+    return false;
+}
+
+// Generate a cross produect normal for a triangle
+Vector3 GenTriNormal(Vector3 t1, Vector3 t2, Vector3 t3) {
+  Vector3 u = t2 - t1;
+  Vector3 v = t3 - t1;
+
+  Vector3 normal = math::CrossV3(u, v);
+
+  return normal;
+}
+
+// Check to see if a Vector3 Point is within a 3 Vector3 Triangle
+bool inTriangle(Vector3 point, Vector3 tri1, Vector3 tri2, Vector3 tri3) {
+  // Test to see if it is within an infinite prism that the triangle outlines.
+  bool within_tri_prisim = SameSide(point, tri1, tri2, tri3) &&
+                           SameSide(point, tri2, tri1, tri3) &&
+                           SameSide(point, tri3, tri1, tri2);
+
+  // If it isn't it will never be on the triangle
+  if (!within_tri_prisim)
+    return false;
+
+  // Calulate Triangle's Normal
+  Vector3 n = GenTriNormal(tri1, tri2, tri3);
+
+  // Project the point onto this normal
+  Vector3 proj = math::ProjV3(point, n);
+
+  // If the distance from the triangle to the point is 0
+  //	it lies on the triangle
+  if (math::MagnitudeV3(proj) == 0)
+    return true;
+  else
+    return false;
+}
+
+// Split a String into a string array at a given token
+inline void split(const std::string &in, std::vector<std::string> &out,
+                  std::string token) {
+  out.clear();
+
+  std::string temp;
+
+  for (int i = 0; i < int(in.size()); i++) {
+    std::string test = in.substr(i, token.size());
+
+    if (test == token) {
+      if (!temp.empty()) {
+        out.push_back(temp);
+        temp.clear();
+        i += (int)token.size() - 1;
+      } else {
+        out.push_back("");
+      }
+    } else if (i + token.size() >= in.size()) {
+      temp += in.substr(i, token.size());
+      out.push_back(temp);
+      break;
+    } else {
+      temp += in[i];
+    }
+  }
+}
+
+// Get tail of string after first token and possibly following spaces
+inline std::string tail(const std::string &in) {
+  size_t token_start = in.find_first_not_of(" \t");
+  size_t space_start = in.find_first_of(" \t", token_start);
+  size_t tail_start = in.find_first_not_of(" \t", space_start);
+  size_t tail_end = in.find_last_not_of(" \t");
+  if (tail_start != std::string::npos && tail_end != std::string::npos) {
+    return in.substr(tail_start, tail_end - tail_start + 1);
+  } else if (tail_start != std::string::npos) {
+    return in.substr(tail_start);
+  }
+  return "";
+}
+
+// Get first token of string
+inline std::string firstToken(const std::string &in) {
+  if (!in.empty()) {
+    size_t token_start = in.find_first_not_of(" \t");
+    size_t token_end = in.find_first_of(" \t", token_start);
+    if (token_start != std::string::npos && token_end != std::string::npos) {
+      return in.substr(token_start, token_end - token_start);
+    } else if (token_start != std::string::npos) {
+      return in.substr(token_start);
+    }
+  }
+  return "";
+}
+
+// Get element at given index position
+template <class T>
+inline const T &getElement(const std::vector<T> &elements, std::string &index) {
+  int idx = std::stoi(index);
+  if (idx < 0)
+    idx = int(elements.size()) + idx;
+  else
+    idx--;
+  return elements[idx];
+}
+} // namespace algorithm
+
+// Class: Loader
+//
+// Description: The OBJ Model Loader
+class Loader {
+public:
+  // Default Constructor
+  Loader() {}
+  ~Loader() { LoadedMeshes.clear(); }
+
+  // Load a file into the loader
+  //
+  // If file is loaded return true
+  //
+  // If the file is unable to be found
+  // or unable to be loaded return false
+  bool LoadFile(std::string Path) {
+    // If the file is not an .obj file return false
+    if (Path.substr(Path.size() - 4, 4) != ".obj")
+      return false;
+
+    std::ifstream file(Path);
+
+    if (!file.is_open())
+      return false;
+
+    LoadedMeshes.clear();
+    LoadedVertices.clear();
+    LoadedIndices.clear();
+
+    std::vector<Vector3> Positions;
+    std::vector<Vector2> TCoords;
+    std::vector<Vector3> Normals;
+
+    std::vector<Vertex> Vertices;
+    std::vector<unsigned int> Indices;
+
+    std::vector<std::string> MeshMatNames;
+
+    bool listening = false;
+    std::string meshname;
+
+    Mesh tempMesh;
+
+#ifdef OBJL_CONSOLE_OUTPUT
+    const unsigned int outputEveryNth = 1000;
+    unsigned int outputIndicator = outputEveryNth;
+#endif
+
+    std::string curline;
+    while (std::getline(file, curline)) {
+#ifdef OBJL_CONSOLE_OUTPUT
+      if ((outputIndicator = ((outputIndicator + 1) % outputEveryNth)) == 1) {
+        if (!meshname.empty()) {
+          std::cout << "\r- " << meshname << "\t| vertices > "
+                    << Positions.size() << "\t| texcoords > " << TCoords.size()
+                    << "\t| normals > " << Normals.size() << "\t| triangles > "
+                    << (Vertices.size() / 3)
+                    << (!MeshMatNames.empty()
+                            ? "\t| material: " + MeshMatNames.back()
+                            : "");
+        }
+      }
+#endif
+
+      // Generate a Mesh Object or Prepare for an object to be created
+      if (algorithm::firstToken(curline) == "o" ||
+          algorithm::firstToken(curline) == "g" || curline[0] == 'g') {
+        if (!listening) {
+          listening = true;
+
+          if (algorithm::firstToken(curline) == "o" ||
+              algorithm::firstToken(curline) == "g") {
+            meshname = algorithm::tail(curline);
+          } else {
+            meshname = "unnamed";
+          }
+        } else {
+          // Generate the mesh to put into the array
+
+          if (!Indices.empty() && !Vertices.empty()) {
+            // Create Mesh
+            tempMesh = Mesh(Vertices, Indices);
+            tempMesh.MeshName = meshname;
+
+            // Insert Mesh
+            LoadedMeshes.push_back(tempMesh);
+
+            // Cleanup
+            Vertices.clear();
+            Indices.clear();
+            meshname.clear();
+
+            meshname = algorithm::tail(curline);
+          } else {
+            if (algorithm::firstToken(curline) == "o" ||
+                algorithm::firstToken(curline) == "g") {
+              meshname = algorithm::tail(curline);
+            } else {
+              meshname = "unnamed";
+            }
+          }
+        }
+#ifdef OBJL_CONSOLE_OUTPUT
+        std::cout << std::endl;
+        outputIndicator = 0;
+#endif
+      }
+      // Generate a Vertex Position
+      if (algorithm::firstToken(curline) == "v") {
+        std::vector<std::string> spos;
+        Vector3 vpos;
+        algorithm::split(algorithm::tail(curline), spos, " ");
+
+        vpos.X = std::stof(spos[0]);
+        vpos.Y = std::stof(spos[1]);
+        vpos.Z = std::stof(spos[2]);
+
+        Positions.push_back(vpos);
+      }
+      // Generate a Vertex Texture Coordinate
+      if (algorithm::firstToken(curline) == "vt") {
+        std::vector<std::string> stex;
+        Vector2 vtex;
+        algorithm::split(algorithm::tail(curline), stex, " ");
+
+        vtex.X = std::stof(stex[0]);
+        vtex.Y = std::stof(stex[1]);
+
+        TCoords.push_back(vtex);
+      }
+      // Generate a Vertex Normal;
+      if (algorithm::firstToken(curline) == "vn") {
+        std::vector<std::string> snor;
+        Vector3 vnor;
+        algorithm::split(algorithm::tail(curline), snor, " ");
+
+        vnor.X = std::stof(snor[0]);
+        vnor.Y = std::stof(snor[1]);
+        vnor.Z = std::stof(snor[2]);
+
+        Normals.push_back(vnor);
+      }
+      // Generate a Face (vertices & indices)
+      if (algorithm::firstToken(curline) == "f") {
+        // Generate the vertices
+        std::vector<Vertex> vVerts;
+        GenVerticesFromRawOBJ(vVerts, Positions, TCoords, Normals, curline);
+
+        // Add Vertices
+        for (int i = 0; i < int(vVerts.size()); i++) {
+          Vertices.push_back(vVerts[i]);
+
+          LoadedVertices.push_back(vVerts[i]);
+        }
+
+        std::vector<unsigned int> iIndices;
+
+        VertexTriangluation(iIndices, vVerts);
+
+        // Add Indices
+        for (int i = 0; i < int(iIndices.size()); i++) {
+          unsigned int indnum =
+              (unsigned int)((Vertices.size()) - vVerts.size()) + iIndices[i];
+          Indices.push_back(indnum);
+
+          indnum = (unsigned int)((LoadedVertices.size()) - vVerts.size()) +
+                   iIndices[i];
+          LoadedIndices.push_back(indnum);
+        }
+      }
+      // Get Mesh Material Name
+      if (algorithm::firstToken(curline) == "usemtl") {
+        MeshMatNames.push_back(algorithm::tail(curline));
+
+        // Create new Mesh, if Material changes within a group
+        if (!Indices.empty() && !Vertices.empty()) {
+          // Create Mesh
+          tempMesh = Mesh(Vertices, Indices);
+          tempMesh.MeshName = meshname;
+          int i = 2;
+          while (1) {
+            tempMesh.MeshName = meshname + "_" + std::to_string(i);
+
+            for (auto &m : LoadedMeshes)
+              if (m.MeshName == tempMesh.MeshName)
+                continue;
+            break;
+          }
+
+          // Insert Mesh
+          LoadedMeshes.push_back(tempMesh);
+
+          // Cleanup
+          Vertices.clear();
+          Indices.clear();
+        }
+
+#ifdef OBJL_CONSOLE_OUTPUT
+        outputIndicator = 0;
+#endif
+      }
+      // Load Materials
+      if (algorithm::firstToken(curline) == "mtllib") {
+        // Generate LoadedMaterial
+
+        // Generate a path to the material file
+        std::vector<std::string> temp;
+        algorithm::split(Path, temp, "/");
+
+        std::string pathtomat = "";
+
+        if (temp.size() != 1) {
+          for (int i = 0; i < temp.size() - 1; i++) {
+            pathtomat += temp[i] + "/";
+          }
+        }
+
+        pathtomat += algorithm::tail(curline);
+
+#ifdef OBJL_CONSOLE_OUTPUT
+        std::cout << std::endl
+                  << "- find materials in: " << pathtomat << std::endl;
+#endif
+
+        // Load Materials
+        LoadMaterials(pathtomat);
+      }
+    }
+
+#ifdef OBJL_CONSOLE_OUTPUT
+    std::cout << std::endl;
+#endif
+
+    // Deal with last mesh
+
+    if (!Indices.empty() && !Vertices.empty()) {
+      // Create Mesh
+      tempMesh = Mesh(Vertices, Indices);
+      tempMesh.MeshName = meshname;
+
+      // Insert Mesh
+      LoadedMeshes.push_back(tempMesh);
+    }
+
+    file.close();
+
+    // Set Materials for each Mesh
+    for (int i = 0; i < MeshMatNames.size(); i++) {
+      std::string matname = MeshMatNames[i];
+
+      // Find corresponding material name in loaded materials
+      // when found copy material variables into mesh material
+      for (int j = 0; j < LoadedMaterials.size(); j++) {
+        if (LoadedMaterials[j].name == matname) {
+          LoadedMeshes[i].MeshMaterial = LoadedMaterials[j];
+          break;
+        }
+      }
+    }
+
+    if (LoadedMeshes.empty() && LoadedVertices.empty() &&
+        LoadedIndices.empty()) {
+      return false;
+    } else {
+      return true;
+    }
+  }
+
+  // Loaded Mesh Objects
+  std::vector<Mesh> LoadedMeshes;
+  // Loaded Vertex Objects
+  std::vector<Vertex> LoadedVertices;
+  // Loaded Index Positions
+  std::vector<unsigned int> LoadedIndices;
+  // Loaded Material Objects
+  std::vector<Material> LoadedMaterials;
+
+private:
+  // Generate vertices from a list of positions,
+  //	tcoords, normals and a face line
+  void GenVerticesFromRawOBJ(std::vector<Vertex> &oVerts,
+                             const std::vector<Vector3> &iPositions,
+                             const std::vector<Vector2> &iTCoords,
+                             const std::vector<Vector3> &iNormals,
+                             std::string icurline) {
+    std::vector<std::string> sface, svert;
+    Vertex vVert;
+    algorithm::split(algorithm::tail(icurline), sface, " ");
+
+    bool noNormal = false;
+
+    // For every given vertex do this
+    for (int i = 0; i < int(sface.size()); i++) {
+      // See What type the vertex is.
+      int vtype;
+
+      algorithm::split(sface[i], svert, "/");
+
+      // Check for just position - v1
+      if (svert.size() == 1) {
+        // Only position
+        vtype = 1;
+      }
+
+      // Check for position & texture - v1/vt1
+      if (svert.size() == 2) {
+        // Position & Texture
+        vtype = 2;
+      }
+
+      // Check for Position, Texture and Normal - v1/vt1/vn1
+      // or if Position and Normal - v1//vn1
+      if (svert.size() == 3) {
+        if (svert[1] != "") {
+          // Position, Texture, and Normal
+          vtype = 4;
+        } else {
+          // Position & Normal
+          vtype = 3;
+        }
+      }
+
+      // Calculate and store the vertex
+      switch (vtype) {
+      case 1: // P
+      {
+        vVert.Position = algorithm::getElement(iPositions, svert[0]);
+        vVert.TextureCoordinate = Vector2(0, 0);
+        noNormal = true;
+        oVerts.push_back(vVert);
+        break;
+      }
+      case 2: // P/T
+      {
+        vVert.Position = algorithm::getElement(iPositions, svert[0]);
+        vVert.TextureCoordinate = algorithm::getElement(iTCoords, svert[1]);
+        noNormal = true;
+        oVerts.push_back(vVert);
+        break;
+      }
+      case 3: // P//N
+      {
+        vVert.Position = algorithm::getElement(iPositions, svert[0]);
+        vVert.TextureCoordinate = Vector2(0, 0);
+        vVert.Normal = algorithm::getElement(iNormals, svert[2]);
+        oVerts.push_back(vVert);
+        break;
+      }
+      case 4: // P/T/N
+      {
+        vVert.Position = algorithm::getElement(iPositions, svert[0]);
+        vVert.TextureCoordinate = algorithm::getElement(iTCoords, svert[1]);
+        vVert.Normal = algorithm::getElement(iNormals, svert[2]);
+        oVerts.push_back(vVert);
+        break;
+      }
+      default: {
+        break;
+      }
+      }
+    }
+
+    // take care of missing normals
+    // these may not be truly acurate but it is the
+    // best they get for not compiling a mesh with normals
+    if (noNormal) {
+      Vector3 A = oVerts[0].Position - oVerts[1].Position;
+      Vector3 B = oVerts[2].Position - oVerts[1].Position;
+
+      Vector3 normal = math::CrossV3(A, B);
+
+      for (int i = 0; i < int(oVerts.size()); i++) {
+        oVerts[i].Normal = normal;
+      }
+    }
+  }
+
+  // Triangulate a list of vertices into a face by printing
+  //	inducies corresponding with triangles within it
+  void VertexTriangluation(std::vector<unsigned int> &oIndices,
+                           const std::vector<Vertex> &iVerts) {
+    // If there are 2 or less verts,
+    // no triangle can be created,
+    // so exit
+    if (iVerts.size() < 3) {
+      return;
+    }
+    // If it is a triangle no need to calculate it
+    if (iVerts.size() == 3) {
+      oIndices.push_back(0);
+      oIndices.push_back(1);
+      oIndices.push_back(2);
+      return;
+    }
+
+    // Create a list of vertices
+    std::vector<Vertex> tVerts = iVerts;
+
+    while (true) {
+      // For every vertex
+      for (int i = 0; i < int(tVerts.size()); i++) {
+        // pPrev = the previous vertex in the list
+        Vertex pPrev;
+        if (i == 0) {
+          pPrev = tVerts[tVerts.size() - 1];
+        } else {
+          pPrev = tVerts[i - 1];
+        }
+
+        // pCur = the current vertex;
+        Vertex pCur = tVerts[i];
+
+        // pNext = the next vertex in the list
+        Vertex pNext;
+        if (i == tVerts.size() - 1) {
+          pNext = tVerts[0];
+        } else {
+          pNext = tVerts[i + 1];
+        }
+
+        // Check to see if there are only 3 verts left
+        // if so this is the last triangle
+        if (tVerts.size() == 3) {
+          // Create a triangle from pCur, pPrev, pNext
+          for (int j = 0; j < int(tVerts.size()); j++) {
+            if (iVerts[j].Position == pCur.Position)
+              oIndices.push_back(j);
+            if (iVerts[j].Position == pPrev.Position)
+              oIndices.push_back(j);
+            if (iVerts[j].Position == pNext.Position)
+              oIndices.push_back(j);
+          }
+
+          tVerts.clear();
+          break;
+        }
+        if (tVerts.size() == 4) {
+          // Create a triangle from pCur, pPrev, pNext
+          for (int j = 0; j < int(iVerts.size()); j++) {
+            if (iVerts[j].Position == pCur.Position)
+              oIndices.push_back(j);
+            if (iVerts[j].Position == pPrev.Position)
+              oIndices.push_back(j);
+            if (iVerts[j].Position == pNext.Position)
+              oIndices.push_back(j);
+          }
+
+          Vector3 tempVec;
+          for (int j = 0; j < int(tVerts.size()); j++) {
+            if (tVerts[j].Position != pCur.Position &&
+                tVerts[j].Position != pPrev.Position &&
+                tVerts[j].Position != pNext.Position) {
+              tempVec = tVerts[j].Position;
+              break;
+            }
+          }
+
+          // Create a triangle from pCur, pPrev, pNext
+          for (int j = 0; j < int(iVerts.size()); j++) {
+            if (iVerts[j].Position == pPrev.Position)
+              oIndices.push_back(j);
+            if (iVerts[j].Position == pNext.Position)
+              oIndices.push_back(j);
+            if (iVerts[j].Position == tempVec)
+              oIndices.push_back(j);
+          }
+
+          tVerts.clear();
+          break;
+        }
+
+        // If Vertex is not an interior vertex
+        float angle = math::AngleBetweenV3(pPrev.Position - pCur.Position,
+                                           pNext.Position - pCur.Position) *
+                      (180 / 3.14159265359);
+        if (angle <= 0 && angle >= 180)
+          continue;
+
+        // If any vertices are within this triangle
+        bool inTri = false;
+        for (int j = 0; j < int(iVerts.size()); j++) {
+          if (algorithm::inTriangle(iVerts[j].Position, pPrev.Position,
+                                    pCur.Position, pNext.Position) &&
+              iVerts[j].Position != pPrev.Position &&
+              iVerts[j].Position != pCur.Position &&
+              iVerts[j].Position != pNext.Position) {
+            inTri = true;
+            break;
+          }
+        }
+        if (inTri)
+          continue;
+
+        // Create a triangle from pCur, pPrev, pNext
+        for (int j = 0; j < int(iVerts.size()); j++) {
+          if (iVerts[j].Position == pCur.Position)
+            oIndices.push_back(j);
+          if (iVerts[j].Position == pPrev.Position)
+            oIndices.push_back(j);
+          if (iVerts[j].Position == pNext.Position)
+            oIndices.push_back(j);
+        }
+
+        // Delete pCur from the list
+        for (int j = 0; j < int(tVerts.size()); j++) {
+          if (tVerts[j].Position == pCur.Position) {
+            tVerts.erase(tVerts.begin() + j);
+            break;
+          }
+        }
+
+        // reset i to the start
+        // -1 since loop will add 1 to it
+        i = -1;
+      }
+
+      // if no triangles were created
+      if (oIndices.size() == 0)
+        break;
+
+      // if no more vertices
+      if (tVerts.size() == 0)
+        break;
+    }
+  }
+
+  // Load Materials from .mtl file
+  bool LoadMaterials(std::string path) {
+    // If the file is not a material file return false
+    if (path.substr(path.size() - 4, path.size()) != ".mtl")
+      return false;
+
+    std::ifstream file(path);
+
+    // If the file is not found return false
+    if (!file.is_open())
+      return false;
+
+    Material tempMaterial;
+
+    bool listening = false;
+
+    // Go through each line looking for material variables
+    std::string curline;
+    while (std::getline(file, curline)) {
+      // new material and material name
+      if (algorithm::firstToken(curline) == "newmtl") {
+        if (!listening) {
+          listening = true;
+
+          if (curline.size() > 7) {
+            tempMaterial.name = algorithm::tail(curline);
+          } else {
+            tempMaterial.name = "none";
+          }
+        } else {
+          // Generate the material
+
+          // Push Back loaded Material
+          LoadedMaterials.push_back(tempMaterial);
+
+          // Clear Loaded Material
+          tempMaterial = Material();
+
+          if (curline.size() > 7) {
+            tempMaterial.name = algorithm::tail(curline);
+          } else {
+            tempMaterial.name = "none";
+          }
+        }
+      }
+      // Ambient Color
+      if (algorithm::firstToken(curline) == "Ka") {
+        std::vector<std::string> temp;
+        algorithm::split(algorithm::tail(curline), temp, " ");
+
+        if (temp.size() != 3)
+          continue;
+
+        tempMaterial.Ka.X = std::stof(temp[0]);
+        tempMaterial.Ka.Y = std::stof(temp[1]);
+        tempMaterial.Ka.Z = std::stof(temp[2]);
+      }
+      // Diffuse Color
+      if (algorithm::firstToken(curline) == "Kd") {
+        std::vector<std::string> temp;
+        algorithm::split(algorithm::tail(curline), temp, " ");
+
+        if (temp.size() != 3)
+          continue;
+
+        tempMaterial.Kd.X = std::stof(temp[0]);
+        tempMaterial.Kd.Y = std::stof(temp[1]);
+        tempMaterial.Kd.Z = std::stof(temp[2]);
+      }
+      // Specular Color
+      if (algorithm::firstToken(curline) == "Ks") {
+        std::vector<std::string> temp;
+        algorithm::split(algorithm::tail(curline), temp, " ");
+
+        if (temp.size() != 3)
+          continue;
+
+        tempMaterial.Ks.X = std::stof(temp[0]);
+        tempMaterial.Ks.Y = std::stof(temp[1]);
+        tempMaterial.Ks.Z = std::stof(temp[2]);
+      }
+      // Specular Exponent
+      if (algorithm::firstToken(curline) == "Ns") {
+        tempMaterial.Ns = std::stof(algorithm::tail(curline));
+      }
+      // Optical Density
+      if (algorithm::firstToken(curline) == "Ni") {
+        tempMaterial.Ni = std::stof(algorithm::tail(curline));
+      }
+      // Dissolve
+      if (algorithm::firstToken(curline) == "d") {
+        tempMaterial.d = std::stof(algorithm::tail(curline));
+      }
+      // Illumination
+      if (algorithm::firstToken(curline) == "illum") {
+        tempMaterial.illum = std::stoi(algorithm::tail(curline));
+      }
+      // Ambient Texture Map
+      if (algorithm::firstToken(curline) == "map_Ka") {
+        tempMaterial.map_Ka = algorithm::tail(curline);
+      }
+      // Diffuse Texture Map
+      if (algorithm::firstToken(curline) == "map_Kd") {
+        tempMaterial.map_Kd = algorithm::tail(curline);
+      }
+      // Specular Texture Map
+      if (algorithm::firstToken(curline) == "map_Ks") {
+        tempMaterial.map_Ks = algorithm::tail(curline);
+      }
+      // Specular Hightlight Map
+      if (algorithm::firstToken(curline) == "map_Ns") {
+        tempMaterial.map_Ns = algorithm::tail(curline);
+      }
+      // Alpha Texture Map
+      if (algorithm::firstToken(curline) == "map_d") {
+        tempMaterial.map_d = algorithm::tail(curline);
+      }
+      // Bump Map
+      if (algorithm::firstToken(curline) == "map_Bump" ||
+          algorithm::firstToken(curline) == "map_bump" ||
+          algorithm::firstToken(curline) == "bump") {
+        tempMaterial.map_bump = algorithm::tail(curline);
+      }
+    }
+
+    // Deal with last material
+
+    // Push Back loaded Material
+    LoadedMaterials.push_back(tempMaterial);
+
+    // Test to see if anything was loaded
+    // If not return false
+    if (LoadedMaterials.empty())
+      return false;
+    // If so return true
+    else
+      return true;
+  }
+};
+} // namespace objl
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b5bf8ae..d049449 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,5 +1,5 @@
-cmake_minimum_required(VERSION 3.9) 
-PROJECT(minus_renderer)
+cmake_minimum_required(VERSION 3.14) 
+project(minus_renderer)
 
 set(CMAKE_CXX_STANDARD 17)
 
@@ -12,4 +12,4 @@ add_subdirectory(${CMAKE_SOURCE_DIR}/src)
 # Third Party
 add_subdirectory(${CMAKE_SOURCE_DIR}/3rdparty/eigen)
 add_subdirectory(${CMAKE_SOURCE_DIR}/3rdparty/spdlog)
-
+add_subdirectory(${CMAKE_SOURCE_DIR}/3rdparty/OBJ_Loader)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ed5fdb7..4483314 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,8 +1,8 @@
-cmake_minimum_required(VERSION 3.9) 
+cmake_minimum_required(VERSION 3.14) 
 PROJECT(renderer)
 
 aux_source_directory(${CMAKE_CURRENT_SOURCE_DIR} CPP_FILES)
 
 add_executable(${CMAKE_PROJECT_NAME} ${CPP_FILES})
 target_include_directories(${CMAKE_PROJECT_NAME} PRIVATE ${OpenCV_INCLUDE_DIRS})
-target_link_libraries(${CMAKE_PROJECT_NAME} LINK_PUBLIC spdlog eigen ${OpenCV_LIBS})
+target_link_libraries(${CMAKE_PROJECT_NAME} LINK_PUBLIC OBJ_Loader spdlog eigen ${OpenCV_LIBS})
diff --git a/src/main.cc b/src/main.cc
index 901b71f..7d64874 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -9,36 +9,42 @@
 #include "camera.h"
 #include "common.h"
 #include "minus_renderer.h"
-#include "triangle.h"
 
 int main(int argc, char *argv[]) {
   spdlog::set_level(spdlog::level::trace);
 
-  const int resolution_width = 320;
-  const int resolution_height = 180;
-  const int far = -1000;
-  const int near = -10;
+  const std::string obj_path{argv[1]};
+
+  const int resolution_width = 640;
+  const int resolution_height = 320;
+  const int far = -100;
+  const int near = -1;
   const float fov = 120. / 180. * M_PI;
   const float aspect_ratio = static_cast<float>(resolution_width) /
                              static_cast<float>(resolution_height);
 
-  std::vector<Triangle> primitives{
-      {Point3d{10, 5, -20}, Point3d{30, 25, -20}, Point3d{55, 5, -20}}};
-
   cv::namedWindow("MinusRenderer", cv::WINDOW_AUTOSIZE);
   const auto start = std::chrono::steady_clock::now();
   Camera camera{Vector3d{0, 1, 0}};
-  MinusRenderer renderer{near, far, fov, aspect_ratio, std::move(primitives)};
+  MinusRenderer renderer{near, far, fov, aspect_ratio};
+  assert(renderer.load_mesh(obj_path));
+  TransformMatrix translate{
+      {1., 0., 0., 3.}, {0., 1., 0., 3.}, {0., 0., 1., -4.}, {0., 0., 0., 1.}};
+  renderer.model_transform(translate);
   do {
     const auto elapse = std::chrono::duration_cast<std::chrono::milliseconds>(
                             std::chrono::steady_clock::now() - start)
                             .count();
-    TransformMatrix m{
+    TransformMatrix move2origin{
+        {1, 0, 0, 0}, {0, 1, 0, -3}, {0, 0, 1, 4}, {0, 0, 0, 1}};
+    TransformMatrix rotation{
         Eigen::Transform<double, 3, Eigen::Affine>(
-            Eigen::AngleAxisd(0.03125 * M_PI, Eigen::Vector3d::UnitZ())
+            Eigen::AngleAxisd(0.0625 * M_PI, Eigen::Vector3d::UnitX())
                 .toRotationMatrix())
             .matrix()};
-    renderer.model_transform(m);
+    TransformMatrix move_back{
+        {1, 0, 0, 0}, {0, 1, 0, 3}, {0, 0, 1, -4}, {0, 0, 0, 1}};
+    renderer.model_transform(move_back * rotation * move2origin);
 
     camera.set_position(Point3d{0, 0, 0});
     camera.set_gaze(Vector3d{0, 0, -1});
diff --git a/src/minus_renderer.cc b/src/minus_renderer.cc
index 5fbf3e9..604cec7 100644
--- a/src/minus_renderer.cc
+++ b/src/minus_renderer.cc
@@ -3,16 +3,18 @@
 
 #include <cmath>
 
-#include "minus_renderer.h"
-#include "rasterizer.h"
+#include "OBJ_Loader.h"
 #include "spdlog/spdlog.h"
 
+#include "minus_renderer.h"
+#include "rasterizer.h"
+
 MinusRenderer::MinusRenderer(float near, float far, float fov,
-                             float aspect_ratio,
-                             std::vector<Triangle> &&primitives)
+                             float aspect_ratio)
     : near_(near), far_(far), fov_(fov), aspect_ratio_(aspect_ratio),
-      primitives_(primitives) {
-  spdlog::info("fov: {}, aspect ratio: {}", fov_, aspect_ratio_);
+      projection_matrix_(orthographic_tranform() * squish_tranform()) {
+  spdlog::info("near: {}, far: {}, fov: {}, aspect ratio: {}", near_, far_,
+               fov_, aspect_ratio_);
 }
 
 float MinusRenderer::calculate_height(const float fov, const float near) {
@@ -61,45 +63,84 @@ TransformMatrix MinusRenderer::view_port_transform(const float width,
                          {0, 0, 0, 1}};
 }
 
-void MinusRenderer::model_transform(const TransformMatrix &m) {
-  for (auto &t : primitives_) {
-    t.affine_transform(m);
+void MinusRenderer::model_transform(const TransformMatrix &mtx) {
+  for (auto &m : meshes_) {
+    for (auto &p : m) {
+      p.affine_transform(mtx);
+    }
   }
 }
 
-void MinusRenderer::view_transform(const TransformMatrix &m) {
-  for (auto &t : primitives_) {
-    t.affine_transform(m);
+void MinusRenderer::view_transform(const TransformMatrix &mtx) {
+  for (auto &m : meshes_) {
+    for (auto &p : m) {
+      p.affine_transform(mtx);
+    }
   }
 }
 
 cv::Mat MinusRenderer::render(const int resolution_width,
                               const int resolution_height) {
-  std::vector<Triangle> primitives = primitives_;
-  for (auto &t : primitives) {
-    t.projective_transform(squish_tranform());
-    t.affine_transform(orthographic_tranform());
-    t.affine_transform(
-        view_port_transform(resolution_width, resolution_height));
-  }
-
+  std::vector<Mesh> meshes = meshes_;
   Rasterizer rasterizer{resolution_width, resolution_height};
-  const auto &pixels = rasterizer.rasterize(primitives);
+  for (auto &primitives : meshes) {
+    for (int i = 0; i < primitives.size(); ++i) {
+      auto &t = primitives[i];
+      t.projective_transform(projection_matrix_);
+      t.affine_transform(
+          view_port_transform(resolution_width, resolution_height));
+    }
+    rasterizer.rasterize(primitives);
+  }
+  const auto &pixels = rasterizer.get_picture();
   assert(!pixels.empty());
   cv::Mat color_image(pixels.size(), (pixels[0]).size(), CV_8UC3);
+  cv::Mat depth_image(pixels.size(), (pixels[0]).size(), CV_8UC1);
 
   for (int i = 0; i < pixels.size(); ++i) {
     const auto &row = pixels[i];
     for (int j = 0; j < row.size(); ++j) {
       auto &pixel_color = color_image.at<cv::Vec3b>(pixels.size() - i - 1, j);
-      pixel_color =
-          cv::Vec3b{static_cast<unsigned char>(
-                        row[j].x() * std::numeric_limits<char>::max()),
-                    static_cast<unsigned char>(
-                        row[j].y() * std::numeric_limits<char>::max()),
-                    static_cast<unsigned char>(
-                        row[j].z() * std::numeric_limits<char>::max())};
+      pixel_color = cv::Vec3b{
+          static_cast<unsigned char>(row[j].x() *
+                                     std::numeric_limits<unsigned char>::max()),
+          static_cast<unsigned char>(row[j].y() *
+                                     std::numeric_limits<unsigned char>::max()),
+          static_cast<unsigned char>(
+              row[j].z() * std::numeric_limits<unsigned char>::max())};
+      auto &pixel_depth =
+          depth_image.at<unsigned char>(pixels.size() - i - 1, j);
+      if (std::isinf(std::fabs(row[j].w()))) {
+        pixel_depth = 0;
+      } else {
+        const float k =
+            static_cast<float>(std::numeric_limits<unsigned char>::max()) * 0.5;
+        const float b = k;
+        pixel_depth = static_cast<unsigned char>(k * row[j].w() + b);
+      }
     }
   }
-  return color_image;
+  return depth_image;
+}
+
+bool MinusRenderer::load_mesh(const std::string &file_path) {
+  objl::Loader loader{};
+  bool load_status = loader.LoadFile(file_path.c_str());
+  if (load_status) {
+    for (const auto &mesh : loader.LoadedMeshes) {
+      spdlog::info("Current mesh has {} vertices.", mesh.Vertices.size());
+      Mesh m{};
+      for (int i = 0; i < mesh.Vertices.size(); i += 3) {
+        const auto objl_v0 = mesh.Vertices[mesh.Indices[i]].Position;
+        Point3d v0(objl_v0.X, objl_v0.Y, objl_v0.Z);
+        const auto objl_v1 = mesh.Vertices[mesh.Indices[i + 1]].Position;
+        Point3d v1(objl_v1.X, objl_v1.Y, objl_v1.Z);
+        const auto objl_v2 = mesh.Vertices[mesh.Indices[i + 2]].Position;
+        Point3d v2(objl_v2.X, objl_v2.Y, objl_v2.Z);
+        m.push_back(Triangle(v0, v1, v2));
+      }
+      meshes_.push_back(m);
+    }
+  }
+  return load_status;
 }
diff --git a/src/minus_renderer.h b/src/minus_renderer.h
index f127097..bb0ff90 100644
--- a/src/minus_renderer.h
+++ b/src/minus_renderer.h
@@ -8,15 +8,20 @@
 
 class MinusRenderer {
 public:
-  MinusRenderer(float near, float far, float fov, float aspect_ratio,
-                std::vector<Triangle> &&primitives);
+  MinusRenderer(float near, float far, float fov, float aspect_ratio);
 
 public:
   cv::Mat render(const int resolution_width, const int resolution_height);
 
 public:
-  void model_transform(const TransformMatrix &m);
-  void view_transform(const TransformMatrix &m);
+  bool load_mesh(const std::string &file_path);
+
+public:
+  void model_transform(const TransformMatrix &mtx);
+  void view_transform(const TransformMatrix &mtx);
+
+private:
+  using Mesh = std::vector<Triangle>;
 
 private:
   static float calculate_height(const float fov, const float near);
@@ -34,5 +39,7 @@ private:
   float fov_; // In radian
   float aspect_ratio_;
 
-  std::vector<Triangle> primitives_;
+  TransformMatrix projection_matrix_;
+
+  std::vector<Mesh> meshes_;
 };
diff --git a/src/rasterizer.cc b/src/rasterizer.cc
index 090104b..7fa582c 100644
--- a/src/rasterizer.cc
+++ b/src/rasterizer.cc
@@ -2,9 +2,11 @@
 // Author: tianlei.richard@qq.com (tianlei.richard)
 
 #include <limits>
+#include <thread>
+
+#include "spdlog/spdlog.h"
 
 #include "rasterizer.h"
-#include "spdlog/spdlog.h"
 
 Rasterizer::Rasterizer(const int width, const int height)
     : width_(width), height_(height),
@@ -25,24 +27,54 @@ Rasterizer::rasterize(const std::vector<Triangle> &primitives) {
     const int y_min = std::min(std::max(0, aabb.y), height_);
     const int y_max = std::min(std::max(0, aabb.y + aabb.height), height_);
 
-    // spdlog::info("Range, ({},{},{},{})", x_min, x_max, y_min, y_max);
-    for (int i = y_min; i < y_max; ++i) {
-      for (int j = x_min; j < x_max; ++j) {
-        auto &property = pixels_[i][j];
-        const auto &[is_inside, z_screen] =
-            inside(Point2d{j + 0.5, i + 0.5}, t);
-        if (is_inside && z_screen > property[3]) {
-          property[0] = 0;
-          property[1] = 1;
-          property[2] = 0;
-          property[3] = z_screen;
+    auto partial_rasterize = [x_min, x_max, &t, this](const int start,
+                                                      const int end) {
+      for (int i = start; i < end; ++i) {
+        for (int j = x_min; j <= x_max; ++j) {
+          auto &property = (this->pixels_)[i][j];
+          const auto &[is_inside, z_screen] =
+              inside(Point2d{j + 0.5, i + 0.5}, t);
+          if (is_inside && z_screen > property[3]) {
+            property[0] = 0;
+            property[1] = 1;
+            property[2] = 0;
+            property[3] = z_screen;
+          }
         }
       }
+    };
+    const int thread_num = 6;
+    std::vector<std::thread> rasterize_threads;
+    for (int start = y_min, offset = (y_max - y_min) / thread_num,
+             remainer = (y_max - y_min) % thread_num;
+         start < y_max;) {
+      int end = start + offset;
+      if (remainer > 0) {
+        end += 1;
+        remainer -= 1;
+      }
+      rasterize_threads.push_back(std::thread(partial_rasterize, start, end));
+      start = end;
+    }
+    for (auto &t : rasterize_threads) {
+      t.join();
     }
   }
   return pixels_;
 }
 
+std::vector<std::vector<PixelProperty>> Rasterizer::get_picture() const {
+  return pixels_;
+}
+
+void Rasterizer::reset() {
+  pixels_ = std::vector<std::vector<PixelProperty>>(
+      height_,
+      std::vector<PixelProperty>(
+          width_,
+          PixelProperty{0., 0., 0., -std::numeric_limits<double>::infinity()}));
+}
+
 std::pair<bool, double> Rasterizer::inside(const Point2d &p_screen,
                                            const Triangle &t) {
   const auto points = t.get_points<3>();
diff --git a/src/rasterizer.h b/src/rasterizer.h
index ebc21e4..36d172d 100644
--- a/src/rasterizer.h
+++ b/src/rasterizer.h
@@ -8,12 +8,6 @@
 #include <vector>
 
 class Rasterizer {
-public:
-  typedef enum RasterizeMode {
-    Fill,
-    Line,
-  } RasterizeMode;
-
 public:
   Rasterizer(const int width, const int height);
 
@@ -21,6 +15,9 @@ public:
   std::vector<std::vector<PixelProperty>>
   rasterize(const std::vector<Triangle> &primitives);
 
+  std::vector<std::vector<PixelProperty>> get_picture() const;
+  void reset();
+
 private:
   static std::pair<bool, double> inside(const Point2d &p_screen,
                                         const Triangle &t);
@@ -29,6 +26,4 @@ private:
   int width_;
   int height_;
   std::vector<std::vector<PixelProperty>> pixels_;
-
-  RasterizeMode mode_ = RasterizeMode::Fill;
 };
diff --git a/src/triangle.cc b/src/triangle.cc
index d084fd5..dcb480f 100644
--- a/src/triangle.cc
+++ b/src/triangle.cc
@@ -19,7 +19,7 @@ BBox Triangle::axis_align_bbox() const {
   const int y_min = std::min({points_[0].y(), points_[1].y(), points_[2].y()});
   const int x_max = std::max({points_[0].x(), points_[1].x(), points_[2].x()});
   const int y_max = std::max({points_[0].y(), points_[1].y(), points_[2].y()});
-  return BBox{x_min, y_min, (x_max - x_min), (y_max - y_min)};
+  return BBox{x_min, y_min, (x_max - x_min + 1), (y_max - y_min + 1)};
 }
 
 Vector3d Triangle::normal_vector() const {