% NOTE1: This code is not efficient. You shall Vectorize it easily
% and get rid of the for-loops, yielding a substantial boost in efficiency.
% NOTE2: boundary rows/columns may be left untouched which needs a simple fix
% The transform that you want to appply is (x,y)=T(n,m)
% The transform that we will use is its inverse (n,m)=invT(x,y) (you shall invert your transform)
% Original image I1 is of size S1xS2, transformed image I2 is of size N1xN2
for x=1:N1    % let N1 = S1 and N2= S2 for this example WAVE1
  for y=1:N2

     % s1- compute back coordinates in the original image domain
     % [n m] = invT(x,y);
     m = y;
     n = x - 20*sin(2*pi*y/128);

     % s2-s3 Apply (bilinear) interpolation to find the image intensity, if
     % n,m does not fall into integer indices
     d1 = n - floor(n);
     d2 = m - floor(m);

     % s3- define the image pixel
    if ( ((n>=1)&(n<S(1))) & ((m>=1)&(m<S(2))) )
     I2(x,y) = (1-d1)*(1-d2) * I1(floor(n),floor(m)) + ...
               (1-d1)*d2 * I1(floor(n),floor(m)+1) + ...
               d1*(1-d2) * I1(floor(n)+1,floor(m)) + ...
               d1*d2 * I1(floor(n)+1,floor(m)+1) ;
    end%if

  end
end