%% Known coordinate positions and corresponding data
x = [40,40,86,86,91,91,110,110,75,75,25,25,35,35,10,10,14,14];
y = [0,0,20,20,70,70,110,110,150,150,150,150,130,130,70,70,30,30];
z = [0,1.9,0.1,1.9,0.1,1.9,0.1,1.9,0.1,1.9,0.1,1.9,0.1,1.9,0.1,1.9,0.1,1.9];
%h = [3.9,6.2 ,   4,6.3,   4.3,6.5,  4.5,6.6,  4.6,  6.6,  4.6, 6.7, 4.5,  6.5, 4.6,  6.3,4.1,6.2];
h = [4.2,6.5,   4.3,6.6,   4.6,6.8,   4.8,6.9,   4.9,  6.9,  4.9,  7,   4.8,6.8,4.6,6.7,4.4,6.5];
%h=[  17.5,20.1, 17.6,19.9,   17.9,19.8,   17.7,20.1,      17.9,20.2,    17.5,20.1,       17.5,19.8,  17.2,20.2,  17.4,20.1  ];
%h=[   15.6,18,  15.5,17.9,   15.9,17.8,   15.6,17.9,      15.9,18.1,    15.5,17.8,       15.4,17.8,   15.1,18.2,  15.3,17.9  ];

%% Create a grid
[Y, X, Z] = meshgrid(0:1:150, 0:1:110, 0.1:0.1:1.9);

% Gaussian radial basis function parameters
epsilon = 0.0003; 
n = length(x);
h_pred = zeros(n, 1); % Initialization of predicted values
errors = zeros(n, 1); % Initialization error

% Calculate the h value over the entire grid
h_grid = zeros(size(X));

