% bare bones

% init models
sigma=20;
m(:,g)=128+50*randn(3,1);
pi(g)=1/k;

% E-step
  
% calculate errors
E{g}=zeros(h,w);
e=(I(:,:,c)-m(c,g))/sigma;
E{g}=E{g}+e.*e;

% calculate un-normalized probabilities
q{g}=pi(g)*exp(-0.5*E{g});
  
% normalize probabilities 
p{g}=q{g}./sumq;
  
% M-step
P=p{g};
R=P.*I(:,:,1);
G=P.*I(:,:,2);
B=P.*I(:,:,3);
pi(g)=sum(P(:));
m(:,g) = [sum(R(:));sum(G(:));sum(B(:))]/pi(g);

pi=pi/sum(pi)