function newM = Interpolate2D_nv(M)
% Perform 2-d interpolation on the real-valued data in nr-by-nc matrix M.
% The interpolated data are added between existing data points so newM is 
%   (2*nr-1)-by-(2*nc-1). Use the simple average as the interpolated value.
%
% NON-VECTORIZED

[nr,nc]= size(M);

% Interpolate one dimension at a time.

% New dimensions after interpolating in x-direction: nr-by-(nc*2-1)
for r= 1:nr
  for c= 1:nc-1
    wideIm(r,c*2-1)= M(r,c);
    wideIm(r,c*2)= (M(r,c)+M(r,c+1))/2;
  end
  wideIm(r,nc*2-1)= M(r,nc);
end

% New dimensions after interpolating in y-direction: (nr*2-c)-by-(nc*2-1)
for c= 1:nc
  for r= 1:nr-1
      newM(r*2-1,c)= wideIm(r,c);
      newM(r*2,c)= (wideIm(r,c)+wideIm(r+1,c))/2;
  end
  newM(nr*2-1,c)= wideIm(nr,c);
end



