Posts Tagged ‘matlab’

MATLAB gotcha inverting a (sparse) diagonal matrix

Thursday, November 2nd, 2017

Just got burned by a silly Matlab gotcha. Suppose you have a diagonal matrix D and you want to compute the inverse square root matrix:

Disqrt = diag(1./sqrt(diag(D))

But this will be dense!

Instead use

Disqrt = diag(sqrt(diag(D).^-1)

Or maybe

Disqrt = diag(diag(D).^-0.5)

Not sure if there’s an accuracy difference (hopefully not).

Paper-worthy rendering in MATLAB

Thursday, July 20th, 2017

MATLAB is not a great tool for creating 3D renderings. However, the learning curves for most commercial rendering tools are quite steep. Other tools like Mitsuba can create beautiful pictures, but can feel quite cumbersome for rendering pure geometry rather than the physical scenes their designed for.

Over the years, I’ve developed a way of creating plots of 3D shapes in MATLAB using a few extra functions in gptoolbox. This started as a way to just make images from research prototypes more palatable, but eventually became the usual way that I render images for papers. If the code for my research is already written in MATLAB then one huge advantage is that every image in my paper can have a *.m script that deterministically generates the result and the corresponding image with user intervention. This helps with reproducibility, editing and sharing between collaborators.

Here’s a “VFX Breakdown” of rendering a 3D shape in MATLAB.

t = tsurf(F,V);
set(gcf,'COlor',0.94*[1 1 1]);
teal = [144 216 196]/255;
pink = [254 194 194]/255;
bg_color = pink;
fg_color = teal;
for pass = 1:10
  switch pass
  case 1
    % blank run
    axis([-209.4       119.38      -181.24       262.67      -247.28 247.38]);
  case 2
    axis equal;
    axis([-209.4       119.38      -181.24       262.67      -247.28 247.38]);
    axis vis3d;
  case 3
    t.EdgeColor = 'none';
  case 4
    set(t,fphong,'FaceVertexCData',repmat(fg_color,size(V,1),1));
  case 5
    set(t,fsoft);
  case 6
    l = light('Position',[0.2 -0.2 1]);
  case 7
    set(gca,'Visible','off');
  case 8
    set(gcf,'Color',bg_color);
  case 9
    s = add_shadow(t,l,'Color',bg_color*0.8,'BackgroundColor',bg_color,'Fade','infinite');
  case 10
    apply_ambient_occlusion(t,'AddLights',false,'SoftLighting',false);
  end

  vidObj = VideoWriter(sprintf('nefertiti-%02d.mp4',pass),'MPEG-4');
  vidObj.Quality = 100;
  vidObj.open;
  thetas = linspace(30,-30,450);
  for theta = thetas(1:end-1)
    view(theta,30);
    drawnow;
    vidObj.writeVideo(getframe(gcf));
  end
  vidObj.close;

end

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.

Offset surface of triangle mesh in matlab

Wednesday, March 22nd, 2017

Here’s a little demonstration of how to use gptoolbox and MATLAB to generate an offset surfaces of a triangle mesh. This takes a mesh in V,F and creates a mesh SV,SF of the isosurface at signed distance iso:

% Extract offset at minus 3% of bounind box diagonal length
iso = -0.03;
% Resolution grid → resolution of extracted offset surface
side = 60;
% Amount of smoothing to apply to distance field
sigma = 1.4;
bbd = norm(max(V)-min(V));
% Pad generously for positive iso values 
[BC,side,r] = voxel_grid([V;max(V)+iso*1;min(V)-iso*1],side);
D = signed_distance(BC,V,F);
D = reshape(D,side([2 1 3]));
% Smooth signed distance field
D = imfilter(D,fspecial('gaussian',9,sigma),'replicate');
BC3 = reshape(BC,[side([2 1 3]) 3]);
% Use matlab's built-in marching cubes iso-surface mesher (works most of the time)
surf = isosurface(BC3(:,:,:,1),BC3(:,:,:,2),BC3(:,:,:,3),D,iso*bbd);
SV = surf.vertices;
SF = surf.faces;

Here’s a blue bunny with a positive offset surface, an orange “cage”:

bunny with positive offset surface

Here’s a blue bunny with a negative offset surface. This is useful for hollowing out objects to 3d print:

bunny with positive offset surface

Because the iso-surface extraction will over tesselate low curvature patches of the output surface, it would make a lot of sense to remesh/decimate this mesh.

(to create these fancy renderings:)

clf;
hold on;
t  = tsurf(F,V,'EdgeColor','none',fsoft,  'FaceVertexCData',repmat(blue,size(V,1),1),'FaceAlpha',1+(iso<0)*(0.35-1),fphong);
ts = tsurf(SF,SV,'EdgeAlpha',0.2+(iso<0)*(0-0.2),fsoft,'FaceVertexCData',repmat(orange,size(SV,1),1),fphong,'FaceAlpha',1+(iso>0)*(0.2-1));
apply_ambient_occlusion(ts);
hold off;
axis equal;
view(-20,20)
camlight;
t.SpecularStrength = 0.04;
l = light('Position',[5 -5 10],'Style','local');
add_shadow(t,l,'Color',0.8*[1 1 1],'Fade','local','Ground',[0 0 -1 min([V(:,3);SV(:,3)])]);
set(gca,'pos',[0 0 1 1])
set(gca,'Visible','off');
set(gcf,'Color','w');
drawnow;

Quad meshing in matlab using gptoolbox

Friday, February 17th, 2017

This morning I played around with a very simple way of quad meshing a triangle mesh. This is more or less the “ideal pipeline” for quad meshing via parameterization. It will likely only work for very simple meshes. I’m pleasantly surprised that it works at all.

% Load a mesh with boundary and no topological (donut) holes
[V,F] = load_mesh('~/Dropbox/models/beetle.off');
V = V*axisangle2matrix([1 0 0],-pi/2);
% Construct the LSCM matrix
[~,Q] = lscm(V,F,[],[]);
M = massmatrix(V,F);
% Use Fiedler vector as parameterization
[EV,ED] = eigs(Q,repdiag(M,2),5,'sm');
U = reshape(EV(:,end-3),[],2);
% try to find "canonical" rotation
[sU,sS,sV] = svd(U'*U);
U = U*sU;
% Lay down a grid in 2D
h = ((max(U(:,1))-min(U(:,1)))/(160-2));
[X,Y] = meshgrid(min(U(:,1))-h:h:max(U(:,1))+h,min(U(:,2))-h:h:max(U(:,2))+h);
% Find all quads that are strictly inside the parameterized mesh
[QQ,QU] = surf2patch(X,Y,0*Y);QU = QU(:,1:2);
[sqrD,I,C] = signed_distance([QU QU(:,1)*0],[U U(:,1)*0],F);
QQ = QQ(all(abs(sqrD(QQ))<1e-10,2),:);
[QU,I] = remove_unreferenced(QU,QQ);
QQ = fliplr(I(QQ));
% Snap boundary vertices of quad mesh to boundary of mesh, and
% solve little Dirichlet problem to diffuse the displacements smoothly
QF = [QQ(:,[1 2 3]);QQ(:,[1 3 4])];
L = cotmatrix(QU,QF);
b = unique(outline(QQ));
[~,~,QU(b,:)] = point_mesh_squared_distance(QU(b,:),U,outline(F));
QU = min_quad_with_fixed(-L,[],b,QU(b,:));
% Locate each quad mesh vertex in the parameterized mesh
[sqrD,I,C] = signed_distance([QU QU(:,1)*0],[U U(:,1)*0],F);
B = barycentric_coordinates(C(:,1:2),U(F(I,1),:),U(F(I,2),:),U(F(I,3),:));
% Interpolate 3D positions at quad mesh vertex locations
QV = V(F(I,1),:).*B(:,1) + V(F(I,2),:).*B(:,2) + V(F(I,3),:).*B(:,3);

To create the image above use:

clf;
hold on;
tq = tsurf(QQ,QV,'FaceVertexCData',repmat(blue,size(QV,1),1),'EdgeColor','none',fsoft,fphong);
tt = tsurf(F,V-[0.8 0 0],'FaceVertexCData',repmat(1-0.8*(1-blue),size(V,1),1),'EdgeAlpha',0.5,fsoft,fphong);
axis equal;
camproj('persp');
view(-38,17);
drawn;
camlight;
plot_edges(QV,[QQ(:,2:3);QQ(:,[4 1])],'Color',orange);
plot_edges(QV,[QQ(:,1:2);QQ(:,3:4)],'w');hold off;
l = light('Position',[-10 0 10],'Style','infinite');
sq = add_shadow(tq,l,'Fade','infinite','Color',[0.8 0.8 0.8]);
st = add_shadow(tt,l,'Fade','infinite','Color',[0.8 0.8 0.8]);
set(gca,'Visible','off');
set(gcf,'Color','w');

There’s a lot to improve about the quad mesh: the boundaries are not aligned with the meshing direction, the meshing directions are not always aligned with principle curvature, etc. Not bad for 30 mins of coding using tools that are just lying around.

How does matlab’s hypot function work?

Saturday, January 7th, 2017

The hypot function in matlab purports to be more numerically stable at computing the hypotenuse of a (right-)triangle in 2D. The example the help hypot gives is:

a = 3*[1e300 1e-300];
b = 4*[1e300 1e-300];
c1 = sqrt(a.^2 + b.^2)
c2 = hypot(a,b)

where you see as output:

c1 =

   Inf     0

and

c2 =

       5e+300       5e-300

this is a compiled built-in function so you can’t just open hypot to find out what its doing. It might just be pre-dividing by the maximum absolute value. There’s probably a better reference for this, but I found it in: “Vector length and normalization difficulties” by Mike Day.

Continuing this example:

m = max(abs([a;b]))
c3 = m.*sqrt((a./m).^2 + (b./m).^2)

produces

c3 =

       5e+300       5e-300

While matlab’s hypot only accepts two inputs, pre-dividing by the maximum obviously extends to any dimension.

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.

Plot a bunch of edges in matlab with per-edge color data

Monday, December 19th, 2016

Suppose you have a list of vertices V (#V by dim) and a list of edges E (#E by 2) indexing V, then one hacky to draw all of these edges with per-vertex color data DV (#V by 1) is to use trisurf:

trisurf(E(:,[1 1 2]),V(:,1),V(:,2),V(:,3),'CData',DV,'EdgeColor','interp');

However, if you try to use this trick to draw per-edge color data DE (#E by 1), you’ll find that 'EdgeColor','flat' does not work. Instead you can explode you compact “mesh” of edges into individual edges and repeat you edge color data for each edge point:

trisurf(reshape([1:size(E,1),1:size(E,1)*2],[],3),V(E(:),1),V(E(:),2),V(E(:),3),'CData',repmat(DE,2,1),'EdgeColor','interp');

Signed polygon area in matlab

Sunday, October 30th, 2016

Suppose you have a polygon’s 2D corners stored in P, then the signed area is given by:

signed_polyarea = @(P) 0.5*(P(1:end,1)'*P([2:end 1],2)-P(1:end,2)'*P([2:end 1],1));