Merge branch 'master' of https://github.com/libigl/libigl

Conflicts:
	include/igl/embree/EmbreeIntersector.h
	include/igl/triangle/triangulate.cpp
	include/igl/viewer/Viewer.cpp

Former-commit-id: e24f8a4848385fbe08e76a08111c20b9649aaa79

CMakeLists.txt

@@ -0,0 +1,238 @@
+cmake_minimum_required(VERSION 2.6)
+project(libigl)
+
+SET(CMAKE_SKIP_RULE_DEPENDENCY 1)
+
+# SET(CMAKE_RUNTIME_OUTPUT_DIRECTORY ../../lib)
+SET(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/tutorial/cmake)
+find_package(EIGEN REQUIRED)
+
+add_definitions(-DIGL_STATIC_LIBRARY)
+
+## Check for GLFW
+find_package(GLFW QUIET)
+if (GLFW_FOUND)
+  include_directories( ${GLFW_INCLUDE_DIR})
+endif(GLFW_FOUND)
+
+find_package(OpenGL QUIET)
+
+## Check for Anttweakbar
+find_package(ANTTWEAKBAR QUIET)
+if (ANTTWEAKBAR_FOUND)
+  include_directories( ${ANT_TWEAK_BAR_INCLUDE_DIR})
+endif(ANTTWEAKBAR_FOUND)
+
+
+## Check for CoMiSo, if not available skip the examples that depends on it
+find_package(LIBCOMISO QUIET)
+if (LIBCOMISO_FOUND)
+  include_directories( ${LIBCOMISO_INCLUDE_DIRS})
+endif(LIBCOMISO_FOUND)
+
+## Check for MATLAB, if not available skip the examples that depends on it
+find_package(MATLAB QUIET)
+if (MATLAB_FOUND)
+  include_directories( ${MATLAB_INCLUDE_DIR})
+endif(MATLAB_FOUND)
+
+## Check for EMBREE, if not available skip the examples that depends on it
+find_package(EMBREE QUIET)
+if (EMBREE_FOUND)
+  include_directories( ${EMBREE_INCLUDE_DIR})
+endif(EMBREE_FOUND)
+
+## Check for CGAL, if not available skip the examples that depends on it
+find_package(CGAL QUIET)
+
+## Check for mosek
+find_package(MOSEK QUIET)
+if(NOT MOSEK_FOUND)
+  add_definitions(-DIGL_NO_MOSEK)
+endif(NOT MOSEK_FOUND)
+
+## Check for CORK
+find_package(CORK QUIET)
+if (NOT CORK_FOUND)
+  add_definitions(-DIGL_NO_CORK)
+else(NOT CORK_FOUND)
+  include_directories( ${CORK_INCLUDE_DIR})
+endif(NOT CORK_FOUND)
+
+## Check for LIM
+find_package(LIM QUIET)
+if(LIM_FOUND)
+  include_directories( ${LIM_INCLUDE_DIR})
+endif(LIM_FOUND)
+
+## Check for SVD3X3
+find_package(SVD3X3 QUIET)
+if(SVD3X3_FOUND)
+  include_directories( ${SVD3X3_INCLUDE_DIR})
+endif(SVD3X3_FOUND)
+
+## Check for TETGEN
+find_package(TETGEN QUIET)
+if(TETGEN_FOUND)
+  include_directories( ${TETGEN_INCLUDE_DIR})
+endif(TETGEN_FOUND)
+
+## Check for TRIANGLE
+find_package(TRIANGLE QUIET)
+if(TRIANGLE_FOUND)
+  include_directories( ${TRIANGLE_INCLUDE_DIR})
+endif(TRIANGLE_FOUND)
+
+## Check for TINYXML2
+find_package(TINYXML2 QUIET)
+if(TINYXML2_FOUND)
+  include_directories( ${TINYXML2_INCLUDE_DIR})
+endif(TINYXML2_FOUND)
+
+## Check for COMISO
+find_package(COMISO QUIET)
+if(COMISO_FOUND)
+  include_directories( ${COMISO_INCLUDE_DIR})
+endif(COMISO_FOUND)
+
+## Use openMP if available
+find_package(OpenMP)
+if (OPENMP_FOUND AND NOT WIN32)
+  set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+  set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+endif()
+
+#### Libigl requires a modern C++ compiler that supports c++11
+SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+
+#### Compile the core library that depends only on EIGEN ####
+include_directories( ${EIGEN_INCLUDE_DIR})
+include_directories( ${PROJECT_SOURCE_DIR}/include/)
+
+file(GLOB SOURCES
+  "${PROJECT_SOURCE_DIR}/include/igl/*.cpp"
+)
+
+### HACKS TO CLEAN
+list(REMOVE_ITEM SOURCES ${ui} ${PROJECT_SOURCE_DIR}/include/igl/cocoa_key_to_anttweakbar_key.cpp)
+list(REMOVE_ITEM SOURCES ${ui} ${PROJECT_SOURCE_DIR}/include/igl/ReAntTweakBar.cpp)
+
+add_library(igl STATIC ${SOURCES})
+
+#### Compile the BBW part
+
+file(GLOB SOURCES_BBW
+  "${PROJECT_SOURCE_DIR}/include/igl/bbw/*.cpp"
+)
+
+add_library(iglbbw STATIC ${SOURCES_BBW})
+
+#### Compile the mosek part (untested)
+if (MOSEK_FOUND)
+  file(GLOB SOURCES_MOSEK
+    "${PROJECT_SOURCE_DIR}/include/igl/mosek/*.cpp"
+  )
+
+add_library(iglmosek STATIC ${SOURCES_MOSEK})
+endif (MOSEK_FOUND)
+
+#### Compile the cgal part
+if (CGAL_FOUND)
+  file(GLOB SOURCES_CGAL
+    "${PROJECT_SOURCE_DIR}/include/igl/cgal/*.cpp"
+  )
+
+add_library(iglcgal STATIC ${SOURCES_CGAL})
+
+#### Compile the boolean part
+if (NOT CORK_FOUND)
+  add_definitions(-DIGL_NO_CORK)
+endif(NOT CORK_FOUND)
+
+file(GLOB SOURCES_BOOLEAN
+  "${PROJECT_SOURCE_DIR}/include/igl/boolean/*.cpp"
+)
+
+add_library(iglboolean STATIC ${SOURCES_BOOLEAN})
+endif (CGAL_FOUND)
+
+#### Compile the embree part
+if (EMBREE_FOUND)
+  file(GLOB SOURCES_EMBREE
+    "${PROJECT_SOURCE_DIR}/include/igl/embree/*.cpp"
+  )
+
+  add_library(iglembree STATIC ${SOURCES_EMBREE})
+endif (EMBREE_FOUND)
+
+#### Compile the lim part
+if (LIM_FOUND)
+  file(GLOB SOURCES_LIM
+    "${PROJECT_SOURCE_DIR}/include/igl/lim/*.cpp"
+  )
+
+  add_library(igllim STATIC ${SOURCES_LIM})
+endif (LIM_FOUND)
+
+#### Compile the matlab part
+if (MATLAB_FOUND)
+  file(GLOB SOURCES_MATLAB
+    "${PROJECT_SOURCE_DIR}/include/igl/matlab/*.cpp"
+  )
+
+  add_library(iglmatlab STATIC ${SOURCES_MATLAB})
+endif (MATLAB_FOUND)
+
+#### Compile the svd3x3 part
+if (SVD3X3_FOUND)
+  file(GLOB SOURCES_SVD3X3
+    "${PROJECT_SOURCE_DIR}/include/igl/svd3x3/*.cpp"
+  )
+
+  add_library(iglsvd3x3 STATIC ${SOURCES_SVD3X3})
+endif (SVD3X3_FOUND)
+
+#### Compile the tetgen part
+if (TETGEN_FOUND)
+  file(GLOB SOURCES_TETGEN
+    "${PROJECT_SOURCE_DIR}/include/igl/tetgen/*.cpp"
+  )
+
+  add_library(igltetgen STATIC ${SOURCES_TETGEN})
+endif (TETGEN_FOUND)
+
+#### Compile the triangle part
+if (TRIANGLE_FOUND)
+  file(GLOB SOURCES_TRIANGLE
+    "${PROJECT_SOURCE_DIR}/include/igl/triangle/*.cpp"
+  )
+
+  add_library(igltriangle STATIC ${SOURCES_TRIANGLE})
+endif (TRIANGLE_FOUND)
+
+#### Compile the xml part
+if (TINYXML2_FOUND)
+  file(GLOB SOURCES_XML
+    "${PROJECT_SOURCE_DIR}/include/igl/xml/*.cpp"
+  )
+
+  add_library(iglxml STATIC ${SOURCES_XML})
+endif (TINYXML2_FOUND)
+
+#### Compile the xml part
+if (LIBCOMISO_FOUND)
+  file(GLOB SOURCES_COMISO
+    "${PROJECT_SOURCE_DIR}/include/igl/comiso/*.cpp"
+  )
+
+  add_library(iglcomiso STATIC ${SOURCES_COMISO})
+endif (LIBCOMISO_FOUND)
+
+#### Compile the viewer
+if (GLFW_FOUND AND ANTTWEAKBAR_FOUND)
+  file(GLOB SOURCES_XML
+    "${PROJECT_SOURCE_DIR}/include/igl/viewer/*.cpp"
+  )
+
+  add_library(iglviewer STATIC ${SOURCES_XML})
+endif (GLFW_FOUND AND ANTTWEAKBAR_FOUND)

README.md

