diff --git a/examples/intersections/example.cpp b/examples/intersections/example.cpp
index aa7dc4c7c..57d40e98e 100644
--- a/examples/intersections/example.cpp
+++ b/examples/intersections/example.cpp
@@ -23,6 +23,7 @@
 #include <igl/ReAntTweakBar.h>
 #include <igl/pathinfo.h>
 #include <igl/Camera.h>
+#include <igl/cat.h>
 #include <igl/get_seconds.h>
 #include <igl/cgal/remesh_self_intersections.h>
 #include <igl/cgal/intersect_other.h>
@@ -123,7 +124,7 @@ float light_pos[4] = {0.1,0.1,-0.9,0};
 // C,D  Colors
 // N,W  Normals
 // mid  combined "centroid"
-Eigen::MatrixXd V,N,C,Z,mid,U,W,D;
+Eigen::MatrixXd V,N,C,Z,mid,U,W,D,VU;
 // F,G  faces
 Eigen::MatrixXi F,G;
 bool has_other = false;
@@ -305,7 +306,7 @@ void display()
   // Draw a nice floor
   glPushMatrix();
   const double floor_offset =
-    -2./bbd*(V.col(1).maxCoeff()-mid(1));
+    -2./bbd*(VU.col(1).maxCoeff()-mid(1));
   glTranslated(0,floor_offset,0);
   const float GREY[4] = {0.5,0.5,0.6,1.0};
   const float DARK_GREY[4] = {0.2,0.2,0.3,1.0};
@@ -633,18 +634,14 @@ int main(int argc, char * argv[])
     {
       return 1;
     }
-    mid = 0.25*(V.colwise().maxCoeff() + V.colwise().minCoeff()) +
-      0.25*(U.colwise().maxCoeff() + U.colwise().minCoeff());
-    bbd = max(
-      (V.colwise().maxCoeff() - V.colwise().minCoeff()).maxCoeff(),
-      (U.colwise().maxCoeff() - U.colwise().minCoeff()).maxCoeff());
+    cat(1,V,U,VU);
     color_intersections(V,F,U,G,C,D);
   }else
   {
-    mid = 0.5*(V.colwise().maxCoeff() + V.colwise().minCoeff());
-    bbd = (V.colwise().maxCoeff() - V.colwise().minCoeff()).maxCoeff();
-    color_selfintersections(V,F,C);
+    VU = V;
   }
+  mid = 0.5*(VU.colwise().maxCoeff() + VU.colwise().minCoeff());
+  bbd = (VU.colwise().maxCoeff() - VU.colwise().minCoeff()).maxCoeff();
 
   // Init glut
   glutInit(&argc,argv);
diff --git a/include/igl/WindingNumberAABB.h b/include/igl/WindingNumberAABB.h
index dac211d61..9b78fdfb7 100644
--- a/include/igl/WindingNumberAABB.h
+++ b/include/igl/WindingNumberAABB.h
@@ -1,3 +1,15 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+// # MUTUAL DEPENDENCY ISSUE FOR HEADER ONLY VERSION
+// MUST INCLUDE winding_number.h first before guard:
+#include "winding_number.h"
+
 #ifndef IGL_WINDINGNUMBERAABB_H
 #define IGL_WINDINGNUMBERAABB_H
 #include "WindingNumberTree.h"
diff --git a/include/igl/WindingNumberTree.h b/include/igl/WindingNumberTree.h
index 8bd19c814..33ece46b7 100644
--- a/include/igl/WindingNumberTree.h
+++ b/include/igl/WindingNumberTree.h
@@ -1,5 +1,12 @@
-#ifndef IGL_BOUNDINGTREE_H
-#define IGL_BOUNDINGTREE_H
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_WINDINGNUMBERTREE_H
+#define IGL_WINDINGNUMBERTREE_H
 #include <list>
 #include <map>
 #include <Eigen/Dense>
