--- a
+++ b/libraries/lib_FastCMeans/FastFCMeans.m
@@ -0,0 +1,124 @@
+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);
+