# Rotation by Shearing

(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(Θ)

1+αβ=cos(Θ)
1+αsin(Θ)=cos(Θ)
α=(cos(Θ)-1)/sin(Θ)
α=-tan(Θ/2)

1+αβ=1+βγ
γ=α

Thus we have α=γ=-tan(Θ/2) and &beta = sin(Θ).

### Example

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:

That's it!

### Implementation

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 `P`, and `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
```

Tobin Fricke