VTK/VTKMatlab

Mini HOWTO contributed by Dominik Szczerba, Computer Vision Lab, ETH

--Domel 04:18, 7 July 2006 (EDT)

It is possible (and highly desirable) to provide VTK bindings for Matlab. For that to happen one has to write a C++ file to be compiled with mex compiler right from under Matlab. Worked example:

/****************************************************
(C) 2005 Dominik Szczerba <domi at vision ee ethz ch>
*****************************************************/

#include <cstdlib>
#include <cmath>

/* MEX HEADERS */
#include <mex.h>

/* VTK HEADERS */
#include <vtkPolyData.h>
#include <vtkPoints.h>
#include <vtkCellArray.h>
#include <vtkDoubleArray.h>
#include <vtkPolyData.h>

#include <vtkCell.h>
#include <vtkCellData.h>
#include <vtkDataSet.h>
#include <vtkDataSetAttributes.h>
#include <vtkProperty.h>

#include <vtkDataSetMapper.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkInteractorStyleTrackballCamera.h>

/* MAIN FUNCTION */
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[]){

  double *v;	// Vertices
  double *f;	// Faces
  
  unsigned int nV, nF, nC=0, i, j;

  if(nrhs<2)
	 mexErrMsgTxt("2+ input args required");

  double* dcolors = 0;
  vtkDoubleArray* dcolarr = 0;
  
  nV = mxGetM(prhs[0]);
  nF = mxGetM(prhs[1]);
  v = mxGetPr(prhs[0]);
  f = (double*)mxGetPr(prhs[1]);
  if(nrhs>=3){
	 //mexPrintf("found color information\n");
	 nC = mxGetM(prhs[2]);
	 if(nC!=nF)
		mexErrMsgTxt("dimension mismatch");
	 dcolors = mxGetPr(prhs[2]);
	 dcolarr = vtkDoubleArray::New();
	 dcolarr->SetName("aCellScalar");
	 dcolarr->SetArray(dcolors,nF,0);
  }
  mexPrintf("got %d nodes, %d faces and %d colors\n",nV,nF,nC);

  // Initialize vtk variables
  vtkPolyData *surf = vtkPolyData::New();
  vtkPoints *sverts = vtkPoints::New();
  vtkCellArray *sfaces = vtkCellArray::New();

//   for (i=0; i<nV; i++) mexPrintf("%d %f %f %f\n",i,v[i],v[nV+i],v[2*nV+i]);
//   for (i=0; i<nF; i++){
// 	 for(j=0;j<3;j++)
// 		mexPrintf("%d ",f[j*nF+i]);
//  	 mexPrintf("\n");
//   }

  // Load the point, cell, and data attributes.
  //indexes in VTK are 0-based, in Matlab 1-based
  const int offset = -1;
  for(i=0; i<nV; i++) sverts->InsertPoint(i,v[i],v[nV+i],v[2*nV+i]);
  for(i=0; i<nF; i++){
	 sfaces->InsertNextCell(3);
	 for(j=0;j<3;j++){
		int idx = (unsigned int)f[j*nF+i]+offset;
		if(idx<0 || idx>=nV)
		//mexPrintf("Warning : wrong index\n");
		  mexErrMsgTxt("wrong index");
		sfaces->InsertCellPoint(vtkIdType(idx));
	 }
  }
  
  // assign the pieces to the vtkPolyData.
  surf->SetPoints(sverts);
  surf->SetPolys(sfaces);
  sverts->Delete();
  sfaces->Delete();

  if(nC>0){
	 surf->GetCellData()->SetScalars(dcolarr);
  }

  // now render
  vtkDataSetMapper::SetResolveCoincidentTopologyToPolygonOffset();
  vtkPolyDataMapper *surfMapper = vtkPolyDataMapper::New();
  vtkActor *surfActor = vtkActor::New();
  vtkRenderer *aRenderer = vtkRenderer::New();
  vtkRenderWindow *renWin = vtkRenderWindow::New();
  vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
  
  surfMapper->SetInput(surf);
  //surfMapper->ScalarVisibilityOff();
  if(nC>0){
	 surfMapper->ScalarVisibilityOn();
	 surfMapper->SetScalarRange(0,1);  
  }

  surfActor->SetMapper(surfMapper);
  surfActor->GetProperty()->SetAmbient(0.1);
  surfActor->GetProperty()->SetDiffuse(0.3);
  surfActor->GetProperty()->SetSpecular(1.0);
  surfActor->GetProperty()->SetSpecularPower(30.0);
  
  aRenderer->AddActor(surfActor);
  aRenderer->SetBackground(0.2,0.2,0.2);
  //needed after camera moves
  //aRenderer->ResetCameraClippingRange ();
  
  iren->SetRenderWindow(renWin);
  
  renWin->AddRenderer(aRenderer);
  renWin->SetSize(500,500);
  renWin->Render();
  
  vtkInteractorStyleTrackballCamera *style = 
 	 vtkInteractorStyleTrackballCamera::New();
  iren->SetInteractorStyle(style);

  // Start interactive rendering
  iren->Start();

  // clean up
  surfMapper->Delete();
  surfActor->Delete();
  surf->Delete();
  aRenderer->Delete();
  renWin->Delete();
  iren->Delete();
  style->Delete();

}

