function varargout=VoronoiBound(varargin)
%% VoronoiBound   Calculates a Voronoi diagram with inner and outer bounds.
%   [VX,VY] = VoronoiBound(PX,PY,OX,OY) produces a 2D Voronoi diagram for
%   the points in vectors PX and PY. The Voronoi diagram will be bounded by
%   the polygon described by the vertices in vectors OX and OY.
%   For the i-th point, the cooresponding Voronoi region is the polygon
%   formed by the vectors in the i-th cell entry of VX and VY. Points
%   outside of the bounds return an empty Voronoi region.
%
%   [...] = VoronoiBound(PX,PY,OX,OY,IX,IY) adds an inner boundary to
%   the calculation of the Voronoi diagram. The inner boundary is a polygon
%   described by vertices in vectors IX and IY.
%
%   [VX,VY,L] = VoronoiBound(...) adds the logical vector L to
%   flag the points that fall within the bounding polygon(s).
%
%   VoronoiBound(...) returns a figure as the only output.
%
%   [...] = VoronoiBound() uses a set of example data.

% Last Updated: May 25, 2018; 12:10 PM EDT
% Brian Phillips; University of Maryland; bphilli@umd.edu

%% Inputs

if nargin==0
    P = gallery('uniformdata',[30 2],0);
%     P=[0.25 0.25; 0.25 0.75; 0.75 0.25; 0.75 0.75;0.5,1];
%     P=[0.25 0.25; 0.5 0.25; 0.75 0.25; .9 0.25; .1 0.25];
    
    x_bound=[0;0;0.5;1;1;nan;0.65;0.65;0.35;0.35];
    y_bound=[0;1;1.5;1;0;nan;0.35;0.85;0.85;0.35]; 
end

if nargin>=2
    % Matrix of (x,y) tap locations
    Px=varargin{1}(:);
    Py=varargin{2}(:);
    P=[Px,Py];
    
    x_bound=[];
    y_bound=[];
end

if nargin>=4
    x_outer=varargin{3}(:);
    y_outer=varargin{4}(:);

    % Clockwise outer boundary (vertices)
    [x_outer_cw,y_outer_cw]=poly2cw(x_outer,y_outer);
    
    x_bound=x_outer_cw;
    y_bound=y_outer_cw;
end

if nargin>=6
    x_inner=varargin{5}(:);
    y_inner=varargin{6}(:);
    
    % Counter-clockwise interior boundary (vertices)
    [x_inner_ccw,y_inner_ccw]=poly2ccw(x_inner,y_inner);
    
    x_bound=[x_outer_cw;nan;x_inner_ccw];
    y_bound=[y_outer_cw;nan;y_inner_ccw];
end
    
%% Main

% Extract taps that fall within boundaries
[in,on]=inpolygon(P(:,1),P(:,2),x_bound,y_bound);
P_inl=in==1|on==1;
P_ind=find(P_inl);
P_in=P(P_ind,:);

% Built-in MATLAB function
warning('off','MATLAB:triangulation:EmptyTri2DWarnId');
DT=delaunayTriangulation(P_in);

