Archive for December, 2013

stl `bind` gotcha

Tuesday, December 31st, 2013

I got burned again by naming a function in c++ bind and then while using the std namespace accidentally calling std::bind instead of my bind. The problem is that std::bind seems to be templated and defined so that it takes any number of any possible type of argument. Therefor it will always compile the arguments meant for my bind function.

One solution would be to avoid naming functions already taken by stl. Instead I’m just ensuring my function is found by calling it using the global namespace ::bind.

Joe Pesci Honks His Own Head

Sunday, December 29th, 2013

Simple color contrasting color to grayscale conversion via PCA

Sunday, December 29th, 2013

I remembered a nice paper about trying to do a better job converting color images to grayscale so that parts of a color image with similar intensities but different hues show up contrasted in the grayscale version. They solve this with a dense global optimization. There have been a few follow ups I haven’t read yet. But I thought I’d try out a really simple, perhaps naive, conversion technique. I simply convert an RGB image to Lab (or ntsc’s yiq) colorspace and then conduct principal component analysis (PCA) an a weighted combination of these vectors. This gives me a poor man’s multi-dimensional scaling of the supposedly perceptually metric LAB colorspace. Actually with the right weighting parameters this works pretty well.

monet sunset pca luma gray

From left to right, compare the original, PCA grayscale, bare “Lightness” L channel, matlab’s rgb2gray.

Here’s the matlab code to compute this:

im = im2double(imread('~/Downloads/sunset.png'));
%lab = rgb2lab(im);
lab = rgb2ntsc(im); % not really lab but similar idea
f = 0.45;
wlab = reshape(bsxfun(@times,cat(3,1-f,f/2,f/2),lab),[],3);
[C,S] = pca(wlab);
S = reshape(S,size(lab));
S = S(:,:,1);
gray = (S-min(S(:)))./(max(S(:))-min(S(:)));

Here’re some various results for different weighting factors f.

monet sunset pca parameter range

Implementing forward kinematics for linear blend skinning

Saturday, December 28th, 2013

I’ve implemented linear blend skinning and forward kinematics many times now. Each time I start out I have to scratch my head a little about what gets stored where and how best to construct the final bone transformations. Here’s a little explanation that will hopefully help in the future.

First let’s recall the linear blend skinning formula:

v' = ∑ wi(v) Ti * v

where v' is the new position of some point v on our shape, and wi(v) is bone i’s weight function evaluated at v and Ti is bone i’s transformation.

We’re specifically dealing with rigid bones, in fact we’ll consider our skeleton a rigid kinematic chain. Then the formula becomes:

v' = ∑ wi(v) (Qi * v + ti)

where Qi is the absolute bone (or joint) rotation (e.g. stored as a quaternion) of bone i. These absolute rotations are just a concatenation down from the root to bone i of the relative rotations we’re storing for each bone. These relative rotations are what we’re changing during an animation. The translation part for each bone ti is what gets a tad confusing. It is a function of the static (rest) skeleton (specifically the rigid bone lengths and relative positions in space) and the center of each bone rotation (i.e. joints positions).

Let’s go through a small example. Suppose we have two bones and we’re storing relative rotations q₀ and q₁.

forward kinematics

First we’ll figure out Q₀ and t₀. We want our rotation q₀ to be centered about the bone’s “joint” or upper endpoint. To do this we first translate so that endpoint is at the origin,

forward kinematics

then rotate by q₀,

forward kinematics

and then translate back,

forward kinematics

This affectively computes Q₀ and t₀. Imagine we apply these operations to some arbitrary point x. First translate so that the endpoint is at the origin: x-x₀, then rotate q₀(x-x₀), then translate back q₀(x-x₀)+x₀. Now separate the rotation and translation and you have q₀x-q₀x₀+x₀, revealing that Q₀ = q₀ and t₀ = -q₀x₀+x₀.

Now, a bit more tricky let’s consider then next bone. Again we want to rotate about the bone’s parent joint, so translate such that it’s at the origin,

forward kinematics

then rotate, but remember we need to rotate not just for this bone but also its parent (and in general all its ancestors),

forward kinematics

then finally translate back, but now we want the center to end up connected with its parent bone, so translate to the transformed parent joint location,

