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

function splinetest2() 

% tests the spline thing ....
MAXTIME = 10;
Np      = 4000;
Nv      = 1000;
timestep = 0.1;

% Generate a simple 2d-lattice
n      = 20;
pointsx = linspace(0,1,n)';
pointsy = linspace(0,1,n)';
[X,Y]   = meshgrid(pointsx,pointsy);
points  = [X(:),Y(:)];

lowerleft = (1:n^2-n)';
lowerleft = lowerleft(mod(lowerleft, n) ~= 0);
cells = [lowerleft, lowerleft+1, lowerleft+n; ...
        lowerleft+1, lowerleft+n+1, lowerleft+n];

lattice = assemble_lattice(points, cells);

% sample a function to be "interpolated as test" 
f     = @(x,t) sin(2*pi*(x(:,1)+t)) + cos(pi*x(:,2)+t);

% generate random point data:
pointid = randi(size(lattice.points, 1),Np,1);
times   = max(0,1.1*MAXTIME*rand(Np, 1)-0.1*MAXTIME);

points  = [pointid,...
           f(lattice.points(pointid,:), times)+0.1*randn(Np,1),...
           times];
           
% generate random volume data:
cellid  = randi(size(lattice.cells, 1),Nv,1);
times   = MAXTIME*rand(Nv, 1);
volumes = [cellid, ...
           (f(lattice.points(lattice.cells(cellid,1),:), times)+0.1*randn(Nv,1))*2/((n-1)^2), ...
           times];

figure(1); clf;
p = patch('FaceColor', [1, 0, 1], 'FaceAlpha', 0.4);
%hold on;
%[~,c] = contourf(X,Y,zeros(size(X)),'zLevel', -3.5, 'zLevelMode', 'manual');
%hold off;

axis([0, 1, 0, 1, -3.5, 3.5]);
view([64, 15]);
set(gca, 'cLimMode', 'manual', 'cLim', [-3.5,3.5]);

i = 0;
for t = 0:timestep:MAXTIME
  i = i+1; % number of picture
  splinedata = spline2(lattice, points, volumes, t);
  
  % solve the system
  % (Ax-b)'D_w(Ax-b) -> min <=> A'D_wAx = A'D_wb
  % Regularization:
  mu = 1e-12;
  x = (splinedata.A'*diag(splinedata.weight)*splinedata.A...
         + mu*eye(size(splinedata.A,2))) \ ...
       splinedata.A'*(splinedata.weight.*splinedata.b);
  % compute point values from coefficients:
  
  set(p, 'Vertices', [lattice.points, x], 'Faces', lattice.cells);
  % plot contour: 
%  delete(c);
%  hold on;
%  [~,c] = contourf(X,Y,reshape(x,n,n),'zLevel', -3.5, 'zLevelMode', 'manual');
%  hold off;
  set(gca, 'cLimMode', 'manual', 'cLim', [-3.5,3.5]);
%  hold on;
%  for i = 50:-1:0
%     idx = ((volumes(:,3)<=(t-i*timestep)) & (volumes(:,3)>t-(i+1)*timestep));
%     plot(lattice.points(lattice.cells(volumes(idx,1),1)),(n-1)*volumes(idx,2), 
%          'Color', [i/50,i/50,1], 'Marker', '+', 'LineStyle', 'none');
%     plot(lattice.points(1+lattice.cells(volumes(idx,1),1)),(n-1)*volumes(idx,2), 
%          'Color', [i/50,i/50,1], 'Marker', '+', 'LineStyle', 'none');     
%  end;
%  for i = 50:-1:0
%     idx = find(points(:,3)<=(t-i*timestep) & points(:,3)>t-(i+1)*timestep);
%     plot(lattice.points(points(idx,1)), points(idx,2),
%     'Color',[1,i/50,i/50],'Marker','*','LineStyle', 'none');
%  end;
%  plot(lattice.points, x, 'm-');
%  hold off;
%  axis([0, 1, -2.5, 2.5]);
  title(sprintf('t=%2.1f', t));
%  print(sprintf('images/pic%d.pdf', i), "-dpdf"); 
  sleep(0.01);
end;