(Or, "how do I rotate a bitmap?")
What follows are my own notes on Alan Paeth's ``A Fast Algorithm for General Raster Rotation,'' as published in the proceedings of Graphics Interface '86. This is a very popular algorithm for image rotation, used by many libraries such as ImageMagick, pnmrotate, etc. Initially I hoped that this scheme would be suitable for use in scientific data processing where arrays must be rotated, but now I believe that this scheme is unsuitable due to the poor interpolation. I also question whether this is any faster than bilinear interpolation when run on general purpose hardware (due to the large amount of data movement); however I do not yet have any quantitative measurements. These concerns aside, the decomposition of a rotation into three shears is interesting in its own right. — Tobin 2002-07-07
A two dimensional shear operation axis has the following matrix representations (one shear matrix for a shear parallel to the X axis, and another for a shear parallel to the Y axis):
Shear-X(α)= [1 α ; 0 1 ]
Shear-Y(β)= [1 0 ; β 1 ]
Thus a shear in the X direction produces y':=y (unchanged) and x':=x+αy.
It turns out that any orthogonal transformation can be realised by a combination of a shear along one axis, then a shear along the other axis, followed by another shear along the first axis: (although that is not proved here)
ShearX(α) ShearY(β) ShearX(γ)= [1 α ; 0 1 ] [1 0; β 1 ] [1 γ ; 0 1 ] = [1+αβ, α+γ+αβγ ; β, 1+βγ ]
Doing a rotation by performing three shear operations might be advantageous, because it's easy to do a shear operation. To do a shear operation on a raster image (that is to say, a bitmap), we just shift all the pixels in a given row (column) by an easy-to-calculate displacement. But in order to do a rotation using shears, we'll have to be able to calculate the necessary values of α, β, and γ from the rotation angle Θ.
Recall the familiar rotation matrix:
Rotate(Θ)= [ cos(Θ), -sin(Θ) ; sin(&Theta), cos(&Theta) ]
Set the rotation matrix equal to the product of the three shears:
[1+αβ, α+γ+αβγ ; β, 1+βγ ] = [ cos(Θ), -sin(Θ) ; sin(&Theta), cos(&Theta) ]
Solve for α, β, & γ in terms of Θ:
&beta = sin(Θ)
Thus we have α=γ=-tan(Θ/2) and &beta = sin(Θ).
Let Θ=20°. Then α~-0.18 and β~0.34. Here you can see what a square looks like as the three shears are applied:
Just for fun we can superimpose the original box, the intermediate results, and the final image:
Alan Paeth presents code in a Pascal-like language to do a shear in an efficient manner, using a simple approximation to do some interpolation (see his paper for the details of the approximation). Here's the code:
PROCEDURE xshear(shear, width, height) FOR y:= 0 TO height-1 DO skew := shear * (y + 0.5); skewi := floor(skew); skewf := frac(skew); /* = skew - skewi */ oleft := 0; FOR x:= 0 TO width-1 DO pixel := P(width-x, y); left := pixmult(pixel, skewf); /* pixel - left = right */ pixel := pixel - left + oleft; P(width-x+skewi, y) := pixel; oleft := left; OD P(skewi, y) := oleft; OD END
You can just replace the
pixmult function with multiplication in most cases; its function is just to weight the given pixel value for the "
skewf" fraction. The array to be sheared is given as the variable
shear is the shear factor (i.e. α, β, or γ). To make this useful you'll have to incorporate a change-of-origin to the center of rotation, and a suitable displacement to deal with negative indicies.
Here's the beginnings of an implementation in matlab: (xshear.m)
function data_out = xshear(data,shear,antialias) % function data = xshear(data,shear) % % $Id: xshear.m,v 1.1 2002/07/06 17:38:37 tobin Exp $ the_size = size(data); width = the_size(1); height = the_size(2); for y=1:height, skew = shear * y; % how much to shift this line skewi = floor(skew); % integer part skewf = skew - skewi; % fractional part oleft = 0; for x=width:-1:1, % from right to left... pixel = data(x, y); % sample the pixel if antialias, left = pixel * skewf; else left = 0; end pixel = (pixel - left) + oleft; if (x+skewi > 0), data_out(x+skewi,y) = pixel; end oleft = left; end if (skewi+1 > 0), data_out(skewi+1,y) = oleft; end end