      1function [newImg,finterp] = MLS_applyDeform(oldImg,f,v,grid_size)
      2%UNTITLED2 Summary of this function goes here
      3%   f(v) - deformation vector
      4%   v - point of evaluation
      5%   grid_size - loop info for v points
      6[imh,imw,imd] = size(oldImg);
      7newImg = oldImg;
      8% newImg = ones(imh,imw,imd);
      9newImg(:,:,1) = 50;
     10newImg(:,:,2) = 205;
     11newImg(:,:,3) = 50;
     12% finterp.x = [];
     13% finterp.y = [];
     14% finterp.zX = [];
     15% finterp.zY = [];
     16finterp = 0;
     17% iterate over evaluated points
     18% +--------> x,i=width
     19% |
     20% |
     21% |
     22% V y,j=height
     23stepX = grid_size(2);
     24stepY = grid_size(1);
     25
     26% go over grid
     27for i = 1:stepX:imw
     28    for j = 1:stepY:imh
     29        % next i,j to create square for bilinear interp
     30        ni = i+stepX;
     31        nj = j+stepY;
     32
     33        % length of square sides
     34        x2 = stepX;
     35        y2 = stepY;
     36
     37        % check image borders
     38        if (ni > imw), ni=imw; x2=imw-i;end
     39        if (nj > imh), nj=imh; y2=imh-j;end
     40
     41        %   y A            ni,nj
     42        %     |            x2,y2
     43        % f01 |-----------+  f11
     44        %     |           |
     45        %     |           |
     46        %     |           |
     47        %     |           |
     48        % f00 +----------------> x
     49        %   x1,y1 = 0      f10
     50        %   i,j            ni,j
     51        % call interpolation
     52        [newImg,finterp] = Interpolation(i,j,ni,nj,x2,y2,f,finterp,newImg,imh,imw,imd,oldImg);
     53
     54    end
     55end
     56
     57end
     58%% bilinear interpolation for 1*1 square,
     59function ret = bilinear_interpol(x,y,f00,f10,f01,f11)
     60% bilinear interpolation for 1*1 square, if you pass x=x/x2 and y=y/y2 you
     61% get bilinear interpolation for rectangle in begining of cordinate system
     62ret = (f00*(1-x) + f10*x)*(1-y) + (f01*(1-x) + f11*x)*y;
     63end
     64
     65function ret = first_half_bilinear_interpolation(x,f00,f10,f01,f11)
     66ret = [(f00*(1-x) + f10*x),(f01*(1-x) + f11*x)];
     67end
     68function ret = second_half_bilinear_interpolation(y,p)
     69ret = p(1)*(1-y)+p(2)*y;
     70end
     71function slope = biggest_slope(A,B,C,D)
     72slope = abs([A-B,A-C,A-D,B-C,B-D,C-D]);
     73end
     74
     75%% Interpolation
     76function [newImg,finterp] = Interpolation(i,j,ni,nj,x2,y2,f,finterp,newImg,imh,imw,imd,oldImg)
     77% calculate step form slopes
     78
     79% try
     80step = min([x2,y2])/max([biggest_slope(f(j,i,1), f(j,ni,1), f(nj,i,1), f(nj,ni,1)),biggest_slope(f(j,i,2), f(j,ni,2), f(nj,i,2), f(nj,ni,2))]);
     81% catch e
     82%     f
     83%     fprintf('Size f: %f×%f×%f, i:%f, ni:%f, j:%f, nj:%f',size(f),i,ni,j,nj)
     84%     rethrow e
     85% end
     86
     87for x = 0:step:x2
     88    % fill deformation function for area between control points (1/2)
     89    py = first_half_bilinear_interpolation(x/x2,f(j,i,1), f(j,ni,1), f(nj,i,1), f(nj,ni,1));
     90    px = first_half_bilinear_interpolation(x/x2,f(j,i,2), f(j,ni,2), f(nj,i,2), f(nj,ni,2));
     91
     92    for y = 0:step:y2
     93        fy = second_half_bilinear_interpolation(y/y2,py);
     94        fx = second_half_bilinear_interpolation(y/y2,px);
     95
     96
     97%          finterp.x = [finterp.x;x+i];
     98%          finterp.y = [finterp.y;y+j];
     99%          finterp.zX = [finterp.zX;fx];
    100%          finterp.zY = [finterp.zY;fy];
    101
    102
    103        % check borders
    104        nx = fx;
    105        ny = fy;
    106        if nx > imw, nx = imw; end
    107        if nx < 1, nx = 1; end
    108
    109        if ny > imh, ny = imh; end
    110        if ny < 1, ny = 1; end
    111
    112        % round dest. pixels
    113        f_nx = round(nx);
    114
    115        f_ny = round(ny);
    116
    117        % float from rounding
    118        fpx =  round(nx-(i+x))-(fx-(i+x));
    119        fpy =  round(ny-(j+y))-(fy-(j+y));
    120
    121        % find closest source pixels 1*1 square
    122        x11 = floor(i+x+fpx);
    123        if x11 < 1
    124            x11 = 1; x22 = 1; fpx = 0;
    125        elseif x11 > imw
    126            x11 = imw; x22 = imw; fpx = 0;
    127        else
    128            x22 = ceil(i+x+fpx);
    129            fpx = i+x+fpx-x11; % distance on unit square for interpolation
    130        end
    131
    132
    133        y11 = floor(j+y+fpy);
    134        if y11 < 1
    135            y11 = 1;y22=1;fpy = 0;
    136        elseif y11 > imh
    137            y11 = imh;y22=imh;fpy = 0;
    138        else
    139            y22 = ceil(j+y+fpy)  ;
    140            fpy = j+y+fpy-y11;
    141        end
    142
    143
    144        newImg = interpolate_pixels(x11,x22,y11,y22,f_nx,f_ny,newImg, fpx, fpy, oldImg,imw,imh,imd);
    145    end
    146end
    147end

    148
    149function newImg = interpolate_pixels(x11,x22,y11,y22,f_nx,f_ny,newImg, fpx, fpy, oldImg,imw,imh,imd)
    150if x22 > imw, x22 = imw;end
    151if x11 > imw, x11 = imw;end
    152
    153if y11 > imh, y11 = imh;end
    154if y22 > imh, y22 = imh;end
    155
    156for layer = 1:imd
    157
    158    % interpolate the pixel at x+x2+float,y+y2+float and
    159    % move it to nx,ny
    160    %try
    161    newImg(f_ny, f_nx, layer) = bilinear_interpol(...
    162        fpx, fpy,...
    163        oldImg(y11, x11,layer),...
    164        oldImg(y11, x22,layer),...
    165        oldImg(y22, x11,layer),...
    166        oldImg(y22, x22,layer));
    167    %catch
    168    %    fprintf('f_ny: %f, f_nx: %f, layer: %f, x11: %f, x22: %f, y11: %f,y22: %f',f_ny,f_nx,x11,x22,y11,y22)
    169    %end
    170    end
    171end