% ititialization of the FMM from given curves
% lambda is the stretch parameter
function T = process(curve1, curve2, lambda)

m = size(curve1, 1)-2;
n = size(curve2, 1)-2;

% compute curvatures k1, k2 (disregard endpoints)
k1 = zeros(m,1);
k2 = zeros(n,1);
for i=1:m,
    delX = curve1(i+2, 1) - curve1(i, 1);
    delY = curve1(i+2, 2) - curve1(i, 2);
    k1(i) = 4*( (curve1(i+2, 2) - 2*curve1(i+1, 2) + curve1(i, 2))*delX - ...
                (curve1(i+2, 1) - 2*curve1(i+1, 1) + curve1(i, 1))*delY ) /...
                (delX*delX + delY*delY)^(3/2);
%   k1(i) = norm(curve1(i, 1:2) - 2*curve1(i+1, 1:2) + curve1(i+2, 1:2));
end

for i=1:n,
    delX = curve2(i+2, 1) - curve2(i, 1);
    delY = curve2(i+2, 2) - curve2(i, 2);
    k2(i) = 4*( (curve2(i+2, 2) - 2*curve2(i+1, 2) + curve2(i, 2))*delX - ...
                (curve2(i+2, 1) - 2*curve2(i+1, 1) + curve2(i, 1))*delY ) /...
                (delX*delX + delY*delY)^(3/2);
%   k2(i) = norm(curve2(i, 1:2) - 2*curve2(i+1, 1:2) + curve2(i+2, 1:2));
end

% compute the speed function F
T = [];
F = [];
dtds = m/n;
for i=1:m,
   for j=1:n,
      F(i,j) = abs(k2(j)*dtds - k1(i)) + lambda*abs(dtds - 1);
%      if F(i,j) == 0
%         F(i,j) = 1;
%      end

   end
end

contour(F);

% start point
initSet = [1 1];

% dispatch
T = FMM(F, initSet);