octave_packages/image-1.0.15/imperspectivewarp.m

   1 ## Copyright (C) 2006  Søren Hauberg
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2, or (at your option)
   6 ## 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 file.  If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} @var{warped} = imperspectivewarp(@var{im}, @var{P}, @var{interp}, @var{bbox}, @var{extrapval})
  18 ## @deftypefnx{Function File} [@var{warped}, @var{valid}] = imperspectivewarp(...)
  19 ## Applies the spatial perspective homogeneous transformation @var{P} to the image @var{im}.
  20 ## The transformation matrix @var{P} must be a 3x3 homogeneous matrix, or 2x2 or 2x3
  21 ## affine transformation matrix.
  22 ##
  23 ## The resulting image @var{warped} is computed using an interpolation method that
  24 ## can be selected through the @var{interp} argument. This must be one
  25 ## of the following strings
  26 ## @table @code
  27 ## @item "nearest"
  28 ## Nearest neighbor interpolation.
  29 ## @item "linear"
  30 ## @itemx "bilinear"
  31 ## Bilinear interpolation. This is the default behavior.
  32 ## @item "cubic"
  33 ## @itemx "bicubic"
  34 ## Bicubic interpolation.
  35 ## @end table
  36 ##
  37 ## By default the resulting image contains the entire warped image. In some situation
  38 ## you only parts of the warped image. The argument @var{bbox} controls this, and can
  39 ## be one of the following strings
  40 ## @table @code
  41 ## @item "loose"
  42 ## The entire warped result is returned. This is the default behavior.
  43 ## @item "crop"
  44 ## The central part of the image of the same size as the input image is returned.
  45 ## @item "same"
  46 ## The size and coordinate system of the input image is keept.
  47 ## @end table
  48 ##
  49 ## All values of the result that fall outside the original image will
  50 ## be set to @var{extrapval}. For images of class @code{double} @var{extrapval}
  51 ## defaults to @code{NA} and for other classes it defaults to 0.
  52 ##
  53 ## The optional output @var{valid} is a matrix of the same size as @var{warped}
  54 ## that contains the value 1 in pixels where @var{warped} contains an interpolated
  55 ## value, and 0 in pixels where @var{warped} contains an extrapolated value.
  56 ## @seealso{imremap, imrotate, imresize, imshear, interp2}
  57 ## @end deftypefn
  58
  59 function [warped, valid] = imperspectivewarp(im, P, interp = "bilinear", bbox = "loose", extrapolation_value = NA)
  60   ## Check input
  61   if (nargin < 2)
  62     print_usage();
  63   endif
  64
  65   [imrows, imcols, imchannels, tmp] = size(im);
  66   if (tmp != 1 || (imchannels != 1 && imchannels != 3))
  67     error("imperspectivewarp: first input argument must be an image");
  68   endif
  69
  70   if (ismatrix(P) && ndims(P) == 2)
  71     if (issquare(P) && rows(P) == 3) # 3x3 matrix
  72       if (P(3,3) != 0)
  73         P /= P(3,3);
  74       else
  75         error("imperspectivewarp: P(3,3) must be non-zero");
  76       endif
  77     elseif (rows(P) == 2 && (columns(P) == 2 || columns(P) == 3)) # 2x2 or 2x3 matrix
  78       P(3,3) = 1;
  79     else # unsupported matrix size
  80       error("imperspectivewarp: transformation matrix must be 2x2, 2x3, or 3x3");
  81     endif
  82   else
  83     error("imperspectivewarp: transformation matrix not valid");
  84   endif
  85
  86   if (!any(strcmpi(interp, {"nearest", "linear", "bilinear", "cubic", "bicubic"})))
  87     error("imperspectivewarp: unsupported interpolation method");
  88   endif
  89   if (any(strcmpi(interp, {"bilinear", "bicubic"})))
  90     interp = interp(3:end); # Remove "bi"
  91   endif
  92   interp = lower(interp);
  93
  94   if (!any(strcmpi(bbox, {"loose", "crop", "same"})))
  95     error("imperspectivewarp: bounding box must be either 'loose', 'crop', or 'same'");
  96   endif
  97
  98   if (!isscalar(extrapolation_value))
  99     error("imperspective: extrapolation value must be a scalar");
 100   endif
 101
 102   ## Do the transformation
 103   [y, x, tmp] = size(im);
 104   ## Transform corners
 105   corners = [1, 1, 1;
 106              1, y, 1;
 107              x, 1, 1;
 108              x, y, 1]';
 109   Tcorners = P*corners;
 110   Tx = Tcorners(1,:)./Tcorners(3,:);
 111   Ty = Tcorners(2,:)./Tcorners(3,:);
 112
 113   ## Do cropping?
 114   x1 = round(min(Tx)); x2 = round(max(Tx));
 115   y1 = round(min(Ty)); y2 = round(max(Ty));
 116   # FIXME: This seems to work fine for rotations, but
 117   # somebody who knows computational geometry should
 118   # be able to come up with a better algorithm.
 119   if (strcmpi(bbox, "crop"))
 120     xl = x2 - x1 + 1;
 121     yl = y2 - y1 + 1;
 122     xd = (xl - x)/2;
 123     yd = (yl - y)/2;
 124     x1 += xd; x2 -= xd;
 125     y1 += yd; y2 -= yd;
 126   elseif (strcmpi(bbox, "same"))
 127     x1 = 1; x2 = x;
 128     y1 = 1; y2 = y;
 129   endif
 130
 131   ## Transform coordinates
 132   [X, Y] = meshgrid(x1:x2, y1:y2);
 133   [sy, sx] = size(X);
 134   D = [X(:), Y(:), ones(sx*sy, 1)]';
 135   PD = inv(P)*D;
 136   XI = PD(1,:)./PD(3,:);
 137   YI = PD(2,:)./PD(3,:);
 138   XI = reshape(XI, sy, sx);
 139   YI = reshape(YI, sy, sx);
 140
 141   clear X Y D PD;
 142
 143   ## Interpolate
 144   [warped, valid] = imremap(im, XI, YI, interp, extrapolation_value);
 145
 146 endfunction
 147
 148 %!demo
 149 %! ## Generate a synthetic image and show it
 150 %! I = tril(ones(100)) + abs(rand(100)); I(I>1) = 1;
 151 %! I(20:30, 20:30) = !I(20:30, 20:30);
 152 %! I(70:80, 70:80) = !I(70:80, 70:80);
 153 %! figure(), imshow(I);
 154 %! ## Resize the image to the double size and show it
 155 %! P = diag([1, 1, 0.5]);
 156 %! warped = imperspectivewarp(I, P);
 157 %! figure(), imshow(warped);
 158
 159 %!demo
 160 %! ## Generate a synthetic image and show it
 161 %! I = tril(ones(100)) + abs(rand(100)); I(I>1) = 1;
 162 %! I(20:30, 20:30) = !I(20:30, 20:30);
 163 %! I(70:80, 70:80) = !I(70:80, 70:80);
 164 %! figure(), imshow(I);
 165 %! ## Rotate the image around (0, 0) by -0.4 radians and show it
 166 %! R = [cos(-0.4) sin(-0.4); -sin(-0.4) cos(-0.4)];
 167 %! warped = imperspectivewarp(I, R, :, :, 0);
 168 %! figure(), imshow(warped);