Răsfoiți Sursa

higher precision writeoff and writedmat, fixed bug in mat_max, more templates, sort rows, unique rows, double area

Former-commit-id: 7a5aff93f46c9e02d5ed98ab9e8c7b4bd1fde609
Alec Jacobson (jalec 12 ani în urmă
părinte
comite
b7fe54e95c

+ 5 - 5
Makefile

@@ -15,7 +15,7 @@ $(info Hello, $(IGL_USERNAME)!)
 all: LFLAGS +=
 OPTFLAGS=-O3 -DNDEBUG
 #debug: OPTFLAGS= -g -Wall -Werror
-debug: OPTFLAGS= -g -Wall -Weffc++
+debug: OPTFLAGS= -g -Wall
 CFLAGS += $(OPTFLAGS)
 
 EXTRA_DIRS=
@@ -39,10 +39,10 @@ endif
 .PHONY: examples
 .PHONY: extras
 debug: lib
-lib: obj lib/libigl.a
-examples: lib
+lib: lib/libigl.a
+examples: lib extras
 	make -C examples
-extras: extras
+extras:
 	for p in  $(EXTRA_DIRS); \
 	do \
 	echo "cd $$p" ; \
@@ -86,7 +86,7 @@ lib/libigl.a: $(OBJ_FILES)
 	rm -f $@
 	ar cqs $@ $(OBJ_FILES)
 
-obj/%.o: include/igl/%.cpp include/igl/%.h obj
+obj/%.o: include/igl/%.cpp include/igl/%.h
 	$(GG) $(CFLAGS) $(AFLAGS) -c -o $@ $< $(INC)
 
 clean:

+ 30 - 0
include/igl/IndexComparison.h

@@ -0,0 +1,30 @@
+#ifndef IGL_INDEXCOMPARISON_H
+#define IGL_INDEXCOMPARISON_H
+// Comparison struct used by sort
+// http://bytes.com/topic/c/answers/132045-sort-get-index
+namespace igl{
+
+  // For use with functions like std::sort
+  template<class T> struct IndexLessThan
+  {
+    IndexLessThan(const T arr) : arr(arr) {}
+    bool operator()(const size_t a, const size_t b) const
+    {
+      return arr[a] < arr[b];
+    }
+    const T arr;
+  };
+
+  // For use with functions like std::unique
+  template<class T> struct IndexEquals
+  {
+    IndexEquals(const T arr) : arr(arr) {}
+    bool operator()(const size_t a, const size_t b) const
+    {
+      return arr[a] == arr[b];
+    }
+    const T arr;
+  };
+}
+
+#endif

+ 60 - 0
include/igl/SortableRow.h

@@ -0,0 +1,60 @@
+#ifndef IGL_SORTABLE_ROW_H
+#define IGL_SORTABLE_ROW_H
+
+// Simple class to contain a rowvector which allows rowwise sorting and
+// reordering
+#include <Eigen/Core>
+
+// Templates:
+//   T  should be a matrix that implments .size(), and operator(int i)
+template <typename T>
+class SortableRow
+{
+  public:
+    T data;
+  public:
+    SortableRow():data(){};
+    SortableRow(const T & data):data(data){};
+    bool operator<(const SortableRow & that) const
+    {
+      // Get reference so that I can use parenthesis
+      const SortableRow<T> & THIS = *this;
+      // Lexicographical
+      int minc = (THIS.data.size() < that.data.size()? 
+          THIS.data.size() : that.data.size());
+      // loop over columns
+      for(int i = 0;i<minc;i++)
+      {
+        if(THIS.data(i) == that.data(i))
+        {
+          continue;
+        }
+        return THIS.data(i) < that.data(i);
+      }
+      // All characters the same, comes done to length
+      return THIS.data.size()<that.data.size();
+    };
+    bool operator==(const SortableRow & that) const
+    {
+      // Get reference so that I can use parenthesis
+      const SortableRow<T> & THIS = *this;
+      if(THIS.data.size() != that.data.size())
+      {
+        return false;
+      }
+      for(int i = 0;i<THIS.data.size();i++)
+      {
+        if(THIS.data(i) != that.data(i))
+        {
+          return false;
+        }
+      }
+      return true;
+    };
+    bool operator!=(const SortableRow & that) const
+    {
+      return !(*this == that);
+    };
+};
+
+#endif

+ 43 - 0
include/igl/doublearea.cpp

@@ -0,0 +1,43 @@
+#include "doublearea.h"
+#include "edge_lengths.h"
+#include <cassert>
+
+template <typename DerivedV, typename DerivedF, typename DeriveddblA>
+IGL_INLINE void igl::doublearea( 
+  const Eigen::PlainObjectBase<DerivedV> & V, 
+  const Eigen::PlainObjectBase<DerivedF> & F, 
+  Eigen::PlainObjectBase<DeriveddblA> & dblA)
+{
+  // Only support triangles
+  assert(F.cols() == 3);
+  // Compute edge lengths
+  Eigen::PlainObjectBase<DerivedV> l;
+  edge_lengths(V,F,l);
+  return doublearea(l,dblA);
+}
+
+template <typename Derivedl, typename DeriveddblA>
+IGL_INLINE void igl::doublearea( 
+  const Eigen::PlainObjectBase<Derivedl> & l, 
+  Eigen::PlainObjectBase<DeriveddblA> & dblA)
+{
+  // Only support triangles
+  assert(l.cols() == 3);
+  // Number of triangles
+  const int m = l.rows();
+  // semiperimeters
+  Eigen::PlainObjectBase<Derivedl> s = l.rowwise().sum()*0.5;
+  assert(s.rows() == m);
+  // resize output
+  dblA.resize(l.rows(),1);
+  // Heron's formula for area
+  for(int i = 0;i<m;i++)
+  {
+    dblA(i) = 2.0*sqrt(s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2)));
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+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> >&);
+#endif

