function segmented_image=EM_test(gray_vals,mask_in)
%Expectation-Maximization image segmentation algorithm
%segmented_image = EM(gray_vals,mask_in)
%An image is segmented using the binary values of mask_in
%as initial guesses for the values of the segmentation.

%k is the number of different layers of segmentation
k=2;
[row,col]=size(gray_vals);
[row_check,col_check]=size(mask_in);
if row~=row_check | col~=col_check
    error('Size of image does not match size of mask!')
end
gray_vals=reshape(gray_vals,1,row*col);
mask_in=reshape(mask_in,1,row*col);
%I is the "support map" described in Forsyth and Ponce
I=zeros(k,row*col);
I(1,:)=mask_in;
mask_in=-mask_in+1;
mask_in=mask_in./(k-1);
for i=2:k
    I(i,:)=mask_in;
end;
I(2:k,:)=I(2:k,:)+rand(size(I(2:k,:)));
for i=2:k
    I(i,:)=I(i,:).*mask_in;
end
for i=1:k
    I(i,:)=I(i,:)./sum(I);
end

firsttime=1;
I_prev=I;
counter=0
%Determine initial alpha, and means and variances
alpha=sum(I,2);
for i=1:k
    %This line takes the mean of each region neglecting pixels belonging to
    %other regions
    mean_vals_temp(i,:)=fliplr(sort(I(i,:).*(gray_vals)));
    mean_vals_temp(i,:)=mean(mean_vals_temp(i:max(find(mean_vals_temp(i,:)))),2);

end

%The mean_vals variable is of smaller length than the mean_vals_temp (this is
%just a book keeping change var_vals_temp has many copies of the same
%information)
mean_vals=mean_vals_temp(:,1);
for i=1:k
    %This line takes the variance of each region neglecting pixels
    %belonging to the other regions
    var_vals_temp(i,:)=fliplr(sort(I(i,:).*(gray_vals)));
    var_vals_temp(i,:)=std(var_vals_temp(i:max(find(var_vals_temp(i,:)))),0,2);

end

%The var_vals variable is of smaller length than the var_vals_temp (this is
%just a book keeping change var_vals_temp has many copies of the same
%information)
var_vals=var_vals_temp(:,1);
mean_past=-1*ones(size(mean_vals));
var_past=-1*ones(size(mean_vals));

flag=0;

%Sometimes the entire region is a person, so it's pointless to try to
%segment further and actually produces errors, so this checks for that
%condition.
for i=1:k
    if isnan(mean_vals(i))
        flag=1;
    end
end

%Stopping condition is the maximum change in the mean value
while (max(abs(mean_past-mean_vals))>.5 & flag==0)
    mean_past=mean_vals;
    var_past=var_vals;
    
    %Update the support map
    for i=1:k
        I(i,:)=alpha(i)/sqrt(2*pi*var_vals(i))*exp(-(gray_vals-mean_vals(i)).^2/(2*var_vals(i)));
    end
    Isum=sum(I);
    for i=1:k
        I(i,:)=I(i,:)./Isum ;
    end

    %Now update the params

    alpha=mean(I,2);
    for i=1:k
        mean_vals(i)=sum(gray_vals.*I(i,:))/sum(I(i,:));
    end
    for i=1:k
        var_vals(i)=sum(((gray_vals-mean_vals(i)).^2).*I(i,:))/sum(I(i,:));
    end

    %If the entire region is one segment or the other, break out.
    if isnan(mean_vals(i))
        break
    end
end

%The output has the segmented image as ordered by 1, 2, 3... for each
%segment
[vals,index]=sort(I);
for i=1:k
    segmented_image(:,:,i)=reshape(index(i,:),row,col);
end
segmented_image=segmented_image(:,:,1);
