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