% eda10_05
% 2D interpolation of simulated pressure data
% this script
%             creates simulated data (Figue 1)
%             creates Delaunay triangle mesh (Figure 2)
%             finds triangle enclosong an arbitrary point (Figure 3)
%             linear interpolation (Figure 4)
%             cubic spline interpolation
%
% Note: this script used the griddata() function, which the MatLab
% documentation says is no longer recommended.  See eda08_05 for an
% alternative.

% rows of grid are x
I = 41;
Dx = 1;
xmin = 0;
xmax = 40;
x = xmin + (xmax-xmin)*[0:I-1]'/(I-1);

% columns of grid are y
J = 41;
Dy = 1;
ymin = 0;
ymax = 40;
y = ymin + (ymax-ymin)*[0:J-1]'/(J-1);

% grid 10% populated with data
N0 = floor(I*J/10);
rowindex0=unidrnd(I,N0,1);
colindex0=unidrnd(J,N0,1);

% remove duplicate points, because the interpolation
% function griddata() doesn't like them, by checking
% eack point against a matrix of ones and zeros
A = zeros(I,J);
rowindex=zeros(N0,1);
colindex=zeros(N0,1);
N=0;
for i = [1:N0]
    if ( A(rowindex0(i),colindex0(i)) == 0 )
        N=N+1;
        rowindex(N)=rowindex0(i);
        colindex(N)=colindex0(i);
        end
    A(rowindex0(i),colindex0(i)) = 1;
end
rowindex=rowindex(1:N);
colindex=colindex(1:N);
xobs = x(rowindex);
yobs = y(colindex);
N=length(xobs);

% create simulated pressure data (solution to Laplace's equation)
kappa = 0.1;
dtrue = 10*sin(kappa*xobs).*exp(-kappa*yobs);
sigmad = 0.1;
dobs = dtrue + random('Normal',0.0,sigmad,N,1);

% populate matrix A with observations
A = zeros(I,J);
for i = [1:N]
    A(rowindex(i),colindex(i)) = dobs(i);
end

% plot data
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(A))-min(min(A));
if( range==0 )
    range=1;
end
imagesc( [0, I], [0, J], (A-min(min(A)))/range);
hold on;
xlabel('y');
ylabel('x');
title('data');

% plot data and triangulation
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(A))-min(min(A));
if( range==0 )
    range=1;
end
imagesc( [0, I], [0, J], (A-min(min(A)))/range);
hold on;
xlabel('y');
ylabel('x');
title('data + triangulation');

% triangulation
mytri = DelaunayTri(xobs,yobs);

% plot triangles
XY=mytri.X;
[NXY, i] = size(XY);
TRI=mytri.Triangulation;
[NTRI, i] = size(TRI);
for i = [1:NTRI]
    v1=TRI(i,1); v2=TRI(i,2); v3=TRI(i,3);
    plot( [XY(v1,2), XY(v2,2)], [XY(v1,1), XY(v2,1)], 'k-', 'LineWidth', 2);
    plot( [XY(v2,2), XY(v3,2)], [XY(v2,1), XY(v3,1)], 'k-', 'LineWidth', 2);
    plot( [XY(v3,2), XY(v1,2)], [XY(v3,1), XY(v1,1)], 'k-', 'LineWidth', 2);
end

% find triangle encl;osing exemplary point and plot them
x0=[20,30];
tri0 = pointLocation(mytri, x0);
plot( x0(2), x0(1), 'ko', 'LineWidth', 2);
if( ~isnan(tri0) )
    v1=TRI(tri0,1); v2=TRI(tri0,2); v3=TRI(tri0,3);
    plot( [XY(v1,2), XY(v2,2)], [XY(v1,1), XY(v2,1)], 'k-', 'LineWidth', 3);
    plot( [XY(v2,2), XY(v3,2)], [XY(v2,1), XY(v3,1)], 'k-', 'LineWidth', 3);
    plot( [XY(v3,2), XY(v1,2)], [XY(v3,1), XY(v1,1)], 'k-', 'LineWidth', 3);
end

% unfold matrix into a vector
xp = zeros(I*J,1);
yp = zeros(I*J,1);
dp = zeros(I*J,1);
for i = [1:I]
for j = [1:J]
    k = j + J*(i-1);
    xp(k) = x(i);
    yp(k) = y(j);
end
end

% interpolate data linearly
dp=griddata(xobs,yobs,dobs,xp,yp,'linear');

B = zeros(I,J);
for i = [1:I]
for j = [1:J]
    k = j + J*(i-1);
    B(i,j) = dp(k);
end
end

% plot interpolated data
figure(3);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(B))-min(min(B));
if( range==0 )
    range=1;
end
imagesc( [0, I], [0, J], (B-min(min(B)))/range);
hold on;
xlabel('y');
ylabel('x');
title('linear interpolation');

% interpolate data linearly
dpp=griddata(xobs,yobs,dobs,xp,yp,'cubic');

C = zeros(I,J);
for i = [1:I]
for j = [1:J]
    k = j + J*(i-1);
    C(i,j) = dpp(k);
end
end

% plot interpolated data
figure(4);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
axis([0, I, 0, J]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(C))-min(min(C));
if( range==0 )
    range=1;
end
imagesc( [0, I], [0, J], (C-min(min(C)))/range);
hold on;
xlabel('y');
ylabel('x');
title('cubic interpolation');

