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;