skinem2.m

Matlab skin detector。运用高斯混合模型训练的到人的皮肤颜色分布。用于皮肤和人脸检测。

% Computes the GMM of skin and background using EM of the passed array of
% color images and starting with Rehg and Jones skin color mixture of
% Gaussian model, from "Statistical Color Models with Application to Skin Detection"
% 
%function labelims = jrmogskin(colorims)
%  colorims is matrix of color images to be EMed with each colorim as the following:
%     
%     colorims(:,:,1, X) = red pixel values of image X, 0 to 255
%     colorims(:,:,2, X) = green pix values of image X, 0 to 255
%     colorims(:,:,3, X) = blue pix values of image X,  0 to 255
%  iterations is the number of times to do EM  
%
%  skinMOG is the EMed skin GMMs
%  nonskinMOG is the EMed background GMMs
%
%  base GMM application by:
%  Robert Collins, Penn State University, 2005
%
% added EM extension authored by:
% Anthony Dotterer

function [skinMOG, nonskinMOG] = skinEM(colorims, iterations)

%format of Mixture of Gaussian models:
%one gaussian component per row
%    rmean gmean bmean rvar gvar bvar weight
skinMOG = [...
  73.53  29.94  17.76   765.40  121.44  112.80  0.0294;
  249.71  233.94  217.49   39.94  154.44  396.05  0.0331;
  161.68  116.25  96.95   291.03  60.48  162.85  0.0654;
  186.07  136.62  114.40   274.95  64.60  198.27  0.0756;
  189.26  98.37  51.18   633.18  222.40  250.69  0.0554;
  247.00  152.20  90.84   65.23  691.53  609.92  0.0314;
  150.10  72.66  37.76   408.63  200.77  257.57  0.0454;
  206.85  171.09  156.34   530.08  155.08  572.79  0.0469;
  212.78  152.82  120.04   160.57  84.52  243.90  0.0956;
  234.87  175.43  138.94   163.80  121.57  279.22  0.0763;
  151.19  97.74  74.59   425.40  73.56  175.11  0.1100;
  120.52  77.55  59.82   330.45  70.34  151.82  0.0676;
  192.20  119.62  82.32   152.76  92.14  259.15  0.0755;
  214.29  136.08  87.24   204.90  140.17  270.19  0.0500;
  99.57  54.33  38.06   448.13  90.18  151.29  0.0667;
  238.88  203.08  176.91   178.38  156.27  404.99  0.0749];

nonskinMOG = [...
  254.37  254.41  253.82   2.77  2.81  5.46  0.0637;
  9.39  8.09  8.52   46.84  33.59  32.48  0.0516;
  96.57  96.95  91.53   280.69  156.79  436.58  0.0864;
  160.44  162.49  159.06   355.98  115.89  591.24  0.0636;
  74.98  63.23  46.33   414.84  245.95  361.27  0.0747;
  121.83  60.88  18.31   2502.24  1383.53  237.18  0.0365;
  202.18  154.88  91.04   957.42  1766.94  1582.52  0.0349;
  193.06  201.93  206.55   562.88  190.23  447.28  0.0649;
  51.88  57.14  61.55   344.11  191.77  433.40  0.0656;
  30.88  26.84  25.32   222.07  118.65  182.41  0.1189;
  44.97  85.96  131.95   651.32  840.52  963.67  0.0362;
  236.02  236.27  230.70   225.03  117.29  331.95  0.0849;
  207.86  191.20  164.12   494.04  237.69  533.52  0.0368;
  99.83  148.11  188.17   955.88  654.95  916.70  0.0389;
  135.06  131.92  123.10   350.35  130.30  388.43  0.0943;
  135.96  103.89  66.88   806.44  642.20  350.36  0.0477];

% make all of the mixture's weights sum to 1 or close to 1 
for i = 1:16
   skinMOG(i, 7) = skinMOG(i, 7) / 2;
   nonskinMOG(i, 7) = nonskinMOG(i, 7) / 2;
end;

% get the dimensions
[nr, nc, nb, ni] = size(colorims);
datasize = nr*nc*ni;

% concatenate all of the images into color vectors 
redvals = double(reshape(colorims(:,:,1, :),datasize,1));
grevals = double(reshape(colorims(:,:,2, :),datasize,1));
bluvals = double(reshape(colorims(:,:,3, :),datasize,1));


