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

    function [c,L,I] = xregfindcenter( x, method, options )
%XREGFINDCENTER find the center of a set of points
%
%   XREGFINDCENTER(X) finds the center of the points whose cartesian
%   coordinates are given as the rows of X. 
%
%   XREGFINDCENTER(X,'LeastSquares') finds the center in the least squares
%   sense. This is the default method. 
%
%   XREGFINDCENTER(X,'MinEllipse') finds the center of the minimal
%   enclosing ellipsoid. 
%
%   [C,L,I] = FIND_CENTER(X,'MinEllipse') also returns the triangular matrix
%   L that defines the ellipsoid and the list I of indices that are on the
%   boundary of the ellipse.
%
%   FIND_CENTER(X, 'MinEllipse', OPTIONS) passes the given OPTIONS to
%   FMINCON.
%
%   XREGFINDCENTER('List') returns a cell array of possible options for the
%   center finding algorithms. 

%  Copyright 2000-2009 The MathWorks, Inc. and Ford Global Technologies, Inc.

if nargin < 2, 
    method = 'LeastSquares';
end
if ischar( x ),
    c = { 'LeastSquares', 'MinEllipse' };
    return
end

if nargin < 3,
    % set up options struct
    options = optimset( fmincon('defaults'),...
        'largescale','off',...
        'Algorithm','active-set',...
        'tolx',   1e-8, ...
        'tolfun', 1e-8, ...
        'tolcon', 1e-8, ...
        'MaxFunEvals', 1000, ...
        'Display', 'off'  );  % 'iter', 'off'
end
options.LargeScale = 'off';
options.MaxSQPIter = 20 * size( x, 2 )^2;

switch lower( method ),
case {'leastsquares','least-square','leastsquare'},
        L = [];
        c = mean( x, 1 );
case {'minellipse','min-ellipse'},
    [c, L, I] = i_EllipseOptimize( x, @i_MinVolumeObj, @i_MinVolumeConstraints, options );
case {'max-ellipse'}, % doesn't work
    [c, L, I] = i_EllipseOptimize( x, @i_MaxVolumeObj, @i_MaxVolumeConstraints, options );
otherwise
    error(message('mbc:xregfindcenter:InvalidArgument'));
end

c=c(:)';

%--------------------------------------------------------------------------
function [center, L, indices] = i_EllipseOptimize( ...
    points, objective, constraint, options )

if size(points,1)==1
   center = points;
   L = diag(Inf(size(points)));
   indices = 1;
   return
end

[junk, d] = size( points );

% use the middle of the extent of the data at the initial guess for the
% center.
cmin = min( points, [], 1 ); 
cmax = max( points, [], 1 ); 
center0 = 0.5 * (cmax + cmin); 

% Ensure that there is sufficient distance between the upper & lower bounds
i = (cmax - cmin) < 2e-8;
if any( i ),
    % We make sure that we are least subtracting/adding eps so that the
    % values are actually changing
    cmin(i) = cmin(i) - max( 1e-7, eps( cmin(i) ) );
    cmax(i) = cmax(i) + max( 1e-7, eps( cmax(i) ) );
end

% Initial guess for ellipse
L0 = diag( 1./(sqrt( 2 )*(cmax - cmin)) ); 

% set lower and upper bounds
lower_L = min( L0, sqrt( eps ) ) * eye( d ) - tril( repmat( 1e+3, d, d ), -1 );
lb = i_EllipseToColumn( cmin, lower_L );
ub = i_EllipseToColumn( cmax, repmat( 1e+3, d, d ) );


% Initial guess as a column vector
x0 = i_EllipseToColumn( center0, L0 );

% perform optimization
[xstar, junkFvel, junkFlag, JunkOutput, lambda] = fmincon( ...
    objective, x0, [], [], [], [], lb, ub, constraint, options, points );

% Get the indices of the points on the edge of the ellipse
indices = find( lambda.ineqnonlin );

% interpret solution
[center, L] = i_ColumnToEllipse( xstar, d );

%--------------------------------------------------------------------------
function f = i_MinVolumeObj( X, points )
% objective function for minimizing volume
[junk, d] = size( points );

% [center, L] = i_ColumnToEllipse( col, d )
% f = prod( diag( L ) ); % = det( L )

% the diag of the matrix is stored in the last d entries of X
f = -prod( X(end-d+1:end) ); 

%--------------------------------------------------------------------------
function [g, geq] = i_MinVolumeConstraints( X, points )
% Constraint function for minimizing volume
%
% The ellipse is given by (x-c)'*M*(x-c) = 1 and all data points must
% satisfy (x-c)'*M*(x-c) - 1 < 0 where x is data point (column), c is the
% center of the ellipse (column) and M = L'*L is a positive definite
% matrix.

[n, d] = size( points );
[center, L] = i_ColumnToEllipse( X, d );

g = (points - center(ones( n, 1 ),:)) * L';
g = sum( g .* g, 2 ) - 1;

% There are no equality constraints
geq = [];

%--------------------------------------------------------------------------
function f = i_MaxVolumeObj( X, points )
% objective function for maximizing volume
f = -i_MinVolumeObj( X, points );

%--------------------------------------------------------------------------
function [g, geq] = i_MaxVolumeConstraints( X, points )
% constraint function for maximizing volume
[g, geq] = i_MinVolumeConstraints( X, points );
g = -g;

%--------------------------------------------------------------------------
function col = i_EllipseToColumn( center, L )
% Puts the ellipse information into a column vector
% Inverse operation to i_ColumnToEllipse
d = size( center, 2 );
col = zeros( d + 0.5 * d * (d + 1), 1 );

col(1:d) = center;

j = d+1;
for i = 1:d,
    col(j) = diag( L, -d + i );
    j = [ j + i, j(end) + i + 1];
end

%--------------------------------------------------------------------------
function [center, L] = i_ColumnToEllipse( col, d )
% Form the ellipse parameters from the column vector
% Inverse operation to i_EllipseToColumn
center = col(1:d)';
L = zeros( d );
j = d+1;
for i = 1:d,
    L = L + diag( col(j), -d + i );
    j = [ j + i, j(end) + i + 1];
end

%--------------------------------------------------------------------------
% EOF
%--------------------------------------------------------------------------