@@ -143,6 +143,27 @@ BibTeX entry:
 }
 ```
 
+## Projects/Universities using libigl
+
+ - [Spine by Esoteric Software](http://esotericsoftware.com/) is an animation tool dedicated to 2D characters.
+ - Columbia University, [Columbia Computer Graphics Group](http://www.cs.columbia.edu/cg/), USA
+ - [Cornell University](http://www.graphics.cornell.edu/), USA
+ - EPF Lausanne, [Computer Graphics and Geometry Laboratory](http://lgg.epfl.ch/people.php), Switzerland
+ - ETH Zurich, [Interactive Geometry Lab](http://igl.ethz.ch/) and [Advanced Technologies Lab](http://ait.inf.ethz.ch/), Swizterland
+ - George Mason University, [CraGL](http://cs.gmu.edu/~ygingold/), USA
+ - [Hong Kong University of Science and Technology](http://www.ust.hk/), USA
+ - [National Institute of Informatics](http://www.nii.ac.jp/en/), Japan 
+ - New York University, [Media Research Lab](http://mrl.nyu.edu/), USA
+ - NYUPoly, [Game Innovation Lab](http://game.engineering.nyu.edu/), USA
+ - [Telecom ParisTech](http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html), Paris, France
+ - [TU Delft](http://www.tudelft.nl/en/), Netherlands
+ - [Universidade Federal de Santa Catarina](http://mtm.ufsc.br/~leo/), Brazil
+ - [Università della Svizzera Italiana](http://www.usi.ch/en), Switzerland
+ - [University College London](http://vecg.cs.ucl.ac.uk/), England
+ - [University of Cambridge](http://www.cam.ac.uk/), England
+ - [University of Pennsylvania](http://cg.cis.upenn.edu/), USA
+
+
 ## Contact
 
 Libigl is a group endeavor led by [Alec

RELEASE_HISTORY.txt

@@ -1,3 +1,4 @@
+1.1.6  Major boolean robustness fix, drop CGAL dependency for AABB/distances
 1.1.5  Bug fix in booleans
 1.1.4  Edge collapsing and linear program solving
 1.1.3  Bug fixes in active set and boundary_conditions

VERSION.txt

@@ -3,4 +3,4 @@
 # Anyone may increment Minor to indicate a small change.
 # Major indicates a large change or large number of changes (upload to website)
 # World indicates a substantial change or release
-1.1.5
+1.1.6

build/Makefile

@@ -9,7 +9,6 @@ $(info Hello, $(IGL_USERNAME)!)
 
 # optimized default settings
 all: LFLAGS +=
-OPTFLAGS+=-O3 -DNDEBUG $(OPENMP)
 #debug: OPTFLAGS= -g -Wall -Werror
 debug: OPTFLAGS= -g -Wall
 debug: DEBUG=debug

build/Makefile.conf

@@ -47,7 +47,7 @@ ifeq ($(IGL_USERNAME),whitinge)
 endif
 
 ifeq ($(IGL_USERNAME),ajx)
-	CFLAGS += -DIGL_STATIC_LIBRARY
+        LIBIGL_USE_STATIC_LIBRARY=1
 	MOSEKPLATFORM=osx64x86
 	MOSEKVERSION=7
 	IGL_WITH_VIEWER=1
@@ -85,6 +85,10 @@ ifeq ($(IGL_USERNAME),alecjaco)
 	OPENMP=-fopenmp
 endif
 
+ifeq ($(IGL_USERNAME),jacobson)
+	CFLAGS+=
+endif
+
 ifeq ($(IGL_USERNAME),sorkineo)
 	MOSEKPLATFORM=osx64x86
 	IGL_WITH_TETGEN=1
@@ -116,15 +120,15 @@ endif
 
 ifeq ($(IGL_USERNAME),olkido)
     IGL_WITH_MATLAB=1
-		IGL_WITH_COMISO=1
-		IGL_WITH_XML=1
-		COMISO=/Users/olkido/Documents/igl/MIQ/src/CoMISo
+	IGL_WITH_COMISO=1
+	IGL_WITH_XML=1
+	COMISO=/Users/olkido/Documents/igl/MIQ/src/CoMISo
 	  IGL_WITH_XML=1
     AFLAGS= -m64
     IGL_WITH_EMBREE=1
     IGL_WITH_PNG=1
-		IGL_WITH_VIEWER=1
-		MATLAB=/Applications/MATLAB_R2014b.app/
+	IGL_WITH_VIEWER=1
+	MATLAB=/Applications/MATLAB_R2014b.app/
 endif
 
 ifeq ($(IGL_USERNAME),chrsch)
@@ -209,7 +213,6 @@ endif
 ifndef LIBIGL
 	LIBIGL=$(THIS_DIR)/../
 endif
-LIBIGL_LIB=-L$(LIBIGL)/lib -ligl
 LIBIGL_INC=-I$(LIBIGL)/include
 
 ifndef ANTTWEAKBAR_INC
@@ -221,24 +224,31 @@ ifndef ANTTWEAKBAR_LIB
 endif
 
 ifndef SINGULAR_VALUE_DECOMPOSITION_INC
-				SINGULAR_VALUE_DECOMPOSITION_INC=\
-				-I$(LIBIGL)/external/Singular_Value_Decomposition/
+	SINGULAR_VALUE_DECOMPOSITION_INC=\
+	-I$(LIBIGL)/external/Singular_Value_Decomposition/
 endif
 
 ifndef TETGEN
 # By default I'm using the libigl version. Adjust accordingly
-				TETGEN=$(LIBIGL)/external/tetgen
-				TETGEN_LIB=-L$(TETGEN) -ligltetgen -ltet 
-				TETGEN_INC=-I$(TETGEN)
+	TETGEN=$(LIBIGL)/external/tetgen
+	TETGEN_LIB=-L$(TETGEN) -ligltetgen -ltet 
+	TETGEN_INC=-I$(TETGEN)
 endif
 
 ifndef EMBREE
-				EMBREE=$(LIBIGL)/external/embree
-				EMBREE_INC=-I$(EMBREE)/ -I$(EMBREE)/include
-				EMBREE_LIB=-L$(EMBREE)/build -lembree -lsys
+	EMBREE=$(LIBIGL)/external/embree
+	EMBREE_INC=-I$(EMBREE)/ -I$(EMBREE)/include
+	EMBREE_LIB=-L$(EMBREE)/build -lembree -lsys
 endif
 ifndef YIMG
-				YIMG=$(LIBIGL)/external/yimg
-				YIMG_LIB=-L$(YIMG) -lyimg -lz -L/usr/X11/lib -L$(DEFAULT_PREFIX)/lib -lpng -bind_at_load
-				YIMG_INC=-I/usr/X11/include -I$(YIMG) 
+	YIMG=$(LIBIGL)/external/yimg
+	YIMG_LIB=-L$(YIMG) -lyimg -lz -L/usr/X11/lib -L$(DEFAULT_PREFIX)/lib -lpng -bind_at_load
+	YIMG_INC=-I/usr/X11/include -I$(YIMG) 
 endif
+
+ifdef LIBIGL_USE_STATIC_LIBRARY
+	CFLAGS += -DIGL_STATIC_LIBRARY
+	LIBIGL_LIB=-L$(LIBIGL)/lib -ligl
+endif
+
+OPTFLAGS+=-O3 -DNDEBUG $(OPENMP)

examples/camera/example.cpp

@@ -14,6 +14,8 @@
 #include <igl/readOFF.h>
 #include <igl/per_face_normals.h>
 #include <igl/draw_floor.h>
+#include <igl/project.h>
+#include <igl/unproject.h>
 
 #include <Eigen/Core>
 #include <Eigen/Geometry>
@@ -24,6 +26,22 @@
 #include <GL/glut.h>
 #endif
 
+#ifndef GLUT_WHEEL_UP
+#define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
+#define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
+#define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
+#define GLUT_WHEEL_LEFT  6
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
+#endif
+
 #include <vector>
 #include <stack>
 #include <iostream>
@@ -229,26 +247,6 @@ void draw_scene(const igl::Camera & v_camera,
 
   glMatrixMode(GL_PROJECTION);
   glPushMatrix();
-  glLoadIdentity();
-  if(v_camera.m_angle > IGL_CAMERA_MIN_ANGLE)
-  {
-    gluPerspective(v_camera.m_angle,v_camera.m_aspect,v_camera.m_near,v_camera.m_far);
-  }else
-  {
-    glOrtho(
-      -0.5*v_camera.m_aspect,
-      0.5*v_camera.m_aspect,
-      -0.5,
-      0.5,
-      v_camera.m_near,
-      v_camera.m_far);
-  }
-  //{
-  //  Matrix4d m;
-  //  glGetDoublev(GL_PROJECTION_MATRIX,m.data());
-  //  cout<<matlab_format(m,"glu")<<endl;
-  //}
-
   glLoadIdentity();
   glMultMatrixd(v_camera.projection().data());
   //{
@@ -258,14 +256,13 @@ void draw_scene(const igl::Camera & v_camera,
   //}
   glMatrixMode(GL_MODELVIEW);
   glPushMatrix();
-  //glLoadIdentity();
-  //gluLookAt(
-  //  v_camera.eye()(0), v_camera.eye()(1), v_camera.eye()(2),
-  //  v_camera.at()(0), v_camera.at()(1), v_camera.at()(2),
-  //  v_camera.up()(0), v_camera.up()(1), v_camera.up()(2));
   glLoadIdentity();
-  glMultMatrixd(v_camera.inverse().matrix().data());
-
+  gluLookAt(
+    v_camera.eye()(0), v_camera.eye()(1), v_camera.eye()(2),
+    v_camera.at()(0), v_camera.at()(1), v_camera.at()(2),
+    v_camera.up()(0), v_camera.up()(1), v_camera.up()(2));
+  //glLoadIdentity();
+  //glMultMatrixd(v_camera.inverse().matrix().data());
 
   for(int c = 0;c<(int)s.cameras.size();c++)
   {
@@ -431,7 +428,6 @@ void display()
 
   TwDraw();
   glutSwapBuffers();
-  glutPostRedisplay();
 }
 
 
@@ -475,6 +471,7 @@ void mouse_wheel(int wheel, int direction, int mouse_x, int mouse_y)
       camera.dolly((wheel==0?Vector3d(0,0,1):Vector3d(-1,0,0))*0.1*direction);
       break;
   }
+  glutPostRedisplay();
 }
 
 void mouse(int glutButton, int glutState, int mouse_x, int mouse_y)
@@ -509,6 +506,7 @@ void mouse(int glutButton, int glutState, int mouse_x, int mouse_y)
         break;
       }
       break;
+    }
 #ifdef GLUT_WHEEL_DOWN
     // Scroll down
     case GLUT_WHEEL_DOWN:
@@ -541,8 +539,8 @@ void mouse(int glutButton, int glutState, int mouse_x, int mouse_y)
       break;
     }
 #endif
-    }
   }
+  glutPostRedisplay();
 }
 
 void mouse_drag(int mouse_x, int mouse_y)
@@ -595,6 +593,7 @@ void mouse_drag(int mouse_x, int mouse_y)
         break;
     }
   }
+  glutPostRedisplay();
 }
 
 void key(unsigned char key, int mouse_x, int mouse_y)
@@ -609,6 +608,12 @@ void key(unsigned char key, int mouse_x, int mouse_y)
     // ^C
     case char(3):
       exit(0);
+    case 'o':
+    case 'O':
+      {
+        s.cameras[0].m_orthographic ^= true;
+        break;
+      }
     case 'z':
     case 'Z':
       if(mod & GLUT_ACTIVE_COMMAND)
@@ -628,6 +633,7 @@ void key(unsigned char key, int mouse_x, int mouse_y)
         cout<<"Unknown key command: "<<key<<" "<<int(key)<<endl;
       }
   }
+  glutPostRedisplay();
 }
 
 int main(int argc, char * argv[])

examples/convertmesh/Makefile

@@ -0,0 +1,21 @@
+.PHONY: all
+
+# Shared flags etc.
+include ../../build/Makefile.conf
+
+all: convertmesh
+
+.PHONY: convertmesh
+
+INC=$(LIBIGL_INC) $(EIGEN3_INC) 
+LIB=$(LIBIGL_LIB) 
+OPTFLAGS+=-O3 -DNDEBUG $(OPENMP)
+
+convertmesh: convertmesh.o
+	g++ $(OPTFLAGS) $(AFLAGS) $(CFLAGS) $(LIB) -o convertmesh convertmesh.o 
+
+convertmesh.o: convertmesh.cpp
+	g++ $(OPTFLAGS) $(AFLAGS) $(CFLAGS) -c convertmesh.cpp -o convertmesh.o $(INC)
+clean:
+	rm -f convertmesh.o
+	rm -f convertmesh

examples/convertmesh/convertmesh.cpp

@@ -0,0 +1,31 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/write_triangle_mesh.h>
+#include <string>
+#include <iostream>
+int main(int argc, char * argv[])
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  MatrixXd V;
+  MatrixXi F;
+  string in,out;
+  switch(argc)
+  {
+    case 3:
+      in = argv[1];
+      out = argv[2];
+      break;
+    default:
+      cerr<<R"(
+USAGE:
+  convertmesh input.[mesh|obj|off|ply|stl|wrl] output.[mesh|obj|off|ply|stl|wrl]
+
+  Note: .ply and .stl outputs are binary.
+)";
+    return EXIT_FAILURE;
+  }
+  return 
+    read_triangle_mesh(in,V,F) && write_triangle_mesh(out,V,F,false) ? 
+    EXIT_SUCCESS : EXIT_FAILURE;
+}

examples/flare-eyes/example.cpp

@@ -35,6 +35,22 @@
 #include <GL/glut.h>
 #endif
 
+#ifndef GLUT_WHEEL_UP
+#define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
+#define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
+#define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
+#define GLUT_WHEEL_LEFT  6
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
+#endif
+
 #include <string>
 #include <vector>
 #include <iomanip>

examples/patches/Makefile

@@ -2,7 +2,9 @@
 
 # Shared flags etc.
 include ../../build/Makefile.conf
-LIBIGL_LIB+=-liglembree -liglboost
+ifdef LIBIGL_USE_STATIC_LIBRARY
+	LIBIGL_LIB+=-liglembree -liglboost
+endif
 
 all: obj example
 
@@ -17,16 +19,16 @@ CPP_FILES=$(wildcard ./*.cpp)
 OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o))) 
 
 example: obj $(OBJ_FILES)
-	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -o example $(OBJ_FILES) $(LIB)
+	g++ $(OPTFLAGS) $(OPENMP) $(AFLAGS) $(CFLAGS) -o example $(OBJ_FILES) $(LIB)
 
 obj:
 	mkdir -p obj
 
 obj/%.o: %.cpp
-	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
+	g++ $(OPTFLAGS) $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
 
 obj/%.o: %.cpp %.h
-	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
+	g++ $(OPTFLAGS) $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
 
 clean:
 	rm -f $(OBJ_FILES)

examples/patches/example.cpp

@@ -43,10 +43,18 @@
 
 #ifndef GLUT_WHEEL_UP
 #define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
 #define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
 #define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
 #define GLUT_WHEEL_LEFT  6
-#define GLUT_ACTIVE_COMMAND 1
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
 #endif
 
 #include <ctime>

examples/patches/temp.rbr

@@ -1,7 +0,0 @@
-camera_rotation: TW_TYPE_QUAT4D 0.152085375459856 -0.590495112202136 -0.114509193857288 0.784266029838584
-center_type: CenterType orbit
-wireframe_visible: TW_TYPE_BOOLCPP 1
-fill_visible: TW_TYPE_BOOLCPP 1
-rotation_type: RotationType two-axis-valuator-fixed-up
-orient_method: OrientMethod ambient-occlusion
-

examples/quicklook-mesh/Info.plist

@@ -2,6 +2,8 @@
 <!DOCTYPE plist PUBLIC "-//Apple//DTD PLIST 1.0//EN" "http://www.apple.com/DTDs/PropertyList-1.0.dtd">
 <plist version="1.0">
 <dict>
+	<key>BuildMachineOSBuild</key>
+	<string>13E28</string>
 	<key>CFBundleDevelopmentRegion</key>
 	<string>English</string>
 	<key>CFBundleDocumentTypes</key>
@@ -11,49 +13,63 @@
 			<string>QLGenerator</string>
 			<key>LSItemContentTypes</key>
 			<array>
-				<string>org.mesh</string>
-				<string>org.obj</string>
-				<string>org.off</string>
-				<string>org.ply</string>
-				<string>org.stl</string>
-				<string>org.wrl</string>
+				<string>com.mesh.mesh</string>
+				<string>com.mesh.obj</string>
+				<string>com.mesh.off</string>
+				<string>com.mesh.ply</string>
+				<string>com.mesh.stl</string>
+				<string>com.mesh.wrl</string>
 			</array>
 		</dict>
 	</array>
 	<key>CFBundleExecutable</key>
 	<string>Mesh</string>
-	<key>CFBundleIconFile</key>
-	<string></string>
 	<key>CFBundleIdentifier</key>
-	<string>org.mesh.quicklook</string>
+	<string>com.mesh.quicklook</string>
 	<key>CFBundleInfoDictionaryVersion</key>
 	<string>6.0</string>
 	<key>CFBundleName</key>
 	<string>Mesh</string>
 	<key>CFBundleShortVersionString</key>
-	<string>1.1</string>
+	<string>1.0.2</string>
 	<key>CFBundleVersion</key>
-	<string>1.2</string>
+	<string>1.0</string>
 	<key>CFPlugInDynamicRegisterFunction</key>
 	<string></string>
 	<key>CFPlugInDynamicRegistration</key>
-	<false/>
+	<string>NO</string>
 	<key>CFPlugInFactories</key>
 	<dict>
-		<key>C00EB589-55DF-45ED-B4B8-8AFFA367CFCF</key>
+		<key>8D269336-626C-4F43-903F-D2F6AE5C4976</key>
 		<string>QuickLookGeneratorPluginFactory</string>
 	</dict>
 	<key>CFPlugInTypes</key>
 	<dict>
 		<key>5E2D9680-5022-40FA-B806-43349622E5B9</key>
 		<array>
-			<string>C00EB589-55DF-45ED-B4B8-8AFFA367CFCF</string>
+			<string>8D269336-626C-4F43-903F-D2F6AE5C4976</string>
 		</array>
 	</dict>
 	<key>CFPlugInUnloadFunction</key>
 	<string></string>
+	<key>DTCompiler</key>
+	<string>com.apple.compilers.llvm.clang.1_0</string>
+	<key>DTPlatformBuild</key>
+	<string>5B1008</string>
+	<key>DTPlatformVersion</key>
+	<string>GM</string>
+	<key>DTSDKBuild</key>
+	<string>13C64</string>
+	<key>DTSDKName</key>
+	<string>macosx10.9</string>
+	<key>DTXcode</key>
+	<string>0511</string>
+	<key>DTXcodeBuild</key>
+	<string>5B1008</string>
+	<key>LSApplicationCategoryType</key>
+	<string></string>
 	<key>NSHumanReadableCopyright</key>
-	<string>Copyright 2013 Alec Jacobson</string>
+	<string>Copyright © 2015 Alec Jacobson</string>
 	<key>QLNeedsToBeRunInMainThread</key>
 	<true/>
 	<key>QLPreviewHeight</key>
@@ -63,18 +79,19 @@
 	<key>QLSupportsConcurrentRequests</key>
 	<true/>
 	<key>QLThumbnailMinimumSize</key>
-	<integer>32</integer>
-	<key>UTImportedTypeDeclarations</key>
+	<integer>10</integer>
+	<key>UTExportedTypeDeclarations</key>
 	<array>
+
 		<dict>
 			<key>UTTypeConformsTo</key>
 			<array>
-				<string>public.image</string>
+				<string>public.data</string>
 			</array>
 			<key>UTTypeDescription</key>
 			<string>3d tetrahedral mesh</string>
 			<key>UTTypeIdentifier</key>
-			<string>org.mesh</string>
+			<string>com.mesh.mesh</string>
 			<key>UTTypeReferenceURL</key>
 			<string>http://www.ann.jussieu.fr/frey/publications/RT-0253.pdf#page=33</string>
 			<key>UTTypeTagSpecification</key>
@@ -85,16 +102,15 @@
 				</array>
 			</dict>
 		</dict>
-
 		<dict>
 			<key>UTTypeConformsTo</key>
 			<array>
-				<string>public.image</string>
+				<string>public.data</string>
 			</array>
 			<key>UTTypeDescription</key>
 			<string>3D model format</string>
 			<key>UTTypeIdentifier</key>
-			<string>org.obj</string>
+			<string>com.mesh.obj</string>
 			<key>UTTypeReferenceURL</key>
 			<string>http://en.wikipedia.org/wiki/Wavefront_.obj_file#File_format</string>
 			<key>UTTypeTagSpecification</key>
@@ -105,16 +121,15 @@
 				</array>
 			</dict>
 		</dict>
-
 		<dict>
 			<key>UTTypeConformsTo</key>
 			<array>
-				<string>public.image</string>
+				<string>public.data</string>
 			</array>
 			<key>UTTypeDescription</key>
 			<string>PLY 3D file</string>
 			<key>UTTypeIdentifier</key>
-			<string>org.ply</string>
+			<string>com.mesh.ply</string>
 			<key>UTTypeReferenceURL</key>
 			<string>http://en.wikipedia.org/wiki/PLY_%28file_format%29</string>
 			<key>UTTypeTagSpecification</key>
@@ -125,16 +140,15 @@
 				</array>
 			</dict>
 		</dict>
-
 		<dict>
 			<key>UTTypeConformsTo</key>
 			<array>
-				<string>public.image</string>
+				<string>public.data</string>
 			</array>
 			<key>UTTypeDescription</key>
 			<string>Stereolithography CAD model</string>
 			<key>UTTypeIdentifier</key>
-			<string>org.stl</string>
+			<string>com.mesh.stl</string>
 			<key>UTTypeReferenceURL</key>
 			<string>http://en.wikipedia.org/wiki/STL_(file_format)</string>
 			<key>UTTypeTagSpecification</key>
@@ -148,12 +162,12 @@
 		<dict>
 			<key>UTTypeConformsTo</key>
 			<array>
-				<string>public.image</string>
+				<string>public.data</string>
 			</array>
 			<key>UTTypeDescription</key>
 			<string>Geomview's polyhedral file format</string>
 			<key>UTTypeIdentifier</key>
-			<string>org.off</string>
+			<string>com.mesh.off</string>
 			<key>UTTypeReferenceURL</key>
 			<string>http://tetgen.berlios.de/fformats.off.html</string>
 			<key>UTTypeTagSpecification</key>
@@ -167,12 +181,12 @@
 		<dict>
 			<key>UTTypeConformsTo</key>
 			<array>
-				<string>public.image</string>
+				<string>public.data</string>
 			</array>
 			<key>UTTypeDescription</key>
 			<string>Virtual Reality Modeling Language</string>
 			<key>UTTypeIdentifier</key>
-			<string>org.wrl</string>
+			<string>com.mesh.wrl</string>
 			<key>UTTypeReferenceURL</key>
 			<string>http://en.wikipedia.org/wiki/VRML#WRL_File_Format</string>
 			<key>UTTypeTagSpecification</key>
@@ -183,6 +197,8 @@
 				</array>
 			</dict>
 		</dict>
+
 	</array>
 </dict>
 </plist>
+

examples/quicklook-mesh/Makefile

@@ -28,6 +28,8 @@ MESA_LIB=-L$(MESA)/lib -lOSMesa -lGL
 OBJC_LIB=-lobjc
 
 all: obj Mesh.qlgenerator
+deploy:
+	dylibbundler -od -b -x ./Mesh.qlgenerator/Contents/MacOS/Mesh -d ./Mesh.qlgenerator/Contents/libs -p @loader_path/../libs/
 install:
 	rm -rf /Library/QuickLook/Mesh.qlgenerator
 	cp -R Mesh.qlgenerator /Library/QuickLook/Mesh.qlgenerator

examples/quicklook-mesh/src/render_to_buffer.cpp

@@ -280,7 +280,24 @@ void display()
       glMaterialfv(GL_BACK, GL_SPECULAR, SILVER_SPECULAR);
       glMaterialf (GL_BACK, GL_SHININESS, 128);
     //}
-    draw_mesh(V,F,N);
+    if(F.rows() == 0)
+    {
+      glPointSize(1.+ 20./log10(V.rows()));
+      const bool has_normals = N.rows() == V.rows();
+      glBegin(GL_POINTS);
+      for(int v = 0;v<V.rows();v++)
+      {
+        if(has_normals)
+        {
+          glNormal3d(N(v,0),N(v,1),N(v,2));
+        }
+        glVertex3d(V(v,0),V(v,1),V(v,2));
+      }
+      glEnd();
+    }else
+    {
+      draw_mesh(V,F,N);
+    }
     //if(invert)
     //{
     //  glFrontFace(GL_CCW);
@@ -355,6 +372,7 @@ bool render_to_buffer(
     // Convert extension to lower case
     if(!igl::readOBJ(filename,vV,vTC,vN,vF,vFTC,vFN))
     {
+      cerr<<"readOBJ failed."<<endl;
       red(width,height,buffer);
       return false;
     }
@@ -363,6 +381,7 @@ bool render_to_buffer(
     // Convert extension to lower case
     if(!igl::readOFF(filename,vV,vF,vN))
     {
+      cerr<<"readOFF failed."<<endl;
       red(width,height,buffer);
       return false;
     }
@@ -408,16 +427,25 @@ bool render_to_buffer(
   {
     if(!list_to_matrix(vV,V))
     {
+      cerr<<"list_to_matrix failed."<<endl;
       red(width,height,buffer);
       return false;
     }
+    if(!list_to_matrix(vN,N))
+    {
+      // silently continue
+      N.resize(0,0);
+    }
     polygon_mesh_to_triangle_mesh(vF,F);
   }
   cout<<"IO: "<<(get_seconds()-ts)<<"s"<<endl;
   ts = get_seconds();
 
-  // Computer already normalized per triangle normals
-  per_face_normals(V,F,N);
+  // Compute already normalized per triangle normals
+  if(N.rows() != V.rows() && N.rows() != F.rows())
+  {
+    per_face_normals(V,F,N);
+  }
   //Vmean = 0.5*(V.colwise().maxCoeff()+V.colwise().minCoeff());
   Vmax = V.colwise().maxCoeff();
   Vmin = V.colwise().minCoeff();

examples/randomly-sample-mesh/example.cpp

@@ -32,6 +32,22 @@
 #include <GL/glut.h>
 #endif
 
+#ifndef GLUT_WHEEL_UP
+#define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
+#define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
+#define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
+#define GLUT_WHEEL_LEFT  6
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
+#endif
+
 #include <string>
 #include <vector>
 #include <stack>

examples/rotate-widget/example.cpp

@@ -33,6 +33,22 @@
 # include <GL/glut.h>
 #endif
 
+#ifndef GLUT_WHEEL_UP
+#define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
+#define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
+#define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
+#define GLUT_WHEEL_LEFT  6
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
+#endif
+
 #include <string>
 #include <vector>
 #include <list>

examples/scene-rotation/example.cpp

@@ -36,6 +36,22 @@
 #include <GL/glut.h>
 #endif
 
+#ifndef GLUT_WHEEL_UP
+#define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
+#define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
+#define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
+#define GLUT_WHEEL_LEFT  6
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
+#endif
+
 #include <string>
 #include <vector>
 #include <stack>

examples/shadow-mapping/example.cpp

@@ -43,10 +43,18 @@
 
 #ifndef GLUT_WHEEL_UP
 #define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
 #define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
 #define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
 #define GLUT_WHEEL_LEFT  6
-#define GLUT_ACTIVE_COMMAND 1
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
 #endif
 
 #include <ctime>

examples/skeleton-builder/Makefile

@@ -4,9 +4,9 @@
 include ../../build/Makefile.conf
 LIBIGL_LIB+=-liglembree
 
-all: obj example
+all: obj skeleton_builder
 
-.PHONY: example
+.PHONY: skeleton_builder
 
 INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(GLUT_INC) ${TETGEN_INC} $(MOSEK_INC) $(EMBREE_INC)
 LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(TETGEN_LIB) $(MOSEK_LIB) $(EMBREE_LIB)
@@ -16,8 +16,8 @@ OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
 
 CFLAGS+=-std=c++11 -g
 
-example: obj $(OBJ_FILES)
-	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -o example $(LIB) $(OBJ_FILES) 
+skeleton_builder: obj $(OBJ_FILES)
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -o skeleton_builder $(LIB) $(OBJ_FILES) 
 
 obj:
 	mkdir -p obj
@@ -30,4 +30,4 @@ obj/%.o: %.cpp %.h
 
 clean:
 	rm -f $(OBJ_FILES)
-	rm -f example
+	rm -f skeleton_builder

examples/skeleton-builder/example.cpp

@@ -51,6 +51,22 @@
 #include <GL/glut.h>
 #endif
 
+#ifndef GLUT_WHEEL_UP
+#define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
+#define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
+#define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
+#define GLUT_WHEEL_LEFT  6
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
+#endif
+
 #include <string>
 #include <vector>
 #include <queue>

examples/skeleton-posing/Makefile → examples/skeleton-poser/Makefile

@@ -4,9 +4,9 @@
 include ../../build/Makefile.conf
 LIBIGL_LIB+=-liglbbw -liglcgal
 
-all: obj example
+all: obj skeleton-poser
 
-.PHONY: example
+.PHONY: skeleton-poser
 
 ifdef IGL_NO_MOSEK
 CFLAGS+=-DIGL_NO_MOSEK
@@ -39,8 +39,8 @@ OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
 
 CFLAGS+=-std=c++11 -g
 
-example: obj $(OBJ_FILES)
-	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -o example $(LIB) $(OBJ_FILES) 
+skeleton-poser: obj $(OBJ_FILES)
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -o skeleton-poser $(LIB) $(OBJ_FILES) 
 
 obj:
 	mkdir -p obj
@@ -53,4 +53,4 @@ obj/%.o: %.cpp %.h
 
 clean:
 	rm -f $(OBJ_FILES)
-	rm -f example
+	rm -f skeleton-poser

examples/skeleton-posing/example.cpp → examples/skeleton-poser/example.cpp

@@ -49,13 +49,26 @@
 
 #ifdef __APPLE__
 #include <GLUT/glut.h>
-#ifndef GLUT_ACTIVE_COMMAND
-#  define GLUT_ACTIVE_COMMAND 9
-#endif
 #else
 #include <GL/glut.h>
 #endif
 
+#ifndef GLUT_WHEEL_UP
+#define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
+#define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
+#define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
+#define GLUT_WHEEL_LEFT  6
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
+#endif
+
 #include <string>
 #include <vector>
 #include <queue>

examples/skeleton/example.cpp

@@ -40,6 +40,22 @@
 #include <GL/glut.h>
 #endif
 
+#ifndef GLUT_WHEEL_UP
+#define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
+#define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
+#define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
+#define GLUT_WHEEL_LEFT  6
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
+#endif
+
 #include <string>
 #include <vector>
 #include <stack>

examples/upright/Makefile

@@ -11,7 +11,7 @@ include ../..//build/Makefile.conf
 INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC)
 LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB)
 
-CFLAGS+=-g
+CFLAGS+=-g $(OPTFLAGS)
 
 upright: upright.o
 	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -o upright upright.o $(LIB)

examples/upright/example.cpp

@@ -8,13 +8,7 @@
 // normals, texture coordinates, etc. in the input file will be ignored and
 // lost in the output file.
 //
-#include <igl/readOBJ.h>
-#include <igl/writeOBJ.h>
-#include <igl/writeOFF.h>
-#include <igl/readWRL.h>
-#include <igl/polygon_mesh_to_triangle_mesh.h>
-#include <igl/readOFF.h>
-#include <igl/readMESH.h>
+#include <igl/read_triangle_mesh.h>
 #include <igl/draw_mesh.h>
 #include <igl/draw_floor.h>
 #include <igl/pathinfo.h>
@@ -23,6 +17,7 @@
 #include <igl/material_colors.h>
 #include <igl/trackball.h>
 #include <igl/snap_to_canonical_view_quat.h>
+#include <igl/write_triangle_mesh.h>
 #include <igl/REDRUM.h>
 
 #include <Eigen/Core>
@@ -34,6 +29,22 @@
 #include <GL/glut.h>
 #endif
 
+#ifndef GLUT_WHEEL_UP
+#define GLUT_WHEEL_UP    3
+#endif
+#ifndef GLUT_WHEEL_DOWN
+#define GLUT_WHEEL_DOWN  4
+#endif
+#ifndef GLUT_WHEEL_RIGHT
+#define GLUT_WHEEL_RIGHT 5
+#endif
+#ifndef GLUT_WHEEL_LEFT
+#define GLUT_WHEEL_LEFT  6
+#endif
+#ifndef GLUT_ACTIVE_COMMAND
+#define GLUT_ACTIVE_COMMAND 8
+#endif
+
 #include <string>
 #include <vector>
 #include <stack>
@@ -301,47 +312,15 @@ bool save()
   using namespace std;
   using namespace igl;
   using namespace Eigen;
-  // dirname, basename, extension and filename
-  string d,b,ext,f;
-  pathinfo(out_filename,d,b,ext,f);
-  // Convert extension to lower case
-  transform(ext.begin(), ext.end(), ext.begin(), ::tolower);
-  if(ext == "obj")
+  if(write_triangle_mesh(out_filename,s.V,s.F))
   {
-    // Convert extension to lower case
-    if(!igl::writeOBJ(out_filename,s.V,s.F))
-    {
-      return false;
-    }
-  }else if(ext == "off")
+    cout<<GREENGIN("Current baked model written to "+out_filename+".")<<endl;
+    return true;
+  }else
   {
-    // Convert extension to lower case
-    if(!igl::writeOFF(out_filename,s.V,s.F))
-    {
-      return false;
-    }
-  //}else if(ext == "wrl")
-  //{
-  //  // Convert extension to lower case
-  //  if(!igl::readWRL(filename,vV,vF))
-  //  {
-  //    return 1;
-  //  }
-  //}else
-  //{
-  //  // Convert extension to lower case
-  //  MatrixXi T;
-  //  if(!igl::readMESH(filename,V,T,F))
-  //  {
-  //    return 1;
-  //  }
-  //  //if(F.size() > T.size() || F.size() == 0)
-  //  {
-  //    boundary_facets(T,F);
-  //  }
+    cout<<REDRUM("Current baked model failed to write to "+out_filename+".")<<endl;
+    return false;
   }
-  cout<<GREENGIN("Current baked model written to "+out_filename+".")<<endl;
-  return true;
 }
 
 
@@ -455,56 +434,11 @@ int main(int argc, char * argv[])
   // Read and prepare mesh
   string filename = argv[1];
   out_filename = argv[2];
-  // dirname, basename, extension and filename
-  string d,b,ext,f;
-  pathinfo(filename,d,b,ext,f);
-  // Convert extension to lower case
-  transform(ext.begin(), ext.end(), ext.begin(), ::tolower);
-  vector<vector<double > > vV,vN,vTC;
-  vector<vector<int > > vF,vFTC,vFN;
-  if(ext == "obj")
-  {
-    // Convert extension to lower case
-    if(!igl::readOBJ(filename,vV,vTC,vN,vF,vFTC,vFN))
-    {
-      return 1;
-    }
-  }else if(ext == "off")
-  {
-    // Convert extension to lower case
-    if(!igl::readOFF(filename,vV,vF,vN))
-    {
-      return 1;
-    }
-  }else if(ext == "wrl")
+  if(!read_triangle_mesh(filename,s.V,s.F))
   {
-    // Convert extension to lower case
-    if(!igl::readWRL(filename,vV,vF))
-    {
-      return 1;
-    }
-  //}else
-  //{
-  //  // Convert extension to lower case
-  //  MatrixXi T;
-  //  if(!igl::readMESH(filename,V,T,F))
-  //  {
-  //    return 1;
-  //  }
-  //  //if(F.size() > T.size() || F.size() == 0)
-  //  {
-  //    boundary_facets(T,F);
-  //  }
+    cout<<"Failed to read from "<<filename<<"."<<endl;
+    return EXIT_FAILURE;
   }
-  if(vV.size() > 0)
-  {
-    if(!list_to_matrix(vV,s.V))
-    {
-      return 1;
-    }
-    polygon_mesh_to_triangle_mesh(vF,s.F);
-  }
-
   init_relative();
 
   // Init glut
@@ -520,5 +454,5 @@ int main(int argc, char * argv[])
   //glutPassiveMotionFunc(mouse_move);
   glutMainLoop();
 
-  return 0;
+  return EXIT_SUCCESS;
 }

file-formats/index.html

@@ -26,7 +26,7 @@
     href=http://tetgen.berlios.de/fformats.node.html>TetGen</a>, and <a
     href=http://www.cs.berkeley.edu/~jrs/stellar/#fileformats>Stellar</a>.
     </li>
-    <li><a href="http://tetgen.berlios.de/fformats.off.html">.off</a> Geomview's polyhedral file format</li>
+    <li><a href="http://wias-berlin.de/software/tetgen/fformats.off.html">.off</a> Geomview's polyhedral file format</li>
     <li><a
     href="http://en.wikipedia.org/wiki/Wavefront_.obj_file#File_format">.obj</a>
     Wavefront object file format. Usually unsafe to assume anything more than

include/igl/AABB.h

@@ -0,0 +1,820 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_AABB_H
+#define IGL_AABB_H
+
+// Implementation of semi-general purpose axis-aligned bounding box hierarchy.
+// The mesh (V,Ele) is stored and managed by the caller and each routine here
+// simply takes it as references (it better not change between calls).
+//
+// It's a little annoying that the Dimension is a template parameter and not
+// picked up at run time from V. This leads to duplicated code for 2d/3d (up to
+// dim).
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <vector>
+namespace igl
+{
+  template <typename DerivedV, int DIM>
+    class AABB 
+    {
+      public:
+        typedef typename DerivedV::Scalar Scalar;
+        typedef Eigen::Matrix<Scalar,1,DIM> RowVectorDIMS;
+        typedef Eigen::Matrix<Scalar,DIM,1> VectorDIMS;
+        typedef Eigen::Matrix<Scalar,Eigen::Dynamic,DIM> MatrixXDIMS;
+        // Shared pointers are slower...
+        AABB * m_left;
+        AABB * m_right;
+        Eigen::AlignedBox<Scalar,DIM> m_box;
+        // -1 non-leaf
+        int m_primitive;
+        AABB():
+          m_left(NULL), m_right(NULL),
+          m_box(), m_primitive(-1)
+      {}
+        // http://stackoverflow.com/a/3279550/148668
+        AABB(const AABB& other):
+          m_left(other.m_left ? new AABB(*other.m_left) : NULL),
+          m_right(other.m_right ? new AABB(*other.m_right) : NULL),
+          m_box(other.m_box),
+          m_primitive(other.m_primitive)
+          {
+          }
+        // copy-swap idiom
+        friend void swap(AABB& first, AABB& second)
+        {
+          // Enable ADL
+          using std::swap;
+          swap(first.m_left,second.m_left);
+          swap(first.m_right,second.m_right);
+          swap(first.m_box,second.m_box);
+          swap(first.m_primitive,second.m_primitive);
+        }
+        // Pass-by-value (aka copy)
+        AABB& operator=(AABB other)
+        {
+          swap(*this,other);
+          return *this;
+        }
+        AABB(AABB&& other):
+          // initialize via default constructor
+          AABB() 
+      {
+        swap(*this,other);
+      }
+        // Seems like there should have been an elegant solution to this using
+        // the copy-swap idiom above:
+        inline void deinit()
+        {
+          m_primitive = -1;
+          m_box = Eigen::AlignedBox<Scalar,DIM>();
+          delete m_left;
+          m_left = NULL;
+          delete m_right;
+          m_right = NULL;
+        }
+        ~AABB()
+        {
+          deinit();
+        }
+        // Build an Axis-Aligned Bounding Box tree for a given mesh and given
+        // serialization of a previous AABB tree.
+        //
+        // Inputs:
+        //   V  #V by dim list of mesh vertex positions. 
+        //   Ele  #Ele by dim+1 list of mesh indices into #V. 
+        //   bb_mins  max_tree by dim list of bounding box min corner positions
+        //   bb_maxs  max_tree by dim list of bounding box max corner positions
+        //   elements  max_tree list of element or (not leaf id) indices into Ele
+        //   i  recursive call index {0}
+        template <typename Derivedbb_mins, typename Derivedbb_maxs>
+          inline void init(
+              const Eigen::PlainObjectBase<DerivedV> & V,
+              const Eigen::MatrixXi & Ele, 
+              const Eigen::PlainObjectBase<Derivedbb_mins> & bb_mins,
+              const Eigen::PlainObjectBase<Derivedbb_maxs> & bb_maxs,
+              const Eigen::VectorXi & elements,
+              const int i = 0);
+        // Wrapper for root with empty serialization
+        inline void init(
+            const Eigen::PlainObjectBase<DerivedV> & V,
+            const Eigen::MatrixXi & Ele);
+        // Build an Axis-Aligned Bounding Box tree for a given mesh.
+        //
+        // Inputs:
+        //   V  #V by dim list of mesh vertex positions. 
+        //   Ele  #Ele by dim+1 list of mesh indices into #V. 
+        //   SI  #Ele by dim list revealing for each coordinate where Ele's
+        //     barycenters would be sorted: SI(e,d) = i --> the dth coordinate of
+        //     the barycenter of the eth element would be placed at position i in a
+        //     sorted list.
+        //   I  #I list of indices into Ele of elements to include (for recursive
+        //     calls)
+        // 
+        inline void init(
+            const Eigen::PlainObjectBase<DerivedV> & V,
+            const Eigen::MatrixXi & Ele, 
+            const Eigen::MatrixXi & SI,
+            const Eigen::VectorXi & I);
+        // Return whether at leaf node
+        inline bool is_leaf() const;
+        // Find the indices of elements containing given point.
+        //
+        // Inputs:
+        //   V  #V by dim list of mesh vertex positions. **Should be same as used to
+        //     construct mesh.**
+        //   Ele  #Ele by dim+1 list of mesh indices into #V. **Should be same as used to
+        //     construct mesh.**
+        //   q  dim row-vector query position
+        //   first  whether to only return first element containing q
+        // Returns:
+        //   list of indices of elements containing q
+        template <typename Derivedq>
+        inline std::vector<int> find(
+            const Eigen::PlainObjectBase<DerivedV> & V,
+            const Eigen::MatrixXi & Ele, 
+            const Eigen::PlainObjectBase<Derivedq> & q,
+            const bool first=false) const;
+
+        // If number of elements m then total tree size should be 2*h where h is
+        // the deepest depth 2^ceil(log(#Ele*2-1))
+        inline int subtree_size() const;
+
+        // Serialize this class into 3 arrays (so we can pass it pack to matlab)
+        //
+        // Outputs:
+        //   bb_mins  max_tree by dim list of bounding box min corner positions
+        //   bb_maxs  max_tree by dim list of bounding box max corner positions
+        //   elements  max_tree list of element or (not leaf id) indices into Ele
+        //   i  recursive call index into these arrays {0}
+        template <typename Derivedbb_mins, typename Derivedbb_maxs>
+          inline void serialize(
+              Eigen::PlainObjectBase<Derivedbb_mins> & bb_mins,
+              Eigen::PlainObjectBase<Derivedbb_maxs> & bb_maxs,
+              Eigen::VectorXi & elements,
+              const int i = 0) const;
+        // Compute squared distance to a query point
+        //
+        // Inputs:
+        //   V  #V by dim list of vertex positions
+        //   Ele  #Ele by dim list of simplex indices
+        //   P  3 list of query point coordinates
+        //   min_sqr_d  current minimum squared distance (only find distances
+        //   less than this)
+        // Outputs:
+        //   I  #P list of facet indices corresponding to smallest distances
+        //   C  #P by 3 list of closest points
+        // Returns squared distance
+        //
+        // Known bugs: currently assumes Elements are triangles regardless of
+        // dimension.
+        inline Scalar squared_distance(
+            const Eigen::PlainObjectBase<DerivedV> & V,
+            const Eigen::MatrixXi & Ele, 
+            const RowVectorDIMS & p,
+            int & i,
+            RowVectorDIMS & c) const;
+      private:
+        inline Scalar squared_distance(
+            const Eigen::PlainObjectBase<DerivedV> & V,
+            const Eigen::MatrixXi & Ele, 
+            const RowVectorDIMS & p,
+            const Scalar min_sqr_d,
+            int & i,
+            RowVectorDIMS & c) const;
+      public:
+        template <
+          typename DerivedP, 
+                   typename DerivedsqrD, 
+                   typename DerivedI, 
+                   typename DerivedC>
+                     inline void squared_distance(
+                         const Eigen::PlainObjectBase<DerivedV> & V,
+                         const Eigen::MatrixXi & Ele, 
+                         const Eigen::PlainObjectBase<DerivedP> & P,
+                         Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+                         Eigen::PlainObjectBase<DerivedI> & I,
+                         Eigen::PlainObjectBase<DerivedC> & C) const;
+      private:
+        // Helper function for leaves: works in-place on sqr_d
+        inline void leaf_squared_distance(
+            const Eigen::PlainObjectBase<DerivedV> & V,
+            const Eigen::MatrixXi & Ele, 
+            const RowVectorDIMS & p,
+            Scalar & sqr_d,
+            int & i,
+            RowVectorDIMS & c) const;
+        inline void set_min(
+            const RowVectorDIMS & p,
+            const Scalar sqr_d_candidate,
+            const int i_candidate,
+            const RowVectorDIMS & c_candidate,
+            Scalar & sqr_d,
+            int & i,
+            RowVectorDIMS & c) const;
+      public:
+        template <int SS>
+          static
+          void barycentric_coordinates(
+              const RowVectorDIMS & p, 
+              const RowVectorDIMS & a, 
+              const RowVectorDIMS & b, 
+              const RowVectorDIMS & c,
+              Eigen::Matrix<Scalar,1,SS> & bary);
+    };
+}
+
+// Implementation
+#include "EPS.h"
+#include "barycenter.h"
+#include "colon.h"
+#include "colon.h"
+#include "doublearea.h"
+#include "matlab_format.h"
+#include "project_to_line_segment.h"
+#include "sort.h"
+#include "volume.h"
+#include <iostream>
+#include <limits>
+
+template <typename DerivedV, int DIM>
+  template <typename Derivedbb_mins, typename Derivedbb_maxs>
+inline void igl::AABB<DerivedV,DIM>::init(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXi & Ele, 
+    const Eigen::PlainObjectBase<Derivedbb_mins> & bb_mins,
+    const Eigen::PlainObjectBase<Derivedbb_maxs> & bb_maxs,
+    const Eigen::VectorXi & elements,
+    const int i)
+{
+  using namespace std;
+  using namespace Eigen;
+  if(bb_mins.size() > 0)
+  {
+    assert(bb_mins.rows() == bb_maxs.rows() && "Serial tree arrays must match");
+    assert(bb_mins.cols() == V.cols() && "Serial tree array dim must match V");
+    assert(bb_mins.cols() == bb_maxs.cols() && "Serial tree arrays must match");
+    assert(bb_mins.rows() == elements.rows() &&
+        "Serial tree arrays must match");
+    // construct from serialization
+    m_box.extend(bb_mins.row(i).transpose());
+    m_box.extend(bb_maxs.row(i).transpose());
+    m_primitive = elements(i);
+    // Not leaf then recurse
+    if(m_primitive == -1)
+    {
+      m_left = new AABB();
+      m_left->init( V,Ele,bb_mins,bb_maxs,elements,2*i+1);
+      m_right = new AABB();
+      m_right->init( V,Ele,bb_mins,bb_maxs,elements,2*i+2);
+    }
+  }else
+  {
+    VectorXi allI = colon<int>(0,Ele.rows()-1);
+    MatrixXDIMS BC;
+    barycenter(V,Ele,BC);
+    MatrixXi SI(BC.rows(),BC.cols());
+    {
+      MatrixXDIMS _;
+      MatrixXi IS;
+      igl::sort(BC,1,true,_,IS);
+      // Need SI(i) to tell which place i would be sorted into
+      const int dim = IS.cols();
+      for(int i = 0;i<IS.rows();i++)
+      {
+        for(int d = 0;d<dim;d++)
+        {
+          SI(IS(i,d),d) = i;
+        }
+      }
+    }
+    init(V,Ele,SI,allI);
+  }
+}
+
+  template <typename DerivedV, int DIM>
+inline void igl::AABB<DerivedV,DIM>::init(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXi & Ele)
+{
+  using namespace Eigen;
+  return init(V,Ele,MatrixXDIMS(),MatrixXDIMS(),VectorXi(),0);
+}
+
+  template <typename DerivedV, int DIM>
+inline void igl::AABB<DerivedV,DIM>::init(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXi & Ele, 
+    const Eigen::MatrixXi & SI,
+    const Eigen::VectorXi & I)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(DIM == V.cols() && "V.cols() should matched declared dimension");
+  const Scalar inf = numeric_limits<Scalar>::infinity();
+  m_box = AlignedBox<Scalar,DIM>();
+  // Compute bounding box
+  for(int i = 0;i<I.rows();i++)
+  {
+    for(int c = 0;c<Ele.cols();c++)
+    {
+      m_box.extend(V.row(Ele(I(i),c)).transpose());
+      m_box.extend(V.row(Ele(I(i),c)).transpose());
+    }
+  }
+  switch(I.size())
+  {
+    case 0:
+      {
+        assert(false);
+      }
+    case 1:
+      {
+        m_primitive = I(0);
+        break;
+      }
+    default:
+      {
+        // Compute longest direction
+        int max_d = -1;
+        m_box.diagonal().maxCoeff(&max_d);
+        // Can't use median on BC directly because many may have same value,
+        // but can use median on sorted BC indices
+        VectorXi SIdI(I.rows());
+        for(int i = 0;i<I.rows();i++)
+        {
+          SIdI(i) = SI(I(i),max_d);
+        }
+        // Since later I use <= I think I don't need to worry about odd/even
+        // Pass by copy to avoid changing input
+        const auto median = [](VectorXi A)->Scalar
+        {
+          size_t n = A.size()/2;
+          nth_element(A.data(),A.data()+n,A.data()+A.size());
+          if(A.rows() % 2 == 1)
+          {
+            return A(n);
+          }else
+          {
+            nth_element(A.data(),A.data()+n-1,A.data()+A.size());
+            return 0.5*(A(n)+A(n-1));
+          }
+        };
+        const Scalar med = median(SIdI);
+        VectorXi LI((I.rows()+1)/2),RI(I.rows()/2);
+        assert(LI.rows()+RI.rows() == I.rows());
+        // Distribute left and right
+        {
+          int li = 0;
+          int ri = 0;
+          for(int i = 0;i<I.rows();i++)
+          {
+            if(SIdI(i)<=med)
+            {
+              LI(li++) = I(i);
+            }else
+            {
+              RI(ri++) = I(i);
+            }
+          }
+        }
+        if(LI.rows()>0)
+        {
+          m_left = new AABB();
+          m_left->init(V,Ele,SI,LI);
+        }
+        if(RI.rows()>0)
+        {
+          m_right = new AABB();
+          m_right->init(V,Ele,SI,RI);
+        }
+      }
+  }
+}
+
+template <typename DerivedV, int DIM>
+inline bool igl::AABB<DerivedV,DIM>::is_leaf() const
+{
+  return m_primitive != -1;
+}
+
+template <typename DerivedV, int DIM>
+template <typename Derivedq>
+inline std::vector<int> igl::AABB<DerivedV,DIM>::find(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXi & Ele, 
+    const Eigen::PlainObjectBase<Derivedq> & q,
+    const bool first) const
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(q.size() == DIM && 
+      "Query dimension should match aabb dimension");
+  assert(Ele.cols() == V.cols()+1 && 
+      "AABB::find only makes sense for (d+1)-simplices");
+  const Scalar epsilon = igl::EPS<Scalar>();
+  // Check if outside bounding box
+  bool inside = m_box.contains(q.transpose());
+  if(!inside)
+  {
+    return std::vector<int>();
+  }
+  assert(m_primitive==-1 || (m_left == NULL && m_right == NULL));
+  if(is_leaf())
+  {
+    // Initialize to some value > -epsilon
+    Scalar a1=1,a2=1,a3=1,a4=1;
+    switch(DIM)
+    {
+      case 3:
+        {
+          // Barycentric coordinates
+          typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
+          const RowVector3S V1 = V.row(Ele(m_primitive,0));
+          const RowVector3S V2 = V.row(Ele(m_primitive,1));
+          const RowVector3S V3 = V.row(Ele(m_primitive,2));
+          const RowVector3S V4 = V.row(Ele(m_primitive,3));
+          a1 = volume_single(V2,V4,V3,(RowVector3S)q);
+          a2 = volume_single(V1,V3,V4,(RowVector3S)q);
+          a3 = volume_single(V1,V4,V2,(RowVector3S)q);
+          a4 = volume_single(V1,V2,V3,(RowVector3S)q);
+          break;
+        }
+      case 2:
+        {
+          // Barycentric coordinates
+          typedef Eigen::Matrix<Scalar,2,1> Vector2S;
+          const Vector2S V1 = V.row(Ele(m_primitive,0));
+          const Vector2S V2 = V.row(Ele(m_primitive,1));
+          const Vector2S V3 = V.row(Ele(m_primitive,2));
+          // Hack for now to keep templates simple. If becomes bottleneck
+          // consider using std::enable_if_t 
+          const Vector2S q2 = q.head(2);
+          Scalar a0 = doublearea_single(V1,V2,V3);
+          a1 = doublearea_single(V1,V2,q2);
+          a2 = doublearea_single(V2,V3,q2);
+          a3 = doublearea_single(V3,V1,q2);
+          break;
+        }
+      default:assert(false);
+    }
+    if(
+        a1>=-epsilon && 
+        a2>=-epsilon && 
+        a3>=-epsilon && 
+        a4>=-epsilon)
+    {
+      return std::vector<int>(1,m_primitive);
+    }else
+    {
+      return std::vector<int>();
+    }
+  }
+  std::vector<int> left = m_left->find(V,Ele,q,first);
+  if(first && !left.empty())
+  {
+    return left;
+  }
+  std::vector<int> right = m_right->find(V,Ele,q,first);
+  if(first)
+  {
+    return right;
+  }
+  left.insert(left.end(),right.begin(),right.end());
+  return left;
+}
+
+template <typename DerivedV, int DIM>
+inline int igl::AABB<DerivedV,DIM>::subtree_size() const
+{
+  // 1 for self
+  int n = 1;
+  int n_left = 0,n_right = 0;
+  if(m_left != NULL)
+  {
+    n_left = m_left->subtree_size();
+  }
+  if(m_right != NULL)
+  {
+    n_right = m_right->subtree_size();
+  }
+  n += 2*std::max(n_left,n_right);
+  return n;
+}
+
+
+template <typename DerivedV, int DIM>
+template <typename Derivedbb_mins, typename Derivedbb_maxs>
+inline void igl::AABB<DerivedV,DIM>::serialize(
+    Eigen::PlainObjectBase<Derivedbb_mins> & bb_mins,
+    Eigen::PlainObjectBase<Derivedbb_maxs> & bb_maxs,
+    Eigen::VectorXi & elements,
+    const int i) const
+{
+  using namespace std;
+  using namespace Eigen;
+  // Calling for root then resize output
+  if(i==0)
+  {
+    const int m = subtree_size();
+    //cout<<"m: "<<m<<endl;
+    bb_mins.resize(m,DIM);
+    bb_maxs.resize(m,DIM);
+    elements.resize(m,1);
+  }
+  //cout<<i<<" ";
+  bb_mins.row(i) = m_box.min();
+  bb_maxs.row(i) = m_box.max();
+  elements(i) = m_primitive;
+  if(m_left != NULL)
+  {
+    m_left->serialize(bb_mins,bb_maxs,elements,2*i+1);
+  }
+  if(m_right != NULL)
+  {
+    m_right->serialize(bb_mins,bb_maxs,elements,2*i+2);
+  }
+}
+
+template <typename DerivedV, int DIM>
+inline typename igl::AABB<DerivedV,DIM>::Scalar 
+igl::AABB<DerivedV,DIM>::squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & p,
+  int & i,
+  RowVectorDIMS & c) const
+{
+  return squared_distance(V,Ele,p,std::numeric_limits<Scalar>::infinity(),i,c);
+}
+
+
+template <typename DerivedV, int DIM>
+inline typename igl::AABB<DerivedV,DIM>::Scalar 
+igl::AABB<DerivedV,DIM>::squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & p,
+  Scalar min_sqr_d,
+  int & i,
+  RowVectorDIMS & c) const
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  Scalar sqr_d = min_sqr_d;
+  assert(DIM == 3 && "Code has only been tested for DIM == 3");
+  assert((Ele.cols() == 3 || Ele.cols() == 2 || Ele.cols() == 1)
+    && "Code has only been tested for simplex sizes 3,2,1");
+
+  assert(m_primitive==-1 || (m_left == NULL && m_right == NULL));
+  if(is_leaf())
+  {
+    leaf_squared_distance(V,Ele,p,sqr_d,i,c);
+  }else
+  {
+    bool looked_left = false;
+    bool looked_right = false;
+    const auto & look_left = [&]()
+    {
+      int i_left;
+      RowVectorDIMS c_left;
+      Scalar sqr_d_left = m_left->squared_distance(V,Ele,p,sqr_d,i_left,c_left);
+      set_min(p,sqr_d_left,i_left,c_left,sqr_d,i,c);
+      looked_left = true;
+    };
+    const auto & look_right = [&]()
+    {
+      int i_right;
+      RowVectorDIMS c_right;
+      Scalar sqr_d_right = m_right->squared_distance(V,Ele,p,sqr_d,i_right,c_right);
+      set_min(p,sqr_d_right,i_right,c_right,sqr_d,i,c);
+      looked_right = true;
+    };
+
+    // must look left or right if in box
+    if(m_left->m_box.contains(p.transpose()))
+    {
+      look_left();
+    }
+    if(m_right->m_box.contains(p.transpose()))
+    {
+      look_right();
+    }
+    // if haven't looked left and could be less than current min, then look
+    Scalar  left_min_sqr_d = m_left->m_box.squaredExteriorDistance(p.transpose());
+    Scalar right_min_sqr_d = m_right->m_box.squaredExteriorDistance(p.transpose());
+    if(left_min_sqr_d < right_min_sqr_d)
+    {
+      if(!looked_left && left_min_sqr_d<sqr_d)
+      {
+        look_left();
+      }
+      if( !looked_right && right_min_sqr_d<sqr_d)
+      {
+        look_right();
+      }
+    }else
+    {
+      if( !looked_right && right_min_sqr_d<sqr_d)
+      {
+        look_right();
+      }
+      if(!looked_left && left_min_sqr_d<sqr_d)
+      {
+        look_left();
+      }
+    }
+  }
+  return sqr_d;
+}
+
+template <typename DerivedV, int DIM>
+template <
+  typename DerivedP, 
+  typename DerivedsqrD, 
+  typename DerivedI, 
+  typename DerivedC>
+inline void igl::AABB<DerivedV,DIM>::squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C) const
+{
+  assert(P.cols() == V.cols() && "cols in P should match dim of cols in V");
+  sqrD.resize(P.rows(),1);
+  I.resize(P.rows(),1);
+  C.resize(P.rows(),P.cols());
+  for(int p = 0;p<P.rows();p++)
+  {
+    RowVectorDIMS Pp = P.row(p), c;
+    int Ip;
+    sqrD(p) = squared_distance(V,Ele,Pp,Ip,c);
+    I(p) = Ip;
+    C.row(p) = c;
+  }
+}
+
+template <typename DerivedV, int DIM>
+inline void igl::AABB<DerivedV,DIM>::leaf_squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & p,
+  Scalar & sqr_d,
+  int & i,
+  RowVectorDIMS & c) const
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+
+  // Simplex size
+  const size_t ss = Ele.cols();
+  // Only one element per node
+  // plane unit normal
+  bool inside_triangle = false;
+  Scalar d_j = std::numeric_limits<Scalar>::infinity();
+  RowVectorDIMS pp;
+  // Only consider triangles, and non-degenerate triangles at that
+  if(ss == 3 && 
+      Ele(m_primitive,0) != Ele(m_primitive,1) && 
+      Ele(m_primitive,1) != Ele(m_primitive,2) && 
+      Ele(m_primitive,2) != Ele(m_primitive,0))
+  {
+    const RowVectorDIMS v10 = (V.row(Ele(m_primitive,1))- V.row(Ele(m_primitive,0)));
+    const RowVectorDIMS v20 = (V.row(Ele(m_primitive,2))- V.row(Ele(m_primitive,0)));
+    const RowVectorDIMS n = v10.cross(v20);
+    Scalar n_norm = n.norm();
+    if(n_norm > 0)
+    {
+      const RowVectorDIMS un = n/n.norm();
+      // vector to plane
+      const RowVectorDIMS bc = 
+        1./3.*
+        ( V.row(Ele(m_primitive,0))+
+          V.row(Ele(m_primitive,1))+
+          V.row(Ele(m_primitive,2)));
+      const auto & v = p-bc;
+      // projected point on plane
+      d_j = v.dot(un);
+      pp = p - d_j*un;
+      // determine if pp is inside triangle
+      Eigen::Matrix<Scalar,1,3> b;
+      barycentric_coordinates(
+            pp,
+            V.row(Ele(m_primitive,0)),
+            V.row(Ele(m_primitive,1)),
+            V.row(Ele(m_primitive,2)),
+            b);
+      inside_triangle = fabs(fabs(b(0)) + fabs(b(1)) + fabs(b(2)) - 1.) <= 1e-10;
+    }
+  }
+  const auto & point_point_squared_distance = [&](const RowVectorDIMS & s)
+  {
+    const Scalar sqr_d_s = (p-s).squaredNorm();
+    set_min(p,sqr_d_s,m_primitive,s,sqr_d,i,c);
+  };
+  if(inside_triangle)
+  {
+    // point-triangle squared distance
+    const Scalar sqr_d_j = d_j*d_j;
+    //cout<<"point-triangle..."<<endl;
+    set_min(p,sqr_d_j,m_primitive,pp,sqr_d,i,c);
+  }else
+  {
+    if(ss >= 2)
+    {
+      // point-segment distance
+      // number of edges
+      size_t ne = ss==3?3:1;
+      for(size_t x = 0;x<ne;x++)
+      {
+        const size_t e1 = Ele(m_primitive,(x+1)%ss);
+        const size_t e2 = Ele(m_primitive,(x+2)%ss);
+        const RowVectorDIMS & s = V.row(e1);
+        const RowVectorDIMS & d = V.row(e2);
+        // Degenerate edge
+        if(e1 == e2 || (s-d).squaredNorm()==0)
+        {
+          // only consider once
+          if(e1 < e2)
+          {
+            point_point_squared_distance(s);
+          }
+          continue;
+        }
+        Matrix<Scalar,1,1> sqr_d_j_x(1,1);
+        Matrix<Scalar,1,1> t(1,1);
+        project_to_line_segment(p,s,d,t,sqr_d_j_x);
+        const RowVectorDIMS q = s+t(0)*(d-s);
+        set_min(p,sqr_d_j_x(0),m_primitive,q,sqr_d,i,c);
+      }
+    }else
+    {
+      // So then Ele is just a list of points...
+      assert(ss == 1);
+      const RowVectorDIMS & s = V.row(Ele(m_primitive,0));
+      point_point_squared_distance(s);
+    }
+  }
+}
+
+
+template <typename DerivedV, int DIM>
+inline void igl::AABB<DerivedV,DIM>::set_min(
+  const RowVectorDIMS & p,
+  const Scalar sqr_d_candidate,
+  const int i_candidate,
+  const RowVectorDIMS & c_candidate,
+  Scalar & sqr_d,
+  int & i,
+  RowVectorDIMS & c) const
+{
+#ifndef NDEBUG
+  const Scalar diff = fabs(sqr_d_candidate - (p-c_candidate).squaredNorm());
+  assert(diff<=1e-10 && "distance should match norm of difference");
+#endif
+  if(sqr_d_candidate < sqr_d)
+  {
+    i = i_candidate;
+    c = c_candidate;
+    sqr_d = sqr_d_candidate;
+  }
+}
+
+
+template <typename DerivedV, int DIM>
+template <int SS>
+inline void
+igl::AABB<DerivedV,DIM>::barycentric_coordinates(
+  const RowVectorDIMS & p, 
+  const RowVectorDIMS & a, 
+  const RowVectorDIMS & b, 
+  const RowVectorDIMS & c,
+  Eigen::Matrix<Scalar,1,SS> & bary)
+{
+  assert(SS==3);
+  // http://gamedev.stackexchange.com/a/23745
+  const RowVectorDIMS v0 = b - a;
+  const RowVectorDIMS v1 = c - a;
+  const RowVectorDIMS v2 = p - a;
+  Scalar d00 = v0.dot(v0);
+  Scalar d01 = v0.dot(v1);
+  Scalar d11 = v1.dot(v1);
+  Scalar d20 = v2.dot(v0);
+  Scalar d21 = v2.dot(v1);
+  Scalar denom = d00 * d11 - d01 * d01;
+  bary(1) = (d11 * d20 - d01 * d21) / denom;
+  bary(2) = (d00 * d21 - d01 * d20) / denom;
+  bary(0) = 1.0f - bary(1) - bary(2);
+};
+
+#endif

include/igl/C_STR.h

@@ -14,5 +14,5 @@
 //   func(C_STR("foo"<<1<<"bar"));
 #include <sstream>
 #include <string>
-#define C_STR(X) static_cast<std::ostringstream&>(std::ostringstream().seekp(0) << X).str().c_str()
+#define C_STR(X) static_cast<std::ostringstream&>(std::ostringstream().flush() << X).str().c_str()
 #endif

include/igl/Camera.h

@@ -35,21 +35,16 @@ namespace igl
       //  m_near  near clipping plane {1e-2}
       //  m_far  far clipping plane {100}
       //  m_at_dist  distance of looking at point {1}
+      //  m_orthographic  whether to use othrographic projection {false}
       //  m_rotation_conj  Conjugate of rotation part of rigid transformation of
       //    camera {identity}. Note: we purposefully store the conjugate because
       //    this is what TW_TYPE_QUAT4D is expecting.
       //  m_translation  Translation part of rigid transformation of camera
       //    {(0,0,1)}
       double m_angle, m_aspect, m_near, m_far, m_at_dist;
+      bool m_orthographic;
       Eigen::Quaterniond m_rotation_conj;
       Eigen::Vector3d m_translation;
-    private:
-      // m_at_dist_min_angle  m_at_dist from last time m_angle set to <= IGL_CAMERA_MIN_ANGLE
-      double m_at_dist_min_angle;
-      double m_angle_min_angle;
-    //  // m_last_positive_m_angle
-    //  // m_last_positive_m_angle_m_at_dist
-    //  double m_last_positive_m_angle,m_last_positive_m_angle_m_at_dist;
     public:
       inline Camera();
       inline virtual ~Camera(){}
@@ -154,10 +149,9 @@ namespace igl
 
 inline igl::Camera::Camera():
   m_angle(45.0),m_aspect(1),m_near(1e-2),m_far(100),m_at_dist(1),
+  m_orthographic(false),
   m_rotation_conj(1,0,0,0),
-  m_translation(0,0,1),
-  m_at_dist_min_angle(m_at_dist),
-  m_angle_min_angle(m_angle)
+  m_translation(0,0,1)
 {
 }
 
@@ -165,41 +159,36 @@ inline Eigen::Matrix4d igl::Camera::projection() const
 {
   Eigen::Matrix4d P;
   using namespace std;
-  using namespace igl;
+  const double far = m_at_dist + m_far;
+  const double near = m_near;
   // http://stackoverflow.com/a/3738696/148668
-  if(m_angle >= IGL_CAMERA_MIN_ANGLE)
-  {
-    const double yScale = tan(PI*0.5 - 0.5*m_angle*PI/180.);
-    // http://stackoverflow.com/a/14975139/148668
-    const double xScale = yScale/m_aspect;
-    const double far = m_at_dist + m_far;
-    const double near = m_near;
-    P<< 
-      xScale, 0, 0, 0,
-      0, yScale, 0, 0,
-      0, 0, -(far+near)/(far-near), -1,
-      0, 0, -2.*near*far/(far-near), 0;
-    P = P.transpose().eval();
-  }else
+  if(m_orthographic)
   {
     const double f = 0.5;
     const double left = -f*m_aspect;
     const double right = f*m_aspect;
     const double bottom = -f;
     const double top = f;
-    const double near = m_near;
-    const double far = m_at_dist + m_far;
     const double tx = (right+left)/(right-left);
     const double ty = (top+bottom)/(top-bottom);
     const double tz = (far+near)/(far-near);
-    const double z_fix = 
-      0.5/(m_at_dist_min_angle * tan(m_angle_min_angle/2./180.*M_PI))+
-      (-m_at_dist+m_at_dist_min_angle)/m_at_dist_min_angle;
+    const double z_fix = 0.5 /m_at_dist / tan(m_angle*0.5 * (M_PI/180.) );
     P<<
       z_fix*2./(right-left), 0, 0, -tx,
       0, z_fix*2./(top-bottom), 0, -ty,
       0, 0, -z_fix*2./(far-near),  -tz,
       0, 0, 0, 1;
+  }else
+  {
+    const double yScale = tan(PI*0.5 - 0.5*m_angle*PI/180.);
+    // http://stackoverflow.com/a/14975139/148668
+    const double xScale = yScale/m_aspect;
+    P<< 
+      xScale, 0, 0, 0,
+      0, yScale, 0, 0,
+      0, 0, -(far+near)/(far-near), -1,
+      0, 0, -2.*near*far/(far-near), 0;
+    P = P.transpose().eval();
   }
   return P;
 }
@@ -244,7 +233,6 @@ inline Eigen::Vector3d igl::Camera::up() const
 
 inline Eigen::Vector3d igl::Camera::unit_plane() const
 {
-  using namespace igl;
   // Distance of center pixel to eye
   const double d = 1.0;
   const double a = m_aspect;
@@ -263,7 +251,6 @@ inline void igl::Camera::dolly(const Eigen::Vector3d & dv)
 inline void igl::Camera::push_away(const double s)
 {
   using namespace Eigen;
-  using namespace igl;
 #ifndef NDEBUG
   Vector3d old_at = at();
 #endif
@@ -276,47 +263,36 @@ inline void igl::Camera::push_away(const double s)
 inline void igl::Camera::dolly_zoom(const double da)
 {
   using namespace std;
-  using namespace igl;
   using namespace Eigen;
 #ifndef NDEBUG
   Vector3d old_at = at();
 #endif
   const double old_angle = m_angle;
-  m_angle += da;
-  m_angle = min(89.,max(0.,m_angle));
-  const double eff_angle = (IGL_CAMERA_MIN_ANGLE > m_angle ? IGL_CAMERA_MIN_ANGLE : m_angle);
-  if(old_angle >= IGL_CAMERA_MIN_ANGLE)
+  if(old_angle + da < IGL_CAMERA_MIN_ANGLE)
+  {
+    m_orthographic = true;
+  }else if(old_angle + da > IGL_CAMERA_MIN_ANGLE)
   {
+    m_orthographic = false;
+  }
+  if(!m_orthographic)
+  {
+    m_angle += da;
+    m_angle = min(89.,max(IGL_CAMERA_MIN_ANGLE,m_angle));
     // change in distance
     const double s = 
       (2.*tan(old_angle/2./180.*M_PI)) /
-      (2.*tan(eff_angle/2./180.*M_PI)) ;
+      (2.*tan(m_angle/2./180.*M_PI)) ;
     const double old_at_dist = m_at_dist;
     m_at_dist = old_at_dist * s;
     dolly(Vector3d(0,0,1)*(m_at_dist - old_at_dist));
-    if(eff_angle == IGL_CAMERA_MIN_ANGLE)
-    {
-      m_at_dist_min_angle = m_at_dist;
-      m_angle_min_angle = eff_angle;
-    }
     assert((old_at-at()).squaredNorm() < DOUBLE_EPS);
-  }else if(old_angle < IGL_CAMERA_MIN_ANGLE && m_angle >= IGL_CAMERA_MIN_ANGLE)
-  {
-    // Restore decent length
-    const double z_fix = 
-      // There should be some factor here based on the incoming angle
-      // (m_angle_min_angle) and outgoing angle (m_angle)... For now I set it
-      // to 1. (assumes equality)
-      //0.5/(m_at_dist_min_angle * tan(m_angle_min_angle/2./180.*M_PI))+
-      1.+(-m_at_dist+m_at_dist_min_angle)/m_at_dist_min_angle;
-    m_at_dist = m_at_dist_min_angle / z_fix;
   }
 }
 
 inline void igl::Camera::turn_eye(const Eigen::Quaterniond & q)
 {
   using namespace Eigen;
-  using namespace igl;
   Vector3d old_eye = eye();
   // eye should be fixed
   //
@@ -330,7 +306,6 @@ inline void igl::Camera::turn_eye(const Eigen::Quaterniond & q)
 inline void igl::Camera::orbit(const Eigen::Quaterniond & q)
 {
   using namespace Eigen;
-  using namespace igl;
   Vector3d old_at = at();
   // at should be fixed
   //
@@ -351,7 +326,6 @@ inline void igl::Camera::look_at(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   // http://www.opengl.org/sdk/docs/man2/xhtml/gluLookAt.xml
   // Normalize vector from at to eye
   Vector3d F = eye-at;

include/igl/InElementAABB.h

@@ -1,407 +0,0 @@
-#ifndef IGL_IN_ELEMENT_AABB_H
-#define IGL_IN_ELEMENT_AABB_H
-
-#include <Eigen/Core>
-#include <memory>
-#include <vector>
-namespace igl
-{
-  class InElementAABB
-  {
-    public:
-      std::shared_ptr<InElementAABB> m_left, m_right;
-      Eigen::RowVectorXd m_bb_min,m_bb_max;
-      // -1 non-leaf
-      int m_element;
-      InElementAABB():
-        m_left(NULL), m_right(NULL),
-        m_bb_min(), m_bb_max(), m_element(-1)
-      {}
-      // Build an Axis-Aligned Bounding Box tree for a given mesh and given
-      // serialization of a previous AABB tree.
-      //
-      // Inputs:
-      //   V  #V by dim list of mesh vertex positions. 
-      //   Ele  #Ele by dim+1 list of mesh indices into #V. 
-      //   bb_mins  max_tree by dim list of bounding box min corner positions
-      //   bb_maxs  max_tree by dim list of bounding box max corner positions
-      //   elements  max_tree list of element or (not leaf id) indices into Ele
-      //   i  recursive call index {0}
-      inline void init(
-          const Eigen::MatrixXd & V,
-          const Eigen::MatrixXi & Ele, 
-          const Eigen::MatrixXd & bb_mins,
-          const Eigen::MatrixXd & bb_maxs,
-          const Eigen::VectorXi & elements,
-          const int i = 0);
-      // Wrapper for root with empty serialization
-      inline void init(
-        const Eigen::MatrixXd & V,
-        const Eigen::MatrixXi & Ele);
-      // Build an Axis-Aligned Bounding Box tree for a given mesh.
-      //
-      // Inputs:
-      //   V  #V by dim list of mesh vertex positions. 
-      //   Ele  #Ele by dim+1 list of mesh indices into #V. 
-      //   SI  #Ele by dim list revealing for each coordinate where Ele's
-      //     barycenters would be sorted: SI(e,d) = i --> the dth coordinate of
-      //     the barycenter of the eth element would be placed at position i in a
-      //     sorted list.
-      //   I  #I list of indices into Ele of elements to include (for recursive
-      //     calls)
-      // 
-      inline void init(
-          const Eigen::MatrixXd & V,
-          const Eigen::MatrixXi & Ele, 
-          const Eigen::MatrixXi & SI,
-          const Eigen::VectorXi & I);
-      // Find the indices of elements containing given point.
-      //
-      // Inputs:
-      //   V  #V by dim list of mesh vertex positions. **Should be same as used to
-      //     construct mesh.**
-      //   Ele  #Ele by dim+1 list of mesh indices into #V. **Should be same as used to
-      //     construct mesh.**
-      //   q  dim row-vector query position
-      //   first  whether to only return first element containing q
-      // Returns:
-      //   list of indices of elements containing q
-      inline std::vector<int> find(
-        const Eigen::MatrixXd & V,
-        const Eigen::MatrixXi & Ele, 
-        const Eigen::RowVectorXd & q, 
-        const bool first=false) const;
-  
-      // If number of elements m then total tree size should be 2*h where h is
-      // the deepest depth 2^ceil(log(#Ele*2-1))
-      inline int subtree_size();
-  
-      // Serialize this class into 3 arrays (so we can pass it pack to matlab)
-      //
-      // Outputs:
-      //   bb_mins  max_tree by dim list of bounding box min corner positions
-      //   bb_maxs  max_tree by dim list of bounding box max corner positions
-      //   elements  max_tree list of element or (not leaf id) indices into Ele
-      //   i  recursive call index into these arrays {0}
-      inline void serialize(
-        Eigen::MatrixXd & bb_mins,
-        Eigen::MatrixXd & bb_maxs,
-        Eigen::VectorXi & elements,
-        const int i = 0);
-  };
-}
-
-// Implementation
-#include <igl/volume.h>
-#include <igl/colon.h>
-#include <igl/doublearea.h>
-#include <igl/matlab_format.h>
-#include <igl/colon.h>
-#include <igl/sort.h>
-#include <igl/barycenter.h>
-#include <iostream>
-
-inline void igl::InElementAABB::init(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & Ele, 
-    const Eigen::MatrixXd & bb_mins,
-    const Eigen::MatrixXd & bb_maxs,
-    const Eigen::VectorXi & elements,
-    const int i)
-{
-  using namespace std;
-  using namespace Eigen;
-  using namespace igl;
-  if(bb_mins.size() > 0)
-  {
-    assert(bb_mins.rows() == bb_maxs.rows() && "Serial tree arrays must match");
-    assert(bb_mins.cols() == V.cols() && "Serial tree array dim must match V");
-    assert(bb_mins.cols() == bb_maxs.cols() && "Serial tree arrays must match");
-    assert(bb_mins.rows() == elements.rows() &&
-      "Serial tree arrays must match");
-    // construct from serialization
-    m_bb_min = bb_mins.row(i);
-    m_bb_max = bb_maxs.row(i);
-    m_element = elements(i);
-    // Not leaf then recurse
-    if(m_element == -1)
-    {
-      m_left = make_shared<InElementAABB>();
-      m_left->init( V,Ele,bb_mins,bb_maxs,elements,2*i+1);
-      m_right = make_shared<InElementAABB>();
-      m_right->init( V,Ele,bb_mins,bb_maxs,elements,2*i+2);
-    }
-  }else
-  {
-    VectorXi allI = colon<int>(0,Ele.rows()-1);
-    MatrixXd BC;
-    barycenter(V,Ele,BC);
-    MatrixXi SI(BC.rows(),BC.cols());
-    {
-      MatrixXd _;
-      MatrixXi IS;
-      igl::sort(BC,1,true,_,IS);
-      // Need SI(i) to tell which place i would be sorted into
-      const int dim = IS.cols();
-      for(int i = 0;i<IS.rows();i++)
-      {
-        for(int d = 0;d<dim;d++)
-        {
-          SI(IS(i,d),d) = i;
-        }
-      }
-    }
-    init(V,Ele,SI,allI);
-  }
-}
-
-inline void igl::InElementAABB::init(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & Ele)
-{
-  using namespace Eigen;
-  return init(V,Ele,MatrixXd(),MatrixXd(),VectorXi(),0);
-}
-
-inline void igl::InElementAABB::init(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & Ele, 
-    const Eigen::MatrixXi & SI,
-    const Eigen::VectorXi & I)
-{
-  using namespace Eigen;
-  using namespace std;
-  using namespace igl;
-  const int dim = V.cols();
-  const double inf = numeric_limits<double>::infinity();
-  m_bb_min.setConstant(1,dim, inf);
-  m_bb_max.setConstant(1,dim,-inf);
-  // Compute bounding box
-  for(int i = 0;i<I.rows();i++)
-  {
-    for(int c = 0;c<Ele.cols();c++)
-    {
-      for(int d = 0;d<dim;d++)
-      {
-        m_bb_min(d) = min(m_bb_min(d),V(Ele(I(i),c),d));
-        m_bb_max(d) = max(m_bb_max(d),V(Ele(I(i),c),d));
-      }
-    }
-  }
-  switch(I.size())
-  {
-    case 0:
-      {
-        assert(false);
-      }
-    case 1:
-      {
-        m_element = I(0);
-        break;
-      }
-    default:
-      {
-        // Compute longest direction
-        int max_d = -1;
-        double max_len = -inf;
-        for(int d = 0;d<dim;d++)
-        {
-          const auto diff = (m_bb_max[d] - m_bb_min[d]);
-          if( diff > max_len )
-          {
-            max_len = diff;
-            max_d = d;
-          }
-        }
-        // Can't use median on BC directly because many may have same value,
-        // but can use median on sorted BC indices
-        VectorXi SIdI(I.rows());
-        for(int i = 0;i<I.rows();i++)
-        {
-          SIdI(i) = SI(I(i),max_d);
-        }
-        // Since later I use <= I think I don't need to worry about odd/even
-        // Pass by copy to avoid changing input
-        const auto median = [](VectorXi A)->double
-        {
-          size_t n = A.size()/2;
-          nth_element(A.data(),A.data()+n,A.data()+A.size());
-          if(A.rows() % 2 == 1)
-          {
-            return A(n);
-          }else
-          {
-            nth_element(A.data(),A.data()+n-1,A.data()+A.size());
-            return 0.5*(A(n)+A(n-1));
-          }
-        };
-        const double med = median(SIdI);
-        VectorXi LI((I.rows()+1)/2),RI(I.rows()/2);
-        assert(LI.rows()+RI.rows() == I.rows());
-        // Distribute left and right
-        {
-          int li = 0;
-          int ri = 0;
-          for(int i = 0;i<I.rows();i++)
-          {
-            if(SIdI(i)<=med)
-            {
-              LI(li++) = I(i);
-            }else
-            {
-              RI(ri++) = I(i);
-            }
-          }
-        }
-        if(LI.rows()>0)
-        {
-          m_left = make_shared<InElementAABB>();
-          m_left->init(V,Ele,SI,LI);
-        }
-        if(RI.rows()>0)
-        {
-          m_right = make_shared<InElementAABB>();
-          m_right->init(V,Ele,SI,RI);
-        }
-      }
-  }
-}
-
-inline std::vector<int> igl::InElementAABB::find(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & Ele, 
-    const Eigen::RowVectorXd & q, 
-    const bool first) const
-{
-  using namespace std;
-  using namespace igl;
-  using namespace Eigen;
-  bool inside = true;
-  const int dim = m_bb_max.size();
-  assert(q.size() == m_bb_max.size());
-  const double epsilon = 1e-14;
-  // Check if outside bounding box
-  for(int d = 0;d<q.size()&&inside;d++)
-  {
-    inside &= (q(d)-m_bb_min(d))>=epsilon;
-    inside &= (m_bb_max(d)-q(d))>=epsilon;
-  }
-  cout<<"searching..."<<endl;
-  if(!inside)
-  {
-    cout<<"not in bb"<<endl;
-    return std::vector<int>();
-  }
-  if(m_element != -1)
-  {
-    // Initialize to some value > -epsilon
-    double a1=1,a2=1,a3=1,a4=1;
-    switch(dim)
-    {
-      case 3:
-        {
-          // Barycentric coordinates
-          const RowVector3d V1 = V.row(Ele(m_element,0));
-          const RowVector3d V2 = V.row(Ele(m_element,1));
-          const RowVector3d V3 = V.row(Ele(m_element,2));
-          const RowVector3d V4 = V.row(Ele(m_element,3));
-          a1 = volume_single(V2,V4,V3,(RowVector3d)q);
-          a2 = volume_single(V1,V3,V4,(RowVector3d)q);
-          a3 = volume_single(V1,V4,V2,(RowVector3d)q);
-          a4 = volume_single(V1,V2,V3,(RowVector3d)q);
-          break;
-        }
-      case 2:
-        {
-          // Barycentric coordinates
-          const Vector2d V1 = V.row(Ele(m_element,0));
-          const Vector2d V2 = V.row(Ele(m_element,1));
-          const Vector2d V3 = V.row(Ele(m_element,2));
-          double a0 = doublearea_single(V1,V2,V3);
-          a1 = doublearea_single(V1,V2,(Vector2d)q);
-          a2 = doublearea_single(V2,V3,(Vector2d)q);
-          a3 = doublearea_single(V3,V1,(Vector2d)q);
-          cout<<
-            a0<<" "<<
-            a1<<" "<<
-            a2<<" "<<
-            a3<<" "<<
-            endl;
-          break;
-        }
-      default:assert(false);
-    }
-    if(
-        a1>=-epsilon && 
-        a2>=-epsilon && 
-        a3>=-epsilon && 
-        a4>=-epsilon)
-    {
-      return std::vector<int>(1,m_element);
-    }else
-    {
-      return std::vector<int>();
-    }
-  }
-  std::vector<int> left = m_left->find(V,Ele,q,first);
-  if(first && !left.empty())
-  {
-    return left;
-  }
-  std::vector<int> right = m_right->find(V,Ele,q,first);
-  if(first)
-  {
-    return right;
-  }
-  left.insert(left.end(),right.begin(),right.end());
-  return left;
-}
-
-inline int igl::InElementAABB::subtree_size()
-{
-  // 1 for self
-  int n = 1;
-  int n_left = 0,n_right = 0;
-  if(m_left != NULL)
-  {
-    n_left = m_left->subtree_size();
-  }
-  if(m_right != NULL)
-  {
-    n_right = m_right->subtree_size();
-  }
-  n += 2*std::max(n_left,n_right);
-  return n;
-}
-
-
-inline void igl::InElementAABB::serialize(
-    Eigen::MatrixXd & bb_mins,
-    Eigen::MatrixXd & bb_maxs,
-    Eigen::VectorXi & elements,
-    const int i)
-{
-  using namespace std;
-  // Calling for root then resize output
-  if(i==0)
-  {
-    const int m = subtree_size();
-    //cout<<"m: "<<m<<endl;
-    bb_mins.resize(m,m_bb_min.size());
-    bb_maxs.resize(m,m_bb_max.size());
-    elements.resize(m,1);
-  }
-  //cout<<i<<" ";
-  bb_mins.row(i) = m_bb_min;
-  bb_maxs.row(i) = m_bb_max;
-  elements(i) = m_element;
-  if(m_left != NULL)
-  {
-    m_left->serialize(bb_mins,bb_maxs,elements,2*i+1);
-  }
-  if(m_right != NULL)
-  {
-    m_right->serialize(bb_mins,bb_maxs,elements,2*i+2);
-  }
-}
-#endif

include/igl/MouseController.h

@@ -214,7 +214,6 @@ inline void igl::MouseController::reshape(const int w, const int h)
 inline bool igl::MouseController::down(const int x, const int y)
 {
   using namespace std;
-  using namespace igl;
   m_down_x = m_drag_x =x;
   m_down_y = m_drag_y =y;
   const bool widget_down = any_selection() && m_widget.down(x,m_height-y);

include/igl/OpenGL_convenience.h

@@ -18,6 +18,7 @@
 #if __APPLE__
 #  include <OpenGL/gl.h>
 #  include <OpenGL/glu.h>
+#  include <OpenGL/glext.h>
 #elif defined(_WIN32)
 #    define NOMINMAX
 #    include <Windows.h>

include/igl/RotateWidget.h

@@ -118,7 +118,6 @@ inline Eigen::Quaterniond igl::RotateWidget::axis_q(const int a)
 inline Eigen::Vector3d igl::RotateWidget::view_direction(const int x, const int y)
 {
   using namespace Eigen;
-  using namespace igl;
   const Vector3d win_s(x,y,0), win_d(x,y,1);
   const Vector3d s = unproject(win_s);
   const Vector3d d = unproject(win_d);
@@ -128,7 +127,6 @@ inline Eigen::Vector3d igl::RotateWidget::view_direction(const int x, const int
 inline Eigen::Vector3d igl::RotateWidget::view_direction(const Eigen::Vector3d & pos)
 {
   using namespace Eigen;
-  using namespace igl;
   const Vector3d ppos = project(pos);
   return view_direction(ppos(0),ppos(1));
 }
@@ -151,7 +149,6 @@ inline Eigen::Vector3d igl::RotateWidget::unproject_onto(
   const int y) const
 {
   using namespace Eigen;
-  using namespace igl;
   // KNOWN BUG: This projects to same depths as pos. I think what we actually
   // want is The intersection with the plane perpendicular to the view
   // direction at pos. If the field of view angle is small then this difference
@@ -179,7 +176,6 @@ inline bool igl::RotateWidget::intersect(
   Eigen::Vector3d & hit) const
 {
   using namespace Eigen;
-  using namespace igl;
   Vector3d view = view_direction(x,y);
   const Vector3d ppos = project(pos);
   Vector3d uxy = unproject(Vector3d(x,y,ppos(2)));
@@ -196,7 +192,6 @@ inline bool igl::RotateWidget::intersect(
 inline double igl::RotateWidget::unprojected_inner_radius() const
 {
   using namespace Eigen;
-  using namespace igl;
   Vector3d off,ppos,ppos_off,pos_off;
   project(pos,ppos);
   ppos_off = ppos;
@@ -207,7 +202,6 @@ inline double igl::RotateWidget::unprojected_inner_radius() const
 inline bool igl::RotateWidget::down(const int x, const int y)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   if(!m_is_enabled)
   {
@@ -289,7 +283,6 @@ inline bool igl::RotateWidget::down(const int x, const int y)
 
 inline bool igl::RotateWidget::drag(const int x, const int y)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   if(!m_is_enabled)
@@ -378,7 +371,6 @@ inline void igl::RotateWidget::draw() const
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   glPushAttrib(GL_ENABLE_BIT | GL_LIGHTING_BIT | GL_DEPTH_BUFFER_BIT | GL_LINE_BIT);
   glDisable(GL_CLIP_PLANE0);
 
@@ -474,7 +466,6 @@ inline void igl::RotateWidget::draw_guide() const
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   glPushAttrib(
     GL_DEPTH_BUFFER_BIT | 
     GL_ENABLE_BIT | 

include/igl/STR.h

@@ -14,5 +14,5 @@
 //   void func(std::string c);
 // Then you can write:
 //   func(STR("foo"<<1<<"bar"));
-#define STR(X) static_cast<std::ostringstream&>(std::ostringstream().seekp(0) << X).str()
+#define STR(X) static_cast<std::ostringstream&>(std::ostringstream().flush() << X).str()
 #endif 

include/igl/WindingNumberAABB.h

@@ -28,10 +28,13 @@ namespace igl
       {
         CENTER_ON_LONGEST_AXIS = 0,
         MEDIAN_ON_LONGEST_AXIS = 1,
-        NUM_SPLIT_METHODS = 3
+        NUM_SPLIT_METHODS = 2
       } split_method;
     public:
-      inline WindingNumberAABB(){}
+      inline WindingNumberAABB():
+        total_positive_area(std::numeric_limits<double>::infinity()),
+        split_method(MEDIAN_ON_LONGEST_AXIS)
+      {}
       inline WindingNumberAABB(
         const Eigen::MatrixXd & V,
         const Eigen::MatrixXi & F);
@@ -83,7 +86,6 @@ inline void igl::WindingNumberAABB<Point>::set_mesh(
 template <typename Point>
 inline void igl::WindingNumberAABB<Point>::init()
 {
-  using namespace igl;
   using namespace Eigen;
   assert(max_corner.size() == 3);
   assert(min_corner.size() == 3);
@@ -124,7 +126,6 @@ inline void igl::WindingNumberAABB<Point>::grow()
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   //cout<<"cap.rows(): "<<this->getcap().rows()<<endl;
   //cout<<"F.rows(): "<<this->getF().rows()<<endl;
 
@@ -300,7 +301,6 @@ inline double igl::WindingNumberAABB<Point>::max_simple_abs_winding_number(const
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   // Only valid if not inside
   if(inside(p))
   {

include/igl/WindingNumberTree.h

@@ -45,7 +45,7 @@ namespace igl
       // (Approximate) center (of mass)
       Point center;
     public:
-      inline WindingNumberTree():V(dummyV){}
+      inline WindingNumberTree();
       // For root
       inline WindingNumberTree(
         const Eigen::MatrixXd & V,
@@ -142,6 +142,20 @@ template <typename Point>
 std::map< std::pair<const igl::WindingNumberTree<Point>*,const igl::WindingNumberTree<Point>*>, double>
   igl::WindingNumberTree<Point>::cached;
 
+template <typename Point>
+inline igl::WindingNumberTree<Point>::WindingNumberTree():
+  method(EXACT_WINDING_NUMBER_METHOD),
+  parent(NULL),
+  V(dummyV),
+  SV(),
+  F(),
+  //boundary(igl::boundary_facets<Eigen::MatrixXi,Eigen::MatrixXi>(F))
+  cap(),
+  radius(std::numeric_limits<double>::infinity()),
+  center(0,0,0)
+{
+}
+
 template <typename Point>
 inline igl::WindingNumberTree<Point>::WindingNumberTree(
   const Eigen::MatrixXd & _V,
@@ -406,7 +420,6 @@ inline double igl::WindingNumberTree<Point>::cached_winding_number(
   const Point & p) const
 {
   using namespace std;
-  using namespace igl;
   // Simple metric for "far".
   //   this             that
   //                   --------

include/igl/active_set.cpp

@@ -48,7 +48,6 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 #ifdef ACTIVE_SET_CPP_DEBUG
 #  warning "ACTIVE_SET_CPP_DEBUG"
 #endif
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
   SolverStatus ret = SOLVER_STATUS_ERROR;

include/igl/angles.h

@@ -8,13 +8,14 @@
 #ifndef IGL_ANGLES_H
 #define IGL_ANGLES_H
 #ifdef _WIN32
- #pragma message ( "Deprecated. Use igl/internal_angles.h instead" )
+#  pragma message ( "Deprecated. Use igl/internal_angles.h instead" )
 #else
- #warning "Deprecated. Use igl/internal_angles.h instead"
+#  warning "Deprecated. Use igl/internal_angles.h instead"
 #endif
 
 
 #include "igl_inline.h"
+#include "deprecated.h"
 #include <Eigen/Core>
 namespace igl
 {
@@ -30,10 +31,12 @@ namespace igl
   typename DerivedV,
   typename DerivedF,
   typename Derivedtheta>
-  IGL_INLINE void angles(
-  const Eigen::PlainObjectBase<DerivedV>& V,
-  const Eigen::PlainObjectBase<DerivedF>& F,
-  Eigen::PlainObjectBase<Derivedtheta>& theta);
+  IGL_INLINE 
+  IGL_DEPRECATED(
+    void angles(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<Derivedtheta>& theta));
 
 }
 

include/igl/angular_distance.cpp

@@ -12,7 +12,6 @@ IGL_INLINE double igl::angular_distance(
   const Eigen::Quaterniond & A,
   const Eigen::Quaterniond & B)
 {
-  using namespace igl;
   assert(fabs(A.norm()-1)<FLOAT_EPS && "A should be unit norm");
   assert(fabs(B.norm()-1)<FLOAT_EPS && "B should be unit norm");
   //// acos is always in [0,2*pi)

include/igl/arap_linear_block.cpp

@@ -18,7 +18,6 @@ IGL_INLINE void igl::arap_linear_block(
   const igl::ARAPEnergyType energy,
   Eigen::SparseMatrix<Scalar> & Kd)
 {
-  using namespace igl;
   switch(energy)
   {
     case ARAP_ENERGY_TYPE_SPOKES:
@@ -44,7 +43,6 @@ IGL_INLINE void igl::arap_linear_block_spokes(
   const int d,
   Eigen::SparseMatrix<Scalar> & Kd)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // simplex size (3: triangles, 4: tetrahedra)
@@ -110,7 +108,6 @@ IGL_INLINE void igl::arap_linear_block_spokes_and_rims(
   const int d,
   Eigen::SparseMatrix<Scalar> & Kd)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // simplex size (3: triangles, 4: tetrahedra)
@@ -193,7 +190,6 @@ IGL_INLINE void igl::arap_linear_block_elements(
   const int d,
   Eigen::SparseMatrix<Scalar> & Kd)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // simplex size (3: triangles, 4: tetrahedra)

include/igl/arap_rhs.cpp

@@ -19,7 +19,6 @@ IGL_INLINE void igl::arap_rhs(
   const igl::ARAPEnergyType energy,
   Eigen::SparseMatrix<double>& K)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // Number of dimensions

include/igl/barycentric_coordinates.h

@@ -46,6 +46,8 @@ namespace igl
   // Outputs:
   //   L  #P by e list of barycentric coordinates
   //   
+  // Known bugs: this code is not tested (and probably will not work) for
+  // triangles and queries in 3D even if the query lives in/on the triangle.
   template <
     typename DerivedP,
     typename DerivedA,

include/igl/basename.h

@@ -13,7 +13,7 @@
 
 namespace igl
 {
-  // Function like PHP's basename
+  // Function like PHP's basename: /etc/sudoers.d --> sudoers.d
   // Input:
   //  path  string containing input path
   // Returns string containing basename (see php's basename)

include/igl/boolean/mesh_boolean.cpp

@@ -1,11 +1,13 @@
 #include "mesh_boolean.h"
-#include <igl/peel_outer_hull_layers.h>
 #include <igl/per_face_normals.h>
+#include <igl/cgal/peel_outer_hull_layers.h>
 #include <igl/cgal/remesh_self_intersections.h>
 #include <igl/remove_unreferenced.h>
 #include <igl/unique_simplices.h>
 #include <iostream>
 
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
 //#define IGL_MESH_BOOLEAN_DEBUG
 
 template <
@@ -30,8 +32,8 @@ IGL_INLINE void igl::mesh_boolean(
     const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
     const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
           Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&, 
-          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&)> 
+          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&)>
     empty_fun;
   return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J);
 }
@@ -57,8 +59,8 @@ IGL_INLINE void igl::mesh_boolean(
     const Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
     const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
           Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3> &,
-          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&, 
-          Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1>&)> 
+          Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
+          Eigen::Matrix<typename DerivedFC::Index, Eigen::Dynamic,1>&)>
     empty_fun;
   return mesh_boolean(VA,FA,VB,FB,type,empty_fun,VC,FC,J);
 }
@@ -82,7 +84,7 @@ IGL_INLINE void igl::mesh_boolean(
     const Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
           Eigen::Matrix<typename DerivedVC::Scalar,Eigen::Dynamic,3>&,
           Eigen::Matrix<typename DerivedFC::Scalar, Eigen::Dynamic,3>&,
-          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&)> 
+          Eigen::Matrix<typename DerivedJ::Scalar, Eigen::Dynamic,1>&)>
     & resolve_fun,
   Eigen::PlainObjectBase<DerivedVC > & VC,
   Eigen::PlainObjectBase<DerivedFC > & FC,
@@ -93,9 +95,12 @@ IGL_INLINE void igl::mesh_boolean(
   using namespace igl;
   MeshBooleanType eff_type = type;
   // Concatenate A and B into a single mesh
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+  typedef Kernel::FT ExactScalar;
   typedef typename DerivedVC::Scalar Scalar;
   typedef typename DerivedFC::Scalar Index;
   typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
+  typedef Matrix<ExactScalar,Dynamic,3> MatrixX3ES;
   typedef Matrix<Index,Dynamic,3> MatrixX3I;
   typedef Matrix<Index,Dynamic,2> MatrixX2I;
   typedef Matrix<Index,Dynamic,1> VectorXI;
@@ -117,7 +122,7 @@ IGL_INLINE void igl::mesh_boolean(
       F.block(0,0,FA.rows(),FA.cols()) = FA.rowwise().reverse();
       F.block(FA.rows(),0,FB.rows(),FB.cols()) = FB.array()+VA.rows();
       //F.block(0,0,FA.rows(),3) = FA;
-      //F.block(FA.rows(),0,FB.rows(),3) = 
+      //F.block(FA.rows(),0,FB.rows(),3) =
       //  FB.rowwise().reverse().array()+VA.rows();
       eff_type = MESH_BOOLEAN_TYPE_INTERSECT;
       break;
@@ -131,11 +136,11 @@ IGL_INLINE void igl::mesh_boolean(
   const auto & libigl_resolve = [](
     const MatrixX3S & V,
     const MatrixX3I & F,
-    MatrixX3S & CV,
+    MatrixX3ES & CV,
     MatrixX3I & CF,
     VectorXJ & J)
   {
-    MatrixX3S SV;
+    MatrixX3ES SV;
     MatrixX3I SF;
     MatrixX2I SIF;
     VectorXI SIM,UIM;
@@ -151,6 +156,7 @@ IGL_INLINE void igl::mesh_boolean(
   cout<<"resolve..."<<endl;
 #endif
   MatrixX3S CV;
+  MatrixX3ES EV;
   MatrixX3I CF;
   VectorXJ CJ;
   if(resolve_fun)
@@ -158,7 +164,12 @@ IGL_INLINE void igl::mesh_boolean(
     resolve_fun(V,F,CV,CF,CJ);
   }else
   {
-    libigl_resolve(V,F,CV,CF,CJ);
+    libigl_resolve(V,F,EV,CF,CJ);
+    CV.resize(EV.rows(), EV.cols());
+    std::transform(EV.data(), EV.data() + EV.rows()*EV.cols(),
+            CV.data(), [&](ExactScalar val) {
+            return CGAL::to_double(val);
+            });
   }
 
   if(type == MESH_BOOLEAN_TYPE_RESOLVE)
@@ -183,7 +194,7 @@ IGL_INLINE void igl::mesh_boolean(
   // peel layers keeping track of odd and even flips
   Matrix<bool,Dynamic,1> odd;
   Matrix<bool,Dynamic,1> flip;
-  peel_outer_hull_layers(CV,CF,CN,odd,flip);
+  peel_outer_hull_layers(EV,CF,CN,odd,flip);
 
 #ifdef IGL_MESH_BOOLEAN_DEBUG
   cout<<"categorize..."<<endl;
@@ -250,7 +261,7 @@ IGL_INLINE void igl::mesh_boolean(
       assert(ug < uG2G.size());
       uG2G[ug].push_back(g);
       // is uG(g,:) just a rotated version of G(g,:) ?
-      const bool consistent = 
+      const bool consistent =
         (G(g,0) == uG(ug,0) && G(g,1) == uG(ug,1) && G(g,2) == uG(ug,2)) ||
         (G(g,0) == uG(ug,1) && G(g,1) == uG(ug,2) && G(g,2) == uG(ug,0)) ||
         (G(g,0) == uG(ug,2) && G(g,1) == uG(ug,0) && G(g,2) == uG(ug,1));
@@ -293,6 +304,18 @@ IGL_INLINE void igl::mesh_boolean(
 }
 
 #ifdef IGL_STATIC_LIBRARY
+
+// This is a hack to discuss
+#include <igl/remove_unreferenced.cpp>
+
+template void igl::remove_unreferenced<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+
+#include <igl/cgal/peel_outer_hull_layers.cpp>
+template unsigned long igl::peel_outer_hull_layers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+
+#include <igl/cgal/outer_hull.cpp>
+template void igl::outer_hull<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+
 // Explicit template specialization
 template void igl::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/boundary_conditions.cpp

@@ -26,7 +26,6 @@ IGL_INLINE bool igl::boundary_conditions(
   Eigen::MatrixXd &       bc )
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
 
   if(P.size()+BE.rows() == 0)

include/igl/boundary_facets.cpp

@@ -21,7 +21,6 @@ IGL_INLINE void igl::boundary_facets(
   std::vector<std::vector<IntegerF> > & F)
 {
   using namespace std;
-  using namespace igl;
 
   if(T.size() == 0)
   {
@@ -111,7 +110,6 @@ IGL_INLINE void igl::boundary_facets(
   assert(T.cols() == 0 || T.cols() == 4 || T.cols() == 3);
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   // Cop out: use vector of vectors version
   vector<vector<typename Eigen::PlainObjectBase<DerivedT>::Scalar> > vT;
   matrix_to_list(T,vT);
@@ -134,6 +132,8 @@ Ret igl::boundary_facets(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
+// generated by autoexplicit.sh
 template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::boundary_facets<int, int>(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);

include/igl/boundary_loop.cpp

@@ -19,7 +19,6 @@ IGL_INLINE void igl::boundary_loop(
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
 
   if(F.rows() == 0)
     return;

include/igl/bounding_box_diagonal.cpp

@@ -13,7 +13,6 @@
 IGL_INLINE double igl::bounding_box_diagonal(
   const Eigen::MatrixXd & V)
 {
-  using namespace igl;
   using namespace Eigen;
   VectorXd maxV,minV;
   VectorXi maxVI,minVI;

include/igl/cat.cpp

@@ -123,9 +123,8 @@ IGL_INLINE Mat igl::cat(const int dim, const Mat & A, const Mat & B)
 }
 
 template <class Mat>
-IGL_INLINE void cat(const std::vector<std::vector< Mat > > & A, Mat & C)
+IGL_INLINE void igl::cat(const std::vector<std::vector< Mat > > & A, Mat & C)
 {
-  using namespace igl;
   using namespace std;
   // Start with empty matrix
   C.resize(0,0);

include/igl/cgal/CGAL_includes.hpp

@@ -8,6 +8,11 @@
 #ifndef IGL_CGAL_INCLUDES_H
 #define IGL_CGAL_INCLUDES_H
 
+// This causes unknown bugs during intersection meshing:
+//// http://www.alecjacobson.com/weblog/?p=4291
+//#define CGAL_INTERSECTION_VERSION 1
+// Use this instead to mute errors resulting from bad CGAL assertions
+#define CGAL_KERNEL_NO_ASSERTIONS
 // Triangle triangle intersection
 #include <CGAL/intersections.h>
 // THIS CANNOT BE INCLUDED IN THE SAME FILE AS <CGAL/intersections.h>

include/igl/cgal/RemeshSelfIntersectionsParam.h

@@ -0,0 +1,23 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_REMESH_SELF_INTERSECTIONS_PARAM_H
+#define IGL_REMESH_SELF_INTERSECTIONS_PARAM_H
+
+namespace igl
+{
+  // Optional Parameters
+  //   DetectOnly  Only compute IF, leave VV and FF alone
+  struct RemeshSelfIntersectionsParam
+  {
+    bool detect_only;
+    bool first_only;
+    RemeshSelfIntersectionsParam():detect_only(false),first_only(false){};
+  };
+}
+
+#endif

include/igl/cgal/SelfIntersectMesh.h

@@ -9,13 +9,14 @@
 #define IGL_SELFINTERSECTMESH_H
 
 #include "CGAL_includes.hpp"
-#include "remesh_self_intersections.h"
+#include "RemeshSelfIntersectionsParam.h"
 
 #include <Eigen/Dense>
 #include <list>
 #include <map>
 #include <vector>
 
+//#define IGL_SELFINTERSECTMESH_DEBUG
 #ifndef IGL_FIRST_HIT_EXCEPTION
 #define IGL_FIRST_HIT_EXCEPTION 10
 #endif
@@ -216,6 +217,11 @@ namespace igl
         SelfIntersectMesh * SIM, 
         const Box &a, 
         const Box &b);
+      // Annoying wrappers to conver from cgal to double or cgal
+      static inline void to_output_type(const typename Kernel::FT & cgal,double & d);
+      static inline void to_output_type(
+        const typename CGAL::Epeck::FT & cgal,
+        CGAL::Epeck::FT & d);
   };
 }
 
@@ -325,18 +331,22 @@ inline igl::SelfIntersectMesh<
   using namespace std;
   using namespace Eigen;
 
-  //const auto & tictoc = []()
-  //{
-  //  static double t_start = igl::get_seconds();
-  //  double diff = igl::get_seconds()-t_start;
-  //  t_start += diff;
-  //  return diff;
-  //};
-  //tictoc();
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  const auto & tictoc = []()
+  {
+    static double t_start = igl::get_seconds();
+    double diff = igl::get_seconds()-t_start;
+    t_start += diff;
+    return diff;
+  };
+  tictoc();
+#endif
 
   // Compute and process self intersections
   mesh_to_cgal_triangle_list(V,F,T);
-  //cout<<"mesh_to_cgal_triangle_list: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"mesh_to_cgal_triangle_list: "<<tictoc()<<endl;
+#endif
   // http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Box_intersection_d/Chapter_main.html#Section_63.5 
   // Create the corresponding vector of bounding boxes
   std::vector<Box> boxes;
@@ -355,7 +365,9 @@ inline igl::SelfIntersectMesh<
       // _1 etc. in global namespace)
       std::placeholders::_1,
       std::placeholders::_2);
-  //cout<<"boxes and bind: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"boxes and bind: "<<tictoc()<<endl;
+#endif
   // Run the self intersection algorithm with all defaults
   try{
     CGAL::box_self_intersection_d(boxes.begin(), boxes.end(),cb);
@@ -368,7 +380,9 @@ inline igl::SelfIntersectMesh<
     }
     // Otherwise just fall through
   }
-  //cout<<"box_self_intersection_d: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"box_self_intersection_d: "<<tictoc()<<endl;
+#endif
 
   // Convert lIF to Eigen matrix
   assert(lIF.size()%2 == 0);
@@ -387,7 +401,9 @@ inline igl::SelfIntersectMesh<
       i++;
     }
   }
-  //cout<<"IF: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"IF: "<<tictoc()<<endl;
+#endif
 
   if(params.detect_only)
   {
@@ -409,6 +425,9 @@ inline igl::SelfIntersectMesh<
   map<typename CDT_plus_2::Vertex_handle,Index> v2i;
   // Loop over offending triangles
   const size_t noff = offending.size();
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  double t_proj_del = 0;
+#endif
   // Unfortunately it looks like CGAL has trouble allocating memory when
   // multiple openmp threads are running. Crashes durring CDT...
 //# pragma omp parallel for if (noff>1000)
@@ -417,7 +436,13 @@ inline igl::SelfIntersectMesh<
     // index in F
     const Index f = offending[o];
     {
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+      const double t_before = get_seconds();
+#endif
       projected_delaunay(T[f],F_objects[f],cdt[o]);
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+      t_proj_del += (get_seconds()-t_before);
+#endif
     }
     // Q: Is this also delaunay in 3D?
     // A: No, because the projection is affine and delaunay is not affine
@@ -527,7 +552,10 @@ inline igl::SelfIntersectMesh<
       }
     }
   }
-  //cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
+#endif
+
   assert(NV_count == (Index)NV.size());
   // Build output
 #ifndef NDEBUG
@@ -565,7 +593,7 @@ inline igl::SelfIntersectMesh<
   }
   // Append vertices
   VV.resize(V.rows()+NV_count,3);
-  VV.block(0,0,V.rows(),3) = V;
+  VV.block(0,0,V.rows(),3) = V.template cast<typename DerivedVV::Scalar>();
   {
     Index i = 0;
     for(
@@ -576,13 +604,8 @@ inline igl::SelfIntersectMesh<
       for(Index d = 0;d<3;d++)
       {
         const Point_3 & p = *nvit;
-        VV(V.rows()+i,d) = CGAL::to_double(p[d]);
-        // This distinction does not seem necessary:
-//#ifdef INEXACT_CONSTRUCTION
-//        VV(V.rows()+i,d) = CGAL::to_double(p[d]);
-//#else
-//        VV(V.rows()+i,d) = CGAL::to_double(p[d].exact());
-//#endif
+        // Don't convert via double if output type is same as Kernel
+        to_output_type(p[d], VV(V.rows()+i,d));
       }
       i++;
     }
@@ -619,7 +642,9 @@ inline igl::SelfIntersectMesh<
       v++;
     }
   }
-  //cout<<"Output + dupes: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"Output + dupes: "<<tictoc()<<endl;
+#endif
 
   // Q: Does this give the same result as TETGEN?
   // A: For the cow and beast, yes.
@@ -1187,5 +1212,57 @@ inline void igl::SelfIntersectMesh<
   }
 }
 
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline 
+void 
+igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::to_output_type(
+    const typename Kernel::FT & cgal,
+    double & d)
+{
+  d = CGAL::to_double(cgal);
+}
+
+template <
+  typename Kernel,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedIF,
+  typename DerivedJ,
+  typename DerivedIM>
+inline 
+void 
+igl::SelfIntersectMesh<
+  Kernel,
+  DerivedV,
+  DerivedF,
+  DerivedVV,
+  DerivedFF,
+  DerivedIF,
+  DerivedJ,
+  DerivedIM>::to_output_type(
+    const typename CGAL::Epeck::FT & cgal,
+    CGAL::Epeck::FT & d)
+{
+  d = cgal;
+}
+
 #endif
 

include/igl/cgal/complex_to_mesh.cpp

@@ -11,7 +11,7 @@
 #include <igl/remove_unreferenced.h>
 
 #include <CGAL/Surface_mesh_default_triangulation_3.h>
-#include <CGAL/Triangulation_cell_base_with_circumcenter_3.h>
+#include <CGAL/Delaunay_triangulation_cell_base_with_circumcenter_3.h>
 #include <set>
 #include <stack>
 
@@ -147,6 +147,5 @@ IGL_INLINE bool igl::complex_to_mesh(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-template bool igl::complex_to_mesh<CGAL::Delaunay_triangulation_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_data_structure_3<CGAL::Surface_mesh_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Triangulation_cell_base_with_circumcenter_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Surface_mesh_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_cell_base_3<void> > > > >, CGAL::Default>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(CGAL::Complex_2_in_triangulation_3<CGAL::Delaunay_triangulation_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_data_structure_3<CGAL::Surface_mesh_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Triangulation_cell_base_with_circumcenter_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Surface_mesh_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_cell_base_3<void> > > > >, CGAL::Default>, void> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template bool igl::complex_to_mesh<CGAL::Delaunay_triangulation_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_data_structure_3<CGAL::Surface_mesh_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Delaunay_triangulation_cell_base_with_circumcenter_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Surface_mesh_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_cell_base_3<void> > > >, CGAL::Sequential_tag>, CGAL::Default, CGAL::Default>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(CGAL::Complex_2_in_triangulation_3<CGAL::Delaunay_triangulation_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_data_structure_3<CGAL::Surface_mesh_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Delaunay_triangulation_cell_base_with_circumcenter_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Surface_mesh_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_cell_base_3<void> > > >, CGAL::Sequential_tag>, CGAL::Default, CGAL::Default>, void> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/cgal/mesh_to_cgal_triangle_list.cpp

@@ -9,10 +9,13 @@
 
 #include <cassert>
 
-template <typename Kernel>
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Kernel>
 IGL_INLINE void igl::mesh_to_cgal_triangle_list(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
   std::vector<CGAL::Triangle_3<Kernel> > & T)
 {
   typedef CGAL::Point_3<Kernel>    Point_3;
@@ -35,7 +38,4 @@ IGL_INLINE void igl::mesh_to_cgal_triangle_list(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::mesh_to_cgal_triangle_list<CGAL::Epeck>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > >&);
-template void igl::mesh_to_cgal_triangle_list<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
-template void igl::mesh_to_cgal_triangle_list<CGAL::Simple_cartesian<double> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >&);
 #endif

include/igl/cgal/mesh_to_cgal_triangle_list.h

@@ -22,10 +22,13 @@ namespace igl
   //   F  #F by 3 list of triangle indices
   // Outputs:
   //   T  #F list of CGAL triangles
-  template <typename Kernel>
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename Kernel>
   IGL_INLINE void mesh_to_cgal_triangle_list(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
     std::vector<CGAL::Triangle_3<Kernel> > & T);
 }
 #ifndef IGL_STATIC_LIBRARY

include/igl/cgal/order_facets_around_edges.cpp

@@ -0,0 +1,247 @@
+#include "order_facets_around_edges.h"
+#include "../sort_angles.h"
+#include <type_traits>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedN,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename uE2oEType,
+    typename uE2CType >
+IGL_INLINE
+typename std::enable_if<!std::is_same<typename DerivedV::Scalar,
+typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
+igl::order_facets_around_edges(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedN>& N,
+        const Eigen::PlainObjectBase<DerivedE>& E,
+        const Eigen::PlainObjectBase<DeriveduE>& uE,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        std::vector<std::vector<uE2oEType> >& uE2oE,
+        std::vector<std::vector<uE2CType > >& uE2C ) {
+
+    typedef Eigen::Matrix<typename DerivedN::Scalar, 3, 1> Vector3F;
+    const typename DerivedV::Scalar EPS = 1e-12;
+    const size_t num_faces = F.rows();
+    const size_t num_undirected_edges = uE.rows();
+
+    auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
+    auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
+
+    uE2oE.resize(num_undirected_edges);
+    uE2C.resize(num_undirected_edges);
+
+    for(size_t ui = 0;ui<num_undirected_edges;ui++)
+    {
+        const auto& adj_edges = uE2E[ui];
+        const size_t edge_valance = adj_edges.size();
+        assert(edge_valance > 0);
+
+        const auto ref_edge = adj_edges[0];
+        const auto ref_face = edge_index_to_face_index(ref_edge);
+        Vector3F ref_normal = N.row(ref_face);
+
+        const auto ref_corner_o = edge_index_to_corner_index(ref_edge);
+        const auto ref_corner_s = (ref_corner_o+1)%3;
+        const auto ref_corner_d = (ref_corner_o+2)%3;
+
+        const typename DerivedF::Scalar o = F(ref_face, ref_corner_o);
+        const typename DerivedF::Scalar s = F(ref_face, ref_corner_s);
+        const typename DerivedF::Scalar d = F(ref_face, ref_corner_d);
+
+        Vector3F edge = V.row(d) - V.row(s);
+        auto edge_len = edge.norm();
+        bool degenerated = edge_len < EPS;
+        if (degenerated) {
+            if (edge_valance <= 2) {
+                // There is only one way to order 2 or less faces.
+                edge.setZero();
+            } else {
+                edge.setZero();
+                Eigen::Matrix<typename DerivedN::Scalar, Eigen::Dynamic, 3>
+                    normals(edge_valance, 3);
+                for (size_t fei=0; fei<edge_valance; fei++) {
+                    const auto fe = adj_edges[fei];
+                    const auto f = edge_index_to_face_index(fe);
+                    normals.row(fei) = N.row(f);
+                }
+                for (size_t i=0; i<edge_valance; i++) {
+                    size_t j = (i+1) % edge_valance;
+                    Vector3F ni = normals.row(i);
+                    Vector3F nj = normals.row(j);
+                    edge = ni.cross(nj);
+                    edge_len = edge.norm();
+                    if (edge_len >= EPS) {
+                        edge.normalize();
+                        break;
+                    }
+                }
+
+                // Ensure edge direction are consistent with reference face.
+                Vector3F in_face_vec = V.row(o) - V.row(s);
+                if (edge.cross(in_face_vec).dot(ref_normal) < 0) {
+                    edge *= -1;
+                }
+
+                if (edge.norm() < EPS) {
+                    std::cerr << "=====================================" << std::endl;
+                    std::cerr << "  ui: " << ui << std::endl;
+                    std::cerr << "edge: " << ref_edge << std::endl;
+                    std::cerr << "face: " << ref_face << std::endl;
+                    std::cerr << "  vs: " << V.row(s) << std::endl;
+                    std::cerr << "  vd: " << V.row(d) << std::endl;
+                    std::cerr << "adj face normals: " << std::endl;
+                    std::cerr << normals << std::endl;
+                    std::cerr << "Very degenerated case detected:" << std::endl;
+                    std::cerr << "Near zero edge surrounded by "
+                        << edge_valance << " neearly colinear faces" <<
+                        std::endl;
+                    std::cerr << "=====================================" << std::endl;
+                }
+            }
+        } else {
+            edge.normalize();
+        }
+
+        Eigen::MatrixXd angle_data(edge_valance, 3);
+        std::vector<bool> cons(edge_valance);
+
+        for (size_t fei=0; fei<edge_valance; fei++) {
+            const auto fe = adj_edges[fei];
+            const auto f = edge_index_to_face_index(fe);
+            const auto c = edge_index_to_corner_index(fe);
+            cons[fei] = (d == F(f, (c+1)%3));
+            assert( cons[fei] ||  (d == F(f,(c+2)%3)));
+            assert(!cons[fei] || (s == F(f,(c+2)%3)));
+            assert(!cons[fei] || (d == F(f,(c+1)%3)));
+            Vector3F n = N.row(f);
+            angle_data(fei, 0) = ref_normal.cross(n).dot(edge);
+            angle_data(fei, 1) = ref_normal.dot(n);
+            if (cons[fei]) {
+                angle_data(fei, 0) *= -1;
+                angle_data(fei, 1) *= -1;
+            }
+            angle_data(fei, 0) *= -1; // Sort clockwise.
+            angle_data(fei, 2) = (cons[fei]?1.:-1.)*(f+1);
+        }
+
+        Eigen::VectorXi order;
+        igl::sort_angles(angle_data, order);
+
+        auto& ordered_edges = uE2oE[ui];
+        auto& consistency = uE2C[ui];
+
+        ordered_edges.resize(edge_valance);
+        consistency.resize(edge_valance);
+        for (size_t fei=0; fei<edge_valance; fei++) {
+            ordered_edges[fei] = adj_edges[order[fei]];
+            consistency[fei] = cons[order[fei]];
+        }
+    }
+}
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedN,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename uE2oEType,
+    typename uE2CType >
+IGL_INLINE 
+typename std::enable_if<std::is_same<typename DerivedV::Scalar,
+typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
+igl::order_facets_around_edges(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedN>& N,
+        const Eigen::PlainObjectBase<DerivedE>& E,
+        const Eigen::PlainObjectBase<DeriveduE>& uE,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        std::vector<std::vector<uE2oEType> >& uE2oE,
+        std::vector<std::vector<uE2CType > >& uE2C ) {
+
+    typedef Eigen::Matrix<typename DerivedN::Scalar, 3, 1> Vector3F;
+    typedef Eigen::Matrix<typename DerivedV::Scalar, 3, 1> Vector3E;
+    const typename DerivedV::Scalar EPS = 1e-12;
+    const size_t num_faces = F.rows();
+    const size_t num_undirected_edges = uE.rows();
+
+    auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
+    auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
+
+    uE2oE.resize(num_undirected_edges);
+    uE2C.resize(num_undirected_edges);
+
+    for(size_t ui = 0;ui<num_undirected_edges;ui++)
+    {
+        const auto& adj_edges = uE2E[ui];
+        const size_t edge_valance = adj_edges.size();
+        assert(edge_valance > 0);
+
+        const auto ref_edge = adj_edges[0];
+        const auto ref_face = edge_index_to_face_index(ref_edge);
+        Vector3F ref_normal = N.row(ref_face);
+
+        const auto ref_corner_o = edge_index_to_corner_index(ref_edge);
+        const auto ref_corner_s = (ref_corner_o+1)%3;
+        const auto ref_corner_d = (ref_corner_o+2)%3;
+
+        const typename DerivedF::Scalar o = F(ref_face, ref_corner_o);
+        const typename DerivedF::Scalar s = F(ref_face, ref_corner_s);
+        const typename DerivedF::Scalar d = F(ref_face, ref_corner_d);
+
+        Vector3E exact_edge = V.row(d) - V.row(s);
+        exact_edge.array() /= exact_edge.squaredNorm();
+        Vector3F edge(
+                CGAL::to_double(exact_edge[0]),
+                CGAL::to_double(exact_edge[1]),
+                CGAL::to_double(exact_edge[2]));
+        edge.normalize();
+
+        Eigen::MatrixXd angle_data(edge_valance, 3);
+        std::vector<bool> cons(edge_valance);
+
+        for (size_t fei=0; fei<edge_valance; fei++) {
+            const auto fe = adj_edges[fei];
+            const auto f = edge_index_to_face_index(fe);
+            const auto c = edge_index_to_corner_index(fe);
+            cons[fei] = (d == F(f, (c+1)%3));
+            assert( cons[fei] ||  (d == F(f,(c+2)%3)));
+            assert(!cons[fei] || (s == F(f,(c+2)%3)));
+            assert(!cons[fei] || (d == F(f,(c+1)%3)));
+            Vector3F n = N.row(f);
+            angle_data(fei, 0) = ref_normal.cross(n).dot(edge);
+            angle_data(fei, 1) = ref_normal.dot(n);
+            if (cons[fei]) {
+                angle_data(fei, 0) *= -1;
+                angle_data(fei, 1) *= -1;
+            }
+            angle_data(fei, 0) *= -1; // Sort clockwise.
+            angle_data(fei, 2) = (cons[fei]?1.:-1.)*(f+1);
+        }
+
+        Eigen::VectorXi order;
+        igl::sort_angles(angle_data, order);
+
+        auto& ordered_edges = uE2oE[ui];
+        auto& consistency = uE2C[ui];
+
+        ordered_edges.resize(edge_valance);
+        consistency.resize(edge_valance);
+        for (size_t fei=0; fei<edge_valance; fei++) {
+            ordered_edges[fei] = adj_edges[order[fei]];
+            consistency[fei] = cons[order[fei]];
+        }
+    }
+}

include/igl/cgal/order_facets_around_edges.h

@@ -0,0 +1,79 @@
+#ifndef ORDER_FACETS_AROUND_EDGES
+#define ORDER_FACETS_AROUND_EDGES
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+namespace igl {
+    // For each undirected edge, sort its adjacent faces.
+    //
+    // Inputs:
+    //   V    #V by 3 list of vertices.
+    //   F    #F by 3 list of faces
+    //   N    #F by 3 list of face normals.
+    //   E    #F*3 by 2 list vertex indices, represents directed edges.
+    //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
+    //  EMAP  #F*3 list of indices that maps E to uE. (a many-to-one map)
+    //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
+    //
+    // Outputs:
+    //   uE2oE #uE list of lists that maps uE to an ordered list of E. (a
+    //         one-to-many map)
+    //   uE2C  #uE list of lists of bools indicates whether each face in
+    //         uE2oE[i] is consistently orientated as the ordering.
+    //
+    template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedN,
+        typename DerivedE,
+        typename DeriveduE,
+        typename DerivedEMAP,
+        typename uE2EType,
+        typename uE2oEType,
+        typename uE2CType >
+    IGL_INLINE
+    typename std::enable_if<!std::is_same<typename DerivedV::Scalar,
+    typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
+    order_facets_around_edges(
+            const Eigen::PlainObjectBase<DerivedV>& V,
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const Eigen::PlainObjectBase<DerivedN>& N,
+            const Eigen::PlainObjectBase<DerivedE>& E,
+            const Eigen::PlainObjectBase<DeriveduE>& uE,
+            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+            const std::vector<std::vector<uE2EType> >& uE2E,
+            std::vector<std::vector<uE2oEType> >& uE2oE,
+            std::vector<std::vector<uE2CType > >& uE2C );
+
+    template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedN,
+        typename DerivedE,
+        typename DeriveduE,
+        typename DerivedEMAP,
+        typename uE2EType,
+        typename uE2oEType,
+        typename uE2CType >
+    IGL_INLINE 
+    typename std::enable_if<std::is_same<typename DerivedV::Scalar,
+    typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
+    order_facets_around_edges(
+            const Eigen::PlainObjectBase<DerivedV>& V,
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const Eigen::PlainObjectBase<DerivedN>& N,
+            const Eigen::PlainObjectBase<DerivedE>& E,
+            const Eigen::PlainObjectBase<DeriveduE>& uE,
+            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+            const std::vector<std::vector<uE2EType> >& uE2E,
+            std::vector<std::vector<uE2oEType> >& uE2oE,
+            std::vector<std::vector<uE2CType > >& uE2C );
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "order_facets_around_edges.cpp"
+#endif
+
+#endif

include/igl/outer_hull.cpp → include/igl/cgal/outer_hull.cpp

@@ -1,18 +1,22 @@
 #include "outer_hull.h"
-#include "outer_facet.h"
-#include "sortrows.h"
-#include "facet_components.h"
-#include "winding_number.h"
-#include "triangle_triangle_adjacency.h"
-#include "unique_edge_map.h"
-#include "barycenter.h"
-#include "per_face_normals.h"
+#include "order_facets_around_edges.h"
+#include "../outer_facet.h"
+#include "../sortrows.h"
+#include "../facet_components.h"
+#include "../winding_number.h"
+#include "../triangle_triangle_adjacency.h"
+#include "../unique_edge_map.h"
+#include "../barycenter.h"
+#include "../per_face_normals.h"
+#include "../writePLY.h"
+#include "../sort_angles.h"
 
 #include <Eigen/Geometry>
 #include <vector>
 #include <map>
 #include <queue>
 #include <iostream>
+#include <CGAL/number_utils.h>
 //#define IGL_OUTER_HULL_DEBUG
 
 template <
@@ -30,9 +34,11 @@ IGL_INLINE void igl::outer_hull(
   Eigen::PlainObjectBase<DerivedJ> & J,
   Eigen::PlainObjectBase<Derivedflip> & flip)
 {
+#ifdef IGL_OUTER_HULL_DEBUG
+  std::cerr << "Extracting outer hull" << std::endl;
+#endif
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   typedef typename DerivedF::Index Index;
   Matrix<Index,DerivedF::RowsAtCompileTime,1> C;
   typedef Matrix<typename DerivedV::Scalar,Dynamic,DerivedV::ColsAtCompileTime> MatrixXV;
@@ -42,16 +48,17 @@ IGL_INLINE void igl::outer_hull(
   typedef Matrix<typename DerivedN::Scalar,1,3> RowVector3N;
   const Index m = F.rows();
 
-  const auto & duplicate_simplex = [&F](const int f, const int g)->bool
-  {
-    return 
-      (F(f,0) == F(g,0) && F(f,1) == F(g,1) && F(f,2) == F(g,2)) ||
-      (F(f,1) == F(g,0) && F(f,2) == F(g,1) && F(f,0) == F(g,2)) ||
-      (F(f,2) == F(g,0) && F(f,0) == F(g,1) && F(f,1) == F(g,2)) ||
-      (F(f,0) == F(g,2) && F(f,1) == F(g,1) && F(f,2) == F(g,0)) ||
-      (F(f,1) == F(g,2) && F(f,2) == F(g,1) && F(f,0) == F(g,0)) ||
-      (F(f,2) == F(g,2) && F(f,0) == F(g,1) && F(f,1) == F(g,0));
-  };
+  // UNUSED:
+  //const auto & duplicate_simplex = [&F](const int f, const int g)->bool
+  //{
+  //  return
+  //    (F(f,0) == F(g,0) && F(f,1) == F(g,1) && F(f,2) == F(g,2)) ||
+  //    (F(f,1) == F(g,0) && F(f,2) == F(g,1) && F(f,0) == F(g,2)) ||
+  //    (F(f,2) == F(g,0) && F(f,0) == F(g,1) && F(f,1) == F(g,2)) ||
+  //    (F(f,0) == F(g,2) && F(f,1) == F(g,1) && F(f,2) == F(g,0)) ||
+  //    (F(f,1) == F(g,2) && F(f,2) == F(g,1) && F(f,0) == F(g,0)) ||
+  //    (F(f,2) == F(g,2) && F(f,0) == F(g,1) && F(f,1) == F(g,0));
+  //};
 
 #ifdef IGL_OUTER_HULL_DEBUG
   cout<<"outer hull..."<<endl;
@@ -62,106 +69,31 @@ IGL_INLINE void igl::outer_hull(
 #endif
   typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixX2I;
   typedef Matrix<typename DerivedF::Index,Dynamic,1> VectorXI;
+  typedef Matrix<typename DerivedV::Scalar, 3, 1> Vector3F;
   MatrixX2I E,uE;
   VectorXI EMAP;
   vector<vector<typename DerivedF::Index> > uE2E;
   unique_edge_map(F,E,uE,EMAP,uE2E);
-
-  // TODO:
-  // uE --> face-edge index, sorted CCW around edge according to normal
-  // uE --> sorted order index 
-  // uE --> bool, whether needed to flip face to make "consistent" with unique
-  //   edge
-  // Place order of each half-edge in its corresponding sorted list around edge
-  VectorXI diIM(3*m);
-  // Whether face's edge used for sorting is consistent with unique edge
-  VectorXI dicons(3*m);
-  // dihedral angles of faces around edge with face of edge in dicons
-  vector<vector<typename DerivedV::Scalar> > di(uE2E.size());
-  // For each list of face-edges incide on a unique edge
-  for(size_t ui = 0;ui<(size_t)uE.rows();ui++)
-  {
-    // Base normal vector to orient against
-    const auto fe0 = uE2E[ui][0];
-    const RowVector3N & eVp = N.row(fe0%m);
-    MatrixXd di_I(uE2E[ui].size(),2);
-
-    const typename DerivedF::Scalar d = F(fe0%m,((fe0/m)+2)%3);
-    const typename DerivedF::Scalar s = F(fe0%m,((fe0/m)+1)%3);
-    // Edge vector
-    const auto & eV = (V.row(d)-V.row(s)).normalized();
-
-    vector<bool> cons(uE2E[ui].size());
-    // Loop over incident face edges
-    for(size_t fei = 0;fei<uE2E[ui].size();fei++)
-    {
-      const auto & fe = uE2E[ui][fei];
-      const auto f = fe % m;
-      const auto c = fe / m;
-      // source should match destination to be consistent
-      cons[fei] = (d == F(f,(c+1)%3));
-      assert( cons[fei] ||  (d == F(f,(c+2)%3)));
-      assert(!cons[fei] || (s == F(f,(c+2)%3)));
-      assert(!cons[fei] || (d == F(f,(c+1)%3)));
-      // Angle between n and f
-      const RowVector3N & n = N.row(f);
-      di_I(fei,0) = M_PI - atan2( eVp.cross(n).dot(eV), eVp.dot(n));
-      if(!cons[fei])
-      {
-        di_I(fei,0) = di_I(fei,0) + M_PI;
-        if(di_I(fei,0)>=2.*M_PI)
-        {
-          di_I(fei,0) = di_I(fei,0) - 2.*M_PI;
-        }
-      }
-      // This signing is very important to make sure different edges sort
-      // duplicate faces the same way, regardless of their orientations
-      di_I(fei,1) = (cons[fei]?1.:-1.)*f;
-    }
-
-    // Despite the effort to get stable normals the atan2 up doesn't
-    // compute (exactly) -θ for -n if it computes θ for n. So just
-    // explicitly check if there's a duplicate face
-    // Shitty O(val^2) implementation
-    for(size_t fei = 0;fei<uE2E[ui].size();fei++)
-    {
-      const auto & fe = uE2E[ui][fei];
-      const auto f = fe % m;
-      for(size_t gei = fei+1;gei<uE2E[ui].size();gei++)
-      {
-        const auto & ge = uE2E[ui][gei];
-        const auto g = ge % m;
-        if(duplicate_simplex(f,g))
-        {
 #ifdef IGL_OUTER_HULL_DEBUG
-          cout<<"Forcing duplicate: "<<(f+1)<<","<<(g+1)<<endl;
-#endif
-          di_I(gei,0) = di_I(fei,0);
-        }
+  for (size_t ui=0; ui<uE.rows(); ui++) {
+      std::cout << ui << ": " << uE2E[ui].size() << " -- (";
+      for (size_t i=0; i<uE2E[ui].size(); i++) {
+          std::cout << uE2E[ui][i] << ", ";
       }
-    }
-    VectorXi IM;
-    //igl::sort(di[ui],true,di[ui],IM);
-    // Sort, but break ties using "signed index" to ensure that duplicates
-    // always show up in same order.
-    MatrixXd s_di_I;
-    igl::sortrows(di_I,true,s_di_I,IM);
-    di[ui].resize(uE2E[ui].size());
-    for(size_t i = 0;i<di[ui].size();i++)
-    {
-      di[ui][i] = s_di_I(i,0);
-    }
-
-    // copy old list
-    vector<typename DerivedF::Index> temp = uE2E[ui];
-    for(size_t fei = 0;fei<uE2E[ui].size();fei++)
-    {
-      uE2E[ui][fei] = temp[IM(fei)];
-      const auto & fe = uE2E[ui][fei];
-      diIM(fe) = fei;
-      dicons(fe) = cons[IM(fei)];
-    }
+      std::cout << ")" << std::endl;
+  }
+#endif
 
+  std::vector<std::vector<typename DerivedF::Index> > uE2oE;
+  std::vector<std::vector<bool> > uE2C;
+  order_facets_around_edges(V, F, N, E, uE, EMAP, uE2E, uE2oE, uE2C);
+  uE2E = uE2oE;
+  VectorXI diIM(3*m);
+  for (auto ue : uE2E) {
+      for (size_t i=0; i<ue.size(); i++) {
+          auto fe = ue[i];
+          diIM[fe] = i;
+      }
   }
 
   vector<vector<vector<Index > > > TT,_1;
@@ -186,6 +118,7 @@ IGL_INLINE void igl::outer_hull(
   vector<MatrixXG> vG(ncc);
   vector<MatrixXJ> vJ(ncc);
   vector<MatrixXJ> vIM(ncc);
+  //size_t face_count = 0;
   for(size_t id = 0;id<ncc;id++)
   {
     vIM[id].resize(counts[id],1);
@@ -222,6 +155,9 @@ IGL_INLINE void igl::outer_hull(
     outer_facet(V,F,N,IM,f,f_flip);
 #ifdef IGL_OUTER_HULL_DEBUG
   cout<<"outer facet: "<<f<<endl;
+  cout << V.row(F(f, 0)) << std::endl;
+  cout << V.row(F(f, 1)) << std::endl;
+  cout << V.row(F(f, 2)) << std::endl;
 #endif
     int FHcount = 1;
     FH[f] = true;
@@ -231,6 +167,8 @@ IGL_INLINE void igl::outer_hull(
     Q.push(f+1*m);
     Q.push(f+2*m);
     flip(f) = f_flip;
+    //std::cout << "face " << face_count++ << ": " << f << std::endl;
+    //std::cout << "f " << F.row(f).array()+1 << std::endl;
     //cout<<"flip("<<f<<") = "<<(flip(f)?"true":"false")<<endl;
 #ifdef IGL_OUTER_HULL_DEBUG
   cout<<"BFS..."<<endl;
@@ -244,6 +182,12 @@ IGL_INLINE void igl::outer_hull(
       const int f = e%m;
       // corner
       const int c = e/m;
+#ifdef IGL_OUTER_HULL_DEBUG
+      std::cout << "edge: " << e << ", ue: " << EMAP(e) << std::endl;
+      std::cout << "face: " << f << std::endl;
+      std::cout << "corner: " << c << std::endl;
+      std::cout << "consistent: " << uE2C[EMAP(e)][diIM[e]] << std::endl;
+#endif
       // Should never see edge again...
       if(EH[e] == true)
       {
@@ -256,9 +200,24 @@ IGL_INLINE void igl::outer_hull(
       const int fd = flip(f)?F(f,(c+1)%3):F(f,(c+2)%3);
       // edge valence
       const size_t val = uE2E[EMAP(e)].size();
+#ifdef IGL_OUTER_HULL_DEBUG
+      std::cout << "vd: " << V.row(fd) << std::endl;
+      std::cout << "vs: " << V.row(fs) << std::endl;
+      std::cout << "edge: " << V.row(fd) - V.row(fs) << std::endl;
+      for (size_t i=0; i<val; i++) {
+          if (i == diIM(e)) {
+              std::cout << "* ";
+          } else {
+              std::cout << "  ";
+          }
+          std::cout << i << ": "
+              << " (e: " << uE2E[EMAP(e)][i] << ", f: "
+              << uE2E[EMAP(e)][i] % m * (uE2C[EMAP(e)][i] ? 1:-1) << ")" << std::endl;
+      }
+#endif
       //// find overlapping face-edges
       //const auto & neighbors = uE2E[EMAP(e)];
-      //// normal after possible flipping 
+      //// normal after possible flipping
       //const auto & fN = (flip(f)?-1.:1.)*N.row(f);
       //// Edge vector according to f's (flipped) orientation.
       ////const auto & eV = (V.row(fd)-V.row(fs)).normalized();
@@ -271,7 +230,7 @@ IGL_INLINE void igl::outer_hull(
       //const auto es = F(fe0%m,((fe0/m)+1)%3);
 
       // is edge consistent with edge of face used for sorting
-      const int e_cons = (dicons(e) ? 1: -1);
+      const int e_cons = (uE2C[EMAP(e)][diIM(e)] ? 1: -1);
       int nfei = -1;
       // Loop once around trying to find suitable next face
       for(size_t step = 1; step<val+2;step++)
@@ -279,25 +238,34 @@ IGL_INLINE void igl::outer_hull(
         const int nfei_new = (diIM(e) + 2*val + e_cons*step*(flip(f)?-1:1))%val;
         const int nf = uE2E[EMAP(e)][nfei_new] % m;
         // Don't consider faces with identical dihedral angles
-        if((di[EMAP(e)][diIM(e)] != di[EMAP(e)][nfei_new]))
+        //if ((di[EMAP(e)][diIM(e)].array() != di[EMAP(e)][nfei_new].array()).any())
+        //if((di[EMAP(e)][diIM(e)] != di[EMAP(e)][nfei_new]))
 //#warning "THIS IS HACK, FIX ME"
-//        if( abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new]) < 1e-16 )
+//        if( abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new]) < 1e-15 )
         {
 #ifdef IGL_OUTER_HULL_DEBUG
-        cout<<"Next facet: "<<(f+1)<<" --> "<<(nf+1)<<", |"<<
-          di[EMAP(e)][diIM(e)]<<" - "<<di[EMAP(e)][nfei_new]<<"| = "<<
-            abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new])
-            <<endl;
+        //cout<<"Next facet: "<<(f+1)<<" --> "<<(nf+1)<<", |"<<
+        //  di[EMAP(e)][diIM(e)]<<" - "<<di[EMAP(e)][nfei_new]<<"| = "<<
+        //    abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new])
+        //    <<endl;
 #endif
-        
+
 
 
           // Only use this face if not already seen
           if(!FH[nf])
           {
             nfei = nfei_new;
+          //} else {
+          //    std::cout << "skipping face " << nfei_new << " because it is seen before"
+          //        << std::endl;
           }
           break;
+        //} else {
+        //    std::cout << di[EMAP(e)][diIM(e)].transpose() << std::endl;
+        //    std::cout << di[EMAP(e)][diIM(nfei_new)].transpose() << std::endl;
+        //    std::cout << "skipping face " << nfei_new << " with identical dihedral angle"
+        //        << std::endl;
         }
 //#ifdef IGL_OUTER_HULL_DEBUG
 //        cout<<"Skipping co-planar facet: "<<(f+1)<<" --> "<<(nf+1)<<endl;
@@ -380,6 +348,8 @@ IGL_INLINE void igl::outer_hull(
         }
 #endif
         FH[nf] = true;
+        //std::cout << "face " << face_count++ << ": " << nf << std::endl;
+        //std::cout << "f " << F.row(nf).array()+1 << std::endl;
         FHcount++;
         // corner of neighbor
         const int nc = max_ne/m;
@@ -399,7 +369,7 @@ IGL_INLINE void igl::outer_hull(
         }
       }
     }
-    
+
     {
       vG[id].resize(FHcount,3);
       vJ[id].resize(FHcount,1);
@@ -427,18 +397,18 @@ IGL_INLINE void igl::outer_hull(
   // Is A inside B? Assuming A and B are consistently oriented but closed and
   // non-intersecting.
   const auto & is_component_inside_other = [](
-    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::MatrixXd & V,
     const MatrixXV & BC,
     const MatrixXG & A,
     const MatrixXJ & AJ,
     const MatrixXG & B)->bool
   {
     const auto & bounding_box = [](
-      const Eigen::PlainObjectBase<DerivedV> & V,
+      const Eigen::MatrixXd & V,
       const MatrixXG & F)->
-      MatrixXV
+      Eigen::MatrixXd
     {
-      MatrixXV BB(2,3);
+      Eigen::MatrixXd BB(2,3);
       BB<<
          1e26,1e26,1e26,
         -1e26,-1e26,-1e26;
@@ -456,8 +426,8 @@ IGL_INLINE void igl::outer_hull(
     };
     // A lot of the time we're dealing with unrelated, distant components: cull
     // them.
-    MatrixXV ABB = bounding_box(V,A);
-    MatrixXV BBB = bounding_box(V,B);
+    Eigen::MatrixXd ABB = bounding_box(V,A);
+    Eigen::MatrixXd BBB = bounding_box(V,B);
     if( (BBB.row(0)-ABB.row(1)).maxCoeff()>0  ||
         (ABB.row(0)-BBB.row(1)).maxCoeff()>0 )
     {
@@ -468,34 +438,31 @@ IGL_INLINE void igl::outer_hull(
     // POTENTIAL ROBUSTNESS WEAK AREA
     ////////////////////////////////////////////////////////////////////////
     //
-    // q could be so close (<~1e-16) to B that the winding number is not a robust way to
+
+    // winding_number_3 expects colmajor
+    // q could be so close (<~1e-15) to B that the winding number is not a robust way to
     // determine inside/outsideness. We could try to find a _better_ q which is
     // farther away, but couldn't they all be bad?
-    MatrixXV q = BC.row(AJ(0));
+    double q[3] = {
+        CGAL::to_double(BC(AJ(0), 0)),
+        CGAL::to_double(BC(AJ(0), 1)),
+        CGAL::to_double(BC(AJ(0), 2)) };
     // In a perfect world, it's enough to test a single point.
     double w;
-
-    // winding_number_3 expects colmajor
-    const typename DerivedV::Scalar * Vdata;
-    Vdata = V.data();
-    Matrix<
-      typename DerivedV::Scalar,
-      DerivedV::RowsAtCompileTime,
-      DerivedV::ColsAtCompileTime,
-      ColMajor> Vcol;
-    if(DerivedV::IsRowMajor)
-    {
-      // copy to convert to colmajor
-      Vcol = V;
-      Vdata = Vcol.data();
-    }
     winding_number_3(
-      Vdata,V.rows(),
+      V.data(),V.rows(),
       B.data(),B.rows(),
-      q.data(),1,&w);
-    return fabs(w)>0.5;
+      q,1,&w);
+    return w > 0.5 || w < -0.5;
   };
 
+  Eigen::MatrixXd Vcol(V.rows(), V.cols());
+  for (size_t i=0; i<(size_t)V.rows(); i++) {
+      for (size_t j=0; j<(size_t)V.cols(); j++) {
+          Vcol(i, j) = CGAL::to_double(V(i, j));
+      }
+  }
+
   // Reject components which are completely inside other components
   vector<bool> keep(ncc,true);
   size_t nG = 0;
@@ -508,7 +475,7 @@ IGL_INLINE void igl::outer_hull(
       {
         continue;
       }
-      const bool inside = is_component_inside_other(V,BC,vG[id],vJ[id],vG[oid]);
+      const bool inside = is_component_inside_other(Vcol,BC,vG[id],vJ[id],vG[oid]);
 #ifdef IGL_OUTER_HULL_DEBUG
       cout<<id<<" is inside "<<oid<<" ? "<<inside<<endl;
 #endif
@@ -556,7 +523,18 @@ IGL_INLINE void igl::outer_hull(
   return outer_hull(V,F,N,G,J,flip);
 }
 
+
 #ifdef IGL_STATIC_LIBRARY
+
+#include <igl/barycenter.cpp>
+template void igl::barycenter<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
+
+#include <igl/outer_facet.cpp>
+template void igl::outer_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
+
+#include <igl/cgal/order_facets_around_edges.cpp>
+template std::__1::enable_if<std::is_same<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>::Scalar, CGAL::Lazy_exact_nt<CGAL::Gmpq> >::value, void>::type igl::order_facets_around_edges<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > > const&, std::__1::vector<std::__1::vector<long, std::__1::allocator<long> >, std::__1::allocator<std::__1::vector<long, std::__1::allocator<long> > > >&, std::__1::vector<std::__1::vector<bool, std::__1::allocator<bool> >, std::__1::allocator<std::__1::vector<bool, std::__1::allocator<bool> > > >&);
+
 // Explicit template specialization
 template void igl::outer_hull<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
 template void igl::outer_hull<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/outer_hull.h → include/igl/cgal/outer_hull.h

@@ -1,6 +1,6 @@
 #ifndef IGL_OUTER_HULL_H
 #define IGL_OUTER_HULL_H
-#include "igl_inline.h"
+#include "../igl_inline.h"
 #include <Eigen/Core>
 namespace igl
 {

include/igl/peel_outer_hull_layers.cpp → include/igl/cgal/peel_outer_hull_layers.cpp

@@ -1,13 +1,13 @@
 #include "peel_outer_hull_layers.h"
-#include "per_face_normals.h"
+#include "../per_face_normals.h"
 #include "outer_hull.h"
 #include <vector>
 #include <iostream>
 //#define IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
 #ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
-#include "writePLY.h"
-#include "writeDMAT.h"
-#include "STR.h"
+#include "../writePLY.h"
+#include "../writeDMAT.h"
+#include "../STR.h"
 #endif
 
 using namespace std;

include/igl/peel_outer_hull_layers.h → include/igl/cgal/peel_outer_hull_layers.h

@@ -1,6 +1,6 @@
-#ifndef PEEL_OUTER_HULL_LAYERS_H
-#define PEEL_OUTER_HULL_LAYERS_H
-#include "igl_inline.h"
+#ifndef IGL_PEEL_OUTER_HULL_LAYERS_H
+#define IGL_PEEL_OUTER_HULL_LAYERS_H
+#include "../igl_inline.h"
 #include <Eigen/Core>
 namespace igl
 {

include/igl/cgal/point_mesh_squared_distance.cpp

@@ -9,7 +9,7 @@
 #include "mesh_to_cgal_triangle_list.h"
 
 template <typename Kernel>
-IGL_INLINE void igl::point_mesh_squared_distance(
+IGL_INLINE void igl::cgal::point_mesh_squared_distance(
   const Eigen::MatrixXd & P,
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
@@ -30,7 +30,7 @@ IGL_INLINE void igl::point_mesh_squared_distance(
 }
 
 template <typename Kernel>
-IGL_INLINE void igl::point_mesh_squared_distance_precompute(
+IGL_INLINE void igl::cgal::point_mesh_squared_distance_precompute(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   CGAL::AABB_tree<
@@ -78,7 +78,7 @@ IGL_INLINE void igl::point_mesh_squared_distance_precompute(
 }
 
 template <typename Kernel>
-IGL_INLINE void igl::point_mesh_squared_distance(
+IGL_INLINE void igl::cgal::point_mesh_squared_distance(
   const Eigen::MatrixXd & P,
   const CGAL::AABB_tree<
     CGAL::AABB_traits<Kernel, 
@@ -120,8 +120,8 @@ IGL_INLINE void igl::point_mesh_squared_distance(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::point_mesh_squared_distance_precompute<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
-template void igl::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
-template void igl::point_mesh_squared_distance<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-template void igl::point_mesh_squared_distance<CGAL::Simple_cartesian<double> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+template void igl::cgal::point_mesh_squared_distance_precompute<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
+template void igl::cgal::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
+template void igl::cgal::point_mesh_squared_distance<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+template void igl::cgal::point_mesh_squared_distance<CGAL::Simple_cartesian<double> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 #endif

include/igl/cgal/point_mesh_squared_distance.h

@@ -5,75 +5,77 @@
 // This Source Code Form is subject to the terms of the Mozilla Public License 
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef IGL_POINT_MESH_SQUARED_DISTANCE_H
-#define IGL_POINT_MESH_SQUARED_DISTANCE_H
+#ifndef IGL_CGAL_POINT_MESH_SQUARED_DISTANCE_H
+#define IGL_CGAL_POINT_MESH_SQUARED_DISTANCE_H
 #include <igl/igl_inline.h>
 #include <Eigen/Core>
 #include <vector>
 #include "CGAL_includes.hpp"
 namespace igl
 {
-  // Compute distances from a set of points P to a triangle mesh (V,F)
-  //
-  // Templates:
-  //   Kernal  CGAL computation and construction kernel (e.g.
-  //     CGAL::Simple_cartesian<double>)
-  // Inputs:
-  //   P  #P by 3 list of query point positions
-  //   V  #V by 3 list of vertex positions
-  //   F  #F by 3 list of triangle indices
-  // Outputs:
-  //   sqrD  #P list of smallest squared distances
-  //   I  #P list of facet indices corresponding to smallest distances
-  //   C  #P by 3 list of closest points
-  //
-  // Known bugs: This only computes distances to triangles. So unreferenced
-  // vertices are ignored.
-  template <typename Kernel>
-  IGL_INLINE void point_mesh_squared_distance(
-    const Eigen::MatrixXd & P,
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
-    Eigen::VectorXd & sqrD,
-    Eigen::VectorXi & I,
-    Eigen::MatrixXd & C);
-  // Probably can do this in a way that we don't pass around `tree` and `T`
-  //
-  // Outputs:
-  //   tree  CGAL's AABB tree
-  //   T  list of CGAL triangles in order of F (for determining which was found
-  //     in computation)
-  template <typename Kernel>
-  IGL_INLINE void point_mesh_squared_distance_precompute(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
-    CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+  namespace cgal
+  {
+    // Compute distances from a set of points P to a triangle mesh (V,F)
+    //
+    // Templates:
+    //   Kernal  CGAL computation and construction kernel (e.g.
+    //     CGAL::Simple_cartesian<double>)
+    // Inputs:
+    //   P  #P by 3 list of query point positions
+    //   V  #V by 3 list of vertex positions
+    //   F  #F by 3 list of triangle indices
+    // Outputs:
+    //   sqrD  #P list of smallest squared distances
+    //   I  #P list of facet indices corresponding to smallest distances
+    //   C  #P by 3 list of closest points
+    //
+    // Known bugs: This only computes distances to triangles. So unreferenced
+    // vertices and degenerate triangles (segments) are ignored.
+    template <typename Kernel>
+    IGL_INLINE void point_mesh_squared_distance(
+      const Eigen::MatrixXd & P,
+      const Eigen::MatrixXd & V,
+      const Eigen::MatrixXi & F,
+      Eigen::VectorXd & sqrD,
+      Eigen::VectorXi & I,
+      Eigen::MatrixXd & C);
+    // Probably can do this in a way that we don't pass around `tree` and `T`
+    //
+    // Outputs:
+    //   tree  CGAL's AABB tree
+    //   T  list of CGAL triangles in order of F (for determining which was found
+    //     in computation)
+    template <typename Kernel>
+    IGL_INLINE void point_mesh_squared_distance_precompute(
+      const Eigen::MatrixXd & V,
+      const Eigen::MatrixXi & F,
+      CGAL::AABB_tree<
+        CGAL::AABB_traits<Kernel, 
+          CGAL::AABB_triangle_primitive<Kernel, 
+            typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+          >
         >
-      >
-    > & tree,
-    std::vector<CGAL::Triangle_3<Kernel> > & T);
-  // Inputs:
-  //  see above
-  // Outputs:
-  //  see above
-  template <typename Kernel>
-  IGL_INLINE void point_mesh_squared_distance(
-    const Eigen::MatrixXd & P,
-    const CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      > & tree,
+      std::vector<CGAL::Triangle_3<Kernel> > & T);
+    // Inputs:
+    //  see above
+    // Outputs:
+    //  see above
+    template <typename Kernel>
+    IGL_INLINE void point_mesh_squared_distance(
+      const Eigen::MatrixXd & P,
+      const CGAL::AABB_tree<
+        CGAL::AABB_traits<Kernel, 
+          CGAL::AABB_triangle_primitive<Kernel, 
+            typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+          >
         >
-      >
-    > & tree,
-    const std::vector<CGAL::Triangle_3<Kernel> > & T,
-    Eigen::VectorXd & sqrD,
-    Eigen::VectorXi & I,
-    Eigen::MatrixXd & C);
-
+      > & tree,
+      const std::vector<CGAL::Triangle_3<Kernel> > & T,
+      Eigen::VectorXd & sqrD,
+      Eigen::VectorXi & I,
+      Eigen::MatrixXd & C);
+  }
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cgal/remesh_self_intersections.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "remesh_self_intersections.h"
 #include "SelfIntersectMesh.h"
@@ -45,7 +45,7 @@ IGL_INLINE void igl::remesh_self_intersections(
 //#endif
 
     typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-    typedef 
+    typedef
       SelfIntersectMesh<
         Kernel,
         DerivedV,
@@ -54,7 +54,7 @@ IGL_INLINE void igl::remesh_self_intersections(
         DerivedFF,
         DerivedIF,
         DerivedJ,
-        DerivedIM> 
+        DerivedIM>
       SelfIntersectMeshK;
     SelfIntersectMeshK SIM = SelfIntersectMeshK(V,F,params,VV,FF,IF,J,IM);
 
@@ -65,7 +65,7 @@ IGL_INLINE void igl::remesh_self_intersections(
   }else
   {
     typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
-    typedef 
+    typedef
       SelfIntersectMesh<
         Kernel,
         DerivedV,
@@ -74,7 +74,7 @@ IGL_INLINE void igl::remesh_self_intersections(
         DerivedFF,
         DerivedIF,
         DerivedJ,
-        DerivedIM> 
+        DerivedIM>
       SelfIntersectMeshK;
     SelfIntersectMeshK SIM = SelfIntersectMeshK(V,F,params,VV,FF,IF,J,IM);
   }

include/igl/cgal/remesh_self_intersections.h

@@ -8,6 +8,7 @@
 #ifndef IGL_REMESH_SELF_INTERSECTIONS_H
 #define IGL_REMESH_SELF_INTERSECTIONS_H
 #include <igl/igl_inline.h>
+#include "RemeshSelfIntersectionsParam.h"
 
 #include <Eigen/Dense>
 
@@ -17,18 +18,9 @@
 #  undef assert
 #  define assert( isOK ) ( (isOK) ? (void)0 : (void) mexErrMsgTxt(C_STR(__FILE__<<":"<<__LINE__<<": failed assertion `"<<#isOK<<"'"<<std::endl) ) )
 #endif
