function [MSTLength, dLength] = mst_length(rr,drr,gamma, qLevels, display)



%function [MSTLength, dLength] = mst_length(rr,drr, gamma, qLevels, disp)

%Finds the total EMST Length using Delaunay triangulation + Kruskal algorithm to compute the EMST

%rr is the feature matrix in the format : rr = [Feature_from1 Feature_from2]

%where Feature_i is a column and represents the i'th dimension 

%drr is the derivative (w.r.t. transformation parameters) of feaures in the format: 

%drr = [dFeauture_from2] <numOfSamples -by- Dimension of samples/2 * numOfTransParams>



%Output:    MSTLength: scalar <double>

%           dLength: 1 -by- numOfTransParams

%           

%Author: Mert Rory Sabuncu

% **********************************
% NEEDS:
% 1) sort_out_samples.dll (or mex) -- generated from sort_out_samples.c
% 2) computeEMST.dll (or mex) -- generated from computeEMST.c
% 3) derivativeEMST.dll (or mex) -- generated from derivativeEMST.c
% 4) plot_mst_2d.m -- to plot a 2D mst's
% **********************************

%Date: June 2005



if (nargin < 5)

    display =0; %Plot the EMST?

end;

if (nargin < 4)

    qLevels = 10;

end

if (nargin < 3)

    gamma = 0.3;

end



if (nargin == 1)

    drr = zeros(size(rr,1),1);

end

dim = size(rr,2);

der_dim = size(drr,2) * 2 / dim;

[qrr_unique, drr_unique] = sort_out_samples(rr, drr, qLevels);

if (sum(qrr_unique(:,1),1) == 0)

    MSTLength = 0;

    return;

end

TRI = delaunay(qrr_unique (:,1)', qrr_unique (:,2)');

EDGE = [TRI(:,1) TRI(:,2); TRI(:,1) TRI(:,3); TRI(:,2) TRI(:,3)];

EDGE = unique(EDGE,'rows');

[MSTLength, E] = computeEMST(qrr_unique', EDGE',gamma); %dll file written by Mert Rory Sabuncu

MSTLength = MSTLength/size(rr,1)^((size(rr,2)-gamma)/size(rr,2));

%E = [I_unique(E(1,:))'; I_unique(E(2,:))'];

if (display)

    h = figure();

    plot_mst_2d(qrr_unique,E',h);

    hold off;

end

[dLength] = derivativeEMST(qrr_unique', drr_unique', E, gamma); %dll file written by Mert Rory Sabuncu

return