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