% Ensure that Delaunay triangulation was successful
if isempty(DT.ConnectivityList)~=1

    % Built-in MATLAB functions
    [V,R]=voronoiDiagram(DT);
    [vx,vy]=voronoi(DT);
    vx=vx.';
    vy=vy.';

    % Eliminate zero-length vectors
    vd=sqrt((vx(:,2)-vx(:,1)).^2+(vy(:,2)-vy(:,1)).^2);
    vx(vd==0,:)=[];
    vy(vd==0,:)=[];

    % Find unique values in vx (and vy) to replace "infinite vertices"
    [u,~,indu]=unique([vx(:),vy(:)],'rows');
    uu=u(accumarray(indu,1)==1,:);
    ux=uu(:,1);
    uy=uu(:,2);

    % Create infinite vectors to help bound infinite regions
    kx=zeros(length(ux),2);
    ky=zeros(length(uy),2);

    for ii=1:length(ux)
        [induu]=find(vx==ux(ii)&vy==uy(ii));
        [indi,~]=ind2sub(size(vx),induu);
        kx(ii,:)=[vx(indi,1);vx(indi,2)].';
        ky(ii,:)=[vy(indi,1);vy(indi,2)].';
    end

    % Extend the infinite line segments
    kd=sqrt((kx(:,2)-kx(:,1)).^2+(ky(:,2)-ky(:,1)).^2);
    kextend=max(abs([x_bound;y_bound]))*3;
    kx(:,2)=kx(:,1)+(kx(:,2)-kx(:,1))./kd*kextend;
    ky(:,2)=ky(:,1)+(ky(:,2)-ky(:,1))./kd*kextend;

    % Sort the infinite line segments in counter-clockwise order
    xmid=(max(P_in(:,1))-min(P_in(:,1)))/2;
    ymid=(max(P_in(:,2))-min(P_in(:,2)))/2;
    [~,kbound]=sort(atan2(ky(:,2)-ymid,kx(:,2)-xmid));
    kbound(end+1)=kbound(1);

    % Initialize variables
    x_trib_in=cell(length(R),1);
    y_trib_in=cell(length(R),1);
    R_new=cell(length(R),1);
    V_new=V;

    % Loop over all tap tributary areas
    for ii=1:size(R,1)
        % Eliminate repeated vertices
        [~,IA,~]=unique([V(R{ii},1),V(R{ii},2)],'rows','stable');
        R{ii}=R{ii}(IA);

        % Fix regions with infinite vertices
        if any(R{ii}==1)
            % Look for infinite vectors that begin closest to each end of unbounded region
            v1_ind=knnsearch([kx(:,1),ky(:,1)],[V(R{ii}(2),1),V(R{ii}(2),2)],'Distance','euclidean','k',2);
            v2_ind=knnsearch([kx(:,1),ky(:,1)],[V(R{ii}(end),1),V(R{ii}(end),2)],'Distance','euclidean','k',2);

            flag1=kx(v1_ind(1),1)==kx(v1_ind(2),1)&&ky(v1_ind(1),1)==ky(v1_ind(2),1);
            flag2=kx(v2_ind(1),1)==kx(v2_ind(2),1)&&ky(v2_ind(1),1)==ky(v2_ind(2),1);

            % Method 1: use if infinite vectors begin in same location
            if flag1==1||flag2==1
                for jj=1:length(kbound)-1
                    % Trial vertices
                    v1=kbound(jj);
                    v2=kbound(jj+1);
                    Vx=[kx(v2,2);kx(v2,1);V(R{ii}(3:end-1),1);kx(v1,1);kx(v1,2)];
                    Vy=[ky(v2,2);ky(v2,1);V(R{ii}(3:end-1),2);ky(v1,1);ky(v1,2)];

                    % Update V if trial vertices enclose tap
                    if inpolygon(P_in(ii,1),P_in(ii,2),Vx,Vy)==1
                        lv=size(V,1);
                        V_new(lv+1,:)=[Vx(1),Vy(1)];
                        V_new(lv+2,:)=[Vx(end),Vy(end)];
                        R_new{ii}=[lv+1,R{ii}(2:end),lv+2];
                    end

                end

            % Method 2: (quicker) use if infinite vectors begin in unique locations
            else
                % Update V using infinite vectors
                lv=size(V,1);
                V_new(lv+1,:)=[kx(v1_ind(1),2),ky(v1_ind(1),2)];
                V_new(lv+2,:)=[kx(v2_ind(1),2),ky(v2_ind(1),2)];
                R_new{ii}=[lv+1,R{ii}(2:end),lv+2];
            end

        else
            R_new{ii}=R{ii};
        end

        % Make voronoi regions clockwise
        [xcw,ycw]=poly2cw(V_new(R_new{ii},1),V_new(R_new{ii},2));

        % Calculate the portion of tributary area bound within limits
        [x_trib_temp,y_trib_temp]=polybool('intersection',xcw,ycw,x_bound,y_bound);
        x_trib_in{ii}=x_trib_temp;
        y_trib_in{ii}=y_trib_temp;
    end

