embree extra and bone visible implementation (for bone heat computation)

Former-commit-id: 64f552164141706d30ff1b97343ffac3217375fb

Makefile

@@ -25,17 +25,20 @@ ifeq ($(IGL_WITH_TETGEN),1)
 	EXTRA_DIRS+=include/igl/tetgen
 endif
 ifeq ($(IGL_WITH_MATLAB),1)
-	# append matlab extra dir
 	EXTRA_DIRS+=include/igl/matlab
 endif
 ifeq ($(IGL_WITH_MOSEK),1)
-	# append mosek extra dir
 	EXTRA_DIRS+=include/igl/mosek
 endif
 ifeq ($(IGL_WITH_PNG),1)
-	# append mosek extra dir
 	EXTRA_DIRS+=include/igl/png
 endif
+ifeq ($(IGL_WITH_XML),1)
+	EXTRA_DIRS+=include/igl/xml
+endif
+ifeq ($(IGL_WITH_EMBREE),1)
+	EXTRA_DIRS+=include/igl/embree
+endif
 
 .PHONY: examples
 .PHONY: extras

Makefile.conf

@@ -23,11 +23,13 @@ endif
 ifeq ($(IGL_USERNAME),ajx)
 	MOSEKPLATFORM=osx64x86
 	IGL_WITH_TETGEN=1
+	IGL_WITH_EMBREE=1
 	IGL_WITH_MATLAB=1
 	IGL_WITH_MOSEK=1
 	IGL_WITH_PNG=1
 	# I don't use llvm
 #AFLAGS = -m64 -march="corei7-avx"
+	# msse4.2 is necessary for me to get embree to compile correctly
 	AFLAGS = -m64 -msse4.2
 	OPENMP=-fopenmp
 #SSE=-mavx

include/igl/embree/EmbreeIntersector.cpp

