function [C, U, LUT, H] = FastFCMeans(im, c, q, opt)
% Partition N-dimensional grayscale image into c classes using a memory 
% efficient implementation of the fuzzy c-means (FCM) clustering algorithm. 
% Computational efficiency is achieved by using the histogram of image
% intensities during clustering instead of the raw image data.
%
% INPUT:
%   - im  : N-dimensional grayscale image in integer format. 
%   - c   : positive integer greater than 1 specifying the number of
%           clusters. c=2 is the default setting. Alternatively, c can be
%           specified as a k-by-1 array of initial cluster (aka prototype)
%           centroids.
%   - q   : fuzzy weighting exponent. q must be a real number greater than
%           1.1. q=2 is the default setting. Increasing q leads to an
%           increased amount of fuzzification, while reducing q leads to
%           crispier class memberships. Note that while in principle
%           setting q==1 is equivalent to using a classical c-means 
%           algorithm, this setting cannot be used in practice because it
%           produces an infinite exponent in the membership update formula.
%   - opt : optional logical argument used to indicate how to initialize
%           cluster centroids. If opt=true {default} then centroids are
%           initialized are by sampling the intensity range at uniform
%           intervals. If opt=false then the initial centroids are set 
%           using the c-means algorithm. 
%
% OUTPUT  :
%   - C   : 1-by-k array of cluster centroids.
%   - U   : L-by-k array of fuzzy class memberships, where k is the number
%           of classes and L is the intensity range of the input image, 
%           such that L=numel(min(im(:)):max(im(:))).
%   - LUT : L-by-1 array that specifies the intensity-class relations,
%           where L is the dynamic intensity range of the input image. 
%           Specifically, LUT(1) corresponds to the (class) label assigned to 
%           min(im(:)) and LUT(L) corresponds to the label assigned 
%           to max(im(:)). LUT is used as input to 'apply_LUT' function to
%           create a label image.
%   - H   : image histogram. If I=min(im(:)):max(im(:)) are the intensities
%           present in the input image, then H(i) is the number of  
%           pixels/voxels with intensity I(i). 
%
% AUTHOR    : Anton Semechko (a.semechko@gmail.com)
%


% Default input arguments
if nargin<2 || isempty(c), c=2; end
if nargin<3 || isempty(q), q=2; end
if nargin<4 || isempty(opt), opt=true; end

% Basic error checking
if nargin<1 || isempty(im)
    error('Insufficient number of input arguments')
end
msg='Revise variable used to specify class centroids. See function documentation for more info.';
if ~isnumeric(c) || ~isvector(c)
    error(msg)
end
if numel(c)==1 && (~isnumeric(c) || round(c)~=c || c<2)
    error(msg)
end
if ~isnumeric(q) || numel(q)~=1 || q<1.1
    error('3rd input argument (q) must be a real number > 1.1')
end
if ~islogical(opt) || numel(opt)>1
    error('4th input argument (opt) must a Boolean')
end    

% Check image format
if isempty(strfind(class(im),'int'))
    error('Input image must be specified in integer format (e.g. uint8, int16)')
end

% Intensity range
Imin=double(min(im(:)));
Imax=double(max(im(:)));
I=(Imin:Imax)';

% Initialize cluster centroids
if numel(c)>1 % user-defined centroids
    C=c;
    opt=true;
else % automatic initialization
    if opt
        dI=(Imax-Imin)/c;
        C=Imin+dI/2:dI:Imax;
    else
        [C,~,H]=FastCMeans(im,c);
    end
end

% Compute intensity histogram
if opt
    H=hist(double(im(:)),I);
    H=H(:);
end
clear im

% Update fuzzy memberships and cluster centroids
dC=Inf;
while dC>1E-3
    
    C0=C;
    
    % Distance to the centroids
    D=abs(bsxfun(@minus,I,C));
    D=D.^(2/(q-1))+eps;
    
    % Compute fuzzy memberships
    U=bsxfun(@times,D,sum(1./D,2));
    U=1./(U+eps);
    
    % Update the centroids
    UH=bsxfun(@times,U.^q,H);
    C=sum(bsxfun(@times,UH,I),1)./sum(UH,1);
    C=sort(C,'ascend'); % enforce natural order
    
    % Change in centroids 
    dC=max(abs(C-C0));
    
end

% Defuzzify
[~,LUT]=max(U,[],2);