function [dd, f1,f12,f41,f1324, g1,g12,g41,g1324, c1,c2,c3] = MLS_applyDeform2(x,y,X,Y,Fx,Fy, f1,f12,f41,f1324, g1,g12,g41,g1324, c1,c2,c3)
% f = deformation in x
% g = deformation in y
% X = matrix containing x position of vertex
% Y = matrix containing y position of vertex
% f00, dfx0, df0x, ddf; g00, dgx0, dg0x, ddg; c1, c2, c3  
% - constans for computation, 10x faster when prealocated
% - for first run zeros(size(X)-1) - done in MLS_getPoints, kept saved in MLS structure 

% Coordinate system f(j,i)
%                                       y,j
%                                          A
%  f1              f2                   f4 |            f3
%     +-----------+----> x,i               +-----------+
%     |           |                        |           |
%     |           |             ==         |           |
%     |           |                        |           |
%     |           |                        |           |
%     +-----------+                        +-----------+----> x,i
%  f4 |            f3                    f1             f2
%     V
%     y,j


%% Preparation

[f1,f12,f41,f1324] = constants(Fx);
[g1,g12,g41,g1324] = constants(Fy);
a = g1324.*f12 - f1324.*g12;
b = g41.*f12 + f1324.*g1 - f41.*g12 - g1324.*f1;
c = f41.*g1 - g41.*f1;


minf = floor(bsxfun(@min,Fx(1:end-1,1:end-1),Fx(2:end,1:end-1)));minf(minf<1)=1;
maxf = ceil(bsxfun(@max,Fx(1:end-1,2:end),Fx(2:end,2:end)));maxf(maxf>x(end))=x(end);

ming = floor(bsxfun(@min,Fy(1:end-1,1:end-1),Fy(1:end-1,2:end)));ming(ming<1)=1;
maxg = ceil(bsxfun(@max,Fy(2:end,1:end-1),Fy(2:end,2:end)));maxg(maxg>y(end))=y(end);

FF = [reshape(f1,[],1),reshape(f12,[],1),reshape(f41,[],1),reshape(f1324,[],1)];
% %%
% 
% % for all squares formed by discrete verteces 
% % redo: for provided vector (pairs of i and j)
% for i=1:length(x)-1
%     for j=1:length(y)-1
% 
%         for k = minf(j,i):maxf(j,i)
%             for m = ming(j,i):maxg(j,i)
%                 q = i+j-1;
%                 xx = sqrt(FF(q,1).*FF(q,2).*FF(q,3).*FF(q,4))*k;
%             end
%         end
%         % if max(abs(F))&max(abs(G))<1 copy to new image directly
%         
%         % compute substitution constants used by all points
%         %fprintf('[%d,%d]-[%d,%d] ',x(i),y(j),x(i+1),y(j+1)) 
%         % for all whole number points in range min(F):max(F), min(G):max(G)
% 
% 
%         % get x,y in original image
% 
%         % get interpolated pixel and write it to newImg
%     end
%     %fprintf('\n')
% end

% for all squares formed by discrete verteces 
% redo: for provided vector (pairs of i and j)
dd = 0;
for i=1:length(x)-1
    for j=1:length(y)-1
        for k = minf(j,i):maxf(j,i)
            for m = ming(j,i):maxg(j,i)
                dd = sqrt(f1(j,i).*f12(j,i).*f41(j,i).*f1324(j,i))*k;
            end
        end
        % if max(abs(F))&max(abs(G))<1 copy to new image directly
        
        % compute substitution constants used by all points
        %fprintf('[%d,%d]-[%d,%d] ',x(i),y(j),x(i+1),y(j+1)) 
        % for all whole number points in range min(F):max(F), min(G):max(G)


        % get x,y in original image

        % get interpolated pixel and write it to newImg
    end
    %fprintf('\n')
end

end
%% helper functions
function [f1,f12,f41,f1324] = constants(Fx)
f1 = Fx(1:end-1,1:end-1);
f12 = conv2(Fx,flip(fliplr([1 -1;0 0])),'valid');
f41 = conv2(Fx,flip(fliplr([-1 0;1 0])),'valid');
f1324 = conv2(Fx,flip(fliplr([1 -1;-1 1])),'valid');

end