% Define boundary points
boundary_points = [
    25,150, 0;75,150,0; 25.5,149,0;76,149, 0; 34.5,131,0;94,131,0;  34,132,0;93,132,0; 33.5,133,0;92,133,0;  33,134,0;91,134,0;  32.5,135,0;90,135,0; 32,136,0;89,136,0; 31.5,137,0;88,137,0; 31,138,0;87,138,0;  30.5,139,0;86,139,0; 85,140,0;30,140,0; 29.5,141,0;84,141,0; 83,142,0;29,142,0; 28.5,143,0;82,143,0; 81,144,0;28,144,0; 80,145,0;27.5,145,0;  27,146,0; 79,146,0; 78,147,0;26.5,147,0; 77,148,0; 26,148,0;     
    35,130, 0; 95,130,0;  34,129,0;96,129, 0; 24.5,111,0;109,111,0;  25,112,0;108,112,0; 25.5,113,0;107,113,0;  26,114,0;106,114,0;  26.5,115,0;105,115,0; 27,116,0;104,116,0; 27.5,117,0;103.5,117,0; 28,118,0;103,118,0;  28.5,119,0;102.5,119,0; 102,120,0;29,120,0; 29.5,121,0;101.5,121,0; 101,122,0;30,122,0; 30.5,123,0;100.5,123,0; 100,124,0;31,124,0; 31.5,125,0;99.5,125,0;  32,126,0; 99,126,0; 32.5,127,0;98,127,0; 33,128,0; 97,128,0;   
    24,110, 0; 110,110,0; 23.5,109,0;109.5,109, 0; 14.5,91,0;96.5,91,0;  15,92,0;97,92,0; 15.5,93,0;97.5,93,0;  16,94,0;98,94,0;  16.5,95,0;98.5,95,0; 17,96,0;99,96,0; 17.5,97,0;99.5,97,0; 18,98,0;100,98,0;  18.5,99,0;100.5,99,0; 19,100,0;101,100,0; 102,101,0;19.5,101,0; 103,102,0;20,102,0; 20.5,103,0;104,103,0; 105,104,0;21,104,0; 21.5,105,0;106,105,0;  22,106,0; 107,106,0; 22.5,107,0;108,107,0; 23,108,0; 109,108,0; 
    14,90, 0; 96,90, 0; 13.96,89,0;95.5,89, 0; 10.25,71,0;86.5,71,0;  10.5,72,0;87,72,0; 10.75,73,0;87.5,73,0;  11,74,0;88,74,0;  11.25,75,0;88.5,75,0; 11.5,76,0;89,76,0; 11.75,77,0;89.5,77,0; 12,78,0;90,78,0;  12.25,79,0;90.5,79,0; 12.5,80,0;91,80,0; 12.75,81,0;91.5,81,0; 13,82,0;92,82,0; 13.25,83,0;92.5,83,0; 13.5,84,0;93,84,0; 13.75,85,0;93.5,85,0;  13.85,86,0; 94,86,0; 13.9,87,0;94.5,87,0; 13.93,88,0; 95,88,0; 
    10,70, 0; 86,70, 0; 9.9,69,0;85.5,69, 0; 8.1,51,0;76.5,51,0;  8.2,52,0;77,52,0; 8.3,53,0;77.5,53,0;  8.4,54,0;78,54,0;  8.5,55,0;78.5,55,0;  8.6,56,0;79,56,0; 8.7,57,0;79.5,57,0; 8.8,58,0;80,58,0;  8.9,59,0;80.5,59,0; 9,60,0;81,60,0; 9.1,61,0;81.5,61,0; 9.2,62,0;82,62,0; 9.3,63,0;82.5,63,0; 9.4,64,0;83,64,0; 9.5,65,0;83.5,65,0;  9.6,66,0; 84,66,0;  9.7,67,0;84.5,67,0; 9.8,68,0; 85,68,0; 
    8,50, 0; 76,50, 0; 8.1,49,0;76,49, 0; 13.6,31,0;70.3,31,0;  13.3,32,0;70.6,32,0; 12.9,33,0;70.9,33,0; 12.6,34,0;71.3,34,0;  12.3,35,0;71.6,35,0;  11.9,36,0;71.9,36,0;  11.6,37,0;72.3,37,0; 11.3,38,0;72.6,38,0;  10.9,39,0;72.9,39,0; 10.6,40,0;73.3,40,0;  10.3,41,0;73.6,41,0; 9.9,42,0;73.9,42,0; 9.6,43,0;74.3,43,0; 9.3,44,0;74.6,44,0; 8.9,45,0;74.9,45,0;  8.6,46,0; 75.6,46,0;  8.7,47,0;75.8,47,0; 8.2,48,0; 75.9,48,0; 
    14,30, 0; 70,30, 0; 14.5,29,0;69,29, 0; 22.9,11,0;60,11,0;  22.7,12,0;60.5,12,0; 22.3,13,0;61,13,0; 22,14,0;61.5,14,0;  21.5,15,0;62,15,0;  21,16,0;62.5,16,0;  20.5,17,0;63,17,0; 20,18,0;63.5,18,0;  19.5,19,0;64,19,0; 19,20,0;64.5,20,0;  18.5,21,0;65,21,0; 18,22,0;65.5,22,0; 17.5,23,0;66,23,0; 17,24,0;66.5,24,0; 16.5,25,0;67,25,0;  16,26,0;67.5,26,0;  15.5,27,0;68,27,0;  15,28,0;68.5,28,0; 
    23,10, 0; 60,10, 0;  25,9,0;57,9, 0;  39,0.6,0; 38,1.2,0;69,1,0; 37,1.8,0;36,2,0; 35,3,0;34,4,0; 39,4,0;  32,5,0;37,5,0;  30,6,0;35,6,0;  28.5,7,0;33,7,0;  27,8,0;31,8,0;
    40,0, 0;34,3.6,0;33,4.2,0;32,4.8,0;
    
    
    
    25,150, 2;75,150, 2; 25.5,149,2;76,149, 2; 34.5,131,2;94,131,2;  34,132,2;93,132,2; 33.5,133,2;92,133,2;  33,134,2;91,134,2;  32.5,135,2;90,135,2; 32,136,2;89,136,2; 31.5,137,2;88,137,2; 31,138,2;87,138,2;  30.5,139,2;86,139,2; 85,140,2;30,140,2; 29.5,141,2;84,141,2; 83,142,2;29,142,2; 28.5,143,2;82,143,2; 81,144,2;28,144,2; 80,145,2;27.5,145,2;  27,146,2; 79,146,2; 78,147,2;26.5,147,2; 77,148,2; 26,148,2;  
    35,130, 2; 95,130, 2;34,129,2;96,129, 2; 24.5,111,2;109,111,2;  25,112,2;108,112,2; 25.5,113,2;107,113,2;  26,114,2;106,114,2;  26.5,115,2;105,115,2; 27,116,2;104,116,2; 27.5,117,2;103.5,117,2; 28,118,2;103,118,2;  28.5,119,2;102.5,119,2; 102,120,2;29,120,2; 29.5,121,2;101.5,121,2; 101,122,2;30,122,2; 30.5,123,2;100.5,123,2; 100,124,2;31,124,2; 31.5,125,2;99.5,125,2;  32,126,2; 99,126,2; 32.5,127,2;98,127,2; 33,128,2; 97,128,2; 
    24,110, 2; 110,110, 2; 23.5,109,2;109.5,109, 2; 14.5,91,2;96.5,91,2;  15,92,2;97,92,2; 15.5,93,2;97.5,93,2;  16,94,2;98,94,2;  16.5,95,2;98.5,95,2; 17,96,2;99,96,2; 17.5,97,2;99.5,97,2; 18,98,2;100,98,2;  18.5,99,2;100.5,99,2; 19,100,2;101,100,2; 102,101,2;19.5,101,2; 103,102,2;20,102,2; 20.5,103,2;104,103,2; 105,104,2;21,104,2; 21.5,105,2;106,105,2;  22,106,2; 107,106,2; 22.5,107,2;108,107,2; 23,108,2; 109,108,2; 
    14,90, 2; 96,90, 2;13.96,89,2;95.5,89,2; 10.25,71,2;86.5,71,2;  10.5,72,2;87,72,2; 10.75,73,0;87.5,73,0;  11,74,2;88,74,2;  11.25,75,2;88.5,75,2; 11.5,76,2;89,76,2; 11.75,77,2;89.5,77,2; 12,78,2;90,78,2;  12.25,79,2;90.5,79,2; 12.5,80,2;91,80,2; 12.75,81,2;91.5,81,2; 13,82,2;92,82,2; 13.25,83,2;92.5,83,2; 13.5,84,2;93,84,2; 13.75,85,2;93.5,85,2;  13.85,86,2; 94,86,2; 13.9,87,2;94.5,87,2; 13.93,88,2; 95,88,2;
    10,70, 2; 86,70, 2; 9.9,69,2;85.5,69,2; 8.1,51,2;76.5,51,2;  8.2,52,2;77,52,2; 8.3,53,2;77.5,53,2;  8.4,54,2;78,54,2;  8.5,55,2;78.5,55,2;  8.6,56,2;79,56,2; 8.7,57,2;79.5,57,2; 8.8,58,2;80,58,2;  8.9,59,2;80.5,59,2; 9,60,2;81,60,2; 9.1,61,2;81.5,61,2; 9.2,62,2;82,62,2; 9.3,63,2;82.5,63,2; 9.4,64,2;83,64,2; 9.5,65,2;83.5,65,2;  9.6,66,2; 84,66,2;  9.7,67,2;84.5,67,2; 9.8,68,2; 85,68,2; 
    8,50, 2; 76,50, 2; 8.1,49,2;76,49, 2; 13.6,31,0;70.3,31,2;  13.3,32,2;70.6,32,2; 12.9,33,2;70.9,33,2; 12.6,34,2;71.3,34,2;  12.3,35,2;71.6,35,2;  11.9,36,2;71.9,36,2;  11.6,37,2;72.3,37,2; 11.3,38,2;72.6,38,2;  10.9,39,2;72.9,39,2; 10.6,40,2;73.3,40,2;  10.3,41,2;73.6,41,2; 9.9,42,2;73.9,42,2; 9.6,43,2;74.3,43,2; 9.3,44,2;74.6,44,2; 8.9,45,2;74.9,45,2;  8.6,46,2; 75.6,46,2;  8.7,47,2;75.8,47,2; 8.2,48,2; 75.9,48,2; 
    14,30, 2; 70,30, 2; 14.5,29,2;69,29,2; 22.9,11,2;60,11,2;  22.7,12,2;60.5,12,2; 22.3,13,2;61,13,2; 22,14,2;61.5,14,2;  21.5,15,2;62,15,2;  21,16,2;62.5,16,2;  20.5,17,2;63,17,2; 20,18,2;63.5,18,2;  19.5,19,2;64,19,2; 19,20,2;64.5,20,2;  18.5,21,2;65,21,2; 18,22,2;65.5,22,2; 17.5,23,2;66,23,2; 17,24,2;66.5,24,2; 16.5,25,2;67,25,2;  16,26,2;67.5,26,2;  15.5,27,2;68,27,2;  15,28,2;68.5,28,2; 
    23,10, 2; 60,10, 2;25,9,2;57,9,2;  39,0.6,2;38,1.2,2;69,1,2; 37,1.8,2; 36,2,2; 35,3,2; 34,4,2;39,4,2;   32,5,2;37,5,2; 30,6,2;35,6,2;  28.5,7,2;33,7,2;  27,8,2;31,8,2;
    40,0, 2;34,3.6,2;33,4.2,2;32,4.8,2
];

