mathr / blog / #

In my work-in-progress et project for escape-time fractals, I currently represent the viewing transform (from pixel grid to complex plane coordinates) by a translation (high precision coordinates of the center of the view), a scaling (high range scale factor), and a 2×2 matrix that accounts for non-uniform scaling and rotation (4 low precision numbers, defaulting to the identity [1,0;0,1]) - this matrix should have determinant 1 as any global scaling belongs in the scale factor.

However, editing raw matrix values is not very user friendly, so I plan to add a friendly user interface based on marking points in the GUI before moving them around with the mouse (eventually: multi-touch support for tablets etc). An intermediate step might be to represent the matrix in a more human-relevant way, decomposing it into rotations and non-uniform scaling (a shear is a non-uniform scaling conjugated by a rotation, no need to handle it separately). It so happens that this is a well-known linear algebra problem, called polar decomposition. The matrix M is decomposed into a rotation V and a stretch P, such that M = V P, and further the stretch P can be decomposed into a rotation U and a diagonal matrix D, such that P = U D U^-1 (though this last decomposition is not unique).

A good description of the problem and examples in higher dimensions is:

Matrix Animation and Polar Decomposition

Ken Shoemake and Tom Duff

Abstract General 3×3 linear or 4×4 homogenous matrices can be formed by composing primitive matrices for translation, rotation, scale, shear, and perspective. Current 3-D computer graphics systems manipulate and interpolate parametric forms of these primitives to generate scenes and motion. For this and other reasons, decomposing a composite matrix in a meaningful way has been a long-standing challenge. This paper presents a theory and method for doing so, proposing that the central issue is rotation extraction, and that the best way to do that is Polar Decomposition. This method also is useful for renormalizing a rotation matrix containing excessive error.

For the 2D case there is a simple explicit formula given in:

Explicit polar decomposition and a near-characteristic polynomial: The 2×2 case

Frank Uhlig

Abstract Explicit algebraic formulas for the polar decomposition of a nonsingular real 2×2 matrix A are given, as well as a classification of all integer 2×2 matrices that admit a rational polar decomposition. These formulas lead to a functional identity which is satisfied by all nonsingular real 2×2 matrices A as well as by exactly one type of exceptional matrix A_n, for each n > 2.

Translated into Octave code (presumably Matlab compatible), assuming that M is real with det(M) > 0:

M = [1,-3;2,2]
scale = sqrt(det(M))
A = M / scale;
B = A + inv(A');
b = sqrt(abs(det(B)));
V = B / b;
P = (A' * A + eye(2)) / b;
[U,D] = eig(P);
stretch = D(1);
stretchAngle = atan2(U(2,1), U(1,1));
if (stretch < 1)
  stretch = 1 / stretch
  stretchAngle = stretchAngle + pi / 2;
endif
stretchA = mod(stretchAngle, pi)
rotation = atan2(V(2,1), V(1,1))
R = [ cos(rotation), -sin(rotation); sin(rotation), cos(rotation) ];
S = [ stretch, 0; 0, 1/stretch ];
T = [ cos(stretchA), -sin(stretchA); sin(stretchA), cos(stretchA) ];
N = scale * R * T * S * T';
error = norm(M - N)

Example output:

M =

   1  -3
   2   2

scale =  2.8284
stretch =  1.2808
stretchA =  1.4483
rotation =  1.0304
error =    1.0721e-15

Note that eigenvalues and eigenvectors can be found explicitly in 2D, see for example Eigenvalues and eigenvectors of 2x2 matrices.

Final things to note: currently et uses only the view scale factor to choose which number type to use for calculations. This might lead to pixelation artifacts from insufficient precision in highly stretched images near number type thresholds, so the stretch factor should be taken into account too when determining the minimal pixel spacing. Handling reflection (det < 0) is left for future investigation.