% function [clustID] = cutDendrogram(den, data, threshold, startingLevel)
% cutDendrogram cuts dendrogram den based on the threshold parameter, 
% to extract the corresponding clustering solution. Here a method based 
% on the average intra-cluster standard deviation is used: this value is 
% computed at each level of the dendrogram starting from the leaves, 
% averaging over all existent clusters for that level. When the computed
% value exceeds the set threshold, the process stops, and the 
% clustering solution at the previous level of the dendrogram is returned.
% The function also requires the original data that have been clustered,
% and the level (startingLevel) at which to start climbing the dendrogram
% (= 0 when starting from the leaves).
% clustID lists, for each data point, the ID of the cluster containing that
% point.
%
% This file is part of Cluster4MetaAnalysis.
%     Cluster4MetaAnalysis is free software: you can redistribute it 
%     and/or modify it under the terms of the GNU General Public License 
%     as published by the Free Software Foundation, either version 3 of 
%     the License, or (at your option) any later version.
% 
%     Cluster4MetaAnalysis is distributed in the hope that it will be 
%     useful, but WITHOUT ANY WARRANTY; without even the implied warranty 
%     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with Cluster4MetaAnalysis.  If not, see
%     <http://www.gnu.org/licenses/>.


function [clustID] = cutDendrogram(den, data, threshold, startingLevel)

% Try all the cuts until the stopping criterion is met...
% If startingLevel = 0, or it is greater than the total number of data 
% points, let us just start from the lowest level. 
% Otherwise, let us take it as a starting level for cutting (for speed sake)

% number of data points
N = size(data, 1);

if (startingLevel == 0 || startingLevel >= N)
    level = N - 1;          % starting from N - 1 clusters
else
    level = startingLevel;
end

% preparing the data
merge = zeros(N, 4);
merge(:, 1:3) = data;

avgStdX = 0;
avgStdY = 0;
avgStdZ = 0;

% climbing the dendrogram to compute average standard deviation at each
% level, until the stopping criterion is met
while (avgStdX <= threshold && avgStdY <= threshold && avgStdZ <= threshold)
    
    if (level == 0)
        disp('Warning: the threshold you set might be too high! Returning one cluster.');
        clustID = cluster(den, 'MaxClust', 1);
        return;
    else
        avgStdX = 0;
        avgStdY = 0;
        avgStdZ = 0;
        % find clusters
        cl = cluster(den, 'MaxClust', level);
        merge(:, 4) = cl;

        % Compute average standard deviations over all clusters at this
        % level
        for i=1:level
            avgStdX = avgStdX + std(merge(find(merge(:, 4) == i), 1));
            avgStdY = avgStdY + std(merge(find(merge(:, 4) == i), 2));
            avgStdZ = avgStdZ + std(merge(find(merge(:, 4) == i), 3));
        end
        avgStdX = avgStdX/level;
        avgStdY = avgStdY/level;
        avgStdZ = avgStdZ/level;

        % next level
        level = level - 1;
    end

end
level = level + 1;


% last cycle execution made us go over threshold... let us take the
% clustering solution at the previous level.
level = level + 1;
clustID = cluster(den, 'MaxClust', level);














