function [bestIdx, bestCenter, bestSumMinDist] = kmeans(X, K)

    max_iter = 100;
    n = size( X, 1 );
    
%    center = X( randsample( n, K ), : );
    center = X( ceil( n * rand( K, 1 ) ), : );    
    sqrDist = ml_sqrDist( center', X' );
    [ minDist idx ] = min( sqrDist );
    count = hist( idx, 1:K );
    for k = 1 : K
        center( k, : ) = mean( X( idx == k, : ), 1 );
    end
    sqrDist = ml_sqrDist( center', X' );
    
    bestSumMinDist = sum( minDist );
    bestCenter = center;
    bestIdx = idx;    
    
    iter = 0;
    while( 1 )
        
        % first phrase
        while( 1 )
            % reindexing
            regmat = repmat( count' ./ ( count' + 1 ), [ 1 n ] );
            index_in = sub2ind( [ K n ], idx, 1 : n );
            regmat( index_in ) = 1;
%{            
            for k = 1 : n
                if count( idx( k ) ) > 1
                    regmat( idx( k ), k ) = count( idx( k ) ) / ( count( idx( k ) ) - 1 );
                end
            end
%}            
            coff_in = count' ./ max( ( count' - 1 ), 1 );
            index_in = sub2ind( [ K, n ], idx, 1 : n );
            regmat( index_in ) = coff_in( idx );

            sqrDist = sqrDist .* regmat;
            [ minDist idx ] = min( sqrDist );

            % recentering
            count = hist( idx, 1:K );
            for k = 1 : K
                if count( k ) == 0
                    while( 1 )
                        [ maxDist maxPosi ] = max( minDist );
                        minDist( maxPosi ) = -1;
                        if count( idx( maxPosi ) ) > 1 break
                        end
                    end
                    count( k ) = count( k ) + 1;
                    count( idx( maxPosi ) ) = count( idx( maxPosi ) ) - 1;
                    idx( maxPosi ) = k;
                end
            end
            for k = 1 : K
                center( k, : ) = mean( X( idx == k, : ), 1 );
            end
            
            sqrDist = ml_sqrDist( center', X' );
            [ minDist idx_next ] = min( sqrDist );        
            sumMinDist = sum( minDist );
            if sumMinDist < bestSumMinDist
                bestSumMinDist = sumMinDist;
                bestCenter = center;
                bestIdx = idx;
                iter = iter + 1;
                if ( iter >= max_iter ) break;
                end
            else
                break;
            end
            idx = idx_next;
        end
        
        if ( iter >= max_iter ) break;
        end
        
        % second Phrase
        updated = false;
        sqrDist = ml_sqrDist( bestCenter', X' );
        [ minDist idx ] = min( sqrDist );
        clusterDist = zeros( K, 1 );
        for k = 1 : K
            clusterDist( k ) = sum( minDist( idx == k ) );
        end
        for k = 1 : K
            if sum( idx == k ) <= 1 continue;
            end
            sqrDistrow = sqrDist( k, : );
            sqrDistrow( idx ~= k ) = -1;
            [ maxDist maxPosi ] = max( sqrDistrow, [], 2 );
            sqrDist( k, maxPosi ) = Inf;
            [ minDist minPosi ] = min( sqrDist( :, maxPosi ) );
            idx( maxPosi ) = minPosi;
            newCenter1 = mean( X( idx == k, : ), 1 );
            newCenter2 = mean( X( idx == minPosi, : ), 1 );
            newDist1 = sum( ml_sqrDist( newCenter1', X( idx == k, : )' ) );
            newDist2 = sum( ml_sqrDist( newCenter2', X( idx == minPosi, : )' ) );
            if newDist1 + newDist2 < clusterDist( k ) + clusterDist( minPosi ) - 0.0001
                clusterDist( k ) = newDist1;
                clusterDist( minPosi ) = newDist2;
                updated = true;
            else
                idx( maxPosi ) = k;
            end
        end
        if ~updated break;
        end
        
        count = hist( idx, 1:K );
        for k = 1 : K
            center( k, : ) = mean( X( idx == k, : ), 1 );
        end
        sqrDist = ml_sqrDist( center', X' );
        [ minDist next_idx ] = min( sqrDist );
        sumMinDist = sum( minDist );
        bestSumMinDist = sumMinDist;
        bestCenter = center;
        bestIdx = idx;
    end
end