function [zTilde] = f(x,y,data, order ,h);
% Locally weighted polynomial regression
% Tri-cube weight-function
% Stig B. Mortensen, 2006
% 
% zTilde = f(x,y,data, order ,h)
% 
% Input:
% Data:   data(i, [x y z])
% x,y:    vectors of coordinates for evaluation
% order:  order of fitted polynomium
% h:      bandwidth 0<h<1
%



% number of nearest neighbors
NN = round(h*length(data));

dataSc = data(:,1:2);

for p=1:length(x)

    % find distance from current point x,y to data
    dist = sqrt((dataSc(:,1)-x(p)).^2 + (dataSc(:,2)-y(p)).^2);
    [sortDist sortIdx] = sort(dist);
    
    %extract nearest neighbors points
    sortIdx = sortIdx(1:NN);
    distNN = dist(sortIdx); 
    distMax(p) = max(distNN); %store size of local neighborhood
    xNN = dataSc(sortIdx,1);
    yNN = dataSc(sortIdx,2);
    zNN = data(sortIdx,3);
    
    % designmatrice og eval punkt
    U = ones(NN,1);
    ep = 1;
    switch order
        case 1
            U = [U xNN yNN];
            ep = [ep x(p) y(p)];
        case 2
            U = [U xNN yNN xNN.^2 yNN.^2 xNN.*yNN];
            ep = [ep x(p) y(p) x(p)^2 y(p)^2 x(p)*y(p)];
    end
    
    % weights tri-cube
    distNN = distNN/max(distNN); %truncate distNN to [0,1] for proper wheigths
    W = (1-distNN.^3).^3; 
    W = spdiags(W,0,length(W),length(W));
    
    % solve normal equation
    % theta = (U'*W*U)^(-1)*U'*W*zNN;
    theta = (U'*W*U)\(U'*W*zNN);
    
    % evaluate estimated locally wheighted polynomial
    zTilde(p) = ep*theta; %theta(1); %ep*theta;

    
end