function vc = vcMatrix(M,N,Gi)
% vcMatrix(M,N,Gi)
%	Since this function may take hours to generate a large Void and
%	Cluster matrix, this function checks Matlab's search path for
%	previously built matrices.  If a previously built matrix of the
%	appropriate size is found then it will load that.  If one cannot
%	be found then it will generate one.  This newly generated Void and
%	Cluster matrix will be saved for future use.
% 
%	Returns a MxN bluenoise matrix generated using the Void and 
%	Cluster method described by Ulichney.  The value of Gi specifies 
%	the seed value used.  The seed value Gi defaults to 1/9 if not 
%	specified.
%
% This function must execute three phases before completion.  Since this
% may take a considerable amount of time (a few hours) for large matrices
% the current phase is displayed as the function executes. 
%
% Functions called
%	gauss2()
%	cirConv2()
%	random()
%	scroll2()
%	scale()
%
% Rick Anthony
% Dan S. Lee

if(nargin < 2) N = M; end
if(nargin < 3) Gi = 1/9; end

str = ['vc' int2str(M) 'x' int2str(N) '.mat'];
if(exist(str))
   disp('Void and Cluster matrix found!  Loading');
   eval(['load ' str]);
   return;
end 

disp('Void and Cluster matrix not found.  Generating ...');
B = rand(M,N) < Gi*ones(M,N);
D = zeros(M,N);
g = gauss2(3.5,M,3.5,N);
[gx, gy] = find(g == max(max(g)));

b1 = B;

tempMax = cirConv2((b1==1),g);
tempMin = cirConv2((b1==0),g);

disp('vcMatrix:  Phase 1');
while 1,
    [maxYindex, maxXindex] = find(tempMax == max(max(tempMax.*(b1==1))));
    maxXindex = maxXindex(random(length(maxXindex)));
    maxYindex = maxYindex(random(length(maxYindex)));

    b1(maxYindex,maxXindex) = 0;
    tempMax = tempMax - scroll2(g, maxYindex-gx, maxXindex-gy);
    tempMin = tempMin + scroll2(g, maxYindex-gx, maxXindex-gy);

    [minYindex, minXindex] = find(tempMin == max(max(tempMin.*(b1==0))));
    minXindex = minXindex(random(length(minXindex)));
    minYindex = minYindex(random(length(minYindex)));

    b1(minYindex,minXindex) = 1;
    tempMin = tempMin - scroll2(g, minYindex-gx, minXindex-gy);
    tempMax = tempMax + scroll2(g, minYindex-gx, minXindex-gy);

    if((maxXindex == minXindex) & (maxYindex == minYindex))
        break;
    end
end


b2 = b1;
Rank = sum(sum(b2))-1;

disp('vcMatrix:  Phase 2');
while Rank >= 0,
    [maxYindex, maxXindex] = find(tempMax == max(max(tempMax.*(b2==1))));
    maxXindex = maxXindex(random(length(maxXindex)));
    maxYindex = maxYindex(random(length(maxYindex)));

    b2(maxYindex,maxXindex) = 0;
    tempMax = tempMax - scroll2(g, maxYindex-gx, maxXindex-gy);

    D(maxYindex,maxXindex) = Rank;

    Rank = Rank-1;
end


b3 = b1;
Rank = sum(sum(b3));

disp('vcMatrix:  Phase 3');
while Rank < M*N,
    [minYindex, minXindex] = find(tempMin == max(max(tempMin.*(b1==0))));
    minXindex = minXindex(random(length(minXindex)));
    minYindex = minYindex(random(length(minYindex)));

    b1(minYindex,minXindex) = 1;
    tempMin = tempMin - scroll2(g, minYindex-gx, minXindex-gy);

    D(minYindex,minXindex) = Rank;

    Rank = Rank+1;
end


vc = scale(D,0,1,-0.5,M*N-0.5);
disp('Saving Void and Cluster matrix for future use.');
eval(['save ' str ' vc']);