Archive for March, 2012

Load workspace in background “thread” in MATLAB

Tuesday, March 27th, 2012

Here’s a proof of concept that you can run in your matlab IDE that loads a workspace from a file in a background, worker thread. I imagine this is useful if you want to keep your Matlab session open, but have a 3rd party program be able to talk to matlab, (not just be called from matlab (mex) or call matlab functions (matlab engine).

First create some dummy workspaces.


X = 1;
save('workspace.mat','X');
X = 2;
save('workspace2.mat','X');

Now run this code to start the background “thread”:


% load the global workspace from 'workspace.mat' every 1 second, print value of X
t = timer('TimerFcn','load(''workspace.mat'');disp(X)','Period',1,'ExecutionMode','fixedDelay');
start(t);

You should see that matlab is spitting out:


X =
     1

But in the meantime you could issue matlab calls like:


Y = 'banana'

Now pretend you’re a third party program that wants to send a new value of X to the matlab session. Do this by copying our dummy workspace2.mat onto the workspace.mat that’s being polled:


!cp workspace2.mat workspace.mat

Now we see that it has worked as matlab is spitting out:


X = 
    2

Finally stop the polling and clean up the timer:


stop(t);
delete(t);

AJX MassMailer source

Monday, March 26th, 2012

ajx mass mailer splash
I decided to release the source code of an old project of mine. It’s a simple Java app that handles batch/scripted emails to large numbers of recipients using your gmail/yahoo/etc. account.
The idea is you load/fill in a table where each row is a recipient and each column is a piece of corresponding data about that recipient. Then you can write a template email that gets sent to each recipient, pulling the relevant data from each row of the table.
ajx mass mailer screen shot
Project page
Source as zip

STL unique != MATLAB unique

Monday, March 26th, 2012

Got burned on the STL unique algorithm not actually implementing what I inferred its name to mean. In MATLAB unique is finds unique elements in an unordered set.

unique([1 2 3 3 2 1]) --> [1 2 3]

In STL, unique reorders an ordered set so that repeated equal entries are replaced with a single entry:

unique([1 2 3 3 2 1]) --> [1 2 3 2 1]

To get the behavior of the matlab unique you’ll need to sort first.

#!/bin/bash
/*/../bin/ls > /dev/null
# BEGIN BASH SCRIPT
printf "//" | cat - $0 | g++ -I/opt/local/include/eigen3 -I$IGL_LIB/include -o .main -x c++ - && ./.main
rm -f .main
# END BASH SCRIPT
exit
*/
#include <iostream>
#include <vector>
#include <algorithm>
#include <iterator>
int main(int argc,char * argv[])
{
  using namespace std;
  vector<int> V;
  V.push_back(1);
  V.push_back(2);
  V.push_back(3);
  V.push_back(3);
  V.push_back(2);
  V.push_back(1);
  vector<int> SV = V;

  cout<<"V=";
  copy(V.begin(), V.end(), ostream_iterator<int>(cout, " "));
  cout<<endl;

  V.erase(unique(V.begin(), V.end()), V.end());
  cout<<"unique(V)=";
  copy(V.begin(), V.end(), ostream_iterator<int>(cout, " "));
  cout<<endl;

  sort(SV.begin(),SV.end());
  SV.erase(unique(SV.begin(), SV.end()), SV.end());
  cout<<"unique(sort(V))=";
  copy(SV.begin(), SV.end(), ostream_iterator<int>(cout, " "));
  cout<<endl;

  return 0;
}

which produces:


V=1 2 3 3 2 1 
unique(V)=1 2 3 2 1 
unique(sort(V))=1 2 3 

Eigen gotcha: aliased copy of SparseMatrix results in empty matrix

Saturday, March 24th, 2012

Found out the hard way today that you can’t copy an Eigen sparse matrix onto itself (i.e. you cannot issue A=B, where A and B are references to the same sparse matrix. The result is that both A and B end up empty.

Here’s a self-compiling snippet that shows what happens:

#!/bin/bash
/*/../bin/ls >/dev/null
# BEGIN BASH SCRIPT
printf "//" | cat - $0 | g++ -g -o .main -x c++ - -I/opt/local/include/eigen3/ -I/opt/local/include/eigen3/upsupported  && ./.main
rm -f .main
# END BASH SCRIPT
exit
*/
#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#include <Eigen/Sparse>

int main(int argc, char** argv)
{
  using namespace std;
  using namespace Eigen;
  SparseMatrix<double> I(2,2);
  I.insert(0,0) = 1;
  I.insert(1,1) = 1;
  I.finalize();
  cout<<"I.nonZeros(): "<<I.nonZeros()<<endl;
  SparseMatrix<double> Icopy;
  cout<<"Icopy=I"<<endl;
  Icopy = I;
  cout<<"I.nonZeros(): "<<I.nonZeros()<<endl;
  cout<<"Icopy.nonZeros(): "<<Icopy.nonZeros()<<endl;
  cout<<"I=I"<<endl;
  I = I;
  cout<<"I.nonZeros(): "<<I.nonZeros()<<endl;
}

Which produces:


I.nonZeros(): 2
Icopy=I
I.nonZeros(): 2
Icopy.nonZeros(): 2
I=I
I.nonZeros(): 0

the shallow heart – gary crane’s biography of gary crane

Tuesday, March 20th, 2012

50¢ compact disc from 2005
album cover for shallow heart gary crane's biography of gary crane
tiny liner notes for shallow heart gary crane's biography of gary crane

mp3
mp3
mp3
Bonus:
mp3
Whole Album + bonus track

mnt bnncrt

Monday, March 19th, 2012

mnt bnncrt

http://alecjacobson.com/art/digital/
http://alecjacobson.com/art/

Isoline plots on triangle meshes in Matlab

Friday, March 9th, 2012

I recently posted how to plot nice looking iso intervals of scalar fields on triangle meshes using matlab. Along with the intervals I gave an ad hoc way of also showing isolines. This works ok if your data has fairly uniform gradient but they’re not really isolines.

Here’s a little matlab function to output real isolines. Save it in a file called isoline.m:


function [LS,LD,I] = isolines(V,F,S,iso)
  % ISOLINES compute a list of isolines for a scalar field S defined on the
  % mesh (V,F)
  %
  % [LS,LD,I] = isolines(V,F,S,iso)
  %
  % Inputs:
  %   V  #V by dim list of vertex positions
  %   F  #F by 3 list of triangle indices
  %   S  #V list of scalar values defined on V
  %   iso #iso list of iso values
  % Outputs:
  %   LS  #L by dim list of isoline start positions
  %   LD  #L by dim list of isoline destination positions
  %   I  #L list of indices into iso revealing corresponding iso value
  %
  % alternative: tracing isolines should result in smaller plot data (every
  % vertex only appears once
  %

  % make sure iso is a ROW vector
  iso = iso(:)';
  % number of isolines
  niso = numel(iso);

  % number of domain positions
  n = size(V,1);
  % number of dimensions
  dim = size(V,2);

  % number of faces
  m = size(F,1);

  % Rename for convenience
  S1 = S(F(:,1),:);
  S2 = S(F(:,2),:);
  S3 = S(F(:,3),:);

  % t12(i,j) reveals parameter t where isovalue j falls on the line from
  % corner 1 to corner 2 of face i
  t12 = bsxfun(@rdivide,bsxfun(@minus,iso,S1),S2-S1);
  t23 = bsxfun(@rdivide,bsxfun(@minus,iso,S2),S3-S2);
  t31 = bsxfun(@rdivide,bsxfun(@minus,iso,S3),S1-S3);

  % replace values outside [0,1] with NaNs
  t12( (t12<-eps)|(t12>(1+eps)) ) = NaN;
  t23( (t23<-eps)|(t23>(1+eps)) ) = NaN;
  t31( (t31<-eps)|(t31>(1+eps)) ) = NaN;

  % masks for line "parallel" to 12
  l12 = ~isnan(t23) & ~isnan(t31);
  l23 = ~isnan(t31) & ~isnan(t12);
  l31 = ~isnan(t12) & ~isnan(t23);

  % find non-zeros (lines) from t23 to t31
  [F12,I12,~] = find(l12);
  [F23,I23,~] = find(l23);
  [F31,I31,~] = find(l31);
  % indices directly into t23 and t31 corresponding to F12 and I12
  %ti12 = sub2ind(size(l12),F12,I12);
  %ti23 = sub2ind(size(l23),F23,I23);
  %ti31 = sub2ind(size(l31),F31,I31);
  % faster sub2ind
  ti12 = F12+(I12-1)*size(l12,1);
  ti23 = F23+(I23-1)*size(l23,1);
  ti31 = F31+(I31-1)*size(l31,1);

  % compute actual position values
  LS = [ ...
    ... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t23(ti12), V(F(F12,2),:)) + ...
    bsxfun(@times,  t23(ti12), V(F(F12,3),:)) ...
    ; ... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t31(ti23), V(F(F23,3),:)) + ...
    bsxfun(@times,  t31(ti23), V(F(F23,1),:)) ...
    ;... % average of vertex positions between 2 and 3
    bsxfun(@times,1-t12(ti31), V(F(F31,1),:)) + ...
    bsxfun(@times,  t12(ti31), V(F(F31,2),:)) ...
    ];
  LD = [...
    ... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t31(ti12), V(F(F12,3),:)) + ...
    bsxfun(@times,  t31(ti12), V(F(F12,1),:)) ...
    ;... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t12(ti23), V(F(F23,1),:)) + ...
    bsxfun(@times,  t12(ti23), V(F(F23,2),:)) ...
    ;... % hcat with average of vertex positions between 3 and 1
    bsxfun(@times,1-t23(ti31), V(F(F31,2),:)) + ...
    bsxfun(@times,  t23(ti31), V(F(F31,3),:)) ...
    ];

  % I is just concatenation of each I12
  I = [I12;I23;I31];