diff --git a/include/igl/cgal/point_mesh_squared_distance.cpp b/include/igl/cgal/point_mesh_squared_distance.cpp
index 6fa9c7c89..2bee4138b 100644
--- a/include/igl/cgal/point_mesh_squared_distance.cpp
+++ b/include/igl/cgal/point_mesh_squared_distance.cpp
@@ -18,26 +18,72 @@ IGL_INLINE void igl::point_mesh_squared_distance(
   Eigen::MatrixXd & C)
 {
   using namespace std;
-  typedef CGAL::Point_3<Kernel>    Point_3;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  Tree tree;
+  vector<Triangle_3> T;
+  point_mesh_squared_distance_precompute(V,F,tree,T);
+  return point_mesh_squared_distance(P,tree,T,sqrD,I,C);
+}
+
+template <typename Kernel>
+IGL_INLINE void igl::point_mesh_squared_distance_precompute(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  std::vector<CGAL::Triangle_3<Kernel> > & T)
+{
+  using namespace std;
+
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+
+  // Must be 3D
+  assert(V.cols() == 3);
+  // Must be triangles
+  assert(F.cols() == 3);
+  // Make list of cgal triangles
+  mesh_to_cgal_triangle_list(V,F,T);
+  tree.clear();
+  tree.insert(T.begin(),T.end());
+  tree.accelerate_distance_queries();
+}
+
+template <typename Kernel>
+IGL_INLINE void igl::point_mesh_squared_distance(
+  const Eigen::MatrixXd & P,
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const std::vector<CGAL::Triangle_3<Kernel> > & T,
+  Eigen::VectorXd & sqrD,
+  Eigen::VectorXi & I,
+  Eigen::MatrixXd & C)
+{
   typedef CGAL::Triangle_3<Kernel> Triangle_3; 
   typedef typename std::vector<Triangle_3>::iterator Iterator;
   typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
   typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
   typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
   typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
-
-  // Must be 3D
-  assert(V.cols() == 3);
+  typedef CGAL::Point_3<Kernel>    Point_3;
   assert(P.cols() == 3);
-  // Must be triangles
-  assert(F.cols() == 3);
-  // Make list of cgal triangles
-  Tree tree;
-  vector<Triangle_3> T;
-  mesh_to_cgal_triangle_list(V,F,T);
-  tree.insert(T.begin(),T.end());
-
-  tree.accelerate_distance_queries();
   const int n = P.rows();
   sqrD.resize(n,1);
   I.resize(n,1);
@@ -59,5 +105,6 @@ IGL_INLINE void igl::point_mesh_squared_distance(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
+template void igl::point_mesh_squared_distance_precompute<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
 
 #endif
diff --git a/include/igl/cgal/point_mesh_squared_distance.h b/include/igl/cgal/point_mesh_squared_distance.h
index 964baef84..2d8aafc44 100644
--- a/include/igl/cgal/point_mesh_squared_distance.h
+++ b/include/igl/cgal/point_mesh_squared_distance.h
@@ -9,6 +9,7 @@
 #define IGL_POINT_MESH_SQUARED_DISTANCE_H
 #include <igl/igl_inline.h>
 #include <Eigen/Core>
+#include <vector>
 #include "CGAL_includes.hpp"
 namespace igl
 {
@@ -36,7 +37,45 @@ namespace igl
     Eigen::VectorXd & sqrD,
     Eigen::VectorXi & I,
     Eigen::MatrixXd & C);
+  // Probably can do this in a way that we don't pass around `tree` and `T`
+  //
+  // Outputs:
+  //   tree  CGAL's AABB tree
+  //   T  list of CGAL triangles in order of F (for determining which was found
+  //     in computation)
+  template <typename Kernel>
+  IGL_INLINE void point_mesh_squared_distance_precompute(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    std::vector<CGAL::Triangle_3<Kernel> > & T);
+  // Inputs:
+  //  see above
+  // Outputs:
+  //  see above
+  template <typename Kernel>
+  IGL_INLINE void point_mesh_squared_distance(
+    const Eigen::MatrixXd & P,
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const std::vector<CGAL::Triangle_3<Kernel> > & T,
+    Eigen::VectorXd & sqrD,
+    Eigen::VectorXi & I,
+    Eigen::MatrixXd & C);
+
 }
+
 #ifndef IGL_STATIC_LIBRARY
 #  include "point_mesh_squared_distance.cpp"
 #endif
diff --git a/include/igl/doublearea.cpp b/include/igl/doublearea.cpp
index 71d632561..b580db54d 100644
--- a/include/igl/doublearea.cpp
+++ b/include/igl/doublearea.cpp
@@ -139,6 +139,7 @@ IGL_INLINE void igl::doublearea(
 template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::doublearea<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
diff --git a/include/igl/edge_lengths.cpp b/include/igl/edge_lengths.cpp
index f58d6039f..772d7748c 100644
--- a/include/igl/edge_lengths.cpp
+++ b/include/igl/edge_lengths.cpp
@@ -72,4 +72,5 @@ template void igl::edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen:
 template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
 template void igl::edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif
diff --git a/include/igl/exterior_edges.h b/include/igl/exterior_edges.h
index 9a1f8135e..b003f86ee 100644
--- a/include/igl/exterior_edges.h
+++ b/include/igl/exterior_edges.h
@@ -20,7 +20,7 @@ namespace igl
   Eigen::MatrixXi exterior_edges( const Eigen::MatrixXi & F);
 }
 #ifndef IGL_STATIC_LIBRARY
-#  include "exterior_edges.h"
+#  include "exterior_edges.cpp"
 #endif
 
 #endif
diff --git a/include/igl/internal_angles.cpp b/include/igl/internal_angles.cpp
index 9e9c5f07b..62808182f 100644
--- a/include/igl/internal_angles.cpp
+++ b/include/igl/internal_angles.cpp
@@ -45,4 +45,5 @@ IGL_INLINE void igl::internal_angles(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::internal_angles<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::internal_angles<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif
diff --git a/include/igl/per_edge_normals.cpp b/include/igl/per_edge_normals.cpp
index 702cde41a..ccac7cf43 100644
--- a/include/igl/per_edge_normals.cpp
+++ b/include/igl/per_edge_normals.cpp
@@ -14,6 +14,7 @@ template <
 IGL_INLINE void igl::per_edge_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
+  const PerEdgeNormalsWeightingType weighting,
   Eigen::PlainObjectBase<DerivedN> & N,
   Eigen::PlainObjectBase<DerivedE> & E,
   Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
@@ -34,20 +35,51 @@ IGL_INLINE void igl::per_edge_normals(
   MatrixXd FN;
   per_face_normals(V,F,FN);
 
-  VectorXd dblA;
-  doublearea(V,F,dblA);
+  Eigen::VectorXd W(F.rows());
+  switch(weighting)
+  {
+    case PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM:
+      W.setConstant(1.);
+      break;
+    default:
+      assert(false && "Unknown weighting type");
+    case PER_EDGE_NORMALS_WEIGHTING_TYPE_DEFAULT:
+    case PER_EDGE_NORMALS_WEIGHTING_TYPE_AREA:
+    {
+      doublearea(V,F,W);
+      break;
+    }
+  }
+
   N.setConstant(E.rows(),3,0);
   for(int f = 0;f<m;f++)
   {
     for(int c = 0;c<3;c++)
     {
-      N.row(EMAP(f+c*m)) += dblA(f) * FN.row(f);
+      N.row(EMAP(f+c*m)) += W(f) * FN.row(f);
     }
   }
   N.rowwise().normalize();
   
 }
 
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedN,
+  typename DerivedE,
+  typename DerivedEMAP>
+IGL_INLINE void igl::per_edge_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedN> & N,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  return 
+    per_edge_normals(V,F,PER_EDGE_NORMALS_WEIGHTING_TYPE_DEFAULT,N,E,EMAP);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciation
 template void igl::per_edge_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
diff --git a/include/igl/per_edge_normals.h b/include/igl/per_edge_normals.h
index 4a5a6eb2f..4b12aa9c2 100644
--- a/include/igl/per_edge_normals.h
+++ b/include/igl/per_edge_normals.h
@@ -11,15 +11,39 @@
 #include <Eigen/Core>
 namespace igl
 {
+  enum PerEdgeNormalsWeightingType
+  {
+    // Incident face normals have uniform influence on edge normal
+    PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM = 0,
+    // Incident face normals are averaged weighted by area
+    PER_EDGE_NORMALS_WEIGHTING_TYPE_AREA = 1,
+    // Area weights
+    PER_EDGE_NORMALS_WEIGHTING_TYPE_DEFAULT = 2,
+    NUM_PER_EDGE_NORMALS_WEIGHTING_TYPE = 3
+  };
   // Compute face normals via vertex position list, face list
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F  #F by 3 eigen Matrix of face (triangle) indices
+  //   weighting  weighting type
   // Output:
   //   N  #2 by 3 matrix of mesh edge 3D normals per row
   //   E  #E by 2 matrix of edge indices per row
   //   EMAP  #E by 1 matrix of indices from all edges to E
   //
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedN,
+    typename DerivedE,
+    typename DerivedEMAP>
+  IGL_INLINE void per_edge_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const PerEdgeNormalsWeightingType weight,
+    Eigen::PlainObjectBase<DerivedN> & N,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
   template <
     typename DerivedV, 
     typename DerivedF, 
diff --git a/include/igl/per_face_normals.cpp b/include/igl/per_face_normals.cpp
index 013f18ada..91b78fd89 100644
--- a/include/igl/per_face_normals.cpp
+++ b/include/igl/per_face_normals.cpp
@@ -19,7 +19,7 @@ IGL_INLINE void igl::per_face_normals(
   N.resize(F.rows(),3);
   // loop over faces
   int Frows = F.rows();
-#pragma omp parallel for
+#pragma omp parallel for if (Frows>10000)
   for(int i = 0; i < Frows;i++)
   {
     const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v1 = V.row(F(i,1)) - V.row(F(i,0));
diff --git a/include/igl/per_vertex_normals.cpp b/include/igl/per_vertex_normals.cpp
index 9e3a08890..0b2ae8e08 100644
--- a/include/igl/per_vertex_normals.cpp
+++ b/include/igl/per_vertex_normals.cpp
@@ -8,28 +8,62 @@
 #include "per_vertex_normals.h"
 
 #include "per_face_normals.h"
+#include "doublearea.h"
+#include "internal_angles.h"
 
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::per_vertex_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
+  const igl::PerVertexNormalsWeightingType weighting,
   Eigen::PlainObjectBase<DerivedV> & N)
 {
   Eigen::PlainObjectBase<DerivedV> PFN;
   igl::per_face_normals(V,F,PFN);
-  return igl::per_vertex_normals(V,F,PFN,N);
+  return per_vertex_normals(V,F,weighting,PFN,N);
 }
 
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::per_vertex_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedV> & N)
+{
+  return per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT,N);
+}
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::per_vertex_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const igl::PerVertexNormalsWeightingType weighting,
   const Eigen::PlainObjectBase<DerivedV>& FN,
   Eigen::PlainObjectBase<DerivedV> & N)
 {
   // Resize for output
   N.setZero(V.rows(),3);
 
+  Eigen::MatrixXd W(F.rows(),3);
+  switch(weighting)
+  {
+    case PER_VERTEX_NORMALS_WEIGHTING_TYPE_UNIFORM:
+      W.setConstant(1.);
+      break;
+    default:
+      assert(false && "Unknown weighting type");
+    case PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT:
+    case PER_VERTEX_NORMALS_WEIGHTING_TYPE_AREA:
+    {
+      Eigen::VectorXd A;
+      doublearea(V,F,A);
+      W = A.replicate(1,3);
+      break;
+    }
+    case PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE:
+      internal_angles(V,F,W);
+      break;
+  }
+
   // loop over faces
   const int Frows = F.rows();
 //// Minimum number of iterms per openmp thread
@@ -44,15 +78,29 @@ IGL_INLINE void igl::per_vertex_normals(
     {
       // Does this need to be critical?
 //#pragma omp critical
-      N.row(F(i,j)) += FN.row(i);
+      N.row(F(i,j)) += W(i,j)*FN.row(i);
     }
   }
+  // take average via normalization
   N.rowwise().normalize();
 }
 
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::per_vertex_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedV>& FN,
+  Eigen::PlainObjectBase<DerivedV> & N)
+{
+  return
+    per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT,FN,N);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::per_vertex_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::PerVertexNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::per_vertex_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::per_vertex_normals<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
 template void igl::per_vertex_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
