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

m = size(curve1, 1)-2;
n = size(curve2, 1)-2;
CURVATURE_THRESHOLD = 0.08;
SMALL = 0.000001;

% compute curvatures k1, k2 (disregard endpoints)
k1 = zeros(m,1);
k2 = zeros(n,1);
for i=1:m,
   k1(i) = norm(curve1(i, 1:2) - 2*curve1(i+1, 1:2) + curve1(i+2, 1:2));
   if k1(i) > CURVATURE_THRESHOLD
       k1(i) = 0;
   end
end

for i=1:n,
   k2(i) = norm(curve2(i, 1:2) - 2*curve2(i+1, 1:2) + curve2(i+2, 1:2));
   if k2(i) > CURVATURE_THRESHOLD
       k2(i) = 0;
   end
end

% compute the speed function F
T = [];
F = [];
SMALL = 0.001;
dtds = abs(m - n)/(m + n);
for i=1:m,
   for j=1:n,
       if (k2(j) + k1(i)) > SMALL,
           F(i,j) = abs(k2(j) - k1(i))/(k2(j) + k1(i)) + lambda*dtds;
       else
           % both curvatures are zero
           F(i,j) = dtds;
       end
%      F(i,j) = abs(k2(j)*dtds - k1(i)) + lambda*abs(dtds - 1);
   end
end

contour(F);

% start point
initSet = [1 1];

% dispatch
T = FMM(F, initSet);

% compute level set normals on T 
[dTdx dTdy] = gradient(T);
% update F adding the stretching term based on normals
for i=1:m,
   for j=1:n,
       normVec = [dTdx(i,j) dTdy(i,j)];
       F(i,j) = F(i,j) + 0.005*abs(normVec(2)/normVec(1) - 1);
   end
end

% recompute T with updated F
T = FMM(F, initSet);

