Saturday, March 10, 2012

2D Affine Transformation as Action on a Triangle

Okay, in this post we'll move on from 2D linear transformations to 2D affine transformations.

The main motivation here is that we want a to be working with a transform algebra that has all of the "nice" behaviors of linear transforms (preservation of co-linearity of points and ratios of distances along lines, composition properties, matrix representation, etc.) but dispenses with the inconvenient and anti-intuitive fixed origin. That is, we'd like to include translations as well as the rotations, dilations, shears, and reflections that we already had at hand with linear transformations, per the previous post.

We can proceed algebraically as follows: we get a transform of the desired type (termed affine) by taking a linear transformation matrix \(A\) and following it up with translation by a vector \(\mathbf{b}\):\[\mathbf{T(x)} = A\mathbf{x} + \mathbf{b}\]It turns out that we can accomplish this computationally with single matrix multiplication via the following trick: we augment all vectors with a "1" at the bottom, and augment all matrices with a row of zeros at the bottom and an extra column to the right for the translation vector. This renders the above equation as:\[\left(\begin{array}{cc}\mathbf{T(x)}\\1\end{array}\right) = \left(\begin{array}{cc}A&\mathbf{b}\\0,...,0&1\end{array}\right)\left(\begin{array}{cc}\mathbf{x}\\1\end{array}\right)\]This checks out if you grind through the matrix multiplication; as you run across each row of the augmented matrix and down the augmented column vector, you end up with the usual computation for the corresponding component of \(\mathbf{T(x)}\), plus a final term of \(b_i\cdot1\) which accomplishes the translation. The final row computes as \(0 + 0 ... + 1\), providing the "1" at the bottom of the computed result vector.

All well and good, but there is a deeper geometric interpretation lurking here as well. All of our vectors and matrices have grown by an extra component, so it appears we are now operating in a space of one higher dimension. All the matrix coefficients are still real, so an affine transformation in, say, \(\mathbb{R}^2\) seems to be uniquely identified with a linear transformation in \(\mathbb{R}^3\). Why is this so?