end

Then if you have a mesh in (V,F) and so scalar data in S. You can display isolines using:


colormap(jet(nin));
trisurf(F,V(:,1),V(:,2),V(:,3),'CData',S,'FaceColor','interp','FaceLighting','phong','EdgeColor','none');
axis equal;
[LS,LD,I] = isolines(V,F,S,linspace(min(S),max(S),nin+1));
hold on;
plot([LS(:,1) LD(:,1)]',[LS(:,2) LD(:,2)]','k','LineWidth',3);
hold off
set(gcf, 'Color', [1,1,1]);
set(gca, 'visible', 'off');
%O = outline(F);
%hold on;
%plot([V(O(:,1),1) V(O(:,2),1)]',[V(O(:,1),2) V(O(:,2),2)]','-k','LineWidth',2);
%hold off

Uncomment the last lines to display the outline of the mesh. This produces:
woody matlab isolines
If you use the my anti-alias script you can get pretty good results, almost good enough for a camera-ready paper:
woody matlab isolines anti-aliased

“Free” two-line font

Thursday, March 8th, 2012

Writing SIGGRAPH rebuttals, we’re all paranoid about word count (we get 1000 words of ASCII text to refute the reviews). The usual word count tricks are often employed: contraction “do not” to “don’t”, dropping articles, dropping “el al.”s from citations. I’ve even heard of even riskier tricks like combining “Reviewer #01” into “Rev01”. Most of these tricks come at the expense of clarity and professionalism, but it got me thinking about the word count game.

