http://alecjacobson.com/art/digital/
http://alecjacobson.com/art/
Archive for February, 2011
judith one and one half
Monday, February 28th, 2011tell her or the terror will tear her
Thursday, February 24th, 2011
“Element-wise” matrix vector multiplication in matlab
Wednesday, February 23rd, 2011or
Multiply each matrix in a list of matrices against corresponding vector in list of vectors
Here’s a code snippet to take a stacked list of m by n matrices, S, and multiply each m by n matrix by a corresponding n-length row vector in a list of row vectors, V. Resulting in a list of m-length row vectors:
First of all, make your vertical stack of matrices S. So that S(i:(i+m-1),:) = Mi and so forth. And be sure that your list of row vectors is s by n, in V, where s is the number of matrices and vectors.
Then you may issue:
I = repmat([1:(m*s)],1,n);
J = repmat(reshape(repmat(1:n:(s*n),m,1),1,s*m),1,n)+reshape(repmat(0:(n-1),s*m,1),1,s*m*n);
SV = reshape(sparse(I,J,S(:))*reshape(V',s*n,1),m,s)';
Often I find that I also have a bunch of matrices stored in a 3d-array, say, T. Where T(:,:,i) = Mi and so forth. I can easily convert T, an m by n by s 3d-array into S, an m*s by n array, using:
S = reshape(permute(T,[2,1,3]),[n,s*m])';
Update:
Here’s working code that uses the above:
s = 4;
m = 2;
n = 3;
% generate some bogus matrices
S = repmat(eye(m,n),s,1).*repmat(reshape(repmat(1:s,m,1),s*m,1),1,n);
% generate some bogus row vectors
V = reshape(1:s*n,n,s)';
% compute "element-wise" matrix-vector multiplication
I = repmat([1:(m*s)],1,n);
J = repmat(reshape(repmat(1:n:(s*n),m,1),1,s*m),1,n)+reshape(repmat(0:(n-1),s*m,1),1,s*m*n);
SV = reshape(sparse(I,J,S(:))*reshape(V',s*n,1),m,s)';
% view result
SV
% convert vertical stack into 3d array stack
T = permute(reshape(S,[m,s,n]),[1,3,2]);
% convert 3D array stack to vertical stack
TS = reshape(permute(T,[2,1,3]),[n,s*m])';
Update: Here’s the for loop version, which is much more straight forward, but starts to slow down in comparison for large s:
SV = zeros(s,m);
for ii = 0:(s-1)
SV(ii+1,:) = (S((m*ii+1):(m*ii+m),:)*(V(ii+1,:)'))';
end
List of installed applications on mac
Wednesday, February 23rd, 2011Seems like there are several hacky ways to achieve this. And it depends on what you mean by “installed”. If you just mean, list all application bundles (files ending with .app) then you can achieve it the following command in your terminal:
find / -name "*.app"
That command returns a list where each line is the full path to the application. If you just want the name of the application then use:
find / -name "*.app" | sed -e "s/\(.*\)\/\([^\/]*\).app/\2/g"
If you’re really a stickler and want to see all applications even ones the current user does not have permissions to see then use:
sudo find / -name "*.app" | sed -e "s/\(.*\)\/\([^\/]*\).app/\2/g"
On my computer this takes around ∞ minutes.
But maybe you don’t mean to find every application anywhere on the computer, but say just the ones in the /Applications folder. Then you can slightly change the above command to:
find /Applications -name "*.app" | sed -e "s/\(.*\)\/\([^\/]*\).app/\2/g"
On my computer this takes around 3 minutes.
Notice that this is still taking a considerable amount of time and its returning a bunch of weird applications. These are Utility applications or example applications buried deep in the sub folders of the /Applications folder. You can speed up the previous command by providing a maximum depth for the search:
find /Applications -name -maxdepth 1 "*.app" | sed -e "s/\(.*\)\/\([^\/]*\).app/\2/g"
This only finds .app files directly in the /Applications folder. So it doesn’t find the Microsoft Office or Adobe CS5 apps I have because those are grouped into a folder. But it only takes 0.050 seconds.
I settled on using -maxdepth 2:
find /Applications -name -maxdepth 2 "*.app" | sed -e "s/\(.*\)\/\([^\/]*\).app/\2/g"
This finds all of my meaningful applications and in only 0.060 seconds. Higher maxdepths found a few more apps but they were weirdo utility/example apps and made the command take over a second.
If you need this list as an applescript list just use the following:
set list_of_installed_apps to paragraphs of (do shell script "find /Applications/ -name \"*.app\" -maxdepth 2| sed -e \"s/\\(.*\\)\\/\\([^\\/]*\\).app/\\2/g\"")
Rotate a point around another point
Tuesday, February 22nd, 2011Composing rotations and translations is one of the most important operations in computer graphics. One useful application is the ability to compose rotations and translations to rotate a point around another point.
Here’s how I learned to do this. I find this method the most intuitive.
- Translate so that the other point is at the origin
- Rotate about the origin by however much you want
- Translate back so that the other point is back to its original position
This works out nicely with matrix math since rotations around the origin are easily stored as 2×2 matrices:
/ cos θ -sin θ\
Rotation around origin by θ: | |
\ sin θ cos θ/
If we call that matrix, R, then we can write the whole operation that rotates a point, a, around another point, b,, as:
R*(a-b) + b. Be careful to note the order of operations: (a-b) corresponds to step 1, then left multiply with R to step 2, and finally adding back b is step 3.
One convenient fact is that when we look at the transformation as a whole where T(a): a → R*(a-b) + b. Because T is just the composition of rotations and translations it can be decomposed into a single rotation followed by a single translation. Namely, T(a): a → R*(a-b) + b may be re-written as T(a): a → S(a) + t, for some rotation matrix S and some translation vector t. To see this just use the distributive law: R*(a-b) + b = R*a – R*b + b, then S = R and t = R*b + b.
So that gives a new derivation of the transformation that rotates a point, a, around another point, b,: R*a – R*b + b. Building an intuition as to why this works is a little tricky. We have just seen how it can be derived using linear algebra but actually seeing this version as each step is elusive. Turns out it helps to make a subtle change. Reverse the last two terms, so that you have: R*a + (b – R*b). Now intuitively we can follow the order of operations and build an intuition:
- Rotate first the point about the origin (since the other point is not the origin we’ve “missed” where we should have been rotated by a certain error amount)
- Rotate the other point about the origin by the same amout
- Translate by the difference between the original other point and the rotated other point (this is the error amount, because we know that rotating other pointabout itself shouldn’t change its position)
Update: Here’s an image comparing these two compositions:
Map grayscale to color using colormap
Saturday, February 19th, 2011In matlab you can view a grayscale image with:
imshow(im)
Which for my image im shows:
And you can also view this grayscale image using pseudocolors from a given colormap with something like:
imshow(im,'Colormap',jet(255))
Which shows:
But it’s not obvious how to use the colormap to actually retrieve the RGB values we see in the plot. Here’s a simple way to convert a grayscale image to a red, green, blue color image using a given colormap:
rgb = ind2rgb(gray2ind(im,255),jet(255));
Replace the 255 with the number of colors in your grayscale image. If you don’t know the number of colors in your grayscale image you can easily find out with:
n = size(unique(reshape(im,size(im,1)*size(im,2),size(im,3))),1);
It’s a little overly complicated to handle if im is already a RGB image.
If you don’t mind if the rgb image comes out as a uint8 rather than double you can use the following which is an order of magnitude faster:
rgb = label2rgb(gray2ind(im,255),jet(255));
Then with your colormaped image stored in rgb you can do anything you normally would with a rgb color image, like view it:
imshow(rgb);
which shows the same as above:
Possible function names include real2rgb, gray2rgb.
Read local file and display contents using javascript
Friday, February 18th, 2011The upcoming File API for web applications means that HTML 5 pages will be able to ask user for a local file (a file on the user’s computer) and use the data in that file directly. It’s amazing that this hasn’t really existed before. Web applications either used http requests to send the file to the server then have the server send it back, or use super hacky solutions that typically involved using Java or Flash applets to piggy-back into the local file.
Unfortunately, this is not fully supported by browsers. It seems Firefox and Chrom(e|ium) support this File API already, but Safari and Webkit don’t yet. IE may never.
Check it out:
Progress:
File contents:
Here's the source as a full HTML page:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN"
"http://www.w3.org/TR/html4/strict.dtd">
<html>
<head>
<meta http-equiv='Content-type' content='text/html;charset=UTF-8' >
<script>
function startRead()
{
// obtain input element through DOM
var file = document.getElementById('file').files[0];
if(file)
{
getAsText(file);
}
}
function getAsText(readFile)
{
var reader;
try
{
reader = new FileReader();
}catch(e)
{
document.getElementById('output').innerHTML =
"Error: seems File API is not supported on your browser";
return;
}
// Read file into memory as UTF-8
reader.readAsText(readFile, "UTF-8");
// Handle progress, success, and errors
reader.onprogress = updateProgress;
reader.onload = loaded;
reader.onerror = errorHandler;
}
function updateProgress(evt)
{
if (evt.lengthComputable)
{
// evt.loaded and evt.total are ProgressEvent properties
var loaded = (evt.loaded / evt.total);
if (loaded < 1)
{
// Increase the prog bar length
// style.width = (loaded * 200) + "px";
document.getElementById("bar").style.width = (loaded*100) + "%";
}
}
}
function loaded(evt)
{
// Obtain the read file data
var fileString = evt.target.result;
document.getElementById('output').innerHTML = fileString;
document.getElementById("bar").style.width = 100 + "%";
}
function errorHandler(evt)
{
if(evt.target.error.code == evt.target.error.NOT_READABLE_ERR)
{
// The file could not be read
document.getElementById('output').innerHTML = "Error reading file..."
}
}
</script>
</head>
<body>
<input id="file" type="file" multiple onchange="startRead()">
<h3>Progress:</h3>
<div style="width:100%;height:20px;border:1px solid black;">
<div id="bar" style="background-color:#45F;width:0px;height:20px;"></div>
</div>
<h3>File contents:</h3>
<pre>
<code id="output">
</code>
</pre>
</body>
</html>
Note: This WILL NOT WORK if you load your html page as a local file. Meaning the address starts with file:/// instead of http://, for some reason the page must be served to the browser from a server.
Determine more recent of two files in bash
Friday, February 18th, 2011Here’s a simple bash script that determines the more recent of two files”:
#!/bin/bash
# Usage:
# morerecent file1 file2
if [[ `stat -f %c "$1"` > `stat -f %c "$2"` ]];
then
echo "$1"
else
echo "$2"
fi
Note: It seems stat is infamously implementation dependent so the format/parameters may be different for your machine.
Temporarily remove mosek from matlab path to access overwritten optimization toolbox quadprog
Saturday, February 12th, 2011I’ve been using mosek from matlab, which redefines a number of functions (quadprog, optimget, optimset, etc.) that were originally functions from matlab’s optimization toolbox. I wanted to compare the functionality and quality of mosek’s quadprog vs matlab’s quadprog, but this turned out to be a pain to do. Since I have all of the mosek paths at the top of my matlabpath, its versions are seen first thus the matlab versions are hidden. Even worse if all I want to do is call the matlab version of quadprog I can’t simply get a handle to the matlab quadprog and call it since internally matlab’s quadprog calls optimget and optimset which have also been overwritten by mosek.
If mosek is at the top of your path the you should be able to issue this:
which quadprog
and see something like:
/opt/local/mosek/6/toolbox/r2007a/quadprog.m
Thus the easiest solution was to temporarily eliminate mosek from the matlab path. This is possible programatically with:
% store original path so we can reset to it after we're done
old_path = matlabpath;
% split based on colon to obtain cell array of individual directories
dirs=regexp(old_path,':','split');
% only keep directories that don't contain mosek
new_dirs = dirs(cellfun(@isempty,regexp(dirs,'mosek')));
% append colons to each directory name and concatenate
new_path = cell2mat(strcat(new_dirs,':'));
% set matlab path to new path
matlabpath(new_path);
Now if you issue:
which quadprog
you should see something like:
/Applications/MATLAB_R2009aSV.app/toolbox/optim/optim/quadprog.m
Using std::map and std::pair to store attributes at graph edges in C++
Thursday, February 10th, 2011
I’ve been really battling with my old code in which I stored a graph by keeping a list of vertices and each vertex kept track of its neighbors. This is fine until you want to store attributes at the graph’s edges.
My solution is to store the edges in a std::map that is indexed by a pair of vertices. Here’s a boiled-down example. For know my “vertices” are just represented as ints (their indices, for example). And I just want to store a single string at each edge (i,j). But you could change these to be their own structs or classes. I’ve also tried not to get too fancy with the overloading and extending the std::pair class which I use as a key. It’s possible to do all sorts of tricks but this is the bare bones.
Save this in a file called edge_map.cpp:
#include <map>
#include <string>
#include <cstdio>
// Make a short name for the edge map's key
typedef std::pair<int,int> EMK;
// Make a short name for the type stored at each edge, the edge map's value
typedef std::string EMV;
// Make a short name for the edge map
typedef std::map<EMK,EMV> EdgeMap;
int main()
{
// Store some edges with values
EdgeMap edge_map;
edge_map[EMK(1,2)] = "Purple";
edge_map[EMK(2,3)] = "Blue";
edge_map[EMK(2,4)] = "Yellow";
edge_map[EMK(2,5)] = "Red";
edge_map[EMK(5,3)] = "Green";
printf("Edge map contents:\n");
// Loop over edges in edge_map
for(
EdgeMap::const_iterator emit = edge_map.begin();
emit != edge_map.end();
emit++)
{
printf("(%d,%d) %s\n",
emit->first.first,
emit->first.second,
emit->second.c_str());
}
printf("\n");
printf("Test for edges:\n");
// Edge queries
printf("(%d,%d): %s\n",1,3,(edge_map.count(EMK(1,3)) == 0 ? "No" : "Yes"));
printf("(%d,%d): %s\n",1,5,(edge_map.count(EMK(1,5)) == 0 ? "No" : "Yes"));
printf("(%d,%d): %s\n",1,2,(edge_map.count(EMK(1,2)) == 0 ? "No" : "Yes"));
return 0;
}
Compile with:
g++ edge_map.cpp -o edge_map
And the running ./edge_map should print:
Edge map contents:
(1,2) Purple
(2,3) Blue
(2,4) Yellow
(2,5) Red
(5,3) Green
Test for edges:
(1,3): No
(1,5): No
(1,2): Yes
Note: With fancy operator overloading I think its possible to store an undirected graph so that queries in either direction return correctly (the above code is good for directed graphs or where the caller manages the “undirectedness”).