meanshift.m

快速高斯变换的程序

function [x, idx] = meanshift(xi, h, method, p, Kc, e, itmax, K)
%MEANSHIFT is an implementation of mean shift method
%
%   input parameters:
%      xi:  source points
%       h:  scale of the kernel
%  method:  'steepdesc', steepest descent method
%           'bfgs', BFGS method
%     dim:  dim = 1, vector is column of xi;
%           dim = 2, vector is row of xi;
%   itmax:  max iterations
%       K:  number of clusters
%   output parameters:
%       x:  output mode points(stationary points)
%     idx:  index to the centers
%       c:  centers

%   $Revision: 1.0 $
%   Changjiang Yang     Nov. 2002

%x = xi;
switch method
    case 'steepdesc'
        xx=xi+randn(size(xi))*.025*range(range(xi));
        figure; plot3(xx(1,:),xx(2,:),xx(3,:),'.','MarkerSize',4);
        hold on;
        tic
        x = kctfgt(xi,h,p,Kc,e,itmax);
        toc
        xx=x+randn(size(x))*.015*range(range(x));
        plot3(xx(1,:),xx(2,:),xx(3,:),'r.','MarkerSize',4);
        hold off;
    case 'bfgs'
        OPTIONS=optimset('LargeScale','off');
        OPTIONS=optimset(OPTIONS,'HessUpdate','bfgs');
        OPTIONS=optimset(OPTIONS,'MaxFunEvals',300);
        OPTIONS=optimset(OPTIONS,'MaxIter',itmax);
        %OPTIONS=optimset(OPTIONS,'Display','iter','TolX',1E-4);
        OPTIONS=optimset(OPTIONS,'Display','off');
        OPTIONS = optimset(OPTIONS,'gradobj','on');
        for i = 1:n,
            x0 = x(:,i);
            x(:,i) = fminunc(@funcgrad,x0,OPTIONS,xi,h); 
        end
    otherwise
        error('Invalid method');
end

if nargout > 1,
    idx = kmeans(x', K, 'rep',3, 'EmptyAction','singleton','disp','final');
    %idx = kcenters(x,K)+1;
end

return;

%%--------------------------------------------------------------%%
% function [f, g] = funcgrad(x, datapts, bandwid)
% %FUNCGRAD gives function vaule and gradient at position x
% % x must be a a column vector
% 
% [DIM_DATA, NUM_DATA] = size(datapts);
% 
% dxih = (x(:,ones(NUM_DATA,1)) - datapts) / bandwid;
% expx = exp(-sum(dxih.^2, 1)/2);
% 
% f = -mean(expx);    % Cost function
% g = sum(repmat(expx,DIM_DATA,1).*dxih, 2); % Compute the gradient
% 
% return;