Posts Tagged ‘mex’

Mex wrapper for graph segmentation

Thursday, May 4th, 2017

I wrote a small mex wrapper for the graph segmentation part of the “Graph Based Image Segmentation” code. Most of the previously matlab implementations/wrappers worked on images. I want to apply this to geometry so I needed access to the graph segmentation directly. Here’s the wrapper (soon to be part of gptoolbox):

// mexopts = gptoolbox_mexopts('Static',false,'Debug',true);
// mex('segment_graph.cpp',mexopts{:});
#ifdef MEX
#  include <mex.h>
#  include <igl/C_STR.h>
#  include <igl/matlab/mexErrMsgTxt.h>
#  undef assert
#  define assert( isOK ) ( (isOK) ? (void)0 : (void) ::mexErrMsgTxt(C_STR(__FILE__<<":"<<__LINE__<<": failed assertion `"<<#isOK<<"'"<<std::endl) ) )
#endif
#include "segment-graph.h"
#include <igl/matlab/mexErrMsgTxt.h>
#include <igl/matlab/parse_rhs.h>
#include <igl/unique.h>
#include <igl/matlab/prepare_lhs.h>
#include <igl/matlab/requires_arg.h>
#include <igl/matlab/validate_arg.h>
#include <igl/matlab/MexStream.h>
#include <Eigen/Sparse>
void mexFunction(
 int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  using namespace igl::matlab;
  using namespace Eigen;
  using namespace std;
  igl::matlab::MexStream mout;        
  std::streambuf *outbuf = std::cout.rdbuf(&mout);


  mexErrMsgTxt(nrhs>0,"Too few inputs");
  mexErrMsgTxt(mxIsSparse(prhs[0]),"Matrix should be sparse");
  const mxArray * mx_data = prhs[0];
  const int m = mxGetM(mx_data);
  const int n = mxGetN(mx_data);
  mexErrMsgTxt(n == mxGetM(prhs[0]), "Matrix should be square");
  assert(mxIsSparse(mx_data));
  assert(mxGetNumberOfDimensions(mx_data) == 2);
  // TODO: It should be possible to directly load the data into the sparse
  // matrix without going through the triplets
  // Copy data immediately
  double * pr = mxGetPr(mx_data);
  mwIndex * ir = mxGetIr(mx_data);
  mwIndex * jc = mxGetJc(mx_data);
  const int num_edges = mxGetNzmax(mx_data);
  edge * edges = new edge[num_edges];
  int k = 0;
  for(int j=0; j<n;j++)
  {
    // Iterate over inside
    while(k<(int)jc[j+1])
    {
      //cout<<ir[k]<<" "<<j<<" "<<pr[k]<<endl;
      assert((int)ir[k]<m);
      assert((int)j<n);
      edges[k].a = ir[k];
      edges[k].b = j;
      edges[k].w = pr[k];
      k++;
    }
  }

  // defaults 
  int min_size = 0;
  // Threshold 
  int c = sqrt((double)n);
  {
    int i = 1;
    while(i<nrhs)
    {
      mexErrMsgTxt(mxIsChar(prhs[i]),"Parameter names should be strings");
      // Cast to char
      const char * name = mxArrayToString(prhs[i]);
      if(strcmp("Threshold",name) == 0)
      {
        requires_arg(i,nrhs,name);
        validate_arg_scalar(i,nrhs,prhs,name);
        validate_arg_double(i,nrhs,prhs,name);
        c = (double)*mxGetPr(prhs[++i]);
      }else if(strcmp("MinSize",name) == 0)
      {
        requires_arg(i,nrhs,name);
        validate_arg_scalar(i,nrhs,prhs,name);
        validate_arg_double(i,nrhs,prhs,name);
        min_size = (int)((double)*mxGetPr(prhs[++i]));
      }
      i++;
    }
  }

  universe *u = segment_graph(n, num_edges, edges, c);

  // post process small components
  for (int i = 0; i < num_edges; i++) {
    int a = u->find(edges[i].a);
    int b = u->find(edges[i].b);
    if ((a != b) && ((u->size(a) < min_size) || (u->size(b) < min_size)))
      u->join(a, b);
  }

  switch(nlhs)
  {
    case 1:
    {
      plhs[0] = mxCreateDoubleMatrix(m,1, mxREAL);
      Eigen::VectorXi C(m);
      for(int i = 0;i<m;i++)
      {
        C(i) = u->find(i);
      }
      Eigen::VectorXi uC,I,J;
      igl::unique(C,uC,I,J);
      prepare_lhs_index(J,plhs);
    }
    default: break;
  }

  delete[] edges;
  delete u;
  std::cout.rdbuf(outbuf);
}

