octave_packages/image-1.0.15/imremap.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} = imremap(@var{im}, @var{XI}, @var{YI})
  18 ## @deftypefnx{Function File} @var{warped} = imremap(@var{im}, @var{XI}, @var{YI}, @var{interp}, @var{extrapval})
  19 ## @deftypefnx{Function File} [@var{warped}, @var{valid} ] = imremap(...)
  20 ## Applies any geometric transformation to the image @var{im}.
  21 ##
  22 ## The arguments @var{XI} and @var{YI} are lookup tables that define the resulting
  23 ## image
  24 ## @example
  25 ## @var{warped}(y,x) = @var{im}(@var{YI}(y,x), @var{XI}(y,x))
  26 ## @end example
  27 ## where @var{im} is assumed to be a continuous function, which is achieved
  28 ## by interpolation. Note that the image @var{im} is expressed in a (X, Y)-coordinate
  29 ## system and not a (row, column) system.
  30 ##
  31 ## The argument @var{interp} selects the used interpolation method, and most be one
  32 ## of the following strings
  33 ## @table @code
  34 ## @item "nearest"
  35 ## Nearest neighbor interpolation.
  36 ## @item "linear"
  37 ## @itemx "bilinear"
  38 ## Bilinear interpolation. This is the default behavior.
  39 ## @item "cubic"
  40 ## @itemx "bicubic"
  41 ## Bicubic interpolation.
  42 ## @end table
  43 ##
  44 ## All values of the result that fall outside the original image will
  45 ## be set to @var{extrapval}. For images of class @code{double} @var{extrapval}
  46 ## defaults to @code{NA} and for other classes it defaults to 0.
  47 ##
  48 ## The optional output @var{valid} is a matrix of the same size as @var{warped}
  49 ## that contains the value 1 in pixels where @var{warped} contains an interpolated
  50 ## value, and 0 in pixels where @var{warped} contains an extrapolated value.
  51 ## @seealso{imperspectivewarp, imrotate, imresize, imshear, interp2}
  52 ## @end deftypefn
  53
  54 function [warped, valid] = imremap(im, XI, YI, interp = "bilinear", extrapval = NA)
  55   ## Check input
  56   if (nargin < 3)
  57     print_usage();
  58   endif
  59
  60   [imrows, imcols, imchannels, tmp] = size(im);
  61   if (tmp != 1 || (imchannels != 1 && imchannels != 3))
  62     error("imremap: first input argument must be an image");
  63   endif
  64
  65   if (!size_equal(XI, YI) || !ismatrix(XI) || ndims(XI) != 2)
  66     error("imremap: XI and YI must be matrices of the same size");
  67   endif
  68
  69   if (!any(strcmpi(interp, {"nearest", "linear", "bilinear", "cubic", "bicubic", "spline"})))
  70     error("imremap: unsupported interpolation method");
  71   endif
  72   if (any(strcmpi(interp, {"bilinear", "bicubic"})))
  73     interp = interp(3:end); # Remove "bi"
  74   endif
  75   interp = lower(interp);
  76
  77   if (!isscalar(extrapval))
  78     error("imremap: extrapolation value must be a scalar");
  79   endif
  80
  81   ## Interpolate
  82   if (imchannels == 1) # Gray
  83     warped = grayinterp(im, XI, YI, interp, NA);
  84   else # rgb image
  85     for i = 3:-1:1
  86       warped(:,:,i) = grayinterp(im(:,:,i), XI, YI, interp, NA);
  87     endfor
  88   endif
  89   valid = !isna(warped);
  90   warped(!valid) = extrapval;
  91
  92   ## Change the class of the results according to the class of the image
  93   c = class(im);
  94   if (strcmpi(c, "uint8"))
  95     warped = uint8(warped);
  96   elseif (strcmpi(c, "uint16"))
  97     warped = uint16(warped);
  98   endif
  99
 100 endfunction
 101
 102 function [warped, valid] = grayinterp(im, XI, YI, interp, extrapval)
 103   if (strcmp(interp, "cubic"))
 104     warped = graybicubic(double(im), XI, YI, NA);
 105   else
 106     warped = interp2(double(im), XI, YI, interp, NA);
 107   endif
 108   valid = !isna(warped);
 109   warped(!valid) = extrapval;
 110 endfunction
 111
 112 ## -*- texinfo -*-
 113 ## @deftypefn {Function File} {@var{zi}=} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
 114 ## Reference:
 115 ## Image Processing, Analysis, and Machine Vision, 2nd Ed.
 116 ## Sonka et.al.
 117 ## Brooks/Cole Publishing Company
 118 ## ISBN: 0-534-95393-X
 119 ## @seealso{interp2}
 120 ## @end deftypefn
 121
 122 function ZI = graybicubic (Z, XI, YI, extrapval = NA)
 123
 124   ## Allocate output
 125   [X, Y] = meshgrid(1:columns(Z), 1:rows(Z));
 126   [Zr, Zc] = size(XI);
 127   ZI = zeros(Zr, Zc);
 128
 129   ## Find inliers
 130   inside = !( XI < X(1) | XI > X(end) | YI < Y(1) | YI > Y(end) );
 131
 132   ## Scale XI and YI to match indices of Z (not needed when interpolating images)
 133   #XI = (columns(Z)-1) * ( XI - X(1) ) / (X(end)-X(1)) + 1;
 134   #YI = (rows(Z)-1)    * ( YI - Y(1) ) / (Y(end)-Y(1)) + 1;
 135
 136   ## Start the real work
 137   K = floor(XI);
 138   L = floor(YI);
 139
 140   ## Coefficients
 141   AY1  = bc((YI-L+1)); AX1  = bc((XI-K+1));
 142   AY0  = bc((YI-L+0)); AX0  = bc((XI-K+0));
 143   AY_1 = bc((YI-L-1)); AX_1 = bc((XI-K-1));
 144   AY_2 = bc((YI-L-2)); AX_2 = bc((XI-K-2));
 145
 146   ## Perform interpolation
 147   sz = size(Z);
 148   %ZI(inside) = AY_2 .* AX_2 .* Z(sym_sub2ind(sz, L+2, K+2)) ...
 149   ZI = AY_2 .* AX_2 .* Z(sym_sub2ind(sz, L+2, K+2)) ...
 150      + AY_2 .* AX_1 .* Z(sym_sub2ind(sz, L+2, K+1))    ...
 151      + AY_2 .* AX0  .* Z(sym_sub2ind(sz, L+2, K)) ...
 152      + AY_2 .* AX1  .* Z(sym_sub2ind(sz, L+2, K-1)) ...
 153      + AY_1 .* AX_2 .* Z(sym_sub2ind(sz, L+1, K+2)) ...
 154      + AY_1 .* AX_1 .* Z(sym_sub2ind(sz, L+1, K+1))    ...
 155      + AY_1 .* AX0  .* Z(sym_sub2ind(sz, L+1, K)) ...
 156      + AY_1 .* AX1  .* Z(sym_sub2ind(sz, L+1, K-1)) ...
 157      + AY0  .* AX_2 .* Z(sym_sub2ind(sz, L,   K+2)) ...
 158      + AY0  .* AX_1 .* Z(sym_sub2ind(sz, L,   K+1))    ...
 159      + AY0  .* AX0  .* Z(sym_sub2ind(sz, L,   K)) ...
 160      + AY0  .* AX1  .* Z(sym_sub2ind(sz, L,   K-1)) ...
 161      + AY1  .* AX_2 .* Z(sym_sub2ind(sz, L-1, K+2)) ...
 162      + AY1  .* AX_1 .* Z(sym_sub2ind(sz, L-1, K+1))    ...
 163      + AY1  .* AX0  .* Z(sym_sub2ind(sz, L-1, K)) ...
 164      + AY1  .* AX1  .* Z(sym_sub2ind(sz, L-1, K-1));
 165   ZI(!inside) = extrapval;
 166
 167 endfunction
 168
 169 ## Checks if data is meshgrided
 170 function b = isgriddata(X)
 171   D = X - repmat(X(1,:), rows(X), 1);
 172   b = all(D(:) == 0);
 173 endfunction
 174
 175 ## Checks if data is equally spaced (assumes data is meshgrided)
 176 function b = isequallyspaced(X)
 177   Dx = gradient(X(1,:));
 178   b = all(Dx == Dx(1));
 179 endfunction
 180
 181 ## Computes the interpolation coefficients
 182 function o = bc(x)
 183   x = abs(x);
 184   o = zeros(size(x));
 185   idx1 = (x < 1);
 186   idx2 = !idx1 & (x < 2);
 187   o(idx1) = 1 - 2.*x(idx1).^2 + x(idx1).^3;
 188   o(idx2) = 4 - 8.*x(idx2) + 5.*x(idx2).^2 - x(idx2).^3;
 189 endfunction
 190
 191 ## This version of sub2ind behaves as if the data was symmetrically padded
 192 function ind = sym_sub2ind(sz, Y, X)
 193   Y(Y<1) = 1 - Y(Y<1);
 194   while (any(Y(:)>2*sz(1)))
 195     Y(Y>2*sz(1)) = round( Y(Y>2*sz(1))/2 );
 196   endwhile
 197   Y(Y>sz(1)) = 1 + 2*sz(1) - Y(Y>sz(1));
 198   X(X<1) = 1 - X(X<1);
 199   while (any(X(:)>2*sz(2)))
 200     X(X>2*sz(2)) = round( X(X>2*sz(2))/2 );
 201   endwhile
 202   X(X>sz(2)) = 1 + 2*sz(2) - X(X>sz(2));
 203   ind = sub2ind(sz, Y, X);
 204 endfunction
 205
 206 %!demo
 207 %! ## Generate a synthetic image and show it
 208 %! I = tril(ones(100)) + abs(rand(100)); I(I>1) = 1;
 209 %! I(20:30, 20:30) = !I(20:30, 20:30);
 210 %! I(70:80, 70:80) = !I(70:80, 70:80);
 211 %! figure, imshow(I);
 212 %! ## Resize the image to the double size and show it
 213 %! [XI, YI] = meshgrid(linspace(1, 100, 200));
 214 %! warped = imremap(I, XI, YI);
 215 %! figure, imshow(warped);
 216
 217 %!demo
 218 %! ## Generate a synthetic image and show it
 219 %! I = tril(ones(100)) + abs(rand(100)); I(I>1) = 1;
 220 %! I(20:30, 20:30) = !I(20:30, 20:30);
 221 %! I(70:80, 70:80) = !I(70:80, 70:80);
 222 %! figure, imshow(I);
 223 %! ## Rotate the image around (0, 0) by -0.4 radians and show it
 224 %! [XI, YI] = meshgrid(1:100);
 225 %! R = [cos(-0.4) sin(-0.4); -sin(-0.4) cos(-0.4)];
 226 %! RXY = [XI(:), YI(:)] * R;
 227 %! XI = reshape(RXY(:,1), [100, 100]); YI = reshape(RXY(:,2), [100, 100]);
 228 %! warped = imremap(I, XI, YI);
 229 %! figure, imshow(warped);