octave_packages/image-1.0.15/imrotate.m

   1 ## Copyright (C) 2004-2005 Justus H. Piater
   2 ##
   3 ## This program is free software; you can redistribute it and/or
   4 ## modify it under the terms of the GNU General Public License
   5 ## as published by the Free Software Foundation; either version 2
   6 ## of the License, or (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but
   9 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 ## General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {} imrotate(@var{imgPre}, @var{theta}, @var{method}, @var{bbox}, @var{extrapval})
  18 ## Rotation of a 2D matrix about its center.
  19 ##
  20 ## Input parameters:
  21 ##
  22 ##   @var{imgPre}   a gray-level image matrix
  23 ##
  24 ##   @var{theta}    the rotation angle in degrees counterclockwise
  25 ##
  26 ##   @var{method}
  27 ##     @itemize @w
  28 ##       @item "nearest" neighbor: fast, but produces aliasing effects (default).
  29 ##       @item "bilinear" interpolation: does anti-aliasing, but is slightly slower.
  30 ##       @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.
  31 ##       @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.
  32 ##     @end itemize
  33 ##
  34 ##   @var{bbox}
  35 ##     @itemize @w
  36 ##       @item "loose" grows the image to accommodate the rotated image (default).
  37 ##       @item "crop" rotates the image about its center, clipping any part of the image that is moved outside its boundaries.
  38 ##     @end itemize
  39 ##
  40 ##   @var{extrapval} sets the value used for extrapolation. The default value
  41 ##      is @code{NA} for images represented using doubles, and 0 otherwise.
  42 ##      This argument is ignored of Fourier interpolation is used.
  43 ##
  44 ## Output parameters:
  45 ##
  46 ##   @var{imgPost}  the rotated image matrix
  47 ##
  48 ##   @var{H}        the homography mapping original to rotated pixel
  49 ##                   coordinates. To map a coordinate vector c = [x;y] to its
  50 ##           rotated location, compute round((@var{H} * [c; 1])(1:2)).
  51 ##
  52 ##   @var{valid}    a binary matrix describing which pixels are valid,
  53 ##                  and which pixels are extrapolated. This output is
  54 ##                  not available if Fourier interpolation is used.
  55 ## @end deftypefn
  56
  57 ## Author: Justus H. Piater  <Justus.Piater@ULg.ac.be>
  58 ## Created: 2004-10-18
  59 ## Version: 0.7
  60
  61 function [imgPost, H, valid] = imrotate(imgPre, thetaDeg, interp="nearest", bbox="loose", extrapval=NA)
  62   ## Check input
  63   if (nargin < 2)
  64     error("imrotate: not enough input arguments");
  65   endif
  66   [imrows, imcols, imchannels, tmp] = size(imgPre);
  67   if (tmp != 1 || (imchannels != 1 && imchannels != 3))
  68     error("imrotate: first input argument must be an image");
  69   endif
  70   if (!isscalar(thetaDeg))
  71     error("imrotate: the angle must be given as a scalar");
  72   endif
  73   if (!any(strcmpi(interp, {"nearest", "linear", "bilinear", "cubic", "bicubic", "Fourier"})))
  74     error("imrotate: unsupported interpolation method");
  75   endif
  76   if (any(strcmpi(interp, {"bilinear", "bicubic"})))
  77     interp = interp(3:end); # Remove "bi"
  78   endif
  79   if (!any(strcmpi(bbox, {"loose", "crop"})))
  80     error("imrotate: bounding box must be either 'loose' or 'crop'");
  81   endif
  82   if (!isscalar(extrapval))
  83     error("imrotate: extrapolation value must be a scalar");
  84   endif
  85
  86   ## Input checking done. Start working
  87   thetaDeg = mod(thetaDeg, 360); # some code below relies on positive angles
  88   theta = thetaDeg * pi/180;
  89
  90   sizePre = size(imgPre);
  91
  92   ## We think in x,y coordinates here (rather than row,column), except
  93   ## for size... variables that follow the usual size() convention. The
  94   ## coordinate system is aligned with the pixel centers.
  95
  96   R = [cos(theta) sin(theta); -sin(theta) cos(theta)];
  97
  98   if (nargin >= 4 && strcmp(bbox, "crop"))
  99     sizePost = sizePre;
 100   else
 101     ## Compute new size by projecting zero-base image corner pixel
 102     ## coordinates through the rotation:
 103     corners = [0, 0;
 104                (R * [sizePre(2) - 1; 0             ])';
 105                (R * [sizePre(2) - 1; sizePre(1) - 1])';
 106                (R * [0             ; sizePre(1) - 1])' ];
 107     sizePost(2) = round(max(corners(:,1)) - min(corners(:,1))) + 1;
 108     sizePost(1) = round(max(corners(:,2)) - min(corners(:,2))) + 1;
 109     ## This size computation yields perfect results for 0-degree (mod
 110     ## 90) rotations and, together with the computation of the center of
 111     ## rotation below, yields an image whose corresponding region is
 112     ## identical to "crop". However, we may lose a boundary of a
 113     ## fractional pixel for general angles.
 114   endif
 115
 116   ## Compute the center of rotation and the translational part of the
 117   ## homography:
 118   oPre  = ([ sizePre(2);  sizePre(1)] + 1) / 2;
 119   oPost = ([sizePost(2); sizePost(1)] + 1) / 2;
 120   T = oPost - R * oPre;         # translation part of the homography
 121
 122   ## And here is the homography mapping old to new coordinates:
 123   H = [[R; 0 0] [T; 1]];
 124
 125   ## Treat trivial rotations specially (multiples of 90 degrees):
 126   if (mod(thetaDeg, 90) == 0)
 127     nRot90 = mod(thetaDeg, 360) / 90;
 128     if (mod(thetaDeg, 180) == 0 || sizePre(1) == sizePre(2) ||
 129         strcmpi(bbox, "loose"))
 130       imgPost = rot90(imgPre, nRot90);
 131       return;
 132     elseif (mod(sizePre(1), 2) == mod(sizePre(2), 2))
 133       ## Here, bbox is "crop" and the rotation angle is +/- 90 degrees.
 134       ## This works only if the image dimensions are of equal parity.
 135       imgRot = rot90(imgPre, nRot90);
 136       imgPost = zeros(sizePre);
 137       hw = min(sizePre) / 2 - 0.5;
 138       imgPost   (round(oPost(2) - hw) : round(oPost(2) + hw),
 139                  round(oPost(1) - hw) : round(oPost(1) + hw) ) = ...
 140           imgRot(round(oPost(1) - hw) : round(oPost(1) + hw),
 141                  round(oPost(2) - hw) : round(oPost(2) + hw) );
 142       return;
 143     else
 144       ## Here, bbox is "crop", the rotation angle is +/- 90 degrees, and
 145       ## the image dimensions are of unequal parity. This case cannot
 146       ## correctly be handled by rot90() because the image square to be
 147       ## cropped does not align with the pixels - we must interpolate. A
 148       ## caller who wants to avoid this should ensure that the image
 149       ## dimensions are of equal parity.
 150     endif
 151   end
 152
 153   ## Now the actual rotations happen
 154   if (strcmpi(interp, "Fourier"))
 155     c = class (imgPre);
 156     imgPre = im2double (imgPre);
 157     if (isgray(imgPre))
 158       imgPost = imrotate_Fourier(imgPre, thetaDeg, interp, bbox);
 159     else # rgb image
 160       for i = 3:-1:1
 161         imgPost(:,:,i) = imrotate_Fourier(imgPre(:,:,i), thetaDeg, interp, bbox);
 162       endfor
 163     endif
 164     valid = NA;
 165
 166     switch (c)
 167       case "uint8"
 168         imgPost = im2uint8 (imgPost);
 169       case "uint16"
 170         imgPost = im2uint16 (imgPost);
 171       case "single"
 172         imgPost = single (imgPost);
 173     endswitch
 174   else
 175     [imgPost, valid] = imperspectivewarp(imgPre, H, interp, bbox, extrapval);
 176   endif
 177 endfunction
 178
 179 %!test
 180 %! ## Verify minimal loss across six rotations that add up to 360 +/- 1 deg.:
 181 %! methods = { "nearest", "bilinear", "bicubic", "Fourier" };
 182 %! angles     = [ 59  60  61  ];
 183 %! tolerances = [ 7.4 8.5 8.6     # nearest
 184 %!                3.5 3.1 3.5     # bilinear
 185 %!                2.7 2.0 2.7     # bicubic
 186 %!                2.7 1.6 2.8 ]/8;  # Fourier
 187 %!
 188 %! # This is peaks(50) without the dependency on the plot package
 189 %! x = y = linspace(-3,3,50);
 190 %! [X,Y] = meshgrid(x,y);
 191 %! x = 3*(1-X).^2.*exp(-X.^2 - (Y+1).^2) \
 192 %!      - 10*(X/5 - X.^3 - Y.^5).*exp(-X.^2-Y.^2) \
 193 %!      - 1/3*exp(-(X+1).^2 - Y.^2);
 194 %!
 195 %! x -= min(x(:));            # Fourier does not handle neg. values well
 196 %! x = x./max(x(:));
 197 %! for m = 1:(length(methods))
 198 %!   y = x;
 199 %!   for i = 1:5
 200 %!     y = imrotate(y, 60, methods{m}, "crop", 0);
 201 %!   end
 202 %!   for a = 1:(length(angles))
 203 %!     assert(norm((x - imrotate(y, angles(a), methods{m}, "crop", 0))
 204 %!                 (10:40, 10:40)) < tolerances(m,a));
 205 %!   end
 206 %! end
 207
 208
 209 %!#test
 210 %! ## Verify exactness of near-90 and 90-degree rotations:
 211 %! X = rand(99);
 212 %! for angle = [90 180 270]
 213 %!   for da = [-0.1 0.1]
 214 %!     Y = imrotate(X,   angle + da , "nearest", :, 0);
 215 %!     Z = imrotate(Y, -(angle + da), "nearest", :, 0);
 216 %!     assert(norm(X - Z) == 0); # exact zero-sum rotation
 217 %!     assert(norm(Y - imrotate(X, angle, "nearest", :, 0)) == 0); # near zero-sum
 218 %!   end
 219 %! end
 220
 221
 222 %!#test
 223 %! ## Verify preserved pixel density:
 224 %! methods = { "nearest", "bilinear", "bicubic", "Fourier" };
 225 %! ## This test does not seem to do justice to the Fourier method...:
 226 %! tolerances = [ 4 2.2 2.0 209 ];
 227 %! range = 3:9:100;
 228 %! for m = 1:(length(methods))
 229 %!   t = [];
 230 %!   for n = range
 231 %!     t(end + 1) = sum(imrotate(eye(n), 20, methods{m}, :, 0)(:));
 232 %!   end
 233 %!   assert(t, range, tolerances(m));
 234 %! end
 235