function segment01(k,filename)
% segment01

if nargin<1,k=3;end
if nargin<2,filename='kenzo.jpg';end
I = imread(filename);
figure(1);title('original');imshow(I)
I=double(I);

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

for t=1:50
  h=size(I,1);w=size(I,2);

  % E-step
  
  % calculate errors
  sumE=zeros(h,w); % constant for making exp behave later
  figure(2);title('errors');
  for g=1:k
    % calculate errors (RGB)
    E{g}=zeros(h,w);
    for c=1:3
      e=(I(:,:,c)-m(c,g))/sigma;
      E{g}=E{g}+e.*e;
    end
    sumE=sumE+E{g};
    subplot(1,k,g);
    imagesc(E{g});axis image;colormap(hot)
  end
  cte=sumE/k;
  waitforbuttonpress

  % calculate un-normalized probabilities
  sumq=zeros(h,w);
  for g=1:k
    q{g}=pi(g)*exp(-0.5*(E{g}-cte));
    sumq=sumq+q{g};
  end
  
  % normalize probabilities and show them
  figure(3);title('probabilities');
  for g=1:k
    p{g}=q{g}./sumq;
    subplot(1,k,g);
    imagesc(p{g});axis image;colormap(gray)
  end
  waitforbuttonpress
  
  % create and show index image
  figure(4);title('index');
  best=p{1};
  index=ones(h,w);
  for g=2:k
    P=p{g};
    better=P>best;
    best(better)=P(better);
    index(better)=g;
  end
  imshow(index,lines(k));
  waitforbuttonpress
  
  % M-step
  for g=1:k
    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);
  end
  figure(5);title('average color');
  imshow(index,m'/255);
  
  pi=pi/sum(pi)
  m
  waitforbuttonpress
end