% Plot commands for the 3D voronoi polyhedra.
% Arguments : M, indexed individuals
% S, a matrix of the referenced individuals to be highlighted.
% for example, S = [1 10 19] would highlight the first, tenth, and
% nineteenth individual (ranked in order of ascending distance from center)

function vorPlot3dV03(M,S)
% First, there's no point in running this function if the selected cells
% don't exist.  If the user want us to highlight cell #75, and our
% individual matrix (M) only has 20 points, this whole thing is going to
% break, and the user won't know why.  So we check to make sure that that's
% not going to happen, and we send a little polite message to the user.  The
% figure screen will still pop up, but it will be blank.
cantwork = 0;
for i = 1:size(S,2)
    if (S(1,i)>size(M,1)&(cantwork==0))
        'Please limit your cell selections (S) to cells within the array.'
        cantwork = 1;
    end
end

% So, if we can continue, we do so!
if (cantwork == 0)
    
    % find the center of mass, using the center_of_mass function
    
    center = [];
    center = center_of_mass(M);
    centerX = center(1,1);
    centerY = center(1,2);
    centerZ = center(1,3);
    
    % sort the array, with those closest to the center at the top of the array
    
    M = sortbydistance(M);
    
    % find the greatest distance between any two points, which is going to
    % be our group diameter
    
    currentMax = 0;
    j=0;
    
    for i=1:(size(M,1))
        for j=1:(size(M,1))
            if dist3d(M(i,2), M(i,3), M(i,4), M(j,2), M(j,3), M(j,4))>currentMax
                currentMax = dist3d(M(i,2), M(i,3), M(i,4), M(j,2), M(j,3), M(j,4));
            end
        end
    end
    
    diameter = currentMax;
    
    
    % define a regular dodecahedron, centered at the origin with radius of
    % one.  This was produced in Maple, and just transplanted here.
    
    BOUND = [[.3568220896, -.9341723582, 0.]; [-.3568220896, -.9341723582, 0.]; [-.5773502686, -.5773502686, .5773502686]; [0., -.3568220896, .9341723586]; [.5773502686, -.5773502686, .5773502686]; [0., .3568220896, -.9341723582]; [-.5773502686, .5773502686, -.5773502686]; [-.3568220896, .9341723586, 0.]; [.3568220896, .9341723586, 0.]; [.5773502686, .5773502686, -.5773502686]; [0., -.3568220896, -.9341723582]; [-.5773502686, -.5773502686, -.5773502686]; [-.9341723582, 0., -.3568220896]; [-.9341723582, 0., .3568220896]; [-.5773502686, .5773502686, .5773502686]; [0., .3568220896, .9341723586]; [.5773502686, .5773502686, .5773502686]; [.9341723586, 0., .3568220896]; [.9341723586, 0., -.3568220896]; [.5773502686, -.5773502686, -.5773502686]];
    
    % explode BOUND about the origin.  This makes the bounding dodecahedron
    % just a bit bigger than the diameter of the group, so that all
    % points of the dodecahedron's surface will be outside of the group
    
    BOUND = BOUND*diameter*1.5; 
    
    % this removes our index from M, so that it and BOUND will have the
    % same number of columns (three, in fact, with the form [x y z]).
    M = remindex(M);
    
    % Now we shift the bounding dodecahedron to be centered around our
    % group's center of mass, which we found earlier.  After this step, the
    % group will be contained by the dodecahedron, as planned.
    for i = 1:size(BOUND,1)
        BOUND(i,1)=BOUND(i,1)+centerX;
        BOUND(i,2)=BOUND(i,2)+centerY;
        BOUND(i,3)=BOUND(i,3)+centerZ;
    end
    
    % We build a full matrix (F), which is our original individual points
    % and the bounding vertices of the dodecahedron.  It's important that
    % the original matrix (M) is at the top, as we shall see.
    F = [M; BOUND];
    
    % Now we actually do the voronoi calculation.  V holds the vertices of
    % the voronoi cells, C is the index of the cells themselves.  Voronoin
    % is a default Matlab function, so see "help voronoin" for a good
    % explanation of what's going on here.
    [V,C] = voronoin(F);
    
    % Here we start the graphics calls.
    hold on
    figure(1)
    
    % We perform a convex hull calculation on the cells from our original
    % group (which is why we set hull equal to the size of M, instead of
    % the size of F).  You can check "help convhulln" in Matlab to find out
    % about the convex hull function.
    
    for hull = 1:size(M,1);
        
        % For each volume 'hull'
        X = V(C{hull},:);      % View hull Voronoi cell.
        K = convhulln(X);
        d = [1:size(K,2), 1];
        % d = [1 2 3 1]; Index into K
        % For each face i;
        
        for i = 1:size(K,1)
            j = K(i,d);
            % Plot a polygon that is almost completely transparent green,
            % on each face i, with slightly less opaque edges
            h(i)=patch(X(j,1),X(j,2),X(j,3),'w');
            set(h(i),'facecolor', 'green');
            set(h(i),'facealpha', .15);
            set(h(i),'edgealpha', .05);
            
        end
        for q=1:size(S,2);
            
            if (hull==S(1,q))
                
                for i = 1:size(K,1)
                    j = K(i,d);
                    % Plot a polygon that is white with black lines, on each face i
                    h(i)=patch(X(j,1),X(j,2),X(j,3),'w');
                    
                    % set the faces of this to blue
                    set(h(i),'facecolor', 'blue');
                    set(h(i),'edgealpha', .05);
                    
                    
                    
                end
            end
        end
    end
end
hold off
view(3)
axis equal
title('highlighted voronoi cells')
rotate3d;