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

Conflicts:
	build/Makefile_embree
	build/Makefile_png
	include/igl/cut_mesh_from_singularities.cpp

Former-commit-id: d4c7a9e29883a56fdb660f7ed6bf4d9f07ba3f8c

README.md

@@ -12,27 +12,32 @@ mesh-viewing utilities for OpenGL and GLSL, and many core functions for matrix
 manipulation which make [Eigen](http://eigen.tuxfamily.org) feel a lot more
 like MATLAB.
 
-It is first and foremost a header library. Each header file contains a single
-function. Most are tailored to operate on a generic triangle mesh stored in an
-n-by-3 matrix of vertex positions V and an m-by-3 matrix of triangle indices F.
-The library may also be [compiled](build/) into a statically linked
-library, for faster compile times with your projects.
+It is **a header-only library**. You do not need to compile anything to use,
+just include igl headers (e.g. `#include <igl/cotmatrix.h>`) and run.  Each
+header file contains a single function (e.g. `igl/cotmatrix.h` contains
+`igl::cotmatrix()`). Most are tailored to operate on a generic triangle mesh
+stored in an n-by-3 matrix of vertex positions V and an m-by-3 matrix of
+triangle indices F. 
+
+_Optionally_ the library may also be [pre-compiled](build/) into a statically
+linked library, for faster compile times with your projects. This only effects
+compile time (run-time performance and behavior is identical). If in doubt, use
+the header-only default mode: (i.e. just include the headers you want to use).
 
 We use the [Eigen](http://eigen.tuxfamily.org) library heavily in our code. Our
-group prototypes a lot in MATLAB, and we have a useful [conversion
-table](matlab-to-eigen.html) from
-MATLAB to libigl/Eigen.
+group prototypes a lot in MATLAB, and we have a useful [MATLAB to libigl+Eigen
+conversion table](matlab-to-eigen.html).
 
 ## Tutorial
 
 As of version 1.0, libigl includes an introductory
-[tutorial](tutorial/tutorial.html) that covers
-its basic functionalities.
+[tutorial](tutorial/tutorial.html) that covers many functionalities.
 
 ## Installation
-Libigl is a *header* library. You do **not** need to build anything to install.
-Simply add `igl/` to your include path and include relevant headers. Here is a
-small "Hello, World" program:
+
+Libigl is a **header-only** library. You do **not** need to build anything to
+install.  Simply add `libigl/include` to your include path and include relevant
+headers.  Here is a small "Hello, World" program:
 
 ```cpp
 #include <igl/cotmatrix.h>
@@ -57,10 +62,10 @@ int main()
 ```
 
 If you save this in `hello.cpp`, then you could compile this with (assuming
-Eigen is installed in /opt/local/include/eigen3):
+Eigen is installed in `/usr/local/include/eigen3`):
 
 ```bash
-gcc -I/opt/local/include/eigen3 -I./igl/ hello.cpp -o hello
+gcc -I/usr/local/include/eigen3 -I./libigl/include/ hello.cpp -o hello
 ```
 
 Running `./hello` would then produce
@@ -86,10 +91,9 @@ GCC 4.7 and clang will work correctly.
 
 ### OpenMP and Windows
 Some of our functions will take advantage of OpenMP if available. However, it
-has come to our attention that Visual Studio + Eigen does not work properly
-with OpenMP. Since OpenMP only improves performance without affecting
-functionality we recommend avoiding OpenMP on Windows or proceeding with
-caution.
+has come to our attention that Visual Studio + Eigen + OpenMP does not work
+properly. Since we use OpenMP only to improve performance, we recommend
+avoiding OpenMP on Windows or proceeding with caution.
 
 ## Download
 You can keep up to date by cloning a read-only copy of our GitHub
@@ -97,9 +101,9 @@ You can keep up to date by cloning a read-only copy of our GitHub
 
 ## Known Issues
 We really heavily on Eigen. Nearly all inputs and outputs are Eigen matrices of
-some kind. However, we currently _only_ support Eigen's default column-major
-ordering. That means, we **do not** expect our code to work for matrices using
-the `Eigen::RowMajor` flag. If you can, change definitions like:
+some kind. However, we currently _only_ officially support Eigen's default
+column-major ordering. That means, we **do not** expect our code to work for
+matrices using the `Eigen::RowMajor` flag. If you can, change definitions like:
 
 ```cpp
 Eigen::Matrix<double, Eigen::Dynamic, 3, Eigen::RowMajor> A;

RELEASE_HISTORY.txt

@@ -1,12 +1,18 @@
-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
-1.1.1  PLY file format support
-1.1.0  Mesh boolean operations using CGAL and cork, implementing [Attene 14]
-1.0.3  Bone heat method
-1.0.2  Bug fix in winding number code
-1.0.1  Bug fixes and more CGAL support
-1.0.0  Major beta release: many renames, tutorial, triangle wrapper, org. build
+# Version tracking
+
+Version | Short description
+--------|----------------------------------------------------------------------
+1.1.7   | Switch build for static library to cmake.
+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
+1.1.1   | PLY file format support
+1.1.0   | Mesh boolean operations using CGAL and cork, implementing [Attene 14]
+1.0.3   | Bone heat method
+1.0.2   | Bug fix in winding number code
+1.0.1   | Bug fixes and more CGAL support
+1.0.0   | Major beta release: many renames, tutorial, triangle, org. build
 
 ## Version 1.0 Changes ##
 Our beta release marks our confidence that this library can be used outside of

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.7

build/Makefile

@@ -1,127 +0,0 @@
-.PHONY: all
-all: lib extras examples
-framework: lib extras lib/igl.framework/
-
-# Shared flags etc.
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-$(info Hello, $(IGL_USERNAME)!)
-
-# optimized default settings
-all: LFLAGS +=
-#debug: OPTFLAGS= -g -Wall -Werror
-debug: OPTFLAGS= -g -Wall
-debug: DEBUG=debug
-CFLAGS += $(OPTFLAGS)
-#CFLAGS += -DIGL_NO_OPENGL -DIGL_NO_ANTTWEAKBAR
-# We use well-supported features of c++11
-
-EXTRAS=
-ifeq ($(IGL_WITH_BBW),1)
-	EXTRAS += bbw
-endif
-ifeq ($(IGL_WITH_BOOLEAN),1)
-	EXTRAS += boolean
-endif
-ifeq ($(IGL_WITH_BOOST),1)
-	EXTRAS += boost
-endif
-ifeq ($(IGL_WITH_CGAL),1)
-	EXTRAS += cgal
-endif
-ifeq ($(IGL_WITH_EMBREE),1)
-	EXTRAS += embree
-endif
-ifeq ($(IGL_WITH_COMISO),1)
-	EXTRAS += comiso
-endif
-ifeq ($(IGL_WITH_MATLAB),1)
-	EXTRAS += matlab
-endif
-ifeq ($(IGL_WITH_MOSEK),1)
-	EXTRAS += mosek
-endif
-ifeq ($(IGL_WITH_PNG),1)
-	EXTRAS += png
-endif
-ifeq ($(IGL_WITH_SVD3X3),1)
-	EXTRAS += svd3x3
-endif
-ifeq ($(IGL_WITH_TETGEN),1)
-	# append tetgen extra dir
-	EXTRAS += tetgen
-endif
-ifeq ($(IGL_WITH_VIEWER),1)
-	EXTRAS += viewer
-endif
-ifeq ($(IGL_WITH_XML),1)
-	EXTRAS += xml
-endif
-
-.PHONY: examples
-.PHONY: extras
-debug: lib extras
-lib: ../lib/libigl.a
-examples: lib extras
-	make -C ../examples
-extras:
-	for p in  $(EXTRAS); \
-	do \
-	$(MAKE) -f Makefile_$$p $(DEBUG); \
-	done
-
-
-#############################################################################
-# SOURCE
-#############################################################################
-CPP_FILES=$(wildcard ../include/igl/*.cpp)
-H_FILES=$(wildcard ../include/igl/*.h)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-#############################################################################
-# DEPENDENCIES
-#############################################################################
-INC+=$(OPENGL_INC)
-INC+=$(EIGEN3_INC)
-INC+=$(ANTTWEAKBAR_INC)
-
-.PHONY: obj
-obj:
-	mkdir -p obj
-
-../lib/libigl.a: obj $(OBJ_FILES)
-	mkdir -p ../lib
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: ../include/igl/%.cpp ../include/igl/%.h
-	$(GG) $(CFLAGS) $(AFLAGS) -c -o $@ $< $(INC)
-
-#lib/igl.framework/:
-#	mkdir -p $@
-#	cp lib/*.a $@
-#	mv $@/libigl.a $@/igl
-#	mkdir -p $@/Libraries
-#	mv $@/*.a $@/Libraries
-#	mkdir -p $@/Headers
-#	cp $(H_FILES) $@/Headers
-#	for p in $(EXTRAS); \
-#	do \
-#	mkdir $@/Headers/$$p; \
-#	cp include/igl/$$p/*.h $@/Headers/$$p; \
-#	done
-
-
-clean:
-	rm -rf ../lib/igl.framework/
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libigl.a
-	make -C ../examples clean
-	for p in  $(EXTRA_DIRS); \
-	do \
-	echo "cd $$p" ; \
-	$(MAKE) -C $$p clean; \
-	done

build/Makefile_bbw

@@ -1,59 +0,0 @@
-.PHONY: all
-all: libiglbbw
-debug: libiglbbw
-
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: OPTFLAGS += -O3 -DNDEBUG $(OPENMP)
-debug: OPTFLAGS += -g -Wall
-CFLAGS += $(OPTFLAGS)
-
-.PHONY: libiglbbw
-libiglbbw: obj ../lib/libiglbbw.a
-
-SRC_DIR=../include/igl/bbw/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-ifdef IGL_NO_MOSEK
-CFLAGS+=-DIGL_NO_MOSEK
-else
-# mosek dependency
-# TODO: linux, 32 bit etc
-ifndef MOSEKPLATFORM
-  MOSEKPLATFORM=osx64x86
-endif
-ifndef MOSEKVERSION
-  MOSEKVERSION=6
-endif
-IGLMOSEK=../mosek/
-IGLMOSEK_INC=-I$(IGLMOSEK)/
-INC+=${IGLMOSEK_INC}
-MOSEK=/usr/local/mosek
-MOSEK_INC=-I$(MOSEK)/$(MOSEKVERSION)/tools/platform/$(MOSEKPLATFORM)/h
-MOSEK_LIB=-L$(MOSEK)/$(MOSEKVERSION)/tools/platform/$(MOSEKPLATFORM)/bin -lmosek64
-INC+=$(MOSEK_INC)
-endif
-
-obj: 
-	mkdir -p obj
-
-../lib/libiglbbw.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	$(GG) $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglbbw.a

build/Makefile_boolean

@@ -1,54 +0,0 @@
-include Makefile.conf
-
-.PHONY: all
-all: libiglboolean
-debug: libiglboolean
-
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG -std=c++11
-debug: CFLAGS += -g -Wall -std=c++11
-
-.PHONY: libiglboolean
-libiglboolean: obj ../lib/libiglboolean.a
-
-SRC_DIR=../include/igl/boolean/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# cork dependency
-#CORK=../external/cork
-#CORK_INC=-I$(cork)/common -I$(cork)/rtcore/
-CORK=../external/cork/
-CORK_INC=-I$(CORK)/include/
-INC+=$(CORK_INC)
-
-# CGAL dependency
-CGAL=$(DEFAULT_PREFIX)
-CGAL_INC=-I$(CGAL)/include
-CGAL_FLAGS=-frounding-math -fsignaling-nans 
-CFLAGS+=$(CGAL_FLAGS)
-INC+=$(CGAL_INC)
-
-obj: 
-	mkdir -p obj
-
-../lib/libiglboolean.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	g++ $(AFLAGS) $(OPENMP) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglboolean.a

build/Makefile_boost

@@ -1,45 +0,0 @@
-include Makefile.conf
-
-.PHONY: all
-all: libiglboost
-debug: libiglboost
-
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG
-debug: CFLAGS += -g -Wall 
-
-.PHONY: libiglboost
-libiglboost: obj ../lib/libiglboost.a
-
-SRC_DIR=../include/igl/boost/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# boost dependency
-BOOST=/opt/local/
-BOOST_INC=-I$(BOOST)/include
-INC+=$(BOOST_INC)
-
-obj: 
-	mkdir -p obj
-
-../lib/libiglboost.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	g++ $(AFLAGS) $(OPENMP) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglboost.a

build/Makefile_cgal

@@ -1,46 +0,0 @@
-
-.PHONY: all
-all: libiglcgal
-debug: libiglcgal
-
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG 
-debug: CFLAGS += -g -Wall -Werror
-
-.PHONY: libcgal
-libiglcgal: obj ../lib/libiglcgal.a
-
-SRC_DIR=../include/igl/cgal/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# CGAL dependency
-CGAL=$(DEFAULT_PREFIX)
-CGAL_INC=-I$(CGAL)/include
-CGAL_FLAGS=-frounding-math -fsignaling-nans ${OPENMP}
-CFLAGS+=$(CGAL_FLAGS)
-INC+=$(CGAL_INC)
-
-obj: 
-	mkdir -p obj
-
-../lib/libiglcgal.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglcgal.a

build/Makefile_comiso

@@ -1,45 +0,0 @@
-include Makefile.conf
-
-.PHONY: all
-all: libiglcomiso
-debug: libiglcomiso
-
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG -std=c++11
-debug: CFLAGS += -g -Wall -std=c++11
-CFLAGS += -DINCLUDE_TEMPLATES
-
-.PHONY: libiglcomiso
-libiglcomiso: obj ../lib/libiglcomiso.a
-
-SRC_DIR=../include/igl/comiso/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# comiso dependency
-COMISO_INC=-I$(COMISO)/ -I$(COMISO)/gmm/include -I$(COMISO)/../
-INC+=$(COMISO_INC)
-
-obj:
-	mkdir -p obj
-
-../lib/libiglcomiso.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	g++ $(AFLAGS) $(OPENMP) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglcomiso.a

build/Makefile_lim

@@ -1,45 +0,0 @@
-include Makefile.conf
-
-.PHONY: all
-all: libigllim
-debug: libigllim
-
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG
-debug: CFLAGS += -g -Wall 
-
-.PHONY: libigllim
-libigllim: obj ../lib/libigllim.a
-
-SRC_DIR=../include/igl/lim/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# LIM dependency
-LIM=../external/lim
-LIM_INC=-I$(LIM)
-INC+=$(LIM_INC)
-
-obj: 
-	mkdir -p obj
-
-../lib/libigllim.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	g++ $(AFLAGS) $(OPENMP) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libigllim.a

build/Makefile_matlab

@@ -1,47 +0,0 @@
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG
-debug: CFLAGS += -g -Wall -Werror
-
-.PHONY: all
-all: libiglmatlab
-debug: libiglmatlab
-
-.PHONY: libmatlab
-libiglmatlab: obj ../lib/libiglmatlab.a
-
-SRC_DIR=../include/igl/matlab/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# Matlab dependency
-ifndef MATLAB
-  MATLAB=/Applications/MATLAB_R2014b.app/
-  $(warning MATLAB undefined. Setting MATLAB=${MATLAB})
-endif
-MATLAB_INC=-I$(MATLAB)/extern/include/
-MATLAB_LIB=-L$(MATLAB)/bin/maci64 -lmx -leng
-INC+=$(MATLAB_INC)
-
-obj:
-	mkdir -p obj
-
-../lib/libiglmatlab.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglmatlab.a

build/Makefile_mosek

@@ -1,54 +0,0 @@
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-
-.PHONY: all
-all: libiglmosek
-debug: libiglmosek
-
-include Makefile.conf
-all: OPTFLAGS += -O3 -DNDEBUG $(OPENMP)
-debug: OPTFLAGS += -g -Wall
-CFLAGS += $(OPTFLAGS)
-
-.PHONY: libiglmosek
-libiglmosek: obj ../lib/libiglmosek.a
-
-SRC_DIR=../include/igl/mosek/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# mosek dependency
-# TODO: linux, 32 bit etc
-ifndef MOSEKPLATFORM
-  MOSEKPLATFORM=osx64x86
-endif
-ifndef MOSEKVERSION
-  MOSEKVERSION=6
-endif
-MOSEK=/usr/local/mosek
-MOSEK_INC=-I$(MOSEK)/$(MOSEKVERSION)/tools/platform/$(MOSEKPLATFORM)/h
-MOSEK_LIB=-L$(MOSEK)/$(MOSEKVERSION)/tools/platform/$(MOSEKPLATFORM)/bin -lmosek64
-INC+=$(MOSEK_INC)
-
-obj: 
-	mkdir -p obj
-
-../lib/libiglmosek.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglmosek.a

build/Makefile_svd3x3

@@ -1,53 +0,0 @@
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-
-.PHONY: all
-all: libiglsvd3x3
-debug: libiglsvd3x3
-
-include Makefile.conf
-all: OPTFLAGS += -O3 -DNDEBUG $(OPENMP)
-debug: OPTFLAGS += -g -Wall
-CFLAGS += $(OPTFLAGS) -std=c++11
-
-.PHONY: libiglsvd3x3
-libiglsvd3x3: obj ../lib/libiglsvd3x3.a
-
-SRC_DIR=../include/igl/svd3x3/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
-
-# SVD 
-SINGULAR_VALUE_DECOMPOSITION_INC=\
-	-I../external/Singular_Value_Decomposition/
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC) $(SINGULAR_VALUE_DECOMPOSITION_INC)
-
-# Matlab dependency
-ifndef MATLAB
-  MATLAB=/Applications/MATLAB_R2013b.app/
-  $(warning MATLAB undefined. Setting MATLAB=${MATLAB})
-endif
-MATLAB_INC=-I$(MATLAB)/extern/include/
-MATLAB_LIB=-L$(MATLAB)/bin/maci64 -lmx -leng
-INC+=$(MATLAB_INC)
-
-obj: 
-	mkdir -p obj
-
-../lib/libiglsvd3x3.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	$(GG) $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglsvd3x3.a

build/Makefile_tetgen

@@ -1,45 +0,0 @@
-
-.PHONY: all
-all: libigltetgen
-debug: libigltetgen
-
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG 
-debug: CFLAGS += -g -Wall -Werror
-
-.PHONY: libtetgen
-libigltetgen: obj ../lib/libigltetgen.a
-
-SRC_DIR=../include/igl/tetgen/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# Tetgen dependency
-TETGEN=../external/tetgen
-TETGEN_INC=-I$(TETGEN)
-INC+=$(TETGEN_INC)
-TETGEN_STATIC_LIB=$(TETGEN)/libtet.a
-
-obj: 
-	mkdir -p obj
-
-../lib/libigltetgen.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libigltetgen.a

build/Makefile_triangle

@@ -1,45 +0,0 @@
-
-.PHONY: all
-all: libigltriangle
-debug: libigltriangle
-
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG 
-debug: CFLAGS += -g -Wall -Werror
-
-.PHONY: libtriangle
-libigltriangle: obj ../lib/libigltriangle.a
-
-SRC_DIR=../include/igl/triangle/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# triangle dependency
-TRIANGLE=../external/triangle
-TRIANGLE_INC=-I$(TRIANGLE)
-INC+=$(TRIANGLE_INC)
-TRIANGLE_STATIC_LIB=$(TRIANGLE)/libtriangle.a
-
-obj: 
-	mkdir -p obj
-
-../lib/libigltriangle.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libigltriangle.a

build/Makefile_viewer

@@ -1,51 +0,0 @@
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG -fopenmp
-debug: CFLAGS += -g -Wall -Werror -fopenmp
-
-.PHONY: all
-all: libiglviewer
-debug: libiglviewer
-
-.PHONY: libviewer
-libiglviewer: obj ../lib/libiglviewer.a
-
-SRC_DIR=../include/igl/viewer/
-CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY
-
-# Eigen dependency
-EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-INC+=$(EIGEN3_INC)
-
-# GLFW dependency
-ifndef GLFW
-	GLFW=../external/glfw/
-endif
-GLFW_INC=-I$(GLFW)/include
-INC+=$(GLFW_INC)
-ifndef ANTTWEAKBAR
-	ANTTWEAKBAR=../external/AntTweakBar/
-endif
-# Viewer also uses source files from AntTweakBar for font rendering
-ANTTWEAKBAR_INC=-I$(ANTTWEAKBAR)/include -I$(ANTTWEAKBAR)/src
-INC+=$(ANTTWEAKBAR_INC)
-
-obj:
-	mkdir -p obj
-
-../lib/libiglviewer.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: $(SRC_DIR)/%.cpp $(SRC_DIR)/%.h
-	$(GG) $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglmatlab.a

build/Makefile_xml

@@ -1,61 +0,0 @@
-CFLAGS += -DIGL_STATIC_LIBRARY
-include Makefile.conf
-
-.PHONY: all
-#all:
-#debug:
-all: libiglxml
-debug: libiglxml
-
-include Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG
-debug: CFLAGS += -g -Wall
-
-.PHONY: libiglxml
-libiglxml: obj ../lib/libiglxml.a
-
-#SRC_DIR=../include/igl/xml/
-#CPP_FILES=$(wildcard $(SRC_DIR)*.cpp)
-CPP_FILES=$(wildcard ../include/igl/xml/*.cpp)
-OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
-
-# include igl headers
-INC+=-I../include/
-
-# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY
-
-# Eigen dependency
-ifndef EIGEN3_INC
-	EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
-endif
-INC+=$(EIGEN3_INC)
-
-#AntTweakbar dependency
-ANTTWEAKBAR_INC=-I../external/AntTweakBar/include
-
-# xml dependency
-# TODO: linux, 32 bit etc
-TINYXML2=../external/tinyxml2
-TINYXML2_INC=-I$(TINYXML2)
-#TINYXML2_LIB=-L$(TINYXML2) -ltinyxml2
-INC+=$(TINYXML2_INC) $(ANTTWEAKBAR_INC)
-
-# AntTweakBar dependency
-ANTTWEAKBAR_INC=-I../external/AntTweakBar/include
-INC+=$(ANTTWEAKBAR_INC)
-
-#CFLAGS+=-std=c++11
-
-obj:
-	mkdir -p obj
-
-../lib/libiglxml.a: $(OBJ_FILES)
-	rm -f $@
-	ar cqs $@ $(OBJ_FILES)
-
-obj/%.o: ../include/igl/xml/%.cpp 
-	$(GG) $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
-
-clean:
-	rm -f $(OBJ_FILES)
-	rm -f ../lib/libiglxml.a

examples/MatlabWorkspace/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 LIBIGL_LIB+=-liglmatlab
 
 MATLAB_INC=-I$(MATLAB)/extern/include/

examples/ReAntTweakBar/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/affine/Makefile

@@ -3,7 +3,7 @@
 all: example
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 .PHONY: example
 

examples/ambient-occlusion/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 LIBIGL_LIB+=-liglembree
 
 all: example

examples/arap/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 LIBIGL_LIB+=-liglsvd3x3
 
 all: example

examples/bbw/Makefile

@@ -1,4 +1,4 @@
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 CXX=g++
 

examples/beach-balls/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: obj example
 

examples/camera/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: obj example
 

examples/colored-mesh/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/components/Makefile

@@ -0,0 +1,33 @@
+.PHONY: all
+
+# Shared flags etc.
+include ../Makefile.conf
+LIBIGL_LIB+=-liglembree -liglboost
+
+all: obj example
+
+.PHONY: example
+
+CFLAGS+=-g -std=c++11
+
+INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(GLUT_INC) $(EMBREE_INC)
+LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(EMBREE_LIB)
+
+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)
+
+obj:
+	mkdir -p obj
+
+obj/%.o: %.cpp
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
+
+obj/%.o: %.cpp %.h
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
+
+clean:
+	rm -f $(OBJ_FILES)
+	rm -f example

examples/components/example.cpp

@@ -0,0 +1,1095 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/init_render_to_texture.h>
+#include <igl/draw_floor.h>
+#include <igl/report_gl_error.h>
+#include <igl/per_face_normals.h>
+#include <igl/trackball.h>
+#include <igl/snap_to_canonical_view_quat.h>
+#include <igl/REDRUM.h>
+#include <igl/Camera.h>
+#include <igl/ReAntTweakBar.h>
+#include <igl/get_seconds.h>
+#include <igl/jet.h>
+#include <igl/rgb_to_hsv.h>
+#include <igl/hsv_to_rgb.h>
+#include <igl/randperm.h>
+#include <igl/boost/components.h>
+#include <igl/C_STR.h>
+#include <igl/write_triangle_mesh.h>
+#include <igl/two_axis_valuator_fixed_up.h>
+#include <igl/snap_to_fixed_up.h>
+#include <igl/create_shader_program.h>
+
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+#ifdef __APPLE__
+#include <GLUT/glut.h>
+#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 <ctime>
+#include <string>
+#include <vector>
+#include <stack>
+#include <iostream>
+
+int cc_hover = -1;
+
+Eigen::MatrixXd V;
+Eigen::VectorXd Vmid,Vmin,Vmax;
+double bbd = 1.0;
+Eigen::MatrixXi F;
+Eigen::VectorXi CC;
+Eigen::MatrixXd N;
+struct State
+{
+  igl::Camera camera;
+  Eigen::VectorXf I;
+  Eigen::Matrix<GLubyte,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> selected;
+  GLuint mask_id;
+} s;
+std::string out_filename;
+
+GLuint pick_tex = 0;
+GLuint pick_fbo = 0;
+GLuint pick_dfbo = 0;
+
+// See README for descriptions
+enum RotationType
+{
+  ROTATION_TYPE_IGL_TRACKBALL = 0,
+  ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP = 1,
+  NUM_ROTATION_TYPES = 2,
+} rotation_type;
+
+enum CenterType
+{
+  CENTER_TYPE_ORBIT = 0,
+  CENTER_TYPE_FPS  = 1,
+  NUM_CENTER_TYPES = 2,
+} center_type = CENTER_TYPE_ORBIT;
+
+
+std::stack<State> undo_stack;
+std::stack<State> redo_stack;
+
+bool wireframe_visible = false;
+bool fill_visible = true;
+
+bool is_rotating = false;
+int down_x,down_y;
+igl::Camera down_camera;
+
+bool is_animating = false;
+double animation_start_time = 0;
+double ANIMATION_DURATION = 0.5;
+Eigen::Quaterniond animation_from_quat;
+Eigen::Quaterniond animation_to_quat;
+
+int width,height;
+Eigen::Vector4f light_pos(-0.1,-0.1,0.9,0);
+
+#define REBAR_NAME "temp.rbr"
+igl::ReTwBar rebar;
+
+// Forward
+void init_components();
+void init_relative();
+
+void push_undo()
+{
+  undo_stack.push(s);
+  // Clear
+  redo_stack = std::stack<State>();
+}
+
+void TW_CALL set_rotation_type(const void * value, void * clientData)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  const RotationType old_rotation_type = rotation_type;
+  rotation_type = *(const RotationType *)(value);
+  if(rotation_type == ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP &&
+    old_rotation_type != ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP)
+  {
+    animation_from_quat = s.camera.m_rotation_conj;
+    snap_to_fixed_up(animation_from_quat,animation_to_quat);
+    // start animation
+    animation_start_time = get_seconds();
+    is_animating = true;
+  }
+}
+void TW_CALL get_rotation_type(void * value, void *clientData)
+{
+  RotationType * rt = (RotationType *)(value);
+  *rt = rotation_type;
+}
+
+void reshape(int width, int height)
+{
+  ::width = width;
+  ::height = height;
+  glViewport(0,0,width,height);
+  // Send the new window size to AntTweakBar
+  TwWindowSize(width, height);
+  s.camera.m_aspect = (double)width/(double)height;
+  igl::init_render_to_texture(width,height, pick_tex, pick_fbo, pick_dfbo);
+  igl::report_gl_error("init_render_to_texture: ");
+  glutPostRedisplay();
+}
+
+void push_scene()
+{
+  using namespace igl;
+  using namespace std;
+  glMatrixMode(GL_PROJECTION);
+  glPushMatrix();
+  glLoadIdentity();
+  auto & camera = s.camera;
+  glMultMatrixd(camera.projection().data());
+  glMatrixMode(GL_MODELVIEW);
+  glPushMatrix();
+  glLoadIdentity();
+  gluLookAt(
+    camera.eye()(0), camera.eye()(1), camera.eye()(2),
+    camera.at()(0), camera.at()(1), camera.at()(2),
+    camera.up()(0), camera.up()(1), camera.up()(2));
+  glScaled(2./bbd,2./bbd,2./bbd);
+  glTranslated(-Vmid(0),-Vmid(1),-Vmid(2));
+}
+
+void pop_scene()
+{
+  glMatrixMode(GL_PROJECTION);
+  glPopMatrix();
+  glMatrixMode(GL_MODELVIEW);
+  glPopMatrix();
+}
+
+void draw_mesh(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & N,
+  const Eigen::VectorXf & S,
+  const GLuint & S_loc)
+{
+  using namespace Eigen;
+  using namespace std;
+  static Matrix<float,Dynamic,3,RowMajor> VR,NR;
+  static Matrix<int,Dynamic,3,RowMajor> FR;
+  static Matrix<float,Dynamic,1,ColMajor> SR;
+  static GLuint ibo,vbo,sbo,nbo;
+  static bool scene_dirty = true;
+  if(scene_dirty)
+  {
+    VR.resize(F.rows()*3,3);
+    NR.resize(F.rows()*3,3);
+    SR.resize(F.rows()*3,1);
+    FR.resize(F.rows(),3);
+    for(int f = 0;f<F.rows();f++)
+    {
+      for(int c = 0;c<3;c++)
+      {
+        VR.row(3*f+c) = V.row(F(f,c)).cast<float>();
+        SR(3*f+c) = S(F(f,c));
+        NR.row(3*f+c) = N.row(f).cast<float>();
+        FR(f,c) = 3*f+c;
+      }
+    }
+
+    glGenBuffers(1,&ibo);
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER,ibo);
+    glBufferData(GL_ELEMENT_ARRAY_BUFFER,sizeof(GLuint)*FR.size(),FR.data(),GL_STATIC_DRAW);
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER,0);
+    glGenBuffers(1,&vbo);
+    glGenBuffers(1,&nbo);
+    glGenBuffers(1,&sbo);
+
+    glBindBuffer(GL_ARRAY_BUFFER,vbo);
+    glBufferData(GL_ARRAY_BUFFER,sizeof(float)*VR.size(),VR.data(),GL_STATIC_DRAW);
+    glBindBuffer(GL_ARRAY_BUFFER,nbo);
+    glBufferData(GL_ARRAY_BUFFER,sizeof(float)*NR.size(),NR.data(),GL_STATIC_DRAW);
+
+    glBindBuffer(GL_ARRAY_BUFFER,sbo);
+    glBufferData(GL_ARRAY_BUFFER,sizeof(float)*SR.size(),SR.data(),GL_STATIC_DRAW);
+    igl::report_gl_error("glBindBuffer: ");
+
+    scene_dirty = false;
+  }
+
+  glEnableClientState(GL_VERTEX_ARRAY);
+  glBindBuffer(GL_ARRAY_BUFFER,vbo);
+  glVertexPointer(3,GL_FLOAT,0,0);
+  glEnableClientState(GL_NORMAL_ARRAY);
+  glBindBuffer(GL_ARRAY_BUFFER,nbo);   
+  glNormalPointer(GL_FLOAT,0,0);
+
+  glBindBuffer(GL_ARRAY_BUFFER,sbo);   
+  glVertexAttribPointer(S_loc, 1, GL_FLOAT, GL_FALSE, 0, 0);
+  glEnableVertexAttribArray(S_loc);
+
+  glBindBuffer(GL_ELEMENT_ARRAY_BUFFER,ibo);
+  glDrawElements(GL_TRIANGLES,FR.size(),GL_UNSIGNED_INT,0);
+  glBindBuffer(GL_ARRAY_BUFFER,0);
+}
+
+// Set up double-sided lights
+void lights()
+{
+  using namespace std;
+  using namespace Eigen;
+  glEnable(GL_LIGHTING);
+  glLightModelf(GL_LIGHT_MODEL_TWO_SIDE,GL_TRUE);
+  glEnable(GL_LIGHT0);
+  float WHITE[4] = {1,1,1,1.};
+  float BLACK[4] = {0.,0.,0.,1.};
+  Vector4f pos = light_pos;
+  glLightfv(GL_LIGHT0,GL_AMBIENT,BLACK);
+  glLightfv(GL_LIGHT0,GL_DIFFUSE,WHITE);
+  glLightfv(GL_LIGHT0,GL_SPECULAR,BLACK);
+  glLightfv(GL_LIGHT0,GL_POSITION,pos.data());
+  //glEnable(GL_LIGHT1);
+  //pos(0) *= -1;
+  //pos(1) *= -1;
+  //pos(2) *= -1;
+  //glLightfv(GL_LIGHT1,GL_AMBIENT,BLACK);
+  //glLightfv(GL_LIGHT1,GL_DIFFUSE,NEAR_BLACK);
+  //glLightfv(GL_LIGHT1,GL_SPECULAR,BLACK);
+  //glLightfv(GL_LIGHT1,GL_POSITION,pos.data());
+}
+
+template <int Rows, int Cols>
+GLuint generate_1d_texture(
+  const Eigen::Matrix<GLubyte,Rows,Cols,Eigen::RowMajor> & colors)
+{
+  assert(colors.cols() == 3 && "Seems colors.cols() must be 3");
+  GLuint tex_id = 0;
+  glGenTextures(1,&tex_id);
+  glBindTexture(GL_TEXTURE_1D,tex_id);
+  glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
+  glTexImage1D(GL_TEXTURE_1D, 0, colors.cols(),colors.rows(),
+    0,GL_RGB, GL_UNSIGNED_BYTE,
+    colors.data());
+  igl::report_gl_error("glTexImage1D: ");
+  glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP);
+  glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_T, GL_CLAMP);
+  glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER,GL_NEAREST);
+  igl::report_gl_error("texture: ");
+  return tex_id;
+}
+ 
+GLuint color_shader(const size_t max_ids, GLuint & scalar_loc, GLuint & tex_id)
+{
+  std::string vertex_shader = R"(
+#version 120
+attribute float scalar_in;
+varying float scalar_out;
+void main()
+{
+  scalar_out = scalar_in;
+  gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;
+}
+)";
+  std::string fragment_shader = R"(
+#version 120
+varying float scalar_out;
+uniform float cmin;
+uniform float cmax;
+uniform sampler1D color_map;
+void main()
+{
+  float scalar_normalized = max(min((scalar_out-cmin)/(cmax-cmin),1.0),0.0);
+  gl_FragColor = texture1D(color_map,scalar_normalized);
+}
+)";
+  Eigen::Matrix<GLubyte,Eigen::Dynamic,3,Eigen::RowMajor> colors(max_ids,3);
+  for(size_t id = 0;id<max_ids;id++)
+  {
+    size_t index = id;
+    size_t re = (index)%(256*256);
+    colors(id,0) = (index-re)/(256*256);
+    index = re;
+    re = index%(256);
+    colors(id,1) = (index-re)/(256);
+    colors(id,2) = re;
+  }
+  tex_id = generate_1d_texture(colors);
+  return igl::create_shader_program(
+    vertex_shader.c_str(), 
+    fragment_shader.c_str(),
+    {{"scalar_in",scalar_loc}}
+    );
+}
+
+
+void display()
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  glClearColor(0.8,0.8,0.8,0);
+  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+  if(is_animating)
+  {
+    double t = (get_seconds() - animation_start_time)/ANIMATION_DURATION;
+    if(t > 1)
+    {
+      t = 1;
+      is_animating = false;
+    }
+    Quaterniond q = animation_from_quat.slerp(t,animation_to_quat).normalized();
+    auto & camera = s.camera;
+    switch(center_type)
+    {
+      default:
+      case CENTER_TYPE_ORBIT:
+        camera.orbit(q.conjugate());
+        break;
+      case CENTER_TYPE_FPS:
+        camera.turn_eye(q.conjugate());
+        break;
+    }
+  }
+
+  glEnable(GL_DEPTH_TEST);
+  glEnable(GL_NORMALIZE);
+  lights();
+  push_scene();
+
+
+  const auto & color_components_shader = [](
+     const GLuint scalar_loc,
+     GLuint & tex_id)->GLuint
+  {
+  std::string vertex_shader = R"(
+#version 120
+attribute float scalar_in;
+varying vec3 normal;
+varying float scalar_out;
+void main()
+{
+  gl_Position = gl_ModelViewProjectionMatrix * gl_Vertex;
+  normal = normalize(gl_NormalMatrix * gl_Normal);
+  scalar_out = scalar_in;
+}
+)";
+  std::string fragment_shader = R"(
+#version 120
+varying vec3 normal;
+varying float scalar_out;
+uniform float cmin;
+uniform float cmax;
+uniform float cc_hover;
+uniform sampler1D color_map;
+uniform sampler1D selected_mask;
+void main()
+{
+  float scalar_normalized = max(min((scalar_out-cmin)/(cmax-cmin),1.0),0.0);
+  vec4 texture_color = texture1D(color_map,scalar_normalized);
+  bool is_selected = texture1D(selected_mask,scalar_normalized).x > 0.5;
+  const vec4 selected_color = vec4(1,0.2,0.2,1);
+  if(scalar_out==cc_hover)
+  {
+    texture_color = 0.5*(texture_color + selected_color);
+  }
+  if(is_selected)
+  {
+    texture_color = selected_color;
+  }
+  const float num_lights = 1.0;
+  vec4 diffuse = (1.0/num_lights)*(gl_LightSource[0].diffuse);
+  vec4 ambient = vec4(0,0,0,0);
+  ambient += (1.0/num_lights)*(gl_FrontMaterial.ambient * gl_LightSource[0].ambient);
+  ambient += (1.0/num_lights)*(gl_LightModel.ambient * gl_FrontMaterial.ambient);
+  vec4 color = ambient;
+  // Phong
+  vec3 lightDir = normalize(vec3(gl_LightSource[0].position));
+  vec3 halfVector = gl_LightSource[0].halfVector.xyz;
+  vec3 n = normalize(normal);
+  float NdotL = max(abs(dot(n.xyz,lightDir)), 0.0);  
+  vec4 specular = vec4(0.0,0.0,0.0,0.0);
+  if (NdotL > 0.0) {
+      color += diffuse * NdotL;
+      vec3 halfV = normalize(halfVector);
+      float NdotHV = max(abs(dot(n,halfV)),0.0);
+      specular += gl_FrontMaterial.specular * gl_LightSource[0].specular * pow(NdotHV, gl_FrontMaterial.shininess);
+  }
+  gl_FragColor = color * texture_color + specular;
+}
+
+)";
+
+    typedef Matrix<GLubyte,64,3,RowMajor> Matrix64_3_R_ubyte;
+    typedef Matrix<float,64,3,RowMajor> Matrix64_3_R_float;
+    Matrix64_3_R_ubyte colors;
+    {
+      Matrix64_3_R_float rgb = (Matrix64_3_R_ubyte()<<
+        255,   0,   0,
+        255,  24,   0,
+        255,  48,   0,
+        255,  72,   0,
+        255,  96,   0,
+        255, 120,   0,
+        255, 143,   0,
+        255, 167,   0,
+        255, 191,   0,
+        255, 215,   0,
+        255, 239,   0,
+        247, 255,   0,
+        223, 255,   0,
+        199, 255,   0,
+        175, 255,   0,
+        151, 255,   0,
+        128, 255,   0,
+        104, 255,   0,
+         80, 255,   0,
+         56, 255,   0,
+         32, 255,   0,
+          8, 255,   0,
+          0, 255,  16,
+          0, 255,  40,
+          0, 255,  64,
+          0, 255,  88,
+          0, 255, 112,
+          0, 255, 135,
+          0, 255, 159,
+          0, 255, 183,
+          0, 255, 207,
+          0, 255, 231,
+          0, 255, 255,
+          0, 231, 255,
+          0, 207, 255,
+          0, 183, 255,
+          0, 159, 255,
+          0, 135, 255,
+          0, 112, 255,
+          0,  88, 255,
+          0,  64, 255,
+          0,  40, 255,
+          0,  16, 255,
+          8,   0, 255,
+         32,   0, 255,
+         56,   0, 255,
+         80,   0, 255,
+        104,   0, 255,
+        128,   0, 255,
+        151,   0, 255,
+        175,   0, 255,
+        199,   0, 255,
+        223,   0, 255,
+        247,   0, 255,
+        255,   0, 239,
+        255,   0, 215,
+        255,   0, 191,
+        255,   0, 167,
+        255,   0, 143,
+        255,   0, 120,
+        255,   0,  96,
+        255,   0,  72,
+        255,   0,  48,
+        255,   0,  24).finished().cast<float>()/255.f;
+
+      Matrix64_3_R_float H;
+      rgb_to_hsv(rgb,H);
+      H.col(1) *= 0.1;
+      H.col(2) = (H.col(2).array() + 0.1*(1.-H.col(2).array())).eval();
+      hsv_to_rgb(H,rgb);
+      colors = (rgb*255.).cast<GLubyte>();
+    }
+
+    tex_id = generate_1d_texture(colors);
+
+    GLuint prog_id = igl::create_shader_program(
+      vertex_shader.c_str(), 
+      fragment_shader.c_str(),
+      {{"scalar_in",scalar_loc}}
+      );
+    igl::report_gl_error("create_shader_program: ");
+    return prog_id;
+  };
+  static GLuint scalar_loc = 1;
+  static GLuint tex_id = 0;
+  static GLuint color_components_prog = 
+    color_components_shader(scalar_loc,tex_id);
+
+  // Set material properties
+  glEnable(GL_COLOR_MATERIAL);
+  glColorMaterial(GL_FRONT_AND_BACK,GL_AMBIENT_AND_DIFFUSE);
+  glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,(const GLfloat[]){1,1,1,1});
+  if(wireframe_visible)
+  {
+    glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
+    if(fill_visible)
+    {
+      glColor3f(0,0,0);
+      glUseProgram(0);
+      draw_mesh(V,F,N,s.I,scalar_loc);
+    }else
+    {
+      glUseProgram(color_components_prog);
+      igl::report_gl_error("UseProgram: ");
+      draw_mesh(V,F,N,s.I,scalar_loc);
+    }
+  }
+
+  glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
+
+  glPushAttrib(GL_ALL_ATTRIB_BITS);
+  glUseProgram(color_components_prog);
+    igl::report_gl_error("use: ");
+  glUniform1f(glGetUniformLocation(color_components_prog,"cmin"),s.I.minCoeff());
+  glUniform1f(glGetUniformLocation(color_components_prog,"cmax"),s.I.maxCoeff());
+  //glUniform1f(glGetUniformLocation(color_components_prog,"cc_selected"),cc_selected);
+  glUniform1f(glGetUniformLocation(color_components_prog,"cc_hover"),cc_hover);
+  glActiveTexture(GL_TEXTURE0);
+  glBindTexture(GL_TEXTURE_1D, tex_id);
+  glUniform1i(glGetUniformLocation(color_components_prog,"color_map"),0);
+  glActiveTexture(GL_TEXTURE1);
+  glBindTexture(GL_TEXTURE_1D, s.mask_id);
+  glUniform1i(glGetUniformLocation(color_components_prog,"selected_mask"),1);
+
+    igl::report_gl_error("unif: ");
+  if(fill_visible)
+  {
+    glEnable(GL_POLYGON_OFFSET_FILL); // Avoid Stitching!
+    glPolygonOffset(1.0, 0);
+  }
+  draw_mesh(V,F,N,s.I,scalar_loc);
+  glPopAttrib();
+  glUseProgram(0);
+
+  // Draw a nice floor
+  glPushMatrix();
+  const double floor_offset =
+    -2./bbd*(V.col(1).maxCoeff()-Vmid(1));
+  glTranslated(0,floor_offset,0);
+  const float GREY[4] = {0.5,0.5,0.6,1.0};
+  const float DARK_GREY[4] = {0.2,0.2,0.3,1.0};
+  draw_floor(GREY,DARK_GREY);
+  glPopMatrix();
+
+  pop_scene();
+
+  TwDraw();
+  glutSwapBuffers();
+  if(is_animating)
+  {
+    glutPostRedisplay();
+  }
+}
+
+void mouse_wheel(int wheel, int direction, int mouse_x, int mouse_y)
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  GLint viewport[4];
+  glGetIntegerv(GL_VIEWPORT,viewport);
+  if(wheel == 0 && TwMouseMotion(mouse_x, viewport[3] - mouse_y))
+  {
+    static double mouse_scroll_y = 0;
+    const double delta_y = 0.125*direction;
+    mouse_scroll_y += delta_y;
+    TwMouseWheel(mouse_scroll_y);
+    return;
+  }
+
+  auto & camera = s.camera;
+  switch(center_type)
+  {
+    case CENTER_TYPE_ORBIT:
+      if(wheel==0)
+      {
+        // factor of zoom change
+        double s = (1.-0.01*direction);
+        //// FOV zoom: just widen angle. This is hardly ever appropriate.
+        //camera.m_angle *= s;
+        //camera.m_angle = min(max(camera.m_angle,1),89);
+        camera.push_away(s);
+      }else
+      {
+        // Dolly zoom:
+        camera.dolly_zoom((double)direction*1.0);
+      }
+      break;
+    default:
+    case CENTER_TYPE_FPS:
+      // Move `eye` and `at`
+      camera.dolly((wheel==0?Vector3d(0,0,1):Vector3d(-1,0,0))*0.1*direction);
+      break;
+  }
+  glutPostRedisplay();
+}
+
+bool pick(const int x, const int y, int & cc_selected)
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+  static GLuint scalar_loc = 1;
+  static GLuint tex_id = 0;
+  static const size_t max_ids = s.I.maxCoeff()+1;
+  static GLuint color_shader_prog = color_shader(max_ids,scalar_loc,tex_id);
+  const int pick_s = 0;
+  const int pick_w = pick_s;
+  GLint old_vp[4];
+  glGetIntegerv(GL_VIEWPORT,old_vp);
+  const double pick_ratio = double(pick_w)/double(old_vp[2]);
+  // ceil, cause might otherwise round down to 0
+  const int pick_h = ceil(double(old_vp[3])*pick_ratio);
+  glViewport(
+    x-pick_w,
+    old_vp[3]-y-pick_h,2*pick_w+1,2*pick_h+1);
+  glMatrixMode(GL_PROJECTION);
+  Matrix4d proj;
+  glGetDoublev(GL_PROJECTION_MATRIX,proj.data());
+  glPushMatrix();
+  glLoadIdentity();
+  gluPickMatrix(
+    x,
+    old_vp[3]-y, 
+    pick_w*2+1,
+    pick_h*2+1,
+    old_vp);
+  glMultMatrixd(proj.data());
+  glMatrixMode(GL_MODELVIEW);
+  // Activate color shader
+  glUseProgram(color_shader_prog);
+  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT,pick_fbo);
+  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT,pick_dfbo);
+  // Clear screen
+  glClearColor(0,0,0,0);
+  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+  glPushAttrib(GL_ALL_ATTRIB_BITS);
+  glEnable(GL_TEXTURE_1D);
+  glBindTexture(GL_TEXTURE_1D, tex_id);
+  glUniform1f(glGetUniformLocation(color_shader_prog,"cmin"),s.I.minCoeff());
+  glUniform1f(glGetUniformLocation(color_shader_prog,"cmax"),s.I.maxCoeff());
+  draw_mesh(V,F,N,s.I,scalar_loc);
+  glPopAttrib();
+
+  glMatrixMode(GL_PROJECTION);
+  glPopMatrix();
+  glMatrixMode(GL_MODELVIEW);
+  glViewport(old_vp[0],old_vp[1],old_vp[2],old_vp[3]);
+
+  Matrix<GLubyte,1,4> pixel;
+  glReadPixels(x,old_vp[3]-y,1,1,GL_RGBA,GL_UNSIGNED_BYTE,pixel.data());
+
+  glUseProgram(0);
+  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT,0);
+  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT,0);
+
+  if(pixel(3) == 0)
+  {
+    cc_selected = -1;
+    return false;
+  }
+  cc_selected = pixel(0)*256*256+pixel(1)*256+pixel(2);
+  return true;
+}
+
+void regenerate_mask()
+{
+  if(glIsTexture(s.mask_id))
+  {
+    glDeleteTextures(1,&s.mask_id);
+  }
+  s.mask_id = generate_1d_texture(s.selected);
+}
+
+void mouse(int glutButton, int glutState, int mouse_x, int mouse_y)
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  bool tw_using = TwEventMouseButtonGLUT(glutButton,glutState,mouse_x,mouse_y);
+  int mod = glutGetModifiers();
+  switch(glutButton)
+  {
+    case GLUT_RIGHT_BUTTON:
+    {
+      switch(glutState)
+      {
+        case 1:
+          // up
+          glutSetCursor(GLUT_CURSOR_INHERIT);
+          is_rotating = false;
+          break;
+        case 0:
+          glutSetCursor(GLUT_CURSOR_CYCLE);
+          // collect information for trackball
+          is_rotating = true;
+          down_camera = s.camera;
+          down_x = mouse_x;
+          down_y = mouse_y;
+        break;
+      }
+      break;
+    }
+    case GLUT_LEFT_BUTTON:
+    {
+      switch(glutState)
+      {
+        case 1:
+          // up
+          glutSetCursor(GLUT_CURSOR_INHERIT);
+          is_rotating = false;
+          break;
+        case 0:
+          if(!tw_using)
+          {
+            push_scene();
+            int cc_selected=-1;
+            if(pick(mouse_x,mouse_y,cc_selected))
+            {
+              push_undo();
+              if(!(mod & GLUT_ACTIVE_SHIFT))
+              {
+                s.selected.setConstant(0);
+              }
+              s.selected(cc_selected,0) = 255;
+              regenerate_mask();
+            }else
+            {
+              glutSetCursor(GLUT_CURSOR_CYCLE);
+              // collect information for trackball
+              is_rotating = true;
+              down_camera = s.camera;
+              down_x = mouse_x;
+              down_y = mouse_y;
+            }
+            pop_scene();
+          }
+        break;
+      }
+      break;
+    }
+    // Scroll down
+    case GLUT_WHEEL_DOWN:
+    {
+      mouse_wheel(0,-1,mouse_x,mouse_y);
+      break;
+    }
+    // Scroll up
+    case GLUT_WHEEL_UP:
+    {
+      mouse_wheel(0,1,mouse_x,mouse_y);
+      break;
+    }
+    // Scroll left
+    case GLUT_WHEEL_LEFT:
+    {
+      mouse_wheel(1,-1,mouse_x,mouse_y);
+      break;
+    }
+    // Scroll right
+    case GLUT_WHEEL_RIGHT:
+    {
+      mouse_wheel(1,1,mouse_x,mouse_y);
+      break;
+    }
+  }
+  glutPostRedisplay();
+}
+
+void mouse_move(int mouse_x, int mouse_y)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  bool tw_using = TwMouseMotion(mouse_x,mouse_y);
+  push_scene();
+  pick(mouse_x,mouse_y,cc_hover);
+  pop_scene();
+  glutPostRedisplay();
+}
+
+void mouse_drag(int mouse_x, int mouse_y)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  bool tw_using = TwMouseMotion(mouse_x,mouse_y);
+
+  if(is_rotating)
+  {
+    glutSetCursor(GLUT_CURSOR_CYCLE);
+    Quaterniond q;
+    auto & camera = s.camera;
+    switch(rotation_type)
+    {
+      case ROTATION_TYPE_IGL_TRACKBALL:
+      {
+        // Rotate according to trackball
+        igl::trackball<double>(
+          width,
+          height,
+          2.0,
+          down_camera.m_rotation_conj.coeffs().data(),
+          down_x,
+          down_y,
+          mouse_x,
+          mouse_y,
+          q.coeffs().data());
+          break;
+      }
+      case ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP:
+      {
+        // Rotate according to two axis valuator with fixed up vector
+        two_axis_valuator_fixed_up(
+          width, height,
+          2.0,
+          down_camera.m_rotation_conj,
+          down_x, down_y, mouse_x, mouse_y,
+          q);
+        break;
+      }
+      default:
+        break;
+    }
+    switch(center_type)
+    {
+      default:
+      case CENTER_TYPE_ORBIT:
+        camera.orbit(q.conjugate());
+        break;
+      case CENTER_TYPE_FPS:
+        camera.turn_eye(q.conjugate());
+        break;
+    }
+  }
+  glutPostRedisplay();
+}
+
+void init_relative()
+{
+  using namespace Eigen;
+  using namespace igl;
+  per_face_normals(V,F,N);
+  Vmax = V.colwise().maxCoeff();
+  Vmin = V.colwise().minCoeff();
+  Vmid = 0.5*(Vmax + Vmin);
+  bbd = (Vmax-Vmin).norm();
+}
+
+void init_components()
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+  components(F,CC);
+  s.I = CC.cast<float>();
+  s.selected = Matrix<GLubyte,Dynamic,Dynamic>::Zero(s.I.maxCoeff()+1,3);
+  cout<<"s.selected: "<<s.selected.rows()<<endl;
+  regenerate_mask();
+}
+
+void undo()
+{
+  using namespace std;
+  if(!undo_stack.empty())
+  {
+    redo_stack.push(s);
+    s = undo_stack.top();
+    undo_stack.pop();
+  }
+  regenerate_mask();
+}
+
+void redo()
+{
+  using namespace std;
+  if(!redo_stack.empty())
+  {
+    undo_stack.push(s);
+    s = redo_stack.top();
+    redo_stack.pop();
+  }
+  regenerate_mask();
+}
+
+bool save(const std::string & out_filename)
+{
+  using namespace std;
+  using namespace igl;
+  if(write_triangle_mesh(out_filename,V,F))
+  {
+    cout<<GREENGIN("Saved mesh to `"<<out_filename<<"` successfully.")<<endl;
+    return true;
+  }else
+  {
+    cout<<REDRUM("Failed to save mesh to `"<<out_filename<<"`.")<<endl;
+    return false;
+  }
+}
+
+void TW_CALL saveCB(void * /*clientData*/)
+{
+  save(out_filename);
+}
+
+void key(unsigned char key, int mouse_x, int mouse_y)
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  int mod = glutGetModifiers();
+  switch(key)
+  {
+    // ESC
+    case char(27):
+      rebar.save(REBAR_NAME);
+    // ^C
+    case char(3):
+      exit(0);
+    case 'z':
+    case 'Z':
+      if(mod & GLUT_ACTIVE_COMMAND)
+      {
+        if(mod & GLUT_ACTIVE_SHIFT)
+        {
+          redo();
+        }else
+        {
+          undo();
+        }
+      }else
+      {
+        Quaterniond q;
+        snap_to_canonical_view_quat(s.camera.m_rotation_conj,1.0,q);
+        switch(center_type)
+        {
+          default:
+          case CENTER_TYPE_ORBIT:
+            s.camera.orbit(q.conjugate());
+            break;
+          case CENTER_TYPE_FPS:
+            s.camera.turn_eye(q.conjugate());
+            break;
+        }
+      }
+      break;
+    case 'u':
+        mouse_wheel(0, 1,mouse_x,mouse_y);
+        break;
+    case 'j':
+        mouse_wheel(0,-1,mouse_x,mouse_y);
+        break;
+    default:
+      if(!TwEventKeyboardGLUT(key,mouse_x,mouse_y))
+      {
+        cout<<"Unknown key command: "<<key<<" "<<int(key)<<endl;
+      }
+  }
+
+  glutPostRedisplay();
+}
+
+int main(int argc, char * argv[])
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  string filename = "../shared/truck.obj";
+  switch(argc)
+  {
+    case 3:
+      out_filename = argv[2];
+    case 2:
+      // Read and prepare mesh
+      filename = argv[1];
+      break;
+    default:
+      cerr<<"Usage:"<<endl<<"    ./example input.obj (output.obj)"<<endl;
+      cout<<endl<<"Opening default mesh..."<<endl;
+      break;
+  }
+
+  // print key commands
+  cout<<"[Click] and [drag]  Rotate model using trackball."<<endl;
+  cout<<"[Z,z]               Snap rotation to canonical view."<<endl;
+  cout<<"[Command+Z]         Undo."<<endl;
+  cout<<"[Shift+Command+Z]   Redo."<<endl;
+  cout<<"[^C,ESC]            Exit."<<endl;
+
+  read_triangle_mesh(filename,V,F);
+
+
+  // Init glut
+  glutInit(&argc,argv);
+  if( !TwInit(TW_OPENGL, NULL) )
+  {
+    // A fatal error occured
+    fprintf(stderr, "AntTweakBar initialization failed: %s\n", TwGetLastError());
+    return 1;
+  }
+  // Create a tweak bar
+  rebar.TwNewBar("bar");
+  TwDefine("bar label='Components' size='200 550' text=light alpha='200' color='68 68 68'");
+  rebar.TwAddVarRW("camera_rotation", TW_TYPE_QUAT4D,
+    s.camera.m_rotation_conj.coeffs().data(), "open readonly=true");
+  TwType RotationTypeTW = ReTwDefineEnumFromString("RotationType",
+    "igl_trackball,two-axis-valuator-fixed-up");
+  rebar.TwAddVarCB( "rotation_type", RotationTypeTW,
+    set_rotation_type,get_rotation_type,NULL,"keyIncr=] keyDecr=[");
+  TwType CenterTypeTW = ReTwDefineEnumFromString("CenterType","orbit,fps");
+  rebar.TwAddVarRW("center_type", CenterTypeTW,&center_type,
+    "keyIncr={ keyDecr=}");
+
+  rebar.TwAddVarRW("wireframe_visible",TW_TYPE_BOOLCPP,&wireframe_visible,"key=l");
+  rebar.TwAddVarRW("fill_visible",TW_TYPE_BOOLCPP,&fill_visible,"key=f");
+  if(out_filename != "")
+  {
+    rebar.TwAddButton("save",
+      saveCB,NULL,
+      C_STR("label='Save to `"<<out_filename<<"`' "<<
+      "key=s"));
+  }
+  rebar.load(REBAR_NAME);
+
+
+  animation_from_quat = Quaterniond(1,0,0,0);
+  s.camera.m_rotation_conj = animation_from_quat;
+  animation_start_time = get_seconds();
+
+  glutInitDisplayString( "rgba depth double samples>=8");
+  glutInitWindowSize(glutGet(GLUT_SCREEN_WIDTH)/2.0,glutGet(GLUT_SCREEN_HEIGHT)/2.0);
+  glutCreateWindow("components");
+  glutDisplayFunc(display);
+  glutReshapeFunc(reshape);
+  glutKeyboardFunc(key);
+  glutMouseFunc(mouse);
+  glutMotionFunc(mouse_drag);
+  glutPassiveMotionFunc(mouse_move);
+
+  init_components();
+  init_relative();
+  regenerate_mask();
+
+  std::cout<<"OpenGL version: "<<glGetString(GL_VERSION)<<std::endl;
+  glutMainLoop();
+
+  return 0;
+}

examples/convertmesh/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: convertmesh
 

examples/dmat/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/eigen-gotchas/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/embree/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/example_fun/Makefile

@@ -3,7 +3,7 @@
 all: example_static example_header_only
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 CFLAGS+=-g
 INC=$(LIBIGL_INC)

examples/file_contents_as_string/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/flare-eyes/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: obj example
 

examples/get_seconds/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/glslversion/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/glut_speed_test/Makefile

@@ -3,7 +3,7 @@
 all: example
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 .PHONY: example
 

examples/harwell_boeing/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/intersections/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 LIBIGL_LIB=+liglcgal
 
 all: example

examples/is_dir/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/marching_cubes/Makefile

@@ -5,7 +5,7 @@ all: example
 .PHONY: example
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 CFLAGS+=-g
 inc=$(LIBIGL_INC) $(EIGEN3_INC)

examples/meshio/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/mode/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/multi-viewport/Makefile

@@ -6,7 +6,7 @@ all: example
 
 .PHONY:  example
 
-include ../..//build/Makefile.conf
+include ../Makefile.conf
 
 CFLAGS+=-std=c++11 -g -Wno-deprecated-declarations
 

examples/patches/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 ifdef LIBIGL_USE_STATIC_LIBRARY
 	LIBIGL_LIB+=-liglembree -liglboost
 endif

examples/patches/example.cpp

@@ -1,4 +1,5 @@
 #include <igl/readOBJ.h>
+#include <igl/readPLY.h>
 #include <igl/writeOBJ.h>
 #include <igl/writeOFF.h>
 #include <igl/readWRL.h>
@@ -750,6 +751,13 @@ int main(int argc, char * argv[])
     {
       return 1;
     }
+  }else if(ext == "ply")
+  {
+    // Convert extension to lower case
+    if(!igl::readPLY(filename,vV,vF,vN,vTC))
+    {
+      return 1;
+    }
   }else if(ext == "wrl")
   {
     // Convert extension to lower case

examples/path_tests/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/pathinfo/Makefile

@@ -3,7 +3,7 @@
 all: example
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 .PHONY: example
 

examples/randomly-sample-mesh/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: obj example
 

examples/render_to_png/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 LIBIGL_LIB+=-liglpng
 
 # YIMG dependency

examples/rotate-widget/Makefile

@@ -2,7 +2,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: obj example
 

examples/scene-rotation/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: obj example
 

examples/shadow-mapping/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: obj example
 

examples/skeleton-builder/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 LIBIGL_LIB+=-liglembree
 
 all: obj skeleton_builder

examples/skeleton-poser/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 LIBIGL_LIB+=-liglbbw -liglcgal
 
 all: obj skeleton-poser

examples/skeleton/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 LIBIGL_LIB+=-liglbbw -liglmosek
 
 all: obj example

examples/slice/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/sort/Makefile

@@ -2,7 +2,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/sortrows/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/stdin_to_temp/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/svd/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/textured-mesh/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 LIBIGL_LIB+=-liglpng
 
 all: example

examples/trackball/Makefile

@@ -3,7 +3,7 @@
 all: example
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 .PHONY: example
 

examples/transparency/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/transpose_blocks/Makefile

@@ -1,7 +1,7 @@
 .PHONY: all
 
 # Shared flags etc.
-include ../../build/Makefile.conf
+include ../Makefile.conf
 
 all: example
 

examples/upright/Makefile

@@ -6,7 +6,7 @@ all: upright
 
 .PHONY: upright
 
-include ../..//build/Makefile.conf
+include ../Makefile.conf
 
 INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC)
 LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB)

include/igl/AABB.h

@@ -0,0 +1,1142 @@
+// 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;
+      //Scalar m_max_sqr_d;
+      //int m_depth;
+      AABB():
+        m_left(NULL), m_right(NULL),
+        m_box(), m_primitive(-1)
+        //m_max_sqr_d(std::numeric_limits<double>::infinity()),
+        //m_depth(0)
+    {}
+      // 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)
+        //m_max_sqr_d(other.m_max_sqr_d),
+        //m_depth(std::max(
+        //   m_left ? m_left->m_depth + 1 : 0,
+        //   m_right ? m_right->m_depth + 1 : 0))
+        {
+        }
+      // 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);
+        //swap(first.m_max_sqr_d,second.m_max_sqr_d);
+        //swap(first.m_depth,second.m_depth);
+      }
+      // 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;
+
+      template < 
+        typename Derivedother_V,
+        typename DerivedsqrD, 
+        typename DerivedI, 
+        typename DerivedC>
+      inline void squared_distance(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const AABB<Derivedother_V,DIM> & other,
+        const Eigen::PlainObjectBase<Derivedother_V> & other_V,
+        const Eigen::MatrixXi & other_Ele, 
+        Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+        Eigen::PlainObjectBase<DerivedI> & I,
+        Eigen::PlainObjectBase<DerivedC> & C) const;
+private:
+      template < 
+        typename Derivedother_V,
+        typename DerivedsqrD, 
+        typename DerivedI, 
+        typename DerivedC>
+      inline Scalar squared_distance_helper(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const AABB<Derivedother_V,DIM> * other,
+        const Eigen::PlainObjectBase<Derivedother_V> & other_V,
+        const Eigen::MatrixXi & other_Ele, 
+        const Scalar min_sqr_d,
+        Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+        Eigen::PlainObjectBase<DerivedI> & I,
+        Eigen::PlainObjectBase<DerivedC> & C) const;
+      // 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
+      inline void barycentric_coordinates(
+        const RowVectorDIMS & p, 
+        const RowVectorDIMS & a, 
+        const RowVectorDIMS & b, 
+        const RowVectorDIMS & c,
+        Eigen::Matrix<Scalar,1,SS> & bary);
+public:
+      EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+    };
+}
+
+// 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 <iomanip>
+#include <limits>
+#include <list>
+
+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);
+      //m_depth = std::max( m_left->m_depth, m_right->m_depth)+1;
+    }
+  }else
+  {
+    VectorXi allI = colon<int>(0,Ele.rows()-1);
+    MatrixXDIMS BC;
+    if(Ele.cols() == 1)
+    {
+      // points
+      BC = V;
+    }else
+    {
+      // Simplices
+      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);
+            }
+          }
+        }
+        //m_depth = 0;
+        if(LI.rows()>0)
+        {
+          m_left = new AABB();
+          m_left->init(V,Ele,SI,LI);
+          //m_depth = std::max(m_depth, m_left->m_depth+1);
+        }
+        if(RI.rows()>0)
+        {
+          m_right = new AABB();
+          m_right->init(V,Ele,SI,RI);
+          //m_depth = std::max(m_depth, m_right->m_depth+1);
+        }
+      }
+  }
+}
+
+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 = c;
+      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 = c;
+      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>
+template < 
+  typename Derivedother_V,
+  typename DerivedsqrD, 
+  typename DerivedI, 
+  typename DerivedC>
+inline void igl::AABB<DerivedV,DIM>::squared_distance(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const AABB<Derivedother_V,DIM> & other,
+  const Eigen::PlainObjectBase<Derivedother_V> & other_V,
+  const Eigen::MatrixXi & other_Ele, 
+  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C) const
+{
+  assert(other_Ele.cols() == 1 && 
+    "Only implemented for other as list of points");
+  assert(other_V.cols() == V.cols() && "other must match this dimension");
+  sqrD.setConstant(other_Ele.rows(),1,std::numeric_limits<double>::infinity());
+  I.resize(other_Ele.rows(),1);
+  C.resize(other_Ele.rows(),other_V.cols());
+  // All points in other_V currently think they need to check against root of
+  // this. The point of using another AABB is to quickly prune chunks of
+  // other_V so that most points just check some subtree of this.
+
+  // This holds a conservative estimate of max(sqr_D) where sqr_D is the
+  // current best minimum squared distance for all points in this subtree
+  double min_sqr_d = std::numeric_limits<double>::infinity();
+  squared_distance_helper(
+    V,Ele,&other,other_V,other_Ele,min_sqr_d,sqrD,I,C);
+}
+
+template <typename DerivedV, int DIM>
+template < 
+  typename Derivedother_V,
+  typename DerivedsqrD, 
+  typename DerivedI, 
+  typename DerivedC>
+inline typename igl::AABB<DerivedV,DIM>::Scalar igl::AABB<DerivedV,DIM>::squared_distance_helper(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const AABB<Derivedother_V,DIM> * other,
+  const Eigen::PlainObjectBase<Derivedother_V> & other_V,
+  const Eigen::MatrixXi & other_Ele, 
+  const Scalar min_sqr_d,
+  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C) const
+{
+  using namespace std;
+  using namespace Eigen;
+
+  // This implementation is a bit disappointing. There's no major speed up. Any
+  // performance gains seem to come from accidental cache coherency and
+  // diminish for larger "other" (the opposite of what was intended).
+
+  // Base case
+  if(other->is_leaf() && this->is_leaf())
+  {
+    Scalar sqr_d = sqrD(other->m_primitive);
+    int i = I(other->m_primitive);
+    RowVectorDIMS c = C.row(      other->m_primitive);
+    RowVectorDIMS p = other_V.row(other->m_primitive);
+    leaf_squared_distance(V,Ele,p,sqr_d,i,c);
+    sqrD( other->m_primitive) = sqr_d;
+    I(    other->m_primitive) = i;
+    C.row(other->m_primitive) = c;
+    //cout<<"leaf: "<<sqr_d<<endl;
+    //other->m_max_sqr_d = sqr_d;
+    return sqr_d;
+  }
+
+  if(other->is_leaf())
+  {
+    Scalar sqr_d = sqrD(other->m_primitive);
+    int i = I(other->m_primitive);
+    RowVectorDIMS c = C.row(      other->m_primitive);
+    RowVectorDIMS p = other_V.row(other->m_primitive);
+    sqr_d = squared_distance(V,Ele,p,sqr_d,i,c);
+    sqrD( other->m_primitive) = sqr_d;
+    I(    other->m_primitive) = i;
+    C.row(other->m_primitive) = c;
+    //other->m_max_sqr_d = sqr_d;
+    return sqr_d;
+  }
+
+  //// Exact minimum squared distance between arbitary primitives inside this and
+  //// othre's bounding boxes
+  //const auto & min_squared_distance = [&](
+  //  const AABB<DerivedV,DIM> * A,
+  //  const AABB<Derivedother_V,DIM> * B)->Scalar
+  //{
+  //  return A->m_box.squaredExteriorDistance(B->m_box);
+  //};
+
+  if(this->is_leaf())
+  {
+    //if(min_squared_distance(this,other) < other->m_max_sqr_d)
+    if(true)
+    {
+      this->squared_distance_helper(
+        V,Ele,other->m_left,other_V,other_Ele,0,sqrD,I,C);
+      this->squared_distance_helper(
+        V,Ele,other->m_right,other_V,other_Ele,0,sqrD,I,C);
+    }else
+    {
+      // This is never reached...
+    }
+    //// we know other is not a leaf
+    //other->m_max_sqr_d = std::max(other->m_left->m_max_sqr_d,other->m_right->m_max_sqr_d);
+    return 0;
+  }
+
+  // FORCE DOWN TO OTHER LEAF EVAL
+  //if(min_squared_distance(this,other) < other->m_max_sqr_d)
+  if(true)
+  {
+    if(true)
+    {
+      this->squared_distance_helper(
+        V,Ele,other->m_left,other_V,other_Ele,0,sqrD,I,C);
+      this->squared_distance_helper(
+        V,Ele,other->m_right,other_V,other_Ele,0,sqrD,I,C);
+    }else // this direction never seems to be faster
+    {
+      this->m_left->squared_distance_helper(
+        V,Ele,other,other_V,other_Ele,0,sqrD,I,C);
+      this->m_right->squared_distance_helper(
+        V,Ele,other,other_V,other_Ele,0,sqrD,I,C);
+    }
+  }else
+  {
+    // this is never reached ... :-(
+  }
+  //// we know other is not a leaf
+  //other->m_max_sqr_d = std::max(other->m_left->m_max_sqr_d,other->m_right->m_max_sqr_d);
+
+  return 0;
+#if false
+
+  // _Very_ conservative approximation of maximum squared distance between
+  // primitives inside this and other's bounding boxes
+  const auto & max_squared_distance = [](
+    const AABB<DerivedV,DIM> * A,
+    const AABB<Derivedother_V,DIM> * B)->Scalar
+  {
+    AlignedBox<Scalar,DIM> combo = A->m_box;
+    combo.extend(B->m_box);
+    return combo.diagonal().squaredNorm();
+  };
+
+  //// other base-case
+  //if(other->is_leaf())
+  //{
+  //  double sqr_d = sqrD(other->m_primitive);
+  //  int i = I(other->m_primitive);
+  //  RowVectorDIMS c = C.row(m_primitive);
+  //  RowVectorDIMS p = other_V.row(m_primitive);
+  //  leaf_squared_distance(V,Ele,p,sqr_d,i,c);
+  //  sqrD(other->m_primitive) = sqr_d;
+  //  I(other->m_primitive) = i;
+  //  C.row(m_primitive) = c;
+  //  return;
+  //}
+  std::vector<const AABB<DerivedV,DIM> * > this_list;
+  if(this->is_leaf())
+  {
+    this_list.push_back(this);
+  }else
+  {
+    assert(this->m_left);
+    this_list.push_back(this->m_left);
+    assert(this->m_right);
+    this_list.push_back(this->m_right);
+  }
+  std::vector<AABB<Derivedother_V,DIM> *> other_list;
+  if(other->is_leaf())
+  {
+    other_list.push_back(other);
+  }else
+  {
+    assert(other->m_left);
+    other_list.push_back(other->m_left);
+    assert(other->m_right);
+    other_list.push_back(other->m_right);
+  }
+
+  //const std::function<Scalar(
+  //  const AABB<Derivedother_V,DIM> * other)
+  //    > max_sqr_d = [&sqrD,&max_sqr_d](const AABB<Derivedother_V,DIM> * other)->Scalar
+  //  {
+  //    if(other->is_leaf())
+  //    {
+  //      return sqrD(other->m_primitive);
+  //    }else
+  //    {
+  //      return std::max(max_sqr_d(other->m_left),max_sqr_d(other->m_right));
+  //    }
+  //  };
+
+  //// Potentially recurse on all pairs, if minimum distance is less than running
+  //// bound
+  //Eigen::Matrix<Scalar,Eigen::Dynamic,1> other_max_sqr_d =
+  //  Eigen::Matrix<Scalar,Eigen::Dynamic,1>::Constant(other_list.size(),1,min_sqr_d);
+  for(size_t child = 0;child<other_list.size();child++)
+  {
+    auto other_tree = other_list[child];
+
+    Eigen::Matrix<Scalar,Eigen::Dynamic,1> this_max_sqr_d(this_list.size(),1);
+    for(size_t t = 0;t<this_list.size();t++)
+    {
+      const auto this_tree = this_list[t];
+      this_max_sqr_d(t) = max_squared_distance(this_tree,other_tree);
+    }
+    if(this_list.size() ==2 &&
+      ( this_max_sqr_d(0) > this_max_sqr_d(1))
+      )
+    {
+      std::swap(this_list[0],this_list[1]);
+      //std::swap(this_max_sqr_d(0),this_max_sqr_d(1));
+    }
+    const Scalar sqr_d = this_max_sqr_d.minCoeff();
+
+
+    for(size_t t = 0;t<this_list.size();t++)
+    {
+      const auto this_tree = this_list[t];
+
+      //const auto mm = max_sqr_d(other_tree);
+      //const Scalar mc = other_max_sqr_d(child);
+      //assert(mc == mm);
+      // Only look left/right in this_list if can possible decrease somebody's
+      // distance in this_tree.
+      const Scalar min_this_other = min_squared_distance(this_tree,other_tree); 
+      if(
+          min_this_other < sqr_d && 
+          min_this_other < other_tree->m_max_sqr_d)
+      {
+        //cout<<"before: "<<other_max_sqr_d(child)<<endl;
+        //other_max_sqr_d(child) = std::min(
+        //  other_max_sqr_d(child),
+        //  this_tree->squared_distance_helper(
+        //    V,Ele,other_tree,other_V,other_Ele,other_max_sqr_d(child),sqrD,I,C));
+        //cout<<"after: "<<other_max_sqr_d(child)<<endl;
+          this_tree->squared_distance_helper(
+            V,Ele,other_tree,other_V,other_Ele,0,sqrD,I,C);
+      }
+    }
+  }
+  //const Scalar ret = other_max_sqr_d.maxCoeff();
+  //const auto mm = max_sqr_d(other);
+  //assert(mm == ret);
+  //cout<<"non-leaf: "<<ret<<endl;
+  //return ret;
+  if(!other->is_leaf())
+  {
+    other->m_max_sqr_d = std::max(other->m_left->m_max_sqr_d,other->m_right->m_max_sqr_d);
+  }
+  return 0;
+#endif
+}
+
+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)
+  {
+    cout<<"pp"<<endl;
+    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
+  //std::cout<<matlab_format(c_candidate,"c_candidate")<<std::endl;
+  const Scalar pc_norm = (p-c_candidate).squaredNorm();
+  const Scalar diff = fabs(sqr_d_candidate - pc_norm);
+  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/InElementAABB.h

@@ -1,404 +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;
-  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;
-  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 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/ReAntTweakBar.h

@@ -276,9 +276,10 @@ namespace igl
 //TW_API int      TW_CALL TwRemoveVar(TwBar *bar, const char *name);
 //TW_API int      TW_CALL TwRemoveAllVars(TwBar *bar);
 
-#ifndef IGL_STATIC_LIBRARY
+// Until AntTweakBar dependency folder exists, this is header-only
+//#ifndef IGL_STATIC_LIBRARY
 #  include "ReAntTweakBar.cpp"
-#endif
+//#endif
 
 #endif
 #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);

include/igl/WindingNumberMethod.h

@@ -1,3 +1,10 @@
+// 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_WINDINGNUMBERMETHOD_H
 #define IGL_WINDINGNUMBERMETHOD_H
 namespace igl

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,

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/bone_parents.cpp

@@ -1,3 +1,10 @@
+// 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/.
 #include "bone_parents.h"
 
 template <typename DerivedBE, typename DerivedP>

include/igl/boolean/MeshBooleanType.h

@@ -1,3 +1,10 @@
+// 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 MESH_BOOLEAN_TYPE_H
 #define MESH_BOOLEAN_TYPE_H
 

include/igl/boolean/from_cork_mesh.cpp

@@ -1,3 +1,10 @@
+// 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_NO_CORK
 #include "from_cork_mesh.h"
 

include/igl/boolean/from_cork_mesh.h

@@ -1,3 +1,10 @@
+// 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_FROM_CORK_MESH_H
 #define IGL_FROM_CORK_MESH_H
 #ifndef IGL_NO_CORK

include/igl/boolean/mesh_boolean.cpp

@@ -1,11 +1,20 @@
+// 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/.
 #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 +39,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 +66,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 +91,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 +102,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 +129,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 +143,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 +163,7 @@ IGL_INLINE void igl::mesh_boolean(
   cout<<"resolve..."<<endl;
 #endif
   MatrixX3S CV;
+  MatrixX3ES EV;
   MatrixX3I CF;
   VectorXJ CJ;
   if(resolve_fun)
@@ -158,7 +171,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 +201,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 +268,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 +311,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/boolean/mesh_boolean.h

@@ -1,3 +1,10 @@
+// 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 MESH_BOOLEAN_H
 #define MESH_BOOLEAN_H
 

include/igl/boolean/mesh_boolean_cork.cpp

@@ -1,3 +1,10 @@
+// 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_NO_CORK
 #include "mesh_boolean_cork.h"
 #include "to_cork_mesh.h"

include/igl/boolean/mesh_boolean_cork.h

@@ -1,3 +1,10 @@
+// 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 MESH_BOOLEAN_CORK_H
 #define MESH_BOOLEAN_CORK_H
 #ifndef IGL_NO_CORK

include/igl/boolean/to_cork_mesh.cpp

@@ -1,3 +1,10 @@
+// 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_NO_CORK
 #include "to_cork_mesh.h"
 template <

include/igl/boolean/to_cork_mesh.h

@@ -1,3 +1,10 @@
+// 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_TO_CORK_MESH_H
 #define IGL_TO_CORK_MESH_H
 #ifndef IGL_NO_CORK

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

@@ -38,4 +38,8 @@ IGL_INLINE void igl::mesh_to_cgal_triangle_list(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+template void igl::mesh_to_cgal_triangle_list<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > >&);
+template void igl::mesh_to_cgal_triangle_list<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
+template void igl::mesh_to_cgal_triangle_list<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epeck>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<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<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
 #endif

include/igl/cgal/order_facets_around_edges.cpp

@@ -0,0 +1,254 @@
+// 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/.
+#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,86 @@
+// 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 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,30 @@
+// 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/.
 #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 <type_traits>
+#include <CGAL/number_utils.h>
 //#define IGL_OUTER_HULL_DEBUG
 
 template <
@@ -30,6 +42,9 @@ 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;
   typedef typename DerivedF::Index Index;
@@ -41,16 +56,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;
@@ -61,118 +77,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));
 #ifdef IGL_OUTER_HULL_DEBUG
-      if(di_I(fei,0) != di_I(fei,0) )
-      {
-        cout<<"NaN from face: "<<(f+1)<<endl;
-        cout<<"  n: "<<n<<endl;
-        cout<<"  eVp: "<<eVp<<endl;
-        cout<<"  eV: "<<eV<<endl;
-        cout<<"  eVp x n . eV: "<<(eVp.cross(n).dot(eV))<<endl;
-        cout<<"  eVp . n: "<<(eVp.dot(n))<<endl;
+  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] << ", ";
       }
+      std::cout << ")" << std::endl;
+  }
 #endif
-      assert(di_I(fei,0) == di_I(fei,0) && "NaN Alert!");
-      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);
-        }
+  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;
       }
-    }
-    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)];
-    }
-
   }
 
   vector<vector<vector<Index > > > TT,_1;
@@ -197,6 +126,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);
@@ -233,6 +163,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;
@@ -242,6 +175,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;
@@ -255,6 +190,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)
       {
@@ -267,9 +208,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();
@@ -282,7 +238,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++)
@@ -290,25 +246,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;
@@ -391,6 +356,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;
@@ -410,7 +377,7 @@ IGL_INLINE void igl::outer_hull(
         }
       }
     }
-    
+
     {
       vG[id].resize(FHcount,3);
       vJ[id].resize(FHcount,1);
@@ -438,18 +405,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;
@@ -467,8 +434,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 )
     {
@@ -479,34 +446,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;
@@ -519,7 +483,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
@@ -567,8 +531,14 @@ IGL_INLINE void igl::outer_hull(
   return outer_hull(V,F,N,G,J,flip);
 }
 
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+#undef IGL_STATIC_LIBRARY
+#include <igl/barycenter.cpp>
+#include <igl/outer_facet.cpp>
+#include <igl/cgal/order_facets_around_edges.cpp>
+#define IGL_STATIC_LIBRARY
 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> >&);
 template void igl::outer_hull<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::Matrix<bool, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);

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

@@ -1,6 +1,13 @@
+// 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_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,20 @@
+// 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/.
 #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,13 @@
-#ifndef PEEL_OUTER_HULL_LAYERS_H
-#define PEEL_OUTER_HULL_LAYERS_H
-#include "igl_inline.h"
+// 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_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/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