ATR.zip_ATR_K._atr 算法

function [idx,ctrs,Nc] = my_isodata(X,k,Nc,Tn,Ts,Tc,L,MaxIter) [m,n] = size(X); % X is an m-by-n matrix ctrs = X(1:m,1:Nc); % ctrs is an m-by-k matrix, Nc for initialization idx = zeros(n,1); % idx is an n-dimension vector flag = 1; while(flag) % step 2 dist = zeros(Nc,1); for i = 1:n for j = 1:Nc a = X(:,i);b = ctrs(:,j); c = [a b]'; dist(j) = pdist(c,'cityblock'); end d = find(dist==min(dist)); idx(i) = d(1); end %step 3 i = 1; while(i<=Nc) if length(find(idx==i))<Tn idx(idx==i) = 0; idx(idx>i) = idx(idx>i)-1; Nc = Nc-1; else i = i+1; end end % step 4 % Nc ctrs = zeros(m,Nc); for i = 1:Nc a = X(:,idx==i); ctrs(:,i) = mean(a,2); end % step 5 avgd = zeros(Nc,1);totd = avgd; for i = 1:Nc a = X(:,idx==i);b = ctrs(:,i); len = length(find(idx==i)); dist = zeros(len,1); for j = 1:len c = [a(:,j) b]'; dist(i) = pdist(c,'cityblock'); end avgd(i) = mean(dist); totd(i) = sum(dist); end % step 6 avgD = sum(totd)/n; % step 7 if flag==MaxIter Tc = 0;break; % % step 11 & 12 % Dist = inf*ones(Nc); % for i = 1:Nc-1 % for j = i+1:Nc % a = ctrs(:,i);b = ctrs(:,j); % c = [a b]';d = pdist(c,'cityblock'); % if d<Tc % Dist(i,j) = d; % end % end % end % a = reshape(Dist,1,Nc^2);a = sort(a); % ref = a(L); % % step 13 % [row,col] = find(Dist<=ref);ctrs_new = zeros(m,L); % for i = 1:L % len1 = length(find(idx==row(i))); % len2 = length(find(idx==col(i))); % ctrs_new(:,i) = (ctrs(:,row(i))*len1+ctrs(:,col(i))*len2)... % /(len1+len2); % end % a = unique([row;col]); % ctrs(:,a) = [];ctrs = [ctrs ctrs_new]; % Nc = Nc+L-length(a); elseif Nc<=ceil(k/2) % step 8 & 9 sig = zeros(m,Nc);sigmax = zeros(1,Nc); for i = 1:Nc a = X(:,idx==i)'; sig(:,i) = std(a,1)'; sigmax(i) = max(sig(:,i)); end % step 10 j = 1; for i = 1:Nc len = length(find(idx==i)); if (sigmax(i)>Ts) if ((avgd(i)>avgD)&&(len>2*(Tn+1)))||(Nc<=ceil(k/2)) ctrs = [ctrs zeros(m,1)]; ctrs(:,i) = ctrs(:,i)+0.5*sigmax(i); ctrs(:,Nc+j) = ctrs(:,i)-0.5*sigmax(i); j = j+1; end end end Nc = Nc+j-1; elseif (mod(flag,2)==0)||(Nc>2*k) % step 11 & 12 Dist = inf*ones(Nc); for i = 1:Nc-1 for j = i+1:Nc a = ctrs(:,i);b = ctrs(:,j); c = [a b]';d = pdist(c,'cityblock'); if d<Tc Dist(i,j) = d; end end end a = reshape(Dist,1,Nc^2);a = sort(a); ref = a(L); % step 13 [row,col] = find(Dist<=ref);ctrs_new = zeros(m,L); for i = 1:L len1 = length(find(idx==row(i))); len2 = length(find(idx==col(i))); ctrs_new(:,i) = (ctrs(:,row(i))*len1+ctrs(:,col(i))*len2)... /(len1+len2); end a = unique([row;col]); ctrs(:,a) = [];ctrs = [ctrs ctrs_new]; Nc = Nc+L-length(a); else % step 8 & 9 sig = zeros(m,Nc);sigmax = zeros(1,Nc); for i = 1:Nc a = X(:,idx==i)'; sig(:,i) = std(a,1)'; sigmax(i) = max(sig(:,i)); end % step 10 j = 1;len = length(find(idx==i)); for i = 1:Nc if (sigmax(i)>Ts) if ((avgd(i)>avgD)&&(len>2*(Tn+1)))||(Nc<=ceil(k/2)) ctrs = [ctrs zeros(m,1)]; ctrs(:,i) = ctrs(:,i)+0.5*sigmax(i); ctrs(:,Nc+j) = ctrs(:,i)-0.5*sigmax(i); j = j+1; end end end Nc = Nc+j-1; end flag = flag+1; end end