% GriddataExample0306NEW.m  L Braile 03/05/06, revised 3/26/06
% Modified from example in Matlab Help, Griddata
% Pauses removed (3/26/06)
rand('seed',0) % set random number seed 
% -----------------------------------------------------------
% Find 100 random x,y locations between -2 to 2
% -----------------------------------------------------------
x = rand(100,1)*4-2; y = rand(100,1)*4-2;
% -----------------------------------------------------------
% Calculate function z(x,y) at x,y locations
% -----------------------------------------------------------
z = x.*exp(-x.^2-y.^2); % function definition
ti = [-2:.25:2]; % set grid parameters
[XI,YI] = meshgrid(ti,ti); % make x,y grid
% The gridding using nearest is to obtain values at and near
% the edges of the grid.  Near the edges, the gridding using
% cubic interpolation will result in Z(x,y) values will be NaN.
% Away from the edges, the interpolated values using the cubic
% interpolation method will be used.
ZI = griddata(x,y,z,XI,YI,'cubic'); % grid (2-D interpolate)
ZI2 = griddata(x,y,z,XI,YI,'nearest'); % grid (2-D interpolate)
% -----------------------------------------------------------
% The next section of the code replaces the Z(x,y) = NaN values
% from the cubic interpolation with values obtained from the
% nearest interpolation.
n = size(ZI);
nx = n(1,2);
ny = n(1,1);
for ix = 1:nx
    for iy = 1:ny
        iNaN = isnan(ZI(iy,ix));
        if iNaN == 1
            ZI(iy,ix) = ZI2(iy,ix);
        end
    end
end
% -----------------------------------------------------------
mesh(XI,YI,ZI), hold on % Make a mesh of gridded surface
plot3(x,y,z,'o'), hold off % Show original x,y,z data values
% -----------------------------------------------------------
% Plot color image and overlay contours
figure
pcolor(XI,YI,ZI)
shading interp
colormap jet
colorbar
hold on
% plot the x-y locations of the original data
plot(x,y,'+k')
% Contour the gridded data
v = [-0.5:0.1:0.5];
[c,h] = contour(XI,YI,ZI,v,'w-');
clabel(c,h,'color','w','fontweight','bold','fontsize',15)
axis square
hold off
%
% -----------------------------------------------------------
% use interp2 to regrid for increased resolution and smoothness
% -----------------------------------------------------------
figure
[X2,Y2] = meshgrid(-2:0.05:2); % set finer grid interval
Z2 = interp2(XI,YI,ZI,X2,Y2,'cubic');
% -----------------------------------------------------------
% replot color display and contour map for interpolated data
% -----------------------------------------------------------
pcolor(X2,Y2,Z2)
shading interp
colorbar
hold on
% -----------------------------------------------------------
% Contour the gridded data
% -----------------------------------------------------------
v = [-0.5:0.1:0.5];
[c,h] = contour(X2,Y2,Z2,v,'w-');
clabel(c,h,'color','w','fontweight','bold','fontsize',15)
axis square
hold off
% -----------------------------------------------------------
% Calculate and plot exact 2-D function for comparison
% -----------------------------------------------------------
t3 = -2:.05:2; % set grid parameters
[X3,Y3] = meshgrid(t3,t3); % make x,y grid
Z3 = X3.*exp(-X3.^2-Y3.^2); % function definition
figure
pcolor(X3,Y3,Z3)
shading interp
colormap jet
colorbar
hold on
% -----------------------------------------------------------
% Contour the gridded data
% -----------------------------------------------------------
v = [-0.5:0.1:0.5];
[c,h] = contour(X3,Y3,Z3,v,'w-');
clabel(c,h,'color','w','fontweight','bold','fontsize',15)
axis square
hold off
% -----------------------------------------------------------
% Calculate the residual surface (errors) and plot
% -----------------------------------------------------------
ZR = Z3 - Z2;
figure
pcolor(X3,Y3,ZR)
shading interp
colormap jet
colorbar
hold on
% plot the x-y locations of the original data
plot(x,y,'+k')
% -----------------------------------------------------------
% Contour the gridded data
% -----------------------------------------------------------
v = [-0.05 -0.02 -0.01 0.01 0.02 0.05];
[c,h] = contour(X3,Y3,ZR,v,'w-');
clabel(c,h,'color','w','fontweight','bold','fontsize',15)
axis square
set(gca,'fontsize',16,'linewidth',2)
xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
title('100 random x,y data, z = x.*exp(-x.^2-y.^2), residuals','fontsize',16)
hold off
% -----------------------------------------------------------