Mercurial > hg > octave-image

changeset 137:abe4962005a4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[for Justus Piater] fix half-pixel alignment problem and add optimization for n*90 degree rotations.
author pkienzle
date Wed, 16 Feb 2005 04:07:59 +0000
parents 78b9878c50f0
children 4d4f17064410
files imrotate.m imrotate_Fourier.m
diffstat 2 files changed, 374 insertions(+), 47 deletions(-) [+]
line wrap: on
line diff
--- a/imrotate.m
+++ b/imrotate.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2004 Justus Piater
+## Copyright (C) 2004-2005 Justus Piater
 ##
 ## This program is free software; you can redistribute it and/or
 ## modify it under the terms of the GNU General Public License
@@ -21,19 +21,27 @@
 ##
 ## Input parameters:
 ##
-##   @var{imgPre}   an input image matrix
+##   @var{imgPre}   a gray-level image matrix
 ##
 ##   @var{theta}    the rotation angle in degrees counterclockwise
 ##
-##   @var{method}   "nearest" neighbor (default; faster, but produces
-##                  aliasing effects) or "bilinear" interpolation
-##                     (preferred; does anti-aliasing)
+##   @var{method}
+##     @itemize @w
+##       @item "nearest" neighbor: fast, but produces aliasing effects (default).
+##       @item "bilinear" interpolation: does anti-aliasing, but is slightly slower.
+##       @item "bicubic" interpolation: does anti-aliasing, preserves edges better than bilinear interpolation, but gray levels may slightly overshoot at sharp edges. This is probably the best method for most purposes, but also the slowest.
+##       @item "Fourier" uses Fourier interpolation, decomposing the rotation matrix into 3 shears. This method often results in different artifacts than homography-based methods.  Instead of slightly blurry edges, this method can result in ringing artifacts (little waves near high-contrast edges).  However, Fourier interpolation is better at maintaining the image information, so that unrotating will result in an image closer to the original than the other methods.
+##     @end itemize
 ##
-##   @var{bbox}     "loose" (default) or "crop"
+##   @var{bbox}
+##     @itemize @w
+##       @item "loose" grows the image to accommodate the rotated image (default).
+##       @item "crop" rotates the image about its center, clipping any part of the image that is moved outside its boundaries.
+##     @end itemize
 ##
 ## Output parameters:
 ##
-##   @var{imgPos}t  the rotated image matrix
+##   @var{imgPost}  the rotated image matrix
 ##
 ##   @var{H}        the homography mapping original to rotated pixel
 ##                   coordinates. To map a coordinate vector c = [x;y] to its
@@ -42,71 +50,138 @@
 
 ## Author: Justus H. Piater  <Justus.Piater@ULg.ac.be>
 ## Created: 2004-10-18
-## Version: 0.1
+## Version: 0.7
 
-function [imgPost, H] = imrotate(imgPre, theta, method, bbox)
-  if (nargin < 2)
-    help("imrotate");
-    return;
+function [imgPost, H] = imrotate(imgPre, thetaDeg, method, bbox)
+  if (nargin < 4)
+    bbox = "loose";
+    if (nargin < 3)
+      method = "nearest";
+      if (nargin < 2)
+	usage("imrotate(img, angle [, method [, bbox]]");
+      endif
+    endif
   endif
 
-  theta = theta * pi/180;
+  thetaDeg = mod(thetaDeg, 360); # some code below relies on positive angles
+  theta = thetaDeg * pi/180;
 
-  sizePre  = size(imgPre);
+  sizePre = size(imgPre);
 
-  ## We think in x,y coordinates here (rather than row,column).
+  ## We think in x,y coordinates here (rather than row,column), except
+  ## for size... variables that follow the usual size() convention. The
+  ## coordinate system is aligned with the pixel centers.
 
   R = [cos(theta) sin(theta); -sin(theta) cos(theta)];
 
-  if (nargin == 4 && strcmp(bbox, "crop"))
+  if (nargin >= 4 && strcmp(bbox, "crop"))
     sizePost = sizePre;
   else
