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