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

    function [X, Z, A] = removeDuplicatePts1(X, Z, tol)
%REMOVEDUPLICATEPTS1 Remove duplicate points from a list
%
%   [X, Z] = REMOVEDUPLICATEPTS1(OBJ, X, Z, TOL) removes duplicates from
%   the 1-d points, X. All duplicates within the tolerance, TOL, are
%   replaced by one point centered at the mean of the duplicates. The
%   output values at the 1-d points, Z, are also modified to reflect the
%   new 1-d points.
%
%   [X, Z, A] = REMOVEDUPLICATEPTS1(OBJ, X, Z, TOL) returns a matrix,
%   A, which maps [X, Z] with duplicates onto the unique [X, Z].
%
%   See also CGMATHSOBJECT/REMOVEDUPLICATEPTS2

%   Copyright 2006-2007 The MathWorks, Inc.
 
% Return immediately if no data is specified
if isempty(X) 
    return
end

% Sort the data. Perform a sort on the index so we can restore the data to
% the original order at the end.
[Xsort, idx] = sort(X);
[unused, invidx] = sort(idx);
Zsort = Z(idx);

% Initialisation for linear map return
if nargout > 2
    szX = size(X);
    nX = max(szX);
    A = speye(nX);
end

% Perform one pass through the sorted list, clustering points that are
% within tolerance
keep = true(size(Xsort));
AvX = Xsort(1);
AvZ = Zsort(1);
nAv = 1;
currentUniquePt = 1;
for m = 2:length(Xsort)
    Xd = Xsort(m)-AvX;

    if Xd<tol && Xd>-tol

        % This point is too close to the current cluster's average
        keep(m) = false;

        % Update cluster averages
        AvX = nAv*AvX;
        AvZ = nAv*AvZ;

        nAv = nAv+1;

        AvX = (AvX + Xsort(m))/nAv;
        AvZ = (AvZ + Zsort(m))/nAv;

    else
        if nargout > 2 && m-currentUniquePt > 1
            A(currentUniquePt, currentUniquePt:m-1) = 1/nAv;
        end
        Xsort(currentUniquePt) = AvX;
        Zsort(currentUniquePt) = AvZ;
        currentUniquePt = m;
        AvX = Xsort(m);
        AvZ = Zsort(m);
        nAv = 1;
    end
end

% Deal with the case where the last point is a duplicate
if currentUniquePt ~= m
    if nargout > 2 
        A(currentUniquePt, currentUniquePt:m) = 1/nAv;
    end
    Xsort(currentUniquePt) = AvX;
    Zsort(currentUniquePt) = AvZ;
end
    
% Rearrange X, Z and A back into the order they were given to the function
X = Xsort(invidx);
X = X(keep(invidx));
Z = Zsort(invidx);
Z = Z(keep(invidx));

% Return linear map if required
if nargout > 2
    A = A(invidx(keep(invidx)), invidx);
    if szX(1) == 1
        % Transpose A for row vectors
        A = A';
    end
end