function [skin1, skin2, optimalThreshold] = SegmentSkin(filename, bmean, rmean, brcov)
% Assume the skinmodel.m is run
% Produce two images, skinlikelihood greyscale image, skin1
% and skin segment binary image, skin2

im = imread(filename);
imycbcr = rgb2ycbcr(im);

dim = size(im);
skin1 = zeros(dim(1), dim(2));
for i = 1:dim(1)
   for j = 1:dim(2)
      cb = double(imycbcr(i,j,2));
      cr = double(imycbcr(i,j,3));
      x = [(cb-bmean); (cr-rmean)];
      skin1(i,j) = exp(-0.5* x'*inv(brcov)* x);
   end
end

lpf= 1/9*ones(3);
skin1 = filter2(lpf,skin1);
skin1 = skin1./max(max(skin1));

% Adaptive Thresholding
previousSkin2 = zeros(i,j);
changelist = [];
for threshold = 0.55:-0.1:0.05
   skin2 = zeros(i,j);
   skin2(find(skin1>threshold)) = 1;
   change = sum(sum(skin2 - previousSkin2));
   changelist = [changelist change];
   previousSkin2 = skin2;
end

% Finding the optimal threshold
[C, I] = min(changelist);
optimalThreshold = (7-I)*0.1
skin2 = zeros(i,j);
skin2(find(skin1>optimalThreshold)) = 1;

figure(1)
imshow(skin1, [0 1]);
figure(2)
imshow(skin2, [0 1]);