function [center, U, obj_fun] = hcm(data, cluster_n, options)
%HCM finds clusters by hard C-means clustering.

%Carson KAO, May 12, 1997

%======================
%stop_impro = 0;
%======================
if nargin ~= 2 & nargin ~= 3,
	error('Too many or too few input arguments!');
end

% Change the following to set default options
default_options = [1000;	% max. number of iteration
		1e-5;	% min. amount of improvement
		1];	% info display during iteration 

if nargin == 2,
	options = default_options;
else
	% If "options" is not fully specified, pad it with default values.
	if length(options) < 3,
		tmp = default_options;
		tmp(1:length(options)) = options;
		options = tmp;
	end
	% If some entries of "options" are nan's, replace them with defaults.
	nan_index = find(isnan(options)==1);
	options(nan_index) = default_options(nan_index);
end

max_iter = options(1);		% Max. iteration
min_impro = options(2);		% Min. improvement
display = options(3);		% Display info or not

obj_fun = zeros(max_iter, 1);	% Array for objective function
%generates initial cluster centers by randomly selecting c points
center = inithcm(cluster_n,data);
%main loop 
for i = 1:max_iter,
    dij = disthcm(center,data);	%compute distance between jpt with center i
    %Determine the membership matrix U of hard C-means
    U = mfuhcm(center,data,dij);
    obj_fun(i) = sum(sum(U.*dij.*dij));	%Compute the cost function value
    if display, 
	fprintf('Iteration count = %d, obj. fun = %f\n', i, obj_fun(i));
    end
    % check termination condition
    if i > 1,
%===============================
	if abs(obj_fun(i) - obj_fun(i-1)) < min_impro, break; end,
%if abs(obj_fun(i) - obj_fun(i-1)) == stop_impro, break; end,
%===============================
    end
    %Update cluster centers
    Gi = (sum(U'))'*ones(1,size(center,2));
    center = U*data./Gi;
end

iter_n = i;	% Actual number of iterations 
obj_fun(iter_n+1:max_iter) = [];