-    ## Compute new size by projecting image corners through the rotation:
+    ## Compute new size by projecting zero-base image corner pixel
+    ## coordinates through the rotation:
     corners = [0, 0;
-	       (R * [sizePre(2); 1])';
-	       (R * [sizePre(2); sizePre(1)])';
-	       (R * [1; sizePre(1)])'          ];
-    sizePost(2) = ceil(max(corners(:,1))) - floor(min(corners(:,1)));
-    sizePost(1) = ceil(max(corners(:,2))) - floor(min(corners(:,2)));
+	       (R * [sizePre(2) - 1; 0             ])';
+	       (R * [sizePre(2) - 1; sizePre(1) - 1])';
+	       (R * [0             ; sizePre(1) - 1])' ];
+    sizePost(2) = round(max(corners(:,1)) - min(corners(:,1))) + 1;
+    sizePost(1) = round(max(corners(:,2)) - min(corners(:,2))) + 1;
+    ## This size computation yields perfect results for 0-degree (mod
+    ## 90) rotations and, together with the computation of the center of
+    ## rotation below, yields an image whose corresponding region is
+    ## identical to "crop". However, we may lose a boundary of a
+    ## fractional pixel for general angles.
   endif
 
-  ## Compute the translation part of the homography:
-  oPre  = ([sizePre(2) ; sizePre(1) ] + 1) / 2;
+  ## Compute the center of rotation and the translational part of the
+  ## homography:
+  oPre  = ([ sizePre(2);  sizePre(1)] + 1) / 2;
   oPost = ([sizePost(2); sizePost(1)] + 1) / 2;
-  T = oPost - R * oPre;
+  T = oPost - R * oPre;		# translation part of the homography
 
   ## And here is the homography mapping old to new coordinates:
   H = [[R; 0 0] [T; 1]];
 
-  Hinv = inv(H);
-
-  ## "Pre"  variables hold pre -rotation values;
-  ## "Post" variables hold post-rotation values.
-
-  ## Target coordinates:
-  [xPost, yPost] = meshgrid(1:(sizePost(2)), 1:(sizePost(1)));
+  ## Treat trivial rotations specially (multiples of 90 degrees):
+  if (mod(thetaDeg, 90) == 0)
+    nRot90 = mod(thetaDeg, 360) / 90;
+    if (mod(thetaDeg, 180) == 0 || sizePre(1) == sizePre(2) ||
+	strcmp(bbox, "loose"))
+      imgPost = rot90(imgPre, nRot90);
+      return;
+    elseif (mod(sizePre(1), 2) == mod(sizePre(2), 2))
+      ## Here, bbox is "crop" and the rotation angle is +/- 90 degrees.
+      ## This works only if the image dimensions are of equal parity.
+      imgRot = rot90(imgPre, nRot90);
+      imgPost = zeros(sizePre);
+      hw = min(sizePre) / 2 - 0.5;
+      imgPost   (round(oPost(2) - hw) : round(oPost(2) + hw),
+		 round(oPost(1) - hw) : round(oPost(1) + hw) ) = ...
+	  imgRot(round(oPost(1) - hw) : round(oPost(1) + hw),
+		 round(oPost(2) - hw) : round(oPost(2) + hw) );
+      return;
+    else
+      ## Here, bbox is "crop", the rotation angle is +/- 90 degrees, and
+      ## the image dimensions are of unequal parity. This case cannot
+      ## correctly be handled by rot90() because the image square to be
+      ## cropped does not align with the pixels - we must interpolate. A
+      ## caller who wants to avoid this should ensure that the image
+      ## dimensions are of equal parity.
+    endif
+  end
 
-  ## Compute corresponding source coordinates:
-  xPre = Hinv(1,1) * xPost + Hinv(1,2) * yPost + Hinv(1,3);
-  yPre = Hinv(2,1) * xPost + Hinv(2,2) * yPost + Hinv(2,3);
-  ## zPre is guaranteed to be 1, since the last row of H (and thus of
-  ## Hinv) is [0 0 1].
+  ## For better readability of this spaghetti implementation, I keep the
+  ## branches pertaining to the various 'method's all at the first
+  ## level, even though this causes a slight redundancy in the if
+  ## statements.
+
+  imgPost = [];
+
+  if (strcmp(method, "Fourier"))
+    imgPost = imrotate_Fourier(imgPre, thetaDeg, method, bbox);
+  else
+    ## This section pertains to all non-Fourier methods.
+
+    ## "Pre"  variables hold pre -rotation values;
+    ## "Post" variables hold post-rotation values.
 
-  ## Now map the image, either by nearest neighbor or by bilinear
-  ## interpolation:
-  if (nargin < 3 || !size(method) || strcmp(method, "nearest"))
+    ## General rotation: map pixel coordinates back from the Post to the
+    ## Pre img
+    Hinv = inv(H);
+
+    ## Target coordinates:
+    [xPost, yPost] = meshgrid(1:(sizePost(2)), 1:(sizePost(1)));
+
+    ## Compute corresponding source coordinates:
+    xPre = Hinv(1,1) * xPost + Hinv(1,2) * yPost + Hinv(1,3);
+    yPre = Hinv(2,1) * xPost + Hinv(2,2) * yPost + Hinv(2,3);
+    ## zPre is guaranteed to be 1, since the last row of H (and thus of
+    ## Hinv) is [0 0 1].
+  endif
+
+  ## Now map the image using the coordinates computed in the else branch above:
+  if (strcmp(method, "nearest"))
     ## nearest-neighbor: simply round Pre coordinates
     xPre = round(xPre);
     yPre = round(yPre);
     valid = find(1 <= xPre & xPre <= sizePre(2) &
 		 1 <= yPre & yPre <= sizePre(1)  );
+    if (!length(valid))
+      warning("input image too small");
+      imgPost = 0;
+      return;
+    endif
+    
     iPre  = sub2ind(sizePre , yPre (valid), xPre (valid));
     iPost = sub2ind(sizePost, yPost(valid), xPost(valid));
 
     imgPost = zeros(sizePost);
     imgPost(iPost) = imgPre(iPre);
-  else
-    ## bilinear interpolation between the four floor and ceiling coordinates
+  elseif(!strcmp(method, "Fourier"))
+    ## This section pertains to "bilinear" and "bicubic" methods.
+
+    ## With interpolation, one unavoidably loses up to one or two pixel
+    ## rows or columns at the image boundaries.
+
     xPreFloor = floor(xPre);
     xPreCeil  = ceil (xPre);
     yPreFloor = floor(yPre);
@@ -114,6 +189,11 @@
 
     valid = find(1 <= xPreFloor & xPreCeil <= sizePre(2) &
 		 1 <= yPreFloor & yPreCeil <= sizePre(1)  );
+    if (!length(valid))
+      warning("input image too small");
+      imgPost = 0;
+      return;
+    endif
 
     xPreFloor = xPreFloor(valid);
     xPreCeil  = xPreCeil (valid);
@@ -131,12 +211,89 @@
     yPreFrac = yPre(valid) - yPreFloor;
 
     iPost = sub2ind(sizePost, yPost(valid), xPost(valid));
+  endif
 
+  if (strcmp(method, "bilinear"))
     imgPost = zeros(sizePost);
-    imgPost(iPost) = ...
-	round(imgPre(iPreFF) .* (1 - xPreFrac) .* (1 - yPreFrac) + ...
-	      imgPre(iPreCF) .*      xPreFrac  .* (1 - yPreFrac) + ...
-	      imgPre(iPreCC) .*      xPreFrac  .*      yPreFrac  + ...
-	      imgPre(iPreFC) .* (1 - xPreFrac) .*      yPreFrac       );
+    ## bilinear interpolation between the four floor and ceiling coordinates
+    imgPost(iPost) = (imgPre(iPreFF) .* (1 - xPreFrac) .* (1 - yPreFrac) +
+		      imgPre(iPreCF) .*      xPreFrac  .* (1 - yPreFrac) +
+		      imgPre(iPreCC) .*      xPreFrac  .*      yPreFrac  +
+		      imgPre(iPreFC) .* (1 - xPreFrac) .*      yPreFrac   );
+  elseif (strcmp(method, "bicubic"))
+    ## bicubic interpolation (see Numerical Recipes)
+    ## This code, together with the prerequisites above, is not limited
+    ## to this particular use but applies to generic bicubic
+    ## interpolation in the following scenario:
+    ## - source data are stored in a matrix,
+    ## - interpolated coordinates may lie anywhere, no regularity is assumed.
+
+    ## precompute the required derivatives at the source image pixels:
+    imgPreDx  = conv2(imgPre  , [ 0.5 0 -0.5] , "same");
+    imgPreDy  = conv2(imgPre  , [-0.5 0  0.5]', "same");
+    imgPreDxy = conv2(imgPreDx, [-0.5 0  0.5]', "same");
+
+    ## Interpolation is done on a square of pixels and their
+    ## derivatives along x, y, and xy. The square is indexed as:
+    ##   43                        FF CF
+    ##   12  which corresponds to  FC CC
+
+    ## Coefficient matrix W
+    ## C11 12         21            31            41                p deriv
+    W = [1  0 -3  2    0  0  0  0   -3  0  9 -6    2  0 -6  4    ## 1
+	 0  0  0  0    0  0  0  0    3  0 -9  6   -2  0  6 -4	 ## 2
+	 0  0  0  0    0  0  0  0    0  0  9 -6    0  0 -6  4	 ## 3
+	 0  0  3 -2    0  0  0  0    0  0 -9  6    0  0  6 -4	 ## 4
+	 
+	 0  0  0  0    1  0 -3  2   -2  0  6 -4    1  0 -3  2	 ## 1 x
+	 0  0  0  0    0  0  0  0   -1  0  3 -2    1  0 -3  2	 ## 2 x
+	 0  0  0  0    0  0  0  0    0  0 -3  2    0  0  3 -2	 ## 3 x
+	 0  0  0  0    0  0  3 -2    0  0 -6  4    0  0  3 -2	 ## 4 x
+	 
+	 0  1 -2  1    0  0  0  0    0 -3  6 -3    0  2 -4  2	 ## 1 y
+	 0  0  0  0    0  0  0  0    0  3 -6  3    0 -2  4 -2	 ## 2 y
+	 0  0  0  0    0  0  0  0    0  0 -3  3    0  0  2 -2	 ## 3 y
+	 0  0 -1  1    0  0  0  0    0  0  3 -3    0  0 -2  2	 ## 4 y
+	 
+	 0  0  0  0    0  1 -2  1    0 -2  4 -2    0  1 -2  1	 ## 1 xy
+	 0  0  0  0    0  0  0  0    0 -1  2 -1    0  1 -2  1	 ## 2 xy
+	 0  0  0  0    0  0  0  0    0  0  1 -1    0  0 -1  1	 ## 3 xy
+	 0  0  0  0    0  0 -1  1    0  0  2 -2    0  0 -1  1];  ## 4 xy
+
+    u = 1 - yPreFrac;
+    values = zeros(size(valid));
+    for ci = 4:-1:1
+      ## compute ci'th row of matrix C:
+
+      col = 4*(ci - 1) + 1;
+      c{1} = (W( 1,col) * imgPre   (iPreFC) + W( 2,col) * imgPre   (iPreCC) +
+	      W( 5,col) * imgPreDx (iPreFC) + W( 6,col) * imgPreDx (iPreCC)  );
+
+      col++;
+      c{2} = (W( 9,col) * imgPreDy (iPreFC) + W(10,col) * imgPreDy (iPreCC) +
+	      W(13,col) * imgPreDxy(iPreFC) + W(14,col) * imgPreDxy(iPreCC)  );
+
+      for cii = 3:4
+	col++;
+	c{cii} = ...
+	    (W( 1,col) * imgPre   (iPreFC) + W( 2,col) * imgPre   (iPreCC) +
+	     W( 3,col) * imgPre   (iPreCF) + W( 4,col) * imgPre   (iPreFF) +
+	     W( 5,col) * imgPreDx (iPreFC) + W( 6,col) * imgPreDx (iPreCC) +
+	     W( 7,col) * imgPreDx (iPreCF) + W( 8,col) * imgPreDx (iPreFF) +
+	     W( 9,col) * imgPreDy (iPreFC) + W(10,col) * imgPreDy (iPreCC) +
+	     W(11,col) * imgPreDy (iPreCF) + W(12,col) * imgPreDy (iPreFF) +
+	     W(13,col) * imgPreDxy(iPreFC) + W(14,col) * imgPreDxy(iPreCC) +
+	     W(15,col) * imgPreDxy(iPreCF) + W(16,col) * imgPreDxy(iPreFF)  );
+      endfor
+
+      values .*= xPreFrac;
+      values  += ((c{4} .* u + c{3}) .* u + c{2}) .* u + c{1};
+    endfor
+    imgPost = zeros(sizePost);
+    imgPost(iPost) = values;
+  endif
+
+  if (!prod(size(imgPost)))
+    error(sprintf("Interpolation method %s not implemented", method));
   endif
 endfunction
new file mode 100644
--- /dev/null
+++ b/imrotate_Fourier.m
@@ -0,0 +1,170 @@
+## Copyright (C) 2002 Jeff Orchard <jorchard@cs.uwaterloo.ca>
+##
+## This program is free software; you can redistribute it and/or
+## modify it under the terms of the GNU General Public License
+## as published by the Free Software Foundation; either version 2
+## of the License, or (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, write to the Free Software
+## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+## -*- 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 <jorchard@cs.uwaterloo.ca>
+## Created: Oct. 14, 2002
+
+function fs = imrotate_Fourier(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);
+
+	# This finds the index coords of the centre of the image (indices are base-1)
+	#   eg. if xdim_orig=8, then xcentre_orig=4.5 (half-way between 1 and 8)
+	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), 'loose');
+		f2 = imshear(f1, 'y', S2(2,1), 'loose');
+		fs = real( imshear(f2, 'x', S1(1,2), 'loose') );
+		%fs = f2;
+		xcentre_new = (size(fs,2)+1) / 2;
+		ycentre_new = (size(fs,1)+1) / 2;
+	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) );
+		fs = fs( (ymin_loose+1):(ymax_loose-1) , (xmin_loose+1):(xmax_loose-1) );
+
+	endif
+
+
+endfunction