@@ -0,0 +1,80 @@
+#include "EmbreeIntersector.h"
+
+template <typename RowVector3>
+inline embree::Vec3f toVec3f(const RowVector3 &p) { return embree::Vec3f((float)p[0], (float)p[1], (float)p[2]); }
+
+template <
+typename PointMatrixType,
+typename FaceMatrixType,
+typename RowVector3>
+EmbreeIntersector < PointMatrixType, FaceMatrixType, RowVector3>
+::EmbreeIntersector(const PointMatrixType & V, const FaceMatrixType & F)
+{
+  static bool inited = false;
+  if(!inited)
+  {
+	//embree::TaskScheduler::start();//init();
+    inited = true;
+  }
+  
+  size_t numVertices = 0;
+  size_t numTriangles = 0;
+  
+  triangles = (embree::BuildTriangle*) embree::rtcMalloc(sizeof(embree::BuildTriangle) * F.rows());
+  vertices = (embree::BuildVertex*) embree::rtcMalloc(sizeof(embree::BuildVertex) * V.rows());
+
+  for(int i = 0; i < (int)V.rows(); ++i)
+  {
+    vertices[numVertices++] = embree::BuildVertex((float)V(i,0),(float)V(i,1),(float)V(i,2));
+  }
+  
+  for(int i = 0; i < (int)F.rows(); ++i)
+  {
+    triangles[numTriangles++] = embree::BuildTriangle((int)F(i,0),(int)F(i,1),(int)F(i,2),i);
+  }
+  
+  _accel = embree::rtcCreateAccel("default", "default", triangles, numTriangles, vertices, numVertices);
+  _intersector = _accel->queryInterface<embree::Intersector>();
+}
+
+template <
+typename PointMatrixType,
+typename FaceMatrixType,
+typename RowVector3>
+EmbreeIntersector < PointMatrixType, FaceMatrixType, RowVector3>
+::~EmbreeIntersector()
+{
+	embree::rtcFreeMemory();
+}
+
+template <
+typename PointMatrixType,
+typename FaceMatrixType,
+typename RowVector3>
+bool 
+EmbreeIntersector < PointMatrixType, FaceMatrixType, RowVector3>
+::intersectRay(const RowVector3& origin, const RowVector3& direction, embree::Hit &hit) const
+{
+	embree::Ray ray(toVec3f(origin), toVec3f(direction), 1e-4f);
+	_intersector->intersect(ray, hit);
+	return hit ; 
+}
+
+template <
+typename PointMatrixType,
+typename FaceMatrixType,
+typename RowVector3>
+bool 
+EmbreeIntersector < PointMatrixType, FaceMatrixType, RowVector3>
+::intersectSegment(const RowVector3& a, const RowVector3& ab, embree::Hit &hit) const
+{
+	embree::Ray ray(toVec3f(a), toVec3f(ab), embree::zero, embree::one);
+	_intersector->intersect(ray, hit);
+	return hit ; 
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template instanciation
+#include <Eigen/Core>
+template class EmbreeIntersector<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >;
+#endif

include/igl/embree/EmbreeIntersector.h

@@ -0,0 +1,32 @@
+#ifndef EMBREE_INTERSECTOR_H
+#define EMBREE_INTERSECTOR_H
+
+#undef interface
+#undef near
+#undef far
+
+#include "common/intersector.h"
+#include "common/accel.h"
+//#include "types.h"
+
+template <
+typename PointMatrixType,
+typename FaceMatrixType,
+typename RowVector3>
+class EmbreeIntersector
+{
+public:
+  EmbreeIntersector(const PointMatrixType & V, const FaceMatrixType & F);
+  virtual ~EmbreeIntersector();
+  
+  bool intersectRay(const RowVector3& origin, const RowVector3& direction, embree::Hit &hit) const;
+  bool intersectSegment(const RowVector3& a, const RowVector3& ab, embree::Hit &hit) const;
+  
+private:
+  embree::BuildTriangle *triangles;
+  embree::BuildVertex *vertices;
+  embree::Ref<embree::Accel> _accel;
+  embree::Ref<embree::Intersector> _intersector;
+};
+
+#endif //EMBREE_INTERSECTOR_H

include/igl/embree/Makefile

@@ -0,0 +1,44 @@
+include ../../../Makefile.conf
+
+.PHONY: all
+all: libiglembree
+debug: libiglembree
+
+include ../../../Makefile.conf
+all: CFLAGS += -O3 -DNDEBUG
+debug: CFLAGS += -g -Wall -Werror
+
+.PHONY: libiglembree
+libiglembree: obj ../../../lib/libiglembree.a
+
+CPP_FILES=$(wildcard *.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)
+
+# embree dependency
+EMBREE=../../../external/embree
+EMBREE_INC=-I$(EMBREE)/common -I$(EMBREE)/rtcore/
+#EMBREE_LIB=-L$(EMBREE)/bin -lrtcore -lsys
+INC+=$(EMBREE_INC)
+
+obj: 
+	mkdir -p obj
+
+../../../lib/libiglembree.a: $(OBJ_FILES)
+	rm -f $@
+	ar cqs $@ $(OBJ_FILES)
+
+obj/%.o: %.cpp %.h
+	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
+
+clean:
+	rm -f obj/*.o
+	rm -f ../../../lib/libiglembree.a

include/igl/embree/bone_visible.cpp

@@ -0,0 +1,111 @@
+#include "bone_visible.h"
+#include "EmbreeIntersector.h"
+#include <igl/project_to_line.h>
+#include <igl/EPS.h>
+#include <iostream>
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedSD,
+  typename Derivedflag>
+void bone_visible(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedSD> & s,
+  const Eigen::PlainObjectBase<DerivedSD> & d,
+  Eigen::PlainObjectBase<Derivedflag>  & flag)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  flag.resize(V.rows());
+  // Initialize intersector
+  const EmbreeIntersector<MatrixXd,MatrixXi,Vector3d> ei = 
+        EmbreeIntersector<MatrixXd,MatrixXi,Vector3d>(V,F);
+  const double sd_norm = (s-d).norm();
+  // loop over mesh vertices
+  // TODO: openmp?
+  for(int v = 0;v<V.rows();v++)
+  {
+    Vector3d Vv = V.row(v);
+    // Project vertex v onto line segment sd
+    //embree.intersectSegment
+    double t,sqrd;
+    Vector3d projv;
+    // degenerate bone, just snap to s
+    if(sd_norm < DOUBLE_EPS)
+    {
+      t = 0;
+      sqrd = (Vv-s).array().pow(2).sum();
+      projv = s;
+    }else
+    {
+      // project onto (infinite) line
+      project_to_line(
+        Vv(0),Vv(1),Vv(2),s(0),s(1),s(2),d(0),d(1),d(2),
+        projv(0),projv(1),projv(2),t,sqrd);
+      // handle projections past endpoints
+      if(t<0)
+      {
+        t = 0;
+        sqrd = (Vv-s).array().pow(2).sum();
+        projv = s;
+      } else if(t>1)
+      {
+        t = 1;
+        sqrd = (Vv-d).array().pow(2).sum();
+        projv = d;
+      }
+    }
+    embree::Hit hit;
+    // perhaps 1.0 should be 1.0-epsilon, or actually since we checking the
+    // incident face, perhaps 1.0 should be 1.0+eps
+    const Vector3d dir = (Vv-projv)*1.0;
+    if(ei.intersectSegment(projv,dir, hit))
+    {
+      //if(v == 1228-1)
+      //{
+      //  int fi = hit.id0;
+      //  Vector3d bc,P;
+      //  bc << 1 - hit.u - hit.v, hit.u, hit.v; // barycentric
+      //  P = V.row(F(fi,0))*bc(0) + 
+      //      V.row(F(fi,1))*bc(1) + 
+      //      V.row(F(fi,2))*bc(2);
+      //  cout<<(fi+1)<<endl;
+      //  cout<<bc.transpose()<<endl;
+      //  cout<<P.transpose()<<endl;
+      //  cout<<Vv.transpose()<<endl;
+      //}
+
+      const int fi = hit.id0; // face id
+      // Assume hit is valid, so not visible
+      flag(v) = false;
+      // loop around corners of triangle
+      for(int c = 0;c<F.cols();c++)
+      {
+        if(F(fi,c) == v)
+        {
+          // hit self, so no hits before, so vertex v is visible
+          flag(v) = true;
+          break;
+        }
+      }
+      // Hit is actually past v
+      if(!flag(v) && (hit.t*hit.t*dir.squaredNorm())>sqrd)
+      {
+        flag(v) = true;
+      }
+    }else
+    {
+      // no hit so vectex v is visible
+      flag(v) = true;
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template instanciation
+template void bone_visible<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 1, 0, 3, 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, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+template void bone_visible<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<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<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/embree/bone_visible.h

@@ -0,0 +1,29 @@
+#ifndef BONE_VISIBLE_H
+#define BONE_VISIBLE_H
+#include <Eigen/Core>
+//
+// BONE_VISIBLE  test whether vertices of mesh are "visible" to a given bone,
+// where "visible" is defined as in [Baran & Popovic 07].
+//
+// [flag] = bone_visible(V,F,s,d);
+//
+// Input:
+//    s  row vector of position of start end point of bone
+//    d  row vector of position of dest end point of bone
+//    V  #V by 3 list of vertex positions
+//    F  #F by 3 list of triangle indices
+// Output:
+//    flag  #V by 1 list of bools (true) visible, (false) obstructed
+//
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedSD,
+  typename Derivedflag>
+void bone_visible(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedSD> & s,
+  const Eigen::PlainObjectBase<DerivedSD> & d,
+  Eigen::PlainObjectBase<Derivedflag>  & flag);
+#endif

readme.txt

@@ -25,6 +25,8 @@ GLEW  Windows only
 
   = Optional (included in external/ =
   TetGen  libigltetgen extra only
+  Embree  libiglembree extra only
+  tinyxml2  libiglxml extra only
 
 = Header only =
 libigl is designed to work "out-of-the-box" as a headers only library. To
@@ -77,7 +79,7 @@ Then run this example with:
     = Extras =
     Once you've set up an IGL_USERNAME and enabled extras within Makefile.conf.
     You can build the extra libraries (into lib/ligiglpng.a, lib/libiglmatlab.a,
-    lib/libigltetgen.a and lib/libiglmosek.a) by issuing:
+    lib/libigltetgen.a, lib/libiglmosek.a, etc.) by issuing:
   
     make extras