% If Delaunay triangulation was unsuccessful (e.g., due to collinear points), add four bounding points
else
  
    % Append list of points with four new bounding points
    Pextend=max(abs([x_bound;y_bound]))*10;
    P_add=[P_in;[1 1; -1 1; -1 -1; 1 -1]*Pextend];
    
    % Built-in MATLAB functions
    DT=delaunayTriangulation(P_add);
    [V,R]=voronoiDiagram(DT);
    
    % Remove Voronoi regions created by additional points
    R_new=R(1:size(P_in,1));
    V_new=V;
    
    % Initialize variables
    x_trib_in=cell(length(R_new),1);
    y_trib_in=cell(length(R_new),1);

    % Loop over all tap tributary areas
    for ii=1:size(R_new,1)
        % Make Voronoi regions clockwise
        [xcw,ycw]=poly2cw(V_new(R_new{ii},1),V_new(R_new{ii},2));

        % Calculate the portion of tributary area bound within limits
        [x_trib_temp,y_trib_temp]=polybool('intersection',xcw,ycw,x_bound,y_bound);
        x_trib_in{ii}=x_trib_temp;
        y_trib_in{ii}=y_trib_temp;
    end
    
end

% Initialize variables
x_trib=cell(size(P,1),1);
y_trib=cell(size(P,1),1);

% Expand to all taps (including those outside bounds)
for ii=1:length(P_ind)
    x_trib{P_ind(ii)}=x_trib_in{ii};
    y_trib{P_ind(ii)}=y_trib_in{ii};
end

%% Outputs

% Define outputs of function
if nargout==3
    varargout={x_trib,y_trib,P_inl};
elseif nargout==2
    varargout={x_trib,y_trib};
else
    varargout={};
end

%% Plot

