% GriddataExample0305NEW.m  L Braile 
% Modified from example in Matlab Help, Griddata
%
rand('seed',0)
% 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, grid spacing  = 0.25
[XI,YI] = meshgrid(ti,ti); % make x,y grid
ZI = griddata(x,y,z,XI,YI,'cubic'); % grid (2-D interpolate)
% This section interpolates values at the edges of the map area 
% and within the map area where data points are sparse and the
% griddata code returns 'NaN' values.  This section assumes 
% that the grid area is slightly larger in each dimension than
% the x,y space of the input data.
%
% Find number of NaN values not on outer edge
[nx ny] = size(ZI); % nx is num of rows, ny is num of col
ny1 = ny-1;  nx1 = nx-1;
ic = 0;
for ix = 2:nx1
    for iy = 2:ny1
        iNaN = isnan(ZI(ix,iy));
        if iNaN == 1
            ic = ic +1;
        else
        end
    end
end
d = zeros(ic,1);
inx = zeros(ic,1);
iny = zeros(ic,1);
xm = mean(XI(:));
ym = mean(YI(:));
ic = 0;
for ix = 2:nx1
    for iy = 2:ny1
        iNaN = isnan(ZI(ix,iy));
        if iNaN == 1
            ic = ic +1;
            d(ic,1) = sqrt((XI(ix,iy) - xm)^2 + (YI(ix,iy) - ym)^2);
            inx(ic,1) = ix;
            iny(ic,1) = iy;
        else
        end
    end
end
A = [d inx iny];
% sort A so that NaN values are interpolated beginning with those
% closest to the center of the map and working outward
B = sortrows(A,1);
inx = B(:,2);
iny = B(:,3);
for i = 1:ic
    it = isnan(ZI(inx(i),:));
    xt = [];
    yt = [];
    zt = [];
    for iy = 2:ny1
        if it(iy) == 0
            xt = [xt XI(inx(i),iy)];
            zt = [zt ZI(inx(i),iy)];
        else
        end
    end
    ZI(inx(i),iny(i)) = interp1(xt,zt,XI(inx(i),iny(i)),'spline','extrap');
end
% Extrapolate to outer edge of map
% left and right edge
for ix = 2:nx1
    ZI(ix,1) = interp1(XI(ix,2:ny1),ZI(ix,2:ny1),XI(ix,1),'nearest','extrap');
    ZI(ix,ny) = interp1(XI(ix,2:ny1),ZI(ix,2:ny1),XI(ix,ny),'nearest','extrap');
end
% top and bottom and corners
for iy = 1:ny
    ZI(1,iy) = interp1(YI(2:nx1,iy),ZI(2:nx1,iy),YI(1,iy),'nearest','extrap');
    ZI(nx,iy) = interp1(YI(2:nx1,iy),ZI(2:nx1,iy),YI(ny,iy),'nearest','extrap');
end
%
figure
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
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), Mesh','fontsize',16)
%
% 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
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), coarse grid','fontsize',16)
hold off
%
% use interp2 to regrid for increased resolution and smoothness
figure
[X2,Y2] = meshgrid(-2:0.05:2); % set finer grid interval, grid = 0.05
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
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), fine grid','fontsize',16)
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
set(gca,'fontsize',16,'linewidth',2)
xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
title('z = x.*exp(-x.^2-y.^2), exact function','fontsize',16)
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