ambient occlusion mex, wrapper for ambient occlusion with just mesh

Former-commit-id: 855b89c2af9603bbc1e447e941325302a36e4163

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
-0.3.3
+0.3.4

examples/ambient-occlusion-mex/README

@@ -0,0 +1,10 @@
+Small example showing how to build a mex function to call the
+igl/embree/ambient_occlusion function from matlab.
+
+To compile open matlab, `cd` to this directory and issue:
+
+    compile
+
+To run the example, add the IGL_TOOLBOX to your path and issue:
+
+    example

examples/ambient-occlusion-mex/ambient_occlusion.m

@@ -0,0 +1,14 @@
+  % AMBIENT_OCCLUSION Compute ambient occlusion per given point
+  %
+  % S = ambient_occlusion(V,F,P,N,num_samples)
+  %
+  % Inputs:
+  %    V  #V by 3 list of mesh vertex positions
+  %    F  #F by 3 list of mesh triangle facet indices into V
+  %    P  #P by 3 list of origin points
+  %    N  #P by 3 list of origin normals
+  %    num_samples  number of samples
+  % Outputs:
+  %    S  #P list of ambient occlusion values between 1 (fully occluded) and 0
+  %      (not occluded)
+  %

examples/ambient-occlusion-mex/compile.m

@@ -0,0 +1,9 @@
+mex -v -o ambient_occlusion -DMEX -largeArrayDims ...
+  -I/opt/local/include/eigen3 -I/usr/local/igl/libigl/include ...
+  -L/usr/local/igl/libigl/lib -ligl -liglmatlab -liglembree ...
+  CXXFLAGS="\$CXXFLAGS -m64 -msse4.2 -fopenmp" ...
+  -I/usr/local/igl/libigl/external/embree/rtcore ...
+  -I/usr/local/igl/libigl/external/embree/common ...
+  -L/usr/local/igl/libigl/external/embree/bin -lrtcore -lsys ... 
+  mexFunction.cpp parse_rhs.cpp;
+

examples/ambient-occlusion-mex/example.m

@@ -0,0 +1,30 @@
+[V,F] = load_mesh('../shared/cheburashka.obj');
+N = per_vertex_normals(V,F);
+S = ambient_occlusion(V,F,V,N,1000);
+L = cotmatrix(V,F);
+M = massmatrix(V,F,'voronoi');
+[EV,~] = eigs(-L,M,15,'sm');
+Z = EV(:,7);
+qZ = round(255*(Z - min(Z))/(max(Z)-min(Z)));
+C = ind2rgb(qZ,jet(256));
+
+nsp = 3;
+fs = 20;
+subplot(nsp,1,1);
+set(tsurf(F,V),'CData',Z,'EdgeColor','none',fphong);
+title('Pseudo-color  ','FontSize',fs);
+axis equal;
+view(2);
+colormap(jet(256));
+
+subplot(nsp,1,2);
+set(tsurf(F,V),'CData',permute(1-[S S S],[1 3 2]),'EdgeColor','none',fphong);
+title('Inverted ambient occlusion  ','FontSize',fs);
+axis equal;
+view(2);
+
+subplot(nsp,1,3);
+set(tsurf(F,V),'CData',bsxfun(@times,1-S,C),'EdgeColor','none',fphong);
+title('Product  ','FontSize',fs);
+axis equal;
+view(2);

examples/ambient-occlusion-mex/mexFunction.cpp

@@ -0,0 +1,56 @@
+#include "parse_rhs.h"
+
+#include <igl/matlab/mexStream.h>
+#include <igl/matlab/MatlabWorkspace.h>
+#include <igl/embree/EmbreeIntersector.h>
+#include <igl/embree/ambient_occlusion.h>
+#include <igl/matlab_format.h>
+
+#include <igl/read.h>
+#include <igl/per_vertex_normals.h>
+
+#include <mex.h>
+
+#include <iostream>
+#include <string>
+
+void mexFunction(int nlhs, mxArray *plhs[], 
+    int nrhs, const mxArray *prhs[])
+{
+  // This is useful for debugging whether Matlab is caching the mex binary
+  //mexPrintf("%s %s\n",__TIME__,__DATE__);
+  igl::mexStream mout;
+  std::streambuf *outbuf = std::cout.rdbuf(&mout);
+
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+
+  MatrixXd V,P,N;
+  VectorXd S;
+  MatrixXi F;
+  int num_samples;
+  parse_rhs(nrhs,prhs,V,F,P,N,num_samples);
+  // Prepare left-hand side
+  nlhs = 1;
+
+  //read("../shared/cheburashka.obj",V,F);
+  //P = V;
+  //per_vertex_normals(V,F,N);
+  EmbreeIntersector<Eigen::MatrixXd,Eigen::MatrixXi,Eigen::Vector3d> ei;
+  ei = EmbreeIntersector<MatrixXd,MatrixXi,Vector3d>(V,F);
+  ambient_occlusion(ei,P,N,num_samples,S);
+  MatlabWorkspace mw;
+  mw.save(V,"V");
+  mw.save(P,"P");
+  mw.save(N,"N");
+  mw.save_index(F,"F");
+  mw.save(S,"S");
+  mw.write("out.mat");
+
+  plhs[0] = mxCreateDoubleMatrix(S.rows(),S.cols(), mxREAL);
+  copy(S.data(),S.data()+S.size(),mxGetPr(plhs[0]));
+
+  // Restore the std stream buffer Important!
+  std::cout.rdbuf(outbuf);
+}