It takes the graph as a sparse matrix and outputs the component ids:

C = segment_graph(A);

Using glfw background window in matlab mex thread

Monday, April 3rd, 2017

A while ago, I tried to get glfw playing nicely with matlab’s mex files. I didn’t quite succeed.

My current project requires “GP-GPU” programming. I’m currently using glfw’s “background window” feature to create an OpenGL context. My first tests show that matlab’s mex will play nicely if glfwTerminate() is never called:

#include <igl/opengl/glfw/background_window.h>
#include <mex.h>
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  GLFWwindow * window;
  igl::opengl::glfw::background_window(window);
  glfwDestroyWindow(window);
  //glfwTerminate();
}

You can compile this (on mac os x) with:

mex( ...
  '-largeArrayDims','-DMEX','CXXFLAGS=$CXXFLAGS -std=c++11', ...
  'LDFLAGS=\$LDFLAGS -framework Foundation -framework AppKit -framework Accelerate -framework OpenGL -framework AppKit -framework Carbon -framework QuartzCore -framework IOKit ', ...
  '-I/usr/local/igl/libigl//external/nanogui/ext/eigen/', ...
  '-I/usr/local/igl/libigl/include', ...
  '-I/usr/local/igl/libigl/external/nanogui/ext/glfw/include/', ...
  '-L/usr/local/igl/libigl/lib/','-lglfw3', ...
  'background_glfw.cpp');

I don’t know why including glfwTerminate() causes matlab to sometimes crash. The error report is impossible to read and seems to have something to do with Quartz.

Friday, January 6th, 2017

Like clockwork, my matlab updates have gotten out of sync with Xcode updates. It seems like fixing this SDK error always requires a different hack each year. This year I got the error:

No supported compiler or SDK was found. For options, visit
http://www.mathworks.com/support/compilers/current_release/. 

To fix it, I replaced all occurrences of 10.9 with 10.11 in /Applications/MATLAB_R2017a.app/bin/maci64/mexopts/clang{++,}_maci64.xml

I’m still getting linker warnings:

ld: warning: object file was built for newer OSX version (10.11) than being linked (10.9)

For now, I’m assuming that I can ignore them. We’ll see how far that gets me.

Fix dyld linker errors when installing new mosek toolbox

Monday, October 3rd, 2016

Each time I upgrade my mosek library, matlab panics and can’t find it. The old solution was to monkey with the DYLD_LIBRARY_PATH in all sorts of funny places: ~/.profile, /etc/launch.d, ~/Library/LaunchAgents/environment.plist

These all worked at some point, but as far as I can tell, no longer do. They’re also the wrong way to install libraries.

Fortunately, mosek has made life easy. Just

cd /usr/local/mosek/8/tools/platform/osx64x86/bin/
python install.py

This will actually fix all of the binaries and mex files using otool and install_name_tool to find dynamic libraries in their installed locations.

dyld symbol not found, matlab c++

Sunday, December 27th, 2015

I upgrade to the newest Matlab 2015b and suddenly got some unexpected runtime errors in my command line programs compiled with matlab:

dyld: Symbol not found: _mxDestroyArray
  Referenced from: /Users/ajx/./.main
  Expected in: /usr/lib/libSystem.B.dylib
 in /Users/ajx/./.main

I finally tracked this down to a problem in my DYLD_FALLBACK_LIBRARY_PATH. The path to matlab was there but after /usr/lib/ and for some reason that order was throwing off the dynamic linker. I put the matlab path first and that fixed the problem.

GLFW c++ mesh viewer from matlab with Core OpenGL

Monday, June 15th, 2015

Here’s a proof of concept mex function for matlab that passes a triangle mesh and renders it in a new window using glfw and native opengl:

// make sure the modern opengl headers are included before any others
#include <OpenGL/gl3.h>
#define __gl_h_

