Installing and running fmmlib from MATLAB on Mac OS X

December 26th, 2017

No better way to finish 2017 than trying to compile some old fortran code. This time it’s the fast multipole method library from nyu, fmmlib.

Here’s how I did it on a MacBook Pro running OS X 10.12 with MATLAB r2018a pre-release.

First install gfortran using homebrew via:

brew install gcc

Now we need to convince to use find and use gfortran.

Download gfortran.xml


Download gfortran.xml

In MATLAB issue:

copyfile('~/Downloads/xcode7_mexopts/gfortran.xml',fullfile( matlabroot, 'bin', 'maci64', 'mexopts' ))

On the command line, install gfortrani8 in the same directory as gfortran. I did this via:

cp ~/Downloads/gfortrani8 /usr/local/Cellar/gcc/7.2.0/bin/gfortrani8

Now, the following should succeed in MATLAB to setup the gfortran compiler:

mex -setup -v FORTRAN

For posterity: If you get “sdk” related issues, you may need to copy all the 10.13 lines in gfortran.xml and create analogous lines for 10.14 etc.

It seems like the '-compatibleArrayDims' flag should be used with fmmlib, but honestly I’m not sure. Something to revisit if it crashes on giant input I guess.

Finally, I had trouble locating the libgfortran.a library and convincing MATLAB’s Mex to use it. Therefore, I added:

[~,libgfortran] = system('gfortran -print-file-name=libgfortran.dylib');
libgfortran = libgfortran(1:end-1);

to my Mex script so that I can issue lines like this:


This is enough to compile all of the fortran source files one by one and then link them during compilation and linking of the fmm3d_r2012a.c file.

I’ve added a matlab/compile.m script to my GitHub mirror/fork of the fmmlib library. So if you do the steps above you should be to just issue:

cd ~/Downloads
git clone 
cd ~/Downloads/fmmlib3d-1.2/matlab/


Eigen Gotcha: sizeof Eigen::Index depends on system

December 19th, 2017

Part of the magic behind libigl’s header-only / static-library duality is making use of explicit template instantiations. Usually I can just literally copy the compiler “missing symbol” error lines, stick template in front of them and a semicolon ; at the end to add a symbol to the static library.

Unfortunately this falls apart for any templates that are defined via Eigen::Index (e.g., via Eigen::MatrixXi::Index).

My compiler (clang on Mac OS X) translates these to long in the output lines. And on Mac OSX with clang a long is 64 bits long.

Unfortunately, long is only 32 bits long on Windows and there Eigen::Index is translated to __int64 (which has 64 bits). So then I still get a missing symbol.

This is probably more of a clang compiler gotcha. I’m not sure Eigen could do anything to help me here.

Update: I take that back. Eigen could define EIGEN_DEFAULT_DENSE_INDEX_TYPE to a fixed size value (e.g., int64_t) rather than the system dependent std::ptrdiff_t, but there’s an obvious advantage to using the best type on this system rather than a fixed size that might be too big (or too small).

Accessing Free Guangzhou Airport Wi-Fi on MacBook Pro

December 2nd, 2017

Guangzhou airport has a crazy way of providing free wifi. If you have a Chinese phone you can receive a login account via text, but if you don’t have a Chinese number, then you have to go to a physical ticket machine, scan your passport, and quickly jot down (or photograph) the code that flashes on the screen for 15 seconds.

Once you have this code, you have to click “Ticket Cert” on the login page and enter the info.

This worked fine for me on my iPhone.

But when I tried to follow the same steps on my laptop, the “Ticket Cert” option did not appear on my browser login prompt. I saw a slightly different page that only had the SMS option (and a bunch of half-loaded CSS and javascript errors).

I tried many things including spoofing the UserAgent on my browser. Nothing seemed to work.

Finally, I changed the MAC address on my laptop to match my already-online iphone’s