If you consider the "1" added to the bottom of each column vector, what has been done is that the x-y plane in \(\mathbb{R}^2\) has been mapped to the z=1 plane in \(\mathbb{R}^3\).  So our original origin in \(\mathbb{R}^2\) has been moved off the problematic origin of \(\mathbb{R}^3\) to (0, 0, 1) and can now be moved about freely by a linear transform in \(\mathbb{R}^3\); in fact, a z shear in \(\mathbb{R}^3\) will accomplish a translation of the entire z=1 plane, as illustrated here (forgive the odd transparency artifacts; order-independent transparency is a tough problem and three.js doesn't go there in order to keep complexity down).

For those interested, coordinates transformed in this way to a fixed subspace of a higher dimensional space are known as homogeneous coordinates, more generally belonging to the study of projective geometry.  That general study is beyond the scope of what I'd like to cover here, so for now I'm just cherry-picking some useful techniques and points-of-view.

Okay, so, next, let's consider the higher dimensional analogue of the computational technique described in the previous post for divining a transformation matrix from a set of input and output points. In the previous (linear 2D) case, two input points and two output points allowed us to calculate a unique 2x2 transformation matrix. Two points to two points admitted a characterization of 2D linear transforms by their action on a line segment. To do the same trick for a linear transformation matrix in \(\mathbb{R}^3\) we would require 3 input points and 3 output points. Thus, a 3D linear transform may be characterized by its action on a triangle.

Furthermore, if we restrict the points of our input and output triangles to the z=1 plane, the transformation characterized by the resulting calculated matrix will be closed in the z=1 plane (any point in the z=1 plane transformed by the resulting matrix will be carried to some other point in the z=1 plane). This follows from the fact that the transformation implemented by our result matrix must be at least linear in 3D, combined with the co-linearity preservation quality of all linear transformations (to see this, characterize any input point in the z=1 plane as a weighted vector sum of two of the sides of the input triangle).

Is such a calculated transformation actually affine in the z=1 plane, however? As it turns out, yes. This is easiest to see (well, I think!) by examination of the decomposition taxonomy of linear transformations to scales, shears, reflections, stretches, dilations, and rotations. Consider the effects of each of these restricted to cases were the z=1 plane is closed, and we arrive at a decomposition taxonomy equivalent to a linear transformation in the z=1 plane followed by a translation in the z=1 plane.

Okay, so here is a web toy for playing around with 2D affine transformations. You can click in the image to drag out a couple of triangles, then repeatedly apply the transformation whose action maps the first triangle to the second one, operating on one of several synthetic images. Click the "+" to go offsite and investigate the code.

Compare this to the linear transformation web toy in the previous post to get a feel for the difference between linear and affine transformations. With affine transformations at our disposal, we are almost ready for some geodesic dome programs! More on this in an upcoming post.

Tuesday, January 3, 2012

2D Linear Transformation as Action on a Line Segment

Okay, this post is about 2D linear transformations as a stepping stone to 3D affine transformations for geodesic dome programs. I think a lot of people miss a useful bit of understanding about affine transformations: that they are just linear transformations in a higher-dimensional space. Some familiarity with the properties and "feel" of linear transformations first is helpful, and there is less bookkeeping and it is easier to see what is going on in 2D.

We are going to be working with linear transformations (maps) from the 2D real plane (\(\mathbb{R}^2\) viewed as a vector space) to itself. Given a transformation \(\mathbf{T}\), vectors \(\mathbf{a}\), \(\mathbf{b}\) and scalar \(\alpha\), we have the usual general linearity conditions:

  • \(\mathbf{T(a} + \mathbf{b)} = \mathbf{T(a)} + \mathbf{T(b)}\)
  • \(\mathbf{T(}\alpha\mathbf{a)} = \alpha\mathbf{T(a)}\)

From these, we can derive in particular over \(\mathbb{R}^2\) (though these also generalize to higher dimensional spaces):

  • The origin (0,0) is a fixed point.
  • Compositions of linear transformations are also linear.
  • Inverses of linear transformations (if they exist) are also linear.
  • Linear transformations on \(\mathbb{R}^2\) are representable as 2x2 matrices with real coefficients (and any such 2x2 matrix implements a linear transformation), operating on column vectors in \(\mathbb{R}^2\).
  • Linear transformations preserve co-linearity (lines map to lines) and proportional distance between colinear points.
  • A linear transformation on \(\mathbb{R}^2\) is completely specified by its action on two points, or, equivalently, its action on a single line segment.

This last property will be of interest to us -- it will turn out to generalize nicely to affine transformations and higher dimensional spaces. This, and the co-linearity property, will be directly useful for geodesic dome programs.

The big pain/limitation/disappointment with linear transformations is that they leave the origin fixed. This is what will be overcome by moving to affine transformations.

Here is a web toy for playing around with 2D linear transformations. You can click in the image to specify endpoints for a couple of line segments, then repeatedly apply the transformation whose action maps the first line segment to the second one, operating on one of several synthetic images. Click the "+" to go offsite and investigate the code (computation explanation follows). Some patience required; depending on your browser it can be a bit slow (it's all client side JavaScript).

So, if you play around with this a bit, most of the qualities of linear transformations are noticeable. Note particularly the "inconvenience" of having the origin fixed. Transformations defined by segments that are co-linear or close to co-linear with the origin are singular or nearly so.

You can also get a feel here for that fact that linear transformations have a pretty simple taxonomy/decomposition -- there are only reflections, rotations, scales, stretches (unequal scale in x and y) shears (flattening of squares to parallelograms), and combinations of these.

The above toy works computationally by representing transformations as 2x2 real matrices. The matrix operations are "unrolled" into primitive algebraic operations in the code in lieu of using a prefab linear algebra JavaScript package, mostly for performance reasons. A matrix is generated for the inverse of the transformation that we want; then each desired point in the output image is transformed back through that matrix to find its corresponding point in the synthetic test subject.

The most interesting calculation here is the one that yields the matrix of a transform, given the endpoints for the source and transformed line segments. This is a useful computational trick if you haven't seen it before. It works like this:

Suppose we are looking for the transformation that maps vectors \(\mathbf{a},\mathbf{b}\) to vectors \(\mathbf{c},\mathbf{d}\), that is:\[\begin{align*}(a_x,a_y) &\mapsto (c_x,c_y)\\(b_x,b_y) &\mapsto (d_x,d_y)\end{align*}\]
First, consider the matrix formed from adjoining the two source points \(\mathbf{a},\mathbf{b}\), and it's inverse. We have:\[\left(\begin{array}{cc}a_x&b_x\\a_y&b_y\end{array}\right)^{-1}\left(\begin{array}{cc}a_x&b_x\\a_y&b_y\end{array}\right) = \left(\begin{array}{cc}1&0\\0&1\end{array}\right)\]
Since a column of the result matrix on right is the result of the leftmost matrix operating on the corresponding column of the middle matrix (definition of matrix multiplication), we can see the the leftmost matrix here sends each of our source points to an orthonormal basis vector. That is:\[\begin{align*}\mathbf{a}&\mapsto(1,0)\\\mathbf{b}&\mapsto(0,1)\end{align*}\]
Likewise, we have:\[\left(\begin{array}{cc}c_x&d_x\\c_y&d_y\end{array}\right)\left(\begin{array}{cc}1&0\\0&1\end{array}\right) = \left(\begin{array}{cc}c_x&d_x\\c_y&d_y\end{array}\right)\]
so the leftmost matrix here sends each basis vector on to one of our destination points:\[\begin{align*}(1,0)&\mapsto\mathbf{c}\\(0,1)&\mapsto\mathbf{d}\end{align*}\]
Composing these two matrices (destination point matrix multiplied by the inverse of source point matrix) gives the matrix result which will carry the source points to the destination points and all other points along for the ride, obeying the constraints/behaviors of linear transformations. This same method is perfectly generalizable to higher dimensional spaces, and will be one of the principal techniques used in the upcoming geodesic dome programs.

Saturday, December 3, 2011

Synthetic Test Subjects

Been quite busy lately, so not much time to write here.  Recently, though, a work colleague brought up a hobby project that touched on geodesic dome geometry, which is something I've played around with a bit and had been intending to explore here.  That's prodded me to take find a little time to write another post or two.

The programs I've written dealing with geodesic dome geometry make a lot of use of projection via affine transform.  Basic stuff, but I think interesting and visual in itself and often poorly explained.  So I might write up a post or two and put together a couple of web toys just dealing with that as some groundwork, before moving on to anything having to do with geodesic domes directly.

A web toy dealing with projections/mapping/warping needs a visual test subject.  A usual approach is to take a sample bitmap and show the result of running it through the transform.  This can produce less than appealing results, because of the inherent limited resolution of the source bitmap.  Depending on the transform, very small (subpixel) areas of the source image can end up taking up very large areas of the result image.  Intermediate pixels then need to be synthesized/interpolated, and the results get bumpy/blurry.

I find it more convenient to deal with a synthetic, functional, test subject.  Given a point (x,y) in the real image plane, a function is invoked to return the color of that point in test subject.  This sort of test subject can have effectively infinite (within the limits of floating point) resolution.

Since the test subject has effectively infinite resolution, it is also easy to oversample the source with respect to the output bitplane and then average the sub-pixel results to obtain a greater effective visual resolution.

Code for a synthetic test image need not be complicated.  Here are some coded synthetic test subjects in just a dozen or so lines of JavaScript each, and also a couple functions to oversample and render the results (hosted at; click the "+" button to go offsite and explore the code).

I'll be using these test subjects in subsequent articles to illustrate the effects of transform math.

Thursday, September 8, 2011

Three.js and jsFiddle

So, the previous post touched on the fact that all of the free or open-source math and computer algebra environments that I've been looking at are a bit weak in the general 3D polygon mesh visualization department.  Since the purpose of the blog is to put up things people can immediately see and play with, I thought I'd see what was at hand these days for 3D visualization directly in the web browser.  Without too much effort, I found Three.js and jsFiddle.  Taken together, they can make toys like this:

Here, Three.js provides a nice browser-neutral 3D JavaScript library.  It can render via any of WebGL, HTML5 Canvas, or SVG back ends, to various levels of capability.  I've chosen SVG here, so this example should even work in Internet Eplorer :-)

JsFiddle provides a surrounding environment for web toys.  If you click on the "+" button above, you'll go off site to a full multi-pane view of the toy, where you can inspect and modify and play with the various parts of the HTML, CSS, and JavaScript.

I did hit one wrinkle to do with Internet Explorer 9, which I'll mention here in case anybody else runs up against it.  Originally, I was referencing Three.js from its canonical location on GitHub.  GitHub was exporting this with a MIME type of "text/plain", however, and it seems IE9 will no longer load <script> tags whose specified type does not match the MIME type of their source.  In the end, I needed to copy the script files I wanted to use to a different server that exports them as "application/javascript" to make things work correctly under IE9.

MATLAB alternatives

Okay, so I like MATLAB -- it's interactive, expressive, handles a lot of tedious math, and draws nice pictures. It's a little arcane (which can be fun in a sort of puzzling way) and it makes you think about everything in a sort of vectorized way to get decent performance (which can be fun in a sort of challenging way).

I used to have an old crusty copy for 32-bit Windows on some previous laptop from some previous job.  These days I'm using a MacBook (and 64-bit Windows when I have to) so I figured I'd go see about getting an updated copy for personal use...

$2,100.00 for a single, personal use license.  Are you f-ing kidding me?  I could buy a really nice bass for $2,100.00.  Y'know, for $200.00 I would have just bought MATLAB right away.  For anything up to $500.00 I probably would have just bought it after some serious thought...  For $2,100.00 -- no way!

I want to post snippets of mathematical code here from time to time, and I'd like students and non-millionaire-hobbyists to be able to conceivably grab those and play with them.  So to heck with MATLAB.  Started looking into free and/or open-source alternatives, and few that I played with were not that bad:

SciLab (

Octave (

NumPy / SciPy (

I'll comment more on these soon after I've had a little more of a chance to play around with them.  I'm probably leaning towards NumPy / SciPy, since so many people seem to be familiar with Python these days.

The most immediate weakness of all three of these seems to be the lack of a fully developed 3D polygon patch visualization system like MATLAB has (damn you, MATLAB!) so I guess I'll be looking elsewhere for that.

Wednesday, August 31, 2011

Moebius transformation animated GIFs

Here are some animated GIFs that I created a few years ago with MATLAB.  These characterize the action of the four classes of Moebius transformations, mapping the complex plane to itself.  They were inspired by the illustrations and analysis in a section of Tristan Needham's excellent book Visual Complex Analysis:


So what's it all about?  A Moebius transformation is a mapping on the complex numbers of the form \[ M(z) = \frac{az + b}{cz + d}, \] where \(a, b, c, d\) are complex constants.  Moebius tranformations have many nice features: they are one-to-one and onto on the complex plane (extended with the addition of the "point at \(\infty\)"), they form a group under composition, they are conformal (angle preserving).  Interestingly, though conformality only implies that they map infinitesimal circles to other infinitesimal circles, Moebius transformations actually map circles of any size in the complex plane to other circles in the plane.  For a really nice exposition and proofs of these qualities, see Needham.

Fixed points on the plane under a Moebius transform will simply be solutions of \[ z = \frac{az+b}{cz+d}, \] which is just quadratic in \(z\).  So a (non-identity) Moebius transform will have at most two fixed points, or one when there is a repeated root.  Following Needham, for the remainder of this discussion we'll refer to the fixed points of a given Moebius transformation as \(\xi_+\), \(\xi_-\) in the two root case, or just \(\xi\) in the repeated root case.

Most of the qualities above are readily observable in the animations: one or two fixed points, circles mapping to circles, one-to-one, easy to imagine extension to the whole plane.  Conformality is a little harder to see, but if you look closely you can see that the angles at each of the corners in these crazed checkerboards always remain 90°.

The GIFs were generated using a transform method.  Given a Moebius transformation with two fixed points (we will revisit the single fixed point case shortly), consider the additional Moebius transformation defined by \[ F(z) = \frac{z - \xi_+}{z - \xi_-}. \]This will send \(\xi_+\) to \(0\), and \(\xi_-\) to \(\infty\).  We can now construct \[ \widetilde{M} = F \circ M \circ F^{-1}, \] which is itself a Moebius transformation (since the Moebius transformations are a group under composition) and which will have fixed points at \(0\) and \(\infty\).  We can consider \(M\) and \(\widetilde{M}\) to be a transform pair; for any \(M\) there is a corresponding \(\widetilde{M}\) with fixed points at \(0\) and \(\infty\).

Now it turns out that Moebius transformations with fixed points at both \(0\) and \(\infty\) reduce to a particularly simple form — that of a single complex multiplication, i.e. just a dilation and/or rotation about the origin of the complex plane.  See Needham for exposition of this.  The elliptic, hyperbolic, and loxodromic classes of Moebius transformations turn out to be those whose corresponding \(\widetilde{M}\) is a pure rotation, pure dilation, or general combination of the two, respectively.

To generate these GIFs, we decorate the complex plane on the \(\widetilde{M}\) side like a circular checkerboard or dart board, and represent the action of a given class of \(\widetilde{M}\) transformations by rotating it, dilating it, or a combination of both, e.g.:

elliptic \(\widetilde{M}\)
hyperbolic \(\widetilde{M}\)

On the \(M\) side, we can pick fixed points \(\xi_+\), \(\xi_-\) wherever we like, and then derive the corresponding \(F\) as above.  For each frame, we take each point on the \(M\) side, map it through \(F\) to find the corresponding point on the \(\widetilde{M}\) side, and check that against a dynamic checkerboard model to see if the point should be black or white in this frame.

Returning to the repeated root, single fixed point case: we treat this similarly, but set \[ F = \frac{1}{z - \xi}, \] which sends \(\xi\) to \(\infty\).  As before, it turns out that Moebius transformations of this form (repeated fixed point at \(\infty\)) reduce to a very simple form: this time, a pure translation.  The parabolic class of Moebius transformations are those whose corresponding \(\widetilde{M}\) transformation is a pure translation.  To illustrate these, we use a translating dynamic checkerboard on the \(\widetilde{M}\) side that looks like this:

parabolic \(\widetilde{M}\)

All of this can be done quite succinctly in MATLAB, once it is understood what needs to be done.  Below is the MATLAB snippet which was used to generate these GIFs.  Commenting and uncommenting various lines chooses different \(F\) and \(\widetilde{M}\) side checkerboard models, resulting in the different outputs.

The first few lines create a 512x512 matrix of complex numbers, ranging from -3 to 3 on both the real and imaginary axes, to represent a portion of the complex plane.  This is then pre-mapped through an appropriate \(F\).  The for loop iterates on each frame.  The dynamic checkerboard result is calculated as the xor of various versions of functions g1 and g2 operating over the premapped points.  Each frame is downsampled and written to disk as a separate file, then at the end all the frames are stitched into a movie.  I must then have used some non-matlab utility to convert the movies to animated GIFs, but I'm not sure what that was...

Z = -3:6/512:3-(6/512); Z = Z(ones(1,length(Z)),:); Z = complex(Z,Z'); %Z = 1 ./ Z; % one finite fixed point Z = (Z+1)./(Z-1); % two finite fixed points im = []; g1 = []; g2 = []; for frame = 1:30 g1 = mod(log(abs(Z))*4,2)<1; % radial, static % g1 = mod(log(abs(Z))*4+(frame-1)/15,2)<1; % radial, dynamic % g1 = mod(real(Z)+(frame-1)*.8/30,.8)<.4; % vertical, dynamic % g2 = mod(angle(Z), pi/6)<(pi/12); % circumferential, static g2 = mod(angle(Z)+(frame-1)*pi/180, pi/6)<(pi/12); % circumferential, dynamic % g2 = mod(imag(Z),.8)<.4; % horizontal, static im(:,:,1,frame) = imresize(xor(g1,g2), .25, 'bilinear'); imwrite(im(:,:,1,frame), sprintf('frame%.2d.bmp', frame)); end mov = immovie(255*im,gray(255)); movie(mov,50);

Sunday, August 28, 2011


I'll be wanting a lot of mathematical content here, so took a look around to see what the current state of that art is.  Was delighted to find MathJax:
"No more setup for readers.  No more browser plugins.  No more font installations...  It just works."
Nice; the above tag-line is pretty much the exact opposite of MathML (ha!)  Basically, they've got client-side LaTeX typesetting via JavaScript.  I'm a big fan of LaTeX.  A few years ago I started working on my own server-side LaTeX-on-the-web solution via JSP, but it was rather involved and a bit of a pain.  MathJax seems to perform much better than my prototype system did, and has additional client-side niceties (dynamic scaling, cut and paste, etc.)

To get this up and running in a system like Blogger, one just sources a single script from MathJax' CDN in the blog template (see here for details).  Then type LaTeX source in your blog posts; MathJax will scrape it and re-render it right in the client's browser.  Here's a sample: \[ \int\limits_M{\partial\omega} = \int\limits_{\partial M}{\omega}. \]

There are a few niceties for poor suffering MathML users thrown in as well.