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

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

%   Copyright 2006-2007 The MathWorks, Inc.
    
if nargout > 3
    szX = size(X);
    nX = max(szX);
    A = speye(nX);
end

keep = true(size(X));
for n = 1:length(X)
    if keep(n)
        AvX = X(n);
        AvY = Y(n);
        AvZ = Z(n);
        nAv = 1;

        % Going to use a logical vector to indicate duplicates. This is
        % because it is quicker to build up a logical vector and index
        % into the sparse A at the end of the inner loop than indexing into
        % A inside the loop.
        if nargout > 3
            idxDup = false(1, length(X));
            idxDup(n) = true;
        end
        
        % Look for close points that haven't been clustered already
        for m = (n+1):length(X)
            Xd = X(m)-AvX;
            Yd = Y(m)-AvY;

            if keep(m) ...
                    && Xd<tol && Xd>-tol...
                    && Yd<tol && Yd>-tol

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

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

                nAv = nAv+1;

                AvX = (AvX + X(m))/nAv;
                AvY = (AvY + Y(m))/nAv;
                AvZ = (AvZ + Z(m))/nAv;
                
                  if nargout > 3
                      idxDup(m) = true;
                  end
            end
        end

        X(n) = AvX;
        Y(n) = AvY;
        Z(n) = AvZ;
        
        if nargout > 3
           A(n, idxDup) = 1/nAv;
        end
    end
end

X = X(keep);
Y = Y(keep);
Z = Z(keep);
if nargout > 3
    A = A(keep, :);
    if szX(1) == 1
        % Transpose A for row vectors
        A = A';
    end
end