www.gusucode.com > mbcdesign 工具箱 matlab 源码程序 > mbcdesign/@conconvexhull/findBoundaryPoints.m

    function [c, bp, ok] = findBoundaryPoints(c, opts, X)
%FINDBOUNDARYPOINTS -- Find the boundary points for a data set and a constraint
%
%  [CON, BP] = FINDBOUNDARYPOINTS(CON, OPTS, X)
%
%  Note that for a convex hull constraint, finding the boundary points is
%  equivalent to find the constraint.
%
%  See also CONCONVEXHULL.

%  Copyright 2005-2012 The MathWorks, Inc.
Xc = pFilterFactors( c, X );

% normalise Xc to be in [-1, 1]
[scale, shift, XcNorm] = scaleShift(Xc);

% Compute the convex hull
try
    k = convhulln( XcNorm );
catch
    % CONVHULLN can error is if there is a problem with the data, e.g., if
    % all the data lies on a line. This corresponds to a "failed to fit"
    k = [];
end

% Data size
nFaces = size( k, 1 );
nDim = size( XcNorm, 2 ); 

% Allocate space for the matrix A and right-hand side b
A = zeros( nFaces, nDim );
b = zeros( nFaces, 1 );

% Compute the contents of A and b
xmean = mean( XcNorm ); 
j1 = 2:nDim;
j2 = ones( nDim - 1, 1 );
for i = 1:nFaces,
    % Xi is a matrix with one row for each point on face i
    Xi = XcNorm(k(i,:),:);
    % 
    Ai = Xi(j1,:) - Xi(j2,:);
    % The null space of Ai is orthogonal to the row space of Ai
    ni = null( Ai );
    if size( ni, 2 ) ~= 1,
        % The null space spans more than one dimension
        ni = Inf( nDim, 1 );
    end
    
    b(i) = Xi(1,:) * ni;

    if xmean * ni > b(i),
        % The middle point is more outside than the edge point. This means
        % that the normal "ni" points in the wrong direction.
        ni   = -ni;
        b(i) = -b(i);
    end
    
    A(i,:) = ni.';
end

% simplify conv hull
if ~isnull(opts) && ~get(opts, 'KeepAllFacets')
    tol = get(opts, 'Tolerance');
    
    % find unique rows within Inf-norm tolerance
    uniqueX = curvefitlib.internal.uniqueWithinTol([A,b], tol);

    % set A,b values 
    A = uniqueX(:,1:end-1);
    b = uniqueX(:,end);
end

if ~isempty(A)
    % un-shift and scale
    % refer to scaleShift function
    b = b+A*shift';
    A = bsxfun(@times, A, scale);
end

i = isfinite( b );

% Set the output
ok = ~isempty( k );
c.A = A(i,:);
c.b = b(i);
c = pUpdateScaleFactor( c );
bp = unique( k(:) );
if ok
    c.CenterPoint = mean( Xc(bp,:) );
    % make sure all X are inside the boundary
    md = max(constraintDistance(c,X));
    if md>0
        c.b = c.b + 2*md*c.ScaleFactor;
    end
else
    c.CenterPoint = zeros( 1, size( Xc, 2 ) );
end

function [scale,shift,XcNorm] = scaleShift(Xc)
if isempty(Xc)
    XcNorm = Xc;
    scale = eye(size(Xc, 1));
    shift = zeros(1,size(Xc,1));
else
    % XcNorm = Xc*scale - shift
    a = min(Xc, [], 1);
    b = max(Xc, [], 1);
    scale = 2./(b-a);
    shift = (2*a)./(b-a)+1;
    XcNorm = bsxfun(@minus, bsxfun(@times, Xc, scale), shift);
end