## Copyright (C) 2015 Motherboard
## 
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.

## Author: Motherboard <cantfind@gmail.com>
## Created: 2015-12-1

## -*- texinfo -*-
## affine2d takes a transpose of the affine matrix as described in standard literature
## it performs the transformation as follows:
## v = (u*T)(1:2), where u = [x y 1] and T = [a b 0; c d 0; e f 1] 
## [a b; c d] is a transposed rotation\shear matrix, [e f] is the translation vector, where e = dx, f = dy.

classdef affine2d
  properties
    T                     # Transformation matrix.
    Dimensionality = 2;
  end
  methods
## constructor, if no argment is given, default is the identity matrix.
    function obj = affine2d(tform)
      if (nargin == 0)
        obj.T = eye(3);
      elseif (ndims(tform) == 2 && size(tform) == [3 3] && tform(:,3) == [0; 0; 1]) 
##    check tform is a valid affine transform
        if (det(tform(1:2,1:2)) ~= 0)
          obj.T = tform;
        else
          error("given transform has singular matrix")
        endif	         
      else
        error("given transform is not affine! should be of the form [a b 0; c d 0; e f 1]")	
      endif
    endfunction
##  invert transform
    function obj = invert(tform)
      inv_T = inv(tform.T);
      inv_T(:,3) = [0; 0; 1];
      obj = affine2d (inv_T);
    endfunction
    function ans = isRigid(tform) 
##  check if transform is only rotation or translation
      S = tform.T(1:2,1:2);
      ans = abs(S(1,1)*S(2,2)-S(1,2)*S(2,1)-1) < eps;
    endfunction
    function ans = isSimilarity(tform) 
##  check if transform is only homogeneous scaling, rotation, or translation
      submatrix = tform.T(1:2,1:2);
      s2 = submatrix*submatrix';
      ans = all((abs(s2 - s2(1,1)*eye(2)) < ones(2)*eps)(:));
    endfunction
    function ans = isTranslation(tform) 
##  check if transform is only a translation
      submatrix = tform.T(1:2,1:2);
      ans = all((abs(submatrix - eye(2)) < ones(2)*eps)(:));
    endfunction
    function [limitsX limitsY] = outputLimits(tform, xlims, ylims) 
## given a bounding box corner coordinated in xlims and ylims (top left, right bottom) - return the new bounding box after transformation.
      xlims2 = [xlims(:); xlims(:)];
      ylims2 = [ylims(:); ylims(end:-1:1)(:)];
      temp = [xlims2 ylims2 ones(4,1)];
      ans = temp*tform.T;
      limitsX = [min(ans(1:4,1)), max(ans(1:4,1))];
      limitsY = [min(ans(1:4,2)), max(ans(1:4,2))];
    endfunction
    function [x,y] = transformPointsForward(tform,u,v) 
## apply transformation on the set of u, v points (1xn vectors) or on U (2xn matrix)
      N = size(u,1);
      if exist("v","var")
        u = [u(:) v(:)];
      elseif size(u) ~= [N 2]
        error("provided matrix to transform should be 2 by N in size");
      endif
      U = [u ones(N,1)];
      x = (U*tform.T)(:,1:2);
      if isargout(2)
        y = x(:,2);
        x = x(:,1);
      endif
    endfunction
    function [x,y] = transformPointsInverse(tform,u,v) 
##  apply inverse transformation on the set of u, v points (1xn vectors) or on U (2xn matrix)
      N = size(u,1);
      if exist("v","var")
        u = [u(:) v(:)];
      elseif size(u) ~= [N 2]
        error("provided matrix to transform should be 2 by N in size");
      endif
      U = [u ones(N,1)];
      x = (U*inv(tform.T))(:,1:2);
      if isargout(2)
        y = x(:,2);
        x = x(:,1);
      endif
    endfunction
  end
end

%!test
%! theta = 10;
%! tform = affine2d([cosd(theta) -sind(theta) 0; sind(theta) cosd(theta) 0; 0 0 1]);
%! [X,Y] = transformPointsForward(tform,5,10);
%! assert(abs([X Y] - [6.6605 8.9798]) < e^-4) 
%! [U,V] = transformPointsInverse(tform,X,Y);
%! assert(abs([U V] - [5 10]) < e^-10) 
%! assert(isRigid(tform))
%! assert(~isTranslation(tform))
%! assert(isSimilarity(tform));

%!test
%! tform = affine2d([1 0 0; 0 1 0; 5 10 1]);
%! [X Y] = transformPointsForward(tform,[1 2; 3 4; 5 6; 7 8]);
%! assert(round(X) == [6;8;10;12])
%! assert(round(Y) == [12;14;16;18])
%! [U,V] = transformPointsInverse(tform,X,Y);
%! assert(round(U) == [1;3;5;7])
%! assert(round(V) == [2;4;6;8])
%! assert(isRigid(tform))
%! assert(isTranslation(tform))
%! assert(isSimilarity(tform));

%!test
%! tform = affine2d([1 1e-16 0; 1e-16 1 0; 5 10 1]);
%! assert(isRigid(tform))
%! tform = affine2d([2 1e-16 0; 1e-16 1 0; 5 10 1]);
%! assert(~isRigid(tform))

%!test 
%! theta = 10;
%! tform = affine2d([cosd(theta) -sind(theta) 0; sind(theta) cosd(theta) 0; 0 0 1]);
%! [xlim ylim] = outputLimits(tform,[1 240],[1 291]);
%! assert(abs(xlim - [1.1585  286.8855]) < [1 1]*e^-4)
%! assert(abs(ylim - [-40.6908  286.4054]) < [1 1]*e^-4)

%!test
%!error affine2d ([0 0 0; 0 0 0])
%!error affine2d ([0 0 0  0 0 0  0 0 1])
%!error affine2d ([0 0 0; 0 0 0; 0 0 0])
%!error affine2d ([0 0 0; 0 0 0; 0 0 1])