-
+  
 namespace igl
 {
-  // Optional Parameters
-  //   DetectOnly  Only compute IF, leave VV and FF alone
-  struct RemeshSelfIntersectionsParam
-  {
-    bool detect_only;
-    bool first_only;
-    RemeshSelfIntersectionsParam():detect_only(false),first_only(false){};
-  };
-  
   // Given a triangle mesh (V,F) compute a new mesh (VV,FF) which is the same
   // as (V,F) except that any self-intersecting triangles in (V,F) have been
   // subdivided (new vertices and face created) so that the self-intersection

include/igl/cgal/signed_distance.cpp

@@ -1,292 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-#include "signed_distance.h"
-#include "../per_vertex_normals.h"
-#include "../per_edge_normals.h"
-#include "../per_face_normals.h"
-#include "../get_seconds.h"
-#include "point_mesh_squared_distance.h"
-
-
-template <typename Kernel>
-IGL_INLINE void igl::signed_distance(
-  const Eigen::MatrixXd & P,
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  const SignedDistanceType sign_type,
-  Eigen::VectorXd & S,
-  Eigen::VectorXi & I,
-  Eigen::MatrixXd & C,
-  Eigen::MatrixXd & N)
-{
-  using namespace Eigen;
-  using namespace std;
-  assert(V.cols() == 3 && "V should have 3d positions");
-  assert(P.cols() == 3 && "P should have 3d positions");
-  assert(F.cols() == 3 && "F should have triangles");
-
-  typedef typename Kernel::FT FT;
-  typedef typename Kernel::Point_3 Point_3;
-  typedef typename CGAL::Triangle_3<Kernel> Triangle_3;
-  typedef typename std::vector<Triangle_3>::iterator Iterator;
-  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
-  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
-  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
-  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
-
-  // Prepare distance computation
-  Tree tree;
-  vector<Triangle_3 > T;
-  point_mesh_squared_distance_precompute(V,F,tree,T);
-
-  Eigen::MatrixXd FN,VN,EN;
-  Eigen::MatrixXi E;
-  Eigen::VectorXi EMAP;
-  WindingNumberAABB<Eigen::Vector3d> hier;
-  switch(sign_type)
-  {
-    default:
-      assert(false && "Unknown SignedDistanceType");
-    case SIGNED_DISTANCE_TYPE_DEFAULT:
-    case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
-      hier.set_mesh(V,F);
-      hier.grow();
-      break;
-    case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
-      // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
-      // [Bærentzen & Aanæs 2005]
-      per_face_normals(V,F,FN);
-      per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
-      per_edge_normals(
-        V,F,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);
-      N.resize(P.rows(),3);
-      break;
-  }
-
-  S.resize(P.rows(),1);
-  I.resize(P.rows(),1);
-  C.resize(P.rows(),3);
-  for(int p = 0;p<P.rows();p++)
-  {
-    const Point_3 q(P(p,0),P(p,1),P(p,2));
-    double s,sqrd;
-    Point_and_primitive_id pp;
-    switch(sign_type)
-    {
-      default:
-        assert(false && "Unknown SignedDistanceType");
-      case SIGNED_DISTANCE_TYPE_DEFAULT:
-      case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
-        signed_distance_winding_number(tree,hier,q,s,sqrd,pp);
-        break;
-      case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
-      {
-        Vector3d n(0,0,0);
-        signed_distance_pseudonormal(tree,T,F,FN,VN,EN,EMAP,q,s,sqrd,pp,n);
-        N.row(p) = n;
-        break;
-      }
-    }
-    I(p) = pp.second - T.begin();
-    S(p) = s*sqrt(sqrd);
-    C(p,0) = pp.first.x();
-    C(p,1) = pp.first.y();
-    C(p,2) = pp.first.z();
-  }
-}
-
-
-template <typename Kernel>
-IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
-  const CGAL::AABB_tree<
-    CGAL::AABB_traits<Kernel, 
-      CGAL::AABB_triangle_primitive<Kernel, 
-        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
-      >
-    >
-  > & tree,
-  const std::vector<CGAL::Triangle_3<Kernel> > & T,
-  const Eigen::MatrixXi & F,
-  const Eigen::MatrixXd & FN,
-  const Eigen::MatrixXd & VN,
-  const Eigen::MatrixXd & EN,
-  const Eigen::VectorXi & EMAP,
-  const typename Kernel::Point_3 & q)
-{
-  typename CGAL::AABB_tree<
-    CGAL::AABB_traits<Kernel, 
-      CGAL::AABB_triangle_primitive<Kernel, 
-        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
-        > 
-      >::Point_and_primitive_id pp;
-  typename Kernel::FT s,sqrd;
-  Eigen::Vector3d n(0,0,0);
-  signed_distance_pseudonormal<Kernel>(tree,T,F,FN,VN,EN,EMAP,q,s,sqrd,pp,n);
-  return s*sqrt(sqrd);
-}
-
-template <typename Kernel>
-IGL_INLINE void igl::signed_distance_pseudonormal(
-  const CGAL::AABB_tree<
-    CGAL::AABB_traits<Kernel, 
-      CGAL::AABB_triangle_primitive<Kernel, 
-        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
-      >
-    >
-  > & tree,
-  const std::vector<CGAL::Triangle_3<Kernel> > & T,
-  const Eigen::MatrixXi & F,
-  const Eigen::MatrixXd & FN,
-  const Eigen::MatrixXd & VN,
-  const Eigen::MatrixXd & EN,
-  const Eigen::VectorXi & EMAP,
-  const typename Kernel::Point_3 & q,
-  typename Kernel::FT & s,
-  typename Kernel::FT & sqrd,
-  typename CGAL::AABB_tree<
-    CGAL::AABB_traits<Kernel, 
-      CGAL::AABB_triangle_primitive<Kernel, 
-        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
-        > 
-      >::Point_and_primitive_id & pp,
-   Eigen::Vector3d & n)
-{
-  using namespace Eigen;
-  using namespace std;
-  typedef typename Kernel::FT FT;
-  typedef typename Kernel::Point_3 Point_3;
-  typedef typename CGAL::Triangle_3<Kernel> Triangle_3;
-  typedef typename std::vector<Triangle_3>::iterator Iterator;
-  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
-  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
-  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
-  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
-
-  pp = tree.closest_point_and_primitive(q);
-  Point_3 & p = pp.first;
-  const auto & qp = q-p;
-  sqrd = qp.squared_length();
-  Vector3d v(qp.x(),qp.y(),qp.z());
-  const int f = pp.second - T.begin();
-  const Triangle_3 & t = *pp.second;
-  // barycentric coordinates
-  const auto & area = [&p,&t](const int i, const int j)->FT
-  {
-    return sqrt(Triangle_3(p,t.vertex(i),t.vertex(j)).squared_area());
-  };
-  Vector3d b(area(1,2),area(2,0),area(0,1));
-  b /= b.sum();
-  // Determine which normal to use
-  const double epsilon = 1e-12;
-  const int type = (b.array()<=epsilon).cast<int>().sum();
-  switch(type)
-  {
-    case 2:
-      // Find vertex
-      for(int c = 0;c<3;c++)
-      {
-        if(b(c)>epsilon)
-        {
-          n = VN.row(F(f,c));
-          break;
-        }
-      }
-      break;
-    case 1:
-      // Find edge
-      for(int c = 0;c<3;c++)
-      {
-        if(b(c)<=epsilon)
-        {
-          n = EN.row(EMAP(F.rows()*c+f));
-          break;
-        }
-      }
-      break;
-    default:
-      assert(false && "all barycentric coords zero.");
-    case 0:
-      n = FN.row(f);
-      break;
-  }
-  s = (v.dot(n) >= 0 ? 1. : -1.);
-}
-
-template <typename Kernel>
-IGL_INLINE typename Kernel::FT igl::signed_distance_winding_number(
-  const CGAL::AABB_tree<
-    CGAL::AABB_traits<Kernel, 
-      CGAL::AABB_triangle_primitive<Kernel, 
-        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
-      >
-    >
-  > & tree,
-  const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
-  const typename Kernel::Point_3 & q)
-{
-  typename CGAL::AABB_tree<
-    CGAL::AABB_traits<Kernel, 
-      CGAL::AABB_triangle_primitive<Kernel, 
-        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
-        > 
-      >::Point_and_primitive_id pp;
-  typename Kernel::FT s,sqrd;
-  signed_distance_winding_number<Kernel>(tree,hier,q,s,sqrd,pp);
-  return s*sqrt(sqrd);
-}
-
-
-template <typename Kernel>
-IGL_INLINE void igl::signed_distance_winding_number(
-  const CGAL::AABB_tree<
-    CGAL::AABB_traits<Kernel, 
-      CGAL::AABB_triangle_primitive<Kernel, 
-        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
-      >
-    >
-  > & tree,
-  const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
-  const typename Kernel::Point_3 & q,
-  typename Kernel::FT & s,
-  typename Kernel::FT & sqrd,
-  typename CGAL::AABB_tree<
-    CGAL::AABB_traits<Kernel, 
-      CGAL::AABB_triangle_primitive<Kernel, 
-        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
-        > 
-      >::Point_and_primitive_id & pp)
-{
-  typedef typename Kernel::FT FT;
-  typedef typename Kernel::Point_3 Point_3;
-  typedef typename CGAL::Triangle_3<Kernel> Triangle_3;
-  typedef typename std::vector<Triangle_3>::iterator Iterator;
-  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
-  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
-  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
-  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
-  using namespace Eigen;
-  using namespace std;
-  using namespace igl;
-  pp = tree.closest_point_and_primitive(q);
-  Point_3 & p = pp.first;
-  const auto & qp = q-p;
-  const Vector3d eq(q.x(),q.y(),q.z()); 
-  sqrd = qp.squared_length();
-  const double w = hier.winding_number(eq);
-  s = 1.-2.*w;
-}
-
-#ifdef IGL_STATIC_LIBRARY
-// This template is necessary for the others to compile with clang
-// http://stackoverflow.com/questions/27748442/is-clangs-c11-support-reliable
-template void igl::signed_distance<CGAL::Epick>( const Eigen::MatrixXd & , const Eigen::MatrixXd & , const Eigen::MatrixXi & , const SignedDistanceType , Eigen::VectorXd & , Eigen::VectorXi &, Eigen::MatrixXd & , Eigen::MatrixXd & );
-template CGAL::Epick::FT igl::signed_distance_winding_number<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, CGAL::Epick::Point_3 const&);
-template CGAL::Epick::FT igl::signed_distance_pseudonormal<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, CGAL::Epick::Point_3 const&);
-template void igl::signed_distance_pseudonormal<CGAL::Simple_cartesian<double> >(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Simple_cartesian<double>, CGAL::AABB_triangle_primitive<CGAL::Simple_cartesian<double>, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >::iterator, CGAL::Boolean_tag<false> > > > const&, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, CGAL::Simple_cartesian<double>::Point_3 const&, CGAL::Simple_cartesian<double>::FT&, CGAL::Simple_cartesian<double>::FT&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Simple_cartesian<double>, CGAL::AABB_triangle_primitive<CGAL::Simple_cartesian<double>, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >::iterator, CGAL::Boolean_tag<false> > > >::Point_and_primitive_id&, Eigen::Matrix<double, 3, 1, 0, 3, 1>&);
-template void igl::signed_distance_winding_number<CGAL::Simple_cartesian<double> >(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Simple_cartesian<double>, CGAL::AABB_triangle_primitive<CGAL::Simple_cartesian<double>, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >::iterator, CGAL::Boolean_tag<false> > > > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, CGAL::Simple_cartesian<double>::Point_3 const&, CGAL::Simple_cartesian<double>::FT&, CGAL::Simple_cartesian<double>::FT&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Simple_cartesian<double>, CGAL::AABB_triangle_primitive<CGAL::Simple_cartesian<double>, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >::iterator, CGAL::Boolean_tag<false> > > >::Point_and_primitive_id&);
-#endif