#include <igl/frustum.h>
#include <igl/read_triangle_mesh.h>
#include <igl/get_seconds.h>
#include <igl/matlab/mexStream.h>
#include <igl/matlab/parse_rhs.h>
#include <igl/matlab/mexErrMsgTxt.h>
#include <Eigen/Core>
#define GLFW_INCLUDE_GLU
#include <GLFW/glfw3.h>
#include <mex.h>

#include <string>
#include <chrono>
#include <thread>
#include <iostream>

std::string vertex_shader = R"(
#version 330 core
uniform mat4 proj;
uniform mat4 model;
in vec3 position;
void main()
{
  gl_Position = proj * model * vec4(position,1.);
}
)";
std::string fragment_shader = R"(
#version 330 core
out vec3 color;
void main()
{
  color = vec3(0.8,0.4,0.0);
}
)";

// width, height, shader id, vertex array object
int w=800,h=600;
double highdpi=1;
GLuint prog_id=0;
GLuint VAO;
// Mesh data: RowMajor is important to directly use in OpenGL
Eigen::Matrix< float,Eigen::Dynamic,3,Eigen::RowMajor> V;
Eigen::Matrix<GLuint,Eigen::Dynamic,3,Eigen::RowMajor> F;
void mexFunction(int nlhs, mxArray *plhs[], 
   int nrhs, const mxArray *prhs[])
