Offscreen meshes rendering using embree (#1652)

* Added simple raytracing renderer which allows off-screen rendering of meshes.
Using embree (already included in libigl). Addresses #1458, #601
Does not depend on OpenGL.
Tutorial example renders `fertility` example with average gaussian curvature field, in 0.1 sec.

* Separated EmbreeDevice

* Changed singleton pattern, updated variable name, removed EmbreeDevice.cpp

include/igl/embree/EmbreeDevice.h

@@ -0,0 +1,86 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Vladimir Fonov <[email protected]>
+//               2013 Alec Jacobson <[email protected]>
+//               2014 Christian Schüller <[email protected]>
+//
+// 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_EMBREE_EMBREE_DEVICE_H
+#define IGL_EMBREE_EMBREE_DEVICE_H
+#include <embree3/rtcore.h>
+#include <iostream>
+
+namespace igl
+{
+  namespace embree
+  {
+     // keep track of embree device
+     struct EmbreeDevice
+     {
+        RTCDevice embree_device;
+        int       embree_device_cntr;
+
+        static EmbreeDevice & instance()
+        {
+            static EmbreeDevice s;
+            return s;
+        } // instance
+
+        EmbreeDevice(const EmbreeDevice &) = delete;
+        EmbreeDevice & operator = (const EmbreeDevice &) = delete;
+
+        static RTCDevice get_device(const char *config=nullptr)
+        {
+            return instance().get(config);
+        }
+
+        static void release_device(void)
+        {
+            instance().release();
+        }
+
+    private:
+
+        EmbreeDevice():embree_device(nullptr),embree_device_cntr(0) {}
+
+        ~EmbreeDevice() 
+        { 
+            if(embree_device)
+                rtcReleaseDevice(embree_device);
+        }
+
+        RTCDevice get(const char *config=nullptr) 
+        {
+            if(!embree_device) 
+            {
+                embree_device = rtcNewDevice (config);
+                if(rtcGetDeviceError (embree_device) != RTC_ERROR_NONE)
+                    std::cerr << "Embree: An error occurred while initializing embree core!" << std::endl;
+            #ifdef IGL_VERBOSE
+                else
+                    std::cerr << "Embree: core initialized." << std::endl;
+            #endif
+            }
+            ++embree_device_cntr;
+            return embree_device;
+        }
+
+        void release()
+        {
+            if(!--embree_device_cntr) {
+                rtcReleaseDevice (embree_device);
+                embree_device = nullptr;                
+            #ifdef IGL_VERBOSE
+                    std::cerr << "Embree: core released." << std::endl;
+            #endif
+            }
+        }
+     };
+  }
+}
+
+#endif // IGL_EMBREE_EMBREE_DEVICE_H

include/igl/embree/EmbreeIntersector.h

@@ -19,27 +19,20 @@
 #include "../Hit.h"
 #include <Eigen/Geometry>
 #include <Eigen/Core>
-#include <Eigen/Geometry>
 
 #include <embree3/rtcore.h>
 #include <embree3/rtcore_ray.h>
 #include <iostream>
 #include <vector>
 
+#include "EmbreeDevice.h"
+
 namespace igl
 {
   namespace embree
   {
     class EmbreeIntersector
     {
-    public:
-      // Initialize embree engine. This will be called on instance `init()`
-      // calls. If already inited then this function does nothing: it is harmless
-      // to call more than once.
-      static inline void global_init();
-    private:
-      // Deinitialize the embree engine.
-      static inline void global_deinit();
     public:
       typedef Eigen::Matrix<float,Eigen::Dynamic,3> PointMatrixType;
       typedef Eigen::Matrix<int,Eigen::Dynamic,3> FaceMatrixType;
@@ -176,6 +169,8 @@ namespace igl
       Triangle* triangles;
       bool initialized;
 
+      RTCDevice device;
+
       inline void createRay(
         RTCRayHit& ray,
         const Eigen::RowVector3f& origin,
@@ -189,39 +184,7 @@ namespace igl
 
 // Implementation
 #include <igl/EPS.h>
-// This unfortunately cannot be a static field of EmbreeIntersector because it
-// would depend on the template and then we might end up with initializing
-// embree twice. If only there was a way to ask embree if it's already
-// initialized...
-namespace igl
-{
-  namespace embree
-  {
-    // Keeps track of whether the **Global** Embree intersector has been
-    // initialized. This should never been done at the global scope.
-    static RTCDevice g_device = nullptr;
-  }
-}
 
-inline void igl::embree::EmbreeIntersector::global_init()
-{
-  if(!g_device)
-  {
-    g_device = rtcNewDevice (NULL);
-    if(rtcGetDeviceError (g_device) != RTC_ERROR_NONE)
-      std::cerr << "Embree: An error occurred while initializing embree core!" << std::endl;
-#ifdef IGL_VERBOSE
-    else
-      std::cerr << "Embree: core initialized." << std::endl;
-#endif
-  }
-}
-
-inline void igl::embree::EmbreeIntersector::global_deinit()
-{
-  rtcReleaseDevice (g_device);
-  g_device = nullptr;
-}
 
 inline igl::embree::EmbreeIntersector::EmbreeIntersector()
   :
@@ -229,7 +192,8 @@ inline igl::embree::EmbreeIntersector::EmbreeIntersector()
   geomID(0),
   vertices(NULL),
   triangles(NULL),
-  initialized(false)
+  initialized(false),
+  device(igl::embree::EmbreeDevice::get_device())
 {
 }
 
@@ -278,8 +242,7 @@ inline void igl::embree::EmbreeIntersector::init(
     deinit();
 
   using namespace std;
-  global_init();
-
+  
   if(V.size() == 0 || F.size() == 0)
   {
     std::cerr << "Embree: No geometry specified!";
@@ -289,14 +252,14 @@ inline void igl::embree::EmbreeIntersector::init(
   RTCBuildQuality buildQuality = isStatic ? RTC_BUILD_QUALITY_HIGH : RTC_BUILD_QUALITY_MEDIUM;
 
   // create a scene
-  scene = rtcNewScene(g_device);
+  scene = rtcNewScene(device);
   rtcSetSceneFlags(scene, RTC_SCENE_FLAG_ROBUST);
   rtcSetSceneBuildQuality(scene, buildQuality);
 
   for(int g=0;g<(int)V.size();g++)
   {
     // create triangle mesh geometry in that scene
-    RTCGeometry geom_0 = rtcNewGeometry (g_device, RTC_GEOMETRY_TYPE_TRIANGLE);
+    RTCGeometry geom_0 = rtcNewGeometry (device, RTC_GEOMETRY_TYPE_TRIANGLE);
     rtcSetGeometryBuildQuality(geom_0,buildQuality);
     rtcSetGeometryTimeStepCount(geom_0,1);
     geomID = rtcAttachGeometry(scene,geom_0);
@@ -328,7 +291,7 @@ inline void igl::embree::EmbreeIntersector::init(
 
   rtcCommitScene(scene);
 
-  if(rtcGetDeviceError (g_device) != RTC_ERROR_NONE)
+  if(rtcGetDeviceError (device) != RTC_ERROR_NONE)
       std::cerr << "Embree: An error occurred while initializing the provided geometry!" << endl;
 #ifdef IGL_VERBOSE
   else
@@ -343,15 +306,16 @@ igl::embree::EmbreeIntersector
 {
   if(initialized)
     deinit();
+  igl::embree::EmbreeDevice::release_device();
 }
 
 void igl::embree::EmbreeIntersector::deinit()
 {
-  if(g_device && scene)
+  if(device && scene)
   {
     rtcReleaseScene(scene);
 
-    if(rtcGetDeviceError (g_device) != RTC_ERROR_NONE)
+    if(rtcGetDeviceError (device) != RTC_ERROR_NONE)
     {
         std::cerr << "Embree: An error occurred while resetting!" << std::endl;
     }
@@ -386,7 +350,7 @@ inline bool igl::embree::EmbreeIntersector::intersectRay(
     ray.hit.Ng_z = -ray.hit.Ng_z;
   }
 #ifdef IGL_VERBOSE
-  if(rtcGetDeviceError (g_device) != RTC_ERROR_NONE)
+  if(rtcGetDeviceError (device) != RTC_ERROR_NONE)
       std::cerr << "Embree: An error occurred while resetting!" << std::endl;
 #endif
 

include/igl/embree/EmbreeRenderer.cpp

@@ -0,0 +1,423 @@
+#include "EmbreeRenderer.h"
+// Implementation
+
+//IGL viewing parts
+#include "../unproject.h"
+#include "../look_at.h"
+#include "../frustum.h"
+#include "../ortho.h"
+
+// color map
+#include "../jet.h"
+
+#include "../PI.h"
+
+
+IGL_INLINE void igl::embree::EmbreeRenderer::init_view()
+{
+  camera_base_zoom = 1.0f;
+  camera_zoom = 1.0f;
+
+  camera_view_angle = 45.0;
+  camera_dnear = 1.0;
+  camera_dfar = 100.0;
+  camera_base_translation << 0, 0, 0;
+  camera_translation << 0, 0, 0;
+  camera_eye << 0, 0, 5;
+  camera_center << 0, 0, 0;
+  camera_up << 0, 1, 0;
+  
+  rot_matrix = Eigen::Matrix3f::Identity();
+
+  view = Eigen::Matrix4f::Identity();
+  proj = Eigen::Matrix4f::Identity();
+  norm = Eigen::Matrix4f::Identity();
+
+  orthographic = false;
+  
+
+  uC << 1,0,0;
+}
+
+IGL_INLINE igl::embree::EmbreeRenderer::EmbreeRenderer()
+  :
+  scene(NULL),
+  geomID(0),
+  initialized(false),
+  device(igl::embree::EmbreeDevice::get_device())
+{
+  init_view();
+}
+
+IGL_INLINE igl::embree::EmbreeRenderer::EmbreeRenderer(
+  const EmbreeRenderer &)
+  :// To make -Weffc++ happy
+  scene(NULL),
+  geomID(0),
+  initialized(false)
+{
+  assert(false && "Embree: Copying EmbreeRenderer is not allowed");
+}
+
+IGL_INLINE igl::embree::EmbreeRenderer & igl::embree::EmbreeRenderer::operator=(
+  const EmbreeRenderer &)
+{
+  assert(false && "Embree: Assigning an EmbreeRenderer is not allowed");
+  return *this;
+}
+
+
+IGL_INLINE void igl::embree::EmbreeRenderer::init(
+  const PointMatrixType& V,
+  const FaceMatrixType& F,
+  bool isStatic)
+{
+  std::vector<const PointMatrixType*> Vtemp;
+  std::vector<const FaceMatrixType*> Ftemp;
+  std::vector<int> masks;
+  Vtemp.push_back(&V);
+  Ftemp.push_back(&F);
+  masks.push_back(0xFFFFFFFF);
+  init(Vtemp,Ftemp,masks,isStatic);
+}
+
+IGL_INLINE void igl::embree::EmbreeRenderer::init(
+  const std::vector<const PointMatrixType*>& V,
+  const std::vector<const FaceMatrixType*>&  F,
+  const std::vector<int>& masks,
+  bool isStatic)
+{
+  if(initialized)
+    deinit();
+
+  using namespace std;
+
+  if(V.size() == 0 || F.size() == 0)
+  {
+    std::cerr << "Embree: No geometry specified!";
+    return;
+  }
+  RTCBuildQuality buildQuality = isStatic ? RTC_BUILD_QUALITY_HIGH : RTC_BUILD_QUALITY_MEDIUM;
+
+  // create a scene
+  scene = rtcNewScene(device);
+  rtcSetSceneFlags(scene, RTC_SCENE_FLAG_ROBUST);
+  rtcSetSceneBuildQuality(scene, buildQuality);
+
+  for(int g=0;g<(int)V.size();g++)
+  {
+    // create triangle mesh geometry in that scene
+    RTCGeometry geom_0 = rtcNewGeometry (device, RTC_GEOMETRY_TYPE_TRIANGLE);
+    rtcSetGeometryBuildQuality(geom_0, buildQuality);
+    rtcSetGeometryTimeStepCount(geom_0,1);
+    geomID = rtcAttachGeometry(scene,geom_0);
+    rtcReleaseGeometry(geom_0);
+
+    // fill vertex buffer, have to be 16 byte wide( sizeof(float)*4 )
+    Eigen::Map<Eigen::Matrix<float,-1,4,Eigen::RowMajor>> vertices(
+        (float*)rtcSetNewGeometryBuffer(geom_0,RTC_BUFFER_TYPE_VERTEX,0,RTC_FORMAT_FLOAT3,4*sizeof(float),V[g]->rows()),
+        V[g]->rows(),4
+    );
+    vertices.block(0,0,V[g]->rows(),3) = V[g]->cast<float>();
+
+    // fill triangle buffer
+    Eigen::Map<Eigen::Matrix<unsigned int,-1,3,Eigen::RowMajor>> triangles(
+      (unsigned int*) rtcSetNewGeometryBuffer(geom_0,RTC_BUFFER_TYPE_INDEX,0,RTC_FORMAT_UINT3,3*sizeof(unsigned int), F[g]->rows()),
+      F[g]->rows(),3
+    );
+    triangles = F[g]->cast<unsigned int>();
+    //TODO: store vertices and triangles in array for whatever reason?
+
+    rtcSetGeometryMask(geom_0, masks[g]);
+    rtcCommitGeometry(geom_0);
+  }
+
+  rtcCommitScene(scene);
+
+  if(rtcGetDeviceError (device) != RTC_ERROR_NONE)
+      std::cerr << "Embree: An error occurred while initializing the provided geometry!" << endl;
+#ifdef IGL_VERBOSE
+  else
+    std::cerr << "Embree: geometry added." << endl;
+#endif
+
+  initialized = true;
+}
+
+IGL_INLINE igl::embree::EmbreeRenderer::~EmbreeRenderer()
+{
+  if(initialized)
+    deinit();
+  
+  igl::embree::EmbreeDevice::release_device();
+}
+
+IGL_INLINE void igl::embree::EmbreeRenderer::deinit()
+{
+  if(scene)
+  {
+    rtcReleaseScene(scene);
+
+    if(rtcGetDeviceError (device) != RTC_ERROR_NONE)
+    {
+        std::cerr << "Embree: An error occurred while resetting!" << std::endl;
+    }
+#ifdef IGL_VERBOSE
+    else
+    {
+      std::cerr << "Embree: geometry removed." << std::endl;
+    }
+#endif
+  }
+}
+
+IGL_INLINE bool igl::embree::EmbreeRenderer::intersect_ray(
+  const Eigen::RowVector3f& origin,
+  const Eigen::RowVector3f& direction,
+  Hit&  hit,
+  float tnear,
+  float tfar,
+  int mask) const
+{
+  RTCRayHit ray;
+  ray.ray.flags = 0;
+  create_ray(ray, origin,direction,tnear,tfar,mask);
+
+  // shot ray
+  {
+    RTCIntersectContext context;
+    rtcInitIntersectContext(&context);
+    rtcIntersect1(scene, &context, &ray);
+    ray.hit.Ng_x = -ray.hit.Ng_x; // EMBREE_FIXME: only correct for triangles,quads, and subdivision surfaces
+    ray.hit.Ng_y = -ray.hit.Ng_y;
+    ray.hit.Ng_z = -ray.hit.Ng_z;
+  }
+#ifdef IGL_VERBOSE
+  if(rtcGetDeviceError (device) != RTC_ERROR_NONE)
+      std::cerr << "Embree: An error occurred while resetting!" << std::endl;
+#endif
+
+  if((unsigned)ray.hit.geomID != RTC_INVALID_GEOMETRY_ID)
+  {
+    hit.id = ray.hit.primID;
+    hit.gid = ray.hit.geomID;
+    hit.u = ray.hit.u;
+    hit.v = ray.hit.v;
+    hit.t = ray.ray.tfar;
+    hit.N = Vec3f(ray.hit.Ng_x, ray.hit.Ng_y, ray.hit.Ng_z);
+    return true;
+  }
+
+  return false;
+}
+
+IGL_INLINE void 
+igl::embree::EmbreeRenderer
+::create_ray(RTCRayHit& ray, const Eigen::RowVector3f& origin, 
+ const Eigen::RowVector3f& direction, float tnear, float tfar, int mask) const
+{
+  ray.ray.org_x = origin[0];
+  ray.ray.org_y = origin[1];
+  ray.ray.org_z = origin[2];
+  ray.ray.dir_x = direction[0];
+  ray.ray.dir_y = direction[1];
+  ray.ray.dir_z = direction[2];
+  ray.ray.tnear = tnear;
+  ray.ray.tfar = tfar;
+  ray.ray.id = RTC_INVALID_GEOMETRY_ID;
+  ray.ray.mask = mask;
+  ray.ray.time = 0.0f;
+
+  ray.hit.geomID = RTC_INVALID_GEOMETRY_ID;
+  ray.hit.instID[0] = RTC_INVALID_GEOMETRY_ID;
+  ray.hit.primID = RTC_INVALID_GEOMETRY_ID;
+}
+
+
+IGL_INLINE void
+igl::embree::EmbreeRenderer
+::set_mesh(const Eigen::Matrix<double,Eigen::Dynamic,3> & MV,
+           const Eigen::Matrix<int,  Eigen::Dynamic,3>  & MF,
+           bool is_static)
+{
+  V = MV.cast<float>();
+  F = MF;
+  this->init(V,F,is_static);
+
+  auto min_point = V.colwise().minCoeff();
+  auto max_point = V.colwise().maxCoeff();
+  auto centroid  = (0.5*(min_point + max_point)).eval();
+
+  camera_base_translation.setConstant(0);
+  camera_base_translation.head(centroid.size()) = -centroid.cast<float>();
+  camera_base_zoom = 2.0 / (max_point-min_point).array().abs().maxCoeff();
+}
+
+IGL_INLINE void
+igl::embree::EmbreeRenderer
+::render_buffer(PixelMatrixType& R, PixelMatrixType&G, PixelMatrixType &B,PixelMatrixType &A)
+{
+  assert(R.rows()==G.rows());assert(R.rows()==B.rows());assert(R.rows()==A.rows());
+  assert(R.cols()==G.cols());assert(R.cols()==B.cols());assert(R.cols()==A.cols());
+
+  Eigen::Vector4f viewport(0,0,R.rows(),R.cols());
+
+  float width  = R.rows();
+  float height = R.cols();
+
+  // update view matrix
+  igl::look_at( camera_eye, camera_center, camera_up, view);
+
+  view = view
+    * (rot_matrix * Eigen::Scaling(camera_zoom * camera_base_zoom)
+    * Eigen::Translation3f(camera_translation + camera_base_translation)).matrix();
+  
+  if (orthographic)
+  {
+    float length = (camera_eye - camera_center).norm();
+    float h = tan(camera_view_angle/360.0 * igl::PI) * (length);
+    igl::ortho(-h*width/height, h*width/height, -h, h, camera_dnear, camera_dfar, proj);
+  } else {
+    float fH = tan(camera_view_angle / 360.0 * igl::PI) * camera_dnear;
+    float fW = fH * (double)width/(double)height;
+    igl::frustum(-fW, fW, -fH, fH, camera_dnear, camera_dfar, proj);
+  }
+
+  // go over all pixels in the "view"
+  for(int x=0;x<(int)width;++x) 
+  {
+    for(int y=0;y<(int)height;++y) 
+    {
+      Vec3f s,d,dir;
+      igl::embree::EmbreeRenderer::Hit hit;
+
+      // look into the screen
+      Vec3f win_s(x,y,0);
+      Vec3f win_d(x,y,1);
+      // Source, destination and direction in world
+      igl::unproject(win_s,this->view,this->proj,viewport,s);
+      igl::unproject(win_d,this->view,this->proj,viewport,d);
+      dir = d-s;
+      dir.normalize();
+
+      auto clamp=[](float x)->unsigned char {return (unsigned char)(x<0?0:x>1.0?255:x*255);};
+
+      if(this->intersect_ray(s,dir,hit))
+      {
+        if ( dir.dot(hit.N) > 0.0f)
+        {
+          // TODO: interpolate normals ?
+          hit.N.normalize();
+          // cos between ray and face normal
+          float face_proj=dir.dot(hit.N);
+
+          Eigen::RowVector3f c;
+
+          if(this->uniform_color)
+          {
+            // same color for the whole mesh
+            c=uC;
+          } else if(this->face_based) {
+            // flat color per face
+            c=this->C.row(hit.id);
+          } else { //use barycentric coordinates to interpolate colour
+              c=this->C.row(F(hit.id,1))*hit.u+
+                this->C.row(F(hit.id,2))*hit.v+
+                this->C.row(F(hit.id,0))*(1.0-hit.u-hit.v);
+          }
+
+          R(x,y) = clamp(face_proj*c(0));
+          G(x,y) = clamp(face_proj*c(1));
+          B(x,y) = clamp(face_proj*c(2));
+        }
+        // give the same alpha to all points with something behind
+        A(x,y)=255; 
+      } else {
+        R(x,y)=0;
+        G(x,y)=0;
+        B(x,y)=0;
+        A(x,y)=0;
+      }
+
+    }
+  }
+}
+
+
+IGL_INLINE void
+igl::embree::EmbreeRenderer
+::set_colors(const Eigen::MatrixXd & C)
+{
+  if(C.rows()==V.rows()) // per vertex color
+  {
+    face_based = false;
+    this->C = C.cast<float>();
+    this->uniform_color=false;
+  } else if (C.rows()==F.rows()) {
+    face_based = true;
+    this->C = C.cast<float>();
+    this->uniform_color=false;
+  } else if (C.rows()==1) {
+    face_based = true;
+    this->uC = C.cast<float>();
+    this->uniform_color=true;
+  }else {
+    // don't know what to do
+    this->uniform_color=true;
+    assert(false); //?
+  }
+}
+
+IGL_INLINE void
+igl::embree::EmbreeRenderer
+::set_data(const Eigen::VectorXd & D, igl::ColorMapType cmap)
+{
+  const double caxis_min = D.minCoeff();
+  const double caxis_max = D.maxCoeff();
+  return set_data(D,caxis_min,caxis_max,cmap);
+}
+
+
+IGL_INLINE void igl::embree::EmbreeRenderer::set_data(
+  const Eigen::VectorXd & D,
+  double caxis_min,
+  double caxis_max,
+  igl::ColorMapType cmap)
+{
+    Eigen::MatrixXd C;
+    igl::colormap(cmap,D,caxis_min,caxis_max,C);
+    set_colors(C);
+}
+
+IGL_INLINE void 
+igl::embree::EmbreeRenderer::set_rot(const Eigen::Matrix3d &r)
+{
+  this->rot_matrix = r.cast<float>();
+}
+
+IGL_INLINE void 
+igl::embree::EmbreeRenderer::set_zoom(double zoom)
+{
+  this->camera_zoom=zoom;
+}
+
+IGL_INLINE void 
+igl::embree::EmbreeRenderer::set_translation(const Eigen::Vector3d &tr)
+{
+  this->camera_translation=tr.cast<float>();
+}
+
+IGL_INLINE void 
+igl::embree::EmbreeRenderer::set_face_based(bool _f)
+{
+  this->face_based=_f;
+}
+
+IGL_INLINE void 
+igl::embree::EmbreeRenderer::set_orthographic(bool o)
+{
+  this->orthographic=o;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+#endif //IGL_STATIC_LIBRARY

include/igl/embree/EmbreeRenderer.h

@@ -0,0 +1,255 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// 
+// Copyright (C) 2020 Vladimir Fonov <[email protected]>
+//               2013 Alec Jacobson <[email protected]>
+//               2014 Christian Schüller <[email protected]>
+//
+// 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_EMBREE_EMBREE_RENDERER_H
+#define IGL_EMBREE_EMBREE_RENDERER_H
+
+#include <igl/colormap.h>
+
+#include <Eigen/Geometry>
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+#include <embree3/rtcore.h>
+#include <embree3/rtcore_ray.h>
+#include <iostream>
+#include <vector>
+
+#include "EmbreeDevice.h"
+
+
+namespace igl
+{
+  namespace embree
+  {
+    // embree-based mesh renderer
+    class EmbreeRenderer
+    {
+    public:
+      typedef Eigen::RowVector3f Vec3f;
+
+      struct Hit
+      {
+        int id;    // primitive id
+        int gid;   // geometry id
+        float u,v; // barycentric coordinates
+        float t;   // distance = direction*t to intersection
+        Vec3f N;    // element normal
+      };
+
+      // Initialize embree engine. This will be called on instance `init()`
+      // calls. If already inited then this function does nothing: it is harmless
+      // to call more than once.
+      static IGL_INLINE void global_init();
+    private:
+      // Deinitialize the embree engine.
+      static IGL_INLINE void global_deinit();
+    public:
+      typedef Eigen::Matrix<float,Eigen::Dynamic,3> PointMatrixType;
+      typedef Eigen::Matrix<float,Eigen::Dynamic,3> ColorMatrixType;
+      typedef Eigen::Matrix<int,  Eigen::Dynamic,3> FaceMatrixType;
+      
+      typedef Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> PixelMatrixType;
+
+    public:
+      IGL_INLINE EmbreeRenderer();
+    private:
+      // Copying and assignment are not allowed.
+      IGL_INLINE EmbreeRenderer(const EmbreeRenderer & that);
+      IGL_INLINE EmbreeRenderer & operator=(const EmbreeRenderer &);
+    public:
+      virtual inline ~EmbreeRenderer();
+      
+      // Specify mesh, this call reinitializes embree structures
+      // Inputs:
+      //   V  #V x dim matrix of vertex coordinates
+      //   F  #F x simplex_size  matrix of indices of simplex corners into V
+      //   is_static - optimize for static thene (HQ rendering)
+      IGL_INLINE void set_mesh(const Eigen::Matrix<double,Eigen::Dynamic,3> & V,
+                    const Eigen::Matrix<int,  Eigen::Dynamic,3>  & F,
+                    bool is_static=true);
+      // Specify per-vertex or per-face color
+      // Inputs:
+      //   C  #V x 3 matrix of vertex colors
+      //    or #F x 3 matrix of face colors
+      //    or 1 x 3 matrix of uniform color
+      IGL_INLINE void set_colors(const Eigen::MatrixXd & C);
+
+
+      // Use min(D) and max(D) to set caxis.
+      IGL_INLINE void set_data(const Eigen::VectorXd & D,
+                    igl::ColorMapType cmap = igl::COLOR_MAP_TYPE_VIRIDIS);
+  
+      // Specify per-vertex or per-face scalar field
+      //   that will be converted to color using jet color map
+      // Inputs:
+      //   caxis_min  caxis minimum bound
+      //   caxis_max  caxis maximum bound
+      //   D  #V by 1 list of scalar values
+      //   cmap colormap type
+      //   num_steps number of intervals to discretize the colormap
+      IGL_INLINE void set_data(
+        const Eigen::VectorXd & D,
+        double caxis_min,
+        double caxis_max,
+        igl::ColorMapType cmap = igl::COLOR_MAP_TYPE_VIRIDIS);
+
+      // Specify mesh rotation
+      // Inputs:
+      //   r  3 x 3 rotaton matrix
+      IGL_INLINE void set_rot(const Eigen::Matrix3d &r);
+
+      // Specify mesh magnification
+      // Inputs:
+      //   z  magnification ratio
+      IGL_INLINE void set_zoom(double z);
+
+      // Specify mesh translation
+      // Inputs:
+      //   tr  translation vector 
+      IGL_INLINE void set_translation(const Eigen::Vector3d &tr);
+
+      // Specify that color is face based
+      // Inputs:
+      //    f - face or vertex colours
+      IGL_INLINE void set_face_based(bool f);
+
+      // Use orthographic projection
+      // Inputs:
+      //    f - orthographic or perspective projection
+      IGL_INLINE void set_orthographic(bool f );
+
+      // render full buffer
+      // Outputs:
+      //   all outputs should have the same size (size of the output picture)
+      //     area outside of the visible object will have zero alpha component (transparant)
+      //   R - red channel
+      //   G - green channel
+      //   B - blue channel
+      //   A - alpha channel
+      IGL_INLINE void render_buffer(PixelMatrixType &R,
+                         PixelMatrixType &G, 
+                         PixelMatrixType &B,
+                         PixelMatrixType &A);
+
+      // Given a ray find the first hit
+      //
+      // Inputs:
+      //   origin     3d origin point of ray
+      //   direction  3d (not necessarily normalized) direction vector of ray
+      //   tnear      start of ray segment
+      //   tfar       end of ray segment
+      //   mask      a 32 bit mask to identify active geometries.
+      // Output:
+      //   hit        information about hit
+      // Returns true if and only if there was a hit
+      IGL_INLINE bool intersect_ray(
+        const Eigen::RowVector3f& origin,
+        const Eigen::RowVector3f& direction,
+        Hit& hit,
+        float tnear = 0,
+        float tfar = std::numeric_limits<float>::infinity(),
+        int mask = 0xFFFFFFFF) const;
+
+    private:
+
+      // Initialize with a given mesh.
+      //
+      // Inputs:
+      //   V  #V by 3 list of vertex positions
+      //   F  #F by 3 list of Oriented triangles
+      //   isStatic  scene is optimized for static geometry
+      // Side effects:
+      //   The first time this is ever called the embree engine is initialized.
+      IGL_INLINE void init(
+        const PointMatrixType& V,
+        const FaceMatrixType& F,
+        bool isStatic = false);
+
+      // Initialize embree with a given mesh.
+      //
+      // Inputs:
+      //   V  vector of #V by 3 list of vertex positions for each geometry
+      //   F  vector of #F by 3 list of Oriented triangles for each geometry
+      //   masks  a 32 bit mask to identify active geometries.
+      //   isStatic  scene is optimized for static geometry
+      // Side effects:
+      //   The first time this is ever called the embree engine is initialized.
+      IGL_INLINE void init(
+        const std::vector<const PointMatrixType*>& V,
+        const std::vector<const FaceMatrixType*>& F,
+        const std::vector<int>& masks,
+        bool isStatic = false);
+
+
+      // Deinitialize embree datasctructures for current mesh.  Also called on
+      // destruction: no need to call if you just want to init() once and
+      // destroy.
+      IGL_INLINE void deinit();
+      // initialize view parameters
+      IGL_INLINE void init_view();
+
+      // scene data
+      PointMatrixType V; // vertices
+      FaceMatrixType  F; // faces
+      ColorMatrixType C; // colours
+
+      Eigen::RowVector3f uC; // uniform color
+
+      bool face_based;
+      bool uniform_color;
+
+      // Camera parameters
+      float camera_base_zoom;
+      float camera_zoom;
+      
+      Eigen::Vector3f camera_base_translation;
+      Eigen::Vector3f camera_translation;
+      Eigen::Vector3f camera_eye;
+      Eigen::Vector3f camera_up;
+      Eigen::Vector3f camera_center;
+      float camera_view_angle;
+      float camera_dnear;
+      float camera_dfar;
+
+      // projection matrixes      
+      Eigen::Matrix4f view;
+      Eigen::Matrix4f proj;
+      Eigen::Matrix4f norm;
+
+      Eigen::Matrix3f rot_matrix;
+
+      bool orthographic;
+
+      // embree data
+      RTCScene scene;
+      unsigned geomID;
+      bool initialized;
+
+      RTCDevice device;
+
+      IGL_INLINE void create_ray(
+        RTCRayHit& ray,
+        const Eigen::RowVector3f& origin,
+        const Eigen::RowVector3f& direction,
+        float tnear,
+        float tfar,
+        int mask) const;
+
+    };
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "EmbreeRenderer.cpp"
+#endif
+#endif //IGL_EMBREE_EMBREE_RENDERER_H

include/igl/unproject.cpp

@@ -72,4 +72,5 @@ template Eigen::Matrix<double, 3, 1, 0, 3, 1> igl::unproject<double>(Eigen::Matr
 template void igl::unproject<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 1, 0, 4, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 1, 0, 4, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::unproject<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 1, 0, 4, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 1, 0, 4, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&);
 template void igl::unproject<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 1, 0, 4, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 1, 0, 4, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);
+template void igl::unproject<Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 1, 0, 4, 1>, Eigen::Matrix<float, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 4, 1, 0, 4, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> >&);
 #endif

tutorial/608_RayTrace/CMakeLists.txt

@@ -0,0 +1,5 @@
+get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+project(${PROJECT_NAME})
+
+add_executable(${PROJECT_NAME} main.cpp)
+target_link_libraries(${PROJECT_NAME} igl::core igl::embree igl::png tutorials)

tutorial/608_RayTrace/main.cpp

@@ -0,0 +1,70 @@
+#include <igl/gaussian_curvature.h>
+#include <igl/massmatrix.h>
+#include <igl/invert_diag.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/png/writePNG.h>
+#include <igl/PI.h>
+#include <Eigen/Geometry>
+// embree
+#include <igl/embree/EmbreeRenderer.h>
+
+#include "tutorial_shared_path.h"
+
+#include <iostream>
+
+int main(int argc, char *argv[])
+{
+  int width=640;
+  int height=480;
+
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+
+  const char *mesh_file=argc > 1 ? argv[1] : TUTORIAL_SHARED_PATH "/fertility.off";
+  const char *png_file=argc > 2 ? argv[2]: "fertility_curvature.png";
+
+  // Load mesh 
+  igl::read_triangle_mesh(mesh_file, V,F);
+  
+  Eigen::VectorXd K;
+  // Compute integral of Gaussian curvature
+  igl::gaussian_curvature(V,F,K);
+  // Compute mass matrix
+  Eigen::SparseMatrix<double> M,Minv;
+  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,M);
+  igl::invert_diag(M,Minv);
+  // Divide by area to get integral average
+  K = (Minv*K).eval();
+
+  // embree object
+  igl::embree::EmbreeRenderer er;
+  er.set_mesh(V,F,true);
+
+  //er.set_uniform_color(Eigen::RowVector3d(0.3,0.3,0.3));
+  er.set_data(K,igl::COLOR_MAP_TYPE_JET);
+
+  Eigen::Matrix3d rot_matrix;
+
+  // specify rotation
+  rot_matrix =  Eigen::AngleAxisd( 10*igl::PI/180.0, Eigen::Vector3d::UnitX())
+              * Eigen::AngleAxisd(  5*igl::PI/180.0, Eigen::Vector3d::UnitY())
+              * Eigen::AngleAxisd(  4*igl::PI/180.0, Eigen::Vector3d::UnitZ());
+  er.set_rot(rot_matrix);
+
+  er.set_zoom(1.5);
+  er.set_orthographic(false);
+     
+  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> R(width, height);
+  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> G(width, height);
+  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> B(width, height);
+  Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> A(width, height);
+
+  // render view using embree
+  er.render_buffer(R,G,B,A);
+
+  std::cout<<"Rendered scene saved to "<<png_file<<std::endl;
+
+  // save to PNG file
+  igl::png::writePNG(R,G,B,A,png_file);
+  return 0;
+}

tutorial/CMakeLists.txt

@@ -134,6 +134,9 @@ if(TUTORIALS_CHAPTER6)
     add_subdirectory("111_MatCap")
     add_subdirectory("607_ScreenCapture")
   endif()
+  if(LIBIGL_WITH_PNG AND LIBIGL_WITH_EMBREE)
+    add_subdirectory("608_RayTrace")
+  endif()
   if(LIBIGL_WITH_CGAL)
     add_subdirectory("609_Boolean")
     add_subdirectory("610_CSGTree")