VTK/VTKMatlab

Mini HOWTO contributed by Dominik Szczerba, Computer Vision Lab, ETH --Domel 04:18, 7 July 2006 (EDT)

Users may also find this link interesting Matlab ITK interface

---

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,1);//save dcolors
  }
  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();

  // 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);
  
  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/vtk-5.1-CVS/include/vtk-5.1/ vtkPolyDataRenderer.cpp -L/home/domi/local/pack/vtk-5.1-CVS/lib -lvtkFiltering -lvtkRendering

You can easilly create an appropriate config file, e.g., .matlab/mexopts_vtk.sh, put all include/lib stuff there, and call mex as e.g.:

>> mex -O -f /home/domi/.matlab/mexopts_vtk.sh vtkPolyDataRenderer.cpp

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);

To exit the viewer hit 'q' over the window (don't close the window itself).

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 has an annoying property of not outputting its content immediately, only after the main() returns. The Mathworks somehow remains blind and deaf to these problems.

---

VTK: [Welcome | Site Map]