diff --git a/include/igl/per_vertex_normals.h b/include/igl/per_vertex_normals.h
index b8a844fe0..3cc1b0deb 100644
--- a/include/igl/per_vertex_normals.h
+++ b/include/igl/per_vertex_normals.h
@@ -9,19 +9,37 @@
 #define IGL_PER_VERTEX_NORMALS_H
 #include "igl_inline.h"
 #include <Eigen/Core>
-// Note: So for this only computes normals per vertex as uniformly weighted
-// averages of incident triangle normals. It would be nice to support more or
-// all of the methods here:
+// Note: It would be nice to support more or all of the methods here:
 // "A comparison of algorithms for vertex normal computation"
 namespace igl
 {
+  enum PerVertexNormalsWeightingType
+  {
+    // Incident face normals have uniform influence on vertex normal
+    PER_VERTEX_NORMALS_WEIGHTING_TYPE_UNIFORM = 0,
+    // Incident face normals are averaged weighted by area
+    PER_VERTEX_NORMALS_WEIGHTING_TYPE_AREA = 1,
+    // Incident face normals are averaged weighted by incident angle of vertex
+    PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE = 2,
+    // Area weights
+    PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT = 3,
+    NUM_PER_VERTEX_NORMALS_WEIGHTING_TYPE = 4
+  };
   // Compute vertex normals via vertex position list, face list
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   weighting  Weighting type
   // Output:
   //   N  #V by 3 eigen Matrix of mesh vertex 3D normals
   template <typename DerivedV, typename DerivedF>
+  IGL_INLINE void per_vertex_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const igl::PerVertexNormalsWeightingType weighting,
+    Eigen::PlainObjectBase<DerivedV> & N);
+  // Without weighting
+  template <typename DerivedV, typename DerivedF>
   IGL_INLINE void per_vertex_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