examples/ambient-occlusion-mex/parse_rhs.cpp

@@ -0,0 +1,59 @@
+#ifdef MEX
+#include "parse_rhs.h"
+#include <algorithm>
+#include <functional>
+
+void parse_rhs(
+  const int nrhs, 
+  const mxArray *prhs[], 
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F,
+  Eigen::MatrixXd & P,
+  Eigen::MatrixXd & N,
+  int & num_samples)
+{
+  using namespace std;
+  if(nrhs < 5)
+  {
+    mexErrMsgTxt("nrhs < 5");
+  }
+
+  const int dim = mxGetN(prhs[0]);
+  if(dim != 3)
+  {
+    mexErrMsgTxt("Mesh vertex list must be #V by 3 list of vertex positions");
+  }
+  if(dim != (int)mxGetN(prhs[1]))
+  {
+   mexErrMsgTxt("Mesh facet size must be 3");
+  }
+  if(mxGetN(prhs[2]) != dim)
+  {
+    mexErrMsgTxt("Point list must be #P by 3 list of origin locations");
+  }
+  if(mxGetN(prhs[3]) != dim)
+  {
+    mexErrMsgTxt("Normal list must be #P by 3 list of origin normals");
+  }
+  if(mxGetN(prhs[4]) != 1 || mxGetM(prhs[4]) != 1)
+  {
+    mexErrMsgTxt("Number of samples must be scalar.");
+  }
+
+
+  V.resize(mxGetM(prhs[0]),mxGetN(prhs[0]));
+  copy(mxGetPr(prhs[0]),mxGetPr(prhs[0])+V.size(),V.data());
+  F.resize(mxGetM(prhs[1]),mxGetN(prhs[1]));
+  copy(mxGetPr(prhs[1]),mxGetPr(prhs[1])+F.size(),F.data());
+  F.array() -= 1;
+  P.resize(mxGetM(prhs[2]),mxGetN(prhs[2]));
+  copy(mxGetPr(prhs[2]),mxGetPr(prhs[2])+P.size(),P.data());
+  N.resize(mxGetM(prhs[3]),mxGetN(prhs[3]));
+  copy(mxGetPr(prhs[3]),mxGetPr(prhs[3])+N.size(),N.data());
+  if(*mxGetPr(prhs[4]) != (int)*mxGetPr(prhs[4]))
+  {
+    mexErrMsgTxt("Number of samples should be non negative integer.");
+  }
+  num_samples = (int) *mxGetPr(prhs[4]);
+}
+#endif

examples/ambient-occlusion-mex/parse_rhs.h

@@ -0,0 +1,27 @@
+#ifndef PARSE_RHS_H
+#define PARSE_RHS_H
+
+#include <mex.h>
+#include <Eigen/Core>
+// Parse right hand side arguments for a matlab mex function.
+//
+// Inputs:
+//   nrhs  number of right hand side arguments
+//   prhs  pointer to right hand side arguments
+// Outputs:
+//   V  n by dim list of mesh vertex positions
+//   F  m by dim list of mesh face indices
+//   P  #P by 3 list of origin points
+//   N  #P by 3 list of origin normals
+//   num_samples  number of samples
+// Throws matlab errors if dimensions are not sane.
+void parse_rhs(
+  const int nrhs, 
+  const mxArray *prhs[], 
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F,
+  Eigen::MatrixXd & P,
+  Eigen::MatrixXd & N,
+  int & num_samples);
+
+#endif 

include/igl/embree/ambient_occlusion.cpp

@@ -49,7 +49,31 @@ void igl::ambient_occlusion(
   }
 }
 
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+void igl::ambient_occlusion(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  using namespace igl;
+  using namespace Eigen;
+  EmbreeIntersector<
+    PlainObjectBase<DerivedV>,
+    PlainObjectBase<DerivedF>,
+    Matrix<typename DerivedP::Scalar,3,1> > ei(V,F);
+  return ambient_occlusion(ei,P,N,num_samples,S);
+}
+
 #ifndef IGL_HEADER_ONLY
 // Explicit template instanciation
 template void igl::ambient_occlusion<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, const int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::ambient_occlusion<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/embree/ambient_occlusion.h

@@ -33,8 +33,21 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedN> & N,
     const int num_samples,
     Eigen::PlainObjectBase<DerivedS> & S);
-
-
+  // Wrapper which builds new EmbreeIntersector for (V,F). That's expensive so
+  // avoid this if repeatedly calling.
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  void ambient_occlusion(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
 };
 #ifdef IGL_HEADER_ONLY
 #  include "ambient_occlusion.cpp"