% LOOCV
for i = 1:n
    % Current test point
    xtest = x(i);
    ytest = y(i);
    ztest = z(i);
    htest = h(i);

    % Training data
    x_train = x;
    y_train = y;
    z_train = z;
    h_train = h;

    % Remove the test point
    x_train(i) = [];
    y_train(i) = [];
    z_train(i) = [];
    h_train(i) = [];

    % Calculate the RBF matrix for the training data
    distances_train = sqrt((x_train - x_train').^2 + (y_train - y_train').^2 + (z_train - z_train').^2);
    phi_train = rbf_gaussian(distances_train, epsilon);

    % Calculate the weights
    W = phi_train \ h_train';

    % Calculate the distance of the test point
    distances_test = sqrt((x_train - xtest).^2 + (y_train - ytest).^2 + (z_train - ztest).^2);
    phi_test = rbf_gaussian(distances_test, epsilon);

    % Predict the h-value of the test point
    h_pred(i) = phi_test * W;
    
  
    errors(i) = htest - h_pred(i);
end

for i = 1:numel(X)
    xtest = X(i);
    ytest = Y(i);
    ztest = Z(i);
    
    % 
    distances_test = sqrt((x_train - xtest).^2 + (y_train - ytest).^2 + (z_train - ztest).^2);
    phi_test = rbf_gaussian(distances_test, epsilon);

    % 
    h_grid(i) = phi_test * W;
end

% Constrain h_grid to within the boundary points' range
valid_mask = inpolygon(X(:), Y(:), boundary_points(:, 1), boundary_points(:, 2));
h_grid(~valid_mask) = NaN; % When outside the range, the value is set to  NaN

% Calculation Error Metrics
% MSE = mean(errors.^2);
% RMSE = sqrt(MSE);
% MAE = mean(abs(errors));

% Printing Error Index
% fprintf('MSE: %.4f\n', MSE);
% fprintf('RMSE: %.4f\n', RMSE);
% fprintf('MAE: %.4f\n', MAE);


figure;
% Plot a colored 3D graph
figure;
slice(Y,X, Z, h_grid, [], [], 0:0.05:2); % Select a slice along the Z-axis
colorbar;
xlabel('Length(m)');
ylabel('Width(m)');
zlabel('Water depth(m)');
title('DO concentration(mg/L)');
%title('Water Temperature concentration(℃)');
%colormap(jet);
shading interp;
view(3); 
colorbar; 
%xlabel('Length(m)');
%ylabel('Width(m)');
%zlabel('Water depth(m)');
%title('DO concentration(mg/L)');

colormap(parula);
shading interp;

function phi = rbf_gaussian(distances, epsilon)
    phi = exp(-epsilon * distances.^2);
end