% Construct a parameter vector that will hold each EM iteration
lamda = zeros(32, 7, 1);
lamda(1:16, :, 1) = skinMOG(:, :);
lamda(17:32, :, 1) = nonskinMOG(: ,:);

% let it iterate for a bit (just to see if this works at all)
for t = 1:iterations
   
   % compute the propabilities of all datapoints for each class
   props = zeros(32, datasize);
   for i = 1:32
      tmp = (redvals - lamda(i, 1, t)).^2 / lamda(i, 4, t) + (grevals - lamda(i, 2, t)).^2 / lamda(i, 5, t) + (bluvals - lamda(i, 3, t)).^2 / lamda(i, 6, t);
      tmp = lamda(i, 7, t) * exp(- tmp / 2) / ((2 * pi)^(3/2) * sqrt(lamda(i, 4, t)*lamda(i, 5, t)*lamda(i, 6, t)));
      props(i, :) = tmp(:, 1)';
   end;
   
   % compute the sum of the propabilities of all datapoints for each class
   propssum = sum(props, 1);
   
   % compute expected classes for all datapoints for each class
   Estep = zeros(32, datasize);
   for i = 1:32
      Estep(i, :) = props(i, :) ./ propssum;
   end;
   
   % compute the sum of the expected classes for all datapoints
   Estepsum = sum(Estep, 2);
   
   % compute the new weights for each class
   newWeights = [];
   for i = 1:32
      newWeights(i) = Estepsum(i) / datasize;
   end;
   
   % compute the new means for each class for each color
   newMeans = zeros(32, 3);
   for i = 1:32
      if (Estepsum(i) ~= 0)
         tmp = sum(Estep(i, :)' .*  redvals);
         newMeans(i, 1) = sum(Estep(i, :)' .*  redvals) / Estepsum(i);
         newMeans(i, 2) = sum(Estep(i, :)' .*  grevals) / Estepsum(i);
         newMeans(i, 3) = sum(Estep(i, :)' .*  bluvals) / Estepsum(i);
      else
         newMeans(i, :) = lamda(i, 1:3, t); % 0 estep sum means 0 weight just pass previous value
      end;
   end;
   
   % compute the new variances for each class for each color
   % NOTE: if the variance drop to zero, because it is too small then use the smallest double value possible
   newVar = zeros(32, 3);
   for i = 1:32
      if (Estepsum(i) ~= 0) 
         newVar(i, 1) = max(sum(Estep(i, :)' .*  (redvals - newMeans(i, 1)).^2) / Estepsum(i), eps);
         newVar(i, 2) = max(sum(Estep(i, :)' .*  (grevals - newMeans(i, 2)).^2) / Estepsum(i), eps);
         newVar(i, 3) = max(sum(Estep(i, :)' .*  (bluvals - newMeans(i, 3)).^2) / Estepsum(i), eps);
      else
         newVar(i, :) = lamda(i, 4:6, t); % 0 estep sum means 0 weight just pass previous value
      end;      
   end;
   
   % setup the next run in the lamda vector
   newLamda = zeros(32, 7, 1);
   for i = 1:32
      newLamda(i, 1, 1) = newMeans(i, 1);
      newLamda(i, 2, 1) = newMeans(i, 2);
      newLamda(i, 3, 1) = newMeans(i, 3);
      newLamda(i, 4, 1) = newVar(i, 1);
      newLamda(i, 5, 1) = newVar(i, 2);
      newLamda(i, 6, 1) = newVar(i, 3);
      newLamda(i, 7, 1) = newWeights(i);
   end;
   lamda(:, :, t+1) = newLamda;
   %lamda(:, :, t+1) = [ newMeans, newVar, newWeights' ];
   
end;

% replace the values of the skin and non-skin models
skinMOG = lamda(1:16, :, t+1);
nonskinMOG = lamda(17:32, :, t+1);


% try to normalize the skin's total probalities to 1 and background's to 1
totalSkinProp = sum(skinMOG(:, 7));
totalNonskinProp = sum(nonskinMOG(:, 7));
for i = 1:16
   skinMOG(i, 7) = skinMOG(i, 7) / totalSkinProp;
   nonskinMOG(i, 7) = nonskinMOG(i, 7) / totalNonskinProp;
end;

return;