include/igl/cgal/signed_distance.h

@@ -1,163 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef IGL_SIGNED_DISTANCE_H
-#define IGL_SIGNED_DISTANCE_H
-#include <igl/igl_inline.h>
-#include <igl/WindingNumberAABB.h>
-#include <Eigen/Core>
-#include <vector>
-#include "CGAL_includes.hpp"
-namespace igl
-{
-  enum SignedDistanceType
-  {
-    SIGNED_DISTANCE_TYPE_PSEUDONORMAL   = 0,
-    SIGNED_DISTANCE_TYPE_WINDING_NUMBER = 1,
-    SIGNED_DISTANCE_TYPE_DEFAULT        = 2,
-    SIGNED_DISTANCE_TYPE_IGNORE         = 3,
-    NUM_SIGNED_DISTANCE_TYPE            = 4
-  };
-  // Computes signed distance to a mesh
-  //
-  // Templates:
-  //   Kernal  CGAL computation and construction kernel (e.g.
-  //     CGAL::Simple_cartesian<double>)
-  // Inputs:
-  //   P  #P by 3 list of query point positions
-  //   V  #V by 3 list of vertex positions
-  //   F  #F by 3 list of triangle indices
-  //   sign_type  method for computing distance _sign_ (see signed_distance.h)
-  // Outputs:
-  //   S  #P list of smallest signed distances
-  //   I  #P list of facet indices corresponding to smallest distances
-  //   C  #P by 3 list of closest points
-  //   N  #P by 3 list of closest normals (only set if
-  //     sign_type=SIGNED_DISTANCE_TYPE_PSEUDONORMAL)
-  //
-  // Known bugs: This only computes distances to triangles. So unreferenced
-  // vertices are ignored.
-  template <typename Kernel>
-  IGL_INLINE void signed_distance(
-    const Eigen::MatrixXd & P,
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
-    const SignedDistanceType sign_type,
-    Eigen::VectorXd & S,
-    Eigen::VectorXi & I,
-    Eigen::MatrixXd & C,
-    Eigen::MatrixXd & N);
-  // Computes signed distance to mesh
-  //
-  // Templates:
-  //   Kernal  CGAL computation and construction kernel (e.g.
-  //     CGAL::Simple_cartesian<double>)
-  // Inputs:
-  //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
-  //   T  #F list of CGAL triangles (see point_mesh_squared_distance.h)
-  //   F  #F by 3 list of triangle indices
-  //   FN  #F by 3 list of triangle normals 
-  //   VN  #V by 3 list of vertex normals (ANGLE WEIGHTING)
-  //   EN  #E by 3 list of edge normals (UNIFORM WEIGHTING)
-  //   EMAP  #F*3 mapping edges in F to E
-  //   q  Query point
-  // Returns signed distance to mesh
-  //
-  template <typename Kernel>
-  IGL_INLINE typename Kernel::FT signed_distance_pseudonormal(
-    const CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
-        >
-      >
-    > & tree,
-    const std::vector<CGAL::Triangle_3<Kernel> > & T,
-    const Eigen::MatrixXi & F,
-    const Eigen::MatrixXd & FN,
-    const Eigen::MatrixXd & VN,
-    const Eigen::MatrixXd & EN,
-    const Eigen::VectorXi & EMAP,
-    const typename Kernel::Point_3 & q);
-  // Outputs:
-  //   s  sign
-  //   sqrd  squared distance
-  //   pp  closest point and primitve
-  //   n  normal
-  template <typename Kernel>
-  IGL_INLINE void signed_distance_pseudonormal(
-    const CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
-        >
-      >
-    > & tree,
-    const std::vector<CGAL::Triangle_3<Kernel> > & T,
-    const Eigen::MatrixXi & F,
-    const Eigen::MatrixXd & FN,
-    const Eigen::MatrixXd & VN,
-    const Eigen::MatrixXd & EN,
-    const Eigen::VectorXi & EMAP,
-    const typename Kernel::Point_3 & q,
-    typename Kernel::FT & s,
-    typename Kernel::FT & sqrd,
-    typename CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
-          > 
-        >::Point_and_primitive_id & pp,
-   Eigen::Vector3d & n);
-
-  // Inputs:
-  //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
-  //   hier  Winding number evaluation hierarchy
-  //   q  Query point
-  // Returns signed distance to mesh
-  template <typename Kernel>
-  IGL_INLINE typename Kernel::FT signed_distance_winding_number(
-    const CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
-        >
-      >
-    > & tree,
-    const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
-    const typename Kernel::Point_3 & q);
-  // Outputs:
-  //   s  sign
-  //   sqrd  squared distance
-  //   pp  closest point and primitve
-  template <typename Kernel>
-  IGL_INLINE void signed_distance_winding_number(
-    const CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
-        >
-      >
-    > & tree,
-    const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
-    const typename Kernel::Point_3 & q,
-    typename Kernel::FT & s,
-    typename Kernel::FT & sqrd,
-    typename CGAL::AABB_tree<
-      CGAL::AABB_traits<Kernel, 
-        CGAL::AABB_triangle_primitive<Kernel, 
-          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
-          > 
-        >::Point_and_primitive_id & pp);
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "signed_distance.cpp"
-#endif
-
-#endif
-

include/igl/cgal/signed_distance_isosurface.cpp

@@ -8,21 +8,20 @@
 #include "signed_distance_isosurface.h"
 #include "point_mesh_squared_distance.h"
 #include "complex_to_mesh.h"
-#include "signed_distance.h"
 
-#include <igl/per_face_normals.h>
-#include <igl/per_edge_normals.h>
-#include <igl/per_vertex_normals.h>
-#include <igl/centroid.h>
-#include <igl/WindingNumberAABB.h>
-#include <igl/matlab_format.h>
-#include <igl/remove_unreferenced.h>
+#include "../AABB.h"
+#include "../per_face_normals.h"
+#include "../per_edge_normals.h"
+#include "../per_vertex_normals.h"
+#include "../centroid.h"
+#include "../WindingNumberAABB.h"
+#include "../matlab_format.h"
+#include "../remove_unreferenced.h"
 
 #include <CGAL/Surface_mesh_default_triangulation_3.h>
 #include <CGAL/Complex_2_in_triangulation_3.h>
 #include <CGAL/make_surface_mesh.h>
 #include <CGAL/Implicit_surface_3.h>
-#include <CGAL/Polyhedron_3.h>
 #include <CGAL/IO/output_surface_facets_to_polyhedron.h>
 // Axis-aligned bounding box tree for tet tri intersection
 #include <CGAL/AABB_tree.h>
@@ -42,6 +41,7 @@ IGL_INLINE bool igl::signed_distance_isosurface(
   Eigen::MatrixXi & F)
 {
   using namespace std;
+  using namespace Eigen;
 
   // default triangulation for Surface_mesher
   typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
@@ -53,18 +53,9 @@ IGL_INLINE bool igl::signed_distance_isosurface(
   typedef GT::FT FT;
   typedef std::function<FT (Point_3)> Function;
   typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;
-  typedef CGAL::Polyhedron_3<GT> Polyhedron;
-  typedef GT::Kernel Kernel;
-  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
-  typedef typename std::vector<Triangle_3>::iterator Iterator;
-  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
-  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
-  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
-  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
 
-  Tree tree;
-  vector<Triangle_3 > T;
-  point_mesh_squared_distance_precompute(IV,IF,tree,T);
+  AABB<Eigen::MatrixXd,3> tree;
+  tree.init(IV,IF);
 
   Eigen::MatrixXd FN,VN,EN;
   Eigen::MatrixXi E;
@@ -74,6 +65,9 @@ IGL_INLINE bool igl::signed_distance_isosurface(
   {
     default:
       assert(false && "Unknown SignedDistanceType");
+    case SIGNED_DISTANCE_TYPE_UNSIGNED:
+      // do nothing
+      break;
     case SIGNED_DISTANCE_TYPE_DEFAULT:
     case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
       hier.set_mesh(IV,IF);
@@ -105,39 +99,36 @@ IGL_INLINE bool igl::signed_distance_isosurface(
   {
     default:
       assert(false && "Unknown SignedDistanceType");
+    case SIGNED_DISTANCE_TYPE_UNSIGNED:
+      fun = 
+        [&tree,&IV,&IF,&level](const Point_3 & q) -> FT
+        {
+          int i;
+          RowVector3d c;
+          const double sd = tree.squared_distance(
+            IV,IF,RowVector3d(q.x(),q.y(),q.z()),i,c);
+          return sd-level;
+        };
     case SIGNED_DISTANCE_TYPE_DEFAULT:
     case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
       fun = 
-        [&tree,&hier,&level](const Point_3 & q) -> FT
+        [&tree,&IV,&IF,&hier,&level](const Point_3 & q) -> FT
         {
-          return signed_distance_winding_number(tree,hier,q)-level;
+          const double sd = signed_distance_winding_number(
+            tree,IV,IF,hier,RowVector3d(q.x(),q.y(),q.z()));
+          return sd-level;
         };
       break;
     case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
-      fun = [&tree,&T,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) -> FT
+      fun = [&tree,&IV,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) -> FT
         {
-          return 
-            igl::signed_distance_pseudonormal(tree,T,IF,FN,VN,EN,EMAP,q) - 
-            level;
+          const double sd = 
+            igl::signed_distance_pseudonormal(
+              tree,IV,IF,FN,VN,EN,EMAP,RowVector3d(q.x(),q.y(),q.z()));
+          return sd- level;
         };
       break;
   }
-    //[&tree,&hier,&T,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) ->FT
-    //{
-    //  const FT p = signed_distance_pseudonormal(tree,T,IF,FN,VN,EN,EMAP,q);
-    //  const FT w = signed_distance_winding_number(tree,hier,q);
-    //  if(w*p < 0 && (fabs(w) > 0.1 || fabs(p) > 0.1))
-    //  {
-    //    cout<<"q=["<<q.x()<<","<<q.y()<<","<<q.z()<<"];"<<endl;
-    //    cout<<matlab_format(n.transpose().eval(),"n")<<endl;
-    //    cout<<matlab_format(b.transpose().eval(),"b")<<endl;
-    //    cout<<"Sig difference: "<<type<<endl;
-    //    cout<<"w: "<<w<<endl;
-    //    cout<<"p: "<<p<<endl;
-    //    exit(1);
-    //  }
-    //  return w;
-    //},
   Sphere_3 bounding_sphere(cmid, (r+level)*(r+level));
   Surface_3 surface(fun,bounding_sphere);
   CGAL::Surface_mesh_default_criteria_3<Tr> 

include/igl/cgal/signed_distance_isosurface.h

@@ -7,8 +7,8 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_SIGNED_DISTANCE_ISOSURFACE_H
 #define IGL_SIGNED_DISTANCE_ISOSURFACE_H
-#include <igl/igl_inline.h>
-#include "signed_distance.h"
+#include "../igl_inline.h"
+#include "../signed_distance.h"
 #include <Eigen/Core>
 namespace igl
 {
@@ -18,11 +18,12 @@ namespace igl
   // Inputs:
   //   IV  #IV by 3 list of input mesh vertex positions
   //   IF  #IF by 3 list of input triangle indices
-  //   level  iso-level to contour (e.g. 0)
+  //   level  iso-level to contour in world coords, negative is inside.
   //   angle_bound  lower bound on triangle angles (mesh quality) (e.g. 28)
   //   radius_bound  upper bound on triangle size (mesh density?) (e.g. 0.02)
   //   distance_bound  cgal mysterious parameter (mesh density?) (e.g. 0.01)
-  //   sign_type  method for computing distance _sign_ (see signed_distance.h)
+  //   sign_type  method for computing distance _sign_ (see
+  //     ../signed_distance.h)
   // Outputs:
   //   V  #V by 3 list of input mesh vertex positions
   //   F  #F by 3 list of input triangle indices

include/igl/collapse_edge.cpp

@@ -29,7 +29,6 @@ IGL_INLINE bool igl::collapse_edge(
   // never get collapsed to anything else since it is the smallest index)
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   const int eflip = E(e,0)>E(e,1);
   // source and destination
   const int s = eflip?E(e,1):E(e,0);

include/igl/comb_cross_field.cpp

@@ -144,4 +144,5 @@ IGL_INLINE void igl::comb_cross_field(const Eigen::PlainObjectBase<DerivedV> &V,
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::comb_cross_field<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::comb_cross_field<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/comb_frame_field.cpp

@@ -72,4 +72,5 @@ IGL_INLINE void igl::comb_frame_field(const Eigen::PlainObjectBase<DerivedV> &V,
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::comb_frame_field<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::comb_frame_field<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/comiso/miq.cpp.REMOVED.git-id

@@ -1 +1 @@
-7417ec66154d8aa8468ea9b36c103ed6e7899ec7
+53de49ce1ce4a0b841b852e3e973c71387a7bcae

include/igl/comiso/nrosy.cpp

@@ -6,6 +6,8 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 
+#include "nrosy.h"
+
 #include <igl/comiso/nrosy.h>
 #include <igl/triangle_triangle_adjacency.h>
 #include <igl/edge_topology.h>
@@ -22,6 +24,7 @@
 #include <CoMISo/Solver/ConstrainedSolver.hh>
 #include <CoMISo/Solver/MISolver.hh>
 #include <CoMISo/Solver/GMM_Tools.hh>
+
 namespace igl
 {
 class NRosyField

include/igl/compile_and_link_program.cpp

@@ -0,0 +1,33 @@
+#include "compile_and_link_program.h"
+#include "compile_shader.h"
+#include "report_gl_error.h"
+#include <iostream>
+#include <cassert>
+
+IGL_INLINE GLuint igl::compile_and_link_program(
+  const char * v_str, const char * f_str)
+{
+  GLuint vid = compile_shader(GL_VERTEX_SHADER,v_str);
+  GLuint fid = compile_shader(GL_FRAGMENT_SHADER,f_str);
+
+  GLuint prog_id = glCreateProgram();
+  assert(prog_id != 0 && "Failed to create shader.");
+  glAttachShader(prog_id,vid);
+  igl::report_gl_error("glAttachShader (vid): ");
+  glAttachShader(prog_id,fid);
+  igl::report_gl_error("glAttachShader (fid): ");
+
+  glLinkProgram(prog_id);
+  igl::report_gl_error("glLinkProgram: ");
+
+  GLint status;
+  glGetProgramiv(prog_id, GL_LINK_STATUS, &status);
+  if (status != GL_TRUE)
+  {
+    char buffer[512];
+    glGetProgramInfoLog(prog_id, 512, NULL, buffer);
+    std::cerr << "Linker error: " << std::endl << buffer << std::endl;
+    prog_id = 0;
+  }
+  return prog_id;
+}

include/igl/compile_and_link_program.h

@@ -0,0 +1,22 @@
+#ifndef IGL_COMPILE_AND_LINK_PROGRAM_H
+#define IGL_COMPILE_AND_LINK_PROGRAM_H
+#include "igl_inline.h"
+#include "OpenGL_convenience.h"
+namespace igl
+{
+  // Compile and link very simple vertex/fragment shaders
+  //
+  // Inputs:
+  //   v_str  string of vertex shader contents
+  //   f_str  string of fragment shader contents
+  // Returns id of program
+  //
+  // Known bugs: this seems to duplicate `create_shader_program` with less
+  // functionality.
+  IGL_INLINE GLuint compile_and_link_program(
+    const char * v_str, const char * f_str);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "compile_and_link_program.cpp"
+#endif
+#endif

include/igl/compile_shader.cpp

@@ -0,0 +1,28 @@
+#include "compile_shader.h"
+#include "report_gl_error.h"
+#include <iostream>
+
+IGL_INLINE GLuint igl::compile_shader(const GLint type, const char * str)
+{
+  GLuint id = glCreateShader(type);
+  igl::report_gl_error("glCreateShader: ");
+  glShaderSource(id,1,&str,NULL);
+  igl::report_gl_error("glShaderSource: ");
+  glCompileShader(id);
+  igl::report_gl_error("glCompileShader: ");
+
+  GLint status;
+  glGetShaderiv(id, GL_COMPILE_STATUS, &status);
+  if (status != GL_TRUE)
+  {
+    char buffer[512];
+    if (type == GL_VERTEX_SHADER)
+      std::cerr << "Vertex shader:" << std::endl;
+    else if (type == GL_FRAGMENT_SHADER)
+      std::cerr << "Fragment shader:" << std::endl;
+    std::cerr << str << std::endl << std::endl;
+    glGetShaderInfoLog(id, 512, NULL, buffer);
+    std::cerr << "Error: " << std::endl << buffer << std::endl;
+  }
+  return id;
+}

include/igl/compile_shader.h

@@ -0,0 +1,28 @@
+#ifndef IGL_COMPILE_SHADER_H
+#define IGL_COMPILE_SHADER_H
+#include "OpenGL_convenience.h"
+#include "igl_inline.h"
+namespace igl
+{
+  // Compile a shader given type and string of shader code
+  // 
+  // Inputs:
+  //   type  either GL_VERTEX_SHADER or GL_FRAGMENT_SHADER
+  //   str  contents of shader code
+  // Returns result of glCreateShader (id of shader)
+  //
+  // Example:
+  //     GLuint vid = compile_shader(GL_VERTEX_SHADER,vertex_shader.c_str());
+  //     GLuint fid = compile_shader(GL_FRAGMENT_SHADER,fragment_shader.c_str());
+  //     GLuint prog_id = glCreateProgram();
+  //     glAttachShader(prog_id,vid);
+  //     glAttachShader(prog_id,fid);
+  //     glLinkProgram(prog_id);
+  //
+  // Known bugs: seems to be duplicate of `load_shader`
+  IGL_INLINE GLuint compile_shader(const GLint type, const char * str);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "compile_shader.cpp"
+#endif
+#endif

include/igl/compute_frame_field_bisectors.cpp

@@ -75,4 +75,5 @@ IGL_INLINE void igl::compute_frame_field_bisectors(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::compute_frame_field_bisectors<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::compute_frame_field_bisectors<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/cotmatrix.cpp

@@ -22,7 +22,6 @@ IGL_INLINE void igl::cotmatrix(
   const Eigen::PlainObjectBase<DerivedF> & F, 
   Eigen::SparseMatrix<Scalar>& L)
 {
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
 

include/igl/cotmatrix_entries.cpp

@@ -21,7 +21,6 @@ IGL_INLINE void igl::cotmatrix_entries(
   const Eigen::PlainObjectBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedC>& C)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // simplex size (3: triangles, 4: tetrahedra)

include/igl/covariance_scatter_matrix.cpp

@@ -20,7 +20,6 @@ IGL_INLINE void igl::covariance_scatter_matrix(
   const ARAPEnergyType energy,
   Eigen::SparseMatrix<double>& CSM)
 {
-  using namespace igl;
   using namespace Eigen;
   // number of mesh vertices
   int n = V.rows();

include/igl/create_shader_program.cpp

@@ -10,20 +10,22 @@
 
 #include "load_shader.h"
 #include "print_program_info_log.h"
+#include <iostream>
 #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,
+  const std::string & geom_source,
+  const std::string & vert_source,
+  const std::string & frag_source,
+  const std::map<std::string,GLuint> & attrib,
   GLuint & id)
 {
+  using namespace std;
   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");
+    cerr<<
+      "create_shader_program() could not create shader program,"
+      " both .vert and .frag source given were empty"<<endl;
     return false;
   }
 
@@ -31,18 +33,29 @@ IGL_INLINE bool igl::create_shader_program(
   id = glCreateProgram();
   if(id == 0)
   {
-    fprintf(
-      stderr,
-      "Error: create_shader_program() could not create shader program.\n");
+    cerr<<"create_shader_program() could not create shader program."<<endl;
     return false;
   }
 
+  if(geom_source != "")
+  {
+    // load vertex shader
+    GLuint g = igl::load_shader(geom_source.c_str(),GL_GEOMETRY_SHADER_EXT);
+    if(g == 0)
+    {
+      cerr<<"geometry shader failed to compile."<<endl;
+      return false;
+    }
+    glAttachShader(id,g);
+  }
+
   if(vert_source != "")
   {
     // load vertex shader
     GLuint v = igl::load_shader(vert_source.c_str(),GL_VERTEX_SHADER);
     if(v == 0)
     {
+      cerr<<"vertex shader failed to compile."<<endl;
       return false;
     }
     glAttachShader(id,v);
@@ -54,6 +67,7 @@ IGL_INLINE bool igl::create_shader_program(
     GLuint f = igl::load_shader(frag_source.c_str(),GL_FRAGMENT_SHADER);
     if(f == 0)
     {
+      cerr<<"fragment shader failed to compile."<<endl;
       return false;
     }
     glAttachShader(id,f);
@@ -78,4 +92,36 @@ IGL_INLINE bool igl::create_shader_program(
 
   return true;
 }
+
+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 & prog_id)
+{
+  return create_shader_program("",vert_source,frag_source,attrib,prog_id);
+}
+
+
+IGL_INLINE GLuint igl::create_shader_program(
+  const std::string & geom_source,
+  const std::string & vert_source,
+  const std::string & frag_source,
+  const std::map<std::string,GLuint> & attrib)
+{
+  GLuint prog_id = 0;
+  create_shader_program(geom_source,vert_source,frag_source,attrib,prog_id);
+  return prog_id;
+}
+
+IGL_INLINE GLuint igl::create_shader_program(
+  const std::string & vert_source,
+  const std::string & frag_source,
+  const std::map<std::string,GLuint> & attrib)
+{
+  GLuint prog_id = 0;
+  create_shader_program(vert_source,frag_source,attrib,prog_id);
+  return prog_id;
+}
+
 #endif

include/igl/create_shader_program.h

@@ -21,6 +21,8 @@ namespace igl
   // source strings and vertex attributes assigned from a map before linking the
   // shaders to the program, making it ready to use with glUseProgram(id)
   // Inputs:
+  //   geom_source  string containing source code of geometry shader (can be
+  //     "" to mean use default pass-through)
   //   vert_source  string containing source code of vertex shader
   //   frag_source  string containing source code of fragment shader
   //   attrib  map containing table of vertex attribute strings add their
@@ -34,10 +36,25 @@ namespace igl
   //
   // See also: destroy_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,
+    const std::string &geom_source,
+    const std::string &vert_source,
+    const std::string &frag_source,
+    const std::map<std::string,GLuint> &attrib,
     GLuint & id);
+  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);
+  IGL_INLINE GLuint create_shader_program(
+    const std::string & geom_source,
+    const std::string & vert_source,
+    const std::string & frag_source,
+    const std::map<std::string,GLuint> &attrib);
+  IGL_INLINE GLuint create_shader_program(
+    const std::string & vert_source,
+    const std::string & frag_source,
+    const std::map<std::string,GLuint> &attrib);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cross_field_missmatch.cpp

@@ -133,5 +133,6 @@ IGL_INLINE void igl::cross_field_missmatch(const Eigen::PlainObjectBase<DerivedV
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::cross_field_missmatch<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::cross_field_missmatch<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 
 #endif

include/igl/crouzeix_raviart_massmatrix.cpp

@@ -39,7 +39,6 @@ void igl::crouzeix_raviart_massmatrix(
     const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
     Eigen::SparseMatrix<MT> & M)
 {
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
   assert(F.cols() == 3);

include/igl/cut_mesh_from_singularities.cpp

@@ -198,4 +198,6 @@ IGL_INLINE void igl::cut_mesh_from_singularities(const Eigen::PlainObjectBase<De
 template void igl::cut_mesh_from_singularities<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::cut_mesh_from_singularities<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::cut_mesh_from_singularities<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::cut_mesh_from_singularities<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::cut_mesh_from_singularities<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/dated_copy.cpp

@@ -19,7 +19,6 @@
 IGL_INLINE bool igl::dated_copy(const std::string & src_path, const std::string & dir)
 {
   using namespace std;
-  using namespace igl;
   // Get time and date as string
   char buffer[80];
   time_t rawtime;

include/igl/deprecated.h

@@ -0,0 +1,21 @@
+#ifndef IGL_DEPRECATED_H
+#define IGL_DEPRECATED_H
+// Macro for marking a function as deprecated.
+// 
+// http://stackoverflow.com/a/295229/148668
+#ifdef __GNUC__
+#define IGL_DEPRECATED(func) func __attribute__ ((deprecated))
+#elif defined(_MSC_VER)
+#define IGL_DEPRECATED(func) __declspec(deprecated) func
+#else
+#pragma message("WARNING: You need to implement IGL_DEPRECATED for this compiler")
+#define IGL_DEPRECATED(func) func
+#endif
+// Usage:
+//
+//     template <typename T> IGL_INLINE void my_func(Arg1 a);
+//
+// becomes 
+//
+//     template <typename T> IGL_INLINE IGL_DEPRECATED(void my_func(Arg1 a));
+#endif

include/igl/dirname.h

@@ -13,7 +13,7 @@
 
 namespace igl
 {
-  // Function like PHP's dirname
+  // Function like PHP's dirname: /etc/passwd --> /etc, 
   // Input:
   //  path  string containing input path
   // Returns string containing dirname (see php's dirname)

include/igl/dot_row.cpp

@@ -22,4 +22,5 @@ IGL_INLINE Eigen::PlainObjectBase<DerivedV> igl::dot_row(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > igl::dot_row<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&);
 #endif

include/igl/doublearea.cpp

@@ -194,4 +194,5 @@ template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::M
 template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template Eigen::Matrix<double, 2, 1, 0, 2, 1>::Scalar igl::doublearea_single<Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&);
 template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/dqs.cpp

@@ -22,7 +22,6 @@ IGL_INLINE void igl::dqs(
   Eigen::PlainObjectBase<DerivedU> & U)
 {
   using namespace std;
-  using namespace igl;
   assert(V.rows() <= W.rows());
   assert(W.cols() == (int)vQ.size());
   assert(W.cols() == (int)vT.size());

include/igl/draw_mesh.h

@@ -20,7 +20,7 @@ namespace igl
   //
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V|#F by 3 eigen Matrix of 3D normals
   IGL_INLINE void draw_mesh(
     const Eigen::MatrixXd & V,
@@ -32,7 +32,7 @@ namespace igl
   //
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V|#F by 3 eigen Matrix of 3D normals
   //   C  #V|#F|1 by 3 eigen Matrix of RGB colors
   IGL_INLINE void draw_mesh(
@@ -42,7 +42,7 @@ namespace igl
     const Eigen::MatrixXd & C);
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V|#F by 3 eigen Matrix of 3D normals
   //   C  #V|#F|1 by 3 eigen Matrix of RGB colors
   //   TC  #V|#F|1 by 3 eigen Matrix of Texture Coordinates
@@ -58,7 +58,7 @@ namespace igl
   //
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V by 3 eigen Matrix of mesh vertex 3D normals
   //   C  #V by 3 eigen Matrix of mesh vertex RGB colors
   //   TC  #V by 3 eigen Matrix of mesh vertex UC coorindates between 0 and 1
@@ -84,14 +84,14 @@ namespace igl
   //
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V by 3 eigen Matrix of mesh vertex 3D normals
-  //   NF  #F by 3 eigen Matrix of face (triangle) normal indices, <0 means no
-  //     normal
+  //   NF  #F by 3 eigen Matrix of face (triangle/quad) normal indices, <0
+  //     means no normal
   //   C  #V by 3 eigen Matrix of mesh vertex RGB colors
   //   TC  #V by 3 eigen Matrix of mesh vertex UC coorindates between 0 and 1
-  //   TF  #F by 3 eigen Matrix of face (triangle) texture indices, <0 means no
-  //     texture
+  //   TF  #F by 3 eigen Matrix of face (triangle/quad) texture indices, <0
+  //     means no texture
   //   W  #V by #H eigen Matrix of per mesh vertex, per handle weights
   //   W_index  Specifies the index of the "weight" vertex attribute: see
   //     glBindAttribLocation, if W_index is 0 then weights are ignored

include/igl/draw_rectangular_marquee.cpp

@@ -15,7 +15,6 @@ IGL_INLINE void igl::draw_rectangular_marquee(
   const int to_x,
   const int to_y)
 {
-  using namespace igl;
   using namespace std;
   int l;
   glGetIntegerv(GL_LIGHTING,&l);

include/igl/edge_collapse_is_valid.cpp

@@ -23,7 +23,6 @@ IGL_INLINE bool igl::edge_collapse_is_valid(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   // For consistency with collapse_edge.cpp, let's determine edge flipness
   // (though not needed to check validity)
   const int eflip = E(e,0)>E(e,1);

include/igl/edge_lengths.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "edge_lengths.h"
 #include <iostream>
@@ -88,4 +88,10 @@ template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Ma
 template void igl::edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
 template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/embree/EmbreeIntersector.h

@@ -16,16 +16,15 @@
 #ifndef IGL_EMBREE_INTERSECTOR_H
 #define IGL_EMBREE_INTERSECTOR_H
 
+#include "Hit.h"
 #include <Eigen/Geometry>
 #include <Eigen/Core>
 #include <Eigen/Geometry>
 
-#include <vector>
 #include <embree2/rtcore.h>
 #include <embree2/rtcore_ray.h>
 #include <iostream>
-#include "Hit.h"
-#include <iostream>
+#include <vector>
 
 namespace igl
 {
@@ -279,7 +278,7 @@ inline void igl::EmbreeIntersector::init(
   // create a scene
   scene = rtcNewScene(RTC_SCENE_ROBUST | RTC_SCENE_HIGH_QUALITY,RTC_INTERSECT1);
 
-  for(int g=0;g<V.size();g++)
+  for(int g=0;g<(int)V.size();g++)
   {
     // create triangle mesh geometry in that scene
     geomID = rtcNewTriangleMesh(scene,RTC_GEOMETRY_STATIC,F[g]->rows(),V[g]->rows(),1);
@@ -328,14 +327,21 @@ igl::EmbreeIntersector
 
 void igl::EmbreeIntersector::deinit()
 {
-  rtcDeleteScene(scene);
+  if(scene)
+  {
+    rtcDeleteScene(scene);
 
-  if(rtcGetError() != RTC_NO_ERROR)
-      std::cerr << "Embree: An error occured while resetting!" << std::endl;
+    if(rtcGetError() != RTC_NO_ERROR)
+    {
+        std::cerr << "Embree: An error occured while resetting!" << std::endl;
+    }
 #ifdef IGL_VERBOSE
-  else
-    std::cerr << "Embree: geometry removed." << std::endl;
+    else
+    {
+      std::cerr << "Embree: geometry removed." << std::endl;
+    }
 #endif
+  }
 }
 
 inline bool igl::EmbreeIntersector::intersectRay(
@@ -347,7 +353,7 @@ inline bool igl::EmbreeIntersector::intersectRay(
   int mask) const
 {
   RTCRay ray;
-  createRay(ray, origin,direction,tnear,std::numeric_limits<float>::infinity(),mask);
+  createRay(ray, origin,direction,tnear,tfar,mask);
   
   // shot ray
   rtcIntersect(scene,ray);
@@ -356,7 +362,7 @@ inline bool igl::EmbreeIntersector::intersectRay(
       std::cerr << "Embree: An error occured while resetting!" << std::endl;
 #endif
   
-  if(ray.geomID != RTC_INVALID_GEOMETRY_ID)
+  if((unsigned)ray.geomID != RTC_INVALID_GEOMETRY_ID)
   {
     hit.id = ray.primID;
     hit.gid = ray.geomID;
@@ -388,8 +394,7 @@ inline bool igl::EmbreeIntersector::intersectBeam(
   else
     bestHit.t = 0;
 
-  hasHit = (intersectRay(origin,direction,hit,tnear,tfar,mask) && (hit.gid == geoId || geoId == -1));
-  if(hasHit)
+  if(hasHit = (intersectRay(origin,direction,hit,tnear,tfar,mask) && (hit.gid == geoId || geoId == -1)))
     bestHit = hit;
   
   // sample points around actual ray (conservative hitcheck)
@@ -448,7 +453,7 @@ igl::EmbreeIntersector
     ray.instID = RTC_INVALID_GEOMETRY_ID;
     num_rays++;
     rtcIntersect(scene,ray);
-    if(ray.geomID != RTC_INVALID_GEOMETRY_ID)
+    if((unsigned)ray.geomID != RTC_INVALID_GEOMETRY_ID)
     {
       // Hit self again, progressively advance
       if(ray.primID == last_id0 || ray.tfar <= min_t)
@@ -523,7 +528,7 @@ igl::EmbreeIntersector
   
   rtcIntersect(scene,ray);
 
-  if(ray.geomID != RTC_INVALID_GEOMETRY_ID)
+  if((unsigned)ray.geomID != RTC_INVALID_GEOMETRY_ID)
   {
     hit.id = ray.primID;
     hit.gid = ray.geomID;