      1function [newImg,f] = 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% iterate over evaluated points
      9% +--------> x,i=width
     10% |
     11% |
     12% |
     13% V y,j=height
     14stepX = grid_size(2);
     15stepY = grid_size(1);
     16
     17% for i = 1:stepX:imw+stepX
     18%
     19%     if i >= imw; i = imw; end
     20%
     21%     for j = 1:stepY:imh+stepY
     22%         if j == 1; continue;end % skip the [1,1] point (already created)
     23%         if j >= imh; j = imh; end % hit the end row of column
     24%         v = [v;[j,i]] ;           % append the eval point
     25%         if j == imh; break;end    % next colimn
     26%     end
     27%
     28%     if i == imw; break;end        % end after last columm
     29% end
     30
     31
     32%%%
     33for i = 1:stepX:imw
     34    for j = 1:stepY:imh
     35        % next i,j to create square for bilinear interp
     36        ni = i+stepX;
     37        nj = j+stepY;
     38
     39        % length of square sides
     40        x2 = stepX;
     41        y2 = stepY;
     42
     43        % check image borders
     44        if (ni > imw), ni=imw; x2=imw-i;end
     45        if (nj > imh), nj=imh; y2=imh-j;end
     46
     47        %     A
     48        %     |            x2,y2
     49        % f01 |-----------+  f11
     50        %     |           |
     51        %     |           |
     52        %     |           |
     53        %     |           |
     54        % f00 +---------------->
     55        %   x1,y1 = 0      f10
     56
     57        for x = 0:x2-1
     58            for y = 0:y2-1
     59                % fill deformation function for area between control points
     60                % bilinear_interpol(x,y,f00,f10,f01,f11),      f[y,x,[y,x]]
     61                fy = bilinear_interpol(x/x2, y/y2, f(j,i,1), f(j,ni,1), f(nj,i,1), f(nj,ni,1));
     62                fx = bilinear_interpol(x/x2, y/y2, f(j,i,2), f(j,ni,2), f(nj,i,2), f(nj,ni,2));
     63
     64                % get new coordinates for this point
     65                %newx = i + x + fx;
     66                %newy = j + y + fy;
     67
     68                nx = fx;
     69                ny = fy;
     70
     71                % check borders
     72                if nx > imw, nx = imw; end
     73                if nx < 1, nx = 1; end
     74
     75                if ny > imh, ny = imh; end
     76                if ny < 1, ny = 1; end
     77
     78                f(j+y,i+x,[1,2]) = [ny,nx];
     79
     80                % do interpolation with new and old pixel
     81
     82                % vector from interpolated p to q
     83                vx = nx - (i+x);
     84                vy = ny - (j+y);
     85
     86                % float part of this vector and x2,y2 for interpolation
     87                if vx > 0
     88                    fpx = abs(vx - floor(vx)); x22 = i+x+1;
     89                elseif vx < 0
     90                    fpx = abs(vx - ceil(vx)); x22 = i+x-1;
     91                else
     92                    fpx = 0; x22 = i+x;
     93                end
     94
     95                if vy > 0
     96                    fpy = abs(vy - floor(vy)); y22 = j+y+1;
     97                elseif vy < 0
     98                    fpy = abs(vy - ceil(vy)); y22 = j+y-1;
     99                else
    100                    fpy = 0; y22 = j+y;
    101                end
    102
    103                x11 = i+x;
    104                y11 = j+y;
    105                f_nx = floor(nx);
    106                f_ny = floor(ny);
    107
    108
    109                for layer = 1:imd
    110                    % interpolate the pixel at x+x2+float,y+y2+float and
    111                    % move it to nx,ny
    112                    newImg(f_ny, f_nx, layer) = bilinear_interpol(...
    113                            fpx, fpy,...
    114                            oldImg(y11, x11,layer),...
    115                            oldImg(y11, x22,layer),...
    116                            oldImg(y22, x11,layer),...
    117                            oldImg(y22, x22,layer));
    118                end
    119            end
    120        end
    121    end
    122end
    123
    124
    125
    126end
    127

    128function f = bilinear_interpol(x,y,f00,f10,f01,f11)
    129% bilinear interpolation for 1*1 square, if you pass x=x/x2 and y=y/y2 you
    130% get bilinear interpolation for rectangle in begining of cordinate system
    131f = (f00*(1-x) + f10*x)*(1-y) + (f01*(1-x) + f11*x)*y;
    132end
    133