forward kinematics

Again, this computes Q₁ and t₁. Let’s go through it again for some arbitrary point x. First translate so that the parent joint is at the origin x-x₁, then rotate by the parent’s rotation and this bone’s q₁q₀(x-x₁), then translate to meet the parent’s transformation of the joint q₁q₀(x-x₁) + (q₀x₁-q₀x₀+x₀). Separating the rotation from the translation we have q₁q₀x - q₁q₀x₁ + q₀x₁-q₀x₀+x₀, revealing that Q₁ = q₁q₀ and t₁ = - q₁q₀x₁ + q₀x₁-q₀x₀+x₀.

Finally, putting both bones together we see the fully transformed skeleton.

forward kinematics

In general, we can write the absolute rotations and translations determined by forward kinematics as a recursive function of the relative rotations and rest joint locations:

Qi = qi Qp 
ti = - Qi xi + Qp xi - Qp xp + xp

where bone p is the parent of bone i.

Note: If you’re storing your skeleton as a tree where each bone i knows its parent p then one affective and efficient way to compute the result of this recursive formula is using dynamic programming. Alternatively, if each bone p knows each of its children i then one can just use breadth or depth first traversal.

Note: This is just for determining the absolute rigid transformations corresponding to each bone. It does not depend on the fact that we’re using linear blend skinning. The same construction applies if fancier skinning formulae are used like dual quaternion skinning.

Free wifi at Charles de Gaulle airport on laptop

Wednesday, December 25th, 2013

The Charles de Gaulle airport in Paris offers limited free wifi for mobile devices. To use this on your laptop just open your browser and go to the mobile start page: mobile site

I found this by switching my browser to use the iPhone user agent and opening a random page.

3D printable cookie cutters from any input shape (curve)

Monday, December 23rd, 2013

Here’s my pipeline for automatically converting a drawing of a shape into a 3d printed cookie cutter.

Today I printed a cookie-shaped cookie cutter. I photographed a drawing of the shape:

cookie cookie cutter drawing

Then I converted this to a curve using MATLAB, robustly created an offset curve, extruded it with a small lip.

cookie cookie cutter cutter

I also extrude an slightly more shrunken offset of the curve to make a “stamper” to help push the cookie out of the “cutter” above:

cookie cookie cutter stamper

Notice that these are oriented so that they are height fields in the y-direction, making them easy to print with the makerbot. The stamper fits inside, under the lip of cutter.

After 2 hours the printing is done. There are a few errors, but it more or less looks and functions like it should.

Pictures of the cookies to come.

To try it out yourself, download this zip of matlab files and run:

cookie_cookie_cutter

You’ll need to grab the winding number (a bit of overkill) code, if you don’t already. Then you can try out this completely procedure pipeline on your own shape (curve).

Update: Here’s the edible proof:

cookie-shaped cookie

Cubic Bezier curve in matlab

Friday, December 20th, 2013

Here’s a nasty anonymous function to compute points along a cubic bezier curve.

For example, the formula for a cubic is pretty compact. If you have some parametric samples t and knot points in the rows of P then:

bez = @(t,P) ...
  bsxfun(@times,(1-t).^3,P(1,:)) + ...
  bsxfun(@times,3*(1-t).^2.*t,P(2,:)) + ...
  bsxfun(@times,3*(1-t).^1.*t.^2,P(3,:)) + ...
  bsxfun(@times,t.^3,P(4,:));

builds an Anonymous Function so that X = bez(t,P) gives you points along the curve.

For a concrete example try:

t = linspace(0,1)';
P = [0 0;0.6 0.1;0.9 0.4;1 1];
X = bez(t,P);
plot(X(:,1),X(:,2));
hold on;
plot(P(:,1),P(:,2),'o-r')
hold off;

which displays the curve in blue and the knots in red:

cubic bezier curve

Reshape 3D array to match SPSS long format

Thursday, December 19th, 2013

The input format for repeated measures with two-conditions in SPSS is a little funny. In MATLAB, I would to store this data in a 3D m-by-n-by-p array where m is the number of variants across condition B, n is the number of variants across the condition in the rows, m across columns and p the number of repetitions (trials) across tubes. However, the wide format of SPSS expects repetitions to all be on the same row with flags indicating which row and column this tube corresponds to.

