boundarySearch(X0,lb,ub,limitFunction, varargin)

function [xb,P, a] = boundarySearch(X0,lb,ub,limitFunction, varargin)
%BOUNDARYSEARCH searches a 2-D space to identify the failure region(s)
%  BOUNDARYSEARCH(X, LIMITFUNCTION) searches over a 2D space to identify
%  regions where the limit function returns (True) 1 if failure is present.
%  X0 (nx2) is initial search starting points of (x,y) pairs. LIMITFUNCTION
%  is a handle to the limit state function that that takes X as an input.
%
%  Optional input are:
%    'showplot'     show plot update of algorithm (voronoi cells and
%                   next points to update)
%    'ptol'         Perimeter change tolerance, followed by numeric value.
%
%  SYNTAX
%     [xb,P]   = boundarySearch(X0, limitFunction)
%     [xb,P,A] = boundarySearch(X0, limitFunction,'showplot','ptol',val)
%
%  Example
%    [x,y,z] = peaks(10);
%    X0 = [x(:) y(:)];
%    ub = [3 3];
%    lb = [-3 -3];
%    [xb P a] = boundarySearch(X0,lb,ub,@peaksLimitFunction,'showplot','ptol',0.1)
%
%  See also clusterdata, DelaunayTri, voronoi, peaksLimitFunction
%

%   Copyright 2009 The MathWorks, Inc.

%% Input Parsing
debug = false; showplot = false;
ptol = 1e-3; vin = 1; maxregion = 5;
if nargin > 2
    while vin <= length(varargin)
        switch lower(varargin{vin})
            case 'debug'
                debug = true;
            case 'showplot'
                showplot = true;
            case 'ptol'
                vin = vin+1;
                ptol = varargin{vin};
            case 'maxregion'
                vin = vin+1;
                maxregion = varargin{vin};
            otherwise
                warning('BOUNDARYSEARCH:InvalidInputs',...
                    'Ignoring additional input:  %s', varargin{vin});
        end
        vin = vin+1;
    end
end

%% Initialization
X = [];
failed = logical([]);
nextX = X0;
P = [];
Plast = 0;
deltaP = 100;
iter = 1;
while deltaP > ptol && iter < 1000
    X = [X; nextX];
    %% Call LimitFunction with New Inputs
    nextFail = limitFunction(nextX);
    failed = [failed; nextFail];
    nextX = [];
    
    %% Use Delaunay Triangulation
    % Later we'll use the convex hull and voronoi diagrams to update search
    % region
    DT = DelaunayTri(X(:,1),X(:,2));
    [V,R] = voronoiDiagram(DT);
    
    if debug || showplot
        voronoi(X(:,1),X(:,2)), drawnow
        axis([lb(1) ub(1) lb(2) ub(2)])
    end
    
    %% Test for failure point detected
    if sum(failed) > 0 % Found failure point(s)
        
        %% Use clustering to identify number of differnt regions
        zones = clusterdata([X,failed],maxregion);
        if debug || showplot
            %unique(zones)
            hold on
            scatter(X(:,1),X(:,2),50,zones,'filled')
            axis square equal
            axis([lb(1) ub(1) lb(2) ub(2)])
            hold off
            drawnow
        end
        
        %% Determine Number of Failure Regions
        zoneNumbers = unique(zones.*failed);
        zoneNumbers = zoneNumbers( zoneNumbers>0);
        
        %% Find Failure Boundary
        col = 'rgbmcy';
        for zn = 1:length(zoneNumbers)
            newDT = DelaunayTri(X(:,1),X(:,2));
            [V,R] = voronoiDiagram(newDT);
            % Remove other zones from newDT
            indx = zones == zoneNumbers(zn);
            if sum(indx) <= 2
                % if only one/two points, grab voronoi cell points and continue
                pts = find(indx > 0);
                for i = 1:length(pts)
                    nextX = [nextX; V(R{pts(i)},:)];
                end
                continue
            else
                indx = find(~indx); % convert to integer index
                newDT.X(indx,:) = [];
            end
            %% Find Each Failure Region's Convex Hull
            try
                [ch a(zn)] = convexHull(newDT);
            catch eDT
                % Check for empty triangulation (collinear points likely)
                if strcmpi(eDT.identifier,'MATLAB:DelaunayTri:TriIsEmptyErrId')
                    pts = find(indx > 0);
                    for i = 1:length(pts)
                        nextX = [nextX; V(R{pts(i)},:)];
                    end
                    continue
                else
                    rethrow(eDT)
                end
            end
            
            % Calculate the perimeter of the convex hull
            xb{zn} = newDT.X(ch,:);
            dX = diff(newDT.X(ch,:));
            P(zn)  = sum( sqrt( sum(dX.^2,2) ) );
            
            % Make note of any points that are in the convex hull or on it
            in = inpolygon(X(:,1),X(:,2),newDT.X(ch,1),newDT.X(ch,2));
            
            % Use the point in the convex hull to find the voronoi cells
            vc = R(in);
            
            % Use the voronoi cell points as the next points to evaluate
            idxNewPts = [];
            for j = 1:length(vc)
                idxNewPts = [idxNewPts, vc{j}];
            end
            idxNewPts = unique(idxNewPts);
            newX = V(idxNewPts,:);
            
            % Remove any points in/on the convex hull
            in = inpolygon(newX(:,1),newX(:,2),newDT.X(ch,1),newDT.X(ch,2));
            newX = newX(~in,:);
            
            % Update next points array
            nextX = [nextX;newX];
            
            if showplot || debug
                hold on
                plot(newDT.X(ch,1),newDT.X(ch,2),col(zn))
                plot(nextX(:,1),nextX(:,2),'r+')
                axis([lb(1) ub(1) lb(2) ub(2)])
                drawnow
                if debug, pause, end
                hold off
            end
            iter = iter+1;
        end %zn loop
    end
    %% Remove INF and repeat (x,y) pairs
    if isempty(nextX)      % No failure points found, use voronoi diagram points to update
        nextX = V;
        deltaP = deltaP - 1;
    else
        if isempty(P)
            deltaP = deltaP -1;
        else
            deltaP = sum(P) - Plast;
            Plast = sum(P);
        end
    end
    nextX = unique(nextX,'rows');
    remove = (nextX(:,1) == Inf | nextX(:,2) == Inf);
    nextX(remove,:) = [];
    
    %% Remove any Points Not In Bounds
    % Remove x's
    nextX = nextX(nextX(:,1) <= ub(1),:);
    nextX = nextX(nextX(:,1) >= lb(1),:);
    nextX = nextX(nextX(:,2) <= ub(2),:);
    nextX = nextX(nextX(:,2) >= lb(2),:);
    
end %while
%repeat