www.gusucode.com > mbcdata 工具箱 matlab 源码程序 > mbcdata/@cgmathsobject/private/pExtrapolateSpecialCases.m

    function [V, Done] = pExtrapolateSpecialCases(X, Y, Z, Xi, Yi)
%PEXTRAPOLATESPECIALCASES  Detect and perform special case extrapolations
%
%  [T, DONE] = PEXTRAPOLATESPECIALCASES(X, Y, Z, XI, YI) extrapolates from
%  a defined set of values (X, Y, Z) to generate the values at the points
%  (XI, YI).  This routine is a helper function that only looks for and
%  performs 6 special cases:
%
%  1. X, Y, Z specified as empty, in which case we set the output to be
%     empty.
%
%  2. Only one point specified, in which case we set the output to a
%     constant matrix with this value.
%
%  3. Degenerate case, where all (X, Y) points are the same. In this case 
%     we set the output to a constant matrix with the value being the
%     average of the Z-values.
%
%  4. Colinear data is given, in which case we will produce a surface
%     obtained by translating this line along a vector perpendicular to the
%     direction of the line and the z axis.
%
%  5. Three non colinear points are given, in which case they define a
%     plane and we use this to provide our values.
%
%  6. The same number of points as there are evaluation points is
%     given, and every point lies exactly on a different point, i.e. every
%     point has a single exact data value.

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


% Initialize output
V = [];
Done = false;

% Number of points to be evaluated at
Vsize = numel(Xi);

% First up, if we have not been supplied with data, then return empty.
if isempty(X)
    V = [];
    Done = true;
    return
end

% Now, what if we only have 1 trustworthy point - then initialise the
% matrix to this number.
if length(X)==1
    V = ones(Vsize, 1)*Z;
    Done = true;
    return
end

% Next, what if N points are specified but they are all duplicates of the
% same point. Initialise the matrix to the average of the Z-values at this
% point.
uniqueX = X(1);
uniqueY = Y(1);
degenerateCase = true;
for i = 2:length(X)
    if X(i) ~= uniqueX || Y(i) ~= uniqueY
        degenerateCase = false;
        break
    end
end
if degenerateCase
    Z = mean(Z);
    V = ones(Vsize, 1)*Z;
    Done = true;
    return
end

% Now to handle the colinear case when all the trust worthy data is on a
% line.  We assume the data is linear, and then hit it with some tests to
% see if this is true
linearflag = true;
if length(X)>2
    % Test for colinearity by checking singular values from svd
    % decomposition.  Scaled values are used for numerical stability.
    sPts = [i_scale(X), i_scale(Y)].';
    linearflag = isColinear(sPts);
end

if linearflag % Colinear data
    % For colinear data we will translate the line defined by the data in a
    % direction normal to both the line and the z axis to generate a
    % surface. In practice what this means is that we will have a
    % coordinate s which represents distance along the line, we generate an
    % s value for every point in the matrix, and then we use extrinterp to
    % fill out the table.

    % This extrapolation uses the first points to define the line in (X,
    % Y). We should pick the first two distinct points.
    XY = [X Y];
    tol = sqrt(eps)*max(mean(abs(XY),1),1);
    
    XY = curvefitlib.internal.uniqueWithinTol(XY,tol);
    X1 = [XY(1, 1); XY(1, 2); 0];
    X2 = [XY(2, 1); XY(2, 2); 1];

    % plane is of the form Ax+By+Cz = d, where z gives us our s coordinate, 0 at X1, and 1 at X2;
    dX = X1-X2;
    A = -dX(1)*dX(3);
    B = -dX(2)*dX(3);
    C = dX(1)^2+dX(2)^2;
    D = [A B C]*X1;

    % This gives us the s coordinates
    Seval = (D-A*Xi-B*Yi)/C;
    Sin = (D-A*X-B*Y)/C;
    % make sure that new coordindates are sorted so linear extrapolation
    % works
    [Sin,ind]= unique(Sin);
    V = extinterp1(Sin,Z(ind),Seval); 
    Done = true;
    return
end


% Now to handle the coplanar case when there are 3 points that are not on a
% line
if length(X) == 3
    Xmat = [X,Y,ones(3,1)]';
    A = Z'/Xmat;

    V = A(1)*Xi+A(2)*Yi+A(3);
    Done = true;
    return
end


% Now handle case when all output values have had a data value specified
% exactly on them.  For this to be the case, all [X, Y] values must be
% [Xi, Yi] values, there must be no repeated Z values and there must be
% exactly Vsize Z values.
if numel(X)==Vsize
    [uniqueXY, idxData] = unique([X, Y], 'rows');
    [uniqueXiYi, idxPts] = unique([Xi, Yi], 'rows');  
    if size(uniqueXY, 1) == Vsize && isequal(uniqueXY, uniqueXiYi)
        V = zeros(Vsize, 1);
        V(idxPts) = Z(idxData);
        Done = true;
    end
end



% Helper functions that do scaling of the X/Y values
function Out = i_scale(In, mn,mx)
if nargin<2
    [mn,mx] = i_getScale(In);
end
Out = (In-mn)/(mx-mn);

function [mn,mx] = i_getScale(In)
mx = max(In);
mn = min(In);
if mx == mn
    mn = mn-1;
    mx = mx+1;
end