% ititialization of the FMM from given curves using shape context histograms and chi-squared 
% histogram distance measure
%function T = process(curve1, curve2, lambda)

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

disp('Computing the shape contexts...');

NUM_BINS_DIST = 5;
NUM_BINS_ORIENT = 12;
h1 = zeros(NUM_BINS_DIST, NUM_BINS_ORIENT, m);
h2 = zeros(NUM_BINS_DIST, NUM_BINS_ORIENT, n);
% curve 1
% compute the maximum distance between a pair of points
maxDist = 0;
for i=1:m,
	for j=1:m,
        distVec = curve1(j,:) - curve1(i,:);
        distVecLen = norm(distVec);
        if distVecLen > maxDist
            maxDist = distVecLen;
        end
	end
end
logBase = maxDist^(1/(NUM_BINS_DIST-1));
denom = log(logBase);
TWOPI = 6.2832;
SMALL = 1e-6;
ORIENT_BIN_SIZE = TWOPI/NUM_BINS_ORIENT;
for i=1:m,
	for j=1:m,
        if i == j
            % i-th point is not in its own histogram
            continue;
        end
        % distance bin
        distVec = curve1(j,:) - curve1(i,:);
        distVecLen = norm(distVec);
        if (distVecLen - maxDist) < SMALL
            % farthest point - count in last bin
            distBin = NUM_BINS_DIST;
        elseif distVecLen < 1
            distBin = 1;
        else
            distBin = 2 + floor(log(distVecLen)/denom);
        end
        % orientation bin
        angleVec = atan2(distVec(2), distVec(1));
        if angleVec < 0
            angleVec = TWOPI - abs(angleVec);
        end
        orientBin = floor(angleVec/ORIENT_BIN_SIZE) + 1;
        % count the point in its bin
        h1(distBin, orientBin, i) = h1(distBin, orientBin, i) + 1;
	end
end
% normalize the histogram
h1 = h1/(m-1);
% curve 2
% compute the maximum distance between a pair of points
maxDist = 0;
for i=1:n,
	for j=1:n,
        distVec = curve2(j,:) - curve2(i,:);
        distVecLen = norm(distVec);
        if distVecLen > maxDist
            maxDist = distVecLen;
        end
	end
end
logBase = maxDist^(1/(NUM_BINS_DIST-1));
denom = log(logBase);
TWOPI = 6.2832;
ORIENT_BIN_SIZE = TWOPI/NUM_BINS_ORIENT;
for i=1:n,
	for j=1:n,
        if i == j
            % i-th point is not in its own histogram
            continue;
        end
        % distance bin
        distVec = curve2(j,:) - curve2(i,:);
        distVecLen = norm(distVec);
        if (distVecLen - maxDist) < SMALL
            % farthest point - count in last bin
            distBin = NUM_BINS_DIST;
        elseif distVecLen < 1
            distBin = 1;
        else
            distBin = 2 + floor(log(distVecLen)/denom);
        end
        % orientation bin
        angleVec = atan2(distVec(2), distVec(1));
        if angleVec < 0
            angleVec = TWOPI - abs(angleVec);
        end
        orientBin = floor(angleVec/ORIENT_BIN_SIZE) + 1;
        % count the point in its bin
        h2(distBin, orientBin, i) = h2(distBin, orientBin, i) + 1;
	end
end
% normalize the histogram
h2 = h2/(n-1);

disp('Computing the cost function F...'); 

% compare histograms of i and j using chi-squared test statistic
F = zeros(m,n);
SMALL = 0.0001;
for i=1:m,
    for j=1:n,
        for k=1:NUM_BINS_DIST,
            for l=1:NUM_BINS_ORIENT,
                histSum = (h1(k,l,i) + h2(k,l,j));
                if histSum ~= 0
                    histDist = h1(k,l,i) - h2(k,l,j);
                    F(i,j) = F(i,j) + histDist*histDist/histSum;
                end
            end
        end 
    end
end
F = F/2;

% adaptive lambda
valLambda = sum(sum(F))/(m*n);
F = F + valLambda;

% assumes closed curves as input - duplicate curve2 points along x-axis
F = [F F; F F; F F];

disp('Computing the weighted distance function T...');

% compute grid spacing array for each curve
spacing1 = zeros(m, 1);
spacing2 = zeros(n, 1);
spacing1(1) = norm(curve1(1,:) - curve1(m,:));
for i = 2:m,
    spacing1(i) = norm(curve1(i,:) - curve1(i-1,:));
end
spacing1 = [spacing1; spacing1; spacing1];
spacing2(1) = norm(curve2(1,:) - curve2(n,:));
for j = 2:n,
    spacing2(j) = norm(curve2(j,:) - curve2(j-1,:));
end
spacing2 = [spacing2; spacing2];

% start point
initSet = [1 1];
% begin 
T = FMM(initSet, F);

