diff --git a/.clang-format b/.clang-format
new file mode 100644
index 0000000000000000000000000000000000000000..503f80e217906435071464b965eaac25b2e6a48f
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,7 @@
+BasedOnStyle: Google
+IndentWidth: 4
+ColumnLimit: 100
+BreakConstructorInitializersBeforeComma: true
+AccessModifierOffset: -4
+KeepEmptyLinesAtTheStartOfBlocks: true
+SortIncludes: false
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..e7a3c00a2db611fd2a15350ecda449e93ce9da1e
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,20 @@
+.DS_Store
+.idea
+.vscode
+
+build
+cmake-build-debug
+cmake-build-minsizerel
+cmake-build-release
+cmake-build-relwithdebinfo
+
+*.html
+# Build files from LaTeX
+*.aux
+*.bbl
+*.blg
+*.fdb_latexmk
+*.fls
+*.log
+*-blx.bib
+*.run.xml
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000000000000000000000000000000000000..7394ceb69f88c1e75b761fd93fe3bc8921c049f9
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,19 @@
+# Use the official gcc image
+image: gcc
+
+stages:
+    - build
+
+before_script:
+    # Install build dependencies
+    - apt-get update && apt-get -y install cmake freeglut3 freeglut3-dev libglew-dev libwxgtk3.0-dev
+    - gcc --version
+
+MoA:
+    # Build the application
+    stage: build
+    script:
+        - mkdir build
+        - cd build
+        - cmake ..
+        - make
\ No newline at end of file
diff --git a/Decimation/QuadricDecimationMesh.cpp b/Decimation/QuadricDecimationMesh.cpp
index 0d156f46fbe3d71cea32b5ac2146a02b86a5d5dc..af08eb099ab0d124a24c73fe0a3d55c3ab837430 100644
--- a/Decimation/QuadricDecimationMesh.cpp
+++ b/Decimation/QuadricDecimationMesh.cpp
@@ -11,33 +11,32 @@
  *************************************************************************************************/
 #include "QuadricDecimationMesh.h"
 
-const QuadricDecimationMesh::VisualizationMode
-    QuadricDecimationMesh::QuadricIsoSurfaces =
-        NewVisualizationMode("Quadric Iso Surfaces");
+const QuadricDecimationMesh::VisualizationMode QuadricDecimationMesh::QuadricIsoSurfaces =
+    NewVisualizationMode("Quadric Iso Surfaces");
 
 void QuadricDecimationMesh::Initialize() {
-  // Allocate memory for the quadric array
-  size_t numVerts = mVerts.size();
-  mQuadrics.reserve(numVerts);
-  std::streamsize width = std::cerr.precision(); // store stream precision
-  for (size_t i = 0; i < numVerts; i++) {
+    // Allocate memory for the quadric array
+    size_t numVerts = mVerts.size();
+    mQuadrics.reserve(numVerts);
+    std::streamsize width = std::cerr.precision();  // store stream precision
+    for (size_t i = 0; i < numVerts; i++) {
 
-    // Compute quadric for vertex i here
-    mQuadrics.push_back(createQuadricForVert(i));
+        // Compute quadric for vertex i here
+        mQuadrics.push_back(createQuadricForVert(i));
 
-    // Calculate initial error, should be numerically close to 0
+        // Calculate initial error, should be numerically close to 0
 
-    Vector3<float> v0 = mVerts[i].pos;
-    Vector4<float> v(v0[0], v0[1], v0[2], 1);
-    Matrix4x4<float> m = mQuadrics.back();
+        Vector3<float> v0 = mVerts[i].pos;
+        Vector4<float> v(v0[0], v0[1], v0[2], 1);
+        Matrix4x4<float> m = mQuadrics.back();
 
-    float error = v * (m * v);
-    // std::cerr << std::scientific << std::setprecision(2) << error << " ";
-  }
-  std::cerr << std::setprecision(width) << std::fixed; // reset stream precision
+        float error = v * (m * v);
+        // std::cerr << std::scientific << std::setprecision(2) << error << " ";
+    }
+    std::cerr << std::setprecision(width) << std::fixed;  // reset stream precision
 
-  // Run the initialize for the parent class to initialize the edge collapses
-  DecimationMesh::Initialize();
+    // Run the initialize for the parent class to initialize the edge collapses
+    DecimationMesh::Initialize();
 }
 
 /*! \lab2 Implement the computeCollapse here */
