a = imread('PSWebImageBitMap0001.bmp',  'bmp');

b = double(a(:,:,2));
figure(1)
clf
imshow(b)

[pos1, pos2] = find(b > 228);
X = [pos1, pos2];

k = 3;
[gIdx,c]=k_means(X, k);
rc = round(c);

S = zeros(size(b));
for i = 1:size(c,1);
    S(rc(i,2), rc(i,1)) = 1;
end
[g, NR, SI, TI] = regiongrow(b, S, 11);
[x1,y1] = find(g == 1);
[x2,y2] = find(g == 2);
[x3,y3] = find(g == 3);

figure(4)
clf
plot(y1,x1,'r.',y2,x2,'g.', y3,x3,'b.', rc(:,1),rc(:,2),'kx');
figure(3)
clf
imshow(g)