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, skin2im = 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);   endendlpf= 1/9*ones(3);skin1 = filter2(lpf,skin1);skin1 = skin1./max(max(skin1));% Adaptive ThresholdingpreviousSkin2 = 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.1skin2 = zeros(i,j);skin2(find(skin1>optimalThreshold)) = 1;figure(1)imshow(skin1, [0 1]);figure(2)imshow(skin2, [0 1]);