//#define mexErrMsgTxt(X) std::cerr<<X<<std::endl;exit(EXIT_FAILURE);
//int main()
{
  using namespace std;
  igl::MexStream mout;        
  std::streambuf *outbuf = cout.rdbuf(&mout);
  // Init glfw and create window + OpenGL context
  if(!glfwInit())
  {
    mexErrMsgTxt("Could not initialize glfw");
  }
  const auto & error = []
    (int error, const char* description){
      mexErrMsgTxt(description);
    };
  glfwSetErrorCallback(error);
  glfwWindowHint(GLFW_SAMPLES, 4);
  glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
  glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 2);
  glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
  glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
  GLFWwindow* window = glfwCreateWindow(w, h, "OpenGL in matlab? ! ? !", NULL, NULL);
  if(!window)
  {
    glfwTerminate();
    mexErrMsgTxt("Could not create glfw window");
  }

  glfwMakeContextCurrent(window);

      int major, minor, rev;
      major = glfwGetWindowAttrib(window, GLFW_CONTEXT_VERSION_MAJOR);
      minor = glfwGetWindowAttrib(window, GLFW_CONTEXT_VERSION_MINOR);
      rev = glfwGetWindowAttrib(window, GLFW_CONTEXT_REVISION);
      printf("OpenGL version recieved: %d.%d.%d\n", major, minor, rev);
      printf("Supported OpenGL is %s\n", (const char*)glGetString(GL_VERSION));
      printf("Supported GLSL is %s\n", (const char*)glGetString(GL_SHADING_LANGUAGE_VERSION));

  glfwSetInputMode(window,GLFW_CURSOR,GLFW_CURSOR_NORMAL);
  const auto & reshape = [] (GLFWwindow* window, int w, int h)
  {
    ::w=w,::h=h;
  };
  glfwSetWindowSizeCallback(window,reshape);

  {
    int width, height;
    glfwGetFramebufferSize(window, &width, &height);
    int width_window, height_window;
    glfwGetWindowSize(window, &width_window, &height_window);
    highdpi = width/width_window;
    reshape(window,width_window,height_window);
  }


  // Compile each shader
  const auto & compile_shader = [](const GLint type,const char * str) -> GLuint
  {
    GLuint id = glCreateShader(type);
    glShaderSource(id,1,&str,NULL);
    glCompileShader(id);
    return id;
  };
  GLuint vid = compile_shader(GL_VERTEX_SHADER,vertex_shader.c_str());
  GLuint fid = compile_shader(GL_FRAGMENT_SHADER,fragment_shader.c_str());
  // attach shaders and link
  prog_id = glCreateProgram();
  glAttachShader(prog_id,vid);
  glAttachShader(prog_id,fid);
  glLinkProgram(prog_id);
  GLint status;
  glGetProgramiv(prog_id, GL_LINK_STATUS, &status);
  glDeleteShader(vid);
  glDeleteShader(fid);

  // Read input mesh from rhs
  const int dim = mxGetN(prhs[1]);
  igl::mexErrMsgTxt(dim == 3,
    "Mesh vertex list must be #V by 3 list of vertex positions");
  igl::mexErrMsgTxt(dim == mxGetN(prhs[0]),
    "Mesh \"face\" simplex size must equal dimension");
  igl::parse_rhs_double(prhs+1,V);
  igl::parse_rhs_index(prhs+0,F);

  V.rowwise() -= V.colwise().mean();
  V /= (V.colwise().maxCoeff()-V.colwise().minCoeff()).maxCoeff();
  V /= 1.2;

  // Generate and attach buffers to vertex array
  glGenVertexArrays(1, &VAO);
  GLuint VBO, EBO;
  glGenBuffers(1, &VBO);
  glGenBuffers(1, &EBO);
  glBindVertexArray(VAO);
  glBindBuffer(GL_ARRAY_BUFFER, VBO);
  glBufferData(GL_ARRAY_BUFFER, sizeof(float)*V.size(), V.data(), GL_STATIC_DRAW);
  glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, EBO);
  glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(GLuint)*F.size(), F.data(), GL_STATIC_DRAW);
  glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(GLfloat), (GLvoid*)0);
  glEnableVertexAttribArray(0);
  glBindBuffer(GL_ARRAY_BUFFER, 0); 
  glBindVertexArray(0);

  // Main display routine
  while (!glfwWindowShouldClose(window))
  {
      double tic = igl::get_seconds();
      // clear screen and set viewport
      glClearColor(0.0,0.4,0.7,0.);
      glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
      glViewport(0,0,w,h);

      // Projection and modelview matrices
      Eigen::Matrix4f proj = Eigen::Matrix4f::Identity();
      float near = 0.01;
      float far = 100;
      float top = tan(35./360.*M_PI)*near;
      float right = top * (double)::w/(double)::h;
      igl::frustum(-right,right,-top,top,near,far,proj);
      Eigen::Affine3f model = Eigen::Affine3f::Identity();
      model.translate(Eigen::Vector3f(0,0,-1.5));
      // spin around
      static size_t count = 0;
      model.rotate(Eigen::AngleAxisf(0.005*count++,Eigen::Vector3f(0,1,0)));

      // select program and attach uniforms
      glUseProgram(prog_id);
      GLint proj_loc = glGetUniformLocation(prog_id,"proj");
      glUniformMatrix4fv(proj_loc,1,GL_FALSE,proj.data());
      GLint model_loc = glGetUniformLocation(prog_id,"model");
      glUniformMatrix4fv(model_loc,1,GL_FALSE,model.matrix().data());

      // Draw mesh as wireframe
      glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
      glBindVertexArray(VAO);
      glDrawElements(GL_TRIANGLES, F.size(), GL_UNSIGNED_INT, 0);
      glBindVertexArray(0);

      glfwSwapBuffers(window);

      {
        glfwPollEvents();
        // In microseconds
        double duration = 1000000.*(igl::get_seconds()-tic);
        const double min_duration = 1000000./60.;
        if(duration<min_duration)
        {
          std::this_thread::sleep_for(std::chrono::microseconds((int)(min_duration-duration)));
        }
      }
  }
  glfwDestroyWindow(window);
  glfwTerminate();
  // Restore the std stream buffer Important!
  std::cout.rdbuf(outbuf);
}

Compile with something like:

mex('LDOPTIMFLAGS=','CXX=clang++','LD=clang++','-DMEX -v -largeArrayDims','CXXFLAGS=$CXXFLAGS -std=c++11','-I/usr/local/igl/libigl/include','-I/usr/local/igl/libigl/external/glfw/include','-L/usr/local/igl/libigl/external/glfw/lib -lglfw3_clang','-I/opt/local/include/eigen3/','LDFLAGS=\$LDFLAGS -framework Foundation -framework AppKit -framework Accelerate -framework GLUT -framework OpenGL -framework IOKit -framework Cocoa -framework CoreVideo', 'tsurf_opengl.cpp');

And then if you issue:

[F,V] = surf2patch(surface(16*membrane(1,10)),'triangles');
V = V*axisangle2matrix([1 0 0],pi*0.5)*axisangle2matrix([0 1 0],-0.25*pi)*axisangle2matrix([1 0 0],-pi*0.1);
tsurf_opengl(F,V);