@@ -29,6 +47,14 @@ namespace igl
   // Inputs:
   //   FN  #F by 3 matrix of face (triangle) normals
   template <typename DerivedV, typename DerivedF>
+  IGL_INLINE void per_vertex_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const PerVertexNormalsWeightingType weighting,
+    const Eigen::PlainObjectBase<DerivedV>& FN,
+    Eigen::PlainObjectBase<DerivedV> & N);
+  // Without weighting
+  template <typename DerivedV, typename DerivedF>
   IGL_INLINE void per_vertex_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
diff --git a/include/igl/triangle_fan.h b/include/igl/triangle_fan.h
index d6dcb45e7..2ebc9b3da 100644
--- a/include/igl/triangle_fan.h
+++ b/include/igl/triangle_fan.h
@@ -18,6 +18,6 @@ namespace igl
   IGL_INLINE Eigen::MatrixXi triangle_fan( const Eigen::MatrixXi & E);
 }
 #ifndef IGL_STATIC_LIBRARY
-#  include "triangle_fan.h"
+#  include "triangle_fan.cpp"
 #endif
 #endif
diff --git a/include/igl/winding_number.h b/include/igl/winding_number.h
index 409bad10e..f3eb9e32f 100644
--- a/include/igl/winding_number.h
+++ b/include/igl/winding_number.h
@@ -61,7 +61,7 @@ namespace igl
 }
 
 #ifndef IGL_STATIC_LIBRARY
-#  include "winding_number.h"
+#  include "winding_number.cpp"
 #endif
 
 #endif