So if X is your m-by-n-by-p array then create a reshaped version of X and also create some index arrays:

sX = reshape(permute(X,[2 1 3]),size(X,1)*size(X,2),size(X,3));
srows = reshape(permute(rows,[2 1 3]),size(X,1)*size(X,2),1);
scols = reshape(permute(cols,[2 1 3]),size(X,1)*size(X,2),1);

Then you can combine these into one matrix:

SPSS = [srows scols sX];

So if your original data is:

>> X

X(:,:,1) =

   111   121
   211   221
   311   321


X(:,:,2) =

   112   122
   212   222
   312   322


X(:,:,3) =

   113   123
   213   223
   313   323


X(:,:,4) =

   114   124
   214   224
   314   324

Then this outputs:

>> SPSS

SPSS =

     1     1   111   112   113   114
     1     2   121   122   123   124
     2     1   211   212   213   214
     2     2   221   222   223   224
     3     1   311   312   313   314
     3     2   321   322   323   324

I guess this is a tensor equivalent of IJV format often used for sparse matrices.

Update: Here’s a nasty anonymous function to convert directly:

to_spss = @(X) [reshape(permute(repmat((1:size(X,1))',[1 size(X,2)]),[2 1 3]),size(X,1)*size(X,2),1) reshape(permute(repmat((1:size(X,2)),[size(X,1) 1 1]),[2 1 3]),size(X,1)*size(X,2),1) reshape(permute(X,[2 1 3]),size(X,1)*size(X,2),size(X,3))];

Then you can just call

SPSS = to_spss(X);

Reshape 3D array to match expected anova2 input

Thursday, December 19th, 2013

The anova2 function of MATLAB allows you to load in observations across two conditions with repetitions. The input format is a little funny though. I would expect to store this data in a 3D m-by-n-by-p array where m is the number of variants across condition B, n is the number of variants across condition A and p is the number of repetitions (actually I would expect the rows to be called “A” and columns “B” but help anova2 labels them this way, so I will, too). However, MATLAB expects repetitions to come as rows under each corresponding variant of condition B.

So to reshape a m-by-n-by-p matrix X use:

Y = reshape(permute(X,[3 1 2]),size(X,1)*size(X,3),size(X,2));

then you can call:

p = anova2(Y,size(X,3));

Update: Here’s a nasty Anonymous Function:

to_anova = @(X) reshape(permute(X(:,:,:),[3 1 2]),size(X,1)*size(X,3)*size(X,4),size(X,2));

then you can call:

p = anova2(to_anova(X),size(X,3));

Append a progress bar to a video or image sequence and show frame reordering

Thursday, December 19th, 2013

For a project we’re reordering the frames of a video. Here’s how to append a progress bar that visualizes where the frames are coming from in the original video.

Load in a video, here I use a video object:

video = VideoReader('rhinos.avi');

Let’s define a vector of indices I to be a remapping of the video frames so that frame f of the new video will be frame I(f) of the original video. We’ll use the original ordering

I = 1:video.NumberOfFrames;

or a random ordering:

I = randperm(video.NumberOfFrames);

Then append a progress bar and store the composite as an Height-by-Width-by-3-by-NumberOfFrames image sequence matrix S:

bar_w = video.Width;
bar_h = 20;
S = zeros(video.Height+bar_h,video.Width,3,video.NumberOfFrames);
for f = 1:video.NumberOfFrames
  S(1:video.Height,:,:,f) = read(video,f);
  col_If = round((I(f)-1)/(video.NumberOfFrames-1) * (bar_w-1)) + 1;
  col_f = round((f-1)/(video.NumberOfFrames-1) * (bar_w-1)) + 1;
  S(video.Height+(1:bar_h),1:col_f,:,f) = 0.2;
  S(video.Height+(1:bar_h),col_If,1,f) = 1;
end

Here’s what the result looks like for the original ordering:

original rhino ordering with progress bar

And with the random ordering:

remapped rhino ordering with progress bar

Note: This assumes that the remapped video and the original video are the same length.