2D Matrix Decomposition

This demonstrates the polar decomposition of 2D matrices into their angle, scale, and shear components. The polarDecomp algorithm is taken from Shoemake & Duff "Matrix Animation and Polar Decomposition", found at <http://www.cs.wisc.edu/graphics/Courses/838-s2002/Papers/polar-decomp.pdf>, 1992.

I found that the formula for the polar decomposition, on page 3, produces rotation matrices that don't account for diagonal axis flipping (swapping x with y). The article seems to factor out the flipping, but I didn't want to have a fourth component. Removing the sign(det(M)) factor eliminated this problem by incorporating the flipping into the rotation matrix.

It's very desirable to decompose a matrix generated by {angle, scale, shear} back into the same part values, but this happens only when scale is positive and shear is zero. If not, the rotation is affected. At least, the parts are stable once derived from the first matrix. This may help in UI, and certainly in interpolation.

Matrix generators

All the matrix ops are Vec . Mat

Rotation

rmat[a_]:={{Cos[a],-Sin[a]},{Sin[a],Cos[a]}}

Scaling

smat[xs_,ys_]:={{xs,0},{0,ys}}

smat[{xs_,ys_}]:=smat[xs,ys]

Shear

hmat[xh_,yh_]:={{1,xh},{yh,1}}

hmat[{xh_,yh_}]:=hmat[xh,yh]

A handy pre-defined symbolic matrix

(mat={{a,b},{c,d}})//MatrixForm

Returns Q (the rotation matrix) and S, the stretch matrix

Unlike the version in Shoemake & Duff's polar_decomp paper, this one always returns a normal rotation matrix, and thereby accounts for axis flipping

polarDecomp[m:{{m11_,m12_},{m21_,m22_}}]:=
With[{q=m+{{m22,-m21},{-m12,m11}}},
With[{qs=(q/(Sqrt[((q[[1,1]])^2)+((q[[2,1]])^2)]))},
{Simplify[qs],Simplify[Transpose[qs].m]}]]

Returns the angle of a rotation matrix

angle[{{a_,b_},{c_,d_}}]:=ArcTan[d,c]

Returns the {x, y} scale of a stretch matrix

scale[{{a_,b_},{c_,d_}}]:={a,d}

Returns the {x, y} shear of a stretch matrix

shear[{{a_,b_},{c_,d_}}]:={(b/a),(c/d)}

Converts the Q and S matrices into {angle, scale, shear}

decompParts[{qMat_,sMat_}]:={angle[qMat],scale[sMat],shear[sMat]}

Composition to produce the parts directly from the original matrix

polarParts[mat_]:=decompParts[polarDecomp[mat]]

Recreate the original matrix from the parts

compParts[{ang_,scl_,shr_}]:=rmat[ang].smat[scl].hmat[shr]

Testing

This loses precision enough to require the SetPrecision

testDecomp[mat_]:=
With[{qs=polarDecomp[mat]},
With[{parts=decompParts[qs]},
SetPrecision[Simplify[mat],10]!=SetPrecision[Simplify[compParts[parts]],10]]]

Test that matrix decomposition and recomposition is equal to the original

genRandMats[n_]:=Table[Partition[RandomReal[{-10.0,10.0},4],2],{n}]

MatrixForm/@(failedMats=Select[genRandMats[100],testDecomp])

Test that the same components are recovered from matrices generated by random components

Note that they aren't the same unless scale is positive and shear is zero

genRandComponents[n_]:=Table[{RandomReal[{-π,π}],RandomReal[{0.01,2},2],{0,0}},{n}]

testRandComponents[parts_]:=
With[{m=compParts[parts]},
With[{p=polarParts[m]},
With[{cp1=Chop[parts],cp2=Chop[p]},
cp1≠cp2]]]

failedParts=Select[genRandComponents[10],testRandComponents]

Test that components are equal after two passes of decomp/comp. The 1st pass output is "canonical"

This always succeeds, so long as the precision is kept down

genRandComponents2Pass[n_]:=Table[{RandomReal[{-π,π}],RandomReal[{-10,10},2],RandomReal[{-1,1},2]},{n}]

test2PassComponents[parts_]:=
With[{m=compParts[parts]},
With[{p1=polarParts[m]},
With[{m2=compParts[p1]},
With[{p2=polarParts[m2]},
With[{cp1=SetPrecision[p1,10],cp2=SetPrecision[p2,10]},
cp1≠cp2]]]]]

failedParts=Select[genRandComponents2Pass[100],test2PassComponents]

Examples

Note that Geometric transform does M . V, but that doesn't matter so long as it's consistent

Draw a transformed unit square, colored to make it easy to see axis flipping

transformRect[xform_,opts:OptionsPattern[Graphics]]:=Graphics[{Thick,GeometricTransformation[Line[{{0,0},{0,1},{1,1},{1,0},{0,0}},VertexColors->{Red,Green,Blue,Yellow,Red}],xform]},opts,Axes->True,AxesOrigin->{0,0},AspectRatio->Automatic]

Draw some random examples

With[{randMats=genRandMats[4]},
GraphicsGrid[{transformRect/@randMats,transformRect/@compParts/@polarParts/@randMats}]]