function [clusts, centers, numClust, affinity, dist] = ...
  doSimilarityComponentClustering(affinity, dist, ...
  numClusters, features, positivePatches, patchIndexes, relLocs)
% Generate affinity matrix.
if isempty(affinity)
  disp('Generating affinity matrix');
  tic;
  [affinity, dist] = getL2DistanceAffinityLoc(features, patchIndexes, ...
    relLocs);
  % [affinity, dist] = getL2DistanceAffinity(features);
  % [affinity, dist] = getHistogramIntersectionAffinity(features);
  toc;
end
disp('Doing Segmentation');
tic;
clusts = doNCutSegmentation(affinity, numClusters);

% [clusts, affinity] = doComponentSegmentation2(affinity, positivePatches);
toc;
disp('Done');

numClust = length(unique(clusts));
centers = getClusterCenters(features, clusts, numClust);
end

function centers = getClusterCenters(features, clusts, numClusts)
centers = zeros(numClusts, size(features, 2));
for i = 1 : numClusts
  inds = find(clusts==i);
  feats = features(inds, :);
  centers(i, :) = sum(feats) ./ length(inds);
end
end

function clusts = doNCutSegmentation(affinity, numClusters)
disp('N-Cutting');
members = ncutW(affinity, numClusters);
clusts = sum(members .* repmat(1:numClusters, size(members, 1), 1), 2);
uclust = unique(clusts);
index(uclust) = 1:length(uclust);
clusts = index(clusts);
end

function members = doComponentSegmentation(affinity, numClusters)
disp('Componenting');
maxIter = 10;
thresh = 0.1;
acceptableDeviation = numClusters * 0.1;
for i = 1 : maxIter
  aff = affinity;
  aff(aff<thresh) = 0;
  fprintf('.');
  [numComps, members] = graphconncomp(sparse(aff), 'Directed', false);
  fprintf('[%d, %f]\n', numComps, thresh);
  if abs(numClusters - numComps) < acceptableDeviation
    break;
  end
  thresh = thresh + thresh * ...
    min(abs(numClusters - numComps) / numClusters, 0.2) * ...
    sign(numClusters - numComps);
end
end

function [members, affNew] = doComponentSegmentation2(affinity, ...
  positivePatches)
affut = triu(affinity);
allAff = affut(:);
[~, inds] = sort(allAff, 'descend');

degreeTally = zeros(size(affinity, 1), 1);
wbar = waitbar(0, '');
[R C] = ind2sub(size(affinity), inds);
cutoff = 3000;
maxDegree = 2;
affNew = zeros(size(affinity));
includedEdges = 0;
for i = 1 : length(inds)
  if rem(includedEdges, 10) == 0
    waitbar(includedEdges/cutoff, wbar);
  end
  if allAff(inds(i)) < eps
    break;
  end
  r = R(i);
  c = C(i);
  D1 = positivePatches(r);
  D2 = positivePatches(c);
  if degreeTally(r) >= maxDegree || ...
      degreeTally(c) >= maxDegree || ...
      r == c || ...
      strcmp(D1.im, D2.im)
    % Do nothing.
  else
    includedEdges = includedEdges + 1;
    degreeTally(r) = degreeTally(r) + 1;
    degreeTally(c) = degreeTally(c) + 1;
    affNew(c, r) = affut(r, c);
    affNew(r, c) = affut(r, c);
  end
  if includedEdges == cutoff
    break;
  end
end
close(wbar);
disp('Finding connected components');
[~, members] = graphconncomp(sparse(affNew), 'Directed', false);
disp('Done');
end

function affinity = getAffinityMatrix1(features)
% Affinity by the weird distance.
means = sum(features, 2) / size(features, 2);
normFeat = features - repmat(means, 1, size(features, 2));
affinity = normFeat * normFeat';
% affinity = affinity .* affinity';
maxVal = max(max(affinity));
minVal = min(min(affinity));
affinity = (affinity - minVal) / (maxVal - minVal);
end

function affinity = getAffinityMatrix3(features)
% Affinity weighted by average value.
wts = exp(-features);
fwts = features.*wts;
ata = zeros(size(features, 1));
parfor i = 1 : size(features, 1)
  ata(i, :) = fwts(i, :) * ...
    (repmat(features(i, :), size(features, 1), 1) .* wts)';
end
btb = ata';
atb = fwts * fwts';
bta = atb';
d = ata + btb - bta - atb;
affinity = exp(-d);
end
