function Umap=FM2map(im,U,H)

% Unpack fuzzy-membership functions to produce membership maps.

% 

% INPUT:

%   - im    : N-dimensional grayscale image in integer format. 

%   - U     : L-by-c array of fuzzy class memberships, where c is the 

%             number of classes and L is the intensity range of the input 

%             image, such that L=numel(min(im(:)):max(im(:))). See 

%             'FastFCMeans' function for more info.

%   - H     : image histogram returned by 'FastFCMeans' function.

%

% OUTPUT:

%   - Umap  : membership maps. Umap has the same size as the input image

%             plus an additional dimension to account for c classes. For

%             example, if im is a 2D M-by-N image then U will be 

%             M-by-N-by-c array where U(:,:,i) is a membership map for the

%             i-th class.

%

% AUTHOR    : Anton Semechko (a.semechko@gmail.com)

%

if nargin<3 || isempty(H)

    % Intensity range

    Imin=double(min(im(:)));

    Imax=double(max(im(:)));

    I=(Imin:Imax)';

    % Intensity histogram

    H=hist(double(im(:)),I);

    H=H(:);

end

% Unpack memberships

Umap=zeros(sum(H),size(U,2));

i1=1; i2=0;

for i=1:numel(H)

    i2=i2+H(i);

    Umap(i1:i2,:)=repmat(U(i,:),[H(i) 1]);

    i1=i2+1;

end

% Find positional mapping

[~,idx]=sort(im(:), 'ascend');

[~,idx]=sort(idx(:),'ascend');

% Reshape membership maps to match image dimensions

Umap=reshape(Umap(idx,:),[size(im) size(U,2)]);