+ 49 - 0
include/igl/doublearea.h

@@ -0,0 +1,49 @@
+#ifndef IGL_DOUBLEAREA_H
+#define IGL_DOUBLEAREA_H
+#include "igl_inline.h"
+#include <Eigen/Dense>
+namespace igl
+{
+  // DOUBLEAREA computes twice the area for each input triangle
+  //
+  // Templates:
+  //   DerivedV  derived type of eigen matrix for V (e.g. derived from
+  //     MatirxXd)
+  //   DerivedF  derived type of eigen matrix for F (e.g. derived from
+  //     MatrixXi)
+  //   DeriveddblA  derived type of eigen matrix for dblA (e.g. derived from
+  //     MatirxXd)
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by simplex_size list of mesh faces (must be triangles)
+  // Outputs:
+  //   dblA  #F list of triangle double areas
+  template <typename DerivedV, typename DerivedF, typename DeriveddblA>
+  IGL_INLINE void doublearea( 
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedF> & F, 
+    Eigen::PlainObjectBase<DeriveddblA> & dblA);
+  // Same as above but use instrinsic edge lengths rather than (V,F) mesh
+  // Templates:
+  //   DerivedV  derived type of eigen matrix for V (e.g. derived from
+  //     MatirxXd)
+  //   DerivedF  derived type of eigen matrix for F (e.g. derived from
+  //     MatrixXi)
+  //   DeriveddblA  derived type of eigen matrix for dblA (e.g. derived from
+  //     MatirxXd)
+  // Inputs:
+  //   l  #F by dim list of edge lengths using 
+  //     for triangles, columns correspond to edges 23,31,12
+  // Outputs:
+  //   dblA  #F list of triangle double areas
+  template <typename Derivedl, typename DeriveddblA>
+  IGL_INLINE void doublearea( 
+    const Eigen::PlainObjectBase<Derivedl> & l, 
+    Eigen::PlainObjectBase<DeriveddblA> & dblA);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "doublearea.h"
+#endif
+
+#endif

+ 1 - 0
include/igl/list_to_matrix.cpp

@@ -96,4 +96,5 @@ template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<d
 template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2> >&);
 template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >&);
 template bool igl::list_to_matrix<int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::vector<int, std::allocator<int> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
+template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

+ 4 - 3
include/igl/mat_max.cpp