The word count algorithm employed by SIGGRAPH’s SIS system follows a view rules revealed through their source. The gist is that “non-alphanumerics” get replaced by spaces, then the count of words are just the number of tokens separated by whitespace. In reg-ex form non-alphanumerics are defined here to be:


[^A-Za-z0-9']

Note:The only interesting thing here is that the apostrophe is OK meaning, where as “light-blue” counts as two words “don’t” counts as one.
Note:There’s another slight subtlety that if an apostrophe occurs alone it is OK. This just means ” ‘ ” is 0 words and ” -‘- ” is also 0 words.

But if non-alphanumerics show up alone, that is, never neighboring an alphanumeric character, then you get 0 words.

This leads to a ridiculous “trick” for gaming the word count system by employing a “two-line” font composed entirely of non-alphanumeric words. Here’s quick prototype:


/'`  /_\  |\ |   \ / /'\ | |   |_) /_  /_\  |'\
\_. /   \ | \|    |  \_/ |_|   | \ \_ /   \ |_/


''|'' |_| |  ('   [' /'\ |\ | ''|''   )
  |   | | |  _)   |  \_/ | \|   |     .

Note: The obvious danger (besides being embarrassed by actually submitting something written in this font) is that the reviewers may not see it in a correctly line-separated or monospaced, original font. In that case it just looks like garbage.