ifconfig en0 ether [iphone's "Wi-Fi Address"]

This worked.

PS: So far, I have not managed to get any sort of vpn, proxy or ssh tunnelling to work.

Convincing maple to solve an ODE with Neumann conditions at a symbolic valued location

November 17th, 2017

I can use maple to solve a 1D second-order ODE with Dirichlet boundary conditions at symbolic-valued locations:

# Z'' = 0, Z(a)=0, Z(b) = 1
dsolve({diff(Z(r),r,r) = 0,Z(a)=0,Z(b)=1}); 

This correctly returns

                  r       a
       Z(r) = - ----- + -----
                a - b   a - b

I can also easily convince maple to solve this ODE with some Neumann (normal derivative) boundary conditions at at fixed-value, numeric location:

# Z'' = 0, Z(a) = 1, Z'(0) = 0
dsolve({diff(Z(r),r,r) = 0,Z(a)=1,eval(diff(Z(r),r),r=0)=0});


                 Z(r) = 1

But if I try naively to use a Neumann condition at a symbolic value location

# Z'' = 0, Z(a) = 1, Z'(b) = 0
dsolve({diff(Z(r),r,r) = 0,Z(a)=1,eval(diff(Z(r),r),r=b)=0});

then I get an error:

Error, (in dsolve) found differentiated functions with same name but depending on different arguments in the given DE system: {Z(b), Z(r)}

After a long hunt, I found the solution. dsolve takes an optional second argument that can tell it what the dependent variable actually is. So the correct call is:

# Z'' = 0, Z(a) = 1, Z'(b) = 0
dsolve({diff(Z(r),r,r) = 0,Z(a)=1,eval(diff(Z(r),r),r=b)=0});

and this gives the correct answer

                 Z(r) = 1

MATLAB gotcha inverting a (sparse) diagonal matrix

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

Eigen performance gotcha calling non-templated function from templated one

July 25th, 2017

I just spent a while tracking down a rather surprising performance bug in my code.

Here’s a minimal example:

#include <Eigen/Dense>
#include <iostream>

int simple_size(const Eigen::MatrixXi & Z)
  return Z.size();

template <typename T> int templated_size(const Eigen::MatrixBase<T> & Y)
  return simple_size(Y);

int main(int argc, const char * argv[])
  const int s = 40000;
  Eigen::MatrixXi X = Eigen::MatrixXi::Zero(40000,40000);
  std::cout<<(X.size()         ?"done":"")<<std::endl;
  std::cout<<(simple_size(X)          ?"done":"")<<std::endl;


Running this, it will show that the last call to templated_size is taking way too long. Inspection will show that a copy of Y is being created to create a Eigen::MatrixXi & reference.

Now, clearly it’s poor design to call a function expecting a Eigen::MatrixXi & reference with a generic templated type Eigen::MatrixBase<T> &, but unfortunately this happens quite often with legacy libigl functions. My expectation was that since T is Eigen::MatrixXi in this case a simple reference would be passed.

It’s worth noting that const is actually creating/hiding the problem. Because simple_size takes a const reference, the compiler is happy to construct a Eigen::MatrixXi on the fly to create a valid reference. Without the consts the compiler stops at an error.

Paper-worthy rendering in MATLAB

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
  case 5
  case 6
    l = light('Position',[0.2 -0.2 1]);
  case 7
  case 8
  case 9
    s = add_shadow(t,l,'Color',bg_color*0.8,'BackgroundColor',bg_color,'Fade','infinite');
  case 10

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


Inflate Wire Mesh in libigl C++ or gptoolbox MATLAB

July 12th, 2017

For a visualization and 3D printing, it’s often useful to “inflate” a edge-network into a thickened surface mesh. One method to do this is described “Sculptural Forms from Hyperbolic Tessellations” by George C Hart. This method works by adding rotated polygons at the ends of each edge offset a bit from the vertices. Then for each vertex the convex hull of incident edges’ polygons is computed and unioned with the convex hull of the polygons at either end of each edge. Hart writes that polygons shared by “edge hulls” and “vertex hulls” can simply be discarded. This is unfortunately not true, in general. It’s not super easier to categorize which faces can be discarded (even in general position) since the answer depends on the thickness, the number of sides of the polygons, their rotations, their offsets, and the angle between neighbouring edges. Fortunately, libigl is very good at conducting unions. We can just conduct the union explicitly and exactly using libigl.

I’ve written a new function for libigl igl::wire_mesh that takes in a wire network and spits out a solid (i.e., closed, watertight, manifold) mesh of a the inflated surface.

I’ve also wrapped this up in a Matlab Mex function in gptooolbox wire_mesh.

Read animated gif and convert to rgb frames

June 14th, 2017

Matlab’s built in imread (as of 2017) doesn’t load animated gifs correctly. You can fix this by changing the line:

map = info.ColorTable;

in /Applications/ with

map = reshape([info(:).ColorTable],[],3,n);

For a single frame indexed image you can use ind2rgb to convert to an rgb image. To do this on entire animated gif. You can use an arrayfun for-loop hack:

[X,M] = imread('input.gif');
Y = cell2mat(permute(arrayfun(@(C) ind2rgb(X(:,:,:,C),M(:,:,C)),1:size(X,4),'UniformOutput',false),[1 4 3 2]));


Use this instead:

map = zeros(0,3,n);
for j = 1:n
  map(1:size(info(j).ColorTable,1),:,j) = info(j).ColorTable;

Project page for “Generalized Matryoshka: Computational Design of Nesting Objects”

June 14th, 2017

This April I had fun working on a little project that’s been tickling my mind for a while. Can you make any shape into a Matryoshka doll?

I’ll be presenting my paper on this at SGP 2017. It’s entitled
Generalized Matryoshka: Computational Design of Nesting Objects.