You should see something like:

opengl matlab glfw There are still quite a few kinks, but it’s at least a start.

Paste directly from clipboard into matlab image array

Sunday, October 12th, 2014

I find myself often wanting to experiment with images from papers I’m reading. To load the image into matlab properly I should extract it from the pdf, save it in a file and load the file via matlab. Often I skip the first step and just take a screengrab. But I still need to create a dummy file just to load it into matlab.

Here’s a mex function to load images directly from the clipboard. It even maintains all color channels (RGB, CMYK, RGBA, etc.)

Save this in a file paste.h:

#include <cstddef>
// Get image from clipboard as an RGB image
//
// Outputs:
//   IM  h*w*c list of rgb pixel color values as uint8, running over height,
//     then width, then rgb channels
//   h  height
//   w  width
//   c  channels
bool paste(unsigned char *& IM,size_t & h, size_t & w, size_t & c);

and in paste.mm:

#import "paste.h"
#import <Foundation/Foundation.h>
#import <Cocoa/Cocoa.h>
#import <unistd.h>

bool paste(unsigned char *& IM,size_t & h, size_t & w, size_t & c)
{
  NSPasteboard *pasteboard = [NSPasteboard generalPasteboard];
  NSArray *classArray = [NSArray arrayWithObject:[NSImage class]];
  NSDictionary *options = [NSDictionary dictionary];
  BOOL ok = [pasteboard canReadObjectForClasses:classArray options:options]; 
  if(!ok)
  {
    //printf("Error: clipboard doesn't seem to contain image...");
    return false;
  }
  NSArray *objectsToPaste = [pasteboard readObjectsForClasses:classArray options:options];
  NSImage *image = [objectsToPaste objectAtIndex:0];
  NSBitmapImageRep* bitmap = [NSBitmapImageRep imageRepWithData:[image TIFFRepresentation]];
  // http://stackoverflow.com/a/19649616/148668
  w = [bitmap pixelsWide];
  h = [bitmap pixelsHigh];
  size_t rowBytes = [bitmap bytesPerRow];
  c = rowBytes/w;
  unsigned char* pixels = [bitmap bitmapData];
  IM = new unsigned char[w*h*c];
  for(size_t y = 0; y < h ; y++)
  {
    for(size_t x = 0; x < w; x++)
    {
      for(size_t d = 0;d<c;d++)
      {
        // For some reason colors are invertex
        IM[y+h*(x+d*w)] = pixels[y*rowBytes + x*c + d];
      }
    }
  }
  [image release];
  return true;
}

and in impaste.cpp

#ifdef MEX
#include <mex.h>
#include "paste.h"
#include <iostream>

void mexFunction(
  int nlhs, mxArray *plhs[], 
  int nrhs, const mxArray *prhs[])
{
  unsigned char * IM;
  size_t h,w,c;
  if(!paste(IM,h,w,c))
  {
    mexErrMsgTxt("Clipboard doesn't contain image.");
  }
  switch(nlhs)
  {
    default:
    {
      mexErrMsgTxt("Too many output parameters.");
    }
    case 1:
    {
      size_t dims[] = {h,w,c};
      plhs[0] = mxCreateNumericArray(3,dims,mxUINT8_CLASS,mxREAL);
      unsigned char * IMp = (unsigned char *)mxGetData(plhs[0]);
      std::copy(IM,IM+w*h*c,IMp);
      // Fall through
    }
    case 0: break;
  }
}

#endif

Then compile on Mac OS X using something like:

mex -v -largeArrayDims -DMEX CXX='/usr/bin/clang++' LD='/usr/bin/clang++' LDOPTIMFLAGS='-O ' LDFLAGS='\$LDFLAGS -framework Foundation -framework AppKit' -o impaste impaste.cpp paste.mm

Notice, I’m forcing the mex compiler to use clang and include the relevant frameworks.

Then you can call from matlab using something like:

IM = impaste();

or

imshow(impaste());

Compiling pardiso matlab interface finding wrong gfortran dylib

Wednesday, May 28th, 2014

I tried to compile the pardiso mex interface for matlab and got a very annoying linker error upon running:

Invalid MEX-file '/Users/ajx/.../pardisoinit.mexmaci64':
dlopen(/Users/ajx/.../pardisoinit.mexmaci64, 6): Symbol
not found: __gfortran_transfer_array_write
  Referenced from: /usr/local/lib/libpardiso500-MACOS-X86-64.dylib
  Expected in: /Applications/MATLAB_R2014a.app/sys/os/maci64/libgfortran.3.dylib
 in /usr/local/lib/libpardiso500-MACOS-X86-64.dylib

The problem seems to be that matlab is finding gfortran.3.dylib in its own directory path rather than the one that pardisoinit was linked to. These libraries are also unfortunately different. I think there’s an otool/install_name_tool trick for this, but so far I’m using this hack. I replace the libgfortran.3.dylib in Matlab with the one from gcc:

mv /Applications/MATLAB_R2014a.app/sys/os/maci64/libgfortran.3.dylib /Applications/MATLAB_R2014a.app/sys/os/maci64/libgfortran.3.dylib.off
cp /opt/local/lib/libgcc/libgfortran.3.dylib /Applications/MATLAB_R2014a.app/sys/os/maci64/libgfortran.3.dylib

So far this hasn’t messed up anything else in MATLAB, but I guess, without knowing better, it’s only a matter of luck. Too bad pardiso doesn’t seem to support mac and isn’t open source.

Matlab’s mex can’t find std headers

Friday, April 18th, 2014

Trying to compile the ecos mex files I ran into the following errors:

In file included from /opt/local/lib/gcc47/gcc/x86_64-apple-darwin13/4.7.3/include-fixed/syslimits.h:7:0,
                 from /opt/local/lib/gcc47/gcc/x86_64-apple-darwin13/4.7.3/include-fixed/limits.h:34,
                 from ../external/SuiteSparse_config/SuiteSparse_config.h:45,
                 from ../external/ldl/include/ldl.h:11,
                 from ../external/ldl/src/ldl.c:157:
/opt/local/lib/gcc47/gcc/x86_64-apple-darwin13/4.7.3/include-fixed/limits.h:169:61: error: no include path in which to search for limits.h
In file included from ../external/ldl/include/ldl.h:11:0,
                 from ../external/ldl/src/ldl.c:157:
../external/SuiteSparse_config/SuiteSparse_config.h:46:20: fatal error: stdlib.h: No such file or directory
compilation terminated.

    mex: compile of ' "../external/ldl/src/ldl.c"' failed.

Turns out this has nothing to do with ecos, rather I had upgraded my OS since installing matlab. I need to replace 10.7 with 10.9 in my /Applications/MATLAB_R2013b.app/bin/mexopts.sh file. Then everything compiled fine.

The offending flags were:

-isysroot /Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.7.sdk/ -mmacosx-version-min=10.7

And the correct versions were:

-isysroot /Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.9.sdk/ -mmacosx-version-min=10.9

Ambient occlusion for matlab pseudo-color plots

Sunday, October 13th, 2013

I’m often plotting scalar functions on surfaces using the trisurf routine of matlab. This has basic per-vertex, per-face coloring. And with some hassle you could even set up lights. Instead I’ve settled on a better way to convey the shape of the surface and the function by multiply the colors with the inverse of the mesh’s ambient occlusion. I wrote a matlab wrapper mex function for our libigl function. Find it in libigl/examples/ambient-occlusion-mex (version ≥0.3.4). Once compiled you can call from matlab it like this:

S = ambient_occlusion(V,F,P,N,num_samples);

Where (V,F) is your triangle mesh, P are the evaluation points and N are the normals (defining a hemisphere).

Here’s an example of multiplying it against some pseudo-color values:

ambient occlusion over pseudo color plot matlab

For computing ambient occlusion per-vertex and then applying it to per-vertex colors, use:

AO = ambient_occlusion(V,F,V,per_vertex_normals(V,F),1000);
tsurf(F,V,'FaceVertexCData',bsxfun(@times,(1-AO),C),fphong,'EdgeColor','none');

Update: If you have a scalar field and you’d like to visualize it using ambient occlusion on top of a pseudo-color, try:

tsurf(F,V,'FaceVertexCData',bsxfun(@times,1-AO,squeeze(ind2rgb(floor(matrixnormalize(Wbh(:,1))*256),jet(256)))),fphong,'EdgeColor','none');