if nargout==0
    
    % Outer limits for plot based on outer boundaries
    xmax=max(x_bound);
    xmin=min(x_bound);
    ymax=max(y_bound);
    ymin=min(y_bound);
    
    % Figure 1: Taps and bounds
    
    figure
    hold on
    % Plot taps
    plot(P(P_inl,1),P(P_inl,2),'ok')
    plot(P(~P_inl,1),P(~P_inl,2),'xk')
    % Plot surface bounds
    if any(isnan(x_bound))
        ind=find(isnan(x_bound));
        plot([x_bound(1:ind-1);x_bound(1)],[y_bound(1:ind-1);y_bound(1)],'-k','linewidth',2)
        plot([x_bound(ind+1:end);x_bound(ind+1)],[y_bound(ind+1:end);y_bound(ind+1)],'-k','linewidth',2)
    else
        plot([x_bound;x_bound(1)],[y_bound;y_bound(1)],'-k','linewidth',2)
    end      
    hold off
    axis([xmin,xmax,ymin,ymax])
    title('Taps and Bounds')
    
    % If Delaunay Triangulation is not empty
    if isempty(DT.ConnectivityList)~=1
      
        % Figure 2: Delaunay Triangulation
        
        figure
        hold on
        for ii=1:size(DT.ConnectivityList,1)
        patch(DT.Points(DT.ConnectivityList(ii,:),1),DT.Points(DT.ConnectivityList(ii,:),2),[0 0.1 1],'EdgeColor',[1 0 0],'linewidth',2,'FaceAlpha',0.05)
        end
        % Plot taps
        plot(P(P_inl,1),P(P_inl,2),'ok')
        plot(P(~P_inl,1),P(~P_inl,2),'xk')
        % Plot surface bounds
        if any(isnan(x_bound))
            ind=find(isnan(x_bound));
            plot([x_bound(1:ind-1);x_bound(1)],[y_bound(1:ind-1);y_bound(1)],'-k','linewidth',2)
            plot([x_bound(ind+1:end);x_bound(ind+1)],[y_bound(ind+1:end);y_bound(ind+1)],'-k','linewidth',2)
        else
            plot([x_bound;x_bound(1)],[y_bound;y_bound(1)],'-k','linewidth',2)
        end      
        hold off
        axis([xmin,xmax,ymin,ymax])
        title('Delaunay Triangulation')
        
        % Figure 3: Perpendicular Bisectors
        
        figure
        hold on
        for ii=1:size(DT.ConnectivityList,1)
        patch(DT.Points(DT.ConnectivityList(ii,:),1),DT.Points(DT.ConnectivityList(ii,:),2),[0 0.1 1],'EdgeColor',[1 0 0],'linewidth',2,'FaceAlpha',0.05)
        end
        % Plot taps
        plot(P(P_inl,1),P(P_inl,2),'ok')
        plot(P(~P_inl,1),P(~P_inl,2),'xk')
        % Plot surface bounds
        if any(isnan(x_bound))
            ind=find(isnan(x_bound));
            plot([x_bound(1:ind-1);x_bound(1)],[y_bound(1:ind-1);y_bound(1)],'-k','linewidth',2)
            plot([x_bound(ind+1:end);x_bound(ind+1)],[y_bound(ind+1:end);y_bound(ind+1)],'-k','linewidth',2)
        else
            plot([x_bound;x_bound(1)],[y_bound;y_bound(1)],'-k','linewidth',2)
        end      
        for ii=1:size(vx,1)
            plot(vx(ii,:),vy(ii,:),'-b','linewidth',2)
        end
        for ii=1:size(kx,1)
            plot(kx(ii,:),ky(ii,:),'-b','linewidth',2)
        end
        hold off
        axis([xmin,xmax,ymin,ymax])
        title('Perpendicular Bisectors')
        
    end
    
    % Figure 4: Voronoi Regions
    
    N_taps=size(P,1);
    c_map=jet(N_taps);
    figure
    for ii=1:N_taps
        % Plot tributary areas with holes
        if any(isnan(x_trib{ii}))&&P_inl(ii)
            [patch_f,patch_v]=poly2fv(x_trib{ii},y_trib{ii});
            patch('Faces',patch_f,'Vertices',patch_v,'FaceColor',c_map(ii,:),'LineStyle','none','FaceAlpha',0.5);

            ind=find(isnan(x_trib{ii}));
            patch(x_trib{ii}(1:ind-1),y_trib{ii}(1:ind-1),c_map(ii,:),'FaceColor','none','linewidth',2);
            patch(x_trib{ii}(ind+1:end),y_trib{ii}(ind+1:end),c_map(ii,:),'FaceColor','none','linewidth',2);
        % Plot tributary areas without holes
        elseif P_inl(ii)
            patch(x_trib{ii},y_trib{ii},c_map(ii,:),'linewidth',2,'FaceAlpha',0.5);   
        end    
    end
    hold on
    % Plot taps
    plot(P(P_inl,1),P(P_inl,2),'ok')
    plot(P(~P_inl,1),P(~P_inl,2),'xk')
    % Plot surface bounds
    if any(isnan(x_bound))
        ind=find(isnan(x_bound));
        plot([x_bound(1:ind-1);x_bound(1)],[y_bound(1:ind-1);y_bound(1)],'-k','linewidth',2)
        plot([x_bound(ind+1:end);x_bound(ind+1)],[y_bound(ind+1:end);y_bound(ind+1)],'-k','linewidth',2)
    else
        plot([x_bound;x_bound(1)],[y_bound;y_bound(1)],'-k','linewidth',2)
    end
    hold off
    axis([xmin,xmax,ymin,ymax])
    title('Voronoi Regions')

end

warning('on','MATLAB:triangulation:EmptyTri2DWarnId');

end