@@ -18,9 +18,9 @@ IGL_INLINE void igl::mat_max(
   // loop over dimension opposite of dim
   for(int j = 0;j<n;j++)
   {
-    int PHONY;
-    int i;
-    int m;
+    typename Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::Index PHONY;
+    typename Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::Index i;
+    T m;
     if(dim==1)
     {
       m = X.col(j).maxCoeff(&i,&PHONY);
@@ -35,4 +35,5 @@ IGL_INLINE void igl::mat_max(
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
+template void igl::mat_max<double>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, int, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
 #endif

+ 0 - 2
include/igl/mat_min.cpp

@@ -1,7 +1,5 @@
 #include "mat_min.h"
 
-#include "verbose.h"
-
 template <typename T>
 IGL_INLINE void igl::mat_min(
   const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,

+ 8 - 0
include/igl/reorder.cpp

@@ -1,4 +1,8 @@
 #include "reorder.h"
+#include "SortableRow.h"
+#ifndef IGL_NO_EIGEN
+#include <Eigen/Core>
+#endif
 
 // This implementation is O(n), but also uses O(n) extra memory
 template< class T >
@@ -21,4 +25,8 @@ IGL_INLINE void igl::reorder(
 // Explicit template specialization
 // generated by autoexplicit.sh
 template void igl::reorder<double>(std::vector<double, std::allocator<double> > const&, std::vector<unsigned long, std::allocator<unsigned long> > const&, std::vector<double, std::allocator<double> >&);
+template void igl::reorder<int>(std::vector<int, std::allocator<int> > const&, std::vector<unsigned long, std::allocator<unsigned long> > const&, std::vector<int, std::allocator<int> >&);
+#ifndef IGL_NO_EIGEN
+template void igl::reorder<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > > const&, std::vector<unsigned long, std::allocator<unsigned long> > const&, std::vector<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > >&);
+#endif
 #endif

+ 9 - 16
include/igl/sort.cpp

@@ -1,7 +1,11 @@
 #include "sort.h"
 
-#include <algorithm>
+#include "SortableRow.h"
 #include "reorder.h"
+#include "IndexComparison.h"
+
+#include <cassert>
+#include <algorithm>
 
 template <typename DerivedX, typename DerivedIX>
 IGL_INLINE void igl::sort(
@@ -56,20 +60,6 @@ IGL_INLINE void igl::sort(
   }
 }
 
-// Comparison struct used by sort
-// http://bytes.com/topic/c/answers/132045-sort-get-index
-namespace igl{
-  template<class T> struct index_cmp
-  {
-    index_cmp(const T arr) : arr(arr) {}
-    bool operator()(const size_t a, const size_t b) const
-    {
-      return arr[a] < arr[b];
-    }
-    const T arr;
-  };
-}
-
 template <class T>
 IGL_INLINE void igl::sort(
   const std::vector<T> & unsorted,
@@ -87,7 +77,7 @@ IGL_INLINE void igl::sort(
   std::sort(
     index_map.begin(),
     index_map.end(),
-    igl::index_cmp<const std::vector<T>& >(unsorted));
+    igl::IndexLessThan<const std::vector<T>& >(unsorted));
 
   // if not ascending then reverse
   if(!ascending)
@@ -104,4 +94,7 @@ IGL_INLINE void igl::sort(
 // Explicit template specialization
 // generated by autoexplicit.sh
 template void igl::sort<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&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > > const&, bool, std::vector<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
+template void igl::sort<int>(std::vector<int, std::allocator<int> > const&, bool, std::vector<int, std::allocator<int> >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
 #endif

+ 1 - 0
include/igl/sort.h

@@ -46,6 +46,7 @@ namespace igl
       const bool ascending,
       std::vector<T> &sorted,
       std::vector<size_t> &index_map);
+
 }
 
 #ifdef IGL_HEADER_ONLY

+ 43 - 0
include/igl/sortrows.cpp

@@ -0,0 +1,43 @@
+#include "sortrows.h"
+
+#include "SortableRow.h"
+#include "sort.h"
+
+#include <vector>
+
+template <typename DerivedX, typename DerivedIX>
+IGL_INLINE void igl::sortrows(
+  const Eigen::PlainObjectBase<DerivedX>& X,
+  const bool ascending,
+  Eigen::PlainObjectBase<DerivedX>& Y,
+  Eigen::PlainObjectBase<DerivedIX>& IX)
+{
+  using namespace std;
+  using namespace Eigen;
+  typedef Eigen::Matrix<typename DerivedX::Scalar, Eigen::Dynamic, 1> RowVector;
+  vector<SortableRow<RowVector> > rows;
+  rows.resize(X.rows());
+  // Loop over rows
+  for(int i = 0;i<X.rows();i++)
+  {
+    RowVector ri = X.row(i);
+    rows[i] = SortableRow<RowVector>(ri);
+  }
+  vector<SortableRow<RowVector> > sorted;
+  std::vector<size_t> index_map;
+  // Perform sort on rows
+  igl::sort(rows,ascending,sorted,index_map);
+  // Resize output
+  Y.resize(X.rows(),X.cols());
+  IX.resize(X.rows(),1);
+  // Convert to eigen
+  for(int i = 0;i<X.rows();i++)
+  {
+    Y.row(i) = sorted[i].data;
+    IX(i) = index_map[i];
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+template void igl::sortrows<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

+ 35 - 0
include/igl/sortrows.h

@@ -0,0 +1,35 @@
+#ifndef IGL_SORTROWS_H
+#define IGL_SORTROWS_H
+#include "igl_inline.h"
+
+#include <vector>
+#include <Eigen/Core>
+namespace igl
+{
+  // Act like matlab's [Y,I] = sortrows(X)
+  //
+  // Templates:
+  //   DerivedX derived scalar type, e.g. MatrixXi or MatrixXd
+  //   DerivedIX derived integer type, e.g. MatrixXi
+  // Inputs:
+  //   X  m by n matrix whose entries are to be sorted
+  //   ascending  sort ascending (true, matlab default) or descending (false)
+  // Outputs:
+  //   Y  m by n matrix whose entries are sorted
+  //   IX  m by n matrix of indices so that if dim = 1, then in matlab notation
+  //     for j = 1:n, Y(:,j) = X(I(:,j),j); end
+  template <typename DerivedX, typename DerivedIX>
+  IGL_INLINE void sortrows(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    const bool ascending,
+    Eigen::PlainObjectBase<DerivedX>& Y,
+    Eigen::PlainObjectBase<DerivedIX>& IX);
+
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "sortrows.cpp"
+#endif
+
+#endif
+

+ 98 - 0
include/igl/unique.cpp

@@ -0,0 +1,98 @@
+#include "unique.h"
+#include "sort.h"
+#include "IndexComparison.h"
+#include "SortableRow.h"
+#include "list_to_matrix.h"
+
+#include <algorithm>
+#include <iostream>
+
+template <typename T>
+IGL_INLINE void igl::unique(
+  const std::vector<T> & A,
+  std::vector<T> & C,
+  std::vector<size_t> & IA,
+  std::vector<size_t> & IC)
+{
+  using namespace std;
+  std::vector<size_t> IM;
+  std::vector<T> sortA;
+  igl::sort(A,true,sortA,IM);
+  // Original unsorted index map
+  IA.resize(sortA.size());
+  for(int i=0;i<(int)sortA.size();i++)
+  {
+    IA[i] = i;
+  }
+  IA.erase(
+    std::unique(
+    IA.begin(),
+    IA.end(),
+    igl::IndexEquals<const std::vector<T>& >(sortA)),IA.end());
+
+  IC.resize(A.size());
+  {
+    int j = 0;
+    for(int i = 0;i<(int)sortA.size();i++)
+    {
+      if(sortA[IA[j]] != sortA[i])
+      {
+        j++;
+      }
+      IC[IM[i]] = j;
+    }
+  }
+
+  C.resize(IA.size());
+  // Reindex IA according to IM
+  for(int i = 0;i<(int)IA.size();i++)
+  {
+    IA[i] = IM[IA[i]];
+    C[i] = A[IA[i]];
+  }
+}
+
+template <typename DerivedA, typename DerivedIA, typename DerivedIC>
+IGL_INLINE void igl::unique_rows(
+  const Eigen::PlainObjectBase<DerivedA>& A,
+  Eigen::PlainObjectBase<DerivedA>& C,
+  Eigen::PlainObjectBase<DerivedIA>& IA,
+  Eigen::PlainObjectBase<DerivedIC>& IC)
+{
+  using namespace std;
+
+  typedef Eigen::Matrix<typename DerivedA::Scalar, Eigen::Dynamic, 1> RowVector;
+  vector<SortableRow<RowVector> > rows;
+  rows.resize(A.rows());
+  // Loop over rows
+  for(int i = 0;i<A.rows();i++)
+  {
+    RowVector ri = A.row(i);
+    rows[i] = SortableRow<RowVector>(ri);
+  }
+  vector<SortableRow<RowVector> > vC;
+
+  // unique on rows
+  vector<size_t> vIA;
+  vector<size_t> vIC;
+  unique(rows,vC,vIA,vIC);
+
+  // Convert to eigen
+  C.resize(vC.size(),A.cols());
+  IA.resize(vIA.size(),1);
+  IC.resize(vIC.size(),1);
+  for(int i = 0;i<C.rows();i++)
+  {
+    C.row(i) = vC[i].data;
+    IA(i) = vIA[i];
+  }
+  for(int i = 0;i<A.rows();i++)
+  {
+    IC(i) = vIC[i];
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+template void igl::unique<int>(std::vector<int, std::allocator<int> > const&, std::vector<int, std::allocator<int> >&, std::vector<unsigned long, std::allocator<unsigned long> >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
+template void igl::unique_rows<Eigen::Matrix<int, -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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -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> >&);
+#endif

+ 53 - 0
include/igl/unique.h

@@ -0,0 +1,53 @@
+#ifndef IGL_UNIQUE_H
+#define IGL_UNIQUE_H
+#include "igl_inline.h"
+
+#include <vector>
+#include <Eigen/Core>
+namespace igl
+{
+  // Act like matlab's [C,IA,IC] = unique(X)
+  //
+  // Templates:
+  //   T  comparable type T
+  // Inputs:
+  //   A  #A vector of type T
+  // Outputs:
+  //   C  #C vector of unique entries in A
+  //   IA  #A index vector so that C = A(IA);
+  //   IC  #C index vector so that A = C(IC);
+  template <typename T>
+  IGL_INLINE void unique(
+    const std::vector<T> & A,
+    std::vector<T> & C,
+    std::vector<size_t> & IA,
+    std::vector<size_t> & IC);
+
+  // Act like matlab's [C,IA,IC] = unique(X,'rows')
+  //
+  // Templates:
+  //   DerivedA derived scalar type, e.g. MatrixXi or MatrixXd
+  //   DerivedIA derived integer type, e.g. MatrixXi
+  //   DerivedIC derived integer type, e.g. MatrixXi
+  // Inputs:
+  //   A  m by n matrix whose entries are to unique'd according to rows
+  // Outputs:
+  //   C  #C vector of unique rows in A
+  //   IA  #A index vector so that C = A(IA,:);
+  //   IC  #C index vector so that A = C(IC,:);
+  template <typename DerivedA, typename DerivedIA, typename DerivedIC>
+  IGL_INLINE void unique_rows(
+    const Eigen::PlainObjectBase<DerivedA>& A,
+    Eigen::PlainObjectBase<DerivedA>& C,
+    Eigen::PlainObjectBase<DerivedIA>& IA,
+    Eigen::PlainObjectBase<DerivedIC>& IC);
+
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "unique.cpp"
+#endif
+
+#endif
+
+

+ 3 - 1
include/igl/writeDMAT.cpp

@@ -22,7 +22,7 @@ IGL_INLINE bool igl::writeDMAT(const std::string file_name, const Mat & W)
     // loop over rows (down columns) quickly
     for(int i = 0;i < W.rows();i++)
     {
-      fprintf(fp,"%lg\n",(double)W(i,j));
+      fprintf(fp,"%0.17lg\n",(double)W(i,j));
     }
   }
   fclose(fp);
@@ -68,4 +68,6 @@ IGL_INLINE bool igl::writeDMAT(
 // generated by autoexplicit.sh
 template bool igl::writeDMAT<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&);
 template bool igl::writeDMAT<double>(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >);
+template bool igl::writeDMAT<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&);
+template bool igl::writeDMAT<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&);
 #endif

+ 28 - 21
include/igl/writeOFF.cpp

@@ -4,29 +4,36 @@
 // write mesh to an ascii off file
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE bool igl::writeOFF(
-                              const std::string fname,
-                              const Eigen::PlainObjectBase<DerivedV>& V,
-                              const Eigen::PlainObjectBase<DerivedF>& F)
+  const std::string fname,
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F)
 {
-    FILE *fp = fopen (fname.c_str(), "w");
+  FILE *fp = fopen (fname.c_str(), "w");
   
-    
-    if (!fp)
-    {
-        fprintf (stderr, "writeOFF(): could not open file %s", fname.c_str());
-      return false;
-    }
-    
-    fprintf (fp, "OFF\n%d %d 0\n",  (int) V.rows(), (int) F.rows());
-    
-    for (int i = 0; i < V.rows(); i++)
-        fprintf (fp, "%f %f %f\n", V(i,0), V(i,1), V(i,2));
-    
-    for (int i = 0; i < F.rows(); i++)
-        fprintf (fp, "3 %d %d %d\n", F(i,0), F(i,1), F(i,2));
-    
-    fclose (fp);
-    return true;
+  
+  if (!fp)
+  {
+      fprintf (stderr, "writeOFF(): could not open file %s", fname.c_str());
+    return false;
+  }
+  
+  fprintf (fp, "OFF\n%d %d 0\n",  (int) V.rows(), (int) F.rows());
+  
+  for (int i = 0; i < V.rows(); i++)
+  {
+    fprintf(
+      fp,
+      "%0.17g %0.17g %0.17g\n",
+      (double)V(i,0), 
+      (double)V(i,1), 
+      (double)V(i,2));
+  }
+  
+  for (int i = 0; i < F.rows(); i++)
+      fprintf (fp, "3 %d %d %d\n", F(i,0), F(i,1), F(i,2));
+  
+  fclose (fp);
+  return true;
 }
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization