Interpolation of sparse pointwise and integral geographical Data over Time/spline2.m

function splinedata = spline2(lattice, points, volumes, time, eta)

% Gives the "best fitting" spline in the sense of least squares optimization.
% 
% lattice has two parts, lattice.points is an array of points (rowwise)
%                        lattice.cells  is an array of cells  (rowwise pointids)
% points contains the prescribed point values in the form 
%                      (pointid, value, timeofmeasurement)
% volumes contains the prescribed volume values in the form 
%                      (cellid,  value, timeofmeasurement)
% time is the currecnt time 
% eta is the speed of decay of the weighting function \phi
%
% We return a system matrix in splinedata.A and a right hand side in 
% splinedata.b, and a weight vector, splinedata.weight

% Assembling of 2d-degree1-spline matrix: 

weightfunction = @(t) (t>=0).*exp(-eta*t);
% kill "future" and irrelevant values": 
idx = find(weightfunction(time-points(:,3))>0);
points = points(idx, :);
idx = find(weightfunction(time-volumes(:,3))>0);
volumes = volumes(idx,:);
 
numberOfEquations = size(points, 1) + size(volumes, 1);
numberOfVariables = size(lattice.points,1);
 
splinedata.A = sparse(numberOfEquations, numberOfVariables);
splinedata.b = zeros(numberOfEquations, 1);
splinedata.weight = zeros(numberOfEquations, 1);

 
% for each point value, we have an equation of the form 
%            x_{pointid} = value

%for i = 1:size(points, 1)
%  splinedata.A(i, points(i,1)) = 1;
  % "bigger" basis functions 
%  splinedata.A(i, neighbours(lattice, points(i,1))) = 1/2;
%  splinedata.b(i)              = points(i,2);
%  splinedata.weight(i)         = weightfunction(time - points(i,3));
%end;

splinedata.A(1:size(points,1), points(:,1)) = eye(size(points,1));
splinedata.b(1:size(points,1)) = points(:,2);
splinedata.weight(1:size(points,1)) = weightfunction(time - points(:,3));



% for each cell value, we have an equation of the form 
%        (sum of x_{points})/3 *volume = value

% hence compute the cellvolumes 
p1 = lattice.points(lattice.cells(:,1),:);
p2 = lattice.points(lattice.cells(:,2),:);
p3 = lattice.points(lattice.cells(:,3),:);

vols = 1/2*abs((p1(:,1)-p3(:,1)).*(p2(:,2)-p1(:,2)) ...
            -  (p1(:,1)-p2(:,1)).*(p3(:,2)-p1(:,2)));

p = size(points, 1);
for i = 1:size(volumes, 1)
  cell              = volumes(i, 1);
  splinedata.A(p+i, lattice.cells(cell, :)) = 1/3*vols(cell);
  splinedata.b(p+i)                         = volumes(i, 2);
  splinedata.weight(p+i) = weightfunction(time - volumes(i,3));
end;
 
%%% THE FOLLOWING ENSURES smoothness, but is a bit too hard
 
%% for each point, we want its value to be the mean 
%% of its neighbour points
p = size(points,1) + size(volumes, 1);
for i = 1:size(lattice.points,1)
  nbhd = neighbours(lattice, i);
  splinedata.A(p+i, nbhd)        = 1;
  splinedata.A(p+i,    i)        = -length(nbhd);
  splinedata.b(p+i)              = 0;
  splinedata.weight(p+i)         = 0.1;  
end;