Get function handle to matlab’s quadprog, while shadowed by mosek

Monday, March 5th, 2012

Here’s a small function I use to retrieve a function handle to MATLAB’s quadprog which is shadowed when I have installed the mosek toolbox:


function [quadprog_matlab,matlab_path] = get_quadprog_matlab()
  % GET_QUADPROG_MATLAB returns function handle to matlab's builtin quadprog
  % (useful if, say, mosek has shadowed it)
  %
  % [quadprog_matlab,matlab_path] = get_quadprog_matlab();
  %
  % Outputs:
  %  quadprog_matlab  function handle to matlab's quadprog
  %  matlab_path  path to matlab's quadprog
  %
  % [quadprog_matlab,matlab_path] = get_quadprog_matlab();
  % help(matlab_path);
  % X = quadprog_matlab(H,f,A,b);

  % get list of all quadprogs
  paths = which('quadprog.m','-ALL');
  % remember current directory
  old_dir = pwd;
  matlab_path = paths(~cellfun(@isempty,regexp(paths,'MATLAB')));
  matlab_path = matlab_path{1};
  % move to that directory
  cd(regexprep(matlab_path,'quadprog.m',''));
  quadprog_matlab = @quadprog;
  cd(old_dir);
end

Recover Lagrange multiplier values for known values in quadratic energy minimization

Monday, March 5th, 2012

Often I’m minimizing energies of the form:


E = ½ xT yT Q x + xTyT b
              y

where Q is a quadratic coefficient matrix, b is a vector of linear coefficients, x are unknown, and y are known. Now normally if I just want to minimize this energy I treat the known values (y) as constants and take the derivative with respect to the unknowns (x) and set it to zero. Then I just move all the terms involving y to the right hand side.

However if I’d like to know how setting these known values affected the energy minimization then its useful to treat the known values not as constant but as constraints. So I have the same energy as above, where now both x and y are variables, but subject to the constraints that


y = y*

where y* are the known values of y.
I can enforce these via Lagrange multipliers. Then once I’ve found my solution the values of the Lagrange multipliers corresponding to each y will tell me the effect of that y‘s constraint on the energy minimization. The problem now becomes finding the saddle point (equilibrium) of the following Lagrangian:


Λ = ½ xT yT Q x + xTyT b + λT y - λT y*
              y

Or equivalently by repeating the λTy terms:


Λ = ½ xT yT λT Qxx Qxy 0 x + xTyT λT bx
              Qyx Qyy I y           by
              0   I   0 λ          -y*

Now, let’s start to take derivatives with respect to each of the sets of variables. We begin with λ. Setting ∂Λ/∂λT = 0 gives us:


∂Λ    0 I 0 x + -y*
--- =       y       = 0
∂λT         λ

which reduces simply to revealing our known values:


y = y*

Now look at setting ∂Λ/∂xT = 0 which gives us:


∂Λ    Qxx Qxy 0 x + bx
--- =           y       = 0
∂xT             λ

which after substituting our result y=y* reduces to the system of equations we are used to seeing when minimizing quadratic energies with fixed values:


Qxx x = -bx - Qxy y*

So we can invert Qxx and solve for x:


x* = Qxx-1 ( -bx - Qxy y* )

Now x and y are known, all that’s left is to solve for λ by taking the last set of derivatives. So we set ∂Λ/∂yT = 0 which gives us:


∂Λ    Qyx Qyy I x + by
--- =           y       = 0
∂yT             λ

Plugging in our known values and moving them to the right-hand side reveals the values of λ:


λ* = - Qyx x* - Qyy y* - by

The great thing about all of this is that you don’t have to do anything extra. If you’re already set up to immediately solve for x treating y as constant you’ve got everything you need to determine the values of λ without every actually implementing the equivalent Lagrangian.