Compile it under matlab with:

> mex -I<includesToVtkHeaders> fileToCompile.cpp -L<libDir> -lAppropriateLibraries

Under linux this is normal compiling stuff, in doubt just look up your own running VTK projects' makefiles to provide includes and linking directives, see also GCC documentation. In my case (latest VTK-5 via CVS built with shared libs) the following worked (change paths where appropriate):

>> mex -O -I/home/domi/local/pack/matlab-7.2.R2006a/extern/include -I/home/domi/local/pack/vtk-5.1-CVS/include/vtk-5.1/ vtkPolyDataRenderer.cpp -L/home/domi/local/pack/vtk-5.1-CVS/lib -lvtkFiltering -lvtkRendering

Under Windows this is more trouble. First off, you need a reasonable C++ compiler - I use the well known MinGW (http://www.mingw.org/). Install, learn how to use it, it will take a few moments. Download CMake and VTK sources and compile VTK yourself (version 5.0.0 compiled silently). The linux way I am used to work with cmake (i.e. cmake src_dir called from a build_dir) failed rather violently, with Microsoft apologizing for inconvenience. You have to do something like:

/c/Program\ Files/CMake\ 2.4/bin/CMakeSetup

from the MSYS console (in the build directory) and configure everything there.

Then you need something to compile the above given VTK mex example - not surprisingly, the Matlab C compiler didnt work out for me. I used gnumex (http://gnumex.sourceforge.net/) - a way to use Mingw with Matlab through mex. Install and follow documentation therein. I eventually had to change the mexopts to use g++ instead of gcc, otherwise the C++ standard stuff (eg. cout etc.) was not available. Also, beware of spaces in filenames everywhere. Avoid wherever possible, otherwise use 'quotes'. Under Windows I compiled VTK with static libs and could find no way of using -L/-l directives. If you find out let me know. I ended up doing it explicitely:

>> mex -ID:\MyProjects\VTK-5.0\include\vtk-5.0 vtkPolyDataRenderer.cxx D:\MyProjects\VTK-5.0\lib\libvtkRendering.a D:\MyProjects\VTK-5.0\lib\libvtkIO.a D:\MyProjects\VTK-5.0\lib\libvtkGraphics.a D:\MyProjects\VTK-5.0\lib\libvtkImaging.a D:\MyProjects\VTK-5.0\lib\libvtkftgl.a D:\MyProjects\VTK-5.0\lib\libvtkfreetype.a C:\mingw\lib\libopengl32.a D:\MyProjects\VTK-5.0\lib\libvtkFiltering.a D:\MyProjects\VTK-5.0\lib\libvtkCommon.a D:\MyProjects\VTK-5.0\lib\libvtksys.a D:\MyProjects\VTK-5.0\lib\libvtkDICOMParser.a D:\MyProjects\VTK-5.0\lib\libvtkpng.a D:\MyProjects\VTK-5.0\lib\libvtktiff.a D:\MyProjects\VTK-5.0\lib\libvtkzlib.a D:\MyProjects\VTK-5.0\lib\libvtkjpeg.a D:\MyProjects\VTK-5.0\lib\libvtkexpat.a C:\mingw\lib\libvfw32.a D:\MyProjects\VTK-5.0\lib\libvtkMPEG2Encode.a C:\mingw\lib\libgdi32.a

which looks quite hairy, but actually works fine. Change paths where appropriate, think of smarter syntax (get -L -l stuff to work or modify your mexopts.bat) Refer to mex documentation for more information.

And finally, here is an example how to use the above given mex file, once compiled:

[x y z v] = flow;
q = z./x.*y.^3;
mesh=isosurface(x, y, z, q, -.08, v);
[N,M]=size(mesh.faces);
scalars=[0:1:N-1]'/(N-1);
vtkPolyDataRenderer(mesh.vertices,mesh.faces,scalars);

Possible problems:

1) Matlab may die or even your computer may crash on a buggy opengl driver. It was a case with me under Debian/stable, update to Debain/testing fixed the problem (apparently the driver for my ATI gfx card was updated)

2) Matlab may crash on errors in the mex file, e.g. errors in pointer stuff. Where you get a normal VTK application crash you will get the whole Matlab crashing.

3) Matlab mex files cannot be easilly interrupted (e.g. with CTRL+C). On linux with matlab run in terminal (with -nodesktop) this can be remedied by intercepting signals, but isint easy and limited not only to linux but also -nodesktop usage. Also, printf() will not work, you need the funny mexPrintf() instead, which doesnt output its text immediately, only after the main() returns. This is annoying and the Mathworks somehow remains blind and deaf to these problems.