1## Copyright (C) 2012 Pantxo Diribarne 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 3 of the License, or 6## (at your option) any later version. 7## 8## This program is distributed in the hope that it will be useful, 9## but WITHOUT ANY WARRANTY; without even the implied warranty of 10## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11## GNU General Public License for more details. 12## 13## You should have received a copy of the GNU General Public License 14## along with Octave; see the file COPYING. If not, see 15## <http://www.gnu.org/licenses/>. 16 17## -*- texinfo -*- 18## @deftypefn {Function File} {@var{IM_OUT} =} imtransform (@var{IM_IN}, @var{T}) 19## @deftypefnx {Function File} {@var{IM_OUT} =} imtransform (@var{IM_IN}, @var{T}, @var{interp}) 20## @deftypefnx {Function File} {@var{IM_OUT} =} imtransform (@dots{}, @var{prop}, @var{val}) 21## @deftypefnx {Function File} {[@var{IM_OUT}, @var{xdata}, @var{ydata}] =} imtransform (@dots{}) 22## Transform image. 23## 24## Given an image @var{IM_IN}, return an image @var{IM_OUT} 25## resulting from the forward transform defined in the transformation 26## structure @var{T}. An additional input argument @var{interp}, 27## 'bicubic', 'bilinear' (default) or 'nearest', 28## specifies the interpolation method to be used. Finally, the 29## transform can be tuned using @var{prop}/@var{val} pairs. The 30## following properties @var{prop} are supported: 31## 32## @table @asis 33## @item "udata" 34## Specifies the input space horizontal limits. @var{val} must be a two 35## elements vector [minval maxval]. Default: [1 columns(@var{IM_IN})] 36## 37## @item "vdata" 38## Specifies the input space vertical limits. @var{val} must be a two 39## elements vector [minval maxval]. Default: [1 rows(@var{IM_IN})] 40## 41## @item "xdata" 42## Specifies the required output space horizontal limits. @var{val} must 43## be a two elements vector [minval maxval]. Default: estimated using 44## udata, vdata and findbounds function. 45## 46## @item "ydata" 47## Specifies the required output space vertical limits. @var{val} must 48## be a two elements vector [minval maxval]. Default: estimated using 49## udata, vdata and findbounds function. 50## 51## @item "xyscale" 52## Specifies the size of pixels in the output space. If a scalar 53## is provided, both vertical and horizontal dimensions are scaled the 54## same way. If @var{val} is a two element vector, it must indicate 55## consecutively width and height of the output pixels. The default is 56## to use the width and height of input pixels provided 57## that it does not lead to a too large output image. 58## 59## @item "size" 60## Size of the output image (1-by-2 vector). Overrides the effect of 61## "xyscale" property. 62## 63## @item "fillvalues" 64## Color of the areas where no interpolation is possible, e.g. when 65## coordinates of points in the output space are out of the limits of 66## the input space. @var{val} must be coherent with the input image 67## format: for grayscale and indexed images (2D) @var{val} must be 68## scalar, for RGB (n-by-m-by-3) @var{val} must be a 3 element vector. 69## 70## @end table 71## 72## The actual output limits, @var{xdata} and @var{ydata} vectors, 73## are returned respectively as second and third output variables. 74## @seealso{maketform, cp2tform, tforminv, tformfwd, findbounds} 75## @end deftypefn 76 77## Author: Pantxo Diribarne <pantxo@dibona> 78 79function [varargout] = imtransform (im, T, varargin) 80 81 if (nargin < 2) 82 print_usage (); 83 elseif (! istform (T)) 84 error ("imtransform: T must be a transformation structure (see `maketform')"); 85 endif 86 87 ## Parameters 88 interp = "linear"; 89 imdepth = size (im, 3); 90 maximsize = [20000 20000]; 91 92 udata = [1; columns(im)]; 93 vdata = [1; rows(im)]; 94 xdata = ydata = []; 95 xyscale = []; 96 imsize = []; 97 fillvalues = zeros (1, imdepth); 98 99 if (isempty (varargin)) 100 xydata = findbounds (T, [udata vdata]); 101 xdata = xydata(:,1); 102 ydata = xydata(:,2); 103 else 104 ## interp 105 if (floor (numel (varargin)/2) != numel (varargin)/2) 106 allowed = {"bicubic", "bilinear", "nearest"}; 107 tst = strcmp (varargin{1}, allowed); 108 if (!any (tst)) 109 error ("imtransform: expect one of %s as interp method", disp (allowed)); 110 else 111 interp = {"pchip", "linear", "nearest"}{find (tst)}; 112 endif 113 varargin = varargin(2:end); 114 endif 115 116 ## options 117 allowed = {"udata", "vdata", "xdata", "ydata", ... 118 "xyscale", "size", "fillvalues"}; 119 props = varargin(1:2:end); 120 vals = varargin(2:2:end); 121 np = numel (props); 122 if (!all (cellfun (@ischar, props))) 123 error ("imtransform: expect property/value pairs."); 124 endif 125 126 props = tolower (props); 127 tst = cellfun (@(x) any (strcmp (x, allowed)), props); 128 if (!all (tst)) 129 error ("imtransform: unknown property %s", disp (props{!tst})); 130 endif 131 132 ## u(vxy)data 133 iolims = allowed(1:4); 134 for ii = 1:numel (iolims) 135 tst = cellfun (@(x) any (strcmp (x, iolims{ii})), props); 136 if (any (tst)) 137 prop = props{find (tst)(1)}; 138 val = vals{find (tst)(1)}; 139 if (isnumeric (val) && numel (val) == 2) 140 if (isrow (val)) 141 val = val'; 142 endif 143 eval (sprintf ("%s = val;", prop), 144 "error (\"imtransform: %s\n\", lasterr ());"); 145 else 146 error ("imtransform: expect 2 elements real vector for %s", prop) 147 endif 148 endif 149 endfor 150 if (isempty (xdata) && isempty (ydata)) 151 xydata = findbounds (T, [udata vdata]); 152 xdata = xydata(:,1); 153 ydata = xydata(:,2); 154 elseif (isempty (xdata)) 155 xydata = findbounds (T, [udata vdata]); 156 xdata = xydata(:,1); 157 elseif (isempty (ydata)) 158 xydata = findbounds (T, [udata vdata]); 159 ydata = xydata(:,2); 160 endif 161 162 ## size and xyscale 163 tst = strcmp ("size", props); 164 if (any (tst)) 165 val = vals{find (tst)(1)}; 166 if (isnumeric (val) && numel (val) == 2 && 167 all (val > 0)) 168 imsize = val; 169 else 170 error ("imtransform: expect 2 elements real vector for size"); 171 endif 172 elseif (any (tst = strcmp ("xyscale", props))) 173 val = vals{find (tst)(1)}; 174 if (isnumeric (val) && all (val > 0)) 175 if (numel (val) == 1) 176 xyscale(1:2) = val; 177 elseif (numel (val) == 2) 178 xyscale = val; 179 else 180 error ("imtransform: expect 1 or 2 element(s) real vector for xyscale"); 181 endif 182 else 183 error ("imtransform: expect 1 or 2 elements real vector for xyscale"); 184 endif 185 endif 186 ## Fillvalues 187 tst = strcmp ("fillvalues", props); 188 if (any (tst)) 189 val = vals{find (tst)(1)}; 190 if (isnumeric (val) && numel (val) == 1) 191 fillvalues(1:end) = val; 192 elseif (isnumeric (val) && numel (val) == 3) 193 fillvalues = val; 194 else 195 error ("imtransform: expect 1 or 3 elements real vector for `fillvalues'"); 196 endif 197 endif 198 endif 199 200 ## Output/Input pixels 201 if (isempty (imsize)) 202 if (isempty (xyscale)) 203 nx = columns (im); 204 ny = rows (im); 205 xyscale = [(diff (udata) / (nx - 1)), ... 206 (diff (vdata) / (ny - 1))]; 207 endif 208 xscale = xyscale(1); 209 yscale = xyscale(2); 210 xsize = ceil (diff (xdata) / xscale) + 1; 211 ysize = ceil (diff (ydata) / yscale) + 1; 212 ## xyscale takes precedence: recompute x/ydata considering roundoff errors 213 xdata(2) = xdata(1) + (xsize - 1) * xscale; 214 ydata(2) = ydata(1) + (ysize - 1) * yscale; 215 if (xsize > maximsize(2) || ysize > maximsize(1)) 216 if (xsize >= ysize) 217 scalefactor = (diff (xdata) / maximsize(2)) / xscale; 218 else 219 scalefactor = (diff (ydata) / maximsize(1)) / yscale; 220 endif 221 xscale *= scalefactor; 222 yscale *= scalefactor; 223 xsize = floor (diff (xdata) / xscale); 224 ysize = floor (diff (ydata) / yscale); 225 warning ("imtransform: output image two large, adjusting the largest dimension to %d", maximsize); 226 endif 227 imsize = [ysize xsize]; 228 endif 229 [xx yy] = meshgrid (linspace (xdata(1), xdata(2), imsize(2)), 230 linspace (ydata(1), ydata(2), imsize(1))); 231 232 [uu vv] = meshgrid (linspace (udata(1), udata(2), size(im)(2)), 233 linspace (vdata(1), vdata(2), size(im)(1))); 234 235 ## Input coordinates 236 [uui, vvi] = tforminv (T, reshape (xx, numel (xx), 1), 237 reshape (yy, numel (yy), 1)); 238 uui = reshape (uui, size (xx)); 239 vvi = reshape (vvi, size (yy)); 240 ## Interpolation 241 for layer = 1:imdepth 242 imout(:,:,layer) = interp2 (uu, vv, im(:,:,layer), ... 243 uui, vvi, interp, fillvalues(layer)); 244 endfor 245 if (nargout != 0) 246 varargout = {imout, xdata(:).', ydata(:).'}; 247 endif 248endfunction 249 250%!demo 251%! ## Various linear transforms 252%! figure (); 253%! im = [checkerboard(20, 2, 4); checkerboard(40, 1, 2)]; 254%! %input space corners 255%! incp = [1 1; 160 1; 160 160; 1 160]; 256%! udata = [min(incp(:,1)) max(incp(:,1))]; 257%! vdata = [min(incp(:,2)) max(incp(:,2))]; 258%! subplot (2,3,1); 259%! imshow (im) 260%! hold on 261%! plot (incp(:,1), incp(:,2), 'ob') 262%! axis on 263%! xlabel ('Original') 264%! 265%! % Translation and scaling 266%! outcp = incp * 2; 267%! outcp(:,1) += 200; 268%! outcp(:,2) += 500; 269%! T = maketform ('affine', incp(1:3,:), outcp(1:3,:)); 270%! subplot (2,3,2); 271%! [im2 xdata ydata] = imtransform (im, T, 'udata', udata, 272%! 'vdata', vdata, 'fillvalues', 1); 273%! imh = imshow (im2); set (imh, 'xdata', xdata, 'ydata', ydata) 274%! set (gca, 'xlim', xdata, 'ylim', ydata) 275%! axis on, hold on, xlabel ('Translation / Scaling'); 276%! plot (outcp(:,1), outcp(:,2), 'or') 277%! 278%! % Shear 279%! outcp = [1 1; 160 1; 140 160; -19 160]; % affine only needs 3 control points 280%! T = maketform ('affine', incp(1:3,:), outcp(1:3,:)); 281%! subplot (2,3,3); 282%! [im2 xdata ydata] = imtransform (im, T, 'udata', udata, 283%! 'vdata', vdata, 'fillvalues', 1); 284%! imh = imshow (im2); set (imh, 'xdata', xdata, 'ydata', ydata) 285%! set (gca, 'xlim', xdata, 'ylim', ydata) 286%! axis on, hold on, xlabel ('Shear'); 287%! plot (outcp(:,1), outcp(:,2), 'or') 288%! 289%! % Rotation 290%! theta = pi/4; 291%! T = maketform ('affine', [cos(theta) -sin(theta); ... 292%! sin(theta) cos(theta); 0 0]); 293%! outcp = tformfwd (T, incp); 294%! subplot (2,3,4); 295%! [im2 xdata ydata] = imtransform (im, T, 'udata', udata, 296%! 'vdata', vdata, 'fillvalues', 1 ); 297%! imh = imshow (im2); set (imh, 'xdata', xdata, 'ydata', ydata) 298%! set (gca, 'xlim', xdata, 'ylim', ydata) 299%! axis on, hold on, xlabel ('Rotation'); 300%! plot (outcp(:,1), outcp(:,2), 'or') 301%! 302%! % Reflection around x axis 303%! outcp = incp; 304%! outcp(:,2) *= -1; 305%! T = cp2tform (incp, outcp, 'similarity'); 306%! subplot (2,3,5); 307%! [im2 xdata ydata] = imtransform (im, T, 'udata', udata, 308%! 'vdata', vdata, 'fillvalues', 1 ); 309%! imh = imshow (im2); set (imh, 'xdata', xdata, 'ydata', ydata) 310%! set (gca, 'xlim', xdata, 'ylim', ydata) 311%! axis on, hold on, xlabel ('Reflection'); 312%! plot (outcp(:,1), outcp(:,2), 'or') 313%! 314%! % Projection 315%! outcp = [1 1; 160 -40; 220 220; 12 140]; 316%! T = maketform ('projective', incp, outcp); 317%! subplot (2,3,6); 318%! [im2 xdata ydata] = imtransform (im, T, 'udata', udata, 319%! 'vdata', vdata, 'fillvalues', 1 ); 320%! imh = imshow (im2); set (imh, 'xdata', xdata, 'ydata', ydata) 321%! set (gca, 'xlim', xdata, 'ylim', ydata) 322%! axis on, hold on, xlabel ('Projection'); 323%! plot (outcp(:,1), outcp(:,2), 'or') 324 325%!demo 326%! ## Streched image 327%! rad = 2; % minimum value: 4/pi 328%! [uu vv] = meshgrid ((-2:2)/rad, (-2:2)/rad); 329%! rescfactor = sin ((uu.^2 + vv.^2).^.5); 330%! inpts = [(reshape (uu, numel (uu), 1)), (reshape (vv, numel (uu), 1))]; 331%! xx = rescfactor .* sign(uu); 332%! yy = rescfactor .* sign(vv); 333%! outpts = [reshape(xx, numel (xx), 1) reshape(yy, numel (yy), 1)]; 334%! 335%! T = cp2tform (inpts, outpts, "polynomial", 4); 336%! figure; 337%! subplot (1,2,1) 338%! im = zeros (800, 800, 3); 339%! im(:,:,1) = checkerboard (100) > 0.2; 340%! im(:,:,3) = checkerboard (100) < 0.2; 341%! [im2 xdata ydata] = imtransform (im, T, 'udata', [-2 2], 342%! 'vdata', [-2 2], 'fillvalues', 343%! [0 1 0]); 344%! imh = imshow (im2); 345%! set (imh, 'xdata', xdata, 'ydata', ydata) 346%! set (gca, 'xlim', xdata, 'ylim', ydata) 347%! [im cmap] = imread ('default.img'); 348%! subplot (1,2,2) 349%! [im2 xdata ydata] = imtransform (im, T, 'udata', [-1 1], 350%! 'vdata', [-1 1], 'fillvalues', 351%! round (length (cmap) / 2)); 352%! imh = imshow (im2, cmap); 353 354%!test 355%! im = checkerboard (); 356%! incp = [0 0; 0 1; 1 1]; 357%! scl = 10; 358%! outcp = scl * incp; 359%! T = maketform ('affine', incp, outcp); 360%! [im2 xdata ydata] = imtransform (im, T, 'udata', [0 1], 361%! 'vdata', [0 1], 'size', [500 500]); 362%! assert (xdata, scl * ([0 1])) 363%! assert (ydata, scl * ([0 1])) 364%! assert (size (im2), [500 500]) 365 366%!test 367%! im = checkerboard (); 368%! incp = [0 0; 0 1; 1 1]; 369%! scl = 10; 370%! outcp = scl * incp; 371%! xyscale = scl; 372%! T = maketform ('affine', incp, outcp); 373%! [im2 xdata ydata] = imtransform (im, T, 'xyscale', xyscale); 374%! assert (size (im2), size (im), 1) 375 376%!test 377%! im = checkerboard (100, 10, 4); 378%! theta = 2 * pi; 379%! T = maketform ("affine", [cos(theta) -sin(theta); ... 380%! sin(theta) cos(theta); 0 0]); 381%! im2 = imtransform (im, T, "nearest", "xdata", [1 800], "ydata", [1 2000]); 382%! im = im(2:end-1, 2:end-1); %avoid boundaries 383%! im2 = im2(2:end-1, 2:end-1); 384%! assert (im, im2) 385 386%!test 387%! im = checkerboard (20, 10, 4); 388%! theta = pi/6; 389%! T = maketform ('affine', [cos(theta) -sin(theta); ... 390%! sin(theta) cos(theta); 0 0]); 391%! [im2, xdata] = imtransform (im, T); 392%! nu = columns(im); 393%! nv = rows(im); 394%! nx = xdata(2); 395%! diag = sqrt (nu^2 + nv^2); 396%! ang = atan (nv / nu); 397%! assert (nx, diag * abs (cos (theta - ang)), 398%! diag * 1 / size (im2, 2)) 399 400## Test default fillvalues 401%!test 402%! im = rand (2); 403%! tmat = [eye(2); 0 0]; 404%! T = maketform ("affine", tmat); 405%! im2 = imtransform (im, T, "xdata", [1 3]); 406%! assert (im2(:,3), zeros (2,1)) 407 408## Test image size when forcing x/ydata 409%!test 410%! im = rand (2); 411%! tmat = [eye(2); 0 0]; 412%! T = maketform ('affine', tmat); 413%! im2 = imtransform (im, T, "xdata", [1 3]); 414%! assert (size (im2), [2 3]) 415 416## Test image size when forcing xyscale 417%!test 418%! im = rand (2); 419%! tmat = [eye(2); 0 0]; 420%! T = maketform ('affine', tmat); 421%! im2 = imtransform (im, T, "xyscale", [0.5 0.5]); 422%! assert (size (im2), [3 3]) 423