new file mode 100644
--- /dev/null
+++ b/MakeShears.m
@@ -0,0 +1,7 @@
+function [S1, S2] = MakeShears(theta)
+
+S1 = eye(2);
+S2 = eye(2);
+
+S1(1,2) = -tan(theta/2);
+S2(2,1) = sin(theta);

new file mode 100644
--- /dev/null
+++ b/imrotate.m
@@ -0,0 +1,148 @@
+## -*- texinfo -*-
+## @deftypefn {Function File} {} imrotate(@var{M}, @var{theta}, @var{method}, @var{bbox})
+## Rotation of a 2D matrix.
+##
+## Applies a rotation of @var{THETA} degrees to matrix @var{M}.
+##
+## The @var{method} argument is not implemented, and is only included for compatibility with Matlab.
+## This function uses Fourier interpolation,
+## decomposing the rotation matrix into 3 shears.
+##
+## @var{bbox} can be either 'loose' or 'crop'.
+## 'loose' allows the image to grow to accomodate the rotated image.
+## 'crop' keeps the same size as the original, clipping any part of the image
+## that is moved outside the bounding box.
+## @end deftypefn
+
+## Author: Jeff Orchard <jjo@cs.sfu.ca>
+## Created: Oct. 14, 2002
+
+function fs = imrotate(f,theta,method,bbox)
+
+	if ( nargin == 2 )
+		method = "fourier";
+		bbox = "loose";
+	elseif ( nargin == 3 )
+		bbox = "loose";
+	endif
+
+	# Get original dimensions.
+	[ydim_orig, xdim_orig] = size(f);
+	xcentre_orig = (xdim_orig+1) / 2;
+	ycentre_orig = (ydim_orig+1) / 2;
+
+	# Pre-process the angle ===========================================================
+	# Whichever 90 degree multiple theta is closest to, that multiple of 90 will
+	# be implemented by rot90. The remainder will be done by shears.
+
+	# This ensures that 0 <= theta < 360.
+	theta = rem( rem(theta,360) + 360, 360 );
+
+	# This is a flag to keep track of 90-degree rotations.
+	perp = 0;
+
+	if ( theta>=0 && theta<=45 )
+		phi = theta;
+	elseif ( theta>45 && theta<=135 )
+		phi = theta - 90;
+		f = rot90(f,1);
+		perp = 1;
+	elseif ( theta>135 && theta<=225 )
+		phi = theta - 180;
+		f = rot90(f,2);
+	elseif ( theta>225 && theta<=315 )
+		phi = theta - 270;
+		f = rot90(f,3);
+		perp = 1;
+	else
+		phi = theta;
+	endif
+
+
+
+	if ( phi == 0 )
+		fs = f;
+		if ( strcmp(bbox,"loose") == 1 )
+			return;
+		else
+			xmax = xcentre_orig;
+			ymax = ycentre_orig;
+			if ( perp == 1 )
+				xmax = max([xmax ycentre_orig]);
+				ymax = max([ymax xcentre_orig]);
+				[ydim xdim] = size(fs);
+				xpad = ceil( xmax - (xdim+1)/2 );
+				ypad = ceil( ymax - (ydim+1)/2 );
+				fs = impad(fs, [xpad,xpad], [ypad,ypad], "zeros");
+			endif
+			xcentre_new = (size(fs,2)+1) / 2;
+			ycentre_new = (size(fs,1)+1) / 2;
+		endif
+	else
+
+		# At this point, we can assume -45<theta<45 (degrees)
+
+		phi = phi * pi / 180;
+		theta = theta * pi / 180;
+		R = [ cos(theta) -sin(theta) ; sin(theta) cos(theta) ];
+
+		# Find max of each dimension... this will be expanded for "loose" and "crop"
+		xmax = xcentre_orig;
+		ymax = ycentre_orig;
+
+		# If we don't want wrapping, we have to zeropad.
+		# Cropping will be done later, if necessary.
+		if ( strcmp(bbox, "wrap") == 0 )
+			corners = ( [ xdim_orig xdim_orig -xdim_orig -xdim_orig ; ydim_orig -ydim_orig ydim_orig -ydim_orig ] + 1 )/ 2;
+			rot_corners = R * corners;
+			xmax = max([xmax rot_corners(1,:)]);
+			ymax = max([ymax rot_corners(2,:)]);
+
+			# If we are doing a 90-degree rotation first, we need to make sure our
+			# image is large enough to hold the rot90 image as well.
+			if ( perp == 1 )
+				xmax = max([xmax ycentre_orig]);
+				ymax = max([ymax xcentre_orig]);
+			endif
+
+			[ydim xdim] = size(f);
+			xpad = ceil( xmax - xdim/2 );
+			ypad = ceil( ymax - ydim/2 );
+			f = impad(f, [xpad,xpad], [ypad,ypad], "zeros");
+			xcentre_new = (size(f,2)+1) / 2;
+			ycentre_new = (size(f,1)+1) / 2;
+		endif
+
+		size(f)
+		[S1 S2] = MakeShears(phi);
+
+		tic;
+		f1 = imshear(f, 'x', S1(1,2), 'wrap');
+		f2 = imshear(f1, 'y', S2(2,1), 'wrap');
+		fs = real( imshear(f2, 'x', S1(1,2), 'crop') );
+		toc
+	endif
+
+	if ( strcmp(bbox, "crop") == 1 )
+
+		# Translate the current centre to centre_orig
+		fs = imtranslate(fs, xcentre_orig-xcentre_new, -ycentre_orig+ycentre_new, "wrap");
+
+		# Crop to original dimensions
+		fs = fs(1:ydim_orig, 1:xdim_orig);
+
+	elseif ( strcmp(bbox, "loose") == 1 )
+
+		# Find tight bounds on size of rotated image
+		# These should all be positive, or 0.
+		xmax_loose = ceil( xcentre_new + max(rot_corners(1,:)) );
+		xmin_loose = floor( xcentre_new - max(rot_corners(1,:)) );
+		ymax_loose = ceil( ycentre_new + max(rot_corners(2,:)) );
+		ymin_loose = floor( ycentre_new - max(rot_corners(2,:)) );
+
+		fs = fs( (ymin_loose+1):(ymax_loose-1) , (xmin_loose+1):(xmax_loose-1) );
+
+	endif
+
+
+endfunction

