www.gusucode.com > mbcdesign 工具箱 matlab 源码程序 > mbcdesign/@constar/private/ray_casting.m

    function [Ind, R] = ray_casting( X, R, S, center )
%RAY_CASTING boundary identification using the ray casting technique
%   RAY_CASTING(X) is a list of the indices of the points on the boundary of 
%   a superset of those data points given in the rows of X. 
% 
%   RAY_CASTING(X,R) uses spheres of radius R.
%   RAY_CASTING(X,'Auto') automatically estimates the required radius for 
%   the spheres by using the separation distance.
%   [I,R] = RAY_CASTING(X,'Auto') returns the estimated radius in R.
%
%   RAY_CASTING(X,R,S), where S is a scalar, uses S directions. 
%   RAY_CASTING(X,R,S), where S is a matrix with as many columns as X, uses 
%   the directions given by the rows of S. 
%   RAY_CASTING(X,R,'Data'), uses the data to generate the directions.
%
%   RAY_CASTING(X,R,S,C), where C is a row vector with the same number of 
%   columns as X, specifies that C be used as the center of the sphere for 
%   the rays. 
%   RAY_CASTING(X,R,S,'LeastSquares') determines the center this sphere by 
%   the point that approximates the data best in the least square sense. 
%   This is the default. 
%   RAY_CASTING(X,R,S,'MinEllipse') uses as the center of the sphere the 
%   center of the minimal volume bounding ellipsoid.
%
%   See also: SEPARATION, FIND_CENTER.
                
%  Copyright 2000-2011 The MathWorks, Inc. and Ford Global Technologies, Inc.



narginchk(1,4);
if nargin < 2,
    R = 'Auto';
end
if nargin < 3,
    S = 'Data';
end
if nargin < 4,
    center = 'LeastSquares';
end

if strcmpi( R, 'auto' ),
    R = separation( X );
end

[n,d] = size( X );

% if this is a 1-d problem, return min and max indices
if d == 1,
    [~, i] = min( X );
    [~, j] = max( X );
    Ind = [ i; j ];
    return
end
    
% find a center
if ischar( center )
    center = xregfindcenter( X, center );
else
    % check the size of center
    if any( size( center ) ~= [1, d] ),
        error(message('mbc:constar:InvalidArgument1'));
    end
end

% subtract center from data points
X = X - center(ones(n,1),:);

% get directions
if strcmpi( S, 'data' ),
    S = X;
elseif numel( S ) == 1,
    switch d,
    case 2,
        t = linspace( -pi, pi, S+1 );
        t = t(1:S)';
        S = [ sin( t ), cos( t ) ];
    %case 3,
    %    S = partsphere( S );
    otherwise,
        % for the arbitrary case, use random sample-- 
        %     need a more general version of partsphere
        S = 1 - 2*rand( S, d );
    end
end

% inner products and norms
XX = sum( X.^2, 2 );
VV = sum( S.^2, 2 );

% If any of the "directions" coincide with the center then we remove them
% from the analysis. These "directions" don't provide us with an actual
% direction.
i = (VV <= eps );
S(i,:) = [];
VV(i) = [];
N = size( S, 1 );

% Main computation
t_max(1:N)= -Inf;
t_min(1:N)= Inf;
i_min= zeros(1,N);
i_max= zeros(1,N);

R2= R^2;
for j = 1:N
    % rows of S
    sj= 2*S(j,:).';
    vj= 4*VV(j);
    for i=1:n
        % rows of X
        p= X(i,:)*sj;
        
        D = p*p - vj*(XX(i)-R2);
        if D >= 0
            tmpMax = ( p + sqrt(D) )./ (2*VV(j)) ;
            tmpMin = ( p - sqrt(D) )./ (2*VV(j)) ;
            if tmpMax > t_max(j)
                t_max(j)= tmpMax;
                i_max(j)= i;
            end
            if tmpMin < t_min(j)
                t_min(j) = tmpMin;
                i_min(j) = i;
            end
        end
    end
end

% do the output thing
Ind = union( i_max, i_min );
if ~isempty(Ind) && Ind(1)==0
    % remove zero point
    Ind= Ind(2:end);
end

return




function D= separation(X)

D= 0;
N= size(X,1);
for i=1:N-1
    x= X(i,:);
    h= Inf;
    for j= i+1:N
        e= X(j,:)-x;
        h= min( h , e*e');
    end
    D= max(D,h);
end

D= sqrt(D);

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