%function computes a binary (pixel-wise) mask of an input image using the EM
%algorithm, median filter and image closing, which's size is given as an input value. The
%mask should highlight the foreground. medFiltSize can be scalar, defining
%a square median filter or a 2-element vector defining a rectangular median
%filter. seSize can be scalar, defining a square structuring element or a
%2-element vector defining a rectangular structuring element. The outputs
%mask and maskedI have the same amount of rows and columns as the input
%image. Note that for the 4-by-4 downsampled versions of the images used in
%our experiments, a good choice for medFiltSize=[20, 20]=seSize.

function [mask,maskedI]=segmentForeBackEmMedClosing(I,medFiltSize,seSize)

%two clusters: foreground and background
noClusters=2;
if numel(seSize)==1
   seSize=[1 1]*seSize; 
elseif numel(seSize)~=2
   error('Size of structuring element for closing is a vector with more than two elements.')
end
if numel(medFiltSize)==1
   medFiltSize=[1 1]*medFiltSize; 
elseif numel(medFiltSize)~=2
   error('Size of median filter is a vector with more than two elements.')
end
if length(size(I))==3
    %if yes, compute luminance
    Iy=rgb2ycbcr(double(I));
    Iy=Iy(:,:,1);
elseif length(size(I))==2
    %if not, use original image as luminance
    Iy=I;
else
    error('Input image has more than three or less than two dimensions.')
end
[M,N]=size(Iy);
%gaussian mixtures, image
classyData=double(Iy(:));
errFlag=1;
while errFlag==1
        objIm = gmdistribution.fit(classyData,noClusters) ;
        gaussMixIxIm = cluster(objIm,classyData);
        if objIm.mu(1)~=objIm.mu(2)
            errFlag=0;
        else
            errFlag=1;
            disp('Warning: EM has not converged and will be started again. Risk of an infinite loop.')
        end
end
% objIm = gmdistribution.fit(classyData,noClusters) ;
% gaussMixIxIm = cluster(objIm,classyData);
%the smaller sigma refers to background
[~,minIx]=min(objIm.Sigma);
%for mask, reduce labels to 0 and 1
mask = abs(gaussMixIxIm-minIx);
%reshape to original image size
mask=reshape(mask,[M,N]);
% figure; imshow(mask);
sE = strel('rectangle',seSize);
%do image closing on mask
mask=double(imclose(mask,sE));
%do median filter on mask
mask=double(medfilt2(uint8(mask),medFiltSize));
%multiply mask with image, resulting in a foreground-only image
maskedI=uint8(mask).*uint8(Iy);