new file mode 100644
--- /dev/null
+++ b/imshear.m
@@ -0,0 +1,114 @@
+## -*- texinfo -*-
+## @deftypefn {Function File} {} imshear (@var{M}, @var{axis}, @var{alpha}, @var{bbox})
+## Applies a shear to @var{M}.
+##
+## The argument @var{M} is either a matrix or an RGB image.
+##
+## @var{axis} is the axis along which the shear is to be applied, and can
+## be either 'x' or 'y'.
+## For example, to shear sideways is to shear along the 'x' axis. Choosing
+## 'y' causes an up/down shearing.
+##
+## @var{alpha} is the slope of the shear. For an 'x' shear, it is the 
+## horizontal shift (in pixels) applied to the pixel above the 
+## center. For a 'y' shear, it is the vertical shift (in pixels)
+## applied to the pixel just to the right of the center pixel.
+##
+## NOTE: @var{alpha} does NOT need to be an integer.
+##
+## @var{bbox} can be one of 'loose', 'crop' or 'wrap'.
+## 'loose' allows the image to grow to accomodate the new transformed image.
+## 'crop' keeps the same size as the original, clipping any part of the image
+## that is moved outside the bounding box.
+## 'wrap' keeps the same size as the original, but does not clip the part of the
+## image that is outside the bounding box. Instead, it wraps it back into the
+## image.
+##
+## If called with only 3 arguments, @var{bbox} is set to 'loose' by default.
+## @end deftypefn
+
+## Author: Jeff Orchard <jjo@sfu.ca>
+## Created: June 2002
+
+function g = imshear(m, axis, alpha, bbox, noshift)
+
+	# The code below only does y-shearing. This is because of
+	# the implementation of fft (operates on columns, but not rows).
+	# So, transpose first for x-shearing.
+	if ( strcmp(axis, "x")==1 )
+		m = m';
+	endif
+
+	if ( nargin < 4 )
+		bbox = "loose";
+		noshift = 0;
+	elseif ( nargin < 5 )
+		noshift = 0;
+	endif
+
+	[ydim_orig xdim_orig] = size(m);
+	if ( strcmp(bbox, "wrap") == 0 )
+		ypad = ceil( (xdim_orig+1)/2 * abs(alpha) );
+		m = impad(m, [0,0], [ypad,ypad]);
+	endif
+
+	[ydim_new xdim_new] = size(m);
+	xcentre = ( xdim_new + 1 ) / 2;
+	ycentre = ( ydim_new + 1 ) / 2;
+
+	# This applies FFT to columns of m (x-axis remains a spatial axis).
+	# Because the way that fft and fftshift are implemented, the origin
+	# will move by 1/2 pixel, depending on the polarity of the image
+	# dimensions.
+	#
+	# If dim is even (=2n), then the origin of the fft below is located
+	# at the centre of pixel (n+1). ie. if dim=16, then centre is at 9.
+	#
+	# If dim is odd (=2n+1), then the origin of the fft below is located
+	# at the centre of pixel (n). ie. if dim=15, then centre is at 8.
+	if ( noshift==1 )
+		M = fft(m);
+	else
+		#M = imtranslate(fft(imtranslate(m, -xcentre, ycentre, "wrap")), xcentre, -ycentre, "wrap");
+		M = fftshift(fft(fftshift(m)));
+	endif
+
+	[ydim xdim] = size(m);
+	x = zeros(ydim, xdim);
+
+	# Find coords of the origin of the image.
+	if ( noshift==1 )
+		xc_coord = 1;
+		yc_coord = 1;
+		l = (1:ydim)' - yc_coord;
+		r = (1:xdim) - xc_coord;
+		if ( strcmp(bbox, "wrap")==1 )
+			l((ydim/2):ydim) = l((ydim/2):ydim) - ydim;
+			r((xdim/2):xdim) = r((xdim/2):xdim) - xdim;
+		endif
+	else
+		xc_coord = (xdim+1)/2;
+		yc_coord = (ydim+1)/2;
+		l = (1:ydim)' - yc_coord;
+		r = (1:xdim) - xc_coord;
+	endif
+	x = l * r;
+
+	Ms = M.* exp(2*pi*I*alpha/ydim * x);
+
+	if ( noshift==1 )
+		g = abs(ifft(Ms));
+	else
+		#g = abs(imtranslate( ifft( imtranslate(Ms, -xcentre, ycentre, "wrap") ), xcentre, -ycentre, "wrap"));
+		g = abs( fftshift(ifft(ifftshift(Ms))) );
+	endif
+
+	if ( strcmp(bbox, "crop")==1 )
+		g = g(ypad+1:ydim_orig+ypad, :);
+	endif
+
+	# Un-transpose if x-shearing was wanted
+	if ( strcmp(axis, "x")==1 )
+		g = g';
+	endif
+endfunction

