% Capacity-Constrained Voronoi Tesselation - 04/29/11
% Jean-Marc Berthommé
%
% - 04/29/11 :
%   . implementation of "Capacity-Constrained Voronoï Diagrams in Finite
%     Spaces" by M. Balzer and D. Heck 2008

function ccvt
n = 3;                      % nb of sites
r = 30;                     % resolution   -  10 -   200
m = r^2;                    % nb of points - 100 - 40000
r = sqrt(m);                % resolution
k = 0;                      % iteration counter
stable = false;             % loop flag
scr = get(0,'ScreenSize');  % screen size

RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(100*clock)));

% sites coordinates
sx = floor(r*[0.40, 0.70, 0.30]); % n = 3
sy = floor(r*[0.25, 0.50, 0.80]);
% sx = 1+floor(r*rand(1,n));
% sy = 1+floor(r*rand(1,n));

% m points equidistributed on n sites
S = 1:n;
I = 1+floor(n*rand(r,r));
provide_capacities(n, m, I, 'CCVT');

% display I
f1 = figure(1); set(gcf,'Color',[0.2,0.2,0.2]);
set(f1,'Position',[2*scr(3)/3, scr(4)/2, scr(3)/3, scr(4)/2-45*2]);
image(uint8(I)); colormap(gray(1+n)); axis image;
hold on; plot(sx, sy, 'ko'); hold off;

% distance map
[X Y]   = meshgrid(1:r,1:r); % image grid
[SX XX] = ndgrid(sx, X);
[SY YY] = ndgrid(sy, Y);
D2 = (SX-XX).^2+(SY-YY).^2;

[~,MIN] = min(D2);
MIN = reshape(MIN, size(I));
provide_capacities(n, m, MIN, 'Voronoï');

f2 = figure(2); set(gcf,'Color',[0.2,0.2,0.2]);
set(f2,'Position',[2*scr(3)/3, 0, scr(3)/3, scr(4)/2-50*2]);
image(uint8(MIN)); colormap(gray(1+n)); axis image;
hold on; plot(sx, sy, 'ko'); hold off;
title('Voronoï tesselation','Color','w');

p = init_pair(S);

while ~stable
    stable = true;
    k = k+1; tit = sprintf('CCVT - k = %d', k);

    for i=1:nchoosek(n,2); % nb of remaining increments
        H1 = zeros(size(I(I == p(1))));
        H2 = zeros(size(I(I == p(2))));

        X1 = X(I == p(1)); Y1 = Y(I == p(1));
        X2 = X(I == p(2)); Y2 = Y(I == p(2));
        for j=1:length(H1)
            H1(j) = delta(X1(j), Y1(j), sx, sy, p(1), p(2));
        end
        for j=1:length(H2)
            H2(j) = delta(X2(j), Y2(j), sx, sy, p(2), p(1));
        end

        while ~isempty(H1) && ~isempty(H2) && max(H1)+max(H2) > 0
            [~, Id1] = max(H1);

            I(r*(X1(Id1)-1) + Y1(Id1)) = p(2);
            H1(Id1) = []; X1(Id1) = []; Y1(Id1) = [];

            [~, Id2] = max(H2);

            I(r*(X2(Id2)-1) + Y2(Id2)) = p(1);
            H2(Id2) = []; X2(Id2) = []; Y2(Id2) = [];

            figure(1);
            image(uint8(I)); colormap(gray(1+n)); axis image;
            hold on; plot(sx, sy, 'ko'); hold off; title(tit,'Color','w');

            stable = false;
        end

        p = inc_pair(S, p);
    end
end

disp('Press any key to continue...');
pause; clear all; close all;

%*************************************************************************%
function delta = delta(x, y, sx, sy, p1, p2)
delta = ((x-sx(p1)).^2+(y-sy(p1)).^2) - ((x-sx(p2)).^2+(y-sy(p2)).^2);

% function E = energy_map(n,m,I,D2)
% Id = ones(n,1)*(I(:)');
% [~, Yd] = meshgrid(1:m,1:n);
% Idx = (Yd == Id);
% E = D2(Idx);             % energy vector
% E = reshape(E, size(I)); % energy matrix

function p = init_pair(S)
% Initialize pair
p = zeros(1,2);
p(1) = S(1); p(2) = S(2);

function p = inc_pair(S, p)
% Increment pair
if p(2) + 1 <= max(S)
    p(2) = p(2) + 1;
elseif p(1) + 1 <= max(S) - 1
    p(1) = p(1) + 1;
    p(2) = p(1) + 1;
else
    p = init_pair(S);
end

function provide_capacities(n, m, I, str)
C = zeros(1,n); % capacities vector
str = sprintf('%s%s',str,'\n');
fprintf(str);
for i=1:n
    C(i) = sum(sum(I==i))/m;
    fprintf('C[%d] = %0.1f %%\n', i, 100*C(i));
end
N = 1/n*ones(1,n);
fprintf('-> # ~ %0.2f %%\n\n', 100*mean(abs(C-N)./N));

if abs(sum(C) - 1) > eps % check the sum
    error('Sum of capacities not equal to 1');
end

CCVT
C[1] = 36.2 %
C[2] = 32.2 %
C[3] = 31.6 %
-> # ~ 5.78 %

Voronoï
C[1] = 30.1 %
C[2] = 40.4 %
C[3] = 29.4 %
-> # ~ 14.22 %

Press any key to continue...