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