function options = getWarpingFunction(options)
% this function contains all of the known warping functions to be used. if
% a new type is needed please add the warping function here, the forward
% warping function, and the compositional function 

switch options.WarpType
    case 'affine2d'
        options.warpingFcn = @(P) affine2d(P);
        options.fwdWarpingFcn = @(P,x) affine2dFwd(P,x);
        options.composition = @(P,dP) affine2dComposition(P,dP);
end
end

function W = affine2d(P)
W = [1+P(1)     P(3)    P(5)
       P(2)   1+P(4)    P(6)
       0        0       1];
end


function x = affine2dFwd(X,inputs)
narginchk(2,2) % check to make sure 2 inputs are passed
P = inputs{1}; % set P = to the first argin
if numel(inputs)==1
    K = eye(3);
else
    K = inputs{2};
end
W = affine2d(P);
x = K.inv*W*K.fwd*X;
end

function P = affine2dComposition(P,dP)
% This composition is the same as building a warp for P and dP:
% W(x,P):                           W(x,dP):
% [1+P(1)  P(3)    P(5) ;           [1+dP(1)  dP(3)    dP(5) ;
%  P(2)    1+P(4)  P(6) ;            dP(2)    1+dP(4)  dP(6) ;
%  0       0       1   ];            0       0       1      ];
% from the parameters P and dP, then muliplying the resulting matrices:
% W(x,P)*W(x,dP)^-1
%
% derived from equation 16
%
%Here it is calculated explicity for speed:

% den = ( (1+dP(1))*(1+dP(4))- dP(2)*dP(3) ) ^-1;
% 
% dP = [-dP(1) - dP(1)*dP(4) + dP(2)*dP(3)
%                     -dP(2)
%                     -dP(3)
%       -dP(4) - dP(1)*dP(4) + dP(2)*dP(3)
%       -dP(5) - dP(4)*dP(5) + dP(3)*dP(6)
%       -dP(6) - dP(1)*dP(6) + dP(2)*dP(5)];
%   
% dP = den*dP;
%   
% P =[P(1) + dP(1) + P(1)*dP(1) + P(3)*dP(2) ...
%     P(2) + dP(2) + P(2)*dP(1) + P(4)*dP(2) ...
%     P(3) + dP(3) + P(1)*dP(3) + P(3)*dP(4) ...
%     P(4) + dP(4) + P(2)*dP(3) + P(4)*dP(4) ...
%     P(5) + dP(5) + P(1)*dP(5) + P(3)*dP(6) ...
%     P(6) + dP(6) + P(2)*dP(5) + P(4)*dP(6)  ];

PW =  affine2d(P);
dPW = affine2d(dP);
W = PW/dPW;
W = W-eye(3);
W(end,:) = [];
P = W(:)';

end


   