@@ -46,110 +45,114 @@ void QuadricDecimationMesh::Initialize() {
  * DecimationMesh::EdgeCollapse
  */
 void QuadricDecimationMesh::computeCollapse(EdgeCollapse *collapse) {
-  // Compute collapse->position and collapse->cost here
-  // based on the quadrics at the edge endpoints
-
-    size_t edge1 = collapse->halfEdge;
-    size_t edge2 = e(edge1).pair;
-    size_t v0 = e(edge1).vert;
-    size_t v1 = e(edge2).vert;
-
-    Matrix4x4<float> Q_saved = mQuadrics[v0] + mQuadrics[v1]; //createQuadricForVert(v0) + createQuadricForVert(v1) ;
-    Matrix4x4<float> Q = Q_saved;
-    Q(3,0) = 0.0f;
-    Q(3,1) = 0.0f;
-    Q(3,2) = 0.0f;
-    Q(3,3) = 1.0f;
-
-    Vector4<float> v_;
-
-    //std::cout << "5" << std::endl;
-   if(!Q.IsSingular()){
-        v_ = Q.Inverse() * Vector4<float>(0.0f,0.0f,0.0f,1.0f);
-        collapse->position = Vector3<float>(v_[0], v_[1], v_[2]);
-        collapse->cost = v_* ( Q_saved * v_);
-
-    }
-    else {
+    // Compute collapse->position and collapse->cost here
+    // based on the quadrics at the edge endpoints
+
+    size_t v0 = e(collapse->halfEdge).vert;
+    size_t v1 = e(e(collapse->halfEdge).pair).vert;
+
+    Matrix4x4<float> Q1 = createQuadricForVert(v0);
+    Matrix4x4<float> Q2 = createQuadricForVert(v1);
+    Matrix4x4<float> Q = Q1 + Q2;
+    Matrix4x4<float> Q_saved = Q;
+    Q(3, 0) = 0.0f;
+    Q(3, 1) = 0.0f;
+    Q(3, 2) = 0.0f;
+    Q(3, 3) = 1.0f;
+    Vector4<float> v;
+
+    // std::cout << "5" << std::endl;
+    if (!Q.IsSingular()) {
+        v = Q.Inverse() * Vector4<float>(0.0f, 0.0f, 0.0f, 1.0f);
+        collapse->position = Vector3<float>(v[0], v[1], v[2]);
+       
+        collapse->cost = v * (Q_saved * v);
+        // std::cout << "6" << std::endl;
+    } else {
         collapse->position = (mVerts[v0].pos + mVerts[v1].pos) * 0.5;
         collapse->cost = (collapse->position - mVerts[v0].pos).Length();
-
+        // std::cout << "7" << std::endl;
     }
+
+    // std::cerr << "computeCollapse in QuadricDecimationMesh not implemented.\n";
 }
 
 /*! After each edge collapse the vertex properties need to be updated */
 void QuadricDecimationMesh::updateVertexProperties(size_t ind) {
-  DecimationMesh::updateVertexProperties(ind);
-  mQuadrics[ind] = createQuadricForVert(ind);
+    DecimationMesh::updateVertexProperties(ind);
+    mQuadrics[ind] = createQuadricForVert(ind);
 }
 
 /*!
  * \param[in] indx vertex index, points into HalfEdgeMesh::mVerts
  */
-Matrix4x4<float>
-QuadricDecimationMesh::createQuadricForVert(size_t indx) const {
+Matrix4x4<float> QuadricDecimationMesh::createQuadricForVert(size_t indx) const {
     float q[4][4] = {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};
-  Matrix4x4<float> Q(q);
-
-  // The quadric for a vertex is the sum of all the quadrics for the adjacent
-  // faces Tip: Matrix4x4 has an operator +=
+    Matrix4x4<float> Q(q);
+    // std::cout << "3" << std::endl;
+    // The quadric for a vertex is the sum of all the quadrics for the adjacent
+    // faces Tip: Matrix4x4 has an operator +=
+    // size_t start = v(indx).edge;
+    // size_t temp = start;
+    // do {
+    //    size_t face = e(temp).face;
+    //    Q += createQuadricForFace(face);
+    //    temp = e(e(temp).prev).pair;
+    //} while (start != temp);
     std::vector<size_t> faces = FindNeighborFaces(indx);
-    for(size_t face : faces) {
+    for (size_t face : faces) {
         Q += createQuadricForFace(face);
     }
-  return Q;
+    // std::cout << "4" << std::endl;
+
+    return Q;
 }
 
 /*!
  * \param[in] indx face index, points into HalfEdgeMesh::mFaces
  */
-Matrix4x4<float>
-QuadricDecimationMesh::createQuadricForFace(size_t indx) const {
-
-  // Calculate the quadric (outer product of plane parameters) for a face
-  // here using the formula from Garland and Heckbert
-
-  //1. Get two vectors (B-A) and (C-A)
-  //2. (B-A)x(C-A) is the normal. Normalize: a, b, c are done
-  float a, b, c, d;
-  Vector3<float> normal = f(indx).normal;
-  normal = normal.Normalize();
-  Vector3<float> vertex = v(e(f(indx).edge).vert).pos;
-
-  a = normal[0];
-  b = normal[1];
-  c = normal[2];
-
-    //3. solve for d: d = -n1*x - n2*y - n3*z = - (n1*x + n2*y + n3*z) = - (n*vertex)
-  d = -1.0f*(normal*vertex);
-
-  float m[4][4] =
-    {
-        {a*a, a*b, a*c, a*d},
-        {a*b, b*b, b*c, b*d},
-        {a*c, b*c, c*c, c*d},
-        {a*d, b*d, c*d, d*d}
-    };
-  //Calculate the matrix
-  Matrix4x4<float> Kp(m);
-
-  return Kp;
+Matrix4x4<float> QuadricDecimationMesh::createQuadricForFace(size_t indx) const {
+
+    // Calculate the quadric (outer product of plane parameters) for a face
+    // here using the formula from Garland and Heckbert
+
+    // 1. Get two vectors (B-A) and (C-A)
+    // 2. (B-A)x(C-A) is the normal. Normalize: a, b, c are done
+    float a, b, c, d;
+    Vector3<float> normal = f(indx).normal;
+    Vector3<float> vertex = v(e(f(indx).edge).vert).pos;
+    normal.Normalize();
+    a = normal[0];
+    b = normal[1];
+    c = normal[2];
+    // std::cout << "1" << std::endl;
+    // 3. solve for d: d = -n1*x - n2*y - n3*z = - (n1*x + n2*y + n3*z) = - (n*vertex)
+    d = -1.0f * (normal * vertex);
+
+    float m[4][4] = {{a * a, a * b, a * c, a * d},
+                     {a * b, b * b, b * c, b * d},
+                     {a * c, b * c, c * c, c * d},
+                     {a * d, b * d, c * d, d * d}};
+    // Calculate the matrix
+    Matrix4x4<float> Kp(m);
+    // std::cout << "2" << std::endl;
+    return Kp;
 }
 
 void QuadricDecimationMesh::Render() {
-  DecimationMesh::Render();
+    DecimationMesh::Render();
 
-  glEnable(GL_LIGHTING);
-  glMatrixMode(GL_MODELVIEW);
+    glEnable(GL_LIGHTING);
+    glMatrixMode(GL_MODELVIEW);
 
-  if (mVisualizationMode == QuadricIsoSurfaces) {
-    // Apply transform
-    glPushMatrix(); // Push modelview matrix onto stack
+    if (mVisualizationMode == QuadricIsoSurfaces) {
+        // Apply transform
+        glPushMatrix();  // Push modelview matrix onto stack
 
-    // Implement the quadric visualization here
-    std::cout << "Quadric visualization not implemented" << std::endl;
+        // Implement the quadric visualization here
+        std::cout << "Quadric visualization not implemented" << std::endl;
 
-    // Restore modelview matrix
-    glPopMatrix();
-  }
+        // Restore modelview matrix
+        glPopMatrix();
+    }
 }
diff --git a/Decimation/SimpleDecimationMesh.cpp b/Decimation/SimpleDecimationMesh.cpp
index 1fee1cd517da9654babe93b12551e04c04586f01..3613898aaed85a1afaf762989c36b945bfb8a79d 100644
--- a/Decimation/SimpleDecimationMesh.cpp
+++ b/Decimation/SimpleDecimationMesh.cpp
@@ -6,7 +6,7 @@
  *   Ken Museth (ken.museth@itn.liu.se)
  *   Michael Bang Nielsen (bang@daimi.au.dk)
  *   Ola Nilsson (ola.nilsson@itn.liu.se)
- *   Andreas S�derstr�m (andreas.soderstrom@itn.liu.se)
+ *   Andreas S�derstr�m (andreas.soderstrom@itn.liu.se)
  *
  *************************************************************************************************/
 #include "SimpleDecimationMesh.h"
@@ -20,26 +20,5 @@ void SimpleDecimationMesh::computeCollapse(EdgeCollapse *collapse) {
   const Vector3<float> &v1 =
       mVerts[mEdges[mEdges[collapse->halfEdge].pair].vert].pos;
   collapse->position = (v0 + v1) * 0.5;
-  //collapse->cost = (collapse->position - v0).Length();
-
-    Vector3<float> CameraPos(0.0f,0.0f,0.0f);
-    Vector3<float> CameraForw(0.0f,0.0f,-1.0f);
-    //Distance to camera
-    //collapse->cost = log(1/(collapse->position - CameraPos).Length());
-
-    //Visible to camera
-    size_t vert = mEdges[collapse->halfEdge].vert;
-
-   // float cost = v(vert).normal * CameraForw;
-    std::vector<size_t> faces = FindNeighborFaces(vert);
-    std::vector<size_t> verts  = FindNeighborVertices(vert);
-
-    float cost{};
-    for(size_t vert : verts) {
-         cost = v(vert).normal * CameraForw;
-         if(cost > 0 ) {
-             collapse->cost = (collapse->position - v0).Length();
-         }
-    }
-
+  collapse->cost = (collapse->position - v0).Length();
 }
diff --git a/GUI/FrameMain.cpp b/GUI/FrameMain.cpp
index eea4161168dc0f7a3ea81dee4151249b6a759574..5c7fdff6754011719f8d0502ca9d65872f581b18 100644
--- a/GUI/FrameMain.cpp
+++ b/GUI/FrameMain.cpp
@@ -741,6 +741,10 @@ void FrameMain::AddObjectImplicitMesh(wxCommandEvent &event) {
 void FrameMain::AddObjectQuadricPlane(wxCommandEvent &event) {
   Matrix4x4<float> M;
   // Construct the quadric matrix here
+  M(0,0) = M(1,1) = M(2,2) = M(3,3) = 0.0f;
+  M(0,3) = M(3,0) = 0.5f;
+  M(1,3) = M(3,1) = 0.5f;
+  M(2,3) = M(3,2) = 0.5f;
 
   Quadric *Q = new Quadric(M);
   Q->SetBoundingBox(Bbox(-1, 1));
@@ -755,6 +759,8 @@ void FrameMain::AddObjectQuadricPlane(wxCommandEvent &event) {
 void FrameMain::AddObjectQuadricCylinder(wxCommandEvent &event) {
   Matrix4x4<float> M;
   // Construct the quadric matrix here
+  M(2,2) = 0.0f;
+  M(3,3) = -1.0f;
 
   Quadric *Q = new Quadric(M);
   Q->SetBoundingBox(Bbox(-1, 1));
@@ -769,6 +775,8 @@ void FrameMain::AddObjectQuadricCylinder(wxCommandEvent &event) {
 void FrameMain::AddObjectQuadricEllipsoid(wxCommandEvent &event) {
   Matrix4x4<float> M;
   // Construct the quadric matrix here
+  M(1,1) = 1.0f/0.25f;
+  M(3,3) = -1.0f;
 
   Quadric *Q = new Quadric(M);
   Q->SetBoundingBox(Bbox(-1, 1));
@@ -783,6 +791,8 @@ void FrameMain::AddObjectQuadricEllipsoid(wxCommandEvent &event) {
 void FrameMain::AddObjectQuadricCone(wxCommandEvent &event) {
   Matrix4x4<float> M;
   // Construct the quadric matrix here
+  M(2,2) = -1.0f;
+  M(3,3) = 0.0f;
 
   Quadric *Q = new Quadric(M);
   Q->SetBoundingBox(Bbox(-1, 1));
@@ -797,6 +807,15 @@ void FrameMain::AddObjectQuadricCone(wxCommandEvent &event) {
 void FrameMain::AddObjectQuadricParaboloid(wxCommandEvent &event) {
   Matrix4x4<float> M;
   // Construct the quadric matrix here
+  //Version 1
+  /* M(0,0) = M(1,1)= 1.0f;
+  M(2,3) = M(3,2) = -1.0f;
+  M(2,2) = M(3,3) = 0.0f; */
+
+  //Version 2
+  M(0,0) = 1.0f;
+  M(1,1) = M(2,3) = M(3,2) = -1.0f;
+  M(2,2) = M(3,3) = 0.0f;
 
   Quadric *Q = new Quadric(M);
   Q->SetBoundingBox(Bbox(-2, 2));
@@ -811,6 +830,20 @@ void FrameMain::AddObjectQuadricParaboloid(wxCommandEvent &event) {
 void FrameMain::AddObjectQuadricHyperboloid(wxCommandEvent &event) {
   Matrix4x4<float> M;
   // Construct the quadric matrix here
+  //Hyperboloid 1
+  /*M(0, 0) = M(1, 1) = 2.0f;
+  M(2, 2) = -2.0f;
+  M(3, 3) = 0.0f;*/
+
+  //Hyperboloid 2
+  /*M(0, 0) = M(1, 1) = 1.0f;
+  M(2, 2) = -1.0f;
+  M(3, 3) = 1.0f;*/
+
+  //Hyperboloid 3
+  M(0, 0) = M(1, 1) = 1.0f;
+  M(2, 2) = -1.0f;
+  M(3, 3) = -1.0f;
 
   Quadric *Q = new Quadric(M);
   Q->SetBoundingBox(Bbox(-2, 2));
diff --git a/Geometry/CSG.h b/Geometry/CSG.h
index 5741fbc68a8f8b9877de6b1c769b8379f33102a5..ec84b5fcb236bc4de56ca584f726f1e90d051a75 100644
--- a/Geometry/CSG.h
+++ b/Geometry/CSG.h
@@ -14,7 +14,7 @@
 #define __CSG_H__
 
 #include "Geometry/Implicit.h"
-
+#include <cmath>
 /*! \brief CSG Operator base class */
 class CSG_Operator : public Implicit {
 
@@ -43,7 +43,19 @@ public:
     // and can be transformed like all implicit surfaces.
     // Then, get values from left and right children and perform the
     // boolean operation.
-    return 0;
+    float p = 7.0f;
+    TransformW2O(x, y, z);
+
+    float A = left->GetValue(x, y, z);
+    float B = right->GetValue(x, y, z);
+    
+    float DA = exp(-A);
+    float DB = exp(-B);
+
+    float DAUB = pow((pow(DA, p) + pow(DB, p)), (1.0f / p));
+
+
+    return -log(DAUB);
   }
 };
 
@@ -54,7 +66,21 @@ public:
     mBox = BoxIntersection(l->GetBoundingBox(), r->GetBoundingBox());
   }
 
-  virtual float GetValue(float x, float y, float z) const { return 0; }
+  virtual float GetValue(float x, float y, float z) const { 
+      
+    float p = 7.0f;
+    TransformW2O(x, y, z);
+
+    float A = left->GetValue(x, y, z);
+    float B = right->GetValue(x, y, z);
+
+    float DA = exp(-A);
+    float DB = exp(-B);
+
+    float DAUB = pow((pow(DA, -p) + pow(DB, -p)), (-1.0f / p));
+
+    return -log(DAUB);
+  }
 };
 
 /*! \brief Difference boolean operation */
@@ -64,7 +90,20 @@ public:
     mBox = l->GetBoundingBox();
   }
 
-  virtual float GetValue(float x, float y, float z) const { return 0; }
+  virtual float GetValue(float x, float y, float z) const {
+    float p = 7.0f;
+    TransformW2O(x, y, z);
+
+    float A = left->GetValue(x, y, z);
+    float B = right->GetValue(x, y, z);
+
+    float DA = exp(-A);
+    float DB = exp(B);
+
+    float DAUB = pow((pow(DA, -p) + pow(DB, -p)), (-1.0f / p));
+
+    return -log(DAUB);
+  }
 };
 
 /*! \brief BlendedUnion boolean operation */
diff --git a/Geometry/HalfEdgeMesh.cpp b/Geometry/HalfEdgeMesh.cpp
index 995bf832f3cf5f599051bfa2debb4ac44bd7e5cf..91c4a0a8dd541f5dfb90a6311afbe33d45d3638d 100644
--- a/Geometry/HalfEdgeMesh.cpp
+++ b/Geometry/HalfEdgeMesh.cpp
@@ -12,8 +12,7 @@
 #include "Geometry/HalfEdgeMesh.h"
 
 const size_t HalfEdgeMesh::BORDER = (std::numeric_limits<size_t>::max)();
-const size_t HalfEdgeMesh::UNINITIALIZED =
-    (std::numeric_limits<size_t>::max)() - 1;
+const size_t HalfEdgeMesh::UNINITIALIZED = (std::numeric_limits<size_t>::max)() - 1;
 
 HalfEdgeMesh::HalfEdgeMesh() {}
 
@@ -26,49 +25,43 @@ HalfEdgeMesh::~HalfEdgeMesh() {}
  * \param[in] v3 vertex 3, Vector3<float>
  */
 bool HalfEdgeMesh::AddFace(const std::vector<Vector3<float> > &verts) {
-  // Add your code here
-  //std::cerr << "ADD TRIANGLE NOT IMPLEMENTED. ";
+    // Add your code here
+    //std::cerr << "ADD TRIANGLE NOT IMPLEMENTED. ";
 
-  // Add the vertices of the face/triangle
+    // Add the vertices of the face/triangle
     size_t ind1, ind2, ind3;
     AddVertex(verts.at(0), ind1);
     AddVertex(verts.at(1), ind2);
     AddVertex(verts.at(2), ind3);
 
-  // Add all half-edge pairs
-    size_t half_edge_out1, half_edge_out2, half_edge_out3;
-    size_t half_edge_in1, half_edge_in2, half_edge_in3;
-
-    AddHalfEdgePair(ind1, ind2, half_edge_in1, half_edge_out1);
-    AddHalfEdgePair(ind2, ind3, half_edge_in2, half_edge_out2);
-    AddHalfEdgePair(ind3, ind1, half_edge_in3, half_edge_out3);
-
-  // Connect inner ring
-    e(half_edge_in1).next = half_edge_in2;
-    e(half_edge_in2).next = half_edge_in3;
-    e(half_edge_in3).next = half_edge_in1;
-
-    e(half_edge_in1).prev = half_edge_in3;
-    e(half_edge_in2).prev = half_edge_in1;
-    e(half_edge_in3).prev = half_edge_in2;
-
-  // Finally, create the face, don't forget to set the normal (which should be
-  // normalized)
-    Face tri{};
-    tri.edge = half_edge_in1;
+    // Add all half-edge pairs
+    size_t halfEdgeIn1, halfEdgeIn2, halfEdgeIn3, halfEdgeOut1, halfEdgeOut2, halfEdgeOut3;
+
+    AddHalfEdgePair(ind1, ind2, halfEdgeIn1, halfEdgeOut1);
+    AddHalfEdgePair(ind2, ind3, halfEdgeIn2, halfEdgeOut2);
+    AddHalfEdgePair(ind3, ind1, halfEdgeIn3, halfEdgeOut3);
+    // Connect inner ring
+    e(halfEdgeIn1).next = halfEdgeIn2;
+    e(halfEdgeIn2).next = halfEdgeIn3;
+    e(halfEdgeIn3).next = halfEdgeIn1;
+
+    e(halfEdgeIn1).prev = halfEdgeIn3;
+    e(halfEdgeIn2).prev = halfEdgeIn1;
+    e(halfEdgeIn3).prev = halfEdgeIn2;
+    // Finally, create the face, don't forget to set the normal (which should be
+    // normalized)
+    Face tri;
+    tri.edge = halfEdgeIn1;
     mFaces.push_back(tri);
-    // Compute and assign a normal
     size_t faceIndex = mFaces.size() - 1;
     mFaces.back().normal = FaceNormal(faceIndex);
-
-  // All half-edges share the same left face (previously added)
-    e(half_edge_in1).face = faceIndex;
-    e(half_edge_in2).face = faceIndex;
-    e(half_edge_in3).face = faceIndex;
-
-  // Optionally, track the (outer) boundary half-edges
-  // to represent non-closed surfaces
-  return true;
+    // All half-edges share the same left face (previously added)
+    e(halfEdgeIn1).face = faceIndex;
+    e(halfEdgeIn2).face = faceIndex;
+    e(halfEdgeIn3).face = faceIndex;
+    // Optionally, track the (outer) boundary half-edges
+    // to represent non-closed surfaces
+    return true;
 }
 
 /*!
@@ -78,19 +71,18 @@ bool HalfEdgeMesh::AddFace(const std::vector<Vector3<float> > &verts) {
  * inserted (true) or already existed (false)
  */
 bool HalfEdgeMesh::AddVertex(const Vector3<float> &v, size_t &indx) {
-  std::map<Vector3<float>, size_t>::iterator it = mUniqueVerts.find(v);
-  if (it != mUniqueVerts.end()) {
-    indx = (*it).second; // get the index of the already existing vertex
-    return false;
-  }
-
-  mUniqueVerts[v] = indx =
-      GetNumVerts(); // op. [ ] constructs a new entry in map
-  Vertex vert;
-  vert.pos = v;
-  mVerts.push_back(vert); // add it to the vertex list
-
-  return true;
+    std::map<Vector3<float>, size_t>::iterator it = mUniqueVerts.find(v);
+    if (it != mUniqueVerts.end()) {
+        indx = (*it).second;  // get the index of the already existing vertex
+        return false;
+    }
+
+    mUniqueVerts[v] = indx = GetNumVerts();  // op. [ ] constructs a new entry in map
+    Vertex vert;
+    vert.pos = v;
+    mVerts.push_back(vert);  // add it to the vertex list
+
+    return true;
 }
 
 /*!
@@ -102,47 +94,44 @@ bool HalfEdgeMesh::AddVertex(const Vector3<float> &v, size_t &indx) {
  * half-edge from v2 to v1, size_t \return a bool indicating whether the
  * half-edge pair was successfully inserted (true) or already existed (false)
  */
-bool HalfEdgeMesh::AddHalfEdgePair(size_t v1, size_t v2, size_t &indx1,
-                                   size_t &indx2) {
-  std::map<OrderedPair, size_t>::iterator it =
-      mUniqueEdgePairs.find(OrderedPair(v1, v2));
-  if (it != mUniqueEdgePairs.end()) {
-    indx1 = it->second;
-    indx2 = e(it->second).pair;
-    if (v1 != e(indx1).vert) {
-      std::swap(indx1, indx2); // sort correctly
+bool HalfEdgeMesh::AddHalfEdgePair(size_t v1, size_t v2, size_t &indx1, size_t &indx2) {
+    std::map<OrderedPair, size_t>::iterator it = mUniqueEdgePairs.find(OrderedPair(v1, v2));
+    if (it != mUniqueEdgePairs.end()) {
+        indx1 = it->second;
+        indx2 = e(it->second).pair;
+        if (v1 != e(indx1).vert) {
+            std::swap(indx1, indx2);  // sort correctly
+        }
+        return false;
     }
-    return false;
-  }
 
-  // If not found, calculate both half-edges indices
-  indx1 = mEdges.size();
-  indx2 = indx1 + 1;
+    // If not found, calculate both half-edges indices
+    indx1 = mEdges.size();
+    indx2 = indx1 + 1;
 
-  // Create edges and set pair index
-  HalfEdge edge1, edge2;
-  edge1.pair = indx2;
-  edge2.pair = indx1;
+    // Create edges and set pair index
+    HalfEdge edge1, edge2;
+    edge1.pair = indx2;
+    edge2.pair = indx1;
 
-  // Connect the edges to the verts
-  edge1.vert = v1;
-  edge2.vert = v2;
+    // Connect the edges to the verts
+    edge1.vert = v1;
+    edge2.vert = v2;
 
-  // Connect the verts to the edges
-  v(v1).edge = indx1;
-  v(v2).edge = indx2;
+    // Connect the verts to the edges
+    v(v1).edge = indx1;
+    v(v2).edge = indx2;
 
-  // Store the edges in mEdges
-  mEdges.push_back(edge1);
-  mEdges.push_back(edge2);
+    // Store the edges in mEdges
+    mEdges.push_back(edge1);
+    mEdges.push_back(edge2);
 
-  // Store the first edge in the map as an OrderedPair
-  OrderedPair op(v1, v2);
-  mUniqueEdgePairs[op] =
-      indx1; // op. [ ] constructs a new entry in map, ordering not important
-  // sorting done when retrieving
+    // Store the first edge in the map as an OrderedPair
+    OrderedPair op(v1, v2);
+    mUniqueEdgePairs[op] = indx1;  // op. [ ] constructs a new entry in map, ordering not important
+    // sorting done when retrieving
 
-  return true;
+    return true;
 }
 
 /*! \lab1 HalfEdgeMesh Implement the MergeAdjacentBoundaryEdge */
@@ -151,13 +140,13 @@ bool HalfEdgeMesh::AddHalfEdgePair(size_t v1, size_t v2, size_t &indx1,
  * \param [in] indx the index of the INNER half-edge, size_t
  */
 void HalfEdgeMesh::MergeOuterBoundaryEdge(size_t innerEdge) {
-  // Add your code here
-  // 1. Merge first loop (around innerEdge->vert)
-  // 2. Find leftmost edge, last edge counter clock-wise
-  // 3. Test if there's anything to merge
-  // 3a. If so merge the gap
-  // 3b. And set border flags
-  // 4. Merge second loop (around innerEdge->pair->vert)
+    // Add your code here
+    // 1. Merge first loop (around innerEdge->vert)
+    // 2. Find leftmost edge, last edge counter clock-wise
+    // 3. Test if there's anything to merge
+    // 3a. If so merge the gap
+    // 3b. And set border flags
+    // 4. Merge second loop (around innerEdge->pair->vert)
 }
 
 /*! Proceeds to check if the mesh is valid. All indices are inspected and
@@ -166,78 +155,67 @@ void HalfEdgeMesh::MergeOuterBoundaryEdge(size_t innerEdge) {
  * findNeighbourFaces method.
  */
 void HalfEdgeMesh::Validate() {
-  std::vector<HalfEdge>::iterator iterEdge = mEdges.begin();
-  std::vector<HalfEdge>::iterator iterEdgeEnd = mEdges.end();
-  while (iterEdge != iterEdgeEnd) {
-    if ((*iterEdge).face == UNINITIALIZED ||
-        (*iterEdge).next == UNINITIALIZED ||
-        (*iterEdge).pair == UNINITIALIZED ||
-        (*iterEdge).prev == UNINITIALIZED || (*iterEdge).vert == UNINITIALIZED)
-      std::cerr << "HalfEdge " << iterEdge - mEdges.begin()
-                << " not properly initialized" << std::endl;
-
-    iterEdge++;
-  }
-  std::cerr << "Done with edge check (checked " << GetNumEdges() << " edges)"
-            << std::endl;
-
-  std::vector<Face>::iterator iterTri = mFaces.begin();
-  std::vector<Face>::iterator iterTriEnd = mFaces.end();
-  while (iterTri != iterTriEnd) {
-    if ((*iterTri).edge == UNINITIALIZED)
-      std::cerr << "Tri " << iterTri - mFaces.begin()
-                << " not properly initialized" << std::endl;
-
-    iterTri++;
-  }
-  std::cerr << "Done with face check (checked " << GetNumFaces() << " faces)"
-            << std::endl;
-
-  std::vector<Vertex>::iterator iterVertex = mVerts.begin();
-  std::vector<Vertex>::iterator iterVertexEnd = mVerts.end();
-  while (iterVertex != iterVertexEnd) {
-    if ((*iterVertex).edge == UNINITIALIZED)
-      std::cerr << "Vertex " << iterVertex - mVerts.begin()
-                << " not properly initialized" << std::endl;
-
-    iterVertex++;
-  }
-  std::cerr << "Done with vertex check (checked " << GetNumVerts()
-            << " vertices)" << std::endl;
-
-  std::cerr << "Looping through triangle neighborhood of each vertex... ";
-  iterVertex = mVerts.begin();
-  iterVertexEnd = mVerts.end();
-  int emptyCount = 0;
-  std::vector<size_t> problemVerts;
-  while (iterVertex != iterVertexEnd) {
-    std::vector<size_t> foundFaces =
-        FindNeighborFaces(iterVertex - mVerts.begin());
-    std::vector<size_t> foundVerts =
-        FindNeighborVertices(iterVertex - mVerts.begin());
-    if (foundFaces.empty() || foundVerts.empty())
-      emptyCount++;
-    std::set<size_t> uniqueFaces(foundFaces.begin(), foundFaces.end());
-    std::set<size_t> uniqueVerts(foundVerts.begin(), foundVerts.end());
-    if (foundFaces.size() != uniqueFaces.size() ||
-        foundVerts.size() != uniqueVerts.size())
-      problemVerts.push_back(iterVertex - mVerts.begin());
-    iterVertex++;
-  }
-  std::cerr << std::endl
-            << "Done: " << emptyCount << " isolated vertices found"
-            << std::endl;
-  if (problemVerts.size()) {
+    std::vector<HalfEdge>::iterator iterEdge = mEdges.begin();
+    std::vector<HalfEdge>::iterator iterEdgeEnd = mEdges.end();
+    while (iterEdge != iterEdgeEnd) {
+        if ((*iterEdge).face == UNINITIALIZED || (*iterEdge).next == UNINITIALIZED ||
+            (*iterEdge).pair == UNINITIALIZED || (*iterEdge).prev == UNINITIALIZED ||
+            (*iterEdge).vert == UNINITIALIZED)
+            std::cerr << "HalfEdge " << iterEdge - mEdges.begin() << " not properly initialized"
+                      << std::endl;
+
+        iterEdge++;
+    }
+    std::cerr << "Done with edge check (checked " << GetNumEdges() << " edges)" << std::endl;
+
+    std::vector<Face>::iterator iterTri = mFaces.begin();
+    std::vector<Face>::iterator iterTriEnd = mFaces.end();
+    while (iterTri != iterTriEnd) {
+        if ((*iterTri).edge == UNINITIALIZED)
+            std::cerr << "Tri " << iterTri - mFaces.begin() << " not properly initialized"
+                      << std::endl;
+
+        iterTri++;
+    }
+    std::cerr << "Done with face check (checked " << GetNumFaces() << " faces)" << std::endl;
+
+    std::vector<Vertex>::iterator iterVertex = mVerts.begin();
+    std::vector<Vertex>::iterator iterVertexEnd = mVerts.end();
+    while (iterVertex != iterVertexEnd) {
+        if ((*iterVertex).edge == UNINITIALIZED)
+            std::cerr << "Vertex " << iterVertex - mVerts.begin() << " not properly initialized"
+                      << std::endl;
+
+        iterVertex++;
+    }
+    std::cerr << "Done with vertex check (checked " << GetNumVerts() << " vertices)" << std::endl;
+
+    std::cerr << "Looping through triangle neighborhood of each vertex... ";
+    iterVertex = mVerts.begin();
+    iterVertexEnd = mVerts.end();
+    int emptyCount = 0;
+    std::vector<size_t> problemVerts;
+    while (iterVertex != iterVertexEnd) {
+        std::vector<size_t> foundFaces = FindNeighborFaces(iterVertex - mVerts.begin());
+        std::vector<size_t> foundVerts = FindNeighborVertices(iterVertex - mVerts.begin());
+        if (foundFaces.empty() || foundVerts.empty()) emptyCount++;
+        std::set<size_t> uniqueFaces(foundFaces.begin(), foundFaces.end());
+        std::set<size_t> uniqueVerts(foundVerts.begin(), foundVerts.end());
+        if (foundFaces.size() != uniqueFaces.size() || foundVerts.size() != uniqueVerts.size())
+            problemVerts.push_back(iterVertex - mVerts.begin());
+        iterVertex++;
+    }
+    std::cerr << std::endl << "Done: " << emptyCount << " isolated vertices found" << std::endl;
+    if (problemVerts.size()) {
+        std::cerr << std::endl
+                  << "Found " << problemVerts.size() << " duplicate faces in vertices: ";
+        std::copy(problemVerts.begin(), problemVerts.end(),
+                  std::ostream_iterator<size_t>(std::cerr, ", "));
+        std::cerr << "\n";
+    }
     std::cerr << std::endl
-              << "Found " << problemVerts.size()
-              << " duplicate faces in vertices: ";
-    std::copy(problemVerts.begin(), problemVerts.end(),
-              std::ostream_iterator<size_t>(std::cerr, ", "));
-    std::cerr << "\n";
-  }
-  std::cerr << std::endl
-            << "The mesh has genus " << Genus() << ", and consists of "
-            << Shells() << " shells.\n";
+              << "The mesh has genus " << Genus() << ", and consists of " << Shells()
+              << " shells.\n";
 }
 
 /*! \lab1 Implement the FindNeighborVertices */
@@ -245,24 +223,25 @@ void HalfEdgeMesh::Validate() {
  * counter clockwise. \param [in] vertexIndex  the index to vertex, size_t
  * \return a vector containing the indices to all the found vertices.
  */
-std::vector<size_t>
-HalfEdgeMesh::FindNeighborVertices(size_t vertexIndex) const {
-  // Collected vertices, sorted counter clockwise!
-  std::vector<size_t> oneRing;
-
-  size_t start_edge = v(vertexIndex).edge;
-  size_t tmp_edge = start_edge;
-  size_t next_edge;
-  size_t prev_edge;
-
-   do {
-      next_edge = e(tmp_edge).next;
-      prev_edge = e(tmp_edge).prev;
-      oneRing.push_back(e(next_edge).vert);
-      tmp_edge = e(prev_edge).pair;
-  } while (start_edge != tmp_edge);
-
-  return oneRing;
+std::vector<size_t> HalfEdgeMesh::FindNeighborVertices(size_t vertexIndex) const {
+    // Collected vertices, sorted counter clockwise!
+    std::vector<size_t> oneRing;
+
+    size_t start_edge = v(vertexIndex).edge;
+    size_t tmp_edge = start_edge;
+    size_t next_edge;
+    size_t prev_edge;
+
+    do {
+        next_edge = e(tmp_edge).next;
+        prev_edge = e(tmp_edge).prev;
+        oneRing.push_back(e(next_edge).vert);
+        tmp_edge = e(prev_edge).pair;
+    } while (start_edge != tmp_edge);
+
+    // Add your code here
+
+    return oneRing;
 }
 
 /*! \lab1 Implement the FindNeighborFaces */
@@ -271,10 +250,9 @@ HalfEdgeMesh::FindNeighborVertices(size_t vertexIndex) const {
  * \return a vector containing the indices to all the found faces.
  */
 std::vector<size_t> HalfEdgeMesh::FindNeighborFaces(size_t vertexIndex) const {
-  // Collected faces, sorted counter clockwise!
-  std::vector<size_t> foundFaces;
+    // Collected faces, sorted counter clockwise!
+    std::vector<size_t> foundFaces;
 
-  // Add your code here
     size_t start_edge = v(vertexIndex).edge;
     size_t tmp_edge = start_edge;
     size_t prev_edge;
@@ -283,20 +261,46 @@ std::vector<size_t> HalfEdgeMesh::FindNeighborFaces(size_t vertexIndex) const {
         prev_edge = e(tmp_edge).prev;
         foundFaces.push_back(e(tmp_edge).face);
         tmp_edge = e(prev_edge).pair;
-    } while(start_edge != tmp_edge);
-
-        return foundFaces;
+    } while (start_edge != tmp_edge);
+    // Add your code here
+    return foundFaces;
 }
 
 /*! \lab1 Implement the curvature */
 float HalfEdgeMesh::VertexCurvature(size_t vertexIndex) const {
+    //// Copy code from SimpleMesh or compute more accurate estimate
+    //  std::vector<size_t> oneRing = FindNeighborVertices(vertexIndex);
+    //  assert(oneRing.size() != 0);
+
+    //  size_t curr, next;
+    //  const Vector3<float> &vi = mVerts.at(vertexIndex).pos;
+    //  float angleSum = 0;
+    //  float area = 0;
+    //  for (size_t i = 0; i < oneRing.size(); i++) {
+    //      // connections
+    //      curr = oneRing.at(i);
+    //      if (i < oneRing.size() - 1)
+    //          next = oneRing.at(i + 1);
+    //      else
+    //          next = oneRing.front();
+
+    //      // find vertices in 1-ring according to figure 5 in lab text
+    //      // next - beta
+    //      const Vector3<float> &nextPos = mVerts.at(next).pos;
+    //      const Vector3<float> &vj = mVerts.at(curr).pos;
+
+    //      // compute angle and area
+    //      angleSum += acos((vj - vi) * (nextPos - vi) / ((vj - vi).Length() * (nextPos -
+    //      vi).Length())); area += Cross((vi - vj), (nextPos - vj)).Length() * 0.5f;
+    //  }
+    //  return (2.0f * static_cast<float>(M_PI) - angleSum) / area;
 
     size_t vi = vertexIndex;
     std::vector<size_t> vjs = FindNeighborVertices(vertexIndex);
-    //Alpha & Beta
+    // Alpha & Beta
     Vector3<float> v_alpha_vi, v_beta_vi;
     Vector3<float> v_alpha_vj, v_beta_vj;
-    double alpha_j ,beta_j, cot_alpha_plus_beta;
+    double alpha_j, beta_j, cot_alpha_plus_beta;
 
     size_t vj, v_alpha, v_beta;
 
@@ -305,270 +309,204 @@ float HalfEdgeMesh::VertexCurvature(size_t vertexIndex) const {
 
     for (int i = 0; i < vjs.size(); i++) {
 
-        if(i == 0) {
-            v_alpha = vjs[vjs.size()-1];
+        if (i == 0) {
+            v_alpha = vjs[vjs.size() - 1];
             vj = vjs[i];
-            v_beta = vjs[i+1];
-        }
-        else if( i == vjs.size() - 1) {
+            v_beta = vjs[i + 1];
+        } else if (i == vjs.size() - 1) {
             v_alpha = vjs[i - 1];
             vj = vjs[i];
             v_beta = vjs[0];
-        }
-        else {
-            v_alpha = vjs[i-1];
-            vj =  vjs[i];
+        } else {
+            v_alpha = vjs[i - 1];
+            vj = vjs[i];
             v_beta = vjs[i + 1];
         }
 
-        //Getting alpha
+        // Getting alpha
         /* v_alpha_vi = v(vi).pos - v(v_alpha).pos;
         v_alpha_vj = v(vj).pos - v(v_alpha).pos;
         alpha_j = acos((v_alpha_vi * v_alpha_vj) / (v_alpha_vi.Length() * v_alpha_vj.Length())); */
         alpha_j = Cotangent(v(vi).pos, v(v_alpha).pos, v(vj).pos);
 
-        //Getting beta
+        // Getting beta
         /* v_beta_vi = v(vi).pos - v(v_beta).pos;
         v_beta_vj = v(vj).pos - v(v_beta).pos;
         beta_j = acos((v_beta_vi * v_beta_vj) / (v_beta_vi.Length() * v_beta_vj.Length())); */
-        beta_j = Cotangent(v(vj).pos,v(v_beta).pos, v(vi).pos);
+        beta_j = Cotangent(v(vj).pos, v(v_beta).pos, v(vi).pos);
 
-        //cot alpha + cot beta
-        //cot_alpha_plus_beta = (cos(alpha_j) / sin(alpha_j) + cos(beta_j) / sin(beta_j));
+        // cot alpha + cot beta
+        // cot_alpha_plus_beta = (cos(alpha_j) / sin(alpha_j) + cos(beta_j) / sin(beta_j));
         cot_alpha_plus_beta = alpha_j + beta_j;
-        Tj +=  cot_alpha_plus_beta * (v(vi).pos - v(vj).pos);
-        area += cot_alpha_plus_beta * (v(vi).pos - v(vj).pos).Length() * (v(vi).pos - v(vj).pos).Length();
+        Tj += cot_alpha_plus_beta * (v(vi).pos - v(vj).pos);
+        area += cot_alpha_plus_beta * (v(vi).pos - v(vj).pos).Length() *
+                (v(vi).pos - v(vj).pos).Length();
     }
     area /= 8.0;
-    return Tj.Length() * 0.25/area;
-/*
- *  Gaussian Curvature
-    std::vector<size_t> oneRing = FindNeighborVertices(vertexIndex);
-    assert(oneRing.size() != 0);
-
-    size_t curr, next;
-    const Vector3<float> &vi = mVerts.at(vertexIndex).pos;
-    float angleSum = 0;
-    float area = 0;
-    for (size_t i = 0; i < oneRing.size(); i++) {
-        // connections
-        curr = oneRing.at(i);
-        if (i < oneRing.size() - 1)
-            next = oneRing.at(i + 1);
-        else
-            next = oneRing.front();
-
-        // find vertices in 1-ring according to figure 5 in lab text
-        // next - beta
-        const Vector3<float> &nextPos = mVerts.at(next).pos;
-        const Vector3<float> &vj = mVerts.at(curr).pos;
-
-        // compute angle and area
-        angleSum += acos((vj - vi) * (nextPos - vi) /
-                         ((vj - vi).Length() * (nextPos - vi).Length()));
-        area += Cross((vi - vj), (nextPos - vj)).Length() * 0.5f;
-    }
-    return (2.0f * static_cast<float>(M_PI) - angleSum) / area;
-    */
+    return Tj.Length() * 0.25 / area;
 }
 
 float HalfEdgeMesh::FaceCurvature(size_t faceIndex) const {
-  // NB Assumes vertex curvature already computed
-  size_t indx = f(faceIndex).edge;
-  const EdgeIterator it = GetEdgeIterator(indx);
+    // NB Assumes vertex curvature already computed
+    size_t indx = f(faceIndex).edge;
+    const EdgeIterator it = GetEdgeIterator(indx);
 
-  const Vertex &v1 = v(it.GetEdgeVertexIndex());
-  const Vertex &v2 = v(it.Next().GetEdgeVertexIndex());
-  const Vertex &v3 = v(it.Next().GetEdgeVertexIndex());
+    const Vertex &v1 = v(it.GetEdgeVertexIndex());
+    const Vertex &v2 = v(it.Next().GetEdgeVertexIndex());
+    const Vertex &v3 = v(it.Next().GetEdgeVertexIndex());
 
-  return (v1.curvature + v2.curvature + v3.curvature) / 3.f;
+    return (v1.curvature + v2.curvature + v3.curvature) / 3.f;
 }
 
 Vector3<float> HalfEdgeMesh::FaceNormal(size_t faceIndex) const {
-  size_t indx = f(faceIndex).edge;
-  const EdgeIterator it = GetEdgeIterator(indx);
+    size_t indx = f(faceIndex).edge;
+    const EdgeIterator it = GetEdgeIterator(indx);
 
-  const Vector3<float> &p1 = v(it.GetEdgeVertexIndex()).pos;
-  const Vector3<float> &p2 = v(it.Next().GetEdgeVertexIndex()).pos;
-  const Vector3<float> &p3 = v(it.Next().GetEdgeVertexIndex()).pos;
+    const Vector3<float> &p1 = v(it.GetEdgeVertexIndex()).pos;
+    const Vector3<float> &p2 = v(it.Next().GetEdgeVertexIndex()).pos;
+    const Vector3<float> &p3 = v(it.Next().GetEdgeVertexIndex()).pos;
 
-  const Vector3<float> e1 = p2 - p1;
-  const Vector3<float> e2 = p3 - p1;
-  return Cross(e1, e2).Normalize();
+    const Vector3<float> e1 = p2 - p1;
+    const Vector3<float> e2 = p3 - p1;
+    return Cross(e1, e2).Normalize();
 }
 
 Vector3<float> HalfEdgeMesh::VertexNormal(size_t vertexIndex) const {
 
-  Vector3<float> n(0, 0, 0);
-  std::vector<size_t> neighbourFaces = FindNeighborFaces(vertexIndex);
+    Vector3<float> n(0, 0, 0);
+
+    std::vector<size_t> neighbourFaces = FindNeighborFaces(vertexIndex);
 
-  for(size_t face : neighbourFaces) {
-      n += f(face).normal;
-  }
+    for (size_t face : neighbourFaces) {
+        n += f(face).normal;
+    }
 
-  n /= n.Length();
+    n /= n.Length();
 
-  return n;
+    return n;
 }
 
 void HalfEdgeMesh::Initialize() {
-  Validate();
-  Update();
+    Validate();
+    Update();
 }
 
 void HalfEdgeMesh::Update() {
-  // Calculate and store all differentials and area
-
-  // First update all face normals and triangle areas
-  for (size_t i = 0; i < GetNumFaces(); i++) {
-    f(i).normal = FaceNormal(i);
-  }
-  // Then update all vertex normals and curvature
-  for (size_t i = 0; i < GetNumVerts(); i++) {
-    // Vertex normals are just weighted averages
-    mVerts.at(i).normal = VertexNormal(i);
-  }
-
-  // Then update vertex curvature
-  for (size_t i = 0; i < GetNumVerts(); i++) {
-    mVerts.at(i).curvature = VertexCurvature(i);
-    //    std::cerr <<   mVerts.at(i).curvature << "\n";
-  }
-
-  // Finally update face curvature
-  for (size_t i = 0; i < GetNumFaces(); i++) {
-    f(i).curvature = FaceCurvature(i);
-  }
-
-  std::cerr << "Area: " << Area() << ".\n";
-  std::cerr << "Volume: " << Volume() << ".\n";
-
-  // Update vertex and face colors
-  if (mVisualizationMode == CurvatureVertex) {
-    std::vector<Vertex>::iterator iter = mVerts.begin();
-    std::vector<Vertex>::iterator iend = mVerts.end();
-    float minCurvature = (std::numeric_limits<float>::max)();
-    float maxCurvature = -(std::numeric_limits<float>::max)();
-    while (iter != iend) {
-      if (minCurvature > (*iter).curvature)
-        minCurvature = (*iter).curvature;
-      if (maxCurvature < (*iter).curvature)
-        maxCurvature = (*iter).curvature;
-      iter++;
+    // Calculate and store all differentials and area
+
+    // First update all face normals and triangle areas
+    for (size_t i = 0; i < GetNumFaces(); i++) {
+        f(i).normal = FaceNormal(i);
     }
-    std::cerr << "Mapping color based on vertex curvature with range ["
-              << minCurvature << "," << maxCurvature << "]" << std::endl;
-    iter = mVerts.begin();
-    while (iter != iend) {
-      (*iter).color =
-          mColorMap->Map((*iter).curvature, minCurvature, maxCurvature);
-      iter++;
+    // Then update all vertex normals and curvature
+    for (size_t i = 0; i < GetNumVerts(); i++) {
+        // Vertex normals are just weighted averages
+        mVerts.at(i).normal = VertexNormal(i);
     }
-  } else if (mVisualizationMode == CurvatureFace) {
-    std::vector<Face>::iterator iter = mFaces.begin();
-    std::vector<Face>::iterator iend = mFaces.end();
-    float minCurvature = (std::numeric_limits<float>::max)();
-    float maxCurvature = -(std::numeric_limits<float>::max)();
-    while (iter != iend) {
-      if (minCurvature > (*iter).curvature)
-        minCurvature = (*iter).curvature;
-      if (maxCurvature < (*iter).curvature)
-        maxCurvature = (*iter).curvature;
-      iter++;
+
+    // Then update vertex curvature
+    for (size_t i = 0; i < GetNumVerts(); i++) {
+        mVerts.at(i).curvature = VertexCurvature(i);
+        //    std::cerr <<   mVerts.at(i).curvature << "\n";
     }
-    std::cerr << "Mapping color based on face curvature with range ["
-              << minCurvature << "," << maxCurvature << "]" << std::endl;
-    iter = mFaces.begin();
-    while (iter != iend) {
-      (*iter).color =
-          mColorMap->Map((*iter).curvature, minCurvature, maxCurvature);
-      iter++;
+
+    // Finally update face curvature
+    for (size_t i = 0; i < GetNumFaces(); i++) {
+        f(i).curvature = FaceCurvature(i);
+    }
+
+    std::cerr << "Area: " << Area() << ".\n";
+    std::cerr << "Volume: " << Volume() << ".\n";
+
+    // Update vertex and face colors
+    if (mVisualizationMode == CurvatureVertex) {
+        std::vector<Vertex>::iterator iter = mVerts.begin();
+        std::vector<Vertex>::iterator iend = mVerts.end();
+        float minCurvature = (std::numeric_limits<float>::max)();
+        float maxCurvature = -(std::numeric_limits<float>::max)();
+        while (iter != iend) {
+            if (minCurvature > (*iter).curvature) minCurvature = (*iter).curvature;
+            if (maxCurvature < (*iter).curvature) maxCurvature = (*iter).curvature;
+            iter++;
+        }
+        std::cerr << "Mapping color based on vertex curvature with range [" << minCurvature << ","
+                  << maxCurvature << "]" << std::endl;
+        iter = mVerts.begin();
+        while (iter != iend) {
+            (*iter).color = mColorMap->Map((*iter).curvature, minCurvature, maxCurvature);
+            iter++;
+        }
+    } else if (mVisualizationMode == CurvatureFace) {
+        std::vector<Face>::iterator iter = mFaces.begin();
+        std::vector<Face>::iterator iend = mFaces.end();
+        float minCurvature = (std::numeric_limits<float>::max)();
+        float maxCurvature = -(std::numeric_limits<float>::max)();
+        while (iter != iend) {
+            if (minCurvature > (*iter).curvature) minCurvature = (*iter).curvature;
+            if (maxCurvature < (*iter).curvature) maxCurvature = (*iter).curvature;
+            iter++;
+        }
+        std::cerr << "Mapping color based on face curvature with range [" << minCurvature << ","
+                  << maxCurvature << "]" << std::endl;
+        iter = mFaces.begin();
+        while (iter != iend) {
+            (*iter).color = mColorMap->Map((*iter).curvature, minCurvature, maxCurvature);
+            iter++;
+        }
     }
-  }
 }
 
 /*! \lab1 Implement the area */
 float HalfEdgeMesh::Area() const {
-  float area = 0;
-  size_t v1;
-  size_t v2;
-  size_t v3;
-
-  for(Face face : mFaces) {
-      v1 = e(face.edge).vert;
-      v2 = e(e(face.edge).next).vert;
-      v3 = e(e(face.edge).prev).vert;
-      Vector3<float> crossProduct = Cross((v(v2).pos - v(v1).pos), v(v3).pos - v(v1).pos);
-      area += crossProduct.Length() * 0.5;
-  }
-
-  return area;
+    float area = 0;
+    size_t v1;
+    size_t v2;
+    size_t v3;
+
+    for (Face face : mFaces) {
+        v1 = e(face.edge).vert;
+        v2 = e(e(face.edge).next).vert;
+        v3 = e(e(face.edge).prev).vert;
+        Vector3<float> crossProduct = Cross((v(v2).pos - v(v1).pos), v(v3).pos - v(v1).pos);
+        area += crossProduct.Length() * 0.5;
+    }
+
+    return area;
 }
 
 /*! \lab1 Implement the volume */
 float HalfEdgeMesh::Volume() const {
-  float volume = 0;
-  float area{};
-  size_t v1;
-  size_t v2;
-  size_t v3;
-
-    for(Face face : mFaces) {
-        //Vertices
+    float volume = 0;
+    // Add code here
+    float area{};
+    size_t v1;
+    size_t v2;
+    size_t v3;
+
+    for (Face face : mFaces) {
+        // Vertices
         v1 = e(face.edge).vert;
         v2 = e(e(face.edge).next).vert;
         v3 = e(e(face.edge).prev).vert;
 
-        //Centroid vertex
+        // Centroid vertex
         Vector3<float> centroid = (v(v1).pos + v(v2).pos + v(v3).pos) * 0.333333;
 
-        //Area for a face
+        // Area for a face
         Vector3<float> crossProduct = Cross((v(v2).pos - v(v1).pos), v(v3).pos - v(v1).pos);
         area = crossProduct.Length() * 0.5;
 
         volume += (centroid * face.normal) * area;
     }
 
-  return volume * 0.333333;
+    return volume * 0.333333;
 }
 
 /*! \lab1 Calculate the number of shells  */
-int HalfEdgeMesh::Shells() const {
-/*
-    std::set<size_t> vertexQueueSet;
-    std::map<size_t, int> vertexTaggedSet;
-    size_t v_ = 0;
-    int shells{};
-
-    vertexQueueSet.insert(v_);
-    while(vertexQueueSet.size() != 0) {
-        vertexQueueSet.erase(v_);
-        vertexTaggedSet.insert(v_);
-        for(size_t vi : FindNeighborVertices(v_)) {
-            if(vertexTaggedSet.find(vi) == vertexTaggedSet.end() ) {
-                vertexQueueSet.insert(vi);
-            }
-        }
-    }
-*/
-
-
-   /* push first vertex → vertexQueueSet
-    while vertexQueueSet ̸= empty do
-        v ← vertexQueueSet.pop
-    push v → vertexTaggedSet
-    for all vi in findNeighborVertices(v) do
-        if vi ∈/ vertexTaggedSet then push vi → vertexQueueSet
-    end if end for
-    end while
-    */
-
-    return 1;
-}
+int HalfEdgeMesh::Shells() const { return 1; }
 
 /*! \lab1 Implement the genus */
 size_t HalfEdgeMesh::Genus() const {
-
     // Add code here
     size_t genus;
     genus = -1 * ((static_cast<int>(mVerts.size()) - static_cast<int>(mEdges.size()) + static_cast<int>(mFaces.size()) - 2) * 0.5);
@@ -577,138 +515,136 @@ size_t HalfEdgeMesh::Genus() const {
 }
 
 void HalfEdgeMesh::Dilate(float amount) {
-  std::vector<Vertex>::iterator iter = mVerts.begin();
-  std::vector<Vertex>::iterator iend = mVerts.end();
-  while (iter != iend) {
-    (*iter).pos += amount * (*iter).normal;
-    iter++;
-  }
-
-  Initialize();
-  Update();
+    std::vector<Vertex>::iterator iter = mVerts.begin();
+    std::vector<Vertex>::iterator iend = mVerts.end();
+    while (iter != iend) {
+        (*iter).pos += amount * (*iter).normal;
+        iter++;
+    }
+
+    Initialize();
+    Update();
 }
 
 void HalfEdgeMesh::Erode(float amount) {
-  std::vector<Vertex>::iterator iter = mVerts.begin();
-  std::vector<Vertex>::iterator iend = mVerts.end();
-  while (iter != iend) {
-    (*iter).pos -= amount * (*iter).normal;
-    iter++;
-  }
-
-  Initialize();
-  Update();
+    std::vector<Vertex>::iterator iter = mVerts.begin();
+    std::vector<Vertex>::iterator iend = mVerts.end();
+    while (iter != iend) {
+        (*iter).pos -= amount * (*iter).normal;
+        iter++;
+    }
+
+    Initialize();
+    Update();
 }
 
 void HalfEdgeMesh::Smooth(float amount) {
-  std::vector<Vertex>::iterator iter = mVerts.begin();
-  std::vector<Vertex>::iterator iend = mVerts.end();
-  while (iter != iend) {
-    (*iter).pos -= amount * (*iter).normal * (*iter).curvature;
-    iter++;
-  }
-
-  Initialize();
-  Update();
+    std::vector<Vertex>::iterator iter = mVerts.begin();
+    std::vector<Vertex>::iterator iend = mVerts.end();
+    while (iter != iend) {
+        (*iter).pos -= amount * (*iter).normal * (*iter).curvature;
+        iter++;
+    }
+
+    Initialize();
+    Update();
 }
 
 void HalfEdgeMesh::Render() {
-  glEnable(GL_LIGHTING);
-  glMatrixMode(GL_MODELVIEW);
+    glEnable(GL_LIGHTING);
+    glMatrixMode(GL_MODELVIEW);
 
-  // Apply transform
-  glPushMatrix(); // Push modelview matrix onto stack
+    // Apply transform
+    glPushMatrix();  // Push modelview matrix onto stack
 
-  // Convert transform-matrix to format matching GL matrix format
-  // Load transform into modelview matrix
-  glMultMatrixf(mTransform.ToGLMatrix().GetArrayPtr());
+    // Convert transform-matrix to format matching GL matrix format
+    // Load transform into modelview matrix
+    glMultMatrixf(mTransform.ToGLMatrix().GetArrayPtr());
 
-  if (mWireframe)
-    glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
+    if (mWireframe) glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
 
-  // Draw geometry
-  glBegin(GL_TRIANGLES);
-  const auto numTriangles = GetNumFaces();
-  for (size_t i = 0; i < numTriangles; i++) {
+    // Draw geometry
+    glBegin(GL_TRIANGLES);
+    const auto numTriangles = GetNumFaces();
+    for (size_t i = 0; i < numTriangles; i++) {
 
-    Face &face = f(i);
+        Face &face = f(i);
 
-    HalfEdge *edge = &e(face.edge);
+        HalfEdge *edge = &e(face.edge);
 
-    Vertex &v1 = v(edge->vert);
-    edge = &e(edge->next);
+        Vertex &v1 = v(edge->vert);
+        edge = &e(edge->next);
 
-    Vertex &v2 = v(edge->vert);
-    edge = &e(edge->next);
+        Vertex &v2 = v(edge->vert);
+        edge = &e(edge->next);
 
-    Vertex &v3 = v(edge->vert);
+        Vertex &v3 = v(edge->vert);
 
-    if (mVisualizationMode == CurvatureVertex) {
-      glColor3fv(v1.color.GetArrayPtr());
-      glNormal3fv(v1.normal.GetArrayPtr());
-      glVertex3fv(v1.pos.GetArrayPtr());
-
-      glColor3fv(v2.color.GetArrayPtr());
-      glNormal3fv(v2.normal.GetArrayPtr());
-      glVertex3fv(v2.pos.GetArrayPtr());
-
-      glColor3fv(v3.color.GetArrayPtr());
-      glNormal3fv(v3.normal.GetArrayPtr());
-      glVertex3fv(v3.pos.GetArrayPtr());
-    } else {
-      glColor3fv(face.color.GetArrayPtr());
-      glNormal3fv(face.normal.GetArrayPtr());
-
-      glVertex3fv(v1.pos.GetArrayPtr());
-      glVertex3fv(v2.pos.GetArrayPtr());
-      glVertex3fv(v3.pos.GetArrayPtr());
+        if (mVisualizationMode == CurvatureVertex) {
+            glColor3fv(v1.color.GetArrayPtr());
+            glNormal3fv(v1.normal.GetArrayPtr());
+            glVertex3fv(v1.pos.GetArrayPtr());
+
+            glColor3fv(v2.color.GetArrayPtr());
+            glNormal3fv(v2.normal.GetArrayPtr());
+            glVertex3fv(v2.pos.GetArrayPtr());
+
+            glColor3fv(v3.color.GetArrayPtr());
+            glNormal3fv(v3.normal.GetArrayPtr());
+            glVertex3fv(v3.pos.GetArrayPtr());
+        } else {
+            glColor3fv(face.color.GetArrayPtr());
+            glNormal3fv(face.normal.GetArrayPtr());
+
+            glVertex3fv(v1.pos.GetArrayPtr());
+            glVertex3fv(v2.pos.GetArrayPtr());
+            glVertex3fv(v3.pos.GetArrayPtr());
+        }
     }
-  }
-  glEnd();
+    glEnd();
 
-  // Mesh normals by courtesy of Richard Khoury
-  if (mShowNormals) {
-    glDisable(GL_LIGHTING);
-    glBegin(GL_LINES);
-    const auto numTriangles = GetNumFaces();
-    for (size_t i = 0; i < numTriangles; i++) {
+    // Mesh normals by courtesy of Richard Khoury
+    if (mShowNormals) {
+        glDisable(GL_LIGHTING);
+        glBegin(GL_LINES);
+        const auto numTriangles = GetNumFaces();
+        for (size_t i = 0; i < numTriangles; i++) {
 
-      Face &face = f(i);
+            Face &face = f(i);
 
-      HalfEdge *edge = &e(face.edge);
+            HalfEdge *edge = &e(face.edge);
 
-      Vertex &v1 = v(edge->vert);
-      edge = &e(edge->next);
+            Vertex &v1 = v(edge->vert);
+            edge = &e(edge->next);
 
-      Vertex &v2 = v(edge->vert);
-      edge = &e(edge->next);
+            Vertex &v2 = v(edge->vert);
+            edge = &e(edge->next);
 
-      Vertex &v3 = v(edge->vert);
+            Vertex &v3 = v(edge->vert);
 
-      Vector3<float> faceStart = (v1.pos + v2.pos + v3.pos) / 3.0f;
-      Vector3<float> faceEnd = faceStart + face.normal * 0.1f;
+            Vector3<float> faceStart = (v1.pos + v2.pos + v3.pos) / 3.0f;
+            Vector3<float> faceEnd = faceStart + face.normal * 0.1f;
 
-      glColor3f(1, 0, 0); // Red for face normal
-      glVertex3fv(faceStart.GetArrayPtr());
-      glVertex3fv(faceEnd.GetArrayPtr());
+            glColor3f(1, 0, 0);  // Red for face normal
+            glVertex3fv(faceStart.GetArrayPtr());
+            glVertex3fv(faceEnd.GetArrayPtr());
 
-      glColor3f(0, 1, 0); // Vertex normals in Green
-      glVertex3fv(v1.pos.GetArrayPtr());
-      glVertex3fv((v1.pos + v1.normal * 0.1f).GetArrayPtr());
-      glVertex3fv(v2.pos.GetArrayPtr());
-      glVertex3fv((v2.pos + v2.normal * 0.1f).GetArrayPtr());
-      glVertex3fv(v3.pos.GetArrayPtr());
-      glVertex3fv((v3.pos + v3.normal * 0.1f).GetArrayPtr());
+            glColor3f(0, 1, 0);  // Vertex normals in Green
+            glVertex3fv(v1.pos.GetArrayPtr());
+            glVertex3fv((v1.pos + v1.normal * 0.1f).GetArrayPtr());
+            glVertex3fv(v2.pos.GetArrayPtr());
+            glVertex3fv((v2.pos + v2.normal * 0.1f).GetArrayPtr());
+            glVertex3fv(v3.pos.GetArrayPtr());
+            glVertex3fv((v3.pos + v3.normal * 0.1f).GetArrayPtr());
+        }
+        glEnd();
+        glEnable(GL_LIGHTING);
     }
-    glEnd();
-    glEnable(GL_LIGHTING);
-  }
 
-  if (mWireframe)
-    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+    if (mWireframe) glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
 
-  // Restore modelview matrix
-  glPopMatrix();
+    // Restore modelview matrix
+    glPopMatrix();
 
-  GLObject::Render();
+    GLObject::Render();
 }
diff --git a/Geometry/HalfEdgeMesh.h b/Geometry/HalfEdgeMesh.h
index f10244a6d6c1be40403810eaca10244b55b8c7fb..fbbd2b4096bfa6879ad59a92fe2978533f652f9d 100644
--- a/Geometry/HalfEdgeMesh.h
+++ b/Geometry/HalfEdgeMesh.h
@@ -148,6 +148,7 @@ protected:
   struct Face : public Mesh::Face {
     Face() : edge(UNINITIALIZED) {}
     size_t edge; //!< index into mEdges
+    
   };
 
   //! The edges of the mesh
diff --git a/Geometry/Quadric.cpp b/Geometry/Quadric.cpp
index 85f10dbd77bf86ae2c291ba54f9a355d2571e253..cdfa350c9f9aac66c4e7989ec8e66f4bb53c7ff6 100644
--- a/Geometry/Quadric.cpp
+++ b/Geometry/Quadric.cpp
@@ -20,11 +20,20 @@ Quadric::~Quadric() {}
  * of the world-coordinates by mWorld2Obj, or transformation of the quadric
  * coefficient matrix by GetTransform() ONCE (see Section 2.2 in lab text).
  */
-float Quadric::GetValue(float x, float y, float z) const { return 0; }
+float Quadric::GetValue(float x, float y, float z) const { 
+    TransformW2O(x,y,z);
+    Vector4<float> p = {x, y, z, 1.0f};
+    
+    return p*(mQuadric*p);
+}
 
 /*!
  * Use the quadric matrix to evaluate the gradient.
  */
 Vector3<float> Quadric::GetGradient(float x, float y, float z) const {
-  return Vector3<float>(0, 0, 0);
+  TransformW2O(x,y,z);
+  Vector4<float> p = {x, y, z, 1.0f};
+  Vector4<float> vec = mQuadric*p;
+  Vector3<float> grad =  {vec[0], vec[1], vec[2]};
+  return 2.0*grad;
 }
diff --git a/Math/Vector4.h b/Math/Vector4.h
index ae0f1ca92ff955d80cf5bdbab2c01c84cdf42bee..72baef7fae2febfc310bcad634689fa51d1ebcba 100644
--- a/Math/Vector4.h
+++ b/Math/Vector4.h
@@ -6,7 +6,7 @@
  *   Ken Museth (ken.museth@itn.liu.se)
  *   Michael Bang Nielsen (bang@daimi.au.dk)
  *   Ola Nilsson (ola.nilsson@itn.liu.se)
- *   Andreas S�derstr�m (andreas.soderstrom@itn.liu.se)
+ *   Andreas S�derstr�m (andreas.soderstrom@itn.liu.se)
  *
  *************************************************************************************************/
 #ifndef _vector4_h
@@ -60,9 +60,9 @@ public:
                    v[3] * vec4.v[3]);
   }
 
-  template<typename Real2>
-  friend std::ostream& std::operator << (std::ostream &os, const
-   Vector4<Real2> &m);
+  // template<typename Real2>
+  // friend std::ostream& std::operator << (std::ostream &os, const
+  // Vector4<Real2> &m);
 
 protected:
   Real v[4];
diff --git a/README.md b/README.md
index 7c8d875b43b96d30d094d2c222bcb22421de1f71..ecec9df49c194ae9b5e89f05eb45ed2142868890 100644
--- a/README.md
+++ b/README.md
@@ -72,6 +72,7 @@ wxWidgets Windows libraries compiled from git branch ```WX_2_8_BRANCH``` using t
 and then choosing one of the color maps.
 - In the output log, it should say that 0 isolated vertices were found.
 - If you want to compute the genus, remember that in the halfedge data structure, there are two halfedges stored for every real edge.
+- When calculating the the genus, remember that `size_t` (e.g., returned by `GetNumVerts()`) is an unsigned type. This is problematic when intermediate results are negative (such as `V-E`). Make sure to use a signed type instead.
 
 ### Lab 2
 
diff --git a/Subdivision/LoopSubdivisionMesh.cpp b/Subdivision/LoopSubdivisionMesh.cpp
index 4ebd1ccd57e9d215dfe57c0bbb81e17dd0f56aec..c38dd9f909ea720eadd6a58ae2c93fa6c23294f7 100644
--- a/Subdivision/LoopSubdivisionMesh.cpp
+++ b/Subdivision/LoopSubdivisionMesh.cpp
@@ -16,121 +16,149 @@
 /*! Subdivides the mesh uniformly one step
  */
 void LoopSubdivisionMesh::Subdivide() {
-  // Create new mesh and copy all the attributes
-  HalfEdgeMesh subDivMesh;
-  subDivMesh.SetTransform(GetTransform());
-  subDivMesh.SetName(GetName());
-  subDivMesh.SetColorMap(GetColorMap());
-  subDivMesh.SetWireframe(GetWireframe());
-  subDivMesh.SetShowNormals(GetShowNormals());
-  subDivMesh.SetOpacity(GetOpacity());
-  if (IsHovering())
-    subDivMesh.Hover();
-  if (IsSelected())
-    subDivMesh.Select();
-  subDivMesh.mMinCMap = mMinCMap;
-  subDivMesh.mMaxCMap = mMaxCMap;
-  subDivMesh.mAutoMinMax = mAutoMinMax;
-
-  // loop over each face and create 4 new ones
-  for (size_t i = 0; i < mFaces.size(); i++) {
-    // subdivide face
-    std::vector<std::vector<Vector3<float>>> faces = Subdivide(i);
-
-    // add new faces to subDivMesh
-    for (size_t j = 0; j < faces.size(); j++) {
-      subDivMesh.AddFace(faces.at(j));
+    // Create new mesh and copy all the attributes
+    HalfEdgeMesh subDivMesh;
+    subDivMesh.SetTransform(GetTransform());
+    subDivMesh.SetName(GetName());
+    subDivMesh.SetColorMap(GetColorMap());
+    subDivMesh.SetWireframe(GetWireframe());
+    subDivMesh.SetShowNormals(GetShowNormals());
+    subDivMesh.SetOpacity(GetOpacity());
+    if (IsHovering()) subDivMesh.Hover();
+    if (IsSelected()) subDivMesh.Select();
+    subDivMesh.mMinCMap = mMinCMap;
+    subDivMesh.mMaxCMap = mMaxCMap;
+    subDivMesh.mAutoMinMax = mAutoMinMax;
+
+    // loop over each face and create 4 new ones
+    for (size_t i = 0; i < mFaces.size(); i++) {
+        // subdivide face
+        std::vector<std::vector<Vector3<float>>> faces = Subdivide(i);
+
+        // add new faces to subDivMesh
+        for (size_t j = 0; j < faces.size(); j++) {
+            subDivMesh.AddFace(faces.at(j));
+        }
     }
-  }
 
-  // Assigns the new mesh
-  *this = LoopSubdivisionMesh(subDivMesh, ++mNumSubDivs);
-  Update();
+    // Assigns the new mesh
+    *this = LoopSubdivisionMesh(subDivMesh, ++mNumSubDivs);
+    Update();
 }
 
 /*! Subdivides the face at faceindex into a vector of faces
  */
-std::vector<std::vector<Vector3<float>>>
-LoopSubdivisionMesh::Subdivide(size_t faceIndex) {
-  std::vector<std::vector<Vector3<float>>> faces;
-  EdgeIterator eit = GetEdgeIterator(f(faceIndex).edge);
-
-  // get the inner halfedges
-  size_t e0, e1, e2;
-  // and their vertex indices
-  size_t v0, v1, v2;
-
-  e0 = eit.GetEdgeIndex();
-  v0 = eit.GetEdgeVertexIndex();
-  eit.Next();
-  e1 = eit.GetEdgeIndex();
-  v1 = eit.GetEdgeVertexIndex();
-  eit.Next();
-  e2 = eit.GetEdgeIndex();
-  v2 = eit.GetEdgeVertexIndex();
-
-  // Compute positions of the vertices
-  Vector3<float> pn0 = VertexRule(v0);
-  Vector3<float> pn1 = VertexRule(v1);
-  Vector3<float> pn2 = VertexRule(v2);
-
-  // Compute positions of the edge vertices
-  Vector3<float> pn3 = EdgeRule(e0);
-  Vector3<float> pn4 = EdgeRule(e1);
-  Vector3<float> pn5 = EdgeRule(e2);
-
-  // add the four new triangles to new mesh
-  std::vector<Vector3<float>> verts;
-  verts.push_back(pn0);
-  verts.push_back(pn3);
-  verts.push_back(pn5);
-  faces.push_back(verts);
-  verts.clear();
-  verts.push_back(pn3);
-  verts.push_back(pn4);
-  verts.push_back(pn5);
-  faces.push_back(verts);
-  verts.clear();
-  verts.push_back(pn3);
-  verts.push_back(pn1);
-  verts.push_back(pn4);
-  faces.push_back(verts);
-  verts.clear();
-  verts.push_back(pn5);
-  verts.push_back(pn4);
-  verts.push_back(pn2);
-  faces.push_back(verts);
-  return faces;
+std::vector<std::vector<Vector3<float>>> LoopSubdivisionMesh::Subdivide(size_t faceIndex) {
+    std::vector<std::vector<Vector3<float>>> faces;
+    EdgeIterator eit = GetEdgeIterator(f(faceIndex).edge);
+
+    // get the inner halfedges
+    size_t e0, e1, e2;
+    // and their vertex indices
+    size_t v0, v1, v2;
+
+    e0 = eit.GetEdgeIndex();
+    v0 = eit.GetEdgeVertexIndex();
+    eit.Next();
+    e1 = eit.GetEdgeIndex();
+    v1 = eit.GetEdgeVertexIndex();
+    eit.Next();
+    e2 = eit.GetEdgeIndex();
+    v2 = eit.GetEdgeVertexIndex();
+
+    // Compute positions of the vertices
+    Vector3<float> pn0 = VertexRule(v0);
+    Vector3<float> pn1 = VertexRule(v1);
+    Vector3<float> pn2 = VertexRule(v2);
+
+    // Compute positions of the edge vertices
+    Vector3<float> pn3 = EdgeRule(e0);
+    Vector3<float> pn4 = EdgeRule(e1);
+    Vector3<float> pn5 = EdgeRule(e2);
+
+    // add the four new triangles to new mesh
+    std::vector<Vector3<float>> verts;
+    verts.push_back(pn0);
+    verts.push_back(pn3);
+    verts.push_back(pn5);
+    faces.push_back(verts);
+    verts.clear();
+    verts.push_back(pn3);
+    verts.push_back(pn4);
+    verts.push_back(pn5);
+    faces.push_back(verts);
+    verts.clear();
+    verts.push_back(pn3);
+    verts.push_back(pn1);
+    verts.push_back(pn4);
+    faces.push_back(verts);
+    verts.clear();
+    verts.push_back(pn5);
+    verts.push_back(pn4);
+    verts.push_back(pn2);
+    faces.push_back(verts);
+    return faces;
 }
 
 /*! Computes a new vertex, replacing a vertex in the old mesh
  */
 Vector3<float> LoopSubdivisionMesh::VertexRule(size_t vertexIndex) {
-  // Get the current vertex
-  Vector3<float> vtx = v(vertexIndex).pos;
+    // Get the current vertex
+    Vector3<float> vtx = v(vertexIndex).pos;
 
-  return vtx;
+    int valence = 0;
+    double weight = 0.5f;
+    std::vector<size_t> n_verts = FindNeighborVertices(vertexIndex);
+    Vector3<float> sum = Vector3<float>(0,0,0);
+
+    for (size_t n : n_verts) {
+        sum += v(n).pos;
+    }
+
+    valence = n_verts.size();
+    weight = Beta(valence);
+
+    //std::cout << weight << std::endl;
+
+    vtx = vtx * (1 - valence * weight) + (sum*weight);
+
+    return vtx;
 }
 
 /*! Computes a new vertex, placed along an edge in the old mesh
  */
 Vector3<float> LoopSubdivisionMesh::EdgeRule(size_t edgeIndex) {
 
-  // Place the edge vertex halfway along the edge
-  HalfEdge &e0 = e(edgeIndex);
-  HalfEdge &e1 = e(e0.pair);
-  Vector3<float> &v0 = v(e0.vert).pos;
-  Vector3<float> &v1 = v(e1.vert).pos;
-  return (v0 + v1) * 0.5f;
+    // Place the edge vertex halfway along the edge
+    HalfEdge &e0 = e(edgeIndex);
+    HalfEdge &e1 = e(e0.pair);
+    Vector3<float> &v0 = v(e0.vert).pos;
+    Vector3<float> &v1 = v(e1.vert).pos;
+
+    // size_t temp = e(edgeIndex).pair;
+    // int valence = 0;
+    HalfEdge &e2 = e(e0.prev);
+    HalfEdge &e3 = e(e1.prev);
+    Vector3<float> &v2 = v(e2.vert).pos;
+    Vector3<float> &v3 = v(e3.vert).pos;
+
+    double weight = 0.5f;
+
+    // valence = FindNeighborVertices(e0.vert).size();
+    // weight = Beta(valence);
+
+    //std::cout << weight << std::endl;
+    //return (v0 + v1) * weight;
+    return (((3.0 / 8.0) * (v0 + v1)) + ((1.0 / 8.0) * (v2 + v3)));
 }
 
 //! Return weights for interior verts
 float LoopSubdivisionMesh::Beta(size_t valence) {
-  if (valence == 6) {
-    return 1.0f / 16.0f;
-  } else if (valence == 3) {
-    return 3.0f / 16.0f;
-  } else {
-    return 3.0f / (8.0f * valence);
-  }
+    if (valence == 6) {
+        return 1.0f / 16.0f;
+    } else if (valence == 3) {
+        return 3.0f / 16.0f;
+    } else {
+        return 3.0f / (8.0f * valence);
+    }
 }
diff --git a/Subdivision/UniformCubicSplineSubdivisionCurve.cpp b/Subdivision/UniformCubicSplineSubdivisionCurve.cpp
index 1efaa0bf3b559a1e044d20c9c3b0e0cb6d9c928d..405b81fbef1250e426196a84698d38738fc3b87e 100644
--- a/Subdivision/UniformCubicSplineSubdivisionCurve.cpp
+++ b/Subdivision/UniformCubicSplineSubdivisionCurve.cpp
@@ -2,56 +2,69 @@
 #include "UniformCubicSplineSubdivisionCurve.h"
 
 UniformCubicSplineSubdivisionCurve::UniformCubicSplineSubdivisionCurve(
-    const std::vector<Vector3<float>> &joints, Vector3<float> lineColor,
-    float lineWidth)
+    const std::vector<Vector3<float>> &joints, Vector3<float> lineColor, float lineWidth)
     : mCoefficients(joints), mControlPolygon(joints) {
-  this->mLineColor = lineColor;
-  this->mLineWidth = lineWidth;
+    this->mLineColor = lineColor;
+    this->mLineWidth = lineWidth;
 }
 
 void UniformCubicSplineSubdivisionCurve::Subdivide() {
-  // Allocate space for new coefficients
-  std::vector<Vector3<float>> newc;
+    // Allocate space for new coefficients
+    std::vector<Vector3<float>> newc;
 
-  assert(mCoefficients.size() > 4 && "Need at least 5 points to subdivide");
+    assert(mCoefficients.size() > 4 && "Need at least 5 points to subdivide");
 
-  // Implement the subdivision scheme for a natural cubic spline here
+    // Implement the subdivision scheme for a natural cubic spline here
+    Vector3<float> updated_coeff;
+    Vector3<float> subdiv_coeff;
 
+    newc.push_back(mCoefficients[0]);
+    newc.push_back(0.125 * ((4 * mCoefficients[0]) + (4*mCoefficients[1])));
 
-  // If 'mCoefficients' had size N, how large should 'newc' be? Perform a check
-  // here!
-  assert(true && "Incorrect number of new coefficients!");
 
-  mCoefficients = newc;
+    for (int i = 1; i < mCoefficients.size() - 1; ++i) {
+        updated_coeff = 0.125 * (mCoefficients[i - 1] + (6 * mCoefficients[i]) + mCoefficients[i + 1]);
+        subdiv_coeff = 0.125 * ((4 * mCoefficients[i]) + (4*mCoefficients[i + 1]));
+
+        newc.push_back(updated_coeff);
+        newc.push_back(subdiv_coeff);
+    }
+
+    newc.push_back(mCoefficients[mCoefficients.size()-1]);
+    // If 'mCoefficients' had size N, how large should 'newc' be? Perform a check
+    // here!
+    assert(((mCoefficients.size()*2)-1) == newc.size());
+
+    mCoefficients = newc;
 }
 
 void UniformCubicSplineSubdivisionCurve::Render() {
-  // Apply transform
-  glPushMatrix(); // Push modelview matrix onto stack
+    // Apply transform
+    glPushMatrix();  // Push modelview matrix onto stack
 
-  // Convert transform-matrix to format matching GL matrix format
-  // Load transform into modelview matrix
-  glMultMatrixf(mTransform.ToGLMatrix().GetArrayPtr());
+    // Convert transform-matrix to format matching GL matrix format
+    // Load transform into modelview matrix
+    glMultMatrixf(mTransform.ToGLMatrix().GetArrayPtr());
 
-  mControlPolygon.Render();
+    mControlPolygon.Render();
 
-  // save line point and color states
-  glPushAttrib(GL_POINT_BIT | GL_LINE_BIT | GL_CURRENT_BIT);
+    // save line point and color states
+    glPushAttrib(GL_POINT_BIT | GL_LINE_BIT | GL_CURRENT_BIT);
 
-  // draw segments
-  glLineWidth(mLineWidth);
-  glColor3fv(mLineColor.GetArrayPtr());
-  glBegin(GL_LINE_STRIP);
-  // just draw the spline as a series of connected linear segments
-  for (size_t i = 0; i < mCoefficients.size(); i++) {
-    glVertex3fv(mCoefficients.at(i).GetArrayPtr());
-  }
-  glEnd();
+    // draw segments
+    glLineWidth(mLineWidth);
+    glColor3fv(mLineColor.GetArrayPtr());
+    glBegin(GL_LINE_STRIP);
+    // just draw the spline as a series of connected linear segments
+    for (size_t i = 0; i < mCoefficients.size(); i++) {
+        glVertex3fv(mCoefficients.at(i).GetArrayPtr());
+    }
+    glEnd();
 
-  // restore attribs
-  glPopAttrib();
+    // restore attribs
+    glPopAttrib();
 
-  glPopMatrix();
+    glPopMatrix();
 
-  GLObject::Render();
+    GLObject::Render();
 }