Mauch MIREX 2010

function song = segmentation_auto(song, param, varargin)

%% calculating similarity matrix

% exp of minus the cosine distance

if any(strcmp(fieldnames(song.syncchroma),'wide'))

	v = song.syncchroma.wide;

else

	v = song.syncchroma.treble + song.syncchroma.bass(:,1:12) * 0.8;

end

nBeat = song.nBeat;

flatPCP = std(v')<.000001;

v(flatPCP,1) = 1000; % arbitrary

simmat0 = squareform(pdist(v,'correlation'));

% simmat0(flatPCP,:) = 2;

% simmat0(:,flatPCP) = 2;

% simmat = qnormalise(exp(-simmat0),inf,1);

% simmat(flatPCP,:) = 0;

% simmat(:,flatPCP) = 0;

simmat0(flatPCP,:) = 1;

simmat0(:,flatPCP) = 1;

simmat = 1-simmat0/2;

simmat(flatPCP,:) = 0;

simmat(:,flatPCP) = 0;

%%

median_simmat = zeros(size(simmat));

for  i = 1:nBeat

    temp1 = diag(medfilt1(diag(simmat,i-1),param.seg.medfilt_length,[],1),i-1);

    median_simmat = median_simmat + temp1;

end

median_simmat = median_simmat + median_simmat' - diag(diag(median_simmat));

median_simmat(isnan(median_simmat)) = 0;

%%median_simmat = (median_simmat - mean(median_simmat(:)))/std(median_simmat(:));

%%%%%%%% ATTENTION!!!!!!!! using MEDIAN %%%%%%%%%%%%%

% median_simmat = (median_simmat - median(median_simmat(:)))/std(median_simmat(:));

if param.seg.standardise

    med_median_simmat = repmat(median(median_simmat),nBeat,1);

    std_median_simmat = repmat(std(median_simmat),nBeat,1);

    median_simmat = (median_simmat - med_median_simmat) ./ std_median_simmat;

end

%% get bar boundaries!!!

potential_duplicates = triu(median_simmat > param.seg.thresh_beat);

%%

partlengths = param.seg.minlength:4:param.seg.maxlength;

nPartlengths = length(partlengths);

% initialise arrays

simArray = zeros(nBeat,nBeat,nPartlengths);

decisionArray2 = zeros(nBeat,nBeat,nPartlengths);

fprintf(1,'\n');

fprintf(2,'    ');

for iLength = 1:nPartlengths

    len = partlengths(iLength);

    nUsedBeat = nBeat - len + 1; % number of potential rep beginnings: they can't overlap at the end of the song

    fprintf(1,'\rLength: %1.0d',len)

    fprintf(2,'\b\b\b\b%3.0f%%',iLength/nPartlengths*100);

    for iBeat = 1:nUsedBeat % looping over all columns (arbitrarily chosen columns)

        help2 = find(potential_duplicates(1:nUsedBeat,iBeat));

        for kBeat = help2' % check only potential duplicates

            % measure how well two length len segments go together

            simArray(iBeat,kBeat,iLength) = ...

                quantile(diag(median_simmat(iBeat+(0:len-1),kBeat+(0:len-1))),param.seg.quantile);

        end

    end

    simArray(1:nUsedBeat,1:nUsedBeat,iLength) = simArray(1:nUsedBeat,1:nUsedBeat,iLength) ...

        + simArray(1:nUsedBeat,1:nUsedBeat,iLength)' ...

        - eye(nUsedBeat) .* simArray(1:nUsedBeat,1:nUsedBeat,iLength);

    simArray(:,:,iLength) = conv2(simArray(:,:,iLength),[.01 .98 .01], 'same');

    % take only over-average bars that do not overlap

    for iBeat = 1:nUsedBeat

        temp = simArray(:,iBeat,iLength) > param.seg.thresh_seg;

        decisionArray2(temp,iBeat,iLength) = ...

            simArray(temp,iBeat,iLength);

    end

    decisionArray2(:,:,iLength) = decisionArray2(:,:,iLength) .* (decisionArray2(:,:,iLength) >= maxfilt1(decisionArray2(:,:,iLength),len-1));

    decisionArray2(:,:,iLength) = decisionArray2(:,:,iLength) .* decisionArray2(:,:,iLength)';

%     potential_duplicates = triu(simArray(:,:,iLength) > param.seg.thresh_seg * 0.7);

    potential_duplicates = potential_duplicates .* (simArray(:,:,iLength) > param.seg.thresh_seg);

end

%% milk the data

bestvals = [];

for iLength = 1:nPartlengths

    currLogicSum = sum(decisionArray2(:,:,iLength)>0,2);

    for iBeat = 1:nBeat

        if currLogicSum(iBeat) > 1.0

            currSum = (mean(decisionArray2(decisionArray2(:,iBeat,iLength)>0,iBeat,iLength)))/2;

            bestvals = [bestvals; (currLogicSum(iBeat)-1) * partlengths(iLength), currSum, iLength, iBeat, currLogicSum(iBeat)];

        end

    end

end

%%

%%

% make a table of all valid sets of parts

partletters = {'A','B','C','D','E','F','G', 'H','I','J','K','L','M','N','O','P','Q','R','S'};

valid_sets = ones(1,size(bestvals,1));

parts = [];

if ~isempty(bestvals)

    bestvals(:,2) = bestvals(:,2) / max(bestvals(:,2)*2);

    bestvals(:,1) = bestvals(:,1) + bestvals(:,2);

    bestvals(:,2) = [];

    for kSeg = 1:7

        currbestvals = bestvals(logical(valid_sets),:);

        [m,i] = max(currbestvals(:,1));

        if isempty(m)

            break

        end

        bestLength = partlengths(currbestvals(i,2));

        bestIndices = decisionArray2(currbestvals(i,3),:,currbestvals(i,2));

        islands = conv2(1*(bestIndices > 0), [zeros(1,bestLength-1) ones(1,bestLength)], 'same');

        %%

        if param.seg.newcheck

            iBetter = i;

            bestSubval = 0;

            subm1 = find(currbestvals(:,4) > currbestvals(i,4)); % has to have more part instances

            for kSub = 1:length(subm1)

                iSub = subm1(kSub);

                subIndices = decisionArray2(currbestvals(iSub,3),:,currbestvals(iSub,2));

                subLength = partlengths(currbestvals(iSub,2));

                subIslands = conv2(1*(subIndices > 0), [zeros(1,subLength-1) ones(1,subLength)], 'same');

                if sum(islands.*subIslands) == currbestvals(i,4) * subLength ...

                        && currbestvals(iSub,1) > bestSubval

                    iBetter = iSub;

                    bestSubval = currbestvals(iSub,1);

                end

            end

            i = iBetter;

            bestLength = partlengths(currbestvals(i,2));

            bestIndices = decisionArray2(currbestvals(i,3),:,currbestvals(i,2));

            islands = conv2(1*(bestIndices > 0), [zeros(1,bestLength-1) ones(1,bestLength)], 'same');

        end

        %%

        part.n = bestLength;

        part.indices = find(bestIndices);

        part.letter = partletters{kSeg};

        part.level = kSeg;

        parts = [parts part];

        for iSet = find(valid_sets)

            currislands = conv2(2*(decisionArray2(bestvals(iSet,3),:,bestvals(iSet,2)) > 0), ...

                [zeros(1,partlengths(bestvals(iSet,2))-1) ones(1,partlengths(bestvals(iSet,2)))], 'same');

            if any(islands .* currislands > 0)

                valid_sets(iSet) = false;

            else

                %             stem(islands+currislands*.5)

                %             aaa = 1;

            end

        end

    end

else

    part.n = nBeat;

    part.indices = 1;

    part.letter = 'A';

    part.level = 1;

    parts = [parts part];

end

parts = [parts, nullpart(parts,1:nBeat)];

if param.seg.editor

    [pa, ta] = partarray(parts);

    parts = editorssearch(pa, ta, parts);

    parts = [parts, nullpart(parts,1:nBeat)];

end

parts = mergenulls(parts);

song.parts = parts;

% song.parts