new file mode 100644
--- /dev/null
+++ b/imtranslate.m
@@ -0,0 +1,51 @@
+# -*- texinfo -*-
+# @deftypefn {Function File} {@var{Y}} = imtranslate (@var{M}, @var{x}, @var{y} [, @var{bbox}])
+# Translate a 2D image by (x,y) using Fourier interpolation.
+#
+# @var{M} is a matrix, and is translated to the right by @var{X} pixels
+# and translated up by @var{Y} pixels.
+#
+# @var{bbox} can be either 'crop' or 'wrap' (default).
+#
+# @end deftypefn
+
+function Y = imtranslate(X, a, b, bbox_in)
+
+	bbox = "wrap";
+	if ( nargin > 3 )
+		bbox = bbox_in;
+	endif
+
+	if ( strcmp(bbox, "crop")==1 )
+
+		xpad = [0,0];
+		if (a>0)
+			xpad = [0,ceil(a)];
+		elseif (a<0)
+			xpad = [-ceil(a),0];
+		endif
+
+		ypad = [0,0];
+		if (b>0)
+			ypad = [ceil(b),0];
+		elseif (b<0)
+			ypad = [0,-ceil(b)];
+		endif
+
+		X = impad(X, xpad, ypad, 'zeros')
+	endif
+
+
+	[dimy, dimx] = size(X);
+	x = ifftshift(fft2(fftshift(X)));
+	px = exp(-2*pi*i*a*(0:dimx-1)/dimx);
+	py = exp(-2*pi*i*b*(0:dimy-1)/dimy)';
+	P = py * px;
+	y = x .* P;
+	Y = abs( ifftshift(ifft2(fftshift(y))) );
+	#Y = ifftshift(ifft2(fftshift(y)));
+
+	if ( strcmp(bbox, "crop")==1 )
+		Y = Y(  ypad(1)+1:dimy-ypad(2) , xpad(1)+1:dimx-xpad(2));
+	endif
+endfunction