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

produces

                 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<<"Compare:"<<std::endl;
  std::cout<<(X.size()         ?"done":"")<<std::endl;
  std::cout<<(simple_size(X)          ?"done":"")<<std::endl;
  std::cout<<(templated_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
    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

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/MATLAB_R2017a.app/toolbox/matlab/imagesci/private/readgif.m 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]));

Update:

Use this instead:

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

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.

Pause (and then resume) Battery-Guzzling programs

June 7th, 2017

My laptop battery dies quickly these days. Certain apps (cough, cough, Slack), have very high idle CPU-usage. You can pause these programs with

killall -STOP Slack

And later you can resume the application with

killall -CONT Slack

Convincing LatexIt and Illustrator to use the new SIGGRAPH fonts

May 20th, 2017

The SIGGRAPH Latex style changed to the Libertine font. Here’re the steps to convince Latexit to use the new stylesheet and then to convince Illustrator to use the libertine font for drag and drop math.

mkdir ~/Library/texmf/tex/latex/local/acmart.cls/
cp ~/Dropbox/boundary/Paper/acmart.cls ~/Library/texmf/tex/latex/local/acmart.cls

In Latexit, open up Preferences, add a new SIGGRAPH “Template” containing:

\documentclass[sigconf, review]{acmart}
\pagenumbering{gobble}

If you try to drag and drop these into illustrator you’ll see that illustrator has replaced the nice math font with Myriad or something silly.

Download the TTF libertine font pack

Drag this into FontBook.app

cp /usr/local/texlive/2015/texmf-dist/fonts/type1/public/libertine/*.pfb ~/Library/Application\ Support/Adobe/Fonts/

Update: I also had to issue:

cp /usr/local/texlive/2015/texmf-dist/fonts/type1/public/txfonts/*.pfb ~/Library/Application\ Support/Adobe/Fonts/
cp /usr/local/texlive/2015/texmf-dist/fonts/type1/public//newtx/*.pfb  ~/Library/Application\ Support/Adobe/Fonts/

If you see boxes with X’s replacing symbols after dragging and dropping from LaTeXit, then drag into Finder instead (to create a .pdf file), then open this directly and Illustrator will give a warning and tell you which font it’s (still) missing.

Restart Adobe Illustrator

Mex wrapper for graph segmentation

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