split most .h files into .h and .cpp, most templated .cpp are missing any explicit instanciations so libigl.a builds but is rather useless

Former-commit-id: b328d179fbe47be83f8876aeab384788623f5613

EPS.cpp

@@ -0,0 +1,15 @@
+#include "EPS.h"
+
+template <> IGL_INLINE float igl::EPS()
+{
+  return igl::FLOAT_EPS;
+}
+
+template <> IGL_INLINE double igl::EPS()
+{
+  return igl::DOUBLE_EPS;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

EPS.h

@@ -1,5 +1,6 @@
 #ifndef IGL_EPS_H
 #define IGL_EPS_H
+#include "igl_inline.h"
 // Define a standard value for double epsilon
 namespace igl
 {
@@ -8,20 +9,14 @@ namespace igl
   const float FLOAT_EPS    = 1.0e-7;
   const float FLOAT_EPS_SQ = 1.0e-14;
   // Function returning EPS for corresponding type
-  template <typename S_type> inline S_type EPS();
+  template <typename S_type> IGL_INLINE S_type EPS();
   // Template specializations for float and double
-  template <> inline float EPS<float>();
-  template <> inline double EPS<double>();
+  template <> IGL_INLINE float EPS<float>();
+  template <> IGL_INLINE double EPS<double>();
 }
 
-// Implementation
-template <> inline float igl::EPS()
-{
-  return igl::FLOAT_EPS;
-}
+#ifdef IGL_HEADER_ONLY
+#  include "EPS.cpp"
+#endif
 
-template <> inline double igl::EPS()
-{
-  return igl::DOUBLE_EPS;
-}
 #endif

ReAntTweakBar.cpp

@@ -0,0 +1,679 @@
+#include "ReAntTweakBar.h"
+
+#include <cstdio>
+#include <sstream>
+#include <map>
+
+#define MAX_CB_VAR_SIZE 10
+// Max line size for reading files
+#define MAX_LINE 1000
+#define MAX_WORD 100
+
+// GLOBAL WRAPPERS
+namespace igl
+{
+  std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > > ReTw_custom_types;
+}
+
+TwType igl::ReTwDefineEnum(
+  const char *name, 
+  const TwEnumVal *enumValues, 
+  unsigned int nbValues)
+{
+  // copy enum valus into vector
+  std::vector<TwEnumVal> enum_vals;
+  enum_vals.resize(nbValues);
+  for(unsigned int j = 0; j<nbValues;j++)
+  {
+    enum_vals[j] = enumValues[j];
+  }
+  TwType type = TwDefineEnum(name,enumValues,nbValues);
+
+  ReTw_custom_types[type] = 
+    std::pair<const char *,std::vector<TwEnumVal> >(name,enum_vals);
+  return type;
+}
+
+namespace igl
+{
+  struct ReTwTypeString
+  {
+    TwType type;
+    const char * type_str;
+  };
+
+  #define RETW_NUM_DEFAULT_TYPE_STRINGS 23
+  ReTwTypeString ReTwDefaultTypeStrings[RETW_NUM_DEFAULT_TYPE_STRINGS] = 
+  {
+    {TW_TYPE_UNDEF,"TW_TYPE_UNDEF"},
+    {TW_TYPE_BOOLCPP,"TW_TYPE_BOOLCPP"},
+    {TW_TYPE_BOOL8,"TW_TYPE_BOOL8"},
+    {TW_TYPE_BOOL16,"TW_TYPE_BOOL16"},
+    {TW_TYPE_BOOL32,"TW_TYPE_BOOL32"},
+    {TW_TYPE_CHAR,"TW_TYPE_CHAR"},
+    {TW_TYPE_INT8,"TW_TYPE_INT8"},
+    {TW_TYPE_UINT8,"TW_TYPE_UINT8"},
+    {TW_TYPE_INT16,"TW_TYPE_INT16"},
+    {TW_TYPE_UINT16,"TW_TYPE_UINT16"},
+    {TW_TYPE_INT32,"TW_TYPE_INT32"},
+    {TW_TYPE_UINT32,"TW_TYPE_UINT32"},
+    {TW_TYPE_FLOAT,"TW_TYPE_FLOAT"},
+    {TW_TYPE_DOUBLE,"TW_TYPE_DOUBLE"},
+    {TW_TYPE_COLOR32,"TW_TYPE_COLOR32"},
+    {TW_TYPE_COLOR3F,"TW_TYPE_COLOR3F"},
+    {TW_TYPE_COLOR4F,"TW_TYPE_COLOR4F"},
+    {TW_TYPE_CDSTRING,"TW_TYPE_CDSTRING"},
+    {TW_TYPE_STDSTRING,"TW_TYPE_STDSTRING"},
+    {TW_TYPE_QUAT4F,"TW_TYPE_QUAT4F"},
+    {TW_TYPE_QUAT4D,"TW_TYPE_QUAT4D"},
+    {TW_TYPE_DIR3F,"TW_TYPE_DIR3F"},
+    {TW_TYPE_DIR3D,"TW_TYPE_DIR3D"}
+  };
+}
+
+
+// BAR WRAPPERS
+void igl::ReTwBar::TwNewBar(const char *barName)
+{
+  // double colon without anything in front of it means look for this in the
+  // global namespace... I hope...
+  this->bar = ::TwNewBar(barName);
+}
+
+int igl::ReTwBar::TwAddVarRW(
+  const char *name, 
+  TwType type, 
+  void *var, 
+  const char *def,
+  const bool record)
+{
+  int ret = ::TwAddVarRW(this->bar,name,type,var,def);
+  if(ret && record)
+  {
+    rw_items.push_back(ReTwRWItem(name,type,var));
+  }
+  return ret;
+}
+
+int igl::ReTwBar::TwAddVarCB(
+  const char *name, 
+  TwType type, 
+  TwSetVarCallback setCallback, 
+  TwGetVarCallback getCallback, 
+  void *clientData, 
+  const char *def,
+  const bool record)
+{
+  int ret = 
+    ::TwAddVarCB(this->bar,name,type,setCallback,getCallback,clientData,def);
+  if(ret && record)
+  {
+    cb_items.push_back(ReTwCBItem(name,type,setCallback,getCallback,clientData));
+  }
+  return ret;
+}
+
+int igl::ReTwBar::TwAddVarRO(
+  const char *name, 
+  TwType type, 
+  void *var, 
+  const char *def)
+{
+  int ret = ::TwAddVarRO(this->bar,name,type,var,def);
+  // Read only variables are not recorded
+  //if(ret)
+  //{
+  //  rw_items.push_back(ReTwRWItem(name,type,var));
+  //}
+  return ret;
+}
+
+int igl::ReTwBar::TwAddButton(
+  const char *name, 
+  TwButtonCallback buttonCallback, 
+  void *clientData, 
+  const char *def)
+{
+  int ret = 
+    ::TwAddButton(this->bar,name,buttonCallback,clientData,def);
+  // buttons are not recorded
+  //if(ret)
+  //{
+  //  cb_items.push_back(ReTwCBItem(name,type,setCallback,getCallback,clientData));
+  //}
+  return ret;
+}
+
+int igl::ReTwBar::TwSetParam(
+  const char *varName, 
+  const char *paramName, 
+  TwParamValueType paramValueType, 
+  unsigned int inValueCount, 
+  const void *inValues)
+{
+  // For now just pass these along
+  return 
+    ::TwSetParam(
+      this->bar,
+      varName,
+      paramName,
+      paramValueType,
+      inValueCount,
+      inValues);
+}
+
+int igl::ReTwBar::TwGetParam(
+  const char *varName, 
+  const char *paramName, 
+  TwParamValueType paramValueType, 
+  unsigned int outValueMaxCount,
+  void *outValues)
+{
+  return 
+    ::TwGetParam(
+      this->bar,
+      varName,
+      paramName,
+      paramValueType,
+      outValueMaxCount,
+      outValues);
+}
+
+int igl::ReTwBar::TwRefreshBar()
+{
+  return ::TwRefreshBar(this->bar);
+}
+bool igl::ReTwBar::save(const char *file_name)
+{
+  FILE * fp;
+  if(file_name == NULL)
+  {
+    fp = stdout;
+  }else
+  {
+    fp = fopen(file_name,"w");
+  }
+
+  if(fp == NULL)
+  {
+    printf("ERROR: not able to open %s for writing...\n",file_name);
+    return false;
+  }
+
+  // Print all RW variables
+  for(
+    std::vector<ReTwRWItem>::iterator it = rw_items.begin(); 
+    it != rw_items.end(); 
+    it++)
+  {
+    const char * name = (*it).name;
+    TwType type = (*it).type;
+    void * var = (*it).var;
+    fprintf(fp,"%s: %s\n",
+      name,
+      get_value_as_string(var,type).c_str());
+  }
+
+  char var[MAX_CB_VAR_SIZE];
+  // Print all CB variables
+  for(
+    std::vector<ReTwCBItem>::iterator it = cb_items.begin(); 
+    it != cb_items.end(); 
+    it++)
+  {
+    const char * name = it->name;
+    TwType type = it->type;
+    //TwSetVarCallback setCallback = it->setCallback;
+    TwGetVarCallback getCallback = it->getCallback;
+    void * clientData = it->clientData;
+    // I'm not sure how to do what I want to do. getCallback needs to be sure
+    // that it can write to var. So var needs to point to a valid and big
+    // enough chunk of memory
+    getCallback(var,clientData);
+    fprintf(fp,"%s: %s\n",
+      name,
+      get_value_as_string(var,type).c_str());
+  }
+
+  fprintf(fp,"\n");
+
+  if(file_name != NULL)
+  {
+    fclose(fp);
+  }
+  // everything succeeded
+  return true;
+}
+
+std::string igl::ReTwBar::get_value_as_string(
+  void * var, 
+  TwType type)
+{
+  std::stringstream sstr;
+  switch(type)
+  {
+    case TW_TYPE_BOOLCPP:
+      {
+        sstr << "TW_TYPE_BOOLCPP" << " ";
+        sstr << *(static_cast<bool*>(var));
+        break;
+      }
+    case TW_TYPE_QUAT4F:
+      {
+        sstr << "TW_TYPE_QUAT4F" << " ";
+        // Q: Why does casting to float* work? shouldn't I have to cast to
+        // float**?
+        float * q = static_cast<float*>(var);
+        sstr << q[0] << " " << q[1] << " " << q[2] << " " << q[3];
+        break;
+      }
+    case TW_TYPE_COLOR4F:
+      {
+        sstr << "TW_TYPE_COLOR4F" << " ";
+        float * c = static_cast<float*>(var);
+        sstr << c[0] << " " << c[1] << " " << c[2] << " " << c[3];
+        break;
+      }
+    case TW_TYPE_COLOR3F:
+      {
+        sstr << "TW_TYPE_COLOR3F" << " ";
+        float * c = static_cast<float*>(var);
+        sstr << c[0] << " " << c[1] << " " << c[2];
+        break;
+      }
+    case TW_TYPE_DIR3F:
+      {
+        sstr << "TW_TYPE_DIR3F" << " ";
+        float * d = static_cast<float*>(var);
+        sstr << d[0] << " " << d[1] << " " << d[2];
+        break;
+      }
+    case TW_TYPE_BOOL32:
+      {
+        sstr << "TW_TYPE_BOOL32" << " ";
+        sstr << *(static_cast<int*>(var));
+        break;
+      }
+    case TW_TYPE_INT32:
+      {
+        sstr << "TW_TYPE_INT32" << " ";
+        sstr << *(static_cast<int*>(var));
+        break;
+      }
+    case TW_TYPE_FLOAT:
+      {
+        sstr << "TW_TYPE_FLOAT" << " ";
+        sstr << *(static_cast<float*>(var));
+        break;
+      }
+    case TW_TYPE_DOUBLE:
+      {
+        sstr << "TW_TYPE_DOUBLE" << " ";
+        sstr << *(static_cast<double*>(var));
+        break;
+      }
+    default:
+      {
+        std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > >::iterator iter = 
+          ReTw_custom_types.find(type);
+        if(iter != ReTw_custom_types.end())
+        {
+          sstr << (*iter).second.first << " ";
+          int enum_val = *(static_cast<int*>(var));
+          // try find display name for enum value
+          std::vector<TwEnumVal>::iterator eit = (*iter).second.second.begin();
+          bool found = false;
+          for(;eit<(*iter).second.second.end();eit++)
+          {
+            if(enum_val == eit->Value)
+            {
+              sstr << eit->Label;
+              found = true;
+              break;
+            }
+          }
+          if(!found)
+          {
+            sstr << "ERROR_ENUM_VALUE_NOT_DEFINED";
+          }
+        }else
+        {
+          sstr << "ERROR_TYPE_NOT_SUPPORTED";
+        }
+        break;
+      }
+  }
+  return sstr.str();
+}
+
+bool igl::ReTwBar::load(const char *file_name)
+{
+  FILE * fp;
+  fp = fopen(file_name,"r");
+
+  if(fp == NULL)
+  {
+    printf("ERROR: not able to open %s for reading...\n",file_name);
+    return false;
+  }
+
+  // go through file line by line
+  char line[MAX_LINE];
+  bool still_comments;
+  char name[MAX_WORD];
+  char type_str[MAX_WORD];
+  char value_str[MAX_WORD];
+
+
+  // line number
+  int j = 0;
+  bool finished = false;
+  while(true)
+  {
+    // Eat comments
+    still_comments = true;
+    while(still_comments)
+    {
+      if(fgets(line,MAX_LINE,fp) == NULL)
+      {
+        finished = true;
+        break;
+      }
+      // Blank lines and lines that begin with # are comments
+      still_comments = (line[0] == '#' || line[0] == '\n');
+      j++;
+    }
+    if(finished)
+    {
+      break;
+    }
+
+    sscanf(line,"%[^:]: %s %[^\n]",name,type_str,value_str);
+    //printf("%s: %s %s\n",name, type_str,value_str);
+
+    TwType type;
+    if(!type_from_string(type_str,type))
+    {
+      printf("ERROR: %s type not found... Skipping...\n",type_str);
+      continue;
+    }
+    set_value_from_string(name,type,value_str);
+
+  }
+
+  fclose(fp);
+  
+  // everything succeeded
+  return true;
+}
+
+bool igl::ReTwBar::type_from_string(const char *type_str, TwType & type)
+{
+  // first check default types
+  for(int j = 0; j < RETW_NUM_DEFAULT_TYPE_STRINGS; j++)
+  {
+    if(strcmp(type_str,ReTwDefaultTypeStrings[j].type_str) == 0)
+    {
+      type = ReTwDefaultTypeStrings[j].type;
+      return true;
+      break;
+    }
+  }
+
+  // then check custom types
+  std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > >::iterator iter = 
+    ReTw_custom_types.begin();
+  for(;iter != ReTw_custom_types.end(); iter++)
+  {
+    if(strcmp((*iter).second.first,type_str)==0)
+    {
+      type = (*iter).first;
+      return true;
+    }
+  }
+  return false;
+}
+
+bool igl::ReTwBar::set_value_from_string(
+  const char * name, 
+  TwType type, 
+  const char * value_str)
+{
+  void * value = NULL;
+  // possible value slots
+  int i;
+  float v;
+  double dv;
+  float f[4];
+  bool b;
+
+  // First try to get value from default types
+  switch(type)
+  {
+    case TW_TYPE_BOOLCPP:
+      {
+        int ib;
+        if(sscanf(value_str," %d",&ib) == 1)
+        {
+          b = ib!=0;
+          value = &b;
+        }else
+        {
+          printf("ERROR: Bad value format...\n");
+          return false;
+        }
+        break;
+      }
+    case TW_TYPE_QUAT4F:
+    case TW_TYPE_COLOR4F:
+      {
+        if(sscanf(value_str," %f %f %f %f",&f[0],&f[1],&f[2],&f[3]) == 4)
+        {
+          value = &f;
+        }else
+        {
+          printf("ERROR: Bad value format...\n");
+          return false;
+        }
+        break;
+      }
+    case TW_TYPE_COLOR3F:
+    case TW_TYPE_DIR3F:
+      {
+        if(sscanf(value_str," %f %f %f",&f[0],&f[1],&f[2]) == 3)
+        {
+          value = &f;
+        }else
+        {
+          printf("ERROR: Bad value format...\n");
+          return false;
+        }
+        break;
+      }
+    case TW_TYPE_BOOL32:
+    case TW_TYPE_INT32:
+      {
+        if(sscanf(value_str," %d",&i) == 1)
+        {
+          value = &i;
+        }else
+        {
+          printf("ERROR: Bad value format...\n");
+          return false;
+        }
+        break;
+      }
+    case TW_TYPE_FLOAT:
+      {
+        if(sscanf(value_str," %f",&v) == 1)
+        {
+          value = &v;
+        }else
+        {
+          printf("ERROR: Bad value format...\n");
+          return false;
+        }
+        break;
+      }
+    case TW_TYPE_DOUBLE:
+      {
+        if(sscanf(value_str," %lf",&dv) == 1)
+        {
+          value = &dv;
+        }else
+        {
+          printf("ERROR: Bad value format...\n");
+          return false;
+        }
+        break;
+      }
+    default:
+      // Try to find type in custom enum types
+      std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > >::iterator iter = 
+        ReTw_custom_types.find(type);
+      if(iter != ReTw_custom_types.end())
+      {
+        std::vector<TwEnumVal>::iterator eit = (*iter).second.second.begin();
+        bool found = false;
+        for(;eit<(*iter).second.second.end();eit++)
+        {
+          if(strcmp(value_str,eit->Label) == 0)
+          {
+            i = eit->Value;
+            value = &i;
+            found = true;
+            break;
+          }
+        }
+        if(!found)
+        {
+          printf("ERROR_ENUM_VALUE_NOT_DEFINED");
+        }
+      }else
+      {
+        printf("ERROR_TYPE_NOT_SUPPORTED\n");
+      }
+
+      break;
+  }
+
+
+  // Find variable based on name
+  // First look in RW items
+  bool item_found = false;
+  for(
+    std::vector<ReTwRWItem>::iterator it = rw_items.begin(); 
+    it != rw_items.end(); 
+    it++)
+  {
+    if(strcmp(it->name,name) == 0)
+    {
+      void * var = it->var;
+      switch(type)
+      {
+        case TW_TYPE_BOOLCPP:
+          {
+            bool * bvar = static_cast<bool*>(var);
+            bool * bvalue = static_cast<bool*>(value);
+            *bvar = *bvalue;
+            break;
+          }
+        case TW_TYPE_QUAT4F:
+        case TW_TYPE_COLOR4F:
+          {
+            float * fvar = static_cast<float*>(var);
+            float * fvalue = static_cast<float*>(value);
+            fvar[0] = fvalue[0];
+            fvar[1] = fvalue[1];
+            fvar[2] = fvalue[2];
+            fvar[3] = fvalue[3];
+            break;
+          }
+        case TW_TYPE_COLOR3F:
+        case TW_TYPE_DIR3F:
+          {
+            float * fvar = static_cast<float*>(var);
+            float * fvalue = static_cast<float*>(value);
+            fvar[0] = fvalue[0];
+            fvar[1] = fvalue[1];
+            fvar[2] = fvalue[2];
+            break;
+          }
+        case TW_TYPE_BOOL32:
+        case TW_TYPE_INT32:
+          {
+            int * ivar = static_cast<int*>(var);
+            int * ivalue = static_cast<int*>(value);
+            *ivar = *ivalue;
+            break;
+          }
+        case TW_TYPE_FLOAT:
+          {
+            float * fvar = static_cast<float*>(var);
+            float * fvalue = static_cast<float*>(value);
+            *fvar = *fvalue;
+            break;
+          }
+        case TW_TYPE_DOUBLE:
+          {
+            double * dvar =   static_cast<double*>(var);
+            double * fvalue = static_cast<double*>(value);
+            *dvar = *fvalue;
+            break;
+          }
+        default:
+          // Try to find type in custom enum types
+          std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > >::iterator iter = 
+            ReTw_custom_types.find(type);
+          if(iter != ReTw_custom_types.end())
+          {
+            int * ivar = static_cast<int*>(var);
+            std::vector<TwEnumVal>::iterator eit = (*iter).second.second.begin();
+            bool found = false;
+            for(;eit<(*iter).second.second.end();eit++)
+            {
+              if(strcmp(value_str,eit->Label) == 0)
+              {
+                *ivar = eit->Value;
+                found = true;
+                break;
+              }
+            }
+            if(!found)
+            {
+              printf("ERROR_ENUM_VALUE_NOT_DEFINED");
+            }
+          }else
+          {
+            printf("ERROR_TYPE_NOT_SUPPORTED\n");
+          }
+          break;
+      }
+      item_found = true;
+      break;
+    }
+  }
+
+  // Try looking in CB items
+  if(!item_found)
+  {
+    for(
+      std::vector<ReTwCBItem>::iterator it = cb_items.begin(); 
+      it != cb_items.end(); 
+      it++)
+    {
+      if(strcmp(it->name,name) == 0)
+      {
+        it->setCallback(value,it->clientData);
+        item_found = true;
+        break;
+      }
+    }
+  }
+
+  if(!item_found)
+  {
+    printf("ERROR: item not found\n");
+  }
+  return true;
+}

ReAntTweakBar.h

@@ -1,5 +1,6 @@
 #ifndef IGL_REANTTWEAKBAR_H
 #define IGL_REANTTWEAKBAR_H
+#include "igl_inline.h"
 // ReAntTweakBar is a minimal wrapper for the AntTweakBar library that allows
 // "bars" to be saved and load from disk. Changing your existing app that users
 // AntTweakBar to use ReAntTweakBar is trivial.
@@ -188,683 +189,8 @@ namespace igl
 //TW_API int      TW_CALL TwRemoveVar(TwBar *bar, const char *name);
 //TW_API int      TW_CALL TwRemoveAllVars(TwBar *bar);
 
-// Implementation
-
-#include <cstdio>
-#include <sstream>
-#include <map>
-
-#define MAX_CB_VAR_SIZE 10
-// Max line size for reading files
-#define MAX_LINE 1000
-#define MAX_WORD 100
-
-// GLOBAL WRAPPERS
-namespace igl
-{
-  std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > > ReTw_custom_types;
-}
-
-TwType igl::ReTwDefineEnum(
-  const char *name, 
-  const TwEnumVal *enumValues, 
-  unsigned int nbValues)
-{
-  // copy enum valus into vector
-  std::vector<TwEnumVal> enum_vals;
-  enum_vals.resize(nbValues);
-  for(unsigned int j = 0; j<nbValues;j++)
-  {
-    enum_vals[j] = enumValues[j];
-  }
-  TwType type = TwDefineEnum(name,enumValues,nbValues);
-
-  ReTw_custom_types[type] = 
-    std::pair<const char *,std::vector<TwEnumVal> >(name,enum_vals);
-  return type;
-}
-
-namespace igl
-{
-  struct ReTwTypeString
-  {
-    TwType type;
-    const char * type_str;
-  };
-
-  #define RETW_NUM_DEFAULT_TYPE_STRINGS 23
-  ReTwTypeString ReTwDefaultTypeStrings[RETW_NUM_DEFAULT_TYPE_STRINGS] = 
-  {
-    {TW_TYPE_UNDEF,"TW_TYPE_UNDEF"},
-    {TW_TYPE_BOOLCPP,"TW_TYPE_BOOLCPP"},
-    {TW_TYPE_BOOL8,"TW_TYPE_BOOL8"},
-    {TW_TYPE_BOOL16,"TW_TYPE_BOOL16"},
-    {TW_TYPE_BOOL32,"TW_TYPE_BOOL32"},
-    {TW_TYPE_CHAR,"TW_TYPE_CHAR"},
-    {TW_TYPE_INT8,"TW_TYPE_INT8"},
-    {TW_TYPE_UINT8,"TW_TYPE_UINT8"},
-    {TW_TYPE_INT16,"TW_TYPE_INT16"},
-    {TW_TYPE_UINT16,"TW_TYPE_UINT16"},
-    {TW_TYPE_INT32,"TW_TYPE_INT32"},
-    {TW_TYPE_UINT32,"TW_TYPE_UINT32"},
-    {TW_TYPE_FLOAT,"TW_TYPE_FLOAT"},
-    {TW_TYPE_DOUBLE,"TW_TYPE_DOUBLE"},
-    {TW_TYPE_COLOR32,"TW_TYPE_COLOR32"},
-    {TW_TYPE_COLOR3F,"TW_TYPE_COLOR3F"},
-    {TW_TYPE_COLOR4F,"TW_TYPE_COLOR4F"},
-    {TW_TYPE_CDSTRING,"TW_TYPE_CDSTRING"},
-    {TW_TYPE_STDSTRING,"TW_TYPE_STDSTRING"},
-    {TW_TYPE_QUAT4F,"TW_TYPE_QUAT4F"},
-    {TW_TYPE_QUAT4D,"TW_TYPE_QUAT4D"},
-    {TW_TYPE_DIR3F,"TW_TYPE_DIR3F"},
-    {TW_TYPE_DIR3D,"TW_TYPE_DIR3D"}
-  };
-}
-
-
-// BAR WRAPPERS
-void igl::ReTwBar::TwNewBar(const char *barName)
-{
-  // double colon without anything in front of it means look for this in the
-  // global namespace... I hope...
-  this->bar = ::TwNewBar(barName);
-}
-
-int igl::ReTwBar::TwAddVarRW(
-  const char *name, 
-  TwType type, 
-  void *var, 
-  const char *def,
-  const bool record)
-{
-  int ret = ::TwAddVarRW(this->bar,name,type,var,def);
-  if(ret && record)
-  {
-    rw_items.push_back(ReTwRWItem(name,type,var));
-  }
-  return ret;
-}
-
-int igl::ReTwBar::TwAddVarCB(
-  const char *name, 
-  TwType type, 
-  TwSetVarCallback setCallback, 
-  TwGetVarCallback getCallback, 
-  void *clientData, 
-  const char *def,
-  const bool record)
-{
-  int ret = 
-    ::TwAddVarCB(this->bar,name,type,setCallback,getCallback,clientData,def);
-  if(ret && record)
-  {
-    cb_items.push_back(ReTwCBItem(name,type,setCallback,getCallback,clientData));
-  }
-  return ret;
-}
-
-int igl::ReTwBar::TwAddVarRO(
-  const char *name, 
-  TwType type, 
-  void *var, 
-  const char *def)
-{
-  int ret = ::TwAddVarRO(this->bar,name,type,var,def);
-  // Read only variables are not recorded
-  //if(ret)
-  //{
-  //  rw_items.push_back(ReTwRWItem(name,type,var));
-  //}
-  return ret;
-}
-
-int igl::ReTwBar::TwAddButton(
-  const char *name, 
-  TwButtonCallback buttonCallback, 
-  void *clientData, 
-  const char *def)
-{
-  int ret = 
-    ::TwAddButton(this->bar,name,buttonCallback,clientData,def);
-  // buttons are not recorded
-  //if(ret)
-  //{
-  //  cb_items.push_back(ReTwCBItem(name,type,setCallback,getCallback,clientData));
-  //}
-  return ret;
-}
-
-int igl::ReTwBar::TwSetParam(
-  const char *varName, 
-  const char *paramName, 
-  TwParamValueType paramValueType, 
-  unsigned int inValueCount, 
-  const void *inValues)
-{
-  // For now just pass these along
-  return 
-    ::TwSetParam(
-      this->bar,
-      varName,
-      paramName,
-      paramValueType,
-      inValueCount,
-      inValues);
-}
-
-int igl::ReTwBar::TwGetParam(
-  const char *varName, 
-  const char *paramName, 
-  TwParamValueType paramValueType, 
-  unsigned int outValueMaxCount,
-  void *outValues)
-{
-  return 
-    ::TwGetParam(
-      this->bar,
-      varName,
-      paramName,
-      paramValueType,
-      outValueMaxCount,
-      outValues);
-}
-
-int igl::ReTwBar::TwRefreshBar()
-{
-  return ::TwRefreshBar(this->bar);
-}
-bool igl::ReTwBar::save(const char *file_name)
-{
-  FILE * fp;
-  if(file_name == NULL)
-  {
-    fp = stdout;
-  }else
-  {
-    fp = fopen(file_name,"w");
-  }
-
-  if(fp == NULL)
-  {
-    printf("ERROR: not able to open %s for writing...\n",file_name);
-    return false;
-  }
-
-  // Print all RW variables
-  for(
-    std::vector<ReTwRWItem>::iterator it = rw_items.begin(); 
-    it != rw_items.end(); 
-    it++)
-  {
-    const char * name = (*it).name;
-    TwType type = (*it).type;
-    void * var = (*it).var;
-    fprintf(fp,"%s: %s\n",
-      name,
-      get_value_as_string(var,type).c_str());
-  }
-
-  char var[MAX_CB_VAR_SIZE];
-  // Print all CB variables
-  for(
-    std::vector<ReTwCBItem>::iterator it = cb_items.begin(); 
-    it != cb_items.end(); 
-    it++)
-  {
-    const char * name = it->name;
-    TwType type = it->type;
-    //TwSetVarCallback setCallback = it->setCallback;
-    TwGetVarCallback getCallback = it->getCallback;
-    void * clientData = it->clientData;
-    // I'm not sure how to do what I want to do. getCallback needs to be sure
-    // that it can write to var. So var needs to point to a valid and big
-    // enough chunk of memory
-    getCallback(var,clientData);
-    fprintf(fp,"%s: %s\n",
-      name,
-      get_value_as_string(var,type).c_str());
-  }
-
-  fprintf(fp,"\n");
-
-  if(file_name != NULL)
-  {
-    fclose(fp);
-  }
-  // everything succeeded
-  return true;
-}
-
-std::string igl::ReTwBar::get_value_as_string(
-  void * var, 
-  TwType type)
-{
-  std::stringstream sstr;
-  switch(type)
-  {
-    case TW_TYPE_BOOLCPP:
-      {
-        sstr << "TW_TYPE_BOOLCPP" << " ";
-        sstr << *(static_cast<bool*>(var));
-        break;
-      }
-    case TW_TYPE_QUAT4F:
-      {
-        sstr << "TW_TYPE_QUAT4F" << " ";
-        // Q: Why does casting to float* work? shouldn't I have to cast to
-        // float**?
-        float * q = static_cast<float*>(var);
-        sstr << q[0] << " " << q[1] << " " << q[2] << " " << q[3];
-        break;
-      }
-    case TW_TYPE_COLOR4F:
-      {
-        sstr << "TW_TYPE_COLOR4F" << " ";
-        float * c = static_cast<float*>(var);
-        sstr << c[0] << " " << c[1] << " " << c[2] << " " << c[3];
-        break;
-      }
-    case TW_TYPE_COLOR3F:
-      {
-        sstr << "TW_TYPE_COLOR3F" << " ";
-        float * c = static_cast<float*>(var);
-        sstr << c[0] << " " << c[1] << " " << c[2];
-        break;
-      }
-    case TW_TYPE_DIR3F:
-      {
-        sstr << "TW_TYPE_DIR3F" << " ";
-        float * d = static_cast<float*>(var);
-        sstr << d[0] << " " << d[1] << " " << d[2];
-        break;
-      }
-    case TW_TYPE_BOOL32:
-      {
-        sstr << "TW_TYPE_BOOL32" << " ";
-        sstr << *(static_cast<int*>(var));
-        break;
-      }
-    case TW_TYPE_INT32:
-      {
-        sstr << "TW_TYPE_INT32" << " ";
-        sstr << *(static_cast<int*>(var));
-        break;
-      }
-    case TW_TYPE_FLOAT:
-      {
-        sstr << "TW_TYPE_FLOAT" << " ";
-        sstr << *(static_cast<float*>(var));
-        break;
-      }
-    case TW_TYPE_DOUBLE:
-      {
-        sstr << "TW_TYPE_DOUBLE" << " ";
-        sstr << *(static_cast<double*>(var));
-        break;
-      }
-    default:
-      {
-        std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > >::iterator iter = 
-          ReTw_custom_types.find(type);
-        if(iter != ReTw_custom_types.end())
-        {
-          sstr << (*iter).second.first << " ";
-          int enum_val = *(static_cast<int*>(var));
-          // try find display name for enum value
-          std::vector<TwEnumVal>::iterator eit = (*iter).second.second.begin();
-          bool found = false;
-          for(;eit<(*iter).second.second.end();eit++)
-          {
-            if(enum_val == eit->Value)
-            {
-              sstr << eit->Label;
-              found = true;
-              break;
-            }
-          }
-          if(!found)
-          {
-            sstr << "ERROR_ENUM_VALUE_NOT_DEFINED";
-          }
-        }else
-        {
-          sstr << "ERROR_TYPE_NOT_SUPPORTED";
-        }
-        break;
-      }
-  }
-  return sstr.str();
-}
-
-bool igl::ReTwBar::load(const char *file_name)
-{
-  FILE * fp;
-  fp = fopen(file_name,"r");
-
-  if(fp == NULL)
-  {
-    printf("ERROR: not able to open %s for reading...\n",file_name);
-    return false;
-  }
-
-  // go through file line by line
-  char line[MAX_LINE];
-  bool still_comments;
-  char name[MAX_WORD];
-  char type_str[MAX_WORD];
-  char value_str[MAX_WORD];
-
-
-  // line number
-  int j = 0;
-  bool finished = false;
-  while(true)
-  {
-    // Eat comments
-    still_comments = true;
-    while(still_comments)
-    {
-      if(fgets(line,MAX_LINE,fp) == NULL)
-      {
-        finished = true;
-        break;
-      }
-      // Blank lines and lines that begin with # are comments
-      still_comments = (line[0] == '#' || line[0] == '\n');
-      j++;
-    }
-    if(finished)
-    {
-      break;
-    }
-
-    sscanf(line,"%[^:]: %s %[^\n]",name,type_str,value_str);
-    //printf("%s: %s %s\n",name, type_str,value_str);
-
-    TwType type;
-    if(!type_from_string(type_str,type))
-    {
-      printf("ERROR: %s type not found... Skipping...\n",type_str);
-      continue;
-    }
-    set_value_from_string(name,type,value_str);
-
-  }
-
-  fclose(fp);
-  
-  // everything succeeded
-  return true;
-}
-
-bool igl::ReTwBar::type_from_string(const char *type_str, TwType & type)
-{
-  // first check default types
-  for(int j = 0; j < RETW_NUM_DEFAULT_TYPE_STRINGS; j++)
-  {
-    if(strcmp(type_str,ReTwDefaultTypeStrings[j].type_str) == 0)
-    {
-      type = ReTwDefaultTypeStrings[j].type;
-      return true;
-      break;
-    }
-  }
-
-  // then check custom types
-  std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > >::iterator iter = 
-    ReTw_custom_types.begin();
-  for(;iter != ReTw_custom_types.end(); iter++)
-  {
-    if(strcmp((*iter).second.first,type_str)==0)
-    {
-      type = (*iter).first;
-      return true;
-    }
-  }
-  return false;
-}
-
-bool igl::ReTwBar::set_value_from_string(
-  const char * name, 
-  TwType type, 
-  const char * value_str)
-{
-  void * value = NULL;
-  // possible value slots
-  int i;
-  float v;
-  double dv;
-  float f[4];
-  bool b;
-
-  // First try to get value from default types
-  switch(type)
-  {
-    case TW_TYPE_BOOLCPP:
-      {
-        int ib;
-        if(sscanf(value_str," %d",&ib) == 1)
-        {
-          b = ib!=0;
-          value = &b;
-        }else
-        {
-          printf("ERROR: Bad value format...\n");
-          return false;
-        }
-        break;
-      }
-    case TW_TYPE_QUAT4F:
-    case TW_TYPE_COLOR4F:
-      {
-        if(sscanf(value_str," %f %f %f %f",&f[0],&f[1],&f[2],&f[3]) == 4)
-        {
-          value = &f;
-        }else
-        {
-          printf("ERROR: Bad value format...\n");
-          return false;
-        }
-        break;
-      }
-    case TW_TYPE_COLOR3F:
-    case TW_TYPE_DIR3F:
-      {
-        if(sscanf(value_str," %f %f %f",&f[0],&f[1],&f[2]) == 3)
-        {
-          value = &f;
-        }else
-        {
-          printf("ERROR: Bad value format...\n");
-          return false;
-        }
-        break;
-      }
-    case TW_TYPE_BOOL32:
-    case TW_TYPE_INT32:
-      {
-        if(sscanf(value_str," %d",&i) == 1)
-        {
-          value = &i;
-        }else
-        {
-          printf("ERROR: Bad value format...\n");
-          return false;
-        }
-        break;
-      }
-    case TW_TYPE_FLOAT:
-      {
-        if(sscanf(value_str," %f",&v) == 1)
-        {
-          value = &v;
-        }else
-        {
-          printf("ERROR: Bad value format...\n");
-          return false;
-        }
-        break;
-      }
-    case TW_TYPE_DOUBLE:
-      {
-        if(sscanf(value_str," %lf",&dv) == 1)
-        {
-          value = &dv;
-        }else
-        {
-          printf("ERROR: Bad value format...\n");
-          return false;
-        }
-        break;
-      }
-    default:
-      // Try to find type in custom enum types
-      std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > >::iterator iter = 
-        ReTw_custom_types.find(type);
-      if(iter != ReTw_custom_types.end())
-      {
-        std::vector<TwEnumVal>::iterator eit = (*iter).second.second.begin();
-        bool found = false;
-        for(;eit<(*iter).second.second.end();eit++)
-        {
-          if(strcmp(value_str,eit->Label) == 0)
-          {
-            i = eit->Value;
-            value = &i;
-            found = true;
-            break;
-          }
-        }
-        if(!found)
-        {
-          printf("ERROR_ENUM_VALUE_NOT_DEFINED");
-        }
-      }else
-      {
-        printf("ERROR_TYPE_NOT_SUPPORTED\n");
-      }
-
-      break;
-  }
-
-
-  // Find variable based on name
-  // First look in RW items
-  bool item_found = false;
-  for(
-    std::vector<ReTwRWItem>::iterator it = rw_items.begin(); 
-    it != rw_items.end(); 
-    it++)
-  {
-    if(strcmp(it->name,name) == 0)
-    {
-      void * var = it->var;
-      switch(type)
-      {
-        case TW_TYPE_BOOLCPP:
-          {
-            bool * bvar = static_cast<bool*>(var);
-            bool * bvalue = static_cast<bool*>(value);
-            *bvar = *bvalue;
-            break;
-          }
-        case TW_TYPE_QUAT4F:
-        case TW_TYPE_COLOR4F:
-          {
-            float * fvar = static_cast<float*>(var);
-            float * fvalue = static_cast<float*>(value);
-            fvar[0] = fvalue[0];
-            fvar[1] = fvalue[1];
-            fvar[2] = fvalue[2];
-            fvar[3] = fvalue[3];
-            break;
-          }
-        case TW_TYPE_COLOR3F:
-        case TW_TYPE_DIR3F:
-          {
-            float * fvar = static_cast<float*>(var);
-            float * fvalue = static_cast<float*>(value);
-            fvar[0] = fvalue[0];
-            fvar[1] = fvalue[1];
-            fvar[2] = fvalue[2];
-            break;
-          }
-        case TW_TYPE_BOOL32:
-        case TW_TYPE_INT32:
-          {
-            int * ivar = static_cast<int*>(var);
-            int * ivalue = static_cast<int*>(value);
-            *ivar = *ivalue;
-            break;
-          }
-        case TW_TYPE_FLOAT:
-          {
-            float * fvar = static_cast<float*>(var);
-            float * fvalue = static_cast<float*>(value);
-            *fvar = *fvalue;
-            break;
-          }
-        case TW_TYPE_DOUBLE:
-          {
-            double * dvar =   static_cast<double*>(var);
-            double * fvalue = static_cast<double*>(value);
-            *dvar = *fvalue;
-            break;
-          }
-        default:
-          // Try to find type in custom enum types
-          std::map<TwType,std::pair<const char *,std::vector<TwEnumVal> > >::iterator iter = 
-            ReTw_custom_types.find(type);
-          if(iter != ReTw_custom_types.end())
-          {
-            int * ivar = static_cast<int*>(var);
-            std::vector<TwEnumVal>::iterator eit = (*iter).second.second.begin();
-            bool found = false;
-            for(;eit<(*iter).second.second.end();eit++)
-            {
-              if(strcmp(value_str,eit->Label) == 0)
-              {
-                *ivar = eit->Value;
-                found = true;
-                break;
-              }
-            }
-            if(!found)
-            {
-              printf("ERROR_ENUM_VALUE_NOT_DEFINED");
-            }
-          }else
-          {
-            printf("ERROR_TYPE_NOT_SUPPORTED\n");
-          }
-          break;
-      }
-      item_found = true;
-      break;
-    }
-  }
-
-  // Try looking in CB items
-  if(!item_found)
-  {
-    for(
-      std::vector<ReTwCBItem>::iterator it = cb_items.begin(); 
-      it != cb_items.end(); 
-      it++)
-    {
-      if(strcmp(it->name,name) == 0)
-      {
-        it->setCallback(value,it->clientData);
-        item_found = true;
-        break;
-      }
-    }
-  }
+#ifdef IGL_HEADER_ONLY
+#  include "ReAntTweakBar.cpp"
+#endif
 
-  if(!item_found)
-  {
-    printf("ERROR: item not found\n");
-  }
-  return true;
-}
 #endif

adjacency_list.cpp

@@ -0,0 +1,124 @@
+#include "adjacency_list.h"
+
+#include "verbose.h"
+
+template <typename T, typename M>
+IGL_INLINE void igl::adjacency_list(
+    const M & F, 
+    std::vector<std::vector<T> >& A,
+    bool sorted)
+{
+  A.clear(); 
+  A.resize(F.maxCoeff()+1);
+  
+  // Loop over faces
+  for(int i = 0;i<F.rows();i++)
+  {
+    // Loop over this face
+    for(int j = 0;j<F.cols();j++)
+    {
+      // Get indices of edge: s --> d
+      int s = F(i,j);
+      int d = F(i,(j+1)%F.cols());
+      A.at(s).push_back(d);
+      A.at(d).push_back(s);
+    }
+  }
+  
+  // Remove duplicates
+  for(int i=0; i<(int)A.size();++i)
+  {
+    std::sort(A[i].begin(), A[i].end());
+    A[i].erase(std::unique(A[i].begin(), A[i].end()), A[i].end());
+  }
+  
+  // If needed, sort every VV
+  if (sorted)
+  {
+    // Loop over faces
+    
+    // for every vertex v store a set of ordered edges not incident to v that belongs to triangle incident on v.
+    std::vector<std::vector<std::vector<int> > > SR; 
+    SR.resize(A.size());
+    
+    for(int i = 0;i<F.rows();i++)
+    {
+      // Loop over this face
+      for(int j = 0;j<F.cols();j++)
+      {
+        // Get indices of edge: s --> d
+        int s = F(i,j);
+        int d = F(i,(j+1)%F.cols());
+        // Get index of opposing vertex v
+        int v = F(i,(j+2)%F.cols());
+        
+        std::vector<int> e(2);
+        e[0] = d;
+        e[1] = v;
+        SR[s].push_back(e);
+      }
+    }
+    
+    for(int v=0; v<(int)SR.size();++v)
+    {
+      std::vector<int>& vv = A.at(v);
+      std::vector<std::vector<int> >& sr = SR[v];
+      
+      std::vector<std::vector<int> > pn = sr;
+      
+      // Compute previous/next for every element in sr
+      for(int i=0;i<(int)sr.size();++i)
+      {
+        int a = sr[i][0];
+        int b = sr[i][1];
+        
+        // search for previous
+        int p = -1;
+        for(int j=0;j<(int)sr.size();++j)
+          if(sr[j][1] == a)
+            p = j;
+        pn[i][0] = p;
+        
+        // search for next
+        int n = -1;
+        for(int j=0;j<(int)sr.size();++j)
+          if(sr[j][0] == b)
+            n = j;
+        pn[i][1] = n;
+        
+      }
+      
+      // assume manifoldness (look for beginning of a single chain)
+      int c = 0;
+      for(int j=0; j<=(int)sr.size();++j)
+        if (pn[c][0] != -1)
+          c = pn[c][0];
+      
+      if (pn[c][0] == -1) // border case
+      {
+        // finally produce the new vv relation
+        for(int j=0; j<(int)sr.size();++j)
+        {
+          vv[j] = sr[c][0];
+          if (pn[c][1] != -1)
+            c = pn[c][1];
+        }
+        vv.back() = sr[c][1];
+      }
+      else
+      {
+        // finally produce the new vv relation
+        for(int j=0; j<(int)sr.size();++j)
+        {
+          vv[j] = sr[c][0];
+          
+          c = pn[c][1];
+        }
+      }
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

adjacency_list.h

@@ -1,5 +1,6 @@
 #ifndef IGL_ADJACENCY_LIST_H
 #define IGL_ADJACENCY_LIST_H
+#include "igl_inline.h"
 
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
@@ -23,130 +24,14 @@ namespace igl
   //
   // See also: edges, cotmatrix, diag
   template <typename T, typename M>
-  inline void adjacency_list(
+  IGL_INLINE void adjacency_list(
     const M & F, 
     std::vector<std::vector<T> >& A,
     bool sorted = false);
 }
 
-// Implementation
-#include "verbose.h"
-
-template <typename T, typename M>
-inline void igl::adjacency_list(
-    const M & F, 
-    std::vector<std::vector<T> >& A,
-    bool sorted)
-{
-  A.clear(); 
-  A.resize(F.maxCoeff()+1);
-  
-  // Loop over faces
-  for(int i = 0;i<F.rows();i++)
-  {
-    // Loop over this face
-    for(int j = 0;j<F.cols();j++)
-    {
-      // Get indices of edge: s --> d
-      int s = F(i,j);
-      int d = F(i,(j+1)%F.cols());
-      A.at(s).push_back(d);
-      A.at(d).push_back(s);
-    }
-  }
-  
-  // Remove duplicates
-  for(int i=0; i<(int)A.size();++i)
-  {
-    std::sort(A[i].begin(), A[i].end());
-    A[i].erase(std::unique(A[i].begin(), A[i].end()), A[i].end());
-  }
-  
-  // If needed, sort every VV
-  if (sorted)
-  {
-    // Loop over faces
-    
-    // for every vertex v store a set of ordered edges not incident to v that belongs to triangle incident on v.
-    std::vector<std::vector<std::vector<int> > > SR; 
-    SR.resize(A.size());
-    
-    for(int i = 0;i<F.rows();i++)
-    {
-      // Loop over this face
-      for(int j = 0;j<F.cols();j++)
-      {
-        // Get indices of edge: s --> d
-        int s = F(i,j);
-        int d = F(i,(j+1)%F.cols());
-        // Get index of opposing vertex v
-        int v = F(i,(j+2)%F.cols());
-        
-        std::vector<int> e(2);
-        e[0] = d;
-        e[1] = v;
-        SR[s].push_back(e);
-      }
-    }
-    
-    for(int v=0; v<(int)SR.size();++v)
-    {
-      std::vector<int>& vv = A.at(v);
-      std::vector<std::vector<int> >& sr = SR[v];
-      
-      std::vector<std::vector<int> > pn = sr;
-      
-      // Compute previous/next for every element in sr
-      for(int i=0;i<(int)sr.size();++i)
-      {
-        int a = sr[i][0];
-        int b = sr[i][1];
-        
-        // search for previous
-        int p = -1;
-        for(int j=0;j<(int)sr.size();++j)
-          if(sr[j][1] == a)
-            p = j;
-        pn[i][0] = p;
-        
-        // search for next
-        int n = -1;
-        for(int j=0;j<(int)sr.size();++j)
-          if(sr[j][0] == b)
-            n = j;
-        pn[i][1] = n;
-        
-      }
-      
-      // assume manifoldness (look for beginning of a single chain)
-      int c = 0;
-      for(int j=0; j<=(int)sr.size();++j)
-        if (pn[c][0] != -1)
-          c = pn[c][0];
-      
-      if (pn[c][0] == -1) // border case
-      {
-        // finally produce the new vv relation
-        for(int j=0; j<(int)sr.size();++j)
-        {
-          vv[j] = sr[c][0];
-          if (pn[c][1] != -1)
-            c = pn[c][1];
-        }
-        vv.back() = sr[c][1];
-      }
-      else
-      {
-        // finally produce the new vv relation
-        for(int j=0; j<(int)sr.size();++j)
-        {
-          vv[j] = sr[c][0];
-          
-          c = pn[c][1];
-        }
-      }
-    }
-  }
-}
+#ifdef IGL_HEADER_ONLY
+#  include "adjacency_list.cpp"
+#endif
 
 #endif

adjacency_matrix.cpp

@@ -0,0 +1,33 @@
+#include "adjacency_matrix.h"
+
+#include "verbose.h"
+
+template <typename T>
+IGL_INLINE void igl::adjacency_matrix(
+  const Eigen::MatrixXi & F, 
+  Eigen::SparseMatrix<T>& A)
+{
+  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> 
+    dyn_A(F.maxCoeff()+1, F.maxCoeff()+1);
+  dyn_A.reserve(6*(F.maxCoeff()+1));
+
+  // Loop over faces
+  for(int i = 0;i<F.rows();i++)
+  {
+    // Loop over this face
+    for(int j = 0;j<F.cols();j++)
+    {
+      // Get indices of edge: s --> d
+      int s = F(i,j);
+      int d = F(i,(j+1)%F.cols());
+      dyn_A.coeffRef(s, d) = 1;
+      dyn_A.coeffRef(d, s) = 1;
+    }
+  }
+
+  A = Eigen::SparseMatrix<T>(dyn_A);
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

adjacency_matrix.h

@@ -1,5 +1,6 @@
 #ifndef IGL_ADJACENCY_MATRIX_H
 #define IGL_ADJACENCY_MATRIX_H
+#include "igl_inline.h"
 
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
@@ -31,38 +32,13 @@ namespace igl
   //
   // See also: edges, cotmatrix, diag
   template <typename T>
-  inline void adjacency_matrix(
+  IGL_INLINE void adjacency_matrix(
     const Eigen::MatrixXi & F, 
     Eigen::SparseMatrix<T>& A);
 }
 
-// Implementation
-#include "verbose.h"
-
-template <typename T>
-inline void igl::adjacency_matrix(
-  const Eigen::MatrixXi & F, 
-  Eigen::SparseMatrix<T>& A)
-{
-  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> 
-    dyn_A(F.maxCoeff()+1, F.maxCoeff()+1);
-  dyn_A.reserve(6*(F.maxCoeff()+1));
-
-  // Loop over faces
-  for(int i = 0;i<F.rows();i++)
-  {
-    // Loop over this face
-    for(int j = 0;j<F.cols();j++)
-    {
-      // Get indices of edge: s --> d
-      int s = F(i,j);
-      int d = F(i,(j+1)%F.cols());
-      dyn_A.coeffRef(s, d) = 1;
-      dyn_A.coeffRef(d, s) = 1;
-    }
-  }
-
-  A = Eigen::SparseMatrix<T>(dyn_A);
-}
+#ifdef IGL_HEADER_ONLY
+#  include "adjacency_matrix.cpp"
+#endif
 
 #endif

all_pairs_distances.cpp

@@ -0,0 +1,29 @@
+#include "all_pairs_distances.h"
+
+template <typename Mat>
+IGL_INLINE void igl::all_pairs_distances(
+  const Mat & V,
+  const Mat & U,
+  const bool squared,
+  Mat & D)
+{
+  // dimension should be the same
+  assert(V.cols() == U.cols());
+  // resize output
+  D.resize(V.rows(),U.rows());
+  for(int i = 0;i<V.rows();i++)
+  {
+    for(int j=0;j<U.rows();j++)
+    {
+      D(i,j) = (V.row(i)-U.row(j)).array().pow(2).sum();
+      if(!squared)
+      {
+        D(i,j) = sqrt(D(i,j));
+      }
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

all_pairs_distances.h

@@ -1,5 +1,6 @@
 #ifndef IGL_ALL_PAIRS_DISTANCES_H
 #define IGL_ALL_PAIRS_DISTANCES_H
+#include "igl_inline.h"
 
 namespace igl
 {
@@ -19,37 +20,15 @@ namespace igl
   //     squareed distance between V(i,:) and U(j,:)
   // 
   template <typename Mat>
-  inline void all_pairs_distances(
+  IGL_INLINE void all_pairs_distances(
     const Mat & V,
     const Mat & U,
     const bool squared, 
     Mat & D);
 }
 
-// Implementation
-
-template <typename Mat>
-inline void igl::all_pairs_distances(
-  const Mat & V,
-  const Mat & U,
-  const bool squared,
-  Mat & D)
-{
-  // dimension should be the same
-  assert(V.cols() == U.cols());
-  // resize output
-  D.resize(V.rows(),U.rows());
-  for(int i = 0;i<V.rows();i++)
-  {
-    for(int j=0;j<U.rows();j++)
-    {
-      D(i,j) = (V.row(i)-U.row(j)).array().pow(2).sum();
-      if(!squared)
-      {
-        D(i,j) = sqrt(D(i,j));
-      }
-    }
-  }
-}
+#ifdef IGL_HEADER_ONLY
+#  include "all_pairs_distances.cpp"
+#endif
 
 #endif

axis_angle_to_quat.cpp

@@ -0,0 +1,29 @@
+#include "axis_angle_to_quat.h"
+
+// http://www.antisphere.com/Wiki/tools:anttweakbar
+template <typename Q_type>
+IGL_INLINE void igl::axis_angle_to_quat(
+  const Q_type *axis, 
+  const Q_type angle,
+  Q_type *out)
+{
+    Q_type n = axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2];
+    if( fabs(n)>igl::EPS<Q_type>())
+    {
+        Q_type f = 0.5*angle;
+        out[3] = cos(f);
+        f = sin(f)/sqrt(n);
+        out[0] = axis[0]*f;
+        out[1] = axis[1]*f;
+        out[2] = axis[2]*f;
+    }
+    else
+    {
+        out[3] = 1.0;
+        out[0] = out[1] = out[2] = 0.0;
+    }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

axis_angle_to_quat.h

@@ -1,5 +1,6 @@
 #ifndef IGL_AXIS_ANGLE_TO_QUAT_H
 #define IGL_AXIS_ANGLE_TO_QUAT_H
+#include "igl_inline.h"
 
 #include "EPS.h"
 #include <cmath>
@@ -14,35 +15,14 @@ namespace igl
   // Outputs:
   //   quaternion
   template <typename Q_type>
-  inline void axis_angle_to_quat(
+  IGL_INLINE void axis_angle_to_quat(
     const Q_type *axis, 
     const Q_type angle,
     Q_type *out);
 }
 
-// Implementation
-// http://www.antisphere.com/Wiki/tools:anttweakbar
-template <typename Q_type>
-inline void igl::axis_angle_to_quat(
-  const Q_type *axis, 
-  const Q_type angle,
-  Q_type *out)
-{
-    Q_type n = axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2];
-    if( fabs(n)>igl::EPS<Q_type>())
-    {
-        Q_type f = 0.5*angle;
-        out[3] = cos(f);
-        f = sin(f)/sqrt(n);
-        out[0] = axis[0]*f;
-        out[1] = axis[1]*f;
-        out[2] = axis[2]*f;
-    }
-    else
-    {
-        out[3] = 1.0;
-        out[0] = out[1] = out[2] = 0.0;
-    }
-}
+#ifdef IGL_HEADER_ONLY
+#  include "axis_angle_to_quat.cpp"
+#endif
 
 #endif

basename.cpp

@@ -0,0 +1,31 @@
+#include "basename.h"
+
+#include <algorithm>
+
+IGL_INLINE std::string igl::basename(const std::string & path)
+{
+  if(path == "")
+  {
+    return std::string("");
+  }
+  // http://stackoverflow.com/questions/5077693/dirnamephp-similar-function-in-c
+  std::string::const_reverse_iterator last_slash =
+    std::find(
+      path.rbegin(), 
+      path.rend(), '/');
+  if( last_slash == path.rend() )
+  {
+    // No slashes found
+    return path;
+  }else if(1 == (last_slash.base() - path.begin()))
+  {
+    // Slash is first char
+    return std::string(path.begin()+1,path.end());
+  }else if(path.end() == last_slash.base() )
+  {
+    // Slash is last char
+    std::string redo = std::string(path.begin(),path.end()-1);
+    return igl::basename(redo);
+  }
+  return std::string(last_slash.base(),path.end());
+}

basename.h

@@ -1,5 +1,6 @@
 #ifndef IGL_BASENAME_H
 #define IGL_BASENAME_H
+#include "igl_inline.h"
 
 #include <string>
 
@@ -11,38 +12,11 @@ namespace igl
   // Returns string containing basename (see php's basename)
   //
   // See also: dirname, pathinfo
-  inline std::string basename(const std::string & path);
+  IGL_INLINE std::string basename(const std::string & path);
 }
 
-// Implementation
-#include <algorithm>
-
-inline std::string igl::basename(const std::string & path)
-{
-  if(path == "")
-  {
-    return std::string("");
-  }
-  // http://stackoverflow.com/questions/5077693/dirnamephp-similar-function-in-c
-  std::string::const_reverse_iterator last_slash =
-    std::find(
-      path.rbegin(), 
-      path.rend(), '/');
-  if( last_slash == path.rend() )
-  {
-    // No slashes found
-    return path;
-  }else if(1 == (last_slash.base() - path.begin()))
-  {
-    // Slash is first char
-    return std::string(path.begin()+1,path.end());
-  }else if(path.end() == last_slash.base() )
-  {
-    // Slash is last char
-    std::string redo = std::string(path.begin(),path.end()-1);
-    return igl::basename(redo);
-  }
-  return std::string(last_slash.base(),path.end());
-}
+#ifdef IGL_HEADER_ONLY
+#  include "basename.cpp"
 #endif
 
+#endif

canonical_quaternions.cpp

@@ -0,0 +1,14 @@
+#include "canonical_quaternions.h"
+
+template <> IGL_INLINE const float igl::CANONICAL_VIEW_QUAT<float>(int i, int j)
+{
+  return igl::CANONICAL_VIEW_QUAT_F[i][j];
+}
+template <> IGL_INLINE const double igl::CANONICAL_VIEW_QUAT<double>(int i, int j)
+{
+  return igl::CANONICAL_VIEW_QUAT_D[i][j];
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

canonical_quaternions.h

@@ -1,5 +1,6 @@
 #ifndef IGL_CANONICAL_QUATERNIONS_H
 #define IGL_CANONICAL_QUATERNIONS_H
+#include "igl_inline.h"
 // Define some canonical quaternions for floats and doubles
 // A Quaternion, q, is defined here as an arrays of four scalars (x,y,z,w),
 // such that q = x*i + y*j + z*k + w
@@ -94,8 +95,7 @@ namespace igl
       {             0,-SQRT_2_OVER_2, SQRT_2_OVER_2,             0},
       {           0.5,          -0.5,           0.5,          -0.5}
     };
-
-  const size_t NUM_CANONICAL_VIEW_QUAT = 24;
+#define NUM_CANONICAL_VIEW_QUAT 24
 
   // NOTE: I want to rather be able to return a Q_type[][] but C++ is not
   // making it easy. So instead I've written a per-element accessor
@@ -107,25 +107,18 @@ namespace igl
   //   j  index of coordinate in quaternion i
   // Returns values of CANONICAL_VIEW_QUAT_*[i][j]
   template <typename Q_type> 
-  inline const Q_type CANONICAL_VIEW_QUAT(size_t i, size_t j);
+  IGL_INLINE const Q_type CANONICAL_VIEW_QUAT(int i, int j);
   // Template specializations for float and double
   template <> 
-  inline const float CANONICAL_VIEW_QUAT<float>(size_t i, size_t j);
+  IGL_INLINE const float CANONICAL_VIEW_QUAT<float>(int i, int j);
   template <> 
-  inline const double CANONICAL_VIEW_QUAT<double>(size_t i, size_t j);
+  IGL_INLINE const double CANONICAL_VIEW_QUAT<double>(int i, int j);
 
 #  undef SQRT_2_OVER_2
 }
 
-// Implementation
-
-template <> inline const float igl::CANONICAL_VIEW_QUAT<float>(size_t i, size_t j)
-{
-  return igl::CANONICAL_VIEW_QUAT_F[i][j];
-}
-template <> inline const double igl::CANONICAL_VIEW_QUAT<double>(size_t i, size_t j)
-{
-  return igl::CANONICAL_VIEW_QUAT_D[i][j];
-}
+#ifdef IGL_HEADER_ONLY
+#  include "canonical_quaternions.cpp"
+#endif
 
 #endif

cat.cpp

@@ -0,0 +1,136 @@
+#include "cat.h"
+
+// Sparse matrices need to be handled carefully. Because C++ does not 
+// Template:
+//   Scalar  sparse matrix scalar type, e.g. double
+template <typename Scalar>
+IGL_INLINE void igl::cat(
+    const int dim, 
+    const Eigen::SparseMatrix<Scalar> & A, 
+    const Eigen::SparseMatrix<Scalar> & B, 
+    Eigen::SparseMatrix<Scalar> & C)
+{
+  assert(dim == 1 || dim == 2);
+  using namespace Eigen;
+  // Special case if B or A is empty
+  if(A.size() == 0)
+  {
+    C = B;
+    return;
+  }
+  if(B.size() == 0)
+  {
+    C = A;
+    return;
+  }
+
+  DynamicSparseMatrix<Scalar, RowMajor> dyn_C;
+  if(dim == 1)
+  {
+    assert(A.cols() == B.cols());
+    dyn_C.resize(A.rows()+B.rows(),A.cols());
+  }else if(dim == 2)
+  {
+    assert(A.rows() == B.rows());
+    dyn_C.resize(A.rows(),A.cols()+B.cols());
+  }else
+  {
+    fprintf(stderr,"cat.h: Error: Unsupported dimension %d\n",dim);
+  }
+
+  dyn_C.reserve(A.nonZeros()+B.nonZeros());
+
+  // Iterate over outside of A
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+    {
+      dyn_C.coeffRef(it.row(),it.col()) += it.value();
+    }
+  }
+
+  // Iterate over outside of B
+  for(int k=0; k<B.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
+    {
+      int r = (dim == 1 ? A.rows()+it.row() : it.row());
+      int c = (dim == 2 ? A.cols()+it.col() : it.col());
+      dyn_C.coeffRef(r,c) += it.value();
+    }
+  }
+
+  C = SparseMatrix<Scalar>(dyn_C);
+}
+
+template <typename Derived, class MatC>
+IGL_INLINE void igl::cat(
+  const int dim,
+  const Eigen::MatrixBase<Derived> & A, 
+  const Eigen::MatrixBase<Derived> & B,
+  MatC & C)
+{
+  assert(dim == 1 || dim == 2);
+  // Special case if B or A is empty
+  if(A.size() == 0)
+  {
+    C = B;
+    return;
+  }
+  if(B.size() == 0)
+  {
+    C = A;
+    return;
+  }
+
+  if(dim == 1)
+  {
+    assert(A.cols() == B.cols());
+    C.resize(A.rows()+B.rows(),A.cols());
+    C << A,B;
+  }else if(dim == 2)
+  {
+    assert(A.rows() == B.rows());
+    C.resize(A.rows(),A.cols()+B.cols());
+    C << A,B;
+  }else
+  {
+    fprintf(stderr,"cat.h: Error: Unsupported dimension %d\n",dim);
+  }
+}
+
+template <class Mat>
+IGL_INLINE Mat igl::cat(const int dim, const Mat & A, const Mat & B)
+{
+  assert(dim == 1 || dim == 2);
+  Mat C;
+  igl::cat(dim,A,B,C);
+  return C;
+}
+
+template <class Mat>
+IGL_INLINE void cat(const std::vector<std::vector< Mat > > & A, Mat & C)
+{
+  using namespace igl;
+  using namespace std;
+  // Start with empty matrix
+  C.resize(0,0);
+  for(typename vector<vector< Mat > >::const_iterator rit = A.begin(); rit != A.end(); rit++)
+  {
+    // Concatenate each row horizontally
+    // Start with empty matrix
+    Mat row(0,0);
+    for(typename vector<vector< Mat > >::iterator cit = A.begin(); rit != A.end(); rit++)
+    {
+      row = cat(2,row,*cit);
+    }
+    // Concatenate rows vertically
+    C = cat(1,C,row);
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

cat.h

@@ -1,5 +1,6 @@
 #ifndef IGL_CAT_H
 #define IGL_CAT_H
+#include "igl_inline.h"
 
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Sparse>
@@ -27,20 +28,20 @@ namespace igl
   //   C  output matrix
   //   
   template <typename Scalar>
-  inline void cat(
+  IGL_INLINE void cat(
       const int dim, 
       const Eigen::SparseMatrix<Scalar> & A, 
       const Eigen::SparseMatrix<Scalar> & B, 
       Eigen::SparseMatrix<Scalar> & C);
   template <typename Derived, class MatC>
-  inline void cat(
+  IGL_INLINE void cat(
     const int dim,
     const Eigen::MatrixBase<Derived> & A, 
     const Eigen::MatrixBase<Derived> & B,
     MatC & C);
   // Wrapper that returns C
   template <class Mat>
-  inline Mat cat(const int dim, const Mat & A, const Mat & B);
+  IGL_INLINE Mat cat(const int dim, const Mat & A, const Mat & B);
 
   // Note: Maybe we can autogenerate a bunch of overloads D = cat(int,A,B,C),
   // E = cat(int,A,B,C,D), etc. 
@@ -53,140 +54,11 @@ namespace igl
   // Output:
   //   C
   template <class Mat>
-  inline void cat(const std::vector<std::vector< Mat > > & A, Mat & C);
+  IGL_INLINE void cat(const std::vector<std::vector< Mat > > & A, Mat & C);
 }
 
-// Implementation
-
-// Sparse matrices need to be handled carefully. Because C++ does not 
-// Template:
-//   Scalar  sparse matrix scalar type, e.g. double
-template <typename Scalar>
-inline void igl::cat(
-    const int dim, 
-    const Eigen::SparseMatrix<Scalar> & A, 
-    const Eigen::SparseMatrix<Scalar> & B, 
-    Eigen::SparseMatrix<Scalar> & C)
-{
-  assert(dim == 1 || dim == 2);
-  using namespace Eigen;
-  // Special case if B or A is empty
-  if(A.size() == 0)
-  {
-    C = B;
-    return;
-  }
-  if(B.size() == 0)
-  {
-    C = A;
-    return;
-  }
-
-  DynamicSparseMatrix<Scalar, RowMajor> dyn_C;
-  if(dim == 1)
-  {
-    assert(A.cols() == B.cols());
-    dyn_C.resize(A.rows()+B.rows(),A.cols());
-  }else if(dim == 2)
-  {
-    assert(A.rows() == B.rows());
-    dyn_C.resize(A.rows(),A.cols()+B.cols());
-  }else
-  {
-    fprintf(stderr,"cat.h: Error: Unsupported dimension %d\n",dim);
-  }
-
-  dyn_C.reserve(A.nonZeros()+B.nonZeros());
-
-  // Iterate over outside of A
-  for(int k=0; k<A.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
-    {
-      dyn_C.coeffRef(it.row(),it.col()) += it.value();
-    }
-  }
-
-  // Iterate over outside of B
-  for(int k=0; k<B.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
-    {
-      int r = (dim == 1 ? A.rows()+it.row() : it.row());
-      int c = (dim == 2 ? A.cols()+it.col() : it.col());
-      dyn_C.coeffRef(r,c) += it.value();
-    }
-  }
-
-  C = SparseMatrix<Scalar>(dyn_C);
-}
-
-template <typename Derived, class MatC>
-inline void igl::cat(
-  const int dim,
-  const Eigen::MatrixBase<Derived> & A, 
-  const Eigen::MatrixBase<Derived> & B,
-  MatC & C)
-{
-  assert(dim == 1 || dim == 2);
-  // Special case if B or A is empty
-  if(A.size() == 0)
-  {
-    C = B;
-    return;
-  }
-  if(B.size() == 0)
-  {
-    C = A;
-    return;
-  }
-
-  if(dim == 1)
-  {
-    assert(A.cols() == B.cols());
-    C.resize(A.rows()+B.rows(),A.cols());
-    C << A,B;
-  }else if(dim == 2)
-  {
-    assert(A.rows() == B.rows());
-    C.resize(A.rows(),A.cols()+B.cols());
-    C << A,B;
-  }else
-  {
-    fprintf(stderr,"cat.h: Error: Unsupported dimension %d\n",dim);
-  }
-}
-
-template <class Mat>
-inline Mat igl::cat(const int dim, const Mat & A, const Mat & B)
-{
-  assert(dim == 1 || dim == 2);
-  Mat C;
-  igl::cat(dim,A,B,C);
-  return C;
-}
-
-template <class Mat>
-inline void cat(const std::vector<std::vector< Mat > > & A, Mat & C)
-{
-  using namespace igl;
-  using namespace std;
-  // Start with empty matrix
-  C.resize(0,0);
-  for(typename vector<vector< Mat > >::const_iterator rit = A.begin(); rit != A.end(); rit++)
-  {
-    // Concatenate each row horizontally
-    // Start with empty matrix
-    Mat row(0,0);
-    for(typename vector<vector< Mat > >::iterator cit = A.begin(); rit != A.end(); rit++)
-    {
-      row = cat(2,row,*cit);
-    }
-    // Concatenate rows vertically
-    C = cat(1,C,row);
-  }
-}
+#ifdef IGL_HEADER_ONLY
+#  include "cat.cpp"
+#endif
 
 #endif

cocoa_key_to_anttweakbar_key.cpp

@@ -0,0 +1,72 @@
+#include "cocoa_key_to_anttweakbar_key.h"
+
+#include <AntTweakBar.h>
+
+IGL_INLINE int igl::cocoa_key_to_anttweakbar_key(int key)
+{
+  // I've left commented the AntTweakBar key codes that correspond to keys I
+  // don't have on my keyboard. Please fill this in if you have those keys
+  switch(key)
+  {
+    case 127:
+      return TW_KEY_BACKSPACE;
+    case 9:
+      return TW_KEY_TAB;
+  //TW_KEY_CLEAR        = 0x0c,
+    case 3://ENTER
+    case 13:
+      return TW_KEY_RETURN;
+    case 27:
+      return TW_KEY_ESCAPE;
+    case 32:
+      return TW_KEY_SPACE;
+    case 63272:
+      return TW_KEY_DELETE;
+    case 63232:
+      return TW_KEY_UP;
+    case 63233:
+      return TW_KEY_DOWN;
+    case 63235:
+      return TW_KEY_RIGHT;
+    case 63234:
+      return TW_KEY_LEFT;
+  //TW_KEY_INSERT,
+  //TW_KEY_HOME,
+  //TW_KEY_END,
+  //TW_KEY_PAGE_UP,
+  //TW_KEY_PAGE_DOWN,
+    case 63236:
+      return TW_KEY_F1;
+    case 63237:
+      return TW_KEY_F2;
+    case 63238:
+      return TW_KEY_F3;
+    case 63239:
+      return TW_KEY_F4;
+    case 63240:
+      return TW_KEY_F5;
+    case 63241:
+      return TW_KEY_F6;
+    case 63242:
+      return TW_KEY_F7;
+    case 63243:
+      return TW_KEY_F8;
+    case 63244:
+      return TW_KEY_F9;
+    case 63245:
+      return TW_KEY_F10;
+    case 63246:
+      return TW_KEY_F11;
+    case 63247:
+      return TW_KEY_F12;
+    case 63248:
+      return TW_KEY_F13;
+    case 63249:
+      return TW_KEY_F14;
+    case 63250:
+      return TW_KEY_F15;
+    default:
+      break;
+  }
+  return key;
+}

cocoa_key_to_anttweakbar_key.h

@@ -1,5 +1,6 @@
 #ifndef IGL_COCOA_KEY_TO_ANTTWEAKBAR_KEY_H
 #define IGL_COCOA_KEY_TO_ANTTWEAKBAR_KEY_H
+#include "igl_inline.h"
 
 
 namespace igl
@@ -10,78 +11,11 @@ namespace igl
   // Inputs:
   //   key  unsigned char key from keyboard
   // Returns int of new key code 
-  inline int cocoa_key_to_anttweakbar_key(int key);
+  IGL_INLINE int cocoa_key_to_anttweakbar_key(int key);
 }
 
-// Implementation
-#include <AntTweakBar.h>
+#ifdef IGL_HEADER_ONLY
+#  include "cocoa_key_to_anttweakbar_key.cpp"
+#endif
 
-inline int igl::cocoa_key_to_anttweakbar_key(int key)
-{
-  // I've left commented the AntTweakBar key codes that correspond to keys I
-  // don't have on my keyboard. Please fill this in if you have those keys
-  switch(key)
-  {
-    case 127:
-      return TW_KEY_BACKSPACE;
-    case 9:
-      return TW_KEY_TAB;
-  //TW_KEY_CLEAR        = 0x0c,
-    case 3://ENTER
-    case 13:
-      return TW_KEY_RETURN;
-    case 27:
-      return TW_KEY_ESCAPE;
-    case 32:
-      return TW_KEY_SPACE;
-    case 63272:
-      return TW_KEY_DELETE;
-    case 63232:
-      return TW_KEY_UP;
-    case 63233:
-      return TW_KEY_DOWN;
-    case 63235:
-      return TW_KEY_RIGHT;
-    case 63234:
-      return TW_KEY_LEFT;
-  //TW_KEY_INSERT,
-  //TW_KEY_HOME,
-  //TW_KEY_END,
-  //TW_KEY_PAGE_UP,
-  //TW_KEY_PAGE_DOWN,
-    case 63236:
-      return TW_KEY_F1;
-    case 63237:
-      return TW_KEY_F2;
-    case 63238:
-      return TW_KEY_F3;
-    case 63239:
-      return TW_KEY_F4;
-    case 63240:
-      return TW_KEY_F5;
-    case 63241:
-      return TW_KEY_F6;
-    case 63242:
-      return TW_KEY_F7;
-    case 63243:
-      return TW_KEY_F8;
-    case 63244:
-      return TW_KEY_F9;
-    case 63245:
-      return TW_KEY_F10;
-    case 63246:
-      return TW_KEY_F11;
-    case 63247:
-      return TW_KEY_F12;
-    case 63248:
-      return TW_KEY_F13;
-    case 63249:
-      return TW_KEY_F14;
-    case 63250:
-      return TW_KEY_F15;
-    default:
-      break;
-  }
-  return key;
-}
 #endif

colon.cpp

@@ -0,0 +1,71 @@
+#include "colon.h"
+
+#include <cstdio>
+
+template <typename L,typename S,typename H,typename T>
+IGL_INLINE void igl::colon(
+  const L low, 
+  const S step, 
+  const H hi, 
+  Eigen::Matrix<T,Eigen::Dynamic,1> & I)
+{
+  if(low < hi)
+  {
+    if(step < 0)
+    {
+      I.resize(0);
+      fprintf(stderr,"Error: colon() low(%g)<hi(%g) but step(%g)<0\n",
+        (double)low,
+        (double)hi,
+        (double)step);
+      return;
+    }
+  }
+  if(low > hi)
+  {
+    if(step > 0)
+    {
+      I.resize(0);
+      fprintf(stderr,"Error: colon() low(%g)>hi(%g) but step(%g)<0\n",
+        (double)low,
+        (double)hi,
+        (double)step);
+      return;
+    }
+  }
+  // resize output
+  int n = floor(double((hi-low)/step))+1;
+  I.resize(n);
+  int i = 0;
+  T v = (T)low;
+  while((low<hi && (H)v<=hi) || (low>hi && (H)v>=hi))
+  {
+    I(i) = v;
+    v = v + (T)step;
+    i++;
+  }
+  assert(i==n);
+}
+
+template <typename L,typename H,typename T>
+IGL_INLINE void igl::colon(
+  const L low, 
+  const H hi, 
+  Eigen::Matrix<T,Eigen::Dynamic,1> & I)
+{
+  return igl::colon(low,(T)1,hi,I);
+}
+
+template <typename T,typename L,typename H>
+IGL_INLINE Eigen::Matrix<T,Eigen::Dynamic,1> igl::colon(
+  const L low, 
+  const H hi)
+{
+  Eigen::Matrix<T,Eigen::Dynamic,1> I;
+  igl::colon(low,hi,I);
+  return I;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

colon.h

@@ -1,5 +1,6 @@
 #ifndef IGL_COLON_H
 #define IGL_COLON_H
+#include "igl_inline.h"
 #include <Eigen/Dense>
 namespace igl
 {
@@ -26,90 +27,26 @@ namespace igl
   // Output:
   //   I  list of values from low to hi with step size step
   template <typename L,typename S,typename H,typename T>
-  inline void colon(
+  IGL_INLINE void colon(
     const L low, 
     const S step, 
     const H hi, 
     Eigen::Matrix<T,Eigen::Dynamic,1> & I);
   // Same as above but step == (T)1
   template <typename L,typename H,typename T>
-  inline void colon(
+  IGL_INLINE void colon(
     const L low, 
     const H hi, 
     Eigen::Matrix<T,Eigen::Dynamic,1> & I);
   // Return output rather than set in reference
   template <typename T,typename L,typename H>
-  inline Eigen::Matrix<T,Eigen::Dynamic,1> colon(
+  IGL_INLINE Eigen::Matrix<T,Eigen::Dynamic,1> colon(
     const L low, 
     const H hi);
 }
 
-// Implementation
-#include <cstdio>
-
-template <typename L,typename S,typename H,typename T>
-inline void igl::colon(
-  const L low, 
-  const S step, 
-  const H hi, 
-  Eigen::Matrix<T,Eigen::Dynamic,1> & I)
-{
-  if(low < hi)
-  {
-    if(step < 0)
-    {
-      I.resize(0);
-      fprintf(stderr,"Error: colon() low(%g)<hi(%g) but step(%g)<0\n",
-        (double)low,
-        (double)hi,
-        (double)step);
-      return;
-    }
-  }
-  if(low > hi)
-  {
-    if(step > 0)
-    {
-      I.resize(0);
-      fprintf(stderr,"Error: colon() low(%g)>hi(%g) but step(%g)<0\n",
-        (double)low,
-        (double)hi,
-        (double)step);
-      return;
-    }
-  }
-  // resize output
-  int n = floor(double((hi-low)/step))+1;
-  I.resize(n);
-  int i = 0;
-  T v = (T)low;
-  while((low<hi && (H)v<=hi) || (low>hi && (H)v>=hi))
-  {
-    I(i) = v;
-    v = v + (T)step;
-    i++;
-  }
-  assert(i==n);
-}
-
-template <typename L,typename H,typename T>
-inline void igl::colon(
-  const L low, 
-  const H hi, 
-  Eigen::Matrix<T,Eigen::Dynamic,1> & I)
-{
-  return igl::colon(low,(T)1,hi,I);
-}
-
-template <typename T,typename L,typename H>
-inline Eigen::Matrix<T,Eigen::Dynamic,1> igl::colon(
-  const L low, 
-  const H hi)
-{
-  Eigen::Matrix<T,Eigen::Dynamic,1> I;
-  igl::colon(low,hi,I);
-  return I;
-}
-
+#ifdef IGL_HEADER_ONLY
+#  include "colon.cpp"
 #endif
 
+#endif

concat.cpp

@@ -0,0 +1,42 @@
+#include "concat.h"
+
+#include <cstdio>
+
+template <typename T>
+IGL_INLINE void igl::concat(
+                   const T A, 
+                   const T B,
+                   const bool horiz,                 
+                   T& O)
+{
+  if (horiz)
+  {
+    // O = [A,B]
+    assert(A.rows() == B.rows());
+    O = T(A.rows(),A.cols()+B.cols());
+    O << A,B;
+  }
+  else
+  {
+    // O = [A;B]
+    assert(A.cols() == B.cols());
+    O = T(A.rows()+B.rows(),A.cols());
+    O << A,B;
+  }
+}
+
+template <typename T>
+IGL_INLINE T igl::concat(
+                const T A, 
+                const T B,
+                bool horiz
+                )
+{
+  T O = T(1,1);
+  concat(A,B,horiz,O);
+  return O;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

concat.h

@@ -1,5 +1,6 @@
 #ifndef IGL_CONCAT_H
 #define IGL_CONCAT_H
+#include "igl_inline.h"
 #include <Eigen/Dense>
 namespace igl
 {
@@ -13,14 +14,14 @@ namespace igl
   // Output:
   //   O if horiz = false return [A;B] else [A,B]
   template <typename T>
-  inline void concat(
+  IGL_INLINE void concat(
                      const T A, 
                      const T B,
                      const bool horiz,                 
                      T& O);
   
   template <typename T>
-  inline T concat(
+  IGL_INLINE T concat(
                   const T A, 
                   const T B,
                   bool horiz = false
@@ -28,43 +29,8 @@ namespace igl
 
 }
 
-// Implementation
-#include <cstdio>
-
-template <typename T>
-inline void igl::concat(
-                   const T A, 
-                   const T B,
-                   const bool horiz,                 
-                   T& O)
-{
-  if (horiz)
-  {
-    // O = [A,B]
-    assert(A.rows() == B.rows());
-    O = T(A.rows(),A.cols()+B.cols());
-    O << A,B;
-  }
-  else
-  {
-    // O = [A;B]
-    assert(A.cols() == B.cols());
-    O = T(A.rows()+B.rows(),A.cols());
-    O << A,B;
-  }
-}
-
-template <typename T>
-inline T igl::concat(
-                const T A, 
-                const T B,
-                bool horiz
-                )
-{
-  T O = T(1,1);
-  concat(A,B,horiz,O);
-  return O;
-}
-
+#ifdef IGL_HEADER_ONLY
+#  include "concat.cpp"
 #endif
 
+#endif

cotangent.cpp

@@ -0,0 +1,129 @@
+#include "cotangent.h"
+
+#include "verbose.h"
+#include <Eigen/Dense>
+
+template <class MatV, class MatF, class MatC>
+IGL_INLINE void igl::cotangent(const MatV & V, const MatF & F, MatC & C)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  // simplex size (3: triangles, 4: tetrahedra)
+  int simplex_size = F.cols();
+  // Number of elements
+  int m = F.rows();
+
+  if(simplex_size == 3)
+  {
+    // Triangles
+    // edge lengths numbered same as opposite vertices
+    Matrix<typename MatC::Scalar,Dynamic,3> l(m,3);
+    // loop over faces
+    for(int i = 0;i<m;i++)
+    {
+      l(i,0) = sqrt((V.row(F(i,1))-V.row(F(i,2))).array().pow(2).sum());
+      l(i,1) = sqrt((V.row(F(i,2))-V.row(F(i,0))).array().pow(2).sum());
+      l(i,2) = sqrt((V.row(F(i,0))-V.row(F(i,1))).array().pow(2).sum());
+    }
+    // semiperimeters
+    Matrix<typename MatC::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
+    assert(s.rows() == m);
+    // Heron's forumal for area
+    Matrix<typename MatC::Scalar,Dynamic,1> dblA(m);
+    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)));
+    }
+    // cotangents and diagonal entries for element matrices
+    // correctly divided by 4 (alec 2010)
+    C.resize(m,3);
+    for(int i = 0;i<m;i++)
+    {
+      C(i,0) = (l(i,1)*l(i,1) + l(i,2)*l(i,2) - l(i,0)*l(i,0))/dblA(i)/4.0;
+      C(i,1) = (l(i,2)*l(i,2) + l(i,0)*l(i,0) - l(i,1)*l(i,1))/dblA(i)/4.0;
+      C(i,2) = (l(i,0)*l(i,0) + l(i,1)*l(i,1) - l(i,2)*l(i,2))/dblA(i)/4.0;
+    }
+  }else if(simplex_size == 4)
+  {
+    // Tetrahedra
+    typedef Matrix<typename MatC::Scalar,3,1> Vec3;
+    typedef Matrix<typename MatC::Scalar,3,3> Mat3;
+    typedef Matrix<typename MatC::Scalar,3,4> Mat3x4;
+    typedef Matrix<typename MatC::Scalar,4,4> Mat4x4;
+
+    // preassemble right hand side
+    // COLUMN-MAJOR ORDER FOR LAPACK
+    Mat3x4 rhs;
+    rhs <<
+      1,0,0,-1,
+      0,1,0,-1,
+      0,0,1,-1;
+
+    bool diag_all_pos = true;
+    C.resize(m,6);
+
+    // loop over tetrahedra
+    for(int j = 0;j<F.rows();j++)
+    {
+      // points a,b,c,d make up the tetrahedra
+      size_t a = F(j,0);
+      size_t b = F(j,1);
+      size_t c = F(j,2);
+      size_t d = F(j,3);
+      //const std::vector<double> & pa = vertices[a];
+      //const std::vector<double> & pb = vertices[b];
+      //const std::vector<double> & pc = vertices[c];
+      //const std::vector<double> & pd = vertices[d];
+      Vec3 pa = V.row(a);
+      Vec3 pb = V.row(b);
+      Vec3 pc = V.row(c);
+      Vec3 pd = V.row(d);
+
+      // Following definition that appears in the appendix of: ``Interactive
+      // Topology-aware Surface Reconstruction,'' by Sharf, A. et al
+      // http://www.cs.bgu.ac.il/~asharf/Projects/InSuRe/Insure_siggraph_final.pdf
+
+      // compute transpose of jacobian Jj
+      Mat3 JTj;
+      JTj.row(0) = pa-pd;
+      JTj.row(1) = pb-pd;
+      JTj.row(2) = pc-pd;
+
+      // compute abs(determinant of JTj)/6 (volume of tet)
+      // determinant of transpose of A equals determinant of A
+      double volume = fabs(JTj.determinant())/6.0;
+      //printf("volume[%d] = %g\n",j+1,volume);
+
+      // solve Jj' * Ej = [-I -1], for Ej
+      // in other words solve JTj * Ej = [-I -1], for Ej
+      Mat3x4 Ej = JTj.inverse() * rhs;
+      // compute Ej'*Ej
+      Mat4x4 EjTEj = Ej.transpose() * Ej;
+
+      // Kj =  det(JTj)/6 * Ej'Ej 
+      Mat4x4 Kj = EjTEj*volume;
+      diag_all_pos &= Kj(0,0)>0 & Kj(1,1)>0 & Kj(2,2)>0 & Kj(3,3)>0;
+      C(j,0) = Kj(1,2);
+      C(j,1) = Kj(2,0);
+      C(j,2) = Kj(0,1);
+      C(j,3) = Kj(3,0);
+      C(j,4) = Kj(3,1);
+      C(j,5) = Kj(3,2);
+    }
+    if(diag_all_pos)
+    {
+      verbose("cotangent.h: Flipping sign of cotangent, so that cots are positive\n");
+      C *= -1.0;
+    }
+  }else
+  {
+    fprintf(stderr,
+      "cotangent.h: Error: Simplex size (%d) not supported\n", simplex_size);
+    assert(false);
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

cotangent.h

@@ -1,5 +1,6 @@
 #ifndef IGL_COTANGENT_H
 #define IGL_COTANGENT_H
+#include "igl_inline.h"
 namespace igl
 {
   // COTANGENT compute the cotangents of each angle in mesh (V,F)
@@ -16,131 +17,11 @@ namespace igl
   //     for triangles, columns correspond to edges 23,31,12
   //     for tets, columns correspond to edges 23,31,12,41,42,43
   template <class MatV, class MatF, class MatC>
-  inline void cotangent(const MatV & V, const MatF & F, MatC & C);
+  IGL_INLINE void cotangent(const MatV & V, const MatF & F, MatC & C);
 }
 
-// Implementation
-#include <verbose.h>
-#include <Eigen/Dense>
-
-template <class MatV, class MatF, class MatC>
-inline void igl::cotangent(const MatV & V, const MatF & F, MatC & C)
-{
-  using namespace igl;
-  using namespace std;
-  using namespace Eigen;
-  // simplex size (3: triangles, 4: tetrahedra)
-  int simplex_size = F.cols();
-  // Number of elements
-  int m = F.rows();
-
-  if(simplex_size == 3)
-  {
-    // Triangles
-    // edge lengths numbered same as opposite vertices
-    Matrix<typename MatC::Scalar,Dynamic,3> l(m,3);
-    // loop over faces
-    for(int i = 0;i<m;i++)
-    {
-      l(i,0) = sqrt((V.row(F(i,1))-V.row(F(i,2))).array().pow(2).sum());
-      l(i,1) = sqrt((V.row(F(i,2))-V.row(F(i,0))).array().pow(2).sum());
-      l(i,2) = sqrt((V.row(F(i,0))-V.row(F(i,1))).array().pow(2).sum());
-    }
-    // semiperimeters
-    Matrix<typename MatC::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
-    assert(s.rows() == m);
-    // Heron's forumal for area
-    Matrix<typename MatC::Scalar,Dynamic,1> dblA(m);
-    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)));
-    }
-    // cotangents and diagonal entries for element matrices
-    // correctly divided by 4 (alec 2010)
-    C.resize(m,3);
-    for(int i = 0;i<m;i++)
-    {
-      C(i,0) = (l(i,1)*l(i,1) + l(i,2)*l(i,2) - l(i,0)*l(i,0))/dblA(i)/4.0;
-      C(i,1) = (l(i,2)*l(i,2) + l(i,0)*l(i,0) - l(i,1)*l(i,1))/dblA(i)/4.0;
-      C(i,2) = (l(i,0)*l(i,0) + l(i,1)*l(i,1) - l(i,2)*l(i,2))/dblA(i)/4.0;
-    }
-  }else if(simplex_size == 4)
-  {
-    // Tetrahedra
-    typedef Matrix<typename MatC::Scalar,3,1> Vec3;
-    typedef Matrix<typename MatC::Scalar,3,3> Mat3;
-    typedef Matrix<typename MatC::Scalar,3,4> Mat3x4;
-    typedef Matrix<typename MatC::Scalar,4,4> Mat4x4;
-
-    // preassemble right hand side
-    // COLUMN-MAJOR ORDER FOR LAPACK
-    Mat3x4 rhs;
-    rhs <<
-      1,0,0,-1,
-      0,1,0,-1,
-      0,0,1,-1;
-
-    bool diag_all_pos = true;
-    C.resize(m,6);
-
-    // loop over tetrahedra
-    for(int j = 0;j<F.rows();j++)
-    {
-      // points a,b,c,d make up the tetrahedra
-      size_t a = F(j,0);
-      size_t b = F(j,1);
-      size_t c = F(j,2);
-      size_t d = F(j,3);
-      //const std::vector<double> & pa = vertices[a];
-      //const std::vector<double> & pb = vertices[b];
-      //const std::vector<double> & pc = vertices[c];
-      //const std::vector<double> & pd = vertices[d];
-      Vec3 pa = V.row(a);
-      Vec3 pb = V.row(b);
-      Vec3 pc = V.row(c);
-      Vec3 pd = V.row(d);
-
-      // Following definition that appears in the appendix of: ``Interactive
-      // Topology-aware Surface Reconstruction,'' by Sharf, A. et al
-      // http://www.cs.bgu.ac.il/~asharf/Projects/InSuRe/Insure_siggraph_final.pdf
-
-      // compute transpose of jacobian Jj
-      Mat3 JTj;
-      JTj.row(0) = pa-pd;
-      JTj.row(1) = pb-pd;
-      JTj.row(2) = pc-pd;
-
-      // compute abs(determinant of JTj)/6 (volume of tet)
-      // determinant of transpose of A equals determinant of A
-      double volume = fabs(JTj.determinant())/6.0;
-      //printf("volume[%d] = %g\n",j+1,volume);
-
-      // solve Jj' * Ej = [-I -1], for Ej
-      // in other words solve JTj * Ej = [-I -1], for Ej
-      Mat3x4 Ej = JTj.inverse() * rhs;
-      // compute Ej'*Ej
-      Mat4x4 EjTEj = Ej.transpose() * Ej;
+#ifdef IGL_HEADER_ONLY
+#  include "cotangent.cpp"
+#endif
 
-      // Kj =  det(JTj)/6 * Ej'Ej 
-      Mat4x4 Kj = EjTEj*volume;
-      diag_all_pos &= Kj(0,0)>0 & Kj(1,1)>0 & Kj(2,2)>0 & Kj(3,3)>0;
-      C(j,0) = Kj(1,2);
-      C(j,1) = Kj(2,0);
-      C(j,2) = Kj(0,1);
-      C(j,3) = Kj(3,0);
-      C(j,4) = Kj(3,1);
-      C(j,5) = Kj(3,2);
-    }
-    if(diag_all_pos)
-    {
-      verbose("cotangent.h: Flipping sign of cotangent, so that cots are positive\n");
-      C *= -1.0;
-    }
-  }else
-  {
-    fprintf(stderr,
-      "cotangent.h: Error: Simplex size (%d) not supported\n", simplex_size);
-    assert(false);
-  }
-}
 #endif

cotmatrix.cpp

@@ -0,0 +1,71 @@
+#include "cotmatrix.h"
+
+// For error printing
+#include <cstdio>
+#include "cotangent.h"
+
+template <typename DerivedV, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L)
+{
+  using namespace igl;
+  using namespace Eigen;
+
+  DynamicSparseMatrix<Scalar, RowMajor> dyn_L (V.rows(), V.rows());
+  Matrix<int,Dynamic,2> edges;
+  int simplex_size = F.cols();
+  // 3 for triangles, 4 for tets
+  assert(simplex_size == 3 || simplex_size == 4);
+  if(simplex_size == 3)
+  {
+    // This is important! it could decrease the comptuation time by a factor of 2
+    // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
+    // row
+    dyn_L.reserve(7*V.rows());
+    edges.resize(3,2);
+    edges << 
+      1,2,
+      2,0,
+      0,1;
+  }else if(simplex_size == 4)
+  {
+    dyn_L.reserve(17*V.rows());
+    edges.resize(6,2);
+    edges << 
+      1,2,
+      2,0,
+      0,1,
+      3,0,
+      3,1,
+      3,2;
+  }else
+  {
+    return;
+  }
+  // Gather cotangents
+  Matrix<Scalar,Dynamic,Dynamic> C;
+  cotangent(V,F,C);
+  
+  // Loop over triangles
+  for(int i = 0; i < F.rows(); i++)
+  {
+    // loop over edges of element
+    for(int e = 0;e<edges.rows();e++)
+    {
+      int source = F(i,edges(e,0));
+      int dest = F(i,edges(e,1));
+      dyn_L.coeffRef(source,dest) += C(i,e);
+      dyn_L.coeffRef(dest,source) += C(i,e);
+      dyn_L.coeffRef(source,source) += -C(i,e);
+      dyn_L.coeffRef(dest,dest) += -C(i,e);
+    }
+  }
+    // Corner indices of this triangle
+  L = SparseMatrix<Scalar>(dyn_L);
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

cotmatrix.h

@@ -1,5 +1,6 @@
 #ifndef IGL_COTMATRIX_H
 #define IGL_COTMATRIX_H
+#include "igl_inline.h"
 
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
@@ -34,77 +35,14 @@ namespace igl
   // Known bugs: off by 1e-16 on regular grid. I think its a problem of
   // arithmetic order in cotangent.h: C(i,e) = (arithmetic)/dblA/4
   template <typename DerivedV, typename DerivedF, typename Scalar>
-  inline void cotmatrix(
+  IGL_INLINE void cotmatrix(
     const Eigen::MatrixBase<DerivedV> & V, 
     const Eigen::MatrixBase<DerivedF> & F, 
     Eigen::SparseMatrix<Scalar>& L);
 }
 
-// Implementation
-
-// For error printing
-#include <cstdio>
-#include "cotangent.h"
-
-template <typename DerivedV, typename DerivedF, typename Scalar>
-inline void igl::cotmatrix(
-  const Eigen::MatrixBase<DerivedV> & V, 
-  const Eigen::MatrixBase<DerivedF> & F, 
-  Eigen::SparseMatrix<Scalar>& L)
-{
-  using namespace igl;
-  using namespace Eigen;
+#ifdef IGL_HEADER_ONLY
+#  include "cotmatrix.cpp"
+#endif
 
-  DynamicSparseMatrix<Scalar, RowMajor> dyn_L (V.rows(), V.rows());
-  Matrix<int,Dynamic,2> edges;
-  int simplex_size = F.cols();
-  // 3 for triangles, 4 for tets
-  assert(simplex_size == 3 || simplex_size == 4);
-  if(simplex_size == 3)
-  {
-    // This is important! it could decrease the comptuation time by a factor of 2
-    // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
-    // row
-    dyn_L.reserve(7*V.rows());
-    edges.resize(3,2);
-    edges << 
-      1,2,
-      2,0,
-      0,1;
-  }else if(simplex_size == 4)
-  {
-    dyn_L.reserve(17*V.rows());
-    edges.resize(6,2);
-    edges << 
-      1,2,
-      2,0,
-      0,1,
-      3,0,
-      3,1,
-      3,2;
-  }else
-  {
-    return;
-  }
-  // Gather cotangents
-  Matrix<Scalar,Dynamic,Dynamic> C;
-  cotangent(V,F,C);
-  
-  // Loop over triangles
-  for(int i = 0; i < F.rows(); i++)
-  {
-    // loop over edges of element
-    for(int e = 0;e<edges.rows();e++)
-    {
-      int source = F(i,edges(e,0));
-      int dest = F(i,edges(e,1));
-      dyn_L.coeffRef(source,dest) += C(i,e);
-      dyn_L.coeffRef(dest,source) += C(i,e);
-      dyn_L.coeffRef(source,source) += -C(i,e);
-      dyn_L.coeffRef(dest,dest) += -C(i,e);
-    }
-  }
-    // Corner indices of this triangle
-  L = SparseMatrix<Scalar>(dyn_L);
-}
 #endif

create_index_vbo.cpp

@@ -0,0 +1,39 @@
+#include "create_index_vbo.h"
+
+// http://www.songho.ca/opengl/gl_vbo.html#create
+IGL_INLINE void igl::create_index_vbo(
+  const Eigen::MatrixXi & F,
+  GLuint & F_vbo_id)
+{
+  // Generate Buffers
+  glGenBuffersARB(1,&F_vbo_id);
+  // Bind Buffers
+  glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB,F_vbo_id);
+  // Copy data to buffers
+  // We expect a matrix with each vertex position on a row, we then want to
+  // pass this data to OpenGL reading across rows (row-major)
+  if(F.Options & Eigen::RowMajor)
+  {
+    glBufferDataARB(
+      GL_ELEMENT_ARRAY_BUFFER_ARB,
+      sizeof(int)*F.size(),
+      F.data(),
+      GL_STATIC_DRAW_ARB);
+  }else
+  {
+    // Create temporary copy of transpose
+    Eigen::MatrixXi FT = F.transpose();
+    // If its column major then we need to temporarily store a transpose
+    glBufferDataARB(
+      GL_ELEMENT_ARRAY_BUFFER_ARB,
+      sizeof(int)*F.size(),
+      FT.data(),
+      GL_STATIC_DRAW);
+  }
+  // bind with 0, so, switch back to normal pointer operation
+  glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

create_index_vbo.h

@@ -1,5 +1,6 @@
 #ifndef IGL_CREATE_INDEX_VBO_H
 #define IGL_CREATE_INDEX_VBO_H
+#include "igl_inline.h"
 // NOTE: It wouldn't be so hard to template this using Eigen's templates
 
 #include <Eigen/Core>
@@ -25,46 +26,13 @@ namespace igl
   // Outputs:
   //   F_vbo_id  buffer id for face indices
   //
-  inline void create_index_vbo(
+  IGL_INLINE void create_index_vbo(
     const Eigen::MatrixXi & F,
     GLuint & F_vbo_id);
 }
 
-// Implementation
-
-// http://www.songho.ca/opengl/gl_vbo.html#create
-inline void igl::create_index_vbo(
-  const Eigen::MatrixXi & F,
-  GLuint & F_vbo_id)
-{
-  // Generate Buffers
-  glGenBuffersARB(1,&F_vbo_id);
-  // Bind Buffers
-  glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB,F_vbo_id);
-  // Copy data to buffers
-  // We expect a matrix with each vertex position on a row, we then want to
-  // pass this data to OpenGL reading across rows (row-major)
-  if(F.Options & Eigen::RowMajor)
-  {
-    glBufferDataARB(
-      GL_ELEMENT_ARRAY_BUFFER_ARB,
-      sizeof(int)*F.size(),
-      F.data(),
-      GL_STATIC_DRAW_ARB);
-  }else
-  {
-    // Create temporary copy of transpose
-    Eigen::MatrixXi FT = F.transpose();
-    // If its column major then we need to temporarily store a transpose
-    glBufferDataARB(
-      GL_ELEMENT_ARRAY_BUFFER_ARB,
-      sizeof(int)*F.size(),
-      FT.data(),
-      GL_STATIC_DRAW);
-  }
-  // bind with 0, so, switch back to normal pointer operation
-  glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
-}
-
+#ifdef IGL_HEADER_ONLY
+#  include "create_index_vbo.cpp"
 #endif
 
+#endif

create_mesh_vbo.cpp

@@ -0,0 +1,36 @@
+#include "create_mesh_vbo.h"
+
+#include "create_vector_vbo.h"
+#include "create_index_vbo.h"
+
+// http://www.songho.ca/opengl/gl_vbo.html#create
+IGL_INLINE void igl::create_mesh_vbo(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  GLuint & V_vbo_id,
+  GLuint & F_vbo_id)
+{
+  // Create VBO for vertex position vectors
+  create_vector_vbo(V,V_vbo_id);
+  // Create VBO for face index lists
+  create_index_vbo(F,F_vbo_id);
+}
+
+// http://www.songho.ca/opengl/gl_vbo.html#create
+IGL_INLINE void igl::create_mesh_vbo(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & N,
+  GLuint & V_vbo_id,
+  GLuint & F_vbo_id,
+  GLuint & N_vbo_id)
+{
+  // Create VBOs for faces and vertices
+  create_mesh_vbo(V,F,V_vbo_id,F_vbo_id);
+  // Create VBO for normal vectors
+  create_vector_vbo(N,N_vbo_id);
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

create_mesh_vbo.h

@@ -1,5 +1,6 @@
 #ifndef IGL_CREATE_MESH_VBO_H
 #define IGL_CREATE_MESH_VBO_H
+#include "igl_inline.h"
 // NOTE: It wouldn't be so hard to template this using Eigen's templates
 
 #include <Eigen/Core>
@@ -31,7 +32,7 @@ namespace igl
   // NOTE: when using glDrawElements VBOs for V and F using MatrixXd and
   // MatrixXi will have types GL_DOUBLE and GL_UNSIGNED_INT respectively
   //
-  inline void create_mesh_vbo(
+  IGL_INLINE void create_mesh_vbo(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     GLuint & V_vbo_id,
@@ -45,7 +46,7 @@ namespace igl
   //   V_vbo_id  buffer id for vertex positions
   //   F_vbo_id  buffer id for face indices
   //   N_vbo_id  buffer id for vertex positions
-  inline void create_mesh_vbo(
+  IGL_INLINE void create_mesh_vbo(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     const Eigen::MatrixXd & N,
@@ -55,36 +56,8 @@ namespace igl
 
 }
 
-// Implementation
-#include "create_vector_vbo.h"
-#include "create_index_vbo.h"
-
-// http://www.songho.ca/opengl/gl_vbo.html#create
-inline void igl::create_mesh_vbo(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  GLuint & V_vbo_id,
-  GLuint & F_vbo_id)
-{
-  // Create VBO for vertex position vectors
-  create_vector_vbo(V,V_vbo_id);
-  // Create VBO for face index lists
-  create_index_vbo(F,F_vbo_id);
-}
-
-// http://www.songho.ca/opengl/gl_vbo.html#create
-inline void igl::create_mesh_vbo(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  const Eigen::MatrixXd & N,
-  GLuint & V_vbo_id,
-  GLuint & F_vbo_id,
-  GLuint & N_vbo_id)
-{
-  // Create VBOs for faces and vertices
-  create_mesh_vbo(V,F,V_vbo_id,F_vbo_id);
-  // Create VBO for normal vectors
-  create_vector_vbo(N,N_vbo_id);
-}
+#ifdef IGL_HEADER_ONLY
+#  include "create_mesh_vbo.cpp"
+#endif
 
 #endif

create_shader_program.cpp

@@ -0,0 +1,72 @@
+#include "create_shader_program.h"
+
+#include "load_shader.h"
+#include "print_program_info_log.h"
+#include <cstdio>
+
+IGL_INLINE bool igl::create_shader_program(
+  const std::string vert_source,
+  const std::string frag_source,
+  const std::map<std::string,GLuint> attrib,
+  GLuint & id)
+{
+  if(vert_source == "" && frag_source == "")
+  {
+    fprintf(
+      stderr,
+      "Error: create_shader_program() could not create shader program,"
+      " both .vert and .frag source given were empty\n");
+    return false;
+  }
+
+  // create program
+  id = glCreateProgram();
+  if(id == 0)
+  {
+    fprintf(
+      stderr,
+      "Error: create_shader_program() could not create shader program.\n");
+    return false;
+  }
+
+  if(vert_source != "")
+  {
+    // load vertex shader
+    GLuint v = igl::load_shader(vert_source.c_str(),GL_VERTEX_SHADER);
+    if(v == 0)
+    {
+      return false;
+    }
+    glAttachShader(id,v);
+  }
+
+  if(frag_source != "")
+  {
+    // load fragment shader
+    GLuint f = igl::load_shader(frag_source.c_str(),GL_FRAGMENT_SHADER);
+    if(f == 0)
+    {
+      return false;
+    }
+    glAttachShader(id,f);
+  }
+
+  // loop over attributes
+  for(
+    std::map<std::string,GLuint>::const_iterator ait = attrib.begin();
+    ait != attrib.end();
+    ait++)
+  {
+    glBindAttribLocation(
+      id,
+      (*ait).second,
+      (*ait).first.c_str());
+  }
+  // Link program
+  glLinkProgram(id);
+
+  // print log if any
+  igl::print_program_info_log(id);
+
+  return true;
+}

create_shader_program.h

@@ -1,5 +1,6 @@
 #ifndef IGL_CREATE_SHADER_PROGRAM_H
 #define IGL_CREATE_SHADER_PROGRAM_H
+#include "igl_inline.h"
 #include <string>
 #include <map>
 
@@ -32,82 +33,15 @@ namespace igl
   // leaking a shader (since it will be overwritten)
   //
   // See also: destroy_shader_program
-  inline bool create_shader_program(
+  IGL_INLINE bool create_shader_program(
     const std::string vert_source,
     const std::string frag_source,
     const std::map<std::string,GLuint> attrib,
     GLuint & id);
 }
 
-// Implementation
-#include "load_shader.h"
-#include "print_program_info_log.h"
-#include <cstdio>
-
-inline bool igl::create_shader_program(
-  const std::string vert_source,
-  const std::string frag_source,
-  const std::map<std::string,GLuint> attrib,
-  GLuint & id)
-{
-  if(vert_source == "" && frag_source == "")
-  {
-    fprintf(
-      stderr,
-      "Error: create_shader_program() could not create shader program,"
-      " both .vert and .frag source given were empty\n");
-    return false;
-  }
-
-  // create program
-  id = glCreateProgram();
-  if(id == 0)
-  {
-    fprintf(
-      stderr,
-      "Error: create_shader_program() could not create shader program.\n");
-    return false;
-  }
-
-  if(vert_source != "")
-  {
-    // load vertex shader
-    GLuint v = igl::load_shader(vert_source.c_str(),GL_VERTEX_SHADER);
-    if(v == 0)
-    {
-      return false;
-    }
-    glAttachShader(id,v);
-  }
-
-  if(frag_source != "")
-  {
-    // load fragment shader
-    GLuint f = igl::load_shader(frag_source.c_str(),GL_FRAGMENT_SHADER);
-    if(f == 0)
-    {
-      return false;
-    }
-    glAttachShader(id,f);
-  }
-
-  // loop over attributes
-  for(
-    std::map<std::string,GLuint>::const_iterator ait = attrib.begin();
-    ait != attrib.end();
-    ait++)
-  {
-    glBindAttribLocation(
-      id,
-      (*ait).second,
-      (*ait).first.c_str());
-  }
-  // Link program
-  glLinkProgram(id);
-
-  // print log if any
-  igl::print_program_info_log(id);
+#ifdef IGL_HEADER_ONLY
+#  include "create_shader_program.cpp"
+#endif
 
-  return true;
-}
 #endif

create_vector_vbo.cpp

@@ -0,0 +1,45 @@
+#include "create_vector_vbo.h"
+
+#include <cassert>
+
+// http://www.songho.ca/opengl/gl_vbo.html#create
+template <typename T>
+IGL_INLINE void igl::create_vector_vbo(
+  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & V,
+  GLuint & V_vbo_id)
+{
+  //// Expects that input is list of 3D vectors along rows
+  //assert(V.cols() == 3);
+
+  // Generate Buffers
+  glGenBuffers(1,&V_vbo_id);
+  // Bind Buffers
+  glBindBuffer(GL_ARRAY_BUFFER,V_vbo_id);
+  // Copy data to buffers
+  // We expect a matrix with each vertex position on a row, we then want to
+  // pass this data to OpenGL reading across rows (row-major)
+  if(V.Options & Eigen::RowMajor)
+  {
+    glBufferData(
+      GL_ARRAY_BUFFER,
+      sizeof(T)*V.size(),
+      V.data(),
+      GL_STATIC_DRAW);
+  }else
+  {
+    // Create temporary copy of transpose
+    Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> VT = V.transpose();
+    // If its column major then we need to temporarily store a transpose
+    glBufferData(
+      GL_ARRAY_BUFFER,
+      sizeof(T)*V.size(),
+      VT.data(),
+      GL_STATIC_DRAW);
+  }
+  // bind with 0, so, switch back to normal pointer operation
+  glBindBuffer(GL_ARRAY_BUFFER, 0);
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

create_vector_vbo.h

@@ -1,5 +1,6 @@
 #ifndef IGL_CREATE_VECTOR_VBO_H
 #define IGL_CREATE_VECTOR_VBO_H
+#include "igl_inline.h"
 // NOTE: It wouldn't be so hard to template this using Eigen's templates
 
 #include <Eigen/Core>
@@ -28,51 +29,13 @@ namespace igl
   //   V_vbo_id  buffer id for vectors
   //
   template <typename T>
-  inline void create_vector_vbo(
+  IGL_INLINE void create_vector_vbo(
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & V,
     GLuint & V_vbo_id);
 }
 
-// Implementation
-#include <cassert>
-
-// http://www.songho.ca/opengl/gl_vbo.html#create
-template <typename T>
-inline void igl::create_vector_vbo(
-  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & V,
-  GLuint & V_vbo_id)
-{
-  //// Expects that input is list of 3D vectors along rows
-  //assert(V.cols() == 3);
-
-  // Generate Buffers
-  glGenBuffers(1,&V_vbo_id);
-  // Bind Buffers
-  glBindBuffer(GL_ARRAY_BUFFER,V_vbo_id);
-  // Copy data to buffers
-  // We expect a matrix with each vertex position on a row, we then want to
-  // pass this data to OpenGL reading across rows (row-major)
-  if(V.Options & Eigen::RowMajor)
-  {
-    glBufferData(
-      GL_ARRAY_BUFFER,
-      sizeof(T)*V.size(),
-      V.data(),
-      GL_STATIC_DRAW);
-  }else
-  {
-    // Create temporary copy of transpose
-    Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> VT = V.transpose();
-    // If its column major then we need to temporarily store a transpose
-    glBufferData(
-      GL_ARRAY_BUFFER,
-      sizeof(T)*V.size(),
-      VT.data(),
-      GL_STATIC_DRAW);
-  }
-  // bind with 0, so, switch back to normal pointer operation
-  glBindBuffer(GL_ARRAY_BUFFER, 0);
-}
-
+#ifdef IGL_HEADER_ONLY
+#  include "create_vector_vbo.cpp"
 #endif
 
+#endif

cross.cpp

@@ -0,0 +1,12 @@
+#include "cross.h"
+
+// http://www.antisphere.com/Wiki/tools:anttweakbar
+IGL_INLINE void igl::cross(
+  const double *a, 
+  const double *b,
+  double *out)
+{
+  out[0] = a[1]*b[2]-a[2]*b[1];
+  out[1] = a[2]*b[0]-a[0]*b[2];
+  out[2] = a[0]*b[1]-a[1]*b[0];
+}

cross.h

@@ -1,5 +1,6 @@
 #ifndef IGL_CROSS_H
 #define IGL_CROSS_H
+#include "igl_inline.h"
 namespace igl
 {
   // Computes out = cross(a,b)
@@ -8,21 +9,14 @@ namespace igl
   //   b  right 3d vector
   // Outputs:
   //   out  result 3d vector
-  inline void cross(
+  IGL_INLINE void cross(
     const double *a, 
     const double *b,
     double *out);
 }
 
-// Implementation
-// http://www.antisphere.com/Wiki/tools:anttweakbar
-inline void igl::cross(
-  const double *a, 
-  const double *b,
-  double *out)
-{
-  out[0] = a[1]*b[2]-a[2]*b[1];
-  out[1] = a[2]*b[0]-a[0]*b[2];
-  out[2] = a[0]*b[1]-a[1]*b[0];
-}
+#ifdef IGL_HEADER_ONLY
+#  include "cross.cpp"
+#endif
+
 #endif

destroy_shader_program.cpp

@@ -0,0 +1,32 @@
+#include "destroy_shader_program.h"
+#include <cstdio>
+
+IGL_INLINE bool igl::destroy_shader_program(const GLuint id)
+{
+  // Don't try to destroy id == 0 (no shader program)
+  if(id == 0)
+  {
+    fprintf(stderr,"Error: destroy_shader_program() id = %d"
+      " but must should be positive\n",id);
+    return false;
+  }
+  // Get each attached shader one by one and detach and delete it
+  GLsizei count;
+  // shader id
+  GLuint s;
+  do
+  {
+    // Try to get at most *1* attached shader
+    glGetAttachedShaders(id,1,&count,&s);
+    // Check that we actually got *1*
+    if(count == 1)
+    {
+      // Detach and delete this shader
+      glDetachShader(id,s);
+      glDeleteShader(s);
+    }
+  }while(count > 0);
+  // Now that all of the shaders are gone we can just delete the program
+  glDeleteProgram(id);
+  return true;
+}

destroy_shader_program.h

@@ -1,5 +1,6 @@
 #ifndef IGL_DESTROY_SHADER_PROGRAM_H
 #define IGL_DESTROY_SHADER_PROGRAM_H
+#include "igl_inline.h"
 
 #ifdef __APPLE__
 #  include <OpenGL/gl.h>
@@ -24,38 +25,11 @@ namespace igl
   // to use id as if it still contains a program
   // 
   // See also: create_shader_program
-  inline bool destroy_shader_program(const GLuint id);
+  IGL_INLINE bool destroy_shader_program(const GLuint id);
 }
 
-// Implementation
-inline bool igl::destroy_shader_program(const GLuint id)
-{
-  // Don't try to destroy id == 0 (no shader program)
-  if(id == 0)
-  {
-    fprintf(stderr,"Error: destroy_shader_program() id = %d"
-      " but must should be positive\n",id);
-    return false;
-  }
-  // Get each attached shader one by one and detach and delete it
-  GLsizei count;
-  // shader id
-  GLuint s;
-  do
-  {
-    // Try to get at most *1* attached shader
-    glGetAttachedShaders(id,1,&count,&s);
-    // Check that we actually got *1*
-    if(count == 1)
-    {
-      // Detach and delete this shader
-      glDetachShader(id,s);
-      glDeleteShader(s);
-    }
-  }while(count > 0);
-  // Now that all of the shaders are gone we can just delete the program
-  glDeleteProgram(id);
-  return true;
-}
+#ifdef IGL_HEADER_ONLY
+#  include "destroy_shader_program.cpp"
+#endif
 
 #endif

diag.cpp

@@ -0,0 +1,48 @@
+#include "diag.h"
+
+#include "verbose.h"
+
+template <typename T>
+IGL_INLINE void igl::diag(
+  const Eigen::SparseMatrix<T>& X, 
+  Eigen::SparseVector<T>& V)
+{
+  // Get size of input
+  int m = X.rows();
+  int n = X.cols();
+  V = Eigen::SparseVector<T>((m>n?n:m));
+
+  // Iterate over outside
+  for(int k=0; k<X.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
+    {
+      if(it.col() == it.row())
+      {
+        V.coeffRef(it.col()) += it.value();
+      }
+    }
+  }
+}
+
+template <typename T>
+IGL_INLINE void igl::diag(
+  const Eigen::SparseVector<T>& V,
+  Eigen::SparseMatrix<T>& X)
+{
+  // clear and resize output
+  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(V.size(),V.size());
+
+  // loop over non-zeros
+  for(typename Eigen::SparseVector<T>::InnerIterator it(V); it; ++it)
+  {
+    dyn_X.coeffRef(it.index(),it.index()) += it.value();
+  }
+
+  X = Eigen::SparseMatrix<T>(dyn_X);
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

diag.h

@@ -1,5 +1,7 @@
 #ifndef IGL_DIAG_H
 #define IGL_DIAG_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Sparse>
 
 namespace igl
@@ -15,7 +17,7 @@ namespace igl
   // Outputs:
   //   V  a min(m,n) sparse vector
   template <typename T>
-  inline void diag(
+  IGL_INLINE void diag(
     const Eigen::SparseMatrix<T>& X, 
     Eigen::SparseVector<T>& V);
   // Templates:
@@ -25,52 +27,13 @@ namespace igl
   // Outputs:
   //  X  a m by m sparse matrix
   template <typename T>
-  inline void diag(
+  IGL_INLINE void diag(
     const Eigen::SparseVector<T>& V,
     Eigen::SparseMatrix<T>& X);
 }
 
-// Implementation
-#include "verbose.h"
-
-template <typename T>
-inline void igl::diag(
-  const Eigen::SparseMatrix<T>& X, 
-  Eigen::SparseVector<T>& V)
-{
-  // Get size of input
-  int m = X.rows();
-  int n = X.cols();
-  V = Eigen::SparseVector<T>((m>n?n:m));
-
-  // Iterate over outside
-  for(int k=0; k<X.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
-    {
-      if(it.col() == it.row())
-      {
-        V.coeffRef(it.col()) += it.value();
-      }
-    }
-  }
-}
-
-template <typename T>
-inline void igl::diag(
-  const Eigen::SparseVector<T>& V,
-  Eigen::SparseMatrix<T>& X)
-{
-  // clear and resize output
-  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(V.size(),V.size());
-
-  // loop over non-zeros
-  for(typename Eigen::SparseVector<T>::InnerIterator it(V); it; ++it)
-  {
-    dyn_X.coeffRef(it.index(),it.index()) += it.value();
-  }
+#ifdef IGL_HEADER_ONLY
+#  include "diag.cpp"
+#endif
 
-  X = Eigen::SparseMatrix<T>(dyn_X);
-}
 #endif

dirname.cpp

@@ -0,0 +1,32 @@
+#include "dirname.h"
+
+#include <algorithm>
+#include "verbose.h"
+
+IGL_INLINE std::string igl::dirname(const std::string & path)
+{
+  if(path == "")
+  {
+    return std::string("");
+  }
+  // http://stackoverflow.com/questions/5077693/dirnamephp-similar-function-in-c
+  std::string::const_reverse_iterator last_slash =
+    std::find(
+      path.rbegin(), 
+      path.rend(), '/');
+  if( last_slash == path.rend() )
+  {
+    // No slashes found
+    return std::string(".");
+  }else if(1 == (last_slash.base() - path.begin()))
+  {
+    // Slash is first char
+    return std::string("/");
+  }else if(path.end() == last_slash.base() )
+  {
+    // Slash is last char
+    std::string redo = std::string(path.begin(),path.end()-1);
+    return igl::dirname(redo);
+  }
+  return std::string(path.begin(),last_slash.base()-1);
+}

dirname.h

@@ -1,5 +1,6 @@
 #ifndef IGL_DIRNAME_H
 #define IGL_DIRNAME_H
+#include "igl_inline.h"
 
 #include <string>
 
@@ -11,39 +12,11 @@ namespace igl
   // Returns string containing dirname (see php's dirname)
   //
   // See also: basename, pathinfo
-  inline std::string dirname(const std::string & path);
+  IGL_INLINE std::string dirname(const std::string & path);
 }
 
-// Implementation
-#include <algorithm>
-#include "verbose.h"
-
-inline std::string igl::dirname(const std::string & path)
-{
-  if(path == "")
-  {
-    return std::string("");
-  }
-  // http://stackoverflow.com/questions/5077693/dirnamephp-similar-function-in-c
-  std::string::const_reverse_iterator last_slash =
-    std::find(
-      path.rbegin(), 
-      path.rend(), '/');
-  if( last_slash == path.rend() )
-  {
-    // No slashes found
-    return std::string(".");
-  }else if(1 == (last_slash.base() - path.begin()))
-  {
-    // Slash is first char
-    return std::string("/");
-  }else if(path.end() == last_slash.base() )
-  {
-    // Slash is last char
-    std::string redo = std::string(path.begin(),path.end()-1);
-    return igl::dirname(redo);
-  }
-  return std::string(path.begin(),last_slash.base()-1);
-}
+#ifdef IGL_HEADER_ONLY
+#  include "dirname.cpp"
+#endif
 
 #endif

dot.cpp

@@ -0,0 +1,9 @@
+#include "dot.h"
+
+// http://www.antisphere.com/Wiki/tools:anttweakbar
+IGL_INLINE double igl::dot(
+  const double *a, 
+  const double *b)
+{
+  return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
+}

dot.h

@@ -1,5 +1,6 @@
 #ifndef IGL_DOT_H
 #define IGL_DOT_H
+#include "igl_inline.h"
 namespace igl
 {
   // Computes out = dot(a,b)
@@ -7,17 +8,13 @@ namespace igl
   //   a  left 3d vector
   //   b  right 3d vector
   // Returns scalar dot product
-  inline double dot(
+  IGL_INLINE double dot(
     const double *a, 
     const double *b);
 }
 
-// Implementation
-// http://www.antisphere.com/Wiki/tools:anttweakbar
-inline double igl::dot(
-  const double *a, 
-  const double *b)
-{
-  return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
-}
+#ifdef IGL_HEADER_ONLY
+#  include "dot.cpp"
+#endif
+
 #endif

edges.cpp

@@ -0,0 +1,31 @@
+#include "edges.h"
+
+#include <map>
+#include "adjacency_matrix.h"
+
+IGL_INLINE void igl::edges( const Eigen::MatrixXi& F, Eigen::MatrixXi& E)
+{
+  // build adjacency matrix
+  Eigen::SparseMatrix<int> A;
+  igl::adjacency_matrix(F,A);
+  // Number of non zeros should be twice number of edges
+  assert(A.nonZeros()%2 == 0);
+  // Resize to fit edges
+  E.resize(A.nonZeros()/2,2);
+  int i = 0;
+  // Iterate over outside
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(Eigen::SparseMatrix<int>::InnerIterator it (A,k); it; ++it)
+    {
+      // only add edge in one direction
+      if(it.row()<it.col())
+      {
+        E(i,0) = it.row();
+        E(i,1) = it.col();
+        i++;
+      }
+    }
+  }
+}

edges.h

@@ -1,5 +1,6 @@
 #ifndef IGL_EDGES_H
 #define IGL_EDGES_H
+#include "igl_inline.h"
 
 #include <Eigen/Dense>
 
@@ -16,38 +17,11 @@ namespace igl
   //   E #E by 2 list of edges in no particular order
   //
   // See also: adjacency_matrix
-  inline void edges( const Eigen::MatrixXi& F, Eigen::MatrixXi& E);
+  IGL_INLINE void edges( const Eigen::MatrixXi& F, Eigen::MatrixXi& E);
 }
 
-// Implementation
-#include <map>
-#include <adjacency_matrix.h>
-
-inline void igl::edges( const Eigen::MatrixXi& F, Eigen::MatrixXi& E)
-{
-  // build adjacency matrix
-  Eigen::SparseMatrix<int> A;
-  igl::adjacency_matrix(F,A);
-  // Number of non zeros should be twice number of edges
-  assert(A.nonZeros()%2 == 0);
-  // Resize to fit edges
-  E.resize(A.nonZeros()/2,2);
-  int i = 0;
-  // Iterate over outside
-  for(int k=0; k<A.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(Eigen::SparseMatrix<int>::InnerIterator it (A,k); it; ++it)
-    {
-      // only add edge in one direction
-      if(it.row()<it.col())
-      {
-        E(i,0) = it.row();
-        E(i,1) = it.col();
-        i++;
-      }
-    }
-  }
-}
+#ifdef IGL_HEADER_ONLY
+#  include "edges.cpp"
+#endif
 
 #endif

edgetopology.h

@@ -27,7 +27,7 @@ namespace igl
     const Eigen::MatrixXi& EF);
 }
 
-// Implementation
+// Broken Implementation
 #include <algorithm>
 #include "is_manifold.h"
 

faces_first.cpp

@@ -0,0 +1,79 @@
+#include "faces_first.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+template <typename MatV, typename MatF, typename VecI>
+IGL_INLINE void igl::faces_first(
+  const MatV & V, 
+  const MatF & F, 
+  MatV & RV, 
+  MatF & RF, 
+  VecI & IM)
+{
+  using namespace std;
+  using namespace Eigen;
+  vector<bool> in_face(V.rows());
+  for(int i = 0; i<F.rows(); i++)
+  {
+    for(int j = 0; j<F.cols(); j++)
+    {
+      in_face[F(i,j)] = true;
+    }
+  }
+  // count number of vertices not in faces
+  int num_in_F = 0;
+  for(int i = 0;i<V.rows();i++)
+  {
+    num_in_F += (in_face[i]?1:0);
+  }
+  // list of unique vertices that occur in F
+  VectorXi U(num_in_F);
+  // list of unique vertices that do not occur in F
+  VectorXi NU(V.rows()-num_in_F);
+  int Ui = 0;
+  int NUi = 0;
+  // loop over vertices
+  for(int i = 0;i<V.rows();i++)
+  {
+    if(in_face[i])
+    {
+      U(Ui) = i;
+      Ui++;
+    }else
+    {
+      NU(NUi) = i;
+      NUi++;
+    }
+  }
+  IM.resize(V.rows());
+  // reindex vertices that occur in faces to be first
+  for(int i = 0;i<U.size();i++)
+  {
+    IM(U(i)) = i;
+  }
+  // reindex vertices that do not occur in faces to come after those that do
+  for(int i = 0;i<NU.size();i++)
+  {
+    IM(NU(i)) = i+U.size();
+  }
+  RF.resize(F.rows(),F.cols());
+  // Reindex faces
+  for(int i = 0; i<F.rows(); i++)
+  {
+    for(int j = 0; j<F.cols(); j++)
+    {
+      RF(i,j) = IM(F(i,j));
+    }
+  }
+  RV.resize(V.rows(),V.cols());
+  // Reorder vertices
+  for(int i = 0;i<V.rows();i++)
+  {
+    RV.row(IM(i)) = V.row(i);
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

faces_first.h

@@ -1,5 +1,6 @@
 #ifndef IGL_FACES_FIRST_H
 #define IGL_FACES_FIRST_H
+#include "igl_inline.h"
 namespace igl
 {
   // FACES_FIRST Reorder vertices so that vertices in face list come before
@@ -24,7 +25,7 @@ namespace igl
   //    and RV(IM,:) = V
   //
   template <typename MatV, typename MatF, typename VecI>
-  inline void faces_first(
+  IGL_INLINE void faces_first(
     const MatV & V, 
     const MatF & F, 
     MatV & RV, 
@@ -32,79 +33,8 @@ namespace igl
     VecI & IM);
 }
 
-// Implementation
-#include <vector>
-#include <Eigen/Dense>
-
-template <typename MatV, typename MatF, typename VecI>
-inline void igl::faces_first(
-  const MatV & V, 
-  const MatF & F, 
-  MatV & RV, 
-  MatF & RF, 
-  VecI & IM)
-{
-  using namespace std;
-  using namespace Eigen;
-  vector<bool> in_face(V.rows());
-  for(int i = 0; i<F.rows(); i++)
-  {
-    for(int j = 0; j<F.cols(); j++)
-    {
-      in_face[F(i,j)] = true;
-    }
-  }
-  // count number of vertices not in faces
-  int num_in_F = 0;
-  for(int i = 0;i<V.rows();i++)
-  {
-    num_in_F += (in_face[i]?1:0);
-  }
-  // list of unique vertices that occur in F
-  VectorXi U(num_in_F);
-  // list of unique vertices that do not occur in F
-  VectorXi NU(V.rows()-num_in_F);
-  int Ui = 0;
-  int NUi = 0;
-  // loop over vertices
-  for(int i = 0;i<V.rows();i++)
-  {
-    if(in_face[i])
-    {
-      U(Ui) = i;
-      Ui++;
-    }else
-    {
-      NU(NUi) = i;
-      NUi++;
-    }
-  }
-  IM.resize(V.rows());
-  // reindex vertices that occur in faces to be first
-  for(int i = 0;i<U.size();i++)
-  {
-    IM(U(i)) = i;
-  }
-  // reindex vertices that do not occur in faces to come after those that do
-  for(int i = 0;i<NU.size();i++)
-  {
-    IM(NU(i)) = i+U.size();
-  }
-  RF.resize(F.rows(),F.cols());
-  // Reindex faces
-  for(int i = 0; i<F.rows(); i++)
-  {
-    for(int j = 0; j<F.cols(); j++)
-    {
-      RF(i,j) = IM(F(i,j));
-    }
-  }
-  RV.resize(V.rows(),V.cols());
-  // Reorder vertices
-  for(int i = 0;i<V.rows();i++)
-  {
-    RV.row(IM(i)) = V.row(i);
-  }
-}
+#ifdef IGL_HEADER_ONLY
+#  include "faces_first.cpp"
+#endif
 
 #endif

file_contents_as_string.cpp

@@ -0,0 +1,25 @@
+#include "file_contents_as_string.h"
+
+#include <fstream>
+#include <cstdio>
+
+IGL_INLINE bool igl::file_contents_as_string(
+  const std::string file_name,
+  std::string & content)
+{
+  std::ifstream ifs(file_name.c_str());
+  // Check that opening the stream worked successfully
+  if(!ifs.good())
+  {
+    fprintf(
+      stderr,
+      "IOError: file_contents_as_string() cannot open %s\n",
+      file_name.c_str());
+    return false;
+  }
+  // Stream file contents into string
+  content = std::string(
+    (std::istreambuf_iterator<char>(ifs)),
+    (std::istreambuf_iterator<char>()));
+  return true;
+}

file_contents_as_string.h

@@ -1,5 +1,6 @@
 #ifndef IGL_FILE_CONTENTS_AS_STRING_H
 #define IGL_FILE_CONTENTS_AS_STRING_H
+#include "igl_inline.h"
 
 #include <string>
 namespace igl
@@ -10,33 +11,13 @@ namespace igl
   // Outputs:
   //   content  output string containing contents of the given file
   // Returns true on succes, false on error
-  inline bool file_contents_as_string(
+  IGL_INLINE bool file_contents_as_string(
     const std::string file_name,
     std::string & content);
 }
 
-// Implementation
-#include <fstream>
-#include <cstdio>
+#ifdef IGL_HEADER_ONLY
+#  include "file_contents_as_string.cpp"
+#endif
 
-inline bool igl::file_contents_as_string(
-  const std::string file_name,
-  std::string & content)
-{
-  std::ifstream ifs(file_name.c_str());
-  // Check that opening the stream worked successfully
-  if(!ifs.good())
-  {
-    fprintf(
-      stderr,
-      "IOError: file_contents_as_string() cannot open %s\n",
-      file_name.c_str());
-    return false;
-  }
-  // Stream file contents into string
-  content = std::string(
-    (std::istreambuf_iterator<char>(ifs)),
-    (std::istreambuf_iterator<char>()));
-  return true;
-}
 #endif

file_exists.cpp

@@ -0,0 +1,9 @@
+#include "file_exists.h"
+
+#include <sys/stat.h>
+
+IGL_INLINE bool igl::file_exists(const char* filename)
+{
+  struct stat status;
+  return (stat(filename,&status)==0);
+}

file_exists.h

@@ -1,5 +1,6 @@
 #ifndef IGL_FILE_EXISTS_H
 #define IGL_FILE_EXISTS_H
+#include "igl_inline.h"
 namespace igl
 {
   // Check if a file or directory exists like PHP's file_exists function:
@@ -8,16 +9,11 @@ namespace igl
   //   filename  path to file
   // Returns true if file exists and is readable and false if file doesn't
   // exist or *is not readable*
-  inline bool file_exists(const char * filename);
+  IGL_INLINE bool file_exists(const char * filename);
 }
 
+#ifdef IGL_HEADER_ONLY
+#  include "file_exists.cpp"
+#endif
 
-// Implementation
-#include <sys/stat.h>
-
-inline bool igl::file_exists(const char* filename)
-{
-  struct stat status;
-  return (stat(filename,&status)==0);
-}
 #endif

find.cpp

@@ -0,0 +1,54 @@
+#include "find.h"
+
+#include "verbose.h"
+  
+template <typename T>
+IGL_INLINE void igl::find(
+  const Eigen::SparseMatrix<T>& X,
+  Eigen::Matrix<int,Eigen::Dynamic,1> & I,
+  Eigen::Matrix<int,Eigen::Dynamic,1> & J,
+  Eigen::Matrix<T,Eigen::Dynamic,1> & V)
+{
+  // Resize outputs to fit nonzeros
+  I.resize(X.nonZeros());
+  J.resize(X.nonZeros());
+  V.resize(X.nonZeros());
+
+  int i = 0;
+  // Iterate over outside
+  for(int k=0; k<X.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
+    {
+      V(i) = it.value();
+      I(i) = it.row();
+      J(i) = it.col();
+      i++;
+    }
+  }
+}
+  
+template <typename T>
+IGL_INLINE void igl::find(
+  const Eigen::SparseVector<T>& X,
+  Eigen::Matrix<int,Eigen::Dynamic,1> & I,
+  Eigen::Matrix<T,Eigen::Dynamic,1> & V)
+{
+  // Resize outputs to fit nonzeros
+  I.resize(X.nonZeros());
+  V.resize(X.nonZeros());
+
+  int i = 0;
+  // loop over non-zeros
+  for(typename Eigen::SparseVector<T>::InnerIterator it(X); it; ++it)
+  {
+    I(i) = it.index();
+    V(i) = it.value();
+    i++;
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

find.h

@@ -1,5 +1,6 @@
 #ifndef IGL_FIND_H
 #define IGL_FIND_H
+#include "igl_inline.h"
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
@@ -18,7 +19,7 @@ namespace igl
   //   V  nnz vector of type T non-zeros entries in X
   //
   template <typename T>
-  inline void find(
+  IGL_INLINE void find(
     const Eigen::SparseMatrix<T>& X,
     Eigen::Matrix<int,Eigen::Dynamic,1> & I,
     Eigen::Matrix<int,Eigen::Dynamic,1> & J,
@@ -36,60 +37,14 @@ namespace igl
   //   I  nnz vector of indices of non zeros entries in X
   //   V  nnz vector of type T non-zeros entries in X
   template <typename T>
-  inline void find(
+  IGL_INLINE void find(
     const Eigen::SparseVector<T>& X,
     Eigen::Matrix<int,Eigen::Dynamic,1> & I,
     Eigen::Matrix<T,Eigen::Dynamic,1> & V);
 }
 
-// Implementation
-#include "verbose.h"
-  
-template <typename T>
-inline void igl::find(
-  const Eigen::SparseMatrix<T>& X,
-  Eigen::Matrix<int,Eigen::Dynamic,1> & I,
-  Eigen::Matrix<int,Eigen::Dynamic,1> & J,
-  Eigen::Matrix<T,Eigen::Dynamic,1> & V)
-{
-  // Resize outputs to fit nonzeros
-  I.resize(X.nonZeros());
-  J.resize(X.nonZeros());
-  V.resize(X.nonZeros());
-
-  int i = 0;
-  // Iterate over outside
-  for(int k=0; k<X.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
-    {
-      V(i) = it.value();
-      I(i) = it.row();
-      J(i) = it.col();
-      i++;
-    }
-  }
-}
-  
-template <typename T>
-inline void igl::find(
-  const Eigen::SparseVector<T>& X,
-  Eigen::Matrix<int,Eigen::Dynamic,1> & I,
-  Eigen::Matrix<T,Eigen::Dynamic,1> & V)
-{
-  // Resize outputs to fit nonzeros
-  I.resize(X.nonZeros());
-  V.resize(X.nonZeros());
-
-  int i = 0;
-  // loop over non-zeros
-  for(typename Eigen::SparseVector<T>::InnerIterator it(X); it; ++it)
-  {
-    I(i) = it.index();
-    V(i) = it.value();
-    i++;
-  }
-}
+#ifdef IGL_HEADER_ONLY
+#  include "find.cpp"
+#endif
 
 #endif

full.cpp

@@ -0,0 +1,22 @@
+#include "full.h"
+
+template <typename T>
+IGL_INLINE void igl::full(
+  const Eigen::SparseMatrix<T> & A,
+  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B)
+{
+  B = Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::Zero(A.rows(),A.cols());
+  // Iterate over outside
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
+    {
+      B(it.row(),it.col()) = it.value();
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

full.h

@@ -1,5 +1,6 @@
 #ifndef IGL_FULL_H
 #define IGL_FULL_H
+#include "igl_inline.h"
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
@@ -13,27 +14,13 @@ namespace igl
   // Output:
   //   B  m by n dense/full matrix
   template <typename T>
-  inline void full(
+  IGL_INLINE void full(
     const Eigen::SparseMatrix<T> & A,
     Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B);
 }
 
-// Implementation
+#ifdef IGL_HEADER_ONLY
+#  include "full.cpp"
+#endif
 
-template <typename T>
-inline void igl::full(
-  const Eigen::SparseMatrix<T> & A,
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B)
-{
-  B = Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::Zero(A.rows(),A.cols());
-  // Iterate over outside
-  for(int k=0; k<A.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
-    {
-      B(it.row(),it.col()) = it.value();
-    }
-  }
-}
 #endif

get_seconds.cpp

@@ -0,0 +1,20 @@
+#include "get_seconds.h"
+
+#if _WIN32
+#  include <ctime>
+IGL_INLINE double igl::get_seconds()
+{
+  // This does not work on mac os x with glut in the main loop
+  return double(clock())/CLOCKS_PER_SEC;
+}
+#else
+#  include <sys/time.h>
+IGL_INLINE double igl::get_seconds()
+{
+  timeval time;
+  gettimeofday(&time, NULL);
+  return time.tv_sec + time.tv_usec / 1e6;
+  // This does not work on mac os x with glut in the main loop
+  //return double(clock())/CLOCKS_PER_SEC;
+}
+#endif

get_seconds.h

@@ -1,31 +1,16 @@
 #ifndef IGL_GET_SECONDS_H
 #define IGL_GET_SECONDS_H
+#include "igl_inline.h"
 
 namespace igl
 {
   // Return the current time in seconds since program start
-  inline double get_seconds();
+  IGL_INLINE double get_seconds();
 
 }
 
-// Implementation
-
-#if _WIN32
-#  include <ctime>
-inline double igl::get_seconds()
-{
-  // This does not work on mac os x with glut in the main loop
-  return double(clock())/CLOCKS_PER_SEC;
-}
-#else
-#  include <sys/time.h>
-inline double igl::get_seconds()
-{
-  timeval time;
-  gettimeofday(&time, NULL);
-  return time.tv_sec + time.tv_usec / 1e6;
-  // This does not work on mac os x with glut in the main loop
-  //return double(clock())/CLOCKS_PER_SEC;
-}
+#ifdef IGL_HEADER_ONLY
+#  include "get_seconds.cpp"
 #endif
+
 #endif

get_seconds_hires.cpp

@@ -0,0 +1,21 @@
+#include "get_seconds_hires.h"
+
+#if _WIN32
+#  include <windows.h>
+IGL_INLINE double igl::get_seconds_hires()
+{
+  LARGE_INTEGER li_freq, li_current;
+  const bool ret = QueryPerformanceFrequency(&li_freq);
+  const bool ret2 = QueryPerformanceCounter(&li_current);
+  assert(ret && ret2);
+  assert(li_freq.QuadPart > 0);
+  return double(li_current.QuadPart) / double(li_freq.QuadPart);
+}
+#else
+#  include "get_seconds.h"
+IGL_INLINE double igl::get_seconds_hires()
+{
+  // Sorry I've no idea how performance counters work on Mac...
+  return igl::get_seconds();
+}
+#endif

get_seconds_hires.h

@@ -1,31 +1,15 @@
 #ifndef IGL_GET_SECONDS_HIRES_H
 #define IGL_GET_SECONDS_HIRES_H
+#include "igl_inline.h"
 
 namespace igl
 {
   // Return the current time in seconds using performance counters
-  inline double get_seconds_hires();
+  IGL_INLINE double get_seconds_hires();
 }
 
-// Implementation
-
-#if _WIN32
-#  include <windows.h>
-inline double igl::get_seconds_hires()
-{
-  LARGE_INTEGER li_freq, li_current;
-  const bool ret = QueryPerformanceFrequency(&li_freq);
-  const bool ret2 = QueryPerformanceCounter(&li_current);
-  assert(ret && ret2);
-  assert(li_freq.QuadPart > 0);
-  return double(li_current.QuadPart) / double(li_freq.QuadPart);
-}
-#else
-#  include "get_seconds.h"
-inline double igl::get_seconds_hires()
-{
-  // Sorry I've no idea how performance counters work on Mac...
-  return igl::get_seconds();
-}
+#ifdef IGL_HEADER_ONLY
+#  include "get_seconds_hires.cpp"
 #endif
+
 #endif

gl_type_size.cpp

@@ -0,0 +1,23 @@
+#include "gl_type_size.h"
+#include <cassert>
+
+IGL_INLINE int igl::gl_type_size(const GLenum type)
+{
+  switch(type)
+  {
+    case GL_DOUBLE:
+      return 8;
+      break;
+    case GL_FLOAT:
+      return 4;
+      break;
+    case GL_INT:
+      return 4;
+      break;
+    default:
+      // should handle all other GL_[types]
+      assert(false);
+      break;
+  }
+  return -1;
+}

gl_type_size.h

@@ -1,5 +1,6 @@
 #ifndef IGL_GL_TYPE_SIZE_H
 #define IGL_GL_TYPE_SIZE_H
+#include "igl_inline.h"
 
 #if __APPLE__
 #  include <OpenGL/gl.h>
@@ -18,31 +19,11 @@ namespace igl
   // Inputs:
   //   type  enum value of opengl type
   // Returns size in bytes of type
-  inline int gl_type_size(const GLenum type);
+  IGL_INLINE int gl_type_size(const GLenum type);
 }
 
-// Implementation
-
-inline int igl::gl_type_size(const GLenum type)
-{
-  switch(type)
-  {
-    case GL_DOUBLE:
-      return 8;
-      break;
-    case GL_FLOAT:
-      return 4;
-      break;
-    case GL_INT:
-      return 4;
-      break;
-    default:
-      // should handle all other GL_[types]
-      assert(false);
-      break;
-  }
-  return -1;
-}
-
-#endif 
+#ifdef IGL_HEADER_ONLY
+#  include "gl_type_size.cpp"
+#endif
 
+#endif

grad.cpp

@@ -0,0 +1,51 @@
+#include "grad.h"
+
+template <typename T>
+IGL_INLINE void igl::grad(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
+  const Eigen::MatrixXi &F,
+  const Eigen::Matrix<T, Eigen::Dynamic, 1>&X,
+  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &G )
+{
+  G = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>::Zero(F.rows(),3);
+  for (int i = 0; i<F.rows(); ++i)
+  {
+    // renaming indices of vertices of triangles for convenience
+    int i1 = F(i,0);
+    int i2 = F(i,1);
+    int i3 = F(i,2);
+    
+    // #F x 3 matrices of triangle edge vectors, named after opposite vertices
+    Eigen::Matrix<T, 1, 3> v32 = V.row(i3) - V.row(i2);
+    Eigen::Matrix<T, 1, 3> v13 = V.row(i1) - V.row(i3);
+    Eigen::Matrix<T, 1, 3> v21 = V.row(i2) - V.row(i1);
+    
+    // area of parallelogram is twice area of triangle
+    // area of parallelogram is || v1 x v2 || 
+    Eigen::Matrix<T, 1, 3> n  = v32.cross(v13); 
+    
+    // This does correct l2 norm of rows, so that it contains #F list of twice
+    // triangle areas
+    double dblA = std::sqrt(n.dot(n));
+    
+    // now normalize normals to get unit normals
+    Eigen::Matrix<T, 1, 3> u = n / dblA;
+    
+    // rotate each vector 90 degrees around normal
+    double norm21 = std::sqrt(v21.dot(v21));
+    double norm13 = std::sqrt(v13.dot(v13));
+    Eigen::Matrix<T, 1, 3> eperp21 = u.cross(v21);
+    eperp21 = eperp21 / std::sqrt(eperp21.dot(eperp21));
+    eperp21 *= norm21;
+    Eigen::Matrix<T, 1, 3> eperp13 = u.cross(v13);
+    eperp13 = eperp13 / std::sqrt(eperp13.dot(eperp13));
+    eperp13 *= norm13;
+    
+    G.row(i) = ((X[i2] -X[i1]) *eperp13 + (X[i3] - X[i1]) *eperp21) / dblA;
+  };
+}
+  
+  
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

grad.h

@@ -8,6 +8,7 @@
 
 #ifndef IGL_GRAD_H
 #define IGL_GRAD_H
+#include "igl_inline.h"
 
 #include <Eigen/Core>
 
@@ -33,57 +34,14 @@ namespace igl {
   // 90 degrees
   //
   template <typename T>
-  inline void grad(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
+  IGL_INLINE void grad(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
     const Eigen::MatrixXi &F,
     const Eigen::Matrix<T, Eigen::Dynamic, 1>&X,
     Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &G );
 }
 
-// Implementation
+#ifdef IGL_HEADER_ONLY
+#  include "grad.cpp"
+#endif
 
-template <typename T>
-inline void igl::grad(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
-  const Eigen::MatrixXi &F,
-  const Eigen::Matrix<T, Eigen::Dynamic, 1>&X,
-  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &G )
-{
-  G = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>::Zero(F.rows(),3);
-  for (int i = 0; i<F.rows(); ++i)
-  {
-    // renaming indices of vertices of triangles for convenience
-    int i1 = F(i,0);
-    int i2 = F(i,1);
-    int i3 = F(i,2);
-    
-    // #F x 3 matrices of triangle edge vectors, named after opposite vertices
-    Eigen::Matrix<T, 1, 3> v32 = V.row(i3) - V.row(i2);
-    Eigen::Matrix<T, 1, 3> v13 = V.row(i1) - V.row(i3);
-    Eigen::Matrix<T, 1, 3> v21 = V.row(i2) - V.row(i1);
-    
-    // area of parallelogram is twice area of triangle
-    // area of parallelogram is || v1 x v2 || 
-    Eigen::Matrix<T, 1, 3> n  = v32.cross(v13); 
-    
-    // This does correct l2 norm of rows, so that it contains #F list of twice
-    // triangle areas
-    double dblA = std::sqrt(n.dot(n));
-    
-    // now normalize normals to get unit normals
-    Eigen::Matrix<T, 1, 3> u = n / dblA;
-    
-    // rotate each vector 90 degrees around normal
-    double norm21 = std::sqrt(v21.dot(v21));
-    double norm13 = std::sqrt(v13.dot(v13));
-    Eigen::Matrix<T, 1, 3> eperp21 = u.cross(v21);
-    eperp21 = eperp21 / std::sqrt(eperp21.dot(eperp21));
-    eperp21 *= norm21;
-    Eigen::Matrix<T, 1, 3> eperp13 = u.cross(v13);
-    eperp13 = eperp13 / std::sqrt(eperp13.dot(eperp13));
-    eperp13 *= norm13;
-    
-    G.row(i) = ((X[i2] -X[i1]) *eperp13 + (X[i3] - X[i1]) *eperp21) / dblA;
-  };
-}
-  
-  
 #endif

is_border_vertex.cpp

@@ -0,0 +1,27 @@
+#include "is_border_vertex.h"
+#include <vector>
+
+#include "tt.h"
+
+template<typename T>
+IGL_INLINE std::vector<bool> igl::is_border_vertex(const T& V, const Eigen::MatrixXi& F)
+{
+  Eigen::MatrixXi FF;
+  igl::tt(V,F,FF);
+  std::vector<bool> ret(V.rows());
+  for(int i=0; i<ret.size();++i)
+    ret[i] = false;
+  
+  for(int i=0; i<F.rows();++i)
+    for(int j=0;j<F.cols();++j)
+      if(FF(i,j) == -1)
+      {
+        ret[F(i,j)]       = true;
+        ret[F(i,(j+1)%3)] = true;
+      }
+  return ret;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

is_border_vertex.h

@@ -5,6 +5,7 @@
 
 #ifndef IGL_IS_BORDER_VERTEX_H
 #define IGL_IS_BORDER_VERTEX_H
+#include "igl_inline.h"
 
 #include <Eigen/Core>
 #include <vector>
@@ -12,30 +13,11 @@
 namespace igl 
 {
   template<typename T>
-  inline std::vector<bool> is_border_vertex(const T& V, const Eigen::MatrixXi& F)
-}
-
-// Implementation
-#include "tt.h"
-
-template<typename T>
-inline std::vector<bool> igl::is_border_vertex(const T& V, const Eigen::MatrixXi& F)
-{
-  Eigen::MatrixXi FF;
-  igl::tt(V,F,FF);
-  vector<bool> ret(V.rows());
-  for(int i=0; i<ret.size();++i)
-    ret[i] = false;
-  
-  for(int i=0; i<F.rows();++i)
-    for(int j=0;j<F.cols();++j)
-      if(FF(i,j) == -1)
-      {
-        ret[F(i,j)]       = true;
-        ret[F(i,(j+1)%3)] = true;
-      }
-  return ret;
+  IGL_INLINE std::vector<bool> is_border_vertex(const T& V, const Eigen::MatrixXi& F);
 }
 
+#ifdef IGL_HEADER_ONLY
+#  include "is_border_vertex.cpp"
+#endif
 
 #endif

is_dir.cpp

@@ -0,0 +1,23 @@
+#include "is_dir.h"
+
+#include <sys/stat.h>
+
+#ifndef S_ISDIR
+#define S_ISDIR(mode)  (((mode) & S_IFMT) == S_IFDIR)
+#endif
+
+#ifndef S_ISREG
+#define S_ISREG(mode)  (((mode) & S_IFMT) == S_IFREG)
+#endif
+
+IGL_INLINE bool igl::is_dir(const char * filename)
+{
+  struct stat status;
+  if(stat(filename,&status)!=0)
+  {
+    // path does not exist
+    return false;
+  }
+  // Tests whether existing path is a directory
+  return S_ISDIR(status.st_mode);
+}

is_dir.h

@@ -1,5 +1,6 @@
 #ifndef IGL_IS_DIR_H
 #define IGL_IS_DIR_H
+#include "igl_inline.h"
 namespace igl
 {
   // Act like php's is_dir function
@@ -10,31 +11,12 @@ namespace igl
   //     be checked relative to the current working directory. 
   // Returns TRUE if the filename exists and is a directory, FALSE
   // otherwise.
-  inline bool is_dir(const char * filename);
+  IGL_INLINE bool is_dir(const char * filename);
 
 }
 
-// Implementation
-#include <sys/stat.h>
-
-#ifndef S_ISDIR
-#define S_ISDIR(mode)  (((mode) & S_IFMT) == S_IFDIR)
-#endif
-
-#ifndef S_ISREG
-#define S_ISREG(mode)  (((mode) & S_IFMT) == S_IFREG)
+#ifdef IGL_HEADER_ONLY
+#  include "is_dir.cpp"
 #endif
 
-inline bool igl::is_dir(const char * filename)
-{
-  struct stat status;
-  if(stat(filename,&status)!=0)
-  {
-    // path does not exist
-    return false;
-  }
-  // Tests whether existing path is a directory
-  return S_ISDIR(status.st_mode);
-}
-
 #endif

is_file.cpp

@@ -0,0 +1,14 @@
+#include "is_file.h"
+
+#include <sys/stat.h>
+IGL_INLINE bool igl::is_file(const char * filename)
+{
+  struct stat status;
+  if(stat(filename,&status)!=0)
+  {
+    // path does not exist
+    return false;
+  }
+  // Tests whether existing path is a regular file
+  return S_ISREG(status.st_mode);
+}

is_file.h

@@ -1,5 +1,6 @@
 #ifndef IGL_IS_FILE_H
 #define IGL_IS_FILE_H
+#include "igl_inline.h"
 namespace igl
 {
   // Act like php's is_file function
@@ -10,22 +11,12 @@ namespace igl
   //     be checked relative to the current working directory. 
   // Returns TRUE if the filename exists and is a regular file, FALSE
   // otherwise.
-  inline bool is_file(const char * filename);
+  IGL_INLINE bool is_file(const char * filename);
 
 }
 
-// Implementation
-#include <sys/stat.h>
-inline bool igl::is_file(const char * filename)
-{
-  struct stat status;
-  if(stat(filename,&status)!=0)
-  {
-    // path does not exist
-    return false;
-  }
-  // Tests whether existing path is a regular file
-  return S_ISREG(status.st_mode);
-}
+#ifdef IGL_HEADER_ONLY
+#  include "is_file.cpp"
+#endif
 
 #endif

is_manifold.cpp

@@ -0,0 +1,40 @@
+#include "is_manifold.h"
+
+#include <algorithm>
+
+template<typename T>
+IGL_INLINE bool igl::is_manifold(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& V, const Eigen::MatrixXi& F)
+  {
+    std::vector<std::vector<int> > TTT;
+    for(int f=0;f<F.rows();++f)
+      for (int i=0;i<3;++i)
+      {
+        // v1 v2 f ei 
+        int v1 = F(f,i);
+        int v2 = F(f,(i+1)%3);
+        if (v1 > v2) std::swap(v1,v2);
+        std::vector<int> r(4);
+        r[0] = v1; r[1] = v2;
+        r[2] = f;  r[3] = i;
+        TTT.push_back(r);
+      }
+    std::sort(TTT.begin(),TTT.end());
+    
+    for(int i=2;i<(int)TTT.size();++i)
+    {
+      std::vector<int>& r1 = TTT[i-2];
+      std::vector<int>& r2 = TTT[i-1];
+      std::vector<int>& r3 = TTT[i];
+      if ( (r1[0] == r2[0] && r2[0] == r3[0]) 
+          && 
+          (r1[1] == r2[1] && r2[1] == r3[1]) )
+      {
+        return false;
+      }
+    }
+    return true;
+  }
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

is_manifold.h

@@ -5,6 +5,7 @@
 
 #ifndef IGL_IS_MANIFOLD_H
 #define IGL_IS_MANIFOLD_H
+#include "igl_inline.h"
 
 #include <Eigen/Core>
 #include <vector>
@@ -18,44 +19,11 @@ namespace igl
   // Known Bugs:
   //  Does not check for non-manifold vertices
   template<typename T>
-  inline bool is_manifold(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& V, const Eigen::MatrixXi& F);
+  IGL_INLINE bool is_manifold(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& V, const Eigen::MatrixXi& F);
 }
 
-// Implementation
-#include <algorithm>
-
-template<typename T>
-inline bool igl::is_manifold(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& V, const Eigen::MatrixXi& F);
-  {
-    std::vector<std::vector<int> > TTT;
-    for(int f=0;f<F.rows();++f)
-      for (int i=0;i<3;++i)
-      {
-        // v1 v2 f ei 
-        int v1 = F(f,i);
-        int v2 = F(f,(i+1)%3);
-        if (v1 > v2) std::swap(v1,v2);
-        std::vector<int> r(4);
-        r[0] = v1; r[1] = v2;
-        r[2] = f;  r[3] = i;
-        TTT.push_back(r);
-      }
-    std::sort(TTT.begin(),TTT.end());
-    
-    for(int i=2;i<(int)TTT.size();++i)
-    {
-      std::vector<int>& r1 = TTT[i-2];
-      std::vector<int>& r2 = TTT[i-1];
-      std::vector<int>& r3 = TTT[i];
-      if ( (r1[0] == r2[0] && r2[0] == r3[0]) 
-          && 
-          (r1[1] == r2[1] && r2[1] == r3[1]) )
-      {
-        return false;
-      }
-    }
-    return true;
-  }
-}
+#ifdef IGL_HEADER_ONLY
+#  include "is_manifold.cpp"
+#endif
 
 #endif

is_readable.cpp

@@ -0,0 +1,47 @@
+#include "is_readable.h"
+
+#ifdef _WIN32
+#  include <cstdio>
+IGL_INLINE bool igl::is_readable(const char* filename)
+{
+  FILE * f = fopen(filename,"r");
+  if(f == NULL)
+  {
+    return false;
+  }
+  fclose(f);
+  return true;
+}
+#else
+#  include <sys/stat.h>
+#  include <unistd.h>
+IGL_INLINE bool igl::is_readable(const char* filename)
+{
+  // Check if file already exists
+  struct stat status;
+  if(stat(filename,&status)!=0)
+  {
+    return false;
+  }
+
+  // Get current users uid and gid
+  uid_t this_uid = getuid();
+  gid_t this_gid = getgid();
+
+  // Dealing with owner
+  if( this_uid == status.st_uid )
+  {
+    return S_IRUSR & status.st_mode;
+  }
+
+  // Dealing with group member
+  if( this_gid == status.st_gid )
+  {
+    return S_IRGRP & status.st_mode;
+  }
+
+  // Dealing with other
+  return S_IROTH & status.st_mode;
+
+}
+#endif

is_readable.h

@@ -1,5 +1,6 @@
 #ifndef IGL_IS_READABLE_H
 #define IGL_IS_READABLE_H
+#include "igl_inline.h"
 namespace igl
 {
   // Check if a file is reabable like PHP's is_readable function:
@@ -10,55 +11,11 @@ namespace igl
   // exist or *is not readable*
   //
   // Note: Windows version will not check user or group ids
-  inline bool is_readable(const char * filename);
+  IGL_INLINE bool is_readable(const char * filename);
 }
 
-
-// Implementation
-#ifdef _WIN32
-#  include <cstdio>
-inline bool igl::is_readable(const char* filename)
-{
-  FILE * f = fopen(filename,"r");
-  if(f == NULL)
-  {
-    return false;
-  }
-  fclose(f);
-  return true;
-}
-#else
-#  include <sys/stat.h>
-#  include <unistd.h>
-inline bool igl::is_readable(const char* filename)
-{
-  // Check if file already exists
-  struct stat status;
-  if(stat(filename,&status)!=0)
-  {
-    return false;
-  }
-
-  // Get current users uid and gid
-  uid_t this_uid = getuid();
-  gid_t this_gid = getgid();
-
-  // Dealing with owner
-  if( this_uid == status.st_uid )
-  {
-    return S_IRUSR & status.st_mode;
-  }
-
-  // Dealing with group member
-  if( this_gid == status.st_gid )
-  {
-    return S_IRGRP & status.st_mode;
-  }
-
-  // Dealing with other
-  return S_IROTH & status.st_mode;
-
-}
+#ifdef IGL_HEADER_ONLY
+#  include "is_readable.cpp"
 #endif
 
 #endif

is_symmetric.cpp

@@ -0,0 +1,21 @@
+#include "is_symmetric.h"
+
+template <typename T>
+IGL_INLINE bool igl::is_symmetric(const Eigen::SparseMatrix<T>& A)
+{
+  if(A.rows() != A.cols())
+  {
+    return false;
+  }
+  Eigen::SparseMatrix<T> AT = A.transpose();
+  Eigen::SparseMatrix<T> AmAT = A-AT;
+  //// Eigen screws up something with LLT if you try to do
+  //SparseMatrix<T> AmAT = A-A.transpose();
+  //// Eigen crashes at runtime if you try to do
+  // return (A-A.transpose()).nonZeros() == 0;
+  return AmAT.nonZeros() == 0;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

is_symmetric.h

@@ -1,5 +1,8 @@
 #ifndef IGL_IS_SYMMETRIC_H
 #define IGL_IS_SYMMETRIC_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Sparse>
 namespace igl
 {
   // Returns true if the given matrix is symmetric
@@ -7,24 +10,11 @@ namespace igl
   //   A  m by m matrix
   // Returns true if the matrix is square and symmetric
   template <typename T>
-  inline bool is_symmetric(const Eigen::SparseMatrix<T>& A);
+  IGL_INLINE bool is_symmetric(const Eigen::SparseMatrix<T>& A);
 }
 
-// Implementation
+#ifdef IGL_HEADER_ONLY
+#  include "is_symmetric.cpp"
+#endif
 
-template <typename T>
-inline bool igl::is_symmetric(const Eigen::SparseMatrix<T>& A)
-{
-  if(A.rows() != A.cols())
-  {
-    return false;
-  }
-  Eigen::SparseMatrix<T> AT = A.transpose();
-  Eigen::SparseMatrix<T> AmAT = A-AT;
-  //// Eigen screws up something with LLT if you try to do
-  //SparseMatrix<T> AmAT = A-A.transpose();
-  //// Eigen crashes at runtime if you try to do
-  // return (A-A.transpose()).nonZeros() == 0;
-  return AmAT.nonZeros() == 0;
-}
 #endif

is_writable.cpp

@@ -0,0 +1,51 @@
+#include "is_writable.h"
+
+#ifdef _WIN32
+#include <sys/stat.h>
+#ifndef S_IWUSR
+#  define S_IWUSR S_IWRITE
+#endif
+IGL_INLINE bool is_writable(const char* filename)
+{
+  // Check if file already exists
+  struct stat status;
+  if(stat(filename,&status)!=0)
+  {
+    return false;
+  }
+
+  return S_IWUSR & status.st_mode;
+}
+#else
+#include <sys/stat.h>
+#include <unistd.h>
+
+IGL_INLINE bool igl::is_writable(const char* filename)
+{
+  // Check if file already exists
+  struct stat status;
+  if(stat(filename,&status)!=0)
+  {
+    return false;
+  }
+
+  // Get current users uid and gid
+  uid_t this_uid = getuid();
+  gid_t this_gid = getgid();
+
+  // Dealing with owner
+  if( this_uid == status.st_uid )
+  {
+    return S_IWUSR & status.st_mode;
+  }
+
+  // Dealing with group member
+  if( this_gid == status.st_gid )
+  {
+    return S_IWGRP & status.st_mode;
+  }
+
+  // Dealing with other
+  return S_IWOTH & status.st_mode;
+}
+#endif

is_writable.h

@@ -1,5 +1,6 @@
 #ifndef IGL_IS_WRITABLE_H
 #define IGL_IS_WRITABLE_H
+#include "igl_inline.h"
 namespace igl
 {
   // Check if a file exists *and* is writable like PHP's is_writable function:
@@ -10,58 +11,11 @@ namespace igl
   // exist or *is not writable*
   //
   // Note: Windows version will not test group and user id
-  inline bool is_writable(const char * filename);
+  IGL_INLINE bool is_writable(const char * filename);
 }
 
-
-// Implementation
-#ifdef _WIN32
-#include <sys/stat.h>
-#ifndef S_IWUSR
-#  define S_IWUSR S_IWRITE
+#ifdef IGL_HEADER_ONLY
+#  include "is_writable.cpp"
 #endif
-inline bool is_writable(const char* filename)
-{
-  // Check if file already exists
-  struct stat status;
-  if(stat(filename,&status)!=0)
-  {
-    return false;
-  }
-
-  return S_IWUSR & status.st_mode;
-}
-#else
-#include <sys/stat.h>
-#include <unistd.h>
-
-inline bool igl::is_writable(const char* filename)
-{
-  // Check if file already exists
-  struct stat status;
-  if(stat(filename,&status)!=0)
-  {
-    return false;
-  }
-
-  // Get current users uid and gid
-  uid_t this_uid = getuid();
-  gid_t this_gid = getgid();
-
-  // Dealing with owner
-  if( this_uid == status.st_uid )
-  {
-    return S_IWUSR & status.st_mode;
-  }
 
-  // Dealing with group member
-  if( this_gid == status.st_gid )
-  {
-    return S_IWGRP & status.st_mode;
-  }
-
-  // Dealing with other
-  return S_IWOTH & status.st_mode;
-}
 #endif
-#endif

limit_faces.cpp

@@ -0,0 +1,56 @@
+#include "limit_faces.h"
+
+#include <vector>
+#include <Eigen/Dense>
+
+template <typename MatF, typename VecL>
+IGL_INLINE void igl::limit_faces(
+  const MatF & F, 
+  const VecL & L, 
+  const bool exclusive,
+  MatF & LF)
+{
+  using namespace std;
+  using namespace Eigen;
+  vector<bool> in(F.rows(),false);
+  int num_in = 0;
+  // loop over faces
+  for(int i = 0;i<F.rows();i++)
+  {
+    bool all = true;
+    bool any = false;
+    for(int j = 0;j<F.cols();j++)
+    {
+      bool found = false;
+      // loop over L
+      for(int l = 0;l<L.size();l++)
+      {
+        if(F(i,j) == L(l))
+        {
+          found = true;
+          break;
+        }
+      }
+      any |= found;
+      all &= found;
+    }
+    in[i] = (exclusive?all:any);
+    num_in += (in[i]?1:0);
+  }
+
+  LF.resize(num_in,F.cols());
+  // loop over faces
+  int lfi = 0;
+  for(int i = 0;i<F.rows();i++)
+  {
+    if(in[i])
+    {
+      LF.row(lfi) = F.row(i);
+      lfi++;
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

limit_faces.h

@@ -1,5 +1,6 @@
 #ifndef IGL_LIMIT_FACES_H
 #define IGL_LIMIT_FACES_H
+#include "igl_inline.h"
 namespace igl
 {
   // LIMIT_FACES limit given faces F to those which contain (only) indices found
@@ -21,62 +22,15 @@ namespace igl
   //   in  #F list of whether given face was included
   //
   template <typename MatF, typename VecL>
-  inline void limit_faces(
+  IGL_INLINE void limit_faces(
     const MatF & F, 
     const VecL & L, 
     const bool exclusive,
     MatF & LF);
 }
 
-// Implementation
-#include <vector>
-
-template <typename MatF, typename VecL>
-inline void igl::limit_faces(
-  const MatF & F, 
-  const VecL & L, 
-  const bool exclusive,
-  MatF & LF)
-{
-  using namespace std;
-  using namespace Eigen;
-  vector<bool> in(F.rows(),false);
-  int num_in = 0;
-  // loop over faces
-  for(int i = 0;i<F.rows();i++)
-  {
-    bool all = true;
-    bool any = false;
-    for(int j = 0;j<F.cols();j++)
-    {
-      bool found = false;
-      // loop over L
-      for(int l = 0;l<L.size();l++)
-      {
-        if(F(i,j) == L(l))
-        {
-          found = true;
-          break;
-        }
-      }
-      any |= found;
-      all &= found;
-    }
-    in[i] = (exclusive?all:any);
-    num_in += (in[i]?1:0);
-  }
-
-  LF.resize(num_in,F.cols());
-  // loop over faces
-  int lfi = 0;
-  for(int i = 0;i<F.rows();i++)
-  {
-    if(in[i])
-    {
-      LF.row(lfi) = F.row(i);
-      lfi++;
-    }
-  }
-}
+#ifdef IGL_HEADER_ONLY
+#  include "limit_faces.cpp"
+#endif
 
 #endif

list_to_matrix.cpp

@@ -0,0 +1,48 @@
+#include "list_to_matrix.h"
+
+#include <cassert>
+#include <cstdio>
+
+#include "max_size.h"
+#include "min_size.h"
+#define VERBOSE
+#include "verbose.h"
+
+template <typename T, class Mat>
+IGL_INLINE bool igl::list_to_matrix(const std::vector<std::vector<T > > & V,Mat & M)
+{
+  // number of columns
+  int m = V.size();
+  if(m == 0)
+  {
+    fprintf(stderr,"Error: list_to_matrix() list is empty()\n");
+    return false;
+  }
+  // number of rows
+  int n = igl::min_size(V);
+  if(n != igl::max_size(V))
+  {
+    fprintf(stderr,"Error: list_to_matrix()"
+      " list elements are not all the same size\n");
+    return false;
+  }
+  assert(n != -1);
+  // Resize output
+  M.resize(m,n);
+
+  // Loop over rows
+  for(int i = 0;i<m;i++)
+  {
+    // Loop over cols
+    for(int j = 0;j<n;j++)
+    {
+      M(i,j) = V[i][j];
+    }
+  }
+
+  return true;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

list_to_matrix.h

@@ -1,5 +1,6 @@
 #ifndef IGL_LIST_TO_MATRIX_H
 #define IGL_LIST_TO_MATRIX_H
+#include "igl_inline.h"
 #include <vector>
 namespace igl
 {
@@ -15,50 +16,11 @@ namespace igl
   //   M  an m by n matrix
   // Returns true on success, false on errors
   template <typename T, class Mat>
-  inline bool list_to_matrix(const std::vector<std::vector<T > > & V,Mat & M);
+  IGL_INLINE bool list_to_matrix(const std::vector<std::vector<T > > & V,Mat & M);
 }
 
-// Implementation
-#include <cassert>
-#include <cstdio>
-
-#include "max_size.h"
-#include "min_size.h"
-#define VERBOSE
-#include "verbose.h"
-
-template <typename T, class Mat>
-inline bool igl::list_to_matrix(const std::vector<std::vector<T > > & V,Mat & M)
-{
-  // number of columns
-  int m = V.size();
-  if(m == 0)
-  {
-    fprintf(stderr,"Error: list_to_matrix() list is empty()\n");
-    return false;
-  }
-  // number of rows
-  int n = igl::min_size(V);
-  if(n != igl::max_size(V))
-  {
-    fprintf(stderr,"Error: list_to_matrix()"
-      " list elements are not all the same size\n");
-    return false;
-  }
-  assert(n != -1);
-  // Resize output
-  M.resize(m,n);
-
-  // Loop over rows
-  for(int i = 0;i<m;i++)
-  {
-    // Loop over cols
-    for(int j = 0;j<n;j++)
-    {
-      M(i,j) = V[i][j];
-    }
-  }
+#ifdef IGL_HEADER_ONLY
+#  include "list_to_matrix.cpp"
+#endif
 
-  return true;
-}
 #endif

load_shader.cpp

@@ -0,0 +1,20 @@
+#include "load_shader.h"
+
+// Copyright Denis Kovacs 4/10/08
+#include "print_shader_info_log.h"
+#include <cstdio>
+IGL_INLINE GLuint igl::load_shader(const char *src,const GLenum type)
+{
+  GLuint s = glCreateShader(type);
+  if(s == 0)
+  {
+    fprintf(stderr,"Error: load_shader() failed to create shader.\n");
+    return 0;
+  }
+  // Pass shader source string
+  glShaderSource(s, 1, &src, NULL);
+  glCompileShader(s);
+  // Print info log (if any)
+  igl::print_shader_info_log(s);
+  return s;
+}

load_shader.h

@@ -1,5 +1,6 @@
 #ifndef IGL_LOAD_SHADER_H 
-#define IGL_LOAD_SHADER_H 
+#define IGL_LOAD_SHADER_H
+#include "igl_inline.h" 
 
 #ifdef __APPLE__
 #  include <OpenGL/gl.h>
@@ -22,29 +23,11 @@ namespace igl
   //     GL_FRAGMENT_SHADER
   //     GL_GEOMETRY_SHADER
   // Returns  index id of the newly created shader, 0 on error
-  inline GLuint load_shader(const char *src,const GLenum type);
-}
-
-// Implementation
-
-// Copyright Denis Kovacs 4/10/08
-#include "print_shader_info_log.h"
-#include <cstdio>
-inline GLuint igl::load_shader(const char *src,const GLenum type)
-{
-  GLuint s = glCreateShader(type);
-  if(s == 0)
-  {
-    fprintf(stderr,"Error: load_shader() failed to create shader.\n");
-    return 0;
-  }
-  // Pass shader source string
-  glShaderSource(s, 1, &src, NULL);
-  glCompileShader(s);
-  // Print info log (if any)
-  igl::print_shader_info_log(s);
-  return s;
+  IGL_INLINE GLuint load_shader(const char *src,const GLenum type);
 }
 
+#ifdef IGL_HEADER_ONLY
+#  include "load_shader.cpp"
 #endif
 
+#endif

lu_lagrange.cpp

@@ -0,0 +1,137 @@
+#include "lu_lagrange.h"
+
+// Cholesky LLT decomposition for symmetric positive definite
+#include <Eigen/SparseExtra>
+#include <cassert>
+#include "find.h"
+#include "sparse.h"
+
+template <typename T>
+IGL_INLINE bool igl::lu_lagrange(
+  const Eigen::SparseMatrix<T> & ATA,
+  const Eigen::SparseMatrix<T> & C,
+  Eigen::SparseMatrix<T> & L,
+  Eigen::SparseMatrix<T> & U)
+{
+  // number of unknowns
+  int n = ATA.rows();
+  // number of lagrange multipliers
+  int m = C.cols();
+
+  assert(ATA.cols() == n);
+  if(m != 0)
+  {
+    assert(C.rows() == n);
+    if(C.nonZeros() == 0)
+    {
+      // See note above about empty columns in C
+      fprintf(stderr,"Error: lu_lagrange() C has columns but no entries\n");
+      return false;
+    }
+  }
+
+  // Check that each column of C has at least one entry
+  std::vector<bool> has_entry; has_entry.resize(C.cols(),false);
+  // Iterate over outside
+  for(int k=0; k<C.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (C,k); it; ++it)
+    {
+      has_entry[it.col()] = true;
+    }
+  }
+  for(int i=0;i<(int)has_entry.size();i++)
+  {
+    if(!has_entry[i])
+    {
+      // See note above about empty columns in C
+      fprintf(stderr,"Error: lu_lagrange() C(:,%d) has no entries\n",i);
+      return false;
+    }
+  }
+
+
+
+  // Cholesky factorization of ATA
+  //// Eigen fails if you give a full view of the matrix like this:
+  //Eigen::SparseLLT<SparseMatrix<T> > ATA_LLT(ATA);
+  Eigen::SparseMatrix<T> ATA_LT = ATA.template triangularView<Eigen::Lower>();
+  Eigen::SparseLLT<Eigen::SparseMatrix<T> > ATA_LLT(ATA_LT);
+
+  Eigen::SparseMatrix<T> J = ATA_LLT.matrixL();
+
+  //if(!ATA_LLT.succeeded())
+  if(!((J*0).eval().nonZeros() == 0))
+  {
+    fprintf(stderr,"Error: lu_lagrange() failed to factor ATA\n");
+    return false;
+  }
+
+  if(m == 0)
+  {
+    // If there are no constraints (C is empty) then LU decomposition is just L
+    // and L' from cholesky decomposition
+    L = J;
+    U = J.transpose();
+  }else
+  {
+    // Construct helper matrix M
+    Eigen::SparseMatrix<T> M = C;
+    J.template triangularView<Eigen::Lower>().solveInPlace(M);
+
+    // Compute cholesky factorizaiton of M'*M
+    Eigen::SparseMatrix<T> MTM = M.transpose() * M;
+
+    Eigen::SparseLLT<Eigen::SparseMatrix<T> > MTM_LLT(MTM.template triangularView<Eigen::Lower>());
+
+    Eigen::SparseMatrix<T> K = MTM_LLT.matrixL();
+
+    //if(!MTM_LLT.succeeded())
+    if(!((K*0).eval().nonZeros() == 0))
+    {
+      fprintf(stderr,"Error: lu_lagrange() failed to factor MTM\n");
+      return false;
+    }
+
+    // assemble LU decomposition of Q
+    Eigen::Matrix<int,Eigen::Dynamic,1> MI;
+    Eigen::Matrix<int,Eigen::Dynamic,1> MJ;
+    Eigen::Matrix<T,Eigen::Dynamic,1> MV;
+    igl::find(M,MI,MJ,MV);
+
+    Eigen::Matrix<int,Eigen::Dynamic,1> KI;
+    Eigen::Matrix<int,Eigen::Dynamic,1> KJ;
+    Eigen::Matrix<T,Eigen::Dynamic,1> KV;
+    igl::find(K,KI,KJ,KV);
+
+    Eigen::Matrix<int,Eigen::Dynamic,1> JI;
+    Eigen::Matrix<int,Eigen::Dynamic,1> JJ;
+    Eigen::Matrix<T,Eigen::Dynamic,1> JV;
+    igl::find(J,JI,JJ,JV);
+
+    int nnz = JV.size()  + MV.size() + KV.size();
+
+    Eigen::Matrix<int,Eigen::Dynamic,1> UI(nnz);
+    Eigen::Matrix<int,Eigen::Dynamic,1> UJ(nnz);
+    Eigen::Matrix<T,Eigen::Dynamic,1> UV(nnz);
+    UI << JJ,                        MI, (KJ.array() + n).matrix();
+    UJ << JI, (MJ.array() + n).matrix(), (KI.array() + n).matrix(); 
+    UV << JV,                        MV,                     KV*-1;
+    igl::sparse(UI,UJ,UV,U);
+
+    Eigen::Matrix<int,Eigen::Dynamic,1> LI(nnz);
+    Eigen::Matrix<int,Eigen::Dynamic,1> LJ(nnz);
+    Eigen::Matrix<T,Eigen::Dynamic,1> LV(nnz);
+    LI << JI, (MJ.array() + n).matrix(), (KI.array() + n).matrix();
+    LJ << JJ,                        MI, (KJ.array() + n).matrix(); 
+    LV << JV,                        MV,                        KV;
+    igl::sparse(LI,LJ,LV,L);
+  }
+
+  return true;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

lu_lagrange.h

@@ -1,5 +1,6 @@
 #ifndef IGL_LU_LAGRANGE_H
 #define IGL_LU_LAGRANGE_H
+#include "igl_inline.h"
 
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
@@ -42,7 +43,7 @@ namespace igl
   // so you will get an error.
   //
   template <typename T>
-  inline bool lu_lagrange(
+  IGL_INLINE bool lu_lagrange(
     const Eigen::SparseMatrix<T> & ATA,
     const Eigen::SparseMatrix<T> & C,
     Eigen::SparseMatrix<T> & L,
@@ -50,137 +51,8 @@ namespace igl
 
 }
 
-// Implementation
-// Cholesky LLT decomposition for symmetric positive definite
-#include <Eigen/SparseExtra>
-#include <cassert>
-#include "find.h"
-#include "sparse.h"
-
-template <typename T>
-inline bool igl::lu_lagrange(
-  const Eigen::SparseMatrix<T> & ATA,
-  const Eigen::SparseMatrix<T> & C,
-  Eigen::SparseMatrix<T> & L,
-  Eigen::SparseMatrix<T> & U)
-{
-  // number of unknowns
-  int n = ATA.rows();
-  // number of lagrange multipliers
-  int m = C.cols();
-
-  assert(ATA.cols() == n);
-  if(m != 0)
-  {
-    assert(C.rows() == n);
-    if(C.nonZeros() == 0)
-    {
-      // See note above about empty columns in C
-      fprintf(stderr,"Error: lu_lagrange() C has columns but no entries\n");
-      return false;
-    }
-  }
-
-  // Check that each column of C has at least one entry
-  std::vector<bool> has_entry; has_entry.resize(C.cols(),false);
-  // Iterate over outside
-  for(int k=0; k<C.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (C,k); it; ++it)
-    {
-      has_entry[it.col()] = true;
-    }
-  }
-  for(int i=0;i<(int)has_entry.size();i++)
-  {
-    if(!has_entry[i])
-    {
-      // See note above about empty columns in C
-      fprintf(stderr,"Error: lu_lagrange() C(:,%d) has no entries\n",i);
-      return false;
-    }
-  }
-
-
-
-  // Cholesky factorization of ATA
-  //// Eigen fails if you give a full view of the matrix like this:
-  //Eigen::SparseLLT<SparseMatrix<T> > ATA_LLT(ATA);
-  Eigen::SparseMatrix<T> ATA_LT = ATA.template triangularView<Eigen::Lower>();
-  Eigen::SparseLLT<Eigen::SparseMatrix<T> > ATA_LLT(ATA_LT);
-
-  Eigen::SparseMatrix<T> J = ATA_LLT.matrixL();
-
-  //if(!ATA_LLT.succeeded())
-  if(!((J*0).eval().nonZeros() == 0))
-  {
-    fprintf(stderr,"Error: lu_lagrange() failed to factor ATA\n");
-    return false;
-  }
-
-  if(m == 0)
-  {
-    // If there are no constraints (C is empty) then LU decomposition is just L
-    // and L' from cholesky decomposition
-    L = J;
-    U = J.transpose();
-  }else
-  {
-    // Construct helper matrix M
-    Eigen::SparseMatrix<T> M = C;
-    J.template triangularView<Eigen::Lower>().solveInPlace(M);
-
-    // Compute cholesky factorizaiton of M'*M
-    Eigen::SparseMatrix<T> MTM = M.transpose() * M;
-
-    Eigen::SparseLLT<Eigen::SparseMatrix<T> > MTM_LLT(MTM.template triangularView<Eigen::Lower>());
-
-    Eigen::SparseMatrix<T> K = MTM_LLT.matrixL();
-
-    //if(!MTM_LLT.succeeded())
-    if(!((K*0).eval().nonZeros() == 0))
-    {
-      fprintf(stderr,"Error: lu_lagrange() failed to factor MTM\n");
-      return false;
-    }
-
-    // assemble LU decomposition of Q
-    Eigen::Matrix<int,Eigen::Dynamic,1> MI;
-    Eigen::Matrix<int,Eigen::Dynamic,1> MJ;
-    Eigen::Matrix<T,Eigen::Dynamic,1> MV;
-    igl::find(M,MI,MJ,MV);
-
-    Eigen::Matrix<int,Eigen::Dynamic,1> KI;
-    Eigen::Matrix<int,Eigen::Dynamic,1> KJ;
-    Eigen::Matrix<T,Eigen::Dynamic,1> KV;
-    igl::find(K,KI,KJ,KV);
-
-    Eigen::Matrix<int,Eigen::Dynamic,1> JI;
-    Eigen::Matrix<int,Eigen::Dynamic,1> JJ;
-    Eigen::Matrix<T,Eigen::Dynamic,1> JV;
-    igl::find(J,JI,JJ,JV);
-
-    int nnz = JV.size()  + MV.size() + KV.size();
-
-    Eigen::Matrix<int,Eigen::Dynamic,1> UI(nnz);
-    Eigen::Matrix<int,Eigen::Dynamic,1> UJ(nnz);
-    Eigen::Matrix<T,Eigen::Dynamic,1> UV(nnz);
-    UI << JJ,                        MI, (KJ.array() + n).matrix();
-    UJ << JI, (MJ.array() + n).matrix(), (KI.array() + n).matrix(); 
-    UV << JV,                        MV,                     KV*-1;
-    igl::sparse(UI,UJ,UV,U);
-
-    Eigen::Matrix<int,Eigen::Dynamic,1> LI(nnz);
-    Eigen::Matrix<int,Eigen::Dynamic,1> LJ(nnz);
-    Eigen::Matrix<T,Eigen::Dynamic,1> LV(nnz);
-    LI << JI, (MJ.array() + n).matrix(), (KI.array() + n).matrix();
-    LJ << JJ,                        MI, (KJ.array() + n).matrix(); 
-    LV << JV,                        MV,                        KV;
-    igl::sparse(LI,LJ,LV,L);
-  }
-
-  return true;
-}
+#ifdef IGL_HEADER_ONLY
+#  include "lu_lagrange.cpp"
+#endif
 
 #endif

mat_max.cpp

@@ -0,0 +1,38 @@
+#include "mat_max.h"
+
+template <typename T>
+IGL_INLINE void igl::mat_max(
+  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
+  const int dim,
+  Eigen::Matrix<T,Eigen::Dynamic,1> & Y,
+  Eigen::Matrix<int,Eigen::Dynamic,1> & I)
+{
+  assert(dim==1||dim==2);
+
+  // output size
+  int n = (dim==1?X.cols():X.rows());
+  // resize output
+  Y.resize(n);
+  I.resize(n);
+
+  // loop over dimension opposite of dim
+  for(int j = 0;j<n;j++)
+  {
+    int PHONY;
+    int i;
+    int m;
+    if(dim==1)
+    {
+      m = X.col(j).maxCoeff(&i,&PHONY);
+    }else
+    {
+      m = X.row(j).maxCoeff(&PHONY,&i);
+    }
+    Y(j) = m;
+    I(j) = i;
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

mat_max.h

@@ -1,5 +1,6 @@
 #ifndef IGL_MAT_MAX_H
 #define IGL_MAT_MAX_H
+#include "igl_inline.h"
 #include <Eigen/Dense>
 
 namespace igl
@@ -22,47 +23,15 @@ namespace igl
   //   I  vector the same size as Y containing the indices along dim of maximum
   //     entries
   template <typename T>
-  inline void mat_max(
+  IGL_INLINE void mat_max(
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
     const int dim,
     Eigen::Matrix<T,Eigen::Dynamic,1> & Y,
     Eigen::Matrix<int,Eigen::Dynamic,1> & I);
 }
 
-// Implementation
-
-template <typename T>
-inline void igl::mat_max(
-  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
-  const int dim,
-  Eigen::Matrix<T,Eigen::Dynamic,1> & Y,
-  Eigen::Matrix<int,Eigen::Dynamic,1> & I)
-{
-  assert(dim==1||dim==2);
-
-  // output size
-  int n = (dim==1?X.cols():X.rows());
-  // resize output
-  Y.resize(n);
-  I.resize(n);
-
-  // loop over dimension opposite of dim
-  for(int j = 0;j<n;j++)
-  {
-    int PHONY;
-    int i;
-    int m;
-    if(dim==1)
-    {
-      m = X.col(j).maxCoeff(&i,&PHONY);
-    }else
-    {
-      m = X.row(j).maxCoeff(&PHONY,&i);
-    }
-    Y(j) = m;
-    I(j) = i;
-  }
-}
-
+#ifdef IGL_HEADER_ONLY
+#  include "mat_max.cpp"
 #endif
 
+#endif

mat_min.cpp

@@ -0,0 +1,40 @@
+#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,
+  const int dim,
+  Eigen::Matrix<T,Eigen::Dynamic,1> & Y,
+  Eigen::Matrix<int,Eigen::Dynamic,1> & I)
+{
+  assert(dim==1||dim==2);
+
+  // output size
+  int n = (dim==1?X.cols():X.rows());
+  // resize output
+  Y.resize(n);
+  I.resize(n);
+
+  // loop over dimension opposite of dim
+  for(int j = 0;j<n;j++)
+  {
+    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).minCoeff(&i,&PHONY);
+    }else
+    {
+      m = X.row(j).minCoeff(&PHONY,&i);
+    }
+    Y(j) = m;
+    I(j) = i;
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

mat_min.h

@@ -1,5 +1,6 @@
 #ifndef IGL_MAT_MIN_H
 #define IGL_MAT_MIN_H
+#include "igl_inline.h"
 #include <Eigen/Dense>
 
 namespace igl
@@ -24,48 +25,15 @@ namespace igl
   //
   // See also: mat_max
   template <typename T>
-  inline void mat_min(
+  IGL_INLINE void mat_min(
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
     const int dim,
     Eigen::Matrix<T,Eigen::Dynamic,1> & Y,
     Eigen::Matrix<int,Eigen::Dynamic,1> & I);
 }
 
-// Implementation
-#include "verbose.h"
-
-template <typename T>
-inline void igl::mat_min(
-  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
-  const int dim,
-  Eigen::Matrix<T,Eigen::Dynamic,1> & Y,
-  Eigen::Matrix<int,Eigen::Dynamic,1> & I)
-{
-  assert(dim==1||dim==2);
-
-  // output size
-  int n = (dim==1?X.cols():X.rows());
-  // resize output
-  Y.resize(n);
-  I.resize(n);
-
-  // loop over dimension opposite of dim
-  for(int j = 0;j<n;j++)
-  {
-    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).minCoeff(&i,&PHONY);
-    }else
-    {
-      m = X.row(j).minCoeff(&PHONY,&i);
-    }
-    Y(j) = m;
-    I(j) = i;
-  }
-}
-
+#ifdef IGL_HEADER_ONLY
+#  include "mat_min.cpp"
 #endif
 
+#endif

max_size.cpp

@@ -0,0 +1,22 @@
+#include "max_size.h"
+
+
+template <typename T>
+IGL_INLINE int igl::max_size(const std::vector<T> & V)
+{
+  int max_size = -1;
+  for(
+    typename std::vector<T>::const_iterator iter = V.begin();
+    iter != V.end(); 
+    iter++)
+  {
+    int size = (int)iter->size();
+    max_size = (max_size > size ? max_size : size);
+  }
+  return max_size;
+}
+
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

max_size.h

@@ -1,5 +1,6 @@
 #ifndef IGL_MAX_SIZE_H
 #define IGL_MAX_SIZE_H
+#include "igl_inline.h"
 #include <vector>
 
 namespace igl
@@ -11,23 +12,11 @@ namespace igl
   //   V  vector of list types T
   // Returns max .size() found in V, returns -1 if V is empty
   template <typename T>
-  inline int max_size(const std::vector<T> & V);
+  IGL_INLINE int max_size(const std::vector<T> & V);
 }
 
-// Implementation 
+#ifdef IGL_HEADER_ONLY
+#  include "max_size.cpp"
+#endif
 
-template <typename T>
-inline int igl::max_size(const std::vector<T> & V)
-{
-  int max_size = -1;
-  for(
-    typename std::vector<T>::const_iterator iter = V.begin();
-    iter != V.end(); 
-    iter++)
-  {
-    int size = (int)iter->size();
-    max_size = (max_size > size ? max_size : size);
-  }
-  return max_size;
-}
 #endif

min_quad_dense.cpp

@@ -0,0 +1,84 @@
+#include "min_quad_dense.h"
+
+#include <Eigen/Core>
+#include <Eigen/LU>
+#include "EPS.h"
+
+template <typename T>
+IGL_INLINE void igl::min_quad_dense_precompute(
+  const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& A,
+  const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& Aeq,    
+  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& S)
+{
+  typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> Mat;
+        // This threshold seems to matter a lot but I'm not sure how to
+        // set it
+  const T treshold = igl::FLOAT_EPS;
+  //const T treshold = igl::DOUBLE_EPS;
+
+  const int n = A.rows();
+  assert(A.cols() == n);
+  const int m = Aeq.rows();
+  assert(Aeq.cols() == n);
+
+  // Lagrange multipliers method:
+  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> LM(n + m, n + m);
+  LM.block(0, 0, n, n) = A;
+  LM.block(0, n, n, m) = Aeq.transpose();
+  LM.block(n, 0, m, n) = Aeq;
+  LM.block(n, n, m, m).setZero();
+
+//#define USE_LU_DECOMPOSITION
+#ifdef USE_LU_DECOMPOSITION    
+  // if LM is close to singular, use at your own risk :)
+  Mat LMpinv = LM.inverse();
+#else // use SVD
+  typedef Eigen::Matrix<T, Eigen::Dynamic, 1> Vec; 
+  Vec singValues;
+  Eigen::JacobiSVD<Mat> svd;
+  svd.compute(LM, Eigen::ComputeFullU | Eigen::ComputeFullV );
+  const Mat& u = svd.matrixU();
+  const Mat& v = svd.matrixV();
+  const Vec& singVals = svd.singularValues();
+
+  Vec pi_singVals(n + m);
+  int zeroed = 0;
+  for (int i=0; i<n + m; i++)
+  {
+    T sv = singVals(i, 0);
+    assert(sv >= 0);      
+               // printf("sv: %lg ? %lg\n",(double) sv,(double)treshold);
+    if (sv > treshold) pi_singVals(i, 0) = T(1) / sv;
+    else 
+    {
+      pi_singVals(i, 0) = T(0);
+      zeroed++;
+    }
+  }
+
+  printf("min_quad_dense_precompute: %i singular values zeroed (threshold = %e)\n", zeroed, treshold);
+  Eigen::DiagonalMatrix<T, Eigen::Dynamic> pi_diag(pi_singVals);
+
+  Mat LMpinv = v * pi_diag * u.transpose();
+#endif
+  S = LMpinv.block(0, 0, n, n + m);
+
+  //// debug:
+  //mlinit(&g_pEngine);
+  //
+  //mlsetmatrix(&g_pEngine, "A", A);
+  //mlsetmatrix(&g_pEngine, "Aeq", Aeq);
+  //mlsetmatrix(&g_pEngine, "LM", LM);
+  //mlsetmatrix(&g_pEngine, "u", u);
+  //mlsetmatrix(&g_pEngine, "v", v);
+  //MatrixXd svMat = singVals;
+  //mlsetmatrix(&g_pEngine, "singVals", svMat);
+  //mlsetmatrix(&g_pEngine, "LMpinv", LMpinv);
+  //mlsetmatrix(&g_pEngine, "S", S);
+
+  //int hu = 1;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

min_quad_dense.h

@@ -1,15 +1,13 @@
 #ifndef IGL_MIN_QUAD_DENSE_H
 #define IGL_MIN_QUAD_DENSE_H
+#include "igl_inline.h"
 
-#include <Eigen/Core>
 #include <Eigen/Dense>
-#include <Eigen/LU>
 
 //// debug
 //#include <matlabinterface.h>
 //Engine *g_pEngine;
 
-#include <EPS.h>
 
 namespace igl
 {
@@ -27,86 +25,14 @@ namespace igl
   //   S  n by (n + m) "solve" matrix, such that S*[B', Beq'] is a solution
   // Returns true on success, false on error
   template <typename T>
-  inline void min_quad_dense_precompute(
+  IGL_INLINE void min_quad_dense_precompute(
     const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& A,
     const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& Aeq,    
     Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& S);
 }
 
-// Implementation
-template <typename T>
-inline void igl::min_quad_dense_precompute(
-  const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& A,
-  const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& Aeq,    
-  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& S)
-{
-  typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> Mat;
-        // This threshold seems to matter a lot but I'm not sure how to
-        // set it
-  const T treshold = igl::FLOAT_EPS;
-  //const T treshold = igl::DOUBLE_EPS;
-
-  const int n = A.rows();
-  assert(A.cols() == n);
-  const int m = Aeq.rows();
-  assert(Aeq.cols() == n);
-
-  // Lagrange multipliers method:
-  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> LM(n + m, n + m);
-  LM.block(0, 0, n, n) = A;
-  LM.block(0, n, n, m) = Aeq.transpose();
-  LM.block(n, 0, m, n) = Aeq;
-  LM.block(n, n, m, m).setZero();
-
-//#define USE_LU_DECOMPOSITION
-#ifdef USE_LU_DECOMPOSITION    
-  // if LM is close to singular, use at your own risk :)
-  Mat LMpinv = LM.inverse();
-#else // use SVD
-  typedef Eigen::Matrix<T, Eigen::Dynamic, 1> Vec; 
-  Vec singValues;
-  Eigen::JacobiSVD<Mat> svd;
-  svd.compute(LM, Eigen::ComputeFullU | Eigen::ComputeFullV );
-  const Mat& u = svd.matrixU();
-  const Mat& v = svd.matrixV();
-  const Vec& singVals = svd.singularValues();
-
-  Vec pi_singVals(n + m);
-  int zeroed = 0;
-  for (int i=0; i<n + m; i++)
-  {
-    T sv = singVals(i, 0);
-    assert(sv >= 0);      
-               // printf("sv: %lg ? %lg\n",(double) sv,(double)treshold);
-    if (sv > treshold) pi_singVals(i, 0) = T(1) / sv;
-    else 
-    {
-      pi_singVals(i, 0) = T(0);
-      zeroed++;
-    }
-  }
-
-  printf("min_quad_dense_precompute: %i singular values zeroed (threshold = %e)\n", zeroed, treshold);
-  Eigen::DiagonalMatrix<T, Eigen::Dynamic> pi_diag(pi_singVals);
-
-  Mat LMpinv = v * pi_diag * u.transpose();
+#ifdef IGL_HEADER_ONLY
+#  include "min_quad_dense.cpp"
 #endif
-  S = LMpinv.block(0, 0, n, n + m);
-
-  //// debug:
-  //mlinit(&g_pEngine);
-  //
-  //mlsetmatrix(&g_pEngine, "A", A);
-  //mlsetmatrix(&g_pEngine, "Aeq", Aeq);
-  //mlsetmatrix(&g_pEngine, "LM", LM);
-  //mlsetmatrix(&g_pEngine, "u", u);
-  //mlsetmatrix(&g_pEngine, "v", v);
-  //MatrixXd svMat = singVals;
-  //mlsetmatrix(&g_pEngine, "singVals", svMat);
-  //mlsetmatrix(&g_pEngine, "LMpinv", LMpinv);
-  //mlsetmatrix(&g_pEngine, "S", S);
-
-  //int hu = 1;
-}
 
 #endif

min_quad_with_fixed.cpp

@@ -0,0 +1,251 @@
+#include "min_quad_with_fixed.h"
+
+#include <Eigen/SparseExtra>
+#include <cassert>
+#include <cstdio>
+#include <iostream>
+
+#include "slice.h"
+#include "is_symmetric.h"
+#include "find.h"
+#include "sparse.h"
+#include "repmat.h"
+#include "lu_lagrange.h"
+#include "full.h"
+
+template <typename T>
+IGL_INLINE bool igl::min_quad_with_fixed_precompute(
+  const Eigen::SparseMatrix<T>& A,
+  const Eigen::Matrix<int,Eigen::Dynamic,1> & known,
+  const Eigen::SparseMatrix<T>& Aeq,
+  const bool pd,
+  igl::min_quad_with_fixed_data<T> & data
+  )
+{
+  // number of rows
+  int n = A.rows();
+  // cache problem size
+  data.n = n;
+
+  int neq = Aeq.rows();
+  // default is to have 0 linear equality constraints
+  if(Aeq.size() != 0)
+  {
+    assert(n == Aeq.cols());
+  }
+
+  assert(A.rows() == n);
+  assert(A.cols() == n);
+
+  // number of known rows
+  int kr = known.size();
+
+  assert(kr == 0 || known.minCoeff() >= 0);
+  assert(kr == 0 || known.maxCoeff() < n);
+  assert(neq <= n);
+
+  // cache known
+  data.known = known;
+  // get list of unknown indices
+  data.unknown.resize(n-kr);
+  std::vector<bool> unknown_mask;
+  unknown_mask.resize(n,true);
+  for(int i = 0;i<kr;i++)
+  {
+    unknown_mask[known(i)] = false;
+  }
+  int u = 0;
+  for(int i = 0;i<n;i++)
+  {
+    if(unknown_mask[i])
+    {
+      data.unknown(u) = i;
+      u++;
+    }
+  }
+  // get list of lagrange multiplier indices
+  data.lagrange.resize(neq);
+  for(int i = 0;i<neq;i++)
+  {
+    data.lagrange(i) = n + i;
+  }
+  // cache unknown followed by lagrange indices
+  data.unknown_lagrange.resize(data.unknown.size()+data.lagrange.size());
+  data.unknown_lagrange << data.unknown, data.lagrange;
+
+  Eigen::SparseMatrix<T> Auu;
+  igl::slice(A,data.unknown,data.unknown,Auu);
+
+  // determine if A(unknown,unknown) is symmetric and/or positive definite
+  data.Auu_sym = igl::is_symmetric(Auu);
+  // Positive definiteness is *not* determined, rather it is given as a
+  // parameter
+  data.Auu_pd = pd;
+
+  // Append lagrange multiplier quadratic terms
+  Eigen::SparseMatrix<T> new_A;
+  Eigen::Matrix<int,Eigen::Dynamic,1> AI;
+  Eigen::Matrix<int,Eigen::Dynamic,1> AJ;
+  Eigen::Matrix<T,Eigen::Dynamic,1> AV;
+  igl::find(A,AI,AJ,AV);
+  Eigen::Matrix<int,Eigen::Dynamic,1> AeqI(0,1);
+  Eigen::Matrix<int,Eigen::Dynamic,1> AeqJ(0,1);
+  Eigen::Matrix<T,Eigen::Dynamic,1> AeqV(0,1);
+  if(neq > 0)
+  {
+    igl::find(Aeq,AeqI,AeqJ,AeqV);
+  }
+  Eigen::Matrix<int,Eigen::Dynamic,1> new_AI(AV.size()+AeqV.size()*2);
+  Eigen::Matrix<int,Eigen::Dynamic,1> new_AJ(AV.size()+AeqV.size()*2);
+  Eigen::Matrix<T,Eigen::Dynamic,1>   new_AV(AV.size()+AeqV.size()*2);
+  new_AI << AI, (AeqI.array()+n).matrix(), AeqJ;
+  new_AJ << AJ, AeqJ, (AeqI.array()+n).matrix();
+  new_AV << AV, AeqV, AeqV;
+  igl::sparse(new_AI,new_AJ,new_AV,n+neq,n+neq,new_A);
+
+  // precompute RHS builders
+  if(kr > 0)
+  {
+    Eigen::SparseMatrix<T> Aulk;
+    igl::slice(new_A,data.unknown_lagrange,data.known,Aulk);
+    Eigen::SparseMatrix<T> Akul;
+    igl::slice(new_A,data.known,data.unknown_lagrange,Akul);
+
+    //// This doesn't work!!!
+    //data.preY = Aulk + Akul.transpose();
+    Eigen::SparseMatrix<T> AkulT = Akul.transpose();
+    data.preY = Aulk + AkulT;
+  }else
+  {
+    data.preY.resize(data.unknown_lagrange.size(),0);
+  }
+
+  // Create factorization
+  if(data.Auu_sym && data.Auu_pd)
+  {
+    data.sparse = true;
+    Eigen::SparseMatrix<T> Aequ(0,0);
+    if(neq>0)
+    {
+      Eigen::Matrix<int,Eigen::Dynamic,1> Aeq_rows;
+      Aeq_rows.resize(neq);
+      for(int i = 0;i<neq;i++)
+      {
+        Aeq_rows(i) = i;
+      }
+      igl::slice(Aeq,Aeq_rows,data.unknown,Aequ);
+    }
+    // Get transpose of Aequ
+    Eigen::SparseMatrix<T> AequT = Aequ.transpose();
+    // Compute LU decomposition
+    bool lu_success = igl::lu_lagrange(Auu,AequT,data.L,data.U);
+    if(!lu_success)
+    {
+      return false;
+    }
+  }else
+  {
+    Eigen::SparseMatrix<T> NA;
+    igl::slice(new_A,data.unknown_lagrange,data.unknown_lagrange,NA);
+    // Build LU decomposition of NA
+    data.sparse = false;
+    fprintf(stderr,
+      "Warning: min_quad_with_fixed_precompute() resorting to dense LU\n");
+    Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NAfull;
+    igl::full(NA,NAfull);
+    data.lu = 
+      Eigen::FullPivLU<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >(NAfull);
+    if(!data.lu.isInvertible())
+    {
+      fprintf(stderr,
+        "Error: min_quad_with_fixed_precompute() LU not invertible\n");
+      return false;
+    }
+  }
+  return true;
+}
+
+
+template <typename T>
+IGL_INLINE bool igl::min_quad_with_fixed_solve(
+  const igl::min_quad_with_fixed_data<T> & data,
+  const Eigen::Matrix<T,Eigen::Dynamic,1> & B,
+  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Y,
+  const Eigen::Matrix<T,Eigen::Dynamic,1> & Beq,
+  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Z)
+{
+  // number of known rows
+  int kr = data.known.size();
+  if(kr!=0)
+  {
+    assert(kr == Y.rows());
+  }
+  // number of columns to solve
+  int cols = Y.cols();
+
+  // number of lagrange multipliers aka linear equality constraints
+  int neq = data.lagrange.size();
+
+  // append lagrange multiplier rhs's
+  Eigen::Matrix<T,Eigen::Dynamic,1> BBeq(B.size() + Beq.size());
+  BBeq << B, (Beq*-2.0);
+
+  // Build right hand side
+  Eigen::Matrix<T,Eigen::Dynamic,1> BBequl;
+  igl::slice(BBeq,data.unknown_lagrange,BBequl);
+  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> BBequlcols;
+  igl::repmat(BBequl,1,cols,BBequlcols);
+  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NB;
+  if(kr == 0)
+  {
+    NB = BBequlcols;
+  }else
+  {
+    NB = data.preY * Y + BBequlcols;
+  }
+
+  // resize output
+  Z.resize(data.n,cols);
+
+  // Set known values
+  for(int i = 0;i < kr;i++)
+  {
+    for(int j = 0;j < cols;j++)
+    {
+      Z(data.known(i),j) = Y(i,j);
+    }
+  }
+
+  //std::cout<<"NB=["<<std::endl<<NB<<std::endl<<"];"<<std::endl;
+
+  if(data.sparse)
+  {
+    //std::cout<<"data.LIJV=["<<std::endl;print_ijv(data.L,1);std::cout<<std::endl<<"];"<<
+    //  std::endl<<"data.L=sparse(data.LIJV(:,1),data.LIJV(:,2),data.LIJV(:,3),"<<
+    //  data.L.rows()<<","<<data.L.cols()<<");"<<std::endl;
+    //std::cout<<"data.UIJV=["<<std::endl;print_ijv(data.U,1);std::cout<<std::endl<<"];"<<
+    //  std::endl<<"data.U=sparse(data.UIJV(:,1),data.UIJV(:,2),data.UIJV(:,3),"<<
+    //  data.U.rows()<<","<<data.U.cols()<<");"<<std::endl;
+    data.L.template triangularView<Eigen::Lower>().solveInPlace(NB);
+    data.U.template triangularView<Eigen::Upper>().solveInPlace(NB);
+  }else
+  {
+    NB = data.lu.solve(NB);
+  }
+  // Now NB contains sol/-0.5
+  NB *= -0.5;
+  // Now NB contains solution
+  // Place solution in Z
+  for(int i = 0;i<(NB.rows()-neq);i++)
+  {
+    for(int j = 0;j<NB.cols();j++)
+    {
+      Z(data.unknown_lagrange(i),j) = NB(i,j);
+    }
+  }
+  return true;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

min_quad_with_fixed.h

@@ -1,5 +1,6 @@
 #ifndef IGL_MIN_QUAD_WITH_FIXED_H
 #define IGL_MIN_QUAD_WITH_FIXED_H
+#include "igl_inline.h"
 
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Core>
@@ -31,7 +32,7 @@ namespace igl
   //     using min_quad_with_fixed_solve
   // Returns true on success, false on error
   template <typename T>
-  inline bool min_quad_with_fixed_precompute(
+  IGL_INLINE bool min_quad_with_fixed_precompute(
     const Eigen::SparseMatrix<T>& A,
     const Eigen::Matrix<int,Eigen::Dynamic,1> & known,
     const Eigen::SparseMatrix<T>& Aeq,
@@ -46,7 +47,7 @@ namespace igl
   //   Z  n by cols solution
   // Returns true on success, false on error
   template <typename T>
-  inline bool min_quad_with_fixed_solve(
+  IGL_INLINE bool min_quad_with_fixed_solve(
     const min_quad_with_fixed_data<T> & data,
     const Eigen::Matrix<T,Eigen::Dynamic,1> & B,
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Y,
@@ -54,20 +55,6 @@ namespace igl
     Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Z);
 }
 
-// Implementation
-#include <Eigen/SparseExtra>
-#include <cassert>
-#include <cstdio>
-#include <iostream>
-
-#include "slice.h"
-#include "is_symmetric.h"
-#include "find.h"
-#include "sparse.h"
-#include "repmat.h"
-#include "lu_lagrange.h"
-#include "full.h"
-
 template <typename T>
 struct igl::min_quad_with_fixed_data
 {
@@ -97,236 +84,8 @@ struct igl::min_quad_with_fixed_data
   Eigen::FullPivLU<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > lu;
 };
 
-template <typename T>
-inline bool igl::min_quad_with_fixed_precompute(
-  const Eigen::SparseMatrix<T>& A,
-  const Eigen::Matrix<int,Eigen::Dynamic,1> & known,
-  const Eigen::SparseMatrix<T>& Aeq,
-  const bool pd,
-  igl::min_quad_with_fixed_data<T> & data
-  )
-{
-  // number of rows
-  int n = A.rows();
-  // cache problem size
-  data.n = n;
-
-  int neq = Aeq.rows();
-  // default is to have 0 linear equality constraints
-  if(Aeq.size() != 0)
-  {
-    assert(n == Aeq.cols());
-  }
-
-  assert(A.rows() == n);
-  assert(A.cols() == n);
-
-  // number of known rows
-  int kr = known.size();
-
-  assert(kr == 0 || known.minCoeff() >= 0);
-  assert(kr == 0 || known.maxCoeff() < n);
-  assert(neq <= n);
-
-  // cache known
-  data.known = known;
-  // get list of unknown indices
-  data.unknown.resize(n-kr);
-  std::vector<bool> unknown_mask;
-  unknown_mask.resize(n,true);
-  for(int i = 0;i<kr;i++)
-  {
-    unknown_mask[known(i)] = false;
-  }
-  int u = 0;
-  for(int i = 0;i<n;i++)
-  {
-    if(unknown_mask[i])
-    {
-      data.unknown(u) = i;
-      u++;
-    }
-  }
-  // get list of lagrange multiplier indices
-  data.lagrange.resize(neq);
-  for(int i = 0;i<neq;i++)
-  {
-    data.lagrange(i) = n + i;
-  }
-  // cache unknown followed by lagrange indices
-  data.unknown_lagrange.resize(data.unknown.size()+data.lagrange.size());
-  data.unknown_lagrange << data.unknown, data.lagrange;
-
-  Eigen::SparseMatrix<T> Auu;
-  igl::slice(A,data.unknown,data.unknown,Auu);
-
-  // determine if A(unknown,unknown) is symmetric and/or positive definite
-  data.Auu_sym = igl::is_symmetric(Auu);
-  // Positive definiteness is *not* determined, rather it is given as a
-  // parameter
-  data.Auu_pd = pd;
-
-  // Append lagrange multiplier quadratic terms
-  Eigen::SparseMatrix<T> new_A;
-  Eigen::Matrix<int,Eigen::Dynamic,1> AI;
-  Eigen::Matrix<int,Eigen::Dynamic,1> AJ;
-  Eigen::Matrix<T,Eigen::Dynamic,1> AV;
-  igl::find(A,AI,AJ,AV);
-  Eigen::Matrix<int,Eigen::Dynamic,1> AeqI(0,1);
-  Eigen::Matrix<int,Eigen::Dynamic,1> AeqJ(0,1);
-  Eigen::Matrix<T,Eigen::Dynamic,1> AeqV(0,1);
-  if(neq > 0)
-  {
-    igl::find(Aeq,AeqI,AeqJ,AeqV);
-  }
-  Eigen::Matrix<int,Eigen::Dynamic,1> new_AI(AV.size()+AeqV.size()*2);
-  Eigen::Matrix<int,Eigen::Dynamic,1> new_AJ(AV.size()+AeqV.size()*2);
-  Eigen::Matrix<T,Eigen::Dynamic,1>   new_AV(AV.size()+AeqV.size()*2);
-  new_AI << AI, (AeqI.array()+n).matrix(), AeqJ;
-  new_AJ << AJ, AeqJ, (AeqI.array()+n).matrix();
-  new_AV << AV, AeqV, AeqV;
-  igl::sparse(new_AI,new_AJ,new_AV,n+neq,n+neq,new_A);
-
-  // precompute RHS builders
-  if(kr > 0)
-  {
-    Eigen::SparseMatrix<T> Aulk;
-    igl::slice(new_A,data.unknown_lagrange,data.known,Aulk);
-    Eigen::SparseMatrix<T> Akul;
-    igl::slice(new_A,data.known,data.unknown_lagrange,Akul);
-
-    //// This doesn't work!!!
-    //data.preY = Aulk + Akul.transpose();
-    Eigen::SparseMatrix<T> AkulT = Akul.transpose();
-    data.preY = Aulk + AkulT;
-  }else
-  {
-    data.preY.resize(data.unknown_lagrange.size(),0);
-  }
-
-  // Create factorization
-  if(data.Auu_sym && data.Auu_pd)
-  {
-    data.sparse = true;
-    Eigen::SparseMatrix<T> Aequ(0,0);
-    if(neq>0)
-    {
-      Eigen::Matrix<int,Eigen::Dynamic,1> Aeq_rows;
-      Aeq_rows.resize(neq);
-      for(int i = 0;i<neq;i++)
-      {
-        Aeq_rows(i) = i;
-      }
-      igl::slice(Aeq,Aeq_rows,data.unknown,Aequ);
-    }
-    // Get transpose of Aequ
-    Eigen::SparseMatrix<T> AequT = Aequ.transpose();
-    // Compute LU decomposition
-    bool lu_success = igl::lu_lagrange(Auu,AequT,data.L,data.U);
-    if(!lu_success)
-    {
-      return false;
-    }
-  }else
-  {
-    Eigen::SparseMatrix<T> NA;
-    igl::slice(new_A,data.unknown_lagrange,data.unknown_lagrange,NA);
-    // Build LU decomposition of NA
-    data.sparse = false;
-    fprintf(stderr,
-      "Warning: min_quad_with_fixed_precompute() resorting to dense LU\n");
-    Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NAfull;
-    igl::full(NA,NAfull);
-    data.lu = 
-      Eigen::FullPivLU<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >(NAfull);
-    if(!data.lu.isInvertible())
-    {
-      fprintf(stderr,
-        "Error: min_quad_with_fixed_precompute() LU not invertible\n");
-      return false;
-    }
-  }
-  return true;
-}
-
-
-template <typename T>
-inline bool igl::min_quad_with_fixed_solve(
-  const igl::min_quad_with_fixed_data<T> & data,
-  const Eigen::Matrix<T,Eigen::Dynamic,1> & B,
-  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Y,
-  const Eigen::Matrix<T,Eigen::Dynamic,1> & Beq,
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Z)
-{
-  // number of known rows
-  int kr = data.known.size();
-  if(kr!=0)
-  {
-    assert(kr == Y.rows());
-  }
-  // number of columns to solve
-  int cols = Y.cols();
-
-  // number of lagrange multipliers aka linear equality constraints
-  int neq = data.lagrange.size();
-
-  // append lagrange multiplier rhs's
-  Eigen::Matrix<T,Eigen::Dynamic,1> BBeq(B.size() + Beq.size());
-  BBeq << B, (Beq*-2.0);
-
-  // Build right hand side
-  Eigen::Matrix<T,Eigen::Dynamic,1> BBequl;
-  igl::slice(BBeq,data.unknown_lagrange,BBequl);
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> BBequlcols;
-  igl::repmat(BBequl,1,cols,BBequlcols);
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NB;
-  if(kr == 0)
-  {
-    NB = BBequlcols;
-  }else
-  {
-    NB = data.preY * Y + BBequlcols;
-  }
-
-  // resize output
-  Z.resize(data.n,cols);
-
-  // Set known values
-  for(int i = 0;i < kr;i++)
-  {
-    for(int j = 0;j < cols;j++)
-    {
-      Z(data.known(i),j) = Y(i,j);
-    }
-  }
-
-  //std::cout<<"NB=["<<std::endl<<NB<<std::endl<<"];"<<std::endl;
+#ifdef IGL_HEADER_ONLY
+#  include "min_quad_with_fixed.cpp"
+#endif
 
-  if(data.sparse)
-  {
-    //std::cout<<"data.LIJV=["<<std::endl;print_ijv(data.L,1);std::cout<<std::endl<<"];"<<
-    //  std::endl<<"data.L=sparse(data.LIJV(:,1),data.LIJV(:,2),data.LIJV(:,3),"<<
-    //  data.L.rows()<<","<<data.L.cols()<<");"<<std::endl;
-    //std::cout<<"data.UIJV=["<<std::endl;print_ijv(data.U,1);std::cout<<std::endl<<"];"<<
-    //  std::endl<<"data.U=sparse(data.UIJV(:,1),data.UIJV(:,2),data.UIJV(:,3),"<<
-    //  data.U.rows()<<","<<data.U.cols()<<");"<<std::endl;
-    data.L.template triangularView<Eigen::Lower>().solveInPlace(NB);
-    data.U.template triangularView<Eigen::Upper>().solveInPlace(NB);
-  }else
-  {
-    NB = data.lu.solve(NB);
-  }
-  // Now NB contains sol/-0.5
-  NB *= -0.5;
-  // Now NB contains solution
-  // Place solution in Z
-  for(int i = 0;i<(NB.rows()-neq);i++)
-  {
-    for(int j = 0;j<NB.cols();j++)
-    {
-      Z(data.unknown_lagrange(i),j) = NB(i,j);
-    }
-  }
-  return true;
-}
 #endif

min_size.cpp

@@ -0,0 +1,28 @@
+#include "min_size.h"
+
+template <typename T>
+IGL_INLINE int igl::min_size(const std::vector<T> & V)
+{
+  int min_size = -1;
+  for(
+    typename std::vector<T>::const_iterator iter = V.begin();
+    iter != V.end(); 
+    iter++)
+  {
+    int size = (int)iter->size();
+    // have to handle base case
+    if(min_size == -1)
+    {
+      min_size = size;
+    }else{
+      min_size = (min_size < size ? min_size : size);
+    }
+  }
+  return min_size;
+}
+
+
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif